From terra at gnome.org  Mon Oct 25 18:04:08 2004
From: terra at gnome.org (Morten Welinder)
Date: Mon Oct 25 18:04:31 2004
Subject: [Rd] pgamma discontinuity (PR#7307)
In-Reply-To: <x2fz46h1lo.fsf@biostat.ku.dk> (message from Peter Dalgaard on 22
	Oct 2004 23:41:55 +0200)
References: <20041022184744.457E61094C@slim.kubism.ku.dk>
	<x2fz46h1lo.fsf@biostat.ku.dk>
Message-ID: <20041025160408.BDF7B1422D5E@darter.rentec.com>


A little code study, formula study and experimentation reveals that the
situation is mostly fixable:

1. Get rid of the explicit alpha limit.  (A form of it is implicit in
   (2) and (3) below.)

2. Use the series formula when

       (x < alph + 10 && x < 0.99 * (alph + 1000))

   This guarantees that the sum converges reasonably fast.  (For extremely
   large x and alpha that will take about -53/log2(0.99) iterations for
   53 significant bits, i.e., about 3700 iterations.)

3. Use the continued fraction formula when

       (alph < x && alph - 100 < 0.99 * x)

   Aka, you don't want to use the formula either near the critical point
   where alpha/x ~ 1 unless the numbers are small.

4a. Go to a library and figure out how Temme does it for alpha near x,
    both large.  In this case the 0.99 from above could probably be
    lowered a lot for faster convergence.

or

4b. Use the pnorm approximation.  It seems to do a whole lot better for
    alpha near x than it did for the 10:1 case I quoted.

Comments, please.

Morten


[Rd] pgamma discontinuity (PR#7307)
Morten Welinder