From: Morten Welinder <terra_at_gnome.org>

Date: Sat 23 Oct 2004 - 13:51:28 EST

static gnm_float

calculate_loggos (gnm_float traffic, gnm_float circuits) {

R-devel@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sat Oct 23 13:57:10 2004

Date: Sat 23 Oct 2004 - 13:51:28 EST

> Make that 30400 orders of magnitude (natural logs y'know)...

Right. (/me raises hands showing 2.7 fingers.)

> What the devil are you calculating? The probability that a random

*> configuration of atoms would make up the known universe?
*

Not quite. Where you see a cdf for the gamma distribution I see the incomplete gamma function. Same function, different hat. I am using it to compute the Erlang B function ("Grade Of Service"), see

http://www.dcss.mcmaster.ca/~qiao/publications/erlang/newerlang.html

And here is my code for the log version of this. (Link's c==circuit; link's rho==traffic; link's p is a typo for rho.)

static gnm_float

calculate_loggos (gnm_float traffic, gnm_float circuits) {

double f;

if (traffic < 0 || circuits < 1) return gnm_nan; if (traffic == 0) return gnm_ninf; #ifdef CANCELLATION /* Calculated this way we get cancellation. */ f = circuits * loggnum (traffic) - lgamma1p (circuits) - traffic; #else f = (circuits - traffic) + (1 - loggnum (sqrtgnum (2 * M_PIgnum))) - loggnum (circuits + 1) / 2.0 - logfbit (circuits) + circuits * (loggnum (traffic / (circuits + 1)));#endif

return f - pgamma (traffic, circuits + 1, 1, FALSE, TRUE); }

The two #ifdef branches calculate the same thing, but the bottom version suffers a lot less from cancellation. I might still need to consider cancellation in the final subtraction. (Read "double" where the above says "gnm_float" and forget the "gnum" suffixes.)

In the traffic=1e6,circuits=1e5 case I quoted I could use the second formula from the link above instead, but that won't work when the two are close to each other. Sadly I need it there too.

Googling suggests that the canonical reference for this problem is

Temme N M (1987) On the computation of the incomplete gamma functions for large values of the parameters Algorithms for Approximation J C Mason and M G Cox (ed) Oxford University Press.

(This reference from nag's manual.)

Morten

R-devel@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sat Oct 23 13:57:10 2004

*
This archive was generated by hypermail 2.1.8
: Wed 03 Nov 2004 - 22:45:23 EST
*