ESL-chapter8-EM算法介绍1-混合高斯的例子

最新推荐文章于 2020-12-02 11:24:39 发布

果然好吃

最新推荐文章于 2020-12-02 11:24:39 发布

阅读量832

点赞数

本文链接：https://blog.csdn.net/u010198460/article/details/38117565

版权

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)