65.9K
CodeProject 正在变化。 阅读更多。
Home

Bernoulli 工厂算法

Peter Occil
starIconstarIconstarIconstarIconstarIcon

5.00/5 (12投票s)

2020年9月23日

公共领域

142分钟阅读

viewsIcon

33513

downloadIcon

201

将有偏“硬币抛掷”变成有偏“硬币抛掷”的算法，以及如何对其进行编码。

引言

给定一个正面概率为 λ 的未知硬币，从中抽取概率 f(λ)。换句话说，将一个朝一个方向偏斜（λ）的硬币变成朝另一个方向偏斜（f(λ)）的硬币。这就是Bernoulli 工厂问题。

本页面将对解决各种函数形式的 Bernoulli 工厂问题进行分类，这些算法被称为Bernoulli 工厂。

其中许多算法在 (Flajolet et al., 2010)(1) 中提出，但很多情况下没有提供分步说明。本页面提供了这些说明，以帮助程序员实现他们所描述的 Bernoulli 工厂。Python 模块 bernoulli.py 包含了几个 Bernoulli 工厂的实现。

本页面还包含仅使用随机比特来精确采样无理数概率的算法，这与 Bernoulli 工厂问题有关。（无理数是不能写成两个整数之比的数。）同样，其中许多算法也在 (Flajolet et al., 2010)(1) 中提出。

本页面侧重于精确采样概率的方法，不引入除输入中已存在的误差之外的舍入误差或其他误差（并假设我们有一个独立且无偏的随机比特源）。

补充说明可在以下位置找到： Bernoulli 工厂算法的补充说明

关于本文档

本文档为开源文档；如需最新版本，请参阅 源代码 或其在 GitHub 上的 渲染版本。您可以在 GitHub issues 页面 上就本文档发表评论。请参阅“请求和公开问题”以获取我正在寻求答案的本文档相关问题列表。

我鼓励读者实现本页面提供的任何算法，并报告他们的实现经验。这有助于改进本页面。

目录

关于 Bernoulli 工厂

Bernoulli 工厂（Keane 和 O'Brien 1994）(2) 是一种算法，它接收一个输入硬币（一种以未知概率返回 1（正面），否则返回 0（反面）的方法），并以一个取决于输入硬币正面概率的概率返回 0 或 1。未知的正面概率在此文档中称为λ。例如，Bernoulli 工厂算法可以接收一个正面概率为 λ 的硬币，并生成一个正面概率为 exp(−λ) 的硬币。

工厂函数是关联旧概率和新概率的已知函数。其定义域为闭区间 [0, 1] 或其子集，并返回 [0, 1] 中的概率。

许多 Bernoulli 工厂还可以访问外部随机性（例如无偏随机比特源）。如果一个工厂函数可以实现一个仅使用输入硬币且不使用外部随机性的 Bernoulli 工厂，则称该工厂函数是强可模拟的。

Keane 和 O'Brien (1994)(2) 表明，当且仅当—

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

以下显示了一些是工厂函数和不是工厂函数的函数。在下表中，ϵ 是大于 0 且小于 1/2 的数。

函数 f(λ)定义域f 能成为工厂函数吗？
0[0, 1]是；常数。
1[0, 1]是；常数。
1/2(0, 1)是；常数。
floor(λ/2)*3+1/4(0, 1)否；不连续。
2*λ[0, 1] 或 [0, 1/2)否；由于其图形在定义域的 (0, 1) 区间内某处触及 1，因此不是多项式有界的。(3)。
1−2*λ[0, 1] 或 [0, 1/2)否；由于其图形在定义域的 (0, 1) 区间内某处触及 0，因此不是多项式有界的。
2*λ[0, 1/2−ϵ]是；在定义域上连续且多项式有界（Keane 和 O'Brien 1994）(2)。
min(2 * λ, 1 − ϵ)[0, 1]是；在定义域上连续且多项式有界（Huber 2014，引言）(4)。
0 若 λ = 0，否则为 exp(−1/λ)(0, 1)否；不是多项式有界的，因为它比任何多项式都慢地远离 0。
ϵ 若 λ = 0，否则为 exp(−1/λ) + ϵ(0, 1)是；连续且远离 0 和 1。

如果 f 的定义域包含 0 和/或 1（因此输入硬币允许每次返回 1 或每次返回 0，分别），那么 f 只能是工厂函数，如果—

  1. f 在其定义域上是常数，或在定义域上是连续且多项式有界的，并且
  2. 当 0 在 f 的定义域中时，f(0) 等于 0 或 1，并且
  3. 当 1 在 f 的定义域中时，f(1) 等于 0 或 1，

除非可以使用外部随机性（除了输入硬币）。

算法

本节将展示一些工厂函数的算法，允许从输入硬币中采样不同类型的概率。

此处描述的算法并不总是带来最佳性能。实现可以修改这些算法，只要它们产生与此处描述的算法相同的结果即可。

这些算法假设除了输入硬币之外，还可以获得独立且无偏的随机比特源。但在许多情况下，它们可以仅使用这些硬币作为随机性来源来实现。有关详细信息，请参阅附录。

采样概率 f(λ) 的 Bernoulli 工厂算法充当 f(λ) 的无偏估计量。有关详细信息，请参阅附录。

实现说明

本节介绍适用于本文中算法的实现说明。应遵循这些说明，以避免在算法中引入错误。

在以下算法中

  • λ 是输入硬币的未知正面概率。
  • choose(n, k) = (1*2*3*...*n)/((1*...*k)*(1*...*(n−k))) = n!/(k! * (n − k)!) 是一个二项式系数，即从 n 个已标记项中选择 k 个项的方式数。例如，可以通过计算 Manolopoulos (2002)(5) 中对整数 i 在 [n−k+1, n] 范围内的 i/(n−i+1) 的结果，然后将这些结果相乘来计算。请注意，对于每个 m>0，choose(m, 0) = choose(m, m) = 1 且 choose(m, 1) = choose(m, m−1) = m；此外，在此文档中，当 k 小于 0 或大于 n 时，choose(n, k) 为 0。
  • n! = 1*2*3*...*n 也称为n 的阶乘。
  • “生成一个 uniform(0, 1) 随机变量”的指令可以实现为—
  • “生成指数随机变量”的指令可以实现为—
    • 通过创建一个空的指数 PSRN（最精确），或
    • 通过获取（在我关于 PSRN 的文章中描述的）**ExpRand** 或 **ExpRand2** 算法的结果，其速率为 1，或
    • 通过生成 -ln(1/X)，其中 X 是半开区间 (0, 1] 中的均匀随机变量（例如，“随机化和采样方法”中的 RNDRANGEMinExc(0, 1)）（精度较低）。

  • “按权重选择[整数]”的指令可以按以下方式之一实现：

    例如，“以与权重 [A, B, C] 成比例的概率选择 0、1 或 2”意味着以随机方式选择 0、1 或 2，以便 0 的选择概率为 A/(A+B+C)，1 的选择概率为 B/(A+B+C)，2 的选择概率为 C/(A+B+C)。

  • 采样一个数 u 意味着生成一个数，该数以概率 u 返回 1，否则返回 0。
    • 如果该数是一个均匀 PSRN，则调用 **SampleGeometricBag** 算法并传入 PSRN，然后取该调用的结果（将是 0 或 1）（最精确）。（**SampleGeometricBag** 在我关于 PSRN 的文章中有描述。）
    • 否则，可以通过生成另一个 uniform(0, 1) 随机变量 v 来实现，如果 v 小于 u 则生成 1，否则生成 0（精度较低）。
  • 当算法说“如果 a 小于 b”，其中 a 和 b 是随机变量时，意味着对这两个数字运行 **RandLess** 算法（如果它们都是 PSRN），或根据情况对 a 和 b 执行小于操作。（**RandLess** 在我关于 PSRN 的文章中有描述。）
  • 当算法说“如果 a 小于（或等于）b”，其中 a 和 b 是随机变量时，意味着对这两个数字运行 **RandLess** 算法（如果它们都是 PSRN），或根据情况对 a 和 b 执行小于等于操作。
  • 当算法中的一个步骤说“以概率 x”来引用一个可能发生也可能不发生的事件时，可以通过以下方式之一来实现：
    • 生成一个 uniform(0, 1) 随机变量 v（见上文）。如果 v 小于 x（见上文），则事件发生。
    • 将 x 转换为有理数 y/z，然后调用 ZeroOrOne(y, z)。如果调用返回 1，则事件发生。例如，如果指令是“以 3/5 的概率返回 1”，则实现为“调用 ZeroOrOne(3, 5)。如果调用返回 1，则返回 1。”ZeroOrOne 在我关于随机采样方法的文章中有描述。请注意，如果 x 不是有理数，则会产生舍入误差。
  • 为获得最佳结果，应使用精确的有理数算术（例如 Python 中的 Fraction 或 Ruby 中的 Rational）来实现这些算法。不鼓励使用浮点算术，因为它可能由于固定精度计算（如舍入和抵消）而引入错误。

通用函数 λ 的算法

本节介绍了用于采样多项式、有理函数或一般函数的概率的通用算法。

特定多项式

任何多项式都可以写成Bernstein 形式，即 ∑i = 0, ..., n choose(n, i) * λi * (1 − λ)n − i * a[i]，其中 n 是多项式的次数，a[i] 是其 n+1 个系数。

但是，唯一可以实现 Bernoulli 工厂的多项式是其系数（一旦多项式以 Bernstein 形式写出）全部在区间 [0, 1] 内的多项式，而这些多项式也是唯一可以用固定数量的硬币抛掷来模拟的函数（Goyal 和 Sigman 2012(6)；Qian et al. 2011）(7)；另见 Wästlund 1999，第 4 节(8)）。Goyal 和 Sigman 提供了一种模拟这些多项式的算法，该算法如下。

  1. 抛掷输入硬币 n 次，令 j 为此过程中硬币返回 1 的次数。
  2. 以 a[j] 的概率返回 1。否则，返回 0。

对于具有重复系数的某些多项式，以下是该算法的优化版本，并非由 Goyal 和 Sigman 提供。

  1. 将 j 设置为 0，将 i 设置为 0。如果 n 为 0，则返回 0。
  2. 如果 i 大于或等于 n，或者系数 a[k]，其中 k 在区间 [j, j+(n−i)] 内，全部相等，则返回一个以 a[j] 的概率为 1，否则为 0 的数。
  3. 抛掷输入硬币。如果返回 1，则将 j 加 1。
  4. 将 i 加 1 并转到步骤 2。

以下是另一个优化算法

  1. 将 j 设置为 0，将 i 设置为 0。如果 n 为 0，则返回 0。否则，生成一个 uniform(0, 1) 随机变量，称之为 u。
  2. 如果 u 小于最低系数的下界，则返回 1。否则，如果 u 小于（或等于）最高系数的上界，则转到下一步。否则，返回 0。
  3. 如果 i 大于或等于 n，或者系数 a[k]，其中 k 在区间 [j, j+(n−i)] 内，全部相等，则返回一个数，如果 u 小于 a[j]，则为 1，否则为 0。
  4. 抛掷输入硬币。如果返回 1，则将 j 加 1。
  5. 将 i 加 1 并转到步骤 3。

注释:

  1. 每个 a[i] 都作为一维贝塞尔曲线的控制点，其中 λ 是该曲线上的相对位置，曲线从 a[0] 开始，到 a[n] 结束。例如，给定控制点 0.2、0.3 和 0.6，当 λ = 0 时，曲线为 0.2，当 λ = 1 时，曲线为 0.6。（曲线在 λ = 1/2 时不一定在 0.3 处；通常，贝塞尔曲线除了第一个和最后一个外，不会穿过其控制点。）
  2. 模拟 Bernstein 形式的多项式问题与随机逻辑相关，它涉及模拟布尔函数（仅使用 AND、OR、NOT 和 XOR 操作的函数）产生的概率，这些函数以固定数量的比特作为输入，每个比特有单独的 1 而非 0 的概率，并输出单个比特（有关更多讨论，请参见 (Qian et al. 2011)(7)，Qian 和 Riedel 2008(9)）。
  3. 这些算法可以作为近似模拟任何连续函数 f（甚至任何将区间 [0, 1] 映射到 [0, 1] 的函数，即使它不连续）的方法。在这种情况下，a[j] 计算为 f(j/n)，从而使生成的多项式紧密逼近函数。事实上，如果 f 是连续的，可以选择足够大的 n 来达到给定的最大误差（这是所谓的“Weierstrass 近似定理”的结果），但在 f 不连续的情况下通常不是这样。有关更多信息，请参阅我的关于 Bernoulli 工厂的补充说明。

示例：考虑 (Thomas 和 Blanchet 2012)(10) 中讨论的以下抛物线函数：(1−4*(λ−1/2)2)*c，其中 c 在区间 (0, 1) 内。这是一个 2 次多项式，可以重写为 −4*c*λ2+4*c*λ，因此该幂形式的系数为 (0, 4*c, −4*c)，次数（n）为 2。通过将多项式重写为 Bernstein 形式（例如，通过 Ray 和 Nataraj (2012)(11) 的矩阵方法），我们得到系数 (0, 2*c, 0)。因此，对于此多项式，a[0] 为 0，a[1] 为 2*c，a[2] 为 0。因此，如果 c 在区间 (0, 1/2] 内，我们可以如下模拟此函数：“抛掷输入硬币两次。如果恰好一次抛掷返回 1，则返回一个以 2*c 的概率为 1，否则为 0 的数。否则，返回 0。”对于其他 c 值，该算法需要将多项式重写为 Bernstein 形式，然后多次提高重写多项式的次数以使其系数在 [0, 1] 内；随着 c 接近 1，所需次数趋于无穷大。(12)

 

多个硬币。 Niazadeh et al. (2021)(13) 描述了涉及一个或多个硬币的单项式，形式为 λ[1]a[1] * (1−λ[1])b[1]*λ[2]a[2] * (1−λ[2])b[2]* ... * λ[n]a[n] * (1−λ[n])b[n]，其中有 n 个硬币，λ[i] 是硬币 i 的正面概率，a[i] ≥ 0 和 b[i] ≥ 0 是硬币 i 的参数（特别是，在 a+b 次抛掷中，前 a 次必须是正面，其余必须是反面才能成功）。

  1. 对于 [1, n] 中的每个 i
    1. 抛掷 λ[i] 输入硬币 a[i] 次。如果任何一次抛掷返回 0，则返回 0。
    2. 抛掷 λ[i] 输入硬币 b[i] 次。如果任何一次抛掷返回 1，则返回 0。
  2. 返回 1。

同一篇论文还描述了这些单项式的加权和形式的多项式，即形式为 P = ∑j = 1, ..., k c[j]*M[j](λ) 的多项式，其中有 k 个单项式，M[j](.) 标识单项式 j，λ 标识硬币的正面概率，c[j] ≥ 0 是单项式 j 的权重。

令 C 为所有 c[j] 的总和。要模拟概率 P/C，请以与其权重成比例的概率选择一个单项式（参见“带权重选择（有放回）”），然后运行上述算法处理该单项式（另见“凸组合”，稍后）。

以下是一个特例

  • 如果只有一个硬币，多项式 P 是 Bernstein 形式的，当且仅当 c[j] = α[j]*choose(k−1, j−1)，其中 α[j] 是区间 [0, 1] 内的系数，并且对于每个单项式 j，a[1] = j−1 且 b[1] = k−j。

特定有理函数

有理函数是多项式的比。

根据 Mossel 和 Peres (2005)(14)，一个映射开区间 (0, 1) 到 (0, 1) 的函数可以通过有限状态机来模拟，当且仅当该函数可以写成系数为有理数的有理函数。

以下算法基于 Mossel 和 Peres 的论文以及 (Thomas 和 Blanchet 2012)(10)。它假设有理函数的形式为 D(λ)/E(λ)，其中—

  • D(λ) = ∑i = 0, ..., n λi * (1 − λ)n − i * d[i],
  • E(λ) = ∑i = 0, ..., n λi * (1 − λ)n − i * e[i],
  • 每个 d[i] 小于或等于对应的 e[i]，并且
  • 每个 d[i] 和每个 e[i] 是区间 [0, choose(n, i)] 内的整数或有理数，其中上界是具有 i 个 1 的 n 位字的总体数。

此处，d[i] 类似于具有 i 个 1 的 n 位字中“通过”的数量，而 e[i] 类似于该数量加上具有 i 个 1 的“失败” n 位字的数量。

算法如下。

  1. 抛掷输入硬币 n 次，令 j 为此过程中硬币返回 1 的次数。
  2. 以这些权重成比例地选择 0、1 或 2：[e[j] − d[j], d[j], choose(n, j) − e[j]]。如果以此方式选择了 0 或 1，则返回它。否则，转到步骤 1。

注释

  1. 在上面的公式中—

    • d[i] 可以替换为 δ[i] * choose(n,i)，其中 δ[i] 是区间 [0, 1] 内的有理数（从而表示给定字是具有 i 个 1 的所有 n 位字中“通过”字的概率），并且
    • e[i] 可以替换为 η[i] * choose(n,i)，其中 η[i] 是区间 [0, 1] 内的有理数（从而表示给定字是“通过”或“失败”字 among 所有具有 i 个 1 的 n 位字中的概率），

    然后 δ[i] 和 η[i] 可以被视为两个不同一维贝塞尔曲线的控制点，其中 δ 曲线始终位于 η 曲线的“上方”或“下方”。对于每条曲线，λ 是该曲线上的相对位置，曲线从 δ[0] 或 η[0] 开始，到 δ[n] 或 η[n] 结束。另请参阅下一节。

  2. 此算法可以修改为避免除输入硬币抛掷之外的额外随机性，方法是将硬币抛掷打包成一个 n 位字，并查找该字是否在所有具有 j 个 1 的 n 位字中属于某个特定的字类，但这样做可能不切实际（通常，大小为 2n 的查找表必须在设置步骤中先构建；随着 n 的增长，表的大小呈指数增长）。此外，这种方法仅在 d[i] 和 e[i] 是整数时才有效（或者如果 d[i] 被替换为 floor(d[i]) 且 e[i] 被替换为 ceil(e[i]) (Nacu 和 Peres 2005)(15)，但当然，这样做时会产生舍入误差）。另请参阅 (Thomas 和 Blanchet 2012)(10)。

  3. 与多项式一样，此算法（或稍后给出的算法）可以作为近似模拟任何工厂函数的方法，通过一个有理函数紧密逼近该函数。n 越高，近似越好，通常，n 次有理函数比 n 次多项式更好地逼近给定函数。但是，要用有理函数达到给定的误差容限，必须优化 n 以及 d[i] 和 e[i]。这与多项式情况不同，后者仅需优化次数 n。

示例：考虑函数 f(λ) = 1/(λ−2)2。这是一个有理函数，在这种情况下，它是由两个多项式的比值构成的，这两个多项式在区间 [0, 1] 上都非负。一种模拟 f 的算法是
(1) 抛掷输入硬币两次，令 heads 为此过程中硬币返回 1 的次数。
(2) 根据 heads，以与以下权重成比例的概率选择 0、1 或 2：heads=0 → [3, 1, 0]，heads=1 → [1, 1, 2]，heads=2 → [0, 1, 3]；如果以此方式选择了 0 或 1，则返回它；否则，转到步骤 1。
以下是为推导出此算法而准备 f 的方式
(1) 取分子 1，以及分母 (λ−2)2。将分母重写为 1*λ2 − 4*λ + 4。
(2) 将分子和分母重写为 2 次的齐次多项式（项具有相同次数的多项式）；请参阅“准备有理函数”中的“齐次化”部分。结果分别为 (1, 2, 1) 和 (4, 4, 1)。
(3) 将两个多项式都除以 4（最大系数），使它们都为 1 或更小。结果是 d = (1/4, 1/2, 1/4)，e = (1, 1, 1/4)
(4) 按照原始算法步骤 2 中的方式准备权重。结果是 [3/4, 1/4, 0]，[1/2, 1/2, 1]，以及 [0, 1/4, 3/4]，分别对应不同的 heads 计数。在此情况下，这些权重可以简化为整数而不影响算法：[3, 1, 0]，[1, 1, 2]，[0, 1, 3]。

“骰子企业”特例。 以下算法实现了 Morina et al. (2019)(16) 的“骰子企业”方法的一个特例。该算法以概率 PX(λ) / (P0(λ) + P1(λ) + ... + Pm−1(λ)) 返回 m 个结果之一（即 X，一个在 [0, m) 内的整数），其中 λ 是输入硬币的正面概率，m 大于或等于 2。具体来说，概率是一个有理函数，即多项式的比值。这里，所有 Pk(λ) 都具有以下形式的多项式：

  • 这些多项式是齐次的，即它们写成 ∑i = 0, ..., n λi * (1 − λ)n − i * a[i]，其中 n 是多项式的次数，a[i] 是系数。
  • 这些多项式具有相同的次数（即 n），并且所有 a[i] 都大于或等于 0。
  • 对于从 0 到 n 的每个 j，jth 系数的总和大于 0，除了总和列表可能以零开头和/或以零结尾。称此列表为 R。例如，如果 R 是 (2, 4, 4, 2) 或 (0, 2, 4, 0)，则此条件成立，但如果 R 是 (2, 0, 4, 3) 则不成立。

任何可以转化为上述形式的有理函数都可以通过附录中的“准备有理函数”进行处理。在此算法中，令 R[j] 为多项式jth 系数的总和（j 从 0 开始）。首先，定义以下操作

  • 给定 state、b、u 和 n，获取新状态:
    1. 如果 state > 0 且 b 为 0，则返回 state−1（如果 u 小于（或等于）PA）或 state（否则），其中 PA 是 R[state−1]/max(R[state], R[state−1])。
    2. 如果 state < n 且 b 为 1，则返回 state+1（如果 u 小于（或等于）PB）或 state（否则），其中 PB 是 R[state+1]/max(R[state], R[state+1])。
    3. 返回 state。

然后算法如下：

  1. 创建两个空列表：blist 和 ulist。
  2. 将 state1 设置为 R 中第一个非零项的位置。将 state2 设置为 R 中最后一个非零项的位置。在这两种情况下，位置都从 0 开始。如果 R 中的所有项都为零，则返回 0。
  3. 抛掷输入硬币并将结果（0 或 1）追加到 blist 的末尾。生成一个 uniform(0, 1) 随机变量并将其追加到 ulist 的末尾。
  4. (来自过去的单调耦合 (Morina et al., 2019)(16), (Propp 和 Wilson 1996)(17)。) 将 i 设置为 blist 的项数减 1，然后当 i 大于或等于 0 时
    1. 令 b 为 blist 中位置 i（从 0 开始）的项，令 u 为 ulist 中该位置的项。
    2. 给定 state1、b、u 和 n，获取新状态，并将 state1 设置为新状态。
    3. 给定 state2、b、u 和 n，获取新状态，并将 state2 设置为新状态。
    4. 将 i 减 1。
  5. 如果 state1 和 state2 不相等，则转到步骤 2。
  6. 令 b(j) 为多项式 j 的系数 a[state1]。以这些权重成比例地选择一个 [0, m) 中的整数：[b(0), b(1), ..., b(m−1)]。然后返回所选的整数。

注释

  1. 如果只有两个结果，那么这是特例 Bernoulli 工厂；算法将以 P1(λ) / (P0(λ) + P1(λ)) 的概率返回 1。
  2. 如果 R[j] = choose(n, j)，则步骤 1 到 5 的效果与计算 n 次输入硬币抛掷中的 1 的数量相同（在这种情况下，这将存储在 state1 中），但可惜的是，这些步骤效率不高。在这种情况下，PA 等同于“如果 state 大于 floor(n/2)，则为 1，否则为 state/(n+1−state)”，而 PB 等同于“如果 state 小于 floor(n/2)，则为 1，否则为 (n−state)/(state+1)”。

示例：令 P0(λ) = 2*λ*(1−λ) 且 P1(λ) = (4*λ*(1−λ))2/2。目标是以 P1(λ) / (P0(λ) + P1(λ)) 的概率生成 1。在准备此函数后（并注意到最大次数为 n = 4），我们得到系数总和 R = (0, 2, 12, 2, 0)。由于 R 以 0 开头和结尾，因此步骤 2 将 state1 和 state2 分别设置为第一个或最后一个非零项的位置，即 1 或 3。（或者，因为 R 以 0 开头和结尾，我们可以包含第三个多项式，即常数 P2(λ) = 0.001，这样新的系数总和将是 R′ = (0.001, 10.004, 12.006, 2.006, 0.001) [通过将系数 0.001*choose(n, i) 添加到 i 处的总和（从 i = 0 开始）得到]。现在我们将使用 R′ 运行算法，如果它返回 2 [意味着选择了常数多项式]，我们将重试，直到算法不再返回 2。）

特定代数函数

(Flajolet et al., 2010)(1) 展示了如何通过生成比特串并确定该比特串是否属于特定类别的比特串来模拟某些函数。确定比特串是否属于该类别的规则称为二进制随机文法，它只使用两种“字母”的字母表，或更一般的随机文法。这些函数属于称为代数函数（可以作为多项式方程的解的函数）的类别。

根据 (Mossel 和 Peres 2005)(14)，如果一个函数可以作为系数为有理数的多项式方程的解，则可以通过下推自动机（带有栈的状态机）来模拟一个工厂函数。

以下算法模拟以下代数函数

  • ∑k = 0, 1, 2, ... (λk * (1 − λ) * W(k) / βk)，或者
  • (1 − λ) * OGF(λ/β)，

其中—

  • W(k) 是区间 [0, βk] 内的一个数，在许多情况下是随机文法所能产生的 k 字母单词的数量，
  • β 是字母表大小，即字母表中的“字母”数（例如，Flajolet 论文中讨论的情况为 2），并且是一个大于或等于 2 的整数，
  • 普通生成函数 OGF(x) = W(0) + W(1) * x + W(2) * x2 + W(3) * x3 + ..., 并且
  • 第二个公式结合了该论文定理 3.2 的修正(18)。

(此处，OGF(x) 的 kth 系数对应于 W(k)。) 算法如下。

  1. 反复抛掷输入硬币，直到返回 0。令 g 为此过程中硬币返回 1 的次数。
  2. 返回一个概率为 W(g)/βg 的数（1 或 0）。（在 Flajolet 论文中，这是通过生成一个 g 字母单词并使用二进制随机文法“解析”该单词来确定的，以判断该单词是否能由该文法生成。请注意，这可以通过生成单词的每个“字母”来完成。）

该算法的一个扩展（Flajolet et al. 论文中未提及）是使用比两个字母更大的字母表的随机文法。例如，对于三元随机文法，字母表大小为 3，算法中 β 为 3。通常，对于β-元随机文法，字母表大小为 β，可以是任何大于或等于 2 的整数。

注意：OGF 的收敛半径是使得函数在区间 [0, ρ) 上定义的最大的数 ρ。在此算法中，收敛半径在区间 [1/β, 1] 内（Flajolet 1987）(19)。例如，平方根构造中涉及的 OGF（下文示例中给出）的收敛半径为 1/2。

示例

  1. 以下是 Flajolet et al. 论文中的一个示例。可以如下“解析”一个 g 字母的二进制单词以确定该单词是否编码三元树：“2. 如果 g 为 0，则返回 0。否则，将 i 设置为 1，将 d 设置为 1。; 2a. 生成一个无偏随机比特（即，以相等概率选择 0 或 1），如果该比特为 0，则将 d 减 1，否则加 2。; 2b. 将 i 加 1。然后，如果 i < g 且 d > 0，则转到步骤 3a。; 2c. 如果 d 为 0 且 i 等于 g，则返回 1，否则返回 0。”

  2. 如果 W(g)，即随机文法在问题中可以生成的 g 字母单词的数量，形式为—

    • choose(g, g/t) * (β−1)g−g/t（具有恰好 g/t 个 A 的 g 字母单词的数量，对于字母表大小为 β），如果 g 可被 t 整除(20)，并且
    • 否则为 0，

    其中 t 是大于或等于 2 的整数，β 是字母表大小且大于或等于 2 的整数，则步骤 2 的算法可以如下进行：“2. 如果 g 不能被 t 整除，则返回 0。否则，生成 [0, β) 中的 g 个均匀随机整数（例如，如果 β 是 2，则生成 g 个无偏随机比特），然后如果恰好有 g/t 个零被生成，则返回 1，否则返回 0。”如果 β = 2，则这再现了 Flajolet 论文中的另一个示例。

  3. 如果 W(g) 的形式为—
        choose(g * α, g) * (β−1)g*α−g / βg*α−g,
    其中 α 是大于或等于 1 的整数，β 是字母表大小且大于或等于 2 的整数(21)，步骤 2 的算法可以如下进行：“2. 生成 [0, β) 中的 g * α 个均匀随机整数（例如，如果 β 是 2，则生成 g * α 个无偏随机比特），然后如果恰好生成了 g 个零，则返回 1，否则返回 0。”如果 α = 2 且 β = 2，则这表示 Flajolet et al. 论文中提到的平方根构造 sqrt(1 − λ)。如果 α 为 1，则修改后的算法模拟以下概率：(β*(λ−1))/(λ−β)。有趣的是，我发现如果 α 为 2 或更大，概率将简化为涉及超几何函数。具体来说，概率变为—

    • (1 − λ) * α−1Fα−2(1/α, 2/α, ..., (α−1)/α; 1/(α−1), ..., (α−2)/(α−1); λ * αα/((α−1)α−1 * 2α)) 如果 β = 2，或者更一般地，
    • (1 − λ) * α−1Fα−2(1/α, 2/α, ..., (α−1)/α; 1/(α−1), ..., (α−2)/(α−1); λ*αα*(β−1)α−1/((α−1)α−1 * βα)).

    此修改后算法的普通生成函数因此是—
        OGF(z) = α−1Fα−2(1/α, 2/α, ..., (α−1)/α; 1/(α−1), ..., (α−2)/(α−1); z*αα*(β−1)α−1/((α−1)α−1 * βα−1)).

  4. 示例 2 中涉及的概率同样涉及超几何函数

    • (1 − λ) * t−1Ft−2(1/t, 2/t, ..., (t−1)/t; 1/(t−1), ..., (t−2)/(t−1); λt*tt*(β−1)t−1/((t−1)t−1 * βt)).

特定幂级数

Mendo (2019)(22) 给出了一种 Bernoulli 工厂算法，用于可以将形式为—

    1 − (c[0] * (1 − λ) + ... + c[i] * (1 − λ)i + 1 + ...),

其中 c[i] ≥ 0 是级数的系数，并且总和为 1。（根据 Mendo 的说法，这意味着该级数是可微的 — 其图形没有“尖角” — 并且当 λ 分别接近 0 或 1 时，其值接近 0 或 1。）算法如下。

  1. 令 v 为 1，令 result 为 1。
  2. 将 dsum 设置为 0，将 i 设置为 0。
  3. 抛掷输入硬币。如果返回 v，则返回 result。
  4. 如果 i 等于或大于系数的数量，则将 ci 设置为 0。否则，将 ci 设置为 c[i]。
  5. 以 ci/(1 − dsum) 的概率，返回 1 减去 result。
  6. 将 ci 添加到 dsum，将 1 添加到 i，然后转到步骤 3。

正如 Mendo (2019)(22) 所指出的，该算法的变体适用于形式为—

  1. (c[0] * (1 − λ) + ... + c[i] * (1 − λ)i + 1 + ...)，或
  2. (c[0] * λ + ... + c[i] * λi + 1 + ...)，或
  3. 1 − (c[0] * λ + ... + c[i] * λi + 1 + ...)。

在前两种情况下，将算法中的“令 result 为 1”替换为“令 result 为 0”。在后两种情况下，将“令 v 为 1”替换为“令 v 为 0”。此外，正如 Mendo 所指出的，c[i] 的总和也可以小于 1，在这种情况下，如果算法会返回 1，则它返回一个数，该数以等于所有 c[i] 总和的概率为 1，否则为 0。

(Łatuszyński et al. 2009/2011)(23) 给出了一种算法，适用于一类广泛的级数和其他结构，这些结构从上方和下方收敛到所需的概率。

其中一种情况是当 f(λ) 可以写成—

f(λ) = d[0] − d[1] * λ + d[2] * λ2 − ...,

这是一个交错级数，其中 d[i] 都在区间 [0, 1] 内，形成一个非递增的系数序列，并且 f(1) 必须收敛到一个在半开区间 [0, 1) 内的数。

以下是这种级数的通用算法，称为**通用鞅算法**。它接受一个系数列表和一个输入硬币，并返回 1（概率为上述级数给出的值），否则返回 0。

  1. 令 d[0], d[1] 等为交错级数的第一、第二等系数。将 u 设置为 d[0]，将 w 设置为 1，将 ℓ 设置为 0，将 n 设置为 1。
  2. 生成一个 uniform(0, 1) 随机变量，称之为 ret。
  3. 如果 w 不为 0，则抛掷输入硬币并将结果（0 或 1）乘以 w。
  4. 如果 n 是偶数，则将 u 设置为 ℓ + w * d[n]。否则，将 ℓ 设置为 u − w * d[n]。
  5. 如果 ret 小于（或等于）ℓ，则返回 1。如果 ret 小于 u，则转到下一步。如果都不是，则返回 0。（如果 ret 是均匀 PSRN，这些比较应通过 **URandLessThanReal 算法**进行，该算法在我关于 PSRN 的文章中有所描述。）
  6. 将 n 加 1 并转到步骤 3。

另一种情况是当 f(λ) 可以写成—

f(λ) = d[0] − d[1] * λ2 + d[2] * λ4 − ...,

这同样是一个交错级数，其中 d[i] 都在区间 [0, 1] 内，形成一个非递增的系数序列，并且 f(1) 必须收敛到一个在半开区间 [0, 1) 内的数。在这种情况下，通过在步骤 3 之后添加以下内容来修改通用鞅算法：“3a. 重复步骤 3 一次。”（这种级数的例子可以在 sin(λ) 和 cos(λ) 中找到。）

(Nacu 和 Peres 2005，命题 16)(15)。以下算法模拟形式为—

    d[0] + d[1] * λ + d[2] * λ2 + ...,

其中每个 d[i] 大于或等于 0，并接受两个参数 t 和 ϵ，其中 t 必须选择使得 t 在 (0, 1] 内，f(t) < 1，并且 λ < t − 2*ϵ。

  1. 创建一个 ν 输入硬币，执行以下操作：“(1) 将 n 设置为 0。(2) 以 ϵ/t 的概率转到下一步。否则，将 n 加 1 并重复此子步骤。(3) 以 1 − d[n]*tn 的概率返回 0。(4) 使用 (λ) 输入硬币，x/y = 1/(t − ϵ)，i = 1（对于 2019 算法），以及 ϵ = ϵ，调用下面描述的**2014 算法**、**2016 算法**或**2019 算法**n 次。如果其中任何一个调用返回 0，则返回 0。否则，返回 1。”
  2. 使用前面描述的 ν 输入硬币，x/y = t/ϵ，i = 1（对于 2019 算法），以及 ϵ = ϵ，调用**2014 算法**、**2016 算法**或**2019 算法**一次，并返回结果。

通用工厂函数

具有未知正面概率 λ 的硬币可以转化为具有正面概率 f(λ) 的硬币，其中 f 是任何工厂函数，通过一个算法，该算法基于硬币抛掷的结果构建 f(λ) 的随机边界。这些随机边界来自两个多项式序列

  • 一个多项式序列从上方收敛到 f，另一个从下方收敛。
  • 对于每个序列，多项式必须具有递增的次数。
  • 多项式以Bernstein 形式写成（参见“特定多项式”）。
  • 对于每个序列，次数为 n 的多项式的系数必须位于前一个上方多项式和前一个下方多项式（一旦多项式提升到次数 n）的“内部”。这也称为一致性要求。

本节提出了两种通过多项式模拟工厂函数的算法。在两种算法中

  • fbelow(n, k) 是 Bernstein 形式的次数为 n 的多项式逼近 f 从下方收敛的第 k 个系数的下界，其中 k 在区间 [0, n] 内。例如，它可以是 f(k/n) 减去一个取决于 n 的常数。（见下文注 3。）
  • fabove(n, k) 是 Bernstein 形式的次数为 n 的多项式逼近 f 从上方收敛的第 k 个系数的上界。例如，它可以是 f(k/n) 加上一个取决于 n 的常数。（见下文注 3。）

第一个算法实现了 Łatuszyński et al. (2009/2011)(23) 中的反向时间鞅框架（算法 4）和 Flegal 和 Herbei (2012)(24) 的算法 I 中的次数加倍建议，尽管下面指出了算法 I 中的一个错误。第一个算法如下。

  1. 生成一个 uniform(0, 1) 随机变量，称之为 ret。
  2. 将 ℓ 和 ℓt 设置为 0。将 u 和 ut 设置为 1。将 lastdegree 设置为 0，将 ones 设置为 0。
  3. 将 degree 设置为使得第一对多项式的次数等于 degree 并且所有系数都在 [0, 1] 内。例如，可以如下完成：令 fbound(n) 为 fbelow(n, k) 的最小值和 fabove(n,k) 的最大值（对于 k 在区间 [0, n] 内）；然后设置 degree 为 1；然后当 fbound(degree) 返回小于 0 或大于 1 的上界或下界时，将 degree 乘以 2；然后转到下一步。
  4. 将 startdegree 设置为 degree。
  5. (剩余步骤现在将重复进行，直到算法通过返回值结束。) 抛掷输入硬币 t 次，其中 t 为 degree − lastdegree。对于每次硬币返回 1 的情况，将 ones 加 1。
  6. 计算 ℓ 和 u 如下：
    1. 定义 FB(a, b) 如下：令 c 为 choose(a, b)。计算 fbelow(a, b) 作为下界和上界 LB 和 UB，它们足够精确，使得 floor(LB*c) = floor(UB*c)，然后返回 floor(LB*c)。
    2. 定义 FA(a, b) 如下：令 c 为 choose(a, b)。计算 fabove(a, b) 作为下界和上界 LB 和 UB，它们足够精确，使得 ceil(LB*c) = ceil(UB*c)，然后返回 ceil(LB*c)。
    3. 令 c 为 choose(degree, ones)。设置 ℓ 为 (FB(degree, ones))/c，设置 u 为 (FA(degree, ones))/c。
  7. (此步骤和下一步找到当前硬币抛掷给定下的先前 ℓ 和 u 的期望值。) 如果 degree 等于 startdegree，则将 ℓs 设置为 0，将 us 设置为 1。（Flegal 和 Herbei 2012 的算法 I 未考虑这一点。）
  8. 如果 degree 大于 startdegree：令 nh 为 choose(degree, ones)，令 od 为 degree/2。设置 ℓs 为 ∑j=0,...,ones FB(od,j)*choose(degree−od, ones−j)/nh，并设置 us 为 ∑j=0,...,ones FA(od,j)*choose(degree−od, ones−j)/nh。
  9. 令 m 为 (ut−ℓt)/(us−ℓs)。设置 ℓt 为 ℓt+(ℓ−ℓs)*m，并设置 ut 为 ut−(us−u)*m。
  10. 如果 ret 小于（或等于）ℓt，则返回 1。如果 ret 小于 ut，则转到下一步。如果都不是，则返回 0。（如果 ret 是均匀 PSRN，这些比较应通过 **URandLessThanReal 算法**进行，该算法在我关于 PSRN 的文章中有所描述。）
  11. (找到下一对多项式并重新启动循环。) 将 degree 增加，使得下一对多项式的次数等于更高的 degree 值并更接近目标函数（例如，将 degree 乘以 2）。然后，转到步骤 5。

第二个算法由 Thomas 和 Blanchet (2012)(10) 给出；它假定与前一个算法相同的多项式序列可用。下面给出了一个等同于该算法的算法。

  1. 将 ones 设置为 0，将 lastdegree 设置为 0。
  2. 将 degree 设置为使得第一对多项式的次数等于 degree 并且所有系数都在 [0, 1] 内。例如，可以如下完成：令 fbound(n) 为 fbelow(n, k) 的最小值和 fabove(n,k) 的最大值（对于 k 在区间 [0, n] 内）；然后设置 degree 为 1；然后当 fbound(degree) 返回小于 0 或大于 1 的上界或下界时，将 degree 乘以 2；然后转到下一步。
  3. 将 startdegree 设置为 degree。
  4. (剩余步骤现在将重复进行，直到算法通过返回值结束。) 抛掷输入硬币 t 次，其中 t 为 degree − lastdegree。对于每次硬币返回 1 的情况，将 ones 加 1。
  5. 将 c 设置为 choose(degree, ones)。可选地，将 c 乘以 2degree（见下文注 3）。
  6. 如下找到 acount 和 bcount：
    1. 计算 fbelow(degree, ones) 作为下界和上界 LB 和 UB，它们足够精确，使得 floor(LB*c) = floor(UB*c)。然后设置 a[degree,ones] 和 acount 为 floor(LB*c)。
    2. 计算 1−fabove(degree, ones) 作为下界和上界 LB 和 UB，它们足够精确，使得 floor(LB*c) = floor(UB*c)。然后设置 b[degree,ones] 和 bcount 为 floor(LB*c)。
    3. 从 c 中减去 (acount + bcount)。
  7. 如果 degree 大于 startdegree，则
    1. 令 diff 为 degree−lastdegree，令 u 为 max(0, ones−lastdegree)，令 v 为 min(ones, diff)。 （以下子步骤从 acount 和 bcount 中移除本应更早终止算法的结果。此过程与论文第 3 节的步骤 (f) 不同，该步骤似乎不正确，并且该过程是从我请求 A. C. Thomas 提供的补充源代码派生的。）
    2. 令 g 为 choose(lastdegree, ones−u)。令 h 为 1。如果 c 如步骤 5 中那样乘以 2，则将 h 乘以 2lastdegree（见下文注 3）。
    3. 对于区间 [u, v] 内的每个整数 k
      1. 令 d 为 choose(diff, k)。令 ω 为 ones−k。
      2. 如果尚未计算：计算 fbelow(lastdegree, ω) 作为下界和上界 LB 和 UB，它们足够精确，使得 floor(LB*g*h) = floor(UB*g*h)。然后设置 a[lastdegree, ω] 为 floor(LB*g*h)。
      3. 如果尚未计算：计算 1−fabove(lastdegree, ω) 作为下界和上界 LB 和 UB，它们足够精确，使得 floor(LB*g*h) = floor(UB*g*h)。然后设置 b[lastdegree, ω] 为 floor(LB*g*h)。
      4. 从 acount 中减去 (a[lastdegree, ω]*d)。
      5. 从 bcount 中减去 (b[lastdegree, ω]*d)。
      6. 将 g 乘以 ω，然后将 g 除以 (lastdegree+1−ω)。 (将 g 设置为 choose(lastdegree, (ones−k)−1)。)
  8. 以这些权重成比例地选择 0、1 或 2：[acount, bcount, c]。
  9. 如果上一步选择的数字是 0，则返回 1。如果上一步选择的数字是 1，则返回 0。
  10. (找到下一对多项式并重新启动循环。) 将 lastdegree 设置为 degree，然后增加 degree，使得下一对多项式的次数等于更高的 degree 值并更接近目标函数（例如，将 degree 乘以 2）。然后，转到步骤 4。

注释

  1. 这两个算法的效率取决于许多因素，包括 f 的“平滑度”以及计算 fbelow 和 fabove 的适当值有多容易。为给定函数 f 实现 fbelow 和 fabove 的最佳方法需要对该函数进行深入的数学分析。有关更多信息，请参阅我的关于 Bernoulli 工厂的补充说明。
  2. 在某些情况下，一对多项式序列可能无法快速收敛到所需函数 f，特别是当 f 不够“平滑”时。Thomas 和 Blanchet (2012)(10) 的一个有趣建议是使用多对收敛到 f 的多项式序列，其中每对都针对 λ 的特定范围进行了优化：首先抛掷输入硬币几次以获得 λ 的粗略估计，然后选择针对估计 λ 优化的对，并在该对上运行本节中的任一算法。
  3. 第二个算法（在 Thomas 和 Blanchet (2012)(10) 中给出）基于 Nacu 和 Peres (2005)(15) 的算法。在这两篇论文中，该算法仅在 λ 在区间 (0, 1) 内时才有效。如果 λ 可以是 0 或 1（意味着输入硬币允许每次返回 1 或每次返回 0），那么根据 Holtz et al. (2011)(25) 的建议，步骤 5 中的 c 可以乘以 2degree，步骤 7 子步骤 2 中的 h 乘以 2lastdegree，以确保对所有 λ 值都正确。

通用无理数常量的算法

本节展示了通用算法，用于生成正面概率等于一个无理数（不是两个整数之比的数）的概率，该数已知其数字或级数展开、连分式或连对数。

但另一方面，有理数常量的概率很容易模拟。如果提供了公平硬币，则应使用我关于随机采样方法的文章中描述的 ZeroOrOne 方法。如果提供了正面概率未知的硬币，则应使用随机性提取方法将其转换为公平硬币。

数字展开

概率可以表示为数字展开（形式为 0.dddddd...）。以下算法返回 1 的概率为 p，否则返回 0，其中 p 是半开区间 [0, 1) 中的概率。请注意，数字 0 也是零的无限数字展开，数字 1 也是负一的无穷大数字展开。无理数始终具有无限的数字展开，必须“即时”计算。

在算法中（另见 (Brassard et al., 2019)(26), (Devroye 1986, p. 769)(27)），BASE 是数字基数，例如二进制为 2，十进制为 10。

  1. 将 u 设置为 0，将 k 设置为 1。
  2. 将 u 设置为 (u * BASE) + v，其中 v 是 [0, BASE) 区间内的均匀随机整数（如果 BASE 为 2，则 v 只是一个无偏随机比特）。计算 pa，它是 p 的近似值，使得 abs(p−pa) ≤ BASE−k。将 pk 设置为 pa 的数字展开，小数点后最多 k 位。例如：如果 p 是 π/4，BASE 是 10，k 是 5，则 pk = 78539。
  3. 如果 pk + 1 <= u，则返回 0。如果 pk - 2 >= u，则返回 1。如果都不是，则将 1 添加到 k 并转到步骤 2。

连分式

以下算法模拟了以下形式的简单连分数表示的概率：0 + 1 / (a[1] + 1 / (a[2] + 1 / (a[3] + ... )))。a[i] 是部分分母，它们的绝对值都不能小于 1。受 (Flajolet et al., 2010, "Finite graphs (Markov chains) and rational functions")(1) 的启发，我开发了以下算法。

算法 1。此算法仅当每个 a[i] 的绝对值大于或等于 1 且 a[1] 大于 0 时才有效，但否则，每个 a[i] 可以是负数和/或非整数。算法从 pos 等于 1 开始。然后执行以下步骤。

  1. 将 k 设置为 a[pos]。
  2. 如果 pos 处的部分分母是最后一个，则返回一个概率为 1/k 的数（1 或 0）。
  3. 如果 a[pos] 小于 0，则将 kp 设置为 k − 1，将 s 设置为 0。否则，将 kp 设置为 k，将 s 设置为 1。（此步骤处理负的部分分母。）
  4. 执行以下过程，直到此算法运行返回一个值
    1. 以 kp/(1+kp) 的概率，返回一个概率为 1/kp 的数（1 或 0）。
    2. 运行当前正在运行的算法的一个独立运行，但 pos = pos + 1。如果独立运行返回 s，则返回 0。

广义连分数的形式为 0 + b[1] / (a[1] + b[2] / (a[2] + b[3] / (a[3] + ... )))。a[i] 与之前相同，但 b[i] 是部分分子。以下是模拟广义连分数形式的概率的两种算法。

算法 2。此算法仅当每个比率 b[i]/a[i] 小于或等于 1 时才有效，但否则，每个 b[i] 和每个 a[i] 可以是负数和/或非整数。此算法采用一种等价变换，将广义连分数转换为简单连分数。算法从 pos 和 r 都等于 1 开始。然后执行以下步骤。

  1. 将r设为 1 / (r * b[pos])，然后将k设为a[pos] * r。（k是等价的简单连分数的部分分母。）
  2. 如果pos处的部分分子/分母对是最后一个，则返回一个概率为 1/abs(k) 为 1，否则为 0 的数字。
  3. 将kp设为abs(k)，将s设为 1。
  4. 将r2设为 1 / (r * b[pos + 1])。如果a[pos + 1] * r2小于 0，将kp减 1，将s设为 0。（此步骤处理负数部分分子和分母。）
  5. 执行以下过程，直到此算法运行返回一个值
    1. 以 kp/(1+kp) 的概率，返回一个概率为 1/kp 的数（1 或 0）。
    2. 对当前正在运行的算法进行单独运行，但将pos设为pos + 1，将r设为r。如果单独运行返回s，则返回 0。

算法 3。此算法仅在每个比率b[i]/a[i]小于等于 1 且每个b[i]和每个a[i]大于 0 时有效，但否则，每个b[i]和每个a[i]都可以是非整数。该算法从等于 1 的pos开始。然后执行以下步骤。

  1. 如果pos处的部分分子/分母对是最后一个，则返回一个概率为 b[pos]/a[pos] 为 1，否则为 0 的数字。
  2. 执行以下过程，直到此算法运行返回一个值
    1. 以 a[pos]/(1 + a[pos]) 的概率，返回一个概率为 b[pos]/a[pos] 为 1，否则为 0 的数字。
    2. 对当前正在运行的算法进行单独运行，但将pos设为pos + 1。如果单独运行返回 1，则返回 0。

算法 3 的正确性证明请参见附录。

注释

  • 如果这些算法中的任何一个遇到一个不在 [0, 1] 区间内的概率，则整个算法将针对该连分数失败。

  • 如果——，这些算法将适用于“1 − ...”形式（而不是“0 + ...”）的连分数

    • 在运行算法之前，删除第一个部分分子和分母的符号，并且
    • 在运行算法之后，取 1 减去结果（而不是仅结果）。

  • 这些算法旨在允许“即时”计算部分分子和分母。

  • 以下是算法 1 的另一种写法，它更好地展示了其灵感，因为它展示了如何使用“偶数校验构造”（或双硬币特殊情况）以及“1 − x”构造来开发与连分数的展开一样大的有理数模拟器，正如 Flajolet 论文的引用部分所建议的。但是，仅当连分数展开的大小（此处为size）提前知道时，它才有效。
    1. 将i设为size。
    2. 创建一个输入硬币，该硬币执行以下操作：“返回一个概率为 1/a[size] 为 1，否则为 0 的数字”。
    3. 当i等于 1 或更大时
      1. 将k设为a[i]。
      2. 创建一个输入硬币，该硬币接收之前的输入硬币和k，并执行以下操作：“(a) 以 k/(1+k) 的概率，返回一个概率为 1/k 为 1，否则为 0 的数字；(b) 翻转之前的输入硬币。如果结果是 1，则返回 0。否则，转到步骤 (a)”。（概率 k/(1+k) 与 λ/(1+λ) = 1 − 1/(1+λ) 相关，后者涉及 1/(1+λ) 的偶数校验构造——或双硬币特殊情况——以及“1 − x”的补集。）
      3. 将 i 减 1。
    4. 翻转此算法创建的最后一个输入硬币，并返回结果。

连对数

一个数在开区间 (0, 1) 中的连对数（Gosper 1978）(28)，（Borwein 等人，2016）(29) 具有以下连分数形式：0 + (1 / 2c[1]) / (1 + (1 / 2c[2]) / (1 + ...))，其中c[i]是连对数的系数，且都为 0 或更大。我提出了以下算法，该算法模拟一个表示为连分数展开的概率。

该算法从等于 1 的pos开始。然后执行以下步骤。

  1. 如果pos处的系数是最后一个，则返回一个概率为 1/(2c[pos]) 为 1，否则为 0 的数字。
  2. 执行以下过程，直到此算法运行返回一个值
    1. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则返回一个概率为 1/(2c[pos]) 为 1，否则为 0 的数字。
    2. 对当前正在运行的算法进行单独运行，但将pos设为pos + 1。如果单独运行返回 1，则返回 0。

有关正确性证明，请参见附录。

特定收敛级数

Mendo (2020)(30) 给出了一个通用算法，该算法可以模拟开区间 (0, 1) 中的任何概率，只要它可以重写为收敛级数——

  • 该级数的形式为 a[0] + a[1] + ...，其中a[n]均为大于 0 的有理数，并且
  • 存在一个序列err[0], err[1], ...，该序列是非递增的并收敛到 0，其中err[n]是截断级数a的前n+1项的和所产生的误差的上界。

算法如下。

  1. 将ϵ设为 1，然后将n, lamunq, lam, s和k全部设为 0。
  2. 将 1 加到k，然后将s/(2k) 加到lam。
  3. 如果lamunq+ϵ ≤ lam + 1/(2k)，则转到步骤 8。
  4. 如果lamunq > lam + 1/(2k)，则转到步骤 8。
  5. 如果lamunq > lam + 1/(2k+1) 且lamunq+ϵ < 3/(2k+1)，则转到步骤 8。
  6. 将a[n]加到lamunq，并将ϵ设为err[n]。
  7. 将 1 加到n，然后转到步骤 3。
  8. 令bound为lam+1/(2k)。如果lamunq+ϵ ≤ bound，则将s设为 0。否则，如果lamunq > bound，则将s设为 2。否则，将s设为 1。
  9. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则转到步骤 2。否则，返回一个数字：如果s为 0，则返回 0；如果s为 2，则返回 1；否则，返回一个无偏随机比特（0 或 1 的概率相等）。

如果上面的a是一个收敛到开区间 (0, 1) 中概率的基-2 对数的序列，那么我开发的以下算法可以模拟该概率。为简单起见，即使这种概率的对数是负数，所有a[i]都必须为 0 或更大（因此是已为负的对数近似值的负值），并且必须形成一个非递减序列，所有err[i]都必须为 0 或更大。

  1. 将intinf设为 floor(max(0, abs(a[0])))。（这是级数第一项的绝对整数部分，或 0，取较大者。）
  2. 如果intinf大于 0，则生成无偏随机比特，直到生成零比特或生成了intinf个比特。如果生成了零比特，则返回 0。
  3. 生成一个速率为 ln(2) 的指数随机变量E。例如，可以使用 "任意精度采样中的更多算法" 中给出的算法。（我们利用了指数分布的无记忆性：给定一个指数随机变量E大于intinf，E减去intinf具有相同的分布。）
  4. 将n设为 0。
  5. 重复执行以下过程，直到算法返回一个值
    1. 将inf设为 max(0, a[n])，然后将sup设为 min(0, inf+err[n])。
    2. 如果E小于inf+intinf，则返回 0。如果E小于sup+intinf，则转到下一步。如果都不是，则返回 1。
    3. 将n设为 1。

当a收敛到自然对数而不是基-2 对数时，情况就微不足道了。同样，对于此算法，所有a[i]都必须为 0 或更大，并且形成一个非递减序列，所有err[i]都必须为 0 或更大。

  1. 生成一个指数随机变量E（速率为 1）。
  2. 将n设为 0。
  3. 重复执行以下过程，直到算法返回一个值
    1. 将inf设为 max(0, a[n])，然后将sup设为 min(0, inf+err[n])。
    2. 如果E小于inf+intinf，则返回 0。如果E小于sup+intinf，则转到下一步。如果都不是，则返回 1。
    3. 将n设为 1。

示例:

  • 令f(λ) = cosh(1)−1。本节的第一个算法可以模拟此常数，方法是修改步骤 6，使其读作：“令m为((n+1)*2)，令α为 1/(m!)（泰勒级数的一项）。将α加到lamunq，并将ϵ设为 2/((m+1)!)（误差项）”。(31)
  • 当n可以非常大（例如，高达 230）时，对数可以作为模拟概率z = choose(n, k)/2n的高效算法的基础，而无需依赖浮点算术。在此示例中，choose(n, k) 的二项式系数的平凡算法通常需要随n和k增长的存储空间。另一方面，任何常数都可以平均使用最多两个无偏随机比特来模拟，对于此处涉及的常数，甚至少于这个数量（Kozen 2014）(32)。与其直接计算二项式系数，不如计算一个级数，该级数收敛到该系数的对数，例如 ln(choose(n, k))，即使对于较大的n和k，它在空间上也经济高效。然后可以使用该级数与上述算法来模拟概率z。类似的方法也已实现（参见 interval.py 和 betadist.py）。另请参见（Bringmann 等人 2014）(33) 中的附录。

其他通用算法

 

凸组合

假设我们有一个或多个输入硬币hi(λ)，它们以取决于λ的概率返回正面。（硬币数量可以是无限的。）以下算法随机选择其中一个硬币，然后翻转该硬币。具体来说，该算法以等于以下加权和的概率生成 1：g(0) * h0(λ) + g(1) * h1(λ) + ...，其中g(i)是选择硬币i的概率，hi是硬币i模拟的函数，所有g(i)的总和为 1。参见（Wästlund 1999，定理 2.7）(8)。（或者，该算法可以看作是以降“hX(λ)”的期望或平均值，即X是识别随机选择的硬币的数字的期望或平均值，来返回正面。）

  1. 以某种方式生成一个随机整数X。例如，它可以是 [1, 6] 中的均匀随机整数，或者它可以是泊松随机变量。（具体来说，数字X以g(X)的概率生成。）
  2. 翻转由X表示的硬币并返回结果。

注释

  1. 构建凸组合。 假设我们有一个函数f(λ) = ∑n=0,1,... wn(λ)，其中wn是连续函数，其在域 [0, 1] 中的最大值总和为 1 或更少。令g(n)为随机选择的数字X为n的概率。那么通过**生成X并以wX(λ)/g(X)的概率翻转硬币**，我们可以将概率f(λ)模拟为凸组合——

    f(λ) = ∑n=0,1,... g(n) * (wn(λ) / g(n)),

    但这仅在以下两个条件对每个整数n≥0都满足时才有效

    • g(n) ≥ wn(λ) ≥ 0 对区间 [0, 1] 中的每个λ都成立（大致意思是wn有界于上或被g(n)“支配”）。
    • 函数wn(λ)/g(n)允许一个伯努利工厂（如果它在区间 (0, 1) 内触及 0 或 1，但不是常数，例如，它不会）。

    另请参见 Mendo (2019)(22)。

  2. 具有非负级数展开的常数。 注意 1 的特殊情况。令g如注意 1 所示。假设我们有一个常数，其级数展开为：c = a0 + a1 + a2 + ...，其中——

    • an均为 0 或更大，且总和为 1 或更少，并且
    • g(n) ≥ an ≥ 0 对每个整数n≥0都成立。

    那么通过**生成X并以aX/g(X)的概率翻转硬币**，我们可以将概率c模拟为凸组合——

    c = ∑n=0,1,... g(n) * (an / g(n)).

示例

  1. 生成一个泊松(μ)随机变量X，然后翻转输入硬币。以 1/(1+X) 的概率，返回翻转硬币的结果；否则，返回 0。这对应于g(i)是泊松(μ)概率，而硬币hi以 1/(1+i) 的概率返回 1，否则返回 0。此方法返回 1 的概率为 **E**[1/(1+X)]，即 (exp(μ)−1)/(exp(μ)*μ)。
  2. 生成一个泊松(μ)随机变量X，如果X为 0，则返回 1，否则返回 0。这是一个先前提到的 exp(−μ) 的伯努利工厂，对应于g(i)是泊松(μ)概率，而硬币hi在i为 0 时返回 1，否则返回 0。
  3. 生成一个泊松(μ)随机变量X，运行**exp(−z)的算法**，其中z = X，并返回结果。返回 1 的概率为 **E**[exp(−X)]，即 exp(μ*exp(−1)−μ)。下面的 Python 代码使用计算机代数库 SymPy 来查找此概率：from sympy.stats import *; E(exp(-Poisson('P', x))).simplify()。
  4. 伯努利竞赛 (Dughmi 等人 2017)(34)：假设我们有n个硬币，然后随机均匀地选择其中一个并翻转该硬币。如果翻转返回 1，则返回X；否则，重复此算法。此算法根据其正面概率选择一个随机硬币。每次迭代对应于g(i)为 1/n，而hi()为相应硬币i的概率。
  5. (Wästlund 1999)(8)：生成一个均值为 1 的泊松随机变量X，然后将输入硬币翻转X次。如果任何一次翻转返回 1，则返回 0；否则返回 1。这是一个 exp(−λ) 的伯努利工厂，对应于g(i)是泊松概率，即 1/(i!*exp(1))，而hi()是 (1−λ)i。
  6. 多元伯努利工厂 (Huber 2016)(35)，形式为 R = C0*λ0 + C1*λ1 + ... + Cm−1*λm−1，其中Ci是已知的常数，大于 0，且对于任何ϵ > 0，R ≤ 1 − ϵ：随机均匀地选择 [0, m) 中的一个整数，称之为i，然后运行 (m*Ci)*λi 的线性伯努利工厂。这不同于 Huber 的“稀疏化”泊松过程的建议，该过程由多个输入硬币驱动。
  7. 概率生成函数 (PGF) (Dughmi 等人 2017)(34)。生成正面，其概率为 **E**[λX]，即λX的期望或平均值。 **E**[λX]是X分布的 PGF。算法如下：(1) 以某种方式生成一个随机整数X；(2) 翻转输入硬币，直到翻转返回 0 或翻转了X次硬币，以先到者为准。如果所有翻转（包括最后一次）都返回 1（或如果X为 0），则返回 1；否则返回 0。
  8. 假设X是显示 0 的无偏随机比特的数量，直到生成第一个 1。那么g(n) = 1/(2n+1)。

积分

粗略地说，函数f(x)在区间 [a, b] 上的积分是该函数图形下的面积，当函数限制在该区间内时。

算法 1。 (Flajolet 等人，2010)(1) 展示了如何将模拟f(λ)的算法转换为模拟以下概率的算法

  • (1/λ) * ∫[0, λ] f(u) du，或者等价地，
  • ∫[0, 1] f(u * λ) du（一个积分），

使用以下算法

  1. 生成一个均匀(0, 1)随机变量u。
  2. 创建一个输入硬币，执行以下操作：“翻转原始输入硬币，然后从数字u采样。如果调用和翻转都返回 1，则返回 1，否则返回 0。”
  3. 运行原始伯努利工厂算法，使用步骤 2 中描述的输入硬币而不是原始输入硬币。返回该运行的结果。

算法 2。 算法 1 的一个特殊情况是积分 ∫[0, 1] f(u) du，当原始输入硬币始终返回 1 时

  1. 生成一个均匀(0, 1)随机变量u。
  2. 创建一个输入硬币，执行以下操作：“从数字u采样并返回结果。”
  3. 运行原始伯努利工厂算法，使用步骤 2 中描述的输入硬币而不是原始输入硬币。返回该运行的结果。

算法 3。 我发现可以模拟以下积分——

  • ∫[a, b] f(u) du,

其中 [a, b] 是 [0, 1] 或其中的一个闭区间，使用以下算法

  1. 生成一个均匀(0, 1)随机变量u。然后如果u小于a或大于b，则重复此步骤。（如果u是均匀的PSR N，这些比较应该通过 **URandLessThanReal** 算法完成。）
  2. 创建一个输入硬币，执行以下操作：“从数字u采样并返回结果。”
  3. 运行原始伯努利工厂算法，使用步骤 2 中描述的输入硬币。如果运行返回 0，则返回 0。否则，生成一个均匀(0, 1)随机变量v并返回一个数字：如果v小于a或大于b，则返回 0；否则返回 1。

注意：如果a为 0 且b为 1，则此算法模拟的概率将是单调递增的（将不断上升），斜率不超过 1，并且在点 0 处等于 0。

广义 Bernoulli 竞赛

以下算法（Schmon 等人 2019）(36) 推广了“凸组合”部分中的伯努利竞赛。它以概率返回i——

  • φi = g(i)*hi(λ) / ∑k=0,...,r g(k)*hk(λ)，

其中

  • r 是一个大于 0 的整数。该算法可以从中选择 r+1 个值。
  • g(i) 接受一个整数i并返回一个 0 或更大的数字。这充当“硬币”标记为i的权重；权重越高，“硬币”被“翻转”的可能性越大。
  • hi(λ) 接收一个数字i和输入硬币的正面概率，并返回一个在 [0, 1] 区间内的数字。这代表了 r+1 个选择中的一个“硬币”。

算法如下。

  1. 随机选择 [0, r] 中的一个整数，其概率与以下权重成比例：[g(0), g(1), ..., g(r)]。称选择的整数为i。
  2. 运行一个模拟 hi(λ) 的伯努利工厂算法。如果运行返回 0，则转到步骤 1。
  3. 返回i。

注释

  1. 伯努利竞赛（参见“凸组合”）是此算法的一个特殊情况，其中对于每个k，g(k) = 1。
  2. 如果我们定义S为 [0, r] 中的整数子集，并将步骤 3 替换为“如果i在集合S中，则返回 1。否则，返回 0。”，则该算法以概率 ∑k in S φk 返回 1，否则返回 0。在这种情况下，修改后的算法具有 Agrawal 等人 (2021, 附录 D)(37) 所谓的“骰子硬币算法”，其中——
    g(k) = ck*dr−k，
    hk(λ, μ) = λk*μr−k （对于以下算法：翻转λ硬币k次和μ硬币r−k次；如果所有翻转都返回 1，则返回 1，否则返回 0），并且
    S是闭区间 [1, r]，
    其中c≥0, d≥0，而λ和μ是两个输入硬币的正面概率。在该论文中，c, d, λ, 和μ分别对应于cy, cx, py, 和px，在“门户”算法中，或者分别对应于cy, cx, py, 和px，在“翻转门户”算法中。
  3. 虽然 Schmon 论文中未提及，但算法中的r可以是无穷大（另见 Wästlund 1999，定理 2.7(8)）。在这种情况下，步骤 1 更改为“随机选择一个大于等于 0 的整数，其概率为g(k)，对于整数k。称选择的整数为i。”作为一个例子，步骤 1 可以从泊松分布采样，该分布可以取任何大于等于 0 的整数。

特定函数 λ 的算法

本节描述了特定函数的算法，特别是当它们比通用算法具有更方便的模拟时。

exp(−λ)

此算法改编自通用鞅算法（在“某些幂级数”中，参见上面），并利用了 exp(−λ) 可以重写为 1 − λ + λ2/2 − λ3/6 + λ4/24 − ...的事实，这是一个交错级数，其系数为 1, 1, 1/(2!), 1/(3!), 1/(4!), ...。此算法在开区间 (0, 1) 的任何地方都快速收敛。（换句话说，该算法是均匀快速的，意味着对于任何λ和其他参数的选择，平均运行时间都是有限的（Devroye 1986，特别是第 717 页）(27)。(38)）

  1. 将u设为 1，将w设为 1，将ℓ设为 0，并将n**设为 1**。
  2. 生成一个 uniform(0, 1) 随机变量，称之为 ret。
  3. 重复执行以下过程，直到此算法返回一个值
    1. 如果w不为 0，则翻转输入硬币，将w乘以翻转的结果，然后将w除以n。（这与通用鞅算法不同，是为了更有效地考虑第二项及以后系数中的阶乘。）
    2. 如果n是偶数，则将u设为ℓ + w。否则，将ℓ设为u − w。
    3. 如果ret小于（或等于）ℓ，则返回 1。如果ret小于u，则转到下一个子步骤。如果都不是，则返回 0。（如果ret是均匀的PSR N，这些比较应该通过 **URandLessThanReal** 算法完成，该算法在我的关于 PSRN 的文章中有所描述。）
    4. 将 1 加到n。

(exp(λ)−1)/exp(−λ)

此算法使用通用鞅算法；此函数可以重写为 λ*(1 − λ/2 + λ2/6 − λ3/24 + ...)，其中包含一个交错级数，其系数为 1, 1/(2!), 1/(3!), 1/(4!), ...。

  1. 翻转输入硬币。如果它返回 0，则返回 0。
  2. 将u设为 1，将w设为 1，将ℓ设为 0，并将n**设为 2**。
  3. 生成一个 uniform(0, 1) 随机变量，称之为 ret。
  4. 重复执行以下过程，直到此算法返回一个值
    1. 如果w不为 0，则翻转输入硬币，将w乘以翻转的结果，然后将w除以n。（这与通用鞅算法不同，是为了更有效地考虑第二项及以后系数中的阶乘。）
    2. 如果n是偶数，则将u设为ℓ + w。否则，将ℓ设为u − w。
    3. 如果ret小于（或等于）ℓ，则返回 1。如果ret小于u，则转到下一个子步骤。如果都不是，则返回 0。（如果ret是均匀的PSR N，这些比较应该通过 **URandLessThanReal** 算法完成，该算法在我的关于 PSRN 的文章中有所描述。）
    4. 将 1 加到n。

exp(−(λk * c))

在以下算法中，k是一个大于等于 0 的整数，c是一个实数。

第一个算法在 **c大于等于 0** 时有效。

  1. 特殊情况：如果c为 0，则返回 1。如果k为 0，则运行（本页稍后给出）**exp(−z)的算法**，其中z = c，并返回结果。
  2. 生成一个泊松(c)随机整数，称之为N。（关于生成此整数的精确方法，请参见关于 von Neumann 模式的附录。）
  3. 将i设为 0，然后当i < N时
    1. 翻转输入硬币，直到翻转返回 0 或翻转了k次硬币，以先到者为准。如果所有硬币翻转（包括最后一次）都返回 1，则返回 0。
    2. 将 1 加到i。
  4. 返回 1。

第二个算法应用了通用鞅算法，但仅在 **c是区间 [0, 1] 中的有理数** 时有效。目标函数表示为级数 1 − λk*c + λ2*k*c/2! − λ3*k*x/3!, ...，系数为 1, c, c/(2!), c/(3!), ...。

  1. 特殊情况：如果c为 0，则返回 1。如果k为 0，则运行（本页稍后给出）**exp(−x/y)的算法**，其中x/y = c，并返回结果。
  2. 将u设为 1，将w设为 1，将ℓ设为 0，并将n设为 1。
  3. 生成一个 uniform(0, 1) 随机变量，称之为 ret。
  4. 如果w不为 0，则翻转输入硬币k次或直到翻转返回 0。如果任何一次翻转返回 0，则将w设为 0；如果所有翻转都返回 1，则将w除以n。然后，将w乘以一个概率为c为 1，否则为 0 的数字。
  5. 如果n是偶数，则将u设为ℓ + w。否则，将ℓ设为u − w。
  6. 如果 ret 小于（或等于）ℓ，则返回 1。如果 ret 小于 u，则转到下一步。如果都不是，则返回 0。（如果 ret 是均匀 PSRN，这些比较应通过 **URandLessThanReal 算法**进行，该算法在我关于 PSRN 的文章中有所描述。）
  7. 将 1 加到n并转到步骤 4。

第三个算法基于第二个算法，并且在 **c是有理数且大于等于 0** 时有效。

  1. 令m为 floor(c)。调用第二个算法m次，其中k = k，c = 1。如果任何一次调用返回 0，则返回 0。
  2. 如果c是整数，则返回 1。
  3. 调用第二个算法一次，其中k = k，c = c − floor(c)。返回该调用的结果。

示例： **exp(−((1−λ)1 * c))** ((Dughmi 等人 2017)(34)；对输入硬币应用指数权重——此处为c)：“(1) 如果c为 0，则返回 1；(2) 生成一个泊松(c)随机整数，称之为N；(3) 翻转输入硬币，直到翻转返回 0 或翻转了N次硬币，以先到者为准，然后返回一个数字：如果N为 0 或所有硬币翻转（包括最后一次）都返回 1，则返回 1，否则返回 0。”

exp(−(λ + m)k)

在以下算法中，m和k均为大于等于 0 的整数，除非另有说明。

  1. 如果k为 0，则运行（本页稍后给出）**exp(−x/y)的算法**，其中x/y = 1/1，并返回结果。
  2. 如果k为 1 且m为 0，则运行**exp(−λ)的算法**并返回结果。
  3. 如果k为 1 且m大于 0（在此情况下，m可以是任何有理数）
    • 运行（本页稍后给出）**exp(−x/y)的算法**，其中x/y = m。如果算法返回 0，则返回 0。否则，返回**exp(−λ)的算法**的结果。
  4. 运行（本页稍后给出）**exp(−x/y)的算法**，其中x/y = mk / 1。如果算法返回 0，则返回 0。
  5. 运行**exp(−(λk * c))的算法**，其中k = k，x = 1。如果算法返回 0，则返回 0。
  6. 如果m为 0，则返回 1。
  7. 将i设为 1，然后当i < k时
    1. 令z为 choose(k, i) * mk − i。
    2. 运行**exp(−(λk * c))的算法**z次，其中k = i，x = 1。如果任何一次调用返回 0，则返回 0。
    3. 将 1 加到i。
  8. 返回 1。

exp(λ)*(1−λ)

(Flajolet 等人，2010)(1)

  1. 将k和w都设为 0。
  2. 翻转输入硬币。如果它返回 0，则返回 1。
  3. 生成一个均匀(0, 1)随机变量U。
  4. 如果k > 0 且w小于U，则返回 0。
  5. 将w设为U，将 1 加到k，然后转到步骤 2。

1/(2k + λ) 或 exp(−(k + λ)*ln(2))

这个新算法使用了基-2 对数k + λ，其中k是一个大于等于 0 的整数，并且当这个对数非常大时很有用。

  1. 如果k > 0，则生成无偏随机比特，直到生成零比特或生成了k个比特，以先到者为准。如果生成了零比特，则返回 0。
  2. 创建一个输入硬币μ，执行以下操作：“翻转输入硬币，然后运行（稍后给出）**ln(1+y/z)的算法**，其中y/z = 1/1。如果调用和翻转都返回 1，则返回 1。否则，返回 0。”
  3. 使用μ输入硬币运行**exp(−μ)的算法**，并返回结果。

1/(2m*(k + λ)) 或 1/((2m)*(k + λ)) 或 exp(−(k + λ)*ln(2m))

前面算法的扩展。这里，m是一个大于 0 的整数。

  1. 如果k > 0，则生成无偏随机比特，直到生成零比特或生成了k*m个比特，以先到者为准。如果生成了零比特，则返回 0。
  2. 创建一个输入硬币μ，执行以下操作：“翻转输入硬币，然后运行（稍后给出）**ln(1+y/z)的算法**，其中y/z = 1/1。如果调用和翻转都返回 1，则返回 1。否则，返回 0。”
  3. 使用μ输入硬币运行**exp(−μ)的算法**m次。如果任何一次调用返回 0，则返回 0。否则，返回 1。

1/(1+λ)

此算法是（Gonçalves 等人，2017）(39) 的双硬币伯努利工厂的一个特殊情况，并且对于所有λ参数都是均匀快速的。在本文档中，它将被调用**双硬币特殊情况**。(40)

  1. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则返回 1。
  2. 翻转输入硬币。如果它返回 1，则返回 0。否则，转到步骤 1。

λ/(1+λ)

返回 1 减去**1/(1+λ)的算法**的结果。

1/(2 − λ)

从双硬币特殊情况修改。

  1. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则返回 1。
  2. 翻转输入硬币。**如果它返回 0**，则返回 0。否则，转到步骤 1。

c * λ * β / (β * (c * λ + d * μ) − (β − 1) * (c + d))

这是（Gonçalves 等人，2017）(39) 和（Vats 等人 2020）(41) 的通用**双硬币算法**。它需要两个输入硬币，它们分别以概率λ或μ输出正面（1）。它还接受参数c和d，每个参数都大于等于 0，以及一个在 [0, 1] 区间的β，这是一个所谓的“门户”或早期拒绝参数（当β = 1 时，公式简化为c * λ / (c * λ + d * μ)）。在 Vats 等人 (2020)(41) 中，β, c, d, λ 和μ分别对应于“门户”算法中的β, cy, cx, py, 和px，或者分别对应于“翻转门户”算法中的β, c̃x, c̃y, p̃x, 和p̃y。

  1. 以β的概率，转到步骤 2。否则，返回 0。（例如，使用β的分子和分母调用ZeroOrOne，如果该调用返回 0，则返回 0，否则转到步骤 2。ZeroOrOne在我的随机采样方法的文章中有描述。）
  2. 以c / (c + d)的概率，翻转λ输入硬币。否则，翻转μ输入硬币。如果λ输入硬币返回 1，则返回 1。如果μ输入硬币返回 1，则返回 0。如果相应的硬币返回 0，则转到步骤 1。

c * λ / (c * λ + d) 或 (c/d) * λ / (1 + (c/d) * λ))

此算法，也称为**logistic Bernoulli 工厂** (Huber 2016)(35), (Morina 等人., 2019)(16)，是上面双硬币算法的一个特殊情况，但这次仅使用一个输入硬币。

  1. 以d / (c + d)的概率，返回 0。
  2. 翻转输入硬币。如果翻转返回 1，则返回 1。否则，转到步骤 1。

（注意，Huber [2016] 在泊松点过程中指定了这个伯努利工厂，这似乎平均需要更多的随机性。）

(d + λ) / c

此算法目前仅在d和c为整数且 0 ≤ d < c时有效。

  1. 随机均匀地生成 [0, c) 中的一个整数，称之为i。
  2. 如果i < d，则返回 1。如果i = d，则翻转输入硬币并返回结果。如果都不是，则返回 0。

d / (c + λ)

在此算法中，c和d必须是有理数，1 ≤ c，且 0 ≤ d ≤ c。另请参见连分数的算法。（例如，当d = 1 时，此算法可以模拟形式为 1 / z 的概率，其中z大于 0 且由整数部分（c）和可以通过伯努利工厂模拟的分数部分（λ）组成。）

  1. 以c / (1 + c)的概率，返回一个概率为 d/c 为 1，否则为 0 的数字。
  2. 翻转输入硬币。如果翻转返回 1，则返回 0。否则，转到步骤 1。

示例： **1 / (c + λ)**：运行d / (c + λ)的算法，其中d = 1。

注意：快速证明此算法有效：令x为所需概率。那么——
x = (c / (1 + c)) * (d/c) +
(1−c / (1 + c)) * (λ*0 + (1−λ)*x)，
解出x得到x=d/(c+λ)。

(d + μ) / (c + λ)

合并了前两个章节的算法。此算法目前仅在d和c为整数且 0 ≤ d < c时有效。

  1. 以c / (1 + c)的概率，执行以下操作
    1. 随机均匀地生成 [0, c) 中的一个整数，称之为i。
    2. 如果i < d，则返回 1。如果i = d，则翻转μ输入硬币并返回结果。如果都不是，则返回 0。
  2. 翻转λ输入硬币。如果翻转返回 1，则返回 0。否则，转到步骤 1。

(d + μ) / ((d + μ) + (c + λ))

在此算法中，c和d是大于等于 0 的整数，而λ和μ是两个不同输入硬币的正面概率。在此算法的预期用途中，λ和μ由两个均匀 PSRN 的小数部分支持，而c和d分别是它们的整数部分。

  1. 运行下面的子算法，使用μ输入硬币，并令a = d+c，b = 1+d+c。如果它返回 1
    1. 如果c为 0，则返回 1。
    2. 使用μ输入硬币运行子算法，并令a = d，b = d + c。如果它返回 1，则返回 1。否则，返回 0。
  2. 翻转λ输入硬币。如果翻转返回 1，则返回 0。否则，转到步骤 1。

下面的子算法模拟 (a+μ) / (b+μ)。

  1. 以b / (1+b)的概率，执行以下操作
    1. 随机均匀地生成 [0, b) 中的一个整数，称之为i。
    2. 如果i < a，则返回 1。如果i = a，则翻转输入硬币并返回结果。如果都不是，则返回 0。
  2. 翻转μ输入硬币。如果翻转返回 1，则返回 0。否则，转到步骤 1。

dk / (c + λ)k, 或 (d / (c + λ))k

在此算法中，c必须大于等于 1，d必须在 [0, c] 区间内，而k必须是一个大于等于 0 的整数。

  1. 将i设为 0。
  2. 如果k为 0，则返回 1。
  3. 以c / (1 + c)的概率，执行以下操作
    1. 以d/c的概率，将 1 加到i，然后要么返回 1（如果i现在大于等于k），要么中止这些子步骤并转到步骤 2，否则。
    2. 返回 0。
  4. 翻转输入硬币。如果翻转返回 1，则返回 0。否则，转到步骤 2。

1 / (1 + (c/d)*λ)

此算法是双硬币算法的一个特殊情况。在此算法中，c/d必须大于等于 0。

  1. 如果c为 0，则翻转μ输入硬币并返回结果。
  2. 以d/(c+d)的概率，翻转μ输入硬币并返回结果。
  3. 翻转输入硬币。如果翻转返回 1，则返回 0。否则，转到步骤 2。

示例： **μ / (1 + (c/d)*λ)**（需要两个模拟λ或μ的输入硬币）：运行 1 / (1 + (c/d)*λ)的**算法**，使用λ输入硬币。如果它返回 0，则返回 0。否则，翻转μ输入硬币并返回结果。

λ + μ

(Nacu 和 Peres 2005，命题 14(iii))(15)。此算法需要两个模拟λ或μ的输入硬币，以及一个参数ϵ，该参数必须大于 0 且选择得使得 λ + μ < 1 − ϵ。

  1. 创建一个ν输入硬币，执行以下操作：“生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则翻转λ输入硬币并返回结果。否则，翻转μ输入硬币并返回结果。”
  2. 调用**2014 算法**、**2016 算法**或**2019 算法**，使用ν输入硬币，x/y = 2/1，i = 1（对于 2019 算法），和ϵ = ϵ，并返回结果。

λ − μ

(Nacu 和 Peres 2005，命题 14(iii-iv))(15)。此算法需要两个模拟λ或μ的输入硬币，以及一个参数ϵ，该参数必须大于 0 且选择得使得 λ − μ > ϵ（并且应选择得使得ϵ略小于λ − μ）。

  1. 创建一个ν输入硬币，执行以下操作：“生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则翻转λ输入硬币并返回**1 减去结果**。否则，翻转μ输入硬币并返回结果。”
  2. 调用**2014 算法**、**2016 算法**或**2019 算法**，使用ν输入硬币，x/y = 2/1，i = 1（对于 2019 算法），和ϵ = ϵ，并返回 1 减去结果。

1 − λ

(Flajolet 等人，2010)(1)：翻转λ输入硬币，如果结果是 1，则返回 0，否则返回 1。

ν * λ + (1 − ν) * μ

(Flajolet 等人，2010)(1)：翻转ν输入硬币。如果结果是 0，则翻转λ输入硬币并返回结果。否则，翻转μ输入硬币并返回结果。

λ + μ − (λ * μ)

(Flajolet 等人，2010)(1)：翻转λ输入硬币和μ输入硬币。如果任一翻转返回 1，则返回 1；否则返回 0。

(λ + μ) / 2

(Nacu 和 Peres 2005，命题 14(iii))(15)；(Flajolet 等人，2010)(1)：生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则翻转λ输入硬币并返回结果。否则，翻转μ输入硬币并返回结果。

λx/y

在下面的算法中，x/y在开区间 (0, 1) 中的情况是 Mendo (2019)(22) 的最新工作。该算法仅在x/y大于等于 0 时有效。

  1. 如果x/y为 0，则返回 1。
  2. 如果x/y等于 1，则翻转输入硬币并返回结果。
  3. 如果x/y大于 1
    1. 将ipart设为 floor(x/y)，将fpart设为 `rem(x, y)`。
    2. 如果fpart大于 0，则从ipart中减去 1，然后递归调用此算法，令x = floor(fpart/2) 和y = y，然后再次递归调用此算法，令x = fpart − floor(fpart/2) 和y = y。如果任何一次调用返回 0，则返回 0。（这样做是为了避免调用此算法时小数部分非常接近 0，因为算法的运行速度比小数部分接近 1 时要慢得多。）
    3. 如果ipart大于等于 1，则将输入硬币翻转ipart次。如果任何一次翻转返回 1，则返回 0。
    4. 返回 1。
  4. x/y小于 1，因此将i设为 1。
  5. 翻转输入硬币；如果它返回 1，则返回 1。
  6. 以x/(y*i)的概率，返回 0。
  7. 将 1 加到i并转到步骤 5。

注意：当x/y小于 1 时，此算法平均需要的最小硬币翻转次数将随着λ接近 0 而无界增长。实际上，在没有关于λ的额外信息的情况下，任何快速的伯努利工厂算法都无法避免这种无界增长（Mendo 2019）(22)。

λμ

此算法基于前面的算法，但修改为接受第二个输入硬币（其正面概率为μ），而不是固定的指数值。据我所知，我不知道任何其他作者的文章或论文提出了这个特定的伯努利工厂。

  1. 将i设为 1。
  2. 翻转模拟基数的输入硬币λ；如果它返回 1，则返回 1。
  3. 翻转模拟指数的输入硬币μ；如果它返回 1，则以 1/i 的概率返回 0。
  4. 将 1 加到i并转到步骤 1。

sqrt(λ)

使用 λ1/2 的算法。

λ * μ

(Flajolet 等人，2010)(1)：翻转λ输入硬币和μ输入硬币。如果两次翻转都返回 1，则返回 1；否则返回 0。

λ * x/y (线性 Bernoulli 工厂)

一般来说，当x/y > 1 时，此函数将在开区间 (0, 1) 内触及 0 或 1。这使得在此情况下模拟该函数相对不平凡。

Huber 多年来为该函数提出了几种算法。

第一个算法在本文件中称为**2014 算法** (Huber 2014)(4)。它使用三个参数：x, y, 和ϵ，其中x/y > 0 且ϵ > 0。当x/y大于 1 时，必须选择ϵ参数，使得 λ * x/y < 1 − ϵ，以便将函数限制在 0 和 1 之外。因此，算法必须能够获得对λ的一些了解。（事实上，如模拟结果所示，ϵ的选择对该算法的性能至关重要；为获得最佳效果，应选择ϵ使得 λ * x/y 略小于 1 − ϵ。）下面描述的算法还包括某些 Huber 未提及的特殊情况，以使其更通用。

  1. 特殊情况：如果x为 0，则返回 0。否则，如果x等于y，则翻转输入硬币并返回结果。否则，如果x小于y，则：(a) 以 x/y 的概率翻转输入硬币并返回结果；否则 (b) 返回 0。
  2. 将c设为 x/y，将k设为 23 / (5 * ϵ)。
  3. 如果ϵ大于 644/1000，则将ϵ设为 644/1000。
  4. 将i设为 1。
  5. 翻转输入硬币。如果它返回 0，则生成数字，每个数字的概率为 (c − 1) / c 为 1，否则为 0，直到这样生成 0 为止，然后将这样生成的数字数量（包括最后一个）加到i。
  6. 将 1 从i中减去，然后如果i为 0，则返回 1。
  7. 如果i小于k，则转到步骤 5。
  8. 如果i等于k或更大
    1. 生成i个数字，每个数字的概率为 2 / (ϵ + 2) 为 1，否则为 0。如果其中任何一个数字是 0，则返回 0。
    2. 将c乘以 2 / (ϵ + 2)，将ϵ除以 2，并将k乘以 2。
  9. 如果i为 0，则返回 1。否则，转到步骤 5。

第二个算法称为**2016 算法** (Huber 2016)(35)，它使用相同的参数x, y, 和ϵ，并且其描述使用了相同的特殊情况。不同之处在于它涉及所谓的“logistic Bernoulli 工厂”，在本文件中将其替换为模拟相同函数的另一个工厂。当x/y大于 1 时，必须选择ϵ参数，使得 λ * x/y ≤ 1 − ϵ。

  1. 适用于 2014 算法的相同特殊情况。
  2. 将m设为 ceil(1 + 9 / (2 * ϵ))。
  3. 将β设为 1 + 1 / (m − 1)。
  4. **算法 A** 是 Huber 对此步骤的称呼。将s设为 1，然后当s大于 0 且小于m时
    1. 运行c/d = β * x/y 的**logistic Bernoulli 工厂**算法。
    2. 令s = s − z * 2 + 1，其中z是 logistic Bernoulli 工厂的结果。
  5. 如果s不等于 0，则返回 0。
  6. 以 1/β 的概率，返回 1。
  7. 对当前正在运行的算法进行单独运行，令x/y = β * x/y 和ϵ = 1 − β * (1 − ϵ)。如果单独运行返回 0，则返回 0。
  8. **高功率 logistic Bernoulli 工厂** 是 Huber 对此步骤的称呼。将s设为 1，然后当s大于 0 且小于等于m减 2 时
    1. 运行c/d = β * x/y 的**logistic Bernoulli 工厂**算法。
    2. 令s = s + z * 2 − 1，其中z是 logistic Bernoulli 工厂的结果。
  9. 如果s等于m减 1，则返回 1。
  10. 将 1 从m中减去，然后转到步骤 7。

提出 2016 算法的论文还包括第三种算法，如下所述，该算法仅在 λ * x / y 已知小于 1/2 时有效。第三种算法接受三个参数：x, y, 和m，并且必须选择m使得 λ * x / y ≤ m < 1/2。

  1. 适用于 2014 算法的相同特殊情况。
  2. 运行c/d = (x/y) / (1 − 2 * m) 的**logistic Bernoulli 工厂**算法。如果它返回 0，则返回 0。
  3. 以 1 − 2 * m 的概率，返回 1。
  4. 运行 2014 算法或 2016 算法，令x/y = (x/y) / (2 * m) 和ϵ = 1 − m。

注意：有关模拟 λ*(x/y) 的近似方法，请参见 "伯努利工厂算法的补充说明" 页面。

(λ * x/y)i

(Huber 2019)(42)。此算法在本文件中称为**2019 算法**，它使用四个参数：x, y, i, 和ϵ，其中x/y > 0，i ≥ 0 是一个整数，并且ϵ > 0。当x/y大于 1 时，必须选择ϵ参数，使得 λ * x/y < 1 − ϵ。它还具有 Huber 2019 中未提及的特殊情况。

  1. 特殊情况：如果i为 0，则返回 1。如果x为 0，则返回 0。否则，如果x等于y且i等于 1，则翻转输入硬币并返回结果。
  2. 特殊情况：如果x小于y且i = 1，则：(a) 以 x/y 的概率翻转输入硬币并返回结果；否则 (b) 返回 0。
  3. 特殊情况：如果x小于y，则创建一个辅助硬币μ，执行以下操作：“(a) 以 x/y 的概率翻转输入硬币并返回结果；否则 (b) 返回 0”，然后使用此辅助硬币运行**（μi/1)的算法**（前面给出）。
  4. 将t设为 355/100，将c设为 x/y。
  5. 如果i为 0，则返回 1。
  6. 当i = t / ϵ时
    1. 令β为 (1 − ϵ / 2) / (1 − ϵ)。
    2. 运行**（1/β）i的算法**（在无理数常数部分给出，标记为**xy**）。如果它返回 0，则返回 0。
    3. 将c乘以β，然后将ϵ除以 2。
  7. 运行**logistic Bernoulli 工厂**，其中c/d = c，然后将结果设为z。令i = i + 1 − z * 2，然后转到步骤 5。

ϵ / λ

(Lee 等人 2014)(43) 此算法除了输入硬币外，还有一个参数ϵ，它必须大于 0 且选择得使得 ϵ 小于λ。

  1. 将β设为 max(ϵ, 1/2)，将γ设为 1 − (1 − β) / (1 − (β / 2))。
  2. 创建一个μ输入硬币，该硬币翻转输入硬币并返回 1 减去结果。
  3. 以ϵ的概率，返回 1。
  4. 调用**2014 算法**、**2016 算法**或**2019 算法**，使用μ输入硬币，x/y = 1 / (1 − ϵ)，i = 1（对于 2019 算法），和ϵ = γ。如果结果为 0，则返回 0。否则，转到步骤 3。（注意，以这种方式运行算法模拟的概率是 (λ − ϵ)/(1 − ϵ) 或 1 − (1 − λ)/(1 − ϵ)）。

arctan(λ) /λ

(Flajolet 等人，2010)(1)

  1. 生成一个均匀(0, 1)随机变量u。
  2. 对数字u采样两次，并翻转输入硬币两次。如果任何一个调用或翻转返回 0，则返回 1。
  3. 对数字u采样两次，并翻转输入硬币两次。如果任何一个调用或翻转返回 0，则返回 0。否则，转到步骤 2。

注意到 Flajolet 论文中使用的偶数校验构造等价于双硬币特殊情况，对于每个λ参数都是均匀快速的，因此可以通过以下方式使上述算法均匀快速

  1. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则返回 1。
  2. 生成一个均匀(0, 1)随机变量u，如果它尚未生成。
  3. 对数字u采样两次，并翻转输入硬币两次。如果所有这些调用和翻转都返回 1，则返回 0。否则，转到步骤 1。

arctan(λ)

(Flajolet 等人，2010)(1)：调用 **arctan(λ) /λ 的算法** 并翻转输入硬币。如果调用和翻转都返回 1，则返回 1；否则返回 0。

cos(λ)

此算法改编自通用鞅算法，用于此函数的级数展开。事实上，这是（Łatuszyński 等人 2009/2011）(23) 的算法 3 的一个特殊情况（该算法比命题 3.4，即通用鞅算法更通用）。cos(λ) 的级数展开为 1 − λ2/(2!) + λ4/(4!) − ...，这是一个交错级数，但指数每项增加 2（而不是 1）。因此，系数为 1, 1/(2!), 1/(4!), ...。下截断级数是结尾为负项的截断级数，而相应的上截断是相同截断但没有最后一项减号。此级数展开满足算法 3 的要求，因为概率为 1——

  • 下截断小于或等于其相应的上截断，
  • 下截断和上截断在 [0, 1] 区间内，
  • 每个下截断大于或等于前一个下截断，
  • 每个上截断小于或等于前一个上截断，并且
  • 下截断和上截断的期望值从下方和上方接近λ。

模拟 cos(λ) 的算法如下。

  1. 将u设为 1，将w设为 1，将ℓ设为 0，将n设为 1，并将fac设为 2。
  2. 生成一个 uniform(0, 1) 随机变量，称之为 ret。
  3. 如果w不为 0，则翻转输入硬币。如果翻转返回 0，则将w设为 0。再次执行此步骤。（注意，在通用鞅算法中，此步骤只翻转一个硬币。这里最多翻转两个硬币，因为指数增加了 2 而不是 1。）
  4. 如果n是偶数，则将u设为 ℓ + w / fac。否则，将ℓ设为u − w / fac。（此处我们除以 2 乘以 n 的阶乘。）
  5. 如果 ret 小于（或等于）ℓ，则返回 1。如果 ret 小于 u，则转到下一步。如果都不是，则返回 0。（如果 ret 是均匀 PSRN，这些比较应通过 **URandLessThanReal 算法**进行，该算法在我关于 PSRN 的文章中有所描述。）
  6. 将 1 加到n，然后将fac乘以 (n * 2 − 1) * (n * 2)，然后转到步骤 3。

sin(λ)

此算法同样是（Łatuszyński 等人 2009/2011）(23) 的算法 3 的一个特殊情况。sin(λ) 可以重写为 λ * (1 − λ2/(3!) + λ4/(5!) − ...)，它包含一个交错级数，其中指数每项增加 2（而不是 1）。因此，系数为 1, 1/(3!), 1/(5!), ...。这个级数展开出于与 cos(λ) 级数相同的理由，满足算法 3 的要求。

模拟 sin(λ) 的算法如下。

  1. 翻转输入硬币。如果它返回 0，则返回 0。
  2. 将u设为 1，将w设为 1，将ℓ设为 0，将n设为 1，并将fac设为 6。
  3. 生成一个 uniform(0, 1) 随机变量，称之为 ret。
  4. 如果w不为 0，则翻转输入硬币。如果翻转返回 0，则将w设为 0。再次执行此步骤。
  5. 如果n是偶数，则将u设为 ℓ + w / fac。否则，将ℓ设为u − w / fac。
  6. 如果 ret 小于（或等于）ℓ，则返回 1。如果 ret 小于 u，则转到下一步。如果都不是，则返回 0。（如果 ret 是均匀 PSRN，这些比较应通过 **URandLessThanReal 算法**进行，该算法在我关于 PSRN 的文章中有所描述。）
  7. 将 1 加到n，然后将fac乘以 (n * 2) * (n * 2 + 1)，然后转到步骤 4。

(1−λ)/cos(λ)

(Flajolet 等人，2010)(1)

  1. 翻转输入硬币，直到翻转返回 0。然后将返回 1 的次数设为G。
  2. 如果G是**奇数**，则返回 0。
  3. 生成一个均匀(0, 1)随机变量U，然后将i设为 1。
  4. 当i小于G时
    1. 生成一个均匀(0, 1)随机变量V。
    2. 如果i是奇数且V小于U，则返回 0。
    3. 如果i是偶数且U小于V，则返回 0。
    4. 将 1 加到i，然后将U设为V。
  5. 返回 1。

(1−λ) * tan(λ)

(Flajolet 等人，2010)(1)

  1. 翻转输入硬币，直到翻转返回 0。然后将返回 1 的次数设为G。
  2. 如果G是**偶数**，则返回 0。
  3. 生成一个均匀(0, 1)随机变量U，然后将i设为 1。
  4. 当i小于G时
    1. 生成一个均匀(0, 1)随机变量V。
    2. 如果i是奇数且V小于U，则返回 0。
    3. 如果i是偶数且U小于V，则返回 0。
    4. 将 1 加到i，然后将U设为V。
  5. 返回 1。

ln(1+λ)

(Flajolet 等人，2010)(1)

  1. 生成一个均匀(0, 1)随机变量u。
  2. 翻转输入硬币。如果它返回 0，则再次翻转硬币并返回结果。
  3. 对数字u采样。如果结果为 0，则翻转输入硬币并返回结果。
  4. 翻转输入硬币。如果它返回 0，则返回 0。
  5. 对数字u采样。如果结果为 0，则返回 0。否则，转到步骤 2。

注意到 Flajolet 论文中使用的偶数校验构造等价于双硬币特殊情况，对于每个λ参数都是均匀快速的，因此可以通过以下方式使上述算法均匀快速

  1. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则翻转输入硬币并返回结果。
  2. 生成一个均匀(0, 1)随机变量u，如果u尚未生成。
  3. 对数字u采样，然后翻转输入硬币。如果调用和翻转都返回 1，则返回 0。否则，转到步骤 1。

1 − ln(1+λ)

返回 1 减去**ln(1+λ)的算法**的结果。(44)

arcsin(λ) + sqrt(1 − λ2) − 1

(Flajolet 等人，2010)(1)。此处给出的算法使用双硬币特殊情况而不是偶数校验构造。

  1. 生成一个均匀(0, 1)随机变量u。
  2. 创建一个辅助硬币μ，执行以下操作：“对数字u采样两次，并翻转输入硬币两次。如果所有这些调用和翻转都返回 1，则返回 0。否则，返回 1。”
  3. 使用辅助硬币μ调用 **μ1/2 的算法**。如果它返回 0，则返回 0。
  4. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则翻转输入硬币并返回结果。
  5. 对数字u采样一次，并翻转输入硬币一次。如果调用和翻转都返回 1，则返回 0。否则，转到步骤 4。

arcsin(λ) / 2

Flajolet 论文没有详细解释 arcsin(λ)/2 如何通过伯努利工厂构造从 arcsin(λ) + sqrt(1 − λ2) − 1 产生，但这里有一个算法。(45) 但是，预期的输入硬币翻转次数会随着λ接近 1 而无界增长。

  1. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则运行**arcsin(λ) + sqrt(1 − λ2) − 1 的算法**并返回结果。
  2. 创建一个辅助硬币μ，执行以下操作：“翻转输入硬币两次。如果两次翻转都返回 1，则返回 0。否则，返回 1。”（该硬币模拟 1 − λ2。）
  3. 使用辅助硬币μ调用 **μ1/2 的算法**。如果它返回 0，则返回 1；否则返回 0。（此步骤有效地消除了 sqrt(1 − λ2) − 1 部分并将结果除以 2。）

涉及多对数表达式

以下算法模拟表达式 Lir(λ) * (1 / λ − 1)，其中r是一个大于等于 1 的整数。然而，即使对于相对较小的r（如 6），该表达式也很快接近一条直线。

如果λ为 1/2，则此表达式简化为 Lir(1/2)。另请参见 (Flajolet 等人，2010)(1)。另请参见 "凸组合"（1/2 的情况通过将形成多对数常数的级数分解为g(i) = (1/2)i（总和为 1）和hi() = 1/ir（其中i ≥ 1）来工作）。

  1. 翻转输入硬币，直到它返回 0，令t为 1 加上硬币返回 1 的次数。
  2. 返回一个概率为 1/tr 为 1，否则为 0 的数字。

特定常量的算法

本节展示了模拟特定种类无理数概率的算法。

1 / φ （1 除以黄金比例）

此算法使用连分数部分所述的算法来模拟 1 除以黄金比例，其连分数的局部分母为 1, 1, 1, 1, ...。

  1. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则返回 1。
  2. 对当前正在运行的算法进行单独运行。如果单独运行返回 1，则返回 0。否则，转到步骤 1。

另一种模拟此概率的算法，在法语中（Penaud 和 Roques 2002）(46) 中出现，它是迭代的但需要算术运算。

sqrt(2) − 1

另一个连分数的例子是 2 的平方根的小数部分，其中局部分母为 2, 2, 2, 2, ...。模拟此数字的算法如下

  1. 以 2/3 的概率，生成一个无偏随机比特并返回该比特。
  2. 对当前正在运行的算法进行单独运行。如果单独运行返回 1，则返回 0。否则，转到步骤 1。

1/sqrt(2)

第三个连分数示例展示了如何模拟概率 1/z，其中z > 1 具有已知的简单连分数展开。在这种情况下，局部分母如下：floor(z), a[1], a[2], ..., 其中a[i]是z的局部分母（不包括z的整数部分）。在 1/sqrt(2) 的例子中，局部分母为 1, 2, 2, 2, ..., 其中 1 排在最前面，因为 floor(sqrt(2)) = 1。模拟 1/sqrt(2) 的算法如下

该算法从等于 1 的pos开始。然后执行以下步骤。

  1. 如果pos为 1，则以 1/2 的概率返回 1。如果pos大于 1，则以 2/3 的概率，生成一个无偏随机比特并返回该比特。
  2. 对当前正在运行的算法进行单独运行，但将pos设为pos + 1。如果单独运行返回 1，则返回 0。否则，转到步骤 1。

tanh(1/2) 或 (exp(1) − 1) / (exp(1) + 1)

算法从等于 2 的k开始。然后执行以下步骤。

  1. 以 k/(1+k) 的概率，返回一个概率为 1/k 为 1，否则为 0 的数字。
  2. 对当前正在运行的算法进行单独运行，但将k设为k + 4。如果单独运行返回 1，则返回 0。否则，转到步骤 1。

arctan(x/y) * y/x

(Flajolet 等人，2010)(1)

  1. 生成一个均匀(0, 1)随机变量u。
  2. 生成一个概率为 x * x/(y * y) 为 1，否则为 0 的数字。如果该数字为 0，则返回 1。
  3. 对数字u采样两次。如果这两个调用中的任何一个返回 0，则返回 1。
  4. 生成一个概率为 x * x/(y * y) 为 1，否则为 0 的数字。如果该数字为 0，则返回 0。
  5. 对数字u采样两次。如果这两个调用中的任何一个返回 0，则返回 0。否则，转到步骤 2。

注意到 Flajolet 论文中使用的偶数校验构造等价于双硬币特殊情况，它是均匀快速的，因此可以通过以下方式使上述算法均匀快速

  1. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则返回 1。
  2. 生成一个均匀(0, 1)随机变量u，如果它尚未生成。
  3. 以 x * x/(y * y) 的概率，对数字u采样两次。如果这两个调用都返回 1，则返回 0。
  4. 转到步骤 1。

π / 12

两种算法

  • 第一种算法：使用**arcsin(1/2) / 2 的算法**。在算法指示“翻转输入硬币”的地方，而是生成一个无偏随机比特。
  • 第二种算法：以 2/3 的概率，返回 0。否则，运行**π / 4 的算法**并返回结果。

π / 4

(Flajolet 等人，2010)(1)

  1. 随机生成一个 [0, 6) 区间中的整数，称之为n。
  2. 如果n小于 3，则返回**arctan(1/2) * 2 的算法**的结果。否则，如果n等于 3，则返回 0。否则，返回**arctan(1/3) * 3 的算法**的结果。

1 / π

(Flajolet 等人，2010)(1)

  1. 将t设为 0。
  2. 以 1/4 的概率，将 1 加到t并重复此步骤。否则，转到步骤 3。
  3. 以 1/4 的概率，将 1 加到t并重复此步骤。否则，转到步骤 4。
  4. 以 5/9 的概率，将 1 加到t。
  5. 生成 2*t 个无偏随机比特（即，0 或 1，以相等的概率选择），并且如果生成的零比特多于一个比特或一个比特多于零比特，则返回 0。（注意，这个条件可以在所有比特生成之前检查。）重复此步骤两次。
  6. 返回 1。

有关此算法如何推导的草图，请参见附录。

(a/b)x/y

在下面的算法中，a, b, x, 和y是整数，并且x/y在开区间 (0, 1) 中的情况是 Mendo (2019)(22) 的最新工作。该算法仅在以下情况下有效——

  • x/y 大于等于 0 且 a/b 在区间 [0, 1] 内，或者
  • x/y 小于 0 且 a/b 大于等于 1。

算法如下。

  1. 如果 x/y 小于 0，交换 a 和 b，并去掉 x/y 的符号。如果此时 a/b 不再位于区间 [0, 1] 内，则返回错误。
  2. 如果 x/y 等于 1，则以 a/b 的概率返回 1，否则返回 0。
  3. 如果 x 为 0，则返回 1。否则，如果 a 为 0，则返回 0。否则，如果 a 等于 b，则返回 1。
  4. 如果x/y大于 1
    1. 将ipart设为 floor(x/y)，将fpart设为 `rem(x, y)`。
    2. 如果fpart大于 0，则从ipart中减去 1，然后递归调用此算法，令x = floor(fpart/2) 和y = y，然后再次递归调用此算法，令x = fpart − floor(fpart/2) 和y = y。如果任何一次调用返回 0，则返回 0。（这样做是为了避免调用此算法时小数部分非常接近 0，因为算法的运行速度比小数部分接近 1 时要慢得多。）
    3. 如果 ipart 大于等于 1，则随机生成一个数，该数为 1 的概率为 aipart/bipart，否则为 0。（或者，随机生成 ipart 个数，每个数以 a/b 的概率为 1 或 0，然后将它们全部相乘得到一个数。）如果该数为 0，则返回 0。
    4. 返回 1。
  5. 将i设为 1。
  6. 以 a/b 的概率返回 1。
  7. 否则，以 x/(y*i) 的概率返回 0。
  8. 将 i 加 1 并转到步骤 6。

exp(−x/y)

该算法接受非负整数 x ≥ 0 和正整数 y > 0，并以 exp(-x/y) 的概率输出 1，否则输出 0。它源自 (Canonne et al. 2020)(47)。

  1. 特殊情况：如果 x 为 0，则返回 1。（这是因为概率变为 exp(0) = 1。）
  2. 如果 x > y（因此 x/y 大于 1），则使用 x = y = 1 调用此算法（递归地）floor(x/y) 次，并使用 x = x − floor(x/y) * y 和 y = y 调用一次。如果所有这些调用都返回 1，则返回 1；否则，返回 0。
  3. 将 r 设置为 1，将 i 设置为 1。
  4. 以 (y * i − x) / (y * i) 的概率返回 r。
  5. 将 r 设置为 1 − r，将 i 加 1，然后转到步骤 4。

exp(−z)

该算法与之前的算法类似，不同之处在于指数 z 可以是任何大于等于 0 的实数，只要 z 可以写成若干个分量之和，而每个分量的小数部分都可以通过一个伯努利工厂算法来模拟，该算法以等于该小数部分的概率输出正面。（这利用了恒等式 exp(−a) = exp(−b) * exp(−c)。）

更具体地说，

  1. 将 z 分解为 n > 0 个大于 0 的分量之和。例如，如果 z = 3.5，它只能分解为一个分量 3.5（其小数部分易于模拟）；如果 z = π，它可以分解为四个分量，每个分量都是 (π / 4)（其模拟如本页前面所述，并不简单）。
  2. 对于通过此方式找到的每个分量 LC[i]，令 LI[i] = floor(LC[i])，令 LF[i] = LC[i] − floor(LC[i]）（即 LC[i] 的小数部分）。

然后算法如下：

  • 对于每个分量 LC[i]，调用 **exp(− LI[i]/1) 的算法**，并调用适用于 **exp(−λ)** 的 **通用鞅算法**，使用模拟 LF[i] 的输入硬币。如果其中任何一个调用返回 0，则返回 0；否则，返回 1。（另请参见 (Canonne et al. 2020)(47)。）

(a/b)z

该算法与之前的幂运算算法类似，不同之处在于指数 z 可以是任何大于等于 0 的实数，只要 z 可以写成若干个分量之和，而每个分量的小数部分都可以通过一个伯努利工厂算法来模拟，该算法以等于该小数部分的概率输出正面。此算法利用了与 exp 类似的恒等式，并且仅在 z ≥ 0 且 a/b 在区间 [0, 1] 内时才有效。

像 **exp(− z) 算法** 一样，将 z 分解为 LC[i], LI[i], 和 LF[i]。然后算法如下。

  • 如果 z 为 0，则返回 1。否则，如果 a 为 0，则返回 0。否则，对于每个分量 LC[i]（直到算法返回一个数字）
    1. 调用 **(a/b)LI[i]/1 的算法**。如果它返回 0，则返回 0。
    2. 将 j 设置为 1。
    3. 随机生成一个数，该数为 1 的概率为 a/b，否则为 0。如果该数为 1，则中止这些步骤并继续处理下一个分量，或者，如果没有更多分量，则返回 1。
    4. 翻转模拟 LF[i]（即指数）的输入硬币；如果它返回 1，则以 1/j 的概率返回 0。
    5. 将 j 加 1 并转到子步骤 2。

1 / (1 + exp(x / (y * 2prec)) (LogisticExp)

这是指数随机变量（速率为 x/y）的 prec 位（小数点后第 prec 位）设置为 1 的概率。该算法是 **logistic Bernoulli factory** 的特例。

  1. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则返回 1。
  2. 调用 **exp(− x/(y * 2prec)) 的算法**。如果调用返回 1，则返回 1。否则，转到步骤 1。

1 / (1 + exp(z / 2prec)) (LogisticExp)

这与之前的算法类似，不同之处在于 z 可以是 **exp(−z) 算法** 中描述的任何实数。

像 **exp(−z) 算法** 一样，将 z 分解为 LC[i], LI[i], 和 LF[i]。然后算法如下。

  1. 对于每个分量 LC[i]，创建一个输入硬币，该硬币执行以下操作：（a）以 1/(2prec) 的概率返回 1，如果模拟 LF[i] 的输入硬币返回 1；（b）返回 0。
  2. 以 1/2 的概率返回 0。
  3. 使用 x = ∑i LI[i] 和 y = 2prec 调用 **exp(− x/y) 的算法**。如果此调用返回 0，则转到步骤 2。
  4. 对于每个分量 LC[i]，调用 **exp(−λ) 的算法**，使用步骤 1 中创建的相应输入硬币。如果其中任何一个调用返回 0，则转到步骤 2。否则，返回 1。

ζ(3) * 3 / 4 和其他与 Zeta 相关的常数

(Flajolet et al., 2010)(1)。它可以看作是函数 1/(1 + a * b * c) 的三重积分，其中 a、b 和 c 是 uniform(0, 1) 的随机变量。该算法如下给出，但使用双硬币特殊情况而不是偶校验构造。请注意，论文第 5 节中的三重积分是 ζ(3) * 3 / 4，而不是 ζ(3) * 7 / 8。（此处，ζ(x) 是黎曼 zeta 函数。）

  1. 生成三个 uniform(0, 1) 的随机变量。
  2. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则返回 1。
  3. 对步骤 1 中生成的三个数字中的每一个进行采样。如果所有三个调用都返回 1，则返回 0。否则，转到步骤 2。（这实现了涉及均匀随机变量的三重积分。）

这可以扩展到覆盖任何形式为 ζ(k) * (1 − 2−(k − 1)) 的常数，其中 k ≥ 2 是整数，正如 Flajolet 论文在提到 ζ(5) * 31 / 32 时（可能应为 ζ(5) * 15 / 16）所暗示的，使用以下算法。

  1. 生成 k 个 uniform(0, 1) 的随机变量。
  2. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则返回 1。
  3. 对步骤 1 中生成的 k 个数字中的每一个进行采样。如果所有 k 个调用都返回 1，则返回 0。否则，转到步骤 2。

erf(x)/erf(1)

在以下算法中，x 是区间 [0, 1] 内的一个实数。

  1. 生成一个 uniform(0, 1) 随机变量，称之为 ret。
  2. 将 u 设置为指向与 ret 相同的值，并将 k 设置为 1。
  3. (在此及下一步中，我们创建 v，它是两个 uniform [0, 1] 随机变量的最大值。) 生成两个 uniform(0, 1) 随机变量，称为 a 和 b。
  4. 如果 a 小于 b，则将 v 设置为 b。否则，将 v 设置为 a。
  5. 如果 v 小于 u，则将 u 设置为 v，然后将 k 加 1，然后转到步骤 3。
  6. 如果 k 是奇数，如果 ret 小于 x，则返回 1，否则返回 0。（如果 ret 实现为均匀 PSRN，则此比较应通过 **URandLessThanReal 算法** 完成，该算法在我的 关于 PSRN 的文章 中有描述。）
  7. 转到步骤 1。

实际上，该算法利用了与 Forsythe 随机抽样方法相关的定理 (Forsythe 1972)(48)。有关更多信息，请参阅附录中的“由某些排列引起的概率”部分。

注意：如果算法的最后一步是“返回 0”而不是“转到步骤 1”，那么该算法模拟的概率是 erf(x)*sqrt(π)/2。

2 / (1 + exp(2)) 或 (1 + exp(0)) / (1 + exp(1))

该算法利用了附录中“由某些排列引起的概率”部分提到的公式 2。这里，相关概率被重写为 1 − (∫(−∞, 1) (1 − exp(−max(0, min(1, z)))) * exp(−z) dz) / (∫(−∞, ∞) (1 − exp(−max(0, min(1, z))) * exp(−z) dz)。

  1. 生成一个 **指数** 随机变量 ex，然后将 k 设置为 1。
  2. 将 u 设置为指向与 ex 相同的值。
  3. 生成一个 **uniform(0,1)** 随机变量 v。
  4. 如果 u 小于 v，则将 stop 设置为 1，否则设置为 0。
  5. 如果 stop 为 1 且 k **为偶数**，则返回一个数：如果 ex **小于 1**，则为 0，否则为 1。否则，如果 stop 为 1，则转到步骤 1。
  6. 将 u 设置为 v，然后将 k 加 1，然后转到步骤 3。

(1 + exp(1)) / (1 + exp(2))

该算法利用了附录中“由某些排列引起的概率”部分提到的定理。这里，相关概率被重写为 1 − (∫(−∞, 1/2) exp(−max(0, min(1, z))) * exp(−z) dz) / (∫(−∞, ∞) exp(−max(0, min(1, z)) * exp(−z) dz)。

  1. 生成一个 **指数** 随机变量 ex，然后将 k 设置为 1。
  2. 将 u 设置为指向与 ex 相同的值。
  3. 生成一个 **uniform(0,1)** 随机变量 v。
  4. 如果 u 小于 v，则将 stop 设置为 1，否则设置为 0。
  5. 如果 stop 为 1 且 k **为奇数**，则返回一个数：如果 ex **小于 1/2**，则为 0，否则为 1。否则，如果 stop 为 1，则转到步骤 1。
  6. 将 u 设置为 v，然后将 k 加 1，然后转到步骤 3。

(1 + exp(k)) / (1 + exp(k + 1))

该算法通过计算 exp(1) 的下界和上界来模拟此概率，这些界限随着计算的数字位数的增加而提高。这些界限是通过 Citterio 和 Pavani (2016)(49) 的算法计算的。注意此算法使用了 (Łatuszyński et al. 2009/2011, algorithm 2)(23) 中的方法。在此算法中，k 必须是大于等于 0 的整数。

  1. 如果 k 为 0，则运行 **2 / (1 + exp(2)) 的算法** 并返回结果。如果 k 为 1，则运行 **(1 + exp(1)) / (1 + exp(2)) 的算法** 并返回结果。
  2. 生成一个 uniform(0, 1) 随机变量，称之为 ret。
  3. 如果 k 大于等于 3，如果 ret 大于 38/100，则返回 0，如果 ret 小于 36/100，则返回 1。（这是一个提前返回步骤。如果 ret 实现为均匀 PSRN，则这些比较应通过 **URandLessThanReal 算法** 完成，该算法在我的 关于 PSRN 的文章 中有描述。）
  4. 将 d 设置为 2。
  5. 使用 Citterio 和 Pavani 算法计算 exp(1) 的下界和上界（分别为 LB 和 UB），形式为分子最多有 d 位数字的有理数。有关详细信息，请参阅附录。
  6. 将 rl 设置为 (1+LBk) / (1+UBk + 1)，将 ru 设置为 (1+UBk) / (1+LBk + 1)；这两个数字都应使用有理数算术计算。
  7. 如果 ret 大于 ru，则返回 0。如果 ret 小于 rl，则返回 1。（如果 ret 实现为均匀 PSRN，则这些比较应通过 **URandLessThanReal** 完成。）
  8. 将 d 加 1 并转到步骤 5。

欧拉-马斯克若尼常数 γ

模拟欧拉-马斯刻若尼常数 γ 的算法是 Mendo (2020)(30) 提出的。这解决了 (Flajolet et al., 2010)(1) 中提出的一个开放性问题。使用的级数由 Sondow (2005)(50) 给出。尽管 γ 是否为无理数尚不确定，但这里仍然给出了一个 γ 的算法。

  1. 将 ϵ 设置为 1，然后将 n、lamunq、lam、s、k 和 prev 全部设置为 0。
  2. 将 1 加到k，然后将s/(2k) 加到lam。
  3. 如果lamunq+ϵ ≤ lam + 1/(2k)，则转到步骤 8。
  4. 如果lamunq > lam + 1/(2k)，则转到步骤 8。
  5. 如果lamunq > lam + 1/(2k+1) 且lamunq+ϵ < 3/(2k+1)，则转到步骤 8。
  6. (此步骤将 γ 级数的一项加到 lamunq，并将 ϵ 设置为如果级数在此项和之前的项之后截断所产生的误差的上限。) 如果 n 为 0，则将 1/2 加到 lamunq 并将 ϵ 设置为 1/2。否则，将 B(n)/(2*n*(2*n+1)*(2*n+2)) 加到 lamunq 并将 ϵ 设置为 min(prev, (2+B(n)+(1/n))/(16*n*n))，其中 B(n) 是存储 n 所需的最小位数（或 b≥1 且 n < 2b 的最小整数 b）。
  7. 将 n 加 1，然后将 prev 设置为 ϵ，然后转到步骤 3。
  8. 令bound为lam+1/(2k)。如果lamunq+ϵ ≤ bound，则将s设为 0。否则，如果lamunq > bound，则将s设为 2。否则，将s设为 1。
  9. 生成一个无偏随机比特。如果该比特是 1（这种情况发生的概率为 1/2），则转到步骤 2。否则，返回一个数字：如果s为 0，则返回 0；如果s为 2，则返回 1；否则，返回一个无偏随机比特（0 或 1 的概率相等）。

exp(−x/y) * z/t

该算法同样基于 Mendo (2020)(30) 的算法。在此算法中，x、y、z 和 t 是大于 0 的整数，其中 x 和/或 z 可以为 0，并且必须满足 exp(−x/y) * z/t 在区间 [0, 1] 内。

  1. 如果 z 为 0，则返回 0。如果 x 为 0，则返回一个数：以 z/t 的概率为 1，否则为 0。
  2. 将ϵ设为 1，然后将n, lamunq, lam, s和k全部设为 0。
  3. 将 1 加到k，然后将s/(2k) 加到lam。
  4. 如果 lamunq+ϵ ≤ lam + 1/(2k)，则转到步骤 9。
  5. 如果 lamunq > lam + 1/(2k)，则转到步骤 9。
  6. 如果lamunq > lam + 1/(2k+1) 且lamunq+ϵ < 3/(2k+1)，则转到步骤 8。
  7. (此步骤将 exp(−x/y) 交错级数的两项乘以 z/t 加到 lamunq，并将 ϵ 设置为当前总和与期望概率的接近程度的上限。) 令 m = n*2。设置 ϵ = z*xm/(t*(m!)*ym)。如果 m 为 0，则将 z*(y−x)/(t*y) 加到 lamunq。否则，将 z*xm*(m*y−x+y) / (t*ym+1*((m+1)!)) 加到 lamunq。
  8. 将 1 加到n并转到步骤 4。
  9. 令bound为lam+1/(2k)。如果lamunq+ϵ ≤ bound，则将s设为 0。否则，如果lamunq > bound，则将s设为 2。否则，将s设为 1。
  10. 生成一个无偏随机比特。如果该比特为 1（以 1/2 的概率发生），则转到步骤 3。否则，返回一个数：如果 s 为 0，则为 0；如果 s 为 2，则为 1；否则为无偏随机比特（0 或 1 的概率相等）。

ln(1+y/z)

另请参阅前面给出的 ln(1+λ) 算法。在此算法中，y/z 是区间 [0, 1] 内的有理数。（因此，当 y/z = 1/1 时，得到特殊情况 ln(2)。）

  1. 如果 y/z 为 0，则返回 0。
  2. 生成一个无偏随机比特。如果该比特为 1（以 1/2 的概率发生），则返回一个数：以 y/z 的概率为 1，否则为 0。
  3. 生成一个均匀(0, 1)随机变量u，如果u尚未生成。
  4. 对数字 u 进行采样，然后生成一个数：以 y/z 的概率为 1，否则为 0。如果调用返回 1 且生成的数字为 1，则返回 0。否则，转到步骤 2。

请求和公开问题

  1. 除了本页上的算法之外，还有哪些“相对简单”且以 0 和 1 之间的无理数概率成功的模拟算法？对于“相对简单”的伯努利工厂算法呢？这里的“相对简单”意味着算法

    • 应仅使用均匀随机整数（或比特）和整数算术。
    • 不使用浮点算术，也不直接使用平方根或超越函数。
    • 不直接计算基-n 展开。
    • 除非万不得已，否则不使用有理数算术或越来越复杂的近似。

    另请参阅 Flajolet et al. (2010)(1)。有许多方法可以描述无理数概率或工厂函数。我正在寻找描述无理数常数或工厂函数的论文或书籍，方法有以下几种：

    • 对于无理数常数
      • 简单的 连分数 展开。
      • 单位正方形内的封闭形状，其面积为无理数。（包括判断一个盒子位于形状内部、外部，或部分内部/外部的算法。） 示例。
      • 在封闭形状内生成一个 uniform (x, y) 点，然后以 x 的概率返回 1。对于哪些形状，x 的期望值是无理数？ 示例。
      • 将 [0, 1] 映射到 [0, 1] 的函数，其积分（曲线下的面积）为无理数。
    • 对于伯努利工厂函数
      • 具有以下任一级数展开的函数，仅使用有理数算术：
        • 幂级数，其中 f(0) 为 0 且 f(1) 为有理数，反之亦然（参见“某些幂级数”）。
        • 项非负的级数，可以“嵌入”离散概率质量函数之下（参见“凸组合”）。
        • 交错幂级数，其系数都在 [0, 1] 区间内，并且构成一个非递增序列（参见“某些幂级数”）。
        • 项非负且截断误差有界的级数（参见“某些收敛级数”）。
      • 一种计算两种多项式序列的方法，这些序列以 Bernstein 形式表示，并从上方和下方收敛到工厂函数，如下所示：（a）每个序列的多项式系数必须位于 [0, 1] 内，并且度数递增；（b）度为 n 的多项式的系数必须位于先前上层多项式和先前下层多项式（一旦多项式提升到度 n）的系数之上或“内部”。有关这些多项式的正式说明，请参见我在 Mathematics Stack Exchange 上的问题。

        注意：“补充说明”包括计算这些多项式的公式，适用于大类工厂函数，但没有一个能确保一般的期望硬币翻转次数有限，并且有人怀疑除非工厂函数是 C2 连续的（有两个或多个连续的“斜率”函数），否则有限次翻转是不可能的。因此，一个问题是：给定一个 C2 连续的工厂函数，是否有实际算法来构建收敛于该函数的多项式，其方式对伯努利工厂问题是必需的，并且期望的翻转次数是有限的（除了本文或补充说明中的算法之外）？

  2. 给定一个排列类（例如降序排列）和两个连续概率分布 D 和 E。考虑以下算法：生成一个独立随机变量序列（第一个分布为 D，其余为 E），直到序列不再遵循排列类，然后返回 n，即生成的数字个数减 1。在这种情况下

    1. 返回 n 的概率是多少？
    2. n 为奇数或偶数，或属于某类数字的概率是多少？
    3. 给定 n 为奇数或 n 为偶数，第一个生成数字的分布函数（CDF）是多少？

    显然，这些答案取决于特定的排列类和/或分布 D 和 E。因此，只适用于特定类和/或分布的答案是受欢迎的。另请参阅我的 Stack Exchange 问题 由排列引起的概率。

  3. 是否有更简单或更快的方法来实现二项系数的底 2 或自然对数？参见“某些收敛级数”部分中的示例。

  4. Łatuszyński et al. (2009/2011)(23) 的逆时间鞅算法（参见“通用工厂函数”）的一部分是为了模拟工厂函数 f(λ)，如下所示。对于从 1 开始的每个 n

    1. 翻转输入硬币，并计算 f 的第 n 个上界和下界（给定迄今为止的正面数），称为 L 和 U。
    2. 计算 f 的 (n−1)th 上界和下界（给定迄今为止的正面数），称为 L′ 和 U′。（这些界限必须与未来硬币翻转的结果无关，并且区间 [L′, U′] 必须等于或完全包含区间 [L, U]。）

    这些算法部分适用于任何两个收敛于 f 的函数序列（不仅仅是多项式），其中 L 或 L′ 和 U 或 U′ 是它们的下界和上界近似值。通用工厂函数部分展示了如何为多项式实现此算法。但是，当近似函数（收敛于 f 的函数）是系数为整数的有理函数时，这些步骤如何工作？系数为有理数的有理函数？任意近似函数？

  5. 推入式自动机 是一种状态机，它持有一个符号堆栈。Mossel 和 Peres (2005)(3) 研究了这些机器在给定概率为 λ 的硬币的无限“磁带”翻转时可以模拟哪些函数 (f(λ))。他们证明了推入式自动机只能模拟代数函数，但可能并非所有代数函数。（参见“某些代数函数”。）问题是：推入式自动机可以模拟的代数函数的精确类别是什么？我写了一篇 文章附录 来展示我的进展，但关于这个问题还有其他结果吗？

正确性和性能图表

展示这些算法的正确性和性能的图表可在 单独页面 上找到。

致谢

我感谢 Luis Mendo（他回应了我的一项开放性问题）以及 C. Karney。

注释

  • (1) Flajolet, P., Pelletier, M., Soria, M., "On Buffon machines and numbers", arXiv:0906.5560 [math.PR], 2010.
  • (2) Keane, M. S., and O'Brien, G. L., "A Bernoulli factory", ACM Transactions on Modeling and Computer Simulation 4(2), 1994.
  • (3) 有一个与伯努利工厂问题类似的叫做量子伯努利工厂的问题，目标是模拟未知概率的函数，但这次使用涉及量子力学运算的算法（与不使用此类运算的经典算法不同）。然而，量子力学编程的普及程度远不及经典编程，并且在可预见的未来很可能仍将如此。因此，量子伯努利工厂超出了本文档的范围，但应指出的是，使用量子力学运算可以构造比经典算法更多的工厂函数。例如，在 [0, 1] 上定义的工厂函数必须满足 Keane 和 O'Brien 证明的要求，但它可以在定义域的有限个点处触摸 0 和/或 1（Dale, H., Jennings, D. and Rudolph, T., 2015, "Provable quantum advantage in randomness processing", Nature communications 6(1), pp.1-4）。
  • (4) Huber, M., "Nearly optimal Bernoulli factories for linear functions", arXiv:1308.1562v2 [math.PR], 2014.
  • (5) Yannis Manolopoulos. 2002. "Binomial coefficient computation: recursion or iteration?", SIGCSE Bull. 34, 4 (December 2002), 65–67. DOI: https://doi.org/10.1145/820127.820168.
  • (6) Goyal, V. and Sigman, K., 2012. On simulating a class of Bernstein polynomials. ACM Transactions on Modeling and Computer Simulation (TOMACS), 22(2), pp.1-5.
  • (7) Weikang Qian, Marc D. Riedel, Ivo Rosenberg, "Uniform approximation and Bernstein polynomials with coefficients in the unit interval", European Journal of Combinatorics 32(3), 2011, https://doi.org/10.1016/j.ejc.2010.11.004 http://www.sciencedirect.com/science/article/pii/S0195669810001666
  • (8) Wästlund, J., "Functions arising by coin flipping", 1999.
  • (9) Qian, W. and Riedel, M.D., 2008, June. The synthesis of robust polynomial arithmetic with stochastic logic. In 2008 45th ACM/IEEE Design Automation Conference (pp. 648-653). IEEE.
  • (10) Thomas, A.C., Blanchet, J., "A Practical Implementation of the Bernoulli Factory", arXiv:1106.2508v3 [stat.AP], 2012.
  • (11) S. Ray, P.S.V. Nataraj, "A Matrix Method for Efficient Computation of Bernstein Coefficients", Reliable Computing 17(1), 2012.
  • (12) 并且这表明，如果允许 c 等于 1，则该多项式无法模拟，因为所需的次数将是无穷大；事实上，在这种情况下，该多项式将在 0.5 点处触及 1，从而排除了任何算法的模拟（参见前面“关于伯努利工厂”）。
  • (13) Niazadeh, R., Leme, R.P., Schneider, J., "Combinatorial Bernoulli Factories: Matchings, Flows, and Polytopes", in Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pp. 833-846, June 2021; also at https://arxiv.org/abs/2011.03865.pdf.
  • (14) Mossel, Elchanan, and Yuval Peres. New coins from old: computing with unknown bias. Combinatorica, 25(6), pp.707-724, 2005.
  • (15) Nacu, Şerban, and Yuval Peres. "Fast simulation of new coins from old", The Annals of Applied Probability 15, no. 1A (2005): 93-115.
  • (16) Morina, G., Łatuszyński, K., et al., "From the Bernoulli Factory to a Dice Enterprise via Perfect Sampling of Markov Chains", arXiv:1912.09229 [math.PR], 2019/2020.
  • (17) Propp, J.G., Wilson, D.B., "Exact sampling with coupled Markov chains and applications to statistical mechanics", 1996.
  • (18) Flajolet 论文中定理 3.2 中给出的概率，即“∑ k = 0, 1, 2, ... (W(k) * (λ/2)k)”，似乎与该论文的图 4 不符。
  • (19) Flajolet, Ph., "Analytic models and ambiguity of context-free languages", Theoretical Computer Science 49, pp. 283-309, 1987
  • (20) 这里，“choose(g, g/t)”表示在 g 个字母中，必须有 g/t 个是 A，而“(β−1)g−g/t”是在其余字母为 A 的情况下，具有 g−g/t 个非 A 字母的词的数量。
  • (21) 在这个公式中，它与示例 2 类似，除以 βg*α−g 将 W(g) 从区间 [0, βg*α]（(g*α)-字母词）带到区间 [0, βg]（g-字母词），如主算法所需。
  • (22) Mendo, Luis. "An asymptotically optimal Bernoulli factory for certain functions that can be expressed as power series." Stochastic Processes and their Applications 129, no. 11 (2019): 4366-4384.
  • (23) Łatuszyński, K., Kosmidis, I., Papaspiliopoulos, O., Roberts, G.O., "Simulating events of unknown probabilities via reverse time martingales", arXiv:0907.4018v2 [stat.CO], 2009/2011.
  • (24) Flegal, J.M., Herbei, R., "Exact sampling from intractible probability distributions via a Bernoulli factory", Electronic Journal of Statistics 6, 10-37, 2012.
  • (25) Holtz, O., Nazarov, F., Peres, Y., "New Coins from Old, Smoothly", Constructive Approximation 33 (2011).
  • (26) Brassard, G., Devroye, L., Gravel, C., "Remote Sampling with Applications to General Entanglement Simulation", Entropy 2019(21)(92), https://doi.org/10.3390/e21010092 .
  • (27) Devroye, L., Non-Uniform Random Variate Generation, 1986.
  • (28) Bill Gosper, "Continued Fraction Arithmetic", 1978.
  • (29) Borwein, J. et al. “Continued Logarithms and Associated Continued Fractions.” Experimental Mathematics 26 (2017): 412 - 429.
  • (30) Mendo, L., "Simulating a coin with irrational bias using rational arithmetic", arXiv:2010.14901 [math.PR], 2020.
  • (31) 误差项，根据泰勒定理得出，其分子为 2，因为 2 大于 f 的斜率、斜率的斜率等函数在点 1 处可达到的最大值（在 cosh(1) 中）。
  • (32) Kozen, D., "Optimal Coin Flipping", 2014.
  • (33) 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.
  • (34) 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).
  • (35) Huber, M., "Optimal linear Bernoulli factories for small mean problems", arXiv:1507.00843v2 [math.PR], 2016.
  • (36) Schmon, S.M., Doucet, A. and Deligiannidis, G., 2019, April. Bernoulli race particle filters. In The 22nd International Conference on Artificial Intelligence and Statistics (pp. 2350-2358).
  • (37) Agrawal, S., Vats, D., Łatuszyński, K. and Roberts, G.O., 2021. "Optimal Scaling of MCMC Beyond Metropolis", arXiv:2104.02020.
  • (38) exp(−λ) 还有另一个算法，但当 λ 接近 1 时，它收敛缓慢。一个是通用鞅算法，因为当 λ 在 [0, 1] 内时，此函数是形式为 1 - x + x^2 - x^3 + ... 的交错级数，其系数为 1, 1, 1, 1, ...。另一个是 Flajolet et al. 2010 中所谓的“偶校验”构造：“(1) 翻转输入硬币。如果它返回 0，则返回 1。(2) 翻转输入硬币。如果它返回 0，则返回 0。否则，转到步骤 1。”
  • (39) Gonçalves, F. B., Łatuszyński, K. G., Roberts, G. O. (2017). Exact Monte Carlo likelihood-based inference for jump-diffusion processes.
  • (40) 这个函数还有另外两种算法，但当 λ 非常接近 1 时，它们都收敛得非常慢。一种是通用鞅算法，因为当 λ 在 [0, 1] 内时，该函数是形式为 1 - x + x^2 - x^3 + ... 的交错级数，其系数为 1, 1, 1, 1, ...。另一种是 Flajolet et al. 2010 中所谓的“偶校验”构造：“(1) 翻转输入硬币。如果它返回 0，则返回 1。(2) 翻转输入硬币。如果它返回 0，则返回 0。否则，转到步骤 1。”
  • (41) Vats, D., Gonçalves, F. B., Łatuszyński, K. G., Roberts, G. O., "Efficient Bernoulli factory MCMC for intractable posteriors", Biometrika, 2021 (也收录于 arXiv:2004.07471 [stat.CO])。
  • (42) Huber, M., "Designing perfect simulation algorithms using local correctness", arXiv:1907.06748v1 [cs.DS], 2019.
  • (43) Lee, A., Doucet, A. and Łatuszyński, K., 2014. "Perfect simulation using atomic regeneration with application to Sequential Monte Carlo", arXiv:1407.5770v1 [stat.CO].
  • (44) 这个函数还有一个算法，它使用通用鞅算法，但当 λ 接近 1 时，平均需要更多的比特。这里，交错级数是 1 - x + x^2/2 - x^3/3 + ...，其系数为 1, 1, 1/2, 1/3, ...
  • (45) 我能找到的为数不多的（如果不是唯一的）实现之一是 Haskell 实现。
  • (46) Penaud, J.G., Roques, O., "Tirage à pile ou face de mots de Fibonacci", Discrete Mathematics 256, 2002.
  • (47) Canonne, C., Kamath, G., Steinke, T., "The Discrete Gaussian for Differential Privacy", arXiv:2004.00010 [cs.DS], 2020.
  • (48) Forsythe, G.E., "Von Neumann's Comparison Method for Random Sampling from the Normal and Other Distributions", Mathematics of Computation 26(120), October 1972.
  • (49) Citterio, M., Pavani, R., "A Fast Computation of the Best k-Digit Rational Approximation to a Real Number", Mediterranean Journal of Mathematics 13 (2016).
  • (50) Sondow, Jonathan. “New Vacca-Type Rational Series for Euler's Constant and Its 'Alternating' Analog ln 4/π.”, 2005.
  • (51) von Neumann, J., "Various techniques used in connection with random digits", 1951.
  • (52) Pae, S., "Random number generation using a biased source", dissertation, University of Illinois at Urbana-Champaign, 2005.
  • (53) Peres, Y., "Iterating von Neumann's procedure for extracting random bits", Annals of Statistics 1992,20,1, p. 590-597.
  • (54) 将 λ 估计为 λ′，然后找到 f(λ′)，不一定是 f(λ) 的无偏估计量，即使 λ′ 是无偏估计量。
  • (55) Glynn, P.W., "Exact simulation vs exact estimation", Proceedings of the 2016 Winter Simulation Conference, 2016.
  • (56) Flajolet, P., Sedgewick, R., Analytic Combinatorics, Cambridge University Press, 2009.
  • (57) Monahan, J.. "Extensions of von Neumann's method for generating random variables." Mathematics of Computation 33 (1979): 1065-1069.
  • (58) Tsai, Yi-Feng, Farouki, R.T., "Algorithm 812: BPOLY: An Object-Oriented Library of Numerical Algorithms for Polynomials in Bernstein Form", ACM Trans. Math. Softw. 27(2), 2001.

附录

 

随机算法与非随机算法

非随机化算法 是一种模拟算法，它除了输入硬币作为随机性来源外，不使用任何其他东西（与使用其他随机性来源的随机化算法相反）（Mendo 2019）(22)。随机化算法不是生成外部随机性，而是可以实现一个 随机性提取 程序，使用输入硬币本身来生成该随机性。这样，算法就变成了一个非随机化算法。例如，如果一个算法实现了**双硬币特殊情况**，通过在步骤 1 中生成一个随机比特，那么它可以替换生成该比特的操作，改为翻转输入硬币两次，直到第一次翻转得到 0 然后是 1，或第一次翻转得到 1 然后是 0，然后将结果分别作为 0 或 1（von Neumann 1951）(51)。

事实上，将一个头出现概率为 (λ) 的硬币变成一个头出现概率为 (τ = f(λ)) 的硬币所需的平均硬币翻转次数存在一个下限。它被称为熵界限（参见，例如，(Pae 2005)(52)，(Peres 1992)(53)），计算如下——

  • ((τ − 1) * ln(1 − τ) − τ * ln(τ)) / ((λ − 1) * ln(1 − λ) − λ * ln(λ))。

例如，如果 f(λ) 是一个常数，非随机化算法通常需要不断增加的硬币翻转次数来模拟该常数，如果输入硬币强烈倾向于正面或反面（正面出现概率为 λ）。请注意，此公式仅在除了硬币翻转外不允许任何其他随机性时才有效。

对于某些 λ 值，Kozen (2014)(32) 证明了一个更严格的此类下界，但这个下界通常是非平凡的，并假定 λ 是已知的。然而，如果 λ 为 1/2（输入硬币无偏），这个下界很简单：平均而言，至少需要 2 次输入硬币翻转才能模拟一个已知的常数 τ，除非 τ 是 1/(2n)（对于任何整数 n）的倍数。

一个函数 f(λ) 被称为强可模拟（Keane and O'Brien 1994）(22)，如果它允许一个非随机化伯努利工厂算法。强可模拟的一个显著例子是，如果输入硬币的正面概率 (λ) 既不是 0 也不是 1。

模拟概率与估计概率

如果一个算法——

  • 翻转一个具有未知正面概率 (λ) 的硬币，并且
  • 以一个取决于 λ 的概率 (f(λ)) 产生正面，

该算法是 f(λ) 的一个无偏估计量，它以概率 1 生成 [0, 1] 内的估计值（Łatuszyński et al. 2009/2011）(23)。因此，理论上可以通过以下方式模拟概率 f(λ)——

  1. 以某种方式找到 f(λ) 的无偏估计量；(54)
  2. 生成一个 [0,1] 内的均匀随机变量，称为 u；并且
  3. 如果 u 小于 v，则返回 1，否则返回 0。

然而，在实践中，这种方法容易出现许多错误，其中包括由于步骤 1 和 2 中使用固定精度（例如舍入和抵消）而导致的错误。出于这个原因，并且因为“精确采样”是本页的重点，本页不涵盖直接估计 λ 或 f(λ) 的算法。另请参见 (Mossel and Peres 2005, section 4)(14)。

只有工厂函数（如“关于伯努利工厂”中所定义）才能具有无偏估计算法，该算法以概率 1 输出 [0, 1] 内的估计值（Łatuszyński et al. 2009/2011）(23)。

Glynn (2016)(55) 区分了——

  • 精确模拟，或模拟与 g(X) 具有相同分布（相同的“形状”、“位置”和“尺度”）的方法，并以概率 1 结束；
  • 精确估计，或模拟与 g(X) 具有相同期望值（无偏估计量，不仅仅是一致或渐近无偏估计量）的方法，并以概率 1 结束；

其中 g(X) 是一个遵循期望分布的随机值，基于随机变量 X。同样，本页的重点是“精确采样”（精确模拟），而不是“精确估计”，但是正面概率为 λ 的输入硬币可以是任何 λ 的“精确估计量”（如上定义），它输出 0 或 1。

注意：偏差和方差是随机估计算法中错误的两个来源。无偏估计量没有偏差，但并非没有误差。在当前情况下，伯努利工厂 f(λ) 的方差等于 f(λ) * (1−f(λ))，并且可能高达 1/4。有一些方法可以减小这种方差，但它们超出了本文档的范围。估计算法的均方误差等于方差加上偏差的平方。

连对数模拟算法的正确性证明

定理。如果“连分数”中的算法以概率 1 终止，则它以恰好等于连分数 c 所代表的数字的概率返回 1，否则返回 0。

证明。此正确性证明利用了 Huber 的“完美模拟基本定理”（Huber 2019）(42)。使用 Huber 的定理需要证明两件事：

  • 根据假设，算法以概率 1 终止。
  • 其次，当循环中的递归调用被一个模拟正确“连续子对数”的“黑盒”替换时，我们证明该算法是局部正确的。如果步骤 1 达到最后一个系数，则算法显然以正确的概率通过。否则，我们将模拟概率 (1 / 2c[i]) / (1 + x)，其中 x 是“连续子对数”，并且根据构造最多为 1。步骤 2 定义了一个循环，该循环将概率空间分为三个部分：第一部分占一半，第二部分（在第二个子步骤中）占另一半的一部分（此处等于 x/2），最后一部分是“拒绝部分”，它会重新运行循环。由于此循环不更改会影响后续迭代的任何变量，因此每次迭代都像一个接受/拒绝算法，该算法已被证明是 Huber 的完美模拟器。算法将在第一个子步骤中以概率 p = (1 / 2c[i]) / 2 通过，并在循环的第一个子步骤中以概率 f1 = (1 − 1 / 2c[i]) / 2 或在第二个子步骤中以概率 f2 = x/2 失败（所有这些概率相对于整个迭代）。最后，将通过数除以通过数和失败数的总和（p / (p + f1 + f2)）得到 (1 / 2c[i]) / (1 + x)，这就是我们想要的概率。

由于 Huber 定理的两个条件都得到满足，因此证明完成。□

连分式模拟算法 3 的正确性证明

定理。假设一个广义连分数的偏分子为 b[i]，全部大于 0，其偏分母为 a[i]，全部大于等于 1，并且进一步假设每个 b[i]/a[i] 小于等于 1。那么，算法“连分数”中的算法 3 以恰好等于该连分数所代表的数字的概率返回 1，否则返回 0。

证明。我们在正确性证明中再次使用了 Huber 的“完美模拟基本定理”。

  • 该算法以概率 1 终止，因为随着每次递归，该方法永远不会比不递归地进行更多的递归；请注意，a[i] 永远不能大于 1，因此在每次迭代中完成运行的概率 a[i]/(1+a[i]) 始终是 1/2 或更大。
  • 如果循环中的递归调用被一个模拟正确“子分数”的“黑盒”替换，则算法在局部是正确的。如果步骤 1 达到连分数的最后一个元素，则算法显然以正确的概率通过。否则，我们将模拟概率 b[i] / (a[i] + x)，其中 x 是“连续子分数”，根据假设最多为 1。步骤 2 定义了一个循环，该循环将概率空间分为三个部分：第一部分占 h = a[i]/(a[i] + 1) 的一部分，第二部分（在第二个子步骤中）占剩余部分的一部分（此处等于 x * (1 − h)），最后一部分是“拒绝部分”。算法将在第一个子步骤中以概率 p = (b[i] / a[pos]) * h 通过，并在循环的第一个子步骤中以概率 f1 = (1 − b[i] / a[pos]) * h 或在第二个子步骤中以概率 f2 = x * (1 − h) 失败（所有这些概率相对于整个迭代）。最后，将通过数除以通过数和失败数的总和得到 b[i] / (a[i] + x)，这就是我们想要的概率，因此 Huber 的两个条件都得到满足，证明完成。□

冯·诺依曼方案

(Flajolet et al., 2010)(1) 描述了他们称之为冯·诺依曼模式（第 2 节）。尽管冯·诺依曼模式用于此处给出的几个伯努利工厂，但它本身并不是一个伯努利工厂，因为它可能产生大于 0 和 1 的随机变量，因此本节出现在附录中。给定一个排列类和一个输入硬币，冯·诺依曼模式会生成一个随机整数 n，大于等于 0，其概率等于——

  • (λn * V(n) / n!) / EGF(λ)，

其中—

  • EGF(λ) = ∑k = 0, 1, ... (λk * V(k) / k!)（指数生成函数或 EGF，它完全确定了一个排列类），并且
  • V(n) 是一个在 [0, n!] 范围内的数字，并且是大小为 n 的排列中满足所讨论的排列类要求的数量。

实际上，通过翻转硬币直到它返回 0 并计算 1 的数量（该论文称 G 为几何(λ) 随机变量，但在本文档中避免使用此术语，因为它有多个冲突的含义），然后以 V(G)/(G!) 的概率接受，否则拒绝，来生成一个随机变量 G。r 次变量被拒绝的概率是 p*(1 − p)r，其中 p = (1 − λ) * EGF(λ)。

排列类的例子包括——

  • 单循环排列（EGF(λ) = Cyc(λ) = ln(1/(1 − λ)); V(n) = (n − 1)!）
  • 排序排列，或数字按降序排列的排列（EGF(λ) = Set(λ) = exp(λ); V(n) = 1），
  • 所有排列（EGF(λ) = Seq(λ) = 1/(1 − λ); V(n) = n!），
  • 偶数大小的交替排列（EGF(λ) = 1/cos(λ); 从 n = 0 开始的 V(n) 是在线整数序列百科全书中的A000364），以及
  • 奇数大小的交替排列（EGF(λ) = tan(λ); 从 n = 0 开始的 V(n) 是A000182），

使用“Analytic Combinatorics”（Flajolet and Sedgewick 2009）(56) 中的符号。

以下算法生成一个遵循冯·诺依曼模式的随机变量。

  1. 将 r 设置为 0。（这是算法拒绝随机变量的次数。）
  2. 翻转输入硬币，直到翻转返回 0。然后将返回 1 的次数设为G。
  3. 以 V(G)/G! 的概率返回 G（如果需要，也可以返回 r）。（实际上，概率检查是通过生成 G 个 uniform(0, 1) 随机变量并确定这些数字是否满足给定的排列类来完成的。这尤其因为 G!（G 的阶乘）很容易变得非常大。）
  4. 将 r 加 1 并转到步骤 2。

根据 EGF 和 Bernoulli 工厂算法处理的 G 和/或 r 的值，冯·诺依曼模式可以产生各种伯努利工厂概率函数。以下 Python 函数使用 SymPy 计算机代数库来查找给定排列类的 EGF 的概率和其他有用信息，以应用冯·诺依曼模式。

def coeffext(f, x, power):
    # Extract a coefficient from a generating function
    # NOTE: Can also be done with just the following line:
    # return diff(f,(x,power)).subs(x,0)/factorial(power)
    px = 2
    for i in range(10):
      try:
        poly=Poly(series(f, x=x, n=power+px).removeO())
        return poly.as_expr().coeff(x, power)
      except:
        px+=2
    # Failed, assume 0
    return 0

def number_n_prob(f, x, n):
    # Probability that the number n is generated
    # for the von Neumann schema with the given
    # exponential generating function (e.g.f.)
    # Example: number_n_prob(exp(x),x,1) --> x**exp(-x)
    return (x**n*coeffext(f, x, n))/f

def r_rejects_prob(f, x, r):
    # Probability that the von Neumann schema
    # with the given e.g.f. will reject r random variates
    # before accepting the next one
    p=(1-x)*f
    return p*(1-p)**r

def valid_perm(f, x, n):
    # Number of permutations of size n that meet
    # the requirements of the permutation class
    # determined by the given e.g.f. for the
    # von Neumann schema
    return coeffext(f, x, n)*factorial(n)

注意：冯·诺依曼模式可以模拟任何幂级数分布（例如泊松分布、负二项分布、几何分布和对数级数），只要有一个合适的指数生成函数。然而，由于步骤 2，当 λ 接近 1 时，模式所需的输入硬币翻转次数会无界增长。

示例：使用排序排列类，我们可以通过冯·诺依曼模式生成一个 Poisson(λ) 随机变量，其中 λ 是输入硬币的正面概率。这将导致 exp(−λ) 的算法——如果 Poisson(λ) 随机变量为 0，则返回 1，否则返回 0——但由于上述原因，当 λ 接近 1 时，此算法收敛缓慢。此外，如果 c > 0 是一个实数，将 Poisson(floor(c)) 随机变量添加到 Poisson(c−floor(c)) 随机变量会生成 Poisson(c) 随机变量。

冯·诺依曼模式的一个变体发生在 G 的生成方式不同于步骤 2 中给出的情况，但仍然是通过翻转输入硬币生成的。在这种情况下，算法将返回 n 的概率为——

  • (κ(n; λ)*V(n)/(n!)) / p，

其中 p = ( ∑k=0,1,... (κ(k; λ)*V(k)/(k!)) )，并且 κ(n; λ) 是 G 为 n 的概率，参数为 λ 或输入硬币的正面概率。同样，修改后的算法拒绝 r 个随机变量的概率是 p*(1 − p)r。

示例：如果 G 是 Poisson(z2/4) 随机变量且使用排序排列类，则该算法将以 1/I0(z) 的概率返回 0，其中 I0(.) 是第一类修正贝塞尔函数。

某些排列产生的概率

某些有趣的概率函数可以从排列中产生，例如排序排列或最高数字出现在第一个的排列。受前面附录中冯·诺依曼模式的启发，我们可以描述以下算法：

给定一个排列类（例如降序排列）和两个连续概率分布 D 和 E。考虑以下算法：生成一个独立随机变量序列（第一个分布为 D，其余为 E），直到序列不再遵循排列类，然后返回 n，即生成的数字个数减 1。

那么算法的行为如下表所示。

排列类分布 D 和 E算法返回 n 的概率n 为...的概率
降序排列两者均为 uniform(0,1)n / ((n + 1)!).奇数为 1−exp(−1)。
偶数为 exp(−1)。
降序排列任意(∫(−∞,∞) DPDF(z) * (ECDF(z)n−1/((n−1)!) − ECDF(z)n/(n!)) dz)，对于每个 n > 0（另见 (Devroye 1986, Chapter IV)(27) 定理 2.1 的证明。DPDF 和 ECDF 稍后定义。奇数为下面公式 1 的分母。
交替排列两者均为 uniform(0,1)(an * (n + 1) − an + 1) / (n + 1)!，其中 ai 是序列 A000111 在在线整数序列百科全书中位置 i（从 0 开始）处的整数。奇数为 1−cos(1)/(sin(1)+1)；偶数为 cos(1)/(sin(1)+1)。
任意两者均为 uniform(0,1)(∫[0, 1] 1 * (zn−1*V(n)/((n−1)!) − zn*V(n+1)/(n!)) dz)，对于每个 n > 0。V(n) 是大小为 n 的排列中满足排列类要求的数量。对于此算法，V(n) 必须在 (0, n!] 区间内；例如，如果有 0 个奇数大小的排列，此算法将无法工作。奇数为 1 − 1 / EGF(1)；偶数为 1/EGF(1)。
小于 k 为 (V(0) − V(k)/(k!)) / V(0)。
排列类分布 D 和 E给定 n 为...时，序列中的第一个数字小于 x 的概率
降序排列任意奇数为 ψ(x) = (∫(−∞, x) exp(−ECDF(z)) * DPDF(z) dz) / (∫(−∞, ∞) exp(−ECDF(z)) * DPDF(z) dz) （公式 1；参见 (Devroye 1986, Chapter IV)(27) 定理 2.1(iii)；另参见 Forsythe 1972(48)）。此处，DPDF 是 D 的概率密度函数 (PDF)，ECDF 是 E 的累积分布函数 (CDF)。
如果 x 是 uniform(0, 1)，则此概率变为 ∫[0, 1] ψ(z) dz。
降序排列任意偶数为 (∫(−∞, x) (1 − exp(−ECDF(z))) * DPDF(z) dz) / (∫(−∞, ∞) (1 − exp(−ECDF(z))) * DPDF(z) dz) （公式 2；另参见 Monahan 1979(57)）。DPDF 和 ECDF 如上所述。
降序排列两者均为 uniform(0,1)奇数为 ((1−exp(−x))−exp(1))/(1−exp(1))。因此，序列中的第一个数字分布为 exponential(1) 并“截断”到区间 [0, 1]（von Neumann 1951）(51)。
降序排列D 是 uniform(0,1)；E 是两个 uniform(0,1) 的最大值奇数为 erf(x)/erf(1) （使用公式 1，其中 z 在 [0, 1] 内时 DPDF(z) = 1 且 ECDF(z) = z2；另参见 erf(x)/erf(1)）。

注释

  1. 公式 1 和 2 中所有可能的函数都是非递减函数。这两个公式都表达了累积分布函数 FD(x | n 为奇数) 或 FD(x | n 为偶数)，分别。
  2. EGF(z) 是算法中涉及的排列类型的指数生成函数（EGF）。例如，交替排列类（数字在低和高之间交替的排列，即 X1 > X2 < X3 > ...）使用 EGF tan(λ)+1/cos(λ)。在冯·诺依曼模式部分给出了 EGF 的其他示例。

开放性问题：如何为其他排列类和不同的分布组合 D 和 E 填充上表？

1 / π 算法推导草图

Flajolet 论文描述了一个模拟 1 / π 的算法，但未提供推导。以下是该算法如何工作的草图。

该算法是凸组合技术的一个应用。即，1 / π 可以看作是两个分量的凸组合

  • g(n): 26 * n * (6 * n + 1) / 28 * n + 2 = 2−2 * n * (6 * n + 1) / 4 = (6 * n + 1) / (22 * n + 2)，这是以下独立随机变量之和等于 n 的概率

    • 两个随机变量，每个变量都表示第一次成功之前的失败次数，其中成功的几率为 1−1/4（该论文称这两个数字为几何(1/4) 随机变量，但在本文档中避免使用此术语，因为它有多个冲突的含义）。
    • 一个 Bernoulli(5/9) 随机变量。

    这对应于凸组合算法的步骤 1 和 1 / π 算法的步骤 2 到 4。（这也表明 Flajolet 论文中给出的 1 / π 的恒等式存在错误：“8 n + 4”应为“8 n + 2”。）

  • hn(): (choose(n * 2, n) / 2n * 2)3，这是编号为 n 的“硬币”的正面概率。这对应于凸组合算法的步骤 2 和 1 / π 算法的步骤 5。

注释

  1. 9 * (n + 1) / (22 * n + 4) 是两个独立随机变量之和等于 n 的概率（负二项分布），其中每个数字表示第一次成功之前的失败次数，而成功的几率为 p。
  2. pm * (1 − p)n * choose(n + m − 1, m − 1) 是 m 个独立随机变量之和等于 n 的概率（负二项分布），其中每个 m 个数字表示第一次成功之前的失败次数，而成功的几率为 p。
  3. f(z) * (1 − p) + f(z − 1) * p 是两个独立随机变量（一个 Bernoulli(p) 数和一个概率质量函数为 f(.) 的整数 z）之和等于 z 的概率。

计算 exp(1) 的界

以下实现了 Citterio 和 Pavani 算法 (2016)(49) 中计算 exp(1) 的有理数形式的下界和上界所需的部分。

定义以下运算

  • 设置：将 p 设置为列表 [0, 1]，将 q 设置为列表 [1, 0]，将 a 设置为列表 [0, 0, 2]（两个零，后跟 exp(1) 的整数部分），将 v 设置为 0，将 av 设置为 0。
  • 确保 n：当 v 小于等于 n 时
    1. (确保从 0 开始的部分分母 v 可用。) 如果 v + 2 大于等于 a 的大小，则按顺序将 1、av 和 1 追加到列表 a，然后将 av 加 2。
    2. (计算从 0 开始的收敛项 v。) 将 a[n+2] * p[n+1]+p[n] 追加到列表 p，并将 a[n+2] * q[n+1]+q[n] 追加到列表 q。（列表中的位置从 0 开始。例如，p[0] 表示 p 中的第一项；p[1] 表示第二项；依此类推。）
    3. 将 v 加 1。
  • 获取收敛项 n 的分子：确保 n，然后返回 p[n+2]。
  • 获取收敛项 n：确保 n，然后返回 p[n+2]/q[n+2]。
  • 给定 d 的半收敛项 n
    1. 确保 n，然后将 m 设置为 floor(((10d)−1−p[n+1])/p[n+2])。
    2. 返回 (p[n+2] * m +p[n+1]) / (q[n+2] * m +q[n+1])。

然后计算 exp(1) 的下界和上界的算法如下，给定 d：

  1. 将 i 设置为 0，然后运行 **设置**。
  2. **获取收敛项 i 的分子**，称之为 c。如果 c 小于 10d，则将 i 加 1 并重复此步骤。否则，转到下一步。
  3. **获取收敛项 i − 1** 和 **获取给定 d 的半收敛项 i − 1**，分别称为 conv 和 semi。
  4. 如果 (i − 1) 是奇数，则返回 semi 作为下界，返回 conv 作为上界。否则，返回 conv 作为下界，返回 semi 作为上界。

准备有理函数

本节描述了如何将单变量有理函数（多项式之比）转换为应用**“Dice Enterprise”特殊情况**（在“某些有理函数”中描述）所需的系列多项式。简而言之，执行此操作的步骤可描述为分离、齐次化和增强。

分离。 如果有理函数的分子（D）和分母（E）被写成——

  • 形如 z*λi*(1−λ)j 的项之和，其中 z 是实数，i≥0 和 j≥0 是整数（在本节中称为形式 1），

则该函数可以分离成两个相加等于分母的多项式。（此处，i+j 是该项的次数，多项式的次数是其项中的最高次数。）为了进行此分离，从分母中减去分子得到一个新多项式（G），使得 G = E − D（或 D + G = E）。（然后 D 和 G 将是我们使用的两个多项式。）类似地，如果我们有多个具有共同分母的有理函数，即 (D1/E), ..., (DN/E)，其中 D1, ..., DN 和 E 以形式 1 表示，则可以通过从分母中减去分子来将它们分离成 N + 1 个多项式，使得 G = E − D1 − ... − DN。（然后 D1, ..., DN 和 G 是我们将使用的多项式。）然而，为了在算法中使用这些多项式，需要对它们进行齐次化，然后增强，如下面所述。

示例：我们有有理函数 (4*λ1*(1−λ)2) / (7 − 5*λ1*(1−λ)2)。从分母中减去分子得到：7 − 1*λ1*(1−λ)2。

齐次化。 下一步是齐次化这些多项式，使它们具有相同的次数和特定的形式。对于此步骤，选择 n 为一个不小于多项式中最高次数的整数。

假设一个多项式——

  • 对于 [0, 1] 区间内的每个 λ，其值大于等于 0，
  • 次数小于等于 n，并且
  • 如上所示，以形式 1 表示。

那么该多项式可以通过以下方式转换为一个*齐次多项式*（次数为 n，即其所有项的次数都为 n）。

  • 对于 [0, n] 范围内的每个整数 m，新齐次多项式中 m 的系数计算如下：
    1. 将 r 设置为 0。
    2. 对于（旧多项式中）形如 z*λi*(1−λ)j 的每个项：
      • 如果 i ≤ m，且 (n−m) ≥ j，且 i + j ≤ n，则将 z*choose(n−(i+j), (n−m)−j) 加到 r 中。
    3. 现在，r 就是新的系数（对应于 r* λm*(1−λ)n−m 的项）。

如果多项式以所谓的“幂形式”表示为 c[0] + c[1]*λ + c[2]*λ2 + ... + c[n]*λn，则方法如下：

  • 对于 [0, n] 范围内的每个整数 m，新齐次多项式中 m 的系数计算如下：
    1. 将 r 设置为 0。
    2. 对于 [0, m] 范围内的每个整数 i，如果存在系数 c[i]，则将 c[i]*choose(n−i, n−m) 加到 r 中。
    3. 现在，r 就是新的系数（对应于 r* λm*(1−λ)n−m 的项）。

示例：我们有以下多项式：3*λ2 + 10*λ1*(1−λ)2。这是一个 3 次多项式，我们希望将其转换为 5 次齐次多项式。结果将是各项的总和：

  • 0 * λ0*(1−λ)5；
  • 10*choose(2, 2) * λ1*(1−λ)4 = 10* λ1*(1−λ)4；
  • (3*choose(3, 3) + 10*choose(2, 1)) * λ2*(1−λ)3 = 23* λ2*(1−λ)3；
  • (3*choose(3, 2) + 10*choose(2, 0)) * λ3*(1−λ)2 = 19* λ3*(1−λ)2；
  • 3*choose(3, 1) * λ4*(1−λ)1 = 9* λ4*(1−λ)1；和
  • 3*choose(3, 0) * λ5*(1−λ)0 = 3* λ5*(1−λ)0，

最终得到新齐次多项式的系数（0, 10, 23, 19, 9, 3）。

增强。如果我们有一个同次数的齐次单变量多项式数组，它们就准备好用于Dice Enterprise 特例，前提是：

  • 这些多项式具有相同的次数，即 n，
  • 它们的系数都大于或等于 0，并且
  • 对于从 0 到 n 的每个 j，第 j 个系数的总和大于 0，但系数之和列表可以以零开始和/或结束。

如果未满足这些条件，则每个多项式都可以根据需要进行增强以满足条件（Morina 等人，2019）(16)。对于这里相关的多项式，增强一个多项式相当于进行次数提升，这与 Bernstein 形式的多项式类似（另见 Tsai 和 Farouki 2001(58)）。实现方法如下：

  • 设 n 为多项式的旧次数。对于 [0, n+1] 范围内的每个 k，新多项式中 k 的系数计算如下：
    • 设 c[j] 为旧多项式的第 j 个系数（从 0 开始）。对于 j 在区间 [max(0, k−1), min(n, k)] 内的每个 j，计算 c[j] * choose(1, k−j)，然后将它们相加。总和即为新系数。

根据 Morina 的论文，对每个多项式进行 n 次增强足以使整个数组满足上述条件（尽管少于 n 次通常也足够）。

注意：为获得最佳结果，输入多项式的系数应为有理数。如果不是，则需要采用特殊方法来确保精确结果，例如区间算术，用于计算下界和上界。

© . All rights reserved.