c++ - 範囲内で偏りのないランダムな整数を生成するための最適なアルゴリズムは何ですか?

Question

このStackOverflowの質問では：

範囲からランダムな整数を生成する

受け入れられた答えは、与えられたminとの間のランダムな整数を生成するための次の式を示唆しています。maxminmax

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);

多数の「バケツ」にどれだけ均等に「ボール」が「落ち」、分散されているかを確認します。1 つは浮動小数点外挿のテストで、もう 1 つはwhile-loop アルゴリズムのテストです。しかし、「ボール」（および「バケツ」）の数によって結果が異なることが判明したため、勝者を簡単に選ぶことはできませんでした。作業コードは、この Ideone ページにあります。たとえば、バケットが 10 個でボールが 100 個の場合、バケット間の理想的な確率からの最大偏差は、whileループ アルゴリズムよりも浮動小数点外挿の方が小さくなります (それぞれ 0.04 と 0.05) が、ボールが 1000 個の場合、while-loop アルゴリズムはより小さく (0.024 および 0.011)、10000 個のボールを使用すると、浮動小数点外挿は再びより良く (0.0034 および 0.0053) 実行され、一貫性があまりありません。whileどのアルゴリズムも一貫して他のアルゴリズムよりも優れた均一分布を生成しない可能性を考えると、 -loop アルゴリズムよりも高速に実行されるように見えるため、浮動小数点外挿に傾倒します。では、浮動小数点外挿アルゴリズムを選択しても問題ありませんか、それとも私のテスト/結論は完全に正しくありませんか?

Answer 1

問題は、モジュロ演算を行っていることです。RAND_MAXモジュラスで割り切れる場合は問題ありませんが、通常はそうではありません。非常に不自然な例として、RAND_MAX11 でモジュラスが 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 よりもわずかに確率が高くなります。

これを解決する 1 つのオプションは、棄却サンプリングです。上記の 9 と 10 の数字を許可しないことで、結果の分布を再び均一にすることができます。トリッキーな部分は、それを効率的に行う方法を考え出すことです。非常に優れた例 (なぜ機能するのかを理解するのに 2 日かかった例) は、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 メソッドはまったく完全ではありません。バイアスされた数値を別の方法で分散しているだけです。

Answer 2

この問題は、乱数ジェネレーターからの出力数 (RAND_MAX+1) が目的の範囲 (最大-最小+1) で割り切れない場合に発生します。乱数から出力への一貫したマッピングがあるため、一部の出力は他よりも多くの乱数にマッピングされます。これは、マッピングがどのように行われるかに関係なく、モジュロ、除算、浮動小数点への変換など、思いつくブードゥーを使用できますが、基本的な問題は残ります。

問題の規模は非常に小さく、要求の厳しいアプリケーションでは、通常、無視して問題を解決できます。範囲が小さく、RAND_MAX が大きいほど、効果は目立たなくなります。

私はあなたのサンプルプログラムを取り、少し調整しました。最初にrand、効果をよりよく示すために、範囲が 0 ～ 255 しかない の特別なバージョンを作成しました。にいくつかの微調整を行いましたrangeRandomAlg2。最後に、一貫性を向上させるために「ボール」の数を 1000000 に変更しました。ここで結果を見ることができます: http://ideone.com/4P4HY

浮動小数点バージョンは、0.101 または 0.097 に近い 2 つの密接にグループ化された確率を生成することに注意してください。これがバイアスの作用です。

これを「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;
}

Answer 3

このアルゴリズムが偏ったサンプルを生成する理由は簡単にわかります。rand()関数が set から一様整数を返すとします{0, 1, 2, 3, 4}。これを使用してランダム ビット0またはを生成したい場合は1、 と言うでしょうrand() % 2。セット{0, 2, 4}は私0に を与え、セット{1, 3}は私に与えます1- 明らかに、私0は 60% と140% の可能性でサンプリングし、まったく均一ではありません!

これを修正するには、目的の範囲が乱数ジェネレーターの範囲を分割することを確認するか、そうでなければ、乱数ジェネレーターがターゲット範囲の可能な最大倍数よりも大きい数値を返すたびに結果を破棄する必要があります。

上記の例では、ターゲット範囲は 2 で、ランダム生成範囲に収まる最大の倍数は 4 であるため、セットに含まれていないサンプルを破棄して、{0, 1, 2, 3}もう一度ロールします。

Answer 4

最も簡単な解決策は、std::uniform_int_distribution<int>(min, max).

Answer 5

Without loss of generality, the problem of generating random integers on [a, b] can be reduced to the problem of generating random integers on [0, s). The state of the art for generating random integers on a bounded range from a uniform PRNG is represented by the following recent publication:

Daniel Lemire,"Fast Random Integer Generation in an Interval." ACM Trans. Model. Comput. Simul. 29, 1, Article 3 (January 2019) (ArXiv draft)

Lemire shows that his algorithm provides unbiased results, and motivated by the growing popularity of very fast high-quality PRNGs such as Melissa O'Neill's PCG generators, shows how to the results can be computed fast, avoiding slow division operations almost all of the time.

彼のアルゴリズムの 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;
}