Editing Stirling's approximation (section)

== Alternative derivations ==
An alternative formula for <math>n!</math> using the [[gamma function]] is
<math display=block> n! = \int_0^\infty x^n e^{-x}\,{\rm d}x.</math>
(as can be seen by repeated integration by parts). Rewriting and changing variables {{math|''x'' {{=}} ''ny''}}, one obtains
<math display=block> n! = \int_0^\infty e^{n\ln x-x}\,{\rm d}x = e^{n \ln n} n \int_0^\infty e^{n(\ln y -y)}\,{\rm d}y.</math>
Applying [[Laplace's method]] one has
<math display=block>\int_0^\infty e^{n(\ln y -y)}\,{\rm d}y \sim \sqrt{\frac{2\pi}{n}} e^{-n},</math>
which recovers Stirling's formula:
<math display=block>n! \sim e^{n \ln n} n \sqrt{\frac{2\pi}{n}} e^{-n}
= \sqrt{2\pi n}\left(\frac{n}{e}\right)^n.
</math>

=== Higher orders ===
In fact, further corrections can also be obtained using Laplace's method. From previous result, we know that <math>\Gamma(x) \sim x^x e^{-x}</math>, so we "peel off" this dominant term, then perform two changes of variables, to obtain:<math display="block">x^{-x}e^x\Gamma(x) = \int_\R e^{x(1+t-e^t)}dt</math>To verify this: <math>\int_\R e^{x(1+t-e^t)}dt \overset{t \mapsto \ln t}{=} e^x \int_0^\infty t^{x-1} e^{-xt} dt \overset{t \mapsto t/x}{=} x^{-x} e^x \int_0^\infty e^{-t} t^{x-1} dt = x^{-x} e^x \Gamma(x)</math>.

Now the function <math>t \mapsto 1+t - e^t</math> is unimodal, with maximum value zero. Locally around zero, it looks like <math>-t^2/2</math>, which is why we are able to perform Laplace's method. In order to extend Laplace's method to higher orders, we perform another change of variables by <math>1+t-e^t = -\tau^2/2</math>. This equation cannot be solved in closed form, but it can be solved by serial expansion, which gives us <math>t = \tau - \tau^2/6 + \tau^3/36 + a_4 \tau^4 + O(\tau^5) </math>. Now plug back to the equation to obtain<math display="block">x^{-x}e^x\Gamma(x) = \int_\R e^{-x\tau^2/2}(1-\tau/3 + \tau^2/12 + 4a_4 \tau^3 + O(\tau^4)) d\tau = \sqrt{2\pi}(x^{-1/2} + x^{-3/2}/12) + O(x^{-5/2})</math>notice how we don't need to actually find <math>a_4</math>, since it is cancelled out by the integral. Higher orders can be achieved by computing more terms in <math>t = \tau + \cdots</math>, which can be obtained programmatically.{{NoteTag|note=For example, a program in Mathematica: <syntaxhighlight lang="mathematica">
series = tau - tau^2/6 + tau^3/36 + tau^4*a + tau^5*b;
(*pick the right a,b to make the series equal 0 at higher orders*)
Series[tau^2/2 + 1 + t - Exp[t] /. t -> series, {tau, 0, 8}]

(*now do the integral*)
integral = Integrate[Exp[-x*tau^2/2] * D[series /. a -> 0 /. b -> 0, tau], {tau, -Infinity, Infinity}];
Simplify[integral/Sqrt[2*Pi]*Sqrt[x]]

</syntaxhighlight>|name=mathematica-program|content=content|text=text}}

Thus we get Stirling's formula to two orders:<math display="block"> n! = \sqrt{2\pi n}\left(\frac{n}{e}\right)^n \left(1 + \frac{1}{12 n}+O\left(\frac{1}{n^2}\right) \right).
</math>

=== Complex-analytic version ===
A complex-analysis version of this method{{r|flajolet-sedgewick}} is to consider <math>\frac{1}{n!}</math> as a [[Taylor series|Taylor coefficient]] of the exponential function <math>e^z = \sum_{n=0}^\infty \frac{z^n}{n!}</math>, computed by [[Cauchy's integral formula]] as
<math display="block">\frac{1}{n!} = \frac{1}{2\pi i} \oint\limits_{|z|=r} \frac{e^z}{z^{n+1}} \, \mathrm dz. </math>

This line integral can then be approximated using the [[Method of steepest descent|saddle-point method]] with an appropriate choice of contour radius <math>r = r_n</math>. The dominant portion of the integral near the saddle point is then approximated by a real integral and Laplace's method, while the remaining portion of the integral can be bounded above to give an error term.

=== Using the Central Limit Theorem and the Poisson distribution ===
An alternative version uses the fact that the [[Poisson distribution]] converges to a [[normal distribution]] by the [[Central limit theorem|Central Limit Theorem]].<ref>{{Cite book |last=MacKay |first=David J. C. |title=Information theory, inference, and learning algorithms |date=2019 |publisher=Cambridge University Press |isbn=978-0-521-64298-9 |edition=22nd printing |location=Cambridge}}</ref>

Since the Poisson distribution with parameter <math>\lambda</math> converges to a normal distribution with mean <math>\lambda</math> and variance <math>\lambda</math>, their [[Probability density function|density functions]] will be approximately the same:

<math>\frac{\exp(-\mu)\mu^x}{x!}\approx \frac{1}{\sqrt{2\pi\mu}}\exp(-\frac{1}{2}(\frac{x-\mu}{\sqrt{\mu}}))</math>

Evaluating this expression at the mean, at which the approximation is particularly accurate, simplifies this expression to:

<math>\frac{\exp(-\mu)\mu^\mu}{\mu!}\approx \frac{1}{\sqrt{2\pi\mu}}</math>

Taking logs then results in:

<math>-\mu+\mu\ln\mu-\ln\mu!\approx -\frac{1}{2}\ln 2\pi\mu</math>

which can easily be rearranged to give:

<math>\ln\mu!\approx  \mu\ln\mu - \mu + \frac{1}{2}\ln 2\pi\mu</math>

Evaluating at <math>\mu=n</math> gives the usual, more precise form of Stirling's approximation.