EM算法R代码及疑问

学习了EM的理论后想写一下练练手。在网上找到了混合高斯的R代码模拟数据是样本集5000个，前2000个是以3为均值，1为方差的高斯分布，后3000个是以-2为均值，2为方差的高斯分布。1. #模拟数据2. miu1-33. miu2-24. sigma1-15. sigma2-26. alpha1-0.47. alpha2-0.68. #生成两种高斯分布的样本9. n-500010. x-rep(0,n)11. n1-floor(n*alpha1)12. n2-n-n113. x1:n1-rnorm(n1)*sigma1+miu114. x(n1+1):n-rnorm(n2)*sigma2+miu215. hist(x,freq=F)16. lines(density(x),col=red)17. #设置初始值18. m-219. miu-runif(m)20. sigma-runif(m)21. alpha-c(0.2,0.8)22. prob-matrix(rep(0,n*m),ncol=m)23. #迭代24. for(stepin1:100)25. #E步，已知参数求概率26. for(jin1:m)27. prob,j-sapply(x,dnorm,miuj,sigmaj)28. 29. sumprob-rowSums(prob)30. prob-prob/sumprob#求出样本属于各类的概率31. 32. 33. oldmiu-miu34. oldsigma-sigma35. oldalpha-alpha36. 37. 38. #M步更新参数39. for(jin1:m)40. p1-sum(prob,j)41. p2-sum(prob,j*x)42. miuj-p2/p143. alphaj-p1/n44. p3-sum(prob,j*(x-miuj)2)45. sigmaj-sqrt(p3/p1)46. 47. 48. 49. #结束迭代的条件50. epsilo-1e-451. if(sum(abs(miu-oldmiu)epsilo&52. sum(abs(sigma-oldsigma)epsilo&53. sum(abs(alpha-oldalpha)epsilo)break54. cat(step,step,miu,miu,sigma,sigma,alpha,alpha,n)55. 56. 57. 有疑问的是在E步里面求样本属于各个类的概率的步骤，上面给出的代码是用P(Xi|z=j)/P(xi|zj)计算。但是在看理论的时候E步计算的应该是P(z=j|xi)，也就是z=j的后验概率。用贝叶斯公式，后验概率可以从先验概率计算得到。（这玩意怎么编辑公式？只能截图了）其中P(x|z)是高斯分布的概率密度，P(z)就是代码中的alpha。于是我修改了一下上面代码中E步的prob参数计算的方法：1. for(jin1:m)2. prob,j-sapply(x,dnorm,miuj,sigmaj)3. 4. 5. for(iin(1:dim(prob)1)6. px-alpha%*%probi,7. probi,-(probi,*alpha)/px8. 这样计算的prob才是z=j的后验概率。分别贴出来修改前后的计算结果：修改前：step 35 miu 2. -2. sigma 1. 1. alpha 0. 0.到第35次迭代时收敛，得到的均值是2.，-2. （真实值是3，-2），方差1.，1.（真实值是1，2），alpha0.，0. （真实值0.4,0.6）。修改后的：step 57 miu -1. 3. sigma 2. 0. alpha 0. 0.到第57次迭代时收敛，得到的均值是3.，-1. （真实值是3，-2），方差0.，2.0384（真实值是1，2），alpha 0.3971，0.6029 （真实值0.4,0.6）。显然修改后的参数估计更接近实际值。其实如果alpha=（0.5,0.5）的话，修改后的计算公式就是等于修改前的。也就是说，如果每个类别的样本数相差不大的话，修改前的计算结果应该是接近修改后的，反之则不然。好了，下面开始喜闻乐见的实验时间。之前模拟数据一组是2000，一组3000个。两者相差不大，现在吧数量差距加大，一个1000，一个4000。为了使数据分得开，把sigma2调小了一点。1. miu1-32. miu2-23. sigma1-14. sigma2-1.55. alpha1-0.26. alpha2-0.87. #生成两种高斯分布的样本8. n-50009. x-rep(0,n)10. n1-floor(n*alpha1)11. n2-n-n112. x1:n1-rnorm(n1)*sigma1+miu113. x(n1+1):n-rnorm(n2)*sigma2+miu214. hist(x,freq=F)15. lines(density(x),col=red)然后计算，修改前的代码计算的结果：step 30 miu 0. -2. sigma 2. 1. alpha 0. 0.可以看到偏离真实值非常夸张，根本不可用了。修改后的结果呢：step 64 miu 2.95897 -1. sigma 1