[Rd] [Fwd: Re: pt inaccurate when x is close to 0 (PR#9945)]

From: Duncan Murdoch <murdoch_at_stats.uwo.ca>
Date: Thu, 11 Oct 2007 11:10:49 -0400


Here's a contribution from Ian Smith that got bounced from the list.

Duncan,

I tried sending the rest of this to R-devel but it was rejected as spam, hence the personal e-mail.

R calculates the pt value from

nx = 1 + (x/n)*x;
val = pbeta(1./nx, n / 2., 0.5, /*lower_tail*/1, log_p);

whereas Gnumeric calculates the value as

val =? (n > x * x)
? pbeta (x * x / (n + x * x), 0.5, n / 2, /*lower_tail*/0, log_p) : pbeta (n / (n + x * x), n / 2.0, 0.5, /*lower_tail*/1, log_p);

thus avoiding the loss of accuracy in the pbeta routine when 1-1./nx is calculated.

It also makes the

if (n > 4e5) { /*-- Fixme(?): test should depend on `n' AND `x' ! */
??? /* Approx. from? Abramowitz & Stegun 26.7.8 (p.949) */
??? val = 1./(4.*n);
??? return pnorm(x*(1. - val)/sqrt(1. + x*x*2.*val), 0.0, 1.0,
???????? lower_tail, log_p);
}

code unneccessary.

Ian Smith

Personally, I think the code should also guard against the possible overflow of the x * x expressions.



Get a FREE AOL Email account with unlimited storage. Plus, share and store photos and experience exclusively recorded live music Sessions from your favourite artists. Find out more at http://info.aol.co.uk/joinnow/?ncid=548.

R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu 11 Oct 2007 - 15:10:59 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 25 Oct 2007 - 11:37:10 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.