（八）统计学习方法 | EM算法及其推广

---

1. EM算法

1.1 简介与定义

$\mathrm{E}\mathrm{M}$算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。$\mathrm{E}\mathrm{M}$算法的每次迭代由两步组成：$\mathrm{E}$步：求期望，$\mathrm{M}$步：求极大。首先根据一个例子引入$\mathrm{E}\mathrm{M}$算法。

例题 假设有$3$枚硬币，分别记作$\mathrm{A},\mathrm{B},\mathrm{C}$。这些硬币正面出现的概率分别是$\pi ,p,q$。进行如下掷硬币实验：先掷硬币$\mathrm{A}$，根据其结果选出硬币$\mathrm{B}$或$\mathrm{C}$，正面选择$\mathrm{B}$，反面选择$\mathrm{C}$；然后掷选出的硬币，掷硬币的结果出现正面，记为$1$，出现反面记为$0$。现独立地重复$n=10$次实验，观测结果为$1,1,0,1,0,0,1,0,1,1$。现在问如何估计三硬币正面出现的概率，即三硬币模型的参数。

解 三硬币模型可以写作：$$\begin{array}{rl}P(y|\theta )& =\sum _{z}P(y,z|\theta )=\sum _{z}P(z|\theta )P(y|z,\theta )\\ & =P(z=1|\theta )P(y|z=1,\theta )+P(z=0|\theta )P(y|z=0,\theta )\\ & =\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}\end{array}\text{\hspace{1em}(1)}$$

这里的$y$表示观测变量，表示一次实验观测的结果是$1$或$0$；随机变量$z$是隐变量，表示为观测到$\mathrm{A}$的结果；$\theta =(\pi ,p,q)$是模型参数。将观测数据表示为 Y = { Y 1 , Y 2 , . . . , Y n } T Y=\{Y_1,Y_2,...,Y_n\}^{\rm T} Y={Y1​,Y2​,...,Yn​}T，未观测数据表示为 Z = { Z 1 , Z 2 , . . . , Z n } T Z=\{Z_1,Z_2,...,Z_n\}^{\rm T} Z={Z1​,Z2​,...,Zn​}T。则观测数据的似然函数为：$$P(Y|\theta )=\sum _{Z}P(Z|\theta )P(Y|Z,\theta )\text{\hspace{1em}(2)}$$

即：$$P(Y|\theta )=\prod _{j=1}^n[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}]\text{\hspace{1em}(3)}$$

考虑求模型$\theta =(\pi ,p,q)$的极大似然估计，即：$$\hat{\theta }=\arg \underset{\theta }{\max }\log P(Y|\theta )$$

这个问题没有解析解，只有通过迭代的方法求解。$\mathrm{E}\mathrm{M}$算法就是可以用求解这个问题的一种迭代算法。$\mathrm{E}\mathrm{M}$算法首先选取参数的初值，记作${\theta }^{(0)}=({\pi }^{(0)},p^{(0)},q^{(0)})$，然后通过下面的步骤迭代计算估计值，直到收敛为止。第$i$次迭代的估计值为${\theta }^{(i)}=({\pi }^{(i)},p^{(i)},q^{(i)})$，$\mathrm{E}\mathrm{M}$算法第$i+1$次迭代如下。

E步： 计算在模型参数${\pi }^{(i)},p^{(i)},q^{(i)}$下观测数据$y_j$来自掷硬币$\mathrm{B}$的概率：$${\mu }_j^{(i+1)}=\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}}\text{\hspace{1em}(4)}$$

M步： 计算模型参数的新估计值：$${\pi }^{(i+1)}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}\text{\hspace{1em}(5)}$$

$$p^{(i+1)}=\frac{\sum _{j=1}^n{\mu }_j^{(i+1)}{y}_j}{\sum _{j=1}^n{\mu }_j^{(i+1)}}\text{\hspace{1em}(6)}$$

$$q^{(i+1)}=\frac{\sum _{j=1}^n(1-{\mu }_j^{(i+1)})y_j}{\sum _{j=1}^n(1-{\mu }_j^{(i+1)})}\text{\hspace{1em}(7)}$$

现在进行数值计算。假设模型的参数初始取值为：$${\pi }^{(0)}=0.5,p^{(0)}=0.5,q^{(0)}=0.5$$

由式（4），对于$y_i=1$和$y_i=0$均有${\mu }_j^{(1)}=0.5$。则更新的参数为：$${\pi }^{(1)}=0.5,p^{(1)}=0.6,q^{(1)}=0.6$$

下一次更新的参数为：$${\pi }^{(2)}=0.5,p^{(2)}=0.6,q^{(2)}=0.6$$

于是得到模型参数$\theta$的极大似然估计：$$\hat{\pi }=0.5,\hat{p}=0.6,\hat{q}=0.6$$

如果取初值为${\pi }^{(0)}=0.4,p^{(0)}=0.6,q^{(0)}=0.7$，那么得到的模型参数的极大似然估计是$\hat{\pi }=0.4064,\hat{p}=0.5368,\hat{q}=0.6432$。即，$\mathrm{E}\mathrm{M}$算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。

EM算法

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

输出 模型参数$\theta$。

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

（2）$\mathrm{E}$步，记${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计值，在第$i+1$次迭代的$\mathrm{E}$步，计算$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =E_Z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}]\\ & =\sum _Z\log P(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})\end{array}\text{\hspace{1em}(8)}$$

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

（3）$\mathrm{M}$步，求使$Q(\theta ,{\theta }^{(i)})$极大化的$\theta$，确定第$i+1$次迭代的参数的估计值${\theta }^{(i+1)}$：$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\text{\hspace{1em}(9)}$$

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

式（8）的函数$Q(\theta ,{\theta }^{(i)})$是$\mathrm{E}\mathrm{M}$算法的核心，称为Q函数。

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

关于EM算法的几点：

• 参数的初值可以任意选择，但$\mathrm{E}\mathrm{M}$算法对初值敏感；
• $\mathrm{E}$步求$Q(\theta ,{\theta }^{(i)})$。$Q$函数中$Z$是未观测数据，$Y$是观测数据。$Q(\theta ,{\theta }^{(i)})$的第一个元表示要极大化的参数，第二个元表示参数的当前估计值。每次迭代实际在求$Q$函数及其极大；
• $\mathrm{M}$步求$Q(\theta ,{\theta }^{(i)})$的极大化，得到${\theta }^{(i+1)}$，完成一次迭代${\theta }^{(i)}\to {\theta }^{(i+1)}$；
• 给出停止迭代的条件，一般是对较小的正数${\epsilon }_1,{\epsilon }_2$，若满足：$$||{\theta }^{(i+1)}-{\theta }^{(i)}||<{\epsilon }_1||||Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})||<{\epsilon }_2$$
  则停止迭代。

1.2 EM算法的收敛性

定理 设$P(Y|\theta )$为观测数据的似然函数，${\theta }^{(i)}(i=1,2,...)$为$\mathrm{E}\mathrm{M}$算法得到的参数估计序列，$P(Y|{\theta }^{(i)})$为对应的似然函数序列，则$P(Y|{\theta }^{(i)})(i=1,2,...)$是单调递增的，即：$$P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)})\text{\hspace{1em}(10)}$$

定理 设$L(\theta )=\log P(Y|\theta )$为观测数据的对数似然函数，${\theta }^{(i)}(i=1,2,...)$为$\mathrm{E}\mathrm{M}$算法得到的参数估计序列，$L({\theta }^{(i)})(i=1,2,...)$为对应的对数似然函数序列。则：

（1）如果$P(Y|\theta )$有上界，则$L({\theta }^{(i)})=\log P(Y|{\theta }^{(i)})$收敛到某一值$L^{\ast }$；

（2）在函数$Q(\theta ,{\theta }^{\prime })$与$L(\theta )$满足一定条件下，由$\mathrm{E}\mathrm{M}$算法得到的参数估计序列${\theta }^{(i)}$的收敛值${\theta }^{\ast }$是$L(\theta )$的稳定点。

上述定理只能保证参数估计序列收敛到对数似然函数的稳定点，不能保证收敛到极大值点。所以在实际应用中，初值的选择变得非常重要，常用的办法是选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。

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

1.3.1 高斯混合模型

高斯混合模型 高斯混合模型是指具有如下形式的概率分布模型：$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)\text{\hspace{1em}(11)}$$

其中${\alpha }_k$是系数，${\alpha }_k\ge 0$，$\sum _{k=1}^K{\alpha }_k=1$；$\varphi (y|{\theta }_k)$是高斯分布密度，${\theta }_k=({\mu }_k,{\sigma }_k^2)$：$$\varphi (y|{\theta }_k)=\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)\text{\hspace{1em}(12)}$$

称为第$k$个分模型。

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

假设观测数据$y_1,y_2,...,y_N$由高斯混合模型生成：$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)\text{\hspace{1em}(13)}$$

其中，$\theta =({\alpha }_1,{\alpha }_2,...,{\alpha }_K;{\theta }_1,{\theta }_2,...,{\theta }_K)$。使用$\mathrm{E}\mathrm{M}$算法估计高斯混合模型的参数$\theta$。

明确隐变量，写出完全数据的对数似然函数 设想观测数据$y_i,j=1,2,...,N$是这样产生的：首先依概率${\alpha }_k$选择第$k$个高斯分布模型$\varphi (y|{\theta }_k)$，然后依第$k$个分模型的概率分布$\varphi (y|{\theta }_k)$生成观测数据$y_j$。这是观测数据$y_i,j=1,2,...,N$是已知的；反映观测数据$y_j$来自第$k$个分模型的数据是未知的，$k=1,2,...,K$，以隐变量${\gamma }_{jk}$表示，其定义如下：
$${\gamma }_{jk}=\{\begin{array}{rl}1,& 第j个观测来自第k个分模型\\ 0,& 否则\end{array}j=1,2,...,N;k=1,2,...,K\text{\hspace{1em}(14)}$$

有了观测数据$y_j$及未观测数据${\gamma }_{jk}$，那么完全数据是：$$(y_j,{\gamma }_{j2},,...,{\gamma }_{jK}),j=1,2,...,N$$

于是，完全数据的对数似然函数为：$$\begin{array}{rl}P(y,\gamma |\theta )& =\prod _{j=1}^{N}P(y_j,{\gamma }_{j2},,...,{\gamma }_{jK}|\theta )\\ & =\prod _{k=1}^K\prod _{j=1}^N[{\alpha }_k\varphi (y_j|{\theta }_k){]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N[\varphi (y_j|{\theta }_k){]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left[\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)\right]}^{{\gamma }_{jk}}\end{array}$$

式中：$$n_k=\sum _{j=1}^N{\gamma }_{jk},\sum _{k=1}^K{n}_k=N$$

那么，完全数据的对数似然函数为：$$\log P(y,\gamma |\theta )=\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\right]\right\}$$

EM算法的E步：确定Q函数

$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =E[\log P(y,\gamma |\theta )y,{\theta }^{(i)}]\\ & =E\left\{\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\right]\right\}\right\}\\ & =\sum _{k=1}^K\left\{\sum _{j=1}^N(E_{{\gamma }_{jk}})\log {\alpha }_k+\sum _{j=1}^N(E_{{\gamma }_{jk}})\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\right]\right\}\end{array}\text{\hspace{1em}(15)}$$

这里需要计算$E({\gamma }_{jk}|y,\theta )$，记为${\hat{\gamma }}_{jk}$：$$\begin{array}{rl}{\hat{\gamma }}_{jk}& =E({\gamma }_{jk}|y,\theta )=P({\gamma }_{jk}=1|y,\theta )\\ & =\frac{P({\gamma }_{jk}=1,y_j|\theta )}{\sum _{k=1}^{K}P({\gamma }_{jk}=1,y_j\theta )}\\ & =\frac{P({\gamma }_{jk}=1,y_j|\theta )P({\gamma }_{jk}=1|\theta )}{\sum _{k=1}^{K}P({\gamma }_{jk}=1,y_j\theta )P({\gamma }_{jk}=1|\theta )}\\ & =\frac{{\alpha }_k\varphi (y_j|{\theta }_k)}{\sum _{k=1}^K{\alpha }_k\varphi (y_j|{\theta }_k)}\end{array}$$

${\hat{\gamma }}_{jk}$表示当前模型参数下第$j$个观测数据来自第$k$个分模型的概率，称为分模型$k$对观测数据$y_j$的响应度。则综合上式，最后得到的$Q$函数为：$$Q(\theta ,{\theta }^{(i)})=\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\hat{\gamma }}_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\right]\right\}\text{\hspace{1em}(16)}$$

确定EM算法的M步

迭代的$\mathrm{M}$步是求函数$Q(\theta ,{\theta }^{(i)})$对$\theta$的极大值，即求新一轮迭代的模型参数：$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$

用${\hat{\mu }}_k,{\hat{\sigma }}_k^2$和${\hat{\alpha }}_k$，$k=1,2,...,K$表示${\theta }^{(i+1)}$的各参数。求${\hat{\mu }}_k,{\hat{\sigma }}_k^2$只需将式（16）分别对${\mu }_k,{\sigma }_k^2$求偏导数并令其为$0$即可得到；求${\hat{\alpha }}_k$是在$\sum _{k=1}^K{\alpha }_k=1$条件下求偏导数并令其为$0$得到的。求得的结果如下：$${\hat{\mu }}_k=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{\sum _{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,...,K\text{\hspace{1em}(17)}$$

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

$${\hat{\alpha }}_k=\frac{n_k}{N}=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}}{N}\text{\hspace{1em}(19)}$$

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

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

输入 观测数据$y_1,y_2,...,y_N$，高斯混合模型；

输出 高斯混合模型参数。

（1）$\mathrm{E}$步，根据当前模型参数，计算分模型$k$对观测数据$y_j$的响应度：
$${\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi (y_j|{\theta }_k)}{\sum _{k=1}^K{\alpha }_k\varphi (y_j|{\theta }_k)},j=1,2,...,N;k=1,2,...,K$$

（2）$\mathrm{M}$步，计算新一轮迭代的模型参数：$${\hat{\mu }}_k=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{\sum _{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,...,K$$

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

$${\hat{\alpha }}_k=\frac{n_k}{N}=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}}{N}$$

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

---

2. EM算法的推广

2.1 F函数的极大-极大算法

F函数 假设隐变量数据$Z$的概率分布$\tilde{P}(Z)$，定义分布$\tilde{P}$与参数$\theta$的函数$F(\tilde{P},\theta )$如下：$$F(\tilde{P},\theta )=E_{\tilde{P}}[\log P(Y,Z|\theta )]+H(\tilde{P})\text{\hspace{1em}(20)}$$

称为$F$函数，式中$H(\tilde{P})=-E_{\tilde{P}}\log \tilde{P}(Z)$是分布$\tilde{P}(Z)$的熵。

引理 对于固定的$\theta$，存在唯一的分布${\tilde{P}}_{\theta }$极大化$F(\tilde{P},\theta )$，这时${\tilde{P}}_{\theta }$由以下式子给出：$${\tilde{P}}_{\theta }(Z)=P(Z|Y,\theta )\text{\hspace{1em}(21)}$$

并且${\tilde{P}}_{\theta }$随$\theta$连续变化。

引理 若${\tilde{P}}_{\theta }(Z)=P(Z|Y,\theta )$，则：$$F(\tilde{P},\theta )=\log P(Y|\theta )\text{\hspace{1em}(22)}$$

定理 设$L(\theta )=\log P(Y|\theta )$为观测数据的对数似然函数，${\theta }^{(i)}(i=1,2,...)$为$\mathrm{E}\mathrm{M}$算法得到的参数估计序列，函数$F(\tilde{P},\theta )$由式（20）定义。如果$F(\tilde{P},\theta )$在${\tilde{P}}^{\ast }$和${\theta }^{\ast }$有局部极大值，那么$L(\theta )$在${\theta }^{\ast }$也有局部极大值。类似地，如果$F(\tilde{P},\theta )$在${\tilde{P}}^{\ast }$和${\theta }^{\ast }$达到全局最大值，那么$L(\theta )$在${\theta }^{\ast }$也达到全局最大值。

定理 $\mathrm{E}\mathrm{M}$算法的一次迭代可由$F$函数的极大-极大算法实现。设${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计，${\tilde{P}}^{(i)}$为第$i$次迭代函数$\tilde{P}$的估计。在第$i+1$次迭代的两步为：

（1）对固定的${\theta }^{(i)}$，求${\tilde{P}}^{(i+1)}$使$F(\tilde{P},{\theta }^{(i)})$极大化；

（2）对固定的${\tilde{P}}^{(i+1)}$，求${\theta }^{(i+1)}$使$F({\tilde{P}}^{(i+1)},\theta )$极大化。

2.2 GEM算法

GEM算法1

输入 观测数据，$F$函数；

输出 模型参数。

（1）初始化参数${\theta }^{(0)}$，开始迭代；

（2）第$i+1$次迭代，第$1$步：记${\theta }^{(i)}$为参数$\theta$的估计值，${\tilde{P}}^{(i)}$为函数$\tilde{P}$的估计，求${\tilde{P}}^{(i+1)}$使$\tilde{P}$极大化$F(\tilde{P},{\theta }^{(i)})$；

（3）第$2$步：求${\theta }^{(i+1)}$使$F({\tilde{P}}^{(i+1)},\theta )$极大化；

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

GEM算法2

输入 观测数据，$F$函数；

输出 模型参数。

（1）初始化参数${\theta }^{(0)}$，开始迭代；

（2）第$i+1$次迭代，第$1$步：记${\theta }^{(i)}$为参数$\theta$的估计值，计算：$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =E_Z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}]\\ & =\sum _Z\log P(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})\end{array}$$

（3）第$2$步：求${\theta }^{(i+1)}$使：$$Q({\theta }^{(i+1)},{\theta }^{(i)})>Q({\theta }^{(i)},{\theta }^{(i)})$$

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

GEM算法3

输入 观测数据，$F$函数；

输出 模型参数。

（1）初始化参数${\theta }^{(0)}=({\theta }_1^{(0)},{\theta }_2^{(0)},...,{\theta }_d^{(0)})$，开始迭代；

（2）第$i+1$次迭代，第$1$步：记${\theta }^{(i)}=({\theta }_1^{(i)},{\theta }_2^{(i)},...,{\theta }_d^{(i)})$为参数$\theta =({\theta }_1,{\theta }_2,...,{\theta }_d)$的估计值，计算：$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =E_Z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}]\\ & =\sum _Z\log P(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})\end{array}$$

（3）第$2$步：进行$d$次条件极大化。首先，在${\theta }_2^{(i)},...,{\theta }_d^{(i)}$保持不变的条件下求使$Q(\theta ,{\theta }^{(i)})$达到极大的${\theta }_1^{(i+1)}$；然后在${\theta }_1={\theta }_1^{(i+1)},{\theta }_j={\theta }_j^{(i)},j=3,4,...,d$的条件下求使$Q(\theta ,{\theta }^{(i)})$达到极大的${\theta }_2^{(i+1)}$；如此继续，经过$d$次条件极大化，求得${\theta }^{(i+1)}=({\theta }_1^{(i+1)},{\theta }_2^{(i+1)},...,{\theta }_d^{(i+1)})$使得：$$Q({\theta }^{(i+1)},{\theta }^{(i)})>Q({\theta }^{(i)},{\theta }^{(i)})$$

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

---

3. Python实现EM算法

这里实现$\mathrm{E}\mathrm{M}$算法在高斯混合模型中的应用。首先，我们需要自定义伪造的满足高斯分布的输入数据集：

def load_data(mu0, sigma0, alpha0, mu1, sigma1, alpha1):
    # 使用高斯分布生成伪造数据集，0和1分别代表两个高斯分布的参数且数据集长度为1000
    len_data = 1000
    # 第一个高斯分布
    data0 = np.random.normal(mu0, sigma0, int(len_data * alpha0))
    # 第二个高斯分布
    data1 = np.random.normal(mu1, sigma1, int(len_data * alpha1))
    # 总数据集
    data_set = []
    data_set.extend(data0)
    data_set.extend(data1)
    # 返回打乱后的数据
    random.shuffle(data_set)
    return data_set

根据高斯函数的概率密度形式计算：

def cal_gauss(data_set_arr, mu, sigma):
    # 根据高斯函数的概率密度函数形式，式(11)
    result = 1 / (math.sqrt(2 * math.pi) * sigma) * \
             np.exp(-1 * (data_set_arr - mu) ** 2 / 2 * sigma ** 2)
    # 返回
    return result

$\mathrm{E}\mathrm{M}$算法的$\mathrm{E}$步：

def e_step(data_set_arr, alpha0, mu0, sigma0, alpha1, mu1, sigma1):
    # EM算法的E步，首先计算分模型k对观测数据yi的响应度
    gamma0 = alpha0 * cal_gauss(data_set_arr, mu0, sigma0)
    gamma1 = alpha1 * cal_gauss(data_set_arr, mu1, sigma1)
    # 根据公式计算各自的响应度
    su = gamma0 + gamma1
    gamma0 = gamma0 / su
    gamma1 = gamma1 / su
    # 返回
    return gamma0, gamma1

$\mathrm{E}\mathrm{M}$算法的$\mathrm{M}$步：

def m_step(data_set_arr, mu0, mu1, gamma0, gamma1):
    # 根据高斯混合模型参数估计的EM算法中的M步，计算新的mu值
    new_mu_0 = np.dot(gamma0, data_set_arr) / np.sum(gamma0)
    new_mu_1 = np.dot(gamma1, data_set_arr) / np.sum(gamma1)
    # 计算新的sigma值
    new_sigma_0 = math.sqrt(np.dot(gamma0, (data_set_arr - mu0) ** 2) / np.sum(gamma0))
    new_sigma_1 = math.sqrt(np.dot(gamma1, (data_set_arr - mu1) ** 2) / np.sum(gamma1))
    # 计算新的alpha值
    new_alpha_0 = np.sum(gamma0) / len(gamma0)
    new_alpha_1 = np.sum(gamma1) / len(gamma1)
    # 返回
    return new_mu_0, new_mu_1, new_sigma_0, new_sigma_1, new_alpha_0, new_alpha_1

模型训练并返回预测结果：

def model_train(data_set, it=500):
    # 转换数据格式
    data_set = np.array(data_set)
    # 给定初始值
    alpha0 = 0.5
    mu0 = 0
    sigma0 = 1
    alpha1 = 0.5
    mu1 = 1
    sigma1 = 1
    # 开始迭代
    s = 0
    while s < it:
        s += 1
        # E步，计算响应度
        gamma0, gamma1 = e_step(data_set, alpha0, mu0, sigma0, alpha1, mu1, sigma1)
        # M步，更新参数
        mu0, mu1, sigma0, sigma1, alpha0, alpha1 = m_step(data_set, mu0, mu1, gamma0, gamma1)
    # 返回迭代后的结果
    return alpha0, mu0, sigma0, alpha1, mu1, sigma1

---

4. EM算法总结

$\mathrm{E}\mathrm{M}$算法是含有隐变量的概率模型极大似然估计或极大后验概率估计迭代算法。含有隐变量的概率模型的数据表示为$P(Y,Z|\theta )$，其中$Y$是观测变量的数据，$Z$是隐变量的数据，$\theta$是模型参数。$\mathrm{E}\mathrm{M}$算法的每次迭代包括两步：$\mathrm{E}$步：求期望，$\mathrm{M}$步：求极大。$\mathrm{E}\mathrm{M}$算法的应用及其广泛，高斯混合模型的参数估计是$\mathrm{E}\mathrm{M}$算法的一个重要应用，后面介绍的隐马尔可夫的无监督学习也是$\mathrm{E}\mathrm{M}$算法的一个重要应用。

---

参考

1. 统计学习方法/李航著。—2版。—北京：清华大学出版社，2019（2019.6重印）.
2. https://github.com/Dod-o/Statistical-Learning-Method_Code（EM算法代码）.

---

完

---