ガンマ関数 - 丸井綜研

非負整数nの階乗n!はn×(n-1)×(n-2)×(n-3)×...×2×1と定義されています。また、0!=1とします。例えば6!であれば、6×5×4×3×2×1=720です。これを、整数以外にも拡張したのがガンマ関数で、nが整数の時にはΓ(n+1) = n!となります。例えばΓ(6)=120、Γ(7)=720となり、その中間のΓ(6.5)は約287.885になります。整数のところでしか値がない階乗を、非整数のところも補間していったというイメージです。

このガンマ関数は統計に出てくる様々な関数を計算する際にも基本となるものです。例えば組み合わせの数（二項係数）には以下のような性質があります。大規模なABテストのようにnやkが大きくなる時には、愚直に階乗を計算するよりもガンマ関数の近似値で計算したほうが良い場合がありそうです。

$\displaystyle \binom{n}{k} = \frac{n!}{(n-k)!k!} = \frac{\Gamma\left(n+1\right)}{\Gamma\left(k+1\right)\Gamma\left(n-k+1\right)}$

というわけで、ガンマ関数をのちのち使えるように、Common Lispで実装してみます。正規分布の関数を実装したとき同様にAbramowitz and StegunやNumerical Recipesなどに頼ってもいいのですが、たまたまarXivで見つけた以下の論文があったので、そちらを参考にしてみます。

arxiv.org

これはプレプリントですが、査読を通過したものは Causley, M.F. The gamma function via interpolation. Numer Algor 90, 687–707 (2022). https://doi.org/10.1007/s11075-021-01204-8 のようです。ペイウォールに阻まれて読めなかったので、実装はarXivのものを参考にします。

この論文はAdaptive Antoulous Andersonアルゴリズムという方法を使って、倍精度で従来法と同じ精度を保ちつつも10～15パーセントの高速化を実現したというものっぽいです。紹介されていた従来法のMatlabコードのほうがAAA版よりも簡単かつ精度は同じということだったので、そちらを実装しました。（論文を参照した意味がないという……）

;;; Gamma Function
;;;
;;; Reference:
;;; Causley, M.F. The gamma function via interpolation.
;;; Numer Algor 90, 687–707 (2022).
;;; https://doi.org/10.1007/s11075-021-01204-8
;;; preprint: https://arxiv.org/abs/2104.00697
;;;
(defun gamma (z)
  "gamma function using a shifted Stiring series"
  (let* ((k 5.0d0)
         (n 16.0d0)
         (c (/ (log (* 2 pi)) 2))
         (a '(1/12 -1/360 1/1260 -1/1680 1/1188))
         (j (<= (realpart z) 1/2))      ;reflection
         (zz (if j (- 1 z) z))
         (tt (+ zz n))
         (gam (/ (car a) tt))
         (u 1.0d0)
         (lg 0.0d0)
         (g 0.0d0))
    (loop for i from 1 to (1- k)
          do (setq tt (* tt (expt (+ zz n) 2)))
             (setq gam (+ gam (/ (nth i a) tt))))
    (loop for i from 0 to (1- n)
          do (setq u (* u (+ zz i))))
    (setq lg (+ c (- (+ zz n)) (* (+ zz n -1/2) (log (+ zz n))) (- (log u)) gam))
    (setq g (exp lg))
    (when j                             ;reflection
      (let* ((m (truncate (realpart zz))))
        (setq g (/ (* (expt -1 m) pi) (* g (sin (* pi (- zz m))))))))
    g))

もともとMatlabのコードが与えられていたので、それを移植するのは難しくありませんでしたが、つまづいたのがnやkの定数の指定でした。最初は(let ((k 5) (n 16)...のように整数を指定していましたが、そのままだとMatlabで実行した結果と食い違いが大きいのです。割り算があっても有理数のまま扱われるので精度が高くなるのかと思ったのですが、そもそも倍精度で数値計算するためのアルゴリズムなので、ちゃんと倍精度の値だと指定しないといけなかったようです。あるいは有理数を浮動小数点数に変換するときに明示的に倍精度にしないといけなかったのかもしれません。

MatlabにもRにもgamma関数がありますが、どちらも引数は実数に限られており、複素数までは拡張されていないようです。上記のコードであれば複素数領域のガンマ関数も計算できるので、ちょっとお得な気分です。

計算結果を見てみましょう。まずは上記のコードから。

CL-USER> (mapcar #'gamma (linspace 1 20 20))
(1.0000000000000049d0 0.9999999999999994d0 2.0000000000000027d0
 5.999999999999997d0 24.000000000000153d0 119.99999999999997d0
 720.0000000000047d0 5040.00000000002d0 40319.99999999981d0
 362880.00000000047d0 3628799.999999996d0 3.991680000000017d7
 4.790016000000061d8 6.2270208000000105d9 8.717829120000078d10
 1.307674367999989d12 2.092278988800013d13 3.5568742809600825d14
 6.402373705727994d15 1.2164510040883208d17)

そしてMatlab。

>> for n=1:20
disp(gamma(n))
end
     1
     1
     2
     6
    24
   120
   720
        5040
       40320
      362880
     3628800
    39916800
   479001600
     6.227020800000000e+09
     8.717829120000000e+10
     1.307674368000000e+12
     2.092278988800000e+13
     3.556874280960000e+14
     6.402373705728000e+15
     1.216451004088320e+17

だいたい一致してますね。論文には13桁くらいの精度はあると書かれているので、僕の用途には十分すぎます。