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

Question

我已经看到这个问题了很多但从未见过真正的具体答案.所以我将在这里发布一个,希望能帮助人们理解为什么在使用随机数生成器时会出现"模数偏差",就像rand()在C++中一样.

Answer 1

所以rand()是选择0之间的自然数和伪随机数发生器RAND_MAX,它是在定义的常数cstdlib(见本文章有关的一般概述rand()).

如果你想在0到2之间生成一个随机数,会发生什么？为了便于解释,假设RAND_MAX为10,我决定通过调用生成0到2之间的随机数rand()%3.但是,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_MAX/n调用rand().

一种更有效的公式方法是采用一个可被整除的长度的大范围n,例如RAND_MAX - RAND_MAX % n,保持生成随机数,直到得到一个位于该范围内,然后取模数:

int x;

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

x %= n;

对于较小的值n,这将很少需要多次调用rand().

作品引用和进一步阅读:

• CPlusPlus参考

• 永远的困惑

Answer 2

继续选择随机是消除偏见的好方法.

更新

如果我们搜索可被整除的范围内的x,我们可以快速编写代码n.

// 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次迭代.

Answer 3

@ user1413793关于这个问题是正确的.我不打算进一步讨论,除了提出一点:是的,对于小值n和大值RAND_MAX,模偏差可能非常小.但是使用偏置诱导模式意味着每次计算随机数时都必须考虑偏差,并为不同情况选择不同的模式.如果你做出错误的选择,它引入的错误是微妙的,几乎不可能进行单元测试.与仅使用适当的工具(例如arc4random_uniform)相比,这是额外的工作,而不是更少的工作.做更多工作并获得更糟糕的解决方案是一项糟糕的工程,尤其是在大多数平台上每次都很容易做到这一点.

不幸的是,解决方案的实现都是错误的或效率低于应有的.(每个解决方案都有各种解释问题的注释,但没有解决任何解决方案来解决这些问题.)这可能会使偶然的答案者感到困惑,所以我在这里提供了一个已知良好的实现.

同样,最好的解决方案只是arc4random_uniform在提供它的平台上使用,或者为您的平台使用类似的远程解决方案(例如Random.nextInt在Java上).它会做正确的事情,无需代码成本.这几乎总是正确的召唤.

如果你没有arc4random_uniform,那么你可以使用opensource的强大功能来确切地了解它是如何在更广泛的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-upper_bound%upper_bound''.通过使用32位余数而不是64位余数,简化代码并使其在ILP32和LP64架构上都相同,并且在LP64架构上也略快一些.

Jorden Verwer在tech @ ok 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;
 }

Answer 4

定义

模数偏差是使用模运算将输出集减少到输入集子集的固有偏差.通常,只要输入和输出集之间的映射不是均匀分布就存在偏差,如在输出集的大小不是输入集大小的除数时使用模运算的情况.

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

使用更多位:而不是3位,使用4.这产生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%.

这引出了另外两个问题:

1. 如果我们添加足够的位,是否可以保证丢弃的概率会减少？
2. 一般情况下多少位就足够了？

一般解决方案

值得庆幸的是,第一个问题的答案是肯定的.6的问题是2 ^ x mod 6在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;
    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价值观,看看在大多数情况下实际发生了多少次重拍.持怀疑态度的人也可能希望将计算值保存到文件并验证分布是否正常.

Answer 5

使用模数有两种常见的抱怨.

• 一个对所有发电机都有效.在极限情况下更容易看到.如果您的生成器的RAND_MAX为2(不符合C标准)并且您只需要0或1作为值,则使用modulo将生成0两次(当生成器生成0和2时),因为它将生成1(当生成器生成1时).请注意,只要不删除值,就会发生这种情况,无论您使用从生成器值到所需的映射的映射,其中一个的频率将是另一个频率的两倍.

• 某种类型的发生器具有较少的有效位,其随机性较低,至少对于它们的一些参数而言,但遗憾的是这些参数具有其他有趣的特性(例如,能够使RAND_MAX小于2的幂).这个问题是众所周知的,并且很长一段时间库实现可能会避免这个问题(例如C标准中的示例rand()实现使用这种生成器,但是丢弃了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溢出)

BTW,C++ 11引入了减少和其他生成器的标准方法,而不是rand().

Answer 6

Mark的解决方案(已接受的解决方案)几乎是完美的.

int x;

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

x %= n;

于2016年3月25日23:16编辑

Mark Amery 39k21170211

但是,它有一个警告,在RAND_MAX(RM)小于N的倍数(其中N =可能的有效结果的数量)的任何情况下,丢弃1个有效的结果集.

即,当"丢弃的值的数量"(D)等于N时,它们实际上是有效集合(V),而不是无效集合(I).

使用Mark的解决方案,在以下情况下丢弃值:X => RM - RM%N

EG: 

Ran Max Value (RM) = 255
Valid Outcome (N) = 4

When X => 252, Discarded values for X are: 252, 253, 254, 255

So, if Random Value Selected (X) = {252, 253, 254, 255}

Number of discarded Values (I) = RM % N + 1 == N

 IE:

 I = RM % N + 1
 I = 255 % 4 + 1
 I = 3 + 1
 I = 4

   X => ( RM - RM % N )
 255 => (255 - 255 % 4) 
 255 => (255 - 3)
 255 => (252)

 Discard Returns $True

正如您在上面的示例中所看到的,当X的值(我们从初始函数得到的随机数)为252,253,254或255时,我们会丢弃它,即使这四个值包含一组有效的返回值.

IE:当值的计数Discarded(I)= 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%

由于Rerolls所需的百分比增加,N越接近RM,这可能是许多不同值的有效关注点,这取决于运行代码的系统的约束和所寻找的值.

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

 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.

广义版本1:

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语句放入最终的乘数.但也许你有一些代码可能有正当理由在n = 1时调用函数时返回1 ...

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

广义版本2:

int x;

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

    x %= n;
} else {
    x = rand();
}

当RM + 1是n的乘积时,这两种解决方案都解决了不必要地丢弃的有效结果的问题.

当您需要n等于RAND_MAX中包含的总可能值集时,第二个版本还涵盖了边缘情况.

两者中的修改方法是相同的,并且允许更一般地解决提供有效随机数和最小化丢弃值的需要.

重申:

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

 int x;

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

 x %= n;

扩展通用解决方案允许RAND_MAX + 1 = n的另一个场景:

int x;

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

    x %= n;
} else {
    x = rand();
}