EM 算法

MLE 与 EM算法在参数估计里应用真是很多， PLSA就是用 EM 来求解的 ，估计这些都是概率图模型中会涉及到的，以后有机会再去系统的学习下概率图模型。

Maximum Likelihood Estimate

极大似然估计（MLE）是给定数据集后用来求解模型参数的方法，其问题形式是这样的，给定来自随机变量 $X$ 的观测数据集合 $left {  x_i ight }_{i=1}^N$ , 随机变量 $X$ 的概率密度函数 $f(x| heta)$ ,其中 $ heta$ 是为概率密度的未知参数，现在即可根据 MLE 求最优的参数 $ heta^*$ ，步骤如下：

1） 根据密度函数，列出观测数据的联合分布：

[f(x_1,x_2,...,x_N| heta) = f(x_1| heta)f(x_2| heta)...f(x_N| heta) = prod_{i=1}^{N}f(x_i| heta)]

2）求使得联合分布取值最大的参数 $ heta^*$，可先两边同时取 log ，然后求 log 函数的极大值即可，即：

[ heta^* = argmax_{ heta}L( heta)= argmax_{ heta} logprod_{i=1}^{N}f(x_i| heta) =  argmax_{ heta}sum_{i=1}^{N}log f(x_i| heta) ]

显然对 $L( heta)$ 求导，另导数得 0 ，求得的解即为最优的 $ heta^*$了， 而且对 MLE 来说数据量越多，所得到的模型会越能反映数据的真实分布。

Expectation Maximization Algorithm

期望最大化算法（Expectation Maximization Algorithm，EM）用于含有隐变量概率模型的 MLE 或 MAP，所谓的隐变量就是每个可见随机变量的值都对应着一个隐藏的随机变量，比如说 PLSA 里的主题或者 GMM 里的类别均为隐变量，用 $X$ 表示可见变量集合，第 $n$ 行 $x_n^T$ 用于表示第 $n$ 个观测数据， 对应的隐变量为 $z_n^T$，把所有的隐变量一起记作 $Z$ ， 整个概率模型的参数记做 $ heta$ ，这时 log 似然函数的形式如下：

[L( heta) = log P(X| heta) = log left { sum_Z P(X,Z| heta) ight }]

对数似然中的 $log left { sum ight }$ 部分涉及到 log里边求和的形式，很难直接计算结果，所以一般采用 EM 算法求解。EM 通过迭代逐步极大化 $L( heta)$ ，假设第 $i$ 次迭代后 $ heta$ 的估计值是 $ heta^{(i)}$ ， $ heta^{(i)}$ 已知后，下一次迭代需要使得 $L( heta)$ 更大，迭代前后的差可表示为：

[L( heta) - L( heta^{(i)}) =log left[sum_Z  P(X|Z, heta)P(Z| heta) ight ] - log P(X| heta^{(i)})  ]

计算过程中，会涉及到 $log left { sum ight }$  ，解决方法便是使用 Jensen 不等式来去掉这种难以计算的形式，其形式如下：

Jensen's Inequality

当 $f$ 为凹函数且 $sum_i lambda_i  =1 , lambda_i ge 0$ 时，有 $f(sum_ilambda_ix_i) ge sum_i lambda_i  f(x_i)$.

现在来计算迭代前后的差值：

egin{aligned}
L( heta)-L( heta^{(i)}) &= log left [sum_Z P(X|Z, heta)P(Z| heta) ight] - log P(X| heta^{(i)}) \&= log left [sum_Z P(X|Z, heta)P(Z| heta) cdot frac{P(Z|X, heta^{(i)})}{P(Z|X, heta^{(i)}}) ight] - log P(X| heta^{(i)}) \
&= logleft[sum_Z P(Z|X, heta^{(i)})frac{P(X|Z, heta)P(Z| heta)}{P(Z|X, heta^{(i)})} ight ]- log P(X| heta^{(i)}) \
& mathbf{ Using Jensen's Inequality. }\
&ge sum_Z P(Z|X, heta^{(i)}) logleft [frac{P(X|Z, heta)P(Z| heta)}{P(Z|X, heta^{(i)})} ight ]- log P(X| heta^{(i)}) \
&= sum_Z P(Z|X, heta^{(i)}) log left [frac{P(X|Z, heta)P(Z| heta)}{P(Z|X, heta^{(i)})P(X| heta^{(i)})} ight ]
end{aligned}

因为 $log$ 函数为凹函数,所以可以直接使用 Jensen's Inequality ，使用之后，还能得到似然函数的下界，令：

egin{aligned}
B( heta , heta^{(i)})=L( heta^{(i)}) +sum_ZP(Z|X, heta^{(i)})logleft [frac{P(X|Z, heta)P(Z| heta)}{P(Z|X, heta^{(i)})P(X| heta^{(i)})} ight ]
end{aligned}

于是可得下界为 ：

[L( heta) ge B( heta, heta^{(i)})]

若下次迭代还取 $ heta^{(i)}$ ，会导致 $B( heta^{(i)} , heta^{(i)})  = L( heta^{(i)})$ ：

egin{aligned}
B( heta^{(i)} , heta^{(i)})
&= L( heta^{(i)}) +sum_ZP(Z|X, heta^{(i)})logleft [frac{P(X|Z, heta^{(i)})P(Z| heta^{(i)})}{P(Z|X, heta^{(i)})P(X| heta^{(i)})} ight ]\
&= L( heta^{(i)}) +sum_ZP(Z|X, heta^{(i)})logleft [frac{P(X,Z | heta^{(i)})}{P(X,Z | heta^{(i)})} ight ]\
&= L( heta^{(i)}) +sum_ZP(Z|X, heta^{(i)})log 1\
&=L( heta^{(i)})
end{aligned}

下界 $B( heta, heta^{(i)})$ 是一个关于 $ heta$ 的函数，为了使 $L( heta)$ 尽可能增加，需要最大化下界，就是说找到满足条件的 $ heta$ 使 $B( heta, heta^{(i)}$) 达到极大，记下一次迭代的值为 $ heta^{(i+1)}$:

egin{aligned}
heta^{(i+1)} &= argmax_{ heta}B( heta, heta^{(i)})\
&=argmax_{ heta} left {L( heta^{(i)}) +sum_ZP(Z|X, heta^{(i)}) log frac{P(X|Z, heta)P(Z| heta)}{P(Z|X, heta^{(i)})P(X| heta^{(i)})} ight } \
& mathbf{ Now drop terms which are constant w.r.t. heta} \
&=argmax_{ heta} left {sum_ZP(Z|X , heta^{(i)})log P(X|Z, heta)P(Z| heta) ight } \
&=argmax_{ heta} left {sum_ZP(Z|X, heta^{(i)}) log frac{ P(X,Z, heta)}{P(Z, heta)}frac{ P(Z, heta)}{P( heta)} ight } \
&=argmax_{ heta} left {sum_ZP(Z|X, heta^{(i)})log P(X,Z| heta) ight} \
&=argmax_{ heta} Q( heta, heta^{(i)}) \
end{aligned}

这里 $P(Z|X, heta^{(i)})$ 是给观测数据 $X$ 和当前参数 $ heta^{(i)}$ 时 $Z$ 的后验概率分布， $Q$ 函数就是完全数据的对数似然函数 $log P(X,Z| heta)$ 在给定下关于隐变量 $Z$ 分布的期望。

[Q( heta, heta^{(i)})=E_{Z|X, heta^{(i)}}left[log P(X,Z| heta) ight]=sum_Z log P(X,Z| heta) P(Z|X, heta^{(i)})]

综上给出 EM 算法：

输入：观测数据 $X$，隐变量数据 $Z$ ，联合分布 $P(X,Z| heta)$；

输出：极大似然参数 $ heta$.

（1）选择初始参数 $ heta^{(0)}$ ；

（2）E Step：计算隐变量 $Z$ 在参数 $ heta^{(i)}$ 下的后验分布：[P(Z|X, heta^{(i)})]

（3）M Step：估计 $ heta^{(i+1)}$ 的值：

[ heta^{(i+1)} = arg max_{ heta} Q( heta, heta^{(i)})]（4）重复（2）至（3），直到收敛.

关于 EM 算法有几点需要注意：

• 关于 $P(Z|X, heta)$ 的计算可采用贝叶斯公式计算；
• 有的 EM 版本写法不一样，比如统计学习方法中在 E 步计算的是 $Q( heta , heta^{(i)})$ ，本书同 PRML 的写法；
• 初值 $ heta^{(0)}$ 可以任意选择 ，不同初值，结果差距很大，因为 EM 是对初值敏感的；
• 这里的 EM 是 MLE 版本，如果对参数 $ heta$ 引入一个先验 $P( heta)$，就变成了 MAP 版本，计算时 E 步不变，M 步极大化  $Q( heta, heta^{(i)})+P( heta)$ 即可；
• 迭代停止的条件，给定较小的正数 $varepsilon _1$ $varepsilon _2$， 迭代的停止条件为以下两者任一个均可：

[ || heta^{(i+1)}- heta^{(i)}||< varepsilon_1 , || Q( heta^{(i+1)}, heta^{(i)})-Q( heta^{(i)}, heta^{(i)})||<varepsilon _2 ]

下图形象的展示了 EM 算法的执行过程，上方曲线为 $L( heta)$ ,下方为 $B( heta, heta^{(i)})$ ，刚才分析已经得到 $B( heta, heta^{(i)})$ 是 $L( heta)$ 的下界，两个函数在点 $ heta = heta^{(i)}$ 处相等， EM 可以寻找下一个点 $ heta^{(i+1)}$ 使得 $B( heta, heta^{(i)})$ 与 $Q( heta, heta^{(i)})$ 极大化，因为 $B( heta, heta^{(i)})$ 的增加会导致 $L( heta)$ 也是增加的，接下来在方才寻找的 $ heta^{(i+1)}$ 的基础上继续进行下一次迭代，整个迭代过程似然函数 不断增大，从而达到估计参数的目的，但需要注意 EM 找到的并非全部最优解。为了得到更好的解，EM 通常可以选取不同的初值点 $ heta^{(0)}$ ,迭代完成后，比较各个结果选出最优的一个即可。

EM 算法是收敛的，设 $P(X| heta)$ 为观测数据的似然函数， $ heta^{(i)}$ 为 EM 算法得到的参数估计序列，$P(X| heta^{(i)})$  为对应的似然函数序列，则有：

[P(X| heta^{(i+1)}) ge P(X| heta^{(i)}) ]

证明从略，参考统计学习方法即可！最后来看关于 EM 算法的示例，也就是三硬币模型，三个硬币分别记做 A,B,C ,这些硬币正面朝上的概率分别为 $pi,p,q$ 即混合模型的参数 $ heta = (pi,p,q)$ ，实验是这样的，首先选择抛硬币 A ，若 A 为正，则选择 B ，否则选择 C ，然后抛出所选硬币进行观测，正面记做 1 ，反面记做 0 ，重复 10 次得到如下的结果：

根据以上数据，直接进行三个二项分布的 MLE 即可得到 $ heta$，如果没有给出可见隐变量 $Z$ ，这时 MLE 便无法求解了，采用 EM 算法即可。

数据集中每个 $x_i$ 对应着一个隐变量 $z_i$ ，这里隐变量的取值为 $Zin(B,C)$， 这时 $X$ 的联合分布为：

[P(X| heta) = sum_ZP(Z| heta)P(X|Z, heta) = prod_{i=1}^N left{sum_{j=1}^M P(z_j| heta)P(x_i|z_j, heta) ight }]

这里有 $N = 10 ,M = 2$ ，令 $P(Z = B) = pi, P(Z = C) = 1 - pi$ ，则完全数据的似然函数可写作：

[L( heta) =log P(X,Z| heta) = sum_{i=1}^N log left { pi p^{x_i}(1-p)^{1-x_i} + (1-pi)q^{x_i}(1-q)^{1-x_i} ight }]

接下来便是 EM 算法，首先初始化参数:

[ heta^{old} = left { pi^{old},p^{old},q^{old}  ight }]

E 步：根据 $ heta^{old}$,使用贝叶斯公式，计算 $Z$ 的后验概率分布：

[P(Z = B|X, heta^{old})
=  frac{P(Z=B)P(X|Z=B, heta^{old})}{P(Z=B)P(X|Z=B, heta^{old})+P(Z=C)P(X|Z=C, heta^{old})} ]

对于样本 $x_i$，将其隐变量记做 $mu_i = P(z_i = B|x_i, heta^{old})$，则其隐变量的后验概率为：

[mu_i=frac{pi^{old}(p^{old})^{y_i}(1-p^{old})^{1-y_i}}{pi^{old}(p^{old})^{y_i}(1-p^{old})^{1-y_i} + (1-pi^{old})(q^{old})^{y_i}(1-q^{old})^{1-y_i}} ]

M 步：极大化 $Q( heta, heta^{old})$ 即可，这里 Q 函数的形式如下：

egin{aligned}
Q( heta, heta^{old})
&= sum_Z P(Z|X, heta^{old})log P(X,Z| heta)\
&=sum_{i=1}^N P(z_i = B|x_i, heta^{old})log P(x_i, z_i=B| heta)+P(z_i = C|x_i, heta^{old})log P(x_i, z_i=C| heta)\
& mathbf{ Import  Sign } mu_i  =P(z_i = B|x_i, heta^{old})  \
&=sum_{i=1}^N mu_i log left{pi p^{y_i}(1-p)^{1-y_i} ight } + (1-mu_i)log left{(1-pi) q^{y_i} (1-q)^{1-y_i} ight}
end{aligned}

分别对 $pi,p,q$ 求导，即可得到其极大值，比如对于 $pi$ 求导 ，令导数得 0 ，即可得到 $pi^{new}$：
[frac{partial Q}{partial pi} = sum_{i=1}^Nleft (frac{mu_i}{pi}  -frac{1-mu_i}{1-pi} ight ) = frac{sum_{i=1}^N mu_i -Npi}{pi(1-pi)} = 0 Rightarrow pi^{new} = frac{1}{N}sum_{i=1}^N mu_i]

计算得到参数值 $ heta^{new}$ 分别为：

egin{aligned}
pi^{new}&= frac{1}{N}sum_{i=1}^N mu_i\
p^{new}&= frac{sum_{i=1}^N mu_ix_i}{sum_{i=1}^N mu_i}\
q^{new}&= frac{sum_{i=1}^N (1-mu_i)x_i}{sum_{i=1}^N (1-mu_i)}
end{aligned}

不断迭代这个过程直到收敛便可得到模型参数，这是一个完整的 EM 算法的例子，算法非常经典，并且适用于一系列带有隐变量的混合模型求解，典型的有 GMM,HMM,PLSA 等。

 

参考文献:

1. 书： 统计学习方法 | PRML

2. http://blog.tomtung.com/2011/10/em-algorithm/#id8

3. http://www.thebigdata.cn/upload/2015-02/15020913247095.pdf

4. http://alexkong.net/2014/07/GMM-and-EM/

5. http://freemind.pluskid.org/machine-learning/regularized-gaussian-covariance-estimation/

6. https://www.cs.utah.edu/~piyush/teaching/EM_algorithm.pdf(主要参考)

7. http://chenrudan.github.io/blog/2015/12/02/emexample.html