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

# はじめに

[「値をランダムに分配するにはどうすればいいのか」](http://qiita.com/khsk/items/4dde8384832a405e6e5d)が元ネタです。数学には疎いものの、**「あれ？正規分布ってまさにこの用途じゃね？」**と思いついたのでググりながらいろいろ試してみました。

# 正規分布って何？

Wikipedia先生に聞きましょう。

- [Wikipedia - 正規分布](https://ja.wikipedia.org/wiki/%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83)

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

![720px-Normal_Distribution_PDF.svg.png](https://qiita-image-store.s3.amazonaws.com/0/25060/e097073f-56c5-b737-a5fa-ebd7480c39b5.png)

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

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

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

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

- [PHP Manual - stats_rand_gen_normal](http://php.net/manual/ja/function.stats-rand-gen-normal.php)

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

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

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

- [Wikipedia - ボックス＝ミュラー法](https://ja.wikipedia.org/wiki/%E3%83%9C%E3%83%83%E3%82%AF%E3%82%B9%EF%BC%9D%E3%83%9F%E3%83%A5%E3%83%A9%E3%83%BC%E6%B3%95)
- [Wikipedia - Box-Muller transform](https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform)

# 関数の実装

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

```php
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以上のときのみ計算を実行する<br />(中央値が1以上のときのみ計算を実行する)
- 分散が大きい時には前に大きな数が偏るので`shuffle`を最後に実行する<br />(必要に応じて[`mt_shuffle`](http://qiita.com/mpyw/items/c644dab7d067626fc629#mt_array_rand-mt_shuffle)に替えてください)

```php:分配用関数の定義
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));` で実行してみます。 

```php:分散が0.2のとき
Array
(
    [0] => 100
    [1] => 100
    [2] => 99
    [3] => 100
    [4] => 100
    [5] => 101
    [6] => 100
    [7] => 100
    [8] => 100
    [9] => 100
)
```

```php:分散が1.0のとき
Array
(
    [0] => 101
    [1] => 100
    [2] => 100
    [3] => 102
    [4] => 101
    [5] => 100
    [6] => 98
    [7] => 98
    [8] => 99
    [9] => 101
)
```

```php:分散が5.0のとき
Array
(
    [0] => 93
    [1] => 98
    [2] => 100
    [3] => 100
    [4] => 103
    [5] => 90
    [6] => 114
    [7] => 97
    [8] => 106
    [9] => 99
)
```

```php:分散が20.0のとき
Array
(
    [0] => 105
    [1] => 68
    [2] => 109
    [3] => 136
    [4] => 133
    [5] => 90
    [6] => 96
    [7] => 87
    [8] => 123
    [9] => 53
)
```

```php:分散が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_rand`は`rand`に比べるとマシになっているものの、暗号学的に安全なレベルの乱数では無いようです。この用途で使うには、従来は以下の`openssl_rand`のように自分で関数を作る必要がありました。

- [Qiita - array_randやshuffleの乱択精度を改善しよう](http://qiita.com/mpyw/items/c644dab7d067626fc629)

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

- [PHP Manual - random_int](http://php.net/manual/ja/function.random-int.php)