EM算法

概率模型有时既含有观测变量（observable vriable），又含有隐变量或潜在变量（latent variable）。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数。

EM算法[Dempster, 1977]是一种迭代算法，用于含有隐变量（hidden variable）的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望（expectation）；M步，求极大（maximization）。所以这一算法称为期望极大算法（expectation maximization algorithm，EM）

1 EM 算法

import numpy as np

def e_step(y, pi, p, q):
    
    mu_1 = pi * p ** y * (1 - p) ** (1 - y)
    mu_2 = (1 - pi) * q ** y * (1 - q) ** (1 - y)
    
    mu = mu_1 / (mu_1 + mu_2)
    
    return mu

def m_step(y, mu):
    
    n = len(y)
    pi = np.sum(mu) / n
    p = sum(y * mu) / sum(mu)
    q = sum(y * (1 - mu)) / sum(1 - mu)
    
    return pi, p, q

def diff(pi, p, q, pi_, p_, q_):
    
    return np.sum(np.abs([pi - pi_, p - p_, q - q_]))

def em(y, pi, p, q):
    cnt = 1
    while True:

        print("-" * 10)
        print("iter %d:" % cnt)
        pi_ = pi
        p_ = p
        q_ = q

        mu = e_step(y, pi, p, q)
        print(mu)
        pi, p, q = m_step(y, mu)
        print(pi, p, q)

        if diff(pi, p, q, pi_, p_, q_) < 0.001:
            break

        cnt += 1
        
    return pi, p, q

y = np.array([1, 1, 0, 1, 0, 0, 1, 0, 1, 1])

print("*" * 10)
pi = 0.5
p = 0.5
q = 0.5

pi, p, q = em(y, pi, p, q)

print("*" * 10)
pi = 0.4
p = 0.6
q = 0.7

pi, p, q = em(y, pi, p, q)

print("*" * 10)
pi = 0.46
p = 0.55
q = 0.67

pi, p, q = em(y, pi, p, q)

**********
----------
iter 1:
[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5]
0.5 0.6 0.6
----------
iter 2:
[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5]
0.5 0.6 0.6
**********
----------
iter 1:
[0.36363636 0.36363636 0.47058824 0.36363636 0.47058824 0.47058824
 0.36363636 0.47058824 0.36363636 0.36363636]
0.40641711229946526 0.5368421052631579 0.6432432432432431
----------
iter 2:
[0.36363636 0.36363636 0.47058824 0.36363636 0.47058824 0.47058824
 0.36363636 0.47058824 0.36363636 0.36363636]
0.40641711229946526 0.5368421052631579 0.6432432432432431
**********
----------
iter 1:
[0.41151594 0.41151594 0.53738318 0.41151594 0.53738318 0.53738318
 0.41151594 0.53738318 0.41151594 0.41151594]
0.461862835113919 0.5345950037850112 0.6561346417857326
----------
iter 2:
[0.41151594 0.41151594 0.53738318 0.41151594 0.53738318 0.53738318
 0.41151594 0.53738318 0.41151594 0.41151594]
0.46186283511391907 0.5345950037850112 0.6561346417857326

通常情况， $Y$表示观测随机变量的数据， $Z$表示隐随机变量的数据。 $Y$和 $Z$均已知称为完全数据（complete-data），仅有观测数据 $Y$称为不完全数据（incomplete-data）。假设给定观测数据 $Y$，其概率分布是 $P(Y; \theta)$，其中 $\theta$是需要估计的模型参数，那么不完全数据 $Y$的似然函数是 $P(Y; \theta)$，对数似然函数 $L(\theta) = \log P(Y; \theta)$；假设 $Y$和 $Z$的联合概率分布是 $P(Y, Z; \theta)$，那么完全数据的对数似然函数是 $L(\theta) = \log P(Y, Z; \theta)$。

EM算法通过迭代求解 $L(\theta) = \log P(Y, Z ; \theta)$的极大似然估计。每次迭代包含两步：E步，求期望；M步，求极大化。

算法9.1（EM算法）

输入：观测变量数据 $Y$，隐变量数据 $Z$，联合分布 $P(Y, Z; \theta)$，条件分布 $P(Z | Y; \theta)$；

输出：模型参数 $\theta$。

1. 选择参数初值 $\theta^{(0)}$，开始迭代；

2. E步：记 $\theta^{(i)}$为第 $i$次迭代参数 $\theta$的估计值，在第 $i + 1$次迭代的E步，计算

$\begin{aligned} Q(\theta, \theta^{(i)}) & = \text{E}_{Z} \left[ \log P(Y, Z; \theta) | Y; \theta^{(i)} \right] \\ & = \sum_{Z} P(Z | Y; \theta^{(i)}) \log P(Y, Z; \theta) \end{aligned} \tag {9}$

其中， $P(Z; Y, \theta^{(i)})$是在给定观测数据 $Y$和当前的参数估计 $\theta^{(i)}$下隐变量数据 $Z$的条件概率分布；

3. M步：求使 $Q(\theta, \theta^{(i)})$极大化的 $\theta$，确定第 $i + 1$次迭代的参数的估计值 $\theta^{(i + 1)}$

$\theta^{(i + 1)} = \argmax_{\theta} Q(\theta, \theta^{(i)}) \tag {10}$

4. 重复第2步和第3步，直到收敛。

方程（9）的函数 $Q(\theta, \theta^{(i)})$是EM算法的核心，称为 $Q$函数（ $Q$ function）。

定义9.1（Q函数） 完全数据的对数似然函数 $P(Y, Z | \theta)$是关于给定观测数据 $Y$和当前参数 $\theta^{(i)}$下，对未观测数据 $Z$条件概率分布 $P(Z | Y, \theta^{(i)})$的期望，称为 $Q$函数，即

$Q(\theta, \theta^{(i)}) = \text{E}_{Z} \left[ \log P(Y, Z; \theta) | Y; \theta^{(i)} \right] \tag {11}$

关于EM算法的几点说明：

步骤1：参数初值可以任意选择，但EM算法对初值敏感；

步骤2：E步求 $Q(\theta, \theta^{(i)})$， $Q$函数中 $Z$是未观测数据， $Y$是观测数据。 $Q(\theta, \theta^{(i)})$的第1个变元表示要极大化的参数，第2个变元表示参数的当前估计值。每次迭代实际在求 $Q$函数及其极大。

步骤3：M步极大化 $Q(\theta, \theta^{(i)})$，得到 $\theta^{(i + 1)}$，完成一次迭代 $\theta^{(i)} \rightarrow \theta^{(i + 1)}$。每次迭代使似然函数增大或达到局部极值。

步骤4：停止迭代的条件一般是对较小的正数 $\epsilon_{1}$、 $\epsilon_{2}$，若满足

$\| \theta^{(i + 1)} - \theta^{(i)} \| \lt \epsilon_{1}$

或

$\| Q(\theta^{(i + 1)}, \theta^{(i)}) - Q(\theta^{(i)}, \theta^{(i)}) \| \lt \epsilon_{2}$

则停止迭代。

推导

EM算法可通过近似求解观测数据的对数似然函数极大化问题导出：考虑一个含有隐变量的概率模型，目标是极大化观测数据（不完全数据） $Y$关于参数 $\theta$的对数似然函数，即极大化

$L(\theta) = \log P(Y; \theta) = \log \sum_{Z} P(Y, Z; \theta) = \log \sum_{Z} P(Y | Z; \theta) P(Z; \theta) \tag {12}$

假设在第 $i$次迭代后 $\theta$的估计值为 $\theta^{(i)}$。迭代求解要求 $\theta$的新估计值使 $L(\theta)$增加，即 $L(\theta) \gt L(\theta^{(i)})$，并逐步达到极大值。考虑两者差值：

$\begin{aligned} L(\theta) - L(\theta^{(i)}) = \log \sum_{Z} P(Y | Z; \theta) P(Z; \theta) - \log P(Y; \theta^{(i)}) \end{aligned}$

由Jensen不等式（Jensen inequality），其下界为：

$\begin{aligned} L(\theta) - L(\theta^{(i)}) & = \log \left( \sum_{Z} P(Z | Y; \theta^{(i)}) \frac{P(Y | Z; \theta) P(Z; \theta)}{P(Z | Y; \theta^{(i)})} \right) - \log P(Y; \theta^{(i)}) \\ & \geq \sum_{Z} P(Z | Y; \theta^{(i)})\log \left( \frac{P(Y | Z; \theta) P(Z; \theta)}{P(Z | Y; \theta^{(i)})} \right) - \log P(Y; \theta^{(i)}) \\ & = \sum_{Z} P(Z | Y; \theta^{(i)})\log \left( \frac{P(Y | Z; \theta) P(Z; \theta)}{P(Z | Y; \theta^{(i)}) P(Y; \theta^{(i)})} \right) \end{aligned}$

令

$B(\theta, \theta^{(i)}) \triangleq L(\theta^{(i)}) + \sum_{Z} P(Z | Y; \theta^{(i)})\log \left( \frac{P(Y | Z; \theta) P(Z; \theta)}{P(Z | Y; \theta^{(i)}) P(Y; \theta^{(i)})} \right) \tag {13}$

则

$L(\theta) \geq B(\theta, \theta^{(i)}) \tag {14}$

即函数 $B(\theta, \theta^{(i)})$是 $L(\theta)$的一个下界，由方程（13）可知，

$L(\theta^{(i)}) = B(\theta^{(i)}, \theta^{(i)}) \tag {15}$

为使 $L(\theta)$尽可能大的增长， $\theta^{(i + 1)}$应选择

$\theta^{(i + 1)} = \argmax_{\theta} B(\theta, \theta^{(i)}) \tag {16}$

由方程（10）、（13）和（16）可得

$\begin{aligned} \theta^{(i + 1)} & = \argmax_{\theta} B(\theta, \theta^{(i)}) \\ & = \argmax_{\theta} \left( L(\theta^{(i)}) + \sum_{Z} P(Z | Y; \theta^{(i)}) \log \frac{ P(Y | Z; \theta) P(Z; \theta) }{ P(Z | Y; \theta^{(i)}) P(Y; \theta^{(i)}) } \right) \\ & = \argmax_{\theta} \left( \sum_{Z} P(Z | Y; \theta^{(i)}) \log P(Y | Z; \theta) P(Z; \theta) \right) \\ & = \argmax_{\theta} \left( \sum_{Z} P(Z | Y; \theta^{(i)}) \log P(Y, Z; \theta) \right) \\ & = \argmax_{\theta} Q(\theta, \theta^{(i)}) \\ \end{aligned} \tag {16}$

EM算法的直观解释：图中上方曲线为 $L(\theta)$、下方曲线为 $B(\theta, \theta^{(i)})$。 $B(\theta, \theta^{(i)})$为 $L(\theta)$的下界，由方程（15）， $B(\theta, \theta^{(i)})$和 $L(\theta)$在 $\theta = \theta^{(i)}$处相等。由方程（16）、（17）可知，EM算法寻找的下一个点 $\theta^{(i + 1)}$使 $B(\theta, \theta^{(i)})$（即 $Q(\theta, \theta^{(i)})$）极大化。EM算法在 $\theta^{(i + 1)}$处重新计算函数 $Q$的值，进行下一轮迭代。在迭代过程中，对数似然函数 $L(\theta)$不断增大，但EM算法不能保证找到全局最优解。

非监督学习中的应用

EM算法可用于生成模型的非监督学习，生成模型由联合概率分布 $P(X, Y)$表示，可以认为非监督学习训练数据是联合概率分布产生的数据，其中， $X$为观测数据， $Y$为未观测数据。

2 EM算法的收敛性

EM算法提供一种近似计算含有隐变量概率模型的极大似然估计的方法，其最大优点是简单性和普适性。

定理9.1 设 $P(Y; \theta)$为观测数据的似然函数， $\theta^{(i)}$（ $i = 1, 2, \cdots$）为EM算法得到的参数估计序列， $P(Y; \theta^{(i)})$为对应的似然函数序列，则 $P(Y; \theta^{(i)})$为单调递增的，即

$P(Y; \theta^{(i + 1)}) \geq P(Y; \theta^{(i)}) \tag {18}$

证明：

由于

$P(Y; \theta) = \frac{P(Y, Z; \theta)}{P(Z | Y; \theta)}$

取对数有

$\log P(Y; \theta) = \log P(Y, Z; \theta) - \log P(Z | Y; \theta)$

由方程（11）

$Q(\theta, \theta^{(i)}) = \sum_{Z} P(Z | Y; \theta^{(i)}) \log P(Y, Z; \theta)$

令

$H(\theta, \theta^{(i)}) = \sum_{Z} P(Z | Y; \theta^{(i)}) \log P(Z | Y; \theta) \tag {19}$

则对数似然函数改写为

$\log P(Y; \theta) = Q(\theta, \theta^{(i)}) - H(\theta, \theta^{(i)}) \tag {20}$

可知

$\begin{aligned} \log P(Y; \theta^{(i + 1)}) - \log P(Y; \theta^{(i)}) & = [Q(\theta^{(i + 1)}, \theta^{(i)}) - Q(\theta^{(i)}, \theta^{(i)})] - [H(\theta^{(i + 1)}, \theta^{(i)}) - H(\theta^{(i)}, \theta^{(i)})] \end{aligned} \tag {21}$

因此，需证明方程（21）右端非负。对于第一项，由M步定义可知：

$Q(\theta^{(i + 1)}, \theta^{(i)}) - Q(\theta^{(i)}, \theta^{(i)}) \geq 0 \tag {22}$

对于第二项，由Jensen不等式：

$\begin{aligned} H(\theta^{(i + 1)}, \theta^{(i)}) - H(\theta^{(i)}, \theta^{(i)}) & = \sum_{Z} P(Z | Y; \theta^{(i)}) \log \frac{P(Z | Y; \theta^{(i + 1)})}{P(Z | Y; \theta^{(i)})} \\ & \leq \log \left( \sum_{Z} P(Z | Y; \theta^{(i)}) \frac{P(Z | Y; \theta^{(i + 1)})}{P(Z | Y; \theta^{(i)})} \right) \\ & = \log \left( \sum_{Z} P(Z | Y; \theta^{(i + 1)}) \right) \\ & = 0 \end{aligned} \tag {23}$

即

$P(Y; \theta^{(i + 1)}) \geq P(Y; \theta^{(i)})$

得证。

定理9.2 设 $L(\theta) = \log P(Y; \theta)$为观测数据的对数似然函数， $\theta^{(i)}$（ $i = 1, 2, \cdots$）为EM算法得到的参数估计序列， $L(\theta^{(i)})$为对应的对数似然函数序列，

1. 如果 $P(Y; \theta)$有上界，则 $L(\theta^{(i)}) = \log P(Y; \theta^{(i)})$收敛到某一值 $L^{\ast}$；

2. 在函数 $Q(\theta, \theta^{\prime})$与 $L(\theta)$满足一定条件下，由EM算法得到的参数估计序列 $\theta^{(i)}$的收敛值 $\theta^{\ast}$是 $L(\theta)$的稳定点。

定理9.2关于函 $Q(\theta, \theta^{\prime})$与 $L(\theta)$的条件在大多数情况下都是满足的，EM算法的收敛性包含关于对数似然函数序列 $L(\theta^{(i)})$的收敛性和关于参数估计序列 $\theta^{(i)}$的收敛性，前者并不蕴涵后者。此外，该定理只能保证参数估计序列收敛到对数似然函数序列的稳定点，不能保证收敛到极大值点。所以在应用中，初值的选择非常重要，常用的办法是选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，择优选取。