Understanding a particular approximation for Stirling's number of the second kind

Question

I noticed and employed (without a problem) an approximation for Stirling's number of the second kind found on Wikipedia (http://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind), in particular this approximate expression: (http://upload.wikimedia.org/math/8/a/f/8afeed15dd40295320cf974418895cc2.png):

$n \brace k$$\space \approx \frac{\sqrt{n-k}}{\sqrt{n(1-G)}G^k(v-G)^{n-k}} (\frac{n-k}{e})^{n-k}$ ${n}\choose{k}$$,\space \space \space \forall{k}$ s.t. $1<k<n$

Where:

$G = -W_0(-ve^{-v}) \space \space \space$ ($W_0$ being the main branche of Lambert W function)

$v = (\frac{n}{k})$

However, upon looking up the stated references [Ref. 13 and 14 on the Wikipedia page]:

1. W. E. Bleick and Peter C. C. Wang, Asymptotics of Stirling Numbers of the Second Kind, AMS Vol.42 No.2, 1974.

2. N. M. Temme, Asymptotic Estimates of Stirling Numbers, STUDIES IN APPLIED MATHEMATICS 89:233-243 (1993), Elsevier Science Publishing.

It's not clear to me where this exact expression comes from (though approximate expressions in the same vein are presented in both papers that rely on the Lambert W function). Can anyone help me out in terms of understanding where the exact form of the expression on Wikipedia comes from and/or what liberties are being taken by the Wikipedia author of the expression?

For convenience, here's the expression (http://upload.wikimedia.org/math/8/a/f/8afeed15dd40295320cf974418895cc2.png) in Mathematica format:

v = n/k;

G = -ProductLog[0, -v*E^(-v)];

stirlingApprox = (n - k)^(1/2)/((n*(1 - G))^(1/2)G^k(v - G)^(n - k))*((n - k)/E)^(n - k)*Binomial[n, k];

Note that "StirlingS2[n, k]" directly gives the Stirling number of the second kind in Mathematica.

The approximation appears to satisfy the stated relevant error of $\approx (\frac{0.06}{n})$, however I haven't checked this carefully yet.

---

As an update --- "The approximation appears to satisfy the stated relevant error of $\approx (\frac{0.06}{n})$, however I haven't checked this carefully yet."

I've now checked $n \brace k$ for $n \in [2, 1000]$ and $\forall k < n$ for each $n$. I found: (1) that the value for $\delta = |(1 - \frac{(exact)}{(approximation)})|$ decreases monotonically with increasing $n$ (regardless of the value of $k < n$); and (2) that for each value of $n$, the value of $\delta$ depends on $k$ in a particular, though consistent manner for most values of $n$. More specifically, there seems to be one local maxima for the smaller values of $k$, a global maxima appears near the mid-range values of $1<k<n$, and then we see that $\delta$ decreases monotonically after the global maxima with increasing $k$.

The largest value of the error we encounter is $\delta_{max} = 0.040932861725805439857312709933506101086741032506081...$ for $(n,k) = (2,1)$.

Now, in terms of whether $\approx (\frac{0.06}{n})$ is a good upper bound for $\delta$, actually it's not conservative enough -- we can notice that all values of $\delta$ very slightly exceed this bound. So perhaps $\approx (\frac{0.066}{n})$ is a more appropriate upperbound, as it is satisfied by all values of $\delta$ for $n \geq 6$.

Answer 1

I cannot provide you with a reference, but I can sketch with my own derivation (I have a draft somewhere...).

First, let $X$ be a zero-truncated multinomial $(n,k)$ random variable (we throw $n$ ball in $k$ urns and discard the results having empty urns). Its probability mass function is

$$P_X({\bf x})=\frac{1}{k! \, {n \brace k}} \frac{n!}{\prod_{i=1}^k x_i! } \left[ x_i \ge 1 \right] \left[ \sum x_i = n\right] \tag{1}$$

Further, let $Z$ be $k$ iid random variables with truncated Poisson distribution: $$P_z(z_i)=(e^\lambda-1)^{-1} \frac{\lambda^{z_i}}{z_i!} [z_i \ge 1] \tag{2}$$ $$P_Z({\bf z})=\prod_{i=1}^k P_z(z_i)$$

This second model is equivalent to the first if conditioned to the event $\mathcal{E}\equiv \sum z_i=n$. That is:

$$ P_X = P_{Z \mid \mathcal{E}} = \frac{P_{\mathcal{E} \mid Z} }{P_\mathcal{E}} P_Z \tag{3}$$

Then, by picking any configuration inside the event set $\mathcal{E}$, we can write:

$${n \brace k} \frac{k!}{n!} = P_\mathcal{E} \frac{(e^\lambda-1)^k}{\lambda^n} \tag{4}$$

Notice that this equation is exact (see also Springael, Johan & Nieuwenhuyse, Inneke. (2006). On the sum of independent zero-truncated Poisson random variables, for a different derivation).

We now approximate $P(\mathcal{E})$ by using the CLT (we could also get higher terms using an Edgeworth expansion, but a first order approximation is enough for this question).

Let's set $E[z_i]=E[x_i]= \frac{n}{k} \triangleq v$ [*]. Then we get the equation $$\frac{\lambda}{1-e^{-\lambda}}=v\tag{5}$$ which implicitly defines $\lambda$ as function of $v$, say $$\lambda= h(v)=v+ W(-v \, e^{-v}) \tag{6}$$

Equivalently, we can set $g=v-\lambda$, with $g$ being the solution in $(0,1)$ of $ \displaystyle g \, e^{-g} =v \, e^{-v} $.

Because $E(\sum z_i)=n$ and the variance of a Poisson truncated is $v \left(1+\lambda-v\right) $, the CLT gives:

$$\begin{align} {n \brace k} &\approx \frac{n!}{k!} \frac{1}{\sqrt{2\pi n (1+\lambda - v )}}\frac{(e^\lambda -1)^k}{\lambda^n} \\ &={n \choose k} \frac{(n-k)!}{\sqrt{2\pi n (1+\lambda - v)}} \frac{1}{\lambda^{n-k} (v-\lambda)^k} \\ &={n \choose k} \frac{(n-k)!}{\sqrt{2\pi n (1-g)}} \frac{1}{ g^k (v-g)^{n-k}} \\ \tag{7} \end{align}$$

Replacing in the last one $(n-k)!$ by its Stirling approximation we get the same expression you posted.

This also shows the order of the approximation: for increasing $k$, keeping $v=n/k$ constant, we have

$$ \frac{{\rm approx}}{\rm exact}=1+O(k^{-1})$$

Again, using the first term of the Edgeworth expansion (involving third and fourth moments of the Truncated Poisson, see eg here) we can obtain the coefficient of the $1/k$ term (or equivalently $1/n$) and attain a more refined $O(k^{-2})$ approximation. The $1/k$ coefficient would depend on $v=n/k$, though.

[*] Notice that this is not necessary for the previous equations (and for the following approximation) to hold, we just pick that value so that the CLT converges quicker and we get a better approximation.

---

Edit: A sample code in JS for computing $\lambda = h(v)$:

function lambda_stirling(v) {
  // solves v = lambda/(1-exp(-lambda)) with newton raphson
  // we use u = v - lambda , so that u = v * exp(u-v)
  let u = 0.0;
  let du = v * Math.exp(-v);
  do {
    u += du;
    du = (v * Math.exp(u - v) - u) / (1.0 - u);
  } while (Math.abs(du) > 1e-15);
  return v - u;
}