共计 2121 个字符，预计需要花费 6 分钟才能阅读完成。
原文链接：http://tecdat.cn/ 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
require(mvtnorm)
N = nrow(X)
mu = params$mu
var = params$var
probs = params$probs
ri = matrix(0, ncol=clusters, nrow=N)
ll = 0
it = 0
converged = FALSE
if (showits)
cat(paste(“Iterations of EM:”, “\n”))
while (!converged & it < maxits) {
probsOld = probs
llOld = ll
riOld = ri
# Compute responsibilities
for (k in 1:clusters){
ri[,k] = probs[k] * dmvnorm(X, mu[k,], sigma = var[[k]], log=F)
}
ri = ri/rowSums(ri)
rk = colSums(ri)
probs = rk/N
for (k in 1:clusters){
varmat = matrix(0, ncol=ncol(X), nrow=ncol(X))
for (i in 1:N){
varmat = varmat + ri[i,k] * X[i,]%*%t(X[i,])
}
mu[k,] = (t(X) %*% ri[,k]) / rk[k]
var[[k]] = varmat/rk[k] – mu[k,]%*%t(mu[k,])
ll[k] = -.5 * sum(ri[,k] * dmvnorm(X, mu[k,], sigma = var[[k]], log=T) )
}
ll = sum(ll)
parmlistold = c(llOld, probsOld)
parmlistcurrent = c(ll, probs)
it = it + 1
if (showits & it == 1 | it%%5 == 0)
cat(paste(format(it), “…”, “\n”, sep = “”))
converged = min(abs(parmlistold – parmlistcurrent)) <= tol
}
clust = which(round(ri)==1, arr.ind=T)
clust = clust[order(clust[,1]), 2]
out = list(probs=probs, mu=mu, var=var, resp=ri, cluster=clust, ll=ll)
}
后果，能够用图像化来示意：
qplot(x=xs, y=ys, data=Y)
ggplot(aes(x=xs, y=ys), data=Y) +
geom_point(aes(color=factor(test$cluster)))
相似的其余状况这里不出现了，另外 r 语言提供了 EMCluster 包能够比拟不便的实现 EM 进行参数估计和后果的误差剖析。
ret <- init.EM(Y, nclass = 2)
em.aic(x=Y,emobj=list(pi = ret$pi, Mu = ret$Mu, LTSigma = ret$LTSigma))# 计算结果的 AIC
通过比拟不同状况的 AIC，咱们能够筛选出适宜的聚类数参数值。
