qgamma precision (PR#2214)

[terra@diku.dk]

Full_Name: Morten Welinder
Version: 1.5.1
OS: Solaris
Submission from: (NULL) (65.213.85.136)


I was having problems with qgamma's precision in the sense that
pgamma(qpgamma(x))
was not as close to the identity function as I would like.  I was seeing
relative
errors with random input of about 1e-8.  This fits nicely with the code'd EPS2
value of 5e-7.

To solve this I added a newton step at the end.  This seems far cheaper than
the
current iterations.  Here's the code I added (to replace the final return):


    /* Special Gnumeric patch to improve precision.  */
    {
	gnum_float x0 = 0.5*scale*ch;
	gnum_float e0 = pgamma (x0, alpha, scale, lower_tail, log_p) - p;

	if (e0 != 0 && lower_tail && !log_p) {
	    gnum_float d0 = dgamma (x0, alpha, scale, log_p);
	    if (d0) {
		gnum_float x1 = x0 - e0 / d0;
		gnum_float e1 = pgamma (x1, alpha, scale, lower_tail, log_p) - p;
		if (gnumabs (e1) < gnumabs (e0))
		    x0 = x1;
	    }
	}

	return x0;
    }

(Just read "gnum_float" as "double" and "gnumabs" as "fabs".)

The test for lower_tail is because (a) that is where I use it, and (b) the d0
calculation is probably wrong otherwise.

With this change, the relative precision improves to something like 1e-15.


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._