蒙特卡洛应用讲义 完整版

蒙特卡罗方法 在核技术中的应用 林谦 目 录 第一章 蒙特卡罗方法概述 第二章 随机数 第三章 由已知分布的随机抽样 第四章 蒙特卡罗方法解粒子输运问题 教材 蒙特卡罗方法在实验核物理中的应用 许淑艳 编著 原子能出版社 蒙特卡罗方法 清华大学 参考书 蒙特卡罗方法及其在粒子输运问题中的应用 裴鹿成 张孝泽 编著 科学出版社 蒙特卡罗方法 徐钟济 编著 上海科学技术出版社 联系方式 电话 83918 电子邮件 第一章 蒙特卡罗方法概述 蒙特卡罗方法的基本思想 蒙特卡罗方法的收敛性，误差 蒙特卡罗方法的特点 蒙特卡罗方法的主要应用范围 作 业 第一章 蒙特卡罗方法概述 蒙特卡罗方法 又称 随机抽样技巧 或 统计试验方法 。半个多世纪以来，由于科学技术的发展和电子计算机的发明 ，这种方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。 1. 蒙特卡罗方法的基本思想 二十世纪四十年代中期，由于科学技术的发展和电子计算机的发明，蒙特卡罗方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。但其基本思想并非新颖，人们在生产实践和科学试验中就已发现，并加以利用。 两个例子 例 1. 蒲丰氏问题 例 2. 射击问题（打靶游戏） 基本思想 计算机模拟试验过程 例 1. 蒲丰氏问题 为了求得圆周率 值，在十九世纪后期，有很多人作了这样的试验：将长为 2l的一根针任意投到地面上，用针与一组相间距离为 2a（ l a）的平行线相交的频率代替概率 P，再利用准确的关系式： 求出 值 其中 为投计次数， n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。 alP2)(22 nNa laP l 一些人进行了实验，其结果列于下表 ： 实验者 年份 投计次数 的实验值 沃尔弗 (Wolf) 1850 5000 3.1596 斯密思 (Smith) 1855 3204 3.1553 福克斯 (Fox) 1894 1120 3.1419 拉查里尼(Lazzarini) 1901 3408 3.1415929 例 2. 射击问题（打靶游戏） 设 r表示射击运动员的弹着点到靶心的距离， (r)表示击中 r处相应的得分数（环数）， f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射击水平。该运动员的射击成绩为 用概率语言来说， 是随机变量 (r)的数学期望，即 )( rgEg 0 )()( drrfrgg 现假设该运动员进行了 次射击 ， 每次射击的弹着点依次为 r1 ， r2 ， ， rN ， 则 次得分 g(r1) ，g(r2)， ， g(rN)的算术平均值 代表了该运动员的成绩 。 换言之 ， 为积分 的估计值 ， 或近似值 。 在该例中，用 次试验所得成绩的算术平均值作为数学期望 的估计值（积分近似值）。 NiiN rgNg1)(1 基本思想 由以上两个例子可以看出 ， 当所求问题的解是某个事件的概率 ， 或者是某个随机变量的数学期望 ， 或者是与概率 、 数学期望有关的量时 ， 通过某种试验的方法 ， 得出该事件发生的频率 ， 或者该随机变量若干个具体观察值的算术平均值 ， 通过它得到问题的解 。这就是蒙特卡罗方法的基本思想 。 当随机变量的取值仅为 1或 0时 ， 它的数学期望就是某个事件的概率 。 或者说 ， 某种事件的概率也是随机变量 （ 仅取值为 1或 0） 的数学期望 。 因此，可以通俗地说，蒙特卡罗方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数 f(r)的随机变量 (r)的数学期望 通过某种试验 ， 得到 个观察值 r1， r2， ， rN（ 用概率语言来说 ， 从分布密度函数 f(r)中抽取 个子样 r1，r2， ， rN， ） ， 将相应的 个随机变量的值 g(r1)，g(r2)， ， g(rN)的算术平均值 作为积分的估计值（近似值）。 NiiN rgNg1)(1 0 )()( drrfrgg 为了得到具有一定精确度的近似解，所需试验的次数是很多的，通过人工方法作大量的试验相当困难，甚至是不可能的。因此，蒙特卡罗方法的基本思想虽然早已被人们提出，却很少被使用。本世纪四十年代以来，由于电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验交由计算机完成，使得蒙特卡罗方法得以广泛地应用，在现代化的科学技术中发挥应有的作用。 计算机模拟试验过程 计算机模拟试验过程 ， 就是将试验过程 （ 如投针 ，射击 ） 化为数学问题 ， 在计算机上实现 。 以上述两个问题为例 ， 分别加以说明 。 例 1. 蒲丰氏问题 例 2. 射击问题（打靶游戏） 由上面两个例题看出 ， 蒙特卡罗方法常以一个“ 概率模型 ” 为基础 ， 按照它所描述的过程 ， 使用由已知分布抽样的方法 ， 得到部分试验结果的观察值 ，求得问题的近似解 。 例蒲丰氏问题 设针投到地面上的位置可以用一组参数（ x,）来描述， x为针中心的坐标， 为针与平行线的夹角，如图所示。 任意投针，就是意味着 x与都是任意取的，但 x的范围限于 0， a，夹角 的范围限于 0， 。在此情况下，针与平行线相交的数学条件是 针在平行线间的位置 s in 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,)， 为 如果投针 次 ， 则 是针与平行线相交概率 的估计值。事实上， 于是有 其他当,0s i n,1),(lxxsNiiiN xsNs1),(1 aladxdd xdfxfxsPl 2)()(),(s i n0021NsalaPl 22 例射击问题 设射击运动员的弹着点分布为 用计算机作随机试验（射击）的方法为，选取一个随机数 ，按右边所列方法判断得到成绩。 这样，就进行了一次随机试验（射击），得到了一次成绩 (r)，作 次试验后，得到该运动员射击成绩的近似值 环数 7 8 9 10 概率 0.1 0.1 0.3 0.5 环中命环命中环命中环命中1095.082.071.0 NiiN rgNg1)(12. 蒙特卡罗方法的收敛性，误差 蒙特卡罗方法作为一种计算方法 ， 其收敛性与误差是普遍关心的一个重要问题 。 收敛性 误差 减小方差的各种技巧 效率 收敛性 由前面介绍可知，蒙特卡罗方法是由随机变量 X的简单子样 X1， X2， ， XN的算术平均值： 作为所求解的近似值。由大数定律可知， 如 X1， X2， ， XN独立同分布，且具有有限期望值（ E(X)0 。对该分布的直接抽样方法如下： !)(nePnxPnn n0ii1n0iii!i!,当 enX F例 3. 掷骰子点数的抽样 掷骰子点数 X=n的概率为： 选取随机数 ，如 则 在等概率的情况下，可使用如下更简单的方法： 其中表示取整数。 61)( nXP661 nn nX F 16 FX例 4. 碰撞核种类的确定 中子或光子在介质中发生碰撞时，如介质是由多种元素组成，需要确定碰撞核的种类。假定介质中每种核的宏观总截面分别为 1， 2， ， n，则中子或光子与每种核碰撞的概率分别为： 其中 t 1 2 n。碰撞核种类的确定方法为：产生一个随机数 ，如果 则中子或光子与第 I种核发生碰撞。 niPtii ,2,1 I1ii1I1ii PP例 5. 中子与核的反应类型的确定 假设中子与核的反应类型有如下几种 :弹性散射，非弹性散射，裂变，吸收，相应的反应截面分别为 el，in， f， a。则发生每一种反应类型的概率依次为 ： 其中反应总截面 t el in f a。 taatfftinintelel PPPP 反应类型的确定方法为：产生一个随机数 收吸裂变非弹性散射弹性散射 finelinelelPPPPPP2) 连续型分布的直接抽样方法 对于连续型分布 ， 如果分布函数 F(x) 的反函数 F 1(x)存在 ， 则直接抽样方法是 ： )(1 FX F例 6. 在 a， b上均匀分布的抽样 在 a， b上均匀分布的分布函数为： 则 bxbxaabaxaxxF当当当10)( )( abaX F例 7. 分布 分布为连续型分布，作为它的一个特例是： 其分布函数为： 则 FX10,2)()( 20 xxtdtdttfxF xx10,2)( xxxf例 8. 指数分布 指数分布为连续型分布，其一般形式如下： 其分布函数为： 则 因为 1 也是随机数，可将上式简化为 0,1)()( 0 xedteadttfxF axx atx0,)( xeaxf ax)1ln (1 aX Fln1aX F 连续性分布函数的直接抽样方法对于分布函数的反函数存在且容易实现的情况 ， 使用起来是很方便的 。但是对于以下几种情况 ， 直接抽样法是不合适的 。 分布函数无法用解析形式给出 ， 因而其反函数也无法给出 。 分布函数可以给出其解析形式 ， 但是反函数给不出来 。 分布函数即使能够给出反函数 ， 但运算量很大 。 下面叙述的挑选抽样方法是克服这些困难的比较好的方法 。 3. 挑选抽样方法 为了实现从己知分布密度函数 f(x)抽样 ， 选取与f(x)取值范围相同的分布密度函数 h(x)， 如果 则挑选抽样方法为： )()(s u pxhxfMxhfhhXXXhMXf)()( 即从 h(x)中抽样 xh， 以 的概率接受它 。 下面证明 xf 服从分布密度函数 f(x)。 证明：对于任意 x )()(hhxhM xf)()()()(,)()()(hhhhhhhhfXhMXfPXhMXfdxxXxPXhMXfdxxXxPdxxXxP dxxfdXXfdXXfdXXhXhMXfdXXhXhMXfddXXhddXXhhhdxxxhhhhhhdxxxhhhhXhMXfhhdxxxXhMXfhhhhhh)()()()()()()()()()()()()(0)()(0 使用挑选抽样方法时 ， 要注意以下两点：选取 h(x)时要使得 h(x)容易抽样且 M的值要尽量小 。 因为 M小能提高抽样效率 。 抽样效率是指在挑选抽样方法中进行挑选时被选中的概率 。 按此定义 ， 该方法的抽样效率 E为： 所以 ， M越小 ， 抽样效率越高 。 MdXXhXhMXfXhMXfPEhhhhhh1)()()()()( 当 f(x) 在 0， 1 上定义时 ， 取 h(x)=1， Xh= ， 此时挑选抽样方法为 )(s u p10xfMx fXMf )(例 9. 圆内均匀分布抽 样 令圆半径为 R0，点到圆心的距离为 r，则 r的分布密度函数为 分布函数为 容易知道，该分布的直接抽样方法是 其它当002)( 020RrRrrf202)(RrrF 0Rr f 由于开方运算在计算机上很费时间 ， 该方法不是好方法 。 下面使用挑选抽样方法：取 则抽样框图为 00022)()(1)( RrMRrrhrfRrh h， 2021Rr f 显然 ， 没有必要舍弃 1 2的情况 ， 此时 ， 只需取 就可以了 ， 亦即 另一方面，也可证明 与 具有相同的分布 。 10 Rr f),m ax ( 210 Rr f ),m a x( 21 2)( rrF 4. 复合抽样方法 在实际问题中 ， 经常有这样的随机变量 ， 它服从的分布与一个参数有关 ， 而该参数也是一个服从确定分布的随机变量 ， 称这样的随机变量服从复合分布 。例如 ， 分布密度函数 是一个复合分布 。 其中 Pn0， n=1， 2， ， 且 fn(x)为与参数 n有关的分布密度函数 ， n=1， 2， ， 参数 n服从如下分布 1)()(nnn xfPxf11nnPynnPyF )( 复合分布的一般形式为： 其中 f2(x/y)表示与参数 y有关的条件分布密度函数 ， F1(y)表示分布函数 。 复合分布的抽样方法为 :首先由分布函数 F1(y) 或分布密度函数 f1(y)中抽样 YF1或 Yf1，然后再由分布密度函数 f2(x/ YF1)中抽样确定 Xf2 (x/YF) 证明： 所以 ， Xf所服从的分布为 f (x)。 )()()( 12 ydFyxfxf )/( 12 FYxff XX dxxfYd x dFYxfdxxXxpdxxXxpFYxff)()()()()(12)/( 12例 10. 指数函数分布的抽 样 指数函数分布的一般形式为： 引入如下两个分布密度函数： 其它当00)( 1xdyyenxE nxyn 其它当其它当00)(01)(211xeyyxfyynyfxyn 则 使用复合抽样方法 ， 首先从 f1(y)中抽取 y 再由 f2(x/ YF1)中抽取 x 1 12 )()()( dyyfyxfxE n),m a x (11211nnfY 1211ln),m a x (ln1nnfnfYX5. 复合挑选抽样方法 考虑另一种形式的复合分布如下： 其中 0H(x,y)M， f2(x/y)表示与参数 y有关的条件分布密度函数 ， F1(y)表示分布函数 。 抽样方法如下： )()(),()( 12 ydFyxfyxHxf )/()/(12112),(FFYxffFYxfXXMYXH 证明： 抽样效率为： E=1/M dxxfdxydFyxfyxHyd x d FyxfMyxHyd x d FyxfMyxHdyd x d Fyxfdyd x d FyxfMYXHPMYXHdxxXxPMYXHdxxXxPdxxXxPdxxxMyxHdxxxMyxHFfFffFfff)()()(),()()(),()()(),()()()()(),(),(,),()(121212),(012),(01212122122 为了实现某个复杂的随机变量 y 的抽样，将其表示成若干个简单的随机变量 x1， x2， ， xn 的函数 得到 x1， x2， ， xn 的抽样后，即可确定 y 的抽样，这种方法叫作替换法抽样。即 6. 替换抽样方法 ),( 21 nxxxgy ),( 21 nf XXXgY 例 11. 散射方位角余弦分布的抽 样 散射方位角 在 0,2上均匀分布，则其正弦和余弦 sin和 cos服从如下分布： 直接抽样方法为： 其它当011111)( 2xxxf 2co sco s2s ins in 令 =2， 则 在 0, 上均匀分布 ， 作变换 其中 01， 0， 则 (x,y) 表示上半个单位圆内的点 。 如果 (x,y) 在上半个单位圆内均匀分布 ， 则 在 0, 上均匀分布 ， 由于 s inc o syx2222s inc o syxyyxx222222222co ss i n22s i ns i ns i nco s2co sco syxxyyxyx 因此抽样 sin和 cos的问题就变成在上半个单位圆内均匀抽样 (x,y) 的问题 。 为获得上半个单位圆内 的均匀点 ， 采用挑选法 ， 在 上半个单位圆的外切矩形内 均匀投点 （ 如图 ） 。 舍弃圆外的点 ， 余下的就是所要求的点 。 抽样方法为： 抽样效率 E=/40.785 21 yx2221212221222122212s i n,c o s1 为实现散射方位角余弦分布抽样 ， 最重要的是在上半个单位圆内产生均匀分布点 。 下面这种方法 ， 首先在单位圆的半个外切正六边形内产生均匀分布点 ，如图所示 。 于是便有了抽样效率更高的抽样方法： 抽样效率 222121222122212221221121332s i n,33c o s131,123 906.032 E例 12. 正态分布的抽 样 标准正态分布密度函数为： 引入一个与标准正态随机变量 X独立同分布的随机变量 Y，则（ X,Y）的联合分布密度为： 作变换 2221)( xexf 2)( 2221),( yxeyxf s inco syx 则 （ ,） 的联合分布密度函数为： 由此可知 ， 与 相互独立 ， 其分布密度函数分别为 分别抽取 ， ： 222),( ef 21)()(2212 ff e212ln2 从而得到一对服从标准正态分布的随机变量 X和 Y： 对于一般的正态分布密度函数 N(,2) 的抽样 ， 其抽样结果为： ffffYYXX)2s i n(ln2)2c os (ln22121ffYX例 13. 分布的抽 样 分布密度函数的一般形式为： 其中 n， k为整数。为了实现 分布的抽样，将其看作一组简单的相互独立随机变量的函数，通过这些简单随机变量的抽样，实现 分布的抽样。设 x1， x2， ，xn 为一组相互独立、具有相同分布 F(x) 的随机变量，k为 x1， x2， ， xn 按大小顺序排列后的第 k个，记为： 10)1()!()!1(!)( 1 xxxknknxf knk),( 21 nkk xxxR 则 k的分布函数为： 当 F(x)=x 时 ， 不难验证 ， k的分布密度函数为 分布 。 因此 ， 分布的抽样可用如下方法实现： 选取 n个随机数 ， 按大小顺序排列后取第 k个 ， 即 ininkiin xFxFCxF k )(1)()(ininkiin xxCxF k )1()(),( 21 nkf RX 7. 随机抽样的一般方法 加抽样方法 减抽样方法 乘抽样方法 乘加抽样方法 乘减抽样方法 对称抽样方法 积分抽样方法 1) 加抽样方法 加抽样方法是对如下加分布给出的一种抽样方法： 其中 Pn0， ，且 fn(x)为与参数 n有关的分布密度函数， n=1， 2， 。 由复合分布抽样方法可知，加分布的抽样方法为 :首先抽样确定 n，然后由 fn(x)中抽样 x，即： 1)()(nnn xfPxf11nnPn1nn1n1nn PP,当 nffXX例 14. 多项式分布抽 样 多项式分布密度函数的一般形式为： 将 f(x) 改写成如下形式： 则该分布的抽样方法为： 0)(iii xaxf00)()1(1)(iiiiii xfPxiiaxf n0ii1n0ii11 PP),m ax (当 nfX 例 15. 球壳内均匀分布抽 样 设球壳内半径为 R0，外半径为 R1，点到球心的距离为 r，则 r的分布密度函数为 分布函数为 该分布的直接抽样方法是 其它当03)( 1030312RrRRRrrf3031303)(RRRrrF 31303031 )( RRRr f 为避免开立方根运算 ， 作变换： 则 x 0,1， 其分布密度函数为： 其中 001 )( RxRRr 132)(33)()(200102201 RxRRRxRRxf211020 RRRR 则 x及 r的抽样方法为： 001432322101201)(),m a x (),m a x (33RXRRrXXXRRRfffff 2) 减抽样方法 减抽样方法是对如下形式的分布密度所给出的一种抽样方法： 其中 A1、 A2为非负实数， f1(x) 、 f2(x)均为分布密度函数。 减抽样方法分为以下两种形式： )()()( 2211 xfAxfAxf 以上两种形式的抽样方法，究竟选择哪种好，要看 f1(x) 、 f2(x)哪一个容易抽样，如相差不多，选用第一种方法抽样效率高。 （ 1） 将 f (x)表示为 令 m表示 f2(x) f1(x)的下界 ， 使用挑选法 ， 从 f1(x)中抽取 Xf1 抽样效率为： )()()()(12211 xfxfAAxfxf 111)()(12212211ffffXXXfXfmAAAmAAA211mAAE （ 2） 将 f (x)表示为 使用挑选法 ， 从 f2(x)中抽取 Xf2 抽样效率为： 22112 )()()()( AxfxfAxfxfmEmAAmE 21 22221221211)()(ffffXXmAAmAXfXfmAAmA例 16. 分布抽样 分布的一个特例： 取 A1 2， A2 1， f1(x) 1， f2(x) 2x，此时 m 0，则根据第一种形式的减抽样方法，有 或 10),1(2)( xxxf 221 1fX 212 1fX 由于 1 1可用 1代替 ， 该抽样方法可简化为： 对于 2 1的情况 ， 可取 Xf 1 ， 因此 与 分布的推论相同 。 212fX ),m in ( 21 fX 如下形式的分布称为乘分布： 其中 H(x)为非负函数， f1(x)为任意分布密度函数。 令 M为 H(x)的上界，乘抽样方法如下： 抽样效率为： 3) 乘抽样方法 )()()( 1 xfxHxf ME111)(fffXXMXH 例 17. 倒数分布抽样 倒数分布密度函数为： 其直接抽样方法为： 下面采用乘抽样方法，考虑如下分布族： 其中 i = 1， 2， ，该分布的直接抽样方法为： axxaxf 1,1ln1)( af eaX lniif aXi 1)1(1 axxxaixfiii 1,)1(1)( 11 利用这一分布族 ， 将倒数分布 f(x) 表示成 : 其中 ， 乘法分布的抽样方法如下 : 该分布的抽样效率为： )()()( xfxHxf i,1)(,ln)1(,ln)1()(1111iiiixMxHaaiMxaaixH iifiaXa1)1(11)1(21211 )1(ln1 iaiaE例 18. 麦克斯韦 (Maxwell)分布抽样 麦克斯韦分布密度函数的一般形式为： 使用乘抽样方法，令 该分布的直接抽样方法为： 0,2)(23 xexxf x2ln231fX0,32)( 321 xexf x 此时 则麦克斯韦分布的抽样方法为 : 该分布的抽样效率为： eMxexxHx 2270,3)( 3121 22221ln23lnfXe 7 9 5.0272 eE 在实际问题中，经常会遇到如下形式的分布： 其中 Hn(x)为非负函数， fn(x) 为任意分布密度函数，n=1， 2， 。不失一般性，只考虑 n=2的情况： 将 f(x) 改写成如下的加分布形式： 4) 乘加抽样方法 1)()()(nnn xfxHxf)()()()()( 2211 xfxHxfxHxf )()()()()()()(*22*1122221111xfPxfPxfPxHPxfPxHPxf 其中 )()()()()()()()()()(222*2111*1222111xfPxHxfxfPxHxfdxxfxHPdxxfxHP 乘加抽样方法为： 该方法的抽样效率为： 22222)(fffXXMXH 11112)(fffXXMXH 11 P 2221212221111 MPMPMPPMPPE 这种方法需要知道 P1的值 (P2=1 P1） ， 这对有些分布是很困难的 。 下面的方法可以不用计算 P1 ： 对于任意小于 1的正数 P1 ， 令 P2=1 P1 ； ynnPyFyxfyxfyxfyPxHyPxHyxH)(2),(21),()(2)(21,)(),(12122211),m in (),m a x (22112211MPMPEPMPMM则采用复合挑选抽样方法，有： 当取 时 ， 抽样效率最高 这时 ， 乘加抽样方法为： 22222)(fffXXMXH 11112)(fffXXMXH 2111 MMM 2121MME21222111 MMMPMMMP 由于 可知第一种方法比第二种方法的抽样效率高 。 0)()()(2)()()()()(121212211221212121222121222121212221212221212221212122211212122122212121MMMMPMPMMMMMPPMMPMPMMMMMMMPPMMPMPMMMMMMMPMMMPMMMMMMPMPEE例 19. 光子散射后能量分布的抽样 令光子散射前后的能量分别为 和 （以 m0c2 为单位， m0为电子静止质量， c 为光速）， ， 则 x 的分布密度函数为： 该分布即为光子散射能量分布，它是由著名的 KlinNishina 公式确定的。其中 K() 为归一因子： 211,1111)(1)(322 xxxxxxKxf x22 )21(21421)21l n ()1(21)( K 把光子散射能量分布改写成如下形式： 在 1, 1+2 上定义如下函数 : 3222)1(111)(1)(xxxxKxf32221221)1()(2)(11)21)(2)(21)(1221)(xxKxHxKxHxfxxf 则有 使用乘加抽样方法： )()()()()( 2211 xfxHxfxHxf 221122)(278)()21)(4)(21212121MKxHMKxHXXff 光子散射能量分布的抽样方法为： 该方法的抽样效率为： 222211132321232)1(427212942711212121fffffffffXXXXXXXXX )294(4)()21(27121 KMME 乘减分布的形式为： 其中 H1(x) 、 H2(x)为非负函数， f1(x)、 f2(x) 为任意分布密度函数。 与减抽样方法类似，乘减分布的抽样方法也分为两种。 5) 乘减抽样方法 )()()()()( 2211 xfxHxfxHxf （ 1） 将 f (x) 表示为 令 H1(x)的上界为 M1， 的下界为 m， 使用乘抽 样方法得到如下乘减抽样方法： )()()()(1)()()(112211 xfxHxfxHxHxfxf 11111)()()()()1(112211ffffffXXXfXfXHXHmM)()()()(1122xfxHxfxH （ 2） 将 f (x) 表示为 令 H2(x)的上界为 M2， 使用乘抽样方法 ， 得到另一种乘减抽样方法： 1)()()()()()()(221122 xfxHxfxHxHxfxf 22222)()()()()1(22112ffffffXXXHXfXfXHmMm例 20. 裂变中子谱分布抽样 裂变中子谱分布的一般形式为： 其中 A， B， C， Emin， Emax 均为与元素有关的量。令 其中 为归一因子， 为任意参数。 m a xm i n,sh)( EEEBEeCEfAE m a xm i n21 ,)()( EEEeEfEfE 相应的 H1(E)， H2(E) 为： 于是裂变中子谱分布可以表示成乘减分布形式 : 容易确定 H1(E) 的上界为： 为提高抽样效率 ， 应取 使得 M1 达到最小 ， 此时 1141e xp21 ABCMBEEACEHBEEACEH11e x p2)(11e x p2)(21)()()()()( 2211 EfEHEfEHEf 116181ABABA 取 m 0， 令 则裂变中子谱分布的抽样方法为 : 抽样效率 4,11 BA1111221m i n)( l nex plnfffffEEBEEEE eCME 211 对称分布的一般形式为： 其中 f1(x) 为任意分布密度函数，满足偶函数对称条件，H(x) 为任意奇函数，即对任意 x满足： 对称分布的抽样方法如下：取 =2 1 6) 对称抽样方法 )()()( 1 xHxfxf )()()()( 11xHxHxfxf1ff XX 1ff XX )()(111 ffXfXH 证明： 因为 =2 1， x 相当于 ， 因此 xffxfffxfffxfffxfffxfffxffffxfffffffffffffffffdXXfdXXHXfdXXHXfdXXHXfdXXHXfdXXHXfdXXfXfXHdXXfXfXHXfXHxXPXfXHxXPXfXHxXPXfXHxXPxXP11111111111111111111111111111111111)()()()()(21)()(21)()(21)()(21)()()(121)()()(121)()(121,)()(121,)()(121,)()(121,)(1111111111111 x121例 21. 质心系各向同性散射角余弦分布抽样 在质心系各向同性散射的假设下，为得到实验室系散射角余弦，需首先抽样确定质心条散射角余弦： 再利用下面转换公式 : 得到实验室系散射角余弦 L。其中 A为碰撞核质量，C、 L 分别为质心系和实验室系散射角。 12co s11,21)(CCCCfCCLLAAA211c os2 为避免开方运算 ， 可以使用对称分布抽样 。 根据转换公式可得： 依照质心系散射各向同性的假定 ， 可得到实验室系散射角余弦 L 的分布如下： 该密度函数中的第一项为偶函数 ， 第二项为奇函数 ，因而是对称分布 。 其中 11,212121)(2222 LLLLLAAAf 11,12121)(22221 LLLLAAAf 111 222 LLLC AA 从 f1(L) 的抽样可使用挑选法 然后再以 的概率决定接受或取负值 。 上述公式涉及开方运算 ， 需要进一步简化 。 22221212LLLAA1221221221121LAAAA 注意以下事实：对于任意 0a1 令 则上述挑选抽样中的挑选条件简化为： 另一方面 ， 在 即 的条件下 ， 2/a 在 1, 1 上均匀分布 ， 故可令 2/a， 则最终决定取正负值的条件简化为： aaPaP )()( 2222 AAAAa 1121 2212212 21222122221)21(1 AAAA222 a 11 2 a12221 AA 于是 ， 得到质心系各向同性散射角余弦分布的抽样方法为： 112221222122221211)21(1 LLAAAAAA 如下形式的分布密度函数 称为积分分布密度函数，其中 f0(x， y) 为任意二维分布密度函数， H(x)为任意函数。该分布密度函数的抽样方法为： 7) 积分抽样方法 )(0)(0),(),()(xHxHd xd yyxfdyyxfxf000)(ffffXXXHY 证明：对于任意 x dxxfd x d yyxfd x d yyxfXHYPXHYdxxXxPXHYdxxXxPdxxXxPxHxHfffffffff)(),(),()()(,)()()(0)(000000000 例 22. 各向同性散射方向的抽样 为了确定各向同性散射方向 ，根据公式： 对于各向同性散射， cos在 1, 1上均匀分布， 在 0, 2上均匀分布。由于 直接抽样需要计算三角函数和开方。 kji wvu co ss ins inco ss inwvu22 1c o s1s i n w 定义两个随机变量： 可以证明 ， 当 时 ， 随机变量 x 和 y 服从如下分布 : 定义区域为： xxAyxxyAyxf12,12m i n0,111221),(22232222212322222123222221AAyAAAAx12322 则 w cos 的分布可以用上述分布表示成积分分布的形式： 令 ， 则属于上述积分限内的 y 一定满足 条件 。 3 163A12322 21),(),()(31213121)()( xxd x d yyxfdyyxfxf 各向同性散射方向的抽样方法为： 抽样效率为： 2322222123222221232222213123222221211223222221c o s2s i ns i n2c o ss i n)(AAAAwAAAvAAAuAA555.012 2 AE 8. 随机抽样的其它方法 偏倚抽样方法 近似抽样方法 近似 -修正抽样方法 多维分布抽样方法 指数分布的抽样 使用蒙特卡罗方法计算积分 时，可考虑将积分 I改写为 其中 f *(x) 为一个与 f (x) 有相同定义域的新的分布密度函数。于是可以这样计算积分 I： 这里 Xi 是从 f *(x) 中抽取的第 i 个子样。 1) 偏移抽样方法 dxxfxgdxxfxf xfxgI )()()()( )()( * dxxfxgI )()(NiiiiNi iiNiiN XgXWNXgXfXfNXgNI11*1* )()(1)()()(1)(1 由此可以看出 ， 原来由 f (x) 抽样 ， 现改为由另一个分布密度函数 f *(x) 抽样 ， 并附带一个权重纠偏因子 这种方法称为偏倚抽样方法 。 从 f (x) 中抽取的 Xf ， 满足 而对于偏倚抽样 ， 有 一般情况下 ， Xf 是具有分布 f (x) 总体的简单子样的个体 ， 只代表一个 。 Xf* 是具有分布 f *(x) 总体的简单子样的个体 ， 但不代表一个 ， 而是代表 W(Xf*) 个 ，这时 Xf*是带权 W(Xf*)服从分布 f (x) 。 )()()( * xfxfxW dxxfdxdXxP f )()( dxxfdxxfxWdxdXxPXW ff )()()()()( * 在实际问题中，分布密度函数的形式有时是非常复杂的，有些甚至不能用解析形式给出，只能用数据或曲线形式给出。如中子散射角余弦分布多数是以曲线形式给出的。对于这样的分布，需要用近似分布密度函数代替原来的分布密度函数， 用近似分布密度函数的抽样代替原分布密度函数