【学习笔记】统计学习方法——EM算法

1、摘要

EM算法是一种迭代算法，1977年由Dempster等人总结提出，用于含有隐变量（hidden variable）的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望（expectation）；M步，求极大（maximization）。所以这一算法称为期望极大算法（expectation maximization algorithm），简称EM算法。

2、预备知识

2.1、联合概率、边缘概率、条件概率

1、联合概率

联合概率指的是包含多个条件且所有条件同时成立的概率，记作P(X=a,Y=b)或P(a,b)。 一定要注意是所有条件同时成立。

2、边缘概率

边缘概率是与联合概率对应的，P(X=a)或P(Y=b)，这类仅与单个随机变量有关的概率称为边缘概率。

1&2、联合概率和边缘概率的关系
$$\begin{gathered} P(X=a)=\sum _{b}P(X=a,Y=b) \\ P(Y=b)=\sum _{a}P(Y=b,X=a) \end{gathered}$$
求和符号表示穷举所有Y（或X）所能取得b（或a）后，所有对应值相加得到的和。

注意：上面X和Y谁写在前面是没有关系的，我只是为了更容易理解，所以改了下位置。

3、条件概率

条件概率表示在条件Y=b成立的情况下，X=a的概率，记作P(X=a|Y=b)或P(a|b)，它具有如下性质：

“在条件Y=b下X的条件分布”也是一种“X的概率分布”，因此穷举X的可取值之后，所有这些值对应的概率之和为1即：
$$\sum _{a}P(X=a|Y=b)=1$$
1&2&3、联合概率、边缘概率和条件概率之间的关系
$$P(X=a|Y=b)=\frac{P(X=a,Y=b)}{P(Y=b)}$$
当然，如果X和Y相互独立的话公式就变成这样了，$P(X=a)=\frac{P(X=a,Y=b)}{P(Y=b)}$

2.2、极大似然估计

在说极大似然估计前，先来看看贝叶斯决策。

贝叶斯决策

首先来看贝叶斯分类，我们都知道经典的贝叶斯公式：
$$P(w|x)=\frac{P(x|w)P(w)}{P(x)}$$
其中：$P(w)$为先验概率，表示每种类别分布的概率；$P(x|w)$为类条件概率，表示在某种类别前提下，某事发生的概率。而$P(w|x)$为后验概率，表示某事发生了，并且它属于某一类别的概率，了这个后验概率，我们就可以对样本进行分类。后验概率越大，说明某事物属于这个类别的可能性越大。

我们来看一个直观的例子：**已知：**在夏季，某公园男性穿凉鞋的概率为1/2，女性穿凉鞋的概率为2/3，并且该公园中男女比例通常为2:1，**问题：**若你在公园中随机遇到一个穿凉鞋的人，请问他的性别为男性或女性的概率分别为多少？

从问题来看，就是上面讲的，某事发生了，它属于某一类别的概率是多少？即后验概率。

设 $w_1=男性$，$w_2=女性$，$x=穿凉鞋$。

由已知可以得：
$$\begin{gathered} 先验概率P(w_1)=\frac{2}{3}，P(w_2)=\frac{1}{3} \\ 类条件概率P(x|w_1)=\frac{1}{2}，P(x|w_2)=\frac{2}{3} \end{gathered}$$
男性和女性穿凉鞋独立（至少这题没特别说明，两者不独立），所以在公园里遇到一个人穿凉鞋的概率为：
$$P(x)=P(w_1)P(x|w_1)+P(w_2)P(x|w_2)=\frac{2}{3}\times \frac{1}{2}+\frac{1}{3}\times \frac{2}{3}=\frac{1}{3}+\frac{2}{9}=\frac{5}{9}$$
有贝叶斯公式算出：
$$\begin{gathered} P(w_1|x)=\frac{P(x|w_1)P(w_1)}{P(x)}=\frac{\frac{1}{2}\times \frac{2}{3}}{\frac{5}{9}}=\frac{3}{5} \\ P(w_2|x)=\frac{P(x|w_2)P(w_2)}{P(x)}=\frac{\frac{2}{3}\times \frac{1}{3}}{\frac{5}{9}}=\frac{2}{5} \end{gathered}$$
所以遇到一个穿凉鞋的人，判断这个人是男性的概率为$\frac{3}{5}$，是女性的概率为$\frac{2}{5}$。当然贝叶斯决策并不是我们要说的重点。

问题引出

​ 但是在实际问题中并不都是这样幸运的，我们能获得的数据可能只有有限数目的样本数据，而先验概率$P(w_i)$和类条件概率(各类的总体分布)$P(x|w_i)$都是未知的。根据仅有的样本数据进行分类时，一种可行的办法是我们需要先对先验概率和类条件概率进行估计，然后再套用贝叶斯分类器。

​ 先验概率的估计较简单，1、每个样本所属的自然状态都是已知的（有监督学习）；2、依靠经验；3、用训练样本中各类出现的频率估计。

​ 类条件概率的估计（非常难），原因包括：概率密度函数包含了一个随机变量的全部信息；样本数据可能不多；特征向量x的维度可能很大等等。总之要直接估计类条件概率的密度函数很难。解决的办法就是，把估计完全未知的概率密度$P(x|w_i)$转化为估计参数。这里就将概率密度估计问题转化为参数估计问题，极大似然估计就是一种参数估计方法。当然了，概率密度函数的选取很重要，模型正确，在样本区域无穷时，我们会得到较准确的估计值，如果模型都错了，那估计半天的参数，肯定也没啥意义了。

极大似然估计

极大似然估计的原理，用一张图片来说明，如下图所示：

总结起来，最大似然估计的目的就是： 利用已知的样本结果，反推最有可能（最大概率）导致这样结果的参数值。

原理：极大似然估计是建立在极大似然原理的基础上的一个统计方法，是概率论在统计学中的应用。极大似然估计提供了一种给定观察数据来评估模型参数的方法，即：“模型已定，参数未知”。通过若干次试验，观察其结果，利用试验结果得到某个参数值能够使样本出现的概率为最大，则称为极大似然估计。

由于样本集中的样本都是独立同分布，可以只考虑一类样本集D，来估计参数向量θ。记已知的样本集为：
D = { x 1 , x 2 , x 3 , ⋯   , x N } D=\lbrace x_1,x_2,x_3,\cdots,x_N \rbrace D={x1​,x2​,x3​,⋯,xN​}
似然函数（likehood function）：联合概率密度$P(D|\theta )$称为相对于 { x 1 , x 2 , x 3 , ⋯   , x N } \lbrace x_1,x_2,x_3,\cdots,x_N \rbrace {x1​,x2​,x3​,⋯,xN​}的$\theta$的似然幻术。
$$记l(\theta )=P(D|\theta )=p(x_1,x_2,x_3,\cdots ,x_N|\theta )=\prod _{i=1}^{N}P(x_i|\theta )$$
如果$\hat{\theta }$是参数空间中能使似然函数$l(\theta )$最大的$\theta$值，则$\hat{\theta }$应该是最有可能的参数值，那么$\hat{\theta }$就是$\theta$的极大似然估计量。它是样本集的函数，记作：
$$\begin{gathered} \hat{\theta }=d(x_1,x_2,\cdots ,x_N)=d(D) \\ \hat{\theta }(x_1,x_2,\cdots ,x_N)称作极大似然函数估计值 \end{gathered}$$
ML估计：求使得出现改组样本的概率最大的$\theta$值
$$\hat{\theta }=arg{\max }_{\theta }l(\theta )=arg{\max }_{\theta }\prod _{i=1}^{N}P(x_i|\theta )$$
为了便于计算，通常把$\prod$，利用对数可以转换成$\log$。因此就定义了对数似然函数（通常取e作为底数）：
$$\begin{gathered} 记H(\theta )=logl(\theta ) \\ \hat{\theta }=arg{\max }_{\theta }H(\theta )=arg{\max }_{\theta }logl(\theta )=arg{\max }_{\theta }\sum _{i=1}^{N}log(P(x_i|\theta )) \end{gathered}$$
然后通过，求导使式子等于0就可以求出$\theta$的值了。

2.3、Jensen不等式

注：若函数 $f$ 是凹函数，Jensen不等式符号相反。

3、EM算法

3.1、EM算法的引入

首先通过书上的三硬币模型来导入，这里我就直接将书上的题目拍下来了。

$$\{\begin{array}{l}\pi p+(1-\pi )q，y=1\\ \pi (1-p)+(1-\pi )(1-q)，y=0\end{array}\text{\hspace{1em}(1)}$$
将(1)式合并一下就得到：

$$\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}\text{\hspace{1em}(2)}$$
接下来来推导书上的三硬币模型也就是公式(9.1)。

$$P(y|\theta )=\sum _{z}P(y,z|\theta )$$
这一步就是前面说过的联合概率和边缘概率的关系。下面来证明
$$\begin{gathered} \sum _{z}P(y,z|\theta )=\sum _{z}P(z|\theta )P(y|z,\theta ) \\ 左边=\sum _z\frac{P(y,z,\theta )}{P(\theta )} \\ 右边=\sum _z\frac{P(z,\theta )}{P(\theta )}P(y|z,\theta )=\sum _z\frac{P(y,z,\theta )}{P(\theta )} \end{gathered}$$
因此该等式成立。然后再结合公式(2)就可以得出书上的公式(9.1)。

将观测数据表示为$Y=(Y_1,Y_2,Y_3,\cdots ,Y_n{)}^T$，未观测数据表示为$Z=(Z_1,Z_2,Z_3,\cdots ,Z_n{)}^T$，所以观测数据的似然函数为
$$P(Y|\theta )=\sum _{Z}P(Y,Z|\theta )=\sum _{Z}P(Z|\theta )P(Y|Z,\theta ))\text{\hspace{1em}(9.2)}$$
即
$$P(Y|\theta )=\prod _{j=1}^n[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}]\text{\hspace{1em}(9.3)}$$
考虑求模型参数$\theta =(\theta ,p,q)$的极大似然估计，即
$$\hat{\theta }=arg{\max }_{\theta }logP(Y|\theta )\text{\hspace{1em}(9.4)}$$
由于这个问题没有解析解，只有通过迭代的方法求解，EM算法就是可以用于求解这个问题的一种迭代算法。每次迭代包含两步：E步，求期望；M步，求极大化。下面来介绍EM算法。

3.2、EM算法的导出

我们面对一个含有隐变量的概率模型，目标是极大化观测数据（不完全数据）Y关于参数$\theta$的对数似然函数，即极大化
$$L(\theta )=logP(Y|\theta )=log\sum _{Z}P(Y,Z|\theta )=log(\sum _{Z}P(Y|Z,\theta )P(Z|\theta ))\text{\hspace{1em}(9.12)}$$
然后我们要来极大化$L(\theta )$，也就是要对$\theta$求导咯。我们先来展开看看$L(\theta )$长什么样子。
$$L(\theta )=log[(Y|Z_1,\theta )P(Z_1|\theta )+(Y|Z_2,\theta )P(Z_2|\theta )+\cdots +(Y|Z_n,\theta )P(Z_n|\theta )]$$
这个式子求导后的形式自己脑补一下就知道很复杂了吧。所以很难通过求导来求解得到$\theta$。因此我们想一下可不可以把$\sum$符号从$log$当中提出来呢。因此，这里我们前面说过的Jensen不等式就派上用场了。得到了它的下界。
$$L(\theta )=log(\sum _{Z}P(Y|Z,\theta )P(Z|\theta ))=log(\sum _{Z}Q(Z)\frac{P(Y|Z,\theta )P(Z|\theta )}{Q(Z)})\text{\hspace{1em}(3)}$$

$$>=\sum _{Z}Q(Z)log\frac{P(Y|Z,\theta )P(Z|\theta )}{Q(Z)}\text{\hspace{1em}(4)}$$

上面公式(3)中引入了一个新的未知的分布Q(Z)，满足：
$$\sum _{z}Q(z)=1,其中，0<=Q(z)<=1$$
公式(4)用到了Jensen不等式（对数函数是凹函数）：
$$log(E(x))>=E(log(x))$$
其中：
$$\begin{gathered} E(x)=\sum _i{\lambda }_i{x}_i,{\lambda }_i>=0,\sum _i{\lambda }_i=1 \\ {x}_i=\frac{P(Y,Z^i|\theta )}{Q(Z^i)} \\ {\lambda }_i=Q(Z^i) \end{gathered}$$
由 Jensen 不等式可知，等式成立的条件是随机变量是常数，则有：
$$\frac{P(Y,Z^i|\theta )}{Q(Z^i)}=c$$
其中c为常数，对于任意i，我们得到
$$P(Y,Z^i|\theta )=cQ(Z^i)\text{\hspace{1em}(5)}$$
方程(5)两边同时累加和：
$$\sum _{z}P(Y,Z|\theta )=c\sum _{z}Q(z)$$
由于${\sum }_{Z}Q(Z)=1$。从上面式子中得到${\sum }_{Z}P(Y,Z|\theta )=c$

所以：$Q_{(}{Z}^i)=\frac{P(Y,Z^i|\theta )}{c}=\frac{P(Y,Z^i|\theta )}{{\sum }_{Z}P(Y,Z|\theta )}=\frac{P(Y,Z|\theta )}{P(Y|\theta )}=P(Z|Y,\theta )$。这就是书中公式(9.13)上面为什么分子分母同时乘以$P(Z|Y,{\theta }^{(i)})$。

4、EM算法案例-两硬币模型

假设有两枚硬币A、B，以相同的概率随机选择一个硬币，进行如下的掷硬币实验：共做 5 次实验，每次实验独立的掷十次，结果如图中 a 所示，例如某次实验产生了H、T、T、T、H、H、T、H、T、H (H代表正面朝上)。a 是在知道每次选择的是A还是B的情况下进行，b是在不知道选择的是A还是B的情况下进行，问如何估计两个硬币正面出现的概率？

CASE a

已知每个实验选择的是硬币A 还是硬币 B，重点是如何计算输出的概率分布，这其实也是极大似然求导所得。

上面这个式子求导之后发现，5 次实验中A正面向上的次数再除以总次数作为即为 ${\hat{\theta }}_A$，5次实验中B正面向上的次数再除以总次数作为即为 ，即:

CASE b

由于并不知道选择的是硬币 A 还是硬币 B，因此采用EM算法。

计算出每个实验为硬币 A 和硬币 B 的概率，然后进行加权求和。

针对L函数求导来对参数求导，例如对 ${\theta }_A$求导：

求导等于 0 之后就可得到图中的第一次迭代之后的参数值:

当然，基于Case a 我们也可以用一种更简单的方法求得：

class EM:
    def __init__(self, prob):
        self.pro_A, self.pro_B, self.pro_C = prob

    # e_step
    def pmf(self, i):
        # 硬币B产生y的概率
        pro_1 = self.pro_A * math.pow(self.pro_B, data[i]) * math.pow(
            (1 - self.pro_B), 1 - data[i])
        # 硬币A产生y的概率
        pro_2 = (1 - self.pro_A) * math.pow(self.pro_C, data[i]) * math.pow(
            (1 - self.pro_C), 1 - data[i])
        # y来自硬币B的概率，也就是公式(9.5)
        return pro_1 / (pro_1 + pro_2)

    # m_step
    def fit(self, data):
        count = len(data)
        print('init prob:{}, {}, {}'.format(self.pro_A, self.pro_B,
                                            self.pro_C))
        for d in range(count):
            _ = yield
            _pmf = [self.pmf(k) for k in range(count)]
            # 公式(9.6)
            pro_A = 1 / count * sum(_pmf)
            # 公式(9.7)
            pro_B = sum([_pmf[k] * data[k] for k in range(count)]) / sum(
                [_pmf[k] for k in range(count)])
            # 公式(9.8)
            pro_C = sum([(1 - _pmf[k]) * data[k]
                         for k in range(count)]) / sum([(1 - _pmf[k])
                                                        for k in range(count)])
            print('{}/{}  pro_a:{:.3f}, pro_b:{:.3f}, pro_c:{:.3f}'.format(
                d + 1, count, pro_A, pro_B, pro_C))
            self.pro_A = pro_A
            self.pro_B = pro_B
            self.pro_C = pro_C

# 我这里直接使用书上的数据
data = [1,1,0,1,0,0,1,0,1,1]
# 初始化的概率也选择跟书上一样
em = EM(prob=[0.5, 0.5, 0.5])
f = em.fit(data)
next(f)
# 输出 init prob:0.5, 0.5, 0.5

# 第一次迭代
f.send(1)
# 输出 1/10  pro_a:0.500, pro_b:0.600, pro_c:0.600

# 第二次
f.send(2)
# 输出 2/10  pro_a:0.500, pro_b:0.600, pro_c:0.600

可以看出跟书上得出的结果是一样一样的。

# 书中另一个初始化的概率为[0.4, 0.6, 0.7]
em = EM(prob=[0.4, 0.6, 0.7])
f2 = em.fit(data)
next(f2)
# 输出 init prob:0.4, 0.6, 0.7

# 第一次迭代
f2.send(1)
# 输出 1/10  pro_a:0.406, pro_b:0.537, pro_c:0.643

# 第二次迭代
f2.send(2)
# 输出 2/10  pro_a:0.406, pro_b:0.537, pro_c:0.643

# 第三次迭代
f2.send(3)
# 输出 3/10  pro_a:0.406, pro_b:0.537, pro_c:0.643

一直迭代下去也就是书上的那个概率。至此EM算法的一些笔记到这里就结束了，如果跟书本有出入的，还是以书本为准。溜了~

Reference

• https://blog.csdn.net/tick_tock97/article/details/79885868
• https://blog.csdn.net/zengxiantao1994/article/details/72787849
• 李航《统计学习方法》第二版
• https://www.zhihu.com/collection/582310071
• https://zhuanlan.zhihu.com/p/35698329
• http://ai.stanford.edu/~chuongdo/papers/em_tutorial.pdf
• https://github.com/wzyonggege/statistical-learning-method