利用python对混合高斯分布的研究

EM算法是一种迭代算法，用于含有隐变量的概率模型的极大似然估计，或极大后验概率估计。EM算法有两步，E步和M步，E步求期望，M步求最大化。合称EM算法，本文主要以the element of statistic learning 第八章关于EM算法的介绍过程为主线撰写，对难以理解的地方尽量加以解释。

    先来看个例子。有一堆数据，如下表。

在R中，我们直接把它写到list中。代码:y=c(-0.39, 0.12,0.94 ,1.67, 1.76, 2.44 ,3.72, 4.28, 4.92 ,5.53,0.06 ,0.48, 1.01, 1.68, 1.80 ,3.25, 4.12, 4.60, 5.28 ,6.22 )

hist(y,breaks=15,col="red")#画出书上图8.5所示的直方图。

从这个图中可以看出，这堆数据有两个峰值，用单个正态分布难以描述，所以可以用两个正态分布的混合来建模。令

其中 Δ 属于（0，1）。这个生成过程解释如下：

先以概率π生成一个参数Δ∈（0，1），然后把这个Δ代入Y中，则Y不是等于Y1 就是等于Y2。现在我们想用Y分布函数来建模上面这组数据。其中要确定的数据为

基于N个训练数据的似然函数是

直接求这个似然函数比较难，因为其中log函数里面有一个加号。这里有一个简单的策略可以避开在log里面求和。因为Δ不是取0就是取1.当Δ取0时，样本来自Y1，反之，当Δ取1时，样本来自Y2.假定我们对于每个样本，已经知道Δ的值了。于是对数似然函数可以写成

其中的π，mu和sigma都可以根据各自的样本类分别估算出。π可根据Δ为1的比例得出。但现在的问题是我们并不知道Δ的值。因此，我们采用迭代的方法，每次采用估算出来的Δ的期望值。

这个期望值是根据theta和已知样本计算出来的，theta一开始赋予一个初值，后面在迭代中持续更新。具体的赋值方法见下面的算法流程。注意对于每一个样本n，我们都要计算一个γ，对于k类问题，每一类都要计算一个γ。所以有NK个γ。这里因为是两类问题，所以只要计算一个，因为另一个只要1-γ即可。伽马可以看做是Δ的期望值。

     在EM算法中，我们已有的输入是：不完全数据Y，隐变量Z，隐变量和Y的联合分布-这里是混合高斯分布。隐变量自身的分布，这里是伯努利分布（取0，1）。

EM算法就是希望能够根据隐变量的后验分布（根据给定初值和样本计算）来计算完全数据（包括样本和隐变量）的期望值。并且使这个期望最大。

下面贴出该算法的R代码：以及相应的计算结果。

mu1 <-y[sample(1:20,1)];mu2 <- y[sample(1:20,1)];pi <- 0.5
 sigma1 <- sd(y);sigma2 <-sd(y);
 loglikelihood <- c()

 for (i in 1:50)
     {
         gama <- c()#E-step
         for(j in 1:length(y))
             {
               fhtheta1 <- dnorm(y[j],mean=mu1,sd=sigma1)
               fhtheta2 <- dnorm(y[j],mean=mu2,sd=sigma2)
               gamaTemp <- pi*fhtheta2/((1-pi)*fhtheta1+pi*fhtheta2)
               gama <- c(gama,gamaTemp)
             }
         gamasum <- sum(gama) #M-step
         gamasumY <- sum(gama*y)
         gamasuminverse <- sum(1-gama)
         gamasuminverseY <- sum((1-gama)*y)
         mu1 <- gamasuminverseY/gamasuminverse
         mu2 <- gamasumY/gamasum
         sigma1temp <- sum((1-gama)*(y-mu1)^2)
         sigma1 <- sqrt(sigma1temp/gamasuminverse)
         sigma2temp <- sum(gama*(y-mu2)^2)
         sigma2 <- sqrt(sigma2temp/gamasum)
           pi <- gamasum/length(y)
         loglikelihoodTemp <- sum(log((1-pi)*dnorm(y,mean=mu1,sd=sigma1)+pi*dnorm(y,mean=mu2,sd=sigma2)))
         loglikelihood <- c(loglikelihood,loglikelihoodTemp)}

输出结果为

pi;mu1;mu2;sigma1;sigma2
0.5545902
 [1] 4.655913
 [1] 1.083162
 [1] 0.9048721
 [1] 0.9007611

和书上的结果略有不同，可能是由于初始值的原因。不过大致能够说明情况。下一篇博客介绍一般性的EM算法。

plot(loglikelihood[1:50],col="green")#画出书上的图8.6

大致趋势是一样的。