第三章 由已知分布的随机抽样 随机抽样及其特点 直接抽样方法 挑选抽样方法 复合抽样方法 复合挑选抽样方法 替换抽样方法 随机抽样的一般方法 第三章 由已知分布的随机抽样 随机抽样及其特点 直接抽样方法 挑选抽样方法 复合抽样方法 复合挑选抽样方法 替换抽样方法 随机抽样的一般方法 随机抽样的其它方法 作 业

第三章 由已知分布的随机抽样 本章叙述由己知分布抽样的各主要方法，并给出在粒子输运问题中经常用到的具体实例。

随机抽样及其特点 由巳知分布的随机抽样指的是由己知分布的总体中抽取简单子样。随机数序列是由单位均匀分布的总体中抽取的简单子样，属于一种特殊的由已知分布的随机抽样问题。本章所叙述的由任意已知分布中抽取简单子样，是在假设随机数为已知量的前提下，使用严格的数学方法产生的。 为方便起见，用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, Measure theory, N.Y.Von 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为：

当 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）的联合分布密度为： 例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. 倒数分布抽样 倒数分布密度函数为： 其直接抽样方法为： 下面采用乘抽样方法，考虑如下分布族： 例17. 倒数分布抽样 倒数分布密度函数为： 其直接抽样方法为： 下面采用乘抽样方法，考虑如下分布族： 其中 i = 1，2，…，该分布的直接抽样方法为：

利用这一分布族，将倒数分布 f(x) 表示成: 其中， 乘法分布的抽样方法如下: 该分布的抽样效率为： ＞ ≤

例18. 麦克斯韦(Maxwell)分布抽样 麦克斯韦分布密度函数的一般形式为： 使用乘抽样方法，令 该分布的直接抽样方法为：

此时 则麦克斯韦分布的抽样方法为: 该分布的抽样效率为： ＞ ≤

乘加抽样方法 在实际问题中，经常会遇到如下形式的分布： 其中Hn(x)为非负函数，fn(x) 为任意分布密度函数，n=1，2，…。不失一般性，只考虑n=2的情况： 将 f(x) 改写成如下的加分布形式：

其中

乘加抽样方法为： 该方法的抽样效率为： > ≤

这种方法需要知道P1的值(P2=1－P1），这对有些分布是很困难的。下面的方法可以不用计算P1 ： 则采用复合挑选抽样方法，有：

当取 时，抽样效率最高 这时，乘加抽样方法为： > ≤

由于 可知第一种方法比第二种方法的抽样效率高。

例19. 光子散射后能量分布的抽样 令光子散射前后的能量分别为 和 （以 m0c2 为单位，m0为电子静止质量，c 为光速）， ， 例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 均为与元素有关的量。令 例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. 质心系各向同性散射角余弦分布抽样 在质心系各向同性散射的假设下，为得到实验室系散射角余弦，需首先抽样确定质心条散射角余弦： 例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 均为与元素有关的量。 例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 。 二维分布密度函数的抽样方法是: 首先由 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 Y 所用随机数的平均个数 N＝e2 / ( e－1)≈4.3

方法二 ＞ ≤ N Y

作 业 光子散射后能量分布的抽样 把光子散射能量分布改写成如下形式进行抽样：

在［1, 1+2α］上定义如下函数:

＞ ≤