在一个范围内生成无偏随机整数的最佳算法是什么？答案

【问题标题】：What is the optimal algorithm for generating an unbiased random integer within a range?在一个范围内生成无偏随机整数的最佳算法是什么？
【发布时间】：2012-07-30 07:45:34
【问题描述】：

在这个 StackOverflow 问题中：

Generating random integer from a range

接受的答案建议使用以下公式在给定的min 和max 之间生成一个随机整数，其中min 和max 包含在该范围内：

output = min + (rand() % (int)(max - min + 1))

但它也这么说

这仍然略微偏向于较低的数字......它也是可以扩展它以消除偏见。

但它没有解释为什么它偏向于较低的数字或如何消除这种偏见。所以，问题是：这是在（有符号）范围内生成随机整数的最佳方法，而不依赖任何花哨的东西，只是rand() 函数，如果它是最优的，如何消除偏差?

编辑：

我刚刚针对浮点外推法测试了@Joey 建议的while-loop 算法：

static const double s_invRandMax = 1.0/((double)RAND_MAX + 1.0);
return min + (int)(((double)(max + 1 - min))*rand()*s_invRandMax);

查看有多少“球”均匀地“落入”并分布在多个“桶”中，一个测试用于浮点外推，另一个测试用于while-loop 算法。但是结果会根据“球”（和“桶”）的数量而有所不同，所以我无法轻易选出获胜者。工作代码可以在this Ideone page 找到。例如，对于 10 个桶和 100 个球，浮点外推法与桶之间理想概率的最大偏差小于while-loop 算法（分别为 0.04 和 0.05），但对于 1000 个球，最大偏差为while-loop 算法较小（0.024 和 0.011），并且对于 10000 个球，浮点外推再次做得更好（0.0034 和 0.0053），依此类推，没有太大的一致性。考虑到没有一种算法能够始终如一地产生比其他算法更好的均匀分布的可能性，这让我倾向于浮点外推，因为它似乎比while-loop 算法执行得更快。那么选择浮点外推算法是否可以，否则我的测试/结论并不完全正确？

【问题讨论】：

它偏向于较小的数字，因为 rand() 产生一定范围内的数字（由 RAND_MAX 定义），而 RAND_MAX 通常不能被除数整除，所以从 0 到所有数字RAND_MAX % 除数 - 1 将有更高的机会被选中。
假设你有一个随机数生成器，它给出 0、1 或 2 的机会均等。如果你对它应用模 2 得到 0 或 1，你可以看到 0 有两倍被选中的几率。我想这就是你引用的陈述所暗示的。
eternallyconfuzzled.com/arts/jsw_art_rand.aspx
输入是整数，输出是整数 - 在转换期间转换为浮点数不会改变某些数字比其他数字更有可能的事实，尽管它可能会改变哪个这些数字将是。
Generating random integer from a range的可能重复

标签： c++ c random uniform

【解决方案1】：

问题是您正在执行模运算。如果RAND_MAX 可以被您的模数整除，这将没有问题，但通常情况并非如此。作为一个非常人为的示例，假设 RAND_MAX 为 11，您的模数为 3。您将获得以下可能的随机数和以下结果余数：

0 1 2 3 4 5 6 7 8 9 10
0 1 2 0 1 2 0 1 2 0 1

如您所见，0 和 1 的可能性略高于 2。

解决此问题的一个选项是拒绝抽样：通过禁止上面的数字 9 和 10，您可以使结果分布再次均匀。棘手的部分是弄清楚如何有效地做到这一点。在 Java 的 java.util.Random.nextInt(int) 方法中可以找到一个非常好的示例（我花了两天时间才理解为什么它有效）。

Java 的算法有点棘手的原因是它们避免了像乘法和除法这样的慢操作来进行检查。如果你不太在乎，你也可以用幼稚的方式来做：

int n = (int)(max - min + 1);
int remainder = RAND_MAX % n;
int x, output;
do {
  x = rand();
  output = x % n;
} while (x >= RAND_MAX - remainder);
return min + output;

编辑： 更正了上述代码中的一个栅栏错误，现在它可以正常工作了。我还创建了一个小示例程序（C#；为 0 到 15 之间的数字采用统一的 PRNG，并通过各种方式为 0 到 6 之间的数字构建一个 PRNG）：

using System;

class Rand {
    static Random r = new Random();

    static int Rand16() {
        return r.Next(16);
    }

    static int Rand7Naive() {
        return Rand16() % 7;
    }

    static int Rand7Float() {
        return (int)(Rand16() / 16.0 * 7);
    }

    // corrected
    static int Rand7RejectionNaive() {
        int n = 7, remainder = 16 % n, x, output;
        do {
            x = Rand16();
            output = x % n;
        } while (x >= 16 - remainder);
        return output;
    }

    // adapted to fit the constraints of this example
    static int Rand7RejectionJava() {
        int n = 7, x, output;
        do {
            x = Rand16();
            output = x % n;
        } while (x - output + 6 > 15);
        return output;
    }

    static void Test(Func<int> rand, string name) {
        var buckets = new int[7];
        for (int i = 0; i < 10000000; i++) buckets[rand()]++;
        Console.WriteLine(name);
        for (int i = 0; i < 7; i++) Console.WriteLine("{0}\t{1}", i, buckets[i]);
    }

    static void Main() {
        Test(Rand7Naive, "Rand7Naive");
        Test(Rand7Float, "Rand7Float");
        Test(Rand7RejectionNaive, "Rand7RejectionNaive");
    }
}

结果如下（粘贴到Excel中，并添加了单元格的条件着色，使差异更加明显）：

现在我修正了我在上述拒绝抽样中的错误，它可以正常工作（在它偏向 0 之前）。如您所见，float 方法并不完美，它只是以不同的方式分配有偏的数字。

【讨论】：

您是否认为一种明显的基于浮点的方法旨在将 RAND_MAX 范围映射到自定义范围（几乎）而无需分支，其 CPU 效率会低于多个分区 (%) 和分支while 循环暗示？
如果你改用(double)rand() / RAND_MAX * n，你会遇到同样的问题，只是不是较低的数字更有可能，而是你在整个范围内分配了一些偏差。使用该方法根本无法消除偏见。您仍然存在将 10 个输入数字均匀地拟合为 3 个输出数字的问题，这是不可能的。
我不明白distribute the bias somewhat over the whole range这句话——似乎自相矛盾。另外，ITYM ... evenly fitting *11* input numbers into 3 output numbers ...
例如，您可能会偏向 0 和 2，而不是像模数那样在前面设置更多偏向的数字。或者 1 和 2。这有点不可预测，因为它取决于浮点数的样子以及它们在哪里切换到另一个数字。你仍然有偏见。是的，这是 11 个数字。请原谅我无法再编辑的评论中的一个小错字。要点保持不变。
@Joey 请查看我的帖子的编辑内容，了解我使用this code 运行的一些测试。

【解决方案2】：

很容易看出为什么该算法会产生有偏差的样本。假设您的rand() 函数从集合{0, 1, 2, 3, 4} 返回统一整数。如果我想用它来生成随机位0 或1，我会说rand() % 2。集合{0, 2, 4} 给了我0，集合{1, 3} 给了我1——很明显，我对0 的采样率为60%，1 的采样率为40%，根本不统一！

要解决此问题，您必须确保所需的范围除以随机数生成器的范围，否则只要随机数生成器返回的数字大于最大的数字，就丢弃结果目标范围的可能倍数。

在上面的例子中，目标范围是2，适合随机生成范围的最大倍数是4，所以我们丢弃任何不在集合{0, 1, 2, 3}中的样本并再次滚动。

【讨论】：

【解决方案3】：

当随机数生成器 (RAND_MAX+1) 的输出数量不能被所需范围 (max-min+1) 整除时，就会出现问题。由于从随机数到输出会有一致的映射，因此某些输出将比其他输出映射到更多的随机数。这与映射如何完成无关 - 您可以使用模数、除法、转换为浮点数，无论您能想出什么巫术，基本问题仍然存在。

问题的严重性非常小，要求不高的应用程序通常可以忽略它而侥幸成功。范围越小，RAND_MAX越大，效果越不明显。

我采用了您的示例程序并对其进行了一些调整。首先我创建了一个特殊版本的rand，它的范围只有0-255，以更好地展示效果。我对rangeRandomAlg2 做了一些调整。最后我将“球”的数量更改为 1000000 以提高一致性。你可以在这里看到结果：http://ideone.com/4P4HY

请注意，浮点版本会产生两个紧密组合的概率，接近 0.101 或 0.097，两者之间没有任何关系。这就是行动中的偏见。

我认为将此称为“Java 的算法”有点误导 - 我确信它比 Java 古老得多。

int rangeRandomAlg2 (int min, int max)
{
    int n = max - min + 1;
    int remainder = RAND_MAX % n;
    int x;
    do
    {
        x = rand();
    } while (x >= RAND_MAX - remainder);
    return min + x % n;
}

【讨论】：

嗯，在阅读 Java 的源代码时，我发现了它们的实现（尽管有点不明显，我仍然觉得它很好）。当然，这不是此处任何答案中所写的内容，因为您从我的答案中复制的幼稚检查更容易理解。
我想知道它是否应该正确处理n大于RAND_MAX的情况，到目前为止它是一个无限循环。
嗯，这永远不会终止（或需要很长时间），所有可能的 PRNG？我想 RNG 最终必须返回 [0, RAND_MAX] 中的所有值才能使其始终有效。
@DesmondHume，不，它不能处理 n 大于 RAND_MAX 的情况——当某些值不可能时，你怎么能希望得到均匀分布？您可以检查该条件的代码。
@AmbrozBizjak，最坏的情况是当 n 为 RAND_MAX/2+1 时，平均会丢弃一半的随机数。当然，最坏的情况可能会更糟，但如果运行足够多的时间，它应该会回落到平均值。没有产生其范围内所有可能值的 RNG 将无法满足大多数随机定义。

【解决方案4】：

到目前为止，最简单的解决方案是std::uniform_int_distribution<int>(min, max)。

【讨论】：

【解决方案5】：

如果你真的想得到一个完美的生成器，假设你拥有的 rand() 函数是完美的，你需要应用下面解释的方法。

我们将创建一个随机数r，从0到max-min=b-1，然后很容易移动到你想要的范围，只需取r+min

我们将创建一个 b

程序：

以原始 RAND_MAX 大小取一个随机数 r，不进行任何截断
以 b 为底显示此数字
取这个数的前 m=floor(log_b(RAND_MAX)) 个数字，得到从 0 到 b-1 的 m 个随机数
按 min （即 r+min）将它们移动到所需的范围 (min,max) 中

由于 log_b(RAND_MAX) 不一定是整数，所以表示中的最后一位数字被浪费了。

原来只使用 mod (%) 的方法被

弄错了

(log_b(RAND_MAX) - floor(log_b(RAND_MAX)))/ceil(log_b(RAND_MAX))

你可能会同意这并不多，但如果你坚持要精确，那就是程序。

【讨论】：

【解决方案6】：

不失一般性，在[a,b]上生成随机整数的问题可以简化为在[0,s)上生成随机整数的问题。以下最近的出版物代表了从统一 PRNG 生成有界范围内随机整数的最新技术：

Daniel Lemire，“区间内的快速随机整数生成”。 ACM 翻译。模型。计算。 Simul. 29, 1, Article 3 (January 2019) (ArXiv draft)

Lemire 表明他的算法提供了无偏的结果，并且受到越来越受欢迎的非常快速的高质量 PRNG（例如 Melissa O'Neill 的PCG generators）的推动，展示了如何快速计算结果，几乎避免了缓慢的除法运算所有的时间。

下面的randint() 显示了他的算法的示例性 ISO-C 实现。在这里，我结合 George Marsaglia 的旧版KISS64 PRNG 来演示它。出于性能原因，所需的 64×64→128 位无符号乘法通常最好通过机器特定的内在函数或直接映射到适当硬件指令的内联汇编来实现。

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>

/* PRNG state */
typedef struct Prng_T *Prng_T;
/* Returns uniformly distributed integers in [0, 2**64-1] */
uint64_t random64 (Prng_T);
/* Multiplies two 64-bit factors into a 128-bit product */
void umul64wide (uint64_t, uint64_t, uint64_t *, uint64_t *);

/* Generate in bias-free manner a random integer in [0, s) with Lemire's fast
   algorithm that uses integer division only rarely. s must be in [0, 2**64-1].

   Daniel Lemire, "Fast Random Integer Generation in an Interval," ACM Trans.
   Model. Comput. Simul. 29, 1, Article 3 (January 2019)
*/
uint64_t randint (Prng_T prng, uint64_t s) 
{
    uint64_t x, h, l, t;
    x = random64 (prng);
    umul64wide (x, s, &h, &l);
    if (l < s) {
        t = (0 - s) % s;
        while (l < t) {
            x = random64 (prng);
            umul64wide (x, s, &h, &l);
        }
    }
    return h;
}

#define X86_INLINE_ASM (0)

/* Multiply two 64-bit unsigned integers into a 128 bit unsined product. Return
   the least significant 64 bist of the product to the location pointed to by
   lo, and the most signfiicant 64 bits of the product to the location pointed
   to by hi.
*/
void umul64wide (uint64_t a, uint64_t b, uint64_t *hi, uint64_t *lo)
{
#if X86_INLINE_ASM
    uint64_t l, h;
    __asm__ (
        "movq  %2, %%rax;\n\t"  // rax = a
        "mulq  %3;\n\t"         // rdx:rax = a * b
        "movq  %%rax, %0;\n\t"  // l = (a * b)<31:0>
        "movq  %%rdx, %1;\n\t"  // h = (a * b)<63:32>
        : "=r"(l), "=r"(h)
        : "r"(a), "r"(b)
        : "%rax", "%rdx");
    *lo = l;
    *hi = h;
#else // X86_INLINE_ASM
    uint64_t a_lo = (uint64_t)(uint32_t)a;
    uint64_t a_hi = a >> 32;
    uint64_t b_lo = (uint64_t)(uint32_t)b;
    uint64_t b_hi = b >> 32;

    uint64_t p0 = a_lo * b_lo;
    uint64_t p1 = a_lo * b_hi;
    uint64_t p2 = a_hi * b_lo;
    uint64_t p3 = a_hi * b_hi;

    uint32_t cy = (uint32_t)(((p0 >> 32) + (uint32_t)p1 + (uint32_t)p2) >> 32);

    *lo = p0 + (p1 << 32) + (p2 << 32);
    *hi = p3 + (p1 >> 32) + (p2 >> 32) + cy;
#endif // X86_INLINE_ASM
}

/* George Marsaglia's KISS64 generator, posted to comp.lang.c on 28 Feb 2009
   https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
*/
struct Prng_T {
    uint64_t x, c, y, z, t;
};

struct Prng_T kiss64 = {1234567890987654321ULL, 123456123456123456ULL,
                        362436362436362436ULL, 1066149217761810ULL, 0ULL};

/* KISS64 state equations */
#define MWC64 (kiss64->t = (kiss64->x << 58) + kiss64->c,            \
               kiss64->c = (kiss64->x >> 6), kiss64->x += kiss64->t, \
               kiss64->c += (kiss64->x < kiss64->t), kiss64->x)
#define XSH64 (kiss64->y ^= (kiss64->y << 13), kiss64->y ^= (kiss64->y >> 17), \
               kiss64->y ^= (kiss64->y << 43))
#define CNG64 (kiss64->z = 6906969069ULL * kiss64->z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)
uint64_t random64 (Prng_T kiss64)
{
    return KISS64;
}

int main (void)
{
    int i;
    Prng_T state = &kiss64;

    for (i = 0; i < 1000; i++) {
        printf ("%llu\n", randint (state, 10));
    }
    return EXIT_SUCCESS;
}

【讨论】：

【解决方案7】：

您谈到了涉及随机整数算法的两点：它是最优，它是无偏吗？

最佳

有很多方法可以定义“最佳”算法。在这里，我们根据平均使用的随机位数来查看“最佳”算法。从这个意义上说，rand 对于随机生成的数字来说是一种糟糕的方法，部分原因是它不一定需要产生随机位（因为没有完全指定 RAND_MAX）*。相反，我们假设我们有一个“真正的”随机生成器，它可以产生无偏且独立的随机位。

1976 年，DE Knuth 和 AC Yao 表明，任何仅使用随机位以给定概率生成随机整数的算法都可以表示为二叉树，其中随机位指示遍历树和每个叶子的方式（端点）对应于一个结果。（Knuth 和 Yao，“非均匀随机数生成的复杂性”，Algorithms and Complexity，1976 年。）他们还给出了给定算法平均需要的比特数来完成这项任务。在这种情况下，在[0, n) 中统一生成整数的最佳算法将需要平均至少需要log2(n) 和最多log2(n) + 2 位。

在这个意义上有很多最优算法的例子。其中一个是 J. Lumbroso (2013) 的 Fast Dice Roller（在下面实现），也许另一个例子是 2004 年 Math Forum 中给出的算法。另一方面，所有算法 surveyed by M. O'Neill 都不是最优的，因为它们依赖于一次生成随机位块。另请参阅我在 integer generating algorithms 上的说明。

下面显示了快速骰子滚子的实现；尽管它是在 JavaScript 中，而不是在 C 或 C++ 中，但它很容易适应任何一种语言，并且其想法是表明以最佳方式从位生成整数并不复杂。在代码中，(Math.random() < 0.5 ? 0 : 1) 是 JavaScript 生成无偏随机位的方式。

function randomInt(minInclusive, maxExclusive) {
  var maxInclusive = (maxExclusive - minInclusive) - 1
  var x = 1
  var y = 0
  while(true) {
    x = x * 2
    var randomBit = (Math.random() < 0.5 ? 0 : 1)
    y = y * 2 + randomBit
    if(x > maxInclusive) {
      if (y <= maxInclusive) { return y + minInclusive }
      // Rejection
      x = x - maxInclusive - 1
      y = y - maxInclusive - 1
    }
  }
}

无偏见

然而，任何最优整数生成器也是无偏的，一般来说，在最坏的情况下会永远运行，正如 Knuth 和 Yao 所展示的那样。回到二叉树，n 结果标签中的每一个都留在二叉树中，因此 [0, n) 中的每个整数都可以以 1/n 的概率出现。但是如果 1/n 有一个非终止的二元展开式（如果 n 不是 2 的幂就是这种情况），这棵二叉树必然要么——

具有“无限”深度，或
在树的末端包含“拒绝”叶子，

在任何一种情况下，算法都不会在恒定时间内运行，并且在最坏的情况下会永远运行。（另一方面，当n 是 2 的幂时，最优二叉树将具有有限深度且没有拒绝节点。）快速骰子滚轮是使用“拒绝”事件来确保其不偏不倚；请参阅上面代码中的注释。

对于一般n，没有办法在不引入偏差的情况下“修复”这种最坏情况的时间复杂度。例如，模减少（包括您问题中的min + (rand() % (int)(max - min + 1))）相当于一棵二叉树，其中拒绝叶子被标记的结果替换 - 但由于可能的结果比拒绝叶子多，因此只有一些结果可以采用拒绝离开的地方，引入了偏见。如果您在一定次数的迭代后停止拒绝，则会产生相同类型的二叉树 - 以及相同类型的偏差。（但是，根据应用程序，这种偏差可能可以忽略不计。随机整数生成也有安全方面的问题，这太复杂了，无法在此答案中讨论。）

注意

* 还有other problems with rand()。也许这里最严重的是，C 标准没有为rand() 返回的数字指定特定的分布。

【讨论】：