《李航 统计学习方法》学习笔记——第九章EM算法及其推广

EM算法及其推广

• 9.1 EM算法

• • 9.1.1 EM算法介绍
  • 9.1.2 EM算法推导

• 9.2 EM算法的收敛性

• 9.3 EM算法在高斯混合模型学习下的应用

• • 9.3.1 高斯混合模型GMM
  • 9.3.2 高斯混合模型参数估计的EM算法

• 9.4 EM算法的推广

• • 9.4.1 F函数的极大-极大算法
  • 9.4.2 GEM算法

• 习题9.1python实现

• 习题9.2

• 习题9.3 python实现

• 习题9.4 混合贝叶斯模型

• 参考

9.1 EM算法

适用于无监督学习 。训练数据只有输入没有对应的输出，可以认为无监督学习的训练数据是联合概率分布产生的数据。

应用场景：引自博主
如果我们要求解的是一个混合模型，只知道混合模型中各个类的分布模型（譬如都是高斯分布）和对应的采样数据，而不知道这些采样数据分别来源于哪一类（隐变量 ），那这时候就可以借鉴EM算法。EM算法可以用于解决数据缺失的参数估计问题（隐变量的存在实际上就是数据缺失问题，缺失了各个样本来源于哪一类的记录）

9.1.1 EM算法介绍

EM算法：含有隐变量的概率模型参数的极大似然估计法或极大后验概率估计法。每次迭代分为E步（极大期望）和M步（求极大）。

三硬币模型

在这个实例中，抛硬币A的结果（这里记为Z ）是无法观测的，所以该结果称之为隐变量 。

在该例子中，也可以称（Y,Z）为完全数据 ，Y为不完全数据

设随意变量Y为观测到的抛硬币的结果（即观测变量），\theta=(\pi,p,q)是模型参数，则可以得到模型的似然函数(因为这里是离散变量，所以是求和，连续变量则为求积分)：
P(Y|\theta)=\sum\limits_ZP(Z|\theta)P(Y|Z,\theta)\\ =\prod\limits_{j=1}^n\left [ \pi p^{y_j}(1-p)^{1-y_j}+(1-\pi)q^{y_j}(1-q)^{1-y_j} \right ]
按照正常的求解过程，对参数\theta进行极大似然估计，即：
\theta = \arg\max\limits_{\theta}\log P(Y|\theta)
这个问题没有解析解，因为模型中既有未知变量Z也有模型变量 \theta。只能通过迭代计算参数的估计值。EM算法就是用于求解这种问题的迭代算法。

算法9.1 EM算法
输入：观测数据变量Y,隐变量数据Z,联合分布P(Y,Z|\theta),条件分布P(Z|Y,\theta)
输出：模型参数\theta
（1）选择参数的初值\theta^{(0)}，开始迭代；
（2）E步:记\theta^{(i)}为第i此迭代参数\theta的估计值，在第i+1此迭代的E步，计算
Q(\theta,\theta^{(i)})=E_Z\left [ \log(P(Y,Z|\theta)|Y,\theta^{(i)} \right]\\ = \sum_Z\log(P(Y,Z|\theta)P(Z|Y,\theta^{(i)})

上式子的思想主要为：由于Zi是隐变量，为了能够计算\theta是排除其他可变量，所以通过利用固定值 Z的期望代替变量 Zi求近似极大化，即 Z_i \rightarrow E(Z) Q函数的意义为在给定观测数据Y和当前参数\theta^{(i)}下对为观测数据Z的条件概率分布P(Z|Y,\theta^{(i)})的期望。Q中的第一个变元表示要极大化的参数，第二个变元表示当前该参数的估计值。每次迭代的最终目的为求Q函数的极大 。

（3）M步：求使得Q(\theta,\theta^{(i)})极大化的\theta,确定第i+1次的迭代的参数估计值\theta^{(i+1)}:
\theta^{(i+1)} =\arg\max_\theta Q(\theta,\theta^{(i)})

每一次迭代产生\theta^{(i+1)},会使得似然函数增大或达到局部极值。

（4）重复（2）（3）直到收敛。（较小的正数\varepsilon_1,\varepsilon_2,若满足：
||\theta^{(i+1)}-\theta^{(i)}||<\varepsilon_1\\ or\\ ||Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})|| < \varepsilon_2
则停止迭代。）

9.1.2 EM算法推导

模型目标：极大化观测数据Y关于参数\theta的对数似然函数,即极大化：
L(\theta)=\log P(Y|\theta)=\log\sum_Z P(Y,Z|\theta)\\ =\log(\sum_ZP(Y|Z,\theta)P(Z|\theta))
在极大化的过程中，因为涉及到隐变量Z未知，所以极大化求解困难。
所以EM的主要思想是迭代逐步近似极大化 L(\theta)。
假设新估计值\theta,上一次的估计值\theta^{(i)},为了每次迭代能够使得极大值尽可能增加，所以这里考虑每次迭代对数似然函数的增加量 ：
L(\theta)-L(\theta^{(i)})=\log(\sum_ZP(Y|Z,\theta)P(Z|\theta))-\log P(Y|\theta^{(i)})\\ =\log(\sum_ZP(Z|Y,\theta^{(i)})\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})})-\log P(Y|\theta^{(i)})

Jensen不等式：
\log\sum_j\lambda_jy_j\geqslant\sum_j\lambda_j\log y_i \\ \lambda_j \geqslant 0 \ \ and \ \ \sum_j\lambda_j=1

利用以上不等式可以变形为
L(\theta)-L(\theta^{(i)}) \geqslant \sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}-\log P(Y|\theta^{(i)})\\ =\sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}
把L(\theta^{(i)})移向到等号右侧，令：
B(\theta,\theta^{(i)}) = L(\theta^{(i)})\ + \ \sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}
得到下界L(\theta)的一个下界B(\theta,\theta^{(i)})：
L(\theta)\geqslant B(\theta,\theta^{(i)})
因此，任意一个使得下界增大的\theta,也可以使似然函数增大，所以为了似然函数尽可能大的增长，应选择\theta^{(i+1)}使得下界达到极大：
\theta^{(i+1)} =\arg\max_\theta B(\theta,\theta^{(i)})
在求解该式子中，省去一些不包含\theta变量得常数项，得：
\theta^{(i+1)} = \arg\max_\theta(L(\theta^{(i)})\ + \ \sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})})\\ =\arg\max_\theta(\sum_ZP(Z|Y,\theta^{(i)})\log{P(Y|Z,\theta)P(Z|\theta)}) 保留了P(Z|Y,\theta^{(i)})该常数项\\ =\arg\max_\theta(\sum_ZP(Z|Y,\theta^{(i)})\log(P(Y,Z|\theta))\\ =\arg\max_\theta Q(\theta,\theta^{(i)})
到这里完成了EM算法得一次迭代，即Q函数得极大化。方式：通过不断求解下界得极大化逼近求解对数似然函数极大化得算法 。

通过上图可知，在每一步迭代中，先固定\theta_i确保似然函数改变参数后依旧增长,调整\theta^{(i+1)}使得下界取到极值，下一次迭代中，固定\theta^{(i+1)}，寻找下一个参数使得下界取到极值，循环往复，直到收敛到最优值\hat{\theta}。但是也可以看出EM算法不能保证找到全局最优值。

9.2 EM算法的收敛性

定理9.1

这个定理表明，EM算法每次迭代求得最大化完全数据模型似然函数的参数\theta会使观测数据的似然函数增大，这是一个很好的性质，因为最终目的就是使得观测数据似然函数最大。
定理9.2 （收敛性的说明——局部最优）

表明EM的收敛性包含关于对数似然函数序列的收敛性与关于参数估计序列的收敛性两层意思；且参数估计序列收敛到对数似然函数序列的稳定点，不能保证收敛到极大值点，所以EM算法会收敛到局部的最大值 。因此：初值的选择相对重要 ，可尝试不同初值带入，并加以比较。

• 若L(\theta，\theta^{(i)})是凸函数，则可以收敛到全局最大值，这点和梯度下降算法相似。

9.3 EM算法在高斯混合模型学习下的应用

9.3.1 高斯混合模型GMM

高斯分布就是俗称的正态分布 。
高斯混合模型的具体含义（以一维为例）：
观测数据Y不来自于同一高斯分布，而是以\alpha_k的概率来自于参数\theta_k=(\mu_k,\sigma _k)不同的k个高斯分布模型\phi (y|\theta_k)。在这个过程之中，来自哪个高斯分布是无法观测，因此可以称该变量叫做隐变量，因此适用于EM求解。
定义9.2 高斯混合模型
概率分布模型为：
P(y|\theta)=\sum\limits_{k=1}^K\alpha_k\phi (y|\theta_k)
其中：
\alpha_k \geqslant 0\\ \sum_{k=1}^K\alpha_k = 1\\ \phi (y|\theta_k) = \frac{1}{\sqrt{2\pi}\sigma_k}\exp(-\frac{(y-\mu_k)^2}{2\sigma_k^2})

类比三硬币模型的形式： P(Y|\theta)=\sum\limits_ZP(Z|\theta)P(Y|Z,\theta)
其中P(Z|\theta)看作是\alpha,当参数确定时，为Z=Z_i的概率；P(Y|Z,\theta)看作是不同的高斯概率分布，当参数和Z=Z_i确定时，该分布也确定了。

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

1. 确定隐变量 ，写出完全数据的对数似然函数 。
   假定观测数据产生如下：
   \alpha_k,\phi (y|\theta_k)\Rightarrow y_j
   但是y_i来自于第k个分模型中的哪一个无法观测，定义隐变量
   \gamma_{jk}=\left\{\begin{matrix} 1,\ \ \text{第j个观测来自于第k个分模型}& \\ 0,否则& \end{matrix}\right.\\ j=1,2,...,N;\ \ \ k=1,2,...,K

例如y_1来自于第2个分模型，则隐变量为 \gamma_{1,2}=(0,1,0,...,0)

由此得到的完全数据为：
(y_i,\gamma_{j1},\gamma_{j2}...,\gamma_{jK}) , \ \ j=1,2,...,N
由此得到完全函数的对数似然函数：
P(y,\gamma|\theta) = \prod_{j=1}^NP(y_i,\gamma_{j1},\gamma_{j2}...,\gamma_{jK}|\theta)\\ =\prod_{k=1}^{K}\prod_{j=1}^N\left [ \alpha_k,\phi (y|\theta_k)\right ]^{\gamma_{jk}}\\ 令n_k=\sum\limits_{j=1}^N\gamma_{jk}表示来自k分模型中样本的个数,\ \sum_{k=1}^Kn_k=N\\ = \prod_{k=1}^K \alpha_k^{n_k}\prod\limits_{j=1}^N\left [\phi (y|\theta_k)\right ]^{\gamma_{jk}}\\ =\prod_{k=1}^K \alpha_k^{n_k}\prod\limits_{j=1}^N\left[ \frac{1}{\sqrt{2\pi}\sigma_k}\exp(-\frac{(y-\mu_k)^2}{2\sigma_k^2})\right]^{r_{jk}}
取对数，得对数似然函数：
\log P(y,\gamma|\theta) = \begin{Bmatrix} \sum_{k=1}^Kn_k\log\alpha_k+\sum_{j=1}^Nr_{jk}[ \log(\frac{1}{\sqrt{2\pi}})-\log\sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2] \end{Bmatrix}

2.E步：确定Q函数
Q(\theta,\theta^{(i)})=E_{r_{jk}}[\log P(y,\gamma|\theta)|y,\theta^{(i)}]\\ =E_{r_{jk}}\begin{Bmatrix} \sum_{k=1}^Kn_k\log\alpha_k+\sum_{j=1}^Nr_{jk}[ \log(\frac{1}{\sqrt{2\pi}})-\log\sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2] \end{Bmatrix}\\ \text{} \\ \text用E(r_{jk})取代r_{jk}\\ \\ = \sum_{k=1}^K\begin{Bmatrix} \sum_{j=1}^N(E(r_{jk}))\log\alpha_k+\sum_{j=1}^NE(r_{jk})[ \log(\frac{1}{\sqrt{2\pi}})-\log\sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2]\end{Bmatrix}
这里涉及到E(r_{jk}|y,\theta)得求解：
\hat{\gamma_{jk}}=E({r_{jk}}|y,\theta)=P(\gamma_{jk}=1|y,\theta)\\ = \frac{P(\gamma_{jk}=1,y_j|\theta)}{\sum_{k=1}^KP(\gamma_{jk}=1,y_j|\theta)}\\ =\frac {P(y_j|\gamma_{jk}=1,\theta)\color{Red}P(\gamma_{jk}=1|\theta)}{\sum_{k=1}^KP(y_j|\gamma_{jk}=1,\theta){\color{Red}P(\gamma_{jk}=1|\theta)}}\\ =\frac{ \alpha_k\phi (y|\theta_k)}{\sum_{k=1}^K\alpha_k\phi (y|\theta_k)}
这里称\hat{\gamma_{jk}}是当前模型参数下第j个观测数据来自第k个分模型的概率，称为分模型对观测数据y_j得响应度。
将得到的E_{\gamma_{jk}}带入Q函数中(带入过程中n_k也相应的变成\sum_{j=1}^NE(r_{jk}))，得到：
Q(\theta,\theta^{(i)})=\sum_{k=1}^K\begin{Bmatrix} {\color{Red}n_k}\log\alpha_k+\sum_{j=1}^N{\color{Red}\hat{\gamma_{jk}}}[ \log(\frac{1}{\sqrt{2\pi}})-\log\sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2]\end{Bmatrix}

3. M步：极大化Q函数
   求函数Q对\theta的最大值作为新一轮迭代模型的参数：
   \theta^{(i+1)}=\arg\max Q(\theta,\theta^{(i)})
   模型参数\theta^{(i+1)}包括(\alpha_k,\mu_k,\sigma_k^2),通过极大似然估计分别对\alpha_k,\mu_k,\sigma_k^2求偏导数令其为0,其中\hat{\alpha_k}是在\sum_{k=1}^K\alpha_k=1限制条件下求解的。通过该过程求得的解如下：
   \hat{\mu_k}=\frac{\sum_{j=1}^N\hat{\gamma_{jk}}y_j}{\sum_{j=1}^N\hat{\gamma_{jk}}} \ \ \sigma_k^2=\frac{\sum_{j=1}^N\hat{\gamma_{jk}}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat{\gamma_{jk}}} \ \ \alpha_k=\frac{n_k}{N}=\frac{\sum_{j=1}^N\hat{\gamma_{jk}}}{N}
   到此一次EM迭代结束，重复以上步骤，直到实现停止条件（对数似然函数不明显变化或参数不明显变化。

算法9.2 高斯混合模型参数估计的EM算法
输入：观测数据Y，高斯混合模型；
输出：高斯混合模型参数。
（1）去参数的初始值开始迭代（由于EM算法对初始值敏感，可以通过赋值不同的初始值比较模型 ）
（2）E步，根据当前的模型参数，计算分模型k对观测数据的响应度 。
\hat{\gamma_{jk}} = \frac{ \alpha_k\phi (y|\theta_k)}{\sum_{k=1}^K\alpha_k\phi (y|\theta_k)}\ \ j=1,2,...,N;\ \ k=1,2,...K
（3）计算新一轮的迭代模型参数\theta：
\hat{\mu_k}=\frac{\sum_{j=1}^N\hat{\gamma_{jk}}y_j}{\sum_{j=1}^N\hat{\gamma_{jk}}} \ \ \sigma_k^2=\frac{\sum_{j=1}^N\hat{\gamma_{jk}}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat{\gamma_{jk}}} \ \ \alpha_k=\frac{n_k}{N}=\frac{\sum_{j=1}^N\hat{\gamma_{jk}}}{N}
（4）重复(2)(3)直到收敛。

9.4 EM算法的推广

EM算法的另外一种解释（F函数）与应用（GEM）。
首先通过图形的方式理解EM的算法原理。

图引自《The Elements of Statistical Learning》

红色曲线为P(Z|Y,\theta)，横坐标为隐变量，纵坐标为模型参数，等高线可以理解为Q函数或F函数，如图所示，EM算法可以理解，E步为了寻找最优的隐变量分布（即红色线上的点），M步为了寻找在当前隐变量分布下，使得Q函数或F函数极大的模型参数 ，并不断迭代，最终至收敛。

9.4.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})
其中H(\tilde{P})=-E_{\tilde{P}}\log\tilde{P}(Z)，是分布{\tilde{P}(Z)}的熵。

EM函数中分别通过E和M步骤求解隐函数分布和参数，而F函数将这两个步骤融合在一个函数中，通过固定其一，调整另一个的策略完成求解。

• 第一个极大 ：固定\theta，关于Z的分布求F极大值：
  
  证明：
  拉格朗日方法求解关于{\tilde{P}(Z)}的极大值问题：
  L=E_{\tilde{P}}[\log P(Y,Z|\theta)]\ -E_{\tilde{P}}\log\tilde{P}(Z) +\lambda[1-\sum_Z\tilde{P}(Z)]
  其中\sum_z\tilde{P}(Z)=1为约束条件。
  将其对\tilde{P}(Z)求导并令导数等于0求得：
  \lambda=\log P(Y,Z|\theta) - \log\tilde{P_\theta}(Z)+1
  变换该式子得：
  \frac{P(Y,Z|\theta)}{\tilde{P_\theta}(Z)}=e^{\lambda+1}
  通过该式子可以得知，因为P(Y,Z|\theta)是连续的，所以\tilde{P_\theta}(Z)也随\theta连续变化。
  又因为\sum_z\tilde{P}(Z)=1，所以可以得出：
  \tilde{P}_\theta(Z)=P(Z|Y,\theta)
  这个性质确保了第一个极大步骤的可行性，固定\theta，寻找最佳的隐变量分布时，目的是使得当前的F_\theta(\tilde{P},\theta)最大。该性质表明，通过该方式寻找到的最佳隐变量唯一且正好取在{P(Z|Y,\theta)}这个点上，结果很理想。

• 第二个极大 ：固定{\tilde{P}(Z)}，关于\theta求F极大值。
  
  证明：
  F(\tilde{P},\theta)=E_{\tilde{P}}[\log P(Y,Z|\theta)]\ +\ H(\tilde{P}) \\ =E_{\tilde{P}}[\log P(Y,Z|\theta)]\ -E_{\tilde{P}}\log\tilde{P}(Z)\\ =\tilde{P}_\theta(Z)\sum_Z\log P(Y,Z|\theta)\ -\tilde{P}(Z)\sum_Z\log\tilde{P}(Z)\\ 若\tilde{P}_\theta(Z)=P(Z|Y,\theta),带入上式\\ =P(Z|Y,\theta)\sum_Z\log P(Y,Z|\theta)\ -P(Z|Y,\theta)\sum_Z\log P(Z|Y,\theta)\\ =[\log P(Y,Z|\theta)\ - \log P(Z|Y,\theta)]\sum_ZP(Z|Y,\theta)\\ =[\log \frac {P(Y,Z|\theta)}{P(Z|Y,\theta)}]\sum_ZP(Z|Y,\theta)\\ =\log P(Y|\theta)\sum_ZP(Z|Y,\theta)\\ 因为\sum_ZP(Z|Y,\theta)=1\\ =\log P(Y|\theta)
  这个性质表明：当我们把\tilde{P}_\theta(Z)设为P(Z|Y,\theta)时，恰好目标函数\log P(Y|\theta)等于F函数，而且去除了隐变量，这为第二个极大提供了保证，关于\theta求极大。
  这步正是EM算法要解决的问题，该步骤等价于EM中的：
  \theta = \arg\max_\theta\log P(Y|\theta)

通过以上两个极大步骤的迭代求得得\tilde{P}^*,\theta^*与EM算法的关系
这个性质说明了，两个方法求得的最优解相同 。

F函数极大-极大算法具体的步骤：

证明过程
可以得知EM算法中极大化Q函数就等价于极大化F函数 。

通过该算法得知，F函数的极大-极大算法求得的参数估计序列与EM算法求得的参数估计序列相同 ，很好的解释了EM算法。

9.4.2 GEM算法

GEM：广义EM算法。
算法9.3 GEM算法1
输入：观测数据，F函数；
输出：模型参数。
（1）初始化参数\theta^{(0)},开始迭代；
（2）第（i+1）次迭代：
第一步：记\theta^{(i)}为参数\theta的估计值，\tilde{P}^{(i)}为函数\tilde{P}的估计，求\tilde{P}^{(i+1)}使\tilde{P}极大化F(\tilde{P},\theta^{(i)})；
第二步：求\theta^{(i+1)}使得F(\tilde{P}^{(i+1)},\theta)极大化。
（3）重复（2）的迭代，直到收敛。

GEM算法1中有时Q(\theta,\theta^{(i)})(或F(\tilde{P}^{(i+1)},\theta))过于复杂，直接求它的极大化很困难。GEM算法2和GEM算法3为了解决这个问题，不直接求使得Q函数达到极大的 \theta^{(i+1)}，而是找一个\theta^{(i+1)}使得Q(\theta^{(i+1)},(\theta^{(i)})>Q((\theta^{(i)},(\theta^{(i)})

算法9.4 GEM算法2 （模型参数维数为1维）
输入：观测数据，Q函数；
输出：模型参数。
（1）初始化参数\theta^{(0)},开始迭代；
（2）第（i+1）次迭代：
第一步：记\theta^{(i)}为参数\theta的估计值，计算：
Q(\theta,\theta^{(i)})=E_Z\left [ \log(P(Y,Z|\theta)|Y,\theta^{(i)} \right]\\ = \sum_Z\log(P(Y,Z|\theta)P(Z|Y,\theta^{(i)})
第二步：求\theta^{(i+1)}使得
Q(\theta^{(i+1)},(\theta^{(i)})>Q((\theta^{(i)},(\theta^{(i)})
（3）重复（2）的迭代，直到收敛。

算法9.5 GEM算法3 （模型参数维数为d维,d>=2）
基本思想：将EM算法的M步分解为d次条件极大化，每次只改变参数向量的一个分量，其余分量不变。
输入：观测数据，Q函数;
输出：模型参数。
（1）初始化参数\theta^{(0)}=(\theta_1^{(0)},\theta_2^{(0)},...,\theta_d^{(0)}),开始迭代；
（2）第（i+1）次迭代：

• 第一步：记\theta^{(i)}=(\theta_1^{(i)},\theta_2^{(i)},...,\theta_d^{(i)})为参数\theta的估计值，计算：
  Q(\theta,\theta^{(i)})=E_Z\left [ \log(P(Y,Z|\theta)|Y,\theta^{(i)} \right]\\ = \sum_Z\log(P(Y,Z|\theta)P(Z|Y,\theta^{(i)})

• 第二步：进行d次条件极大化
  首先，在\theta_2^{(i)},\theta_3^{(i)},...,\theta_d^{(i)}保持不变的条件下求使Q函数达到极大的\theta_1^{(i+1)};
  然后在\theta_1=\theta_1^{(i+1)}，\theta_j=\theta_j^{(i)}，j=3,4,...d的条件下求使Q函数达到极大的\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)})
  （3）重复（2）的迭代，直到收敛。

习题9.1python实现

    def cal_u(pi, p, q, xi):
    """
    u值计算
    :param pi: 下一次迭代开始的 pi
    :param p:  下一次迭代开始的 p
    :param q:  下一次迭代开始的 q
    :param xi: 观察数据第i个值，从0开始
    :return:
    """
    return pi * math.pow(p, xi) * math.pow(1 - p, 1 - xi) / \
           float(pi * math.pow(p, xi) * math.pow(1 - p, 1 - xi) +
                 (1 - pi) * math.pow(q, xi) * math.pow(1 - q, 1 - xi))
    
    def e_step(pi,p,q,x):
    """
        e步计算
    :param pi: 下一次迭代开始的 pi
    :param p:  下一次迭代开始的 p
    :param q:  下一次迭代开始的 q
    :param x: 观察数据
    :return:
    """
    return [cal_u(pi,p,q,xi) for xi in x]

    def m_step(u,x):
    """
     m步计算
    :param u:  m步计算的u
    :param x:  观察数据
    :return:
    """
    pi1=sum(u)/len(u)
    p1=sum([u[i]*x[i] for i in range(len(u))]) / sum(u)
    q1=sum([(1-u[i])*x[i] for i in range(len(u))]) / sum([1-u[i] for i in range(len(u))])
    return [pi1,p1,q1]

终止条件：

    def run(observed_x, start_pi, start_p, start_q, iter_num):
    """
    
    :param observed_x:  观察数据
    :param start_pi:  下一次迭代开始的pi $\pi$
    :param start_p:  下一次迭代开始的p
    :param start_q:  下一次迭代开始的q
    :param iter_num:  迭代次数
    :return:
    """
    for i in range(iter_num):
        u = e_step(start_pi, start_p, start_q, observed_x)
        print (i,[start_pi,start_p,start_q])
        if [start_pi,start_p,start_q]== m_step(u, observed_x):#收敛条件
            break
        else:
            [start_pi,start_p,start_q]= m_step(u, observed_x)

调用：

    # 观察数据
    x = [1, 1, 0, 1, 0, 0, 1, 0, 1, 1]
    # 初始化 pi，p q
    [pi, p, q] = [0.45, 0.55, 0.67]
    # 迭代计算
    run(x,pi,p,q,100)

习题9.2

证明引理9.2

证明：
F(\tilde{P},\theta)=E_{\tilde{P}}[\log P(Y,Z|\theta)]\ +\ H(\tilde{P}) \\ =E_{\tilde{P}}[\log P(Y,Z|\theta)]\ -E_{\tilde{P}}\log\tilde{P}(Z)\\ =\tilde{P}_\theta(Z)\sum_Z\log P(Y,Z|\theta)\ -\tilde{P}(Z)\sum_Z\log\tilde{P}(Z)\\ 若\tilde{P}_\theta(Z)=P(Z|Y,\theta),带入上式\\ =P(Z|Y,\theta)\sum_Z\log P(Y,Z|\theta)\ -P(Z|Y,\theta)\sum_Z\log P(Z|Y,\theta)\\ =[\log P(Y,Z|\theta)\ - \log P(Z|Y,\theta)]\sum_ZP(Z|Y,\theta)\\ =[\log \frac {P(Y,Z|\theta)}{P(Z|Y,\theta)}]\sum_ZP(Z|Y,\theta)\\ =\log P(Y|\theta)\sum_ZP(Z|Y,\theta)\\ 因为\sum_ZP(Z|Y,\theta)=1\\ =\log P(Y|\theta)

习题9.3 python实现

    import numpy as np
    from scipy.stats import norm
    
    y = np.array([-67, -48, 6, 8, 14, 16, 23, 24, 28, 29, 41, 49, 56, 60, 75])
    K = 2  # 分模型数
    N = 15  # 观测值个数
    
    # 参数初始化
    mu = np.array([0.5, 0.5])
    sg = np.array([1.0, 1.0]) * 10 #sigma
    al = np.array([0.5, 0.5])#alpha
    
    for i in range(10):
    gm = np.zeros((N, K))#隐变量gamma矩阵
    
    # E 步
    for j in range(N):
        for k in range(K):
            gm[j, k] = al[k] * norm(mu[k], sg[k]).pdf(y[j]) #计算gamma矩阵
        gm[j, :] /= sum(gm[j, :])
    
    # M 步
    mu2 = y.dot(gm) / sum(gm)
    al2 = sum(gm) / N
    sg2 = np.zeros((2,))
    sg2[0] = sum(gm[:, 0] * (y - mu[0]) ** 2) / sum(gm[:, 0])
    sg2[1] = sum(gm[:, 1] * (y - mu[1]) ** 2) / sum(gm[:, 1])
    
    if sum((mu - mu2) ** 2 + (sg - sg2) ** 2 + (al - al2) ** 2) < 0.01:#设置阈值，参数平方误差和小于0.01
        break
    
    mu = mu2
    sg = sg2
    al = al2
    
    print("迭代次数： ", i + 1, "\n", "mu = ", mu, "\n", "sg = ", sg, "\n", "al = ", al, "\n")

习题9.4 混合贝叶斯模型

首先拓展EM的一般形式：
这里提到的Q和\omega^{(i)}都指隐函数Z的一种概率分布。
具体的证明过程

按照一般形式，对混合贝叶斯模型进行求解。

参考

[1] 《统计学习导论》李航
[2] 《The Elements of Statistical Learning》
[3]
[4] https://www.zhihu.com/question/47830123
[5]
[6]
[7]

笔者刚刚入门学习机器学习，因为水平有限，李航老师的书对入门不是特别友好，还在生啃阶段，如果有错误还请之处。