原理

极大似然估计

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

比如说，

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

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

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

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

（1）写出似然函数；

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

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

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

EM算法

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

概率模型有时既含有观测变量（observable variable），又含有隐变量或潜在变量（latent variable）。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数。但是，当模型含有隐变量时，就不能简单地使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法。

EM算法和极大似然估计的前提是一样的，都要假设数据总体的分布，如果不知道数据分布，是无法使用EM算法的)

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

三硬币模型

假设有3枚硬币，分别记作A,B,C。这些硬币正面出行的概率分别是\(\pi\),p和q。进行如下抛硬币试验：先抛硬币A，根据其结果选出硬币B或硬币C，正面选硬币B，反面选硬币C；然后抛选出的硬币，抛硬币的结果，出现正面记作1，出现反面记作0；独立地重复n次试验（这里，n=10），观测结果如下：

\[1,1,0,1,0,0,1,0,1,1\]

假设只能观测到抛硬币的结果，不能观测抛硬币的过程。问如何估计三硬币正面出现的概率，即三硬币模型的参数。

解 三硬币模型可以写作

\[\begin{align} P(y \mid \theta) &= \sum_z P(y,z \mid \theta) = \sum_z P(z \mid \theta) P(y \mid z, \theta) \\ &= \pi p^y (1-p)^{1-y}) + (1-\pi)q^y(1-q)^{1-y} \\ \end{align}\]

这里，随机变量y是观测变量，表示一次试验观测的结果是1或0；随机变量z是隐变量，表示未观测到的抛硬币A的结果；\(\theta=(\pi,p,q)\)是模型参数。这一模型是以上数据的生成模型。注意，随机变量y的数据可以观测，随机变量z的数据不可观测。

将观测数据表示为\(Y=(Y_1,Y_2,...,Y_n)^T\)，未观测数据表示为\(Z=(Z_1,Z_2,...,Z_n)^T\)，则观测数据的似然函数为

\[P(Y \mid \theta) = \sum_Z P(Z \mid \theta) P(Y \mid Z, \theta)\]

即

\[P(Y \mid \theta) = \prod_{j=1}^n [\pi p^{y_j} (1-p)^{1-y_j}) + (1-\pi)q^{y_j} (1-q)^{1-y_j}]\]

考虑模型参数\(\theta = (\pi, p,q)\)的极大似然估计，即

\[\hat{\theta} = arg max_\theta log P(Y \mid \theta)\]

这个问题没有解析解，只有通过迭代的方法求解。EM算法就是可以用于求解这个问题的一种迭代算法。

EM算法首先选取参数的初值，记作\(\theta^{(0)}=(\pi^{(0)}, p^{(0)},q^{(0)})\)，然后通过下面的步骤迭代计算参数的估计值，直到收敛为止。第i次迭代参数估计值为\(\theta^{(i)}=(\pi^{(i)}, p^{(i)},q^{(i)})\)。EM算法的第i+1次迭代如下。

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

\[\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)}}\]

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

\[\pi^{(i+1)} = \frac{1}{n} \sum_{j=1}^n \mu_j^{(i+1)}\] \[p^{(i+1)} = \frac{\sum_{j=1}^n \mu_j^{(i+1)}y_j}{\sum_{j=1}^n \mu_j^{(i+1)}}\] \[q^{(i+1)} = \frac{\sum_{j=1}^n (1-\mu_j^{(i+1)})y_j}{\sum_{j=1}^n (1-\mu_j^{(i+1)})}\]

进行数字计算，假设模型参数的初值取为

\[\pi^{(0)} = 0.5, p^{(0)}=0.5, q^{(0)}=0.5\]

由E步得到，\(\mu_j^{(1)} = 0.5, j=1,2,...,10\)

由M步得到，

\[\pi^{(1)} = 0.5, p^{(1)}=0.6, q^{(1)}=0.6\]

迭代，由E步得到，\(\mu_j^{(2)} = 0.5, j=1,2,...,10\)

由M步得到，

\[\pi^{(2)} = 0.5, p^{(2)}=0.6, q^{(2)}=0.6\]

于是得到模型参数\(\theta\)的极大似然估计：

\[\hat{\pi} = 0.5, \hat{p} = 0.6, \hat{q}=0.6\]

\(\pi=0.5\)表示硬币A是均匀的，这一结果容易理解。

如果取初值\(\pi^{(0)} = 0.4, p^{(0)}=0.6, q^{(0)}=0.7\)，那么得到的模型参数的极大似然估计是\(\hat{\pi} = 0.4064, \hat{p} = 0.5368, \hat{q}=0.6432\)。这就是说，EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。

一般地，用Y表示观测随机变量的数据，Z表示隐随机变量的数据。Y和Z连在一起称为完全数据（complete-data），观测数据Y又称为不完全数据（incomplete-data）。假设给定观测数据Y，其概率分布是\(P(Y \mid \theta)\)，其中\(\theta\)是需要估计的模型参数，那么不完全数据Y的似然函数是\(P(Y \mid \theta)\)，对数似然函数\(L(\theta) = log P(Y \mid \theta)\)；假设Y和Z的联合概率分布是\(P(Y,Z \mid \theta)\)，那么完全数据的对数似然函数是\(log P(Y,Z \mid \theta)\)。

EM算法通过迭代求\(L(\theta) = log P(Y \mid \theta)\)的极大似然估计。每次迭代包括两步：E步，求期望；M步，求极大化。

EM算法：

输入：观测变量数据Y，隐变量数据Z，联合分布\(P(Y,Z \mid \theta)\)，条件分布\(P(Z \mid Y, \theta)\);

输出：模型参数\(\theta\)

(1) 选择参数的初值\(\theta^{(0)}\)，开始迭代；

(2) E步：记\(\theta^{(i)}\)为第i次迭代参数\(\theta\)的估计值，在第i+1次迭代的E步，计算

\[\begin{align} Q(\theta, \theta^{(i)}) &= E_Z[log P(Y,Z \mid \theta) \mid Y, theta^{(i)}] \\ & \sum_Z log P(Y,Z \mid \theta) P(Z \mid Y, \theta^{(i)}) \\ \end{align}\]

这里，\(P(Z \mid Y, \theta^{(i)})\)是在给定观测数据Y和当前的参数估计\(\theta^{(i)}\)下的隐变量数据Z的条件概率分布；

(3) M步：求使\(Q(\theta, \theta^{(i)})\)极大化的\(\theta\)，确定第i+1次迭代的参数的估计值\(\theta^{(i+1)}\)

\[\theta^{(i+1)} = arg max_\theta Q(\theta, \theta^{(i)})\]

(4) 重复第(2)步和第(3)步，直到收敛。

其中，第(2)步的函数\(Q(\theta,\theta^{(i)})\)是EM算法的核心，成为Q函数。

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

\[Q(\theta, \theta^{(i)}) = E_Z[log P(Y,Z \mid \theta) \mid Y, \theta^{(i)}]\]

下面关于EM算法作几点说明：

步骤(1) 参数的初值可以任意选择，但需注意EM算法对初值是敏感的。

步骤(2) E步求\(Q(\theta,\theta^{(i)})\)。Q函数式中Z是未观测数据，Y是观测数据。[log P(Y,Z \mid \theta) \mid Y, \theta^{(i)}]即类似三硬币模型中E步的\(\mu_j^{(i+1)}\)，而对其求期望，则类似于三硬币模型M步中的，\(\frac{1}{n}\sum_{j=1}^n \mu_j^{(i+1)}\)，我们可知得到期望，即能够算出三硬币模型的新估计值。注意，\(Q(\theta,\theta^{(i)})\)第1个变元表示要极大化的参数，第2个变元表示参数的当前估计值。每次迭代实际在求Q函数及其极大。E-step得名于Expectation(期望), 其实也就是我们通常没必要显示计算出数据补全的概率分布, 而只需要计算期望的显著统计量.

步骤(3) M步求\(Q(\theta, \theta^{(i)})\)的极大化，得到\(\theta^{(i+1)}\)，完成一次迭代\(\theta^{(i)} \to \theta^{(i+1)}\)。后面将正面每次迭代使似然函数增大或达到局部极值。在某些情况下，三硬币模型中，可以选择不更新\(\pi\)，来不断迭代计算新的\(p,q\)。如参考文献(Do C B, Batzoglou S. 2008)中所举的例子。

步骤(4) 给出停止条件，一般是对较小的正数\(\varepsilon_1, \varepsilon_2\)，若满足

\[\|\theta^{(i+1)} - \theta^{(i)} \| < \varepsilon_1 或 \| Q(\theta^{(i+1)},\theta^{(i)}) - Q(\theta^{(i)},\theta^{(i)})\| < \varepsilon_2\]

则停止迭代。

如果我们从算法思想的角度来思考EM算法，我们可以发现我们的算法里已知的是观察数据，未知的是隐含数据和模型参数，在E步，我们所做的事情是固定模型参数的值，优化隐含数据的分布，而在M步，我们所做的事情是固定隐含数据分布，优化模型参数的值。比较下其他的机器学习算法，其实很多算法都有类似的思想。比如SMO算法（支持向量机原理(四)SMO算法原理），坐标轴下降法(Lasso回归算法： 坐标轴下降法与最小角回归法小结), 都使用了类似的思想来求解问题。

参考文献