机器学习入门：第十章 EM算法

引子

做完了课件间操，给了一项任务，测量同学的身高，整个学校的同学太多了，我们准备采用随机抽样的形式，留下了学号尾随为1的同学来估计学生$(x_1,x_2,...,x_n)$的身高情况。从以往的经验，身高是符合高斯分布的，概率密度函数为：

$$f(x)=\frac{1}{\delta \sqrt{2\pi }}{e}^{-\frac{(x-\mu {)}^2}{2{\delta }^2}}$$

其中，参数$\mu$和$\delta$是我们通过实验需要估计的参数。

这是很常见的概率参数的估计，用极大似然很容易可以得到结果：

$$\hat{\mu }=\frac{1}{n}\sum _{i=1}^n{x}_i=\overline{x},{\hat{\delta }}^2=\frac{1}{n}\sum _{i=1}^n(x_i-\overline{x}{)}^2$$

很容易就能得到学生身高的概率密度分布。
随着我们对身高研究的深入，我们知道，男女同学的身高是有区别的，我们需要了解男女身高的分布，即男学生身高的概率密度函数，女学生的概率密度函数，当然，男女学生的概率是同属于一种分布，只是参数不同罢了，一样是高斯分布。但目前我们是不能通过直接的计算就能得到的，因为我们在采集数据的时候，并没有把男女同学的身高分别开来。事后才知道对样本数据的采样时，属性的记录也同样重要，但我们还是要找到一种方式，能够把参数估计出来，也就是估计四个参数，${\mu }_{male}$和${\delta }_{male}$，${\mu }_{female}$和${\delta }_{female}$。
其实从分类的角度来看，就是要把男生和女生分开来，这也是为什么EM算法是一个重要的机器学习算法，他是一个有效的分类算法。

三硬币模型

再来看一个问题，通过这个问题，我们能把这一类问题，归到一种数学模型下。也就是大家熟悉的三硬币问题，三硬币模型。

假设有三枚硬币A、B、C，这些硬币正面出现的概率分别是 π 、p 、q 。进行如下掷硬币试验： 先掷A，如果A是正面则再掷B，如果A是反面则再掷C。对于B或C的结果，如果是正面则记为1，如果是反面则记为0。进行 NN 次独立重复实验，得到结果。现在只能观测到结果，不能观测到掷硬币的过程，估计模型参数 θ=(π,p,q) 。

在这个问题中，实验结果是可观测数据 Y=(y1,…,yn) ，硬币A的结果是不可观测数据 Z=(z1,…,zn) 且 z只有两种可能取值1和0。

这个模型简化了刚才测身高的例子，他是一种伯努利分布。而如果我们这样来看待测身高，其实这两个模型是一样样的，这里多了一个硬币A，他决定了B，C的出场，如果非要套到测量身高的例子，可以这样来看，我们在选择学生的时候，其实也是一个随机的过程，大概的流程跟那一个硬币A是一样的，他决定了男女的出场。他们的共同点，其实隐含着一个隐藏变量，也就是这里的硬币A。也是这个隐藏变量，使我们没办法用之前的极大似然方法，对概率进行估计。

EM算法

这里直接给出EM算法的步骤，证明的话，网上有许多的证明过程，可以自己查询。
上式完整表明了EM算法中的一步迭代中所需要的两个步骤：E-step，求期望；M-step，求极大值。EM算法的流程：

  输入：观测数据 Y，不可观测数据 Z；
  输出：参数 θ；

步骤：

• 给出参数初始化值 ${\theta }_0$

• E步：记 ${\theta }_0$为第 n次迭代后参数的估计值。在第 n+1 次迭代的E步，求期望（ Q 函数）

$$\begin{gathered} Q(\theta |thet{a}_n)=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}] \\ =\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)}) \end{gathered}$$

• M步：求 Q 函数的极大值点，来作为第 n+1 次迭代所得到的参数估计值 ${\theta }_{n+1}$

• 重复上面两步，直至达到停机条件。

三硬币模型的算法实现

来撸起袖子，我们始终没有忘记我们关键的问题，但是现在还不是时候，还是来看三硬币模型，用EM算法来求解这个模型。三枚硬币为A，B，C，出现正面的概率分别为$\pi ,q,p$，随机变量$z$是隐藏变量，也就是A到底是出现正面还是反面，对于第 j次试验，可以分两步求得，一种A是正面的情况，一种A是反面的情况。

$$\begin{gathered} P(y_i|\theta )=\sum _{z}P(y_i,z|\theta ) \\ =\sum _{z}P(z|\theta )P(y_i|z,\theta ) \\ =P(z=1|\theta )P(y_j|z=1,\theta )+P(z=0|\theta )P(y_j|z=0,\theta ) \\ =f(n)=\{\begin{array}{ll}\pi p+(1-\pi )q& if{y}_j=1\\ \pi (1-p)+(1-\pi )(1-q)& if{y}_j=1\end{array} \\ =\pi y^{y_i}(1-p{)}^{1-y_i}+(1-\pi )q^{y_i}(1-q{)}^{1-y_i} \end{gathered}$$

所以有

$$P(Y|\theta )=\prod _{j=1}^{N}P(y_i|\theta )=\prod (\pi p^{y_i}(1-p{)}^{1-y_i}+(1-\pi )q^{y_i}(1-q{)}^{1-y_i})$$

• E-step，求期望（Q函数）：

$$\begin{gathered} Q(\theta |thet{a}_n)=\sum _{z}P(z|Y,{\theta }_n)lnP(Y,z|\theta ) \\ =\sum _{j=1}^N(\sum _{z}P(z|y_i,{\theta }_n)lnP(y_i,z|\theta )) \\ =\sum _{j=1}^N(P(z=1|y_i,{\theta }_n)lnP(y_i,z=1|{\theta }_n)lnP(y_i,z=0|\theta ) \end{gathered}$$

先求$P(z|y_i,{\theta }_n)$，

因为他是伯努利分布，概率分布可以简化为：

$$P(z|y_i,{\theta }_n)=P(z|y_i,{\theta }_n)=\{\begin{array}{l}\frac{{\pi }_n{p}_n^{y_j}(1-p_n{)}^{1-y_j}}{{\pi }_n{p}_n^{y_j}(1-p_n{)}^{1-y_j}+(1-{\pi }_n)q_n^{y_j}(1-q_n{)}^{1-y_j}}=\mu ifz=11-{\pi }_{j,n},ifz=0\end{array}$$

因此 $P(y_i,z|{\theta }_n)=P(z|\theta )P(y_j|z,\theta )$

$$P(y_i,z|{\theta }_n)\{\begin{array}{l}\pi p^{y_j}(1-p{)}^{1-y_j}ifz=1\\ (1-\pi )q^{y_j}(1-q{)}^{1-y_j}ifz=0\end{array}$$

因此，Q函数表达式为：

$$Q(\theta |{\theta }_n)=\sum _{j=1}^N({\mu }_{j,n}ln[\pi p^{y_j}(1-p{)}^{1-y_j}]+(1-{\mu }_{j,n})ln[\pi q^{y_j}(1-q{)}^{1-y_j}])$$

• M-step，求解Q函数的极大值，来就求解$\pi$，p，q：
  要求$\pi$，就对$\pi$求导，令他等于0：

$$\begin{gathered} \frac{\partial Q(\theta |{\theta }_n)}{\partial \pi }=\sum _{j=1}^N(\frac{{\mu }_{j,n}ln[\pi p^{y_j}(1-p{)}^{1-y_j}]+(1-{\mu }_{j,n})ln[\pi q^{y_j}(1-q{)}^{1-y_j}]}{\partial \pi }) \\ =\sum _{j=1}^N({\mu }_{j,n}\frac{p^{y_j}(1-p{)}^{1-y_j}}{\pi p^{y_j}(1-p{)}^{1-y_j}}+(1-{\mu }_{j,n})\frac{-q^{y_j}(1-q{)}^{1-y_j}}{(1-\pi )q^{y_j}(1-q{)}^{1-y_j}}) \\ =\sum _{j=1}^N(\frac{{\mu }_{j,n}}{\pi (1-\pi )})=\frac{{\sum }_{j=1}^N{\mu }_{j,n}-n\pi }{\pi (1-\pi )} \end{gathered}$$

令$\frac{\partial Q(\theta |{\theta }_n)}{\partial \pi }=0$，可得：
$\pi =\frac{1}{n}{\sum }_{j=1}^N{\mu }_{j,n}$，可得${\pi }_{n+1}=\frac{1}{n}{\sum }_{j=1}^N{\mu }_{j,n}$
同样的过程，可以得到

$$p_{n+1}=\frac{{\sum }_{j=1}^N{\mu }_{j,n}{y}_i}{{\sum }_{j=1}^N{\mu }_{j,n}{y}_i+{\sum }_{j=1}^N{\mu }_{j,n}(1-y_i)},q_{n+1}=\frac{{\sum }_{j=1}^N(1-{\mu }_{j,n})y_i}{{\sum }_{j=1}^N(1-{\mu }_{j,n})y_i+{\sum }_{j=1}^N(1-{\mu }_{j,n})(1-y_i)}$$

注意这里公式的推导过程跟书上的过程是一样的，但是可能是表示方法的不同，表现的形式不同。这里$y_i$的意思为，第i次试验中，产生的正面的次数，$1-y_i$的意思是第i次试验中出现负面的次数。$y_i$和$1-y_i$分别代表一个数字，不可以分开。

公式和推导比较多，关键是下面的几个公式：

$$P(z|y_i,{\theta }_n)=\{\begin{array}{l}\frac{{\pi }_n{p}_n^{y_j}(1-p_n{)}^{1-y_j}}{{\pi }_n{p}_n^{y_j}(1-p_n{)}^{1-y_j}+(1-{\pi }_n)q_n^{y_j}(1-q_n{)}^{1-y_j}}=\mu ifz=11-{\pi }_{j,n},ifz=0\end{array}$$

${\pi }_{n+1}=\frac{1}{n}{\sum }_{j=1}^N{\mu }_{j,n}$
和

$$p_{n+1}=\frac{{\sum }_{j=1}^N{\mu }_{j,n}{y}_i}{{\sum }_{j=1}^N{\mu }_{j,n}{y}_i+{\sum }_{j=1}^N{\mu }_{j,n}(1-y_i)},q_{n+1}=\frac{{\sum }_{j=1}^N(1-{\mu }_{j,n})y_i}{{\sum }_{j=1}^N(1-{\mu }_{j,n})y_i+{\sum }_{j=1}^N(1-{\mu }_{j,n})(1-y_i)}$$

通过不断的迭代这个过程，p,q就能收敛。
具体的来看一个简单的例子，这个例子是三硬币模型。
假设有两枚硬币A、B，以相同的概率随机选择一个硬币，进行如下的掷硬币实验：共做5次实验，每次实验独立的掷十次，结果如图中a所示，例如某次实验产生了H、T、T、T、H、H、T、H、T、H，H代表证明朝上。a是在知道每次选择的是A还是B的情况下进行，b是在不知道选择的硬币情况下进行，问如何估计两个硬币正面出现的概率？

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-oDl3S0fn-1649987967520)(/images/kuxue/machine_learning/EM.png)]

第一种算法假设知道实验的时候知道哪个是A，哪个是B，属于极大似然法则，不在详细的介绍。因为图片的符号跟我们的稍微有点不同，为了统一改变一下，
看第二种算法，这里一开始假设

$$p_0=0.60,q_0=0.50,{\pi }_0=0.50$$

这里可以带入公式了：

$${\mu }_{j,n}=\frac{{\pi }_n{p}_n^{y_j}(1-p_n{)}^{1-y_j}}{{\pi }_n{p}_n^{y_j}(1-p_n{)}^{1-y_j}+(1-{\pi }_n)q_n^{y_j}(1-q_n{)}^{1-y_j}}$$
，这里第一次迭代n取为0，程序简化为：

$${\mu }_{j,0}=\frac{p_0^{y_j}(1-p_0{)}^{1-y_j}}{p_0^{y_j}(1-p_0{)}^{1-y_j}+q_0^{y_j}(1-q_0{)}^{1-y_j}}$$

，第一次试验，也就是产生了H、T、T、T、H、H、T、H、T、H的那次，这里的意思为，第i次试验中，产生的正面的次数，的意思是第i次试验中出现负面的次数，第一次试验中，，可以得到：

$${\mu }_{1,0}=\frac{p_0^5(1-p_0{)}^5}{p_0^5(1-p_0{)}^5+q_0^5(1-q_0{)}^5}=\frac{0.6^5(1-0.6{)}^5}{0.6^5(1-0.6{)}^5+0.5^5(1-0.5{)}^5}=0.45$$

$${\mu }_{2,0}=\frac{p_0^9(1-p_0{)}^1}{p_0^9(1-p_0{)}^1+q_0^9(1-q_0{)}^1}=\frac{0.6^9(1-0.6{)}^1}{0.6^9(1-0.6{)}^1+0.5^9(1-0.5{)}^1}=0.8$$

可以算到 ${\mu }_{3,0}=0.73,{\mu }_{4,0}=0.35,{\mu }_{5,0}=0.65$

$${\pi }_i=\frac{1}{5}\sum _{j=1}^5{\mu }_{j,n}=\frac{1}{5}(0.45+0.8+0.73+0.35+0.65)=0.596$$

$$\begin{gathered} p_i=\frac{{\sum }_{j=1}^5{\mu }_{j,0}{y}_i}{{\sum }_{j=1}^5{\mu }_{j,0}+{\sum }_{j=1}^5{\mu }_{j,0}(1-y_j)} \\ =\frac{0.45\bullet 5+0.8\bullet 9+0.73\bullet 8+0.35\bullet 4+0.65\bullet 7}{0.45\bullet 5+0.8\bullet 9+0.73\bullet 8+0.35\bullet 4+0.65\bullet 7+0.45\bullet 5+0.8\bullet 1+0.73\bullet 2+0.35\bullet 6+0.65\bullet 3} \\ =\frac{21.3}{21.3+8.6}=0.71 \end{gathered}$$

同样得到：

$q_1=0.58$

这样就完成了一次有效的迭代过程。继续重复上面的过程。

这里计算的是三硬币模型的情况，好像里我们测量身高的还有一定的距离，像前面所说的，测量身高的模型，跟三硬币模型的区别是概率模型的不同，三硬币属于伯努利分布，而身高模型属于高斯分布，推导的过程相似，这里直接给出计算公式。

高斯混合模型

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

输入：测量数据$y_1,y_2,...,y_N$，高斯混合模型；
输出：告示混合模型参数。

• 取参数的初始值开始迭代

• E步：依据当前模型参数，计算分模型k对观测数据$y_i$的响应度：

$${\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi (y_i|{\theta }_k)}{{\sum }_{k=1}^K{\alpha }_k\varphi (y_i|{\theta }_k)}$$

• M步：计算新一轮迭代的模型参数

$${\hat{\mu }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}{y}_i}{{\sum }_{k=1}^N}{\hat{\gamma }}_{jk},k=1,2,...,K$$

$${\hat{\sigma }}_k^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}(y_i-{\mu }_k{)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,...,K$$

$${\hat{\alpha }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N},k=1,2,...,K$$

• 重复(2)步和第(3)步，直到收敛。

用python代码来表示这两个过程，如下：
E步：

M步

重复的迭代这个两个过程，训练参数，两个高斯分布的参数就越来越明显了。

通过上边的迭代过程就解决了上面测量身高的过程，一开始测量身高的过程，主要是测量四个参数和，和，因为里面有一个隐含的参数，也就是男女的选取概率，这也是应用EM算法的特征了。
EM算法的全名是Expectation Maximization，中文名叫期望最大化算法。它是一个在含有隐变量的模型中常用的算法，在最大似然估计（MLE）和最大后验估计（MAP）中常用。在GMM、HMM、PCFG、IBM 5个对齐模型以及K-Means聚类方法中均有它的影子。

总结

EM算法，跟之前参数估计的方法最大的不同是多了一个隐藏的参数，这个隐藏的参数阻碍了我们利用极大似然估计。既然这个隐藏的参数那么难办，有没有可能绕过这个参数，进而利用已知的方法进行求解呢？那就假设这个隐藏的参数是已知的吧，我们通过慢慢的调这个隐藏的参数，只要这个隐藏的参数慢慢的达到正确的估值，那么我们所求的两个分布函数自然而然的可以用极大似然的方法进行求解了。而求解这个隐藏参数的方法，变成了寻找一个值，使得他最符合我们的观察值，也就是期望最大了。
EM算法可以通过纯数学的方法进行证明，EM算法是具有收敛性的，这就是为何我们可以通过一次一次的迭代，能够慢慢的靠近真值的理论依据了。

参考文献：1.LDA数学八卦
2.EM算法学习笔记与三硬币模型推导
3.EM算法实例分析