为什么 rand()%6 有偏差?

Why is rand()%6 biased?

在阅读如何使用 std::rand 时,我在 cppreference.com

上找到了这段代码
int x = 7;
while(x > 6) 
    x = 1 + std::rand()/((RAND_MAX + 1u)/6);  // Note: 1+rand()%6 is biased

右边的表达有什么问题?试过了,效果很好。

这里有深藏不露:

  1. RAND_MAX + 1u中小u的使用。 RAND_MAX 被定义为 int 类型,并且通常是最大可能的 int。在 signed 类型溢出的情况下,RAND_MAX + 1 的行为将是 undefined。写入 1u 会强制将 RAND_MAX 类型转换为 unsigned,从而避免溢出。

  2. % 6的使用(但在std::rand的每个实现中我都看到没有't) 在所提供的替代方案之外引入任何额外的统计偏差。 % 6 是危险的这种情况是数字生成器在低阶位中具有相关平原的情况,例如我认为 1970 年代相当著名的 rand 的 IBM 实现(C 语言)将高位和低位翻转为 "a final flourish"。进一步的考虑是 6 非常小 cf。 RAND_MAX,因此如果 RAND_MAX 不是 6 的倍数(可能不是),影响将很小。

总而言之,最近,由于它的易处理性,我会使用 % 6。除了生成器本身引入的统计异常之外,它不太可能引入任何统计异常。如果您仍然有疑问,请测试您的生成器以查看它是否具有适合您的用例的统计属性。

rand() % 6 有两个问题(1+ 不影响任何一个问题)。

首先,正如几个答案所指出的,如果rand()的低位不适当统一,则余数运算符的结果也不统一。

其次,如果 rand() 产生的不同值的数量不是 6 的倍数,那么余数将产生比高值更多的低值。即使 rand() returns 完美分布值也是如此。

作为一个极端的例子,假设 rand() 产生 [0..6] 范围内均匀分布的值。如果您查看这些值的余数,当 rand() returns 范围内的值 [0..5] 时,余数会在 [0..5] 范围内产生均匀分布的结果。当 rand() returns 6, rand() % 6 returns 0 时,就好像 rand() 返回了 0。所以你得到的 0 的分布是其他分布的两倍值。

第二个是真正的问题rand() % 6

避免该问题的方法是丢弃 会产生不均匀重复的值。您计算小于或等于 RAND_MAX 的 6 的最大倍数,并且每当 rand() returns 大于或等于该倍数的值时,您将拒绝它并调用 `rand()再次,需要多少次。

所以:

int max = 6 * ((RAND_MAX + 1u) / 6)
int value = rand();
while (value >= max)
    value = rand();

这是相关代码的不同实现方式,旨在更清楚地展示正在发生的事情。

此示例代码说明 std::rand 是一个遗留的货物崇拜胡言乱语的案例,每次看到它都会让您大吃一惊。

这里有几个问题:

人们通常假设的合同——即使是可怜的不幸的灵魂,他们不知道更好,也不会准确地用这些术语来思考——是 rand 来自 的样本在 0、1、2、...、RAND_MAX 中的整数上均匀分布 ,并且每次调用都会产生一个 独立 样本。

第一个问题是假设的合同,每次调用中的独立统一随机样本,实际上并不是文档所说的那样——在实践中,历史上的实现甚至无法提供最独立的模拟. 例如,C99 §7.20.2.1 'rand 函数'说,没有详细说明:

The rand function computes a sequence of pseudo-random integers in the range 0 to RAND_MAX.

这是一个无意义的句子,因为伪随机性是函数(或函数族)的属性,不是一个整数,但这并不能阻止甚至 ISO 官僚滥用该语言。毕竟,只有 reader 的人会因为害怕他们的脑细胞衰退而不去阅读 rand 的文档。

C 中典型的历史实现是这样的:

static unsigned int seed = 1;

static void
srand(unsigned int s)
{
    seed = s;
}

static unsigned int
rand(void)
{
    seed = (seed*1103515245 + 12345) % ((unsigned long)RAND_MAX + 1);
    return (int)seed;
}

这有不幸的属性,即使单个样本可能在均匀随机种子下均匀分布(这取决于[=15的具体值=]), 它在连续调用中在偶数和奇数整数之间交替——在

之后
int a = rand();
int b = rand();

表达式 (a & 1) ^ (b & 1) 以 100% 的概率产生 1,对于 independent 支持偶数和奇数整数的任何分布的随机样本,情况并非如此。因此,出现了一种货物崇拜,人们应该丢弃低位比特以追逐难以捉摸的“更好的随机性”野兽。 (剧透警告:这不是一个技术术语。这表明您正在阅读谁的散文要么不知道他们在说什么,要么认为 毫无头绪,必须屈尊。)

第二个问题是 即使每个调用都独立于均匀随机分布 在 0, 1, 2, …, RAND_MAXrand() % 6 的结果不会像掷骰子一样均匀分布在 0、1、2、3、4、5 中,除非 RAND_MAX 与 -1 模 6 一致。 简单的反例:如果RAND_MAX = 6,那么从rand(),所有结果的概率都是1/7,但是从rand() % 6,结果0的概率是2/7,而所有其他结果有概率 1/7.

正确的做法是拒绝抽样:重复从0抽取一个独立的均匀随机样本s , 1, 2, …, RAND_MAX, 和 拒绝 (例如)结果 0, 1, 2, …, ((RAND_MAX + 1) % 6) - 1——如果你得到其中之一那些,重新开始;否则,yield s % 6.

unsigned int s;
while ((s = rand()) < ((unsigned long)RAND_MAX + 1) % 6)
    continue;
return s % 6;

这样,我们接受的 rand() 的结果集可以被 6 整除,并且 s % 6 的每个可能结果都由相同数量的 接受 来自 rand() 的结果,因此如果 rand() 是均匀分布的,那么 s 也是均匀分布的。试验次数没有bound,但期望次数小于2,成功概率随次数呈指数增长试验。

您拒绝 rand() 哪个 结果的选择是 immaterial,前提是您将相同数量的结果映射到 6 以下的每个整数. cppreference.com 处的代码做出了 不同的 选择,因为上面的第一个问题——rand() 的输出的分布或独立性没有任何保证,并且在实践中,低阶位表现出的模式“看起来不够随机”(不要介意下一个输出是前一个输出的确定性函数)。

reader 的练习:如果 rand() 在 0, 1, 2, …, RAND_MAX.

reader 的练习:为什么您可能更愿意拒绝一个或其他子集?两种情况下每次试验需要什么计算?

第三个问题是种子 space 太小了,即使种子是均匀分布的,一个知道你的程序和一个结果但不知道种子的对手也可以很容易预测种子和后续结果,这使得它们毕竟看起来不是那么随机。 所以甚至不要考虑将其用于密码学。

您可以走花哨的过度设计路线和 C++11 的 std::uniform_int_distribution class,使用适当的随机设备和您最喜欢的随机引擎,例如广受欢迎的梅森扭曲器 std::mt19937和你四岁的堂兄一起玩骰子,但即使这样也不适合生成加密密钥 material——而且梅森扭曲器也是一个可怕的 space 猪,有几千字节state wreaking havoc on your CPU's cache with an obscene setup time,所以即使对于 e.g.,并行 Monte Carlo 具有可重现的子计算树的模拟也是不好的;它的受欢迎程度可能主要来自其朗朗上口的名字。但是你可以像这个例子一样用它来玩玩具掷骰子!

另一种方法是使用具有小状态的简单密码伪随机数生成器,例如简单的 fast key erasure PRNG,或者如果您有信心(例如,在自然科学研究的Monte Carlo模拟中)如果状态受到损害,预测过去的结果不会产生不利后果。

我无论如何都不是经验丰富的 C++ 用户,但有兴趣看看是否有关于 std::rand()/((RAND_MAX + 1u)/6)1+std::rand()%6 更少的偏见实际上是正确的。所以我写了一个测试程序来列出两种方法的结果(我已经很久没有写过 C++,请检查一下)。 link 运行 找到代码 here。同时转载如下:

// Example program
#include <cstdlib>
#include <iostream>
#include <ctime>
#include <string>

int main()
{
    std::srand(std::time(nullptr)); // use current time as seed for random generator

    // Roll the die 6000000 times using the supposedly unbiased method and keep track of the results

    int results[6] = {0,0,0,0,0,0};

    // roll a 6-sided die 20 times
    for (int n=0; n != 6000000; ++n) {
        int x = 7;
        while(x > 6) 
            x = 1 + std::rand()/((RAND_MAX + 1u)/6);  // Note: 1+rand()%6 is biased

        results[x-1]++;
    }

    for (int n=0; n !=6; n++) {
        std::cout << results[n] << ' ';
    }

    std::cout << "\n";


    // Roll the die 6000000 times using the supposedly biased method and keep track of the results

    int results_bias[6] = {0,0,0,0,0,0};

    // roll a 6-sided die 20 times
    for (int n=0; n != 6000000; ++n) {
        int x = 7;
        while(x > 6) 
            x = 1 + std::rand()%6;

        results_bias[x-1]++;
    }

    for (int n=0; n !=6; n++) {
        std::cout << results_bias[n] << ' ';
    }
}

然后我得到这个输出并使用 R 中的 chisq.test 函数进行 运行 卡方检验,看看结果是否与预期有显着差异。这个 stackexchange 问题更详细地介绍了使用卡方检验来测试模具公平性:How can I test whether a die is fair?。以下是几 运行 秒的结果:

> ?chisq.test
> unbias <- c(100150, 99658, 100319, 99342, 100418, 100113)
> bias <- c(100049, 100040, 100091, 99966, 100188, 99666 )

> chisq.test(unbias)

Chi-squared test for given probabilities

data:  unbias
X-squared = 8.6168, df = 5, p-value = 0.1254

> chisq.test(bias)

Chi-squared test for given probabilities

data:  bias
X-squared = 1.6034, df = 5, p-value = 0.9008

> unbias <- c(998630, 1001188, 998932, 1001048, 1000968, 999234 )
> bias <- c(1000071, 1000910, 999078, 1000080, 998786, 1001075   )
> chisq.test(unbias)

Chi-squared test for given probabilities

data:  unbias
X-squared = 7.051, df = 5, p-value = 0.2169

> chisq.test(bias)

Chi-squared test for given probabilities

data:  bias
X-squared = 4.319, df = 5, p-value = 0.5045

> unbias <- c(998630, 999010, 1000736, 999142, 1000631, 1001851)
> bias <- c(999803, 998651, 1000639, 1000735, 1000064,1000108)
> chisq.test(unbias)

Chi-squared test for given probabilities

data:  unbias
X-squared = 7.9592, df = 5, p-value = 0.1585

> chisq.test(bias)

Chi-squared test for given probabilities

data:  bias
X-squared = 2.8229, df = 5, p-value = 0.7273

在我所做的三个 运行 中,两种方法的 p 值始终大于用于检验显着性的典型 alpha 值 (0.05)。这意味着我们不会认为它们中的任何一个有偏见。有趣的是,所谓的无偏方法始终具有较低的 p 值,这表明它实际上可能更偏。需要注意的是我只做了 3 运行s.

更新:在我写答案时,Konrad Rudolph 发布了一个采用相同方法的答案,但得到的结果却截然不同。我没有资格评论他的回答,所以我将在这里解决。首先,最主要的是他使用的代码在每次 运行 时都为随机数生成器使用相同的种子。如果你改变种子,你实际上会得到各种各样的结果。第二,如果你不改变种子,而是改变试验次数,你也会得到各种各样的结果。尝试增加或减少一个数量级以了解我的意思。第三,在预期值不太准确的地方有一些整数 t运行cation 或四舍五入。它可能不足以产生影响,但它就在那里。

基本上,总而言之,他只是碰巧得到了正确的种子和试验次数,他可能会得到错误的结果。

人们可以将随机数生成器想象成处理二进制数字流。生成器通过将流切成块将其转换为数字。如果 std:rand 函数使用的 RAND_MAX 为 32767,则它在每个切片中使用 15 位。

当取 0 到 32767(含)之间的数字的模块时,会发现 5462 个“0”和“1”,但只有 5461 个“2”、“3”、“4”和“5”。因此结果是有偏差的。 RAND_MAX值越大,偏差越小,但不可避免。

不偏的是[0..(2^n)-1]范围内的数字。您可以通过提取 3 位,将它们转换为 0..7 范围内的整数并拒绝 6 和 7 来生成 0..5 范围内的(理论上)更好的数字。

人们希望比特流中的每一位都有相同的机会成为“0”或“1”,而不管它在流中的位置或其他比特的值。这在实践中异常困难。许多不同的软件 PRNG 实现在速度和质量之间提供了不同的折衷。诸如 std::rand 之类的线性同余生成器以最低的质量提供最快的速度。密码生成器以最低速度提供最高质量。