对于一般的分布的采样，在很多的编程语言中都有实现，如最基本的满足均匀分布的随机数，但是对于复杂的分布，要想对其采样，却没有实现好的函数，在这里，可以使用马尔可夫链蒙特卡罗(Markov Chain Monte Carlo, MCMC)方法，其中Metropolis-Hastings采样和Gibbs采样是MCMC中使用较为广泛的两种形式。
MCMC的基础理论为马尔可夫过程，在MCMC算法中，为了在一个指定的分布上采样，根据马尔可夫过程，首先从任一状态出发，模拟马尔可夫过程，不断进行状态转移，最终收敛到平稳分布。
一、马尔可夫链
1、马尔可夫链
设 Xt 表示随机变量 X 在离散时间 t 时刻的取值。若该变量随时间变化的转移概率仅仅依赖于它的当前取值，即
P(Xt+1=sj∣X0=s0,X1=s1,⋯,Xt=si)=P(Xt+1=sj∣Xt=si)
也就是说状态转移的概率只依赖于前一个状态。称这个变量为马尔可夫变量，其中， s0,s1,⋯,si,sj∈Ω 为随机变量 X 可能的状态。这个性质称为马尔可夫性质，具有马尔可夫性质的随机过程称为马尔可夫过程。
马尔可夫链指的是在一段时间内随机变量 X 的取值序列 (X0,X1,⋯,Xm) ，它们满足如上的马尔可夫性质。
2、转移概率
马尔可夫链是通过对应的转移概率定义的，转移概率指的是随机变量从一个时刻到下一个时刻，从状态 si 转移到另一个状态 sj 的概率，即：
P(i→j):=Pi,j=P(Xt+1=sj∣Xt=si)
记 π(t)k 表示随机变量 X 在时刻 t 的取值为 sk 的概率，则随机变量 X 在时刻 t+1 的取值为 si 的概率为：
π(t+1)i=P(Xt+1=si)=∑kP(Xt+1=si∣Xt=sk)⋅P(Xt=sk)=∑kPk,i⋅π(t)k
假设状态的数目为 n ，则有：
(π(t+1)1,⋯,π(t+1)n)=(π(t)1,⋯,π(t)n)⎡⎣⎢⎢⎢⎢⎢P1,1P2,1⋮Pn,1P1,2P2,2⋮Pn,2⋯⋯⋯P1,nP2,n⋮Pn,n⎤⎦⎥⎥⎥⎥⎥
3、马尔可夫链的平稳分布
对于马尔可夫链，需要注意以下的两点：
- 1、周期性：即经过有限次的状态转移，又回到了自身；
- 2、不可约：即两个状态之间相互转移；
如果一个马尔可夫过程既没有周期性，又不可约，则称为各态遍历的。
对于一个各态遍历的马尔可夫过程，无论初始值 π(0) 取何值，随着转移次数的增多，随机变量的取值分布最终都会收敛到唯一的平稳分布 π∗ ，即：
limt→∞π(0)Pt=π∗
且这个平稳分布 π∗ 满足：
π∗P=π∗
其中， P=(pi,j)n×n 为转移概率矩阵。
二、马尔可夫链蒙特卡罗方法
1、基本思想
对于一个给定的概率分布 P(X) ，若是要得到其样本，通过上述的马尔可夫链的概念，我们可以构造一个转移矩阵为 P 的马尔可夫链，使得该马尔可夫链的平稳分布为 P(X) ，这样，无论其初始状态为何值，假设记为 x0 ，那么随着马尔科夫过程的转移，得到了一系列的状态值，如： x0,x1,x2,⋯,xn,xn+1,⋯, ，如果这个马尔可夫过程在第 n 步时已经收敛，那么分布 P(X) 的样本即为 xn,xn+1,⋯ 。
2、细致平稳条件
对于一个各态遍历的马尔可夫过程，若其转移矩阵为 P ，分布为 π(x) ，若满足：
π(i)Pi,j=π(j)Pj,i
则 π(x) 是马尔可夫链的平稳分布，上式称为细致平稳条件。
3、Metropolis采样算法
Metropolis采样算法是最基本的基于MCMC的采样算法。
3.1、Metropolis采样算法的基本原理
假设需要从目标概率密度函数 p(θ) 中进行采样，同时， θ 满足 −∞<θ<∞ 。Metropolis采样算法根据马尔可夫链去生成一个序列：
θ(1)→θ(2)→⋯θ(t)→
其中， θ(t) 表示的是马尔可夫链在第 t 代时的状态。
在Metropolis采样算法的过程中，首先初始化状态值 θ(1) ，然后利用一个已知的分布 q(θ∣θ(t−1)) 生成一个新的候选状态 θ(∗) ，随后根据一定的概率选择接受这个新值，或者拒绝这个新值，在Metropolis采样算法中，概率为：
α=min⎛⎝⎜1,p(θ(∗))p(θ(t−1))⎞⎠⎟
这样的过程一直持续到采样过程的收敛，当收敛以后，样本 θ(t) 即为目标分布 p(θ) 中的样本。
3.2、Metropolis采样算法的流程
基于以上的分析，可以总结出如下的Metropolis采样算法的流程：
- 初始化时间 t=1
- 设置 u 的值，并初始化初始状态 θ(t)=u
- 重复一下的过程：
- 令 t=t+1
- 从已知分布 q(θ∣θ(t−1)) 中生成一个候选状态 θ(∗)
- 计算接受的概率： α=min(1,p(θ(∗))p(θ(t−1)))
- 从均匀分布 Uniform(0,1) 生成一个随机值 a
- 如果 a⩽α ，接受新生成的值： θ(t)=θ(∗) ；否则： θ(t)=θ(t−1)
- 直到 t=T
3.3、Metropolis算法的解释
要证明Metropolis采样算法的正确性，最重要的是要证明构造的马尔可夫过程满足如上的细致平稳条件，即：
π(i)Pi,j=π(j)Pj,i
对于上面所述的过程，分布为 p(θ) ，从状态 i 转移到状态 j 的转移概率为：
Pi,j=αi,j⋅Qi,j
其中， Qi,j 为上述已知的分布。
对于选择该已知的分布，在Metropolis采样算法中，要求该已知的分布必须是对称的，即 Qi,j=Qj,i ，即
q(θ=θ(t)∣θ(t−1))=q(θ=θ(t−1)∣θ(t))
常用的符合对称的分布主要有：正态分布，柯西分布以及均匀分布等。
接下来，需要证明在Metropolis采样算法中构造的马尔可夫链满足细致平稳条件。
p(θ(i))Pi,j=p(θ(i))⋅αi,j⋅Qi,j=p(θ(i))⋅min⎧⎩⎨⎪⎪1,p(θ(j))p(θ(i))⎫⎭⎬⎪⎪⋅Qi,j=min{p(θ(i))Qi,j,p(θ(j))Qi,j}=p(θ(j))⋅min⎧⎩⎨⎪⎪p(θ(i))p(θ(j)),1⎫⎭⎬⎪⎪⋅Qj,i=p(θ(j))⋅αj,i⋅Qj,i=p(θ(j))Pj,i
因此，通过以上的方法构造出来的马尔可夫链是满足细致平稳条件的。
3.4、实验
假设需要从柯西分布中采样数据，我们利用Metropolis采样算法来生成样本，其中，柯西分布的概率密度函数为：
f(θ)=1π(1+θ2)
那么，根据上述的Metropolis采样算法的流程，接受概率 α 的值为：
α=min⎛⎝⎜⎜1,1+[θ(t)]21+[θ(∗)]2⎞⎠⎟⎟
代码如下：
''' Date:20160629 @author: zhaozhiyong '''
import random
from scipy.stats import norm
import matplotlib.pyplot as plt
def cauchy(theta):
y = 1.0 / (1.0 + theta ** 2)
return y
T = 5000
sigma = 1
thetamin = -30
thetamax = 30
theta = [0.0] * (T+1)
theta[0] = random.uniform(thetamin, thetamax)
t = 0
while t < T:
t = t + 1
theta_star = norm.rvs(loc=theta[t - 1], scale=sigma, size=1, random_state=None)
#print theta_star
alpha = min(1, (cauchy(theta_star[0]) / cauchy(theta[t - 1])))
u = random.uniform(0, 1)
if u <= alpha:
theta[t] = theta_star[0]
else:
theta[t] = theta[t - 1]
ax1 = plt.subplot(211)
ax2 = plt.subplot(212)
plt.sca(ax1)
plt.ylim(thetamin, thetamax)
plt.plot(range(T+1), theta, 'g-')
plt.sca(ax2)
num_bins = 50
plt.hist(theta, num_bins, normed=1, facecolor='red', alpha=0.5)
plt.show()实验的结果：
对于Metropolis采样算法，其要求选定的分布必须是对称的，为了弥补这样的一个缺陷，在下一篇中，介绍一下Metropolis-Hastings采样算法，其是Metropolis采样算法的推广形式。
参考文献
- 1、马尔可夫链蒙特卡罗算法
- 2、受限玻尔兹曼机（RBM）学习笔记（一）预备知识
- 3、LDA数学八卦
