【メモ帳】ガンマ関数とスターリングの公式

ガンマ関数

正の整数値に対してのみ定義されている階乗を、定義域を実数に広げて、微分可能な連続関数として拡張した物がガンマ関数 $\Gamma(n) = \int_0^\infty x^{n-1}e^{-x}dx$ ですが、階乗とは $n!=\Gamma(n+1)$ で結びついています。*1

ぱっと見、なんでこれが階乗になるのか意味不明な式ですが、 $n!$ の場合の被積分関数を少し見みてみることにして、これを $f(x)=x^ne^{-x}$ と置いておきます。

この被積分関数 $f(x)$ は、無限大に発散してゆく $x$ のべき $x^n$ と尻すぼみに減衰させる指数関数 $\exp(-x)$ の積になっています。べき関数よりも減衰関数の方が強いので、 $x$ を０から増やしてゆくと、被積分関数 $f(x)$ は初めは勢いよく増えますが、どこかで天井を迎えてその後０に向かって減衰してゆくことになります。値は常に正になっています。

その天井の x 座標は、 $f(x)$ の微係数が 0 になる所ですが、それは $f'(x) = e^{-x}x^{n-1}(n-x)$ なので $x = n$ と分かります。ここで、整数 $n$ がちゃんと形をとどめて出てくるので少し納得です。

fortran66.hatenablog.com

これが半無間区間で積分されるわけですが、どうせ被積分関数のピーク位置以外の寄与はカスみたいに小さいというのがパターンなので、今回もそんなところだろうと予想して、被積分関数の値 $f(n)=e^{-n}n^n$ を数値的に計算してガンマ関数と比較して見ることにします。

プログラム１

被積分関数のピーク値と積分した値に対応するガンマ関数の比を計算します。ガンマ関数は $(n-1)!$ と $n!$ 相当の二つの場合を採っておきます。 $e^{-n}n^n/\Gamma(n)$ と $e^{-n}n^n/\Gamma(n+1)$

    program gam1
      use, intrinsic :: iso_fortran_env 
      implicit none
      integer, parameter :: kq = real128
      real(kq) :: x
      integer  :: i
      do i = 1, 1000, 100
         x = i 
         print '(i3,2(x,g0))', i, exp(-x) * x**x / gamma(x), exp(-x) * x**x / gamma(x + 1)  
      end do    
    end program gam1

計算結果１

  1 .367879441171442321595523770161461 .367879441171442321595523770161461
101 4.00601365311257804891572563285707 .396635015159661192961953032956148E-001
201 5.65363854114232006542228612995131 .281275549310563187334442096017479E-001
301 6.91947392501882452722897920740193 .229882854651788190273388013534928E-001
401 7.98715292661699106820555066172287 .199180870987954889481435178596563E-001
501 8.92805158365303910433045178190687 .178204622428204373339929177283547E-001
601 9.77883421755885529035262693536473 .162709387979348673716349865813030E-001
701 10.5613016450838296424350939106298 .150660508489070323001927159923414E-001
801 11.2896672383911749737148370845054 .140944659655320536500809451741626E-001
901 11.9738080975208380196990875864777 .132894651470819511872353913279397E-001

計算結果を見ると、被積分関数のピーク値は $(n-1)!$ よりは大きく、積分後の値 $(n)!$ よりも（当たり前ですが）小さくなっています。そこで、戯れにその中間の $(n-1/2)!$ 相当の値と比較してみることにします。

プログラム２

    program gam2
      use, intrinsic :: iso_fortran_env 
      implicit none
      integer, parameter :: kq = real128
      real(kq) :: x
      integer  :: i
      do i = 1, 1000, 100
         x = i 
         print '(i3,x,g0)', i, exp(-x) * x**x / gamma(x + 0.5_kq)  
      end do    
    end program gam2

計算結果２

  1 .415107497420594703340268249441337
101 .399106893561339198388218895228938
201 .399024988331960555536648899830374
301 .398997508756508333925078382047327
401 .398983735395490147411070127245623
501 .398975460605820122182191678946810
601 .398969939583726224779254994017756
701 .398965993792399175599028757808330
801 .398963033242696722186803273419373
901 .398960729877773834988078478722955

すると、うまい具合に一定比　0.39896... に収束しているように見えます。0.39896 という値を電卓で適当に叩いて、逆数を取って二乗すると 6.2826... という値になったので多分 $2\pi$ であろうと思われます。

すなわち数値的には、ガンマ関数の階乗相当地点での被積分関数のピーク値は $f(n) / \Gamma(n+{1\over2}) \sim {1\over\sqrt{2\pi}}$ と求まりました。少しずれてはいますが、確かに被積分関数のピーク値が階乗と関係深いことが見て取れます。

スターリングの公式

スターリングの公式は、俗には階乗の近似式 $n!\sim n(\ln(n)-1)$ もしくはこれにもう少し係数などが付いた物を指します。

今、求めた式を少し移項した上で両辺の対数を取ってみると、たしかに $\ln\Gamma(n+{1\over2})=(n-{1\over2})!\sim n\ln(n)-n+\ln(\sqrt{2\pi})$ となってよく似た式になります。

あるいは $1/2$ だけずらすことにすると $n!=\Gamma(n+1)$ なので、 $n\rightarrow n+{1\over2}$ と置き換えて、 $n!=\sqrt{2\pi}f(n+1/2)\sim\sqrt{2\pi}e^{-(n+1/2)}(n+1/2)^{(n+1/2)}$ となります。

数値的に $\ln\Gamma(n+{1\over2})\sim n\ln(n)-n+\ln(\sqrt{2\pi})$ の差を取って比較してみます。

プログラム３

差 $\ln\Gamma(n+{1\over2})-\left({1\over2}\ln{2\pi}-n+n\ln(n)\right)$ とその差に $n$ を掛けたものを出力します。 $n$ を掛けるのはスターリングの公式の近似の高次項からの類比による戯れです。

    program gam3
      use, intrinsic :: iso_fortran_env 
      implicit none
      integer, parameter :: kq = real128
      real(kq), parameter :: pi = 4 * atan(1.0_kq)
      real(kq) :: x
      integer  :: i
      do i = 1, 1000, 100
         x = i 
         print '(i3,2(x,g0))', i, log(gamma(x + 0.5_kq)) - ( log(2 * pi) / 2 - x + x * log(x) ), & 
                                ( log(gamma(x + 0.5_kq)) - ( log(2 * pi) / 2 - x + x * log(x) ) ) * i 
      end do    
    end program gam3

計算結果

  1 -.397207708399179641258481821872649E-001 -.397207708399179641258481821872649E-001
101 -.412538895125284184356880268682438E-003 -.416664284076537026200449071369262E-001
201 -.207296549782888548831458473423111E-003 -.416666065063605983151231531580453E-001
301 -.138427374882831965831612061863506E-003 -.416666398397324217153152306209154E-001
401 -.103906861724177446060175670243060E-003 -.416666515513951558701304437674669E-001
501 -.831669800064454703178056020480413E-004 -.416666569832291806292206066260687E-001
601 -.693288850874847850447310043197465E-004 -.416666599375783558118833335961676E-001
701 -.594388897581945996700061322962287E-004 -.416666617204944143686742987396564E-001
801 -.520183057158577580123800189897250E-004 -.416666628784020641679163952107697E-001
901 -.462449097365583306083968409790626E-004 -.416666636726390558781655537221354E-001

計算結果をみると、差は結構小さくよく近似されていることが分かります。また高次の補正項が $-0.04166666... / n = -1/24n$ であることも分かります。

そこで調子に乗って、もう一段同じ手続きを進めます。

プログラム４

差 $\ln\Gamma(n+{1\over2})-\left({1\over2}\ln{2\pi}-n+n\ln(n)-{1\over24n}\right)$ とその差に $n^3$ を掛けたものを出力します。 $n^3$ は、計算結果からの逆算です。

    program gam4
      use, intrinsic :: iso_fortran_env 
      implicit none
      integer, parameter :: kq = real128
      real(kq), parameter :: pi = 4 * atan(1.0_kq)
      real(kq) :: x
      integer  :: i
      do i = 1, 1000, 100
         x = i 
         print '(i3,2(x,g0))', i, log(gamma(x + 0.5_kq)) - ( log(2 * pi) / 2 - x + x * log(x) - 1.0_kq / (24 * x) ), & 
                                ( log(gamma(x + 0.5_kq)) - ( log(2 * pi) / 2 - x + x * log(x) - 1.0_kq / (24 * x) ) ) * i**3 
      end do    
    end program gam4

計算結果４

  1 .194589582674870254081848447940181E-002 .194589582674870254081848447940181E-002
101 .235900012835689724514385849088463E-008 .243048019124623958856896254701693E-002
201 .299305005315181808524918841287283E-009 .243053652546747070948926446747635E-002
301 .891260273918649549369722899779121E-010 .243054706952683739145563255973093E-002
401 .376939439171983446954764391899663E-010 .243055077420759382117630006617249E-002
501 .193282185350048823554048332146143E-010 .243055249243288499852057323939240E-002
601 .111964863741344148635793726845698E-010 .243055342696905859466697952881662E-002
701 .705588053109556685876677566004617E-011 .243055399095148574762536148061177E-002
801 .472941897939918885290992188350704E-011 .243055435722780067273990289109438E-002
901 .332300511741232087996684894263910E-011 .243055460846571599191284911257351E-002

計算結果をみると、差は $10^{-11}$ 程度とかなり小さくなっていることが分かります。また高次の補正項は電卓をいじっていると $-0.00230554... / n^3 = -7/2880n^3$ であることが分かります。

すなわち $\ln{(n-{1\over2})!}=\ln\Gamma(n+{1\over2})\sim\left({1\over2}\ln{2\pi}-n+n\ln(n)-{1\over24n}+{7\over2880n^3}\right)$ となっていると考えられます。

ここで近似は 10^-17 程度と倍精度を超えるところまで来ており、これ以上は同じやり方を繰り返しても数値誤差のためにうまくゆきませんでした。これを普通のスターリングの公式の補正項と比較してみると、べきも二次の（逆数）項がないなど、こちらの方が補正項が小さいことが分かります。

結論

以上、ガンマ関数が被積分関数のピーク値の寄与しかないと考えることで、おおむね直観的に理解でき、しかもそこから出てくるスターリングの公式に似た式は、~~元のスターリングの公式よりも補正項が小さくなっており、より良い近似になっていることも分かりました~~。

ここでは数値的に補正項を類推しましたが、被積分関数をピーク値での１の幅の短冊積分になっていると思ってやれば、オイラー＝マクローリンの公式で導出できると思われます。

補足[H31.1.26]　休みの徒然に逆さまに迎えに行くと

スターリングの公式より、ベルヌーイ数を使った展開から $\log\Gamma(z)\simeq{1\over2}\log2\pi-z+(z-{1\over2})\log z +{1\over12z}-{1\over360z^3}+\cdots$ なので、 $z$ を $z+{1\over2}$ にまとめ直すとここで求めた係数が出る。

$(z-{1\over2})\log(z-{1\over2})=(z-{1\over2})\log z(1-{1\over2z})=(z-{1\over2})(\log z + \log(1-{1\over2z}))$

対数関数をテイラー展開してやれば
$(z-{1\over2})\log(1-{1\over2z})=(z-{1\over2})(\log z - ( {1\over2z}+{1\over2(2z)^2}+{1\over3(2z)^3}+\cdots))$

開いて移項すれば
$(z-{1\over2})\log z=(z-{1\over2})\log(z-{1\over2})+{1\over2}+({1\over8z}-{1\over4z})+({1\over24z^2}-{1\over16z^2})+({1\over64z^3}-{1\over48z^3})+\cdots$

これを元の式に代入すると、
$\log\Gamma(z)\simeq{1\over2}\log2\pi-(z-{1\over2})+(z-{1\over2})\log(z-{1\over2})+({1\over8}-{1\over4}+{1\over12}){1\over z}+({1\over24}-{1\over16}){1\over z^2}+({1\over64}-{1\over48}-{1\over360}){1\over z^3}+\cdots$

係数を整理して、
$\log\Gamma(z)\simeq{1\over2}\log2\pi-(z-{1\over2})+(z-{1\over2})\log(z-{1\over2})-{1\over24z}-{1\over48z^2}-{23\over2880z^3}+\cdots$

ここで ${1\over24z}$ の項の $z$ を $z-{1\over2}$ の形にすると、
${1\over24(z-{1\over2})}={1\over24z}(1+{1\over2z}+{1\over4z^2})+\cdots={1\over24z}+{1\over48z^2}+{1\over96z^3}+\cdots$
なので、これによって置き換えると二乗の項は消え、また ${1\over96}={30\over2880}$ なので、
$\log\Gamma(z)\simeq{1\over2}\log2\pi-(z-{1\over2})+(z-{1\over2})\log(z-{1\over2})-{1\over24(z-{1\over2})}+{7\over2880z^3}+\cdots$
となる。 $z$ の三次以上の項を無視するなら ${1\over z^3}$ を ${1\over (z-{1\over2})^3}$ に置き換えてもいい。