一，思考问题

1. L(x; θ)是似然函数，那么θ代表什么？

   θ代表的是概率，参考《图解EM算法投硬币的例子good.pdf》
   寻找具有给定观测值(x值)的最大可能性的θ值，参考《为什么要使用最大似然估计good.docx》
   似然函数实质上是样本的分布律或分布密度。

2. 最大似然估计是干什么的？

   就是估计一个最大最接近的真实情况的概率。

3. 最大似然估计法的一般步骤？

   写似然函数L；取对数；求导数，得驻点，最大值点；作结论。

4. 最大似然估计和EM算法的区别？

   EM算法中需要用到最大似然估计的思想。即已知某个參数能使这个样本出现的概率最大，我们当然不会再去选择其它小概率的样本，所以干脆就把这个參数作为预计的真实值。
   所以M步骤就用到了最大似然去估计Q去优化θ。

5. EM算法的终结条件是什么？

   用迭代法的话，最后的终止条件可以自定义，当迭代到一定次数，或者算法收敛到一定精度则停止。参考《图解EM算法投硬币的例子good.pdf》大概第10次的时候停止了。一开始的概率θ是估计的，然后求出期望来代替实际值。就又可以求概率θ。依次循环的一个过程。

6. Jasen不等式和EM算法有什么关系？

   EM算法需要利用Jasen不等式来计算和推导。

7. 推导的公式与实际工程计算有什么关系？

   推导的目的是证明了了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

https://blog.csdn.net/u011300443/article/details/46763743