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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
