EM算法

文章目录

一、EM算法函数

1.1 EM算法问题导入

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

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

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

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

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

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

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

若给定一组样本$x_1,x_2…x_n，$已知它们来自于高斯分布$N(μ,σ)$，试估计参数 $μ,σ$。

$f(x) = \frac{1}{\sqrt{2\pi }\sigma }e^{-\frac{(x-\mu )^2}{2\sigma^2 }}~~~~~~~（高斯分布的概率密度函数）$

此时我们已经获得符合高斯分布的概率密度函数：

对 $μ$ 最大化，我们得到最大似然解（即样本均值）：

$\mu = \frac 1 N \sum_{n=1}^N x_n$

对 $σ^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{aligned} \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_mx_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{aligned}$

用到了：

• 当$m=n$时，$\mathbb E[x_mx_n]=\mathbb E[x_n^2] = \mu^2 + \sigma^2$
• 当$m\neq n$时，$\mathbb E[x_mx_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$

附 混合高斯模型

随机变量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^*,$若$\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{aligned} \ell (\theta)&=\sum_{i=1} \ln p(x^{(i)};\theta) \\ &=\sum_{i} \ln \sum_{z^{(i)}} p(x^{(i)},z^{(i)};\theta) \end{aligned}$
  将$p(x^{(i)})$
  这里，$z$ 是隐随机变量，直接找到参数的估计是很困难的。我们建立 $l(\theta)$的下界，并求该下界的最大值，重复这个过程，直到收敛到局部最大值。
• 定义隐含变量$Z$ 的分布 $Q_i$：
  令 $Q_i$ 是z的某一个分布，满足$\mathbf{\sum_zQ_iz^{(i)}=1}$, $\mathbf{Q_i \geqslant 0}$，有：
  $\begin{aligned} \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)})} \\ &\geqslant \sum_{i}\sum_{z^{(i)}} Q_i(z^{(i)}) \ln \frac{p(x^{(i)},z^{(i)};\theta)}{ Q_i(z^{(i)})} \end{aligned}$
  上述不等式是根据Jensen不等式的来，其中$\ln x$为凹函数。
  为了使上式等号成立：
  $\begin{aligned} \frac{p(x^{(i)},z^{(i)};\theta)}{ Q_i(z^{(i)})} &= C \\ &\\ p(x^{(i)},z^{(i)};\theta)&=C(Q_i(z^{(i)})) \\ &\\ \sum_zp(x^{(i)},z^{(i)};\theta)&=C(\sum_zQ_i(z^{(i)})) \end{aligned}$
  满足$\mathbf{\sum_zQ_iz^{(i)}=1}$, $\mathbf{Q_i \geqslant 0}$，由上式推导出：
  $\sum_zp(x^{(i)},z^{(i)};\theta)=C$
  进一步推导：
  $\begin{aligned} Q_i(z^{(i)}) &= \frac{p(x^{(i)},z^{(i)};\theta)}{\sum_zp(x^{(i)},z^{(i)};\theta)} \\ &\\ &=\frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)} \\ &\\ &=p(z^{(i)}|x^{(i)};\theta) \end{aligned}$
  最终我们得到了EM算法的一般形式，一般化形式的思想是：
  在 E-step：找到对于当前参数 $\theta$ 使不等式成立的 $Q$分布。
  在 M-step：对似然函数下界进行极大似然估计，得到新的参数，形式化表述为：

$\begin{aligned} \textbf{E-step} ~~~~~~~& Q_i(z^{(i)}) = p(z^{(i)}|x^{(i)};\theta)\\ \textbf{M-step} ~~~~~~~& \theta =arg \max_{\theta} \sum_{i=1}^n\sum_{z^{(i)}} Q_i(z^{(i)}) \ln \frac{p(x^{(i)}),z^{(i)};\theta}{Q_i(z^{(i)})} \end{aligned}$

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
● 现在我们抹去每轮投掷时使用的硬币标记，如上右表（表二（测试数据））。

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

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

概率计算

如果是硬币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

然后按照最大似然概率法则来估计新的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轮，使用硬币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步。

总结