EM(expectation and maximization)算法学习笔记

Prerequisite

• 好心情

模型

Driving Example

扔钢镚儿：
钢镚一号 钢镚二号
图片来自：哔哩哔哩

设：
一号正面和二号正面出现的概率分别为：$pq$
并且有奇数次投掷一号，偶数次投掷二号，反复十次得到结果：
$1,1,0,1,0,0,1,0,1,1$
现在试问用极大似然估计的方法求 $pq$：

解：
由题那么有：
似然函数为：$L(p,q)=P(\mathbf{x}|p,q)=p(x_1|p,q)\cdot p(x_2|p,q)\cdots p(x_{10}|p,q)$
即有：
$L(p,q)=p\ast q\ast (1-p)\ast q\ast (1-p)\ast (1-q)\ast p\ast (1-q)\ast p\ast q$
$=p^3\ast q^3\ast (1-p{)}^2\ast (1-q{)}^2$
若要 $maxL(p,q)$, 只需分别 $max(p^3\ast (1-p{)}^2)$ 和 $max(q^3\ast (1-q{)}^2)$
显然 $p=q=0.6$

上面的例子很简单，但相传有一好事者，他是这样子做的：
好事者扔钢镚儿：
钢镚一号钢镚二号钢镚三号
设：
一号正面，二号正面，三号正面出现的概率分别为：$\pi ,p,q$
并且有，先投掷一号，然后根据它的正反，来决定投掷二号（一号为正）或者三号（一号为负）。反复十次投掷有：$1,1,0,1,0,0,1,0,1,1$ 。
好事者说：
假设只通过观测到掷硬币的结果，不能观测掷硬币的过程，问如何估计三个硬币出现正面的概率？
解：
比起上面的奇偶投掷法，这里无非就是把奇偶选择硬币的方式，换成了用另外一个硬币来做选择。但尴尬的是，好事者不让我们知道一号的结果，因此我们无法直接用$p(x_1|p,q)\cdot p(x_2|p,q)\cdots$的方法来建立似然函数了，因为我们没有一号硬币的信息。那么就让我们来建立一个三参数模型，引入一号硬币的信息。
$$\begin{array}{rl}P(y|\theta )& =\sum _{z}P(y,z|\theta )=\sum _{z}P(z|\theta )P(y|z,\theta )\\ & =\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}\end{array}$$
y是观测变量，表示一次观测结果是1或0，z是隐藏变量(hidden variable)，表示掷硬币一号的结果，这个是观测不到结果的，$\theta =(\pi ,p,q)$表示模型参数，将观测数据表示为$Y=(Y_1,Y_2,...,Y_n{)}^T$，未观测的数据表示为$Z=(Z_1,Z_2,...,Z_n{)}^T$，则观测函数的似然函数是：
$$\begin{array}{rl}P(Y|\theta )& =\sum _{Z}P(Z|\theta )P(Y|Z,\theta )\\ & or\\ & =\prod _{i=0}^{n}P(y_i|\theta )\\ & =\prod _{i=0}^n(\pi p^{y_i}(1-p{)}^{1-y_i}+(1-\pi )q^{y_i}(1-q{)}^{1-y_i})\end{array}$$
注意$Z=(Z_1,Z_2,...,Z_n{)}^T$是向量，不是取其中一个值！！
那么对$\theta =(\pi ,p,q)$求极大似然估计有：
$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}logP(Y|\theta )$
就这？

是不是lingo, mathlab启动就完事了？！！
然而这个问题没有解析解，所以我们就只能另寻他路。

EM算法：

E:expectation M:maximization

$L(\theta )-L({\theta }^{(i)})$:

上文我们提到有似然函数：
$$\begin{array}{rl}L(\theta )& =logP(Y|\theta )=log\sum _{Z}P(Y,Z|\theta )\\ & =log(\sum _{Z}P(Y|Z,\theta )P(Z|\theta ))\end{array}$$
假设我们有随机初始值：${\theta }^{(i)}$, 现在我们的任务是找到一个$\theta$有：$L(\theta )>L({\theta }^{(i)})$。那么也即是有：
$$L(\theta )-L({\theta }^{(i)})=log(\sum _{Z}P(Y|Z,\theta )P(Z|\theta ))-logP(Y|{\theta }^{(i)})>0$$
现在我们希望对这个式子进行变形，从而找到$L(\theta )$的下界，如果有下界大于$L({\theta }^{(i)})$，那么就达到了差值大于零的目的。并且我们尽量是下界大，那么肯定$L(\theta )$也会尽量变大。但要注意，下界最大时，$L(\theta )$并不一定就取得了最大。下文会提到这一点。
$$\begin{array}{rl}L(\theta )-L({\theta }^{(i)})& =log(\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-logP(Y|{\theta }^{(i)})\\ & 由log凹函数性质有:log\sum _j{\lambda }_j{y}_j\ge \sum _j{\lambda }_{j}log{y}_j,{\lambda }_j\ge 0,\sum _j{\lambda }_j=1，则\\ & \ge \sum _{Z}P(Z|Y,{\theta }^{(i)}))log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-logP(Y|{\theta }^{(i)})\\ & Known:\sum _{Z}P(Z|Y,\theta )=1,then\\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{logP(Y|{\theta }^{(i)})}{\sum _{Z}P(Z|Y,{\theta }^{(i)})}\\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}\end{array}$$

$B(\theta ,{\theta }^{(i)})$

令：
$$B(\theta ,{\theta }^{(i)})\widehat{=}L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}$$
则：
$$L(\theta )\ge B(\theta ,{\theta }^{(i)})$$
当令$\theta ={\theta }^{(i)}$时，我们有$L({\theta }^{(i)})=B({\theta }^{(i)},{\theta }^{(i)})$
上文说到极大$B({\theta }^{(i)},{\theta }^{(i)})$，那么也会尽量增大$L(\theta )$。
那么我们的任务就变成了：
$${\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}B(\theta ,{\theta }^{(i)})$$
去除一些关于$\theta$的常数：
$$\begin{array}{rl}{\theta }^{(i+1)}=& \underset{\theta }{\mathrm{argmax}}(L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})})\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Y|Z,\theta )P(Z|\theta ))-\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Y|{\theta }^{(i)})P(Z|Y,{\theta }^{(i)}))\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Y|Z,\theta )P(Z|\theta ))\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )\\ & =\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^{(i)})--划重点：Q函数！\end{array}$$

最优性：

回到EM算法能否找到最优解的问题上，看下图：

图片来源：李航-统计学习方法

可以发现，虽然说我们能使$B(\theta ,{\theta }^{(i)})$最大，但是由于推导中的不等式，此时$L({\theta }^{(i+1)})$不一定最大。因而EM算法不能保证找到全局最优解。

Q函数

定义：

完全数据的对称似然函数$logP(Y,Z|\theta )$关于在给定观测数据Y和当前参数${\theta }^{(i)}$下对未观测数据Z的条件概率分布$P(Z|Y,{\theta }^{(i)})$的期望称为Q函数。（摘自李航统计学习）
即：
$$Q(\theta ,{\theta }^{(i)})=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]$$

伪代码：

$初始化\theta ,需要注意一点是EM算法对初始值比较敏感$
$fori=1tok:$
$E步：Q(\theta ,{\theta }^{(i)})=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]$
$M步:{\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^{(i)})$
$until||{\theta }^{(i+1)}-{\theta }^{(i)}||<{ϵ}_{1}or||Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})||<{ϵ}_2$
通过上面的流程也可以看得出来为什么叫EM算法。

收敛性：

please jump: https://blog.csdn.net/sinat_26376671/article/details/44313885?depth_1-utm_source=distribute.pc_relevant_right.none-task-blog-BlogCommendFromBaidu-3&utm_source=distribute.pc_relevant_right.none-task-blog-BlogCommendFromBaidu-3

高斯混合模型的参数估计

定义：

高斯混合模型是指有如下形式的概率分布模型：
$$\begin{gathered} P(y|\theta )=\sum _{k=1}^K{a}_k\varphi (y|{\theta }_k) \\ {a}_k>0,\sum a_k=1,{\theta }_k=(u_k,{\delta }_k^2) \\ \varphi (y|{\theta }_k)=\frac{1}{\sqrt{2\pi }{\delta }_k}exp(-\frac{(y-{\mu }_k{)}^2}{2{\delta }_k^2}),称为第k个分模型 \end{gathered}$$
考虑：
让高斯混合模型生成样本$Y=(y_1,y_2,\cdots y_j{)}^T,j=1,2,3,\cdots ,N$，那么我们如何来理解这个生成过程呢？
由于$\sum {\alpha }_k=1$，我们可以认为，第 i 次生成时，根据${\alpha }_k$的相对大小来选择第K个模型来生成$y_i$。
这跟扔钢镚很像，先确定一个前置量（无法观测），再用相应的后置量来生成样本。因此，我们的隐变量（所谓的前置量）也就确定了。
在这里我们设隐变量为：$Z=(z_1,z_2,\cdots z_j{)}^T,j=1,2,3,\cdots ,N$，它的每一个分量表征的是每一次采样所用的是哪一个分模型$\varphi (y|{\theta }_k)$,即 $z_j=1,2,3,\cdots ,k$

E步：

由上文，我们知道了，观察变量$Y$，和隐变量$Z$。又由极大似然估计的EM算法，我们有：
$$\begin{array}{rl}{\theta }^{(i+1)}& =\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^{(i)})\\ & =\underset{\theta }{\mathrm{argmax}}{E}_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )\\ 由样本间的独立性有：& \\ & =\underset{\theta }{\mathrm{argmax}}\sum _{z_1,z_2,\cdots ,z_j}P(z_1|y_1,{\theta }^{(i)})\cdots P(z_j|y_j,{\theta }^{(i)})logP(y_1,z_1|\theta )\cdots P(y_j,z_j|\theta )\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{z_1,z_2,\cdots ,z_j}\prod _{j=1}^{N}P(z_j|y_j,{\theta }^{(i)})\sum _{j=1}^{N}logP(y_j,z_j|\theta )\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{z_1,z_2,\cdots ,z_j}\prod _{j=1}^{N}P(z_j|y_j,{\theta }^{(i)})(logP(y_1,z_1|\theta )+\sum _{j=2}^{N}log(y_j,z_j|\theta ))\\ 取前一加项有：& \\ & =\sum _{z_1,z_2,\cdots ,z_j}\prod _{j=1}^{N}P(z_j|y_j,{\theta }^{(i)})logP(y_1,z_1|\theta )\\ & =\sum _{z_1,z_2,\cdots ,z_j}logP(y_1,z_1|\theta )P(z_1|y_1,{\theta }^{(i)})\prod _{j=2}^{N}P(z_j|y_j,{\theta }^{(i)})\\ & =\sum _{z_1}logP(y_1,z_1|\theta )P(z_1|y_1,{\theta }^{(i)})\sum _{z_2,\cdots ,z_j}\prod _{j=2}^{N}P(z_j|y_j,{\theta }^{(i)})\\ & =\sum _{z_1}logP(y_1,z_1|\theta )P(z_1|y_1,{\theta }^{(i)})\sum _{z_2}\sum _{z_3,\cdots ,z_j}P(z_2|y_2,{\theta }^{(i)})\prod _{j=3}^{N}P(z_j|y_j,{\theta }^{(i)})\\ & =\sum _{z_1}logP(y_1,z_1|\theta )P(z_1|y_1,{\theta }^{(i)})\sum _{z_2}P(z_2|y_2,{\theta }^{(i)})\sum _{z_3,\cdots ,z_j}\prod _{j=3}^{N}P(z_j|y_j,{\theta }^{(i)})\\ 可以发现\sum _{z_2}P(z_2|y_2,{\theta }^{(i)})=1，所以有：& \\ & =\sum _{z_1}logP(y_1,z_1|\theta )P(z_1|y_1,{\theta }^{(i)})\\ 综上讨论有：& \\ & \sum _{z_1,z_2,\cdots ,z_j}\prod _{j=1}^{N}P(z_j|y_j,{\theta }^{(i)})logP(y_1,z_1|\theta )\\ & =\sum _{z_1}logP(y_1,z_1|\theta )P(z_1|y_1,{\theta }^{(i)})\\ 即是说有：& \\ & \sum _{z_1,z_2,\cdots ,z_j}\prod _{j=1}^{N}P(z_j|y_j,{\theta }^{(i)})logP(y_j,z_j|\theta )\\ & =\sum _{z_j}logP(y_j,z_j|\theta )P(z_j|y_j,{\theta }^{(i)})\\ 那么：& \\ & \underset{\theta }{\mathrm{argmax}}\sum _{z_1,z_2,\cdots ,z_j}\prod _{j=1}^{N}P(z_j|y_j,{\theta }^{(i)})(logP(y_1,z_1|\theta )+\sum _{j=2}^{N}log(y_j,z_j|\theta ))\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{j=1}^N\sum _{z_j}logP(y_j,z_j|\theta )P(z_j|y_j,{\theta }^{(i)})\\ 又由于：& \\ 令：& \\ & {\alpha }_k={\alpha }_{z_j},{\alpha }_k的取值取决于z_j的取值\\ & P(y_j,z_j|\theta )=P(z_j|\theta )P(y_j|z_j\theta )={\alpha }_k\varphi (y_j|{\theta }_k)={\alpha }_{z_j}\varphi (y_j|{\theta }_k)\\ & P(z_j|y_j,{\theta }^{(i)})=\frac{P(z_j,y_j|{\theta }^{(i)})}{P(y_j|{\theta }^{(i)})}=\frac{{\alpha }_k\varphi (y_j|{\theta }_k^{(i)})}{P(y_j|{\theta }^{(i)})}=\frac{{\alpha }_k\varphi (y_j|{\theta }_k^{(i)})}{\sum _{k=1}^K{\alpha }_k\varphi (y_j|{\theta }_k^{(i)})}=\frac{{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}{\sum _{k=1}^K{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}\\ 终于有：& \\ & \underset{\theta }{\mathrm{argmax}}\sum _{j=1}^N\sum _{z_j}logP(y_j,z_j|\theta )P(z_j|y_j,{\theta }^{(i)})\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{j=1}^N\sum _{z_j}log{\alpha }_{z_j}\varphi (y_j|{\theta }_k)\frac{{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}{\sum _{k=1}^K{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}\end{array}$$

M步：

由上面的推证，我们得到：
$$\begin{array}{rl}{\theta }^{(i+1)}& =\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^{(i)})\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{j=1}^N\sum _{z_j}log{\alpha }_{z_j}\varphi (y_j|{\theta }_k)\frac{{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}{\sum _{k=1}^K{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{j=1}^N\sum _{z_j}(log{\alpha }_{z_j}+log\varphi (y_j|{\theta }_k))\frac{{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}{\sum _{k=1}^K{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{j=1}^N\sum _{z_j}log{\alpha }_{z_j}\frac{{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}{\sum _{k=1}^K{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}+log\varphi (y_j|{\theta }_k)\frac{{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}{\sum _{k=1}^K{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}\end{array}$$
可以发现，若要取到最大，只需要加号两边分别最大，即是说：
$$\begin{array}{r}{\alpha }_k=\underset{{\alpha }_{z_j}}{\mathrm{argmax}}\sum _{j=1}^N\sum _{z_j}log{\alpha }_{z_j}\frac{{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}{\sum _{k=1}^K{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}\\ \theta =\underset{\theta }{\mathrm{argmax}}\sum _{j=1}^N\sum _{z_j}log\varphi (y_j|{\theta }_k)\frac{{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}{\sum _{k=1}^K{\alpha }_{z_j}\varphi (y_j|{\theta }_k^{(i)})}\end{array}$$
这样一来，分模型的概率值就可以和高斯分布的参数分开来求了。由于求解过程涉及到宇宙起源，篇幅有限，就不再赘述了。

for details please jump: https://blog.csdn.net/qq_37334135/article/details/85493330
总结：

• EM算法，简单来说就是列出Q函数，求使其最大化的参数，并迭代直到达到理想效果。
• EM算法，不仅仅可以用极大似然估计法，还可以用最大后验估计法。
• 除开高斯分布，其他分布也可以用来建模
• sklearn的mixture模块的讲解的蛮好的，贝叶斯估计也有，可以去看看
• https://scikit-learn.org/stable/modules/mixture.html#gmm

参考书籍：
李航：《统计学习方法第二版》

参考链接：
https://www.bilibili.com/video/BV1aE411o7qd?p=68
https://github.com/datawhalechina/team-learning/tree/master/%E6%9C%BA%E5%99%A8%E5%AD%A6%E4%B9%A0%E7%AE%97%E6%B3%95%E5%9F%BA%E7%A1%80