学习笔记——EM算法
EM算法是一种迭代算法，用于含有隐变量(hidden variable)的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望(expectation)；M步，求极大（Maximization）。
EM算法的引入
给一些观察数据，可以使用极大似然估计法，或贝叶斯估计法估计模型参数。但是当模型含有隐变量时，就不能简单地使用这些方法。有些时候，参数的极大似然估计问题没有解析解，只能通过迭代的方法求解，EM算法就是可以用于求解这个问题的一种迭代算法。
- EM算法
- 输入：观测变量数据Y，隐变量数据Z，联合分布\(P(Y, X | \theta)\)，条件分布\(P(Z|Y, \theta)\)；
- 输出：模型参数\(\theta\)
- 选择参数的初值\(\theta^{(0)}\)，开始迭代；
- E步：记\(\theta^{(0)}\)为第\(i\)次迭代参数\(\theta\)的估计值，在第\(i+1\)次迭代的E步，计算$$Q(\theta, \theta^{(i)}) = E_Z[\log P(Y, Z|\theta) | Y, \theta^{(i)}]$$ $$= \sum_Z \log P(Y, Z | \theta) P (Z|Y, \theta^{(i)})$$ 这里，\(P(Z|Y, \theta^{(i)})\)是在给定观测数据Y和当前的参数估计\(\theta^{(i)}\)下隐变量数据Z的条件概率分布；
- M步：求使\(Q(\theta, \theta^{(i)})\)极大化的\(\theta\)，确定第\(i+1\)次迭代的参数的估计值\(\theta^{(i+1)}\) $$\theta^{(i+1)} = \arg \max_{\theta} Q(\theta, \theta^{(i)})$$
- 重复第(2)步和第(3)步，直到收敛。
其中的函数\(Q(\theta, \theta^{(i)})\)是EM算法的核心。称为Q函数。
-
Q函数：完全数据的对数似然函数\(\log P(Y, Z | \theta)\)关于在给定观测数据Y和当前参数\(\theta^{(i)}\)下对未观测数据Z的条件概率分布\(P(Z|Y, \theta^{i})\)的期望，$$Q(\theta, \theta^{(i)}) = E_Z[\log P(Y, Z|\theta) | Y, \theta^{(i)}]$$ $$= \sum_Z P(Z | Y, \theta^{(i)}) \log P(Y, Z | \theta)$$
-
可以证明，每次迭代使似然函数增大或达到局部极值。
-
EM算法可以用于生成模型的非监督学习，生成模型由联合概率分布\(P(X, Y)\)表示，可以认为非监督学习训练数据是联合概率分布产生的数据，X为观测数据，Y为未观测数据。
-
在应用中，初值的选择变得非常重要。常用的办法是选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。
EM算法在高斯混合模型学习中的应用
- 高斯混合模型是指具有以下形式的概率分布模型：$$P(y | \theta) = \sum_{k = 1}^K \alpha_k \phi(y|\theta_k)$$ 其中，\(\alpha_k\)是系数，\(\alpha_k \geq 0\)，\(\sum_{k=1}^K \alpha_k = 1\)；\(\phi\)是高斯分布密度，\(\theta_k = (\mu_k, \sigma_k^2)\)，$$\phi(y|\theta_k) = \frac{1}{\sqrt{2\pi} \sigma_k} \exp (-\frac{(y - \mu_k)^2}{2 \sigma_k^2})$$称为第k个分模型。
使用EM算法进行高斯混合模型参数估计
假设观测数据\(y_1,y_2,...,y_N\)由高斯混合模型生成，$$P(y | \theta) = \sum_{k = 1}^K \alpha_k \phi(y|\theta_k)$$其中，\(\theta = (\alpha_1, \alpha_2, ..., \alpha_K; \theta_1, \theta_2,...,\theta_K)\)。我们用EM算法估计高斯混合模型的参数\(\theta\)。
-
明确隐变量，写出完全数据的对数似然函数
可以认为观测数据是这样产生的：首先依概率\(\alpha_k\)选择第k个分模型，然后根据这个分模型的概率分布生成一个数据。计算出完全数据的对数似然函数为$$\log P(y, \gamma | \theta) = \sum_{k = 1}^K n_k \log \alpha_k + \sum_{j = 1}^N \gamma_{jk} [\log(\frac{1}{\sqrt{2\pi}}) - \log \sigma_k - \frac{1}{2\sigma^2} (y_i - \mu_k)^2]$$其中，\(\gamma_{jk}\)为1代表第j个观测来自第k个分模型，否则为0，\(n_k = \sum_{j=1}^N \gamma_{jk}\)，\(\sum_{k=1}^K n_k = N\)。
观测数据\(y_j\)及未观测数据\(\gamma_{jk}\)，完全数据是\((y_j,\gamma_{j1},\gamma_{j2},...,\gamma_{jK}, j = 1,2,...,N\) -
EM算法的E步：确定Q函数
计算出Q函数为：$$Q(\theta, \theta^{(i)}) = \sum_{k=1}^K n_k \log \alpha_k + \sum_{k=1}^N \hat{\gamma}{jk} [\log (\frac{1}{\sqrt{2\pi}}) - \log \sigma_k - \frac{1}{2\sigma_k^2} (y_i - \mu_k)^2]$$ 其中，$$\hat{\gamma} = E(\gamma_{jk} | y, \theta) = \frac{\alpha_k \phi (y_i | \theta_k)}{\sum_{k=1}^K \alpha_k \phi (y_i | \theta_k)}$$ $$n_k = \sum_{j = 1}^N E\gamma_{jk}$$ -
确定EM算法中的M步
迭代的M步是求函数\(Q(\theta, \theta^{(i)})\)对\(\theta\)的极大值，即求新一轮迭代的模型参数：$$\theta^{(i+1)} = \arg \max_{\theta} Q(\theta, \theta^{(i)})$$
通过求偏导可以得到计算相应参数的表达式。
- 高斯混合模型参数估计的EM算法
- 输入：观测数据\(y_1, y_2,...,y_N\)，高斯混合模型
- 输出：高斯混合模型参数。
- 取参数的初始值开始迭代
- E步：依据当前模型参数，计算分模型k对观测数据\(y_i\)的响应度$$\hat{\gamma}{jk} = \frac{\alpha_k \phi (y_i | \theta_k)}{\sum^K \alpha_k \phi (y_i | \theta_k)}, j = 1,2,...,N; k = 1,2,...,K$$
- M步：计算新一轮迭代的模型参数：$$\hat{\mu_k} = \frac{\sum_{j=1}^N \hat{\gamma}{jk} y_j}{\sum^N \hat{\gamma}{jk}}, k = 1,2,...,K$$ $$\hat{\sigma}k^2 = \frac{\sum^N \hat{\gamma} (y_j - \mu_k)2}{\sum_{k=1}N \hat{\gamma}{jk}}, k = 1,2,...,K$$ $$\hat{\alpha}k = \frac{\sum^N \hat{\gamma}{jk}}{N}, k = 1,2,...,K$$
- 重复第（2）步和第（3）步，直到收敛。
- EM算法还可以解释为F函数的极大-极大算法，基于这个解释有若干变形与推广，如广义期望极大（GEM）算法。
(注：本文为读书笔记与总结，侧重算法原理，来源为[《统计学习方法》](http://book.douban.com/subject/10590856/)一书第九章)
作者：[rubbninja](http://www.cnblogs.com/rubbninja/) 出处：[http://www.cnblogs.com/rubbninja/](http://www.cnblogs.com/rubbninja/) 关于作者：目前主要研究领域为机器学习与无线定位技术，欢迎讨论与指正！
