《统计学习方法》 李航著 第九章 EM算法

我是小白一个；本文代码转载地址文末有注释；除代码和部分注释外大部分自己书写，有问题请多指教

模拟课本第一个例子，即用EM算法估计三个硬币模型的参数。

问题引入：有三个硬币：A、B、C。现在先投掷A，根据A的结果进行选择，A为正面选B，A为反面选C。然后用选出的硬币进行投掷。结果[1, 1, 0, 1, 0, 0, 1, 0, 1, 1]。问题：先只能观测到抛硬币结果，不能观测过程的前提下，如何估计三个硬币正面出现的概率，即三个硬币模型的参数。

EM算法解决代码：

import numpy as np

import math

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 # 关于yield强烈推荐这篇： https://developer.ibm.com/zh/articles/os-cn-python-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

data = [1, 1, 0, 1, 0, 0, 1, 0, 1, 1]

# 由于EM算法与初始值有关，这里有两组不同的初始值，以便体会

# 第一组

# em = EM(prob=[0.5, 0.5, 0.5])

# f = em.fit(data)

# 如果用yield

# 需要下边代码

# next(f) # 执行到上面的yield卡主，即只完成了初始化

# next(f) # 继续执行，到下一次yield之前。也就是完成了一次迭代

# 第二组

em = EM(prob=[0.4, 0.6, 0.7])

f2 = em.fit(data)

结果：

没有用yield时：

init prob:0.4, 0.6, 0.7

1/10 pro_a:0.406, pro_b:0.537, pro_c:0.643

2/10 pro_a:0.406, pro_b:0.537, pro_c:0.643

3/10 pro_a:0.406, pro_b:0.537, pro_c:0.643

4/10 pro_a:0.406, pro_b:0.537, pro_c:0.643

5/10 pro_a:0.406, pro_b:0.537, pro_c:0.643

6/10 pro_a:0.406, pro_b:0.537, pro_c:0.643

7/10 pro_a:0.406, pro_b:0.537, pro_c:0.643

8/10 pro_a:0.406, pro_b:0.537, pro_c:0.643

9/10 pro_a:0.406, pro_b:0.537, pro_c:0.643

10/10 pro_a:0.406, pro_b:0.537, pro_c:0.643

用yield的时候一条一条输出上面的语句