我看到很多人问过这个问题,但从未见过一个真正具体的答案。所以我打算在这里发布一个,希望能帮助人们理解为什么在使用随机数生成器时,比如c++中的rand(),会有“模偏置”。
当前回答
定义
模偏置是使用模算术将输出集缩减为输入集的子集时的固有偏置。一般来说,只要输入和输出集之间的映射不是均匀分布的,就会存在偏置,例如当输出集的大小不是输入集大小的除数时使用模算术。
这种偏差在计算中尤其难以避免,在计算中,数字被表示为比特串: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()是一个伪随机数生成器,它在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
RAND_MAX值为3(实际上它应该比这个值高得多,但偏差仍然存在),从这些计算中可以看出存在偏差:
1% 2 = 1 2% 2 = 0 3% 2 = 1 Random_between(1,3) % 2 =更可能是1
在本例中,当您想要0到1之间的随机数时,不应该使用% 2。你可以通过% 3得到一个0到2之间的随机数,因为在这种情况下:RAND_MAX是3的倍数。
另一种方法
有更简单的方法,但要加上其他答案,这是我的解,得到一个0到n - 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
正如公认的答案所示,“模偏置”的根源在于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;
}
}
}
不断随机选取是去除偏差的好方法。
更新
如果我们在能被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次迭代。
推荐文章
- decltype(auto)的一些用途是什么?
- Shared_ptr转换为数组:应该使用它吗?
- Printf与std::字符串?
- 禁用复制构造函数
- 如何在Ruby中生成a和b之间的随机数?
- 只接受特定类型的c++模板
- c#和Java中的泛型有什么不同?和模板在c++ ?
- c++ 11中的递归lambda函数
- 在c++中指针使用NULL或0(零)吗?
- 在c++中,如何将int值附加到字符串中?
- java.util.Random真的那么随机吗?我怎么能生成52!(阶乘)可能的序列?
- 从数组中随机选择一个元素
- 就性能而言,使用std::memcpy()还是std::copy()更好?
- 为什么布尔值是1字节而不是1位?
- 四舍五入到一个数字的最接近倍数