随机采样：基于Markov链的优化方法

Markov链的一个重要应用是：产生服从给定概率分布w的样本。假设我们只能对服从均匀分布的随机变量进行采样，也就是说，1)每次都随机地在1,2,3,···,N中选取一个数，每个数被选取的概率都一样，或者，2)每次都从区间[0,1]中随机选取一个数，相应的概率分布为均匀分布。我们的目标是：通过控制上述两种采样过程，来生成一些服从给定概率分布w的样本。

事实上，式(15.19)是一个算法迭代过程：

每一步的迭代过程为式(15.16)。我们要解决的第一个问题是：将每一步的迭代过程转换为一个采样过程，假设第n步时的样本为xk，在第n+1步时，首先，通过均匀分布的随机采样得到新的样本xl，然后，采用“接受/拒绝”采样进行筛选，判断是否更新xk。具体过程为：首先，生成一个随机数u∈[0,1]，如果u＜αlk，则将样本更新为xl；反之，样本不变，仍为xk。样本由xk更新为xl的条件概率为plk，因此，采样过程应该通过条件u＜plk来控制是否将xk更新为xl。需要注意的是，不将xk更新为xl，并不等于说xk继续保持不变，还包括将xk更新为其他xj的情况。“xk继续保持不变”的概率pkk远小于“xk不被更新为xl”的概率1plk。一种修正方法为：将plk做同比例放大，也就是说，

但是此时转移矩阵P发生了变化，无法保障其平衡状态仍然是w。另一种方法是：等比例地保留“xk持续保持不变”过程中采集到的m个样本中的pkkm/(1pk l)个样本，而不是保留全部m个样本xk。于是，大量的样本被舍弃掉了。此外，对于一些极端情况，例如式(15.29)和(15.31)中的矩阵P，由于大量的pkk=0，随机采样将难以进行下去。

值得庆幸的是：我们的任务是根据平衡状态w来设计转移矩阵P，进而实现式(15.49)所对应的采样过程。我们在设计P的过程中，要将上述问题考虑进去，这也是Hasting对Metropolis 1954年提出的原创算法进行改进的根本出发点。

我们现在来介绍Metropolis-Hasting采样算法（简称M-H采样），也称为Markov链Monte Carlo方法^[2]（简称MCMC方法）。正如我们前面所谈到的，对于平衡状态w，每一个节点的流入总量等于流出总量，参见式(15.20)。当然，这并不等价于：沿着每条边的流入量等于沿着相同的边的流出量，但是，如果沿着每条边的流入量等于沿着相同边的流出量，那么每一个节点的流入总量一定等于流出总量。这为我们设计转移矩阵P提供了依据，具体地说，就是令：

上式中，wi和wj是已知的，而pj,i和pj,i是未知的。我们只需令：

即可使式(15.50)成立，其中，ci,j＞0是任意常数。此时，pj,i已经进行了比例放大，因此，不再是条件概率。在比例放大的过程中，始终保障了“流入等于流出”，因此，平衡状态w始终保持不变。

剩下的一个问题是：ci,j应该取多大。正如式(15.49)所描述的，算法中将“xk不更新为xl”替换成了“xk继续保持不变”，使得pkk被放大成了1plk，因此，应该对pi,j做相应的放大，但是，采样过程是通过对比随机数u∈[0,1]和pi,j=ci,jwj（或pj,i=ci,jwi）来进行的，因此，pi,j=ci,jwj和pj,i=ci,jwi中的一个最大可以被放大到1，也就是说，ci,j应该选为^[3]

此时，相应的pj,i和pi,j分别为：

于是，我们得到了（基于均匀分布采样的）MCMC方法（也称为M-H采样）：

1.随机选取初始值x(0)=xk。

2.进行迭代采样，令n=1,2,3,···(www.xing528.com)

(a)假设x(n)=xi，通过均匀分布随机采样，得到一个1到N之间的正整数j。令xj作为x(n+1)的备选采样结果。

(b)生成一个随机数u∈[0,1]，如果u＜pj,i，则令x(n+1)=xj，其中pj,i=m in{wj/wi,1}（参见式(15.2)）。

(c)如果u≤pj,i，就令x(n+1)=xi。

3.选取一个足够大的m，将x(m)后续生成的M个采样结果x(m+1),x(m+2),x(m+3),···,x(m+M)作为最终的样本集。

图15.2给出了一个具体的例子。我们令xk=k，其中k=1,2,3,···,11。假设随机变量X n服从θ=0.4的二项式分布，也就是说，P(Xn=xk)=，如图15.2(c)中的“黑圈”所示。采用上面介绍的（基于均匀分布采样的）MCMC方法进行20000次采样，结果（包含20000个样本点）如图15.2(a)所示。我们分别对：1)第5000次到第10000次的采样结果（共5001个样本点）和2)第10000次到第20000次的采样结果（共10001个样本点）做统计，所得到的统计直方图如图15.2(b)所示。进一步对图15.2(b)中的统计结果做“归一化”处理，所得到的统计分布如图15.2(c)所示。可以看到，统计分布非常接近随机变量Xn的概率分布P(Xn=k)=（其中k=1,2,3,···,11），说明了M-H采样方法的有效性。

图15.3给出了另外一个实验结果。我们令xk=k，其中k=1,2,3,···,40，随机生成一个概率分布，如图15.3(c)中的“黑圈”所示。采用上面介绍的（基于均匀分布采样的）MCMC方法进行40000次采样，结果（包含40000个样本点）如图15.3(a)所示。我们分别对：1)第5000次到第10000次的采样结果（共5001个样本点）和2)第10000次到第20000次的采样结果（共10001个样本点）做统计，所得到的统计直方图如图15.3(b)所示。进一步，对图15.3(b)中的统计结果做“归一化”处理，所得到的统计分布如图15.3(c)所示。可以看到，统计分布非常接近随机变量Xn的概率分布。这个实验进一步验证了M-H采样方法的有效性。

虽然M-H采样算法在理论上存在一些缺陷：将“xk不被更新为xl”替换成“xk继续保持不变”的过程中，忽略了“xk被更新为其他xi”的情况，但是，MCMC方法有效地实现了随机采样。图15.4所示的是：直接采用“接受/拒绝”采样方法所得到的结果。

图15.2　我们令xk=k，其中k=1,2,3,···,11，将随机变量X n的概率分布设置为一个二项式分布，也就是说，P(X n=xk)= 0.4k-1 0.6(n+1-k)。(a)采用（基于均匀分布采样）的MCMC方法进行20000次采样，所得到的20000个样本点。(b)分别对：1)第5000次到第10000次的采样结果（共5001个样本点）和2)第10000次到第20000次的采样结果（共10001个样本点）做统计，所得到的统计直方图。(c)对统计直方图做“归一化”处理，所得到的统计分布。这两个统计分布非常接近随机变量的概率分布。

为了与图15.3进行对比，我们也令xk=k，其中k=1,2,3,···,40，然后，随机生成一个概率分布，如图15.4(c)中的“黑圈”所示。直接采用“接受/拒绝”采样方法进行40000次实验，每次实验都随机生成：一个1到40之间的整数k和一个0到1之间的随机数u，如果u＜P(Xn=k)，则保留该采样结果，否则，直接舍弃该采样结果。在40000次实验中，共舍弃了38989次采样结果（约97.47%），只保留了1011个采样结果（约2.53%），图15.4(a)中给出了相应的实验结果（即保留下来的1011个样本点）。对“接受/拒绝”采样方法所得到的1011个样本点做统计，所得到的统计直方图如图15.4(b)所示。进一步，对图15.4(b)中的统计结果做“归一化”处理，所得到的统计分布如图15.4(c)所示。可以看到，由于有效样本点数目不足（共1011个，约占采样次数的2.53%），统计分布与概率分布之间存在一定的偏差，参见图15.4(c)。

图15.3　我们令xk=k，其中k=1,2,3,···,40，随机生成一个概率分布，然后，采用（基于均匀分布采样）的MCMC方法进行随机采样。(a)MCMC采样所得到的40000个样本点。(b)分别对：1)第5000次到第10000次的采样结果（共5001个样本点）和2)第30000次到第40000次的采样结果（共10001个样本点）做统计，所得到的统计直方图。(c)对统计直方图做“归一化”处理，所得到的统计分布。这两个统计分布非常接近随机变量的概率分布。

随机变量的取值xk可以是任意的（甚至可以是一张图像）。此时，我们仍然对k=1,2,3,···,N进行随机采样，再根据样本点k找到相应的xk。通过随机采样的方式，我们可以根据一个小的样本集生成一个大的数据集，例如：

图15.4　直接采用“接受/拒绝”方法得到的采样结果。经过40000次采样共保留了1011个有效样本（约占采样次数的2.53%），由于有效样本数目不足，统计分布与概率分布之间存在一定偏差。(a)在40000次实验中，共舍弃了38989次采样结果，只保留了1011个有效样本。(b)对“接受/拒绝”采样方法所得到的1011个样本点做统计，所得到的统计直方图。(c)根据1011个有效采样结果所得到的统计分布，对比给定的概率分布，两者之间存在一定的偏差。

当然，我们也可以直接对小样本集做统计，然后，根据统计直方图来进行MCMC采样，从而生成大的样本集。

需要指出的是，Markov链的应用远不止于此。此外，我们还可以通过设置停止时间，来对Markov过程进行优化控制。在习题15.1中，我们详细讨论了一个基于Markov链的最优控制策略问题。在后续章节中，我们还会介绍一些相关的实际应用案例。