本篇博客主要介绍李航《统计学习方法(第2版)》中讲解EM算法涉及到的三硬币模型案例，原文中该模型的推导过程被省略了。本篇博客主要是将该模型的具体推导过程。

1 三硬币模型

  假设有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$$目前只能观测到掷硬币的结果，不能观测掷硬币的过程。EM算法要解决的问题就是在缺少掷硬币的过程信息的情况下如何估计$\pi$，$p$和$q$的值。

2 推导过程

  假设观测数据记为$Y=(y_1,y_2,\dots ,y_n{)}^T$，其中$y_i$表示第$i$次试验的观测结果是1或0；隐变量数据记为$Z=(z_1,z_2,\dots ,z_n{)}^T$，其中$z_i$表示第$i$次未观测到的掷硬币A的结果。$\theta =(\pi ,p,q)$记为模型参数。对于任意一次试验中的观测变量$y_i$，其概率为$$\begin{array}{rl}P(y_i|\theta )& =\sum _{z_i}P(y_i,z_i|\theta )=\sum _{z_i}p(z_i|\theta )P(y_i|z_i,\theta )\\ & =\pi p^{y_i}(1-p{)}^{1-y_i}+(1-\pi )q^{y_i}(1-q{)}^{1-y_i}\end{array}$$所以，观测数据$Y$的似然函数为$$\begin{array}{rl}P(Y|\theta )& =\sum _{Z}P(Z|\theta )P(Y|Z,\theta )\\ & =\prod _{i=1}^n[\pi p^{y_i}(1-p_i{)}^{1-y_i}+(1-\pi )q^{y_i}(1-q{)}^{1-y_i}]\end{array}$$若使用传统的极大似然方法求解参数$\theta$，其解为$$\hat{\theta }=arg\underset{\theta }{\max }logP(Y|\theta )$$但该公式并不存在解析解，该方法行不通。这时就需要用到EM算法了。EM算法是一个不断迭代的过程，直到满足终止条件则结束。以下过程主要展示如何利用第$i$次的参数结果推导出第$i+1$次的参数。

• 首先给模型参数指定初始值。这里记初值为${\theta }^{(0)}=({\pi }^{(0)},p^{(0)},q^{(0)})$，
• 假设我们现在已经得到第$i$次迭代后的参数值，记为${\theta }^{(i)}=({\pi }^{(i)},p^{(i)},q^{(i)})$。现在我们要做的是：${\theta }^{(i)}\to {\theta }^{(i+1)}$。
• 根据${\theta }^{(i)}$估计隐藏数据的值并计算$Q$函数(E步)；
  这里要估计的是隐藏数据$Z$出现的概率，也就是每次试验中选中硬币B或硬币C的概率。在得到${\theta }^{(i)}$后，观测数据$y_j$来自硬币$B$的概率如下：$$u_j^{(i+1)}=P(z_j=B|y_j,{\theta }^{(i)})=\frac{{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}}{{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}+(1-{\pi }^{(i)})(q^{(i)}{)}^{y_j}(1-q^{(i)}{)}^{1-y_j}}$$观测数据$y_j$来自硬币C的概率为:$1-u_j^{(i+1)}$. 计算$Q$函数(这里不讲$Q$函数的推导，书上讲的很详细)，如下：$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]\\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )\end{array}$$这里要介绍一下$Q$函数代表的意义：完全数据的对数似然函数$logP(Y,Z|\theta )$在给定观测数据$Y$和当前参数${\theta }^{(i)}$下对未观测数据$Z$的条件概率分布$P(Z|Y,{\theta }^{(i)})$的期望。而所求的${\theta }^{(i+1)}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$。
  这里补充一点关于条件概率期望的知识。在$X=x$下随机变量$Y$的期望的计算公式如下：$$E(Y|X=x)=\{\begin{array}{ll}\sum y_{j}P(Y=y_j|X=x)& \text{(X,Y)为二维离散随机变量}\\ {\int }_{-\infty }^{\infty }yp(y|x)dy& \text{(X,Y)为二维连续随机变量}\end{array}$$接下来，为了能够完整地展现$Q(\theta ,{\theta }^{(i)})$的计算过程，我们先将这个三硬币模型的试验次数调整为2。
  假设现在在2次随机试验过程下得到的观测数据为$Y=(y_1,y_2)$。那么$Q(\theta ,{\theta }^{(i)})$所代表的就是每一种可能的观测数据集合$Z$出现的条件概率与完全数据的似然函数$logP(Y,Z|\theta )$的乘积之和。
  现在取$Z=(1,0)$，那么此时$P(Z|Y,{\theta }^{(i)})=u_1^{(i+1)}(1-u_2^{(i+1)})$，而$$\begin{array}{rl}logP(Y,Z|\theta )& =log[\pi p^{y_1}(1-p{)}^{1-y_1}(1-\pi )q^{y_2}(1-q{)}^{1-y_2}]\\ & =log\pi +y_{1}logp+(1-y_1)log(1-p)+log(1-\pi )+y_{2}logq+(1-y_2)log(1-q)\end{array}$$补充上$Z=(1,1)$、$Z=(0,1)$、$Z=(0,0)$的情况下的$P(Z|Y,{\theta }^{(i)})$和$logP(Y,Z|\theta )$后即可得到对应的$Q(\theta ,{\theta }^{(i)})$。这里记为$Q_2(\theta ,{\theta }^{(i)})$，那么$$\begin{array}{rl}{Q}_2(\theta ,{\theta }^{(i)})& =u_1^{i+1}(1-u_2^{(i+1)})(log\pi +log(1-\pi )+y_{1}logp+(1-y_1)log(1-p)+y_{2}logq+(1-y_2)log(1-q))\\ & +u_1^{(i+1)}{u}_2^{(i+1)}(log\pi +log\pi +y_{1}logp+(1-y_1)log(1-p)+y_{2}logp+(1-y_2)log(1-p))\\ & +(1-u_1^{(i+1)})u_2^{(i+1)}(log(1-\pi )+log\pi +y_{1}logq+(1-y_1)log(1-q)+y_{2}logp+(1-y_2)log(1-p))\\ & +(1-u_1^{(i+1)})(1-u_2^{(i+1)})(log(1-\pi )+log(1-\pi )+y_{1}logq+(1-y_1)log(1-q)+y_{2}logq+(1-y_2)log(1-q))\\ & =(u_1^{(i+1)}+u_2^{(i+1)})log\pi +(2-u_1^{(i+1)}-u_2^{(i+1)})log(1-\pi )\\ & +u_1^{(i+1)}[y_{1}logp+(1-y_1)log(1-p)]+u_2^{(i+1)}[y_{2}logp+(1-y_2)log(1-p)]\\ & +(1-u_1^{(i+1)})[y_{1}logq+(1-y_1)log(1-q)+(1-u_2^{(i+1)})[y_{2}logq+(1-y_2)log(1-q)]\end{array}$$对上述的结果进行分析，可以推测出进行$n$次数随机试验时，其$Q$函数可以表示为：$$Q(\theta ,{\theta }^{(i)})=\sum _{j=1}^n[u_j^{(i+1)}[log\pi +y_{j}logp+(1-y_j)log(1-p)]+(1-u_j^{(i+1)})[log(1-\pi )+y_{j}logq+(1-y_j)log(1-q)]]$$
• 接着，为了得到新一轮的参数${\theta }^{(i+1)}$, 需要使用$Q(\theta ,{\theta }^{(i)})$函数依次对$\pi$,$p$和$q$求偏导。具体如下：
  对$\pi$求偏导得到如下公式：$$\frac{\partial Q}{\partial \pi }=\sum _{j=1}^n[u_j^{(i+1)}\cdot \frac{1}{\pi }-\frac{1-u_j^{(i+1)}}{1-\pi }]$$令该计算式为0，则得到：$${\pi }^{(i+1)}=\frac{1}{n}\cdot \sum _{j=1}^n{u}_j^{(i+1)}$$对$p$求偏导得到如下公式：$$\frac{\partial Q}{\partial p}=\sum _{j=1}^n{u}_j^{(i+1)}(\frac{y_j}{p}+\frac{y_j-1}{1-p})$$令该计算式为0，得到:$$p^{(i+1)}=\frac{{\sum }_{j=1}^n{u}_j^{(i+1)}\cdot y_j}{{\sum }_{j=1}^n{u}_j^{(i+1)}}$$对$q$求偏导得到如下公式：$$\frac{\partial Q}{\partial q}=\sum _{j=1}^n(1-u_j^{(i+1)})(\frac{y_j}{q}+\frac{y_j-1}{1-q})$$同理，可以得到:$$q^{(i+1)}=\frac{{\sum }_{j=1}^n(1-u_j^{(i+1)})y_j}{{\sum }_{j=1}^n(1-u_j^{(i+1)})}$$

至此，三硬币模型的EM算法推导完毕。

参考资料

1. https://blog.csdn.net/weixin_41566471/article/details/106219019
2. 《统计学习方法》