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