经典机器学习模型(九)EM算法的推导

1 相关数据基础

1.1 数学期望

1.1.1 数学期望的定义

根据定义，我们可以求得掷骰子对应的期望：
$$\begin{gathered} E(X)=X_1\ast p(X_1)+X_2\ast p(X_2)+...+X_6\ast p(X_6) \\ =1\ast \frac{1}{6}+2\ast \frac{1}{6}+1\ast \frac{1}{6}+3\ast \frac{1}{6}+4\ast \frac{1}{6}+5\ast \frac{1}{6}+6\ast \frac{1}{6} \\ =3.5 \end{gathered}$$

• 要注意区分平均值和期望。平均值是一个统计量(对观察样本的统计)，期望是一种概率论概念，是一个数学特征。比如我们进行掷骰子，掷了六次，点数分别为2，2，2，4，4，4，这六次的观察就是我们的样本，于是我们可以说平均值为(2+2+2+4+4+4)/6=3，但是千万不能说期望是3。
• 平均值和期望的联系也是大数定理联系起来的。如果说概率是频率随样本趋于无穷的极限 ，期望就是平均数随样本趋于无穷的极限。
• 可以用加权平均值来理解期望。

1.1.2 函数期望公式

1.1.3 常见分布的期望

具体推导可以参考：

数学期望及常见分布的期望计算与推导

1.2 极大似然估计

• 极大似然估计法的出发点是已知被观测对象的分布，但不知道其参数。
• 极大似然法用得到观测值（样本）最高概率的那些参数的值来估计该分布的参数，即产生该样本概率最大的原则。
• 设$f(y,\theta )$ 是随机变量 $Y$ 的密度函数，其中$\theta$是该分布的未知参数，若有一随机样本 $Y_1,Y_2,\cdots ,Y_n$，则$\theta$的极大似然估计值是具有产生该观测样本的最高概率的那个$\theta$值，或者换句话说，$\theta$的极大似然估计值是使密度函数 $f(y,\theta )$ 达到最大的$\theta$值。
• 由于总体有离散型和连续型两种分，离散型分布通过分布律来构造似然函数，而连续型分布通过概率密度函数来构造似然函数。

求解极大似然估计问题步骤：

1. 写出似然函数；
2. 对似然函数取对数，并整理；
3. 求导数，令导数为0，得到似然方程；
4. 解似然方程，得到的参数即为所求；

1.2.1 离散型随机变量的极大似然原理

• 若总体为离散型分布，其分布律为 $P(Y=y)=p(y,\theta )$，分布律的形式已知。
  • 其中， $\hat{\theta }=({\theta }_1,{\theta }_2,\cdots ,{\theta }_k)\prime$ 是待估参数向量
  • 其中一维随机变量离散型分布主要有：

设 $Y_1,Y_2,\cdots ,Y_n$ 表示总体 $Y$ 的一个样本，它们独立同分布。 $y_1,y_2,\cdots ,y_n$ 是相应于样本 $Y_1,Y_2,\cdots ,Y_n$ 的一组样本值，容易求得从 $Y_1,Y_2,\cdots ,Y_n$ 取到观察值$ y_1,y_2,⋯,y_n$ 的概率，即事件 {$Y_1=y_1,Y_2=y_2,\cdots ,Y_n=y_n$}发生的概率，这个概率为
$$L(\hat{\theta })=L(y_1,y_2,...,y_n;\hat{\theta })=\prod _{i=1}^{n}p(y_i,\hat{\theta })$$

• 这一概率随$\hat{\theta }$的取值而变化，它是$\hat{\theta }$的函数，$L(\hat\theta) $称为样本的似然函数。

• 极大似然估计法就是在$\hat{\theta }$取值的可能范围内挑选使似然函数 $L(y_1,y_2,\cdots ,y_n;\hat{\theta })$ 达到最大的参数值$\hat{\theta }$

• 一般求解的步骤就是：对该似然函数取对数，然后求导数，令导数为0，得到似然方程，最后解似然方程。取对数是因为对数似然函$lnL(\hat{\theta })$ 可以把乘积形式转为和的形式，从而为简化运算提供了方便。

1.2.2 离散型分布的极大似然估计的举例

首先来看1次抛硬币，假设参数正面向上的概率为$\theta$，满足伯努利分布(也称0-1分布)，可能的事件有2个（正面向上的次数可能为0、1次），其正面的条件概率为：
$$p(x)={\theta }^x(1-\theta {)}^{1-x},x为正面向上的次数$$
现在我们朝空中扔$N$次，其中有$x$次显示的是正面，有$N-x$次显示的是反面，那么它所对应的**「正面的条件概率」**就可以写成下式（满足二项分布）：
$$P(x|\theta )=C_N^x{\theta }^x(1-\theta {)}^{N-x},x为正面向上的次数$$

现在，我们已经知道硬币正面向上的概率$\theta$、一共抛了$N$次，$x$次显示的是正面，那么代入上式，很容易就能求出该次事件出现的概率。

假如，现在知道硬币一共抛了$N$次，$x$次显示的是正面，但是不知道硬币朝上的概率$\theta$，那么该怎么办呢？这时，咱们就可以让条件概率最大化来找到对应的$\theta$，即
$$arg\underset{\theta }{\max }P(x|\theta )$$
此时，我们可以把它写成似然函数的形式$L(\theta )$，当然，由于原来的条件概率函数都是指数乘积的形式，为了计算方便，我们接着把似然函数写成 【对数似然函数】。
$$\begin{gathered} \hat{\theta }=arg\underset{\theta }{\max }ln(L(\theta )) \\ =arg\underset{\theta }{\max }ln({\theta }^x(1-\theta {)}^{N-x})(忽略常数项) \\ =arg\underset{\theta }{\max }(xln\theta +(N-x)ln(1-\theta )) \\ 求最大值，我们对\theta 求导，并令其为0，那么 \\ \frac{\partial L(\theta )}{\partial \theta }=\frac{\partial (xln\theta +(N-x)ln(1-\theta ))}{\partial \theta }=\frac{x}{\theta }-\frac{N-x}{1-\theta }=0 \\ 可以求得\theta =\frac{x}{N} \end{gathered}$$

• 频率学派相信概率是确定的，或者说，$\theta$是个常量，采样数据$x$则是基于这个参数为$\theta$的分布中随机采样的，因此通过采样数据可以求得$\theta$（通过极大似然估计MLE），而采样数据越多，$\theta$越准确。

• 而贝叶斯学派则认为θ并非是个未知的常量，而是个满足某种分布的随机变量，而对于这个θ，会有一个最初始的信仰，即一个先验假设（比如抛硬币中，θ可以被视为一个均值为0.5的正态分布）。

• 具体可参考：概率学派和贝叶斯学派的区别

• 可以看出，极大似然估计可以看成是对应于一组完全数据的情况，但是当出现不完全的数据时，比如未被观测到或者是缺失的数据时，这时用极大似然估计来求解就相当复杂了。

1.2.3 连续型随机变量的极大似然原理

若总体为连续型分布，其概率密度函数为$f(y,\hat{\theta })$，密度函数的形式已知。

• 其中， $\hat{\theta }=({\theta }_1,{\theta }_2,\cdots ,{\theta }_k)\prime$ 是待估参数向量
• 其中一维随机变量连续型分布主要有：

解法和离散型分布一致，对该似然函数取对数，然后求导数，令导数为0，得到似然方程，最后解似然方程。

1.2.4 连续型分布的极大似然估计举例

我们假设样本服从高斯分布：

1.3 Jensen不等式

• 这里简单介绍下结论，感兴趣的可以详细了解下该不等式。

• 如果$f$是凸函数(如下图)，X是随机变量，那么有$E[f(X)]>=f(E[X])$，也就是函数的期望大于等于期望的函数。

• 对于凹函数，不等号方向反向，即$E[f(X)]<=f(E[X])$。

2 EM算法举例

如下图，我们抛两枚硬币A和B，一共抛了5轮，每轮抛10次。

如果知道每次抛的是A还是B，那么根据之前讲的极大似然估计，就直接可以估计每种硬币的参数${\theta }_A,{\theta }_B$（正面朝上的概率）。

假如此时我们并不知道每轮抛掷的是A硬币还是B硬币，只能知道每组实验的10次结果。这时候，我们就需要EM算法了，这时每组未知的硬币就是隐变量。

• EM算法的核心就是猜数+迭代。

• 对于第一轮抛掷，使用硬币 A 的概率是 0.45，使用硬币 B 的概率是 0.55。同理其他轮。这一步我们实际上是估计出了 Z 的概率分布，这步就是 E-Step。

到隐变量$z$后(即每轮是A硬币，还是B硬币)，我们可以去进行M步计算极大似然估计求得更好的θ

3 EM算法的推导

3.1 EM算法流程

我们先看下《统计学习方法》中EM算法的流程，然后我们再去进行推导。

3.2 EM算法中E步的推导

对$m$个样本观察数据$y=(y^{(1)},y^{(2)},...y^{(m)})$中，找出样本的模型参数θ，最大化模型分布的对数似然函数如下：
$$\begin{gathered} \hat{\theta }=arg\underset{\theta }{\max }log(L(\theta )) \\ =arg\underset{\theta }{\max }\sum _{i=1}^{m}log(P(y^i|\theta )) \\ =arg\underset{\theta }{\max }log(P(Y|\theta )) \\ 它表示参数\theta 的条件下对应的观测变量Y的概率对数化 \\ Y是离散变量时，对应的是概率;Y是连续变量，对应的是概率密度。 \\ 后文以概率为代表解释，概率密度类似。 \end{gathered}$$

1）因为这里包含缺失数据和隐变量，无法直接得到对应的概率，于是把对应的隐变量$Z$添进来，以全概率公式展开。
$$\begin{gathered} L(\theta )=log(P(Y|\theta ))=log\sum _Z(P(Y,Z|\theta )) \\ 这里的隐变量Z是对样本空间的分割，就是在Z的所有取值下对应的概率。 \end{gathered}$$
2）接下来，我们就要用到极大似然估计的思想了。

我们令每次迭代后的$L(\theta )$都比上轮的$L({\theta }^{(i)})$大，目的就是为了让更新之后的似然函数更大，即：
$$\begin{gathered} L(\theta )-L({\theta }^{(i)})=log\sum _Z(P(Y,Z|\theta ))-log(P(Y|{\theta }^{(i)}))\ge 0 \\ 注意：上式中{\theta }^{(i)}已知(上一轮推测出来的)，\theta 是未知的 \end{gathered}$$

• 等号右侧前面的式子是将所有隐变量考虑在内的全概率密度函数，也就是在EM例子中，既然不知道每轮是A硬币还是B硬币，那就把A和B的概率都分别计算出来。
• 而后面的式子中，参数${\theta }^{(i)}$是根据上轮结果推测出来A硬币/B硬币正面朝上的概率，我们认为它是已知的，目标就是用它来推测出新的一轮隐变量对应的概率，并且找到新一轮的似然函数。

3）接着，用贝叶斯公式，来把前面的式子进行改造，通过观测结果$Y$和上一轮的参数${\theta }^{(i)}$去猜出隐变量 对应的概率$Z$，将，上式就可以继续写成：
$$\begin{gathered} L(\theta )-L({\theta }^{(i)})=log\sum _Z(P(Y,Z|\theta ))-log(P(Y|{\theta }^{(i)})) \\ =log\sum _Z[P(Z|Y,{\theta }^{(i)})\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^{(i)}}]-log(P(Y|{\theta }^{(i)})) \\ 分子分母同乘以一个数，原式保持不变，这样就引入了隐变量z的概率分布。 \\ 另外对数函数中有一个求和式，直接计算是十分复杂的。 \end{gathered}$$

4）此时，我们需要借助Jensen不等式了。

我们已经知道，期望$E(X)=\sum xp(x)$，函数期望则为：$E(f(x))=\sum f(x)p(x)$

$\sum _{Z}P(Z|Y,{\theta }^{(i)})=1$，这是因变量$z$的概率分布，相当于上式的$p(x)$。

我们$\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^{(i)})}$把看作一个整体，相当于上式中的$f(x)$，因此$\sum _Z[P(Z|Y,{\theta }^{(i)})\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^{(i)})}]$可以看作一个期望。

根据Jensen不等式，对于凸函数，函数的期望大于等于期望的函数；凹函数相反

$log\sum _Z[P(Z|Y,{\theta }^{(i)})\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^{(i)}}]$就是期望的函数，此函数为log，而log属于凹函数，取相反结论，因此函数的期望小于等于期望的函数。

我们可以得到下式：
$$\begin{gathered} L(\theta )-L({\theta }^{(i)})=log\sum _Z(P(Y,Z|\theta ))-log(P(Y|{\theta }^{(i)})) \\ =log\sum _Z[P(Z|Y,{\theta }^{(i)})\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^{(i)})}]-log(P(Y|{\theta }^{(i)})) \\ \ge \sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-log(P(Y|{\theta }^{(i)})) \\ 由于\sum _{Z}P(Z|Y,{\theta }^{(i)})=1，因此在log(P(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|Y,{\theta }^{(i)})}-\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Y|{\theta }^{(i)})) \\ =\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^{(i)})(P(Y|{\theta }^{(i)}))} \\ 接着，对分母的式子进行简化，根据乘法公式P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})=P(Y,Z|{\theta }^{(i)})，它可以继续写成 \\ =\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^{(i)})} \\ \ge 0 \end{gathered}$$
5)可以看出，对数中的分子是在未知参数$\theta$下的联合概率，而分母是在已知参数${\theta }^{(i)}$的联合概率，都是完全概率，我们的目的是为了让每轮迭代后的概率比上轮的大，自然写成：
$$\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^{(i)})}\ge 0$$
到这一步，我们的似然函数基本上就敲定了。

我们来验证下：令似然函数最大就意味着令每轮迭代后的隐变量概率都比上轮大。

很简单，只要对数部分大于零，不等式自然成立，也就是对数里面的部分大于等于1，继续写成
$$\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^{(i)})}\ge 1$$
验证了每轮迭代的概率都比上轮大，等同于令似然函数越来越大。

6)接着，我们的目标就是要求出对应似然函数下的缺失变量$\theta$的概率值。
$$\begin{gathered} 我们已经知道L(\theta )-L({\theta }^{(i)})\ge \sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^{(i)})} \\ 即L(\theta )\ge L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^{(i)})} \\ 我们令L(\theta )的下界B(\theta ,{\theta }^{(i)})=L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^{(i)})} \\ 我们极大化似然函数L(\theta )等同于最大化下界B(\theta ,{\theta }^{(i)}) \\ 而L({\theta }^{(i)})为常数，就等同于最大化\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^{(i)})} \end{gathered}$$
我们将式子拆开：
$$\begin{gathered} \sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^{(i)})} \\ =\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )-\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|{\theta }^{(i)}) \end{gathered}$$
在减号后面的部分还是个常数，不好含未知参数$\theta$。因此，我们只需要对前面的部分求出当它最大时对应的$\theta$ 。
$$\begin{gathered} 即Q(\theta ,{\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta ) \\ \sum _{Z}P(Z|Y,{\theta }^{(i)})=1，是一个条件概率分布，因此可以化简为期望形式，即： \\ =E_{Z|Y,{\theta }^{(i)}}logP(Y,Z|\theta ) \end{gathered}$$
这就是EM算法的E步的推导。

通过EM算法的E步，我们得到了Q函数
$$\begin{gathered} Q(\theta ,{\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta ) \\ =E_{Z|Y,{\theta }^{(i)}}logP(Y,Z|\theta ) \end{gathered}$$
接下来，我们只需要求使$Q(\theta ,{\theta }^{(i)})$极大化的$\theta$，确定第$i+1$次的迭代参数的估计值${\theta }^{(i+1)}$，即M步。

3.3 EM算法的图形解释

$$\begin{gathered} 在推导EM算法过程中我们知道，L(\theta )\ge L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^{(i)})} \\ 我们令下界B(\theta ,{\theta }^{(i)})=L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^{(i)})} \end{gathered}$$

• 我们想让似然函数越来越大，等同于是让下界$B(\theta ,{\theta }^{(i)})$越来越大。

• 注意，在每轮的迭代中都它所对应的下界函数$B(\theta ,{\theta }^{(i)})$都是在不断更新的。

• 我们说明一下，下图中$L(\theta )$和$B(\theta ,{\theta }^{(i)})$在${\theta }^{(i)}$处相等这句话。

$$\begin{gathered} L(\theta )\ge L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^{(i)})} \\ 当\theta ={\theta }^{(i)}时，\frac{P(Y,Z|{\theta }^{(i)})}{P(Y,Z|{\theta }^{(i)})}=1,那么log\frac{P(Y,Z|{\theta }^{(i)})}{P(Y,Z|{\theta }^{(i)})}=0 \\ 因此L({\theta }^{(i)})\ge L({\theta }^{(i)})+0,即L(\theta )和B(\theta ,{\theta }^{(i)})在{\theta }^{(i)}处相等 \end{gathered}$$

• 现在，我们再看3.1EM算法的流程，就应该有更深刻的理解了。

有关EM算法收敛的证明可以参考：

学会EM算法

接下来会介绍下《统计学习方法》中EM算法的三硬币模型例子，以及EM算法在高斯混合模型中的应用。