为什么人们说使用随机数生成器时存在模偏差

Why do people say there is modulo bias when using a random number generator?

本文关键字：存在随机数生成器人们说为什么更新时间：2023-10-16

我见过这个问题问了很多，但从未见过真正的具体答案。因此，我将在这里发布一个，希望能帮助人们理解为什么在使用随机数生成器时存在"模偏差"，例如C++中的rand()。

所以rand()是一个伪随机数生成器，它在0和RAND_MAX之间选择一个自然数，这是一个在cstdlib中定义的常量（有关rand()的一般概述，请参阅本文）。

现在，如果您想生成一个介于 0 和 2 之间的随机数，会发生什么？为了解释起见，假设RAND_MAX是 10，我决定通过调用 rand()%3 生成一个介于 0 和 2 之间的随机数。但是，rand()%3不会以相等的概率产生 0 到 2 之间的数字！

当rand()返回 0、3、6 或 9 时，rand()%3 == 0 .因此，P（0） = 4/11

当rand()返回 1、4、7 或 10 时，rand()%3 == 1 .因此，P（1） = 4/11

当rand()返回 2、5 或 8 时，rand()%3 == 2 .因此，P（2） = 3/11

这不会以相等的概率生成 0 到 2 之间的数字。当然，对于小范围，这可能不是最大的问题，但对于较大的范围，这可能会扭曲分布，使较小的数字偏向。

那么rand()%n什么时候以相等的概率返回从 0 到 n-1 的数字范围呢？当RAND_MAX%n == n - 1.在这种情况下，除了我们之前的假设rand()确实以相等的概率返回一个介于 0 和 RAND_MAX 之间的数字外，n 的模类也将平均分布。

那么我们如何解决这个问题呢？一种粗略的方法是不断生成随机数，直到获得所需范围内的数字：

int x; 
do {
    x = rand();
} while (x >= n);

但是对于低值的n，这是低效的，因为您只有n/RAND_MAX的机会获得范围内的值，因此您需要执行RAND_MAX/n调用才能平均rand()。

一种更有效的公式方法是取一些长度可被n整除的大范围，如RAND_MAX - RAND_MAX % n，不断生成随机数，直到得到一个位于范围内的随机数，然后取模：

int x;
do {
    x = rand();
} while (x >= (RAND_MAX - RAND_MAX % n));
x %= n;

对于n的小值，这很少需要多次调用rand()。

引用作品及延伸阅读：

CPlusPlus参考
永远的迷惑

继续选择随机是消除偏见的好方法。

更新

如果我们在可被 n 整除的范围内搜索 x，我们可以使代码快速。

// Assumptions
// rand() in [0, RAND_MAX]
// n in (0, RAND_MAX]
int x; 
// Keep searching for an x in a range divisible by n 
do {
    x = rand();
} while (x >= RAND_MAX - (RAND_MAX % n)) 
x %= n;

上面的循环应该非常快，比如平均 1 次迭代。

@user1413793对

这个问题是正确的。我不打算进一步讨论这个问题，只是要指出一点：是的，对于n的小值和RAND_MAX的大值，模偏差可能非常小。但是使用偏差诱导模式意味着每次计算随机数并为不同情况选择不同的模式时都必须考虑偏差。如果你做出了错误的选择，它引入的错误是微妙的，几乎不可能进行单元测试。与仅使用适当的工具（例如arc4random_uniform）相比，这是额外的工作，而不是更少的工作。做更多的工作并获得更糟糕的解决方案是可怕的工程，尤其是在大多数平台上每次都做对的时候。

不幸的是，该解决方案的实现都是不正确的或效率低于应有的效率。（每个解决方案都有各种注释来解释问题，但尚未修复任何解决方案来解决这些问题。这可能会使随意的答案寻求者感到困惑，因此我在这里提供了一个已知良好的实现。

同样，最好的解决方案只是在提供它的平台上使用arc4random_uniform，或者为您的平台使用类似的范围解决方案（例如 Java 上的 Random.nextInt）。它将做正确的事情，而无需花费代码。这几乎总是正确的选择。

如果你没有arc4random_uniform，那么你可以利用开源的力量来确切地看到它是如何在更广泛的RNG之上实现的（在这种情况下ar4random，但类似的方法也可以在其他RNG之上工作）。

下面是 OpenBSD 的实现：

/*
 * Calculate a uniformly distributed random number less than upper_bound
 * avoiding "modulo bias".
 *
 * Uniformity is achieved by generating new random numbers until the one
 * returned is outside the range [0, 2**32 % upper_bound).  This
 * guarantees the selected random number will be inside
 * [2**32 % upper_bound, 2**32) which maps back to [0, upper_bound)
 * after reduction modulo upper_bound.
 */
u_int32_t
arc4random_uniform(u_int32_t upper_bound)
{
    u_int32_t r, min;
    if (upper_bound < 2)
        return 0;
    /* 2**32 % x == (2**32 - x) % x */
    min = -upper_bound % upper_bound;
    /*
     * This could theoretically loop forever but each retry has
     * p > 0.5 (worst case, usually far better) of selecting a
     * number inside the range we need, so it should rarely need
     * to re-roll.
     */
    for (;;) {
        r = arc4random();
        if (r >= min)
            break;
    }
    return r % upper_bound;
}

值得注意的是，对于那些需要实现类似事情的人来说，关于这段代码的最新提交注释：

更改 arc4random_uniform（）以将2**32 % upper_bound计算为 -upper_bound % upper_bound . 简化代码并使其成为在 ILP32 和 LP64 架构上相同，在使用 32 位余数而不是 64 位余数的 LP64 体系结构剩余。
Jorden Verwer在tech@上指出好的德拉特;DJM或OTTO没有异议

Java 实现也很容易找到（请参阅上一个链接）：

public int nextInt(int n) {
   if (n <= 0)
     throw new IllegalArgumentException("n must be positive");
   if ((n & -n) == n)  // i.e., n is a power of 2
     return (int)((n * (long)next(31)) >> 31);
   int bits, val;
   do {
       bits = next(31);
       val = bits % n;
   } while (bits - val + (n-1) < 0);
   return val;
 }

定义

模

偏置是使用模算术将输出集减少为输入集子集的固有偏差。通常，只要输入集和输出集之间的映射分布不均匀，就会存在偏差，例如当输出集的大小不是输入集大小的除数时，使用模算术的情况。

这种偏差在计算中特别难以避免，其中数字表示为位字符串：0 和 1。找到真正的随机性来源也非常困难，但超出了本次讨论的范围。对于此答案的其余部分，假设存在无限的真正随机位源。

问题示例

让我们考虑使用这些随机位模拟骰子（0 到 5）。有 6 种可能性，因此我们需要足够的位来表示数字 6，即 3 位。不幸的是，3 个随机位会产生 8 种可能的结果：

000 = 0, 001 = 1, 010 = 2, 011 = 3
100 = 4, 101 = 5, 110 = 6, 111 = 7

我们可以通过取模 6 来将结果集的大小减少到 6，但这会带来模偏差问题：110产生 0，111产生 1。此模具已加载。

潜在解决方案

方法 0：

理论上，人们可以雇佣一小支军队整天掷骰子并将结果记录在数据库中，然后每个结果只使用一次，而不是依赖随机比特。这听起来很实用，而且很可能不会产生真正的随机结果（双关语）。

方法1：

一个幼稚但数学上正确的解决方案不是使用模数，而是丢弃产生110和111的结果，只需使用 3 个新位重试。不幸的是，这意味着每次滚动都有 25% 的机会需要重新滚动，包括每个重新滚动本身。除了最微不足道的用途之外，这显然是不切实际的。

方法2：

使用更多位：使用 4 位而不是 3 位。这会产生 16 种可能的结果。当然，只要结果大于 5，重新滚动就会让事情变得更糟（10/16 = 62.5%），因此仅凭这一点无济于事。

请注意，2 * 6 = 12 <16，因此我们可以安全地获取任何小于 12 的结果，并减少该模 6 以均匀分布结果。其他 4 个结果必须丢弃，然后像前面的方法一样重新滚动。

起初听起来不错，但让我们检查一下数学：

4 discarded results / 16 possibilities = 25%

在这种情况下，额外的 1 位根本没有帮助！

这个结果很不幸，但让我们再试一次 5 位：

32 % 6 = 2 discarded results; and
2 discarded results / 32 possibilities = 6.25%

这是一个明确的改进，但在许多实际案例中还不够好。好消息是，添加更多位永远不会增加需要丢弃和重新滚动的机会。这不仅适用于骰子，而且适用于所有情况。

但是，如前所述，添加 1 个额外的位可能不会改变任何东西。事实上，如果我们将滚动增加到 6 位，概率仍然是 6.25%。

这引出了另外两个问题：

如果我们添加足够的位，是否可以保证丢弃的概率会降低？
在一般情况下，多少位才足够？

通用解决方案

值得庆幸的是，第一个问题的答案是肯定的。6 的问题在于 2^x mod 6 在 2 和 4 之间翻转，巧合的是 2 的倍数，因此对于偶数 x> 1，

[2^x mod 6] / 2^x == [2^(x+1) mod 6] / 2^(x+1)

因此，6 是一个例外而不是规则。有可能找到更大的模量，以相同的方式产生 2 的连续幂，但最终这必须环绕，并且丢弃的概率将降低。

没有提供进一步的证据，通常使用双倍的数字所需的位将提供一个较小的，通常微不足道的，丢弃的机会。

概念验证

这是一个使用OpenSSL的libcrypo提供随机字节的示例程序。编译时，请务必链接到大多数人都应该可用的-lcrypto库。

#include <iostream>
#include <assert.h>
#include <limits>
#include <openssl/rand.h>
volatile uint32_t dummy;
uint64_t discardCount;
uint32_t uniformRandomUint32(uint32_t upperBound)
{
    assert(RAND_status() == 1);
    uint64_t discard = (std::numeric_limits<uint64_t>::max() - upperBound) % upperBound;
    RAND_bytes((uint8_t*)(&randomPool), sizeof(randomPool));
    while(randomPool > (std::numeric_limits<uint64_t>::max() - discard)) {
        RAND_bytes((uint8_t*)(&randomPool), sizeof(randomPool));
        ++discardCount;
    }
    return randomPool % upperBound;
}
int main() {
    discardCount = 0;
    const uint32_t MODULUS = (1ul << 31)-1;
    const uint32_t ROLLS = 10000000;
    for(uint32_t i = 0; i < ROLLS; ++i) {
        dummy = uniformRandomUint32(MODULUS);
    }
    std::cout << "Discard count = " << discardCount << std::endl;
}

我鼓励使用MODULUS和ROLLS值，看看在大多数情况下实际发生了多少次重滚。持怀疑态度的人可能还希望将计算值保存到文件中，并验证分布是否正常。

马克的解决方案（公认的解决方案）是近乎完美的。

int x;
do {
    x = rand();
} while (x >= (RAND_MAX - RAND_MAX % n));
x %= n;
编辑于 Mar 25 '16 at 23：16
马克·阿默里 39k21170211

但是，它有一个警告，即在RAND_MAX（RM）小于N的倍数（其中N=可能的有效结果数）的任何情况下，丢弃1个有效结果集。

即，当"丢弃的值计数"（D）等于N时，那么它们实际上是一个有效的集合（V)，而不是一个无效的集合（I）。

导致这种情况的原因是在某些时候马克忽视了N和Rand_Max之间的区别。

N 是一个有效成员仅由正整数组成的集合，因为它包含有效的响应计数。（例如：设置N = {1, 2, 3, ... n }）

Rand_max 然而，它是一个集合（为我们的目的而定义）包括任意数量的非负整数。

在最通用的形式中，这里定义为Rand Max是所有有效结果的集合，理论上可以包括负数或非数值。

因此，Rand_Max最好定义为"可能的响应"集合。

但是，N对有效响应集中的值计数进行操作，因此即使在我们的特定情况下定义，Rand_Max的值也将比它包含的总数少一个。

使用 Mark 的解决方案，在以下情况下丢弃值： X => RM - RM % N

EG: 
Ran Max Value (RM) = 255
Valid Outcome (N) = 4
When X => 252, Discarded values for X are: 252, 253, 254, 255
So, if Random Value Selected (X) = {252, 253, 254, 255}
Number of discarded Values (I) = RM % N + 1 == N
 IE:
 I = RM % N + 1
 I = 255 % 4 + 1
 I = 3 + 1
 I = 4
   X => ( RM - RM % N )
 255 => (255 - 255 % 4) 
 255 => (255 - 3)
 255 => (252)
 Discard Returns $True

正如你在上面的例子中看到的，当 X（我们从初始函数获得的随机数）的值为 252、253、254 或 255 时，即使这四个值包含一组有效的返回值，我们也会丢弃它。

IE：当丢弃的值计数（I） = N（有效结果的数量）时，原始函数将丢弃一组有效的返回值。

如果我们将值 N 和 RM 之间的差异描述为 D，即：

D = (RM - N)

然后，随着 D 的值变小，由于此方法导致的不需要的重新滚动的百分比在每个自然乘法时都会增加。（当RAND_MAX不等于质数时，这是有效的问题）

例如：

RM=255 , N=2 Then: D = 253, Lost percentage = 0.78125%
RM=255 , N=4 Then: D = 251, Lost percentage = 1.5625%
RM=255 , N=8 Then: D = 247, Lost percentage = 3.125%
RM=255 , N=16 Then: D = 239, Lost percentage = 6.25%
RM=255 , N=32 Then: D = 223, Lost percentage = 12.5%
RM=255 , N=64 Then: D = 191, Lost percentage = 25%
RM=255 , N= 128 Then D = 127, Lost percentage = 50%

由于所需的重滚百分比随着 N 越接近 RM，因此在许多不同的值下，这可能是有效的关注点，具体取决于运行代码的系统约束和要查找的值。

为了否定这一点，我们可以做一个简单的修正，如下所示：

 int x;
 
 do {
     x = rand();
 } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) );
 
 x %= n;

这提供了公式的更通用版本，该公式考虑了使用模量定义最大值的额外特性。

对 N 的乘法RAND_MAX使用小值的示例。

标记'原始版本：

RAND_MAX = 3, n = 2, Values in RAND_MAX = 0,1,2,3, Valid Sets = 0,1 and 2,3.
When X >= (RAND_MAX - ( RAND_MAX % n ) )
When X >= 2 the value will be discarded, even though the set is valid.

通用版本 1：

RAND_MAX = 3, n = 2, Values in RAND_MAX = 0,1,2,3, Valid Sets = 0,1 and 2,3.
When X > (RAND_MAX - ( ( RAND_MAX % n  ) + 1 ) % n )
When X > 3 the value would be discarded, but this is not a vlue in the set RAND_MAX so there will be no discard.

此外，在 N 应该是 RAND_MAX 中的值数的情况下;在这种情况下，您可以设置 N = RAND_MAX +1，除非 RAND_MAX = INT_MAX。

循环方面，你可以只使用 N = 1，但是，X 的任何值都将被接受，并为最终乘数输入一个 IF 语句。但也许您的代码可能有正当理由在调用函数时返回 n = 1...

因此，当您

希望 n = RAND_MAX+1 时，最好使用 0，它通常会提供 Div 0 错误

通用版本 2：

int x;
if n != 0 {
    do {
        x = rand();
    } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) );
    x %= n;
} else {
    x = rand();
}

这两种解决方案都解决了当 RM+1 是 n 的乘积时会发生的不必要丢弃有效结果的问题。

第二个版本还涵盖了边缘情况，当您需要 n 等于 RAND_MAX 中包含的总可能值集时。

两者的修改方法是相同的，并且允许提供更通用的解决方案，以满足提供有效随机数和最小化丢弃值的需求。

重申：

扩展 mark 示例的基本通用解决方案：

// Assumes:
//  RAND_MAX is a globally defined constant, returned from the environment.
//  int n; // User input, or externally defined, number of valid choices.
 int x;
 
 do {
     x = rand();
 } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) ) );
 
 x %= n;

扩展通用解决方案，允许一个 RAND_MAX+1 = n 的附加方案：

// Assumes:
//  RAND_MAX is a globally defined constant, returned from the environment.
//  int n; // User input, or externally defined, number of valid choices.
int x;
if n != 0 {
    do {
        x = rand();
    } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) ) );
    x %= n;
} else {
    x = rand();
}

在某些语言（特别是解释型语言）中，在while条件之外进行比较操作的计算可能会导致更快的结果，因为无论需要多少次重试，这都是一次性计算。扬子晚报！

// Assumes:
//  RAND_MAX is a globally defined constant, returned from the environment.
//  int n; // User input, or externally defined, number of valid choices.
int x; // Resulting random number
int y; // One-time calculation of the compare value for x
y = RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) 
if n != 0 {
    do {
        x = rand();
    } while (x > y);
    x %= n;
} else {
    x = rand();
}

使用模有两个常见的抱怨。

一个对所有生成器都有效。在极限情况下更容易看到。如果您的生成器的RAND_MAX为 2（不符合 C 标准），并且您只需要 0 或 1 作为值，则使用模数生成 0 的频率（当生成器生成 0 和 2 时）是生成 1（当生成器生成 1 时）的两倍。请注意，只要您不删除值，这是正确的，无论您使用从生成器值到所需值的映射，一个值的发生频率将是另一个的两倍。
某种生成器的不那么有效比另一个的随机位少，但遗憾的是，这些参数具有其他有趣的特征（例如能够RAND_MAX小于 2 的幂）。这个问题是众所周知的，并且在很长一段时间内，库实现可能会避免这个问题（例如，C 标准中的示例 rand（）实现使用这种生成器，但删除了 16 个不太有效的位），但有些人喜欢抱怨这一点，你可能会运气不好

使用类似的东西

int alea(int n){ 
 assert (0 < n && n <= RAND_MAX); 
 int partSize = 
      n == RAND_MAX ? 1 : 1 + (RAND_MAX-n)/(n+1); 
 int maxUsefull = partSize * n + (partSize-1); 
 int draw; 
 do { 
   draw = rand(); 
 } while (draw > maxUsefull); 
 return draw/partSize; 
}

生成一个介于 0 和 n 之间的随机数将避免这两个问题（并且避免溢出 RAND_MAX == INT_MAX）

顺便说一句，C++11 引入了除 rand（）之外的归约和其他生成器的标准方法。

当RAND_MAX值为 3 时（实际上它应该比这高得多，但偏差仍然存在），从这些计算中可以理解存在偏差：

1 % 2 = 1 2 % 2 = 0 3 % 2 = 1 random_between(1, 3) % 2 = more likely a 1

在这种情况下，当您想要 0 和 1 之间的随机数时，% 2是你不应该做的。不过，您可以通过执行% 3来获得 0 到 2 之间的随机数，因为在这种情况下：RAND_MAX 是 3 的倍数。

另一种方法

还有更简单的，但要添加到其他答案中，这是我的解决方案，可以在 0 和 n - 1 之间获取一个随机数，因此n不同的可能性，没有偏见。

编码可能性数所需的位数（不是字节数）是您需要的随机数据位数
对随机位中的数字进行编码
如果此数字为 >= n ，则重新启动（无模数）。

真正的随机数据并不容易获得，所以为什么要使用比需要的更多的位。

下面是 Smalltalk 中的一个示例，使用来自伪随机数生成器的位缓存。我不是安全专家，所以使用风险自负。

next: n
    | bitSize r from to |
    n < 0 ifTrue: [^0 - (self next: 0 - n)].
    n = 0 ifTrue: [^nil].
    n = 1 ifTrue: [^0].
    cache isNil ifTrue: [cache := OrderedCollection new].
    cache size < (self randmax highBit) ifTrue: [
        Security.DSSRandom default next asByteArray do: [ :byte |
            (1 to: 8) do: [ :i |    cache add: (byte bitAt: i)]
        ]
    ].
    r := 0.
    bitSize := n highBit.
    to := cache size.
    from := to - bitSize + 1.
    (from to: to) do: [ :i |
        r := r bitAt: i - from + 1 put: (cache at: i)
    ].
    cache removeFrom: from to: to.
    r >= n ifTrue: [^self next: n].
    ^r

模约简是一种常见的方法，可以使随机整数生成器避免永远运行的最坏情况。

然而，当可能的整数范围未知时，通常没有办法在不引入偏差的情况下"修复"这种永远运行的最坏情况。不仅仅是模约简（rand() % n，在公认的答案中讨论）会以这种方式引入偏差，还有Daniel Lemire的"乘移"约简，或者如果你在一定次数的迭代后停止拒绝一个结果。（需要明确的是，这并不意味着没有办法解决伪随机生成器中存在的偏差问题。例如，即使模和其他约简通常是有偏差的，但如果可能的整数范围是 2 的幂，并且随机生成器产生无偏的随机位或块，它们也不会有偏差问题。

我的以下答案讨论了随机生成器中的运行时间和偏差之间的关系，假设我们有一个"真正的"随机生成器，可以产生无偏差和独立的随机位。答案甚至不涉及 C 中的 rand() 函数，因为它有很多问题。也许这里最严重的是C标准没有明确地为rand()返回的数字指定一个特定的分布，甚至没有一个均匀的分布。

如何在不浪费位的情况下从随机位流中生成 [0，n] 范围内的随机整数？

正如公认的答案所表明的那样，"模偏置"的根源在于RAND_MAX的低值。他使用极小的值 RAND_MAX （10）来表明，如果RAND_MAX是 10，那么您尝试使用 % 生成一个介于 0 和 2 之间的数字，结果如下：

rand() % 3   // if RAND_MAX were only 10, gives
output of rand()   |   rand()%3
0                  |   0
1                  |   1
2                  |   2
3                  |   0
4                  |   1
5                  |   2
6                  |   0
7                  |   1
8                  |   2
9                  |   0

因此，有 4 个 0 输出（4/10 机会

），只有 3 个输出 1 和 2（每个 3/10 机会）。

所以这是有偏见的。数字越低，出局的机会越大。

但这只有在RAND_MAX很小的时候才会如此明显地表现出来。或者更具体地说，当您修改的数字与RAND_MAX相比很大时。

比循环（效率极低，甚至不应该被建议）更好的解决方案是使用具有更大输出范围的 PRNG。 Mersenne Twister 算法的最大输出为 4，294，967，295。因此，出于所有意图和目的，做MersenneTwister::genrand_int32() % 10将平均分布，模偏置效应将几乎消失。

我想要各种软件的随机双倍。如果我使用（（double）rand（）/RAND_MAX，我发现范围更"随机"。所以我猜如果你乘以你的数字范围，你可以得到一个偏差较小的随机数？

即（（

双）兰德（）/RAND_MAX） * 3.

我读到了一个关于从 2 中做一个随机数的答案。 isodd（rand（））？

我刚刚为冯·诺依曼的无偏硬币翻转方法编写了一段代码，理论上应该消除随机数生成过程中的任何偏差。更多信息可在（http://en.wikipedia.org/wiki/Fair_coin）中找到

int unbiased_random_bit() {    
    int x1, x2, prev;
    prev = 2;
    x1 = rand() % 2;
    x2 = rand() % 2;
    for (;; x1 = rand() % 2, x2 = rand() % 2)
    {
        if (x1 ^ x2)      // 01 -> 1, or 10 -> 0.
        {
            return x2;        
        }
        else if (x1 & x2)
        {
            if (!prev)    // 0011
                return 1;
            else
                prev = 1; // 1111 -> continue, bias unresolved
        }
        else
        {
            if (prev == 1)// 1100
                return 0;
            else          // 0000 -> continue, bias unresolved
                prev = 0;
        }
    }
}