At some point in our mathematical careers, we encounter the Stirling approximation: $n! \sim \sqrt{2\pi n} \left( \frac ne\right)^n\;\;\text{as }n\to\infty$ The $\left(\frac ne\right)^n$ makes some amount of sense, because heuristically, $n!$ should be similar to, but less than, $n^n$. It also corresponds to the more intuitive-looking logarithmic form of the Stirling approximation, $\ln n! \sim n\ln n - n$. But then there's the $\sqrt{2\pi n}$. What? Huh? Why is there a $\sqrt n$ there? What's $\pi$ doing there? ### Integral version The chief idea to approximating $\ln n!$ is to express it as an integral and try to approximate the integral. This is not very difficult; for integer $n$ we know that $n!$ can be expressed in terms of the Gamma function as $n! = \Gamma(n+1) = \int_0^\infty e^{-x} x^n dx$ For reasons that will be clear in a moment, we will write this as $n! = \int_0^\infty e^{-x + n\ln x} dx$ If we look at the graph of $n\ln x - x$, it has a hump before it crashes down into the negatives. This means that exponentiating it gives a reasonable value only in a small positive region, before the strongly negative exponent forces the integrand to be zero. This is the idea behind Laplace's approximation ### Laplace's approximation Suppose we want to approximate an integral of the form $\int_a^b e^{M f(x)} dx$ where $M$ is a large number. The basic heuristic is that this is going to behave basically like $e^{M \max{f(x)}}$ as $M$ grows without bound. Let's make this idea more concrete. Let $f(x)$ be a function with a global maximum at $x_0$. Then $f'(x_0) = 0$ and around that global maximum the function behaves like $f(x) = f(x_0) + \frac12 f''(x_0)(x_0 - x)^2$ (by the Taylor expansion). Since $x_0$ is a global maximum, the value of $f''(x_0)$ is negative, so we can write everything positively and say $f(x) = f(x_0) - \frac 12 f''(x_0) |x-x_0|^2$ Then our integral $\int_a^b e^{M f(x)} dx = e^{Mf(x_0)}\int_a^be^{-\frac12 Mf''(x_0)|x-x_0|^2} dx$which looks like a Gaussian integral! The main difference is that for us to apply the Gaussian formula, the limits of integration should be $\pm \infty$, but we can make those the limits with little error because the exponential decays very fast for large positive and negative $x$. Now we can use the Gaussian integral formula that $\int_{-\infty}^{\infty} e^{-ax^2}dx = \sqrt{\frac\pi a}$ and get $\int_{-\infty}^\infty e^{-\frac12 Mf''(x_0)|x-x_0|^2} dx = \sqrt{\frac{2\pi}{M|f''(x_0)|}}$so that the overall integral goes to $\int_a^b e^{M f(x)} dx = e^{Mf(x_0)}\int_a^be^{-\frac12 Mf''(x_0)|x-x_0|^2} dx=\sqrt{\frac{2\pi}{M|f''(x_0)|}}Me^{f(x_0)}$ ### Applying the approximation We can now apply this approximation to $e^{-x + n\ln x}$, where $n$ takes the role of the large number $n$. Thus $f(x) = -\frac xn + \ln x$, and its first derivative is $f'(x) = -\frac 1n + \frac 1x$, and the second derivative $f''(x) = -\frac1{x^2}$ which means that the maximum occurs at $x=n$ with the values $f(x_0) = -1 + \ln n$ and $f''(x_0) = -\frac 1{n^2}$. Plugging this into the formula gives $n! = \sqrt{\frac{2\pi}{n\left|\frac1{n^2}\right|}}e^{n(-1 + \ln n)} = \sqrt{2\pi n}\left(\frac ne\right)^n$ So _that's_ where the $\pi$ comes from! It comes from the Gaussian integral hidden in approximating the Gamma function. We could keep tracing this back and ask where $\pi$ comes from in the Gaussian integral, but [other people have done this better than me](https://www.youtube.com/watch?v=cy8r7WSuT1I). ### References See the [Wikipedia page on the Laplace approximation](https://en.wikipedia.org/wiki/Laplace%27s_method#Example:_Stirling's_approximation) for some helpful visualizations on this.