The n-Category Café

Nice! Now that you say it, I suppose it’s inevitable that the $\sqrt{\pi}$ in Stirling’s formula comes from the integral of $e^{-x^2}$.

Typos: just before the heading “Laplace’s method”, there are a couple of equals signs that should be $\sim$.

Regarding probabilistic interpretations, Dan Piponi’s comment here gives me a small scrap of an idea. Stirling’s formula can be rearranged as

$$\frac{n!}{n^n} \sim \sqrt{2\pi n}\,e^{-n},$$

the left-hand side being the probability that a random endomorphism of an $n$-element set is an automorphism. But why should the right-hand side be what it is? I have no idea. We need a Gaussian!

This post contains a quick probabilistic proof of Stirling’s formula using the central limit theorem. From the perspective of this proof, the deeper explanation is very roughly that:

1. $n!$ is a value of the $\Gamma$ function,

2. which is used to describe the $\Gamma$ distribution,

3. which is the distribution of a sum of independent exponential random variables,

4. which (by the CLT) is approximately Gaussian in distribution.

Another simple way to connect $n!$ to a Gaussian distribution is that the Euclidean norm squared of a two-dimensional Gaussian random vector has an exponential distribution, which has factorial moments. (The former of these two observations is closely related to the computation John mentioned of the Gaussian integral by squaring it and using polar coordinates.)

Therefore, if $X$ and $Y$ are independent Gaussian random variables with mean $0$ and variance $1/2$, then $$n! = \E (X^2 + Y^2)^n.$$ I don’t immediately see a way to get from this to Stirling’s formula, though.

I am reminded of this post on Joyal’s isomorphism between doubly rooted labeled trees and endomorphisms. https://golem.ph.utexas.edu/category/2019/12/avisualtellingofjoyals_pro.html

A theorem of Harer and Zagier asserts that the orbifold Euler characteristic of the moduli space $\M_{g,1}$ of genus $g$ smooth curves with one marked point equals $$ζ(1 − 2g) = − B_{2g}/2g ,$$ where $ζ$ is the Riemann zeta function and $B_{2g}$ is a Bernoulli number (indexing conventions vary). On the other hand Stirling’s formula comes from the leading term in an asymptotic expansion $$\log \Gamma (z) = z \log z - z + 1/2 \log (2\pi/z) - \sum_{g \geq 1} \zeta(1 - 2g) \frac{z^{1-2g}}{2g-1} + \dots$$ of the Gamma function. Kontsevich [Comm. Math. Phys. 92] and Mirzakhani [JAMS 07] and others have made good use of this by interpreting the orbifold Euler characteristic as a volume. Homotopy theorists will recognize $\zeta(1- 2g)/2$ mod $Z$ as the image of the $J$-homomorphism in the stable homotopy group $\pi^S_{4g-1}$ of spheres. These are widely thought not to be coincidences.

Also check out Section 3 here, which discusses many different articles on Stirling’s formula in the American Mathematical Monthly:

• Jonathan M. Borwein and Robert M. Corless, Gamma and factorial in the Monthly.

Judging by the number of articles in the Monthly on the subject, Stirling’s formula approximating $n!$ for large $n$ is by far the most popular aspect of the $\Gamma$ function. There are “some remarks”, “notes”, more “remarks”; there are “simple proofs”, “direct proofs”, “new proofs”, “corrections”, “short proofs”, “very short proofs”, “elementary” proofs, “probabilistic” proofs, “new derivations”, and (our favourite title) “The (n+1)th proof”.

That should be “(n+1)st”. But I want to check out the probabilistic proofs!

Don’t forget the regularized formula “$\infty!=\sqrt{2\pi}$” coming from $\zeta'(0)=-\log\sqrt{2\pi}$.

(FWIW this can be thought of as the determinant of an operator with eigenvalues 1, 2, 3, … so it arises when regularising certain path integrals…)

H. Suyari (2004) published a paper on a q-algebraic derivation of Stirling’s formula based on a q-version of the factorial (ie. Tsallis statistics in non-extensive versions of thermodynamics).. This generalization may provide you with some additional insights relating to q-generalizations of the Gaussian distribution.

I’d been thinking about Stirling’s formula recently, and then remembered this Café posting. I’ll take this opportunity to record another reasonably simple proof of Stirling’s formula. The stuff in the beginning is really well-known and appears in many calculus textbooks. The end is not something I’ve seen put quite this way before, although the basic plan of attack seems fairly natural to me.

Trapezoidal sum

My own favorite way of remembering Stirling’s formula is by working through a trapezoidal sum estimate for the integral

$$I_n = \int_1^n \log x\; dx = \left. x\log x - x\right|_1^n = n \log n - n + 1 = 1 + \log \frac{n^n}{e^n}$$ where $\log$ is of course the natural logarithm. Partitioning the interval $[1, n]$ into subintervals of unit length, the corresponding trapezoidal estimate is

$$T_n = \frac{\log 1 + \log 2}{2} + \frac{\log 2 + \log 3}{2} + \cdots + \frac{\log (n-1) + \log n}{2} = \log (n!) - \frac1{2} \log n = \log \frac{n!}{\sqrt{n}}.$$

$T_n$ is an underestimate of $I_n$ because the graph of $\log$ is concave down. In other words, the error terms

$$\int_k^{k+1} \log x\; dx - \frac{\log k + \log (k+1)}{2} \qquad (1)$$ are positive, and it follows that the sum of these error terms from $k=1$ to $n-1$, which is $I_n - T_n$, increases with $n$. However, using the error term estimate for the trapezoidal rule, the term (1) is on the order of $1/k^2$ (the maximum value of $|(\log)''(x)| = 1/x^2$ over $[k, k+1]$). Since $\sum 1/k^2$ converges, we see that the increasing sequence $I_n - T_n$ has a uniform upper bound and therefore a limit. Hence

$$\underset{n \to \infty}{\lim} (I_n - 1) - T_n = \underset{n \to \infty}{\lim} \log \frac{n^n \sqrt{n}}{e^n \cdot n!}$$ exists.

We have thus shown that $n! \sim C\frac{n^n \sqrt{n}}{e^n}$ for some constant $C$.

Connection with Wallis’s formula

It remains to see that $C = \sqrt{2\pi}$. By some scaling analysis, $n! \sim C\frac{n^n \sqrt{n}}{e^n}$ leads to

$$\frac{(2n)!}{n! \cdot n!} \sim \frac{2^{2n} \sqrt{2n}}{C n}.$$ We have

$$\frac1{C} \sim \frac{(2n)!}{2^{2n} n! \cdot n!} \cdot \frac{n}{\sqrt{2n}} = \frac{1 \cdot 3 \cdot 5 \cdot \ldots \cdot (2n-1)}{2 \cdot 4 \cdot 6 \cdot \ldots \cdot 2n} \cdot \frac{n}{\sqrt{2n}}\; , \qquad (2)$$

$$\frac{2}{C} \sim \frac{3 \cdot 5 \cdot \ldots \cdot (2n-1)(2n+1)}{2 \cdot 4 \cdot 6 \cdot \ldots \cdot 2n} \cdot \frac1{\sqrt{2n}}\; ,$$ and now by multiplying the previous two lines,

$$\frac{2}{C^2} \sim \frac{1 \cdot 3}{2 \cdot 2} \cdot \frac{3 \cdot 5}{4 \cdot 4} \cdot \ldots \cdot \frac{(2n-1)(2n+1)}{2n \cdot 2n} \cdot \frac1{2}.$$

This last is recognizable as a product appearing in a ‘Wallis formula’ for $\pi$. One way of evaluating its limit is first by rewriting it as

$$\frac1{2} \prod_{k=1}^n \frac{(2k-1)(2k+1)}{(2k)^2} = \frac1{2}\prod_{k=1}^n \left(1 - \frac1{(2k)^2}\right)$$ and then by using the product formula for the sine function,

$$\frac{\sin \pi x}{\pi x} = \prod_{k = 1}^\infty \left(1 - \frac{x^2}{k^2}\right).$$ But the product formula is a whole other story, involving things like the reflection formula for the Gamma function. Let’s not get too sidetracked: remember that the original goal was to connect $C$ directly to the Gaussian integral

$$\int_{-\infty}^\infty e^{-x^2/2}\; dx$$ and that’s what we should do here as well.

Calculation with Gaussian integrals and the Gamma function

Let’s start afresh. The Gaussian integral is a somewhat protean object, and is inevitably connected with the Gamma function, for example with the value $\Gamma(\frac1{2})$. Indeed,

$$\Gamma\left(\frac1{2}\right) = \int_0^\infty x^{-1/2} e^{-x}\; dx = 2\int_0^\infty e^{-u^2} \; du = \int_{-\infty}^\infty e^{-u^2}\; du$$ (substitute $u^2$ for $x$) and here I think we may take it for granted that we know the classic and beautiful proof that

$$\int_{-\infty}^\infty e^{-u^2}\; du = \sqrt{\pi}$$ by squaring both sides and changing to polar coordinates. We can also take for granted the functional equation

$$\Gamma(x+1) = x\Gamma(x).$$ Thinking asymptotically, one may then guess the interpolation

$$\Gamma(x + \frac1{2}) \approx \sqrt{x} \cdot \Gamma(x).$$ Let’s see where this takes us. For large integers $n$, consider

$$\frac{\Gamma(n + \frac1{2})}{\Gamma(n)\sqrt{n}} = \frac{(n-\frac1{2})(n-\frac{3}{2}) \ldots \frac1{2}\Gamma(1/2)}{(n-1)!\sqrt{n}} = \frac{(2n-1)(2n-3) \ldots 1 \cdot \sqrt{\pi}}{(2n-2)(2n-4)\ldots 2 \cdot \sqrt{n}}. \qquad (3)$$ We know the limit of (3) as $n \to \infty$ exists, based on our similar expression for $\frac1{C}$:

$$\frac1{C} = \underset{n \to \infty}{\lim} \frac{1 \cdot 3 \cdot 5 \cdot \ldots \cdot (2n-1)}{2 \cdot 4 \cdot 6 \cdot \ldots \cdot 2n} \cdot \frac{n}{\sqrt{2n}}. \qquad (2)$$

All that remains is to show the limit of (3) is $1$ (as we suspect from our asymptotic guess), because this is equivalent to $C = \sqrt{2\pi}$. For that, all we need to do is show that we get the same limit

$$L = \underset{x \to \infty}{\lim} \frac{\Gamma(x+ \frac1{2})}{\sqrt{x} \cdot \Gamma(x)} \qquad (4)$$ for half-integers $x$ as we do for whole integers $x$, for then

$$\frac{\Gamma(n+1)}{\sqrt{n+(1/2)}\cdot \Gamma(n + (1/2))} \cdot \frac{\Gamma(n + (1/2))}{\sqrt{n}\Gamma(n)} = \frac{n\Gamma(n)}{\sqrt{(n+1/2)n}\cdot \Gamma(n)} = \frac{n}{\sqrt{n^2 + (n/2)}} \to 1$$ implies $L^2 = 1$, whence $L = 1$: exactly what we want.

The crucial fact we need here is the log-convexity of $\Gamma$, i.e., the fact that the function $\log \Gamma(x)$ is convex over $0\lt x \lt \infty$. For the moment we take this for granted. The immediate consequence of log-convexity is that for $x \gt 0$,

$$\Gamma(x+ 1/2)^2 \leq \Gamma(x)\Gamma(x+1)$$ so that

$$\frac{\Gamma(x-1/2)}{\Gamma(x - 1)} \leq \frac{\Gamma(x)}{\Gamma(x-1/2)} \leq \frac{\Gamma(x+1/2)}{\Gamma(x)}.$$ This means that these values for half-integers $x$ are nestled between values for whole integers $x$ on either side. It follows easily that we get the same limit (4) for half-integers $x$ as we do for whole integers. And now we’re done.

Log-convexity of Gamma function

This is actually not very difficult, but the basic fact on which it rests is in some sense profound:

Theorem: Log-convex functions $(a, b) \to [0, \infty)$ form a rig over the rig of nonnegative real numbers, closed under sups.

That log-convex functions are closed under pointwise multiplication is very easy. Another easy fact that can be left as an exercise is that log-convex functions are closed under taking pointwise suprema. But the real kicker is that log-convex functions are closed under pointwise addition!

I’ll just point to the nLab for that. In essence, it boils down to an application of Hölder’s inequality.

Note that for fixed positive $x$, the function $t \mapsto x^{t-1}$ is log-convex (it is log-linear after all). The integral

$$\Gamma(t) = \int_0^\infty x^{t-1} e^{-x}\;dx$$ can be expressed as a sup over Riemann sums

$$\sum_{i=1}^n x_i^{t-1} e^{x_i} \Delta x_i$$ which are positive linear combinations of log-convex functions. It follows from the theorem that $\Gamma(t)$ is log-convex.

One last thing I’ll recall is that $\Gamma$ is uniquely characterized as a real-valued function $f(x)$ defined for all real $x$ except $x = 0, -1, -2, \ldots$ and satisfying log-convexity, the functional equation $f(x+1) = x f(x)$, and $f(1) = 1$. This is the Bohr-Mollerup theorem. Various convexity arguments are put to use with a kind of terrifying vise-like effect, to squeeze down $f$ to a single possibility. Emil Artin made the Bohr-Mollerup characterization the centerpiece of his famous notes on the Gamma function, which may be found here.

Afterword

While I was thinking about this and researching this, I came across some interesting articles. One of them, “From the Wallis formula to the Gaussian distribution” by Mark Reeder, appears to be a write-up for students that is rich in historical notes, for example following the trail of Wallis in his amazing efforts toward “quadrature of the circle”, and connecting up with Euler, Gauss, and even Stieltjes. I didn’t know, for instance, that the symbol $\pi$ was first introduced by a Welsh mathematician named William Jones in 1706, after Wallis had died. (Sheesh, really? What did the old codgers like Newton and Leibniz use?) Wallis himself denoted the quantity we’d call $4/\pi$ by a symbol which, if I had to name it in TeX, I’d call \boxdot! A box with a dot inside it.

Another is a very short proof of Stirling’s formula due to Dan Romik that appeared in the American Mathematical Monthly. Normally I’m not a huge fan of “World’s Shortest Proofs”, like Don Zagier’s “one-line proof” of the two-squares theorem, because (a) usually they look a little show-offy, and (b) they often leave a lot of work to the reader to actually grok the damned thing. Rabbits out of hats and so forth. [Incidentally, a much kinder explanation of Zagier’s proof can be found somewhere on YouTube, maybe by Mathologer; I forget. IIRC it’s about 30 minutes long, but it really gives the intuition, which is indeed lovely.] But anyway, in the case of Romik, I’ll admit that my jaw dropped. To me his proof has a very Euler-MacLaurin feel to it, but very slickly arranged/optimized.

Finally, I’ll mention A Probabilistic Proof of Wallis’s Formula for $\pi$ by Stephen J. Miller, which involves Student’s t- and seems to have something in common with Romik’s set-up, although I’d need to think longer to be clear on that.