スターリングの近似を試してみた

スターリングの近似とは

スターリングの近似（スターリングの公式）とは$n!$の近似です。
$n!$の計算は、プログラミングでは再帰的な計算が必要になり$n$が大きくなるにつれ計算コストが増え大変になります。
そこでスターリングの近似式を利用することで大幅に計算コストを削減できます。

簡単なスターリング近似は(1)のようになります。

n! = \exp(n\log n - n + \varepsilon)\tag{1}

もしくは対数をとった(2)になります。

\log n! = n\log n - n + \varepsilon\tag{2}

$\varepsilon$は近似誤差であり$n$の値によって変わります（nによらない定数値ではありません）。
$\varepsilon$を0に近づけるため、もっと複雑にしたものもあります。
wikipedia等では(2)で説明されているため、ここでもスターリングの近似式は(2)として扱います。

スターリングの近似式の簡単な導出

※詳しい導出はwikipediaを参考にしてください
まず(2)を式変形し(3)のように総和を取ることができます。

\log n! = \sum_{k=1}^{n} \log k \tag{3}

その後(4)のように、$\log x$の$k-\frac{1}{2}$から$k+\frac{1}{2}$の積分が$\log k$と近似します。
（$\log x$の$k-\frac{1}{2}$から$k+\frac{1}{2}$の間が線形ではないため誤差が発生します。）

\sum_{k=1}^{n} \log k \approx \sum_{k=1}^{n} \int_{k-\frac{1}{2}}^{k+\frac{1}{2}} \log x dx \tag{4} 

\sum_{k=1}^{n} \int_{k-\frac{1}{2}}^{k+\frac{1}{2}} \log x dx = \int_{\frac{1}{2}}^{n+\frac{1}{2}} \log x dx 

= \left( n + \frac{1}{2} \right) \log \left( n+\frac{1}{2} \right) -\left( n+\frac{1}{2} \right)
 -\frac{1}{2} \log \frac{1}{2} +\frac{1}{2} \tag{5} 

最終的に(5)のようになり$n$を無限にした場合、定数項は無視され(6)になります

= n\log n - n　+ \varepsilon \tag{6}

pythonで検証

※あまり良いコードではありません

import matplotlib.pyplot as plt
import numpy as np
import scipy.misc as sm

plt.subplot(2, 1, 1)
x_size = 100
x_max = 10000
x = np.linspace(1, x_max, x_size)

# スターリング近似
appr_y = x * np.log(x) - x 
plt.plot(x, appr_y, label='Approximate')

# 真値（オーバーフローを避けるためややこしくなっています）
tv_y = np.zeros(x_size)
y_count = 0
for i in x:
    for j in range(int(i)):
        tv_y[y_count] += np.log(j+1)   
    y_count += 1
plt.plot(x, tv_y, label='True value')
plt.legend()


plt.subplot(2, 1, 2)
# 誤差
err_y = tv_y - appr_y
plt.plot(x, err_y, label='Error')

plt.legend()
plt.show()

図上がスターリング近似と真値になっておりほぼほぼ重なっており、うまく近似できていることが確認できます。真値から近似式を引いた近似誤差が図下になっています。
$\varepsilon$を0に近づけるため、もっと複雑にしたものがwikipediaにあるので必要に応じて変更してください。