我看到很多人问过这个问题,但从未见过一个真正具体的答案。所以我打算在这里发布一个,希望能帮助人们理解为什么在使用随机数生成器时,比如c++中的rand(),会有“模偏置”。


因此rand()是一个伪随机数生成器,它在0和RAND_MAX之间选择一个自然数,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()执行RAND_MAX/n次调用。

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

int x;

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

x %= n;

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


引用作品及进一步阅读:

CPlusPlus参考 永远Confuzzled


不断随机选取是去除偏差的好方法。

更新

如果我们在能被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次迭代。

对于模的使用,有两种常见的抱怨。

one is valid for all generators. It is easier to see in a limit case. If your generator has a RAND_MAX which is 2 (that isn't compliant with the C standard) and you want only 0 or 1 as value, using modulo will generate 0 twice as often (when the generator generates 0 and 2) as it will generate 1 (when the generator generates 1). Note that this is true as soon as you don't drop values, whatever the mapping you are using from the generator values to the wanted one, one will occurs twice as often as the other. some kind of generator have their less significant bits less random than the other, at least for some of their parameters, but sadly those parameter have other interesting characteristic (such has being able to have RAND_MAX one less than a power of 2). The problem is well known and for a long time library implementation probably avoid the problem (for instance the sample rand() implementation in the C standard use this kind of generator, but drop the 16 less significant bits), but some like to complain about that and you may have bad luck

使用类似于

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的低值。他使用一个非常小的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很小的时候才会很明显。或者更具体地说,当你modding的数字比RAND_MAX大的时候。

一个比循环更好的解决方案(循环效率非常低,甚至不应该被建议使用)是使用输出范围大得多的PRNG。梅森Twister算法的最大输出为4,294,967,295。这样做MersenneTwister::genrand_int32() % 10,将是均匀分布的,模偏效应将几乎消失。

@user1413793 is correct about the problem. I'm not going to discuss that further, except to make one point: yes, for small values of n and large values of RAND_MAX, the modulo bias can be very small. But using a bias-inducing pattern means that you must consider the bias every time you calculate a random number and choose different patterns for different cases. And if you make the wrong choice, the bugs it introduces are subtle and almost impossible to unit test. Compared to just using the proper tool (such as arc4random_uniform), that's extra work, not less work. Doing more work and getting a worse solution is terrible engineering, especially when doing it right every time is easy on most platforms.

不幸的是,解决方案的实现都是不正确的,或者效率低于应有的水平。(每个解决方案都有各种解释问题的评论,但没有一个解决方案被修复以解决这些问题。)这可能会让那些随意寻求答案的人感到困惑,所以我在这里提供了一个已知的良好实现。

同样,最好的解决方案是在提供arc4random_uniform的平台上使用它,或者为您的平台使用类似的远程解决方案(如Random。nextInt在Java)。它将在没有代码成本的情况下做正确的事情。这几乎总是正确的选择。

如果你没有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;
}

对于那些需要实现类似事情的人来说,值得注意这段代码上的最新commit注释:

更改arc4random_uniform()计算2** 32% upper_bound为 -upper_bound % upper_bound。简化代码并使之成为 在ILP32和LP64架构上都是一样的,而且速度也略快 LP64架构使用32位余数而不是64位余数 余数。 由Jorden Verwer在tech@上指出 好的deraadt;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;
 }

我刚刚为冯·诺依曼无偏抛硬币法写了一段代码,理论上应该可以消除随机数生成过程中的任何偏差。更多信息请访问(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;
        }
    }
}

定义

模偏置是使用模算术将输出集缩减为输入集的子集时的固有偏置。一般来说,只要输入和输出集之间的映射不是均匀分布的,就会存在偏置,例如当输出集的大小不是输入集大小的除数时使用模算术。

这种偏差在计算中尤其难以避免,在计算中,数字被表示为比特串: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和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值来查看在大多数情况下实际发生了多少次重新滚动。持怀疑态度的人也可能希望将计算值保存到文件中,并验证分布是否正常。

RAND_MAX值为3(实际上它应该比这个值高得多,但偏差仍然存在),从这些计算中可以看出存在偏差:

1% 2 = 1 2% 2 = 0 3% 2 = 1 Random_between(1,3) % 2 =更可能是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

马克的解决方案(公认的解决方案)近乎完美。

int x; {做 X = rand(); } while (x >= (RAND_MAX - RAND_MAX % n)); X %= n; 编辑于2016年3月25日23:16 Mark Amery 39k21170211

然而,它有一个警告,在RAND_MAX (RM)小于N的倍数(其中N =可能有效结果的数量)的任何场景中,它会丢弃1个有效结果集。

也就是说,当'count of values discarded' (D)等于N时,那么它们实际上是一个有效的集合(V),而不是一个无效的集合(I)。

造成这种情况的原因是Mark在某些时候忽略了N和Rand_Max之间的差异。

N是一个集合,它的有效成员仅由正整数组成,因为它包含了有效响应的计数。(例如:Set N ={1,2,3,…N})

Rand_max然而是一个集合,它(根据我们的目的定义)包括任意数量的非负整数。

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

因此,Rand_Max最好被定义为“可能的响应”的集合。

然而,N是针对有效响应集合中的值的计数进行操作的,因此即使在我们的特定情况下定义,Rand_Max也将是一个比它所包含的总数小1的值。

使用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不等于质数时,这是有效的关注)

EG:

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,所需的rerroll的百分比就越高,因此根据运行代码的系统的约束条件和所寻找的值,在许多不同的值上,这可能是值得关注的问题。

要否定这一点,我们可以做一个简单的修正,如下所示:

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

这提供了一个更通用的公式版本,说明了使用模量定义最大值的附加特性。

使用小值RAND_MAX的例子,它是N的乘法。

Mark’original版:

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调用时,返回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条件之外执行比较操作的计算可能会导致更快的结果,因为无论需要重试多少次,这都是一次性计算。YMMV !

// 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();
}

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

When the range of possible integers is unknown, however, there is no way in general to "fix" this worst case of running forever without introducing bias. It's not just modulo reduction (rand() % n, discussed in the accepted answer) that will introduce bias this way, but also the "multiply-and-shift" reduction of Daniel Lemire, or if you stop rejecting an outcome after a set number of iterations. (To be clear, this doesn't mean there is no way to fix the bias issues present in pseudorandom generators. For example, even though modulo and other reductions are biased in general, they will have no issues with bias if the range of possible integers is a power of 2 and if the random generator produces unbiased random bits or blocks of them.)

这个答案的其余部分将显示随机生成器中运行时间和偏差之间的关系。从这里开始,我们将假设我们有一个“真正的”随机生成器,可以产生无偏和独立的随机比特

In 1976, D. E. Knuth and A. C. Yao showed that any algorithm that produces random integers with a given probability, using only random bits, can be represented as a binary tree, where random bits indicate which way to traverse the tree and each leaf (endpoint) corresponds to an outcome. In this case, we're dealing with algorithms that generate random integers in [0, n), where each integer is chosen with probability 1/n. The algorithm is unbiased if the same number of leaves appear in the tree for all outcomes. But if 1/n has a non-terminating binary expansion (which will be the case if n is not a power of 2), the algorithm will be unbiased only if—

二叉树具有“无限”深度,或者 二叉树的末端包含“排斥”叶,

无论哪种情况,算法都不会在常数时间内运行在最坏的情况下会一直运行下去。(另一方面,当n是2的幂时,最优二叉树的深度是有限的,并且没有拒绝节点。)

The binary tree concept also shows that any way to "fix" this worst-case time complexity will lead to bias in general. (Again, this doesn't mean there is no way to fix the bias issues present in pseudorandom generators.) For instance, modulo reductions are equivalent to a binary tree in which rejection leaves are replaced with labeled outcomes — but since there are more possible outcomes than rejection leaves, only some of the outcomes can take the place of the rejection leaves, introducing bias. The same kind of binary tree — and the same kind of bias — results if you stop rejecting after a set number of iterations. (However, this bias may be negligible depending on the application. There are also security aspects to random integer generation, which are too complicated to discuss in this answer.)

为了说明这一点,下面的JavaScript代码实现了J. Lumbroso(2013)提出的名为Fast Dice Roller的随机整数算法。请注意,它包括一个拒绝事件和一个循环,这是在一般情况下使算法无偏倚所必需的。

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
    }
 }
}

Note

*这个答案不会涉及到C中的rand()函数,因为它有很多问题。这里最严重的问题可能是,C标准没有明确地为rand()返回的数字指定特定的分布,甚至没有统一的分布。