Skip to content
Prev 4690 / 63424 Next
Original Raw

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

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?

Thread (3 messages)

Peter Dalgaard non-centrality parameter in pf() (PR#752) Dec 5 Peter Dalgaard non-centrality parameter in pf() (PR#752) Dec 5 Peter Dalgaard non-centrality parameter in pf() (PR#752) Dec 5