随机化和采样方法
4.90/5 (82投票s)
提供了许多应用程序可以使用底层（伪）随机数生成器进行采样的方法，并包含其中许多方法的伪代码。
引言
本页目录随机化方法和采样方法。随机化或采样方法由“随机数源”驱动，并生成称为随机变量的数字或其他值。这些变量是随机化的结果。（“随机数源”在实践中通常通过所谓的伪随机数生成器或PRNG来模拟。）本文档涵盖了许多方法，包括—
本页重点介绍精确从描述的分布中采样的随机化和采样方法，而不会引入除输入中已存在的错误之外的任何其他错误（并假设存在理想的“随机数源”）。如果可供选择的值数量有限，情况将是如此。但对于正态分布和其他具有无限多值的分布，总会涉及某种程度的近似；在这种情况下，本页的重点是最大程度地减少其引入的错误的方法。
本文档展示了许多方法的伪代码，并且可以下载示例 Python 代码，该代码实现了本文档中的许多方法，以及代码文档。
本页介绍的随机化方法假定我们有一个无限的、独立随机选择的、均匀分布的数字源。有关更多信息，请参阅附录中的“随机数源”。
总的来说，以下内容不在本文档的讨论范围之内
- 本文档不涵盖如何为特定应用程序选择底层PRNG（或模拟“随机数源”的设备或程序），包括在安全性、性能和质量方面。我在另一份文档中写了更多关于建议的内容。
- 本文档不包含针对特定PRNG的算法，例如Mersenne Twister、PCG、xorshift、线性同余生成器或基于哈希函数的生成器。
- 本文档不涵盖如何测试PRNG的正确性或充分性，同样也适用于模拟“随机数源”的其他设备和程序。测试涵盖在例如“测试PRNG以获得高质量随机性”中。
- 本文档不解释如何为PRNG指定或生成“种子”。这在其他地方有详细介绍。
- 本文档不展示如何生成随机安全参数，例如加密密钥。
- 本文档不涵盖随机性提取（也称为去偏、去抖动或白化）。请参阅我的关于随机性提取的说明。
- “方差缩减”方法，例如重要性采样或公共随机数，不在本文档的讨论范围之内。
此外，本页不侧重于用于计算机图形渲染的采样方法（例如泊松盘采样、多重重要性采样、蓝噪声和梯度噪声），因为此应用程序倾向于将性能和视觉可接受性置于准确性和精确采样之上。
关于本文档
本文档是开源的；更新版本请参阅 源代码 或其在 GitHub 上的 渲染版本。您可以就本文档通过 CodeProject 或 GitHub issues 页面 发送评论。
本文的受众是具备数学知识但对微积分不熟悉或非常不熟悉的计算机程序员。
我鼓励读者实现本文中的任何算法，并报告他们的实现体验。特别是，我寻求有关以下方面的评论
- 本文中的算法是否易于实现？是否有人在阅读本文后可以为该算法编写代码？
- 本文是否有应纠正的错误？
- 是否有方法可以使本文对目标受众更有用？
欢迎对本文档的其他方面提出意见。
目录
符号
本文档中
- 本文档使用伪代码约定。
- 区间符号表示如下：[
a,b) 表示“a或更大，但小于b”。(a,b) 表示“大于a，但小于b”。(a,b] 表示“大于a且小于或等于b”。[a,b] 表示“a或更大且b或更小”。 log1p(x)等同于ln(1 + x)，并且是ln(1 + x)的“稳健”替代方案，其中x是浮点数（Pedersen 2018）[^1]。MakeRatio(n, d)创建一个分子为n、分母为d的有理数。Sum(list)计算给定整数或有理数列表中数字的总和。（朴素地对浮点数求和可能会引入舍入误差。）
均匀随机整数
本节展示了如何借助“随机数源”（或模拟该源的设备或程序）来生成独立的均匀随机整数。
本节描述了四种方法：RNDINT、RNDINTEXC、RNDINTRANGE、RNDINTEXCRANGE。其中，接下来描述的 RNDINT 可以作为其余方法的依据。
RNDINT：[0, N] 中的随机整数
在本文档中，RNDINT(maxInclusive) 是本文档的核心方法；它从“随机数源”生成[0, maxInclusive] 区间内的独立均匀整数[^2]。在下面的伪代码中，它实现了 RNDINT
NEXTRAND()读取“（均匀）随机数源”生成的下一个数字，该源在附录中定义（一个无限的数字源，其中每个数字都是独立选择的，且分布均匀）。如附录中所述，伪随机数生成器在实践中可以模拟此源。对于此方法，源可以是均匀分布，也可以是非均匀分布。MODULUS是该源可能产生的不同结果的数量。
具体来说
| 如果底层源产生 | 那么 NEXTRAND() 是 | 而 MODULUS 是 |
|---|---|---|
| 非均匀数字[^3]。 | 将底层源的输出作为随机性提取技术的输入，以生成独立的无偏随机比特（零或一）的新源的下一个比特。 | 2. |
| 下面未描述的均匀数字。 | 同上。 | 2n。 |
| 均匀的 32 位非负整数。 | 源的下一个数字。 | 232. |
| 均匀的 64 位非负整数。 | 源的下一个数字。 | 264. |
| 区间 [0, n) 内的均匀整数。 | 源的下一个数字。 | n. |
| 区间 [0, 1) 内的均匀数字，已知它们以 p 的倍数间隔分布（例如 dSFMT）。 | 源的下一个数字，乘以 p。 | 1/p。 |
区间 [0, 1) 内的均匀数字，其中 [0, 0.5) 中的数字和 [0.5, 1) 中的数字发生的概率相等（例如 Java 的 Math.random()）。 | 如果源输出的数字小于 0.5，则为 0，否则为 1。 | 2. |
METHOD RndIntHelperNonPowerOfTwo(maxInclusive)
if maxInclusive <= MODULUS - 1:
// NOTE: If the programming language implements
// division with two integers by discarding the result's
// fractional part, the division can be used as is without
// using a "floor" function.
nPlusOne = maxInclusive + 1
maxexc = floor((MODULUS - 1) / nPlusOne) * nPlusOne
while true // until a value is returned
ret = NEXTRAND()
if ret < nPlusOne: return ret
if ret < maxexc: return rem(ret, nPlusOne)
end
else
cx = floor(maxInclusive / MODULUS) + 1
while true // until a value is returned
ret = cx * NEXTRAND()
// NOTE: The addition operation below should
// check for integer overflow and should reject the
// number if overflow would result.
ret = ret + RNDINT(cx - 1)
if ret <= maxInclusive: return ret
end
end
END METHOD
METHOD RndIntHelperPowerOfTwo(maxInclusive)
// NOTE: Finds the number of bits minus 1 needed
// to represent MODULUS (in other words, the number
// of random bits returned by NEXTRAND() ). This will
// be a constant here, though.
modBits = ln(MODULUS)/ln(2)
// Lumbroso's Fast Dice Roller.
x = 1
y = 0
nextBit = modBits
rngv = 0
while true // until a value is returned
if nextBit >= modBits
nextBit = 0
rngv = NEXTRAND()
end
x = x * 2
y = y * 2 + rem(rngv, 2)
rngv = floor(rngv / 2)
nextBit = nextBit + 1
if x > maxInclusive
if y <= maxInclusive: return y
x = x - maxInclusive - 1
y = y - maxInclusive - 1
end
end
END METHOD
METHOD RNDINT(maxInclusive)
// maxInclusive must be 0 or greater
if maxInclusive < 0: return error
if maxInclusive == 0: return 0
if maxInclusive == MODULUS - 1: return NEXTRAND()
// NOTE: Finds the number of bits minus 1 needed
// to represent MODULUS (if it's a power of 2).
// This will be a constant here, though.
modBits=ln(MODULUS)/ln(2)
if floor(modBits) == modBits // Is an integer
return RndIntHelperPowerOfTwo(maxInclusive)
else
return RndIntHelperNonPowerOfTwo(maxInclusive)
end
END METHOD
下表显示了为均匀随机选择整数而提出的算法。有些是无偏的（精确的），有些是有偏的，但总的来说，除非算法在最坏情况下永远运行，否则该算法只能是无偏的。所列算法以 n 作为参数，其中 n = maxInclusive + 1，因此从区间 [0, n) 中采样。（“无偏？”列表示算法是否生成无偏的随机整数，即使 n 不是 2 的幂。）
| 算法 | 最优？ | 无偏？ | 时间复杂度 |
|---|---|---|---|
| 拒绝采样：在更大的区间内采样，直到采样的数字适合较小的区间。 | 不总是 | 是 | 最坏情况下永远运行 |
乘法-移位约简：生成 bignumber，一个由 k 个无偏比特组成的整数，其中 k 远大于 n，然后计算 (bignumber * n) >> k（参见（Lemire 2016）[^4]、（Lemire 2018）[^5]以及 M. O'Neill 调查的“整数乘法”算法）。 | 否 | 否 | 恒定 |
模约简：生成 bignumber 如上所述，然后计算 rem(bignumber, n)。 | 否 | 否 | 恒定 |
| 快速骰子滚筒 (Lumbroso 2013)[^6]（参见上面的伪代码）。 | 是 | 是 | 最坏情况下永远运行 |
| Math Forum (2004)[^7] 或 (Mennucci 2018)[^8]（批处理/回收随机比特）。 | 是 | 是 | 最坏情况下永远运行 |
| M. O'Neill 调查的“FP Multiply”。 | 否 | 否 | 恒定 |
| O'Neill 在“结论”部分介绍的算法。 | 否 | 是 | 最坏情况下永远运行 |
| M. O'Neill 调查的“去偏”和“带拒绝的位掩码”。 | 否 | 是 | 最坏情况下永远运行 |
注释
RNDINT作为二叉树遍历器。 Donald E. Knuth 和 Andrew C. Yao (1976)[^9] 表明，任何使用随机无偏比特（每个比特是 0 或 1 的概率相等）生成随机整数的算法都可以描述为二叉树（也称为DDG 树或离散分布生成树）。（这也适用于RNDINT算法。）随机无偏比特沿着此树的路径，每个叶子（终端节点）代表一个结果。对于RNDINT(maxInclusive)，有n = maxInclusive + 1个结果，每个结果的概率为1/n。
Knuth 和 Yao 表明，任何最优 DDG 树算法平均需要至少log2(n)比特，最多需要log2(n) + 2比特（其中log2(x) = ln(x)/ln(2)）。[^10] 但他们也表明，为了使算法无偏（精确），它必须在最坏情况下永远运行，即使它平均使用很少的随机比特（也就是说，一般情况下无法在保持无偏性的同时“修复”此最坏情况）。这是因为1/n将具有无限的二进制展开，除非n是 2 的幂，因此生成的 DDG 树将需要无限深，或在树的末尾包含“拒绝叶”。
例如，模约简方法可以由一个 DDG 树表示，其中拒绝叶被标记的结果替换，但该方法是有偏的，因为只有某些结果可以通过这种方式替换拒绝叶。出于同样的原因，在固定次数后停止拒绝采样器也会导致偏差。但是，哪些结果会以这种方式产生偏差取决于算法。- 减少“比特浪费”。 如上所述，任何
RNDINT实现平均每次选择的整数至少需要log2(n)比特，但大多数实现使用的比特远多于此。有多种方法可以将算法的比特数接近log2(n)。它们包括批处理、比特回收和随机性提取，并且例如在 Math Forum 页面以及上面引用的 Lumbroso 和 Mennucci 的论文以及 Devroye 和 Gravel (2020, section 2.3)[^11] 中进行了讨论。*批处理示例*：要生成 0 到 9 之间的三个数字，我们可以调用RNDINT(999)来生成 [0, 999] 范围内的整数，然后将返回的数字分解为三个数字。- 模拟骰子。 如果我们有一个（虚拟的）公平的 *p* 面骰子，我们如何用它来模拟 *k* 面骰子的掷骰？除非“能整除 k 的每个素数也能整除 p”（参见 B. Kloeckner 的“用骰子模拟骰子”，2008），否则无法在不“浪费”随机性的情况下完成。（参见我的关于随机性提取的说明）。但是，随机性提取可以将骰子掷骰变成无偏比特，因此前面的讨论适用于此。
RNDINTRANGE：[N, M] 中的随机整数
生成*区间* [`minInclusive`, `maxInclusive`] 内的随机整数的朴素方法如下所示，对于非负整数和任意精度整数效果很好。
METHOD RNDINTRANGE(minInclusive, maxInclusive)
// minInclusive must not be greater than maxInclusive
if minInclusive > maxInclusive: return error
return minInclusive + RNDINT(maxInclusive - minInclusive)
END METHOD
然而，朴素方法的效果不佳，如果整数格式可以表示负数和非负整数，并且 `maxInclusive` 和 `minInclusive` 之间的差值超过了该格式的最高可能整数。对于可以表示—
- 范围 [-1 -
MAXINT,MAXINT] 内的每个整数（例如 Javaint、short或long），或 - 范围 [-
MAXINT,MAXINT] 内的每个整数（例如 Javafloat和double以及 .NET 的System.Decimal实现），
其中 MAXINT 是一个大于 0 的整数，可以使用以下 RNDINTRANGE 伪代码。
METHOD RNDINTRANGE(minInclusive, maxInclusive)
// minInclusive must not be greater than maxInclusive
if minInclusive > maxInclusive: return error
if minInclusive == maxInclusive: return minInclusive
if minInclusive==0: return RNDINT(maxInclusive)
// Difference does not exceed maxInclusive
if minInclusive > 0 or minInclusive + MAXINT >= maxInclusive
return minInclusive + RNDINT(maxInclusive - minInclusive)
end
while true // until a value is returned
ret = RNDINT(MAXINT)
// NOTE: For case 1, use the following line:
if RNDINT(1) == 0: ret = -1 - ret
// NOTE: For case 2, use the following three lines
// instead of the preceding line; these lines
// avoid negative zero
// negative = RNDINT(1) == 0
// if negative: ret = 0 - ret
// if negative and ret == 0: continue
if ret >= minInclusive and ret <= maxInclusive: return ret
end
END METHOD
RNDINTEXC：[0, N) 中的随机整数
RNDINTEXC(maxExclusive)，它生成*区间* [0, maxExclusive) 内的随机整数，可以如下实现[^12]
METHOD RNDINTEXC(maxExclusive)
if maxExclusive <= 0: return error
return RNDINT(maxExclusive - 1)
END METHOD
注意：
RNDINTEXC未被指定为核心随机生成方法，因为它在此方法中更难用随机比特填充流行的整数格式中的整数。
RNDINTEXCRANGE：[N, M) 中的随机整数
RNDINTEXCRANGE 返回*区间* [minInclusive, maxExclusive) 内的随机整数。它可以使用RNDINTRANGE 实现，如下面的伪代码所示。
METHOD RNDINTEXCRANGE(minInclusive, maxExclusive)
if minInclusive >= maxExclusive: return error
if minInclusive >=0: return RNDINTRANGE(
minInclusive, maxExclusive - 1)
while true // until a value is returned
ret = RNDINTRANGE(minInclusive, maxExclusive)
if ret < maxExclusive: return ret
end
END METHOD
均匀随机比特
惯用语 RNDINT((1 << b) - 1) 是生成*均匀随机* b 位整数（最大值为 2b - 1）的朴素方法。
在实践中，内存通常被分成*字节*，或区间 [0, 255] 内的 8 位整数。在这种情况下，可以填充一个内存块的随机比特—
- 将块中的每个字节设置为
RNDINT(255)，或 - 通过 PRNG（或模拟“随机数源”的其他设备或程序），如果它一次输出一个或多个 8 位块。
使用 RNDINT 系列的示例
- 要以相等的概率选择 −1 或 1（*Rademacher 分布*），可以使用以下惯用语之一：
(RNDINT(1) * 2 - 1)或(RNDINTEXC(2) * 2 - 1)。 - 要生成一个可被正整数 (
DIV) 整除的随机整数，请使用任何方法（例如RNDINT）生成整数，令X为该整数，然后计算X - rem(X, DIV)（如果X >= 0），或X - (DIV - rem(abs(X), DIV))（否则）。（根据方法，结果整数可能超出范围，在这种情况下，该过程应重复。） - N×M 网格上的随机二维点可以表示为单个整数，如下所示
- 要生成随机的 N×M 点
P，请生成P = RNDINT(N * M - 1)（因此P在区间 [0,N * M) 内）。 - 要将点
P转换为其二维坐标，请生成[rem(P, N), floor(P / N)]。（每个坐标从 0 开始。） - 要将二维坐标
coord转换为 N×M 点，请生成P = coord[1] * N + coord[0]。
- 要生成随机的 N×M 点
- 模拟掷 N 面骰子（N > 1）：
RNDINTRANGE(1, N)，它以相等的概率选择区间 [1, N] 中的数字。 - 生成一个具有一个十进制数字的随机整数：
RNDINTRANGE(0, 9)。 - 生成一个具有 N 个十进制数字（N ≥ 2）的随机整数，其中第一个数字不能为 0：
RNDINTRANGE(pow(10, N-1), pow(10, N) - 1)。 - 以
step的增量，以相等的概率选择区间 [mn,mx) 中的数字：mn+step*RNDINTEXC(ceil((mx-mn)/(1.0*step)))。 - 在区间 [0,
X) 中随机选择一个整数- 并倾向于中间的数字：
floor((RNDINTEXC(X) + RNDINTEXC(X)) / 2)。 - 并倾向于较大的数字：
max(RNDINTEXC(X), RNDINTEXC(X))。 - 并倾向于较小的数字：
min(RNDINTEXC(X), RNDINTEXC(X))。 - 并强烈倾向于较大的数字：
max(RNDINTEXC(X), RNDINTEXC(X), RNDINTEXC(X))。 - 并强烈倾向于较小的数字：
min(RNDINTEXC(X), RNDINTEXC(X), RNDINTEXC(X))。
- 并倾向于中间的数字：
随机化技术
本节介绍常用的随机化技术，如洗牌、选择多个唯一项以及创建随机文本字符串。
布尔（真/假）条件
要生成具有指定概率的真条件，请在 if 条件中使用以下惯用语
- 以相等的概率真或假：
RNDINT(1) == 0。 - 以 X% 的概率为真：
RNDINTEXC(100) < X。 - 以 X/Y 的概率为真（*伯努利试验*）：
RNDINTEXC(Y) < X。 - 以 X 比 Y 的几率（odds）为真：
RNDINTEXC(X + Y) < X。
以下辅助方法生成 1 的概率为 x/y，否则为 0
METHOD ZeroOrOne(x,y)
if RNDINTEXC(y)<x: return 1
return 0
END METHOD
该方法也可以按以下方式实现（如 Lumbroso (2013, Appendix B)[^6] 所指出）
// NOTE: Modified from Lumbroso
// Appendix B to add 'z==0' and error checks
METHOD ZeroOrOne(x,y)
if y <= 0: return error
if x==y: return 1
z = x
while true // until a value is returned
z = z * 2
if z >= y
if RNDINT(1) == 0: return 1
z = z - y
else if z == 0 or RNDINT(1) == 0: return 0
end
END METHOD
注意： 概率可以是无理数或有理数。有理数概率的形式为
n/d，可以使用上面的ZeroOrOne进行模拟。无理数概率（如exp(-x/y)或ln(2)）具有无限的数字展开（0.ddddd....），它们需要特殊的算法来模拟；请参阅我关于伯努利工厂算法的页面上的“无理数常量算法”。示例:
- 以 3/8 的概率为真：
RNDINTEXC(8) < 3。- 以 100 比 1 的几率（odds）为真：
RNDINTEXC(101) < 1。- 以 20% 的概率为真：
RNDINTEXC(100) < 20。- 使用更少的随机比特而不是朴素方法，在 [0,
y) 中生成随机整数，或者如果该数字小于x则返回 -1：if ZeroOrOne(x, y) == 1: return -1; else: return RNDINTEXCRANGE(x, y)。
随机采样
本节包含从集合中选择一项或多项的方法，其中集合中的每个项都有与其他项相同的被选中的机会。这称为*随机采样*，可以*带放回*或*无放回*地进行。
有放回采样：从列表中随机选择一项
有放回采样本质上意味着随机抽取一项并将其放回。要从列表中选择一项—
- 其大小可预知，使用惯用语
list[RNDINTEXC(size(list))]；或 - 其大小不可预知，生成
RandomKItemsFromFile(file, 1)，在稍后给出的伪代码中（结果将是一个 1 项列表，或如果没有任何项则为空列表）。
无放回采样：选择多个唯一项
无放回采样本质上意味着不放回地抽取一项。有几种方法可以从 n 个可用项或值中均匀随机地选择 k 个唯一项或值，具体取决于 n 是否已知以及 n 和 k 的大小。
- 如果
n未预知：使用水库抽样方法；参见稍后伪代码中的RandomKItemsFromFile方法。 - 如果
n相对较小（例如，有 200 个可用项，或有 0 到 200 的数字范围可供选择）- 如果必须从列表中按*相对（索引）顺序*选择项，或者如果
n为 1，则使用RandomKItemsInOrder（稍后给出）。 - 否则，如果采样项的顺序不重要，并且每个项都可以通过其*索引*（项的位置，从 0 开始的整数）来派生，而无需从列表中查找：使用
RandomKItemsFromFile方法。[^13] - 否则，如果
k远小于n，则按项目 3 的方式进行。 - 否则，以下任何一种方法都可以随机顺序选择
k项- 将所有项存储在列表中，洗牌该列表，然后从该列表中选择前
k项。 - 如果项已存储在列表中且列表顺序可以更改，则洗牌该列表并从洗牌后的列表中选择前
k项。 - 如果项已存储在列表中且列表顺序不能更改，则将这些项的索引存储在另一个列表中，洗牌后者列表，然后从后者列表中选择前
k个索引（或与这些索引对应的项）。
- 将所有项存储在列表中，洗牌该列表，然后从该列表中选择前
- 如果必须从列表中按*相对（索引）顺序*选择项，或者如果
- 如果
k远小于n且项的顺序是随机的或不重要- 如果项存储在顺序可以更改的列表中： 对该列表进行*部分洗牌*，然后选择该列表的*最后*
k项。*部分洗牌*按“洗牌”部分所述进行，不同之处在于部分洗牌在进行k次交换后停止（将一个项与其自身交换算作一次交换）。 - 否则：创建另一个空列表
newlist，并创建一个键/值数据结构，例如哈希表。然后，对于区间 [0, k−1] 中的每个整数i，执行j = RNDINTEXC(n-i); AddItem(newlist, HGET(j,j)); HSET(j,HGET(n-i-1,n-i-1))，其中HSET(k,v)将哈希表中键为k的项设置为v，HGET(k,v)获取表中键为k的项，如果没有该项则返回v（Ting 2021）[^14]。新列表存储选定项的索引，按随机顺序。
- 如果项存储在顺序可以更改的列表中： 对该列表进行*部分洗牌*，然后选择该列表的*最后*
- 如果
n - k远小于n，项存储在列表中，并且采样项的顺序不重要- 如果列表顺序可以更改： 对该列表进行*部分洗牌*，不同之处在于进行了
n-k次交换而不是k次，然后选择该列表的*前*k项。（注意 5 中的“洗牌”不能使用。） - 否则，如果
n不是很大（例如，小于 5000）： 将这些项的索引存储在另一个列表中，对后者列表进行*部分洗牌*，不同之处在于进行了n-k次交换而不是k次，然后选择该列表的*前*k个索引（或与这些索引对应的项）。（注意 5 中的“洗牌”不能使用。） - 否则：按项目 5 的方式进行。
- 如果列表顺序可以更改： 对该列表进行*部分洗牌*，不同之处在于进行了
-
否则（例如，如果选择 32 位或更大整数，则
n为 232，或者n否则非常大）- 如果必须按*相对（索引）顺序*选择项：令
n2 = floor(n/2)。生成h = PolyaEggenberger(k, n2, n, -1)。通过递归运行此算法（项目 1 到 5），从列表[0, 1, ..., n2 - 1]中按相对顺序采样h个整数，然后通过递归运行此算法，从列表[n2, n2 + 1, ..., n - 1]中按相对顺序采样k - h个整数。这样选择的整数是所需项在相对（索引）顺序中的索引（Sanders et al. 2019）[^15]。 -
否则，创建一个数据结构来存储已选索引。当随机选择一个新项的索引时，如果它不在数据结构中，则添加它，否则选择一个新的随机索引。重复此过程，直到向该数据结构添加了
k个索引。合适的结构示例包括—许多应用程序出于各种原因需要生成“唯一随机”值来标识数据库记录或其他共享资源。有关生成此类值的方法，请参阅我的推荐文档。
- 如果必须按*相对（索引）顺序*选择项：令
洗牌
Fisher–Yates 洗牌方法可以洗牌列表（将其项置于随机顺序），使得该列表的所有排列（排列）都以相同的概率发生。但是，该方法也很容易出错 — 请参阅（Atwood 2007）[^16]。下面的伪代码旨在洗牌列表的内容。
METHOD Shuffle(list)
// NOTE: Check size of the list early to prevent
// `i` from being less than 0 if the list's size is 0 and
// `i` is implemented using a nonnegative integer
// type available in certain programming languages.
if size(list) >= 2
// Set i to the last item's index
i = size(list) - 1
while i > 0
// Choose an item ranging from the first item
// up to the item given in `i`. Observe that the item
// at i+1 is excluded.
j = RNDINTEXC(i + 1)
// Swap item at index i with item at index j;
// in this case, i and j may be the same
tmp = list[i]
list[i] = list[k]
list[k] = tmp
// Move i so it points to the previous item
i = i - 1
end
end
// NOTE: An implementation can return the
// shuffled list, as is done here, but this is not required.
return list
END METHOD
注释
j = RNDINTEXC(i + 1)不能用j = RNDINTEXC(size(list))替换，因为它会引入偏差。如果该行被替换为j = RNDINTEXC(i)，则结果是 Sattolo 算法（它从具有循环的排列中选择），而不是 Fisher–Yates 洗牌。- 在洗牌方面，伪随机数生成器（或任何模拟“随机数源”的设备）的选择很重要；请参阅我的关于洗牌的推荐文档。
- （Bacher et al., 2015）[^17] 描述了一种可以并行进行的洗牌算法。
- 错排是每个元素都移动到不同位置的排列。可以使用以下方法生成随机错排（Merlini et al. 2007）[^18]：(1) 修改
Shuffle，在j = RNDINTEXC(i + 1)之后添加以下行：if i == list[j]: return nothing，并将while i > 0改为while i >= 0；(2) 使用以下伪代码以及修改后的Shuffle方法：while True; list = []; for i in 0...n: AddItem(list, i); s=Shuffle(list); if s!=nothing: return s; end。- Ting (2021)[^14] 展示了如何通过*哈希表*来降低洗牌的空间复杂度：修改
Shuffle如下：(1) 在方法开始时创建一个哈希表；(2) 不是执行交换tmp = list[i]; list[i] = list[j]; list[j] = tmp，而是使用以下方法：list[i] = HGET(j,j); HSET(k,HGET(i,i)); if k==i: HDEL(i)，其中HSET(k,v)将哈希表中键为k的项设置为v；HGET(k,v)获取表中键为k的项，如果没有该项则返回v；HDEL(k)从表中删除键为k的项。（哈希表也可以是任何键/值映射结构，包括红黑树。这可以与注 4 结合使用来生成错排，只是将注 4 中的list[j]替换为HGET(j,j)。）
随机字符字符串
要生成随机字符字符串
- 准备字符串可以包含的字母、数字和其他字符的列表。示例稍后在本节中提供。
- 构建一个新字符串，其字符从该字符列表中随机选择。下面的伪代码演示了该方法。该方法以有放回的方式随机采样字符，并返回采样字符的列表。（如何将此列表转换为文本字符串取决于编程语言，超出了本页的范围。）该方法接受两个参数：
characterList是步骤 1 中的列表，stringSize是随机字符的数量。
METHOD RandomString(characterList, stringSize)
i = 0
newString = NewList()
while i < stringSize
// Choose a character from the list
randomChar = characterList[RNDINTEXC(size(characterList))]
// Add the character to the string
AddItem(newString, randomChar)
i = i + 1
end
return newString
END METHOD
以下是字符列表的示例
- 对于*字母数字字符串*，或字母和数字的字符串，字符可以是 Unicode 标准中提供的基本数字“0”到“9”（U+0030-U+0039，号 48-57）、基本大写字母“A”到“Z”（U+0041-U+005A，号 65-90）和基本小写字母“a”到“z”（U+0061-U+007A，号 96-122）。
- 对于十进制数字字符串，字符只能是基本数字。
- 对于十六进制数字（十六进制）字符串，字符可以是基本数字以及基本字母“A”到“F”或“a”到“f”（不能同时是）。
注释
随机采样伪代码
以下伪代码实现了两种方法
RandomKItemsFromFile实现水库抽样；它从无限大小的文件（file）中选择最多k个随机项。虽然伪代码引用了文件和行，但该技术适用于从大小未知的的数据集或列表中一次检索一项。在伪代码中，ITEM_OUTPUT(item, thisIndex)是一个占位符，用于返回要存储在列表中的项的代码；这可以包括项的值、从 0 开始的索引，或两者兼有。RandomKItemsInOrder返回一个列表，其中最多包含k个来自给定列表（list）的随机项，且顺序与它们在列表中的出现顺序相同。它基于（Devroye 1986）[^19]，第 620 页介绍的技术。
METHOD RandomKItemsFromFile(file, k)
list = NewList()
j = 0
index = 0
while true // until a value is returned
// Get the next line from the file
item = GetNextLine(file)
thisIndex = index
index = index + 1
// If the end of the file was reached, break
if item == nothing: break
// NOTE: The following line is OPTIONAL
// and can be used to choose only random lines
// in the file that meet certain criteria,
// expressed as MEETS_CRITERIA below.
// ------
// if not MEETS_CRITERIA(item): continue
// ------
if j < k // phase 1 (fewer than k items)
AddItem(list, ITEM_OUTPUT(item, thisIndex))
j = j + 1
else // phase 2
j = RNDINT(thisIndex)
if j < k: list[j] = ITEM_OUTPUT(item, thisIndex)
end
end
// NOTE: Shuffle at the end in case k or
// fewer lines were in the file, since in that
// case the items would appear in the same
// order as they appeared in the file
// if the list weren't shuffled. This line
// can be removed, however, if the order of
// the items in the list is unimportant.
if size(list)>=2: Shuffle(list)
return list
end
METHOD RandomKItemsInOrder(list, k)
n = size(list)
// Special case if k is 1
if k==1: return [list[RNDINTEXC(n)]]
i = 0
kk = k
ret = NewList()
while i < n and size(ret) < k
u = RNDINTEXC(n - i)
if u <= kk
AddItem(ret, list[i])
kk = kk - 1
end
i = i + 1
end
return ret
END METHOD
示例
- 从
n个项的列表中移除k个随机项（list）等同于通过RandomKItemsInOrder(list, n - k)生成一个新列表。- 过滤：如果应用程序需要反复采样同一个列表（带放回或无放回），但仅限于该列表的某些项，它可以创建一个要从中采样的项列表（或指向这些项的索引列表），然后从新列表中采样。[^20] 然而，对于大小无限或非常大的列表，这效果不佳。
拒绝采样
拒绝采样是一种生成满足某些要求的随机内容的简单而灵活的方法。实现拒绝采样
- 通过任何方法和任何分布和范围生成随机内容（例如数字或文本字符串）。
- 如果内容不满足预定标准，则转到步骤 1。
示例标准包括检查—
- 以这种方式生成的数字—
- 不小于最小阈值（*左截断*），
- 不大于最大阈值（*右截断*），
- 是素数，
- 可被某些数字整除或不被某些数字整除，
- 不是采样方法最近选择的数字之一，
- 尚未被选中（借助哈希表、红黑树或类似结构），
- 未连续出现超过期望的次数，或
- 不在数字“黑名单”中，
- 以此方式生成的点是否与之前的随机点足够远（借助 KD 树或类似结构），
- 以此方式生成的点是否位于简单或复杂的形状内，
- 以此方式生成的文本字符串是否与正则表达式匹配，或
- 上述两项或多项标准。
（KD 树、哈希表、红黑树、素数测试算法和正则表达式超出了本文档的范围。）
注释
- 拒绝采样的运行时间取决于接受率，即采样器接受采样输出的频率，而不是拒绝它。总的来说，该比率是可接受的输出数量除以总输出数量。
- 所有拒绝采样策略都有被拒绝数据的可能性，因此它们都有*可变运行时间*（事实上，它们可能无限期运行）。但是，图形处理单元（GPU）和其他同时运行多个任务的设备在所有任务都完成工作时工作效果更好。如果它们都实现了拒绝采样技术，由于其可变运行时间，这种情况是不可能发生的。如果拒绝采样器的每次迭代拒绝率很低，一种解决方案是让每个任务运行一次采样器，并使用自己的“随机数源”（例如，来自其自己的 PRNG 生成的数字），然后采用第一个未被任务如此拒绝的采样数字（这可能以非常低的速率失败）。[^21]
随机游走
随机游走是随时间变化的随机过程。随机游走的简单形式涉及随机选择一个改变游走状态的数字。下面的伪代码生成一个具有 *n* 步的随机游走，其中 STATEJUMP() 是添加到当前状态的下一个数字（参见本节后面的示例）。
METHOD RandomWalk(n)
// Create a new list with an initial state
list=[0]
// Add 'n' new numbers to the list.
for i in 0...n: AddItem(list, list[i] + STATEJUMP())
return list
END METHOD
注释
- 白噪声过程是通过创建一组以相同方式生成的独立随机变量的列表来模拟的。此类过程通常模拟不依赖于时间或当前状态的时间行为。一个例子是
ZeroOrOne(px,py)（用于模拟*伯努利过程*，其中每个数字要么以px/py的概率为 1，要么为 0）。- 一个有用的参考是 De Bruyne et al. (2021)[^22]。
示例
- 如果
STATEJUMP()是RNDINT(1) * 2 - 1，则随机游走生成的数字每次都比上一个数字增加 -1 或 1，随机选择。- 如果
STATEJUMP()是ZeroOrOne(px,py) * 2 - 1，则随机游走生成的数字，每次增加 1 的概率为px/py，否则为 -1。- 二项过程：如果
STATEJUMP()是ZeroOrOne(px,py)，则随机游走以px/py的概率推进状态。
随机日期和时间
可以使用类似以下的伪代码，在两个日期-时间（date1、date2）之间生成*随机日期-时间*：dtnum1 = DATETIME_TO_NUMBER(date1); dtnum2 = DATETIME_TO_NUMBER(date2); num = RNDINTRANGE(date1, date2); result = NUMBER_TO_DATETIME(num)。
在该伪代码中，DATETIME_TO_NUMBER 和 NUMBER_TO_DATETIME 分别将日期-时间转换为数字或从数字转换回日期-时间，精确到所需粒度，例如月、日或小时粒度（此类转换的细节取决于日期-时间格式，超出了本文档的范围）。可以使用任何其他随机选择策略代替 RNDINTRANGE(date1, date2)。
统计检验中的随机化
统计检验使用洗牌和*自助法*来帮助通过随机化对数据得出结论。
- 洗牌用于数据集中每个项属于多个互斥组之一的情况。这里，通过洗牌原始数据集并按顺序对洗牌后的数据集中的每个项重新分组来生成一个或多个*模拟数据集*，使得模拟数据集的每个组中的项数与原始数据集相同。
- 自助法是一种创建现有数据集的一个或多个随机样本（模拟数据集）的方法，其中每个样本中的项是*带放回*随机选择的（参见从列表中随机选择一项）。（这样，每个随机样本可以包含重复项。）另请参阅（Brownlee 2018）[^23]。
创建模拟数据集后，将计算一个或多个统计量（例如均值）以及原始数据集的每个模拟数据集，然后将模拟数据集的统计量与原始数据集的统计量进行比较（此类比较超出了本文档的范围）。
马尔可夫链
马尔可夫链模拟一个或多个*状态*（例如，单个字母或音节），并存储从一个状态转移到另一个状态的概率（例如，“b”到“e”的概率为 20%，或“b”到“b”的概率为 1%）。因此，每个状态都可以看作是具有其自身的、针对每个相关状态转换的*权重*列表（参见“有放回加权选择”）。例如，用于生成*“可发音”单词*或与自然语言单词类似的单词的马尔可夫链可以包含单词开始和结束的“start”和“stop”状态。
一种称为*过去耦合*（Propp 和 Wilson 1996）[^24] 的算法可以从马尔可夫链的*平稳分布*（即链的稳态）中采样状态，通过从不同状态启动多个链并运行它们直到它们在同一时间达到同一状态。但是，除非采取预防措施（Fill 1998）[^25]，否则过早停止算法可能会引入偏差。如果链具有有限数量的状态且不是周期的，并且每个状态都可以从任何其他状态到达，则该算法可以正确工作。
以下伪代码实现了过去耦合。在该方法中，StateCount 是马尔可夫链的状态数，UPDATE(chain, state, random) 在给定当前状态和随机变量的情况下将马尔可夫链转移到下一个状态，RANDOM() 生成 UPDATE 所需的一个或多个随机变量。
METHOD CFTP(chain)
states=[]
numstates=StateCount(chain)
done=false
randoms=[]
while not done
// Start multiple chains at different states. NOTE:
// If the chain is monotonic (meaning the states
// are ordered and, whenever state A is less
// than state B, A's next state is never higher than
// B's next state), then just two chains can be
// created instead, starting
// at the first and last state, respectively.
for i in 0...numstates: AddItem(states, i)
// Update each chain with the same randomness
AddItem(randoms, RANDOM())
for k in 0...size(randoms):
for i in 0...numstates: states[i]=
UPDATE(chain, states[i], randoms[size(randoms)-1-k])
end
// Stop when all states are the same
fs=states[0]
done=true
for i in 1...numstates: done=(done and states[i]==fs)
end
return states[0] // Return the steady state
END METHOD
注意： 离散相型分布由马尔可夫链、起始状态和结束状态组成。它模拟了马尔可夫链从起始状态转移到结束状态所需的（随机）步数减 1。
随机图
图是点及其之间连接的列表。点称为*顶点*，连接称为*边*。
表示图的便捷方法是*邻接矩阵*。这是一个 *n*×*n* 矩阵，有 *n* 行和 *n* 列（表示图中有 *n* 个顶点）。对于简单图，邻接矩阵仅包含 1 和 0 — 对于第 *r* 行和第 *c* 列的单元格，1 表示存在一条从顶点 *r* 指向顶点 *c* 的边，0 表示不存在。
在本节中，Zeros(n) 创建一个 *n*×*n* 的零矩阵（例如，一个由 n 个列表组成的列表，每个列表包含 n 个零）。
以下方法生成一个遵循 G(*n*, *p*) 模型（也称为*Gilbert 模型*（Gilbert 1959）[^26]）的随机 *n* 顶点图，其中每条边以 px/py 的概率绘制（Batagelj 和 Brandes 2005）[^27]。
METHOD GNPGraph(n, px, py)
graph=Zeros(n)
for i in 2...n
j = i
while j > 0
j = j - 1 - min(NegativeBinomialInt(1, px, py), j - 1)
if j > 0
// Build an edge
graph[i][j]=1
graph[j][i]=1
end
end
end
return graph
END METHOD
还有其他类型的图，包括*Erdős–Rényi 图*（给定一个 *n*×*n* 邻接矩阵，无放回地均匀选择 m 条随机边）、Chung–Lu 图、优先连接图等。例如，*网格图*是矩形网格形式的图，其中顶点是网格矩形的角，边是其边。随机*迷宫*是网格图的随机生成树（Costantini 2020）[^28]。Penschuck et al. (2020)[^29] 提供了随机图生成技术的调查。
关于随机变量排序的说明
总的来说，对随机变量进行排序与对任何其他数据进行排序没有区别。（排序算法超出了本文档的范围。）[^30]
通用非均匀分布
一些应用需要选择随机值，使得其中一些值被选中的几率更大（*非均匀*分布）。本节中的大多数技术都展示了如何使用均匀随机整数方法来生成此类随机值。
加权选择
加权选择方法从集合中生成随机项或值，其中每个项或值被选中的概率不同。加权选择有几种类型。
有放回加权选择
第一种称为加权选择*带放回*（可以看作是抽取一个球，然后放回）或*分类分布*，其中选择每个项的概率在选择项时不会改变。在下面的伪代码中
WeightedChoice接受一个权重列表weights（非负整数）并返回该列表中权重的*索引*。权重越大，其索引被选中的几率越大。CumulativeWeightedChoice接受一个 N 个*累积权重*的列表；它们从 0 开始，下一个权重不小于前一个。返回一个区间 [0, N - 1) 中的数字。NormalizeRatios计算一个整数列表，其比例与给定有理数列表（形式为x/y的数字）的比例相同。这对于将有理数权重转换为整数权重以供WeightedChoice使用非常有用。gcd(a, b)是两个数之间的最大公约数（其中gcd(0, a) = gcd(a, 0) = a，只要a >= 0）。
METHOD WChoose(weights, value)
// Choose the index according to the given value
runningValue = 0
for i in 0...size(weights) - 1
if weights[i] > 0
newValue = runningValue + weights[i]
// NOTE: Includes start, excludes end
if value < newValue: break
runningValue = newValue
end
end
// Should not happen with integer weights
return error
END METHOD
METHOD WeightedChoice(weights)
return WChoose(weights,
RNDINTEXC(Sum(weights)))
END METHOD
METHOD CumulativeWeightedChoice(weights)
if size(weights)==0 or weights[0]!=0: return error
value = RNDINTEXC(weights[size(weights) - 1])
// Choose the index according to the given value
for i in 0...size(weights) - 1
// Choose only if difference is positive
if weights[i] < weights[i+1] and
weights[i]>=value: return i
end
return 0
END METHOD
METHOD NormalizeRatios(ratios)
prod=1; gc=0
for r in ratios: prod*=r[1] // Multiply denominators
weights=[]
for r in ratios
rn = floor(r * prod)
gc = gcd(rn, gc); AddItem(weights, rn)
end
if gc==0: return weights
for i in 0...size(weights): weights[i]=floor(weights[i]/gc)
return weights
END METHOD
以下是实现 WeightedChoice 的各种方法。其中许多需要使用特殊的数据结构。
| 算法 | 注释 |
|---|---|
| 线性搜索 | 参见上面 WeightedChoice 的伪代码。 |
| 带累积权重的线性搜索 | 计算累积权重列表（也称为*累积分布表*或*CDT*），然后随机生成一个小于（原始）权重总和（如果权重是整数，则应为整数）的数字，然后对新列表进行线性扫描以查找第一个累积权重超过生成数字的项。 |
| 快速加载骰子滚筒 (Saad et al., 2020a)[^31]。 | 仅使用整数权重，并使用随机比特（“公平硬币”）进行采样。该采样器每次采样所需的随机比特数平均接近最优值（相差最多 6 比特）。 |
| （Saad et al., 2020b）[^32] 中介绍的采样器 | 仅使用整数权重，并使用随机比特进行采样。只要权重之和为 2k 或 2k - 2m 的形式，采样器每次采样所需的随机比特数平均接近最优值（相差最多 2 比特）。 |
| 拒绝采样 | 给定一个权重列表（weights）包含 n 个权重：(1) 找到最高权重并称之为*max*；(2) 将 i 设置为 RNDINT(n - 1)；(3) 以 weights[i]/max 的概率（例如，对于整数权重，如果 ZeroOrOne(weights[i], max)==1），返回 i，否则转到步骤 2。（参见，例如，快速加载骰子滚筒论文的第 4 节，或 Tang 或 Klundert 的论文。weights[i] 也可以是一个“即时”计算权重的函数；在这种情况下，max 是所有 i 的 weights[i] 的最大值。）如果权重是 log-weight（即，每个权重是 ln(x)/ln(b)，其中 x 是原始权重，b 是对数底），则步骤 3 变为：“(3) 如果 Expo(ln(b)) > max - weights[i]（这以 pow(b, -(max - weights[i])) 的概率发生），返回 i，否则转到步骤 2。”如果权重是“硬币”，每枚硬币都有一个单独但未知的正面概率，该算法也称为*伯努利竞赛*（Dughmi et al. 2017）[^33]：(1) 将 i 设置为 RNDINT(n - 1)；(2) 抛掷硬币 i（第一个硬币是 0，第二个是 1，依此类推），如果返回 1 或正面，则返回 i，否则转到步骤 1。 |
| Bringmann 和 Panagiotou (2012/2017)[^34]。 | 展示了一种专为对排序的权重列表进行操作而设计的采样器。 |
| 别名方法（Walker 1977）[^35] | Michael Vose 的别名方法版本（Vose 1991）[^36] 在“飞镖、骰子和硬币：离散分布采样”中进行了描述。权重应该是无理数。 |
| （Klundert 2019）[^37] | 各种数据结构，重点介绍了它们如何支持权重更改。 |
| Bringmann–Larsen 紧凑型数据结构（Bringmann 和 Larsen 2013）[^38] | 如果权重之和很大，则使用拒绝采样，否则使用压缩结构。 |
| Hübschle-Schneider 和 Sanders (2019)[^39]。 | 并行加权随机采样器。 |
| （Tang 2019）[^40]。 | 提出了各种算法，包括二维和多级搜索、二分搜索（带累积权重）以及新的“平面”方法。 |
| “用有偏硬币掷出的骰子” | 给定一个求和为 1 且应为有理数的概率列表 probs：(1) 设置 cumu 为 1，i 为 0；(2) 以 probs[i]/cumu 的概率返回 i；(3) 从 cumu 中减去 probs[i]，然后将 i 加 1，然后转到步骤 2。正确性证明请参阅“飞镖、骰子和硬币”。如果 probs 中的每个概率都是“即时”计算的，这也称为顺序搜索；请参阅 Devroye (1986)[^19] 的第 10 章（但这通常不精确，除非所有涉及的概率都是有理数）。 |
| Knuth 和 Yao (1976)[^9] | 从概率的二进制展开（即它们是 0.bbbbbb... 的形式，其中 b 是 0 或 1）生成 DDG 二进制树。每次采样所需的随机比特数平均接近最优值（相差最多 2 比特）。这在 Devroye (1986)[^19] 的第 15 章练习 3.4.2 中有建议，在 *randomgen.py* 中实现为 discretegen 方法，并在（Devroye 和 Gravel 2020）[^11] 中进行了描述。只要能“即时”计算二进制展开，discretegen 就可以处理无理数概率（具有无限二进制展开）。 |
| Han 和 Hoshi (1997)[^41] | 使用累积概率作为输入，每次采样所需的随机比特数平均接近最优值（相差最多 3 比特）。也在（Devroye 和 Gravel 2020）[^11] 中进行了描述。 |
注释
- 加权选择算法作为二叉树遍历器。 就像
RNDINT算法一样（参见RNDINT部分），加权选择算法都可以描述为在 DDG 二叉树上随机行走。在这种情况下，概率不一定均匀，平均而言，该算法需要至少等于所有概率的*二元熵*之和的无偏随机比特数。例如，假设我们为四个整数 1、2、3、4 分配以下权重：3、15、1、2。这些权重的二元熵加起来为 0.4010... + 0.3467... + 0.2091... + 0.3230... = 1.2800...（因为权重之和为 21，3/21 的二元熵为0 - (3/21) * log2(3/21) = 0.4010...，其中log2(x) = ln(x)/ln(2)，依此类推）。因此，任何加权采样器平均需要至少 1.2800... 比特来生成概率与这些权重成比例的数字。[^10]RNDINT部分的“减少‘比特浪费’”说明也适用于此处。- 为获得最佳结果，传递给上表中算法的权重应首先*转换为整数*（例如，使用上面伪代码中的
NormalizeRatios），*或在指明时转换为有理数*。（显然，如果权重最初是无理数，则此转换将是有损的，算法将不精确，除非表中另有说明。）此外，在算法中使用浮点数可能会引入不必要的舍入误差，因此应避免使用此类数字。- Python 示例代码包含
WeightedChoice伪代码的一个变体，用于一次生成多个随机点。示例
- 假设我们有以下列表：
["apples", "oranges", "bananas", "grapes"]，以及以下weights：[3, 15, 1, 2]。“oranges”的权重为 15，“apples”的权重为 3。由于“oranges”的权重高于“apples”，因此使用WeightedChoice方法时，“oranges”的索引（1）被选中的概率高于“apples”的索引（0）。以下惯用语实现了使用该方法从列表中随机选择一项的方法：item = list[WeightedChoice(weights)]。- 如果
weights是以下累积权重列表：[0, 3, 18, 19, 21]，则可以使用CumulativeWeightedChoice而不是WeightedChoice来实现示例 1。- 分段常数分布。 假设使用示例 1 的权重，并且列表包含以下内容：
[0, 5, 10, 11, 13]（比权重多一项）。这表示四个区间：[0, 5)、[5, 10) 等。使用index = WeightedChoice(weights)选择一个随机索引（从而选择一个区间）。然后独立地，在该区间内均匀随机选择一个数字（例如，以下代码以这种方式选择一个随机整数：number = RNDINTEXCRANGE(list[index], list[index + 1])）。
无放回加权选择
加权选择*无放回*可以看作是*不*放回地抽取一个球。
以下是实现加权选择无放回的方法，其中每个项*最多只能被选中一次*。权重具有这样的性质：权重较高的项首先出现的几率更大。
- 使用
WeightedChoice选择随机索引。每次选择一个索引时，将该索引的权重设置为 0，以防止再次选中它。或者... - 为每个索引分配一个随机指数随机变量（速率等于该索引的权重，该权重必须是整数 1 或更大）。创建一个将每个数字分配给索引的对列表，然后按升序（从小到大）对该列表进行排序（参见下面的注 2）。排序后的索引列表将对应于无放回的加权选择。参见“无放回采样算法”，另请参阅（Efraimidis 2015）[^42]。（Efraimidis 和 Spirakis 2005）[^43] 描述了如何使用水库抽样（本质上是具有关联的
k个最低随机变量的*优先级队列*）在不存储每个项的情况下实现此方法。
注释
一些应用程序（特别是某些游戏）希望控制哪些随机结果出现，以使这些结果对用户来说“更公平”（例如，避免连续出现好结果或坏结果）。只要在所有结果都选完后重新填充列表，就可以使用无放回加权选择列表结果来实现此目的。但是，在关心信息安全的应用中应避免此类采样，包括当玩家或用户可以轻易猜测结果时，如果结果易于猜测，则可能会获得显著且不公平的优势。
ExpoNew(weight)使用速率weight创建一个“空的”指数随机变量，其比特尚未确定（参见“部分采样随机数”）。对变量进行排序依赖于通过ExpoLess(a, b)比较两个变量，它返回一个指数数（a）是否小于另一个（b），并根据需要构建两个变量的比特。对于不太精确的算法，将ExpoNew(weight)替换为ExpoApprox(weight)，并将ExpoLess(a, b)替换为a < b，其中ExpoApprox具有以下伪代码，该伪代码使用 von Neumann 的（1951）[^44] 算法：METHOD ExpoApprox(w); b=1000000; count=0; while true; y1=RNDINTEXC(b); y=y1; accept=true; while true; z=RNDINTEXC(b); if y<=z: break; if accept: accept=false; else: accept=true; y=z; end; if accept: count=count+y1; else: count=count+b; if accept: return MakeRatio(count, b*w); end; END METHOD。
不等概率采样
对于上一节中给出的方法，权重具有这样的性质：权重较高的项首先被选中，但每个项的权重不一定是从给定 n 项样本中包含该项的几率（*包含概率*）。以下方法从项列表中选择 n 个索引的随机样本（其权重是存储在列表 weights 中的整数），使得索引 k 在样本中的几率为 weights[k]*n/Sum(weights)。选定的索引不一定按随机顺序。该方法实现了“拆分方法”（Deville 和 Tillé 1998）[^45]。
METHOD InclusionSelect(weights, n)
if n>size(weights): return error
if n==0: return []
ws=Sum(weights)
wts=[]
items=[]
// Calculate inclusion probabilities
for i in 0...size(weights):
AddItem(wts,[MakeRatio(weights[i],ws)*n, i])
Sort(wts)
// Check for invalid inclusion probabilities
if wts[size(wts)-1][0]>1: return error
last=size(wts)-n
if n==size(wts)
for i in 0...n: AddItem(items,i)
return items
end
while true // until a value is returned
lamda=min(MakeRatio(1,1)-wts[last-1][0],wts[last][0])
if lamda==0: return error
if ZeroOrOne(lamda[0],lamda[1])
for k in 0...size(wts)
if k+1>ntotal-n:AddItem(items,wts[k][1])
end
return items
end
newwts=[]
for k in 0...size(wts)
newwt=(k+1<=last) ?
wts[k][0]/(1-lamda) : (wts[k][0]-lamda)/(1-lamda)
AddItem(newwts,[newwt,wts[k][1]])
end
wts=newwts
Sort(wts)
end
END METHOD
对于项列表大小未知且其权重可以“即时”计算的情况，请参阅（Chao 1982）[^46]。
分布混合
混合由两个或多个概率分布组成，每个分布都有单独的采样概率。要从混合物生成随机内容—
- 生成
index = WeightedChoice(weights)，其中weights是每个分布在混合物中被采样的相对概率列表，然后 - 根据
index的值，从相应的分布生成随机内容。
示例
- 一种混合物由三个六面虚拟骰子点数之和以及一个六面骰子点数组成，但有 80% 的几率掷一个六面虚拟骰子而不是三个。以下伪代码显示了如何采样这种混合物：
index = WeightedChoice([80, 20]); number = 0; if index==0: number = RNDINTRANGE(1,6); else: number = RNDINTRANGE(1,6) + RNDINTRANGE(1,6) + RNDINTRANGE(1,6)。- 从复杂形状（任意维度）中选择一个独立的均匀随机点，等同于从此复杂形状组成的更简单形状的混合物中进行采样（此处，
weights列表保存每个更简单形状的 n 维“体积”）。例如，一个简单的闭合二维多边形可以被三角剖分，或分解为三角形，然后可以对这些三角形的混合物进行采样。[^47]取一组不重叠的整数区间（例如 [0, 5]，[7, 8]，[20, 25]）。要从这些区间中选择一个独立的均匀随机整数
- 创建一个权重列表（
weights），为每个区间分配一个权重。每个区间赋予权重(mx - mn) + 1，其中mn是该区间的最小值，mx是其最大值。- 使用
WeightedChoice(weights)选择一个索引，然后生成RNDINTRANGE(mn, mx)，其中mn是相应区间的最小值，mx是其最大值。此方法可以应用于具有共同分母的有理数，通过将涉及的整数视为这些有理数的分子。例如，[0/100, 5/100]，[7/100, 8/100]，[20/100, 25/100]，其中分子与前面的示例相同。
- 在伪代码
index = WeightedChoice([80, 20]); list = [[0, 5], [5, 10]]; number = RNDINTEXCRANGE(list[index][0], list[index][1])中，以 80% 的几率选择 [0, 5) 中的随机整数，以 20% 的几率选择 [5, 10) 中的随机整数。- 超指数分布是*指数分布*的混合物，每个指数分布都有单独的权重和单独的速率参数。
随机变量的变换
可以通过组合和/或转换一个或多个随机变量以及/或丢弃其中一些来生成随机变量。
例如，Red Blob Games 的“概率与游戏：伤害掷骰”包含交互式图形，显示了最低值、最高值、丢弃最低值和重掷的掷骰机制的分数分布。[^48] 这些以及类似的分布可以概括如下。
生成一个或多个随机变量（数字），每个变量具有不同的概率分布，然后
- 最高值：选择最高生成的变量。
- 丢弃最低值：加总所有生成的变量，除了最低的。
- 重掷最低值：加总所有生成的变量，除了最低的，然后添加一个由不同概率分布随机生成的数字。
- 最低值：选择最低生成的数字。
- 丢弃最高值：加总所有生成的变量，除了最高的。
- 重掷最高值：加总所有生成的变量，除了最高的，然后添加一个由不同概率分布随机生成的数字。
- 总和：将所有生成的随机变量相加。
- 平均值：将所有生成的随机变量相加，然后将总和除以随机变量的数量。
- 几何变换：将随机变量视为一个n维点，然后对该点应用几何变换，例如旋转或其他仿射变换[^49]。
如果概率分布相同，则策略 1 到 3 会使较大的数字具有更高的概率，而策略 4 到 6 则使较小的数字具有更高的概率。
注意：策略 4 的变体——例如，选择第二、第三或第 n 个最小的数字——分别称为第二、第三或第 n 个阶统计量分布。
示例
- 习语
min(RNDINTRANGE(1, 6), RNDINTRANGE(1, 6))获取两个六面骰子结果中的最小值（策略 4）。由于这种方法，1 出现的几率大于 6。- 习语
RNDINTRANGE(1, 6) + RNDINTRANGE(1, 6)获取两个六面骰子的结果（另请参阅“骰子”）（策略 7）。- 二项分布对由
ZeroOrOne(px,py)（策略 7）生成的n个数字的总和进行建模（参见“二项分布”）。- 泊松二项分布对由每个数字都有单独的概率为 1 而不是 0 的
n个数字的总和进行建模（策略 7）。给定probs，即n个概率作为有理数的列表，伪代码如下：for i in 0...n: x=x+ZeroOrOne(probs[i][0],probs[i][1]); return x。截断随机变量。这是变换随机变量的一个示例。要生成一个截断随机变量，请像往常一样随机生成一个数字，然后——
- 如果该数字小于最小阈值，则使用最小阈值（左截尾），和/或
- 如果该数字大于最大阈值，则使用最大阈值（右截尾）。
截断随机变量的一个例子是
min(200, RNDINT(255))。- 复合泊松分布对由相同方式随机选择的 n 个数字的总和进行建模，其中 n 服从泊松分布（例如，
n = PoissonInt(10, 1)表示平均 10 个数字）（策略 7，总和）。- Pólya–Aeppli 分布是一种复合泊松分布，其中数字由
NegativeBinomial(1, 1-p)+1生成，其中p是固定的。
特定非均匀分布
本节包含有关一些最常见的非均匀概率分布的信息。
骰子
以下方法生成虚拟骰子掷出的随机结果。它接受三个参数：骰子数量 (dice)、每个骰子的面数 (sides) 以及添加到结果中的数字 (bonus)（该数字可以是负数，但如果结果大于 0，则该方法的返回值为 0）。另请参阅 Red Blob Games 的“概率与游戏：伤害掷骰”。
METHOD DiceRoll(dice, sides, bonus)
if dice < 0 or sides < 1: return error
ret = 0
for i in 0...dice: ret=ret+RNDINTRANGE(1, sides)
return max(0, ret + bonus)
END METHOD
示例：掷出——
- 四个六面虚拟骰子（“4d6”）的结果是
DiceRoll(4,6,0)，- 三个十面虚拟骰子，加上 4（“3d10 + 4”）的结果是
DiceRoll(3,10,4)，以及- 两个六面虚拟骰子，减去 2（“2d6 - 2”）的结果是
DiceRoll(2,6,-2)。
二项分布
二项分布使用两个参数：trials 和 p。此分布对固定数量的独立试验（等于 trials）中的成功次数进行建模，每次试验的成功概率相同（等于 p，其中 p <= 0 表示从不成功，p >= 1 表示总是成功，p = 1/2 表示成功或失败的几率相等）。在本文件中，Binomial(trials, p) 是具有给定参数的二项随机变量。
此分布有一个简单的实现：count = 0; for i in 0...trials: count=count+ZeroOrOne(px, py)。但对于大量的试验，这可能会非常慢。
以下伪代码实现了一个精确的该分布采样器，并根据（Farach-Colton and Tsai 2015）[^50] 进行了一些优化。（另一个精确采样器在（Bringmann et al. 2014）[^51] 中给出，并在我的“随机化杂项观察”中进行了描述。）这里，参数 p 表示为比率 px/py。
METHOD BinomialInt(trials, px, py)
if trials < 0: return error
if trials == 0: return 0
// Always succeeds
if mx: return trials
// Always fails
if p <= 0.0: return 0
count = 0
ret = 0
recursed = false
if py*2 == px // Is half
if i > 200
// Divide and conquer
half = floor(trials / 2)
return BinomialInt(half, 1, 2) + BinomialInt(trials - half, 1, 2)
else
if rem(trials,2)==1
count=count+RNDINT(1)
trials=trials-1
end
// NOTE: This step can be made faster
// by precalculating an alias table
// based on a list of n + 1 binomial(1/2)
// weights, which consist of n-choose-i
// for every i in [0, n], and sampling based on
// that table (see Farach-Colton and Tsai).
for i in 0...trials: count=count+RNDINT(1)
end
else
// Based on proof of Theorem 2 in Farach-Colton and Tsai.
// Decompose px/py into its base-2 digits.
pw = MakeRatio(px, py)
pt = MakeRatio(1, 2)
while trials>0 and pw>0
c=BinomialInt(trials, 1, 2)
if pw>=pt
count=count+c
trials=trials-c
pw=pw-pt
else
trials=c
end
pt=pt/2 // NOTE: Not rounded
end
end
if recursed: return count+ret
return count
END METHOD
注意：如果
px/py为1/2，则二项分布模拟“抛 N 枚硬币，然后计算正面朝上的数量”的任务，随机总和被称为汉明距离（将每次试验视为一个“位”，成功时设为 1，失败时设为 0）。如果px为1，则此分布模拟“掷n枚py面骰子，然后计算显示数字 1 的骰子数量”的任务。
负二项分布
在本文件中，负二项分布模拟在固定数量的成功试验（successes）之前发生的失败试验的数量。每次试验都是独立的，成功概率为 px/py（其中 0 表示从不成功，1 表示总是成功）。以下是一个朴素的实现；另请参阅几何分布的注意事项，它是此分布的一个特例。
METHOD NegativeBinomialInt(successes, px, py)
// Needs to be 0 or greater; px must not be 0
if successes < 0 or px == 0: return error
if successes == 0 or px >= py: return 0
total = 0
count = 0
while total < successes
if ZeroOrOne(px, py) == 1: total = total + 1
else: count = count + 1
end
return count
END METHOD
如果 successes 是非整数，则该分布通常称为Pólya 分布。在这种情况下，可以使用以下伪代码进行采样（Heaukulani and Roy 2019）[^52]
METHOD PolyaInt(sx, sy, px, py)
isinteger=rem(sx,sy)==0
sxceil=ceil(sx/sy)
while true // until a value is returned
w=NegativeBinomialInt(sxceil, px, py)
if isinteger or w==0: return w
tmp=MakeRatio(sx,sy)
anum=tmp
for i in 1...w: anum=anum*(tmp+i)
tmp=sxceil
aden=tmp
for i in 1...w: aden=aden*(tmp+i)
a=anum/aden
if ZeroOrOne(a[0], a[1])==1: return w
end
END METHOD
几何分布
几何分布是 successes = 1 的负二项分布。在本文件中，几何随机变量是指在发生一次成功之前发生的失败次数。例如，如果 p 为 1/2，则几何分布模拟“抛硬币直到出现反面，然后计算正面朝上的次数”的任务。作为几何分布的一个独特属性，给定 n 次试验失败，新的失败试验次数具有相同的分布（其中 n 是大于 0 的整数）。
注释
- 负二项分布和几何分布在不同的文献中有不同的定义。例如，Mathematica 的定义排除了最后一次成功，而（Devroye 1986, p. 498）[^19] 中的定义则包含了它。有些文献可能将负二项数定义为 N 次失败之前的成功次数，而不是反之。
- 有界几何随机变量是 n（大于 0 的整数）或几何随机变量，取两者中较小者。几何分布和有界几何分布的精确且高效的采样器在（Bringmann and Friedrich 2013）[^53] 中给出，并在我的“随机化杂项观察”中进行了描述。）
指数分布
指数分布使用一个称为 λ 的参数，即率或逆尺度。通常，λ 是在给定时间跨度（例如，给定的一天或一年）内发生的给定类型独立事件的概率，随机结果是直到该事件发生的那个时间跨度的数量。通常，λ 等于 1，或 1/1。1/λ 是尺度（均值），通常是同一类型两个独立事件之间的平均等待时间。
在本文件中，Expo(lamda) 是具有率 lamda 的指数分布随机变量。有关任意精度采样指数随机变量的算法，请参阅“部分采样随机数”。
泊松分布
泊松分布使用参数 mean（也称为 λ）。λ 是每固定单位时间或空间（例如，每天、每小时或每平方公里）的给定类型独立事件的平均数量。泊松分布的数字是该单元内此类事件的数量。
在本文件中，如果 mean 大于 0，则 Poisson(mean) 是泊松分布的数字，如果 mean 为 0，则为 0。
方法 PoissonInt 生成均值为 mx/my 的泊松随机变量，其中 Poisson1 方法使用 Duchon 和 Duvignau (2016) 的算法[^54]。
METHOD Poisson1()
ret=1; a=1; b=0
while true // until this method returns
j=RNDINT(a)
if j<a and j<b: return ret
if j==a: ret=ret+1
else
ret=ret-1; b=a+1
end
a=a+1
end
END METHOD
METHOD PoissonInt(mx, my)
if my == 0: return error
if mx == 0 or (mx < 0 and my < 0) or (mx > 0 and my < 0): return 0
r=0
while mx>=my
r=r+Poisson1(); mx=mx-my
end
if mx>0
// At this point, mx/my < 1, so sum a Poisson number
// of coin flips with heads prob. mx/my; see Devroye 1986, p. 487
r=r+BinomialInt(Poisson1(), mx, my)
end
return r
END METHOD
注意：要生成具有单独均值的
n个独立泊松随机变量的总和，请生成均值为这些均值总和的泊松随机变量（参见（Devroye 1986）[^19]，第 501 页）。例如，要生成均值为 1/1000000 的 1000 个独立泊松随机变量的总和，只需生成PoissonInt(1, 1000)（因为 1/1000000 * 1000 = 1/1000）。
Pólya–Eggenberger 分布
假设从标记为 1 或 0 的项目集合中随机抽取项目，并且在抽取一个项目后，将其放回并将相同标签的 m 个项目添加进去。那么
- Pólya–Eggenberger 分布模拟以这种方式抽取的标记为
1的项目的数量。 - 逆 Pólya–Eggenberger 分布模拟在抽取
successes个标记为1的项目之前抽取的标记为0的项目的数量。
（Johnson and Kotz 1969）[^55]。在下面的方法中，trials 是随机抽取的项目数，ones 是集合中标记为 1 的项目数，count 是集合中标记为 1 或 0 的项目数，m 是每次抽取后添加的项目数（或不放回抽样为 −1），successes 是抽取到的标记为 1 的项目的数量。
METHOD PolyaEggenberger(trials, ones, count, m)
if ones < 0 or count < 0 or trials < 0 or
ones > count or trials > count
return error
end
if ones == 0: return 0
zeros=count-ones
ret=0
for i in 0...trials
if zeros==0 or ZeroOrOne(ones,zeros)==1
ones=ones+m
ret=ret+1
else: zeros=zeros+m
end
return ret
END METHOD
METHOD InversePolyaEggenberger(successes, ones, count, m)
if ones <= 0 or count < 0 or successes < 0 or
ones > count or successes > count
return error
end
zeros=count-ones
ret=0; trials=0
while ret<successes
if zeros==0 or ZeroOrOne(ones,zeros)==1
ones=ones+m
ret=ret+1
else: zeros=zeros+m
trials=trials+1
end
return trials-successes
END METHOD
注释
- 超几何分布是一种 Pólya–Eggenberger 分布，其中
m=-1。例如，在 52 张盎格鲁-美式扑克牌的牌堆中，有 12 张是人头牌（J、Q 或 K）。洗牌并抽取七张牌后，以这种方式抽取的人头牌数量遵循超几何分布，其中trials为 7，ones为 12，count为 52，m为 -1。- 负超几何分布是一种逆 Pólya–Eggenberger 分布，其中
m=-1。
具有给定正和的随机整数
以下伪代码展示了如何生成具有给定正总和的 n 个随机整数，并以随机顺序排列（具体来说，是该总和的均匀随机选择的分区，允许重复且顺序随机）。（此算法在（Smith and Tromble 2004）[^56] 中介绍。）在以下伪代码中——
- 方法
PositiveIntegersWithSum返回n个大于 0 且总和为total的整数，并以随机顺序排列。 - 方法
IntegersWithSum返回n个大于等于 0 且总和为total的整数，并以随机顺序排列。 Sort(list)按升序对list中的项进行排序（但排序算法超出本文档范围）。
METHOD PositiveIntegersWithSum(n, total)
if n <= 0 or total <=0: return error
ls = [0]
ret = NewList()
while size(ls) < n
c = RNDINTEXCRANGE(1, total)
found = false
for j in 1...size(ls)
if ls[j] == c
found = true
break
end
end
if found == false: AddItem(ls, c)
end
Sort(ls)
AddItem(ls, total)
for i in 1...size(ls): AddItem(ret,
ls[i] - ls[i - 1])
return ret
END METHOD
METHOD IntegersWithSum(n, total)
if n <= 0 or total <=0: return error
ret = PositiveIntegersWithSum(n, total + n)
for i in 0...size(ret): ret[i] = ret[i] - 1
return ret
END METHOD
注释
- 要以随机顺序生成具有给定正平均值
avg（整数）的N个随机整数，请生成IntegersWithSum(N, N * avg)。- 要生成具有给定正总和
sum（整数）的min或更大的N个随机整数，并以随机顺序排列，请生成IntegersWithSum(N, sum - N * min)，然后将min添加到每个以此方式生成的数字中。Python 示例代码实现了一种有效的方法来生成此类整数，前提是每个整数都不能超过给定的最大值；该算法归功于 John McClane 在 Stack Overflow 上的一个答案（questions/61393463）。- 要生成总和为
tx/ty的N个随机有理数，请调用IntegersWithSum(N, tx * ty * x)或PositiveIntegersWithSum(N, tx * ty * x)（取决于情况）（其中x是所需的精度，表示为整数，例如pow(2, 32)或pow(2, 53)，以便结果的精度达到1/x或更低），然后对于结果中的每个数字c，将其转换为MakeRatio(c, tx * ty * x) * MakeRatio(tx, ty)。
多项分布
多项分布使用两个参数：trials 和 weights。它模拟在给定数量的试验（trials）中，多个互斥事件发生的次数，其中每个事件都有单独的发生概率（以 weights 列表的形式给出）。
一个朴素的实现是，将一个列表填充为与 weights 相同数量的零，然后对于每次试验，选择 index = WeightedChoice(weights)，并将列表中的选定 index 处的项加 1。生成的列表遵循多项分布。下面的伪代码显示了（Durfee et al., 2018, Corollary 45）[^57] 中建议的一个优化，但它假设所有权重都是整数。
METHOD Multinomial(trials, weights)
if trials < 0: return error
// create a list of successes
list = []
ratios = []
sum=Sum(weights)
for i in 0...size(weights): AddItem(ratios,
MakeRatio(weights[i], sum))
end
for i in 0...size(weights)
r=ratios[i]
b=BinomialInt(t,r[0],r[1])
AddItem(list, b)
trials=trials-b
if trials>0: for j in range(i+1,
len(weights)): ratios[j]=ratios[j]/(1-r)
end
return list
END METHOD
实数随机化
本节描述了使用随机实数而非仅随机整数的随机化方法。这些包括随机有理数、定点数和浮点数。
但是，只要有可能，应用程序应使用随机整数而不是其他随机实数。这是因为
- 任何计算机都无法在两个数字之间的所有实数中进行选择，因为它们有无穷多个。
- 处理整数的算法比处理其他实数（尤其是浮点数）的算法更具可移植性。[^58] 整数算法更容易控制其精度级别。
- 对于可能关心可重现的“随机”数字（单元测试、模拟、机器学习等）的应用程序，使用非整数（尤其是浮点数）可能会使方法从一次运行到另一次运行或跨计算机进行重现的任务复杂化。
本节中的方法不应用于信息安全目的的随机抽样，即使有安全的“随机数源”。请参阅附录中的“安全注意事项”。
均匀随机实数
本节定义了一种方法，即 RNDRANGEMinMaxExc(a, b)，用于生成在 a 和 b 之间（不包含端点）的独立“均匀”随机实数。[^59]
本节展示了如何为定点、有理数和浮点数实现此方法。然而，这三种格式都使用预先确定且固定的精度。其他随机实数格式没有这个限制，包括部分采样随机数和“构造实数”或“递归实数”（Boehm 2020）[^60]。
对于定点数格式
对于表示 1/n 倍数的定点数格式，此方法很简单。以下代码返回一个表示定点的整数。在下面的方法中（以及在注释中），fpa 和 fpb 是生成的定点的边界，它们是表示定点的整数（使得 fpa = a * n 和 fpb = b * n）。例如，如果 n 为 100，要生成开区间 (6.35, 9.96) 中的一个数字，请生成 RNDRANGEMinMaxExc(6.35, 9.96) 或 RNDINTRANGE(635 + 1, 996 - 1)。
RNDRANGEMinMaxExc(a, b):RNDINTRANGE(fpa + 1, fpb - 1)，如果fpa >= fpb 或 a == fpb - 1则报错。但如果a为 0 且b为 1：(RNDINT(n - 2) + 1)或(RNDINTEXC(n - 1) + 1)。
注意：下面提供了在不同区间采样定点的其他方法，但本文档其余部分未使用。
RNDRANGE(a, b)，区间 [a,b]:RNDINTRANGE(fpa, fpb)。但如果a为 0 且b为 1：RNDINT(n)。RNDRANGEMinExc(a, b)，区间 (a,b]:RNDINTRANGE(fpa + 1, fpb)，如果fpa >= fpb则报错。但如果a为 0 且b为 1：(RNDINT(n - 1) + 1)或(RNDINTEXC(n) + 1)。RNDRANGEMaxExc(a, b)，区间 [a,b):RNDINTEXCRANGE(fpa, fpb)。但如果a为 0 且b为 1：RNDINTEXC(n)或RNDINT(n - 1)。
对于有理数格式
有理数是有理数的比值。如果该有理数的分子为 n（必须大于等于 1），则使用上一节生成其分子，以便该有理数是 1/n 的倍数。
对于浮点数格式
对于表示 FPSign * s * FPRADIXe 形式的数字的浮点数格式 [^61]，以下伪代码实现了 RNDRANGEMinMaxExc(lo, hi)。在伪代码中
MINEXP是数字在浮点格式中可以具有的最低指数。对于 IEEE 754 binary64 格式（Javadouble），MINEXP = -1074。对于 IEEE 754 binary32 格式（Javafloat），MINEXP = -149。FPPRECISION是浮点格式中的有效数字位数，无论格式是否以此方式存储。对于 binary64 为 53，对于 binary32 为 24。FPRADIX是浮点格式的数字基数。对于 binary64 和 binary32 为 2。FPExponent(x)返回数字x的e值，使得s中的位数等于FPPRECISION。如果x = 0或e小于MINEXP，则返回MINEXP。FPSignificand(x)返回s，即数字x的尾数。如果x = 0，则返回 0。除非FPExponent(x) == MINEXP，否则具有FPPRECISION位数字。FPSign(x)返回 -1 或 1，表示数字是正数还是负数。即使s为 0，也可以为 -1。
另请参阅（Downey 2007）[^62] 和Rademacher 浮点库。
METHOD RNDRANGEMinMaxExc(lo, hi)
if mn >= mx: return error
return RNDRANGEHelper(lo, hi)
END METHOD
METHOD RNDRANGEHelper(lo, hi)
losgn = FPSign(lo)
hisgn = FPSign(hi)
loexp = FPExponent(lo)
hiexp = FPExponent(hi)
losig = FPSignificand(lo)
hisig = FPSignificand(hi)
if lo > hi: return error
if losgn == 1 and hisgn == -1: return error
if losgn == -1 and hisgn == 1
// Straddles negative and positive ranges
// NOTE: Changes negative zero to positive
mabs = max(abs(lo),abs(hi))
while true // until a value is returned
ret=RNDRANGEHelper(0, mabs)
neg=RNDINT(1)
if neg==0: ret=-ret
if ret>=lo and ret<=hi: return ret
end
end
if lo == hi: return lo
if losgn == -1
// Negative range
return -RNDRANGEHelper(abs(lo), abs(hi))
end
// Positive range
expdiff=hiexp-loexp
if loexp==hiexp
// Exponents are the same
// NOTE: Automatically handles
// subnormals
s=RNDINTRANGE(losig, hisig)
return s*1.0*pow(FPRADIX, loexp)
end
while true // until a value is returned
ex=hiexp
while ex>MINEXP
v=RNDINTEXC(FPRADIX)
if v==0: ex=ex-1
else: break
end
s=0
if ex==MINEXP
// Has FPPRECISION or fewer digits
// and so can be normal or subnormal
s=RNDINTEXC(pow(FPRADIX,FPPRECISION))
else if FPRADIX != 2
// Has FPPRECISION digits
s=RNDINTEXCRANGE(
pow(FPRADIX,FPPRECISION-1),
pow(FPRADIX,FPPRECISION))
else
// Has FPPRECISION digits (bits), the highest
// of which is always 1 because it's the
// only nonzero bit
sm=pow(FPRADIX,FPPRECISION-1)
s=RNDINTEXC(sm)+sm
end
ret=s*1.0*pow(FPRADIX, ex)
if ret>=lo and ret<=hi: return ret
end
END METHOD
注释
下面提供了在不同区间采样“均匀”浮点数的其他方法，但本文档其余部分未使用。
RNDRANGE(mn, mx)，区间 [mn,mx]: 生成RNDRANGEHelper(mn, mx)。RNDRANGEMaxExc(mn, mx)，区间 [mx,mx): 如果mn >= mx，则返回错误。否则，在循环中生成RNDRANGEHelper(mn, mx)，直到生成一个非mx的数字为止。RNDRANGEMinExc(mn, mx)，区间 (mn,mx]: 如果mn >= mx，则返回错误。否则，在循环中生成RNDRANGEHelper(mn, mx)，直到生成一个非mn的数字为止。许多软件库通过将均匀随机整数乘以或除以常数来采样“均匀”实数。例如，从半开区间 [0, 1) 中“均匀”采样的函数通常实现为
RNDINTEXC(X) * (1.0/X)或RNDINTEXC(X) / X，其中 X 取决于软件库。[^63] 缺点是这样做并不一定能覆盖浮点格式在 (Goualard 2020)[^64] 范围内可表示的所有数字。另一个例子是，从半开区间 [a,b) 中“均匀”采样的函数通常实现为a + Math.random() * (b - a)，其中Math.random()是 [0, 1) 中的“均匀”随机浮点数；然而，这不仅有同样的缺点，而且在涉及浮点数时还有许多其他问题（Monahan 1985）[^65]。
蒙特卡罗采样：期望值、积分和优化
需要随机实数。
随机化是蒙特卡洛采样的核心。蒙特卡洛采样有三种主要用途：估计、积分和优化。
-
估计期望值。蒙特卡洛采样有助于估计采样分布的期望值（均值或“长期平均值”），或该分布采样值的函数。在本节中，此函数称为
EFUNC(x)，其中x是样本中的值之一。估计期望值的算法称为估计量。一个这样的估计量是采样n个值，将EFUNC(x)应用于每个采样值x，将值相加，然后除以n（见下文注释）。然而，这个估计量并不适用于所有分布，因为它们可能具有无限的期望值，并且它也不允许控制估计的误差。这个估计量被称为- 如果
EFUNC(x)是pow(x, n)，则为第n个样本原始矩（原始矩是n次幂的均值）。 - 如果
EFUNC(x)是x或pow(x, 1)，则为样本均值。 - 如果
EFUNC(x)是pow(x-m, n)，其中m是样本均值，则为第n个样本中心矩（中心矩是关于均值的矩）。 - （有偏）样本方差，即第二个样本中心矩。
- 如果满足某个条件则为
1，否则为0，则为样本概率。
蒙特卡洛估计量存在两个误差源：偏差和方差。如果一个估计量的期望值等于分布的期望值，那么它就是无偏的（偏差为 0）。例如，任何第
n个样本原始矩都是无偏估计量，前提是样本量至少为n，但样本方差不是无偏的，对于除第一个以外的任何样本中心矩，它都不是无偏的（Halmos 1946）[^66]。（“方差缩减”方法超出本文档范围。）估计量的均方误差等于方差加上偏差的平方。有关具有精度保证的蒙特卡洛估计量，请参阅“随机估计算法”。
- 如果
-
蒙特卡洛积分。这通常是蒙特卡洛估计的一个特例，它近似采样域上的多维积分；这里，
EFUNC(z)是要积分的函数，其中z是采样域中的随机选择点（例如，如果点在真实体积内则为 1，否则为 0）。 -
随机优化。它使用随机性来帮助找到具有一个或多个变量的函数的最小值或最大值；例如包括模拟退火和同时扰动随机逼近（另请参阅（Spall 1998）[^67]）。
注意：假设总体具有有限的均值和方差，则样本均值是均值的无偏估计量，但样本方差通常是除整个总体外的所有样本的方差的有偏估计量。以下伪代码返回一个包含样本均值和方差的无偏估计量（按此顺序）的两项列表，使用 J. D. Cook 提出的Welford 方法。此处计算的方差的平方根是许多 API 所称的标准差（例如，Python 的
statistics.stdev）。对于常规（有偏）样本方差，请将以下伪代码中的(size(list)-1)替换为size(list)。伪代码如下：if size(list)==0: return [0, 0]; if size(list)==1: return [list[0], 0]; xm=list[0]; xs=0; i=1; while i < size(list); c = list[i]; i = i + 1; cxm = (c - xm); xm = xm + cxm *1.0/ i; xs = xs + cxm * (c - xm); end; return [xm, xs*1.0/(size(list)-1)]。
点采样选择
需要随机实数。
已开发出各种方法来选择均匀分布的样本点，尤其用于蒙特卡洛方法。
在这些方法中，基数 N 的低差异序列（或拟随机点集或伪蒙特卡洛点集）是一个确定性点序列，与该分布的独立点相比，它对盒子上的均匀分布具有低差异。以下是示例
- 基数 N 的van der Corput 序列生成如下：对于序列中的每个非负整数索引，将其作为基数 N 的数字，然后将最低有效基数 N 位除以 N，下一个位除以 N2，再下一个除以 N3，依此类推，并将这些除法结果相加。
- Halton 序列是一组两个或多个具有不同素数基的 van der Corput 序列；给定索引处的 Halton 点具有等于该索引处 van der Corput 序列点的坐标。
- Roberts, M. 在“拟随机序列的非凡有效性”中，提出了一种基于黄金比例“广义”版本的低差异序列。
- Sobol 序列的解释请参见 S. Joe 和 F. Kuo 的“Sobol 序列生成器”。
低差异序列的点可以通过伪随机数生成器（或其他模拟“随机数源”的设备或程序）进行“打乱”。在蒙特卡洛采样中，低差异序列通常用于实现更有效的“随机”采样，但一般来说，只有在不跳过其任何点的情况下才能安全使用（Owen 2020）[^68]。
其他产生均匀分布的点样本的方法包括以下几种。
- 分层抽样将 N 维盒子分成大小相同的较小盒子，并在每个盒子中随机选择一个或多个点。
- 拉丁超立方抽样可以使用以下伪代码实现一个
n个数字的序列：lhs = []; for i in 0...n: AddItem(RNDRANGEMinMaxExc(i*1.0/n,(i+1)*1.0/n)); lhs = Shuffle(lhs)。 - 伪随机数生成器的特殊版本。一个例子是模
m、具有全周期和“良好晶格结构”的线性同余生成器；然后，对于区间 [1,m] 中的每个整数i，一个n维点序列是[MLCG(i), MLCG(i+1), ..., MLCG(i+n-1)]（L'Ecuyer 1999）[^69]。一个例子是MLCG(seed)：rem(92717*seed,262139)/262139.0。另一个例子是某些线性反馈移位寄存器生成器（Harase 2020）[^70]。 - 如果低差异序列输出区间 [0, 1] 中的数字，则序列的Baker 映射是
2 * (MakeRatio(1,2)-abs(x - MakeRatio(1,2)))，其中x是序列中的每个数字。 - Kamath (2022)[^71] 描述了其他随机点采样方法，包括泊松圆盘采样、“最佳候选算法”和 N-最远点算法。
关于实数随机化的说明
需要随机实数。
随机游走：其他示例
- 白噪声过程的一个例子是
Normal(0, 1)数字列表（高斯白噪声）。 - 如果
STATEJUMP()是RNDRANGEMinMaxExc(-1, 1)，则随机状态会推进一个 -1 到 1 之间的随机实数。 - 连续时间过程模拟的是在每个时刻的随机行为，而不仅仅是在离散时间。有两个流行的例子：
- 维纳过程（也称为布朗运动）具有随机状态和服从正态分布的跳跃。对于遵循维纳过程的随机游走，
STATEJUMP()是Normal(mu * timediff, sigma * sqrt(timediff))，其中mu是漂移（或每单位时间的平均值），sigma是波动率，timediff是样本之间的时间差。布朗桥（Revuz and Yor 1999）[^72] 修改维纳过程如下：对于每个时间 X，计算W(X) - W(E) * (X - S) / (E - S)，其中S和E分别是过程的开始和结束时间，W(X)和W(E)分别是时间 X 和 E 的状态。 - 在泊松点过程中，每次事件之间的时间是其自身的指数随机变量，具有其自身的速率参数（例如，
Expo(rate)）（参见“指数分布”），对 N 个随机RNDRANGEMinMaxExc(x, y)进行排序可以得到区间[x, y]中的 N 个到达时间。如果所有速率都相同，则该过程是齐次的；如果速率是每个事件跳跃之前“时间戳”（风险率函数）的函数，则该过程是非齐次的；要在这里生成到达时间，会以最大可能速率（maxrate）生成潜在到达时间，并在RNDRANGEMinMaxExc(0, maxrate) < thisrate时接受每个时间，其中thisrate是给定到达时间的速率（Lewis and Shedler 1979）[^73]。
- 维纳过程（也称为布朗运动）具有随机状态和服从正态分布的跳跃。对于遵循维纳过程的随机游走，
变换：其他示例
- Bates 分布：求给定范围内 n 个均匀随机变量的均值（例如，由
RNDRANGEMinMaxExc(minimum, maximum)生成的变量）（策略 8，均值）。 - 随机点 (
x,y) 可以通过变换（策略 9，几何变换）得到具有相关随机坐标（旧x，新x）的点，如下所示（参见（Saucier 2000）[^74]，第 3.8 节）：[x, y*sqrt(1 - rho * rho) + rho * x]，其中x和y是以相同方式随机选择的独立数字，rho是区间 [-1, 1] 中的相关系数（如果rho为 0，则x和y不相关）。 - 谈论 N 个随机变量的总和或平均值是合理的，其中 N 具有小数部分。在这种情况下，将生成
ceil(N)个随机变量，并将最后一个变量乘以该小数部分。例如，要采样 2.5 个随机变量的总和，请生成三个随机变量，将最后一个乘以 0.5（2.5 的小数部分），然后将这三个变量相加。 - 次指数分布模拟 n 个服从指数分布且每个都有单独速率参数的随机变量的总和（参见“指数分布”）。
- P. Jacob 提到的最大耦合方法[^XX] 生成来自两个分布 P 和 Q 的相关随机变量，这两个分布具有已知的概率密度函数或 PDF（分别为
PPDF和QPDF）；这仅在每个 PDF 下的面积为 1 时有效：从分布 P 中随机采样一个数字x，如果RNDRANGEMinMaxExc(0, PPDF(x)) < QPDF(x)，则返回[x, x]。否则，从分布 Q 中随机采样一个数字y，直到PPDF(y) < RNDRANGEMinMaxExc(0, QPDF(y))，然后返回[x, y]。
从数据点分布中采样
需要随机实数。
基于数据点列表的分布生成随机数据点涉及机器学习领域：拟合数据模型到数据点，然后基于该模型预测新的数据点，并加入随机性。下面将描述的三种数据模型可用于此目的。（拟合工作原理超出本页范围。）
-
密度估计模型。密度估计模型旨在描述给定数据集中的数据点的分布，其中点数较多的区域有更大的几率被采样。[^75] 以下是示例。
- 直方图是一组一个或多个非重叠的箱，通常大小相等。直方图是混合，每个箱的权重是该箱中的数据点数。随机选择一个箱后，生成一个可能适合该箱的随机数据点（该点不必是现有数据点）。
- 高斯混合模型也是混合，在这种情况下，是多个高斯（正态）分布的混合。
- 核密度分布是采样分布的混合，每个数据点一个。估计核密度分布称为核密度估计。要从核密度分布采样
- 从列表中选择一个数字或点，并进行均匀随机有放回抽样。
- 为选定的数字或点添加随机化的“抖动”；例如，将单独生成的
Normal(0, sigma)添加到选定的数字或选定点的每个分量，其中sigma是带宽[^76]。
- 随机插值在（Saucier 2000）[^74]，第 5.3.4 节中有描述。
- 拟合已知分布（例如正态分布），具有未知参数，可以使用最大似然估计等方法。
-
回归模型。回归模型是一种将数据概括为公式和误差项的模型。如果应用程序具有输入和输出形式的数据（例如，月销售额），并且想在给定已知输入点的情况下采样一个随机但合理的输出（例如，未来一个月的销售额），则应用程序可以拟合并采样该数据的回归模型。例如，线性回归模型模拟给定已知输入
a和b时的y值，可以如下采样：y = c1 * a + c2 * b + c3 + Normal(0, sqrt(mse))，其中mse是均方误差，c1、c2和c3是模型的系数。（这里，Normal(0, sqrt(mse))是误差项。） -
生成模型。这些是机器学习模型，它们以随机变量作为输入，并生成与它们已经见过的数据相似的输出（例如图像或声音）。
注释
- 通常，可以选择一种或多种数据模型和/或机器学习模型来拟合给定数据集（例如，多种密度估计模型、回归模型、参数分布和/或决策树）。如果存在多种模型作为拟合选择，则应选择显示对数据集可接受的预测准确性（例如，信息准则、精确度、召回率）的最简单的模型。
- 如果现有的数据点每个都属于几个类别中的一个，则可以通过选择一个随机数，其概率与每个类别中的数据点数量成正比来选择一个随机类别（参见“加权选择”）。
- 如果现有的数据点每个都属于几个类别中的一个，则可以通过——
- 选择一个随机数据点，然后找到其类别（例如，通过称为分类模型的机器学习模型），或者
- 选择一个如上给出的随机类别，然后仅根据该类别的现有数据点选择一个随机数据点。
从任意分布中采样
需要随机实数。
许多概率分布可以根据以下任一方式定义：
- 累积分布函数，或 CDF，
CDF(x)，是随机选择小于或等于x的数字的概率。概率在 [0, 1] 区间内。 - 离散分布[^77] 具有概率质量函数，或 PMF，它给出每个数字被随机选择的概率。
- 绝对连续分布具有概率密度函数，或 PDF，
PDF(x)，它是 CDF 的“斜率”函数，或随机选择“接近”x的数字的相对概率。相对概率为 0 或更大，PDF 下方的面积为 1。 - 分位数函数（也称为反累积分布函数或反 CDF）将区间 (0, 1) 中的数字映射到分布中的数字，从低到高。
在本节中，类 PDF 函数是 PDF、PMF 或任何函数乘以（可能未知的）正常数。
以下各节展示了根据分布生成随机变量的不同方法，具体取决于对该分布的了解程度。
注意：CDF、类 PDF 函数或分位数函数的列表超出了本文档的范围。
离散分布采样
如果分布是离散的，则可以通过选择覆盖分布的全部或几乎全部区域的点，找到它们的权重或累积权重，并根据这些权重选择一个随机点来采样接近该分布的数字。
如果——
- 离散分布具有已知的类 PDF 函数
PDF(x)，其中x必须是整数， - 区间 [
mini,maxi] 覆盖了该分布的全部，并且 - 函数的取值都是有理数（形式为
y/z的数字，其中y和z是整数），
则以下方法可以从该分布中精确采样。
METHOD SampleDiscrete(mini, maxi)
// Setup
ratios=[]
for i in mini..maxi: AddItem(ratios, PDF(i))
ratios=NormalizeRatios(ratios)
// Sampling
return mini + WeightedChoice(ratios)
END METHOD
如果——
- 离散分布具有已知的 CDF
CDF(x)，其中x必须是整数， - 区间 [
mini,maxi] 覆盖了该分布的全部，并且 - CDF 的取值都是有理数，
则以下方法可以从该分布中精确采样。
METHOD SampleDiscreteCDF(mini, maxi)
// Setup
ratios=[MakeRatio(0,1)]
for i in mini..maxi: AddItem(ratios, CDF(i))
ratios=NormalizeRatios(ratios)
// Sampling
value=ratios[size(ratios) - 1]
for i in 0...size(ratios) - 1
if ratios[i] < ratios[i+1] and
ratios[i]>=value: return mini + i
end
return mini
END METHOD
在其他情况下，离散分布仍然可以近似采样。以下情况将导致近似采样器，除非 CDF 或类 PDF 函数的值覆盖了全部分布并且被精确计算（无误差）。
- 如果 CDF 或类 PDF 函数的值作为形式为
FPSignificand*FPRadixFPExponent的浮点数计算（包括 Java 的double和float）[^78]，则有多种方法可以将这些数字转换为有理数或整数。- 一种方法是使用
FPRatio(x)（在下面的伪代码中），这是无损的，并且可以计算给定浮点数x的有理数。 - 另一种方法是将值缩放并四舍五入为整数（例如，
floor(x * mult)如果floor(x * mult) < 0.5否则ceil(x * mult)，其中mult是一个大的整数）；这是有损的。 - 第三种方法是以限制误差的方式将类 PDF 函数的值近似为整数，例如（Saad et al., 2020）[^79] 中给出的方法；这是有损的，并且仅适用于类 PDF 函数。
- 一种方法是使用
- 如果 CDF 或类 PDF 函数的值作为有理数计算，则可以使用
NormalizeRatios（无损）或上述 (2) 或 (3)（有损）将这些数字转换为整数权重。 - 如果分布取无限多个值，则可以选择一个合适的区间 [
mini,maxi] 来覆盖该分布的几乎全部。
METHOD FPRatio(fp)
expo=FPExponent(fp)
sig=FPSignificand(fp)
radix=FPRadix(fp)
if expo>=0: return MakeRatio(sig * pow(radix, expo), 1)
return MakeRatio(sig, pow(radix, abs(expo)))
END METHOD
注意：此外，一些分布仅通过一个生成该分布随机变量的神谕（或“黑箱”）来得知。算法可以使用此神谕来生成遵循不同分布的新随机变量。当新变量的均值（“长期平均值”）等于神谕的均值时，该算法称为无偏估计量（因为从某种意义上说，它“估计”了神谕数字的某个方面）。一个例子是伯努利工厂（请参阅我的文章“伯努利工厂算法”），它接收具有不同正面概率（神谕）的硬币翻转，并生成具有不同正面概率的新“硬币”的翻转。另一个例子是加权选择中描述的“伯努利比赛”。
逆变换采样
逆变换采样（或简单地说反演）是从概率分布中采样数字的最通用方法。
如果分布具有已知的分位数函数，则生成 (0, 1) 中的一个均匀随机变量（如果该数字尚未预先生成），并取该数字的分位数。然而
- 在大多数情况下，无法获得分位数函数。因此，它必须被近似。
- 即使分位数函数可用，朴素的分位数计算（例如，
ICDF(RNDRANGEMinMaxExc(0, 1))）也可能意味着均匀数字的微小变化会导致分位数发生巨大变化，从而在采样覆盖中产生差距（Monahan 1985，第 4 和 6 节）[^65]。
以下方法通过反演从分布中采样，精度为 1/BASEprecision（（Devroye and Gravel 2020）[^11]，但扩展到任何基数；另请参阅（Bringmann and Friedrich 2013，附录 A）[^53]）。在该方法中，ICDF(u, ubits, prec) 返回一个两项列表，分别包含一个数字的上限和下限，该数字在 true 分位数 u/BASEubits 的 1/BASEprec 范围内，而 BASE 是数字基数（例如，二进制为 2，十进制为 10）。
METHOD Inversion(precision)
u=0
ubits=0
threshold=MakeRatio(1,pow(BASE, precision))*2
incr=8
while true // until a value is returned
incr=8
if ubits==0: incr=precision
// NOTE: If a uniform number (`n`) is already pregenerated,
// use the following instead:
// u = rem(floor(n*pow(BASE, ubits+incr)), pow(BASE, incr))
u=u*pow(BASE,incr)+RNDINTEXC(pow(BASE,incr))
ubits=ubits+incr
// Get upper and lower bound
bounds=ICDF(u,ubits,precision)
if lower>upper: return error
diff=bounds[1]-bounds[0]
if diff<=threshold: return bounds[1]+diff/2
end
end
某些应用程序需要通过分位数将 [0, 1] 中的预生成数字（通常是从均匀分布中采样的数字）u01 转换为非均匀分布。显著的例子包括 copula 方法、阶统计量和涉及低差异序列的蒙特卡洛方法。以下计算分位数的方法在理论上是精确的
- 分布是离散的，具有已知的 PMF（并且分布取整数）：顺序搜索（Devroye 1986，第 85 页）[^19]：
i = 0; p = PMF(i); while u01 > p; u01 = u01 - p; i = i + 1; p = PMF(i); end; return p，但这并不总是很快。（这仅在PMF的值总和为 1 时有效，这就是为什么此处只允许使用 PMF 而不是类 PDF 函数。）
此外，以下方法近似分位数
- 分布是离散的，具有已知的类 PDF 函数（并且分布取整数）：如果区间 [a, b] 覆盖了该分布的全部或几乎全部，则应用程序可以将该区间内的类 PDF 函数的值存储在一个列表中，然后调用
WChoose：wsum = 0; for i in a..b: wsum=wsum+PDF(i); for i in a..b: AddItem(weights, PDF(i)); return a + WChoose(weights, u01 * wsum)。但是，基于CDF而非类 PDF 函数查找分位数可能会引入更多误差（Walter 2019）[^80]。另请参阅Python 示例代码中的integers_from_u01。 - 分布是绝对连续的，具有已知的类 PDF 函数：
ICDFFromContPDF(u01, mini, maxi, step)，下面找到基于类 PDF 函数在 [mini,maxi] 中的分段线性逼近的近似分位数，其中分段宽度最多为step。（Devroye and Gravel 2020）[^11]。另请参阅Python 示例代码中的DensityInversionSampler、numbers_from_u01和numbers_from_dist_inversion（Derflinger et al. 2010）[^81]、（Devroye and Gravel 2020）[^11]。 - 分布是绝对连续的，具有已知的 CDF：请参阅Python 示例代码中的
numbers_from_u01。
METHOD ICDFFromContPDF(u01, mini, maxi, step)
pieces=[]
areas=[]
// Setup
lastvalue=i
lastweight=PDF(i)
cumuarea=0
i = mini+step; while i <= maxi
weight=i; value=PDF(i)
cumuarea=cumuarea+abs((weight + lastweight) * 0.5 *
(value - lastvalue))
AddItem(pieces,[lastweight,weight,lastvalue,value])
AddItem(areas,cumuarea)
lastweight=weight;lastvalue=value
if i==maxi: break
i = min(i + step, maxi)
end
for i in 0...size(areas): areas[i]=areas[i]/cumuarea
// Sampling
prevarea=0
for i in 0...size(areas)
cu=areas[i]
if u01<=cu
p=pieces[i]; u01=(u01-prevarea)/(cu-prevarea)
s=p[0]; t=p[1]; v=u01
if s!=t: v=(s-sqrt(t*t*u01-s*s*u01+s*s))/(s-t)
return p[2]+(p[3]-p[2])*v
end
prevarea=cu
end
return error
END METHOD
注释
- 如果仅有数据的百分位数（例如中位数或第 50 百分位数、最小值或第 0 百分位数、最大值或第 100 百分位数）可用，则可以通过这些百分位数近似分位数函数。第 N 百分位数对应于
N/100.0的分位数。然后可以通过插值（例如样条拟合）来填充分位数函数的缺失值。如果原始数据点可用，请参阅“从数据点分布中采样”。- 取
n个相同分布的随机变量的第k小值，等同于取n个 [0, 1) 区间内的均匀随机变量的第k小值（也称为第k个阶统计量；例如，BetaDist(k, n+1-k)）并找到其分位数（Devroye 2006）[^83]；（Devroye 1986，第 30 页）[^19]。
带 PDF 状函数的拒绝采样
如果分布具有已知的类 PDF 函数 (PDF)，并且该函数可以更容易地通过另一个具有其自身类 PDF 函数 (PDF2) 的分布进行采样，该分布“支配”PDF，即在每个有效 x 处 PDF2(x) >= PDF(x)，则生成后者的分布的随机变量，直到生成满足 r < PDF(n) 的变量（称为 n），其中 r = RNDRANGEMinMaxExc(0, PDF2(n))（即，在 PDF2 中采样点，直到一个点落在 PDF 范围内）。
拒绝采样的变体是挤压原理，其中选择第三个类 PDF 函数 (PDF3)，它被第一个函数 (PDF)“支配”且比 PDF 更容易采样。这里，当 r < PDF3(n) 或 r < PDF(n) 时接受一个数字，其中 r = RNDRANGEMinMaxExc(0, PDF2(n))（Devroye 1986，第 53 页）[^19]。
另请参阅（von Neumann 1951）[^44]；（Devroye 1986）[^19]，第 41-43 页；“拒绝采样”；以及“生成伪随机数”。
示例
- 要从 PDF 峰值不大于
peak的区间 [low,high) 中的一个 PDF 类函数采样随机变量，请生成x = RNDRANGEMinMaxExc(low, high)和y = RNDRANGEMinMaxExc(0, peak)直到y < PDF(x)，然后取最后生成的x。（另请参阅 Saucier 2000，第 6-7 页。）如果分布是离散的且为整数值，则生成x，其中x = RNDINTEXCRANGE(low, high)。- 自定义分布的类 PDF 函数
exp(-abs(x*x*x))，而指数分布的PDF2是exp(-x)。指数类 PDF 函数PDF2“支配”x大于等于 0 的地方），如果我们将其乘以 1.5，则PDF2变为1.5 * exp(-x)。现在我们可以通过对指数点进行采样直到一个点落在n = Expo(1)直到RNDRANGEMinMaxExc(0, PDF2(n)) < PDF(n)来完成的。- 正态分布的倒置钟形曲线的类 PDF 函数为
1-exp(-(x*x))，该函数在该区间内的最高点为peak = max(1-exp(-(low*low)), 1-exp(-(high*high)))。采样此分布然后使用示例 1 中的算法。注意：在 Python 示例代码中，moore.py 和
numbers_from_dist通过拒绝采样从分布中采样（Devroye and Gravel 2020）[^11]、（Sainudiin and York 2013）[^84]。
交错级数
如果目标分布的类 PDF 函数不精确知道，但可以用两个级数展开从上方和下方近似，这两个级数随着更多项的加入会收敛到该函数，则可以使用交替级数法。这仍然需要一个“支配”的类 PDF 函数 (PDF2(x)) 作为“易于采样”的分布。将级数展开称为 UPDF(x, n) 和 LPDF(x, n)，其中 n 是级数要相加的项数。要使用此方法（Devroye 2006）[^83] 对该分布进行采样：（1）从“支配”分布中采样，并将采样到的数字设为 x；（2）将 n 设置为 0；（3）如果 r < LPDF(x, n) 则接受 x，或者如果 r >= UPDF(x, n) 则转到步骤 1，或者如果都不是，则通过将 n 加 1 来重复此步骤，其中 r = RNDRANGEMinMaxExc(0, PDF2(n))。
马尔可夫链蒙特卡罗
马尔可夫链蒙特卡洛（MCMC）是一系列算法，用于通过构建随机值马尔可夫链来从概率分布中采样，该链随着链的步数增加而接近给定分布。但总的来说，MCMC 是近似的，并且由给定链生成的值彼此之间存在统计依赖性（这就是为什么通常采用“抽样”（仅保留每 N 个样本）等技术的原因）。[^85]
MCMC 算法[^86] 包括Metropolis–Hastings、slice sampling 和Gibbs sampling（另请参阅Python 示例代码）。后者比较特别，因为它不使用类 PDF 函数，而是使用两个或多个分布，每个分布使用从前一个分布（条件分布）中随机采样的数字，这些分布收敛到联合分布。
示例：在一个 Gibbs 采样器中，先选择一个初始
y值，然后生成多个x,y对的随机变量，其中x = BetaDist(y, 5)然后y = Poisson(x * 10)。
分段线性分布
需要随机实数。
分段线性分布描述了一种绝对连续分布，该分布具有已知点的权重，其他权重由线性插值（平滑）确定。PiecewiseLinear 方法（在下面的伪代码中）接受两个列表，如下所示（另请参阅（Kscischang 2019）[^87]）：
values是一个有理数列表。数字应按升序排列。weights是给定数字的有理值权重列表（其中每个数字及其权重在两个列表中具有相同的索引）。数字的权重越大，选择接近该数字的数字的概率就越大。每个权重应大于等于 0。
METHOD PiecewiseLinear(values, weights)
if size(values)!=size(weights) or size(values)==0: return error
if size(values)==1: return values[0]
areas=[]
for i in 1...size(values)
area=abs((weights[i] + weights[i-1]) *
(values[i] - values[i-1]) / 2) // NOTE: Not rounded
AddItem(areas,area)
end
// NOTE: If values and weights are rational
// numbers, use `areas=NormalizeRatios(areas)` instead
// of finding `areas` as given below.
ratios=[]
for w in areas: AddItem(ratios, FPRatio(w))
areas=NormalizeRatios(ratios)
index=WeightedChoice(areas)
w=values[index+1]-values[index]
if w==0: return values[index]
m=(weights[index+1]-weights[index])/w
h2=(weights[index+1]+weights[index])
ww=w/2.0; hh=h2/2.0
x=RNDRANGEMinMaxExc(-ww, ww)
if RNDRANGEMinMaxExc(-hh, hh)>x*m: x=-x
return values[index]+x+ww
END METHOD
注意：Python 示例代码包含该方法的变体，用于一次返回多个随机变量。
示例：假设
values如下：[0, 1, 2, 2.5, 3]，weights如下：[0.2, 0.8, 0.5, 0.3, 0.1]。2 的权重是 0.5，2.5 的权重是 0.3。由于 2 的权重大于 2.5，使用PiecewiseLinear方法，接近 2 的数字被选择的概率大于接近 2.5 的数字。
特定分布
在单独页面中提供了采样其他分布的方法。它们涵盖了正态、伽马、beta、von Mises、稳定和多元正态分布以及 copulas。
非均匀分布索引
许多分布需要随机实数。
分布旁边的 † 符号表示该分布的样本可以通过位置参数（mu）进行平移，然后通过大于 0 的尺度参数（sigma）进行缩放。示例：num * sigma + mu。
分布旁边的 ⬦ 符号表示样本可以缩放到任何范围，该范围由最小值 mini 和最大值 maxi 给出。示例：mini + (maxi - mini) * num。
有关进一步的示例和分布，请参阅（Devroye 1996）[^88] 和（Crooks 2019）[^89]。
最常用的
- Beta 分布⬦：请参阅Beta 分布。
- 二项分布：请参阅二项分布。
- 双正态分布：请参阅多元正态（多正态）分布。
- Cauchy（洛伦兹）分布†：
Stable(1, 0)。此分布类似于正态分布，但具有“更肥厚的”尾部。基于（McGrath and Irving 1975）[^90] 中提到的算法的替代算法：生成x = RNDRANGEMinMaxExc(0,1)和y = RNDRANGEMinMaxExc(0,1)直到x * x + y * y <= 1，然后生成(RNDINT(1) * 2 - 1) * y / x。 - 卡方分布：
GammaDist(df * 0.5 + Poisson(sms * 0.5))*2，其中df是自由度，sms是均方和（其中sms非 0 表示非中心分布）。 - 骰子：请参阅骰子。
- 指数分布：请参阅指数分布。朴素的实现
-ln(1-RNDRANGEMinMaxExc(0, 1)) / lamda有几个问题，例如在大值处因分布的右侧尾部而条件病态（Pedersen 2018）[^1]。应用程序可以通过应用 Pedersen 的建议来减少其中一些问题，即使用-ln(RNDRANGEMinMaxExc(0, 0.5))或-log1p(-RNDRANGEMinMaxExc(0, 0.5))（而不是-ln(1-RNDRANGEMinMaxExc(0, 1))），每次随机均匀选择；另一种选择是（Devroye 2006）[^83] 中提到的ln(1/RNDRANGEMinMaxExc(0,1))。 - 极值分布：请参阅广义极值分布。
- Gamma 分布：请参阅Gamma 分布。广义伽马分布包括Stacy 分布（
pow(GammaDist(a), 1.0 / c) * b，其中c是另一个形状参数）和Amoroso 分布（Crooks 2015）[^91]、（pow(GammaDist(a), 1.0 / c) * b + d，其中d是最小值）。 - 高斯分布：请参阅正态（高斯）分布。
- 几何分布：请参阅几何分布。假设计算机可以“精确”地处理实数，以下是“精确”的：
floor(-Expo(1)/ln(1-p))（Devroye 1986，第 500 页）[^19]（将 ceil 替换为 floor，因为本页对几何分布的定义不同）。 - Gumbel 分布：请参阅广义极值分布。
- 逆 Gamma 分布：
b / GammaDist(a)，其中a和b的含义与 Gamma 分布相同。或者，1.0 / (pow(GammaDist(a), 1.0 / c) / b + d)，其中c和d分别是形状和位置参数。 - Laplace（双指数）分布†：
(Expo(1) - Expo(1))。另外，Normal(0,1) * Normal(0, 1) - Normal(0, 1) * Normal(0, 1)（Kotz et al. 2012）[^92]。 - 对数分布⬦：
RNDRANGEMinMaxExc(0, 1) * RNDRANGEMinMaxExc(0, 1)（Saucier 2000，第 26 页）。在此分布中，较低的数字比较高的数字指数级更可能出现。 - 对数正态分布：
exp(Normal(mu, sigma))，其中mu和sigma是底层正态分布的参数。 - 多正态分布：请参阅多元正态分布。
- 多元正态分布：请参阅多元正态（多正态）分布。
- 正态分布：请参阅正态（高斯）分布。
- 泊松分布：请参阅“泊松分布”。假设计算机可以“精确”地处理实数，以下是“精确”的（Devroye 1986，第 504 页）[^19]：
c = 0; s = 0; while true; sum = sum + Expo(1); if sum>=mean: return c; else: c = c + 1; end；此外，还可以使用（Devroye 1991）[^93] 中的以下优化：while mean > 20; n=ceil(mean-pow(mean,0.7)); g=GammaDist(n); if g>=mean: return c+(n-1-Binomial(n-1,(g-mean)/g)); mean = mean - g; c = c + n; end。对于mean大于 50，以下近似由（Giammatteo and Di Mascio 2020）[^92] 建议：floor(1.0/3 + pow(max(0, Normal(0, 1)*pow(mean,1/6.0)*2/3 + pow(mean, 2.0/3)), 3.0/2))。 - Pareto 分布：
pow(RNDRANGEMinMaxExc(0, 1), -1.0 / alpha) * minimum，其中alpha是形状参数，minimum是最小值。 - Rayleigh 分布†：
sqrt(Expo(0.5))。如果尺度参数（sigma）服从对数正态分布，则结果是Suzuki 分布。 - 标准正态分布†：
Normal(0, 1)。另请参阅正态（高斯）分布。 - Student's t 分布：
Normal(cent, 1) / sqrt(GammaDist(df * 0.5)*2 / df)，其中df是自由度，cent 是正态分布随机变量的均值。cent非 0 表示非中心分布。或者，cos(RNDRANGEMinMaxExc(0, pi * 2)) * sqrt((pow(RNDRANGEMinMaxExc(0, 1),-2.0/df)-1) * df)（Bailey 1994）[^94]。 - 三角分布†（Stein and Keblis 2009）[^95]：
(1-alpha) * min(a, b) + alpha * max(a, b)，其中alpha在 [0, 1] 中，a = RNDRANGEMinMaxExc(0, 1)，b = RNDRANGEMinMaxExc(0, 1)。 - Weibull 分布：请参阅广义极值分布。
杂项
- Archimedean copulas：请参阅高斯和其他 copulas。
- Arcsine 分布⬦：
BetaDist(0.5, 0.5)（Saucier 2000，第 14 页）。 - Bates 分布：请参阅随机变量变换：附加示例。
- Beckmann 分布：请参阅多元正态（多正态）分布。
- Beta 二项分布：
Binomial(trials, BetaDist(a, b))，其中a和b是 beta 分布的两个参数，trials是二项分布的参数。 - Beta 负二项分布：
NegativeBinomial(successes, BetaDist(a, b))，其中a和b是 beta 分布的两个参数，successes是负二项分布的参数。如果 successes 为 1，结果是Waring–Yule 分布。如果 successes 和 b 都为 1（例如，在 Mathematica 中），则结果为Yule–Simon 分布；如果 successes 和 a 都为 1（在其他文献中）。 - Beta-PERT 分布：
startpt + size * BetaDist(1.0 + (midpt - startpt) * shape / size, 1.0 + (endpt - midpt) * shape / size)。该分布从startpt开始，在midpt处达到峰值，并在endpt结束，size是endpt - startpt，shape是一个大于等于 0 的形状参数，但通常为 4。如果已知均值（mean）而不是峰值，则midpt = 3 * mean / 2 - (startpt + endpt) / 4。 - Beta 质数分布†：
pow(GammaDist(a), 1.0 / alpha) / pow(GammaDist(b), 1.0 / alpha)，其中a、b和alpha是形状参数。如果 a 为 1，则结果是Singh–Maddala 分布；如果 b 为 1，则是Dagum 分布；如果 a 和 b 都为 1，则是对数logistic 分布。 - Birnbaum–Saunders 分布：
pow(sqrt(4+x*x)+x,2)/(4.0*lamda)，其中x = Normal(0,gamma)，gamma是形状参数，lamda是尺度参数。 - Borel 分布（Borel 1942）[^108]：
r=0; q=1; while q>=1; q+=Poisson(la); q-=1; r+=1; end; return r。la，平均到达数，应在区间 (0, 1) 内。 - Chi 分布：卡方随机变量的平方根。请参阅卡方分布。
- 复合泊松分布：请参阅随机变量变换：附加示例。
- 余弦分布⬦：
atan2(x, sqrt(1 - x * x)) / pi，其中x = (RNDINT(1) * 2 - 1) * RNDRANGEMinMaxExc(0, 1)（Saucier 2000，第 17 页；将反正弦替换为atan2等效项）。 - CUB 分布（Piccolo 2003）[^107]：
if ZeroOrOne(px,py)==1: return 1+BinomialInt(m-1, zy-zx, zy); else: return RNDINTRANGE(1, m)，其中m>3，px/py在 [0, 1] 内，zx/zy在 [0, 1] 内。 - Dagum 分布：请参阅 beta 质数分布。
- Dirichlet 分布：假设我们（1）生成 n+1 个随机Gamma 分布的随机变量，每个都有单独的参数；（2）求它们的和；（3）将每个数除以该总和；然后（4）将每个数乘以大于 0 的实数
x。那么- 在步骤 (4) 之后，如果
x为 1，则Dirichlet 分布（例如，（Devroye 1986）[^23]，第 593-594 页）模拟这些数字的前 n 个。 - 如果步骤 (1) 中的数字每个都生成为
Expo(1)（Gamma 分布的一个特例），则步骤 (4) 之后的结果是总和为x的 n+1 个数字的均匀分布的和（另请参阅上面的链接文章）。
- 在步骤 (4) 之后，如果
- 双对数分布⬦：
(0.5 + (RNDINT(1) * 2 - 1) * RNDRANGEMinMaxExc(0, 0.5) * RNDRANGEMinMaxExc(0, 1))（另请参阅 Saucier 2000，第 15 页，其中显示了错误的 X 轴）。 - Erlang 分布：
GammaDist(n)/lamda，其中n是大于 0 的整数。返回一个模拟具有给定lamda参数的n个指数随机变量总和的数字。 - Estoup 分布：请参阅 zeta 分布。
- 指数幂分布（广义正态分布版本 1）：
(RNDINT(1) * 2 - 1) * pow(GammaDist(1.0/a), a)，其中a是形状参数。 - Fréchet 分布：请参阅广义极值分布。
- Fréchet–Hoeffding 下界 copula：请参阅高斯和其他 copulas。
- Fréchet–Hoeffding 上界 copula：请参阅高斯和其他 copulas。
- 高斯 copula：请参阅高斯和其他 copulas。
- 广义极值（Fisher–Tippett 或广义最大值）分布 (
GEV(c))†：如果c != 0则为(pow(Expo(1), -c) - 1) / c，否则为-ln(Expo(1))，其中c是形状参数。特殊情况- 结果的负值表示广义最小值。在这种情况下，参数
c = 0产生Gumbel 分布。 - 参数
c = 0产生极值分布。 - Weibull 分布：
1 - 1.0/a * GEV(-1.0/a)（或pow(Expo(1), 1.0/a)），其中a是形状参数。 - Fréchet 分布：
1 + 1.0/a * GEV(1.0/a)（或pow(Expo(1), -1.0/a)），其中a是形状参数。
- 结果的负值表示广义最小值。在这种情况下，参数
- 广义 Tukey lambda 分布：
(s1 * (pow(x, lamda1)-1.0)/lamda1 - s2 * (pow(1.0-x, lamda2)-1.0)/lamda2) + loc，其中x是RNDRANGEMinMaxExc(0, 1)，lamda1和lamda2是形状参数，s1和s2是尺度参数，loc是位置参数。 - 半正态分布。参数化包括
- Mathematica：
abs(Normal(0, sqrt(pi * 0.5) / invscale)))，其中invscale是半正态分布的参数。 - MATLAB：
abs(Normal(mu, sigma)))，其中mu和sigma是底层正态分布的参数。
- Mathematica：
- 超指数分布：请参阅分布混合。
- 超几何分布：请参阅Polya–Eggenberger 分布。
- 次指数分布：请参阅随机变量变换。
- 逆卡方分布†：
df / (GammaDist(df * 0.5)*2)，其中df是自由度。尺度参数（sigma）通常为1.0 / df。 - 逆高斯分布（Wald 分布）：生成
n = mu + (mu*mu*y/(2*lamda)) - mu * sqrt(4 * mu * lamda * y + mu * mu * y * y) / (2 * lamda)，其中y = pow(Normal(0, 1), 2)，然后以mu / (mu + n)的概率返回n（例如，如果RNDRANGEMinMaxExc(0, 1) <= mu / (mu + n)），否则返回mu * mu / n。mu是均值，lamda是尺度；两个参数都大于 0。基于 (Devroye 1986)[^19] 发表的方法。 k阶统计量：BetaDist(k, n+1-k)。返回n个 [**0, 1) 区间内的均匀随机变量中的第k小值。另见 (Devroye 1986, p. 210)[^19]。- Kumaraswamy 分布⬦：
pow(BetaDist(1, b), 1.0 / a)，其中a和b是形状参数。 - Landau 分布：见稳定分布。
- Lévy 分布†：
0.5 / GammaDist(0.5)。尺度参数（sigma）也称为离散度。 - 对数Logistic 分布：见 Beta Prime 分布。
- 对数级数分布：生成
n = NegativeBinomialInt(1, py - px, py)+1（其中px/py是 (0,1) 区间内的参数），然后如果ZeroOrOne(1, n) == 1则返回n，否则重复此过程（Flajolet 等人，2010）[^96]。以下是“精确”的，假设计算机可以“精确”地处理实数：floor(1.0 - Expo(log1p(-pow(1.0 - p, RNDRANGEMinMaxExc(0, 1)))))，其中p是 (0, 1) 区间内的参数；见 (Devroye 1986)[^19]。 - Logistic 分布†：
(ln(x)-log1p(-x))(logit 函数)，其中x是RNDRANGEMinMaxExc(0, 1)。 - Log-multinormal 分布：见 多变量正态（Multinormal）分布。
- 最大均匀分布：
BetaDist(n, 1)。返回一个模拟 [**0, 1) 区间内n个均匀随机变量中最大值的数字。另见 (Devroye 1986, p. 675)[^19]。 - Maxwell 分布†：
sqrt(GammaDist(1.5)*2)。 - 最小均匀分布：
BetaDist(1, n)。返回一个模拟 [**0, 1) 区间内n个均匀随机变量中最小值的数字。另见 (Devroye 1986, p. 210)[^19]。 - Moyal 分布：见 Python 示例代码。
- Multinomial 分布：见 Multinomial 分布。
- Multivariate Poisson 分布：见 Python 示例代码。
- Multivariate t-copula：见 Python 示例代码。
- Multivariate t-distribution：见 Python 示例代码。
- 负二项分布（
NegativeBinomial(successes, p)）：见 负二项分布。以下是“精确”的，假设计算机可以“精确”地处理实数：Poisson(GammaDist(successes)*(1 - p) / p)（即使successes不是整数也能工作）。 - 负多项分布：见 Python 示例代码。
- 非中心 Beta 分布⬦：
BetaDist(a + Poisson(nc), b)，其中nc（非中心性）、a和b都大于 0。 - 抛物线分布⬦：
BetaDist(2, 2)(Saucier 2000, p. 30)。 - Pascal 分布：
NegativeBinomial(successes, p) + successes，其中successes和p的含义与负二项分布相同，但successes始终是整数。 - Pearson VI 分布：
GammaDist(v) / GammaDist(w)，其中v和w是大于 0 的形状参数 (Saucier 2000, p. 33；其中定义了一个额外的b参数，但在源代码中该参数已被抵消)。 - 分段常数分布：见 有放回加权选择。
- 分段线性分布：见 分段线性分布。
- Pólya–Aeppli 分布：见 随机变量的变换：附加示例。
- 幂分布：
BetaDist(alpha, 1) / b，其中alpha是形状，b是域。名义上在区间 (0, 1) 内。 - 幂律分布：
pow(RNDRANGEMinMaxExc(pow(mn,n+1),pow(mx,n+1)), 1.0 / (n+1))，其中n是指数，mn是最小值，mx是最大值。参考。 - 幂对数正态分布：见 Python 示例代码。
- 幂正态分布：见 Python 示例代码。
- 积分布（Product copula）：见 高斯和其他 Copula。
- Rice 分布：见 多变量正态（Multinormal）分布。
- Rice–Norton 分布：见 多变量正态（Multinormal）分布。
- Singh–Maddala 分布：见 Beta Prime 分布。
- sin^k 分布：生成
x = BetaDist(k+1, k+1) * pi，然后如果0-Expo(k) <= ln(pi*pi*sin(x) / ((4*x*(pi - x)))则返回x，否则重复此过程（Makalic 和 Schmidt 2018）[^97]。 - Skellam 分布：
Poisson(mean1) - Poisson(mean2)，其中mean1和mean2是Poisson方法中使用的均值。 - 偏态正态分布† (Ghorbanzadeh 等人 2014)[^98]：生成
c*max(a, b) + (1-c)*min(a, b)，其中 a =Normal(0, 1)且独立地，b = Normal(0, 1)，以及c = (1+th)/sqrt(2.0*(1+th))，且th是 [0, 1] 区间内的实数。特殊情况：如果th=0，生成Normal(0, 1)；如果th=1，生成max(a, b)；如果th=1，生成min(a, b)。 - Snedecor (Fisher) F 分布：
GammaDist(m * 0.5)*n / (GammaDist(n * 0.5 + Poisson(sms * 0.5)) * m)，其中m和n是两个卡方分布随机变量的自由度，如果sms非零，则其中一个分布是*非中心*的，均方和等于sms。 - 稳定分布：见 稳定分布。四参数稳定分布：
Stable(alpha, beta) * sigma + mu，其中mu是均值，sigma是尺度；如果alpha和beta为 1，则结果是*Landau 分布*。“Type 0”稳定分布：Stable(alpha, beta) * sigma + (mu - sigma * beta * x)，其中x是ln(sigma)*2.0/pi（如果alpha为 1），否则是tan(pi*0.5*alpha)。 - 标准复正态分布：见 多变量正态（Multinormal）分布。
- Suzuki 分布：见 Rayleigh 分布。
- Tukey lambda 分布：
(pow(x, lamda)-pow(1.0-x,lamda))/lamda，其中x是RNDRANGEMinMaxExc(0, 1)且lamda是形状参数。 - Twin-t 分布 (Baker 和 Jackson 2018)[^99]：生成
x，一个随机的学生 t 分布数（非中心）。以z = pow((1 + y) / ((1 + y * y) + y), (df + 1) * 0.5)的概率接受x（例如，如果RNDRANGEMinMaxExc(0, 1) < z），其中y = x * x / df且df是生成该数时使用的自由度；否则重复此过程。 - von Mises 分布：见 von Mises 分布。
- Waring–Yule 分布：见 Beta Negative Binomial 分布。
- Wigner（半圆）分布†：
(BetaDist(1.5, 1.5)*2-1)。尺度参数（sigma）是半圆半径。 - Yule–Simon 分布：见 Beta Negative Binomial 分布。
- Zeta 分布：生成
n = floor(pow(RNDRANGEMinMaxExc(0, 1), -1.0 / r))，如果d / pow(2, r) < RNDRANGEMinMaxExc((d - 1) * n / (pow(2, r) - 1.0))，其中d = pow((1.0 / n) + 1, r)，则重复此过程。参数r大于 0。基于 (Devroye 1986)[^19] 中描述的方法。通过拒绝大于某个大于 0 的整数的随机值来*截断*的 zeta 分布称为*Zipf 分布*或*Estoup 分布*。（Devroye 使用“Zipf 分布”来指代未截断的 zeta 分布。） - Zipf 分布：见 Zeta 分布。
几何采样
需要随机实数。
本节介绍在几何形状内部或表面选择独立均匀随机点的方法。
单纯形内的随机点
以下伪代码生成一个*n*维单纯形（最简单的凸形，如线段、三角形或四面体）内的随机点。它需要一个参数*points*，一个包含单纯形 n+1 个顶点的列表，所有顶点都在 n 维或更高维度。3 个顶点的特殊情况来自 Osada 等人（2002）[^100]。
METHOD VecAddProd(a, b, c)
for j in 0...size(a): a[j]=a[j]+b[j]*c
END METHOD
METHOD RandomPointInSimplex(points):
ret=NewList()
if size(points) > size(points[0])+1: return error
if size(points)==1 // Return a copy of the point
for i in 0...size(points[0]): AddItem(ret,points[0][i])
return ret
end
if size(points)==3
// Equivalent to sqrt(RNDRANGEMinMaxExc(0, 1))
rs=max(RNDRANGEMinMaxExc(0, 1), RNDRANGEMinMaxExc(0, 1))
r2=RNDRANGEMinMaxExc(0, 1)
ret=[0,0,0]
VecAddProd(ret,points[0],1.0-rs)
VecAddProd(ret,points[1],(1.0-r2)*rs)
VecAddProd(ret,points[2],r2*rs)
return ret
end
gammas=NewList()
// Sample from the simplex
for i in 0...size(points): AddItem(gammas, Expo(1))
tsum=0 // Will store sum of all gammas
for i in 0...size(gammas): tsum=tsum+gammas[i]
for i in 0...size(gammas): gammas[i] = gammas[i] / tsum
gammas[size(gammas)-1]=0 // To omit last gamma in sum
tot = 1.0 // Will store 1 minus the sum of all gammas
for i in 0...size(gammas): tot=tot - gammas[i]
// Build the final point
for i in 0...size(points[0]): AddItem(ret, points[0][i]*tot)
for i in 1...size(points): VecAddProd(
ret, points[i], gammas[i-1])
return ret
END METHOD
球体上的随机点
以下伪代码展示了如何生成一个以原点为中心、具有以下参数的球体（球的表面）上的随机点
dims，球体的维度（以及随机点的维度）。radius，球体的半径（如果radius为 1，则结果也可以作为dims维空间中的单位向量）。p大于 0，或为无穷大，它描述了球体的形状（如果p为 2，则球体是常规的）。
参见 Schechtmann 和 Zinn (1990)[^101]。这里，EPD 生成一个*指数幂*随机变量（Devroye 1986, pp. 174-175）[^19]。
METHOD PNorm(vec, p)
ret=0
if p==infinity
for i in 0...size(vec): ret=max(ret,abs(vec[i]))
return ret
else
for i in 0...size(vec): ret=ret+pow(abs(vec[i]),p)
return pow(ret,1.0/p)
end
END METHOD
METHOD EPD(p)
# Infinity case is uniform in (-1,1) to be
# appropriate for this section's purposes
if p==infinity: return RNDRANGEMinMaxExc(-1,1)
if p==2: return Normal(0,1)
return (RNDINT(1) * 2 - 1)*pow(GammaDist(1/p),1/p)
END METHOD
METHOD RandomPointOnSphere(dims, radius, p)
x=0
while x==0
ret=[]
for i in 0...dims: AddItem(ret, EPD(p))
x=PNorm(ret, p)
end
invnorm=radius/x
for i in 0...dims: ret[i]=ret[i]*invnorm
return ret
END METHOD
注释
PNorm(vec, p)，也称为 ℓp范数，是一种广义的距离概念。PNorm(vec, 2)是“常规”距离，例如，它构成了“常规”球体的基础。p可以是任何大于或等于 0 的数，也可以是无穷大。- Python 示例代码（Python sample code）包含了一种针对圆（2 维球体，
p=2）上点的优化方法。示例： 要生成圆柱体（沿 Z 轴）表面上的随机点，请在圆上生成随机的 X 和 Y 坐标，并通过
RNDRANGEMinMaxExc(mn, mx)生成随机的 Z 坐标，其中mn和mx是可能的最高和最低 Z 坐标。
盒子、球体、壳或圆锥内的随机点
生成内部的随机点——
- N 维*盒子*，为每个坐标生成
RNDRANGEMinMaxExc(mn, mx)，其中mn和mx是该坐标的下界和上界。例如——- 要生成一个矩形内的随机点，该矩形沿 X 轴在 [0, 2) 区间内，沿 Y 轴在 [3, 6) 区间内，则生成
[RNDRANGEMinMaxExc(0,2), RNDRANGEMinMaxExc(3,6)]，并且 - 生成一个实部和虚部都限制在 [0, 1] 内的*复数*，则生成
[RNDRANGEMinMaxExc(0, 1), RNDRANGEMinMaxExc(0, 1)]。
- 要生成一个矩形内的随机点，该矩形沿 X 轴在 [0, 2) 区间内，沿 Y 轴在 [3, 6) 区间内，则生成
- 以给定半径为中心的*N*维*球体*，遵循
RandomPointOnSphere中的伪代码，但将PNorm(ret, p)替换为pow(pow(PNorm(ret, p),p)+Expo(1),1.0/p)(Barthe 等人 2005)[^102]。[^103] - 以原点为中心、内半径为 A、外半径为 B（其中 A 小于 B）的*N*维*球面壳*（空心球），生成一个半径为
pow(RNDRANGEMinMaxExc(pow(A, N), pow(B, N)), 1.0 / N)的 N 维球体表面上的随机点[^104]。 - 一个高为
H、底面半径为R且沿 Z 轴的*圆锥*，通过Z = max(max(RNDRANGEMinMaxExc(0, H), RNDRANGEMinMaxExc(0, H)), RNDRANGEMinMaxExc(0, H))生成一个随机 Z 坐标，然后生成半径等于max(RNDRANGEMinMaxExc(0,Z*(R/H)), RNDRANGEMinMaxExc(0,Z*(R/H)))的圆盘（2 维球体）内的随机 X 和 Y 坐标[^105]。
示例： 要生成一个圆柱体（沿 Z 轴）内部的随机点，请在圆盘（2 维球体）内生成随机的 X 和 Y 坐标，并通过
RNDRANGEMinMaxExc(mn, mx)生成随机的 Z 坐标，其中mn和mx是可能的最高和最低 Z 坐标。注释
- Python 示例代码（Python sample code）包含了一种用于生成模拟地球的椭球体表面上的随机点的方法。
- 采样半球、半球体或半壳，可以通过采样完整球体或球壳，然后将结果的一个维度替换为其绝对值来完成。
- Lacko 和 Harman (2012)[^106] 定义了一族*非均匀*分布的球内点：生成
RandomPointOnSphere(dims, r*pow(BetaDist(dims/p, d/p), 1.0/p),p)，其中r>0是半径，dims和p与RandomPointOnSphere中的相同，d>=0是形状参数。如果d = p，则分布在球体内是均匀的。
随机纬度和经度
生成球体表面上的随机点（以纬度和经度表示，以弧度为单位，西方和南方坐标为负）[^107]——
- 生成经度
RNDRANGEMinMaxExc(-pi, pi)，其中经度在 [-π, π) 区间内，并且 - 生成纬度
atan2(sqrt(1 - x * x), x) - pi / 2，其中x = RNDRANGEMinMaxExc(-1, 1)且纬度在 [-π/2, π/2] 区间内（该区间排除了两极，两极有多种等效形式；如果不需要两极，则生成x直到以此方式不生成 -1 也不生成 1）。
致谢
我感谢该页面 CodeProject 版本的评论者，包括 George Swan，他向我推荐了蓄水池抽样方法。
我还感谢 Christoph Conrads，他在本文的某些部分提出了建议。
其他文档
以下是我撰写的一些关于随机化和伪随机变量生成的附加文章。它们都是开源的。
- 随机数生成器应用建议
- 更多随机采样方法
- 离散分布的代码生成器
- 涉及随机化最常见的主题
- 用于连续分布精确采样的部分采样随机数
- Bernoulli 工厂算法
- 测试 PRNG 的高质量随机性
- 高质量 PRNG 的示例
注释
[^1]: Pedersen, K., "Reconditioning your quantile function", arXiv:1704.07949v3 [stat.CO], 2018。
[^2]: 对于通过 RNDINT 伪代码部分解决的练习，请参见 A. Koenig 和 B. E. Moo 的《Accelerated C++》，2000 年；另请参见 Johnny Chan 的博客文章。
[^3]: 这种源的一个例子是高斯噪声发生器。这种源通常被称为*熵源*。
[^4]: Lemire, D., "A fast alternative to the modulo reduction", Daniel Lemire's blog, 2016。
[^5]: Lemire, D., "Fast Random Integer Generation in an Interval", arXiv:1805.10941v4 [cs.DS], 2018。
[^6]: Lumbroso, J., "Optimal Discrete Uniform Generation from Coin Flips, and Applications", arXiv:1304.1916 [cs.DS]
[^7]: "Probability and Random Numbers", 2004 年 2 月 29 日。
[^8]: Mennucci, A.C.G., "Bit Recycling for Scaling Random Number Generators", arXiv:1012.4290 [cs.IT], 2018。
[^9]: Knuth, Donald E. 和 Andrew Chi-Chih Yao。"The complexity of nonuniform random number generation"，载于《Algorithms and Complexity: New Directions and Recent Results》，1976 年。
[^10]: 这是因为 p = 1/n 的*二进制熵*为 p * log2(1/p) = log2(n) / n，而 n 个概率为 1/n 的结果的 n 个二进制熵之和为 log2(n) = ln(n)/ln(2)。
[^11]: Devroye, L., Gravel, C., "Random variate generation using only finitely many unbiased, independently and identically distributed random bits", arXiv:1502.02539v6 [cs.IT], 2020。
[^12]: 在某些语言（如 JavaScript）中常见的朴素 RNDINTEXC 实现是 floor(Math.random() * maxExclusive) 这个惯用语，其中 Math.random() 是任何表现得像 [0, 1) 区间内独立均匀随机变量的浮点数输出方法。然而，没有 Math.random() 的实现可以从 [0, 1) 中的所有实数中选择，因此根据 maxExclusive 的值，此惯用语可能会导致某些结果比其他结果有偏差。例如，如果 Math.random() 实现为 RNDINT(X - 1)/X 且 X 不能被 maxExclusive 整除，则结果将有偏差。此外，实现可能会在 floor 之前对 Math.random() * maxExclusive 进行四舍五入（到它可以表示的最近数字）；在极少数情况下，根据四舍五入模式，这可能是 maxExclusive。如果应用程序担心这些问题，它应该将 Math.random() 的实现视为 RNDINT 的“随机数源”，并通过 RNDINT 来实现 RNDINTEXC。
[^13]: 用户“BVtp”来自Stack Overflow社区，他让我有了这个见解。
[^14]: Daniel Ting, "Simple, Optimal Algorithms for Random Sampling Without Replacement", arXiv:2104.05091, 2021。
[^15]: Sanders, P., Lamm, S., et al., "Efficient Parallel Random Sampling – Vectorized, Cache-Efficient, and Online", arXiv:1610.0514v2 [cs.DS], 2019。
[^16]: Jeff Atwood, "The danger of naïveté", 2007 年 12 月 7 日。
[^17]: Bacher, A., Bodini, O., et al., "MergeShuffle: A Very Fast, Parallel Random Permutation Algorithm", arXiv:1508.03167 [cs.DS], 2015。
[^18]: Merlini, D., Sprugnoli, R., Verri, M.C., "An Analysis of a Simple Algorithm for Random Derangements", 2007。
[^19]: Devroye, L., Non-Uniform Random Variate Generation, 1986。
[^20]: 另见Stack Overflow问题 "Random index of a non zero value in a numpy array"。
[^21]: S. Linderman, "A Parallel Gamma Sampling Implementation", Laboratory for Independent Probabilistic Systems Blog, 2013 年 2 月 21 日，举例说明了一个 GPU 实现的 gamma 分布随机变量采样器。
[^22]: De Bruyne, B., et al., "Generating discrete-time constrained random walks and Lévy flights, arXiv:2104.06145 (2021)。
[^23]: Brownlee, J. "A Gentle Introduction to the Bootstrap Method", Machine Learning Mastery, 2018 年 5 月 25 日。
[^24]: Propp, J.G., Wilson, D.B., "Exact sampling with coupled Markov chains and applications to statistical mechanics", 1996。
[^25]: Fill, J.A., "An interruptible algorithm for perfect sampling via Markov chains", Annals of Applied Probability 8(1), 1998。
[^26]: E. N. Gilbert, "Random Graphs", Annals of Mathematical Statistics 30(4), 1959。
[^27]: V. Batagelj 和 U. Brandes, "Efficient generation of large random networks", Phys.Rev. E 71:036113, 2005。
[^28]: Costantini, Lucia. "Algorithms for sampling spanning trees uniformly at random." 硕士论文，Universitat Politècnica de Catalunya，2020 年。https://upcommons.upc.edu/bitstream/handle/2117/328169/memoria.pdf
[^29]: Penschuck, M., et al., "Recent Advances in Scalable Network Generation", arXiv:2003.00736v1 [cs.DS], 2020。
[^30]: Jon Louis Bentley 和 James B. Saxe, "Generating Sorted Lists of Random Numbers", ACM Trans. Math. Softw. 6 (1980), pp. 359-364，描述了一种生成特定类型随机变量的排序列表的方法，但此处未给出，因为它依赖于生成 [0, 1] 区间内的实数，这本身就是不完美的，因为计算机无法选择 0 和 1 之间的所有实数，而且它们是无限的。
[^31]: (2020a) Saad, F.A., Freer C.E., et al., "The Fast Loaded Dice Roller: A Near-Optimal Exact Sampler for Discrete Probability Distributions", arXiv:2003.03830v2 [stat.CO]，也发表在《AISTATS 2020: Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research》108 期，巴勒莫，西西里岛，意大利，2020 年。
[^32]: Feras A. Saad, Cameron E. Freer, Martin C. Rinard, 和 Vikash K. Mansinghka, "Optimal Approximate Sampling From Discrete Probability Distributions", arXiv:2001.04555v1 [cs.DS]，也发表在 Proc. ACM Program. Lang. 4, POPL, Article 36 (2020 年 1 月), 33 页。
[^33]: Shaddin Dughmi, Jason D. Hartline, Robert Kleinberg, and Rad Niazadeh. 2017. Bernoulli Factories and Black-Box Reductions in Mechanism Design. In Proceedings of 49th Annual ACM SIGACT Symposium on the Theory of Computing, Montreal, Canada, June 2017 (STOC’17)。
[^34]: K. Bringmann 和 K. Panagiotou, "Efficient Sampling Methods for Discrete Distributions." Algorithmica 79 (2007)，也发表在 Proc. 39th International Colloquium on Automata, Languages, and Programming (ICALP'12), 2012 年。
[^35]: A.J. Walker, "An efficient method for generating discrete random variables with general distributions", ACM Transactions on Mathematical Software 3, 1977。
[^36]: Vose, Michael D. "A linear algorithm for generating random numbers with a given distribution." IEEE Transactions on software engineering 17, no. 9 (1991): 972-975。
[^37]: Klundert, B. van de, "Efficient Generation of Discrete Random Variates", Faculty of Science Theses, Universiteit Utrecht, 2019。
[^38]: K. Bringmann 和 K. G. Larsen, "Succinct Sampling from Discrete Distributions", In: Proc. 45th Annual ACM Symposium on Theory of Computing (STOC'13), 2013。
[^39]: L. Hübschle-Schneider 和 P. Sanders, "Parallel Weighted Random Sampling", arXiv:1903.00227v2 [cs.DS], 2019。
[^40]: Y. Tang, "An Empirical Study of Random Sampling Methods for Changing Discrete Distributions", 硕士论文，University of Alberta，2019。
[^41]: T. S. Han 和 M. Hoshi, "Interval algorithm for random number generation", IEEE Transactions on Information Theory 43(2), 1997 年 3 月。
[^42]: Efraimidis, P. "Weighted Random Sampling over Data Streams", arXiv:1012.0256v2 [cs.DS], 2015。
[^43]: Efraimidis, P. 和 Spirakis, P. "Weighted Random Sampling (2005; Efraimidis, Spirakis)", 2005。
[^44]: von Neumann, J., "Various techniques used in connection with random digits", 1951。
[^45]: Deville, J.-C. 和 Tillé, Y. Unequal probability sampling without replacement through a splitting method. Biometrika 85 (1998)。
[^46]: Chao, M-T., "A general purpose unequal probability sampling plan", Biometrika 69 (1982)。
[^47]: Python 示例代码（Python sample code）包含一个 ConvexPolygonSampler 类，该类实现了此类采样（用于凸多边形）；与一般多边形不同，凸多边形可以轻松分解为三角形。
[^48]: 该文章还提到了一个暴击分布，实际上它是两种分布的*混合*：一次掷骰子和两次掷骰子的总和。
[^49]: *仿射变换*是保持直线不变且平行线保持平行的变换。
[^50]: Farach-Colton, M. 和 Tsai, M.T., 2015. Exact sublinear binomial sampling. Algorithmica 73(4), pp. 637-651。
[^51]: K. Bringmann, F. Kuhn, et al., “Internal DLA: Efficient Simulation of a Physical Growth Model.” In: Proc. 41st International Colloquium on Automata, Languages, and Programming (ICALP'14), 2014。
[^52]: Heaukulani, C., Roy, D.M., "Black-box constructions for exchangeable sequences of random multisets", arXiv:1908.06349v1 [math.PR], 2019。但请注意，此参考文献将负二项分布定义为 N 次失败前的成功次数（而非反之）。
[^53]: Bringmann, K., and Friedrich, T., 2013, July. Exact and efficient generation of geometric random variates and random graphs, in International Colloquium on Automata, Languages, and Programming (pp. 267-278)。
[^54]: Duchon, P., Duvignau, D., "Preserving the number of cycles of length k in a growing uniform permutation", Electronic Journal of Combinatorics 23(4), 2016。
[^55]: Johnson 和 Kotz, "Discrete Distributions", 1969。
[^56]: Smith, Noah A., 和 Roy W. Tromble. "Sampling uniformly from the unit simplex." Johns Hopkins University, Tech. Rep 29 (2004)。
[^57]: Durfee, et al., "l1 Regression using Lewis Weights Preconditioning and Stochastic Gradient Descent", Proceedings of Machine Learning Research 75(1), 2018。
[^58]: NVIDIA 的白皮书 "Floating Point and IEEE 754 Compliance for NVIDIA GPUs"，以及 Bruce Dawson 的 "Floating-Point Determinism"，更详细地讨论了浮点数问题。
[^59]: “Uniform”加引号表示，对于该数字格式来说，尽可能接近均匀分布。排除两个边界是因为在数学上，来自均匀分布的任何特定实数出现的概率为 0。
[^60]: Boehm, Hans-J. "Towards an API for the real numbers." In Proceedings of the 41st ACM SIGPLAN Conference on Programming Language Design and Implementation, pp. 562-576. 2020。
[^61]: 这包括整数（如果 e 限制为 0），以及定点数（如果 e 限制为小于 0 的单个指数）。
[^62]: Downey, A. B. "Generating Pseudo-random Floating Point Values", 2007。
[^63]: 理想情况下，X 是最大的整数 p，使得 [0, 1] 区间内 1/p 的所有倍数在该数字格式中都可表示。例如，对于 binary64，X 为 2^53 (9007199254740992)，对于 binary32，X 为 2^24 (16777216)。
[^64]: Goualard F. (2020) Generating Random Floating-Point Numbers by Dividing Integers: A Case Study. In: Krzhizhanovskaya V. et al. (eds) Computational Science – ICCS 2020. ICCS 2020. Lecture Notes in Computer Science, vol 12138. Springer, Cham。https://doi.org/10.1007/978-3-030-50417-5_2
[^65]: Monahan, J.F., "Accuracy in Random Number Generation", Mathematics of Computation 45(172), 1985。
[^66]: Halmos, P.R., "The theory of unbiased estimation", Annals of Mathematical Statistics 17(1), 1946。
[^67]: Spall, J.C., "An Overview of the Simultaneous Perturbation Method for Efficient Optimization", Johns Hopkins APL Technical Digest 19(4), 1998, pp. 482-492。
[^68]: Owen, Art B. "On dropping the first Sobol’point." In International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing, pp. 71-86. Springer, Cham, 2022. "Preprint": arXiv:2008.08051。
[^69]: P. L'Ecuyer, "Tables of Linear Congruential Generators of Different Sizes and Good Lattice Structure", Mathematics of Computation 68(225), 1999 年 1 月，附有勘误。
[^70]: Harase, S., "A table of short-period Tausworthe generators for Markov chain quasi-Monte Carlo", arXiv:2002.09006 [math.NA], 2020。
[^71]: Kamath, Chandrika. "Intelligent sampling for surrogate modeling, hyperparameter optimization, and data analysis." Machine Learning with Applications (2022): 100373, https://doi.org/10.1016/j.mlwa.2022.100373.
[^72]: D. Revuz, M. Yor, "Continuous Martingales and Brownian Motion", 1999。
[^73]: Lewis, P.W., Shedler, G.S., "Simulation of nonhomogeneous Poisson processes by thinning", Naval Research Logistics Quarterly 26(3), 1979。
[^74]: Saucier, R. "Computer Generation of Statistical Distributions", 2000 年 3 月。
[^75]: 其他关于密度估计的参考资料包括关于多变量核密度估计的维基百科文章，以及 M. Kay 的博客文章。
[^76]: “Jitter”（在此步骤中使用）遵循一种称为*核*的分布，其中正态分布是其中一种示例。带宽应设置为使估计的分布适合数据并保持平滑。一种更复杂的“抖动”（用于多分量数据点）包括从均值全为 0 的*多变量正态分布*生成的点，以及一个*协方差矩阵*，在此上下文中用作*带宽矩阵*。“Jitter”和带宽在本文件中不再讨论。
[^77]: *离散分布*是一种将可数数量的项目与每个项目关联的概率的分布。这些项目通常是整数，但不必是。例如，项目可以是负整数（例如，x/y 的概率为 x/(1+y)），只要它们可以转换为整数或从整数转换而来。两个示例：- 一个最简分数可以表示为整数，方法是将分子和分母的位交错。- 整数量化数字（在“深度学习”神经网络中很受欢迎）占用相对较少的位数（通常是 8 位甚至更少）。8 位量化数字格式有效地映射 256 个整数到实数。
[^78]: 这包括整数（如果 FPExponent 限制为 0），以及定点数（如果 FPExponent 限制为小于 0 的单个指数）。
[^79]: Saad, F.A., et al., "Optimal Approximate Sampling from Discrete Probability Distributions", arXiv:2001.04555 [cs.DS], 2020。另见相关的源代码。
[^80]: Walter, M., "Sampling the Integers with Low Relative Error", In: International Conference on Cryptology in Africa, Jul. 2019, pp. 157-180。
[^81]: Gerhard Derflinger, Wolfgang Hörmann, 和 Josef Leydold, "Random variate generation by numerical inversion when only the density is known", ACM Transactions on Modeling and Computer Simulation 20(4) article 18, October 2010。
[^82]: numbers_from_u01 的一部分使用了 Arnas, D., Leake, C., Mortari, D. 的算法，“Random Sampling using k-vector”，Computing in Science & Engineering 21(1) pp. 94-107, 2019，以及 Mortari, D., Neta, B., "k-Vector Range Searching Techniques"。
[^83]: Devroye, L., "Non-Uniform Random Variate Generation". In Handbooks in Operations Research and Management Science: Simulation, Henderson, S.G., Nelson, B.L. (eds.), 2006, p.83。
[^84]: Sainudiin, Raazesh, 和 Thomas L. York. "An Auto-Validating, Trans-Dimensional, Universal Rejection Sampler for Locally Lipschitz Arithmetical Expressions," Reliable Computing 18 (2013): 15-54。
[^85]: 许多马尔可夫链收敛到*平稳分布*。链的*混合时间*是收敛的速度；它是经过多少步后下一次采样将遵循一个接近平稳分布的分布。然后通过独立随机地初始化几个马尔可夫链，运行每个链直到其混合时间（“预热”），然后取每个链的下一次采样来采样近似分布。有关更多信息，请参见 Levin 和 Peres 的《Markov chains and mixing times》，第二版，2017 年。
[^86]: Tran, K.H., "A Common Derivation for Markov Chain Monte Carlo Algorithms with Tractable and Intractable Targets", arXiv:1607.01985v5 [stat.CO], 2018，它提供了一个描述许多 MCMC 算法（包括 Metropolis–Hastings、slice sampling 和 Gibbs sampling）的通用框架。
[^87]: Kschischang, Frank R. "A Trapezoid-Ziggurat Algorithm for Generating Gaussian Pseudorandom Variates." (2019)。
[^88]: Devroye, L., 1996, December, "Random variate generation in one line of code" In Proceedings Winter Simulation Conference (pp. 265-272). IEEE。
[^89]: Crooks, G.E., Field Guide to Continuous Probability Distributions, 2019。
[^90]: McGrath, E.J., Irving, D.C., "Techniques for Efficient Monte Carlo Simulation, Volume II", Oak Ridge National Laboratory, April 1975。
[^91]: Crooks, G.E., "The Amoroso Distribution", arXiv:1005.3274v2 [math.ST], 2015。
[^92]: Kotz, Samuel, Tomasz Kozubowski, 和 Krzystof Podgorski. The Laplace distribution and generalizations: a revisit with applications to communications, economics, engineering, and finance. Springer Science & Business Media, 2012。
[^93]: Devroye, L., "Expected Time Analysis of a Simple Recursive Poisson Random Variate Generator", Computing 46, pp. 165-173, 1991。
[^94]: Bailey, R.W., "Polar generation of random variates with the t distribution", Mathematics of Computation 62 (1994)。
[^95]: Stein, W.E. 和 Keblis, M.F., "A new method to simulate the triangular distribution", Mathematical and Computer Modelling 49(5-6), 2009, pp.1143-1147。
[^96]: Flajolet, P., Pelletier, M., Soria, M., "On Buffon machines and numbers", arXiv:0906.5560v2 [math.PR], 2010。
[^97]: Makalic, E., Schmidt, D.F., "An efficient algorithm for sampling from sin^k(x) for generating random correlation matrices", arXiv:1809.05212v2 [stat.CO], 2018。
[^98]: Ghorbanzadeh, D., Jaupi, L., Durand, P., "A Method to Simulate the Skew Normal Distribution", Applied Mathematics 5(13), 2014。
[^99]: Baker, R., Jackson, D., "A new distribution for robust least squares", arXiv:1408.3237 [stat.ME], 2018。
[^100]: Osada, R., Funkhouser, T., et al., "Shape Distributions", _ACM Transactions on Graphics 21(4), Oct. 2002。
[^101]: Schechtman, G., Zinn, J., On the volume of intersection of two Lp^n balls. 1990。
[^102]: Barthe, F., Guédon, O., et al., "A probabilistic approach to the geometry of the lP^N-ball", Annals of Probability 33(2), 2005。
[^103]: 或者，如果 p 是大于 0 的整数，则生成一个具有给定半径的 N+p 维球体表面的随机点（例如，使用 RandomPointOnSphere(N+p,radius,p)），然后丢弃该点的最后 p 个坐标（Lacko, V., & Harman, R. (2012) 的推论 1。A conditional distribution approach to uniform sampling on spheres and balls in Lp spaces. Metrika, 75(7), 939-951）。
[^104]: 参见Mathematics Stack Exchange问题 "Random multivariate in hyperannulus"，questions/1885630。
[^105]: 参见Stack Overflow问题 "Uniform sampling (by volume) within a cone"，questions/41749411。平方和立方根被最大值替换。
[^106]: Lacko, V., & Harman, R. (2012). A conditional distribution approach to uniform sampling on spheres and balls in Lp spaces. Metrika, 75(7), 939-951。
[^107]: 参考：MathWorld 中的 "Sphere Point Picking"（用 atan2 等效项替换了反余弦）。
[^108]: 在隐私背景下，例如，参见 Awan, J. 和 Rao, V., 2021. "Privacy-Aware Rejection Sampling", arXiv:2108.00965。
[^109]: 例如，参见 Balcer, V., Vadhan, S., "Differential Privacy on Finite Computers", 2018 年 12 月 4 日；以及 Micciancio, D. 和 Walter, M., "Gaussian sampling over the integers: Efficient, generic, constant-time", 在 Annual International Cryptology Conference, August 2017 (pp. 455-485)。
附录
随机数源
本页面上介绍的所有随机化方法都假设我们有一个无尽的数字源，其中——
- 这些数字遵循*均匀分布*，并且
- 每个数字都*独立于任何其他选择*。
也就是说，这些方法假定我们已经有一个“（均匀）随机数源”。（因此，严格来说，这些方法本身不*生成*随机数，而是假设我们已经有了一个随机数源。）
然而，这是一个理想化的假设，在实践中很难实现，甚至不可能。
事实上，大多数应用程序都使用*伪随机数生成器*（PRNGs），它们是产生*随机行为*数字的算法，也就是说，它们模拟了上面提到的理想“随机数源”的数字。因此，即使理论上这些方法与 PRNG 的质量无关，但在实践中，它们的效果和质量将取决于 PRNG（或其他随机行为数字生成器）的质量。
“随机数源”可以通过各种设备和程序来模拟，包括 PRNG、所谓的“真随机数生成器”以及为应用程序提供均匀随机行为数字的应用程序编程接口（API）。应用程序应选择能够充分模拟“随机数源”的设备或程序，包括在统计质量、“不可猜测性”或两者兼顾方面。但是，进一步的建议超出本文档的范围。
本文档中的随机化方法是确定性的（也就是说，在给定相同状态和输入的情况下，它们会产生相同的值），而不管是什么在模拟“随机数源”（例如 PRNG 或“真随机数生成器”）。例外情况如下：
- 这些方法“不知道”由“随机数源”（或模拟该源的任何东西）接下来会产生什么数字。
- 一些方法从大小未知的文件中读取行；在读取它们之前，它们不会“知道”这些行的内容。
实现注意事项
- Shell 脚本和 Microsoft Windows 批处理文件旨在运行其他程序，而不是通用编程。然而，批处理文件和
bash（一个 shell 脚本解释器）可能支持一个变量，该变量返回 [0, 32767] 区间内均匀分布的“随机”整数（分别称为%RANDOM%或$RANDOM）；这两个变量都不是为信息安全设计的。只要有可能，本文档中的方法就不应在 shell 脚本或批处理文件中实现，尤其是在以信息安全为目标时。 - SQL 等查询语言没有循环和分支等过程式元素。此外，标准 SQL 没有随机选择数字的方法，但流行的 SQL 方言经常有——其行为是特有的——并且描述 SQL 方言之间的差异超出了本文档的范围。只要有可能，本文档中的方法就不应在 SQL 中实现，尤其是在以信息安全为目标时。
- 无状态 PRNG。 目前普遍使用的伪随机数生成器（PRNG）设计大多维护一个内部状态，并在每次生成伪随机数时更新该状态。但对于*无状态* PRNG 设计（包括所谓的“可分割”PRNG），本文档中的
RNDINT()、NEXTRAND()和其他随机采样方法可能需要相应调整（通常是通过添加一个额外的参数）。 - 多线程。 多线程可以作为一次生成多个随机变量的快速方法；它没有在本页给出的伪代码中体现。通常，这涉及将一块内存分成多个块，将每个块分配给一个线程，为每个线程提供自己的伪随机数生成器实例（或其他模拟“随机数源”的程序），并让每个线程用随机变量填充其分配的块。有关示例，请参见“Multithreaded Generation”。
- 固定数量的“随机性”。 给定一个 k 位整数 n（位于 [0, 2k) 区间内，并随机均匀选择），可以通过找到 (n+1)/2k+2 的分位数，或者使用 n 来帮助初始化本地 PRNG 并使用 PRNG 从该分布生成样本来生成近似于概率分布（例如，
Poisson、Normal）的值。应用程序应仅在希望确保最终绘制的每个采样结果具有固定数量的“随机性”时才使用此建议，因为通过这种方式，采样方法只能返回 2k 种不同的结果中的一种。（一般来说，除非 n 的不同结果数量是 2 的幂，否则无法使用*固定*数量的随机选择的位来均匀随机地选择 n。）
安全注意事项
如果应用程序出于信息安全目的进行随机采样，例如随机生成密码或加密密钥，则以下规定适用：
- “加密生成器”。 应用程序必须使用一个设备或程序，该设备或程序能够生成难以从信息安全角度猜测的随机行为数字（所谓的“加密生成器”）。选择此类设备或程序超出了本文档的范围。
- 时序攻击。 某些安全和隐私攻击利用时序和其他差异来恢复明文、加密密钥或其他秘密或私有数据。因此，已开发出安全算法，使其在内存访问模式等方面没有时序差异，从而不会泄露任何关于秘密或私有输入（如密钥、密码或伪随机数生成器的“种子”）的信息。（此类算法包括所谓的“恒定时间”算法。）但是，即使算法的运行时间可变（例如，拒绝采样），它也可能（也可能不）具有安全相关的时序差异，特别是如果它不重用秘密[^108]。
- 安全算法范围之外。 使用随机秘密生成随机安全参数（如加密密钥、公私钥对、椭圆曲线或椭圆曲线上的点）的安全算法超出了本文档的范围。
- 浮点数。 为安全目的随机选择的数字几乎总是整数（极少数情况下是定点数）。即使在少数安全应用中使用浮点数（尤其是差分隐私和基于格的密码学），也有方法可以避免使用这些浮点数[^109]。