R-alpha: New pbeta.c : Still problems for BIG 'df' (and an S-bug)

[Martin Maechler]

I've only made a few experiments using  qt(.)  with the new pbeta.c:

At least the problem of  "BIG df" --> 'long/infinite' CPU  remains even
though it has been pushed slightly

For 0.49 'plain', 
	qt(.99, 10^k) already takes quite some time for k=5, and 
			takes for ever for k=7
Try
	for(ee in 1:6) cat(10^ee,":", unix.time(qt(.99,df=10^ee))[1],"\n")

For the new  pbeta.c,  the problem appears about two exponents later.

 	qt(.99, df = 1e8)	## 'never' returns

The following code shows you a bit of what is going on.

BTW, S-plus seems to use a different approximation which remains fast even
when df = BIG.
However, S-plus results then become COMPLETELY RUBBISH!

Here are my few tests:

p <- c(.975, .995, .9995, 1 - 5e-6)
df1 <- c(1:5,10^c(1:8,10,15,10*(2:10)))
df0 <- df1[1:12] ## df0  works for  R 0.49 -patch2-

dig <- 8
for(df in df0)
  cat(formatC(df,w=6),":", formatC(unix.time(qq <- qt(p, df))[1],form='f'), "|",
      formatC(qq, form='f', dig=dig,wid=dig+6),"\n")
cat(" _Inf_ :  ---   |", formatC(qnorm(p), form='f', dig=dig,wid=dig+6),"\n")

## R 0.49  with  pbeta.c  of 30 April 1997:

   df  :  CPU[s]|     q  .975           .995         .9995       .999995=1-5e-6

     1 : 0.0100 |    12.70620474    63.65674116   636.61924876 63661.96977709
     2 : 0.0200 |     4.30265273     9.92484320    31.59905458   316.22539430
     3 : 0.0000 |     3.18244631     5.84090931    12.92397864    60.39682404
     4 : 0.0000 |     2.77644511     4.60409487     8.61030158    27.77164766
     5 : 0.0100 |     2.57058184     4.03214298     6.86882663    17.89686614
    10 : 0.0100 |     2.22813885     3.16927267     4.58689386     8.15028656
   100 : 0.0900 |     1.98397152     2.62589052     3.39049131     4.65424006
  1000 : 0.0300 |     1.96233908     2.58075470     3.30028265     4.43992647
 1e+04 : 0.0400 |     1.96020124     2.57632105     3.29149997     4.41943950
 1e+05 : 0.1300 |     1.95998771     2.57587847     3.29062403     4.41739993
 1e+06 : 0.2800 |     1.95996636     2.57583422     3.29053646     4.41719606
 1e+07 : 0.3800 |     1.95996422     2.57582979     3.29052770     4.41717568
 _Inf_ :  ---   |     1.95996398     2.57582930     3.29052672     4.41718074

## S-plus (3.4): From df ~= 1e4 (p=1-5e-6) results are COMPLETELY wrong !
## ----- You get quick results for even higher 'df' -- but they are wrong...

     1 : 0.0100 |    12.70620474    63.65674116   636.61924876 63661.96770658
     2 : 0.0000 |     4.30265273     9.92484320    31.59905458   316.22539430
     3 : 0.0000 |     3.18244631     5.84090931    12.92397864    60.39682404
     4 : 0.0100 |     2.77644511     4.60409487     8.61030158    27.77164766
     5 : 0.0000 |     2.57058184     4.03214298     6.86882663    17.89686614
    10 : 0.0000 |     2.22813885     3.16927267     4.58689386     8.15028656
   100 : 0.0100 |     1.98397152     2.62589052     3.39049131     4.65424006
  1000 : 0.0200 |     1.96233908     2.58075470     3.30028265     4.43992647
 1e+04 : 0.0200 |     1.96020124     2.57632105     3.29966223     4.70204286
 1e+05 : 0.0200 |     1.95998771     2.57587848     5.14482712    12.48628360
 1e+06 : 0.0100 |     1.95996636     2.57583569    14.70818229    38.67299101
 1e+07 : 0.0200 |     1.95996422     2.57600108    45.98857869   122.03496539
 _Inf_ :  ---   |     1.95996398     2.57582930     3.29052673     4.41717341

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
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
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-