一、Gibbs采样概述

前面介绍的Metropolis-Hastings采样为从指定分布中进行采样提供了一个统一的框架，但是采样的效率依赖于指定的分布的选择，若是选择的不好，会使得接受率比较低，大量的采样被拒绝，影响到整体的收敛速度。

Gibbs采样是Metropolis-Hastings采样算法的特殊形式，即找到一个已知的分布，使得接受率α=1<script type="math/tex" id="MathJax-Element-1">\alpha =1</script>。这样，每次的采样都会被接受，可以提高MCMC的收敛速度。

二、Gibbs采样算法的流程

在这部分，先直接给出Gibbs采样算法的流程，对于Gibbs采样算法的有效性将在第三部分给出论述，Gibbs采样算法的具体流程如下所述：

• 初始化时间t=1<script type="math/tex" id="MathJax-Element-197">t=1</script>
• 设置u=(u1,u2,⋯,uN)<script type="math/tex" id="MathJax-Element-198">\mathbf{u}=\left ( u_1,u_2,\cdots ,u_N \right )</script>的值，并初始化初始状态Θ(t)=u<script type="math/tex" id="MathJax-Element-199">\Theta ^{\left (t \right )}=\mathbf{u}</script>
• 重复以下的过程：
  
  • 令t=t+1<script type="math/tex" id="MathJax-Element-200">t=t+1</script>
  • 对每一维：i=1,2,⋯N<script type="math/tex" id="MathJax-Element-201">i=1,2,\cdots N</script>
    
    • θ(t)1∼p(θ1∣θ(t−1)2,⋯,θ(t−1)N)<script type="math/tex" id="MathJax-Element-202"> \theta _1^{\left ( t \right )}\sim p\left ( \theta _1\mid \theta _2^{\left ( t-1 \right )},\cdots ,\theta _N^{\left ( t-1 \right )} \right )</script>
    • θ(t)2∼p(θ2∣θ(t)1,⋯,θ(t−1)N)<script type="math/tex" id="MathJax-Element-203">\theta _2^{\left ( t \right )}\sim p\left ( \theta _2\mid \theta _1^{\left ( t \right )},\cdots ,\theta _{N}^{\left ( t-1 \right )} \right )</script>
    • ⋯<script type="math/tex" id="MathJax-Element-204">\cdots </script>
    • θ(t)N−1∼p(θN−1∣θ(t)1,⋯,θ(t−1)N)<script type="math/tex" id="MathJax-Element-205">\theta _{N-1}^{\left ( t \right )}\sim p\left ( \theta _{N-1}\mid \theta _1^{\left ( t \right )},\cdots ,\theta _{N}^{\left ( t-1 \right )} \right )</script>
    • θ(t)N∼p(θN∣θ(t)1,⋯,θ(t)N−1)<script type="math/tex" id="MathJax-Element-206">\theta _N^{\left ( t \right )}\sim p\left ( \theta _N\mid \theta _1^{\left ( t \right )},\cdots ,\theta _{N-1}^{\left ( t \right )} \right )</script>
• 直到t=T<script type="math/tex" id="MathJax-Element-207">t=T</script>

Gibbs采样有一个缺陷，必须知道条件分布。

三、上述过程满足细致平稳条件

为简单起见，我们假设所需采样的分布为一个二元分布f(x,y)<script type="math/tex" id="MathJax-Element-13">f\left ( x,y \right )</script>，假设两个状态为(x1,y1)<script type="math/tex" id="MathJax-Element-14">\left ( x_1,y_1 \right )</script>和(x1,y2)<script type="math/tex" id="MathJax-Element-15">\left ( x_1,y_2 \right )</script>。已知：

<script type="math/tex; mode=display" id="MathJax-Element-16">p\left ( x_1,y_1 \right )\cdot p\left ( y_2\mid x_1 \right )=p\left ( x_1 \right )\cdot p\left ( y_1\mid x_1 \right )\cdot p\left ( y_2\mid x_1 \right )</script>

<script type="math/tex; mode=display" id="MathJax-Element-17">p\left ( x_1,y_2 \right )\cdot p\left ( y_1\mid x_1 \right )=p\left ( x_1 \right )\cdot p\left ( y_2\mid x_1 \right )\cdot p\left ( y_1\mid x_1 \right )</script>

所以有：

<script type="math/tex; mode=display" id="MathJax-Element-18">p\left ( x_1,y_1 \right )\cdot p\left ( y_2\mid x_1 \right )=p\left ( x_1,y_2 \right )\cdot p\left ( y_1\mid x_1 \right )</script>

由此可见，Gibbs采样的过程是满足细致平稳条件的。这里直接取p(y2∣x1)<script type="math/tex" id="MathJax-Element-19">p\left ( y_2\mid x_1 \right )</script>为转移概率，则α=1<script type="math/tex" id="MathJax-Element-20">\alpha =1</script>，可见Gibbs采样算法是Metropolis-Hastings采样的特殊形式。

四、实验

4.1、前提

假设从二项正态分布中进行采样，假设Θ=(θ1,θ2)<script type="math/tex" id="MathJax-Element-274">\Theta =\left ( \theta _1,\theta _2 \right )</script>，且：

<script type="math/tex; mode=display" id="MathJax-Element-241">\Theta \sim Norm\left (\mu ,\Sigma \right )</script>

其中

<script type="math/tex; mode=display" id="MathJax-Element-246">\mu =\left ( \mu _1,\mu _2 \right )</script>

<script type="math/tex; mode=display" id="MathJax-Element-254">\Sigma =\begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}</script>

已知：

<script type="math/tex; mode=display" id="MathJax-Element-259">\theta _1\sim Norm\left ( \mu _1+\rho \left ( \theta _2-\mu _2 \right ), \sqrt{1-\rho ^2} \right )</script>

<script type="math/tex; mode=display" id="MathJax-Element-276">\theta _2\sim Norm\left ( \mu _2+\rho \left ( \theta _1-\mu _1 \right ), \sqrt{1-\rho ^2} \right )</script>

4.2、流程

• 初始化时间t=1<script type="math/tex" id="MathJax-Element-358">t=1</script>
• 设置u=(u1,u2)<script type="math/tex" id="MathJax-Element-359">\mathbf{u}=\left ( u_1,u_2 \right )</script>的值，并初始化初始状态Θ(t)=u<script type="math/tex" id="MathJax-Element-360">\Theta ^{\left (t \right )}=\mathbf{u}</script>
• 重复以下的过程：
  
  • 令t=t+1<script type="math/tex" id="MathJax-Element-361">t=t+1</script>
  • 对每一维：i=1,2<script type="math/tex" id="MathJax-Element-362">i=1,2</script>
    
    • θ(t)1∼Norm(μ1+ρ(θ2−μ2),1−ρ2−−−−−√)<script type="math/tex" id="MathJax-Element-363"> \theta _1^{\left ( t \right )}\sim Norm\left ( \mu _1+\rho \left ( \theta _2-\mu _2 \right ), \sqrt{1-\rho ^2} \right )</script>
    • θ(t)2∼Norm(μ2+ρ(θ1−μ1),1−ρ2−−−−−√)<script type="math/tex" id="MathJax-Element-364">\theta _2^{\left ( t \right )}\sim Norm\left ( \mu _2+\rho \left ( \theta _1-\mu _1 \right ), \sqrt{1-\rho ^2} \right )</script>
• 直到t=T<script type="math/tex" id="MathJax-Element-365">t=T</script>

4.3、实验代码

'''
Date:20160704
@author: zhaozhiyong
'''
import random
import math
import matplotlib.pyplot as plt

def p_ygivenx(x, m1, m2, s1, s2):
    return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))

def p_xgiveny(y, m1, m2, s1, s2):
    return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))

N = 5000
K = 20
x_res = []
y_res = []
m1 = 10
m2 = -5
s1 = 5
s2 = 2

rho = 0.5
y = m2

for i in xrange(N):
    for j in xrange(K):
        x = p_xgiveny(y, m1, m2, s1, s2)
        y = p_ygivenx(x, m1, m2, s1, s2)
        x_res.append(x)
        y_res.append(y)

num_bins = 50
plt.hist(x_res, num_bins, normed=1, facecolor='green', alpha=0.5)
plt.hist(y_res, num_bins, normed=1, facecolor='red', alpha=0.5)
plt.title('Histogram')
plt.show()

4.4、实验结果