JavaScriptでt分布をそれなりにコンパクトに実装してみた

JavaScript組み込みの機能だけでt分布の累積分布関数(CDF)、確率密度関数(PDF)とCDFの逆関数である分位点関数を実装してみたのでそのメモです。
ちなみにCDFとPDFを同時に求めるために双対数を利用しているので、双対数の簡単な実用例としてもお読みいただけるかもしれません。

##実装方針
自由度$\nu$のt分布のPDFは以下です。

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

CDFはこの関数を$t$について$-\infty$から積分したものになります。
ここで、$\Gamma(x)$はガンマ関数です。
これを素直に実装しようと思うとガンマ関数がネックなわけです。
JavaScript組み込みのMathにはガンマ関数がありません。

まあ、$\nu$は自由度なので自然数しかとらないとしてしまえば

$$\Gamma(x+1) = x\Gamma(x)$$

の関係と$\Gamma(1)=1$、$\Gamma(1/2)=\sqrt{\pi}$から単純な階乗計算でガンマ関数の値を求めることもできます。
しかしながらこの方針で実装するとすぐに以下のようなかなしみを背負います。

• JavaScriptのNumberは64bit浮動小数点なので階乗なんか素直に実装したらあっという間に$2^{53}-1$を超えてしまう。計算結果が信用できない
• PDFはいいとしてもCDFをどうするか。ここからさらに数値積分実装するとコード長くなりそう

どうしたものでしょうか。

###CDFの実装
そんな僕に光を与えてくれたのが奥村晴彦先生のC言語による最新アルゴリズム事典です。
1991年に初版発行で2018年にも改訂されているすごい本です。みんなも買おう！

こちらで紹介されている方法は$f(t)$の積分が$t^{2}/\nu=\tan^{2}\theta$と置けば$\cos ^{v-1}\theta$の積分になることと
$$
\int \cos ^{m} \theta d \theta=\frac{1}{m} \sin \theta \cos ^{m-1} \theta+\frac{m-1}{m} \int \cos ^{m-2} \theta d \theta
$$
を利用して$\cos\theta$の次数を下げていくと途中の係数とガンマ関数がうまく打ち消しあうというエレガントなものです。
結果として$\sin\theta,\cos\theta$をそれぞれ$s,c$とすると、t分布のCDFは$\nu$が奇数の時
$$
\int_{0}^{t} f(u) d u=\frac{1}{\pi}\left(\theta+s c+\frac{2}{3} s c^{3}+\cdots+\frac{2 \cdot 4 \ldots(v-3)}{3 \cdot 5 \ldots(v-2)} s c^{v-2}\right)
$$
偶数の時
$$
\int_{0}^{t} f(u) d u=\frac{s}{2}\left(1+\frac{1}{2} c^{2}+\cdots+\frac{1 \cdot 3 \ldots(v-3)}{2 \cdot 4 \ldots(v-2)} c^{v-2}\right)
$$
となります。すごい！

…って、あれ？
CDFはこれでいいとして、PDFはどうすれば…？

###PDFの実装
PDFはCDFの導関数なのでした。
というわけで、CDFの計算をしながら同時にその微分も求めてしまえばいいのです。
でも中心差分で数値微分するのではステップ幅を$h$としたときに$O(h^{2})$の誤差が乗ってしまいます。それを嫌って$h$を小さくしすぎれば今度は桁落ちの問題が…。

いいえ、追加的な誤差はほとんどなく微分値を求めることができます。双対数ならね。
双対数自体は数学的にとても難解なのですが、ここでは利用するためだけの大雑把な説明をします。(間違いがありましたらご指摘ください。数学つよつよの方々…)

双対数は複素数と同じように実部と非実部がある2次元の値です。
$$d = a+b\epsilon$$この$\epsilon$が複素数でいうところの$i$のようなものです。
$i^2=-1$だったわけですが双対数では$\epsilon^2=0$と定義します。

定義を守って計算すれば双対数同士の四則演算は以下になることが簡単にわかります。
〇足し算
$$(a_1+b_1\epsilon)+(a_2+b_2\epsilon) = (a_1+a_2)+(b_1+b_2)\epsilon$$
〇引き算
$$(a_1+b_1\epsilon)-(a_2+b_2\epsilon) = (a_1-a_2)+(b_1-b_2)\epsilon$$
〇掛け算
$$\left(a_{1}+b_{1} \epsilon\right) \cdot\left(a_{2}+b_{2} \epsilon\right)=a_{1} a_{2}+\left(a_{1} b_{2}+a_{2} b_{1}\right) \epsilon$$
〇割り算
$$\frac{\left(a_{2}+b_{2} \epsilon\right)}{\left(a_{1}+b_{1} \epsilon\right)} =
\frac{a_{2}}{a_{1}}+\frac{a_{1} b_{2}-a_{2} b_{1}}{a_{1}^{2}} \epsilon
$$
何かお気づきでしょうか？
掛け算と割り算の$\epsilon$の係数に注目してください。
この形、どこかで見覚えがありませんか？
実は積の微分、商の微分公式と同じ形になっています。
つまり、非実部の係数が実部の導関数になっていれば四則演算後もその関係が保たれるのです。

もうひとつ、双対数の便利な性質があります。
あるTaylor展開な関数$f$に対して以下が成り立ちます。
$$f(x+\epsilon)=f(x)+f'(x)\epsilon$$
なぜかというと、刻み幅を$\epsilon$とした$f$の$x$まわりのTaylor展開を考えると
$$f(x+\epsilon)=f(x)+f^{\prime}(x) \epsilon+\frac{1}{2 !} f^{\prime \prime}(x) \epsilon^{2}+\cdots$$
となり、定義から$\epsilon^2=0$なので右辺3項目以降は全て0になるからです。

なんとも不思議な話ですが、以上より双対数は

• 非実部が実部の導関数という関係が成り立っているもの同士なら四則演算した後も非実部が実部の導関数であり続ける
• Taylor展開可能な関数に双対数を代入すれば非実部がその関数の導関数になる

という非常に便利な性質を持っていることが分かります。

実際にこの双対数を用いて導関数を求めるには、双対数同士の四則演算と計算に用いる全ての関数の導関数をあらかじめ実装しておけばOKです。

というわけで双対数を実装してCDFの計算に用いれば、非実部がそのままPDFになるわけです。魔法みたい。
今回t分布のCDFの計算で用いるあらかじめ導関数を与えておく必要がある関数は平方根とarctanのみです。それぞれ以下のようになります。
$$(\sqrt{x})' = \frac{1}{2\sqrt{x}}$$$$(\arctan x)' = \frac{1}{1+x^{2}}$$
とっても簡単。

###分位点関数(CDFの逆関数)
最後に信頼区間計算などに必要な分位点関数を実装します。
できるだけ簡単に実装したいので、こちらを参考にしてNewton法で求めることにします。
Newton法は言わずと知れた微分可能な非線形関数の求根アルゴリズムです。
すなわち根を求めたい関数を$f$として
$$f(x)=0$$となる$x$を求めます。
今回のケースではt分布のCDFを$CDF$として$f$を以下のように考えればいいわけですね。
$$f(t) = CDF(t)-p$$
与えられた$p$に対して右辺を0にするような$t$を求めればいいので、
$$t^{(i+1)} = t^{(i)} - \frac{CDF(t^{(i)})-p}{PDF(t^{(i)})}$$
に従って$t$を適当な精度になるまで更新すればOKです。

##実装コード
これまでの説明をコードに落とし込んだのが以下です。
双対数はきちんと実装するならDual Numberクラスを用意すべきですが、今回は横着してArrayの0番目の要素を実部、1番目の要素を非実部として計算してます。

tdist.js

// t分布のCDFとPDF
function TDISTFUNC(t,df) {
  // 双対数の四則演算を定義
  const dadd = (d1,d2) => [d1[0]+d2[0], d1[1]+d2[1]];
  const dsub = (d1,d2) => [d1[0]-d2[0], d1[1]-d2[1]];
  const dmult = (d1,d2) => [d1[0]*d2[0], d1[0]*d2[1]+d1[1]*d2[0]];
  const ddiv = (d1,d2) => [d1[0]/d2[0], (d2[0]*d1[1]-d1[0]*d2[1])/d2[0]/d2[0]];
  // 双対数の平方根とarctanを定義
  const dsqrt = (d) => [Math.sqrt(d[0]), d[1]*0.5/Math.sqrt(d[0])];
  const darctan = (d) => [Math.atan(d[0]), d[1]/(1+d[0]*d[0])];
  // ここからt分布のCDFとPDFを求める計算
  const c2 = ddiv([df,0],dadd([df,0],dmult([t,1],[t,1])));
  let s = dmult([Math.sign(t),0],dsqrt(dsub([1,0],c2)));
  let p = [0,0];
  for(let i=df%2+2; i<=df; i+=2){
    p = dadd(p,s); 
    s = dmult(s,dmult([((i-1)/i),0],c2));
  }
  const theta = darctan(ddiv([t,1],dsqrt([df,0])));
  // 自由度が奇数か偶数かで計算結果を変える
  const result = df%2===0? 
    dmult([0.5,0],dadd([1,0],p)) : 
    dadd([0.5,0], ddiv(dadd(dmult(p,dsqrt(c2)),theta),[Math.PI,0]));
  return {'CDF':result[0] , 'PDF':result[1]};
}
// t分布の分位点関数(CDFの逆関数) <- t分布のCDFとPDF
function TQUANTILE(p,df,tol=1e-10,maxiter=1000) {
  let t, gain, iter;
  [t, gain, iter] = [1, Infinity, 0];
  // tをNewton法で更新
  while(Math.abs(gain)>tol && iter<=maxiter){
    gain = -(TDISTFUNC(t,df).CDF-p)/TDISTFUNC(t,df).PDF;
    t += gain;
    iter++;
  }
  return t;
}

以上！それなりにコンパクトに実装できたと思います。