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