MLE在R中的实现

Theory Part

In statistics, maximum-likelihood estimation (MLE) is a method of estimating the parameters of a statistical model. When applied to a data set and given a statistical model, maximum-likelihood estimation provides estimates for the model’s parameters.

MLE(Most likelihood estimation)极大似然估计,是是一个用来估计一个统计模型中参数的方法.在给定数据和统计的模型(或者是数据服从的某种分布),MLE体统了对于模型(分布中未知参数)的估计.

MLE是数据统计中无法避开不谈的问题之一,MLE的思想是极大化似然函数从而得到系数的估计. MLE如此基础的原因之一是在于它的使用面很广,从估计分布的未知参数到估计模型的参数,无处不在.而且,基于MLE的一些推广,比如EM算法等也在统计学中占有举足轻重的地位.

先来看下MLE的基本流程:

假设(X_1,X_2,…,X_n)独立同分布于一个连续分布函数(f),对于一组观测值(x_1,x_2,…,x_n),产生这组观测值的概率为:

[P(x_1<X_1\le x_1+dx,,…,x_n<X_n\le x_n+dx)] [= \prod_{i=1}^n P(x_i<X_i\le x_i+dx)=\prod_{i=1}^nf(x_i)]

在这个式子中假设随机变量取值在一个长度为(dx)的小区间内,在观测值确定的情况下(dx) 趋近于0.因此CDF函数可以转化为PDF函数,从而得到最后的(\prod_{i=1}^nf(x_i)dx).接下来做的是:

[\hat{\theta} = Argmax{\prod_{i=1}^nf(x_i)dx}] Argmax代表的是极大化(\prod_{i=1}^nf(x_i)),从而得到参数的估计,就是(\hat{\theta}).(\prod_{i=1}^nf(x_i)dx) 就是似然函数,不过因为它是和未知参数(\theta)有关的,所以一般会写作: [L(\theta)=f(x_1,x_2,\ldots,x_n|\theta)= \prod_{i=1}^nf(x_i)] 总的来说,就是MLE是一个最大化似然函数来计算(\theta)的过程,而极大化这个似然函数的意义在于使得在我们估计出(\theta)的情况下,产生这组数据的概率最大.

Simulate Part

如果你没能看懂上面这堆东西是啥,恭喜你,你的生活一定很快乐..

我其实是想讲MLE在R中的实现.

要在R中实现,其实是把MLE中理论推导的过程撇开,不讨论求偏导,解方程,求极值从而极大化一个函数,而是改用数值的方法来实现似然函数的极大化.

数据:采用1697-1999年 北美1720个位置的降雨信息.

对于降雨数据,一般采用Gamma分布来进行建模比较理想,这里,使用矩估计的方法来对gamma分布的参数进行估计.

#install.packages("spuRs")
library(spuRs)
data(kew)
names(kew)

##  [1] "year" "jan"  "feb"  "mar"  "apr"  "may"  "jun"  "jul"  "aug"  "sep"
## [11] "oct"  "nov"  "dec"

kew[, 2:13] <- kew[, 2:13]/10
kew.mean <- apply(kew[-1], 2, mean)
kew.var <- apply(kew[-1], 2, var)
lambda.mm <- kew.mean/kew.var
m.mm <- kew.mean^2/kew.var

比较直方图和通过矩估计得到的参数曲线:

c(m.mm[7],lambda.mm[7])

##    jul    jul
## 3.0750 0.5235

hist(kew$jun,breaks=25,freq=F)
t=seq(0,length.out = 200,15)
lines(t,dgamma(t,m.mm[7],lambda.mm[7]),lty=2)

x <- kew$jul
x[x == 0] <- 0.1


newtonraphson <- function(ftn, x0, tol = 1e-9, max.iter = 100, …) {
 # find a root of ftn(x, …) near x0 using Newton-Raphson
 # initialise
 x <- x0
 fx <- ftn(x, …)
 iter <- 0
 # continue iterating until stopping conditions are met
 while ((abs(fx[1]) > tol) && (iter < max.iter)) {
 x <- x - fx[1]/fx[2]
 fx <- ftn(x, …)
 iter <- iter + 1
 }
 # output depends on success of algorithm
 if (abs(fx[1]) > tol) {
 stop("Algorithm failed to converge\n")
 } else {
 return(x)
 }
 }

dl <- function(m, a) {
 return(c(log(m) - digamma(m) - a, 1/m - trigamma(m)))
 }

m.ml <- newtonraphson(dl, m.mm[7], a = log(mean(x)) - mean(log(x)))
 lambda.ml <- m.ml/mean(x)


hist(kew$jun,breaks=25,freq=F)
lines(t,dgamma(t,m.mm[7],lambda.mm[7]),lty=2)
lines(t, dgamma(t, m.ml, lambda.ml))

方法2

• optim函数完成最大化

ell <- function(theta, x) {
 alpha <- theta[1]
 beta <- theta[2]
 sum(log(dgamma(x,alpha,beta)))
}
alpha=1;beta=1
ell.optim <- optim(c(alpha, beta), ell, x = x,
 control = list(fnscale = -1,ndeps=1e-9))

## Warning: NaNs produced

## Parameter:
ell.optim$par

## [1] 2.8544 0.4859

hist(kew$jun,breaks=25,freq=F)
lines(t,dgamma(t,m.mm[7],lambda.mm[7]),lty=2)
lines(t, dgamma(t,ell.optim$par[1], ell.optim$par[2]))

参数对比: