EM算法公式推导和范例详解

在研究kallisto的时候,遇到了EM算法,为了能够理解其代码和自己进行适当的运用,在百度谷歌上进行了一干搜索,最后逐步理解了EM算法在做什么,怎么做,局限性等等。本文将从EM算法的推导讲起,最后用两个范例来详细地演示EM算法是如何运作的。

如何理解EM算法?

通常来讲,一堆由已知分布得到的数据,如果模型中的变量都是可以观测到的,为了求其中的参数$\theta$,可以直接使用极大似然估计。然而有这么一种情况,当模型含有隐含变量的时候,极大似然估计的计算过程就会显得极为复杂,甚至不可解。这时候就需要用到EM算法了,那么,EM算法的思想究竟是什么呢?
看这里,有一篇写的很好的博客————如何感性地理解EM算法

用一种非数学公式推导,而且较为直观的弄懂EM算法思想后,下面,进行数学推导。

EM算法的数学公式推导

公式推导过程

假设,有一个含有隐含变量$Z$的模型,其概率密度函数为$P(y,z|\theta)$。现在,我们希望得到$$\hat{\theta}=\mathop{argmax}_{\theta}logL(Y|\theta)$$
而$logL(Y|\theta)$有,$$ \begin{align} logL(Y|\theta) & = \mathop{\sum}_{j}logP(y_j|\theta) \\ & = \mathop{\sum}_{j}log\mathop{\sum}_{z^{(i)}}P(y_j,z^{(i)}|\theta) \end{align}$$

这一步,由边缘分布$P(y_j|\theta) = \mathop{\sum}_{z^{(i)}}P(y_j,z^{(i)}|\theta)$得到,接下来,有,$$\begin{align} logL(Y|\theta) & = \mathop{\sum}_{j}log\mathop{\sum}_{z^{(i)}}Q(z^{(i)})\frac{P(y_j,z^{(i)}|\theta)}{Q(z^{(i)})}\end{align}$$
这一步,$Q(z^{(i)})$,被认为是$z^{(i)}$的概率密度函数,因此有,$$Q(z^{(i)})\ge0,\mathop{\sum}_{z^{(i)}}Q(z^{(i)})=1$$
或许你会问,$Q(z^{(i)})$怎么来的,李航的《统计学习方法》中,认为Q函数是一个定义;斯坦福大学在网易上的公开课中,这一步中Q函数也是直接写出来;更别说某些多而无实的博客了。根据师兄的讲解,这一步Q函数的推导是为了找到一个关于$Z$随机变量的分布,读者请自行体会一下。

我们往下继续推导,$$logL(Y|\theta) = \mathop{\sum}_{j}logE(\frac{P(y_j,z^{(i)}|\theta)}{Q(z^{(i)})})$$
这一步,来自于$E(g(x))=\int g(x)f(x)dx$,$f(x)$是概率密度函数。继续往下,由于$f(x)=logx$是凹函数[1],根据Jense不等式,必有$$logE(\frac{P(y_j,z^{(i)}|\theta)}{Q(z^{(i)})})\ge E(log\frac{P(y_j,z^{(i)}|\theta)}{Q(z^{(i)})})$$因而,有$$logL(Y|\theta)\ge \mathop{\sum}_{j}E(log\frac{P(y_j,z^{(i)}|\theta)}{Q(z^{(i)})})=\mathop{\sum}_{j}\mathop{\sum}_{z^{(i)}}log\frac{P(y_j,z^{(i)}|\theta)}{Q(z^{(i)})}*Q(z^{(i)})$$
做到这一步,可能就不知道我们在干啥了,其实是这么一回事。假设一个对数似然函数$L(Y|\theta)$,其曲线如图1所示,那么我们初始化一个$\theta^{(0)}$,根据推导出来的上式,对数似然函数$L(Y|\theta^{(0)})$可以在图1上如此表示[2],再求出其极大值下的$\theta$,作为$\theta^{(1)}$,依次迭代,如此以来,可以紧逼下界,最终得到局部最优值[3]

在上述求解局部最优值过程中,拿到$L(Y|\theta^{(i)})$被称为E步,而求解$argmaxlogL(Y|\theta^{(i)})$为M步。

Q函数

Q函数的版本有很多,查阅不同的博客,教材,PPT,网课,可能会有不同的数学表达式,数学符号,但无一例外,Q函数一定和$P(Z|Y|\theta)$这个条件分布有关。下面是斯坦福大学网课中的Q函数,$$Q(z^{(i)})=P(z^{(i)}|y^{(i)},\theta^{(i)})$$

E step and M step

E步:$$Q(z^{(i)})=P(z^{(i)}|y^{(i)},\theta^{(i)})$$
M步:$$\hat{\theta} = \mathop{argmax}_{\theta}\mathop{\sum}_{j}\mathop{\sum}_{z^{(i)}}Q(z^{(i)})log\frac{P(y_j,z^{(i)}|\theta)}{Q(z^{(i)})}$$
网易公开课:斯坦福大学公开课——机器学习之Kmeans算法

EM算法范例1:双硬币,伯努利分布

题目:假设有两枚硬币,A和B。随机选一枚进行抛投,抛投的过程中不知道是A还是B,如此进行10次,最后得到1,1,0,1,0,0,1,0,1,1。估算A和B抛出正面的概率。
求解:

  1. 模型建立
    假设随机选择硬币过程中,选A设为$z=1$,概率为$m$,那么选B的概率是$1-m$,显然,因为无法观测到每一次是选择了A还是B,所以这是一个隐变量,设为$Z$。认为A和B抛出正面的概率分别是$a$和$b$,服从伯努利分布。采用EM算法对$\theta=(m,a,b)$进行估计。
  2. E步
    有,$$\begin{align} Q(z^{(i)}) & =P(z^{(i)}|y^{(i)},\theta^{(i)}) \\ \\ & = {\frac{P(y^{(i)},z^{(i)}|\theta^{(i)})}{P(y^{(i)})}} \end{align}$$
    step1,很明显有,$$P(z)_Z = \begin{cases} m, & z=1 \\ 1-m, & z=0 \\ \end{cases}$$,$$P(y|z)_{Y|Z} = \begin{cases} a^y{(1-a)}^{1-y}, & z=1 \\ b^y{(1-b)}^{1-y}, & z=0 \\ \end{cases}$$
    由条件分布公式,有,$$\begin{align} P(y,z)_{Y,Z} & = P(y|z)_{Y|Z}P(z)_Z \\ \\ & = \begin{cases} ma^y{(1-a)}^{1-y}, & z=1 \\ (1-m)b^y{(1-b)}^{1-y}, & z=0 \\ \end{cases} \end{align}$$
    由边缘分布公式,有,$$\begin{align} P(y)_Y & = \mathop{\sum}_{z}P(y,z)_{Y,Z} \\ & = \mathop{\sum}_{z}P(y|z)_{Y|Z}P(z)_Z \\ & = ma^y{(1-a}^{1-y})+(1-m)b^y{(1-b)}^{1-y} \end{align}$$
    综上,有,$$ Q(z^{(i)}) = \begin{cases} {\frac{ma^y{(1-a)}^{1-y}}{ma^y{(1-a)}^{1-y} + (1-m)b^y{(1-b)}^{1-y}}}, & z=1 \\ {\frac{(1-m)b^y{(1-b)}^{1-y}}{ma^y{(1-a)}^{1-y} + (1-m)b^y{(1-b)}^{1-y}}}, & z=0 \\ \end{cases}$$
  3. M步
    $$\begin{align} \hat{\theta} & = \mathop{argmax}_{\theta}\mathop{\sum}_{j}\mathop{\sum}_{z^{(i)}}Q(z^{(i)})log\frac{P(y_j,z^{(i)}|\theta)}{Q(z^{(i)})} \\ & = \mathop{argmax}_{\theta}\mathop{\sum}_{j}\mathop{\sum}_{z^{(i)}}Q(z^{(i)})log\frac{P(y_j|z=z^{(i)})P(z=z^{(i)})}{Q(z^{(i)})} \end{align}$$
    显然,$z^{(i)}$的取值只有0和1,代入上式,有
    $$\begin{align} \hat{\theta} & = \mathop{argmax}_{\theta}\mathop{\sum}_{j}(Q(1)log\frac{ma^{y_j}(1-a)^{1-y_j}}{Q(1)}+Q(0)log\frac{b^{y_j}(1-b)^{1-y_j}}{Q(0)}) \\ & = \mathop{argmax}_{\theta}\mathop{\sum}_{j}(Q(1)[logm+y_jloga+(1-y_j)log(1-a)-log(Q(1))] \\ & +Q(0)[log(1-m)+y_jlogb+(1-y_j)logb+(1-y_j)log(1-b)-log(Q(0))])\end{align}$$
    对$m$求导,有,$$\begin{align} & \frac{dL}{dm^{(i+1)}} = \mathop{\sum}_{j}\frac{Q(1)}{m} - \mathop{\sum}_{j}\frac{Q(0)}{1-m} = 0 \\ \Rightarrow & m^{(i+1)} = \frac{\mathop{\sum}_{j}^{n}Q(1)}{\mathop{\sum}_{j}^{n}[Q(1)+Q(0)]}= \frac{\mathop{\sum}_{j}^{n}Q(1)}{n} \end{align}$$
    对$a$求导,有$$a^{(i+1)} = \frac{\sum_{j}^{n}Q(1)y_j}{\sum_{j}^{n}Q(1)}$$
    对$b$求导,有$$b^{(i+1)} = \frac{\sum_{j}^{n}Q(0)y_j}{\sum_{j}^{n}Q(0)}$$
  4. 初始化$\theta^{(0)}=(m=0.5,a=0.5,b=0.5)$
  5. 迭代计算
    E步:i=0$$Q(z^{(0)}) = \begin{cases} 0.5, & z=1 \\ 0.5, & z=0 \\ \end{cases}$$
    M步:i=0$$m^{(1)} = 0.5, a^{(1)} = 0.6, b^{(1)} = 0.6$$
    E步:i=1
    M步:i=1
    直到收敛
  6. 收敛判别
    一般情况下,找到两个较小的正数$\xi_1, \xi_2$使得$$|\theta^{(i+1)}-\theta^{(i)}|<\xi_1$$或$$|Q(z^{(i+1)})-Q(z^{(i)})|<\xi_2$$

EM算法范例2:双硬币,二项分布

题目:假设有两组硬币,A和B,每组5个。随机选一组进行抛投,抛投的过程中不知道是A还是B,进行5次,得到以下结果:

coin value sum
U 11010 3T-2F
U 00110 2T-3F
U 10000 1T-4F
U 10011 3T-2F
U 01100 2T-3F

求解:

  1. 模型建立
    假设随即选择一组硬币时,选A设为$z=1$,概率为$\pi$,那么选B的概率是$1-\pi$,显然,因为无法观测到每一次是选择了A还是B,所以这是一个隐变量,设为$Z$。认为每组抛投硬币的结果均服从二项分布$B(n~p)$,其中,A组每次成功的概率为$p$,B组每次成功的概率为$q$。
  2. E步
    有,$$\begin{align} Q(z^{(i)}) & =P(z^{(i)}|y^{(i)},\theta^{(i)}) \\ \\ & = \frac{P(y|z^{(i)})P(z^{(i)})}{\sum_{j}^{n}P(y|z^{(i)})P(z^{(i)})} \\ \\ & = \frac{\pi\begin{pmatrix} n \\ k \end{pmatrix}p^k{(1-p)}^{n-k}}{\pi\begin{pmatrix} n \\ k \end{pmatrix}p^k{(1-p)}^{n-k}+(1-\pi)\begin{pmatrix} n \\ k \end{pmatrix}q^k{(1-q)}^{n-k}}, z=1 \\ or & = \frac{(1-\pi)\begin{pmatrix} n \\ k \end{pmatrix}q^k{(1-q)}^{n-k}}{\pi\begin{pmatrix} n \\ k \end{pmatrix}p^k{(1-p)}^{n-k}+(1-\pi)\begin{pmatrix} n \\ k \end{pmatrix}q^k{(1-q)}^{n-k}}, z=0 \end{align}$$
  3. M步
    $$\begin{align} \hat{\theta} & = \mathop{argmax}_{\theta}\mathop{\sum}_{j}\mathop{\sum}_{z^{(i)}}Q(z^{(i)})log\frac{P(y_j,z^{(i)}|\theta)}{Q(z^{(i)})} \\ & = \mathop{argmax}_{\theta}\mathop{\sum}_{j}\mathop{\sum}_{z^{(i)}}Q(z^{(i)})log\frac{P(y_j|z=z^{(i)})P(z=z^{(i)})}{Q(z^{(i)})} \end{align}$$
    显然,$z^{(i)}$的取值只有0和1,代入上式,有
    $$\begin{align} \hat{\theta} & = \mathop{argmax}_{\theta}\mathop{\sum}_{j}(Q(1)log\frac{\pi\begin{pmatrix} n \\ k_j \end{pmatrix}p^{k_j}{(1-p)}^{n-k_j}}{Q(1)}+Q(0)log\frac{(1-pi)\begin{pmatrix} n \\ k_j \end{pmatrix}q^{k_j}{(1-q)}^{n-k_j}}{Q(0)}) \\ & = \mathop{argmax}_{\theta}\mathop{\sum}_{j}(Q(1)[log\pi+k_jlogp+(n-k_j)log(1-p)+log\begin{pmatrix} n \\ k_j \end{pmatrix}-log(Q(1))] \\ & +Q(0)[log(1-\pi)+k_jlogq+(n-k_j)logq+(n-k_j)log(1-q)+log\begin{pmatrix} n \\ k_j \end{pmatrix} -log(Q(0))])\end{align}$$
    对$m$求导,有,$$\pi^{(i+1)} = \frac{\mathop{\sum}_{j}^{M}Q(1)}{M} $$
    对$a$求导,有$$p^{(i+1)} = \frac{\sum_{j}^{M}Q(1)k_j}{\sum_{j}^{M}Q(1)n}$$
    对$b$求导,有$$q^{(i+1)} = \frac{\sum_{j}^{M}Q(0)k_j}{\sum_{j}^{M}Q(0)n}$$

总结

在上述计算过程中,可以发现M步中log下的$Q(z^{(i)})$会被当作常数消去,因此M步公式可以省略这个log下的$Q(z^{(i)})$。(感谢评论里“牛牛龙”的指出)
不难发现,这里的EM算法,隐含变量$Z$都是服从于伯努利分布,非彼即此;另外一种较为复杂的是隐含变量$Z$服从于多项分布。再难一点,或者更深入的隐含变量的分布,在当前阶段,我还并未遇到。不过我也是表示祈祷,kallisto里最好没有太过复杂的EM应用。即使说思想容易理解了,但具体的步骤实现,收敛判别的划定,都很麻烦或者说是复杂的。
以后如果对EM算法理解更深了,就再开一篇博客做补充,毕竟这个MathJax已经有些超负荷了。


  1. 1.这里很有意思,国内某些教材和国外定义凹凸函数不同,详见维基百科——凸函数
  2. 2.有下界的凹函数,可以求极大值下的参数$\theta^{(i)}$
  3. 3.局部最优值,如图1,没有得到全局最优值
打赏还是得开着的,万一有人打赏呢?