# Erlang C for Large Numbers

The formula for Erlang C gives the probability of not having your call answered immediately at a call center, so a number between 0 and 1:

$${E_C} = {1 \over {1 + \left( {1 - p} \right){{m!} \over {{u^m}}}\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}} }}$$

p (agent occupancy) is a number between 0 and 1. With m (agents on hand) and u (call intensity, m>u) large (e.g., u = 143, m = 144), 64-bit floating-point arithmetic overflows -- so I'm trying to figure out how to calculate the denominator in log space, i.e.,

$${E_C} = {1 \over {1 + \left( {1 - p} \right)\exp \left( {\ln \left( {{{m!} \over {{u^m}}}\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}} } \right)} \right)}}$$

... and got this far:

$$\displaylines{ \ln \left( {{{m!} \over {{u^m}}}\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}} } \right) = \ln m! - \ln {u^m} + \ln \left( {\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}} } \right) \cr = \ln m! - m\ln u + \ln \left( {\sum\limits_{k = 0}^{m - 1} {{u^k}} } \right) - \ln \left( {\sum\limits_{k = 0}^{m - 1} {k!} } \right) \cr}$$

The first two terms are no problem (Stirling's Approximation for the first), but am stuck on the last two. I looked at the log identity for a sum

$$\ln \sum\limits_{i = 0}^N {{a_i}} = \ln {a_0} + \ln \left( {1 + \sum\limits_{i = 1}^N {{{{a_i}} \over {{a_0}}}} } \right)$$

... but it doesn't get me anywhere.

This is an exercise for fun, and I am not a mathematician, but would welcome a suggestion.

In searching, I saw this same question asked at https://forums.adobe.com/thread/841586, but there was a deathly silence. I'm optimistic the person was just asking on the wrong forum.

EDIT: In the alternative, might there be a way to restructure the summation so that the m!/u^m term (a little tiny number) can be calculated incrementally to moderate the gargantuan u^k/k! terms? In that context, I reckon this is more of a numerical analysis or programming problem than math. Would Stack Overflow be a more appropriate forum?

• There is a mistake in those equations: $\ln \left( {\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}} }\right)$ is not $\ln \left( {\sum\limits_{k = 0}^{m - 1} {{u^k}} } \right) - \ln \left( {\sum\limits_{k = 0}^{m - 1} {k!} } \right)$ in general. Computing $\sum_{k=0}^{m-1}\exp(\ln(m!)-\ln(k!)+(k-m)\ln(u))$ seems like a good approach to me. Depending on your environment you can probably use a built-in log gamma like Matlab's gammaln; $\ln(m!)=\mathrm{gammaln}(m+1)$. (Though in practice you'd probably use a lookup table.) – Dap Sep 27 '17 at 14:37
• @Dap: D'oh, of course it's not, thank you. Looking at your suggestion now. – shg Sep 27 '17 at 17:08
• I've been unsuccessful in getting that to work – shg Sep 27 '17 at 20:52
• Oh, happy dance, thank you! I would be grateful if you could briefly explain how you derived that. Also, if anyone is interested into how to implement that in an Excel formula, I would be happy to post it. – shg Sep 27 '17 at 21:03
• Ah: I see it:$$\ln \left( {{{m!} \over {{u^m}}}\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}} } \right) = \ln \left( {m!\sum\limits_{k = 0}^{m - 1} {{{{u^{k - m}}} \over {k!}}} } \right) = \ln m! - \ln k! + \left( {k - m} \right)\ln u$$And Excel also has a GAMMALN function, so that was perfect. Thank you again. – shg Sep 27 '17 at 22:20

The question author shg has come up with a practical solution and posted a VBA function in the Adobe thread linked above.

This answer describes an alternate solution in terms of built-in functions not requiring iteration. Unfortunately it doesn't seem to work accurately in Excel, but it should work in R or Matlab.

The sum can be expressed in terms of a cumulative distribution function: \begin{align*} \frac{m!}{u^m}\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}} &=\frac{m!}{u^m}e^uF_{\mathrm{Poisson}}(m-1;u)\\ &=\exp(\ln(\Gamma(m+1))-m\ln(u)+u)F_{\mathrm{Poisson}}(m-1;u) \end{align*}

where $F_{\mathrm{Poisson}}(m-1;u)$ is the probability that a $\mathrm{Poisson}(u)$ variable is at most $m-1$, and $\Gamma(m+1)$ denotes the Gamma function. The $F_{\mathrm{Poisson}}(m-1;u)$ factor is numerically nice because it is not too small, in fact it looks like $$\tfrac 1 2 < F_{\mathrm{Poisson}}(m-1;u)<1\text{ for u<m.}$$

(I don't know a rigorous proof for the lower bound, but by the normal approximation for Poisson variables, for large $u$ the median will be approximately the mean, giving $F_{\mathrm{Poisson}}(m-1;u)\geq F_{\mathrm{Poisson}}(u;u)\to \tfrac 1 2$.)

Excel's POISSON function can be used to evaluate $F_{\mathrm{Poisson}}$, so the whole expression could be expressed as:

EXP(GAMMALN(m+1)-m*LN(u)+u)*POISSON(m-1,u,TRUE)

Unfortunately Excel's POISSON can't be trusted here; my copy of Excel gives POISSON(1e6,1e6,TRUE)=0.863245255, where R gives the more plausible ppois(1e6,lambda=1e6)=0.500266.

• My primary goal was to create a formula to replace my current (agent-limited) formula, and the result was = 1 / (1 + (1 - u / m) * SUMPRODUCT(EXP(GAMMALN(m + 1) - GAMMALN(RowVec(1, m)) + (RowVec(0, m - 1) - m) * LN(u)))), which works great. The RowVec function just generates a simple series in lieu of Excel's awkward INDIRECT(ROW(...)) construct. The workbook is at app.box.com/s/7suleld2atsv7z5z526byvxxk9m9t4zl – shg Sep 28 '17 at 15:22
• In Excel 2003, I get the result you show for Excel's POISSON function. In Excel 2010, both POISSON and POISSON.DIST return the same value as R's function. – shg Sep 28 '17 at 15:37
• So the modified formula works great: = 1 / (1 + (1 - u / m) * EXP(GAMMALN(m + 1) - m * LN(u) + u) * POISSON(m - 1, u, TRUE)) – shg Sep 28 '17 at 15:52
• Updated the function at the linked Adobe thread for this formulation. Dap, thank you again. – shg Sep 28 '17 at 16:29

Since the summation is over $k$, you cannot have $k$ in the result. $$\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}}=\frac{ \Gamma (m,u)}{\Gamma (m)}\,e^u$$ where appears the incomplete gamma function. $$\frac {m!}{u^m}\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}}=m \,e^u\, u^{-m} \,\Gamma (m,u)$$

Have a look here for the calculation of the incomplete gamma function using Excel.

• This may still result in unnecessary overflows since $\Gamma(m,u)$ will be on the scale of $(m-1)!$ for $m\approx u$ – Dap Sep 28 '17 at 11:09
• I understand your comment, and don't doubt that the expression shown in my prior comment is incorrect, but the summation is outside the expression in the resulting formula. Thank you. – shg Sep 28 '17 at 16:34