先日、学生から「Octaveではrandn()で正規分布に従う乱数が作れるが、C/C++にはそのようなものが用意されていないようで困っている」との相談がありました。聞いてみると、厳密に正規分布に従う必要はなく、ゼロ近傍の発生確率が高く、ゼロから離れるに従って発生確率が低くなればいいとのこと。試行錯誤して、sin(rand())など、乱数を生成した後に非線形の関数で変換してやればなんとかなりそうだ、というアイデアもすでに持っていたので、こんな方法もあるよと、二つの方法を紹介しました。（学生はC++のコードを書いていましたが、ここではMatlabコードで書いています。）
まずは一様乱数ふたつをかけ算するもの。これは-log(x)に従う乱数になります。（図の上側パネル参照）
N = 10000; x = rand(N,1) .* rand(N,1); hist(x, sqrt(N));
次に一様乱数ふたつを足し算するもの。rand()は[0, 1)の範囲の一様乱数なので、2個加算すると、1.0付近に山のある[0, 2)の範囲の乱数になります。そこから1.0を引くと、[-1, +1)の範囲の乱数になります。
N = 10000; x = rand(N,1) + rand(N,1) - 1; hist(x, sqrt(N));
足し算する回数を増やすと、徐々に正規乱数に近づいてゆきます。以下の例では12回加算して6を引いているので、[-6, +6)の乱数となります*1。（図の下側パネル参照）
N = 10000; x = rand(N,1) + rand(N,1) + rand(N,1) + rand(N,1) + rand(N,1) + rand(N,1) + ... rand(N,1) + rand(N,1) + rand(N,1) + rand(N,1) + rand(N,1) + rand(N,1) - 6; hist(x, sqrt(N));
このように一様乱数を足していって正規乱数に近づくのは、中心極限定理の応用例と言えます。
【追記】12個の一様乱数を加算する理由として、[0,1)の一様乱数の分散は1/12なので、標準正規分布の条件である分散1にするためには12個が必要だとのことでした。
*1:一様乱数からの正規乱数の生成については、『C言語による最新アルゴリズム事典
