得られた情報は一般に有限次元であるが、これを無限次元にして一般化して分析することがよくある。その際に、有限次元で一般化した無限次元との違いを知りたいときがよくある。その際に階乗の計算が必要になる。そこでスターリングの近似の精度が長い間気になっていたので試してみた。簡単な検証だが、驚くような結果になった。計算方法によりこんなに精度が違う。
最も簡単な近似は
$ \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で試してみる。
階乗
%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$の桁数を上げてみる。
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であればもっと先まで計算できる。
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
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$まで計算ができる。
組み合わせ
同様の分析を組み合わせでやってみた。
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
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
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
近似式の威力が実感できる。
参考：
ウィキ：スターリングの近似