深度学习与自然语言处理第二次作业——EM算法解决硬币问题

深度学习与自然语言处理第二次作业——EM算法解决硬币问题

利用EM算法解决三类硬币抛掷问题

---

---

一、解题背景

一个袋子中三种硬币的混合比例为：s1, s2与1-s1-s2 (0<=si<=1), 三种硬币掷出正面的概率分别为：p, q, r。 （1）自己指定系数s1, s2, p, q, r，生成N个投掷硬币的结果（由01构成的序列，其中1为正面，0为反面），利用EM算法来对参数进行估计并与预先假定的参数进行比较。

---

二、解题原理

最大期望算法（Expectation-maximization algorithm，又译为期望最大化算法），是在概率模型中寻找参数最大似然估计或者最大后验估计的算法，其中概率模型依赖于无法观测的隐性变量。
最大期望算法经过两个步骤交替进行计算：
第一步是计算期望（E），利用对隐藏变量的现有估计值，计算其最大似然估计值；
第二步是最大化（M），最大化在E步上求得的最大似然值来计算参数的值。M步上找到的参数估计值被用于下一个E步计算中，这个过程不断交替进行，直到算法收敛。
具体收敛性证明，可见前人分析。

---

三、实验分析

该模型可视为三个二项分布组合，定义1为正面，0为反面，由数学知识，其二项分布的概率密度函数为$P(X=1)=p;P(X=0)=1-p$。在模型中，$\hat{\theta }$表示对于当前的模型求得的对应的期望值$\hat{\theta }=<s1,s2,p,q,r>$，其中分别为硬币s1的概率，硬币s2的概率，抛s1硬币为正面的概率，抛s2硬币为正面的概率，抛s3硬币为正面的概率。
在模型中，通过实验已知变量$X$为硬币的正反面，隐藏变量$U$为硬币的种类。

Step1：初始化

在实验开始时，生成N个硬币投掷结果并输入假设的$\theta$值。

Step2：E-step

对于第j次迭代出来的第i个$x$，算出来的隐含变量$u_1$,$u_2$表达式如下：

$$u_e^1(x^{(i)})=\frac{p_{(e-1)}^{x^{(i)}}(1-p_{(e-1)}{)}^{1-x^{(i)}}{s}_{(e-1)}^1}{p_{(e-1)}^{x^{(i)}}(1-p_{(e-1)}{)}^{1-x^{(i)}}{s}_{(e-1)}^1+q_{(e-1)}^{x^{(i)}}(1-q_{(e-1)}{)}^{1-x^{(i)}}{s}_{(e-1)}^2+r_{(e-1)}^{x^{(i)}}(1-r_{(e-1)}{)}^{1-x^{(i)}}(1-s_{(e-1)}^1-s_{(e-1)}^2)}$$
$$u_e^2(x^{(i)})=\frac{q_{(e-1)}^{x^{(i)}}(1-q_{(e-1)}{)}^{1-x^{(i)}}{s}_{(e-1)}^2}{p_{(e-1)}^{x^{(i)}}(1-p_{(e-1)}{)}^{1-x^{(i)}}{s}_{(e-1)}^1+q_{(e-1)}^{x^{(i)}}(1-q_{(e-1)}{)}^{1-x^{(i)}}{s}_{(e-1)}^2+r_{(e-1)}^{x^{(i)}}(1-r_{(e-1)}{)}^{1-x^{(i)}}(1-s_{(e-1)}^1-s_{(e-1)}^2)}$$

Step3：M-step

对于第e次迭代，根据极大似然估计的${\hat{\theta }}_{(e)}$为
$${\hat{\theta }}_{(e)}[0]=s_{(e)}^1=\frac{{\sum }_{i=1}^{i=N}{u}_{(e)}^1(x^{(i)})}{N}$$
$${\hat{\theta }}_{(e)}[1]=s_{(e)}^2=\frac{{\sum }_{i=1}^{i=N}{u}_{(e)}^2(x^{(i)})}{N}$$
$${\hat{\theta }}_{(e)}[2]=p_{(e)}=\frac{{\sum }_{i=1}^{i=N}{x}^{(i)}\ast u_{(e)}^1(x^{(i)})}{{\sum }_{i=1}^{i=N}{u}_{(e)}^1(x^{(i)})}$$
$${\hat{\theta }}_{(e)}[3]=q_{(e)}=\frac{{\sum }_{i=1}^{i=N}{x}^{(i)}\ast u_{(e)}^2(x^{(i)})}{{\sum }_{i=1}^{i=N}{u}_{(e)}^2(x^{(i)})}$$
$${\hat{\theta }}_{(e)}[4]=r_{(e)}=\frac{{\sum }_{i=1}^{i=N}{x}^{(i)}\ast (1-u_{(e)}^1(x^{(i)})-u_{(e)}^2(x^{(i)}))}{N-{\sum }_{i=1}^{i=N}{u}_{(e)}^1(x^{(i)})-{\sum }_{i=1}^{i=N}{u}_{(e)}^2(x^{(i)})}$$

step4：迭代

不断进行Step2：E-Step和Step3：M-Step直到算法收敛

---

四、实验代码

1.随机投掷硬币数据生成模块

根据题目要求，需要自己指定系数s1, s2, p, q, r，生成N个投掷硬币的结果。
代码中的true_alphaL 、 true_thetaL分别表示抽到某类硬币的概率和某类硬币正面向上的概率；N表示投掷硬币的数量，size表示每枚硬币每轮的投掷次数，生成的投掷硬币结果为$N=N\ast size$
在实验过程中发现，起初size为1时，算法收敛得比较慢，且结果可能不收敛。因此设置size，增加每枚硬币的投掷次数，加快收敛速度。

代码如下（数据生成）：

	true_alphaL = [0.2, 0.3, 0.5]
    true_thetaL = [0.5, 0.7, 0.9]
    N = 2500    #拥有的硬币数量
    size = 100  #单枚硬币每轮的投掷次数
    data = []
    for alpha, theta_1 in zip(true_alphaL, true_thetaL):
        for _ in range(int(alpha * N)):
            data.append(np.random.binomial(1, theta_1, size=size))
    O = data

2.单次EM模块

根据前一次返回的$\theta$值依次进行E-Step、M-step，并返回这一轮得到的新$\hat{\theta }$值。

代码如下（单次EM）：

def em_single(theta, O):  # 对于当前的模型求对应的期望值(估算步骤)
    # 改为5个值分别为硬币s1的概率、抛s2硬币的概率、抛s1硬币为正面的概率、抛s2硬币为正面的概率，抛s3硬币为正面的概率
    pi_1 = theta[0]
    pi_2 = theta[1]
    h1 = theta[2]
    h2 = theta[3]
    h3 = theta[4]
    
    new_pi_1 = 0
    new_pi_2 = 0
    new_h1 = 0
    new_h2 = 0
    new_h3 = 0
    PA = []
    PB = []
    head = []

    for o in O:
        t = len(o)
        cnt = o.sum()  # 正面的数量
        head.append(cnt)
    for o in O:
        t = len(o)
        heads = o.sum()  # 正面的数量
        tails = t - heads  # 反面的数量
#e_step
        u_1 = (pi_1 * math.pow(h1, heads) * math.pow(1 - h1, tails)) /\
              (pi_1 * math.pow(h1, heads) * math.pow(1 - h1, tails) + (pi_2) * math.pow(h2, heads) * math.pow(1 - h2, tails)+(1-pi_1-pi_2) * math.pow(h3, heads) * math.pow(1 - h3, tails))
        u_2 = (pi_2 * math.pow(h2, heads) * math.pow(1 - h2, tails)) / \
              (pi_1 * math.pow(h1, heads) * math.pow(1 - h1, tails) + (pi_2) * math.pow(h2, heads) * math.pow(1 - h2, tails)+(1-pi_1-pi_2) * math.pow(h3, heads) * math.pow(1 - h3, tails))
        PA.append(u_1)
        PB.append(u_2)
        
    l = len(O)
#m_step
    new_pi_1 = sum(PA) / N
    new_pi_2 = sum(PB) / N

    for i in range(l):
        new_h1 += PA[i] * head[i] / t
        new_h2 += PB[i] * head[i] / t
        new_h3 += (1 - PA[i]-PB[i]) * head[i] / t
        # new_h1 += PA[i] * head[i]
        # new_h2 += PB[i] * head[i]
        # new_h3 += (1 - PA[i] - PB[i]) * head[i]
    new_h1 = new_h1 / (sum(PA))
    new_h2 = new_h2 / (sum(PB))
    new_h3 = new_h3 / (l - sum(PB) - sum(PA))
    # 因为有c1 c2两种可能，总观测数减去c1的即为c2的
    return [new_pi_1, new_pi_2, new_h1, new_h2,  new_h3]

3.迭代模块

进行最大化步骤，当迭代后一次的结果与前一次的结果变化值小于tol时或者进行1000次迭代后，则得到最终的收敛结果。

代码如下（最大化步骤）

def em(theta, Y, tol):  # 最大化步骤
    j = 0
    while j < 1000:
        new_theta = em_single(theta, Y)
        change = np.abs(new_theta[0] - theta[0])
        if change < tol:
            break
        else:
            theta = new_theta
            j = j + 1
        print("迭代后的参数为:")
        print(new_theta)
        print(j)
    return new_theta

五、实验结果与分析

进行了多组实验，发现实验收敛结果受到多种因素的影响。

1、硬币投掷方法的影响

当N=25000，size=1时，迭代结果如下（实验未完全进行），且结果收敛速度较慢：

增加size值，即N=2500，size=100时，经过39次迭代即得到与真实值收敛的参数结果。

因此，从中我们可以知道：

• 系统的收敛速度与硬币的投掷方法有关，单个硬币投掷多次的收敛速度明显优于多次抽取硬币进行投掷，即size的值越大，收敛速度越优；
• 但是应注意到的是，size的大小不能无穷大，当size大小过大时，会出现系统无可行解的情况。例如，当N=250，size=1000时，没有可行解产生。
  

2、初始参数$\theta$的影响

本实验的最终输出$\theta =<s1,s2,p,q,r>$受到初始参数设置的影响，根据EM模型的原理我们可知，最终得到的$\hat{\theta }$值向着最接近估计值靠近。
在本实验中首先根据真实参数值$\theta$产生的N个硬币投掷结果，再根据EM算法计算出估计的参数值。若设置的初始参数不同，则最终得到的$\hat{\theta }$与真实参数值$\theta$可能存在相应差异。

1. 例如，针对真实参数值$\theta =<s1,s2,p,q,r>=<0.2,0.3,0.5,0.7,0,9>$产生一组$N=250000$的硬币投掷数据，当初始参数$\theta =<0.1,0.3,0.1,0.6,0.8>$时经过37次迭代，产生以下结果：
   
   此时$\hat{\theta }$与真实参数值$\theta$近似相等。
2. 当初始尝试值改为$\theta =<0.8,0.1,0.3,0.6,0.4>$时，最终得到$\hat{\theta }$与真实参数值$\theta$存在着对应关系的差异，即真实的s2的值应为$<0.3,0.7>$，此时对应的硬币应为硬币s3，即此时s2与s3硬币的概率发生改变：
   
   考虑原因应为初始参数设计时，s2的初始参数更靠近s3的真实参数值，进而经过迭代后，s2与s3的参数值发生互换，但对实验结果没有影响，只是这两者顺序发生改变。
3. 初始尝试值改为$\theta =<0.8,0.1,0.6,0.6,0.9>$时，最终得到$\hat{\theta }$与真实参数值$\theta$存在着极大的差异：
   
   考虑原因应为由于初始参数设计时，s1和s2的正面向上概率相似，此时EM算法无法识别s1与s2 两种硬币，而将这两种硬币视为同一种硬币。

实验结果表格如下所示：

参数值	真实参数值	初始参数为<0.1,0.3,0.1,0.6,0.8>	初始参数为<0.8,0.1,0.3,0.6,0.4>	初始参数为<0.8, 0.1, 0.6, 0.6, 0.9>
s1	0.2	0.20046934103514671	0.20127147448954955	0.4300894611937678
s2	0.3	0.29995687083883144	0.4981815495712238	0.05376118264922097
p	0.5	0.4980744170221924	0.4995940693133312	0.6141768099240248
q	0.7	0.7003422897528382	0.9004644692018112	0.6141768099240248
r	0.9	0.8996446145996789	0.6991625178254812	0.8978140717545678

---

六、实验总结

本实验使用EM算法对三种硬币的参数进行估计，实验结果看出：

• 算法的收敛速度与硬币投掷方式相关
• 最终参数$\hat{\theta }$收敛性与初始参数$\theta$的设定相关

因此，在投掷阶段，我们可以对一枚硬币进行多次投掷；在实验使用EM算法开始前，我们可先针对袋子中的硬币性质进行估计，对初始参数进行设计，使之能得到更准确的结果。

---

实验总体代码

实验代码使用python实现

import numpy as np
import math
import random

def em_single(theta, O):  # 对于当前的模型求对应的期望值(估算步骤)
    # 改为5个值分别为硬币s1的概率、抛s2硬币的概率、抛s1硬币为正面的概率、抛s2硬币为正面的概率，抛s3硬币为正面的概率
    pi_1 = theta[0]
    pi_2 = theta[1]
    h1 = theta[2]
    h2 = theta[3]
    h3 = theta[4]

    new_pi_1 = 0
    new_pi_2 = 0
    new_h1 = 0
    new_h2 = 0
    new_h3 = 0
    PA = []
    PB = []
    head = []

    for o in O:
        t = len(o)
        cnt = o.sum()  # 正面的数量
        head.append(cnt)
    for o in O:
        t = len(o)
        heads = o.sum()  # 正面的数量
        tails = t - heads  # 反面的数量

        u_1 = (pi_1 * math.pow(h1, heads) * math.pow(1 - h1, tails)) /\
              (pi_1 * math.pow(h1, heads) * math.pow(1 - h1, tails) + (pi_2) * math.pow(h2, heads) * math.pow(1 - h2, tails)+(1-pi_1-pi_2) * math.pow(h3, heads) * math.pow(1 - h3, tails))
        u_2 = (pi_2 * math.pow(h2, heads) * math.pow(1 - h2, tails)) / \
              (pi_1 * math.pow(h1, heads) * math.pow(1 - h1, tails) + (pi_2) * math.pow(h2, heads) * math.pow(1 - h2, tails)+(1-pi_1-pi_2) * math.pow(h3, heads) * math.pow(1 - h3, tails))
        PA.append(u_1)
        PB.append(u_2)


    l = len(O)


    new_pi_1 = sum(PA) / N
    new_pi_2 = sum(PB) / N

    for i in range(l):
        new_h1 += PA[i] * head[i] / t
        new_h2 += PB[i] * head[i] / t
        new_h3 += (1 - PA[i]-PB[i]) * head[i] / t
        # new_h1 += PA[i] * head[i]
        # new_h2 += PB[i] * head[i]
        # new_h3 += (1 - PA[i] - PB[i]) * head[i]
    new_h1 = new_h1 / (sum(PA))
    new_h2 = new_h2 / (sum(PB))
    new_h3 = new_h3 / (l - sum(PB) - sum(PA))
    # 因为有c1 c2两种可能，总观测数减去c1的即为c2的
    return [new_pi_1, new_pi_2, new_h1, new_h2,  new_h3]


def em(theta, Y, tol):  # 最大化步骤
    j = 0
    while j < 1000:
        new_theta = em_single(theta, Y)
        change = np.abs(new_theta[0] - theta[0])
        if change < tol:
            break
        else:
            theta = new_theta
            j = j + 1
        print("迭代后的参数为:")
        print(new_theta)
        print(j)
    return new_theta



if __name__ == "__main__":
    # 硬币投掷结果

    true_alphaL = [0.2, 0.3, 0.5]
    true_thetaL = [0.5, 0.7, 0.9]
    N = 2500   #拥有的硬币数量
    size = 100  #单枚硬币每轮的投掷次数
    data = []
    for alpha, theta_1 in zip(true_alphaL, true_thetaL):
        for _ in range(int(alpha * N)):
            data.append(np.random.binomial(1, theta_1, size=size))
    O = data
    #print(O)

    theta = [0.8, 0.1, 0.6, 0.6, 0.9]
    #theta = [0.1, 0.3, 0.1, 0.6, 0.8]
    #theta = [0.001, 0.001, 0.001, 0.001, 0.8]
    t = len(theta)
    print("初始参数为:")
    print(theta)
    theta = em(theta, O, 1e-15)