I tried the accuracy of three Stirling's approximations in python

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.

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 ) It has a precision of 8 decimal places. Also, as a simpler approximation

\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) there is.

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

I tried the accuracy of three Stirling's approximations in python
I tried simulating the "birthday paradox" in Python
I tried the least squares method in Python
I tried to graph the packages installed in Python
I tried to summarize the string operations of Python
I tried "gamma correction" of the image with Python + OpenCV
I tried programming the chi-square test in Python and Java.
[Python] I tried to visualize the follow relationship of Twitter
I tried to implement the mail sending function in Python
I wrote the queue in Python
[Python] I tried collecting data using the API of wikipedia
I tried Line notification in Python
I tried to implement blackjack of card game in Python
I wrote the stack in Python
I tried to create a Python script to get the value of a cell in Microsoft Excel
I wrote a doctest in "I tried to simulate the probability of a bingo game with Python"
I tried scraping the ranking of Qiita Advent Calendar with Python
I compared the calculation time of the moving average written in Python
Movement that changes direction in the coordinate system I tried Python 3
I wrote the code to write the code of Brainf * ck in python
I tried to improve the efficiency of daily work with Python
I tried to implement PLSA in Python
I tried the asynchronous server of Django 3.0
Check the behavior of destructor in Python
I tried to implement permutation in Python
I tried to implement PLSA in Python 2
I tried using Bayesian Optimization in Python
I didn't know the basics of Python
The result of installing python in Anaconda
I tried to implement ADALINE in Python
The basics of running NoxPlayer in Python
I tried to implement PPO in Python
Python: I tried the traveling salesman problem
The Python project template I think of.
In search of the fastest FizzBuzz in Python
I tried the Python Tornado Testing Framework
I tried to summarize the contents of each package saved by Python pip in one line
I tried how to improve the accuracy of my own Neural Network
I want to batch convert the result of "string" .split () in Python
I want to explain the abstract class (ABCmeta) of Python in detail.
I tried to get the authentication code of Qiita API with Python.
(Python) I tried to analyze 1 million hands ~ I tried to estimate the number of AA ~
I tried to verify and analyze the acceleration of Python by Cython
I tried to streamline the standard role of new employees with Python
I tried to get the movie information of TMDb API with Python
I made a program to check the size of a file in Python
I tried to display the altitude value of DTM in a graph
I tried to implement a card game of playing cards in Python
I tried touching touch related methods in the scene module of pythonista
[Python / DynamoDB / boto3] List of operations I tried
I tried "smoothing" the image with Python + OpenCV
I tried hundreds of millions of SQLite with python
Output the number of CPU cores in Python
I tried the pivot table function of pandas
[Python] I tried substituting the function name for the function name
I tried cluster analysis of the weather map
[Python] Sort the list of pathlib.Path in natural sort
vprof --I tried using the profiler for Python
I tried "differentiating" the image with Python + OpenCV
Get the caller of a function in Python
Match the distribution of each group in Python