[R] remarks on function pnbeta

From: Ali Baharev <ali.baharev_at_gmail.com>
Date: Fri, 25 Apr 2008 16:14:12 +0200


Dear Developers,

The pnbeta function has been reviewed recently in the article of A. Baharev, S. Kemény, 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



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Fri 25 Apr 2008 - 14:20:25 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 25 Apr 2008 - 14:30:32 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive