The following table shows the peak errors (in units of epsilon)
found on various platforms with various floating point types. No comparison
to the R-2.5.1 Math library,
or to the FORTRAN implementations of AS226 or AS310 are given since these
appear to only guarantee absolute error: this would causes our test harness
to assign an "infinite" error to these libraries
for some of our test values when measuring relative error.
Unless otherwise specified any floating-point type that is narrower than
the one shown will have effectively
zero error.

Table 2.4. Errors In CDF of the Noncentral Beta

Significand Size

Platform and Compiler

α, β,λ < 200

α,β,λ > 200

53

Win32, Visual C++ 8

Peak=620 Mean=22

Peak=8670 Mean=1040

64

RedHat Linux IA32, gcc-4.1.1

Peak=825 Mean=50

Peak=2.5x104 Mean=4000

64

Redhat Linux IA64, gcc-3.4.4

Peak=825 Mean=30

Peak=1.7x104 Mean=2500

113

HPUX IA64, aCC A.06.06

Peak=420 Mean=50

Peak=9200 Mean=1200

Error rates for the PDF, the complement of the CDF and for the quantile
functions are broadly similar.

There are two sets of test data used to verify this implementation: firstly
we can compare with a few sample values generated by the R
library. Secondly, we have tables of test data, computed with this
implementation and using interval arithmetic - this data should be accurate
to at least 50 decimal digits - and is the used for our accuracy tests.

First we determine which of the two values (the CDF or its complement)
is likely to be the smaller, the crossover point is taken to be the mean
of the distribution: for this we use the approximation due to: R. Chattamvelli
and R. Shanmugam, "Algorithm AS 310: Computing the Non-Central Beta
Distribution Function", Applied Statistics, Vol. 46, No. 1. (1997),
pp. 146-156.

Then either the CDF or its complement is computed using the relations:

The summation is performed by starting at i = λ/2, and then recursing in
both directions, using the usual recurrence relations for the Poisson PDF
and incomplete beta functions. This is the "Method 2" described
by:

Of these, the Posten reference provides the most complete overview, and
includes the modification starting iteration at λ/2.

The main difference between this implementation and the above references
is the direct computation of the complement when most efficient to do so,
and the accumulation of the sum to -1 rather than subtracting the result
from 1 at the end: this can substantially reduce the number of iterations
required when the result is near 1.

The PDF is computed using the methodology of Benton and Krishnamoorthy
and the relation:

Quantiles are computed using a specially modified version of bracket_and_solve_root,
starting the search for the root at the mean of the distribution. (A Cornish-Fisher
type expansion was also tried, but while this gets quite close to the root
in many cases, when it is wrong it tends to introduce quite pathological
behaviour: more investigation in this area is probably warranted).