C C :------------------------: C : : C : 'CLABROC' PROGRAM : JAN 94 C : : C :------------------------: C C C******************************************************************************* C * C * C THIS PROGRAM IS APPLICABLE TO CORRELATED, CONTINUOUSLY-DISTRIBUTED * C TEST-RESULT DATA. FOR CORRELATED BUT INHERENTLY CATEGORICAL DATA, * C THE STAND-ALONE PROGRAM ENTITLED 'CORROC2' (FROM WHICH 'CLABROC' IS * C ADAPTED) SHOULD BE USED INSTEAD. * C * C THE PURPOSES OF THIS PROGRAM ARE: * C * C (I) TO CALCULATE MAXIMUM LIKELIHOOD ESTIMATES OF THE PARAMETERS * C OF A MODEL FOR CORRELATED, CONTINUOUSLY-DISTRIBUTED * C TEST-RESULT DATA AND THE BINORMAL ROC CURVES IMPLIED BY THOSE * C DATA. * C AND: * C * C (II) TO CALCULATE THE STATISTICAL SIGNIFICANCE OF THE DIFFERENCE * C BETWEEN THE TWO ESTIMATED ROC CURVES BY ANY ONE OF THREE * C DISTINCT STATISTICAL TESTS FOR CORRELATED ROC DATA: * C * C (1) A BIVARIATE CHI-SQUARE TEST OF THE SIMULTANEOUS * C DIFFERENCES BETWEEN THE 'A' PARAMETERS AND BETWEEN * C THE 'B' PARAMETERS OF THE TWO ROC CURVES. (NULL * C HYPOTHESIS: THE DATA SETS AROSE FROM THE SAME * C BINORMAL ROC CURVE). * C * C (2) A UNIVARIATE Z-SCORE TEST OF THE DIFFERENCE BETWEEN * C THE AREAS UNDER THE TWO ROC CURVES. (NULL * C HYPOTHESIS: THE DATA SETS AROSE FROM BINORMAL ROC * C CURVES WITH EQUAL AREAS BENEATH THEM); OR * C * C (3) A UNIVARIATE Z-SCORE TEST OF THE DIFFERENCE BETWEEN * C THE TRUE-POSITIVE-FRACTIONS ON THE TWO ROC CURVES * C AT A SELECTED FALSE-POSITIVE-FRACTION. (NULL * C HYPOTHESIS: THE DATA SETS AROSE FROM BINORMAL ROC * C CURVES HAVING THE SAME TPF AT THE SELECTED FPF). * C * C THE REQUIRED DATA ARE TWO SETS OF PAIRED, CONTINUOUSLY-DISTRIBUTED * C TEST RESULTS -- ONE SET FROM ACTUALLY NEGATIVE (E.G., 'NORMAL') CASES * C AND THE OTHER SET FROM ACTUALLY POSITIVE (E.G., 'ABNORMAL') CASES. * C THE TWO TEST RESULTS ASSOCIATED WITH EACH CASE COULD ARISE, * C FOR EXAMPLE, FROM TWO DIFFERENT DIAGNOSTIC TESTS PERFORMED ON THE * C SAME PATIENT. * C * C *** NOTE *** THE CASES USED AS INPUT MUST BE DISTINCT, BECAUSE * C THIS PROGRAM ASSUMES THAT EACH PAIR OF CONTINUOUSLY- * C DISTRIBUTED TEST RESULTS WAS SAMPLED INDEPENDENTLY * C FROM ONE OF TWO BIVARIATE DISTRIBUTIONS (CORRESPONDING TO * C 'ACTUALLY NEGATIVE' AND 'ACTUALLY POSITIVE' CASES). * C HENCE, THE INPUT DATA MUST NOT INCLUDE MORE THAN ONE * C PAIR OF TEST RESULTS FROM EACH CASE (E.G., PATIENT). * C IF MULTIPLE TEST-RESULT PAIRS FROM EACH CASE ARE POOLED * C IN THE INPUT, THE PROGRAM WILL OVERESTIMATE THE * C STATISTICAL SIGNIFICANCE OF ANY APPARENT DIFFERENCE * C BETWEEN THE ROC CURVES, THEREBY INVALIDATING THE * C STATISTICAL TEST. WHEN MULTIPLE DIAGNOSTIC TEST-RESULT * C PAIRS ARE AVAILABLE FOR EACH CASE (FOR EXAMPLE, FROM * C REPLICATION IN EACH OF THE TWO CONDITIONS), THE GENERAL * C APPROACH DESCRIBED BY SWETS AND PICKETT ('EVALUATION * C OF DIAGNOSTIC SYSTEMS: METHODS FROM SIGNAL DETECTION * C THEORY.' ACADEMIC PRESS, NEW YORK, 1982, CHAPTER 4) * C SHOULD BE EMPLOYED INSTEAD OF THIS PROGRAM. * C * C * C IN A PRELIMINARY STEP, CLABROC AUTOMATICALLY CATEGORIZES THE * C CONTINUOUSLY-DISTRIBUTED INPUT DATA IN AN ARBITRARY BUT EFFECTIVE * C WAY. NEXT, THE MARGINAL CATEGORICAL DATA SETS THUS CREATED * C ARE ANALYZED INDEPENDENTLY BY A MODIFIED DORFMAN PROGRAM * C (J. MATH. PSYCH., VOL.6, 1969, PP.487-496, AND VOL.9, 1972, * C PP. 128-139) TO OBTAIN MAXIMUM LIKELIHOOD ESTIMATES OF THE * C 'A', 'B', AND CATEGORY-BOUNDARY PARAMETERS SEPARATELY FOR EACH * C ROC CURVE. TOGETHER WITH PRODUCT-MOMENT CORRELATION COEFFICIENTS * C CALCULATED DIRECTLY FROM THE JOINT CATEGORIZED DATA, THESE RESULTS * C ARE THEN USED AS INITIAL ESTIMATES IN A LATER PART OF THE ALGORITHM * C TO COMPUTE (BY THE METHOD OF SCORING) MAXIMUM LIKELIHOOD ESTIMATES * C OF THE PARAMETERS OF A BIVARIATE-BINORMAL MODEL THAT IS ASSUMED * C TO UNDERLIE THE CORRELATED TEST RESULTS. * C * C THE BIVARIATE-BINORMAL MODEL USED BY THIS PROGRAM AND THE STATISTICAL * C TESTS PERFORMED ARE DESCRIBED IN THE PAPER 'A NEW APPROACH FOR * C TESTING THE SIGNIFICANCE OF DIFFERENCES BETWEEN ROC CURVES MEASURED * C FROM CORRELATED DATA' BY METZ, WANG AND KRONMAN IN 'INFORMATION * C PROCESSING IN MEDICAL IMAGING' (F. DECONINCK, ED.). THE HAGUE: * C MARTINUS NIJHOFF, 1984, PP. 432-445. THE MODEL IS BASED ON AN * C ASSUMPTION THAT THE TEST RESULTS -- OR SOME UNKNOWN MONOTONIC * C TRANSFORMATIONS THEREOF -- FOLLOW BIVARIATE NORMAL DISTRIBUTIONS; * C THIS ASSUMPTION MAY BE VALID EVEN WHEN THE DISTRIBUTIONS OF THE RAW * C DATA ARE SKEWED AND/OR MULTI-MODAL. * C * C PORTIONS OF THIS PROGRAM (SUBROUTINES 'ESTM1', 'DORALF', AND * C THE SUBROUTINES THAT THEY CALL) ASSOCIATED WITH CALCULATING * C MAXIMUM LIKELIHOOD ESTIMATES OF THE ROC PARAMETERS OF THE * C MARGINAL DATA SETS -- WHICH ARE USED AS INITIAL ESTIMATES FOR * C THE BIVARIATE MODEL CALCULATION -- WHERE TAKEN (WITH SOME * C MODIFICATION) FROM THE PROGRAM 'RSCORE II' BY DONALD D. DORFMAN, * C DEPARTMENT OF PSYCHOLOGY, THE UNIVERSITY OF IOWA. * C * C * C ALGORITHM DESIGNED BY CHARLES E. METZ, JONG-HER SHEN, * C HELEN B. KRONMAN, AND PU-LAN WANG. * C PROGRAM WRITTEN BY JONG-HER SHEN, HELEN B. KRONMAN, * C AND PU-LAN WANG. * C * C DEPARTMENT OF RADIOLOGY AND THE FRANKLIN MCLEAN * C MEMORIAL RESEARCH INSTITUTE * C THE UNIVERSITY OF CHICAGO, CHICAGO, ILLINOIS 60637 * C * C INQUIRIES OR COMMENTS CONCERING THIS PROGRAM SHOULD BE DIRECTED * C TO C.E. METZ AT THE ABOVE ADDRESS. * C * C THIS WORK WAS SUPPORTED BY THE US 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******************************************************************************* C * C REQUIRED INPUT AND DESCRIPTION OF OUTPUT: * C * C------------------------------------------------------------------------------* C * C * C INPUT --- * C * C 1. FIRST A CODE TO INDICATE WHETHER THE ENTIRE INPUT SHOULD BE * C PRINTED; IF THE CODE BEGINS WITH 'Y' (E.G., 'YES'), THEN * C THE PROGRAM WILL ECHO-PRINT ALL OF THE INPUT DATA. * C 2. TWO LINES FOR EACH OF TWO 'CONDITIONS', X AND Y, WITH THE * C FOLLOWING INFORMATION: * C (NOTE: THE TWO 'CONDITIONS' (X AND Y) CAN REPRESENT DIFFERENT * C DIAGNOSTIC TESTS, TWO DIFFERENT QUANTITIES PROVIDED * C BY A SINGLE DIAGNOSTIC TEST, INITIAL AND SUBSEQUENT * C RESULTS FROM A SINGLE TEST, ETC.) * C LINE #2 AND #4 : A FREE-TEXT DESCRIPTION OF THE CONDITION * C (UP TO 80 CHARACTERS). * C LINE #3 AND #5 : AN ALPHABETIC CODE WORD ('SMALL' OR * C 'LARGE') INDICATING THAT SMALLER OR * C LARGER TEST RESULTS, RESPECTIVELY, ARE * C ASSOCIATED WITH STRONGER EVIDENCE OF * C POSITIVITY (E.G., 'ABNORMALITY'). * C (NOTE: THE COMPUTER READS ONLY THE FIRST * C CHARACTER OF THIS CODE WORD.) * C FORMAT : FREE INPUT FORMAT. INPUT DATA CAN START IN ANY * C COLUMN BETWEEN COLUMN 1 AND COLUMN 80. * C 3. TEST-RESULT PAIRS FOR ACTUALLY NEGATIVE CASES AND THEN FOR * C ACTUALLY POSITIVE CASES. * C FORMAT : ON A LINE FOR EACH CASE, ENTER THE TEST RESULT FOR * C THE X-CONDITION, ONE OR MORE BLANK SPACES, THE TEST * C RESULT FOR THE Y-CONDITION, ONE OR MORE BLANK SPACES, * C AND FINALLY ANY DESIRED FREE-TEXT DESCRIPTION OF THE * C CASE (THE DESCRIPTION IS OPTIONAL). ENTER ALL * C OF THE ACTUALLY NEGATIVE (E.G., 'NORMAL') CASES FIRST * C AND THEN ALL OF THE ACTUALLY POSITIVE (E.G., 'ABNORMAL') * C CASES. EACH OF THESE TWO INPUT SEQUENCES MUST BE * C TERMINATED BY A SEPARATE LINE WITH AN ASTERISK (*). * C (NOTE: THIS PROGRAM ACCEPTS UP TO 1000 ACTUALLY * C NEGATIVE TEST-RESULT PAIRS AND UP TO 1000 * C ACTUALLY POSITIVE TEST-RESULT PAIRS. FOR EACH * C TEST RESULT, THE PROGRAM ACCEPTS UP TO 6 DIGITS * C TO THE LEFT OF THE DECIMAL POINT AND UP TO 8 * C DIGITS TO THE RIGHT OF THE DECIMAL POINT.) * C 4. AN ALPHABETIC CODE WORD INDICATING THE STATISTICAL TEST DESIRED: * C FOR THE BIVARIATE TEST, ENTER: BIVARIATE * C FOR THE AREA TEST, ENTER: AREA * C FOR THE 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 TEST OR THE AREA TEST WAS REQUESTED. * C FORMAT : FREE INPUT FORMAT. * C * C * C INPUT EXAMPLE FOR TPF TEST: * C --------------------------------------------------------- * C YES * C SAMPLE INPUT DATA --- X * C LARGE * C SAMPLE INPUT DATA --- Y * C LARGE * C 0.12 -0.2 NORMAL CASE #1 * C 1.1 .29 NORMAL CASE #2 * C 0.99 1.2 NORMAL CASE #3 * C ... (AND SO ON) * C * * C 1.1 2.2 ABNORMAL CASE #1 * C 1.3 -0.5 ABNORMAL CASE #2 * C 2.3 1.4 ABNORMAL CASE #3 * C ... (AND SO ON) * C * * C TPF * C 0.10 * C --------------------------------------------------------- * C * C REMARKS CONCERNING THE INPUT: * C THIS PROGRAM ALLOWS THE USE OF SEVERAL DATA SETS AS INPUT; * C SIMPLY ENTER THE DATA SETS SEQUENTIALLY, BEGINNING WITH * C LINE #1. * C * C * C * C OUTPUT --- * C * C 1. ECHO-PRINT ALL INPUT DATA, IF REQUESTED. * C 2. HEADING, WITH A DESCRIPTION OF THE STATISTICAL TEST TO * C BE EMPLOYED AND THEN, FOR EACH OF THE TWO CONDITIONS, * C A DESCRIPTION OF THE CONDITION AND OF THE 'DIRECTION' OF * C THE TEST-RESULT SCALE FOR EACH CONDITION. * C 3. TOTAL NUMBER OF ACTUALLY NEGATIVE CASES AND TOTAL NUMBER OF * C ACTUALLY POSITIVE CASES. * C 4. ROC OPERATING POINTS OBTAINED DIRECTLY FROM THE MARGINAL * C DATA FROM THE TWO CONDITIONS, RESPECTIVELY, AFTER CATEGORIZATION * C BUT BEFORE CURVE FITTING. * C 5. THE NUMBER OF ITERATIONS EMPLOYED. * C 6. FINAL ESTIMATES OF THE BINORMAL ROC PARAMETERS FOR EACH * C CONDITION, AND OF THE EFFECTIVE INTER-CONDITION CORRELATIONS * C FOR ACTUALLY NEGATIVE CASES AND ACTUALLY POSITIVE CASES. FOR * C EACH CONDITION, THE PARAMETER 'A' REPRESENTS THE VERTICAL * C INTERCEPT AND 'B' REPRESENTS THE SLOPE OF THE FITTED ROC * C CURVE WHEN IT IS PLOTTED AS A STRAIGHT LINE ON 'NORMAL- * C DEVIATE' AXES. THE OUTPUT CORRELATION VALUES FOR ACTUALLY * C NEGATIVE AND ACTUALLY POSITIVE CASES REPRESENT, THE INTER- * C CONDITION CORRELATION COEFFICIENTS OF THE CORRESPONDING * C BIVARIATE-NORMAL DISTRIBUTIONS THAT WOULD BE OBTAINED BY * C APPROPRIATE MONOTONIC TRANSFORMATIONS OF THE TEST RESULTS FROM * C CONDITIONS 'X' AND 'Y'. * C 7. ESTIMATED STANDARD DEVIATIONS OF THE FINAL ESTIMATES OF * C THE ROC PARAMETERS AND OF THE INTER-CONDITION CORRELATION. * C 8. TPF VALUES FOR CONDITION X AND FOR CONDITION Y ON THE FITTED * C ROC CURVES, AT SELECTED FPF VALUES. * C 9. THE AREAS ('A SUB Z') UNDER THE FITTED ROC CURVES FOR THE TWO * C CONDITIONS, WITH ESTIMATED STANDARD DEVIATIONS AND CORRELATIONS. * C 10. ONE OF THREE TEST-STATISTIC VALUES (DEPENDING ON THE TEST * C SELECTED), WITH ASSOCIATED SIGNIFICANCE LEVEL OR LEVELS: * C -- THE COMPUTED 'BIVARIATE CHI-SQUARE TEST' STATISTIC * C VALUE WITH ITS CORRESPONDING P-LEVEL. * C -- THE COMPUTED 'AREA TEST' STATISTIC VALUE WITH * C CORRESPONDING TWO-TAILED AND ONE-TAILED P-LEVELS. * C -- THE COMPUTED 'TRUE-POSITIVE FRACTION TEST' STATISTIC * C VALUE AT THE SPECIFIED FALSE-POSITIVE FRACTION, * C WITH CORRESPONDING TWO-TAILED AND ONE-TAILED * C P-LEVELS. * C * C REMARKS CONCERNING THE OUTPUT: * C 1. BEFORE MAXIMUM LIKELIHOOD ESTIMATION OF THE ROC CURVES IS * C ATTEMPTED, THE CONTINUOUSLY-DISTRIBUTED INPUT DATA ARE * C AUTOMATICALLY SORTED INTO 'K' CATEGORIES, SUCH THAT THE * C FRACTION OF NORMAL CASES PLUS THE FRACTION OF ABNORMAL CASES * C IN EACH CATEGORY IS APPROXIMATLY CONSTANT. THIS CAUSES THE * C VERTICAL DISTANCE PLUS THE HORIZONTAL DISTANCE BETWEEN * C ADJACENT OPERATING POINTS TO BE APPROXIMATELY CONSTANT ON * C LINEAR AXES, THUS PROVIDING A ROUGHLY UNIFORM SPREAD OF * C OPERATING POINTS WITH POINTS SLIGHTLY CLOSER TOGETHER WHERE * C CURVATURE OF THE ROC IS GREATER. IF BOTH THE NUMBER OF * C ACTUALLY NEGATIVE CASES AND THE NUMBER OF ACTUALLY * C POSITIVE CASES IN THE INPUT ARE GREATER THAN OR EQUAL * C 100, THEN K=20; OTHERWISE, K=10. HOWEVER, DATA SETS CONTAINING * C LARGE NUMBERS OF TIED TEST-RESULT VALUES SOMETIMES DO NOT * C ALLOW THESE NUMBERS OF CATEGORIES TO BE ACHIEVED; IN SUCH * C SITUATIONS, THE PROGRAM PRODUCES AS MANY CATEGORIES AS * C POSSIBLE. * C 2. IT IS CRUCIALLY IMPORTANT TO NOTE THAT THIS PROGRAM * C IS NOT -- REPEAT, IS NOT -- APPROPRIATE FOR POOLED * C TEST-RESULT PAIRS (SEE NOTE ON FIRST PAGE). * C 3. THIS PROGRAM WILL CHECK BOTH MARGINAL DATA SETS FOR * C DEGENERACY. DEGENERATE DATA SETS ARE THOSE FOR WHICH * C AN EXACT FIT TO THE OBSERVED OPERATING POINTS IS PROVIDED * C BY A BINORMAL ROC CURVE AND ASSOCIATED CUT-OFFS WITH * C ONE OR MORE INFINITE PARAMETER VALUES, SO THAT THE * C ITERATIVE CALCULATIONS CANNOT CONVERGE. IF * C EITHER MARGINAL DATA SET IS DEGENERATE, THIS PROGRAM * C WILL OUTPUT A MESSAGE DESCRIBING THE KIND OF DEGENERACY * C FOUND AND THEN ABORT EXECUTION OF THE CORRESPONDING * C BIVARIATE DATA SET. IN GENERAL, DEGENERACY SHOULD BE FOUND * C ONLY IN VERY SMALL DATA SETS. * C 4. THE ROC CURVE THAT THIS PROGRAM ESTIMATES FOR EACH CONDITION * C FROM PAIRED RATING DATA MAY BE SLIGHTLY DIFFERENT FROM THE * C CURVE THAT THE 'ROCFIT' PROGRAM ESTIMATES FROM THE RATING * C DATA OF THE CONDITION ALONE. MORE GENERALLY, THE ROC CURVE * C THAT 'CLABROC' ESTIMATES FOR A GIVEN CONDITION CAN DEPEND TO * C SOME DEGREE UPON THE DATA WITH WHICH THE CONDITION IS PAIRED. * C FOR EXAMPLE, IF TRIPLETS OF RATINGS ARE OBTAINED FROM THREE * C CONDITIONS (X, Y, AND Z), THEN THE ROC THAT 'CLABROC' * C ESTIMATES FOR CONDITION X MAY DEPEND ON WHETHER THE RATINGS * C FROM CONDITION X ARE PAIRED WITH THOSE FROM CONDITION Y OR * C THOSE FROM CONDITION Z. IN OUR EXPERIENCE, THESE EFFECTS * C USUALLY ARE SMALL AND CAN BE IGNORED. IN SITUATIONS WHERE * C DIFFERENCES AMONG MULTIPLE ESTIMATES OF A PARTICULAR ROC * C ARE MORE SUBSTANTIAL, WE RECOMMEND THAT AN OVERALL ESTIMATE * C OF THE ROC BE OBTAINED BY AVERAGING THE INDIVIDUAL ESTIMATES * C OF THE 'A' PARAMETER AND THE INDIVIDUAL ESTIMATES OF THE 'B' * C PARAMETER OF THAT CURVE. THEN THE OVERALL ROC IS GIVEN BY * C THE EXPRESSIONS FPF=PHI(V) AND TPF=PHI(A+BV), WHERE 'V' TAKES * C ON VALUES FROM -3 TO +3 AND PHI(Z) REPRESENTS THE CUMULATIVE * C STANDARD-NORMAL DISTRIBUTION FUNCTION. * C * C * C******************************************************************************* COMMON/BLK1/VV(21,21),WW(21,21),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) COMMON/E/IEND,IERROR,IECHO COMMON/DECIML/MAXDX,MAXDY COMMON/OP/ZZX(20),ZZY(20) INTEGER TCODE DIMENSION ACTNI(1000),ACTSNI(1000),ACTNJ(1000),ACTSNJ(1000), 1 IDN(1000,80),IDSN(1000,80) IERROR=0 IEND=0 OPEN(UNIT=5,FILE='~/data.dat',STATUS='OLD') OPEN(UNIT=6,FILE='CLABROC.OUT',STATUS='NEW') C C READ IN ECHO PRINT REQUEST CODE C 10 CALL RDECHO IF(IEND.EQ.1)GO TO 199 C C READ DATA DESCRIPTION, CODE WORD AND TWO SEQUENCES OF C QUANTITATIVE DATA FOR ACTUALLY NEGATIVE CASES AND C ACTUALLY POSITIVE CASES C CALL READHL(KEYI,IDENTI) IF(IERROR.GT.0)GO TO 199 CALL READHL(KEYJ,IDENTJ) IF(IERROR.GT.0)GO TO 199 MAXDX=0 MAXDY=0 CALL READQ(MN,ACTNI,ACTNJ,IDN) IF(IERROR.GT.0)GO TO 199 IF(IECHO.EQ.1)CALL PPRINT(MN,ACTNI,ACTNJ,IDN) CALL READQ(MS,ACTSNI,ACTSNJ,IDSN) IF(IERROR.GT.0)GO TO 199 IF(IECHO.EQ.1)CALL PPRINT(MS,ACTSNI,ACTSNJ,IDSN) C C READ STATISTICAL REQUEST CODE WORD C CALL READT(TCODE,FP) C C INITIALIZE VARIABLES C EMN=MN EMS=MS C C CONVERT QUANTITATIVE DATA INTO CATEGORICAL DATA C CALL CATGRZ(1,KEYI,ACTNI,ACTSNI) CALL CATGRZ(2,KEYJ,ACTNJ,ACTSNJ) C C SETUP RATING DATA MATRIX C CALL SETMTX(MN,ACTNI,ACTNJ,VV) CALL SETMTX(MS,ACTSNI,ACTSNJ,WW) C C RUN THROUGH "CORROC2" C NNI=KATI+1 NNJ=KATJ+1 100 CALL CORROC2 GO TO 10 199 STOP END C-------------------------------------------- SUBROUTINE RDECHO C-------------------------------------------- C C READ ECHO PRINT CODE C COMMON/E/IEND,IERROR,IECHO COMMON/PASS/LSTRING(80),LINE(80),LENGTH DATA NECHO1/1HN/,NECHO2/1Hn/,KECHO1/1HY/,KECHO2/1Hy/ C READ(5,100,END=300)(LINE(I),I=1,80) 100 FORMAT(80A1) WRITE(6,110) 110 FORMAT(1H1,13X,'CLABROC(JUNE 1993 VERSION):'/ 1 2X,'MAXIMUM LIKELIHOOD ESTIMATION OF BINORMAL ROC CURVES'/ 2 2X,'FROM CORRELATED CONTINUOUSLY-DISTRIBUTED TEST RESULTS'/ 3 2X,' AND CALCULATION OF THE STATISTICAL SIGNIFICANCE OF '/ 4 2X,' ANY APPARENT DIFFERENCE BETWEEN THE CURVES.'/) CALL GETWRD(NOTGET) IF(NOTGET.EQ.0)GO TO 120 WRITE(6,115) 115 FORMAT(2X,'--- NO INPUT DATA FOR ECHO-PRINT REQUEST CODE ---'/ + 5X,'ASSIGN BY DEFAULT: ECHO-PRINT REQUEST CODE = YES') RETURN C 120 IF((LSTRING(1).EQ.NECHO1) .OR. (LSTRING(1).EQ.NECHO2) .OR. + (LSTRING(1).EQ.KECHO1) .OR. (LSTRING(1).EQ.KECHO2)) + GO TO 130 C C INVALID ECHO-PRINTED REQUEST CODE DETECTED C IECHO=1 WRITE(6,125)(LSTRING(I),I=1,LENGTH) 125 FORMAT(2X,80A1) WRITE(6,126) 126 FORMAT(1X,'--- INVALID ECHO-PRINTED REQUEST CODE ---'/ + 5X,'ASSIGN BY DEFAULT = YES'/ + 2X,'ECHO-PRINT OF THE INPUT DATA THAT WILL BE USED IN ', + 'THE ANALYSIS THAT FOLLOWS:') RETURN C C NO ECHO-PRINTED REQUEST C 130 IF((LSTRING(1).EQ.KECHO1) .OR. (LSTRING(1).EQ.KECHO2))GO TO 140 IECHO=0 WRITE(6,135) 135 FORMAT(2X,'THE INPUT DATA WILL NOT BE ECHO-PRINTED ') RETURN C C ECHO-PRINTED IS REQUIRED C 140 WRITE(6,145) 145 FORMAT(2X,'ECHO-PRINT OF THE INPUT DATA THAT WILL BE USED IN ', + 'THE ANALYSIS THAT FOLLOWS:') IECHO=1 WRITE(6,155)(LSTRING(I),I=1,LENGTH) 155 FORMAT(2X,80A1) RETURN 300 IEND=1 RETURN END C-------------------------------------------- SUBROUTINE READHL(KEY,ITITLE) C-------------------------------------------- C C THIS SUBROUTINE READS IN A FREE-TEXT DESCRIPTION OF THE DATA C AND A CODE WHICH INDICATES WHETHER LARGE VALUES OR C SMALL VALUES ARE ASSOCIATED WITH ACTUALLY POSITIVE CASES C COMMON/E/IEND,IERROR,IECHO COMMON/PASS/LSTRING(80),LINE(80),LENGTH DIMENSION ITITLE(*) DATA IBLANK/1H /,ISMALL1/1HS/,ISMALL2/1Hs/,ILARGE1/1HL/, + ILARGE2/1Hl/ C READ(5,100)(LINE(I),I=1,80) 100 FORMAT(80A1) DO 120 I=1,80 IF(LINE(I).EQ.IBLANK)GO TO 120 IPNT=I GO TO 130 120 CONTINUE WRITE(6,125) 125 FORMAT(2X,'--- NO INPUT FOR DATA DESCRIPTION ---') IPNT=1 130 DO 150 I=IPNT,80 ITITLE(I-IPNT+1)=LINE(I) 150 CONTINUE IF(IECHO.EQ.1)WRITE(6,155)(ITITLE(I),I=1,80) 155 FORMAT(2X,80A1) C 200 READ(5,100)(LINE(I),I=1,80) CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 210 IF(LSTRING(1).EQ.ISMALL1 .OR. LSTRING(1).EQ.ISMALL2)KEY=-1 IF(LSTRING(1).EQ.ILARGE1 .OR. LSTRING(1).EQ.ILARGE2)KEY=1 IF(IECHO.EQ.1)WRITE(6,205)(LSTRING(I),I=1,LENGTH), + (LINE(I),I=1,80) 205 FORMAT(2X,80A1) RETURN 210 IERROR=1 WRITE(6,215)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 215 FORMAT(2X,80A1/1X,'--- INVALID CODE WORD DETECTED ---') RETURN END C-------------------------------------------- SUBROUTINE READQ(NUM,F1,F2,IDENT) C-------------------------------------------- C C THIS SUBROUTINE READS IN A SEQUENCE OF INPUT DATA PAIR IN FREE C FORMAT. THE ONLY FORMAT REQUIREMENTS ARE THAT (1) C ANY TWO INPUT VALUES MUST BE SEPARATED BY AT LEAST ONE C IBLANK COLUMN AND (2) THE INPUT DATA SEQUENCE MUST BE TERMINATED C BY AN ASTERISK (*). C NEGATIVE VALUES ARE ACCEPTABLE. C COMMON/E/IEND,IERROR,IECHO COMMON/PASS/LSTRING(80),LINE(80),LENGTH COMMON/DECIML/MAXDX,MAXDY DIMENSION F1(*),F2(*),IDENT(1000,*) DATA IBLANK/1H /,ISTAR/1H*/ NUM=0 DO 50 I=1,1000 DO 50 J=1,80 IDENT(I,J)=IBLANK 50 CONTINUE C 100 READ(5,110)(LINE(I),I=1,80) 110 FORMAT(80A1) C C READ TEST VALUE FOR X-CONDITION C 120 CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 130 IF(LSTRING(1).EQ.ISTAR)RETURN CALL TONUM(MAXDX,LENGTH,LSTRING,RVALUE,IERROR) IF(IERROR.GT.0)GO TO 130 NUM=NUM+1 F1(NUM)=RVALUE C C READ TEST VALUE FOR Y-CONDITION C CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 130 CALL TONUM(MAXDY,LENGTH,LSTRING,RVALUE,IERROR) IF(IERROR.GT.0)GO TO 130 F2(NUM)=RVALUE C C READ TEXT-FREE DATA DESCRIPTION FOR EACH TEST-RESULT PAIR C DO 125 I=1,80 IF(LINE(I).EQ.IBLANK)GO TO 125 IPNT=I GO TO 128 125 CONTINUE GO TO 100 128 DO 129 I=IPNT,80 IDENT(NUM,I-IPNT+1)=LINE(I) 129 CONTINUE GO TO 100 C 130 WRITE(6,135)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 135 FORMAT(2X,80A1/1X,'--- INVALID TEST RESULT PAIR DETECTED ---') RETURN END C-------------------------------------------- SUBROUTINE READT(TCODE,FP) C-------------------------------------------- C C READ IN STATISTICAL TEST CODE WORD AND C FPF (IF 'TPF' WAS ENTERED) C COMMON/E/IEND,IERROR,IECHO COMMON/PASS/LSTRING(80),LINE(80),LENGTH DATA IB1/1HB/,IB2/1Hb/,IA1/1HA/,IA2/1Ha/,IT1/1HT/,IT2/1Ht/ INTEGER TCODE C TCODE=0 READ(5,100)(LINE(I),I=1,80) 100 FORMAT(80A1) CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 180 TCODE=99 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.NE.1 .AND. TCODE.NE.2 .AND. TCODE.NE.3)GO TO 180 IF(IECHO.EQ.1)WRITE(6,120)(LSTRING(I),I=1,LENGTH), + (LINE(I),I=1,80) 120 FORMAT(2X,80A1) IF(TCODE.NE.3)RETURN C C READ IN FPF VALUE IF TCODE=3 C MAXD=0 READ(5,100)(LINE(I),I=1,80) CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 190 CALL TONUM(MAXD,LENGTH,LSTRING,RVALUE,IERROR) IF(IERROR.GT.0)GO TO 190 FP=RVALUE IF(FP.GT.1.OR.FP.LT.0)GO TO 190 CALL GETWRD(NOTGET) IF(NOTGET.EQ.0)GO TO 190 IF(IECHO.EQ.1)WRITE(6,135)FP 135 FORMAT(2X,F5.3) RETURN 180 IERROR=2 WRITE(6,185)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 185 FORMAT(2X,80A1/ + 1X,'--- INVALID STATISTICAL TEST CODE WORD DETECTED ---') RETURN 190 IERROR=3 WRITE(6,195)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 195 FORMAT(2X,80A1/1X,'--- INVALID FPF VALUE DETECTED ---') 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(MAXD,LEN,LINPUT,RNUM,IERR) C----------------------------------------------------------- C C CONVERT CHARACTER STRING TO REAL NUMBER C DIMENSION LINPUT(*) DATA IDEC/1H./,IZERO/1H0/,NEG/1H-/ C IERR= 0 IF(LEN .EQ. 0) GO TO 120 C C CHECK INVALID CHARACTER (I.E., ANY CHARACTER EXCEPT 0-9, - OR .) C IP= 0 IN= 0 DO 20 I= 1,LEN IK= IORD(LINPUT(I)) - IORD(IZERO) IF(IK .GE. 0 .AND. IK .LE. 9) GO TO 20 IF(LINPUT(I) .NE. IDEC) GO TO 10 IP= IP+1 IF(IP.GT.1)GO TO 120 GO TO 20 10 IF(LINPUT(I) .NE. NEG) GO TO 120 IN= IN+1 IF(IN.GT.1)GO TO 120 20 CONTINUE IF(IP.EQ.1.AND.LEN.EQ.1)GO TO 120 IF(IN.EQ.1.AND.LEN.EQ.1)GO TO 120 C C CONVERT TO NUMERICAL STRING C IU= 0 TOTAL= 0 K= 1 IF(LINPUT(K) .EQ. NEG) K=K+1 40 IF(LINPUT(K) .EQ. IDEC) GO TO 70 IK= IORD(LINPUT(K)) - IORD(IZERO) IF(IK .LT. 0 .OR. IK .GT. 9) GO TO 120 TOTAL= TOTAL*10 + IK K= K+1 IF(K.GT.LEN) GO TO 100 GO TO 40 C 70 K= K+1 IF(K.GT.LEN) GO TO 100 IK= IORD(LINPUT(K)) - IORD(IZERO) IF(IK .LT. 0 .OR. IK .GT. 9) GO TO 120 TOTAL= TOTAL*10 + IK IU= IU+1 GO TO 70 C 100 RNUM = TOTAL/(10.0**IU) IF(MAXD.LT.IU)MAXD=IU IF(IN.EQ.1)RNUM=-RNUM RETURN C 120 IERR= 1 200 RETURN END C---------------------------------------------- INTEGER FUNCTION IORD(LCH) C---------------------------------------------- C C CHANGE FROM CHARACTER TO ASCII CODE C C Changed 6/8/93; put \\ for \ in Holl const. 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 PPRINT(M,F1,F2,iTEXT) C------------------------------------------- C C ECHO PRINT C DIMENSION F1(*),F2(*),itext(1000,*) COMMON/DECIML/MAXDX,MAXDY IF(MAXDX.EQ.0)GO TO 10 IF(MAXDX.EQ.1)GO TO 11 IF(MAXDX.EQ.2)GO TO 12 IF(MAXDX.EQ.3)GO TO 13 IF(MAXDX.EQ.4)GO TO 14 IF(MAXDX.EQ.5)GO TO 15 IF(MAXDX.EQ.6)GO TO 16 IF(MAXDX.EQ.7)GO TO 17 IF(MAXDY.EQ.0)WRITE(6,180)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.1)WRITE(6,181)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.2)WRITE(6,182)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.3)WRITE(6,183)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.4)WRITE(6,184)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.5)WRITE(6,185)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.6)WRITE(6,186)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.7)WRITE(6,187)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.GE.8)WRITE(6,188)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) WRITE(6,200) RETURN C 10 IF(MAXDY.EQ.0)WRITE(6,100)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.1)WRITE(6,101)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.2)WRITE(6,102)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.3)WRITE(6,103)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.4)WRITE(6,104)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.5)WRITE(6,105)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.6)WRITE(6,106)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.7)WRITE(6,107)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.GE.8)WRITE(6,108)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) WRITE(6,200) RETURN C 11 IF(MAXDY.EQ.0)WRITE(6,110)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.1)WRITE(6,111)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.2)WRITE(6,112)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.3)WRITE(6,113)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.4)WRITE(6,114)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.5)WRITE(6,115)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.6)WRITE(6,116)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.7)WRITE(6,117)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.GE.8)WRITE(6,118)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) WRITE(6,200) RETURN C 12 IF(MAXDY.EQ.0)WRITE(6,120)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.1)WRITE(6,121)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.2)WRITE(6,122)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.3)WRITE(6,123)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.4)WRITE(6,124)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.5)WRITE(6,125)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.6)WRITE(6,126)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.7)WRITE(6,127)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.GE.8)WRITE(6,128)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) WRITE(6,200) RETURN C 13 IF(MAXDY.EQ.0)WRITE(6,130)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.1)WRITE(6,131)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.2)WRITE(6,132)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.3)WRITE(6,133)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.4)WRITE(6,134)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.5)WRITE(6,135)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.6)WRITE(6,136)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.7)WRITE(6,137)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.GE.8)WRITE(6,138)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) WRITE(6,200) RETURN C 14 IF(MAXDY.EQ.0)WRITE(6,140)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.1)WRITE(6,141)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.2)WRITE(6,142)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.3)WRITE(6,143)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.4)WRITE(6,144)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.5)WRITE(6,145)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.6)WRITE(6,146)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.7)WRITE(6,147)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.GE.8)WRITE(6,148)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) WRITE(6,200) RETURN C 15 IF(MAXDY.EQ.0)WRITE(6,150)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.1)WRITE(6,151)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.2)WRITE(6,152)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.3)WRITE(6,153)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.4)WRITE(6,154)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.5)WRITE(6,155)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.6)WRITE(6,156)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.7)WRITE(6,157)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.GE.8)WRITE(6,158)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) WRITE(6,200) RETURN C 16 IF(MAXDY.EQ.0)WRITE(6,160)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.1)WRITE(6,161)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.2)WRITE(6,162)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.3)WRITE(6,163)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.4)WRITE(6,164)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.5)WRITE(6,165)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.6)WRITE(6,166)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.7)WRITE(6,167)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.GE.8)WRITE(6,168)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) WRITE(6,200) RETURN C 17 IF(MAXDY.EQ.0)WRITE(6,170)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.1)WRITE(6,171)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.2)WRITE(6,172)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.3)WRITE(6,173)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.4)WRITE(6,174)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.5)WRITE(6,175)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.6)WRITE(6,176)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.EQ.7)WRITE(6,177)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) IF(MAXDY.GE.8)WRITE(6,178)(F1(I),F2(I),(itext(I,J),J=1,80),I=1,M) WRITE(6,200) RETURN 100 FORMAT(3X,F8.0,2X,F8.0,5X,80A1) 101 FORMAT(3X,F8.0,2X,F9.1,5X,80A1) 102 FORMAT(3X,F8.0,2X,F10.2,5X,80A1) 103 FORMAT(3X,F8.0,2X,F11.3,5X,80A1) 104 FORMAT(3X,F8.0,2X,F12.4,5X,80A1) 105 FORMAT(3X,F8.0,2X,F13.5,5X,80A1) 106 FORMAT(3X,F8.0,2X,F14.6,5X,80A1) 107 FORMAT(3X,F8.0,2X,F15.7,5X,80A1) 108 FORMAT(3X,F8.0,2X,F16.8,5X,80A1) C 110 FORMAT(3X,F9.1,2X,F8.0,5X,80A1) 111 FORMAT(3X,F9.1,2X,F9.1,5X,80A1) 112 FORMAT(3X,F9.1,2X,F10.2,5X,80A1) 113 FORMAT(3X,F9.1,2X,F11.3,5X,80A1) 114 FORMAT(3X,F9.1,2X,F12.4,5X,80A1) 115 FORMAT(3X,F9.1,2X,F13.5,5X,80A1) 116 FORMAT(3X,F9.1,2X,F14.6,5X,80A1) 117 FORMAT(3X,F9.1,2X,F15.7,5X,80A1) 118 FORMAT(3X,F9.1,2X,F16.8,5X,80A1) C 120 FORMAT(3X,F10.2,2X,F8.0,5X,80A1) 121 FORMAT(3X,F10.2,2X,F9.1,5X,80A1) 122 FORMAT(3X,F10.2,2X,F10.2,5X,80A1) 123 FORMAT(3X,F10.2,2X,F11.3,5X,80A1) 124 FORMAT(3X,F10.2,2X,F12.4,5X,80A1) 125 FORMAT(3X,F10.2,2X,F13.5,5X,80A1) 126 FORMAT(3X,F10.2,2X,F14.6,5X,80A1) 127 FORMAT(3X,F10.2,2X,F15.7,5X,80A1) 128 FORMAT(3X,F10.2,2X,F16.8,5X,80A1) C 130 FORMAT(3X,F11.3,2X,F8.0,5X,80A1) 131 FORMAT(3X,F11.3,2X,F9.1,5X,80A1) 132 FORMAT(3X,F11.3,2X,F10.2,5X,80A1) 133 FORMAT(3X,F11.3,2X,F11.3,5X,80A1) 134 FORMAT(3X,F11.3,2X,F12.4,5X,80A1) 135 FORMAT(3X,F11.3,2X,F13.5,5X,80A1) 136 FORMAT(3X,F11.3,2X,F14.6,5X,80A1) 137 FORMAT(3X,F11.3,2X,F15.7,5X,80A1) 138 FORMAT(3X,F11.3,2X,F16.8,5X,80A1) C 140 FORMAT(3X,F12.4,2X,F8.0,5X,80A1) 141 FORMAT(3X,F12.4,2X,F9.1,5X,80A1) 142 FORMAT(3X,F12.4,2X,F10.2,5X,80A1) 143 FORMAT(3X,F12.4,2X,F11.3,5X,80A1) 144 FORMAT(3X,F12.4,2X,F12.4,5X,80A1) 145 FORMAT(3X,F12.4,2X,F13.5,5X,80A1) 146 FORMAT(3X,F12.4,2X,F14.6,5X,80A1) 147 FORMAT(3X,F12.4,2X,F15.7,5X,80A1) 148 FORMAT(3X,F12.4,2X,F16.8,5X,80A1) C 150 FORMAT(3X,F13.5,2X,F8.0,5X,80A1) 151 FORMAT(3X,F13.5,2X,F9.1,5X,80A1) 152 FORMAT(3X,F13.5,2X,F10.2,5X,80A1) 153 FORMAT(3X,F13.5,2X,F11.3,5X,80A1) 154 FORMAT(3X,F13.5,2X,F12.4,5X,80A1) 155 FORMAT(3X,F13.5,2X,F13.5,5X,80A1) 156 FORMAT(3X,F13.5,2X,F14.6,5X,80A1) 157 FORMAT(3X,F13.5,2X,F15.7,5X,80A1) 158 FORMAT(3X,F13.5,2X,F16.8,5X,80A1) C 160 FORMAT(3X,F14.6,2X,F8.0,5X,80A1) 161 FORMAT(3X,F14.6,2X,F9.1,5X,80A1) 162 FORMAT(3X,F14.6,2X,F10.2,5X,80A1) 163 FORMAT(3X,F14.6,2X,F11.3,5X,80A1) 164 FORMAT(3X,F14.6,2X,F12.4,5X,80A1) 165 FORMAT(3X,F14.6,2X,F13.5,5X,80A1) 166 FORMAT(3X,F14.6,2X,F14.6,5X,80A1) 167 FORMAT(3X,F14.6,2X,F15.7,5X,80A1) 168 FORMAT(3X,F14.6,2X,F16.8,5X,80A1) C 170 FORMAT(3X,F15.7,2X,F8.0,5X,80A1) 171 FORMAT(3X,F15.7,2X,F9.1,5X,80A1) 172 FORMAT(3X,F15.7,2X,F10.2,5X,80A1) 173 FORMAT(3X,F15.7,2X,F11.3,5X,80A1) 174 FORMAT(3X,F15.7,2X,F12.4,5X,80A1) 175 FORMAT(3X,F15.7,2X,F13.5,5X,80A1) 176 FORMAT(3X,F15.7,2X,F14.6,5X,80A1) 177 FORMAT(3X,F15.7,2X,F15.7,5X,80A1) 178 FORMAT(3X,F15.7,2X,F16.8,5X,80A1) C 180 FORMAT(3X,F16.8,2X,F8.0,5X,80A1) 181 FORMAT(3X,F16.8,2X,F9.1,5X,80A1) 182 FORMAT(3X,F16.8,2X,F10.2,5X,80A1) 183 FORMAT(3X,F16.8,2X,F11.3,5X,80A1) 184 FORMAT(3X,F16.8,2X,F12.4,5X,80A1) 185 FORMAT(3X,F16.8,2X,F13.5,5X,80A1) 186 FORMAT(3X,F16.8,2X,F14.6,5X,80A1) 187 FORMAT(3X,F16.8,2X,F15.7,5X,80A1) 188 FORMAT(3X,F16.8,2X,F16.8,5X,80A1) C 200 FORMAT(2X,'*') END C-------------------------------------------------- SUBROUTINE CONVRT(N,DATA,K,CUTOFF) C-------------------------------------------------- C C CONVERT QUANTITATIVE DATA INTO CATEGORICAL DATA C DIMENSION DATA(*),CUTOFF(*) DO 1040 I=1,N IF(DATA(I) .GT. CUTOFF(1)) GO TO 1010 DATA(I)=1 GO TO 1040 1010 DO 1030 J=1,K IF(.NOT.(DATA(I).GT.CUTOFF(J).AND.DATA(I).LE.CUTOFF(J+1))) + GO TO 1030 DATA(I)=J+1 GO TO 1040 1030 CONTINUE IF(DATA(I) .GT. CUTOFF(K)) DATA(I)=K+1 1040 CONTINUE RETURN END C------------------------------------------------------ SUBROUTINE SETMTX(NUM,AI,AJ,QQ) C------------------------------------------------------ C C SET UP THE RATING DATA MATRIX C DIMENSION QQ(21,21),AI(*),AJ(*) DO 110 I=1,21 DO 110 J=1,21 QQ(I,J)=0. 110 CONTINUE DO 120 I=1,NUM IX=AI(I) IY=AJ(I) QQ(IY,IX)=QQ(IY,IX)+1. 120 CONTINUE RETURN END C-------------------------------------------- SUBROUTINE SORT(NUM,F2) C-------------------------------------------- C C SORT A VECTOR C DIMENSION F2(2,*) C NN = NUM-1 DO 30 I= 1, NN SMALL= F2(1,I) JJ= I C DO 20 J= I,NUM IF (F2(1,J) .GE. SMALL) GO TO 20 SMALL = F2(1,J) JJ = J 20 CONTINUE C IF (JJ .EQ. I) GO TO 30 TEMP = F2(1,I) F2(1,I)= F2(1,JJ) F2(1,JJ)= TEMP TEMP = F2(2,I) F2(2,I)= F2(2,JJ) F2(2,JJ)= TEMP 30 CONTINUE RETURN END C----------------------------------------------- SUBROUTINE CORROC2 C----------------------------------------------- C C THIS PROGRAM USE THE ALGORITHM FROM 'LABROC1' TO CATEGORIZE C QUANTITATIVE DATA AND RUN THROUGH CORROC2. C COMMON/BLK1/VV(21,21),WW(21,21),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) COMMON/BLK3/TX(19),UY(19),TEXPX(21),TEXPY(21), 1 TXPTUR(21,21),TXPTUP(21,21),FYXP(21,21),FXYP(21,21), 2 FUTR(19,21),FTUR(21,19),ELR(21,21),ELP(21,21), 3 CPELR(20,20),CPELP(20,20),SKP,SKR,RADP,RADR,FUNLIK COMMON/BLK4/FOP(44),SOPNEG(44,44),VCOV(44,44),ESTCOR(44), 1 XXDUM(2000),ITER,CRIT,LSTAT COMMON/BLK5/FPVAL(26),TPVALX(26),TPVALY(26) COMMON/BLK6/ICLASS,JCLASS,IERX,IERY,ICON,JCON,IARY(21),JARY(21) COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(21),U(21) COMMON/E/IEND,IERROR,IECHO INTEGER TCODE REAL XI(21),XJ(21) IERX=0 IERY=0 ICLASS=0 JCLASS=0 CALL HEADLN(IERROR) 10 CALL CLAPSE CALL TPFFPF crit=1.0 C C CLASSIFY DATA SET DEGENERACY, IF ANY. IF DATA ARE NOT DEGENERATE, C GET INITIAL ESTIMATES OF THE AX,BX, T(I), AY, BY, AND U(J) C PARAMETERS FOR THE BIVARIATE CALCULATION BY RUNNING SUBROUTINE C 'DORALF' ON THE MARGINAL DATA SETS. C CALL CHECK(ICLASS,KATI,FPFX,TPFX,ICON,IARY) CALL CHECK(JCLASS,KATJ,FPFY,TPFY,JCON,JARY) KEYDEG=0 IF(ICLASS.NE.0.OR.JCLASS.NE.0)CALL DEGMSG(KEYDEG) IF(KEYDEG.EQ.1)RETURN CALL ESTM1(AX,BX,KATI,FPFX,TPFX,XI,EMN,EMS) CALL DORALF(AX,BX,KATI,EMN,EMS,SUMNT,SUMST,XI,VAX,VBX,VABX,IERX) IF(IERX.NE.0)THEN CALL ITRMSG RETURN ENDIF CALL ESTM1(AY,BY,KATJ,FPFY,TPFY,XJ,EMN,EMS) CALL DORALF(AY,BY,KATJ,EMN,EMS,SUMNU,SUMSU,XJ,VAY,VBY,VABY,IERY) IF (IERY.NE.0)THEN CALL ITRMSG RETURN ENDIF C C THIS PART CALCULATES MAXIMUM LIKELIHOOD ESTIMATES OF THE PARAMETERS C OF THE BIVARIATE MODEL FROM CORRELATED PAIRS OF RATINGS. C INITIAL ESTIMATES OF AX, BX, AY, BY AND THE T'S AND U'S WERE C OBTAINED ABOVE BY USING SUBROUTINE 'DORALF' ON THE MARGINAL DATA C SETS. INITIAL ESTIMATES OF THE CORRELATION COEFFICIENTS OF THE C UNDERLYING BIVARIATE NOISE AND SIGNAL-PLUS-NOISE DISTRIBUTIONS C WILL BE OBTAINED BY CALCULATING A PRODUCT MOMENT CORRELATION C COEFFICIENT FROM THE CORRESPONDING RATING-DATA MATRIX. C 20 ITER=0 LSTAT=0 IU=0 C C THE LSTAT VARIABLE IS A DUMMY VARIABLE FOR THE STATUS OF THE C ITERATIVE PROCEDURE. 0-INCOMPLETE , 1-COMPLETE C T(1)=0. U(1)=0. T(NNI)=0. U(NNJ)=0. DO 30 I=1,KKI T(I+1)=-XI(I) 30 CONTINUE DO 40 I=1,KKJ U(I+1)=-XJ(I) 40 CONTINUE CALL CORCOE AX0=AX BX0=BX AY0=AY BY0=BY R0=R RHO0=RHO CALL INIVAR 50 ITER=ITER+1 CALL TERMS CALL FIRST CALL SECOND CALL ITERAT(IU) IF(IU.EQ.1)GO TO 10 IF(IU.EQ.2)THEN CALL OUTINI(AX0,BX0,AY0,BY0,R0,RHO0) RETURN ENDIF IF (ITER.GT.100)THEN WRITE(6,60) CRIT 60 FORMAT(1X,'PROCEDURE DOES NOT CONVERGE AFTER 100 ITERATIONS'/ 1 1X,'ON THE LAST ITERATION CRIT=',F9.5///) RETURN ENDIF IF (LSTAT.NE.2) GO TO 50 CALL PRTOPS CALL TPFVAL CALL TEST(IERROR) CALL OUTRSL CALL INFORM(IERROR) RETURN END C-------------------------------------- BLOCK DATA C-------------------------------------- COMMON/BLK5/FPVAL(26),TPVALX(26),TPVALY(26) DATA FPVAL/0.005,0.01,0.02,0.03,0.04,0.05, 1 0.06,0.07,0.08,0.09,0.10,0.11, 2 0.12,0.13,0.14,0.15,0.20,0.25, 3 0.30,0.40,0.50,0.60,0.70,0.80, 4 0.90,0.95/ END C------------------------------------------------------------------ SUBROUTINE HEADLN(IERROR) C------------------------------------------------------------------ C C CALCULATE THE MARGINAL DATA FOR THE TWO EXPERIMENTAL CONDITIONS, C RESPECTIVELY, BY SUMMING THE 'NOISE' AND 'SINGLE-PLUS-NOISE' MATRICES C IN ROWS (FOR MARGINAL DATA OF CONDITION X) AND IN COLUMNS (FOR C MARGINAL DATA OF CONDITION Y). C COMMON/BLK1/VV(21,21),WW(21,21),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) INTEGER TCODE DO 40 I=1,KATI SUMNT(I)=0. SUMST(I)=0. 40 CONTINUE DO 45 J=1,KATJ SUMNU(J)=0. SUMSU(J)=0. 45 CONTINUE DO 50 I=1,KATI DO 50 J=1,KATJ SUMNT(I)=SUMNT(I)+VV(J,I) SUMST(I)=SUMST(I)+WW(J,I) 50 CONTINUE DO 55 J=1,KATJ DO 55 I=1,KATI SUMNU(J)=SUMNU(J)+VV(J,I) SUMSU(J)=SUMSU(J)+WW(J,I) 55 CONTINUE C C PRINT PAGE HEADING, IDENTIFICATION NAMES, AND RATING-DATA MATRICES C FOR NOISE CASES AND SIGNAL-PLUS-NOISE CASES C WRITE(6,70) 70 FORMAT(1H1,26X,'CLABROC (JUNE 1993 VERSION):'/ +15X,'MAXIMUM LIKELIHOOD ESTIMATION OF BINORMAL ROC CURVES'/ +15X,'FROM CORRELATED CONTINUOUSLY-DISTRIBUTED TEST RESULTS'/ +15X,' AND CALCULATION OF THE STATISTICAL SIGNIFICANCE OF'/ +15X,' ANY APPARENT DIFFERENCE BETWEEN THE CURVES.') WRITE(6,71) 71 FORMAT(///10X,'STATISTICAL TEST TO BE EMPLOYED:') IF(TCODE.EQ.1)WRITE(6,72) IF(TCODE.EQ.2)WRITE(6,73) IF(TCODE.EQ.3)WRITE(6,74)FP IF(IERROR.GT.1)WRITE(6,75) 72 FORMAT(16X,'CHI-SQUARE TEST FOR THE DIFFERENCE BETWEEN ', 1 '(AX,BX) AND (AY,BY)') 73 FORMAT(16X,'AREA (A SUB Z) TEST') 74 FORMAT(16X,'TPF TEST AT FPF = ',F7.4) 75 FORMAT(16X,'INVALID TEST CODE INPUT, BUT PROCEEDING WITH ', 1 'MAXIMUM LIKELIHOOD ESTIMATION') WRITE(6,80)(IDENTI(I),I=1,80) 80 FORMAT(//10X,'CONDITION X: ',80A1) IF(KSNX.EQ.KATI)WRITE(6,81) IF(KSNX.EQ.1)WRITE(6,82) 81 FORMAT(10X,'LARGER VALUES OF THE TEST RESULT REPRESENT ', + 'STRONGER EVIDENCE THAT THE CASE IS ACTUALLY ', + 'POSITIVE'/10X,'(E.G., THAT THE PATIENT IS ACTUALLY ', + 'ABNORMAL)') 82 FORMAT(10X,'SMALLER VALUES OF THE TEST RESULT REPRESENT ', + 'STRONGER EVIDENCE THAT THE CASE IS ACTUALLY ', + 'POSITIVE'/10X,'(E.G., THAT THE PATIENT IS ACTUALLY ', + 'ABNORMAL)') WRITE(6,85)(IDENTJ(I),I=1,80) 85 FORMAT(/10X,'CONDITION Y: ',80A1) IF(KSNY.EQ.KATJ)WRITE(6,81) IF(KSNY.EQ.1)WRITE(6,82) WRITE(6,90)EMN,EMS 90 FORMAT(/10X,'TOTAL NUMBER OF ACTUALLY NEGATIVE CASES = ',F4.0, + 7X,'TOTAL NUMBER OF ACTUALLY POSITIVE CASES = ',F4.0) RETURN END C----------------------------------------------------------------- SUBROUTINE CLAPSE C----------------------------------------------------------------- C C COLLAPSE DATA C COMMON/BLK1/VV(21,21),WW(21,21),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) INTEGER TCODE DO 40 I=1,KATI SUMNT(I)=0. SUMST(I)=0. 40 CONTINUE DO 45 J=1,KATJ SUMNU(J)=0. SUMSU(J)=0. 45 CONTINUE DO 50 I=1,KATI DO 50 J=1,KATJ SUMNT(I)=SUMNT(I)+VV(J,I) SUMST(I)=SUMST(I)+WW(J,I) 50 CONTINUE DO 55 J=1,KATJ DO 55 I=1,KATI SUMNU(J)=SUMNU(J)+VV(J,I) SUMSU(J)=SUMSU(J)+WW(J,I) 55 CONTINUE C KATI1=KATI DO 5005 I=1,KATI IT=KATI-I+1 IF (SUMNT(IT).NE.0.OR.SUMST(IT).NE.0) GO TO 5005 KATI1=KATI1-1 IF (IT.GT.KATI1) GO TO 5005 DO 5000 J=IT,KATI1 SUMNT(J)=SUMNT(J+1) SUMST(J)=SUMST(J+1) DO 5003 IX=1,KATJ VV(IX,J)=VV(IX,J+1) WW(IX,J)=WW(IX,J+1) 5003 CONTINUE 5000 CONTINUE 5005 CONTINUE KATI=KATI1 IF(KSNX.NE.1)KSNX=KATI C KATJ1=KATJ DO 5015 I=1,KATJ IT=KATJ-I+1 IF (SUMNU(IT).NE.0.OR.SUMSU(IT).NE.0) GO TO 5015 KATJ1=KATJ1-1 IF (IT.GT.KATJ1) GO TO 5015 DO 5010 J=IT,KATJ1 SUMNU(J)=SUMNU(J+1) SUMSU(J)=SUMSU(J+1) DO 5007 IY=1,KATI VV(J,IY)=VV(J+1,IY) WW(J,IY)=WW(J+1,IY) 5007 CONTINUE 5010 CONTINUE 5015 CONTINUE KATJ=KATJ1 IF(KSNY.NE.1)KSNY=KATJ KKI=KATI-1 KKJ=KATJ-1 NNI=KKI+2 NNJ=KKJ+2 KKIL=KKI-1 KKJL=KKJ-1 RETURN END C--------------------------------------------------------------- SUBROUTINE TPFFPF C--------------------------------------------------------------- C C CALCULATE OBSERVED TRUE POSITIVE FRACTIONS AND C FALSE POSITIVE FRACTIONS AND OUTPUT THE OBSERVED C OPERATING POINTS C COMMON/BLK1/VV(21,21),WW(21,21),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) INTEGER TCODE FPFX(NNI)=0. TPFX(NNI)=0. FPFX(KATI)=SUMNT(KATI) TPFX(KATI)=SUMST(KATI) DO 5020 I=2,KATI FPFX(KATI-I+1)=FPFX(KATI-I+2)+SUMNT(KATI-I+1) TPFX(KATI-I+1)=TPFX(KATI-I+2)+SUMST(KATI-I+1) 5020 CONTINUE FPFY(NNJ)=0. TPFY(NNJ)=0. FPFY(KATJ)=SUMNU(KATJ) TPFY(KATJ)=SUMSU(KATJ) DO 5025 I=2,KATJ FPFY(KATJ-I+1)=FPFY(KATJ-I+2)+SUMNU(KATJ-I+1) TPFY(KATJ-I+1)=TPFY(KATJ-I+2)+SUMSU(KATJ-I+1) 5025 CONTINUE DO 5030 I=1,KATI FPFX(I)=FPFX(I)/EMN TPFX(I)=TPFX(I)/EMS 5030 CONTINUE DO 5035 I=1,KATJ FPFY(I)=FPFY(I)/EMN TPFY(I)=TPFY(I)/EMS 5035 CONTINUE RETURN END C-------------------------------------------- SUBROUTINE PRTOPS C-------------------------------------------- C C PRINT OUT OPERATING POINTS C COMMON/BLK1/VV(21,21),WW(21,21),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) INTEGER TCODE WRITE(6,5090) 5090 FORMAT(//20X,'OPERATING POINTS CORRESPONDING TO THE ', + 'CATEGORIZED INPUT DATA FOR CONDITION X:'/) IF(NNI.GE.11)GO TO 5092 WRITE(6,5095) (FPFX(NNI-I+1),I=1,NNI) WRITE(6,5096) (TPFX(NNI-I+1),I=1,NNI) WRITE(6,5098) GO TO 6000 5092 WRITE(6,5095) (FPFX(NNI-I+1),I=1,10) WRITE(6,5096) (TPFX(NNI-I+1),I=1,10) IF(NNI.GE.21)GO TO 5094 WRITE(6,5095) (FPFX(NNI-I+1),I=11,NNI) WRITE(6,5096) (TPFX(NNI-I+1),I=11,NNI) WRITE(6,5098) GO TO 6000 5094 WRITE(6,5095) (FPFX(NNI-I+1),I=11,20) WRITE(6,5096) (TPFX(NNI-I+1),I=11,20) IF(NNI.EQ.21) WRITE(6,5097) FPFX(1),TPFX(1) 5095 FORMAT(10X,' FPF: ',10(F6.4,2X)) 5096 FORMAT(10X,' TPF: ',10(F6.4,2X)/) 5097 FORMAT(10X,' FPF: ',F6.4/10X,' TPF: ',F6.4//) 5098 FORMAT(//) 6000 WRITE(6,6010) 6010 FORMAT(20X,'OPERATING POINTS CORRESPONDING TO THE ', + 'CATEGORIZED INPUT DATA FOR CONDITION Y:'/) IF(NNJ.GE.11)GO TO 6012 WRITE(6,5095) (FPFY(NNJ-I+1),I=1,NNJ) WRITE(6,5096) (TPFY(NNJ-I+1),I=1,NNJ) WRITE(6,5098) RETURN 6012 WRITE(6,5095) (FPFY(NNJ-I+1),I=1,10) WRITE(6,5096) (TPFY(NNJ-I+1),I=1,10) IF(NNJ.GE.21)GO TO 6014 WRITE(6,5095) (FPFY(NNJ-I+1),I=11,NNJ) WRITE(6,5096) (TPFY(NNJ-I+1),I=11,NNJ) WRITE(6,5098) RETURN 6014 WRITE(6,5095) (FPFY(NNJ-I+1),I=11,20) WRITE(6,5096) (TPFY(NNJ-I+1),I=11,20) IF(NNJ.EQ.21) WRITE(6,5097) FPFY(1),TPFY(1) RETURN END C------------------------------------------------------------------ SUBROUTINE DEGMSG(KEYDEG) C------------------------------------------------------------------ C C WRITE OUT DEGENERACY MESSAGE FOR CONDITION X OR CONDITION Y C COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) COMMON/BLK6/ICLASS,JCLASS,IERX,IERY,ICON,JCON,IARY(21),JARY(21) 5100 IF (ICLASS.EQ.0) GO TO 5120 IF (ICLASS.GE.5.AND.ICLASS.LE.8) WRITE(6,5105) IF (ICLASS.LT.5.OR.ICLASS.GT.8) WRITE(6,5110) 5105 FORMAT(///,15X,'MARGINAL DATA FOR CONDITION X ARE DEGENERATE.') 5110 FORMAT(///,15X,'MARGINAL DATA FOR CONDITION X ARE DEGENERATE') CALL DEGENE(ICLASS,IARY,FPFX,TPFX) 5120 IF (JCLASS.EQ.0) GO TO 5135 IF (JCLASS.GE.5.AND.JCLASS.LE.8) WRITE(6,5125) IF (JCLASS.LT.5.OR.JCLASS.GT.8) WRITE(6,5130) 5125 FORMAT(///,15X,'MARGINAL DATA FOR CONDITION Y ARE DEGENERATE.') 5130 FORMAT(///,15X,'MARGINAL DATA FOR CONDITION Y ARE DEGENERATE') CALL DEGENE(JCLASS,JARY,FPFY,TPFY) 5135 WRITE(6,5140) 5140 FORMAT(//,15X,'EXECUTION OF THIS BIVARIATE DATA SET ABORTED.') KEYDEG=1 RETURN END C-------------------------------------------------------------------- SUBROUTINE ITRMSG C-------------------------------------------------------------------- C C WRITE OUT MESSAGE DESCRIBING ANY MATRIX INVERSION PROBLEM OR C REPORTING THAT THE NUMBER OF ITERATIONS EXCEEDED MAXIMUM C RESTRICTION FOR MARGINAL DATA SETS C COMMON/BLK6/ICLASS,JCLASS,IERX,IERY,ICON,JCON,IARY(21),JARY(21) IF (IERX.EQ.0) GO TO 5205 IF (IERX.GT.100) WRITE(6,5203) IF (IERX.LE.100) WRITE(6,5204) IERX 5203 FORMAT(1X,'PROCEDURE DOES NOT CONVERGE -- 100 ITERATIONS --', 1 'THIS WAS WITH MARGINAL DATA SET FOR CONDITION X') 5204 FORMAT(1X,'PROBLEM WITH MATRIX INVERSION ON ',I3,'TH ITERATION', 2 ' WITH MARGINAL DATA SET FOR CONDITION X') IF (IERY.EQ.0) RETURN 5205 IF (IERY.GT.100) WRITE(6,5206) IF (IERY.LE.100) WRITE(6,5207) IERY 5206 FORMAT(1X,'PROCEDURE DOES NOT CONVERGE -- 100 ITERATIONS --', 1 'THIS WAS WITH MARGINAL DATA SET FOR CONDITION Y') 5207 FORMAT(1X,'PROBLEM WITH MATRIX INVERSION ON ',I3,'TH ITERATION', 2 ' WITH MARGINAL DATA SET FOR CONDITION Y') RETURN END C----------------------------------------------------------------- SUBROUTINE CORCOE C----------------------------------------------------------------- C C CALCULATE INITIAL ESTIMATES FOR R=R(NOISE) AND RHO=R(SIGNAL-PLUS-NOISE) C BY COMPUTING PEARSON PRODUCT MOMENT CORRELATION COEFFICIENTS FOR THE C NOISE RATING-DATA MATRIX AND THE SIGNAL-PLUS-NOISE RATING-DATA MATRIX, C RESPECTIVELY. C COMMON/BLK1/VV(21,21),WW(21,21),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(21),U(21) INTEGER TCODE SUMRX=0. SUMPX=0. SUMRY=0. SUMPY=0. SUMRX2=0. SUMPX2=0. SUMRY2=0. SUMPY2=0. SUMRXY=0. SUMPXY=0. DO 9065 I=1,KATI SUMRX=SUMRX+SUMNT(I)*I SUMPX=SUMPX+SUMST(I)*I SUMRX2=SUMRX2+SUMNT(I)*I**2 SUMPX2=SUMPX2+SUMST(I)*I**2 9065 CONTINUE DO 9070 J=1,KATJ SUMRY=SUMRY+SUMNU(J)*J SUMPY=SUMPY+SUMSU(J)*J SUMRY2=SUMRY2+SUMNU(J)*J**2 SUMPY2=SUMPY2+SUMSU(J)*J**2 9070 CONTINUE DO 9075 I=1,KATI DO 9075 J=1,KATJ SUMRXY=SUMRXY+VV(J,I)*J*I SUMPXY=SUMPXY+WW(J,I)*J*I 9075 CONTINUE R=(EMN*SUMRXY-SUMRX*SUMRY)/SQRT((EMN*SUMRX2-SUMRX**2)*(EMN*SUMRY2 1-SUMRY**2)) RHO=(EMS*SUMPXY-SUMPX*SUMPY)/SQRT((EMS*SUMPX2-SUMPX**2)*(EMS* 1SUMPY2-SUMPY**2)) RETURN END C--------------------------------------------------------------------- SUBROUTINE INIVAR C--------------------------------------------------------------------- C C INITIALIZE VARIABLES C COMMON/BLK1/VV(21,21),WW(21,21),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) COMMON/BLK3/TX(19),UY(19),TEXPX(21),TEXPY(21), 1 TXPTUR(21,21),TXPTUP(21,21),FYXP(21,21),FXYP(21,21), 2 FUTR(19,21),FTUR(21,19),ELR(21,21),ELP(21,21), 3 CPELR(20,20),CPELP(20,20),SKP,SKR,RADP,RADR,FUNLIK COMMON/BLK4/FOP(44),SOPNEG(44,44),VCOV(44,44),ESTCOR(44), 1 XXDUM(2000),ITER,CRIT,LSTAT INTEGER TCODE C TEXPX(1)=0. TEXPX(NNI)=0. TEXPY(1)=0. TEXPY(NNI)=0. DO 9100 J=1,NNJ FYXP(1,J)=0. FYXP(NNI,J)=0. ELR(1,J)=0. ELP(1,J)=0. TXPTUR(1,J)=0. TXPTUR(NNI,J)=0. TXPTUP(1,J)=0. TXPTUP(NNI,J)=0. 9100 CONTINUE DO 9105 I=2,KATI FYXP(I,1)=0. FYXP(I,NNJ)=1. TXPTUR(I,1)=0. TXPTUR(I,NNJ)=0. TXPTUP(I,1)=0. TXPTUP(I,NNJ)=0. FUTR(I-1,NNJ)=1. FUTR(I-1,1)=0. 9105 CONTINUE DO 9110 I=1,NNI FXYP(I,1)=0. FXYP(I,NNJ)=0. 9110 CONTINUE DO 9115 J=2,KATJ FXYP(1,J)=0. FXYP(NNI,J)=1. FTUR(NNI,J-1)=1. FTUR(1,J-1)=0. 9115 CONTINUE DO 9120 I=2,NNI ELP(I,1)=0. ELR(I,1)=0. 9120 CONTINUE ELP(NNI,NNJ)=1. ELR(NNI,NNJ)=1. MTRX=NNI+NNJ+2 DO 9125 I=1,MTRX DO 9125 J=1,MTRX SOPNEG(I,J)=0. VCOV(I,J)=0. 9125 CONTINUE RETURN END C--------------------------------------------------------------------- SUBROUTINE TERMS C--------------------------------------------------------------------- C C COMPUTE THE MAJOR TERMS THAT OCCUR IN THE DERIVATIVE EXPRESSIONS. C COMMON/BLK1/VV(21,21),WW(21,21),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) COMMON/BLK3/TX(19),UY(19),TEXPX(21),TEXPY(21), 1 TXPTUR(21,21),TXPTUP(21,21),FYXP(21,21),FXYP(21,21), 2 FUTR(19,21),FTUR(21,19),ELR(21,21),ELP(21,21), 3 CPELR(20,20),CPELP(20,20),SKP,SKR,RADP,RADR,FUNLIK COMMON/BLK4/FOP(44),SOPNEG(44,44),VCOV(44,44),ESTCOR(44), 1 XXDUM(2000),ITER,CRIT,LSTAT COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(21),U(21) INTEGER TCODE IF (RHO.GT.0.9999) RHO=0.9999 IF (R.GT.0.9999) R=0.9999 IF (RHO.LT.-0.9999) RHO=-0.9999 IF (R.LT.-0.9999) R=-0.9999 SKP=1.-RHO**2 SKR=1.-R**2 RADP=SQRT(SKP) RADR=SQRT(SKR) DO 9135 I=2,KATI TX(I-1)=T(I)*BX-AX TEXPX(I)=EXP(-.5*TX(I-1)**2) 9135 CONTINUE DO 9140 J=2,KATJ UY(J-1)=U(J)*BY-AY TEXPY(J)=EXP(-.5*UY(J-1)**2) 9140 CONTINUE DO 9145 I=1,KKI DO 9145 J=2,KATJ ARG=(U(J)-R*T(I+1))/RADR CALL NDTR(ARG,FUTR(I,J),D) 9145 CONTINUE DO 9150 J=1,KKJ DO 9150 I=2,KATI ARG=(T(I)-R*U(J+1))/RADR CALL NDTR(ARG,FTUR(I,J),D) 9150 CONTINUE DO 9155 I=2,KATI DO 9155 J=2,KATJ TXPTUR(I,J)=EXP(-1./(2.*SKR)*(T(I)**2-2.*R*T(I)*U(J)+U(J)**2)) TXPTUP(I,J)=EXP(-1./(2.*SKP)*(TX(I-1)**2-2.*RHO*TX(I-1)*UY(J-1)+ 1UY(J-1)**2)) ARG1=(UY(J-1)-RHO*TX(I-1))/RADP CALL NDTR(ARG1,FYXP(I,J),D) ARG2=(TX(I-1)-RHO*UY(J-1))/RADP CALL NDTR(ARG2,FXYP(I,J),D) 9155 CONTINUE IF (ITER.GT.1) RETURN DO 9156 I=2,KATI CALL NDTR(TX(I-1),ELP(I,NNJ),D) 9156 CALL NDTR (T(I),ELR(I,NNJ),D) DO 9157 J=2,KATJ CALL NDTR(UY(J-1),ELP(NNI,J),D) 9157 CALL NDTR(U(J),ELR(NNI,J),D) DO 9158 I=2,KATI DO 9158 J=2,KATJ CALL MDBNOR(TX(I-1),UY(J-1),RHO,ELP(I,J),IER) 9158 CALL MDBNOR(T(I),U(J),R,ELR(I,J),IER) DO 9160 I=1,KATI DO 9160 J=1,KATJ CPELR(I,J)=ELR(I+1,J+1)+ELR(I,J)-ELR(I+1,J)-ELR(I,J+1) CPELP(I,J)=ELP(I+1,J+1)+ELP(I,J)-ELP(I+1,J)-ELP(I,J+1) IF(CPELR(I,J).LT.1.0E-07) CPELR(I,J)=1.0E-07 IF(CPELP(I,J).LT.1.0E-07) CPELP(I,J)=1.0E-07 9160 CONTINUE FUNLIK=VLOGL(KATI,KATJ,VV,WW,CPELR,CPELP) RETURN END C--------------------------------------------------------------------- SUBROUTINE FIRST C--------------------------------------------------------------------- C C COMPUTE THE VALUE OF EACH OF THE FIRST-ORDER PARTIAL DERIVATIVE C EXPRESSIONS C COMMON/BLK1/VV(21,21),WW(21,21),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) COMMON/BLK3/TX(19),UY(19),TEXPX(21),TEXPY(21), 1 TXPTUR(21,21),TXPTUP(21,21),FYXP(21,21),FXYP(21,21), 2 FUTR(19,21),FTUR(21,19),ELR(21,21),ELP(21,21), 3 CPELR(20,20),CPELP(20,20),SKP,SKR,RADP,RADR,FUNLIK COMMON/BLK4/FOP(44),SOPNEG(44,44),VCOV(44,44),ESTCOR(44), 1 XXDUM(2000),ITER,CRIT,LSTAT COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(21),U(21) INTEGER TCODE GPI=3.14159265 C C WITH RESPECT TO AX C 9162 SUM=0. DO 9163 I=1,KATI DO 9163 J=1,KATJ SUM=SUM+WW(J,I)*(TEXPX(I)*(FYXP(I,J+1)-FYXP(I,J))-TEXPX(I+1)* 1(FYXP(I+1,J+1)-FYXP(I+1,J)))/CPELP(I,J) 9163 CONTINUE FOP(1)=SUM/SQRT(2.*GPI) C C WITH RESPECT TO BX C SUM=0. DO 9165 I=1,KATI DO 9165 J=1,KATJ SUM=SUM+WW(J,I)*(-T(I)*TEXPX(I)*(FYXP(I,J+1)-FYXP(I,J))+T(I+1)* 1TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J)))/CPELP(I,J) 9165 CONTINUE FOP(2)=SUM/SQRT(2.*GPI) C C WITH RESPECT TO AY C SUM=0. DO 9170 I=1,KATI DO 9170 J=1,KATJ SUM=SUM+WW(J,I)*(TEXPY(J)*(FXYP(I+1,J)-FXYP(I,J))-TEXPY(J+1)* 1(FXYP(I+1,J+1)-FXYP(I,J+1)))/CPELP(I,J) 9170 CONTINUE FOP(3)=SUM/SQRT(2.*GPI) C C WITH RESPECT TO BY C SUM=0. DO 9175 I=1,KATI DO 9175 J=1,KATJ SUM=SUM+WW(J,I)*(-U(J)*TEXPY(J)*(FXYP(I+1,J)-FXYP(I,J))+U(J+1)* 1TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1)))/CPELP(I,J) 9175 CONTINUE FOP(4)=SUM/SQRT(2.*GPI) C C WITH RESPECT TO R=R(NOISE) C SUM=0. DO 9180 I=1,KATI DO 9180 J=1,KATJ SUM=SUM+VV(J,I)*(TXPTUR(I+1,J+1)+TXPTUR(I,J)-TXPTUR(I+1,J)- 1TXPTUR(I,J+1))/CPELR(I,J) 9180 CONTINUE FOP(5)=SUM/(2.*GPI*RADR) C C WITH RESPECT TO RHO=R(SIGNAL-PLUS-NOISE) C SUM=0. DO 9185 I=1,KATI DO 9185 J=1,KATJ SUM=SUM+WW(J,I)*(TXPTUP(I+1,J+1)+TXPTUP(I,J)-TXPTUP(I+1,J)- 1TXPTUP(I,J+1))/CPELP(I,J) 9185 CONTINUE FOP(6)=SUM/(2.*GPI*RADP) C C WITH RESPECT TO EACH OF THE T(I) C DO 9195 I=1,KKI SUMN=0. SUMS=0. DO 9190 J=1,KATJ SUMN=SUMN+EXP(-.5*T(I+1)**2)*(FUTR(I,J+1)-FUTR(I,J))*(VV(J,I)/ 1CPELR(I,J)-VV(J,I+1)/CPELR(I+1,J)) SUMS=SUMS+TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J))*(WW(J,I)/ 1CPELP(I,J)-WW(J,I+1)/CPELP(I+1,J)) 9190 CONTINUE FOP(I+6)=(SUMN+BX*SUMS)/SQRT(2.*GPI) 9195 CONTINUE C C WITH RESPECT TO EACH OF THE U(J) C DO 9205 J=1,KKJ SUMN=0. SUMS=0. DO 9200 I=1,KATI SUMN=SUMN+EXP(-.5*U(J+1)**2)*(FTUR(I+1,J)-FTUR(I,J))*(VV(J,I)/ 1CPELR(I,J)-VV(J+1,I)/CPELR(I,J+1)) SUMS=SUMS+TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1))*(WW(J,I)/ 1CPELP(I,J)-WW(J+1,I)/CPELP(I,J+1)) 9200 CONTINUE FOP(J+KKI+6)=(SUMN+BY*SUMS)/SQRT(2.*GPI) 9205 CONTINUE RETURN END C------------------------------------------------------------------ SUBROUTINE SECOND C------------------------------------------------------------------ C C COMPUTE THE EXPECTED VALUE OF EACH OF THE MINUS SECOND-ORDER C PARTIAL DERIVATIVE EXPRESSIONS. C COMMON/BLK1/VV(21,21),WW(21,21),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) COMMON/BLK3/TX(19),UY(19),TEXPX(21),TEXPY(21), 1 TXPTUR(21,21),TXPTUP(21,21),FYXP(21,21),FXYP(21,21), 2 FUTR(19,21),FTUR(21,19),ELR(21,21),ELP(21,21), 3 CPELR(20,20),CPELP(20,20),SKP,SKR,RADP,RADR,FUNLIK COMMON/BLK4/FOP(44),SOPNEG(44,44),VCOV(44,44),ESTCOR(44), 1 XXDUM(2000),ITER,CRIT,LSTAT COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(21),U(21) INTEGER TCODE GPI=3.14159265 C C WITH RESPECT TO AX C SUM=0. DO 9210 I=1,KATI DO 9210 J=1,KATJ SUM=SUM+(TEXPX(I)*(FYXP(I,J+1)-FYXP(I,J))-TEXPX(I+1)* 1(FYXP(I+1,J+1)-FYXP(I+1,J)))**2/CPELP(I,J) 9210 CONTINUE SOPNEG(1,1)=EMS*SUM/(2.*GPI) C C WITH RESPECT TO BX C SUM=0. DO 9215 I=1,KATI DO 9215 J=1,KATJ SUM=SUM+(-T(I)*TEXPX(I)*(FYXP(I,J+1)-FYXP(I,J))+T(I+1)*TEXPX(I+1)* 1(FYXP(I+1,J+1)-FYXP(I+1,J)))**2/CPELP(I,J) 9215 CONTINUE SOPNEG(2,2)=EMS*SUM/(2.*GPI) C C WITH RESPECT TO AY C SUM=0. DO 9220 I=1,KATI DO 9220 J=1,KATJ SUM=SUM+(TEXPY(J)*(FXYP(I+1,J)-FXYP(I,J))-TEXPY(J+1)* 1(FXYP(I+1,J+1)-FXYP(I,J+1)))**2/CPELP(I,J) 9220 CONTINUE SOPNEG(3,3)=EMS*SUM/(2.*GPI) C C WITH RESPECT TO BY C SUM=0. DO 9225 I=1,KATI DO 9225 J=1,KATJ SUM=SUM+(-U(J)*TEXPY(J)*(FXYP(I+1,J)-FXYP(I,J))+U(J+1)*TEXPY(J+1)* 1(FXYP(I+1,J+1)-FXYP(I,J+1)))**2/CPELP(I,J) 9225 CONTINUE SOPNEG(4,4)=EMS*SUM/(2.*GPI) C C WITH RESPECT TO R=R(NOISE) C SUM=0. DO 9230 I=1,KATI DO 9230 J=1,KATJ SUM=SUM+(TXPTUR(I+1,J+1)+TXPTUR(I,J)-TXPTUR(I+1,J)- 1TXPTUR(I,J+1))**2/CPELR(I,J) 9230 CONTINUE SOPNEG(5,5)=EMN*SUM/(4.*GPI**2*SKR) C C WITH RESPECT TO RHO=R(SIGNAL-PLUS-NOISE) C SUM=0. DO 9235 I=1,KATI DO 9235 J=1,KATJ SUM=SUM+(TXPTUP(I+1,J+1)+TXPTUP(I,J)-TXPTUP(I+1,J)- 1TXPTUP(I,J+1))**2/CPELP(I,J) 9235 CONTINUE SOPNEG(6,6)=EMS*SUM/(4.*GPI**2*SKP) C C WITH RESPECT TO EACH OF THE T(I) C DO 9245 I=1,KKI SUMN=0. SUMS=0. DO 9240 J=1,KATJ SUMN=SUMN+(EXP(-.5*T(I+1)**2)*(FUTR(I,J+1)-FUTR(I,J)))**2*(1./ 1CPELR(I,J)+1./CPELR(I+1,J)) SUMS=SUMS+(TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J)))**2*(1./ 1CPELP(I,J)+1./CPELP(I+1,J)) 9240 CONTINUE SOPNEG(I+6,I+6)=(EMN*SUMN+EMS*BX**2*SUMS)/(2.*GPI) 9245 CONTINUE C C WITH RESPECT TO EACH OF THE U(J) C DO 9255 J=1,KKJ SUMN=0. SUMS=0. DO 9250 I=1,KATI SUMN=SUMN+(EXP(-.5*U(J+1)**2)*(FTUR(I+1,J)-FTUR(I,J)))**2*(1./ 1CPELR(I,J)+1./CPELR(I,J+1)) SUMS=SUMS+(TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1)))**2*(1./ 1CPELP(I,J)+1./CPELP(I,J+1)) 9250 CONTINUE SOPNEG(J+KKI+6,J+KKI+6)=(EMN*SUMN+EMS*BY**2*SUMS)/(2.*GPI) 9255 CONTINUE C C MIXED WITH RESPECT TO AX AND BX C SUM=0. DO 9260 I=1,KATI DO 9260 J=1,KATJ SUM=SUM+(TEXPX(I)*(FYXP(I,J+1)-FYXP(I,J))-TEXPX(I+1)* 1(FYXP(I+1,J+1)-FYXP(I+1,J)))*(-T(I)*TEXPX(I)*(FYXP(I,J+1)- 2FYXP(I,J))+T(I+1)*TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J)))/ 3CPELP(I,J) 9260 CONTINUE SOPNEG(1,2)=EMS*SUM/(2.*GPI) SOPNEG(2,1)=SOPNEG(1,2) C C MIXED WITH RESPECT TO AX AND AY C SUM=0. DO 9265 I=1,KATI DO 9265 J=1,KATJ SUM=SUM+(TEXPX(I)*(FYXP(I,J+1)-FYXP(I,J))-TEXPX(I+1)* 1(FYXP(I+1,J+1)-FYXP(I+1,J)))*(TEXPY(J)*(FXYP(I+1,J)-FXYP(I,J))- 2TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1)))/CPELP(I,J) 9265 CONTINUE SOPNEG(1,3)=EMS*SUM/(2.*GPI) SOPNEG(3,1)=SOPNEG(1,3) C C MIXED WITH RESPECT TO AX AND BY C SUM=0. DO 9270 I=1,KATI DO 9270 J=1,KATJ SUM=SUM+(TEXPX(I)*(FYXP(I,J+1)-FYXP(I,J))-TEXPX(I+1)* 1(FYXP(I+1,J+1)-FYXP(I+1,J)))*(-U(J)*TEXPY(J)*(FXYP(I+1,J)- 2FXYP(I,J))+U(J+1)*TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1)))/ 3CPELP(I,J) 9270 CONTINUE SOPNEG(1,4)=EMS*SUM/(2.*GPI) SOPNEG(4,1)=SOPNEG(1,4) C C MIXED WITH RESPECT TO AX AND RHO=R(SIGNAL-PLUS-NOISE) C SUM=0. DO 9275 I=1,KATI DO 9275 J=1,KATJ SUM=SUM+(TEXPX(I)*(FYXP(I,J+1)-FYXP(I,J))-TEXPX(I+1)* 1(FYXP(I+1,J+1)-FYXP(I+1,J)))*(TXPTUP(I+1,J+1)+TXPTUP(I,J)- 2TXPTUP(I+1,J)-TXPTUP(I,J+1))/CPELP(I,J) 9275 CONTINUE SOPNEG(1,6)=EMS*SUM/SQRT(8.*GPI**3*SKP) SOPNEG(6,1)=SOPNEG(1,6) C C MIXED WITH RESPECT TO AX AND EACH OF THE T(I) C DO 9285 I=1,KKI SUM=0. DO 9280 J=1,KATJ SUM=SUM+TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J))*((TEXPX(I)* 1(FYXP(I,J+1)-FYXP(I,J))-TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J)))/ 2CPELP(I,J)-(TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J))-TEXPX(I+2)* 3(FYXP(I+2,J+1)-FYXP(I+2,J)))/CPELP(I+1,J)) 9280 CONTINUE SOPNEG(1,I+6)=EMS*BX*SUM/(2.*GPI) SOPNEG(I+6,1)=SOPNEG(1,I+6) 9285 CONTINUE C C MIXED WITH RESPECT TO AX AND EACH OF THE U(J) C DO 9295 J=1,KKJ SUM=0. DO 9290 I=1,KATI SUM=SUM+TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1))*((TEXPX(I)* 1(FYXP(I,J+1)-FYXP(I,J))-TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J)))/ 2CPELP(I,J)-(TEXPX(I)*(FYXP(I,J+2)-FYXP(I,J+1))-TEXPX(I+1)* 3(FYXP(I+1,J+2)-FYXP(I+1,J+1)))/CPELP(I,J+1)) 9290 CONTINUE SOPNEG(1,J+KKI+6)=EMS*BY*SUM/(2.*GPI) SOPNEG(J+KKI+6,1)=SOPNEG(1,J+KKI+6) 9295 CONTINUE C C MIXED WITH RESPECT TO BX AND AY C SUM=0. DO 9300 I=1,KATI DO 9300 J=1,KATJ SUM=SUM+(-T(I)*TEXPX(I)*(FYXP(I,J+1)-FYXP(I,J))+T(I+1)*TEXPX(I+1)* 1(FYXP(I+1,J+1)-FYXP(I+1,J)))*(TEXPY(J)*(FXYP(I+1,J)-FXYP(I,J))- 2TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1)))/CPELP(I,J) 9300 CONTINUE SOPNEG(2,3)=EMS*SUM/(2.*GPI) SOPNEG(3,2)=SOPNEG(2,3) C C MIXED WITH RESPECT TO BX AND BY C SUM=0. DO 9305 I=1,KATI DO 9305 J=1,KATJ SUM=SUM+(-T(I)*TEXPX(I)*(FYXP(I,J+1)-FYXP(I,J))+T(I+1)*TEXPX(I+1)* 1(FYXP(I+1,J+1)-FYXP(I+1,J)))*(-U(J)*TEXPY(J)*(FXYP(I+1,J)- 2FXYP(I,J))+U(J+1)*TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1)))/ 3CPELP(I,J) 9305 CONTINUE SOPNEG(2,4)=EMS*SUM/(2.*GPI) SOPNEG(4,2)=SOPNEG(2,4) C C MIXED WITH RESPECT TO BX AND RHO=R(SIGNAL-PLUS-NOISE) C SUM=0. DO 9310 I=1,KATI DO 9310 J=1,KATJ SUM=SUM+(-T(I)*TEXPX(I)*(FYXP(I,J+1)-FYXP(I,J))+T(I+1)*TEXPX(I+1)* 1(FYXP(I+1,J+1)-FYXP(I+1,J)))*(TXPTUP(I+1,J+1)+TXPTUP(I,J)- 2TXPTUP(I+1,J)-TXPTUP(I,J+1))/CPELP(I,J) 9310 CONTINUE SOPNEG(2,6)=EMS*SUM/SQRT(8.*GPI**3*SKP) SOPNEG(6,2)=SOPNEG(2,6) C C MIXED WITH RESPECT TO BX AND EACH OF THE T(I) C DO 9320 I=1,KKI SUM=0. DO 9315 J=1,KATJ SUM=SUM+TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J))*((-T(I)*TEXPX(I)* 1(FYXP(I,J+1)-FYXP(I,J))+T(I+1)*TEXPX(I+1)*(FYXP(I+1,J+1)- 2FYXP(I+1,J)))/CPELP(I,J)-(-T(I+1)*TEXPX(I+1)*(FYXP(I+1,J+1)- 3FYXP(I+1,J))+T(I+2)*TEXPX(I+2)*(FYXP(I+2,J+1)-FYXP(I+2,J)))/ 4CPELP(I+1,J)) 9315 CONTINUE SOPNEG(2,I+6)=EMS*BX*SUM/(2.*GPI) SOPNEG(I+6,2)=SOPNEG(2,I+6) 9320 CONTINUE C C MIXED WITH RESPECT TO BX AND EACH OF THE U(J) C DO 9330 J=1,KKJ SUM=0. DO 9325 I=1,KATI SUM=SUM+TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1))*((-T(I)*TEXPX(I)* 1(FYXP(I,J+1)-FYXP(I,J))+T(I+1)*TEXPX(I+1)*(FYXP(I+1,J+1)- 2FYXP(I+1,J)))/CPELP(I,J)-(-T(I)*TEXPX(I)*(FYXP(I,J+2)-FYXP(I,J+1)) 3+T(I+1)*TEXPX(I+1)*(FYXP(I+1,J+2)-FYXP(I+1,J+1)))/CPELP(I,J+1)) 9325 CONTINUE SOPNEG(2,J+KKI+6)=EMS*BY*SUM/(2.*GPI) SOPNEG(J+KKI+6,2)=SOPNEG(2,J+KKI+6) 9330 CONTINUE C C MIXED WITH RESPECT TO AY AND BY C SUM=0. DO 9335 I=1,KATI DO 9335 J=1,KATJ SUM=SUM+(TEXPY(J)*(FXYP(I+1,J)-FXYP(I,J))-TEXPY(J+1)* 1(FXYP(I+1,J+1)-FXYP(I,J+1)))*(-U(J)*TEXPY(J)*(FXYP(I+1,J)- 2FXYP(I,J))+U(J+1)*TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1)))/ 3CPELP(I,J) 9335 CONTINUE SOPNEG(3,4)=EMS*SUM/(2.*GPI) SOPNEG(4,3)=SOPNEG(3,4) C C MIXED WITH RESPECT TO AY AND RHO=R(SIGNAL-PLUS-NOISE) C SUM=0. DO 9340 I=1,KATI DO 9340 J=1,KATJ SUM=SUM+(TEXPY(J)*(FXYP(I+1,J)-FXYP(I,J))-TEXPY(J+1)* 1(FXYP(I+1,J+1)-FXYP(I,J+1)))*(TXPTUP(I+1,J+1)+TXPTUP(I,J)- 2TXPTUP(I+1,J)-TXPTUP(I,J+1))/CPELP(I,J) 9340 CONTINUE SOPNEG(3,6)=EMS*SUM/SQRT(8.*GPI**3*SKP) SOPNEG(6,3)=SOPNEG(3,6) C C MIXED WITH RESPECT TO AY AND EACH OF THE T(I) C DO 9350 I=1,KKI SUM=0. DO 9345 J=1,KATJ SUM=SUM+TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J))*((TEXPY(J)* 1(FXYP(I+1,J)-FXYP(I,J))-TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1)))/ 2CPELP(I,J)-(TEXPY(J)*(FXYP(I+2,J)-FXYP(I+1,J))-TEXPY(J+1)* 3(FXYP(I+2,J+1)-FXYP(I+1,J+1)))/CPELP(I+1,J)) 9345 CONTINUE SOPNEG(3,I+6)=EMS*BX*SUM/(2.*GPI) SOPNEG(I+6,3)=SOPNEG(3,I+6) 9350 CONTINUE C C MIXED WITH RESPECT TO AY AND EACH OF THE U(J) C DO 9360 J=1,KKJ SUM=0. DO 9355 I=1,KATI SUM=SUM+TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1))*((TEXPY(J)* 1(FXYP(I+1,J)-FXYP(I,J))-TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1)))/ 2CPELP(I,J)-(TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1))-TEXPY(J+2)* 3(FXYP(I+1,J+2)-FXYP(I,J+2)))/CPELP(I,J+1)) 9355 CONTINUE SOPNEG(3,J+KKI+6)=EMS*BY*SUM/(2.*GPI) SOPNEG(J+KKI+6,3)=SOPNEG(3,J+KKI+6) 9360 CONTINUE C C MIXED WITH RESPECT TO BY AND RHO=R(SIGNAL-PLUS-NOISE) C SUM=0. DO 9365 I=1,KATI DO 9365 J=1,KATJ SUM=SUM+(-U(J)*TEXPY(J)*(FXYP(I+1,J)-FXYP(I,J))+U(J+1)*TEXPY(J+1)* 1(FXYP(I+1,J+1)-FXYP(I,J+1)))*(TXPTUP(I+1,J+1)+TXPTUP(I,J)- 2TXPTUP(I+1,J)-TXPTUP(I,J+1))/CPELP(I,J) 9365 CONTINUE SOPNEG(4,6)=EMS*SUM/SQRT(8.*GPI**3*SKP) SOPNEG(6,4)=SOPNEG(4,6) C C MIXED WITH RESPECT TO BY AND EACH OF THE T(I) C DO 9375 I=1,KKI SUM=0. DO 9370 J=1,KATJ SUM=SUM+TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J))*((-U(J)*TEXPY(J)* 1(FXYP(I+1,J)-FXYP(I,J))+U(J+1)*TEXPY(J+1)*(FXYP(I+1,J+1)- 2FXYP(I,J+1)))/CPELP(I,J)-(-U(J)*TEXPY(J)*(FXYP(I+2,J)-FXYP(I+1,J)) 3+U(J+1)*TEXPY(J+1)*(FXYP(I+2,J+1)-FXYP(I+1,J+1)))/CPELP(I+1,J)) 9370 CONTINUE SOPNEG(4,I+6)=EMS*BX*SUM/(2.*GPI) SOPNEG(I+6,4)=SOPNEG(4,I+6) 9375 CONTINUE C C MIXED WITH RESPECT TO BY AND EACH OF THE U(J) C DO 9385 J=1,KKJ SUM=0. DO 9380 I=1,KATI SUM=SUM+TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1))*((-U(J)*TEXPY(J)* 1(FXYP(I+1,J)-FXYP(I,J))+U(J+1)*TEXPY(J+1)*(FXYP(I+1,J+1)- 2FXYP(I,J+1)))/CPELP(I,J)-(-U(J+1)*TEXPY(J+1)*(FXYP(I+1,J+1)- 3FXYP(I,J+1))+U(J+2)*TEXPY(J+2)*(FXYP(I+1,J+2)-FXYP(I,J+2)))/ 4CPELP(I,J+1)) 9380 CONTINUE SOPNEG(4,J+KKI+6)=EMS*BY*SUM/(2.*GPI) SOPNEG(J+KKI+6,4)=SOPNEG(4,J+KKI+6) 9385 CONTINUE C C MIXED WITH RESPECT TO R=R(NOISE) AND EACH OF THE T(I) C DO 9395 I=1,KKI SUM=0. DO 9390 J=1,KATJ SUM=SUM+EXP(-.5*T(I+1)**2)*(FUTR(I,J+1)-FUTR(I,J))* 1((TXPTUR(I+1,J+1)+TXPTUR(I,J)-TXPTUR(I+1,J)-TXPTUR(I,J+1))/ 2CPELR(I,J)-(TXPTUR(I+2,J+1)+TXPTUR(I+1,J)-TXPTUR(I+2,J)- 3TXPTUR(I+1,J+1))/CPELR(I+1,J)) 9390 CONTINUE SOPNEG(5,I+6)=EMN*SUM/SQRT(8.*GPI**3*SKR) SOPNEG(I+6,5)=SOPNEG(5,I+6) 9395 CONTINUE C C MIXED WITH RESPECT TO R=R(NOISE) AND EACH OF THE U(J) C DO 9405 J=1,KKJ SUM=0. DO 9400 I=1,KATI SUM=SUM+EXP(-.5*U(J+1)**2)*(FTUR(I+1,J)-FTUR(I,J))* 1((TXPTUR(I+1,J+1)+TXPTUR(I,J)-TXPTUR(I+1,J)-TXPTUR(I,J+1))/ 2CPELR(I,J)-(TXPTUR(I+1,J+2)+TXPTUR(I,J+1)-TXPTUR(I+1,J+1)- 3TXPTUR(I,J+2))/CPELR(I,J+1)) 9400 CONTINUE SOPNEG(5,J+KKI+6)=EMN*SUM/SQRT(8.*GPI**3*SKR) SOPNEG(J+KKI+6,5)=SOPNEG(5,J+KKI+6) 9405 CONTINUE C C MIXED WITH RESPECT TO RHO=R(SIGNAL-PLUS-NOISE) AND EACH OF THE T(I) C DO 9415 I=1,KKI SUM=0. DO 9410 J=1,KATJ SUM=SUM+TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J))*((TXPTUP(I+1,J+1)+ 1TXPTUP(I,J)-TXPTUP(I+1,J)-TXPTUP(I,J+1))/CPELP(I,J)- 2(TXPTUP(I+2,J+1)+TXPTUP(I+1,J)-TXPTUP(I+2,J)-TXPTUP(I+1,J+1))/ 3CPELP(I+1,J))