开启掘金成长之旅！这是我参与「掘金日新计划 · 2 月更文挑战」的第 32 天，点击查看活动详情
（本文是第35篇活动文章）
6.2.1 示例1：三硬币模型
Step1：观测变量是掷硬币的结果，隐变量是使用哪一枚硬币，而模型参数是硬币的出现正面的概率。
重复以上计算， 直到对数似然函数值不再有明显的变化为止。
完整算法步骤：
代码实现：
import numpy as np
import math
def em_single(theta,Y):
"""每一轮EM算法"""
u_group=[] #Y_i：丢的是第一枚硬币的概率
pi=theta[0]
p=theta[1]
q=theta[2]
#E：通过计算E得到隐变量
for y in Y:
u=(pi*math.pow(p,y)*math.pow(1-p,1-y))/(pi*math.pow(p,y)*math.pow(1-p,1-y)+(1-pi)*math.pow(q,y)*math.pow(1-q,1-y))
u_group.append(u)
#M：计算使Q极大化的参数
new_pi=sum(u_group)/len(Y)
new_p=sum(u_group[i]*Y[i] for i in range(len(Y)))/sum(u_group)
new_q=sum((1-u_group[i])*Y[i] for i in range(len(Y)))/sum(1-u_group[i] for i in range(len(Y)))
return [new_pi,new_p,new_q]
def threecoins(theta,Y):
j=0 #步数
while j<1000:
new_theta=em_single(theta,Y)
if new_theta==theta: #收敛（这里设置一个阈值可能会更好）
break
else:
theta = new_theta
j=j+1
print(new_theta)
return new_theta
if __name__=="__main__":
Y = [1, 1, 0, 1, 0, 0, 1, 0, 1, 1] #其实这里写成X可能会更好
theta=[0.4,0.6,0.7] #初始参数
result=threecoins(theta,Y)
print("模型参数的极大似然估计是",result)
