最小実装を読む — vanilla JavaScript で Poisson を書く

式をそのまま関数へ落とす

このコースの JavaScript は、数式を隠さずにそのまま関数へ移した最小実装です。外部ライブラリは使わず、前提を満たさない入力は例外を送出します。

最小実装の JavaScript

ポアソン分布の計算とシミュレータで実際に使っている最小実装は次のとおりです。

function assertFiniteNonNegative(name, value) {
  if (!Number.isFinite(value) || value < 0) {
    throw new RangeError(`${name} must be a finite non-negative number`);
  }
}

function assertNonNegativeInteger(name, value) {
  if (!Number.isInteger(value) || value < 0) {
    throw new RangeError(`${name} must be a non-negative integer`);
  }
}

function factorialInt(n) {
  assertNonNegativeInteger("n", n);
  let acc = 1;
  for (let i = 2; i <= n; i += 1) {
    acc *= i;
  }
  return acc;
}

function poissonPmf(lambda, k) {
  assertFiniteNonNegative("lambda", lambda);
  assertNonNegativeInteger("k", k);
  return Math.exp(-lambda) * (lambda ** k) / factorialInt(k);
}

function poissonCdf(lambda, k) {
  assertFiniteNonNegative("lambda", lambda);
  assertNonNegativeInteger("k", k);
  let sum = 0;
  for (let i = 0; i <= k; i += 1) {
    sum += poissonPmf(lambda, i);
  }
  return sum;
}

function samplePoissonKnuth(lambda, rng = Math.random) {
  assertFiniteNonNegative("lambda", lambda);
  if (lambda === 0) return 0;

  const limit = Math.exp(-lambda);
  let product = 1;
  let count = 0;

  while (product > limit) {
    product *= rng();
    count += 1;
  }
  return count - 1;
}

poissonPmf は 1 本の棒、poissonCdf は左から順に足した累積、samplePoissonKnuth は乱数から 1 サンプル引く関数です。

コードと式の対応

入門では、まず poissonPmf と poissonCdf が読めれば十分です。サンプラは、あとでシミュレーションやモンテカルロへ進むときの入口になります。

Knuth サンプラはなぜ動くのか

samplePoissonKnuth は、ポアソン分布から乱数を 1 つ引くための古典的アルゴリズム（Donald Knuth, 1969）です。コードだけ見るとブラックボックスに感じますが、その背景には「ポアソン過程の到着間隔は指数分布に従う」という事実があります。

• 発生率 λ（単位区間あたり）のポアソン過程では、連続するイベントの間隔は Exponential(λ) に従う。
• 区間 [0, 1] に到着するイベント数 X は Poisson(λ) に従う。
• つまり「指数分布の和を区間 1 まで足し続けて何件入るか」を数えれば、Poisson 乱数が得られます。

指数乱数 E_i は、一様乱数 U_i ∈ (0, 1] から E_i = -ln(U_i) / λ で作れます。「E_1 + E_2 + ... + E_n ≤ 1 を満たす最大の n」が Poisson(λ) 乱数です。両辺を -λ 倍して exp を取ると次のように変形できます。

この同値変形がアルゴリズムの肝です。コードでは limit = e^(−λ) を最初に固定し、product に一様乱数を掛け続けて、product が limit を下回るまでのループ回数を数えます。具体的には次の対応です。

• product: それまでに引いた一様乱数の積 U_1 × U_2 × ... × U_count。
• while (product > limit): 積がまだ e^(−λ) を上回っている、つまり指数乱数の累積和が 1 を超えていない間ループを続ける。
• count - 1: 「ループを抜ける直前まで条件を満たしていた回数」を返す。count は product *= rng() の直後に + 1 しているので、最後の 1 回（条件を破った乱数）を取り除く必要があり、戻り値は count - 1。

lambda === 0 のときに早期 return しているのは、e^0 = 1 となって最初の product = 1 がそのまま比較され、ループが回らない（または微妙な丸めで不安定になる）ことを避けるためです。

このアルゴリズムは λ が小さいときに直感的で実装も短いという長所がありますが、λ が大きい（おおむね 30 以上）と平均で λ 回程度の乗算が必要になり遅くなります。本格的な用途では PTRS や PA など、大きい λ 向けのアルゴリズムを使い分けるのが普通です。

実装上の注意点

この実装はあくまで教育用の最小コードです。実用に向けて発展させるときに留意したい点を 4 つ挙げておきます。

• lambda ** k と ES2016: 指数演算子 ** は ES2016 で導入されました。Node.js や近年のブラウザでは問題ありませんが、古い環境を意識する場合は Math.pow(lambda, k) の方が互換性が高くなります。本講座では数式 λ^k との見た目の対応を優先して ** を採用しています。
• factorialInt のオーバーフロー: factorialInt は素直なループ実装です。JavaScript の Number は 53 bit 精度で、22! あたりから誤差が入り、171! 以上は Infinity になります。λ が大きく k も大きいケースを扱うなら、P(X = k+1) = P(X = k) × λ / (k+1) の漸化式で計算する、対数空間で計算する、などの工夫が必要です。本実装は λ が小さい教材用途を想定しています。
• Knuth サンプラの戻り値 count - 1: 上の節で説明したとおり、ループを抜けた時点では「条件を破った 1 回ぶん」が count に含まれているため、count - 1 を返すことで Poisson(λ) サンプルになります。count から 0, 1, 2, ... と数え始めれば自然なオフセットです。
• 例外の種類: RangeError: lambda < 0 や k が非整数のような「型は正しいが値が範囲外」というケースは RangeError が適切です。TypeError は引数の型自体が間違っているとき（例: 数値ではなく文字列が渡された）に使うのが慣習です。教材としても、暗黙に丸めたり NaN を返したりするより、明示的に例外を中断させたほうが原因を追いやすくなります。

第 6 章のまとめ

• Poisson の pmf / cdf は、数式をそのまま JavaScript へ写せる。
• Knuth サンプラは「指数分布の累積和が 1 を超えるまでの件数」を一様乱数の積で書き換えたもの。
• 不正な入力は RangeError で例外を送出すると、教材としても実装としても追いやすい。
• シミュレータの数字は、この最小実装と同じ計算で出せる。