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

Apr 10, 2019 / Mar 2, 2026 · 4 分で読了 · 統計学 javascript 正規分布 ·

共有する:

概要

Javascriptで正規分布の乱数発生（rnorm）、確率密度関数（dnorm）、累積分布関数（pnorm）、累積分布の逆関数（qnorm）を実装する。すべて標準正規分布を想定。 Javascriptに限らず使えるアルゴリズムだが、日本語でまとまっている情報があまりないのと、ブラウザ上でA/Bテストなど有意性をみる検定などできたら面白いということでJSでやってみる。

なお、実務で手軽に使いたい場合は stdlib-js や jStat といったライブラリも検討するとよい。本記事はアルゴリズムの中身を理解する目的で、ライブラリを使わずスクラッチで実装する。

正規乱数の生成（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());
}

javascript

後ろの係数は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);
}

javascript

累積分布関数（pnorm）

Abramowitz and Stegun, Handbook of Mathematical Functions (1964)から。 https://personal.math.ubc.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;
}

javascript

これと同じ https://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/ 別の近似式を使って実装したものの例 https://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;
}

javascript

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

7.1.25 と 26.2.16
7.1.26 と 26.2.17
7.1.27 と 26.2.18
7.1.28 と 26.2.19

がそれぞれ対応している。

累積分布関数の逆関数（qnorm）

Peter J. Acklamのアルゴリズムを用いる。中央領域と裾領域で別々の有理ミニマックス近似を適用するもので、相対誤差の絶対値が全域で 1.15 × 10⁻⁹ 未満という高精度な近似である。

アルゴリズムの原典は以下（元サイトは閉鎖されているためWayback Machineのアーカイブ）。 https://web.archive.org/web/20151030215612/http://home.online.no/~pjacklam/notes/invnorm/

アルゴリズムの解説記事: https://stackedboxes.org/2017/05/01/acklams-normal-quantile-function/

Javascriptで実装

// Peter J. Acklam's algorithm for the inverse normal CDF
// precision: abs(relative err) < 1.15e-9

function qnorm(p) {

    // coefficients for rational approximation
    var a1 = -3.969683028665376e+01;
    var a2 =  2.209460984245205e+02;
    var a3 = -2.759285104469687e+02;
    var a4 =  1.383577518672690e+02;
    var a5 = -3.066479806614716e+01;
    var a6 =  2.506628277459239e+00;

    var b1 = -5.447609879822406e+01;
    var b2 =  1.615858368580409e+02;
    var b3 = -1.556989798598866e+02;
    var b4 =  6.680131188771972e+01;
    var b5 = -1.328068155288572e+01;

    var c1 = -7.784894002430293e-03;
    var c2 = -3.223964580411365e-01;
    var c3 = -2.400758277161838e+00;
    var c4 = -2.549732539343734e+00;
    var c5 =  4.374664141464968e+00;
    var c6 =  2.938163982698783e+00;

    var d1 =  7.784695709041462e-03;
    var d2 =  3.224671290700398e-01;
    var d3 =  2.445134137142996e+00;
    var d4 =  3.754408661907416e+00;

    // break-points
    var pLow  = 0.02425;
    var pHigh = 1 - pLow;

    var q, r;

    if (p < pLow) {
        // rational approximation for lower region
        q = Math.sqrt(-2 * Math.log(p));
        return (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) /
               ((((d1 * q + d2) * q + d3) * q + d4) * q + 1);
    } else if (p <= pHigh) {
        // rational approximation for central region
        q = p - 0.5;
        r = q * q;
        return (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q /
               (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1);
    } else {
        // rational approximation for upper region
        q = Math.sqrt(-2 * Math.log(1 - p));
        return -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) /
                ((((d1 * q + d2) * q + d3) * q + d4) * q + 1);
    }
}

javascript

入力 p は 0 < p < 1 の範囲で有効。qnorm(0.975) が約 1.96 を返すことで動作確認できる。