stirlerr.m

function delta = stirlerr(n)
%STIRLERR Error of Stirling-De Moivre approximation to n factorial.
%    DELTA=STIRLERR(N) computes
%        gammaln(n+1) - (0.5*log(2*pi*n)+n*log(n)-n)
%    This function is used in C. Loader, "Fast and Accurate Calculations of
%    Binomial Probabilities", July 9, 2000, to compute binomial and Poisson
%    probabilities.
%
%    DELTA is approximated as
%        delta(n) = 1/(12*n) - 1/(360*n^3) + 1/(1260*n^5) + O(1/n^7)
%    For small values of "n", the pre-calculated values for delta(n) are
%    used.

%   Copyright 2010 The MathWorks, Inc.


delta = zeros(size(n),'like',n);
nn = n.*n;

% Define S0=1/12 S1=1/360 S2=1/1260 S3=1/1680 S4=1/1188
S0 = 8.333333333333333e-02;
S1 = 2.777777777777778e-03;
S2 = 7.936507936507937e-04;
S3 = 5.952380952380952e-04;
S4 = 8.417508417508418e-04;

% Define delta(n) for n<0:0.5:15
sfe=[                    0;...
     1.534264097200273e-01;...
     8.106146679532726e-02;...
     5.481412105191765e-02;...
     4.134069595540929e-02;... 
     3.316287351993629e-02;...
     2.767792568499834e-02;...
     2.374616365629750e-02;...
     2.079067210376509e-02;...
     1.848845053267319e-02;...
     1.664469118982119e-02;...
     1.513497322191738e-02;...
     1.387612882307075e-02;...
     1.281046524292023e-02;...
     1.189670994589177e-02;...
     1.110455975820868e-02;...
     1.041126526197210e-02;...
     9.799416126158803e-03;...
     9.255462182712733e-03;...
     8.768700134139385e-03;...
     8.330563433362871e-03;...
     7.934114564314021e-03;...
     7.573675487951841e-03;...
     7.244554301320383e-03;...
     6.942840107209530e-03;...
     6.665247032707682e-03;...
     6.408994188004207e-03;...
     6.171712263039458e-03;...
     5.951370112758848e-03;...
     5.746216513010116e-03;...
     5.554733551962801e-03];

k = find(n<=15);
if any(k)
    n1 = n(k);
    n2 = 2*n1;
    if all(n2==round(n2))
        delta(k) = sfe(n2+1);
    else
        lnsr2pi = 0.9189385332046728;
        delta(k) = gammaln(n1+1)-(n1+0.5).*log(n1)+n1-lnsr2pi;
    end
end
k = find(n>15 & n<=35);
if any(k)
    delta(k) = (S0-(S1-(S2-(S3-S4./nn(k))./nn(k))./nn(k))./nn(k))./n(k);
end
k = find(n>35 & n<=80);
if any(k)
    delta(k) = (S0-(S1-(S2-S3./nn(k))./nn(k))./nn(k))./n(k);
end
k = find(n>80 & n<=500);
if any(k)
    delta(k) = (S0-(S1-S2./nn(k))./nn(k))./n(k);
end
k = find(n>500);
if any(k)
    delta(k) = (S0-S1./nn(k))./n(k);
end

end