EM算法(Expectation maximization algorithm)

• • 极大似然估计
  • 极大似然估计定义
    • 求解过程
  • EM算法
    • EM算法和极大似然估计的区别
    • 鸡生蛋，蛋生鸡问题
    • EM算法思想
    • 三硬币模型
    • jensen不等式
    • Q函数
    • EM算法的推导
  • K-means中EM思想
  • 参考资料

纠结了好几天，总算搞清楚了EM算法的大概。因此写下这篇博客做个笔记，由于这方面懂得不是很多，可能存在理解错误的地方，欢迎大家指正，好了闲话不多说。

极大似然估计

极大似然估计定义

在正式介绍EM算法之前，我们需要先来了解一下最大似然估计。这个我们应该都在概率论中学过，其实思想比较简单，而且我们在生活中也经常用到，举一个简单的例子：
$$ 某位同学与一位猎人一起外出打猎，一只野兔从前方窜过。只听一声枪响，野兔应声到下，如果要你推测，这一发命中的子弹是谁打的？你就会想，只发一枪便打中，由于猎人命中的概率一般大于这位同学命中的概率，看来这一枪是猎人射中的。这便是最大似然的思想，看起来是不是非常简单。下面来看看极大似然的定义。

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

求解过程

现在我们以一个正态分布为例，假设有一组样本

D=(x1,x2,x3....xn) $$D=(x_1,x_2,x_3....x_n)$$

我们知道样本 x1,x2,...xn $x_1,x_2,...x_n$独立同分布于一个正太分布函数：

f(x)=12π−−√δe−(x−μ)22δ2 $$f(x)=\frac{1}{\sqrt{2\pi }\delta }{e}^{-\frac{(x-\mu {)}^2}{2{\delta }^2}}$$

并且 δ $\delta$是已知的，而 μ $\mu$未知，那么最大似然需要最什么呢？求出当均值为多少时，产生这种采样数据的概率最大。
我们令似然函数 l(θ)=P(D|θ) $l(\theta )=P(D|\theta )$，这里我们解释下什么是 θ $\theta$，其实 θ $\theta$表示的是是当前概率最大的最大似然函数的模型，什么意思呢？即取到当前样本最大的时候对应的函数参数 (μ,δ2) $(\mu ,{\delta }^2)$。在这里其实等同于 μ $\mu$。

由于样本之间是独立同分布所以

l(θ)=∏i=1nf(xi|θ) $$l(\theta )=\underset{i=1}{\overset{n}{\prod }}f(x_i|\theta )$$

下一步要做的便是找到一个和是的 θ^ $\hat{\theta }$使得 l(θ) $l(\theta )$最大，即

θ^=argmaxl(θ) $$\hat{\theta }=argmaxl(\theta )$$

具体做法，对 θ $\theta$求导，然后解出最大值。

∂L(θ)∂θ=0 $$\frac{\mathrm{\partial }L(\theta )}{\mathrm{\partial }\theta }=0$$

然后求解 θ $\theta$。

具体计算过程可以看看极大似然估计详解，可以更好的理解 θ $\theta$以及整个流程。

所以计算过程总结下来就是

（1）写出似然函数；

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

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

（4）解似然方程，得到的参数即为所求；

EM算法

既然我们前面讲到极大似然估计，那么EM算法到底和他有什么关系呢？以生活中送快递为例子

EM算法和极大似然估计的区别

极大似然估计面临的情况
$$一个快递员给你送货。若他到你家只有一条路（结果的实现依赖一个概率分布），但却不知道这条路今天修不修路（不知道该概率分布的参数），修路的话今天快递员就没法送货，若结果是快递员送到货了，那这条路修路了没？答案很明显：没修路。

EM算法面临的情况
$$快递员到你家的路有N条（结果的实现依赖多个概率分布），但快递员只会选择一条路，即，今天他不会选择第二条路，若他选择的路修路，那他就不给你送货了，即使这会而让你暴跳如雷。问：如果今天快递员送到货了，则他选择的哪条路？那条路修路了吗？对于这个，因为你不知道他选择的哪条路（他把货送到就走了，根本不给你问他话的时间），所以你唯一能做的就是估计出这N条路被他选择的N个概率（即：每个概率分布的权值），然后在根据极大似然估计来得出：这条路没修路（求出每个概率分布的参数）。

一句话总结就是极大似然估计是知道概率分布，不知道参数，现在需要通过求解参数使得当前观测值的可能性最大，而EM算法是知道观测值属于哪一个概率分布，也不知道参数

鸡生蛋，蛋生鸡问题

对于EM算法，我们既不知道样本属于那个概率分布，也不知道具体的参数，要求在什么参数下会使观测概率最大？这样就带来了一个问题，要知道属于哪个分布，必须要知道具体的参数。然而参数是未知的。要求极大参数，知道概率分布是前提。两个相互依赖，但是又都是未知的，就造成了鸡生蛋蛋生鸡的问题。

EM算法思想

面对上面的问题，EM算法怎么做的呢？举例说明：现在你有一堆糖果，现在有两个盘子，你需要将糖果均分到两个盘子中，你又不想一个一个的数，嫌太麻烦了。那可以先随机把糖果分成两堆，分别放到两个盘子当中，然后用手分别拿起盘子掂量一下分量，判断哪一个重一些，然后将重的那一盘糖果中那一部分出来放入轻的当中，重复这个过程直到两边分量感觉起来差不多。

EM算法就是这样，假设我们想估计知道A和B两个参数，在开始状态下二者都是未知的，但如果知道了A的信息就可以得到B的信息，反过来知道了B也就得到了A。可以考虑首先赋予A某种初值，以此得到B的估计值，然后从B的当前值出发，重新估计A的取值，这个过程一直持续到收敛为止。

上面第一步赋初值对应的EM算法中的E步，求期望(expcetation)，后面的迭代表示的M步，求极大值。

三硬币模型

假设有三枚硬币，分别记作A,B,C。这些硬币正面出现的概率分别为 π,p,q $\pi ,p,q$。现在进行抛硬币实现：先抛A，根据A的结果抛出B或C，正面选B，反面选C。然后掷出所选的硬币。出现正面记为1，反面为0。独立重复n次(n=10)，观测结果如下:

1,1,0,1,0,0,1,0,1,1 $$1,1,0,1,0,0,1,0,1,1$$

假设只能观测到B和C的结果，不能观测A的结果，求三枚硬币出现正面的概率 π,p,q $\pi ,p,q$
三硬币模型可以写作

p(y|θ)=∑zp(y,z|θ)=∑zp(y|z,θ)p(z|θ)(1) $$p(y|\theta )=\underset{z}{\sum }p(y,z|\theta )=\underset{z}{\sum }p(y|z,\theta )p(z|\theta )(1)$$

p(y|θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y(2) $$p(y|\theta )=\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}(2)$$

其中y表示观测变量，即B、C的结果：1或0。随机变量z是隐含变量，即A最后的结果，我们是无法观测。 θ=(π,p,q) $\theta =(\pi ,p,q)$是模型参数。对公示(2)的解释， p(y|θ) $p(y|\theta )$表示在 θ $\theta$下 y $y$的概率，假设 y=1 $y=1$那么公式则为 p(y|θ)=πpy+(1−π)qy $p(y|\theta )=\pi p^y+(1-\pi )q^y$,当 y=0 $y=0$的时候， p(y|θ)=π(1−p)y+(1−π)(1−q)y $p(y|\theta )=\pi (1-p{)}^y+(1-\pi )(1-q{)}^y$，由于 y=1or0 $y=1or0$，将两种结果同一即得到(2)

将上述观测数据表示为 Y=(y1,y2...yn) $Y=(y_1,y_2...y_n)$，为观测数据表示为 Z=(Z1,Z2,....Zn) $Z=(Z_1,Z_2,....Z_n)$。则观测数据的似然函数为：

P(Y|θ)=∑zp(Y,Z|θ)p(Z|θ) $$P(Y|\theta )=\underset{z}{\sum }p(Y,Z|\theta )p(Z|\theta )$$

即：

P(Y|θ)=∏j=1n[πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj] $$P(Y|\theta )=\underset{j=1}{\overset{n}{\prod }}[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}]$$

模型参数 θ=(π,p,q) $\theta =(\pi ,p,q)$的极大似然估计为:

θ^=argmaxlogP(Y|θ) $$\hat{\theta }=argmaxlogP(Y|\theta )$$

这个问题就没有解析解，只有通过迭代方式求解。EM算法就是可以用于求解这个问题的迭代算法。EM算法首先选取参数的初始值，记作 θ(0)=(π(0),p(0),q(0)) ${\theta }^{(0)}=({\pi }^{(0)},p^{(0)},q^{(0)})$，然后通过下面的步骤迭代计算参数的估计值，直到收敛为止。 θ(i)=(π(i),p(i),q(i)) ${\theta }^{(i)}=({\pi }^{(i)},p^{(i)},q^{(i)})$标识的是第i次迭代后的模型参数。第 i+1 $i+1$次的迭代我们可以这样表示：

E步：计算模型参数 π(i),p(i),q(i) ${\pi }^{(i)},p^{(i)},q^{(i)}$下观测数据 yj $y_j$来自硬币B的概率：

μ(i+1)j=π(i)(p(i))yj(1−p(i))1−yjπ(i)p(i))yj(1−p(i))1−yj+(1−π(i))(q(i))yj(1−q(i))(1−yj)(3) $${\mu }_j^{(i+1)}=\frac{{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}}{{\pi }^{(i)}{p}^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}+(1-{\pi }^{(i)})(q^{(i)}{)}^{y_j}(1-q^{(i)}{)}^{(1-y_j)}}(3)$$

这一步其实就是前面公式(2)的第i次迭代，计算出来自B的概率之后接下来需要对参数重新估值。

M步：计算参数模型的新估值。

π(i+1)=1n∑j=1nμ(i+1)j(4) $${\pi }^{(i+1)}=\frac{1}{n}\underset{j=1}{\overset{n}{\sum }}{\mu }_j^{(i+1)}(4)$$

p(i+1)=1n∑nj=1μ(i+1)yj∑nj=1μ(i+1)(5) $$p^{(i+1)}=\frac{1}{n}\frac{\underset{j=1}{\overset{n}{\sum }}{\mu }^{(i+1)}{y}_j}{\underset{j=1}{\overset{n}{\sum }}{\mu }^{(i+1)}}(5)$$

q(i+1)=∑nj=1(1−μ(i+1))yj∑nj=1(1−μ(i+1))(6) $$q^{(i+1)}=\frac{\underset{j=1}{\overset{n}{\sum }}(1-{\mu }^{(i+1)})y_j}{\underset{j=1}{\overset{n}{\sum }}(1-{\mu }^{(i+1)})}(6)$$

这四个公式看起来有点复杂，其实理解起来没有那么难。公式(3)和前面公式(2)是一样的，这里就不在赘述。现在看看公式(4), π $\pi$其实就是观测是 yj $y_j$来自B的概率，所以只需要将公式(3)求均值就行了。公式(5)表示的观测值 yj $y_j$来自B并且为正面的概率，即求在观测值来自B的条件下，观测为正面的条件概率。关于这个具体计算，可以参考李航统计学习方法里面。需要注意的一点是：EM算法的参数估计值与选取的初始值有关。

jensen不等式

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

Jensen不等式表述如下：

如果f是凸函数，X是随机变量，那么： E[f(X)]≥f(E[X]) $E[f(X)]\ge f(E[X])$

特别地，如果f是严格凸函数，当且仅当X是常量时，上式取等号。

Q函数

定义： 完全数据的对数似然函数 logP(Y,Z|θ) $logP(Y,Z|\theta )$关于在给定的观测数据 Y $Y$和当前参数 θ(i) ${\theta }^{(i)}$下对为观测数据Z的条件概率分布 P(Z|Y,θ(i)) $P(Z|Y,{\theta }^{(i)})$的期望，即

Q(θ,θ(i))=EZ[logP(Y,Z|θ)|Y,θ(i)] $$Q(\theta ,{\theta }^{(i)})=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]$$

=∑ZP(Z|Y,θ(i))logP(Y,Z|θ) $$=\underset{Z}{\sum }P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )$$

什么意思呢？其实简化之后及为 Q(θ,θ(i)) $$=\underset{Z}{\sum }P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )$$

什么意思呢？其实简化之后及为 Q(θ,θ(i))=∑ZP(Z)logP(Y,Z) $Q(\theta ,{\theta }^{(i)})=\underset{Z}{\sum }P(Z)logP(Y,Z)$

EM算法的推导

现在我们来看看EM算法的一般推导，一般的，用 Y $Y$表示观测随机变量的数据，Z表示隐随机变量的数据， YY $Y$和 Z $Z$连在一起被称为完全数据，观测数据 Y $Y$又称为不完全观测数。给定观测数据Y，其概率分布为 P(Y|θ) $P(Y|\theta )$， θ $\theta$为模型参数，那么不完全数据 Y $Y$的对数似然函数为 Y