一，思考问题
-
L(x; θ)是似然函数，那么θ代表什么？
θ代表的是概率，参考《图解EM算法投硬币的例子good.pdf》
寻找具有给定观测值(x值)的最大可能性的θ值，参考《为什么要使用最大似然估计good.docx》
似然函数实质上是样本的分布律或分布密度。 -
最大似然估计是干什么的？
就是估计一个最大最接近的真实情况的概率。
-
最大似然估计法的一般步骤？
写似然函数L；取对数；求导数，得驻点，最大值点；作结论。
-
最大似然估计和EM算法的区别？
EM算法中需要用到最大似然估计的思想。即已知某个參数能使这个样本出现的概率最大，我们当然不会再去选择其它小概率的样本，所以干脆就把这个參数作为预计的真实值。
所以M步骤就用到了最大似然去估计Q去优化θ。 -
EM算法的终结条件是什么？
用迭代法的话，最后的终止条件可以自定义，当迭代到一定次数，或者算法收敛到一定精度则停止。参考《图解EM算法投硬币的例子good.pdf》大概第10次的时候停止了。一开始的概率θ是估计的，然后求出期望来代替实际值。就又可以求概率θ。依次循环的一个过程。
-
Jasen不等式和EM算法有什么关系？
EM算法需要利用Jasen不等式来计算和推导。
-
推导的公式与实际工程计算有什么关系？
推导的目的是证明了了l(θ)会单调增加。一种收敛方法是l（θ）不再变化，还有一种就是变化幅度很小。EM可以看作是J的坐标上升法，E步固定θ，优化Q，M步固定Q优化θ。
直白的说，就是证明每次调参，趋势是向极大似然值走的，并且结果收敛。
二，Python动手实验
import numpy
import scipy.stats
#硬币投掷结果
observations = numpy.array([[1,0,0,0,1,1,0,1,0,1],
[1,1,1,1,0,1,1,1,0,1],
[1,0,1,1,1,1,1,0,1,1],
[1,0,1,0,0,0,1,1,0,0],
[0,1,1,1,0,1,1,1,0,1]])
def em_single(priors,observations):
"""
EM算法的单次迭代
Arguments
------------
priors:[theta_A,theta_B]
observation:[m X n matrix]
Returns
---------------
new_priors:[new_theta_A,new_theta_B]
:param priors:
:param observations:
:return:
"""
counts = {'A': {'H': 0, 'T': 0}, 'B': {'H': 0, 'T': 0}}
theta_A = priors[0]
theta_B = priors[1]
#E step
for observation in observations:
len_observation = len(observation)
num_heads = observation.sum()
num_tails = len_observation-num_heads
#二项分布求解公式
#第一组的10个数据中，包括5正5反，对A投出5正5反的概率为
contribution_A = scipy.stats.binom.pmf(num_heads,len_observation,theta_A)
# 第一组的10个数据中，包括5正5反，对B投出5正5反的概率为
contribution_B = scipy.stats.binom.pmf(num_heads,len_observation,theta_B)
#print(contribution_A,contribution_B)
#第一组实验选择的硬币是来自A的概率为
weight_A = contribution_A / (contribution_A + contribution_B)
#第一组实验选择的硬币是来自B的概率为
weight_B = contribution_B / (contribution_A + contribution_B)
#print(weight_A,weight_B)
#更新在当前参数下A，B硬币产生的正反面次数
counts['A']['H'] += weight_A * num_heads
counts['A']['T'] += weight_A * num_tails
counts['B']['H'] += weight_B * num_heads
counts['B']['T'] += weight_B * num_tails
# M step
new_theta_A = counts['A']['H'] / (counts['A']['H'] + counts['A']['T'])
new_theta_B = counts['B']['H'] / (counts['B']['H'] + counts['B']['T'])
print(new_theta_A,new_theta_B)
return [new_theta_A,new_theta_B]
def em(observations,prior,tol = 1e-6,iterations=10000):
"""
EM算法
：param observations :观测数据
：param prior：模型初值
：param tol：迭代结束阈值
：param iterations：最大迭代次数
：return：局部最优的模型参数
"""
iteration = 0;
while iteration < iterations:
new_prior = em_single(prior,observations)
delta_change = numpy.abs(prior[0]-new_prior[0])
if delta_change < tol:
break
else:
prior = new_prior
iteration +=1
#print(new_prior,iteration)
return [new_prior,iteration]
print("start")
print (em(observations,[0.6,0.5]))
print("end")
三，参考
https://www.cnblogs.com/zfyouxi/p/4297500.html
http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html
