极大似然估计与牛顿拉普森极值算法

03 May 2013

---

---

题目有点大了，其实就想说一个简单的事。前两天，一同学找我，问我能不能用极大似然估计来求解一个分布参数，我想了良久，最后从网上找到了答案。极大似然的关键是写出密度函数的乘积形式，然后取对数，剩下的工作就是求这个对数似然函数的极值点。如何求解最值对应的参数点最为棘手，但R中已提供了求解函数最值得函数–nlm.下面举一例子说明其应用。

现有一负二项分布样本，容量为100。

>set.seed(1)
> x <- rnbinom(100, 9, 0.6)
> x
  [1]  7  5  4  5  5  5  5  8  5  6  5  7  2  5  6  3 10  5  2 11  8 10 11  3
 [25]  5  4 10  6  5  8 15  6  1  5  5  4  3  3 12  6  4  6  5  5  2 11  4  5
 [49]  7  6  5  6  8  5  4  5  3  2  2 10  9  4  6  7  6  5  6  8  4 14  3  8
 [73]  4  5  1  1  3  6  4  2  9  2  4  3 14  7  7  4  6 17  9  9  1 12  5  4
 [97]  3 13  7  8

分布参数r,p未知，使用极大似然方法进行估计。首先写出对数似然函数。

fun <- function(par, x) {
  r <- par[1]
  p <- par[2]
  r <-   sum(log(dnbinom(x, r, p)))
  r1 <- -r
  return(r1)
}

注意到，函数fun的返回值为对数似然函数的相反数，这是因为nlm默认是求函数极小值点，如果我们要求极大值，就需要取其相反数（加负号）。

我们把此函数带入nlm中如下，其中，c(9, 0.2)为初始值，可以自由设置，只要别太离谱就行。

> es <- nlm(fun, c(9,0.2), x)
> es
$minimum
[1] 249.6542

$estimate
[1] 9.1897195 0.6065934

$gradient
[1] -3.555327e-06  1.262791e-04

$code
[1] 1

$iterations
[1] 20

函数的极小值为249.6542，极小值点为r=9.1897195,p= 0.6065934,意味着对数似然函数极大值为-249.6542，参数估计值为r=9.1897195,p= 0.6065934。可见参数估计值与我们最初生成样本时使用的参数很接近，也就是说，这种方法很不坏！

---

图形会使枯燥的数字生动。把对数似然函数用等高线表示在一个二维平面上。

> fun1 <- function(par, x) {
+   r <- par[1]
+   p <- par[2]
+   r <-   sum(log(dnbinom(x, r, p)))
+   r1 <- -r
+   return(r)
+ }
> r <- seq(1, 15, length = 40)
> p <- seq(0.4, 0.9, length = 40)
> m.s <- expand.grid(r,p)
> z.matrix1 <- apply(m.s, 1, fun1, x)
> 
> grids <- cbind(m.s, z.matrix1)
> names(grids) <- c("r", "p", "z")
> z.matrix <- matrix(z.matrix1, nrow = 40)
> image(r, p, z.matrix, col=cm.colors(30))
> contour(r,p, z.matrix, nlevels = 100, add = T)
> points(model$es[1],model$es[2] , pch = 15, col = 2, cex = 1)

其中，红色方块表示极大值点位置。

也可以把对数似然函数表示在一个三维立体图中,如下。

> library(lattice)       
> wireframe(z~r*p, data=grids, xlab=expression(r),
+           ylab=expression(p), zlab="Log-\nlikelihood",
+           scales = list(arrows = FALSE), drape=TRUE,
+           par.settings=list(par.zlab.text=list(cex=.75),
+           axis.text=list(cex=.75)), screen = list(z = -30, x = -60))

---

参考：

• R中的极大似然估计
• nlm,R帮助文档

---

---