二项分布的对数似然面:算法与R语言实现

对数似然面是极大似然估计的结果之一，这里介绍在R语言中如何绘制对数似然面（Log Likelihood Surface）,这里用到的对数似然函数是二项分布的，对于泊松分布，正态分布，负二项分布等，只需对对数似然函数进行适当修改即可。

绘制二项分布的对数似然面举例：
主要过程

  1. 写出对数似然函数
  2. 用nlm求对数似然函数的估计值
  3. 若绘制似然曲线，则固定一个参数值，并变化另外一个参数，计算Loglikelihood，用lines添加曲线。
  4. 若绘制趋势面，则先要计算参数组合矩阵，并求各栅格处的 LogLikelihood，用persp, image, contour等函数即可。也可以用lattice函数包的wireframe()

以下的例子给出全部的R代码和相应说明。

该例子引自:
http://www.unc.edu/courses/2010fall/ecol/563/001/docs/lectures/lecture10.htm#Rcode

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
## 实例： 假设有以下数据, 记录每个枝条上蚜虫的数量，该数据出自
##Krebs, Charles. 1999. Ecological Methodology.
## Menlo Park, CA: Addison-Wesley. 130页

###########################################################
#######################原始数据############################
###########################################################

### 数据导入

num.stems <- c(6,8,9,6,6,2,5,3,1,4)
aphids_count <- c(0,1,2,3,4,5,6,7,8,9)
aphids <- data.frame(aphids_count, num.stems)

### Raw Data,整理出二项分布拟合所需要的格式
aphids_raw <- rep(aphids_count, num.stems)

###########################################################
####################写出对数似然函数#######################
###########################################################

## 二项分布的对数似然函数，用来绘制对数似然曲线和趋势面

nnominb.loglik <- function(p) {
    res <- sum(log(dnbinom(aphids_raw, mu = p[1], size = p[2])))
    return(res)
}
## 二项分布的负对数似然函数，nlm()函数只能求最小值，因此极大

## 似然值时必须用到此函数
neg.nnominb.loglik <- function(p) {
    res <- sum(log(dnbinom(aphids_raw, mu = p[1], size = p[2])))
    return(-res)

}

###########################################################
###################估计对数似然函数极值####################
###########################################################

#利用nlm估计极值
out.NB <- nlm(neg.nnominb.loglik, c(4,4))
# estimate 是对两个参数的估计，即对neg.nnominb.loglik函数参数p的估计
out.NB

###########################################################
###################图 1 ##绘制似然曲线#####################
###########################################################

 

### x 序列
x <- seq(0.1, 10, by = 0.1)

### 计算y
y_max <- sapply(x, function(x) nnominb.loglik(p = c(out.NB$estimate[1], x)))
plot(y_max ~ x, type = "l", ylab = "Log Likelihood", xlab = "Size")

### 当mu = 2时
y_2 <- sapply(x, function(x) nnominb.loglik(p = c(2, x)))
lines(y_2 ~ x, col = 2)

### 当mu = 5时
y_5 <- sapply(x, function(x) nnominb.loglik(p = c(5, x)))
lines(y_5 ~ x, col = 3)

### 添加图例
legend('bottomright', c(expression(mu=='MLE (3.46)'),
        expression(mu==2), expression(mu==5)),
        col=c(1,2,3), lty=c(1,1,1), cex=.9, bty='n')

###########################################################
##################图 2#绘制栅格图##########################


### 绘制栅格图

## 按照值设色
image(seq(2,5,.1), seq(1,4,.1), z.matrix, col = cm.colors(30),
      xlab = expression(mu), ylab = "Number of Trials")
## 添加等高线
contour(seq(2,5,.1), seq(1,4,.1), z.matrix, add = TRUE)
## 极大似然值
points(x = out.NB$estimate[1], y = out.NB$estimate[2], pch = 3)

###########################################################
###############图 3绘制3D 对数似然面#######################
###########################################################


### 计算趋势面原始数据
### 构建网格
grids <- expand.grid(mu=seq(2,5,.1), theta=seq(1,4,.1))
### 计算z值
grids$z <- apply(grids, 1, function(x) nnominb.loglik(p = c(x[1], x[2])))
### 趋势面原始数据
z.matrix <- matrix(grids$z, nrow=length(seq(2,5,.1)))


### 绘制三维趋势面

persp(seq(2,5,.1), seq(1,4,.1), z.matrix,
      xlab=expression(mu), ylab=expression(theta), zlab='Log Likelihood',
      ticktype='detailed', theta=30, phi=20, d = 4)
###########################################################

##########图 4##绘制对数似然面的三维设色图#################
library(lattice)       
wireframe(z~mu*theta, data=grids, xlab=expression(mu),
          ylab=expression(theta), zlab="Log-nlikelihood",
          scales = list(arrows = FALSE), drape=TRUE,
          par.settings=list(par.zlab.text=list(cex=.75),
          axis.text=list(cex=.75)))