LoginSignup
21
18

More than 5 years have passed since last update.

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

Last updated at Posted at 2015-07-02

はじめに

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

正規分布って何?

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を使う必要性は薄れていくかもしれませんね。

21
18
3

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
21
18