为什么 rand()%6 有偏见

Why is rand()%6 biased?

本文关键字:有偏见 rand 为什么      更新时间:2023-10-16

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

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

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

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

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

这里有隐藏的深度:

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

  2. 使用% 6可以(但在我见过的每个std::rand实现中都没有)引入任何超出所提出替代方案的额外统计偏差。% 6危险的例子是数字生成器在低阶位中具有相关平原的情况,例如我认为在 1970 年代相当著名的 IBM 实现(C 语言)rand,它将高位和低位翻转为"最终繁荣"。进一步的考虑是 6 非常小 cf。RAND_MAX,所以如果RAND_MAX不是 6 的倍数,效果将很小,而它可能不是。

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

此示例代码说明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,则从rand()开始,所有结果的概率相等 1/7,但从rand() % 6开始,结果 0 的概率为 2/7,而所有其他结果的概率为 1/7。

正确的方法是使用拒绝抽样:0、1、2、...、RAND_MAX重复绘制一个独立的均匀随机样本s,并拒绝(例如)结果 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的每个可能结果都是通过相同数量的rand()接受结果获得的,所以如果rand()均匀分布,那么s也是如此。 试验次数没有限制,但预期数量小于 2,成功的概率随着试验次数呈指数级增长。

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

读者练习:证明 cppreference.com 处的代码在模具辊上产生均匀分布,如果rand()在 0、1、2、...、RAND_MAX上产生均匀分布。

读者练习:为什么你更喜欢一个或另一个子集来拒绝? 在这两个案件中,每个试验需要什么计算?

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

你可以走花哨的过度设计的路线和C++11的std::uniform_int_distribution类,使用适当的随机设备和你最喜欢的随机引擎,如广受欢迎的梅森捻线机std::mt19937和你四岁的表弟玩骰子,但即使这样也不适合生成加密密钥材料——梅森捻线机也是一个可怕的太空猪,一个数千字节的状态对你的CPU缓存造成严重破坏,一个淫秽的设置时间,因此即使对于具有可重现的子计算树的并行蒙特卡罗模拟也是如此;它的受欢迎程度可能主要源于其朗朗上口的名字。 但是你可以用它来掷玩具骰子,就像这个例子一样!

另一种方法是使用具有小状态的简单加密伪随机数生成器,例如简单的快速密钥擦除PRNG,或者只是流密码,例如AES-CTR或ChaCha20,如果您确信(例如,在自然科学研究的蒙特卡罗模拟中)如果状态受到损害,预测过去的结果不会产生不利影响。

可以将随机数生成器视为处理二进制数字流。生成器通过将流切成块来将流转换为数字。如果std:rand函数使用 32767 的RAND_MAX,则在每个切片中使用 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这样的线性全等发电机以最低质量提供最快的速度。加密生成器以最低速度提供最高质量。

无论如何,我都不是一个有经验的C++用户,但有兴趣看看其他答案是否关于std::rand()/((RAND_MAX + 1u)/6)1+std::rand()%6实际的偏见更少。所以我写了一个测试程序来列出两种方法的结果(我已经很久没有写C++了,请检查一下)。此处提供了运行代码的链接。还复制如下:

// 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函数运行卡方检验,以查看结果是否与预期明显不同。这个堆栈交换问题更详细地介绍了使用卡方检验来测试骰子公平性:如何测试骰子是否公平?以下是几次运行的结果:

> ?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 次。

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

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