pythonで３つのスターリングの近似式の精度を試してみた

得られた情報は一般に有限次元であるが、これを無限次元にして一般化して分析することがよくある。その際に、有限次元で一般化した無限次元との違いを知りたいときがよくある。その際に階乗の計算が必要になる。そこでスターリングの近似の精度が長い間気になっていたので試してみた。簡単な検証だが、驚くような結果になった。計算方法によりこんなに精度が違う。

最も簡単な近似は

$ \log n! = n \log n - n +O(\log n)$
である。

ガンマ関数$\Gamma (z)$について正の整数$n$に対して$\gamma(n)=(n-1)!$が成り立つ。ガンマ関数の計算機向け近似として次の式が成り立つ。

$2\ln \gamma(z) \sim \ln 2\pi -\ln z + z \left (  2\ln z+\ln  \left(z \sinh \frac{1}{z}+\frac{1}{810z^6} \right)-2\right )$
これは小数点以下８桁の精度をもつ。また、より単純な近似式として

$\ln \gamma(z) \sim 0.5\left(\ln 2\pi -\ln z\right) + z \left ( \ln \left(z+  \frac{1}{12z-\frac{1}{10z}}-1\right )\right)$
がある。

これらの精度をPythonで試してみる。

# 階乗

```python
%matplotlib inline
from scipy.special import factorial
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import statsmodels.api as sm
import itertools
for i in range(1,10,1):
    f=factorial(i)
    ff=1
    for ii in range(1,i+1,1):
        ff=ff*ii
    fff=np.exp(i*np.log(i)-i)
    j=i+1
    ffff=0.5*(np.log(2*np.pi)-np.log(j))+j*(np.log(j+1/(12*j-1/10/j))-1)
    fffff=0.5*(np.log(2*np.pi)-np.log(j)+j*(2*np.log(j)+np.log(j*np.sinh(1/j)+1/810/j**6)-2))
    print("{0:10d}{1:10.0f}{2: 10.0f}{3:15.8f}{4:18.8f}{5:18.8f}".format(i,f,ff,fff,np.exp(ffff),np.exp(fffff)))

         1         1         1     0.36787944        0.99998309        1.00000267
         2         2         2     0.54134113        1.99999509        2.00000037
         3         6         6     1.34425085        5.99999637        6.00000016
         4        24        24     4.68880356       23.99999516       24.00000014
         5       120       120    21.05608437      119.99999017      120.00000020
         6       720       720   115.64866155      719.99997253      720.00000040
         7      5040      5040   750.97400956     5039.99990097     5040.00000112
         8     40320     40320  5628.12896825    40319.99955910    40320.00000396
         9    362880    362880 47811.48664666   362879.99765201   362880.00001709

```
これらの結果の精度は驚くべきものだ.　さらに$n$の桁数を上げてみる。

```Python
for i in range(100,1000,100):
    f=factorial(i)
    ff=1
    for ii in range(1,i+1,1):
        ff=ff*ii
    fff=np.exp(i*np.log(i)-i)
    j=i+1
    ffff=0.5*(np.log(2*np.pi)-np.log(j))+j*(np.log(j+1/(12*j-1/10/j))-1)
    fffff=0.5*(np.log(2*np.pi)-np.log(j)+j*(2*np.log(j)+np.log(j*np.sinh(1/j)+1/810/j**6)-2))
    print(i,f,fff,np.exp(ffff),np.exp(fffff))

100 9.332621544394415e+157 3.720075976020918e+156 9.332621544394225e+157 9.332621544394755e+157
200 inf inf inf inf
300 inf inf inf inf
400 inf inf inf inf
500 inf inf inf inf
600 inf inf inf inf
700 inf inf inf inf
800 inf inf inf inf
900 inf inf inf inf
```
ffとfffについては適切に出力できないために割愛。ここまでで近似式の計算限界が示される。しかし、ln Factorialであればもっと先まで計算できる。

```Python
for i in range(100,1000,100):
    f=factorial(i)
    ff=1
    for ii in range(1,i+1,1):
        ff=ff*ii
    fff=i*np.log(i)-i
    j=i+1
    ffff=0.5*(np.log(2*np.pi)-np.log(j))+j*(np.log(j+1/(12*j-1/10/j))-1)
    fffff=0.5*(np.log(2*np.pi)-np.log(j)+j*(2*np.log(j)+np.log(j*np.sinh(1/j)+1/810/j**6)-2))
    print("{0:10d}{1: 18.8f}{2:18.8f}{2:18.8f}".format(i,fff,ffff,fffff))

       100      360.51701860      363.73937556      363.73937556
       200      859.66347331      863.23198719      863.23198719
       300     1411.13474240     1414.90584995     1414.90584995
       400     1996.58581884     2000.50069798     2000.50069798
       500     2607.30404921     2611.33045846     2611.33045846
       600     3238.15779313     3242.27533538     3242.27533538
       700     3885.75623453     3889.95083228     3889.95083228
       800     4547.68938213     4551.95073070     4551.95073070
       900     5222.15528699     5226.47551550     5226.47551550
```

```Python
for i in range(1000,10000,1000):
    f=factorial(i)
    ff=1
    for ii in range(1,i+1,1):
        ff=ff*ii
    fff=i*np.log(i)-i
    j=i+1
    ffff=0.5*(np.log(2*np.pi)-np.log(j))+j*(np.log(j+1/(12*j-1/10/j))-1)
    fffff=0.5*(np.log(2*np.pi)-np.log(j)+j*(2*np.log(j)+np.log(j*np.sinh(1/j)+1/810/j**6)-2))
    print("{0:10d}{1: 18.8f}{2:18.8f}{2:18.8f}".format(i,fff,ffff,fffff))

      1000     5907.75527898     5912.12817849     5912.12817849
      2000    13201.80491908    13206.52435051    13206.52435051
      3000    21019.10270295    21024.02485305    21024.02485305
      4000    29176.19856041    29181.26454459    29181.26454459
      5000    37585.96595708    37591.14350888    37591.14350888
      6000    46197.08848926    46202.35719906    46202.35719906
      7000    54975.65799626    54981.00377941    54981.00377941
      8000    63897.57456530    63902.98711266    63902.98711266
      9000    72944.81870687    72950.29014459    72950.29014459
```
階乗の対数をとった計算であればかなり大きな$n$まで計算ができる。

# 組み合わせ

同様の分析を組み合わせでやってみた。

```Python
def lnfactorial(i):
    j=i+1
    return 0.5*(np.log(2*np.pi)-np.log(j)+j*(2*np.log(j)+np.log(j*np.sinh(1/j)+1/810/j**6)-2))

n=10
for i in range(1,10,1):
    cnk=factorial(n)/factorial(i)/factorial(n-i)
    p=i/n
    cnk2=np.exp((-p*np.log(p)-(1-p)*np.log(1-p))*n)
    cnk3=np.exp(lnfactorial(n)-lnfactorial(i)-lnfactorial(n-i))
    print("{0:18.0f}{1:18.0f}{2:18.0f}".format(cnk,cnk2,cnk3))

                10                26                10
                45               149                45
               120               450               120
               210               837               210
               252              1024               252
               210               837               210
               120               450               120
                45               149                45
                10                26                10
```

```Python
n=100
for i in range(10,100,10):
    cnk=factorial(n)/factorial(i)/factorial(n-i)
    p=i/n
    cnk2=np.exp((-p*np.log(p)-(1-p)*np.log(1-p))*n)
    cnk3=np.exp(lnfactorial(n)-lnfactorial(i)-lnfactorial(n-i))
    print("{0:38.0f}{1:38.0f}{2:38.0f}".format(cnk,cnk2,cnk3))

                        17310309456440                       131272619177803                        17310309456022
                 535983370403809656832                5397605346934002810880                 535983370403677536256
            29372339821610946892136448           338453926949027754838851584            29372339821610843812921344
         13746234145802808216380243968        169248693106310399347688210432         13746234145802504751170977792
        100891344545564166887342866432       1267650600228231653296516890624        100891344545564518731063754752
         13746234145802806017356988416        169248693106310399347688210432         13746234145802700464240721920
            29372339821610942597169152           338453926949027754838851584            29372339821611264719716352
                 535983370403809656832                5397605346934002810880                 535983370403666067456
                        17310309456440                       131272619177803                        17310309456021
```

```Python
n=1000
for i in range(100,1000,100):
    cnk=factorial(n)/factorial(i)/factorial(n-i)
    p=i/n
    cnk2=np.exp((-p*np.log(p)-(1-p)*np.log(1-p))*n)
    cnk3=np.exp(lnfactorial(n)-lnfactorial(i)-lnfactorial(n-i))
    print("{0:*<10} {1:*<10} {2:*<10}".format(cnk,cnk2,cnk3))

nan******* 1.5196427575205235e+141 6.385051192635767e+139
nan******* 2.0989943697369585e+217 6.61715556065951e+215
nan******* 1.9724026353324775e+265 5.428250046407123e+263
nan******* 1.9286491874506976e+292 4.965272386250766e+290
nan******* 1.0715086071861939e+301 2.7028824094538324e+299
nan******* 1.9286491874506976e+292 4.965272386250766e+290
nan******* 1.9724026353324775e+265 5.428250046405889e+263
nan******* 2.0989943697369585e+217 6.617155560660263e+215
nan******* 1.5196427575205235e+141 6.385051192633227e+139
```

近似式の威力が実感できる。

参考：
[ウィキ：スターリングの近似](https://ja.wikipedia.org/wiki/%E3%82%B9%E3%82%BF%E3%83%BC%E3%83%AA%E3%83%B3%E3%82%B0%E3%81%AE%E8%BF%91%E4%BC%BC)