统计机器学习笔记——EM算法及其应用(1)

EM算法的适用场景：

EM算法用于估计含有隐变量的概率模型参数的极大似然估计，或者极大后验概率估计。当概率模型既含有观测值，又含有隐变量或潜在变量时，就可以使用EM算法来求解概率模型的参数。当概率模型只含有观测值时，直接使用极大似然估计法，或者贝叶斯估计法估计模型参数就可以了。

EM算法的入门简单例子：

已知有三枚硬币A，B，C，假设抛掷A，B，C出现正面的概率分别为$\pi$，$p$，$q$。
单次实验的过程是:

1. 首先抛掷硬币A，如果A出现正面选择硬币B，否则，选择硬币C。
2. 抛掷所选择的硬币，正面输出1，反面输出0。

重复上述单词实验n次，需要估计抛掷硬币A，B，C出现正面的概率$\pi$，$p$，$q$。其中每次实验步骤1的抛掷结果不可见，可见的是所挑选硬币的抛掷结果。

解：

$P(y|\theta)=\sum\limits_{z}{P(y,z|\theta)}=\sum\limits_{z}{P(z|\theta)P(y|z,\theta)}$

注：

1. $\theta$ 表示的是整个模型的参数，也就是我们学习的目标。
2. $y$是二元随机变量(取值为0或者1,也可以对应于C或者B)，为观测随机变量。
3. $z$是二元随机变量(取值为0或者1,也可以对应于C或者B), 为隐变量（不可观测）
4. 根据全概率公式可得$P(y|\theta)=\sum\limits_{z}{P(y,z|\theta)}$
5. 而求和中的每一项根据条件概率公式得$P(y,z|\theta)=P(z|\theta)P(y|z,\theta)$
6. $\theta=(\pi, p, q)$

$P(y|\theta)=\pi p^y (1-p)^{1-y} + (1-\pi) q^y (1-q)^{1-y}$

也可做如下理解：

1. $P(y|\theta)=P(y,B|\theta) + P(y,C|\theta)=P(C|\theta)P(y|C,\theta) + P(B|\theta)P(y|B,\theta)$
2. $P(y|\theta)=(1-\pi)P(y|C,\theta) + \pi P(y|B,\theta)$
3. $P(y|\theta)=(1-\pi)q^{y}(1-q)^{1-y} + \pi p^{y}(1-q)^{1-y}$

设$Z=(Z_1,Z_2,\cdots, Z_n)$为n次实验的隐状态（不可观测）序列，$Y=(Y_1,Y_2,\cdots, Y_n)$为n次实验的观测序列

那么每次实验对应观测值的概率为：$P(Y_i|\theta)=\sum\limits_{Z}{P(Z|\theta)P(Y|Z,\theta)}$

这个观测序列的概率为$P(Y|\theta)=\prod\limits_{i=1}^{n}{P(Y_i|\theta)}=\prod\limits_{i=1}^{n}{\{\pi p^y (1-p)^{1-y} + (1-\pi) q^y (1-q)^{1-y}\}}$

那么我们的学习目标就是让$P(Y|\theta)$出现的概率近可能的大，即$\hat{\theta}=\arg\max\limits_{\theta}{\log{P(Y|\theta)}}$

这个问题没有解析解（未知量的个数大于样本数，多出来一些隐变量），只能通过迭代的方法求解。EM算法就是用来解决这类问题的。

求解步骤如下：

1. 初始化参数$\theta ^{(0)}$

2. 然后通过下面的步骤计算参数的估计值,直至收敛:

2.1 E步骤（对隐变量进行估计,是在每个观测值上都要计算的）：计算在参数$\pi ^ {(i)}$, $p ^ {(i)}$, $q ^ {(i)}$下观测数据$y_j$来自投掷硬币B的概率：

$P(B|y_j,\theta)=\frac{P(B,y_j|\theta)}{P(y_j|\theta)}$

$P(y_j|\theta)=P(y_j,B|theta) + P(y_j,C|theta)$

$P(y_j|\theta)=\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)}$

$P(B,y_j|\theta)=P(B|\theta)P(y_j|\theta)$

$P(B,y_j|\theta)=\pi ^ {(i)}(p ^ {(i)})^{y_j}(1-p^{(i)})^{(1-y_j)}$

$P(B|y_j,\theta)=\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)}}=\mu ^{(i+1)} _ j$

2.2 M步骤（是在每个观测值计算E后进行的）：更新模型的权重参数

$P(B_j|\theta)=\pi _ j=\frac{1}{n}\sum\limits_{j=1}^{n}{P(B_j|y_j,\theta)}$

$\pi ^ {(i+1)}=\frac{1}{n}\sum\limits_{j=1}^{n}P(B_j|y_j,\theta)^{(i+1)}=\frac{1}{n}\sum\limits_{j=1}^{n}{\mu ^{(i+1)}_{j}}$

$P(y_j=1|B_j,\theta)=p_j=\frac{P(y_j=1,B_j|\theta)}{P(B_j|\theta)}$

$p ^ {(i+1)}=\frac{\sum\limits_{j=1}^{n}{\mu ^{(i+1)}_j y_j}}{\sum\limits_{j=1}^{n}{\mu ^{(i+1)}_{j}}}$

$P(y_j=1|C_j,\theta)=q_j=\frac{P(y_j=1,C_j|\theta)}{P(C_j|\theta)}=\frac{P(y_j=1,C_j|\theta)}{P(C_j|\theta)}$

$p ^ {(i+1)}=\frac{\sum\limits_{j=1}^{n}{(1-\mu ^{(i+1)}_j) y_j}}{\sum\limits_{j=1}^{n}{(1-\mu ^{(i+1)}_{j})}}$

使用具体数值进行运算：

1. 设初值为：$\pi^{(0)}=0.5$，$p^{(0)}=0.5$，$q^{(0)}=0.5$，观测序列为1，1，0，1，0，0，1，0，1，1
2. $\mu^{(1)}=0.5$
3. $\pi^{(1)}=0.5$，$p^{(1)}=0.6$，$q^{(1)}=0.6$
4. 继续迭代，得$\pi^{(2)}=0.5$，$p^{(2)}=0.6$，$q^{(2)}=0.6$
5. 于是最终的模型参数$\theta$的极大似然估计：$\hat{\pi}=0.5$，$\hat{p}=0.6$，$\hat{q}=0.6$

当换一组初始化权重参数：

1. 设初值为：$\pi^{(0)}=0.4$，$p^{(0)}=0.6$，$q^{(0)}=0.7$，观测序列为1，1，0，1，0，0，1，0，1，1
2. 于是最终的模型参数$\theta$的极大似然估计：$\hat{\pi}=0.4064$，$\hat{p}=0.5368$，$\hat{q}=0.6432$

对应的程序代码为：

import numpy as np


def generate_observe_sequence(n):
    return (np.random.random(size=n)> 0.35).astype(np.int)

def Estep(observe_list, theta):

    def sample_mu(y):
        up_1 = theta[0] * np.power(theta[1], y) * np.power((1-theta[1]),(1-y))
        up_2 = (1-theta[0]) * np.power(theta[2], y) * np.power((1-theta[2]),(1-y))
        return up_1/(up_1 + up_2)

    return [sample_mu(y) for y in observe_list]

def MStep(observe_list, mus):
    p = [0.0, 0.0, 0.0]
    p[0] = sum(mus)/len(mus)
    p[1] = sum([mus[i] * observe_list[i] for i in range(len(observe_list))])/sum(mus)
    p[2] = sum([(1-mus[i]) * observe_list[i] for i in range(len(observe_list))])/sum([1-mu for mu in mus])
    return p[:]

if __name__ == "__main__":
    records = []
    theta = [0.4, 0.6, 0.7]
    m = 1e-7
    records.append(theta)
    observe_list = [1,1,0,1,0,0,1,0,1,1]
    #observe_list = generate_observe_sequence(5)
    print theta
    while True:
        mus = Estep(observe_list, theta)
        new_theta = MStep(observe_list, mus)
        print new_theta
        records.append(new_theta)
        err = 0
        for old, new in zip(theta, new_theta):
            err += np.abs(old-new)
        print err
        if err < m:
            break
        theta = new_theta[:]
    print "###########################"
    for record in records:
        print record

通过上面的例子可以发现，EM算法受初值影响明显。

参考书目：《统计机器学习》李航