3.2　蒙特卡洛采样

在3.1节中，通过求解定积分这个具体的例子来演示了蒙特卡洛方法的基本思想，而其中的核心就是使用随机数。当所求解问题可以转化为某种随机分布的特征数（如随机事件出现的概率，或者随机变量的期望值等）时，往往就可以考虑使用蒙特卡洛方法。通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解。这种方法多用于求解复杂的高维积分问题。

实际应用中，所要面对的第一个问题就是如何抽样？注意，在计算机模拟时，这里所说的抽样其实是指从一个概率分布中生成观察值（observations）的方法。而这个分布通常是由其概率密度函数来表示的。前面曾经提过，即使在已知PDF的情况下，让计算机自动生成观测值也不是一件容易的事情。从本质上来说，计算机只能实现对均匀分布的采样。幸运的是，仍然可以在此基础上对更为复杂的分布进行采样。

3.2.1　逆采样

比较简单的一种情况是，可以通过PDF与CDF之间的关系，求出相应的CDF。或者根本就不知道PDF，但是知道CDF。此时就可以使用CDF反函数（及分位数函数）的方法来进行采样。这种方法又称为逆变换采样（Inverse Transform Sampling）。

假设已经得到了CDF的反函数F^-1（u），如果想得到m个观察值，则重复下面的步骤m次：

（1）从U（0，1）中随机生成一个值（计算机可以实现从均匀分布中采样），用u表示；

（2）计算F^-1（u）的值x，则x就是从目标分布f（x）中得出的一个采样点。

面对一个具有复杂表达式的函数，逆变换采样法真有效吗？来看一个例子，假设现在希望从具有下面这个PDF的分布中采样

可以算得相应的CDF为

对于u∈［0，1］，上述CDF的反函数为

下面在R中利用上面的方法采样10000个点，并以此来演示抽样的效果。

将所得结果与真实的PDF函数图形进行对照，如图3-5所示。可见由逆变换采样法得到的点所呈现处理的分布与目标分布非常吻合。

下面再举一个稍微复杂一点的例子，已知分布的PDF如下

可以算得相应的CDF为

对于u∈［0，1］，它的反函数为

同样，给出R中的示例代码如下。

下面这段代码利用R中提供的一些内置函数实现了已知PDF时基于逆变换方法的采样，将新定义的函数命名为samplepdf（）。当然，对于那些过于复杂的PDF函数（例如很难积分的），samplepdf（）确实有力所不及的情况。但是对于标准的常规PDF，该函数的效果还是不错的。

下面就用samplepdf（）函数来对上面给定的h（x）进行采样，然后再跟之前所得结果进行对比。

代码执行结果如图3-6所示。

3.2.2　博克斯-穆勒变换

博克斯-穆勒变换（Box-Muller Transform）最初由乔治·博克斯（George Box）与默文·穆勒（Mervin Muller）在1958年共同提出。博克斯是统计学的一代大师，统计学中的很多名词术语都以其名字命名。博克斯与统计学的家学渊源相当深厚，他的导师是统计学开山鼻祖皮尔逊的儿子，英国统计学家埃贡·皮尔逊（Egon Pearson），博克斯还是统计学的另外一位巨擘级奠基人费希尔的女婿。统计学中的名言“所有模型都是错的，但其中一些是有用的”也出自博克斯之口。

本质上来说，计算机只能生产符合均匀分布的采样。如果要生成其他分布的采样，就需要借助一些技巧性的方法。而在众多的“其他分布”中，正态分布无疑占据着相当重要的地位。下面这个定理，就为生成符合正态分布的采样（随机数）提供了一种方法，而且这也是很多软件或者编程语言的库函数中生成正态分布随机数时所采样的方法。

定理（Box-Muller变换）：如果随机变量U₁和U₂是独立同分布的，且U₁，U₂~U［0，1］，则

其中，Z₀和Z₁独立且服从标准正态分布。

如何来证明这个定理呢？这需要用到一些微积分中的知识，首先回忆一下二重积分化为极坐标下累次积分的方法

假设现在有两个独立的标准正态分布X~N（0，1）和Y~N（0，1），由于两者相互独立，则联合概率密度函数为

做极坐标变换，则x=Rcosθ，y=Rsinθ，则有

这个结果可以看成是两个概率分布的密度函数的乘积，其中一个可以看成是［0，2π］上均匀分布，将其转换为标准均匀分布则有θ~U（0，2π）=2πU₂。

另外一个的密度函数为

则其累计分布函数CDF为

这个CDF函数的反函数可以写成

根据逆变换采样的原理，如果有个PDF为P（R）的分布，那么对齐CDF的反函数进行均匀采样所得的样本分布将符合P（R）的分布，而如果u是均匀分布的，那么U₁=1-u也将是均匀分布的，于是用U₁替换1-u，最后可得

结论得证。最后来总结一下利用Box-Muller变换生成符合高斯分布的随机数的方法：

（1）产生两个随机数U₁，U₂~U［0，1］；

（2）用它们来创造半径和和夹角θ=2πU₂；

（3）将（R，θ）从极坐标转换到笛卡儿坐标：（Rcosθ，Rsinθ）。

3.2.3　拒绝采样与自适应拒绝采样

读者已经看到逆变换采样的方法确实有效。但其实它的缺点也是很明显的，那就是有些分布的CDF可能很难通过对PDF的积分得到，再或者CDF的反函数也很不容易求。这时可能需要用到另外一种采样方法，这就是下面即将要介绍的拒绝采样（Reject Sampling）。

图3-7阐释了拒绝采样的基本思想。假设对PDF为p（x）的函数进行采样，但是由于种种原因（例如这个函数很复杂），对其进行采样是相对困难的。但是另外一个PDF为q（x）的函数则相对容易采样，例如采用逆变换方法可以很容易对对它进行采样，甚至q（x）就是一个均匀分布（别忘了计算机可以直接进行采样的分布就只有均匀分布）。那么，当我们将q（x）与一个常数M相乘之后，可以实现如图3-7所示关系，即M·q（x）将p（x）完全“罩住”。

然后重复如下步骤，直到获得m个被接受的采样点：

（1）从q（x）中获得一个随机采样点x_i；

（2）对于x_i计算接受概率（acceptanceprobability）

（3）从U（0，1）中随机生成一个值，用u表示；

（4）如果α≥u，则接受x_i作为一个来自p（x）的采样值，否则就拒绝x_i并回到第一步。

当然可以采用严密的数学推导来证明拒绝采样的可行性。但它的原理从直观上来解释也是相当容易理解的。可以想象一下在图2-7的例子中，从哪些位置抽出的点会比较容易被接受。显然，红色曲线（位于上方）和绿色曲线（位于下方）所示之函数更加接近的地方接受概率较高，即是更容易被接受，所以在这样的地方采到的点就会比较多，而在接受概率较低（即两个函数差距较大）的地方采到的点会比较少，这也就保证了这个方法的有效性。

还是以本章前面给出的那个分段函数f（x）为例来演示拒绝采样方法。如图3-8所示，所选择的参考分布是均匀分布（当然也可以选择其他的分布，但采用均匀分布显然是此处最简单的一种处理方式）。而且令常数M=3。

下面给出R中的示例代码，易见此处采样点数目为10 000。

上述代码的执行结果如图3-9所示，可见采样结果是非常理想的。

下面的例子演示了对（表达式非常复杂的）beta（3，6）分布进行拒绝采样的效果。这里采用均匀分布作为参考分布。而且这里的Mq（x）所取之值就是beta（3，6）分布的极大值，它的函数图形应该是与beta（3，6）的极值点相切的一条水平直线。

图3-10给出了采样50 000个点后的密度分布情况，可见采样分布与目标分布beta（3，6）非常吻合。

拒绝采样的方法确实可以解决我们的问题。但是它的一个不足涉及其采样效率的问题。针对上面给出的例子而言，我们选择了离目标函数最近的参考函数，就均匀分布而言，已经不能有更进一步的方法了。但即使这种，在这个类似钟形的图形两侧其实仍然会拒绝掉很多很多采样点，这种开销相当浪费。最理想的情况下，参考分布应该跟目标分布越接近越好，从图形上来看就是包裹的越紧实越好。但是这种情况的参考分布往往又不那么容易得到。在满足某些条件的时候也确实可以采用所谓的改进方法，即自适应的拒绝采样（Adaptive Rejection Sampling）。

拒绝采样的弱点在于当被拒绝的点很多时，采样的效率会非常不理想。同时我们也知道，如果能够找到一个跟目标分布函数非常接近的参考函数，那么就可以保证被接受的点占大多数（被拒绝的点很少）。这样一来便克服了拒绝采样效率不高的弱点。如果函数是log-concave的话，那么就可以采用自适应的拒绝采样方法。什么是log-concave呢？还是回到之前介绍过的贝塔分布，用下面的代码来绘制beta（2，3）的概率密度函数图像，以及将beta（2，3）的函数取对数之后的图形。

上述代码的执行结果如图3-11所示，其中（a）是beta（2，3）的概率密度函数图形，（b）是将beta（2，3）的函数取对数之后的图形，你可以发现结果是一个凹函数（concave）。那么beta（2，3）就满足log-concave的要求。

然后在对数图像上找一些点做图像的切线，如图3-12所示。因为取对数后的函数是凹函数，所以每个切线都相当于一个超平面，而且对数图像只会位于超平面的一侧。

同时给出用以绘制图3-12的代码，要知道R语言的一个强项就是绘图。

再把这些切线转换回原始的beta（2，3）图像中，显然原来的线性函数会变成指数函数，它们将对应图3-13中的一些曲线，这些曲线会被原函数的图形紧紧包裹住。特别是当这些的指数函数变得很多很稠密时，以彼此的交点作为分界线，其实相当于得到了一个分段函数。这个分段函数是原函数的一个逼近。用这个分段函数来作为参考函数再执行拒绝采样，自然就完美地解决了之前的问题。

下面是用来画出图3-13的R语言代码。

这无疑是一种绝妙的想法。而且这种想法，在前面其实已经暗示过。在上一部分最后一个例子中，我们其实就是选择了一个与原函数相切的均匀分布函数来作为参考函数。我们当然会想去选择更多与原函数相切的函数，然后用这个函数的集合来作为新的参考函数。只是由于原函数的凹凸性无法保证，所以直线并不是一种好的选择。而自适应拒绝采样（Adaptive Rejection Sampling，ARS）所采用的策略则非常巧妙地解决了我们的问题。当然函数是log-concave的条件必须满足，否则就不能使用ARS。

下面给出一个在R中进行自适应拒绝采样的例子。显然，该例子要比之前的代码简单许多。因为R中ars包已经提供了一个现成的用于执行自适应拒绝采样的函数，即ars（）。关于这个函数在用法上的一些细节，读者还可以进一步参阅R的帮助文档，这里不再赘言。此次我们需要指出：ars（）函数中两个重要参数，一个是对原分布的PDF取对数，另外一个则是对PDF的对数形式再进行求导（在求导时我们忽略了前面的系数项），其实也就是为了确定切线。

上述代码的执行结果如图3-14所示。