为什么人们说在使用随机数生成器时存在模偏置?

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

65647 次浏览

所以rand()是一个伪随机数生成器,它在0和RAND_MAX之间选择一个自然数,RAND_MAXcstdlib中定义的常量(有关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 < em > < / em >

这不会以相等的概率生成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()进行多次调用。


引用作品及进一步阅读:


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

更新

如果我们在能被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次迭代。

对于模的使用,有两种常见的抱怨。

  • 1对所有生成器有效。在极限情况下更容易看出来。如果你的发电机有RAND_MAX 2(不符合C标准),你希望只有0或1的值,使用模将生成0 2倍(当生成器生成0和2),它将生成1(当生成器生成1)。请注意,这是真的只要你不下降值,无论您使用的是发电机的值映射到想要的,人会出现另一个的两倍。

  • 少一些发电机有明显比另一方更随机,至少他们的一些参数,但遗憾的是这些参数有其他有趣的特点(RAND_MAX能够有一个小于2的乘方)。问题是众所周知的,可能很长一段时间内库实现避免问题(例如C标准的样本rand()实现使用这种发电机,但下降16不那么重要部分),但有些人喜欢抱怨,你可能会有坏运气

使用类似于

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()以外的其他生成器。

正如接受的答案所示,“模偏置”的根源在于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很小时,才会如此明显地显示。或者更具体地说,当你正在修改的数字比RAND_MAX大时。

一个比循环更好的解决方案(它效率低得离谱,甚至不应该被建议)是使用一个输出范围大得多的PRNG。梅森素数捻线机算法的最大输出为4,294,967,295。这样一来,无论出于何种目的,执行MersenneTwister::genrand_int32() % 10都将是均匀分布的,模偏效应将几乎消失。

@user1413793是正确的。我不打算进一步讨论这个问题,除了说明一点:是的,对于n的小值和RAND_MAX的大值,模偏置可以非常小。但是使用偏差诱导模式意味着每次计算随机数时都必须考虑偏差,并为不同的情况选择不同的模式。如果你做出了错误的选择,它所带来的错误就会非常微妙,几乎不可能进行单元测试。与仅仅使用适当的工具(如arc4random_uniform)相比,这是额外的工作,而不是减少工作。做更多的工作却得到更糟糕的解决方案是糟糕的工程,尤其是在大多数平台上,每次都很容易做到正确的时候。

不幸的是,解决方案的实现都是不正确的,或者效率低于应有的水平。(每个解决方案都有各种解释问题的评论,但没有一个解决方案被修复以解决这些问题。)这可能会让那些随意寻求答案的人感到困惑,所以我在这里提供了一个已知的良好实现。

同样,最好的解决方案是在提供arc4random_uniform的平台上使用它,或者为您的平台使用类似的远程解决方案(例如Java上的Random.nextInt)。它将在没有代码成本的情况下做正确的事情。这几乎总是正确的选择。

如果你没有arc4random_uniform,那么你可以使用开源的力量来查看它是如何在更广泛的RNG上实现的(在这种情况下是ar4random,但类似的方法也可以在其他RNG上工作)。

下面是OpenBSD的实现:

/*
* Calculate a uniformly distributed random number less than upper_bound
* avoiding "modulo bias".
*
* Uniformity is achieved by generating new random numbers until the one
* returned is outside the range [0, 2**32 % upper_bound).  This
* guarantees the selected random number will be inside
* [2**32 % upper_bound, 2**32) which maps back to [0, upper_bound)
* after reduction modulo upper_bound.
*/
u_int32_t
arc4random_uniform(u_int32_t upper_bound)
{
u_int32_t r, min;


if (upper_bound < 2)
return 0;


/* 2**32 % x == (2**32 - x) % x */
min = -upper_bound % upper_bound;


/*
* This could theoretically loop forever but each retry has
* p > 0.5 (worst case, usually far better) of selecting a
* number inside the range we need, so it should rarely need
* to re-roll.
*/
for (;;) {
r = arc4random();
if (r >= min)
break;
}


return r % upper_bound;
}

对于那些需要实现类似事情的人来说,值得注意这段代码上的最新commit注释:

更改arc4random_uniform()计算2**32 % upper_bound-upper_bound % upper_bound。简化代码并使之成为 在ILP32和LP64架构上都是一样的,而且速度也略快 LP64架构使用32位余数而不是64位余数 余数。< / p >

由Jorden Verwer在tech@上指出 好的deraadt;DJM或otto

无异议

Java实现也很容易找到(见之前的链接):

public int nextInt(int n) {
if (n <= 0)
throw new IllegalArgumentException("n must be positive");


if ((n & -n) == n)  // i.e., n is a power of 2
return (int)((n * (long)next(31)) >> 31);


int bits, val;
do {
bits = next(31);
val = bits % n;
} while (bits - val + (n-1) < 0);
return val;
}

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

定义

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

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

而不是使用模量,一个简单但数学上正确的解决方案是丢弃产生110111的结果,并简单地重新尝试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%

确实有了改进,但在许多实际情况下还不够好。好消息是,增加更多的比特永远不会增加需要丢弃和重新滚动的机会。这不仅适用于骰子,而且适用于所有情况。

事实上,如果我们将摇到6位,概率仍然是6.25%。

这就引出了另外两个问题:

  1. 如果我们添加足够多的比特,是否能保证丢弃的概率会降低?
  2. 多少位才够呢在一般情况下?

通解

幸运的是,第一个问题的答案是肯定的。6的问题在于,2^x mod 6在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;
}

我鼓励使用MODULUSROLLS值,看看在大多数情况下实际发生了多少次重新滚动。持怀疑态度的人也可能希望将计算值保存到文件中,并验证分布是否正常。

对于RAND_MAX3(实际上它应该比这个值高得多,但偏差仍然存在),从这些计算中可以理解存在偏差:

< p > 1 % 2 = 1 2 % 2 = 0 3 % 2 = 1 random_between(1, 3) % 2 = more likely a 1 < / p >

在这种情况下,当你想要一个介于01之间的随机数时,不应该使用% 2。你可以通过% 3得到一个介于02之间的随机数,因为在这种情况下:RAND_MAX3的倍数。

另一种方法

有一个更简单的答案,但要添加到其他答案,这里是我的解决方案,以获得0n - 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

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

int x;


do {
x = rand();
} while (x >= (RAND_MAX - RAND_MAX % n));


x %= n;

编辑于2016年3月25日23:16

Mark Amery 39k21170211

然而,它有一个警告,在任何场景中,如果RAND_MAX (RM)比N的倍数小1(其中N =可能的有效结果的数量),它会丢弃1组有效结果。

例如,当'count of values discarded' (D)等于N时,它们实际上是一个有效集(V),而不是一个无效集(I)。

造成这种情况的原因是Mark在某些时候忽略了NRand_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不等于质数时,这是有效的关注)

例如:

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的乘法。

马克'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() % n,在已接受的答案中讨论过)会这样引入偏差,而且“乘-移”也会这样引入偏差。Daniel Lemire的还原,或者如果您在一组迭代之后不再拒绝某个结果。(需要明确的是,这并不意味着没有办法修复伪随机生成器中存在的偏差问题。例如,即使模和其他约简通常是有偏差的,如果可能的整数的范围是2 而且的幂,如果随机生成器产生无偏置的随机比特或块,则它们不会有偏差问题。)

这个答案的其余部分将显示随机生成器中运行时间和偏差之间的关系。从这里开始,我们将假设我们有一个“真实的”;随机发生器,可以产生无偏和独立的随机位

1976年,D. E. Knuth和a . C. Yao表明,任何产生具有给定概率的随机整数的算法,只使用随机位,都可以表示为二叉树,其中随机位表示遍历树的路径,每个叶子(端点)对应一个结果。在这种情况下,我们处理的是生成[0,n)随机整数的算法,其中每个整数的选择概率为1/n。如果所有结果的树中出现相同数量的叶子,则算法为无偏见的。但是如果1/n有一个不终止的二进制展开(当n不是2的幂时就是这种情况),算法只有在-时才会是无偏的

  • 二叉树有一个“无穷”;深度,或
  • 二叉树包括“拒绝”;叶子在最后,

无论哪种情况,算法都不会在常数时间内运行在最坏的情况下会一直运行下去。(另一方面,当n是2的幂时,最优二叉树将具有有限深度并且没有拒绝节点。)

二叉树的概念还表明,任何“修正”的方法;这种最坏情况下的时间复杂度通常会导致偏差。(再次强调,这并不意味着没有办法修复伪随机生成器中存在的偏差问题。)例如,模约相当于一个二叉树,其中拒绝叶被标记的结果取代——但由于可能的结果比拒绝叶更多,只有一些结果可以取代拒绝叶,从而引入偏差。如果您在一组迭代之后停止拒绝,则会得到相同类型的二叉树和相同类型的偏差。(然而,根据应用程序的不同,这种偏差可以忽略不计。随机整数生成也有安全方面的问题,这太复杂了,不能在这个答案中讨论。)

为了说明这一点,下面的JavaScript代码实现了一个随机整数算法快速掷骰子机,由J. Lumbroso(2013)编写。请注意,它包括一个拒绝事件和一个循环,这是在一般情况下使算法无偏倚所必需的。

function randomInt(minInclusive, maxExclusive) {
var maxInclusive = (maxExclusive - minInclusive) - 1
var x = 1
var y = 0
while(true) {
x = x * 2
var randomBit = (Math.random() < 0.5 ? 0 : 1)
y = y * 2 + randomBit
if(x > maxInclusive) {
if (y <= maxInclusive) { return y + minInclusive }
// Rejection
x = x - maxInclusive - 1
y = y - maxInclusive - 1
}
}
}

请注意

*这个答案不会涉及C语言中的rand()函数,因为它是有很多问题。也许这里最严重的是,C标准没有显式地为rand()返回的数字指定特定的分布,甚至不是均匀的分布。