Skip to content

non-centrality parameter in pf() (PR#752)

3 messages · Peter Dalgaard

Original
Peter Dalgaard
# 
jerome@stat.ubc.ca writes:
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?
Peter Dalgaard
# 
Peter Dalgaard BSA <p.dalgaard@biostat.ku.dk> writes:
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.
Peter Dalgaard
# 
Peter Dalgaard BSA <p.dalgaard@biostat.ku.dk> writes:
^^^^
Whoops. Read: "large"
I have now committed a fix to appear in R v. 1.2