利用EM算法解决三类硬币抛掷问题
文章目录
- 深度学习与自然语言处理第二次作业——EM算法解决硬币问题
- 一、解题背景
- 二、解题原理
- 三、实验分析
- Step1：初始化
- Step2：E-step
- Step3：M-step
- step4：迭代
- 四、实验代码
- 1.随机投掷硬币数据生成模块
- 2.单次EM模块
- 3.迭代模块
- 五、实验结果与分析
- 1、硬币投掷方法的影响
- 2、初始参数 θ \theta θ的影响
- 六、实验总结
- 实验总体代码
一、解题背景
一个袋子中三种硬币的混合比例为：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
P(X=1)=p;P(X=0)=1-p
P(X=1)=p;P(X=0)=1−p。在模型中，
θ
^
\hat\theta
θ^表示对于当前的模型求得的对应的期望值
θ
^
=
<
s
1
,
s
2
,
p
,
q
,
r
>
\hat\theta=
在模型中，通过实验已知变量
X
X
X为硬币的正反面，隐藏变量
U
U
U为硬币的种类。
在实验开始时，生成N个硬币投掷结果并输入假设的 θ \theta θ值。
Step2：E-step对于第j次迭代出来的第i个 x x x，算出来的隐含变量 u 1 u_1 u1, u 2 u_2 u2表达式如下：
u
e
1
(
x
(
i
)
)
=
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^1 _e(x^{(i)} )=\frac{p^{x^{(i)}}_{(e-1)}(1-p_{(e-1)})^{1-x^{(i)}}s^1_{(e-1)}}{p^{x^{(i)}}_{(e-1)}(1-p_{(e-1)})^{1-x^{(i)}}s^1_{(e-1)}+q^{x^{(i)}}_{(e-1)}(1-q_{(e-1)})^{1-x^{(i)}}s^2_{(e-1)}+r^{x^{(i)}}_{(e-1)}(1-r_{(e-1)})^{1-x^{(i)}}(1-s^1_{(e-1)}-s^2_{(e-1)})}
ue1(x(i))=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)p(e−1)x(i)(1−p(e−1))1−x(i)s(e−1)1
u
e
2
(
x
(
i
)
)
=
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
)
u^2 _e(x^{(i)} )=\frac{q^{x^{(i)}}_{(e-1)}(1-q_{(e-1)})^{1-x^{(i)}}s^2_{(e-1)}} {p^{x^{(i)}}_{(e-1)}(1-p_{(e-1)})^{1-x^{(i)}}s^1_{(e-1)}+q^{x^{(i)}}_{(e-1)}(1-q_{(e-1)})^{1-x^{(i)}}s^2_{(e-1)}+r^{x^{(i)}}_{(e-1)}(1-r_{(e-1)})^{1-x^{(i)}}(1-s^1_{(e-1)}-s^2_{(e-1)})}
ue2(x(i))=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)q(e−1)x(i)(1−q(e−1))1−x(i)s(e−1)2
对于第e次迭代，根据极大似然估计的
θ
^
(
e
)
{\hat{\theta}}_{(e)}
θ^(e)为
θ
^
(
e
)
[
0
]
=
s
(
e
)
1
=
∑
i
=
1
i
=
N
u
(
e
)
1
(
x
(
i
)
)
N
{\hat{\theta}}_{(e)}[0]=s^1_{(e)}=\frac{\sum^{i=N}_{i=1}u^1_{(e)}(x^{(i)})}N
θ^(e)[0]=s(e)1=N∑i=1i=Nu(e)1(x(i))
θ
^
(
e
)
[
1
]
=
s
(
e
)
2
=
∑
i
=
1
i
=
N
u
(
e
)
2
(
x
(
i
)
)
N
{\hat{\theta}}_{(e)}[1]=s^2_{(e)}=\frac{\sum^{i=N}_{i=1}u^2_{(e)}(x^{(i)})}N
θ^(e)[1]=s(e)2=N∑i=1i=Nu(e)2(x(i))
θ
^
(
e
)
[
2
]
=
p
(
e
)
=
∑
i
=
1
i
=
N
x
(
i
)
∗
u
(
e
)
1
(
x
(
i
)
)
∑
i
=
1
i
=
N
u
(
e
)
1
(
x
(
i
)
)
{\hat{\theta}}_{(e)}[2]=p_{(e)}=\frac{\sum^{i=N}_{i=1}x^{(i)}*u^1_{(e)}(x^{(i)})}{\sum^{i=N}_{i=1}u^1_{(e)}(x^{(i)})}
θ^(e)[2]=p(e)=∑i=1i=Nu(e)1(x(i))∑i=1i=Nx(i)∗u(e)1(x(i))
θ
^
(
e
)
[
3
]
=
q
(
e
)
=
∑
i
=
1
i
=
N
x
(
i
)
∗
u
(
e
)
2
(
x
(
i
)
)
∑
i
=
1
i
=
N
u
(
e
)
2
(
x
(
i
)
)
{\hat{\theta}}_{(e)}[3]=q_{(e)}=\frac{\sum^{i=N}_{i=1}x^{(i)}*u^2_{(e)}(x^{(i)})}{\sum^{i=N}_{i=1}u^2_{(e)}(x^{(i)})}
θ^(e)[3]=q(e)=∑i=1i=Nu(e)2(x(i))∑i=1i=Nx(i)∗u(e)2(x(i))
θ
^
(
e
)
[
4
]
=
r
(
e
)
=
∑
i
=
1
i
=
N
x
(
i
)
∗
(
1
−
u
(
e
)
1
(
x
(
i
)
)
−
u
(
e
)
2
(
x
(
i
)
)
)
N
−
∑
i
=
1
i
=
N
u
(
e
)
1
(
x
(
i
)
)
−
∑
i
=
1
i
=
N
u
(
e
)
2
(
x
(
i
)
)
{\hat{\theta}}_{(e)}[4]=r_{(e)}=\frac{\sum^{i=N}_{i=1}x^{(i)}*(1-u^1_{(e)}(x^{(i)})-u^2_{(e)}(x^{(i)}))}{N-\sum^{i=N}_{i=1}u^1_{(e)}(x^{(i)})-\sum^{i=N}_{i=1}u^2_{(e)}(x^{(i)})}
θ^(e)[4]=r(e)=N−∑i=1i=Nu(e)1(x(i))−∑i=1i=Nu(e)2(x(i))∑i=1i=Nx(i)∗(1−u(e)1(x(i))−u(e)2(x(i)))
不断进行Step2：E-Step和Step3：M-Step直到算法收敛
四、实验代码 1.随机投掷硬币数据生成模块
根据题目要求，需要自己指定系数s1, s2, p, q, r，生成N个投掷硬币的结果。
代码中的true_alphaL 、 true_thetaL分别表示抽到某类硬币的概率和某类硬币正面向上的概率；N表示投掷硬币的数量，size表示每枚硬币每轮的投掷次数，生成的投掷硬币结果为
N
=
N
∗
s
i
z
e
N=N*size
N=N∗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
根据前一次返回的 θ \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]
进行最大化步骤，当迭代后一次的结果与前一次的结果变化值小于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时，没有可行解产生。
本实验的最终输出
θ
=
<
s
1
,
s
2
,
p
,
q
,
r
>
\theta=
在本实验中首先根据真实参数值
θ
\theta
θ产生的N个硬币投掷结果，再根据EM算法计算出估计的参数值。若设置的初始参数不同，则最终得到的
θ
^
\hat\theta
θ^与真实参数值
θ
\theta
θ可能存在相应差异。
- 例如，针对真实参数值
θ
=
<
s
1
,
s
2
,
p
,
q
,
r
>
=
<
0.2
,
0.3
,
0.5
,
0.7
,
0
,
9
>
\theta=
此时 θ ^ \hat\theta θ^与真实参数值 θ \theta θ近似相等。 - 当初始尝试值改为
θ
=
<
0.8
,
0.1
,
0.3
,
0.6
,
0.4
>
\theta=<0.8,0.1,0.3,0.6,0.4>
θ=<0.8,0.1,0.3,0.6,0.4>时，最终得到
θ
^
\hat\theta
θ^与真实参数值
θ
\theta
θ存在着对应关系的差异，即真实的s2的值应为
<
0.3
,
0.7
>
<0.3,0.7>
<0.3,0.7>，此时对应的硬币应为硬币s3，即此时s2与s3硬币的概率发生改变：
考虑原因应为初始参数设计时，s2的初始参数更靠近s3的真实参数值，进而经过迭代后，s2与s3的参数值发生互换，但对实验结果没有影响，只是这两者顺序发生改变。 - 初始尝试值改为
θ
=
<
0.8
,
0.1
,
0.6
,
0.6
,
0.9
>
\theta=<0.8,0.1,0.6,0.6,0.9>
θ=<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)
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）