em算法系列1：python实现三硬币模型

三硬币模型

假设有3枚硬币，分别记做A，B，C。这些硬币正面出现的概率分别是π,p和q。进行如下掷硬币实验：先掷硬币A，根据其结果选出硬币B或C，正面选B，反面选硬币C；然后投掷选重中的硬币，出现正面记作1，反面记作0；独立地重复n次（n=10)，结果为1111110000
我们只能观察投掷硬币的结果，而不知其过程，估计这三个参数π,p和q。

EM算法

可以看到投掷硬币时到底选择了B或者C是未知的。我们设隐藏变量Z 来指示来自于哪个硬币， Z = { z 1 , z 2 , … , z n } Z=\{z_1,z_2,…,z_n\} Z={z1​,z2​,…,zn​}，令 θ = { π , p , q } θ=\{π,p,q\} θ={π,p,q}，观察数据$X=x_1,x_2,\dots ,x_n$。

写出生成一个硬币时的概率：
$$\begin{array}{rl}P(x|\theta )& =\sum zP(x,z|\theta )=\sum zP(z|\pi )P(x|z,\theta )\\ & =\pi p^x(1-p{)}^{1-x}+(1-\pi )q^x(1-q{)}^{1-x}\end{array}$$
有了一个硬币的概率，我们就可以写出所有观察数据的log似然函数：
$$\begin{array}{r}L(\theta |X)=logP(X|\theta )=\sum _{j=1}^{n}log[\pi p^{x_j}(1-p{)}^{1-x_j}+(1-\pi )q^{x_j}(1-q{)}^{1-x_j}]\end{array}$$
然后求似然函数最大值
$${\theta }^{\ast }=argmaxL(\theta |X)$$
其中$L(\theta |X)=logP(X|\theta )=log{\sum }_{Z}P(X,Z|\theta )$。因为$log$里面带着加和所以这个极大似然是求不出解析解的。 如果我们知道隐藏变量Z的话，求以下的似然会容易很多:
L(θ|X,Z)=logP(X,Z|θ)
但是隐藏变量Z的值谁也不知道，所以我们转用EM算法求它的后验概率期望$E_{Z|X,\theta }(L(\theta |X,Z))$的最大值。

• E步：：假设当前模型的参数为$\pi ,p,q$时，隐含变量来自于硬币B的后验概率
  $$\begin{array}{r}\mu =P(Z|X,\theta )=\frac{\pi p^{x_i}(1-p{)}^{1-x_i}}{\pi p^{x_i}(1-p{)}^{1-x_i}+(1-\pi )q^{x_i}(1-q{)}^{1-x_i}}\end{array}$$
  那么隐含变量来自于硬币C的后验概率自然为$1-\mu$。

def cal_u(pi1,p1,q1,xi):
    return pi1*math.pow(p1,xi) * math.pow(1-p1,1-xi)/\
           float(pi1* math.pow(p1,xi) * math.pow(1-p1,1-xi)+
                 (1-pi1) * math.pow(q1,xi) * math.pow(1-q1,1-xi))

• M步： 先写出似然关于后验概率的期望，它是似然期望下界函数的最大值，
  $$E_{Z|X,\theta }(L(\theta |X,Z))=\sum _{j=1}^n[\mu log\frac{\pi p^{x_i}(1-p{)}^{1-x_i}}{\mu }+(1-\mu )\frac{log(1-\pi )q^{x_i}(1-q{)}^{1-x_i}}{1-\mu }]$$
  要注意这里把$\mu$看做固定值。然后我们分别求偏导，获得参数$\pi ,p,q$的新估计值
  $$\begin{array}{rl}\pi & =\frac{1}{n}\sum _{j=1}^n{\mu }_j\\ p& =\frac{{\sum }_{j=1}^n{\mu }_j{x}_j}{{\sum }_{j=1}^n{\mu }_j}\\ q& =\frac{{\sum }_{j=1}^n(1-{\mu }_j)xj}{{\sum }_{j=1}^n(1-{\mu }_j)}\end{array}$$

def cal_u(pi, p, q, xi):
    """
      u值计算
    :param pi: 下一次迭代开始的 pi
    :param p:  下一次迭代开始的 p
    :param q:  下一次迭代开始的 q
    :param xi: 观察数据第i个值，从0开始
    :return: 
    """
    return pi * math.pow(p, xi) * math.pow(1 - p, 1 - xi) / \
           float(pi * math.pow(p, xi) * math.pow(1 - p, 1 - xi) +
                 (1 - pi) * math.pow(q, xi) * math.pow(1 - q, 1 - xi))

def e_step(pi,p,q,x):
    """
        e步计算
    :param pi: 下一次迭代开始的 pi
    :param p:  下一次迭代开始的 p
    :param q:  下一次迭代开始的 q
    :param x: 观察数据
    :return:
    """
    return [cal_u(pi,p,q,xi) for xi in x]

算法首先选取参数初始值$\theta (0)=\pi (0),p(0),q(0)$，然后迭代到收敛为止

python实现

# -*- coding: utf-8 -*-
import math

def cal_u(pi, p, q, xi):
    """
      u值计算
    :param pi: 下一次迭代开始的 pi
    :param p:  下一次迭代开始的 p
    :param q:  下一次迭代开始的 q
    :param xi: 观察数据第i个值，从0开始
    :return: 
    """
    return pi * math.pow(p, xi) * math.pow(1 - p, 1 - xi) / \
           float(pi * math.pow(p, xi) * math.pow(1 - p, 1 - xi) +
                 (1 - pi) * math.pow(q, xi) * math.pow(1 - q, 1 - xi))

def e_step(pi,p,q,x):
    """
        e步计算
    :param pi: 下一次迭代开始的 pi
    :param p:  下一次迭代开始的 p
    :param q:  下一次迭代开始的 q
    :param x: 观察数据
    :return:
    """
    return [cal_u(pi,p,q,xi) for xi in x]

def m_step(u,x):
    """
     m步计算
    :param u:  m步计算的u
    :param x:  观察数据
    :return:
    """
    pi1=sum(u)/len(u)
    p1=sum([u[i]*x[i] for i in range(len(u))]) / sum(u)
    q1=sum([(1-u[i])*x[i] for i in range(len(u))]) / sum([1-u[i] for i in range(len(u))])
    return [pi1,p1,q1]

def run(observed_x, start_pi, start_p, start_q, iter_num):
    """

    :param observed_x:  观察数据
    :param start_pi:  下一次迭代开始的pi $\pi$
    :param start_p:  下一次迭代开始的p
    :param start_q:  下一次迭代开始的q
    :param iter_num:  迭代次数
    :return:
    """
    for i in range(iter_num):
        u=e_step(start_pi, start_p, start_q, observed_x)
        print (i,[start_pi,start_p,start_q])
        if [start_pi,start_p,start_q]==m_step(u, observed_x):
            break
        else:
            [start_pi,start_p,start_q]=m_step(u, observed_x)
if __name__ =="__main__":
    # 观察数据
    x = [1, 1, 0, 1, 0, 0, 1, 0, 1, 1]
    # 初始化 pi，p q
    [pi, p, q] = [0.4, 0.6, 0.7]
    # 迭代计算
    run(x,pi,p,q,100)

运行结果如下

0 [0.4, 0.6, 0.7]
1 [0.40641711229946526, 0.5368421052631579, 0.6432432432432431]
2 [0.40641711229946537, 0.5368421052631579, 0.6432432432432431]
3 [0.40641711229946537, 0.536842105263158, 0.6432432432432431]