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

From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>
Date: Sat 23 Oct 2004 - 07:41:55 EST

terra@gnome.org writes:

> for (x = 99990; x <= 100009; x++)
> printf ("%.14g --> %.14g\n", x, pgamma (1000000.0, x, 1.0, 0, 1));
....
> and have no further local changes.
>
> I get...
> src/nmath/standalone> LD_LIBRARY_PATH=. ./test
...
> 99999 --> -669752.66592471
> 100000 --> -669750.36332851
> 100001 --> -599731.36192347 <--- Look at me jump!
> 100002 --> -599729.89768519

> i.e., the result changes 70000 orders of magnitude right here. Ugh.
> This is affecting some erlang calculations of mine pretty badly.

Make that 30400 orders of magnitude (natural logs y'know)... On something thats about 300000 orders of magnitude below 1, mind you! What the devil are you calculating? The probability that a random configuration of atoms would make up the known universe?  

> Judging from comments in pgamma someone else might have gotten hit by this
> around R1.8.1.

Beginning to look like a method bug. Recipe is to read the original paper, check for transcription errors from the Fortran algorithm, then check for errors going from theory to Fortran, then check the theory...

PS: It happens on the R command line too:

> pgamma(1000000,100000.001,1.0,,0,1) - pgamma(1000000,100000,1.0,,0,1)
[1] 70017.54

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907

______________________________________________
R-devel@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Sat Oct 23 07:50:43 2004

This archive was generated by hypermail 2.1.8 : Fri 18 Mar 2005 - 09:00:56 EST