R 程式設計/最大似然

簡介

最大似然估計只是一個最佳化問題。您需要寫下您的對數似然函式並使用一些最佳化技術。有時您還需要編寫您的得分函式（對數似然的導數）或 Hessian 矩陣（對數似然的二階導數）。

一維

如果只有一個引數，我們可以使用optimize()來最佳化對數似然函式。

型別 1 Pareto 分佈示例

我們提供了一個型別 1 Pareto 分佈的示例。請注意，在這個示例中，我們將最小值視為已知，並且不進行估計。因此，這是一個一維問題。

我們使用rpareto1()（actuar）函式從形狀為 1，最小值為 500 的型別 1 Pareto 分佈中生成一個隨機向量。我們使用dpareto1()（actuar）函式，並設定log = TRUE選項來編寫對數似然函式。然後，我們只需要使用optimize()，並設定maximum=TRUE。我們使用interval選項為引數提供最小值和最大值。

> library(actuar)
> y <- rpareto1(1000, shape = 1, min = 500)
> ll <- function(mu, x) { 
+    sum(dpareto1(x,mu[1],min = min(x),log = TRUE)) 
+   } 
> optimize(f = ll, x = y, interval = c(0,10), maximum = TRUE)

多維

• fitdistr()（MASS 包）透過最大似然法擬合單變數分佈。它是optim()的包裝器。
• 如果您需要自己編寫最大似然估計器（MLE），則必須使用內建最佳化器，例如nlm(), optim(). R 還包含以下最佳化器 
• stats4 包中的mle()
• maxLik 包

Logistic 分佈示例

例如，我們從 Logistic 分佈中抽取樣本，並使用 . 來估計引數。

> # draw from a gumbel distribution using the inverse cdf simulation method
> e.1 <- -log(-log(runif(10000,0,1))) 
> e.2 <- -log(-log(runif(10000,0,1)))
> u <- e.2 - e.1  # u follows a logistic distribution (difference between two gumbels.)
> fitdistr(u,densfun=dlogis,start=list(location=0,scale=1))

Cauchy 分佈示例

例如，我們可以使用nlm()最佳化器為 Cauchy 分佈編寫一個簡單的最大似然估計器。我們首先從 Cauchy 分佈中抽取一個向量 x。然後我們定義對數似然函式，然後使用nlm()函式進行最佳化。請注意nlm()是極小值點，而不是極大值點。

> n <- 100
> x <- rcauchy(n)
> mlog.1 <- function(mu, x) { 
+   - sum(dcauchy(x, location = mu, log = TRUE)) 
+   } 
> mu.start <- median(x)
> out <- nlm(mlog.1, mu.start, x = x)

Beta 分佈示例

這是另一個使用 Beta 分佈和optim()函式的示例。

> y <- rbeta(1000,2,2)
> loglik <- function(mu, x) { 
+    sum(-dbeta(x,mu[1],mu[2],log = TRUE)) 
+    } 
> 
> out <- optim(par = c(1,1), fn=loglik,x=y,method = "L-BFGS-B",lower=c(0,0))

檢驗

似然比檢驗

• lrtest()在lmtest包中。

一些具體案例

• gum.fit()（ismev 包）提供 Gumbel 分佈的 MLE

資源

• Charles Geyer：R 中的最大似然法 (pdf，9 頁)
• Marco Steenbergen：R 中的最大似然法程式設計 (pdf，7 頁)