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

因为近期在分析数据时用到了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]

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

install.packages("MASS")

library(MASS)

install.packages("EMCluster")

library(EMCluster)

install.packages("ggplot2")

library(ggplot2)

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

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

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

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

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

首先输入初始参数：

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

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

probs = c(.01, .99)

然后编写em.norm函数，注意其中的clusters值需要根据不同的初始参数进行修改，em.norm = function(X,mustart,covstart,probs){

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

clusters = 2

tol=.00001

maxits=100

showits=T