机器学习 | EM 算法

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