EM算法学习记录

本文主要内容来自于《统计学习方法》，主要分为以下部分：

1. 极大似然估计
2. EM算法

EM算法是用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计，每次迭代由两步组成：E步，求期望；M步，求极大。
所以，本文首先讲一下极大似然估计。
极大似然估计，其实就是已知样本的观测结果，反推最有可能产生这个结果的参数值。
对于式子P(x|θ），θ为参数，x为观测结果，如果已知参数，则是概率问题，如果已知观测结果求参数则是似然问题。
举个例子
假设掷十次硬币，如果已知这枚硬币是正面的概率是0.5，求掷出正正反反正正正正反正的概率，则是统计问题，即θ=0.5，而x=[正，正，反，正，反，反，正，反，正，正]，$P(x|θ)=0.5^{10})$
但反过来，并不知道这枚硬币正面的概率是多少，只知道这次试验掷出了正正反反正正正正反正，因此，$P(x|θ)=θ^6(1-θ)^4$，是个关于θ的函数，当θ取不同值时，会有不同的概率，当函数取最大值时，θ是多少，这就是极大似然估计。
那么EM算法又是什么呢？
最开始就说了，EM算法其实是包含隐变量的极大似然估计（此处只讲极大似然估计，不讲极大后验概率估计）
那么什么是隐变量呢，再举个掷硬币的例子（以下例子来自统计学习方法）。
假设有三枚硬币，分别记作A,B,C，这些硬币正面出现的概率分别是π，p，q。进行如下试验：先掷硬币A，根据其结果选出硬币B或C，正面选硬币B，反面选硬币C；然后掷选出的硬币，掷硬币的结果，出现正面记作1，出现反面记作0；独立地重复n次试验（这里，n=10)，预测结果如跟上面那个十次结果一样，也就是：
【1,1,0,1,0，0,1,0,1,1】
这时我们假设只能观测到最后所掷硬币的结果，不能观测掷硬币的过程，也就是说，先掷的硬币A的结果我们是不知道的。
这时，我们的概率模型如下：
$p(y|θ)=\sum_{z}P(y,z|θ)=\sum_{z}P(z|θ)P(y|z,θ)=πp^y(1-p)^{1-y}+(1-π)q^y(1-q)^{1-y}$
上述模型中，y为观测结果，即最后硬币B或硬币C是正还是反，而z是隐变量，就是中间没有观测到的硬币A的结果，而此处的θ则是(π,p,q)，是模型中的参数。
将n次的观测结果表示为$Y=(Y_1,Y_2,...,Y_n)^T$,而n次未观测数据表示为$Z=(Z_1,Z_2,...,Z_n)^T$
则似然函数为：
$P(Y|θ)=\sum_{Z}P(Z|θ)P(Y|Z,θ)=\prod_{j=1}^n[πp^{y_j}(1-p)^{1-{y_j}}+(1-π)q^{y_j}(1-q)^{1-y_j}]$
如上面例子来说，第一次观测结果为1，则有两种情况，一种是硬币A为正，选硬币B，硬币B为正，第二种为硬币A为反，选硬币C，硬币C为正，因此第一次观测结果为1（即为正面）的概率为P=硬币A为正的概率硬币B为正的概率+硬币A为反的概率硬币C为正的概率=$π*p+(1-π)q$，将观测结果1,0代入，则概率=$πp^y(1-p)^{1-y}+(1-π)q^y(1-q)^{1-y}$，即y为1时取p，y为0时取1-p,q同理。
每一次结果各自独立，所以把每一次的概率都相乘，则就是上处似然函数，其中n为10。
上述问题，就是要对P（Y|θ）进行极大似然估计。
EM算法，就是可以用于求解这个问题的一种迭代算法。
针对上面例子，应用EM算法。
首先，对参数取一个初始值，假设$π^0=p^0=q^0=0.5$
那么，不管观测结果是1还是0，都有观测结果$y_j$来自于硬币B的概率为$μ_j^1=（0.5*0.5）/（0.5*0.5+0.5*0.5）=0.5$，来自硬币C的概率为$1-μ_j^1=0.5$。
由观测结果，重新求参数，则硬币A为正的概率$π^1=0.5*10/10=0.5$
硬币B为正的概率$p^1=(0.5*6)/(0.5*10)=0.6$
硬币C为正的概率$q^1=(0.5*6)/(0.5*10)=0.6$
继续，则当观测结果是1时，此时观测结果来自硬币B的概率为$μ_j^2=（0.5*0.6）/（0.5*0.6+0.5*0.6）=0.5$，来自硬币C的概率为$1-μ_j^2=0.5$，当观测结果是0时，$μ_j^2=（0.5*0.4）/（0.5*0.4+0.5*0.4）=0.5$，来自硬币C的概率为$1-μ_j^2=0.5$。
继续由观测结果重新求参数，则硬币A为正的概率$π^2=0.5*10/10=0.5$
硬币B为正的概率$p^2=(0.5*6)/(0.5*10)=0.6$
硬币C为正的概率$q^2=(0.5*6)/(0.5*10)=0.6$
收敛，因此，模型参数θ的极大似然估计就是（0.5,0.6,0.6），即当π=0.5，p=0.6，q=0.6时，出现该观测结果的可能性最大。
现在我们来换一下初始参数，假设$π^0=0.4，p^0=0.6，q^0=0.7$
那么，第一步，当观测结果为1时，来自硬币B的概率为，$μ_j^1=（0.4*0.6）/（0.4*0.6+0.6*0.7）=0.3636$，来自硬币C的概率为，$1-μ_j^1=0.6364$；当观测结果为0时，来自硬币B的概率为，$μ_j^1=（0.4*0.4）/（0.4*0.4+0.6*0.3）=0.4706$，来自硬币C的概率为0.5294.
在此基础上，由观测结果重新计算参数，则硬币A为正的概率（即观测结果来自硬币B）为$π^1=（0.3636*6+0.4706*4）/10=0.4064$
硬币B为正的概率$p^1=(0.3636*6)/(0.3636*6+0.4706*4)=0.5368$
硬币C为正的概率$q^1=(0.6364*6)/(0.6364*6+0.5294*4)=0.6433$
继续计算，显然会与之前的最终解不同。
上述过程就是EM算法，其中在已知参数模型下，求观测结果来自于硬币B的概率为E步，是在计算不同隐变量的期望。
而在此基础上，重新估算参数值的过程就是M步，意在求极大化。
不停迭代，得到的最后收敛的值便是EM算法的最终解。
由上述可知，EM算法的最后结果与初始值有较大的关系。
对于上述三硬币模型来说。
E步：计算在模型参数$π^i,p^i,q^i$下观测数据来自硬币B的概率：
$μ_j^{i+1}=\frac{π^i(p^i)^{y_j}(1-p^i)^{1-y_j}}{π^i(p^i)^{y_j}(1-p^i)^{1-y_j}+(1-π)^i(q^i)^{y_j}(1-q^i)^{1-y_j}}$
M步：计算模型参数的新估计值
$π^{i+1}=\frac{\sum_{j=1}^nμ_j^{i+1}}{n}$

$p^{i+1}=\frac{\sum_{j=1}^nμ_j^{i+1}y_j}{\sum_{j=1}^nμ_j^{i+1}}$

$q^{i+1}=\frac{\sum_{j=1}^n(1-μ_j^{i+1})y_j}{\sum_{j=1}^n(1-μ_j^{i+1})}$
综合以上的例子，引出EM算法。
输入：观测变量数据Y，隐变量数据Z，联合分布P(Y,Z|θ),条件分布P(Z|Y,θ）；
输出：模型参数θ
（1）选择参数的初值$θ^{(0)}$，开始迭代；
（2）E步：记$θ^{(i)}$为第i次迭代参数θ的估计值，在第i+1次迭代的E步，计算
$Q(θ,θ^{(i)})=E_z[logP(Y,Z|θ)|Y,θ^{(i)}]=\sum_zlogP(Y,Z|θ)P(Z|Y,θ^{(i)})$
这里，$P(Z|Y,θ^{(i)})$是在给定观测数据Y和当前的参数估计$θ^{(i)}$下隐变量数据Z的条件概率分布，$logP(Y,Z|θ)$是完全数据的对数似然函数；
（3）M步：求使$Q(θ,θ^{(i)})$极大化的θ，确定第i+1次迭代的参数的估计值$θ^{(i+1)}$
$θ^{(i+1)}=argmax_θQ(θ,θ^{(i)})$
（4）重复（2）（3）步，直至收敛
其中，$Q(θ,θ^{(i)})$是EM算法的核心，称为Q函数，其表示完全数据(Y,Z)的对数似然函数$logP(Y,Z|θ)$关于在给定观测数据Y和当前参数$θ^{(i)}$下对未观测数据Z的条件概率分布$P(Z|Y,θ^{(i)})$的期望。