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