为什么人们说在使用随机数生成器时会出现模数偏差?

329

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

11个回答

463

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_MAX/nrand()调用。

更有效的方法是使用一些长度可被n整除的大范围,比如RAND_MAX - RAND_MAX % n,持续生成随机数直到获得一个在该范围内的数,然后取模即可:

int x;

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

x %= n;

对于小的n值,这很少需要多次调用rand()函数。


参考文献和进一步阅读:



12
另一个理解RAND_MAX%n == n - 1的方法是(RAND_MAX + 1) % n == 0。在阅读代码时,我更容易将% something == 0理解为“可以被整除”,而不是其他计算方式。当然,如果你的C ++ stdlib将RAND_MAX设置为与INT_MAX相同的值,则(RAND_MAX + 1)肯定不起作用。因此,Mark的计算仍然是最安全的实现。 - Slipp D. Thompson
1
我可能有点挑剔,但如果目标是减少浪费的位数,我们可以针对RAND_MAX(RM)仅比N等分低1个的边缘条件稍作改进。在这种情况下,不需要通过执行X> =(RM-RM%N)来浪费任何位,这对于小值N的价值很小,但对于大值N的价值变得更大。如Slipp D. Thompson所述,有一个解决方案,只有当INT_MAX(IM)> RAND_MAX时才起作用,但当它们相等时会出现问题。然而,我们可以通过以下简单的方法修改计算X> =(RM-RM%N)来解决这个问题: - Ben Personick
2
X >= RM - (((RM % N) + 1) % N) - Ben Personick
我发布了一个额外的答案,详细解释了问题并提供了示例代码解决方案。 - Ben Personick
在这种情况下,使用循环是否会引入侧信道攻击的可能性? - Rodolfo Carvalho
关键弱点:当 n > RAND_MAX 时,do { ... } while (x >= (RAND_MAX - RAND_MAX % n)); 将会进入无限循环。 - chux - Reinstate Monica

37

随机选择是消除偏差的一种好方法。

更新

如果在范围内查找可被 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次。


2
呸 :-P 将其转换为双精度,然后乘以MAX_UPPER_LIMIT/RAND_MAX更加简洁,性能也更好。 - boycy
23
你错过了重点。如果 rand() 可返回的值的数量不是 n 的倍数,那么无论你做什么,都会出现“取模偏差”的问题,除非你丢弃其中一些值。user1413793 解释得很好(尽管该答案提出的解决方案确实有些麻烦)。 - TonyK
6
@TonyK,非常抱歉,我错过了重点。没有好好思考,认为偏差只会出现在使用显式模操作的方法中。感谢您纠正我 :-) - boycy
4
如果 RAND_MAX 等于整型最大值 *(在大多数系统上是这样)*,这个方法就行不通了。请参考我在评论区给 @user1413793 的第二条评论。 - BlueRaja - Danny Pflughoeft
1
@BlueRaja-DannyPflughoeft 在大多数系统上?我从未见过 RAND_MAX 不是 32767 的 libc 实现 - 包括 Microsoft 的 Visual libc、GLibC、BSD libc,甚至跨架构的实现。 - cat
显示剩余5条评论

23

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

不幸的是,所有解决方案的实现都不正确或不如应该的效率高。(每个解决方案都有各种说明问题的注释,但是没有一个解决方案已经被修复以解决它们。)这很可能会使普通的寻答者感到困惑,因此我在这里提供一个已知好的实现。

再次强调,最好的解决方案就是在提供这个功能的平台上使用arc4random_uniform,或者使用你的平台提供的类似范围解决方案(例如Java中的Random.nextInt)。这将在不增加代码成本的情况下正确地执行。这几乎总是正确的选择。

如果你没有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;
}

对于那些需要实现类似功能的人,值得注意这段代码的最新提交评论:

将arc4random_uniform()更改为计算2 ** 32%upper_bound, 以-upper_bound%upper_bound形式简化代码,并使其在ILP32和LP64架构上相同, 并使用32位余数而不是64位余数,在LP64架构上略微更快。

由Jorden Verwer指出技术@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
请注意,如果 arcfour_random() 函数在其实现中实际使用了真正的 RC4 算法,则输出将肯定存在一些偏差。希望您的库作者已经转而使用更好的 CSPRNG 实现相同的接口。我记得其中一个 BSD 现在实际上使用 ChaCha20 算法来实现 arcfour_random()。有关 RC4 输出偏差的更多信息,这使其在安全或其他关键应用程序(如视频扑克)方面变得无用:http://blog.cryptographyengineering.com/2013/03/attack-of-week-rc4-is-kind-of-broken-in.html?m=1 - rmalayter
5
在iOS和OS X上,arc4random从/dev/random读取数据,这是系统中最高质量的熵源。(名称中的“arc4”是历史遗留问题,并为了兼容性而保留。) - Rob Napier
2
@Rob_Napier 很好知道,但是 /dev/random 在过去的某些平台上也使用了 RC4(Linux 使用计数器模式下的 SHA-1)。不幸的是,我通过搜索找到的 man 页面表明,在提供 arc4random 的各种平台上仍在使用 RC4(尽管实际代码可能不同)。 - rmalayter
2
我有点困惑。-upper_bound % upper_bound == 0 不是成立的吗? - Jon McClung
2
@JonMcClung 如果 int 的宽度大于32位,则 -upper_bound%upper_bound 的确会为0。假设 u_int32_tuint32_t 的 BSD 特定类型,则应为 (u_int32_t)-upper_bound%upper_bound - Ian Abbott
显示剩余4条评论

20

定义

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

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

与其使用模数,一个朴素但数学上正确的解决方案是丢弃产生110111的结果,并只需用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;
    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的值,以查看在大多数情况下实际发生了多少次重投。一个怀疑的人也许还想将计算出的值保存到文件中,并验证分布是否正常。


1
我真的希望没有人盲目地复制了你的均匀随机实现。由于断言,randomPool = RAND_bytes(...)行将始终导致randomPool == 1。这总是导致丢弃和重新滚动。我认为你想在另一行上声明。因此,这导致RNG在每次迭代中返回1 - Qix - MONICA WAS MISTREATED
需要明确的是,根据OpenSSL对RAND_bytes()文档randomPool将始终评估为1,因为它将始终成功,这要归功于RAND_status()断言。 - Qix - MONICA WAS MISTREATED

11

马克的解决方案(被接受的解决方案)几乎完美。

int x;

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

x %= n;

然而,它有一个警告,在任何RAND_MAX (RM) 比 N(其中 N = 可能的有效结果数量) 少 1 的倍数少 1 的情况下,会丢弃一套有效的结果。

也就是说,当“值被丢弃的数量”(D)等于 N 时,它们实际上是一个有效的集合(V),而不是无效的集合(I)。

这是由于某些时候 Mark 忽视了 NRand_Max 之间的差异所引起的。

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

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

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

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

然而,N 作用于有效响应集合内的值的计数,因此即使在我们特定的情况下定义了它,Rand_Max 也将是一个比它包含的总数少 1 的值。

使用 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时,即使这四个值构成了一个有效的返回值集合,我们也会将其丢弃。

例如: 当被丢弃的值的数量 (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%

由于所需的重新投掷百分比随着N越接近RM而增加,因此这可能是一个有效的关注点,具体取决于运行代码的系统约束和寻找的值。

为了消除这个问题,我们可以进行简单的修改,如下所示:

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

这提供了一个更为通用的公式,考虑了使用模数来定义最大值时的额外特殊情况。

以N的乘法因子作为RAND_MAX的小值使用的示例。

Mark的原始版本:

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中的总可能值集合时。

这两个修改后的方法相同,并允许提供有效的随机数字并最小化丢弃的值的更通用的解决方案。

再次强调:

基本通用解决方案扩展了马克的示例:

// 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循环条件之外可能会导致更快的结果,因为这是一次性的计算,无论需要多少次重试。 个人经验因素!

// 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();
}

2
难道可以说Mark的解决方案存在问题在于他将RAND_MAX和n视为同一“度量单位”,而实际上它们代表两个不同的概念?其中n代表可能的“数量”,而RAND_MAX仅代表原始可能性的最大值,其中RAND_MAX + 1将是原始可能性的数量。我很惊讶他没有得出你的结论,因为他似乎已经承认n和RAND_MAX不是同一个东西,其等式为:RAND_MAX%n = n-1 - Danilo Souza Morães
@DaniloSouzaMorães 谢谢你,Danilo。你已经非常简洁地阐述了问题。我试图演示他正在做什么以及为什么和如何做,但我认为我从未能够清楚地陈述他在做什么错误,因为我过于沉迷于逻辑细节上的问题,没有清晰地说明问题所在。你介意我修改我的答案,使用你在这里写的一些内容作为我自己对问题的总结,以及需要在顶部解决的接受的解决方案在做什么? - Ben Personick
最后一次编辑(2020年)在我看来是错误的,@BenPersonick。y没有在n!= 0分支之外使用,并且由于零除法(... % n),在分支之外毫无意义。 - Palec
@palec y 可以避免在每次 rhencode 迭代时运行静态计算,从而减少 CPU 循环等待的负担。我通常在新年晚宴上,但这只是加速代码的示例。Y 必须在每次运行时计算一次,这会创建 6 个内存空间使用,但意味着它将成为一个缓存内存调用,可能在 CPU 缓存中进行比较,而不是实际的 CPU 计算。但是,CPU 比较也可能完全来自缓存,因此可能没有区别,或者哪个更快可能不同。YMMV - Ben Personick
@BenPersonick,我理解为什么需要y,即一些编译器不会将其提升出循环,需要手动提升。我只是认为y的定义应该在do-while循环之前,而不是更早。想想当n == 0时。新年快乐! :-) - Palec
@palec 啊,好的,我明白了。把它放在第一个比较里是公平的,这个边缘情况不需要计算。但这只是一个优化的例子。 - Ben Personick

10

使用取模运算有两个常见问题。

  • 其中一个问题适用于所有生成器。在一个极限情况下更容易看出来。如果您的生成器的RAND_MAX为2(这不符合C标准),而您只想要值为0或1,则使用取模将使得0的出现次数比1多一倍(当生成器生成0和2时会生成0两次)。请注意,只要您不放弃任何值,不管您从生成器值到所需值的映射是什么,其中一个值将比另一个值出现两次频率都高。

  • 某些类型的生成器其低位更不随机,至少对于其中一些参数是这样,但不幸的是这些参数具有其他有趣的特性(例如能够使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内的偏差,但我无法精确地将代码映射到这个想法上。 - zinking
1
天真的版本应该是(RAND_MAX+1)/(n+1),因为有RAND_MAX+1个值要分成n+1个桶。为了避免计算RAND_MAX+1时出现溢出,可以将其转换为1+(RAND_MAX-n)/(n+1)。为了避免计算n+1时出现溢出,首先检查n==RAND_MAX的情况。 - AProgrammer
+加号,除法似乎比重新生成数字的成本更高。 - zinking
4
取模和除法的成本相同。有些指令集甚至只提供一个指令来同时执行这两个操作。重新生成数字的成本将取决于n和RAND_MAX。如果n相对于RAND_MAX很小,那么重建数字的成本可能很高。显然,您可以决定偏差对于您的应用程序不重要;我只是提供了一种避免它们的方法。 - AProgrammer

1

RAND_MAX 的值为 3 时(实际上应该比这个高得多,但偏差仍然存在),从这些计算中可以看出有偏差:

1 % 2 = 1 2 % 2 = 0 3 % 2 = 1 random_between(1, 3) % 2 = 更可能是 1

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

另一种方法

有一种更简单的方法来避免偏差,以补充其他答案,这是我获取介于 0n - 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

1
模数约简是一种常见的方法,可以使随机整数生成器避免永远运行的最坏情况。
然而,当可能整数范围未知时,通常没有办法“修复”这个最坏情况,而不引入偏差。这不仅适用于模数约简(rand() % n,在接受的答案中讨论),还适用于Daniel Lemire的“乘法和移位”约简,或者如果您在一定次数的迭代后停止拒绝结果。(要明确的是,这并不意味着没有办法解决伪随机生成器存在的偏差问题。例如,即使模数和其他约简通常是有偏差的,但如果可能整数范围是2的幂次方,并且随机生成器产生无偏随机位或块,则它们将不会有偏差问题。)
以下是我的回答,讨论了在假设有一个“真正”的随机生成器可以产生无偏和独立的随机位的情况下,运行时间和偏差之间的关系。这个答案甚至不涉及C语言中的rand()函数,因为它存在很多问题。也许最严重的问题在于C标准并没有明确指定rand()返回的数字分布,甚至没有均匀分布。


除了处理一些与 OP 的问题无关的转移范围之外,(包括本答案在内的所有答案中的 IMP 似乎只会使人们对正在完成的事情产生困惑)。话虽如此,这段代码似乎只是在解决模数偏差本身的相同根本原因,即 RAND_MAX 总是2的幂,因此当 SET 不是2的幂时,您必须丢弃掉落到错误集合中的值。这在我的答案和被接受的答案中得到了解决,但您似乎认为它没有。 - Ben Personick
@BenPersonick:我的回答说没有办法在不引入偏差的情况下“修复”永远运行的最坏情况,并不是说伪随机生成器存在偏差问题就没有办法解决。当整数范围未知时,偏差问题通常只能通过拒绝抽样来解决,例如你的答案或这个答案中提供的技术,并且拒绝抽样具有无上限的最坏情况运行时间。我将澄清这个答案。 - Peter O.
啊,我明白了,你的意思是要提出我们所有代码都存在的隐含问题,这一点对我来说并不十分清楚。尽管如此,在实际情况下,除非基础的伪随机数生成具有显著偏差,否则它永远运行的可能性相当小。每一轮都有可能被丢弃,从未真正达到50%。 - Ben Personick
2^(N-1)-1 是最大的丢弃值(其中 N 是表示结果集 RAND_MAX 的 2 的幂次方,即 2^N 是随机函数可能返回的值的集合计数,而 RAND_MAX2^N-1)。因此,为了便于审查,我们将每轮的最大丢弃概率称为 1/2。这种情况会一直持续下去吗?是的,这是可能的,但是会吗?这是极不可能的。 - Ben Personick
@BenPersonick:是的,拒绝抽样可以像你提到的那样在常数_期望_时间内实现。 - Peter O.
因此,对于我们使用1/2的机会(也称为50%的概率)的参数,我们有(惊喜)另一个2的幂。您在丢弃集合X次的概率是1 /(2 ^ X)(例如,1个丢弃1 /(2 ^ 1) 50% - 4个丢弃是1 / 16 6.25% - 10个丢弃是1 / 1024 0.00098% - 连续被丢弃100次的概率是1 /(2 ^ 100),即7.88860905E-31%,实际上您可以将每次投掷的百分比相加,以查看“93.75%的尝试不会超过4次丢弃”,而“96.87%”不会超过5次丢弃。“98.43%”在6次尝试中,没有超过10个丢弃的机会为99.88%。这是最坏情况。 - Ben Personick

-1
正如被接受的答案所指出的那样,“模数偏差”源于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(10 中 4 的机会),只有 3 个输出是 1 和 2(每个 10 中 3 的机会)。

所以这是有偏差的。较低的数字更有可能出现。

但仅当 RAND_MAX 很小时才会很明显。或更具体地说,当您要进行取模运算的数字与 RAND_MAX 相比较大时。

循环 更好的解决方案(极其低效,甚至不应建议使用)是使用具有更大输出范围的 PRNG。Mersenne Twister 算法的最大输出为 4,294,967,295。因此,对于所有目的,执行 MersenneTwister::genrand_int32() % 10 将是等分布的,而模数偏差效应几乎会消失。


3
你的方法更高效,如果RAND_MAX明显比你取模的数字大得多,则可能是正确的。然而,即使如此,你的方法仍然会存在偏差。尽管这些都是伪随机数发生器,本身也是一个不同的话题,但是如果你假设完全的随机数生成器,你的方法仍然会对较小的值产生偏差。 - user1413793
因为最高值是奇数,MT::genrand_int32()%2 会在50 + 2.3e-8%的时间内选取0,在50 - 2.3e-8%的时间内选取1。除非你在构建赌场的随机数发生器(对于这种情况,你可能需要使用更大范围的RGN),否则任何用户都不会注意到额外2.3e-8%的时间差。在这里讨论的是微不足道的数字。 - bobobobo
8
循环是最好的解决方案。它并非“极其低效”; 在最坏的平均情况下,所需的迭代次数不到两倍。使用高RAND_MAX值将减少模数偏差,但无法消除它。只有循环才能实现。 - Jared Nielsen
6
如果RAND_MAX比你对其取模的数大得足够多,那么需要重新生成随机数的次数将非常少,并且不会影响效率。我建议保留循环,只要你在测试中使用的是最大的n的倍数而不仅仅是像被采纳的答案建议的那样使用n - Mark Ransom

-1

我很喜欢使用各种软件生成随机双精度数。如果我使用((double)rand() / RAND_MAX),我发现范围更加“随机”。所以我猜想,如果你将这个乘以你的数值范围,你可以得到一个更少偏倚的随机数?

例如,((double)rand() / RAND_MAX) * 3。

我读到一个关于从2中得到一个随机数的答案。isodd(rand())?


这并没有真正回答问题。如果你有其他问题,可以点击提问来提出。如果你想在这个问题有新的回答时收到通知,你可以关注此问题。一旦你拥有足够的声望,你还可以添加悬赏以吸引更多关注。- 来自评论 - YesThatIsMyName

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接