*DECK MDCHI C IMSL ROUTINE NAME - MDCHI C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - INVERSE CHI-SQUARED PROBABILITY DISTRIBUTION C FUNCTION C C USAGE - CALL MDCHI (P,DF,X,IER) C C ARGUMENTS P - INPUT PROBABILITY IN THE EXCLUSIVE RANGE C (0,1) C DF - INPUT NUMBER OF DEGREES OF FREEDOM. DF MUST BE C IN THE EXCLUSIVE RANGE (.5,200000.). C X - OUTPUT CHI-SQUARED VALUE, SUCH THAT A RANDOM C VARIABLE, DISTRIBUTED AS CHI-SQUARED WITH C DF DEGREES OF FREEDOM, WILL BE LESS THAN OR C EQUAL TO X WITH PROBABILITY P. C IER - ERROR PARAMETER. (OUTPUT) C TERMINAL ERROR C IER = 129 INDICATES THAT THE BOUNDS WHICH C ENCLOSED P COULD NOT BE FOUND WITHIN 20 C (NCT) ITERATIONS C IER = 130 INDICATES THAT AN ERROR OCCURRED C IN IMSL ROUTINE MDNRIS (P IS NOT A VALID C VALUE) C IER = 131 INDICATES THAT AN ERROR OCCURRED C IN IMSL ROUTINE MDCH C IER = 132 INDICATES THAT THE VALUE X COULD C NOT BE FOUND WITHIN 50 (ITMAX) ITERATIONS, C SO THAT THE ABSOLUTE VALUE OF P1-P WAS C LESS THAN OR EQUAL TO EPS. (P1 IS THE C CALCULATED PROBABILITY AT X, EPS = .0001) C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - H32/MDCH,MDNOR,MDNRIS,MERFI,MERRC=ERFC, C MGAMAD=DGAMMA,UERTST,UGETIO C - H36,H48,H60/MDCH,MDNOR,MDNRIS,MERFI, C MERRC=ERFC,MGAMA=GAMMA,UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE MDCHI (P,DF,X,IER) C DATA EPS/.0001/,ITMAX/50/,NCT/20/,NSIG/5/ C FIRST EXECUTABLE STATEMENT IC = 0 IER = 0 C ESTIMATE STARTING X CALL MDNRIS (P,XP,IER) IF (IER .NE. 0) GO TO 80 DFF = .22222222222222 DFF = DFF/DF D = SQRT(DFF) X = DF*(1. -DFF + XP*D)**3 C IS THE CASE ASYMPTOTICALLY NORMAL IF (DF .GE. 40.) GO TO 9005 C FIND BOUNDS (IN X) WHICH ENCLOSE P NCNT = 0 IST = 0 ISW = 0 DX = X * .125 5 IF (X) 10,15,20 10 X = 0.0 DX = -DX GO TO 20 15 DX = .1 20 CALL MDCH (X,DF,P1,IER) DX = DX + DX IF (IER .NE. 0) GO TO 85 CSS = X NCNT = NCNT + 1 IF (NCNT .GT. NCT) GO TO 90 IF (P1 - P) 25,9005,30 25 X = X + DX ISW = 1 IF (IST .EQ. 0) GO TO 5 GO TO 35 30 X = X - DX IST = 1 IF (ISW .EQ. 0) GO TO 5 XR = CSS XL = X GO TO 40 35 XL = CSS XR = X C PREPARE FOR ITERATION TO FIND X 40 EPSP = 10.**(-NSIG) IF (XL .LT. 0.) XL = 0.0 CALL MDCH (XL,DF,P1,IER) IF (IER .NE. 0) GO TO 85 FXL = P1 - P CALL MDCH (XR,DF,P1,IER) IF (IER .NE. 0) GO TO 85 FXR = P1-P IF (FXL*FXR .NE. 0.0) GO TO 45 X = XR IF (FXL .EQ. 0.0) X = XL GO TO 9005 45 IF (DF .LE. 2. .OR. P .GT. .98) GO TO 50 C REGULA FALSI METHOD X = XL + FXL*(XR-XL)/(FXL-FXR) GO TO 55 C BISECTION METHOD 50 X = (XL+XR) * .5 55 CALL MDCH (X,DF,P1,IER) IF (IER .NE. 0) GO TO 85 FCS = P1-P IF (ABS(FCS) .GT. EPS) GO TO 60 GO TO 9005 60 IF (FCS * FXL .GT. 0.0) GO TO 65 XR = X FXR = FCS GO TO 70 65 XL = X FXL = FCS 70 IF (XR-XL .GT. EPSP*ABS(XR)) GO TO 75 GO TO 9005 75 IC = IC+1 IF (IC .LE. ITMAX) GO TO 45 IER = 132 GO TO 9000 C ERROR RETURNED FROM MDNRIS 80 IER = 130 GO TO 9000 C ERROR RETURNED FROM MDCH 85 IER = 131 GO TO 9000 90 IER = 129 9000 CONTINUE CALL UERTST (IER,6HMDCHI ) 9005 RETURN END *DECK MDNOR C IMSL ROUTINE NAME - MDNOR C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - NORMAL OR GAUSSIAN PROBABILITY DISTRIBUTION C FUNCTION C C USAGE - CALL MDNOR (Y,P) C C ARGUMENTS Y - INPUT VALUE AT WHICH FUNCTION IS TO BE C EVALUATED. C P - OUTPUT PROBABILITY THAT A RANDOM VARIABLE C HAVING A NORMAL (0,1) DISTRIBUTION WILL BE C LESS THAN OR EQUAL TO Y. C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - MERRC=ERFC C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE MDNOR (Y,P) C SPECIFICATIONS FOR ARGUMENTS REAL P,Y C SPECIFICATIONS FOR LOCAL VARIABLES REAL SQR1D2 DATA SQR1D2/.70710678118655/ C FIRST EXECUTABLE STATEMENT P = .5 * ERFC(-Y*SQR1D2) RETURN END *DECK MERRC C IMSL ROUTINE NAME - MERRC=ERFC C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JUNE 1, 1981 C C PURPOSE - EVALUATE THE COMPLEMENTED ERROR FUNCTION C C USAGE - RESULT = ERFC(Y) C C ARGUMENTS Y - INPUT ARGUMENT OF THE COMPLEMENTED ERROR C FUNCTION. C ERFC - OUTPUT VALUE OF THE COMPLEMENTED ERROR C FUNCTION. C C PRECISION/HARDWARE - SINGLE/ALL C NOTE - ERFC MAY NOT BE SUPPLIED BY IMSL IF IT C RESIDES IN THE MATHEMATICAL SUBPROGRAM C LIBRARY SUPPLIED BY THE MANUFACTURER. C C REQD. IMSL ROUTINES - NONE REQUIRED C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C REAL FUNCTION ERFC(Y) C SPECIFICATIONS FOR ARGUMENTS REAL Y C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER ISW,I DIMENSION P(5),Q(3),P1(8),Q1(7),P2(5),Q2(4) REAL P,Q,P1,Q1,P2,Q2,XMIN,XLARGE,SSQPI,X, * RES,XSQ,XNUM,XDEN,XI C COEFFICIENTS FOR 0.0 .LE. Y .LT. C .477 DATA P(1)/-.44422647396874/, 1 P(2)/10.731707253648/, 2 P(3)/15.915606197771/, 3 P(4)/374.81624081284/, 4 P(5)/2.5612422994823E-02/ DATA Q(1)/17.903143558843/, 1 Q(2)/124.82892031581/, 2 Q(3)/332.17224470532/ C COEFFICIENTS FOR .477 .LE. Y C .LE. 4.0 DATA P1(1)/7.2117582508831/, 1 P1(2)/43.162227222057/, 2 P1(3)/152.98928504694/, 3 P1(4)/339.32081673434/, 4 P1(5)/451.91895371187/, 5 P1(6)/300.45926102016/, 6 P1(7)/-1.3686485738272E-07/, 7 P1(8)/.56419551747897/ DATA Q1(1)/77.000152935229/, 1 Q1(2)/277.58544474399/, 2 Q1(3)/638.98026446563/, 3 Q1(4)/931.35409485061/, 4 Q1(5)/790.95092532790/, 5 Q1(6)/300.45926095698/, 6 Q1(7)/12.782727319629/ C COEFFICIENTS FOR 4.0 .LT. Y DATA P2(1)/-.22695659353969/, 1 P2(2)/-4.9473091062325E-02/, 2 P2(3)/-2.9961070770354E-03/, 3 P2(4)/-2.2319245973418E-02/, 4 P2(5)/-2.7866130860965E-01/ DATA Q2(1)/1.0516751070679/, 1 Q2(2)/.19130892610783/, 2 Q2(3)/1.0620923052847E-02/, 3 Q2(4)/1.9873320181714/ C CONSTANTS DATA XMIN/1.0E-8/,XLARGE/5.6875E0/ C ERFC(XBIG) .APPROX. SETAP DATA XBIG/25.90625/ DATA SSQPI/.56418958354776/ C FIRST EXECUTABLE STATEMENT X = Y ISW = 1 IF (X.GE.0.0E0) GO TO 5 ISW = -1 X = -X 5 IF (X.LT..477E0) GO TO 10 IF (X.LE.4.0E0) GO TO 30 IF (ISW .GT. 0) GO TO 40 IF (X.LT.XLARGE) GO TO 45 RES = 2.0E0 GO TO 65 C ABS(Y) .LT. .477, EVALUATE C APPROXIMATION FOR ERFC 10 IF (X.LT.XMIN) GO TO 20 XSQ = X*X XNUM = P(5) DO 15 I = 1,4 XNUM = XNUM*XSQ+P(I) 15 CONTINUE XDEN = ((Q(1)+XSQ)*XSQ+Q(2))*XSQ+Q(3) RES = X*XNUM/XDEN GO TO 25 20 RES = X*P(4)/Q(3) 25 IF (ISW.EQ.-1) RES = -RES RES = 1.0E0-RES GO TO 65 C .477 .LE. ABS(Y) .LE. 4.0 C EVALUATE APPROXIMATION FOR ERFC 30 XSQ = X*X XNUM = P1(7)*X+P1(8) XDEN = X+Q1(7) DO 35 I=1,6 XNUM = XNUM*X+P1(I) XDEN = XDEN*X+Q1(I) 35 CONTINUE RES = XNUM/XDEN GO TO 55 C 4.0 .LT. ABS(Y), EVALUATE C MINIMAX APPROXIMATION FOR ERFC 40 IF (X.GT.XBIG) GO TO 60 45 XSQ = X*X XI = 1.0E0/XSQ XNUM = P2(4)*XI+P2(5) XDEN = XI+Q2(4) DO 50 I = 1,3 XNUM = XNUM*XI+P2(I) XDEN = XDEN*XI+Q2(I) 50 CONTINUE RES = (SSQPI+XI*XNUM/XDEN)/X 55 RES = RES*EXP(-XSQ) IF (ISW.EQ.-1) RES = 2.0E0-RES GO TO 65 60 RES = 0.0E0 65 ERFC = RES RETURN END *DECK MGAMA C IMSL ROUTINE NAME - MGAMA=GAMMA C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JUNE 1, 1981 C C PURPOSE - EVALUATE THE GAMMA FUNCTION C C USAGE - RESULT = GAMMA(X) C C ARGUMENTS X - INPUT ARGUMENT. C GAMMA IS SET TO MACHINE INFINITY, WITH THE C PROPER SIGN, WHENEVER C X IS ZERO, C X IS A NEGATIVE INTEGER, C ABS(X) .LE. XMIN, OR C ABS(X) .GE. XMAX. C XMIN IS OF THE ORDER OF 10**(-39) AND C XMAX IS AT LEAST 34.8. THE EXACT VALUES C OF XMIN AND XMAX MAY ALLOW LARGER RANGES C FOR X ON SOME COMPUTERS. C SEE THE PROGRAMMING NOTES IN THE MANUAL C FOR THE EXACT VALUES. C GAMMA - OUTPUT SINGLE PRECISION VALUE OF THE GAMMA C FUNCTION. C C PRECISION/HARDWARE - SINGLE/ALL C NOTE - GAMMA MAY NOT BE SUPPLIED BY IMSL IF C IT RESIDES IN THE MATHEMATICAL SUBPROGRAM C LIBRARY SUPPLIED BY THE MANUFACTURER. C C REQD. IMSL ROUTINES - UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C REMARKS AN ERROR MESSAGE PRINTED BY UERTST FROM GAMMA SHOULD C BE INTERPRETED AS FOLLOWS C IER - ERROR INDICATOR C TERMINAL ERROR C IER = 129 INDICATES THAT THE ABSOLUTE VALUE C OF THE INPUT ARGUMENT X WAS SPECIFIED C GREATER THAN OR EQUAL TO XMAX. GAMMA C IS SET TO MACHINE INFINITY. C IER = 130 INDICATES THAT THE INPUT ARGUMENT C X WAS SPECIFIED AS ZERO OR A NEGATIVE C INTEGER OR THAT THE ABSOLUTE VALUE OF THE C INPUT ARGUMENT WAS LESS THAN OR EQUAL TO C XMIN. GAMMA IS SET TO MACHINE INFINITY C WITH THE PROPER SIGN FOR THE GAMMA C FUNCTION. IF X IS ZERO OR AN EVEN C NEGATIVE INTEGER, GAMMA HAS A NEGATIVE C SIGN. OTHERWISE IT HAS A POSITIVE SIGN. C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C REAL FUNCTION GAMMA (X) C SPECIFICATIONS FOR ARGUMENTS REAL X C SPECIFICATIONS FOR LOCAL VARIABLES REAL P,Q,P4,BIG1,PI,XMIN,SIGN,Y,T,R,A,TOP, * DEN,B INTEGER I,IEND,IEND1,IEND2,IER,J LOGICAL MFLAG DIMENSION P(7),Q(6),P4(5) C COEFFICIENTS FOR MINIMAX C APPROXIMATION TO GAMMA(X), C 2.0 .LE. X .LE. 3.0 DATA P(1)/3.4109112397125E+01/, 1 P(2)/ -4.8341273405598E+01/, 2 P(3)/4.3005887829594E+02/, 3 P(4)/-5.5688734338586E+01/, 4 P(5)/2.0585220673644E+03/, 5 P(6)/7.7192407739800E-01/, 6 P(7)/-3.1721064346240E+00/ DATA Q(1)/2.4455138217658E+02/, 1 Q(2)/-1.0174768492818E+03/, 2 Q(3)/1.1615998272754E+03/, 3 Q(4)/2.0512896777972E+03/, 4 Q(5)/6.8080353498091E-01/, 5 Q(6)/-2.5386729086746E+01/ C COEFFICIENTS FOR MINIMAX C APPROXIMATION TO LN(GAMMA(X)), C 12.0 .LE. X DATA P4(1)/9.1893853320467E-01/, 1 P4(2)/8.3333333333267E-02/, 2 P4(3)/-2.7777776505151E-03/, 3 P4(4)/7.9358449768E-04/, 4 P4(5)/-5.82399983E-04/ DATA IEND/7/,IEND1/6/,IEND2/5/ DATA PI/3.1415926535898/ C GAMMA(XMIN) .APPROX. R1MACH(1) C GAMMA(BIG1) .APPROX. R1MACH(1) DATA XMIN/0.0/ DATA BIG1/177.803/ C FIRST EXECUTABLE STATEMENT IER = 0 MFLAG = .FALSE. T = X IF (ABS(T).GT.XMIN) GO TO 5 IER = 130 GAMMA = R1MACH(1) IF (T.LE.0.0) GAMMA = -R1MACH(1) GO TO 9000 5 IF (ABS(T).LT.BIG1) GO TO 10 IER = 129 GAMMA = R1MACH(1) GO TO 9000 10 IF (T.GT.0.0) GO TO 25 C ARGUMENT IS NEGATIVE MFLAG = .TRUE. T = -T R = AINT(T) SIGN = 1.0 IF (AMOD(R,2.0) .EQ. 0.0) SIGN = -1. R = T-R IF (R.NE.0.0) GO TO 20 IER = 130 GAMMA = R1MACH(1) IF (SIGN.EQ.-1.0) GAMMA = -R1MACH(1) GO TO 9000 C ARGUMENT IS NOT A NEGATIVE INTEGER 20 R = PI/SIN(R*PI)*SIGN T = T+1.0 C EVALUATE APPROXIMATION FOR GAMMA(T) C T .GT. XMIN 25 IF (T.GT.12.0) GO TO 60 I = T A = 1.0 IF (I.GT.2) GO TO 40 I = I+1 GO TO (30,35,50),I C 0.0 .LT. T .LT. 1.0 30 A = A/(T*(T+1.0)) T = T+2.0 GO TO 50 C 1.0 .LE. T .LT. 2.0 35 A = A/T T = T+1.0 GO TO 50 C 3.0 .LE. T .LE. 12.0 40 DO 45 J=3,I T = T-1.0 A = A*T 45 CONTINUE C 2.0 .LE. T .LE. 3.0 50 TOP = P(IEND1)*T+P(IEND) DEN = T+Q(IEND1) DO 55 J=1,IEND2 TOP = TOP*T+P(J) DEN = DEN*T+Q(J) 55 CONTINUE Y = (TOP/DEN)*A IF (MFLAG) Y = R/Y GAMMA = Y GO TO 9005 C T .GT. 12.0 60 TOP = ALOG(T) TOP = (T-1.5)*(TOP-1.0)+TOP-1.5 T = 1.0/T B = T*T Y = (((P4(5)*B+P4(4))*B+P4(3))*B+P4(2))*T+P4(1)+TOP Y = EXP(Y) IF (MFLAG) Y = R/Y GAMMA = Y GO TO 9005 9000 CONTINUE CALL UERTST(-IER,6H MGAMA) CALL UERTST(IER,6HGAMMA ) 9005 RETURN END