Skip to main content

4.4 马尔可夫链蒙特卡洛方法

到目前为止,前面讲述的对分布采样的方法称为直接采样法,在这种情况下概率分布通常是已知的,有些时候也能很容易地计算出归一化常数。然而直接采样方法在某些情况下会遇到困难甚至不适用,因此本章介绍另一种(也是非常重要的一种)基于马尔可夫链的梅特罗波利斯算法。

直接采样算法的遇到的第一个困难是,当随机变量的维度变得越来越高时,其采样效率变得越来越低,因此不适合于对高维分布进行采样。考虑随机向量: Xf(x1,x2,,xs)X\sim f(x_1,x_2,\cdots,x_s ),其中ss为随机向量XX的维度,每个分量的定义域为:biaib_i-a_i,则直接采样方法中取舍采样的效率为:

η=1Li=1s(biai) \eta= \cfrac{1}{L}\prod^{s}_{i=1}(b_i-a_i)

(式1

其中,L=maxf(x1,x2,,xs)L=\max{f(x_1,x_2,\cdots,x_s)},如果LL值很大,或者维数ss很高,则分母值很大,取舍算法的接受概率很低,拒绝率很高,因此采样效率非常低。直观上看,这是由于维数越高,用于包围目标概率分布的建议概率分布与目标概率分布之间存在更多的间隙空间,这些间隙空间即是被拒绝的采样点。

直接采样方法的另一个困难来自于概率分布为不完全已知概率分布,此时直接采样方法完全不适应。不完全已知概率分布是指,虽然概率分布的形式已知,但是其中包含部分项或本身就是很难进行数值计算的,例如为很难进行积分计算的高维积分,因此相当于是未知的。

一个直接的例子是渲染中每一帧图像的分布函数,我们可以将一张图像看成是某个值的概率分布,由于图像的每个像素的值是由RGB三个颜色值组成,因此这个值可以是这三个分量某种形式的函数,或者说我们可以以一个图像的亮度作为这个分布函数的值。在这个分布中,每一个像素位置的值是光线通过无数表面反射/折射传播的结果,因此每一个值本身就是一个高维积分计算,我们将在下一章看到渲染方程的路径形式,每一条路径都是一个高维积分。

本节要讨论的梅特罗波利斯算法是一种能够从任意概率分布p(x)p(x)中进行采样的方法,它提供一种机制使我们可以通过计算一个正比于目标分布函数p(x)p(x)的分布函数f(x)f(x)的值来获取采样值。由于f(x)f(x)只需要正比于目标分布,而不需要精确计算目标分布的归一化常数,因此梅特罗波利斯算法可以用于从任意分布采样。

马尔可夫链

梅特罗波利斯算法可以看做是取舍算法的一种推广,取舍算法每一步从建议分布产生一个全局的随机数,并与目标概率分布进行比较,以此来决定对新随机数的接受与拒绝,然而取舍算法的建议概率分布和目标概率分布之间的间隙随着随机变量维数的增加而增加,因此其效率极低。

而梅特罗波利斯算法以马尔可夫链为基础,整个取样的过程产生一条马尔可夫链,它将目标概率分布当做一个稳态的系统,每个新产生的随机数根据该系统的动态平衡关系进行取舍,所以最终整条马尔可夫链上的随机数序列就是服从目标概率分布的,并且梅特罗波利斯算法只需要能够计算目标概率分布每个定义域处的值,而不需要计算其归一化常数(注意:计算p(x)p(x)的每个值和计算p(x)p(x)的归一化常数是不一样的,后者实际上是积分p(x)dx{\rm \int} p(x){\rm d}x的计算,这个积分计算可能非常困难。)。

所以为了理解梅特罗波利斯算法,我们需要首先了解马尔可夫链(Markov chain)的概念及其属性,这是推导和理解梅特罗波利斯算法的重要基础。基于马尔可夫链的蒙特卡洛方法称为马尔可夫蒙特卡洛方法(Markov chain Monte Carlo, MCMC),而梅特罗波利斯算法是最重要的一种马尔可夫蒙特卡洛方法。

设一个系统状态序列为x0,x1,x2,x_0,x_1,x_2,\cdots,每个状态可以以一定的概率转移到另一个状态,假设我们从其中任意一个起始状态开始行走,每个时刻由一个状态转移到另一个状态(也可以转移到自身),如果对任何时刻nn,系统状态的转移条件概率为:

P(Xn+1=xX0=x0,X1=x1,,Xn=xn)=P(Xn+1=xXn=xn) P(X_{n+1}=x|X_0=x_0,X_1=x_1,\cdots,X_n=x_n)=P(X_{n+1}=x|X_n=x_n)

(式2

则此转移过程形成的状态序列X0,X1,,Xn+1X_0,X_1,\cdots,X_{n+1}称为马尔可夫链(注意这里XnX_n表示nn时刻产生的一个随机数,而xnx_n表示nn时刻的一个系统状态值,其中两个相邻时刻的系统状态值可能是相同的。)。马尔可夫链的定义说明,其将来时刻n+1n+1的状态,只依赖于当前时刻nn的状态,与以前任何时刻的状态都无关。据此,系统状态序列x0,x1,,xn,xn+1x_0,x_1,\cdots,x_n,x_{n+1}发生的概率可分解为:

P(x0,x1,,xn,xn+1)=P(X0=x0)P(X1=x1X0=x0)P(Xn+1=xn+1Xn=xn) P(x_0,x_1,\cdots,x_n,x_{n+1})=P(X_0=x_0)P(X_1=x_1|X_0=x_0)\cdots P(X_{n+1}=x_{n+1}|X_n=x_n)

(式3

式中,P(X0=x0)P(X_0=x_0)为初始状态的概率,通常是从初始状态x0x_0出发,所以P(X0=x0)=1P(X_0=x_0)=1;式中条件概率P(Xn+1=xXn=y)P(X_{n+1}=x|X_n=y)称为一步转移概率,简称为转移概率,表示在n+1n+1时刻由yy状态转移到xx状态,转移概率通常与时刻nn有关。若转移概率P(Xn+1=xXn=y)P(X_{n+1}=x|X_n=y)与时刻nn无关,即P(Xn+1=xn+1Xn=xn)=P(Xn=xXn1=y)P(X_{n+1}=x_{n+1}|X_n=x_n)=P(X_n=x|X_{n-1}=y),则称马尔可夫链是均匀的,或者时齐的(time-homogeneous),通常我们研究的都是时齐的马尔可夫链,因此后面的表述我们将去掉时间nn相关信息,yyxx的状态条件转移概率直接表述为: P(xy)P(x|y)

马尔可夫链将一个系统的定义域当做一个状态空间,因此马尔可夫链又分为离散时间马尔可夫链(discrete-time Markov chain, DTMC)和连续时间马尔可夫链(continuous-time Markov chain, CTMC)。

马尔可夫链是在一个状态空间中从一个状态到另一个状态的不断转换的随机过程,这个过程中经过的每个状态值就是我们获得的服从这个状态空间概率分布的采样值。在一个平衡的系统状态中,每个状态到另一个状态的转移概率是固定的(我们可以将系统状态表述为一个状态图,如图(1)是一个具有两个状态的状态图),这些转移概率的结果就是使得状态空间的每个状态服从系统的分布,因此马尔可夫链是服从系统分布的一个随机采样序列。

初学者在理解马尔可夫链的时候往往比较吃力,这里尝试做一些解释,请读者细细体会。统计方法是关于用概率来解决确定性问题的方法,它的实际操作方式就是用大量的重复试验以获取足够的样本,而这些样本是满足目标概率分布的,所以能够用来解决确定性问题。具体地,以前面讲述的蒙特卡洛方法为例,当我们从一个概率分布p(x)p(x)采样产生一个随机数序列时,这个随机数序列里的值并不是唯一的,而是包含大量重复的值,每个随机数对应的概率越高,则其重现的次数越多,最终每个唯一随机数xix_i在这个随机序列中出现的次数正比于它的概率pip_i

所以在马尔可夫链中也是一样,给定一个初始状态,马尔可夫链按照各个状态的转移概率在系统状态中随机行走,由于这些转移概率是按照系统状态的概率决定的,因此马尔可夫链就是一条状态链,这里马尔可夫链里出现的每个状态就是一个随机值,那些具有较高概率的系统状态,则其在马尔可夫链中出现的次数越多。因此,马尔可夫链本质上和前面讲述的其他采样方法是一样的,其目的都是产生服从目标概率分布的采样,但是操作方法不一样。

图(1):一个具有两个转换状态的马尔可夫链(图片来自Wikipedia)

在一个给定平衡状态概率分布的系统中,如图(1)就是一个包含两个状态的系统,每个状态到其他状态的转移概率是固定的(实际上对于一个有限状态空间,整个系统状态的转移可以表述为一个转移矩阵,每个状态的矢量只需要乘以这个矩阵就能按照相关的概率进行转移,这里的例子参见Examples_of_Markov_chains;在本书后面的辐射度(radiosity)方法里,其思路就是将整个场景分成多个有限的面积块,然后由几何关系计算出每个面积块到其他面积块的转移比例,即是概率,然后给定任意一个光源,经过一定时间的反射,就能得到一个稳定的光照分布。),因此给定任意一个初始状态,例如晴天为100%100\%,雨天0%0\%,则等到一定时间tt后系统处于平衡状态时,期间产生的马尔可夫链中晴天和雨天的分布一定服从平衡状态下的概率分布,即:约83.3%83.3\%的时间为晴天,剩下为雨天。

所以,利用马尔可夫链对目标概率分布进行采样的过程,实际上就是设一个状态系统服从目标概率分布,然后给定一个任意初始值,经历一定时间等系统趋于平衡之后,产生的马尔可夫链就是服从目标概率分布的,因此可作为服从该目标概率分布的采样值。由此也可以看出马尔可夫蒙特卡洛方法的一个问题,那就是开始一段趋近稳定状态的过程存在很大的方差,实践中往往需要通过某种方法省略掉这部分采样值。

而在马尔可夫蒙特卡洛方法中,梅特罗波利斯算法主要就是以一种简单的方式从已知概率分布中产生每个状态下正确的转移概率,要理解梅特罗波利斯算法的推导过程,我们还需要了解马尔可夫链的一些属性。

马尔可夫链属性

马尔可夫链有许多属性,为了简化概念,这里只讨论跟后面推导梅特罗波利斯算法相关的属性。

不可约性

如果一个系统从状态ii开始,存在非零的概率经过n1n\geq 1步之后可以达到状态jj,则称ii可达到(accessible)jj,记为iji\rightarrow j,表述为:

P(Xnij=jX0=i)=pij(nij)>0 P(X_{n_{ij}}=j|X_0=i)=p^{(n_{ij})}_{ij}>0

(式4

如果iji\rightarrow jjij\rightarrow i同时存在(nijn_{ij}njin_{ji}可能不相等),则称iijj是联通的(communicating)。一个联通集合(communicating class)是指存在一个最大状态集合CC,使得其中的每对状态相互都是联通的。

如果一个状态空间就是一个单一的联通集合,则称该状态空间内的马尔可夫链是不可约的(irreducible),也就是说,在该状态空间中,可以由任何一个状态ii出发,经过nijn_{ij}步之后,达到另外一个任意状态jj

非周期性

假设从一个状态ii出发,经过一定的步数niin_{ii}后回到状态ii自身,如果niin_{ii}kk的倍数,则称状态ii具有周期(period)kk,称为周期态,一个状态的周期表述为:

k=gcd{n>0:P(Xn=i)X0=i)>0} k=gcd\{n>0:P(X_n=i)|X_0=i)>0\}

(式5

这里gcdgcd表示取最大公约数(注意:即使状态ii拥有最大公约数kk,但并不意味着从ii出发经过kk步之后一定可以回到状态ii,例如可能回到原状态的步数为{6,8,10,12,}\{6,8,10,12,\cdots\},这里kk为2,但是2并没有出现在返回步数列表中。)(greatest common divisor),如果k=1k=1,则称该状态是非周期态(aperiodic state),它意味着回到状态自身的步数可以是任意的;如果至少有一个状态是非周期态,则称为非周期的马尔可夫链;若状态ii的周期数k>1k>1,则称其为周期态。

回返性

假设从一个状态ii出发,经过一定步数nn之后回到状态ii的概率记为:

pii(n)=P(Ti=n) p^{(n)}_{ii}=P(T_i=n)

(式6

概率pii(n)p^{(n)}_{ii}又称为首达概率,如果平均返回时间Ti<T_i<\infty,则称状态ii为正回返(positive recurrent)状态,如果Ti=T_i=\infty,则称其为零回返状态,或者暂态(transient)。有限马尔可夫链的回返状态必为正回返状态。

遍历性

如果状态ii是正回返,不可约和非周期状态,则称状态ii是各态遍历(ergodic)状态,或遍历态。换句话说,如果一个状态具有周期为1,并且能够在有限的平均步数内回返,则为遍历态。如果一个不可约的马尔可夫链中所有的状态都是遍历的,则该马尔可夫链是遍历的。

更一般地,如果一个马尔可夫链是遍历的,则始终存在一个步数nn,使得从任意其他非ii状态经过nn步之后可以达到ii状态。也就是说,马尔可夫链达到ii状态的概率与初始状态无关,其实遍历性就是初始状态独立性。遍历性是马尔可夫链收敛到平稳分布的定量条件。

平稳分布

设马尔可夫链在一个状态空间中从某个初始状态出发经历一定的时间后,各个状态jj的离散分布概率为πj\pi_{j},如果对任意状态jj存在概率密度函数πj>0\pi_j>0,并使得:

i=1nπipij=πj \sum^{n}_{i=1}\pi_i p_{ij}=\pi_j

(式7

即从所有其他状态转移到jj状态的转移概率之和为状态jj的概率分布,记为:

πP=π \mathbf{\pi} \mathbf{P}=\mathbf{\pi}

(式8

则称马尔可夫链达到总体平衡,这里P\mathbf{P}是一个转移矩阵,π\mathbf{\pi}为所有状态的概率构成的矢量。一个正回返和非周期的马尔可夫链,当且仅当它是一个遍历链时,经历足够多的状态转移步数,这个马尔可夫链具有一个不变的概率分布πj>0\pi_j>0,这时无论初始分布如何,最终都会趋向于平稳概率分布πj\pi_j,于是有:

limmpijn=πj \lim_{m\rightarrow \infty}p^{n}_{ij}=\pi_j

(式9

这个唯一的近似不变的概率分布就是平稳分布(stationary distribution),它是一个极限分布。这就表明,最终状态不再发生变化,系统达到平稳状态。马尔可夫链具有平稳分布,是马尔可夫链的一个重要特征,在马尔可夫链蒙特卡洛方法中具有重要意义。

我们可以以前面的天气预报的状态分布为例来分析马尔可夫链从任意一个状态出发,经历一定时间后逐渐趋近于平稳分布的过程。如图(1)所示,这里状态空间的转移矩阵为:

P=[0.90.10.50.5] \mathbf{P}=\begin{bmatrix} 0.9 & 0.1\\0.5 & 0.5 \end{bmatrix}

(式10

设初始状态为晴天概率为100%100\%,雨天为0%0\%,用一个矢量表述为:

x(0)=[1.00.0] \mathbf{x}^{(0)}=\begin{bmatrix} 1.0&0.0 \end{bmatrix}

(式11

则经过一次随机行走之后状态概率分布为:

x(1)=x(0)P=[1.00.0][0.90.10.50.5]=[0.90.1] \mathbf{x}^{(1)}=\mathbf{x}^{(0)}\mathbf{P}=\begin{bmatrix} 1.0&0.0 \end{bmatrix}\begin{bmatrix} 0.9&0.1\\0.5&0.5 \end{bmatrix}=\begin{bmatrix} 0.9&0.1 \end{bmatrix}

(式12

即依照这样的概率分布,第二天为晴天的概率为90%90\%,为雨天的概率相应为10%10\%,以此类推我们可以得到第nn天的状态概率分布,如果nn足够大,使系统趋近于稳定分布,则最终的天气预报的概率为:

[q1q2]=[0.8330.167] \begin{bmatrix} q_1&q_2 \end{bmatrix}=\begin{bmatrix} 0.833&0.167 \end{bmatrix}

(式13

细致平衡

在一个状态空间中,如果存在概率分布π\pi,使得:

πipij=πjpji \pi_i p_{ij}=\pi_j p_{ji}

(式14

该条件称为细致平衡(detailed balance),又称为局部平衡(local balance)。细致平衡是微观平衡,比上述的平稳分布(宏观平衡)具有更多的约束条件。

由细致平衡可以得到平稳分布,反之则不成立,平稳分布并不能保证细致平衡。满足细致平衡条件的马尔可夫链是可逆的(细致平衡本就是源于在一个平衡的动力学系统中,每个基础过程与其逆过程保持平衡。),称为可逆链,因此细致平衡与可逆性是等价的,可逆性保证了平稳分布。

梅特罗波利斯算法

有了平稳分布的条件,对于使用马尔可夫链的蒙特卡洛方法的思路则比较容易理解:它首先把目标概率分布的作用域看做一个状态空间,然后从任意一个状态出发随机行走以建立马尔可夫链,这个随机行走过程中的每一步遵循状态空间的转移概率分布,则经过一定时间之后,马尔可夫链中的状态概率分布会趋向于目标概率分布。在这个过程中,最重要的是怎样根据已知的目标概率求出转移概率分布,梅特罗波利斯算法(Metropolis algorithm)的核心就是使用了一种简单求解转移概率的方法。

为了更好地理解梅特罗波利斯算法及其最关键的概念,我们有必要区分一下统计学中利用马尔可夫链解决积分问题的两种不同的思路:第一种是已知目标概率分布pip_i,使用蒙特卡洛方法对这个目标概率分布进行采样,然后计算积分,此种方法的关键是要根据已知目标概率求出状态转移概率pijp_{ij},这样使得马尔可夫链能够逼近目标分布,这即是本节讨论的梅特罗波利斯算法;第二种是已知转移概率分布,然后给定任意一个初始概率分布,则马尔可夫链会根据转移概率逐渐趋向于稳定分布,这种思路在计算机图形学中被用于本书后面第(8)章会讨论的辐射度方法(radiosity)中;当然后者是一种迭代法(而不是蒙特卡洛方法),但是它与马尔可夫链性质的运用其实有很深的联系,在学习辐射度理论的时候可以细细体会这种联系。

梅特罗波利斯算法通过逐渐产生一序列的随机数的方式工作,这个随机数序列是以迭代的方式产生的,在这个迭代的过程中,下一个随机数的产生仅依赖于上一个随机数(因此整个随机数序列是一个马尔可夫链),随着产生越来越多的随机数,这个随机数序列会逐渐逼近于目标概率分布。

特别地,在每一次迭代中,该算法基于当前随机数产生一个新的候选随机数,然后这个候选随机数根据一定的概率被接受(这种情况下新的候选随机数会被用于下一个迭代中)或拒绝(这种情况下新的候选随机数被丢弃,当前随机数继续被用在下一个迭代中),这个接受概率被适当选取以使其能够近似目标概率的分布。

f(x)f(x)为一个正比于目标概率分布p(x)p(x)的近似分布,梅特罗波利斯算法的基本步骤如下:

  1. 初始化阶段:选取一个任意点x0x_0作为起始采样点,并选择一个任意概率密度函数g(xy)g(x|y)用来根据当前采样随机数yy产生下一个候选随机数xx,在梅特罗波利斯算法中,gg必须是对称的,即满足g(xy)=g(yx)g(x|y)=g(y|x)。一个常用的选择是高斯分布(Gaussian distribution),也即正态分布(normal distribution),这种情况下,靠近yy的点会被以更高的概率成为候选采样点。gg被称为建议密度函数(proposal density function)或跳跃分布(jumping distribution)。
  2. 迭代阶段,对于每一次迭代nn: 3.1. 首先从建议密度函数g(xxn)g(x^{'}|x_n)产生一个候选采样点xx^{'}。 3.2 计算接受率α=f(x)f(xn)\alpha= \cfrac{f(x^{'})}{f(x_n)},这个接受率用于决定新的候选采样点被接受或拒绝。因为f(x)f(x)正比于目标概率密度函数p(x)p(x),所以α=f(x)f(xn)=p(x)p(xn)\alpha= \cfrac{f(x^{'})}{f(x_n)}= \cfrac{p(x^{'})}{p(x_n)}。 3.3 如果α1\alpha\geq 1,新的候选点直接被接受:xn+1=xx_{n+1}=x^{'};否则,按照概率α\alpha来接受新的候选点,如果新的候选点被拒绝,则设置:xn+1=xnx_{n+1}=x_n

梅特罗波利斯算法的过程是通过在采样空间随机行走实现,在随机行走的每一步,新的随机采样点可能被接受,也可能被拒绝,这完全由接受率α\alpha来决定。接受率α\alpha表述的是根据目标概率分布p(x)p(x),新产生的候选随机数相对于当前随机数有多大的概率被接受,例如如果α1\alpha\geq 1,说明新的候选采样点xx^{'}p(x)p(x)中处于更高概率区域,即是它比xnx_n具有更高的概率,因此新的候选点总是被接受;如果α<1\alpha<1,则说明新的候选点在p(x)p(x)中处于较低概率区域,此时它可能会被接受,可能会被拒绝,α\alpha越小,新的候选点拥有更小的几率被接受。因此,通过这样的机制,马尔可夫链中的随机数倾向于停留在目标概率分布p(x)p(x)的高概率区域,而只有偶尔很少的机会会访问低概率区域。直观上看,这就是梅特罗波利斯算法能够工作的原因,也使它最终返回的采样点是符合目标概率分布的。

收敛性证明

梅特罗波利斯算法的目标是根据目标概率分布p(x)p(x)产生一个状态序列,为了达到这个目的,它使用一个马尔可夫过程来渐渐地逼近一个唯一的稳定分布π(x)\pi(x),使得π(x)=p(x)\pi(x)=p(x)。那么怎样证明该算法是收敛于p(x)p(x)的呢?

马尔可夫链蒙特卡洛方法的各种算法的收敛性是基于马尔可夫链的收敛性理论。一个马尔可夫过程由其转移概率P(xx)P(x^{'}|x)唯一决定,要使该马尔可夫过程存在一个唯一的稳定分布π(x)\pi(x),必须满足以下两个条件:

  1. 存在稳定分布。
  2. 稳定分布是唯一的。

要使马尔可夫链存在稳定分布,由前面的内容可知,一个充分不必要条件是满足细致平衡(detailed balance)条件,细致平衡条件要求每个转移过程是可逆的,对于每一对状态xxxx^{'},处于状态xx并转移到状态xx^{'}的概率必须等于处于状态xx^{'}并转移到状态xx的概率,即:π(x)P(xx)=π(x)P(xx)\pi(x)P(x^{'}|x)=\pi(x^{'})P(x|x^{'})

稳定分布的唯一性可以由遍历性(ergodicity)来保证,通过前面的马尔可夫链属性可知,马尔可夫链的遍历性必须要求非周期性和正回返性。非周期性其实暗指每个状态不是只能从某些特定路径返回,而是可以沿任意路径返回,其意义就是每个状态可以与任意状态联通;而正回返性保证这种联通的概率小于\infty。由不可约性可知,如果状态空间内所有状态都是联通的,从每个状态出发经过一定步数之后都可以达到另一个状态,则马尔可夫链就是与初始状态无关的,所以遍历性保证了稳定分布的唯一性。

梅特罗波利斯算法就是通过设计一个满足上述两个条件的马尔可夫链,以使稳定分布π\pi逼近p(x)p(x)。该算法的推导首先从满足细致平衡条件开始:

P(xx)p(x)=P(xx)p(x) P(x^{'}|x)p(x)=P(x|x^{'})p(x^{'})

(式15

上式可以改写为:

P(xx)P(xx)=p(x)p(x) \cfrac{P(x^{'}|x)}{P(x|x^{'})}= \cfrac{p(x^{'})}{p(x)}

(式16

梅特罗波利斯算法将转移概率的计算分成两步:建议概率g(xx)g(x^{'}|x)和接受率α(xx)\alpha(x^{'}|x),并使得:

P(xx)=g(xx)α(xx) P(x^{'}|x)=g(x^{'}|x)\alpha(x^{'}|x)

(式17

将上式代入前面的细致平衡关系,得到:

α(xx)α(xx)=p(x)p(x)g(xx)g(xx) \cfrac{\alpha(x^{'}|x)}{\alpha(x|x^{'})}= \cfrac{p(x^{'})}{p(x)} \cfrac{g(x|x^{'})}{g(x^{'}|x)}

(式18

接下来就是怎么选择接受概率α\alpha以满足上面的条件,由前面的梅特罗波利斯算法可知其选择为:

α(xx)=min(1,p(x)p(x)g(xx)g(xx))=min(1,p(x)p(x)) \alpha(x^{'}|x)=\min \bigg( 1, \cfrac{p(x^{'})}{p(x)} \cfrac{g(x|x^{'})}{g(x^{'}|x)}\bigg)=\min\bigg(1, \cfrac{p(x^{'})}{p(x)}\bigg)

(式19

上式后半部分的等式成立是因为梅特罗波利斯算法中建议概率gg是对称的,从此也可以看出,梅特罗波利斯算法中中的建议概率gg是可以任意选择的,因为它并不会参与到最终转移概率的计算上。梅特罗波利斯算法中的转移概率是满足细致平衡条件的,而细致平衡条件可以保证稳定分布,也就是总体平衡,因此梅特罗波利斯算法是收敛的。

细致平衡条件

对于好奇的读者可以由以下证明过程证明其满足细致平衡条件:

p(x)P(yx)=p(x)g(yx)α(yx)=p(x)g(yx)min(1,p(y)/p(x))=min(p(x)g(yx),g(yx)p(y))=p(y)g(xy)min(p(x)/p(y),1)=p(y)g(xy)α(xy)=p(y)P(xy) \begin{aligned} p(x)P(y|x)&=p(x)g(y|x)\alpha(y|x)=p(x)g(y|x)\min(1,p(y)/p(x))\\ &=\min(p(x)g(y|x),g(y|x)p(y))=p(y)g(x|y)\min(p(x)/p(y),1)\\ &=p(y)g(x|y)\alpha(x|y)=p(y)P(x|y) \end{aligned}

(式20

梅特罗波利斯算法的不足

尽管梅特罗波利斯算法中是收敛的,但是马尔可夫链蒙特卡洛方法的各种算法收敛性证明只能证明其样本值遵从已知概率分布,但不能证明样本值是独立无关的,也没有给出收敛速度是多少。本节介绍梅特罗波利斯算法主要的两个方面的不足。

首先,尽管梅特罗波利斯算法的建议概率gg可以是任意选取的,并且能够保证收敛性,但是它的选取却对效率和方差由一定的影响。直观上,为了使马尔可夫链收敛的速度更快,以及更快地遍历整个状态空间,我们趋向于选取步幅更大的建议概率。但是对于用于高频率分布的目标分布函数,过大的建议步幅在高频部分将导致p(x)p(x)p(x^{'})\gg p(x),因此使得接受率α\alpha非常低,其结果就是大量的随机数停留在尖锐的区域,导致方差增大,并且收敛速度变慢(随机数被卡在高频部分很长时间),如图(2)所示;相反,过小的建议步幅则会直接导致更慢的收敛速度,所以合适地选取建议概率函数非常重要,我们将在后面第(7)章学习MLT算法的时候更具体地讨论一些建议概率函数选择的策略。

图(2):如果分布函数的某些区域比较尖锐,或者由于建议概率产生的候选随机数的步幅过大,则会导致大部分候选随机数由于接受率α\alpha太低而被拒绝,因此在尖锐的部分聚集大量相同的随机数,导致方差增大

其次,虽然梅特罗波利斯算法最终会收敛于目标概率分布,但是在这个收敛的过程中,一些初始的随机数序列则可能完全不是服从目标概率分布的,这部分随机数通常需要从马尔可夫链中去除掉,然而并没有一些好的算法用于预测这个收敛的长度需要多少初始随机数。

尽管如此,梅特罗波利斯算法能够轻松应对高维概率分布的采样是其他方法不可比拟的,有时候几乎梅特罗波利斯算法是唯一的选择。在本章中我们仅介绍梅特罗波利斯算法的一些基本概念,在第(7)章还会有更多更具体的内容来讨论梅特罗波利斯算法的一些取舍以及更多的实现细节。