为什么人们说在使用随机数发生器时存在模数偏差？
我已经看到这个问题很多，但从来没有看到一个真正的具体答案。 所以我打算在这里发布一个，希望能够帮助人们理解为什么在使用随机数生成器（如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，而不是使用模数，丢弃结果产生110和111 ，只是再试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％。
这引出了另外两个问题：
- 如果我们添加足够的比特，是否可以保证丢弃的概率会降低？
- 一般情况下有多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; }
我鼓励使用MODULUS和ROLLS值来查看在大多数情况下实际发生了多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是当你想要一个0到1之间的随机数时你不应该做的。 您可以通过执行% 3来获得0到2之间的随机数，因为在这种情况下： RAND_MAX是3的倍数。
另一种方法
有更简单的，但添加到其他答案，这里是我的解决scheme获得0和n - 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。;
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; } } }
