EM算法

EM算法引入

EM算法：含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计法。

假设有3枚硬币A，B和C，正面出现的概率分别为$\pi$，$p$和$q$。进行如下的掷币实验：先掷硬币A，根据结果选出硬币B或C，正面选B，反面选C；然后掷选出的硬币，掷硬币的结果，出现正面记作1，出现反面记作0；独立地重复$n$次实验($n=10$)，观测结果如下：

$$
1,1,0,1,0,0,1,0,1,1
$$

设随机变量$y$为观测变量，表示一次实验的观测结果是0或1,；随机变量$z$是隐变量，表示未观测到的掷硬币A的结果；$\theta=(\pi, p, q)$是模型参数。

三硬币模型可以写作：

$$
P(y|\theta) = \sum\limits_z P(y,z|\theta)
= \sum\limits_z P(z|\theta) P(y|\theta, z)
=\pi p^y (1-p)^{1-y} + (1-\pi) q^y (1-q)^{1-y}
$$

类似的，如果将观测数据写作$Y=(Y_1, \dots, Y_n)^T$，未观测数据写作$Z=(Z_1, \dots, Z_n)^T$，则观测数据的似然函数是：

$$
P(Y|\theta)=\sum\limits_Z P(Z|\theta) P(Y|Z, \theta)
$$

即：

$$
P(Y|\theta) = \prod \limits_{j=1}^n \left[ \pi p^{y_j} (1-p)^{1-y_j} + (1-\pi) q ^{y_j} (1-q) ^{1-y_j} \right]
$$

求$\theta$的极大似然估计，即

$$
\hat\theta = \arg \max\limits_\theta \log P(Y|\theta)
$$

E步：计算模型参数$\pi^{(i)},p^{(i)},\epsilon^{(i)}$下观测数据来自掷硬币B的概率：

$$
\mu^{(i+1)}_j = \frac{\pi^{(i)} (p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}}{\pi^{(i)} (p^{(i)})^{1-y_j} + (1-\pi^{(i)}) (q^{(i)})^{y_j} (1-q^{(i)}) ^ {1-y_j}}
$$

M步：更新模型参数的估计值

$$
\pi^{(i+1)} = \frac{1}{n} \sum\limits_{j=1}^n \mu_j^{i+1}
$$

$$
p^{(i+1)} = \frac{ \sum\limits_{j=1}^n \mu_j^{(i+1)} y_j }{ \sum\limits_{j=1}^n \mu_j^{(i+1)} }
$$

$$
q^{(i+1)} = \frac {\sum\limits_{j=1}^n \left( 1-\mu_j^{(i+1)} \right) y_j} {\sum\limits_{j=1}^n \left( 1-\mu_j^{(i+1)} \right)}
$$

EM算法

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

输出：模型参数$\theta$

（1）选择参数的初值$\theta^{(0)}$，开始迭代。

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

$$
Q(\theta, \theta^{(i+1)}) = E_Z \left[ \log P(Y,Z|\theta) | Y, \theta^{(i)} \right]
= \sum\limits_Z \log P(Y,Z|\theta) P(Z|Y, \theta^{(i)})
$$

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

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

$$
\theta^{(i+1)} = \arg \max \limits_\theta Q(\theta, \theta^{(i)})
$$

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

上述算法中的函数$Q(\theta, \theta^{(i)})$是EM算法的核心，称为Q函数。

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

$$
Q(\theta, \theta^{(i+1)}) = E_Z \left[ \log P(Y,Z|\theta) | Y, \theta^{(i)} \right]
$$

说明：

1. 参数$\theta$的初值$\theta^{(0)}$可随意选取，但EM算法对初值敏感。
2. E步求$Q(\theta, \theta^{(i)})$，Q函数式中$Z$是未观测数据，$Y$是观测数据，$\theta$是要极大化的参数，$\theta^{(i)}$是参数的当前估计值。每次迭代求Q及其极大。
3. M步求Q的极大化得到$\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+1)}) - Q(\theta^{(i)}, \theta^{(i+1)}) || \lt \epsilon_2
$$

EM算法的导出

目标是极大化观测数据$Y$关于参数$\theta$的对数似然函数，即极大化

$$
\begin {split}
L(\theta) &= \log P(Y|\theta) \\
&= \log \sum\limits_Z P(Y,Z|\theta) \\
&= \log \left( \sum\limits_Z P(Y|Z, \theta) P(Z|\theta) \right)
\end {split}
$$

极大化上式的困难在于含有未观测数据$Z$及和的对数。

我们希望通过迭代的方式逐步近似地去极大化$L(\theta)$，即经过一次迭代$\theta^{(i)} \rightarrow \theta^{(i+1)}$后，$L(\theta)$能够增加，即$L(\theta) \gt L(\theta^{(i)})$。

考虑两者的差，

$$
L(\theta) - L(\theta^{(i)}) = \log \left( \sum\limits_Z P(Y|Z, \theta) P(Z|\theta) \right) - \log P(Y|\theta^{(i)})
$$

利用Jensen不等式：

$$
\log \sum\limits_j \lambda_j y_j \ge \sum\limits_j \lambda_j \log y_j, \lambda_j \ge 0, \sum\limits_j \lambda_j = 1
$$

得到上式的下届：

$$
\begin {split}
L(\theta) - L(\theta^{(i)}) &= \log( \sum\limits_Z P(Y|Z, \theta^{(i)}) \frac{P(Y|Z, \theta)P(Z|\theta)}{ P(Y|Z, \theta^{(i)}) } ) - \log P(Y|\theta^{(i)}) \\
& \ge \sum\limits_Z P(Z|Y, \theta^{(i)}) \log \frac{ P(Y|Z, \theta)P(Z|\theta) }{ P(Y|Z, \theta^{(i)}) } - \log P(Y|\theta^{(i)}) \\
& = \sum\limits_Z P(Z|Y, \theta^{(i)}) \log \frac {P(Y|Z, \theta) P(Z, \theta)} { P(Y|Z, \theta^{(i)}) P(Y|\theta^{(i)}) }
\end{split}
$$

令

$$
B(\theta, \theta^{(i)}) = L(\theta^{(i)}) + \sum\limits_Z P(Z|Y, \theta^{(i)}) \log \frac {P(Y|Z, \theta) P(Z, \theta)} {P(Y|Z, \theta^{(i)}) P(Y|\theta^{(i)})}
$$

则有$L(\theta) \ge B(\theta, \theta^{(i)})$， 且有$L(\theta^{(i)}) = B(\theta^{(i)}, \theta^{(i)})$。

因此，任何可以使$B(\theta, \theta^{(i)})$增大的$\theta$，也可以使$L(\theta)$增大。为了使$L(\theta)$尽可能增大，选择$\theta^{(i+1)}$使$B(\theta, \theta^{(i)})$达到极大，即

$$
\begin{split}
\theta^{(i+1)} &= \arg \max \limits_\theta B(\theta, \theta^{(i)}) \\
&= \arg \max \limits_\theta \left( L(\theta^{(i)}) + \sum\limits_Z P(Z|Y, \theta^{(i)}) \log \frac {P(Y|Z, \theta) P(Z, \theta)} {P(Y|Z, \theta^{(i)}) P(Y|\theta^{(i)})} \right) \\
&= \arg \max \limits_\theta \left( \sum\limits_Z P(Z|Y,\theta^{(i)}) \log (P(Y|Z,\theta)P(Z|\theta)) \right) \\
&= \arg \max \limits_\theta \left( \sum\limits_Z P(Z|Y, \theta^{(i)}) \log P(Y, Z|\theta) \right) \\
&= \arg \max \limits_\theta Q(\theta, \theta^{(i)})
\end{split}
$$

EM算法的收敛性

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

$$
P(Y|\theta^{(i)}) \ge P(Y|\theta^{(i)})
$$

证明：

由于$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)
$$

由于：

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

令：

$$
H(\theta, \theta^{(i)}) = \sum\limits_Z \log P(Y, Z|\theta) P(Z|Y, \theta^{(i)})
$$

于是对数似然函数可以写成：

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

上式的$\theta$分别取$\theta^{(i+1)}$和$\theta^{(i)}$并相减，有

$$
\begin{split}
& \log P(Y|\theta^{(i+1)}) - \log P(Y|\theta^{(i)}) \\
=& \left[Q(\theta^{(i+1)}, \theta^{(i)}) - Q(\theta^{(i)}, \theta^{(i)}) \right] - \left[ H(\theta^{(i+1)}, \theta^{(i)}) - H(\theta^{(i)}, \theta^{(i)}) \right]
\end{split}
$$

由于$\theta^{(i+1)}$使$Q(\theta, \theta^{(i)})$达到最大，于是有

$$
Q(\theta^{(i+1)}, \theta^{(i)}) - Q(\theta^{(i)}, \theta^{(i)}) \ge 0
$$

而

$$
\begin{split}
& H(\theta^{(i+1)}, \theta^{(i)}) - H(\theta^{(i)}, \theta^{(i)}) \\
=& \sum\limits_Z \left( \log \frac{P(Z|Y,\theta^{(i+1)})} {P(Z|Y,\theta^{(i)})} \right)P(Z|Y,\theta^{(i)}) \\
\le & \log \left( \sum\limits_Z \frac{P(Z|Y, \theta^{(i+1)})} {P(Z|Y, \theta^{(i)})} P(Z|Y, \theta^{(i)}) \right) （Jensen不等式） \\
=& \log \left( \sum\limits_Z P(Z|Y, \theta^{(i+1)}) \right) \\
=& 0
\end{split}
$$

于是有$\log P(Y|\theta^{(i+1)}) - \log P(Y|\theta^{(i)}) \le 0$，证毕。

EM算法在高斯混合模型中的应用

高斯混合模型的定义

具有如下形式的概率分布模型：

$$
P(y|\theta)=\sum\limits_{k=1}^K \alpha_k \phi(y|\theta_k)
$$

其中，$\alpha_k$是系数，$\alpha_k \gt 0$，$\sum\limits_{k=1}^K \alpha_k = 1$。

$\phi(u|\theta_k)$是高斯分布密度，$\theta_k=(\mu_k, \sigma_k^2)$，

$$
\phi(y|\theta_k)=\frac{1}{\sqrt{2\pi}\sigma_k} \exp \left(- \frac{(y-\mu_k)^2}{2\sigma_k^2}\right)
$$

称为第$k$个分模型。

一般混合模型可以由任意概率分布密度代替上式中的高斯分布密度。

高斯混合模型参数估计中的EM算法

假设观测数据$y_1, y_2, \dots, y_N$由高斯混合模型生成：

$$
P(y|\theta)=\sum\limits_{k=1}^K \alpha_k \phi(y|\theta_k)
$$

其中，$\theta=(\alpha_1, \alpha_2, \dots, \alpha_K; \theta_1, \theta_2, \dots, \theta_K)$。

下面用EM算法估计高斯混合模型的参数$\theta$。

• 1 明确隐变量，写出完全数据的对数似然函数。

设隐变量$\gamma_{jk}$反映观测数据$y_j$来自第$k$个分模型。

$$
\gamma_{jk} =
\begin{cases}
1, 第j个观测数据来自第k个分模型 \\
0, 否则
\end{cases}
$$

$\gamma$是$0-1$随机变量。

有了观测数据$y_j$以及未观测数据$\gamma_{jk}$，那么完全数据是

$$
(y_j, \gamma_{j1}, \gamma_{j2}, \dots, \gamma_{jk}), j=1,2,\dots,N
$$

于是可以写出完全数据的似然函数：

$$
\begin{split}
P(y, \gamma | \theta) =& \prod \limits_{j=1}^N P(y_j, \gamma_{j1}, \gamma_{j2}, \dots, \gamma_{jk} | \theta) \\
=& \prod \limits_{k=1}^K \prod\limits_{j=1}^N \left[\alpha_k \phi (j_j|\theta_k) \right] ^ { \gamma _ {jk} } \\
=& \prod \limits_{k=1}^{K} \alpha_{k}^{n_k} \prod \limits_{j=1}^N \left[ \phi(y_j|\theta_k) \right]^{\gamma_jk} \\
=& \prod \limits_{k=1}^{K} \alpha_{k}^{n_k} \prod \limits_{j=1}^N \left[ \frac{1}{\sqrt{2\pi}\sigma_k} \exp \left(- \frac{(y_j-\mu_k)^2}{2\sigma_k^2}\right) \right]^{\gamma_{jk}}
\end{split}
$$

其中$n_k = \sum\limits_{j=1}^N \gamma_{jk}$， $\sum\limits_{k=1}^K n_k = N$。

于是，完全数据的对数似然函数为：

$$
\log P(y, \gamma |\theta) = \sum\limits_{k=1}^K \left( n_k \log \alpha_k + \sum\limits_{j=1}^N \gamma_{jk} \left[\log (\frac{1}{\sqrt{2\pi}}) - \log \sigma_k - \frac{1}{2\sigma_k^2} (y_j - \mu_k)^2\right] \right)
$$

• 2 E步：确定Q函数。

$$
\begin{split}
Q(\theta, \theta^{(i)}) =& E[\log P(y, \gamma | \theta) | y, \theta^{(i)}] \\
=& \sum\limits_{k=1}^K \left( n_k \log \alpha_k + \sum\limits_{j=1}^N \gamma_{jk} \left[\log (\frac{1}{\sqrt{2\pi}}) - \log \sigma_k - \frac{1}{2\sigma_k^2} (y_j - \mu_k)^2\right] \right) \\
=& \sum\limits_{k=1}^K \left( \sum\limits_{j=1}^N (E\gamma_{jk})\log\alpha_k + \sum\limits_{j=1}^N(E\gamma_{jk})\left( \log(\frac{1}{\sqrt{2\pi}}) - \log\sigma_k - \frac{(y_j-\mu_k)^2}{2\sigma^2} \right) \right)
\end {split}
$$

这里需要计算$E(\gamma_{jk} | y,\theta)$，记为$\hat\gamma_{jk}$。

$$
\begin{split}
\hat\gamma_{jk} =& E(\gamma_{jk} | y,\theta) \\
=& P(\gamma_{jk}=1|y,\theta) \\
=& \frac {P(\gamma_{jk}=1,y_j|\theta)} {\sum\limits_{k=1}^K P(\gamma_{jk}=1,y_j | \theta)} \\
=& \frac {P(y_j|\gamma_{jk}=1,\theta) P(\gamma_{jk}=1|\theta)} {\sum\limits_{k=1}^K P(y_j | \gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)} \\
=& \frac{\alpha_k \phi(y_j | \theta_k)} {\sum\limits_{k=1}^K \alpha_k \phi(y_j | \theta_k)} , j=1,2,\dots,N;k=1,2,\dots,K
\end {split}
$$

$\hat\gamma_{jk}$是当前模型参数下第$j$个观测数据来自第$k$个分模型的概率，称为分模型$k$对观测数据$y_j$的响应度。

由$\hat\gamma_{jk}=E(\gamma_{jk})$及$n_k = \sum\limits_{j=1}^N \gamma_{jk}$，有

$$
Q(\theta, \theta^{(i)}) = \sum\limits_{k=1}^K \left( n_k \log \alpha_k + \sum\limits_{j=1}^N \left[\log (\frac{1}{\sqrt{2\pi}}) - \log \sigma_k - \frac{1}{2\sigma_k^2} (y_j - \mu_k)^2\right] \right)
$$

• 3 M步，求Q函数对$\theta$的极大值。

用$\hat\mu_k$，$\hat\sigma_k^2$及$\hat\alpha_k,k=1,2,\dots,K$表示$\theta^{(i+1)}$的各参数。求$\hat\mu_k,\hat\sigma_k^2$只需对上式求偏导令为0；求$\hat\alpha_k$是在$\sum\limits_{k=1}^K \alpha_k=1$的条件下求偏导令等于0.

于是：

$$
\begin{split}
\hat\mu_k &= \frac{\sum\limits_{j=1}^N \hat\gamma_{jk} y_j} {\sum\limits_{j=1}^N \hat\gamma_{jk}} \\
\sigma_k^2 &= \frac{\sum\limits_{j=1}^N \hat\gamma_{jk} (y_j - \mu_k)^2} { \sum\limits_{j=1}^N \hat\gamma_{jk}} \\
\hat\alpha_k &= \frac{n_k}{N} = \frac {\sum\limits_{j=1}^N \gamma_{jk}} {N} \\
k &= 1,2,\dots,K
\end{split}
$$

重复以上计算，直到对数似然函数不再有明显变化为止。