I've been worried about the accuracy of Stirling's approximation for a long time, but I tried it. It's a simple verification, but the results are amazing. The accuracy is so different depending on the calculation method.

The simplest approximation is

$ \log n! = n \log n - n +O(\log n)$ It will be.

For the gamma function $ \ Gamma (z) $, $ \ gamma (n) = (n-1)! $ holds for the positive integer $ n $. The following equation holds as a computer approximation of the gamma function.

Try these precisions in 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
```

The accuracy of these results was amazing. Let's increase the number of digits of $ 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,np.exp(ffff),np.exp(fffff))
100 9.332621544394415e+157 9.332621544394225e+157 9.332621544394755e+157
200 inf inf inf
300 inf inf inf
400 inf inf inf
500 inf inf inf
600 inf inf inf
700 inf inf inf
800 inf inf inf
900 inf inf inf
```

ff and fff are omitted because they cannot be output properly. Up to this point, the calculation limits of the calculation approximation formula are shown. But with ln Factorial, you can calculate farther.

```
for i in range(1000,10000,1000):
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: 18.8f}{2:18.8f}".format(i,ffff,fffff))
1000 5912.12817849 5912.12817849
2000 13206.52435051 13206.52435051
3000 21024.02485305 21024.02485305
4000 29181.26454459 29181.26454459
5000 37591.14350888 37591.14350888
6000 46202.35719906 46202.35719906
7000 54981.00377941 54981.00377941
8000 63902.98711266 63902.98711266
9000 72950.29014459 72950.29014459
```

Calculations that take the logarithm of the factorial can be calculated up to a fairly large $ n $. This can greatly extend the computational limit, such as when using multiplicity with very small probabilities.

reference: Wiki: Stirling's approximation

Recommended Posts