Box-Muller法と正規乱数生成器

PRML上巻の序章を読んだ頃に、ガウス分布に肄う乱数発生器を作りたいと思っていたので作ってみた。つまり、平均値と標準偏差を引数に取る乱数生成関数である。Web検索するとBox-Muller法を使うのと、中心極限定理を使う方法があるのでそれぞれ調べてみた。

Box-Muller法を使う

Box-Muller法は次の式によって作られる2つの独立変数を使うらしい。

$\sqrt{-2 * log(α)} * sin(2πβ)$
$\sqrt{-2 * log(α)} * cos(2πβ)$

式だけ見てもイメージが湧かないので、Octaveを使ってサンプル数2000でテストしてみた。

check_box_muller.matlab

clear;
alpha = rand(1, 1000);
beta = rand(1, 1000);
val1 = sqrt(-2 * log(alpha)) .* sin(2 * pi * beta);
val2 = sqrt(-2 * log(alpha)) .* cos(2 * pi * beta);
hist([val1,val2], 20);

ちゃんと正規分布っぽい。Box-Muller法で生成した標準正規分布に標準偏差を掛けて、平均値を足せば完成。

"use strict";

/**
 * @param {Number} mean 平均値
 * @param {Number} std 標準偏差
 */
function gaussianRandom(mean, std) {
  return randByBoxMullerTransform() * std + mean;
}

// Box Muller法を使って標準正規分布N(0,1)を得る
var randByBoxMullerTransform = (function() {
  var vals = [];

  function calc() {
    var alpha = Math.random(),
        beta = Math.random();
    return [
      Math.sqrt(-2 * Math.log(alpha)) * Math.sin(2 * Math.PI * beta),
      Math.sqrt(-2 * Math.log(alpha)) * Math.cos(2 * Math.PI * beta)
    ];
  }

  return function() {
    vals = vals.length == 0 ? calc() : vals;
    return vals.pop();
  }
})();

// Module export
module.exports = gaussianRandom;

中心極限定理を使う

中心極限定理を利用すると、複数の一様乱数の和を使うだけで標準正規分布N(0,1)が得られる。

"use strict";

/**
 * @param {Number} mean 平均
 * @param {Number} std 標準偏差
 */
function gaussianRandom(mean, std) {
  return randByCentralLimitTheorem() * std + mean;
}
 
// 中心極限定理を使って標準正規分布N(0,1)を得る
function randByCentralLimitTheorem() {
  var v = 0;
  for (var i = 0; i < 12; i++) {
    v += Math.random();
  }
  return v - 6;
}
 
// Module export
module.exports = gaussianRandom;

こっちの方が圧倒的に簡単でびびる。