C C :------------------------: C : : C : 'CORROC2' PROGRAM : Jan 1994 C : : C :------------------------: C C C******************************************************************************* C * C THIS PROGRAM IS APPLICABLE TO CORRELATED, INHERENTLY CATEGORICAL * C RATING-SCALE DATA. FOR CORRELATED BUT CONTINUOUSLY-DISTRIBUTED DATA, * C THE 'CLABROC' PROGRAM SHOULD BE USED INSTEAD. * C * C THE PURPOSES OF THIS PROGRAM ARE: * C (I) TO CALCULATE MAXIMUM LIKELIHOOD ESTIMATES OF THE PARAMETERS * C OF A MODEL FOR CORRELATED ROC RATING-SCALE DATA BASED ON * C AN EFFECTIVE PAIR OF UNDERLYING BIVARIATE-NORMAL DECISION- * C VARIABLE DISTRIBUTIONS; * 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 ROC CURVE -- THAT IS AX=AY AND BX=BY); * 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 ROC CURVES WITH * C 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 ROC CURVES * C HAVING THE SAME TPF AT THE SELECTED FPF). * C * C THE REQUIRED DATA ARE TWO SETS OF CASES -- A SET OF ACTUALLY * C NEGATIVE CASES AND A SET OF ACTUALLY POSITIVE CASES -- WITH A PAIR * C OF CONFIDENCE RATINGS FOR EACH CASE. THE TWO RATINGS ASSOCIATED WITH * C EACH CASE COULD ARISE, FOR EXAMPLE, FROM A SINGLE OBSERVER IN RESPONSE * C TO A STIMULUS PROCESSED OR PRESENTED UNDER TWO CONDITIONS, FROM TWO * C OBSERVERS IN RESPONSE TO THE SAME STIMULUS, OR FROM OBSERVATIONS OF * C TWO DIFFERENT DIAGNOSTIC TESTS PERFORMED ON THE SAME PATIENT. * C * C *** NOTE *** THE CASES USED AS INPUT MUST BE DISTINCT, BECAUSE * C THIS PROGRAM ASSUMES THAT EACH PAIR OF CONFIDENCE- * C RATING DATA AROSE INDEPENDENTLY FROM ONE OF TWO * C MULTINOMIAL DISTRIBUTIONS (CORRESPONDING TO 'ACTUALLY * C NEGATIVE' AND 'ACTUALLY POSITIVE' CASES). HENCE, THE INPUT * C DATA MUST NOT INCLUDE MORE THAN ONE PAIR OF RATINGS FROM * C EACH CASE. IF MULTIPLE CONFIDENCE-RATING PAIRS * C FROM EACH CASE ARE POOLED IN THE INPUT, THE PROGRAM * C WILL OVERESTIMATE THE STATISTICAL SIGNIFICANCE * C OF ANY APPARENT DIFFERENCE BETWEEN THE ROC CURVES, * C THEREBY INVALIDATING THE STATISTICAL TEST. WHEN * C MULTIPLE CONFIDENCE-RATING PAIRS ARE AVAILABLE FROM * C EACH CASE (FOR EXAMPLE, FROM MULTIPLE OBSERVERS WHO * C READ IN BOTH OF TWO OBSERVATION CONDITIONS), THE * C DATA SETS SHOULD BE RUN SEPARATELY (FOR EXAMPLE, FOR * C EACH OBSERVER), OR THE GENERAL APPROACH DESCRIBED BY * C SWETS AND PICKETT ('EVALUATION OF DIAGNOSTIC SYSTEMS: * C METHODS FROM SIGNAL DETECTION THEORY.' ACADEMIC PRESS, * C NEW YORK, 1982, CHAPTER 4) SHOULD BE EMPLOYED. * C * C THIS 'CORROC2' PROGRAM DIFFERS FROM THE EARLIER 'CORROC' PROGRAM * C (WHICH IT REPLACES) IN THAT 'CORROC2' AUTOMATICALLY CREATES THE TWO * C DATA MATRICES THAT ARE REQUIRED BY THE MAXIMUM LIKELIHOOD * C ESTIMATION ALGORITHM. THOSE MATRICES HAD TO BE TABULATED BY THE * C USER OF 'CORROC'. * C * C * C CORROC2 (LIKE CORROC) FIRST TREATS THE MARGINAL DATA SETS INDEPENDENTLY * C AND USES A MODIFIED DORFMAN PROGRAM (J. MATH. PSYCH., VOL.6, 1969, * C PP.487-496, AND VOL.9, 1972, PP. 128-139) TO OBTAIN MAXIMUM * C LIKELIHOOD ESTIMATES OF THE 'A', 'B', AND CATEGORY-BOUNDARY PARAMETERS * C SEPARATELY FOR EACH ROC CURVE. TOGETHER WITH PRODUCT-MOMENT CORRELATION * C COEFFICIENTS CALCULATED DIRECTLY FROM THE TWO DATA MATRICES, THESE * 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 THE BIVARIATE-BINORMAL MODEL THAT IS ASSUMED TO * C UNDERLIE THE CORRELATED RATING DATA. * 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. * C * C THE 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 * C AND HELEN B. KRONMAN, 1979. * C PROGRAM WRITTEN BY HELEN B. KRONMAN, APRIL 1980. * C SUBSEQUENT REVISIONS BY PU-LAN WANG AND JONG-HER SHEN. * 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. A CODE TO INDICATE WHETHER THE ENTIRE INPUT SHOULD BE * C PRINTED. IF THIS 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: EXAMPLES OF THE TWO 'CONDITIONS' ARE GIVEN IN REMARK (1) * C BELOW.) * C LINE #2 AND #4 : A FREE-TEXT DESCRIPTION OF THE CONDITION * C (UP TO 80 CHARACTERS). * C LINE #3 AND #5 : THE NUMBER OF CATEGORIES FOR THE RATING DATA * C FROM THAT CONDITION, AND THE CATEGORY * C REPRESENTING THE STRONGEST EVIDENCE OF * C POSITIVITY (E.G., "ABNORMALITY" OR "SIGNAL * C PRESENCE"). (SEE REMARKS (2) AND (3) BELOW.) * C FORMAT : FREE INPUT FORMAT, I.E., INPUT DATA CAN START * C IN ANY COLUMN BETWEEN COLUMN 1 AND COLUMN 80, * C BUT MULTIPLE ENTRIES ON A SINGLE LINE MUST BE * C SEPARATED BY AT LEAST ONE SPACE. * C 3. RATING-DATA PAIRS FOR ACTUALLY NEGATIVE CASES AND THEN FOR * C ACTUALLY POSITIVE CASES. * C FORMAT : ON A LINE FOR EACH PAIR, ENTER THE X-CONDITION RATING, * C ONE OR MORE BLANK SPACES, THE Y-CONDITION RATING, * C ONE OR MORE BLANK SPACES, AND THEN ANY FREE-TEXT * C DESCRIPTION OF THE PAIR DESIRED. ENTER ALL PAIRS FOR * C ACTUALLY NEGATIVE CASES FIRST AND THEN ALL PAIRS FOR * C ACTUALLY POSITIVE CASES. EACH OF THESE TWO INPUT * C SEQUENCES MUST BE TERMINATED BY A FREE-TEXT LINE * C SUCH THAT THE FIRST CHARACTER IS AN ASTERISK (*). * 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 WILL READ ONLY * C THE FIRST CHARACTER OF THE INPUT 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 TYPICAL EXAMPLE OF 5-CATEGORY DATA--X * C 5 5 * C TYPICAL EXAMPLE OF 5-CATEGORY DATA--Y * C 5 5 * C 3 1 ACTUALLY NEGATIVE CASE #1 * C 2 1 ACTUALLY NEGATIVE CASE #2 * C 3 5 ACTUALLY NEGATIVE CASE #3 * C ... (AND SO ON) * C *END OF ACTUALLY NEGATIVE CASES * C 5 2 ACTUALLY POSITIVE CASE #1 * C 1 5 ACTUALLY POSITIVE CASE #2 * C 3 4 ACTUALLY POSITIVE CASE #3 * C ... (AND SO ON) * C *END OF ACTUALLY POSITIVE CASES * C TPF * C 0.10 * C --------------------------------------------------------- * C * C REMARKS CONCERNING INPUT: * C 1. THE TWO 'CONDITIONS' (X AND Y) CAN REPRESENT DIFFERENT * C READINGS OF THE SAME STIMULI (E.G., IMAGES) BY A SINGLE * C OBSERVER (PERHAPS UNDER DIFFERENT READING, DISPLAY, * C OR DATA-MANIPULATION CONDITIONS); DIFFERENT * C OBSERVERS WHO READ THE SAME SET OF IMAGES; READINGS * C OF DIFFERENT IMAGES OF THE SAME SET OF PATIENTS; ETC. * C 2. THE NUMBERS OF CATEGORIES (FIRST INPUT VALUES ON LINE #3 * C AND LINE #5) CAN BE 3 TO 11 AND NEED NOT BE THE SAME FOR * C THE TWO CONDITIONS. * C 3. THE SECOND INPUT VALUES ON LINE #3 AND LINE #5 INDICATE * C WHETHER THE HIGHEST OR LOWEST NUMBERED CATEGORY REPRESENTS * C THE STRONGEST EVIDENCE OF POSITIVITY (E.G., 'ABNORMALITY') * C IN THE CORRESPONDING CONDITION. IF HIGHER-NUMBERED CATEGORIES * C REPRESENT STRONGER EVIDENCE OF POSITIVITY (E.G., * C 'ABNORMALITY'), THEN THIS VALUE SHOULD BE EQUAL TO * C THE TOTAL NUMBER OF CATEGORIES INPUT FIRST ON THE LINE. * C ALTERNATIVELY, IF LOWER-NUMBERED CATEGORIES REPRESENT * C STRONGER EVIDENCE OF POSITIVITY, THEN THIS VALUE SHOULD * C BE '1'. * C 4. THIS PROGRAM ALLOWS THE USE OF SEVERAL DATA SETS AS INPUT; * C SIMPLY ENTER DATA SETS SEQUENTIALLY. * C * 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 RATING SCALE * C USED FOR THAT CONDITION. * C 3. RATING-DATA MATRICES FOR THE ACTUALLY NEGATIVE AND ACTUALLY * C POSITIVE CASES, RESPECTIVELY. DIFFERENT COLUMNS IN EACH DATA * C MATRIX INDICATE DIFFERENT RATINGS FROM CONDITION X, WHILE * C DIFFERENT ROWS INDICATE DIFFERENT RATINGS FROM CONDITION Y. * C IN EITHER MATRIX, EACH ELEMENT REPRESENTS THE NUMBER OF PAIRED * C READINGS (E.G., FROM A COMMON IMAGE OR A COMMON PATIENT) * C THAT PRODUCED THE PAIR OF RATINGS INDICATED BY THE * C ELEMENT'S POSITION. * C 4. OBSERVED OPERATING POINTS FOR THE TWO CONDITIONS, RESPECTIVELY. * C 5. INITIAL ESTIMATES OF THE PARAMETER VALUES, WITH THE * C LOG-LIKELIHOOD OF THE DATA. * C FOR EACH CONDITION, 'A' REPRESENTS THE 'Y-INTERCEPT' * C AND 'B' REPRESENTS THE SLOPE OF THAT CONDITION'S ROC CURVE WHEN * C IT IS PLOTTED ON NORMAL-DEVIATE AXES. R(NEGATIVE CASES) AND * C R(POSITIVE CASES) REPRESENT THE CORRELATION COEFFICIENTS OF * C THE EFFECTIVE UNDERLYING BIVARIATE-NORMAL DECISION-VARIABLE * C DISTRIBUTIONS FOR ACTUALLY NEGATIVE AND ACTUALLY POSITIVE * C CASES, RESPECTIVELY. T(I) AND U(J) REPRESENT THE CUT-OFFS ON * C THE DECISION AXES FOR CONDITION X AND CONDITION Y, RESPECTIVELY, * C EXPRESSED IN STANDARD DEVIATIONS OF THE EFFECTIVE ACTUALLY * C NEGATIVE DISTRIBUTION; THEY ARE LISTED FROM LEFT TO RIGHT AND * C FROM LOWER TO HIGHER ON THE TWO DECISION-VARIABLE AXES, * C RESPECTIVELY, WITH THE MEAN OF THE ACTUALLY POSITIVE * C DISTRIBUTION ASSUMED TO BE ABOVE AND TO THE RIGHT OF THE MEAN * C OF THE ACTUALLY NEGATIVE DISTRIBUTION. * C 6. THE NUMBER OF ITERATIONS EMPLOYED. * C 7. FINAL ESTIMATES OF THE PARAMETER VALUES, WITH * C THE CORRESPONDING LOG-LIKELIHOOD OF THE DATA. * C 8. VARIANCE-COVARIANCE AND CORRELATION MATRICES FOR THE FINAL * C PARAMETER ESTIMATES. * C 9. TPF VALUES FOR CONDITION X AND CONDITION Y ON THE FITTED ROC * C CURVES, AT SELECTED FPF VALUES. * C 10. THE AREAS ('A SUB Z') UNDER THE FITTED ROC CURVES FOR THE TWO * C CONDITIONS, WITH ESTIMATED STANDARD DEVIATIONS, AND CORRELATION. * C 11. ONE OF THREE TEST-STATISTIC VALUES, WITH ASSOCIATED * C SIGNIFICANCE 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 OUTPUT: * C 1. IT IS CRUCIALLY IMPORTANT TO NOTE THAT THIS PROGRAM * C IS NOT -- REPEAT, IS NOT -- APPROPRIATE FOR POOLED * C CONFIDENCE-RATING PAIRS (SEE NOTE ON FIRST PAGE). * C 2. 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 THE CORRESPONDING EXACT-FIT ROC, AND THEN ABORT * C EXECUTION OF THE BIVARIATE DATA SET. * C 3. SOME DATA SETS MAY CAUSE UNDERFLOW DURING THE MAXIMUM- * C LIKELIHOOD CALCULATIONS. TO CIRCUMVENT THIS PROBLEM, THE * C PROGRAM MAY MERGE TWO OR MORE RATING CATEGORIES IN ONE OR * C BOTH CONDITIONS BEFORE PROCEEDING. THE PROGRAM OUTPUTS BOTH * C THE ORIGINAL AND COLLAPSED DATA MATRICES IN SUCH SITUATIONS. * C IN OUR EXPERIENCE, THIS NEED OCCURS ONLY WITH VERY HIGHLY * C CORRELATED 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 'CORROC2' 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 'CORROC2' * 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(11,11),WW(11,11),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) COMMON/BLK3/TX(10),UY(10),TEXPX(12),TEXPY(12), 1 TXPTUR(12,12),TXPTUP(12,12),FYXP(12,12),FXYP(12,12), 2 FUTR(10,12),FTUR(12,10),ELR(12,12),ELP(12,12), 3 CPELR(11,11),CPELP(11,11),SKP,SKR,RADP,RADR,FUNLIK COMMON/BLK4/FOP(26),SOPNEG(26,26),VCOV(26,26),ESTCOR(26), 1 CORR(26,26),XXDUM(496),ITER,CRIT,LSTAT COMMON/BLK5/FPVAL(26),TPVALX(26),TPVALY(26) COMMON/BLK6/ICLASS,JCLASS,IERX,IERY,ICON,JCON,IARY(11),JARY(11) COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(12),U(12) COMMON/E/IEND,IERROR,IECHO INTEGER TCODE REAL XI(10),XJ(10),T0(12),U0(12) IEND=0 5 IERX=0 IERY=0 ICLASS=0 JCLASS=0 IUNDER=0 CALL READIN IF (IEND.NE.0 .OR. IERROR.EQ.1) GO TO 100 CALL HEADLN 10 CALL CLAPSE CALL TPFFPF crit=1.0 IF(IUNDER.EQ.0)CALL PRTOP 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)GO TO 5 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.EQ.0) GO TO 15 CALL ITRMSG GO TO 5 15 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.EQ.0) GO TO 20 CALL ITRMSG GO TO 5 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 AX0=AX BX0=BX AY0=AY BY0=BY 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) T0(I+1)=T(I+1) 30 CONTINUE DO 40 I=1,KKJ U(I+1)=-XJ(I) U0(I+1)=U(I+1) 40 CONTINUE CALL CORCOE R0=R RHO0=RHO CALL INIVAR 50 ITER=ITER+1 CALL TERMS IF(ITER.EQ.1)FUN0=FUNLIK CALL FIRST CALL SECOND CALL ITERAT(IU) IF(IU.EQ.1)THEN IUNDER=1 GO TO 10 ENDIF IF(IU.EQ.2)GO TO 80 IF (ITER.GT.100)THEN WRITE(6,60) CRIT 60 FORMAT(1X,'PROCEDURE DOES NOT CONVERGE AFTER 100 ITERATIONS'/ 1 'ON THE LAST ITERATION CRIT=',F9.5///) GO TO 5 ENDIF IF (LSTAT.NE.2) GO TO 50 C C GET CORRELATION MATRIX ON FINAL ITERATION C DO 70 I=1,MTRX DO 70 J=1,MTRX CORR(I,J)=VCOV(I,J)/SQRT(VCOV(I,I)*VCOV(J,J)) 70 CONTINUE 80 IF(IUNDER.NE.1)GO TO 90 CALL PRTMTX CALL PRTOP 90 CALL OUTINI(AX0,BX0,AY0,BY0,R0,RHO0,T0,U0,FUN0, + VAX,VBX,VABX,VAY,VBY,VABY) IF(IU.EQ.2)THEN WRITE(6,95) 95 FORMAT(///6X,'**** NO FINAL ESTIMATES CAN BE DERIVED FOR THIS ', + 'DATA SET ****') GO TO 5 ENDIF CALL TPFVAL CALL TEST CALL OUTRSL CALL INFORM GO TO 5 100 STOP 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 READIN C----------------------------------------------- C C READ INPUT DATA C COMMON/BLK1/VV(11,11),WW(11,11),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) COMMON/E/IEND,IERROR,IECHO INTEGER TCODE C C READ IN ECHO PRINT CODE C 10 CALL RDECHO IF(IEND.EQ.1)RETURN C C READ DATA DESCRIPTION, TOTAL NUMBER OF CATEGORY, AND C THE CATEGORY REPRESENTING THE STRONGEST EVIDENCE OF C SIGNAL C CALL READHL(IDENTI,KATI,KSNX) IF(IERROR.EQ.1)GO TO 199 CALL READHL(IDENTJ,KATJ,KSNY) IF(IERROR.EQ.1)GO TO 199 C C READ RATING-DATA PAIRS AND GENERATE RATING-DATA MATRICES C FOR NOISE AND SIGNAL-PLUS-NOISE CASES C CALL READQ(MN,KATI,KATJ,VV) IF(IERROR.EQ.1)GO TO 199 CALL READQ(MS,KATI,KATJ,WW) IF(IERROR.EQ.1)GO TO 199 80 NNI=KATI+1 NNJ=KATJ+1 C C READ "TCODE" C CALL READT(TCODE,FP) 199 RETURN 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=200)(LINE(I),I=1,80) 100 FORMAT(80A1) WRITE(6,110) 110 FORMAT(1H1,'CORROC2:') CALL GETWRD(NOTGET) IF(NOTGET.EQ.0)GO TO 120 IECHO=1 WRITE(6,115) 115 FORMAT(2X,'--- NO INPUT DATA FOR ECHO-PRINTED REQUEST CODE ', + 5X,'ASSIGN BY DEFAULT: ECHO-PRINTED 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 IECHO=1 WRITE(6,125)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 125 FORMAT(2X,80A1) WRITE(6,126) 126 FORMAT(2X,'--- INVALID ECHO-PRINTED REQUEST CODE DETECTED ---'/ + 5X,'ASSIGN BY DEFAULT: YES'/ + 2X,'ECHO-PRINT OF INPUT DATA THAT WILL BE USED IN ', + 'THE ANALYSIS THAT FOLLOWS:') RETURN C 130 IF(LSTRING(1).EQ.KECHO1 .OR. LSTRING(1).EQ.KECHO2)GO TO 140 IECHO=0 WRITE(6,135) 135 FORMAT(2X,'INPUT DATA WILL NOT BE ECHO-PRINTED ') RETURN C 140 WRITE(6,145) 145 FORMAT(2X,'ECHO-PRINT OF INPUT DATA THAT WILL BE USED IN ', + 'THE ANALYSIS THAT FOLLOWS:') IECHO=1 WRITE(6,155)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 155 FORMAT(2X,80A1) RETURN 200 IEND=1 RETURN END C-------------------------------------------- SUBROUTINE READHL(ITITLE,KAT,KSN) C-------------------------------------------- C C THIS SUBROUTINE READS IN A FREE-TEXT DESCRIPTION OF THE DATA, C TOTAL NUMBER OF CATEGORY AND A CODE WHICH INDICATES WHETHER C CATEGORY 1 OR CATEGORY KAT IS ASSOCIATED WITH ACTUALLY C POSITIVE CASES. C COMMON/E/IEND,IERROR,IECHO COMMON/PASS/LSTRING(80),LINE(80),LENGTH DIMENSION ITITLE(*) DATA IBLANK/1H / C KAT=0 KSN=0 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 J=IPNT,80 ITITLE(J-IPNT+1)=LINE(J) 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.EQ.1)GO TO 210 CALL TONUM(LENGTH,LSTRING,RVALUE,IERROR) IF(IERROR.EQ.1)GO TO 210 DO 202 I=3,11 IF(RVALUE.EQ.I)GO TO 203 202 CONTINUE GO TO 210 203 KAT=RVALUE C CALL GETWRD(IERROR) IF(IERROR.EQ.1)GO TO 220 CALL TONUM(LENGTH,LSTRING,RVALUE,IERROR) IF(IERROR.EQ.1)GO TO 220 IF(RVALUE.NE.1.AND.RVALUE.NE.KAT)GO TO 220 KSN=RVALUE IF(IECHO.EQ.1)WRITE(6,205)KAT,KSN 205 FORMAT(2X,I2,3X,I2) RETURN 210 IERROR=1 WRITE(6,215)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 215 FORMAT(2X,80A1/2X,'--- INVALID TOTAL NUMBER OF CATEGORY ', + 'DETECTED ---') RETURN 220 IERROR=1 WRITE(6,225)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 225 FORMAT(2X,80A1/2X,'--- INVALID KSN DETECTED ---') RETURN END C-------------------------------------------- SUBROUTINE READQ(NUM,MAXX,MAXY,QQ) 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 NON NUMERICAL CHARACTER C COMMON/E/IEND,IERROR,IECHO COMMON/PASS/LSTRING(80),LINE(80),LENGTH DIMENSION QQ(11,11) DATA IBLANK/1H /,ISTAR/1H*/ C NUM=0 DO 50 I=1,11 DO 50 J=1,11 QQ(I,J)=0 50 CONTINUE C 100 READ(5,110)(LINE(I),I=1,80) 110 FORMAT(80A1) C C READ RATING-DATA FOR X-CONDITION C 120 CALL GETWRD(IERROR) IF(IERROR.EQ.1)GO TO 130 IF(LSTRING(1).EQ.ISTAR)GO TO 150 CALL TONUM(LENGTH,LSTRING,RVALUE,IERROR) IF(IERROR.EQ.1)GO TO 130 DO 121 I=1,MAXX IF(RVALUE.EQ.I)GO TO 122 121 CONTINUE GO TO 130 122 NUM=NUM+1 IX=RVALUE C C READ RATING-DATA FOR Y-CONDITION C CALL GETWRD(IERROR) IF(IERROR.EQ.1)GO TO 130 CALL TONUM(LENGTH,LSTRING,RVALUE,IERROR) IF(IERROR.EQ.1)GO TO 130 DO 123 I=1,MAXY IF(RVALUE.EQ.I)GO TO 124 123 CONTINUE GO TO 130 124 IY=RVALUE DO 125 I=1,80 IF(LINE(I).EQ.IBLANK)GO TO 125 IPNT=I GO TO 126 125 CONTINUE IPNT=80 126 IF(IECHO.EQ.1)WRITE(6,127)IX,IY,(LINE(I),I=IPNT,80) 127 FORMAT(2X,I2,5X,I2,8X,80A1) QQ(IY,IX)=QQ(IY,IX)+1 GO TO 100 C 130 IERROR=1 WRITE(6,135)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 135 FORMAT(2X,80A1/1X,'--- INVALID RATING-DATA PAIR DETECTED ---') RETURN 150 IF(IECHO.EQ.1)WRITE(6,155)(LSTRING(I),I=1,LENGTH), + (LINE(I),I=1,80) 155 FORMAT(2X,80A1) RETURN END C-------------------------------------------- SUBROUTINE READT(TCODE,FP) C-------------------------------------------- C C READ STATISTICAL TEST CODE 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 READ(5,100)(LINE(I),I=1,80) 100 FORMAT(80A1) TCODE=99 CALL GETWRD(IERROR) IF(IERROR.EQ.1)GO TO 180 IF(LSTRING(1).EQ.IB1 .OR. LSTRING(1).EQ.IB2)TCODE=1 IF(LSTRING(1).EQ.IA1 .OR. LSTRING(1).EQ.IA2)TCODE=2 IF(LSTRING(1).EQ.IT1 .OR. LSTRING(1).EQ.IT2)TCODE=3 IF(TCODE.LT.1 .OR. TCODE.GT.3)GO TO 180 IF(IECHO.EQ.1)WRITE(6,120)(LSTRING(I),I=1,LENGTH) 120 FORMAT(2X,80A1) IF(TCODE.NE.3)RETURN C C READ IN FPF VALUE IF TCODE=3 C READ(5,100)(LINE(I),I=1,80) CALL GETWRD(IERROR) IF(IERROR.EQ.1)GO TO 190 CALL TONUM(LENGTH,LSTRING,RVALUE,IERROR) IF(IERROR.EQ.1)GO TO 190 FP=RVALUE IF(FP.GE.1.OR.FP.LE.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 ---') 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 CONTINUE 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 C 90 NOTGET=1 RETURN END C----------------------------------------------------------- SUBROUTINE TONUM(LEN,LINPUT,RNUM,IERR) C----------------------------------------------------------- C C CONVERT CHARACTER STRING TO REAL NUMBER C DIMENSION LINPUT(*) DATA IZERO/1H0/,IDEC/1H./ C IERR= 0 IF(LEN .EQ. 0) GO TO 120 C C CHECK INVALID CHARACTER (I.E., ANY CHARACTER EXCEPT 0-9, AND PERIOD) C IP= 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 120 IP= IP+1 IF(IP .GT. 1 .OR. (IP.EQ.1.AND.LEN.EQ.1)) GO TO 120 20 CONTINUE C C CONVERT TO NUMERICAL STRING C IU= 0 TOTAL= 0 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) RETURN C 120 IERR= 1 200 RETURN END C---------------------------------------------- INTEGER FUNCTION IORD(LCH) C---------------------------------------------- C C CHANGE FROM CHARACTER TO ASCII CODE C 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 HEADLN 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 REAL RVV(11,11),RWW(11,11),RSUMNT(11),RSUMST(11),RSUMNU(11), 1 RSUMSU(11) COMMON/E/IEND,IERROR,IECHO COMMON/BLK1/VV(11,11),WW(11,11),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) 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 EMN=0. EMS=0. DO 60 I=1,KATI EMN=EMN+SUMNT(I) EMS=EMS+SUMST(I) 60 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,36X,'C O R R O C 2 (JUNE 1993 VERSION) :'// + 25X,'M A X I M U M L I K E L I H O O D E S T I ' 1,'M A T I O N'/39X,'O F T H E P A R A M E T E R S'/25X, 2'O F T H E B I V A R I A T E B I N O R M A L ', 3'M O D E L'/29X,'F O R C O R R E L A T E D R A T I N G ', 4'D A T A') 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),KATI,KSNX,(IDENTJ(J),J=1,80), 1KATJ,KSNY,(I,I=1,KATI) 80 FORMAT(//10X,'CONDITION X: ',80A1,/, 110X,'DATA COLLECTED IN ',I2,' CATEGORIES',/, 210X,'WITH CATEGORY ',I2,' REPRESENTING STRONGEST EVIDENCE ', 3'OF POSITIVITY (E.G., THAT ABNORMALITY IS PRESENT).'// 410X,'CONDITION Y: ',80A1,/, 510X,'DATA COLLECTED IN ',I2,' CATEGORIES',/ 610X,'WITH CATEGORY ',I2,' REPRESENTING STRONGEST EVIDENCE ', 7'OF POSITIVITY (E.G., THAT ABNORMALITY IS PRESENT).', 8///35X,'RATING-DATA MATRIX FOR ACTUALLY NEGATIVE CASES:'// 953X,'CONDITION X RATING'/36X,'CONDITION',3X,11(I2,4X)) WRITE(6,81) 81 FORMAT(36X,'Y RATING') DO 86 J=1,KATJ JOPP=NNJ-J WRITE(6,85)JOPP,(int(VV(JOPP,I)),I=1,KATI),int(SUMNU(JOPP)) 85 FORMAT(40X,I2,4X,12(I4,2X)/) 86 CONTINUE WRITE(6,87) (int(SUMNT(I)),I=1,KATI),int(EMN) 87 FORMAT(/39X,'SUMX',3X,12(I4,2X)/) WRITE(6,88) (I,I=1,KATI) 88 FORMAT(//,35X,'RATING-DATA MATRIX FOR ACTUALLY POSITIVE CASES:', 1//,53X,'CONDITION X RATING'/36X,'CONDITION',3X,11(I2,4X)) WRITE(6,81) DO 89 J=1,KATJ JOPP=NNJ-J WRITE(6,85)JOPP,(int(WW(JOPP,I)),I=1,KATI),int(SUMSU(JOPP)) 89 CONTINUE WRITE(6,87) (int(SUMST(I)),I=1,KATI),int(EMS) IF(KSNX.EQ.KATI)GO TO 200 DO 110 I=1,KATI DO 110 J=1,KATJ RVV(J,I)=VV(J,KATI+1-I) 110 RWW(J,I)=WW(J,KATI+1-I) DO 120 I=1,KATI DO 120 J=1,KATJ VV(J,I)=RVV(J,I) 120 WW(J,I)=RWW(J,I) DO 130 I=1,KATI RSUMNT(KATI-I+1)=SUMNT(I) 130 RSUMST(KATI-I+1)=SUMST(I) DO 140 I=1,KATI SUMNT(I)=RSUMNT(I) 140 SUMST(I)=RSUMST(I) C 200 IF(KSNY.EQ.KATJ)GO TO 300 DO 210 I=1,KATI DO 210 J=1,KATJ RVV(J,I)=VV(KATJ+1-J,I) 210 RWW(J,I)=WW(KATJ+1-J,I) DO 220 I=1,KATI DO 220 J=1,KATJ VV(J,I)=RVV(J,I) 220 WW(J,I)=RWW(J,I) DO 230 J=1,KATJ RSUMNU(KATJ-J+1)=SUMNU(J) 230 RSUMSU(KATJ-J+1)=SUMSU(J) DO 240 J=1,KATJ SUMNU(J)=RSUMNU(J) 240 SUMSU(J)=RSUMSU(J) 300 RETURN END C------------------------------------------------------------------ SUBROUTINE CLAPSE C----------------------------------------------------------------- C C COLLAPSE DATA C COMMON/BLK1/VV(11,11),WW(11,11),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) 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 EMN=0. EMS=0. DO 60 I=1,KATI EMN=EMN+SUMNT(I) EMS=EMS+SUMST(I) 60 CONTINUE C KATI1=KATI DO 1005 I=1,KATI IT=KATI-I+1 IF (SUMNT(IT).NE.0.OR.SUMST(IT).NE.0) GO TO 1005 KATI1=KATI1-1 IF (IT.GT.KATI1) GO TO 1005 DO 1004 J=IT,KATI1 SUMNT(J)=SUMNT(J+1) SUMST(J)=SUMST(J+1) DO 1003 IX=1,KATJ VV(IX,J)=VV(IX,J+1) WW(IX,J)=WW(IX,J+1) 1003 CONTINUE 1004 CONTINUE 1005 CONTINUE KATI=KATI1 KATJ1=KATJ DO 1015 I=1,KATJ IT=KATJ-I+1 IF (SUMNU(IT).NE.0.OR.SUMSU(IT).NE.0) GO TO 1015 KATJ1=KATJ1-1 IF (IT.GT.KATJ1) GO TO 1015 DO 1010 J=IT,KATJ1 SUMNU(J)=SUMNU(J+1) SUMSU(J)=SUMSU(J+1) DO 1007 IY=1,KATI VV(J,IY)=VV(J+1,IY) WW(J,IY)=WW(J+1,IY) 1007 CONTINUE 1010 CONTINUE 1015 CONTINUE KATJ=KATJ1 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(11,11),WW(11,11),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) INTEGER TCODE FPFX(NNI)=0. TPFX(NNI)=0. FPFX(KATI)=SUMNT(KATI) TPFX(KATI)=SUMST(KATI) DO 1020 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) 1020 CONTINUE FPFY(NNJ)=0. TPFY(NNJ)=0. FPFY(KATJ)=SUMNU(KATJ) TPFY(KATJ)=SUMSU(KATJ) DO 1025 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) 1025 CONTINUE DO 1030 I=1,KATI FPFX(I)=FPFX(I)/EMN TPFX(I)=TPFX(I)/EMS 1030 CONTINUE DO 1035 I=1,KATJ FPFY(I)=FPFY(I)/EMN TPFY(I)=TPFY(I)/EMS 1035 CONTINUE 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(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) COMMON/BLK6/ICLASS,JCLASS,IERX,IERY,ICON,JCON,IARY(11),JARY(11) C IF (ICLASS.EQ.0) GO TO 1120 IF (ICLASS.GE.5.AND.ICLASS.LE.8) WRITE(6,1105) IF (ICLASS.LT.5.OR.ICLASS.GT.8) WRITE(6,1110) 1105 FORMAT(///,15X,'MARGINAL DATA FOR CONDITION X ARE DEGENERATE.') 1110 FORMAT(///,15X,'MARGINAL DATA FOR CONDITION X ARE DEGENERATE') CALL DEGENE(ICLASS,IARY,FPFX,TPFX) 1120 IF (JCLASS.EQ.0) GO TO 1135 IF (JCLASS.GE.5.AND.JCLASS.LE.8) WRITE(6,1125) IF (JCLASS.LT.5.OR.JCLASS.GT.8) WRITE(6,1130) 1125 FORMAT(///,15X,'MARGINAL DATA FOR CONDITION Y ARE DEGENERATE.') 1130 FORMAT(///,15X,'MARGINAL DATA FOR CONDITION Y ARE DEGENERATE') CALL DEGENE(JCLASS,JARY,FPFY,TPFY) 1135 WRITE(6,1140) 1140 FORMAT(//,15X,'EXECUTION OF THIS BIVARIATE DATA SET ABORTED.') KEYDEG=1 RETURN END C-------------------------------------------------------------------- SUBROUTINE ITRMSG C-------------------------------------------------------------------- C C WRITE OUT MESSAGE OF MATRIX INVERSION PROBLEM OR ITERATIONS C EXCEEDING MAXIMUM RESTRICTION FOR MARGINAL DATA SETS C COMMON/BLK6/ICLASS,JCLASS,IERX,IERY,ICON,JCON,IARY(11),JARY(11) IF (IERX.EQ.0) GO TO 1205 IF (IERX.GT.100) WRITE(6,1203) IF (IERX.LE.100) WRITE(6,1204) IERX 1203 FORMAT(1X,'PROCEDURE DOES NOT CONVERGE -- 100 ITERATIONS --', 1 'THIS WAS WITH MARGINAL DATA SET FOR CONDITION X') 1204 FORMAT(1X,'PROBLEM WITH MATRIX INVERSION ON ',I3,'TH ITERATION', 2 ' WITH MARGINAL DATA SET FOR CONDITION X') IF (IERY.EQ.0) RETURN 1205 IF (IERY.GT.100) WRITE(6,1206) IF (IERY.LE.100) WRITE(6,1207) IERY 1206 FORMAT(1X,'PROCEDURE DOES NOT CONVERGE -- 100 ITERATIONS --', 1 'THIS WAS WITH MARGINAL DATA SET FOR CONDITION Y') 1207 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(11,11),WW(11,11),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(12),U(12) INTEGER TCODE SUMRX=0. SUMPX=0. SUMRY=0. SUMPY=0. SUMRX2=0. SUMPX2=0. SUMRY2=0. SUMPY2=0. SUMRXY=0. SUMPXY=0. DO 1065 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 1065 CONTINUE DO 1070 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 1070 CONTINUE DO 1075 I=1,KATI DO 1075 J=1,KATJ SUMRXY=SUMRXY+VV(J,I)*J*I SUMPXY=SUMPXY+WW(J,I)*J*I 1075 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 OUTINI(AX,BX,AY,BY,R,RHO,T,U,FUNLIK, + VAX,VBX,VABX,VAY,VBY,VABY) C------------------------------------------------------------------- C C PRINT INITIAL ESTIMATES OF ALL OF THE MODEL PARAMETERS TO BE C USED NEXT IN CALCULATING PARAMETERS OF THE BIVARIATE MODEL BY C THE 'METHOD OF SCORING'. C DIMENSION T(*),U(*) COMMON/BLK1/VV(11,11),WW(11,11),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) INTEGER TCODE C WRITE(6,1095) AX,BX,AY,BY,R,RHO 1095 FORMAT(//42X,'INITIAL ESTIMATES OF THE PARAMETERS:'/ 132X,'AX=',F7.4,5X,'BX=',F7.4,5X,'AY=',F7.4,5X,'BY=',F7.4/ 232X,'R(NEGATIVE CASES)=',F7.4,5X,'R(POSITIVE CASES)=',F7.4) WRITE(6,1096) (T(I),I=2,KATI) 1096 FORMAT(34X,'T(I)',10(2X,F6.3)) WRITE(6,1097) (U(I),I=2,KATJ) 1097 FORMAT(34X,'U(J)',10(2X,F6.3)) CORX=VABX/SQRT(VAX*VBX) CORY=VABY/SQRT(VAY*VBY) WRITE(6,1098) VAX,VBX,VABX,CORX,VAY,VBY,VABY,CORY 1098 FORMAT(7X,'(VAX=',F6.4,2X,'VBX=',F6.4,2X,'COVABX=',F7.4,2X, 1'CORABX=',F7.4,4X,'VAY=',F6.4,2X,'VBY=',F6.4,2X,'COVABY=',F7.4,2X, 2'CORABY=',F7.4,')') WRITE(6,1099)FUNLIK 1099 FORMAT(/36X,'LOGL =',F10.3) RETURN END C--------------------------------------------------------------------- SUBROUTINE INIVAR C--------------------------------------------------------------------- C C INITIALIZE VARIABLES C COMMON/BLK1/VV(11,11),WW(11,11),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) COMMON/BLK3/TX(10),UY(10),TEXPX(12),TEXPY(12), 1 TXPTUR(12,12),TXPTUP(12,12),FYXP(12,12),FXYP(12,12), 2 FUTR(10,12),FTUR(12,10),ELR(12,12),ELP(12,12), 3 CPELR(11,11),CPELP(11,11),SKP,SKR,RADP,RADR,FUNLIK COMMON/BLK4/FOP(26),SOPNEG(26,26),VCOV(26,26),ESTCOR(26), 1 CORR(26,26),XXDUM(496),ITER,CRIT,LSTAT INTEGER TCODE C TEXPX(1)=0. TEXPX(NNI)=0. TEXPY(1)=0. TEXPY(NNI)=0. DO 1100 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. 1100 CONTINUE DO 1105 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. 1105 CONTINUE DO 1110 I=1,NNI FXYP(I,1)=0. FXYP(I,NNJ)=0. 1110 CONTINUE DO 1115 J=2,KATJ FXYP(1,J)=0. FXYP(NNI,J)=1. FTUR(NNI,J-1)=1. FTUR(1,J-1)=0. 1115 CONTINUE DO 1120 I=2,NNI ELP(I,1)=0. ELR(I,1)=0. 1120 CONTINUE ELP(NNI,NNJ)=1. ELR(NNI,NNJ)=1. MTRX=NNI+NNJ+2 DO 1125 I=1,MTRX DO 1125 J=1,MTRX SOPNEG(I,J)=0. VCOV(I,J)=0. CORR(I,J)=0. 1125 CONTINUE RETURN END C--------------------------------------------------------------------- SUBROUTINE TERMS C--------------------------------------------------------------------- C C COMPUTE THE MAJOR TERMS THAT OCCUR IN THE DERIVATIVE EXPRESSIONS. C COMMON/BLK1/VV(11,11),WW(11,11),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) COMMON/BLK3/TX(10),UY(10),TEXPX(12),TEXPY(12), 1 TXPTUR(12,12),TXPTUP(12,12),FYXP(12,12),FXYP(12,12), 2 FUTR(10,12),FTUR(12,10),ELR(12,12),ELP(12,12), 3 CPELR(11,11),CPELP(11,11),SKP,SKR,RADP,RADR,FUNLIK COMMON/BLK4/FOP(26),SOPNEG(26,26),VCOV(26,26),ESTCOR(26), 1 CORR(26,26),XXDUM(496),ITER,CRIT,LSTAT COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(12),U(12) 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 1135 I=2,KATI TX(I-1)=T(I)*BX-AX TEXPX(I)=EXP(-.5*TX(I-1)**2) 1135 CONTINUE DO 1140 J=2,KATJ UY(J-1)=U(J)*BY-AY TEXPY(J)=EXP(-.5*UY(J-1)**2) 1140 CONTINUE DO 1145 I=1,KKI DO 1145 J=2,KATJ ARG=(U(J)-R*T(I+1))/RADR CALL NDTR(ARG,FUTR(I,J),D) 1145 CONTINUE DO 1150 J=1,KKJ DO 1150 I=2,KATI ARG=(T(I)-R*U(J+1))/RADR CALL NDTR(ARG,FTUR(I,J),D) 1150 CONTINUE DO 1155 I=2,KATI DO 1155 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) 1155 CONTINUE IF (ITER.GT.1) RETURN DO 1156 I=2,KATI CALL NDTR(TX(I-1),ELP(I,NNJ),D) 1156 CALL NDTR (T(I),ELR(I,NNJ),D) DO 1157 J=2,KATJ CALL NDTR(UY(J-1),ELP(NNI,J),D) 1157 CALL NDTR(U(J),ELR(NNI,J),D) DO 1158 I=2,KATI DO 1158 J=2,KATJ CALL MDBNOR(TX(I-1),UY(J-1),RHO,ELP(I,J),IER) 1158 CALL MDBNOR(T(I),U(J),R,ELR(I,J),IER) DO 1160 I=1,KATI DO 1160 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 1160 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(11,11),WW(11,11),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) COMMON/BLK3/TX(10),UY(10),TEXPX(12),TEXPY(12), 1 TXPTUR(12,12),TXPTUP(12,12),FYXP(12,12),FXYP(12,12), 2 FUTR(10,12),FTUR(12,10),ELR(12,12),ELP(12,12), 3 CPELR(11,11),CPELP(11,11),SKP,SKR,RADP,RADR,FUNLIK COMMON/BLK4/FOP(26),SOPNEG(26,26),VCOV(26,26),ESTCOR(26), 1 CORR(26,26),XXDUM(496),ITER,CRIT,LSTAT COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(12),U(12) INTEGER TCODE GPI=3.14159265 C C WITH RESPECT TO AX C 1162 SUM=0. DO 1163 I=1,KATI DO 1163 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) 1163 CONTINUE FOP(1)=SUM/SQRT(2.*GPI) C C WITH RESPECT TO BX C SUM=0. DO 1165 I=1,KATI DO 1165 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) 1165 CONTINUE FOP(2)=SUM/SQRT(2.*GPI) C C WITH RESPECT TO AY C SUM=0. DO 1170 I=1,KATI DO 1170 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) 1170 CONTINUE FOP(3)=SUM/SQRT(2.*GPI) C C WITH RESPECT TO BY C SUM=0. DO 1175 I=1,KATI DO 1175 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) 1175 CONTINUE FOP(4)=SUM/SQRT(2.*GPI) C C WITH RESPECT TO R=R(NOISE) C SUM=0. DO 1180 I=1,KATI DO 1180 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) 1180 CONTINUE FOP(5)=SUM/(2.*GPI*RADR) C C WITH RESPECT TO RHO=R(SIGNAL-PLUS-NOISE) C SUM=0. DO 1185 I=1,KATI DO 1185 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) 1185 CONTINUE FOP(6)=SUM/(2.*GPI*RADP) C C WITH RESPECT TO EACH OF THE T(I) C DO 1195 I=1,KKI SUMN=0. SUMS=0. DO 1190 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)) 1190 CONTINUE FOP(I+6)=(SUMN+BX*SUMS)/SQRT(2.*GPI) 1195 CONTINUE C C WITH RESPECT TO EACH OF THE U(J) C DO 1205 J=1,KKJ SUMN=0. SUMS=0. DO 1200 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)) 1200 CONTINUE FOP(J+KKI+6)=(SUMN+BY*SUMS)/SQRT(2.*GPI) 1205 CONTINUE RETURN END C------------------------------------------------------------------ SUBROUTINE SECOND C------------------------------------------------------------------ C C COMPUTE THE EXPECTED VALUE OF EACH OF THE MINUS SECOND-ORDER C PARTIAL DERIVATIVE EQUATIONS. C COMMON/BLK1/VV(11,11),WW(11,11),KATI,KATJ,TCODE,FP,NNI,NNJ,KKI, 1 KSNX,KSNY,KKJ,KKIL,KKJL,MTRX,IDENTI(80),IDENTJ(80) COMMON/BLK2/EMN,EMS,SUMNT(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) COMMON/BLK3/TX(10),UY(10),TEXPX(12),TEXPY(12), 1 TXPTUR(12,12),TXPTUP(12,12),FYXP(12,12),FXYP(12,12), 2 FUTR(10,12),FTUR(12,10),ELR(12,12),ELP(12,12), 3 CPELR(11,11),CPELP(11,11),SKP,SKR,RADP,RADR,FUNLIK COMMON/BLK4/FOP(26),SOPNEG(26,26),VCOV(26,26),ESTCOR(26), 1 CORR(26,26),XXDUM(496),ITER,CRIT,LSTAT COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(12),U(12) INTEGER TCODE GPI=3.14159265 C C WITH RESPECT TO AX C SUM=0. DO 1210 I=1,KATI DO 1210 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) 1210 CONTINUE SOPNEG(1,1)=EMS*SUM/(2.*GPI) C C WITH RESPECT TO BX C SUM=0. DO 1215 I=1,KATI DO 1215 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) 1215 CONTINUE SOPNEG(2,2)=EMS*SUM/(2.*GPI) C C WITH RESPECT TO AY C SUM=0. DO 1220 I=1,KATI DO 1220 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) 1220 CONTINUE SOPNEG(3,3)=EMS*SUM/(2.*GPI) C C WITH RESPECT TO BY C SUM=0. DO 1225 I=1,KATI DO 1225 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) 1225 CONTINUE SOPNEG(4,4)=EMS*SUM/(2.*GPI) C C WITH RESPECT TO R=R(NOISE) C SUM=0. DO 1230 I=1,KATI DO 1230 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) 1230 CONTINUE SOPNEG(5,5)=EMN*SUM/(4.*GPI**2*SKR) C C WITH RESPECT TO RHO=R(SIGNAL-PLUS-NOISE) C SUM=0. DO 1235 I=1,KATI DO 1235 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) 1235 CONTINUE SOPNEG(6,6)=EMS*SUM/(4.*GPI**2*SKP) C C WITH RESPECT TO EACH OF THE T(I) C DO 1245 I=1,KKI SUMN=0. SUMS=0. DO 1240 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)) 1240 CONTINUE SOPNEG(I+6,I+6)=(EMN*SUMN+EMS*BX**2*SUMS)/(2.*GPI) 1245 CONTINUE C C WITH RESPECT TO EACH OF THE U(J) C DO 1255 J=1,KKJ SUMN=0. SUMS=0. DO 1250 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)) 1250 CONTINUE SOPNEG(J+KKI+6,J+KKI+6)=(EMN*SUMN+EMS*BY**2*SUMS)/(2.*GPI) 1255 CONTINUE C C MIXED WITH RESPECT TO AX AND BX C SUM=0. DO 1260 I=1,KATI DO 1260 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) 1260 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 1265 I=1,KATI DO 1265 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) 1265 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 1270 I=1,KATI DO 1270 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) 1270 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 1275 I=1,KATI DO 1275 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) 1275 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 1285 I=1,KKI SUM=0. DO 1280 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)) 1280 CONTINUE SOPNEG(1,I+6)=EMS*BX*SUM/(2.*GPI) SOPNEG(I+6,1)=SOPNEG(1,I+6) 1285 CONTINUE C C MIXED WITH RESPECT TO AX AND EACH OF THE U(J) C DO 1295 J=1,KKJ SUM=0. DO 1290 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)) 1290 CONTINUE SOPNEG(1,J+KKI+6)=EMS*BY*SUM/(2.*GPI) SOPNEG(J+KKI+6,1)=SOPNEG(1,J+KKI+6) 1295 CONTINUE C C MIXED WITH RESPECT TO BX AND AY C SUM=0. DO 1300 I=1,KATI DO 1300 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) 1300 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 1305 I=1,KATI DO 1305 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) 1305 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 1310 I=1,KATI DO 1310 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) 1310 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 1320 I=1,KKI SUM=0. DO 1315 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)) 1315 CONTINUE SOPNEG(2,I+6)=EMS*BX*SUM/(2.*GPI) SOPNEG(I+6,2)=SOPNEG(2,I+6) 1320 CONTINUE C C MIXED WITH RESPECT TO BX AND EACH OF THE U(J) C DO 1330 J=1,KKJ SUM=0. DO 1325 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)) 1325 CONTINUE SOPNEG(2,J+KKI+6)=EMS*BY*SUM/(2.*GPI) SOPNEG(J+KKI+6,2)=SOPNEG(2,J+KKI+6) 1330 CONTINUE C C MIXED WITH RESPECT TO AY AND BY C SUM=0. DO 1335 I=1,KATI DO 1335 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) 1335 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 1340 I=1,KATI DO 1340 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) 1340 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 1350 I=1,KKI SUM=0. DO 1345 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)) 1345 CONTINUE SOPNEG(3,I+6)=EMS*BX*SUM/(2.*GPI) SOPNEG(I+6,3)=SOPNEG(3,I+6) 1350 CONTINUE C C MIXED WITH RESPECT TO AY AND EACH OF THE U(J) C DO 1360 J=1,KKJ SUM=0. DO 1355 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)) 1355 CONTINUE SOPNEG(3,J+KKI+6)=EMS*BY*SUM/(2.*GPI) SOPNEG(J+KKI+6,3)=SOPNEG(3,J+KKI+6) 1360 CONTINUE C C MIXED WITH RESPECT TO BY AND RHO=R(SIGNAL-PLUS-NOISE) C SUM=0. DO 1365 I=1,KATI DO 1365 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) 1365 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 1375 I=1,KKI SUM=0. DO 1370 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)) 1370 CONTINUE SOPNEG(4,I+6)=EMS*BX*SUM/(2.*GPI) SOPNEG(I+6,4)=SOPNEG(4,I+6) 1375 CONTINUE C C MIXED WITH RESPECT TO BY AND EACH OF THE U(J) C DO 1385 J=1,KKJ SUM=0. DO 1380 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)) 1380 CONTINUE SOPNEG(4,J+KKI+6)=EMS*BY*SUM/(2.*GPI) SOPNEG(J+KKI+6,4)=SOPNEG(4,J+KKI+6) 1385 CONTINUE C C MIXED WITH RESPECT TO R=R(NOISE) AND EACH OF THE T(I) C DO 1395 I=1,KKI SUM=0. DO 1390 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)) 1390 CONTINUE SOPNEG(5,I+6)=EMN*SUM/SQRT(8.*GPI**3*SKR) SOPNEG(I+6,5)=SOPNEG(5,I+6) 1395 CONTINUE C C MIXED WITH RESPECT TO R=R(NOISE) AND EACH OF THE U(J) C DO 1405 J=1,KKJ SUM=0. DO 1400 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)) 1400 CONTINUE SOPNEG(5,J+KKI+6)=EMN*SUM/SQRT(8.*GPI**3*SKR) SOPNEG(J+KKI+6,5)=SOPNEG(5,J+KKI+6) 1405 CONTINUE C C MIXED WITH RESPECT TO RHO=R(SIGNAL-PLUS-NOISE) AND EACH OF THE T(I) C DO 1415 I=1,KKI SUM=0. DO 1410 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)) 1410 CONTINUE SOPNEG(6,I+6)=EMS*BX*SUM/SQRT(8.*GPI**3*SKP) SOPNEG(I+6,6)=SOPNEG(6,I+6) 1415 CONTINUE C C MIXED WITH RESPECT TO RHO=R(SIGNAL-PLUS-NOISE) AND EACH OF THE U(J) C DO 1425 J=1,KKJ SUM=0. DO 1420 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)) 1420 CONTINUE SOPNEG(6,J+KKI+6)=EMS*BY*SUM/SQRT(8.*GPI**3*SKP) SOPNEG(J+KKI+6,6)=SOPNEG(6,J+KKI+6) 1425 CONTINUE C C MIXED WITH RESPECT TO T(I) AND T(I+1) FOR I=1,KKI-1 C DO 1435 I=1,KKIL SUMN=0. SUMS=0. DO 1430 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) 1430 CONTINUE SOPNEG(I+6,I+7)=(EMN*SUMN+EMS*BX**2*SUMS)/(2.*GPI) SOPNEG(I+7,I+6)=SOPNEG(I+6,I+7) 1435 CONTINUE C C MIXED WITH RESPECT TO EACH OF THE T(I) AND EACH OF THE U(J) C DO 1440 I=1,KKI DO 1440 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) 1440 CONTINUE C C MIXED WITH RESPECT TO U(J) AND U(J+1)FOR J=1,KKJ-1 C DO 1450 J=1,KKJL SUMN=0. SUMS=0. DO 1445 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) 1445 CONTINUE