Table 1

Computing times on the CDC Cyber 180/855 for the

c.d.f. of t'"(v,6,A) using CDFDNT for selected

parameter values.

€=10 No. beta
+ CPU c.d.f.
v 6 a x P{t"<x} sec values
1 1 1 0.7071 0.433771 0.01 128
1 1 100 0.0995 0.498015 0.02 1,120.
1 1 10000 0.0100 0.500000 0.10 11,248
1 10 I 7.0711 0.349271 0.02 1,128
1 10 100 0.9950 0.485863 0.08 9,870
1 10 10000 0.1000 0.500000 0.67 99,123
1 100 1 70.7107 0.347264 0.11 11,256
1 100 100 9.9504 0.480221 0.66 98,490
1 100 10000 1.0000 0.500000 6.47 989,121
10 1 1 0.9535 0.490326 0.01 128
10 1 100 0.3015 0.498251 0.01 1,120
10 1 10000 0.0316 0.499892 0.10 11,248
10 10 1 9.5346 0.448390 0.02 1,128
10 10 100 3.0151 0.487089 0.08 9,870
10 10 10000 0.3161 0.500181 0.75 99,123
10 100 1 95.3463 0.441153 0.12 11,256
10 100 100 30.1511 0.480930 0.65 98,490
10 100 10000 3.1607 0.498611 6.48 989,121
100 1 1 0.9950 0.498990 0.01 128
100 1 100 0.7071 0.499248 0.01 1,120
100 1 10000 0.0995 0.499965 0.12 11,248
100 10 1 9.9504 0.490966 0.02 1,128
100 10 100 7.0711 0.493307 0.08 9,870
100 10 10000 0.9950 0.499656 0.73 99,123
100 100 1 99.5037 0.481469 0.13 11,256
100 100 100 70-7107 0.485762 0.76 98,490
100 100 10000 9.9504 0.498682 6.71 989,121

t x = 6//1 + A/v rounded to four decimal places

+ example in Carey [4] with large values of v, 6, and d


V(I-INC) = V(I)
3@ CONTINUE
NSPAN = NK-INC

RETURN

END
Cc

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

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

DIMENSION BFK(+) ,POI(*) ,POJ(*)
DOUBLE PRECISION DARG, DEUFLO,GAMLOG
DATA DEUFLO / ~69.0D0 /
FD = FC-1.0
K = MAX@(L,MIN@(NK, INT((GC-1.)*XX/YY-FD) ))
FK = FD+REAL(K)
CALL CDFBET (XX, FK,GC,EPS3, IFLAG,BFK(K))
IF (IFLAG.NE.@) RETURN
IF (L.EQ.1) BFK(K) = 1.0-BFK(K)
IF _(NK.EQ.1) GO TO 40
DARG = DBLE(FK) *DLOG(DBLE(XX) )+DBLE(GC) «DLOG(DBLE(YY))—-
+ — DLOG(DBLE(FK) )+GAMLOG (FK+GC)-GAMLOG(FK)—GAMLOG(GC
IF (DARG.LT.DEUFLO) THEN
OK = 0.0

SE
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 1@ I = K+1, NK-1
IF (KFLAG.EQ.1) THEN
DI = DI+(FD+GC+REAL(I-1)) #XX/(FD+REAL(I))
IF (DK+DI.EQ.DK) THEN
KFLAG = @
DI = 0.0
ENDIF
ENDIF
BFK(I+1) = BFK(I)-DI
1@ CONTINUE
2@ DI = DK
KFLAG = 1
DO 3@ 1 = K-1, L, -1
IF (KFLAG.EQ.1) THEN
DI = DIs(FC+REAL(1))/( (FD+GC+REAL(I))#XX)
IF (DK+DI.£Q.DK) THEN
KFLAG = @
DI = 0.0
ENDIF

ENDIF
BFK(I) = BFK(I+1)+DI
3@ CONTINUE
40 DO 50 I = L, NK
CDFX = CDFX+POI(I)#POJ (1) *BFK(1)
5@ CONTINUE
RETURN
END

Eu

SUBROUTINE GRID (NI,NJ,FC,GC,BFI,BFU,POI,POU,XX, YY, EPS3,CDFX, IFLAG

C—— COMPUTE DOUBLE SUMMATION OF COMPONENTS OF THE T" C.D.F. OVER THE
C—— GRID I=IMIN TO IMAX AND J=JMIN TO JMAX

DIMENSION BFI(*) ,BFJ(*),POI(*) ,POU(*)
C-—- COMPUTE BETA C.D.F. BY RECURRENCE WHEN I=IMIN, J=JMIN TO JMAX +

ore EDGET (NU,GC,FC,YY,XX,BFu ,CDFX,POJ,POI,EPS3, IFLAG, 1)
IF (NI.LE.1.0R.IFLAG.NE.@) RETURN

Cc
C—— COMPUTE BETA C.D.F. BY RECURRENCE WHEN J=JMIN, I=IMIN TO IMAX
c

BFI(1) = BFU(1)
CALL EDGET (NI,FC,GC,XX,YY,BFI ,CDFX,POI ,POJ,EPS3, IFLAG,2)
IF (NJ.LE.1.0R. IFLAG.NE.Q) RETURN


*CDFDNT

SUBROUTINE CDFDNT (X,DF,DELTA,ALAMB, EPS, IFLAG,CDFX)

CDFONT 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 T DISTRIBUTION TO A SPECIFIED ACCURACY (TRUNCATION
ERROR IN THE INFINITE SERIES REPRESENTATION GIVEN BY EQUATION
4 IN REFERENCE 1 BELOW). WHEN X<@ THE C.D.F. IS COMPUTED
FROM CDF(X,DF,DELTA,ALAMB) = 1 — CDF(-X,DF,-DELTA,ALAMB) .

THE BETA C.D.F. ROUTINE IS CALLED AT MOST FOUR 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 T 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 100
FOR DELTA AND 10,00@ FOR LAMBDA CAN BE HANDLED WITH THE
CURRENT DIMENSION LIMITS. THE COMPUTED VALUE CDFX IS VALID
ONLY IF IFLAG=@ ON RETURN.

NOTE: IN SUBROUTINE EDGET THE DOUBLE PRECISION CONSTANT DEUFLO IS

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

SUBPROGRAMS CALLED: CDFBET (BETA C.D.F.)

GAMLOG (DOUBLE PRECISION LOG OF GAMMA FUNCTION)
POISST, EDGET, GRID (ATTACHED)

CURRENT VERSION COMPLETED SEPTEMBER 29, 1988

REFERENCES:

1.

KRISHNAN, MARAKATHA, ‘SERIES REPRESENTATIONS OF THE DOUBLY
NONCENTRAL T DISTRIBUTION’, JOURNAL OF THE AMERICAN STATISTICAL
ASSOCIATION, SEPTEMBER 1968, VOLUME 63, NO. 323, PP. 1004-1012.

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 T C.D.F. TO A SPECIFIED ACCURACY’, STATISTICAL
ENGINEERING DIVISION NOTE 86-5, DECEMBER 1986.

0000000000000000000000000090000090090990990902099998999999GH9999A9D99999G9G99E0

DEFINITION OF PASSED PARAMETERS:

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

* DELTA = THE NONCENTRALITY PARAMETER FOR THE NUMERATOR (REAL)

(EQUAL TO ZERO FOR THE CENTRAL T DISTRIBUTION]

* ALAMB = THE NONCENTRALITY PARAMETER (>=@) FOR THE DENOMINATOR

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

* EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (REAL)

[1 >= EPS >= 10++#(-10)]

IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION:

@ —> NO ERRORS DETECTED

2 -—> ERROR FLAGS FROM SUBROUTINE CDFBET

3 -—> ALAMB IS < @

4—-> DF IS <= @

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

6 —> VECTOR DIMENSIONS ARE TOO SMALL - INCREASE NX

1,

CDFX = THE DOUBLY NONCENTRAL T C.D.F. EVALUATED AT X (REAL)


8.

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. and Amos, D.E., "A Note on Representations of the Doubly
Non-Central t Distribution", Journal of the American Statistical
Association, Vol. 63, No. 323, September 1968, pp- 1013-1019.

Bulgren, W.G., "Probability Integral of the Doubly Noncentral t-
Distribution with Degrees of Freedom n and Non-Centrality Parameters 6
and A", Selected Tables in Mathematical Statistics, Vol. II, 1974,
pp. 1-138.

Carey, Michele B., "Evaluation of the Doubly Noncentral t Cumulative
Distribution Function", Computer Science and Statistics - The Interface,
James E. Gentle (ed.), North Holland Publishing Company, 1983,

Pp» 339-343.

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

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

Krishnan, Marakatha, "Series Representations of the Doubly Noncentral t-
Distribution", Journal of the American Statistical Association, Vol. 63,
No. 323, September 1968, pp. 1004-1012.

Mudholkar, Govind S. and Chaubey, Yogendra P., "A Simple Approximation for
the Doubly Noncentral t-Distribution", Communications in Statistics -
Simulation and Computation, Vol. B5, Nos. 2&3, 1976, pp» 85-92.

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

Patnaik, P.B., "Hypotheses Concerning the Means of Observations in Normal
Samples", Sankhya, Vol. 15, 1955, pp» 343-372.

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.

Robbins, Herbert, "The Distribution of Student's t When the Population
Means are Unequal", Annals of Mathematical Statistics, Vol. 19, No. 3,
September 1948, pp. 406-410.



L = K-IMIN+1
POI(L) = SIGN(SNGL(DEXP(DARG)) ,SDELTA)

SUM = POI(L)

DO 20 I = K-1, IMIN, ~1
L=tL-1
POI(L) = POI(L+1)#(REAL(I)+1.5)/FA
SUM = SUM+PO1(L)

28 CONTINUE
L = K-IMIN+1
DO 30 I = K+1, IMIN+NI-1

L = Lt
POI(L) = POI(L-1)#FA/(REAL(I)+0.5)
SUM = SUM+POI(L)
30 CONTINUE
IF (YY.GE.1.0) GO TO 40
FC = 1.0+REAL(IMIN)
c
C-—— SUM THE TERMS CORRESPONDING TO ODD’ VALUES OF INDEX I
c

CALL GRID (NI,NJ,FC,GC,8FI,8FJ,POI,POJ,xX, YY, EPS3,CDFX, IFLAG)
IF (IFLAG.NE.@) RETURN
ENDIF

Cc
C—— COMPUTE THE NORMAL C.D.F. AT —SDELTA
c

4@ PHI = @.5*(1.@-SUM)

c
C—— COMPUTE THE DOUBLY NONCENTRAL T C.D.F. AT X, USING AN IDENTITY
c— IF x<0
Cc

CDFX = @.5*CDFX+PHI

5@ IF (LL) CDFX = 1.0-CDFX
RETURN
END

SUBROUTINE POISST (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/3, 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
DLIMIT = 1.0D@—2.@D@*DBLE(EPS)/3.@00

K = INT(ALAMB)

Leet
IF (ALAMB.EQ.0.0) THEN
PL=1.
ELSE
DAL = DBLE(ALAMB)
QBLECK)
SNGL(DEXP (DK*DLOG(DAL)—DAL-GAMLOG(REAL(K+1))))

K = ALAMB+PL/REAL(L)
NK = NV/2
NL = NK#1
DSUM = 0.0
10 IF (PL.LT.PK) THEN
NK = NKH
IF (NK.GT.NV) THEN
IFLAG = 6
RETURN
ENDIF
V(NK) = PK
DSUM = DSUM+DBLE(PK)
K = K+
IF (DSUM.GE.DLIMIT) GO TO 20
PK = ALAMBsPK/REAL (KH)
ELSE
NL = NL-1
V(NL) = PL
DSUM = DSUM+DBLE(PL)
b= tt
IF (DSUM.GE.DLIMIT) GO TO 20
PL = REAL(L)*PL/ALAMB
ENDIF
GO To 10
20 INC = NL-1
DO 30 I = NL, NK



National Bureau of Standards N-86-5
Statistical Engineering Division December 1986
Gaithersburg, MD 20899

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

Charles P. Reeve

Let Z be a normally distributed random variable with mean 8 and variance 1,
and X be a noncentral chi-squared random variable, independent of Z, with
degrees of freedom v>0 and noncentrality parameter 430. Then the random
variable

Y= 2/VX/v qd)

has the doubly noncentral t distribution, indicated by Y~ t"(v,6,A). This
distribution was introduced by Robbins [14] as the distribution of Student's t
s istic when the fions have unequal po ulation ns. It was later
used by Pati (11] in testing hypotheses concerning the standardized means
of nonhomogeneous normal populations.

Krishnan [7] gives a series representation for the t" cumulative distribution’
function (c.d.f.) in terms of incomplete beta functions. Alternative series
representations are given by Bulgren and Amos [2], Bulgren [3], Carey [4], ana
[7]. Approximations are given by Johnson and Kotz [6] and Mulholkar and
Chaubey [8]. Numerical examples of usage are given in [3] and {7].

The author has been unable to find any published algorithms for computing
exact vy. £' ced.f. although computer Programs have obviously been
used in generating published tables in [2,3,4,7,8]. The Purpose of this note
is to present an efficient algorithm for computing the t" c.d.f. to a
specified accuracy using exact formulas.

The algorithm uses the series representation in eq- (4) of [7] which can be
re-written

Fy) = JF A.B, I(u,1/2+4/2,v/24j)/2 + Loa. J p.c-1y4/2 (2)
: 32 ji . J: 1
j=0 i=0 jr0 i=0

where A, = (4/2)4e/2 pr 543), Bes (De 2/065 aay, a= x2/(x24),
and x>0. When x<0 the c.d.f. is computed from the relation Fy (x5v,8,A) =
1 - Fy (-x3v,-6,A). The A; are Poisson Probabilities, and

u
Ku,a,b) = f oat) ae ya¢a,b) is the c.d.f. of the beta distribution
(also called the incomplete beta ratio) where O<u<1, a>0, b>0, and

1 *
B(a,b) = ff eT ety lac. The quantity (8//2)7 in Be is erroneously given
0


Cc
C—— COMPUTE BETA C.D.F. BY RECURRENCE WHEN I>IMIN, J>JMIN
Cc

DO 20 I = 2, NI
BFU(1) = BFI(I)
DO 10 J = 2, NU
BFU(J) = XX*BFU(J)+YY«BFu(J-1)
CDFX = CDFX+POI(I)*POU(J)*BFU(J)
10 CONTINUE
2@ CONTINUE
RETURN
END


* INDICATES PARAMETERS REQUIRING INPUT VALUES

2900

Cc
c—
Cc

Cc

1

PARAMETER (NX=1000)

DIMENSION BFI(NX) ,BFJ(NX) ,POI (NX) ,POJ(NX)
DOUBLE PRECISION DARG,DFA

LOGICAL LL

COFX = 0.0

CHECK VALIDITY OF ARGUMENTS.

IF (ALAMB.LT.0.@) THEN
IFLAG = 3
RETURN

ENDIF

IF (DF.LE.@.0) THEN
IFLAG = 4
RETURN

ENDIF

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

ENDIF

IFLAG = @

SET ERROR CRITERION FOR THE BETA C.D.F. (PECULIAR TO CDFBET)
EPS3 = @.001*EPS.

DELSQ = DELTA++2
FA = @.5*DELSQ
@.5*ALAMB
@.54DF
DF/(DF+X*X)
XX = 1.0-YY

IF X<@ SET LL=.TRUE., REVERSE SIGN OF DELTA, AND USE THE
IDENTITY DESCRIBED UP FRONT FOR COMPUTING THE C.D.F.

LL = X.LT.0.0
IF (XX.GE.1.0) THEN

CDFX = 1.0
GO TO 5@
ENDIF

SDELTA = DELTA
IF (LL) SDELTA = —DELTA

COMPUTE POISSON PROBABILITIES IN VECTOR POI

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

IF (YY.GE.1.0) GO To 10

FC = @.5+REAL(IMIN)

COMPUTE POISSON PROBABILITIES IN VECTOR POJ

CALL POISST (GA,EPS,JMIN,NJ,POJ,NX, IFLAG)
IF (IFLAG.NE.@) RETURN
GC = GB+REAL(JMIN)

SUM THE TERMS CORRESPONDING TO *EVEN’ VALUES OF INDEX I

CALL GRID (NI,NJ,FC,GC,BFI,BFJ,POI ,POJ,xXX, YY, EPS3,CDFX, IFLAG)
IF (IFLAG,NE.@) RETURN
@ IF (DELTA.EQ.0.0) THEN
NI =@
SUM = 0.0
IF (YY.GE.1.@) GO To 40
LSE

C—— COMPUTE ’POISSON-LIKE’ PROBABILITIES IN VECTOR POI
Cc

K = INT(FA)
IF (IMIN.GT.@) THEN
IMIN = IMIN-1
NI = NIH
ENDIF yy

DFA = DBLE(FA) ———
LE(K)+@. 5D) *DLOG(DFA)-DFA-GAMLOG(REAL(K)+1.5)



The above algorithm has been incorporated into the FORTRAN subroutine CDFDNT.
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 CDFNDT these
routines are the subroutine CDFBET and the function GAMLOG as described in
Reeve [12,13]. Other routines which can be substituted for CDFBET and GAMLOG
are those in [5,9,10]. For all practical purposes the absolute error
criterion ¢ will be met if the beta c.d.f. routine is accurate to two or three
digits beyond e.

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 6 and A become large. The computing time and storage also
increase as € becomes small, but are unaffected by v. In table 1 the t"
c.d.f. is computed for selected parameter values. The computing time in CPU
seconds and the number of beta c.d.f.'s computed are included. The case
v=6=A=100 corresponds to the example in Carey [4] who defines A a bit
differently. i eri i ul ears we
sui for utatio n_ § and/o alues. Note that CDFDNT
required only 0.76 CPU seconds in this case. Were all 98,490 beta c.d.f.
evaluations done by separate calls to the beta c.d.f. routine, the CPU time
would have been at least 100 times greater. The computations in table 1 were
done on the CDC Cyber 180/855 computer at NBS.

The current dimension limits in CDFDNT allow values of 6 up to 100 and 4 up to

10,000 with € as small as 10°, but these limits could easily be increased by
the user. The limiting factor in using CDFDNT 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 figures 1 and 2 are actually stored.

The Poisson probabilities Agi sees Ayn and Bie eeey Bin are also stored as
e

°

are B° wey Bins

rt?
°

If A=0 then the doubly noncentral t reduces to the (singly) noncentral t, and
if 6=A=0 it reduces to the central t. In either case, CDFDNT will run almost
as efficiently as routines designed for those specific cases.

Portions of tables in [2,3,4,7,8] were reproduced by CDFDNT and agreed to
within roundoff error in each case.

A listing of CDFDNT is an appendix to this note. It is invoked by
CALL CDFDNT(X, DF, DELTA, ALAMB, EPS, IFLAG, CDFX)

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 CDFDNT 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.


recurrence relations

I(x,a,b) = Ux,a,b+l) - x°(1-x)>/ppB(a,b)], (5a)
I(x,a,b) = U(x,atl,b) + x2(1-x)?/[aB(a,b)], and (5b)
I(x,a,b) = xI(x,a-l,b) + (1-x)I(x,a,b-1) (5e)

as found in Abramowitz and Stegun [1]. Subject to the restrictions J's sess"
and IjSIksI", Je and Ik are chosen to maximize the magnitudes of the rightmost

terms in (5a) and (5b) respectively. In applying each of these two recurrence
relations only one direct evaluation of B(a,b) is necessary, computed by

InI'(a) + Inf(b) - Inf(a+b) an} it

o

where T(a) =f t dt . Further
0

values are easily computed from the identities B(a+l,b) = aB(a,b)/(a+b) and

B(a,b+l) = bB(a,b)/(a+b). In a similar fashion I(u,1+1) ,v/2+I%) and

B(a,b) =e

I(u, 1+1%,v/2+5") are computed directly, and the same procedure is followed
with similar restrictions on Je and Tk. These computations are illustrated in
in figure 2. The double summations in (4) are accumulated as the beta c.d.f.'s

are computed.

In figures 1 and 2 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 "c" being done last.

jes" -| aoc e coe ec jes" -| aoc e ce ec oe
-| aoc ee ¢ e oc -| ac ec ¢ ec oc
-| aoc ec ¢ e oc -| ac e ce ec oc
jad* -| x C oee © C Corre © iss -| Kk CGC eee tO C Corre c¢
-| ac e cc e oc -| ac ec ¢ eoc
-| ac coc ce e oc -| ac e coc e oc
jes' -| a b boy b b b j=s' -| a b boy b b b
isI* isi" isi! isi"
e °
Figure 1 Figure 2

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


as (82/2) 4/2 in [7]. In the latter form the quantity would be (incorrectly)
positive when 6 is negative and i is an odd integer.

Each summation over i in (2) can be split into two summations over even and

4 -87/2

2 i :
odd values of i. For 4=0,1,2,++- let By = B,; = (8 /2)"e /T(i+1) and

: 2
Bo = 8... = (8/12)(82/2) 22? /?

4 Qi+l /T(i+3/2). Then (2) takes the form

FG = J I ape I(u,1/2+i,v/2+j)/2

j=0 i=0
(3)
oe
+ fy ABE S1(u,L+i,v/24j)/2 + {1 -} Be} /2
j=0 i=0 i=0
~
since } Ay = } BY = 1. For computational purposes the infinite series
ja0 i=0
must be truncated, thus (3) is re-expressed as
yuo
F(x) = Jf A,BeM(u,1/2+i,v/2+5)/2
¥ szyt geq! Jt
jes3' isl
e
(4)
yu" aa
+ J J Apgr(u,isi,w/24j)/2 + {1 - By} /2 +R
jest dst) J ist!

where q, Tj, I", J', and J" are non-negative integers and R is the remainder.

If the beta c.d.f.'s are computed without error, it can easily be shown that

su i"
choosing I}, I}, I", J', and Js" such that J A p> 1-26/3, J By > 1-2e/3,
jes! ist}

and 1) = max{I}-1,0} yields 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 j and i over the largest of the Poisson probabilities

A, and BF respectively. It then follows that i also indexes over the Be which

are largest in absolute value.

The final task is to compute the (21N-1)-1)+2)(J"-5' 41) beta c.d.f.'s and the

summations. An efficient procedure for doing this is to first compute only
I(u,1/2+1! ,v/2+I*) and I(u,1/24+1*,v/2+5') directly, indicated by the symbols

"x" and "y" in figure 1. The remaining beta c-d.f's are computed using the