引例:硬币模型

单硬币抛掷

假设某个盒子里有一枚硬币,抛掷该硬币只有两种结果,即正面(HeadH)和背面(TailT)。对抛硬币进行NN 次独立重复实验,其中正面朝上的概率设为pp 。记正面朝上为1,背面朝上为0,就不难得出每一次抛硬币的结果X{0,1}X\in\{0,1\} 且服从于二项分布,即Pr{X=x}=px(1p)1x\Pr\{X=x\}=p^x(1-p)^{1-x},其中pp 是未知参数。

换句话说,我们将问题建模出了一个含参的概率模型,我们的目的就是求出这个参数。统计学知识告诉我们对参数pp 的估计p^\hat p 可以通过极大似然估计(Maximum Likelihood Estimate,MLE)求得。对上述问题来说可以通过最大化似然函数求得参数:

p^=argmaxpi=1Npxi(1p)1xi\hat p=\arg\max_p\prod_{i=1}^Np^{x_i}(1-p)^{1-x_i}

通常我们可以采用取对数的方法简化求解过程,从而可以记目标函数为l(p)\mathscr l(p)

l(p)=logpi=1Nxi+log(1p)i=1N(1xi)=Nxˉlogp+N(1xˉ)log(1p)\begin{aligned} \mathscr l(p)&=\log p\sum_{i=1}^Nx_i+\log(1-p)\sum_{i=1}^N(1-x_i)\\ &=N\bar x\log p+N(1-\bar x)\log(1-p) \end{aligned}

dl(p)/dp=0\mathrm d\mathscr l(p)/\mathrm d p=0 解得p^=xˉ\hat p=\bar x ,是正面朝上占总试验次数的比例。

双硬币抛掷

现假设盒子中有A,BA,B 两枚硬币,它们正面朝上的概率分别记为θA,θB\theta_A,\theta_B。以相同的概率随机选择一个硬币,进行抛硬币实验,共做N=5N=5 次,每次实验独立的抛十次,结果如下图所示。

2coinTest

显然,在明确知道每一次实验抛的是AA 还是BB 的情况下,我们可以分别对A,BA,B 利用单硬币抛掷的方法求得参数(图中a所示)。但是,如果每次实验时并不知道抛的是哪一个硬币,那么传统的 MLE 方法就不再适用了。解决这种问题我们可以使用 EM 算法(图中b所示),这就是本文的重点内容。

EM算法及其原理

EM算法,全称为 Expectation Maximization Algorithm,译作最大期望化算法或期望最大算法。它是一种从不完全数据或有数据丢失的数据集或者说含有隐变量(hidden variable)的概率参数模型中求解参数估计迭代算法。

原理推导

对于观测数据Y\mathbf Y 和概率模型的参数Θ\mathbf\Theta ,为求参数原本可通过最大化似然函数求得,或者是最大化通常所使用的对数似然函数(log-likelihood function):

L(Θ;Y)=logP(Y;Θ)\mathcal L(\mathbf\Theta;\mathbf Y)=\log P(\mathbf Y;\mathbf\Theta)

其中分号 (;) 分隔自变量和参数,P(Y;Θ)P(\mathbf Y;\mathbf\Theta) 是全体观测样本yYy\in\mathbf Y 的联合分布,也是全体样本的似然函数。当每个样本都独立时,P(Y;Θ)=yYP(y;Θ)P(\mathbf Y;\mathbf\Theta)=\prod_{y\in\mathbf Y} P(y;\mathbf\Theta).

而引入隐变量Z\mathbf Z 之后,根据概率论的知识我们有:

P(Y;Θ)=ZP(Y,Z;Θ)=Z[P(YZ;Θ)P(Z;Θ)]P(\mathbf Y;\mathbf\Theta)=\sum_{\mathbf Z}P(\mathbf {Y,Z};\mathbf\Theta)=\sum_{\mathbf Z}\left[P(\mathbf Y|\mathbf Z;\mathbf\Theta)\cdot P(\mathbf Z;\mathbf\Theta)\right]

系联合分布P(Y,Z;Θ)P(\mathbf {Y,Z};\mathbf\Theta)Y\mathbf Y 的边缘分布(离散版本),也是边缘似然(marginal likelihood)。

接下来,我们利用 Jensen不等式,引入一个函数Q(Z)Q(\mathbf Z) 来重审目标优化函数:

L=logZP(Y,Z;Θ)=logZ(Q(Z)P(Y,Z;Θ)Q(Z))=logEQ[P(Y,Z;Θ)Q(Z)]Z(Q(Z)logP(Y,Z;Θ)Q(Z))=EQ[logP(Y,Z;Θ)Q(Z)]\begin{aligned} \mathcal L&=\log\sum_{\mathbf Z}P(\mathbf {Y,Z};\mathbf\Theta)\\ &=\log\sum_{\mathbf Z}\left(Q(\mathbf Z)\cdot\frac{P(\mathbf {Y,Z};\mathbf\Theta)}{Q(\mathbf Z)}\right) \color{grey}{=\log\mathbb E_Q\left[\frac{P(\mathbf {Y,Z};\mathbf\Theta)}{Q(\mathbf Z)}\right]}\\ &\geq \sum_{\mathbf Z}\left(Q(\mathbf Z)\cdot\log\frac{P(\mathbf {Y,Z};\mathbf\Theta)}{Q(\mathbf Z)}\right) \color{grey}{=\mathbb E_Q\left[\log\frac{P(\mathbf {Y,Z};\mathbf\Theta)}{Q(\mathbf Z)}\right]} \end{aligned}

其中保证ZQ(Z)=1\sum_{\mathbf Z}Q(\mathbf Z)=1 . 事实上,我们可以把Q(Z)Q(\mathbf Z) 视为隐变量Z\mathbf Z 的概率分布,也就是 Jensen 不等式的概率版本(上式灰色部分,可适用于连续型分布的情况,用积分表示)。

显然不等式右边是目标函数的下界,在取等条件下,极大化这个下界就等同于极大化目标函数。并且还可以进一步将与Θ\mathbf\Theta 无关的项摘掉,这并不影响极大化操作:

Θ=argmaxΘL=argmaxΘEQ[logP(Y,Z;Θ)Q(Z)]=argmaxΘEQ[logP(Y,Z;Θ)]=argmaxΘZQ(Z)logP(Y,Z;Θ)\begin{aligned} \mathbf\Theta^*=\arg\max_{\mathbf\Theta}\mathcal L &=\arg\max_{\mathbf\Theta}\mathbb E_Q\left[\log\frac{P(\mathbf {Y,Z};\mathbf\Theta)}{Q(\mathbf Z)}\right]\\ &=\arg\max_{\mathbf\Theta}\mathbb E_Q\left[\log P(\mathbf {Y,Z};\mathbf\Theta)\right]\\ &=\arg\max_{\mathbf\Theta}\sum_{\mathbf Z}Q(\mathbf Z)\log P(\mathbf {Y,Z};\mathbf\Theta) \end{aligned}

值得注意的是,此处的QQ 函数是我们人工引入的Z\mathbf Z 有关的函数,唯有这样上述的推导才是正确的!

而根据 Jensen 不等式,当且仅当P(Y,Z;Θ)Q(Z)=c\frac{P(\mathbf {Y,Z};\mathbf\Theta)}{Q(\mathbf Z)}=c 为定值时,可以取等。从而我们可以推导确定QQ 函数的表达式:

Q(Z)=P(Y,Z;Θ)ZP(Y,Z;Θ)=P(Y,Z;Θ)P(Y;Θ)=P(ZY;Θ)\begin{aligned} Q(\mathbf Z) = \frac{P(\mathbf {Y,Z};\mathbf\Theta)}{\sum_{\mathbf Z}P(\mathbf {Y,Z};\mathbf\Theta)} = \frac{P(\mathbf {Y,Z};\mathbf\Theta)}{P(\mathbf Y;\mathbf\Theta)} = P(\mathbf Z|\mathbf Y;\mathbf\Theta) \end{aligned}

这个表达式与Θ\mathbf\Theta 有关,所以 EM 算法作为一个迭代算法,提出了一个 2-step 操作。在每次迭代时分别对Z\mathbf Z 的分布和Θ\mathbf\Theta 的选取进行计算:

  • E step (Expectation):固定当前迭代次数tt 已得到的参数Θ(t)\mathbf\Theta^{(t)} 计算P(ZY;Θ(t))P(\mathbf Z|\mathbf Y;\mathbf\Theta^{(t)})
  • M step (Maximization):固定Z\mathbf Z 的分布,最优化目标函数求得当前最优的参数Θ(t+1)\mathbf\Theta^{(t+1)} ,即Θ(t+1)=argmaxΘL(Θ,Θ(t))\mathbf\Theta^{(t+1)}=\arg\max_{\mathbf\Theta}\mathcal L(\mathbf\Theta,\mathbf\Theta^{(t)})

此外算法的收敛性同样需要我们思考:

  1. EM算法能保证收敛吗?
  2. EM算法如果收敛,那么能否保证收敛到全局最大值?

详见:机器学习-白板推导系列(十)-EM算法

算法流程

虽然推导看似很复杂,但是我们可以总结出这样一个流程:

  1. 输入观测样本数据Y\mathbf Y ,隐变量数据Z\mathbf Z,联合分布(或联合密度函数)P(Y,Z;Θ)P(\mathbf {Y,Z};\mathbf\Theta) 和条件分布(或条件密度函数)P(ZY;Θ)P(\mathbf Z|\mathbf Y;\mathbf\Theta) 以及迭代次数TT;
  2. 随机初始化参数的初值Θ(0)\mathbf\Theta^{(0)}
  3. E步→M步;
  4. 判断参数是否收敛,否则继续上一步操作。

双硬币抛掷问题

2coinTest

回看双硬币抛掷问题,我们使用 EM 算法来进行求解(图b所示)。

首先是初始化欲估参数θA(0)=0.6,θB(0)=0.5\theta_A^{(0)}=0.6,\theta_B^{(0)}=0.5。规定隐变量zz 表示当前抛掷的硬币是AA 还是BB。接下来计算联合分布:

  • 第一行观测数据(5个正面5个背面)是AA 的概率P(Y1z=A)=P(Y_1|z=A)=

非标准多项式分布

ECM 算法

参考

  1. EM算法 - LeeLIn。 - 博客园
  2. EM算法原理总结 - 刘建平Pinard - 博客园
  3. EM算法及其推广-码农场
  4. 机器学习-白板推导系列(十)-EM算法(Expectation Maximization)学习笔记 - 知乎