统计学习方法中算法实现（基于Python） --- EM算法

1．EM算法是含有隐变量的概率模型极大似然估计或极大后验概率估计的迭代算法。含有隐变量的概率模型的数据表示为$\theta$ 。这里，$Y$是观测变量的数据，$Z$是隐变量的数据，$\theta$ 是模型参数。EM算法通过迭代求解观测数据的对数似然函数$L(\theta )=\log P(\mathrm{Y}|\theta )$的极大化，实现极大似然估计。每次迭代包括两步：

$E$步，求期望，即求$logP\left(Z|Y,\theta \right)$ )关于$P\left(Z|Y,{\theta }^{(i)}\right)$)的期望：

$$Q\left(\theta ,{\theta }^{(i)}\right)=\sum _Z\log P(Y,Z|\theta )P\left(Z|Y,{\theta }^{(i)}\right)$$
称为$Q$函数，这里${\theta }^{(i)}$是参数的现估计值；

$M$步，求极大，即极大化$Q$函数得到参数的新估计值：

$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q\left(\theta ,{\theta }^{(i)}\right)$$
在构建具体的EM算法时，重要的是定义$Q$函数。每次迭代中，EM算法通过极大化$Q$函数来增大对数似然函数$L(\theta )$。

2．EM算法在每次迭代后均提高观测数据的似然函数值，即

$$P\left(Y|{\theta }^{(i+1)}\right)⩾P\left(Y|{\theta }^{(i)}\right)$$
在一般条件下EM算法是收敛的，但不能保证收敛到全局最优。

3．EM算法应用极其广泛，主要应用于含有隐变量的概率模型的学习。高斯混合模型的参数估计是EM算法的一个重要应用，下一章将要介绍的隐马尔可夫模型的非监督学习也是EM算法的一个重要应用。

4．EM算法还可以解释为$F$函数的极大-极大算法。EM算法有许多变形，如GEM算法。GEM算法的特点是每次迭代增加$F$函数值（并不一定是极大化$F$函数），从而增加似然函数值。

在统计学中，似然函数（likelihood function，通常简写为likelihood，似然）是一个非常重要的内容，在非正式场合似然和概率（Probability）几乎是一对同义词，但是在统计学中似然和概率却是两个不同的概念。概率是在特定环境下某件事情发生的可能性，也就是结果没有产生之前依据环境所对应的参数来预测某件事情发生的可能性，比如抛硬币，抛之前我们不知道最后是哪一面朝上，但是根据硬币的性质我们可以推测任何一面朝上的可能性均为50%，这个概率只有在抛硬币之前才是有意义的，抛完硬币后的结果便是确定的；而似然刚好相反，是在确定的结果下去推测产生这个结果的可能环境（参数），还是抛硬币的例子，假设我们随机抛掷一枚硬币1,000次，结果500次人头朝上，500次数字朝上（实际情况一般不会这么理想，这里只是举个例子），我们很容易判断这是一枚标准的硬币，两面朝上的概率均为50%，这个过程就是我们运用出现的结果来判断这个事情本身的性质（参数），也就是似然。

$$P(Y|\theta )=\prod [\pi p^{y_i}(1-p{)}^{1-y_i}+(1-\pi )q^{y_i}(1-q{)}^{1-y_i}]$$
E step:
$${\mu }^{i+1}=\frac{\pi (p^i{)}^{y_i}(1-(p^i){)}^{1-y_i}}{\pi (p^i{)}^{y_i}(1-(p^i){)}^{1-y_i}+(1-\pi )(q^i{)}^{y_i}(1-(q^i){)}^{1-y_i}}$$

M step:
$${\pi }^{i+1}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{i+1}$$$$p^{i+1}=\frac{{\sum }_{j=1}^n{\mu }_j^{i+1}{y}_i}{{\sum }_{j=1}^n{\mu }_j^{i+1}}$$$$q^{i+1}=\frac{{\sum }_{j=1}^n(1-{\mu }_j^{i+1}{y}_i)}{{\sum }_{j=1}^n(1-{\mu }_j^{i+1})}$$

class EM:
    def __init__(self, prob):
        self.pro_A, self.pro_B, self.pro_C = prob

    # e_step
    def pmf(self, i):
        pro_1 = self.pro_A * math.pow(self.pro_B, data[i]) * math.pow(
            (1 - self.pro_B), 1 - data[i])
        pro_2 = (1 - self.pro_A) * math.pow(self.pro_C, data[i]) * math.pow(
            (1 - self.pro_C), 1 - data[i])
        return pro_1 / (pro_1 + pro_2)

    # m_step
    def fit(self, data):
        count = len(data)
        print('init prob:{}, {}, {}'.format(self.pro_A, self.pro_B,
                                            self.pro_C))
        for d in range(count):
            _ = yield
            _pmf = [self.pmf(k) for k in range(count)]
            pro_A = 1 / count * sum(_pmf)
            pro_B = sum([_pmf[k] * data[k] for k in range(count)]) / sum(
                [_pmf[k] for k in range(count)])
            pro_C = sum([(1 - _pmf[k]) * data[k]
                         for k in range(count)]) / sum([(1 - _pmf[k])
                                                        for k in range(count)])
            print('{}/{}  pro_a:{:.3f}, pro_b:{:.3f}, pro_c:{:.3f}'.format(
                d + 1, count, pro_A, pro_B, pro_C))
            self.pro_A = pro_A
            self.pro_B = pro_B
            self.pro_C = pro_C