[R] Problems with pf() with certain noncentral values/degrees of freedom combinations
Thomas Lumley
tlumley at u.washington.edu
Thu Nov 3 20:37:25 CET 2005
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
More information about the R-help
mailing list