我看到很多人问过这个问题,但从未见过一个真正具体的答案。所以我打算在这里发布一个,希望能帮助人们理解为什么在使用随机数生成器时,比如c++中的rand(),会有“模偏置”。
当前回答
@user1413793 is correct about the problem. I'm not going to discuss that further, except to make one point: yes, for small values of n and large values of RAND_MAX, the modulo bias can be very small. But using a bias-inducing pattern means that you must consider the bias every time you calculate a random number and choose different patterns for different cases. And if you make the wrong choice, the bugs it introduces are subtle and almost impossible to unit test. Compared to just using the proper tool (such as arc4random_uniform), that's extra work, not less work. Doing more work and getting a worse solution is terrible engineering, especially when doing it right every time is easy on most platforms.
不幸的是,解决方案的实现都是不正确的,或者效率低于应有的水平。(每个解决方案都有各种解释问题的评论,但没有一个解决方案被修复以解决这些问题。)这可能会让那些随意寻求答案的人感到困惑,所以我在这里提供了一个已知的良好实现。
同样,最好的解决方案是在提供arc4random_uniform的平台上使用它,或者为您的平台使用类似的远程解决方案(如Random。nextInt在Java)。它将在没有代码成本的情况下做正确的事情。这几乎总是正确的选择。
如果你没有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位余数 余数。 由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;
}
其他回答
@user1413793 is correct about the problem. I'm not going to discuss that further, except to make one point: yes, for small values of n and large values of RAND_MAX, the modulo bias can be very small. But using a bias-inducing pattern means that you must consider the bias every time you calculate a random number and choose different patterns for different cases. And if you make the wrong choice, the bugs it introduces are subtle and almost impossible to unit test. Compared to just using the proper tool (such as arc4random_uniform), that's extra work, not less work. Doing more work and getting a worse solution is terrible engineering, especially when doing it right every time is easy on most platforms.
不幸的是,解决方案的实现都是不正确的,或者效率低于应有的水平。(每个解决方案都有各种解释问题的评论,但没有一个解决方案被修复以解决这些问题。)这可能会让那些随意寻求答案的人感到困惑,所以我在这里提供了一个已知的良好实现。
同样,最好的解决方案是在提供arc4random_uniform的平台上使用它,或者为您的平台使用类似的远程解决方案(如Random。nextInt在Java)。它将在没有代码成本的情况下做正确的事情。这几乎总是正确的选择。
如果你没有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位余数 余数。 由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;
}
因此rand()是一个伪随机数生成器,它在0和RAND_MAX之间选择一个自然数,RAND_MAX是cstdlib中定义的一个常量(有关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
这不会以相等的概率生成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()。
引用作品及进一步阅读:
CPlusPlus参考 永远Confuzzled
模约化是一种常见的方法,可以使随机整数生成器避免永远运行的最坏情况。
When the range of possible integers is unknown, however, there is no way in general to "fix" this worst case of running forever without introducing bias. It's not just modulo reduction (rand() % n, discussed in the accepted answer) that will introduce bias this way, but also the "multiply-and-shift" reduction of Daniel Lemire, or if you stop rejecting an outcome after a set number of iterations. (To be clear, this doesn't mean there is no way to fix the bias issues present in pseudorandom generators. For example, even though modulo and other reductions are biased in general, they will have no issues with bias if the range of possible integers is a power of 2 and if the random generator produces unbiased random bits or blocks of them.)
这个答案的其余部分将显示随机生成器中运行时间和偏差之间的关系。从这里开始,我们将假设我们有一个“真正的”随机生成器,可以产生无偏和独立的随机比特
In 1976, D. E. Knuth and A. C. Yao showed that any algorithm that produces random integers with a given probability, using only random bits, can be represented as a binary tree, where random bits indicate which way to traverse the tree and each leaf (endpoint) corresponds to an outcome. In this case, we're dealing with algorithms that generate random integers in [0, n), where each integer is chosen with probability 1/n. The algorithm is unbiased if the same number of leaves appear in the tree for all outcomes. But if 1/n has a non-terminating binary expansion (which will be the case if n is not a power of 2), the algorithm will be unbiased only if—
二叉树具有“无限”深度,或者 二叉树的末端包含“排斥”叶,
无论哪种情况,算法都不会在常数时间内运行在最坏的情况下会一直运行下去。(另一方面,当n是2的幂时,最优二叉树的深度是有限的,并且没有拒绝节点。)
The binary tree concept also shows that any way to "fix" this worst-case time complexity will lead to bias in general. (Again, this doesn't mean there is no way to fix the bias issues present in pseudorandom generators.) For instance, modulo reductions are equivalent to a binary tree in which rejection leaves are replaced with labeled outcomes — but since there are more possible outcomes than rejection leaves, only some of the outcomes can take the place of the rejection leaves, introducing bias. The same kind of binary tree — and the same kind of bias — results if you stop rejecting after a set number of iterations. (However, this bias may be negligible depending on the application. There are also security aspects to random integer generation, which are too complicated to discuss in this answer.)
为了说明这一点,下面的JavaScript代码实现了J. Lumbroso(2013)提出的名为Fast Dice Roller的随机整数算法。请注意,它包括一个拒绝事件和一个循环,这是在一般情况下使算法无偏倚所必需的。
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
}
}
}
Note
*这个答案不会涉及到C中的rand()函数,因为它有很多问题。这里最严重的问题可能是,C标准没有明确地为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很小的时候才会很明显。或者更具体地说,当你modding的数字比RAND_MAX大的时候。
一个比循环更好的解决方案(循环效率非常低,甚至不应该被建议使用)是使用输出范围大得多的PRNG。梅森Twister算法的最大输出为4,294,967,295。这样做MersenneTwister::genrand_int32() % 10,将是均匀分布的,模偏效应将几乎消失。
我刚刚为冯·诺依曼无偏抛硬币法写了一段代码,理论上应该可以消除随机数生成过程中的任何偏差。更多信息请访问(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;
}
}
}
推荐文章
- 未定义对静态constexpr char的引用[]
- 什么是ORM,它是如何工作的,我应该如何使用它?
- 我能在服务器端应用程序(PHP、Ruby、Python等)上读取URL的哈希部分吗?
- 在c++中,restrict关键字是什么意思?
- c++中类似于java的instanceof
- 多少个参数是太多?
- include_directories和target_include_directories在CMake中的区别是什么?
- std::make_pair与std::pair的构造函数的目的是什么?
- 如何追加一个字符到std::字符串?
- 为什么要在c++中使用嵌套类?
- 如何处理11000行c++源文件?
- 使用g++编译多个.cpp和.h文件
- 如何在c++中追加文本到文本文件?
- 在c++中使用"super
- Mmap () vs.读取块