不写R包的分析师不是好全栈

MLE在R中的实现

    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)

plot of chunk unnamed-chunk-2



第一个方法:


Newton-Raphen


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/mean(x)来计算(不全是MLE估计,使用了一部分矩估计的性质)


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

plot of chunk unnamed-chunk-5


参数对比:






















methodalphalambda
Moment3.0750.5235
mle-Newton2.8550.4860




方法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]))

plot of chunk unnamed-chunk-7


参数对比:



























methodalphalambda
Moment3.0750.5235
MLE-Newton2.8550.4860
MLE-optim2.8540.4859

page PV:  ・  site PV:  ・  site UV: