EM算法

一、EM算法函数

1.1 EM算法问题导入

调查我们学校的男生和女生的身高分布。 假设你在校园里随便找了100个男生和100个女生。调查学生的身高分布。将他们按照性别划分为两组，然后先统计抽样得到的100个男生的身高。假设他们的身高是服从正态分布的。但是这个分布的均值μ和方差 ${\sigma }^2$我们不知道，这两个参数就是我们要估计的。
记作 $\theta =[\mu ,\sigma {]}^T。$

问题分析：设样本集$X=x_1,x_2,\dots ,x_N,$(N=100), $p(x_i|\theta )$为概率密度函数，表示抽到男生$x_i$的概率。

由于100个样本之间独立同分布，所以我同时抽到这100个男生的概率就是他们各自概率的乘积，也就是样本集X中各个样本的联合概率，用 似然函数 $\mathbf{L}(\theta )$ 表示：
$$L(\theta )=L(x_1,x_2,\dots ,x_n;\theta )=\prod _{i=1}^{n}p(x_i;\theta )\theta \in \Theta$$

找到一个参数 $\theta$，使得抽到 $X$ 这组样本的概率最大，也就是说需要其对应的似然函数 $L(\theta )$最大。满足条件的$\theta$ 叫做 $\theta$的最大似然估计量:
$$\hat{\theta }=arg\max L(\theta )$$

1.2 求解最大似然函数 估计值

1. 对 似然函数$L(\theta )$ 取对数：
   $$H(\theta )=\ln (\theta )=\ln \prod _{i=1}^{n}p(x_i;\theta )=\sum _{i=1}^n\ln p(x_i;\theta )$$
2. 对上式求导，令导数为 0 ，得到似然方程。
3. 求解似然方程，得到的参数 $\theta$ 即为所求。

附 符合高斯分布 概率密度函数

若给定一组样本$x_1,x_2\dots x_n，$已知它们来自于高斯分布$N(\mu ,\sigma )$，试估计参数 $\mu ,\sigma$。
$$f(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}（高斯分布的概率密度函数）$$

此时我们已经获得符合高斯分布的概率密度函数：
对 $\mu$ 最大化，我们得到最大似然解（即样本均值）：
$$\mu =\frac{1}{N}\sum _{n=1}^N{x}_n$$

对 ${\sigma }^2$ 最大化，我们得到（样本方差。）：
$${\sigma }^2=\frac{1}{N}\sum _{n=1}^N(x_n-\mu {)}^2$$

但是这个解不是无偏的，由 $\mathbb{E}[x^2]=\mathbb{E}[x{]}^2+\text{var}[x]={\mu }^2+{\sigma }^2$ 我们可以计算它们的期望：
$$\begin{array}{rl}\mathbb{E}[\mu ]& =\frac{1}{N}\sum _{n=1}^N\mathbb{E}[x_n]=\mu \\ \mathbb{E}[{\sigma }^2]& =\mathbb{E}\left[\frac{1}{N}\sum _{n=1}^N(x_n-\frac{1}{N}\sum _{m=1}^N{x}_m)\right]\\ & =\frac{1}{N}\sum _{i=1}^N\mathbb{E}\left[x_n^2-\frac{2}{N}{x}_n\sum _{m=1}{x}_m+\frac{1}{N^2}\sum _{m=1}^N\sum _{l=1}^N{x}_m{x}_l\right]\\ & =({\mu }^2+{\sigma }^2)-2({\mu }^2+\frac{1}{N}{\sigma }^2)+{\mu }^2+\frac{1}{N}{\sigma }^2\\ & =\left(\frac{N-1}{N}\right){\sigma }^2\end{array}$$

用到了：

• 当 $m=n$ 时，$\mathbb{E}[x_m{x}_n]=\mathbb{E}[x_n^2]={\mu }^2+{\sigma }^2$
• 当 $m\ne n$ 时，$\mathbb{E}[x_m{x}_n]=\mathbb{E}[x_n]\mathbb{E}[x_m]={\mu }^2$

因此，方差的一个无偏估计为：
$${\tilde{\sigma }}^2=\frac{N}{N-1}{\sigma }^2=\frac{1}{N-1}\sum _{n=1}^N(x_n-\mu {)}^2$$

随着 $N$ 的增大，方差估计的误差也随之增大。

附 混合高斯模型

直观理解猜测GMM的参数估计：

随机变量X是有K个高斯分布混合而成，取各个高斯分布的概率为π1π2… πK，
第i个高斯分布的均值为μi，方差为Σi。
若观测到随机变量X的一系列样本x1,x2,…,xn，试估计参数π，μ，Σ。

在学习混合高斯模型的推论时，使我想起蒙特卡洛方法。虽然两者相差很大，但有一点相通。
混合高斯模型：通过数据点的分布来反推公式中的参数
蒙特卡洛方法：通过落入特定形状数据点的概率来反推落入特定形状面积。
关于混合高斯模型（GMM）的参数估计 请点击

二、Jensen不等式

lazy Statistician规则：

设 y 是随机变量 x 的函数，$y=g(x)(g是连续函数)$，那么：

(1) x 是离散随机变量，分布为$p(x=x_k)=p_k,k\in N^{\ast },$若${\sum }_{k=1}^{\infty }g(x)p_k$绝对收敛，则：
$$E(y)=E|g(x)|=\sum _{k=1}^{\infty }g(x_k)p_k$$
(2) x 是连续随机变量，分布概率为$f(x)$,若${\int }_{-\infty }^{\infty }g(x)f(x)dx$绝对收敛，则：
$$E(y)=E|g(x)|={\int }_{-\infty }^{\infty }g(x)f(x)dx$$

三、算法流程

3.1 数据及参数

观测数据：观测到的随机变量X的样本：$X=(x_1,...,x_n)$
隐含变量：未观测到的随机变量Z的值：$Z=(z_1,...,z_n)$
完整数据：包含观测到的随机变量 $X$ 和隐含变量 $Z$ 的数据：

$Y=(X,Z)$
$Y=((x_1,z_1),...,(x_n,z_n))$

3.2 EM算法的推导

EM算法是从含有隐含变量的数据(完整数据)中计算极大似然估计。Z为隐含变量，则从可观测数据入手，对参数进行极大似然估计。

• 根据边缘分布列的定义：
  $$\sum _{j=1}^{+\infty }P(X=x_i,Y=y_j)=P(X=x_i)$$

  改写成 $L(\theta )$ 并取对数似然函数：$$\begin{array}{rl}\ell (\theta )& =\sum _{i=1}\ln p(x^{(i)};\theta )\\ & =\sum _i\ln \sum _{z^{(i)}}p(x^{(i)},z^{(i)};\theta )\end{array}$$

  将$p(x^{(i)})$ 用边缘分布列反向拆分为联合分布。

  这里，$z$ 是隐随机变量，直接找到参数的估计是很困难的。我们建立 $l(\theta )$的下界，并求该下界的最大值，重复这个过程，直到收敛到局部最大值。

• 定义隐含变量 $Z$ 的分布 $Q_i$：
  令 $Q_i$ 是z的某一个分布，满足${\sum }_{\mathbf{z}}{\mathbf{Q}}_{\mathbf{i}}{\mathbf{z}}^{(\mathbf{i})}=1$, ${\mathbf{Q}}_{\mathbf{i}}⩾0$，有：
  $$\begin{array}{rl}\sum _i\ln p(x^{(i)};\theta )& =\sum _i\ln \sum _{z^{(i)}}p(x^{(i)},z^{(i)};\theta )\\ & =\sum _i\ln \sum _{z^{(i)}}{Q}_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}\\ & ⩾\sum _i\sum _{z^{(i)}}{Q}_i(z^{(i)})\ln \frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}\end{array}$$
  上述不等式是根据Jensen不等式的来，其中$\ln x$为凹函数。

  为了使上式等号成立：
  $$\begin{array}{rl}\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}& =C\\ & \\ p(x^{(i)},z^{(i)};\theta )& =C(Q_i(z^{(i)}))\\ & \\ \sum _{z}p(x^{(i)},z^{(i)};\theta )& =C(\sum _z{Q}_i(z^{(i)}))\end{array}$$
  满足${\sum }_{\mathbf{z}}{\mathbf{Q}}_{\mathbf{i}}{\mathbf{z}}^{(\mathbf{i})}=1$, ${\mathbf{Q}}_{\mathbf{i}}⩾0$，由上式推导出：
  $$\sum _{z}p(x^{(i)},z^{(i)};\theta )=C$$
  进一步推导：
  $$\begin{array}{rl}{Q}_i(z^{(i)})& =\frac{p(x^{(i)},z^{(i)};\theta )}{{\sum }_{z}p(x^{(i)},z^{(i)};\theta )}\\ & \\ & =\frac{p(x^{(i)},z^{(i)};\theta )}{p(x^{(i)};\theta )}\\ & \\ & =p(z^{(i)}|x^{(i)};\theta )\end{array}$$
  最终我们得到了EM算法的一般形式，一般化形式的思想是：
  在 E-step：找到对于当前参数 $\theta$ 使不等式成立的 $Q$分布。
  在 M-step：对似然函数下界进行极大似然估计，得到新的参数，形式化表述为：

$$\begin{array}{rl}\text{E-step}& Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\theta )\\ \text{M-step}& \theta =arg\underset{\theta }{\max }\sum _{i=1}^n\sum _{z^{(i)}}{Q}_i(z^{(i)})\ln \frac{p(x^{(i)}),z^{(i)};\theta }{Q_i(z^{(i)})}\end{array}$$

3.1 算法思路

直线式迭代优化的路径：

可以看到每一步都会向最优值前进一步，而且前进路线是平行于坐标轴的，因为每一步只优化一个变量。
这犹如在x-y坐标系中找一个曲线的极值，然而曲线函数不能直接求导，因此什么梯度下降方法就不适用了。

但固定一个变量后，另外一个可以通过求导得到，因此可以使用坐标上升法，一次固定一个变量，
对另外的求极值，最后逐步逼近极值。
对应到EM上， E步：固定θ，优化Q；
M步：固定Q，优化θ；
交替将极值推向最大。 （如图）

3.2 EM 算法示例：

● 假设现在有两枚硬币1和2,随机抛掷后正面朝上概率分别为$P_1，P_2。$
为了估计这两个概率，做实验，每次取一枚硬币，连掷5下，记录下结果，
如下(表一（实验数据）)：

表一（实验数据）				表二（测试数据）
硬币	结果	统计		硬币	结果	统计
1	正正反正反	3正-2反		UnKnow	正正反正反	3正-2反
2	反反正正反	2正-3反		UnKnow	反反正正反	2正-3反
1	正反反反反	1正-4反		UnKnow	正反反反反	1正-4反
2	正反反正正	3正-2反		UnKnow	正反反正正	3正-2反
1	反正正反反	2正-3反		UnKnow	反正正反反	2正-3反

表一：可以很容易地估计出P1和P2，如下：
P1 = （3+1+2）/ 15 = 0.4
P2= （2+3）/10 = 0.5
● 现在我们抹去每轮投掷时使用的硬币标记，如上右表（表二（测试数据））。

目标没变，还是估计P1和P2:

此时多了一个隐变量z，可以把它认为是一个5维的向量（z1,z2,z3,z4,z5)，
代表每次投掷时所使用的硬币，比如z1，就代表第一轮投掷时使用的硬币是1还是2。
但是，这个变量z不知道，就无法去估计P1和P2。

我们可以先随机初始化一个P1和P2，用它来估计z，然后基于z，还是按照最大似然概率法则去估计新的P1和P2
当与我们初始化的P1和P2一样时，说明是P1和P2很有可能就是真实的值。这里面包含了两个交互的最大似然估计。

概率计算

当假设P1 = 0.2，P2 = 0.8时： （看看第一轮抛掷最可能是哪个硬币。）

如果是硬币1，得出3正2反的概率为：
0.2×0.2×0.2×0.8×0.8= 0.00512
如果是硬币2，得出3正2反的概率为：
0.8×0.8×0.8×0.2×0.2=0.03087
依次求出其他4轮中的相应概率。（如：表三）

按照最大似然法则：（表三）	轮数	若是硬币1	若是硬币2
第1轮中最有可能的是硬币2	1	0.00512	0.03087
第2轮中最有可能的是硬币1	2	0.02048	0.01323
第3轮中最有可能的是硬币1	3	0.08192	0.00567
第4轮中最有可能的是硬币2	4	0.00512	0.02048
第5轮中最有可能的是硬币1	5	0.02048	0.01323

我们就把上面的值作为z的估计值。

然后按照最大似然概率法则来估计新的P1和P2。（由表三，查到 表二的正反数据计算P值）
P1 = （2+1+2）/15 = 0.33
P2=（3+3）/10 = 0.6
将算的的P1P2反复迭代：
设想我们知道每轮抛掷时的硬币就是如标示的那样，那么，P1和P2的最大似然估计就是0.4和0.5
（下文中将这两个值称为P1和P2的真实值）。那么对比下我们初始化的P1和P2和新估计出的P1和P2：

初始化P1	估计出P1	真实的P1		初始化P2	估计出P2	真实的P2
0.2	0.33	0.4		0.7	0.6	0.5

3.3 示例优化：

新估计出的 P1和P2 一定会更接近真实的 P1和P2
迭代不一定会收敛到真实的P1和P2。取决于P1和P2的初始化值。

我们是用最大似然概率法则估计出的z值，然后再用z值按照最大似然概率法则估计新的P1和P2。
也就是说，我们使用了一个最可能的z值，而不是所有可能的z值。
所有的z值有多少个呢？显然，有2^5=32种，需要我们进行32次估值？？？

用期望来简化运算

再次利用表三：我们可以算出每轮抛掷中使用硬币1或者使用硬币2的概率。

比如第1轮，使用硬币1的概率是：0.00512/(0.00512+0.03087)=0.14
使用硬币2的概率是：1-0.14=0.86

我们按照期望最大似然概率的法则来估计新的P1和P2：
以P1估计为例，（表二）第1轮的3正2反相当于
0.14*3=0.42正
0.14*2=0.28反
依次算出其他四轮，表四如下：new_P1=4.22/(4.22+7.98)=0.35

轮数	正面	反面
1	0.42	0.28
2	1.22	1.83
3	0.94	3.76
4	0.42	0.28
5	1.22	1.83
总计	4.22	7.98

可以看到，改变了z值的估计方法后，新估计出的P1要更加接近0.4。原因就是我们使用了所有抛掷的数据，而不是之前只使用了部分的数据。我们根据E步中求出的z的概率分布，依据最大似然概率法则去估计P1和P2，被称作M步。

总结