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

我已经看到这个问题了很多但从未见过真正的具体答案 . 所以我将在这里发布一个,希望能帮助人们理解为什么在使用随机数生成器时确实存在"modulo bias",如C中的 rand() .

回答(9)

3 years ago

因此 rand() 是一个伪随机数生成器,它选择0到 RAND_MAX 之间的自然数,这是 cstdlib 中定义的常量(有关 rand() 的一般概述,请参阅此article) .

如果你想在0到2之间生成一个随机数,会发生什么?为了便于解释,假设 RAND_MAX 为10,我决定通过调用 rand()%3 生成0到2之间的随机数 . 但是, rand()%3 不会以相同的概率产生0到2之间的数字!

When rand() returns 0, 3, 6, or 9, rand()%3 == 0 . 因此,P(0)= 4/11

When rand() returns 1, 4, 7, or 10, rand()%3 == 1 . 因此,P(1)= 4/11

When rand() returns 2, 5, or 8, rand()%3 == 2 . 因此,P(2)= 3/11

这不会以相等的概率生成0到2之间的数字 . 当然对于小范围,这可能不是最大的问题,但是对于更大的范围,这可能会扭曲分布,从而偏向较小的数字 .

那么 rand()%n 何时以相同的概率返回从0到n-1的数字范围?当 RAND_MAX%n == n - 1 . 在这种情况下,与我们之前的假设 rand() 确实以相等的概率返回介于0和_3008613之间的数字,n的模数类也将是均等分布的 .

那么我们如何解决这个问题呢?粗略的方法是保持生成随机数,直到得到所需范围内的数字:

int x; 
do {
    x = rand();
} while (x >= n);

但是对于低值的 n 来说这是低效的,因为你只有 n/RAND_MAX 的机会获得你的范围内的值,所以你需要平均对 RAND_MAX/n 执行 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() .


作品引用和进一步阅读:


3 years ago

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

Update

如果我们搜索可被 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次迭代 .

3 years ago

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

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

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

如果你没有 arc4random_uniform ,那么你可以使用opensource的强大功能来确切地了解它是如何在更宽范围的RNG上实现的(在这种情况下 ar4random ,但类似的方法也可以在其他RNG之上工作) .

这是OpenBSD implementation

/*
 * 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'' . 简化代码并使其在ILP32和LP64架构上都相同,并且通过使用32位在LP64架构上也略快一些余数而不是64位余数 . 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;
 }

3 years ago

定义

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

在计算中特别难以避免这种偏差,其中数字表示为位串:0和1 . 找到真正随机的随机源也非常困难,但超出了本讨论的范围 . For the remainder of this answer, assume that there exists an unlimited source of truly random bits.

问题示例

让我们考虑使用这些随机位来模拟掷骰子(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 . This die is loaded.

潜在解决方案

方法0:

理论上,人们可以雇佣一支小军队整天掷骰子并将结果记录在数据库中,然后只使用一次结果,而不是依赖随机比特 . 这听起来和实际情况一样实际,并且很可能不会产生真正随机的结果(双关语) .

方法1:

而不是使用模数,一个天真但数学上正确的解决方案是丢弃产生 110111 的结果,并简单地再尝试3个新位 . 不幸的是,这意味着自己有一个 25% chance on each roll that a re-roll will be required, including each of the re-rolls . 除了最微不足道的用途之外,这显然是不切实际的 .

方法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%

一个明显的改进,但在许多实际情况下还不够好 . 好消息是, adding more bits will never increase the chances of needing to discard and re-roll . 这不仅适用于骰子,也适用于所有情况 .

正如所示 however, adding an 1 extra bit may not change anything. 事实上,如果我们将滚动增加到6位,概率仍为6.25% .

这引出了另外两个问题:

  • 如果我们添加足够的位,是否可以保证丢弃的概率会减少?

  • 一般情况下多少位就足够了?

一般解决方案

值得庆幸的是,第一个问题的答案是肯定的 . 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;
}

我鼓励玩 MODULUSROLLS 值,看看在大多数情况下实际发生了多少次重拍 . 持怀疑态度的人也可能希望将计算值保存到文件并验证分布是否正常 .

3 years ago

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

  • 一个对所有发电机都有效 . 在极限情况下更容易看到 . 如果您的生成器的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() .

3 years ago

Mark的解决方案(公认的解决方案)几乎是完美的 .

int x;

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: When the count of the values Discarded (I) = N (The number of valid outcomes) then a Valid set of return values will be discarded by the original function.

如果我们将值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,这可能是许多不同值的有效关注点,这取决于运行代码的系统的约束和所寻找的值 .

To negate this we can make a simple amendment As shown here:

int x;

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

 x %= n;

这提供了一个更通用的公式版本,它解释了使用模数来定义最大值的额外特性 .

Examples of using a small value for RAND_MAX which is a multiplicative of 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错误

Generalized Version 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中包含的总可能值集时,第二个版本还涵盖了边缘情况 .

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

重申:

The Basic General Solution which extends mark's example:

int x;

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

 x %= n;

The Extended General Solution which Allows one additional scenario of 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();
}

3 years ago

RAND_MAX 的值为 3 (实际上它应该远高于此值,但偏差仍然存在)从这些计算中可以看出存在偏差:

1 % 2 = 1 2 % 2 = 0 3 % 2 = 1 random_between(1, 3) % 2 = more likely a 1

在这种情况下, % 2 是您想要 01 之间的随机数时不应该做的事情 . 你可以通过 % 3 获得 02 之间的随机数,因为在这种情况下: RAND_MAX3 的倍数 .

Another method

有更简单的但是要添加到其他答案,这里是我的解决方案,以获得 0n - 1 之间的随机数,所以 n 不同的可能性,没有偏见 .

  • 编码可能性数量所需的位数(不是字节数)是您需要的随机数据的位数

  • 从随机位编码数字

  • 如果这样number是 >= 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

3 years ago

正如accepted answer所示,"modulo bias"的根源是 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个机会) .

所以它有偏见 . 较低的数字有更好的机会出来 .

But that only shows up so obviously when RAND_MAX is small . 或者更具体地说,当您修改的数字与 RAND_MAX 相比较大时 .

looping (这是非常低效且甚至不应该被建议)更好的解决方案是使用具有更大输出范围的PRNG . Mersenne Twister算法的最大输出为4,294,967,295 . 因此,对于所有意图和目的而言,将平均分配并且模数偏差效应将几乎消失 .

3 years ago

我刚刚编写了Von Neumann的无偏硬币翻转方法的代码,理论上应该消除随机数生成过程中的任何偏差 . 更多信息可以在(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;
        }
    }
}