This post is mostly for my own amusement, mainly showcasing some ways to have fun with calculus.
I’ll begin with an argument of my own, for half of Stirling approximation of \(n!\) on \(\mathbb{N}\) (recovering the shape of its asymptotic form easily, and without integration, but missing the precise constant); then I’ll work out the details in the standard Laplace-method argument for estimating \(\Gamma(t + 1)\) as \(t \to \infty\) in \(\mathbb{R}\), in full rigor (but with no originality claimed here).
Last fall, I was helping two of my friends with some homework for MATH 131AH (the analysis I class here at UCLA). As part of one weekly assignment, they were given the following problem:
- Determine the convergence of the series $$\sum_{n \geq 1} \frac{n!}{n^n}.$$
It’s obvious what this problem is trying to suggest, but let’s hold off for a brief moment. This question can be answered in increasing levels of sophistication, starting with analyzing the sequence \(a_n = n! / n^n\) and its variants. For a start, we have $$\frac{a_{n + 1}}{a_n} = \left( \frac{n}{n + 1} \right)^n.$$ The task, of course, is to determine the limiting behavior of this object. Just by replacing \(n\) with \(2\) (or with any fixed \(k\)), we may get an \(a_n \lesssim 1 / n^2\) polynomial decay bound (or \(1 / n^k\)), which is already enough for convergence. But let’s go further.
We may attempt to get an exponential decay bound, by first noting that $$\left( \frac{n}{n + 1} \right)^n \to e^{- 1} < 1$$ as \(n \to \infty\), which means that the \(a_n\) will display exponential decay with \(n\).
Of course, as the 131AH class is an introductory real analysis course, and things like the exponential function, logarithms, and calculus have all not yet been introduced (nor the number \(e\), even), this would require some additional work to convert into an acceptable solution for that class.
It can be shown, however, that \(n \mapsto n^n / (n + 1)^n\) is a decreasing map from \(\mathbb{N}\) to \(\mathbb{Q} \cap (0, \infty)\) (by means of some elementary inequalities), and this shows the exponential decay by completely-legitimate means.
After working on this, I realized this approach of analyzing the sequence of ratios seemed potentially-useful for considering the factorial, and I wanted to see how far this line of argument could go. As it turns out, quite far. Specifically, I found that I could recover this weak form of Stirling’s, with nothing more than some basic differential calculus:
Lemma: There exists a real number \(c > 0\) such that $$\lim_{n \to \infty} \frac{n!}{n^{1 / 2} (n / e)^n} = c,$$ where the limit is being taken in the sense of sequences (maps defined on \(\mathbb{N}\)).
Proof: We begin with the sequence $$a_n = \frac{e^n n!}{n^n n^{1 / 2}},$$ as one might expect after being inspired by the above. We then have that $$\frac{a_{n + 1}}{a_n} = e \left( \frac{n}{n + 1} \right)^{n + 1 / 2},$$ after some algebra.
As a consequence, we have $$a_{N + 1} = \exp(c_N) a_1,$$ where $$c_N = \sum_{n = 1}^{N} 1 + (n + 1 / 2) \ln \left( \frac{n}{n + 1} \right) = \sum_{n = 1}^{N} \ln \left( \frac{a_{n + 1}}{a_n} \right),$$ and we now recall the power series of \(x \mapsto \ln(1 \ – \, x)\) in \((- 1, 1)\): $$- \ln(1 \ – \, x) = \sum_{n = 1}^{\infty} \frac{x^n}{n}.$$
For this, we only need a Taylor expansion to third order: $$- \ln(1 \ – \, x) = x + \frac{x^2}{2} + O(x^3),$$ for \(0 < x < \frac{1}{3}\), and some implicit constant \(< \frac{1}{3} + \frac{1}{4} \frac{1 / 3}{1 \ – \, 1 / 3}\) from casually estimating the tail.
After more algebra, we can reduce the $$1 + (n + 1 / 2) \ln \left( 1 \ – \, \frac{1}{n + 1} \right)$$ summands in the definition of \(c_N\) to be exactly equal to $$\frac{1}{4 (n + 1)^2} + O \left( \frac{n + 1 / 2}{(n + 1)^3} \right),$$ for all \(n \geq 3\).
As a consequence, we get that $$\lim_{N \to \infty} c_N = L = \sum_{n = 1}^{\infty} \ln(a_{n + 1} / a_n)$$ exists, as a real number; it is given by an absolutely-convergent series. Consequently, if we define \(c = a_1 e^L = e^{L + 1}\), it follows that we have $$\lim_{n \to \infty} \frac{n!}{n^{1 / 2} (n / e)^n} = c,$$ as desired.
Finally, we note that by the Cauchy condensation test (summing over dyadic blocks) to bound the tail of the Basel series \(\sum_n n^{- 2}\), we can conclude a more quantitative estimate $$c_N = L + O(1 / N),$$ so that we have the more precise asymptotic $$n! = c e^{O(1 / n)} n^{1 / 2} \left( \frac{n}{e} \right)^n = c (1 + O(1 / n)) n^{1 / 2} \left( \frac{n}{e} \right)^n$$ by calculus, as \(n \to \infty\) in \(\mathbb{N}\).
Remark 1: A drawback of our proof is that the simplicity lies in the determination of the asymptotic form of the factorial, but without any hint as to the value of the constant. The constant \(c\) can be determined by using a certain family of trigonometric integrals, as explained in the article by M. R. Murty and K. Sampath, or from consulting the classic little book of M. Rosenlicht, which uses the Wallis product; it turns out, of course, to be \(c = (2 \pi)^{1 / 2}\). Dan Romik has a separate, independent argument in Am. Math. Monthly for determining the constant, which I personally find to be very neat and the most elegant of all the proofs I’ve seen for this fact.
Remark 2: At the moment, I can’t see a clear way of extending the simple differential calculus argument to noninteger values (i.e., the \(\Gamma\)-function statement of Stirling, where the limit is taken as \(s \to \infty\) in \(\mathbb{R}\) instead of \(n \to \infty\) in \(\mathbb{N}\)). There is also the complex variant of the Stirling approximation, taking \(z \to \infty\) in \(\mathbb{C}\), but that is even farther beyond me.
Remark 3: One can probably argue this is a trivial result in the grand scheme of things, and something like this would normally not go on the blog. I just felt compelled to write up this argument because, first, even with the restriction to natural numbers and a lack of an exact \(c\), I believe the weak form is still useful;1 and, second, simply because many proofs of the Stirling approximation, just to reach the same weak fragment we show here, give more-complicated arguments.
Remark 4: The first proof I ever saw was in Baby Rudin, and that argument, with its slick change-of-variables and manipulations of integrals, felt exceptionally unmotivated. (The proofs of A. J. Coleman, W. Feller, P. Diaconis & D. Freedman, J. M. Patin, C. L. Frenzen, and R. Michel should also be mentioned; cf. the bibliography and discussion in Borwein & Corless. The Am. Math. Monthly seems to have an inordinate fondness for “elementary” proofs of the Stirling approximation.) Other arguments I have seen involve Euler-Maclaurin summation, or else Taylor series approximations for \(\ln\) combined with a large amount of estimation and unusually-many inequalities for the logarithm, simply to recover the same result we did, with constant evaluated separately. Few of these arguments have felt very intuitive.
Remark 5: Putting my neck out, I’d like to cautiously express the possibility that this argument might, perhaps, be slightly new, as it uses no integration whatsoever, and no inequalities except the remainder from a quadratic Taylor approximation. Indeed, the standard way I have seen the \(\mathbb{N}\)-Stirling approximation proved is by discretizing and approximating $$\ln(N!) = \sum_{n = 1}^{N} \ln(n) = \int_{1}^{N + 1} \ln(t) \, dt + \epsilon_N$$ as an integral, and finding a way to control the error (through Euler-Maclaurin or by estimating the Riemann sums rigorously, often with lots of inequalities on the function \(\ln\)). By contrast, we are simply reducing matters to the tail of a sum, which just requires a suitably-precise evaluation of the logarithm near \(1\) and some algebra; only differentiation is used. In particular, we make no substitutions and draw on no specialized inequalities, or make any decisions on how to discretize as Riemann sums.
A separate argument, common for explaining the presence of the Gaussian, is to argue as follows: for \(t > 0\) in \(\mathbb{R}\), we have $$\Gamma(t + 1) = \int_{0}^{\infty} y^t e^{- y} \, dy = \int_{0}^{\infty} e^{- t (y / t \ – \, \ln(y))} \, dy = \int_{0}^{\infty} e^{- t (y \ – \, \ln(y))} e^{t \ln(t)} t \, dy,$$ and we reduce this to $$t (t / e)^t \int_{0}^{\infty} e^{- t (y \ – \, \ln(y) \ – \, 1)} \, dy = t (t / e)^t I(t),$$ with the understanding that if we can show the integral behaves asymptotically $$I(t) = (1 + o(1)) (2 \pi / t)^{1 / 2}$$ as \(t \to \infty\), then we will have obtained the desired Stirling approximation.
At this point, many blogs and treatments I’ve seen will invoke Laplace’s method for exponential integrals; some will mention the intuition that at its critical points, a nondegenerate phase function \(\Phi\) means the integrand in \(\int e^{- \lambda \Phi}\) appears like a quadratic. While this forms a good basis for intuition behind the square-root asymptotic (the value of the parameter is the inverse width of a Gaussian), it should be made formal. The Wikipedia proof is rather muddled on handling the infinite interval case, as well as generally not very concise. So, let’s do this rigorously: we’ll give a full, precise argument that we can trust.
(Also, in 247B, our main focus was on oscillatory integrals and sums, where the task was extracting cancellation from a rapidly-changing imaginary phase rather than decay from a real-valued negative phase. This task gives me an excuse to work out some of the details of Laplace’s method on my own.)
We begin with the phase function $$\Phi(y) = y \ – \, \ln(y) \ – \, 1,$$ which has derivatives \(\Phi'(y) = 1 \ – \, 1 / y\), \(\Phi^{\prime \prime}(y) = 1 / y^2\). Thus we see that \(y = 1\) is the unique global minimum for \(\Phi\) on \((0, \infty)\): it is strictly decreasing for \(y < 1\), and strictly increasing for \(y > 1\). Moreover, we have that \(\Phi(1) = 0\), \(\Phi'(1) = 0\), and \(\Phi^{\prime \prime}(1) = 1\) (of course, \(\Phi\) is in fact convex).
Then at \(y = 1\), we perform the exact second-order Taylor expansion $$\Phi(1 + h) = \Phi(1) + h \Phi'(1) + \frac{h^2}{2} \int_{0}^{1} 2 \Phi^{\prime \prime}(1 + \theta h) (1 \ – \, \theta) \, d\theta,$$ for \(|h| \ll 1\) small, and note that we have that $$\Phi(1 + h) = \frac{1}{2} (v(h))^2, \qquad v(h) = h \left( \int_{0}^{1} 2 \Phi^{\prime \prime}(1 + \theta h) (1 \ – \, \theta) \, d\theta \right)^{1 / 2}.$$
We have $$v'(h) = \left( \int_{0}^{1} 2 \Phi^{\prime \prime}(1 + \theta h) (1 \ – \, \theta) \, d\theta \right)^{1 / 2} + h \frac{\int_{0}^{1} 2 \Phi^{\prime \prime \prime}(1 + \theta h) \theta (1 \ – \, \theta) \, d\theta}{2 \left( \int_{0}^{1} 2 \Phi^{\prime \prime}(1 + \theta h) (1 \ – \, \theta) \, d\theta \right)^{1 / 2}},$$ and the point is that we can read off from this expression that \(v\) is well-defined and \(C^1\) in a small neighborhood of \(0\). Moreover, substituting \(h = 0\) gives \(v'(0) = 1\).
By the most elementary variant of the inverse function theorem, we then have that \(v\) is a \(C^1\) diffeomorphism of open intervals \((- \epsilon, \epsilon) \to (- c_{\epsilon}^{-}, c_{\epsilon}^{+})\) for some \(0 < \epsilon \ll 1\), and \((v^{- 1})'(v(0)) = (v^{- 1})'(0) = 1\).
We use this to split the exponential integral \(I(c)\) into three pieces: we first take a smooth cutoff \(\eta \in C^{\infty}(\mathbb{R}; [0, 1])\) supported in \((- 1, 1)\), identically \(1\) on \((- \frac{1}{2}, \frac{1}{2})\). Then we define the bulk term $$\int_{\mathbb{R}} \eta((y \ – \, 1) / \delta) e^{- t \Phi(y)} \, dy = \int_{- \delta}^{\delta} \eta(y / \delta) e^{- t \Phi(y + 1)} \, dy,$$ and the remaining terms $$\int_{0}^{\infty} (1 \ – \, \eta((y \ – \, 1) / \delta)) \eta_0(y) e^{- t \Phi(y)} \, dy,$$ and $$\int_{0}^{\infty} (1 \ – \, \eta((y \ – \, 1) / \delta)) (1 \ – \, \eta_0(y)) e^{- t \Phi(y)} \, dy,$$ where \(\eta_0\) is a bump function that helps us localize to \((0, 1 \ – \, \delta / 2)\) (and thus \(1 \ – \, \eta_0\) plays the same role for \((1 + \delta / 2, \infty)\)).
For the bulk term, we will fix a \(\delta > 0\) sufficient to ensure \((- \delta, \delta)\) is contained in the subdomain \((- \epsilon, \epsilon)\) where \(v\) actually acts as a diffeomorphism. This will ensure that we can change variables $$\int_{- \delta}^{\delta} \tilde{\eta}(h) e^{- t v(h)^2 / 2} \, dh = \int_{v(- \delta)}^{v(\delta)} \rho(s) e^{- t s^2 / 2} \, ds,$$ where $$\rho(s) = (\tilde{\eta} \circ v^{- 1})(s) \frac{1}{v'(v^{- 1}(s))} = (\tilde{\eta} \circ v^{- 1})(s) (v^{- 1})'(s),$$ and we will note that where \(s\) is such that \(v^{- 1}(s)\) exceeds \(\delta\) in magnitude, \(\rho\) vanishes; as a consequence, provided \(0 < \delta \ll \epsilon\), depending on the properties of our diffeomorphism, we can ensure that everything is well-defined on \(\mathbb{R}\), and the domain of integration is just the whole real line. Then this is equal to $$\int_{\mathbb{R}} e^{- t s^2 / 2} \, ds + \int_{\mathbb{R}} (\rho(s) \ – \, 1) e^{- t s^2 / 2} \, ds,$$ and we can write this using the compact support of \(\rho\) and a mean-value bound at the origin as $$\left( \frac{2 \pi}{t} \right)^{1 / 2} + O \bigg( \int_{|s| \lesssim 1} |s| e^{- t s^2 / 2} \, ds \bigg).$$
The error term is given by $$\lesssim \int_{|s| \lesssim (t / 2)^{1 / 2}} \frac{1}{t^{1 / 2}} |s| e^{- s^2} \frac{1}{t^{1 / 2}} \, ds = O(1 / t)$$ for \(t \gg 1\), which will suffice. (Though if we wish to be even more precise, we can note that the diffeomorphism \(v\) we constructed above is actually \(C^{\infty}\), meaning \(\rho\) is smooth, and so for any large \(N\), we can take a Taylor expansion of $$\rho(t) \ – \, 1 = \sum_{1 \leq i < N} \frac{1}{i!} t^i \rho^{(i)}(0) + t^N \rho_N(t)$$ at the origin. The integral of any Gaussian times a pure monomial can be computed exactly, with all the odd terms vanishing by symmetry; the coefficients \(\rho^{(i)}(0)\) all follow deterministically from knowing \(v^{- 1}\) and its derivatives at the origin, which depend simply on \(v\) and its derivatives. In principle, we could develop an asymptotic series in \(t^{- 1}\) for the coefficient up to arbitrary degree. More concretely, we have that the above decay bound should actually be \(O(1 / t^{3 / 2})\).)
To deal with the remaining terms, we first note that $$\int_{0}^{\infty} (1 \ – \, \eta((y \ – \, 1) / \delta)) (1 \ – \, \eta_0(y)) e^{- t \Phi(y)} \, dy$$ has decay faster than any polynomial, in the \(t \to \infty\) limit. Indeed, we claim $$\Phi(y) \geq c y$$ for any \(y \geq 1 + \delta / 2\), for some small \(c > 0\) depending on \(\delta > 0\). To see this, we have that $$\frac{d}{d y} \left( \frac{\Phi(y)}{y} \right) = \ – \frac{1}{y^2} + \frac{\ln(y)}{y^2} + \frac{1}{y^2} > 0$$ for \(y \geq 1 + \delta / 2\), and so \(\Phi(y) \geq (\delta / 2 \ – \, \ln(1 + \delta / 2)) y\) for all such \(y\). As a consequence, we can simply bound this by an exponential integral $$\int_{1 + \delta / 2}^{\infty} e^{- c t y} \, dy \lesssim_{\delta} \frac{1}{t} e^{- c t},$$ which certainly decays faster than any power of \(1 / t\).
The remaining term $$\int_{0}^{\infty} (1 \ – \, \eta((y \ – \, 1) / \delta)) \eta_0(y) e^{- t \Phi(y)} \, dy,$$ can be handled in a similar manner. As \(\Phi(y) = y \ – \, \ln(y) \ – \, 1\), we saw that \(\Phi\) is strictly decreasing on \((0, 1)\), and so \(\Phi(y) \geq \Phi(1 \ – \, \delta / 2) = \, – \delta / 2 \ – \, \ln(1 \ – \, \delta / 2) > 0\). This is strictly positive, and so we get a decay bound of \(\leq e^{- c t}\).
Adding all of these contributions together, we get an estimate on \(I(t)\) and after some rewriting, the desired Stirling’s approximation for the \(\Gamma\)-function follows: $$\Gamma(t + 1) = (2 \pi)^{1 / 2} (1 + O(1 / t)) t^{1 / 2} (t / e)^t$$ as \(t \to \infty\) in \(\mathbb{R}\) (no longer restricted to natural numbers). If one prefers, one can divide by \(t\) and modify the RHS to read \(t^{- 1 / 2}\) to get the corresponding asymptotic for \(\Gamma(t)\) proper.
- I think it’s safe to say that in analysis, we often make do with equivalence up to multiplicative constants; more facetiously, we can treat all values \(0 < c < \infty\) as \(c \sim 1\) in many arguments (until, for whatever reason, we can’t). ↩︎
Leave a Reply