首页 > 代码库 > EM算法原理以及高斯混合模型实践

EM算法原理以及高斯混合模型实践

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=http://www.mamicode.com/iris[,1:4], what="classification",dimens=c(3,4))

table(iris$Species, mc$classification)

EM算法原理以及高斯混合模型实践