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

[Peter Dalgaard BSA]

Peter Dalgaard BSA <p.dalgaard@biostat.ku.dk> writes:

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

Heh. This one goes right back to AS R84 (1990). The whole business of
calculating x0 like that is to avoid calculating a bunch of terms that
are practically zero when the noncentrality parameter is zero, and the
paper details what changes to the FORTRAN are needed to initialize the
recurrence relation for successive terms in mid-sequence. However it
forgets to say that the starting point for XJ also needs to be
changed.... 

A simple fix is to replace j=0 with j=x0, but there are a couple of
other things that might require a touchup - in particular, the whole
thing is based on single precision, so the cutoff induces an error
estimated as pnorm(-5) < 3e-7. Also, one should probably ensure that
the max iteration count is kept constant.

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