Directory of image this file is from
This file as a plain text file
C
C ..................................................................
C
C SUBROUTINE CDTR
C
C PURPOSE
C COMPUTES P(X) = PROBABILITY THAT THE RANDOM VARIABLE U,
C DISTRIBUTED ACCORDING TO THE CHI-SQUARE DISTRIBUTION WITH G
C DEGREES OF FREEDOM, IS LESS THAN OR EQUAL TO X. F(G,X), THE
C ORDINATE OF THE CHI-SQUARE DENSITY AT X, IS ALSO COMPUTED.
C
C USAGE
C CALL CDTR(X,G,P,D,IER)
C
C DESCRIPTION OF PARAMETERS
C X - INPUT SCALAR FOR WHICH P(X) IS COMPUTED.
C G - NUMBER OF DEGREES OF FREEDOM OF THE CHI-SQUARE
C DISTRIBUTION. G IS A CONTINUOUS PARAMETER.
C P - OUTPUT PROBABILITY.
C D - OUTPUT DENSITY.
C IER - RESULTANT ERROR CODE WHERE
C IER= 0 --- NO ERROR
C IER=-1 --- AN INPUT PARAMETER IS INVALID. X IS LESS
C THAN 0.0, OR G IS LESS THAN 0.5 OR GREATER
C THAN 2*10**(+5). P AND D ARE SET TO -1.7E38. 0
C IER=+1 --- INVALID OUTPUT. P IS LESS THAN ZERO OR
C GREATER THAN ONE, OR SERIES FOR T1 (SEE
C MATHEMATICAL DESCRIPTION) HAS FAILED TO
C CONVERGE. P IS SET TO 1.7E38. 0
C
C REMARKS
C SEE MATHEMATICAL DESCRIPTION.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C DLGAM
C NDTR
C
C METHOD
C REFER TO R.E. BARGMANN AND S.P. GHOSH, STATISTICAL
C DISTRIBUTION PROGRAMS FOR A COMPUTER LANGUAGE,
C IBM RESEARCH REPORT RC-1094, 1963.
C
C ..................................................................
C
SUBROUTINE CDTR(X,G,P,D,IER)
DOUBLE PRECISION XX,DLXX,X2,DLX2,GG,G2,DLT3,THETA,THP1,
1GLG2,DD,T11,SER,CC,XI,FAC,TLOG,TERM,GTH,A2,A,B,C,DT2,DT3,THPI
C
C TEST FOR VALID INPUT DATA
C
IF(G-(.5-1.E-5)) 590,10,10
10 IF(G-2.E+5) 20,20,590
20 IF(X) 590,30,30
C
C TEST FOR X NEAR 0.0
C
30 IF(X-1.E-8) 40,40,80
40 P=0.0
IF(G-2.) 50,60,70
50 D=1.7E38 0
GO TO 610
60 D=0.5
GO TO 610
70 D=0.0
GO TO 610
C
C TEST FOR X GREATER THAN 1.E+6
C
80 IF(X-1.E+6) 100,100,90
90 D=0.0
P=1.0
GO TO 610
C
C SET PROGRAM PARAMETERS
C
100 XX=DBLE(X)
DLXX=DLOG(XX)
X2=XX/2.D0
DLX2=DLOG(X2)
GG=DBLE(G)
G2=GG/2.D0
C
C COMPUTE ORDINATE
C
CALL DLGAM(G2,GLG2,IOK)
DD=(G2-1.D0)*DLXX-X2-G2*.6931471805599453 -GLG2
IF(DD-1.68D02) 110,110,120
110 IF(DD+1.68D02) 130,130,140
120 D=1.7E38 0
GO TO 150
130 D=0.0
GO TO 150
140 DD=DEXP(DD)
D=SNGL(DD)
C
C TEST FOR G GREATER THAN 1000.0
C TEST FOR X GREATER THAN 2000.0
C
150 IF(G-1000.) 160,160,180
160 IF(X-2000.) 190,190,170
170 P=1.0
GO TO 610
180 A=DLOG(XX/GG)/3.D0
A=DEXP(A)
B=2.D0/(9.D0*GG)
C=(A-1.D0+B)/DSQRT(B)
SC=SNGL(C)
CALL NDTR(SC,P,DUMMY)
GO TO 490
C
C COMPUTE THETA
C
190 K= IDINT(G2)
THETA=G2-DFLOAT(K)
IF(THETA-1.D-8) 200,200,210
200 THETA=0.D0
210 THP1=THETA+1.D0
C
C SELECT METHOD OF COMPUTING T1
C
IF(THETA) 230,230,220
220 IF(XX-10.D0) 260,260,320
C
C COMPUTE T1 FOR THETA EQUALS 0.0
C
230 IF(X2-1.68D02) 250,240,240
240 T1=1.0
GO TO 400
250 T11=1.D0-DEXP(-X2)
T1=SNGL(T11)
GO TO 400
C
C COMPUTE T1 FOR THETA GREATER THAN 0.0 AND
C X LESS THAN OR EQUAL TO 10.0
C
260 SER=X2*(1.D0/THP1 -X2/(THP1+1.D0))
J=+1
CC=DFLOAT(J)
DO 270 IT1=3,30
XI=DFLOAT(IT1)
CALL DLGAM(XI,FAC,IOK)
TLOG= XI*DLX2-FAC-DLOG(XI+THETA)
TERM=DEXP(TLOG)
TERM=DSIGN(TERM,CC)
SER=SER+TERM
CC=-CC
IF(DABS(TERM)-1.D-9) 280,270,270
270 CONTINUE
GO TO 600
280 IF(SER) 600,600,290
290 CALL DLGAM(THP1,GTH,IOK)
TLOG=THETA*DLX2+DLOG(SER)-GTH
IF(TLOG+1.68D02) 300,300,310
300 T1=0.0
GO TO 400
310 T11=DEXP(TLOG)
T1=SNGL(T11)
GO TO 400
C
C COMPUTE T1 FOR THETA GREATER THAN 0.0 AND
C X GREATER THAN 10.0 AND LESS THAN 2000.0
C
320 A2=0.D0
DO 340 I=1,25
XI=DFLOAT(I)
CALL DLGAM(THP1,GTH,IOK)
T11=-(13.D0*XX)/XI +THP1*DLOG(13.D0*XX/XI) -GTH-DLOG(XI)
IF(T11+1.68D02) 340,340,330
330 T11=DEXP(T11)
A2=A2+T11
340 CONTINUE
A=1.01282051+THETA/156.D0-XX/312.D0
B=DABS(A)
C= -X2+THP1*DLX2+DLOG(B)-GTH-3.951243718581427
IF(C+1.68D02) 370,370,350
350 IF (A) 360,370,380
360 C=-DEXP(C)
GO TO 390
370 C=0.D0
GO TO 390
380 C=DEXP(C)
390 C=A2+C
T11=1.D0-C
T1=SNGL(T11)
C
C SELECT PROPER EXPRESSION FOR P
C
400 IF(G-2.) 420,410,410
410 IF(G-4.) 450,460,460
C
C COMPUTE P FOR G GREATER THAN ZERO AND LESS THAN 2.0
C
420 CALL DLGAM(THP1,GTH,IOK)
DT2=THETA*DLXX-X2-THP1*.6931471805599453 -GTH
IF(DT2+1.68D02) 430,430,440
430 P=T1
GO TO 490
440 DT2=DEXP(DT2)
T2=SNGL(DT2)
P=T1+T2+T2
GO TO 490
C
C COMPUTE P FOR G GREATER THAN OR EQUAL TO 2.0
C AND LESS THAN 4.0
C
450 P=T1
GO TO 490
C
C COMPUTE P FOR G GREATER THAN OR EQUAL TO 4.0
C AND LESS THAN OR EQUAL TO 1000.0
C
460 DT3=0.D0
DO 480 I3=2,K
THPI=DFLOAT(I3)+THETA
CALL DLGAM(THPI,GTH,IOK)
DLT3=THPI*DLX2-DLXX-X2-GTH
IF(DLT3+1.68D02) 480,480,470
470 DT3=DT3+DEXP(DLT3)
480 CONTINUE
T3=SNGL(DT3)
P=T1-T3-T3
C
C SET ERROR INDICATOR
C
490 IF(P) 500,520,520
500 IF(ABS(P)-1.E-7) 510,510,600
510 P=0.0
GO TO 610
520 IF(1.-P) 530,550,550
530 IF(ABS(1.-P)-1.E-7) 540,540,600
540 P=1.0
GO TO 610
550 IF(P-1.E-8) 560,560,570
560 P=0.0
GO TO 610
570 IF((1.0-P)-1.E-8) 580,580,610
580 P=1.0
GO TO 610
590 IER=-1
D=-1.7E38 0
P=-1.7E38 0
GO TO 620
600 IER=+1
P= 1.7E38 0
GO TO 620
610 IER=0
620 RETURN
END