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


当前回答

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

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

其他回答

因此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


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

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

定义

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

这种偏差在计算中尤其难以避免,在计算中,数字被表示为比特串: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的低值。他使用一个非常小的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,将是均匀分布的,模偏效应将几乎消失。