JavaScriptにガンマ関数がないので、Lanczos近似で実装した話

Last updated at 2026-04-24Posted at 2026-04-24

この記事は、統計学習ツール StatPlay を開発する中で学んだことをまとめたものです。AIと協業して開発していますが、数学的な実装の正確性は自分なりに理解・検証しています。完璧かと言われると、まだ学びの途中です。

前回の記事で、統計学習ツール StatPlay を作った経緯を書きました。

前回の最後に「正確性をどう担保したかは次回以降で」と書きました。今回はその話の入り口です。

最初に言っておくと、私は数学が得意ではありません。高校数学の対数あたりから怪しくなって、統計の勉強中に何度も教科書に戻りました。そういう人間が、AIが書いたガンマ関数のコードを必死に理解しようとした話です。

t分布を描きたかっただけなのに

事の発端は、StatPlay の「仮説検定」セクションでt分布のグラフを描きたかったことでした。

StatPlay には正規分布のセクションがあって、そこでは normPDF（正規分布の確率密度関数）を使ってグラフを描いています。正規分布の式は比較的シンプルで、JavaScriptの標準的な数学関数だけで書けます。

// 正規分布の確率密度関数 — Math.exp, Math.sqrt, Math.PI だけで書ける
export function normPDF(x, mu = 0, sd = 1) {
  return Math.exp(-0.5 * ((x - mu) / sd) ** 2) / (sd * Math.sqrt(TAU));
}

ところが、t分布を描こうとしたときに事件が起きました。AIに「t分布のグラフも描いて」と指示したら、こういうコードが出てきました。

// t分布 — 自由度dfのt分布の確率密度関数
export function tPDF(x, df) {
  return Math.exp(lgamma((df + 1) / 2) - lgamma(df / 2))
    / Math.sqrt(df * Math.PI)
    * Math.pow(1 + x * x / df, -(df + 1) / 2);
}

動いて、グラフも出ました。でも後から「これ何やってるんだ」とコードを読み返したら、まったく意味がわからない。lgamma とか Math.exp とかが組み合わさった式が並んでいて、何を計算しているのか見当もつかない。

で、このコードが何を実装しているのか調べたら、t分布の確率密度関数の式にたどり着きました。

$$
f(x; \nu) = \frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\sqrt{\nu\pi},\Gamma\left(\frac{\nu}{2}\right)} \left(1 + \frac{x^2}{\nu}\right)^{-\frac{\nu+1}{2}}
$$

$\Gamma$ って何だ。コードにも式にも出てくるけど、自分が何を動かしているのかわからない。これはまずい。

ガンマ関数とは何か

AIに聞きました。「この式に出てくるΓって何ですか」と。

教えてもらった内容を自分なりに噛み砕くと、こういうことでした。

階乗（$n!$）は知っています。$5! = 5 \times 4 \times 3 \times 2 \times 1 = 120$。ただ、これは自然数でしか使えません。$3.5!$ とか $0.5!$ は定義できない。

ガンマ関数は、この階乗を実数（と複素数）に拡張したものです。自然数 $n$ に対して $\Gamma(n) = (n-1)!$ が成り立ちます。つまり $\Gamma(5) = 4! = 24$ です。

階乗（自然数のみ）	ガンマ関数（実数に対応）	関係
1! = 1	Γ(2) = 1! = 1	Γ(n) = (n−1)!
2! = 2	Γ(3) = 2! = 2	ひとつズレる
3! = 6	Γ(4) = 3! = 6
4! = 24	Γ(5) = 4! = 24

ガンマ関数は非整数にも対応: Γ(2.5), Γ(0.5) = √π なども計算できる

なんで $\Gamma(5) = 5!$ じゃなくて $(5-1)!$ なのかというと、歴史的な経緯でそうなっているらしいです。最初に知ったときは「なんでひとつズレてるんだ」と思いましたが、まあそういうものだと受け入れました。

で、t分布の式にガンマ関数が出てくるのは、自由度が整数じゃないケースにも対応するためです。$\Gamma(\nu/2)$ の $\nu$ が奇数だと、$\Gamma(2.5)$ とか $\Gamma(3.5)$ みたいな値が必要になる。階乗では計算できない。だからガンマ関数が要る。

ここまではわかりました。問題は、JavaScriptにはガンマ関数がない ということでした。

Python には math.lgamma() がある。R には lgamma() がある。でも JavaScript の Math オブジェクトにはない。もっとも、AIは「ない」なんて言いません。何も言わずに勝手に実装してきます。問題はそこに何が書いてあるかでした。

AIが書いてきたコードを読む

AIが書いてきたガンマ関数のコードには lgamma という関数名がついていて、中に「Lanczos近似」というコメントがありました。何それ。

調べてみると、Lanczos近似はガンマ関数を近似計算するためのアルゴリズムで、Python の math.lgamma や多くの数値計算ライブラリでも使われている方法らしい。選択肢として妥当だったようですが、私が選んだわけではありません。AIが勝手に選んでいます。

もうひとつわからなかったのが、なぜガンマ関数そのものではなく 対数ガンマ関数（$\ln\Gamma(z)$）を計算しているのか。これもAIに聞きました。ガンマ関数は値がすぐに巨大になるからだそうです。$\Gamma(100)$ はもう Infinity になってしまう。対数を取っておけば桁あふれを防げる。最終的に値が必要なときに Math.exp(lgamma(z)) で戻せばいい。なるほど。

z	gamma(z) 直接計算	lgamma(z) 対数を取る
5	24	≈ 3.18
20	≈ 1.2 × 10¹⁷	≈ 36.4
172	Infinity（破綻）	≈ 701.4

lgamma(z) = ln Γ(z) なら桁あふれしない。必要なときに Math.exp(lgamma(z)) で戻せる

で、これが実際のコードです。

// 対数ガンマ関数 — Lanczos近似
// 定数 g=7 と係数配列 c[] は近似精度を決めるパラメータ
export function lgamma(z) {
  const g = 7;
  const c = [
    0.99999999999980993,    // c[0]: 級数の初項
    676.5203681218851,      // c[1]〜c[8]: 近似の精度を上げる係数
    -1259.1392167224028,
    771.32342877765313,
    -176.61502916214059,
    12.507343278686905,
    -0.13857109526572012,
    9.9843695780195716e-6,
    1.5056327351493116e-7
  ];

  // z < 0.5 のとき：オイラーの反射公式で正の領域に変換する
  // Γ(z)Γ(1-z) = π / sin(πz) を使う
  if (z < 0.5) return Math.log(Math.PI / Math.sin(Math.PI * z)) - lgamma(1 - z);

  z -= 1;

  // Lanczos級数の計算
  let x = c[0];
  for (let i = 1; i < g + 2; i++) x += c[i] / (z + i);

  const t = z + g + 0.5;

  // 最終的な対数ガンマ値を返す
  return 0.5 * Math.log(TAU) + (z + 0.5) * Math.log(t) - t + Math.log(x);
}

// ガンマ関数本体 — 対数ガンマから復元する
export function gamma(z) {
  return Math.exp(lgamma(z));
}

コードの流れを整理するとこうなります。

Lanczos近似は z ≥ 0.5 で精度が出るよう設計されているため、z < 0.5 は 1−z に変換する

最初に読んだとき、係数配列 c の数値が何なのかさっぱりわかりませんでした。これまたAIに聞くと、Lanczos近似の理論から導出される値で、$g=7$ のときにこの9個の係数を使うと倍精度浮動小数点の範囲で十分な精度が出る、と。正直、導出までは追いきれていません。

理解に一番時間がかかったのは z < 0.5 のときの処理でした。これはオイラーの反射公式というもので、$\Gamma(z)\Gamma(1-z) = \pi / \sin(\pi z)$ という関係を使って、負の領域を正の領域の計算に変換しています。再帰呼び出しになっていて、lgamma(1 - z) を呼ぶと今度は引数が 0.5 以上になるので、本体の計算に入る。ここは理解できたときに「なるほど」と思いました。

……と書いていて気づいたんですが、Lanczos近似だのオイラーの反射公式だの、統計検定2級の試験範囲なんてとっくの昔に通り過ぎています。そもそもオイラーって誰ですか。統計の勉強をしていたはずなのに、いつの間にか18世紀の数学者の公式を実装している。局地的にではありますが、明らかにヤバい方向に進んでいます。

t分布・カイ二乗分布・F分布も、もう書いてあった

lgamma を理解し終わる頃には、AIはとっくにt分布・カイ二乗分布・F分布の確率密度関数も書き終えていました。

// t分布 — 自由度dfのt分布の確率密度関数
export function tPDF(x, df) {
  // lgammaのおかげで Γ((df+1)/2) / Γ(df/2) が計算できる
  return Math.exp(lgamma((df + 1) / 2) - lgamma(df / 2))
    / Math.sqrt(df * Math.PI)
    * Math.pow(1 + x * x / df, -(df + 1) / 2);
}

// カイ二乗分布 — 自由度dfのχ²分布の確率密度関数
export function chi2PDF(x, df) {
  if (x <= 0) return 0;
  return Math.exp(
    (df / 2 - 1) * Math.log(x) - x / 2 - (df / 2) * Math.log(2) - lgamma(df / 2)
  );
}

// F分布 — 自由度(d1, d2)のF分布の確率密度関数
export function fPDF(x, d1, d2) {
  if (x <= 0) return 0;
  // ベータ関数 B(d1/2, d2/2) = Γ(d1/2)Γ(d2/2) / Γ((d1+d2)/2)
  const num = d1 / 2 * Math.log(d1 * x) + d2 / 2 * Math.log(d2)
    - (d1 + d2) / 2 * Math.log(d1 * x + d2);
  const beta = lgamma(d1 / 2) + lgamma(d2 / 2) - lgamma((d1 + d2) / 2);
  return Math.exp(num - beta) / x;
}

読んでみると、どの関数も同じ構造でした。計算の途中で lgamma を使って対数の世界で足し引きして、最後に Math.exp で戻す。対数を取れば掛け算が足し算になるから、桁あふれしない。だから lgamma が先に要る。

こういう設計をAIは一瞬で書いてきます。自分で思いつくのは不可能です。

全部は追いつけない

なぜ対数を使うのか。なぜ再帰呼び出しが必要なのか。係数配列の数値は何なのか。ひとつずつAIに聞いて、対数の性質が怪しくなるたびに高校数学から教えてもらって、やっと追いつく。その繰り返しでした。

ただ、全部には追いつけませんでした。Lanczos近似の係数の導出は、正直いまでもわかっていません。反射公式の構造は理解できたけど、なぜそれで精度が保てるのかと聞かれると怪しい。

AIは一瞬で書く。自分は何時間もかけて読み解く。追いつけるところと追いつけないところがある。改めて人間——特に自分の身の丈を知りました。

でも、背後にある数理の世界が可視化される過程は、感じ取ることができました。
統計って、実はめっちゃ面白い。