C C :---------------------: C : : C : 'ROCFIT' PROGRAM : Jan 1994 C : : C :---------------------: C C C******************************************************************************* C * C A PROGRAM FOR MAXIMUM LIKELIHOOD ESTIMATION OF A BINORMAL * C ROC CURVE AND ITS ASSOCIATED PARAMETERS FROM A SET OF CATEGORICAL * C RATING-SCALE DATA. FOR INHERENTLY CONTINUOUS DATA, THE 'LABROC4' * C PROGRAM SHOULD BE USED INSTEAD. * C * C THE MATHEMATICAL BASIS OF THE CURVE-FITTING ALGORITHM IS DESCRIBED * C IN REFERENCES (1) AND (2) BELOW. THIS PROGRAM IS A MODIFIED VERSION * C OF THE PROGRAM 'RSCORE II' BY: * C * C DONALD D. DORFMAN * C DEPARTMENT OF PSYCHOLOGY * C THE UNIVERSITY OF IOWA * C IOWA CITY, IOWA 52242 * C * C AS LISTED IN REFERENCE (3) BELOW. MODIFICATIONS TO 'RSCORE II' * C IMPLEMENTED IN 'ROCFIT' WERE DESIGNED BY: * C * C CHARLES E. METZ, JONG-HER SHEN, PU-LAN WANG, * C 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 * C AND WERE WRITTEN BY JONG-HER SHEN, PU-LAN WANG, AND HELEN B. KRONMAN. * C THESE MODIFICATIONS INCLUDE: * C (A). PROVIDING FOR INPUT OF A FREE-TEXT LINE TO DESCRIBE THE * C DATA SET IN THE OUTPUT HEADING. * C (B). ALLOWING THE USE OF FREE FORMAT FOR DATA INPUT. * C (C). ALLOWING EITHER HIGHER OR LOWER CATEGORY NUMBERS TO * C INDICATE GREATER CONFIDENCE THAT A SIGNAL (E.G., 'ABNORMALITY) * C IS PRESENT. * C (D). CHECKING CONSISTENCY OF INPUT DATA. * C (E). GIVING VARIABLES MORE MNEMONIC NAMES. * C (F). RENUMBERING STATEMENTS IN ORDER OF APPEARANCE. * C (G). TESTING THE RATING-SCALE DATA SET FOR 'DEGENERACY' * C BEFORE ATTEMPTING THE ITERATIVE CALCULATIONS. * C (H). CALCULATION OF 'DIFFSN' AND 'DIFFN' BETWEEN STATEMENTS * C 1060 AND 1065 IN SUBROUTINE 'TERMS', AND USE OF THESE IN * C SUBSEQUENT COMPUTATIONS. * C (I). MODIFICATION OF THE ITERATIVE PROCEDURE TO ENSURE THAT * C TWO CONDITIONS ARE SATISFIED AT EACH STEP: (1) THE MATRIX * C OF SECOND PARTIAL DERIVATIVES MUST NOT BECOME NUMERICALLY * C SINGULAR, AND (2) THE LOG-LIKELIHOOD FUNCTION MUST NOT * C DECREASE BY MORE THAN 10**(-5) OF ITS PREVIOUS ABSOLUTE * C VALUE. IF EITHER OF THESE CONDITIONS IS VIOLATED AT SOME * C STEP, THE ADJUSTMENTS OF ALL PARAMETER ESTIMATES ARE * C REDUCED SUCCESSIVELY BY 50% UNTIL THE CONDITIONS ARE * C SATISFIED. ITERATION THEN PROCEEDS TO THE NEXT STEP. * C (J). MODIFICATION OF THE STOPPING CRITERION TO REFLECT THE * C NUMBER OF CATEGORIES EMPLOYED IN THE DATA (I.E., TO * C REFLECT THE NUMBER OF ADJUSTABLE PARAMETERS). * C (K). CALCULATION AND OUTPUT OF A 'CORRELATION MATRIX' FOR * C THE PARAMETER ESTIMATES IN ADDITION TO THE 'VARIANCE- * C COVARIANCE MATRIX'. * C (L). CALCULATION AND OUTPUT OF 'TRUE POSITIVE FRACTIONS' ON * C THE FITTED ROC CURVE AT A VARIETY OF 'FALSE POSITIVE * C FRACTIONS', TO ALLOW THE FITTED ROC CURVE TO BE PLOTTED * C EASILY ON LINEAR AXES. * C (M). CALCULATION AND OUTPUT OF THE ASYMMETRIC 95% CONFIDENCE * C INTERVAL FOR 'TRUE POSITIVE FRACTION' AT EACH 'FALSE * C POSITIVE FRACTION'. * C (N). CALCULATION AND OUTPUT OF EXPECTED OPERATING POINTS * C ESTIMATED ON THE FITTED ROC CURVE, TOGETHER WITH ASSYMETRIC * C 95% CONFIDENCE INTERVALS FOR THOSE POINTS ALONG THE * C FITTED CURVE. * C * C REFERENCES: * C (1). DORFMAN, D. D. AND ALF, E. - JOURNAL OF MATHEMATICAL * C PSYCHOLOGY, 6: 487-496, 1969. * C (2). GREY, D. R. AND MORGAN, B. J. T. - JOURNAL OF * C MATHEMATICAL PSYCHOLOGY, 9: 128-139, 1972. * C (3). SWETS, J. A. AND PICKETT, R. M. - EVALUATION OF * C DIAGNOSTIC SYSTEMS: METHODS FROM SIGNAL DETECTION * C THEORY. (NEW YORK: ACADEMIC PRESS), 1982. * C * C * C INQUIRIES OR COMMENTS CONCERNING THIS PROGRAM SHOULD BE DIRECTED * C TO C. E. METZ AT THE ADDRESS ABOVE. * C * C * C ACKNOWLEDGEMENT: * C MODIFICATION OF THIS PROGRAM AT THE UNIVERSITY OF CHICAGO * C WAS SUPPORTED BY THE U.S. DEPARTMENT OF ENERGY UNDER * C CONTRACTS DE-AC02-80EV10359 AND DE-AC02-82ER60033 AND * C UNDER GRANT DE-FG02-86ER60418. * C * C * C******************************************************************************* C******************************************************************************* C * C REQUIRED INPUT AND DESCRIPTION OF OUTPUT: * C * C------------------------------------------------------------------------------* C * C INPUT --- * C * C LINE #1 : DESCRIPTION OF DATA SET (UP TO 80 CHARACTERS). * C LINE #2 : TOTAL NUMBER OF CATEGORIES, AND THE CATEGORY REPRESENTING * C THE STRONGEST EVIDENCE OF POSITIVITY (E.G., PRESENCE OF * C ABNORMALITY). * C LINE #3 : TOTAL NUMBER OF ACTUALLY NEGATIVE (E.G., 'ACTUALLY * C NORMAL') CASES, AND NUMBER OF RESPONSES IN EACH * C CATEGORY, IN THE ORDER OF THE CATEGORY NUMBERS. * C LINE #4 : TOTAL NUMBER OF ACTUALLY POSITIVE (E.G., 'ACTUALLY * C ABNORMAL') CASES, AND NUMBER OF RESPONSES IN EACH * C CATEGORY, IN THE ORDER OF THE CATEGORY NUMBERS. * C * C FORMAT : FREE INPUT FORMAT, I.E., DATA CAN START IN ANY * C COLUMN BETWEEN COLUMN 1 AND COLUMN 80. IF THERE IS MORE * C THAN ONE INPUT VALUE PER LINE (LINES #2, 3, AND 4), * C THOSE INPUT VALUES MUST BE SEPARATED BY AT LEAST * C ONE BLANK SPACE. * C * C INPUT EXAMPLE: * C ---------------------------------------------------------------- * C TYPICAL EXAMPLE OF 5-CATEGORY DATA * C 5 5 * C 60 30 19 8 2 1 * C 50 5 6 5 12 22 * C ---------------------------------------------------------------- * C * C REMARKS CONCERNING THE INPUT: * C 1. THE TOTAL NUMBER OF CATEGORIES (THE FIRST VALUE INPUT * C ON LINE #2) CAN BE 3 TO 11. * C 2. THE SECOND INPUT VALUE ON LINE #2 INDICATES WHETHER THE * C HIGHEST OR LOWEST NUMBERED CATEGORY REPRESENTS THE STRONGEST * C EVIDENCE OF POSITIVITY (E.G., PRESENCE OF ABNORMALITY). * C IF HIGHER-NUMBERED CATEGORIES REPRESENT STRONGER EVIDENCE * C OF POSITIVITY, THEN THIS VALUE SHOULD BE EQUAL TO THE TOTAL * C NUMBER OF CATEGORIES INPUT FIRST ON THIS LINE. ALTERNATIVELY, * C IF LOWER-NUMBERED CATEGORIES REPRESENT STRONGER EVIDENCE * C OF POSITIVITY, THEN THIS VALUE SHOULD BE '1'. * C 3. THIS PROGRAM ALLOWS THE USE OF SEVERAL DATA SETS AS INPUT; * C SIMPLY ENTER THE 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 OUTPUT --- * C * C 1. HEADING, WITH DESCRIPTION OF DATA SET, NUMBER OF CATEGORIES, * C 'DIRECTION' OF RATING SCALE, AND NUMBERS OF CASES. * C 2. CONFIDENCE-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 * C ACTUALLY POSITIVE DISTRIBUTION, WHEREAS THE PARAMETER 'A' * C CAN BE INTERPRETED AS THE ABSOLUTE DIFFERENCE BETWEEN THE * C MEANS OF TWO DISTRIBUTIONS, RELATIVE TO THE STANDARD * C DEVIATION OF THE ACTUALLY POSITIVE DISTRIBUTION. THE Z(K) * C ARE THE CUT-OFFS ON THE DECISION-VARIABLE AXIS THAT DEFINE * C THE CATEGORY BOUNDARIES; THEY ARE MEASURED FROM THE MEAN OF * C THE ACTUALLY NEGATIVE DISTRIBUTION TOWARD THE MEAN OF THE * C ACTUALLY POSITIVE DISTRIBUTION AND ARE EXPRESSED IN UNITS OF * C THE STANDARD DEVIATION OF THE EFFECTIVE NORMAL DISTRIBUTION * C ACTUALLY NEGATIVE CASES. * C 5. NUMBER OF ITERATIONS EMPLOYED. * C 6. FINAL PARAMETER ESTIMATES, WITH LOG LIKELIHOOD OF 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 LINEAR * C AXES (OFTEN DENOTED BY '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 THE ASYMMETRIC 95% CONFIDENCE INTERVAL * C THAT IS ESTIMATED BY ASSUMING THAT ERRORS IN THE ROC PARAMETER * C ESTIMATES 'A' AND 'B' ARE JOINTLY NORMAL. * C 10. FINAL ESTIMATES OF THE EXPECTED OPERATING POINTS ON THE FITTED * C ROC CURVE, WITH LOWER AND UPPER BOUNDS (ALONG THE FITTED CURVE) * C OF THE ASYMMETRIC 95% CONFIDENCE INTERVALS FOR THOSE OPERATING- * C POINT ESTIMATES, CALCULATED ASSUMING THAT ERRORS IN THE Z(K) * C ESTIMATES ARE NORMALLY DISTRIBUTED. ESTIMATES OF THE EXPECTED * C OPERATING POINTS ON THE FITTED CURVE, ON NORMAL-DEVIATE AXES, * C HAVE ABSCISSA VALUES = -Z(K) AND ORDINATE VALUES = A-B*Z(K). * C * C REMARKS CONCERNING THE OUTPUT: * C 1. THIS PROGRAM WILL CHECK DEGENERACY OF THE INPUT DATA. * C DEGENERATE DATA SETS ARE THOSE FOR WHICH AN EXACT FIT * C TO THE OBSERVED OPERATING POINTS IS PROVIDED BY A * C BINORMAL ROC CURVE AND ASSOCIATED CUT-OFFS WITH ONE * C OR MORE INFINITE PARAMETER VALUES, SO THAT THE * C ITERATIVE CALCULATION CANNOT CONVERGE. IF THE DATA * C SET IS DEGENERATE, THIS PROGRAM WILL OUTPUT A MESSAGE * C DESCRIBING THE KIND OF DEGENERACY FOUND AND THE CORRESPONDING * C EXACT-FIT ROC, AND THEN ABORT EXECUTION OF THAT DATA SET. * C 2. THE 'CHI-SQUARE' GOODNESS-OF-FIT MEASURE IS NOT OUTPUT IF ANY * C EXPECTED CELL FREQUENCY IS LESS THAN 5, OR IF ONLY 3 * C RATING CATEGORIES WERE EMPLOYED. * C 3. IF THE OUTPUT VALUE OF THE AREA INDEX ('A SUB Z') * C IS SIGNIFICANTLY LESS THAN 0.5, THE 'DIRECTION' OF THE * C CONFIDENCE-RATING SCALE MAY HAVE BEEN REVERSED: SEE * C THE SECOND 'REMARK CONCERNING THE INPUT' ABOVE TO ENSURE THAT * C THE SECOND VALUE INPUT ON LINE #2 WAS CORRECT. * C * C * C******************************************************************************* INTEGER IARY(11) 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) IEND=0 IERROR=0 C C READ IN C 1000 CALL READIN(IERROR,IEND) IF(IEND.EQ.1 .OR. IERROR.GT.0)GO TO 8500 C C PRINT HEADINGS C CALL OUTINI C C CHECK CONSISTENCY OF INPUT DATA C CALL CHKIN(IERROR) IF(IERROR.EQ.1)GO TO 1000 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) 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.GT.100)GO TO 1000 C C PRINT FINAL RESULT C CALL OUTRSL(LL,XP) GO TO 1000 8500 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 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 Changed 6/1/93: Holle. Const. needs \\ for \on C-prepro. systems 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 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 WRITE(6,1030) 1030 FORMAT(1H1,29X,'R O C F I T (JUNE 1993 VERSION) :'// 1 18X,'M A X I M U M L I K E L I H O O D', 2 ' E S T I M A T I O N',/,24X,'O F A B I N O R ', 3 'M A L R O C C U R V E'/30X,'F R O M R A T I N ', 4 'G D A T A'//) WRITE(6,1040)(IDENT(I),I=1,80) 1040 FORMAT(20X,'DATA DESCRIPTION: ',80A1,/) WRITE(6,1041)KAT,KSN 1041 FORMAT(20X,'DATA COLLECTED IN ',I2,' CATEGORIES'/ 1 20X,'WITH CATEGORY ',I2,' REPRESENTING STRONGEST EVIDENCE', 2 ' OF POSITIVITY (E.G., THAT ABNORMALITY IS PRESENT).'/) WRITE(6,1050)EMN,EMS 1050 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)= ',12(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 tnet=1.0 alpha=1.0 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)/tnet) IF (alpha.le.1.0e-5.or.FUNDIF.LE.1.0E-4) 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(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 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)= ',12(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 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 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 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(12),TPF(12) INTEGER IARY(11) 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(12),TPF(12) INTEGER IARY(11) 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