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

蒙特卡罗方法在核技术中的应用 林谦 目录 第一章蒙特卡罗方法概述第二章随机数第三章由已知分布的随机抽样第四章蒙特卡罗方法解粒子输运问题 教材 蒙特卡罗方法在实验核物理中的应用许淑艳编著原子能出版社蒙特卡罗方法清华大学 参考书 蒙特卡罗方法及其在粒子输运问题中的应用裴鹿成张孝泽编著科学出版社蒙特卡罗方法徐钟济编著上海科学技术出版社 联系方式 电话83918电子邮件linqian 第一章蒙特卡罗方法概述 蒙特卡罗方法的基本思想蒙特卡罗方法的收敛性 误差蒙特卡罗方法的特点蒙特卡罗方法的主要应用范围作业 第一章蒙特卡罗方法概述 蒙特卡罗方法又称随机抽样技巧或统计试验方法 半个多世纪以来 由于科学技术的发展和电子计算机的发明 这种方法作为一种独立的方法被提出来 并首先在核武器的试验与研制中得到了应用 蒙特卡罗方法是一种计算方法 但与一般数值计算方法有很大区别 它是以概率统计理论为基础的一种方法 由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程 解决一些数值方法难以解决的问题 因而该方法的应用领域日趋广泛 蒙特卡罗方法的基本思想 二十世纪四十年代中期 由于科学技术的发展和电子计算机的发明 蒙特卡罗方法作为一种独立的方法被提出来 并首先在核武器的试验与研制中得到了应用 但其基本思想并非新颖 人们在生产实践和科学试验中就已发现 并加以利用 两个例子例1 蒲丰氏问题例2 射击问题 打靶游戏 基本思想计算机模拟试验过程 例1 蒲丰氏问题 为了求得圆周率 值 在十九世纪后期 有很多人作了这样的试验 将长为2l的一根针任意投到地面上 用针与一组相间距离为2a l a 的平行线相交的频率代替概率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 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 xi 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也是相互独立的 P R Halmos Measuretheory N Y VonNosrtand 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 0 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个 即 随机抽样的一般方法 加抽样方法减抽样方法乘抽样方法乘加抽样方法乘减抽样方法对称抽样方法积分抽样方法 加抽样方法 加抽样方法是对如下加分布给出的一种抽样方法 其中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 0 965 B 2 29 C 0 453 Emin 0 Emax 若采用乘减抽样方法 其抽样效率约为0 5 令相应的则从fa x 的抽样为从f1 x 的抽样为 参数 的确定 使1 A 0 且使H1 E 的上界M1最小 裂变中子谱的近似修正抽样方法为对于铀 235 m 0 8746 M 0 2678 0 5543 抽样效率E 0 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 中抽样 得到 前面已经介绍了 指数分布的直接抽样为 这不仅需要计算对数 而且由于要使用伪随机数 受精度的限制 该抽样值在小概率处即数值较大处呈现明显得离散性 下面介绍两种抽样方法可以避免这些问题 指数分布的抽样 所用随机数的平均个数N e2 e 1 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 确定粒子的初始状态 实际上就是要从中子源的空间位置 能量和方向分布中抽样 设源分布为则分别从各自的分布中抽样确定初始状态 对于平板情况 抽样得到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 则介质