ChatGPT解决这个技术问题 Extra ChatGPT

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

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


M
Mark Amery

所以 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 参考

永远的困惑


另一种思考_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 )
我发布了一个附加答案,详细解释了问题并给出了示例代码解决方案。
在这种情况下,是否使用循环为侧信道攻击引入了空间?
k
klutt

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

更新

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


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

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

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

同样,最好的解决方案是在提供它的平台上使用 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。通过使用 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;
 }

注意,如果 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 主义)。
D
Dave R.

定义

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

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

潜在的解决方案

方法 0:

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

方法一:

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

方法二:

使用更多位:使用 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%

一个明确的改进,但在许多实际情况下还不够好。好消息是,添加更多位永远不会增加需要丢弃和重新滚动的机会。这不仅适用于骰子,而且适用于所有情况。

然而,正如演示的那样,添加一个额外的位可能不会改变任何事情。事实上,如果我们将滚动增加到 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;
    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() 断言,它始终会成功。
B
Ben Personick

马克的解决方案(公认的解决方案)几乎是完美的。

诠释 x;做 { x = rand(); } 而 (x >= (RAND_MAX - RAND_MAX % n)); x %= n;已于 2016 年 3 月 25 日 23:16 编辑 Mark Amery 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 不等于质数时,这是值得关注的)

例如:

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

因此,当您希望 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 条件之外进行比较操作的计算可能会导致更快的结果,因为无论需要多少次重试,这都是一次性计算。 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();
}

说马克的解决方案的问题在于他将 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。新年快乐! :-)
M
Mark Amery

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

一个适用于所有发电机。在极限情况下更容易看到。如果您的生成器的 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 很小,则可能会花费很多。显然,您可能会认为这些偏差对您的应用程序并不重要;我只是提供一种避免它们的方法。
R
Rivenfall

如果 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 不同的可能性,没有偏见。

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

从随机位编码数字

如果此数字 >= 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

P
Peter O.

模减少是使随机整数生成器避免永远运行的最坏情况的常用方法。

然而,当可能的整数范围未知时,通常没有办法在不引入偏差的情况下“修复”这种永远运行的最坏情况。不仅模约归约(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%,
IE。 2^(N-1)-1 是最大丢弃(其中 N 是 2 的幂,表示我们的集合 RAND_MAX --- i3 2^N 是随机函数可能返回的值集的计数,而 { 3} 是 2^N-1 ) 因此,为了便于查看,我们将每轮丢弃的最大机会称为 1/2。这能一直持续下去吗?是的,这是可能的,但是,会吗?这是极不可能的。
@BenPersonick:是的,正如您提到的,拒绝抽样可以在恒定的预期时间内实施。
C
Community

正如 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 相比较大时。

looping(效率极低,甚至不应该被建议)更好的解决方案是使用具有更大输出范围的 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
Y
Yavuz Koroglu

我刚刚为冯诺依曼的无偏硬币翻转方法编写了一个代码,理论上应该可以消除随机数生成过程中的任何偏差。可在 (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 的幂次方,需要解决模偏置问题。因此,该解决方案无法消除随机数生成过程中的任何偏差。
@里克嗯。冯诺依曼的方法在生成 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;}