用EM算法来估计GMM的参数。
# 设置初始值
m <- 2
miu <- runif(m)
sigma <- runif(m)
alpha <- c(0.2,0.8)
prob <- matrix(rep(0,n*m),ncol=m)
for (step in 1:100){
# E步骤
for (j in 1:m){
prob[,j]<- sapply(x,dnorm,miu[j],sigma[j])
}
sumprob <- rowSums(prob)
prob<- prob/sumprob
oldmiu <- miu
oldsigma <- sigma
oldalpha <- alpha
# M步骤
for (j in 1:m){
p1 <- sum(prob[ ,j])
p2 <- sum(prob[ ,j]*x)
miu[j] <- p2/p1
alpha[j] <- p1/n
p3 <- sum(prob[ ,j]*(x-miu[j])^2)
sigma[j] <- sqrt(p3/p1)
}
# 变化
epsilo <- 1e-4
if (sum(abs(miu-oldmiu))<epsilo &
sum(abs(sigma-oldsigma))<epsilo &
sum(abs(alpha-oldalpha))<epsilo) break
cat('step',step,'miu',miu,'sigma',sigma,'alpha',alpha,'\n')
}
在33次循环之后运算结果趋于稳定，估计的miu为（-2.2，2.8），sigma为（1.82，1.14）
GMM 模型常用于基于模型的聚类分析，GMM中的每一个高斯分布都可以代表数据的一类，整个数据就是多个高斯分布的混合。在R中的mclust包中的Mclust函数可以用来进行基于GMM的聚类分析。下面即是以最常用的iris数据集为例，聚类结果生成的图形就是文章的第一幅图：
library(mclust)
mc <- Mclust(iris[,1:4], 3)
plot(mc, data=iris[,1:4], what="classification",dimens=c(3,4))
table(iris$Species, mc$classification)
扫码加好友，拉您进群
全部版块
我的主页
收藏
