R语言与点估计学习笔记.pdfVIP

R 语言与点估计学习笔记（EM 算法与 Bootstrap 法） 一、EM 算法 EM 算法是一种在观测到数据后，用迭代法估计未知参数的方法。可以证明EM 算法 得到的序列是稳定单调递增的。这种算法对于截尾数据或参数中有一些我们不感兴趣的参 数时特别有效。 EM 算法的步骤为： E-step （求期望）：在给定y 及theta=theta （i）的条件下，求关于完全数据对数似 然关于潜在变量z 的期望 M-step （求极值）：求上述期望关于theta 的最大值theta （i+1） 重复以上两步，直至收敛即可得到theta 的MLE。 从上面的算法我们可以看到对于一个参数的情况，EM 仅仅只是求解MLE 的一个迭代算 法。M-step 做得就是optimize 函数做得事情。对于EM 算法，我们也没有现成的求解函数 （这个是自然的），我们一样可以通过人机交互的办法处理。 先举一个一元的例子： 设一次实验可能有4 个结果，发生概率分别为0.5+theta/4, 0.25-theta/4 ,0.25- theta/4 ,theta/4.其中theta 在0,1 之间。现进行了197 次实验，结果发生的次数分别为： 125,18,20,34，求theta 的MLE。 计算出theta （i+1）= （195theta （i）+68）/ （197theta （i）+144） 为什么是这个结果，请翻阅王兆军《数理统计讲义》p43-p44 我们用简单的循环就可以解决这个问题，程序及结果如下： fun-function(error=1e-7){ +theta-0.5 +k-1 +while(T){ +k-k+1 +theta[k]-(159*theta[k-1]+68)/(197*theta[k-1]+144) +if(abs(theta[k]-theta[k-1])error) break +} +list(theta-theta[k],iter-k) +} fun() [[1]] [1]0.6268215 [[2]] [1]9 我们再看一个二元的简单例子： 幼儿园里老师给a，b，c，d 四个小朋友发糖吃，但老师有点偏心，不同小朋友得 到糖的概率不同，p(a)=0.5,p(b)=miu, p(c)=2*miu, p(d)=0.5-3*miu 如果确定了参数miu， 概率分布就知道了。我们可以通过观察样本数据来推测参数知道c 和d 二人得到的糖果 数，也知道a 与b 二人的糖果数之和为h，如何来估计出参数miu 呢？前面我们知道了， 如果观察到a，b，c，d 就可以用ML 估计出miu。反之，如果miu 已知，根据概率期望 a/b=0.5/miu，又有a+b=h 。由两个式子可得到 a=0.5*h/(0.5+miu)和b=miu*h/(0.5+miu)。 # 已知条件 h = 20 c = 10 d = 10 # 随机初始两个未知量 miu = runif(1,0,1/6) b = round(runif(1,1,20)) iter = 1 nonstop=TRUE while (nonstop) { + # E 步骤，根据假设的miu 来算b + b = c(b,miu[iter]*h/(0.5+miu[iter])) + print(b) + # M 步骤，根据上面算出的b 再来计算miu + miu = c(miu,(b[iter+1] +c)/(6*(b[iter+1]+c+d))) + print(miu) + # 记录循环次数 + iter = iter + 1 + # 如果前后两次的计算结果差距很小则退出 + nonstop =((miu[iter]-miu[iter-1])10^(-10)) +} [1]3.000000 4.450531 [1]00 print(cbind(miu,b)) miu b [1,]03.000000 [2,]04.450531 关于EM 算法，及后续的发展GME 的理论你可以在多数数理统计书上找到相关结论，也 可以用类似办法编写函数处理它。 二、 自助法（bootstrap） Bootst