讲一个EM算法和Gibbs 抽样的小例子，用于加深理解（变分推断版本请见变分推断学习笔记(3)——三硬币问题的变分推断解法）。
题目(引用自参考1）：假设有3枚硬币，分别记做A，B，C。这些硬币正面出现的概率分别是\(\pi\),\(p\)和\(q\)。进行如下掷硬币实验：先掷硬币A，根据其结果选出硬币B或C，正面选B，反面选硬币C；然后投掷选重中的硬币，出现正面记作1，反面记作0；独立地重复\(n\)次（n=10)，结果为 \[1111110000\] 我们只能观察投掷硬币的结果，而不知其过程，估计这三个参数\(\pi\),\(p\)和\(q\)。
EM算法
可以看到投掷硬币时到底选择了B或者C是未知的。我们设隐藏变量Z 来指示来自于哪个硬币，\(Z=\{z_1,z_2,\ldots,z_n \}\)，令\(\theta=\{\pi,p,q\}\)，观察数据\(X=\{x_1,x_2,\ldots,x_n \}\)。
写出生成一个硬币时的概率： \[\begin{split}P(x|\theta) & =\sum_z P(x,z|\theta)=\sum_z P(z|\pi)P(x|z,\theta) \\& =\pi p^x (1-p)^{1-x}+(1-\pi)q^x(1-q)^{1-x} \\\end{split}\] 有了一个硬币的概率，我们就可以写出所有观察数据的log似然函数： \[L(\theta|X)=\log P(X|\theta)=\sum^n_{j=1}\log[\pi p^{x_j} (1-p)^{1-{x_j}}+(1-\pi)q^{x_j}(1-q)^{1-{x_j}}]\] 然后求极大似然 \[\hat{\theta}=\arg \max L(\theta|X)\] 其中\(L(\theta|X)=\log P(X|\theta)=\log \sum_Z P(X,Z|\theta)\)。因为log里面带着加和所以这个极大似然是求不出解析解的。 如果我们知道隐藏变量Z的话，求以下的似然会容易很多: \[L(\theta|X,Z)=\log P(X,Z|\theta)\] 但是隐藏变量Z的值谁也不知道，所以我们转用EM算法求它的后验概率期望\(E_{Z|X,\theta}(L(\theta|X,Z))\)的最大值。
E步：假设当前模型的参数为\(\pi,p,q\)时，隐含变量来自于硬币B的后验概率 \[\mu=P(Z|X,\theta)=\frac{\pi p^{x_i}(1-p)^{1-x_i}}{ \pi p^{x_i}(1-p)^{1-x_i}+(1-\pi)q^{x_i}(1-q)^{1-x_i}}\] 那么隐含变量来自于硬币C的后验概率自然为\(1-\mu\)。
M步： 先写出似然关于后验概率的期望，它是似然期望下界函数的最大值， \[\begin{split}&E_{Z|X,\theta}(L(\theta|X,Z))=\\&\sum^n_{j=1}[ \mu \log{\frac{\pi p^{x_i}(1-p)^{1-x_i}}{\mu}}+(1-\mu) \log{\frac{(1-\pi)q^{x_i}(1-q)^{1-x_i}}{1-\mu}} ]\\\end{split}\] 要注意这里把\(\mu\)看做固定值。然后我们分别求偏导，获得参数\(\pi,p,q\)的新估计值 \[\begin{split}& \pi=\frac{1}{n}\sum^n_{j=1}\mu_j \\& p=\frac{\sum^n_{j=1}\mu_j x_j}{\sum^n_{j=1}\mu_j} \\& q=\frac{\sum^n_{j=1}(1-\mu_j) x_j}{\sum^n_{j=1}(1-\mu_j)} \\\end{split}\]
算法首先选取参数初始值\(\theta^{(0)}=\pi^{(0)},p^{(0)},q^{(0)}\)，然后迭代到收敛为止。
Gibbs Sampling
其实这道题用EM做就可以了，因为最近看了Gibbs Sampling，所以套一下看看，如果有错敬请指出~
先给出Gibbs sampling的抽象的应用步骤，由于目标是对latent varibles逐一sample得到complete data再算充分统计量。所以关键是求 \(P(Z_i|Z_{-i},O,\alpha)\)：
- 根据Markov bucket性质写出latent variable与对应observed variable的joint分布: \(P(Z,O|\alpha)=f(Z,O,\theta,\alpha)\)
- 根据样本对立性和参数条件独立性对joint distribution进行分解
- 积分参数\(\theta\)，得到仅包含超参和充分统计量的joint分布表达式:\(P(Z,O|\alpha)=f(Z,O,\alpha)\)
- 对单个latent variable, 计算full conditional概率，用来采样: \(P(Z_{i}|Z_{-i},O,\alpha)\)
建模过程
首先我们确定建模过程，画出概率图： 
可以看到硬币B还是硬币C的类标签服从贝努利分布 \[Z_j \sim Bernoulli(\pi)\] 然后硬币正反面的出现也服从贝努利分布 \[W_{z_{j}k} \sim Bernoulli(\theta_{Z_j})\] 其中\(\theta=\{p,q\}\)。需要注意这里参数都变成了随机变量，都有先验分布。为了简单和无偏倚起见，Beta分布的先验参数均设为1, \(\alpha=<1,1>\)，\(\beta=<1,1>\)，这样\(p(\alpha)=1,p(\beta)=1\)，在实际运算中就不考虑了。
确定联合分布
然后写出联合分布来应该是（要注意这里的联合分布应该是似然与先验\(p(\theta)\)的乘积，因为先验\(p(\alpha)=1,p(\beta)=1\)所以省略了后项)
\[ \begin{split}P(C,Z,\pi,p,q)&=p(W|Z,\Theta)p(\Theta|\gamma_{\theta})p(Z|\pi)p(\pi|\gamma_{\pi})\\ &= \prod^N_{i=1} \pi^{z_{i}} (1-\pi)^{1-z_{i}} p^{W_{z_{i0}}} (1-p)^{W_{z_{i1}}} q^{W_{(1-z_{i})0}} (1-q)^{W_{(1-z_{i})1}} \\ & =\pi^{C_0} (1-\pi)^{C_1} p^{N_{C_0}(0)} (1-p)^{N_{C_0}(1)} q^{N_{C_1}(0)} (1-q)^{N_{C_1}(1)} \\ \end{split} \]
其中\(z_i=0\)时选中硬币B，\(z_i=1\)时选中硬币C，\(W_{xj}\)是观察到的硬币B或硬币C的正反结果，正面为\(W_{x0}=1,W_{x1}=0\)，反面为\(W_{x1}=1,W_{x0}=0\) ，\(C_0\)是所有结果中选中硬币B的次数，\(C_1\)是选中硬币C的次数，\(N_{C_0}(0)\)是选B 的硬币中正面的硬币个数，\(N_{C_0}(1)\) 是选B的硬币中反面的硬币个数，\(N_{C_1}(0)\)是选C的硬币中正面的硬币个数，\(N_{C_1}(1)\)是选C的硬币中反面的硬币个数。
Standard Gibbs Sampling
对类标签Z采样
然后根据Gibbs Sampling对隐藏变量Z采样，首先，我们设\({Z}^{-j}\)是除了\(Z_j\)外所有的硬币标签，然后\(C^{(-j)}\)是除了\(W_j\)外所有的硬币结果集合。硬币A的投掷有2个结果，\(Z_j=0\)（硬币B）或\(Z_j=1\)（硬币C）。如果\(Z_j=0\)，那么\(C_0^{-j}=C_0-1,C_1^{-j}=C_1\)。反之如果文档\(Z_j=1\)，那么\(C_1^{-j}=C_1-1,C_0^{-j}=C_0\)。
我们构建Gibbs Sampler，采样公式为： \[\begin{split}P(Z_j=0|Z^{-j},C^{-j},\pi,p,q) &=\frac{P(Z,C,\pi,p,q)}{P(Z^{-j},C^{-j},\pi,p,q)} \\&=\pi p^{W_{j0}} (1-p)^{W_{j1}}\\P(Z_j=1|Z^{-j},C^{-j},\pi,p,q) &=\frac{P(Z,C,\pi,p,q)}{P(Z^{-j},C^{-j},\pi,p,q)} \\&=(1-\pi) q^{W_{j0}} (1-q)^{W_{j1}}\end{split}\]
然后根据这两个概率对\(Z_j^{(t+1)}\)做贝努利实验，得到新采样的结果。
对剩下参数的采样
在对\(Z\)采样获得了所有硬币的正反面之后，我们接下来要估计\(\pi,p,q\)。
Gibbs Sampler的公式分别为: \[\begin{split}&P(\pi|Z,C,p,q)=\pi^{C_0} (1-\pi)^{C_1}\\&P(p|Z,C,\pi,q)=p^{N_{C_0}(0)} (1-p)^{N_{C_0}(1)}\\&P(q|Z,C,\pi,p)=q^{N_{C_1}(0)} (1-q)^{N_{C_1}(1)}\\ \end{split}\] 这里\(\pi,p,q\)均服从贝努利分布，根据Beta-Bernoulli共轭(之前先验分布为Beta(1,1))，我们可以把他们看成后验的Beta分布 \[\begin{split}&\pi \sim Beta(1+C_0,1+C_1)\\&p \sim Beta(1+N_{C_0}(0),1+N_{C_0}(1))\\&q \sim Beta(1+N_{C_1}(0),1+N_{C_1}(1))\\\end{split}\] 然后从各自的Beta分布中采样出\(\pi,p,q\)的值，详见这里的求\(\theta\)部分。
Collapsed Gibbs Sampling
对变量积分
根据《Gibbs Sampling for the UniniTiated》上所说,\(\pi\)可以被积分掉来简化联合分布，这种方法被称作“Collapsed Gibbs Sampling”,指通过积分避开了实际待估计的参数，转而对每个硬币的正反面\(z\)进行采样，一旦每个硬币的z确定下来，被积分的值可以在统计频次后计算出来。正巧这里\(\theta=\{p,q\}\)也只生成一个样本，可以被积分，就统统给积分好了。 \[\begin{split}&P(C,Z)= \int\int\int P(C,Z,\pi,p,q) d \pi d p d q \\&=\int\int\int \pi^{C_0} (1-\pi)^{C_1} p^{N_{C_0}(0)} (1-p)^{N_{C_0}(1)} q^{N_{C_1}(0)} (1-q)^{N_{C_1}(1)} d \pi d p d q \\&=\int \pi^{C_0} (1-\pi)^{C_1}d \pi \int p^{N_{C_0}(0)} (1-p)^{N_{C_0}(1)} dp \int q^{N_{C_1}(0)} (1-q)^{N_{C_1}(1)} \\&=\frac{\Gamma(C_0+1)\Gamma(C_1+1)}{\Gamma(C_0+C_1+2)}\frac{\Gamma(N_{C_0}(0)+1)\Gamma(N_{C_0}(1)+1)}{\Gamma(C_0+2)}\frac{\Gamma(N_{C_1}(0)+1)\Gamma(N_{C_1}(1)+1)}{\Gamma(C_1+2)}\end{split}\] 最后一步是根据Beta分布的积分公式得出的。要注意这个式子里还是存在未知数的\(Z\)的，只不过被加和掩盖了而已。
对类标签Z采样
在这里采样公式就变成了 \[\begin{split}&P(Z_j=0|Z^{-j},C^{-j})=\frac{C_0}{C_0+C_1+1}(\frac{N_{C_0}(0)}{C_0+1})^{W_{j0}}(\frac{N_{C_1}(0)}{C_1+1})^{W_{j1}}\\&P(Z_j=1|Z^{-j},C^{-j})=\frac{C_1}{C_0+C_1+1}(\frac{N_{C_0}(1)}{C_0+1})^{W_{j0}}(\frac{N_{C_1}(1)}{C_1+1})^{W_{j1}}\\ \end{split}\] 根据这两个概率对\(Z_j^{(t+1)}\)做贝努利实验，得到新采样的结果。
在\(Z\)的一轮迭代结束后我们还要估计被积分了的变量的值，已知后验概率为 \[\begin{split}&p(\pi|Z)=\prod^n_{i=1}p(z|\pi)p(\pi|r_\pi)=Beta(\pi|1+C_0,1+C_1)\\&p(p|Z)=\prod^n_{i=1}p(z|p)p(p|r_\theta)=Beta(\pi|1+N_{C_0}(0),1+N_{C_0}(1))\\&p(q|Z)=\prod^n_{i=1}p(z|q)p(q|r_\theta)=Beta(\pi|1+N_{C_1}(0),1+N_{C_1}(1))\\ \end{split}\] 根据贝叶斯推断的话，我们应该求整个参数分布上的期望，所以有 \[E(\pi)=\int \pi \cdot \pi^{C_0+1-1}(1-\pi)^{C_1+1-1} d \pi=\frac{C_0+1}{C_0+C_1+2}\] 这里还是利用Beta分布的积分公式，同理剩下的俩为 \[\begin{split}&E(p)=\frac{N_{C_0}(0)+1}{C_0+2}\\&E(q)=\frac{N_{C_1}(0)+1}{C_1+2}\\ \end{split}\] 这里可以看到先验知识（分子式上面的1和分母的2）和观察到的数据很好的结合了起来。
EM算法和Sampling算法（见PRML的11.1.6和参考3)
两者的联系
两者都常用来估计latent variable分布的参数。主要思想一致：先估计生成一些latent variable，然后看成complete data，算参数。 区别：
EM直接估计所有latent varibles的联合概率;
Gibbs估计单个单个latent variable的条件概率;
所以，EM类似梯度下降，而Gibbs类似于coordinate gradient decent(不是conjugate)。
Monte Carlo EM算法
EM算法中的E步估计所有latent varible的联合概率分布，可能很复杂。所以可以用抽样的方法来估计，以下是似然函数的期望表达式： \[Q(\theta,\theta^{old})=\int p(Z|X,\theta^{old})\ln p(Z,X|\theta)d Z\] 我们可以在后验概率\(p(Z|X,\theta^{old})\)的分布上进行有限次数的抽样来接近这个期望： \[Q(\theta,\theta^{old}) \simeq \frac{1}{L}\sum^L_{l=1}\ln p(Z^{(l)},X|\theta)\] 然后像往常一样在M步优化Q函数，这个方法叫做Monte Carlo EM算法。
而每次在E步只生成一个样本的方法叫做stochastic EM，是Monte Carlo EM 算法的一个特例。在E步，我们从后验概率\(p(Z|X,\theta^{old})\)中抽取一个样本，用来近似计算似然的期望，然后在M步更新模型参数的值。
刚才介绍的Collapsed Gibbs Sampling的方法就是stochastic EM的一个例子，只不过，M步骤的参数估计结果在E步骤中没有用到，所以不需要重复多余的M步骤，只需在最后进行一次M步骤，得到所需要的参数即可。
参考文献：
- 《统计学习方法》
- 《Pattern Recognition and Machine Learning》
- http://xhyan.wordpress.com/2012/04/30/%E3%80%90todo%E3%80%91quick-derivation-in-gibbs-sampling/