Re: [Rd] pgamma discontinuity (PR#7307)

From: Morten Welinder <terra_at_gnome.org>
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