贝叶斯统计学习笔记|Bayesian Statistics|Metropolis-Hastings与Gibbs Sampling

（一） Metropolis-Hastings(MH)

现要从目标分布p(θ)∝g(θ)p(\theta)\propto g(\theta)p(θ)∝g(θ)中抽样，MH算法的思想是：构建一个马氏链，其平稳分布就是目标分布。在构建时，我们首先随机选取一个初始值，从另一个更容易抽样的分布中将抽取出的样本作为候选，进行接受或拒绝操作，所得的马氏链收敛于目标分布，如此也就得到了后验分布的样本（posterior samples)。

MH算法的过程如下：

选取初始值θ0\theta_0θ0
对i=1,...,mi=1,...,mi=1,...,m，重复如下操作：
（1）从proposal distribution q(θ∗∣θi−1)q(\theta^{*}| \theta_{i-1})q(θ∗∣θi−1)中抽取一个候选样本θ∗\theta^{*}θ∗
（2）计算接受率α\alphaα：
α=g(θ∗)/q(θ∗∣θi−1)g(θi−1)/q(θi−1∣θ∗)\alpha = \frac{g(\theta^{*}) /q(\theta^{*}| \theta_{i-1})}{g(\theta_{i-1}) /q(\theta_{i-1}| \theta^{*})}α=g(θi−1)/q(θi−1∣θ∗)g(θ∗)/q(θ∗∣θi−1)
（3）接受或拒绝候选样本（规则如下）：
\quad若α≥1\alpha \geq1α≥1，则令θi=θ∗\theta_{i}=\theta^{*}θi=θ∗（即接受候选样本）；
\quad若α<1\alpha <1α<1，则
\qquad θi=θ∗\theta_{i}=\theta^{*}θi=θ∗ with prob α\alphaα
\qquad θi=θi−1\theta_{i}=\theta_{i-1}θi=θi−1 with prob 1−α1-\alpha1−α

上述过程中的proposal distribution是人为选择、更方便抽样的分布，由于分布q不同于目标分布p，因此MH算法中的接受/拒绝操作可以看做是对分布q的一种修正。对于马氏链的每一步，我们或者拒绝、或者接受候选样本，其规则取决于接受率α\alphaα的大小。如果我们跳过这种筛选机制，每一步都接受候选样本，那么显而易见，所得到的最终结果并非是对目标分布p而是分布q的蒙特卡洛模拟。

例子：已知一枚不均匀的硬币出现正面概率为0.7。现掷另一枚未知均匀的硬币共5次，2次正面，3次反面。求硬币不均匀的概率P(loaded∣X=2)P(loaded|X=2)P(loaded∣X=2)？

方法一：贝叶斯公式
θ={fair,loaded}\theta=\{fair, loaded\}θ={fair,loaded}其中fair代表硬币均匀，loaded代表硬币不均匀
由于5次中出现3次反面，因此可以设先验概率P(θ\thetaθ=loaded)=0.6
似然：f(x∣θ)=(5x)125I{θ=fair}+(5x)0.7x×0.35−xI{θ=loaded}f(x|\theta)=\binom{5}{x}\frac{1}{2^5}I_{\{\theta=fair\}}+\binom{5}{x}0.7^{x}\times 0.3^{5-x}I_{\{\theta=loaded\}}f(x∣θ)=(x5)251I{θ=fair}+(x5)0.7x×0.35−xI{θ=loaded}

从而：
f(θ∣X=2)=f(x∣θ)f(θ)f(x)=125×0.4×I{θ=fair}+0.72×0.33×0.6I{θ=loaded}125×0.4+0.72×0.33×0.6=0.0125I{θ=fair}+0.0794I{θ=loaded}0.0125+0.0794=0.612I{θ=fair}+0.388I{θ=loaded}\begin{aligned} f(\theta|X=2)&=\frac{f(x|\theta)f(\theta)}{f(x)}\\ &=\frac{\frac{1}{2^5}\times 0.4\times I_{\{\theta=fair\}}+0.7^2\times 0.3^3\times 0.6I_{\{\theta=loaded\}}}{\frac{1}{2^5}\times 0.4+0.7^2\times 0.3^3\times 0.6}\\ &=\frac{0.0125I_{\{\theta=fair\}}+0.0794I_{\{\theta=loaded\}}}{0.0125+0.0794} \\ &=0.612 I_{\{\theta=fair\}}+0.388I_{\{\theta=loaded\}} \end{aligned} f(θ∣X=2)=f(x)f(x∣θ)f(θ)=251×0.4+0.72×0.33×0.6251×0.4×I{θ=fair}+0.72×0.33×0.6I{θ=loaded}=0.0125+0.07940.0125I{θ=fair}+0.0794I{θ=loaded}=0.612I{θ=fair}+0.388I{θ=loaded}
因此P(loaded∣X=2)=0.388P(loaded|X=2)=0.388P(loaded∣X=2)=0.388

方法二：MH算法

设定初始值θ0=fair\theta_{0}=fairθ0=fair或θ0=loaded\theta_{0}=loadedθ0=loaded
对i=1,...,mi=1,...,mi=1,...,m，重复如下操作：
（1）令候选样本θ∗\theta^{*}θ∗为不同于θi−1\theta_{i-1}θi−1的状态（即，若θi−1\theta_{i-1}θi−1为正面，则θ∗\theta^{*}θ∗就为反面，反之亦然）
（2）由候选样本的选取知q(θ∗∣θi−1)=q(θi−1∣θ∗)=1q(\theta^{*}| \theta_{i-1})=q(\theta_{i-1}| \theta^{*})=1q(θ∗∣θi−1)=q(θi−1∣θ∗)=1（因为二者互斥，已知其中一个，另一个的状态能完全确定）
\qquad 故接受率α=g(θ∗)g(θn−1)=f(x=2∣θ∗)f(θ∗)f(x=2∣θi−1)f(θi−1)\alpha = \frac{g(\theta^{*})}{g(\theta_{n-1})}=\frac{f(x=2|\theta^{*})f(\theta^{*})}{f(x=2|\theta_{i-1})f(\theta_{i-1})}α=g(θn−1)g(θ∗)=f(x=2∣θi−1)f(θi−1)f(x=2∣θ∗)f(θ∗)

因此，当 θ∗=fair\theta^{*}= fairθ∗=fair 时，α=0.1250.0794=1.574\alpha = \frac{0.125}{0.0794}=1.574α=0.07940.125=1.574。此时应接受θ∗\theta^{*}θ∗，即令θi=fair\theta_{i}=fairθi=fair;

当 θ∗=loaded\theta^{*}= loadedθ∗=loaded 时，α=0.07940.0125=0.635\alpha = \frac{0.0794}{0.0125}=0.635α=0.01250.0794=0.635。此时θi=loaded\theta_{i}=loadedθi=loaded （with prob. 0.635），θi=fair\theta_{i}=fairθi=fair （with prob. 0.365）

得到马氏链及其转移矩阵：P=[0.3650.63510]P = \begin{bmatrix} 0.365 & 0.635\\ 1 & 0 \end{bmatrix} P=[0.36510.6350]

由于马氏链的平稳分布π\piπ满足：πP=π\pi P=\piπP=π
可求得π=[0.6120.388]\pi = [0.612\quad 0.388]π=[0.6120.388]
从而待求的后验概率为0.388。

（二） Gibbs Sampling

Gibbs采样可以方便我们求多个参数的后验分布。现有两个参数θ\thetaθ和ϕ\phiϕ的联合后验分布 p(θ,ϕ∣y)∝g(θ,ϕ)p(\theta, \phi | y) \propto g(\theta, \phi)p(θ,ϕ∣y)∝g(θ,ϕ)。
我们首先介绍Full Conditional Distribution：
在p(θ,ϕ∣y)=p(θ∣ϕ,y)p(ϕ∣y)p(\theta, \phi | y)=p(\theta|\phi,y)p(\phi|y)p(θ,ϕ∣y)=p(θ∣ϕ,y)p(ϕ∣y)中，分布p(θ∣ϕ,y)p(\theta|\phi, y)p(θ∣ϕ,y)称为θ\thetaθ的Full Conditional Distribution。在一些情况中，Full Conditional Distribution是一个我们知道如何抽样的标准分布，此时就无需抽取候选样本再计算接受率来决定究竟是接受它还是决绝它了。我们可以将Full Conditional Distribution当成是候选分布（即MH算法里的proposal candidate distribution），MH算法里的接受率α\alphaα为1。
由于p(θ∣ϕ,y)∝p(θ,ϕ∣y)∝g(θ,ϕ)p(\theta|\phi, y)\propto p(\theta, \phi | y)\propto g(\theta, \phi)p(θ∣ϕ,y)∝p(θ,ϕ∣y)∝g(θ,ϕ)，且p(ϕ∣θ,y)∝p(θ,ϕ∣y)∝g(θ,ϕ)p(\phi|\theta, y)\propto p(\theta, \phi | y)\propto g(\theta, \phi)p(ϕ∣θ,y)∝p(θ,ϕ∣y)∝g(θ,ϕ)，因此Full Conditional Distribution的好处在于它们都正比于Full Joint Posterior Distribution（即p(θ,ϕ∣y)p(\theta, \phi | y)p(θ,ϕ∣y)

Gibbs Sampling算法的过程如下：
假设现有参数θ\thetaθ和ϕ\phiϕ的联合后验分布 p(θ,ϕ∣y)p(\theta, \phi | y)p(θ,ϕ∣y),

首先，选取参数的初始值θ0,ϕ0\theta_0,\phi_0θ0,ϕ0
对i=1,...,mi=1,...,mi=1,...,m，重复如下操作：
（1）利用ϕi−1\phi_{i-1}ϕi−1，生成θi\theta_{i}θi~p(θ∣ϕi−1,y)p(\theta | \phi_{i-1},y)p(θ∣ϕi−1,y)
（2）利用θi\theta_iθi，生成ϕi\phi_{i}ϕi~p(ϕ∣θi,y)p(\phi | \theta_{i},y)p(ϕ∣θi,y)

Gibbs Sampling的思想是，更新时每次只采样一个参数，遍历完所有参数后再循环。在更新其中的一个参数时，我们利用的是其它参数当前状态的值。