EM算法总结

EM算法（Expectation-Maximum）算法。

1，要解决的问题

在一些情况下，我们得到的观察数据有未观察到的隐含数据，此时我们有未知的隐含数据和模型参数，因而无法直接用极大化对数似然函数得到模型分布的参数。怎么办呢？这就是EM算法可以派上用场的地方了。

算法思路

假设未知的隐含数据为z，模型参数为$\theta$。
已知的输入数据为:
$x = (x^{(1)},x^{(2)},...x^{(m)})$
联合分布$p(x,z|\theta)$
条件分布$p(z|x,\theta)$

1. 随即初始画模型参数$\theta$的初值$\theta^0$。

2. for j from 1 to J开始EM算法迭代：
   a. E步：

   $$Q_i(z^{(i)}) = P(z^{(i)}|x^{(i)}, \theta^j)\tag{1}$$
   ((1)需证明，该等式为何成立)
   $$L(\theta, \theta^j) = \sum_{I=1}^m\sum_{z^{(i)}}Q_i(z^{(i)})\log P(x^{(i)},z^{(i)}|\theta)\tag{2}$$
   ((2)需证明，该等式为何成立)
   b. M步：极大化$L(\theta, \theta^j)$,得到$\theta^{j+1}$:
   $$\theta^{j+1} = \arg \max_\theta L(\theta, \theta^j)$$
   c. 如果$\theta^{j+1}已收敛，则算法结束，不然回到步骤a重复迭代$
   ((3)需证明，为何可以收敛)

   数学证明

   证明步骤中，先证明了为什么(2)是最大似然函数。再证明(1)等式成立。最后证明收敛性。
   对于m个样本观察数据$x = (x^{(1)},x^{(2)},...x^{(m)})$，m个未观察到的隐含数据$z = (z^{(1)},z^{(2)},...z^{(m)})$,极大化模型分布的对数似然函数如下：
   

   $$\theta = \arg \max_\theta\sum^m_{i =1}\log P(x^{(i)}|\theta)= \arg\max_\theta\sum^m_{i =1}\log \sum_{z^{(i)}} P(x^{(i)},z^{(i)}|\theta)$$
   再对该等式进行缩放，缩放用到Jensen不等式,通过
   $$\log \sum_j\lambda_j y_j \ge \sum_j\lambda_j \log y_j, \lambda_j \ge 0, \sum_j \lambda_j = 1$$
   可以得出：
   $$\begin{equation} \begin{aligned} \sum^m_{i =1}\log \sum_{z^{(i)}} P(x^{(i)},z^{(i)}|\theta)&=\sum^m_{i =1}\log \sum_{z^{(i)}} Q_i(z^{(i)}) \frac {P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})} \\ &\ge \sum^m_{i =1}\sum_{z^{(i)}} Q_i(z^{(i)}) \log \frac {P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})} \end{aligned} \end{equation}$$
   由于$\sum_z Q_i(z^{(i)}) = 1,$当$\frac {P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})} = c$，等号成立
   而从上面两等式中，我们可以得到：
   $$\begin{equation} \begin{aligned} {Q_i(z^{(i)})} &= \frac{P(x^{(i)},z^{(i)}|\theta)}{c}\\ &= \frac{P(x^{(i)},z^{(i)}|\theta)}{c \times \sum_z Q_i(z^{(i)})}\\ &= \frac{P(x^{(i)},z^{(i)}|\theta)} {\sum_z P(x^{(i)},z^{(i)}|\theta)} \\ &= \frac{P(x^{(i)},z^{(i)}|\theta)} { P(x^{(i)}|\theta)}\\ &= P(z^{(i)}|x^{(i)},\theta) \end{aligned} \end{equation}$$
   所以(1)中证明成立
   所以当等式成立时，$\sum^m_{i =1}\sum_{z^{(i)}} Q_i(z^{(i)}) \log \frac {P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})}$是对数似然函数的下界，求该函数的极大下界。有助于极大化$\sum^m_{i =1}\log \sum_{z^{(i)}} P(x^{(i)},z^{(i)}|\theta)$.所以：
   $$\begin{equation} \begin{aligned} &\arg \max_\theta\sum^m_{i =1}\sum_{z^{(i)}} Q_i(z^{(i)}) \log \frac {P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})} \\ = &\arg \max_\theta\sum^m_{i =1}\sum_{z^{(i)}} Q_i(z^{(i)}) (\log {P(x^{(i)},z^{(i)}|\theta)}-\log {Q_i(z^{(i)})}) \\ = &\arg \max_\theta\sum^m_{i =1}\sum_{z^{(i)}} Q_i(z^{(i)}) \log {P(x^{(i)},z^{(i)}|\theta)} \end{aligned} \end{equation}$$
   所以(2)中证明成立
   (3)要求EM算法收敛时，则我们需要证明我们的对数似然函数的值在迭代的过程中一直在增大。 只需证明$\sum_{i=1}^m \log P(x^{(i)}|\theta^{j+1}) \ge \sum_{i=1}^m \log P(x^{(i)}|\theta^j)$，这里不做证明。详细可以参考EM算法总结
   以上就是EM算法的整个过程，和数学证明。

代码

#-*-coding:utf-8-*-
import numpy as np 
from scipy.stats import binom

def em_single(priors, observations):
    counts = {"A": {"H": 0, "T": 0}, "B": {"H": 0, "T": 0}}
    theta_A = priors[0]
    theta_B = priors[1]

    for observation in  observations:
        len_observation = len(observation)
        num_head = sum(observation)
        num_tail = len_observation - num_head

        contribution_A = binom.pmf(num_head, len_observation, theta_A)
        contribution_B = binom.pmf(num_head, len_observation, theta_B)

        weight_A = contribution_A/(contribution_A + contribution_B)
        weight_B = contribution_B/(contribution_A + contribution_B)

        counts["A"]["H"] += weight_A * num_head
        counts["A"]["T"] += weight_A * num_tail
        counts["B"]["H"] += weight_B * num_head
        counts["B"]["T"] += weight_B * num_tail

    new_theta_A = counts["A"]["H"] / (counts["A"]["H"] + counts["A"]["T"])
    new_theta_B = counts["B"]["H"] / (counts["B"]["H"] + counts["B"]["T"])
    return (new_theta_A, new_theta_B)


def em(observations, prior, tol = 1e-6, itertions = 10000):
    itertion = 0
    while itertion < itertions:
        new_prior = em_single(prior, observations)
        delta_change = np.abs(prior[0] - new_prior[0])
        if delta_change < tol:
            break
        else:
            prior = new_prior
            itertion += 1
    return ([new_prior, itertion])


def main(argv = None):
    observations = np.array([[1,0,0,0,1,1,0,1,0,1],
                        [1,1,1,1,0,1,1,1,0,1],
                        [1,0,1,1,1,1,1,0,1,1],
                        [1,0,1,0,0,0,1,1,0,0],
                        [0,1,1,1,0,1,1,1,0,1]])
    print em(observations, [0.8, 0.2])

if __name__ == "__main__":
    main()