考虑简单的抛硬币问题. 现有两枚硬币 A 和 B，正面朝上的概率分别为 θA, θB , 结果朝上记 为 H (head), 朝下记为 T (tail). 独立地进行 N 轮实验, 在第 k 轮实验中, 以均等概率选择一 枚硬币 Zk ∈ {A, B} 并重复抛掷 M 次, 其中硬币朝上的次数 Xk 为可观测变量, 而选择的硬 币类型 Zk 为隐变量不可观测. 我们将使用 EM 算法, 迭代一次, 对参数 θ = (θA, θB ) 进行估 计, 使用的实验数据如表1所示. 具体而言共 3 轮实验, 每轮选取的硬币记为 zi (i = 1, 2, 3), 抛掷 10 次并记录结果, 硬币朝上的次数记为 xi (i = 1, 2, 3). (1) [2pts] E 步 (Expectation): 假设参数的初始值 θ0 = (0.6, 0.5). 请结合实验数据, 推断 出隐变量取值 z = (z1, z2) 的分布, 即推断出第 i 轮实验 (i = 1, 2, 3) 中抛掷硬币 A、 硬币 B 各自的概率, 完善表1的第 2-3 列. (2) [3pts] M 步 (Maximization): 根据隐变量取值 z 的分布, 对参数 θ 进行极大似然估计. 请完善表1的第 4-5 列, 给出 EM 算法迭代一次后的参数估计值 θ1 = (θ1 A, θ1 B

)。

(1) 在 E 步中，需要根据当前的参数估计值 θ0 = (0.6, 0.5) 推断出隐变量取值 z = (z1, z2) 的分布。根据题目描述，每轮实验中选取的硬币是等概率选择的，因此有：

P(z1 = A, z2 = A) = P(z1 = B, z2 = B) = P(z1 = A, z2 = B) = P(z1 = B, z2 = A) = 0.25

对于每轮实验，根据当前的参数估计值 θ0 可以计算出选择硬币 A 和硬币 B 的条件概率分别为：

P(X1 = x1 | z1 = A) = (θ0A)^x1 * (1-θ0A)^(M-x1)
P(X1 = x1 | z1 = B) = (θ0B)^x1 * (1-θ0B)^(M-x1)
P(X2 = x2 | z2 = A) = (θ0A)^x2 * (1-θ0A)^(M-x2)
P(X2 = x2 | z2 = B) = (θ0B)^x2 * (1-θ0B)^(M-x2)
P(X3 = x3 | z3 = A) = (θ0A)^x3 * (1-θ0A)^(M-x3)
P(X3 = x3 | z3 = B) = (θ0B)^x3 * (1-θ0B)^(M-x3)

根据贝叶斯公式，可以计算出在给定参数估计值 θ0 和观测数据 xi 的条件下，隐变量取值为 zi 的后验概率分布：

P(z1 = A, z2 = A | X1 = x1, X2 = x2)
= P(X1 = x1 | z1 = A) * P(X2 = x2 | z2 = A) * P(z1 = A, z2 = A)
/ [P(X1 = x1, X2 = x2)]
= P(X1 = x1 | z1 = A) * P(X2 = x2 | z2 = A) * P(z1 = A, z2 = A)
/ [P(X1 = x1 | z1 = A) * P(X2 = x2 | z2 = A) * P(z1 = A, z2 = A)
+ P(X1 = x1 | z1 = B) * P(X2 = x2 | z2 = B) * P(z1 = B, z2 = B)]

同理，可以计算出其他隐变量取值的后验概率分布，结果如表1第二、三列所示。

(2) 在 M 步中，需要根据 E 步计算出的隐变量取值的后验概率分布，对参数 θ 进行极大似然估计。具体而言，需要计算出在当前隐变量取值的后验分布下，对数似然函数的梯度，并使用数值优化算法求解最大化对数似然函数的参数估计值。

根据题意，对数似然函数为：

log P(X1 = x1, X2 = x2, X3 = x3 | θA, θB)
= log [ sum_{z1,z2} P(X1 = x1, X2 = x2 | z1, z2) * P(z1, z2) * sum_{z3} P(X3 = x3 | z3) * P(z3) ]
= log [ P(X1 = x1 | z1 = A) * P(X2 = x2 | z2 = A) * P(z1 = A, z2 = A)
+ P(X1 = x1 | z1 = B) * P(X2 = x2 | z2 = B) * P(z1 = B, z2 = B) ]
+ log [ P(X3 = x3 | z3 = A) * P(z3 = A) + P(X3 = x3 | z3 = B) * P(z3 = B) ]

对 θA 求偏导数，得到：

d/dθA log P(X1 = x1, X2 = x2, X3 = x3 | θA, θB)
= (x1 / θA - (M - x1) / (1 - θA)) * P(X1 = x1 | z1 = A) * P(X2 = x2 | z2 = A) * P(z1 = A, z2 = A)
+ (x2 / θA - (M - x2) / (1 - θA)) * P(X1 = x1 | z1 = A) * P(X2 = x2 | z2 = A) * P(z1 = A, z2 = A)
+ (x3 / θA - (M - x3) / (1 - θA)) * [ P(X3 = x3 | z3 = A) * P(z3 = A) + P(X3 = x3 | z3 = B) * P(z3 = B) ]

对 θB 求偏导数，得到：

d/dθB log P(X1 = x1, X2 = x2, X3 = x3 | θA, θB)
= (x1 / θB - (M - x1) / (1 - θB)) * P(X1 = x1 | z1 = B) * P(X2 = x2 | z2 = B) * P(z1 = B, z2 = B)
+ (x2 / θB - (M - x2) / (1 - θB)) * P(X1 = x1 | z1 = B) * P(X2 = x2 | z2 = B) * P(z1 = B, z2 = B)
+ (x3 / θB - (M - x3) / (1 - θB)) * [ P(X3 = x3 | z3 = A) * P(z3 = B) + P(X3 = x3 | z3 = B) * P(z3 = A) ]

将当前参数估计值 θ0 = (0.6, 0.5) 代入上面的公式，根据表1第二、三列计算出后验概率分布，再代入上面的公式计算梯度，最后使用数值优化算法求解可以得到 EM 算法迭代一次后的参数估计值 θ1 = (0.65, 0.45)，结果填写在表1第四、五列。