EM算法推导与三硬币模型

EM算法概览

EM算法是一种极大似然法，专门处理含有隐变量的概率模型。隐变量不能直接观测得到，为了得到对隐变量分布的一个近似估计，EM算法分两步进行迭代求解：

E步：使用当前模型参数求出一个期望
M步：对期望求极大
一个使用EM算法求解的例子：
有三枚硬币A，B，C，先投掷A，如果是正面就投掷B，如果是反面就投掷C，若我们只能观测到最后的投掷结果（B或者C的结果），如何估算三颗硬币的正面率？
在这里面，每次A的结果不能从观测变量中直接得到，因此A是一个隐变量。以下部分会先推导EM算法，然后将其应用在三硬币问题上。

EM算法推导

首先考虑最一般化的极大似然，我们在获得观测变量$Y$时，希望最大化$Y$出现的概率，即极大化目标函数

$L(\theta)=\log P(Y | \theta)=\log \sum_{Z} P(Y, Z | \theta)=\log \left(\sum_{Z} P(Y | Z, \theta) P(Z | \theta)\right)$

其中，$\theta$是模型的参数，$Z$是隐变量。上面使用到了全概率公式和条件概率公式。我们使用迭代的方式逐步更新$\theta$的值，试图在迭代的过程中增大$L(\theta)$的值。那么如何更新$\theta$呢？假设我们现在已经进行了i轮的迭代，当前的模型参数为$\theta^{(i)}$,我们来考虑一下迭代后的目标函数$L(\theta)$和迭代前的目标函数$L(\theta^{(i)})$的差值（注:$L(\theta^{(i)})$是指使用当前模型的参数计算出来的目标函数值)

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

我们可以利用Jensen不等式来进行一个不等式放缩,对于一个凸函数， Jensen不等式为：

$f(E[x]) \geq E[f(x)]$

或者

$f\left(\sum_{i=1}^{M} \lambda_{i} x_{i}\right) \geq \sum_{i=1}^{M} \lambda_{i} f\left(x_{i}\right),其中\sum_{i=1}^{M} \lambda_{i}=1$

那么原式

$L(\theta)-L\left(\theta^{(i)}\right)=\log \left(\sum_{Z} P(Y | Z, \theta) P(Z | \theta)\right)-\log P\left(Y | \theta^{(i)}\right)$ $=\log \left(\sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(i)}\right)}\right)-\log P\left(Y | \theta^{(i)}\right)$ $\geqslant \sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \log \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(i)}\right)}-\log P\left(Y | \theta^{(i)}\right)$ $= \sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \log \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(i)}\right)}-\sum_{Z} P\left(Z | Y, \theta^{(i)}\right)\log P\left(Y | \theta^{(i)}\right)$ $= \sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \left[\log \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(i)}\right)}-P\left(Y | \theta^{(i)}\right)\right]$ $=\sum_{z} P\left(Z | Y, \theta^{(i)}\right) \log \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(t)}\right) P\left(Y | \theta^{(i)}\right)}$

其中，第一个等号在原式子额外添加了一个分子和分母，分子部分作为概率分布，分母和原有部分组成一个值，那么整个 $\sum_{Z} P\left(Y | Z, \theta^{(i)}\right) \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Y | Z, \theta^{(i)}\right)}$ 部分变成了一个求解期望的式子,利用Jensen不等式可以得到下一行的前半部分。第三行后半部分可以添加一个乘积的原因是， $\sum_{Z} P\left(Z | Y, \theta^{(i)}\right)$ 是一个概率分布，加起来为1，而$\log P\left(Y | \theta^{(i)}\right)$中不含有$Z$，所以可以和前面的求和乘起来（相当于一个常数）。整理一下上面的结论，我们有

$L(\theta)-L\left(\theta^{(i)}\right)\geqslant\sum_{z} P\left(Z | Y, \theta^{(i)}\right) \log \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(t)}\right) P\left(Y | \theta^{(i)}\right)}$

引入一个辅助函数 $B(\theta,\theta^{(i)})=L\left(\theta^{(i)}\right)+\sum_{z} P\left(Z | Y, \theta^{(i)}\right) \log \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(t)}\right) P\left(Y | \theta^{(i)}\right)}$ ,
那么有

$L(\theta) \geqslant B\left(\theta, \theta^{(i)}\right)$

说明$B(\theta,\theta^{(i)})$是$L(\theta)$的一个下界，并且有

$L\left(\theta^{(i)}\right)=B\left(\theta^{(i)}, \theta^{(i)}\right)$

因此我们的问题变成了，如何最大化$B(\theta,\theta^{(i)})$，只要可以提高$B(\theta,\theta^{(i)})$的值，那么就能提高$L(\theta)$的值。（注：$\theta^{i}$是已知参数或者说常数，代表当前的变量的值，$\theta$是变量，我们需要求解$\theta$的极值点）。记迭代后的值为$\theta^{(i+1)}$,

$\theta^{(i+1)}=\arg \max _{\theta}\left(L\left(\theta^{(i)}\right)+\sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \log \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(i)}\right) P\left(Y | \theta^{(i)}\right)}\right)$

省去和$\theta$无关的常数，变形如下

$原式=\arg \max _{\theta}\left(\sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \log (P(Y | Z, \theta) P(Z | \theta))\right)$ $=\arg \max _{\theta}\left(\sum_{z} P\left(Z | Y, \theta^{(i)}\right) \log P(Y, Z | \theta)\right)$

我们记 $Q\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{(i)}\right)=\sum_{z} P\left(Z | Y, \theta^{(i)}\right) \log P(Y, Z | \theta)$ ，称为Q函数，至此我们就可以根据不同的任务设计好Q函数，然后求导使得Q函数极大化，增加Q函数的值，从而增加极大似然函数的值。

另一方面，Q函数也可以看做是直接对$Y,Z$的极大似然函数求期望得到，即

$Q\left(\theta, \theta^{(i)}\right)=E_{Z}\left[\log P(Y, Z | \theta) | Y, \theta^{(i)}\right]$ $=\sum_{Z} P(Z | Y, \theta^{(j)}) \log P(Y, Z | Y, \theta)$

给定$Y$,$Z$是随机变量，我们求目标极大似然函数的期望，并希望得到的极大似然函数最大，实际上问题也变成了求Q函数的极大。

整个EM算法的框架如下：

输入：观察变量Y，隐变量数据Z，联合分布$P(Y,Z|\theta)$,条件分布$P(Z|Y,\theta)$
模型参数:$\theta$
(1) 初始化参数值$\theta^{(0)}$，开始迭代
(2) E步：记$\theta^{(i)}$为当前模型的参数估计值，在第i+1次迭代的E步，我们计算Q函数
$Q\left(\theta, \theta^{(i)}\right)=E_{Z}\left[\log P(Y, Z | \theta) | Y, \theta^{(i)}\right]$ $=\sum_{Z} P(Z | Y, \theta^{(j)}) \log P(Y, Z | \theta)$
E步的输出是当前的Q函数，其中的$\theta$是待优化的参数
（3）M步：求Q函数的极值点使其极大化，得到新的参数
$\theta^{(i+1)}=\arg \max _{\theta} Q\left(\theta, \theta^{(i)}\right)$
(4) 重复(2)(3)直到收敛
注：EM算法对初始值敏感，不同的初始值可能导致不同的结果，EM算法也无法保证找到全局最优点

下面针对三硬币模型进行一个EM算法的应用讲解

三硬币模型

问题描述

有三枚硬币A，B，C，先投掷A，如果是正面就投掷B，如果是反面就投掷C，若我们只能观测到最后的投掷结果（B或者C的结果）而不能直到投掷的过程，如何估算三颗硬币的正面率？

形式化

我们将A，B，C正面朝上的概率分别设为$\pi,p,q$,最后的观察结果随机变量记为$Y(Y=0,1)$,变量$Z$表示硬币A的结果，$Z=0/1$。

建立概率模型

目标函数

首先我们的极大似然函数为 $L(\theta)=\prod_{j=1}^{n}P(Y|\theta)$
其中$n$为样本数量，$j$为样本编号。取了log之后为

$L(\theta)=\sum^{n}_{j=1}\log P(Y | \theta)=\sum^{n}_{j=1}\log \sum_{Z} P(Y, Z | \theta)=\sum^{n}_{j=1}\log \left(\sum_{Z} P(Y | Z, \theta) P(Z | \theta)\right)$

由上面的推导可以知道，极大化$L(\theta)$可以通过极大化Q函数得到 $Q(\theta,\theta^{(i)})=\sum^{n}_{j=1}\sum_{Z} \log P(Y, Z | \theta) P\left(Z | Y, \theta^{(i)}\right)$

条件概率

我们先来计算$P\left(Z | Y, \theta^{(i)}\right)$,

$P\left(Z | Y, \theta^{(i)}\right)=\frac{P(Y,Z|\theta^{(i)})}{P(Y|\theta^{(i)})}$ $P(Y=y_j | \theta^{(i)})=\pi^{(i)}\left(p^{(i)}\right)^{y_{j}}\left(1-p^{(i)}\right)^{1-y_{j}}+\left(1-\pi^{(i)}\right)\left(q^{(i)}\right)^{y_{j}}\left(1-q^{(i)}\right)^{1-y_{j}}$

联合分布

关于$P(Y, Z | \theta)$,我们有

$P(Y, Z | \theta)=Z\pi p^Y(1-p)^{1-Y}+(1-Z)\pi q^Y(1-q)^{1-Y}$

特别的，我们有

$P(Y, Z | \theta^{(i)})=P(Y, Z=0 | \theta^{(i)})+P(Y, Z=1 | \theta^{(i)})$ $P(Y=y_j, Z = 1 | \theta^{(i)}) =\pi^{(i)}\left(p^{(i)}\right)^{y_{j}}\left(1-p^{(i)}\right)^{1-y_{j}}$

求解Q函数（E步）

由上面的结果得

$P\left(Z=1 | Y=y_j, \theta^{(i)}\right)=\frac{P(Y,Z=1|\theta^{(i)})}{P(Y=y_j|\theta^{(i)})}\\ =\frac{\pi^{(i)}\left(p^{(i)}\right)^{y_{j}}\left(1-p^{(i)}\right)^{1-y_{j}}}{\pi^{(i)}\left(p^{(i)}\right)^{y_{j}}\left(1-p^{(i)}\right)^{1-y_{j}}+\left(1-\pi^{(i)}\right)\left(q^{(i)}\right)^{y_{j}}\left(1-q^{(i)}\right)^{1-y_{j}}}$

记$u_j^{(i)}=P\left(Z=1 | Y=y_j, \theta^{(i)}\right)$,那么
$P\left(Z=0 | Y=y_j, \theta^{(i)}\right) =1-P\left(Z=1 | Y=y_j, \theta^{(i)}\right)=1-u_j^{(i)}$

我们要优化的Q函数为

$Q(\theta,\theta^{(i)})=\sum^{n}_{j=1}\sum_{Z} \log P(Y, Z | \theta) P\left(Z | Y, \theta^{(i)}\right)\\ =\sum^{n}_{j=1}\log P(Y=y_j, Z=1 | \theta) P\left(Z=1 | Y=y_j, \theta^{(i)}\right)\\+\log P(Y=y_j, Z=0 | \theta) P\left(Z=0 | Y=y_j, \theta^{(i)}\right)\\ =\sum^{n}_{j=1}u_j^{(i)}\log \left(\pi p^{y_j}(1-p)^{1-y_j}\right) +(1-u_j^{(i)})log\left((1-\pi)q^{y_j}(1-q)^{1-y_j}\right)\\ =\sum^{n}_{j=1}u_j^{(i)}\left[log\pi+y_j\log p+(1-y_j)\log(1-p)\right]+\\ (1-u_j^{(i)})\left[log(1-\pi)+y_j\log q+(1-y_j)\log(1-q)\right]$

求Q函数极大(M步)

我们分别求参数$\pi,p,q$对Q函数的偏导

$\frac{\partial Q}{\partial \pi}=\sum^{n}_{j=1}u_j^{(i)}\frac{1}{\pi}-(1-u_j^{(i)})\frac{1}{1-\pi}= \sum^{n}_{j=1}\frac{u_j^{(i)}-\pi}{\pi(1-\pi)}$

令$\frac{\partial Q}{\partial \pi}=0$,得 $\pi^{(i+1)} = \frac{1}{n} \sum_{j=1}^{n} u_{j}^{(i)}$

继续求$ p $的偏导

$\frac{\partial Q}{\partial p}=\sum^{n}_{j=1}u_j^{(i)}(\frac{y_j}{p}-\frac{1-y_j}{1-p})=\sum^{n}_{j=1}u_j^{(i)}\frac{y_j-p}{p(1-p)}$

令$\frac{\partial Q}{\partial p}=0$,得

$p^{(i+1)} = \frac{\sum^{n}_{j=1} u_{j}^{(i)} y_{j} } {\sum^{n}_{j=1}u_{j}^{(i)}}$

同理，我们得到

$q^{(i+1)}=\frac{\sum^{n}_{j=1}\left(1-u_{j}^{(i)}\right) y_{j}}{\sum^{n}_{j=1}\left(1-u_{j}^{(i)}\right)}$

至此，我们的所有参数更新完成，进入下一轮迭代，EM算法一轮完成。

参考文献

《统计学习方法》,李航