EM算法学习笔记

引入

设袋子里有红白两种球，比例为 $p:1-p$ 。现在不知道比例值是多少，因此可以随机抽样： $x_1,x_2...x_n$ ，每一个 $x_i$ 都是随机变量的随机抽样结果，并进行如下估计： $p=\frac{|\{x_i=红\}|}{n}$

如果有两个袋子，各自都有红色的球和白色的球，抽取的时候先随机挑选一个袋子，设选第一个袋子的概率为 $w$ ，第一个袋子中 $红球:白球=p:1-p$ ，第二个为 $红球:白球=q:1-q$ 。

如果此时进行抽取，抽取到红球的概率为： $P\{x=红\}=w\times p+(1-w)\times q$ ，抽到白球的概率同理。

如果想要知道从第一个袋子中抽到红球的概率： $P\{第一个|x_i=红\}$ ，可以把概率乘以权重放到矩阵中：

$\begin{pmatrix} wp&w(1-p)\ (1-w)q&(1-w)(1-q)\ \end{pmatrix}$

取矩阵的第一列分析，不难得出 $P\{第一个|x_i=红\}=\frac{wp}{wp+(1-w)q}$ 。

此时模型的 $w,p,q$ 三个参数均未知，可以使用极大似然法对其进行估计。令1代表红色，0代表白色，每个 $x_i$ 的密度函数为 $wp^{x_i}(1-p)^{1-x_i}+(1-w)q^{x_i}(1-q)^{1-x_i}$ ，如果 $x_i=1$ (红色)，带入上式可以发现结果就是 $wp+(1-w)q$ ，也就是矩阵第一列的和。

根据极大似然估计，需要把所有样本的密度函数做乘积，似然函数就是：

$\Pi_{i=1}^n [wp^{x_i}(1-p)^{1-x_i}+(1-w)q^{x_i}(1-q)^{1-x_i}]$ ，需要对其求极大值。如果按照极大似然的思路，会发现这个问题很复杂（中括号里的部分不易处理，不信你试试）。不过这个函数可以转换为如下的形式：

$\Pi_{i=1}^n [(wp+(1-w)q)^{x_i}(w(1-p)+(1-w)(1-q))^{1-x_i}]$ ，这直观上也很容易理解，因为 $x_i$ 只有0和1两种形式（带入 $x_i=1 \ or\ 0$ 试试）。对其处理的话，可以令 $P=wp+(1-w)q$ ， $Q=w(1-p)+(1-w)(1-q)$ ，函数继续变形为 $\Pi_{i=1}^nP^{x_i}Q^{1-x_i}=\Pi_{i=1}^nP^{x_i}(1-P)^{1-x_i}$ ，对其取对数后得 $(\Sigma_{i=1}^nx_i)lnP+(n-\Sigma_{i=1}^nx_i)ln(1-P)$ ，对P求导并令导数为0，可得 $P=\frac{\Sigma_{i=1}^nx_i}{n}=\frac{k}{n}$ ， $k$ 为样本中1（红球）的个数，且有 $w(1-p)+(1-w)(1-q)=\frac{n-k}{n}$ （本质上等于上面的方程）。

总结一下，刚才设计的模型（先随机选袋子，再随机选球）从极大似然估计的角度， $wp+(1-w)q=\frac{k}{n}$ 就是最佳估计，但是一个方程里含有三个未知数，所以只要三个参数满足这个方程就是一个合理的估计，即给出两个参数就可以把另一个参数解出来。

从另一个角度看问题

现在给出三个参数 $p^{(0)}, q^{(0)}, w^{(0)}$ ，对于 $x_i$ ，要根据它的颜色判断它来自于哪个袋子（利用上面的矩阵来计算）：

红球来自第一个袋子的概率： $\frac{wp}{wp+(1-w)q}=a_i$ ，白球来自第一个袋子的概率： $\frac{w(1-p)}{w(1-p)+(1-w)(1-q)}=b_i$ ...同理设红球和白球来自第二个袋子的概率为 $c_i,d_i$ ，那么对于每个球，这四个概率都可以求出来，列成一张表：

	$x_1$	$x_2$	..
第一个-红	$\frac{wp}{wp+(1-w)q}=a_1$
第一个-白	$\frac{w(1-p)}{wp+(1-w)(1-q)}=b_1$
第二个-红	$\frac{(1-w)p}{wp+(1-w)q}=c_1$
第二个-白	$\frac{(1-w)(1-p)}{wp+(1-w)(1-q)}=d_1$

每一列的总和都是一样的（都是1），因为这里考虑的是后验概率，即每个球已知是红色还是白色，所以对于每一列只是把和已知颜色对应的行加起来（而非把所有的行加起来）。

因此来自于第一个袋子的概率 $w$ 可以给出一个重新的估计： $w^{(1)}=\frac{\Sigma_{x_i=1}a_i+\Sigma_{x_i=1}b_i}{n}$ ，即表中前两行相加再除以n。这个过程就是：选定 $w^{(0)}$ 为初值，根据初值计算出一个个概率值，再根据上表的概率值计算出新的 $w^{(1)}$ 。

同理 $p^{(1)}=\frac{\Sigma_{x_i=1}a_i}{\Sigma_{x_i=1}a_i+\Sigma_{x_i=0}b_i}$ ， $q^{(1)}=\frac{\Sigma_{x_i=1}c_i}{\Sigma_{x_i=1}c_i+\Sigma_{x_i=0}d_i}$ 。如此可以继续迭代下去得到新的参数值，能够证明这个算法是收敛的，至少可以收敛到局部极大值（不一定是全局）。

整个过程就叫做EM算法。

推广

对于高斯分布 $N(\mu,\sigma^2)$ ，设有 $K$ 个高斯分布，从中随机抽取到的概率分别为 $w_1,w_2...w_k$ ，且满足 $\Sigma_{i=1}^k=1$ 。

$X$ 的密度函数 $w_1f_1(x)+....+w_kf_k(x)$ ， $f(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$

现在给出 $x_1,...x_n$ ，如何估计出所有的 $\mu_j,\sigma_j^2,w_i$ ？

按照极大似然的过程， $max\ \Pi_{i=1}^n[\Sigma_{j=1}^Kw_jf_j(x_i)]$ 等价于 $max\ \Sigma_{i=1}^nlog(\Sigma_{j=1}^Kw_jf_j(x_i))$ 。但是这样并没有对问题进行化简，所以根据前面的结果，给出一组参数 $\mu_j^{(0)},\sigma_j^{2(0)},w_i^{(0)}$ ，考虑建一个表：

高斯分布	$x_1$	$x_2$	..
1	$\frac{w_1f_1(x_1)}{w_1f_1(x_1)+...+w_kf_k(x_1)}=w_{11}$	$\frac{w_1f_1(x_2)}{w_1f_1(x_2)+...}=w_{21}$	..
2	$\frac{w_2f_2(x_1)}{w_1f_1(x_1)+...}=w_{12}$	..	..
..	..	..	..
k	..	..	..

这样就又可以进行估计了。同样，每一列的总和是1，则对于 $w_1$ 有一个很自然的估计： $w_1^{(1)}=\frac{\Sigma_{i=1}^nw_{i1}}{n}$ ，其余的w值同理。现在看均值和方差：

$\mu_1^{(1)}=\frac{\Sigma_{i=1}^nw_{i1}x_i}{\Sigma_{i=1}^nw_{i1}}$ ，注意这里分布是对每一行的 $w_{i1}$ 求和，并不是1。 $\sigma_1^{(1)}=\frac{\Sigma_{i=1}^nw_{i1}(x_i-\mu_1^{(1)})^2}{\Sigma_{i=1}^nw_{i1}}$