C语言 在一个范围内生成一个无偏随机整数的最优算法是什么?

wpx232ag  于 2023-02-15  发布在  其他
关注(0)|答案(7)|浏览(109)

在此StackOverflow问题中:
Generating random integer from a range
可接受的答案提出了用于生成给定minmax之间的随机整数的以下公式,其中minmax包括在该范围内:

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

但它也说
这仍然是 * 轻微 * 偏向于较低的数字...这也是可能的,以扩大它,使它消除偏见。
但它没有解释为什么它偏向于较低的数字,也没有解释如何消除这种偏向。所以,问题是:这是在(有符号的)范围内生成随机整数而不依赖于任何花哨的东西,仅依赖于rand()函数的最优方法吗?如果是最优的,如何去除偏差?

    • 编辑:**

我刚刚测试了@Joey建议的while-loop算法与浮点外推法:

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

为了了解有多少均匀的"球"落入并分布在多个"桶"中,一个测试用于浮点外推,另一个测试用于while循环算法。(和"buckets"),所以我不能轻易地选出一个赢家。工作代码可以在this Ideone page中找到。例如,在10个桶和100个球的情况下,对于浮点外推,桶之间与理想概率的最大偏差小于对于while循环算法(分别为0.04和0.05),但对于1000个球,while-循环算法的最大偏差较小(0.024和0.011),并且具有10000个球,浮点外推再次表现得更好(0.0034和0.0053)等等,没有太多的一致性。考虑到没有一种算法一致地产生比另一种算法更好的均匀分布的可能性,使我倾向于浮点外推,因为它似乎比while循环算法执行得更快。那么,选择浮点外推算法还是我的测试/结论并不完全正确

pgx2nnw8

pgx2nnw81#

问题是你正在做一个模运算。如果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之前)。正如你所看到的,浮点方法一点也不完美,它只是以不同的方式分布有偏向的数字。

aydmsdu9

aydmsdu92#

当随机数生成器的输出数(兰德_MAX+1)不能被所需范围整除(max-min+1)。由于从随机数到输出的Map是一致的,因此某些输出将比其他输出Map到更多的随机数。这与Map是如何完成的无关-您可以使用模、除法、浮点转换、不管你能想出什么巫术,基本问题依然存在。
问题的严重程度非常小,要求不高的应用程序通常可以忽略它。范围越小,兰德_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;
}
nwnhqdif

nwnhqdif3#

很容易理解为什么这个算法会产生一个有偏差的样本,假设你的rand()函数从集合{0, 1, 2, 3, 4}返回均匀整数,如果我想用它来产生一个随机位01,我会说rand() % 2,集合{0, 2, 4}给我0,集合{1, 3}得到1--很明显,我以60%的可能性采样0,以40%的可能性采样1,一点也不均匀!
要解决这个问题,您必须确保您想要的范围除以随机数生成器的范围,否则每当随机数生成器返回的数字大于目标范围的最大可能倍数时,就 * 丢弃 * 结果。
在上面的例子中,目标范围是2,适合随机生成范围的最大倍数是4,因此我们丢弃不在集合{0, 1, 2, 3}中的任何样本并再次滚动。

fykwrbwg

fykwrbwg4#

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

9o685dep

9o685dep5#

您已经涉及到随机整数算法的两点:它是“最优的”吗?它是“无偏的”吗?

最佳

定义“最优”算法的方法有很多种。这里我们从平均使用的随机位数的Angular 来看待“最优”算法。从这个意义上讲,rand是一种用于随机生成数的较差方法,因为除了rand()的其他问题之外,它不一定需要生成随机位(因为RAND_MAX没有被精确地指定)。相反,我们将假设我们具有能够产生无偏且独立的随机位的“真”随机生成器。
在1976年,D. E. Knuth和A. C. Yao证明了任何只使用随机位以给定概率产生随机整数的算法都可以表示为二叉树,其中随机位指示遍历树和每个叶子的方式(终点)对应于结果。(Knuth和Yao,“非均匀随机数生成的复杂性”,《算法与复杂性》,1976)。他们还给出了给定算法平均需要的位数的界限。在这种情况下,在[0, n)中均匀生成整数的 * 最优 * 算法,将需要平均至少log2(n)位和最多log2(n) + 2
在这个意义上,有许多“最优”算法的例子,请看我的以下回答:

  • 如何在不浪费位的情况下从随机位流生成范围[0,n]内的随机整数?

无偏见

然而,任何也是“无偏”的“最优”整数生成器通常将在最坏情况下永远运行,也如Knuth和Yao所示。回到二叉树,n结果标签中的每一个在二叉树中标记叶,使得[0,但是,如果1/n具有非终止二进制扩展(如果n不是2的幂,则将是这种情况),则该二叉树将必然是-

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

在这两种情况下,算法都不会在恒定的时间内运行,在最坏的情况下,算法将永远运行下去(另一方面,当n是2的幂时,最优二叉树将具有有限的深度,并且没有拒绝节点)。
对于一般的n,没有办法在不引入偏差的情况下“修正”这种最坏情况下的时间复杂度。(包括你的问题中的min + (rand() % (int)(max - min + 1)))等价于二叉树,在二叉树中,拒绝叶被替换为标记的结果-但是由于可能的结果比拒绝叶多,只有一部分结果可以代替拒绝叶,这就引入了偏差。2如果你在一定次数的迭代后停止拒绝,就会产生同样的二叉树--同样的偏差。(然而,根据应用程序的不同,这种偏差可能可以忽略不计。随机整数生成还有安全方面的问题,这些问题太复杂,无法在本答案中讨论。)

6ss1mwsb

6ss1mwsb6#

不失一般性,在[a,b]上生成随机整数的问题可以简化为在[0,s)上生成随机整数的问题。用于从均匀PRNG在有界范围上生成随机整数的技术状态由以下最近的出版物表示:
丹尼尔Lemire,“在区间内快速生成随机整数”。ACM Trans. Model. Comput. Simul. 29,1,Article 3(January 2019)(ArXiv draft
Lemire展示了他的算法提供了无偏差的结果,并受快速高质量PRNG(如Melissa奥尼尔的PCG generators)日益流行的推动,展示了如何快速计算结果,几乎总是避免缓慢的除法运算。
下面的randint()展示了他的算法的一个示例性ISO-C实现。这里我结合乔治Marsaglia的老KISS 64 PRNG来演示它。出于性能原因,所需的64×64→128位无符号乘法通常最好通过机器特定的intrinsic或内联汇编来实现,这些intrinsic或内联汇编直接Map到相应的硬件指令。

#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;
}
rt4zxlrg

rt4zxlrg7#

如果你真的想得到一个完美的生成器,假设你拥有的rand()函数是完美的,你需要应用下面解释的方法。
我们将创建一个随机数r,从0到max-min = b-1,然后很容易移动到您想要的范围,只需取r + min
我们将创建一个随机数,其中b〈RAND_MAX,但该过程可以很容易地采用任何基数的随机数
程序:
1.取原始RAND_MAX大小的随机数r,不进行任何截断
1.以底数b显示此数字
1.对于从0到b-1的m个随机数,取该数的前m = floor(log_b(RAND_MAX))位
1.将每个值移动最小值(即r + min),以使它们进入所需的范围(最小值,最大值
由于log_b(RAND_MAX)不一定是整数,因此浪费了表示中的最后一位。
原来只使用mod(%)的方法错误在于

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

你可能会同意这并不多,但如果你坚持要精确,这就是程序。

相关问题