飞镖、骰子和硬币：从离散分布中抽样

本篇是原文 Darts, Dice, and Coins: Sampling from a Discrete Distribution 的全文中文翻译版本

在被我的导师拉着参与 cloudwego/kitex 的一个负载均衡算法的改进时 #1184 ，他指出可能可以使用 Alias Method（别名方法）来实现更为高效的负载均衡。于是我仔细地研究了本文的英文原文（链接如上）， 起初我想寻求一个中文译版，但貌似没有找到

于是鄙人借助 GPT 翻译与个人的理解进行部分词汇修正，手打大部分 LaTeX 公式与表格后完成本文，如有谬误，欢迎指出

PS：暗黑主题可能不适合本文阅读，可在右下角切换

原文最近更新日期：December 29, 2011

---

今年早些时候，我在 Stack Overflow 上提出了一个关于不公平的骰子的数据结构的问题。具体来说，我感兴趣的问题是：

“你有一个 $n$ 面的骰子，第 $i$ 面被掷出的概率为 $p_{i}$。模拟掷骰子的最有效数据结构是什么？”

这种数据结构可以用于多种目的。首先，你可以用它来模拟一个公平的六面骰子，通过给骰子的每个面分配 $\frac{1}{6}$ 的概率，或者模拟一个公平的硬币，通过模拟一个两面的骰子，每面有 $\frac{1}{2}$ 的概率出现。你还可以使用这个数据结构直接模拟两个公平六面骰子掷出的总和，通过使用一个 11 面的骰子（面数为 2、3、4、…、12），每个面都适当地加权，以表示如果使用两个公平骰子，这个总数出现的概率。然而，你也可以使用这个数据结构来模拟不公平的骰子。例如，如果你在玩花旗骰，而你知道这些骰子不是完全公平的，你可能会使用这个数据结构来模拟掷骰子的许多次，看看最佳策略会是什么。你也可以考虑以同样的方式模拟一个不完美的轮盘赌。

在游戏领域之外，您还可以在机器人模拟中使用这个数据结构，其中传感器有已知的故障率。例如，如果一个测距传感器有 95% 的几率返回正确的值，4% 的几率返回一个太小的值，和 1% 的几率返回一个太大的值，您可以使用这个数据结构来模拟传感器的读数，通过生成一个随机结果，并在那种情况下模拟传感器的读数。

我在 Stack Overflow 上收到的答案给我留下了深刻的印象，原因有两个。首先，解决方案指出了一种强大的技术，称为别名方法，它在对机器模型进行一些合理的假设后，能够在 $O(1)$ 时间内完成骰子的模拟掷出，仅需一个简单的预处理步骤。其次，也许更令人惊讶的是，这种算法已经被知道了几十年，但我之前从未遇到过它！考虑到有多少处理时间用于模拟，我本以为这种技术会更为人所知。几次快速的谷歌搜索找到了关于这种技术的大量信息，但我找不到一个单独的网站，能够汇编这种技术背后的直觉和解释。

这篇文章是我尝试快速概述模拟不公平的骰子的各种方法，从高度不切实际的简单技术到非常优化和高效的别名方法。我在这里的希望是捕捉到问题的不同直觉，以及每种方法如何突出模拟不公平的骰子的某些新方面。对于每种方法，我的目标是探索激发思想、核心算法、正确性证明以及运行时分析（就时间、内存和所需的随机性而言）。

---

引言

在我详细介绍不同技术的具体细节之前，让我们首先统一我们的符号和术语。

在这篇文章的引言中，我使用了“骰子”这个术语来描述一个通常情况，其中有一组有限的结果，每个结果都有一些相关的概率。正式来说，这被称为离散概率分布，模拟不公平的骰子的问题被称为从离散分布中抽样。为了描述我们的离散概率分布（不公平的骰子），我们将假设我们给定了一组 $n$ 个概率 $p_0, p_1, ..., p_{n-1}$ 与结果 $0, 1, ..., n-1$ 相关联。尽管结果可以是任何东西（正面/反面、骰子上的数字、颜色等），为了简单起见，我将假设结果是与给定索引相对应的某个正自然数。

在计算机上处理实数是一种复杂的计算领域。有许多快速算法的速度完全来自于在常数时间内计算任意实数的向下取整能力，而浮点表示中的数值不准确性可以完全破坏某些算法。因此，在我们进入关于概率算法的暗黑世界之前，我们应该界定清楚计算机能做什么和不能做什么。

在接下来的内容中，我将假设所有以下操作都可以在恒定时间内完成：

1. 加法、减法、乘法、除法和任意实数的比较：我们需要能够进行这些操作来处理概率。这看起来可能是一个非常强的假设，但计算机处理实数时的精度（即能够表示的实数的精确程度）受到其基本处理单元的限制（例如，32位机器上的64位双精度数），那么这其实是很合理的。

2. 生成范围在 [0, 1) 内的均匀实数：为了模拟随机性，我们需要某种随机源。我假设我们可以在恒定时间内生成任意精度的实数。这远远超出了计算机实际能力，但出于本讨论的目的，我认为这样做是可以的。如果我们允许使用 $[0, 1]$ 内的任意 IEEE-754 双精度数，我们确实会失去一些精度，但对于大多数应用来说可能已经足够准确了。

3. 计算实数的整数下限：我们使用 IEEE-754 双精度数完成这类操作是合理的，但在更广泛的计算环境中，这样的要求对计算机是不合理的。

值得考虑的是，假设我们能够高效地执行所有这些操作是否合理。在实践中，我们很少有概率遇到 IEEE-754 双精度浮点数固有因舍入误差而引起的重大精度问题，因此， IEEE 双精度数满足以上所有约束。当然，如果我们要求概率精确地*被指定为高精度有理数，那么这些约束可能是不合理的。

模拟公平骰子

在我们概括到模拟任意加载的骰子这一特殊情况之前，让我们从一个更简单的算法开始，这个算法将作为后续算法的基础：模拟一个公平的、$n$ 面的骰子。例如，我们可能对掷出一个公平的六面骰子来玩大富翁或风险游戏感兴趣，或者抛掷一个公平的硬币（一个两面的骰子）等。

在这个特殊情况下，有一个简单、优雅且高效的算法可以模拟结果。假设我们可以生成真正随机的、均匀分布的实数，范围在 $[0,1)$ 内。我们可以将这个范围可视化如下：

要模拟投掷 $n$ 面的骰子，一种方法是将范围 $[0,1)$ 分成 $n$ 个均匀大小的小区间，每个区间的长度为 $\frac{1}{n}$。如下所示：

如果我们生成一个在范围 $[0,1)$ 内随机选择的实数，它将会落入这些小区间中的一个。我们可以通过确定它落在哪个区间来得知骰子投掷的结果。例如，如果我们随机选择的值落在了这个位置：

我们会说骰子的结果是 $2$（假设骰子的编号从 $0$ 开始）。

从图形上看，很容易看出随机值落在哪个区间，但我们如何使用算法实现呢？这体现骰子的公平性。由于所有的区间都具有相同的大小，即 $\frac{1}{n}$，我们可以找出最大的 $i$ 使得 $\frac{i}{n}$ 不大于随机生成的值（称为 $x$）。如果我们试图找到最大的 $i$ 值，使得 $\dfrac{i}{n}\leq x$，这等价于找到最大的 $n$ 值，使得 $i\leq xn$。但是，根据定义，这意味着 $i=⌊xn⌋$，即不大于 $xn$ 的最大自然数。因此，这给我们提供了以下（非常简单的）模拟公平的 $n$ 面骰子的算法：

算法：模拟公平的骰子

1. 生成在范围 $[0,1)$ 内均匀随机的值 $x$。
2. 返回 $⌊xn⌋$。

根据我们上面的计算假设，这个算法的运行时间是 $O(1)$。

这一部分有两个要点。首先，我们可以将范围 $[0,1)$ 分割成若干块，以便范围内的均匀随机实数自然地映射回我们可用的众多离散选择之一。我们将在本文的其余部分广泛利用这一技巧。其次，确定特定随机值落入哪个范围可能会很复杂，但如果我们了解有关分区的一些信息（在这种情况下，它们都具有相同的大小），则可以在数学上轻松确定特定点位于哪个分区。

用公平骰子模拟不公平骰子

在给定模拟公平骰子的算法的基础上，我们是否可以调整该算法以模拟一个不公平的骰子？有趣的是，答案是肯定的，但这将占用更多的空间。

根据前面的经验，为了模拟一个不公平的骰子，我们只需要将范围 $[0,1)$ 分成若干块，然后确定我们落入了哪一块。然而，一般情况下，这可能比看起来更加困难。例如，假设我们有一个四面的骰子，其各面的概率分别为$\frac{1}{2}$、 $\frac{1}{3}$、 $\frac{1}{12}$ 和 $\frac{1}{12}$（我们可以确认这是一个合法的概率分布，因为 $\frac{1}{2} + \frac{1}{3} + \frac{1}{12} + \frac{1}{12} = \frac{6}{12} + \frac{4}{12} + \frac{1}{12} + \frac{1}{12} = \frac{12}{12}$）。如果我们将范围 [0,1) 分成四个具有这些给定大小的区间，我们得到以下结果：

不幸的是，我们在这里陷入了困境。即使我们知道范围 $[0,1)$ 内的一个随机数，也没有简单的数学技巧可以自动确定该数落入哪个分区。这并不是说这是极其困难的 - 正如您将看到的，有许多很棒的技巧可以使用 - 但没有一个一定具有像模拟公平骰子的算法那样的数学简单性。

然而，我们可以调整公平骰子的技巧以适应这种情况。让我们以这个特定情况为例。骰子各面出现的概率分别为 $\frac{1}{2}$、 $\frac{1}{3}$、 $\frac{1}{12}$ 和 $\frac{1}{12}$。如果我们将这些重新写成具有共同分母的形式，我们得到值 $\frac{6}{12}$、 $\frac{4}{14}$、 $\frac{1}{12}$ 和 $\frac{1}{12}$。因此，我们可以这样考虑这个问题：与其掷一个具有加权概率的四面骰子，为什么不掷一个具有重复标签的12面公平骰子？由于我们知道如何模拟公平骰子，这相当于将范围 $[0,1)$ 分成以下十二个区间：

然后将它们分配给不同的结果，如下所示：

现在，模拟一个骰子的投掷非常简单 - 我们只需投掷这个新的公平骰子，然后读取它的值。这个第一步可以使用上面的算法来完成，它将给我们返回一个在范围 $0,1,...,11$内的整数。要将该整数映射回原始的不公平骰子的一面，我们将存储一个大小为 12 的辅助数组，将这些数字中的每一个映射回原始的结果。从图形上看，我们可以如下所示：

为了将这个过程正式化为一个算法，我们将描述初始化步骤（创建映射表）和生成步骤（模拟随机骰子的投掷）。在这个算法和随后的算法中，这两个步骤都是重要的，因为设置时间可能会很长。

在初始化步骤中，我们首先找到我们所给定的骰子各面的概率的最小公倍数（在我们的示例中，这是12）。最小公倍数在这里很有用，因为它对应着我们可以用于所有分数的最小公共分母，因此也是新的公平骰子的面数。一旦我们有了这个最小公倍数（我们称之为 L），我们需要确定新骰子的每个面将分配给原始不公平骰子的哪个面。在我们的示例中，概率为 $\frac{1}{2}$ 的面在新骰子上有 6 个面，因为 $\frac{1}{2} \times 12=6$。同样地，概率为 $\frac{1}{3}$ 的面在新骰子上有 $4$ 个面，因为 $\frac{1}{3} \times 12=4$。更一般地，如果 $L$ 是概率的最小公倍数，$p_{i}$ 是骰子的第 $i$ 面出现的概率，那么我们将把 $L⋅p_{i}$ 个新骰子的面分配给原始不公平骰子的第 $i$ 面。

以下是上述算法的伪代码：

算法：使用公平骰子模拟不公平骰子

• 初始化：
  1. 找到概率 $p_{0}, p_{1}, ..., p_{n-1}$ 的分母的最小公倍数；将其命名为 $L$。
  2. 分配一个大小为 $L$ 的数组 $A$，用于映射公平骰子的结果到原始骰子的结果。
  3. 对于原始骰子的每一面 $i$，以任意顺序执行以下步骤：
     1. 将接下来 $L⋅p_{i}$ 个 $A$ 中元素设置为 $i$。
• 生成：
  1. 从一个 $L$ 面的公平骰子生成一个公平骰子的结果；将其命名为 $S$。
  2. 返回 $A[S]$。

这个算法可能很简单，但它的效率如何？实际的骰子投掷生成非常快速 - 使用之前的算法生成随机骰子投掷只需要 $O(1)$ 的工作量，再加上表格查找额外的 $O(1)$ 工作量。这总共需要的工作量是 $O(1)$。因此，这个算法在时间复杂度上是非常高效的。

然而，初始化步骤可能会非常耗时。为了使这个算法工作，我们需要分配一个与所有输入分数的分母的最小公倍数一样大的数组空间。对于我们的示例（$\frac{1}{2}$、 $\frac{1}{3}$、 $\frac{1}{12}$ 、$\frac{1}{12}$），这个最小公倍数是 $12$，但对于其他输入来说，它可能会变得非常差劲。举个例子，考虑分数 $\frac{999999}{1000000}$ 和 $\frac{1}{1000000}$。分母的最小公倍数是一百万，所以我们的表格将需要一百万个条目！

这意味着在某些情况下，初始化步骤可能会占用大量的空间。在考虑使用这个算法时，需要注意可能的空间复杂度问题。

不幸的是，情况可能比这更糟。在前面的例子中，我们至少可以说我们“预期”算法会使用大量内存，因为每个分数的分母都是一百万。然而，我们可能会遇到一组概率，其中最小公倍数远远大于任何单个分母。例如，考虑概率 $\frac{1}{15}$、$\frac{1}{10}$、$\frac{1}{6}$。在这里，分母的最小公倍数是$30$，大于任何一个分母本身。这个构造之所以有效，是因为$15=3×5$，$10=2×5$，$6=2×3$；换句话说，每个分母都是从三个素数中选择的两个素数的乘积。它们的最小公倍数因此是所有这些素数的乘积，因为每个分母都必须整除最小公倍数。如果我们将这个构造推广，并考虑任意一组 $k$ 个素数，那么如果我们为这些素数的每个两两乘积选择一个分数，那么最小公倍数可能远远大于任何单个分母。事实上，我们可以得出最小公倍数的一个最好的上界是 $O(\prod ^{n}_{i=1}d_{i})$，其中 $d_{i}$ 是第 $i$ 个概率的分母。这意味着，如果概率不是事先已知的，那么这个算法在任何实际场景中的使用都可能因为需要的内存太多而无法实现。

尽管如此，在许多情况下，这个算法的表现是良好的。如果所有的概率都相同，那么我们作为输入给定的所有概率都是 $\frac{1}{n}$，其中 $n$ 是某个整数。然后，分母的最小公倍数是 $n$，因此我们最终掷出的公平骰子将有 $n$ 个面，原始骰子的每个面将对应于公平骰子的一个面。初始化时间因此是 $O(n)$。从图形上看，这将如下所示：

下表总结了该算法的信息

算法	初始化时间（最好 - 最差）	生成时间	内存使用（最好 - 最差）
用公平的骰子模拟不公平的骰子	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$	$Θ(1)$	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$

这个算法的另一个重要细节是，它假设我们以分数的形式表示概率，这些分数具有适当的分母。如果概率以 IEEE-754 双精度浮点数的形式给出，那么这种方法可能会由于舍入而遭受严重误差；想象一下如果我们有 0.25 和 0.250000000001 作为概率！因此，除非在特殊情况下概率已知为良好且以有利于有理数运算的格式指定，否则最好不要尝试这种方法。

模拟有偏向硬币

我们从简单的随机基元（公平骰子）引申了一个简单但潜在效率低下的用于模拟不公平骰子算法。也许探索其他简单的随机基元可能会为解决这个问题提供不同的方法。

一个简单但有用的案例是使用随机生成器来模拟一个有偏向的硬币。假设有一枚硬币，正面朝上的概率为 $p_{heads}$，那么我们如何模拟投掷这枚有偏向的硬币呢？

我们之前形成的一个直觉是，我们可以将范围 $[0,1)$ 分成一系列小区间，以便当我们在该范围内选择一个随机值时，它最终会以与该区间的大小相等的概率进入某个区间。要使用在范围 $[0,1)$ 内均匀随机值来模拟有偏向的硬币，我们可以考虑将范围 $[0,1)$ 分割如下：

然后生成一个在范围 $[0,1)$ 内均匀随机的值，以确定它包含在哪个区间。幸运的是，由于只有一个分割点，确定点包含在哪个区间相当容易；如果值小于 $p_{heads}$，则硬币正面朝上，否则硬币反面朝上。以下是伪代码：

算法：模拟有偏向的硬币
生成一个在范围 $[0,1)$ 内均匀随机的值 $x$。
如果 $x<p_{heads}$，则返回 “正面”。
如果 $x≥p_{heads}$，则返回 “反面”。
由于我们可以在 $O(1)$ 时间内生成范围 $[0,1)$ 内的均匀随机值，并且也可以在 $O(1)$ 时间内进行实数比较，因此这个算法的运行时间为 $O(1)$。

用有偏向硬币模拟公平骰子

根据我们之前的讨论，我们知道可以使用一个公平骰子来模拟一个不公平的骰子，前提是我们愿意承受可能的空间使用增加。由于我们可以将有偏向的硬币看作一个有偏向的两面骰子，这意味着可以使用一个公平骰子来模拟一个有偏向的硬币。有趣的是，反过来也是可能的，我们可以使用有偏向的硬币来模拟公平骰子。这种构造方法简单而优雅，并且可以轻松推广到使用一组有偏向的硬币来模拟不公平的骰子。

有偏向的硬币的构造方法是基于硬币出现正面的概率将范围 $[0,1)$ 分为 “正面” 与 “反面” 区域 。我们已经在上面看到了类似的技巧，用于模拟一个公平的 $n$ 面骰子，通过将范围 $[0,1)$ 分成 $n$ 个大小相等的区域。例如，当掷一个四面骰子时，我们得到了这样的分区：

现在，假设我们有一组有偏向的硬币可供使用，我们有兴趣模拟掷一个公平骰子。我们可以考虑一种方法，从左到右穿越这些区间，每次都询问我们是否要停在当前所在的区间或继续前进。例如，假设我们想随机选择这四个区间中的一个。从最左边的区间开始，我们可以投掷一个有偏向的硬币，以确定我们是否应该停在这个区间或继续前进。由于我们想以概率 $\frac{1}{4}$ 均匀地选择所有这些区间，我们可以通过投掷一个出现正面概率为 $\frac{1}{4}$ 的有偏向硬币来实现这一点。如果它出现正面，我们就停在这个区间。否则，我们继续前进到下一个区间。

如果硬币出现反面，我们就会进入第二个区间，并且再次询问是否应选择这个区间或继续前进。最初，你可能认为我们应该再次投掷一个出现正面概率为 $\frac{1}{4}$ 的硬币来决定这一点，但实际上这是不正确的！看到这种推理的错误之处的一种方法是将其推向极端 - 如果在每个区间中我们都投掷一个出现正面概率为 $\frac{1}{4}$ 的硬币，有很小的概率每个区间都是反面，我们将拒绝每个区间。在某种程度上，当我们穿越这些区间时，我们需要增加硬币出现正面的概率。在极端情况下，如果我们最终进入最后一个区间，我们需要硬币以概率 $1$ 出现正面，因为如果我们拒绝了每一个之前的区间，正确的决定是停在最后一个区间。

为了确定我们的有偏向硬币在跳过第一个区间后出现正面的概率，注意如果我们确实跳过了第一个区间，那么只剩下三个区间。由于我们在掷一个公平骰子，我们希望这三个区间中的每一个都以概率 $\frac{1}{3}$ 被选择。因此，直观地看，我们应该让第二个硬币以概率 $\frac{1}{3}$ 出现正面。使用类似的论点，如果我们在第二个区间投掷反面，那么我们为第三个区间投掷的硬币应该以概率 $\frac{1}{2}$ 出现正面，而我们为最后一个区间投掷的硬币应该以概率 $1$ 出现正面。

这种直观的想法引导我们到以下算法。请注意，我们尚未论证这个算法是否正确；我们将在稍后进行讨论。

算法：使用有偏向硬币模拟公平骰子

1. 对于 $i$ 从 $0$ 到 $n-1$：
   1. 投掷一个有偏向硬币，出现正面的概率为 $\frac{1}{n-i}$ 。
   2. 如果硬币出现正面，返回 $i$。

这个算法简单，并且在最坏情况下的运行时间为 $O(n)$。但是我们如何知道它实际上是正确的呢？为了证明这一点，我们需要以下定理：

定理：上述算法对于任何 $i$ 的选择都以概率 $\frac{1}{n}$ 输出面 $i$。

证明：假设任意常数 $n ≥ 0$。我们通过强归纳来证明每个面都有概率 $\frac{1}{n}$ 被选择。

作为基本情况，我们证明骰子的第 $0$ 面以概率 $\frac{1}{n}$ 被选择。但是这是直接从算法得出的 - 如果一个有偏向硬币以概率 $\frac{1}{n}$ 出现正面，那么我们就选择第 $0$ 面，这意味着我们以概率 $\frac{1}{n}$ 选择它。

对于归纳步骤，假设对于第 $0, 1, 2, ..., k-1$ 面，这些面都以概率 $\frac{1}{n}$ 被选择，然后考虑第 $k$ 面被选择的概率。第 $k$ 面将被选择当且仅当前 $k$ 面都没有被选择，然后一个概率 $\frac{1}{n-k}$ 出现正面的硬币出现正面。因为前 $k$ 面每个都以概率 $\frac{1}{n}$ 被选择，而且只会选择其中一个面，所以前 $k$ 面中有一个被选择的概率为 $\frac{k}{n}$。这意味着算法不选择前 $k$ 面的概率为 $1 - \frac{k}{n} = \frac{n-k}{n}$。这意味着选择第 $k$ 面的概率为 $\frac{n-k}{n} \frac{1}{n-k} = \frac{1}{n}$，符合要求，完成了归纳。因此，每一面都以均匀随机的概率被选择。

当然，这个算法相当低效 - 我们可以使用之前的技术在 $O(1)$ 时间内模拟一个公平骰子的投掷！ - 但这个算法可以作为一种步骤，用于构建一个使用有偏向硬币模拟不公平骰子的相对高效的算法。

用有偏向硬币模拟不公平骰子

上述算法很有趣，因为它提供了一种使用一组硬币模拟骰子的简单框架。我们首先翻转一枚硬币来确定是选择骰子的第一面还是继续选择剩余的面。在这个过程中，我们必须小心地扩大剩余的概率。

让我们看看如何使用这种技术来模拟一个不公平的骰子。让我们使用之前的例子，概率分别为 $\frac{1}{2}$、 $\frac{1}{3}$、 $\frac{1}{12}$ 和 $\frac{1}{12}$。回顾一下，这将范围 $[0,1)$ 分为以下部分：

现在，让我们考虑如何使用有偏向硬币来模拟这个不公平的骰子。我们可以开始翻转一枚硬币，它有 $\frac{1}{2}$ 的概率出现正面，以确定是否应该输出第 $0$ 面。如果这枚硬币出现正面，那太好了！我们完成了。否则，我们需要再翻一枚硬币来确定是否选择下一面。与之前一样，尽管下一面出现的概率为 $\frac{1}{3}$，但我们不想翻一枚出现正面概率为 $\frac{1}{3}$ 的硬币，因为当我们没有选择第 $\frac{1}{2}$ 面时，就已经丢掉了一半的概率。实际上，由于一半的概率已经消失了，如果我们此时重新调整剩余的概率，我们会得到更新后的概率为 $\frac{2}{3}$、$\frac{1}{6}$、$\frac{1}{6}$。因此，第二枚硬币应该以概率 $\frac{2}{3}$ 出现正面。如果这枚硬币也出现反面，那么我们必须在两个 $\frac{1}{12}$ 面之间进行选择。由于在这一点上已经消失了 $\frac{5}{6}$ 的概率质量，我们将重新调整 $\frac{1}{12}$ 面的概率，以使每个面都有 $\frac{1}{2}$ 的概率出现正面，因此第三枚硬币将以概率 $\frac{1}{2}$ 出现正面。如果最后一枚硬币被翻转，那么它必须以概率 $1$ 出现正面，因为它是最后一个区间。

总结一下，硬币的概率如下：

第一次翻转：$\frac{1}{2}$
第二次翻转：$\frac{2}{3}$
第三次翻转：$\frac{1}{2}$
第四次翻转：$1$

尽管这些数字从直觉上看可能有意义，但要将其转化为算法，我们需要为选择概率提供正式的构建。思路如下 - 在每个步骤中，我们都记住了还有多少概率质量待使用。在开始之前，在翻转任何硬币之前，这个值为 $1$。翻转第一枚硬币后，它变为 $1 - p_{0}$。翻转第二枚硬币后，它变为 $1 - p_{0} - p_{1}$。更一般地说，翻转了 $k$ 枚硬币后，剩余的概率质量为 $1-\sum ^{k-1}_{i=0}p_{i}$ 。每当我们翻转一枚硬币来确定是否选择第 $k$ 个区间时，我们实际上是在翻转一枚硬币，其正面出现的概率等于剩余概率中由概率 $p_{k}$ 占据的部分，这个概率由 $\frac{p_{k}}{1-\sum ^{k-1}_{i=0}p_{i}}$ 给出。这为我们提供了使用一组有偏向硬币模拟有偏向骰子的以下算法（同样，我们稍后将证明其正确性和运行时间）：

算法：使用有偏向硬币模拟有偏向骰子

• 初始化：
  存储后续需要使用的概率 $p_{i}$。
• 生成：
  1. 设置 $mass=1$
  2. 对于 $i$ 从 $0$ 到 $n-1$：
     1. 翻转一枚有偏向硬币，正面出现的概率为 $\frac{p_{i}}{mass}$。
     2. 如果正面出现，返回 $i$。
     3. 否则，设置 $mass=mass−p_{i}$。

虽然这可能在直觉上有一定的道理，但从数学上来说，它是否正确呢？幸运的是，答案是肯定的，这要归功于对上面证明的一个推广，以下是具体的推广内容：

定理：上述算法对于任何选择的 $i$，以概率 $p_{i}$ 输出边 $i$。

证明：考虑任意固定的 $n≥0$。我们通过强归纳来证明 $n$ 个边中的每个边都具有概率 $p_{i}$ 被选择的特性。

作为基本情况，我们展示骰子的边 $0$ 具有概率 $p_{0}$ 被选择。如果第一次硬币翻转出现正面，我们选择边 $0$，这以概率 $\frac{p_{0}}{mass}$ 发生。由于 $mass$ 最初为 $1$，这个概率为 $\frac{p_{0}}{1}=p_{0}$，因此边 $0$ 以概率 $p_{0}$ 被选择，符合要求。

对于归纳步骤，假设边 $0、1、...、k-1$ 被选择的概率分别为 $p_{0}、p_{1}、...、p_{k-1}$，然后考虑边 $k$ 被选择的概率。边 $k$ 将被选择当且仅当前 $k$ 个边都未被选择，然后一枚正面概率为 $\frac{p_{k}}{mass}$ 的硬币翻转正面。因为前 $k$ 个边都以正确的概率被选择，而且只有一个边被选择，所以前 $k$ 个边中至少选择一个的概率由 $\sum ^{k-1}_{i=0}p_{i}$ 给出。这意味着算法不选择前 $k$ 个边的总概率为 $1-\sum ^{k-1}_{i=0}p_{i}$。现在，边 k 的硬币翻转正面的概率为 $\frac{p_{k}}{mass}$，经过 $k$ 次迭代后，我们可以快速归纳出 $mass=1-\sum ^{k-1}_{i=0}p_{i}$。这意味着选择边 $k$ 的总概率为 $(1-\sum ^{k-1}_{i=0}p_{i})\frac{p_{k}}{1-\sum ^{k-1}_{i=0}p_{i}}=p_{k}$，符合要求，完成了归纳证明。

现在，让我们考虑这个算法的运行时复杂度。我们知道初始化时间将是 $Θ(1)$（如果我们保留输入概率数组的浅拷贝），但如果我们要存储自己的版本数组（以防调用者以后想要更改它），则可能是 $Θ(n)$。实际生成骰子投掷结果可能在最坏情况下需要翻转 $Θ(n)$ 枚硬币，但在最佳情况下只需要翻转一次硬币。

然而，经过一些思考，可以清楚地看出输入分布会严重影响我们需要翻转多少次硬币。在绝对最佳情况下，我们有一个概率分布，其中所有的概率质量都集中在骰子的第一个面上，而其余的概率都为零。在这种情况下，我们只需要翻转一次硬币。在绝对最坏情况下，所有的概率质量都集中在骰子的最后一个面上，其他地方都是零，这种情况下，我们将需要翻转 n 枚硬币才能找到结果。

我们可以精确而数学地描述这个算法的预期硬币翻转次数。让我们考虑一个随机变量 $X$，表示在某个特定分布上执行此算法时翻转的硬币数量。也就是说，$ℙ[X=1]$ 表示算法在终止之前翻转一枚硬币的概率，$ℙ[X=2]$ 表示算法在终止之前翻转两枚硬币的概率，依此类推。在这种情况下，我们的算法的预期硬币翻转次数由随机变量 X 的期望值给出，表示为 $𝔼[X]$。根据定义，我们有以下公式：

$$𝔼[X]=\sum ^{n}_{i=1}i⋅ℙ[X=i]$$

那么 $ℙ[X=i]$ 的值是多少呢？嗯，算法将在选择了骰子的某一面后终止。如果选择了第 $0$ 面，它将翻一枚硬币来决定是否停在那里。如果选择了第 $1$ 面，它将翻两枚硬币 - 一枚来表示不想选择第 $0$ 面，另一枚来表示想要选择第 $1$ 面。更一般地说，如果算法选择了第 $i$ 面，它将翻 $i+1$ 枚硬币，其中 $i$ 枚用于决定不选择前面的 $i-1$ 面，一枚用于决定选择第 $i$ 面。结合我们知道的第 $i$ 面被选择的概率 $p_{i}$，这意味着：

$$𝔼[X] = \sum_{i=1}^{n} i \cdot ℙ[X = i] = \sum_{i=1}^{n} i \cdot p_{i-1} = \sum_{i=1}^{n} ((i-1)p_{i-1} + p_{i-1}) = \sum_{i=1}^{n} (i-1)p_{i-1} + \sum_{i=1}^{n} p_{i-1}$$

在最后的简化中，注意到第一项等于 $\sum ^{n-1}_{i=0}i⋅p_{i}$，这等于 $𝔼[p]$；即骰子投掷的期望结果！而第二项是 $1$，因为它是总概率的和。这意味着 $𝔼[X]=𝔼[p]+1$。也就是说，硬币翻转的预期次数等于骰子投掷的预期值再加上1！

算法	初始化时间（最好 - 最差）	生成时间（最好 - 最差）	内存使用（最好 - 最差）
用公平骰子模拟不公平骰子	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$	$Θ(1)$	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$
用有偏向硬币模拟不公平骰子	$Θ(n)$	$Θ(1)$ - $Θ(n)$	$Θ(n)$

归纳偏向硬币：模拟不公平骰子

在上面的例子中，我们能够有效地模拟有偏硬币，因为我们只需要考虑一个分割点。那么，当我们将这个思想推广到有偏骰子，其中面的数量可以任意大时，我们可以有多高效呢？

如果您注意到，有偏硬币正好与只有两个面的有偏骰子完全相同。因此，我们可以将有偏硬币看作是我们感兴趣解决的更一般问题的特殊情况。在解决有偏硬币问题时，我们将范围 $[0,1)$ 分成两个区域 - 一个表示"正面"，一个表示"反面" - 然后利用只有一个分割点的事实来找到我们所属的区域。如果我们有一个 $n$ 面的骰子，那么我们将有多个区域和因此多个分割点。例如，假设我们有一个七面的骰子，概率分别为 $\frac{1}{4}$ , $\frac{1}{5}$ , $\frac{1}{8}$, $\frac{1}{8}$, $\frac{1}{10}$, $\frac{1}{10}$, $\frac{1}{10}$。如果我们要将范围 $[0,1)$ 分成七个部分，我们将按如下方式进行分割：

请注意这些分割点的位置。第一个分区从 $0$ 开始，结束于 $\frac{1}{4}$。第二个分区从 $\frac{1}{4}$ 开始，结束于$\frac{1}{4}+\frac{1}{5}=\frac{9}{20}$。更一般地，如果概率是 $p_{0},p_{1},...,p_{n−1}$，则分割将分别在区间 $[0, p_{0})、[p_{0}, p_{0}+p_{1})、[p_{0}+p_{1}, p_{0}+p_{1}+p_{2})$ 等之间。也就是说，第 $i$ 个区间的范围是

$$\left[ \sum_{j=0}^{i-1} p_j, \sum_{j=0}^{i} p_j \right)$$

注意这两个值之间的差异是 $p_{i}$，因此该区域的总面积就是 $p_{i}$，正如所需的那样。

现在我们知道了分割点在哪里，如果我们选择一个在范围 $[0,1)$ 内均匀随机的值 $x$，我们该如何确定它落入哪个范围呢？以我们的有偏硬币算法为起点，一个想法是从第一个桶的端点开始，持续向上穿过分割，直到找到一个大于 $x$ 值的端点。如果这样做，我们将找到包含点 $x$ 的第一个桶，并且我们有了我们的值。例如，如果我们选择了随机值 $x=\frac{27}{40}$ ，我们将执行以下搜索：

从中我们可以得出骰子投掷了一个 $3$ ，如果从零开始计数。

这种线性扫描算法将为查找投掷的骰子面提供一个 $O(n)$ 时间的算法。然而，我们可以通过使用以下观察来显著提高这个运行时间：分割点的序列形成一个升序序列（因为我们总是添加越来越多的概率，其中没有一个可以小于零）。因此，我们正在尝试回答以下问题：给定一个升序值序列和一些测试点，找到范围内第一个严格大于测试点的值。这是使用二分搜索的绝佳场景！例如，以下是对上述数组进行二分搜索以找到值 $x=\frac{39}{40}$ 所属的桶的执行示例：

这为我们提供了一个将范围 $[0,1)$ 内的均匀随机值映射到投掷的骰子面的 $Θ(logn)$ 算法。此外，仅需要 $Θ(n)$ 预处理时间来构建端点表；我们只需计算所有概率的部分和直到顶部。

这个算法有时被称为轮盘赌选择，因为该算法使用类似于轮盘赌的技术来选择一个随机桶 - 把一个球扔进范围并看它落在哪里。在伪代码中，算法如下所示：

算法：轮盘赌选择

• 初始化：
  1. 分配大小为 $n$ 的数组 $A$
  2. 将 $A[0]$ 设置为 $p_{0}$。
  3. 对于每个概率 $i$ 从 $1$ 到 $n−1$：
     1. 将 $A[i]$ 设置为 $A[i−1]+p_{i}$
• 生成：
  1. 生成一个在范围 $[0,1)$ 内的均匀随机值 $x$
  2. 使用二分查找，找到大于 $x$ 的 $A$ 中最小的元素的索引 $i$。
  3. 返回 $i$

将这个算法与之前的算法比较，令人印象深刻：

算法	初始化时间（最好 - 最差）	生成时间（最好 - 最差）	内存使用（最好 - 最差）
用公平骰子模拟不公平骰子	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$	$Θ(1)$	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$
用有偏向硬币模拟不公平骰子	$Θ(n)$	$Θ(1)$ - $Θ(n)$	$Θ(n)$
轮盘赌选择	$Θ(n)$	$Θ(log n)$	$Θ(n)$

显然，我们现在有了一个比初始算法好得多的算法。最初似乎离散化概率是有希望的，但基于连续值和二分搜索的这种新方法似乎要好得多。然而，通过使用一组聪明的混合技术，我们仍然可以改进，我们将在接下来看到。

这个算法有一个有趣的细节，即虽然使用二分搜索可以保证随机生成的最坏情况时间为 $O(logn)$ ，但它也消除了更快查找的可能性；也就是说，生成时间现在也是 $Ω(logn)$。有没有可能做得更好？答案是肯定的。

假设我们不再使用简单的二分搜索来查找累积概率列表，而是使用二叉搜索树。例如，给定上面的概率集，我们可以构建如下的累积分布二叉搜索树：

现在，如果我们想模拟掷骰子，我们可以生成一个在范围 $[0,1)$ 内均匀分布的数字，然后查找它在这个二叉搜索树中的范围。由于这是一个平衡的二叉搜索树，最佳情况下的查找时间是 $O(1)$ ，最坏情况下的查找时间是 $O(logn)$。

然而，假设我们对概率分布了解更多，那么可能会有更好的方法。例如，假设我们的概率分布是 $\frac{99}{100}$,$\frac{1}{600}$,$\frac{1}{600}$,$\frac{1}{600}$,$\frac{1}{600}$,$\frac{1}{600}$,$\frac{1}{600}$。也就是说，概率分布非常倾斜，几乎所有的概率质量都集中在一个面上。我们可以为这些概率构建一个平衡的二叉搜索树，如下所示：

尽管这棵二叉搜索树是完美平衡的，但对于我们的应用来说，它并不是一个很好的二叉搜索树。由于我们知道 99 次中有 100 次随机值将位于范围 $[0,\frac{99}{100})$ 内，将该范围的节点存储在它当前的位置是没有意义的。事实上，这样做意味着几乎所有情况下我们都会进行两次不必要的与蓝色和黄色节点的比较。由于几乎所有情况下我们希望首先检查大范围的节点，因此通过以牺牲一些特殊情况来使平均情况显著改善是有道理的。如下所示：

现在，我们很有可能在第一次尝试中立即找到我们想要的桶并终止搜索。在非常不太可能的情况下，我们要查找的桶包含在小块 $(\frac{99}{100},1]$内，然后我们会优雅地降级到树的其余部分，这部分确实是平衡的。

更一般地说，我们想解决这个问题：

给定一组概率，找到这些概率的二叉搜索树，以最小化期望的查找时间。

幸运的是，这个问题已经被广泛研究，并被称为最优二叉搜索树问题。有许多解决这个问题的算法；已知可以使用动态规划在 $O(n^{2})$ 的时间内找到精确解，还存在可以找到近似解的良好线性时间算法。此外，自平衡二叉搜索树数据结构——伸展树，可以用来得到最优解的一个常数倍近似解。

有趣的是，这些优化二叉搜索树的最佳情况行为出现在概率分布极不均衡的情况下，因为我们可以将包含绝大部分概率质量的节点移到树的根附近，而它们的最坏情况出现在分布均衡的情况下，因为在这种情况下树必须宽而浅。这与早期使用公平骰子模拟有偏骰子的算法相反！

在最佳情况下，我们有一个有偏骰子，其中一面总是出现（即概率为 $1$，而其他每一面的概率为 $0$ ）。这是对我们早期示例的极端夸张，但会导致搜索总是在一次查找后终止。在最坏情况下，所有的概率都是均匀的，我们必须执行标准的二叉搜索树查找。这给出了以下信息：

算法	初始化时间（最好 - 最差）	生成时间（最好 - 最差）	内存使用（最好 - 最差）
用公平骰子模拟不公平骰子	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$	$Θ(1)$	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$
用有偏向硬币模拟不公平骰子	$Θ(n)$	$Θ(1)$ - $Θ(n)$	$Θ(n)$
轮盘赌选择	$Θ(n)$	$Θ(log n)$	$Θ(n)$
最佳的轮盘赌轮选择	$O(n^{2})$	$Θ(1)$ - $O(logn)$	$Θ(n)$

掷飞镖

到目前为止，我们已经看到了两种在某种程度上帮助我们构建模拟有偏骰子算法的基本概念：公平骰子和有偏硬币。仅使用公平骰子，我们提出了一个（虽然不实际）模拟有偏骰子的算法，而从有偏硬币的直觉出发，我们能够创造一个快速的模拟有偏骰子的算法。是否可能将这两种方法结合起来，构建一个既基于公平骰子又基于有偏硬币的算法？答案是肯定的，事实上，得到的算法优于上述两种方法。

到目前为止，我们一直将范围 $[0,1)$ 和骰子面的概率视为一维范围。上述两种算法都是通过在范围 $[0,1)$ 中选择某个点，然后将其映射到长度对应某个概率的线段上工作的。我们使线段越长，被选择的概率就越高。但是，如果我们不是以一维方式思考，而是以二维方式思考呢？如果我们不是将概率 $p_{i}$ 看作线段的长度，而是将其看作矩形的面积呢？

让我们重新回到早期的例子，即概率为 $\frac{1}{2}$、 $\frac{1}{3}$、$\frac{1}{12}$ 、$\frac{1}{12}$ 的情况。让我们将这些概率视为宽度为 $w$（对于某个任意的 $w>0$ ）且高度为 $p_{i}$ 的矩形（因此总面积为 $w⋅p_{i}$ ）：

注意，这些矩形的总面积合在一起是 $w$ ，因为面积是：

$$\sum_{i=0}^{n-1} wp_i = w \sum_{i=0}^{n-1} p_i = w$$

现在，假设我们在这些矩形周围绘制一个边界框，它的宽度为 $4w$ （因为有四个矩形），高度为 $\frac{1}{2}$（因为最高的矩形高度为 $\frac{1}{2}$）：

我们可以将这个矩形看作被分成五个区域 - 四个区域对应不同的概率，一个区域代表未使用的空间。有了这个划分，我们可以将模拟骰子的随机投掷算法看作是一个飞镖游戏。假设我们向这个目标投掷一枚（完全均匀分布的）飞镖。如果它击中了未使用的空间，我们将移除飞镖并再次投掷，重复这个过程，直到我们击中其中一个矩形。由于较高的概率对应于较大的矩形，因此特定骰子面出现的概率越高，我们最终命中其矩形的概率就越高。实际上，如果我们在实际击中某个矩形的情况下进行条件概率计算，我们有以下结论：

$$ℙ[\text{命中骰子面 }i \text{ 对应的矩形} | \text{命中某个矩形}] = \frac{\text{矩形 } i \text{ 的面积}}{\text{总矩形面积}} = \frac{wp_i}{w} = p_i$$

换句话说，一旦我们最终用均匀随机的飞镖命中某个矩形，我们将以概率 $p_{i}$ 选择加载骰子的第 $i$ 个面的矩形，正是我们希望它具有的概率！换句话说，如果我们能够找到一种高效的方法来模拟在这个矩形上随机投掷飞镖，那么我们就可以高效地模拟投掷随机骰子。

有一种思路是将飞镖投向这个矩形，方法是选择两个均匀随机值在范围 $[0,1)$ 内，将它们缩放到适当的宽度和高度，然后检查飞镖下方的区域。然而，这会导致与早前在尝试确定在早前的情况下随机值位于哪个一维桶时遇到的相同问题。然而，我们可以有一系列非常美妙的观察，可以使我们确定我们命中的位置变得简单，甚至可以说是信手拈来。

一眼看上去，注意到我们已经表明这些矩形的宽度可以任意选择，只要它们都有相同的宽度即可。当然，高度取决于骰子面的概率。然而，如果我们通过某个正实数 $h$ 均匀缩放所有的高度，那么所有矩形的相对面积将是相同的。事实上，对于任何正实数 $h$，一旦它们的高度被 $h$ 缩放，所有矩形的总面积都是由以下公式给出的：

$$\sum_{i=0}^{n-1} whp_i = wh \sum_{i=0}^{n-1} p_i = wh$$

现在考虑选择任何单个矩形的概率，前提是我们至少击中了某个矩形。使用与之前类似的数学方法，我们得到以下结果：

$$ℙ[\text{命中骰子面 }i \text{ 对应的矩形} | \text{命中某个矩形}] = \frac{\text{矩形 } i \text{ 的面积}}{\text{总矩形面积}} = \frac{whp_i}{wh} = p_i$$

因此，实际上，在它们被线性和均匀地缩放的情况下，选择任何单个矩形的概率都不会改变。

由于我们可以选择任何我们喜欢的正缩放因子，如果我们缩放这些矩形，使边界框的高度始终为 $1$ 会怎样？由于边界框的高度由输入概率的最大值 $p_{i}$ 决定，我们可以首先将每个矩形的比例缩小 $\frac{1}{p_{max}}$，其中 $p_{max}$ 是所有输入概率的最大概率。这将使矩形的高度为$1$。同样，由于我们可以选择任意的矩形宽度，让我们选择宽度为 $1$ 。这意味着对于 $n$ 个概率，边界框的总宽度为 $n$ ，总高度为 $1$ 。如下所示：

现在，我们已经准备好思考如何在这个矩形上扔一个随机飞镖并确定我们命中了什么。关键的见解是，我们可以将矩形分解，使其不再由几个较小的矩形和一个奇形怪状的开放空间组成，而是将区域划分为 $2n$ 个矩形的集合，每个概率有两个矩形。如下所示：

请注意这个矩形是如何形成的。对于每个有偏向骰子的一面，我们有一列宽度为 $1$、高度为 $1$ 的空间，被切成了两个部分 - 一个“是”半空间对应于该面的矩形，一个“否”半空间对应于该列的其余部分。

现在，让我们考虑如何在这个矩形中投掷飞镖。一个完全均匀的飞镖在这个矩形中的投掷将具有 $x$ 和 $y$ 两个分量。在这里，$x$ 分量必须在范围 $[0,1)$ 内，对应于飞镖落在哪一列。$y$ 分量必须在范围 $[0,1)$ 内，对应于我们在该列中的高度。$x$ 分量的选择影响我们考虑的有偏向骰子的哪一面，$y$ 分量的选择对应于我们是否选择该面或重新投掷。但等一下 - 我们之前见过这两个想法！选择 $x$ 坐标，它对应于一列，相当于掷一个公平骰子来决定选择哪一列。选择y坐标对应于翻转有偏向的硬币来确定是选择该面还是重新投掷！这个观察非常重要，我们应该特别清楚地阐述一下：

在这个矩形中选择一个随机点等同于掷一个公平骰子和翻转一个有偏向的硬币。

实际上，这个结果可以以更加强大的方式来考虑。为了模拟有偏向的骰子，我们构建了一组有偏向的硬币，每个硬币对应于骰子的一面，然后掷一个公平骰子来决定翻转哪个硬币。根据骰子的结果，如果适当的硬币翻转为正面，我们选择对应于所选列的结果，如果硬币翻转为反面，我们重复这个过程。

让我们总结一下到目前为止的重要点。首先，这些矩形的尺寸如下 - 对于每一面，“是” 矩形的高度由 $\frac{p_{i}}{p_{max}}$ 给出，"否"矩形的高度由 $\frac{p_{max}-p_{i}}{p_{max}}$ 给出。这将矩形的总高度归一化为 $1$ 。其次，每个矩形的宽度为 $1$ ，尽管老实说这个值无关紧要。最后，我们的算法如下：在选择某个结果之前，掷公平骰子以确定我们位于哪一列（换句话说，选择翻转哪个有偏向的硬币）。接下来，翻转适当的有偏向的硬币。如果硬币翻转为正面，选择对应于所选列的结果。否则，重复这个过程。

算法：公平骰子/有偏硬币模拟的不公平骰子

• 初始化：

1. 找到 $p_{i}$ 的最大值；记为 $p_{max}$ 。
2. 分配一个长度为 $n$ 的数组 $Coins$ ，对应于每一行中“是”矩形的高度。
3. 对于从 $0$ 到 $n-1$ 的每个概率i：
   1. 设置 $Coins[i]=\frac{p_{i}}{p_{max}}$

• 生成：

1. 直到找到一个值为止：
   1. 掷一个公平的 $n$ 面骰子并得到一个在范围 $[0,n)$ 内的索引i。
   2. 翻转一个有偏硬币，它以概率 $Coins[i]$ 出现正面。
   3. 如果硬币翻转为正面，返回 $i$

让我们分析一下这个算法的复杂性。在初始化步骤中，查找最大概率需要 $O(n)$ 的时间，然后分配和填充数组 $Coins$ 需要额外的 $O(n)$ 时间，因此总初始化时间为 $O(n)$ 。在生成步骤中，在最好的情况下，我们可能会在第一次投掷时翻转正面，以 $O(1)$ 时间终止。但是期望需要多少次迭代？为了找到这个值，让我们计算在单次迭代之后我们实际上选择某一面的概率。由于硬币的正面概率不都相同，这将取决于实际选择的硬币。幸运的是，由于我们以相同的概率（即 $\frac{1}{n}$ ）选择每个硬币，数学变得更容易。而且，由于我们只翻转了一个硬币，每个硬币被选择和出现正面的事件都是互斥的，因此某个硬币被选择、翻转并出现正面的总概率等于选择每个个体硬币并使该个体硬币出现正面的概率之和。由于我们知道选择面i的概率由 $\frac{p_{i}}{p_{max}}$ 给出，所以某一面被选择的总概率可以表示为：

$$\sum ^{n-1}_{i=0}\left( \dfrac{1}{n}\dfrac{p_{i}}{p_{\max }}\right) = \frac{1}{n}\sum ^{n-1}_{i=0} \frac{p_{i}}{p_{\max }}\ = \frac{1}{n · p_{max}}\sum ^{n-1}_{i=0} p_{i} = \frac{1}{n⋅p_{max}}$$

如果这是在任何一个迭代中选择某个硬币的概率，那么可能发生的期望迭代次数是该分数的倒数，即 $n⋅p_{max}$。但这到底意味着什么呢？这很大程度上取决于 $p_{max}$ 的选择。在一个极端情况下，$p_{max}$可能等于 $1$ （也就是说，骰子每次都是相同的）。在这种情况下，期望迭代次数等于 $n$ ，这意味着从期望上讲，我们需要投掷 $n$ 次公平的骰子。这是有道理的，因为我们只有在选择总是出现正面的那一面的有偏硬币时，才有可能选择一面，因为其他每一面都有一个硬币，根本不会出现正面。另一方面，在另一个极端，$p_{max}$ 的最小值为 $\frac{1}{n}$，因为如果小于这个值，所有面的总概率将小于 $1$ 。如果$p_{max}=\frac{1}{n}$，那么期望的翻转次数为 $1$ 。这也是有道理的。如果 $p_{max}=\frac{1}{n}$，那么每一面被选择的概率都相同（即 $\frac{1}{n}$ ），因此当我们将每一面的概率归一化为 $1$ 时，每一面都有 $1$ 的概率被选择。因此，用于选择要翻转的硬币的骰子投掷将有效地确定结果，因为硬币总是出现正面，我们永远不需要重复。

有趣的是，期望的翻转次数仅取决于 $p_{max}$ 的值，而不取决于任何其他涉及的概率，但如果我们回到我们的图形直觉，这是有道理的。我们投掷飞镖的矩形的总面积总是 $n$ ，因为我们将高度归一化为 $1$ 。此外，表示“是”答案的矩形的总面积是 $\frac{1}{p_{max}}$ ，因为每个矩形的宽度为 $1$ ，高度通过乘以 $\frac{1}{p_{max}}$ 来归一化。这意味着 “是” 矩形的总面积与整个矩形的总面积之比为$\frac{1}{n⋅p_{max}}$ 。换句话说，“否” 矩形占用的空间仅取决于 $p_{max}$ 的值。它可以分散或分布在我们选择的任何位置，但最终其面积是相同的，某个飞镖击中它的概率不受其分布方式的影响。

将这个算法与其他算法进行比较，我们可以得到以下信息：

算法	初始化时间（最好 - 最差）	生成时间（最好 - 最差）	内存使用（最好 - 最差）
用公平骰子模拟不公平骰子	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$	$Θ(1)$	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$
用有偏向硬币模拟不公平骰子	$Θ(n)$	$Θ(1)$ - $Θ(n)$	$Θ(n)$
轮盘赌选择	$Θ(n)$	$Θ(log n)$	$Θ(n)$
最佳的轮盘赌轮选择	$O(n^{2})$	$Θ(1)$ - $O(logn)$	$Θ(n)$
公平骰子+有偏硬币 模拟不公平骰子	$Θ(n)$	$Θ(1)$ - $Θ(n)$ (expected)	$Θ(n)$

在最佳情况下，这个算法比上面的二分搜索算法要好，只需要一次硬币翻转。然而，它的最坏情况性能要差得多，指数级别的差。有没有可能消除这种最坏情况的性能呢？

别名方法

前面的技巧在最佳情况下表现出色，只需一次公平骰子的掷骰子和硬币翻转即可生成随机点数。然而，从期望值来看，它的最坏情况性能要差得多，可能需要线性数量的掷骰子和硬币翻转。这是因为与之前的技巧不同，该算法可能会“错过”，并不断迭代，直到做出决定。从图形上来看，这是因为它是通过在一个可能包含大量未分配给任何结果的空白空间的目标上投掷飞镖来工作的。如果有一种方法可以消除所有的空白空间，以使目标的每一部分都由对应于加载骰子的某一面的矩形覆盖，那么我们只需在其上投掷一枚飞镖并读取结果。

我们可能会有一个特别聪明的想法，那就是调整矩形的高度，使其与平均概率相匹配，而不是与最大概率相匹配。让我们考虑上面的示例。在这里，我们的概率是 $\frac{1}{2}$、 $\frac{1}{3}$、$\frac{1}{12}$ 、$\frac{1}{12}$ 。由于有四个概率，平均概率必须是 $\frac{1}{4}$。如果我们尝试将矩形的高度归一化为平均值 $\frac{1}{4}$，而不是最大值 $\frac{1}{2}$，会发生什么呢？让我们看看会发生什么。我们从高度等于初始概率的宽度为 $1$ 的矩形开始：

现在，我们将所有这些矩形缩放，使得概率为 $\frac{1}{4}$ 的高度为 $1$ 。这是通过将每个概率乘以四来实现的，得到以下设置：

此时，让我们在这个图像上面画一个 $1×4$ 的矩形。这将代表我们将要射击的目标：

正如你所看到的，这并不是很完美，因为 $\frac{1}{2}$ 和 $\frac{1}{3}$ 的矩形不太容易放进这个框里。但如果我们允许自己将矩形切割成更小的部分呢？也就是说，如果我们将 $\frac{1}{2}$ 矩形的一部分空间切割并移到 $\frac{1}{12}$ 矩形的上方的空白部分，会得到这种情况，虽然垂直的条子还有点突出，但不太明显：

现在，我们仍然有一些突出，但不会太多。完全消除突出的一种方法是将 $\frac{1}{2}$ 和 $\frac{1}{3}$ 矩形的额外部分移动到空白处，但实际上有一个更好的解决方案。让我们从第一列中移出足够多的 $\frac{1}{2}$ 条，填补第三列中的空隙，这将在第一列留下一个小空隙，但会填充另一个空隙：

最后，我们可以将第二列的额外部分藏入第一列，得到这个最终的矩形：

以下所述的方法具有多个出色的特性。首先，表示加载骰子每一面的矩形的总面积与原始面积保持不变；我们所做的只是将这些矩形切成片段并移动它们。这意味着只要原始矩形的面积按照原始概率分布成比例地分布，那么每一面骰子所分配的总面积都是相同的。其次，请注意这个新矩形中没有空白空间，这意味着每当我们向它投掷飞镖，我们都能保证会命中某个能给我们最终答案的东西，而不是需要再次投掷飞镖的空白空间。这意味着一次飞镖投掷足以生成我们的随机值。最重要的是，每一列最多只包含两个不同的矩形。这意味着我们可以保留之前的直觉——我们掷骰子来确定要投掷哪个有偏硬币，然后投掷硬币。这一次不同之处在于硬币投掷的含义。正面朝上的投掷意味着我们选择骰子的一面，而反面的投掷现在意味着我们应该选择骰子的另一面（而不是再次投掷）。

在高层次上，别名方法的工作原理如下。首先，我们创建代表骰子不同概率的矩形。接下来，我们将这些矩形切成片段并重新排列它们，以便完全填充一个矩形目标，使得每一列都有固定的宽度，并包含来自加载骰子最多两个不同面的矩形。最后，我们通过将飞镖随机投掷到目标上来模拟掷骰子，这可以通过结合公平的骰子和有偏硬币来实现。

但我们怎么知道能够以一种方式将原始矩形切分开，以使每列最多包含两种不同的概率呢？这似乎不是显而易见的，但令人惊讶的是，这总是可能的。而且，不仅可以将矩形切成片段，使得每列最多包含两个不同的矩形，还可以以一种方式进行切割，其中每列中的一个矩形是最初放置在该列的矩形。如果你注意到，在上面的矩形重新排列中，我们总是切割一个矩形的一部分并将其移到另一列，从未完全将一个矩形从其原始列中移除。这意味着最终排列中的每一列都将包含与最初分配给该列的概率相对应的某个矩形，以及（可选）从其他列中拉出的第二个矩形。这第二个矩形通常被称为该列的别名，因为该列的剩余概率被用作某个其他概率的“别名”。这里使用“别名”一词导致了“别名方法”的名称。

在我们进一步讨论总是可以以这种方式分配概率之前，让我们简要概述一下这个算法实际是如何工作的。由于最终排列的每一列始终包含该列原始矩形的一些（可能是零高度的）部分，为了存储占据一列的（可能的）两个不同的矩形，别名方法的实现通常通过存储两个不同的表格来工作：一个概率表 $Prob$ 和一个别名表 $Alias$ 。这两个表格都具有大小 $n$ 。概率表存储了每一列内原始矩形将被选择的概率，而别名表存储了该列中包含的第二个矩形的标识（如果有的话）。这样，生成骰子的随机投掷可以如下进行。首先，使用一个公平的骰子，在随机情况下选择某一列 $i$。接下来，根据硬币概率 $Prob[i]$ 来投掷一个随机硬币。如果硬币正面朝上，输出骰子掷出 $i$ ，否则输出骰子掷出 $Alias[i]$ 。例如，这是上述配置的概率表和别名表：

（这里原文有一个前端可交互模拟，可以前往原文使用）

证明别名表存在

我们现在需要正式证明，总是可以构建上面的 $Alias$ 和 $Prob$ 表。为了证明这总是可能的，我们需要展示可以做到以下几点：

• 为每个概率 $p_{i}$ 构造 $(n⋅p_{i})×1$的矩形，
• 将它们水平切割成小块，
• 将它们分布到 $n$ 列中，
  • 使得每一列的高度为 $1$ ，
  • 没有列包含超过两个矩形，
  • 某个与第 $i$ 面对应的矩形被放置在第 $i$ 列中。

在进入总是可以这样做的证明之前，让我们通过一个例子来详细讨论。假设我们有和之前一样的四个概率 $\frac{1}{2}$、 $\frac{1}{3}$、$\frac{1}{12}$ 、$\frac{1}{12}$ 。这是一个包含四个概率（ $k=n=4$ ）的集合，它们的总和为 $1=\frac{4}{4}$。虽然我们已经通过实验看到了如何填写别名表，但让我们尝试通过更明确地从完全空的表格开始然后进行填充来更详细地走一遍这个构建过程。我们首先将所有这些概率乘以四，得到以下概率和这个空表格：

现在，请注意我们需要分配的四个矩形中，有两个（ $\frac{1}{3}$ , $\frac{1}{3}$ ）小于 $1$ 。这意味着它们无法完全填满一列，需要一些其他概率来填补剩余部分。让我们选择其中一个（比如，黄色的）并将其放入适当的列中：

现在，我们需要想办法弥补列顶部的差异。为了做到这一点，我们将注意到还没有分配的两个矩形的高度大于 $1$（即 $2$ 和 $\frac{4}{3}$ ）。让我们任意选择其中一个；在这里，让我们使用 $\frac{4}{3}$ 。然后，我们分配足够多的 $\frac{4}{3}$ 到列中，以完全填满它；这最终使用了$\frac{4}{3}$ 中的 $\frac{2}{3}$ ，如下所示：

现在，请注意我们的设置是什么样的。我们现在有三个总面积为 $3$ 的矩形和三个空列，所以似乎应该可以将这些矩形分配到三个列中。为了这样做，我们将使用与以前相同的直觉。请注意，至少有一个矩形的高度小于 $1$ ，因此我们会任意选择一个（假设我们拿到了 $\frac{4}{3}$ 的矩形）并将其放入其列中：

现在，我们需要填满列，所以我们将选择一个至少为 $1$ 的概率，并用它来填满列的其余部分。这里只有一个选择（使用 $2$ ），所以我们将从 $2$ 中取出 $\frac{1}{3}$ 并将其放在列的顶部：

现在我们只剩下两个高度总和为 $2$ 的矩形。我们现在重复这个过程，找到一个高度最多为 $1$ 的矩形（这里是 $\frac{1}{3}$ ）并将其放入相应的列中：

然后，我们找到一个高度至少为 $1$ 的矩形，以填满该列，我们唯一的选择是 $\frac{5}{3}$：

现在我们只剩下一个面积为 $1$ 的矩形，所以我们可以通过将该矩形放在自己的列中来完成构建：

大功告成！我们已经完成了表格的填充。

注意这个构建的一般模式如下：

1. 找到高度最多为 $1$ 的矩形，并将其放入自己的列中，将 $Prob$ 表设置为该矩形的高度。
2. 找到至少为 $1$ 的高度的矩形，并使用它来填充列，将 $Alias$ 表设置为表示矩形对应的骰子面。

我们能证明这种一般构建方法总是可行的吗？也就是说，我们不会以这种方式分配概率时陷入"僵局"吗？幸运的是，答案是肯定的。其背后的直觉是，我们已经按比例缩放了所有的概率，使新概率的平均值现在为 $1$（因为原来的平均值是 $\frac{1}{n}$ ，我们将所有东西都乘以 $n$ ）。我们知道所有缩放后的概率的最小值不能大于平均值，所有缩放后的概率的最大值不能小于平均值，所以当我们刚开始时，总会至少有一个元素小于等于 $1$（即，缩放后的概率中最小的一个）和一个元素大于等于 $1$（即，缩放后的概率中最大的一个）。我们可以将这些元素配对。但是当我们一旦移除了这两个元素后呢？嗯，当我们这样做时，我们会从总和中移除一个概率，并将缩放后的概率的总和减小 $1$ 。这意味着新的平均值没有改变，因为平均缩放后的概率为 $1$ 。然后，我们可以一遍又一遍地重复这个过程，直到最终我们将所有元素配对完毕。

我们可以将这个论点形式化为以下定理：

定理：给定高度为 $h_{0},h_{1},...,h_{k−1}$ 的 $k$ 个宽度为 $1$ 的矩形，使得 $\sum_{i=0}^{k-1} h_i = k$ ，存在一种切割这些矩形并将它们分配到 $k$ 个高度为 $1$ 的列中的方法，使得每个列包含至多两个不同的矩形，第 $i$ 列至少包含第 $i$ 个矩形的一部分。

证明：采用归纳法。首先考虑基本情况，如果 $k=1$ ，则只有一个矩形，其高度必须为 $1$ 。因此，我们可以将其分配到第 $0$ 列。因此，每个列的高度都是 $1$ ，包含至多两个矩形，第 $0$ 列至少包含第 $0$ 个矩形的一部分。

对于归纳步骤，假设对于某个自然数 $k$ ，定理成立，并考虑任意 $k+1$ 个宽度为$1$ 且高度为 $h_{0}，h_{1}，...，h_{k}$ 的矩形，使得 $\sum_{i=0}^{k} h_i = k+1$ 。首先，我们要声明存在某个高度 $h_{l}$ ，满足 $h_{l}≤1$ ，并且存在不同的高度 $h_{g}$（其中 $l≠g$ ），满足 $h_{g}≥1$ 。要理解这一点，假设为了推导矛盾，不存在 $h_{l}$，使得 $h_{l}≤1$ ；这将意味着对于范围 $0≤i≤k$ 内的所有自然数 $i$ ，$h_{i}>1$ 。但是，这将导致$k+1=\sum_{i=0}^{k} h_i>\sum_{i=0}^{k} 1 = k+1$，显然是不可能的。因此，存在某个索引 $l$ ，使得 $h_{l}≤1$ 。现在，假设为了推导矛盾，不存在其他高度 $h_{g}$（其中 $l≠g$ ），使得 $h_{g}≥1$ 。然后我们必须有每个其他 $h_{g}<1$ ，这将（通过类似的逻辑）意味着 $\sum_{i=0}^{k} h_i<k+1$ ，这与已知条件矛盾。因此，我们有 $h_{l}≤1$和 $h_{g}≥1$ 。

现在，考虑以下构造。将 $h_l$ 放入第 $l$ 列，并用矩形 $h_{g}$的一部分填充第l列中剩余的 $1−h_{l}$ 空间（这样的空间必定存在，因为 $0≤1−h_{l}≤1$ ，且 $h_{g}≥1$ ）。这将完全填满该列。现在我们剩下 $k$ 个不同矩形的集合，其总和为 $k$ ，因为我们从矩形中移除了总面积为 $1$（其初始总和为 $k+1$ ）。此外，我们已经完全填满了第 $l$ 列，因此我们不会再尝试将任何矩形的一部分放在那里。因此，根据归纳假设，我们可以将其余的 $k$ 个矩形分配到 $k$ 个列中，同时满足上述条件。结合现在已经填满的第 $l$ 列，这意味着我们有一种方法可以填满所有列并满足约束条件。这完成了归纳证明。

这是一个构造性的证明，它不仅说明我们总是可以构建别名表，还说明了上述算法，即找到一个高度最多为1的矩形并将其与一个高度至少为1的矩形配对，总是会成功。从这里开始，我们可以开始设计越来越快的计算别名表的算法。

生成别名表

根据我们上面提到的内容，我们可以得到一个相当不错的算法来使用别名方法模拟加载的骰子投掷。初始化阶段通过反复扫描输入的概率，找到至多为 $1$ 的值和至少为 $1$ 的值，将它们组合在一起填充一个列：

算法：朴素别名方法

• 初始化：
  1. 将每个概率 $p_{i}$ 乘以 $n$。
  2. 创建大小为 $n$ 的数组 $Alias$ 和 $Prob$。
  3. 对于 $j=1$ 到 $n−1$：
     1. 找到一个满足 $p_{l}≤1$ 的概率 $p_{l}$。
     2. 找到一个满足 $p_{g}≥1$ 且 $g ≠ l$ 的概率 $p_{g}$。
     3. 设置 $Prob[l]=p_{l}$。
     4. 设置 $Alias[l]=g$。
     5. 从初始概率列表中删除 $p_{l}$。
     6. 设置 $p_{g}:=p_{g}-(1−p_{l})$。
  4. 找到最后一个剩余的概率 $i$，它的权重必定为 $1$。
  5. 设置 $Prob[i]=1$。
• 生成：
  1. 从一个 $n$ 面骰子生成一个公平的骰子点数；将其称为面 $i$。
  2. 抛一枚带有概率 $Prob[i]$ 的有偏硬币。
  3. 如果硬币出现“正面”，则返回 $i$。
  4. 否则，返回 $Alias[i]$。

这个算法的生成步骤与上面描述的方法完全相同，并且运行时间为 $Θ(1)$ 。生成步骤需要多次迭代，具体描述如下。首先，我们需要花费 $Θ(n)$ 的时间将每个概率乘以 $n$ 的因子，并需要 $O(n)$ 的时间来分配这两个数组。内循环执行 $Θ(n)$ 次，在每次迭代中需要 $O(n)$ 的工作来扫描数组、删除其中一个数组元素并更新概率。这总共需要 $O(n^{2})$ 的初始化工作。如果我们将这个算法放在上下文中考虑，会得到以下结果：

算法	初始化时间（最好 - 最差）	生成时间（最好 - 最差）	内存使用（最好 - 最差）
用公平骰子模拟不公平骰子	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$	$Θ(1)$	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$
用有偏向硬币模拟不公平骰子	$Θ(n)$	$Θ(1)$ - $Θ(n)$	$Θ(n)$
轮盘赌选择	$Θ(n)$	$Θ(log n)$	$Θ(n)$
最佳的轮盘赌轮选择	$O(n^{2})$	$Θ(1)$ - $O(logn)$	$Θ(n)$
公平骰子+有偏硬币 模拟不公平骰子	$Θ(n)$	$Θ(1)$ - $Θ(n)$ (expected)	$Θ(n)$
朴素别名方法	$O(n^{2})$	$Θ(1)$	$Θ(n)$

与其他高效的模拟技术相比，这种朴素的别名方法具有较高的初始化成本，但随后可以非常高效地模拟掷骰子。如果我们能够以更低的初始化成本（比如 $O(n)$）来实现，那么这种技术将明显优于这里使用的所有其他技术。

为了减少初始化成本的一种简单方法是在进行过程中使用更好的数据结构来存储概率值。在朴素版本中，我们使用一个未排序的数组来保存所有的概率值，这意味着要花费 $O(n)$ 的工作来定位我们想要的两个概率值。更好的替代方法是使用平衡二叉搜索树来保存这些值。这样，我们可以在 $O(logn)$ 的时间内找到值 $p_{g}$ 和 $p_{l}$ ，方法是在树中找到最大值和最小值。删除 $p_{l}$ 可以在 $O(logn)$ 的时间内完成，更新 $p_{g}$ 的概率值也可以在 $O(logn)$ 的时间内完成，方法是从树中移除它并重新插入它。这给出了以下算法：

算法：别名方法

• 初始化：
  1. 创建大小为 $n$ 的数组 $Alias$ 和 $Prob$。
  2. 创建一个平衡二叉搜索树 $T$
  3. 对于每个概率值 $i$ ，将 $n⋅p_{i}$ 插入 $T$ 中
  4. 对于 $j=1$ 到 $n−1$：
     1. 找到找到并移除 $T$ 中的最小值；记为 $p_{l}$。
     2. 找到找到并移除 $T$ 中的最大值；记为 $p_{g}$。
     3. 设置 $Prob[l]=p_{l}$。
     4. 设置 $Alias[l]=g$。
     5. 设置 $p_{g}:=p_{g}-(1−p_{l})$。
     6. 将 $p_{g}$ 插入 $T$ 中
  5. 找到最后一个剩余的概率 $i$，它的权重必定为 $1$。
  6. 设置 $Prob[i]=1$。
• 生成：
  1. 从一个 $n$ 面骰子生成一个公平的骰子点数；将其称为面 $i$。
  2. 抛一枚带有概率 $Prob[i]$ 的有偏硬币。
  3. 如果硬币出现“正面”，则返回 $i$。
  4. 否则，返回 $Alias[i]$。

现在，我们的算法初始化速度快得多。创建 $Alias$ 和 $Prob$ 仍然需要 $O(n)$ 的时间，将概率添加到平衡二叉搜索树 $T$ 中将需要 $Θ(nlogn)$ 的时间。从那里开始，我们执行 $Θ(n)$ 次填充表的迭代，每次迭代需要 $O(logn)$ 的工作。这导致了初始化的总运行时间为 $O(nlogn)$，如下所示：

算法	初始化时间（最好 - 最差）	生成时间（最好 - 最差）	内存使用（最好 - 最差）
用公平骰子模拟不公平骰子	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$	$Θ(1)$	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$
用有偏向硬币模拟不公平骰子	$Θ(n)$	$Θ(1)$ - $Θ(n)$	$Θ(n)$
轮盘赌选择	$Θ(n)$	$Θ(log n)$	$Θ(n)$
最佳的轮盘赌轮选择	$O(n^{2})$	$Θ(1)$ - $O(logn)$	$Θ(n)$
公平骰子+有偏硬币 模拟不公平骰子	$Θ(n)$	$Θ(1)$ - $Θ(n)$ (expected)	$Θ(n)$
朴素别名方法	$O(n^{2})$	$Θ(1)$	$Θ(n)$
别名方法	$O(log n)$	$Θ(1)$	$Θ(n)$

然而，有一种算法比这种方法运行得更快。这个算法非常简单，也许是实现别名方法的所有算法中最简洁的。这个算法最初是在Michael Vose的论文《使用给定分布生成随机数的线性算法》中描述的，已经成为实现别名方法的标准算法。

Vose算法的思想是维护两个工作列表，一个包含高度小于1的元素，另一个包含高度至少为1的元素，并反复配对每个工作列表的第一个元素。在每次迭代中，我们使用来自“小”工作列表的元素，并有可能将来自“大”工作列表的元素的剩余部分移入“小”工作列表。该算法维护了几个不变性：

1. “小”工作列表中的元素都小于1。
2. “大”工作列表中的元素都至少为1。
3. 工作列表中元素的总和始终等于总元素数。

为了简化，每个工作列表不存储实际的概率值，而是存储指回原概率列表的某些指针，指示加载的骰子的哪一面被引用。根据这些不变性，以下是算法的实现：

算法：（不稳定的）Vose的别名方法算法

注意：该算法存在数值不准确的问题。稍后给出一个数值上更合理的算法。

• 初始化：

1. 创建大小为 $n$ 的数组 $Alias$ 和 $Prob$。
2. 创建两个工作列表 $Small$ 和 $Large$
3. 将每个概率 $p_{i}$ 乘以 $n$。
4. 对于每个概率值 $p_{i}$
   1. 如果 $p_{i} < 1$，将 $i$ 添加到 $Small$。
   2. 否则（$p_{i} ≥ 1$），将 $i$ 添加到 $Large$。
5. 循环直到 $Small$ 为空
   1. 弹出 $Small$ 的首元素；称为 $l$
   2. 弹出 $Large$ 的首元素；称为 $g$
   3. 设置 $Prob[l]=p_{l}$。
   4. 设置 $Alias[l]=g$。
   5. 设置 $p_{g}:=p_{g}-(1−p_{l})$。
   6. 如果 $p_{g} < 1$，将 $g$ 添加到 $Small$。
   7. 否则（$p_{g} ≥ 1$），将 $g$ 添加到 $Large$。
6. 循环直到 $Large$ 为空
   1. 弹出 $Large$ 的首元素；称为 $g$
   2. 设置 $Prob[g]=1$。

• 生成：

1. 从一个 $n$ 面骰子生成一个公平的骰子点数；将其称为面 $i$。
2. 抛一枚带有概率 $Prob[i]$ 的有偏硬币。
3. 如果硬币出现“正面”，则返回 $i$。
4. 否则，返回 $Alias[i]$

根据上述三个不变性质，该算法的前半部分（除了最后的循环之外的部分）应该是相对容易理解的：我们不断将 $Small$ 中的某个小元素与 $Large$ 中的某个大元素配对，然后将大元素的剩余部分添加到相应的工作列表中。算法中的最后一个循环需要解释一下。一旦我们耗尽了 $Small$ 列表中的所有元素，$Large$ 列表中至少会有一个元素剩余（因为如果每个元素都在Small中，那么元素的总和必须小于剩余元素的数量，违反了最后一个不变性质）。由于 $Large$ 中的每个元素至少为 $1$ ，并且由于 $Large$ 中的 $k$ 个元素的总和必须等于 $k$ ，这意味着 $Large$ 中的每个元素必须恰好等于 $1$，否则总和将会太大。因此，这个最后的循环将每个大元素的概率设置为 $1$，以便包含大元素的列都等于 $1$。

在这个算法中，工作列表的类型并不重要。Vose的原始论文使用堆栈作为工作列表，因为它们可以使用数组有效地实现，但如果愿意，我们也可以使用队列。然而，为了简单起见，我们将使用堆栈。

在对算法进行分析之前，让我们首先通过一个示例来跟踪看看它是如何工作的。我们考虑使用七个概率值的例子：$$\frac{1}{4}, \frac{1}{5}, \frac{1}{8}, \frac{1}{8}, \frac{1}{10}, \frac{1}{10}, \frac{1}{10}$$。为了突出该算法不对概率进行排序或要求它们按顺序排列的事实，让我们任意地对它们进行排序，如下所示：$$\frac{1}{8}, \frac{1}{5}, \frac{1}{10}, \frac{1}{4}, \frac{1}{10}, \frac{1}{10}, \frac{1}{8}。算法开始通过将这些元素添加到两个工作栈中，如下所示：

现在，我们将 $Small$ 栈的顶部放入其位置，将品红色的矩形放入最终位置：

现在，我们使用 $Large$ 栈的顶部（青色矩形）来填充列的其余部分。由于 $$ \frac{7}{4} - \frac{1}{8} = \frac{13}{8} ≥ 1 $$，我们将青色块留在 $Large$ 栈的顶部，如下所示

然后我们重复这个过程。我们将位于 $Small$ 栈顶部的矩形放入其所在的列，然后用 $Large$ 栈顶部的矩形补足差额：

再来一次：

当我们下次重复这个过程时，我们会发现虽然我们可以使用青色块来填充表格中的空白空间，但这样做会使青色块的高度小于 $1$。因此，我们将青色块移到 $Small$ 栈的顶部，如下所示：

现在，当我们处理 $Small$ 工作列表时，我们最终将青色块放在了它的位置上，然后使用黄色块填充了剩余的空白空间：

然后，我们处理 $Small$ 栈，将橙色块放置在其位置上，并用黄色块填充其顶部：

最后，由于 $Small$ 栈为空，我们将黄色块放入自己的列中，工作完成。

现在，我们对这些概率有一个良好的别名表。

一个切实可行的 Vose 版本

很不幸，上面的算法在其当前形式下并不是数值稳定的。在一个可以进行任意精度实数计算的理想机器上，这个算法可能没问题，但如果你尝试在使用IEEE-754双精度浮点数的计算机上运行它，它可能会完全失败。我们需要在继续之前解决两个不准确性的来源：

1. 确定概率应该属于 $Small$ 还是 $Large$ 工作列表的计算可能不准确。具体来说，通过因子 $n$ 缩放概率可能导致原本等于 $\frac{1}{n}$ 的概率略小于 $1$，从而被放入 $Small$ 列表而不是 $Large$ 列表。

2. 从较大的概率中减去适当的概率质量的计算不是数值稳定的，可能引入显著的舍入误差。这可能会导致本应属于 $Large$ 列表的概率被错误地放入 $Small$ 列表。

这两个因素的结合意味着，该算法可能会意外地将所有概率都放入 Small 工作列表而不是 Large 工作列表。因此，该算法可能会失败，因为它期望在 Small 工作列表非空时，Large 工作列表也非空。

幸运的是，修复这个问题并不特别困难。我们将更新算法的内部循环，以便在两个工作列表中的任何一个为空时终止循环，以避免意外查看 Large 工作列表中不存在的元素。其次，当一个工作列表为空时，我们将设置另一个工作列表中元素的剩余概率都为 1，因为从数学上讲，这只有在所有剩余的概率都精确等于 1 时才会发生。最后，我们将替换更新大概率的计算，采用稍微更稳定的计算方法。具体如下所示：

算法：Vose的别名方法算法

• 初始化：

1. 创建大小为 $n$ 的数组 $Alias$ 和 $Prob$。
2. 创建两个工作列表 $Small$ 和 $Large$
3. 将每个概率 $p_{i}$ 乘以 $n$。
4. 对于每个概率值 $p_{i}$
5. 如果 $p_{i} < 1$，将 $i$ 添加到 $Small$。
6. 否则（$p_{i} ≥ 1$），将 $i$ 添加到 $Large$。
7. 循环直到 $Small$ 或 $Large$ 为空（$Largre$ 可能先为空）
8. 弹出 $Small$ 的首元素；称为 $l$。
9. 弹出 $Large$ 的首元素；称为 $g$。
10. 设置 $Prob[l]=p_{l}$。
11. 设置 $Alias[l]=g$。
12. 设置 $p_{g}:=(p_{g}+p_{l})-1$ （这是一个数值上更稳定的选项）。
13. 如果 $p_{g} < 1$，将 $g$ 添加到 $Small$。
14. 否则（$p_{g} ≥ 1$），将 $g$ 添加到 $Large$。
15. 循环直到 $Large$ 为空
16. 弹出 $Large$ 的首元素；称为 $g$。
17. 设置 $Prob[g]=1$。
18. 循环直到 $Small$ 为空 （这只可能是由于数值不稳定造成的）
    1. 弹出 $Small$ 的首元素；称为 $l$。
    2. 设置 $Prob[g]=l$。

• 生成：

1. 从一个 $n$ 面骰子生成一个公平的骰子点数；将其称为面 $i$。
2. 抛一枚带有概率 $Prob[i]$ 的有偏硬币。
3. 如果硬币出现“正面”，则返回 $i$。
4. 否则，返回 $Alias[i]$

现在剩下的工作就是分析算法的复杂性。填充工作列表总共需要 $Θ(n)$ 的时间，因为我们将每个元素添加到工作列表中的一个。内部循环总共需要 $Θ(1)$ 的工作，因为它需要从工作列表中移除两个元素，更新两个数组，并将一个元素添加回工作列表。它不能执行超过 $O(n)$ 次，因为每次迭代都通过消除较小的概率来减少工作列表中的元素数量（总共）。最后两个循环每个都可以执行最多 $O(n)$ 次，因为 $Large$ 和 $Small$ 工作列表中的元素最多有 $O(n)$ 个。这给出了总运行时间 $Θ(n)$，如下所示，这是我们能够获得的最好的运行时间：

算法	初始化时间（最好 - 最差）	生成时间（最好 - 最差）	内存使用（最好 - 最差）
用公平骰子模拟不公平骰子	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$	$Θ(1)$	$Θ(n)$ - $O(∏_{i=0}^n{d_i})$
用有偏向硬币模拟不公平骰子	$Θ(n)$	$Θ(1)$ - $Θ(n)$	$Θ(n)$
轮盘赌选择	$Θ(n)$	$Θ(log n)$	$Θ(n)$
最佳的轮盘赌轮选择	$O(n^{2})$	$Θ(1)$ - $O(logn)$	$Θ(n)$
公平骰子+有偏硬币 模拟不公平骰子	$Θ(n)$	$Θ(1)$ - $Θ(n)$ (expected)	$Θ(n)$
朴素别名方法	$O(n^{2})$	$Θ(1)$	$Θ(n)$
别名方法	$O(log n)$	$Θ(1)$	$Θ(n)$
Vose 的别名方法	$Θ(n)$	$Θ(1)$	$Θ(n)$

总结一下

哇！我们已经涵盖了很多内容！我们探讨了几种模拟加载骰子的方法，从非常简单的技术开始，到最后得出了非常快速和高效的算法。每种方法展示了不同的技术，我认为最终版本（Vose的别名方法）是我曾经遇到的最有趣和优雅的算法之一。

如果您对查看Vose的别名方法的代码感兴趣，包括由于数值不精确性而在实践中出现的一些复杂性的简要总结，我在有趣代码存档中提供了别名方法的Java实现。

如果您有任何问题、意见或更正意见，请随时通过[email protected]与我联系。

希望这对您有所帮助！