EM 算法实验
最大似然估计
在已知每次实验的硬币是 A 和 B 的情况下，我们可以直接使用最大似然估计求解 A、B 硬币正面出现的概率 $ p_A $ 和 $ p_B $。由于每次实验都是独立地抛十次，所以可以知道实验结果服从二项分布 $ B(10, p_A) $ 和 $ B(10,p_B) $。
设 A 硬币的第 $i$ 次实验正面朝上的次数为 $y_A^{(i)}$， B 硬币的第 $i$ 次实验正面朝上的次数为 $y_B^{(i)}$，分别进行了 $n_A$ 次实验和 $n_B$ 次实验，可以得到 A 的似然估计函数为：
L_A=\prod^{n_A}_{i=1}{10\choose y^{(i)}_A}p_A^{y^{(i)}_A}(1-p_A)^{10 – y^{(i)}_A}
$$
B 的似然估计函数为：
L_B=\prod^{n_B}_{i=1}{10\choose y^{(i)}_B}p_B^{y^{(i)}_B}(1-p_B)^{10 – y^{(i)}_B}
$$
分别取对数并对 $p_A$ 求导后得到：
\frac{\partial\ln L_A}{\partial p_A}=\sum^{n_A}_{i=1}(\frac{y_A^{(i)}}{p_A}-\frac{10 -y_A^{(i)}}{1 – p_A})
$$
令导数为 0，得到：
p_A = \frac{\sum_{i=1}^{n_A}y_A^{(i)}}{\sum_{i=1}^{n_A}y_A^{(i)}+\sum_{i=1}^{n_A}(10-y_A^{(i)})}
$$
同理可得：
p_B = \frac{\sum_{i=1}^{n_B}y_B^{(i)}}{\sum_{i=1}^{n_B}y_B^{(i)}+\sum_{i=1}^{n_B}(10-y_B^{(i)})}
$$
代入数据，即求得：
p_A = \frac{9+8+7}{9+8+7+1+2+3}=\frac{24}{30}=0.8\\
p_B = \frac{5+4}{5+4+5+6}=\frac{9}{20}=0.45
$$
EM 算法
而在不知道每次实验选取的是哪个硬币的时候，我们就没有办法直接对 A 和 B 进行最大似然估计。这时候就需要使用 EM 算法，通过引入一个隐变量 $Z$ 来对硬币概率建模，并规定 $Z\in \{0, 1\}$，$Z=0$ 表示当前硬币是硬币 A，$Z=1$ 表示当前硬币是硬币 B，并记 $P(Z=0)=\pi$。此时模型参数 $\theta=(\pi,p_A,p_B)$，可以知道，$Y$ 与 $Z$ 是相关的，要求 $P(Y|\theta)$，就需要求 $P(Y, Z|\theta)$ 的边缘分布。
记 $z_i$ 和 $y_i$ 分别为第 $i$ 次实验选择的硬币（未知的）和正面朝上的次数，那么，对于所有实验而言，似然函数就变成了
P(Y|\theta)=\prod^{n}_{i=1}P(Y=y_i|\theta)=\prod^{n}_{i=1}\sum_{z_i \in {0,1} }P(Y=y_i,Z=z_i|\theta)=\prod_{i=1}^{n}\sum_{z_i\in{0,1}}P(Z=z_i|\theta)P(Y=y_i|Z=z_i,\theta)\\
=\prod_{i=1}^{n}[\pi{10\choose y_i}\ p_A^{y_i}(1-p_A)^{10-y_i}+(1-\pi){10\choose y_i}\ p_B^{y_i}(1-p_B)^{10-y_i}]
$$
而极大似然估计求参数即为：
\hat\theta=\arg\max_\theta\ln P(Y|\theta)
$$
这个问题很难求出解析解，所以需要使用 EM 迭代法求近似解。
因为
\ln P(Y|\theta)=\sum_{i=1}^{n}\ln[\pi{10\choose y_i}\ p_A^{y_i}(1-p_A)^{10-y_i}+(1-\pi){10\choose y_i}\ p_B^{y_i}(1-p_B)^{10-y_i}]
$$
含有“和的对数”，极难求解，通过 Jensen 不等式可得
\ln P(Y|\theta) = \sum^{n}_{i=1}\ln \sum_{z_i \in {0,1}}P(Y=y_i,Z=z_i|\theta)\\
=\sum^{n}_{i=1}\ln \sum_{z_i \in {0,1}}Q_i(z_i,\theta)\frac{P(Y=y_i,Z=z_i|\theta)}{Q_i(z_i,\theta)}\\
\ge\sum^{n}_{i=1} \sum_{z_i \in {0,1}}Q_i(z_i,\theta)\ln \frac{P(Y=y_i,Z=z_i|\theta)}{Q_i(z_i,\theta)}=Q(\theta;\theta)
$$
而
Q_i(z_i,\theta)=P(Z=z_i|Y=y_i, \theta)=\frac{P(Z=z_i|\theta)P(Y=y_i|Z=z_i, \theta)}{P(Y=y_i|\theta)}\\
=\frac{(1-z_i)\pi{10\choose y_i}\ p_A^{y_i}(1-p_A)^{10-y_i}+z_i(1-\pi){10\choose y_i}\ p_B^{y_i}(1-p_B)^{10-y_i}}{\pi{10\choose y_i}\ p_A^{y_i}(1-p_A)^{10-y_i}+(1-\pi){10\choose y_i}\ p_B^{y_i}(1-p_B)^{10-y_i}}
$$
可以看出：$1-Q_i(0)=Q_i(1)$
在 E 步，我们通过统计数据计算出每次实验的 $Q_i(z)$，然后固定参数 $\theta$，改变参数 $\theta’=(\pi’,p_A’,p_B’)$，最大化下面的函数：
Q(\theta’;\theta)=\sum_{i=1}^{n}\sum_{z_i \in {0,1}}Q_i(z_i,\theta)\ln \frac{P(Y=y_i,Z=z_i|\theta’)}{Q_i(z_i,\theta)}\\
=\sum_{i=1}^{n}Q_i(0,\theta)\ln \frac{\pi'{10\choose y_i}\ p_A’^{y_i}(1-p_A’)^{10-y_i}}{Q_i(0,\theta)} + \sum_{i=1}^{n}Q_i(1,\theta)\ln \frac{(1-\pi’){10\choose y_i}\ p_B’^{y_i}(1-p_B’)^{10-y_i}}{Q_i(1,\theta)}
$$
对 $\pi$ 求偏导得：
\frac{\partial Q(\theta’; \theta)}{\partial\pi’}=\frac{\sum_{i=1}^nQ_i(0,\theta)}{\pi’}-\frac{\sum_{i=1}^nQ_i(1,\theta)}{1-\pi’}
$$
令导数为 0 得：
\pi’=\frac{\sum_{i=1}^{n}Q_i(0,\theta)}{n}
$$
同理可得：
p_A’=\frac{\sum_{i=1}^{n}y_iQ_i(0,\theta)}{\sum_{i=1}^{n}10Q_i(0,\theta)}\\
p_B’=\frac{\sum_{i=1}^{n}y_iQ_i(1,\theta)}{\sum_{i=1}^{n}10Q_i(1,\theta)}
$$
这就是 M 步的更新公式。
EM 算法需要我们手动指定一个初值，然后开始迭代，迭代的终止条件为 $||\theta’-\theta|| < \epsilon_1$ 或 $||Q(\theta’; \theta)-Q(\theta; \theta)|| < \epsilon_2$
指定初值 $\pi=0.5$，$p_A=0.6$，$p_B=0.5$ 的情况下，最终迭代过程如下：
Iteration 1 pi = 0.597 p_A = 0.713, p_B = 0.581
Iteration 2 pi = 0.591 p_A = 0.733, p_B = 0.555
Iteration 3 pi = 0.582 p_A = 0.752, p_B = 0.532
Iteration 4 pi = 0.572 p_A = 0.767, p_B = 0.516
Iteration 5 pi = 0.564 p_A = 0.777, p_B = 0.509
Iteration 6 pi = 0.556 p_A = 0.783, p_B = 0.506
Iteration 7 pi = 0.550 p_A = 0.786, p_B = 0.506
Iteration 8 pi = 0.545 p_A = 0.788, p_B = 0.507
Iteration 9 pi = 0.541 p_A = 0.789, p_B = 0.508
Iteration 10 pi = 0.538 p_A = 0.790, p_B = 0.509
Iteration 11 pi = 0.535 p_A = 0.791, p_B = 0.510
Iteration 12 pi = 0.533 p_A = 0.791, p_B = 0.510
Iteration 13 pi = 0.531 p_A = 0.792, p_B = 0.511
Iteration 14 pi = 0.529 p_A = 0.792, p_B = 0.512
Iteration 15 pi = 0.528 p_A = 0.792, p_B = 0.512
Iteration 16 pi = 0.527 p_A = 0.792, p_B = 0.512
no update in params, stop iteration
Result pi = 0.527 p_A = 0.792, p_B = 0.512