Stirling Approximations for Log[n!]

from Statistical Physics using Mathematica
© James J. Kelly, 1996-2001

Initialization

ClearAll["Global`*"];
Off[General::spell, General::spell1]

$DefaultFont = {"Times", 12};
$TextStyle = {FontFamily -> "Times", 
          FontSize -> 12};

Needs["Graphics`Master`"]

Obtain simple approximation using a graphical argument

Recognizing that  Log[n!] can be expanded as

(∑_(k = 1)^n) Log[k]

we can approximate the sum as an integral.  The accuracy of this approximation can be gauged by comparing the histogram below with the area under the smooth curve described by the integrand or summand.  Below we construct a function displayApprox[n] which compares the sum and integral graphically.  By displaying these plots for increasing arguments, we find that the integral is a very good approximation to the sum even for relatively small arguments.

curve[nmax_] := 
        Plot[Log[n], {n, 1, nmax}, 
          PlotRange -> {{0, nmax + 1}, Automatic}, Frame -> True, 
          FrameLabel -> {"n", "Log[n]"}, 
          PlotLabel -> "Stirling Approximation", 
          DisplayFunction -> Identity];
steps[nmax_] := 
        BarChart[Table[Log[n], {n, 1, nmax}], BarSpacing -> -0.2, 
          DisplayFunction -> Identity];
displayApprox[nmax_] := 
        Show[{curve[nmax], steps[nmax]}, 
          DisplayFunction -> $DisplayFunction];

{displayApprox[10], displayApprox[100]};

Thus, it is clear that Stirling's  approximation becomes quite accurate as n becomes large.  To test this accuracy directly, we evaluate the percent difference between the area under the smooth curve and the actual discrete sum.  This difference is plotted as a function of n below.

logStirling[n] = (∫_1^n) Log[x] \[DifferentialD]x

\(\(LogPlot[
        Evaluate[\(100. \ Abs[Log[\(n!\)] - logStirling[n]]\)\/Log[\(n!\)]], {n, 5, 1000}, Frame -> True, PlotLabel -> "log N!", 
        FrameLabel -> {"N", "percent error"}];\)\)

Asymptotic expansions

In this section we develop several improved versions of Stirling's approximation using the method described in Appendix B of the book by Pathria.

To improve upon the simplest approximation to n!, it is useful to examine the integrand of the gamma function in some detail.  First we define that integrand and verify that we do obtain the expected integral from it.

integrand[n_] = (E^(-x)) (x^n)

(∫_0^∞) integrand[n] \[DifferentialD]x

By inspecting this integrand either graphically or analytically, we find that it resembles a Gaussian centered upon n with width (Sqrt[n]).  Therefore, we make a change of variable to perform an expansion about the Gaussian approximation.

Plot[integrand[10], {x, 0, 20}];

Series[Log[integrand[n]], {x, n, 2}] // Simplify

logApprox[n_, 
        m_] := Log[
              integrand[n] /. x -> n + u (Sqrt[n])] // 
            Series[#, {u, 0, m}] & // Simplify // Normal

The lowest  order approximation to the log of the integrand yields the familiar Stirling approximation as follows.  Notice that it is customary to extend the range of integration to include the entire negative range of u.   Although the plot above demonstrates that the contribution made by negative u is negligible for large n, the student might find it worthwhile to investigate this issue quantitatively.

term0 = logApprox[n, 2]

stirling0 = (Sqrt[n]) (∫_(-∞)^∞) (E^term0) \[DifferentialD]u

We can now develop a correction by expanding the rest of the integrand in powers of 1/(Sqrt[n]).  To expedite this operation, we define an intermediate variable, y  =  1/(Sqrt[n]), and expand with respect to y.  Once we have that expansion, we undo the substitution.   However, in performing these expansions, we must be careful to retain all terms of a given order consistently.  We define term[m] to be the correction to the argument of the exponential with m being the highest power of u.   The highest power of y is then  m - 2.

logApprox[n, 12]

term[
        m_] := logApprox[n, m] - term0 /. n -> 
              (1/(y^2)) // PowerExpand // Simplify

term[6]

Next we define a correction factor of order m, where now m is the highest power of y retained when the exponential is expanded in powers of y.  To improve the efficiency of integration over u, it helps to collect terms with common powers of u.

correctionFactor[
        m_] := Series[(E^term)[m + 2], {y, 0, 
                m}] /. y -> (1/(Sqrt[n])) // Normal // 
        Collect[#1, u] &

correctionFactor[3]

Finally, we obtain a Stirling formula of order m by integration of the Gaussian approximation times the expansion of the remaining exponential factor.  Here m is the highest power of 1/n in the asymptotic expansion.  The range of integration is once again extended to ±∞.

stirlingFormula[
        m_] := (Sqrt[n]) (∫_(-∞)^∞) (E^term0) correctionFactor[2 m + 1] \[DifferentialD]u // 
        Collect[#, n] &

For example, consider the (7^th) order formula.

stirling7 = stirling0 (stirlingFormula[7] / stirling0 //Simplify)

The expansion up to (n^(-4)) reproduces the formula quoted by Pathria, but using Mathematica we can easily carry the analysis to higher order.  Finally, we derive an expansion for Log[n!] to the same power of n.  We separate the common factor to assist Mathematica in handling an essential singularity.

PowerExpand[Log[stirling0]] + 
          Normal[Series[Log[stirling7/stirling0], {n, ∞, 7}]] // 
        Simplify // Collect[#, n] &

Obviously, these formulas can  be derived much more quickly and reliably using Mathematica than by pencil and paper!

Evaluate accuracy of Stirling formulas for Log[n!]

The accuracy of several of the lowest-order Stirling formulas is evaluated numerically.  The results are displayed as plots of the percent deviation from the exact formula versus n.

Stirling formulas

logStirling1 = {Log[n_!] -> 
            n Log[n] - n};
logStirling2 = {Log[n_!] -> 
            n Log[n] - n + (1/2) Log[2 π n]};
logStirling3 = {Log[n_!] -> 
            n Log[n] - n + (1/2) Log[2 π n] + 
              (1/(12 n))};

Accuracy for Log[n!]

logN[n_] = 
      Log[n!]; (logN1[n_] = 
          Log[n!] /. logStirling1;)
logN2[n_] = Log[n!] /. logStirling2;
logN3[n_] = 
        Log[n!] /. logStirling3;

errorLogN1[n_] := N[100 (logN1[n]/logN[n] - 1)]
errorLogN2[n_] := N[100 (logN2[n]/logN[n] - 1)]
errorLogN3[n_] := N[100 (logN3[n]/logN[n] - 1)]

LogPlot[
        Evaluate[Abs[{errorLogN1[n], errorLogN2[n], errorLogN3[n]}]], {n, 5, 
          100}, Frame -> True, PlotLabel -> "log N!", 
        FrameLabel -> {"N", "percent error"}];

The third formula, which still only order (1/n), is already quite accurate even for relatively small n and becomes incredibly accurate as n becomes large.

Problems

Effect of extending integration range

Study the error made by extending the integration range for u from to (0,∞) to (-∞,∞).