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 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) 1450 CONTINUE RETURN END 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(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 c ALPHA=1.0 BNDRY=1.0E-07 FUNLI1=FUNLIK CALL MSTR(SOPNEG,XXDUM,26,0,1) CALL SINV(XXDUM,MTRX,.001,IER) IF(IER.LT.0) THEN 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 MSTR(XXDUM,VCOV,26,1,0) 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 C CHECK STOPPINT 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 THE 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.BNDRY .AND. VV(JOP,IOP).NE.0) .OR. + (CPELP(IOP,JOP).LT.BNDRY .AND. WW(JOP,IOP).NE.0))THEN ALPHA=-0.5*ABS(ALPHA) IF(ABS(ALPHA).LT.BNDRY)THEN CALL TRUNCT(IOP,JOP,IU) RETURN ENDIF ENDIF IF(CPELR(IOP,JOP).LT.BNDRY) CPELR(IOP,JOP)=BNDRY IF(CPELP(IOP,JOP).LT.BNDRY) CPELP(IOP,JOP)=BNDRY 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) cc write(*,11)ax,bx,ay,by,r,rho,funlik,funli1 cc11 format(1x,4(f7.4,1x),2(f6.3,1x),2(f9.3,1x)) IF (FUNLIK.GE.FUNLI1) RETURN FUNDIF=ABS((FUNLIK-FUNLI1)/crit) IF (alpha.le.1.0e-5.or.FUNDIF.LE.1.0E-4) RETURN ALPHA=-0.5*ABS(ALPHA) IF(ABS(ALPHA).GT.BNDRY)GO TO 1450 C CALL SELCEL(I,J) C CALL TRUNCT(I,J,IU) 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 IDENTI & J set to 80 to match the rest of the program 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/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 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(11,11),WW(11,11),CPELR(11,11),CPELP(11,11) 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 IDENTI & j set changed from 60 to 80 2/12/93 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) DIMENSION FPF(12),TPF(12),D(3,4),IARY(11) 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 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 RATGET 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 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(11),SUMST(11),SUMNU(11),SUMSU(11), 1 FPFX(12),TPFX(12),FPFY(12),TPFY(12) 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 OUTRSL C------------------------------------------------------------------- C C COMPLETE THE PRINTED OUTPUT PAGE 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/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 WRITE(6,1505)ITER,AX,BX,AY,BY,R,RHO 1505 FORMAT(//9X,'PROCEDURE CONVERGES AFTER ',I3,' ITERATIONS'// 1 45X,'FINAL ESTIMATES OF THE PARAMETERS:'/ 2 32X,'AX=',F7.4,5X,'BX=',F7.4,5X,'AY=',F7.4,5X,'BY=',F7.4/ 3 32X,'R(NEGATIVE CASES)=',F7.4,5X,'R(POSITIVE CASES)=',F7.4) WRITE(6,1506) (T(I),I=2,KATI) WRITE(6,1507) (U(I),I=2,KATJ) 1506 FORMAT(34X,'T(I)',10(2X,F6.3)) 1507 FORMAT(34X,'U(J)',10(2X,F6.3)) FUNLIK=VLOGL(KATI,KATJ,VV,WW,CPELR,CPELP) WRITE(6,1515)FUNLIK 1515 FORMAT(/36X,'LOGL =',F10.3) WRITE(6,1520) 1520 FORMAT(1H1,////,34X,'V A R I A N C E - C O V A R I A N C E M A', 1' T R I X:',//) WRITE(6,1521) (VCOV(1,J),J=1,MTRX) 1521 FORMAT(2X,'AX ',2X,14(1X,F7.4),/,10X,12(1X,F7.4)) WRITE(6,1522) (VCOV(2,J),J=1,MTRX) 1522 FORMAT(2X,'BX ',2X,14(1X,F7.4),/,10X,12(1X,F7.4)) WRITE(6,1523) (VCOV(3,J),J=1,MTRX) 1523 FORMAT(2X,'AY ',2X,14(1X,F7.4),/,10X,12(1X,F7.4)) WRITE(6,1524) (VCOV(4,J),J=1,MTRX) 1524 FORMAT(2X,'BY ',2X,14(1X,F7.4),/,10X,12(1X,F7.4)) WRITE(6,1525) (VCOV(5,J),J=1,MTRX) 1525 FORMAT(2X,'R(-)',2X,14(1X,F7.4),/,10X,12(1X,F7.4)) WRITE(6,1526) (VCOV(6,J),J=1,MTRX) 1526 FORMAT(2X,'R(+)',2X,14(1X,F7.4),/,10X,12(1X,F7.4)) IPRI=6+KKI IPRJ1=IPRI+1 DO 1528 I=7,IPRI II=I-6 WRITE(6,1527) II,(VCOV(I,J),J=1,MTRX) 1527 FORMAT(2X,'T(',I2,')',1X,14(1X,F7.4),/,12X,12(1X,F7.4)) 1528 CONTINUE DO 1530 I=IPRJ1,MTRX II=I-IPRI WRITE(6,1529) II,(VCOV(I,J),J=1,MTRX) 1529 FORMAT(2X,'U(',I2,')',1X,14(1X,F7.4),/,12X,12(1X,F7.4)) 1530 CONTINUE WRITE(6,1531) 1531 FORMAT(////43X,'C O R R E L A T I O N M A T R I X:'//) WRITE(6,1521) (CORR(1,J),J=1,MTRX) WRITE(6,1522) (CORR(2,J),J=1,MTRX) WRITE(6,1523) (CORR(3,J),J=1,MTRX) WRITE(6,1524) (CORR(4,J),J=1,MTRX) WRITE(6,1525) (CORR(5,J),J=1,MTRX) WRITE(6,1526) (CORR(6,J),J=1,MTRX) DO 1532 I=7,IPRI II=I-6 WRITE(6,1527) II,(CORR(I,J),J=1,MTRX) 1532 CONTINUE DO 1533 I=IPRJ1,MTRX II=I-IPRI WRITE(6,1529) II,(CORR(I,J),J=1,MTRX) 1533 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(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/BLK7/AX,BX,AY,BY,R,RHO,T(12),U(12) DO 1605 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) 1605 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE TEST 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/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/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) COMMON/BLK8/AZX,AZY,VAZX,VAZY,COVAZ,RAREA,XSTAT,PLEVEL, 1 ZSTAT,ZLEVEL,TPSTAT,TPLEVL INTEGER TCODE 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 1549 IF (TCODE.EQ.3) GO TO 1550 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 1549 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 1550 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)-2.*Z*VCOV(1,2) VTPY=VCOV(3,3)+Z*Z*VCOV(4,4)-2.*Z*VCOV(3,4) COVTP=VCOV(1,3)-Z*VCOV(1,4)-Z*VCOV(2,3)+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 INFORM C-------------------------------------------------------------------- C C OUTPUT THE TPF VALUES AT SELECTED FPF VALUES, AND ALSO C AREAS, VARIANCES, COVARIANCE, AND CORRELATION. C 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/BLK5/FPVAL(26),TPVALX(26),TPVALY(26) COMMON/BLK7/AX,BX,AY,BY,R,RHO,T(12),U(12) COMMON/BLK8/AZX,AZY,VAZX,VAZY,COVAZ,RAREA,XSTAT,PLEVEL, 1 ZSTAT,ZLEVEL,TPSTAT,TPLEVL INTEGER TCODE WRITE(6,1610) 1610 FORMAT(1H1,14X,1H','PLOTS',1H',' OF THE FITTED BINORMAL', + ' ROC CURVES:', 1 //,5X,'**********************************', 2 '*********************************', 3 ///,22X,'FPF',7X,'TPF FOR',8X,'TPF FOR',/,31X,'CONDITION X', 4 4X,'CONDITION Y',/,21X,'-----',4X,'-------------',2X, 5 '-------------',/) DO 1630 I=1,26 WRITE(6,1625) FPVAL(I),TPVALX(I),TPVALY(I) 1625 FORMAT(21X,F5.3,6X,F6.4,9X,F6.4) 1630 CONTINUE STDAZX=SQRT(VAZX) STDAZY=SQRT(VAZY) WRITE(6,1635) AZX,AZY,STDAZX,STDAZY,RAREA 1635 FORMAT(////,12X,'AREA(X)= ',F7.4,5X,'AREA(Y)= ',F7.4/12X, 1 'STD DEV OF AREA(X)= ',F7.4,5X,'STD DEV OF AREA(Y)= ',F7.4 2 /12X,'CORRELATION OF AREA(X) AND AREA(Y) =',F7.4) IF(IERROR.EQ.0)GO TO 1639 IF(IERROR.EQ.2)WRITE(6,1636) 1636 FORMAT(///5X,'*** INVALID INPUT OF STATISTICAL TEST CODE ', 1 '(TCODE), NO STATISTICAL TEST CALCULATED ***') IF(IERROR.EQ.3)WRITE(6,1637) 1637 FORMAT(///5X,'*** INVALID INPUT OF FPF VALUE AT WHICH THE TPF', 1 1H','S ARE TO BE COMPARED; STATISTICAL TEST ABORTED ***') RETURN C C OUTPUT THE RESULTS OF THE REQUESTED STATISTICAL TEST C 1639 WRITE(6,1640) 1640 FORMAT(/////,23X,'STATISTICAL SIGNIFICANCE OF THE DIFFERENCE',/ 1 16X,'BETWEEN THE TWO CORRELATED ROC CURVES ACCORDING TO', 2 ' THE SELECTED TEST:',//,2X,'*******', 3 '*********************************************************', 3 '**********************') IF (TCODE.EQ.2) GO TO 1642 IF (TCODE.EQ.3) GO TO 1648 WRITE(6,1641) XSTAT,PLEVEL 1641 FORMAT(//,2X,'THE COMPUTED CORRELATED ',1H','BIVARIATE CHI-', 1 'SQUARE TEST',1H',' STATISTIC VALUE IS',F8.4,/,6X, 2 'WITH A CORRESPONDING P-LEVEL OF',F7.4,'.') RETURN 1642 ONEPL=0.5*ZLEVEL WRITE(6,1645) ZSTAT,ZLEVEL,ONEPL 1645 FORMAT(//,2X,'THE COMPUTED CORRELATED ',1H','AREA TEST',1H', 1' STATISTIC VALUE IS',2X,F8.4,/,6X,'WITH A CORRESPONDING TWO-', 2'TAILED P-LEVEL OF',F7.4,' AND ONE-TAILED P-LEVEL OF',F7.4,'.') RETURN 1648 ONEPL=0.5*TPLEVL WRITE(6,1650) FP,TPSTAT,TPLEVL,ONEPL 1650 FORMAT(//,2X,'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,/,6X,'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 PRTMTX C------------------------------------------------------------------ C C PRINT OUT RATING DATA MATRIX AND OBSERVED OPERATING POINTS C REAL RVV(11,11),RWW(11,11),RSUMNT(11),RSUMST(11),RSUMNU(11), 1 RSUMSU(11) 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 RSUMNT(I)=0. RSUMST(I)=0. 40 CONTINUE DO 45 J=1,KATJ RSUMNU(J)=0. RSUMSU(J)=0. 45 CONTINUE DO 50 I=1,KATI DO 50 J=1,KATJ IF(KSNX.NE.1)GO TO 46 RVV(J,I)=VV(J,KATI+1-I) RWW(J,I)=WW(J,KATI+1-I) GO TO 47 46 RVV(J,I)=VV(J,I) RWW(J,I)=WW(J,I) 47 RSUMNT(I)=RSUMNT(I)+RVV(J,I) RSUMST(I)=RSUMST(I)+RWW(J,I) 50 CONTINUE DO 55 J=1,KATJ DO 55 I=1,KATI IF(KSNY.NE.1)GO TO 51 RVV(J,I)=VV(KATJ+1-J,I) RWW(J,I)=WW(KATJ+1-J,I) GO TO 52 51 RVV(J,I)=VV(J,I) RWW(J,I)=WW(J,I) 52 RSUMNU(J)=RSUMNU(J)+RVV(J,I) RSUMSU(J)=RSUMSU(J)+RWW(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(//6X,'DUE TO UNDERFLOW PROBLEMS, IT IS NECESSARY TO ', + 'MERGE TWO OR MORE CATEGORIES IN ONE OR BOTH CONDITIONS.'/ + 6X,'THE PARAMETER ESTIMATES WILL BE BASED ON THE ', + 'FOLLOWING PARTIALLY COLLAPSED DATA MATRICES:'//) WRITE(6,80)(I,I=1,KATI) 80 FORMAT(//,35X,'RATING-DATA MATRIX FOR ACTUALLY NEGATIVE CASES:' + //53X,'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(RVV(JOPP,I)),I=1,KATI),int(RSUMNU(JOPP)) 85 FORMAT(40X,I2,4X,12(I4,2X)/) 86 CONTINUE WRITE(6,87) (int(RSUMNT(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:', + //,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(RWW(JOPP,I)),I=1,KATI),int(RSUMSU(JOPP)) 89 CONTINUE WRITE(6,87) (int(RSUMST(I)),I=1,KATI),int(EMS) RETURN END C---------------------------------------- SUBROUTINE PRTOP C---------------------------------------- C C PRINT OUT OBSERVED 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 WRITE(6,1090) 1090 FORMAT(/,40X,'OBSERVED OPERATING POINTS FOR CONDITION X:') WRITE(6,1091) (FPFX(NNI-I+1),I=1,NNI) 1091 FORMAT(34X,' FPF: ',12(F6.4,1X)) WRITE(6,1092) (TPFX(NNI-I+1),I=1,NNI) 1092 FORMAT(34X,' TPF: ',12(F6.4,1X)) WRITE(6,1093) 1093 FORMAT(/,40X,'OBSERVED OPERATING POINTS FOR CONDITION Y:') WRITE(6,1091) (FPFY(NNJ-I+1),I=1,NNJ) WRITE(6,1092) (TPFY(NNJ-I+1),I=1,NNJ) 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 C FIXed array subscripts for UNIX machines 2/14/93 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(12),TPF(12) 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