[Rd] non-centrality parameter in pf() (PR#752)

[Peter Dalgaard BSA]

jerome@stat.ubc.ca writes:

> Bug Description:
> 
> Problem with the function pf() when the non-centrality
> parameter is large. Here is a sample command. You should
> see a smooth line from 0 to about 55, and then the values
> of pf() go crazy from 55 to 100.
> ############################
> ncp <- seq(0,100,length=200)
> plot(ncp,pf(5,7,2,ncp=ncp))
> ############################

Confirmed. I see it on Solaris as well. This looks like trouble for
power calculations...

I wouldn't be surprised if it were related to this FIXME inside
pbeta_raw (src/nmath/pbeta.c):

     else {
        /*___ FIXME ___:  This takes forever (or ends wrongly)
          when (one or) both p & q  are huge
        */

Another "interesting" display is 

plot(ncp,pf(2,7,3,ncp=ncp))

Apparently, this relates mainly to cases with small denominator DF,
but even 

plot(ncp,log(pf(5,7,200,ncp=ncp)))

looks odd. Judging from the plots, it looks as if a loop is terminated
prematurely (discontinuous jumps suggesting that there's a change to
using a different number of iterations).

On the other hand, ncp= ca.55 seems to "mark the spot" for quite
widely different DFs and the spot where that number is important could
be in pnbeta.c, at

    x0 = floor(fmax2(c - ualpha * sqrt(c), 0.));

which becomes nonzero when ncp=54. This might indicate that the fault
is in this section of code:

    x0 = floor(fmax2(c - ualpha * sqrt(c), 0.));
    a0 = a + x0;
    lbeta = lgammafn(a0) + lgammafn(b) - lgammafn(a0 + b);
    temp = pbeta_raw(x, a0, b, /* lower = */TRUE);
    gx = exp(a0 * log(x) + b * log(1. - x) - lbeta - log(a0));
    if (a0 > a)
        q = exp(-c + x0 * log(c) - lgammafn(x0 + 1.));
    else
        q = exp(-c);

(This code looks a bit whacked: Why is the test a0 > a and not x0 > 0,
which would seem to be the same thing?)

Any numerical analysts on the line? 

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._