EM 算法实验

Question

最大似然估计

在已知每次实验的硬币是 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