简单易学的机器学习算法—Gibbs采样

一、Gibbs采样概述
前面介绍的Metropolis-Hastings采样为从指定分布中进行采样提供了一个统一的框架，但是采样的效率依赖于指定的分布的选择，若是选择的不好，会使得接受率比较低，大量的采样被拒绝，影响到整体的收敛速度。
Gibbs采样是Metropolis-Hastings采样算法的特殊形式，即找到一个已知的分布，使得接受率α=1。这样，每次的采样都会被接受，可以提高MCMC的收敛速度。
二、Gibbs采样算法的流程
在这部分，先直接给出Gibbs采样算法的流程，对于Gibbs采样算法的有效性将在第三部分给出论述，Gibbs采样算法的具体流程如下所述：
初始化时间t=1
设置u=(u1,u2,⋯,uN)的值，并初始化初始状态Θ(t)=u
重复以下的过程：
令t=t+1
对每一维：i=1,2,⋯N

直到t=T
Gibbs采样有一个缺陷，必须知道条件分布。
三、上述过程满足细致平稳条件
为简单起见，我们假设所需采样的分布为一个二元分布f(x,y)，假设两个状态为(x1,y1)和(x1,y2)。已知：

所以有：

由此可见，Gibbs采样的过程是满足细致平稳条件的。这里直接取p(y2∣x1)为转移概率，则α=1，可见Gibbs采样算法是Metropolis-Hastings采样的特殊形式。
四、实验
4.1、前提
假设从二项正态分布中进行采样，假设Θ=(θ1,θ2)，且：

其中

已知：

4.2、流程

初始化时间t=1
设置u=(u1,u2)的值，并初始化初始状态Θ(t)=u
重复以下的过程：
令t=t+1
对每一维：i=1,2

直到t=T
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、实验结果