机器学习笔记-Task03-极大似然估计&EM算法

1.简介

1.1何为EM算法?

EM算法的全称是Expectation-maximization algorithm（期望最大化算法），是在概率模型中寻找参数最大似然估计或者最大后验估计的算法，其中概率模型依赖于无法观测的隐性变量。
其经常被用在ML和CV领域的数据聚类方面。

该算法的实现步骤主要包含两步：

1.计算期望（E步）
利用对隐藏变量的现有估计值，计算其最大似然估计值；

2.最大化（M步）
最大化在E步上求得的最大似然值来计算参数的值。
M步上找到的参数估计值被用于下一个E步计算中，这个过程不断交替进行

1.2似然函数、极大似然估计

从上文对EM算法的简单介绍中，存在一个最大似然估计的名词，如果不懂这个，可能无法理解EM算法。因此接下来对其进行一个解释。
多数情况下我们是根据已知条件来推算结果，而最大似然估计是已经知道了结果，然后寻求使该结果出现的可能性最大的条件，以此作为估计值。

看到没，概率是根据条件推测结果，而似然则是根据结果反推条件在这种意义上，似然函数可以理解为条件概率的逆反。如下图所示。举个例子来理解什么似然估计

1.2.1 问题描述-调查学生身高分布

我们先假设学校所有学生的身高服从正态分布 $N(\mu ,{\sigma }^2)$。(注意：极大似然估计的前提一定是要假设数据总体的分布，如果不知道数据分布，是无法使用极大似然估计的)，这个分布的均值 $\mu$ 和方差 ${\sigma }^2$ 未知。一旦我们知道了这两个参数，那么学校全体学生的身高分布不就轻松得到了。
所以我们的目标是找到参数 $\mu$ 和方差 ${\sigma }^2$。

但是学校有很多学生，不能挨个统计。所以通常而言，我们通过抽样，由样本来估计总体情况。具体做法是：
假设我们随机抽到了200个学生，并得到里他们的身高数据。并根据这200个数据，估计得到均值 $\mu$ 和方差 ${\sigma }^2$。

换用一种数学的语言：
为了统计学校学生的身高分布，我们独立地按照概率密度$p(x|\theta )$ 抽取了 200 个（身高），组成样本集$X=x_1,x_2,...,x_N$ (其中$x_i$ 表示抽到的第 $i$ 个人的身高，这里 N 就是 200，表示样本个数)，我们想通过样本集 X 来估计出总体的未知参数 $\theta$。这里概率密度$p(x|\theta )$ 服从高斯分布 $N(\mu ,{\sigma }^2)$，其中的未知参数是 $\theta =[\mu ,{\sigma }^2{]}^T$ 。
So，如何估计参数$\theta$呢？

1.2.2参数估计

先回答几个小问题：

Q1：抽到这200个人的概率是多少？
由于每个样本都是独立地从$p(x|\theta )$中抽取的，换句话说这 200 个学生随便捉的，他们之间是没有关系的，即他们之间是相互独立的。假如抽到学生 A（的身高）的概率是$p(x_A|\theta )$,抽到学生B的概率是$p(x_B|\theta )$ ，那么同时抽到男生 A 和男生 B 的概率是 $p(x_A|\theta )\ast p(x_B|\theta )$.同理，我同时抽到这 200 个学生的概率就是他们各自概率的乘积了，即为他们的联合概率，用下式表示：
n=200，为抽样样本数。这个概率反映了，在概率密度函数的参数是$\theta$时，得到 X 这组样本的概率。
上式中等式右侧只有$\theta$是未知数，所以 L 是 $\theta$的函数。

这个函数反映的是在不同的参数$\theta$取值下，取得当前这个样本集的可能性，因此称为参数$\theta$相对于样本集 X 的似然函数（likelihood function），记为$L(\theta )$ 。
对 L 取对数，将其变成连加的，称为对数似然函数，如下式：

取对数，Why?

• 取对数之后累积变为累和，求导更加方便
• 概率累积会出现数值非常小的情况，比如1e-30，造成精度损失

Q2:学校那么多学生，为什么就恰好抽到了这 200 个人 ( 身高) 呢？
从如此多的学生中，恰好抽到了这200位学生，说明在服从同一分布$p(x|\theta )$的情况下，这200位学生的$L(\theta )$值最大，也就是该校学生身高所服从的分布$p(x|\theta )$，恰好使得$L(\theta )$值最大。所以，最终真正的、理想的参数值应该满足以下关系:

这个$\theta$即为所求。
Q3:那么怎么极大似然函数呢？
$L(\theta )$是关于$\theta$的函数，则求$L(\theta )$对所有参数的偏导数，然后让这些偏导数为 0，假设有 n个参数，就有n个方程组成的方程组，那么方程组的解就是似然函数的极值点了，从而得到对应的$\theta$ 了。

1.2.3总结

极大似然估计，只是一种概率论在统计学的应用，它是参数估计的方法之一。说的是已知某个随机样本满足某种概率分布，但是其中具体的参数不清楚，通过若干次试验，观察其结果，利用结果推出参数的大概值。

极大似然估计是建立在这样的思想上：已知某个参数能使这个样本出现的概率极大，我们当然不会再去选择其他小概率的样本，所以干脆就把这个参数作为估计的真实值。

1.2.4求极大似然函数估计的一般步骤

• 写出似然函数，并取对数；
• 求导数，令导数等于0，得到似然方程（组）；
• 求解似然方程（组）得到估计值。

2.EM算法

2.1 问题描述

在前文介绍似然估计时，我们给出所有学生身高服从同一分布的假设，但实际情况确实，男女生的身高服从不一样的分布。如男生服从$N({\mu }_1,{\sigma }_1^2)$分布，女生服从$N({\mu }_2,{\sigma }_2^2)$分布。（注意：EM算法和极大似然估计的前提是一样的，都要假设数据总体的分布，如果不知道数据分布，是无法使用EM算法的）
此时依旧可以按照上文介绍的抽样方法分别得到两个分布的参数。但是，如果抽样完成后，统计表格出现故障，导致数据缺少了性别一栏，此时，我们拿到的200个数据（100个男生+100个女生）因为无法确定性别，不知道每个数据是从哪个分布中来的，也就不能按照前面的方式来进行计算。

2.2问题求解

这个时候，跟前文的不同的地方在于：前文只有一个分布的参数需要顾及。而该问题中，需要估计两个问题。1.每一个数据对应的性别,男生还是女生？；2.数据服从的分布参数这两个问题相互依赖，

• 如果知道了数据对应性别，可以通过极大似然估计得到分布参数；
• 如果知道男女生的分布参数，则可以计算得到每一个数据对应的概率，也就确定了数据的性别。

但是现在二者皆无，就出现了“鸡生蛋、蛋生鸡”的问题。为了解决这个你依赖我，我依赖你的循环依赖问题，总得有一方要先打破僵局，不管了，我先随便整一个值出来，看你怎么变，然后我再根据你的变化调整我的变化，然后如此迭代着不断互相推导，最终就会收敛到一个解（草原上的狼和羊，相生相克）。这就是EM算法的基本思想了。

好了，让我们来解决一下这个问题。

• 先设定男生和女生的身高分布参数(初始值)，例如男生的身高分布为 $N({\mu }_1=172,{\sigma }_1^2=5^2)$ ， 女生的身高分布为 $N({\mu }_2=162,{\sigma }_2^2=5^2)$ ，当然了，刚开始肯定没那么准；
• 然后计算出每个人更可能属于第一个还是第二个正态分布中的（例如，这个人的身高是180，那很明显，他极大可能属于男生），这个是属于Expectation 一步；
• 我们已经大概地按上面的方法将这 200 个人分为男生和女生两部分，我们就可以根据之前说的极大似然估计分别对男生和女生的身高分布参数进行估计（这不变成了极大似然估计了吗？极大即为Maximization）这步称为Maximization；此时已经知道所抽取数据的性别，即确定前面两个问题之一，就可以得到另立一个。
• 然后，当我们更新这两个分布的时候，每一个学生属于女生还是男生的概率又变了，那么我们就再需要调整E步；
• ……如此往复，直到参数基本不再发生变化或满足结束条件为止。

2.3总结

2.3.1相关概念

隐变量：上述例子中的男女性别即是隐变量，也称为隐含参数。一般用Z表示；
模型参数：上述男女身高的分布参数；
观测数据：上述例子中，抽样出来的200个学生身高，即是观测数据，通俗讲就是我们可以确定的数据。一般用Y表示；
完全数据：Z+Y就是完全数据；
不完全数据：仅有一个Y。
所以EM算法就是求解这类带有隐变量Z的问题

2.3.2与K-Means算法对比

EM算法的求解思路与K-Means相似。
EM算法
先猜想隐含参数（EM 算法的 E 步），接着基于观察数据和猜测的隐含参数一起来极大化对数似然，求解我们的模型参数（EM算法的M步)。由于我们之前的隐含参数是猜测的，所以此时得到的模型参数一般还不是我们想要的结果。我们基于当前得到的模型参数，继续猜测隐含参数（EM算法的 E 步），然后继续极大化对数似然，求解我们的模型参数（EM算法的M步)。以此类推，不断的迭代下去，直到模型分布参数基本无变化，算法收敛，找到合适的模型参数。
K-Means算法
在 K-Means 聚类时，每个聚类簇的质心是隐含数据。我们会假设 K 个初始化质心，即 EM 算法的 E 步；然后计算得到每个样本最近的质心，并把样本聚类到最近的这个质心，即 EM 算法的 M 步。重复这个 E 步和 M 步，直到质心不再变化为止，这样就完成了 K-Means 聚类。

3.EM算法推导

3.1相关基础

3.1.1期望

• 对于离散型随机变量X而言，其概率分布为 p i = p { X = x i } p_{i}=p\{X=x_{i}\} pi​=p{X=xi​}，其期望$E(X)$为：
  $E(X)={\sum }_i{x}_{i}p(x_i)$
  其中，$p_i$是权值，满足两个条件:
  $0⩽p_i⩽1$
  ${\sum }_i{p}_i=1$
• 若连续型随机变量X的概率密度函数为$f(x)$，则数学期望$E(X)$为：
  $E(X)={\int }_{-\infty }^{+\infty }xf(x)dx$
  设$Y=g(X)$ ,若X是离散型随机变量，则：
  $E(Y)={\sum }_{i}g(x_i)p_i$
  若X是连续型随机变量，则：
  $E(Y)={\int }_{-\infty }^{+\infty }g(x)f(x)dx$

3.1.2凸函数

设是定义在实数域上的函数，如果对于任意的实数，都有：
那么是凸函数。若不是单个实数，而是由实数组成的向量，此时，如果函数的 Hesse 矩阵是半正定的，即
是凸函数。特别地，如果$f^{{}^{\prime \prime }}⩾0$或者$H^{{}^{\prime \prime }}⩾0$，称为严格凸函数。

3.1.3Jensen不等式

如下图，如果函数$f$ 是凸函数，$x$ 是随机变量，有 0.5 的概率是 a，有 0.5 的概率是 b，$x$ 的期望值就是 a 和 b 的中值了那么：
通俗讲就是：函数的期望大于等于期望的函数。等式成立的条件是随机变量为常数

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

3.2公式推导

对于m个相互独立的样本，$x=(x^{(1)},x^{(2)},...,x^{(m)})$，对应的隐变量$z=(z^{(1)},z^{(2)},...,z^{(m)})$，此时$(x,z)$为完全数据，样本的模型参数为$\theta$，则观察数据$x^{(i)}$的概率为$p(x^{(i)}|\theta )$，完全数据的$(x^{(i)},z^{(i)})$的联合概率为$p(x^{(i)},z^{(i)}|\theta )$.
如果没有隐变量$z$，我们仅需要找到合适的$\theta$极大化对数似然函数即可。
增加隐变量$z$之后，我们目标变成找到合适的$\theta$和$z$使对数似然极大。
常规解法
面对上面这个式子，我们很容易想到分别对$\theta$和$z$求偏导，然后就可以得到使对数似然函数最大的$\theta$和$z$。但是，因为上述函数为$log(连续求和)$这样形式的复合函数，求导后会异常复杂，故而pass.

简单解法
我们可以对$log$里面同乘以一个相同的函数（即隐变量Z的概率分布$Q_i(z^{(i)})$，其概率之和等于1，即${\sum }_z{Q}_i(z^{(i)})=1$），即下图中的（1）式
上图中，式（2）的变化由Jensen不等式得出，因为log函数是凹函数，所以符号相反。如下式所示
其中：
其中，c为常数，则对于任意$i$，有如下公式成立：

方程两边同时累加和，得：
因为上式右边是隐变量$z$的概率分布，所以我们可以得到：
联立上面的式子，得：
可以看出$Q(z)$就是隐变量的分布
见下图，我们固定θ，调整Q(z)使下界J(z,Q)上升至与L(θ)在此点θ处相等（绿色曲线到蓝色曲线），然后固定Q(z)，调整θ使下界J(z,Q)达到最大值（θt到θt+1），然后再固定θ，调整Q(z)……直到收敛到似然函数L(θ)的最大值处的θ*。

3.3EM算法流程

综上所述，我们总结下EM算法的流程。
输入：观测数据$x=(x^{(1)},x^{(2)},...,x^{(m)})$,联合分布$p(x,z|\theta )$，条件分布$p(z|x,\theta )$,极大迭代次数J。

1. 随机初始化模型参数$\theta$的初值为${\theta }^0$;
2. for j from 1 to J:
   E步：计算联合分布的条件概率期望：
   
   M步：极大化$L(\theta )$,得到$\theta$:
   
   重复上述步骤，直至收敛。
3. 输出模型参数$\theta$

3.4EM 算法的另一种理解

3.5EM算法的收敛性思考

EM算法的流程并不复杂，但是还有两个问题需要我们思考：

1） EM算法能保证收敛吗？

2） EM算法如果收敛，那么能保证收敛到全局极大值吗？

3.6EM算法的应用

EM的应用包括：

• 支持向量机的SMO算法
• 混合高斯模型
• K-means
• 隐马尔可夫模型

4.EM算法案例-两硬币模型（文末论文）

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

实验A
已知每次选择的是硬币A or 硬币B，仅仅是不知道模型参数$\theta$,因此直接使用极大似然法求导，
对数似然如下图所示。
分别让上市对${\theta }_A$和${\theta }_B$求偏导，并让偏导为0，可以得到
实验B
由于不知道使用的是硬币A还是硬币B，即使用的硬币为隐变量，因此需要采用EM算法。
第一轮迭代
得到每轮试验选择硬币A还是硬币B的概率，如下表所示
|试验 |硬币A |硬币|

实验次数	P(Z=A)	P(Z=B)
1	0.45	0.55
2	0.8	0.2
3	0.73	0.27
4	0.35	0.65
5	0.65	0.35

后续迭代
通过编写简单python程序，可以顺利计算出后面每一轮迭代过程中，${\theta }_A$和${\theta }_B$可以很容易得出。如下图所示。
第一行代表${\theta }_A$每一次迭代时的变化；
第二行代表${\theta }_B$每一次迭代时的变化：

此外，还可以发现，当${\theta }_A$和${\theta }_B$的初值设为相等时，他们会得到另一组稳定解，且稳定的更快。
第一行代表${\theta }_A$每一次迭代时的变化；
第二行代表${\theta }_B$每一次迭代时的变化：

5.代码

6.参考资料

1.如何通俗理解EM算法
2.你真的了解EM算法吗？
3.Do C B, Batzoglou S. What is the expectation maximization algorithm?[J]. Nature biotechnology, 2008, 26(8): 897-899.
4.李航《统计学习方法》的代码实现
5.知乎【机器学习】EM——期望最大（非常详细）
6.知乎【人人都懂EM算法】
7.B站《白板推导之EM算法》