いつからその方法で偏りのない乱数が得られると錯覚していた?

私はつい最近まで勘違いしていました。

ここのページに書いてあるような方法で、一様分布する整数が得られると。

int random(int n)
{
    return (int)(( rand() / (RAND_MAX + 1.0) ) * n);
}


この方法、一見すると実に一様分布が得られそうに見えるんですよね。

どういう思考回路を通っているかというのを自己分析すると、次のような感じです。


1. rand() では 0〜RAND_MAX のランダムな整数が得られる。

2. それを RAND_MAX + 1 で割ると、[0, 1) に一様分布する実数が得られる。

3. [0, 1) の一様な実数を n 倍して小数点以下を切り捨てたら、0 から n-1 に一様分布する整数が得られる。


これの罠なところは、1 と(特に)3 が正しいというところだと思います。

ただ、2 がダウト。

思いっきりダウト。


そもそも、整数を整数で割って一様分布する実数が得られるわけがないのですが。

RAND_MAX というのは環境によって 32767 だったり 21億なんとかだったりして、どちらにしても大きい数字なので感覚がつかみにくい。


ここで、C言語の仕様とかそういうのは一旦置いておいて、RAND_MAX が 5 である世界を考えてみます。

すると、rand() の結果は 0, 1, 2, 3, 4, 5 をそれぞれ 1/6 の確率で返すというものになります。

つまり、サイコロを振って、その目より 1 少ない数を結果とするというのをイメージすればOKです。


さて、この rand() (=サイコロの目から 1 を引いたもの)を使って、0〜3 の範囲に一様分布する整数を得ることはできるでしょうか?

上の手順でいうと、(rand() / (RAND_MAX + 1.0)) * n の要領で、サイコロの目から 1 を引いたものを 6 で割って、それに 4 を掛け、小数点以下を切り捨てれば 0〜3 の整数が得られるということになります。


たとえば、父・母・兄・妹の 4 人家族がいて、毎日皿洗いの当番をこの方法で決めるとします。

最初の日、サイコロを振るとたまたま 1 が出たので、それから 1 を引いて 6 で割って 4 を掛けると (0 / 6) * 4 = 0 → 0 → お父さん(ここではお父さん: 0、お母さん: 1… という順番とする)。

次の日、サイコロを振るとたまたま 2 が出たので、それから 1 を引いて 6 で割って 4 を掛けると (1 / 6) * 4 = 0.6666... → 0 → お父さん。

その次の日、サイコロを振るとたまたま 3 が出たので、それから 1 を引いて 6 で割って 4 を掛けると (2 / 6) * 4 = 1.3333... → 1 → お母さん。


結局、サイコロの目によって次のような割り当てになります。

1, 2 → お父さん(0)
3 → お母さん(1)
4, 5 → お兄さん(2)
6 → 妹 (3)

これって明らかに不公平ですよね。


ちゃんと公平にしようと思ったら、サイコロを振って5 か 6 が出たらやり直すというのが手っ取り早いやり方です。

(運によっては、5 と 6 ばっかり出てサイコロを何回も振り直さないといけないというリスクはありますが。ちなみに振り直しを減らそうと思ったら、複雑になりますが 2 回サイコロを振って数字を組み合わせて 0 〜 35 の数字にして、32〜35 ならやり直し、0〜31 ならそれで 2 回分決めるなどということも考えられます)


で、RAND_MAX が 32767 で、上の方法で 0 〜 5 の乱数を作る場合は、次のような対応になります。

0〜5461 (5462通り) : 0
5462〜10922 (5461通り) : 1
10923〜16383 (5461通り): 2
16384〜21845 (5462通り): 3
21846〜27306 (5461通り): 4
27307〜32767 (5461通り): 5

つまり、0 と 3 の出る確率がそれぞれ 5462/32768 で、1, 2, 4, 5 の出る確率がそれぞれ 5461/32768 になります。



RAND_MAX が 32767 の VC++(VS 2012) で次のようなコードを実行してみました。

#include <cstdio>
#include <cstdlib>
#include <ctime>

#define DICE_MAX 6

int uniform_int(int n) {
    return (int)(((double)rand() / (RAND_MAX + 1)) * DICE_MAX);
}

int main()
{
    int larger_count = 0;
    for (int i = 0; i < 100; ++i) {
        srand((unsigned int)time(0));
        int dice_array[DICE_MAX] = { 0 };
        for (int j = 0; j < 1000000000; ++j) {
            dice_array[uniform_int(DICE_MAX)]++;
        }
        for (int j = 0; j < DICE_MAX; ++j) {
            printf("dice_array[%d]: %d\n", j, dice_array[j]);
        }
        if (dice_array[0] > dice_array[1]) {
            printf("dice_array[0] is larger than dice_array[1]\n");
            ++larger_count;
        } else {
            printf("dice_array[0] is smaller than dice_array[1]\n");
        }
    }
    printf("# of times where dice_array[0] was larger than dice_array[1]: %d\n", larger_count);
}

上に書いたような方法で 0〜5 に一様分布させる(つもりの)整数を10億回生成するということを 100回繰り返して、0 が 1 より多く出た回数を数えるというものです。

(一定間隔で time で srand() しているとか、そもそも乱数の周期が 21億ぐらいじゃないかという筋の悪さは置いといて)

結果は今回は次のようになりました。

...
dice_array[0]: 166699112
dice_array[1]: 166660744
dice_array[2]: 166661006
dice_array[3]: 166673554
dice_array[4]: 166653122
dice_array[5]: 166652462
dice_array[0] is larger than dice_array[1]
dice_array[0]: 166701419
dice_array[1]: 166655374
dice_array[2]: 166667455
dice_array[3]: 166671795
dice_array[4]: 166652118
dice_array[5]: 166651839
dice_array[0] is larger than dice_array[1]
# of times where dice_array[0] was larger than dice_array[1]: 100

100 回中 100 回すべて、0 のほうが 1 よりも多いという結果でした。

このように、10 億回もやると、ほぼ確実に 0 と 3 が多くなります。

まあ、一般的に 5461 : 5462 の偏りが重要になる場面がどれだけあるかはわかりませんが。

いずれにせよ、この方法で偏りのない乱数が得られるというのは勘違いということです。


じゃあ、どうせ偏りが出るなら剰余でもいいんじゃない? という怠惰な発想もあるかと思いますが、そちらは乱数の実装によって下位ビットに偏りがあるという問題が表面化するためおすすめできません。

このあたりは、良い乱数・悪い乱数というページに詳しく書かれています。


というわけで、rand() を使って偏りなく一様分布する整数を求めようと思ったら、ひと工夫必要になります。

上の uniform_int() 関数を次のように書き換えます。

int uniform_int(int n) {
    int adjusted_max = (RAND_MAX + 1) - (RAND_MAX + 1) % n;
    int r;
    do {
       r = rand();
    } while (r >= adjusted_max);
    return (int)(((double)r / adjusted_max) * n);
}

RAND_MAX + 1 が n で割り切れない場合、乱数がその端数分の範囲になったら捨ててやり直すというものです。


「いや、そんな原始的なことしなくても世の中には mt19937(ry」という意見もあるでしょうが、世の中にはいろいろな大人の事情で高度っぽいライブラリや先端っぽい技術が使えなくて血の涙を流している人も多いと思います。

そういう人は、せめてこの方法を使えば偏りをなくすことだけはできるのでおすすめかなぁと。

まあ、高度っぽいライブラリや先端っぽい技術を使えないような職場では、往々にしてこんな何千分の一の乱数の偏りなんかよりずっと大きな問題が山積していて、これだけ解決しても焼け石に水かもしれませんが。


この間違いをしているページは多く見かけます。

私は、id:niam さんのCのrand()よりmt19937の方が速いことがあるという話を読んで気付きました。


上で最初に挙げたページもそうですが、「ランダムの偏り」というこの記事も勘違いしていますね。

乱数生成の方法」では、「NがRAND_MAXにくらべて十分小さいことを仮定しています」とあるので、この偏りを意識されているのだと思います。


そういうわけで、世の中のゲーム中でのサイコロには 1 と 4 が微妙に出やすいものがけっこうありそうなので、それを使ったチートをやろうと思ったらできるかもしれません(難しいとは思いますが)。


追記: 偏りのない方法として書いたもので剰余を使っていて、ビットの偏りによる問題がある可能性があったため修正しました。