16

在这个 StackOverflow 问题中：

从范围生成随机整数

接受的答案建议使用以下公式在给定min和之间生成一个随机整数max，其中min和max包含在范围内：

output = min + (rand() % (int)(max - min + 1))

但它也说

这仍然略微偏向于较低的数字......也可以扩展它以消除偏差。

但它没有解释为什么它偏向于较低的数字或如何消除这种偏见。所以，问题是：这是在（有符号）范围内生成随机整数的最佳方法，而不依赖于任何花哨的东西，只是rand()函数，如果它是最优的，如何消除偏差？

编辑：

我刚刚while针对浮点外推法测试了@Joey 建议的 -loop 算法：

static const double s_invRandMax = 1.0/((double)RAND_MAX + 1.0);
return min + (int)(((double)(max + 1 - min))*rand()*s_invRandMax);

查看有多少均匀的“球”“落入”并分布在多个“桶”中，一个测试用于浮点外推，另一个测试用于while循环算法。但是结果会根据“球”（和“桶”）的数量而有所不同，所以我无法轻易选出获胜者。工作代码可以在这个 Ideone 页面找到。例如，对于 10 个桶和 100 个球，对于浮点外推法，与while-loop 算法（分别为 0.04 和 0.05）相比，桶之间与理想概率的最大偏差较小（分别为 0.04 和 0.05），但是对于 1000 个球，while-loop 算法较小（0.024 和 0.011），并且对于 10000 个球，浮点外推再次做得更好（0.0034 和 0.0053），依此类推，没有太大的一致性。考虑到没有一种算法始终如一地产生比其他算法更好的均匀分布的可能性，这让我倾向于浮点外推，因为它似乎比while-loop 算法执行得更快。那么选择浮点外推算法是否可以，或者我的测试/结论不完全正确？

4

7 回答 7

14

问题是你正在做一个模运算。RAND_MAX如果可以被您的模数整除，这将没有问题，但通常情况并非如此。作为一个非常人为的示例，假设RAND_MAX为 11，您的模数为 3。您将获得以下可能的随机数和以下结果余数：

0 1 2 3 4 5 6 7 8 9 10
0 1 2 0 1 2 0 1 2 0 1

如您所见，0 和 1 的概率略高于 2。

解决此问题的一种选择是拒绝抽样：通过禁止上面的数字 9 和 10，您可以使结果分布再次均匀。棘手的部分是弄清楚如何有效地做到这一点。在 Java 的方法中可以找到一个非常好的示例（我花了两天时间才理解它为什么起作用）。java.util.Random.nextInt(int)

Java的算法有点棘手的原因是它们避免了像乘法和除法这样的慢操作来进行检查。如果你不太在乎，你也可以用幼稚的方式来做：

int n = (int)(max - min + 1);
int remainder = RAND_MAX % n;
int x, output;
do {
  x = rand();
  output = x % n;
} while (x >= RAND_MAX - remainder);
return min + output;

编辑：更正了上述代码中的栅栏错误，现在它可以正常工作了。我还创建了一个小示例程序（C#；为 0 到 15 之间的数字采用统一的 PRNG，并通过各种方式从中构造一个用于 0 到 6 之间的数字的 PRNG）：

using System;

class Rand {
    static Random r = new Random();

    static int Rand16() {
        return r.Next(16);
    }

    static int Rand7Naive() {
        return Rand16() % 7;
    }

    static int Rand7Float() {
        return (int)(Rand16() / 16.0 * 7);
    }

    // corrected
    static int Rand7RejectionNaive() {
        int n = 7, remainder = 16 % n, x, output;
        do {
            x = Rand16();
            output = x % n;
        } while (x >= 16 - remainder);
        return output;
    }

    // adapted to fit the constraints of this example
    static int Rand7RejectionJava() {
        int n = 7, x, output;
        do {
            x = Rand16();
            output = x % n;
        } while (x - output + 6 > 15);
        return output;
    }

    static void Test(Func<int> rand, string name) {
        var buckets = new int[7];
        for (int i = 0; i < 10000000; i++) buckets[rand()]++;
        Console.WriteLine(name);
        for (int i = 0; i < 7; i++) Console.WriteLine("{0}\t{1}", i, buckets[i]);
    }

    static void Main() {
        Test(Rand7Naive, "Rand7Naive");
        Test(Rand7Float, "Rand7Float");
        Test(Rand7RejectionNaive, "Rand7RejectionNaive");
    }
}

结果如下（粘贴到 Excel 中，并为单元格添加条件着色，使差异更加明显）：

在此处输入图像描述

现在我在上面的拒绝抽样中修复了我的错误，它应该可以正常工作（在它偏向 0 之前）。如您所见，float 方法并不完美，它只是以不同的方式分配有偏差的数字。

于 2012-08-01T12:08:52.553 回答
12

当随机数生成器 (RAND_MAX+1) 的输出数量不能被所需范围 (max-min+1) 整除时，就会出现问题。由于从随机数到输出会有一致的映射，因此某些输出将被映射到比其他输出更多的随机数。这与映射是如何完成的无关——您可以使用模数、除法、转换为浮点数，无论您能想出什么巫术，基本问题仍然存在。

问题的严重性非常小，要求不高的应用程序通常可以忽略它而侥幸成功。范围越小，RAND_MAX 越大，效果越不明显。

我采用了您的示例程序并对其进行了一些调整。首先我创建了一个rand只有0-255范围的特殊版本，以更好地展示效果。我对rangeRandomAlg2. 最后我将“球”的数量更改为 1000000 以提高一致性。你可以在这里看到结果：http: //ideone.com/4P4HY

请注意，浮点版本会产生两个紧密组合的概率，接近 0.101 或 0.097，两者之间没有任何关系。这是行动中的偏见。

我认为将此称为“Java 的算法”有点误导——我敢肯定它比 Java 古老得多。

int rangeRandomAlg2 (int min, int max)
{
    int n = max - min + 1;
    int remainder = RAND_MAX % n;
    int x;
    do
    {
        x = rand();
    } while (x >= RAND_MAX - remainder);
    return min + x % n;
}
于 2012-08-01T20:06:48.010 回答
6

很容易看出为什么这个算法会产生有偏差的样本。假设您的rand()函数从 set 返回统一整数{0, 1, 2, 3, 4}。如果我想用它来生成一个随机位0或1，我会说rand() % 2。集合{0, 2, 4}给了我0，集合{1, 3}给了我——很明显，我以 60% 和40% 的可能性1采样，根本不统一！01

要解决此问题，您必须确保所需的范围除以随机数生成器的范围，或者在随机数生成器返回的数字大于目标范围的最大可能倍数时丢弃结果。

在上面的例子中，目标范围是 2，适合随机生成范围的最大倍数是 4，所以我们丢弃任何不在集合中的样本{0, 1, 2, 3}并再次滚动。

于 2012-08-01T12:12:37.480 回答
3

到目前为止，最简单的解决方案是std::uniform_int_distribution<int>(min, max).

于 2012-08-03T15:05:00.303 回答
3

您谈到了涉及随机整数算法的两点：它是最优的，它是无偏的吗？

最佳

有很多方法可以定义“最优”算法。在这里，我们根据平均使用的随机位数来查看“最佳”算法。从这个意义上说，rand对于随机生成的数字来说，这是一种糟糕的方法，部分原因是它不一定需要生成随机位（因为RAND_MAX没有完全指定）*。相反，我们将假设我们有一个“真正的”随机生成器，它可以产生无偏且独立的随​​机位。

1976 年，DE Knuth 和 AC Yao 表明，任何仅使用随机位以给定概率生成随机整数的算法都可以表示为二叉树，其中随机位指示遍历树和每个叶子（端点）的方式对应一个结果。（Knuth 和 Yao，“非均匀随机数生成的复杂性”，在Algorithms and Complexity中，1976 年。）他们还给出了给定算法在此任务中平均需要的位数的界限。在这种情况下，均匀生成整数的最佳[0, n)算法将需要至少log2(n)和最多log2(n) + 2平均位。

在这个意义上，有很多优化算法的例子。其中之一是 J. Lumbroso (2013) 的Fast Dice Roller（在下面实现），也许另一个例子是2004 年数学论坛中给出的算法。另一方面，所有由 M. O'Neill 调查的算法不是最优的，因为它们依赖于一次生成随机位块。另请参阅我关于整数生成算法的说明。

下面显示了快速骰子滚子的实现；尽管它是在 JavaScript 中，而不是在 C 或 C++ 中，但它很容易适应任何一种语言，并且其想法是表明以最佳方式从位生成整数并不复杂。在代码中，(Math.random() < 0.5 ? 0 : 1)是 JavaScript 生成无偏随机位的方式。

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

不偏不倚

然而，任何同样无偏的最优整数生成器通常会在最坏的情况下永远运行，正如 Knuth 和 Yao 所展示的那样。回到二叉树，每个结果标签都留在二叉树中，因此 [0, n) 中的每个整数都可以以 1/n 的概率出现。但是如果 1/n 有一个非终止的二元展开式（如果 n 不是 2 的幂就是这种情况），这棵二叉树必然要么——</p> n

  • 具有“无限”深度，或
  • 在树的末端包括“拒绝”叶子，

在任何一种情况下，算法都不会在恒定时间内运行，并且在最坏的情况下会永远运行。（另一方面，当n是 2 的幂时，最佳二叉树将具有有限深度且没有拒绝节点。）快速骰子滚轮是使用“拒绝”事件来确保其无偏性的算法示例；请参阅上面代码中的注释。

一般来说n，没有办法在不引入偏差的情况下“修复”这种最坏情况的时间复杂度。例如，模减少（包括min + (rand() % (int)(max - min + 1))您的问题中的）相当于一棵二叉树，其中拒绝叶子被标记的结果替换 - 但由于可能的结果比拒绝叶子更多，因此只有一些结果可以代替拒绝离开，引入偏见。如果您在一定次数的迭代后停止拒绝，则会产生相同类型的二叉树 - 以及相同类型的偏差。（但是，根据应用程序，这种偏差可能可以忽略不计。随机整数生成也有安全方面的问题，在这个答案中讨论太复杂了。）

笔记

* 还有其他问题rand()。也许这里最严重的事实是 C 标准没有为rand().

于 2020-07-14T19:20:25.140 回答
1

不失一般性，在[a,b]上生成随机整数的问题可以简化为在[0,s)上生成随机整数的问题。从统一的 PRNG 生成有界范围内的随机整数的最新技术由以下最近的出版物表示：

Daniel Lemire，“区间内的快速随机整数生成”。ACM 翻译。模型。计算。模拟。29, 1, 第 3 条（2019 年 1 月）（ArXiv 草案）

Lemire 表明他的算法提供了无偏的结果，并且受到越来越受欢迎的非常快速的高质量 PRNG（如 Melissa O'Neill 的PCG 生成器）的推动，展示了如何快速计算结果，几乎一直避免缓慢的除法运算.

他的算法的一个示例性 ISO-C 实现randint()如下所示。在这里，我结合 George Marsaglia 的旧KISS64 PRNG 来演示它。出于性能原因，所需的 64×64→128 位无符号乘法通常最好通过机器特定的内在函数或直接映射到适当硬件指令的内联汇编来实现。

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>

/* PRNG state */
typedef struct Prng_T *Prng_T;
/* Returns uniformly distributed integers in [0, 2**64-1] */
uint64_t random64 (Prng_T);
/* Multiplies two 64-bit factors into a 128-bit product */
void umul64wide (uint64_t, uint64_t, uint64_t *, uint64_t *);

/* Generate in bias-free manner a random integer in [0, s) with Lemire's fast
   algorithm that uses integer division only rarely. s must be in [0, 2**64-1].

   Daniel Lemire, "Fast Random Integer Generation in an Interval," ACM Trans.
   Model. Comput. Simul. 29, 1, Article 3 (January 2019)
*/
uint64_t randint (Prng_T prng, uint64_t s) 
{
    uint64_t x, h, l, t;
    x = random64 (prng);
    umul64wide (x, s, &h, &l);
    if (l < s) {
        t = (0 - s) % s;
        while (l < t) {
            x = random64 (prng);
            umul64wide (x, s, &h, &l);
        }
    }
    return h;
}

#define X86_INLINE_ASM (0)

/* Multiply two 64-bit unsigned integers into a 128 bit unsined product. Return
   the least significant 64 bist of the product to the location pointed to by
   lo, and the most signfiicant 64 bits of the product to the location pointed
   to by hi.
*/
void umul64wide (uint64_t a, uint64_t b, uint64_t *hi, uint64_t *lo)
{
#if X86_INLINE_ASM
    uint64_t l, h;
    __asm__ (
        "movq  %2, %%rax;\n\t"  // rax = a
        "mulq  %3;\n\t"         // rdx:rax = a * b
        "movq  %%rax, %0;\n\t"  // l = (a * b)<31:0>
        "movq  %%rdx, %1;\n\t"  // h = (a * b)<63:32>
        : "=r"(l), "=r"(h)
        : "r"(a), "r"(b)
        : "%rax", "%rdx");
    *lo = l;
    *hi = h;
#else // X86_INLINE_ASM
    uint64_t a_lo = (uint64_t)(uint32_t)a;
    uint64_t a_hi = a >> 32;
    uint64_t b_lo = (uint64_t)(uint32_t)b;
    uint64_t b_hi = b >> 32;

    uint64_t p0 = a_lo * b_lo;
    uint64_t p1 = a_lo * b_hi;
    uint64_t p2 = a_hi * b_lo;
    uint64_t p3 = a_hi * b_hi;

    uint32_t cy = (uint32_t)(((p0 >> 32) + (uint32_t)p1 + (uint32_t)p2) >> 32);

    *lo = p0 + (p1 << 32) + (p2 << 32);
    *hi = p3 + (p1 >> 32) + (p2 >> 32) + cy;
#endif // X86_INLINE_ASM
}

/* George Marsaglia's KISS64 generator, posted to comp.lang.c on 28 Feb 2009
   https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
*/
struct Prng_T {
    uint64_t x, c, y, z, t;
};

struct Prng_T kiss64 = {1234567890987654321ULL, 123456123456123456ULL,
                        362436362436362436ULL, 1066149217761810ULL, 0ULL};

/* KISS64 state equations */
#define MWC64 (kiss64->t = (kiss64->x << 58) + kiss64->c,            \
               kiss64->c = (kiss64->x >> 6), kiss64->x += kiss64->t, \
               kiss64->c += (kiss64->x < kiss64->t), kiss64->x)
#define XSH64 (kiss64->y ^= (kiss64->y << 13), kiss64->y ^= (kiss64->y >> 17), \
               kiss64->y ^= (kiss64->y << 43))
#define CNG64 (kiss64->z = 6906969069ULL * kiss64->z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)
uint64_t random64 (Prng_T kiss64)
{
    return KISS64;
}

int main (void)
{
    int i;
    Prng_T state = &kiss64;

    for (i = 0; i < 1000; i++) {
        printf ("%llu\n", randint (state, 10));
    }
    return EXIT_SUCCESS;
}
于 2020-01-26T03:31:15.777 回答
0

如果你真的想得到一个完美的生成器，假设你拥有的 rand() 函数是完美的，你需要应用下面解释的方法。

我们将创建一个随机数 r，从 0 到 max-min=b-1，然后很容易移动到您想要的范围，只需 r+min

我们将创建一个 b < RAND_MAX 的随机数，但是可以很容易地采用该过程来为任何基数生成一个随机数

程序：

  1. 取一个原始 RAND_MAX 大小的随机数 r，不进行任何截断
  2. 以 b 为基数显示此数字
  3. 取这个数字的前 m=floor(log_b(RAND_MAX)) 个数字，用于从 0 到 b-1 的 m 个随机数
  4. 将每个按 min（即 r+min）移动，以使它们进入您想要的范围（min，max）

由于 log_b(RAND_MAX) 不一定是整数，因此表示中的最后一位数字被浪费了。

仅使用 mod (%) 的原始方法完全被

(log_b(RAND_MAX) - floor(log_b(RAND_MAX)))/ceil(log_b(RAND_MAX))

你可能同意的不是那么多，但如果你坚持要精确，那就是程序。

于 2020-01-25T14:33:20.707 回答