对于一般的分布的采样,在很多的编程语言中都有实现,如最基本的满足均匀分布的随机数,但是对于复杂的分布,要想对其采样,却没有实现好的函数,在这里,可以使用马尔可夫链蒙特卡罗(Markov Chain Monte Carlo, MCMC)方法,其中Metropolis-Hastings采样和Gibbs采样是MCMC中使用较为广泛的两种形式。
MCMC的基础理论为马尔可夫过程,在MCMC算法中,为了在一个指定的分布上采样,根据马尔可夫过程,首先从任一状态出发,模拟马尔可夫过程,不断进行状态转移,最终收敛到平稳分布。
设XtX_t表示随机变量XX在离散时间tt时刻的取值。若该变量随时间变化的转移概率仅仅依赖于它的当前取值,即
P(Xt+1=sj∣X0=s0,X1=s1,⋯,Xt=si)=P(Xt+1=sj∣Xt=si)
P\left ( X_{t+1}=s_j\mid X_0=s_{0},X_1=s_{1}, \cdots ,X_t=s_i \right )=P\left ( X_{t+1}=s_j\mid X_t=s_i \right )
也就是说状态转移的概率只依赖于前一个状态。称这个变量为马尔可夫变量,其中,s0,s1,⋯,si,sj∈Ωs_0,s_1,\cdots ,s_i,s_j\in \Omega 为随机变量XX可能的状态。这个性质称为马尔可夫性质,具有马尔可夫性质的随机过程称为马尔可夫过程。
马尔可夫链指的是在一段时间内随机变量XX的取值序列(X0,X1,⋯,Xm)\left ( X_0,X_1,\cdots ,X_m \right ),它们满足如上的马尔可夫性质。
马尔可夫链是通过对应的转移概率定义的,转移概率指的是随机变量从一个时刻到下一个时刻,从状态sis_i转移到另一个状态sjs_j的概率,即:
P(i→j):=Pi,j=P(Xt+1=sj∣Xt=si)
P\left ( i\rightarrow j \right ):=P_{i,j}=P\left ( X_{t+1}=s_j\mid X_t=s_i \right )
记π(t)k\pi _k^{\left ( t \right )}表示随机变量XX在时刻tt的取值为sks_k的概率,则随机变量XX在时刻t+1t+1的取值为sis_i的概率为:
π(t+1)i=P(Xt+1=si)=∑kP(Xt+1=si∣Xt=sk)⋅P(Xt=sk)=∑kPk,i⋅π(t)k
\begin{align*} \pi _i^{\left ( t+1 \right )} &=P\left ( X_{t+1}=s_i \right ) \\ &= \sum_{k}P\left ( X_{t+1}=s_i\mid X_{t}=s_k \right )\cdot P\left ( X_{t}=s_k \right )\\ &= \sum_{k}P_{k,i}\cdot \pi _k^{\left ( t \right )} \end{align*}
假设状态的数目为nn,则有:
(π(t+1)1,⋯,π(t+1)n)=(π(t)1,⋯,π(t)n)⎡⎣⎢⎢⎢⎢⎢P1,1P2,1⋮Pn,1P1,2P2,2⋮Pn,2⋯⋯⋯P1,nP2,n⋮Pn,n⎤⎦⎥⎥⎥⎥⎥
\left ( \pi _1^{\left ( t+1 \right )},\cdots ,\pi _n^{\left ( t+1 \right )} \right )=\left ( \pi _1^{\left ( t \right )},\cdots ,\pi _n^{\left ( t \right )} \right )\begin{bmatrix} P_{1,1} & P_{1,2} & \cdots & P_{1,n}\\ P_{2,1} & P_{2,2} & \cdots & P_{2,n}\\ \vdots & \vdots & & \vdots \\ P_{n,1} & P_{n,2} & \cdots & P_{n,n} \end{bmatrix}
对于马尔可夫链,需要注意以下的两点:
如果一个马尔可夫过程既没有周期性,又不可约,则称为各态遍历的。
对于一个各态遍历的马尔可夫过程,无论初始值π(0)\pi ^{\left ( 0 \right )}取何值,随着转移次数的增多,随机变量的取值分布最终都会收敛到唯一的平稳分布π∗\pi ^{\ast },即:
limt→∞π(0)Pt=π∗
\underset{t\rightarrow \infty }{lim}\pi ^{\left ( 0 \right )}\mathbf{P}^t=\pi ^{\ast }
且这个平稳分布π∗\pi ^{\ast }满足:
π∗P=π∗
\pi ^{\ast }\mathbf{P}=\pi ^{\ast }
其中,P=(pi,j)n×n\mathbf{P}=\left ( p_{i,j} \right )_{n\times n}为转移概率矩阵。
对于一个给定的概率分布P(X)P\left (X \right ),若是要得到其样本,通过上述的马尔可夫链的概念,我们可以构造一个转移矩阵为P\mathbf{P}的马尔可夫链,使得该马尔可夫链的平稳分布为P(X)P\left (X \right ),这样,无论其初始状态为何值,假设记为x0x_0,那么随着马尔科夫过程的转移,得到了一系列的状态值,如:x0,x1,x2,⋯,xn,xn+1,⋯,x_0,x_1,x_2,\cdots ,x_n,x_{n+1},\cdots ,,如果这个马尔可夫过程在第nn步时已经收敛,那么分布P(X)P\left (X \right )的样本即为xn,xn+1,⋯x_n,x_{n+1},\cdots 。
对于一个各态遍历的马尔可夫过程,若其转移矩阵为P\mathbf{P},分布为π(x)\pi \left ( x \right ),若满足:
π(i)Pi,j=π(j)Pj,i
\pi \left ( i \right )P_{i,j}=\pi \left ( j \right )P_{j,i}
则π(x)\pi \left ( x \right )是马尔可夫链的平稳分布,上式称为细致平稳条件。
Metropolis采样算法是最基本的基于MCMC的采样算法。
假设需要从目标概率密度函数p(θ)p\left ( \theta \right )中进行采样,同时,θ\theta 满足−∞<θ<∞-\infty <\theta <\infty 。Metropolis采样算法根据马尔可夫链去生成一个序列:
θ(1)→θ(2)→⋯θ(t)→
\theta ^{\left ( 1 \right )}\rightarrow \theta ^{\left ( 2 \right )}\rightarrow \cdots \theta ^{\left (t \right )}\rightarrow
其中,θ(t) \theta ^{\left (t \right )}表示的是马尔可夫链在第tt代时的状态。
在Metropolis采样算法的过程中,首先初始化状态值θ(1)\theta ^{\left (1 \right )},然后利用一个已知的分布q(θ∣θ(t−1))q\left ( \theta \mid \theta ^{\left ( t-1 \right )} \right )生成一个新的候选状态θ(∗)\theta ^{\left (\ast \right )},随后根据一定的概率选择接受这个新值,或者拒绝这个新值,在Metropolis采样算法中,概率为:
α=min⎛⎝⎜1,p(θ(∗))p(θ(t−1))⎞⎠⎟
\alpha =min\: \left ( 1,\; \frac{p\left ( \theta ^{\left ( \ast \right )} \right )}{p\left ( \theta ^{\left ( t-1 \right )} \right )} \right )
这样的过程一直持续到采样过程的收敛,当收敛以后,样本θ(t)\theta ^{\left (t \right )}即为目标分布p(θ)p\left ( \theta \right )中的样本。
基于以上的分析,可以总结出如下的Metropolis采样算法的流程:
要证明Metropolis采样算法的正确性,最重要的是要证明构造的马尔可夫过程满足如上的细致平稳条件,即:
π(i)Pi,j=π(j)Pj,i
\pi \left ( i \right )P_{i,j}=\pi \left ( j \right )P_{j,i}
对于上面所述的过程,分布为p(θ)p\left ( \theta \right ),从状态ii转移到状态jj的转移概率为:
Pi,j=αi,j⋅Qi,j
P_{i,j} =\alpha _{i,j}\cdot Q_{i,j}
其中,Qi,jQ_{i,j}为上述已知的分布。
对于选择该已知的分布,在Metropolis采样算法中,要求该已知的分布必须是对称的,即Qi,j=Qj,iQ_{i,j}=Q_{j,i},即 q(θ=θ(t)∣θ(t−1))=q(θ=θ(t−1)∣θ(t)) q\left ( \theta =\theta ^{\left ( t \right )}\mid \theta ^{\left ( t-1 \right )} \right )=q\left ( \theta =\theta ^{\left ( t-1 \right )}\mid \theta ^{\left ( t \right )} \right ) 常用的符合对称的分布主要有:正态分布,柯西分布以及均匀分布等。
接下来,需要证明在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
\begin{align*} p\left ( \theta ^{\left ( i \right )} \right )P_{i,j} &=p\left ( \theta ^{\left ( i \right )} \right )\cdot \alpha _{i,j}\cdot Q_{i,j} \\ &= p\left ( \theta ^{\left ( i \right )} \right )\cdot min\; \left \{ 1,\frac{p\left ( \theta ^{\left ( j \right )} \right )}{p\left ( \theta ^{\left ( i \right )} \right )} \right \}\cdot Q_{i,j}\\ &=min\; \left \{ p\left ( \theta ^{\left ( i \right )} \right )Q_{i,j},p\left ( \theta ^{\left ( j \right )} \right )Q_{i,j} \right \}\\ &=p\left ( \theta ^{\left ( j \right )} \right )\cdot min\; \left \{ \frac{p\left ( \theta ^{\left ( i \right )} \right )}{p\left ( \theta ^{\left ( j \right )} \right )}, 1 \right \}\cdot Q_{j,i}\\ &=p\left ( \theta ^{\left ( j \right )} \right )\cdot \alpha _{j,i}\cdot Q_{j,i}\\ &=p\left ( \theta ^{\left ( j \right )} \right )P_{j,i} \end{align*}
因此,通过以上的方法构造出来的马尔可夫链是满足细致平稳条件的。
假设需要从柯西分布中采样数据,我们利用Metropolis采样算法来生成样本,其中,柯西分布的概率密度函数为:
f(θ)=1π(1+θ2)
f\left ( \theta \right )=\frac{1}{\pi \left ( 1+\theta ^2 \right )}
那么,根据上述的Metropolis采样算法的流程,接受概率α\alpha 的值为:
α=min⎛⎝⎜⎜1,1+[θ(t)]21+[θ(∗)]2⎞⎠⎟⎟
\alpha =min\; \left ( 1,\frac{1+\left [ \theta ^{\left ( t \right )} \right ]^2}{1+\left [ \theta ^{\left ( \ast \right )} \right ]^2} \right )
代码如下:
'''
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采样算法的推广形式。