Octaveでの実装を参考にCommon Lispで正規分布の累積分布関数を計算

正規分布確率密度関数は定義通りに

;; Normal Distribution
;; (normpdf 2 0 1) => 0.05399096581690089
;;
(defun normpdf (x mu sigma)
  (* (/ 1 (* (sqrt (* 2 pi)) sigma))
     (exp (- (/ (expt (- x mu) 2)
		(* 2 sigma sigma))))))

で計算できるのですが、累積分布関数は確率密度関数を定積分したものなので、やっかいでなかなかプログラムが作れずにいました。

Octaveではどういう実装になっているのかと、normcdfのMファイルを読み辿っていったら、誤差関数(erf)を使って計算していることが分かりました。ただし、このerfOctaveの内部実装なので、この中身を知るためには利用している数値計算ライブラリのソースを読むことになります。

そこでWikipediaで誤差関数を調べたテイラー展開による近似式が載っていたので、それを使って実装することにしました。Octaveでの1:50の作り方がよくわからなかったので、その部分はM.Hiroi氏のiotaを拝借しました。*1

;; Error Function
;; (erf 1) ==> 0.8427007929497148
;;
(defun iota (m &optional (n 1) (step 1))
  (if (> n m)
      nil
      (cons n (iota m (+ n step) step))))
(defun erf (x)
  (* (/ 2 (sqrt pi))
     (let ((sum 0))
       (dotimes (n 100 sum)
	 (setq sum (+ sum (* (/ x (+ (* 2 n) 1))
			     (prod (mapcar
				    (lambda (k) (/ (- (* x x)) k))
				    (iota n))))))))))
(defun erfc (x)
  (- 1 (erf x)))

dotimes部分ではテイラー展開を第100項まで計算*2しています。erfcは累積分布関数を計算するときに使う相補誤差関数です。そして、以下が累積分布関数です。

;; Cumulative Normal Distribution
;; (normcdf 2 0 1) => 0.05399096581690089
;;
(defun normcdf1 (x)
  (/ (erfc (/ x (- 0 (sqrt 2)))) 2))
(defun normcdf (x mu sigma)
  (normcdf1 (/ (- x mu) sigma)))

(normcdf 1.96 0 1)とやると、ちゃんと0.975が戻ってきます。

Octaveはもちろん、C言語Rubyにも誤差関数が標準で用意されているので、もしかしたらCommon Lispにもあるかもしれませんが・・・。車輪の再発明

*1:Octaveでは1:10で[1 2 3 4 5 6 7 8 9 10]というベクトルが作成されます。R言語にはseqClojureにはrangeがあるので、Common Lispにも同様な関数が用意されているのでしょうが、見つかりませんでした。

*2:本当はテイラー展開で第100項まで求める必要ないんです。自分が欲しい誤差に収まったらそこで計算を打ち切った方がよほど計算時間は短くなりそうです。