EM算法解决三硬币问题

EM算法解决三硬币问题

1.问题介绍

假如有三个硬币，分别记做A,B,C，这些硬币正面出现的概率分别是$\pi$,$p$,$q$。进行如下掷硬币实验：先掷硬币A，根据其结果选出硬币B或硬币C，正面选择硬币B，反面选择硬币C。然后掷选出的硬币，掷硬币的结果，出现正面记作1，反面记作0，独立重复n次实验。

2.EM算法简介

输入观测数据$Y$，隐变量数据$Z$，联合分布$P(Y,Z|\theta )$，条件分布$P(Z|Y,\theta )$，输出模型参数$\theta$.
（1）选择参数的初值${\theta }^{(0)}$，开始迭代
（2）$E$步：记${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计值，在第$i+1$次迭代的$E$步，计算
$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]\\ & =\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})\end{array}$$
这里，$P(Z|Y,{\theta }^{(i)})$是在给定观测数据$Y$和当前的参数估计${\theta }^{(i)}$下隐变量数据$Z$的条件概率分布
（3）$M$步：求使$Q(\theta ,{\theta }^{(i)})$极大化的$\theta$，确定第$i+1$次迭代的参数估计值${\theta }^{(i+1)}$
$${\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^{(i)})$$
（4）重复（2）和（3）直到收敛。

3.EM算法解决三硬币问题

（1）首先选取参数的初值，记作${\theta }^{(0)}=({\pi }^{(0)},p^{(0)},q^{(0)})$
（2）$E$步，计算在模型参数${\pi }^{(i)},p^{(i)},q^{(i)}$下观测数据$y_j$来自掷硬币B的概率
$${\mu }_j^{(i+1)}=\frac{{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}}{{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)})+(1-{\pi }^{(i)})(q^{(i)}{)}^{y_j}(1-q^{(i)}{)}^{1-y_j}}$$
（3）$M$步，计算模型参数的新估计值
$${\pi }^{(i+1)}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}$$
$$p^{(i+1)}=\frac{\sum _{j=1}^n{\mu }_j^{(i+1)}{y}_j}{\sum _{j=1}^n{\mu }_j^{i+1}}$$
$$q^{i+1}=\frac{\sum _{j=1}^n(1-{\mu }_j^{(i+1)})y_j}{\sum _{j=1}^n(1-{\mu }_j^{(i+1)})}$$
（4）重复（2）（3）步的迭代直到收敛

以上内容均参考李航老师的《统计学习方法》一书

4.推导

《统计学习方法》一书略去了三硬币问题的推导。然而，给出的解决算法中E步和M步所计算的内容与EM算法形式上的E步和M步的计算内容有些出入，因此，有些地方需要更清楚的解释。
这个问题中，我们设定的隐变量是什么？我认为，隐变量$Z$应该代表了A硬币是正面或反面，正面是1,反面是0。
为了便于后边的论述，我们假设$\theta =(\pi ,p,q)$
首先理解E步。为了理解E步的内容，我们从$Q$函数入手。注意到：
$$Q(\theta ,{\theta }^{(i)})=\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$
观察EM算法$E$步的公式，我们首先需要得到观测数据与隐变量在$\theta$条件下的联合分布$P(Y,Z|\theta )$。如何求联合概率分布?基于以下公式
$$P(Y,Z|\theta )=P(Z|\theta )P(Y|Z,\theta )$$
假设我们抽取的n个样本相互独立，显然我们每次掷硬币A都是独立的，即$z_j$之间相互独立。每次观测的数据$y_j$的值也仅仅与$z_j$有关，于是我们可以得到
$$P(Z|\theta )=\prod _{j=1}^{n}P(z_j|\theta )=\prod _{j=1}^n{\pi }^{z_j}(1-\pi {)}^{1-z_j}$$
其中
$$P(z_j|\theta )={\pi }^{z_j}(1-\pi {)}^{1-z_j}$$
这个是由于$z$服从二项分布。$P(z_j=0|\theta )=1-\pi$，$P(z_j=1|\theta )=\pi$。
如何求出$P(Y|Z,\theta )$?当$z_j$已经确定,我们也就知道投掷哪枚硬币了。这样我们根据$z_j$的取值可以写出以下的内容。$P(y_j=1|z_j=1,\theta )=p$，$P(y_j=1|z_j=0,\theta )=q,P(y_j=0|z_j=1,\theta )=1-p$，$P(y_j=0|z_j=0,\theta )=1-q$。对于仅有的这四种情况，我们可以用一个公式包含，即
$$P(y_j|z_j,\theta )=(p^{z_j}{q}^{(1-z_j)}{)}^{y_j}((1-p{)}^{z_j}(1-q{)}^{(1-z_j)}{)}^{(1-y_j)}$$
因此，由独立性假设，我们可以得到
$$P(Y|Z,\theta )=\prod _{j=1}^n(p^{z_j}{q}^{(1-z_j)}{)}^{y_j}((1-p{)}^{z_j}(1-q{)}^{(1-z_j)}{)}^{(1-y_j)}$$
所以,由$P(Y,Z|\theta )=P(Z|\theta )P(Y|Z,\theta )$将$P(Z|\theta )$和$P(Y|Z,\theta )$相乘，整理，得到
$$P(Y,Z|\theta )=\prod _{j=1}^n[{\pi }^{z_j}{p}^{z_j{y}_j}(1-p{)}^{z_j(1-y_j)}][(1-\pi {)}^{1-z_j}{q}^{(1-z_j)y_j}(1-q{)}^{(1-z_j)(1-y_j)}]$$
对上边的公式对$Z$求和，也就是对每一个乘积项将$z_j=0$和$z_j=1$的情况加起来，我们可以得到
$$P(Y|\theta )=\sum _{Z}P(Y,Z|\theta )=\prod _{j=1}^n[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}]$$
由于有了$P(Y,Z|\theta )$和$P(Y|\theta )$,我们可以得到
$$P(Z|Y,\theta )=\frac{P(Y,Z|\theta )}{P(Y|\theta )}=\prod _{j=1}^n\frac{[{\pi }^{z_j}{p}^{z_j{y}_j}(1-p{)}^{z_j(1-y_j)}][(1-\pi {)}^{1-z_j}{q}^{(1-z_j)y_j}(1-q{)}^{(1-z_j)(1-y_j)}]}{[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}]}$$
按照EM算法中的需求带入${\theta }^{(i)}$
$$P(Z|Y,{\theta }^{(i)})=\prod _{j=1}^n\frac{[({\pi }^{(i)}{)}^{z_j}(p^{(i)}{)}^{z_j{y}_j}(1-p^{(i)}{)}^{z_j(1-y_j)}][(1-{\pi }^{(i)}{)}^{1-z_j}(q^{(i)}{)}^{(1-z_j)y_j}(1-q^{(i)}{)}^{(1-z_j)(1-y_j)}]}{[{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}+(1-{\pi }^{(i)})(q^{(i)}{)}^{y_j}(1-q^{(i)}{)}^{1-y_j}]}$$
至此，$Q$函数中所需要的组成部分都计算出来了。我们也就完成了$E$步
李航老师书中的${\mu }_j^{(i+1)}$是什么？我们观察上面的公式，可以发现
$${\mu }_j^{(i+1)}=P(z_j=1|y_j,{\theta }^{(i)})$$$$1-{\mu }_j^{(i+1)}=P(z_j=0|y_j,{\theta }^{(i)})$$
因为使用${\mu }_j^{(i+1)}$这个符号,会使得公式变得简单，后边我也会用这个符号说明问题。
接下类我们来理解$M$步。
$$\begin{array}{rl}{\theta }^{(i+1)}& =\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^{(i)})\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{Z}P(Z|Y,{\theta }^{(i)})\sum _{j=1}^{n}logP(y_j,z_j|\theta )\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{Z}P(Z|Y,{\theta }^{(i)})(logP(y_1,z_1|\theta )+\sum _{j=2}^{n}logP(y_j,z_j|\theta ))\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(y_1,z_1|\theta )+\sum _{Z}P(Z|Y,{\theta }^{(i)})\sum _{j=2}^{n}logP(y_j,z_j|\theta ))\\ & =\underset{\theta }{\mathrm{argmax}}(P(z_1=0|Y,{\theta }^{(i)})logP(y_1,z_1=0|\theta )+P(z_1=1|Y,{\theta }^{(i)})logP(y_1,z_1=1|\theta ))\sum _{Z^{\prime }}P(Z^{\prime }|Y^{\prime },\theta )+\sum _{Z}P(Z|Y,{\theta }^{(i)})\sum _{j=2}^{n}logP(y_j,z_j|\theta ))\\ & =\underset{\theta }{\mathrm{argmax}}(P(z_1=0|Y,{\theta }^{(i)})logP(y_1,z_1=0|\theta )+P(z_1=1|Y,{\theta }^{(i)})logP(y_1,z_1=1|\theta ))+\sum _{Z}P(Z|Y,{\theta }^{(i)})\sum _{j=2}^{n}logP(y_j,z_j|\theta ))\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{z_1}(P(z_1|Y,{\theta }^{(i)})logP(y_1|\theta ))+\sum _{Z}P(Z|Y,{\theta }^{(i)})\sum _{j=2}^{n}logP(y_j,z_j|\theta ))\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{z_1}(P(z_1|Y,{\theta }^{(i)})logP(y_1|\theta ))+\sum _{z_2}(P(z_2|Y,{\theta }^{(i)})logP(y_2|\theta ))+\sum _{Z}P(Z|Y,{\theta }^{(i)})\sum _{j=3}^{n}logP(y_j,z_j|\theta ))\\ & =...\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{j=1}^n\sum _{z_j}P(z_j|y_j,{\theta }^{(i)})log(y_j,z_j|\theta )\end{array}$$
其中$Z^{\prime }=(z_2,z_3,...,z_n),Y^{\prime }=(y_2,y_3,...y_n)$。显然，$\sum _{Z^{\prime }}P(Z^{\prime }|Y^{\prime },\theta )=1$，($Z^{\prime }$的条件概率对$Z^{\prime }$的积分一定是1)
（注：参考了 https://blog.csdn.net/zsdust/article/details/100042491）
带入我们上边得到的公式，得到
θ ( i + 1 ) = a r g m a x θ ∑ j = 1 n ∑ z j P ( z j ∣ y j , θ ( i ) ) l o g ( y j , z j ∣ θ ) = a r g m a x θ ∑ j = 1 n { [ l o g π p y j ( 1 − p ) ( 1 − y j ) ] μ j ( i + 1 ) + [ l o g ( 1 − π ) q y j ( 1 − q ) ( 1 − y j ) ] ( 1 − μ j ( i + 1 ) ) } \begin{aligned} \theta^{(i+1)} &= \mathop{argmax}\limits_{\theta}\sum_{j=1}^n\sum_{z_j}P(z_j|y_j,\theta^{(i)})log(y_j,z_j|\theta)\\ &=\mathop{argmax}\limits_{\theta}\sum_{j=1}^n\{[log\pi p^{y_j}(1-p)^{(1-y_j)}]\mu_j^{(i+1)} + [log(1-\pi)q^{y_j}(1-q)^{(1-y_j)}](1-\mu_j^{(i+1)})\} \end{aligned} θ(i+1)​=θargmax​j=1∑n​zj​∑​P(zj​∣yj​,θ(i))log(yj​,zj​∣θ)=θargmax​j=1∑n​{[logπpyj​(1−p)(1−yj​)]μj(i+1)​+[log(1−π)qyj​(1−q)(1−yj​)](1−μj(i+1)​)}​
上式对$\pi$,$p$,$q$求导数，分别让导数等于0即可得到上边的公式。
对$\pi$求导
$$\begin{array}{rl}\frac{1}{{\pi }^{(i+1)}}\sum _{j=1}^n{\mu }_j^{(i+1)}+\frac{1}{1-{\pi }^{(i+1)}}(1-{\mu }_j^{(i+1)})& =0\\ (1-{\pi }^{(i+1)})\sum _{j=1}^n{\mu }_j^{(i+1)}+{\pi }^{(i+1)}\sum _{j=1}^n(1-{\mu }_j^{(i+1)})& =0\\ {\pi }^{(i+1)}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}& \end{array}$$
对$p$求导
$$\begin{array}{rl}\sum _{j=1}^n[y_j(1-p^{(i+1)})+(y_j-1)p^{(i+1)}]{\mu }_j^{(i+1)}& =0\\ \sum _{j=1}^n[y_j-p^{(i+1)}]{\mu }_j^{(i+1)}& =0\\ \sum _{j=1}^n{y}_j{\mu }_j^{(i+1)}-p^{(i+1)}\sum _{j=1}^n{\mu }_j^{(i+1)}& =0\\ p^{(i+1)}=\frac{\sum _{j=1}^n{\mu }_j^{(i+1)}{y}_j}{\sum _{j=1}^n{\mu }_j^{i+1}}& \end{array}$$
对$q$求导
$$\begin{array}{rl}\sum _{j=1}^n{y}_{j}log{q}^{(i+1)}+(1-y_j)log(1-q^{(i+1)})(1-{\mu }_j^{(i+1)})& =0\\ \sum _{j=1}^n{y}_j(1-{\mu }_j^{(i+1)})-q\sum _{j=1}^n(1-{\mu }_j^{(i+1)})& =0\\ q^{i+1}=\frac{\sum _{j=1}^n(1-{\mu }_j^{(i+1)})y_j}{\sum _{j=1}^n(1-{\mu }_j^{(i+1)})}& \end{array}$$