C #QRCHSA V1A 29-JUN-73. C LAST UPDATE: 19-JUL-73 C FUNCTION QRCHSA (RCHSQ, NFREE) C C *PURPOSE. C TO EVALUATE THE PROBABILITY OF A RANDOM VARIABLE WITH C A REDUCED CHI-SQUARED DISTRIBUTION WITH "NFREE" DEGREES OF C FREEDOM EXCEEDING THE REDUCED CHI-SQUARED VALUE "RCHSQ". C NOTE THAT A REDUCED CHI-SQUARED IS THE STANDARD CHI- C SQUARED STATISTIC DIVIDED BY ITS NUMBER OF DEGREES OF FREEDOM. C NOTE ALSO THAT THIS IS A PROBABILITY FUNCTION OF CUMULATIVE C TYPE AND NOT A DENSITY FUNCTION. C C *PARAMETERS: C C RCHSQ - REAL VALUE OF THE REDUCED CHISQUARED FOR WHICH THE C PROBABILITY OF EXCEEDING IS REQUIRED. C NFREE - THE INTEGER NUMBER OF BEGREES OF FREEDOM. C C INPUT PARAMETERS: C RCHSQ, NFREE C OUTPUT PARAMETERS: C (REAL FUNCTION VALUE ONLY) C C *METHOD. C THE METHOD IS BASED ON TWO STANDARD FORMULA FOR THE C CHI-SQUARE STATISTIC. SEE FOR EXAMPLE M.ABRAMOWITZ AND C I.STEGUN, "HANDBOOK OF MATHEMATICAL FUNCTIONS", FORMULA NUMBERS C 26.4.5 WHEN "NFREE" IS EVEN AND 26.4.4 WHEN "NFREE" IS ODD. C THE METHOD IS NOT PARTICULARLY SUITED TO THE CASE OF "NFREE" C LARGE ( E.G. >60 ). C C *ACCURACY. C THE FORMULA USED FOR "NFREE" ODD INVOLVES AN APPROXIMATION C WITH ABSOLUTE ERROR LESS THAN 1.5E-7 IN ADDITION THERE WILL C BE ROUNDING ERRORS WHICH WILL INCREASE WITH "NFREE". C THE ROUTINE HAS BEEN CHECKED AGAINST TABLES FOR "NFREE" C UP TO 30 BY TAKING RANDOM ENTRIES SPREAD OVER THE RANGES. IT C AGREES TO 5 DECIMAL PLACES. C C *RESTRICTIONS. C "NFREE" SHOULD BE GREATER THEN ZERO AND "RCHSQ" POSITIVE C IF NOT THE VALUE ZERO WILL BE RETURNED. C C *ERROR CONDITIONS. C C *NON STANDARD ROUTINES CALLED. C C *TYPICAL TIMES. C C *ORIGIN. C THIS ROUTINE HAS BEEN WRITTEN BY M.R.MANNING AFTER SOME C TROUBLE WAS EXPERIENCED WITH THE ROUTINE GIVEN BY P.R.BEVINGTON C IN "DATA REDUCTION AND ANALYSIS...". C C *COMMENTS. C THE FORMULA 26.4.7 OF ABRAMOWITZ AND STEGUN AS USED BY C BEVINGTON, SHOULD NOT BE SUMMED DIRECTLY. IT HAS LARGE ROUNDING C ERRORS DUE TO THE OSCILLATING SIGN OF THE TERMS. C C #END. C C SQRT(2./PI) DATA R2OPI/ 0.797884561/ C CONSTANTS FOR NORMAL PROBABILITY APPROXIMATION C SEE ABRAMOWITZ AND STEGUN. FORMULA 26.2.19 DATA D1/ 0.0498673470/ DATA D2/ 0.0211410061/ DATA D3/ 0.0032776263/ DATA D4/ 0.0000380036/ DATA D5/ 0.0000488906/ DATA D6/ 0.0000053830/ C QRCHSA= 0. IF ( NFREE .LE. 0) RETURN X2= RCHSQ*FLOAT(NFREE) Z= X2/2. IF ( Z .LT. 0.) RETURN SUM= 0. IMX= NFREE/2 IF (NFREE .NE. 2*IMX) GO TO 40 C C EVEN NUMBER OF DEGREES OF FREEDOM - FINITE SUM. C 20 TRM= 1. DO 30 I= 1, IMX SUM= SUM + TRM 30 TRM= TRM*Z/FLOAT(I) QRCHSA= SUM*EXP(-Z) RETURN C C ODD NUMBER OF DEGREES OF FREEDOM. C C CALCULATE SUM 40 IF (IMX .EQ. 0) GO TO 60 X= SQRT(X2) TRM= X DO 50 I= 1, IMX SUM= SUM + TRM 50 TRM= TRM*X2/FLOAT(2*I+1) SUM= R2OPI*EXP(-Z)*SUM C C CALCULATE NORMAL PROBABILITY TERM 60 Q= 0. IF (X .GT. 5.5) GO TO 70 Q= 1./(1.+X*(D1+X*(D2+X*(D3+X*(D4+X*(D5+X*D6) ) ) ) ) ) Q= Q**16 C 70 QRCHSA= SUM + Q C RETURN END