EM 算法 解：豌豆分类 y  ( ) 黄色 绿色 ,  ( , y y 1 ) 2  (84,16) 属于各类的概率分布为： (2 2    ,(1   2 ) ) 其中表示等位基因 A 的概率，且 (0,1)  ，试估计 取的先验分布 ( ) 为(0,1)上的均匀分布，则的观测后验分布为 p ( |  Y )   ( ) ( p Y | )   (2 2    y ) (1 1  )  22 y  y  1 (2   y ) (1 1  )  22 y 令 L ( , Y )  y  1 (2   y ) (1 1 22 y ，对 ( , )   log ( , L  d  L Y d  ) y 1   2 Y 取对数后再对求导得 ) y 1    22 y 1   解似然函数 用如下 R 程序： S   Y ) d log ( , L  d   0 functi  on x   84 { x uniroot S c 84 2 x  ( , (0,1))  2*16 1 x  } 得到 ˆ  0.5999978 . EM 算法：求解.我们将第一类豌豆分为两小类 ( AA Aa ，落入第一小类 AA 的 ) , 概率为 2 ，落入第二小类 Aa 的概率为 2 (1   ，引进潜在变量 Z 表示豌豆属于第一小类 ) 的种数，则 1y Z 为豌豆属于第二小类的种数，从而得到的添加后验分布为 p ( |  , Y Z )   ( ) ( , p Y Z | )  

 2 Z   ( )) ) (2 (1  y Z  1 ((1   2 ) ) y 2  y Z  1  (1  )  y 1  22 y  Z 在第 1i  次迭代中假设有估计值 ( )i ，则可通过 E 步和 M 步得到的一个新的估计。 在 E 步中，有 ) Q  ( ) i (   Y , |   (   E z  y  1  Z y 1 ( E Z z )log ( ) i |   , Y  y ( 1  ) log   2 y    y 1 Z  2 ( )log(1 2  y 2 ) | ( E Z z ( ) i     , Y  ( ) i ))log(1 , Y  |  )  因 Z | Y 服从 ( ) ,i ( B y 1 , ( ) i  ( ) i  ) 2 ，故 ( E Z z ( ) i  | , Y )  ( ) i y  1 ( ) i 2   在 M 步中，我们将 Q ( ) i (   | , Y ) 对求导并令其为 0，有 ( 1) i     , Y ) ( ) i  ) y 2  y 1 | ( E Z z 2( y  1 168 400 200  ( ) i  0 0.5,  ep  1 e  7  { 第 次迭代得到的估计值 : ',   ) theta k   function theta     0 numeric   根据上式编写迭代程序如下： EM theta 1 k { repeat   1 { if k     theta k      if abs theta k } { else   168 / theta k     if abs theta k } k } cat }   ,' ' ( 1  k k ', 168 / 400 200*    theta 0   0 theta  ep break     400 200* theta k   1 ep break het t    a k    1  

运行得到 0.59 ˆ  999 99 . 上述两种方法运行结果截图： EM 算法运行结果： 完整的源代码如下 S <- function(x){ } uniroot(S,c(0,1)) 84/x - 84/(2-x) - 32/(1-x) #用 uniroot 函数求解一元方程在区间[0,1]上的根 EM <- function(theta0=0.5,ep=1e-7){ theta <- numeric(0) k=1 repeat{ if(k==1){ } cat('第',k,'次迭代得到的估计值:',theta[k]) } theta[k]=168/(400-200*theta0) if(abs(theta[k]-theta0)< ep)break theta[k]= 168/(400-200*theta[k-1]) if(abs(theta[k]-theta[k-1])< ep)break } else{ } k=k+1 

#测试 theta=0.5 ep=1e-7 EM(theta,ep)