(2023-01-23追記:正規分布の計算を見直してみました →正規分布の確率密度関数と累積分布関数の計算 - 丸井綜研)
;; 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
)を使って計算していることが分かりました。ただし、このerf
はOctaveの内部実装なので、この中身を知るためには利用している数値計算ライブラリのソースを読むことになります。
そこでWikipediaで誤差関数を調べたらテイラー展開による近似式が載っていたので、それを使って実装することにしました。Octaveでの1:50
の作り方がよくわからなかったので、その部分はM.Hiroi氏のiotaを拝借しました。((Octaveでは1:10で[1 2 3 4 5 6 7 8 9 10]というベクトルが作成されます。R言語にはseq
、Clojureにはrange
があるので、Common Lispにも同様な関数が用意されているのでしょうが、見つかりませんでした。))
;; 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項まで計算*1しています。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にもあるかもしれませんが・・・。車輪の再発明。