二項分布の確率質量関数（改）

以下のエントリに、二項分布の確率質量関数をCommon Lispで計算する関数を紹介したのですが、どうも意図した通りに動いてくれないときがありました。少し時間を作って調べてみたら、コードに間違いがあるのが明らかになりました。

marui.hatenablog.com

Loaderによるアルゴリズムには問題はなく、僕がC言語のコードをLispに書き換えた時点でミスしていました。1箇所だけでなく複数箇所を間違えていたのですが、単純なミスだったり型によって正常に動かなかったりするものでした。大きな原因は、カッコの対応が見えていなかったことだと思います。

問題1

まず一つ目はstirlerr関数。

(defun stirlerr (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)))))

となっていたところですが、最終行の変なところにnnが入ってしまっていました。正しくは以下のようになります。

(defun stirlerr (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)) n)))))

問題2

次の問題点はbd0関数で、最終行の(* x (log (/ x np)))がif文のelse節になっていますが、本当はifの外になければなりませんでした。しかも、ループ部分でループ変数jを更新していないので、いつまでもループから外に出られません。

(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)))))

また、npの値が特定の整数や分数だと無限ループに陥ってしまうこともありましたので、冒頭部分でnpを浮動小数点型にcoerceすることにしました。loopの中で、s1とsが同値になるという浮動小数点数の精度を使った収束判定をしているようなので、分数型だと収束しなくなってしまうのでした。

(defun bd0 (x np)
  (setf np (coerce np 'long-float))
  (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)
                 (setq j (1+ j)))))
  (* x (log (/ x np))))

修正後のコード

以上の修正により、おそらく問題なく動いてくれるのではないかと思います。以下にコード全体を掲載します。

;; C. Loader (2002) Fast and Accurate Computation of Binomial Probabilities
;; https://www.r-project.org/doc/reports/CLoader-dbinom-2002.pdf
(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)) n)))))

(defun bd0 (x np)
  (setf np (coerce np 'long-float))
  (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)
                 (setq j (1+ j)))))
  (* x (log (/ x np))))

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

試しにRの結果と比較しましたが、良いようです。

> dbinom(0:10, 10, 1/2)
 [1] 0.0009765625 0.0097656250 0.0439453125 0.1171875000 0.2050781250
 [6] 0.2460937500 0.2050781250 0.1171875000 0.0439453125 0.0097656250
[11] 0.0009765625

CL-USER> (mapcar #'(lambda (x)
                     (dbinom x 10 1/2))
                 (linspace 0 10 11))
(9.765625e-4 0.009765625d0 0.04394531249999997d0 0.11718750000000003d0
 0.20507812499999997d0 0.24609375d0 0.20507812499999997d0 0.11718750000000006d0
 0.04394531249999997d0 0.009765625000000003d0 9.765625e-4)