为啥 rand()%6 有偏见？

Posted 2023-02-25

【中文标题】为啥 rand()%6 有偏见？

【英文标题】：Why is rand()%6 biased?为什么 rand()%6 有偏见？

【发布时间】：2018-09-27 11:16:52

【问题描述】：

在阅读如何使用 std::rand 时，我在 cppreference.com 上找到了这段代码

int x = 7;
while(x > 6) 
    x = 1 + std::rand()/((RAND_MAX + 1u)/6);  // Note: 1+rand()%6 is biased

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

【问题讨论】：

请注意，将std::uniform_int_distribution 用于骰子更好

@Caleth 是的，这只是为了理解为什么这段代码是“错误的”..

将“错误”改为“有偏见”

rand() 在典型实现中非常糟糕，您不妨使用xkcd RNG。所以这是错误的，因为它使用了rand()。

我写了这个东西（好吧，不是评论——那是@Cubbi），我当时的想法是Pete Becker's answer 解释的。 （仅供参考，这与 libstdc++ 的uniform_int_distribution 的算法基本相同。）

【参考方案1】：

这里有隐藏的深度：

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

% 6 的使用可以（但在我见过的std::rand 的每个实现中不会）在上面引入任何额外的统计偏差并超越提出的替代方案。 % 6 是危险的这种情况是数字生成器在低位具有相关性的情况，例如我认为在 1970 年代一个相当著名的 rand 的 IBM 实现（用 C 语言）低位作为“最后的繁荣”。进一步的考虑是 6 非常小 cf。 RAND_MAX，所以如果RAND_MAX 不是 6 的倍数，那么影响很小。

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

【讨论】：

% 6 只要rand() 生成的不同值的数量不是 6 的倍数，就会产生有偏差的结果。鸽子洞原理。当然，当RAND_MAX 远大于 6 时，偏差很小，但它就在那里。对于更大的目标范围，效果当然更大。

@PeteBecker：确实，我应该说清楚。但请注意，由于整数除法截断效应，当您的采样范围接近 RAND_MAX 时，您也会陷入困境。

@Bathsheba 截断效果不会导致结果大于 6 从而重复执行整个操作吗？

@Gerhardh：正确。事实上，它完全导致x==7。基本上，您将[0, RAND_MAX] 范围划分为 7 个子范围，其中 6 个大小相同，最后一个较小的子范围。最后一个子范围的结果被丢弃。很明显，您不能以这种方式在末尾有两个较小的子范围。

@MSalters：确实。但请注意，由于截断，另一种方式仍然受到影响。我的假设是，由于统计陷阱更难理解，人们更喜欢后者！

【参考方案2】：

这个示例代码说明std::rand 是一个传统的货物***的例子，每次看到它都会让你的眉毛扬起。

这里有几个问题：

人们通常假设的合同——即使是那些不了解任何事情并且不会用这些术语来思考它的可怜的倒霉蛋——是rand 来自均匀分布的样本对 0、1、2、...、RAND_MAX 中的整数进行处理，每次调用都会产生一个独立样本。

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

rand 函数计算 0 到 RAND_MAX 范围内的伪随机整数序列。

这是一个没有意义的句子，因为伪随机性是一个函数（或函数族）的属性，而不是一个整数，但它甚至不会停止ISO 官员不要滥用该语言。毕竟，唯一会对此感到不安的读者知道最好不要阅读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;

这有一个不幸的特性，即即使单个样本可能均匀分布在一个均匀的随机种子（这取决于RAND_MAX 的具体值）下，它在偶数和奇数之间交替在连续调用中——之后

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

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

第二个问题是，即使每个调用都在 0、1、2、...、RAND_MAX 上独立于均匀随机分布进行采样，rand() % 6 的结果也会不会像掷骰子一样均匀分布在 0、1、2、3、4、5 中，除非RAND_MAX 与 -1 模 6 一致。 简单反例：如果 RAND_MAX = 6，则从 @ 987654337@，所有结果的概率为 1/7，但从rand() % 6 开始，结果 0 的概率为 2/7，而所有其他结果的概率为 1/7。

正确的做法是拒绝抽样： 重复从 0、1、2、...、RAND_MAX 中抽取一个独立的均匀随机样本 s , 和 reject（例如）结果 0、1、2、...、((RAND_MAX + 1) % 6) - 1——如果你得到其中之一，重新开始；否则，产生s % 6。

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

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

您拒绝rand() 的哪些结果的选择无关紧要，只要您将相等数量的结果映射到每个小于 6 的整数。cppreference.com 上的代码生成了一个 不同的选择，因为上面的第一个问题——rand() 的输出的分布或独立性没有任何保证，并且在实践中，低位显示的模式“看起来不够随机”（没关系，下一个输出是前一个输出的确定性函数）。

读者练习：如果rand() 在 0、1、2、...、RAND_MAX 上产生均匀分布，则证明 cppreference.com 上的代码在骰子上产生均匀分布。

读者练习：为什么您更喜欢其中一个或其他子集来拒绝？两种情况下每次试验需要什么计算？

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

您可以走花哨的过度设计路线和 C++11 的 std::uniform_int_distribution 类，使用适当的随机设备和您最喜欢的随机引擎（如广受欢迎的 Mersenne twister std::mt19937）与您四岁的孩子玩骰子表弟，但即使这样也不适合生成加密密钥材料——而且 Mersenne twister 也是一个可怕的空间占用者，它具有数 KB 的状态，以令人讨厌的设置时间对你的 CPU 缓存造成严重破坏，所以即使它也很糟糕对于，例如，具有可重现的子计算树的并行蒙特卡罗模拟；它的受欢迎程度可能主要源于其朗朗上口的名字。但是你可以像这个例子一样用它来掷骰子！

另一种方法是使用具有小状态的简单加密伪随机数生成器，例如简单的fast key erasure PRNG，或者如果您有信心（例如 em>，在自然科学研究的蒙特卡洛模拟中）如果状态受到损害，预测过去的结果不会产生不利后果。

【讨论】：

"一个淫秽的设置时间" 无论如何，您实际上不应该使用多个随机数生成器（每个线程），因此除非您的程序运行时间不长，否则设置时间将被摊销。

Downvote BTW 因为不理解问题中的循环正在执行完全相同的拒绝采样，具有完全相同的 (RAND_MAX + 1 )% 6 值。 如何 细分可能的结果并不重要。您可以从[0, RAND_MAX) 范围内的任何地方拒绝它们，只要接受范围的大小是 6 的倍数。地狱，您可以完全拒绝任何结果x>6，并且您不需要@987654358 @ 了。

我对这个答案不太满意。咆哮可能很好，但你把它带向了错误的方向。例如，您抱怨“更好的随机性”不是一个技术术语，它毫无意义。这是对的一半。是的，这不是一个技术术语，但在上下文中它是一个非常有意义的简写。暗示使用这样一个术语的用户要么无知，要么恶意，这本身就是其中之一。 “良好的随机性”可能很难准确定义，但很容易掌握函数何时产生具有更好或更差随机性属性的结果。

我喜欢这个答案。这有点咆哮，但它有很多很好的背景信息。请记住，真正的专家只使用硬件随机发生器，问题就是这么难。

对我来说正好相反。虽然它确实包含了很好的信息，但除了意见之外，它太过咆哮了。除了有用性。

【参考方案3】：

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

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

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

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

【参考方案4】：

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

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

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

作为一个极端的例子，假设rand() 在[0..6] 范围内产生均匀分布的值。如果查看这些值的余数，当rand() 返回[0..5] 范围内的值时，余数会在[0..5] 范围内产生均匀分布的结果。当rand() 返回 6 时，rand() % 6 返回 0，就像 rand() 返回 0 一样。因此，您得到的分布中 0 的数量是任何其他值的两倍。

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

避免该问题的方法是丢弃会产生不一致重复的值。您计算小于或等于 RAND_MAX 的 6 的最大倍数，并且每当 rand() 返回一个大于或等于该倍数的值时，您拒绝它并再次调用 `rand()，根据需要多次调用。

所以：

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

这是所讨论代码的不同实现，旨在更清楚地显示正在发生的事情。

【讨论】：

我已经承诺至少有一个常客在这个网站上写一篇关于这个的论文，但我认为抽样和拒绝可以让高光时刻;例如过度夸大方差。

如果 rand_max 为 32768（在某些实现中是这样），我做了一张图表，说明该技术引入了多少偏差。 ericlippert.com/2013/12/16/…

@Bathsheba：确实，某些拒绝函数可能会导致这种情况，但这种简单的拒绝会将统一 IID 转换为不同的统一 IID 分布。没有比特结转，如此独立，所有样本都使用相同的拒绝，如此相同，并且微不足道以显示一致性。而均匀积分随机变量的高阶矩完全由其范围定义。

@MSalters：您的第一句话对于 true 生成器来说是正确的，对于伪生成器来说不一定是正确的。等我退休了，我会写一篇关于这个的论文。

@Anthony 用骰子来思考。你想要一个介于 1 和 3 之间的随机数，而你只有一个标准的 6 面骰子。如果您掷出 4-6，您只需减去 3 即可得到。但是，假设您想要一个介于 1 和 5 之间的数字。如果在掷 6 时减去 5，那么最终得到的 1 是任何其他数字的两倍。这基本上就是 cppreference 代码正在做的事情。正确的做法是重新滚动 6s。这就是皮特在这里所做的：将骰子分开，以便有相同数量的方式来滚动每个数字，并重新滚动任何不适合偶数除法的数字

【参考方案5】：

无论如何，我都不是经验丰富的 C++ 用户，但有兴趣看看其他答案是否与 std::rand()/((RAND_MAX + 1u)/6) 比 1+std::rand()%6 更少偏见实际上是正确的。所以我写了一个测试程序来将这两种方法的结果制成表格（我已经很久没有写过C++了，请检查一下）。运行代码的链接位于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 次。

更新：当我写我的答案时，康拉德鲁道夫发布了一个采用相同方法的答案，但得到的结果却截然不同。我没有评论他的答案的声誉，所以我将在这里解决它。首先，主要的是他使用的代码每次运行时都使用相同的随机数生成器种子。如果你改变种子，你实际上会得到各种各样的结果。其次，如果你不改变种子，而是改变试验次数，你也会得到各种各样的结果。尝试增加或减少一个数量级以了解我的意思。第三，在预期值不太准确的情况下，会进行一些整数截断或舍入。可能不足以产生影响，但它就在那里。

总的来说，他只是碰巧获得了正确的种子和试验次数，他可能会得到错误的结果。

【讨论】：

由于您的误解，您的实现包含一个致命缺陷：引用的段落不是将rand()%6 与rand()/(1+RAND_MAX)/6 进行比较。相反，它将直接取余数与 拒绝采样 进行比较（有关解释，请参阅其他答案）。因此，您的第二个代码是错误的（while 循环什么也不做）。你的统计测试也有问题（你不能只是重复你的测试来测试稳健性，你没有进行校正，......）。

@KonradRudolph 我没有代表对您的答案发表评论，因此我将其添加为我的更新。你的也有一个致命的缺陷，它碰巧使用了一个固定的种子和每次运行的试验次数，这会给出错误的结果。如果你用不同的种子重复运行，你可能已经抓住了。但是是的，你是对的，while 循环什么都不做，但它也不会改变那个特定代码块的结果

我确实重复了，实际上。由于使用std::srand 设置了随机种子（并且没有使用<random>）is quite hard to do in a standards conforming way，因此故意不设置种子，并且我不希望它的复杂性影响其余代码。它也与计算无关：在模拟中重复相同的序列是完全可以接受的。当然，不同的种子会产生不同的结果，有些结果并不显着。这完全取决于 p 值的定义方式。

老鼠们，我在重复时犯了一个错误；你是对的，重复运行的第 95 个分位数非常接近 p=0.05——即，正是我们在 then null 时所期望的。总而言之，我的 std::rand 标准库实现在随机种子范围内为 d6 产生了非常好的抛硬币模拟。

统计上的意义只是故事的一部分。您有一个原假设（均匀分布）和一个备择假设（模偏差）——实际上，一系列备择假设，由RAND_MAX 的选择索引，它决定了模的效果大小偏见。统计显着性是在原假设下你错误地拒绝它的概率。 统计功效是多少？在替代假设下，您的检验正确拒绝原假设的概率是多少？当 RAND_MAX = 2^31 - 1 时，你会以这种方式检测 rand() % 6 吗？