トップ > LISP > 二項分布の確率質量関数

二項分布の確率質量関数

LISP

【2023-04-17】このページの最下部に載っているコードでは、特定の引数の時に間違った答えが出力されたり、無限ループに入ってしまうことが分かりました。たとえば(dbinom 4 10 0.5) ;=>0.557となってしまいます（正しくは0.205）。元のアルゴリズムには問題なく、僕が書いたLispコードが間違っていましたので、修正したものを「二項分布の確率質量関数（改） - 丸井綜研」にアップロードしました。

今日の午後に開催されたLispお茶会イベント（Shibuya.lisp lispteatime #11 - connpass）に参加しました。とくに事前にテーマを決めずにもくもく会に参加したのですが、作り途中になっていた二項分布の確率質量関数（Rでのdbinom）の高速化と精度改善をすることにしました。

二項分布の確率質量関数は以下のように定義されています。高校数学II/Bでも登場するので見たことがある人も多いのではないでしょうか。（高校数学だと(n k)ではなくnCkというような表記でした）

\displaystyle {n \choose k}p^{k}(1-p)^{n-k}

確率pで起こる事象があるとき、n回の反復試行したときに事象がk回起きる確率を計算するというものです。たとえば「確率1/2で表が出るコイントスを20回やったとき、16回表が出る確率は0.46パーセント」「立方体のサイコロを5回ふった時に2の目が3回出る確率は3.2パーセント」などと計算できます。

ナイーブな実装

これをナイーブに実装すると

(defun dbinom (x n &optional (p 1/2))
  (* (choose n x) (expt p x) (expt (- 1 p) (- n x))))

となるでしょうか。ここで呼ばれているchooseは階乗の組み合わせで計算することもできますが、以下のように乗法表示で計算すると効率が良いようです。

(defun choose (n k)
  (let ((res 1))
    (loop for i from 1 to k
          do (setq res (* res (/ (+ n 1 (- i)) i))))
    res))

使ってみると、以下のように正しい値が計算されます。

CL-USER> (dbinom 16 20 1/2)
4845/1048576
CL-USER> (float (dbinom 16 20 1/2))
0.004620552

pの値による問題点

pが0に近いときや1に近いときに有限桁数のコンピュータで計算すると、pxや(1-p)n-xの値が極端に0や1に近くなることで正確でない値が出てきそうです。SBCLとRとで実行したときの結果が以下です。

CL-USER> (dbinom 16 20 0.001)
0.0

> dbinom(16, 20, 0.001)
[1] 4.825649e-45

あまりに小さい確率なので実質的には問題にならない気もしますが、対数をとりながら計算することで精度が上がるのではないかと思いました。

(defun dbinom (x n &optional (p 1/2))
  (exp (+ (log (choose n x))
          (* (log p) x)
          (* (log (- 1 p)) (- n x)))))

実行してみると以下のように、Rでの計算結果に近くなっています。

CL-USER> (dbinom 16 20 0.001)
4.2038954e-45

大きなxとnのとき

上記の方法では、大きなxやnのときにはchooseの部分でのループ回数が大きくなりすぎてしまいます。もしガンマ関数あたりが使用できるのであれば、ガンマ関数の対数（gammaln）を使って近似的に計算ができます。gammalnがO(1)であれば、dbinomもO(1)で計算できるので高速です。ただ、計算の精度はガンマ関数の実装によって制限がかかります。（ガンマ関数を使ったchooseの計算については「ガンマ関数 - 丸井綜研」に書きました）

(defun dbinom (x n &optional (p 1/2))
  (exp (+ (gammaln (1+ n))
          (- (gammaln (1+ x)))
          (- (gammaln (- n x -1)))
          (* (log p) x)
          (* (log (- 1 p)) (- n x)))))

以下のように実行すると、上記のchooseを使用したものでは実用的な時間では計算が終わりませんが、ガンマ関数を使用したバージョンでは何とか計算できています。（とはいえ、下記のアルゴリズムとは計算結果が3倍違うというヤバいことが起きています）

CL-USER> (dbinom 1000000 2000000 1/2)
0.0020879126042240128d0

Rで採用されているアルゴリズム

高精度かつ高速なアルゴリズムを求めて、Rではどのようなアルゴリズムで実装しているのかをhelp(dbinom)で調べてみたら、以下の記事が参考文献に記載されていました。

Catherine Loader (2000). Fast and Accurate Computation of Binomial Probabilities; available as https://www.r-project.org/doc/reports/CLoader-dbinom-2002.pdf

この記事によると、ナイーブな実装をしてしまうとdbinom(1000000, 2000000, 1/2)というような大きなxとnのときには誤差が大きすぎて、7桁くらいの精度しか確保できなくなってしまうとのことです。このエントリでは詳しい記事の内容は省いて、C言語で書かれたコードをCommon Lispに移植した結果だけを掲載しておきます。

【2023-04-17】以下のコードでは、特定の引数を与えたときに問題が生じることが分かりました。たとえば(dbinom 4 10 0.5)の結果が0.557となってしまったり（正しくは0.205）、(dbinom 6 10 1/2)は無限ループに陥ります。元のアルゴリズムには問題なく、僕が書いたLispコードが間違っていましたので、修正したものを「二項分布の確率質量関数（改） - 丸井綜研」にアップロードしました。以下のコードは修正前のものをそのまま掲載しています。

(defun stirlerr (n)
  (let ((s0 (/ 1 12d0))
        (s1 (/ 1 360d0))
        (s2 (/ 1 1260d0))
        (s3 (/ 1 1680d0))
        (s4 (/ 1 1188d0))
        (sfe '(0.0d0
               0.081061466795327258219670264d0
               0.041340695955409294093822081d0
               0.0276779256849983391487892927d0
               0.020790672103765093111522771d0
               0.0166446911898211921631948653d0
               0.013876128823070747998745727d0
               0.0118967099458917700950557241d0
               0.010411265261972096497478567d0
               0.0092554621827127329177286366d0
               0.008330563433362871256469318d0
               0.0075736754879518407949720242d0
               0.006942840107209529865664152d0
               0.0064089941880042070684396310d0
               0.005951370112758847735624416d0
               0.0055547335519628013710386899d0))
        (nn (* n n)))
    (cond ((< n 16) (nth n sfe))
          ((> n 500) (/ (- s0 (/ s1 nn)) n))
          ((> n 80) (/ (- s0 (/ (- s1 (/ s2 nn)) nn)) n))
          ((> n 35) (/ (- s0 (/ (- s1 (/ (- s2 (/ s3 nn)) nn)) nn)) n))
          (t (/ (- s0 (/ (- s1 (/ (- s2 (/ (- s3 (/ s4 nn)) nn) nn) nn)) nn)) n)))))


(defun bd0 (x np)
  (if (< (abs (- x np)) (* 0.1 (+ x np)))
      (let* ((v (/ (- x np) (+ x np)))
             (s (* (- x np) v))
             (ej (* 2 x v))
             (s1 0.0d0))
        (loop with j = 1
              do (setq ej (* ej v v))
              (setq s1 (+ s (/ ej (+ (* 2 j) 1))))
              (when (= s1 s) (return s1))
              (setq s s1)))
      (* x (log (/ x np)))))

(defun dbinom (x n &optional (p 1/2))
  (cond ((= p 0.0) (if (= x 0) 1.0 0.0))
        ((= p 1.0) (if (= x n) 1.0 0.0))
        ((= x 0) (exp (* n (log (- 1 p)))))
        ((= x n) (exp (* n (log p))))
        (t (let ((lc (- (stirlerr n)
                        (stirlerr x)
                        (stirlerr (- n x))
                        (bd0 x (* n p))
                        (bd0 (- n x) (* n (- 1.0 p))))))
             (* (exp lc) (sqrt (/ n (* 2 pi x (- n x)))))))))

以下のように実行します。

CL-USER> (dbinom 1000000 2000000 1/2)
5.641895130240628d-4

同じことをRでやると以下のようになりました。

> dbinom(1000000, 2000000, 1/2)
[1] 0.0005641895

同じ計算アルゴリズムを採用しているので結果が同じになるのは当たり前ですが、なんとなく嬉しいですね。

もくもく会だったこともあって、普段以上に集中して作業をすることができました。次に参加できたら音に関係する何かができるといいなぁ。