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 页)