【C#】正規乱数 (1)：Box-Mullar 法, 中心極限定理を利用した方法

正規分布に従った乱数の高速な生成アルゴリズムである Ziggurat 法の詳説記事を見た。
詳説 Ziggurat 法～正規分布・指数分布乱数の高速生成 - 屋根裏工房改

非常に詳しい解説がなされている。
知識として Ziggurat 法が速いらしいというのは聞いていたけれど、試したことはなかった。

知ってる他の方法も整理してから、試し較べてみたい。

ボリューミーだったので分割。この記事は (1)。

Box-Mullar 法

C# には標準に正規乱数が用意されていないので実装が各所で無限に繰り返されていそうな、デファクトスタンダードと思われる手法。

using static System.Math;

public static (double Z1, double Z2) NextBoxMullar(Random rand)
{
    double x = rand.NextDouble();
    double y = rand.NextDouble();
    double r = Sqrt(-2 * Log(1 - x));
    (double sin, double cos) = SinCos(y * (PI * 2));
    return (r * sin, r * cos);
}

ふたつ値が出るが何らかの都合でひとつずつ出力を得たい場合は、状態に持つか、あるいは片方を捨てる / 計算しないというやり方もある。この場合、得る値はそれぞれ独立なので、どちらを選んでもよい
速度面（Sin, Cos, Log といった数学関数を必ず使うため、それらを省く方法には速度が劣る）や精度面（裾の精度がやや悪いらしい ※ 聞きかじり）で別のアルゴリズムが選択されうるらしい
冒頭の参考リンクでは、独自の乱数生成器使用に伴い [0-1) の一様乱数を得るとき符号なし 64bit 整数から変換する hack を用いているが、標準のRandomを用いる場合はNextDouble()でよい

中心極限定理を利用した方法

一様乱数を12回足して6を引くと正規乱数を得ることができるという方法。
足して引く簡単な操作で正規化されている（とみなせる）のは感嘆。

// Random rand;
double r = rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble()
    + rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble()
    + rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble()
    - 6d;

数学関数を用いないが、現代のコンピュータでは Box-Mullar 法の方が速い
Box-Mullar 法も（原理はともかく）利用が容易なため、この方法の選択は限られると思う

コクのある乱数法

中心極限定理を利用した方法のバリエーションとしてコクのある乱数という手法がある。
参考：乱数にコクを出す方法について - Togetter

おおもとの手法は一様乱数を n 回足して n で割る方法（n = 5 が推奨されている）。
そのままだと n によって分散が異なる、n は整数に限られるという性質があり、以下のリンクでは正規化・連続化を行っている。
参考：コクのある乱数を使いやすくする - Qiita

上記について、一部の計算を N 指定時に追い出すためにオブジェクトとし、C# に落とし込んだのが以下。

using static System.Math;

public sealed class Koku
{
    private readonly Random rand;
    private double a, sd;
    private int ni;

    public double N
    {
        get => ni + a;
        set
        {
            if (value < 1d) throw new ArgumentOutOfRangeException("N < 1");

            double n = value;
            ni = (int)n;
            a = n - ni;
            var sigma = (ni / 3d) + (a * a / 3d);
            sd = Sqrt(sigma);
        }
    }

    public Koku()
    {
        rand = new();
        N = 5;
    }

    public double Next()
    {
        double s = 0d;
        for (int i = 0; i < ni; i++)
        {
            s += rand.NextStandardUniform(); // ※ 注
        }
        s += rand.NextStandardUniform() * a; // ※ 注

        return s / sd;
    }
}

※注 NextStandardUniform()は [min: -1, max: 1) の一様乱数を得る手続きとする
- C# では NextDouble() に min, max を受け取るオーバーロードがない
- 基本的な方法は、2 * NextDouble() - 1;
- この範囲を取る方法として、冒頭の参考リンクでは、符号なし 64bit 整数から変換する hack が紹介されている。

n = 5 の場合、Box-Mullar 法より速い（と思う。自分の環境では速い）
正規分布とは分布がそれなりに異なる

ヒストグラムを以下。

正規分布と較べると勾配がやや緩やかな感じがする。

Box-Mullar 法で取ったヒストグラムと重ねてみると、天辺付近と裾が小さく中腹に寄っている感じ。

正規乱数という文脈だと精度がという話になるけれど、性質として捉えるなら話は別。
もともとビジュアルプログラミングで用いる手法ということで、これがいい味に繋がるコクということなんだと思う。

有識者の見解も見るに、以下のような落とし所となるよう。
・正規乱数の代用としては良くない
・コクを求めるならその限りでない

(2) へ。

てくメモ

trivial な notes

【C#】正規乱数 (1)：Box-Mullar 法, 中心極限定理を利用した方法

Box-Mullar 法

中心極限定理を利用した方法

コクのある乱数法