EM算法有很多的应用:
最广泛的就是GMM混合高斯模型、聚类、HMM等等.
The EM Algorithm
高斯混合模型（Mixtures of Gaussians）和EM算法
EM算法
求最大似然函数估计值的一般步骤：
（1）写出似然函数；
（2）对似然函数取对数，并整理；
（3）求导数，令导数为0，得到似然方程；
（4）解似然方程，得到的参数即为所求.
期望最大化算法（EM算法）:
优点：
1、 简单稳定；
2、 通过E步骤和M步骤使得期望最大化，是自收敛的分类算法，既不需要事先设定类别也不需要数据见的两两比较合并等操作.
缺点:
1、迭代速度慢，次数多；
2、对初始化敏感；
3、当所要优化的函数不是凸函数时，容易陷入局部最优；
4、EM可能收敛到参数空间的边界.
#####################R语言：给定一组数据设置参数########################
###EM算法在高斯混合模型GMM(Gaussian Mixture Model )中有很重要的用途.
###简单来讲GMM就是一些高斯分布的组合.如果我们已知观测到的数据的类别，
###则可以根据ML来估计出GMM的参数.反之，对于没有类别信息一堆数据，如果
###我们已知GMM的参数，可以很容易用贝叶斯公式将它们归入不同的类中；但尴尬
###的问题是我们即不知道GMM参数，也不知道观测数据的类别.以下面生成的一维数据为###例，
###我们希望找到这两个高斯分布的参数，同时为这些数据分类.
# 设置模拟参数
if(FALSE){
miu1 <- 3
miu2 <- -2
sigma1 <- 1
sigma2 <- 2
alpha1 <- 0.4
alpha2 <- 0.6
# 生成两种高斯分布的样本
n <- 5000
x <- rep(0,n)
n1 <- floor(n*alpha1)
n2 <- n - n1
x[1:n1] <- rnorm(n1)*sigma1 + miu1
x[(n1+1):n] <- rnorm(n2)*sigma2 + miu2
hist(x,freq=F)
lines(density(x),col='red')
###下面用EM算法来估计GMM的参数.
}
x <- c(-67,-48,6,8,14,16,23,24,28,29,41,49,56,60,75)
# 设置初始值
n <- 15
m <- 2
miu <- runif(m)
sigma <- runif(m)
alpha <- c(0.5,0.5)
prob <- matrix(rep(0,n*m),ncol=m)
for (step in 1:10){
# E步骤
for (j in 1:m){
prob[,j]<- sapply(x,dnorm,miu[j],sigma[j])
}
sumprob <- rowSums(prob)
prob<- prob/sumprob
####做NAN处理
for(i in 1:n)
for(j in 1:m){
{
if(is.nan(prob[i,j])){prob[i,j] <- 0}
}
}
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-3
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')
}
####得出结果
step 1 miu 6.822826 17.40323 sigma 0.9985392 5.880087 alpha 0.08455481 0.3154452
step 2 miu 6.972619 22.93183 sigma 0.9996251 38.57418 alpha 0.1252252 0.8747748
#####
###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)