《蒙特卡洛在核技术中的应用讲义完整版ppt课件.ppt》由会员分享，可在线阅读，更多相关《蒙特卡洛在核技术中的应用讲义完整版ppt课件.ppt（402页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、林谦第一章第一章 蒙特卡罗方法概述蒙特卡罗方法概述第二章第二章 随机数随机数第三章第三章 由已知分布的随机抽样由已知分布的随机抽样第四章第四章 蒙特卡罗方法解粒子输运问题蒙特卡罗方法解粒子输运问题蒙特卡罗方法在实验核物理中的应用蒙特卡罗方法在实验核物理中的应用许淑艳 编著原子能出版社蒙特卡罗方法蒙特卡罗方法清华大学蒙特卡罗方法及其在粒子输运问题中的应用蒙特卡罗方法及其在粒子输运问题中的应用裴鹿成 张孝泽 编著 科学出版社蒙特卡罗方法蒙特卡罗方法徐钟济 编著上海科学技术出版社电话电话83918电子邮件电子邮件蒙特卡罗方法的基本思想蒙特卡罗方法的基本思想蒙特卡罗方法的收敛性，误差蒙特卡罗方法的收敛
2、性，误差蒙特卡罗方法的特点蒙特卡罗方法的特点蒙特卡罗方法的主要应用范围蒙特卡罗方法的主要应用范围作作 业业 蒙特卡罗方法又称随机抽样技巧或统计试验方法。半个多世纪以来，由于科学技术的发展和电子计算机的发明 ，这种方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。 二十世纪四十年代中期，由于科学技术的发展和电子计算机的发明，蒙特卡罗方法作为一种独立的方法被提出
3、来，并首先在核武器的试验与研制中得到了应用。但其基本思想并非新颖，人们在生产实践和科学试验中就已发现，并加以利用。两个例子 例1. 蒲丰氏问题 例2. 射击问题（打靶游戏）基本思想计算机模拟试验过程 为了求得圆周率值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（ la）的平行线相交的频率代替概率P，再利用准确的关系式： 求出值 其中为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。alP2)(22nNalaPl 一些人进行了实验，其结果列于下表 ：实验者年份投计次数的实验值沃尔弗(Wolf)185050003.159
4、6斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929 设r表示射击运动员的弹着点到靶心的距离，(r)表示击中r处相应的得分数（环数），f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射击水平。该运动员的射击成绩为 用概率语言来说，是随机变量(r)的数学期望，即 )(rgEg 0)()(drrfrgg 现假设该运动员进行了次射击，每次射击的弹着点依次为r1，r2，rN，则次得分g(r1)，g(r2)，g(rN)的算术平均值 代表了该运动员的成绩。换言之，为积分的估计值，或近似值。 在
5、该例中，用次试验所得成绩的算术平均值作为数学期望的估计值（积分近似值）。 NiiNrgNg1)(1 由以上两个例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。 当随机变量的取值仅为1或0时，它的数学期望就是某个事件的概率。或者说，某种事件的概率也是随机变量（仅取值为1或0）的数学期望。 因此，可以通俗地说，蒙特卡罗方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数f(r)的随
6、机变量(r)的数学期望 通过某种试验，得到个观察值r1，r2，rN（用概率语言来说，从分布密度函数f(r)中抽取个子样r1，r2，rN，），将相应的个随机变量的值g(r1)，g(r2)，g(rN)的算术平均值 作为积分的估计值（近似值）。 NiiNrgNg1)(10)()(drrfrgg 为了得到具有一定精确度的近似解，所需试验的次数是很多的，通过人工方法作大量的试验相当困难，甚至是不可能的。因此，蒙特卡罗方法的基本思想虽然早已被人们提出，却很少被使用。本世纪四十年代以来，由于电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验交由计算机完成，使得蒙特卡罗方法得
7、以广泛地应用，在现代化的科学技术中发挥应有的作用。 计算机模拟试验过程，就是将试验过程（如投针，射击）化为数学问题，在计算机上实现。以上述两个问题为例，分别加以说明。 例1. 蒲丰氏问题 例2. 射击问题（打靶游戏） 由上面两个例题看出，蒙特卡罗方法常以一个“概率模型”为基础，按照它所描述的过程，使用由已知分布抽样的方法，得到部分试验结果的观察值，求得问题的近似解。 设针投到地面上的位置可以用一组参数（x,）来描述，x为针中心的坐标，为针与平行线的夹角，如图所示。 任意投针，就是意味着x与都是任意取的，但x的范围限于0，a，夹角的范围限于0，。在此情况下，针与平行线相交的数学条件是针在平行线间
8、的位置 sin lx 如何产生任意的（x,）？x在0，a上任意取值，表示x在0，a上是均匀分布的，其分布密度函数为： 类似地，的分布密度函数为： 因此，产生任意的（x,）的过程就变成了由f1(x)抽样x及由f2()抽样的过程了。由此得到： 其中1，2均为（0,1）上均匀分布的随机变量。 其他, 00,/1)(1axaxf其他, 00,/1)(2f21 ax 每次投针试验，实际上变成在计算机上从两个均匀分布的随机变量中抽样得到（x,），然后定义描述针与平行线相交状况的随机变量s(x,)，为 如果投针次，则 是针与平行线相交概率的估计值。事实上， 于是有 其他当, 0sin, 1),(lxxsNi
9、iiNxsNs1),(1aladxddxdfxfxsPl2)()(),(sin0021NsalaPl22 设射击运动员的弹着点分布为 用计算机作随机试验（射击）的方法为，选取一个随机数，按右边所列方法判断得到成绩。 这样，就进行了一次随机试验（射击），得到了一次成绩 (r)，作次试验后，得到该运动员射击成绩的近似值 环数 78910概率 0.10.10.30.5环中命环命中环命中环命中1095 . 082 . 071 . 0NiiNrgNg1)(1 蒙特卡罗方法作为一种计算方法，其收敛性与误差是普遍关心的一个重要问题。收敛性误差减小方差的各种技巧 效率 由前面介绍可知，蒙特卡罗方法是由随机变量
10、X的简单子样X1，X2，XN的算术平均值： 作为所求解的近似值。由大数定律可知， 如X1，X2，XN独立同分布，且具有有限期望值（E(X)），则 即随机变量X的简单子样的算术平均值 ，当子样数充分大时，以概率1收敛于它的期望值E(X)。NiiNXNX111)(limXEXPNNNX 蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，XN独立同分布，且具有有限非零的方差2 ，即 f(X)是X的分布密度函数。则dtexXEXNPxxtNN2/221)(limdxxfXEx)()(022 当N充分大时，有如下的近似式 其中称为置信度，1称为
11、置信水平。 这表明，不等式 近似地以概率 1成立，且误差收敛速度的阶为 。 通常，蒙特卡罗方法的误差定义为 上式中 与置信度是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出 。122)(02/2dteNXEXPtNNXEXN)()(2/1NON 下面给出几个常用的与的数值： 关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。第二，误差中的均方差是未知的，必须使用其估计值 来代替，在计算所求量的同时，可计算出 。 0.50.050.003 0.67451.9632112)1(1NiiNiiXNXN 显然，当给定置信度
12、后，误差由和N决定。要减小，或者是增大N，或者是减小方差2。在固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。 另一方面，如能减小估计的均方差，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果。因此降低方差的各种技巧，引起了人们的普遍注意。后面课程将会介绍一些降低方差的技巧。 一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就 是蒙特卡罗方法中效率的概念。它定义为 ，其中c 是观察一个子样的平均费用。显然
13、 越小，方法越有效。 c2c2优点能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。受几何条件限制小。收敛速度与问题的维数无关。具有同时计算多个方案与多个未知量的能力。误差容易确定。程序结构简单，易于实现。 缺点收敛速度慢。误差具有概率性。在粒子输运问题中，计算结果与系统大小有关。 从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。 在计算s维空间中的任一区域Ds上的积分 时，无论区域Ds的形状多么特殊，只要能给出描述Ds的几何特征的条件，就可以
14、从Ds中均匀产生N个点 ，得到积分的近似值。 其中Ds为区域Ds的体积。这是数值方法难以作到的。 另外，在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。 ssDdxdxdxxxxggs2121),( ),()()(2)(1isiixxxNiisiisNxxxgNDg1)()(2)(1),( 由误差定义可知，在给定置信水平情况下，蒙特卡罗方法的收敛速度为，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计算时间的变化，不影响误差。也就是说，使用蒙特卡罗方法时，抽取的子样总数N与维数s无关。维数的增加，除了增加相应的计算量外，
15、不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。)(2/1NO 对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。 另外，使用蒙特卡罗方法还可以同时
16、得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。 对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。 一般计算方法常存在着有效位数损失问题，而要解决这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。 在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。 如前所述，蒙特卡罗方法的收敛速度为 ，一般不容易得到精确度较高的近似结果。对于维数少（三维以下）的问
17、题，不如其他方法好。 )(2/1NO 由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。 经验表明，只有当系统的大小与粒子的平均自由程可以相比较时（一般在十个平均自由程左右），蒙特卡罗方法计算的结果较为满意。但对于大系统或小概率事件的计算问题，计算结果往往比真值偏低。而对于大系统，数值方法则是适用的。 因此，在使用蒙特卡罗方法时，可以考虑把蒙特卡罗方法与解析（或数值）方法相结合，取长补短，既能解决解析（或数值）方法难以解决的问题，也可以解决单纯使用蒙特卡罗方法难以解决的问题。这样，可以发挥蒙特卡罗方法的特长，使其应用范围更加广泛。 蒙特卡罗方法所特
18、有的优点，使得它的应用范围越来越广。它的主要应用范围包括：粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。随着科学技术的发展，其应用范围将更加广泛。 蒙特卡罗方法在粒子输运问题中的应用范围主要包括：实验核物理，反应堆物理，高能物理等方面。 蒙特卡罗方法在实验核物理中的应用范围主要包括：通量及反应率，中子探测效率，光子探测效率，光子能量沉积谱及响应函数，气体正比计数管反冲质子谱，多次散射与通量衰减修正等方面。 用蒲丰投针法在计算机上计算值，取a=4、l=3。分别用理论计算和计算机模拟计算，求连续掷两颗骰子，点数之和大于6且第一次掷出的点数大于第二次掷出点数的概
19、率。 随机数的定义及产生方法随机数的定义及产生方法 伪随机数伪随机数产生伪随机数的乘同余方法产生伪随机数的乘同余方法产生伪随机数的乘加同余方法产生伪随机数的乘加同余方法产生伪随机数的其他方法产生伪随机数的其他方法伪随机数序列的均匀性和独立性伪随机数序列的均匀性和独立性作作 业业 由具有已知分布的总体中抽取简单子样，在蒙特卡罗方法中占有非常重要的地位。总体和子样的关系，属于一般和个别的关系，或者说属于共性和个性的关系。由具有已知分布的总体中产生简单子样，就是由简单子样中若干个性近似地反映总体的共性。 随机数是实现由已知分布抽样的基本量，在由已知分布的抽样过程中，将随机数作为已知量，用适当的数学方
20、法可以由它产生具有任意已知分布的简单子样。 随机数的定义及性质随机数表物理方法 在连续型随机变量的分布中，最简单而且最基本的分布是单位均匀分布。由该分布抽取的简单子样称，随机数序列，其中每一个体称为随机数。 单位均匀分布也称为0，1上的均匀分布，其分布密度函数为： 分布函数为 ：其他, 010, 1)(xxf1, 110,0, 0)(xxxxxF 由于随机数在蒙特卡罗方法中占有极其重要的位置，我们用专门的符号表示。由随机数序列的定义可知，1，2，是相互独立且具有相同单位均匀分布的随机数序列。也就是说，独立性独立性、均匀性均匀性是随机数必备的两个特点。 随机数具有非常重要的性质：对于任意自然数s
21、，由s个随机数组成的s维空间上的点（n+1，n+2，n+s）在s维空间的单位立方体Gs上均匀分布，即对任意的ai， 如下等式成立：siai, 2 , 110，siiiinasiaP1), 1,( 其中P()表示事件发生的概率。反之，如果随机变量序列1, 2对于任意自然数s，由s个元素所组成的s维空间上的点（n+1，n+s）在Gs上均匀分布，则它们是随机数序列。 由于随机数在蒙特卡罗方法中所处的特殊地位，它们虽然也属于由具有已知分布的总体中产生简单子样的问题，但就产生方法而言，却有着本质上的差别。 为了产生随机数，可以使用随机数表。随机数表是由0，1，9十个数字组成，每个数字以0.1的等概率出现
22、，数字之间相互独立。这些数字序列叫作随机数字序列。如果要得到n位有效数字的随机数，只需将表中每n个相邻的随机数字合并在一起，且在最高位的前边加上小数点即可。例如，某随机数表的第一行数字为7634258910，要想得到三位有效数字的随机数依次为0.763，0.425，0.891。 因为随机数表需在计算机中占有很大内存，而且也难以满足蒙特卡罗方法对随机数需要量非常大的要求，因此，该方法不适于在计算机上使用。 用物理方法产生随机数的基本原理是：利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。用来作为随机数发生器的物理源主要有两种：一种是根据放射
23、性物质的放射性，另一种是利用计算机的固有噪声。 一般情况下，任意一个随机数在计算机内总是用二进制的数表示的： 其中i（i=1，2，m）或者为0，或者为1。 mm2222211 因此，利用物理方法在计算机上产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为0.5。 用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，需要增加随机数发生器和电路联系等附加设备，费用昂贵。因此，该方法也不适合在计算机上使用。伪随机数伪随机数存在的两个问题伪随机数的周期和最大容量 在计算机上产生随机数最实用、最常见的方法是数学方法，即用如下递
24、推公式： 产生随机数序列。对于给定的初始值1,2，k，确定n+k，=1，2，。经常使用的是k=1的情况，其递推公式为： 对于给定的初始值1，确定n+1，=, )(nknT, 2 , 1),(11nTknnnkn 用数学方法产生的随机数，存在两个问题：递推公式和初始值1,2，k确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。由于随机数序列是由递推公式确定的，而在计算机上所能表示的0，1上的数又是有限的，因此，这种方法产生的随机数序列就不可能不出现无限重复。一旦出现这样的n，n (n n )，使得下面等式成立： 随机数序列便出现了周期性的循环现象。对于k=1的情况，只要有一个随机数
25、重复，其后面的随机数全部重复，这与随机数的要求是不相符的。 kiinin, 2 , 1 由于这两个问题的存在，常称用数学方法产生的随机数为伪随机数。对于以上存在的两个问题，作如下具体分析。 关于第一个问题，不能从本质上加以改变，但只要递推公式选得比较好，随机数间的相互独立性是可以近似满足的。至于第二个问题，则不是本质的。因为用蒙特卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所用随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。 用数学方法产生的伪随机数容易在计算机上得到，可以进行复算，而且不受计算机型号的限制。因此，这种方法虽然存在着一些问题，但仍然被广泛地在计算机
26、上使用，是在计算机上产生伪随机数的主要方法。 发生周期性循环现象的伪随机数的个数称为伪随机数的周期。对于前面介绍的情况，伪随机数的周期为nn。 从伪随机数序列的初始值开始，到出现循环现象为止，所产生的伪随机数的个数称为伪随机数的最大容量。前面的例子中，伪随机数的最大容量为n 。 乘同余方法是由Lehmer在1951年提出来的，它的一般形式是：对于任一初始值x1，伪随机数序列由下面递推公式确定： 其中a为常数。 )(mod,1Mxaxii, 2 , 1,11iMxii 对于任意正整数M，根据数论中的标准分解定理，总可以分解成如下形式： 其中P0=2，P1， Pr表示不同的奇素数，0表示非负整数，
27、1，r表示正整数。a无论取什么值，乘同余方法的最大容量的上界为： 的最小公倍数。其中：rrPPPM1010)()(),()(1010rrPPPM222101)(0200000当当或当PriPPPiiiii, 2 , 1),1()(11 如果a与x1满足如下条件： 对于 ， x1与M互素，则乘同余方法产生的伪随机数序列的最大容量达到最大可能值(M)。 2)8(mod532)4(mod31)2(mod1000当或当当a)(0iiPn)(mod1iinPa 为了便于在计算机上使用，通常取 ：=2s 其中s为计算机中二进制数的最大可能有效位数x1= 奇数 a = 52k+1 其中k为使52k+1在计算
28、机上所能容纳的最大整数，即a为计算机上所能容纳的5的最大奇次幂。一般地，s=32时，a=513；s=48，a=515等。伪随机数序列的最大容量(M)=2s-2 。 乘同余方法是使用的最多、最广的方法，在计算机上被广泛地使用。 产生伪随机数的乘加同余方法是由Rotenberg于1960年提出来的，由于这个方法有很多优点，已成为仅次于乘同余方法产生伪随机数的另一主要方法。 乘加同余方法的一般形式是，对任意初始值x1，伪随机数序列由下面递推公式确定： 其中a和c为常数。 )(mod,1Mcxaxii, 2 , 1,11iMxii 关于乘加同余方法的最大容量问题，有如下结论：如果对于正整数M的所有素数
29、因子P，下式均成立： 当M为4的倍数时，还有下式成立： c与M互素，则乘加同余方法所产生的伪随机数序列的最大容量达到最大可能值M。 )(mod1Pa )4(mod1a 为了便于在计算机上使用，通常取M = 2s 其中s为计算机中二进制数的最大可能有效位数。a = 2b + 1(b2)c = 1 这样在计算中可以使用移位和指令加法，提高计算速度。 取中方法加同余方法 判断伪随机数序列是否满足均匀和相互独立的要求，要靠统计检验的方法实现。对于伪随机数的统计检验，一般包括两大类：均匀性检验和独立性检验。 六十年代初，人们开始用定性的方法研究伪随机数序列的均匀性和独立性问题，简要叙述如下。 这里只考虑
30、伪随机数序列1,2，n全体作为子样时的均匀性问题。其中n为伪随机数序列的最大容量。 对于任意的0 x1，令Nn(x)表示伪随机数序列1,2，n中适合不等式i0 。对该分布的直接抽样方法如下： !)(nePnxPnnn0ii1n0iii!i!,当enXF 掷骰子点数X=n的概率为： 选取随机数，如 则 在等概率的情况下，可使用如下更简单的方法： 其中表示取整数。61)( nXP661nnnXF16FX 中子或光子在介质中发生碰撞时，如介质是由多种元素组成，需要确定碰撞核的种类。假定介质中每种核的宏观总截面分别为1，2，n，则中子或光子与每种核碰撞的概率分别为： 其中t12n。碰撞核种类的确定方法
31、为：产生一个随机数，如果 则中子或光子与第I种核发生碰撞。 niPtii, 2 , 1I1ii1I1iiPP 假设中子与核的反应类型有如下几种:弹性散射，非弹性散射，裂变，吸收，相应的反应截面分别为el，in，f，a。则发生每一种反应类型的概率依次为 ： 其中反应总截面telinfa。 taatfftinintelelPPPP 反应类型的确定方法为：产生一个随机数 收吸裂变非弹性散射弹性散射finelinelelPPPPPP 对于连续型分布，如果分布函数F(x) 的反函数 F1(x)存在，则直接抽样方法是 ：)(1 FXF 在a，b上均匀分布的分布函数为： 则 bxbxaabaxaxxF当当当
32、10)()(abaXF 分布为连续型分布，作为它的一个特例是： 其分布函数为： 则 FX10,2)()(20 xxtdtdttfxFxx10,2)(xxxf 指数分布为连续型分布，其一般形式如下： 其分布函数为： 则 因为1也是随机数，可将上式简化为 0,1)()(0 xedteadttfxFaxxatx0,)(xeaxfax)1ln(1aXFln1aXF 连续性分布函数的直接抽样方法对于分布函数的反函数存在且容易实现的情况，使用起来是很方便的。但是对于以下几种情况，直接抽样法是不合适的。分布函数无法用解析形式给出，因而其反函数也无法给出。分布函数可以给出其解析形式，但是反函数给不出来。分布函
33、数即使能够给出反函数，但运算量很大。 下面叙述的挑选抽样方法是克服这些困难的比较好的方法。 为了实现从己知分布密度函数f(x)抽样，选取与f(x)取值范围相同的分布密度函数h(x)，如果 则挑选抽样方法为：)()(supxhxfMxhfhhXXXhMXf)()( 即从h(x)中抽样xh，以 的概率接受它。 下面证明xf 服从分布密度函数f(x)。证明：对于任意x )()(hhxhMxf)()()()(,)()()(hhhhhhhhfXhMXfPXhMXfdxxXxPXhMXfdxxXxPdxxXxPdxxfdXXfdXXfdXXhXhMXfdXXhXhMXfddXXhddXXhhhdxxxhh
34、hhhhdxxxhhhhXhMXfhhdxxxXhMXfhhhhhh)()()()()()()()()()()()()(0)()(0 使用挑选抽样方法时，要注意以下两点：选取h(x)时要使得h(x)容易抽样且M的值要尽量小。因为M小能提高抽样效率。抽样效率是指在挑选抽样方法中进行挑选时被选中的概率。按此定义，该方法的抽样效率E为： 所以，M越小，抽样效率越高。MdXXhXhMXfXhMXfPEhhhhhh1)()()()()( 当 f(x) 在0，1上定义时，取 h(x)=1，Xh=， 此时挑选抽样方法为)(sup10 xfMxfXMf)( 令圆半径为R0，点到圆心的距离为r，则r的分布密度函
35、数为 分布函数为 容易知道，该分布的直接抽样方法是其它当002)(020RrRrrf202)(RrrF0Rrf 由于开方运算在计算机上很费时间，该方法不是好方法。下面使用挑选抽样方法：取 则抽样框图为 00022)()(1)(RrMRrrhrfRrhh，2021Rrf 显然，没有必要舍弃12的情况，此时，只需取 就可以了，亦即 另一方面，也可证明 与 具有相同的分布 。 10 Rrf),max(210 Rrf),max(212)(rrF 在实际问题中，经常有这样的随机变量，它服从的分布与一个参数有关，而该参数也是一个服从确定分布的随机变量，称这样的随机变量服从复合分布。例如，分布密度函数 是一
36、个复合分布。其中Pn0，n=1，2，且 fn(x)为与参数n有关的分布密度函数，n=1，2， 参数n服从如下分布1)()(nnnxfPxf11nnPynnPyF)( 复合分布的一般形式为： 其中f2(x/y)表示与参数y有关的条件分布密度函数， F1(y)表示分布函数。 复合分布的抽样方法为:首先由分布函数F1(y) 或分布密度函数f1(y)中抽样YF1或Yf1，然后再由分布密度函数f2(x/ YF1)中抽样确定Xf2 (x/YF) 证明： 所以，Xf所服从的分布为f (x)。)()()(12ydFyxfxf)/(12FYxffXXdxxfYdxdFYxfdxxXxpdxxXxpFYxff)(
37、)()()()(12)/(12 指数函数分布的一般形式为： 引入如下两个分布密度函数：其它当00)(1xdyyenxEnxyn其它当其它当00)(01)(211xeyyxfyynyfxyn 则 使用复合抽样方法，首先从f1(y)中抽取y 再由f2(x/ YF1)中抽取x 112)()()(dyyfyxfxEn),max(11211nnfY1211ln),max(ln1nnfnfYX 考虑另一种形式的复合分布如下： 其中0H(x,y)M，f2(x/y)表示与参数y有关的条件分布密度函数，F1(y)表示分布函数。抽样方法如下：)()(),()(12ydFyxfyxHxf)/()/(12112),(
38、FFYxffFYxfXXMYXH 证明： 抽样效率为：E=1/MdxxfdxydFyxfyxHydxdFyxfMyxHydxdFyxfMyxHdydxdFyxfdydxdFyxfMYXHPMYXHdxxXxPMYXHdxxXxPdxxXxPdxxxMyxHdxxxMyxHFfFffFfff)()()(),()()(),()()(),()()()()(),(),(,),()(121212),(012),(01212122122 为了实现某个复杂的随机变量 y 的抽样，将其表示成若干个简单的随机变量 x1，x2，xn 的函数得到 x1，x2，xn 的抽样后，即可确定 y 的抽样，这种方法叫作替换法
39、抽样。即),(21nxxxgy),(21nfXXXgY 散射方位角在0,2上均匀分布，则其正弦和余弦sin和cos服从如下分布：直接抽样方法为：其它当011111)(2xxxf2coscos2sinsin令=2，则在0,上均匀分布，作变换其中01，0，则(x,y) 表示上半个单位圆内的点。如果 (x,y) 在上半个单位圆内均匀分布，则在0,上均匀分布，由于sincosyx2222sincosyxyyxx222222222cossin22sinsinsincos2coscosyxxyyxyx 因此抽样sin和cos的问题就变成在上半个单位圆内均匀抽样 (x,y) 的问题。 为获得上半个单位圆内的
40、均匀点，采用挑选法，在上半个单位圆的外切矩形内均匀投点（如图）。舍弃圆外的点，余下的就是所要求的点。抽样方法为：抽样效率E=/40.78521yx2221212221222122212sin,cos1 为实现散射方位角余弦分布抽样，最重要的是在上半个单位圆内产生均匀分布点。下面这种方法，首先在单位圆的半个外切正六边形内产生均匀分布点，如图所示。 于是便有了抽样效率更高的抽样方法：抽样效率222121222122212221221121332sin,33cos131,123906. 032E 标准正态分布密度函数为：引入一个与标准正态随机变量X独立同分布的随机变量Y，则（X,Y）的联合分布密度为
41、：作变换2221)(xexf2)(2221),(yxeyxfsincosyx则（,）的联合分布密度函数为：由此可知，与相互独立，其分布密度函数分别为分别抽取， ：222),(ef21)()(2212ffe212ln2从而得到一对服从标准正态分布的随机变量X和Y： 对于一般的正态分布密度函数 N(,2) 的抽样，其抽样结果为：ffffYYXX)2sin(ln2)2cos(ln22121ffYX 分布密度函数的一般形式为：其中n，k为整数。为了实现分布的抽样，将其看作一组简单的相互独立随机变量的函数，通过这些简单随机变量的抽样，实现分布的抽样。设 x1，x2，xn 为一组相互独立、具有相同分布 F
42、(x) 的随机变量，k为 x1，x2，xn 按大小顺序排列后的第k个，记为：10)1 ()!()!1(!)(1xxxknknxfknk),(21nkkxxxR则k的分布函数为：当 F(x)=x 时，不难验证，k的分布密度函数为分布。因此， 分布的抽样可用如下方法实现：选取n个随机数，按大小顺序排列后取第k个，即 ininkiinxFxFCxFk)(1)()(ininkiinxxCxFk)1 ()(),(21nkfRX加抽样方法 减抽样方法乘抽样方法乘加抽样方法乘减抽样方法对称抽样方法积分抽样方法 加抽样方法是对如下加分布给出的一种抽样方法： 其中Pn0， ，且 fn(x)为与参数n有关的分布密
43、度函数，n=1，2，。 由复合分布抽样方法可知，加分布的抽样方法为:首先抽样确定n，然后由 fn(x)中抽样x，即：1)()(nnnxfPxf11nnPn1nn1n1nnPP,当nffXX 多项式分布密度函数的一般形式为： 将 f(x) 改写成如下形式： 则该分布的抽样方法为：0)(iiixaxf00)() 1(1)(iiiiiixfPxiiaxfn0ii1n0ii11PP),max(当nfX 设球壳内半径为R0，外半径为R1，点到球心的距离为r，则r的分布密度函数为 分布函数为该分布的直接抽样方法是其它当03)(1030312RrRRRrrf3031303)(RRRrrF31303031)(
44、RRRrf为避免开立方根运算，作变换：则 x0,1，其分布密度函数为：其中001)(RxRRr132)(33)()(200102201RxRRRxRRxf211020RRRR则x及r的抽样方法为：001432322101201)(),max(),max(33RXRRrXXXRRRfffff 减抽样方法是对如下形式的分布密度所给出的一种抽样方法： 其中A1、A2为非负实数，f1(x) 、f2(x)均为分布密度函数。 减抽样方法分为以下两种形式：)()()(2211xfAxfAxf 以上两种形式的抽样方法，究竟选择哪种好，要看f1(x) 、f2(x)哪一个容易抽样，如相差不多，选用第一种方法抽样效
45、率高。 （1）将f (x)表示为 令m表示f2(x)f1(x)的下界，使用挑选法，从f1(x)中抽取Xf1 抽样效率为：)()()()(12211xfxfAAxfxf111)()(12212211ffffXXXfXfmAAAmAAA211mAAE （2）将f (x)表示为 使用挑选法，从f2(x)中抽取Xf2 抽样效率为：22112)()()()(AxfxfAxfxfmEmAAmE2122221221211)()(ffffXXmAAmAXfXfmAAmA分布的一个特例： 取A12，A21，f1(x)1，f2(x)2x，此时m0，则根据第一种形式的减抽样方法，有或 10),1 (2)(xxxf2
46、211fX2121fX由于11可用1代替，该抽样方法可简化为：对于21的情况，可取 Xf1 ，因此与分布的推论相同。212fX),min(21fX 如下形式的分布称为乘分布： 其中H(x)为非负函数， f1(x)为任意分布密度函数。 令M为H(x)的上界，乘抽样方法如下： 抽样效率为：)()()(1xfxHxfME111)(fffXXMXH倒数分布密度函数为： 其直接抽样方法为：下面采用乘抽样方法，考虑如下分布族：其中 i = 1，2，该分布的直接抽样方法为：axxaxf1,1ln1)(afeaXlniifaXi 1) 1(1axxxaixfiii1,) 1(1)(11利用这一分布族，将倒数分
47、布 f(x) 表示成:其中，乘法分布的抽样方法如下:该分布的抽样效率为：)()()(xfxHxfi,1)(,ln) 1(,ln) 1()(1111iiiixMxHaaiMxaaixHiifiaXa 1) 1(1 1) 1(21211) 1(ln1iaiaE麦克斯韦分布密度函数的一般形式为： 使用乘抽样方法，令该分布的直接抽样方法为：0,2)(23xexxfx2ln231fX0,32)(321xexfx此时则麦克斯韦分布的抽样方法为:该分布的抽样效率为：eMxexxHx2270,3)(312122221ln23lnfXe795. 0272eE 在实际问题中，经常会遇到如下形式的分布： 其中Hn(
48、x)为非负函数，fn(x) 为任意分布密度函数，n=1，2，。不失一般性，只考虑n=2的情况： 将 f(x) 改写成如下的加分布形式：1)()()(nnnxfxHxf)()()()()(2211xfxHxfxHxf)()()()()()()(*22*1122221111xfPxfPxfPxHPxfPxHPxf 其中)()()()()()()()()()(222*2111*1222111xfPxHxfxfPxHxfdxxfxHPdxxfxHP 乘加抽样方法为： 该方法的抽样效率为：22222)(fffXXMXH11112)(fffXXMXH11P2221212221111MPMPMPPMPPE
49、这种方法需要知道P1的值(P2=1P1），这对有些分布是很困难的。下面的方法可以不用计算P1 ： 对于任意小于1的正数P1 ，令P2=1P1 ；ynnPyFyxfyxfyxfyPxHyPxHyxH)(2),(21),()(2)(21,)(),(12122211),min(),max(22112211MPMPEPMPMM 则采用复合挑选抽样方法，有： 当取 时，抽样效率最高 这时，乘加抽样方法为：22222)(fffXXMXH11112)(fffXXMXH2111MMM2121MME21222111MMMPMMMP 由于 可知第一种方法比第二种方法的抽样效率高。0)()()(2)()()()()
50、(121212211221212121222121222121212221212221212221212122211212122122212121MMMMPMPMMMMMPPMMPMPMMMMMMMPPMMPMPMMMMMMMPMMMPMMMMMMPMPEE令光子散射前后的能量分别为 和 （以 m0c2 为单位，m0为电子静止质量，c 为光速）， ，则 x 的分布密度函数为： 该分布即为光子散射能量分布，它是由著名的KlinNishina 公式确定的。其中 K() 为归一因子：211,1111)(1)(322xxxxxxKxfx22)21 (21421)21ln() 1(21)(K把光子散射能