Javascriptで正規分布の実装まとめ（乱数、累積分布関数など）

Javascriptで正規分布の乱数発生（rnorm）、確率密度関数（dnorm）、累積分布関数（pnorm）、累積分布の逆関数（qnorm）を実装する（逆関数は参照で）。すべて標準正規分布を想定。

Javascriptに限らず使えるアルゴリズムだが、日本語でまとまっている情報があまりないのと、ブラウザ上でA/Bテストなど有意性をみる検定などできたら面白いということでJSでやってみる。

正規乱数の生成（rnorm）

1行でBox-Muller法で。

Box-Muller法とは？

$$X_1, X_2 \stackrel{i.i.d.}{\sim} {\rm Unif} (0, 1) $$

とするとき

$$Y_1 = \sqrt{-2 \log{X_1}} \cos{2 \pi X_2} $$
$$Y_2 = \sqrt{-2 \log{X_1}} \sin{2 \pi X_2} $$

で生成される

$$Y_1, Y_2 \stackrel{i.i.d.}{\sim} {\rm N} (0, 1) $$

というもの。

今回は1個の正規乱数でいいので、$Y_1$か$Y_2$の一方を採用すればいい。

Javascriptで実装

function rnorm(){
    return Math.sqrt(-2 * Math.log(1 - Math.random())) * Math.cos(2 * Math.PI * Math.random());
}

後ろの係数はMath.cos()でもMath.sin()でもどちらでもいい。JavascriptのMath.random()は戻り値の区間が[0,1)なので、$log {0}$で発散しないように1-Math.random()としている。

確率密度関数（dnorm）

$$Z(x) = \frac{ e^{ -\frac{x^2}{2}} }{\sqrt{2 \pi}} $$

そのまんま

function dnorm(x){
    return Math.exp(-x * x / 2) / Math.sqrt(2 * Math.PI);
}

累積分布関数（pnorm）

Abramowitz and Stegun, Handbook of Mathematical Functions (1964)から。
http://people.math.sfu.ca/~cbm/aands/

26.2が正規分布の累積分布関数の項目。
実際はC. Hastings, Jr., Approximations for Digital Computers (1955)に基づいているとのこと。

26.2.17の

$$P(x) = 1 – Z(x) \left( b_1 t + b_2 t^2 + b_3 t^3 + b_4 t^4 + b_5 t^5 \right) + \epsilon(x) $$
$$t = \frac{1}{1+px}, \quad Z(x) = \frac{ e^{ -\frac{x^2}{2}} }{\sqrt{2 \pi}} $$
$$|\epsilon(x)| \lt 7.5 \times 10^{-8} $$
$$p = .23164 19 $$
$$b_1 = .31938 1530 $$
$$b_2 = -.35656 3782 $$
$$b_3 = 1.78147 7937 $$
$$b_4 = -1.82125 5978 $$
$$b_5 = 1.33027 4429 $$

をそのまま実装

// Abramowitz and Stegun (1964) formula 26.2.17
// precision: abs(err) < 7.5e-8

function cdf(x) {

    // constants
    var p  =  0.2316419;
    var b1 =  0.31938153;
    var b2 = -0.356563782;
    var b3 =  1.781477937;
    var b4 = -1.821255978;
    var b5 =  1.330274429;

    var t = 1 / (1 + p * Math.abs(x));
    var Z = Math.exp(-x * x / 2) / Math.sqrt(2 * Math.PI);
    var y = 1 - Z * ((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t;

    return (x > 0) ? y : 1 - y;
}

これと同じ
http://math.stackexchange.com/questions/888165/abramowitz-and-stegun-approximation-for-cumulative-normal-distribution

参考
https://malishoaib.wordpress.com/2014/04/02/python-code-and-normal-distribution-writing-cdf-from-scratch/

別の近似式を使って実装したものの例
http://www.johndcook.com/blog/normal_cdf_inverse/

参考までに前著の7.1が誤差関数の項目で、こちらを使う場合もある。
cdf()が累積分布関数で、erf()が誤差関数。

// Abramowitz and Stegun (1964) formula 7.1.26
// precision: abs(err) < 1.5e-7

function cdf(x) {
    return (1 + erf( x / Math.sqrt(2) )) / 2;
}

function erf(x) {

    // constants
    var p  =  0.3275911;
    var a1 =  0.254829592;
    var a2 = -0.284496736;
    var a3 =  1.421413741;
    var a4 = -1.453152027;
    var a5 =  1.061405429;

    var t = 1 /(1 + p * Math.abs(x));
    var y = 1 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.exp(-x * x);

    return (x > 0) ? y : -y;
}

誤差関数と正規分布の関係は
http://www.fbs.osaka-u.ac.jp/labs/ishijima/error-01.html
http://bio-info.biz/statistics/distribution_normal_distribution.html
あたりが参考になる。

前著の数式では