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)