复现经典：《统计学习方法》第 9 章 EM 算法及其推广
本文是李航老师的《统计学习方法》[1]一书的代码复现。
作者：黄海广[2]
备注：代码都可以在github[3]中下载。
我将陆续将代码发布在公众号“机器学习初学者”，敬请关注。
代码目录
第 1 章 统计学习方法概论 第 2 章 感知机 第 3 章 k 近邻法 第 4 章 朴素贝叶斯 第 5 章 决策树 第 6 章 逻辑斯谛回归 第 7 章 支持向量机 第 8 章 提升方法 第 9 章 EM 算法及其推广 第 10 章 隐马尔可夫模型 第 11 章 条件随机场 第 12 章 监督学习方法总结 代码参考：wzyonggege[4],WenDesi[5],火烫火烫的[6]
第 9 章 EM 算法及其推广
Expectation Maximization algorithm
Maximum likehood function
likehood & maximum likehood[7]
1．EM 算法是含有隐变量的概率模型极大似然估计或极大后验概率估计的迭代算法。含有隐变量的概率模型的数据表示为 )。这里，是观测变量的数据，是隐变量的数据， 是模型参数。EM 算法通过迭代求解观测数据的对数似然函数的极大化，实现极大似然估计。每次迭代包括两步：
步，求期望，即求 )关于)的期望：
称为函数，这里是参数的现估计值；
步，求极大，即极大化函数得到参数的新估计值：
在构建具体的 EM 算法时，重要的是定义函数。每次迭代中，EM 算法通过极大化函数来增大对数似然函数。
2．EM 算法在每次迭代后均提高观测数据的似然函数值，即
在一般条件下 EM 算法是收敛的，但不能保证收敛到全局最优。
3．EM 算法应用极其广泛，主要应用于含有隐变量的概率模型的学习。高斯混合模型的参数估计是 EM 算法的一个重要应用，下一章将要介绍的隐马尔可夫模型的非监督学习也是 EM 算法的一个重要应用。
4．EM 算法还可以解释为函数的极大-极大算法。EM 算法有许多变形，如 GEM 算法。GEM 算法的特点是每次迭代增加函数值（并不一定是极大化函数），从而增加似然函数值。
在统计学中，似然函数（likelihood function，通常简写为 likelihood，似然）是一个非常重要的内容，在非正式场合似然和概率（Probability）几乎是一对同义词，但是在统计学中似然和概率却是两个不同的概念。概率是在特定环境下某件事情发生的可能性，也就是结果没有产生之前依据环境所对应的参数来预测某件事情发生的可能性，比如抛硬币，抛之前我们不知道最后是哪一面朝上，但是根据硬币的性质我们可以推测任何一面朝上的可能性均为 50%，这个概率只有在抛硬币之前才是有意义的，抛完硬币后的结果便是确定的；而似然刚好相反，是在确定的结果下去推测产生这个结果的可能环境（参数），还是抛硬币的例子，假设我们随机抛掷一枚硬币 1,000 次，结果 500 次人头朝上，500 次数字朝上（实际情况一般不会这么理想，这里只是举个例子），我们很容易判断这是一枚标准的硬币，两面朝上的概率均为 50%，这个过程就是我们运用出现的结果来判断这个事情本身的性质（参数），也就是似然。
E step:
import numpy as npimport mathpro_A, pro_B, por_C = 0.5, 0.5, 0.5
def pmf(i, pro_A, pro_B, por_C): pro_1 = pro_A * math.pow(pro_B, data[i]) * math.pow( (1 - pro_B), 1 - data[i]) pro_2 = pro_A * math.pow(pro_C, data[i]) * math.pow( (1 - pro_C), 1 - data[i]) return pro_1 / (pro_1 + pro_2)M step:
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_Cdata=[1,1,0,1,0,0,1,0,1,1]em = EM(prob=[0.5, 0.5, 0.5])f = em.fit(data)next(f)init prob:0.5, 0.5, 0.5
# 第一次迭代f.send(1)1/10 pro_a:0.500, pro_b:0.600, pro_c:0.600
# 第二次f.send(2)
2/10 pro_a:0.500, pro_b:0.600, pro_c:0.600
em = EM(prob=[0.4, 0.6, 0.7])f2 = em.fit(data)next(f2)
init prob:0.4, 0.6, 0.7
f2.send(1)
1/10 pro_a:0.406, pro_b:0.537, pro_c:0.643
f2.send(2)
2/10 pro_a:0.406, pro_b:0.537, pro_c:0.643
参考资料
[2] 黄海广: https://github.com/fengdu78
[3] github: https://github.com/fengdu78/lihang-code
[4] wzyonggege: https://github.com/wzyonggege/statistical-learning-method
[5] WenDesi: https://github.com/WenDesi/lihang_book_algorithm
[6] 火烫火烫的: https://blog.csdn.net/tudaodiaozhale
往期精彩回顾