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)) 9410 CONTINUE SOPNEG(6,I+6)=EMS*BX*SUM/SQRT(8.*GPI**3*SKP) SOPNEG(I+6,6)=SOPNEG(6,I+6) 9415 CONTINUE C C MIXED WITH RESPECT TO RHO=R(SIGNAL-PLUS-NOISE) AND EACH OF THE U(J) C DO 9425 J=1,KKJ SUM=0. DO 9420 I=1,KATI SUM=SUM+TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1))*((TXPTUP(I+1,J+1)+ 1TXPTUP(I,J)-TXPTUP(I+1,J)-TXPTUP(I,J+1))/CPELP(I,J)- 2(TXPTUP(I+1,J+2)+TXPTUP(I,J+1)-TXPTUP(I+1,J+1)-TXPTUP(I,J+2))/ 3CPELP(I,J+1)) 9420 CONTINUE SOPNEG(6,J+KKI+6)=EMS*BY*SUM/SQRT(8.*GPI**3*SKP) SOPNEG(J+KKI+6,6)=SOPNEG(6,J+KKI+6) 9425 CONTINUE C C MIXED WITH RESPECT TO T(I) AND T(I+1) FOR I=1,KKI-1 C DO 9435 I=1,KKIL SUMN=0. SUMS=0. DO 9430 J=1,KATJ SUMN=SUMN+(-EXP(-.5*T(I+1)**2)*(FUTR(I,J+1)-FUTR(I,J)))*(EXP(-.5* 1T(I+2)**2)*(FUTR(I+1,J+1)-FUTR(I+1,J)))/CPELR(I+1,J) SUMS=SUMS+(-TEXPX(I+1)*(FYXP(I+1,J+1)-FYXP(I+1,J)))*(TEXPX(I+2)* 1(FYXP(I+2,J+1)-FYXP(I+2,J)))/CPELP(I+1,J) 9430 CONTINUE SOPNEG(I+6,I+7)=(EMN*SUMN+EMS*BX**2*SUMS)/(2.*GPI) SOPNEG(I+7,I+6)=SOPNEG(I+6,I+7) 9435 CONTINUE C C MIXED WITH RESPECT TO EACH OF THE T(I) AND EACH OF THE U(J) C DO 9440 I=1,KKI DO 9440 J=1,KKJ SUMN=EXP(-.5*U(J+1)**2)*(FTUR(I+1,J)-FTUR(I,J))*(EXP(-.5* 1T(I+1)**2)*(FUTR(I,J+1)-FUTR(I,J))/CPELR(I,J)-EXP(-.5*T(I+1)**2)* 2(FUTR(I,J+2)-FUTR(I,J+1))/CPELR(I,J+1))-EXP(-.5*U(J+1)**2)* 3(FTUR(I+2,J)-FTUR(I+1,J))*(EXP(-.5*T(I+1)**2)*(FUTR(I,J+1)- 4FUTR(I,J))/CPELR(I+1,J)-EXP(-.5*T(I+1)**2)*(FUTR(I,J+2)- 5FUTR(I,J+1))/CPELR(I+1,J+1)) SUMS=BY*TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1))*(BX*TEXPX(I+1)* 1(FYXP(I+1,J+1)-FYXP(I+1,J))/CPELP(I,J)-BX*TEXPX(I+1)* 2(FYXP(I+1,J+2)-FYXP(I+1,J+1))/CPELP(I,J+1))-BY*TEXPY(J+1)* 3(FXYP(I+2,J+1)-FXYP(I+1,J+1))*(BX*TEXPX(I+1)*(FYXP(I+1,J+1)- 4FYXP(I+1,J))/CPELP(I+1,J)-BX*TEXPX(I+1)*(FYXP(I+1,J+2)- 5FYXP(I+1,J+1))/CPELP(I+1,J+1)) SOPNEG(I+6,J+KKI+6)=(EMN*SUMN+EMS*SUMS)/(2.*GPI) SOPNEG(J+KKI+6,I+6)=SOPNEG(I+6,J+KKI+6) 9440 CONTINUE C C MIXED WITH RESPECT TO U(J) AND U(J+1)FOR J=1,KKJ-1 C DO 9450 J=1,KKJL SUMN=0. SUMS=0. DO 9445 I=1,KATI SUMN=SUMN+(-EXP(-.5*U(J+1)**2)*(FTUR(I+1,J)-FTUR(I,J)))*(EXP(-.5* 1U(J+2)**2)*(FTUR(I+1,J+1)-FTUR(I,J+1)))/CPELR(I,J+1) SUMS=SUMS+(-TEXPY(J+1)*(FXYP(I+1,J+1)-FXYP(I,J+1)))*(TEXPY(J+2)* 1(FXYP(I+1,J+2)-FXYP(I,J+2)))/CPELP(I,J+1) 9445 CONTINUE SOPNEG(J+KKI+6,J+KKI+7)=(EMN*SUMN+EMS*BY**2*SUMS)/(2.*GPI) SOPNEG(J+KKI+7,J+KKI+6)=SOPNEG(J+KKI+6,J+KKI+7) 9450 CONTINUE RETURN END C---------------------------------------------------------------- SUBROUTINE SELCEL(IX,IY) C---------------------------------------------------------------- C C SELECT A TARGET CELL WHICH HAS THE SMALLEST PROBABILITY AND C NONZERO COUNT C C Fixed: Identi & identJ = 80 not 60 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/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 C DO 1010 I=1,KATI DO 1010 J=1,KATJ IF(VV(J,I).LT.1)GO TO 1010 IX=I IY=J SMALL=CPELR(IX,IY) GO TO 1020 1010 CONTINUE 1020 DO 1030 I=1,KATI DO 1030 J=1,KATJ IF(WW(J,I).GT.0 .AND. (CPELP(I,J).LE.SMALL))THEN IX=I IY=J SMALL=CPELP(IX,IY) ENDIF IF(VV(J,I).GT.0 .AND. (CPELR(I,J).LE.SMALL))THEN IX=I IY=J SMALL=CPELR(IX,IY) ENDIF 1030 CONTINUE RETURN END C---------------------------------------------------------------- FUNCTION VLOGL(KATI,KATJ,VV,WW,CPELR,CPELP) C---------------------------------------------------------------- C C COMPUTE THE VALUE OF LOG L C DIMENSION VV(21,21),WW(21,21),CPELR(20,20),CPELP(20,20) SUMN=0. SUMS=0. DO 1500 I=1,KATI DO 1500 J=1,KATJ SUMN=SUMN+VV(J,I)*ALOG(CPELR(I,J)) SUMS=SUMS+WW(J,I)*ALOG(CPELP(I,J)) 1500 CONTINUE VLOGL=SUMN+SUMS RETURN END C-------------------------------------------- SUBROUTINE TRUNCT(IX,IY,IU) C-------------------------------------------- C C COLLAPSE THE (I,J)-TH ELEMENT OF THE RATING DATA MATRIX C C 6/15/93 Fixed Identi & identJ =80 not 60 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/OP/ZZX(20),ZZY(20) DIMENSION FPF(21),TPF(21),D(3,4),IARY(21) NN=0 IU=1 C C CALCULATE THE MINIMUM DISTANCE BETWEEN TWO OPERATION POINTS C --- FOR CONDITION X C IF(IX.EQ.1)THEN NN=NN+1 D(1,NN)=DIST(1,2,3) D2=DIST(1,1,2) IF(D2.LT.D(1,NN))D(1,NN)=D2 D(2,NN)=IX+1 D(3,NN)=-1 ELSE IF(IX.EQ.KATI)THEN NN=NN+1 D(1,NN)=DIST(1,KATI,NNI) D2=DIST(1,KKI,KATI) IF(D2.LT.D(1,NN))D(1,NN)=D2 D(2,NN)=IX D(3,NN)=-1 ELSE NN=NN+1 D(1,NN)=DIST(1,IX,IX-1) D(2,NN)=IX D(3,NN)=-1 NN=NN+1 D(1,NN)=DIST(1,IX+1,IX+2) D(2,NN)=IX+1 D(3,NN)=-1 ENDIF C C CALCULATE THE MINIMUN DISTANCE BETWEEN TWO OPERATING POINTS C --- FOR CONDITION Y C IF(IY.EQ.1)THEN NN=NN+1 D(1,NN)=DIST(2,2,3) D2=DIST(2,1,2) IF(D2.LT.D(1,NN))D(1,NN)=D2 D(2,NN)=IY+1 D(3,NN)=1 ELSE IF(IY.EQ.KATJ)THEN NN=NN+1 D(1,NN)=DIST(2,KATJ,NNJ) D2=DIST(2,KKJ,KATJ) IF(D2.LT.D(1,NN))D(1,NN)=D2 D(2,NN)=IY D(3,NN)=1 ELSE NN=NN+1 D(1,NN)=DIST(2,IY,IY-1) D(2,NN)=IY D(3,NN)=1 NN=NN+1 D(1,NN)=DIST(2,IY+1,IY+2) D(2,NN)=IY+1 D(3,NN)=1 ENDIF C C CHECK WHETHER THE COLLAPSED RATING DATA MATRIX IS DEGENERATE C CALL SORT3(NN,D) ICLAPS=0 DO 1400 II=1,NN IF(ICLAPS.EQ.1)RETURN ID=D(2,II) IF(D(3,II).LT.0)THEN C C CHECK FOR CONDITION X C CALL CMPCT(ID,KATI,FPFX,TPFX,FPF,TPF) CALL CHECK(ICLASS,KATI-1,FPF,TPF,ICON,IARY) IF(ICLASS.EQ.0)THEN C C COLLAPSE THE TARGET CELL C DO 1220 J=1,KATJ VV(J,ID-1)=VV(J,ID-1)+VV(J,ID) WW(J,ID-1)=WW(J,ID-1)+WW(J,ID) DO 1210 I=ID,KATI-1 VV(J,I)=VV(J,I+1) WW(J,I)=WW(J,I+1) 1210 CONTINUE VV(J,KATI)=0. WW(J,KATI)=0. 1220 CONTINUE DO 1230 I=ID-1,KATI-1 ZZX(I)=ZZX(I+1) 1230 CONTINUE ZZX(KATI)=0. KATI=KATI-1 ICLAPS=1 CALL CLAPSE CALL TPFFPF ENDIF ELSE C C CHECK FOR CONDITION Y C CALL CMPCT(ID,KATJ,FPFY,TPFY,FPF,TPF) CALL CHECK(ICLASS,KATJ-1,FPF,TPF,ICON,IARY) IF(ICLASS.EQ.0)THEN C C COLLAPSE THE TARGET CELL C DO 1320 I=1,KATI VV(ID-1,I)=VV(ID-1,I)+VV(ID,I) WW(ID-1,I)=WW(ID-1,I)+WW(ID,I) DO 1310 J=ID,KATJ-1 VV(J,I)=VV(J+1,I) WW(J,I)=WW(J+1,I) 1310 CONTINUE VV(KATJ,I)=0. WW(KATJ,I)=0. 1320 CONTINUE DO 1330 I=ID-1,KATJ-1 ZZY(I)=ZZY(I+1) 1330 CONTINUE ZZX(KATI)=0. KATJ=KATJ-1 ICLAPS=1 CALL CLAPSE CALL TPFFPF ENDIF ENDIF 1400 CONTINUE C C NO COLLAPSE WAS PERFORMED BECAUSE OF DEGENERACY C IU=2 RETURN END C------------------------------------------------------------------- FUNCTION DIST(ICON,N1,N2) C------------------------------------------------------------------- C C CALCULATE DISTANCE BETWEEN TWO CONSECUTIVE OPERATING POINTS C COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) C IF(ICON.EQ.1)THEN DIST=SQRT((FPFX(N1)-FPFX(N2))**2+(TPFX(N1)-TPFX(N2))**2) ELSE DIST=SQRT((FPFY(N1)-FPFY(N2))**2+(TPFY(N1)-TPFY(N2))**2) ENDIF RETURN END C-------------------------------------------- SUBROUTINE SORT3(NUM,F2) C-------------------------------------------- C C REARRANGE A 3-D MATRIX IN ACCENDING ORDER C DIMENSION F2(3,4) 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 TEMP = F2(3,I) F2(3,I)= F2(3,JJ) F2(3,JJ)= TEMP 30 CONTINUE RETURN END C-------------------------------------------------------- SUBROUTINE CMPCT(NN,KAT,FPF1,TPF1,FPF2,TPF2) C-------------------------------------------------------- C C DELETE THE NN-TH OPERATING POINT C DIMENSION FPF1(*),TPF1(*),FPF2(*),TPF2(*) DO 1010 I=1,12 FPF2(I)=0. TPF2(I)=0. 1010 CONTINUE DO 1020 I=2,NN FPF2(I-1)=FPF1(I-1) TPF2(I-1)=TPF1(I-1) 1020 CONTINUE DO 1030 I=NN,KAT FPF2(I)=FPF1(I+1) TPF2(I)=TPF1(I+1) 1030 CONTINUE RETURN END C------------------------------------------------------------------- SUBROUTINE TPFVAL C------------------------------------------------------------------- C C CALCULATE TPF VALUES FOR CONDITION X AND CONDITION Y ON THE FITTED C ROC CURVES, AT SELECTED FPF VALUES. C 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/BLK7/AX,BX,AY,BY,R,RHO,T(21),U(21) DO 9605 I=1,26 CALL ZDEV(FPVAL(I),Q,D,IER) DEVX=AX+BX*Q DEVY=AY+BY*Q CALL NDTR(DEVX,TPVALX(I),D) CALL NDTR(DEVY,TPVALY(I),D) 9605 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE TEST(IERROR) C----------------------------------------------------------------------- C C COMPUTE AREA('A SUB Z') FOR CONDITION X AND CONDITION Y RESPECTIVELY; C ALSO CALCULATE THE VARIANCES AND COVARIANCES OF AREAS. 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/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) COMMON/BLK8/AZX,AZY,VAZX,VAZY,COVAZ,RAREA,XSTAT,PLEVEL, 1 ZSTAT,ZLEVEL,TPSTAT,TPLEVL INTEGER TCODE C IF((KSNX.EQ.1 .AND. KSNY.EQ.1) .OR. (KSNX.EQ.KATI .AND. 1 KSNY.EQ.KATJ))GO TO 9510 R=-R RHO=-RHO C 9510 RDATA=((VCOV(1,2)+VCOV(3,4))-(VCOV(1,4)+VCOV(2,3)))/ 1SQRT((VCOV(1,1)+VCOV(3,3)-2.*VCOV(1,3))*(VCOV(2,2)+VCOV(4,4)-2.* 2VCOV(2,4))) WORKX=BX*BX+1 WORKY=BY*BY+1 ADEVX=AX/SQRT(WORKX) ADEVY=AY/SQRT(WORKY) CALL NDTR(ADEVX,AZX,DENX) CALL NDTR(ADEVY,AZY,DENY) VAZX=(DENX*DENX)*(VCOV(1,1)/WORKX+((AX*BX)**2)*VCOV(2,2)/ 1 (WORKX**3)-2.*AX*BX*VCOV(1,2)/(WORKX*WORKX)) VAZY=(DENY*DENY)*(VCOV(3,3)/WORKY+((AY*BY)**2)*VCOV(4,4)/ 1 (WORKY**3)-2.*AY*BY*VCOV(3,4)/(WORKY*WORKY)) COVAZ=(DENX*DENY)*(VCOV(1,3)/SQRT(WORKX*WORKY)-AY*BY*VCOV(1,4)/ 1 SQRT(WORKX*(WORKY**3))-AX*BX*VCOV(2,3)/SQRT(WORKY*(WORKX**3)) 2 +AX*BX*AY*BY*VCOV(2,4)/SQRT((WORKX*WORKY)**3)) RAREA=COVAZ/SQRT(VAZX*VAZY) IF(IERROR.GT.0)RETURN C C OBTAIN THE RESULTS OF REQUESTED STATISTICAL TEST C IF (TCODE.EQ.2) GO TO 9549 IF (TCODE.EQ.3) GO TO 9550 C C OBTAIN THE BIVARIATE CHI-SQUARE TEST STATISTIC VALUE FOR THE C SIGNIFICANCE OF THE DIFFERENCE BETWEEN THE TWO CORRELATED ROC CURVES C XSTAT=1./(1.-RDATA**2)*((AX-AY)**2/(VCOV(1,1)+VCOV(3,3)-2.* 1VCOV(1,3))+(BX-BY)**2/(VCOV(2,2)+VCOV(4,4)-2.*VCOV(2,4))-2.*RDATA* 2(AX-AY)*(BX-BY)/(SQRT(VCOV(1,1)+VCOV(3,3)-2.*VCOV(1,3))* 3SQRT(VCOV(2,2)+VCOV(4,4)-2.*VCOV(2,4)))) PLEVEL=EXP(-XSTAT/2.) RETURN C C CALCULATE THE AREA Z-SCORE TEST STATISTIC VALUE C 9549 ZSTAT=(AZX-AZY)/SQRT(VAZX+VAZY-2.*COVAZ) CALL NDTR(ZSTAT,P,D) IF (P.LT.0.5) ZLEVEL=2.*P IF (P.GE.0.5) ZLEVEL=2.*(1.-P) RETURN C C COMPUTE THE TRUE-POSITIVE-FRACTION TEST STATISTIC VALUE AT C FALSE-POSITIVE FRACTION SPECIFIED BY USER. C 9550 P=1-FP CALL ZDEV(P,Z,D,IER) XTPFX=Z*BX-AX YTPFY=Z*BY-AY VTPX=VCOV(1,1)+Z*Z*VCOV(2,2)- 1 2.*Z*VCOV(1,2) VTPY=VCOV(3,3)+Z*Z*VCOV(4,4)- 1 2.*Z*VCOV(3,4) COVTP=VCOV(1,3)-Z*VCOV(1,4)-Z*VCOV(2,3)+ 1 Z*Z*VCOV(2,4) TPSTAT=(XTPFX-YTPFY)/SQRT(VTPX+VTPY-2.*COVTP) CALL NDTR(TPSTAT,P,D) IF (P.LT.0.5) TPLEVL=2.*P IF (P.GT.0.5) TPLEVL=2.*(1.-P) RETURN END C------------------------------------------------------------------ SUBROUTINE OUTRSL C------------------------------------------------------------------ C C COMPLETE THE PRINTED OUTPUT PAGE C 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 WRITE(6,9505)ITER,AX,BX,AY,BY 9505 FORMAT(20X,'PROCEDURE CONVERGES AFTER ',I3,' ITERATIONS'// 1 10X,'FINAL ESTIMATES OF THE BINORMAL ROC ', 2 'PARAMETERS AND THE INTER-CONDITION CORRELATION ', 3 'COEFFICIENTS:'/ 4 15X,'BINORMAL ROC PARAMETERS FOR CONDITION X:'/ 5 20X,'A = ',F7.4/20X,'B = ',F7.4/ 6 15X,'BINORMAL ROC PARAMETERS FOR CONDITION Y:'/ 7 20X,'A = ',F7.4/20X,'B = ',F7.4/) WRITE(6,9515)R,RHO 9515 FORMAT(15X,'EFFECTIVE CORRELATION OF TEST RESULTS ', 1 'BETWEEN CONDITIONS:'/ 2 20X,'FOR ACTUALLY NEGATIVE CASES = ',F7.4/ 3 20X,'FOR ACTUALLY POSITIVE CASES = ',F7.4//) SDAX=SQRT(VCOV(1,1)) SDBX=SQRT(VCOV(2,2)) SDAY=SQRT(VCOV(3,3)) SDBY=SQRT(VCOV(4,4)) SDR=SQRT(VCOV(5,5)) SDRHO=SQRT(VCOV(6,6)) FUNLIK=VLOGL(KATI,KATJ,VV,WW,CPELR,CPELP) WRITE(6,9520)SDAX,SDBX,SDAY,SDBY,SDR,SDRHO WRITE(6,9521)VCOV(1,2),VCOV(3,4),VCOV(1,3),VCOV(1,4), + VCOV(2,4),VCOV(2,3) 9520 FORMAT(10X,'ESTIMATED STANDARD DEVIATIONS OF THE BINORMAL ', 1 'ROC PARAMETERS:'/15X,'FOR CONDITION X:'/ 2 20X,'S.D.(A) = ',F7.4/20X,'S.D.(B) = ',F7.4/ 3 15X,'FOR CONDITION Y:'/ 4 20X,'S.D.(A) = ',F7.4/20X,'S.D.(B) = ',F7.4// 5 10X,'ESTIMATED STANDARD DEVIATIONS OF THE INTER-', 6 'CONDITION CORRELATION COEFFICIENTS:'/ 7 20X,'S.D.(R) FOR ACTUALLY NEGATIVE CASES = ',F7.4/ 8 20X,'S.D.(R) FOR ACTUALLY POSITIVE CASES = ',F7.4) 9521 FORMAT(5X,6(F10.8,1X)) RETURN END C-------------------------------------------------------------------- SUBROUTINE INFORM(IERROR) C-------------------------------------------------------------------- C C OUTPUT THE TPF VALUES AT SELECTED FPF VALUES, AND ALSO C AREAS, VARIANCES, COVARIANCE, AND CORRELATION. 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/BLK5/FPVAL(26),TPVALX(26),TPVALY(26) COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(21),U(21) COMMON/BLK8/AZX,AZY,VAZX,VAZY,COVAZ,RAREA,XSTAT,PLEVEL, 1 ZSTAT,ZLEVEL,TPSTAT,TPLEVL COMMON/OP/ZZX(20),ZZY(20) COMMON/DECIML/MAXDX,MAXDY INTEGER TCODE C WRITE(6,9610) 9610 FORMAT(1H1,17X,1H','PLOTS',1H',' OF THE FITTED BINORMAL', 1 ' ROC CURVES:', 1 //12X,'**********************************', 2 '***********************', 3 ///,22X,'FPF',7X,'TPF FOR',8X,'TPF FOR',/,31X,'CONDITION X', 4 4X,'CONDITION Y',/,21X,'-----',4X,'-------------',2X, 5 '-------------',/) DO 9630 I=1,26 WRITE(6,9625) FPVAL(I),TPVALX(I),TPVALY(I) 9625 FORMAT(21X,F5.3,6X,F6.4,9X,F6.4) 9630 CONTINUE SDAZX=SQRT(VAZX) SDAZY=SQRT(VAZY) WRITE(6,9635) AZX,AZY,SDAZX,SDAZY,RAREA 9635 FORMAT(////,12X,'AREA(X)= ',F7.4,5X,'AREA(Y)= ',F7.4,/,12X, 1 'S.D.(AREAX)= ',F7.4,5X,'S.D.(AREAY)= ',F7.4,5X/ 2 12X,'CORRELATION OF AREA(X) AND AREA(Y) =',F7.4) C WRITE(6,9600) 9600 FORMAT(1H1,30X,'FOR CONDITION X:') CALL OPTPFS(AX,BX,KKI,T,ZZX,MAXDX) WRITE(6,9601) 9601 FORMAT(1H1,30X,'FOR CONDITION Y:') CALL OPTPFS(AY,BY,KKJ,U,ZZY,MAXDY) C IF(IERROR.EQ.0)GO TO 9639 IF(IERROR.EQ.2)WRITE(6,9636) 9636 FORMAT(///23X,'*** INVALID INPUT OF STATISTICAL TEST CODE ', 1 '(TCODE), NO STATISTICAL TEST CALCULATED ***') IF(IERROR.EQ.3)WRITE(6,9637) 9637 FORMAT(///23X,'*** INVALID INPUT OF FPF VALUE AT WHICH THE TPF', 1 1H','S ARE TO BE COMPARED;'/ 2 23X,'STATISTICAL TEST ABORTED ***') RETURN C C OUTPUT THE RESULTS OF THE REQUESTED STATISTICAL TEST C 9639 WRITE(6,9640) 9640 FORMAT(/////23X,'STATISTICAL SIGNIFICANCE OF THE DIFFERENCE', 1 ' BETWEEN THE'/ 2 22X,'TWO CORRELATED ROC CURVES ACCORDING TO THE', 3 ' SELECTED TEST:' 4 //12X,'******************************************', 5 '***********************************') IF (TCODE.EQ.2) GO TO 9642 IF (TCODE.EQ.3) GO TO 9648 WRITE(6,9641) XSTAT,PLEVEL 9641 FORMAT(//12X,'THE COMPUTED CORRELATED ',1H','BIVARIATE CHI-', 1 'SQUARE TEST',1H',' STATISTIC VALUE IS',F8.4,/16X, 2 'WITH A CORRESPONDING P-LEVEL OF',F7.4,'.') RETURN 9642 ONEPL=0.5*ZLEVEL WRITE(6,9645) ZSTAT,ZLEVEL,ONEPL 9645 FORMAT(//12X,'THE COMPUTED CORRELATED ',1H','AREA TEST',1H', 1' STATISTIC VALUE IS',2X,F8.4,/16X,'WITH A CORRESPONDING TWO-', 2'TAILED P-LEVEL OF',F7.4,' AND ONE-TAILED P-LEVEL OF',F7.4,'.') RETURN 9648 ONEPL=0.5*TPLEVL WRITE(6,9650) FP,TPSTAT,TPLEVL,ONEPL 9650 FORMAT(//12X,'AT FALSE-POSITIVE-FRACTION EQUAL TO ',F4.2,', ', 1 /,2X,'THE COMPUTED CORRELATED ',1H', 'TRUE-POSITIVE-FRACTION', 2 ' TEST',1H',' STATISTIC VALUE IS',F8.4,/16X,'WITH A CORRES', 3 'PONDING TWO-TAILED P-LEVEL OF',F7.4,' AND ONE-TAILED P-', 4 'LEVEL OF',F7.4,'.') RETURN END C---------------------------------------- SUBROUTINE OPTPFS(A,B,KK,X,ZZ,MAXD) C---------------------------------------- C C COMPUTE THE TPF'S AND FPF'S FOR EACH EXPECTED OPERATING POINT C ESTIMATED ON THE FITTED CURVE. C DIMENSION X(*),ZZ(*),ETPF(20),EFPF(20) DO 3305 I=1,KK Q=X(I+1) CALL NDTR(Q,P,D) EFPF(I)=1-P Q=B*X(I+1)-A CALL NDTR(Q,P,D) ETPF(I)=1-P 3305 CONTINUE WRITE(6,3310) 3310 FORMAT(//20X,'ESTIMATED RELATIONSHIP BETWEEN THE CRITICAL', 1 ' TEST-RESULT VALUE',/,20X,'(WHICH SEPARATES ',1H', 2 'POSITIVE',1H',' RESULTS FROM ',1H','NEGATIVE',1H', 3 ' RESULTS)',/,20X,'AND THE CORRESPONDING OPERATING', 4 ' POINT ON THE FITTED BINORMAL'/20X,'ROC CURVE:'/ 5/15X,'******************************************************', 6 '*****************', 7 //,24X,'CRITICAL TEST',11X,'( FPF , TPF )', 8 /24X,'RESULT VALUE'//) DO 3330 I=1,KK II=KK-I+1 IF(MAXD.EQ.0)WRITE(6,3320)ZZ(II),EFPF(II),ETPF(II) IF(MAXD.EQ.1)WRITE(6,3321)ZZ(II),EFPF(II),ETPF(II) IF(MAXD.EQ.2)WRITE(6,3322)ZZ(II),EFPF(II),ETPF(II) IF(MAXD.EQ.3)WRITE(6,3323)ZZ(II),EFPF(II),ETPF(II) IF(MAXD.EQ.4)WRITE(6,3324)ZZ(II),EFPF(II),ETPF(II) IF(MAXD.EQ.5)WRITE(6,3325)ZZ(II),EFPF(II),ETPF(II) IF(MAXD.EQ.6)WRITE(6,3326)ZZ(II),EFPF(II),ETPF(II) IF(MAXD.EQ.7)WRITE(6,3327)ZZ(II),EFPF(II),ETPF(II) IF(MAXD.GE.8)WRITE(6,3328)ZZ(II),EFPF(II),ETPF(II) 3320 FORMAT(26X,F9.1,13X,'(',F6.3,',',1X,F6.3,')') 3321 FORMAT(25X,F10.2,13X,'(',F6.3,',',1X,F6.3,')') 3322 FORMAT(25X,F11.3,12X,'(',F6.3,',',1X,F6.3,')') 3323 FORMAT(24X,F12.4,12X,'(',F6.3,',',1X,F6.3,')') 3324 FORMAT(24X,F13.5,11X,'(',F6.3,',',1X,F6.3,')') 3325 FORMAT(23X,F14.6,11X,'(',F6.3,',',1X,F6.3,')') 3326 FORMAT(23X,F15.7,10X,'(',F6.3,',',1X,F6.3,')') 3327 FORMAT(22X,F16.8,10X,'(',F6.3,',',1X,F6.3,')') 3328 FORMAT(22X,F17.9,9X,'(',F6.3,',',1X,F6.3,')') 3330 CONTINUE RETURN END C--------------------------------------------- SUBROUTINE M2TOM1(N,XX,VXX) C--------------------------------------------- C C THIS SUBROUTIN STORES THE UPPER TRIANGULAR PART OF C A SYMMETRIC (N BY N) MATRIX XX COLUMNWISE AS A VECTOR C VXX WITH LENGTH N*(N+1)/2. C DIMENSION XX(44,44),VXX(*) K=0 DO 100 J=1,N DO 100 I=1,J K=K+1 VXX(K)=XX(I,J) 100 CONTINUE RETURN END C--------------------------------------------- SUBROUTINE M1TOM2(N,VXX,XX) C--------------------------------------------- C C THIS SUBROUTINE CREATES A SYMETRIC (N BY N) MATRIX FROM C THE VECTOR VXX WITH LENGTH N*(N+1)/2, WHICH REPRESENTS C THE UPPER TRIANGULAR PART OF XX BY COLUMNS. C DIMENSION XX(44,44),VXX(*) K=0 DO 100 J=1,N DO 100 I=1,J K=K+1 XX(I,J)=VXX(K) IF(I.NE.J)XX(J,I)=XX(I,J) 100 CONTINUE RETURN END C------------------------------------------------------ SUBROUTINE OUTINI(AX,BX,AY,BY,R,RHO) C------------------------------------------------------ C C PRINT OUT INITIAL ESTIMATES C WRITE(6,9505)AX,BX,AY,BY 9505 FORMAT(6X,'INITIAL ESTIMATES OF THE BINORMAL ROC ', 2 'PARAMETERS AND'/ 3 6X,'THE INTER-CONDITION CORRELATION COEFFICIENTS:'// 4 11X,'BINORMAL ROC PARAMETERS FOR CONDITION X:'/ 5 16X,'A = ',F7.4/16X,'B = ',F7.4/ 6 11X,'BINORMAL ROC PARAMETERS FOR CONDITION Y:'/ 7 16X,'A = ',F7.4/16X,'B = ',F7.4/) WRITE(6,9515)R,RHO 9515 FORMAT(11X,'EFFECTIVE CORRELATION OF TEST RESULTS ', 1 'BETWEEN CONDITIONS:'/ 2 16X,'FOR ACTUALLY NEGATIVE CASES = ',F7.4/ 3 16X,'FOR ACTUALLY POSITIVE CASES = ',F7.4//) WRITE(6,9525) 9525 FORMAT(///6X,'*** NO FINAL ESTIMATES CAN BE DERIVED FOR THIS ', 1 'DATA SET ***') 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 - CLASS NUMBER OF DEGENERACY C KAT - NUMBER OF CATEGORIES OF DATA SET C FPF, TPF - OBSERVED OPERATING POINTS C ICON,IARY - STORAGE SPACE FOR INFORMATION TO OUPUT MESSAGE, C FOR USE IN SUBROUTINE 'DEGENE' IF NECESSARY. C REAL FPF(*),TPF(*) INTEGER IARY(*) ICLASS=0 NN=KAT+1 ICON=0 DO 10 I=2,KAT T=FPF(I)*TPF(I)*(1.-FPF(I))*(1.-TPF(I)) IF (T.EQ.0.) GO TO 10 ICON=ICON+1 IARY(ICON)=I 10 CONTINUE IF (KAT.EQ.1) GO TO 20 IF (ICON.EQ.0) GO TO 40 I0=IARY(1)-1 I1=IARY(ICON)+1 IF (KAT.EQ.2.AND.ICON.NE.0) GO TO 30 IF (ICON.GT.1) GO TO 15 IF ((TPF(I0).NE.1.0.OR.TPF(I1).NE.0.0).AND.(FPF(I0).EQ.1.AND. 1 FPF(I1).EQ.0)) ICLASS=5 IF ((FPF(I0).NE.1.0.OR.FPF(I1).NE.0.0).AND.(TPF(I0).EQ.1.AND. 1 TPF(I1).EQ.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.(I1.EQ.0.OR. 1 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 ESTM1(A,B,KAT,OFPF,OTPF,X,EMN,EMS) C--------------------------------------------------------------------- C C PURPOSE C OBTAIN INITIAL PARAMETER ESTIMATES ASSOCIATED WITH A MARGINAL C DATA SET, FOR USE IN SUBROUTINE 'DORALF'. C C NOTE: C THE CONTENT OF THIS SUBROUTINE IS TAKEN FROM THE PROGRAM C 'RSCORE II' BY DONALD D. DORFMAN (DEPT. OF PSYCHOLOGY, C THE UNIVERSITY OF IOWA) WITH ONLY MINOR CHANGES. C C C DESCRIPTION OF PARAMETERS C A, B - ESTIMATES OF 'A' AND 'B' PARAMETERS C KAT - NUMBER OF CATEGORIES C FPF, TPF - OBSEVED OERATING POINTS C X - ESTIMATES OF 'CUTOFF' PARAMETERS C EMN - TOTAL NUMBER OF NOISE TRIALS C EMS - TOTAL NUMBER OF S+N TRIALS C REAL OFPF(*),OTPF(*),X(*),XS(19),FPF(21),TPF(21) do 5 i=1,kat+1 fpf(i)=ofpf(i) tpf(i)=otpf(i) 5 continue Q=0. DO 10 I=2,KAT IF (FPF(I).EQ.0.) FPF(I)=1./(2.*EMN) IF (ABS(FPF(I)-1.).LT.1.0E-05) FPF(I)=1-(1./(2.*EMN)) P=FPF(I) CALL ZDEV(P,Q,D,IER) X(I-1)=Q IF (TPF(I).EQ.0.) TPF(I)=1./(2.*EMS) IF (ABS(TPF(I)-1.).LT.1.0E-05) TPF(I)=1.-(1./(2.*EMS)) P=TPF(I) CALL ZDEV(P,Q,D,IER) XS(I-1)=Q 10 CONTINUE IZ=KAT-2 DO 20 I=1,IZ IF (X(IZ-I+1).LE.X(IZ-I+2)) X(IZ-I+1)=X(IZ-I+2)+0.1 IF (XS(IZ-I+1).LE.XS(IZ-I+2)) XS(IZ-I+1)=XS(IZ-I+2)+0.1 20 CONTINUE SUMX=0. SUMY=0. SUMXY=0. SUMX2=0. XK=KAT-1 KK=KAT-1 DO 30 I=1,KK SUMX=SUMX+X(I) SUMY=SUMY+XS(I) SUMXY=SUMXY+XS(I)*X(I) SUMX2=SUMX2+X(I)*X(I) 30 CONTINUE XMEAN=SUMX/XK YMEAN=SUMY/XK B=(XK*SUMXY-SUMX*SUMY)/(XK*SUMX2-SUMX**2) A=YMEAN-B*XMEAN RETURN END C--------------------------------------------------------------------- SUBROUTINE DORALF(A,B,KAT,EMN,EMS,CATN,CATS,X,VA,VB,VAB,IER) C--------------------------------------------------------------------- C C PURPOSE C OBTAIN THE FINAL PARAMETER ESTIMATES ASSOCIATED WITH A MARGINAL C DATA SET, FOR USE AS INITIAL PARAMETER ESTIMATES IN THE MAIN C PROGRAM. C C NOTE: C THE CONTENT OF THIS SUBROUTINE (AND THE SUBROUTINES THAT IT CALLS) C ARE TAKEN FROM THE PROGRAM 'RSCORE II' BY DONALD D. DORFMAN, C DEPT. OF PSYCHOLOGY, THE UNIVERSITY OF IOWA, WITH ONLY MINOR C MODIFICATIONS. C C DESCRIPTION OF PARAMETERS C A,B - ESTIMATES OF 'A' AND 'B' PARAMETERS C KAT - NUMBER OF CATEGORIES C EMN - TOTAL NUMBER OF NOISE TRIALS C EMS - TOTAL NUMBER OF S+N TRIALS C CATN - VECTOR: NUMBER OF NOISE TRIAL RESPONSES IN EACH CATEGORY C CATS - VECTOR: NUMBER OF S+N TRIAL RESPONSES IN EACH CATEGORY C X - ESTIMATES OF 'CUTOFF' PARAMETERS C VA,VB,VAB - ESTIMATES OF VARIANCES AND COVARIANCE OF 'A' AND C 'B' (NOT USED IN 'CORROC') C IER - RESULTING ERROR PARAMETER CODED AS FOLLOWING: C =K 'IER' CODE OBTAINED FROM MATRIX INVERSION C SUBROUTINE C =101 EXCEEDS 100 ITERATIONS. C REAL CATN(*),CATS(*),X(*),XX(21,21),XP(21,21), 1 PX(21),PXBA(21),DX(21),DXBA(21),DIFFSN(21),DIFFN(21), 2 PARX(21),ADJX(21),XXDUM(300) KK=KAT-1 NN=KAT+1 SWT=0 LL=0 C C GET INTEGRALS AND DENSITIES C DO 7045 J=1,NN DO 7045 K=1,NN XX(J,K)=0. 7045 XP(J,K)=0. 7050 LL=LL+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. IF (LL.GT.1) GO TO 7064 DO 7055 I=2,KAT Y=-X(I-1)*B-A CALL NDTR(Y,PXBA(I),DXBA(I)) 7055 CONTINUE DO 7060 I=2,KAT Y=-X(I-1) CALL NDTR(Y,PX(I),DX(I)) 7060 CONTINUE DO 7062 I=1,KAT DIFFSN(I)=PXBA(I+1)-PXBA(I) DIFFN(I)=PX(I+1)-PX(I) IF (DIFFSN(I).LE.5.0E-08) DIFFSN(I)=5.0E-08 IF (DIFFN(I).LE.5.0E-08) DIFFN(I)=5.0E-08 7062 CONTINUE FUNLIK=0. DO 7063 I=1,KAT FUNLIK=FUNLIK+CATN(I)*ALOG(DIFFN(I))+CATS(I)*ALOG(DIFFSN(I)) 7063 CONTINUE C C GET FIRST PARTIALS C WITH RESPECT TO A C 7064 AAA=0. DO 7065 I=2,KAT 7065 AAA=AAA-DXBA(I)*(CATS(I-1)/DIFFSN(I-1)-CATS(I)/DIFFSN(I)) C C WITH RESPECT TO B C BBB=0. DO 7070 I=2,KAT 7070 BBB=BBB-DXBA(I)*X(I-1)*(CATS(I-1)/DIFFSN(I-1)-CATS(I)/ 1DIFFSN(I)) C C WITH RESPECT TO Z'S C DO 7075 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)) 7075 PARX(I-1)=Q1+Q2 C C GET EXPECTED SECOND AND MIXED PARTIALS C WITH RESPECT TO A C XX(1,1)=0. DO 7080 I=2,KAT 7080 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 7095 I=2,KAT D=X(I-1) IF(I.EQ.KAT)GO TO 7090 IF(I.EQ.2)GO TO 7085 DD=X(I-2) DDD=X(I) GO TO 7095 7085 DD=0. DDD=X(I) GO TO 7095 7090 DD=X(I-2) DDD=0. 7095 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 7100 I=2,KAT 7100 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 7105 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))) 7105 XX(I+1,1)=XX(1,I+1) C C WITH RESPECT TO B AND Z'S C DO 7120 I=2,KAT XIL2=0. IF(I.GT.2)XIL2=X(I-2) IF(I.EQ.KAT)GO TO 7110 CHNG=X(I) GO TO 7115 7110 CHNG=0.0 7115 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))) 7120 XX(I+1,2)=XX(2,I+1) C C WITH RESPECT TO Z'S AND MIXED WITH RESPECT TO Z'S C DO 7130 I=2,KAT IF(I.EQ.KAT)GO TO 7125 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) 7125 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))) 7130 CONTINUE C C INVERT MATRIX C DO 7135 I=1,NN DO 7135 J=1,NN 7135 XX(I,J)=-XX(I,J) CALL MSTR(XX,XXDUM,21,0,1) CALL SINV(XXDUM,NN,.001,IER) IF(IER.GE.0) GO TO 7137 ALPHA=-0.5*ABS(ALPHA) GO TO 7157 7137 CALL MSTR(XXDUM,XP,21,1,0) IF(SWT.EQ.1)GO TO 7190 C C FORM SOLUTION VECTOR C ADJA=AAA*XP(1,1)+BBB*XP(1,2) DO 7145 I=1,KK 7145 ADJA=ADJA+PARX(I)*XP(1,I+2) ADJB=AAA*XP(2,1)+BBB*XP(2,2) DO 7150 I=1,KK 7150 ADJB=ADJB+PARX(I)*XP(2,I+2) DO 7155 I=1,KK ADJX(I)=AAA*XP(I+2,1)+BBB*XP(I+2,2) DO 7155 K=3,NN 7155 ADJX(I)=ADJX(I)+PARX(K-2)*XP(I+2,K) C C ITERATE C ALPHA=1.0 7157 A=A+ALPHA*ADJA B=B+ALPHA*ADJB DO 7160 I=1,KK 7160 X(I)=X(I)+ALPHA*ADJX(I) DO 7165 I=2,KAT Y=-X(I-1) CALL NDTR(Y,PX(I),DX(I)) Y=-X(I-1)*B-A CALL NDTR(Y,PXBA(I),DXBA(I)) 7165 CONTINUE DO 7170 I=1,KAT DIFFN(I)=PX(I+1)-PX(I) DIFFSN(I)=PXBA(I+1)-PXBA(I) IF (DIFFN(I).LT.5.0E-8) DIFFN(I)=5.0E-8 IF (DIFFSN(I).LT.5.0E-8) DIFFSN(I)=5.0E-8 7170 CONTINUE IF (ALPHA.NE.1.0) GO TO 7175 C C CHECK FOR MAXIMIZATION C IF (LL.EQ.1) GO TO 7175 TNET=ABS(ADJA)+ABS(ADJB) DO 7173 I=1,KK 7173 TNET=TNET+ABS(ADJX(I)) TNET=6.0*TNET/FLOAT(1+KAT) IF(TNET.GE..001) GO TO 7175 SWT=1 GO TO 7050 7175 IF (ALPHA.EQ.1.0) FUNLI1=FUNLIK FUNLIK=0. DO 7180 I=1,KAT FUNLIK=FUNLIK+CATN(I)*ALOG(DIFFN(I))+CATS(I)*ALOG(DIFFSN(I)) 7180 CONTINUE IF (FUNLIK.GE.FUNLI1) GO TO 7185 FUNDIF=ABS((FUNLIK-FUNLI1)/tnet) IF (alpha.le.1.0e-5 .or. FUNDIF.LE.1.0E-4) GO TO 7185 ALPHA=-0.5*ABS(ALPHA) GO TO 7157 7185 IF(LL.GT.100)GO TO 8430 GO TO 7050 C C ASSIGN THE FINAL VALUES OF A, B, AND THE CUTOFFS OBTAINED FROM C THE MARGINAL DATA SET TO THE RESPECTIVE VARIABLE NAMES THAT C WILL BE USED IN THE MAIN PROGRAM TO ACCOUNT FOR CORRELATIONS C BETWEEN THE MARGINAL DATA SETS. THE VALUES OF THESE NEW VARIABLES C ARE FIRST USED AS INITIAL ESTIMATES IN THE MAIN PROGRAM. C C THE VARIANCES AND COVARIANCE OF THE 'A' AND 'B' PARAMETERS C ESTIMATED FROM THE MARGINAL DATA SET ARE RETURNED ALSO, BUT C ARE TO BE USED FOR COMPARATIVE PURPOSES ONLY --- WITH EACH C OTHER AND WITH THE VARIANCES AND COVARIANCE OF 'A' AND 'B' FOR C THE OTHER MARGINAL DATA SET, WHEN THE TWO SETS OF ROC DATA ARE C TREATED AS CORRELATED IN THE MAIN PROGRAM. C 7190 VA=XP(1,1) VB=XP(2,2) VAB=XP(1,2) RETURN C C THE FOLLOWING TWO SETS OF WRITE AND FORMAT STATEMENTS WILL BE C USED ONLY IF THERE ARE ANY PROBLEMS WITH ITERATION ON THE MARGINAL C DATA SETS IN SUBROUTINE 'DORALF'. C 8430 IER=101 RETURN END C--------------------------------------------------------------------- SUBROUTINE DEGENE(ICLASS,IARY,FPF,TPF) C--------------------------------------------------------------------- C C PURPOSE C OUTPUT DEGENERACY FOR EACH DEGENERACY CLASS. C C DESCRIPTION OF PARAMETERS C ICLASS - CLASS NUMBER OF DEGENERACY C IARY - INFORMATION REGARDING DEGENERATE DATA SET, C FROM SUBROUTINE 'CHECK' C FPF, TPF - OBSERVED OPERATING POINTS C INTEGER IARY(*) REAL FPF(*),TPF(*) 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(*,5) 5 FORMAT(15X,'AND IMPLY PERFECT (BUT PERVERSE) DECISION', 1 /15X,'PERFORMANCE.') RETURN 10 write(*,15) 15 FORMAT(15X,'AND PROVIDE NO OPERATING POINTS OFF THE (0,0)', 1 /15X,'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(*,25) DD 25 FORMAT(15X,'AND PROVIDE ONLY A SINGLE OPERATING POINT OFF', 1 /,15X,'THE (0,0) AND (1,1) CORNERS. BINORMAL ROC CURVE', 2 /15X,'CANNOT BE ESTIMATED UNIQUELY. UNIT SLOPE ROC CURVE', 3 /15X,'(B=1) WOULD HAVE A=',3HD'=,F8.4,'.') 30 write(*,35) TPF(IARY(2)),TPF(IARY(1)) 35 FORMAT(15X,'AND PRODUCE NO OPERATING POINTS OFF THE BORDERS', 1 /15X,'OF THE UNIT SQUARE. IMPLIED EXACT-FIT BINORMAL', 2 /15X,'ROC CURVE IS HORIZONTAL AT UNDETERMINED HEIGHT', 3 /15X,'BETWEEN TPF=',F8.4,'AND TPF=',F8.4,'.') RETURN 40 write(*,45) FPF(IARY(2)),FPF(IARY(1)) 45 FORMAT(15X,'AND PRODUCE NO OPERATING POINTS OFF THE BORDERS', 1 /15X,'OF THE UNIT SQUARE. IMPLIED EXACT-FIT BINORMAL', 2 /15X,'ROC CURVE IS VERTICAL AT UNDETERMINED POSITION', 3 /15X,'BETWEEN FPF=',F8.4,'AND FPF=',F8.4,'.') RETURN 50 write(*,55) TPF(IARY(1)) 55 FORMAT(15X,'IMPLIED EXACT-FIT BINORMAL ROC CURVE IS HORIZONTAL', 1 /,15X,'AT CONSTANT TPF=',F8.4,'.') RETURN 60 write(*,65) FPF(IARY(1)) 65 FORMAT(15X,'IMPLIED EXACT-FIT BINORMAL ROC CURVE IS VERTICAL', 1 /,15X,'AT CONSTANT FPF=',F8.4,'.') RETURN 90 write(*,95) 95 FORMAT(15X,'AND IMPLY PERFECT DECISION 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 MDBNOR (X,Y,RHO,P,IER) C------------------------------------------------------------------ C C FUNCTION - BIVARIATE NORMAL PROBABILITY DISTRIBUTION C FUNCTION C USAGE - CALL MDBNOR(X,Y,RHO,P,IER) C PARAMETERS X - INPUT UPPER LIMIT OF INTEGRATION FOR THE C FIRST VARIABLE C Y - INPUT UPPER LIMIT OF INTEGRATION FOR THE C SECOND VARIABLE C RHO - INPUT CORRELATION COEFFICIENT C P - OUTPUT PROBABILITY THAT THE FIRST VARIABLE C IS LESS THAN OR EQUAL TO X AND THAT THE C SECOND VARIABLE IS LESS THAN OR EQUAL TO Y C IER - ERROR PARAMETER C TERMINAL ERROR = 128+N C N=1 INDICATES THE ABSOLUTE VALUE OF RHO IS C GREATER THAN OR EQUAL TO ONE C PRECISION - SINGLE C REQD. IMSL ROUTINES - MDTNF,UERTST C C C0 = SQRT(.5) DATA C0/.7071068/,C1/.999999/,C2/6.5/,C3/1.193E-7/ EPS = 0.0 IER = 0 IF (ABS(RHO) .LT. C1) GO TO 5 C TERMINAL - RHO OUT OF RANGE IER = 129 GO TO 9000 5 CONTINUE C FOR LARGE X OR Y VALUE - USE GAUSSIAN C APPROXIMATION IF (X.GT.C2 .AND. C2.GT.Y)THEN CALL NDTR(Y,P,D) GO TO 8000 ELSE IF(Y.GT.C2 .AND. C2.GT.X)THEN CALL NDTR(X,P,D) GO TO 8000 ELSE IF(X.GT.C2 .AND. Y.GT.C2)THEN P=1.0 GO TO 9000 ELSE IF(X.LT.-C2 .OR. Y.LT.-C2)THEN P=0.0 GO TO 9000 ENDIF C CALCULATION F1 = 1.0/SQRT(1.0 - RHO**2) XY = X*Y IAX = 0 IAY = 0 IND = 0 IF (XY .EQ. 0.0) GO TO 10 AX = F1*(Y/X - RHO) AY = F1*(X/Y - RHO) GO TO 25 10 IF (X .NE. 0.0) GO TO 15 IF (Y .NE. 0.0) GO TO 20 C 2 1/2 C FOR X=Y=0 AX=AY=(1-RHO)/(1-RHO ) AX = F1*(1.0 - RHO) AY = AX GO TO 25 C FOR Y=0,X LESS THAN 0 TY = -1/4 C FOR Y=0,X GREATER THAN 0 TY = 1/4 15 TY = 0.25 IF (X .LT. 0.0) TY = -TY AX = -F1*RHO IND = 1 GO TO 25 C FOR X=0,Y LESS THAN 0 TX = -1/4 C FOR X=0,Y GREATER THAN 0 TX = 1/4 20 TX = 0.25 IF (Y .LT. 0.0) TX = -TX AY = -F1*RHO GO TO 35 25 IF (AX .GE. 0.0) GO TO 30 IAX = 1 AX = -AX 30 XXX=X MMM=0 33 CALL MDTNF(XXX,AX,EPS,TX,NNN) IF(NNN.GT.100)THEN cc write(6,*)'..............xxx=',xxx IF(XXX.GT.0)THEN XXX=XXX-0.000001 ELSE XXX=XXX+0.000001 ENDIF MMM=MMM+1 IF(MMM.GT.30)THEN IER=1000 RETURN ELSE GO TO 33 ENDIF ENDIF IF (IAX .NE. 0) TX = -TX IF (IND .NE. 0) GO TO 45 35 IF (AY .GE. 0.0) GO TO 40 IAY = 1 AY = -AY 40 YYY=Y MMM=0 43 CALL MDTNF(YYY,AY,EPS,TY,NNN) IF(NNN.GT.100)THEN cc write(6,*)'..............yyy=',yyy IF(YYY.GT.0)THEN YYY=YYY-0.000001 ELSE YYY=YYY+0.000001 ENDIF MMM=MMM+1 IF(MMM.GT.30)THEN IER=1000 RETURN ELSE GO TO 43 ENDIF ENDIF IF (IAY .NE. 0) TY = -TY 45 CALL NDTR(X,QX,DEN) CALL NDTR(Y,QY,DEN) C NOW EVALUATE P P = 0.5*(QX + QY) - TX - TY IF (XY .LE. 0.0 .AND.(XY .NE. 0.0 .OR. X+Y .LT. 0.0)) P = P - 0.5 P = AMIN1(AMAX1(0.0,P),1.0) 8000 CONTINUE IF(P .LT. C3) P=0.0 IF((1.0-P) .LT. C3) P=1.0 9000 CONTINUE IF (IER .NE. 0) CALL UERTST(IER,6HMDBNOR) 9005 RETURN END C-------------------------------------------------------------------- SUBROUTINE MDTNF (Y,Z,EPS,T,nnn) C-------------------------------------------------------------------- C C FUNCTION - INTEGRATE T(Y,Z) FOR NON-CENTRAL T USAGE. C USAGE - CALL MDTNF(Y,Z,EPS,T) C PARAMETERS Y - INPUT PARAMETER. SEE DOCUMENTATION FOR C THE DEFINITION. C Z - INPUT. INTEGRATION IS FROM 0 TO Z. C EPS - INPUT. ACCURACY SHOULD NOT BE LESS THAN EPS. C IF EPS=0 IS ENTERED, EPS=.000001 IS USED. C T - RESULTANT VALUE OF THE INTEGRAL. C PRECISION - SINGLE C DATA C/.1591549/,EXPOV/174.673/ NNN=0 EP1 = EPS IF(EPS .EQ. 0.) EP1 = .000001 T = 0.0 B = ABS(Y) A = ABS(Z) IF(A .EQ. 0.) GO TO 35 5 TA = ATAN(A) IF (A*B .LE. 4.0) GO TO 10 CALL NDTR(B,T,DEN) T = C*(TA+ATAN(1.0/A)) - .5*(T-.5) GO TO 30 C APPROXIMATION FOR SMALL Y*Z 10 HSQB = .5*B*B IF (HSQB .GT. EXPOV) GO TO 35 BEXP = EXP(-HSQB) ASQ = A*A A4 = ASQ*ASQ B4 = HSQB * HSQB A4B4 = A4 * B4 AHSQB = A * HSQB AB4 = A*B4*.5 F = 1.0 SUM = 0.0 G = 3.0 C BEGIN SERIES EXPANSION 15 G1 = G BER = 0.0 TER = AB4 20 BER = BER+TER IF(TER .LE. BER*EP1) GO TO 25 C DEVELOP COEFFICIENT SERIES TER = TER*HSQB/G1 G1 = G1+1.0 GO TO 20 25 D1 = (BER+AHSQB)/F D2 = BER*ASQ/(F+2.0) D = D1-D2 SUM = SUM+D T = TA-SUM*BEXP AEPS = EP1*T AHSQB = AHSQB*A4B4/((G-1.0)*G) AB4 = AB4*A4B4/((G +1.0)*G) F = F+4.0 G = G+2.0 C SHOULD SERIES EXPANSION BE TERMINATED IF (D2*BEXP .GE. AEPS)THEN NNN=NNN+1 IF(NNN.GT.100)RETURN GO TO 15 ENDIF T = T * C 30 IF(Z .LT. 0.0) T = -T 35 RETURN END C---------------------------------------------------------------------- SUBROUTINE UERTST (IER,NAME) C---------------------------------------------------------------------- C C FUNCTION - ERROR MESSAGE GENERATION C USAGE - CALL UERTST(IER,NAME) C PARAMETERS IER - ERROR PARAMETER. TYPE + N WHERE C TYPE= 128 IMPLIES TERMINAL ERROR C 64 IMPLIES WARNING WITH FIX C 32 IMPLIES WARNING C N = ERROR CODE RELEVANT TO CALLING ROUTINE C NAME - INPUT VECTOR CONTAINING THE NAME OF THE C CALLING ROUTINE AS A SIX CHARACTER LITERAL C STRING. C DIMENSION ITYP(5,4),IBIT(4) INTEGER*2 NAME(3) INTEGER WARN,WARF,TERM,PRINTR REAL ITYP EQUIVALENCE (IBIT(1),WARN),(IBIT(2),WARF),(IBIT(3),TERM) DATA ITYP /'WARN','ING ',' ',' ',' ', * 'WARN','ING(','WITH',' FIX',') ', * 'TERM','INAL',' ',' ',' ', * 'NON-','DEFI','NED ',' ',' '/, * IBIT / 32,64,128,0/ DATA PRINTR / 6/ IER2=IER IF (IER2 .GE. WARN) GO TO 5 C NON-DEFINED IER1=4 GO TO 20 5 IF (IER2 .LT. TERM) GO TO 10 C TERMINAL IER1=3 GO TO 20 10 IF (IER2 .LT. WARF) GO TO 15 C WARNING(WITH FIX) IER1=2 GO TO 20 C WARNING 15 IER1=1 C EXTRACT 'N' 20 IER2=IER2-IBIT(IER1) C PRINT ERROR MESSAGE WRITE (PRINTR,25) (ITYP(I,IER1),I=1,5),NAME,IER2,IER 25 FORMAT(' *** I M S L(UERTST) *** ',5A4,4X,3A2,4X,I2, * ' (IER = ',I3,')') 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*M 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 INSERT(NRANK,SUMF,sumz) C----------------------------------------------- C C SELECT THE CANDIDATE FROM THE MAXIMUM (delta FPF+delta TPF) C WHICH HAS MORE THAN ONE POINT IN BETWEEN C real sumz(*) integer*2 mto DIMENSION SUMF(2,1999) max1=0 max2=0 IDELTA=0 NUMCUT=0 LESS=0 mto=2 C MTO = More Than One (in this case two or more points in between) C C COUNT TOTAL NUMBER OF TARGET OPERATING POINTS WHICH HAS C FPF+TPF < 1.0 C DO 1010 I=1,NRANK IF(SUMF(2,I) .LT. 0)GO TO 1010 IF(SUMF(1,I) .LT. 1.0) LESS=LESS+1 NUMCUT=NUMCUT+1 1010 CONTINUE C C IF WE HAVE LESS THAN ONE HALF OF THE TARGET OPERATING POINTS C WITH FPF+TPF < 1 START INSERTING FROM THE SIGNAL END; OTHERWISE, C START FROM THE NOISE END C 1012 IF(LESS .GT. (NUMCUT/2))GO TO 1020 I1=1 JUMP=1 NBOUND=NRANK GO TO 2010 1020 I1=NRANK JUMP=-1 NBOUND=1 C C CHECK FOR TWO CONSECUTIVE TARGET OPERATING POINTS C 2010 IF(SUMF(2,I1) .GT. 0) GO TO 2020 I1=I1+JUMP GO TO 2010 2020 I2=I1 2030 I2=I2+JUMP IF(I2 .EQ. (NBOUND+JUMP))GO TO 2060 IF(SUMF(2,I2) .LT. 0)GO TO 2030 C C ROUND OFF TO 5 SIGNIFICANT DIGITS C DIST=ABS(SUMz(I1)-SUMz(I2))*1.0E+05 IDIST=INT(DIST) IF((DIST-FLOAT(IDIST)) .GE. 0.5)IDIST=IDIST+1 C C SELECT THE MAXIMUM (delta FPF + delta TPF) WHICH ALSO HAS C MORE THAN ONE POINT IN BETWEEN C IF(IDIST .LE. IDELTA) GO TO 2050 if (abs(i2-I1).le.mto) goto 2050 IDELTA=IDIST MAX1=I1 MAX2=I2 2050 I1=I2 GO TO 2030 C C CHOOSE THE MOST MIDDLE POINT AS THE NEW TARGET OPERATING POINT C 2060 if (mto.eq.2.and.idelta.eq.0) then mto=1 goto 1012 endif if (idelta.eq.0) return J1=MAX1+JUMP J2=MAX2-JUMP IF(J1 .LE. J2) GO TO 2070 J1=MAX2-JUMP J2=MAX1+JUMP 2070 DELTA=ABS(SUMz(MAX1)-SUMz(MAX2)) minn=int((j1+j2)/2) DO 2080 J=J1,J2 DIFF=ABS(ABS(SUMz(J)-SUMz(MAX1)) - + ABS(SUMz(J)-SUMz(MAX2))) IF(DIFF .GE. DELTA)GO TO 2080 DELTA=DIFF MINN=J 2080 CONTINUE SUMF(2,MINN)=1 RETURN END C---------------------------------------------------- SUBROUTINE CATGRZ(KOND,KEY,ACTN,ACTSN) C---------------------------------------------------- C C CATEGORIZE QUANTITATIVE DATA INTO CATEGORICAL 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) COMMON/OP/ZZX(20),ZZY(20) COMMON/CAT/V(2,2000),SUMF(2,1999),Z(1999) DIMENSION ACTN(*),ACTSN(*),ZZ(20),fpf(21),tpf(21),iary(21) integer recat recat=0 MN=EMN MS=EMS NUM=MN+MS C C IF SMALLER VALUES SHOW STRONGER EVIDENCE OF SIGNAL, CONVERT C THE SIGN BY MULTIPLYING -1 C IF(KEY.GT.0)GO TO 1030 DO 1010 I=1,MN ACTN(I)=-ACTN(I) 1010 CONTINUE DO 1020 I=1,MS ACTSN(I)=-ACTSN(I) 1020 CONTINUE C C POOL ALL THE SIGNAL-PLUS-NOISE AND NOISE QUANTITATIVE DATA C TOGETHER C 1030 DO 1050 I=1,NUM IF(I.GT.MN)GO TO 1040 V(1,I)=ACTN(I) V(2,I)=-1 GO TO 1050 1040 V(1,I)=ACTSN(I-MN) V(2,I)=1 1050 CONTINUE C C CHECK HOW MANY OPERATING POINTS AVAILABLE FOR THE QUANTITATIVE C DATA C CALL SORT(NUM,V) 1051 CALL RANKS(NRANK,recat) C C IF "NRANK" CUTOFFS ARE AVAILABLE, AND C (A) NRANK <= 9 :: USE ALL OF THE OPERATING POINTS AVAILABLE C (B) 9 < NRANK <= 19 & MIN(MN,MS) < 100 :: COLLAPSE INTO 10 C CATEGORY C (C) 9 < NRANK <= 19 & MIN(MN,MS) >= 100 :: USE ALL OF THE C AVAILABLE OPERATING POINTS C (D) NRANK > 19 & MIN(MN,MS) < 100 :: SAME AS (B) C (E) NRANK > 19 & MIN(MN,MS) >= 100 :: COLLAPSE INTO 20 CATEGORY C IF(NRANK .LT. 10) GO TO 1100 IF(NRANK .GE. 20) GO TO 1110 IF(MN.LT.100 .OR. MS.LT.100)GO TO 1110 1100 DO 1105 I=1,NRANK SUMF(2,I)=1 1105 CONTINUE GO TO 1120 1110 KAT=10 IF(MN.GE.100 .AND. MS.GE.100)KAT=20 CALL SETCUT(NRANK,KAT-1) C C USE ALL THE OPERATING POINTS AVAILABLE C 1120 II=0 DO 1140 J=1,NRANK IF(SUMF(2,NRANK-J+1) .LT. 0)GO TO 1140 II=II+1 ZZ(II)=Z(NRANK-J+1) 1140 CONTINUE KAT=II+1 do il=1,21 fpf(il)=0 tpf(il)=0 enddo do ik=1,Mn do ij=0,kat-2 if (actn(ik).gt.zz(II-ij)) then fpf(II-ij+1)=fpf(II-ij+1)+1 endif enddo fpf(1)=fpf(1)+1 enddo do ik=1,Ms do ij=0,kat-2 if (actsn(ik).gt.zz(II-ij)) then tpf(II-ij+1)=tpf(II-ij+1)+1 endif enddo tpf(1)=tpf(1)+1 enddo do il=1,kat fpf(il)=fpf(il)/Mn tpf(il)=tpf(il)/Ms enddo CALL CHECK(ICLASS,KAT,FPF,TPF,ICON,IARY) if (iclass.ne.0.and.recat.ne.1) then write(*,*) 'Data Degeneracy... Re-Categorizing' recat=1 go to 1051 endif C C CONVERT QUANTITATIVE VALUE INTO CATEGORICAL DATA C CALL CONVRT(MN,ACTN,KAT-1,ZZ) CALL CONVRT(MS,ACTSN,KAT-1,ZZ) C C ASSIGN VARIABLES C IF(KOND.EQ.2)GO TO 1170 KATI=KAT KSNX=KATI IF(KEY.LT.0)KSNX=1 DO 1160 I=2,KATI IF(KEY.LT.0)GO TO 1150 ZZX(I-1)=ZZ(I-1) GO TO 1160 1150 ZZX(I-1)=-ZZ(I-1) 1160 CONTINUE RETURN C 1170 KATJ=KAT KSNY=KATJ IF(KEY.LT.0)KSNY=1 DO 1190 I=2,KATJ IF(KEY.LT.0)GO TO 1180 ZZY(I-1)=ZZ(I-1) GO TO 1190 1180 ZZY(I-1)=-ZZ(I-1) 1190 CONTINUE RETURN END C------------------------------------- SUBROUTINE RANKS(NRANK,recat) C------------------------------------- C C COUNT TOTAL NUMBER OF OPERATING POINTS FOR THE QUANTITATIVE C DATA; COUNT ONE ONLY FOR ANY TIE VALUES C COMMON/BLK2/EMN,EMS,SUMNT(21),SUMST(21),SUMNU(21),SUMSU(21), 1 FPFX(21),TPFX(21),FPFY(21),TPFY(21) COMMON/CAT/V(2,2000),SUMF(2,1999),Z(1999) integer ofpf,otpf,recat double precision rmn,rms C C INITIALIZATION C DO 1005 I=1,1999 SUMF(1,I)=0. SUMF(2,I)=-1 1005 CONTINUE NUM=eMN+eMS chng=0 rmn=0 rms=0 ofpf=0 otpf=2 NRANK=1 if (recat.eq.1) then FPF=0. TPF=0. else fpf=1.0E-7 tpf=1.0E-7 endif DO 1030 I=2,NUM IF(V(2,NUM-I+2) .GT. 0)GO TO 1010 rmn=rmn+1.0 FPF=rmn/eMN ofpf=-1 GO TO 1020 1010 rms=rms+1.0 TPF=rms/eMS ofpf=1 1020 IF(V(1,NUM-I+1).EQ.V(1,NUM-I+2)) then IF(V(2,NUM-I+1).ne.V(2,NUM-I+2)) then chng=1 otpf=2 endif GO TO 1030 endif if (chng.eq.1) ofpf=0 if (fpf.gt.0.and.tpf.gt.0.and.ofpf.ne.otpf) then nrank=nrank+1 chng=0 endif if (recat.eq.1) otpf=ofpf SUMF(1,NRANK)=FPF+TPF Z(NRANK)=(V(1,NUM-I+1)+V(1,NUM-I+2))/2. if (recat.eq.1.and. + (tpf.gt.0.9999.or.fpf.gt.0.9999)) go to 2000 1030 CONTINUE 2000 RETURN END C----------------------------------------- SUBROUTINE SETCUT(NRANK,KK) C----------------------------------------- C C COLLAPSE QUANTITATIVE DATA INTO (KK+1) CATEGORY C COMMON/CAT/V(2,2000),SUMF(2,1999),Z(1999) real h,d,div,lastrank,sumz(1999) C C COUNT TOTAL NUMBER OF DISTINCT TARGET OPERATING POINTS C NUMCUT=0 lastrank=0 J=1 do i=1,nrank c call ndtr((1.5*(sumf(1,i)- c + (sumf(1,1)+sumf(1,nrank))/2.0)),sumz(i),den) call ndtr((1.5*(sumf(1,i)- + 1.0)),sumz(i),den) end do if (sumf(1,1).lt.0.15) then h=0 else h=sumz(1) endif if (sumf(1,nrank).gt.1.85) then d=1 else d=sumz(nrank) endif div=KK if (h.ne.0 .and. d.ne.1) div=div-1 if (h.eq.0 .and. d.eq.1) div=div+1 step = (d-h)/div C Note Div here now represents how many points to grab C either starting from this location or the next div=kk-1 istart=0 if (h.eq.0) then istart=1 div=kk endif DO 1030 I=istart,div TARGET=h+I*step 1010 IF(ABS(SUMz(J+1)-TARGET) .GE. ABS(SUMz(J)-TARGET))then if ((sumz(j)-lastrank).ge.(step/8)) goto 1020 endif J=J+1 if (j.gt.nrank) goto 1031 GO TO 1010 1020 IF(SUMF(2,J) .GT. 0)GO TO 1030 SUMF(2,J)=1 lastrank=sumz(j) NUMCUT=NUMCUT+1 1030 CONTINUE C C IF WE ALREADY HAVE ENOUGH OPERATING POINTS, THEN RETURN; C OTHERWISE, INSERT AS MANY OPERATING POINTS AS WE NEED (= KK-NUMCUT), C IF WE HAVE ENOUGH CANDIDATES OR WE USE UP ALL OF THE C AVAILABLE CANDIDATES, WHICHEVER COMES FIRST. C 1031 NN=KK-NUMCUT IF(NUMCUT .EQ. KK) GO TO 1050 IF((NRANK-NUMCUT) .LT. NN)NN=NRANK-NUMCUT DO 1040 I=1,NN CALL INSERT(NRANK,SUMF,sumz) 1040 CONTINUE C 1050 KK=NUMCUT+NN RETURN END C C THE EXPECTED VALUES OF ALL OTHER MINUS SECOND-ORDER PARTIAL C DERIVATIVES ARE ZERO. THESE ARE THE MIXED PARTIALS WITH C RESPECT TO R AND EACH OF AX,BX,AY,BY AND RHO; THE MIXED PARTIALS C WITH RESPECT TO T(I) AND T(I') WHERE ABS(I-I')>1; THE MIXED PARTIALS C WITH RESPECT TO U(J) AND U(J') WHERE ABS(J-J')>1; AND THE MIXED C PARTIALS WITH RESPECT TO R AND RHO. THESE HAVE ALREADY BEEN SET C TO ZERO WHEN ALL MATRIX ELEMENTS WERE INITIALIZED TO ZERO IN THE C DO-LOOP. C C--------------------------------------------------------------------- SUBROUTINE ITERAT(IU) C--------------------------------------------------------------------- C C NOW THAT THE VALUES OF ALL MATRIX ELEMENTS HAVE BEEN COMPUTED, C THE NEXT STEP IS TO INVERT THE MATRIX. WHEN THE ITERATIVE C PROCEDURE HAS CONVERGED, THE INVERTED MATRIX OF MINUS SECOND-ORDER C PARTIAL DERIVATIVES WILL BE EQUIVALENT TO THE VARIANCE-COVARIANCE C MATRIX. 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 C C crit=1 c write(*,*) 'iteration:',iter ALPHA=1.0 BNDRY=1.0E-06 BNDRY2=5.0E-08 FUNLI1=FUNLIK CALL M2TOM1(44,SOPNEG,XXDUM) CALL SINV(XXDUM,MTRX,.001,IER) IF(IER.LT.0)THEN write(*,*) ' Radicand is non-positive, attempting to correct' ALPHA=-0.5*ABS(ALPHA) IF(ABS(ALPHA).GT.BNDRY)GO TO 1450 CALL SELCEL(I,J) CALL TRUNCT(I,J,IU) RETURN ENDIF IF(IER.GT.0)WRITE(*,1420)IER,ITER 1420 FORMAT(1X,'LOSS OF SIGNIFICANCE. STEP ',I5,'+1 --- ON THE ', 1I3,'TH ITERATION'///) CALL M1TOM2(44,XXDUM,VCOV) IF(LSTAT.EQ.1)LSTAT=2 IF(LSTAT.EQ.2)RETURN C C FORM SOLUTION VECTOR C DO 1440 K=1,MTRX ESTCOR(K)=0. DO 1440 L=1,MTRX ESTCOR(K)=ESTCOR(K)+VCOV(K,L)*FOP(L) 1440 CONTINUE C C COMPUTE PARAMETER VALUES FOR NEXT STEP C 1450 AX=AX+ALPHA*ESTCOR(1) BX=BX+ALPHA*ESTCOR(2) AY=AY+ALPHA*ESTCOR(3) BY=BY+ALPHA*ESTCOR(4) R=R+ALPHA*ESTCOR(5) RHO=RHO+ALPHA*ESTCOR(6) DO 1460 I=1,KKI T(I+1)=T(I+1)+ALPHA*ESTCOR(I+6) 1460 CONTINUE DO 1470 J=1,KKJ U(J+1)=U(J+1)+ALPHA*ESTCOR(J+KKI+6) 1470 CONTINUE C this is non-standard code that supresses line feed C this makes pretty dots waltz across the screen C write(*,'(A \)') '.' C C CHECK FOR STOPPING CRITERION C IF (ALPHA.EQ.1.0 .AND. ITER.GT.1)THEN CRIT=0. DO 1475 I=1,MTRX CRIT=CRIT+ABS(ESTCOR(I)) 1475 CONTINUE CRIT=14.0*CRIT/FLOAT(MTRX) IF(CRIT.LT..001)THEN LSTAT=1 RETURN ENDIF ENDIF C C COMPUTE TERMS USED IN CALCULATION OF LIKELIHOOD FUNCTION C DO 1480 I=2,KATI TX(I-1)=T(I)*BX-AX CALL NDTR(TX(I-1),ELP(I,NNJ),D) CALL NDTR(T(I),ELR(I,NNJ),D) 1480 CONTINUE DO 1490 J=2,KATJ UY(J-1)=U(J)*BY-AY CALL NDTR(UY(J-1),ELP(NNI,J),D) CALL NDTR(U(J),ELR(NNI,J),D) 1490 CONTINUE RHO1=RHO R1=R IF (RHO1.GT.0.9999) RHO1=0.9999 IF (RHO1.LT.-0.9999) RHO1=-0.9999 IF (R1.GT.0.9999) R1=0.9999 IF (R1.LT.-0.9999) R1=-0.9999 DO 1500 I=2,KATI DO 1500 J=2,KATJ CALL MDBNOR(TX(I-1),UY(J-1),RHO1,ELP(I,J),IER) CALL MDBNOR(T(I),U(J),R1,ELR(I,J),IER) 1500 CONTINUE DO 1510 I=1,KATI DO 1510 J=1,KATJ IOP=NNI-I JOP=NNJ-J CPELR(IOP,JOP)=ELR(IOP+1,JOP+1)+ELR(IOP,JOP)-ELR(IOP+1,JOP)- + ELR(IOP,JOP+1) CPELP(IOP,JOP)=ELP(IOP+1,JOP+1)+ELP(IOP,JOP)-ELP(IOP+1,JOP)- + ELP(IOP,JOP+1) IF((CPELR(IOP,JOP).LT.BNDRY2 .AND. VV(JOP,IOP).NE.0) .OR. + (CPELP(IOP,JOP).LT.BNDRY2 .AND. WW(JOP,IOP).NE.0))THEN ALPHA=-0.5*ABS(ALPHA) C write(*,*) C +' Correcting CP:',CPELP(IOP,JOP),CPELR(IOP,JOP),VV(JOP,IOP) IF(ABS(ALPHA).LT.(BNDRY*2.0))THEN CALL TRUNCT(IOP,JOP,IU) RETURN ENDIF go to 1450 ENDIF c IF(CPELR(IOP,JOP).LE.0.0)CPELR(IOP,JOP)=1.0e-7 c IF(CPELP(IOP,JOP).LE.0.0)CPELP(IOP,JOP)=1.0e-7 IF(CPELR(IOP,JOP).LT.BNDRY2)CPELR(IOP,JOP)=BNDRY2/100 IF(CPELP(IOP,JOP).LT.BNDRY2)CPELP(IOP,JOP)=BNDRY2/100 1510 CONTINUE C C CHECK IF LIKELIHOOD FUNCTION VALUES ARE INCREASING, OR THE C DIFFERENCE BETWEEN LIKELIHOOD FUNCTION VALUES ARE QUITE C SMALL; IF NOT, HALVE THE CORRECTION VECTOR. C FUNLIK=VLOGL(KATI,KATJ,VV,WW,CPELR,CPELP) c write(*,11)ax,bx,ay,by,iter+crit*alpha, c + (funlik-funli1)/crit,funlik,funli1 c11 format(1x,4(f8.6,1x),2(f8.5,1x),2(f10.4,1x)) C changed 11/20/93: Make the interval halving routine depend on C the step size rather than the height of the liklihood function C (be fare I mean if you get close your close) IF (FUNLIK.GE.FUNLI1)RETURN FUNDIF=ABS((FUNLIK-FUNLI1)/crit) c write(*,*) funlik,' ',alpha,' ',fundif,' ',crit C Don't check Alpha here like we do in Labroc1 \ C here we use the bndry condition to check alpha IF (abs(alpha).le.BNDRY .or. FUNDIF.LE.1.0E-4)RETURN ALPHA=-0.5*ABS(ALPHA) go to 1450 c IF(ABS(ALPHA).GT.BNDRY)GO TO 1450 c CALL SELCEL(I,J) c CALL TRUNCT(I,J,IU) RETURN END