正規分布の確率密度関数と累積分布関数の計算

目次

Common Lispを音だけでなく、統計にも使いたいです。業務上の統計分析についてはRやJuliaを使っていて、Common Lispで置き換える必要も(今のところは)ないので、みんな大好き車輪の再発明をしていきます。まずは正規分布の計算から。

確率分布として欲しいのは、確率密度関数、累積分布関数、累積分布関数の逆関数、の三つです。あとは分布に従った乱数列もあるといいです。正規分布に従った乱数の生成については以前書いたので、今回は他の三つを作ります。

正規分布確率密度関数

標準正規分布N(0, 1)は以下の式です。

 \displaystyle
Z(x) = \frac{1}{\sqrt{2\pi}}\exp\left(\frac{-x^2}{2}\right)

平均値μと標準偏差σも指定できると便利なので、こっちの式を使います。

 \displaystyle
Z(x) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(\frac{-\left(x-\mu\right)^2}{2\sigma^2}\right)

これは素直にプログラムにすればいいです。

(defun dnorm (x &optional (mu 0) (sigma 1))
  (* (/ 1 (* (sqrt (* 2 pi)) sigma))
     (exp (- (/ (expt (- x mu) 2)
        (* 2 sigma sigma))))))

積分布関数の近似式(Abramowitz and Stegun)

確率密度関数は定義から簡単に計算できるのですが、累積分布関数は確率密度関数を-∞〜xの範囲で積分しないといけません。Wikipediaを見ると誤差関数erfを使った式が載っていますが、このerfの計算が面倒です。

 \displaystyle
P(x) = \frac{1}{2}\left(1 + \mathrm{erf}\,\frac{x - \mu}{\sqrt{2\sigma^2}}\right)

erfの計算にはテイラー展開した式を利用する方法があるのですが(実は以前にその方法を使った関数を作っていました)、実行速度が遅くなるような気がしたので、erfを使わない累積分布関数を調べてみたところ、Abramowitz and Stegun (1964) Handbook of Mathematical Functions (https://personal.math.ubc.ca/~cbm/aands/)に近似値の計算式が載っていました(第26.2.17節)。

(defun pnorm-sub (x)
  (let* ((p 0.2316419)
         (t1 (/ 1 (+ 1 (* p x))))
         (t2 (* t1 t1))
         (t3 (* t1 t2))
         (t4 (* t1 t3))
         (t5 (* t1 t4)))
    (- 1 (* (dnorm x)
            (+ (*  0.319381530 t1)
               (* -0.356563782 t2)
               (*  1.781477937 t3)
               (* -1.821255978 t4)
               (*  1.330274429 t5))))))

(defun pnorm (x &optional (mu 0) (sigma 1))
  (cond ((= x 0) 0.5)
        ((> x 0) (pnorm-sub (/ (- x mu) sigma)))
        ((< x 0) (- 1 (pnorm-sub (/ (- (- x) mu) sigma))))))

試しにSBCL 2.2.7で実行してみます。

CL-USER> (pnorm 1.96 0 1)
0.9750021780440332d0

同じことをMatlab 2022bで実行すると

>> normcdf(1.96, 0, 1)
ans =
   0.975002104851780

となり、小数点以下8桁目から差が出ています。あらかじめ精度がわかっていれば実用上は問題ない気がします。

積分布関数の逆関数の近似式(Abramowitz and Stegun)

こちらも同様にAbramowitz and Stegunの第26.2.23節を参考にします。p ∈ [0, 1]を与えるとz値が返ってきます。0と1の時はそれぞれ負の無限大と正の無限大となります。((ここで使っているdouble-float-negative-infinityとかはSBCL特有のものなので、実装依存になっています。))

(defun qnorm-sub (p)
  (assert (and (>= p 0) (<= p 1)) nil "The argument p has to be bounded to [0, 1].")
  (cond ((= p 0) double-float-negative-infinity)
        ((= p 1) double-float-positive-infinity)
        (t (let* ((t1 (sqrt (log (/ 1 (* p p)))))
                  (t2 (* t1 t1))
                  (t3 (* t1 t2)))
             (- t1 (/ (+    2.515517
                            (* 0.802853 t1)
                            (* 0.010328 t2))
                      (+    1
                            (* 1.432788 t1)
                            (* 0.189269 t2)
                            (* 0.001308 t3))))))))

(defun qnorm (p &optional (mu 0) (sigma 1))
  (if (<= p 0.5)
      (+ (* -1 (qnorm-sub p) sigma) mu)
      (+ (* (qnorm-sub (- 1 p)) sigma) mu)))

SBCL 2.2.7だと

CL-USER> (qnorm 0.975 0 1)
1.9603953

Matlab 2022bだと

>> norminv(0.975, 0, 1)
ans =
   1.959963984540054

となり、こちらはけっこう誤差が大きいです。Abramowitz and Stegunにも誤差は最大で4.5e-4とあったので、しょうがないですね。

積分布関数と逆関数を誤差関数から計算する(Numerical Recipes

Press, Teukolsky, Vetterling, and Flannery (2007) Numerical Recipes 3rd editionにはチェビシェフ関数展開を用いたプログラムが載っていたので、第6.14.1節(p.320)を参考にして確率密度関数と累積分布関数を書き直してみました*1。この部分はNumerical Recipesを参考にしたというよりは、Wikipediaにも載っている数式どおりです。

(defun pnorm (x &optional (mu 0) (sigma 1))
  (/ (erfc (* (/ -1 (sqrt 2)) (/ (- x mu) sigma))) 2))

(defun qnorm (p &optional (mu 0) (sigma 1))
  (- mu (* (sqrt 2) sigma (inverfc (* 2 p)))))

erfc(x) = 1 - erf(x)という関係のerfcと、inverfc(erfc(x)) = xという関係になるinverfcの両方が必要です。Numerical Recipesの第6.2.2節(p.264)には誤差関数に関連した関数群がC++で書かれていますので、そちらをCommon Lispに直していきます。

(defun erf (x)
  "Error function."
  (if (> x 0)
      (- 1 (erfccheb x))
      (- (erfccheb x) 1)))

(defun erfc (x)
  "Complementary error function."
  (if (> x 0)
      (erfccheb x)
      (- 2 (erfccheb (- x)))))

(defun inverfc (p)
  "Inverse of complementary error function. Returns x such that erfc(x) = p for argument p between 0 and 2."
  (cond ((>= p 2) double-float-negative-infinity)
        ((<= p 0) double-float-positive-infinity)
        (t
         (let* ((pp (if (< p 1) p (- 2 p)))
                (tt (sqrt (* -2 (log (/ pp 2)))))
                (x (* (- (sqrt 1/2))
                      (- (/ (+ 2.30753 (* tt 0.27061))
                            (+ 1 (* tt (+ 0.99229 (* tt 0.04481)))))
                         tt))))
           (loop for j from 0 to 1
                 do (setq err (- (erfc x) pp))
                    (setq x (+ x (/ err (- (* 1.12837916709551257 (exp (- (* x x))))
                                           (* x err))))))
           (if (< p 1) x (- x))))))

(defun inverf (p)
  "Inverse of the error function. Returns x such that erf(x) = p for argument p between -1 and 1."
  (inverfc (- 1 p)))

とくにinvertcにはマジックナンバーがいくつも入っていて、どこから出てきたのかがわかりません。気持ち悪い。

erfcの中でerfcchebを呼び出しています。ここがNumerical Recipesで紹介されているチェビシェフ法のキモです。すでに関数展開がされていて、係数のリストが与えられていますので、それを使って計算していくだけです。やっていることはAbramowitz and Stegunの方法と同様で、展開された項がとても多いというところが違いです。

(defun erfccheb (z)
  "Approximation of erfc function using Chebyshev method with precalculated coefficients."
  (assert (> z 0) nil "erfccheb requires nonnegative argument")
  (let* ((cof #(-1.3026537197817094d0 6.4196979235649026d-1 1.9476473204185836d-2 -9.561514786808631d-3 -9.46595344482036d-4 3.66839497852761d-4 4.2523324806907d-5 -2.0278578112534d-5 -1.624290004647d-6 1.303655835580d-6 1.5626441722d-8 -8.5238095915d-8 6.529054439d-9 5.059343495d-9 -9.91364156d-10 -2.27365122d-10 9.6467911d-11 2.394038d-12 -6.886027d-12 8.94487d-13 3.13092d-13 -1.12708d-13 3.81d-16 7.106d-15 -1.523d-15 -9.4d-17 1.21d-16 -2.8d-17))
         (d 0)
         (dd 0)
         (tt (/ 2 (+ 2 z)))
         (ty (- (* 4 tt) 2)))
    (loop for j from (1- (length cof)) downto 1
          do (setq tmp d)
             (setq d (+ (* ty d) (* -1 dd) (aref cof j)))
             (setq dd tmp))
    (* tt (exp (+ (* -1 z z) (* 0.5 (+ (aref cof 0) (* ty d))) (- dd))))))

Apple MacBook Pro (2021, M1 Max)での計算の例をまとめてみました。いずれもMatlabの結果とはやや違いますが、Abramowitz and StegunよりはNumerical RecipesのほうがMatlabの計算値に近い結果になっていますし、Rの計算結果とも一致しています。特に逆関数の方は3〜4桁分も改善しました。また、誤差関数が必要になる他のケースもあるので、高速かつ精度の高い誤差関数が手に入ったのはよかったです。

(pnorm 1.96 0 1)
0.9750021780440332 ;Abramowitz and Stegun, SBCL 2.2.7
0.9750021119513617 ;Numerical Recipes, SBCL 2.2.7
0.975002104851780  ;Matlab 2022b
0.9750021          ;R 4.2.2

(qnorm 0.975 0 1)
1.9603953         ;Abramowitz and Stegun, SBCL 2.2.7
1.959964358933026 ;Numerical Recipes, SBCL 2.2.7
1.959963984540054 ;Matlab 2022b
1.959964          ;R 4.2.2

*1:Numerical Recipesに掲載されているコードを使用する際にはライセンスに注意しないといけないようです。ただ、そのライセンスはC++のコードをそのまま使うときに問題になるのか、今回のようなCommon Lispへの移植でもライセンス条項が当てはまるのかがわかりませんでした。