C C :---------------------: C : : C : 'INDROC' PROGRAM : JUNE 1993 C : : C :---------------------: C C C******************************************************************************* C * C THE PURPOSE OF THIS PROGRAM IS TO CALCULATE THE STATISTICAL * C SIGNIFICANCE OF APPARENT DIFFERENCES BETWEEN TWO BINORMAL ROC * C CURVES ESTIMATED FROM INDEPENDENT CONFIDENCE-RATING DATA SETS, * C USING ANY ONE OF THREE DISTINCT STATISTICAL TESTS: * C * C (1) A BIVARIATE CHI-SQUARE TEST OF THE SIMULTANEOUS * C DIFFERENCES BETWEEN THE A PARAMETERS AND BETWEEN THE * C B PARAMETERS OF THE TWO ROC CURVES. (NULL HYPOTHESIS: * C DATA SETS AROSE FROM SAME ROC CURVE-- THAT IS, A1=A2 * C AND B1=B2). SEE REFERENCE (1). * C * C (2) A UNIVARIATE Z-SCORE TEST OF THE DIFFERENCE BETWEEN * C THE AREAS UNDER THE TWO ROC CURVES. (NULL HYPOTHESIS: * C DATA SETS AROSE FROM ROC CURVES WITH EQUAL AREAS BENEATH * C THEM). SEE REFERENCE (2). * C * C (3) A UNIVARIATE Z-SCORE TEST OF THE DIFFERENCE BETWEEN * C THE TRUE-POSITIVE FRACTIONS ON THE TWO ROC CURVES AT * C A SELECTED FALSE-POSITIVE FRACTION. (NULL HYPOTHESIS: * C THE DATA SETS AROSE FROM ROC CURVES HAVING THE SAME TPF * C AT THE SELECTED FPF). SEE REFERENCE (3). * C * C *** NOTES *** ROC DATA SETS ARE POTENTIALLY CORRELATED WHEN THE SAME * C STIMULI ARE PRESENTED TO DIFFERENT OBSERVERS, OR TO * C A SINGLE OBSERVER UNDER DIFFERENT DISPLAY, DATA * C PROCESSING, OR OBSERVATION CONDITIONS; OR WHEN THE * C STIMULI ARE PAIRED WITH SOME COMMON VARIATION, AS WHEN * C TWO DIAGNOSTIC TESTS OR IMAGING PROCEDURES ARE * C PERFORMED ON THE SAME PATIENT SAMPLE. IN SUCH * C SITUATIONS, THIS PROGRAM IS INAPPROPRIATE; THE DATA * C SHOULD BE ANALYZED BY THE 'CORROC2' PROGRAM INSTEAD. * C * C IN ADDITION TO THE REQUIREMENT THAT THE TWO DATA SETS * C MUST NOT BE POTENTIALLY CORRELATED, THE CASES USED * C AS INPUT WITHIN EACH DATA SET MUST BE DISTINCT, * C BECAUSE THIS PROGRAM ASSUMES THAT THE CONFIDENCE * C RATINGS WITHIN EACH DATA SET ARISE FROM ONE OF * C TWO MULTINOMIAL DESTRIBUTIONS (CORRESPONDING * C TO 'ACTUALLY NEGATIVE' AND 'ACTUALLY POSITIVE' * C CASES). HENCE, MULTIPLE CONFIDENCE RATINGS * C FOR EACH CASE (FOR EXAMPLE, FROM REPLICATED READINGS) * C MUST NOT BE POOLED IN THE INPUT DATA. USE OF POOLED * C DATA WILL CAUSE THE PROGRAM TO OVERESTIMATE THE * C STATISTICAL SIGNIFICANCE OF ANY APPARENT DIFFERENCE * C BETWEEN THE ROC CURVES, THEREBY INVALIDATING THE * C STATISTICAL TEST. IF THE DATA SETS FOR THE TWO ROC * C CURVES MAY BE ASSUMED INDEPENDENT FROM EACH OTHER BUT * C MULTIPLE RATINGS OF EACH CASE ARE AVAILABLE FOR * C ESTIMATION OF EACH ROC, SINGLE RATINGS FOR EACH CASE * C (FOR EXAMPLE, FROM A GIVEN OBSERVER) SHOULD BE RUN * C SEPARATELY (FOR EXAMPLE, OBSERVER-BY-OBSERVER FOR THE * C TWO CONDITIONS), OR THE GENERAL APPROACH DESCRIBED * C SWETS AND PICKETT ('EVALUATION OF DIAGNOSTIC SYSTEMS: * C METHODS FROM SIGNAL DETECTION THEORY.' ACADEMIC PRESS, * C NEW YORK, 1982, CHAPTER 4) SHOULD BE EMPLOYED, WITH * C R-SUB-C SET EQUAL TO ZERO. * C * C * C THE APPROACH USED IN THIS PROGRAM (WHICH IS VALID ONLY FOR * C INDEPENDENT DATA SETS) INVOLVES CALCULATING MAXIMUM * C LIKELIHOOD ESTIMATES OF THE 'A' PARAMETER AND THE 'B' PARAMETER * C OF THE ROC CURVE ASSOCIATED WITH EACH RATING-SCALE DATA SET, * C TOGETHER WITH ESTIMATES OF THE VARIANCES AND COVARIANCE OF THESE * C PARAMETERS (SEE REFERENCE (4) BELOW). THE BIVARIATE STATISTICAL * C TEST IS THEN IMPLEMENTED USING THE METHODS DESCRIBED IN * C REFERENCE (1). THE UNIVARIATE AREA TEST IS CARRIED OUT BY * C APPROXIMATING THE VARIANCE OF THE AREA INDEX IN TERMS OF * C FIRST-ORDER PARTIAL DERIVATIVES OF THE INDEX WITH RESPECT TO * C THE 'A' AND 'B' PARAMETERS AND IN TERMS OF THE ESTIMATED * C VARIANCES AND COVARIANCE OF THE 'A' AND 'B' PARAMETERS FOR EACH * C CURVE; AN APPROXIMATE Z-SCORE STATISTIC IS THEN FORMED BY * C DIVIDING THE OBSERVED INDEX DIFFERENCE BY THE SQUARE ROOT OF THE * C SUM OF THE INDEX VARIANCES (SEE REFERENCE (2)). THE UNIVARIATE TPF * C TEST IS PERFORMED IN ESSENTIALLY THE SAME WAY, EXCEPT THAT THE * C DIFFERENCE BETWEEN THE STANDARD DEVIATES CORRESPONDING TO THE TPF * C VALUES -- RATHER THAN BETWEEN THE TPF VALUES THEMSELVES -- IS TESTED * C (SEE REFERENCE (3)). * C * C THIS PROGRAM WAS DESIGNED BY: * C CHARLES E. METZ AND HELEN B. KRONMAN * C DEPARTMENT OF RADIOLOGY * C AND THE FRANKLIN MCLEAN MEMORIAL * C RESEARCH INSTITUTE * C THE UNIVERSITY OF CHICAGO * C CHICAGO, ILLINOIS 60637 * C AND WRITTEN BY HELEN B. KRONMAN, PU-LAN WANG, AND JONG-HER SHEN. * C * C INQUIRIES OR COMMENTS REGARDING THIS PROGRAM SHOULD BE SENT * C TO C. E. METZ AT THE ADDRESS ABOVE. * C * C THE PARTS OF THIS PROGRAM ASSOCIATED WITH MAXIMUM LIKELIHOOD * C ESTIMATION OF THE 'A' AND 'B' PARAMETERS OF A BINORMAL ROC * C CURVE WERE TAKEN FROM THE PROGRAM 'RSCORE II' BY: * C DONALD D. DORFMAN * C DEPARTMENT OF PSYCHOLOGY * C THE UNIVERSITY OF IOWA * C IOWA CITY, IOWA 52242 * C WITH MINOR CHANGES. * C * C REFERENCES: * C (1) METZ, C.E. AND KRONMAN, H.B. - JOURNAL OF MATHEMATICAL * C PSYCHOLOGY, 22: 218-243, 1980. * C (2) SWETS, J. A. AND PICKETT, R. M. - EVALUATION OF DIAGNOSTIC * C SYSTEMS: METHODS FROM SIGNAL DETECTION THEORY. (NEW YORK: * C ACADEMIC PRESS), 1982. * C (3) MCNEIL, B.J. AND HANLEY, J.A. - MEDICAL DECISION MAKING, * C 4:137-150, 1984. * C (4) DORFMAN, D. D. AND ALF, E. - JOURNAL OF MATHEMATICAL * C PSYCHOLOGY, 6: 487-496, 1969. * C * C ACKNOWLEDGEMENT: * C THIS WORK WAS SUPPORTED BY THE U.S. DEPARTMENT * C OF ENERGY UNDER CONTRACTS DE-AC02-80EV10359 AND * C DE-AC02-82ER60033 AND UNDER GRANT DE-FG02-86ER60418. * C * C******************************************************************************* C******************************************************************************* C * C REQUIRED INPUT AND DESCRIPTION OF OUTPUT: * C * C------------------------------------------------------------------------------* C * C * C INPUT --- * C * C 1. TWO DATA SETS ENTERED SEQUENTIALLY, EACH ON FOUR LINES: * C LINE #1 AND #5: DESCRIPTION OF THE DATA SET (UP TO 80 * C CHARACTERS). * C LINE #2 AND #6: TOTAL NUMBER OF CATEGORIES AND CATEGORY * C REPRESENTING THE STRONGEST EVIDENCE THAT A * C SIGNAL IS PRESENT (E.G., "DEFINITELY POSITIVE").* C LINE #3 AND #7: TOTAL NUMBER OF ACTUALLY NEGATIVE CASES; * C THEN NUMBER OF RESPONSES IN EACH CATEGORY, * C IN THE ORDER OF THE CATEGORY NUMBERS, * C FOR ACTUALLY NEGATIVE CASES. * C LINE #4 AND #8: TOTAL NUMBER OF ACTUALLY POSITIVE CASES; * C THEN NUMBER OF RESPONSES IN EACH CATEGORY, * C IN THE ORDER OF THE CATEGORY NUMBERS, * C FOR ACTUALLY POSITIVE CASES. * C FORMAT: FREE INPUT FORMAT, I.E., INPUT DATA CAN START * C IN ANY COLUMN BETWEEN COLUMN 1 AND COLUMN 80. * C IF THERE IS MORE THAN ONE INPUT VALUE PER * C LINE (LINES #2 THROUGH 4 AND #6 THROUGH 8), THEY * C MUST BE SEPARATED BY AT LEAST ONE BLANK SPACE. * C 2. AN ALPHABETIC CODE WORD INDICATING THE STATISTICAL TEST DESIRED: * C FOR BIVARIATE TEST, ENTER: BIVARIATE * C FOR AREA TEST, ENTER: AREA * C FOR TPF TEST, ENTER: TPF * C FORMAT: FREE INPUT FORMAT (NOTE: THE COMPUTER READS ONLY THE * C FIRST CHARACTER OF THIS CODE WORD.) * C IF 'TPF' WAS ENTERED ABOVE, INPUT ON AN ADDITIONAL LINE THE * C FPF VALUE (GREATER THAN 0.0 AND LESS THAN 1.0) AT WHICH THE * C TPF'S ARE TO BE COMPARED. OMIT THIS LAST INPUT LINE IF THE * C BIVARIATE OR AREA TEST WAS REQUESTED. * C * C * C INPUT EXAMPLE FOR AREA TEST: * C ---------------------------------------------------------------- * C TYPICAL 5-CATEGORY DATA SET * C 5 5 * C 53 22 19 6 4 2 * C 49 4 6 3 7 29 * C ANOTHER 5-CATEGORY DATA SET * C 5 5 * C 53 18 22 8 3 2 * C 49 2 3 4 3 37 * C AREA * C ---------------------------------------------------------------- * C * C INPUT EXAMPLE FOR TPF TEST: * C ---------------------------------------------------------------- * C TYPICAL 5-CATEGORY DATA SET * C 5 5 * C 53 22 19 6 4 2 * C 49 4 6 3 7 29 * C ANOTHER 5-CATEGORY DATA SET * C 5 5 * C 53 18 22 8 3 2 * C 49 2 3 4 3 37 * C TPF * C 0.10 * C ---------------------------------------------------------------- * C * C REMARKS CONCERNING INPUT: * C 1. THE NUMBER OF CATEGORIES FOR EACH DATA SET CAN BE 3 TO 11 * C AND NEED NOT BE THE SAME FOR THE TWO DATA SETS. * C 2. THE SECOND INPUT VALUES ON LINE #2 AND LINE #6 INDICATE * C WHETHER THE HIGHEST OR LOWEST NUMBERED CATEGORY REPRESENTS * C THE STRONGEST EVIDENCE OF POSITIVITY (E.G., THAT * C 'ABNORMALITY' IS PRESENT). IF HIGHER-NUMBERED CATEGORIES * C REPRESENT STRONGER EVIDENCE OF POSITIVITY, THEN THIS VALUE * C SHOULD BE EQUAL TO THE TOTAL NUMBER OF CATEGORIES INPUT * C FIRST ON THE LINE. ALTERNATIVELY, IF LOWER-NUMBERED * C CATEGORIES REPRESENT STRONGER EVIDENCE OF POSITIVITY, THEN * C THIS VALUE SHOULD BE '1'. * C 3. THIS PROGRAM ALLOWS THE USE OF SEVERAL PAIRS OF DATA SETS * C AS INPUT; SIMPLY ENTER PAIRS OF DATA SETS SEQUENTIALLY. * C 4. THE PROGRAM CHECKS TO ENSURE THAT THE INPUT NUMBERS OF * C CASES AND THE SUMS OF THE INPUT NUMBERS OF RESPONSES ARE * C CONSISTENT. * C * C * C * C * C * C OUTPUT --- * C * C I. FOR EACH DATA SET: * C 1. HEADING, WITH DESCRIPTION OF DATA SET, NUMBER OF CATEGORIES, * C 'DIRECTION' OF RATING SCALE, AND NUMBERS OF CASES. * C 2. RATING DATA FROM INPUT. * C 3. ROC OPERATING POINTS CALCULATED DIRECTLY FROM THE INPUT DATA. * C 4. INITIAL ESTIMATES OF THE PARAMETER VALUES, WITH THE LOG * C LIKELIHOOD OF THE DATA AND THE 'CHI-SQUARE' GOODNESS-OF-FIT * C (IF APPROPRIATE). 'A' REPRESENTS THE 'Y-INTERCEPT' AND * C 'B' REPRESENTS THE 'SLOPE' OF THE ROC CURVE WHEN IT IS * C PLOTTED ON NORMAL-DEVIATE AXES. IN TERMS OF AN EFFECTIVE * C PAIR OF UNDERLYING NORMAL DISTRIBUTIONS, THE PARAMETER 'B' * C CAN BE INTERPRETED AS THE RATIO OF THE STANDARD DEVIATION * C OF THE ACTUALLY-NEGATIVE DISTRIBUTION TO THAT OF THE ACTUALLY- * C POSITIVE DISTRIBUTION, WHEREAS THE PARAMETER 'A' CAN BE * C INTERPRETED AS THE ABSOLUTE DIFFERENCE BETWEEN THE MEANS OF * C TWO DISTRIBUTIONS, RELATIVE TO THE STANDARD DEVIATION OF THE * C ACTUALLY-POSITIVE DISTRIBUTION. THE Z(K) ARE THE CUT-OFFS ON * C THE DECISION-VARIABLE AXIS THAT DEFINE THE CATEGORY BOUNDARIES; * C THEY ARE MEASURED FROM THE MEAN OF THE ACTUALLY-NEGATIVE * C DISTRIBUTION TOWARD THE MEAN OF THE ACTUALLY-POSITIVE * C DISTRIBUTION AND ARE EXPRESSED IN UNITS OF THE STANDARD * C DEVIATION OF THE EFFECTIVE NORMALLY-DISTRIBUTED ACTUALLY- * C NEGATIVE DISTRIBUTION. * C 5. NUMBER OF ITERATIONS EMPLOYED. * C 6. FINAL PARAMETER ESTIMATES, WITH LOG LIKELIHOOD OF THE DATA * C AND 'CHI-SQUARE' GOODNESS-OF-FIT (IF APPROPRIATE). * C 7. VARIANCE-COVARIANCE AND CORRELATION MATRICES FOR FINAL * C PARAMETER ESTIMATES. * C 8. ESTIMATE AND STANDARD DEVIATION OF THE AREA BENEATH THE FINAL * C ESTIMATE OF THE BINORMAL ROC CURVE WHEN IT IS PLOTTED ON * C LINEAR AXES (OFTEN CALLED 'A SUB Z'). * C 9. TRUE-POSITIVE FRACTION AS A FUNCTION OF FALSE-POSITIVE * C FRACTION FOR THE ESTIMATED ROC CURVE, WITH LOWER AND * C UPPER BOUNDS OF ASYMMETRIC 95% CONFIDENCE INTERVAL ESTIMATED * C BY ASSUMING THAT ROC PARAMETER ESTIMATES 'A' AND 'B' ARE * C JOINTLY NORMAL. * C 10. ESTIMATED EXPECTED OPERATING POINTS ON THE FITTED ROC CURVE, * C WITH LOWER AND UPPER BOUNDS (ALONG THE FITTED CURVE) OF THE * C ASYMMETRIC 95% CONFIDENCE INTERVAL FOR THOSE OPERATING POINT * C ESTIMATES, CALCULATED ASSUMING THAT ERRORS IN THE Z(K) ESTIMATES * C ARE NORMALLY DISTRIBUTED. ESTIMATES OF THE EXPECTED OPERATING * C POINTS ON THE FITTED CURVE, ON NORMAL-DEVIATE AXES, HAVE ABSCISSA * C VALUES = -Z(K) AND ORDINATE VALUES = A-B*Z(K). * C * C II. ONE OF THREE TEST-STATISTIC VALUES, WITH ASSOCIATED SIGNIFICANCE: * C 1. THE COMPUTED 'BIVARIATE CHI-SQUARE TEST' STATISTIC VALUE WITH * C ITS CORRESPONDING P-LEVEL. * C 2. THE COMPUTED 'AREA TEST' STATISTIC VALUE WITH CORRESPONDING * C TWO-TAILED AND ONE-TAILED P-LEVELS. * C 3. THE COMPUTED 'TRUE-POSITIVE FRACTION TEST' STATISTIC VALUE * C WITH CORRESPONDING TWO-TAILED AND ONE-TAILED P-LEVELS AT THE * C SPECIFIED FALSE-POSITIVE FRACTION. * C * C REMARKS CONCERNING OUTPUT: * C 1. IT IS CRUCIALLY IMPORTANT TO NOTE THAT THE LEVEL OF * C STATISTICAL SIGNIFICANCE CALCULATED BY THIS PROGRAM IS * C NOT -- REPEAT, IS NOT -- VALID IF EITHER RATING-DATA * C SET HAS BEEN FORMED BY POOLING REPEATED OBSERVATIONS OF THE * C SAME STIMULI. * C 2. THE PROGRAM WILL CHECK FOR DEGENERACY IN EACH DATA SET. * C DEGENERATE DATA SETS ARE THOSE FOR WHICH AN EXACT FIT TO * C THE OBSERVED OPERATING POINTS IS PROVIDED BY A BINORMAL * C ROC CURVE AND ASSOCIATED CUT-OFFS WITH ONE OR MORE INFINITE * C PARAMETER VALUES, SO THAT THE ITERATIVE CALCULATION CANNOT * C CONVERGE. IF A DATA SET IS DEGENERATE, THIS PROGRAM WILL * C OUTPUT A MESSAGE DESCRIBING THE KIND OF DEGENERACY FOUND * C AND THE CORRESPONDING EXACT-FIT ROC; ABORT EXECUTION OF THE * C CURVE-FITTING ROUTINE FOR THAT DATA SET; AND ABORT CALCULATIONS * C OF THE STATISTICAL SIGNIFICANCE OF DIFFERENCES BETWEEN THE * C ASSOCIATED PAIR OF DATA SETS. * C 3. THE 'CHI-SQUARE' GOODNESS-OF-FIT MEASURE IS NOT OUTPUT * C IF ANY EXPECTED CELL FREQUENCY IS LESS THAN 5, OR IF ONLY 3 * C RATING CATEGORIES WERE EMPLOYED. * C * C * C******************************************************************************* DIMENSION XP(12,12) COMMON/PARA/A,B,X(10) COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN COMMON/BLK2/FPF(12),TPF(12) COMMON/BLK3/PXBA(12),DXBA(12),DIFFN(11),DIFFSN(11),PX(12),DX(12), + FUNLIK COMMON/BLK4/AAA,BBB,PARX(10),XX(12,12) COMMON/BLK5/TCODE,FIXFP,AX,BX,AY,BY,AZX,AZY,VAX,VBX,VAY,VBY, + VABX,VABY,VAZX,VAZY INTEGER TCODE,IARY(11) C IEND=0 1000 IFLAG=0 IERROR=0 DO 1900 NSET=1,2 C C READ IN C CALL READIN(IERROR,IEND) IF(IEND.EQ.1 .OR. IERROR.GT.0)GO TO 2900 C C PRINT HEADINGS C CALL OUTINI(NSET) C C CHECK CONSISTENCY OF INPUT DATA C CALL CHKIN(IFLAG) IF(IFLAG.EQ.1)GO TO 1900 C C COLLAPSE DATA. C CALL CLAPSE C C CALCULATE OBSERVED OPERATING POINTS C CALL OBPNT C C CHECK DEGENERACY OF DATA SET. IF IT IS A DEGENERATE DATA SET, C WRITE OUT MESSAGE AND ABORT THE EXECUTION OF THIS DATA SET. C CALL CHECK(ICLASS,KAT,FPF,TPF,ICON,IARY) IF (ICLASS.EQ.0) GO TO 1250 CALL DEGENE(ICLASS,IARY,FPF,TPF) IFLAG=1 GO TO 1000 C C GET INITIAL ESTIMATE OF PARAMETERS C 1250 CALL INIFIT C C USE "METHOD OF SCORING" TO GET FINAL ESTIMATES C CALL MOSLOP(LL,XP) IF(LL.LE.100)GO TO 1260 IFLAG=1 GO TO 1900 C C PRINT FINAL RESULT C 1260 CALL OUTRSL(NSET,LL,XP) 1900 CONTINUE C C READ IN REQUIRED STATISTICAL TEST CODE WORD C CALL READT(IERROR) IF(IERROR.NE.0)GO TO 2900 C C COMPUTE VALUES OF CHI-SQUARE TEST STATISTIC, OR AREA Z-TEST C STATISTIC OR TRUE-POSITIVE Z-TEST STATISTIC AND CORRESPONDING C P-LEVEL. C IF(IFLAG.EQ.0)GO TO 1920 WRITE(6,1910) 1910 FORMAT(//,15X,'EXECUTION OF STATISTICAL TEST ABORTED.') GO TO 1000 1920 CALL STAOUT GO TO 1000 2900 STOP END C---------------------------------------------- SUBROUTINE READIN(IERROR,IEND) C---------------------------------------------- C C READ INPUT DATA C COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN COMMON/PASS/LSTRING(80),LINE(80),LENGTH DATA KBLANK/1H / C C READ DATA DESCRIPTION C 100 READ(5,101,END=9999)(LINE(I),I=1,80) 101 FORMAT(80A1) DO 1000 I=1,80 IF(LINE(I).EQ.KBLANK)GO TO 1000 IPNT=I GO TO 1010 1000 CONTINUE IPNT=1 1010 DO 1020 I=IPNT,80 IDENT(I-IPNT+1)=LINE(I) 1020 CONTINUE C C READ IN TOTAL NUMBER OF CATEGORIES AND CATEGORY REPRESENTING C STRONGEST EVIDENCE THAT A SIGNAL IS PRESENT C 1100 READ(5,101)(LINE(I),I=1,80) CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 1170 CALL TONUM(LENGTH,LSTRING,NVALUE,IERROR) IF(IERROR.GT.0 .OR. (NVALUE.LT.3.OR.NVALUE.GT.11))GO TO 1170 KAT=NVALUE C CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 1180 CALL TONUM(LENGTH,LSTRING,NVALUE,IERROR) IF(IERROR.GT.0 .OR. (NVALUE.NE.1.AND.NVALUE.NE.KAT))GO TO 1180 KSN=NVALUE CALL GETWRD(NOTGET) IF(NOTGET.EQ.0)GO TO 1180 GO TO 1200 1170 IERROR=1 WRITE(6,1171)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 1171 FORMAT(2X,80A1/1X,'--- INVALID TOTAL NUMBER OF CATEGORY ---') RETURN 1180 IERROR=1 WRITE(6,1181)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 1181 FORMAT(2X,80A1/1X,'--- INVALID CATEGORY NUMBER REPRESENTING', 1 ' STRONGEST EVIDENCE OF SIGNAL ---') RETURN C C READ IN NUMBER OF N, SN CASES AND NUMBER OF N, SN RESPONSES C FOR EACH CATEGORIES C 1200 DO 1250 I1=1,2 READ(5,101)(LINE(I),I=1,80) CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 1270 CALL TONUM(LENGTH,LSTRING,NVALUE,IERROR) IF(IERROR.GT.0)GO TO 1270 IF(I1.EQ.1)EMN=NVALUE IF(I1.EQ.2)EMS=NVALUE C DO 1240 I2=1,KAT CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 1280 CALL TONUM(LENGTH,LSTRING,NVALUE,IERROR) IF(IERROR.GT.0)GO TO 1280 IF(I1.EQ.1)CATN(I2)=NVALUE IF(I1.EQ.2)CATS(I2)=NVALUE 1240 CONTINUE 1250 CONTINUE CALL GETWRD(NOTGET) IF(NOTGET.EQ.0)GO TO 1280 RETURN 1270 WRITE(6,1271)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 1271 FORMAT(2X,80A1/1X,'--- INVALID TOTAL NUMBER OF CASES ---') RETURN 1280 IERROR=1 WRITE(6,1281)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 1281 FORMAT(2X,80A1/1X,'--- INVALID NUMBER OF N, SN RESPONSES ---') RETURN 9999 IEND=1 RETURN END C------------------------------------- SUBROUTINE READT(IERROR) C------------------------------------- C C READ IN STATISTICAL TEST REQUEST CODE WORD AND C FPF VALUE (IF THE CODE WORD INITIALIZED BY "T") C COMMON/BLK5/TCODE,FIXFP,AX,BX,AY,BY,AZX,AZY,VAX,VBX,VAY,VBY, + VABX,VABY,VAZX,VAZY COMMON/PASS/LSTRING(80),LINE(80),LENGTH INTEGER TCODE DIMENSION NLINE(80) DATA IBLANK/1H /,IDEC/1H./,IZERO/1H0/ DATA IB1/1HB/,IB2/1Hb/,IA1/1HA/,IA2/1Ha/,IT1/1HT/,IT2/1Ht/ C READ(5,1010)(LINE(I),I=1,80) 1010 FORMAT(80A1) WRITE(6,1020) 1020 FORMAT(1H1,31X,'STATISTICAL SIGNIFICANCE OF THE DIFFERENCE',/, 1 25X,'BETWEEN THE TWO ROC CURVES ACCORDING TO THE SELECTED', 2 ' TEST,',/,32X,'ASSUMING THE DATA', 3 ' SETS TO BE INDEPENDENT'//12X,'**********************', 4 '********************************************************', 5 '********') WRITE(6,1030) 1030 FORMAT(/12X,'STATISTICAL TEST TO BE EMPLOYED:') TCODE=99 CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 1050 IF(LSTRING(1).EQ.IB1 .OR. LSTRING(1).EQ.IB2)TCODE=1 IF(LSTRING(1).EQ.IA1 .OR. LSTRING(1).EQ.IA2)TCODE=2 IF(LSTRING(1).EQ.IT1 .OR. LSTRING(1).EQ.IT2)TCODE=3 IF(TCODE.LT.1 .OR. TCODE.GT.3)GO TO 1050 IF(TCODE.EQ.3)GO TO 2000 RETURN 1050 IERROR=1 WRITE(6,1055) 1055 FORMAT(///,12X,'*** INVALID INPUT OF STATISTICAL TEST', 1 ' CODE (TCODE), STATISTICAL TEST ABORTED ***') RETURN C C READ IN FPF VALUE AT WHICH THE TPF'S ARE TO BE COMPARED C 2000 READ(5,1010)(NLINE(I),I=1,80) DO 2010 I=1,80 IF(NLINE(I).EQ.IBLANK)GO TO 2010 IPNT=I GO TO 2020 2010 CONTINUE GO TO 2090 2020 IF(NLINE(IPNT).EQ.IZERO)GO TO 2040 IF(NLINE(IPNT).EQ.IDEC)GO TO 2050 GO TO 2090 2040 IPNT=IPNT+1 IF(NLINE(IPNT).NE.IDEC)GO TO 2090 2050 IPNT=IPNT+1 DO 2060 I=IPNT,80 LINE(I-IPNT+1)=NLINE(I) 2060 CONTINUE INIT=80-IPNT+2 DO 2070 I=INIT,80 LINE(I)=IBLANK 2070 CONTINUE CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 2090 CALL TONUM(LENGTH,LSTRING,NUM,IERROR) IF(IERROR.GT.0)GO TO 2090 FIXFP=FLOAT(NUM)/10**LENGTH IF(FIXFP.LE.0 .OR. FIXFP.GE.1.0)GO TO 2090 RETURN 2090 IERROR=1 WRITE(6,2095) 2095 FORMAT(///12X,'*** INVALID INPUT OF FPF VALUE AT WHICH THE TPF', + 1H','S ARE TO BE COMPARED, STATISTICAL TEST ABORTED ***') RETURN END C-------------------------------------------- SUBROUTINE GETWRD(NOTGET) C-------------------------------------------- C C GET A STRING FROM THE INPUT LINE C COMMON /PASS/LSTRING(80),LINE(80),LENGTH DATA IBLANK/1H / C NOTGET=0 LENGTH=0 I=1 DO 10 J=1,80 LSTRING(J)=IBLANK 10 CONTINUE C 20 IF(I .GT. 80) GO TO 90 IF(LINE(I) .NE. IBLANK) GO TO 40 I= I+1 GO TO 20 C 40 II= 1 50 LSTRING(II)= LINE(I) LENGTH=II I= I+1 IF(I .GT. 80) GO TO 60 IF(LINE(I).EQ.IBLANK) GO TO 70 II= II+1 GO TO 50 C 60 DO 65 J=1,80 LINE(J)=IBLANK 65 CONTINUE RETURN C 70 DO 80 J=I,80 LINE(J-I+1)=LINE(J) 80 CONTINUE INIT=80-I+2 DO 85 J=INIT,80 LINE(J)=IBLANK 85 CONTINUE RETURN 90 NOTGET=1 RETURN END C------------------------------------------------ SUBROUTINE TONUM(LEN,LINPUT,NUM,IERR) C------------------------------------------------ C C CONVERT CHARACTER STRING TO POSITIVE INTEGER C DIMENSION LINPUT(*) DATA IZERO/1H0/ C IERR= 0 IF(LEN .EQ. 0) GO TO 120 C C CONVERT TO NUMERICAL STRING C NUM= 0 K= 1 40 IK= IORD(LINPUT(K)) - IORD(IZERO) IF(IK .LT. 0 .OR. IK .GT. 9) GO TO 120 NUM= NUM*10 + IK K= K+1 IF(K .LE. LEN) GO TO 40 RETURN C 120 IERR= 1 RETURN END C------------------------------------ INTEGER FUNCTION IORD(LCH) C------------------------------------ C C CHANGE FROM CHARACTER TO ASCII CODE C C CHANGE : put in \\ to rep. \ for unix machines DIMENSION LSTD(95) DATA LSTD/1H ,1H!,1H",1H#,1H$,1H%,1H&,1H',1H(,1H),1H*,1H+,1H-, + 1H,,1H.,1H/,1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9, + 1H:,1H;,1H<,1H=,1H>,1H?,1H@,1HA,1HB,1HC,1HD,1HE,1HF, + 1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO,1HP,1HQ,1HR,1HS, + 1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H[,1H\\,1H],1H^,1H_,1H`, + 1Ha,1Hb,1Hc,1Hd,1He,1Hf,1Hg,1Hh,1Hi,1Hj,1Hk,1Hl,1Hm, + 1Hn,1Ho,1Hp,1Hq,1Hr,1Hs,1Ht,1Hu,1Hv,1Hw,1Hx,1Hy,1Hz, + 1H{,1H|,1H},1H~/ DO 20 I= 1,95 IF(LSTD(I).NE.LCH) GO TO 20 IORD=I RETURN 20 CONTINUE 50 WRITE(6,60) 60 FORMAT(2X,'--- CHARACTER NOT FOUND ---') RETURN END C--------------------------------------- SUBROUTINE OUTINI(NSET) C--------------------------------------- C C WRITE OUT HEADING, DATA DESCRIPTION, NUMBERS OF C CATEGORIES AND CASES, AND RESPONSE DATA FROM INPUT. C COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN DIMENSION RCATN(11),RCATS(11) c IF(NSET.EQ.2)GO TO 1040 WRITE(6,1030) 1030 FORMAT(1H1,33X,'I N D R O C (JUNE 1993 VERSION) :'// 1 22X,'MAXIMUM LIKELIHOOD ESTIMATION OF BINORMAL ROC CURVE ', 2 'PARAMETERS'/ 3 48X,'AND'/ 4 22X,'STATISTICAL TESTING OF DIFFERENCES BETWEEN FITTED ', 5 'ROC CURVES'/ 6 48X,'FOR'/ 7 22X,' STATISTICALLY INDEPENDENT RATING-DATA SETS'//) WRITE(6,1035)(IDENT(I),I=1,80),KAT,KSN 1035 FORMAT(20X,'CONDITION X: ',80A1/ 1 20X,'DATA COLLECTED IN ',I2,' CATEGORIES'/ 2 20X,'WITH CATEGORY ',I2,' REPRESENTING STRONGEST EVIDENCE', 3 ' OF POSITIVITY (E.G., THAT ABNORMALITY IS PRESENT).'/) GO TO 1050 C 1040 WRITE(6,1041) 1041 FORMAT(1H1,//) WRITE(6,1045)(IDENT(I),I=1,80),KAT,KSN 1045 FORMAT(20X,'CONDITION Y: ',80A1,/, 120X,'DATA COLLECTED IN ' 2,I2,' CATEGORIES'/ 3 20X,'WITH CATEGORY ',I2,' REPRESENTING STRONGEST EVIDENCE', 4 ' OF POSITIVITY (E.G., THAT ABNORMALITY IS PRESENT).'/) C 1050 WRITE(6,1055)EMN,EMS 1055 FORMAT(20X,'NO. OF ACTUALLY NEGATIVE CASES =',F6.0,5X, 1 'NO. OF ACTUALLY POSITIVE CASES =',F6.0) WRITE(6,1060) 1060 FORMAT(/20X,'RESPONSE DATA:') WRITE(6,1070) (I,I=1,KAT) 1070 FORMAT(20X,' CATEGORY',19X,11(I2,5X)) WRITE(6,1080) (CATN(I),I=1,KAT) 1080 FORMAT(20X,' ACTUALLY NEGATIVE CASES ',11(F6.0,1X)) WRITE(6,1090) (CATS(I),I=1,KAT) 1090 FORMAT(20X,' ACTUALLY POSITIVE CASES ',11(F6.0,1X)) IF(KSN.EQ.KAT)RETURN DO 1091 I=1,KAT RCATN(I)=CATN(KAT+1-I) 1091 RCATS(I)=CATS(KAT+1-I) DO 1092 I=1,KAT CATN(I)=RCATN(I) 1092 CATS(I)=RCATS(I) RETURN END C------------------------------------------- SUBROUTINE CHKIN(IERROR) C------------------------------------------- C C CHECK CONSISTENCY OF INPUT DATA C COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN IERROR=0 SUMR=0. SUMRP=0. DO 1100 I=1,KAT SUMR=SUMR+CATN(I) SUMRP=SUMRP+CATS(I) 1100 CONTINUE IF ((SUMR.EQ.EMN).AND.(SUMRP.EQ.EMS)) RETURN IF (SUMRP.EQ.EMS) GO TO 1130 WRITE(6,1110) SUMRP 1110 FORMAT(//,20X,'SUM OF NO. OF S+N CASES IN EACH CATEGORY=', 1 F5.0) IF (SUMR.NE.EMN) GO TO 1150 WRITE(6,1120) 1120 FORMAT(////,22X,'DATA ARE INCONSISTENT. TOTAL NO. OF S+N CASES', 1 /,22X,'IS NOT EQUAL TO THE SUM OF NO. OF RESPONSES IN', 2 /,22X,'EACH CATEGORY.') GO TO 1900 1130 WRITE(6,1140) SUMR 1140 FORMAT(//,20X,'SUM OF NO. OF N CASES IN EACH CATEGORY=', 1 F5.0////,22X,'DATA ARE INCONSISTENT. TOTAL NO. OF', 2 ' N CASES',/22X,'IS NOT EQUAL TO THE SUM OF NO. OF', 3 ' RESPONSES',/22X,'IN EACH CATEGORY.') GO TO 1900 1150 WRITE(6,1160) SUMR 1160 FORMAT(20X,'SUM OF NO. OF N CASES IN EACH CATEGORY=', 1 F5.0,////,22X,'DATA ARE INCONSISTENT. TOTAL NO. OF', 2 ' N CASES',/22X,'IS NOT EQUAL TO THE SUM OF NO. OF', 3 ' RESPONSES',/22X,'IN EACH CATEGORY. TOTAL NO. OF S+N', 4 ' CASES IS',/,22X,'NOT EQUAL TO THE SUM OF NO. OF', 5 ' RESPONSES IN',/22X,'EACH CATEGORY.') 1900 IERROR=1 RETURN END C--------------------------------- SUBROUTINE CLAPSE C--------------------------------- C C COLLAPSE DATA C COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN KAT1=KAT DO 1190 I=1,KAT IT=KAT-I+1 IF (CATN(IT).NE.0.OR.CATS(IT).NE.0) GO TO 1190 KAT1=KAT1-1 IF (IT.GT.KAT1) GO TO 1190 DO 1180 J=IT,KAT1 CATN(J)=CATN(J+1) CATS(J)=CATS(J+1) 1180 CONTINUE 1190 CONTINUE KAT=KAT1 RETURN END C--------------------------------- SUBROUTINE OBPNT C--------------------------------- C C CALCULATES OBSERVED OPERATING POINTS C COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN COMMON/BLK2/FPF(12),TPF(12) NN=KAT+1 FPF(NN)=0. TPF(NN)=0. FPF(KAT)=CATN(KAT) TPF(KAT)=CATS(KAT) DO 1200 I=2,KAT FPF(KAT-I+1)=FPF(KAT-I+2)+CATN(KAT-I+1) TPF(KAT-I+1)=TPF(KAT-I+2)+CATS(KAT-I+1) 1200 CONTINUE DO 1210 I=1,NN FPF(I)=FPF(I)/EMN TPF(I)=TPF(I)/EMS 1210 CONTINUE WRITE(6,1220) 1220 FORMAT(/20X,'OBSERVED OPERATING POINTS:') WRITE(6,1230) (FPF(NN-I+1),I=1,NN) 1230 FORMAT(20X,' FPF: ',12(F6.4,1X)) WRITE(6,1240) (TPF(NN-I+1),I=1,NN) 1240 FORMAT(20X,' TPF: ',12(F6.4,1X)) RETURN END C-------------------------------------- SUBROUTINE INIFIT C-------------------------------------- C C CALCULATES LEAST SQUARE ESTIMATES C COMMON/PARA/A,B,X(10) COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN COMMON/BLK2/FPF(12),TPF(12) DIMENSION XS(10) KK=KAT-1 DO 1260 I=1,KK X(I)=FPF(I+1) XS(I)=TPF(I+1) 1260 CONTINUE C C TEST FOR 1.'S AND CORRECT BY SUBTRACTING 1/2N. TEST FOR 0'S C AND CORRECT BY CHANGING TO 1/2N. THEN CALL SUBROUTINE ZDEV TO C TRANSFORM THE CUMULATIVE PROBABILITIES IN THE ARRAYS TO STANDARD C NORMAL DEVIATES. C Q=0. DO 1310 I=1,KK IF(X(I).EQ.0.)X(I)=1./(2.*EMN) IF(ABS(X(I)-1.).LT.1.0E-05)X(I)=1.-(1./(2.*EMN)) P=X(I) CALL ZDEV(P,Q,D,IER) X(I)=Q IF(XS(I).EQ.0.)XS(I)=1./(2.*EMS) IF(ABS(XS(I)-1.).LT.1.0E-05)XS(I)=1.-(1./(2.*EMS)) P=XS(I) CALL ZDEV(P,Q,D,IER) XS(I)=Q 1310 CONTINUE IZ=KK-1 DO 1320 I=1,IZ IF(X(IZ-I+1).LE.X(IZ-I+2))X(IZ-I+1)=X(IZ-I+2)+.1 1320 IF(XS(IZ-I+1).LE.XS(IZ-I+2))XS(IZ-I+1)=XS(IZ-I+2)+.1 C C CALCULATE LEAST SQUARES SOLUTIONS C SUMX=0. SUMY=0. SUMXY=0. SUMX2=0. XK=FLOAT(KK) DO 1330 I=1,KK SUMX=SUMX+X(I) SUMX2=SUMX2+X(I)**2 SUMY=SUMY+XS(I) SUMXY=SUMXY+XS(I)*X(I) 1330 CONTINUE XMEAN=SUMX/XK YMEAN=SUMY/XK B=(XK*SUMXY-SUMX*SUMY)/(XK*SUMX2-SUMX**2) A=YMEAN-B*XMEAN C C WRITE OUT INITIAL ESTIMATES OF PARAMETER VALUES C WRITE(6,1340)A,B 1340 FORMAT(/,27X,'INITIAL VALUES OF PARAMETERS:',/,20X,'A= ',F7.4, 14X,'B= ',F7.4) DO 1350 I=1,KK X(I)=-X(I) 1350 CONTINUE 1355 WRITE(6,1360) (X(I),I=1,KK) 1360 FORMAT(20X,'Z(K)= ',10(F7.4,3X)) RETURN END C----------------------------------------- SUBROUTINE MOSLOP(LL,XP) C----------------------------------------- DIMENSION XXDUM(200),ADJX(10),XP(12,12) INTEGER SWT COMMON/PARA/A,B,X(10) COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN COMMON/BLK3/PXBA(12),DXBA(12),DIFFN(11),DIFFSN(11),PX(12),DX(12), + FUNLIK COMMON/BLK4/AAA,BBB,PARX(10),XX(12,12) C NN=KAT+1 KK=KAT-1 SWT=0 LL=0 1050 LL=LL+1 C C GET INTEGRALS AND DENSITIES C CALL TERMS IF(LL.EQ.1 .AND. KAT.GT.3)CALL CHISQR C C CALCULATES PARTIAL DERIVATIVES C CALL PDEV C C INVERT MATRIX C DO 1135 I=1,NN DO 1135 J=1,NN 1135 XX(I,J)=-XX(I,J) CALL MSTR(XX,XXDUM,12,0,1) CALL SINV(XXDUM,NN,.001,IER) C C IF PROBLEMS ARE ENCOUNTERED IN MATRIX INVERSION, REDUCE C ADJUSTMENTS TO SOLUTION VECTOR AND GO ON TO TEST CHANGE C IN LOG-LIKELIHOOD FUNCTION C IF(IER.GE.0) GO TO 1137 ALPHA=-0.5*ABS(ALPHA) GO TO 1157 C 1137 IF(IER.GT.0)WRITE(6,1140)IER 1140 FORMAT(' LOSS OF SIGNIFICANCE. STEP ',I5,'+1') CALL MSTR(XXDUM,XP,12,1,0) IF(SWT.EQ.1)RETURN IF(KAT.LE.3)RETURN C C FORM SOLUTION VECTOR C ADJA=AAA*XP(1,1)+BBB*XP(1,2) DO 1145 I=1,KK 1145 ADJA=ADJA+PARX(I)*XP(1,I+2) ADJB=AAA*XP(2,1)+BBB*XP(2,2) DO 1150 I=1,KK 1150 ADJB=ADJB+PARX(I)*XP(2,I+2) DO 1155 I=1,KK ADJX(I)=AAA*XP(I+2,1)+BBB*XP(I+2,2) DO 1155 K=3,NN 1155 ADJX(I)=ADJX(I)+PARX(K-2)*XP(I+2,K) C C ITERATE C 1156 ALPHA=1.0 1157 A=A+ALPHA*ADJA B=B+ALPHA*ADJB DO 1160 I=1,KK 1160 X(I)=X(I)+ALPHA*ADJX(I) CALL TERMS IF (ALPHA.NE.1.0) GO TO 1175 C C CHECK FOR MAXIMIZATION OF LOG-LIKELIHOOD C IF (LL.EQ.1) GO TO 1175 TNET=ABS(ADJA)+ABS(ADJB) DO 1173 I=1,KK TNET=TNET+ABS(ADJX(I)) 1173 CONTINUE TNET=6.*TNET/FLOAT(1+KAT) IF (TNET.GT.0.001) GO TO 1175 SWT=1 GO TO 1050 C C CHECK THAT LOG-LIKELIHOOD FUNCTION VALUE INCREASES, OR AT LEAST C DOESN'T DECREASE BY MORE THAN 0.001%; OTHERWISE, HALVE THE C CORRECTION VECTOR UNTIL REQUIREMENT IS SATISFIED. C 1175 IF (ALPHA.EQ.1.0) FUNLI1=FUNLIK FUNLIK=0. DO 1180 I=1,KAT FUNLIK=FUNLIK+CATN(I)*ALOG(DIFFN(I))+CATS(I)*ALOG(DIFFSN(I)) 1180 CONTINUE IF (FUNLIK.GE.FUNLI1) GO TO 1185 FUNDIF=ABS((FUNLIK-FUNLI1)/FUNLI1) IF (FUNDIF.LE.1.0E-5) GO TO 1185 ALPHA=-0.5*ABS(ALPHA) GO TO 1157 1185 IF(LL.LE.100)GO TO 1050 WRITE(6,1200) 1200 FORMAT(20X,'DOES NOT CONVERGE -- 100 ITERATIONS') RETURN END C------------------------------------ SUBROUTINE TERMS C------------------------------------ C C CALCULATES INTEGRALS AND DENSITIES C COMMON/PARA/A,B,X(10) COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN COMMON/BLK3/PXBA(12),DXBA(12),DIFFN(11),DIFFSN(11),PX(12),DX(12), + FUNLIK C NN=KAT+1 PXBA(1)=0. PXBA(NN)=1. DXBA(1)=0. DXBA(NN)=0. PX(1)=0. PX(NN)=1. DX(1)=0. DX(NN)=0. DO 1055 I=2,KAT Y=X(I-1)*B-A CALL NDTR(Y,PXBA(I),DXBA(I)) 1055 CONTINUE DO 1060 I=2,KAT Y=X(I-1) CALL NDTR(Y,PX(I),DX(I)) 1060 CONTINUE DO 1065 I=1,KAT DIFFSN(I)=PXBA(I+1)-PXBA(I) DIFFN(I)=PX(I+1)-PX(I) IF (DIFFSN(I).LT.5.0E-8) DIFFSN(I)=5.0E-8 IF (DIFFN(I).LT.5.0E-8) DIFFN(I)=5.0E-8 1065 CONTINUE FUNLIK=0. DO 1070 I=1,KAT FUNLIK=FUNLIK+CATN(I)*ALOG(DIFFN(I))+CATS(I)*ALOG(DIFFSN(I)) 1070 CONTINUE RETURN END C--------------------------------------- SUBROUTINE PDEV C--------------------------------------- C C CALCULATES PARTIAL DERIVATIVES C COMMON/PARA/A,B,X(10) COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN COMMON/BLK3/PXBA(12),DXBA(12),DIFFN(11),DIFFSN(11),PX(12),DX(12), + FUNLIK COMMON/BLK4/AAA,BBB,PARX(10),XX(12,12) C C GET FIRST PARTIAL DERIVATIVES C WITH RESPECT TO A C AAA=0. DO 1065 I=2,KAT 1065 AAA=AAA-DXBA(I)*(CATS(I-1)/DIFFSN(I-1)-CATS(I)/DIFFSN(I)) C C WITH RESPECT TO B C BBB=0. DO 1070 I=2,KAT 1070 BBB=BBB+DXBA(I)*X(I-1)*(CATS(I-1)/DIFFSN(I-1)-CATS(I)/DIFFSN(I)) C C WITH RESPECT TO Z'S C DO 1075 I=2,KAT Q1=DXBA(I)*B*(CATS(I-1)/DIFFSN(I-1)-CATS(I)/DIFFSN(I)) Q2=DX(I)*(CATN(I-1)/DIFFN(I-1)-CATN(I)/DIFFN(I)) 1075 PARX(I-1)=Q1+Q2 C C GET EXPECTED SECOND AND MIXED PARTIAL DERIVATIVES C WITH RESPECT TO A C XX(1,1)=0. DO 1080 I=2,KAT 1080 XX(1,1)=XX(1,1)-DXBA(I)*((DXBA(I)-DXBA(I-1))/DIFFSN(I-1)- 1(DXBA(I+1)-DXBA(I))/DIFFSN(I)) XX(1,1)=EMS*XX(1,1) C C WITH RESPECT TO B C XX(2,2)=0. DO 1095 I=2,KAT D=X(I-1) IF(I.EQ.KAT)GO TO 1090 IF(I.EQ.2)GO TO 1085 DD=X(I-2) DDD=X(I) GO TO 1095 1085 DD=0. DDD=X(I) GO TO 1095 1090 DD=X(I-2) DDD=0. 1095 XX(2,2)=XX(2,2)-DXBA(I)*X(I-1)*((DXBA(I)*D-DXBA(I-1)*DD)/ 1DIFFSN(I-1)-(DXBA(I+1)*DDD-DXBA(I)*D)/DIFFSN(I)) XX(2,2)=EMS*XX(2,2) C C WITH RESPECT TO A AND B C XX(1,2)=0. DO 1100 I=2,KAT 1100 XX(1,2)=XX(1,2)+DXBA(I)*X(I-1)*((DXBA(I)-DXBA(I-1))/ 1DIFFSN(I-1)-(DXBA(I+1)-DXBA(I))/DIFFSN(I)) XX(1,2)=EMS*XX(1,2) XX(2,1)=XX(1,2) C C WITH RESPECT TO A AND Z'S C DO 1105 I=2,KAT XX(1,I+1)=(EMS*B*DXBA(I)*((DXBA(I)-DXBA(I-1))/DIFFSN(I-1)- 1(DXBA(I+1)-DXBA(I))/DIFFSN(I))) 1105 XX(I+1,1)=XX(1,I+1) C C WITH RESPECT TO B AND Z'S C DO 1120 I=2,KAT XIL2=0. IF(I.GT.2)XIL2=X(I-2) IF(I.EQ.KAT)GO TO 1110 CHNG=X(I) GO TO 1115 1110 CHNG=0.0 1115 XX(2,I+1)=-(EMS*DXBA(I)*B*((DXBA(I)*X(I-1)-DXBA(I-1)*XIL2)/ 1DIFFSN(I-1)-(DXBA(I+1)*CHNG-DXBA(I)*X(I-1))/DIFFSN(I))) 1120 XX(I+1,2)=XX(2,I+1) C C WITH RESPECT TO Z'S AND MIXED WITH RESPECT TO Z'S C DO 1130 I=2,KAT IF(I.EQ.KAT)GO TO 1125 XX(I+1,I+2)=(EMS*DXBA(I)*DXBA(I+1)*B**2)/DIFFSN(I)+(EMN* 1DX(I)*DX(I+1))/DIFFN(I) XX(I+2,I+1)=XX(I+1,I+2) 1125 XX(I+1,I+1)=-(EMS*(DXBA(I)*B)**2*(1./DIFFSN(I-1)+1./ 1DIFFSN(I))+EMN*DX(I)**2*(1./DIFFN(I-1)+1./DIFFN(I))) 1130 CONTINUE RETURN END C------------------------------------- SUBROUTINE CHISQR C------------------------------------- C C GET VALUE OF LOG L AND GOODNESS-OF-FIT CHI SQUARE C COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN COMMON/BLK3/PXBA(12),DXBA(12),DIFFN(11),DIFFSN(11),PX(12),DX(12), + FUNLIK C NN=KAT+1 1120 SUMN=0. SUMM=0. DO 1130 I=2,NN SUMM=SUMM+CATS(I-1)*ALOG(DIFFSN(I-1)) SUMN=SUMN+CATN(I-1)*ALOG(DIFFN(I-1)) 1130 CONTINUE FLOGL=SUMN+SUMM WRITE(6,1140)FLOGL 1140 FORMAT(20X,'LOGL= ',F10.4) C IF (KAT.GT.3)GO TO 2020 WRITE(6,2010) 2010 FORMAT(20X,'CHI-SQUARE GOODNESS OF FIT NOT CALCULATED'/20X, 1'BECAUSE THERE ARE NO RESIDUAL DEGREES OF FREEDOM WHEN ONLY 3', 2' CATEGORIES ARE EMPLOYED.') RETURN C 2020 DO 2110 I=1,KAT IF(CATN(I).GE.5 .AND. CATS(I).GE.5)GO TO 2110 WRITE(6,2105) 2105 FORMAT(20X,'CHI-SQUARE GOODNESS OF FIT NOT CALCULATED', 1 ' BECAUSE SOME EXPECTED CELL FREQUENCIES ARE LESS THAN 5.') RETURN 2110 CONTINUE C IDCH=KAT-3 SUMM=0. SUMN=0. DO 2130 I=2,NN SUMN=SUMN+EMS*(CATS(I-1)/EMS-DIFFSN(I-1))**2/DIFFSN(I-1) 2130 SUMM=SUMM+EMN*(CATN(I-1)/EMN-DIFFN(I-1))**2/DIFFN(I-1) CHI=SUMM+SUMN DEGCHI=FLOAT(IDCH) CALL MDCH(CHI,DEGCHI,PVAL,IER1) PVAL=1-PVAL WRITE(6,2140) CHI,IDCH,PVAL 2140 FORMAT(20X,'GOODNESS OF FIT --',/, 120X,'CHI SQUARE=',F10.4,2X,'WITH ',I2, 2' DEGREES OF FREEDOM,', 3X,'P=',F7.4) RETURN END C----------------------------------------- SUBROUTINE OUTRSL(NSET,LL,XP) C----------------------------------------- C C WRITE OUT FINAL VALUES C COMMON/PARA/A,B,X(10) COMMON/BLK1/IDENT(80),EMN,EMS,CATN(11),CATS(11),KAT,KSN COMMON/BLK5/TCODE,FIXFP,AX,BX,AY,BY,AZX,AZY,VAX,VBX,VAY,VBY, + VABX,VABY,VAZX,VAZY DIMENSION XP(12,12),CX(12,12), 1 FPVAL(26),TPVAL(26),BDLOW(26),BDUPP(26), 2 EFPF(10),ETPF(10),FPFUPP(10),TPFUPP(10),FPFLOW(10),TPFLOW(10) DATA FPVAL/0.005,0.01,0.02,0.03,0.04,0.05,0.06, 2 0.07,0.08,0.09,0.10,0.11,0.12, 3 0.13,0.14,0.15,0.20,0.25,0.30, 4 0.40,0.50,0.60,0.70,0.80,0.90,0.95/ C NN=KAT+1 KK=KAT-1 IF(KAT.LE.3)LL=0 WRITE(6,1340)LL 1340 FORMAT(//20X,'PROCEDURE CONVERGES AFTER ',I3,' ITERATIONS.') IF(KAT.GT.3)GO TO 1350 WRITE(6,1342) 1342 FORMAT(20X,'(INITIAL ESTIMATES PROVIDE AN EXACT FIT BECAUSE', 1 ' ONLY 3 CATEGORIES WERE EMPLOYED)') 1350 WRITE(6,1360)A,B 1360 FORMAT(//27X,'FINAL VALUES OF PARAMETERS:'/20X,'A= ',F7.4, 1 4X,'B= ',F7.4) WRITE(6,1365)(X(I),I=1,KK) 1365 FORMAT(20X,'Z(K)= ',10(F7.4,3X)) CALL CHISQR C C GET CORRELATION MATRIX ON FINAL ITERATION C DO 1400 I=1,NN DO 1400 J=1,NN CX(I,J)=XP(I,J)/SQRT(XP(I,I)*XP(J,J)) 1400 CONTINUE WRITE(6,1410) 1410 FORMAT(//,27X,'VARIANCE-COVARIANCE MATRIX:') WRITE(6,1420) (XP(1,J),J=1,NN) 1420 FORMAT(20X,'A',5X,12(F7.4,1X)) WRITE(6,1430) (XP(2,J),J=1,NN) 1430 FORMAT(20X,'B',5X,12(F7.4,1X)) DO 1440 I=3,NN II=I-2 1440 WRITE(6,1450) II,(XP(I,J),J=1,NN) 1450 FORMAT(20X,'Z(',I2,')',1X,12(F7.4,1X)) WRITE(6,1460) 1460 FORMAT(//27X,'CORRELATION MATRIX:') WRITE(6,1420) (CX(1,J),J=1,NN) WRITE(6,1430) (CX(2,J),J=1,NN) DO 1470 I=3,NN II=I-2 1470 WRITE(6,1450) II,(CX(I,J),J=1,NN) C C SUMMARY OF ROC CURVE C WORKV=1+B*B ADEV=A/SQRT(WORKV) CALL NDTR(ADEV,AREA,DENS) CA2=XP(1,1) CB2=XP(2,2) CAB=XP(1,2) VAREA=SQRT((DENS*DENS)*(CA2/WORKV+((A*B)**2)*CB2/(WORKV**3)- 1 2*A*B*CAB/(WORKV*WORKV))) WRITE(6,1490)AREA,VAREA 1490 FORMAT(//26X,'AREA = ',F7.4,7X,'STD. DEV.(AREA) = ',F7.4) C IF (NSET.EQ.2) GO TO 1500 AX=A BX=B VAX=XP(1,1) VBX=XP(2,2) VABX=XP(1,2) AZX=AREA VAZX=VAREA GO TO 1600 1500 AY=A BY=B VAY=XP(1,1) VBY=XP(2,2) VABY=XP(1,2) AZY=AREA VAZY=VAREA C C WITH GIVEN FALSE-POSITIVE FRACTIONS, COMPUTE THE TRUE- C POSITIVE FRACTIONS ON THE ESTIMATED ROC CURVE, WITH C LOWER AND UPPER BOUNDS OF THE ASYMMETRIC 95% CONFIDENCE C INTERVAL. C 1600 DO 2205 I=1,26 P=FPVAL(I) CALL ZDEV(P,Q,D,IER) TPDEV=A+B*Q SI=SQRT(Q*Q*CB2+CA2+2.*Q*CAB) BOUND1=TPDEV-1.96*SI BOUND2=TPDEV+1.96*SI CALL NDTR(TPDEV,P,D) TPVAL(I)=P CALL NDTR(BOUND1,P,D) BDLOW(I)=P CALL NDTR(BOUND2,P,D) BDUPP(I)=P 2205 CONTINUE WRITE(6,2210) (FPVAL(I),TPVAL(I),BDLOW(I),BDUPP(I),I=1,26) 2210 FORMAT(1H1,/,20X,'ESTIMATED BINORMAL ROC CURVE, WITH', 1' LOWER AND UPPER',/,20X,'BOUNDS ON ASYMMETRIC 95% CONFIDENCE', 2' INTERVAL FOR',/,20X,'TRUE-POSITIVE FRACTION AT EACH SPECIFIED ', 3/,20X,'FALSE-POSITIVE FRACTION:',///,23X,'FPF',6X,'TPF',8X, 4'(LOWER BOUND, UPPER BOUND)',/, 526(/,22X,F5.3,4X,F6.4,6X,'(',3X,F6.4,2X,',',3X,F6.4,3X,')')) C C COMPUTE THE LOWER AND UPPER BOUNDS (ALONG THE FITTED CURVE) C OF THE ASYMMETRIC 95% CONFIDENCE INTERVAL FOR EACH C EXPECTED OPERATING POINT ESTIMATED ON THE FITTED CURVE. C DO 2305 I=1,KK Q=X(I) CALL NDTR(Q,P,D) EFPF(I)=1-P QU=Q-1.96*SQRT(XP(I+2,I+2)) QL=Q+1.96*SQRT(XP(I+2,I+2)) CALL NDTR(QU,P,D) FPFUPP(I)=1.-P CALL NDTR(QL,P,D) FPFLOW(I)=1.-P Q=B*X(I)-A CALL NDTR(Q,P,D) ETPF(I)=1-P QU=Q-1.96*B*SQRT(XP(I+2,I+2)) QL=Q+1.96*B*SQRT(XP(I+2,I+2)) CALL NDTR(QU,P,D) TPFUPP(I)=1.-P CALL NDTR(QL,P,D) TPFLOW(I)=1.-P 2305 CONTINUE WRITE(6,2310) 2310 FORMAT(/////,20X,'ESTIMATES OF EXPECTED OPERATING POINTS ON FITT', 1 'ED ROC',/,20X,'CURVE, WITH LOWER AND UPPER BOUNDS OF', 2 ' ASYMMETRIC 95%',/,20X,'CONFIDENCE INTERVALS ALONG THE', 3 ' CURVE FOR THOSE POINTS: ', 4 ///13X,'EXPECTED OPERATING POINT',9X,'LOWER', 5 ' BOUND',10X,'UPPER BOUND'/17X,'( FPF , TPF )',11X, 6 '( FPF , TPF )',5X,'( FPF , TPF )',/) DO 2320 I=1,KK II=KK-I+1 WRITE(6,2315) EFPF(II),ETPF(II),FPFLOW(II),TPFLOW(II),FPFUPP(II), 1 TPFUPP(II) 2315 FORMAT(17X,'(',F6.4,',',1X,F6.4,')',11X,'(',F6.4,',',1X,F6.4,')', 1 5X,'(',F6.4,',',1X,F6.4,')') 2320 CONTINUE RETURN END C------------------------------------------- SUBROUTINE STAOUT C------------------------------------------- C C CALCULATE AND PRINT STATISTICAL TEST RESULTS C COMMON/BLK5/TCODE,FIXFP,AX,BX,AY,BY,AZX,AZY,VAX,VBX,VAY,VBY, + VABX,VABY,VAZX,VAZY INTEGER TCODE IF (TCODE.EQ.2) GO TO 1200 IF (TCODE.EQ.3) GO TO 1300 WORK1=VAX+VAY WORK2=VBX+VBY RHO=(VABX+VABY)/SQRT(WORK1*WORK2) CHISTA=(1./(1.-RHO*RHO))*((AX-AY)*(AX-AY)/WORK1+ 1 (BX-BY)*(BX-BY)/WORK2-2.*RHO*(AX-AY)*(BX-BY)/ 2 SQRT(WORK1*WORK2)) X2LEVL=EXP(-CHISTA/2.) WRITE(6,1110) CHISTA,X2LEVL 1110 FORMAT(16X,'CHI-SQUARE TEST FOR THE DIFFERENCE BETWEEN', 1' (AX,BX) AND (AY,BY)' 2 ,//,12X,'THE COMPUTED ',1H','BIVARIATE CHI-SQUAR', 3 'E TEST',1H',' STATISTIC (2DF) VALUE IS',F8.4,/,16X,'WITH', 4 ' A CORRESPONDING P-LEVEL OF',F7.4,'.') RETURN 1200 AZSTAT=(AZX-AZY)/SQRT(VAZX**2+VAZY**2) CALL NDTR(AZSTAT,P,D) IF (P.GT.0.5) P=1-P AZLEVL=2.*P WRITE(6,1210) AZSTAT,AZLEVL,P 1210 FORMAT(16X,'AREA (A SUB Z) TEST', 1 //,12X,'THE COMPUTED ',1H','AREA TEST',1H',' STATISTIC', 2 ' VALUE IS',F8.4,/,16X,'WITH A CORRESPONDING TWO-TAILED ', 3 'P-LEVEL OF',F7.4,' AND ONE-TAILED P-LEVEL OF',F7.4,'.') RETURN 1300 P=1.-FIXFP CALL ZDEV(P,ZQ,D,IER) TPX=BX*ZQ-AX TPY=BY*ZQ-AY VTPX=VAX+ZQ*ZQ*VBX-2.*ZQ*VABX VTPY=VAY+ZQ*ZQ*VBY-2.*ZQ*VABY TPSTAT=(TPX-TPY)/SQRT(VTPX+VTPY) CALL NDTR(TPSTAT,P,D) IF (P.GT.0.5) P=1-P PLEVEL=2.*P WRITE(6,1310) FIXFP,FIXFP,TPSTAT,PLEVEL,P 1310 FORMAT(16X,'TPF TEST AT FPF = ',F7.4, 1 //,12X,'AT FALSE-POSITIVE-FRACTION EQUAL TO ',F4.2,', ', 2 /,12X,'THE COMPUTED ',1H','TRUE-POSITIVE-FRACTION TEST', 3 1H',' STATISTIC VALUE IS',F8.4,/,16X,'WITH A CORRES', 4 'PONDING TWO-TAILED P-LEVEL OF',F7.4,' AND ONE-TAILED', 5 ' P-LEVEL OF',F7.4,'.') RETURN END C------------------------------------------------------- SUBROUTINE CHECK(ICLASS,KAT,FPF,TPF,ICON,IARY) C------------------------------------------------------- C C PURPOSE C CHECK DEGENERACY OF DATA SET. C C DESCRIPTION OF PARAMETERS C ICLASS - OUTPUT CLASS NUMBER OF DEGENERACY. (USED BY C SUBROUTINE 'DEGENE' IF NECESSARY) C KAT - NUMBER OF CATEGORIES OF DATA SET. C FPF,TPF - OBSERVED OPERATING POINTS. C ICON,IARY - STORAGE SPACE FOR INFORMATION IN MESSAGE. (USED BY C SUBROUTINE 'DEGENE' IF NECESSARY) C REAL FPF(*),TPF(*) INTEGER IARY(*) ICLASS=0 NN=KAT+1 ICON=0 DO 10 I=2,KAT T=FPF(I)*TPF(I)*(1.-FPF(I))*(1.-TPF(I)) IF (T.EQ.0.) GO TO 10 ICON=ICON+1 IARY(ICON)=I 10 CONTINUE IF (KAT.EQ.1) GO TO 20 IF (ICON.EQ.0) GO TO 40 I0=IARY(1)-1 I1=IARY(ICON)+1 IF (KAT.EQ.2.AND.ICON.NE.0) GO TO 30 IF (ICON.GE.2) GO TO 15 IF ((TPF(I0).NE.1.0.OR.TPF(I1).NE.0.0).AND. 1 (FPF(I0).EQ.1.0.AND.FPF(I1).EQ.0.0)) ICLASS=5 IF ((FPF(I0).NE.1.0.OR.FPF(I1).NE.0.0).AND. 1 (TPF(I0).EQ.1.0.AND.TPF(I1).EQ.0.0)) ICLASS=6 RETURN 15 IF ((TPF(IARY(1)).EQ.TPF(IARY(ICON))).AND. 1 (FPF(I0).EQ.1.AND.FPF(I1).EQ.0)) ICLASS=7 IF ((TPF(I0).EQ.1.AND.TPF(I1).EQ.0).AND. 1 (FPF(IARY(1)).EQ.FPF(IARY(ICON)))) ICLASS=8 RETURN 20 ICLASS=1 RETURN 30 ICLASS=2 RETURN 40 ICA=0 ICB=0 DO 50 I=1,NN IF (FPF(I).NE.0.AND.FPF(I).NE.1) ICA=ICA+1 IF (TPF(I).NE.0.AND.TPF(I).NE.1) ICB=ICB+1 50 CONTINUE I1=0 I2=0 DO 60 I=2,KAT IF (FPF(I).EQ.1.AND.TPF(I).NE.1) GO TO 55 IF (FPF(I).EQ.0.AND.TPF(I).NE.0) I2=I2+1 IF (I2.EQ.1) IARY(2)=I GO TO 60 55 I1=I1+1 IARY(1)=I 60 CONTINUE II1=0 II2=0 DO 70 I=2,KAT IF (FPF(I).NE.1.AND.TPF(I).EQ.1) GO TO 65 IF (FPF(I).NE.0.AND.TPF(I).EQ.0) II2=II2+1 IF (II2.EQ.1) IARY(4)=I GO TO 70 65 II1=II1+1 IARY(3)=I 70 CONTINUE IF ((ICA.NE.0.AND.ICB.NE.0).OR.((ICA.EQ.0.AND. 1(I1.EQ.0.OR.I2.EQ.0)).OR.(ICB.EQ.0.AND.(II1.EQ.0.OR.II2.EQ.0)))) 2 GO TO 80 IF (ICA.EQ.0.AND.I1.GE.1.AND.I2.GE.1) GO TO 75 IF (ICB.EQ.0.AND.II1.GE.1.AND.II2.GE.1) ICLASS=4 IARY(1)=IARY(3) IARY(2)=IARY(4) RETURN 75 ICLASS=3 RETURN 80 KK1=IARY(3)+1 K1=IARY(1)+1 KK2=IARY(4)-1 K2=IARY(2)-1 IF ((II1.NE.0.AND.FPF(KK1).EQ.0).OR.(I2.NE.0.AND. 1 TPF(K2).EQ.1)) ICLASS=9 IF ((I1.NE.0.AND.TPF(K1).EQ.0).OR.(II2.NE.0.AND. 1 FPF(KK2).EQ.1)) ICLASS=10 RETURN END C--------------------------------------------------------- SUBROUTINE DEGENE(ICLASS,IARY,FPF,TPF) C--------------------------------------------------------- C C PURPOSE C OUTPUT MESSAGE FOR EACH DEGENERACY CLASS. C C DESCRIPTION OF PARAMETERS C ICLASS - INPUT CLASS NUMBER OF DEGENERACY FROM SUBROUTINE C 'CHECK'. C IARY - INPUT INFORMATION FROM SUBROUTINE 'CHECK'. C FPF,TPF - OBSERVED OPERATING POINTS. C REAL FPF(*),TPF(*) INTEGER IARY(*) IF (ICLASS.EQ.1) GO TO 10 IF (ICLASS.EQ.2) GO TO 20 IF (ICLASS.EQ.3) GO TO 30 IF (ICLASS.EQ.4) GO TO 40 IF (ICLASS.EQ.5.OR.ICLASS.EQ.7) GO TO 50 IF (ICLASS.EQ.6.OR.ICLASS.EQ.8) GO TO 60 IF (ICLASS.EQ.9) GO TO 90 WRITE(6,5) 5 FORMAT(///,22X, 1 'DATA ARE DEGENERATE AND IMPLY PERFECT (BUT PERVERSE)' 2 ,/,22X,'DECISION PERFORMANCE.') RETURN 10 WRITE(6,15) 15 FORMAT(///,22X, 1 'DATA ARE DEGENERATE AND PROVIDE NO OPERATING POINTS',/ 2 ,22X,'OFF THE (0,0) AND (1,1) CORNERS.') RETURN 20 P=FPF(IARY(1)) CALL ZDEV(P,QF,D,IER) P=TPF(IARY(1)) CALL ZDEV(P,QT,D,IER) DD=QT-QF WRITE(6,25) DD 25 FORMAT(///,22X, 1 'DATA ARE DEGENERATE AND PROVIDE ONLY A SINGLE ',/, 2 22X,'OPERATING POINT OFF THE (0,0) AND (1,1) CORNERS.', 3 /,22X,'BINORMAL ROC CURVE CANNOT BE ESTIMATED UNIQUELY. ', 4 /,22X,'UNIT SLOPE ROC CURVE (B=1) WOULD HAVE A=D',2H'=, 5 F8.4,'.') RETURN 30 WRITE(6,35) TPF(IARY(2)),TPF(IARY(1)) 35 FORMAT(///,22X, 1 'DATA ARE DEGENERATE AND PRODUCE NO OPERATING POINTS', 2 /,22X,'OFF THE BORDERS OF THE UNIT SQUARE. IMPLIED', 3 /,22X,'EXACT-FIT BINORMAL ROC CURVE IS HORIZONTAL AT', 4 /,22X,'UNDETERMINED HEIGHT BETWEEN TPF=',F8.4,' AND', 5 /,22X,'TPF=',F8.4,' .') RETURN 40 WRITE(6,45) FPF(IARY(2)),FPF(IARY(1)) 45 FORMAT(///,22X, 1 'DATA ARE DEGENERATE AND PRODUCE NO OPERATING POINTS',/ 2 ,22X,'OFF THE BORDERS OF THE UNIT SQUARE. IMPLIED ', 3 /,22X,'EXACT-FIT BINORMAL ROC CURVE IS VERTICAL AT ', 4 /,22X,'UNDETERMINED POSITION BETWEEN FPF=',F8.4, 5 /,22X,'AND FPF=',F8.4,' .') RETURN 50 WRITE(6,55) TPF(IARY(1)) 55 FORMAT(//,22X,'DATA ARE DEGENERATE. IMPLIED EXACT-FIT', 1 /,22X,'BINORMAL ROC CURVE IS HORIZONTAL AT CONSTANT', 2 /,22X,'TPF=',F8.4,' .') RETURN 60 WRITE(6,65) FPF(IARY(1)) 65 FORMAT(//,22X,'DATA ARE DEGENERATE. IMPLIED EXACT-FIT', 1 /,22X,'BINORMAL ROC CURVE IS VERTICAL AT CONSTANT', 2 /,22X,'FPF=',F8.4,' .') RETURN 90 WRITE(6,95) 95 FORMAT(///,22X, 1 'DATA ARE DEGENERATE AND IMPLY PERFECT DECISION', 2 /,22X,'PERFORMANCE.') RETURN END C----------------------------------------- SUBROUTINE ZDEV(P,X,D,IE) C----------------------------------------- C C PURPOSE C COMPUTES X = P**(-1) (Y), THE ARGUMENT X SUCH THAT Y=P(X)= THE C PROBABILITY THAT THE RANDOM VARIABLE U, DISTRIBUTED NORMALLY C N(0,1), IS LESS THAN OR EQUAL TO X. F(X), THE ORDINATE OF THE C NORMAL DENSITY, AT X, IS ALSO COMPUTED. C C DESCRIPTION OF PARAMETERS C P - INPUT PROBABILITY. C X - OUTPUT ARGUMENT SUCH THAT P=Y=THE PROBABILITY THAT U, C THE RANDOM VARIABLE, IS LESS THAN OR EQUAL TO X. C D - OUTPUT DENSITY, F(X). C IER - OUTPUT ERROR CODE C = -1 IF P IS NOT IN THE INTERVAL (0,1), INCLUSIVE. C X=D=.99999E+38 IN THAT CASE. C = 0 IF THERE IS NO ERROR. SEE REMARKS, BELOW. C C REMARKS C MAXIMUM ERROR IS 0.00045. C IF P=0, X IS SET TO -(10)**38. D IS SET TO 0. C IF P=1, X IS SET TO (10)**38. D IS SET TO 0. C NOTE: ORIGINAL PROGRAM SET X TO + OR -(10)**74. C C SUBROUTINES AND SUBPROGRAMS REQUIRED. C NONE C C METHOD C BASED ON APPROXIMATIONS IN C. HASTINGS, APPROXIMATIONS FOR C DIGITAL COMPUTERS, PRINCETON UNIV. PRESS, PRINCETON, N.J., 1955. C SEE EQUATION 26.2.23, HANDBOOK OF MATHEMATICAL FUNCTIONS, C ABRAMOWITZ AND STEGUN, DOVER PUBLICATIONS, INC., NEW YORK. C IE=0 C X=.99999E+74 X=.99999E+38 D=X IF(P)1,4,2 1 IE=-1 GO TO 12 2 IF(P-1.0)7,5,1 C4 X=-.99999E+74 4 X=-.99999E+38 5 D=0.0 GO TO 12 7 D=P IF(D-0.5)9,9,8 8 D=1.0-D 9 T2=ALOG(1.0/(D*D)) T=SQRT(T2) X=T-(2.515517+0.802853*T+0.010328*T2)/(1.0+1.432788*T+0.189269 1*T2+0.001308*T*T2) IF(P-0.5)10,10,11 10 X=-X 11 D=0.3989423*EXP(-X*X/2.0) 12 RETURN END C------------------------------------ SUBROUTINE NDTR(X,P,D) C------------------------------------ C C PURPOSE C COMPUTES Y=P(X)=PROBABILITY THAT THE RANDOM VARIABLE C DISTRIBUTED NORMALLY IN (0,1), IS LESS THAN OR EQUAL C TO X. F(X), THE ORDINATE OF THE NORMAL DENSITY AT X, C IS ALSO COMPUTED. C C DESCRIPTION OF PARAMETERS C X - INPUT SCALAR FOR WHICH P(X) IS COMPUTED. C P - OUTPUT PROBABILITY. C D - OUTPUT DENSITY. C C REMARKS C MAXIMUM ERROR IS .0000007. C C SUBROUTINE AND SUBPROGRAM REQUIRED C NONE C C METHOD C BASED ON APPROXIMATION IN C. HASTINGS, APPROXIMATIONS C FOR DIGITAL COMPUTERS, PRINCETON UNIV. PRESS, PRINCETON, C N.J., 1955. SEE EQUATION 26.2.17, HANDBOOK OF MATHEMATICAL C FUNCTIONS, ABRAMOWITZ AND STEGUN, DOVER PUBLICATION, INC., C NEW YORK. C AX=ABS(X) T=1.0/(1.0+0.2316419*AX) IF (AX.GT.18) GO TO 5 D=0.3989423*EXP(-X*X/2.0) GO TO 6 5 D=0.0 6 P=1.0-D*T*((((1.330274*T-1.821256)*T+1.781478)*T-0.3565638)*T+ 10.3193815) IF(X)1,2,2 1 P=1.0-P 2 RETURN END C----------------------------------------- SUBROUTINE MSTR(A,R,N,MSA,MSR) C----------------------------------------- C C PURPOSE C CHANGE STORAGE MODE OF A MATRIX C C DESCRIPTION OF PARAMETERS C A - NAME OF INPUT MATRIX C R - NAME OF OUTPUT MATRIX C N - NUMBER OF ROWS AND COLUMNS IN A AND R C MSA - ONE DIGIT NUMBER FOR STORAGE MODE OF MATRIX A C 0 - GENERAL C 1 - SYMMETRIC C 2 - DIAGONAL C MSR - SAME AS MSA EXCEPT FOR MATRIX R C C REMARKS C MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRIX A C MATRIX A MUST BE A SQUARE MATRIX C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C LOC C C METHOD C MATRIX A IS RESTRUCTURED TO FORM MATRIX R C MSA MSR C 0 0 MATRIX A IS MOVED TO MATRIX R C 0 1 THE UPPER TRIANGLE OF ELEMENTS OF A GENERAL MATRIX ARE C USED TO FORM A SYMMETRIC MATRIX C 0 2 THE DIAGONAL ELEMENTS OF A GENERAL MATRIX ARE USED TO FORM C A DIAGONAL MATRIX C 1 0 A SYMMETRIC MATRIX IS EXPANDED TO FORM A GENERAL MATRIX C 1 1 MATRIX A IS MOVED TO MATRIX R C 1 2 THE DIAGONAL ELEMENTS OF A SYMMETRIC MATRIX ARE USED TO C FORM A DIAGONAL MATRIX C 2 0 A DIAGONAL MATRIX IS EXPANDED BY INSERTING MISSING ZERO C ELEMENTS TO FORM A GENERAL MATRIX C 2 1 A DIAGONAL MATRIX IS EXPANDED BY INSERTING MISSING ZERO C ELEMENTS TO FORM A SYMMETRIC MATRIX C 2 2 MATRIX A IS MOVED TO MATRIX R C DIMENSION A(*),R(*) DO 20 I=1,N DO 20 J=1,N C C IF R IS GENERAL, FORM ELEMENT C IF(MSR)5,10,5 C C IF IN LOWER TRIANGLE OF SYMMETRIC OR DIAGONAL R, BYPASS C 5 IF(I-J)10,10,20 10 CALL LOC(I,J,IR,N,MSR) C C IF IN UPPER AND OFF DIAGONAL OF DIAGONAL R, BYPASS C IF(IR)20,20,15 C C OTHERWISE FORM R(I,J) C 15 R(IR)=0.0 CALL LOC(I,J,IA,N,MSA) C C IF THERE IS NO A(I,J), LEAVE R(I,J) AT 0 C IF(IA)20,20,18 18 R(IR)=A(IA) 20 CONTINUE RETURN END C-------------------------------------- SUBROUTINE LOC(I,J,IR,N,MS) C-------------------------------------- C C PURPOSE C COMPUTE A VECTOR SUBSCRIPT FOR AN ELEMENT IN A MATRIX OF C SPECIFIED STORAGE MODE. C C DESCRIPTION OF PARAMETERS C I - ROW NUMBER OF ELEMENT C J - COLUMN NUMBER OF ELEMENT C IR - RESULTANT VECTOR SUBSCRIPT C N - NUMBER OF ROWS IN MATRIX C MS - ONE DIGIT NUMBER FOR STORAGE MODE OF MATRIX C 0 - GENERAL C 1 - SYMMETRIC C 2 - DIAGONAL C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C MS=0 SUBSCRIPT IS COMPUTED FOR A MATRIX WITH N*N ELEMENTS IN C STORAGE (GENERAL MATRIX) C MS=1 SUBSCRIPT IS COMPUTED FOR A MATRIX WITH N*(N+1)/2 IN C STORAGE (UPPER TRIANGLE OF SYMMETRIC MATRIX). IF ELEMENT C IS IN LOWER TRIANGULAR PORTION, SUBSCRIPT IS CORRESPONDING C ELEMENT IN UPPER TRIANGLE. C MS=2 SUBSCRIPT IS COMPUTED FOR A MATRIX WITH N ELEMENTS IN C STORAGE (DIAGONAL ELEMENTS OF DIAGONAL MATRIX). IF C ELEMENT IS NOT ON DIAGONAL (AND THEREFORE NOT IN STORAGE), C IR IS SET TO ZERO. C IX=I JX=J IF(MS-1)10,20,30 10 IRX=N*(JX-1)+IX GO TO 36 20 IF(IX-JX)22,24,24 22 IRX=IX+(JX*JX-JX)/2 GO TO 36 24 IRX=JX+(IX*IX-IX)/2 GO TO 36 30 IRX=0 IF(IX-JX)36,32,36 32 IRX=IX 36 IR=IRX RETURN END C--------------------------------------- SUBROUTINE SINV(A,N,EPS,IER) C--------------------------------------- C C PURPOSE C INVERT A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX C C USAGE C CALL SINV(A,N,EPS,IER) C C DESCRIPTION OF PARAMETERS C A - UPPER TRIANGULAR PART OF THE GIVEN SYMMETRIC POSITIVE C DEFINITE N BY N COEFFICIENT MATRIX. ON RETURN A C CONTAINS THE RESULTANT UPPER TRIANGULAR MATRIX. C N - THE NUMBER OF ROW (COLUMNS) IN GIVEN MATRIX. C EPS - AN INPUT CONSTANT WHICH IS USED AS RELATIVE TOLERANCE C FOR TEST ON LOSS OF SIGNIFICANCE. C IER - RESULTING ERROR PARAMETER CODED AS FOLLOWS: C IER=0 - NO ERROR C IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAMETER N OR C BECAUSE SOME RADICAND IS NONPOSITIVE (MATRIX A C IS NOT POSITIVE DEFINITE, POSSIBLY DUE TO LOSS C OF SIGNIFICANCE) C IER=K - WARNING WHICH INDICATES LOSS OF SIGNIFICANCE. C THE RADICAND FORMED AT FACTORIZATION STEP K+1 C WAS STILL POSITIVE BUT NO LONGER GREATER THAN C ABS(EPS*A(K+1,K+1)). C C REMARKS C THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE C STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS. IN C THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGULAR MATRIX C IS STORED COLUMNWISE TOO. C THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL C CALCULATED RADICANDS ARE POSITIVE. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED. C MFSD C C METHOD C SOLUTION IS DONE USING THE FACTORIZATION BY SUBROUTINE MFSD. C DIMENSION A(*) DOUBLE PRECISION DIN,WORK C C FACTORIZE GIVEN MATRIX BY MEANS OF SUBROUTINE MFSD C C A=TRANSPOSE(T)*T C CALL MFSD(A,N,EPS,IER) IF(IER)9,1,1 C C INVERT UPPER TRIANGULAR MATRIX T C PREPARE INVERSION-LOOP C 1 IPIV=N*(N+1)/2 IND=IPIV C C INITIALIZE INVERSION-LOOP C DO 6 I=1,N DIN=1.D0/DBLE(A(IPIV)) A(IPIV)=DIN MIN=N KEND=I-1 LANF=N-KEND IF(KEND)5,5,2 2 J=IND C C INITIALIZE ROW-LOOP C DO 4 K=1,KEND WORK=0.D0 MIN=MIN-1 LHOR=IPIV LVER=J C C START INNER LOOP C DO 3 L=LANF,MIN LVER=LVER+1 LHOR=LHOR+L 3 WORK=WORK+DBLE(A(LVER)*A(LHOR)) C C END OF INNER LOOP C A(J)=-WORK*DIN 4 J=J-MIN C C END OF ROW-LOOP C 5 IPIV=IPIV-MIN 6 IND=IND-1 C C END OF INVERSION LOOP C C CALCULATE INVERSE(A) BY MEANS OF INVERSE(T) C INVERSE(A)=INVERSE(T)*TRANSPOSE(INVERSE(T)) C INITIALIZE MULTIPLICATION LOOP C DO 8 I=1,N IPIV=IPIV+I J=IPIV C C INITIALIZE ROW-LOOP C DO 8 K=I,N WORK=0.D0 LHOR=J C C START INNER LOOP C DO 7 L=K,N LVER=LHOR+K-I WORK=WORK+DBLE(A(LHOR)*A(LVER)) 7 LHOR=LHOR+L C C END OF INNER LOOP C A(J)=WORK 8 J=J+K C C END OF ROW- AND MULTIPLICATION-LOOP C 9 RETURN END C---------------------------------------- SUBROUTINE MFSD(A,N,EPS,IER) C---------------------------------------- C C PURPOSE C FACTOR A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX C C DESCRIPTION OF PARAMETERS C A - UPPER TRIANGULAR PART OF THE GIVEN SYMMETRIC POSITIVE C DEFINITE N BY N COEFFICIENT MATRIX. ON RETURN A C CONTAINS THE RESULTANT UPPER TRIANGULAR MATRIX. C N - THE NUMBER OF ROW (COLUMNS) IN GIVEN MATRIX. C EPS - AN INPUT CONSTANT WHICH IS USED AS RELATIVE TOLERANCE C FOR TEST ON LOSS OF SIGNIFICANCE. C IER - RESULTING ERROR PARAMETER CODED AS FOLLOWS: C IER=0 - NO ERROR C IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAMETER N OR C BECAUSE SOME RADICAND IS NONPOSITIVE (MATRIX A C IS NOT POSITIVE DEFINITE, POSSIBLY DUE TO LOSS C OF SIGNIFICANCE) C IER=K - WARNING WHICH INDICATES LOSS OF SIGNIFICANCE. C THE RADICAND FORMED AT FACTORIZATION STEP K+1 C WAS STILL POSITIVE BUT NO LONGER GREATER THAN C ABS(EPS*A(K+1,K+1)). C C REMARKS C THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE C STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS. IN C THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGULAR MATRIX C IS STORED COLUMNWISE TOO. C THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL C CALCULATED RADICANDS ARE POSITIVE. C THE PRODUCT OF RETURNED DIAGONAL TERMS IS EQUAL TO THE SQUARE C ROOT OF THE DETERMINANT OF THE GIVEN MATRIX. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C SOLUTION IS DONE USING THE SQUARE-ROOT METHOD OF CHOLESKY. C THE GIVEN MATRIX IS REPRESENTED AS THE PRODUCT OF 2 TRIANGULAR C MATRICES, WHERE THE LEFT HAND FACTOR IS THE TRANSPOSE OF THE C RETURNED RIGHT HAND FACTOR. C DIMENSION A(*) DOUBLE PRECISION DPIV,DSUM C C TEST ON WRONG INPUT PARAMETER N C IF(N-1)12,1,1 1 IER=0 C C INITIALIZE DIAGONAL-LOOP C KPIV=0 DO 11 K=1,N KPIV=KPIV+K IND=KPIV LEND=K-1 C C CALCULATE TOLERANCE C TOL=ABS(EPS*A(KPIV)) C C START FACTORIZATION-LOOP OVER K-TH ROW C DO 11 I=K,N DSUM=0.D0 IF(LEND)2,4,2 C C START INNER LOOP C 2 DO 3 L=1,LEND LANF=KPIV-L LIND=IND-L 3 DSUM=DSUM+DBLE(A(LANF)*A(LIND)) C C END OF INNER LOOP C C TRANSFORM ELEMENT A(IND) C 4 DSUM=DBLE(A(IND))-DSUM IF(I-K)10,5,10 C C TEST FOR NEGATIVE PIVOT ELEMENT AND FOR LOSS OF SIGNIFICANCE. C 5 IF(SNGL(DSUM)-TOL)6,6,9 6 IF(DSUM)12,12,7 7 IF(IER)8,8,9 8 IER=K-1 C C COMPUTE PIVOT ELEMENT C 9 DPIV=DSQRT(DSUM) A(KPIV)=DPIV DPIV=1.D0/DPIV GO TO 11 C C CALCULATE TERMS IN ROW C 10 A(IND)=DSUM*DPIV 11 IND=IND+I C C END OF DIAGONAL-LOOP C RETURN 12 IER=-1 RETURN END C------------------------------------------- SUBROUTINE MDCH(CS,DF,P,IER) C------------------------------------------- C C PURPOSE C CHI-SQUARED PROBABILITY DISTRIBUTION FUNCTION C C DESCRIPTION OF PARAMETERS C CS - INPUT VALUE FOR WHICH THE PROBABILITY IS C COMPUTED. CS MUST BE GREATER THAN OR EQUAL C TO ZERO. C DF - INPUT NUMBER OF DEGREES OF FREEDOM OF THE C CHI-SQUARED DISTRIBUTION. DF MUST BE GREATER THAN C OR EQUAL TO .5 AND LESS THAN OR EQUAL TO 200,000. C P - OUTPUT PROBABILITY THAT A RANDOM VARIABLE WHICH C FOLLOWS THE CHI-SQUARED DISTRIBUTION WITH DF C DEGREES OF FREEDOM IS LESS THAN OR EQUAL TO CS. C IER - OUTPUT ERROR CODE C =129 INDICATES THAT CS OR DF WAS SPECIFIED C INCORRECTLY. C =34 INDICATES THAT THE NORMAL PDF WOULD HAVE C PRODUCED AN UNDERFLOW C C REMARKS C FOR DEGREES OF FREEDOM GREATER THAN MAXDF=100, C THE NORMAL APPROXIMATION IS USED. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NDTR, MGAMAD C C METHOD C SEE EQUATION 26.4.14, 6.5.32, 6.5.29 , HANDBOOK OF C MATHEMATICAL FUNCTIONS, ABRAMOWITZ AND STEGUN, C DOVER PUBLICATIONS, INC., NEW YORK. C REAL CS,DF,P REAL PT2 DOUBLE PRECISION A,Z,DGAM,EPS,W,W1,B,Z1,HALF,ONE,THRTEN,THRD DOUBLE PRECISION MGAMAD REAL X DATA EPS/1.0D-6/,HALF/5.D-1/,THRTEN/13.D0/,ONE/1.D0/ DATA PT2/.2222222E0/ DATA THRD/.3333333333333333D0/ FUNC(W,A,Z)=W*DEXP(A*DLOG(Z)-Z) C C FIRST EXECUTABLE STATEMENT C TEST FOR INVALID INPUT VALUES C IF (DF .GE. .5 .AND. DF .LE. 2.E5 .AND. CS .GE. 0.0) GO TO 5 IER=129 P=-1 GO TO 9000 5 IER=0 C C SET P=0. IF CS IS LESS THAN OR EQUAL TO 10.**(-12) C IF (CS .GT. 1.E-12) GO TO 15 10 P=0.0 GO TO 9005 15 IF (DF .LE. 100.) GO TO 20 C C USE NORMAL DISTRIBUTION APPROXIMATION FOR LARGE C DEGREES OF FREEDOM C IF (CS .LE. 2.0) GO TO 10 X=((CS/DF)**THRD-(ONE-PT2/DF))/SQRT(PT2/DF) IF (X .GT. 5.0) GO TO 50 IF (X .LT. -18.8055) GO TO 55 CALL NDTR(X,P,DEN) GO TO 9005 C C INITIALIZATION FOR CALCULATION USING INCOMPLETE C GAMMA FUNCTION C 20 IF (CS .GT. 200.) GO TO 50 A=HALF*DF Z=HALF*CS DGAM=MGAMAD(A) W=DMAX1(HALF*A,THRTEN) IF (Z .GE. W) GO TO 35 IF (DF .GT. 25. .AND. CS .LT. 2.) GO TO 10 C C CALCULATE USING EQUATION NO. 6.5.29 C W=ONE/(DGAM*A) W1=W DO 25 I=1,50 B=I W1=W1*Z/(A+B) IF (W1 .LE. EPS*W) GO TO 30 W=W+W1 25 CONTINUE 30 P=FUNC(W,A,Z) GO TO 9005 C C CALCULATE USING EQUATION NO. 6.5.32 C 35 Z1=ONE/Z B=A-ONE W1=B*Z1 W=ONE+W1 DO 40 I=2,50 B=B-ONE W1=W1*B*Z1 IF (W1. LE. EPS*W) GO TO 45 W=W+W1 40 CONTINUE 45 W=Z1*FUNC(W,A,Z) P=ONE-W/DGAM GO TO 9005 50 P=1.0 GO TO 9005 55 P=0.0 IER=34 9000 CONTINUE 9005 RETURN END C----------------------------------------------- DOUBLE PRECISION FUNCTION MGAMAD(X) C----------------------------------------------- C C PURPOSE C EVALUATE THE GAMMA FUNCTION OF A POSITIVE DOUBLE PRECISION C ARGUMENT C C DESCRIPTION OF PARAMETERS C X -- INPUT DOUBLE PRECISION ARGUMENT (X MUST BE A C POSITIVE REAL NUMBER) C C METHOD C SEE CODY, W. J. AND HILLSTROM, K. E. - CHEBYSHEV APPROXIMATIONS C FOR THE NATURAL LOGARITHM OF THE GAMMA FUNCTION, MATHEMATICS C OF COMPUTATION, 21(89) 1967, 198-208; ALSO SEE EQUATION 6.1.41, C HANDBOOK OF MATHEMATICAL FUNCTIONS, ABRAMOWITZ AND STEGUN, C DOVER PUBLICATIONS, INC., NEW YORK. C DOUBLE PRECISION P(6,3),Q(6,3),PSUM,QSUM,DLNGA,X DATA PI/3.14159265358979D0/ DATA P/-2.6094066054623D0,-4.1509018875434D01, 1 -9.5564117677317D01,-1.1811439967596D01, 2 2.8113744347038D01,3.5173589912443D0, 3 -2.16192292624703D02,-8.27790897809598D02, 4 1.82987822012009D02,7.06543700154966D02, 5 1.49903662709861D02,4.54827477723909D0, 6 -2.4273113085758D06,1.3860869828508D06, 7 1.8537773351564D06,-6.4279927530351D05, 8 -1.5515971577126D05,-4.2105209252847D03/ DATA Q/5.4587504274950D-01,1.6674969701154D01, 1 7.8556098036754D01,8.7289305773548D01, 2 2.3590762639739D01,1.0D0, 3 1.16412659461333D02,1.20459293663292D03, 4 1.85645035686087D03,7.05287069715149D02, 5 6.65573507467416D01,1.0D0, 6 -1.0542482321634D06,-2.4515705199457D06, 7 -8.6274186723037D05,-6.7771258633073D04, 8 -9.4136613234388D02,1.0D0/ I=1 C INITIALIZE FUNCTION VALUE MGAMAD=0.0D0 IF (X.LE.0) RETURN C C IF X>12, USE EQUATION 6.1.41 C IF (X.GT.12.0) GO TO 60 C C USE COEFFICIENTS P(J,I) AND Q(J,I) DEPENDING ON C VALUES OF X C IF ((X.GE.4.0).AND.(X.LE.12.0)) I=3 IF ((X.GE.1.5).AND.(X.LE.4.0)) I=2 IF (X.LE.0.5) X=X+1.0D0 C C CALCULATE LOG-GAMMA C PSUM=0.D0 QSUM=0.D0 DO 20 J=1,5 PSUM=(PSUM+P(6-J+1,I))*X QSUM=(QSUM+Q(6-J+1,I))*X 20 CONTINUE PSUM=PSUM+P(1,I) QSUM=QSUM+Q(1,I) IF (I.EQ.2) DLNGA=(X-2.0D0)*(PSUM/QSUM) IF (I.EQ.3) DLNGA=PSUM/QSUM IF ((I.EQ.1).AND.(X.LE.0.5)) DLNGA=-DLOG(X)+(PSUM/QSUM) IF ((I.EQ.1).AND.(X.GT.0.5)) DLNGA=(X-1.0D0)*PSUM/QSUM GO TO 70 C C IF X>12, CALCULATE LOG-GAMMA USING EQUATION NO. 6.1.41 C 60 DLNGA=(X-0.5D0)*DLOG(X)-X+0.5D0*DLOG(2.D0*PI)+1.D0/(12.D0*X)- 1 1.D0/(360.D0*X**3)+1.D0/(1260.D0*X**5) 70 MGAMAD=DEXP(DLNGA) RETURN END