10.

ll.

12.

13.

14.

References

Abramowitz, Milton and Stegun, Irene A., Handbook of Mathematical
Functions, NBS Special Publication 55, 1970, p. 944.

Bulgren, W.G., "On Representations of the Doubly Non-Central F
Distribution", Journal of the American Statistical Association, Vol. 66,
No. 333, March 1971, pp. 184-6.

IMSL, Inc., Houston, TX. [MDBETA, DLGAMA]

Johnson, Norman L. and Kotz, Samuel, Continuous Univariate
Distributions-2, Houghton Mifflin Company, 1970, pp. 189-200.

Mudholkar, Govind S., Chaubey, Yogendra P., and Lin, Ching-Chuong,
“Approximations for the Doubly Noncentral-F Distribution", Communications
in Statistics - Theory and Methods, Vol. A5, No. 1, 1976, pp-49-63.

NBS Core Math Library (CMLIB). [BETAI, DLNGAM]
Numerical Algorithms Group (NAG), Downers.Grove, IL. [GOIBDE, S14ABF]

Price, Robert, "Some Noncentral F Distributions Expressed in Closed
Form", Biometrika, Vol. 51, 1964, pp. 107-22.

Reeve, Charles P., "An Algorithm for Computing the Beta C.D.F. to a
Specified Accuracy", SED Note 86-3, October 1986.

Reeve, Charles P., "Accurate Computation of the Log of the Gamma
Function", SED Note 86-1, October 1986.

Scheffe, Henry, The Analysis of Variance, John Wiley & Sons, Inc., 1959,
pp. 134-6.

Tiku, M.L., "Series Expansions of the doubly Noncentral F Distributions",
Australian Journal of Statistics, Vol. 7, 1965, pp. 78-89.

Tiku, M.L., "A Note on Approximating to the Doubly Noncentral F
Distribution", Australian Journal of Statistics, Vol. 14, 1972, pp. 37-40.

Tiku, M.L., "Doubly Noncentral F Distribution - Tables and Applications",
Selected Tables in Mathematical Statistics, Vol. II, 1974, pp. 139-49.



SUBROUTINE POISSF (ALAMB, EPS, L,NSPAN,V,NV, IFLAG)

C—— COMPUTE THE POISSON(ALAMB) PROBABILITIES OVER THE RANGE [L,K]
C—— WHERE THE TOTAL TAIL PROBABILITY IS LESS THAN EPS/2, SUM THE
C—— PROBABILITIES IN DOUBLE PRECISION, AND SHIFT THEM TO THE
C—— BEGINNING OF VECTOR V.

DIMENSION V(*)
DOUBLE PRECISION DAL, DK ,DLIMIT ,DSUM,GAMLOG
OLIMIT = 1 .@D@-@.5D@*DBLE(EPS)
K = INT(ALAMB)
L = K+1
IF (ALAMB.£Q.@.@) THEN
PL = 1.0
Eu

SE
DAL = DBLE(ALAMB)
OK = DBLE(K)
PL = SNGL(DEXP(DK*DLOG(DAL)-DAL-GAMLOG (REAL(K+1))))
ENDIF
PK = ALAMB*PL/REAL(L)
NK = NV/2
NL = NK+T
DSUM = 0.0
40 IF (PL.LT.PK) THEN
NK = NK+1
1F (NK.GT.NV) THEN
IFLAG = 6
RETURN
ENDIF
V(NK) = PK
DSUM = DSUM+DBLE(PK)
K = Ket
IF (DSUM.GE.DLIMIT) GO TO 2@
PK = ALAMB*PK/REAL(K+1)
ELSE
NL = NL-1
V(NL) = PL
DSUM = DSUM+DBLE(PL)
t= t-1
IF (DSUM.GE.DLIMIT) GO TO 20
PL = REAL(L)*PL/ALAMB
ENDIF
GO TO 10
2@ INC = NL-1
DO 30 I = NL, NK
V(I-INC) = V(I)
3@ CONTINUE
NSPAN = NK-INC
RETURN
END

SUBROUTINE EDGEF (NK, FC,GC,XX,YY,BFK,CDFX, POI »POJ,EPS3, IFLAG,L)

C—— COMPUTE THE BETA C.D.F.°S BY A RECURRENCE RELATION ALONG THE EDGES
C—— I = IMIN AND J = JMIN OF A GRID. THE CORRESPONDING COMPONENTS OF
C—— THE F" C.D.F. ARE INCLUDED IN THE SUMMATION. TERMS WHICH MIGHT
C—— CAUSE UNDERFLOW ARE SET TO ZERO.

DIMENSION BFK(#) ,POI(*) ,POU(*)
DOUBLE PRECISION DARG,DEUFLO,GAMLOG
DATA DEUFLO / -69.208 /

FD = FC-1.0
K = MAX@(L,MINO(NK, INT((GC~1.0)*XX/YY-FD) ))
FK = FD+REAL(K)

CALL CDFBET (XX,FK,GC,EPS3, 1FLAG,BFK(K))
IF (1FLAG.NE.@) RETURN
IF (L.EQ.1) BFK(K) = 1.@-BFK(K)
IF (NK.E£Q.1) GO TO 40
DARG = DBLE( FK) +DLOG(DBLE (XX) )+DBLE(GC)*DL0G(DBLE( 17} )-
» DLOG(DBLE(FK) )-+GAMLOG ( FK+GC)-GAMLOG(FK)—-GAMLOG (GC
IF (DARG.LT.DEUFLO) THEN
DK = 0.0

LSE
DK = SNGL(DEXP(DARG))*(-1.0)=#L
ENDIF
IF (K.GE.NK) GO TO 20
BFK(K+1) = BFK(K)-DK
DI = DK
KFLAG = 1
DO 10 1 = K+1, NK-1


where the arguments are defined in the program documentation. The returned
value of CDFX is valid only if IFLAG=0 on return. In passing € (variable name
EPS) to CDFDNF the user should realize that accuracy is limited by the number
of digits carried in a single precision variable, and that roundoff error may
affect the last one or two of these digits.

Table 1

Computing times on the CDC Cyber 180/855 for the
c.def. of EN(V)5Va9A AQ) using both CDFDNF and

the "brute force" method (described in the text),
for selected parameter values.

e=10°° CPU seconds No. beta
———TRrube

YO AY de x P{F"<x}  CDEDNF nae orate
3 3 5 5 2.0 0.757918 ~<0.01 0.13 196
3 3 5 25 2.0 0.997561 <0-01 0.30 476
3 3 25 5 2.0 0.190910 <0.01 0.31 476
3 3 25 25 2.0 0.897835 <0.01 0.69 1,156
3 10 5 5 2.0 0.593795 <0-01 0.14 196
3 10 5 25 2.0 0.943093 <0.01 0.30 476
3 10 25 5 2.0 0.026209 <0.01 0.28 476
3 10 25 25 2.0 0.289601 <0.01 0.73 1,156
10 3 5 5 2.0 0.898330 <0.01 0.12 196
10 3 5 25 2.0 0.999879 <0.01 0.26 476
10 3 25 5 2.0 0.657879 <0-01 0.29 476
10 3 25 25 2.0 0.997703 <0.01 0.66 1,156
10 10 5 5 2.0 0.868071 <0.01 0.14 196
10 10 5 25 2.0 0.998234 <0-01 0.30 476
10 1o 25 5 2.0 0.367101 <0.01 0.32 476
10 10 25 25 2.0 0.934321 <0.01 0.71 1,156
14 15 80 80 lel 0.552328 0.03 2.55 3,969
14 15 400 400 lel 0.582507 0.12 16.28 20,164
14 15 2000 2000 lel 0.664981 0-59 108.10 101,124
14 15 10000 10000 ll 0.825080 2.68 =800 505,521
cases I5——30p00 50000 T.1 0.981351 13.76 - 2,528,100

ft required a doubling of current dimension limits


*CDFDNF

SUBROUTINE COFONF (X,OF1,0F2,ALAMB1,ALAMB2, EPS, IFLAG,CDFX)

CDFDNF WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING

DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
MARYLAND 20899

FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE DOUBLY

NONCENTRAL F DISTRIBUTION TO A SPECIFIED ACCURACY (TRUNCATION
ERROR IN THE INFINITE SERIES REPRESENTATION GIVEN BY EQUATION
2.2 IN REFERENCE 1 BELOW). THE BETA C.D.F. ROUTINE IS CALLED
AT MOST TWO TIMES. FURTHER VALUES OF THE BETA C.D.F. ARE
OBTAINED FROM RECURRENCE RELATIONS GIVEN IN REFERENCE 2.
REFERENCE 3 GIVES A DETAILED DESCRIPTION OF THE ALGORITHM
HEREIN.

THIS PROGRAM MAY ALSO BE EFFICIENTLY USED TO COMPUTE THE
CUMULATIVE DISTRIBUTION FUNCTIONS OF THE SINGLY NONCENTRAL
AND CENTRAL F DISTRIBUTIONS BY SETTING THE APPROPRIATE
NONCENTRALITY PARAMETERS EQUAL TO ZERO.

CHECKS ARE MADE TO ASSURE THAT ALL PASSED PARAMETERS ARE
WITHIN VALID RANGES AS GIVEN BELOW. NO UPPER LIMIT IS SET
FOR THE NONCENTRALITY PARAMETERS, BUT VALUES UP TO ABOUT
1@,@@0 CAN BE HANDLED WITH THE CURRENT DIMENSION LIMITS. THE
COMPUTED VALUE CDFX IS VALID ONLY IF IFLAG=@ ON RETURN.

NOTE: IN EQUATION 2.2 OF REFERENCE 1 THE AUTHOR HAS MISTAKENLY

REVERSED THE ARGUMENTS OF THE INCOMPLETE BETA FUNCTION.
THEY SHOULD READ [(M/2)+R,(N/2+S)] WHERE M AND N ARE THE
DEGREES OF FREEDOM ASSOCIATED WITH THE NUMERATOR AND
DENOMINATOR RESPECTIVELY OF THE F STATISTIC. TO FURTHER
CONFUSE THE ISSUE, THE AUTHOR HAS REVERSED THE USAGE OF
M AND N IN SECTION 1 OF THE PAPER.

NOTE: IN SUBROUTINE EDGEF THE DOUBLE PRECISION CONSTANT DEUFLO IS

THE EXPONENTIAL UNDERFLOW LIMIT WHOSE CURRENT VALUE IS SET
AT -69D2. ON A COMPUTER WHERE DEXP(—69D®) CAUSES UNDERFLOW
THIS LIMIT SHOULD BE CHANGED.

SUBPROGRAMS CALLED: CDFBET (Rove C.0.F.)

GAMLOG (DOUBLE PRECISION LOG OF GAMMA FUNCTION)
POISSF, EDGEF (ATTACHED)

CURRENT VERSION COMPLETED SEPTEMBER 29, 1988

REFERENCES:

1.

BULGREN, W.G., ’ON REPRESENTATIONS OF THE DOUBLY NONCENTRAL F
DISTRIBUTION’, JOURNAL OF THE AMERICAN STATISTICAL ASSOCIATION,
MARCH 1971, VOLUME 66, NO. 333, PP. 184-186.

ABRAMOWITZ, MILTON, AND STEGUN, IRENE A., "HANDBOOK OF
MATHEMATICAL FUNCTIONS’, NATIONAL BUREAU OF STANDARDS APPLIED
MATHEMATICS SERIES 55, NOVEMBER 1970, P. 944.

. REEVE, CHARLES P., ‘AN ALGORITHM FOR COMPUTING THE DOUBLY

NONCENTRAL F C.D.F. TO A SPECIFIED ACCURACY’, STATISTICAL
ENGINEERING DIVISION NOTE 86-4, NOVEMBER 1986.

00000000000000000009000000000000000000000099900000009999099999090 990900 999990090

DEFINITION OF PASSED PARAMETERS:

* X = VALUE (>=@) AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL)
* DF1 = DEGREES OF FREEDOM (>0) IN THE NUMERATOR (REAL)
* DF2 = DEGREES OF FREEDOM (>@) IN THE DENOMINATOR (REAL)

* ALAMB1 = THE NONCENTRALITY PARAMETER (>=) FOR THE NUMERATOR

(REAL) [EQUAL TO ZERO FOR THE CENTRAL F DISTRIBUTION]

* ALAMB2 = THE NONCENTRALITY PARAMETER (>=0) FOR THE DENOMINATOR

(REAL) [EQUAL TO ZERO FOR THE SINGLY NONCENTRAL F AND
CENTRAL F DISTRIBUTIONS]

* EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (REAL)
[1 >= EPS >= 10#«(-10)]

IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION:


30
40

50

IF (KFLAG.EQ.1) THEN
DI = DI*(FD+GC+REAL(I-1))#XX/(FD+REAL(I))
IF (DK+DI.EQ.0K) THEN
KFLAG = @
DI = 0.0
ENDIF
ENDIF
BFK(I+1) = BFK(I)-DI
CONTINUE
DI = DK
KFLAG =
DO 30.1 =K-1, L, -1
IF (KFLAG.EQ.1) THEN
DI = DI*(FC+REAL(I))/((FD+GC+REAL(I))*XX)
IF (DK+DI.£0.0K) THEN
KFLAG = @
DI = 0.0
ENDIF
ENDIF
BFK(I) = BFK(I+1)4DI
CONTINUE
DO 50 I =L, NK
CDFX = CDFX+POI(I)*POU(1)*BFK(I)
CONTINUE
RETURN
END



National Bureau of Standards N-86-4
Statistical Engineering Division November 1986
Gaithersburg, MD 20899

An Algorithm for Computing the Doubly Noncentral F
C.D.F. to a Specified Accuracy

Charles P. Reeve

Let x, and x be independent noncentral chi-squared random variables with
degrees of freedom v) and vy (both >0) and noncentrality parameters at and AD
(both >0) respectively. Then the random variable

¥ = (X,/¥)/(K,/¥) qQ)

is said to have the doubly noncentral F distribution, indicated by
Y~ EMV sVardyrAo) > This distribution has been used in the evaluation of the

power function of analysis of variance tests in which interaction or bias

effects occur, as in Scheffe~ [11]. Numerical examples of this-usage—are
iven i i [13]. It has also been used in engineering

problems in the context of information theory as discussed in Price [8].

Exact formulas for the F" cumulative distribution function (c.d.f.) are given
in Tiku [12] and [2] using the beta c.d.f., in Tiku [14] using Laguerre
Polynomials, and in [8] for special cases of yy and Vor Approximations are

iven son and Kotz [4] and Tiku [12,13] using the central F ced.f,, and

in Mudholkar, Chaubey, and Lin [5] using Edgeworth series expansions.

The author has been unable to find any published algorithms for computing

exact values of the F" c.d.f. although computer programs have obviously been
used in generating published tables in [2,5,12,14]., The purpose of this note

is to present an_efficient algorithm for. _bomputing the F" c.d.f. toa
specified accuracy using exact formulas.

The algorithm uses the series representation in eq. (2.2) of [2] which,
unfortunately, contains typographical errors. With these corrected the c.d.f.
of Y can be re-written

a
F(x) = ~ 14 B, I(u,v, /2+4,V,/245) (2)
i=0 j=0
ays -v,/2

. 2 .
where A, = (A,/2)7e IPGL), Be = O,/2y4e 2 rth, us vx/(v x45)

and x20. The Ay and Bs are Poisson probabilities, and

u
1 — »
I(u,a,b) =f t?7*1-t)>laey(a,b) is the c.d.f. of the beta distribution
0
(also called the incomplete beta ratio) where O<u<l, a>0, b>0, and


—> NO ERRORS DETECTED

—> ERROR FLAGS FROM SUBROUTINE CDFBET

—> EITHER ALAMB1 OR ALAMB2 IS < @

—> EITHER DF1 OR DF2 IS <= @

—> EPS IS OUTSIDE THE RANGE [10**(-10),1]

—> VECTOR DIMENSIONS ARE TOO SMALL — INCREASE NX

OUhUNS

CDFX = THE DOUBLY NONCENTRAL F C.D.F. EVALUATED AT X (REAL)
*» INDICATES PARAMETERS REQUIRING INPUT VALUES

2900009000000

PARAMETER (NX=1000)
DIMENSION BFI(NX),BFU(NX) ,POI (NX) ,POJ(NX)
CDFX = 0.0

c
C—— CHECK VALIDITY OF ARGUMENTS

IF (ALAMB1.LT.@.@.O0R.ALAMB2.LT.@.0) THEN
IFLAG = 3
RETURN

ENDIF

IF (DF1.LE.@.0.OR.DF2.LE.0.0) THEN
IFLAG = 4
RETURN

ENDIF

IF (EPS.GT.1.0.OR.EPS.LT.1.@E-10) THEN
IFLAG = 5
RETURN

ENDIF

IFLAG = 0

Cc
C—— SET ERROR CRITERION FOR THE BETA C.D.F. (PECULIAR TO CDFBET)

c
Cc

EPS3 = @.001*EPS

FA = 0.5*ALAMB1

GA = @.5*ALAMB2

@.5*DF1

@.5*DF2

DF2/(DF2+DF 1x)

IF (YY.GE.1.0) RETURN

XX = 1.0-YY

IF (XX.GE.1.@) THEN
COFX = 1.0
RETURN

ENDIF

°
oO
wow

COMPUTE POISSON PROBABILITIES IN VECTORS POI AND POU

CALL POISSF (FA,EPS,IMIN,NI,POI,NX,IFLAG)
IF (IFLAG.NE.@) RETURN

FC = FB+REAL(IMIN)

CALL POISSF (GA,EPS,JMIN,NJ,POU,NX, IFLAG)
IF (IFLAG.NE.@) RETURN

GC = GB+REAL(JMIN)

COMPUTE BETA C.D.F. BY RECURRENCE WHEN I=IMIN AND J=JMIN TD <JMAX

CALL EDGEF (NJ,GC,FC,YY,XX,BFU,CDFX,POJ,POI,EPS3, IFLAG,1)
IF (NI.LE.1.OR.IFLAG.NE.@) RETURN

COMPUTE BETA C.D.F. BY RECURRENCE WHEN J=JMIN AND I=]MEN ‘TOSMAX
BFI(1) = BFU(1)

CALL EDGEF (NI,FC,GC,XX,YY,BFI,CDFX,POI ,POJ,EPS3, IFEAG,2)
IF (NJ.LE.1.0R.IFLAG.NE.@) RETURN

C—— COMPUTE BETA C.D.F. BY RECURRENCE WHEN I>IMIN AND J>JMIN-

DO 20 I = 2, NI
BFJ(1) = BFI(I)
DO 10.3 = 2, NU
BFU(J) = XX*BFU(J)+YY*BFU(J-1)
CDFX = CDFX+POI(1)*POJ(J)*BFU(J)

10 CONT INUE
2@ CONTINUE

RETURN
END


1

B(a,b) = f ce tty hae. For computational purposes the two infinite
0

series must be truncated, thus (2) is re-expressed as

mo
Fy(x) = oy Pe A,B, T(u,v, /2+i,v,/2+j) +R (3)

where I', I", J', and J" are non-negative integers and Ris the remainder. If
the beta c.d.f.'s are computed without error, it can easily be shown that
q" x"
choosing I', I", J', and J" such that J A, > l-e/2 and l By > l-c/2 yields
i=! jeu"
R<e provided €>0. Therefore, € serves as an absolute error bound on Fy (x).
For maximum computational efficiency, the number of terms in each sum is
minimized by indexing i and j over the largest of the Poisson probabilities
A, and By respectively. The final task is to compute the (I"-I'+1)(J"-J'+l1)

beta c.d.f.'s and the double summation. An efficient procedure for doing this
is to compute only I(u,v, /2+1' Vv, /2+I*) and (u,v /2+1*,v,/2+J") directly,

indicated by the symbols "x" and "y" in figure 1. The remaining beta c.d.f's
are computed using the recurrence relations

I(x,ayb) = Ux,a,bel) - x*(1-x)P/[bB(a,b) 1, (4a)
I(x,a,b) = I(x,atl,b) + x?(1-x)?/[aB(a,b)], and (4b)
I(x,a,b) = xI(x,a-l,b) + (1-x)I(x,a,b-1) (4c)

as found in Abramowitz and Stegun [1]. Subject to the restrictions J'<J*<J"

jas" -| aoc ec coe eo
- a c c c c c c
-| ac ec ec ec coc
j=I* -| x c oe. © © C we. Ce
-| ac ec ec coe
-| ac ec coe eoe
j=s' -| ab boy b bob
TT rT rT

isi! isi* isi"

Figure 1

Computation of the I(u,v, /2+i,v,/2+j)


and I'<I*<I", J* and I* are chosen to maximize the magnitudes of the rightmost
terms in (4a) and (4b) respectively. In applying each of these two recurrence
relations only one direct evaluation of B(a,b) is necessary, computed by

Inf(a) + Inf(b) - Inf(a+tb) al -t

o
B(a,b) =e where T(a) = f t° "e “dt . Further
values are easily computed from the identities B(a+l,b) = aB(a,b)/(atb) and
B(a,b+l) = bB(a,b)/(atb).

In figure 1 the symbols "a", "b", and "c" indicate which of the recurrence
relations (4) is used in computing each beta c.d.f., with those indicated by
“e" being done last.

The above algorithm has been incorporated into the FORTRAN subroutine CDFDNF.
If x<0 then 0 is returned. External routines for computing the beta c.d.f.
and the double precision log of the gamma function are required. In the
current version of CDFNDF these routines are the subroutine CDFBET and the
function GAMLOG as described in Reeve [9,10]. Other routines which can be
substituted for CDFBET and GAMLOG are those in [3,6,7]. Farallpractical
purposes_the absolute error criterion e will_be met if the beta ced.f. routine
i ur. o two or three digits beyond c€.

The recursive method of computing the beta c.d.f.'s requires a little extra
computer programming and storage, but results in a tremendous savings in
computing time as y and AD become large. ‘his is illustrated in table 1

where, for selected parameter values, the CPU time for this method is compared
with that of the "brute force" approach in which each beta c.d.f. is directly
evaluated by a subroutine call. The computing time and storage also increase

as € becomes small, but are unaffected by yy and Vo: The computations in

table 1 were done on the CDC Cyber 180/855 computer at NBS.

The current dimensidén limits in CDFDNF allow values of both A, and i, up to

about 10,000 with € as small as 10°, but these limits could easily be
increased by the user. The limiting factor in using CDFDNF is more likely to
be execution time than storage.

Steps were taken to eliminate underflow situations, minimize the effects of

roundoff error, and minimize storage requirements. Only those c.d.f. values

indicated by "x", "y", "a", or "b' in figure 1 are actually stored. The

Poisson probabilities Ales wees Ain and Bye ones Ban are also stored.

If 44-0 then the doubly noncentral F reduces to the (singly) noncentral F, and
if A =A ,=0 it reduces to the central F. In either case, CDFDNF will run

almost as efficiently as routines designed for those specific cases.

Portions 9, les in [2,5,1 roduced and agreed to within

in each case.
A listing of CDFDNF is an appendix to this note. It is invoked by

CALL CDFDNF(X, DF1,DF2,ALAMB1,ALAMB2, EPS, IFLAG, CDFX)