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

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