为啥人们说使用随机数生成器时存在模偏差?

Posted

技术标签:

【中文标题】为啥人们说使用随机数生成器时存在模偏差?【英文标题】:Why do people say there is modulo bias when using a random number generator?为什么人们说使用随机数生成器时存在模偏差? 【发布时间】:2012-06-14 15:29:35 【问题描述】:

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

【问题讨论】:

【参考方案1】:

取模是一种常见的方法,可以使随机整数生成器避免永远运行的最坏情况。

但是,当可能的整数范围未知时,通常没有办法“修复”这种最坏的情况,即在不引入偏差的情况下永远运行。不仅是模减少(rand() % n,在接受的答案中讨论)会以这种方式引入偏差,而且 Daniel Lemire 的“乘法和移位”减少,或者如果你在设定的次数后停止拒绝结果迭代。 (要清楚,这并不意味着没有办法解决伪随机生成器中存在的偏差问题。例如,即使模数和其他归约通常是有偏差的,如果可能的范围内,它们不会有偏差问题integers 是 2 的幂 并且 如果随机生成器产生无偏的随机位或它们的块。)

此答案的其余部分将显示随机生成器中运行时间和偏差之间的关系。从这里开始,我们将假设我们有一个“真正的”随机生成器,它可以产生无偏且独立的随​​机位。*

1976 年,DE Knuth 和 AC Yao 表明,任何仅使用随机位以给定概率生成随机整数的算法都可以表示为二叉树,其中随机位指示遍历树和每个叶子的方式(端点)对应于一个结果。在这种情况下,我们正在处理在 [0, n) 中生成随机整数的算法,其中每个整数的选择概率为 1/n。如果对于所有结果,树中出现相同数量的叶子,则该算法是无偏的。但是,如果 1/n 具有非终止二进制展开式(如果 n 不是 2 的幂,就会出现这种情况),那么只有在以下情况下,该算法才是无偏的——

二叉树具有“无限”深度,或者 二叉树末尾包含“拒绝”叶子,

在任何一种情况下,算法都不会在恒定时间内运行,并且在最坏的情况下会永远运行。 (另一方面,当n 是 2 的幂时,最优二叉树将具有有限深度且没有拒绝节点。)

二叉树的概念还表明,任何“修复”这种最坏情况时间复杂度的方法通常都会导致偏差。 (同样,这并不意味着没有办法解决伪随机生成器中存在的偏差问题。)例如,模约简相当于一棵二叉树,其中拒绝叶被标记的结果替换——但因为有更多可能结果比拒绝叶子,只有一些结果可以代替拒绝叶子,从而引入偏见。如果您在一定次数的迭代后停止拒绝,则会产生相同类型的二叉树 - 以及相同类型的偏差。 (但是,根据应用程序,这种偏差可能可以忽略不计。随机整数生成也有安全方面的问题,这太复杂了,无法在此答案中讨论。)

为了说明,以下 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
    
 

注意

* 这个答案不会涉及 C 中的 rand() 函数,因为它是 has many issues。这里最严重的可能是 C 标准没有明确指定 rand() 返回的数字的特定分布,甚至没有统一分布。

【讨论】:

除了处理与 OP 的问题无关的偏移范围之外,(包括这个问题在内的所有答案中的 IMP 似乎只会使正在完成的事情变得混乱) .也就是说,这段代码似乎只是解决了模数偏差本身的相同根本原因,即 RAND_MAX 将始终是 2 的幂,因此当 SET 不是 2 的幂时,您必须丢弃掉入坏集。这在我和接受的答案中得到了解决,但您似乎认为不是...... @BenPersonick:我的回答是,如果不引入偏见,就无法“解决”最坏的情况永远运行,并不是说没有办法解决偏见问题存在伪随机生成器。当整数的范围未知时,通常只能通过拒绝抽样来解决偏差问题,例如您的答案或这个答案中给出的技术,并且拒绝抽样具有无限的最坏情况运行时间。我会澄清这个答案。 啊,我明白了,我并不清楚你的意思是提出我们所有代码都存在的隐含问题。虽然,实际上,除非潜在的伪随机数生成具有显着偏差,否则它永远运行的机会非常小。每一轮都有可能被丢弃,但实际上从未达到 50%, 即。 2^(N-1)-1 是最大丢弃(其中 N 是 2 的幂,代表我们的集合 RAND_MAX --- i3 2^N 是随机函数可能返回的值集的计数,而 @ 987654332@ 是2^N-1 ) 因此,为了便于审查,我们将每轮的最大丢弃机会称为 1/2。这能永远持续下去吗?是的,这是可能的,但是,会吗?这是极不可能的。 @BenPersonick:是的,正如您提到的,拒绝抽样可以在恒定的预期时间内实施。【参考方案2】:

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

int x;

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

x %= n;

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

马克艾默里 39k21170211

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

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

造成这种情况的原因是,Mark 在某些时候忽略了 NRand_Max 之间的区别。

N 是一个集合,其有效成员仅由正整数组成,因为它包含有效响应的计数。 (例如:设置N = 1, 2, 3, ... n

Rand_max 然而,它是一个集合(根据我们的目的定义)包含任意数量的非负整数。

在它最通用的形式中,这里定义为 Rand Max 是所有有效结果的集合,理论上可以包括负数或非数字值。

因此,Rand_Max 最好定义为“可能的响应”集合。

但是,N 会针对有效响应集中的值计数进行操作,因此即使在我们的特定案例中定义,Rand_Max 的值也会比它包含的总数小一。

使用 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:当丢弃值的计数 (I) = N(有效结果的数量)时,原始函数将丢弃一组有效的返回值。

如果我们将值N和RM之间的差异描述为D,即:

D = (RM - N)

然后,随着 D 的值变小,由于这种方法导致的不需要重投的百分比在每个自然乘法处都会增加。 (当 RAND_MAX 不等于质数时,这是值得关注的)

EG:

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%

由于 N 越接近 RM,所需的 Rerolls 百分比就会增加,因此根据运行代码的系统的约束和正在寻找的值,这可能会在许多不同的值上引起关注。

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

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

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

对 RAND_MAX 使用较小值的示例,该值是 N 的乘积。

标记原版:

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...

所以最好使用 0,当您希望 n = RAND_MAX+1 时,这通常会提供 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 中包含的所有可能值集的边缘情况。

两者中的修改方法是相同的,并且允许更通用的解决方案来满足提供有效随机数和最小化丢弃值的需求。

重申:

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

// Assumes:
//  RAND_MAX is a globally defined constant, returned from the environment.
//  int n; // User input, or externally defined, number of valid choices.

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

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

// Assumes:
//  RAND_MAX is a globally defined constant, returned from the environment.
//  int n; // User input, or externally defined, number of valid choices.

int x;

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

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

在某些语言(尤其是解释性语言)中,在 while 条件之外进行比较操作的计算可能会导致更快的结果,因为无论需要多少次重试,这都是一次性计算。 YMMV!

// Assumes:
//  RAND_MAX is a globally defined constant, returned from the environment.
//  int n; // User input, or externally defined, number of valid choices.

int x; // Resulting random number
int y; // One-time calculation of the compare value for x

y = RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) 

if n != 0 
    do 
        x = rand();
     while (x > y);

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

【讨论】:

难道不能说 Mark 的解决方案的问题在于他将 RAND_MAX 和 n 视为相同的“度量单位”,而实际上它们意味着两个不同的东西?虽然 n 代表最终的“可能性数”,但 RAND_MAX 仅代表原始可能性的最大值,其中 RAND_MAX + 1 将是原始可能性数。我很惊讶他没有得出你的结论,因为他似乎已经承认 n 和 RAND_MAX 与等式不同:RAND_MAX%n = n - 1 @DaniloSouzaMorães 谢谢 Danilo,你把这件事说得很简洁。我去展示他在做什么以及为什么以及如何做,但不要认为我能够雄辩地说明他做错了什么,因为我非常专注于关于如何和如何做的逻辑细节为什么存在问题,我没有清楚地说明问题所在。您是否介意我修改我的答案以使用您在此处写的一些内容作为我自己的总结,以解决接受的解决方案在做什么和在哪里做什么以及需要在顶部附近解决什么问题? 最后一次编辑(2020 年)是 IMO 错误,@BenPersonick。 y 不在n != 0 分支之外使用,由于被零除(... % n),在分支之外没有任何意义。 @palec y 不再需要在每次运行 rhencode 时多次运行静态计算,因为其他解决方案要求它在每次迭代等待 CPU 周期时运行。我在新年每顿晚餐,但这是如何加速代码的一个例子。每次运行时必须始终计算一次 Y,创建 6 次使用内存空间,但这意味着每次比较与实际 CPU 计算相比,这将是一个可能在 CPU 缓存上的内存调用,但 CPU 比较也可能完全从缓存中完成也一样,所以,可能没有不同,或者哪个更有趣可能不同。 YMMV @BenPersonick,我明白为什么需要y,即某些编译器不会将其提升出循环并且需要手动提升。我只是认为y 的定义应该发生在do-while 循环之前,而不是更早。想想什么时候n == 0。新年快乐! :-)【参考方案3】:

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

不幸的是,解决方案的实施都是不正确的,或者效率低于应有的水平。 (每个解决方案都有不同的 cmets 来解释问题,但没有一个解决方案被修复来解决这些问题。)这可能会使随便的答案寻求者感到困惑,所以我在这里提供了一个已知良好的实现。

同样,最好的解决方案是在提供 arc4random_uniform 的平台上使用它,或者为您的平台使用类似的范围解决方案(例如 Java 上的 Random.nextInt)。它会做正确的事情,而不会给您带来任何代码成本。这几乎总是正确的选择。

如果你没有arc4random_uniform,那么你可以使用开源的力量来看看它是如何在更广泛的 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 架构上相同,并且在 ILP32 和 LP64 架构上也稍快 LP64 架构使用 32 位余数而不是 64 位 余数。

Jorden Verwer 在 tech@ 上指出 好的没有来自 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;
 

【讨论】:

注意,如果arcfour_random()在其实现中实际使用了真正的RC4算法,那么输出肯定会有一些偏差。希望您的库作者已经切换到在同一界面后面使用更好的 CSPRNG。我记得现在有一个 BSD 实际上使用 ChaCha20 算法来实现arcfour_random()。更多关于 RC4 输出偏差使其无法用于安全或其他关键应用(如视频扑克)的信息:blog.cryptographyengineering.com/2013/03/… @rmalayter 在 ios 和 OS X 上,arc4random 从 /dev/random 读取,这是系统中最高质量的熵。 (名称中的“arc4”是历史性的,为了兼容性而保留。) @Rob_Napier 很高兴知道,但 /dev/random 过去在某些平台上也使用过 RC4(Linux 在计数器模式下使用 SHA-1)。不幸的是,我通过搜索找到的手册页表明 RC4 仍在提供arc4random 的各种平台上使用(尽管实际代码可能不同)。 我很困惑。不是-upper_bound % upper_bound == 0吗?? @JonMcClung -upper_bound % upper_bound 如果int 大于 32 位,则确实为 0。它应该是(u_int32_t)-upper_bound % upper_bound)(假设u_int32_tuint32_t 的BSD 主义)。【参考方案4】:

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

更新

如果我们在可被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 次迭代。

【讨论】:

Yuck :-P 转换为双精度,然后乘以 MAX_UPPER_LIMIT/RAND_MAX 更简洁,性能更好。 @boycy:你错过了重点。如果rand() 可以返回的值的数量不是n 的倍数,那么无论您做什么,都将不可避免地得到“模偏差”,除非您丢弃其中一些值。 user1413793 很好地解释了这一点(尽管该答案中提出的解决方案确实很糟糕)。 @TonyK 抱歉,我确实错过了重点。考虑得不够努力,并认为偏差仅适用于使用显式模运算的方法。感谢您修复我:-) 如果RAND_MAX == INT_MAX (就像在大多数系统上一样),这将不起作用。请参阅上面我对@user1413793 的第二条评论。 @BlueRaja-DannyPflughoeft 在大多数系统上?我从未见过 RAND_MAX 不是 32767 的 libc 实现——微软的 Visual libc、GLibC、BSD libc,甚至跨架构【参考方案5】:

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

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

在这种情况下,当您想要一个介于 01 之间的随机数时,您不应该使用 % 2。你可以通过% 3 得到一个介于02 之间的随机数,因为在这种情况下:RAND_MAX3 的倍数。

另一种方法

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

编码可能性数量所需的位数(不是字节数)是您需要的随机数据位数 从随机位编码数字 如果这个数字是&gt;= 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.Ds-s-random 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

【讨论】:

【参考方案6】:

使用模数有两个常见的抱怨。

一个对所有生成器都有效。在极限情况下更容易看到。如果您的生成器的 RAND_MAX 为 2(不符合 C 标准)并且您只需要 0 或 1 作为值,则使用模数生成 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 溢出)

顺便说一句,C++11 为归约和其他生成器引入了标准方法,而不是 rand()。

【讨论】:

n == RAND_MAX ? 1 : (RAND_MAX-1)/(n+1): 我理解这里的想法是先将 RAND_MAX 分成相等的页面大小 N,然后返回 N 内的偏差,但我无法将代码精确映射到这个。 天真的版本应该是 (RAND_MAX+1)/(n+1) 因为有 RAND_MAX+1 值可以划分为 n+1 个桶。如果在计算 RAND_MAX+1 时为了避免溢出,可以将其转化为 1+(RAND_MAX-n)/(n+1)。为了避免计算n+1时溢出,先检查n==RAND_MAX的情况。 +plus,与重新生成数字相比,进行除法的成本似乎更高。 取模和除的成本相同。一些 ISA 甚至只提供一条指令,它总是提供两种指令。重新生成数字的成本将取决于 n 和 RAND_MAX。如果 n 相对于 RAND_MAX 很小,则可能会花费很多。显然,您可能会认为这些偏差对您的应用程序并不重要;我只是提供一种避免它们的方法。【参考方案7】:

所以rand() 是一个伪随机数生成器,它选择一个介于 0 和 RAND_MAX 之间的自然数,这是在 cstdlib 中定义的常数(请参阅此 article 以了解有关 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 Reference

Eternally Confuzzled


【讨论】:

另一种思考方式_RAND_MAX%n == n - 1_ 是(RAND_MAX + 1) % n == 0。阅读代码时,我倾向于将% something == 0 理解为比其他计算方法更容易“整除”。 当然,如果你的 C++ 标准库的 RAND_MAXINT_MAX 的值相同,那么 (RAND_MAX + 1) 肯定行不通;所以 Mark 的计算仍然是最安全的实现。 我可能在吹毛求疵,但如果目标是减少浪费的位,我们可以稍微改进 RAND_MAX (RM) 仅比被 N 整除 1 的边缘条件。在这种情况下,不需要通过执行 X >= (RM - RM % N)) 来浪费任何位,这对于 N 的小值来说价值不大,但对于 N 的大值来说变得更大。正如 Slipp D. Thompson 所提到的,有一种解决方案仅在 INT_MAX (IM) > RAND_MAX 时有效,但在它们相等时会中断。但是,有一个简单的解决方案,我们可以将计算 X >= (RM - RM % N) 修改如下: X >= RM - ( ( ( RM % N ) + 1 ) % N ) 我发布了一个附加答案,详细解释了问题并给出了示例代码解决方案。 在这种情况下使用循环是否会为侧信道攻击引入空间?【参考方案8】:

定义

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

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

与其依赖随机位,理论上可以雇佣一支小军队整天掷骰子并将结果记录在数据库中,然后每个结果只使用一次。这与听起来一样实用,而且很可能无论如何都不会产生真正的随机结果(双关语)。

方法一:

不使用模数,一个简单但数学上正确的解决方案是丢弃产生110111 的结果,然后简单地用3 个新位重试。不幸的是,这意味着每次掷骰有 25% 的机会需要重新掷骰,包括每次重新掷骰。除了最微不足道的用途外,这显然对所有用途都是不切实际的。

方法二:

使用更多位:使用 4 位而不是 3 位。这会产生 16 种可能的结果。当然,在结果大于 5 的任何时候重新滚动会使情况变得更糟(10/16 = 62.5%),因此仅凭这一点是无济于事的。

请注意 2 * 6 = 12

一开始听起来不错,但让我们检查一下数学:

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%。

这引出了 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;

我鼓励使用MODULUSROLLS 值来查看在大多数情况下实际发生了多少重新滚动。持怀疑态度的人也可能希望将计算值保存到文件中并验证分布是否正常。

【讨论】:

我真的希望没有人盲目复制你的统一随机实现。由于断言,randomPool = RAND_bytes(...) 行将始终导致 randomPool == 1。这总是会导致弃牌和重掷。我认为您想在单独的行上声明。因此,这导致 RNG 每次迭代都返回 1 需要明确的是,根据 OpenSSL documentation for RAND_bytes()randomPool 将始终评估为 1,因为由于 RAND_status() 断言,它始终会成功。【参考方案9】:

我刚刚为冯诺依曼的无偏硬币翻转方法编写了一个代码,理论上应该可以消除随机数生成过程中的任何偏差。更多信息请访问 (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 到 n 的均匀分布,其中 n 不小于 2 的幂次方,需要解决模偏置问题。因此,该解决方案无法消除随机数生成过程中的任何偏差。 @Rick 嗯。冯诺依曼方法在生成 1 到 100 之间的随机数时消除模偏差的逻辑扩展是: A) 调用 rand() % 100 100 次。 B) 如果所有结果都不同,则取第一个。 C) 否则,转到 A。这将起作用,但预期迭代次数约为 10^42,您必须非常耐心。并且不朽。 @MarkAmery 确实应该可以。尽管没有正确实现,但请查看此算法。第一个 else 应该是:else if(prev==2) prev= x1; else if(prev!=x1) return prev; prev=2;【参考方案10】:

正如accepted answer 所指出的,“模偏差”源于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 很小的时候才会如此明显。或者更具体地说,当您修改的数字与RAND_MAX 相比较大时。

循环(效率极低,甚至不应该推荐)更好的解决方案是使用输出范围更大的 PRNG。 Mersenne Twister 算法的最大输出为 4,294,967,295。因此,出于所有意图和目的这样做 MersenneTwister::genrand_int32() % 10 将平均分配,并且模偏差效应将几乎消失。

【讨论】:

你的效率更高,如果 RAND_MAX 比你修改的数字大得多,这可能是真的,但是你的仍然会有偏差。当然,无论如何这些都是伪随机数生成器,这本身就是一个不同的主题,但如果你假设一个完全随机数生成器,你的方式仍然会偏向较低的值。 因为最大值是奇数,MT::genrand_int32()%2 选择 0 (50 + 2.3e-8)% 的时间和 1 (50 - 2.3e-8)% 的时间。除非您正在构建赌场的 RGN(您可能会使用更大范围的 RGN),否则任何用户都不会注意到额外的 2.3e-8% 的时间。您在这里谈论的数字太小而无关紧要。 循环是最好的解决方案。它不是“非常低效”;在最坏的平均情况下需要不到两倍的迭代。使用较高的RAND_MAX 值将减少模偏差,但不会消除它。循环将。 如果RAND_MAX 比你修改的数字足够大,你需要重新生成随机数的次数非常少,不会影响效率。我说保持循环,只要您测试的是n 的最大倍数,而不是接受答案所建议的n

以上是关于为啥人们说使用随机数生成器时存在模偏差?的主要内容,如果未能解决你的问题,请参考以下文章

随机数与取模结果对比

为啥我使用PyTorch来生成随机数总会失败?

为啥 MongoDB Java 驱动程序在条件中使用随机数生成器?

为啥 MongoDB Java 驱动程序在条件中使用随机数生成器?

C++ 实现随机数生成(WindowsLinux)

为啥我的随机数生成器在 C# 中看起来不是随机的?