我看到很多人问过这个问题,但从未见过一个真正具体的答案。所以我打算在这里发布一个,希望能帮助人们理解为什么在使用随机数生成器时,比如c++中的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次迭代。
其他回答
@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;
}
我刚刚为冯·诺依曼无偏抛硬币法写了一段代码,理论上应该可以消除随机数生成过程中的任何偏差。更多信息请访问(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:
不使用模量,一个简单但在数学上正确的解决方案是丢弃产生110和111的结果,并简单地重新尝试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%
确实有了改进,但在许多实际情况下还不够好。好消息是,添加更多比特永远不会增加需要丢弃和重新滚动的几率。这不仅适用于骰子,而且适用于所有情况。
然而,如前所述,增加1个额外的位可能不会改变任何东西。事实上,如果我们将点数增加到6位,概率仍然是6.25%。
这就引出了另外两个问题:
如果我们添加足够多的比特,是否能保证丢弃的概率会降低? 一般情况下多少位才够呢?
通解
幸运的是,第一个问题的答案是肯定的。6的问题在于,2^x mod 6在2和4之间翻转而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;
}
我鼓励使用MODULUS和ROLLS值来查看在大多数情况下实际发生了多少次重新滚动。持怀疑态度的人也可能希望将计算值保存到文件中,并验证分布是否正常。
正如公认的答案所示,“模偏置”的根源在于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,将是均匀分布的,模偏效应将几乎消失。
对于模的使用,有两种常见的抱怨。
one is valid for all generators. It is easier to see in a limit case. If your generator has a RAND_MAX which is 2 (that isn't compliant with the C standard) and you want only 0 or 1 as value, using modulo will generate 0 twice as often (when the generator generates 0 and 2) as it will generate 1 (when the generator generates 1). Note that this is true as soon as you don't drop values, whatever the mapping you are using from the generator values to the wanted one, one will occurs twice as often as the other. some kind of generator have their less significant bits less random than the other, at least for some of their parameters, but sadly those parameter have other interesting characteristic (such has being able to have RAND_MAX one less than a power of 2). The problem is well known and for a long time library implementation probably avoid the problem (for instance the sample rand() implementation in the C standard use this kind of generator, but drop the 16 less significant bits), but some like to complain about that and you may have bad luck
使用类似于
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()以外的其他生成器。