EM算法
EM算法，即期望极大算法，是一种迭代算法，每一次迭代包含两步：
- E步：求期望；
- M步：求极大。
算法核心内容
输入：观测变量$Y$，隐变量$Z$，联合分布$P(Y,Z|\theta)$，条件分布$P(Z|Y,\theta)$；
输出：模型参数$\theta$。
算法步骤如下：
选择参数的初值$\theta^{(0)}$，开始迭代；
迭代，分为两步，即E步和M步，
E步：记$\theta^{(i)}$为第$i$次迭代参数$\theta$的估计值，则下一次迭代的E步，计算
M步：计算最优参数
重复第2步，直至算法收敛。
具体实例
这个例子来自《统计学习方法》，只不过书中未给出$Q(\theta,\theta^{(i)})$函数具体的计算过程。
（三硬币模型）假设有3枚硬币，分别记作A，B，C，这些硬币正面出现的概率分别是$\pi,p,q$。下面，进行如下掷硬币试验：先掷硬币A，根据其结果选出硬币B或者C：正面则选硬币B，反面选硬币C；然后掷选出的硬币，将这次掷硬币的结果，出现正面记作1，出现反面记作0；独立地重复n次试验，问如何估计三硬币正面出现的概率，即$\pi,p,q$三个参数。具体地，n=10，每次掷硬币结果为
1,1,0,1,0,0,1,0,1,1
下面使用EM算法求解上述问题。
E步：计算期望
设$Y$为观测变量，即最终掷硬币的结果；$Z$为隐变量，即掷硬币A的结果；参数$\theta=(\pi,p,q)$.
则有
按理说，我们需要把上述$Q(\theta,\theta^{(i)})$的每一项计算出来，然后代入得到完整的表达式，才能进行接下来的M步；但是可以发现，在接下来的M步，即优化求解参数$\displaystyle\theta^{(i+1)} = \arg\max_\theta Q(\theta,\theta^{(i)})$中，我们只需要关注与$\theta$有关的项，其余的项便是常数项，可以忽略。因此，设
故
因此，有
于是，只需计算关键的两项，即$\log P(y_j,Z=0|\theta)$和$\log P(y_j,Z=1|\theta)$.
这个概率很简单求得，只需列出分布律即可计算出来。
| $Z$ | 0 | 1 |
|---|---|---|
| $P$ | $1-\pi$ | $\pi$ |
| $Y | Z=0$ | 0 | 1 |
|---|---|---|---|
| $P$ | $1-q$ | $q$ |
| $Y | Z=1$ | 0 | 1 |
|---|---|---|---|
| $P$ | $1-p$ | $p$ |
而
因此，求得
其中，
M步：计算极大
接下来，进行M步，即求解$\displaystyle \theta^{(i+1)} = \arg\max_\theta Q(\theta,\theta^{(i)})$.对各个参数求偏导：
分别令上式等于零，得到
代码实现
基于上述推导结果，下面编写Python代码实现EM算法的迭代计算。
import numpy as np
# 已知观测变量取值
Y = [1, 1, 0, 1, 0, 0, 1, 0, 1, 1]
Y_array = np.array(Y)
# 设置参数初值
pi, p, q = 0.4, 0.6, 0.7定义计算$\mu_j^{(i)}$和对数似然$LL\left(\theta^{(i)}\right)$的函数：
def mu(pi, p, q, y):
'''
计算μ
'''
cb = pi*p**y*(1-p)**(1-y) # 来自硬币B的概率
cc = (1-pi)*q**y*(1-q)**(1-y) # 来自硬币C的概率
return cb / (cb + cc)
def LL(pi, p, q, Y):
'''
计算对数似然
'''
return sum([np.log(q**y*(1-q)**(1-y)*(1-pi)+p**y*(1-p)**(1-y)*pi) for y in Y])开始运行EM算法：
iter_max = 10
LL_value = LL(pi, p, q, Y)
for iter in range(iter_max):
mu_array = np.array([mu(pi, p, q, y) for y in Y])
# 计算当前次迭代的更新参数
pi_temp = mu_array.mean()
p_temp = (mu_array * Y_array).sum() / mu_array.sum()
q_temp = ((1- mu_array) * Y_array).sum() / (1 - mu_array).sum()
# 计算新的对数似然
LL_value_temp = LL(pi_temp, p_temp, q_temp, Y)
print("[{}/{}] pi: {:.4f} -> {:.4f}, p: {:.4f} -> {:.4f}, q: {:.4f} -> {:.4f}, LL value: {:.4f} -> {:.4f}".format(iter+1, iter_max, pi, pi_temp, p, p_temp, q, q_temp, LL_value, LL_value_temp))
pi, p, q, LL_value = pi_temp, p_temp, q_temp, LL_value_temp运行结果如下。
[1/10] pi: 0.4000 -> 0.4064, p: 0.6000 -> 0.5368, q: 0.7000 -> 0.6432, LL value: -6.8083 -> -6.7301
[2/10] pi: 0.4064 -> 0.4064, p: 0.5368 -> 0.5368, q: 0.6432 -> 0.6432, LL value: -6.7301 -> -6.7301
[3/10] pi: 0.4064 -> 0.4064, p: 0.5368 -> 0.5368, q: 0.6432 -> 0.6432, LL value: -6.7301 -> -6.7301
[4/10] pi: 0.4064 -> 0.4064, p: 0.5368 -> 0.5368, q: 0.6432 -> 0.6432, LL value: -6.7301 -> -6.7301
[5/10] pi: 0.4064 -> 0.4064, p: 0.5368 -> 0.5368, q: 0.6432 -> 0.6432, LL value: -6.7301 -> -6.7301
[6/10] pi: 0.4064 -> 0.4064, p: 0.5368 -> 0.5368, q: 0.6432 -> 0.6432, LL value: -6.7301 -> -6.7301
[7/10] pi: 0.4064 -> 0.4064, p: 0.5368 -> 0.5368, q: 0.6432 -> 0.6432, LL value: -6.7301 -> -6.7301
[8/10] pi: 0.4064 -> 0.4064, p: 0.5368 -> 0.5368, q: 0.6432 -> 0.6432, LL value: -6.7301 -> -6.7301
[9/10] pi: 0.4064 -> 0.4064, p: 0.5368 -> 0.5368, q: 0.6432 -> 0.6432, LL value: -6.7301 -> -6.7301
[10/10] pi: 0.4064 -> 0.4064, p: 0.5368 -> 0.5368, q: 0.6432 -> 0.6432, LL value: -6.7301 -> -6.7301可以发现，根据EM算法，待估计的参数$\pi,p,q$一步即收敛，收敛到$\pi^{(1)}=0.4064,p^{(1)}=0.5368,q^{(1)}=0.6432$.
算法的证明
我们的目标是要最大化对数似然函数（极大似然估计）
因为涉及隐变量$Z$，故对对数似然函数做如下变形
接下来，我们来计算一下对数似然和当前迭代步参数下的似然间的差值：
首先是一个小trick，引入因子$P(Z|Y,\theta^{(i)})$，它将成为整个证明过程的关键：
接下来，由琴生不等式，对上式进行放缩，得到：
另外，注意到$\displaystyle \sum_Z P(Z|Y,\theta^{(i)})=1$，得：
更明显地，记
因此，有
故$B(\theta, \theta^{(i)})$是对数似然函数$LL(\theta)$的一个下界。那么，如果我们每次迭代时，求出使得该下界$B(\theta, \theta^{(i)})$最大的参数$\theta^*$作为最新参数，那么对数似然函数$LL(\theta)$的下界就会随着迭代次数增加而不断增大，进而不断地提高对数似然$LL(\theta)$.那么我们求解
接下来，我们只需要关注$B(\theta,\theta^{(i)})$中与参数$\theta$有关的部分，得到
证毕。
