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

蒙特卡罗方法在核技术中的应用,林谦,目录,第一章蒙特卡罗方法概述第二章随机数第三章由已知分布的随机抽样第四章蒙特卡罗方法解粒子输运问题,教材,蒙特卡罗方法在实验核物理中的应用许淑艳编著原子能出版社蒙特卡罗方法清华大学,参考书,蒙特卡罗方法及其在粒子输运问题中的应用裴鹿成张孝泽编著科学出版社蒙特卡罗方法徐钟济编著上海科学技术出版社,联系方式,电话83918电子邮件linqian,第一章蒙特卡罗方法概述,蒙特卡罗方法的基本思想蒙特卡罗方法的收敛性，误差蒙特卡罗方法的特点蒙特卡罗方法的主要应用范围作业,第一章蒙特卡罗方法概述,蒙特卡罗方法又称随机抽样技巧或统计试验方法。半个多世纪以来，由于科学技术的发展和电子计算机的发明，这种方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。,蒙特卡罗方法的基本思想,二十世纪四十年代中期，由于科学技术的发展和电子计算机的发明，蒙特卡罗方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。但其基本思想并非新颖，人们在生产实践和科学试验中就已发现，并加以利用。两个例子例1.蒲丰氏问题例2.射击问题（打靶游戏）基本思想计算机模拟试验过程,例1.蒲丰氏问题,为了求得圆周率值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（la）的平行线相交的频率代替概率P，再利用准确的关系式：求出值其中为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。,一些人进行了实验，其结果列于下表：,例2.射击问题（打靶游戏）,设r表示射击运动员的弹着点到靶心的距离，(r)表示击中r处相应的得分数（环数），f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射击水平。该运动员的射击成绩为用概率语言来说，是随机变量(r)的数学期望，即,现假设该运动员进行了次射击，每次射击的弹着点依次为r1，r2，rN，则次得分g(r1)，g(r2)，g(rN)的算术平均值代表了该运动员的成绩。换言之，为积分的估计值，或近似值。在该例中，用次试验所得成绩的算术平均值作为数学期望的估计值（积分近似值）。,基本思想,由以上两个例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。当随机变量的取值仅为1或0时，它的数学期望就是某个事件的概率。或者说，某种事件的概率也是随机变量（仅取值为1或0）的数学期望。,因此，可以通俗地说，蒙特卡罗方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数f(r)的随机变量(r)的数学期望通过某种试验，得到个观察值r1，r2，rN（用概率语言来说，从分布密度函数f(r)中抽取个子样r1，r2，rN，），将相应的个随机变量的值g(r1)，g(r2)，g(rN)的算术平均值作为积分的估计值（近似值）。,为了得到具有一定精确度的近似解，所需试验的次数是很多的，通过人工方法作大量的试验相当困难，甚至是不可能的。因此，蒙特卡罗方法的基本思想虽然早已被人们提出，却很少被使用。本世纪四十年代以来，由于电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验交由计算机完成，使得蒙特卡罗方法得以广泛地应用，在现代化的科学技术中发挥应有的作用。,计算机模拟试验过程,计算机模拟试验过程，就是将试验过程（如投针，射击）化为数学问题，在计算机上实现。以上述两个问题为例，分别加以说明。例1.蒲丰氏问题例2.射击问题（打靶游戏）由上面两个例题看出，蒙特卡罗方法常以一个“概率模型”为基础，按照它所描述的过程，使用由已知分布抽样的方法，得到部分试验结果的观察值，求得问题的近似解。,例蒲丰氏问题,设针投到地面上的位置可以用一组参数（x,）来描述，x为针中心的坐标，为针与平行线的夹角，如图所示。任意投针，就是意味着x与都是任意取的，但x的范围限于0，a，夹角的范围限于0，。在此情况下，针与平行线相交的数学条件是,针在平行线间的位置,如何产生任意的（x,）？x在0，a上任意取值，表示x在0，a上是均匀分布的，其分布密度函数为：类似地，的分布密度函数为：因此，产生任意的（x,）的过程就变成了由f1(x)抽样x及由f2()抽样的过程了。由此得到：其中1，2均为（0,1）上均匀分布的随机变量。,每次投针试验，实际上变成在计算机上从两个均匀分布的随机变量中抽样得到（x,），然后定义描述针与平行线相交状况的随机变量s(x,)，为如果投针次，则是针与平行线相交概率的估计值。事实上，于是有,例射击问题,设射击运动员的弹着点分布为用计算机作随机试验（射击）的方法为，选取一个随机数，按右边所列方法判断得到成绩。这样，就进行了一次随机试验（射击），得到了一次成绩(r)，作次试验后，得到该运动员射击成绩的近似值,蒙特卡罗方法的收敛性，误差,蒙特卡罗方法作为一种计算方法，其收敛性与误差是普遍关心的一个重要问题。收敛性误差减小方差的各种技巧效率,收敛性,由前面介绍可知，蒙特卡罗方法是由随机变量X的简单子样X1，X2，XN的算术平均值：作为所求解的近似值。由大数定律可知，如X1，X2，XN独立同分布，且具有有限期望值（E(X)），则即随机变量X的简单子样的算术平均值，当子样数充分大时，以概率1收敛于它的期望值E(X)。,误差,蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，XN独立同分布，且具有有限非零的方差2，即f(X)是X的分布密度函数。则,当N充分大时，有如下的近似式其中称为置信度，1称为置信水平。这表明，不等式近似地以概率1成立，且误差收敛速度的阶为。通常，蒙特卡罗方法的误差定义为上式中与置信度是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出。,下面给出几个常用的与的数值：关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。第二，误差中的均方差是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出。,减小方差的各种技巧,显然，当给定置信度后，误差由和N决定。要减小，或者是增大N，或者是减小方差2。在固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。另一方面，如能减小估计的均方差，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果。因此降低方差的各种技巧，引起了人们的普遍注意。后面课程将会介绍一些降低方差的技巧。,效率,一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中效率的概念。它定义为，其中c是观察一个子样的平均费用。显然越小，方法越有效。,蒙特卡罗方法的特点,优点能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。受几何条件限制小。收敛速度与问题的维数无关。具有同时计算多个方案与多个未知量的能力。误差容易确定。程序结构简单，易于实现。,缺点收敛速度慢。误差具有概率性。在粒子输运问题中，计算结果与系统大小有关。,能够比较逼真地描述具有随机性质的事物的特点及物理实验过程,从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。,受几何条件限制小,在计算s维空间中的任一区域Ds上的积分时，无论区域Ds的形状多么特殊，只要能给出描述Ds的几何特征的条件，就可以从Ds中均匀产生N个点，得到积分的近似值。其中Ds为区域Ds的体积。这是数值方法难以作到的。另外，在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。,收敛速度与问题的维数无关,由误差定义可知，在给定置信水平情况下，蒙特卡罗方法的收敛速度为，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计算时间的变化，不影响误差。也就是说，使用蒙特卡罗方法时，抽取的子样总数N与维数s无关。维数的增加，除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。,具有同时计算多个方案与多个未知量的能力,对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。另外，使用蒙特卡罗方法还可以同时得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。,误差容易确定,对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。一般计算方法常存在着有效位数损失问题，而要解决这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。,程序结构简单，易于实现,在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。,收敛速度慢,如前所述，蒙特卡罗方法的收敛速度为，一般不容易得到精确度较高的近似结果。对于维数少（三维以下）的问题，不如其他方法好。,误差具有概率性,由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。,在粒子输运问题中，计算结果与系统大小有关,经验表明，只有当系统的大小与粒子的平均自由程可以相比较时（一般在十个平均自由程左右），蒙特卡罗方法计算的结果较为满意。但对于大系统或小概率事件的计算问题，计算结果往往比真值偏低。而对于大系统，数值方法则是适用的。因此，在使用蒙特卡罗方法时，可以考虑把蒙特卡罗方法与解析（或数值）方法相结合，取长补短，既能解决解析（或数值）方法难以解决的问题，也可以解决单纯使用蒙特卡罗方法难以解决的问题。这样，可以发挥蒙特卡罗方法的特长，使其应用范围更加广泛。,蒙特卡罗方法的主要应用范围,蒙特卡罗方法所特有的优点，使得它的应用范围越来越广。它的主要应用范围包括：粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。随着科学技术的发展，其应用范围将更加广泛。蒙特卡罗方法在粒子输运问题中的应用范围主要包括：实验核物理，反应堆物理，高能物理等方面。蒙特卡罗方法在实验核物理中的应用范围主要包括：通量及反应率，中子探测效率，光子探测效率，光子能量沉积谱及响应函数，气体正比计数管反冲质子谱，多次散射与通量衰减修正等方面。,作业,用蒲丰投针法在计算机上计算值，取a=4、l=3。分别用理论计算和计算机模拟计算，求连续掷两颗骰子，点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率。,第二章随机数,随机数的定义及产生方法伪随机数产生伪随机数的乘同余方法产生伪随机数的乘加同余方法产生伪随机数的其他方法伪随机数序列的均匀性和独立性作业,第二章随机数,由具有已知分布的总体中抽取简单子样，在蒙特卡罗方法中占有非常重要的地位。总体和子样的关系，属于一般和个别的关系，或者说属于共性和个性的关系。由具有已知分布的总体中产生简单子样，就是由简单子样中若干个性近似地反映总体的共性。随机数是实现由已知分布抽样的基本量，在由已知分布的抽样过程中，将随机数作为已知量，用适当的数学方法可以由它产生具有任意已知分布的简单子样。,随机数的定义及产生方法,随机数的定义及性质随机数表物理方法,随机数的定义及性质,在连续型随机变量的分布中，最简单而且最基本的分布是单位均匀分布。由该分布抽取的简单子样称，随机数序列，其中每一个体称为随机数。单位均匀分布也称为0，1上的均匀分布，其分布密度函数为：分布函数为：,由于随机数在蒙特卡罗方法中占有极其重要的位置，我们用专门的符号表示。由随机数序列的定义可知，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上均匀分布，则它们是随机数序列。由于随机数在蒙特卡罗方法中所处的特殊地位，它们虽然也属于由具有已知分布的总体中产生简单子样的问题，但就产生方法而言，却有着本质上的差别。,随机数表,为了产生随机数，可以使用随机数表。随机数表是由0，1，9十个数字组成，每个数字以0.1的等概率出现，数字之间相互独立。这些数字序列叫作随机数字序列。如果要得到n位有效数字的随机数，只需将表中每n个相邻的随机数字合并在一起，且在最高位的前边加上小数点即可。例如，某随机数表的第一行数字为7634258910，要想得到三位有效数字的随机数依次为0.763，0.425，0.891。因为随机数表需在计算机中占有很大内存，而且也难以满足蒙特卡罗方法对随机数需要量非常大的要求，因此，该方法不适于在计算机上使用。,物理方法,用物理方法产生随机数的基本原理是：利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。用来作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。一般情况下，任意一个随机数在计算机内总是用二进制的数表示的：其中i（i=1，2，m）或者为0，或者为1。,因此，利用物理方法在计算机上产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为0.5。用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，需要增加随机数发生器和电路联系等附加设备，费用昂贵。因此，该方法也不适合在计算机上使用。,伪随机数,伪随机数伪随机数存在的两个问题伪随机数的周期和最大容量,伪随机数,在计算机上产生随机数最实用、最常见的方法是数学方法，即用如下递推公式：产生随机数序列。对于给定的初始值1,2，k，确定n+k，=1，2，。经常使用的是k=1的情况，其递推公式为：对于给定的初始值1，确定n+1，=,伪随机数存在的两个问题,用数学方法产生的随机数，存在两个问题：递推公式和初始值1,2，k确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。由于随机数序列是由递推公式确定的，而在计算机上所能表示的0，1上的数又是有限的，因此，这种方法产生的随机数序列就不可能不出现无限重复。一旦出现这样的n，n(nn)，使得下面等式成立：随机数序列便出现了周期性的循环现象。对于k=1的情况，只要有一个随机数重复，其后面的随机数全部重复，这与随机数的要求是不相符的。,由于这两个问题的存在，常称用数学方法产生的随机数为伪随机数。对于以上存在的两个问题，作如下具体分析。关于第一个问题，不能从本质上加以改变，但只要递推公式选得比较好，随机数间的相互独立性是可以近似满足的。至于第二个问题，则不是本质的。因为用蒙特卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所用随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。用数学方法产生的伪随机数容易在计算机上得到，可以进行复算，而且不受计算机型号的限制。因此，这种方法虽然存在着一些问题，但仍然被广泛地在计算机上使用，是在计算机上产生伪随机数的主要方法。,伪随机数的周期和最大容量,发生周期性循环现象的伪随机数的个数称为伪随机数的周期。对于前面介绍的情况，伪随机数的周期为nn。从伪随机数序列的初始值开始，到出现循环现象为止，所产生的伪随机数的个数称为伪随机数的最大容量。前面的例子中，伪随机数的最大容量为n。,产生伪随机数的乘同余方法,乘同余方法是由Lehmer在1951年提出来的，它的一般形式是：对于任一初始值x1，伪随机数序列由下面递推公式确定：其中a为常数。,乘同余方法的最大容量的上界,对于任意正整数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时，a=513；s=48，a=515等。伪随机数序列的最大容量(M)=2s-2。乘同余方法是使用的最多、最广的方法，在计算机上被广泛地使用。,产生伪随机数的乘加同余方法,产生伪随机数的乘加同余方法是由Rotenberg于1960年提出来的，由于这个方法有很多优点，已成为仅次于乘同余方法产生伪随机数的另一主要方法。乘加同余方法的一般形式是，对任意初始值x1，伪随机数序列由下面递推公式确定：其中a和c为常数。,乘加同余方法的最大容量,关于乘加同余方法的最大容量问题，有如下结论：如果对于正整数M的所有素数因子P，下式均成立：当M为4的倍数时，还有下式成立：c与M互素，则乘加同余方法所产生的伪随机数序列的最大容量达到最大可能值M。,M，x1，a，c的取值,为了便于在计算机上使用，通常取M=2s其中s为计算机中二进制数的最大可能有效位数。a=2b+1(b2)c=1这样在计算中可以使用移位和指令加法，提高计算速度。,产生伪随机数的其他方法,取中方法加同余方法,伪随机数序列的均匀性和独立性,判断伪随机数序列是否满足均匀和相互独立的要求，要靠统计检验的方法实现。对于伪随机数的统计检验，一般包括两大类：均匀性检验和独立性检验。六十年代初，人们开始用定性的方法研究伪随机数序列的均匀性和独立性问题，简要叙述如下。,伪随机数的均匀性,这里只考虑伪随机数序列1,2，n全体作为子样时的均匀性问题。其中n为伪随机数序列的最大容量。对于任意的0x1，令Nn(x)表示伪随机数序列1,2，n中适合不等式i0。对该分布的直接抽样方法如下：,例3.掷骰子点数的抽样,掷骰子点数X=n的概率为：选取随机数，如则在等概率的情况下，可使用如下更简单的方法：其中表示取整数。,例4.碰撞核种类的确定,中子或光子在介质中发生碰撞时，如介质是由多种元素组成，需要确定碰撞核的种类。假定介质中每种核的宏观总截面分别为1，2，n，则中子或光子与每种核碰撞的概率分别为：其中t12n。碰撞核种类的确定方法为：产生一个随机数，如果则中子或光子与第I种核发生碰撞。,例5.中子与核的反应类型的确定,假设中子与核的反应类型有如下几种:弹性散射，非弹性散射，裂变，吸收，相应的反应截面分别为el，in，f，a。则发生每一种反应类型的概率依次为：其中反应总截面telinfa。,反应类型的确定方法为：产生一个随机数,连续型分布的直接抽样方法,对于连续型分布，如果分布函数F(x)的反函数F1(x)存在，则直接抽样方法是：,例6.在a，b上均匀分布的抽样,在a，b上均匀分布的分布函数为：则,例7.分布,分布为连续型分布，作为它的一个特例是：其分布函数为：则,例8.指数分布,指数分布为连续型分布，其一般形式如下：其分布函数为：则因为1也是随机数，可将上式简化为,连续性分布函数的直接抽样方法对于分布函数的反函数存在且容易实现的情况，使用起来是很方便的。但是对于以下几种情况，直接抽样法是不合适的。分布函数无法用解析形式给出，因而其反函数也无法给出。分布函数可以给出其解析形式，但是反函数给不出来。分布函数即使能够给出反函数，但运算量很大。下面叙述的挑选抽样方法是克服这些困难的比较好的方法。,挑选抽样方法,为了实现从己知分布密度函数f(x)抽样，选取与f(x)取值范围相同的分布密度函数h(x)，如果则挑选抽样方法为：,即从h(x)中抽样xh，以的概率接受它。下面证明xf服从分布密度函数f(x)。证明：对于任意x,使用挑选抽样方法时，要注意以下两点：选取h(x)时要使得h(x)容易抽样且M的值要尽量小。因为M小能提高抽样效率。抽样效率是指在挑选抽样方法中进行挑选时被选中的概率。按此定义，该方法的抽样效率E为：所以，M越小，抽样效率越高。,当f(x)在0，1上定义时，取h(x)=1，Xh=，此时挑选抽样方法为,例9.圆内均匀分布抽样,令圆半径为R0，点到圆心的距离为r，则r的分布密度函数为分布函数为容易知道，该分布的直接抽样方法是,由于开方运算在计算机上很费时间，该方法不是好方法。下面使用挑选抽样方法：取则抽样框图为,显然，没有必要舍弃12的情况，此时，只需取就可以了，亦即另一方面，也可证明与具有相同的分布。,复合抽样方法,在实际问题中，经常有这样的随机变量，它服从的分布与一个参数有关，而该参数也是一个服从确定分布的随机变量，称这样的随机变量服从复合分布。例如，分布密度函数是一个复合分布。其中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)。,例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.散射方位角余弦分布的抽样,散射方位角在0,2上均匀分布，则其正弦和余弦sin和cos服从如下分布：直接抽样方法为：,令=2，则在0,上均匀分布，作变换其中01，0，则(x,y)表示上半个单位圆内的点。如果(x,y)在上半个单位圆内均匀分布，则在0,上均匀分布，由于,因此抽样sin和cos的问题就变成在上半个单位圆内均匀抽样(x,y)的问题。为获得上半个单位圆内的均匀点，采用挑选法，在上半个单位圆的外切矩形内均匀投点（如图）。舍弃圆外的点，余下的就是所要求的点。抽样方法为：抽样效率E=/40.785,为实现散射方位角余弦分布抽样，最重要的是在上半个单位圆内产生均匀分布点。下面这种方法，首先在单位圆的半个外切正六边形内产生均匀分布点，如图所示。,于是便有了抽样效率更高的抽样方法：抽样效率,例12.正态分布的抽样,标准正态分布密度函数为：引入一个与标准正态随机变量X独立同分布的随机变量Y，则（X,Y）的联合分布密度为：作变换,则（,）的联合分布密度函数为：由此可知，与相互独立，其分布密度函数分别为分别抽取，：,从而得到一对服从标准正态分布的随机变量X和Y：对于一般的正态分布密度函数N(,2)的抽样，其抽样结果为：,例13.分布的抽样,分布密度函数的一般形式为：其中n，k为整数。为了实现分布的抽样，将其看作一组简单的相互独立随机变量的函数，通过这些简单随机变量的抽样，实现分布的抽样。设x1，x2，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，即：,例14.多项式分布抽样,多项式分布密度函数的一般形式为：将f(x)改写成如下形式：则该分布的抽样方法为：,例15.球壳内均匀分布抽样,设球壳内半径为R0，外半径为R1，点到球心的距离为r，则r的分布密度函数为分布函数为该分布的直接抽样方法是,为避免开立方根运算，作变换：则x0,1，其分布密度函数为：其中,则x及r的抽样方法为：,减抽样方法,减抽样方法是对如下形式的分布密度所给出的一种抽样方法：其中A1、A2为非负实数，f1(x)、f2(x)均为分布密度函数。减抽样方法分为以下两种形式：,以上两种形式的抽样方法，究竟选择哪种好，要看f1(x)、f2(x)哪一个容易抽样，如相差不多，选用第一种方法抽样效率高。,（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)的上界，乘抽样方法如下：抽样效率为：,乘抽样方法,例17.倒数分布抽样,倒数分布密度函数为：其直接抽样方法为：下面采用乘抽样方法，考虑如下分布族：其中i=1，2，该分布的直接抽样方法为：,利用这一分布族，将倒数分布f(x)表示成:其中，乘法分布的抽样方法如下:该分布的抽样效率为：,例18.麦克斯韦(Maxwell)分布抽样,麦克斯韦分布密度函数的一般形式为：使用乘抽样方法，令该分布的直接抽样方法为：,此时则麦克斯韦分布的抽样方法为:该分布的抽样效率为：,在实际问题中，经常会遇到如下形式的分布：其中Hn(x)为非负函数，fn(x)为任意分布密度函数，n=1，2，。不失一般性，只考虑n=2的情况：将f(x)改写成如下的加分布形式：,乘加抽样方法,其中,乘加抽样方法为：该方法的抽样效率为：,这种方法需要知道P1的值(P2=1P1），这对有些分布是很困难的。下面的方法可以不用计算P1：对于任意小于1的正数P1，令P2=1P1；,则采用复合挑选抽样方法，有：,当取时，抽样效率最高这时，乘加抽样方法为：,由于可知第一种方法比第二种方法的抽样效率高。,例19.光子散射后能量分布的抽样,令光子散射前后的能量分别为和（以m0c2为单位，m0为电子静止质量，c为光速），则x的分布密度函数为：该分布即为光子散射能量分布，它是由著名的KlinNishina公式确定的。其中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，B，C，Emin，Emax均为与元素有关的量。令其中为归一因子，为任意参数。,相应的H1(E)，H2(E)为：于是裂变中子谱分布可以表示成乘减分布形式:容易确定H1(E)的上界为：为提高抽样效率，应取使得M1达到最小，此时,取m0，令则裂变中子谱分布的抽样方法为:抽样效率,对称分布的一般形式为：其中f1(x)为任意分布密度函数，满足偶函数对称条件，H(x)为任意奇函数，即对任意x满足：对称分布的抽样方法如下：取=21,对称抽样方法,证明：因为=21，x相当于，因此,例21.质心系各向同性散射角余弦分布抽样,在质心系各向同性散射的假设下，为得到实验室系散射角余弦，需首先抽样确定质心条散射角余弦：再利用下面转换公式:得到实验室系散射角余弦L。其中A为碰撞核质量，C、L分别为质心系和实验室系散射角。,为避免开方运算，可以使用对称分布抽样。根据转换公式可得：依照质心系散射各向同性的假定，可得到实验室系散射角余弦L的分布如下：该密度函数中的第一项为偶函数，第二项为奇函数，因而是对称分布。其中,从f1(L)的抽样可使用挑选法然后再以的概率决定接受或取负值。上述公式涉及开方运算，需要进一步简化。,注意以下事实：对于任意0a1令则上述挑选抽样中的挑选条件简化为：另一方面，在即的条件下，2/a在1,1上均匀分布，故可令2/a，则最终决定取正负值的条件简化为：,于是，得到质心系各向同性散射角余弦分布的抽样方法为：,如下形式的分布密度函数称为积分分布密度函数，其中f0(x，y)为任意二维分布密度函数，H(x)为任意函数。该分布密度函数的抽样方法为：,积分抽样方法,证明：对于任意x,例22.各向同性散射方向的抽样,为了确定各向同性散射方向，根据公式：对于各向同性散射，cos在1,1上均匀分布，在0,2上均匀分布。由于直接抽样需要计算三角函数和开方。,定义两个随机变量：可以证明，当时，随机变量x和y服从如下分布:定义区域为：,则wcos的分布可以用上述分布表示成积分分布的形式：令，则属于上述积分限内的y一定满足条件。,各向同性散射方向的抽样方法为：抽样效率为：,随机抽样的其它方法,偏倚抽样方法近似抽样方法近似-修正抽样方法多维分布抽样方法指数分布的抽样,使用蒙特卡罗方法计算积分时，可考虑将积分I改写为其中f*(x)为一个与f(x)有相同定义域的新的分布密度函数。于是可以这样计算积分I：这里Xi是从f*(x)中抽取的第i个子样。,偏移抽样方法,由此可以看出，原来由f(x)抽样，现改为由另一个分布密度函数f*(x)抽样，并附带一个权重纠偏因子这种方法称为偏倚抽样方法。从f(x)中抽取的Xf，满足而对于偏倚抽样，有一般情况下，Xf是具有分布f(x)总体的简单子样的个体，只代表一个。Xf*是具有分布f*(x)总体的简单子样的个体，但不代表一个，而是代表W(Xf*)个，这时Xf*是带权W(Xf*)服从分布f(x)。,在实际问题中，分布密度函数的形式有时是非常复杂的，有些甚至不能用解析形式给出，只能用数据或曲线形式给出。如中子散射角余弦分布多数是以曲线形式给出的。对于这样的分布，需要用近似分布密度函数代替原来的分布密度函数，用近似分布密度函数的抽样代替原分布密度函数的抽样，这种方法称为近似抽样方法。,近似抽样方法,设fa(x)f(x)，即fa(x)是f(x)的一个近似分布密度函数。对于阶梯近似，有其中，x0，x1，xn为任意分点。在此情况下，近似抽样方法为：或,阶梯近似,对于梯形近似，有其中，c为归一因子，fif(xi)，x0，x1，xn为任意分点。根据对称抽样方法，梯形近似抽样方法为：,梯形近似,除了上述这种近似外，近似抽样方法还包括对直接抽样方法中分布函数反函数的近似处理，以及用具有近似分布的随机变量代替原分布的随机变量。,例23.正态分布的近似抽样,我们知道，随机数的期望值为1/2，方差为1/12，则随机变量渐近正态分布，因此，当n足够大时便可用Xn作为正态分布的近似抽样。特别是n12时，有,对于任意分布密度函数f(x)，设fa(x)是f(x)的一个近似分布密度函数，它的特点是抽样简单，运算量小。令则分布密度函数f(x)可以表示为乘加分布形式：其中H1(x)为非负函数，f1(x)为一分布密度函数。对f(x)而言，fa(x)是它的近似分布密度函数，而H1(x)f1(x)正好是这种近似的修正。,近似-修正抽样方法,近似-修正抽样方法如下：抽样效率由上述近似-修正抽样方法可以看出，如果近似分布密度函数fa(x)选得好，m接近1，这时有很大可能直接从fa(x)中抽取Xfa，而只有很少的情况需要计算与f(x)有关的函数H1(Xf1)。在乘抽样方法中，每一次都要计算H(Xfa)f(Xfa)fa(Xfa)。因此，当f(x)比较复杂时，近似-修正抽样方法有很大好处。,例24.裂变中子谱分布的近似-修正抽样,裂变中子谱分布的一般形式为：其中A，B，C，Emin，Emax均为与元素有关的量。对于铀-235，A=0.965，B=2.29，C=0.453，Emin=0，Emax=。若采用乘减抽样方法，其抽样效率约为0.5。,令相应的则从fa(x)的抽样为从f1(x)的抽样为,参数的确定，使1A0，且使H1(E)的上界M1最小。裂变中子谱的近似修正抽样方法为对于铀-235，m0.8746，M0.2678，0.5543，抽样效率E0.9333。而且近似修正抽样方法有0.8746的概率直接用近似分布抽样，只计算一次对数。因此，较之乘减抽样方法大大节省了计算时间，提高了抽样效率。,为方便起见，这里仅讨论二维分布的情况，对于更高维数的分布，可用类似的方法处理。对于任意二维分布密度函数，总可以用其边缘分布密度函数和条件分布密度函数的乘积表示:其中fl(x)，f2(y|x)分别为分布f(x,y)的边缘分布密度函数和条件分布密度函数，即,多维分布抽样方法,二维分布密度函数的抽样方法是:首先由fl(x)中抽取Xf1，再由f2(y|Xf1)中抽样确定Yf2。对于多维分布密度函数，也可直接采用类似于一维分布密度函数的抽样方法。例如，对如下形式的二维分布密度函数：其中H(x,y)为非负函数，f1(x,y)为任意二维分布密度函数。设M为H(x,y)的上界，则有二维分布的乘抽样方法如下：,例25.下面二维分布密度函数的抽样,将f(x,y)写为其中用直接抽样方法分别从fl(x)和f2(y|Xf1)中抽样，得到,前面已经介绍了，指数分布的直接抽样为：这不仅需要计算对数，而且由于要使用伪随机数，受精度的限制，该抽样值在小概率处即数值较大处呈现明显得离散性。下面介绍两种抽样方法可以避免这些问题。,指数分布的抽样,所用随机数的平均个数Ne2/(e1)4.3,方法一,方法二,N,Y,作业,光子散射后能量分布的抽样把光子散射能量分布改写成如下形式进行抽样：,在1,1+2上定义如下函数:,第四章蒙特卡罗方法解粒子输运问题,屏蔽问题模型直接模拟方法简单加权法统计估计法指数变换法蒙特卡罗方法的效率作业,第四章蒙特卡罗方法解辐射屏蔽问题,辐射（光子和中子）屏蔽问题是蒙特卡罗方法最早广泛应用的领域之一。本章主要从物理直观出发，说明蒙特卡罗方法解决这类粒子输运问题的基本方法和技巧。而这些方法和技巧对于诸如辐射传播、多次散射和通量计算等一般粒子输运问题都是适用的。,屏蔽问题模型,在反应堆工程和辐射的测量与应用中，常常要用一些吸收材料做成屏蔽物挡住光子或中子。我们所关心的是经过屏蔽后射线的强度及其能量分布，这就是屏蔽问题。当屏蔽物的形状复杂，散射各向异性，材料介质不均匀,核反应截面与能量、位置有关时，难以用数值方法求解，用蒙特卡罗方法能够得到满意的结果。,粒子的输运问题带有明显的随机性质，粒子的输运过程是一个随机过程。粒子的运动规律是根据大量粒子的运动状况总结出来的，是一种统计规律。蒙特卡罗模拟，实际上就是模拟相当数量的粒子在介质中运动的状况，使粒子运动的统计规律得以重现。不过，这种模拟不是用实验方法，而是利用数值方法和技巧，即利用随机数来实现的。,为方便起见，选用平板屏蔽模型，在厚度为a，长、宽无限的平板左侧放置一个强度已知，具有已知能量、方向分布的辐射源S。求粒子穿透屏蔽概率（穿透率）及其能量、方向分布。穿透率就是由源发出的平均一个粒子穿透屏蔽的数目。同时，假定粒子在两次碰撞之间按直线运动,且粒子之间的相互作用可以忽略。,直接模拟方法,直接模拟方法就是直接从物理问题出发，模拟粒子的真实物理过程。状态参数与状态序列模拟运动过程记录结果,粒子在介质中的运动的状态，可用一组参数来描述，称之为状态参数。它通常包括：粒子的空间位置r，能量E和运动方向，以S(r,E,)表示。有时还需要其他的参数，如粒子的时间t和附带的权重W，这时状态参数为S(r,E,t,W)。状态参数通常要根据所求问题的类型和所用的方法来确定。对于无限平板几何，取S(z,E,cos)其中z为粒子的位置坐标，为粒子的运动方向与Z轴的夹角。对于球对称几何,取S(r,E,cos)其中r表示粒子所在位置到球心的距离，为粒子的运动方向与其所在位置的径向夹角。,状态参数与状态序列,粒子第m次碰撞后的状态参数为或它表示一个由源发出的粒子，在介质中经过m次碰撞后的状态，其中rm：粒子在第m次碰撞点的位置Em：粒子第m次碰撞后的能量m：粒子第m次碰撞后的运动方向tm：粒子到第m次碰撞时所经历的时间Wm：粒子第m次碰撞后的权重有时，也可选为粒子进入第m次碰撞时的状态参数。,一个由源发出的粒子在介质中运动，经过若干次碰撞后，直到其运动历史结束（如逃出系统或被吸收等）。假定粒子在两次碰撞之间按直线运动，其运动方向与能量均不改变，则粒子在介质中的运动过程可用以下碰撞点的状态序列描述：S0，S1，SM-1，SM或者更详细些,用来描述。这里S0为粒子由源出发的状态，称为初态，SM为粒子的终止状态。M称为粒子运动的链长。这样的序列称为粒子随机运动的历史，模拟一个粒子的运动过程，就变成确定状态序列的问题。,为简单起见，这里以中子穿透均匀平板的模型来说明，这时状态参数取S(z,E,cos)。模拟的步骤如下：(1)确定初始状态S0：确定粒子的初始状态，实际上就是要从中子源的空间位置、能量和方向分布中抽样。设源分布为则分别从各自的分布中抽样确定初始状态。对于平板情况，抽样得到z00。,模拟运动过程,(2)确定下一个碰撞点：已知状态Sm-1，要确定状态Sm，首先要确定下一个碰撞点的位置zm。在相邻两次碰撞之间，中子的输运长度l服从如下分布：对于平板模型，l服从分布：其中，t为介质的中子宏观总截面，积分称为粒子输运的自由程数，系统的大小通常就是用系统的自由程数表示的。,显然，粒子输运的自由程数服从指数分布，因此从f(l)中抽样确定l，就是要从积分方程中解出l。对于单一介质则下一个碰撞点的位置如果zma，则中子穿透屏蔽，若zm0,则中子被反射出屏蔽。这两种情况，均视为中子历史终止。,(3)确定被碰撞的原子核：通常介质由几种原子核组成，中子与核碰撞时，要确定与哪一种核碰撞。设介质由A、B、C三种原子核组成，其核密度分别为NA、NB、NC，则介质的宏观总截面为：其中分别为核A、B、C的宏观总截面。其定义如下：分别表示()核的宏观总截面、核密度和微观总截面。,由于中子截面表示中子与核碰撞可能性的大小，因此，很自然地，中子与A、B、C核发生碰撞的几率分别为：利用离散型随机变量的抽样方法，确定碰撞核种类：,(4)确定碰撞类型：确定了碰撞的核(比如B核)后，就要进一步确定碰撞类型。中子与核的反应类型有弹性散射、非弹性散射、(n,2n)反应，裂变和俘获等，它们的微观截面分别为则有各种反应发生的几率分别为,利用离散型随机变量的抽样方法，确定反应类型。在屏蔽问题中，中子与核反应常只有弹性散射和吸收两种类型，吸收截面为：这时，总截面为：发生弹性散射的几率为：若，则为弹性散射；否则为吸收，发生吸收反应意味着中子的历史终止。,(5)确定碰撞后的能量与运动方向：如果中子被碰撞核吸收，则其输运历史结束。如果发生弹性散射，需要确定散射后中子的能量和运动方向。中子能量Em为：A是碰撞