EM演算法原理以及高斯混合模型實踐

標籤：技術   psi   ilo   重要   amp   strong   ima   ati   組合   

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)

EM演算法原理以及高斯混合模型實踐