该例子引自:

http://www.unc.edu/courses/2010fall/ecol/563/001/docs/lectures/lecture10.htm#Rcode

## 实例： 假设有以下数据, 记录每个枝条上蚜虫的数量，该数据出自

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

### 添加图例