Robbins' version of the Stirling approximation

November 2022

The most beautiful formula in mathematics is certainly not the triviality 1+eiπ=01 + e^{i\pi}=0, as claimed by many tasteless mathematicians with no self-esteem, but rather Stirling's formula, stating that n!nnen2πnn!\sim n^n e^{-n}\sqrt{2\pi n}. Quantifying the error of this approximation can be done to any precision using the Euler-MacLaurin complicated formula, but it is often handy to have a simpler estimation; Herbert Robbins[1] found one (original note here) which is surprisingly accurate:

e112n+1<n!nnen2πn<e112n. e^{\frac{1}{12n + 1}} < \frac{n!}{n^n e^{-n} \sqrt{2\pi n}} < e^{\frac{1}{12n}}.

The inequalities are strict and valid for every nn. The proof of this result is given below. For future reference, we can compare (1) with the higher-order expansion: 

n!nnen2πn=exp{112n1360n3+11260n5+O(1n7)}. \frac{n!}{n^n e^{-n}\sqrt{2\pi n}} = \exp \left\lbrace \frac{1}{12n} - \frac{1}{360n^3} + \frac{1}{1260n^5} + O(\frac{1}{n^7}) \right\rbrace.

Proof of (1)

We have ln(n!)=k=1n1ln(k+1)\ln(n!) = \sum_{k=1}^{n-1} \ln(k+1). Robbins' proof consists in a subtle approximation of ln(k+1)\ln(k+1), as seen as the area of the rectangle [k,k+1]×[0,ln(k+1)][k,k+1] \times [0, \ln(k+1)]. This area is equal to: 


Clearly, Tk=(ln(k+1)ln(k))/2T_k = (\ln(k+1) - \ln(k))/2. We also recall that the antiderative of ln(x)\ln(x) is xln(x)xx\ln(x) - x. Consequently, noting Sn=δ2++δn1S_n = \delta_2 + \dotsb + \delta_{n-1},

ln(n!)=I2++In1+T2++Tn1Sn=1nln(x)dx+12ln(k)Sn=nln(n)n+1+12ln(n)Sn.\begin{aligned} \ln(n!) &= I_2 + \dotsb + I_{n-1} + T_2 + \dotsc + T_{n-1} - S_n\\ &= \int_1^n \ln(x)dx + \frac{1}{2}\ln(k) - S_n\\ &= n\ln(n) - n + 1 + \frac{1}{2}\ln(n) - S_n. \end{aligned}

Estimating SnS_n

How big is SnS_n? Well, first we can reckon the δk\delta_k: they are equal to IkI_k minus the area of the trapezoidal approximation of IkI_k, which is (ln(k)+ln(k+1))/2(\ln(k) + \ln(k+1))/2:

δk=kk+1ln(x)dxln(k)+ln(k+1)2=(k+1)ln(k+1)(k+1)kln(k)+kln(k)+ln(k+1)2=1+ln(k+1k)(k+0.5).\begin{aligned}\delta_k &= \int_k^{k+1}\ln(x)dx - \frac{\ln(k) + \ln(k+1)}{2}\\ &= (k+1)\ln(k+1) - (k+1) - k\ln(k) + k - \frac{\ln(k) + \ln(k+1)}{2}\\ &= -1 + \ln\left(\frac{k+1}{k}\right)(k + 0.5). \end{aligned}

Robbin's trick

Now you might want to develop the determinant as ln(1+k1)=k1+k2/2+\ln(1 + k^{-1}) = k^{-1} + k^{-2}/2 + \dots: don't do this. Instead, do the following dark magic: first, note that k+1k=1+x1x\frac{k+1}{k} = \frac{1 + x}{1-x} where x=(2k+1)1x = (2k+1)^{-1}. Then, use the analytic formula

ln((1+x)/(1x))=2=0x2+12+1\ln((1+x)/(1-x)) = 2\sum_{\ell = 0}^\infty \frac{x^{2\ell +1}}{2\ell +1}

so that

δk=1+1x(x+x33+x55+x77+)=x23x45x67=13(2k+1)215(2k+1)417(2k+1)6\begin{aligned} \delta_k &= -1 + \frac{1}{x}\left(x + \frac{x^3}{3} + \frac{x^5}{5} + \frac{x^7}{7} + \dotsb \right)\\ &= - \frac{x^2}{3} - \frac{x^4}{5} - \frac{x^6}{7} - \dotsb \\ &= - \frac{1}{3(2k+1)^2} - \frac{1}{5(2k+1)^4} - \frac{1}{7(2k+1)^6} - \dotsb \end{aligned}

We now use this to bound δk\delta_k below and above.

Upper bound

Since 3,5,7,93,5,7,9\dots are all greater than 3,

δk<13(2k+1)2+13(2k+1)4+13(2k+1)6+=13(2k+1)211(2k+1)2=112k(k+1)=112k(1k1k+1).\begin{aligned} - \delta_k &< \frac{1}{3(2k+1)^2} + \frac{1}{3(2k+1)^4} + \frac{1}{3(2k+1)^6} + \dotsb \\&= \frac{1}{3(2k+1)^2}\frac{1}{1 - (2k+1)^{-2}}\\&= \frac{1}{12k(k+1)}= \frac{1}{12k}\left(\frac{1}{k} - \frac{1}{k+1} \right). \end{aligned}

Lower bound

On the other hand, 3,5,7,93,5,7,9\dots are smaller than 3,32,33,343,3^2, 3^3, 3^4\dots, so

δk>13(2k+1)2+132(2k+1)4+133(2k+1)6+=13(2k+1)211(3(2k+1))2=112((k+1/2)21/12).\begin{aligned} -\delta_k &> \frac{1}{3(2k+1)^2} + \frac{1}{3^2(2k+1)^4} + \frac{1}{3^3(2k+1)^6} + \dotsb \\&= \frac{1}{3(2k+1)^2}\frac{1}{1 - (3(2k+1))^{-2}}\\ &= \frac{1}{12((k+1/2)^2 - 1/12)}. \end{aligned}

It turns out that (k+1/2)21/12<(k+1/12)(k+1/12+1)(k+1/2)^2 -1/12 < (k+1/12)(k+1/12 + 1) (develop both sides then dismiss some terms), so that

δk>112(1k+1/121k+1/12+1).-\delta_k > \frac{1}{12}\left( \frac{1}{k+1/12} - \frac{1}{k+1/12 + 1} \right).

Adding up

From the preceding bounds we see that the series δ1+δ2+\delta_1 + \delta_2 + \dotsb is indeed convergent. If ss is its sum, Sn=sk>nδkS_n = s - \sum_{k>n}\delta_k and using again the precedings bounds (they are telescoping), we can estimate:

s112n<Sn<s112n+1 s - \frac{1}{12n}< S_n < s - \frac{1}{12n+1}

Now go back last line of (1). Take exponentials to get n!/nnenn=e1Snn!/n^n e^{-n}\sqrt{n} = e^{1-S_n}. Then, with c=1sc=1-s and the estimate above,

ec+112n+1<n!nnenn<ec+112n. e^{c + \frac{1}{12n+1}}<\frac{n!}{n^n e^{-n}\sqrt{n}}< e^{c + \frac{1}{12n}}.

What about the constant ?

It's obviously not over, since we didn't get the constant cc. This, however, is usually done by another means, typically with the Wallis integral asymptotics as they here.

[1]  Herbert Robbins does not have the fame he deserves. He's the Robbins of the Robbins-Munro algorithm, the Lai-Robbins bound on bandit algorithms, the backward algorithm for the secretary problem: three essential contributions to different domains of statistics and computational mathematics. His book What is Mathematics ? with Courant is a masterpiece.