第一章 蒙特卡罗方法概述
蒙特卡罗方法的基本思想
蒙特卡罗方法的收敛性，误差
蒙特卡罗方法的特点
蒙特卡罗方法的主要应用范围
作 业
第一章 蒙特卡罗方法概述
蒙特卡罗方法又称随机抽样技巧或统计试验方法。半个多世纪以来，由于科学技术的发展和电子计算机的发明 ，这种方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。
蒙特卡罗方法的基本思想
二十世纪四十年代中期，由于科学技术的发展和电子计算机的发明，蒙特卡罗方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。但其基本思想并非新颖，人们在生产实践和科学试验中就已发现，并加以利用。
两个例子
例1. 蒲丰氏问题
例2. 射击问题（打靶游戏）
基本思想
计算机模拟试验过程
例1. 蒲丰氏问题
为了求得圆周率π值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（ l＜a）的平行线相交的频率代替概率P，再利用准确的关系式：
求出π值
其中Ｎ为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。
一些人进行了实验，其结果列于下表 ：
实验者
年份
投计次数
π的实验值
沃尔弗(Wolf)
1850
5000
斯密思(Smith)
1855
3204
福克斯(Fox)
1894
1120
拉查里尼(Lazzarini)
1901
3408
例2. 射击问题（打靶游戏）
设r表示射击运动员的弹着点到靶心的距离，ｇ(r)表示击中r处相应的得分数（环数），f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射击水平。该运动员的射击成绩为
用概率语言来说，<g>是随机变量ｇ(r)的数学期望，即
现假设该运动员进行了Ｎ次射击，每次射击的弹着点依次为r1，r2，…，rN，则Ｎ次得分g(r1)，g(r2)，…，g(rN)的算术平均值
代表了该运动员的成绩。换言之，为积分<g>的估计值，或近似值。
在该例中，用Ｎ次试验所得成绩的算术平均值作为数学期望<g>的估计值（积分近似值）。
基本思想
由以上两个例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。
当随机变量的取值仅为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)，作Ｎ次试验后，得到该运动员射击成绩的近似值
概率
10
9
8
7
环数
蒙特卡罗方法的收敛性，误差
蒙特卡罗方法作为一种计算方法，其收敛性与误差是普遍关心的一个重要问题。
收敛性
误差
减小方差的各种技巧
效率
收敛性
由前面介绍可知，蒙特卡罗方法是由随机变量X的简单子样X1，X2，…，XN的算术平均值：
作为所求解的近似值。由大数定律可知，
如X1，X2，…，XN独立同分布，且具有有限期望值（E(X)<∞），则
即随机变量X的简单子样的算术平均值 ，当子样数Ｎ充分大时，以概率1收敛于它的期望值E(X)。
误差
蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，…，XN独立同分布，且具有有限非零的方差σ2 ，即
f(X)是X的分布密度函数。则
当N充分大时，有如下的近似式
其中α称为置信度，1－α称为置信水平。
这表明，不等式 近似地以概率
1－α成立，且误差收敛速度的阶为 。
通常，蒙特卡罗方法的误差ε定义为
上式中 与置信度α是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出 。
下面给出几个常用的α与的数值：
关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。第二，误差中的均方差σ是未知的，必须使用其估计值
来代替，在计算所求量的同时，可计算出 。
3
α
减小方差的各种技巧
显然，当给定置信度α后，误差ε由σ和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十个数字组成，每个数字以的等概率出现，数字之间相互独立。这些数字序列叫作随机数字序列。如果要得到n位有效数字的随机数，只需将表中每n个相邻的随机数字合并在一起，且在最高位的前边加上小数点即可。例如，某随机数表的第一行数字为7634258910…，要想得到三位有效数字的随机数依次为，，。
因为随机数表需在计算机中占有很大内存，而且也难以满足蒙特卡罗方法对随机数需要量非常大的要求，因此，该方法不适于在计算机上使用。
物理方法
用物理方法产生随机数的基本原理是：利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。用来作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。
一般情况下，任意一个随机数在计算机内总是用二进制的数表示的：
其中εi（i=1，2，…，m）或者为0，或者为1。
因此，利用物理方法在计算机上产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为。
用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，需要增加随机数发生器和电路联系等附加设备，费用昂贵。因此，该方法也不适合在计算机上使用。
伪随机数
伪随机数
伪随机数存在的两个问题
伪随机数的周期和最大容量
伪随机数
在计算机上产生随机数最实用、最常见的方法是数学方法，即用如下递推公式：
产生随机数序列。对于给定的初始值ξ1,ξ2…，ξk，确定ξn+k，ｎ=1，2，…。经常使用的是k=1的情况，其递推公式为：
对于给定的初始值ξ1，确定ξn+1，ｎ=１,２…
伪随机数存在的两个问题
用数学方法产生的随机数，存在两个问题：
递推公式和初始值ξ1,ξ2…，ξk确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。
由于随机数序列是由递推公式确定的，而在计算机上所能表示的[0，1]上的数又是有限的，因此，这种方法产生的随机数序列就不可能不出现无限重复。一旦出现这样的n＇，n″ (n＇< n″ )，使得下面等式成立：
随机数序列便出现了周期性的循环现象。对于k=1的情况，只要有一个随机数重复，其后面的随机数全部重复，这与随机数的要求是不相符的。
由于这两个问题的存在，常称用数学方法产生的随机数为伪随机数。对于以上存在的两个问题，作如下具体分析。
关于第一个问题，不能从本质上加以改变，但只要递推公式选得比较好，随机数间的相互独立性是可以近似满足的。至于第二个问题，则不是本质的。因为用蒙特卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所用随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。
用数学方法产生的伪随机数容易在计算机上得到，可以进行复算，而且不受计算机型号的限制。因此，这种方法虽然存在着一些问题，但仍然被广泛地在计算机上使用，是在计算机上产生伪随机数的主要方法。
伪随机数的周期和最大容量
发生周期性循环现象的伪随机数的个数称为伪随机数的周期。对于前面介绍的情况，伪随机数的周期为n″－n＇。
从伪随机数序列的初始值开始，到出现循环现象为止，所产生的伪随机数的个数称为伪随机数的最大容量。前面的例子中，伪随机数的最大容量为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 (b≥2)
c = 1
这样在计算中可以使用移位和指令加法，提高计算速度。
产生伪随机数的其他方法
取中方法
加同余方法
伪随机数序列的均匀性和独立性
判断伪随机数序列是否满足均匀和相互独立的要求，要靠统计检验的方法实现。对于伪随机数的统计检验，一般包括两大类：均匀性检验和独立性检验。
六十年代初，人们开始用定性的方法研究伪随机数序列的均匀性和独立性问题，简要叙述如下。
伪随机数的均匀性
这里只考虑伪随机数序列ξ1,ξ2…，ξn全体作为子样时的均匀性问题。其中n为伪随机数序列的最大容量。
对于任意的0≤x≤1，令Nn(x)表示伪随机数序列ξ1,ξ2…，ξn中适合不等式
ξi< x i=1，2，…，n
的个数，则
标志伪随机数序列ξ1,ξ2…，ξn的均匀程度，称为均匀偏度。
将伪随机数序列ξ1,ξ2…，ξn从小至大重新排列
并令 ，则由δ(n)的定义，容易证明
很明显，对于固定的ｎ，δ(n)的值越小越好。它是描述伪随机数序列均匀程度的基本量。对于任意随机数序列，均有如下不等式成立：
当 时，所对应的伪随机数序列为最佳分布。
可以证明，伪随机数序列为最佳分布的充要条件是它取遍序列
的所有值。
对于计算机上使用的乘同余方法，按照前面介绍的方法选取a、x1时，所产生的伪随机数序列的均匀偏度
对于乘加同余方法
对于部分伪随机数的均匀性问题通常用统计检验方法检验。
伪随机数的独立性
对于任意 ，令 表示(ξ1,ξ2)， (ξ2,ξ3)，…， (ξn,ξn+1)中适合不等式
的个数，根据随机变量间相互独立的定义和频率近似概率的方法，令
则ε(n)标志伪随机数序列ξ1,ξ2…，ξn的独立程度，简称为独立偏度。对于固定的n,ε(n)的值越接近于零，伪随机数序列的独立性越好。
对于乘同余方法，
对于乘加同余方法，
因此，这两种方法的独立性都是很好的。
同伪随机数的均匀性问题一样，伪随机数序列的独立性问题也是对它的全体讨论的。若只考虑伪随机数的一部分，在通常情况下给出ε(i)是相当因难的。因此，伪随机数序列的独立性问题的统计检验方法同样是非常重要的。
作 业
证明1—ξ是随机数。
证明 与 同分布 。
第三章 由已知分布的随机抽样
随机抽样及其特点
直接抽样方法
挑选抽样方法
复合抽样方法
复合挑选抽样方法
替换抽样方法
随机抽样的一般方法
随机抽样的其它方法
作 业
第三章 由已知分布的随机抽样
本章叙述由己知分布抽样的各主要方法，并给出在粒子输运问题中经常用到的具体实例。
随机抽样及其特点
由巳知分布的随机抽样指的是由己知分布的总体中抽取简单子样。随机数序列是由单位均匀分布的总体中抽取的简单子样，属于一种特殊的由已知分布的随机抽样问题。本章所叙述的由任意已知分布中抽取简单子样，是在假设随机数为已知量的前提下，使用严格的数学方法产生的。
为方便起见，用XF表示由己知分布F(x)中产生的简单子样的个体。对于连续型分布，常用分布密度函数f(x)表示总体的己知分布，用Xf表示由己知分布密度函数f(x)产生的简单子样的个体。另外，在抽样过程中用到的伪随机数均称随机数。
直接抽样方法
对于任意给定的分布函数F(x)，直接抽样方法如下：
其中，ξ1，ξ2，…，ξN为随机数序列。为方便起见，将上式简化为：
若不加特殊说明，今后将总用这种类似的简化形式表示，ξ总表示随机数。
证明
下面证明用前面介绍的方法所确定的随机变量序列X1，X2，…，XN具有相同分布F(x)。
对于任意的n成立，因此随机变量序列X1，X2，…，XN具有相同分布F(x)。另外，由于随机数序列ξ1，ξ2，…，ξN是相互独立的，而直接抽样公式所确定的函数是波雷尔（Borel）可测的，因此，由它所确定的X1，X2，…，XN也是相互独立的（[, Measure theory, Nosrtand,1950]§45定理2）。
离散型分布的直接抽样方法
对于任意离散型分布：
其中x1，x2，…为离散型分布函数的跳跃点，P1，P2，…为相应的概率，根据前述直接抽样法，有离散型分布的直接抽样方法如下：
该结果表明，为了实现由任意离散型分布的随机抽样，直接抽样方法是非常理想的。
例1. 二项分布的抽样
二项分布为离散型分布，其概率函数为：
其中，P为概率。对该分布的直接抽样方法如下：
例2. 泊松(Possion)分布的抽样
泊松(Possion)分布为离散型分布，其概率函数为：
其中，λ>0 。对该分布的直接抽样方法如下：
例3. 掷骰子点数的抽样
掷骰子点数X=n的概率为：
选取随机数ξ，如
则
在等概率的情况下，可使用如下更简单的方法：
其中［］表示取整数。
例4. 碰撞核种类的确定
中子或光子在介质中发生碰撞时，如介质是由多种元素组成，需要确定碰撞核的种类。假定介质中每种核的宏观总截面分别为Σ1，Σ2，…，Σn，则中子或光子与每种核碰撞的概率分别为：
其中Σt＝Σ1＋Σ2＋…＋Σn。碰撞核种类的确定方法为：产生一个随机数ξ，如果
则中子或光子与第I种核发生碰撞。
例5. 中子与核的反应类型的确定
假设中子与核的反应类型有如下几种:弹性散射，非弹性散射，裂变，吸收，相应的反应截面分别为Σel，Σin，Σf，Σa。则发生每一种反应类型的概率依次为 ：
其中反应总截面Σt＝Σel＋Σin＋Σf＋Σa。
反应类型的确定方法为：产生一个随机数ξ
连续型分布的直接抽样方法
对于连续型分布，如果分布函数F(x) 的反函数 F－1(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的分布密度函数为
分布函数为
容易知道，该分布的直接抽样方法是
由于开方运算在计算机上很费时间，该方法不是好方法。下面使用挑选抽样方法：取
则抽样框图为
＞
≤
显然，没有必要舍弃ξ1＞ξ2的情况，此时，只需取 就可以了，亦即
另一方面，也可证明 与 具有相同的分布 。
复合抽样方法
在实际问题中，经常有这样的随机变量，它服从的分布与一个参数有关，而该参数也是一个服从确定分布的随机变量，称这样的随机变量服从复合分布。例如，分布密度函数
是一个复合分布。其中Pn≥0，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
复合挑选抽样方法
考虑另一种形式的复合分布如下：
其中0≤H(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,π］上均匀分布，作变换
其中0≤ρ≤1，0≤ρ≤π，则
(x,y) 表示上半个单位圆内的点。如果 (x,y) 在上半个单位圆内均匀分布，则θ在［0,π］上均匀分布，由于
因此抽样sinφ和cosφ的问题就变成在上半个单位圆内均匀抽样 (x,y) 的问题。
为获得上半个单位圆内
的均匀点，采用挑选法，在
上半个单位圆的外切矩形内
均匀投点（如图）。
舍弃圆外的点，余下的就是所要求的点。
抽样方法为：
抽样效率
E=π/4≈
>
为实现散射方位角余弦分布抽样，最重要的是在上半个单位圆内产生均匀分布点。下面这种方法，首先在单位圆的半个外切正六边形内产生均匀分布点，如图所示。
于是便有了抽样效率更高的抽样方法：
抽样效率
>
≤
例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个，即
随机抽样的一般方法
加抽样方法
减抽样方法
乘抽样方法
乘加抽样方法
乘减抽样方法
对称抽样方法
积分抽样方法
加抽样方法
加抽样方法是对如下加分布给出的一种抽样方法：
其中Pn≥0， ，且 fn(x)为与参数n有关的分布密度函数，n=1，2，…。
由复合分布抽样方法可知，加分布的抽样方法为:首先抽样确定n’，然后由 fn’(x)中抽样x，即：
例14. 多项式分布抽样
多项式分布密度函数的一般形式为：
将 f(x) 改写成如下形式：
则该分布的抽样方法为：
例15. 球壳内均匀分布抽样
设球壳内半径为R0，外半径为R1，点到球心的距离为r，则r的分布密度函数为
分布函数为
该分布的直接抽样方法是
为避免开立方根运算，作变换：
则 x∈[0,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. β分布抽样
β分布的一个特例：
取A1＝2，A2＝1，f1(x)＝1，f2(x)＝2x，此时m＝0，则根据第一种形式的减抽样方法，有
或
＞
≤
＞
≤
由于1－ξ1可用ξ1代替，该抽样方法可简化为：
对于ξ2＞ξ1的情况，可取 Xf＝ξ1 ，因此
与β分布的推论相同。
＞
≤
乘抽样方法
如下形式的分布称为乘分布：
其中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=1－P1），这对有些分布是很困难的。下面的方法可以不用计算P1 ：
对于任意小于1的正数P1 ，令P2=1－P1 ；
则采用复合挑选抽样方法，有：
当取
时，抽样效率最高
这时，乘加抽样方法为：
>
>
>
≤
由于
可知第一种方法比第二种方法的抽样效率高。
例19. 光子散射后能量分布的抽样
令光子散射前后的能量分别为 和 （以 m0c2 为单位，m0为电子静止质量，c 为光速）， ，
则 x 的分布密度函数为：
该分布即为光子散射能量分布，它是由著名的Klin－Nishina 公式确定的。其中 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 达到最小，此时
取 m＝0，令
则裂变中子谱分布的抽样方法为:
抽样效率
＞
≤
对称抽样方法
对称分布的一般形式为：
其中 f1(x) 为任意分布密度函数，满足偶函数对称条件，H(x) 为任意奇函数，即对任意x满足：
对称分布的抽样方法如下：取η=2ξ－1
>
≤
证明：
因为η=2ξ－1，η≤x 相当于ξ≤ ，因此
例21. 质心系各向同性散射角余弦分布抽样
在质心系各向同性散射的假设下，为得到实验室系散射角余弦，需首先抽样确定质心条散射角余弦：
再利用下面转换公式:
得到实验室系散射角余弦μL。其中A为碰撞核质量，θC、θL 分别为质心系和实验室系散射角。
为避免开方运算，可以使用对称分布抽样。
根据转换公式可得：
依照质心系散射各向同性的假定，可得到实验室系散射角余弦μL 的分布如下：
该密度函数中的第一项为偶函数，第二项为奇函数，因而是对称分布。其中
从 f1(μL) 的抽样可使用挑选法
然后再以
的概率决定接受或取负值。
上述公式涉及开方运算，需要进一步简化。
＞
≤
注意以下事实：对于任意0≤a≤1
令
则上述挑选抽样中的挑选条件简化为：
另一方面，在 即 的条件下，η2/a 在［－1, 1］上均匀分布，故可令η＝η2/a，则最终决定取正负值的条件简化为：
于是，得到质心系各向同性散射角余弦分布的抽样方法为：
＞
≤
＞
≤
积分抽样方法
如下形式的分布密度函数
称为积分分布密度函数，其中 f0(x，y) 为任意二维分布密度函数，H(x)为任意函数。该分布密度函数的抽样方法为：
>
证明：对于任意x
例22. 各向同性散射方向的抽样
为了确定各向同性散射方向 ，根据公式：
对于各向同性散射，cosθ在［－1, 1］上均匀分布，φ在［0, 2π］上均匀分布。由于
直接抽样需要计算三角函数和开方。
定义两个随机变量：
可以证明，当 时，随机变量 x 和 y 服从如下分布:
定义区域为：
则 w＝cosθ 的分布可以用上述分布表示成积分分布的形式：
令 ，则属于上述积分限内的 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 为归一因子， fi ＝ f (xi) ，x0，x1，… ，xn为任意分点。根据对称抽样方法，梯形近似抽样方法为：
＞
≤
除了上述这种近似外，近似抽样方法还包括对直接抽样方法中分布函数反函数的近似处理，以及用具有近似分布的随机变量代替原分布的随机变量。
例23. 正态分布的近似抽样
我们知道，随机数ξ的期望值为 1/2，方差为 1/12，则随机变量
渐近正态分布，因此，当 n 足够大时便可用 Xn 作为正态分布的近似抽样。特别是 n＝12 时，有
近似-修正抽样方法
对于任意分布密度函数 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=，B=，C=，Emin=0，Emax=∞。
若采用乘减抽样方法，其抽样效率约为。
令
相应的
则
从 fa(x) 的抽样为
从 f1(x) 的抽样为
参数λ的确定，使1－Aλ＞0，且使 H1(E) 的上界M1 最小。裂变中子谱的近似修正抽样方法为
对于铀-235，m≈，M≈，λ≈，抽样效率 E≈。而且近似修正抽样方法有的概率直接用近似分布抽样，只计算一次对数。因此，较之乘减抽样方法大大节省了计算时间，提高了抽样效率。
≤
≤
＞
＞
多维分布抽样方法
为方便起见，这里仅讨论二维分布的情况，对于更高维数的分布，可用类似的方法处理。
对于任意二维分布密度函数，总可以用其边缘分布密度函数和条件分布密度函数的乘积表示:
其中 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) 中抽样，得到
指数分布的抽样
前面已经介绍了，指数分布
的直接抽样为：
这不仅需要计算对数，而且由于要使用伪随机数，受精度的限制，该抽样值在小概率处即数值较大处呈现明显得离散性。
下面介绍两种抽样方法可以避免这些问题。
方法一
所用随机数的平均个数 N＝e2 / ( e－1)≈
＞
≤
N
Y
方法二
＞
≤
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 ：
确定粒子的初始状态，实际上就是要从中子源的空间位置、能量和方向分布中抽样。设源分布为
则分别从各自的分布中抽样确定初始状态。
对于平板情况，
抽样得到 z0＝0。
(2) 确定下一个碰撞点 ：
已知状态Sm-1，要确定状态Sm，首先要确定下一个碰撞点的位置 zm。在相邻两次碰撞之间，中子的输运长度 l 服从如下分布：
对于平板模型，l 服从分布：
其中，Σt 为介质的中子宏观总截面，
积分 称为粒子输运的自由程数，
系统的大小通常就是用系统的自由程数表示的。
显然，粒子输运的自由程数服从指数分布，
因此从 f ( l ) 中抽样确定 l，就是要从积分方程
中解出 l。
对于单一介质
则下一个碰撞点的位置
如果 zm≥a，则中子穿透屏蔽，若 zm≤0, 则中子被反射出屏蔽。这两种情况，均视为中子历史终止。
(3) 确定被碰撞的原子核 ：
通常介质由几种原子核组成，中子与核碰撞时，要确定与哪一种核碰撞。设介质由A、B、C 三种原子核组成，其核密度分别为NA、NB、NC，则介质的宏观总截面为：
其中 分别为核A、B、C 的宏观总截面。其定义如下：
分别表示(·)核的宏观总截面、核密度和微观总截面。
由于中子截面表示中子与核碰撞可能性的大小，因此，很自然地，中子与A、B、C 核发生碰撞的几率分别为：
利用离散型随机变量的抽样方法，确定碰撞核种类：
＞
≤
＞
≤
(4) 确定碰撞类型 ：
确定了碰撞的核(比如B核)后，就要进一步确定碰撞类型。中子与核的反应类型有弹性散射、非弹性散射、(n,2n)反应，裂变和俘获等，它们的微观截面分别为
则有
各种反应发生的几率分别为
利用离散型随机变量的抽样方法，确定反应类型。
在屏蔽问题中，中子与核反应常只有弹性散射和吸收两种类型，吸收截面为：
这时，总截面为：
发生弹性散射的几率为：
若 ，则为弹性散射；否则为吸收，发生吸收反应意味着中子的历史终止。
(5) 确定碰撞后的能量与运动方向：
如果中子被碰撞核吸收，则其输运历史结束。如果发生弹性散射，需要确定散射后中子的能量和运动方向。中子能量 Em 为：
A是碰撞核的质量与中子质量之比，一般就取元素的原子量；θC 为质心系中中子散射前后方向间的夹角，即偏转角。 可从质心系中弹性散射角分布
fC(μC) 中抽样产生。实验室系散射角θL的余弦μL为：
如果给出实验室系散射角余弦分布 fL(μL)，可直接从 fL(μL)中抽取μL，此时能量Em与μL的关系式为：
确定了实验室系散射角θL后，
再使用球面三角公式
确定cosαm ：
其中χ为在［0,2π］上均匀分布的方位角。
至此，由Sm-1完全可以确定Sm。 因此，当中子由源出发后，即S0确定后，重复步骤 (2)～(5)，直到中子游动历史终止。于是得到了一个中子的随机游动历史 S0 ，S1 ，…，SM-1 ，SM，即
也就是模拟了一个由源发出的中子的运动过程。
以上模拟过程可分为两大步：第一步确定粒子的初始状态S0，第二步由状态Sm-1来确定状态Sm。这第二步又分为两个过程：第一个过程是确定碰撞点位置zm ，称为输运过程；第二个过程是确定碰撞后粒子的能量及运动方向，称为碰撞过程。对于中子而言，碰撞过程是先确定散射角，进而确定能量和运动方向；而对于光子，碰撞过程是先确定能量，再确定散射角以及运动方向。重复这两个过程，直至粒子的历史终止。
这种模拟过程，是解任何类型的粒子输运问题所共有的，它是蒙特卡罗方法解题的基本手段。
记录结果
在获得中子的随机游动历史后，我们要对所要计算的物理量进行估计。对于屏蔽问题，我们要计算中子的穿透率。考察每个中子的随机游动历史，它可能穿透屏蔽（zM≥a），可能被屏蔽发射回来（zM≤0），或者被吸收。设第 n 个中子对穿透的贡献为ηn ，则
如果我们共跟踪了N 个中子，则穿透屏蔽的中子数为：
则穿透屏蔽概率的近似值为：
它是穿透率的一个无偏估计。
我们称这种直观地模拟过程和估计方法为直接模
拟方法。在置信水平 1－α＝ 时， 的误差为：
其中 为ηn的均方差，由于ηn是一个服从二项分布的随机变量，所以
或
为得到中子穿透屏蔽的能量、角分布，将能量、角度范围分成若干个间隔：
其中Emax，Emin分别表示能量的上、下限，对于穿透屏蔽的中子按其能量、方向分间隔记录。设一穿透屏蔽的中子能量为EM，其运动方向与Z轴夹角为αM，若能量EM属于第 i 个能量间隔ΔEi，角度αM属于第 j 个角度间隔Δαj，则分别在第 i 个能量计数器及第 j 个角度计数器中加 "1"。
跟踪 N 个中子后，则
分别为穿透中子的能量分布和角分布。其中N1,i 和 N2,i 分别为第 i 个能量和第 j 个角度间隔的穿透中 子数。归一后分别为：
简单加权法
从模拟物理过程来说，直接模拟法是最简单、也是最基本的方法。但是，在直接模拟法中，不管中子在屏蔽中经过多少次碰撞，只要在介质中被吸收，对穿透的贡献就为零；因此在所跟踪的粒子中绝大部分都对穿透没有贡献。而在许多屏蔽问题中，穿透率的数量级在10-6到10-8。进一步，如果我们要求穿透率达相对误差小于1％，即
那么，N 要大到惊人的数量级1010到1012。显然，这时用直接模拟法计算不是很有效。
简单加权法
屏蔽物一般是由吸收强的介质组成，因此在每次碰撞时，粒子很有可能被吸收而停止跟踪。现在改变
模拟方法，在判断碰撞类型 时，可以认为粒子
的 部分是弹性散射，而其余部
分被吸收，即人为地把中子分成两部分，一部分弹性散射，一部分吸收。弹性散射这部分继续跟踪；吸收部分则停止跟踪。也就是说，我们利用中子权重的变化来反应继续弹性散射的部分。这就是简单加权法的基本思想。
显然，在加权法中中子的权重W 已成为中子状态参数的组成部分。这时，中子历史成为：
对源中子，取W0=1。经过碰撞中子权重的变化为：
因子 称为尚存因子。
这时，第 n 个中子对穿透的贡献为：
如果我们共跟踪了N个中子，则穿透率P的无偏估计为：
类似地，可以得到穿透中子的能量分布和角分布。只不过在对各计数器进行的加 "1" 操作改为加WM。
简单加权法的方差
简单加权法的方差估计为：
与直接模拟法相比，有
注意到ζn≤1，有
这表明简单加权法的方差小于直接模拟法的方差。这是因为加权法比直接模拟法减少了一次随机抽样。
权重方法的其它应用
加权法的思想在蒙特卡罗方法中用途很广泛。例如，对于具有中子增殖反应，如裂变，(n,2n)，(n,3n) 反应的中子输运问题，一个中子与核发生碰撞后，根据反应的类型会产生不同数量的次级中子，每个次级中子又会产生新的次级中子，这样链锁反应 下去，使得用直接模拟法模拟每一个中子是非常困难的。这种情况可以利用加权法来处理。
中子与核发生碰撞 后，产生的次级中子平均数为：
这里νf 为裂变次级中子数。于是，碰撞后的权重为：
而决定碰撞类型的几率分别为：
其中
加权法的思想，还可以应用到连续分布情况和偏倚抽样的问题
统计估计法
加权法虽然改进了直接模拟法，但它同样只关心中子是否穿透屏蔽这一信息，因此对每一个中子历史的信息利用得很不充分。统计估计法能够较多地利用中子的历史信息，因而能得到更好的结果。
一个中子，可能在介质内不发生碰撞而直接穿透屏蔽，也可能在介质内发生一次碰撞后再穿透屏蔽，或经过二次碰撞穿
透屏蔽，等等，
这些事件是互不相
容的，因此穿透概
率P 可表示为：
其中Pm 是中子恰好经过 m 次碰撞而穿透屏蔽的概率。这表明，可以用求 Pm (m=0,1, … ) 的方法得到P。这样，中子对穿透概率的贡献就不只限于末次碰撞了。
α0
α1
S0
S1
Sm
αm
P1
P0
Pm
Z
a
0
设中子的历史为：
根据该中子的历史，我们可以估计出中子恰好经
过 m 次碰撞后，穿透屏蔽的部分
显然，具有初态 S0＝( 0, E0, cosα0,W0 ) 的中子，未经碰撞直接穿透的部分是：
类似地，在经过了第 m 次碰撞后的中子具有状态 Sm＝( zm, Em, cosαm,Wm ) ，其可能穿透的部分，正好是一个中子恰好经过 m 次碰撞穿透的部分：
这里的这种估计技巧，由于是对每次碰撞后的状态，求其后未经碰撞直接穿透的贡献，因此该方法也称为最后自由飞行估计。
于是得到该中子对穿透的贡献：
如果我们共跟踪了N个中子，则穿透率P的估计为：
其方差估计为：
指数变换法
在直接模拟方法中，相对误差为
其中 为与置信水平 1－α相应的量。
如果构造一个新的概率模型，使得该模型的穿透率P*与原模型的穿透率P之间存在关系：
使用直接模拟方法 , 相对误差为
如果令ε*＝ε，即
这意味着，达到同样的相对误差，跟踪粒子的数目缩小 K 倍，从而减少 K 倍的计算量。指数变换法就是构造一个新的概率模型的一个有效方法。
构造如下伪过程：宏观总截面为
散射截面仍为Σel(E)。其中 Emin、Emax 分别为能量的下限和上限，α为粒子的运动方向与 Z 轴的夹角。可以证明这个伪过程的穿透概率P* 与原过程的穿透概率P之间有如下关系 :
显然， 。因伪过程与原过程的结果相差 e 指数，所以该方法称为指数变换法。
分析一下伪过程的定义 , 可以明显看出P*增大的原因。当 cosα=1 时，粒子运动方向与 Z 轴方向一致，其截面最小，粒子沿 Z 轴方向输运的距离较远；而当 cosα=－1 时，粒子运动方向背向 Z 轴方向，这时其截面最大，粒子向后输运的距离较短。因此，截面变换的结果是加强了粒子向前运动的能力，因而使穿透概率增大。
伪过程的构造与几何形状及所考虑的问题有关。比如，对球形几何，使用指数变换法求穿透概率时 , 所构造的宏观总截面与平板屏蔽的情况不同，粒子的模拟方法也较复杂。
蒙特卡罗方法的效率
衡量一种蒙特卡罗技巧的好坏，除了看其方差大小外，还要看其所需费用（计算时间）多少，即从该技巧的效率 Ef（方差与费用乘积的倒数）全面考虑：
其中σ2 为方差，T 为所需费用。Ef 大时，所用方法的效率高；否则，效率低。
在一般情况下，有些方法虽然减小了方差，却增加了费用。例如，加权法、统计估计法虽然较直接模拟方法减小了方差，却使每个粒子的运动链长增加，或记录贡献的计算时间增加。因 此，不能认为方差小的方法一定好，要从方法的效率全面考虑。在有些情况下，直接模拟方法仍然是一个被广泛使用的方法。
作 业
光子散射后能量分布的抽样
把光子散射能量分布改写成如下形式进行抽样：
在［1, 1+2α］上定义如下函数:
作 业
给出分布密度函数
的抽样方法。
＞
≤
第五章 蒙特卡罗方法在计算机上的实现
源分布抽样过程
空间、能量和运动方向的随机游动过程
记录贡献和分析结果过程
核截面数据的引用
蒙特卡罗程序结构
作 业
第五章 蒙特卡罗方法在计算机上的实现
蒙特卡罗方法是随着计算机的出现和发展而逐步发展起来的。在计算机上能够产生符合要求的随机数，实现对已知分布的抽样，奠定了蒙特卡罗方法在计算机上得以实现的基础。在计算机上使用蒙特卡罗方法解粒子输运问题大致包括三个过程：源分布抽样过程，空间、能量和运动方向的随机游动过程以及记录、分析结果过程 。
源分布抽样过程
源分布抽样的目的是产生粒子的初始状态
。下面我们介绍一些常见的特定
类型的源分布抽样方法。
源粒子的位置常见分布的随机抽样
圆内均匀分布
设圆半径为R0，粒子在圆内均匀分布时，从发射点到中心的距离 r 的分布密度函数为：
r 的抽样方法为：
圆环内均匀分布
设圆环的内半径为R0，外半径为R1，则粒子在该圆环内均匀分布时，从发射点到中心的距离 r 的分布密度函数为：
r 的抽样方法为：
≤
>
球内均匀分布
设球的半径为R，粒子在球内均匀分布时，从发射点到中心的距离 r 的分布密度函数为：
r 的抽样方法为：
在直角坐标系下，抽样方法为：
≤
>
球壳内均匀分布
设球壳的内半径为R0，外半径为R1，在均匀分布时，从发射点到中心的距离 r 的分布密度函数为：
r 的抽样方法为：
≤
≤
>
>
在直角坐标系下，球壳内点的坐标为：
其中，r 由前面的抽样方法确定，θ、φ服从各向同性分布，其抽样方法为：
＞
≤
圆柱内均匀分布
圆柱内均匀分布是指粒子发射点均匀地分布在底半径为 R，高为 2H 的圆柱内。若固定圆柱的中心为原点，圆柱的轴向为 z 轴，则分布密度函数为：
抽样方法为：
≤
>
点源分布
点源分布是指粒子由一固定点 发射，其分布密度函数为：
其中， 为狄拉克δ函数，源粒子的抽样方法为：
在球坐标系中，粒子发射点到球心的距离 r 的分布密度函数为：
其中， 为点源到球心的距离。源粒子的位置抽样为：
球外平行束源分布
球外平行束源分布是指粒子平行入射到半径为 R 的球面上，或球外点源距离球很远，可以近似地看作平行束源。设 r 为粒子发射点到球心的距离 , 其分布密度函数为：
r 的抽样方法为：
在直角坐标系中，抽样方法为：
≤
>
源粒子的能量常见分布的随机抽样
单能源分布
单能源分布是指粒子的发射能量为一固定值 E0 ，其分布密度函数为 ：
源粒子的能量为：
裂变中子谱分布
裂变中子谱分布的一般形式为：
其中A，B，C，Emin，Emax 均为与元素有关的量。
对于铀-235，
A=，B=，C=，Emin=0，Emax=∞。
采用近似修正抽样，抽样方法为：
其中，m≈，M1≈，λ≈。
此外，裂变谱分布也有以数值曲线形式给出的，此时，用数值曲线抽样方法抽取 E 。
≤
≤
＞
＞
麦克斯韦(Maxwell) 谱分布
麦克斯韦谱分布的一般形式为：
该分布的抽样方法为
＞
≤
源粒子运动方向常见分布的随机抽样
各向同性分布
各向同性分布密度函数为：
其中，μ＝cosθ，θ为运动方向与 z 轴的夹角，φ为方位角。
在直角坐标系下，各方向余弦 u，v，w 为：
其抽样方法为：
＞
≤
半面各向同性分布
不妨设在 x≥0 的半面方向上各向同性发射粒子，则在前述各向同性分布的抽样方法中，用ξ2代替η2就能得到所需分布的抽样。对于其它方向的情况，可用类似的方法处理。
球外平行束源分布
令μ＝cosθ，θ为粒子运动方向的径向夹角，则μ分布密度函数为：
μ的抽样方法为：
球外各向同性点源分布
设球外点源 S 到球心的距离为D0。点源 S 到球的最大张角为θ*，
则球外各向同性点源分布的抽样方法是：
先抽样确定 ，再转换成θ。
在直角坐标系下，取 OS 为 z 轴，抽样方法为：
次级粒子的源分布
在有关次级粒子（如裂变中子，中子生成光子，光子生成中子）的输运过程中，次级粒子源分布的抽样方法，主要可分为以下两种：
直接生成法
可将生成的次级粒子的位置、能量、方向、权重等参数直接作为源分布的抽样结果。也就是直接对生成的次级粒子进行跟踪。这种方法比较简单、直观。
离散分布法
将生成的次级粒子的权重，按空间位置、能量、方向分别记录，得到次级粒子的空间、能量、运动方向的离散的近似分布。再根据该分布，利用各种抽样技巧，得到源分布的抽样，对抽样的源粒子进行跟踪、记录。
当一个问题需要用两个以上的蒙卡程序处理时，可采用这种方法。
空间、能量和运动方向的随机游动过程
粒子由状态Sm到状态Sm+1时，需要确定粒子的空间位置 rm+1，能量 Em+1和运动方向Ωm+1。
碰撞点位置的计算公式
设 rm 为粒子第 m 次碰撞点的位置，Ωm 为碰撞后的运动方向，则粒子第 m+1 次碰撞点的位置 rm+1 为：
即
其中 为 的方向余弦，L 为两次碰撞点间的距离。
L 的分布密度函数为：
由 f (L) 抽样确定 L 的方法通常有三种：
直接抽样方法
确定 L 的直接抽样方法是：
首先由自由程分布
中抽取ρ
再由下列关系式解出 L 。
对于均匀介质，有
对于多层介质，如果
则
其中， 为粒子由 rm 出发，沿Ωm 方向在顺序经过的第 i 个介质区域内走过的距离， 为第 i 个介质
区域的宏观总截面 ( i =1，2，…，Imax )。 当
时，意味着粒子穿出系统。
最大截面法
对于多层介质，或其他介质密度与位置有关的问题，在求 ( i =1，2，…，Imax ) 时，如果系统形状复杂，计算是非常烦杂的。在这种情况下，使用最大截面法更方便。最大截面抽样方法为：
其中
≤
＞
限制抽样法
当介质区域很小时，如使用直接抽样法抽取输运长度，粒子很容易穿出介质，此时使用限制抽样法确定自由程个数ρ较好，ρ的分布密度函数为：
其中 Dm 为粒子由rm 出发，沿Ωm 方向到达区域边界的自由程个数。ρ的抽样方法是：
然后用直接抽样法中根据ρ计算 L 的方法计算输运长度 L 。此时，粒子的权重需乘以纠偏因子 。
碰撞后能量Em+1的随机抽样
粒子在介质中发生碰撞后，首先要确定与哪种原子核发生何种反应。粒子发生碰撞后（吸收除外）的能量 Em+1 一般只与其碰撞前后运动方向的夹角（散射角）有关。
粒子碰撞后常见的能量分布有下面几种情况。
裂变中子谱
中子引起原子核裂变反应时，裂变中子的能量服从裂变谱分布。其抽样方法可参考以前的介绍。
中子弹性散射后能量的确定
中子弹性散射后，能量与质心系散射角θC的关系是：
能量与实验室系散射角θL的关系是：
其中，A 为碰撞核的质量， 。
或 确定后，即可求出 Em+1。
中子非弹性散射后能量的确定
中子非弹性散射后，能量与质心系散射角θC的关系是：
其中， 为第 K 个能级的阈能， 为第 K 个能级的激发态能量。
如果确定了实验室系散射角θL，则根据下式
确定 后，再计算出 Em+1。
光子康普顿（Compton）散射后能量的确定
光子发生康普顿散射后，其能量分布密度函数为：
其中， K(α) 为归一因子。
， 和 分别为光子散射前后的能量，以 m0c2 为单位，m0为电子静止质量，c 为光速。
光子康普顿散射能量分布的抽样方法为：
x 的抽样确定后，散射后的能量为：
＞
＞
＞
≤
≤
≤
碰撞后散射角的随机抽样
粒子碰撞后运动方向Ωm+1的确定，一般与散射角有关。由已知分布抽样确定散射角后，再确定Ωm+1。常见的散射角分布有如下几种：
质心系各向同性分布
散射角在质心系服从各向同性分布时，其抽样方
法为 。质心系散射角θC抽样确定后，
需转换成实验室系散射角θL。
在中子弹性散射情况下，转换公式为：
其中 A 为碰撞核质量， 。
在中子非弹性散射情况下，转换公式为：
其中， 为第 K 个能级的阈能。
中子弹性散射勒让德 (Legendre) 多项式分布
中子弹性散射角分布常以勒让德多项式的展开形式给定。散射角余弦 x＝cosθ的分布密 度函数为：
其中 Pl(x) 为 l 阶勒让德多项式。
该分布即为 n 阶勒让德近似展开。
勒让德多项式由以下递推公式确定：
考虑新的分布：
当选取 x0，x1，… xn 为 Pn+1(x)＝0 的根，且
时，fa(x) 依照勒让德多项式展开的前 n 项与 f (x) 的展开形式相同。因此，可以用 fa(x) 作为 f (x) 的近似分布。
在实际问题中，由于勒让德多项式展开项数不够，可能出现某个 为负值的现象。此时可以采用如下近似分布：
其中：
对于该近似分布，可用加抽样方法进行抽样：
此时，由于偏倚抽样而引起的纠偏因子为 wK ，也就是说，粒子的权重要乘上wK。
光子康普顿散射角分布
光子的康普顿散射角与其散射前后的能量有关 , 它的分布密度函数为：
抽样方法为：
碰撞后运动方向Ωm+1的确定
实验室系散射角θL确定后，依据不同的坐标系的表现形式，有不同的确定方法。
确定方向余弦 um+1，vm+1，wm+1
其中，
方位角 在 [0, 2π] 上均匀分布。
当 时，不能使用上述公式，可用下面的简单公式：
确定Ωm+1的球坐标 (θm+1，φm+1)
设Ωm的球坐标分别为 (θm,φm)，其中，θ为粒子运动方向与 z 轴的夹角， φ为粒子运动方向在 x y 平面上投影的方位角。则Ωm+1的球坐标 (θm+1,φm+1) 分别由下式确定：
球形几何的随机游动公式
一般几何的随机游动公式可以应用到球形几何，而对球对称问题，使用特殊形式更为方便。
下次碰撞点的径向位置 rm+1的确定
两次碰撞点间的距离 L 确定之后，下次碰撞点的径向位置 rm+1的计算公式为：
设系统的外半径为R，如 rm+1≥R，则粒子逃出系统。
粒子碰撞后瞬时运动方向的确定
在球对称系统中，粒子运动方向用其与径向夹角余弦来描述。使用球面三角公式，粒子碰撞后瞬时运动方向与径向夹角余弦 cosθm+1的计算公式为：
其中， 为在 [0, 2π] 上均匀分布的方位角， 为在 rm+1 点进入碰撞前瞬时运动方向与 rm+1 径向之间的夹角。
点到给定边界面的距离
在抽样确定输运距离、判断粒子是否穿透系统时， 常遇到求由 rm 出发，沿Ωm 方向到达某个区域表面的距离问题。在记录对结果的贡献时，也常使用类似的量。区域表面通常是平面或二次曲面。 求到达区域表面的距离问题，实际上是求直线（或半射线）与平面或二次曲面的交点问题。这是 蒙特卡罗方法解粒子输运的各种实际问题时 , 所遇到的基本几何问题。
点到平面的距离
点 沿方向 的直线方程为：
该直线到达方程为
的平面的距离为：
当 与平面平行时，即
直线与平面无交点。
如果 l 为负值，直线与平面也无交点。这时，粒子的运动方向是背离平面的。
点到球面的距离
在三维直角坐标系中，设球心为 rc＝(xc, yc, zc) ，球半径为 R，则球面方程为：
将直线方程代入该球面方程，得到点 r0沿Ω0方向到达球面的距离 l ：
其中
当 R0≤R 时，即 r0 点在球内 ，Δ≥δ2，l 只有一个正根：
当 R0＞R 时，即 r0 点在球外，分以下三种情况：
若δ≥0，l 无正实根，直线与球面无交点。
若δ＜0，Δ＜0，l 无实根，直线与球面无交点。
若δ＜0，Δ≥0，l 有两个正实根，直线与球面有两个交点。
在球坐标系中，不失一般性，设球心为 rc＝0，则球面方程为 r＝R。
当 r0≤R 时，即 r0 点在球内 ，有一个交点：
其中θ0为Ω0与 r0 的径向夹角。
当 r0＞R 时，即 r0 点在球外 ，令
当 cosθ0≥0 时，直线与球面无交点。
当 cosθ0＜0 时，若 d≥R，则直线与球面无交点。
若 d＜R，则有两个交点：
点到圆柱面的距离
设圆柱面的方程为：
其中 (xc, yc, 0) 为圆柱的中心，R 为圆柱底半径。
点 r0沿Ω0方向到达圆柱面的距离 l 为：
其中
当 R0≤R 时，r0 点在圆柱内 ，如果 ，则 l 有一个正根：
如果 ，即Ω0平行于圆柱的对称轴，直线与圆柱面无交点。
当 R0＞R 时，r0 点在圆柱外，分以下三种情况：
若δ≥0，l 无正实根，直线与圆柱面无交点。
若δ＜0，Δ＜0，l 无实根，直线与圆柱面无交点。
若δ＜0，Δ≥0 且 ，l 有两个正实根，直线与圆柱面有两个交点。
在 的情况下，直线与圆柱面不相交。
点到圆锥面的距离
设圆锥顶点在原点，以 z 轴为对称轴，则圆锥面的方程为：
点 r0沿Ω0方向到达圆锥面的距离 l 为：
其中
如果Ω0与锥面某一母线平行，即 ，则
空腔处理
在粒子输运问题中，所考虑的系统常有空腔存在，如中空的球壳 , 平板间的空隙等。粒子输运时，可有两种处理空腔的方法：
将空腔作为宏观总截面Σt＝0 的区域 , 按通常的方法输运。
设 分别为由 rm 出发，沿Ωm 方向到空腔区域的
近端和远端的交点，当粒子超过 时，以 为新的起
点，重新开始输运。
显然，这两种方法在统计上是等价的。
等效的边界条件
全反射边界
在反应堆活性区中，元件盒常常按正方形或六角形排列。假定元件盒足够多，每个盒结构相同，那么活性区中每个盒所占的栅胞的物理情况，可以代表整个活性区中的状况。
进一步假定，元件盒是圆对称的，那么每个栅胞中情况，可以用更小的单位（栅元）来反映。比如对六角形栅胞可取其 1/12 的ΔOAB 来做代表；正方形栅胞可用其 1/8 的ΔO'A'B' 来做代表。这样一来问题就大大简化了。
现在的问题是怎样计算直角三角形栅元的物理量（如通量）。用蒙特卡罗方法如何模拟中子在栅元内的运动，反映出整个活性区对它的影响。
我们可把OA'、OB'、A'B' 作为全反射边界来处理。所谓全反射边界，就是当中子打到该边界上时，按镜面反射的方式，从边界 上全部反射回来，中子的能量与权重均不改变。
在这种边界上的反射条件，称之为全反射条件，就是通常的镜面反射条件。
在全反射边界条件下，一条通过活性区若干个区域的中子径迹，可以用栅元ΔO'A'B' 中的一条折线轨 道来反映出来。
反过来，在直角三角形栅元ΔO'A'B' 中任一条反射成的折线轨道，都代表了中子在活性区内一条直线轨道的作用。由于系统的对称性，在活性区内，凡是与栅元内位置相当的地方，都有相同的物理情况，因此栅元内各处的情况，当然代表了整个活性区的情况。
一般曲面全反射条件
对于一般曲面的全反射，设入射方向为Ω，入射点的内法线方向为 n ，则反射方向Ω' 为：
其中
设
则
平面全反射条件
设三角形栅元的横截面ΔOAB 在 X-Y 平面上，∠OAB＝θ。则边界 OA、OB、AB 上的反射都是平面全反射。在任一与 X-Y 平面垂直且与 X 轴成α角的平面上，全反射条件为：
由此就可得到OA、OB 和 AB 边上的全反射条件，对于 OB 边，α=θ；对于 OA 边，α= 0；对于 AB 边，α=π/ 2。
反射层边界条件
对于具有大反射层的系统，如存放，运输和生产裂变物质的仓 库、车厢和车间等，当中子从里面打到四周墙上或反射层时，还要继续对它进行跟踪。这种跟踪常常要花费很大的计算量，并且在结果中引起的方差也比较大。如果在计算这种系统的不同方案中，反 射层条件不变，那么这种大量重复的计算是很不经济的。
中子射入反射层后，一部分被介质吸收，只有一部分返回，由于中子的散射慢化，损失一部分能量，因此反射回来的中子有一个能量方向分布。显然，对这种反射层，不能应用全反射条件。不过，我们仍然可以把它当做边界，在边界上按反射层的物理作用来处理。
比如，如果反射层是一种平板几何，我们可以用数值方法或蒙特卡罗方法，预先算好在各种不同入射能量 E 下的反照率β(E)，反射中子的能量分布 RE(E→E' )。于是代替在反射层中眼踪中子，我们可在反射层边界上作如下处理：
一旦中子打入反射 层，立即返回，反射后权重为
其中，E 为射入反射层中子的能量，W 为中子的权重。反射后的能量 Eβ 由反射能谱 RE(E→Eβ) 中抽样产生。反射后的方向Ωβ 由半平面各向同性分布或余弦分布中抽样。反射后的中子位置为入射时的位置。
计算表明，对于大尺寸的反射层来说，这样的近似，引 起的结果上的误差是可以忽略的，却能带来计算量的大量节省。
记录贡献与分析结果过程
在粒子输运问题中，除了得到某些量的积分结果外，还需要得到这些量的方差、协方差、以及这些量的空间、能量、方向和时间的分布。这些量可以利用分类记录手续同时得到。
记录与结果
为了得到所求量的估计值，在粒子输运过程中需进行记录，即求每个粒子对所求量的贡 献。
设模拟了 N 个粒子，所求量的估计值为：
其中 gi 为第 i 个粒子的总贡献。
记录的贡献由所求量决定。对于同一个所求量，又随所用的蒙特卡罗技巧的不同而不同。 例如，所求量是粒子穿透屏蔽概率，使用直接模拟法时，如粒子穿透屏蔽，在叠加记录单元加“1” ( 初始值为零 )，否则没有贡献。使用加权法时，如粒子穿透屏蔽，在叠加记录单元加粒子的权重，否则没有贡献。使用统计估计法时，粒子每发生一次碰撞 (包括零次碰撞)，都要记录贡献，等等。
方差和协方差的估计
估计量 g 和 g' 的方差和协方差为：
它们可以用下式估计：
因此，要得到 和 的估计，只要对每一个历史记
录结果的 和 进行记录，并加以累加即可。
方差估计值 确定后，可得到误差
其中 为置信限，它随置信水平 而定。在通常
情况下，取 。
位置、能量、方向、时间分布
在前面已经提到，用蒙特卡罗方法求某种量的空间、能量、方向和时间分布，实质上是得到这种分布的阶梯函数近似的估计值。而求这种估计值是很方便的，只要将跟踪过程中所得到的感兴趣量，按其状态的空间、能量、方向、时间特征，分别记录其权 重，最后将这些记录结果适当处理即可。
事先，将问题的空间、能量、方向（常按相对于某个方向的夹角余弦）、时间范围，各分为如下不同间隔：
再用一批存贮单元 {A} 记录相应间隔上阶梯函数近似的累计值。
核截面数据的引用
用蒙特卡罗方法解粒子输运问题，需要介质所包含的各种原子核的核数据。以中子核数据为例，需要各种涉及到的核的微观总截面、弹性散射、非弹性散射、n-2n 反应、裂变、俘获等截面；也需要这些反应的相应能量、角度分布、次级粒子数，以及其它关心的粒子数及其能量、方向分布。从输运方程中可以知道，有了这些数据，问题就完全确定了；反映到蒙特卡罗模拟中，有了这些数据， 就能决定宏观总截面，决定碰撞核的具体形式，就能实现抽样和跟踪。
在蒙特卡罗计算中，引用的核数据有点截面、分段常数截面和群截面三种形式。
点截面形式
在跟踪粒子时，对粒子的每一种能量，先从截面库中取出需要的核数据，再用插值（或其它方式），求出相应能量的各种截面数据。这种方法是直接的，也是比较精确的。不少通用程序就是这样做的。这样做的条件是要有完备的核截面数据库，计算机有大而快的存贮能力。
分段常数形式
将问题的能量范围（Emin，Emax）分成许多间隔，截面数据在每个间隔上看作与能量无关，即截面取分段常数形式。这种引用截面数据的好处是，数据量相对地少。问题是它要根据问题的物理特点来划分能量间隔。而为了保证误差较小，所取的间隔数一般是比较多的。
群截面形式
解中子输运问题，常常采用多群近似方法。在多群近似下，中子能量 E 用群指标 i 代替。为了实现蒙特卡罗跟踪，只需要以下群截面：
显然，这种跟踪过程数据量小，程序简单。
蒙特卡罗程序结构
在电子计算机上，蒙特卡罗方法解粒子输运问题的程序，一般都可分为：源抽样，空间输运过程，碰撞过程，记录过程和结果的处理与输出等部分。
开始
数据预处理, 各记录单元清零
取一个粒子历史
源分布抽样
输运过程
碰撞过程
历史终止否？
统计处理
做完给定历史数否？
结果的处理与输出
终止
记录过程
记录过程
记录过程
记录过程
至于粒子历史终止条件，根据问题的几何条件、物理假定，处理方法，可归纳为以下几种：
粒子从系统逃脱；
粒子经碰撞被吸收；
经俄国轮盘赌后，历史被终止；
粒子能量低于给定能量（阈能） ;
粒子位置越过某一界面；
粒子飞行时间超过给定时间；
粒子权重小于某个小量。
作 业
第六章 蒙特卡罗方法在通量计算中的应用
通量的定义
通量的能谱和角分布
计算体通量的模拟方法
计算面通量的模拟方法
计算点通量的模拟方法
与通量有关的物理量的计算
作 业
第六章 蒙特卡罗方法在通量计算中的应用
通量计算在粒子输运问题中占有非常重要的地位。很多问题，如碰撞率、反应率以及系统逃脱几率等都可以通过通量来计算。通量计算问题，包括点通量、面通量和体通量的计算问题。相对来说，点通量的计算要困难一些。
通量的定义
设 分别表示粒子的位置、能量和运动方向。则通量 的定义为：
在 r 点的体积元 dV 内，能量E 和运动方向Ω属于dE dΩ的粒子平均径迹长度。
点通量的定义
给定点 r0 的点通量为：
点通量的含义为：
在r0点的体积元dV内，粒子的平均径迹长度。
面通量的定义
给定曲面 A0 上的面通量为：
面通量的含义为：
沿曲面A0的法线方向增加厚度ds 所组成的体积元的体积元A0ds中，粒子的平均径迹长度。
体通量的定义
给定体 V0 内的体通量为：
体通量的含义为：
在体V0内，粒子的平均径迹长度。
粒子各次散射对通量的贡献
通量 可用粒子各次散射对通量的贡献和表示：
其中 为粒子 n 次散射后对通量的贡献，其含义为：
粒子在第 n 次散射到第 n＋1 次散射之间，在 r 点的体积元 dV 内，能量E 和运动方向Ω属于dE dΩ的粒子平均径迹长度。
通量的能谱与角分布
用蒙特卡罗方法计算通量的能谱与角分布，所采用的手段与计算其它物理量一样，即把能量和方向分成若干个区间，分别按粒子状态所处的区间累积记录各自的贡献。
现将能量分成 I 区：ΔE1，ΔE2，…，ΔEI；方向分成 J 区：ΔΩ1，ΔΩ2，…，ΔΩI。
则有：
计算体通量的模拟方法
在实际问题中，经常遇到要计算某一区域V0 的体通量。
在通量的定义部分已经介绍过，通量可以表示为粒子各次散射对通量的贡献和。因此，下面要介绍的各种估计方法，只叙述各次散射后的通量计算方法。
计算体通量的方法主要有以下几种。
解析(统计)估计方法
粒子 n 次散射（n＝0 时为源粒子）后的通量贡献为：
其中，s1和s2分别为粒子由点rn出发，沿Ωn方向到达
区 域V0的近端和远端的交点的距离。如果点rn在V0内，
则 s1＝0。如果粒子沿Ωn方向与V0有多段相交，则
为每段相交线段的通量贡献之和。如果粒子沿
Ωn方向与V0不相交，则 。
解析估计方法就是把体通量的贡献表达式直接计算出来。当系统为均匀介质时，
如果只是V0为均匀介质，则
如果V0由多层介质组成，则需分段计算积分。
在解析估计方法中，粒子每发生一次碰撞（包括零次散射），都要记录通量的贡献值。
径迹长度方法
设粒子从第 n 次散射到第 n＋1 次散射之间走过的径迹长度为 s ，则 n 次散射的通量贡献为：
径迹长度方法就是把粒子在V0内走过的径迹长度记录下来。
下面证明，径迹长度估计是无偏的。
碰撞密度方法
设粒子从第 n 次散射到第 n＋1 次散射之间走过的径迹长度为 s ，则 n 次散射的通量贡献为：
碰撞密度方法就是把粒子在V0内发生的碰撞记录下来。
下面证明，碰撞密度估计是无偏的。
均匀径迹长度方法
确定一个定义在 [s1,s2] 上的概率密度函数 fn(s)，从 fn(s)中抽样 s*，则 n 次散射通量贡献的估计为：
fn(s)的最简单形式是均匀分布
这时
点通量代替方法
设 为在V0上定义的任一概率密度函数，则体通量可表示为：
体通量的估计为：
其中，r*为从 中抽取的一个样本值。
几种方法的比较
解析估计方法：直接计算体通量的贡献表达式，因此该方法的方差小，但计算时间长，需要计算指数函数的积分。
径迹长度方法：记录贡献方法简单，可与输运过程同时进行，只要粒子穿过记录区域就有贡献。但该方法方差大些，对于较小的系统（如自由程个数小于2），该方法较好。
碰撞密度方法：由于只在记录区域内发生碰撞才有贡献，因此方差较大，尤其在记录区域较小时更是如此。但该方法省时间，适用于大的记录区域。
均匀径迹长度方法：在记录区域为多层介质时，较解析估计方法容易实现。但在记录贡献时仍需计算指数函数，也费时间。
点通量代替方法：可以较好地解决小区域的体通量计算问题。尤其是记录区域与粒子的输运区域分开时，更是如此。
计算面通量的模拟方法
计算面通量的方法主要有以下几种。
解析估计方法
设经过 n 次散射的粒子，由点rn出发，沿Ωn方向
到达曲面域A0的距离为 s1，与曲面相交处曲面的法线
方向为 n，则 n 次散射粒子对该曲面的通量贡献为：
如果粒子沿Ωn方向与A0有多个交点，则 为每个
交点处的通量贡献之和。如果粒子沿Ωn方向与A0没有
交点，则 。
解析估计方法就是把面通量的贡献表达式直接计算出来。粒子每发生一次碰撞（包括零次散射），都要记录通量的贡献值。
加权（径迹长度）方法
设粒子从第 n 次散射到第 n＋1 次散射之间走过的径迹长度为 s ，则 n 次散射的通量贡献为：
加权方法只有在粒子穿过曲面A0时，才对该曲面有通量贡献。
点通量代替方法
设 为在A0上定义的任一概率密度函数，则面通量可表示为：
面通量的估计为：
其中，r*为从 中抽取的一个样本值。
体通量代替方法
沿曲面A0的法线方向均匀地增加一个厚度Δs，由此构成的体积为 。 的体通量为：
A0的面通量为：
因此，如取得足够小，有如下近似：
计算点通量的模拟方法
与体通量、面通量的计算相比，点通量的计算最困难。这是因为，在大量的模拟粒子中，只能有很少的粒子穿过该点所包含的一个小区域，因此无法使用通常的通量计算方法。
指向概率方法
设 n 次散射后粒子的状态为 ，
进入 n 次碰撞的粒子的状态为 ，
表示粒子的碰撞核，其定义为：
一个粒子在点 r 发生碰撞后，能量由E'变为E的dE内，方向由Ω'变为Ω的dΩ内的粒子平均数。
则 n 次散射的粒子对点 r* 的通量贡献为：
其中
当 n＝0时，用源分布密度函数 代替碰撞
核 。
光子问题的指向概率方法
光子问题的碰撞核为：
其中光子能量E以电子静止能量mec2≈ MeV为单位；K(E'→E/r) 为Klein－Nishina公式，由下式确定
N(r) 表示在 r 处单位立方体内的原子数，z (r) 表示在 r 处元素的原子序数，r0表示电子的经典半径。
其中
中子问题的指向概率方法
中子问题的碰撞核为：
其中下标A和 i 分别表示不同的原子核和不同的反应；
和 分别表示能量为E'的中子与第A种原子核发生第 i 种反应后产生的平均次级中子数和微观截面；NA(r) 表示在 r 处第A种原子核的核密度；
表示能量为E'和方向为Ω'的中子与第A种原子核发生第 i 种反应后的能量E和方向Ω的分布。
则有中子的通量贡献为：
关于估计量无界问题
当 r*点附近不含散射物质时（如真空），也就是说，粒子的输运区域与记录点分开时，指向概率方法的估计量是有界的，因此是一种比较好的计算点通量的方法。不含散射物质的区域越大，指向概率方法的优点越明显。
然而，当 r*点附近含有散射物质时，由于在指向概率方法的估计量中含有无界因子
因此，指向概率方法的估计量一般来说是无界的。
与通量有关的物理量的计算
一些物理量可以通过通量来计算。
系统逃脱概率
状态为 的粒子的系统逃脱概率为：
系统逃脱概率为：
各种反应率
碰撞密度
裂变中子密度
吸收率
反应率
作 业
第七章 蒙特卡罗方法在积分计算中的应用
蒙特卡罗方法求积分
重要抽样
俄国轮盘赌和分裂
半解析方法
系统抽样
分层抽样
第七章 蒙特卡罗方法在积分计算中的应用
计算多重积分是蒙特卡罗方法的重要应用领域之一。本章着重介绍计算定积分的蒙特卡罗方法的各种基本技巧，而这些技巧在粒子输运问题中也是适用的。
蒙特卡罗方法求积分
蒙特卡罗方法求积分的一般规则如下：任何一个积分，都可看作某个随机变量的期望值，因此，可以用这个随机变量的平均值来近似它。
设欲求积分
其中，P＝P(x1，x2，…，xs) 表示 s 维空间的点，Vs表示积分区域。取Vs上任一联合概率密度函数 f (P)，令
则
即θ是随机变量 g(P) 的数学期望，P的分布密度函数为 f (P) 。 现从 f (P) 中抽取随机向量 P 的 N 个样本：Pi，i＝1，2，…，N， 则
就是θ的近似估计。
重要抽样
偏倚抽样和权重因子
取Vs上任一联合概率密度函数 f1(P)，令
则有
现从 f1(P) 中抽样 N 个点：Pi，i＝1，2，…，N， 则
就是θ的又一个无偏估计。
重要抽样和零方差技巧
要使 最小，就是使泛函I[f1] 极小。
利用变分原理，可以得到最优的 f1(P) 为
特别地，当 g(P)≥0 时，有
这时
即 g1的方差为零。实际上，这时有
不管那种情况，我们称从最优分布 fl(P)的抽样为重要抽样，称函数 | g(P) | 为重要函数。
俄国轮盘赌和分裂
分裂
设整数 n≥1，令
则
于是计算θ的问题，可化为计算 n 个θi 的和来得到，而每个 gi(P) 为原来θ的估计 g(P) 的 1/ n ，这就是分裂技巧。
俄国轮盘赌
令 0 ＜ q＜1，
则
于是θ变为一个两点分布的随机变量ζ的期望值，
ζ的特性为：
这样就可以通过模拟这个概率模型来得到θ，这就是俄国轮盘赌。
重要区域和不重要区域
我们往往称对积分θ贡献大的积分区域为重要区域，或感兴趣的区域；称对积分θ贡献小的区域为不重要区域，或不感兴趣的区域。
考虑二重积分
令R是V2上 x 的积分区域，表为 R＝R1+R2，其中R1是重要区域，R2是不重要区域，两者互不相交。又命Q为V2上相应于 y 的积分区域。则
通常蒙特卡罗方法，由f (x,y)抽样 (x,y)的步骤是：从 fl(x) 中抽取 xi，再由 f2(y|xi) 中抽样确定 yi，然后用
作为θ的一个无偏估计。
现在，改变抽样方案如下：
当x∈R1时，定义一个整数n(xi)≥1，对一个xi，抽取
n(xi)个yij，j＝1，2，…，n(xi)。以平均值
代替上述θ估计式中的 g(yi, xi) 。
当 x∈R2时，定义一个函数q(xi)，0＜ q(xi) ＜1，
以抽样值
代替上述θ估计式中的 g(yi, xi) 。这里ξ是随机数。
显然，这种抽样估计技巧，就是对 x∈R1时，利用分裂技巧，而对 x∈R2时，利用俄国轮盘赌，而使估计的期望值不变。由于对重要区域多抽样，对不重要区域少观察，因此能使估计的有效性增高。
半解析（数值）方法
考虑二重积分
令
则θx为θ的无偏估计。
θx 的方差为
而由 f (x,y)抽样 (x,y)，用 g (x,y)作为θ的估计，其方差为
系统抽样
我们知道，由f (x,y)抽样 (x,y)的步骤是：
从 fl(x) 中抽取 xi，
再由 f2(y|xi) 中抽样确定 yi，
现在改变 xi 的抽样方法如下：
yi 的抽样方法不变。
其方差为
与通常蒙特卡罗方法相比，方差减少了约
分层抽样
考虑积分
在（0，1）间插入J－1个点
0＝α0＜ α1＜ …＜ αJ-1＜ αJ＝1
令
则有
现在，用蒙特卡罗方法计算θj ，对每个θj 利用 fj(x)中的nj 个样本xij ，那么有
第八章 蒙特卡罗方法应用程序介绍
蒙特卡罗方法应用软件的特点
常用的通用蒙特卡罗程序简介
MCNP程序输入的描述
第八章 蒙特卡罗方法应用程序介绍
建立完善的通用蒙特卡罗程序可以避免大量的重复性工作，并且可以在程序的基础上，开展对于蒙特卡罗方法技巧的研究以及对于计算结果的改进和修正的研究，而这些研究成果反过来又可以进一步完善蒙特卡罗程序。
蒙特卡罗方法应用软件的特点
通用蒙特卡罗程序通常具有以下特点：
具有灵活的几何处理能力
参数通用化，使用方便
元素和介质材料数据齐全
能量范围广，功能强，输出量灵活全面
含有简单可靠又能普遍适用的抽样技巧
具有较强的绘图功能
常用的通用蒙特卡罗程序简介
MORSE程序
较早开发的通用蒙特卡罗程序，可以解决中子、光子、中子－光子的联合输运问题。采用组合几何结构，使用群截面数据，程序中包括了几种重要抽样技巧，如俄国轮盘赌和分裂技巧，指数变换技巧，统计估计技巧和能量偏移抽样等。程序提供用户程序，用户可根据需要编写源分布以及记录程序。
EGS程序
EGS是Electron-Gamma Shower 的缩写，它是一个用蒙特卡罗方法模拟在任意几何中，能量从几个KeV到几个TeV的电子-光子簇射过程的通用程序包。由美国Stanford Linear Accelerator Center提供。EGS于1979年第一次公开发表，提供使用。EGS4是1986年发表的EGS程序的最新版本。
MCNP程序
MCNP是美国Los Alamos国家实验室开发的大型多功能通用蒙特卡罗程序，可以计算中子、光子和电子的联合输运问题以及临界问题，中子能量范围从10-11MeV至20MeV，光子和电子的能量范围从1KeV至1000MeV 。程序采用独特的曲面组合几何结构，使用点截面数据，程序通用性较强，与其它程序相比，MCNP程序中的减方差技巧是比较多而全的。
MCNP程序输入的描述
MCNP的输入包括几个文件，但主要的一个是由用户编写的INP文件，该文件包括描述问题所必须的全部输入信息。文件采用卡片结构，每行代表一张卡片，文件由一系列卡片组成，对于任一特定的问题，只需用到INP全部输入卡片的一小部分。
MCNP输入文件中物理量的单位
×1023
阿伏加德罗常数
中子质量的倍
原子量
10-24 厘米2
截面
克 / 厘米3
质量密度
1024 个原子 / 厘米3
原子密度
MeV（kT）
温度
10-8 秒
时间
MeV
能量
厘米
长度
输入文件的基本形式
信息块
信息块的卡片放在INP文件中标题卡之前。信息块给出了MCNP的一些运行信息，信息块上各部分的意思和运行行信息是一样的，当运行行信息与信息块中所指定的信息相矛盾时，则忽略信息块中相应的信息，而以运行行信息为准。
信息块是可选的，信息块的第一张卡片，必须在第 1～8 列写上“MESSAGE:”，从第一张卡片的第 9～80 列到后续卡片的第 1～80 列都可填写运行信息。在标题卡之前用一个空行分隔符结束信息块。
初始运行的输入文件
选择项
其它
其它数据，包括问题类型、源描述、材料描述、计数描述，问题截断条件等。
数据卡
…
空行分隔符
定义组成栅元的曲面信息。
曲面卡
…
空行分隔符
定义构成整个系统的各个基本介质单元以及相应的物理信息。
栅元卡
…
空行分隔符
仅一行，占用第 1～80 列。作为输出标题。
标题卡
选择项
信息块
空行分隔符
接续运行的输入文件
接续运行必须在运行行信息或信息块中给出C项选择，即Cm，表示从RUNTPE文件中读出第m次转储的内容接着运算，如果m未指定，则读最后一次转储的数据。如果不需要改变内容，则不需要接续输入文件，仅需运行RUNTPE以及在运行行加上C选择。
选择项
其它
只允许部分数据卡。（FQ，DD，NPS，CTME，IDUM，RDUM，PRDMP，LOST，DBCN，PRINT，KCODE，MPLOT，ZA，ZB，和 ZC）
数据卡
…
空行分隔符
写在第 1～8 列
CONTINUE
选择项
信息块
空行分隔符
卡片格式
INP输入文件的每一行（称之为一张卡片）都限于使用第 1～80 列并构成卡片映象。大部分输入卡片按行填写；然而，对数据卡允许按列填写。 $ 符号为它所在那行数据的结束符，在 $ 符号后面的内容作为注释，它可从 $ 符号后面的任一列开始。
标题卡只占一行，整行都可填入用户需要的信息，也可以是空行。但要注意在其它地方使用空行是作为结束符或者分隔符。
输入文件中，在标题卡之后及最后的空行结束卡之前的任何地方都可插入注释卡。注释卡必须是字母“C”写在 1～5 列中的任意位置，且至少用一个空格隔开后面的注释内容。
行输入格式
栅元卡、曲面卡和数据卡的书写格式是相同的。必须从 1～5 列开始填写这些卡片相应的名字(或编号)和粒子标识符，后面填写用空格分隔的数据项。如果 1～5 列为空，则表示它是前一张卡片的继续卡。如果在一行的末尾有一个用空格隔开的符号“&”，则表示下一行是该行的继续卡，数据可填写在 1～80 列。一个数据项必须在一张卡片上写完，不得跨到下一张卡片上。完全空白的一行则为两组卡片的分隔符。 对任何给定的带有粒子标识符的类型卡只能有一张。需要整数的数据项必须填写整数，其它数据可填写为整数或浮点数以及MCNP能读的数据。
为书写方便，可以使用四项书写功能：
nR功能，表示将它前面的数据重复n次。
例如：2 4R 等同于 2 2 2 2 2
nI功能，表示在与其前后相邻的两个数之间，插入n个线性插值点。对于 X nI Y 的结构，如果X和Y是整数，且X－Y刚好是n+1的整倍数，则产生标准的整数插值，否则产生实数插值，但Y值直接存储。
例如： 2I ＝ 3 可能不精确
而 1 4I 6 ＝ 1 2 3 4 5 6 都是精确定整数
XM功能，它表示的数值为前面的数据乘上X。
例如：1 1 2M 2M 4M 2M ＝ 1 1 2 4 16 32
nJ功能，表示其后n个数据项使用缺省值。
例如：DD .1（缺省值） 1000 ＝ DD J 1000
如果nR、nI、及nJ项中缺省n，则假设n＝1。
列输入格式
列输入块的格式：
Si必须是MCNP卡片名字，它们必须全部是栅元参数、或者全部是曲面参数、或者全部是其它参数。
Dn1 Dn2 … Dnm
Kn
D21 D22 … D2m
K2
…
…
D11 D12 … D1m
K1
S1 S2 … Sm
＃
6～72 列
1～5 列
粒子标识符
几个输入卡片都需要粒子标识符以区别中子、光子和电子的输入数据。这些卡片是：IMP、EXT、FCL、WWN、WWE、WWP、WWGE、DXT、DXC、F、F5X、F5Y、F5Z、PHYS、ELPT、ESPLT、CUT和PERT。粒子标识符由上述卡片名字后面的冒号、字母N、P或E组成。
例如：中子重要性卡为 IMP:N
光子重要性卡为 IMP:P
缺省值
MCNP的许多输入参数都有缺省值，因此用户不需要每次都给出各个输入参量的值。当缺省值符合用户要求时，便可不在输入文件中指定。当省略某张输入卡时，则该卡上的全部参数均使用缺省值。如果只想改变一张卡上的某一个特定参量时，则它前面的参量仍需指明，或者用nJ方式跳过前面那些使用缺省值的参量。
例如：光子截断卡 CUT:P 3J
表示前3个参量使用缺省值，只改变第四项参量的值。
输入错误信息
MCNP对输入文件出现的错误作广泛的检查，如果用户违反了输入说明的规定，将在终端上以及输出文件中打印致命错误信息，MCNP不再进行粒子输运计算，作业中断。
第一个出现的致命错误是真的，而后面的错误可能不一定是真的，这取决于前面出现的致命错误的情况。
若在MCNP运行行上指定FATAL项，则MCNP忽略致命错误，照常运行。
对于MCNP的警告信息，用户不应忽视，应搞清楚它们的含义。
检查几何错误
MCNP在处理输入文件的数据时，不能检查一种非常重要的输入错误。即MCNP无法查出各栅元之间的重叠和空隙，只有当粒子丢失时，才会发现几何错误。即使如此，可能仍然无法准确判断错误性质。
栅元描述卡
另一个栅元的名字（编号）。
n
任选的栅元参数说明。
params
栅元材料号，与材料卡（Mm）中的序号对应。
m ＝0 为真空栅元。
m
或：
格式：
描述栅元j和栅元n之间差别的栅元参数。
list
栅元的几何说明。由一系列带符号的曲面号经过布尔运算组成。
geom
栅元材料密度。正值为原子密度，负值为质量密度。对于真空栅元，该项缺省，不填写。
d
栅元号，1≤ j ≤99999，写在第 1～5 列上。
j
LIKE n BUT list
j
m d geom params
j
在栅元的几何说明中，关于曲面的指向是一个很重要的概念。假定曲面 S 的曲面方程为 f (x,y,z)＝0，则对于f (x,y,z)＞0的区域对于曲面 S 具有正的指向；而对于f (x,y,z)＜0的区域对于曲面 S 具有负的指向。正指向的区域用+S表示，“+”号可不写；负指向的区域用-S表示。栅元用各相关曲面的布尔运算表示，布尔算符包括交（用空格表示 ）、并（用冒号:表示）和非（用#表示）。缺省的运算顺序是先非，其次是交，最后是并 ，使用括号可控制布尔运算的次序。
非运算有两种形式：
#n，n是某个栅元号， #n表示一个由不在栅元n内的点组成的空间区域。
# ( ---)，括号内是对某一个栅元进行描述的曲面——栅元关系组，这一形式定义的几何区域由不属于括号内描述区域的点组成的空间。
例如：3 0 -1 2 -4 $ 定义栅元3
#3 $ 与下行相同
#(-1 2 -4)
在栅元卡上可定义栅元参数以代替在输入文件中数据卡部分定义的栅元参数。格式为：关键词＝值。这儿允许的关键词是：带有粒子标识符的IMP、VOL、PWT、EXT、FCL、WWN、DXC、NONU、PD和TMP，以及关于重复结构的4个栅元参数卡：U卡、TRCL卡、LAT卡和FILL卡。
在LIKE n BUT格式中，还有两个关键词MAT和RHO，分别表示栅元的介质号和密度。
例如：10 16 1 -2 3 IMP:N=4 IMP:P=8
表示栅元10由曲面1的正面、曲面2的负面和 曲面3的正面的交集组成，填充质量密度为 克 / 厘米3 的16号材料。该栅元的中子重要性为4，光子重要性为8。
例如：2 3 -1 IMP:N=2 IMP:P=4
3 LIKE 2 BUT TRCL=1 IMP:N=10
曲面描述卡
由方程定义曲面
曲面助记符。
a
对应坐标变换卡TRn，表示该曲面是在辅助坐标系下描述的，而该辅助坐标系与基本坐标系之间的关系由TRn卡给出。如果没有坐标变换，即曲面是在基本坐标系下描述的，则该项缺省。
n
格式：
曲面方程参数，1～10项，取决于曲面类型。
参见MCNP手册，表。
list
曲面号，1≤ j ≤99999，写在第 1～5 列上。
如果曲面号前有*号，则该曲面为反射面。
j
n a list
j
用点定义轴对称曲面
类型为X、Y或Z的曲面卡是用坐标点描述曲面而不是用方程系数描述。用这些卡描述的曲面必须是分别关于X、Y或Z轴对称的，并且如果该曲面是由多叶组成的，则指定的坐标点必须全都在同一个叶上。
字母X、Y或Z。
a
TRn卡的号，如果没有坐标变换，则该项缺省。
n
格式：
1～3对点的坐标。
list
曲面号，1≤ j ≤99999，写在第 1～5 列上。
j
n a list
j
每一对坐标点定义这个曲面上的一个点。例如在一张Y卡上可以给出：
j Y y1 r1 y2 r2
其中，( ) 是第 i 点的坐标。给出的坐标点对数的不同，描述的曲面类型也不同。
给出一对坐标，则定义一个平面（PX、PY或PZ）。
给出二对坐标，则定义的是线性曲面（PX、PY、PZ、CX、CY、CZ、KX、KY或KZ）。
给出三对坐标，则定义的是二次曲面（PX、PY、PZ、SO、SX、SY、SZ、CX、CY、CZ、KX、KY、KZ或SQ）。
当用两点定义一个锥面时，只生成一个单叶锥面。
曲面的指向与方程指定曲面（SQ除外）是一样的。
对于SQ，远离对称轴的点具有正指向。而方程定义的SQ可以自由选取指向。
由三个点定义一般平面
MCNP对用户指定的 P 型曲面， 将检查所给的数据个数，若是 4项，则作一般斜置平面方程的系数理解，若多于 4 项时，便作为三维空间点的坐标值理解。每三个数定义空间一个点，MCNP将把它们转换成所需要的曲面系数以产生平面：
AX+BY+GZ－D＝0
该曲面卡的助记符。
P
TRn卡的号，如果没有坐标变换，则该项缺省。
n
（Xi,
格式：
定义该平面的点坐标。
Yi, Zi）
曲面号，1≤ j ≤99999，写在第 1～5 列上。
j
n P X1 Y1 Z1 X2 Y2 Z2 X3 Y3 Z3
j
数据卡
在信息卡、栅元描述卡和曲面描述卡之后输入的是数据卡，数据卡可分为10类：
问题类型
几何卡
减方差
源描述
计数描述
材料及截面描述
能量及热处理
问题截断条件
用户数据数组
外围卡
数据卡中，标识符必须从前5列开始填写。
问题类型（MODE）卡
如果不给出MODE卡，则缺省形式是MODE N，即缺省值是中子输运问题。
E，电子输运。
P，光子输运。
格式：
N，中子输运。
xi ＝
x1 … xi
MODE
几何卡
几何卡有以下几类：
填充卡
FILL
Universe
U
栅元变换
TRCL
格子
LAT
坐标变换
TR
曲面面积
AREA
VOL
助记符
栅元体积
卡片类型
坐标变换卡格式：
TRn O1，O2，O3，B1，B2，B3，B4，B5，B6，
B7，B8，B9，M
n ＝变换号， 1≤ n ≤999 ，*TRn表示 Bi是角度而非角度的余弦。
O1，O2，O3 ＝坐标变换向量的位移。
B1至B9 ＝坐标变换的坐标旋转矩阵。
元素 B1， B2，B3，B4， B5，B6，B7，B8，B9
轴 x,x’ y,x’ z,x’ x,y’ y,y’ z,y’ x,z’ y,z’ z,z’
M ＝1，表示位移是辅助坐标系原点
相对于基本坐标系的位移。
＝－1，表示位移是基本坐标系原点
相对于辅助坐标系的位移。
减方差
MCNP运用以下卡片来减小方差：
强迫碰撞
FCL
次级光子权重
PWT
指数变换
EXT
方向矢量定义
VECT
能量分裂和俄国轮盘赌
ESPLT
IMP
助记符
栅元重要性
卡片类型
DXTRAN贡献
DXC
探测器贡献
PD
重叠重要性网格权重生成器
MESH
权重窗的参数
WWP
权重窗生成器
WWG
权重窗生成器的能量或时间间隔
WWGE
韧致辐射偏倚因子
BBREM
权重窗的边界
WWN
WWE
助记符
权重窗的能量或时间间隔
卡片类型
源定义
临界计算的源起始点
KSRC
源的注释
SCn
写曲面源
SSW
读曲面源
SSR
临界源
KCODE
α特征值源
ACODE
源的概率
SPn
源的偏倚
SBn
相关的源
DSn
源的信息
SIn
SDEF
助记符
通用源
卡片类型
通用源卡：
格式
计数描述
下列卡片用来记录计算结果：
计数能量乘子
EMn
剂量能量/剂量函数
DEn/DFn
计数乘子
FMn
计数打印层次
FQn
计数能量间隔
En
计数时间间隔
Tn
计数方向余弦间隔
Cn
计数注释
FCn
Fna
助记符
计数类型
卡片类型
DXTRAN球参数
DXT
计数涨落打印
TFn
探测器和DXTRAN诊断
DD
子程序TALLYX输入
FUn
计数片段的体积/面积
SDn
计数片段划分
FSn
计数余弦乘子
CMn
计数栅元标志
CFn
计数曲面标志
SFn
计数特殊处理
FTn
计数时间乘子
TMn
助记符
卡片类型
计数类型卡Fna格式：
*Fn单位
Fn单位
类型说明
助记符
无
电荷
沉积电荷
+F8: E
MeV
脉冲
F8: (P、E、P,E)
109 J/克
MeV/克
平均裂变沉积能量
F7: N
109 J/克
MeV/克
平均沉积能量
F6: (N、P、N,P)
MeV/cm2
粒子/cm2
点或环探测器通量
F5a: (N、P)
MeV/cm2
粒子/cm2
体通量
F4: (N、P、E)
MeV/cm2
粒子/cm2
面通量
F2: (N、P、E)
MeV
粒子
面流
F1: (N、P、E)
材料描述
这组卡片用于指定在栅元中所使用的材料成分和使用那些截面数据。
多群特征描述
MGOPT
次级光子产生偏倚
PIKMT
否定材料
VOID
截面文件
XSn
总裂变
TOTNU
裂变截断
NONU
原子量
AWTAB
离散反应截面
DRXS
Mm
助记符
材料成分
卡片类型
材料成分卡Mm格式：
Mm ZAID1 fr1 ZAID2 fr2 … keyword=value …
ZAIDi ＝材料中第 i 种成份的截面数据，
或ZZZAAA
ZZZ是元素的原子序号，AAA是原子量， nn截面库标识号，X是数据分类。
fri ＝材料中第 i 种成份的原子的分量
（负值表示重量比例）。
AAA＝000表示自然元素。
能量和热处理方式指定
这组卡片用于控制MCNP的能量以及其它物理状况。
热时间卡
THTME
S(α,β) 材料卡
MTm
自由气体模型热温度
TMP
PHYS
助记符
能量物理截断
卡片类型
问题截断卡
这组卡片在初始运行或接续运行的输入文件中均可使用，用于终止粒子的历史或中断计算。
历史数截断
NPS
计算时间截断
CTIME
栅元能量截断
ELPT
CUT
助记符
截断粒子历史
卡片类型
用户数据数组卡
MCNP在其COMMON变量中定义了两个数组IDUM（整数）和RDUM（浮点数）供用户使用，每个数组可存放50个数据。这组卡片为这两个用户数组提供输入数据。
IDUM，整型数组卡
格式： IDUM I1，I2，…，In （1≤n≤50）
RDUM，实型数组卡
格式： RDUM R1，R2，…，Rn （1≤n≤50）
外围卡
这组卡片为用户提供方便，不影响MCNP的计算。
运行期间打印计数
MPLOT
打印控制
PRINT
粒子径迹输出
PTRAC
建立用户文件
FILES
调试信息
DBCN
微扰卡
PERT
丢失粒子
LOST
PRDMP
助记符
打印及转储周期
卡片类型
蒙特卡罗中心服务器信息
IP地址：
目前安装了MCNP4C，以后将陆续安装Egs4、Geant4、Fluka等蒙特卡罗程序。
中心电话： 62784552
联系人： 范佳锦、武祯
