8 随机变量生成

本章讨论从各种常用的连续分布与离散分布中抽样的方法。此前的讨论和示例表明，概率分布在刻画那些一般无法预测或存在不确定性的建模活动时很有用。例如，队列处的到达间隔时间（interarrival time）和服务时间（service time）以及产品需求，在相当程度上常常是无法预测的。通常，这类变量被建模为服从某个给定概率分布的随机变量（random variable），并且已有标准的统计程序可用于估计所假设分布的参数，以及检验所假定统计模型的有效性。

在本章中，我们假定一个分布已经完全指定，并寻求从该分布生成样本的方法，以用作仿真模型（simulation model）的输入。本章的目的是解释并举例说明若干广泛使用的随机变量生成（random variate generation）技术，而不是要对最高效的技术作一份最新进展的综述。在实践中，大多数仿真建模人员会使用编程库中已有的例程，或者仿真语言内置的例程。然而，有些编程语言并未为所有常用分布提供内置例程，有些计算机配置也没有随机变量生成库；在这些情况下，建模者必须自行构造一个可接受的例程。尽管发生这种情况的可能性很小，但了解随机变量生成究竟是如何进行的，仍然是值得的。

本章讨论逆变换法（inverse-transform technique），并较为简要地讨论接受-拒绝法（acceptance-rejection technique）和特殊性质。另一种技术——组合法（composition method）——由 Devroye [1986]、Dagpunar [1988] 和 Law [2007] 讨论。本章所有技术都假定可以方便地获得 [0, 1] 上的均匀随机数源 $R_{1} , R_{2} , \ldots$，其中每个 $R_{i}$ 的 pdf 为

cdf 为

在本章中，R 与 $R_{1} , R_{2} , \ldots$ 代表在 [0, 1] 上均匀分布的随机数，它们由随机数生成（random-number generation）技术产生，或取自随机数表，例如附录中的表 A.1。本章中的许多算法已被译成 Visual Basic for Applications，可在 www.bcnn.net 上的电子表格 RandomNumberTools.xls 中找到。

1 逆变换法

逆变换法可用于从指数分布、均匀分布、Weibull 分布和三角分布中抽样，也可用于从经验分布中抽样。此外，它还是从大量离散分布中抽样的基本原理。本节先以指数分布为例详细说明该方法，然后再将其应用于其他分布。就计算而言，它是最直接的方法，但并不总是最高效的方法。

1.1 指数分布

指数分布的概率密度函数（pdf）为

累积分布函数（cdf）为

参数 λ 可以解释为单位时间内发生次数的均值。例如，若到达间隔时间 $X_{1} , X_{2} , X_{3} , . . .$ 服从速率为 λ 的指数分布，那么 λ 就可以解释为单位时间内到达次数的均值，即到达率（arrival rate）。注意，对任意 i，

因此 $1 / \lambda$ 就是平均到达间隔时间。本节的目标是构造一个生成值 $X_{1} , X_{2} , X_{3} , . . .$ 的过程，使其服从指数分布。

原则上，逆变换法可以用于任何分布，但只有当 cdf（$F ( x )$）的形式极为简单、其逆（$F^{- 1}$）易于计算时，该方法才最为有用。（记号 $F^{- 1}$ 表示方程 $r = F ( {\boldsymbol {x}} )$ 关于 x 的解；它并不表示 $1 / F$ 。）下面给出逆变换法的一个分步过程，此处以指数分布为例加以说明：

步骤 1. 计算所求随机变量 $X$ 的 cdf。

对于指数分布，cdf 为 $F ( x ) = 1 - e^{- \lambda x} , x \geq 0$ 。

步骤 2. 在 X 的取值范围内令 $F ( X ) = R$。

对于指数分布，它变为在范围 $X \geq 0$ 上的 $1 - e^{- \lambda X} = R$ 。X 是一个随机变量（在本例中服从指数分布），因此 $1 - e^{- \lambda X}$ 也是一个随机变量，这里称之为 R。如后文将要说明的，R 在区间 [0, 1] 上服从均匀分布。

步骤 3. 将方程 $F ( X ) = R$ 就 R 解出 X。

对于指数分布，求解过程如下：

式 (1) 称为指数分布的随机变量生成器（random-variate generator）。一般地，式 (1) 写作 $X = F^{- 1} ( R )$。要生成一串取值，则通过步骤 4 完成。

步骤 4. 生成（按需要）均匀随机数 $R_{1} , R_{2} , R_{3} , . . .$ ，并按下式计算所求的随机变量

对于指数情形，由式 (1) 得 $F^{- 1} ( R ) = ( - 1 / \lambda ) \ln ( 1 - R )$ ，因此

其中 $i = 1 , 2 , 3 , \dots$ 。式 (2) 中通常采用的一项简化是把 $1 - R_{i}$ 替换为 $R_{i}$ ，从而得到

这一替代做法之所以合理，是因为 $R_{i}$ 与 $1 - R_{i}$ 都在 [0, 1] 上服从均匀分布。

例 1

表 1 给出了取自表 A.1 的一串随机数，以及由式 (2)（取 $\lambda = 1$）计算得到的指数随机变量 $X_{i} ,$ 。图 1(a) 是取自均匀分布的 200 个值 $R_{1} , R_{2} , \ldots , R_{200}$ 的直方图，图 1(b) 是由式 (2) 计算得到的 200 个值 $X_{1} , X_{2} , \ldots , X_{200}$ 的直方图。把这些经验直方图与图 1(c) 和 (d) 中的理论密度函数作比较。如图所示，直方图是底层密度函数的一种估计。

表 1 给定随机数 $R_{i}$ 时，均值为 1 的指数随机变量 $X_{i}$ 的生成

图 2 给出了逆变换法的图解解释。图中所示的 cdf 为 $F ( x ) = 1 - e^{- x}$ ，这是一个速率为 $\lambda = 1$ 的指数分布。要生成一个具有 cdf $F ( x )$ 的值 $X_{1}$ ，先在 0 与 1 之间生成一个随机数 $R_{1}$ ，然后从 $R_{1}$ 向 cdf 的图形作一条水平线，再向下作一条垂直线到 x 轴，从而得到所求结果 $X_{1}$ 。注意 $R_{1}$ 与 $X_{1}$ 之间的互逆关系，即

以及

一般地，该关系写作

以及

由这一过程生成的随机变量 $X_{1}$ 为什么具有所期望的分布呢？取一个值 $x_{0}$ 并计算累积概率

要理解式 (4) 中的第一个等号，请参见图 2，其中固定的数 $x_{0}$ 和 $F ( x_{0} )$ 画在各自的坐标轴上。可以看出，$X_{1} \leq x_{0}$ 当且仅当 $R_{1} \leq F ( x_{0} )$ 。由于 $0 \leq F ( x_{0} ) \leq 1$，式 (4) 中的第二个等号立即由 $R_{1}$ 在 [0, 1] 上均匀分布这一事实得到。式 (4) 表明 $X_{1}$ 的 cdf 就是 $F$ ；因此，$X_{1}$ 具有所期望的分布。

图 1 (a) 200 个均匀随机数的经验直方图；(b) 200 个指数随机变量的经验直方图；(c) [0, 1] 上的理论均匀密度；(d) 均值为 1 的理论指数密度。

1.2 均匀分布

考虑一个在区间 [a, b] 上均匀分布的随机变量 X。生成 X 的一个合理设想是

图 2 逆变换法的图解。

（请记住，R 始终是 [0, 1] 上的随机数。）X 的 pdf 为

式 (5) 的推导沿用第 1.1 节的步骤 1 至步骤 3：

步骤 1. cdf 为

步骤 2. 令 $F ( X ) = ( X - a ) / ( b - a ) = R$ 。

步骤 3. 就 R 解出 X 得 $X = a + ( b - a ) R ,$ ，与式 (5) 一致。

1.3 Weibull 分布

Weibull 分布是作为机器或电子元件失效时间的模型而引入的。当位置参数 ν 取 0 时：

其中 $\alpha > 0$ 和 $\beta > 0$ 分别是该分布的尺度参数（scale parameter）和形状参数（shape parameter）。要生成一个 Weibull 随机变量，请沿用第 1.1 节的步骤 1 至步骤 3：

步骤 1. cdf 为 $F ( X ) = 1 - e ^ {- (x / \alpha) ^ {\beta}} , x \geq 0$ 。

步骤 2. 令 $F ( X ) = 1 - e ^ {- ( X / \alpha) ^ {\beta}} = R .$ 。

步骤 3. 就 R 解出 X 得

式 (6) 的推导留给读者作为习题 10。比较式 (6) 与式 (1) 可以看出，若 X 是 Weibull 随机变量，则 $X^{\beta}$ 是均值为 $\alpha^{\beta}$ 的指数随机变量。反过来，若 Y 是均值为 $\mu ,$ 的指数随机变量，则 $Y^{1 / \beta}$ 是形状参数为 $\beta$ 、尺度参数为 $\alpha = \mu^{1 / \beta}$ 的 Weibull 随机变量。

1.4 三角分布

考虑一个 pdf 为

的随机变量 X，如图 3 所示。这个分布称为三角分布（triangular distribution），其端点为 $( 0 , 2 )$ ，众数在 1 处。其 cdf 为

对于 $0 \leq X \leq 1$，

图 3 某一特定三角分布的密度函数。

而对于 $1 \leq X \leq 2$，

在式 $( 7 )$ 中，$0 \leq X \leq 1$ 意味着 $\begin{array} {r} {0 \leq R \leq \frac {1} {2}} \end{array}$ ，此时 $X = {\sqrt {2 R}}$ 。由式 (8)，$1 \leq X \leq 2$ 意味着 $\begin{array} {r} {\frac {1} {2} \leq R \leq 1} \end{array}$ ，此时 $X = 2 - \sqrt {2 ( 1 - R )}$ 。于是，X 由下式生成：

习题 2、3 和 4 让读者练习处理其他三角分布。注意，若随机变量 X 的 pdf 和 cdf 是分段给出的（即在 X 取值范围的不同部分需要不同的公式），那么应用逆变换法生成 X 时，将得到在 R 取值范围的不同部分上各不相同的公式，如式 (9) 所示。

1.5 经验连续分布

若建模者未能找到一个能为输入数据提供良好模型的理论分布，那么就有必要使用数据的经验分布（empirical distribution）。一种可能做法是直接对观测到的数据本身重抽样。这称为使用经验分布，当已知输入过程只取有限多个值时，这样做尤其合理。本章第 1.7 节给出了这类情形的一个例子，以及生成随机输入的一种方法。

另一方面，若数据被认为取自一个连续取值的输入过程，那么在观测数据点之间作插值以填补空缺就是合理的做法。本节描述一种定义连续经验分布并从中生成数据的方法。

例 2

已经收集到五组消防队响应时间（以分钟计）的观测值，用以开展一项仿真研究，考察可能的备选人员配置与班组排班政策。数据为

在收集更多数据之前，希望先建立一个初步的仿真模型，其中的响应时间分布基于这五组观测值。因此，需要一种从响应时间分布中生成随机变量的方法。初步假定响应时间 X 的取值范围为 $0 \leq X \leq c .$ ，其中 c 未知，但将由 $\widehat {c} = \operatorname* {m a x} \{X_{i} : i = 1 , \ldots , n \} = 2 . 76$ 估计，其中 $\{X_{i} , i = 1 , \dots , n \}$ 为原始数据，$n = 5$ 为观测值个数。

把数据由小到大排列，令 $x_{( 1 )} \leq x_{( 2 )} \leq \cdot \cdot \cdot \leq x_{( n )}$ 表示这些排序后的取值。认为可能的最小取值为 0，故定义 $x_{( 0 )} = 0$ 。给每个区间 $x_{( i - 1 )} < x \leq x_{( i )}$ 赋予概率 $1 / n = 1 / 5$ ，如表 2 所示。由此得到的经验 cdf $\widehat F ( x )$ 见图 4。第 i 段线段的斜率为

逆 cdf 由下式计算：

其中 $( i - 1 ) / n < R \leq i / n$ 。例如，若生成随机数 $R_{1} = 0 . 71$ ，可以看出 $R_{1}$ 落在第四个区间内（在 $3 / 5 = 0 . 60$ 与 $4 / 5 = 0 . 80$ 之间）；于是由式 (10)，

生成过程的图解说明请参见图 4。

表 2 消防队响应时间数据汇总

图 4 消防队响应时间的经验 cdf。

在例 2 中，经验 cdf 把每一个数据点都表示了出来。若可以获得很大的样本（借助自动化数据采集，样本量可以达到数百乃至数万），那么先把数据汇总成区间数少得多的频数分布（frequency distribution），再对该频数分布拟合一条连续经验 cdf，可能更为方便（计算效率也更高）。为此只需对式 (10) 作一点点推广。现在第 i 段线段的斜率为：

其中 $c_{i}$ 是频数分布前 i 个区间的累积概率，$x_{( i - 1 )} < x \leq x_{( i )}$ 是第 i 个区间。逆 cdf 计算如下：

其中 $c_{i - 1} < R \leq c_{i}$ 。

例 3

假定已经收集到 100 个机器修理时间。数据按各个区间中的观测值个数汇总于表 3。例如，0 至 0.5 小时之间有 31 个观测值，0.5 至 1 小时之间有 10 个，依此类推。假定已知所有修理至少需要 15 分钟，因而总有 $X \geq 0 . 25$ 小时。于是置 $x_{( 0 )} = 0 . 25$ ，如表 3 和图 5 所示。

表 3 修理时间数据汇总

图 5 由修理时间数据的经验分布函数生成随机变量（$X \geq 0 . 25$）。

例如，假定所生成的第一个随机数为 $R_{1} = 0 . 83$ 。由于 $R_{1}$ 介于 $c_{3} = 0 . 66$ 与 $c_{4} = 1 . 00$ 之间，

再举一个例子，假定 $R_{2} = 0 . 33$ 。由于 $c_{1} = 0 . 31 < R_{2} \leq 0 . 41 = c_{2}$ ，

点 $( R_{2} = 0 . 33 , X_{2} = 0 . 6 )$ 同样示于图 5。

现在重新审视表 3 中的数据。数据的取值被限制在 $0 . 25 \leq X \leq 2 . 0$ 之内，但其背后的分布可能有更宽的取值范围。这给出了一条重要的理由，说明为何要为数据寻找一个理论概率分布（如伽马分布或 Weibull 分布）：这些分布确实允许更宽的取值范围——即 $0 \leq X < \infty$ 。另一方面，经验分布紧贴数据本身所呈现的内容，而数据往往是能够获得的信息中最好的来源。

当数据以频数区间的形式汇总时，建议采用相对较短的区间，因为这样做能更准确地刻画底层的 cdf。例如，对于表 3 中的修理时间数据，其观测值个数为 $n = 100$ ，若使用 10 至 20 个区间（这当然不算过多），本可以得到远为准确的估计，而不是像这里为便于说明而实际使用的四个相当宽的区间。

以下几点值得注意：

1. 这一过程的计算机化版本会随着区间个数 n 的增加而变得低效。系统化的计算机版本常称为查表生成方案（table-lookup generation scheme），因为给定一个 R 的取值，计算机程序必须在一组 $c_{i}$ 值中查找 R 所落入的区间 i，即满足

的区间 i。区间越多，若按这里所述的粗略方式实现，查找所需的时间平均而言就越长。编程实现该过程时，分析人员应当权衡估计 cdf 的准确性与计算效率之间的取舍。若可以获得大量的观测值，分析人员很可能决定把观测值归入 20 至 50 个区间，然后采用例 3 的过程——或者也可以采用更高效的查表过程，例如 Law [2007] 所描述的那种。

2. 在例 2 中，假定响应时间 X 满足 $0 \leq X \leq 2 . 76$ 。这一假定导致图 4 和表 2 中纳入了点 $x_{( 0 )} = 0$ 与 $x_{( 5 )} = 2 . 76$ 。如果事先已知 X 落在另一个范围内——例如，已知响应时间总是在 15 秒至 3 分钟之间，即

——那么就应当用点 $x_{( 0 )} = 0 . 25$ 与 $x_{( 6 )} = 3 . 0$ 来估计响应时间的经验 cdf。注意，由于纳入了新的点 $x_{( 6 )}$ ，现在共有六个区间而不是五个，每个区间赋予概率 $1 / 6 = 0 . 167$ 。习题 12 说明了对这些附加假定的运用。

1.6 无闭式逆函数的连续分布

有不少有用的连续分布，其 cdf 或其逆函数没有闭式表达式；例子包括正态分布、伽马分布和贝塔分布。因此，人们常说逆变换法不能用于这些分布的随机变量生成。但只要我们愿意近似逆 cdf，或者对 cdf 作数值积分并加以搜索，这一方法实际上是可以用的。尽管这种做法听上去不够精确，请注意即使是闭式逆函数，在计算机上求值时也需要作数值近似。例如，通过逆 cdf $X = F^{- 1} ( R ) = - \ln ( 1 - R ) / \lambda$ 生成服从指数分布的随机变量，就需要对对数函数作数值近似。因此，使用近似逆 cdf 与近似求值一个闭式逆函数之间，并没有本质差别。使用近似逆 cdf 的问题在于，其中有些在计算上求值缓慢。

为说明这一思想，考虑 Schmeiser [1979] 提出的标准正态分布逆 cdf 的一个简单近似：

当 $0 . 0013499 \leq R \leq 0 . 9986501$ 时，该近似的精度至少达到小数点后一位。表 4 把该近似与若干 R 取值下通过数值积分得到的精确值（保留四位小数）作了比较。还存在精度高得多、而复杂程度只是略高的近似。Bratley、Fox 和 Schrage [1996] 是若干分布此类近似的一个很好的来源。

表 4 标准正态分布的近似逆函数与精确值（保留四位小数）比较

1.7 离散分布

所有离散分布都可以用逆变换法生成，或者通过查表过程以数值方式生成，或者在某些情况下以代数方式生成（最终的生成方案为一个公式）。对某些分布还会使用其他技术，例如用于二项分布的卷积技术（convolution technique）。其中一些方法将在后面几节讨论。本小节给出若干例子，涵盖经验分布以及两个标准离散分布——（离散）均匀分布和几何分布。Bratley、Fox 和 Schrage [1996] 以及 Ripley [1987] 给出了这些分布和其他分布的高效查表过程。

例 4：经验离散分布

每天结束时，IHW 公司装卸码头上的发货批次要么是 0，要么是 1，要么是 2，其观察到的相对频数分别为 0.50、0.30 和 0.20。公司内部的顾问（consultants）受命开发一个模型，以改进装卸与搬运作业的效率；作为该模型的一部分，他们需要能够生成代表每天结束时码头发货批次数的取值 X。顾问们决定把 X 建模为一个离散随机变量，其分布如表 5 所示、见图 6。

表 5 发货批次 X 的分布

概率质量函数（pmf）$p(x)$ 为

cdf（$F ( x ) = P ( X \leq x )$ ）为

请记住，离散随机变量的 cdf 总是由若干水平线段组成，在该随机变量可能取到的各个点 x 上，线段发生大小为 $p(x)$ 的跳跃。例如，在图 6 中，在 $x = 0$ 处有一个大小为 $p ( 0 ) = 0 . 5$ 的跳跃，在 $x = 1$ 处有一个大小为 $p ( 1 ) = 0 . 3$ 的跳跃，在 $x = 2$ 处有一个大小为 $p ( 2 ) = 0 . 2$ 的跳跃。

图 6 发货批次 X 的 cdf。

对于离散随机变量的生成，逆变换法变成一个查表过程，但与连续变量的情形不同，它不需要插值。为说明这一过程，假定生成了 $R_{1} = 0 . 73$ 。在图上，如图 6 所示，先在纵轴上找到 $R_{1} = 0 . 73$，接着作一条水平线段直到撞上 cdf 中的一个“跳跃”，然后向横轴作垂线以得到所生成的随机变量。这里 $R_{1} = 0 . 73$ 被变换为 $X_{1} = 1$ 。这一过程与第 1.5 节中用于经验连续分布、并示于图 4 的过程类似，只是省去了最后一步——线性插值。

构造类似于表 6 的表格有助于实现这一查表过程。当生成 $R_{1} = 0 . 73$ 时，先找出 $R_{1}$ 所落入的区间。一般地，对 $R = R_{1}$ ，若

则置 $X_{1} = x_{i}$ 。这里，$r_{0} = 0 ; x_{0} = - \infty ;$ $x_{1} , x_{2} , \ldots , x_{n}$ 是该随机变量的可能取值；而 $r_{k} = p ( x_{1} ) + \cdot \cdot \cdot + p ( x_{k} ) , k = 1 , 2 , \cdot \cdot \cdot , n .$ 对于本例，$n = 3 , x_{1} = 0 , x_{2} = 1 ,$ 且 $x_{3} = 2 ;$ 因而 $r_{1} = 0 . 5 , r_{2} = 0 . 8 ,$ 且 $r_{3} = 1 . 0$ 。（注意，在所有情形下都有 $r_{n} = 1 . 0$ 。）

表 6 用于生成离散随机变量 X 的表

由于 $r_{1} = 0 . 5 < R_{1} = 0 . 73 \leq r_{2} = 0 . 8$ ，置 $X_{1} = x_{2} = 1$ 。生成方案汇总如下：

例 4 说明了查表过程；下一个例子说明可用于某些分布的代数方法。

例 5：离散均匀分布

考虑集合 $\{1 , 2 , \ldots , k \}$ 上的离散均匀分布，其 pmf 与 cdf 分别为

以及

令 $x_{i} = i$ 且 $r_{i} = p ( 1 ) + \cdot \cdot \cdot + p ( x_{i} ) = F ( x_{i} ) = i / k$ （$i = 1 , 2 , \dots , k$ ）。那么由不等式 (13) 可以看出，若所生成的随机数 R 满足

则通过令 $X = i$ 生成 X。现在可以对 i 解出不等式 (14)：

记 $\lceil y \rceil$ 为不小于 y 的最小整数。例如，$\left\lceil 7 . 82 \right\rceil = 8 , \left\lceil 5 . 13 \right\rceil = 6$ ，且 $\lceil - 1 . 32 \rceil = - 1$ 。当 $y \geq 0$ 时，$\lceil y \rceil$ 是一个向上取整的函数。这一记号与不等式 (15) 合在一起，给出了生成 X 的一个公式，即

例如，考虑生成一个在 $\{1 , 2 , \ldots , 10 \}$ 上均匀分布的随机变量 X。该随机变量 X 可以表示装上一辆卡车所要装载的托盘数。以表 A.1 作为随机数 R 的来源，并使用式 (16)（取 $k = 10$），得到

这里讨论的过程稍加修改即可生成取值范围为任意一组连续整数的离散均匀随机变量。习题 13 要求学生为这样一个情形设计生成过程。

例 6：几何分布

考虑 pmf 为

的几何分布，其中 $0 < p < 1$ 。其 cdf 为

其中 $x = 0 , 1 , 2 , . . .$ 。使用逆变换法［即不等式 (13)］，请回忆几何随机变量 X 在下述条件成立时取值为 x：

其中 R 是所生成的随机数，假定 $0 < R < 1$ 。对不等式 (17) 解出 x 的过程如下：

但 $1 - p < 1$ 意味着 $\ln \left( 1 - p \right) < 0$ ，因此

于是，对满足不等式 (18) 的那个整数 x，取 $X = x$ 。简言之，利用向上取整函数 $\lceil \cdot \rceil$ ，有

由于 $p$ 是固定参数，令 $\beta = - 1 / \ln ( 1 - p )$ 。则 $\beta > 0$ ，且由式 (19)，$X = \lceil - \beta \ln ( 1 - R ) - 1 \rceil$ 。由式 $( 1 )$ 可知，$- \beta \ln ( 1 - R )$ 是均值为 $\beta ;$ 的指数随机变量；因此，生成参数为 $p$ 的几何随机变量的一种办法是：生成一个（用任何方法）参数为 $\beta^{- 1} = - \ln ( 1 - p )$ 的指数随机变量，减去 1，再向上取整。

有时人们需要的几何随机变量 X 的可能取值为 $\{ q , q + 1 , q + 2 , \ldots \}$ ，其 pmf 为 $p ( x ) = p ( 1 - p ) ^ {x - q}$ （$x = q , q + 1 , \ldots$ ）。这类随机变量 X 可由式 (19) 生成为

最常见的情形之一是 $q = 1$ 。

例 7

生成取值范围为 $\{X \ge 1 \}$ 、均值为 2 的几何分布的三个取值。这样的几何分布其 pmf 为 $p ( x ) = p ( 1 - p ) ^ {x - 1}$ （$x = 1 , 2 , \dots$ ），均值为 $1 / p = 2 ,$ 即 $p = 1 / 2$ 。因此，可以由式 (20)（取 $q = 1 , p = 1 / 2$，且 $1 / \ln ( 1 - p ) = - 1 . 443$）生成 X。取自表 $\mathrm{A} . 1$ 的 $R_{1} = 0 . 932 , R_{2} = 0 . 105 ,$ 和 $R_{3} = 0 . 687$ 给出

习题 15 讨论几何分布的一个应用。

2 接受-拒绝法

假设分析人员需要设计一种方法，生成在 1/4 与 1 之间均匀分布的随机变量 X。可以按如下步骤进行：

步骤 1. 生成一个随机数 R。

步骤 2a. 若 $R \ge 1 / 4$ ，则接受 $X = R$ ，转到步骤 3。

步骤 2b. 若 $R < 1 / 4$ ，则拒绝 R，返回步骤 1。

步骤 3. 若还需要另一个 [1/4, 1] 上的均匀随机变量，则从步骤 1 开始重复上述过程。否则停止。

每次执行步骤 1，都必须生成一个新的随机数 R。在这一接受-拒绝技术中，步骤 2a 是"接受"，步骤 2b 是"拒绝"。概括来说，该技术是：生成具有某种分布（此处为 [0, 1] 上的均匀分布）的随机变量 R，直到某个条件（此处为 $R > 1 / 4$ ）被满足。当该条件最终被满足时，就可以计算出所求的随机变量 X（此处为 [1/4, 1] 上的均匀分布）（$X = R$ ）。这一过程之所以正确，是因为被接受的 R 的取值是条件化取值；也就是说，R 本身并不具有所期望的分布，但在事件 $\{R \ge 1 / 4\}$ 条件下的 R 才具有所期望的分布。要说明这一点，取 $1 / 4 \leq a < b \leq 1 ;$ 则

这正是 [1/4, 1] 上均匀分布的正确概率。式 (21) 说明，在 R 位于 1/4 与 1 之间这一条件下（R 的其他取值全部被舍去），R 的概率分布就是所期望的分布。因此，若 $1 / 4 \le R \le 1$ ，置 $X = R$ 。

接受-拒绝技术的效率在很大程度上取决于能否把拒绝的次数降到最低。在本例中，拒绝的概率为 $P ( R < 1 / 4 ) = 1 / 4 ,$ ，因此拒绝次数是一个服从几何分布的随机变量，其"成功"概率为 $p = 3 / 4$ ，平均拒绝次数为 $( 1 / p - 1 ) = 4 / 3 - 1 = 1 / 3$ 。（例 6 讨论了几何分布。）生成一个随机变量 X 所需的平均随机数 R 个数比拒绝次数多 1，故为 $4 / 3 = 1 . 33 .$ 。换句话说，生成 1000 个 X 的取值大约需要 1333 个随机数 R。

在当前这种情形下，生成 [1/4, 1] 上均匀变量的另一种过程是存在的——即式 (5)，它化为 $X = 1 / 4 + ( 3 / 4 ) R$ 。究竟是接受-拒绝法更高效，还是像逆变换法［式 (5)］这样的另一种过程更高效，取决于若干考虑。所使用的计算机、程序员的技术水平，以及接受-拒绝法为生成额外的（被拒绝的）随机数所带来的相对低效，都应当与另一种过程所需的计算量作比较。在实践中，对生成效率的关注通常交给专门的专家，他们通过广泛的试验来比较各种方法（即，直到某个仿真模型因所用发生器而开始需要过长的计算机运行时间为止）。

对于 [1/4, 1] 上的均匀分布，式 (5) 的逆变换法无疑比接受-拒绝法容易应用得多，效率也高得多。本例的主要目的是解释并阐明接受-拒绝法的基本概念。然而，对于某些重要的分布，如正态分布、伽马分布和贝塔分布，逆 cdf 没有闭式表达式，因而逆变换法难以使用。这些更为高级的技术由 Bratley、Fox 和 Schrage [1996] 以及 Law [2007] 作了总结。

在下面几小节中，将说明如何用接受-拒绝法生成泊松分布、非平稳泊松分布和伽马分布的随机变量。

2.1 泊松分布

一个均值为 $\alpha > 0$ 的泊松随机变量 N，其 pmf 为

然而更为重要的是，N 可以解释为一个单位时间内泊松到达过程（Poisson arrival process）的到达次数相继到达顾客的到达间隔时间 $A_{1} , A_{2} , \ldots$ 服从速率为 α 的指数分布（即 α 是单位时间内到达次数的均值）；此外，指数随机变量可由式 (3) 生成。因此，离散的泊松分布与连续的指数分布之间存在如下关系：

当且仅当

式 (22)（$N = n$ ）说的是在一个单位时间内恰好发生了 n 次到达；而不等式 (23) 说的是第 n 次到达发生在时刻 1 之前，而第 $( n + 1 )$ 次到达发生在时刻 1 之后。显然，这两个陈述是等价的。现在可以通过生成指数到达间隔时间直到某次到达（比如第 $n + 1$ 次）发生在时刻 1 之后为止，然后置 $N = n$ 来进行。

为了高效生成，不等式 (23) 通常先借助式 (3)（$A_{i} = ( - 1 / \alpha )$ ln R ）加以简化，得到

接下来两边同乘 $- \alpha$ ，这会使不等号方向反转，并利用"对数之和等于乘积之对数"这一事实，得到

最后，对任意数 x 利用关系 $e^{\ln x} = x$ ，得到

它与不等式 (23) 等价。生成泊松随机变量 N 的过程由以下步骤给出：

步骤 1. 置 $n = 0 , P = 1$ 。

步骤 2. 生成一个随机数 $R_{n + 1}$ ，并把 P 替换为 $P \cdot R_{n + 1}$ 。

步骤 3. 若 $\textit {P} < \textit {e} ^ {- \alpha}$ ，则接受 $N = n$ 。否则，拒绝当前的 n，把 n 加 1，返回步骤 2。

注意，步骤 2 完成后，P 等于不等式 (24) 中最右端的表达式。拒绝技术的基本思想在这里再次得到体现：若在步骤 3 中 $P \ge e^{- \alpha}$ ，则 n 被拒绝，生成过程必须至少再进行一次试探。

平均而言，生成一个泊松随机变量 N 需要多少个随机数呢？若 $N = n$ ，则需要 $n + 1$ 个随机数，因此平均个数由下式给出：

若泊松分布的均值 α 很大，这个数就相当大。

例 8

生成三个均值为 $\alpha = 0 . 2$ 的泊松随机变量。首先计算 $e^{- \alpha} = e^{- 0 . 2} = 0 . 8187$ 。接下来，从表 A.1 取一串随机数 R，按前面描述的步骤 1 至步骤 3 进行：

步骤 1. 置 $n = 0 , P = 1$ 。步骤 2. $R_{1} = 0 . 4357 , P = 1 \cdot R_{1} = 0 . 4357$ 。步骤 3. 由于 $P = 0 . 4357 < e^{- \alpha} = 0 . 8187$ ，接受 $N = 0$ 。

步骤 1–3.（$R_{1} = 0 . 4146$ 导致 $N = 0$ 。）步骤 1. 置 $n = 0 , P = 1$ 。步骤 2. $R_{1} = 0 . 8353 , P = 1 \cdot R_{1} = 0 . 8353 .$ 步骤 3. 由于 $P \ge e^{- \alpha}$ ，拒绝 $n = 0$ 并带着 $n = 1$ 返回步骤 2。

步骤 2. $R_{2} = 0 . 9952 , P = R_{1} R_{2} = 0 . 8313 .$ 步骤 3. 由于 $P \ge e^{- \alpha}$ ，拒绝 $n = 1$ 并带着 $n = 2$ 返回步骤 2。

步骤 2. $R_{3} = 0 . 8004 , P = R_{1} R_{2} R_{3} = 0 . 6654 .$ 步骤 3. 由于 $P < e^{- \alpha}$ ，接受 $N = 2$ 。

生成这三个泊松随机变量所需的计算汇总如下：

这里生成三个泊松随机变量（$N = 0 , N = 0 ,$ 和 $N = 2$）用了五个随机数 R；但从长期来看，若要生成比如 1000 个均值为 $\alpha = 0 . 2$ 的泊松随机变量，大约需要 $1000 ( \alpha + 1 )$ 即 1200 个随机数。

例 9

在 Peachtree 与 North Avenue 交叉口的公交站，公交车按泊松过程到达，平均每 15 分钟一辆。生成一个随机变量 N，表示某个 1 小时时段内到达的公交车数。现在，N 服从均值为每小时四辆车的泊松分布。首先计算 $e^{- \alpha} = e^{- 4} = 0 . 0183$ 。使用取自表 A.1 的 12 个随机数，得到如下汇总结果：

可以立即看出，α 的取值越大（此处 $\alpha = 4$ ），通常所需的随机数就越多；若需要 1000 个泊松随机变量，大约需要 $1000 ( \alpha + 1 ) = 5000$ 个随机数。

当 α 较大时（比如 $\alpha \geq 15$ ），这里所勾勒的拒绝技术就变得相当昂贵；所幸，基于正态分布的近似技术效果相当好。当均值 α 较大时，

近似服从均值为零、方差为 1 的正态分布；这一观察提示我们采用一种近似技术。首先按第 3.1 节的式 (28) 生成一个标准正态随机变量 Z，然后按下式生成所求的泊松随机变量 N：

其中 $\lceil \cdot \rceil$ 是第 1.7 节所述的向上取整函数。（若 $\alpha + {\sqrt {\alpha}} Z - 0 . 5 < 0$ ，则置 $N = 0$ 。）公式中所用的 $" 0 . 5 '$ 使向上取整函数变成了"四舍五入到最接近整数"的函数。式 (25) 不是一种接受-拒绝技术，但当它被用作接受-拒绝法的替代方案时，可以为生成均值很大的泊松随机变量提供一种相当高效而准确的方法。

2.2 非平稳泊松过程

另一类接受-拒绝方法（称为抽稀法，thinning）可用于生成非平稳泊松过程（nonstationary Poisson process, NSPP）的到达间隔时间，该过程的到达率为 $\lambda ( t ) , 0 \leq t \leq T$ 。NSPP 是到达率随时间而变化的到达过程。

例如，考虑表 7 给出的、每小时发生变化的到达率函数。抽稀法的思路是：以最快的速率（在本例中为每分钟 1/5 位顾客）生成一个平稳泊松到达过程，但只"接受"或保留其中一部分到达，抽稀的幅度恰好足以得到所期望的时变速率。下面给出通用算法，它生成 $\mathcal {T} _ {i}$ 作为第 i 次到达的时刻。请记住，在平稳泊松到达过程中，到达之间的时间服从指数分布。

步骤 1. 令

为到达率函数的最大值，并置 $t = 0$ 和 $i = 1$ 。

步骤 2. 从速率为 $\lambda^{*}$ 的指数分布中生成 E，并令 $t = t + E$ （这就是该平稳泊松过程的到达时刻）。

步骤 3. 从 $U ( 0 , 1 )$ 分布中生成随机数 R。若 $R \leq \lambda ( t ) / \lambda^{*}$ 则令 $\mathcal {T} _ {i} = t$ 且 $i = i + 1 .$ 步骤 4. 转到步骤 2。

若典型到达率与最大到达率之间的差异很大，抽稀算法可能效率不高。不过，抽稀法的优点在于它适用于任何可积的到达率函数，而不只是像本例这样分段为常数的函数。

表 7 NSPP 示例的到达率

例 10

对于表 7 中的到达率函数，生成前两个到达时刻。

步骤 1. $\lambda^{*} = \operatorname* {m a x} _ {0 \leq t \leq T} \lambda ( t ) = 1 / 5 , t = 0$ 且 $i = 1$ 。步骤 2. 对随机数 $R = 0 . 2130$ ，有 $E = - 5 \ln ( 0 . 213 ) = 13 . 13$ ，且 $t = 0 + 13 . 13 = 13 . 13$ 。

步骤 3. 生成 $R = 0 . 8830$ 。由于 $R = 0 . 8830 \not\leq \lambda ( 13 . 13 ) / \lambda^{*} = ( 1 / 15 ) / ( 1 / 5 ) = 1 / 3$ ，故不生成该次到达。

步骤 4. 转到步骤 2。

步骤 2. 对随机数 $R = 0 . 5530$ ，有 $E = - 5 \ln ( 0 . 553 ) = 2 . 96$ ，且 $t = 13 . 13 + 2 . 96 = 16 . 09$ 。步骤 3. 生成 $R = 0 . 0240$ 。由于 $R = 0 . 0240 \leq \lambda ( 16 . 09 ) / \lambda^{*} = ( 1 / 15 ) / ( 1 / 5 ) = 1 / 3$ ，故置 $T_{1} = t = 16 . 09$ 且 $i = i + 1 = 2$ 。

步骤 4. 转到步骤 2。

步骤 2. 对随机数 $R = 0 . 0001$ ，有 $E = - 5 \ln ( 0 . 0001 ) = 46 . 05$ ，且 $t = 16 . 09 + 46 . 05 = 62 . 14$ 。

步骤 3. 生成 $R = 0 . 1443$ 。由于 $R = 0 . 1443 \leq \lambda ( 62 . 14 ) / \lambda^{*} = ( 1 / 12 ) / ( 1 / 5 ) = 5 / 12$ ，故置 $T_{2} = t = 62 . 14$ 且 $i = i + 1 = 3$ 。

步骤 4. 转到步骤 2。

2.3 伽马分布

人们已经开发出若干用于生成伽马随机变量的接受-拒绝技术。（参见 Bratley、Fox 和 Schrage [1996] 以及 Law [2007]。）其中较为高效的一种由 Cheng [1977] 提出；对形状参数 β 的任意取值 $\beta \geq 1$ ，其平均试探次数介于 1.13 与 1.47 之间。若形状参数 β 为整数，比如 $\beta = k ,$ 则可以采用例 12 中的卷积技术，因为厄兰分布（Erlang distribution）是更一般的伽马分布的一个特例。另一方面，这里所描述的接受-拒绝技术对厄兰分布而言将是极为高效的方法，尤其当 $\beta = k$ 较大时。该例程生成具有尺度参数 θ 和形状参数 β 的伽马随机变量——也就是说，均值为 $1 / \theta$ ，方差为 $1 / \beta \theta^{2}$ 。其步骤如下：

步骤 1. 计算 $a = 1 / ( 2 \beta - 1 ) ^ {1 / 2} , b = \beta - \ln 4$ 。

步骤 2. 生成 $R_{1}$ 和 $R_{2}$ 。置 $V = R_{1} / ( 1 - R_{1} )$ 。

步骤 3. 计算 $X = \beta V^{a}$ 。

步骤 4a. 若

则拒绝 X 并返回步骤 2。

步骤 4b. 若

则把 X 作为所求的随机变量。

步骤 4b 所生成的随机变量，其均值与方差都等于 $\beta .$ 。若希望均值为 $1 / \theta$ 且方差为 $1 / \beta \theta^{2}$ ，则还要加入步骤 5。

（步骤 5. 把 X 替换为 $X / ( \beta \theta ) .$ ）

所有接受-拒绝方法的基本思想在这里再次得到说明，但本例的证明超出了本书的范围。在步骤 3 中，$X = \beta V^{a} = \beta [ R_{1} / ( 1 - R_{1} ) ] ^ {a}$ 并不服从伽马分布，但在步骤 4a 中拒绝 X 的某些取值，保证了步骤 4b 中被接受的取值确实服从伽马分布。

例 11

已发现某高产量糖果制造机的停机时间服从伽马分布，均值为 2.2 分钟，方差为 2.10 分钟^2。因此，$1 / \theta = 2 . 2$ 且 $1 / \beta \theta^{2} = 2 . 10$ ，两者合起来意味着 $\beta = 2 . 30$ 且 $\theta = 0 . 4545$ 。

步骤 1. $a = 0 . 53 , b = 0 . 91$ 。

步骤 2. 生成 $R_{1} = 0 . 832 , R_{2} = 0 . 021$ 。置 $V = 0 . 832 / ( 1 - 0 . 832 ) = 4 . 952$ 。

步骤 3. 计算 $X = 2 . 3 ( 4 . 952 ) ^ {0 . 53} = 5 . 37$ 。

步骤 4. $X = 5 . 37 > 0 . 91 + [ 2 . 3 ( 0 . 53 ) + 1 ] \ln ( 4 . 952 ) - \ln [ ( 0 . 832 ) ^ {2} ( 0 . 021 ) ] = 8 . 68$ ，故拒绝 X 并返回步骤 2。

步骤 2. 生成 $R_{1} = 0 . 434 , R_{2} = 0 . 716$ 。置 $V = 0 . 434 / ( 1 - 0 . 434 ) = 0 . 767$ 。步骤 3. 计算 $X = 2 . 3 ( 0 . 767 ) ^ {0 . 53} = 2 . 00$ 。

步骤 4. 由于 $X = 2 . 00 \leq 0 . 91 + [ 2 . 3 ( 0 . 53 ) + 1 ] \ln ( 0 . 767 ) - \ln [ ( 0 . 434 ) ^ {2} ( 0 . 716 ) ] = 2 . 32$ ，接受 X。

步骤 5. 用 $\beta \theta = 1 . 045$ 除 X ，得 $X = 1 . 91$ 。

本例用了两次试探（即一次拒绝）才生成一个可接受的服从伽马分布的随机变量；但平均而言，若要生成比如 1000 个伽马随机变量，该方法将需要 1130 到 1470 次试探，或者等价地需要 2260 到 2940 个随机数。该方法用于手工计算稍嫌繁琐，但在计算机上编程很容易，并且是已知的最有效的伽马发生器之一。

3 特殊性质

“特殊性质"正如其名称所暗示的那样：它们是一些基于概率分布特定族之特征的随机变量生成技术，而不是像逆变换法或接受-拒绝法那样的通用技术。

3.1 正态分布与对数正态分布的直接变换

人们已经开发出多种生成正态分布随机变量的方法。然而，逆变换法难以应用，因为逆 cdf 无法写成闭式形式。标准正态 cdf 为

本节描述一种直观而吸引人的直接变换，它产生一对相互独立的、均值为零方差为 1 的标准正态随机变量。该方法由 Box 与 Muller [1958] 提出。尽管它不如许多更现代的技术高效，但用科学计算语言（如 FORTRAN、C、C++、Visual Basic 或 Java）编写程序十分容易。随后我们说明如何把一个标准正态随机变量变换为均值为 $\mu$ 、方差为 $\sigma^{2}$ 的正态随机变量。一旦我们有了（用本方法或任何其他方法）从 $N ( \mu , \sigma^{2} )$ 分布生成 $X$ 的方法，就可以利用直接变换 $Y = e^{X}$ 生成参数为 $\mu$ 和 $\sigma^{2}$ 的对数正态（lognormal）随机变量 Y 。（请记住，$\mu$ 和 $\sigma^{2}$ 并不是对数正态分布的均值和方差。）