スターリングの公式

外部リンクで済ますのは気持ち悪いので自分の言葉で書いてみる。

[前の記事へ]  [統計力学の目次へ]  [次の記事へ]


まずは普通のスターリングの公式

 統計力学では良く階乗の計算が出てくる。それが『アボガドロ数程度の大きさの整数』の階乗だったりするものだから、定義に従って具体的な値を算出するのは不可能であるし、もう桁が幾つか違っていてもどうでもいいやと思えるほどに大き過ぎる値となるだろう。

 これをそのまま使っていたのでは式変形をしていくことが無理なので、別の変形しやすい形で表しておいた方が良い。そこで利用されるのが、「スターリングの公式」と呼ばれる次のような近似公式である。

\[ \begin{align*} n!\ \kinji\ \left(\frac{n}{e}\right)^n \tag{1} \end{align*} \]

 あるいはこれと全く同じ意味だが、次のように書いてもいい。

\[ \begin{align*} \log_e (n!) \ \kinji\ n \big( \log_e n - 1 \big) \tag{2} \end{align*} \]

 これを求めるのは非常に簡単である。まず、\( n! \)の対数を取ったものは次のように表せるだろう。

\[ \begin{align*} \log_e(n!)\ &=\ \log_e ( 1 \cdot 2 \cdots n) \\ &=\ \log_e 1 \ +\ \log_e 2 \ +\ \cdots \ +\ \log_e n \end{align*} \]
 この右辺の値の意味するものを図で表してやろう。次のように、幅が 1 であるような帯グラフを描いてやれば、その帯の面積の合計が、その値を表していることになる。

 この図に描かれた曲線は\( y = \log_e x \)を表しており、帯グラフを描く参考のために引いてみたものである。しかし良く考えてみれば、この曲線の関数を積分したものは、帯の合計の面積に非常に近いものとなるではないか。\( \log_e x \)の積分をするためには、\( \log_e x = 1 \times \log_e x \)と考えて部分積分を使ってやればいいのだ。それで次のような計算ができる。

\[ \begin{align*} \log_e(n!)\ &\kinji\ \int_{\frac{1}{2}}^{n+\frac{1}{2}} \log_e x \diff x \\ &=\ \Big[x \log_e x \Big]_{\frac{1}{2}}^{n+\frac{1}{2}}\ -\ \int_{\frac{1}{2}}^{n+\frac{1}{2}} \diff x \\ &=\ \left(n+\frac{1}{2} \right) \log_e \left(n+\frac{1}{2} \right)\ -\ \frac{1}{2} \log_e \frac{1}{2}\ -\ \left(n+\frac{1}{2} \right)\ +\ \frac{1}{2} \end{align*} \]
 積分範囲についてはあまり細かいことを気にしてはいけない。どうせ\( n \)が大きくなれば、両端のわずかな差は気にならない程度となるのだから。

 上の式の定数部分や\( \frac{1}{2} \)などは\( n \)が大きくなれば無視できる程度になるので、これらを省いてやれば

\[ \begin{align*} \log_e(n!) \ \kinji\ n \log_e n \ -\ n \end{align*} \]
となり、(2) 式が出来上がるというわけだ。


精度の検証

 かなり大雑把に項を省いたりしたので、この近似式の精度がどの程度信頼して良いものなのか気にしないではいられないだろう。それで次のような表を作ってみた。

nn !(n/e)n誤差( 比 )
5 120 21.06 5.7倍
10 3.63 × 106 4.54 × 105 8.0倍
15 1.31 × 1012 1.34 × 1011 9.8倍
20 2.43 × 1018 2.16 × 1017 11.3倍
1009.33 × 101573.7 × 10156 25.2倍

 最初は誤差をパーセントで表示してやろうと思っていたが、実際にやってみたらパーセントで表すことが意味の無いほどに数値に開きがあることが分かった。n が大きくなるほどに差が減っていくのかと思ったが、予想に反してどんどん増えていくし、しかも何十倍もの差がある。それでもこのどんどん増える桁数に比べれば十分に近い値を保って追随しているのだと、・・・ちょっと言いにくいが、まぁ言えないことはない。

 \( n \)が増えるほど誤差が減って行くだろうと予想したのは先ほどのグラフを見て考えたせいだ。グラフは対数を取ったものだったので、対数で比較してみるのが筋というものだろう。それで次のような表を作ってみた。

nloge(n !) loge (n/e)n 誤差( % )
5 4.79 3.05 36.3
10 15.1 13.0 13.8
15 27.9 25.6 8.2
20 42.3 39.9 5.7
100363.7 360.5 0.9

 こうして見ると、\( n \)が大きくなるにつれて急速に真の値に近付いているのが分かる。この僅かな誤差が、対数を外した時には大きく出てしまうのだ。

 どちらにせよ、この近似はこれくらいのものなのだというのを感覚的に把握した上で使うことが大切であると思う。


精度の高いスターリングの公式

 上で説明した公式は普通に使うにはあまり問題はない精度があるのだが、場合によってはもっと高い精度で成り立つ式が必要となることもある。そういうときには次のような式を使う。

\[ \begin{align*} n!\ \kinji\ \sqrt{2\pi n} \left(\frac{n}{e}\right)^n \tag{3} \end{align*} \]

 あるいはこれと全く同じ意味だが、次のように書いてもいい。

\[ \begin{align*} \log_e (n!) \ \kinji\ n \big( \log_e n - 1 \big) \ +\ \frac{1}{2} \log_e (2\pi n ) \tag{4} \end{align*} \]

 \( \pi \)が一体どこから出てくるのか、と思うかも知れないが、これを求めるのは、そんなに難しい考察が要るわけではない。(1) (2) 式を求めるときに無視した部分を省かないで考慮に入れてやればいい。つまり、定数部分を\( c \)と置いて、

\[ \begin{align*} \log_e(n!) \ \kinji\ \left(n+\frac{1}{2} \right) \log_e \left(n+\frac{1}{2} \right)\ -\ n \ +\ c \end{align*} \]
としてやると、\( n! \)は、
\[ \begin{align*} n!\ &\kinji\ \left(n+\frac{1}{2} \right)^{n+\frac{1}{2}} e^{-n} e^c \\ &=\ \left(n+\frac{1}{2} \right)^n \left(n+\frac{1}{2} \right)^{\frac{1}{2}} e^{-n} e^c \\ &=\ n^n \left(1 + \frac{1}{2n} \right)^n \left(n+\frac{1}{2} \right)^{\frac{1}{2}} e^{-n} e^c \tag{5} \end{align*} \]
と書ける。ここで\( n \)が無限大に近付くときのことを考えてやると、
\[ \begin{align*} \lim_{n\to\infty} \left(n+\frac{1}{2} \right)^{\frac{1}{2}} \ =\ n^{\frac{1}{2}} \tag{6} \end{align*} \]
であるし、自然対数の底\( e \)の定義が
\[ \begin{align*} e \ \equiv\ \lim_{n\to\infty} \left( 1 + \frac{1}{n} \right)^n \end{align*} \]
であることから、
\[ \begin{align*} \lim_{n\to\infty} \left( 1 + \frac{1}{2n} \right)^{n} \ =\ \lim_{n\to\infty} \left( 1 + \frac{1}{n} \right)^{\frac{1}{2}n} \ =\ \sqrt{e} \tag{7} \end{align*} \]
が言えるので、これら (6) (7) 式を (5) 式に適用してやって、
\[ \begin{align*} n! \ \kinji\ a \sqrt{n} \left( \frac{n}{e} \right)^n \tag{8} \end{align*} \]
という形になる。(1) 式では 定数\( a \)\( \sqrt{n} \)という要素が抜けていたというわけだ。

 あとは、なぜ定数\( a \)\( \sqrt{2\pi} \)になるかというところを説明すれば終わりだが、ここで天下りで申し訳ないのだが、「ウォリスの公式」と呼ばれる次の式を使う。

\[ \begin{align*} \frac{\pi}{2} \ =\ \lim_{n\to\infty} \frac{ (2^n n! )^4 }{ (2n+1) ((2n!))^2 } \end{align*} \]
 この式に (8) 式を代入してやれば\( n \)はうまい具合に全部消えてしまって、結局\( a = \sqrt{2\pi} \)であることが分かるというわけだ。


精度の検証

 これについても精度を調べておこう。

nn ! (3)式右辺 誤差( 比 )
5 120 118 0.98倍
10 3.63 × 106 3.60 × 106 0.991倍
15 1.31 × 1012 1.30 × 1012 0.994倍
20 2.43 × 1018 2.42 × 1018 0.996倍
1009.33 × 101579.32 × 10157 0.999倍

 なんと!こちらは、初めからほぼ正確であるばかりか、\( n \)の増加に伴って、より正確さを増して行くのである。対数を取ったものは当然正確だろうが、一応まとめたものを書いておく。

nloge(n !) (4)式右辺 誤差( % )
5 4.79 4.77 0.35
10 15.1 15.1 0.05
15 27.9 27.9 0.02
20 42.3 42.3 0.01
100363.7 363.7 0.0003

 申し分の無い正確さである。これに手を加えたさらに正確な公式もまだあるのだが、ここでやる必要もないだろう。