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

我已经看到这个问题很多,但从来没有看到一个真正的具体答案。 所以我打算在这里发布一个,希望能够帮助人们理解为什么在使用随机数生成器(如C ++中的rand()时确实存在“模偏差”。

所以rand()是一个伪随机数发生器,它select一个介于0和RAND_MAX之间的自然数, RAND_MAX是一个在cstdlib定义的常量(参见本文以获得有关rand()的一般性概述)。

现在如果你想产生一个0到2之间的随机数,会发生什么? 为了说明起见,假设RAND_MAX是10,我决定通过调用rand()%3来生成一个介于0和2之间的随机数。 但是, rand()%3不会以相等的概率产生0到2之间的数字!

rand()返回rand()或9时, rand()%3 == 0 。 因此,P(0)= 4/11

rand()返回rand()或10时, rand()%3 == 1 。 因此,P(1)= 4/11

rand()返回2,5或8时, rand()%3 == 2 。 因此,P(2)= 3/11

这不会以相等的概率产生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小值,这很less需要对rand()多次调用。


作品引用和进一步阅读:

  • CPlusPlus参考

  • 永远Confuzzled


保持随机select是消除偏见的好方法。

更新

如果我们search一个可被n整除的x,我们可以使代码更快。

 // Assumptions // rand() in [0, RAND_MAX] // n in (0, RAND_MAX] int x = rand(); // Keep searching for an x in a range divisible by n while (x >= RAND_MAX - (RAND_MAX % n)) { x = rand(); } x %= n; 

上面的循环应该非常快,平均说1次迭代。

@ user1413793对于这个问题是正确的。 我不打算进一步讨论,除了要说一点:是的,对于小的n值和RAND_MAX大值,模数偏差可能非常小。 但是使用偏差诱导模式意味着每次计算随机数时必须考虑偏差,并针对不同的情况select不同的模式。 如果你做出错误的select,它引入的错误是微妙的,unit testing几乎是不可能的。 相比于只使用适当的工具(如arc4random_uniform ),这是额外的工作,而不是更less的工作。 做更多的工作,并得到更糟的解决scheme是糟糕的工程,尤其是在大多数平台上,每一次都很简单。

不幸的是,解决scheme的实现都是不正确的或效率低于他们应该。 (每个解决scheme都有不同的解释问题的意见,但是没有一个解决scheme已经解决了这些问题。)这可能会使偶然的答案search者感到困惑,所以我在这里提供了一个很好的实现。

同样,最好的解决scheme就是在提供平台的平台上使用arc4random_uniform ,或者为您的平台使用类似的范围解决scheme(例如Java上的Random.nextInt )。 它会做你没有代码成本的正确的事情。 这几乎总是正确的要求。

如果您没有arc4random_uniform ,那么您可以使用opensource的function来查看它是如何在更宽范围的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; } 

对于那些需要实现类似的东西的人来说,值得注意的是这个代码的最新的提交评论:

改变arc4random_uniform()来计算2**32 % upper_bound'' as -upper_bound%upper_bound''。 简化代码,使其在ILP32和LP64体系结构上保持一致,而在LP64体系结构上使用32位余数而不是64位余数,也稍快一些。

由Jorden Verwer在tech @ ok deraadt指出; 没有djm或奥托的反对意见

Java实现也很容易find(请参阅前面的链接):

 public int nextInt(int n) { if (n <= 0) throw new IllegalArgumentException("n must be positive"); if ((n & -n) == n) // ie, 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; } 

定义

模偏置Modulo Bias)是使用模运算将输出设置为input集的子集的固有偏差。 一般来说,只要input和输出集合之间的映射不是均等分布,就存在一个偏差,例如当输出集合的大小不是input集合的大小除数时使用模运算的情况。

这种偏差在计算中特别难以避免,其中数字表示为位串: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产生111产生1

潜在的解决scheme

方法0:

从理论上说,可以雇佣一支小军队整天掷骰子,并将结果logging到数据库中,然后每个结果只使用一次。 这听起来很实际,而且很可能不会产生真正随机的结果(双关意图)。

方法1:

一个天真但在math上正确的解决scheme,而不是使用模数,丢弃结果产生110111 ,只是再试3个新的比特。 不幸的是,这意味着每个卷筒25%的机会需要重卷,包括每个重卷本身。 除了最琐碎的用途之外,这显然是不切实际的。

方法2:

使用更多的位:而不是3位,使用4.这会产生16个可能的结果。 当然,再次滚动结果大于5会使情况变得更糟(10/16 = 62.5%),这样一来就不会有帮助。

请注意,2 * 6 = 12 <16,所以我们可以安全地采取任何小于12的结果,并减less模6来均匀分配结果。 其他4个结果必须被丢弃,然后像以前的方法那样重新滚动。

听起来不错,但让我们来看看math:

 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. 一般情况下有多less位是足够的?

通用解决scheme

谢天谢地,第一个问题的答案是肯定的。 6的问题是在2和4之间的2 ^ x mod 6翻转巧合地是彼此的倍数,因此对于偶数x> 1,

 [2^x mod 6] / 2^x == [2^(x+1) mod 6] / 2^(x+1) 

因此,6是一个例外,而不是规则。 有可能find更大的模数,以相同的方式产生2的连续幂数,但最终必须环绕,并且丢弃的概率将减小。

如果没有提供进一步的证据,通常使用双倍的比特数将提供更小的,通常不重要的抛弃的机会。

概念validation

这是一个使用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; uint64_t randomPool = 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值来查看在大多数情况下实际发生了多less次重播。 持怀疑态度的人可能也希望将计算值保存到文件中,并validation分布是否正常。

有两种常见的使用模数的投诉。

  • 一个对所有发电机都有效。 在极限情况下更容易看到。 如果你的发生器的RAND_MAX是2(不符合C标准),而你只需要0或1作为值,那么使用模数将产生0两次(当发生器产生0和2时),因为它会产生1(当发生器产生1)。 请注意,一旦您不删除值,就会发生这种情况,无论您使用从生成器值到想要的映射的映射,其中一个会出现两次。

  • 某种发生器的低有效位比其他参数less,至less对于它们的一些参数,但可悲的是这些参数具有其他有趣的特性(例如能够使RAND_MAX小于2的幂)。 这个问题是众所周知的,很长一段时间的库实现可能会避免这个问题(例如C标准中的样例rand()实现使用这种types的生成器,但删除了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引入了标准的方式来减less和其他发电机比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

循环更好的解决scheme(这是非常低效的,甚至不应该被build议)是使用具有更大输出范围的PRNG。 Mersenne Twisteralgorithm最大输出量为4,294,967,295。 因此,对于所有意图和目的, MersenneTwister::genrand_int32() % 10将被平均分配,模数偏差效应将几乎消失。

RAND_MAX值为3 (实际上应该比这更高,但偏差仍然存在),从这些计算中可以看出,存在一个偏差:

1 % 2 = 1 2 % 2 = 0 3 % 2 = 1 random_between(1, 3) % 2 = more likely a 1

在这种情况下, % 2是当你想要一个01之间的随机数时你不应该做的。 您可以通过执行% 3来获得02之间的随机数,因为在这种情况下: RAND_MAX3的倍数。

另一种方法

有更简单的,但添加到其他答案,这里是我的解决scheme获得0n - 1之间的随机数,所以n种不同的可能性,没有偏见。

  • 编码可能数量所需的比特数(不是字节)是你需要的随机数据的比特数
  • 从随机位编码数字
  • 如果这个数字是>= n ,重新启动(不模)。

真正的随机数据是不容易获得的,所以为什么要使用比需要更多的比特。

以下是Smalltalk中的一个例子,使用来自伪随机数生成器的位caching。 我不是安全专家,所以请自担风险。

 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 

马克的解决scheme(公认的解决scheme)是几乎完美的。

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

编辑16年3月25日在23:16

Mark Amery 39k21170211

然而,在RAND_MAX(RM)的值小于N的倍数的所有情况下,它丢弃一个有效的结果集被舍弃。

即,当作为无效(I)被丢弃的值的数量等于N时,它们实际上是有效集合,而不是无效集合。

例如:

 RM = 255 N = 4 Discard X => RM - RM % N When X => 252, Discarded values = 252, 253, 254, 255 Number of discarded Values (I) = RM % N + 1 

正如您在示例中所看到的,丢弃值的计数= 4,当丢弃值的计数= 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% 

要否定这一点,我们可以做一个简单的修改如下所示:

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

这提供了一个更一般的公式版本,它说明了使用模量来定义最大值的额外特性。

RAND_MAX是N的一个乘数的小数值的例子

Mark'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. 

修正版本:

 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语句为你的最终乘数。 但也许你有代码,可能有一个有效的理由返回1时,函数调用n = 1 …

所以最好使用0,当你希望有n = RAND_MAX + 1时,通常会提供Div 0错误

IE:

 int x; if n != 0 { do { x = rand(); } while (x > (RAND_MAX - RAND_MAX % n)); x %= n; } else { x = rand(); } 

这两个解决scheme解决了这个问题,当RM + 1是n的乘积时,将会发生不必要的丢弃的有效结果。

第二个版本也包含边缘情况,当你需要n等于RAND_MAX中包含的全部可能值时。

两者的修改方法是相同的,并且允许提供更一般的解决scheme以提供有效的随机数并最小化丢弃的值。

重申:

扩展标记示例的基本通用解决scheme:

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

允许一个额外的RAND_MAX + 1 = n的场景的扩展通用解决scheme:

 int x; if n != 0 { do { x = rand(); } while (x > (RAND_MAX - RAND_MAX % n)); x %= n; } else { x = rand(); } 

我刚刚写了冯·诺依曼的无偏硬币翻转法的代码,理论上应该消除随机数生成过程中的任何偏差。 更多信息可以在http://en.wikipedia.org/wiki/Fair_coinfind&#x3002;

 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; } } }