Dear Developers,
The pnbeta function has been reviewed recently in the article of A.
Baharev, S. Kem=E9ny, On the computation of the noncentral F and
noncentral beta distribution, Statistics and Computing, 2008, in
press. Preprint of the paper is available here:
http://reliablecomputing.eu/publications.html
and the article here
http://dx.doi.org/10.1007/s11222-008-9061-3
I suggest increasing the itrmax to (at least) 2000, since practically
important values cannot be computed with the current value. For
example:
pnbeta(0.99992057, 25.0, 0.5, 34012.98, 1, 0);
needs 1649 iterations on my machine, and returns 0.0752037 instead of
the true value 0.100000. This is not a bug, i got a warning, however i
think it would be better to increase the iteration limit.
According to my experience the users tend to ignore warning messages
and would use the incorrect values that is why i think it would be
better if pnbeta called exit(EXIT_FAILURE) or throw an exception than
returning a completely incorrect value and let the user use that. This
is my personal opinion, i do not know more about the design and
application logic of R.
The source file pnbeta.c refers to the paper of Frick (1990, AS R84)
giving a FORTRAN implementation. That paper contains typos; those were
corrected by Lam (1995, AS R95). Lam also made some remarks that are
worth considering.
The R implementation is a revised C code, and it seems to me that the
pnbeta's implementation is not influenced by the typos in the paper of
Frick.
However it would be better to indicate that the pnbeta is a revised
implementation compared to AS 226 and AS R84, and is not influenced by
AS R95.
Cancellation may occur during the recursion in the do while loop, we
failed to produce an example using practically meaningful parameter
values where serious cancellation occurs. We suspect it may be an
issue where the corresponding probability is small (less then
1.0E-30).
We plan to contribute a C implementation of our algorithm where
cancellation cannot occur, and newton iteration for the noncentrality
parameter is included. The computation of the noncentrality parameter
is need for ANOVA purposes.
Best regards,
Ali Baharev
