R语言实现：混合正态分布EM最大期望估计法

原文链接：http://tecdat.cn/?p=4815

 

因为近期在分析数据时用到了EM最大期望估计法这个算法，在参数估计中也用到的比较多。然而，发现国内在R软件上实现高斯混合分布的EM的实例并不多，大多数是关于1到2个高斯混合分布的实现，不易于推广，因此这里分享一下自己编写的k个高斯混合分布的EM算法实现请大神们多多指教。并结合EMCluster包对结果进行验算。

      本文使用的密度函数为下面格式：

   对应的函数原型为 em.norm(x,means,covariances,mix.prop)

x为原数据，means为初始均值，covariances为数据的协方差矩阵，mix.prop为混合参数初始值。

使用的数据为MASS包里面的synth.te数据的前两列

x <- synth.te[,-3]

首先安装需要的包，并读取原数据。

1. install.packages("MASS")

2. 


3. library(MASS)

4. 


5. install.packages("EMCluster")

6. 


7. library(EMCluster)

8. 


9. install.packages("ggplot2")

10. 


11. library(ggplot2)

12. 


13. Y=synth.te[,c(1:2)]

14. 


15. qplot(x=xs, y=ys, data=Y) 

然后绘制相应的变量相关图：

从图上我们可以大概估计出初始的平均点为(-0.7,0.4) (-0.3,0.8)(0.5,0.6)

当然 为了试验的严谨性，我可以从两个初始均值点的情况开始估计

首先输入初始参数：

1. mustart = rbind(c(-0.5,0.3),c(0.4,0.5))    

2. 


3. covstart = list(cov(Y), cov(Y))

4. 


5. probs = c(.01, .99)

 

然后编写em.norm函数，注意其中的clusters值需要根据不同的初始参数进行修改，

1. em.norm = function(X,mustart,covstart,probs){

2. 


3. 


4. 


5.   params = list(mu=mustart, var=covstart, probs=probs)   

6. 


7.   clusters = 2 

8. 


9.   tol=.00001

10. 


11.   maxits=100

12. 


13.   showits=T

14. 


15.   require(mvtnorm)

16. 


17. 


18. 


19.   N = nrow(X)

20. 


21.   mu = params$mu

22. 


23.   var = params$var

24. 


25.   probs = params$probs

26. 


27. 


28. 


29. 


30. 


31.   ri = matrix(0, ncol=clusters, nrow=N)         

32. 


33.   ll = 0                                        

34. 


35.   it = 0                                         

36. 


37.   converged = FALSE                            

38. 


39. 


40. 


41.   if (showits)                                 

42. 


43.     cat(paste("Iterations of EM:", "\n"))

44. 


45. 


46. 


47.   while (!converged & it < maxits) { 

49.     probsOld = probs

51.     

53.     llOld = ll

55.     riOld = ri

57.     

59.    

61.     # Compute responsibilities

63.     for (k in 1:clusters){

65.       ri[,k] = probs[k] * dmvnorm(X, mu[k,], sigma = var[[k]], log=F)

67.     }

69.     ri = ri/rowSums(ri)

71.     

73.   

75.     rk = colSums(ri)                             

77.     probs = rk/N

79.     for (k in 1:clusters){

81.       varmat = matrix(0, ncol=ncol(X), nrow=ncol(X))         

83.       for (i in 1:N){

85.         varmat = varmat + ri[i,k] * X[i,]%*%t(X[i,])

87.       }

89.       mu[k,] = (t(X) %*% ri[,k]) / rk[k]

91.       var[[k]] =  varmat/rk[k] - mu[k,]%*%t(mu[k,])

93.       ll[k] = -.5 * sum( ri[,k] * dmvnorm(X, mu[k,], sigma = var[[k]], log=T) )

95.     }

97.     ll = sum(ll)

99.     

101.      

103.     parmlistold =  c(llOld, probsOld)            

105.     parmlistcurrent = c(ll, probs)             

107.     it = it + 1

109.     if (showits & it == 1 | it%%5 == 0)         

111.       cat(paste(format(it), "...", "\n", sep = ""))

113.     converged = min(abs(parmlistold - parmlistcurrent)) <= tol

115.   }

117.   

119.   clust = which(round(ri)==1, arr.ind=T)       

121.   clust = clust[order(clust[,1]), 2]           

123.   out = list(probs=probs, mu=mu, var=var, resp=ri, cluster=clust, ll=ll)

125. } 

结果，可以用图像化来表示：

1. qplot(x=xs, y=ys, data=Y) 

2. 


3. ggplot(aes(x=xs, y=ys), data=Y) +

4. 


5.    geom_point(aes(color=factor(test$cluster))) 

 

 

 类似的其他情况这里不呈现了，另外r语言提供了EMCluster包可以比较方便的实现EM进行参数估计和结果的误差分析。

1. ret <- init.EM(Y, nclass = 2)

2. 


3. em.aic(x=Y,emobj=list(pi = ret$pi, Mu = ret$Mu, LTSigma = ret$LTSigma))#计算结果的AIC

通过比较不同情况的AIC，我们可以筛选出适合的聚类数参数值。（欢迎转载，请注明出处。 ）