蒙特卡洛在核技术中的应用讲义完整版ppt课件.ppt

《蒙特卡洛在核技术中的应用讲义完整版ppt课件.ppt》由会员分享，可在线阅读，更多相关《蒙特卡洛在核技术中的应用讲义完整版ppt课件.ppt（402页珍藏版）》请在三一办公上搜索。

1、蒙特卡罗方法在核技术中的应用,林谦,目 录,第一章 蒙特卡罗方法概述第二章 随机数第三章 由已知分布的随机抽样第四章 蒙特卡罗方法解粒子输运问题,教材,蒙特卡罗方法在实验核物理中的应用许淑艳 编著原子能出版社蒙特卡罗方法清华大学,参考书,蒙特卡罗方法及其在粒子输运问题中的应用裴鹿成 张孝泽 编著科学出版社蒙特卡罗方法徐钟济 编著上海科学技术出版社,联系方式,电话83918电子邮件,第一章 蒙特卡罗方法概述,蒙特卡罗方法的基本思想蒙特卡罗方法的收敛性，误差蒙特卡罗方法的特点蒙特卡罗方法的主要应用范围作 业,第一章 蒙特卡罗方法概述,蒙特卡罗方法又称随机抽样技巧或统计试验方法。半个多世纪以来，由于

2、科学技术的发展和电子计算机的发明 ，这种方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。,蒙特卡罗方法的基本思想,二十世纪四十年代中期，由于科学技术的发展和电子计算机的发明，蒙特卡罗方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。但其基本思想并非新颖，人们在生产实践和科学试验中就已发现，并加以利用。两个例子 例1. 蒲丰氏问题 例2

3、. 射击问题（打靶游戏）基本思想计算机模拟试验过程,例1. 蒲丰氏问题,为了求得圆周率值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（ la）的平行线相交的频率代替概率P，再利用准确的关系式： 求出值 其中为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。,一些人进行了实验，其结果列于下表 ：,例2. 射击问题（打靶游戏）,设r表示射击运动员的弹着点到靶心的距离，(r)表示击中r处相应的得分数（环数），f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射击水平。该运动员的射击成绩为 用概率语言来说，是随机变量(

4、r)的数学期望，即,现假设该运动员进行了次射击，每次射击的弹着点依次为r1，r2，rN，则次得分g(r1)，g(r2)，g(rN)的算术平均值 代表了该运动员的成绩。换言之，为积分的估计值，或近似值。 在该例中，用次试验所得成绩的算术平均值作为数学期望的估计值（积分近似值）。,基本思想,由以上两个例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。 当随机变量的取值仅为1或0时，它的数学期望就是某个事件

5、的概率。或者说，某种事件的概率也是随机变量（仅取值为1或0）的数学期望。,因此，可以通俗地说，蒙特卡罗方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数f(r)的随机变量(r)的数学期望 通过某种试验，得到个观察值r1，r2，rN（用概率语言来说，从分布密度函数f(r)中抽取个子样r1，r2，rN，），将相应的个随机变量的值g(r1)，g(r2)，g(rN)的算术平均值 作为积分的估计值（近似值）。,为了得到具有一定精确度的近似解，所需试验的次数是很多的，通过人工方法作大量的试验相当困难，甚至是不可能的。因此，蒙特卡罗方法的基本思想虽然早已被人们提出，却很少被使用。本世

6、纪四十年代以来，由于电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验交由计算机完成，使得蒙特卡罗方法得以广泛地应用，在现代化的科学技术中发挥应有的作用。,计算机模拟试验过程,计算机模拟试验过程，就是将试验过程（如投针，射击）化为数学问题，在计算机上实现。以上述两个问题为例，分别加以说明。 例1. 蒲丰氏问题 例2. 射击问题（打靶游戏） 由上面两个例题看出，蒙特卡罗方法常以一个“概率模型”为基础，按照它所描述的过程，使用由已知分布抽样的方法，得到部分试验结果的观察值，求得问题的近似解。,例蒲丰氏问题,设针投到地面上的位置可以用一组参数（x,）来描述，x为针中

7、心的坐标，为针与平行线的夹角，如图所示。 任意投针，就是意味着x与都是任意取的，但x的范围限于0，a，夹角的范围限于0，。在此情况下，针与平行线相交的数学条件是,针在平行线间的位置,如何产生任意的（x,）？x在0，a上任意取值，表示x在0，a上是均匀分布的，其分布密度函数为： 类似地，的分布密度函数为： 因此，产生任意的（x,）的过程就变成了由f1(x)抽样x及由f2()抽样的过程了。由此得到： 其中1，2均为（0,1）上均匀分布的随机变量。,每次投针试验，实际上变成在计算机上从两个均匀分布的随机变量中抽样得到（x,），然后定义描述针与平行线相交状况的随机变量s(x,)，为 如果投针次，则 是

8、针与平行线相交概率的估计值。事实上， 于是有,例射击问题,设射击运动员的弹着点分布为 用计算机作随机试验（射击）的方法为，选取一个随机数，按右边所列方法判断得到成绩。 这样，就进行了一次随机试验（射击），得到了一次成绩 (r)，作次试验后，得到该运动员射击成绩的近似值,蒙特卡罗方法的收敛性，误差,蒙特卡罗方法作为一种计算方法，其收敛性与误差是普遍关心的一个重要问题。收敛性误差减小方差的各种技巧 效率,收敛性,由前面介绍可知，蒙特卡罗方法是由随机变量X的简单子样X1，X2，XN的算术平均值： 作为所求解的近似值。由大数定律可知， 如X1，X2，XN独立同分布，且具有有限期望值（E(X)），则 即

9、随机变量X的简单子样的算术平均值 ，当子样数充分大时，以概率1收敛于它的期望值E(X)。,误差,蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，XN独立同分布，且具有有限非零的方差2 ，即 f(X)是X的分布密度函数。则,当N充分大时，有如下的近似式 其中称为置信度，1称为置信水平。 这表明，不等式 近似地以概率 1成立，且误差收敛速度的阶为 。 通常，蒙特卡罗方法的误差定义为 上式中 与置信度是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出 。,下面给出几个常用的与的数值： 关于蒙特卡罗方法的误差需说明

10、两点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。第二，误差中的均方差是未知的，必须使用其估计值 来代替，在计算所求量的同时，可计算出 。,减小方差的各种技巧,显然，当给定置信度后，误差由和N决定。要减小，或者是增大N，或者是减小方差2。在固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。 另一方面，如能减小估计的均方差，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果。因此降低方差的各种技巧，引起了人们的普遍注意。后面课程将会介绍一些降低方差的技巧。,效率,一般来说，降低方差的技巧，往往会使观察一个子样的时间

11、增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就 是蒙特卡罗方法中效率的概念。它定义为 ，其中c 是观察一个子样的平均费用。显然 越小，方法越有效。,蒙特卡罗方法的特点,优点能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。受几何条件限制小。收敛速度与问题的维数无关。具有同时计算多个方案与多个未知量的能力。误差容易确定。程序结构简单，易于实现。,缺点收敛速度慢。误差具有概率性。在粒子输运问题中，计算结果与系统大小有关。,能够比较逼真地描述具有随机性质的事物的特点及物理实验过程,从这个意义上讲，蒙特卡罗方法可

12、以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。,受几何条件限制小,在计算s维空间中的任一区域Ds上的积分 时，无论区域Ds的形状多么特殊，只要能给出描述Ds的几何特征的条件，就可以从Ds中均匀产生N个点 ，得到积分的近似值。 其中Ds为区域Ds的体积。这是数值方法难以作到的。 另外，在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。,收敛速度与问题的维数无关,由误差定义可知，在给定置信水平情况下，蒙特卡罗方法的收敛速度

13、为，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计算时间的变化，不影响误差。也就是说，使用蒙特卡罗方法时，抽取的子样总数N与维数s无关。维数的增加，除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。,具有同时计算多个方案与多个未知量的能力,对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计

14、算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。 另外，使用蒙特卡罗方法还可以同时得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。,误差容易确定,对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。 一般计算方法常存在着有效位数损失问题，而要解决

15、这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。,程序结构简单，易于实现,在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。,收敛速度慢,如前所述，蒙特卡罗方法的收敛速度为 ，一般不容易得到精确度较高的近似结果。对于维数少（三维以下）的问题，不如其他方法好。,误差具有概率性,由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。,在粒子输运问题中，计算结果与系统大小有关,经验表明，只有当系统的大小与粒子的平均自由程可以相比较时（一般在十个平均自由程左右），蒙特卡罗方法计算的结果较为满意。但对于大系统或小概率事件的计算问题，计算结果

16、往往比真值偏低。而对于大系统，数值方法则是适用的。 因此，在使用蒙特卡罗方法时，可以考虑把蒙特卡罗方法与解析（或数值）方法相结合，取长补短，既能解决解析（或数值）方法难以解决的问题，也可以解决单纯使用蒙特卡罗方法难以解决的问题。这样，可以发挥蒙特卡罗方法的特长，使其应用范围更加广泛。,蒙特卡罗方法的主要应用范围,蒙特卡罗方法所特有的优点，使得它的应用范围越来越广。它的主要应用范围包括：粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。随着科学技术的发展，其应用范围将更加广泛。 蒙特卡罗方法在粒子输运问题中的应用范围主要包括：实验核物理，反应堆物理，高能物理等方

17、面。 蒙特卡罗方法在实验核物理中的应用范围主要包括：通量及反应率，中子探测效率，光子探测效率，光子能量沉积谱及响应函数，气体正比计数管反冲质子谱，多次散射与通量衰减修正等方面。,作 业,用蒲丰投针法在计算机上计算值，取a=4、l=3。分别用理论计算和计算机模拟计算，求连续掷两颗骰子，点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率。,第二章 随机数,随机数的定义及产生方法 伪随机数产生伪随机数的乘同余方法产生伪随机数的乘加同余方法产生伪随机数的其他方法伪随机数序列的均匀性和独立性作 业,第二章 随机数,由具有已知分布的总体中抽取简单子样，在蒙特卡罗方法中占有非常重要的地位。总体和子样的

18、关系，属于一般和个别的关系，或者说属于共性和个性的关系。由具有已知分布的总体中产生简单子样，就是由简单子样中若干个性近似地反映总体的共性。 随机数是实现由已知分布抽样的基本量，在由已知分布的抽样过程中，将随机数作为已知量，用适当的数学方法可以由它产生具有任意已知分布的简单子样。,随机数的定义及产生方法,随机数的定义及性质随机数表物理方法,随机数的定义及性质,在连续型随机变量的分布中，最简单而且最基本的分布是单位均匀分布。由该分布抽取的简单子样称，随机数序列，其中每一个体称为随机数。 单位均匀分布也称为0，1上的均匀分布，其分布密度函数为： 分布函数为 ：,由于随机数在蒙特卡罗方法中占有极其重要

19、的位置，我们用专门的符号表示。由随机数序列的定义可知，1，2，是相互独立且具有相同单位均匀分布的随机数序列。也就是说，独立性、均匀性是随机数必备的两个特点。 随机数具有非常重要的性质：对于任意自然数s，由s个随机数组成的s维空间上的点（n+1，n+2，n+s）在s维空间的单位立方体Gs上均匀分布，即对任意的ai， 如下等式成立：,其中P()表示事件发生的概率。反之，如果随机变量序列1, 2对于任意自然数s，由s个元素所组成的s维空间上的点（n+1，n+s）在Gs上均匀分布，则它们是随机数序列。 由于随机数在蒙特卡罗方法中所处的特殊地位，它们虽然也属于由具有已知分布的总体中产生简单子样的问题，但

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

21、用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。用来作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。 一般情况下，任意一个随机数在计算机内总是用二进制的数表示的： 其中i（i=1，2，m）或者为0，或者为1。,因此，利用物理方法在计算机上产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为0.5。 用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，需要增加随机数发生器和电路联系等附加设备，费用昂贵。因此，该方法也不适

22、合在计算机上使用。,伪随机数,伪随机数伪随机数存在的两个问题伪随机数的周期和最大容量,伪随机数,在计算机上产生随机数最实用、最常见的方法是数学方法，即用如下递推公式： 产生随机数序列。对于给定的初始值1,2，k，确定n+k，=1，2，。经常使用的是k=1的情况，其递推公式为： 对于给定的初始值1，确定n+1，=,伪随机数存在的两个问题,用数学方法产生的随机数，存在两个问题：递推公式和初始值1,2，k确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。由于随机数序列是由递推公式确定的，而在计算机上所能表示的0，1上的数又是有限的，因此，这种方法产生的随机数序列就不可能不出现无限重复。

23、一旦出现这样的n，n (n n )，使得下面等式成立： 随机数序列便出现了周期性的循环现象。对于k=1的情况，只要有一个随机数重复，其后面的随机数全部重复，这与随机数的要求是不相符的。,由于这两个问题的存在，常称用数学方法产生的随机数为伪随机数。对于以上存在的两个问题，作如下具体分析。 关于第一个问题，不能从本质上加以改变，但只要递推公式选得比较好，随机数间的相互独立性是可以近似满足的。至于第二个问题，则不是本质的。因为用蒙特卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所用随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。 用数学方法产生的伪随机数容易在计算机上得到

24、，可以进行复算，而且不受计算机型号的限制。因此，这种方法虽然存在着一些问题，但仍然被广泛地在计算机上使用，是在计算机上产生伪随机数的主要方法。,伪随机数的周期和最大容量,发生周期性循环现象的伪随机数的个数称为伪随机数的周期。对于前面介绍的情况，伪随机数的周期为nn。 从伪随机数序列的初始值开始，到出现循环现象为止，所产生的伪随机数的个数称为伪随机数的最大容量。前面的例子中，伪随机数的最大容量为n 。,产生伪随机数的乘同余方法,乘同余方法是由Lehmer在1951年提出来的，它的一般形式是：对于任一初始值x1，伪随机数序列由下面递推公式确定： 其中a为常数。,乘同余方法的最大容量的上界,对于任意

25、正整数M，根据数论中的标准分解定理，总可以分解成如下形式： 其中P0=2，P1， Pr表示不同的奇素数，0表示非负整数，1，r表示正整数。a无论取什么值，乘同余方法的最大容量的上界为： 的最小公倍数。其中：,关于a与x1的取值,如果a与x1满足如下条件： 对于 ， x1与M互素，则乘同余方法产生的伪随机数序列的最大容量达到最大可能值(M)。,乘同余方法在计算机上的使用,为了便于在计算机上使用，通常取 ：=2s 其中s为计算机中二进制数的最大可能有效位数x1= 奇数 a = 52k+1 其中k为使52k+1在计算机上所能容纳的最大整数，即a为计算机上所能容纳的5的最大奇次幂。一般地，s=32时，

26、a=513；s=48，a=515等。伪随机数序列的最大容量(M)=2s-2 。 乘同余方法是使用的最多、最广的方法，在计算机上被广泛地使用。,产生伪随机数的乘加同余方法,产生伪随机数的乘加同余方法是由Rotenberg于1960年提出来的，由于这个方法有很多优点，已成为仅次于乘同余方法产生伪随机数的另一主要方法。 乘加同余方法的一般形式是，对任意初始值x1，伪随机数序列由下面递推公式确定： 其中a和c为常数。,乘加同余方法的最大容量,关于乘加同余方法的最大容量问题，有如下结论：如果对于正整数M的所有素数因子P，下式均成立： 当M为4的倍数时，还有下式成立： c与M互素，则乘加同余方法所产生的伪

27、随机数序列的最大容量达到最大可能值M。,M，x1，a，c的取值,为了便于在计算机上使用，通常取M = 2s 其中s为计算机中二进制数的最大可能有效位数。a = 2b + 1(b2)c = 1 这样在计算中可以使用移位和指令加法，提高计算速度。,产生伪随机数的其他方法,取中方法加同余方法,伪随机数序列的均匀性和独立性,判断伪随机数序列是否满足均匀和相互独立的要求，要靠统计检验的方法实现。对于伪随机数的统计检验，一般包括两大类：均匀性检验和独立性检验。 六十年代初，人们开始用定性的方法研究伪随机数序列的均匀性和独立性问题，简要叙述如下。,伪随机数的均匀性,这里只考虑伪随机数序列1,2，n全体作为子

28、样时的均匀性问题。其中n为伪随机数序列的最大容量。 对于任意的0 x1，令Nn(x)表示伪随机数序列1,2，n中适合不等式i x i=1，2，n 的个数，则 标志伪随机数序列1,2，n的均匀程度，称为均匀偏度。,将伪随机数序列1,2，n从小至大重新排列 并令 ，则由(n)的定义，容易证明 很明显，对于固定的，(n)的值越小越好。它是描述伪随机数序列均匀程度的基本量。对于任意随机数序列，均有如下不等式成立： 当 时，所对应的伪随机数序列为最佳分布。,可以证明，伪随机数序列为最佳分布的充要条件是它取遍序列 的所有值。 对于计算机上使用的乘同余方法，按照前面介绍的方法选取a、x1时，所产生的伪随机数

29、序列的均匀偏度 对于乘加同余方法 对于部分伪随机数的均匀性问题通常用统计检验方法检验。,伪随机数的独立性,对于任意 ，令 表示(1,2)， (2,3)， (n,n+1)中适合不等式 的个数，根据随机变量间相互独立的定义和频率近似概率的方法，令 则(n)标志伪随机数序列1,2，n的独立程度，简称为独立偏度。对于固定的n,(n)的值越接近于零，伪随机数序列的独立性越好。,对于乘同余方法， 对于乘加同余方法， 因此，这两种方法的独立性都是很好的。 同伪随机数的均匀性问题一样，伪随机数序列的独立性问题也是对它的全体讨论的。若只考虑伪随机数的一部分，在通常情况下给出(i)是相当因难的。因此，伪随机数序列

30、的独立性问题的统计检验方法同样是非常重要的。,作 业,证明1是随机数。证明 与 同分布 。,第三章 由已知分布的随机抽样,随机抽样及其特点直接抽样方法挑选抽样方法复合抽样方法复合挑选抽样方法替换抽样方法随机抽样的一般方法随机抽样的其它方法作 业,第三章 由已知分布的随机抽样,本章叙述由己知分布抽样的各主要方法，并给出在粒子输运问题中经常用到的具体实例。,随机抽样及其特点,由巳知分布的随机抽样指的是由己知分布的总体中抽取简单子样。随机数序列是由单位均匀分布的总体中抽取的简单子样，属于一种特殊的由已知分布的随机抽样问题。本章所叙述的由任意已知分布中抽取简单子样，是在假设随机数为已知量的前提下，使用

31、严格的数学方法产生的。 为方便起见，用XF表示由己知分布F(x)中产生的简单子样的个体。对于连续型分布，常用分布密度函数f(x)表示总体的己知分布，用Xf表示由己知分布密度函数f(x)产生的简单子样的个体。另外，在抽样过程中用到的伪随机数均称随机数。,直接抽样方法,对于任意给定的分布函数F(x)，直接抽样方法如下： 其中，1，2，N为随机数序列。为方便起见，将上式简化为： 若不加特殊说明，今后将总用这种类似的简化形式表示，总表示随机数。,证明,下面证明用前面介绍的方法所确定的随机变量序列X1，X2，XN具有相同分布F(x)。 对于任意的n成立，因此随机变量序列X1，X2，XN具有相同分布F(x

32、)。另外，由于随机数序列1，2，N是相互独立的，而直接抽样公式所确定的函数是波雷尔（Borel）可测的，因此，由它所确定的X1，X2，XN也是相互独立的（P.R.Halmos, Measure theory, N.Y.Von Nosrtand,195045定理2）。,离散型分布的直接抽样方法,对于任意离散型分布： 其中x1，x2，为离散型分布函数的跳跃点，P1，P2，为相应的概率，根据前述直接抽样法，有离散型分布的直接抽样方法如下： 该结果表明，为了实现由任意离散型分布的随机抽样，直接抽样方法是非常理想的。,例1. 二项分布的抽样,二项分布为离散型分布，其概率函数为： 其中，P为概率。对该分布

33、的直接抽样方法如下：,例2. 泊松(Possion)分布的抽样,泊松(Possion)分布为离散型分布，其概率函数为： 其中，0 。对该分布的直接抽样方法如下：,例3. 掷骰子点数的抽样,掷骰子点数X=n的概率为： 选取随机数，如 则 在等概率的情况下，可使用如下更简单的方法： 其中表示取整数。,例4. 碰撞核种类的确定,中子或光子在介质中发生碰撞时，如介质是由多种元素组成，需要确定碰撞核的种类。假定介质中每种核的宏观总截面分别为1，2，n，则中子或光子与每种核碰撞的概率分别为： 其中t12n。碰撞核种类的确定方法为：产生一个随机数，如果 则中子或光子与第I种核发生碰撞。,例5. 中子与核的反

34、应类型的确定,假设中子与核的反应类型有如下几种:弹性散射，非弹性散射，裂变，吸收，相应的反应截面分别为el，in，f，a。则发生每一种反应类型的概率依次为 ： 其中反应总截面telinfa。,反应类型的确定方法为：产生一个随机数,连续型分布的直接抽样方法,对于连续型分布，如果分布函数F(x) 的反函数 F1(x)存在，则直接抽样方法是 ：,例6. 在a，b上均匀分布的抽样,在a，b上均匀分布的分布函数为： 则,例7. 分布,分布为连续型分布，作为它的一个特例是： 其分布函数为： 则,例8. 指数分布,指数分布为连续型分布，其一般形式如下： 其分布函数为： 则 因为1也是随机数，可将上式简化为,

35、连续性分布函数的直接抽样方法对于分布函数的反函数存在且容易实现的情况，使用起来是很方便的。但是对于以下几种情况，直接抽样法是不合适的。分布函数无法用解析形式给出，因而其反函数也无法给出。分布函数可以给出其解析形式，但是反函数给不出来。分布函数即使能够给出反函数，但运算量很大。 下面叙述的挑选抽样方法是克服这些困难的比较好的方法。,挑选抽样方法,为了实现从己知分布密度函数f(x)抽样，选取与f(x)取值范围相同的分布密度函数h(x)，如果 则挑选抽样方法为：,即从h(x)中抽样xh，以 的概率接受它。 下面证明xf 服从分布密度函数f(x)。证明：对于任意x,使用挑选抽样方法时，要注意以下两点：

36、选取h(x)时要使得h(x)容易抽样且M的值要尽量小。因为M小能提高抽样效率。抽样效率是指在挑选抽样方法中进行挑选时被选中的概率。按此定义，该方法的抽样效率E为： 所以，M越小，抽样效率越高。,当 f(x) 在0，1上定义时，取 h(x)=1，Xh=， 此时挑选抽样方法为,例9. 圆内均匀分布抽样,令圆半径为R0，点到圆心的距离为r，则r的分布密度函数为 分布函数为 容易知道，该分布的直接抽样方法是,由于开方运算在计算机上很费时间，该方法不是好方法。下面使用挑选抽样方法：取 则抽样框图为,显然，没有必要舍弃12的情况，此时，只需取 就可以了，亦即 另一方面，也可证明 与 具有相同的分布 。,复

37、合抽样方法,在实际问题中，经常有这样的随机变量，它服从的分布与一个参数有关，而该参数也是一个服从确定分布的随机变量，称这样的随机变量服从复合分布。例如，分布密度函数 是一个复合分布。其中Pn0，n=1，2，且 fn(x)为与参数n有关的分布密度函数，n=1，2， 参数n服从如下分布,复合分布的一般形式为： 其中f2(x/y)表示与参数y有关的条件分布密度函数， F1(y)表示分布函数。 复合分布的抽样方法为:首先由分布函数F1(y) 或分布密度函数f1(y)中抽样YF1或Yf1，然后再由分布密度函数f2(x/ YF1)中抽样确定Xf2 (x/YF) 证明： 所以，Xf所服从的分布为f (x)。

38、,例10. 指数函数分布的抽样,指数函数分布的一般形式为： 引入如下两个分布密度函数：,则 使用复合抽样方法，首先从f1(y)中抽取y 再由f2(x/ YF1)中抽取x,复合挑选抽样方法,考虑另一种形式的复合分布如下： 其中0H(x,y)M，f2(x/y)表示与参数y有关的条件分布密度函数，F1(y)表示分布函数。抽样方法如下：,证明： 抽样效率为：E=1/M,为了实现某个复杂的随机变量 y 的抽样，将其表示成若干个简单的随机变量 x1，x2，xn 的函数得到 x1，x2，xn 的抽样后，即可确定 y 的抽样，这种方法叫作替换法抽样。即,替换抽样方法,例11. 散射方位角余弦分布的抽样,散射方

39、位角在0,2上均匀分布，则其正弦和余弦sin和cos服从如下分布：直接抽样方法为：,令=2，则在0,上均匀分布，作变换其中01，0，则(x,y) 表示上半个单位圆内的点。如果 (x,y) 在上半个单位圆内均匀分布，则在0,上均匀分布，由于,因此抽样sin和cos的问题就变成在上半个单位圆内均匀抽样 (x,y) 的问题。 为获得上半个单位圆内的均匀点，采用挑选法，在上半个单位圆的外切矩形内均匀投点（如图）。舍弃圆外的点，余下的就是所要求的点。抽样方法为：抽样效率E=/40.785,为实现散射方位角余弦分布抽样，最重要的是在上半个单位圆内产生均匀分布点。下面这种方法，首先在单位圆的半个外切正六边形

40、内产生均匀分布点，如图所示。,于是便有了抽样效率更高的抽样方法：抽样效率,例12. 正态分布的抽样,标准正态分布密度函数为：引入一个与标准正态随机变量X独立同分布的随机变量Y，则（X,Y）的联合分布密度为：作变换,则（,）的联合分布密度函数为：由此可知，与相互独立，其分布密度函数分别为分别抽取， ：,从而得到一对服从标准正态分布的随机变量X和Y： 对于一般的正态分布密度函数 N(,2) 的抽样，其抽样结果为：,例13. 分布的抽样,分布密度函数的一般形式为：其中n，k为整数。为了实现分布的抽样，将其看作一组简单的相互独立随机变量的函数，通过这些简单随机变量的抽样，实现分布的抽样。设 x1，x2

41、，xn 为一组相互独立、具有相同分布 F(x) 的随机变量，k为 x1，x2，xn 按大小顺序排列后的第k个，记为：,则k的分布函数为：当 F(x)=x 时，不难验证，k的分布密度函数为分布。因此， 分布的抽样可用如下方法实现：选取n个随机数，按大小顺序排列后取第k个，即,随机抽样的一般方法,加抽样方法 减抽样方法乘抽样方法乘加抽样方法乘减抽样方法对称抽样方法积分抽样方法,加抽样方法,加抽样方法是对如下加分布给出的一种抽样方法： 其中Pn0， ，且 fn(x)为与参数n有关的分布密度函数，n=1，2，。 由复合分布抽样方法可知，加分布的抽样方法为:首先抽样确定n，然后由 fn(x)中抽样x，即

42、：,例14. 多项式分布抽样,多项式分布密度函数的一般形式为： 将 f(x) 改写成如下形式： 则该分布的抽样方法为：,例15. 球壳内均匀分布抽样,设球壳内半径为R0，外半径为R1，点到球心的距离为r，则r的分布密度函数为 分布函数为该分布的直接抽样方法是,为避免开立方根运算，作变换：则 x0,1，其分布密度函数为：其中,则x及r的抽样方法为：,减抽样方法,减抽样方法是对如下形式的分布密度所给出的一种抽样方法： 其中A1、A2为非负实数，f1(x) 、f2(x)均为分布密度函数。 减抽样方法分为以下两种形式：,以上两种形式的抽样方法，究竟选择哪种好，要看f1(x) 、f2(x)哪一个容易抽样

43、，如相差不多，选用第一种方法抽样效率高。,（1）将f (x)表示为 令m表示f2(x)f1(x)的下界，使用挑选法，从f1(x)中抽取Xf1 抽样效率为：,（2）将f (x)表示为 使用挑选法，从f2(x)中抽取Xf2 抽样效率为：,例16. 分布抽样,分布的一个特例： 取A12，A21，f1(x)1，f2(x)2x，此时m0，则根据第一种形式的减抽样方法，有或,由于11可用1代替，该抽样方法可简化为：对于21的情况，可取 Xf1 ，因此与分布的推论相同。,如下形式的分布称为乘分布： 其中H(x)为非负函数， f1(x)为任意分布密度函数。 令M为H(x)的上界，乘抽样方法如下： 抽样效率为：

44、,乘抽样方法,例17. 倒数分布抽样,倒数分布密度函数为： 其直接抽样方法为：下面采用乘抽样方法，考虑如下分布族：其中 i = 1，2，该分布的直接抽样方法为：,利用这一分布族，将倒数分布 f(x) 表示成:其中，乘法分布的抽样方法如下:该分布的抽样效率为：,例18. 麦克斯韦(Maxwell)分布抽样,麦克斯韦分布密度函数的一般形式为： 使用乘抽样方法，令该分布的直接抽样方法为：,此时则麦克斯韦分布的抽样方法为:该分布的抽样效率为：,在实际问题中，经常会遇到如下形式的分布： 其中Hn(x)为非负函数，fn(x) 为任意分布密度函数，n=1，2，。不失一般性，只考虑n=2的情况： 将 f(x)

45、 改写成如下的加分布形式：,乘加抽样方法,其中,乘加抽样方法为： 该方法的抽样效率为：,这种方法需要知道P1的值(P2=1P1），这对有些分布是很困难的。下面的方法可以不用计算P1 ： 对于任意小于1的正数P1 ，令P2=1P1 ；,则采用复合挑选抽样方法，有：,当取 时，抽样效率最高 这时，乘加抽样方法为：,由于 可知第一种方法比第二种方法的抽样效率高。,例19. 光子散射后能量分布的抽样,令光子散射前后的能量分别为 和 （以 m0c2 为单位，m0为电子静止质量，c 为光速）， ，则 x 的分布密度函数为： 该分布即为光子散射能量分布，它是由著名的KlinNishina 公式确定的。其中

46、K() 为归一因子：,把光子散射能量分布改写成如下形式：在1, 1+2上定义如下函数:,则有使用乘加抽样方法：,光子散射能量分布的抽样方法为：该方法的抽样效率为：,乘减分布的形式为：其中H1(x) 、H2(x)为非负函数，f1(x)、f2(x) 为任意分布密度函数。 与减抽样方法类似，乘减分布的抽样方法也分为两种。,乘减抽样方法,（1）将 f (x) 表示为令H1(x)的上界为M1， 的下界为m，使用乘抽样方法得到如下乘减抽样方法：,（2）将 f (x) 表示为令H2(x)的上界为M2，使用乘抽样方法，得到另一种乘减抽样方法：,例20. 裂变中子谱分布抽样,裂变中子谱分布的一般形式为： 其中A

47、，B，C，Emin，Emax 均为与元素有关的量。令其中为归一因子，为任意参数。,相应的 H1(E)，H2(E) 为：于是裂变中子谱分布可以表示成乘减分布形式:容易确定 H1(E) 的上界为：为提高抽样效率，应取使得M1 达到最小，此时,取 m0，令则裂变中子谱分布的抽样方法为:抽样效率,对称分布的一般形式为：其中 f1(x) 为任意分布密度函数，满足偶函数对称条件，H(x) 为任意奇函数，即对任意x满足：对称分布的抽样方法如下：取=21,对称抽样方法,证明：因为=21，x 相当于 ，因此,例21. 质心系各向同性散射角余弦分布抽样,在质心系各向同性散射的假设下，为得到实验室系散射角余弦，需首

48、先抽样确定质心条散射角余弦：再利用下面转换公式:得到实验室系散射角余弦L。其中A为碰撞核质量，C、L 分别为质心系和实验室系散射角。,为避免开方运算，可以使用对称分布抽样。 根据转换公式可得： 依照质心系散射各向同性的假定，可得到实验室系散射角余弦L 的分布如下：该密度函数中的第一项为偶函数，第二项为奇函数，因而是对称分布。其中,从 f1(L) 的抽样可使用挑选法然后再以的概率决定接受或取负值。 上述公式涉及开方运算，需要进一步简化。,注意以下事实：对于任意0a1令则上述挑选抽样中的挑选条件简化为：另一方面，在 即 的条件下，2/a 在1, 1上均匀分布，故可令2/a，则最终决定取正负值的条件

49、简化为：,于是，得到质心系各向同性散射角余弦分布的抽样方法为：,如下形式的分布密度函数称为积分分布密度函数，其中 f0(x，y) 为任意二维分布密度函数，H(x)为任意函数。该分布密度函数的抽样方法为：,积分抽样方法,证明：对于任意x,例22. 各向同性散射方向的抽样,为了确定各向同性散射方向 ，根据公式： 对于各向同性散射，cos在1, 1上均匀分布，在0, 2上均匀分布。由于直接抽样需要计算三角函数和开方。,定义两个随机变量：可以证明，当 时，随机变量 x 和 y 服从如下分布:定义区域为：,则 wcos 的分布可以用上述分布表示成积分分布的形式：令 ，则属于上述积分限内的 y 一定满足条

50、件 。,各向同性散射方向的抽样方法为：抽样效率为：,随机抽样的其它方法,偏倚抽样方法近似抽样方法近似-修正抽样方法多维分布抽样方法指数分布的抽样,使用蒙特卡罗方法计算积分时，可考虑将积分I改写为其中 f *(x) 为一个与 f (x) 有相同定义域的新的分布密度函数。于是可以这样计算积分I：这里 Xi 是从 f *(x) 中抽取的第 i 个子样。,偏移抽样方法,由此可以看出，原来由 f (x) 抽样，现改为由另一个分布密度函数 f *(x) 抽样，并附带一个权重纠偏因子这种方法称为偏倚抽样方法。 从 f (x) 中抽取的 Xf ，满足而对于偏倚抽样，有 一般情况下，Xf 是具有分布 f (x)