# [R] bug? and New bug. --- patch for pt() only (PR#138)

**maechler@stat.math.ethz.ch
**
maechler@stat.math.ethz.ch

*Wed, 10 Mar 1999 18:43:07 +0100*

The following patch saves pt(), but not the
pf() and pbeta() ones which are harder..
[pbeta(), the incomplete beta ratio is really underlying all these...
needs another asymptotic formula, and it's even harder to decide
WHEN to switch from the (Taylor kind) series to the asymptotic formula
]
This has to wait for a while, unless someone else....
BTW, S-plus also fouls up completely here:
df3 <- c(10^(3:20))
## completely wrong on S; taking forever on R
plot(df3, ptd3 <- pt(1.96,df=df3), log='x', type='b')
ptd3
giving
[1] 9.748634e-01 9.749882e-01 9.750007e-01 9.750020e-01
[5] 9.750021e-01 9.750021e-01 9.750021e-01 9.750023e-01
[9] 9.750054e-01 9.750436e-01 9.748710e-01 9.725377e-01
[13] 8.914032e-01 -4.261217e+06 NA NA
[17] NA NA
where the last before the NA's is -426121.7 !!
-----------------------
Here is the patch
---------------------------------------------------------------------
--- src/nmath/pt.c.~1~ Wed Feb 3 12:21:46 1999
+++ src/nmath/pt.c Wed Mar 10 18:15:42 1999
@@ -41,6 +41,11 @@
if(!finite(n))
return pnorm(x, 0.0, 1.0);
#endif
+ 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);
+ }
val = 0.5 * pbeta(n / (n + x * x), n / 2.0, 0.5);
return (x > 0.0) ? 1 - val : val;
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._