EM算法R代码及疑问.docxVIP

学习了?EM?的理论后想写一下练练手。在网上找到了混合高斯的?R?代码 模拟数据是样本集?5000?个，前?2000?个是以?3?为均值，1?为方差的高斯分布，后 3000?个是以-2?为均值，2?为方差的高斯分布。 1.?#?模拟数据 2.?miu1?-?3 3.?miu2?-?-2 4.?sigma1?-?1 5.?sigma2?-?2 6.?alpha1?-?0.4 7.?alpha2?-?0.6 8.?#?生成两种高斯分布的样本 9.?n?-?5000 10.?x?-?rep(0,n) 11.?n1?-?floor(n*alpha1) 12.?n2?-?n?-?n1 13.?x[1:n1]?-?rnorm(n1)*sigma1?+?miu1 14.?x[(n1+1):n]?-?rnorm(n2)*sigma2?+?miu2 15.?hist(x,freq=F) 16.?lines(density(x),col=red) 17.?#?设置初始值 18.?m?-?2 19.?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?(step?in?1:100){ 25.?#?E?步，已知参数求概率 26.?for?(j?in?1:m){ 27. prob[,j]-?sapply(x,dnorm,miu[j],sigma[j]) 28. } 29. sumprob?-?rowSums(prob) 30. prob-?prob/sumprob ####求出样本属于各类的概率 31. 32. 33. oldmiu?-?miu 34. oldsigma?-?sigma 35. oldalpha?-?alpha 36. 37. 38.?##?M?步?更新参数 39.?for?(j?in?1:m){ 40. p1?-?sum(prob[?,j]) 41. p2?-?sum(prob[?,j]*x) 42. miu[j]?-?p2/p1 43. alpha[j]?-?p1/n 44. p3?-?sum(prob[?,j]*(x-miu[j])^2) 45. sigma[j]?-?sqrt(p3/p1) 46. } 47. 48. 49.?##?结束迭代的条件 50.?epsilo?-?1e-4 51. if?(sum(abs(miu-oldmiu))epsilo? 52. sum(abs(sigma-oldsigma))epsilo? 53. sum(abs(alpha-oldalpha))epsilo)?break 54. 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?(j?in?1:m){ 2. prob[,j]-?sapply(x,dnorm,miu[j],sigma[j]) 3.?} 4. 5.?for(i?in?(1:dim(prob)[1])){ 6. px-alpha%*%prob[i,] 7. prob[i,]-(prob[i,]*alpha)/px 8.?} 这样计算的?prob?才是?z=j?的后验概率。分别贴出来修改前后的计算结果： 修改前：step?35?miu?2.843996?-2.195416?sigma?1.097323?1.849515?alpha 0.4422249?0.5577751 到第?35?次迭代时收敛，得到的均值是?2.843996，-2.195416?（真实值是?3，- 2），方差?1.097323，1.849515（真实值是?1，2）， alpha?0.4422249，0.5577751?（真实值?0.4,0.6）。 修改后的：step?57?miu?-1.924697?3.005561?sigma?2.038401?0.9782731 alpha?0.6028941?0.3971059 到第?57?次迭代时收敛，得到的均值是?3.005561，-1.924697?（真实值是?3，- 2），方差?0