Box-Muller法を用いて正規分布に従う分配を行う

  • 19
    Like
  • 3
    Comment
More than 1 year has passed since last update.

はじめに

「値をランダムに分配するにはどうすればいいのか」が元ネタです。数学には疎いものの、「あれ?正規分布ってまさにこの用途じゃね?」と思いついたのでググりながらいろいろ試してみました。

正規分布って何?

Wikipedia先生に聞きましょう。

正直自分もあんまり分かってないんですが、素人なりに大雑把に言うと「中心にいい感じに偏ってくれる確率分布」だと認識しています。さらに 「左右対称に」 というべきでしょうか。

720px-Normal_Distribution_PDF.svg.png

偏り方は与える $\mu$ および $\sigma$ によって変わってくるようです。

正規分布に従う乱数の生成方法

2通りの方法を挙げてみます。

PHPネイティブの関数を使う

標準関数の手広さがウリのPHPでも、これは残念ながらカバーしてくれないようです…
(Pythonとかなら標準であるんですかね?)

一様乱数を正規乱数に変換する

世の中には頭のいいことを思いつく人もいるんですね!

関数の実装

Wikipedia英語版のC++の例を参考にして書いてみました。意味が分からなかったところは端折ったりしてます。

function generate_norm(float $average = 0.0, float $variance = 1.0) {
    static $z1, $z2, $mt_max, $ready = true;
    if ($mt_max === null) {
        $mt_max = mt_getrandmax();
    }
    $ready = !$ready;
    if ($ready) {
        return $z2 * $variance + $average;
    }
    $u1 = mt_rand(1, $mt_max - 1) / $mt_max;
    $u2 = mt_rand(1, $mt_max - 1) / $mt_max;
    $v1 = sqrt(-2 * log($u1));
    $v2 = 2 * M_PI * $u2;
    $z1 = $v1 * cos($v2); 
    $z2 = $v1 * sin($v2);
    return $z1 * $variance + $average;
}

関数の利用

1000を10人に分配してみましょう。分散は固定値にし、平均だけ1回ごとに書き換えてみます。

  • 最低でも1は全員に割り振れるときのみ計算を実行する
  • 無限ループもしくは長いループを防ぐために $\frac{1}{2\sqrt{2\pi\sigma^2}}$が1以上のときのみ計算を実行する
    (中央値が1以上のときのみ計算を実行する)
  • 分散が大きい時には前に大きな数が偏るのでshuffleを最後に実行する
    (必要に応じてmt_shuffleに替えてください)
分配用関数の定義
function distribute(int $total, int $count, float $variance = 1.0) {
    switch (true) {
        case $count < 1:
        case $total < $count:
        case $variance < 0:
        case 1 / (2 * sqrt(2 * M_PI * $variance)) > 1:
            return false;
    }
    $result = [];
    for ($i = 0; $i < $count - 1; ++$i) {
        do {
            $result[$i] = round(generate_norm($total / ($count - $i), $variance));
        } while ($result[$i] < 1 || $total - $result[$i] < $count - $i);
        $total -= $result[$i];
    }
    $result[$count - 1] = $total;
    shuffle($result);
    return $result;
}

逐次 print_r(distribute(1000, 10, $variance)); で実行してみます。

分散が0.2のとき
Array
(
    [0] => 100
    [1] => 100
    [2] => 99
    [3] => 100
    [4] => 100
    [5] => 101
    [6] => 100
    [7] => 100
    [8] => 100
    [9] => 100
)
分散が1.0のとき
Array
(
    [0] => 101
    [1] => 100
    [2] => 100
    [3] => 102
    [4] => 101
    [5] => 100
    [6] => 98
    [7] => 98
    [8] => 99
    [9] => 101
)
分散が5.0のとき
Array
(
    [0] => 93
    [1] => 98
    [2] => 100
    [3] => 100
    [4] => 103
    [5] => 90
    [6] => 114
    [7] => 97
    [8] => 106
    [9] => 99
)
分散が20.0のとき
Array
(
    [0] => 105
    [1] => 68
    [2] => 109
    [3] => 136
    [4] => 133
    [5] => 90
    [6] => 96
    [7] => 87
    [8] => 123
    [9] => 53
)
分散が2000.0のとき
Array
(
    [0] => 60
    [1] => 2
    [2] => 325
    [3] => 6
    [4] => 4
    [5] => 1
    [6] => 434
    [7] => 157
    [8] => 10
    [9] => 1
)

これはいい感じじゃないですかね!?あとはwhileを最小値および最大値も加えたループ条件にすれば、実用的に使えるものになるんじゃないでしょうか。

備考

mt_randrandに比べるとマシになっているものの、暗号学的に安全なレベルの乱数では無いようです。この用途で使うには、従来は以下のopenssl_randのように自分で関数を作る必要がありました。

ところが、PHP7以降はrandom_intという関数が利用できます。PHP7の普及に伴い、わざわざ精度の悪いmt_randを使う必要性は薄れていくかもしれませんね。