二項分布の累積分布関数（と正則化不完全ベータ関数）

統計に使える関数の実装を続けています。2月のLISP もくもく会では二項分布の確率質量関数（dbinom）を実装しました。今度は累積分布関数（pbinom）が欲しいところです。

marui.hatenablog.com

原理どおりに動くのであれば確率質量関数を累積加算すると累積分布関数が作れるのですが、計算機内で累積計算していくと誤差も累積されて正しい答えから外れてしまうというのが普通です。試しに下記のようにp=1/2の条件下で10試行のうち8回以下の回数で成功する確率（と20万試行のうち10万回以下の回数で成功する確率）を計算してみましょう。

(defun pbinom (x n &optional (p0 1/2))
  (loop for i from 0 to x
        sum (dbinom i n p0)))

(pbinom 8 10 0.5)
;;=> 0.9892578125d0

(pbinom 100000 200000 1/2)
;;=> 0.5008920609429977d0

Rのpbinom(8, 10, 1/2)で計算した結果は0.9892578、pbinom(100000, 200000, 1/2)で計算した結果は0.5008921なので、驚いたことに、なぜか一致しています。すごい。でも、やはりループを使って加算していくのが、無駄にCPUを使っている気がします。だって、後者のときは10万回もループを回しているんです。十分に精度が高ければ近似値でもいいので、ズバッと答えが出せるのが現実的な理想です。

正則化不完全ベータ関数

Wikipediaの二項分布のページ（二項分布 - Wikipedia）を見ると、累積（積分）するのと同等のことが正則化不完全ベータ関数 I_x (a,b) で計算できると書かれていました。正則化不完全ベータ関数は不完全ベータ関数とベータ関数の比

$\displaystyle I_x\left(a, b\right) = \frac{B_x\left(a, b\right)}{B\left(a, b\right)}$

とのことですが、それを近似計算する方法がDLMF (Digital Library of Mathematical Functions)に書かれていました。（DLMFはAbramowitz & Stegunの後継のようです）

dlmf.nist.gov

不完全ベータ関数はベータ分布の累積分布関数なので、定義どおりに計算しようとすると積分が必要になります。積分は解析的にも数値的にも難しいので可能であれば扱いたくありません。セクション8.17(v)に記載されている正則化不完全ベータ関数を連分数を用いて近似してゆく式のほうがプログラムにしやすいので、いっちょやってみました。下記式のような無限に続く連分数になっています。

$\displaystyle I_{x}\left(a,b\right)=\frac{x^{a}(1-x)^{b}}{a\mathrm{B}\left(a,b\right)}\left(\cfrac{1}{1+\cfrac{d_{1}}{1+\cfrac{d_{2}}{1+\cfrac{d_{3}}{1+}}}}\cdots\right)$

ここでdは

$\displaystyle d_{2m}=\frac{m(b-m)x}{(a+2m-1)(a+2m)}$

$\displaystyle d_{2m+1}=-\frac{(a+m)(a+b+m)x}{(a+2m)(a+2m+1)}$

となります。連分数の形によっては最適化できる部分もあるのかもしれませんが、よく分からないのでチカラワザでやってみます。

;;; naive implementation of regularized incomplete beta function
;;; using continued fraction
;;; reference: https://dlmf.nist.gov/8.17#v
(defun ribeta (a b x)
  (labels ((d (a b x m)
             (cond ((= m 0) 1)
                   ((evenp m) (/ (* (/ m 2) (- b (/ m 2)) x)
                                 (* (+ a m -1) (+ a m))))
                   (t (/ (* -1 (+ a (/ (1- m) 2)) (+ a b (/ (1- m) 2)) x)
                         (* (+ a (1- m)) (+ a m))))))
           (contfrac (a b x m)
             (if (>= m (min (* a b) 100))
                 0
                 (let ((dd (d a b x m)))
                   (if (= dd 0)
                       0
                       (/ dd (1+ (contfrac a b x (1+ m))))))))
           (ribeta-sub (a b x)
             (cond ((<= x 0) 0)
                   ((>= x 1) 1)
                   (t (* (/ (* (expt x a) (expt (- 1 x) b))
                            (* a (beta a b)))
                         (contfrac a b x 0))))))
    (if (< x (/ (1+ a) (+ a b 2)))
        (ribeta-sub a b x)
        (- 1 (ribeta-sub b a (- 1 x))))))

labelsで定義した3つの関数内関数があります。d関数は、連分数の分子にくる値を計算するもので、連分数の奇数番目と偶数番目とで計算式が異なります。DLMFでは2mと2m+1というように場合分けされていましたが、上記のコードではmが偶数か奇数かで分けていますので、DLMFの式とはmの扱いが異なります。

contfrac関数はdをまとめるためのもので、再帰を使って連分数の内部に入っていくようになっています（連分数は英語でcontinued fractionなので、そこから名前を取りました）。無限に続く連分数なので、値が収束したら打ち切らないといけません。再帰から帰ってくるまで連分数全体の値が計算できないため、収束判定をどうすれば良いのかが思いつきません。どのくらいの深さまで連分数を再帰で潜ればいいのか分からなかったので、a×bか100のどちらか小さいほうの深さまでで打ち止めにするようにしています。分子がすべて1の連分数（正則連分数）については英語版Wikipediaに収束に関連する内容が記載されていました（Continued fraction - Wikipedia）が、今回は分子がころころ変化するのでそのままは使えません。先日の「二項分布の確率質量関数（改） - 丸井綜研」のように、浮動小数点数の精度を利用した繰り返しができればスマートなのかもしれませんが、どうやればいいのでしょう？　dの結果がゼロになったら終了というのはパッと思いつきましたが……。（後日、連分数をマクロで生成して解くという方法を知りました。参考にして書き直したい……）

xの値によって収束速度が変わるとのことなので、ribeta-subでは条件分岐によって最終的な計算をしています。

この正則化不完全ベータ関数はベータ分布の累積分布関数と同じだということなので、Rのpbetaを使って検証してみます。とりあえず、以下の条件ではうまくいっているようです。

> pbeta(seq(0, 1, length.out=11), 3, 7)
 [1] 0.00000000 0.05297214 0.26180250 0.53716883 0.76821299 0.91015625
 [7] 0.97496525 0.99570911 0.99968614 0.99999700 1.00000000

CL-USER> (mapcar #'(lambda (x) (ribeta 3 7 x)) (linspace 0 1 11))
(0 0.052972137999999655d0 0.26180249599999833d0 0.5371688339999965d0
 0.7682129920000015d0 0.9101562500000006d0 0.9749652480000002d0
 0.9957091060000001d0 0.999686144d0 0.999997002d0 1)

というわけで、二項分布の累積分布関数は以下のようになりました。いまのところdbinomをループで累積したものと違う値は出ませんが、どこかに穴があるかもしれないので注意しつつ使っていこうと思います。

(defun pbinom (x n &optional (p0 1/2))
  (ribeta (- n x) (+ 1 x) p0))

おまけ

さて、（正則化されていない）不完全ベータ関数が欲しい場合には変換してあげれば良さそうですね。

;;; (non-regularized) incomplete beta function
(defun ibeta (a b x)
  (* (ribeta a b x) (beta a b)))

以下の論文には、不完全ベータ関数のFORTRAN実装の擬似コードが乗っていましたので、それを参考にして作り直してもいいかもしれません。

Armido R. Didonato and Alfred H. Morris. “Algorithm 708: Significant digit computation of the incomplete beta function ratios.” ACM Transactions on Mathematical Software, vol. 18, iss. 3, pp. 360–373. Sept. 1992. https://doi.org/10.1145/131766.131776