机器学习之从极大似然估计到最大熵原理以及EM算法详解

一、极大似然估计

极大似然估计是建立在极大似然原理的基础上的一个统计方法，极大似然原理的直观想法是，一个随机试验如有若干个可能的结果A，B，C，… ，若在一次试验中，结果A出现了，那么可以认为实验条件对A的出现有利，也即出现的概率P(A)较大。极大似然原理的直观想法我们用下面例子说明。设甲箱中有99个白球，1个黑球；乙箱中有1个白球．99个黑球。现随机取出一箱，再从抽取的一箱中随机取出一球，结果是黑球，这一黑球从乙箱抽取的概率比从甲箱抽取的概率大得多，这时我们自然更多地相信这个黑球是取自乙箱的。一般说来，事件A发生的概率与某一未知参数 θ $\theta$有关， θ $\theta$取值不同，则事件A发生的概率P(A| θ $\theta$)也不同，当我们在一次试验中事件A发生了，则认为此时的 θ $\theta$值应是t的一切可能取值中使P(A| θ $\theta$)达到最大的那一个，极大似然估计法就是要选取这样的t值作为参数t的估计值，使所选取的样本在被选的总体中出现的可能性为最大。

假设总体分布为f(x, θ $\theta$)， X1,X2,X3...Xn $X_{1},X_{2},X_{3}...X_{n}$为该总体采样取得的样本。因为 X1,X2,X3...Xn $X_{1},X_{2},X_{3}...X_{n}$是独立同分布的，于是，它的联合密度函数为：

$$L(X_{1},X_{2}...X_{n};\theta _{1},\theta _{2}...\theta _{k})=\prod_{i=1}^{n}f(x_{i};\theta _{1},\theta _{2}...\theta _{k})$$

在上面这个式子中， θ $\theta$被看做是确定但是未知的参数，并且应为样本已经存在，所以 x1,x2,x3...xn $x_{1},x_{2},x_{3}...x_{n}$也是固定的。因此 L(x,θ) $L(x,\theta )$是关于 θ $\theta$的函数，即似然函数。求解参数 θ $\theta$的值，使似然函数取得极大值，这就是极大似然估计。

在实践中，由于求导数的需要，往往将似然函数取对数，得到对数似然函数；如果似然函数可导，那么就可以通过求导数的方式得到驻点，从而算出极大值。

logL(θ1,θ2...θk)=∑ni=1logf(xi;θ1,θ2...θk) $logL(\theta _{1},\theta _{2}...\theta _{k})=\sum_{i=1}^{n}logf(x_{i};\theta _{1},\theta _{2}...\theta _{k})$
∂L(θ)∂θi=0,i=1,2....k $\frac{\partial L(\theta )}{\partial \theta _{i}}=0,i=1,2....k$

求极大似然估计量的一般步骤：
（1）写出似然函数；
（2）对似然函数取对数，并整理；
（3）求导数；
（4）解似然方程。

例子：

极大似然估计的特点：

1.比其他估计方法更加简单；
2.收敛性：无偏或者渐近无偏，当样本数目增加时，收敛性质会更好；
3.如果假设的类条件概率模型正确，则通常能获得较好的结果。但如果假设模型出现偏差，将导致非常差的估计结果。

---

二、最大熵原理

最大熵原理是一种选择随机变量统计特性最符合客观情况的准则，也称为最大信息原理。随机量的概率分布是很难测定的，一般只能测得其各种均值（如数学期望、方差等）或已知某些限定条件下的值（如峰值、取值个数等），符合测得这些值的分布可有多种、以至无穷多种，通常，其中有一种分布的熵最大。选用这种具有最大熵的分布作为该随机变量的分布，是一种有效的处理方法和准则。这种方法虽有一定的主观性，但可以认为是最符合客观情况的一种选择。在投资时常常讲不要把所有的鸡蛋放在一个篮子里，这样可以降低风险。在信息处理中，这个原理同样适用。在数学上，这个原理称为最大熵原理。

那么，到底什么是熵呢？简单来说，熵是对平均不确定性的度量：

H(X)=−∑xϵXp(x)lnp(x) $H(X)=-\sum_{x\epsilon X}^{ }p(x)lnp(x)$

由以上公式可知，熵是随机变量不确定性的度量，不确定性越大，熵也越大；当随机变量变成一个确定的值时，熵就变成了0。需要指出的是均匀分布是“最不确定”的分布。

最大熵的一般模型：

max(pϵP)H(Y|X)=−∑(x,y)p(x,y)logp(y|x) $max(p\epsilon P)H(Y|X)=-\sum_{(x,y)}^{ }p(x,y)logp(y|x)$

其中P={p|p是X上满足条件的概率分布}

例子：

假设随机变量X有5个取值{A,B,C,D,E},且满足条件P(A)+P(B)=0.3且P(A)+P(B)+P(C)+P(D)+P(E)=1。求最大熵模型。

为了方便,分别用 y1 $y_{1}$~ y5 $y_{5}$表示A~E，于是最大熵模型的最优化问题是:

min−H(p)=∑5i=1p(yi)logp(yi) $min-H(p)=\sum_{i=1}^{5}p(y_{i})logp(y_{i})$

s.t.p(y1)+p(y2)=p˜(y1)+p˜(y2)=310 $s.t. p(y_{1})+p(y_{2})=\widetilde{p}(y_{1})+\widetilde{p}(y_{2})=\frac{3}{10}$

∑5i=1p(yi)=∑5i=1p˜(yi)=1 $\sum_{i=1}^{5}p(y_{i})=\sum_{i=1}^{5}\widetilde{p}(y_{i})=1$

引进拉格朗日乘子w0和w1,定义拉格朗日函数如下:

$$L(p,w)=p(y_{i})log p(y_{i})+w_{1}(p(y_{1})+p(y_{2})-\frac{3}{10})+w_{0}(\sum_{i=1}^{5}p(y_{i})-1)$$

根据拉格朗日对偶性,可以通过求解对偶最优化问题得到原始最优化问题的解。所以求解max min L(p,w)首先需要求解关于p的极小化问题。为此需要固定w0和w1。求偏导数:

∂L(p,w)∂p(y1)=1+logp(y1)+w1+w0 $\frac{\partial L(p,w)}{\partial p(y_{1})}=1+log p(y_{1})+w_{1}+w_{0}$

∂L(p,w)∂p(y2)=1+logp(y2)+w1+w0 $\frac{\partial L(p,w)}{\partial p(y_{2})}=1+log p(y_{2})+w_{1}+w_{0}$

∂L(p,w)∂p(y3)=1+logp(y3)+w0 $\frac{\partial L(p,w)}{\partial p(y_{3})}=1+log p(y_{3})+w_{0}$

∂L(p,w)∂p(y4)=1+logp(y4)+w0 $\frac{\partial L(p,w)}{\partial p(y_{4})}=1+log p(y_{4})+w_{0}$

∂L(p,w)∂p(y5)=1+logp(y5)+w0 $\frac{\partial L(p,w)}{\partial p(y_{5})}=1+log p(y_{5})+w_{0}$

令上面的五个偏导数都等于0，可求得：

$$p(y_{1})=p(y_{2})=e^{-w_{1}-w_{0}-1}，p(y_{3})=p(y_{4})=p(y_{5})=e^{-w_{0}-1}$$

把 p(y1)，p(y2)，p(y3)，p(y4)，p(y5) $p(y_{1})，p(y_{2})，p(y_{3})，p(y_{4})，p(y_{5})$代入到 L(p,w) $L(p,w)$中，再求 L(p,w) $L(p,w)$关于w的极大化问题:

$$max L(p_{w},w)=-2 e^{-w_{1}-w_{0}-1}-3e^{-w_{0}-1}-\frac{3}{10}w_{1}-w_{0}$$

分别对 w0 $w_{0}$和 w1 $w_{1}$求偏导,并令其等于0,可以得到：

$$e^{-w_{1}-w_{0}-1}=\frac{3}{20}，e^{-w_{0}-1}=\frac{7}{30}$$

解得：

$$p(y_{1})=p(y_{2})=\frac{3}{20}，p(y_{3})=p(y_{4})=p(y_{5})=\frac{7}{30}$$

总结：

最大熵模型在分类方法里算是比较优的模型，但是由于它的约束函数的数目一般来说会随着样本量的增大而增大，导致样本量很大的时候，对偶函数优化求解的迭代过程非常慢，scikit-learn甚至都没有最大熵模型对应的类库。但是理解它仍然很有意义，尤其是它和很多分类方法都有千丝万缕的联系。　

优点

a) 最大熵统计模型获得的是所有满足约束条件的模型中信息熵极大的模型,作为经典的分类模型时准确率较高。

b) 可以灵活地设置约束条件，通过约束条件的多少可以调节模型对未知数据的适应度和对已知数据的拟合程度。

缺点

由于约束函数数量和样本数目有关系，导致迭代过程计算量巨大，实际应用比较难。

---

三、EM算法

首先举一个例子：

现在一个班里有50个男生，50个女生，且男生站左，女生站右。我们假定男生的身高服从正态分布： N(μ1,σ21) $N(\mu _{1},\sigma _{1}^{2})$，女生的身高则服从另一个正态分布： N(μ2,σ22) $N(\mu _{2},\sigma _{2}^{2})$ 。这时候我们可以用极大似然法（MLE），分别通过这50个男生和50个女生的样本来估计这两个正态分布的参数。但现在我们让情况复杂一点，就是这50个男生和50个女生混在一起了。我们拥有100个人的身高数据，却不知道这100个人每一个是男生还是女生。这时候情况就有点尴尬，因为通常来说，我们只有知道了精确的男女身高的正态分布参数我们才能知道每一个人更有可能是男生还是女生。但从另一方面去考量，我们只有知道了每个人是男生还是女生才能尽可能准确地估计男女各自身高的正态分布的参数。这个时候有人就想到我们必须从某一点开始，并用迭代的办法去解决这个问题：我们先设定男生身高和女生身高分布的几个参数（初始值），然后根据这些参数去判断每一个样本（人）是男生还是女生，之后根据标注后的样本再反过来重新估计参数。之后再多次重复这个过程，直至稳定。这个算法也就是EM算法。

假设我们有一个样本集{ x1,x2...xm $x^{1},x^{2}...x^{m}$}，包含m个独立的样本。但每个样本i对应的类别 zi $z^{i}$是未知的（相当于聚类），也即隐含变量。故我们需要估计概率模型p(x,z)的参数θ，但是由于里面包含隐含变量z，所以很难用最大似然求解，但如果z知道了，那我们就很容易求解了。

下面我们通过建立极大似然函数来建立目标函数：

$$l(\theta )=\sum_{i=1}^{m}log p(x;\theta )=\sum_{i=1}^{m}log \sum_{z}^{ }p(x,z;\theta ))$$
进一步计算可得：
$$\begin{eqnarray}\sum_{i}^{ }log p(x^{(i)};\theta)&=&\sum_{i}^{ }log\sum_{z^{(i)}}^{ } p(x^{(i)},z^{(i)};\theta)\quad\quad\qquad\qquad\qquad(1)\\&=&\sum_{i}^{ }log\sum_{z^{(i)}}^{ } Q_{i}(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta )}{Q_{i}(z^{(i)})}\qquad\qquad(2)\\&\geq&\sum_{i}^{ }\sum_{z^{(i)}}^{ }Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta )}{Q_{i}(z^{(i)})}\qquad\qquad(3)\end{eqnarray}$$

本质上我们是需要最大化（1）式（对（1）式，我们回忆下联合概率密度下某个变量的边缘概率密度函数的求解，注意这里z也是随机变量。对每一个样本i的所有可能类别z求等式右边的联合概率密度函数和，也就得到等式左边为随机变量x的边缘概率密度），也就是似然函数，但是可以看到里面有“和的对数”，求导后形式会非常复杂（自己可以想象下log(f1(x)+ f2(x)+ f3(x)+…)复合函数的求导），所以很难求解得到未知参数z和θ。那OK，我们可否对（1）式做一些改变呢？我们看（2）式，（2）式只是分子分母同乘以一个相等的函数，还是有“和的对数”啊，还是求解不了，那为什么要这么做呢？咱们先不管，看（3）式，发现（3）式变成了“对数的和”，那这样求导就容易了。我们注意点，还发现等号变成了不等号，为什么能这么变呢？这就是Jensen不等式的大显神威的地方。

Jensen不等式

设f是定义域为实数的函数，如果对于所有的实数x。如果对于所有的实数x，f(x)的二次导数大于等于0，那么f是凸函数。当x是向量时，如果其hessian矩阵H是半正定的，那么f是凸函数。如果只大于0，不等于0，那么称f是严格凸函数。

Jensen不等式表述如下：

如果f是凸函数，X是随机变量，那么：E[f(X)]>=f(E[X])。
特别地，如果f是严格凸函数，当且仅当X是常量时，上式取等号。Jensen不等式应用于凹函数时，不等号方向反向。

一般意义上的Jensen不等式可以参考：百度词条：Jensen不等式

回到公式（2），因为f(x)=log x为凹函数（其二次导数为 −1x2<0 $-\frac{1}{x^{2}}<0$）。

（2）式中 ∑z(i)Qi(z(i))[p(x(i),z(i);θ)Qi(z(i))] $\sum_{z^{(i)}}^{ } Q_{i}(z^{(i)})[\frac{p(x^{(i)},z^{(i)};\theta )}{Q_{i}(z^{(i)})}]$是 p(x(i),z(i);θ)Qi(z(i)) $\frac{p(x^{(i)},z^{(i)};\theta )}{Q_{i}(z^{(i)})}$的期望，（考虑到 E(X)=∑x∗p(x) $E(X)=\sum x*p(x)$，f(x)是x的函数，则 E(f(X))=∑f(x)∗p(x) $E(f(X))=\sum f(x)*p(x)$，又 ∑zQi(z(i))=1 $\sum_{z}^{ } Q_{i}(z^{(i)})=1$所以就可以得到公式（3）的不等式了。

OK，到这里，现在式（3）就容易地求导了，但是式（2）和式（3）是不等号啊，式（2）的最大值不是式（3）的最大值啊，而我们想得到式（2）的最大值，那怎么办呢？

现在我们就需要一点想象力了，上面的式（2）和式（3）不等式可以写成：似然函数 L(θ) $L(\theta )$>=J(z,Q)，那么我们可以通过不断的最大化这个下界J，来使得 L(θ) $L(\theta )$不断提高，最终达到它的最大值。

见上图，我们固定 θ $\theta$，调整Q(z)使下界J(z,Q)上升至与 L(θ) $L(\theta )$在此点 θ $\theta$处相等（绿色曲线到蓝色曲线），然后固定Q(z)，调整θ使下界J(z,Q)达到最大值（ θt到θt+1 $\theta ^{t}到\theta ^{t+1}$），然后再固定 θ $\theta$，调整Q(z)……直到收敛到似然函数 L(θ) $L(\theta)$的最大值处的 θ∗ $\theta ^{*}$。这里有两个问题：什么时候下界J(z,Q)与 L(θ) $L(\theta)$在此点 θ $\theta$处相等？为什么一定会收敛？

首先第一个问题，在Jensen不等式中说到，当自变量X是常数的时候，等式成立。而在这里，即：

$$\frac{p(x^{(i)},z^{(i)};\theta )}{Q_{i}(z^{(i)})}=c$$

再推导下，由于 ∑zQi(z(i))=1 $\sum_{z}^{ }Q_{i}(z^{(i)})=1$（因为Q是随机变量z(i)的概率密度函数），则可以得到：分子的和等于c（分子分母都对所有z(i)求和：多个等式分子分母相加不变，这个认为每个样例的两个概率比值都是c），则：

$$Q_{i}(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta )}{\sum_{z}^{ }p(x^{(i)},z;\theta )}=\frac{p(x^{(i)},z^{(i)};\theta )}{p(x^{(i)};\theta )}=p(z^{(i)}|x^{(i)};\theta )$$

至此，我们推出了在固定参数后，使下界拉升的Q(z)的计算公式就是后验概率，解决了Q(z)如何选择的问题。这一步就是E步，建立 L(θ) $L(\theta)$的下界。接下来的M步，就是在给定Q(z)后，调整 θ $\theta$，去极大化 L(θ) $L(\theta)$的下界J（在固定Q(z)后，下界还可以调整的更大）。

EM算法整体框架：

E步(第一步)：
如果是首次运行，则根据先验知识给定一个θ；如果不是，则这个θ就是M步求出来的。利用这个θ和样本x，求出隐变量z的条件概率，即Q。

M步(第二步)： 将E步求出的Q带入式1后求出θ的最大值。 重复上面两步，直到收敛。

详细推导过程可以参考：（EM算法）The EM Algorithm

优缺点：

要有一些训练数据，再定义一个最大化函数，采用EM算法，利用计算机经过若干次迭代，就可以得到所需的模型。EM算法是自收敛的分类算法，既不需要事先设定类别也不需要数据见的两两比较合并等操作。缺点是当所要优化的函数不是凸函数时，EM算法容易给出局部最佳解，而不是最优解。

关于怎么更通俗易懂地理解EM算法可以参考以下链接：

https://www.zhihu.com/question/27976634/answer/39132183

https://www.zhihu.com/question/27976634/answer/153567695

---

理解EM算法的九层境界

参考资料：

1. 从最大似然到EM算法浅解

2. 百度文库：极大似然估计