理论参考#

样例介绍#

有c个硬币，每次随机选一个投掷n次，重复执行m次并记录结果。根据投掷结果计算出每个硬币正面的概率。每次选择的硬币未知。

过程介绍#

随机初始化硬币为正的概率 head_p
根据 head_p 求出选择某个硬币的概率 selected_p
根据 selected_p 计算新的硬币概率 head_p
若 head_p 收敛，执行5；否则，执行2
结束

代码实现#

导入库#

1
import random
2
import numpy as np
3
from tqdm import tqdm
4
from collections import Counter
5
import matplotlib.pyplot as plt

设置硬币朝上的真实值#

coin_num：硬币数

coins：硬币正面朝上的真实值

1
# 设置真实值
2
coin_num = 5
3
coins = []
4
for _ in range(coin_num):
5
    coins.append(random.randint(0, 100)/100)
6
coins

模拟投币#

随机抽取一个硬币
投掷 n 次
循环 m 次

1 为正面，0 为反面

【注】：在计算概率时并没有乘以组合数，因此当n过大时，概率会精度丢失变为0导致优化失败

1
0n = 100
2
m = 1000
3
coin_result = np.zeros((m, n), dtype=int)
4

5
c_selected_record = []
6
for i_m in range(m):
7
    # 选择硬币
8
    coin_p = random.choice(coins)
9
    c_selected_record.append(coin_p)
10
    for i_n in range(n):
11
        # 开始投掷
12
        coin_result[i_m, i_n] = 1 if random.random() < coin_p else 0
13

14
coin_result.shape, Counter(c_selected_record)

EM算法#

初始化：随机初始化硬币正面的概率

E步：计算当前硬币的期望
M步：更新硬币参数

为了便于实现，修改概率计算方式，结果不变：

$p^{n_1}\times (1-p)^{n_2} \\ =exp(n_1 \times \log(p) + n_2 \times \log(1-p))$

【注】：预测的硬币并不是一一对应，顺序会变

1
ini_coin_theta = np.array([random.randint(1, 99)/100 for _ in range(coin_num)])
2
# coin_theta = np.array([0.2, 0.9])
3
print('ini coin:', ini_coin_theta)
4

5
def E(coin_theta, coin_result):
6

7
    h_e_sum = np.zeros_like(coin_theta)
8
    t_e_sum = np.zeros_like(coin_theta)
9

10
    h_num = coin_result.sum(1)[:, None]
11
    t_num = coin_result.shape[1] - h_num
12

13
    # 可以评估每个硬币的得分
14
    coin_selected_p = h_num @ np.log(coin_theta[None]) + t_num @ np.log(1 - coin_theta[None])
15
    coin_selected_p = np.exp(coin_selected_p)
16
    coin_selected_p = coin_selected_p / coin_selected_p.sum(1)[:, None]
17

18
    h_e = coin_selected_p  * h_num
19
    t_e = coin_selected_p  * t_num
20

21
    return h_e.sum(0), t_e.sum(0), coin_selected_p
22

23
def M(h_e_sum, t_e_sum):
24
    return h_e_sum / (h_e_sum + t_e_sum)
25

26
max_step=1000
27
coin_result = np.array(coin_result)
28

29
h_e_record = []
30
t_e_record = []
31
theta_record = []
32
delta_record = []
33

34
coin_theta = ini_coin_theta
35
for i in tqdm(range(max_step)):
36
    h_e_sum, t_e_sum, coin_selected_p = E(coin_theta, coin_result)
37

38
    h_e_record.append(h_e_sum)
39
    t_e_record.append(t_e_sum)
40

41
    new_coin_theta = M(h_e_sum, t_e_sum)
42

43
    theta_record.append(new_coin_theta)
44

45
    delta =  ((new_coin_theta - coin_theta)**2).sum()
46

47
    delta_record.append(delta)
48

49
    # print(new_coin_theta)
50
    if delta < 1e-10:
51
        break
52
    coin_theta = new_coin_theta
53

54

55
h_e_record = np.array(h_e_record)
56
t_e_record = np.array(t_e_record)
57
theta_record = np.array(theta_record)
58
delta_record = np.array(delta_record)
59

60
i, coin_theta, coins
61

62
'''
63
(36,
64
 array([0.62988197, 0.43099465, 0.84265886, 0.99086422, 0.53815304]),
65
 [0.84, 0.99, 0.63, 0.44, 0.54])
66
'''

展示参数变化过程#

1
def plot_list(f, x, y, labels, title):
2
    f.set_title(title)
3
    for i in range(y.shape[1]):
4
        f.plot(x, y[:, i], label = labels[i], linestyle='--')
5

6
index = range(0, i+1)
7
labels = list(range(coin_theta.shape[0]))
8
figure, axes = plt.subplots(2, 2, figsize=(12,12), dpi=80)
9

10

11
axes[0, 0].set_title("delta")
12
# 与上一步结果的差别
13
axes[0, 0].plot(index, delta_record, label="delta")
14
# 硬币正面的概率
15
plot_list(axes[0, 1], index, theta_record, labels=labels, title="theta")
16
# 每个硬币正面的加权和
17
plot_list(axes[1, 0], index, h_e_record, labels=labels, title="h_e")
18
# 每个硬币反面的加权和
19
plot_list(axes[1, 1], index, t_e_record, labels=labels, title="t_e")
20

21
for axe in axes:
22
    for a in axe:
23
        a.legend()

在这里插入图片描述

源文件#

HMTTT/EM-

原文链接：CSDN