从黑球、白球到极大似然估计

Likelihood 是生态和进化数据分析中的重要概念， 表示在给定模型和参数时， 某一事件发生的可能性。最大似然估计，就是已知某一事件发生的状况，计算在给定模型时， 模型参数取什么样的值， 才能让当前事件的可能性最大。

这里的事件， 可以是从袋子里抽取一定次数， 所得黑白球的结果，也可以是比对好的DNA碱基序列， 也可以是物种丰富度和降水量对应关系数据等等。 因此说， 所谓的极大似然估计， 其实是对参数取值的估计。

在分析中，如果一件事情已经发生，在计算其Likelihood时，就是将理论上的各种状态相乘，获得联合概率密度。求该联合概率密度，在参数取什么值时， 保证联合概率密度最大， 就是极大似然估计。

例如： 假设从一个袋子里，随机抽取白球和黑球，每一次抽完， 都放回到袋子里，共抽取了100次， 得到80个黑球， 20个白球， 求袋子中， 黑球与白球的比例是多大时， 最有可能得到当前的结果？

分析： 袋子里的球， 非黑即白， 可以简化为， 二项分布， 即一件事情发生还是不发生的概率。 假设抽到黑球的概率为p， 那么抽到白球的概率为1-p。在本例子中，抽取100次， 抽到了80个黑球， 其概率为p^80, 而抽到白球的概率为 1-p, 发生了20次， 因此概率为 (1 - p)^20。 因此， 产生抽样结果的概率为p^80*(1 - p)^20，这就是联合概率密度。为了便于批量计算， 这里写成了R函数。

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
### Likelihood函数， 用来计算在p取值不同时， 当前结果（即80个黑球， 20个白球）发生的可能性。 
likelihood <- function(p) {p^80*(1 - p)^20} 

###生成p的取值序列，p值取0.01， 意思是 100个白球里面， 有1个黑球， 99个黑球。
### 不难想象， 这种情况下， 抽取100次， 抽到80个黑球的可能性一定很小。 
p <- seq(0.01, 0.99, by = 0.001) 
 
### 计算p不同取值时， 抽取到80个黑球， 20个白球的可能性（likelihood）。
res.like <- likelihood(p) 

### 绘图
plot(res.like, type = "l", col = "red")


### 由于likelihood是概率的乘机， 一般都很小， 所以在取自然对数之后，更加方便计算机处理， 所以一般都取自然对数。
### 这就是loglikelihood。
loglikelihood <- function(p) {log(p^80*(1 - p)^20)} 

### 计算log likelihoodres.loglike <- loglikelihood(p)
### 绘图， 可查看loglikelihood与参数p的关系
plot(res.loglike, type = "l", col = "red")


### 用optimize可以通过数值方法，求函数的极值。
optimize(loglikelihood, interval = c(0, 1),maximum = TRUE)
## $maximum
## [1] 0.7999832
## 
## $objective
## [1] -50.04024

由于获得likelihood的最大值的解析解常常十分困难， 实际情况在求maximum likelihood时， 大部分会用到数值方法。而且也不一定是求最大值， 而是求函数的最小值。 具体可参考R软件的 optim, nlm等函数的帮助文档。