R中的极大似然估计

人们经常用样本(sample)来估计总体(population)的参数， 假定一个概率密度函数，如正态分布，并给出所有可能的参数组合，那么每一种参数组合所生成当前样本的可能性分别有多大？将最可能产生当前样本的参数，作为总体参数的估计，这就是极大似然(Maximum Likelihood)思想的体现。

在实际的计算中， 人们往往需要先假定某一概率密度， 如正态分布，

$p.d.f. (y) = \frac{1} {\sqrt{2 \pi} \sigma} e^{-1/2 \frac{(y- \mu)^2} {\sigma^2}}$

1
pdfi <- 1/(sqrt(2*pi)*sigma) * exp(-1/2 * (obs - mu)^2/(sigma^2))

其中,

  • obs为某一次观测值
  • mu 为 总体平均值的任意估计值
  • sigma 为 总体标准差的任意估计值，
    则可以求得观测到该次观测值的概率 pdfi

对所有的样本数据的概率密度相成， 即获得联合概率密度， 即似然函数。
似然函数 L = pdf1 * pdf2 * pdf3 * ... * pdfn

$L = \prod_{i=1}^n \frac{1} {\sqrt{2 \pi} \sigma} e^{-1/2 \frac{(y- \mu)^2} {\sigma^2}}$

那么，什么样的参数， 才能保证Likelihood函数取最大值？
由于联合概率密度的值往往非常小，在计算机中难以处理，一般情况下， 将似然函数取对数， 便于数值操作。 此时，做一系列的变换。 ：
如：

$LogL = \sum_{i=1}^n \log\left ( \frac{1} {\sqrt{2 \pi} \sigma} e^{-1/2 \frac{(y- \mu)^2} {\sigma^2}} \right )$

$LogL= \sum_{i=1}^n \left [ \log\left ( \frac{1} {\sqrt{2 \pi} \sigma} \right ) + \log\left ( e^{-1/2 \frac{(y- \mu)^2} {\sigma^2}} \right ) \right ]$

$LogL= \sum_{i=1}^n \log(2 \pi)^{-1/2} + \sum_{i=1}^n \log\left (\frac{1} {\sigma}\right ) + \sum_{i=1}^n \left ( -1/2 \frac{(y- \mu)^2} {\sigma^2} \right )$

$LogL= -\frac{n} {2} \log(2 \pi) - n \log(\sigma) - \frac{1} {2 \sigma^2} \sum_{i=1}^n ( y - \mu)^2$

若认为LogL为$\mu$的函数, 若保证LogL最大值，令 $\sum_{i=1}^n ( y - \mu)^2$ 取最大值，则

$\mu = \frac{\sum_{i=1}^n y_i}{n}$
另外一方面， 若认为LogL为$\sigma$的函数, 将LogL函数对$\sigma$求导， 令导数为0， 可求得极值。

$\frac{\delta LogL}{\delta \sigma} = -\frac{n}{\sigma} - (-2) \frac{1}{2} \sum_{i=1}^n (y - \mu)^2 \sigma^{-3}$

$\frac{\delta LogL}{\delta \sigma} = -\frac{n}{\sigma} + \frac{1}{\sigma^2} \sum_{i=1}^n (y - \mu)^2$

令导数为0， 得
$\frac{n}{\sigma} = \frac{1}{\sigma^3} \sum_{i=1}^n (y - \mu)^2$

$\sigma^2 = \frac{\sum_{i=1}^n (y - \mu)^2}{n}$

此时， 带入之前$\mu$，求得$\sigma$的估计

$\hat{\sigma}^2 = \frac{\sum_{i=1}^n (y - \bar{y})^2}{n}$

用极大似然估计求得的$\sigma$值为有偏估计， 即样本的最可能取值。

以上是解析解。 在实际情况下， 人们常用多种优化方法， 求得极大似然估计的数值解。 这些方法包括

  • "Nelder-Mead" （Nelder and Mead ，1965),
  • "BFGS" （Broyden, Fletcher, Goldfarb and Shanno， 1970）,
  • "CG" （Fletcher and Reeves，1964),
  • "L-BFGS-B" （Byrd et. al. ，1995) ,
  • "SANN" (Belisle 1992), 等。

由此看来， 极大似然估计， 首先是按照假定的概率密度， 求得联合概率密度， 即似然函数，进一步写成对数似然函数， 再用数值优化的方法， 求得对数似然函数所对应的点(各参数)。即获得参数的极大然估计。

下面提供一个实例：
问题：
假设一批树木中抽取若干个体，胸径分别为 150, 124, 100, 107, 170, 144, 113, 108, 92, 129, 123, 118，试求这批树木胸径的平均值和标准差。

分析：
如果用这批样本数据直接计算平均值和方差， 可以得到总体平均值和标准差的估计。

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
dat <- c(150, 124, 100, 107, 170, 144, 113, 108, 92, 129, 123, 118)

> mean(dat)
[1] 123.1667

### 求样本平均值
mean.sample <- sum(dat)/length(dat)

### 设置空向量
var.sample <- c()
for(i in 1:length(dat)){
   var.sample[i] <- (dat[i] - mean.sample)^2
}
dd <- sum(var.sample)

## 求样本标准差
sqrt(dd/(length(dat)))

## 求总体标准差
sqrt(dd/(length(dat) - 1))
以上是普通的估计方法


下面用极大似然方法估计
###########################################################################
### 联合概率密度

like = function(data, mu, sigma) {
 like = 1
 for(i in 1:length(data)){
   # like = like *( 1/(sqrt(2*pi)*sigma) * exp(-1/2 * (data[i] - mu)^2/(sigma^2)))
like = like * dnorm(x = dat[i], mean = mu, sd = sigma)
 }
 return(like)
}

like(dat, 120, 20)

### 对数似然函数
loglike = function(data, mu, sigma) {
 loglike = 0
   for(i in 1:length(data)){
     loglike = loglike + log(dnorm(x = dat[i], mean = mu, sd = sigma))
   }
   return(loglike)
}
loglike(dat, 122, 20)

### 扩展参数组合
params = expand.grid(mu = seq(50, 200, 1),
                    sigma = seq(10, 40, 1))
LogLik.surf = with(params, loglike(dat, mu, sigma))

dim(LogLik.surf) <- c(length(seq(50, 200, 1)), length(seq(10, 40, 1))) 

library(lattice)
contourplot(LogLik.surf ~ mu*sigma, data = params, cuts = 20)

# Zooming in
params = expand.grid(mu = seq(120, 126, 0.01), sigma = seq(18, 25, 0.1))
LogLik.surf.zoom = with(params, loglike(dat, mu, sigma))

dim(LogLik.surf.zoom) <- c(length(seq(120, 126, 0.01)), length(seq(18, 25, 0.1))) 

library(lattice)
contourplot(LogLik.surf.zoom ~ mu*sigma, data = params, cuts = 20, region = TRUE)

### nlm进行极大似然优化， 需要提供初始值
### 对数似然函数
loglike000 = function(p) {
   loglike = 0
   for(i in 1:length(dat)){
     loglike = loglike + log(dnorm(x = dat[i], mean = p[1], sd = p[2]))
   }
   return(-loglike)
}

(MLEst <- nlm(f = loglike000, p = c(100, 10)))

参考

http://www.quantumforest.com/2011/10/maximum-likelihood/