R help - Oct 2005 - Problems with pf() with certain noncentral values/degrees of freedom combinations

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

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