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


当前回答

定义

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

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

其他回答

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

更新

如果我们在能被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()以外的其他生成器。

定义

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

这种偏差在计算中尤其难以避免,在计算中,数字被表示为比特串: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

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