EM算法原理和应用

AlgorithmDog

发布于 2018-01-08 16:14:39

1.4K0

发布于 2018-01-08 16:14:39

文章被收录于专栏：AlgorithmDog的专栏

EM算法是带隐变量概率模型的推断算法。今天我们介绍 EM 算法的原理和应用。我们先介绍推导出 EM 算法的一般方法，再介绍另一种 EM 算法推导方法，最后将 EM 算法应用于高斯混合模型。

1.EM 算法的推导

我们让 (\pmb{y}) 表示可见变量，(\pmb{x}) 表示隐变量，(\theta) 表示模型参数。EM算法的目标是得到一个参数(\theta)，使得概率模型 (p(\pmb{y}|\theta)) 似然值极大。

*** QuickLaTeX cannot compile formula:
\begin{eqnarray*}
\theta &=& argmax_{\theta} p(\pmb{y}|\theta) \nonumber \
       &=& argmax_{\theta}lnp(\pmb{y}|\theta) \nonumber \
       &=& argmax_{\theta} ln\sum_{\pmb{x}}p(\pmb{x},\pmb{y}|\theta) \nonumber
\end{eqnarray*}

*** Error message:
Extra alignment tab has been changed to \cr.
leading text:        &=&

在通常情况下, 直接优化 (lnp(\pmb{y}|\theta)) 很困难。我们希望找到 (lnp(\pmb{y}|\theta)) 的一个下界函数，然后极大这个下界。EM 算法处于第 i 步时，参数为 (\theta^i)。

*** QuickLaTeX cannot compile formula:
\begin{eqnarray*}
lnp(\pmb{y}|\theta) &=& ln\sum_{\pmb{x}}p(\pmb{x},\pmb{y}|\theta) \nonumber \
&=& ln\sum_{\pmb{x}} \frac{p(\pmb{x},\pmb{y}|\theta)}{p(\pmb{x}|\pmb{y},\theta^i)}p(\pmb{x}|\pmb{y},\theta^i) \nonumber \
&\ge& \sum_{\pmb{x}} p(\pmb{x}|\pmb{y},\theta^i) lnp(\pmb{x},\pmb{y}|\theta)+H(p(\pmb{x}|\pmb{y},\theta^i)) \quad ln是凹函数 \nonumber \
&\ge& E_{\pmb{x}}[lnp(\pmb{x},\pmb{y}|\theta)|\pmb{y},\theta^i]  \nonumber 
\end{eqnarray*}

*** Error message:
Extra alignment tab has been changed to \cr.
leading text: &=&

定义 Q 函数等于这个下界函数。

(3)

有了 Q 函数，EM 算法交替如下 E 步和 M 步。

E步，计算出 Q 函数; M步，求出

以极大化Q函数,即

。

2.EM 算法的另一种推导

2.1. 交替优化 F 函数

假设隐变量的概率分布为 (\bar{p}(\pmb{x})), 我们定义 F 函数

(4)

其中(p(\pmb{x},\pmb{y}|\theta))是(\theta)的连续函数。F函数是 (lnp(\pmb{y}|\theta)) 的下界函数。我们用交替优化算法 (Alternative Method) 优化 F 函数。交替优化方法适用优化两组参数，先固定一组优化另一组，然后固定刚刚优化的那组优化另一组，交替进行下去。故而 EM 算法交替执行如下 E 步和 M 步。

E步，固定参数

, 求隐变量分布以极大化 F 函数，用拉格朗日法求得隐变量分布如下所示 (求解细节见文末附录)。 (5)

M步，固定隐变量分布，极大化 F 函数求得

(6)

2.2. 利用 F 函数推导 EM 算法性质

EM 算法两种推导方法的结果都是极大化 (Q(\theta,\theta^{i}) )。F 函数让 EM 算法的单调性和一致性比较直观。

直观地推导 EM 算法的单调性。 EM算法得到的参数序列为 (\theta^0,\theta^1,…,\theta^k,..) ，则有 (lnp(y|\theta^i) \le lnp(y|\theta^{i+1}), \forall i)。EM 执行第 i 次迭代的 E 步之后，这是F函数为 (F(\theta^i,\bar{p}{\theta^i}))

*** QuickLaTeX cannot compile formula:
\begin{eqnarray*}
F(\theta^i,\bar{p}</em>{\theta^i}) &=& E_{\pmb{x}}[lnp(\pmb{x},\pmb{y}|\theta)|\theta^i,\pmb{y}]+H(p(\pmb{x}|\theta^i,\pmb{y})) \quad \bar{p}<em>{\theta^i} = p(\pmb{x}|\pmb{y},\theta^i)  \nonumber \
&=&  E</em>{\pmb{x}}[ln\frac{p(\pmb{x},\pmb{y}|\theta^i)}{p(\pmb{x}|\theta^i,\pmb{y})}|\theta^i,\pmb{y}] \nonumber \
&=& E_{\pmb{x}}[lnp(Y|\theta^{i})|\theta^i,\pmb{y}]= lnp(\pmb{y}|\theta^i)
\end{eqnarray*}

*** Error message:
Extra alignment tab has been changed to \cr.
leading text: &=&

这说明我们每次执行完E步，F函数等于目标函数(lnp(y|\theta^i))。又由于(\theta^0,\theta^1,…,\theta^k,..)是EM算法得到的参数序列，F函数随着i增大而增大，则容易证明(lnp(y|\theta^i))随着i增大而增大。如下图所示，黑线是F函数，红线是(lnp(\pmb{y}|\theta))。F函数是(lnp(\pmb{y}|\theta))的下界函数，每次完成M和E步骤，F都上升，并和(lnp(\pmb{y}|\theta))交汇一次。

直观地推导 EM 算法的一致性。F 函数在(\theta^{})和(\bar{p}^{})达到局部最大值，则 (lnp(y|\theta)) 也在 (\theta^{}) 局部最大值。(lnp(y|\theta) = F(\bar{p}{\theta}, \theta)) 对于任意 (\theta) 成立。特别地，对于使 (F(\bar{p}{\theta}, \theta)) 极大的参数 (\theta^),有

(8)

我们需证明(lnp(y|\theta^))也是极大值。假设(lnp(y|\theta^))不是极大值，则存在接近（在任意邻域内）(\theta^)的参数(\theta^{})，使得(F(\bar{p}_{\theta^*}, \theta^*) = lnp(y|\theta^*) < lnp(y|\theta^{}) =F(\bar{p}_{\theta^{}}, \theta^{}))。因为(\bar{p_{\theta}})是(\theta)的连续函数，则(\bar{p}_{\theta^{**}})接近(\bar{p}_{\theta^{*}})。这与F函数在(\theta^{*})和(\bar{p}^{*})达到极大值矛盾。

3. EM 算法求解高斯混合模型

高斯混合模型 (Gauss Mixed Model, GMM) 是概率生成模型，其按多项式分布从 K 个类别中选择一个类别，然后按照该类别对应的高斯分布生成样例。让 (\pmb{\pi}) 表示多项分布的参数，(\pi_{i}) 表示第 i 个项目被选择的概率；让(\pmb{u}) 和 (\pmb{\sigma}) 分别表示高斯分布的期望和标准差, 其中(u_{i}) 和 (\sigma_{i}) 分别表示第 i 个高斯分别的期望和标准差。下图是 k = 2, (\pi_0 = 1/3), (\pi_1 = 2/3), (u_0 = -2), (u_1=0),(\sigma_0 = 0.571), (\sigma_1 = 0.418) 时高斯混合模型的概率密度近似图。

用 EM 算法求解高斯混合模型，可见变量 (\pmb{y}) 是生成的 n 个样例，隐变量 (\pmb{x}) 是被选中的类别 (其中 (x_i) 表示生成第 i 个样例时选中的类别), 参数 (\theta) 包括(\pmb{\pi}),(\pmb{u}) 和 (\pmb{\sigma})。

3.1. E 步推导

EM 算法的 E 步是计算 (Q(\theta,\theta^i))。让(Z = \sum_{\pmb{x}} \prod_{s=1}^{n} \pi_{x_s}^i N(u_{x_s}^i,\sigma_{x_s}^i))，我们有

*** QuickLaTeX cannot compile formula:
\begin{eqnarray*}
&&Q(\theta,\theta^i) \nonumber \
&=& \sum_{\pmb{x}} lnp(\pmb{y}|\pmb{x},\theta) p(\pmb{x}|\pmb{y},\theta^{i}) \nonumber \
&=& \sum_{\pmb{x}} {\sum_{t=1}^{n} ln\pi_{x_t}N(y_t;u_{x_t},\sigma_{x_t})}{\frac{\prod_{s=1}^{n} \pi_{x_s}^i N(y_s;u_{x_s}^i,\sigma_{x_s}^i)}{Z}} \nonumber \
&=& \sum_{t=1}^{n}  \sum_{\pmb{x}} { ln\pi_{x_t}N(y_t;u_{x_t},\sigma_{x_t})\frac{\prod_{s=1}^{n} \pi_{x_s}^i N(y_s;u_{x_s}^i,\sigma_{x_s}^i)}{Z}} \nonumber \
&=& \sum_{t=1}^{n}  \sum_{k=1}^{K} { ln \pi_{k}N(y_t;u_{k},\sigma_{k}) \sum_{\pmb{x}:x_t=k}\frac{\prod_{s=1}^{n} \pi_{x_s}^i N(y_s;u_{x_s}^i,\sigma_{x_s}^i)}{Z}}
 \nonumber \
&=& \sum_{t=1}^{n}  \sum_{k=1}^{K} { ln \pi_{k}N(y_t;u_{k},\sigma_{k}) \frac{ \pi_{k}^i N(y_t;u_{k}^i,\sigma_{k}^i)}{\sum_{k=1}^{K}\pi_{k}^i N(y_t;u_{k}^i,\sigma_{k}^i)}}\nonumber \
&=& \sum_{t=1}^{n}  \sum_{k=1}^{K} { ln \pi_{k}N(y_t;u_{k},\sigma_{k}) w_{t,k}} \nonumber 
\end{eqnarray*}

*** Error message:
Extra alignment tab has been changed to \cr.
leading text: &=&

E 步计算 (w_{t,k})，便得 (Q(\theta,\theta^i))，代码 (代码地址:GitHub) 如下所示。

def E_step(y, gmm):
    n = len(y);
    K = len(gmm.pi);
    w = [[0.0 for j in xrange(K)] for i in xrange(n)];
    
    for i in xrange(n):
        sum1 = 0.0;
        for j in xrange(K):
            sum1 += gmm.pi[j] * normal(gmm.u[j], gmm.d[j], y[i]);
        for j in xrange(K):
            w[i][j]  = gmm.pi[j] * normal(gmm.u[j],gmm.d[j], y[i]);
            w[i][j] /= sum1;
    return w;

3.2. M 步推导

EM 算法的 M 步要极大化 (Q(\theta,\theta^i))。在给定 (w_{t,k})的情况下，我们用拉格朗日方法在 (\sum_{k=1}^{K}\pi_{k} = 1) 约束下极大化 (Q(\theta,\theta^i)) ，求得更新公式。

*** QuickLaTeX cannot compile formula:
\begin{eqnarray*}
\pi_{k} &=& \frac{N_k}{N} \nonumber \
u_{k}   &=& \frac{1}{N_k} \sum_{t=1}^{n} w_{t,k} y_t \nonumber \
\sigma_{k} &=& \frac{1}{N_k} \sum_{t=1}^{n} w_{t,k} (y_t - u_k)^2  \nonumber
\end{eqnarray*}

*** Error message:
Extra alignment tab has been changed to \cr.
leading text: u_{k}   &=&

其中 (N = \sum_{k=1}^{K}\sum_{t=1}^{n} w_{t,k}) 和 (N_k = \sum_{t=1}^{n} w_{t,k} )。 M 步的代码 (代码地址:GitHub)。如下所示。

def M_step(y, gmm, w):
    n = len(y);
    K = len(gmm.pi);
    N_k = [0.0 for i in xrange(K)];
    N   = 0.0;
    for j in xrange(K):
        for i in xrange(n):
            N_k[j] += w[i][j];
        N += N_k[j];

    for j in xrange(K):
        gmm.pi[j] = N_k[j] / N;
        gmm.u[j]  = 0.0;
        gmm.d[j]  = 0.0;
        for i in xrange(n):
            gmm.u[j] += w[i][j] * y[i];
            gmm.d[j] += w[i][j] * math.pow(y[i] - gmm.u[j],2)
        gmm.u[j] /= N_k[j];
        gmm.d[j] /= N_k[j];
        gmm.d[j]  = math.sqrt(gmm.d[j]);

    return gmm;

3.2. 一个小实验

下面是我们做一个小实验，实验代码已经放到 GitHub 上了。下面是实验结果，纵坐标是高斯混合模型的对数似然值的近似值，横坐标是 EM 算法的迭代次数。可以看到，随着 EM 算法的推进，对数似然值在上升。

附录

证明公式(\bar{p}{\theta^i}(\pmb{x})=p(\pmb{x}|\theta^i,\pmb{y}))。对固定的(\theta)，极大化F函数是有约束的优化问题：

*** QuickLaTeX cannot compile formula:
\begin{eqnarray*}
&max& \quad  E</em>{\bar{p}}[lnp(\pmb{x},\pmb{y}|\theta)]+H(\bar{p}(\pmb{x})) \nonumber \
&sub& \quad \sum_{\pmb{x}} \bar{p}(\pmb{x}) = 1 \nonumber \
\end{eqnarray*}

*** Error message:
Extra alignment tab has been changed to \cr.
leading text: &sub&

为此，引入拉格朗日系数 (\lambda)，拉格朗日函数如下

(12)

求拉格朗日函数对(\bar{p})的偏导数：

(13)

得(\bar{p}(\pmb{x}) = \frac{p(\pmb{x},\pmb{y}|\theta)}{\exp(1+\lambda)}), 又(\sum_{\pmb{x}} \bar{p}(\pmb{x}) = 1), 所以我们有

(14)

即公式成立。

本文参与腾讯云自媒体同步曝光计划，分享自作者个人站点/博客。

原始发表：2014年9月21日，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

本文分享自作者个人站点/博客前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

编程算法

登录后参与评论

0 条评论

热度