蒙特卡罗方法课件(清华大学林谦)

蒙特卡罗方法

在核技术中的应用林谦2021/3/271CHENLI目录

第一章蒙特卡罗方法概述第二章随机数第三章由已知分布的随机抽样第四章蒙特卡罗方法解粒子输运问题2021/3/272CHENLI教材蒙特卡罗方法在实验核物理中的应用 许淑艳编著 原子能出版社蒙特卡罗方法 清华大学2021/3/273CHENLI参考书蒙特卡罗方法及其在粒子输运问题中的应用 裴鹿成张孝泽编著 科学出版社蒙特卡罗方法 徐钟济编著 上海科学技术出版社2021/3/274CHENLI联系方式电话 83918

电子邮件 linqian@

2021/3/275CHENLI第一章蒙特卡罗方法概述蒙特卡罗方法的基本思想蒙特卡罗方法的收敛性，误差蒙特卡罗方法的特点蒙特卡罗方法的主要应用范围作业2021/3/276CHENLI第一章蒙特卡罗方法概述

蒙特卡罗方法又称随机抽样技巧或统计试验方法。半个多世纪以来，由于科学技术的发展和电子计算机的发明，这种方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。2021/3/277CHENLI蒙特卡罗方法的基本思想

二十世纪四十年代中期，由于科学技术的发展和电子计算机的发明，蒙特卡罗方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。但其基本思想并非新颖，人们在生产实践和科学试验中就已发现，并加以利用。两个例子

例1.蒲丰氏问题

例2.射击问题（打靶游戏）基本思想计算机模拟试验过程2021/3/278CHENLI例1.蒲丰氏问题

为了求得圆周率π值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（

l＜a）的平行线相交的频率代替概率P，再利用准确的关系式：求出π值其中Ｎ为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。2021/3/279CHENLI

一些人进行了实验，其结果列于下表：实验者年份投计次数π的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)1901340833/2710CHENLI例2.射击问题（打靶游戏）

设r表示射击运动员的弹着点到靶心的距离，ｇ(r)表示击中r处相应的得分数（环数），f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射击水平。该运动员的射击成绩为

用概率语言来说，<g>是随机变量ｇ(r)的数学期望，即

2021/3/2711CHENLI

现假设该运动员进行了Ｎ次射击，每次射击的弹着点依次为r1，r2，…，rN，则Ｎ次得分g(r1)，g(r2)，…，g(rN)的算术平均值代表了该运动员的成绩。换言之，为积分<g>的估计值，或近似值。在该例中，用Ｎ次试验所得成绩的算术平均值作为数学期望<g>的估计值（积分近似值）。

2021/3/2712CHENLI基本思想

由以上两个例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。当随机变量的取值仅为1或0时，它的数学期望就是某个事件的概率。或者说，某种事件的概率也是随机变量（仅取值为1或0）的数学期望。

2021/3/2713CHENLI

因此，可以通俗地说，蒙特卡罗方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数f(r)的随机变量ｇ(r)的数学期望

通过某种试验，得到Ｎ个观察值r1，r2，…，rN（用概率语言来说，从分布密度函数f(r)中抽取Ｎ个子样r1，r2，…，rN，），将相应的Ｎ个随机变量的值g(r1)，g(r2)，…，g(rN)的算术平均值作为积分的估计值（近似值）。

2021/3/2714CHENLI

为了得到具有一定精确度的近似解，所需试验的次数是很多的，通过人工方法作大量的试验相当困难，甚至是不可能的。因此，蒙特卡罗方法的基本思想虽然早已被人们提出，却很少被使用。本世纪四十年代以来，由于电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验交由计算机完成，使得蒙特卡罗方法得以广泛地应用，在现代化的科学技术中发挥应有的作用。

2021/3/2715CHENLI计算机模拟试验过程

计算机模拟试验过程，就是将试验过程（如投针，射击）化为数学问题，在计算机上实现。以上述两个问题为例，分别加以说明。例1.蒲丰氏问题例2.射击问题（打靶游戏）由上面两个例题看出，蒙特卡罗方法常以一个“概率模型”为基础，按照它所描述的过程，使用由已知分布抽样的方法，得到部分试验结果的观察值，求得问题的近似解。2021/3/2716CHENLI例１．蒲丰氏问题

设针投到地面上的位置可以用一组参数（x,θ）来描述，x为针中心的坐标，θ为针与平行线的夹角，如图所示。任意投针，就是意味着x与θ都是任意取的，但x的范围限于［0，a］，夹角θ的范围限于［0，π］。在此情况下，针与平行线相交的数学条件是针在平行线间的位置2021/3/2717CHENLI

如何产生任意的（x,θ）？x在［0，a］上任意取值，表示x在［0，a］上是均匀分布的，其分布密度函数为：类似地，θ的分布密度函数为：因此，产生任意的（x,θ）的过程就变成了由f1(x)抽样x及由f2(θ)抽样θ的过程了。由此得到：其中ξ1，ξ2均为（0,1）上均匀分布的随机变量。2021/3/2718CHENLI

每次投针试验，实际上变成在计算机上从两个均匀分布的随机变量中抽样得到（x,θ），然后定义描述针与平行线相交状况的随机变量s(x,θ)，为如果投针Ｎ次，则是针与平行线相交概率Ｐ的估计值。事实上，于是有2021/3/2719CHENLI例２．射击问题

设射击运动员的弹着点分布为用计算机作随机试验（射击）的方法为，选取一个随机数ξ，按右边所列方法判断得到成绩。这样，就进行了一次随机试验（射击），得到了一次成绩ｇ(r)，作Ｎ次试验后，得到该运动员射击成绩的近似值环数78910概率0.10.10.30.52021/3/2720CHENLI蒙特卡罗方法的收敛性，误差

蒙特卡罗方法作为一种计算方法，其收敛性与误差是普遍关心的一个重要问题。收敛性误差减小方差的各种技巧效率2021/3/2721CHENLI收敛性

由前面介绍可知，蒙特卡罗方法是由随机变量X的简单子样X1，X2，…，XN的算术平均值：作为所求解的近似值。由大数定律可知，如X1，X2，…，XN独立同分布，且具有有限期望值（E(X)<∞），则即随机变量X的简单子样的算术平均值，当子样数Ｎ充分大时，以概率1收敛于它的期望值E(X)。2021/3/2722CHENLI误差

蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，…，XN独立同分布，且具有有限非零的方差σ2

，即

f(X)是X的分布密度函数。则2021/3/2723CHENLI

当N充分大时，有如下的近似式其中α称为置信度，1－α称为置信水平。这表明，不等式近似地以概率

1－α成立，且误差收敛速度的阶为。通常，蒙特卡罗方法的误差ε定义为上式中与置信度α是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出。2021/3/2724CHENLI

下面给出几个常用的α与的数值：

关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。第二，误差中的均方差σ是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出。α0.50.050.003

0.67451.9632021/3/2725CHENLI减小方差的各种技巧

显然，当给定置信度α后，误差ε由σ和N决定。要减小ε，或者是增大N，或者是减小方差σ2。在σ固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。另一方面，如能减小估计的均方差σ，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果。因此降低方差的各种技巧，引起了人们的普遍注意。后面课程将会介绍一些降低方差的技巧。2021/3/2726CHENLI效率

一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中效率的概念。它定义为，其中c

是观察一个子样的平均费用。显然越小，方法越有效。2021/3/2727CHENLI蒙特卡罗方法的特点优点能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。受几何条件限制小。收敛速度与问题的维数无关。具有同时计算多个方案与多个未知量的能力。误差容易确定。程序结构简单，易于实现。缺点收敛速度慢。误差具有概率性。在粒子输运问题中，计算结果与系统大小有关。2021/3/2728CHENLI能够比较逼真地描述具有随机性质的事物的特点及物理实验过程

从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。2021/3/2729CHENLI受几何条件限制小

在计算s维空间中的任一区域Ds上的积分时，无论区域Ds的形状多么特殊，只要能给出描述Ds的几何特征的条件，就可以从Ds中均匀产生N个点，得到积分的近似值。其中Ds为区域Ds的体积。这是数值方法难以作到的。另外，在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。2021/3/2730CHENLI收敛速度与问题的维数无关

由误差定义可知，在给定置信水平情况下，蒙特卡罗方法的收敛速度为，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计算时间的变化，不影响误差。也就是说，使用蒙特卡罗方法时，抽取的子样总数N与维数s无关。维数的增加，除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。2021/3/2731CHENLI具有同时计算多个方案与多个未知量的能力

对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。另外，使用蒙特卡罗方法还可以同时得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。2021/3/2732CHENLI误差容易确定

对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。一般计算方法常存在着有效位数损失问题，而要解决这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。2021/3/2733CHENLI程序结构简单，易于实现

在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。2021/3/2734CHENLI收敛速度慢

如前所述，蒙特卡罗方法的收敛速度为，一般不容易得到精确度较高的近似结果。对于维数少（三维以下）的问题，不如其他方法好。2021/3/2735CHENLI误差具有概率性

由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。2021/3/2736CHENLI在粒子输运问题中，计算结果与系统大小有关

经验表明，只有当系统的大小与粒子的平均自由程可以相比较时（一般在十个平均自由程左右），蒙特卡罗方法计算的结果较为满意。但对于大系统或小概率事件的计算问题，计算结果往往比真值偏低。而对于大系统，数值方法则是适用的。因此，在使用蒙特卡罗方法时，可以考虑把蒙特卡罗方法与解析（或数值）方法相结合，取长补短，既能解决解析（或数值）方法难以解决的问题，也可以解决单纯使用蒙特卡罗方法难以解决的问题。这样，可以发挥蒙特卡罗方法的特长，使其应用范围更加广泛。2021/3/2737CHENLI蒙特卡罗方法的主要应用范围

蒙特卡罗方法所特有的优点，使得它的应用范围越来越广。它的主要应用范围包括：粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。随着科学技术的发展，其应用范围将更加广泛。蒙特卡罗方法在粒子输运问题中的应用范围主要包括：实验核物理，反应堆物理，高能物理等方面。蒙特卡罗方法在实验核物理中的应用范围主要包括：通量及反应率，中子探测效率，光子探测效率，光子能量沉积谱及响应函数，气体正比计数管反冲质子谱，多次散射与通量衰减修正等方面。2021/3/2738CHENLI作业

用蒲丰投针法在计算机上计算π值，取a=4、l=3。分别用理论计算和计算机模拟计算，求连续掷两颗骰子，点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率。2021/3/2739CHENLI第二章随机数随机数的定义及产生方法伪随机数产生伪随机数的乘同余方法产生伪随机数的乘加同余方法产生伪随机数的其他方法伪随机数序列的均匀性和独立性作业2021/3/2740CHENLI第二章随机数

由具有已知分布的总体中抽取简单子样，在蒙特卡罗方法中占有非常重要的地位。总体和子样的关系，属于一般和个别的关系，或者说属于共性和个性的关系。由具有已知分布的总体中产生简单子样，就是由简单子样中若干个性近似地反映总体的共性。随机数是实现由已知分布抽样的基本量，在由已知分布的抽样过程中，将随机数作为已知量，用适当的数学方法可以由它产生具有任意已知分布的简单子样。2021/3/2741CHENLI随机数的定义及产生方法随机数的定义及性质随机数表物理方法2021/3/2742CHENLI随机数的定义及性质

在连续型随机变量的分布中，最简单而且最基本的分布是单位均匀分布。由该分布抽取的简单子样称，随机数序列，其中每一个体称为随机数。单位均匀分布也称为［0，1］上的均匀分布，其分布密度函数为：分布函数为：2021/3/2743CHENLI

由于随机数在蒙特卡罗方法中占有极其重要的位置，我们用专门的符号ξ表示。由随机数序列的定义可知，ξ1，ξ2，…是相互独立且具有相同单位均匀分布的随机数序列。也就是说，独立性、均匀性是随机数必备的两个特点。随机数具有非常重要的性质：对于任意自然数s，由s个随机数组成的s维空间上的点（ξn+1，ξn+2，…ξn+s）在s维空间的单位立方体Gs上均匀分布，即对任意的ai，如下等式成立：2021/3/2744CHENLI

其中P(·)表示事件·发生的概率。反之，如果随机变量序列ξ1,ξ2…对于任意自然数s，由s个元素所组成的s维空间上的点（ξn+1，…ξn+s）在Gs上均匀分布，则它们是随机数序列。由于随机数在蒙特卡罗方法中所处的特殊地位，它们虽然也属于由具有已知分布的总体中产生简单子样的问题，但就产生方法而言，却有着本质上的差别。2021/3/2745CHENLI随机数表

为了产生随机数，可以使用随机数表。随机数表是由0，1，…，9十个数字组成，每个数字以0.1的等概率出现，数字之间相互独立。这些数字序列叫作随机数字序列。如果要得到n位有效数字的随机数，只需将表中每n个相邻的随机数字合并在一起，且在最高位的前边加上小数点即可。例如，某随机数表的第一行数字为7634258910…，要想得到三位有效数字的随机数依次为0.763，0.425，0.891。因为随机数表需在计算机中占有很大内存，而且也难以满足蒙特卡罗方法对随机数需要量非常大的要求，因此，该方法不适于在计算机上使用。2021/3/2746CHENLI物理方法

用物理方法产生随机数的基本原理是：利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。用来作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。一般情况下，任意一个随机数在计算机内总是用二进制的数表示的：

其中εi（i=1，2，…，m）或者为0，或者为1。2021/3/2747CHENLI

因此，利用物理方法在计算机上产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为0.5。用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，需要增加随机数发生器和电路联系等附加设备，费用昂贵。因此，该方法也不适合在计算机上使用。2021/3/2748CHENLI伪随机数伪随机数伪随机数存在的两个问题伪随机数的周期和最大容量2021/3/2749CHENLI伪随机数

在计算机上产生随机数最实用、最常见的方法是数学方法，即用如下递推公式：产生随机数序列。对于给定的初始值ξ1,ξ2…，ξk，确定ξn+k，ｎ=1，2，…。经常使用的是k=1的情况，其递推公式为：

对于给定的初始值ξ1，确定ξn+1，ｎ=１,２…2021/3/2750CHENLI伪随机数存在的两个问题

用数学方法产生的随机数，存在两个问题：递推公式和初始值ξ1,ξ2…，ξk确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。由于随机数序列是由递推公式确定的，而在计算机上所能表示的[0，1]上的数又是有限的，因此，这种方法产生的随机数序列就不可能不出现无限重复。一旦出现这样的n＇，n″(n＇<n″)，使得下面等式成立：随机数序列便出现了周期性的循环现象。对于k=1的情况，只要有一个随机数重复，其后面的随机数全部重复，这与随机数的要求是不相符的。2021/3/2751CHENLI

由于这两个问题的存在，常称用数学方法产生的随机数为伪随机数。对于以上存在的两个问题，作如下具体分析。关于第一个问题，不能从本质上加以改变，但只要递推公式选得比较好，随机数间的相互独立性是可以近似满足的。至于第二个问题，则不是本质的。因为用蒙特卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所用随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。用数学方法产生的伪随机数容易在计算机上得到，可以进行复算，而且不受计算机型号的限制。因此，这种方法虽然存在着一些问题，但仍然被广泛地在计算机上使用，是在计算机上产生伪随机数的主要方法。2021/3/2752CHENLI伪随机数的周期和最大容量

发生周期性循环现象的伪随机数的个数称为伪随机数的周期。对于前面介绍的情况，伪随机数的周期为n″－n＇。从伪随机数序列的初始值开始，到出现循环现象为止，所产生的伪随机数的个数称为伪随机数的最大容量。前面的例子中，伪随机数的最大容量为n″。2021/3/2753CHENLI产生伪随机数的乘同余方法

乘同余方法是由Lehmer在1951年提出来的，它的一般形式是：对于任一初始值x1，伪随机数序列由下面递推公式确定：其中a为常数。2021/3/2754CHENLI乘同余方法的最大容量的上界

对于任意正整数M，根据数论中的标准分解定理，总可以分解成如下形式：

其中P0=2，P1，…Pr表示不同的奇素数，α0表示非负整数，α1，…，αr表示正整数。a无论取什么值，乘同余方法的最大容量的上界为：

的最小公倍数。其中：2021/3/2755CHENLI关于a与x1的取值

如果a与x1满足如下条件：

对于

，

x1与M互素，则乘同余方法产生的伪随机数序列的最大容量达到最大可能值λ(M)。2021/3/2756CHENLI乘同余方法在计算机上的使用

为了便于在计算机上使用，通常取： Ｍ=2s

其中s为计算机中二进制数的最大可能有效位数

x1=奇数

a=52k+1

其中k为使52k+1在计算机上所能容纳的最大整数，即a为计算机上所能容纳的5的最大奇次幂。一般地，s=32时，a=513；s=48，a=515等。伪随机数序列的最大容量λ(M)=2s-2。

乘同余方法是使用的最多、最广的方法，在计算机上被广泛地使用。2021/3/2757CHENLI产生伪随机数的乘加同余方法

产生伪随机数的乘加同余方法是由Rotenberg于1960年提出来的，由于这个方法有很多优点，已成为仅次于乘同余方法产生伪随机数的另一主要方法。

乘加同余方法的一般形式是，对任意初始值x1，伪随机数序列由下面递推公式确定：其中a和c为常数。2021/3/2758CHENLI乘加同余方法的最大容量

关于乘加同余方法的最大容量问题，有如下结论：如果对于正整数M的所有素数因子P，下式均成立：

当M为4的倍数时，还有下式成立：

c与M互素，则乘加同余方法所产生的伪随机数序列的最大容量达到最大可能值M。2021/3/2759CHENLI

M，x1，a，c的取值

为了便于在计算机上使用，通常取

M=2s

其中s为计算机中二进制数的最大可能有效位数。

a=2b+1 (b≥2) c=1

这样在计算中可以使用移位和指令加法，提高计算速度。

2021/3/2760CHENLI产生伪随机数的其他方法取中方法加同余方法2021/3/2761CHENLI伪随机数序列的均匀性和独立性

判断伪随机数序列是否满足均匀和相互独立的要求，要靠统计检验的方法实现。对于伪随机数的统计检验，一般包括两大类：均匀性检验和独立性检验。六十年代初，人们开始用定性的方法研究伪随机数序列的均匀性和独立性问题，简要叙述如下。2021/3/2762CHENLI伪随机数的均匀性

这里只考虑伪随机数序列ξ1,ξ2…，ξn全体作为子样时的均匀性问题。其中n为伪随机数序列的最大容量。对于任意的0≤x≤1，令Nn(x)表示伪随机数序列ξ1,ξ2…，ξn中适合不等式

ξi<xi=1，2，…，n

的个数，则标志伪随机数序列ξ1,ξ2…，ξn的均匀程度，称为均匀偏度。2021/3/2763CHENLI

将伪随机数序列ξ1,ξ2…，ξn从小至大重新排列并令，则由δ(n)的定义，容易证明很明显，对于固定的ｎ，δ(n)的值越小越好。它是描述伪随机数序列均匀程度的基本量。对于任意随机数序列，均有如下不等式成立：当时，所对应的伪随机数序列为最佳分布。2021/3/2764CHENLI

可以证明，伪随机数序列为最佳分布的充要条件是它取遍序列的所有值。对于计算机上使用的乘同余方法，按照前面介绍的方法选取a、x1时，所产生的伪随机数序列的均匀偏度对于乘加同余方法对于部分伪随机数的均匀性问题通常用统计检验方法检验。2021/3/2765CHENLI伪随机数的独立性

对于任意，令表示(ξ1,ξ2)，(ξ2,ξ3)，…，(ξn,ξn+1)中适合不等式

的个数，根据随机变量间相互独立的定义和频率近似概率的方法，令则ε(n)标志伪随机数序列ξ1,ξ2…，ξn的独立程度，简称为独立偏度。对于固定的n,ε(n)的值越接近于零，伪随机数序列的独立性越好。2021/3/2766CHENLI

对于乘同余方法，对于乘加同余方法，因此，这两种方法的独立性都是很好的。同伪随机数的均匀性问题一样，伪随机数序列的独立性问题也是对它的全体讨论的。若只考虑伪随机数的一部分，在通常情况下给出ε(i)是相当因难的。因此，伪随机数序列的独立性问题的统计检验方法同样是非常重要的。2021/3/2767CHENLI作业

证明1—ξ是随机数。证明与同分布。2021/3/2768CHENLI第三章由已知分布的随机抽样随机抽样及其特点直接抽样方法挑选抽样方法复合抽样方法复合挑选抽样方法替换抽样方法随机抽样的一般方法随机抽样的其它方法作业2021/3/2769CHENLI第三章由已知分布的随机抽样

本章叙述由己知分布抽样的各主要方法，并给出在粒子输运问题中经常用到的具体实例。2021/3/2770CHENLI随机抽样及其特点

由巳知分布的随机抽样指的是由己知分布的总体中抽取简单子样。随机数序列是由单位均匀分布的总体中抽取的简单子样，属于一种特殊的由已知分布的随机抽样问题。本章所叙述的由任意已知分布中抽取简单子样，是在假设随机数为已知量的前提下，使用严格的数学方法产生的。为方便起见，用XF表示由己知分布F(x)中产生的简单子样的个体。对于连续型分布，常用分布密度函数f(x)表示总体的己知分布，用Xf表示由己知分布密度函数f(x)产生的简单子样的个体。另外，在抽样过程中用到的伪随机数均称随机数。2021/3/2771CHENLI直接抽样方法

对于任意给定的分布函数F(x)，直接抽样方法如下：其中，ξ1，ξ2，…，ξN为随机数序列。为方便起见，将上式简化为：若不加特殊说明，今后将总用这种类似的简化形式表示，ξ总表示随机数。2021/3/2772CHENLI证明

下面证明用前面介绍的方法所确定的随机变量序列X1，X2，…，XN具有相同分布F(x)。对于任意的n成立，因此随机变量序列X1，X2，…，XN具有相同分布F(x)。另外，由于随机数序列ξ1，ξ2，…，ξN是相互独立的，而直接抽样公式所确定的函数是波雷尔（Borel）可测的，因此，由它所确定的X1，X2，…，XN也是相互独立的（[P.R.Halmos,Measuretheory,N.Y.VonNosrtand,1950]§45定理2）。2021/3/2773CHENLI离散型分布的直接抽样方法

对于任意离散型分布：其中x1，x2，…为离散型分布函数的跳跃点，P1，P2，…为相应的概率，根据前述直接抽样法，有离散型分布的直接抽样方法如下：该结果表明，为了实现由任意离散型分布的随机抽样，直接抽样方法是非常理想的。2021/3/2774CHENLI例1.二项分布的抽样

二项分布为离散型分布，其概率函数为：其中，P为概率。对该分布的直接抽样方法如下：2021/3/2775CHENLI例2.泊松(Possion)分布的抽样

泊松(Possion)分布为离散型分布，其概率函数为：其中，λ>0。对该分布的直接抽样方法如下：2021/3/2776CHENLI例3.掷骰子点数的抽样

掷骰子点数X=n的概率为：选取随机数ξ，如则在等概率的情况下，可使用如下更简单的方法：其中［］表示取整数。2021/3/2777CHENLI例4.碰撞核种类的确定

中子或光子在介质中发生碰撞时，如介质是由多种元素组成，需要确定碰撞核的种类。假定介质中每种核的宏观总截面分别为Σ1，Σ2，…，Σn，则中子或光子与每种核碰撞的概率分别为：其中Σt＝Σ1＋Σ2＋…＋Σn。碰撞核种类的确定方法为：产生一个随机数ξ，如果则中子或光子与第I种核发生碰撞。2021/3/2778CHENLI例5.中子与核的反应类型的确定

假设中子与核的反应类型有如下几种:弹性散射，非弹性散射，裂变，吸收，相应的反应截面分别为Σel，Σin，Σf，Σa。则发生每一种反应类型的概率依次为：其中反应总截面Σt＝Σel＋Σin＋Σf＋Σa。2021/3/2779CHENLI

反应类型的确定方法为：产生一个随机数ξ

2021/3/2780CHENLI连续型分布的直接抽样方法

对于连续型分布，如果分布函数F(x)的反函数

F－1(x)存在，则直接抽样方法是：2021/3/2781CHENLI例6.在［a，b］上均匀分布的抽样

在［a，b］上均匀分布的分布函数为：则2021/3/2782CHENLI例7.β分布β分布为连续型分布，作为它的一个特例是：其分布函数为：

则2021/3/2783CHENLI例8.指数分布

指数分布为连续型分布，其一般形式如下：其分布函数为：

则因为1－ξ也是随机数，可将上式简化为2021/3/2784CHENLI

连续性分布函数的直接抽样方法对于分布函数的反函数存在且容易实现的情况，使用起来是很方便的。但是对于以下几种情况，直接抽样法是不合适的。分布函数无法用解析形式给出，因而其反函数也无法给出。分布函数可以给出其解析形式，但是反函数给不出来。分布函数即使能够给出反函数，但运算量很大。下面叙述的挑选抽样方法是克服这些困难的比较好的方法。2021/3/2785CHENLI挑选抽样方法

为了实现从己知分布密度函数f(x)抽样，选取与f(x)取值范围相同的分布密度函数h(x)，如果则挑选抽样方法为：>2021/3/2786CHENLI

即从h(x)中抽样xh，以的概率接受它。下面证明xf

服从分布密度函数f(x)。证明：对于任意x

2021/3/2787CHENLI2021/3/2788CHENLI

使用挑选抽样方法时，要注意以下两点：选取h(x)时要使得h(x)容易抽样且M的值要尽量小。因为M小能提高抽样效率。抽样效率是指在挑选抽样方法中进行挑选时被选中的概率。按此定义，该方法的抽样效率E为：所以，M越小，抽样效率越高。2021/3/2789CHENLI

当f(x)在［0，1］上定义时，取h(x)=1，Xh=ξ，此时挑选抽样方法为>2021/3/2790CHENLI例9.圆内均匀分布抽样

令圆半径为R0，点到圆心的距离为r，则r的分布密度函数为分布函数为容易知道，该分布的直接抽样方法是2021/3/2791CHENLI

由于开方运算在计算机上很费时间，该方法不是好方法。下面使用挑选抽样方法：取则抽样框图为＞≤2021/3/2792CHENLI

显然，没有必要舍弃ξ1＞ξ2的情况，此时，只需取就可以了，亦即另一方面，也可证明与具有相同的分布。2021/3/2793CHENLI复合抽样方法

在实际问题中，经常有这样的随机变量，它服从的分布与一个参数有关，而该参数也是一个服从确定分布的随机变量，称这样的随机变量服从复合分布。例如，分布密度函数是一个复合分布。其中Pn≥0，n=1，2，…，且

fn(x)为与参数n有关的分布密度函数，n=1，2，…，参数n服从如下分布2021/3/2794CHENLI

复合分布的一般形式为：其中f2(x/y)表示与参数y有关的条件分布密度函数，F1(y)表示分布函数。 复合分布的抽样方法为:首先由分布函数F1(y)或分布密度函数f1(y)中抽样YF1或Yf1，然后再由分布密度函数f2(x/YF1)中抽样确定Xf2(x/YF)

证明：所以，Xf所服从的分布为f

(x)。2021/3/2795CHENLI例10.指数函数分布的抽样

指数函数分布的一般形式为：引入如下两个分布密度函数：2021/3/2796CHENLI

则使用复合抽样方法，首先从f1(y)中抽取y

再由f2(x/YF1)中抽取x

2021/3/2797CHENLI复合挑选抽样方法

考虑另一种形式的复合分布如下：其中0≤H(x,y)≤M，f2(x/y)表示与参数y有关的条件分布密度函数，F1(y)表示分布函数。抽样方法如下：>2021/3/2798CHENLI

证明：抽样效率为：E=1/M2021/3/2799CHENLI

为了实现某个复杂的随机变量y的抽样，将其表示成若干个简单的随机变量x1，x2，…，xn

的函数 得到x1，x2，…，xn的抽样后，即可确定y的抽样，这种方法叫作替换法抽样。即替换抽样方法2021/3/27100CHENLI例11.散射方位角余弦分布的抽样

散射方位角φ在［0,2π］上均匀分布，则其正弦和余弦sinφ和cosφ服从如下分布： 直接抽样方法为：2021/3/27101CHENLI

令φ=2θ，则θ在［0,π］上均匀分布，作变换 其中0≤ρ≤1，0≤ρ≤π，则

(x,y)表示上半个单位圆内的点。如果(x,y)在上半个单位圆内均匀分布，则θ在［0,π］上均匀分布，由于2021/3/27102CHENLI

因此抽样sinφ和cosφ的问题就变成在上半个单位圆内均匀抽样(x,y)的问题。 为获得上半个单位圆内 的均匀点，采用挑选法，在 上半个单位圆的外切矩形内 均匀投点（如图）。 舍弃圆外的点，余下的就是所要求的点。 抽样方法为： 抽样效率

E=π/4≈0.785>2021/3/27103CHENLI

为实现散射方位角余弦分布抽样，最重要的是在上半个单位圆内产生均匀分布点。下面这种方法，首先在单位圆的半个外切正六边形内产生均匀分布点，如图所示。2021/3/27104CHENLI

于是便有了抽样效率更高的抽样方法： 抽样效率>≤2021/3/27105CHENLI例12.正态分布的抽样

标准正态分布密度函数为： 引入一个与标准正态随机变量X独立同分布的随机变量Y，则（X,Y）的联合分布密度为： 作变换2021/3/27106CHENLI

则（ρ,φ）的联合分布密度函数为： 由此可知，ρ与φ相互独立，其分布密度函数分别为 分别抽取ρ，φ

：2021/3/27107CHENLI

从而得到一对服从标准正态分布的随机变量X和Y： 对于一般的正态分布密度函数N(μ,σ2)的抽样，其抽样结果为：2021/3/27108CHENLI例13.β分布的抽样

β分布密度函数的一般形式为： 其中n，k为整数。为了实现β分布的抽样，将其看作一组简单的相互独立随机变量的函数，通过这些简单随机变量的抽样，实现β分布的抽样。设x1，x2，…，xn

为一组相互独立、具有相同分布F(x)的随机变量，ζk为x1，x2，…，xn

按大小顺序排列后的第k个，记为：2021/3/27109CHENLI

则ζk的分布函数为： 当F(x)=x

时， 不难验证，ζk的分布密度函数为β分布。因此，β分布的抽样可用如下方法实现： 选取n个随机数，按大小顺序排列后取第k个，即2021/3/27110CHENLI随机抽样的一般方法

加抽样方法

减抽样方法乘抽样方法乘加抽样方法乘减抽样方法对称抽样方法积分抽样方法2021/3/27111CHENLI加抽样方法

加抽样方法是对如下加分布给出的一种抽样方法：其中Pn≥0，

，且

fn(x)为与参数n有关的分布密度函数，n=1，2，…。 由复合分布抽样方法可知，加分布的抽样方法为:首先抽样确定n’，然后由fn’(x)中抽样x，即：2021/3/27112CHENLI例14.多项式分布抽样

多项式分布密度函数的一般形式为：将f(x)改写成如下形式：则该分布的抽样方法为：2021/3/27113CHENLI例15.球壳内均匀分布抽样

设球壳内半径为R0，外半径为R1，点到球心的距离为r，则r的分布密度函数为分布函数为 该分布的直接抽样方法是2021/3/27114CHENLI

为避免开立方根运算，作变换： 则x∈[0,1]，其分布密度函数为： 其中2021/3/27115CHENLI

则x及r的抽样方法为：≤≤>>2021/3/27116CHENLI减抽样方法

减抽样方法是对如下形式的分布密度所给出的一种抽样方法：其中A1、A2为非负实数，f1(x)

、f2(x)均为分布密度函数。 减抽样方法分为以下两种形式：

以上两种形式的抽样方法，究竟选择哪种好，要看f1(x)

、f2(x)哪一个容易抽样，如相差不多，选用第一种方法抽样效率高。2021/3/27117CHENLI

（1）将f

(x)表示为令m表示f2(x)／f1(x)的下界，使用挑选法，从f1(x)中抽取Xf1

抽样效率为：>2021/3/27118CHENLI

（2）将f

(x)表示为使用挑选法，从f2(x)中抽取Xf2

抽样效率为：>2021/3/27119CHENLI例16.β分布抽样 β分布的一个特例： 取A1＝2，A2＝1，f1(x)＝1，f2(x)＝2x，此时m＝0，则根据第一种形式的减抽样方法，有 或＞≤＞≤2021/3/27120CHENLI

由于1－ξ1可用ξ1代替，该抽样方法可简化为： 对于ξ2＞ξ1的情况，可取Xf＝ξ1

，因此 与β分布的推论相同。＞≤2021/3/27121CHENLI

如下形式的分布称为乘分布：其中H(x)为非负函数，

f1(x)为任意分布密度函数。 令M为H(x)的上界，乘抽样方法如下：抽样效率为：乘抽样方法≤＞2021/3/27122CHENLI例17.倒数分布抽样

倒数分布密度函数为： 其直接抽样方法为： 下面采用乘抽样方法，考虑如下分布族： 其中i=1，2，…，该分布的直接抽样方法为：2021/3/27123CHENLI

利用这一分布族，将倒数分布f(x)表示成:

其中， 乘法分布的抽样方法如下:

该分布的抽样效率为：＞≤2021/3/27124CHENLI例18.麦克斯韦(Maxwell)分布抽样

麦克斯韦分布密度函数的一般形式为： 使用乘抽样方法，令 该分布的直接抽样方法为：2021/3/27125CHENLI

此时 则麦克斯韦分布的抽样方法为:

该分布的抽样效率为：＞≤2021/3/27126CHENLI

在实际问题中，经常会遇到如下形式的分布：其中Hn(x)为非负函数，fn(x)为任意分布密度函数，n=1，2，…。不失一般性，只考虑n=2的情况：

将f(x)改写成如下的加分布形式：乘加抽样方法2021/3/27127CHENLI

其中2021/3/27128CHENLI

乘加抽样方法为：该方法的抽样效率为：>>>≤2021/3/27129CHENLI

这种方法需要知道P1的值(P2=1－P1），这对有些分布是很困难的。下面的方法可以不用计算P1

：对于任意小于1的正数P1

，令P2=1－P1

；

则采用复合挑选抽样方法，有：2021/3/27130CHENLI

当取时，抽样效率最高这时，乘加抽样方法为：>>>≤2021/3/27131CHENLI

由于可知第一种方法比第二种方法的抽样效率高。2021/3/27132CHENLI例19.光子散射后能量分布的抽样

令光子散射前后的能量分别为

和（以m0c2

为单位，m0为电子静止质量，c

为光速），， 则x

的分布密度函数为： 该分布即为光子散射能量分布，它是由著名的Klin－Nishina

公式确定的。其中K(α)为归一因子：2021/3/27133CHENLI

把光子散射能量分布改写成如下形式： 在［1,1+2α］上定义如下函数:2021/3/27134CHENLI

则有 使用乘加抽样方法：2021/3/27135CHENLI

光子散射能量分布的抽样方法为： 该方法的抽样效率为：＞＞＞≤≤≤2021/3/27136CHENLI

乘减分布的形式为： 其中H1(x)、H2(x)为非负函数，f1(x)、f2(x)为任意分布密度函数。 与减抽样方法类似，乘减分布的抽样方法也分为两种。乘减抽样方法2021/3/27137CHENLI

（1）将f

(x)表示为 令H1(x)的上界为M1，的下界为m，使用乘抽 样方法得到如下乘减抽样方法：>2021/3/27138CHENLI

（2）将f

(x)表示为 令H2(x)的上界为M2，使用乘抽样方法，得到另一种乘减抽样方法：>2021/3/27139CHENLI例20.裂变中子谱分布抽样

裂变中子谱分布的一般形式为： 其中A，B，C，Emin，Emax

均为与元素有关的量。令 其中λ为归一因子，γ为任意参数。2021/3/27140CHENLI

相应的H1(E)，H2(E)为： 于是裂变中子谱分布可以表示成乘减分布形式:

容易确定H1(E)的上界为： 为提高抽样效率，应取γ使得M1

达到最小，此时2021/3/27141CHENLI

取m＝0，令 则裂变中子谱分布的抽样方法为:

抽样效率＞≤2021/3/27142CHENLI

对称分布的一般形式为： 其中f1(x)为任意分布密度函数，满足偶函数对称条件，H(x)为任意奇函数，即对任意x满足： 对称分布的抽样方法如下：取η=2ξ－1对称抽样方法>≤2021/3/27143CHENLI

证明： 因为η=2ξ－1，η≤x

相当于ξ≤，因此2021/3/27144CHENLI例21.质心系各向同性散射角余弦分布抽样

在质心系各向同性散射的假设下，为得到实验室系散射角余弦，需首先抽样确定质心条散射角余弦： 再利用下面转换公式:

得到实验室系散射角余弦μL。其中A为碰撞核质量，θC、θL

分别为质心系和实验室系散射角。2021/3/27145CHENLI

为避免开方运算，可以使用对称分布抽样。 根据转换公式可得： 依照质心系散射各向同性的假定，可得到实验室系散射角余弦μL

的分布如下： 该密度函数中的第一项为偶函数，第二项为奇函数，因而是对称分布。其中2021/3/27146CHENLI

从f1(μL)的抽样可使用挑选法 然后再以 的概率决定接受或取负值。 上述公式涉及开方运算，需要进一步简化。＞≤2021/3/27147CHENLI

注意以下事实：对于任意0≤a≤1

令 则上述挑选抽样中的挑选条件简化为： 另一方面，在即的条件下，η2/a

在［－1,1］上均匀分布，故可令η＝η2/a，则最终决定取正负值的条件简化为：2021/3/27148CHENLI

于是，得到质心系各向同性散射角余弦分布的抽样方法为：＞≤＞≤2021/3/27149CHENLI

如下形式的分布密度函数 称为积分分布密度函数，其中f0(x，y)为任意二维分布密度函数，H(x)为任意函数。该分布密度函数的抽样方法为：积分抽样方法>2021/3/27150CHENLI

证明：对于任意x

2021/3/27151CHENLI例22.各向同性散射方向的抽样

为了确定各向同性散射方向，根据公式： 对于各向同性散射，cosθ在［－1,1］上均匀分布，φ在［0,2π］上均匀分布。由于 直接抽样需要计算三角函数和开方。2021/3/27152CHENLI

定义两个随机变量： 可以证明，当时，随机变量x

和y

服从如下分布:

定义区域为：2021/3/27153CHENLI

则w＝cosθ的分布可以用上述分布表示成积分分布的形式： 令，则属于上述积分限内的y

一定满足 条件。2021/3/27154CHENLI

各向同性散射方向的抽样方法为： 抽样效率为：＞≤2021/3/27155CHENLI随机抽样的其它方法

偏倚抽样方法近似抽样方法近似-修正抽样方法多维分布抽样方法指数分布的抽样2021/3/27156CHENLI

使用蒙特卡罗方法计算积分 时，可考虑将积分I改写为 其中f*(x)为一个与f(x)有相同定义域的新的分布密度函数。于是可以这样计算积分I： 这里Xi

是从f*(x)中抽取的第i

个子样。偏移抽样方法2021/3/27157CHENLI

由此可以看出，原来由f(x)抽样，现改为由另一个分布密度函数f*(x)抽样，并附带一个权重纠偏因子 这种方法称为偏倚抽样方法。 从f(x)中抽取的Xf

，满足 而对于偏倚抽样，有 一般情况下，Xf

是具有分布f(x)总体的简单子样的个体，只代表一个。Xf*

是具有分布f*(x)总体的简单子样的个体，但不代表一个，而是代表W(Xf*)个，这时Xf*是带权W(Xf*)服从分布f(x)。2021/3/27158CHENLI

在实际问题中，分布密度函数的形式有时是非常复杂的，有些甚至不能用解析形式给出，只能用数据或曲线形式给出。如中子散射角余弦分布多数是以曲线形式给出的。对于这样的分布，需要用近似分布密度函数代替原来的分布密度函数，用近似分布密度函数的抽样代替原分布密度函数的抽样，这种方法称为近似抽样方法。近似抽样方法2021/3/27159CHENLI

设fa(x)≈f(x)，即fa(x)是f(x)的一个近似分布密度函数。对于阶梯近似，有 其中，x0，x1，…，xn为任意分点。在此情况下，近似抽样方法为：或阶梯近似2021/3/27160CHENLI

对于梯形近似，有 其中，c

为归一因子，fi

＝f(xi)，x0，x1，…，xn为任意分点。根据对称抽样方法，梯形近似抽样方法为：梯形近似＞≤2021/3/27161CHENLI

除了上述这种近似外，近似抽样方法还包括对直接抽样方法中分布函数反函数的近似处理，以及用具有近似分布的随机变量代替原分布的随机变量。2021/3/27162CHENLI例23.正态分布的近似抽样

我们知道，随机数ξ的期望值为1/2，方差为1/12，则随机变量 渐近正态分布，因此，当n

足够大时便可用Xn

作为正态分布的近似抽样。特别是n＝12时，有2021/3/27163CHENLI

对于任意分布密度函数f(x)，设fa(x)是f(x)的一个近似分布密度函数，它的特点是抽样简单，运算量小。令 则分布密度函数f(x)可以表示为乘加分布形式： 其中H1(x)为非负函数，f1(x)为一分布密度函数。 对f(x)而言，fa(x)是它的近似分布密度函数，而H1(x)f1(x)正好是这种近似的修正。近似-修正抽样方法2021/3/27164CHENLI

近似-修正抽样方法如下： 抽样效率 由上述近似-修正抽样方法可以看出，如果近似分布密度函数fa(x)选得好，m

接近1，这时有很大可能直接从fa(x)中抽取Xfa

，而只有很少的情况需要计算与f

(x)有关的函数H1(Xf1)。在乘抽样方法中，每一次都要计算H(Xfa)＝f

(Xfa)／fa(Xfa)。因此，当f

(x)比较复杂时，近似-修正抽样方法有很大好处。≤≤＞＞2021/3/27165CHENLI例24.裂变中子谱分布的近似-修正抽样

裂变中子谱分布的一般形式为： 其中A，B，C，Emin，Emax

均为与元素有关的量。 对于铀-235，

A=0.965，B=2.29，C=0.453，Emin=0，Emax=∞。 若采用乘减抽样方法，其抽样效率约为0.5。2021/3/27166CHENLI

令 相应的 则 从fa(x)的抽样为 从f1(x)的抽样为2021/3/27167CHENLI

参数λ的确定，使1－Aλ＞0，且使H1(E)的上界M1

最小。裂变中子谱的近似修正抽样方法为 对于铀-235，m≈0.8746，M≈0.2678，λ≈0.5543，抽样效率E≈0.9333。而且近似修正抽样方法有0.8746的概率直接用近似分布抽样，只计算一次对数。因此，较之乘减抽样方法大大节省了计算时间，提高了抽样效率。≤≤＞＞2021/3/27168CHENLI

为方便起见，这里仅讨论二维分布的情况，对于更高维数的分布，可用类似的方法处理。 对于任意二维分布密度函数，总可以用其边缘分布密度函数和条件分布密度函数的乘积表示:

其中fl(x)，f2(y|x)分别为分布f(x,y)的边缘分布密度函数和条件分布密度函数，即多维分布抽样方法2021/3/27169CH