> 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
______________________________________________
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
