EM算法原理解释及公式推导

本文参考的是人人都懂EM算法 - August的文章 - 知乎这篇文章

一、极大似然概述

二、EM算法

2.1 EM算法描述

2.2 EM公式推导

三、EM算法案例

一、极大似然概述

假设我们需要调查我们学校学生的身高分布。我们先假设学校所有学生的身高服从正态分布 $N(\mu,\sigma^2)$ 。(注意：极大似然估计的前提一定是要假设数据总体的分布，如果不知道数据分布，是无法使用极大似然估计的)，这个分布的均值 $\mu$ 和方差 $\sigma^2$ 未知，如果我们估计出这两个参数，那我们就得到了最终的结果。那么怎样估计这两个参数呢？

学校的学生这么多，我们不可能挨个统计吧？这时候我们需要用到概率统计的思想，也就是抽样，根据样本估算总体。假设我们随机抽到了 200 个人（也就是 200 个身高的样本数据，为了方便表示，下面“人”的意思就是对应的身高）。然后统计抽样这 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$ 。

那么问题来了怎样估算参数 $\theta$ 呢？

问题一：抽到这 200 个人的概率是多少呢？

n 为抽取的样本的个数，本例中 n=200 ，这个概率反映了，在概率密度函数的参数是 $\theta$ 时，得到 X 这组样本的概率。上式中等式右侧只有 $\theta$ 是未知数，所以 L 是 $\theta$ 的函数。

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

对 L 取对数，将其变成连加的，称为对数似然函数，如下式：

Q：这里为什么要取对数？

取对数之后累积变为累和，求导更加方便
概率累积会出现数值非常小的情况，比如1e-30，由于计算机的精度是有限的，无法识别这一类数据，取对数之后，更易于计算机的识别(1e-30以10为底取对数后便得到-30)。

问题二：学校那么多学生，为什么就恰好抽到了这 200 个人 ( 身高) 呢？

在学校那么学生中，我一抽就抽到这 200 个学生（身高），而不是其他人，那是不是表示在整个学校中，这 200 个人（的身高）出现的概率极大啊，也就是其对应的似然函数 $L(\theta)$ 极大，即

$\hat{\theta}$ 这个叫做 $\theta$ 的极大似然估计量，即为我们所求的值。

问题三：那么怎么极大似然函数？

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

极大似然估计总结

极大似然估计你可以把它看作是一个反推。多数情况下我们是根据已知条件来推算结果，而极大似然估计是已经知道了结果，然后寻求使该结果出现的可能性极大的条件，以此作为估计值。

比如说，

假如一个学校的学生男女比例为 9:1 (条件)，那么你可以推出，你在这个学校里更大可能性遇到的是男生 (结果)；
假如你不知道那女比例，你走在路上，碰到100个人，发现男生就有90个 (结果)，这时候你可以推断这个学校的男女比例更有可能为 9:1 (条件)，这就是极大似然估计。

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

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

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

（1）写出似然函数；

（2）对似然函数取对数，并整理；

（3）求导数，令导数为 0，得到似然方程；

（4）解似然方程，得到的参数。

二、EM算法

2.1 EM算法描述

上面我们先假设学校所有学生的身高服从正态分布 $N(\mu,\sigma^2)$ 。实际情况并不是这样的，男生和女生分别服从两种不同的正态分布，即男生 $\in N(\mu_1,\sigma_1^2)$ ，女生 $\in N(\mu_2,\sigma_2^2)$ ，(注意：EM算法和极大似然估计的前提是一样的，都要假设数据总体的分布，如果不知道数据分布，是无法使用EM算法的)。那么该怎样评估学生的身高分布呢？

简单啊，我们可以随便抽 100 个男生和 100 个女生，将男生和女生分开，对他们单独进行极大似然估计。分别求出男生和女生的分布。

假如某些男生和某些女生好上了，纠缠起来了。咱们也不想那么残忍，硬把他们拉扯开。这时候，你从这 200 个人（的身高）里面随便给我指一个人（的身高），我都无法确定这个人（的身高）是男生（的身高）还是女生（的身高）。用数学的语言就是，抽取得到的每个样本都不知道是从哪个分布来的。那怎么办呢？

这个时候，对于每一个样本或者你抽取到的人，就有两个问题需要估计了，一是这个人是男的还是女的，二是男生和女生对应的身高的正态分布的参数是多少。这两个问题是相互依赖的：

当我们知道了每个人是男生还是女生，我们可以很容易利用极大似然对男女各自的身高的分布进行估计。
反过来，当我们知道了男女身高的分布参数我们才能知道每一个人更有可能是男生还是女生。例如我们已知男生的身高分布为，女生的身高分布为，一个学生的身高为180，我们可以推断出这个学生为男生的可能性更大。

但是现在我们既不知道每个学生是男生还是女生，也不知道男生和女生的身高分布。这就成了一个先有鸡还是先有蛋的问题了。鸡说，没有我，谁把你生出来的啊。蛋不服，说，没有我，你从哪蹦出来啊。为了解决这个你依赖我，我依赖你的循环依赖问题，总得有一方要先打破僵局，不管了，我先随便整一个值出来，看你怎么变，然后我再根据你的变化调整我的变化，然后如此迭代着不断互相推导，最终就会收敛到一个解（草原上的狼和羊，相生相克）。这就是EM算法的基本思想了。

EM的意思是“Expectation Maximization”，具体方法为：

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

总结

上面的学生属于男生还是女生我们称之为隐含参数，女生和男生的身高分布参数称为模型参数。

EM 算法解决这个的思路是使用启发式的迭代方法，既然我们无法直接求出模型分布参数，那么我们可以先猜想隐含参数（EM 算法的 E 步），接着基于观察数据和猜测的隐含参数一起来极大化对数似然，求解我们的模型参数（EM算法的M步)。由于我们之前的隐含参数是猜测的，所以此时得到的模型参数一般还不是我们想要的结果。我们基于当前得到的模型参数，继续猜测隐含参数（EM算法的 E 步），然后继续极大化对数似然，求解我们的模型参数（EM算法的M步)。以此类推，不断的迭代下去，直到模型分布参数基本无变化，算法收敛，找到合适的模型参数。

2.2 EM公式推导

2.2.1 凸函数定义：

设是定义在实数域上的函数，如果对于任意的实数，都有：

那么是凸函数。若不是单个实数，而是由实数组成的向量，此时，如果函数的 Hesse 矩阵是半正定的，即：

是凸函数。特别地，如果或者，称为严格凸函数。

2.2.2 Jensen不等式

如下图，如果函数 $f$ 是凸函数， $x$ 是随机变量，有 0.5 的概率是 a，有 0.5 的概率是 b， $x$ 的期望值就是 a 和 b 的中值了那么：

其中，，这里 a 和 b 的权值为 0.5, $f(a)$ 与 a 的权值相等， $f(b)$ 与 b 的权值相等。

特别地，如果函数 $f$ 是严格凸函数，当且仅当： (即随机变量是常量) 时等号成立。

2.2.3 期望

对于离散型随机变量 X 的概率分布为，数学期望 $E(x)$ 为：

$p_i$ 是权值，满足两个条件。

若连续型随机变量X的概率密度函数为 $f(x)$ ，则数学期望 $E(x)$ 为：

设，若 $X$ 是离散型随机变量，则：

若 $X$ 是连续型随机变量，则：

以上讲的是基础知识，如果不熟悉的同学可以再复习一下概率论。

2.2.4 EM算法推导

对于 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)})$ 的似然函数为。

假如没有隐含变量 $z$ ，我们仅需要找到合适的 $\theta$ 极大化对数似然函数即可：

增加隐含变量 $z$ 之后，我们的目标变成了找到合适的 $\theta$ 和 $z$ 让对数似然函数极大：

这里的 $x_i$ 相当于随便挑一位同学，隐含变量 $z$ 相当于这位同学的性别， $\theta$ 是女生和男生的身高分布参数

不就是多了一个隐变量 $z$ 吗？那我们自然而然会想到分别对未知的 $\theta$ 和 $z$ 分别求偏导，这样做可行吗？

理论上是可行的，然而如果对分别对未知的 $\theta$ 和 $z$ 分别求偏导，由于是边缘概率(建议没基础的同学网上搜一下边缘概率的概念)，而转化为求导后形式会非常复杂（可以想象下复合函数的求导) ，所以很难求解得到 $\theta$ 和 $z$ 。那么我们想一下可不可以将加号从 log 中提取出来呢？我们对这个式子进行缩放如下：

上面第(1)式引入了一个未知的新的分布，满足：

第(2)式用到了 Jensen 不等式 (对数函数是凹函数)：

其中：

也就是说为第 i 个样本，为第 i 个样本对应的权重，那么：

上式我实际上是我们构建了 $L(\theta,z)$ 的下界，我们发现实际上就是的加权求和，由于上面讲过权值累积和为1，因此上式是的加权平均，也是我们所说的期望，这就是Expectation的来历啦。下一步要做的就是寻找一个合适的最优化这个下界(M步)。

假设 $\theta$ 已经给定，那么的值就取决于和了。我们可以通过调整这两个概率使下界逼近的真实值，当不等式变成等式时，说明我们调整后的下界能够等价于了。由 Jensen 不等式可知，等式成立的条件是随机变量是常数，则有：
其中 c 为常数，对于任意 $i$ ，我们得到：

由上面（1）方程，两边同时累加和：

由于。从上面两式，我们可以得到：

其中：

边缘概率公式：

条件概率公式：

从上式可以发现 $Q(z)$ 是已知样本和模型参数下的隐变量分布。

如果 , 则第 (2) 式是我们的包含隐藏数据的对数似然的一个下界。如果我们能极大化这个下界，则也在尝试极大化我们的对数似然。即我们需要极大化下式：

至此，我们推出了在固定参数 $\theta$ 后分布的选择问题，从而建立了的下界，这是 E 步，接下来的M 步骤就是固定后，调整 $\theta$ ，去极大化的下界。

去掉上式中常数的部分，将 $log\frac{p}{q} = logp-logq$ ，而是常数，所以 logq 为0，所以我们需要极大化的对数似然下界为：

2.2.5 EM 算法流程

输入：观察数据 $x=(x^{(1)},x^{(2)}...,x^{(m)})$ ，联合分布，条件分布，极大迭代次数 $J$ 。

1) 随机初始化模型参数 $\theta$ 的初值 $\theta^0$

2) ：

E步：计算联合分布的条件概率期望：

M步：极大化 $L(\theta )$ ,得到 $\theta$ :

重复E、M步骤直到 $\theta$ 收敛

输出：模型参数 $\theta$

三、EM算法案例

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

a 情况：

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

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

b情况：

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

E步：初始化和，计算每个实验中选择的硬币是 A 和 B 的概率，例如第一个实验中选择 A 的概率为：

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

M步：求出似然函数下界， $y_j$ 代表第 j 次实验正面朝上的个数， $\mu_j$ 代表第 j 次实验选择硬币 A 的概率， $1-\mu_j$ 代表第 j 次实验选择硬币B的概率。

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

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

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

第二轮迭代：基于第一轮EM计算好的 , 进行第二轮 EM，计算每个实验中选择的硬币是 A 和 B 的概率（E步），然后在计算M步，如此继续迭代......迭代十步之后