以下这个EM算法例子来自我的导师Prof. Yang, 因为利用了指数分布无记忆性, 因而十分简洁有力。例子如下: 从指数总体 exp(lambda) 观察到下述10个数字, 其中第9,10个观察被删失 0.22 0.86 0.70 1.59 0.86 2.47 0.70 0.77 1.2+ 1.0+ 试问如何估计lambda?
############################################################
## EM algorithm for sovling MLE from censored
## observations from exponential distribution
## 10 observations from exponential distribution exp(lambda)
## 0.22 0.86 0.70 1.59 0.86 2.47 0.70 0.77 1.2+ 1.0+
############################################################
x=c(0.22, 0.86, 0.70, 1.59, 0.86, 2.47, 0.70, 0.77, 1.2, 1.0)
# naive estimator lambda.hat= 1/mean(x)=0.964 ##偏大
# solve lambda by EM algorithm:
epsilon=Inf
lambda=lambda.init = 100 # 初值
while (epsilon>1e-5){
#x9.expected = lambda + 1.2
#x10.expected = lambda + 1.0
x.new=c(x[1:8], x[9:10]+lambda)
lambda.new = 1/mean(x.new)
epsilon=abs(lambda-lambda.new)
print(epsilon)
lambda=lambda.new
}
cat("Initial lambda =", lambda.init, ";\t Converged lambda =", lambda, "\n")
############################################################
一个有意思的扩展是, 如果后9个数据都是被删失的, 即数据为 0.22 0.86+ 0.70+ 1.59+ 0.86+ 2.47+ 0.70+ 0.77+ 1.2+ 1.0+ 需要将上述程序中的 x.new=c(x[1:8], x[9:10]+lambda) 改为 x.new=c(x[1], x[2:10]+lambda) 即可, 你会看到此时程序收敛的慢一些。更进一步地问, 如果数据全是删失的, 上述EM算法还能用吗?