史特靈公式

史特靈公式（英語：Stirling's formula）是一條用來取n階乘近似值的數學公式。一般來說，當n很大的時候，n階乘的計算量十分大，所以史特靈公式十分好用，而且，即使在n很小的時候，史特靈公式的取值已經十分準確。這個公式以詹姆斯·史特靈（英語：James Stirling (mathematician)）的名字命名，雖然亞伯拉罕·棣美弗早於史特靈提出了一個類似的公式，但結果較不精確。^[1][2][3]

當n增加時，(ln n!)與o (n ln n − n)之比趨於1

階乘與史特靈公式

${\displaystyle n!}$（藍色）、${\displaystyle {\sqrt {2\pi n}}\left({\frac {n}{e}}\right)^{n}}$（橘色），數字越大${\displaystyle {\sqrt {2\pi n}}\left({\frac {n}{e}}\right)^{n},}$會越趨近${\displaystyle n!}$。但${\displaystyle {\sqrt {2\pi n}}\left({\frac {n}{e}}\right)^{n}}$在負值則會因為出現虛數而無法使用

史特靈公式為：

${\displaystyle n!\approx {\sqrt {2\pi n}}\,\left({\frac {n}{e}}\right)^{n}.}$

這就是說，對於足夠大的整數n，這兩個數互為近似值。更加精確地：

${\displaystyle \lim _{n\rightarrow \infty }{\frac {n!}{{\sqrt {2\pi n}}\,\left({\frac {n}{e}}\right)^{n}}}=1}$

或

${\displaystyle \lim _{n\rightarrow \infty }{\frac {e^{n}\,n!}{n^{n}{\sqrt {n}}}}={\sqrt {2\pi }}.}$

歷史

這個公式是亞伯拉罕·棣莫弗首先發現的，形式為：

${\displaystyle \displaystyle n!\approx cn^{n+{\frac {1}{2}}}e^{-n}}$，其中c為常數。

斯特林（英語：James_Stirling_(mathematician)）證明了公式中${\displaystyle c={\sqrt {2\pi }}}$，約為2.506628274631。更加精確的形式是雅克·比內發現的。

推導

方法一

這個公式，以及誤差的估計，可以推導如下。首先不直接估計n!，而是考慮它的自然對數：

${\displaystyle \ln(n!)=\ln 1+\ln 2+\cdots +\ln n.}$

即:

${\displaystyle \ln(n!)-{\frac {\ln n}{2}}=\ln 1+\ln 2+\cdots +\ln n-{\frac {\ln n}{2}}.}$

這個方程的右面是積分${\displaystyle \int _{1}^{n}\ln(x)\,dx=n\ln n-n+1}$的近似值（利用梯形法則），而它的誤差由歐拉-麥克勞林公式給出：

${\displaystyle \ln(n!)-{\frac {\ln n}{2}}=\ln 1+\ln 2+\cdots +\ln(n-1)+{\frac {\ln n}{2}}=n\ln n-n+1+\sum _{k=2}^{m}{\frac {B_{k}{(-1)}^{k}}{k(k-1)}}\left({\frac {1}{n^{k-1}}}-1\right)+R_{m,n}.}$

其中B_k是伯努利數，R_m,n是歐拉-麥克勞林公式中的餘項。取極限，可得：

${\displaystyle \lim _{n\to \infty }\left(\ln n!-n\ln n+n-{\frac {\ln n}{2}}\right)=1-\sum _{k=2}^{m}{\frac {B_{k}{(-1)}^{k}}{k(k-1)}}+\lim _{n\to \infty }R_{m,n}.}$

把這個極限記為y。由於歐拉-麥克勞林公式中的餘項R_m,n滿足：

${\displaystyle R_{m,n}=\lim _{n\to \infty }R_{m,n}+O\left({\frac {1}{n^{2m-1}}}\right),}$

其中用到了大O符號，與以上的方程結合，便得出對數形式的近似公式：

${\displaystyle \ln n!=n\ln \left({\frac {n}{e}}\right)+{\frac {\ln n}{2}}+y+\sum _{k=2}^{m}{\frac {B_{k}{(-1)}^{k}}{k(k-1)n^{k-1}}}+O\left({\frac {1}{n^{2m-1}}}\right).}$

兩邊取指數，並選擇任何正整數m，便得到了一個含有未知數e^y的公式。當m=1時，公式為：

${\displaystyle n!=e^{y}{\sqrt {n}}~{\left({\frac {n}{e}}\right)}^{n}\left[1+O\left({\frac {1}{n}}\right)\right]}$

將上述表達式代入沃利斯乘積公式，並令n趨於無窮，便可以得出e^y（${\displaystyle e^{y}={\sqrt {2\pi }}}$）。因此，我們便得出斯特林公式：

${\displaystyle n!={\sqrt {2\pi n}}~{\left({\frac {n}{e}}\right)}^{n}\left[1+O\left({\frac {1}{n}}\right)\right]}$

這個公式也可以反覆使用分部積分法來得出，首項可以通過最速下降法得到。把以下的和

${\displaystyle \ln(n!)=\sum _{j=1}^{n}\ln j}$

用積分近似代替，可以得出不含${\displaystyle {\sqrt {2\pi n}}}$的因子的斯特林公式（這個因子通常在實際應用中無關）：

${\displaystyle \sum _{j=1}^{n}\ln j\approx \int _{1}^{n}\ln x\,dx=n\ln n-n+1.}$

方法二

考慮如下近似關係：

${\displaystyle {\begin{aligned}\ln n!&=\sum _{k=1}^{n}\ln k\\&\approx \int _{1}^{n}\ln x\mathrm {d} x=n\ln n-n+1\\&\approx n\ln n-n\end{aligned}}}$

以及階乘的積分形式（Γ函數）：

${\displaystyle n!=\int _{0}^{\infty }\mathrm {e} ^{-x}x^{n}\mathrm {d} x}$

其中被積函數取對數後求導結果可寫成：

${\displaystyle {\frac {\mathrm {d} }{\mathrm {d} x}}\ln(\mathrm {e} ^{-x}x^{n})={\frac {\mathrm {d} }{\mathrm {d} x}}(n\ln x-x)={\frac {n}{x}}-1}$

因此被積函數只有在${\displaystyle x=n}$附近才趨於平穩，現令${\displaystyle x\equiv n+\xi }$且${\displaystyle \xi \ll n}$，並改寫成如下形式：

${\displaystyle \ln(x^{n}\mathrm {e} ^{-x})=n\ln x-x=n\ln(n+\xi )-(n+\xi )}$

於是有：

${\displaystyle {\begin{aligned}\ln(n+\xi )&=\ln \left[n\left(1+{\frac {\xi }{n}}\right)\right]=\ln n+\ln \left(1+{\frac {\xi }{n}}\right)\\&=\ln n+{\frac {\xi }{n}}-{\frac {\xi ^{2}}{2n^{2}}}+\cdots \end{aligned}}}$

因此可以得到：

${\displaystyle {\begin{aligned}\ln(x^{n}\mathrm {e} ^{-x})&=n\ln(n+\xi )-(n+\xi )\\&=n\ln n+\xi -{\frac {\xi ^{2}}{2n}}-n-\xi +\cdots \\&=n\ln n-n-{\frac {\xi ^{2}}{2n}}+\cdots \end{aligned}}}$

兩邊同時作指數函數運算，變成：

${\displaystyle {\begin{aligned}x^{n}\mathrm {e} ^{-x}&\approx \mathrm {e} ^{n\ln n}\mathrm {e} ^{-n}\mathrm {e} ^{-{\frac {\xi ^{2}}{2n}}}\\&=n^{n}\mathrm {e} ^{-n}\mathrm {e} ^{-{\frac {\xi ^{2}}{2n}}}\end{aligned}}}$

代入到階乘的積分表示形式，得到：

${\displaystyle {\begin{aligned}n!&\approx \int _{-n}^{\infty }n^{n}\mathrm {e} ^{-n}\mathrm {e} ^{-{\frac {\xi ^{2}}{2n}}}\mathrm {d} \xi \\&\approx n^{n}\mathrm {e} ^{-n}\int _{-\infty }^{\infty }\mathrm {e} ^{-{\frac {\xi ^{2}}{2n}}}\mathrm {d} \xi \end{aligned}}}$

約等號後側得到的是高斯積分的形式，計算該積分得出斯特林公式：

${\displaystyle n!\approx n^{n}\mathrm {e} ^{-n}{\sqrt {2\pi n}}}$

收斂速率和誤差估計

y軸表示截斷的斯特林級數的相對誤差，x軸表示所使用的項數。

更加精確的近似公式為：

${\displaystyle n!={\sqrt {2\pi n}}\;\left({\frac {n}{e}}\right)^{n}e^{\lambda _{n}}}$

其中：

${\displaystyle {\frac {1}{12n+1}}<\lambda _{n}<{\frac {1}{12n}}.}$

斯特林公式實際上是以下級數（現在稱為斯特林級數）的第一個近似值：

${\displaystyle n!={\sqrt {2\pi n}}\left({n \over e}\right)^{n}\left(1+{1 \over 12n}+{1 \over 288n^{2}}-{139 \over 51840n^{3}}-{571 \over 2488320n^{4}}+\cdots \right).}$

當${\displaystyle n\to \infty }$時，截斷級數的誤差等於第一個省略掉的項。這是漸近展開式的一個例子。它不是一個收斂級數；對於任何特殊值n，級數的準確性只在取有限個項時達到最大，如果再取更多的項，則準確性將變得越來越差。

階乘的對數的漸近展開式也稱為斯特林級數：

${\displaystyle \ln n!=n\ln n-n+{1 \over 2}\ln(2\pi n)+{1 \over 12n}-{1 \over 360n^{3}}+{1 \over 1260n^{5}}-{1 \over 1680n^{7}}+\cdots .}$

在這種情況下，級數的誤差總是與第一個省略掉的項異號，且最多同大小。

伽瑪函數的斯特林公式

對於所有正整數，有：

${\displaystyle n!=\Pi (n)=\Gamma (n+1).}$

然而，伽瑪函數與階乘不一樣，它對於所有複數都有定義。儘管如此，斯特林公式仍然適用。如果${\displaystyle \Re (z)>0}$，那麼：

${\displaystyle \ln \Gamma (z)=\left(z-{\frac {1}{2}}\right)\ln z-z+{\frac {\ln {2\pi }}{2}}+2\int _{0}^{\infty }{\frac {\arctan {\frac {t}{z}}}{\exp(2\pi t)-1}}\,dt.}$

反覆使用分部積分法，可得以下漸近展開式：

${\displaystyle \ln \Gamma (z)=\left(z-{\frac {1}{2}}\right)\ln z-z+{\frac {\ln {2\pi }}{2}}+\sum _{n=1}^{\infty }{\frac {(-1)^{n-1}B_{n}}{2n(2n-1)z^{2n-1}}}}$

其中B_n是第n個伯努利數。當${\displaystyle |\arg z|<\pi -\epsilon }$，其中ε是正數時，這個公式對於絕對值足夠大的z是適用的，當使用了最初m個項時，誤差項為${\displaystyle O\left(z^{-m-{\frac {1}{2}}}\right)}$。對應的近似值可以寫為：

${\displaystyle \Gamma (z)={\sqrt {\frac {2\pi }{z}}}~{\left({\frac {z}{e}}\right)}^{z}\left[1+O\left({\frac {1}{z}}\right)\right].}$

斯特林公式的收斂形式

欲得出斯特林公式的一個收斂形式，我們必須計算：

${\displaystyle \int _{0}^{\infty }{\frac {2\arctan {\frac {t}{z}}}{\exp(2\pi t)-1}}\,dt=\ln \Gamma (z)-\left(z-{\frac {1}{2}}\right)\ln z+z-{\frac {1}{2}}\ln(2\pi ).}$

一種方法是利用含有上升階乘冪的級數。如果${\displaystyle z^{\overline {n}}=z(z+1)\cdots (z+n-1)}$，那麼：

${\displaystyle \int _{0}^{\infty }{\frac {2\arctan {\frac {t}{z}}}{\exp(2\pi t)-1}}\,dt=\sum _{n=1}^{\infty }{\frac {c_{n}}{(z+1)^{\overline {n}}}}}$

其中：

${\displaystyle c_{n}={\frac {1}{n}}\int _{0}^{1}x^{\overline {n}}\left(x-{\frac {1}{2}}\right)\,dx.}$

從中可以得出斯特林級數的一個收斂形式：

${\displaystyle \ln \Gamma (z)=\left(z-{\frac {1}{2}}\right)\ln z-z+{\frac {\ln {2\pi }}{2}}}$

${\displaystyle {}+{\frac {1}{12(z+1)}}+{\frac {1}{12(z+1)(z+2)}}+{\frac {59}{360(z+1)(z+2)(z+3)}}+{\frac {29}{60(z+1)(z+2)(z+3)(z+4)}}+\cdots }$

它在${\displaystyle \Re (z)>0}$時收斂。

適用於計算器的形式

以下的近似值

${\displaystyle \Gamma (z)\approx {\sqrt {\frac {2\pi }{z}}}\left({\frac {z}{e}}{\sqrt {z\sinh {\frac {1}{z}}+{\frac {1}{810z^{6}}}}}\right)^{z},}$

或

${\displaystyle 2\ln \Gamma (z)\approx \ln(2\pi )-\ln z+z\left[2\ln z+\ln \left(z\sinh {\frac {1}{z}}+{\frac {1}{810z^{6}}}\right)-2\right],}$

可以通過把斯特林公式整理，並注意到它的冪級數與雙曲正弦函數的泰勒級數展開式的相似性來得出。當z的實數部分大於8時，這個近似值精確到小數點後8位。2002年，Robert H. Windschitl建議計算器用這個公式來計算伽瑪函數。

Gergő Nemes在2007年提出了一個近似公式，它的精確度與Windschitl的公式相等，但更加簡單：

${\displaystyle \Gamma (z)\approx {\sqrt {\frac {2\pi }{z}}}\left[{\frac {1}{e}}\left(z+{\frac {1}{12z-{\frac {1}{10z}}}}\right)\right]^{z},}$

或

${\displaystyle \ln \Gamma (z)\approx {\frac {1}{2}}\left[\ln(2\pi )-\ln z\right]+z\left[\ln \left(z+{\frac {1}{12z-{\frac {1}{10z}}}}\right)-1\right]}$

參考文獻

• Abramowitz, M. and Stegun, I., Handbook of Mathematical Functions, http://www.math.hkbu.edu.hk/support/aands/toc.htm（頁面存檔備份，存於網際網路檔案館）
• 夏慧異，章家順, 史特靈公式在高等數學中的應用, http://ordostsg.org.cn:1080/KCMS/detail/detail.aspx?filename=XUSJ200904014&dbcode=CJFR&dbname=
• Paris, R. B., and Kaminsky, D., Asymptotics and the Mellin-Barnes Integrals, Cambridge University Press, 2001
• Whittaker, E. T., and Watson, G. N., A Course in Modern Analysis, fourth edition, Cambridge University Press, 1963. ISBN 0-521-58807-3
• Viktor T. Toth. "Programmable Calculators: Calculators and the Gamma Function". 2006 [2018-11-18]. （原始內容存檔於2007-02-23）., modified 2006
• 埃里克·韋斯坦因. Stirling's Approximation. MathWorld.
• Stirling's approximation at PlanetMath.

1. Dutka, Jacques, The early history of the factorial function, Archive for History of Exact Sciences, 1991, 43 (3): 225–249, doi:10.1007/BF00389433
2. Le Cam, L., The central limit theorem around 1935, Statistical Science, 1986, 1 (1): 78–96, JSTOR 2245503, MR 0833276, doi:10.1214/ss/1177013818; see p. 81, "The result, obtained using a formula originally proved by de Moivre but now called Stirling's formula, occurs in his 'Doctrine of Chances' of 1733."
3. Pearson, Karl, Historical note on the origin of the normal curve of errors, Biometrika, 1924, 16 (3/4): 402–404 [p. 403], JSTOR 2331714, doi:10.2307/2331714, I consider that the fact that Stirling showed that De Moivre's arithmetical constant was ${\displaystyle {\sqrt {2\pi }}}$ does not entitle him to claim the theorem, [...]