[R] Problems with pf() with certain noncentral values/degrees of freedom combinations

[Thomas Lumley]

The problem is in src/nmath/pnbeta.c, which has an iteration limit of 
100, not enough for these problems.  Increasing the iteration limit to 
1000 seems to work.

 	-thomas


On Mon, 24 Oct 2005, Ken Kelley wrote:

> Hello all.
>
> It seems that the pf() function when used with noncentral parameters can
> behave badly at times. I've included some examples below, but what is
> happening is that with some combinations of df and ncp parameters,
> regardless of how large the quantile gets, the same probability value is
> returned. Upon first glance noncentral values greater than 200 may seem
> large, but they are in some contexts not large at all. The problems with
> pf() can thus have serious implications (for example, in the context of
> sample size planning).
>
> I noticed that in in 1999 and 2000 issues with large degrees of freedom
> came about (PR#138), but I couldn't find the present issue reported
> anywhere.
>
> Might there be a way to make the algorithm more stable? I'm not sure how
> difficult this issue might be to fix, but hopefully it won't be too bad
> and can be easily done. Any thoughts on a workaround until then?
>
> Thanks,
> Ken Kelley
>
> # <Begin example code>
> X <- seq(10, 600, 10)
>
> # Gets stuck at .99135
> ############################
> round(pf(X, 10, 1000, 225), 5)
> round(pf(X, 10, 200, 225), 5)
>
> round(pf(X, 5, 1000, 225), 5)
> round(pf(X, 5, 200, 225), 5)
>
> round(pf(X, 1, 1000, 225), 5)
> round(pf(X, 1, 200, 225), 5)
>
> # Gets stuck at .97035
> ############################
> round(pf(X, 10, 1000, 250), 5)
> round(pf(X, 10, 200, 250), 5)
>
> round(pf(X, 5, 1000, 250), 5)
> round(pf(X, 5, 200, 250), 5)
>
> round(pf(X, 1, 1000, 250), 5)
> round(pf(X, 1, 200, 250), 5)
>
> # Gets stuck at .93539
> ############################
> round(pf(X, 10, 1000, 275), 5)
> round(pf(X, 10, 200, 275), 5)
>
> round(pf(X, 5, 1000, 275), 5)
> round(pf(X, 5, 200, 275), 5)
>
> round(pf(X, 1, 1000, 275), 5)
> round(pf(X, 1, 200, 275), 5)
> # <end example code>
>
> > version
>          _
> platform i386-pc-mingw32
> arch     i386
> os       mingw32
> system   i386, mingw32
> status
> major    2
> minor    2.0
> year     2005
> month    10
> day      06
> svn rev  35749
> language R
>
> -- 
> Ken Kelley, Ph.D.
> Inquiry Methodology Program
> Indiana University
> 201 North Rose Avenue, Room 4004
> Bloomington, Indiana 47405
> http://www.indiana.edu/~kenkel
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

Thomas Lumley			Assoc. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle