用于连续分布精确采样的部分采样随机数

引言

本页面介绍了一种部分采样随机数（PSRN）的Python实现。尽管PSRN的结构在此工作之前已被大量描述，但本文档统一了先前工作中此类数字的概念，并展示了它们如何用于采样Beta分布（对于大多数参数集）、指数分布（具有任意速率参数）以及许多其他连续分布——

• 同时避免浮点运算，并且
• 以任意精度和用户指定的误差范围（因此，在（Karney 2016）[^1]中定义的意义上是“精确”的）。

例如，这两点将本文档中的Beta采样器与我所知的任何其他专门设计的Beta采样器区分开来。至于指数分布，有论文讨论使用随机位生成指数随机变量（Flajolet and Saheb 1982）[^2]、（Karney 2016）[^1]、（Devroye and Gravel 2020）[^3]、（Thomas and Luk 2008）[^4]，但它们中的大多数（如果不是全部）都没有处理使用任意速率（而不仅仅是1）生成指数PSRN。(Habibizad Navin et al., 2010)[^5]可能是一个例外；然而，该方法似乎涉及预先生成的数字概率表。

此处讨论的采样器也借鉴了与一种称为伯努利工厂的构造相关的研究（Keane and O'Brien 1994）[^6]（Flajolet et al., 2010）[^7]，它可以在已知硬币正面朝上的概率未知的情况下模拟一个新的概率。伯努利工厂的一个重要特性是它们可以精确地模拟给定的概率，而无需手动计算该概率，这在概率可能是计算机无法精确计算的无理数（例如pow(p, 1/2)或exp(-2)）时非常重要。

本页展示了这些采样器的Python代码。

关于本文档

这是一份开源文档；如需更新版本，请参阅源代码或其在GitHub上的渲染。您可以在CodeProject或GitHub问题页面上就本文档发表评论。

目录

• 引言
  • 关于本文档
• 目录
• 关于Beta分布
• 关于指数分布
• 关于部分采样随机数
  • 均匀部分采样随机数
  • 指数部分采样随机数
  • 其他分布
  • 属性
  • 限制
  • 与构造性实数的关系
• 采样均匀和指数PSRN
  • 采样均匀PSRN
  • 采样E-rands
• PSRN的算术和比较
  • 加法和减法
  • 乘法
  • 倒数和除法
  • 使用算术算法
  • 比较
  • 讨论
• 构建块
  • SampleGeometricBag
  • FillGeometricBag
  • kthsmallest
  • 幂次方均匀子算法
• Beta和指数分布的算法
  • Beta分布
  • 指数分布
• 采样器代码
  • 指数采样器：扩展
• 正确性测试
  • Beta采样器
  • ExpRandFill
  • ExpRandLess
• 在0到1之间支持的连续分布的精确模拟
  • 一个例子：连续伯努利分布
• 复杂性
  • 一般原则
  • 特定算法的复杂性
• 在加权水库采样中的应用
• 致谢
• 其他文档
• 注释
• 附录
  • SampleGeometricBag算法的等效性
  • UniformMultiply算法
  • 均匀的均匀生成均匀的乘积
  • Oberhoff的“精确拒绝采样”方法
  • 按数字概率设置数字
• 许可证

关于Beta分布

Beta分布是一个有界域概率分布；它的两个参数alpha和beta都大于0，并描述了分布的形状。根据alpha和beta的不同，形状可以是平滑的峰值或平滑的谷值。Beta分布可以取值在区间[0, 1]内。此区间中的任何值（x）发生的概率与——

pow(x, alpha - 1) * pow(1 - x, beta - 1).               (1)

尽管alpha和beta都可以大于0，但本文档中介绍的采样器仅在以下情况下有效：

• 两个参数都大于或等于1，或
• 在二进制数的情况下，一个参数等于1，另一个参数大于0。

关于指数分布

指数分布有一个参数λ。通俗地说，遵循指数分布的随机变量是两个事件之间的时间单位数，而λ是每单位时间的预期平均事件数。通常，λ等于1。

指数随机变量通常按如下方式生成：-ln(1 - X) / lamda，其中X是区间(0, 1)内的均匀分布随机实数。（然而，这种特定的算法在实践中并不健壮，原因超出本文档的范围，但请参阅（Pedersen 2018）[^8]。）本页介绍了采样指数随机变量的另一种方法。

关于部分采样随机数

在本文档中，部分采样随机数（PSRN）是一种数据结构，用于存储无限精度的实数，但其内容仅在必要时才进行采样。PSRN为以“精确”方式遵循概率分布、具有任意精度且无需浮点运算的算法打开了大门（请参阅本节后面的“属性”）。

此处指定的PSRN由以下三部分组成

• 具有任意位数目的小数部分。这可以实现为数字数组或包含所有数字的压缩整数。有些算法关心这些数字是已采样还是未采样；在这种情况下，如果一个数字未采样，则其未采样状态可以以区别于已采样数字的方式进行标记（例如，在Python中使用None关键字，或数字-1，或通过存储一个单独的位数组来指示哪些位已采样和未采样）。存储所有数字的基数（例如十进制的基数10或二进制的基数2）是任意的。小数部分的数字形成所谓的数字展开（例如，二进制数字情况下的二进制展开）。小数部分之外的数字是未采样的。

  例如，如果小数部分按顺序存储基数10的数字[1, 3, 5]，则它表示区间[0.135, 0.136]中的随机变量，反映了0.135和0.136之间的数字是未知的。

• 一个可选的整数部分（更具体地说，是数字绝对值的整数部分，即floor(abs(x))）。
• 一个可选的符号（正或负）。

如果未存储整数部分，则假定为0。如果未存储符号，则假定为正。例如，一个实现可以通过仅存储小数部分来只关心区间[0, 1]中的PSRN。

PSRN最终表示介于另外两个数字之间的随机变量；其中一个变量的两个边界之一具有以下形式：符号 * (整数部分 + 小数部分)，如果PSRN为正，则为下限；如果为负，则为上限。例如，如果PSRN存储一个正号，整数3，以及小数部分[3, 5, 6]（以10为基数），则PSRN表示区间[3.356, 3.357]中的随机变量。这里，其中一个边界是使用PSRN的符号、整数部分和小数部分构建的，并且由于PSRN是正的，所以这是一个下限。

本节规定了两种PSRN：均匀和指数。

均匀部分采样随机数

PSRN最简单的例子是均匀分布。

• Flajolet et al. (2010)[^7]使用术语几何袋来指代区间[0, 1]中的均匀PSRN，它存储二进制（基数2）数字，其中一些可能未采样。在这种情况下，PSRN可以只包含小数部分，这可以按前面所述实现。
• （Karney 2016）[^1]使用术语u-rand来指代可以存储符号、整数部分和小数部分的均匀PSRN，其中小数部分数字的基数是任意的，但Karney的概念只考虑从左到右无间隙地采样数字。

均匀PSRN小数部分的每个附加数字都简单地通过将其设置为独立的均匀随机数字来采样，这一观察结果可以追溯到 von Neumann (1951)[^9]在二进制情况下的研究。[^10] 具有此属性的PSRN在本文档中称为均匀PSRN，即使它是使用非均匀随机采样算法（如Karney的算法用于正态分布）生成的。（这主要是因为，一般来说，这种PSRN表示给定区间内的均匀随机变量。例如，如果PSRN是3.356...，那么它表示区间[3.356, 3.357]内的均匀分布随机变量。）

指数部分采样随机数

在本文档中，一个指数PSRN（或e-rand，类似于Karney的用于部分采样均匀随机变量的“u-rands”（Karney 2016）[^1]）采样每个位，当与现有位结合时，产生给定速率的指数分布随机变量。此外，由于-ln(1 - X)（其中X是区间[0, 1]中的均匀随机变量）呈指数分布，因此e-rands也可以表示区间(0, 1]中部分采样均匀随机变量的自然对数。这里的区别在于，额外的位不是作为无偏随机位采样的，而是作为具有消失偏差的位采样的。（更具体地说，指数PSRN通常表示给定区间内的指数分布随机变量。）

e-rand的采样算法在“Beta和指数分布算法”一节中给出。

其他分布

其他分布的PSRN可以通过从均匀分布进行拒绝采样来实现。例子包括以下内容

• Beta和连续伯努利分布，如本文档后面所述。
• 标准正态分布，如（Karney 2016）[^1]通过运行Karney的算法N并均匀随机填充未采样数字所示，或如Du et al.（2020）[^11]改进该算法所示。
• 在[0, n)（不仅仅是[0, 1]）中采样均匀分布，将在后面的“采样均匀PSRN”中描述）。

对于所有这些分布，PSRN的未采样尾部数字收敛到均匀分布，这也适用于任何具有连续概率密度函数（或更一般地说，所谓的“绝对连续”[^12]分布）（Oberhoff 2018）[^13]、（Hill and Schürger 2005, 推论4.4）[^14]。

PSRN也可以通过从指数分布进行拒绝采样来实现，尽管这里没有给出具体的例子。

属性

使用PSRN从非离散分布采样的算法具有以下属性

1. 该算法仅依赖于独立且无偏随机位作为随机性来源。
2. 该算法不依赖于浮点运算或无理数或超越数的固定精度近似。（该算法可以计算收敛到无理数的近似值，只要这些近似值使用任意精度。）
3. 该算法可以使用有理数算术（例如Python中的Fraction或Ruby中的Rational），只要算术是精确的。
4. 如果算法输出一个PSRN，则由采样数字表示的数字必须遵循一个分布，该分布与理想分布的距离不超过b^−m，其中b是PSRN的基数或基（例如二进制为2），m是PSRN小数部分最右侧采样数字的位置，从1开始计数。((Devroye and Gravel 2020)[^3]建议使用 Wasserstein 距离，或“土方距离”，作为此目的的距离。)即使算法的调用者随后随机采样该PSRN的未采样数字（例如，在均匀PSRN的情况下均匀随机采样），该数字也必须以这种方式接近。
5. 如果算法随机填充PSRN的未采样小数位（例如，在均匀PSRN的情况下均匀随机填充），以便数字的小数部分有m位，则数字的分布必须保持接近理想分布，距离不超过b^−m。

注释

1. 将非离散分布的采样器转换为满足这些属性的算法并不容易。本文档后面“讨论”一节中给出了一些原因。
2. Oberhoff (2018)[^13]描述的精确拒绝采样算法产生表现得像PSRN的样本；然而，该算法不具备本节中描述的属性。这是因为该方法需要计算概率的最小值，并且在大多数情况下，实际上需要使用浮点算术（参见上面的属性2）。此外，算法的进展取决于先前采样位的数值，而不仅仅是这些位的位置，这与均匀分布和指数分布不同（另请参见（Thomas and Luk 2008）[^4]）。为完整起见，Oberhoff的方法出现在附录中。

限制

因为PSRN在某个区间内存储随机变量，所以PSRN不适合表示零体积集合中的数字。此类集合包括

• 整数或有理数集。
• 单个点集。
• 二维或更高维空间中的曲线。
• 三维或更高维空间中的曲面。

在曲线和曲面的情况下，PSRN不能直接存储曲线上或曲面上的随机点的空间坐标（因为这些坐标的精确值可能是一个计算机无法存储的无理数，并且任何区间都不能“足够紧密地”限制这些精确坐标），但PSRN可以存储间接给出该点在曲线上或曲面上的位置的上限和下限。

示例

1. 为了表示圆边上的一个点，PSRN可以通过稍后给出的2*π的RandUniformFromReal方法，存储区间[0, 2*π)中的随机变量（例如，它可以存储整数部分2和小数部分[1, 3, 5]，从而表示区间[2.135, 2.136]中的一个数字），并且以这种方式存储的数字表示相对于其起始位置的圆弧距离。关心点X和Y坐标的程序可以生成足够多的PSRN数字，以计算cos(P)和sin(P)的近似值，分别达到所需的精度，其中P是PSRN存储的数字。（然而，直接使用cos和sin等数学函数超出本文档的范围。）
2. 例1非常简单，因为区间上的每个点都均匀地映射到圆上的一个点。但这通常不成立：区间或盒子的点通常不会均匀地映射到曲线或曲面上的点。例如，取两个PSRN描述一个三维圆锥曲面的U和V坐标：U为[1.135, 1.136]，V为[0.288, 0.289]，圆锥的坐标为X = U*cos(V)，Y = U*sin(V)，Z = U。在这个例子中，PSRN形成一个映射到圆锥曲面一小部分的盒子。然而，盒子中的点不会以这种方式均匀地映射到圆锥上，因此，如果没有更多的工作（参见Williamson（1987）[^15]的一个解决方案），生成足够的数字来计算X、Y和Z到所需精度将不会从圆锥的该部分进行均匀采样。

与构造性实数的关系

部分采样随机数与处理所谓“构造性实数”或“递归实数”的一系列工作有关，或者对实数进行运算以计算精确结果的近似值到用户指定的小数位数。例如，在Hans-J. Boehm的实现中（Boehm 2020）[^16]、（Boehm 1987）[^17]，每个“构造性实数”上的运算（如加法、乘法、exp、ln等）都与一个函数f(n)相关联（其中n通常为0或更大），该函数返回一个整数m，使得abs(m/pow(2, n) - x) < 1/pow(2, n)，其中x是运算的精确结果。正如Goubault-Larrecq et al. (2021)[^18]所建议的，还可以有一个操作，它采样[0, 1]中均匀随机变量的数字，并提供对该变量近似值的访问，根据需要采样随机数字。类似地，可以定义此类操作以访问PSRN（包括均匀或指数PSRN）中存储的值的近似值，根据需要采样PSRN的数字。

采样均匀和指数PSRN

 

采样均匀PSRN

有几种算法用于采样给定另一个数字的均匀部分采样随机数。

RandUniform算法生成一个均匀分布的PSRN (a)，其大于0且小于另一个PSRN (b)的概率为1。此算法在必要时采样b的小数部分。如果已知b是实数而不是部分采样随机数，则不应使用此算法，因为此算法可能会超出b在算法开始时所具有（或似乎具有）的值；相反，应使用稍后给出的RandUniformFromReal算法。（例如，如果b是3.425...，此算法的一个可能结果是a = 3.42574...，b = 3.42575...。请注意，在此示例中，3.425...不被视为精确数字。）

1. 创建一个空的均匀PSRN a。令β为存储在b小数部分中的数字的基数（或基）（例如，二进制为2或十进制为10）。如果b的整数部分或符号未采样，或者如果b的符号为负，则返回错误。
2. （我们现在设置a的整数部分和符号。）将a的符号设置为正，并将a的整数部分设置为在[0, bi]中均匀随机选择的整数，其中bi是b的整数部分（请注意包含bi）。如果a的整数部分小于bi，则返回a。
3. （我们现在采样a的小数部分。）将i设置为0。
4. 如果b的整数部分为0，并且b的小数部分以采样的0位数字开头，则将i设置为b的小数部分开头的采样0位数字的数量。非零位数字或未采样位数字结束此序列。然后将i个0附加到a的小数部分。（例如，如果b是5.000302或4.000或0.0008，则有三个采样的0位数字开头b的小数部分，因此将i设置为3，并将三个0附加到a的小数部分。）
5. 如果a小数部分中位置i处的数字未采样，则将该位置的数字设置为均匀随机选择的基数β数字（例如，如果β为2，则为无偏随机位）。(位置从0开始，其中0是小数点后最重要的数字，1是下一个，依此类推）。
6. 如果b小数部分中位置i处的数字未采样，则根据PSRN b的类型对该位置的数字进行采样。（例如，如果b是均匀PSRN且β是2，这可以通过将该位置的数字设置为无偏随机位来完成。）
7. 如果a小数部分中位置i处的数字小于b的对应数字，则返回a。
8. 如果该数字更大，则丢弃a，然后创建一个新的空均匀PSRN a，然后转到步骤2。
9. 将i加1，然后转到步骤5。

注释

1. Karney (2014, sec. 4结尾)[^1]讨论了即使整数部分也可以部分采样，而不是像算法步骤2那样生成整个整数。然而，采纳这一建议将给上述算法增加非平凡的复杂性。
2. RandUniform算法等价于生成随机变量 (b) 和区间 [0, 1] 中均匀随机变量的乘积。
3. 如果b是带有正号，整数部分为0，且小数部分为空的均匀PSRN，则RandUniform算法等同于生成两个区间[0, 1]中的均匀随机变量的乘积。

RandUniformInRangePositive算法生成一个均匀分布的PSRN（a），其大于一个非负实数bmin且小于另一个正实数bmax的概率为1。该算法无论bmin或bmax是否已知为有理数都有效（例如，任一数字都可以是exp(-2)或ln(20)等表达式的结果），但算法指出，如果bmin或bmax已知为有理数，则可以更有效地实现。

1. 如果bmin大于或等于bmax，如果bmin小于0，或者如果bmax小于或等于0，则返回错误。
2. 创建一个空的均匀PSRN a。
3. 特殊情况：如果bmax为1且bmin为0，则将a的符号设置为正，将a的整数部分设置为0，并返回a。
4. 特殊情况：如果bmax和bmin是有理数，并且它们的每个分母都是β的幂，包括1（其中β是均匀PSRN所需的数字基数或基，例如十进制为10或二进制为2），则执行以下操作
   1. 令denom为bmax或bmin中较大的分母。
   2. 将c1设置为floor(bmax*denom)，将c2设置为floor((bmax−bmin)*denom)。
   3. 如果c2大于1，则向c1添加一个在[0, c2)中均匀随机选择的整数（请注意不包含c2）。
   4. 设d为denom的以β为底的对数（这等同于找到存储denom所需的最少β进制位数并减去1）。将c1的最低有效位传输到a的小数部分；变量d指示以这种方式传输到每个PSRN的位数。然后将a的符号设置为正，并将a的整数部分设置为floor(c1/β^d)。（例如，如果β是10，d是3，c1是7342，则将a的小数部分设置为[3, 4, 2]，并将a的整数部分设置为7。）最后，返回a。
5. 计算floor(bmax)，并将结果赋给bmaxi。同样，计算floor(bmin)，并将结果赋给bmini。
6. 如果bmini等于bmin且bmaxi等于bmax，则将a的符号设置为正，并将a的整数部分设置为在[bmini, bmaxi)中均匀随机选择的整数（请注意不包含bmaxi），然后返回a。（应该指出，确定一个实数是否等于另一个实数通常是不可判定的。）
7. （我们现在设置a的整数部分和符号。）将a的符号设置为正，并将a的整数部分设置为在区间[bmini, bmaxi]中均匀随机选择的整数（请注意包含bmaxi）。如果bmaxi等于bmax，则改为在区间[bmini, bmaxi−1]中选择整数。如果以下情况成立，则返回a——
   • a的整数部分大于bmini且小于bmaxi，或
   • bmini等于bmin，并且a的整数部分等于bmini且小于bmaxi。
8. （我们现在对a的小数部分进行采样。）将i设置为0，将istart设置为0。（然后，如果已知bmax是有理数：将bmaxf设置为bmax减去bmaxi，并且，如果已知bmin是有理数，将bminf设置为bmin减去bmini。）
9. （此步骤对于正确性并非至关重要，但有助于提高其效率。它将a的小数部分设置为bmin和bmax共享的初始数字。）如果a的整数部分等于bmini和bmaxi，则在循环中执行以下操作：1. 计算bmax和bmin小数部分中位置i处的基数β数字，并将dmax和dmin分别设置为这些数字。（如果已知bmax是有理数：通过将dmax设置为floor(bmaxf*β)并将dmin设置为floor(bminf*β)来执行此步骤。）2. 如果dmin等于dmax，则将dmin附加到a的小数部分，然后将i加1（并且，如果已知bmax和/或bmin是有理数，则将bmaxf设置为bmaxf*β−d，并将bminf设置为bminf*β−d）。否则，跳出此循环并将istart设置为i。
10. （确保小数部分大于bmin的。）将i设置为istart，然后如果a的整数部分等于bmini
    1. 计算bmin小数部分中位置i处的基数β数字，并将dmin设置为该数字。
    2. 如果a小数部分中位置i处的数字未采样，则将该位置的数字设置为均匀随机选择的基数β数字（例如，如果β为2，则为无偏随机位，或二进制）。(位置从0开始，其中0是小数点后最重要的数字，1是下一个，依此类推）。
    3. 设ad为a小数部分中位置i处的数字。如果ad大于dmin，则中止这些子步骤并转到步骤11。
    4. 丢弃a，创建一个新的空均匀PSRN a，如果ad小于dmin，则中止这些子步骤并转到步骤7。
    5. 将i加1，然后转到第一个子步骤。
11. （确保小数部分小于bmax的。）将i设置为istart，然后如果a的整数部分等于bmaxi
    1. 如果bmaxi为0且不等于bmax，并且如果a的小数部分没有数字，则在循环中执行以下操作
       1. 计算bmax小数部分中位置i处的基数β数字，并将d设置为该数字。（如果已知bmax是有理数：通过将d设置为floor(bmaxf*β)来执行此步骤。）
       2. 如果d为0，则在a的小数部分后附加一个0位数字，然后将i加1（并且，如果已知bmax为有理数，则将bmaxf设置为bmaxf*β−d）。否则，跳出此循环。
    2. 计算bmax小数部分中位置i处的基数β数字，并将dmax设置为该数字。（如果已知bmax是有理数：通过将bmaxf乘以β，然后将dmax设置为floor(bmaxf)，然后从bmaxf中减去dmax来执行此步骤。）
    3. 如果a小数部分中位置i处的数字未采样，则将该位置的数字设置为均匀随机选择的基数β数字。
    4. 设ad为a小数部分中位置i处的数字。如果ad小于dmax，则返回a。
    5. 丢弃a，创建一个新的空均匀PSRN a，如果以下情况成立，则中止这些子步骤并转到步骤7——
       • bmax已知不是有理数，并且ad大于dmax或者bmax小数部分中位置i之后的数字全部为零，或者
       • bmax已知是有理数，并且ad大于dmax或者bmaxf为0
    6. 将i加1，然后转到第二个子步骤。
12. 返回a。

RandUniformInRange算法生成一个均匀分布的PSRN (a)，其大于一个实数bmin且小于另一个实数bmax的概率为1。它适用于正实数和负实数，但为了减少混乱，它与RandUniformInRangePositive分开指定。

1. 如果bmin大于或等于bmax，则返回错误。如果bmin和bmax都大于或等于0，则返回RandUniformInRangePositive的结果。
2. 如果bmin和bmax都小于或等于0，则调用RandUniformInRangePositive，其中bmin = abs(bmax)且bmax = abs(bmin)，将结果的小数部分设置为负，并返回结果。
3. （此时，bmin小于0，bmax大于0。）将bmaxi设置为floor(bmax)（如果bmax大于或等于0）或−ceil(abs(bmax))（否则），并将bmini设置为floor(bmin)（如果bmin大于或等于0）或−ceil(abs(bmin))（否则）。（这样描述是为了防止实现者混淆floor与整数部分。）
4. 将ipart设置为在区间[bmini, bmaxi]中均匀随机选择的整数（请注意包含bmaxi）。如果bmaxi等于bmax，则改为在区间[bmini, bmaxi−1]中选择整数。
5. 如果ipart既不是bmini也不是bmaxi，则创建一个小数部分为空的均匀PSRN a；然后将a的符号设置为正（如果ipart大于或等于0）或负（否则）；然后将a的整数部分设置为abs(ipart+1)（如果ipart小于0）或ipart（否则）；然后返回a。
6. 如果ipart为bmini，则创建一个带正号、整数部分为abs(ipart+1)且小数部分为空的均匀PSRN a；然后运行URandLessThanReal，其中a = a，b = abs(bmin)。如果结果为1，则将a的符号设置为负并返回a。否则，转到步骤3。
7. 如果ipart为bmaxi，则创建一个带正号、整数部分为ipart且小数部分为空的均匀PSRN a；然后运行URandLessThanReal，其中a = a，b = bmax。如果结果为1，则返回a。否则，转到步骤3。

RandUniformFromReal算法生成一个均匀分布的PSRN (a)，其大于0且小于实数b的概率为1。它等价于RandUniformInRangePositive算法，其中a = a，bmin = 0，bmax = b。

UniformComplement算法生成1减去均匀PSRN (a) 的值，如下所示

1. 如果a的符号为负或其整数部分不是0，则返回错误。
2. 对于a小数部分中的每个已采样数字，将其设置为base−1−digit，其中digit是该数字，base是PSRN存储数字的基数，例如二进制的2。
3. 返回a。

采样E-rands

采样e-rand（指数PSRN）利用了两个观察结果（基于指数分布的参数λ）

• 当硬币翻转，正面朝上的概率为exp(-λ)时，如果正面朝上，则指数随机变量增加1。
• 如果硬币翻转，正面朝上的概率为1/(1+exp(λ/2^k))时，如果正面朝上，则指数随机变量增加2^-k，其中k > 0是整数。

Devroye和Gravel (2020)[^3]在第3.8节中已经提出了这些观察结果，但仅限于λ = 1。

为了仅使用随机位实现这些概率，采样器使用了两种算法，这两种算法都支持具有有理数值λ参数的e-rand

1. 一种用于模拟形如exp(-x/y)的概率（此处，在“伯努利工厂算法”中描述的exp(−x/y)算法）。
2. 一种用于模拟形如1/(1+exp(x/(y*pow(2, prec))))的概率（此处，在“伯努利工厂算法”中描述的LogisticExp算法）。

PSRN的算术和比较

本节描述了涉及均匀PSRN的加法、减法、乘法、倒数和除法，并讨论了涉及PSRN算术的其他方面。

加法和减法

以下算法 (UniformAdd) 展示了如何将两个在其小数部分存储相同基数（基）数字的均匀PSRN (a 和 b) 相加，并得到一个均匀PSRN作为结果。输入的PSRN可以有正号或负号，并假定它们的整数部分和符号已被采样。实现此算法的Python代码将在本文档后面给出。

1. 如果a在小数部分最后一个已采样数字之前有未采样数字，则将每个未采样数字设置为均匀随机选择的数字。对b也进行同样的操作。
2. 如果a的小数部分数字少于b（反之亦然），则采样足够的数字（通过将它们设置为均匀随机数字，例如如果a和b存储二进制或基数2数字，则为无偏随机位），以便两个PSRN的小数部分具有相同数量的数字。现在，令digitcount为a小数部分中的数字数量。
3. 如果a的符号为负，则令asign为−1，否则为1。如果b的符号为负，则令bsign为−1，否则为1。令afp为a的整数部分和小数部分打包成一个整数，如示例中所述，并令bfp为b的整数部分和小数部分以相同方式打包。（例如，如果a表示数字83.12344...，则afp为8312344。）令base为a和b存储的数字的基数，例如二进制的2或十进制的10。
4. 计算以下四个数字
   • afp*asign + bfp*bsign。
   • afp*asign + (bfp+1)*bsign。
   • (afp+1)*asign + bfp*bsign。
   • (afp+1)*asign + (bfp+1)*bsign。
5. 将minv设置为刚计算的四个数字中的最小值，将maxv设置为最大值。这些是应用区间加法到PSRN a和b结果的下限和上限。（例如，如果a是0.12344...，b是0.38925...，它们的小数部分相加形成c = 0.51269...，或区间[0.51269, 0.51271]。）然而，生成的PSRN在其区间内并非均匀分布，这是本算法其余部分将解决的问题，因为实际上，区间内数字的分布类似于两个均匀随机变量之和的分布。
6. 一旦将四个数字从低到高排序，令midmin为排序后的第二个数字，令midmax为排序后的第三个数字。
7. 将x设置为区间[0, maxv−minv)中的均匀随机整数。如果x < midmin−minv，则将dir设置为0（和密度函数左侧）。否则，如果x > midmax−minv，则将dir设置为1（和密度函数右侧）。否则，执行以下操作
   1. 将s设置为minv + x。
   2. 创建一个新的均匀PSRN，名为ret。如果s小于0，则将s设置为abs(1 + s)并将ret的符号设置为负。否则，将ret的符号设置为正。
   3. 将s的digitcount个最低有效位传输到ret的小数部分。（请注意，ret的小数部分存储从最高到最低有效位的数字。）然后将ret的整数部分设置为floor(s/base^digitcount)，然后返回ret。（例如，如果base是10，digitcount是3，s是34297，那么ret的小数部分设置为[2, 9, 7]，ret的整数部分设置为34。）
8. 如果dir是0（左侧），将pw设置为x，b设置为midmin−minv。否则（右侧），将pw设置为x−(midmax−minv)，b设置为maxv−midmax。
9. 将newdigits设置为0，然后将y设置为区间[0, b)中的均匀随机整数。
10. 如果dir是0，则将lower设置为pw。否则，将lower设置为b−1−pw。
11. 如果y小于lower，则执行以下操作
    1. 如果dir为0，则将s设置为minv，否则设置为midmax，然后将s设置为s*(base^newdigits) + pw。
    2. 创建一个新的均匀PSRN，名为ret。如果s小于0，则将s设置为abs(1 + s)并将ret的符号设置为负。否则，将ret的符号设置为正。
    3. 将s的digitcount + newdigits个最低有效数字传输到ret的小数部分，然后将ret的整数部分设置为floor(s/base^{digitcount + newdigits})，然后返回ret。
12. 如果y大于lower + 1，则转到步骤7。（这是一个拒绝事件。）
13. 将pw、y和b各自乘以base，然后向pw添加一个均匀随机选择的数字，然后向y添加一个均匀随机选择的数字，然后将newdigits加1，然后转到步骤10。

以下算法 (UniformAddRational) 展示了如何将一个均匀PSRN (a) 与一个有理数 b 相加。输入的PSRN可以有正号或负号，并假定其整数部分和符号已被采样。同样，有理数可以是正数、负数或零。实现此算法的Python代码将在本文档后面给出。

1. 令ai为a的整数部分。特殊情况
   • 如果a的符号为正且其小数部分没有采样的数字，并且b是大于或等于0的整数，则返回一个带正号、整数部分等于ai + b且小数部分为空的均匀PSRN。
   • 如果a的符号为负且其小数部分没有采样的数字，并且b是小于0的整数，则返回一个带负号、整数部分等于ai + abs(b)且小数部分为空的均匀PSRN。
   • 如果a的符号为正，整数部分为0，并且小数部分没有采样的数字，并且b是一个整数，则返回一个小数部分为空的均匀PSRN。如果b小于0，则PSRN的符号为负，整数部分为abs(b)−1。如果b大于或等于0，则PSRN的符号为正，整数部分为abs(b)。
   • 如果b为0，则返回a的副本。
2. 如果a在小数部分最后一个已采样数字之前有未采样数字，则将每个未采样数字设置为均匀随机选择的数字。现在，令digitcount为a小数部分中的数字数量。
3. 如果a的符号为负，则令asign为−1，否则为1。令base为存储在a小数部分中的数字的基数（例如，二进制为2或十进制为10）。将absfrac设置为abs(b)，然后将fraction设置为absfrac − floor(absfrac)。
4. 设afp为a的整数部分和小数部分打包成一个整数，如示例中所述。（例如，如果a表示数字83.12344...，则afp为8312344。）如果
5. 将ddc设置为base^dcount，然后将lower设置为((afp*asign)/ddc)+b（使用有理数算术），然后将upper设置为(((afp+1)*asign)/ddc)+b（同样使用有理数算术）。将minv设置为min(lower, upper)，将maxv设置为max(lower, upper)。
6. 将newdigits设置为0，然后将b设置为1，然后将ddc设置为base^dcount，然后将mind设置为floor(abs(minv*ddc))，然后将maxd设置为floor(abs(maxv*ddc))。（外边界）：然后将rvstart设置为mind−1（如果minv小于0），否则设置为mind，然后将rvend设置为maxd（如果maxv小于0），否则设置为maxd+1。
7. 将rv设置为区间[0, rvend−rvstart)中的均匀随机整数，然后将rvs设置为rv + rvstart。
8. （内边界。）如果minv小于0，则将innerstart设置为mind，否则设置为mind+1，然后如果maxv小于0，则将innerend设置为maxd−1，否则设置为maxd。
9. 如果rvs大于innerstart且小于innerend，则算法几乎完成，因此执行以下操作
   1. 创建一个空的均匀PSRN，称之为ret。如果rvs小于0，则将rvs设置为abs(1 + rvs)并将ret的符号设置为负。否则，将ret的符号设置为正。
   2. 将rvs的digitcount + newdigits个最低有效数字传输到ret的小数部分，然后将ret的整数部分设置为floor(rvs/base^{digitcount + newdigits})，然后返回ret。
10. 如果rvs小于或等于innerstart并且(rvs+1)/ddc（使用有理数算术计算）小于或等于minv，则转到步骤6。（这是一个拒绝事件。）
11. 如果rvs/ddc（使用有理数算术计算）大于或等于maxv，则转到步骤6。（这是一个拒绝事件。）
12. 将newdigits加1，然后将ddc、rvstart、rv和rvend各自乘以base，然后将mind设置为floor(abs(minv*ddc))，然后将maxd设置为floor(abs(maxv*ddc))，然后向rv添加一个均匀随机选择的数字，然后将rvs设置为rv+rvstart，然后转到步骤8。

乘法

以下算法 (UniformMultiplyRational) 展示了如何将均匀PSRN (a) 乘以非零有理数 b。输入的PSRN可以有正号或负号，并假定其整数部分和符号已被采样。实现此算法的Python代码将在本文档后面给出。

1. 如果a在小数部分最后一个已采样数字之前有未采样数字，则将每个未采样数字设置为均匀随机选择的数字。现在，令digitcount为a小数部分中的数字数量。
2. 创建一个均匀PSRN，称之为ret。如果a的符号为正且b小于0，或者a的符号为负且b大于或等于0，则将ret的符号设置为−1，否则设置为1，然后将ret的整数部分设置为0。令base为存储在a小数部分中的数字的基数（例如，二进制为2或十进制为10）。将absfrac设置为abs(b)，然后将fraction设置为absfrac − floor(absfrac)。
3. 设afp为a的整数部分和小数部分打包成一个整数，如示例中所述。（例如，如果a表示数字83.12344...，则afp为8312344。）
4. 将dcount设置为digitcount，然后将ddc设置为base^dcount，然后将lower设置为(afp/ddc)*absfrac（使用有理数算术），然后将upper设置为((afp+1)/ddc)*absfrac（同样使用有理数算术）。
5. 将rv设置为区间[floor(lower*ddc), floor(upper*ddc))中的均匀随机整数。
6. 将rvlower设置为rv/ddc（作为有理数），然后将rvupper设置为(rv+1)/ddc（作为有理数）。
7. 如果rvlower大于或等于lower且rvupper小于upper，则算法几乎完成，因此执行以下操作：将rv的dcount个最低有效数字传输到ret的小数部分（请注意，ret的小数部分存储从最高到最低有效位的数字），然后将ret的整数部分设置为floor(rv/base^dcount)，然后返回ret。（例如，如果base是10，dcount是4，rv是342978，那么ret的小数部分设置为[2, 9, 7, 8]，ret的整数部分设置为34。）
8. 如果rvlower大于upper或者rvupper小于lower，则转到步骤4。
9. 将rv和ddc各自乘以base，然后将dcount加1，然后向rv添加一个均匀随机选择的数字，然后转到步骤6。

附录中给出了另一个算法 (UniformMultiply)，用于将两个均匀PSRN (a 和 b) 相乘——该算法相当复杂，可能更简单的方法是将PSRN与实现任意精度乘法的“构造性实数”（前文已述）连接起来。

注意：设b>0，c≥0，d>0为有理数，其中d>c。要生成两个均匀变量的乘积，一个在[0, b]中，另一个在[c, d]中，可以使用以下算法。
(1) 使用参数b*(d−c)生成一个均匀PSRN，称之为K；
(2) 获取参数为K和b*c的UniformAddRational的结果，称之为M；
(3) 使用参数M生成一个均匀PSRN；返回该PSRN。
广义地说：“生成一个均匀(0, b*(d−c))的随机变量X，然后返回一个均匀(0, X+b*c)的随机变量”。有关此算法（至少当c = 0时）有效的证据，请参见附录。

倒数和除法

以下算法 (UniformReciprocal) 生成1/a，其中a是一个均匀PSRN，并生成一个具有该倒数的新均匀PSRN。输入的PSRN可以有正号或负号，并假定其整数部分和符号已被采样。这里提到的所有除法都应使用有理数算术。实现此算法的Python代码将在本文档后面给出。

1. 如果a在小数部分最后一个已采样数字之前有未采样数字，则将每个未采样数字设置为均匀随机选择的数字。现在，令digitcount为a小数部分中的数字数量。
2. 创建一个均匀PSRN，称之为ret。将ret的符号设置为a的符号。令base为存储在a小数部分中的数字的基数（例如，二进制为2或十进制为10）。
3. 如果a的小数部分没有非零数字，并且整数部分为0，则向a的小数部分附加一个均匀随机选择的数字。如果该数字为0，则重复此步骤。（此步骤对于当两个PSRN的区间都覆盖数字0时的正确性至关重要，因为它们的乘积分布与通常情况不同。）
4. 设afp为a的整数部分和小数部分打包成一个整数，如示例中所述。（例如，如果a表示数字83.12344...，则afp为8312344。）
5. （计算1/a的下限和上限，不考虑a的符号。）将dcount设置为digitcount，然后将ddc设置为base^dcount，然后将lower设置为(ddc/(afp+1))，然后将upper设置为(ddc/afp)。
6. 将lowerdc设置为floor(lower*ddc)。如果lowerdc为0，则将dcount加1，将ddc乘以base，然后重复此步骤。（此步骤对于正确性也很重要。）
7. （rv表示下限和上限之间的紧密区间，rv2表示均匀(0, 1)随机变量，用于与倒数的密度函数进行比较。）将rv设置为区间[lowerdc, floor(upper*ddc))中的均匀随机整数。将rv2设置为区间[0, lowerdc)中的均匀随机整数。
8. （获取紧密区间rv的边界。）将rvlower设置为rv/ddc，然后将rvupper设置为(rv+1)/ddc。
9. 如果rvlower大于或等于lower且rvupper小于upper
   1. 将rvd设置为lowerdc/ddc，然后将rvlower2设置为rv2/lowerdc，然后将rvupper2设置为(rv2+1)/lowerdc。（rvlower2和rvupper2是均匀(0, 1)变量的边界。）
   2. （与上限密度比较：y²/x²，其中y是1/a的下限，x在下限和上限之间。）如果rvupper2小于(rvd*rvd)/(rvupper*rvupper)，则算法几乎完成，因此执行以下操作：将rv的dcount个最低有效数字传输到ret的小数部分（请注意，ret的小数部分存储从最高到最低有效位的数字），然后将ret的整数部分设置为floor(rv/base^dcount)，然后返回ret。（例如，如果base是10，dcount是4，rv是342978，那么ret的小数部分设置为[2, 9, 7, 8]，ret的整数部分设置为34。）
   3. （与下限密度比较。）如果rvlower2大于(rvd*rvd)/(rvlower*rvlower)，则中止这些子步骤并转到步骤5。（这是一个拒绝事件。）
10. 如果rvlower大于upper或者rvupper小于lower，则转到步骤5。（这是一个拒绝事件。）
11. 将rv、rv2、lowerdc和ddc各自乘以base，然后将dcount加1，然后向rv添加一个均匀随机选择的数字，然后向rv2添加一个均匀随机选择的数字，然后转到步骤8。

有了这个算法，描述将一个均匀PSRN a 除以另一个均匀PSRN b 的算法 (这里称为 UniformDivide) 现在变得轻而易举。

1. 对b运行UniformReciprocal算法，创建一个新的均匀PSRN c。
2. 如果c的小数部分没有数字或全部为零，则向c附加均匀随机数字，直到附加了一个非零数字。
3. 按顺序对a和b运行UniformMultiply算法（在附录中给出），并返回该算法的结果。

同样，描述将一个均匀PSRN a 乘以一个非零有理数 b 的算法（这里称为 UniformDivideRational）也是轻而易举的。

1. 如果b为0，则返回错误。
2. 按顺序对a和1/b运行UniformMultiplyRational算法，并返回该算法的结果。

使用算术算法

上述用于加法和乘法的算法对于缩放和移位PSRN非常有用。例如，它们可以将正态分布的PSRN转换为具有任意均值和标准差的PSRN（通过首先将PSRN乘以标准差，然后加上均值）。下面是实现此过程的草图，给定两个参数，location和scale，它们都是有理数。

1. 生成一个均匀PSRN，然后通过使用从均匀分布拒绝采样的算法（例如Karney的标准正态分布算法（Karney 2016）[^1]）将其转换为所需分布的变量。此过程不适用于指数PSRN（e-rands）。
2. 运行UniformMultiplyRational算法，将均匀PSRN乘以有理参数scale，得到一个新的均匀PSRN。
3. 运行UniformAddRational算法，将新的均匀PSRN与有理参数location相加，得到第三个均匀PSRN。返回这个第三个PSRN。

另请参阅本文后面的“讨论”部分。

比较

两个不同分布但存储相同基数（基）数字的PSRN可以使用与本节中类似的算法进行精确比较。

RandLess算法比较两个PSRN a和b（并根据需要从它们采样额外的位），如果a以概率1小于b，则返回1，否则返回0（另请参阅（Karney 2016）[^1]）。

1. 如果a的整数部分尚未采样，则根据a的PSRN类型采样a的整数部分。对b也进行同样的操作。
2. 如果a的符号与b的符号不同，则如果a的符号为负，则返回1；如果为正，则返回0。如果a的符号为正，则如果a的整数部分小于b的，则返回1；如果大于，则返回0。如果a的符号为负，则如果a的整数部分小于b的，则返回0；如果大于，则返回1。
3. 将i设置为0。
4. 如果a小数部分中位置i处的数字未采样，则根据a的PSRN类型设置该位置的数字。（位置从0开始，其中0是小数点后最重要的数字，1是下一个，依此类推。）对b也进行同样的操作。
5. 令da为a小数部分中位置i处的数字，令db为b的对应数字。
6. 如果a的符号为正，则如果da小于db，则返回1；如果da大于db，则返回0。
7. 如果a的符号为负，则如果da小于db，则返回0；如果da大于db，则返回1。
8. 将i加1，然后转到步骤4。

URandLess是RandLess的一个版本，它涉及两个均匀PSRN。URandLess算法在步骤4中通过将位置i处的数字设置为均匀随机选择的数字来采样数字i。（例如，如果a是一个存储基数2或二进制数字的均匀PSRN，这可以通过将该位置的数字设置为无偏随机位来完成。）

注意：要从区间 [0, 1] 中抽取两个均匀随机变量的最大值，或抽取区间 [0, 1] 中一个均匀随机变量的平方根：(1) 生成两个具有正号、整数部分为 0 且小数部分为空的均匀 PSRN a 和 b。(2) 按此顺序在 a 和 b 上运行 RandLess。如果调用返回 0，则返回 a；否则，返回 b。

RandLessThanReal 算法比较 PSRN a 与实数 b，如果 a 以概率 1 小于 b，则返回 1，否则返回 0。此算法根据需要对 a 的小数部分的数字进行采样。无论 b 是否已知为有理数（例如，b 可以是 exp(-2) 或 ln(20) 等表达式的结果），此算法都有效，但该算法指出如果 b 已知为有理数，则可以更有效地实现。

1. 如果 a 的整数部分或符号未采样，则返回错误。
2. 如果 b 小于 0，则将 bs 设置为 -1，否则设置为 1。计算 floor(abs(b))，并将 bi 设置为结果。（如果 b 已知为有理数： 则将 bf 设置为 abs(b) 减去 bi。）
3. 如果 a 的符号与 bs 的符号不同，如果 a 的符号为负则返回 1，如果为正则返回 0。如果 a 的符号为正，如果 a 的整数部分小于 bi 则返回 1，如果大于则返回 0。（如果两者相等则继续。）如果 a 的符号为负，如果 a 的整数部分小于 bi 则返回 0，如果大于则返回 1。（如果两者相等则继续。）
4. 将i设置为0。
5. 如果 a 的小数部分位置 i 处的数字未采样，则根据 a 是哪种 PSRN 类型设置该位置的数字。（位置从 0 开始，其中 0 是小数点后最有效数字，1 是下一个，依此类推。）
6. 计算 b 的小数部分位置 i 处的基数 β 数字，并将 d 设置为该数字。（如果 b 已知为有理数： 通过将 bf 乘以 β，然后将 d 设置为 floor(bf)，然后从 bf 中减去 d 来执行此步骤。）
7. 令 ad 为 a 的小数部分位置 i 处的数字。
8. 如果出现以下情况，则返回 1：
   • ad 小于 d 且 a 的符号为正，
   • ad 大于 d 且 a 的符号为负，或者
   • ad 等于 d，a 的符号为负，并且——
     • b 不已知为有理数 并且 b 的小数部分位置 i 之后的所有数字都为零（表示 a 以概率 1 小于 b），或者
     • b 已知为有理数 并且 bf 为 0（表示 a 以概率 1 小于 b）。
9. 如果出现以下情况，则返回 0：
   • ad 小于 d 且 a 的符号为负，
   • ad 大于 d 且 a 的符号为正，或者
   • ad 等于 d，a 的符号为正，并且——
     • b 不已知为有理数 并且 b 的小数部分位置 i 之后的所有数字都为零（表示 a 以概率 1 大于 b），或者
     • b 已知为有理数 并且 bf 为 0（表示 a 以概率 1 大于 b）。
10. 将i加1，然后转到步骤5。

上述算法中步骤 6 到 9 的替代版本如下（另请参阅 (Brassard et al. 2019)[^19]）：

• (6.) 计算 bp，它是 b 的一个近似值，使得 abs(b − bp) <= β^{−i − 1}，并且 bp 具有与 b 相同的符号。令 bk 为 bp 的数字展开，直到小数点后 i + 1 位（忽略其符号）。例如，如果 b 是 π 或 -π，β 是 10，i 是 4，则一种可能性是 bp = 3.14159 且 bk = 314159。
• (7.) 令 ak 为 a 的数字展开，直到小数点后 i + 1 位（忽略其符号）。
• (8.) 如果 ak <= bk − 2，如果 a 的符号为正则返回 1，否则返回 0。
• (9.) 如果 ak >= bk + 1，如果 a 的符号为负则返回 1，否则返回 0。

URandLessThanReal 是 RandLessThanReal 的一个版本，其中 a 是一个均匀 PSRN。URandLessThanReal 的算法在步骤 4 中通过将位置 i 处的数字设置为均匀随机选择的数字来采样数字 i。

以下展示了当 b 是一个已知分子和分母 num/den 的分数时，如何实现 URandLessThanReal。

1. 如果 a 的整数部分或符号未采样，或者如果 den 为 0，则返回错误。然后，如果 num 和 den 都小于 0，则将它们设置为其绝对值。然后，如果 a 的符号为正，其整数部分为 0，并且 num 为 0，则返回 0。然后，如果 a 的符号为正，其整数部分为 0，并且 num 的符号与 den 的符号不同，则返回 0。
2. 如果 num 或 den（但不是两者）小于 0，则将 bs 设置为 -1，否则设置为 1，然后将 den 设置为 abs(den)，然后将 bi 设置为 floor(abs(num)/den)，然后将 num 设置为 rem(abs(num), den)。
3. 如果 a 的符号与 bs 的符号不同，如果 a 的符号为负则返回 1，如果为正则返回 0。如果 a 的符号为正，如果 a 的整数部分小于 bi 则返回 1，如果大于则返回 0。（如果两者相等则继续。）如果 a 的符号为负，如果 a 的整数部分小于 bi 则返回 0，如果大于则返回 1。（如果两者相等则继续。）如果 num 为 0（表示分数为整数），如果 a 的符号为正则返回 0，否则返回 1。
4. 将 pt 设置为 base，并将 i 设置为 0。（base 是 a 的数字的基数，例如二进制为 2，十进制为 10。）
5. 将 d1 设置为 a 的小数部分第 i 个位置（从 0 开始）的数字。如果该位置的数字未采样，则在该位置放置一个均匀随机选择的数字，并将 d1 设置为该数字。
6. 如果 num * pt >= den，则将 c 设置为 1，否则设置为 0。
7. 将 d2 设置为 floor(num * pt / den)。（在基数 2 中，这等效于将 d2 设置为 c。）
8. 如果 d1 小于 d2，如果 a 的符号为正则返回 1，否则返回 0。如果 d1 大于 d2，如果 a 的符号为正则返回 0，否则返回 1。
9. 如果 c 为 1，则将 num 设置为 num * pt - den * d2，然后将 den 乘以 pt。
10. 如果 num 为 0，如果 a 的符号为正则返回 0，否则返回 1。
11. 将 pt 乘以 base，将 i 加 1，然后转到步骤 5。

讨论

本节讨论 PSRN 算术中涉及的问题。

均匀 PSRN 算术通常会产生非均匀分布。正如本节前面的算术算法（例如 UniformAdd 和 UniformMultiplyRational）所示，PSRN 的加法、乘法和其他算术运算（另请参阅 (Brassard et al., 2019)[^19]）并不像对其整数部分和小数部分进行加法、乘法等操作那样简单。均匀 PSRN 最终是一个区间内的均匀随机变量（这是它的性质），但随机变量的算术运算通常不会产生均匀分布。

一个例子说明了这一点。假设我们有两个均匀 PSRN：A = 0.12345... 和 B = 0.38901...。它们分别表示区间 AI = [0.12345, 0.12346] 和 BI = [0.38901, 0.38902] 中的随机变量。添加两个均匀 PSRN 类似于添加它们的区间（使用区间算术），因此在此示例中，结果 C 位于 CI = [0.12345 + 0.38901, 0.12346 + 0.38902] = [0.51246, 0.51248]。然而，结果随机变量不是在 [0.51246, 0.51248] 中均匀分布的，因此简单地在该区间中选择一个均匀随机变量是行不通的。（对于加法之外的其他算术运算，这通常也是正确的。）这可以通过生成许多对在区间 AI 和 BI 中的均匀随机变量，对每对中的数字求和，并使用和构建直方图（所有和都将位于区间 CI 中）来证明。在这种情况下，直方图将显示一个三角形分布，在 0.51247 处达到峰值。

该示例通常适用于加法之外的大多数其他数学运算（包括乘法、除法、log、sin 等）：对区间 AI 和 BI 执行数学运算，并构建位于结果区间中的随机结果（乘积、商等）的直方图，以找出由此形成的分布。

实现其他运算。与加法、乘法和除法相反，某些其他数学运算在 PSRN 中执行起来微不足道。它们包括求反，如 (Karney 2016)[^1] 中所述，以及仅影响 PSRN 整数部分的运算。

将 PSRN 与其他数学运算（例如乘法、ln 和 exp）连接起来的一种有前景的方法是使用“构造实数”或“递归实数”。请参阅前面“与构造实数的关系”一节。

可以创建一个采样器，它使用目标分布下每个数字出现的概率。但如果分布是非离散的

• 除了非常有限的几类分布（包括均匀分布和指数分布）之外，这些概率将取决于前面的数字；详细信息请参阅附录。
• 对于超出该有限类别的分布，采样器在实践中将是有限精度的（而不是任意精度的），因为它只能保留如此多的数字概率。例如，(Habibizad Navin et al., 2007)[^20] 和 (Nezhad et al., 2013)[^21] 的著作指出要构建这样一个数字概率的“树”。[^22]

最后，如果算术结果以已知概率密度函数 (PDF) 分布，则 PSRN 的算术可能是可行的，从而允许实现从均匀或指数分布中拒绝的算法。这方面的一个例子可以在上面的 UniformReciprocal 算法中找到，或者在我的关于均匀随机变量和的任意精度采样器的文章中找到。然而，该 PDF 可能具有无界峰值，从而实际上排除了拒绝采样。例如，如果 X 是区间 [0, 1] 中的均匀 PSRN，则 X³ 的分布具有 PDF (1/3) / pow(X, 2/3)，其在 0 处具有无界峰值。虽然这在实践中排除了 X³ 的简单拒绝采样器，但仍然可以使用 PSRN 对均匀分布的幂进行采样，本文稍后将对此进行描述。

PSRN 的重用。如果同一 PSRN 在这些算法的不同运行中多次使用，本节中的算术算法可能会给出不正确的结果。

这个问题通常发生在原始采样器在算法中将同一随机变量用于不同目的时（一个例子是“W*Y, (1−W)*Y”，其中 W 和 Y 是独立的随机变量 (Devroye 1986, p. 394)[^23]）。在这种情况下，如果一个 PSRN 产生额外的 PSRN（从而它们变得依赖于第一个 PSRN），那么一旦第一个 PSRN 的额外数字被均匀随机采样，这些额外的 PSRN 可能会变得不准确。（这并非总是如此，但很难描述何时这些额外的 PSRN 会以这种方式变得不准确，何时不会。）

对于 UniformAddRational 或 UniformMultiplyRational 算法，当它与相同的 PSRN 和相同的有理数多次调用时，这个问题很容易看到：尽管每次都应该返回相同的随机变量，但在现实中，这种方式将以概率 1 生成不同的变量，尤其是在之后从它们中采样额外的数字时。

人们可能会认为，刚才描述的问题可以通过下面的算法解决

假设我们想用不同的数字乘以同一个 PSRN。令 vec 为要乘以同一个 PSRN 的有理数向量，令 vec[i] 为向量中位置 i 处的有理数（位置从 0 开始）。

1. 将 i 设置为 0，将 a 设置为输入 PSRN，将 num 设置为 vec[i]，将 'output' 设置为空列表。
2. 将 ret 设置为 UniformMultiplyRational 与 PSRN a 和有理数 num 的结果。
3. 将指向 ret 的指针添加到列表 'output'。如果 vec[i] 是向量中的最后一个数字，则停止此算法。
4. 将 a 指向 ret，然后将 i 加 1，然后将 num 设置为 vec[i]/vec[i-1]，然后转到步骤 2。

然而，即使此算法也不能确保输出 PSRN 将与同一随机变量完全成比例。例如：令 a 为 PSRN 0....（或区间 [0.0, 1.0]），然后令 b 为 UniformMultiplyRational(a, 1/2) 的结果，然后令 c 为 UniformMultiplyRational(b, 1/3) 的结果。b 的一个可能结果是 0.41...，c 的一个可能结果是 0.138...。现在我们用均匀随机位填充 a、b 和 c。因此，作为一个可能的结果，a 现在是 0.13328133...，b 现在是 0.41792367...，c 现在是 0.13860371...。然而，这里 c 除以 b 不完全是 1/3，尽管它非常接近，而 b 除以 a 远非 1/2（特别是由于 a 最初非常粗糙）。尽管此示例显示了具有十进制数字的 PSRN，但对于二进制数字，情况更糟。

构建块

本文档依赖于本节中描述的几个构建块。

其中之一是 Flajolet 等人 (2010)[^7] 的“几何袋”技术，该技术以逐位构建的概率生成正面或反面。

SampleGeometricBag

SampleGeometricBag 算法以均匀 PSRN 小数部分构建的概率返回 1。（Flajolet 等人，2010 年）[^7] 描述了基数 2（二进制）情况的算法，但该算法难以应用于其他数字基数。因此，以下是该算法的通用版本，适用于任何数字基数。为方便起见，此算法忽略了 PSRN 的整数部分和符号。

1. 将 i 设置为 0，并将 b 设置为一个具有正号且整数部分为 0 的均匀 PSRN。
2. 如果输入 PSRN 小数部分位置 i 处的项目未采样（即未设置为数字），则将该位置的项目设置为均匀随机选择的数字，必要时增加小数部分的容量（位置从 0 开始，其中 0 是小数点后最有效数字，1 是下一个，依此类推），并将结果附加到该小数部分的数字展开中。对 b 执行相同的操作。
3. 令 da 为输入 PSRN 小数部分位置 i 处的数字，令 db 为 b 的相应数字。如果 da 小于 db 则返回 0，如果 da 大于 db 则返回 1。
4. 将 i 加 1 并转到步骤 2。

对于基数 2，可以使用以下 SampleGeometricBag 算法，它更接近 Flajolet 论文中给出的算法。它同样忽略了输入 PSRN 的整数部分和符号。

1. 将 N 设置为 0。
2. 以 1/2 的概率，转到下一步。否则，将 N 加 1 并重复此步骤。（当算法转到下一步时，N 是 Flajolet 论文中称为几何随机变量（参数为 1/2）的东西，因此得名“几何袋”，但本文避免使用“几何随机变量”这一术语，因为它在学术著作中有几种相互冲突的含义。）
3. 如果均匀 PSRN 小数部分（位置从 0 开始）位置 N 处的项目未设置为数字（例如，基数 2 为 0 或 1），则将该位置的项目设置为均匀随机选择的数字（例如，基数 2 为 0 或 1），必要时增加小数部分的容量。（由于此步骤，均匀 PSRN 中可能存在尚未采样的数字的“间隙”。）
4. 返回位置 N 处的项目。

有关这两种算法为何等效的更多信息，请参阅附录。

SampleGeometricBagComplement 与 SampleGeometricBag 算法相同，除了返回值是原始返回值减 1。结果是如果 SampleGeometricBag 以概率 U 输出 1，则 SampleGeometricBagComplement 以概率 1 - U 输出 1。

FillGeometricBag

FillGeometricBag 接受一个均匀 PSRN 并生成一个小数部分具有 p 位数的数字，如下所示：

1. 对于 [0, p) 中的每个位置，如果均匀 PSRN 小数部分中该位置的项目未采样，则将该项目设置为均匀随机选择的数字（例如，对于二进制，为 0 或 1），必要时增加小数部分的容量。（位置从 0 开始，其中 0 是小数点后最有效数字，1 是下一个，依此类推。另请参阅 (Oberhoff 2018, sec. 8)[^13]。）
2. 如果 PSRN 的符号为负，则令 sign 为 -1，否则为 1；令 ipart 为 PSRN 的整数部分；令 bag 为 PSRN 的小数部分。取 bag 的前 p 位数字并返回 sign * (ipart + ∑_{i=0, ..., p−1} bag[i] * b⁻ⁱ⁻¹)，其中 b 是基数。

在步骤 2 之后，如果 PSRN 小数部分中超过 p 的数字已经采样（即它们已经设置为一个数字），则实现可以选择填充第一个和最后一个设置数字之间所有未采样的数字，并返回完整的数字，可以选择将其四舍五入到小数部分具有 p 位数的数字，并选择一种舍入模式。（例如，如果 p 为 4，b 为 10，PSRN 为 0.3437500... 或 0.3438500...，它可以使用四舍五入模式将 PSRN 分别四舍五入到数字 0.3438 或 0.3439；因为这是一个具有“无限”但未采样数字展开的 PSRN，所以此处不应用诸如“偶数舍入”之类的平局打破。）

kthsmallest

kthsmallest 方法生成区间 [0, 1] 中 'n' 个 'bitcount' 位均匀随机变量中的第 'k' 个最小的（也称为第 'n' 个顺序统计量），此 Beta 采样器也依赖于此方法。当 a 和 b 都是整数时，它被使用，基于一个已知属性，即在这种情况下，Beta 随机变量是区间 [0, 1] 中 'a + b - 1' 个均匀随机变量中的第 'a' 个最小的（Devroye 1986, p. 431）[^23]。

然而，kthsmallest 并非简单地生成 'n' 个 'bitcount' 位数字然后对它们进行排序。相反，它通过 PSRN 逐位构建它们的数字展开。它利用了这样的观察：(在二进制情况下) 区间 [0, 1] 中的每个均匀随机变量小于一半或大于一半的概率相等；因此，小于一半与大于一半的均匀数字数量遵循二项式(n, 1/2) 分布（并且在小于一半的数字中，例如，小于四分之一与大于四分之一的数字遵循相同的分布，依此类推）。由于这一观察，该算法可以“即时”生成一个排序样本。如果使用多项式分布而不是二项式分布，类似的观察也适用于基数 2 以外的其他基数。我不知道有任何其他文章或论文（除了我的一篇）描述此处给出的 kthsmallest 算法。

算法如下

1. 创建 n 个具有正号且整数部分为 0 的均匀 PSRN。
2. 将 index 设置为 1。
3. 如果 index <= k 且 index + n >= k
   1. 生成 v，一个多项式随机向量，其中 b 个概率都等于 1/b，其中 b 是基数（对于二进制情况，b = 2，因此这等效于生成 LC，一个参数为 n 和 0.5 的二项式随机变量，并将 v 设置为 {LC, n - LC}）。
   2. 从 index 开始，将数字 0 附加到前 v[0] 个 PSRN，将数字 1 附加到接下来的 v[1] 个 PSRN，依此类推，直到将数字 b − 1 附加到最后 v[b − 1] 个 PSRN（对于二进制情况，这意味着将 0 位附加到前 LC 个 PSRN，将 1 位附加到接下来的 n - LC 个 PSRN）。
   3. 对于 [0, b) 中的每个整数 i：如果 v[i] > 1，则使用 index = index + ∑_{j=0, ..., i−1} v[j] 和 n = v[i] 重复步骤 3 和这些子步骤。（对于二进制情况，这意味着：如果 LC > 1，则使用相同的 index 和 n = LC 重复步骤 3 和这些子步骤；然后，如果 n - LC > 1，则使用 index = index + LC 和 n = n - LC 重复步骤 3 和这些子步骤）。
4. 取第 k 个 PSRN（从 1 开始），然后根据需要用均匀随机数字填充它，使其小数部分有 bitcount 位数字（类似于上面的 FillGeometricBag），然后返回该数字。（请注意，后面描述的 Beta 采样器选择通过此算法以这种方式填充 PSRN。）

幂次方均匀子算法

幂-均匀子算法用于以下 Beta 采样器的某些情况。它返回 U^px/py，其中 U 是区间 [0, 1] 中的均匀随机变量，px/py 大于 1，但与朴素算法不同，它支持任意精度，仅使用随机位，并避免浮点运算。它还使用一个补码标志来确定是否返回 1 减去结果。

它利用了以下一些算法

• 它使用了一个采样无界单调 PDF的算法，该算法又类似于 (Devroye 1986, ch. 7, sec. 4.4)[^23] 中的逆变换拒绝算法。这是必需的，因为当 px/py 大于 1 时，U^px/py 的分布具有 PDF (py/px) / pow(U, 1-py/px)，其在 0 处具有无界峰值。
• 它使用了许多伯努利工厂算法，包括 SampleGeometricBag 和在“伯努利工厂算法”中描述的一些算法。

然而，此算法目前仅支持生成小数部分为基数 2（二进制）数字的 PSRN。

幂-均匀算法如下

1. 将 i 设置为 1。
2. 调用“伯努利工厂算法”中描述的算法 (a/b)^x/y，参数为 a = 1, b = 2, x = py, y = px。如果调用返回 1 且 i 小于 n，则将 i 加 1 并重复此步骤。如果调用返回 1 且 i 等于或大于 n，如果 complement 标志为 1 则返回 1，否则返回 0（或返回一个具有正号、整数部分为 0 且小数部分分别填充精确 n 个 1 或 0 的均匀 PSRN）。
3. 结果是，我们现在将采样一个位于区间 [2⁻ⁱ, 2^{−(i − 1)}) 中的数字。我们现在必须在该区间中生成一个均匀随机变量 X，然后以概率 (py / (px * 2ⁱ)) / X^{1 − py / px} 接受它；此公式中的 2ⁱ 是为了帮助避免采样目的的极低概率。以下步骤将无需使用浮点运算即可实现这一点。
4. 创建一个正号零整数部分均匀 PSRN，然后创建一个 geobag 输入硬币，该硬币返回在该 PSRN 上调用 SampleGeometricBag 的结果。
5. 创建一个 powerbag 输入硬币，该硬币执行以下操作：“调用‘伯努利工厂算法’中描述的算法 λ^x/y，使用 geobag 输入硬币，且 x/y = 1 − py / px，然后返回结果。”
6. 在 PSRN 的小数部分附加 i − 1 个零位，然后是单个一位。这将允许我们采样一个限制在前面提到的区间内的均匀随机变量。
7. 调用“伯努利工厂算法”中描述的算法 ϵ / λ，使用 powerbag 输入硬币（表示 b），且 ϵ = py/(px * 2ⁱ)（表示 a），从而以概率 a/b 返回 1。如果调用返回 1，则接受 PSRN，因此执行以下操作：
   1. 如果 complement 标志为 1，则将 PSRN 小数部分中的每个零位变为一位，反之亦然。
   2. 可选地，根据需要用均匀随机数字填充 PSRN，使其小数部分有 n 位数字（类似于上面的 FillGeometricBag），其中 n 是精度参数。然后，返回 PSRN。
8. 如果对算法 ϵ / λ 的调用返回 0，则从 PSRN 的小数部分中删除除前 i 位之外的所有数字，然后转到步骤 7。

Beta和指数分布的算法

 

Beta分布

现在所有构建块都已就位，可以描述一种采样 Beta 分布的新算法，描述如下。它有三个参数：a >= 1 和 b >= 1（或在二进制情况下一个参数为 1 且另一个大于 0）是 Beta 分布的参数，p > 0 是一个精度参数。

1. 特殊情况
   • 如果 a = 1 且 b = 1，则返回一个正号零整数部分均匀 PSRN。
   • 如果 a 和 b 都是整数，则返回 kthsmallest 的结果，其中 n = a - b + 1 且 k = a
   • 在二进制情况下，如果 a 为 1 且 b 小于 1，则调用下面描述的幂-均匀子算法，其中 px/py = 1/b，并且 complement 标志设置为 1，并按原样返回该算法的结果（不按该算法的子步骤 7.2 所述进行填充）。
   • 在二进制情况下，如果 b 为 1 且 a 小于 1，则调用下面描述的幂-均匀子算法，其中 px/py = 1/a，并且 complement 标志设置为 0，并按原样返回该算法的结果（不按该算法的子步骤 7.2 所述进行填充）。
2. 如果 a > 2 且 b > 2，则执行以下步骤，将 a 和 b 分成两部分，这两部分模拟速度更快（并实现 (Devroye 1986, page 47 顶部)[^23] 中的广义拒绝策略）
   1. 将 aintpart 设置为 floor(a) − 1，将 bintpart 设置为 floor(b) − 1，将 arest 设置为 a − aintpart，并将 brest 设置为 b − bintpart。
   2. 单独（递归）运行此算法，但 a = aintpart 且 b = bintpart。将 bag 设置为运行创建的 PSRN。
   3. 创建一个输入硬币 geobag，该硬币返回使用给定 PSRN 调用 SampleGeometricBag 的结果。创建另一个输入硬币 geobagcomp，该硬币返回使用给定 PSRN 调用 SampleGeometricBagComplement 的结果。
   4. 调用“伯努利工厂算法”中描述的算法 λ^x/y，使用 geobag 输入硬币，且 x/y = arest/1，然后使用 geobagcomp 输入硬币和 x/y = brest/1 调用相同的算法。如果两次调用都返回 1，则返回 bag。否则，转到子步骤 2。
3. 创建一个正号零整数部分均匀 PSRN。创建一个输入硬币 geobag，该硬币返回使用给定 PSRN 调用 SampleGeometricBag 的结果。创建另一个输入硬币 geobagcomp，该硬币返回使用给定 PSRN 调用 SampleGeometricBagComplement 的结果。
4. 从 PSRN 的小数部分中删除所有数字。这将导致在以下步骤中，区间 [0, 1] 中有一个“空”的均匀随机变量 U，该步骤将以概率 U^a−1*(1−U)^b−1)（Beta 分布的比例概率）接受 U，因为 U 正在构建中。
5. 调用“伯努利工厂算法”中描述的算法 λ^x/y，使用 geobag 输入硬币，且 x/y = (a − 1)/1（因此以概率 U^a−1 返回）。如果结果为 0，则转到步骤 4。
6. 调用使用 geobagcomp 输入硬币和 x/y = (b − 1)/1 的相同算法（因此以概率 (1−U)^b−1 返回 1）。如果结果为 0，则转到步骤 4。（请注意，此步骤和上一步不相互依赖，可以按任何顺序执行，而不会影响正确性，Python 代码中利用了这一点。）
7. U 已被接受，因此返回 FillGeometricBag 的结果。

一旦 PSRN 被上述步骤接受，可选地用均匀随机数字填充 PSRN 小数部分的未采样数字，以使该数字具有 p 位小数部分（类似于上面的 FillGeometricBag），然后返回生成的数字。

注释

• 参数为 1/x 和 1 的 Beta 随机变量与区间 [0, 1] 中一个均匀随机变量的 x 次幂相同。
• 对于 Beta 分布，alpha 或 beta 越大，接受区域越小（随机变量被步骤 5 和 6 拒绝的可能性越大，从而增加其运行时间）。这是因为 PDF 的峰值 max(u^(alpha-1)*(1-u)^(beta-1)) 随着参数的增大而趋近于 0。为了解决这个问题，加入了步骤 2，在某些情况下，它将 PDF 分成两部分，这两部分模拟起来相对简单（在位复杂度方面）。

指数分布

我们还具备了描述如何采样 e-rands 的必要构建块。在 Python 代码中实现，e-rand 由五个数字组成：第一个是 1/(2^x) 的倍数，第二个是 x，第三个是整数部分（最初为 -1，表示整数部分尚未采样），第四个和第五个分别是 λ 参数的分子和分母。（因为指数随机变量总是大于等于 0，所以 e-rand 的符号隐式为正）。

要采样指数随机变量小数点后第 k 位（其中 k = 1 表示小数点后第一位，k = 2 表示第二位，依此类推），调用 LogisticExp 算法，参数为 x = λ 的分子，y = λ 的分母，prec = k。

ExpRandLess 算法是前面给出的通用 RandLess 算法的一个特例。它比较两个 e-rands a 和 b（并根据需要从它们中采样额外的位），如果 a 结果小于 b 则返回 1，否则返回 0。（请注意，a 和 b 允许具有不同的 λ 参数。）

1. 如果 a 的整数部分尚未采样，则调用算法 exp(−x/y)，其中 x = λ 的分子，y = λ 的分母，直到调用返回 0，然后将整数部分设置为此方式返回 1 的次数。对 b 执行相同的操作。
2. 如果 a 的整数部分小于 b 的整数部分，则返回 1，如果 a 的整数部分大于 b 的整数部分，则返回 0。
3. 将i设置为0。
4. 如果 a 的小数部分有 i 位或更少位，则调用 LogisticExp 算法，参数为 x = λ 的分子，y = λ 的分母，以及 prec = i + 1，并将结果附加到该小数部分的二进制展开中。（例如，如果实现将二进制展开存储为打包整数和大小，则实现可以将打包整数左移 1 位，将算法结果添加到该整数，然后将大小加 1。）对 b 执行相同的操作。
5. 如果 a 的小数部分小于 b 的小数部分，则返回 1，如果 a 的小数部分大于 b 的小数部分，则返回 0。
6. 将i加1，然后转到步骤4。

ExpRandFill 算法接受一个 e-rand 并生成一个小数部分具有 p 位数的数字，如下所示：

1. 对于 [0, p) 中的每个位置 i，如果 e-rand 小数部分中该位置的项目未采样，则调用 LogisticExp 算法，参数为 x = λ 的分子，y = λ 的分母，以及 prec = i + 1，并将位置 i 处的项目设置为结果（这将是 0 或 1），必要时增加小数部分的容量。（位位置从 0 开始，其中 0 是小数点后最有效位，1 是下一个，依此类推。另请参阅 (Oberhoff 2018, sec. 8)[^13]。）
2. 如果 e-rand 的符号为负，则令 sign 为 -1，否则为 1；令 ipart 为 e-rand 的整数部分；令 bag 为 PSRN 的小数部分。取 bag 的前 p 位数字并返回 sign * (ipart + ∑_{i=0, ..., p−1} bag[i] * 2⁻ⁱ⁻¹)。

有关如何处理在此算法步骤 2 之后，如果 PSRN 小数部分中超过 p 的位不知何故已采样（即，它们已设置为数字）的情况，请参阅 FillGeometricBag 中的讨论。

这是第三种算法（称为 ExpRand），它生成遵循指数分布的均匀 PSRN，而不是 e-rand。在该算法中，速率 λ 以大于 0 的有理数形式给出。该方法基于冯·诺依曼的算法 (von Neumann 1951)[^9]。

1. 将 recip 设置为 1/λ，并将 highpart 设置为 0。
2. 将 u 设置为 RandUniformFromReal 的结果，参数为 recip。
3. 将 val 指向与 u 相同的值，并将 accept 设置为 1。
4. 将 v 设置为 RandUniformFromReal 的结果，参数为 recip。
5. 按顺序在 u 和 v 上运行 URandLess 算法。如果调用返回 0，则将 u 设置为 v，然后将 accept 设置为 1 减去 accept，然后转到步骤 4。
6. 如果 accept 为 1，则通过前面给出的 UniformAddRational 算法将 highpart 添加到 val 中，然后返回 val。
7. 将 recip 添加到 highpart 并转到步骤 2。

前面算法的以下替代版本（称为 ExpRand2）包含了 Karney 对冯·诺依曼算法的改进 (Karney 2016)[^1]，即所谓的“早期拒绝步骤”。这里的算法允许任意速率参数 (λ)，以大于 0 的有理数形式给出，这与冯·诺依曼和 Karney 算法不同，其中 λ 为 1。

1. 将 recip 设置为 1/λ，并将 highpart 设置为 0。
2. 将 u 设置为 RandUniformFromReal 的结果，参数为 recip。
3. 在 u 上运行 URandLessThanReal 算法，参数为 recip/2。如果调用返回 0，则将 recip/2 添加到 highpart 并转到步骤 2。（这是 Karney 的“早期拒绝步骤”，当 λ 为 1 时，参数为 1/2。然而，Fan 等人 (2019)[^24] 指出 Karney 的“早期拒绝步骤”中的参数 1/2 不是最优的。）
4. 将 val 指向与 u 相同的值，并将 accept 设置为 1。
5. 将 v 设置为 RandUniformFromReal 的结果，参数为 recip。
6. 按顺序在 u 和 v 上运行 URandLess 算法。如果调用返回 0，则将 u 设置为 v，然后将 accept 设置为 1 减去 accept，然后转到步骤 5。
7. 如果 accept 为 1，则通过前面给出的 UniformAddRational 算法将 highpart 添加到 val 中，然后返回 val。
8. 将 recip/2 添加到 highpart 并转到步骤 2。

采样器代码

以下 Python 代码实现了刚刚描述的 Beta 采样器。它依赖于我编写的两个 Python 模块：

• “bernoulli.py”，它收集了许多伯努利工厂，其中一些被下面的代码所依赖。
• “randomgen.py”，它收集了许多随机变量生成方法，包括 kthsmallest 以及 RandomGen 类。

请注意，代码仅使用浮点算术将采样器的结果转换为方便的形式，即浮点数。

然而，至少在 Python 中，这段代码远非快速。

import math
import random
import bernoulli
from randomgen import RandomGen
from fractions import Fraction

def _toreal(ret, precision):
        # NOTE: Although we convert to a floating-point
        # number here, this is not strictly necessary and
        # is merely for convenience.
        return ret*1.0/(1<<precision)

def _urand_to_geobag(bag):
  return [(bag[0]>>(bag[1]-1-i))&1 for i in range(bag[1])]

def _power_of_uniform_greaterthan1(bern, power, complement=False, precision=53):
    return bern.fill_geometric_bag(
        _power_of_uniform_greaterthan1_geobag(bern, power, complement), precision
    )

def _power_of_uniform_greaterthan1_geobag(bern, power, complement=False, precision=53):
   if power<1:
     raise ValueError("Not supported")
   if power==1:
        return []  # Empty uniform random variate
   i=1
   powerfrac=Fraction(power)
   powerrest=Fraction(1) - Fraction(1)/powerfrac
   # Choose an interval
   while bern.zero_or_one_power_ratio(1,2,
         powerfrac.denominator,powerfrac.numerator) == 1:
      i+=1
   epsdividend = Fraction(1)/(powerfrac * 2**i)
   # -- A choice for epsdividend which makes eps_div
   # -- much faster, but this will require floating-point arithmetic
   # -- to calculate "**powerrest", which is not the focus
   # -- of this article.
   # probx=((2.0**(-i-1))**powerrest)
   # epsdividend=Fraction(probx)*255/256
   bag=[]
   gb=lambda: bern.geometric_bag(bag)
   bf =lambda: bern.power(gb, powerrest.numerator, powerrest.denominator)
   while True:
      # Limit sampling to the chosen interval
      bag.clear()
      for k in range(i-1):
         bag.append(0)
      bag.append(1)
      # Simulate epsdividend / x**(1-1/power)
      if bern.eps_div(bf, epsdividend) == 1:
          # Flip all bits if complement is true
          bag=[x if x==None else 1-x for x in bag] if complement else bag
          return bag

def powerOfUniform(b, px, py, precision=53):
        # Special case of beta, returning power of px/py
        # of a uniform random variate, provided px/py
        # is in (0, 1].
        return betadist(b, py, px, 1, 1, precision)

    return b.fill_geometric_bag(
        betadist_geobag(b, ax, ay, bx, by), precision
    )

def betadist_geobag(b, ax=1, ay=1, bx=1, by=1):
    """ Generates a beta-distributed random variate with arbitrary
          (user-defined) precision.  Currently, this sampler only works if (ax/ay) and
          (bx/by) are both 1 or greater, or if one of these parameters is
         1 and the other is less than 1.
         - b: Bernoulli object (from the "bernoulli" module).
         - ax, ay: Numerator and denominator of first shape parameter.
         - bx, by: Numerator and denominator of second shape parameter.
         - precision: Number of bits after the point that the result will contain.
        """
    # Beta distribution for alpha>=1 and beta>=1
    bag = []
    afrac=(Fraction(ax) if ay==1 else Fraction(ax, ay))
    bfrac=(Fraction(bx) if by==1 else Fraction(bx, by))
    bpower = bfrac - 1
    apower = afrac - 1
    # Special case for a=b=1
    if bpower == 0 and apower == 0:
        return bag
    # Special case if a=1
    if apower == 0 and bpower < 0:
        return _power_of_uniform_greaterthan1_geobag(b, Fraction(by, bx), True)
    # Special case if b=1
    if bpower == 0 and apower < 0:
        return _power_of_uniform_greaterthan1_geobag(b, Fraction(ay, ax), False)
    if apower <= -1 or bpower <= -1:
        raise ValueError
    # Special case if a and b are integers
    if int(bpower) == bpower and int(apower) == apower:
        a = int(afrac)
        b = int(bfrac)
        return _urand_to_geobag(randomgen.RandomGen().kthsmallest_psrn(a + b - 1, a))
    # Split a and b into two parts which are relatively trivial to simulate
    if bfrac > 2 and afrac > 2:
        bintpart = int(bfrac) - 1
        aintpart = int(afrac) - 1
        brest = bfrac - bintpart
        arest = afrac - aintpart
        # Generalized rejection method, p. 47
        while True:
           bag = betadist_geobag(b, aintpart, 1, bintpart, 1)
           gb = lambda: b.geometric_bag(bag)
           gbcomp = lambda: b.geometric_bag(bag) ^ 1
           if (b.power(gbcomp, brest)==1 and \
              b.power(gb, arest)==1):
              return bag
    # Create a "geometric bag" to hold a uniform random
    # number (U), described by Flajolet et al. 2010
    gb = lambda: b.geometric_bag(bag)
    # Complement of "geometric bag"
    gbcomp = lambda: b.geometric_bag(bag) ^ 1
    bp1=lambda: (1 if b.power(gbcomp, bpower)==1 and \
            b.power(gb, apower)==1 else 0)
    while True:
        # Create a uniform random variate (U) bit-by-bit, and
        # accept it with probability U^(a-1)*(1-U)^(b-1), which
        # is the unnormalized PDF of the beta distribution
        bag.clear()
        if bp1() == 1:
            # Accepted
            return ret

def _fill_geometric_bag(b, bag, precision):
        ret = 0
        lb = min(len(bag), precision)
        for i in range(lb):
            if i >= len(bag) or bag[i] == None:
                ret = (ret << 1) | b.rndint(1)
            else:
                ret = (ret << 1) | bag[i]
        if len(bag) < precision:
            diff = precision - len(bag)
            ret = (ret << diff) | b.rndint((1 << diff) - 1)
        # Now we have a number that is a multiple of
        # 2^-precision.
        return ret / (1 << precision)

def exprandless(a, b):
        """ Determines whether one partially-sampled exponential number
           is less than another; returns
           true if so and false otherwise.  During
           the comparison, additional bits will be sampled in both numbers
           if necessary for the comparison. """
        # Check integer part of exponentials
        if a[2] == -1:
            a[2] = 0
            while zero_or_one_exp_minus(a[3], a[4]) == 1:
                a[2] += 1
        if b[2] == -1:
            b[2] = 0
            while zero_or_one_exp_minus(b[3], b[4]) == 1:
                b[2] += 1
        if a[2] < b[2]:
            return True
        if a[2] > b[2]:
            return False
        index = 0
        while True:
            # Fill with next bit in a's exponential number
            if a[1] < index:
                raise ValueError
            if b[1] < index:
                raise ValueError
            if a[1] <= index:
                a[1] += 1
                a[0] = logisticexp(a[3], a[4], index + 1) | (a[0] << 1)
            # Fill with next bit in b's exponential number
            if b[1] <= index:
                b[1] += 1
                b[0] = logisticexp(b[3], b[4], index + 1) | (b[0] << 1)
            aa = (a[0] >> (a[1] - 1 - index)) & 1
            bb = (b[0] >> (b[1] - 1 - index)) & 1
            if aa < bb:
                return True
            if aa > bb:
                return False
            index += 1

def zero_or_one(px, py):
        """ Returns 1 at probability px/py, 0 otherwise.
            Uses Bernoulli algorithm from Lumbroso appendix B,
            with one exception noted in this code. """
        if py <= 0:
            raise ValueError
        if px == py:
            return 1
        z = px
        while True:
            z = z * 2
            if z >= py:
                if random.randint(0,1) == 0:
                    return 1
                z = z - py
            # Exception: Condition added to help save bits
            elif z == 0: return 0
            else:
                if random.randint(0,1) == 0:
                   return 0

def zero_or_one_exp_minus(x, y):
        """ Generates 1 with probability exp(-px/py); 0 otherwise.
               Reference: Canonne et al. 2020. """
        if y <= 0 or x < 0:
            raise ValueError
        if x==0: return 1
        if x > y:
            xf = int(x / y)  # Get integer part
            x = x % y  # Reduce to fraction
            if x > 0 and zero_or_one_exp_minus(x, y) == 0:
                return 0
            for i in range(xf):
                if zero_or_one_exp_minus(1, 1) == 0:
                    return 0
            return 1
        r = 1
        ii = 1
        while True:
            if zero_or_one(x, y*ii) == 0:
                return r
            r=1-r
            ii += 1

# Example of use
def exprand(lam):
   return exprandfill(exprandnew(lam),53)*1.0/(1<<53)

在下面的 Python 代码中，add_psrns 是一个分别生成两个均匀 PSRN 相乘或相加结果的方法。

def psrn_reciprocal(psrn1, digits=2):
    """ Generates the reciprocal of a partially-sampled random number.
        psrn1: List containing the sign, integer part, and fractional part
            of the first PSRN.  Fractional part is a list of digits
            after the point, starting with the first.
        digits: Digit base of PSRNs' digits.  Default is 2, or binary. """
    if psrn1[0] == None or psrn1[1] == None:
        raise ValueError
    for i in range(len(psrn1[2])):
        psrn1[2][i] = (
            random.randint(0, digits - 1) if psrn1[2][i] == None else psrn1[2][i]
        )
    digitcount = len(psrn1[2])
    # Perform multiplication
    frac1 = psrn1[1]
    for i in range(digitcount):
        frac1 = frac1 * digits + psrn1[2][i]
    while frac1 == 0:
        # Avoid degenerate cases
        d1 = random.randint(0, digits - 1)
        psrn1[2].append(d1)
        frac1 = frac1 * digits + d1
        digitcount += 1
    while True:
        dcount = digitcount
        ddc = digits ** dcount
        small = Fraction(ddc, frac1 + 1)
        large = Fraction(ddc, frac1)
        if small>large: raise ValueError
        if small==0: raise ValueError
        while True:
           dc = int(small * ddc)
           if dc!=0: break
           dcount+=1
           ddc*=digits
        if dc == 0:
             print(["dc",dc,"dc/ddc",float(Fraction(dc,ddc)),"small",float(small),"dcount",dcount,"psrn",psrn1])
        dc2 = int(large * ddc) + 1
        rv = random.randint(dc, dc2 - 1)
        rvx = random.randint(0, dc - 1)
        # print([count,float(small), float(large),dcount, dc/ddc, dc2/ddc])
        while True:
            rvsmall = Fraction(rv, ddc)
            rvlarge = Fraction(rv + 1, ddc)
            if rvsmall >= small and rvlarge < large:
                rvd = Fraction(dc, ddc)
                rvxf = Fraction(rvx, dc)
                rvxf2 = Fraction(rvx + 1, dc)
                # print(["dcs",rvx,"rvsmall",float(rvsmall),"rvlarge",float(rvlarge),"small",float(small),
                #   "rvxf",float(rvxf),float(rvxf2),"rvd",float(rvd),
                #   "sl",float((rvd*rvd)/(rvlarge*rvlarge)),float((rvd*rvd)/(rvsmall*rvsmall))])
                if rvxf2 < (rvd * rvd) / (rvlarge * rvlarge):
                    cpsrn = [1, 0, [0 for i in range(dcount)]]
                    cpsrn[0] = psrn1[0]
                    sret = rv
                    for i in range(dcount):
                        cpsrn[2][dcount - 1 - i] = sret % digits
                        sret //= digits
                    cpsrn[1] = sret
                    return cpsrn
                elif rvxf > (rvd * rvd) / (rvsmall * rvsmall):
                    break
            elif rvsmall > large or rvlarge < small:
                break
            rv = rv * digits + random.randint(0, digits - 1)
            rvx = rvx * digits + random.randint(0, digits - 1)
            dcount += 1
            ddc *= digits
            dc *= digits

def multiply_psrn_by_fraction(psrn1, fraction, digits=2):
    """ Multiplies a partially-sampled random number by a fraction.
        psrn1: List containing the sign, integer part, and fractional part
            of the first PSRN.  Fractional part is a list of digits
            after the point, starting with the first.
        fraction: Fraction to multiply by.
        digits: Digit base of PSRNs' digits.  Default is 2, or binary. """
    if psrn1[0] == None or psrn1[1] == None:
        raise ValueError
    fraction = Fraction(fraction)
    for i in range(len(psrn1[2])):
        psrn1[2][i] = (
            random.randint(0, digits - 1) if psrn1[2][i] == None else psrn1[2][i]
        )
    digitcount = len(psrn1[2])
    # Perform multiplication
    frac1 = psrn1[1]
    fracsign = -1 if fraction < 0 else 1
    absfrac = abs(fraction)
    for i in range(digitcount):
        frac1 = frac1 * digits + psrn1[2][i]
    while True:
        dcount = digitcount
        ddc = digits ** dcount
        small = Fraction(frac1, ddc) * absfrac
        large = Fraction(frac1 + 1, ddc) * absfrac
        dc = int(small * ddc)
        dc2 = int(large * ddc) + 1
        rv = random.randint(dc, dc2 - 1)
        while True:
            rvsmall = Fraction(rv, ddc)
            rvlarge = Fraction(rv + 1, ddc)
            if rvsmall >= small and rvlarge < large:
                cpsrn = [1, 0, [0 for i in range(dcount)]]
                cpsrn[0] = psrn1[0] * fracsign
                sret = rv
                for i in range(dcount):
                    cpsrn[2][dcount - 1 - i] = sret % digits
                    sret //= digits
                cpsrn[1] = sret
                return cpsrn
            elif rvsmall > large or rvlarge < small:
                break
            else:
                rv = rv * digits + random.randint(0, digits - 1)
                dcount += 1
                ddc *= digits

def add_psrns(psrn1, psrn2, digits=2):
    """ Adds two uniform partially-sampled random numbers.
        psrn1: List containing the sign, integer part, and fractional part
            of the first PSRN.  Fractional part is a list of digits
            after the point, starting with the first.
        psrn2: List containing the sign, integer part, and fractional part
            of the second PSRN.
        digits: Digit base of PSRNs' digits.  Default is 2, or binary. """
    if psrn1[0] == None or psrn1[1] == None or psrn2[0] == None or psrn2[1] == None:
        raise ValueError
    for i in range(len(psrn1[2])):
        psrn1[2][i] = (
            random.randint(0, digits - 1) if psrn1[2][i] == None else psrn1[2][i]
        )
    for i in range(len(psrn2[2])):
        psrn2[2][i] = (
            random.randint(0, digits - 1) if psrn2[2][i] == None else psrn2[2][i]
        )
    while len(psrn1[2]) < len(psrn2[2]):
        psrn1[2].append(random.randint(0, digits - 1))
    while len(psrn1[2]) > len(psrn2[2]):
        psrn2[2].append(random.randint(0, digits - 1))
    digitcount = len(psrn1[2])
    if len(psrn2[2]) != digitcount:
        raise ValueError
    # Perform addition
    frac1 = psrn1[1]
    frac2 = psrn2[1]
    for i in range(digitcount):
        frac1 = frac1 * digits + psrn1[2][i]
    for i in range(digitcount):
        frac2 = frac2 * digits + psrn2[2][i]
    small = frac1 * psrn1[0] + frac2 * psrn2[0]
    mid1 = frac1 * psrn1[0] + (frac2 + 1) * psrn2[0]
    mid2 = (frac1 + 1) * psrn1[0] + frac2 * psrn2[0]
    large = (frac1 + 1) * psrn1[0] + (frac2 + 1) * psrn2[0]
    minv = min(small, mid1, mid2, large)
    maxv = max(small, mid1, mid2, large)
    # Difference is expected to be a multiple of two
    if abs(maxv - minv) % 2 != 0:
        raise ValueError
    vs = [small, mid1, mid2, large]
    vs.sort()
    midmin = vs[1]
    midmax = vs[2]
    while True:
        rv = random.randint(0, maxv - minv - 1)
        if rv < 0:
            raise ValueError
        side = 0
        start = minv
        if rv < midmin - minv:
            # Left side of sum density; rising triangular
            side = 0
            start = minv
        elif rv >= midmax - minv:
            # Right side of sum density; falling triangular
            side = 1
            start = midmax
        else:
            # Middle, or uniform, part of sum density
            sret = minv + rv
            cpsrn = [1, 0, [0 for i in range(digitcount)]]
            if sret < 0:
                sret += 1
                cpsrn[0] = -1
            sret = abs(sret)
            for i in range(digitcount):
                cpsrn[2][digitcount - 1 - i] = sret % digits
                sret //= digits
            cpsrn[1] = sret
            return cpsrn
        if side == 0:  # Left side
            pw = rv
            b = midmin - minv
        else:
            pw = rv - (midmax - minv)
            b = maxv - midmax
        newdigits = 0
        y = random.randint(0, b - 1)
        while True:
            lowerbound = pw if side == 0 else b - 1 - pw
            if y < lowerbound:
                # Success
                sret = start * (digits ** newdigits) + pw
                cpsrn = [1, 0, [0 for i in range(digitcount + newdigits)]]
                if sret < 0:
                    sret += 1
                    cpsrn[0] = -1
                sret = abs(sret)
                for i in range(digitcount + newdigits):
                    idx = (digitcount + newdigits) - 1 - i
                    while idx >= len(cpsrn[2]):
                        cpsrn[2].append(None)
                    cpsrn[2][idx] = sret % digits
                    sret //= digits
                cpsrn[1] = sret
                return cpsrn
            elif y > lowerbound + 1:  # Greater than upper bound
                # Rejected
                break
            pw = pw * digits + random.randint(0, digits - 1)
            y = y * digits + random.randint(0, digits - 1)
            b *= digits
            newdigits += 1

def add_psrn_and_fraction(psrn, fraction, digits=2):
    if psrn[0] == None or psrn[1] == None:
        raise ValueError
    fraction = Fraction(fraction)
    fracsign = -1 if fraction < 0 else 1
    absfrac = abs(fraction)
    origfrac = fraction
    isinteger = absfrac.denominator == 1
    # Special cases
    # positive+pos. integer or negative+neg. integer
    if ((fracsign < 0) == (psrn[0] < 0)) and isinteger and len(psrn[2]) == 0:
        return [fracsign, psrn[1] + int(absfrac), []]
    # PSRN has no fractional part, fraction is integer
    if (
        isinteger
        and psrn[0] == 1
        and psrn[1] == 0
        and len(psrn[2]) == 0
        and fracsign < 0
    ):
        return [fracsign, int(absfrac) - 1, []]
    if (
        isinteger
        and psrn[0] == 1
        and psrn[1] == 0
        and len(psrn[2]) == 0
        and fracsign > 0
    ):
        return [fracsign, int(absfrac), []]
    if fraction == 0:  # Special case of 0
        return [psrn[0], psrn[1], [x for x in psrn[2]]]
    # End special cases
    for i in range(len(psrn[2])):
        psrn[2][i] = random.randint(0, digits - 1) if psrn[2][i] == None else psrn[2][i]
    digitcount = len(psrn[2])
    # Perform addition
    frac1 = psrn[1]
    frac2 = int(absfrac)
    fraction = absfrac - frac2
    for i in range(digitcount):
        frac1 = frac1 * digits + psrn[2][i]
    for i in range(digitcount):
        digit = int(fraction * digits)
        fraction = (fraction * digits) - digit
        frac2 = frac2 * digits + digit
    ddc = digits ** digitcount
    small = Fraction(frac1 * psrn[0], ddc) + origfrac
    large = Fraction((frac1 + 1) * psrn[0], ddc) + origfrac
    minv = min(small, large)
    maxv = max(small, large)
    while True:
        newdigits = 0
        b = 1
        ddc = digits ** digitcount
        mind = int(minv * ddc)
        maxd = int(maxv * ddc)
        rvstart = mind - 1 if minv < 0 else mind
        rvend = maxd if maxv < 0 else maxd + 1
        rv = random.randint(0, rvend - rvstart - 1)
        rvs = rv + rvstart
        if rvs >= rvend:
            raise ValueError
        while True:
            rvstartbound = mind if minv < 0 else mind + 1
            rvendbound = maxd - 1 if maxv < 0 else maxd
            if rvs > rvstartbound and rvs < rvendbound:
                sret = rvs
                cpsrn = [1, 0, [0 for i in range(digitcount + newdigits)]]
                if sret < 0:
                    sret += 1
                    cpsrn[0] = -1
                sret = abs(sret)
                for i in range(digitcount + newdigits):
                    idx = (digitcount + newdigits) - 1 - i
                    cpsrn[2][idx] = sret % digits
                    sret //= digits
                cpsrn[1] = sret
                return cpsrn
            elif rvs <= rvstartbound:
                rvd = Fraction(rvs + 1, ddc)
                if rvd <= minv:
                    # Rejected
                    break
                else:
                    # print(["rvd",rv+rvstart,float(rvd),float(minv)])
                    newdigits += 1
                    ddc *= digits
                    rvstart *= digits
                    rvend *= digits
                    mind = int(minv * ddc)
                    maxd = int(maxv * ddc)
                    rv = rv * digits + random.randint(0, digits - 1)
                    rvs = rv + rvstart
            else:
                rvd = Fraction(rvs, ddc)
                if rvd >= maxv:
                    # Rejected
                    break
                else:
                    newdigits += 1
                    ddc *= digits
                    rvstart *= digits
                    rvend *= digits
                    mind = int(minv * ddc)
                    maxd = int(maxv * ddc)
                    rv = rv * digits + random.randint(0, digits - 1)
                    rvs = rv + rvstart

指数采样器：扩展

上面的代码支持有理数值的 λ 参数。它可以扩展到支持任何大于 0 的实数值 λ 参数，只要 λ 可以重写为一个或多个分量的和，每个分量的小数部分都可以通过伯努利工厂算法模拟，该算法以等于该小数部分的概率输出正面。[^25]

更具体地说

1. 将 λ 分解为 n > 0 个正分量，这些分量之和为 λ。例如，如果 λ = 3.5，它可以分解为一个分量，即 3.5（其小数部分很容易模拟），如果 λ = π，它可以分解为四个分量，所有分量都是 (π / 4)，其模拟过程在“伯努利工厂算法”中描述，并非微不足道。
2. 对于以这种方式找到的每个分量 LC[i]，令 LI[i] 为 floor(LC[i])，令 LF[i] 为 LC[i] − floor(LC[i])（LC[i] 的小数部分）。

然后可以修改上面的代码，如下所示：

• exprandnew 被修改，使其不再接受 lamdanum 和 lamdaden，而是接受上面描述的分量列表。每个分量都存储为 LI[i] 和一个模拟 LF[i] 的算法。

• zero_or_one_exp_minus(a, b) 被“伯努利工厂算法”中描述的算法 exp(− z) 替换，其中 z 是实数值的 λ 参数。

• logisticexp(a, b, index+1) 被“伯努利工厂算法”中描述的算法 1 / (1 + exp(z / 2^{index + 1})) (LogisticExp) 替换，其中 z 是实数值的 λ 参数。

正确性测试

 

Beta采样器

为了测试本文中 Beta 采样器的正确性，使用 SciPy 的 kstest 方法对各种 alpha 和 beta 值以及默认精度 53 应用了 Kolmogorov-Smirnov 检验。测试代码非常简单：kst = scipy.stats.kstest(ksample, lambda x: scipy.stats.beta.cdf(x, alpha, beta))，其中 ksample 是使用上述采样器生成的随机变量样本。请注意，SciPy 默认使用双边 Kolmogorov-Smirnov 检验。

请参阅正确性测试的结果。对于每对参数，抽取了五个样本，每个样本包含 50,000 个数字，结果显示了这五个样本获得的最低和最高 Kolmogorov-Smirnov 统计量和 p 值。请注意，p 值极接近 0 或 1 强烈表明样本并非来自相应的 Beta 分布。

ExpRandFill

为了测试 exprandfill 方法（它实现了 ExpRandFill 算法）的正确性，使用 SciPy 的 kstest 方法对各种 λ 值和默认精度 53 应用了 Kolmogorov-Smirnov 检验。测试代码非常简单：kst = scipy.stats.kstest(ksample, lambda x: scipy.stats.expon.cdf(x, scale=1/lamda))，其中 ksample 是使用上述 exprand 方法生成的随机变量样本。请注意，SciPy 默认使用双边 Kolmogorov-Smirnov 检验。

下表显示了正确性测试的结果。对于每个参数，抽取了五个样本，每个样本包含 50,000 个数字，结果显示了这五个样本获得的最低和最高 Kolmogorov-Smirnov 统计量和 p 值。请注意，p 值极接近 0 或 1 强烈表明样本并非来自相应的指数分布。

λ	统计量	p 值
1/10	0.00233-0.00435	0.29954-0.94867
1/4	0.00254-0.00738	0.00864-0.90282
1/2	0.00195-0.00521	0.13238-0.99139
2/3	0.00295-0.00457	0.24659-0.77715
3/4	0.00190-0.00636	0.03514-0.99381
9/10	0.00226-0.00474	0.21032-0.96029
1	0.00267-0.00601	0.05389-0.86676
2	0.00293-0.00684	0.01870-0.78310
3	0.00284-0.00675	0.02091-0.81589
5	0.00256-0.00546	0.10130-0.89935
10	0.00279-0.00528	0.12358-0.82974

ExpRandLess

为了测试 exprandless 的正确性，对涉及 e-rands 的分数和涉及 Python random.expovariate 方法的分数应用了双独立样本 T 检验。具体来说，分数计算为指数数与另一个指数数比较为小于的次数；对于相同的 λ，这个事件应该与它比较为大于的事件一样可能。（事实上，对于任何一对具有相同非退化分布的独立随机变量都应该如此；请参阅我的关于随机性提取的笔记中的命题 2。）下表后面的 Python 代码计算了 e-rands 和 expovariate 的这个分数。即使在这里，测试代码也非常简单：kst = scipy.stats.ttest_ind(exppyscores, exprandscores)，其中 exppyscores 和 exprandscores 分别是 exppyscore 或 exprandscore 的 20 个结果列表，并且 exppyscores 和 exprandscores 中包含的结果是相互独立生成的。

下表显示了正确性测试的结果。对于每对参数，结果显示了 20 个结果中获得的最低和最高 T 检验统计量和 p 值。请注意，p 值极接近 0 或 1 强烈表明指数随机变量未以预期概率比较为小于或大于。

左 λ	右 λ	统计量	p 值
1/10	1/10	-1.21015 – 0.93682	0.23369 – 0.75610
1/10	1/2	-1.25248 – 3.56291	0.00101 – 0.39963
1/10	1	-0.76586 – 1.07628	0.28859 – 0.94709
1/10	2	-1.80624 – 1.58347	0.07881 – 0.90802
1/10	5	-0.16197 – 1.78700	0.08192 – 0.87219
1/2	1/10	-1.46973 – 1.40308	0.14987 – 0.74549
1/2	1/2	-0.79555 – 1.21538	0.23172 – 0.93613
1/2	1	-0.90496 – 0.11113	0.37119 – 0.91210
1/2	2	-1.32157 – -0.07066	0.19421 – 0.94404
1/2	5	-0.55135 – 1.85604	0.07122 – 0.76994
1	1/10	-1.27023 – 0.73501	0.21173 – 0.87314
1	1/2	-2.33246 – 0.66827	0.02507 – 0.58741
1	1	-1.24446 – 0.84555	0.22095 – 0.90587
1	2	-1.13643 – 0.84148	0.26289 – 0.95717
1	5	-0.70037 – 1.46778	0.15039 – 0.86996
2	1/10	-0.77675 – 1.15350	0.25591 – 0.97870
2	1/2	-0.23122 – 1.20764	0.23465 – 0.91855
2	1	-0.92273 – -0.05904	0.36197 – 0.95323
2	2	-1.88150 – 0.64096	0.06758 – 0.73056
2	5	-0.08315 – 1.01951	0.31441 – 0.93417
5	1/10	-0.60921 – 1.54606	0.13038 – 0.91563
5	1/2	-1.30038 – 1.43602	0.15918 – 0.86349
5	1	-1.22803 – 1.35380	0.18380 – 0.64158
5	2	-1.83124 – 1.40222	0.07491 – 0.66075
5	5	-0.97110 – 2.00904	0.05168 – 0.74398

def exppyscore(ln,ld,ln2,ld2):
        return sum(1 if random.expovariate(ln*1.0/ld)<random.expovariate(ln2*1.0/ld2) \
              else 0 for i in range(1000))

def exprandscore(ln,ld,ln2,ld2):
        return sum(1 if exprandless(exprandnew(ln,ld), exprandnew(ln2,ld2)) \
              else 0 for i in range(1000))

在0到1之间支持的连续分布的精确模拟

本文档中的 Beta 采样器展示了一种通过伯努利工厂模拟在域 [0, 1] 上具有 PDF 的广泛连续分布的通用方法。这种通用方法可以使用以下算法采样遵循这些分布之一的数字。该算法允许任何任意基数 b（例如二进制的 2）。（另请参阅 (Devroye 1986, ch. 2, sec. 3.8, exercise 14)[^23]。）

1. 创建一个具有正号、整数部分为 0 且小数部分为空的均匀 PSRN。创建一个使用该 PSRN 的 SampleGeometricBag 伯努利工厂。

2. 随着 PSRN 构建一个均匀随机变量，以可以通过伯努利工厂算法表示的概率接受 PSRN（该算法将步骤 1 中的 SampleGeometricBag 工厂作为其输入的一部分），否则拒绝它。（这些算法中的一些可以在“伯努利工厂算法”中找到。）令 f(U) 为该伯努利工厂建模的概率密度函数 (PDF)，其中 U 是由 PSRN 构建的均匀随机变量。f 的域等于开区间 (0, 1) 或该区间的一个子集，并在其域中的任何地方返回 [0, 1] 的值。f 是底层连续分布的 PDF，或 PDF 乘以一个（可能未知）常数因子。然而，正如 Keane 和 O'Brien [^6] 所示，此步骤仅在以下情况下才有效——

   • f(λ) 在其域上是常数，或者
   • f(λ) 在其域上是连续且多项式有界的（多项式有界意味着 f(λ) 和 1−f(λ) 都由 min(λⁿ, (1−λ)ⁿ) 从下方有界，其中 n 为某个整数），

   他们表明，域为 [0, 1/2) 的 2 * λ 是一个不允许伯努利工厂的函数。请注意，概率可以是常数，包括无理数；请参阅“无理常数算法”以了解如何模拟常数概率。

3. 如果 PSRN 被接受，可选地用均匀随机数字填充 PSRN 的小数部分，以使其具有 n 位小数部分（类似于上面的 FillGeometricBag），其中 n 是一个精度参数，然后返回 PSRN。

然而，此算法的速度关键取决于 f 在 [0, 1] 中的众数（最高点）。[^26] 随着众数趋近于 0，平均拒绝率增加。实际上，此步骤在空间中 1×1 的区域内均匀随机生成一个点。如果众数接近 0，f 将只覆盖该区域的一小部分，因此生成的点落到 f 区域之外而不得不被拒绝的可能性很高。

Beta 分布在 (1) 处的 PDF 符合 Keane 和 O'Brien 的要求（对于 alpha 和 beta 都大于 1），因此可以通过伯努利工厂模拟，并涵盖在此通用算法中。

该算法可以修改为在区间 [m, m + y]（其中 m 和 y 是有理数且 y 大于 0）而不是 [0, 1] 中生成随机变量，如下所示：

1. 应用上述算法，但使用修改后的函数 f′(x) = f(x * y + m) 而不是 f，其中 x 是由 PSRN 构建的在 [0, 1] 中的数字，并且不选择按该算法步骤 3 中给出的方式填充 PSRN。
2. 通过“乘法”中的第二个算法将结果随机 PSRN 乘以 y。（请注意，如果 y 的形式为 bⁱ，此步骤相对简单。）
3. 通过“加法和减法”中的第二个算法将 m 添加到结果随机 PSRN 中。

请注意，此处，函数 f′ 必须满足 Keane 和 O'Brien 的要求。（例如，取函数 sqrt((x - 4) / 2)，它不是伯努利工厂函数。如果我们现在寻求从区间 [4, 4+2¹] = [4, 6] 中采样，则步骤 2 中使用的 f 现在是 sqrt(x)，它是伯努利工厂函数，因此我们可以应用此算法。）

一个例子：连续伯努利分布

连续伯努利分布 (Loaiza-Ganem and Cunningham 2019)[^27] 旨在显著提高变分自编码器（一种机器学习模型）在对取值在区间 [0, 1] 的连续数据（包括“近似二进制”图像数据）进行建模时的性能。

连续伯努利分布接受一个参数 lamda（一个在 [0, 1] 中的数字），并在区间 [0, 1] 中取值，其概率与以下值成比例——

pow(lamda, x) * pow(1 - lamda, 1 - x).

同样，此函数满足 Keane 和 O'Brien 提出的要求，因此可以通过伯努利工厂进行模拟。因此，可以按照下文所述在 Python 中模拟此分布。

连续伯努利分布的采样算法如下。它使用一个输入硬币，该硬币以概率 lamda 返回 1。

1. 创建一个正号零整数部分均匀 PSRN。
2. 创建一个互补 lambda 伯努利工厂，它返回输入硬币结果的 1 减去该结果。
3. 从均匀 PSRN 的小数部分中删除所有数字。这将导致在以下步骤中，区间 [0, 1] 中有一个“空”的均匀随机变量 U，该步骤将以概率 lamda^U*(1−lamda)^1−U)（Beta 分布的比例概率）接受 U，因为 U 正在构建中。
4. 调用“伯努利工厂算法”中描述的算法 λ^μ，使用输入硬币作为 λ 硬币，并使用 SampleGeometricBag 作为 μ 硬币（这将以概率 lamda^U 返回 1）。如果结果为 0，则转到步骤 3。
5. 调用算法 λ^μ，使用互补 lambda 伯努利工厂作为 λ 硬币，并使用 SampleGeometricBagComplement 算法作为 μ 硬币（这将以概率 (1-lamda)^1−U 返回 1）。如果结果为 0，则转到步骤 3。（请注意，步骤 4 和 5 相互独立，可以按任何顺序执行，而不会影响正确性。）
6. U 已被接受，因此返回 FillGeometricBag 的结果。

以下是采样连续伯努利分布的 Python 代码。

def _twofacpower(b, fbase, fexponent):
    """ Bernoulli factory B(p, q) => B(p^q).
           - fbase, fexponent: Functions that return 1 if heads and 0 if tails.
             The first is the base, the second is the exponent.
             """
    i = 1
    while True:
        if fbase() == 1:
            return 1
        if fexponent() == 1 and \
            b.zero_or_one(1, i) == 1:
            return 0
        i = i + 1

def contbernoullidist(b, lamda, precision=53):
    # Continuous Bernoulli distribution
    bag=[]
    lamda=Fraction(lamda)
    gb=lambda: b.geometric_bag(bag)
    # Complement of "geometric bag"
    gbcomp=lambda: b.geometric_bag(bag)^1
    fcoin=b.coin(lamda)
    lamdab=lambda: fcoin()
    # Complement of "lambda coin"
    lamdabcomp=lambda: fcoin()^1
    acc=0
    while True:
       # Create a uniform random variate (U) bit-by-bit, and
       # accept it with probability lamda^U*(1-lamda)^(1-U), which
       # is the unnormalized PDF of the beta distribution
       bag.clear()
       # Produce 1 with probability lamda^U
       r=_twofacpower(b, lamdab, gb)
       # Produce 1 with probability (1-lamda)^(1-U)
       if r==1: r=_twofacpower(b, lamdabcomp, gbcomp)
       if r == 1:
             # Accepted, so fill up the "bag" and return the
             # uniform number
             ret=_fill_geometric_bag(b, bag, precision)
             return ret
       acc+=1

复杂性

生成随机变量的算法的位复杂度以该算法平均使用的无偏随机位数来衡量。

一般原则

现有工作展示了如何计算任何概率分布的位复杂度

• 对于具有概率密度函数 (PDF) 的一维分布，位复杂度由 DE + prec - 1 个随机位下界，其中 DE 是分布的微分熵，prec 是随机变量小数部分的位数 (Devroye and Gravel 2020)[^3]。
• 对于离散分布（具有单独发生概率的随机整数分布），位复杂度由所涉及所有概率的二进制熵之和下界 (Knuth and Yao 1976)[^28]。（对于给定概率 p，二进制熵是 0 - p*log2(p)，其中 log2(x) = ln(x)/ln(2)。）最优算法平均将在此下界内达到 2 位。

例如，在指数分布的情况下，DE 是 log2(exp(1)/λ)，因此此分布的最小位复杂度是 log2(exp(1)/λ) + prec − 1，因此如果 prec = 20，当 λ = 1 时，此最小值约为 20.443 位，当 λ 增大时减小，当 λ 减小时增大。在任何其他具有 PDF 的分布的情况下，DE 是 f(x) * log2(1/f(x)) 在所有有效值 x 上的积分，其中 f 是分布的 PDF。

尽管现有工作给出了算法平均所需的随机位数下限，但大多数算法在实践中通常不会达到这些下限。

通常，如果一个算法调用其他生成随机变量的算法，则总预期位复杂度为——

• 对这些其他算法的每次调用的预期次数，乘以
• 每次此类调用的位复杂度。

特定算法的复杂性

此处给出的 Beta 采样器和指数采样器平均将比位复杂度的下限使用更多的位，特别是由于它们一次生成一个数字的 PSRN。

由于 zero_or_one 方法作为涉及随机位的伯努利试验的性质，它平均通常使用 2 个随机位，另请参阅 (Lumbroso 2013, Appendix B)[^29]。但是，如果它的两个参数相同，它不使用任何随机位。

对于基数为 2 的 SampleGeometricBag，位复杂度有两个组成部分。

• 一个组件来自从公平硬币中采样正面数量直到第一次反面，如下所示：
  • 最优下界：由于随机变量的二进制熵为 2，因此最优下界为 2 位。
  • 最优上界：4 位。
• 另一个组件来自用随机位填充部分采样随机数的小数部分。这里的复杂度取决于对同一 PSRN 调用 SampleGeometricBag 的次数，称其为 n。那么预期的位数是 n 次调用后以此方式填充的位位置的预期数量。

SampleGeometricBagComplement 具有与 SampleGeometricBag 相同的位复杂度。

FillGeometricBag 的位复杂度很容易找到。对于基数 2，它只使用一位来采样位置小于 p 的每个未填充数字。（对于基数不为 2 的情况，以这种方式采样每个数字可能不是最优的，因为数字是逐个生成的，并且随机位不会在几个数字之间重复使用。）因此，对于使用 SampleGeometricBag 和 FillGeometricBag 且 p 位的算法，这两个算法对复杂度的贡献平均在 p + g * 2 到 p + g * 4 位之间，其中 g 是对 SampleGeometricBag 的调用次数。（如果 FillGeometricBag 使用除简单截断之外的舍入机制实现，则此复杂度可能会增加 1 位。）

在加权水库采样中的应用

加权水库采样（从未知大小的列表中随机选择一个项目）通常通过以下方式实现：

• 为遇到的每个项目分配一个权重（一个大于或等于 0 的整数），称其为 w，
• 为每个项目提供一个指数随机变量，其 λ = w，称其为键，以及
• 选择具有最小键的项目

（另请参阅 (Efraimidis 2015)[^30]）。然而，使用完全采样的指数随机变量作为键（例如，在常见浮点算术中，朴素的习惯用法 -ln(1-X)/w，其中 X 是区间 [0, 1] 中的均匀随机变量）可能导致不精确采样，因为键的精度有限，可能存在多个项目具有相同随机键的情况（这可能导致这些项目的采样取决于它们的顺序而不是随机性），并且最大权重未知。本文中给出的部分采样 e-rands 消除了不精确采样的问题。这尤其因为 exprandless 方法只返回两个答案之一——“小于”或“大于”——并根据需要从两个 e-rands 中采样，以便它们在操作结束时彼此不同。（这不是问题，因为随机生成的实数预计以概率 1 彼此不同。）另一个原因是部分采样 e-rands 具有潜在的任意精度。

致谢

我感谢 Claude Gravel 审阅了本文的早期版本。

其他文档

以下是我撰写的关于随机数和伪随机数生成主题的一些其他文章。所有这些都是开源的。

• 随机数生成器应用建议
• 随机化和采样方法
• 更多随机采样方法
• 离散分布代码生成器
• 涉及随机化最常见的主题
• Bernoulli 工厂算法
• 更多任意精度采样算法
• 高质量 PRNG 的测试
• 高质量 PRNG 的示例

注释

[^1]: Karney, C.F.F., 2016. 从正态分布精确采样。ACM Transactions on Mathematical Software (TOMS), 42(1), pp.1-14。另见：“从正态分布精确采样”，arXiv:1303.6257v2 [physics.comp-ph], 2014。

[^2]: Philippe Flajolet, Nasser Saheb。生成指数分布变量的复杂性。[研究报告] RR-0159, INRIA. 1982. inria-00076400。

[^3]: Devroye, L., Gravel, C.，“仅使用有限数量的无偏、独立且同分布的随机比特进行随机变量生成”，arXiv:1502.02539v6 [cs.IT]，2020。

[^4]: Thomas, D.B. 和 Luk, W.，2008年9月。使用独立伯努利变量从指数分布中抽样。收录于 2008 年国际现场可编程逻辑和应用会议（第 239-244 页）。IEEE。

[^5]: A. Habibizad Navin, R. Olfatkhah 和 M. K. Mirnia，“指数随机变量的数据导向模型”，2010 年第二届先进计算机控制国际会议，沈阳，2010 年，第 603-607 页，doi: 10.1109/ICACC.2010.5487128。

[^6]: Keane, M. S., 和 O'Brien, G. L.，“伯努利工厂”，ACM Transactions on Modeling and Computer Simulation 4(2), 1994。

[^7]: Flajolet, P., Pelletier, M., Soria, M.，“关于布丰机器和数字”，arXiv:0906.5560v2 [math.PR]，2010。

[^8]: Pedersen, K.，“重新调整您的分位数函数”，arXiv:1704.07949 [stat.CO]，2018。

[^9]: von Neumann, J.，“与随机数字相关的各种技术”，1951。

[^10]: 正如冯·诺依曼 (von Neumann) (1951) 所指出的，一个介于 0 和 1 之间的均匀随机变量可以通过“并置足够的随机二进制数字”来生成。从这个意义上讲，该变量是 X1/B¹ + X2/B² + ...（其中 B 是数字基数 2，X1、X2 等是在区间 [0, B) 中独立的均匀随机整数），或许“强制最后一个 [随机比特] 为 1”以“避免任何偏差”。不难看出，这种方法可以应用于生成任何基数的任何数字展开，而不仅仅是 2。

[^11]: Yusong Du, Baoying Fan, 和 Baodian Wei，“一种改进的标准正态分布精确抽样算法”，arXiv:2008.03855 [cs.DS]，2020。

[^12]: 这意味着分布域的每个零体积（勒贝格测度零）子集（例如有限点集）的概率为零。等价地，这意味着分布具有概率密度函数。

[^13]: Oberhoff, Sebastian，“精确抽样和前缀分布”，论文与学位论文，威斯康星大学密尔沃基分校，2018。

[^14]: Hill, T.P. 和 Schürger, K.，2005 年。随机变量的数字和有效数字的规律性。Stochastic processes and their applications, 115(10), pp.1723-1743。

[^15]: J.F. Williamson，“弯曲表面上点的随机选择”，Physics in Medicine & Biology 32(10)，1987。

[^16]: Boehm, Hans-J. “迈向实数 API。”收录于第 41 届 ACM SIGPLAN 编程语言设计与实现会议论文集，第 562-576 页。2020。

[^17]: Hans-J. Boehm. 1987. 数值程序的构造性实数解释。收录于 SIGPLAN ’87 解释器和解释技术研讨会论文集。214-221。

[^18]: Goubault-Larrecq, Jean, Xiaodong Jia, 和 Clément Théron。“统计编程语言的领域理论方法。”arXiv 预印本 arXiv:2106.16190 (2021)（特别是第 12.3 节）。

[^19]: Brassard, G., Devroye, L., Gravel, C.，“远程采样及其在通用纠缠模拟中的应用”，Entropy 2019(21)(92)，doi:10.3390/e21010092。

[^20]: A. Habibizad Navin, Fesharaki, M.N., Teshnelab, M. 和 Mirnia, M.，2007 年。“均匀随机变量的数据导向建模：应用方法”。World Academy Science Engineering Technology, 21, pp.382-385。

[^21]: Nezhad, R.F., Effatparvar, M., Rahimzadeh, M., 2013 年。“设计一种通用数据导向随机数生成器”，International Journal of Modern Education and Computer Science 2013(2), pp. 19-24。

[^22]: 实际上，对于每个支持的整数 n，树给出获取 [n*p, (n+1)*p] 范围内值的概率，其中 p 是树的分辨率，例如 1/100000。

[^23]: Devroye, L.，非均匀随机变量生成，1986。

[^24]: Fan, Baoying 等。“通过早期拒绝生成指数分布变量。”2019 IEEE 第五届计算机与通信国际会议 (ICCC) (2019): 1307-1311。

[^25]: 实际上，多亏了 Flajolet 等人（2010）的“几何袋”技术，那个小数部分甚至可以来自均匀 PSRN。

[^26]: 更具体地说，是本质上确界，即函数在 [0, 1] 中的最高点，忽略零体积或测度为零的集合。然而，这里众数也是正确的，因为不连续的 PDF 不允许伯努利工厂，而这是步骤 2 所要求的。

[^27]: Loaiza-Ganem, G., Cunningham, J.P.，“连续伯努利：修复变分自动编码器中普遍存在的错误”，arXiv:1907.06845v5 [stat.ML]，2019。

[^28]: Knuth, Donald E. 和 Andrew Chi-Chih Yao.“非均匀随机数生成的复杂性”，收录于 Algorithms and Complexity: New Directions and Recent Results，1976。

[^29]: Lumbroso, J.，“从掷硬币生成最优离散均匀分布及其应用”，arXiv:1304.1916 [cs.DS]。

[^30]: Efraimidis, P. “数据流上的加权随机采样”，arXiv:1012.0256v2 [cs.DS]，2015。

[^31]: Glen, A.G., Leemis, L.M. 和 Drew, J.H.，2004 年。计算两个连续随机变量乘积的分布。计算统计与数据分析，44(3)，第 451-464 页。

[^32]: S. Kakutani，“论无限乘积测度的等价性”，Annals of Mathematics 1948。

[^33]: George Marsaglia. “具有独立二进制数字的随机变量。”Ann. Math. Statist. 42 (6) 1922 - 1929，1971 年 12 月。https://doi.org/10.1214/aoms/1177693058。

[^34]: Chatterji, S. D.。“某些诱导测度及其‘支持’的分数维数。”Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete 3 (1964): 184-192。

附录

SampleGeometricBag算法的等效性

对于 SampleGeometricBag，有两个版本：一个用于二进制（基数 2），另一个用于其他基数。以下是为什么这两个版本在二进制情况下是等效的。第一个算法的步骤 2 采样一个临时随机变量 N。这可以通过生成无偏随机位（即，每个位是 0 或 1，以相同概率选择）直到以这种方式生成一个零来实现。这里有三种相关情况。

• 生成的位是 1，这以 50% 的几率发生。这意味着位位置被跳过，算法继续到下一个位置。在算法 3 中，这对应于转到步骤 3，因为 a 的小数部分等于 b 的小数部分，这与小数部分不相等的情况相比也以 50% 的几率发生（因为 a 在算法过程中完全构建）。
• 生成的位是零，算法在位置 N 采样（或检索）一个零位，这以 25% 的几率发生。在算法 3 中，这对应于返回 0，因为 a 的小数部分小于 b 的小数部分，这以相同的概率发生。
• 生成的位是零，算法在位置 N 采样（或检索）一个一位，这以 25% 的几率发生。在算法 3 中，这对应于返回 1，因为 a 的小数部分大于 b 的小数部分，这以相同的概率发生。

UniformMultiply算法

以下算法 (UniformMultiply) 展示了如何将两个在其小数部分存储相同基数（基数）数字的均匀 PSRN (a 和 b) 相乘，并获得一个均匀 PSRN 作为结果。输入的 PSRN 可能有正负号，并且假定它们的整数部分和符号已被采样。

该算法目前仅在每个 PSRN 的小数部分至少有一个非零数字时有效；否则，它可能会产生遵循不正确分布的结果。这种情况可以通过应用“乘法”部分的注释来处理，但这会使下面的算法进一步复杂化。

1. 如果a在小数部分最后一个已采样数字之前有未采样数字，则将每个未采样数字设置为均匀随机选择的数字。对b也进行同样的操作。
2. 如果 a 的小数部分位数少于 b（反之亦然），则采样足够的位数（通过将它们设置为均匀随机数字，例如如果 a 和 b 存储二进制或基数 2 数字，则设置为无偏随机位），以便两个 PSRN 的小数部分具有相同的位数。
3. 如果 a 或 b 的整数部分为 0 且小数部分没有非零数字，则执行以下操作。
   1. 向 a 的小数部分追加一个均匀随机选择的数字。对 b 执行相同的操作。
   2. 如果 a 或 b 的整数部分为 0 且小数部分没有非零数字，请转到上一个子步骤。
4. 令 afp 为 a 的整数部分和小数部分打包成一个整数，如示例中所述；令 bfp 为 b 的整数部分和小数部分以相同方式打包。 (例如，如果 a 代表数字 83.12344...，则 afp 为 8312344。) 令 digitcount 为 a 的小数部分中的数字位数。
5. 计算 n1 = afp*bfp, n2 = afp*(bfp+1), n3 = (afp+1)*bfp, 和 n4 = (afp+1)*(bfp+1)。
6. 将 minv 设置为 n1，将 maxv 设置为 n2。将 midmin 设置为 min(n2, n3)，将 midmax 设置为 max(n2, n3)。
   • 数字 minv 和 maxv 是将区间乘法应用于 PSRN a 和 b 的结果的下限和上限。例如，如果 a 是 0.12344...，b 是 0.38925...，它们的小数部分相加形成 c = 0.51269...，或者区间 [0.51269, 0.51271]。然而，所得的 PSRN 在其区间内不是均匀分布的。在乘法的情况下，分布几乎是一个梯形，其域是区间 [minv, maxv]，其顶部由 midmin 和 midmax 界定。（参见本节末尾的注 1。）
7. 创建一个新的均匀 PSRN，ret。如果 a 的符号为负且 b 的符号为负，反之亦然，则将 ret 的符号设置为负。否则，将 ret 的符号设置为正。
8. 将 z 设置为区间 [0, maxv−minv) 中的均匀随机整数。
9. 如果 z < midmin−minv 或者 z ≥ midmax − minv，我们将分别从“梯形”的左侧或右侧进行采样。在这种情况下，执行以下操作：
   1. 将 x 设置为 minv + z。创建 psrn，一个具有正号和空小数部分的 PSRN。
   2. 如果 z < midmin − minv (左侧)，将 psrn 的整数部分设置为 x − minv，然后运行稍后给出的子算法 1，参数为 minv 和 psrn。(子算法以 ln((minv+psrn)/minv) 的概率返回 1。)
   3. 如果 z ≥ midmin − minv (右侧)，将 psrn 的整数部分设置为 x − midmax，然后运行稍后给出的子算法 2，参数为 maxv、midmax 和 psrn。(子算法以 ln(maxv/(midmax+psrn)) 的概率返回 1。)
   4. 如果子算法 1 或 2 返回 1，则算法成功，因此执行以下操作：
      1. 将 s 设置为 ru。
      2. 将 s 的 n*2 个最低有效数字转移到 ret 的小数部分，其中 n 是 a 的小数部分中的数字位数。（请注意，ret 的小数部分从最高位到最低位存储数字。）
      3. 将 psrn 小数部分中的数字追加到 ret 小数部分的末尾。
      4. 将 ret 的整数部分设置为 floor(s/base^n*2)。 (例如，如果 base 是 10，n*2 是 4，s 是 342978，那么 ret 的小数部分设置为 [2, 9, 7, 8]，ret 的整数部分设置为 34。) 最后，返回 ret。
   5. 如果子算法 1 或 2 返回 0，则中止这些子步骤并转到步骤 8。
10. （如果到达这里，我们已经到达梯形的中间部分，它是平坦且均匀的。）如果 n2 > n3，则运行稍后给出的子算法 3，参数为 afp（以 ln(1+1/afp) 的概率返回 1）。否则，运行子算法 3，参数为 bfp（以 ln(1+1/bfp) 的概率返回 1）。在任何一种情况下，如果子算法返回 0，则转到步骤 8。
11. （算法成功。）将 s 设置为 minv + z，然后将 s 的 (n*2) 个最低有效数字转移到 ret 的小数部分，然后将 ret 的整数部分设置为 floor(s/base^n*2)，然后返回 ret。

以下子算法由 UniformMultiply 使用。它们都涉及相同的底层函数 ln(1+x)，其算法在“伯努利工厂算法”页面中提及。

• 子算法 ln(1+x) 接受一个输入算法，并以 ln(1+x) 的概率返回 1，其中 x 是输入算法返回 1 的概率。
  • 重复执行以下过程，直到此子算法返回一个值：
    1. 生成一个无偏随机比特。如果该比特是 1（以 1/2 的概率发生），则运行输入算法并返回结果。
    2. 如果 u 尚未创建，则创建 u，一个带有正号、整数部分为 0 且小数部分为空的均匀 PSRN。
    3. 在 u 的小数部分上运行 SampleGeometricBag 算法，然后运行输入算法。如果调用和运行都返回 1，则返回 0。
• 子算法 1 接受两个参数 (minv 和 psrn)，并以 ln((minv+psrn)/minv) 的概率返回 1。运行 ln(1+x) 子算法，其输入算法如下：
  1. 设 p 为 psrn 的整数部分。在 [0, minv) 范围内均匀随机生成一个整数，称之为 i。
  2. 如果 i < p，返回 1。如果 i = p，翻转输入硬币并返回结果。如果都不是，返回 0。
• 子算法 2 接受三个参数 (maxv, midmax 和 psrn)，并以 ln(maxv/(midmax+psrn)) 的概率返回 1。运行 ln(1+x) 子算法，其输入算法如下：
  1. 设 p 为 psrn 的整数部分。将 d 设置为 maxv − p − midmax − 1，将 c 设置为 p + midmax。
  2. 以 c / (1 + c) 的概率，执行以下操作：
     • 在 [0, c) 范围内均匀随机生成一个整数，称之为 i。如果 i < d，返回 1。如果 i = d，在 psrn 的小数部分上运行 SampleGeometricBag 并返回 1 减去结果。如果 i > d，返回 0。
  3. 在 psrn 的小数部分上运行 SampleGeometricBag。如果结果是 1，返回 0。否则，转到步骤 2。
• 子算法 3 接受一个参数（这里称为 n），并以 ln(1+1/n) 的概率返回 1。运行 ln(1+x) 子算法，其输入算法如下：“返回一个数字，以 1/n 的概率为 1，否则为 0。”

注意：两个均匀 PSRN 的乘积分布并非精确的梯形，而是遵循一种不那么平凡的分布；当每个 PSRN 都远离 0 时，分布的左右两侧并非精确的“三角形”形状，而是基于对数函数。然而，当分布的“宽度”变小时，这些对数函数会趋近于三角形形状。参见 Glen 等人（2004）[^31] 和一个Stack Exchange 问题。

均匀的均匀生成均匀的乘积

本节包含证据表明，“乘法”一节注释中给出的算法正确地生成了两个均匀随机变量的乘积，一个在 [0, b] 中，另一个在 [c, d] 中，至少在 c = 0 时如此。

均匀(α, β) 随机变量的概率密度函数 (PDF) 在 x 位于 [α, β] 时为 1/(β−α)，否则为 0。此处称之为 UPDF(x, α, β)。

令 K = b*(d−c)。为了展示结果，我们找到两个如下所述的 PDF。

• 为了找到算法的 PDF，我们计算 UPDF(x, 0, Z+b*c) 的期望值，其中 Z 分布为均匀(0, K)。这通过计算积分（曲线下方的面积）实现，即 UPDF(x, 0, z+b*c)*UPDF(z, 0, K) 在区间 [0, K]（Z 可以取的值集）上关于 z 的积分。结果是 PDF1(x) = ln(b**2*c**2 - b**2*c*d + (b*c - b*d)*min(b*(-c + d), max(0, -b*c + x)))/(b*c - b*d) - ln(b**2*c**2 - b**2*c*d + b*(-c + d)*(b*c - b*d))/(b*c - b*d)。
• 第二个 PDF 是两个均匀随机变量乘积的 PDF，一个在 [0, b] 中，另一个在 [c, d] 中。根据 Rohatgi 的公式（另见 (Glen 等人 2004)[^31]），可以通过计算 UPDF(z, 0, b)*UPDF(x/z, c, d)/z 在区间 [0, ∞) 上关于 z 的积分来找到（注意这里 z 从不为负）。结果是 PDF2(x) = (ln(max(c,x/b)) - ln(max(c,d,x/b)))/(b*c-b*d)。

现在必须证明 PDF1 和 PDF2 在 x 位于区间 (0, b*d) 内时相等。将一个 PDF 减去另一个并简化，可以看到：

• 这两个 PDF 至少在 c = 0 时相等（以及当 b、d 和 x 都大于 0 时），并且到目前为止在所有计算中，当 b、c 和 d 被特定值替换时，它们都相等。
• PDFs 之间的简化差值的积分为 0，这强烈表明 PDFs 相等（这并非结论性，因为简化差值可能为负）。

Oberhoff的“精确拒绝采样”方法

以下描述了 Oberhoff 提出的一种算法，用于对取值范围在 [0, 1] 的连续分布进行采样，只要该分布具有概率密度函数 (PDF) 且 PDF “几乎处处”连续并有上界（Oberhoff 2018，第 3 节）[^13]，另见 (Devroye and Gravel 2020)[^3]。（请注意，如果 PDF 的域比 [0, 1] 宽，则函数需要分成一单位长的片段，以其面积成比例的概率随机选择一个片段，并将该片段移动，使其位于 [0, 1] 而不是其通常位置；参见 Oberhoff 第 11-12 页。）

1. 将 pdfmax 设置为 PDF（或 PDF 乘以可能未知的常数因子）在 [0, 1] 域上的上限。设 base 为返回值的数字的基数或底数（例如二进制为 2，十进制为 10）。
2. 将 prefix 设置为 0，将 prefixLength 设置为 0。
3. 将 y 设置为区间 [0, pdfmax] 中的均匀随机变量。
4. 令 pw 为 base^{−prefixLength}。将 lower 和 upper 分别设置为 PDF（或 PDF 乘以可能未知的常数因子）在域 [prefix * pw, prefix * pw + pw] 上的值的下限或上限。
5. 如果 y 结果大于 upper，则前缀被拒绝，所以转到步骤 2。
6. 如果 y 结果小于 lower，则前缀被接受。现在执行以下操作：
   1. 当 prefixLength 小于所需精度时，将 prefix 设置为 prefix * base + r，其中 r 是一个均匀随机数字，然后将 prefixLength 加 1。
   2. 返回 prefix * base^{−prefixLength}。（如果 prefixLength 以某种方式大于所需精度，则算法可以选择将返回值四舍五入到小数部分具有所需位数的数字，并选择舍入模式。）
7. 将 prefix 设置为 prefix * base + r，其中 r 是一个均匀随机数字，然后将 prefixLength 加 1，然后转到步骤 4。

由于此算法需要评估 PDF（或 PDF 的常数倍）并找到其在某个区间内的最大值和最小值（这通常需要浮点运算且通常不简单），因此此算法出现在附录而非正文中。此外，除非 y 是均匀 PSRN，否则从生成具有固定位数 y 中会产生额外的近似误差（另请参见“加权水库采样的应用”）。出于实际目的，步骤 4 中计算的下限和上限应取决于 prefixLength（prefixLength 越高，精度越高）。

Oberhoff 还描述了前缀分布，它们以与盒子面积成比例的概率对覆盖 PDF 的盒子进行采样，但这些分布必须支持固定的最大前缀长度，因此只能近似底层分布。

按数字概率设置数字

原则上，部分采样随机数可以通过找到一系列数字概率并根据这些概率设置该数字的数字来实现。然而，均匀分布和指数分布是此类仅有的实用分布。详细信息如下。

令 X 是 0.bbbbbbb... 形式的随机变量，其中每个 b 是一个独立的随机二进制数字（0 或 1）。

设 a_j 为位置 j 处的数字等于 1 的概率（小数点后第一个数字的 j = 1 开始）。

那么 Kakutani 定理 (Kakutani 1948)[^32] 指出，当且仅当 (a_j − 1/2) 的平方和收敛时，X 具有绝对连续[^12]分布。换句话说，当二进制数字离二进制点越来越远时，它们的偏差越来越小。另见 (Marsaglia 1971)[^33]，(Chatterji 1964)[^34]。

因此，如果我们能找到一个收敛于 1/2 的无限序列 a_j，并使用这些概率设置 X 的二进制数字，就可以构建这种绝对连续分布。然而，正如 Marsaglia (1971)[^33] 所展示的，绝对连续分布只能是以下之一：

1. 分布的概率密度函数 (PDF) 在 (0, 1) 中的每个开区间中的某处为零，而不在 [0, 1] 的所有地方为 0。因此，PDF 不连续。

2. PDF 在 1/2、1/4、1/8 等处为正，因此 PDF 在 (0, 1) 的所有地方都连续且为正，序列的形式是——

   a_j = exp(w/2^j)/(1 + exp(w/2^j))，

   其中 w 是一个常数。

3. 上述情况 2 中未描述 PDF，但它在 (0, 1) 的某个开区间上为正，因此 PDF 将是分段连续的，并且 X 可以乘以 2 的整数幂，使得新变量的分布具有情况 2 中描述的 PDF。

正如 Marsaglia 也表明的那样，当随机数字的基数不是 2（二进制）时，也有类似的结果。另见我的Stack Exchange 问题。

情况 2 有几个特例，包括：

• 均匀分布 (w = 0)。
• 速率为 1 的指数随机变量的小数部分 (w = −1；(Devroye and Gravel 2020)[^3])。
• 更一般地，速率为 λ 的指数变量的小数部分 (w = −λ)。
• 当 w > 0 时，1 减去速率为 w 的指数变量的小数部分。
• a_j = y^v/2j/(1 + y^v/2j)，其中 w = ln(y)*v，且 y > 0 和 v 为常数。