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 , … , x n X={x_1,x_2,…,x_n} X=x1,x2,,xn

写出生成一个硬币时的概率:
P ( x ∣ θ ) = ∑ z P ( x , z ∣ θ ) = ∑ z P ( z ∣ π ) P ( x ∣ z , θ ) = π p x ( 1 − p ) 1 − x + ( 1 − π ) q x ( 1 − q ) 1 − x \begin{aligned}P(x|θ) &=∑zP(x,z|θ) = ∑zP(z|π)P(x|z,θ) \\ &= πp^x(1−p)^{1−x}+(1−π)q^x(1−q)^{1−x} \end{aligned} P(xθ)=zP(x,zθ)=zP(zπ)P(xz,θ)=πpx(1p)1x+(1π)qx(1q)1x
有了一个硬币的概率,我们就可以写出所有观察数据的log似然函数:
L ( θ ∣ X ) = l o g P ( X ∣ θ ) = ∑ j = 1 n l o g [ π p x j ( 1 − p ) 1 − x j + ( 1 − π ) q x j ( 1 − q ) 1 − x j ] \begin{aligned} L(θ|X)=logP(X|θ)=∑_{j=1}^nlog[πp^{x_j}(1−p)^{1−x_j} + (1−π)q^{x_j}(1−q)^{1−x_j}] \end{aligned} L(θX)=logP(Xθ)=j=1nlog[πpxj(1p)1xj+(1π)qxj(1q)1xj]
然后求似然函数最大值
θ ∗ = a r g m a x L ( θ ∣ X ) θ^*= arg\quad maxL(θ|X) θ=argmaxL(θX)
其中 L ( θ ∣ X ) = l o g P ( X ∣ θ ) = l o g ∑ Z P ( X , Z ∣ θ ) L(θ|X)=logP(X|θ)=log∑_ZP(X,Z|θ) L(θX)=logP(Xθ)=logZP(X,Zθ)。因为 l o g log log里面带着加和所以这个极大似然是求不出解析解的。 如果我们知道隐藏变量Z的话,求以下的似然会容易很多:
L(θ|X,Z)=logP(X,Z|θ)
但是隐藏变量Z的值谁也不知道,所以我们转用EM算法求它的后验概率期望 E Z ∣ X , θ ( L ( θ ∣ X , Z ) ) E_{Z|X,θ}(L(θ|X,Z)) EZX,θ(L(θX,Z))的最大值。

  • E步::假设当前模型的参数为 π , p , q π,p,q π,p,q时,隐含变量来自于硬币B的后验概率
    μ = P ( Z ∣ X , θ ) = π p x i ( 1 − p ) 1 − x i π p x i ( 1 − p ) 1 − x i + ( 1 − π ) q x i ( 1 − q ) 1 − x i \begin{aligned} μ=P(Z|X,θ)= \frac{πp^{x_i}(1−p)^{1−x_i}}{πp^{x_i}(1−p)^{1−x_i}+(1−π)q^{x_i}(1−q)^{1−x_i}} \end{aligned} μ=P(ZX,θ)=πpxi(1p)1xi+(1π)qxi(1q)1xiπpxi(1p)1xi
    那么隐含变量来自于硬币C的后验概率自然为 1 − μ 1−μ 1μ
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 , θ ( L ( θ ∣ X , Z ) ) = ∑ j = 1 n [ μ l o g π p x i ( 1 − p ) 1 − x i μ + ( 1 − μ ) l o g ( 1 − π ) q x i ( 1 − q ) 1 − x i 1 − μ ] E_{Z|X,θ}(L(θ|X,Z))=∑_{j=1}^n [μlog \frac{πp^{x_i}(1−p)^{1−x_i}}{μ} + (1−μ)\frac{log(1−π)q^{x_i}(1−q)^{1−x_i}}{1−μ}] EZX,θ(L(θX,Z))=j=1n[μlogμπpxi(1p)1xi+(1μ)1μlog(1π)qxi(1q)1xi]
    要注意这里把 μ μ μ看做固定值。然后我们分别求偏导,获得参数 π , p , q π,p,q π,p,q的新估计值
    π = 1 n ∑ j = 1 n μ j p = ∑ j = 1 n μ j x j ∑ j = 1 n μ j q = ∑ j = 1 n ( 1 − μ j ) x j ∑ j = 1 n ( 1 − μ j ) \begin{aligned} π &=\frac{1}{n} ∑_{j=1}^nμ_j \\ p&=\frac{∑_{j=1}^nμ_jx_j}{∑_{j=1}^nμ_j }\\ q&=\frac{∑^n_{j=1}(1−μ_j)xj}{∑^n_{j=1}(1−μ_j)} \end{aligned} πpq=n1j=1nμj=j=1nμjj=1nμjxj=j=1n(1μj)j=1n(1μj)xj
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]

算法首先选取参数初始值 θ ( 0 ) = π ( 0 ) , p ( 0 ) , q ( 0 ) θ(0)=π(0),p(0),q(0) θ(0)=π(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]

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值