もう忘れないスターリングの公式
以下の形のスターリングの公式は有名です:
この公式はある程度精度が高く数理統計学や統計物理における漸近公式などに使用される一方で, ちゃんと示そうとする方法はいずれもそこそこ骨が折れ,
この記事では, 前提知識をなるべく減らし, 上記の主要項をラフに導出する方法について述べます
前提知識
高校数学に加えて, 微積分の以下の事実を抑えておきます.
- 2次のテイラー展開:
\phi(x) \sim \phi(t) + \phi'(t)(x-t) + \dfrac{\phi''(t)}{2} (x-t)^2 - 階乗のガンマ関数による表示:
n! = \Gamma(n+1) = \int_0^\infty x^n e^{-x}\,dx - 平均
, 分散\mu\in \mathbf{R} に対する正規分布におけるガウス積分:\sigma^2 >0
\int_{-\infty}^\infty \exp{\dfrac{(x-\mu)^2}{2\sigma^2}}\, dx = \sqrt{2\pi \sigma^2} - (形式計算を受け入れるおおらかな心)
形式計算
まず, ガンマ関数を用いて積分により書き直します
ただしここで
により最大値は
この関数を指数関数に含む定積分
という感じになります. このことから
であり,
が従います. ここで, 定積分
を得ます. 以上の議論をまとめると,
が得られました.
積分が漸近する様子の可視化
上記の議論においては,
の定積分は 【
被積分関数を Python で可視化すると, 以下のような gif が得られます:
漸近の仕方は左右対称ではないものの, 形式計算における 【
Google Colaboratory で以下を用いると, 手元の環境でも試すことができます
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
# 初期化
fig, ax = plt.subplots()
line, = ax.plot([], [], lw=2)
vline = ax.axvline(x=0, color='gray', linestyle='--') # 破線の初期化
hline = ax.axhline(y=0, color='gray', linestyle='--') # 破線の初期化
title = ax.text(0.5, 1.05, "", transform=ax.transAxes, ha="center")
ax.set_xlabel("x")
ax.set_ylabel(r"$x^n e^{-x}$")
x_label = ax.text(0, 0, "", ha='center', va='bottom', fontsize=10, color='gray')
y_label = ax.text(0, 0, "", ha='left', va='top', fontsize=10, color='gray')
sigma_minus = ax.axvline(x=0, color='blue', linestyle=':', lw=1)
sigma_plus = ax.axvline(x=0, color='blue', linestyle=':', lw=1)
# 初期化関数
def init():
line.set_data([], [])
vline.set_xdata([0])
hline.set_ydata([0])
title.set_text("")
x_label.set_text("")
x_label.set_position((0, 0))
y_label.set_text("")
y_label.set_position((0, 0))
sigma_minus.set_xdata([0])
sigma_plus.set_xdata([0])
return line, vline, title
# 更新関数
def update(n):
x = np.linspace(0, 3 * n, 400)
y = x**n * np.exp(-x)
line.set_data(x, y)
sigma = np.sqrt(n)
sigma_minus.set_xdata([n - sigma])
sigma_plus.set_xdata([n + sigma])
xmin = max(0, n - 3 * sigma)
xmax = min(3 * n, n + 3 * sigma)
vline.set_xdata([n])
x_label.set_position((n, 0))
x_label.set_text(rf"$x = n$")
hline.set_ydata([(n/np.e)**n])
y_label.set_position((xmin, (n/np.e)**n))
y_label.set_text(rf"$y = (n/e)^{{n}}$")
ax.set_xlim(xmin, xmax)
ax.set_ylim(np.min(y), np.max(y) * 1.1)
title.set_text(rf"$f(x) = x^{{{n}}} e^{{-x}}$")
return line, vline, title
# アニメーション作成
ani = FuncAnimation(fig, update, frames=range(1, 26),
init_func=init, blit=True, interval=500)
# Colab からのダウンロード
from google.colab import files
ani.save("xn_exp_minus_x.gif", writer="pillow", fps=2)
files.download("xn_exp_minus_x.gif")
リンク
-
指数関数に
が乗っていることから,n における積分への寄与が急激に小さくなることが原因だが、形式計算なのであまり意識しなくて良い ↩︎x\nsim n
Discussion