C C :---------------------: C : : C : 'LABROC4' PROGRAM : Oct 1997 C : : C :---------------------: C C C******************************************************************************* C * C * C A PROGRAM FOR MAXIMUM LIKELIHOOD ESTIMATION OF A BINORMAL * C ROC CURVE AND ITS ASSOCIATED PARAMETERS FROM A SET OF * C CONTINUOUSLY-DISTRIBUTED DATA. FOR INHERENTLY CATEGORICAL * C DATA, THE PROGRAM ENTITLED 'ROCFIT' SHOULD BE USED INSTEAD. * C 'LABROC4' REPLACES THE EARLIER 'LABROC1'. * C * C 'LABROC4' IS BASED UPON AND EMPLOYS CODE FROM THE PROGRAM * C 'RSCORE II' BY: * C * C DONALD D. DORFMAN * C DEPARTMENT OF PSYCHOLOGY * C THE UNIVERSITY OF IOWA * C IOWA CITY, IOWA 52242 * C * C AS LISTED IN REFERENCE (1) BELOW. MODIFICATIONS TO 'RSCORE II' * C IMPLEMENTED IN 'LABROC4' WERE DESIGNED BY: * C * C CHARLES E. METZ, Ben Herman, JONG-HER SHEN, * C PU-LAN WANG, AND HELEN B. KRONMAN * C DEPARTMENT OF RADIOLOGY * C AND THE FRANKLIN MCLEAN MEMORIAL * C RESEARCH INSTITUTE -- MC2026 * C THE UNIVERSITY OF CHICAGO * C CHICAGO, ILLINOIS 60637 * C * C AND WERE WRITTEN BY Ben Herman, JONG-HER SHEN, * C PU-LAN WANG, AND HELEN B. KRONMAN. * C C Comments and requests should be directed to Dr. Metz C at the address above. C C Problems and/or bug reports should be directed to: C Ben Herman -- b-herman@uchicago.edu C or Dr. Metz -- c-metz@uchicago.edu C C C THESE MODIFICATIONS INCLUDE: * C (A). PROVIDING AUTOMATIC CATEGORIZATION OF THE CONTINUOUSLY * C DISTRIBUTED DATA. THE CATEGORIZATION ALGORITHM RANKS THE * C INPUT DATA VALUES AND THEN FORMS EFFECTIVE CATEGORIES * C CORRESPONDING TO RUNS OF CONSECUTIVE ACTUALLY NEGATIVE * C (OR ACTUALLY POSITIVE) CASES. * C (B). ALLOWING INPUT OF A FREE-TEXT LINE TO DESCRIBE THE * C DATA SET IN THE OUTPUT HEADING. * C (C). ALLOWING THE USE OF FREE FORMAT FOR DATA INPUT. * C (D). GIVING VARIABLES MORE MNEMONIC NAMES. * C (E). RENUMBERING STATEMENTS IN ORDER OF APPEARANCE. * C (F). TESTING THE CATEGORIZED DATA SET FOR 'DEGENERACY' * C BEFORE ATTEMPTING THE ITERATIVE CALCULATIONS. * C (G). CALCULATION OF 'DIFFSN' AND 'DIFFN' BETWEEN STATEMENT * C 7020 AND 7030 IN SUBROUTINE 'MOSLOP', AND USE OF THESE * C IN SUBSEQUENT COMPUTATIONS. * C (H). MODIFICATION OF THE ITERATIVE PROCEDURE TO ENSURE THAT * C TWO CONDITIONS ARE SATISFIED AT EACH STEP: (1) THE * C MATRIX OF SECOND PARTICAL DERIVATIVES MUST NOT BECOME * C NUMERICALLY SINGULAR, AND (2) THE LOG LIKELIHOOD * C FUNCTION MUST NOT DECREASE BY MORE THAN 10**(-4) OF * C ITS PREVIOUS ABSOLUTE VALUE WRT THE CURRENT STEP SIZE. * C IF EITHER OF THESE * C CONDITIONS IS VIOLATED AT SOME STEP, THE ADJUSTMENTS * C OF ALL PARAMETER ESTIMATES ARE REDUCED SUCCESSIVELY * C BY 50% UNTIL THE CONDITIONS ARE SATISFIED. ITERATION * C THEN PROCEEDS TO THE NEXT STEP. * C (I). MODIFICATION OF THE STOPPING CRITERION TO REFLECT THE * C NUMBER OF CATEGORIES CREATED FROM THE DATA (I.E., TO * C REFLECT THE NUMBER OF ADJUSTABLE PARAMETERS). * C (J). CALCULATION OF A NON-PARAMETRIC 'TRAPEZOIDAL RULE' * C ESTIMATE OF THE AREA UNDER THE ROC CURVE AND AN ESTIMATE * C OF ITS STANDARD DEVIATION, AS DESCRIBED IN REFERENCE (2). * C (K). CALCULATION AND OUTPUT OF 'TRUE POSITIVE FRACTIONS' ON * C THE FITTED BINORMAL ROC CURVE AT A VARIETY OF 'FALSE * C POSITIVE FRACTIONS', TO ALLOW THE FITTED ROC CURVE TO * C BE PLOTTED EASILY ON LINEAR AXES. * C (L). CALCULATION AND OUTPUT OF THE ASYMMETRIC 95% CONFIDENCE * C INTERVAL FOR 'TRUE POSITIVE FRACTION' AT EACH 'FALSE * C POSITIVE FRACTION' ON THE FITTED ROC. * C (M). CALCULATION AND OUTPUT OF AN ESTIMATED RELATIONSHIP * C BETWEEN THE CRITICAL TEST RESULT VALUE (WHICH SEPARATES * C NOMINALLY 'POSITIVE' RESULTS FROM NOMINALLY 'NEGATIVE' * C RESULTS) AND THE CORRESPONDING OPERATING POINT ON THE * C ESTIMATED ROC CURVE. * C * C * C REMARKS: * C (A). THIS PROGRAM ASSUMES THAT THE TRUE ROC CURVE PLOTS AS * C A STRAIGHT LINE ON "NORMAL-DEVIATE" AXES, OR EQUIVALENTLY, * C THAT THE INPUT DATA FOLLOW NORMAL DISTRIBUTIONS AFTER * C SOME UNKNOWN MONOTONIC TRANSFORMATION (SEE REFERENCE 3). * C ROC CURVES MEASURED IN A BROAD VARIETY OF FIELDS * C DEMONSTRATE THIS FORM (SEE REFERENCE 4). THE ASSUMPTION MAY * C BE SATISFIED EVEN WHEN THE RAW DATA HAVE MULTINODAL AND/OR * C SKEWED DISTRIBUTION. * C (B). THE BINORMAL ROC CURVE THAT LABROC4 FITS TO THE INPUT DATA * C IS NOT A MAXIMUM LIKELIHOOD ESTIMATE IN THE STRICT SENSE * C DESCRIBED IN REFERENCES (5) AND (6), BECAUSE THE INTERVAL * C BOUNDARIES USED TO PRODUCE THE EFFECTIVELY CATEGORICAL * C DATA EMPLOYED IN SEVERAL STEPS OF THE ANALYSIS ARE NOT * C FIXED IN ADVANCE, BUT INSTEAD ARE ESTABLISHED POST-HOC * C FROM THE CONTINUOUSLY-DISTRIBUTED INPUT DATA. HENCE, * C CATEGORY BOUNDARIES ON THE EFFECTIVE NORMAL-DEVIATE AXIS * C (I.E., THE Z(K) VARIABLES) ARE NOT REALLY 'PARAMETERS' OF * C THE MULTINOMIAL STATISTICAL MODEL (REFERENCES (5) AND (6)) * C UPON WHICH THE CURVE-FITTING PROCEDURE IS BASED. HOWEVER, * C THE AUTHORS' EXPERIENCE IN APPLYING LABROC4 TO SIMULATED * C DATA SUGGESTS THAT THIS THEORETICAL INCONSISTENCY CAN BE * C IGNORED IN PRACTICE, BECAUSE THE ACCURACIES OF THE ROC * C PARAMETER ESTIMATES ('A' AND 'B') AND OF THE STANDARD * C DEVIATION AND CORRELATION ESTIMATES PRODUCED BY LABROC4 * C ARE ESSENTIALLY THE SAME AS THE ACCURACIES OBTAINED WITH * C THE 'ROCFIT' AND 'RSCORE II' PROGRAMS WHEN THE LATTER * C PROGRAMS ARE APPLIED TO DATA IN WHICH THE CATEGORY * C BOUNDARIES WERE FIXED A PRIORI. * C (C). THIS PROGRAM ACCEPTS UP TO 4000 ACTUALLY NEGATIVE CASES AND * C UP TO 4000 ACTUALLY POSITIVE CASES. WITH LARGE DATA SETS * C (E.G., THOSE CONTAINING MORE THAN SEVERAL HUNDRED CASES IN * C TOTAL), EXECUTION OF 'LABROC4' MAY REQUIRE A SUBSTANTIAL * C AMOUNT OF COMPUTATION TIME, ESPECIALLY ON A MICROCOMPUTER. * C BECAUSE OF THIS, AND BECAUSE EXTREMELY LARGE DATA SETS CAN * C LEAD TO A LOSS OF NUMERICAL ACCURACY IN THE FINAL STAGES OF * C COMPUTATION, THE PROGRAM DEFAULTS TO AN EXPLICIT 20-CATEGORY * C ANALYSIS (RATHER THAN AN EFFECTIVELY CONTINUOUS ANALYSIS) WHEN * C THE NUMBER OF RUNS OF ACTUALLY NEGATIVE PLUS THE NUMBER OF * C RUNS OF ACTUALLY POSITIVE TRIALS IN THE RANKED INPUT DATA * C EXCEEDS 400. TYPICALLY, THIS OCCURS WHEN A TOTAL OF MORE * C THAN ABOUT 1300 CASES IS INPUT. INITIAL EXPERIENCE WITH * C RELATIVELY SMALL DATA SETS (CONTAINING A TOTAL OF ABOUT * C FIFTY TO A FEW HUNDRED CASES) IS RECOMMENDED. * C * C REFERENCES: * C (1). SWETS, J. A. AND PICKETT, R. M. - EVALUATION OF * C DIAGNOSTIC SYSTEMS: METHODS FROM SIGNAL DETECTION * C THEORY. (NEW YORK: ACADEMIC PRESS), 1982. * C (2). HANLEY, J.A. AND MCNEIL, B.J. - RADIOLOGY, 143:29-36, 1982. * C (3). METZ, C. E. - INVESTIGATIVE RADIOLOGY, 21:720-733, 1986. * C (4). SWETS, J. A. - PSYCHOLOGICAL BULLETIN, 99:181-198, 1986. * C (5). DORFMAN, D. D. AND ALF, E. - JOURNAL OF MATHEMATICAL * C PSYCHOLOGY, 6: 487-496, 1969. * C (6). GREY, D. R. AND MORGAN, B. J. T. - JOURNAL OF * C MATHEMATICAL PSYCHOLOGY, 9: 128-139, 1972. * C * C ACKNOWLEDGEMENT: * C DEVELOPMENT OF THIS PROGRAM AT THE UNIVERSITY OF CHICAGO * C WAS SUPPORTED BY THE U.S. DEPARTMENT OF ENERGY UNDER * C CONTRACTS DE-AC02-80EV10359 AND DE-AC02-82ER60033 AND * C UNDER GRANT DE-FG02-86ER60418. * C * C******************************************************************************* C C******************************************************************************* C * C REQUIRED INPUT AND DESCRIPTION OF OUTPUT: * C * C------------------------------------------------------------------------------* C * C * C FOR EACH DATA SET, INPUT THE FOLLOWING ON SUCCESSIVE LINES: * C 1. UP TO 80 CHARACTERS OF FREE TEXT TO DESCRIBE THE DATA SET. * C 2. AN ALPHABETIC CODE WORD INDICATING WHETHER SMALL OR LARGE * C TEST-RESULT VALUES ARE ASSOCIATED WITH ACTUALLY POSITIVE CASES: * C SMALL --> SMALL VALUES ARE ASSOCIATED WITH ACTUALLY POSITIVE * C CASES. * C LARGE --> LARGE VALUES ARE ASSOCIATED WITH ACTUALLY POSITIVE * C CASES. * C (NOTE: COMPUTER READS ONLY THE FIRST CHARACTER OF THE INPUT WORD.)* C 3. TWO SEQUENCES OF CONTINUOUSLY-DISTRIBUTED TEST-RESULT DATA. * C THE DATA FOR ACTUALLY NEGATIVE CASES SHOULD BE ENTERED FIRST, * C FOLLOWED BY THE DATA FOR ACTUALLY POSITIVE CASES. THE VALUES * C IN EACH SEQUENCE CAN BE INPUT IN ANY ORDER USING 'FREE * C FORMAT'; THE ONLY REQUIREMENTS ARE THAT THE VALUES MUST * C BE SEPARATED BY AT LEAST ONE BLANK SPACE, AND EACH * C SEQUENCE MUST BE TERMINATED BY AN ASTERISK (*). EACH OF THE * C TWO SEQUENCES CAN CONTAIN UP TO 4000 ENTRIES. * C (NOTE: THIS PROGRAM CAN ACCEPT UP TO 6 DIGITS TO THE LEFT OF THE * C DECIMAL POINT AND UP TO 8 DIGITS TO THE RIGHT OF THE * C DECIMAL POINT.) * C * C INPUT EXAMPLE-- * C ---------------------------------------------------------- * C SAMPLE INPUT * C LARGE * C 2.5 -3.01 +4.2 -0.21 * C 5.2 9.3 -4.3 -9.9 * C +3.3 4.3 * * C 0.221 4.232 0.453 * C +3.21 -2.12 -5.43 * C -3.2 * C 8.3 7.2 9.1 3.444 8.334 * C * * C --------------------------------------------------------- * C * C REMARK CONCERNING THE INPUT: * C -- THIS PROGRAM ALLOWS THE USE OF SEVERAL DATA SETS AS * C INPUT; SIMPLY ENTER THE DATA SETS SEQUENTIALLY, BEGINNING * C WITH THE FREE-TEXT DESCRIPTION. * C * C * C * C OUTPUT --- * C * C 1. THE FREE-TEXT DESCRIPTION OF THE DATA SET. * C 2. THE CONTINUOUSLY-DISTRIBUTED DATA SEQUENCES THAT WERE INPUT. * C 3. THE 'DIRECTION' OF THE TEST RESULTS AND THE NUMBERS * C OF ACTUALLY NEGATIVE AND ACTUALLY POSITIVE TRIALS IN THE INPUT. * C 4. OPERATING POINTS CORRESPONDING TO THE INPUT DATA. * C 5. FINAL ESTIMATES OF THE BINORMAL ROC CURVE PARAMETERS A AND B. * C THE PARAMETER A REPRESENTS THE 'Y-INTERCEPT' WHEREAS THE * C PARAMETER B REPRESENTS THE 'SLOPE' OF THE FITTED BINORMAL ROC * C CURVE WHEN IT IS PLOTTED ON NORMAL-DEVIATE AXES. * C 6. ESTIMATED STANDARD DEVIATION AND CORRELATION OF THE FINAL * C ESTIMATES OF THE BINORMAL ROC PARAMETERS A AND B. * C 7. ESTIMATE AND ESTIMATED STANDARD DEVIATION OF THE AREA * C (OFTEN DENOTED BY 'A SUB Z') BENEATH THE FITTED BINORMAL * C ROC CURVE WHEN IT IS PLOTTED ON LINEAR AXES. * C 8. ESTIMATE AND ESTIMATED STANDARD DEVIATION OF THE "TRAPEZOIDAL * C RULE" AREA BENEATH THE ROC, WHICH IS FORMED BY USING LINE * C SEGMENTS TO CONNECT THE OPERATING POINTS ASSOCIATED * C WITH THE INPUT DATA (SEE REFERENCE 2). * C 9. TRUE-POSITIVE FRACTION AS A FUNCTION OF FALSE-POSITIVE * C FRACTION ON THE FITTED BINORMAL ROC CURVE, WITH LOWER AND * C UPPER BOUNDS OF THE ASYMMETRIC 95% CONFIDENCE INTERVAL * C ESTIMATED BY ASSUMING THAT THE ROC PARAMETER ESTIMATES 'A' * C AND 'B' ARE JOINTLY NORMAL. * C 10. AN ESTIMATED RELATIONSHIP BETWEEN THE CRITICAL TEST-RESULT * C VALUE (WHICH SEPARATES NOMINALLY 'POSITIVE' FROM NOMINALLY * C 'NEGATIVE' RESULTS) AND THE CORRESPONDING OPERATING POINT * C ON THE FITTED ROC CURVE. * C * C REMARK CONCERNING THE OUTPUT: * C -- THIS PROGRAM WILL CHECK FOR DEGENERACY IN THE CATEGORIZED DATA. * C DEGENERATE DATA SETS ARE THOSE FOR WHICH AN EXACT FIT TO * C THE DATA IS PROVIDED BY A BINORMAL ROC CURVE AND ITS * C ASSOCIATED CUT-OFFS SUCH THAT ONE OR MORE OF THE PARAMETERS * C TAKES ON AN INFINITE VALUE, SO THAT THE ITERATIVE CALCULATIONS * C CANNOT CONVERGE. IF THE DATA SET IS DEGENERATE, THIS * C PROGRAM WILL OUTPUT A MESSAGE DESCRIBING THE KIND OF * C DEGENERACY FOUND AND THEN ABORT EXECUTION OF THAT DATA SET. * C IN GENERAL, DEGENERACY SHOULD BE FOUND ONLY IN VERY SMALL DATA * C SETS. * C April 1995: c Added support for 4000 act + and 4000 act - cases C Nov. 1996: C Added support for MLE estimation in the Critical value table C Installed better Character to Number conversion in TONUM C March 1997: C Use Labroc5 to compute the initial estimates C Clean up output for 80 column screens rather than C The old 132 column DecWriter II (finally) C Oct. 1997: C Added estimation correction to MOSLOP -- ecor fixes a C looping problem effectively damping the damped newtons method C when the method bounces between 2 values. C * C * C******************************************************************************* COMMON/B1/MN,MS,KATBIG COMMON/B2/ACTN(4000),ACTSN(4000),V(2,8000) COMMON/B3/CATN4(400),CATS4(400),KAT4 character*80 ititle COMMON/B5/TRAPA,VTRAPA,ititle,ZZ(399) COMMON/E/IEND,IERROR COMMON/DECIML/MAXDEC real aeff,beff,mun,sun,mup,sup IEND=0 C C READ INPUT DATA C 5 IERROR=0 C read(5,10,end=101) Iseed 10 format(I10) CALL READHL(IKEY,ititle) IF(IEND.EQ.1.OR.IERROR.GT.0)GO TO 101 MAXDEC=0 CALL READIN(MN,ACTN) IF(IERROR.GT.0)GO TO 101 CALL READIN(MS,ACTSN) IF(IERROR.GT.0)GO TO 101 C C PRINT TITLE AND ECHO PRINT TEST RESULTS C WRITE(6,13) 13 FORMAT(//28X,'LABROC4 (Oct 1997 VERSION): ' + //16X,'MAXIMUM LIKELIHOOD ESTIMATION OF A BINORMAL ROC' + /16X,'CURVE FROM CONTINUOUSLY-DISTRIBUTED TEST RESULTS') WRITE(6,25)MN 25 FORMAT(//// +14X, '-----------------------------------------------------'/ +14X, ' ORIGINAL INPUT OF ',I4,' ACTUALLY NEGATIVE CASES'/ +14X, '-----------------------------------------------------'/) CALL PPRINT(MN,ACTN) WRITE(6,30)MS 30 FORMAT(/// +14X, '----------------------------------------------------'/ +14X, ' ORIGINAL INPUT OF ',I4,' ACTUALLY POSITIVE CASES'/ +14X, '----------------------------------------------------'/) CALL PPRINT(MS,ACTSN) C CHANGE SIGN FOR INPUT DATA IF SMALLER VALUES ARE ASSOCIATED C WITH TRUE ABNORMALITY C IF(IKEY.EQ.0)CALL CVSIGN() C C COLLAPSE QUANTITATIVE DATA BASE ON THEIR "RUNS" C KATBIG=0 CALL CLAPSE C C CALCULATE TRAPEZOIDAL(WILCOXON) AREA C CALL TRAREA() C C USE 'LABROC5' TO GET INITIAL ESTIMATES FOR A AND B C CALL LABROC1(IKEY,A0,B0) IF(IERROR.EQ.1.OR.KATBIG.EQ.1)GO TO 5 C C GET FINAL ESTIMATES C CALL LSTFIT(IKEY,A0,B0) GO TO 5 101 STOP END C--------------------------------------- SUBROUTINE READHL(IKEY,ititle) C--------------------------------------- C C THIS SUBROUTINE READS IN A FREE-TEXT DESCRIPTION OF THE C DATA AND AN ALPHABETIC CODE WORD WHICH INDICATES WHETHER C LARGE VALUES OR SMALL VALUES ARE ASSOCIATED WITH ACTUAL C ABNORMALITY C COMMON/E/IEND,IERROR character line,lstring COMMON/PASS/LSTRING(80),LINE(80),LENGTH character*80 ititle C ititle(:)=' ' READ(5,110,END=300)(LINE(I),I=1,80) 110 FORMAT(80A1) DO 120 I=1,80 IF(LINE(I).EQ.' ')GO TO 120 IPNT=I GO TO 130 120 CONTINUE WRITE(6,125) 125 FORMAT(2X,'--- NO INPUT FOR DATA DESCRIPTION ---') IPNT=1 130 DO 150 I=IPNT,80 if (ichar(line(i)).ge.32.and.ichar(line(i)).lt.128) then ititle(I-IPNT+1:I-IPNT+1)=LINE(I) endif 150 CONTINUE C 200 READ(5,110,end=300)(LINE(I),I=1,80) CALL GETWRD(IERROR) IF(IERROR.GT.0)GO TO 210 IKEY=99 IF(LSTRING(1).EQ.'s'.OR. LSTRING(1).EQ.'S'.or. + lstring(1).eq.'0')IKEY=0 IF(LSTRING(1).EQ.'l'.OR. LSTRING(1).EQ.'L'.or. + lstring(1).eq.'1')IKEY=1 IF(IKEY.EQ.1 .OR. IKEY.EQ.0)RETURN 210 Line(1)=line(1) WRITE(6,215)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 215 FORMAT(2X,80A1/2X,'*** INVALID KEY CODE DETECTED ***') IERROR=1 RETURN 300 IEND=1 RETURN END C---------------------------------- SUBROUTINE READIN(NUM,F1) C---------------------------------- C C THIS SUBROUTINE READS IN A SEQUENCE OF INPUT DATA IN FREE C FORMAT. THE ONLY FORMAT REQUIREMENTS ARE THAT (1) C ANY TWO INPUT VALUES MUST BE SEPARATED BY AT LEAST ONE C IBLANK COLUMN AND (2) THE INPUT DATA SEQUENCE MUST BE TERMINATED C BY AN ASTERISK (*). C NEGATIVE VALUES ARE ACCEPTABLE. C COMMON/E/IEND,IERROR character line,lstring COMMON/PASS/LSTRING(80),LINE(80),LENGTH DIMENSION F1(*) C NUM=0 100 READ(5,110)(LINE(I),I=1,80) 110 FORMAT(80A1) 120 CALL GETWRD(NOMORE) IF(NOMORE.GT.0)GO TO 100 IF(LSTRING(1).EQ.'*')GO TO 130 CALL TONUM(LENGTH,LSTRING,RVALUE,IERROR) IF(IERROR.GT.0)GO TO 140 NUM=NUM+1 F1(NUM)=RVALUE GO TO 120 130 IERROR=0 RETURN 140 line(1)=line(1) WRITE(6,135)(LSTRING(I),I=1,LENGTH),(LINE(I),I=1,80) 135 FORMAT(2X,80A1/2X,'*** INVALID TEST RESULT DETECTED ***') RETURN END C--------------------------------------- SUBROUTINE GETWRD(NOMORE) C--------------------------------------- C C GET A STRING FROM THE INPUT LINE C character line,lstring COMMON/PASS/LSTRING(80),LINE(80),LENGTH NOMORE=0 LENGTH=0 I=1 C C SKIP LEADING BLANKS C DO 10 J=1,80 LSTRING(J)=' ' 10 CONTINUE C 20 IF(I .GT. 80) GO TO 90 IF(LINE(I) .NE. ' ') GO TO 40 I= I+1 GO TO 20 C C GET THE FIRST WORD OF THE CURRENT LINE, AND COPY IT C INTO "LSTRING" C 40 II= 1 50 LSTRING(II)= LINE(I) LENGTH=II I= I+1 IF(I .GT. 80) GO TO 60 IF(LINE(I).EQ.' ') GO TO 70 II= II+1 GO TO 50 C C SHIFT THE CURRENT LINE BY DELETING THE FIRST WORD C 60 DO 65 J=1,80 LINE(J)=' ' 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)=' ' 85 CONTINUE RETURN C C NO MORE WORD IN CURRENT LINE C 90 NOMORE=1 RETURN END C----------------------------------------------------------- SUBROUTINE TONUM(LEN,LINPUT,RNUM,IERR) C----------------------------------------------------------- C C CONVERT CHARACTER STRING TO REAL NUMBER C character linput DIMENSION LINPUT(80) character*80 buff,temp COMMON/DECIML/MAXDEC rnum=0 buff(:)=' ' iu=0 2 format(1x,80a1) 5 format(1x,g12.5) write(buff,2) (linput(in),in=1,len) ie=index(buff,'e') if (ie.eq.0) ie=index(buff,'E') id=index(buff,'.') if (id.eq.0.and.ie.eq.0) then buff(len+2:)='.' else if (id.eq.0) then temp=buff(ie:len+1) buff(ie:ie)='.' buff(ie+1:len+2)=temp iu=ie-id else if (ie.eq.0) then iu=1+len-id endif read(buff,5,err=120,iostat=ierr) rnum if (maxdec.lt.iu) maxdec=iu return 120 IERR= 1 200 RETURN END C----------------------------------- SUBROUTINE PPRINT(NUM, F) C----------------------------------- C C PRINT OUT THE VECTOR "NUM" C COMMON/DECIML/MAXDEC DIMENSION F(*) IF(MAXDEC.EQ.0)WRITE(6,10)(F(J),J=1,NUM) IF(MAXDEC.EQ.1)WRITE(6,11)(F(J),J=1,NUM) IF(MAXDEC.EQ.2)WRITE(6,12)(F(J),J=1,NUM) IF(MAXDEC.EQ.3)WRITE(6,13)(F(J),J=1,NUM) IF(MAXDEC.EQ.4)WRITE(6,14)(F(J),J=1,NUM) IF(MAXDEC.EQ.5)WRITE(6,15)(F(J),J=1,NUM) IF(MAXDEC.EQ.6)WRITE(6,16)(F(J),J=1,NUM) IF(MAXDEC.EQ.7)WRITE(6,17)(F(J),J=1,NUM) IF(MAXDEC.GE.8)WRITE(6,18)(F(J),J=1,NUM) 10 FORMAT(5(F8.0,2X)) 11 FORMAT(5(F9.1,2X)) 12 FORMAT(5(F10.2,2X)) 13 FORMAT(5(F11.3,2X)) 14 FORMAT(5(F12.4,2x)) 15 FORMAT(5(f13.5,2x)) 16 FORMAT(5(f14.6,2x)) 17 FORMAT(5(f15.7,2x)) 18 FORMAT(5(f16.8,2x)) RETURN END C---------------------------------- SUBROUTINE CVSIGN() C---------------------------------- C C CONVERTS ALL INPUT DATA WHICH HAS THE SMALLER VALUES C REPRESENT STRONGER EVIDENCE OF ABNORMALITY TO LARGER VALUES C REPRESENT STRONGER EVIDENCE OF ABNORMALITY C COMMON/B1/MN,MS,KATBIG COMMON/B2/ACTN(4000),ACTSN(4000),V(2,8000) DO 1010 I=1,MN ACTN(I)=-ACTN(I) 1010 CONTINUE DO 1020 I=1,MS ACTSN(I)=-ACTSN(I) 1020 CONTINUE RETURN END C-------------------------------- SUBROUTINE CLAPSE C-------------------------------- C C POOL ALL NOISE AND S+N CASES INTO A DATA MATRIX C BY RANKING THE INPUT TEST-RESULT VALUES AND THEN C CREATING CATEGORIES SUCH THAT EACH CONTAINS A RUN C OF 'N' CASES OR A RUN OF 'S+N' CASES AND NO CASES OF C THE OTHER KIND C COMMON/B1/MN,MS,KATBIG COMMON/B2/ACTN(4000),ACTSN(4000),V(2,8000) COMMON/B3/CATN4(400),CATS4(400),KAT4 character*80 ititle COMMON/B5/TRAPA,VTRAPA,ititle,ZZ(399) integer vt,vcount(8000) common/vc/vcount,vt DIMENSION CATN(400),CATS(400) NUM=MN+MS DO 100 I=1,NUM IF(I.LE.MN)THEN V(1,I)=ACTN(I) V(2,I)=-1 ELSE V(1,I)=ACTSN(I-MN) V(2,I)=1 ENDIF 100 CONTINUE C C SORT DATA MATRIX C CALL SSORT(NUM,V) C C COUNT TOTAL RUNS C KAT4=0 COUNTN=0 COUNTS=0 vcount(1)=0 vt=1 DO 2100 J=1,NUM IF(J.GT.1)THEN IF(V(1,J).NE.V(1,J-1))THEN IF(((J.GT.2).AND.(V(1,J-1).EQ.V(1,J-2))) .OR. + ((J.LT.NUM).AND.(V(1,J).EQ.V(1,J+1))) .OR. + (V(2,J).NE.V(2,J-1)))THEN KAT4=KAT4+1 IF(KAT4.Ge.400)THEN KATBIG=1 RETURN ENDIF CATN4(KAT4)=COUNTN CATS4(KAT4)=COUNTS ZZ(KAT4)=(V(1,J)+V(1,J-1))/2 COUNTN=0 COUNTS=0 ENDIF vt=vt+1 vcount(vt)=0 ENDIF ENDIF vcount(vt)=vcount(vt)+1 IF(V(2,J).LT.0)THEN COUNTN=COUNTN+1 ELSE COUNTS=COUNTS+1 ENDIF 2100 CONTINUE C C ASSIGN FOR THE LAST CATEGORY C KAT4=KAT4+1 IF(KAT4.GT.400)THEN KATBIG=1 RETURN ENDIF CATN4(KAT4)=COUNTN CATS4(KAT4)=COUNTS C C ELIMINATE CONSECUTIVE ZEROS C C --- NEGATIVE CASES C TOTALN=0 TOTALS=0 NKAT=0 DO 4010 J=1,KAT4-1 TOTALN=TOTALN+CATN4(J) TOTALS=TOTALS+CATS4(J) IF(CATN4(J).NE.0 .OR. CATN4(J+1).NE.0)THEN NKAT=NKAT+1 CATN(NKAT)=TOTALN CATS(NKAT)=TOTALS TOTALN=0 TOTALS=0 ENDIF 4010 CONTINUE C C LAST CATEGORY C NKAT=NKAT+1 IF(CATN4(KAT4).EQ.0)THEN CATN(NKAT)=TOTALN+CATN4(KAT4) CATS(NKAT)=TOTALS+CATS4(KAT4) ELSE CATN(NKAT)=CATN4(KAT4) CATS(NKAT)=CATS4(KAT4) ENDIF C C --- POSITIVE CASES C TOTALN=0 TOTALS=0 KAT4=0 DO 4020 J=1,400 CATN4(J)=0 CATS4(J)=0 4020 CONTINUE DO 4030 J=1,NKAT-1 TOTALN=TOTALN+CATN(J) TOTALS=TOTALS+CATS(J) IF(CATS(J).NE.0 .OR. CATS(J+1).NE.0)THEN KAT4=KAT4+1 CATN4(KAT4)=TOTALN CATS4(KAT4)=TOTALS TOTALN=0 TOTALS=0 ENDIF 4030 CONTINUE C C LAST CATEGORY C KAT4=KAT4+1 IF(CATS(NKAT).EQ.0)THEN CATN4(KAT4)=TOTALN+CATN(NKAT) CATS4(KAT4)=TOTALS+CATS(NKAT) ELSE CATN4(KAT4)=CATN(NKAT) CATS4(KAT4)=CATS(NKAT) ENDIF RETURN END C----------------------------------- SUBROUTINE TRAREA() C----------------------------------- C C THIS SUBROUTINE COMPUTES THE AREA UNDER THE C ROC CURVE BY USING TRAPEZOIDAL(WILCOXON) METHOD C COMMON/B1/MN,MS,KATBIG COMMON/B2/ACTN(4000),ACTSN(4000),V(2,8000) character*80 ititle COMMON/B5/TRAPA,VTRAPA,ititle,ZZ(399) C TW=0. DO 1100 I=1,MN DO 1100 J=1,MS IF(ACTN(I).LT.ACTSN(J))TW=TW+1. IF(ACTN(I).EQ.ACTSN(J))TW=TW+0.5 1100 CONTINUE TRAPA=TW/FLOAT(MN*MS) Q1=TRAPA/(2-TRAPA) Q2=(2*TRAPA**2)/(1+TRAPA) VTRAPA=SQRT((TRAPA*(1-TRAPA)+(MS-1)*(Q1-TRAPA**2)+(MN-1)* + (Q2-TRAPA**2))/(MN*MS)) RETURN END C----------------------------------- SUBROUTINE SSORT(NUM,F2) C----------------------------------- C C SORT A 2D VECTOR IN ACCENDING ORDER C DIMENSION F2(2,4000) C NN = NUM-1 DO 30 I= 1, NN SMALL= F2(1,I) JJ= I C DO 20 J= I,NUM IF (F2(1,J) .GE. SMALL) GO TO 20 SMALL = F2(1,J) JJ = J 20 CONTINUE C IF (JJ .EQ. I) GO TO 30 TEMP = F2(1,I) F2(1,I)= F2(1,JJ) F2(1,JJ)= TEMP TEMP = F2(2,I) F2(2,I)= F2(2,JJ) F2(2,JJ)= TEMP 30 CONTINUE RETURN END C----------------------------------- SUBROUTINE SSORT1(NUM,F1) C----------------------------------- C C SORT A VECTOR IN ACCENDING ORDER C DIMENSION F1(*) C NN = NUM-1 DO 30 I= 1, NN SMALL= F1(I) JJ= I C DO 20 J= I,NUM IF (F1(J) .GE. SMALL) GO TO 20 SMALL = F1(J) JJ = J 20 CONTINUE C IF (JJ .EQ. I) GO TO 30 TEMP = F1(I) F1(I)= F1(JJ) F1(JJ)= TEMP 30 CONTINUE RETURN END C------------------------------------------- SUBROUTINE LABROC1(IKEY,A0,B0) C------------------------------------------- COMMON/B1/MN,MS,KATBIG COMMON/B2/ACTN(4000),ACTSN(4000),V(2,8000) COMMON/B3/CATN4(400),CATS4(400),KAT4 character*80 ititle COMMON/B5/TRAPA,VTRAPA,ititle,ZZ(399) COMMON/C1/CATN1(21),CATS1(21),KAT1 real aeff,beff,mun,sun,mup,sup C C COMPUTES INITIAL ESTIMATES OF A AND B C C COLLAPSE DATA INTO 'KAT1' CATEGORIES C (KAT1=20 FOR LARGE DATA SETS OR 10 FOR SMALL DATA SETS) C CALL CATGRZ() kat=kat4 c if (katbig.ne.0) kat=-kat1 c WRITE(6,2495)ititle,kat 2495 FORMAT(7X,A60,I4) c write(*,*) Trapa,vtrapa c write(*,*) aeff, beff C C USE METHOD OF SCORING TO ESITMATE A AND B C FROM THE COLLAPSED DATA C CALL INIFIT(A0,B0,IKEY) RETURN END C---------------------------------------------------- SUBROUTINE CATGRZ() C SUBROUTINE CATGRZ(IKEY) C---------------------------------------------------- C C CATEGORIZE QUANTITATIVE DATA C DIMENSION ZTEMP(19) COMMON/B1/MN,MS,KATBIG COMMON/B2/ACTN(4000),ACTSN(4000),V(2,8000) character*80 ititle COMMON/B5/TRAPA,VTRAPA,ititle,ZZ(399) COMMON/C1/CATN1(21),CATS1(21),KAT1 COMMON/CAT/SUMF(2,7999),Z(7999) C NUM=MN+MS C C CHECK HOW MANY OPERATING POINTS AVAILABLE FOR THE QUANTITATIVE C DATA C CALL RANKS(NRANK) C C IF "NRANK" CUTOFFS ARE AVAILABLE, AND C (A) NRANK <= 9 :: USE ALL OF THE OPERATING POINTS AVAILABLE C (B) 9 < NRANK <= 19 & MIN(MN,MS) < 40 :: COLLAPSE INTO 10 C CATEGORIES C (C) 9 < NRANK <= 19 & MIN(MN,MS) >= 40 :: USE ALL OF THE C OPERATING POINTS AVAILABLE C (D) NRANK > 19 & MIN(MN,MS) < 40 :: SAME AS (B) C (E) NRANK > 19 & MIN(MN,MS) >= 40 :: COLLAPSE INTO 20 CATEGORIES C IF(NRANK .LT. 10) GO TO 1100 IF(NRANK .GE. 20) GO TO 1110 IF(MN.LT.40 .OR. MS.LT.40)GO TO 1110 1100 DO 1105 I=1,NRANK SUMF(2,I)=1 1105 CONTINUE GO TO 1120 1110 KAT1=10 IF(MN.GE.40 .AND. MS.GE.40)KAT1=20 if (mn.lt.ms) then do inx=1,nrank/2 temp=sumf(1,inx) sumf(1,inx)=sumf(1,(nrank+1)-inx) sumf(1,(nrank+1)-inx)=temp enddo do inx=1,nrank sumf(1,inx)=2-sumf(1,inx) enddo endif kk=kat1-1 CALL SETCUT(NRANK,Kk) if (mn.lt.ms) then do inx=1,nrank/2 temp=sumf(2,inx) sumf(2,inx)=sumf(2,(nrank+1)-inx) sumf(2,(nrank+1)-inx)=temp temp=sumf(1,inx) sumf(1,inx)=sumf(1,(nrank+1)-inx) sumf(1,(nrank+1)-inx)=temp enddo do inx=1,nrank sumf(1,inx)=2-sumf(1,inx) enddo endif C C USE ALL THE OPERATING POINTS AVAILABLE C 1120 II=0 DO 1140 J=1,NRANK IF(SUMF(2,NRANK-J+1) .LT. 0)GO TO 1140 II=II+1 ZTEMP(II)=Z(NRANK-J+1) 1140 CONTINUE KAT1=II+1 C C CONVERT QUANTITATIVE VALUE INTO CATEGORICAL DATA C CALL CONVRT(MN,ACTN,KAT1-1,ZTEMP,CATN1) CALL CONVRT(MS,ACTSN,KAT1-1,ZTEMP,CATS1) C C ASSIGN VARIABLES C IF(KATBIG.EQ.0)RETURN DO 1160 I=2,KAT1 ZZ(I-1)=ZTEMP(I-1) 1160 CONTINUE RETURN END C------------------------------------- SUBROUTINE RANKS(NRANK) C------------------------------------- C C COUNT TOTAL NUMBER OF OPERATING POINTS FOR THE QUANTITATIVE C DATA; COUNT ONE ONLY FOR ANY TIE VALUES C COMMON/B1/MN,MS,KATBIG COMMON/B2/ACTN(4000),ACTSN(4000),V(2,8000) COMMON/CAT/SUMF(2,7999),Z(7999) integer ofpf,otpf,recat,chng integer i,num double precision rmn,rms real fpf,tpf C C INITIALIZATION C DO 1005 I=1,3999 SUMF(1,I)=0. SUMF(2,I)=-1 1005 CONTINUE NUM=MN+MS recat=1 chng=0 rmn=0 rms=0 ofpf=0 otpf=2 nrank=1 if (recat.eq.1) then fpf=0. tpf=0. else fpf=1.0e-7 tpf=1.0e-7 endif do 1030 i=2,num if(v(2,num-i+2) .gt. 0)go to 1010 rmn=rmn+1.0D0 fpf=REAL(rmn/mn) ofpf=-1 go to 1020 1010 rms=rms+1.0D0 tpf=REAL(rms/ms) ofpf=1 1020 if(v(1,num-i+1).eq.v(1,num-i+2)) then if(v(2,num-i+1).ne.v(2,num-i+2)) then chng=1 otpf=2 endif c if (i.eq.num) go to 1025 go to 1030 endif if (chng.eq.1) ofpf=0 if (fpf.gt.0.and.tpf.gt.0.and.ofpf.ne.otpf) then if (sumf(1,nrank).ne.0) nrank=nrank+1 chng=0 endif 1025 if (recat.eq.1) otpf=ofpf sumf(1,nrank)=fpf+tpf z(nrank)=(v(1,num-i+1)+v(1,num-i+2))/2. if (recat.eq.1.and. + (tpf.gt.0.9999.or.fpf.gt.0.9999)) go to 2000 1030 continue 2000 return end C----------------------------------------- SUBROUTINE SETCUT(NRANK,KK) C----------------------------------------- C C COLLAPSE QUANTITATIVE DATA INTO (KK+1) CATEGORIES C C COMMON/B2/ACTN(4000),ACTSN(4000),V(2,8000) COMMON/CAT/SUMF(2,7999),Z(7999) real h,d,div,lastrank,sumz(8000) integer nn,i,j,numcut,istart real den,step,target c ----------------------------------------------- c Initialize variables c ----------------------------------------------- do i=1,8000 sumz(i)=0.0 enddo numcut=0 lastrank=0 j=1 c ----------------------------------------------- c count total number of distinct target operating points c ----------------------------------------------- do i=1,nrank sf=sumf(1,i)-1.0D0 if (sf.gt.0) then call ndtr((1.96*sf/(sumf(1,nrank)-1.0)),sumz(i),den) else call ndtr((1.96*sf/(1.0-sumf(1,1))),sumz(i),den) endif enddo if (sumf(1,1).lt.0.15) then h=0 else h=sumz(1) endif if (sumf(1,nrank).gt.1.85) then d=1 else d=sumz(nrank) endif if (sumf(1,1).eq.0) j=j+1 div=kk if (h.ne.0 .and. d.ne.1) div=div-1 if (h.eq.0 .and. d.eq.1) div=div+1 step = (d-h)/div c ----------------------------------------------- c note div here now represents how many points to grab c either starting from this location or the next c ----------------------------------------------- div=kk-1 istart=0 if (h.eq.0) then istart=1 div=kk endif do 1030 i=istart,INT(div) target=h+i*step 1010 if(abs(sumz(j+1)-target) .ge. abs(sumz(j)-target))then if ((sumz(j)-lastrank).ge.(2*step/div)) goto 1020 endif j=j+1 if (j.gt.nrank) goto 1031 go to 1010 1020 if(sumf(2,j) .gt. 0)go to 1030 sumf(2,j)=1 lastrank=sumz(j) numcut=numcut+1 1030 continue c ----------------------------------------------- c if we already have enough operating points, then return; c otherwise, insert as many operating points as we need (= kk-numcut), c if we have enough candidates or we use up all of the c available candidates, whichever comes first. c ----------------------------------------------- 1031 nn=kk-numcut if(numcut .eq. kk) go to 1050 if((nrank-numcut) .lt. nn)nn=nrank-numcut do 1040 i=1,nn call insert(nrank,sumf,sumz) 1040 continue c 1050 kk=numcut+nn return end C-------------------------------------------------- SUBROUTINE CONVRT(N,DATA,K,CUTOFF,CAT) C-------------------------------------------------- C C COLLAPSE DATA INTO "K+1" CATEGORIES C DIMENSION DATA(*),CUTOFF(*),CAT(*) DO 1000 I=1,21 CAT(I)=0 1000 CONTINUE C DO 1040 I=1,N IF(DATA(I) .GT. CUTOFF(1)) GO TO 1010 CAT(1)=CAT(1)+1 GO TO 1040 1010 DO 1030 J=1,K-1 IF(.NOT.(DATA(I).GT.CUTOFF(J).AND.DATA(I).LE.CUTOFF(J+1))) + GO TO 1030 CAT(J+1)=CAT(J+1)+1 GO TO 1040 1030 CONTINUE IF(DATA(I) .GT. CUTOFF(K)) CAT(K+1)=CAT(K+1)+1 1040 CONTINUE RETURN END C----------------------------------------------- SUBROUTINE INSERT(NRANK,SUMF,sumz) C----------------------------------------------- C C SELECT THE CANDIDATE FROM THE MAXIMUM (delta FPF+delta TPF) C WHICH HAS MORE THAN ONE POINT IN BETWEEN C real sumz(*) DIMENSION SUMF(2,7999) C C COUNT TOTAL NUMBER OF TARGET OPERATING POINTS WHICH HAS C FPF+TPF < 1.0 C integer mto real delta,diff,dist integer i1,i2,idelta,idist,j,j1,j2,jump,i,nbound,numcut integer minn,max1,max2,less max1=0 max2=0 idelta=0 numcut=0 less=0 mto=1 c mto = more than one (in this case two or more points in between) C LAbroc5 Requires MTO to be 1 for best performance. c c count total number of target operating points which has c fpf+tpf < 1.0 c do 1010 i=1,nrank if(sumf(2,i) .lt. 0)go to 1010 if(sumf(1,i) .lt. 1.0) less=less+1 numcut=numcut+1 1010 continue c c if we have less than one half of the target operating points c with fpf+tpf < 1 start inserting from the signal end; otherwise, c start from the noise end c 1012 if(less .gt. (numcut/2))go to 1020 i1=1 jump=1 nbound=nrank go to 2010 1020 i1=nrank jump=-1 nbound=1 c c check for two consecutive target operating points c 2010 if(sumf(2,i1) .gt. 0) go to 2020 i1=i1+jump go to 2010 2020 i2=i1 2030 i2=i2+jump if(i2 .eq. (nbound+jump))go to 2060 if(sumf(2,i2) .lt. 0)go to 2030 c c round off to 5 significant digits c dist=abs(sumz(i1)-sumz(i2))*1.0e+05 idist=int(dist) if((dist-float(idist)) .ge. 0.5)idist=idist+1 c c select the maximum (delta fpf + delta tpf) which also has c more than one point in between c if(idist .le. idelta) go to 2050 if (abs(i2-i1).le.mto) goto 2050 idelta=idist max1=i1 max2=i2 2050 i1=i2 go to 2030 c c choose the most middle point as the new target operating point c 2060 if (mto.eq.2.and.idelta.eq.0) then mto=1 goto 1012 endif if (idelta.eq.0) return j1=max1+jump j2=max2-jump if(j1 .le. j2) go to 2070 j1=max2-jump j2=max1+jump 2070 delta=abs(sumz(max1)-sumz(max2)) minn=int((j1+j2)/2) do 2080 j=j1,j2 diff=abs(abs(sumz(j)-sumz(max1)) - + abs(sumz(j)-sumz(max2))) if(diff .ge. delta)go to 2080 delta=diff minn=j 2080 continue sumf(2,minn)=1 return end C--------------------------------------------- SUBROUTINE INIFIT(A0,B0,IKEY) C--------------------------------------------- C C USE METHOD OF SCORING TO ESTIMATE A AND B. C INTEGER IARY(401) DIMENSION X(20),XS(20),FPF(401),TPF(401) COMMON/E/IEND,IERROR COMMON/B1/MN,MS,KATBIG character*80 ititle COMMON/B5/TRAPA,VTRAPA,ititle,ZZ(399) COMMON/C1/CATN1(21),CATS1(21),KAT1 COMMON/F1/CATN(400),CATS(400),KAT C COMMON/F2/DIFFN(400),DIFFSN(400) KAT=KAT1 EMN=MN EMS=MS NN=KAT+1 KK=KAT-1 DO 1010 I=1,NN CATN(I)=CATN1(I) CATS(I)=CATS1(I) 1010 CONTINUE C C IF THE INPUT TEST RESULTS HAVE MORE THAN 400 RUNS, C THE ESTIMATES FROM THIS SUBROUTINE WILL BE USED AS C FINAL ESTIMATES C C WRITE OUT HEADING, DATA DESCRIPTION, NUMBERS OF C CATEGORIES AND CASES, AND RESPONSE DATA FROM INPUT. C IF(KATBIG.EQ.0)GO TO 1060 WRITE(6,1020) 1020 FORMAT(/////5X,'DUE TO THE LARGE AMOUNT OF INPUT DATA',/5x, + '(WHICH INVOLVES MORE THAN 400 RUNS OF',/5x, + 'EXCLUSIVELY NEGATIVE OR EXCLUSIVELY POSITIVE CASES),'/5x, + 'LABROC4 WILL NOT ATTEMPT AN EFFECTIVELY CONTINUOUS ANALYSIS.' + ,/5X,'INSTEAD, THE ANALYSIS WILL BE PERFORMED USING', + ' 20 POST-HOC CATEGORIES.',/5x,'(I.E., LABROC4 WILL DEFAULT', + ' TO LABROC5).') WRITE(6,1030) ititle(1:80) c WRITE(6,1030)ititle,kat c 1030 FORMAT(7X,A80,I4) 1030 FORMAT(//17X,'LABROC4/5 : ',A////) IF(IKEY.EQ.1)WRITE(6,1040) 1040 FORMAT(10X,'LARGER VALUES OF THE TEST RESULT REPRESENT STRONGER', + ' EVIDENCE',/10x,' THAT THE CASE IS ACTUALLY POSITIVE'/ + 10X,'(E.G., THAT THE PATIENT IS ACTUALLY ABNORMAL)') IF(IKEY.EQ.0)WRITE(6,1041) 1041 FORMAT(10X,'SMALLER VALUES OF THE TEST RESULT REPRESENT STRONGER', + ' EVIDENCE',/10x,' THAT THE CASE IS ACTUALLY POSITIVE'/ + 10X,'(E.G., THAT THE PATIENT IS ACTUALLY ABNORMAL)') WRITE(6,1050)EMN,EMS 1050 FORMAT(10X,'NO. OF ACTUALLY NEGATIVE CASES=',F6.0,/10X, 1 'NO. OF ACTUALLY POSITIVE CASES=',F6.0) C CHECK INPUT DATA. C 1060 SUMR=0. SUMRP=0. DO 1100 I=1,KAT SUMR=SUMR+CATN(I) SUMRP=SUMRP+CATS(I) 1100 CONTINUE C C CALCULATE OBSERVED OPERATING POINTS C FPF(NN)=0. TPF(NN)=0. FPF(KAT)=CATN(KAT) TPF(KAT)=CATS(KAT) DO 1200 I=2,KAT FPF(KAT-I+1)=FPF(KAT-I+2)+CATN(KAT-I+1) TPF(KAT-I+1)=TPF(KAT-I+2)+CATS(KAT-I+1) 1200 CONTINUE DO 1210 I=1,NN FPF(I)=FPF(I)/EMN TPF(I)=TPF(I)/EMS 1210 CONTINUE IF(KATBIG.EQ.0)GO TO 1245 WRITE(6,1220) 1220 FORMAT(/10X,'"OPERATING POINTS" CORRESPONDING TO THE ', + 'INPUT DATA:'/) WRITE(6,1230) (FPF(NN-I+1),I=1,10),(TPF(NN-I+1),I=1,10) 1230 FORMAT(2X,' FPF: ',10(F6.4,1X)/2X,' TPF: ',10(F6.4,1X)/) IF(KAT.GT.10)WRITE(6,1230) (FPF(NN-I+1),I=11,20), + (TPF(NN-I+1),I=11,20) WRITE(6,1240) FPF(1),TPF(1) 1240 FORMAT(2X,' FPF: ',F6.4/2X,' TPF: ',F6.4/) C C CHECK DEGENERACY OF DATA SET. IF IT IS A DEGENERATE DATA SET, C WRITE OUT MESSAGE AND ABORT THE EXECUTION OF THIS DATA SET. C 1245 CALL CHECK(ICLASS,KAT,FPF,TPF,ICON,IARY) IF (ICLASS.EQ.0) GO TO 1250 CALL DEGENE(ICLASS,ICON,IARY,FPF,TPF) IERROR=1 RETURN 1250 DO 1260 I=1,KK X(I)=FPF(I+1) XS(I)=TPF(I+1) 1260 CONTINUE X(kat)=0 Xs(kat)=0 C C TEST FOR 1.'S AND CORRECT BY SUBTRACTING 1/2N. TEST FOR 0'S C AND CORRECT BY CHANGING TO 1/2N. THEN CALL SUBROUTINE ZDEV TO C TRANSFORM THE CUMULATIVE PROBABILITIES IN THE ARRAYS TO STANDARD C NORMAL DEVIATES. C Q=0. DO 1310 I=1,KK IF(X(I).EQ.0.)X(I)=1./(2.*EMN) IF(ABS(X(I)-1.).LT.1.0E-05)X(I)=1.-(1./(2.*EMN)) P=X(I) CALL ZDEV(P,Q,D,IER) X(I)=Q IF(XS(I).EQ.0.)XS(I)=1./(2.*EMS) IF(ABS(XS(I)-1.).LT.1.0E-05)XS(I)=1.-(1./(2.*EMS)) P=XS(I) CALL ZDEV(P,Q,D,IER) XS(I)=Q 1310 CONTINUE C C ADJUST THE SEQUENCE OF CUTOFFS C DO 1320 I=2,KK IF(X(I).GE.X(I-1))X(I)=X(I-1)-0.1 1320 IF(XS(I).GE.XS(I-1))XS(I)=XS(I-1)-0.1 C C CALCULATE LEAST SQUARES SOLUTIONS C SUMX=0. SUMY=0. SUMXY=0. SUMX2=0. XK=FLOAT(KK) DO 1330 I=1,KK SUMX=SUMX+X(I) SUMX2=SUMX2+X(I)**2 SUMY=SUMY+XS(I) SUMXY=SUMXY+XS(I)*X(I) 1330 CONTINUE XMEAN=SUMX/XK YMEAN=SUMY/XK B=(XK*SUMXY-SUMX*SUMY)/(XK*SUMX2-SUMX**2) A=YMEAN-B*XMEAN 1345 DO 1350 I=1,KK X(I)=-X(I) 1350 CONTINUE C C USE METHOD OF SCORING TO GET FINAL ESTIMATE OF A AND B C CALL CUTOFF(A,B,kat,catn,cats,X) CALL MOSLOP(LL,A,B,X) IF(IERROR.EQ.0)GO TO 1370 WRITE(6,1366) 1366 FORMAT(3X,'UNDERFLOW OCCURS IN LABROC5') IERROR=1 RETURN 1370 IF(LL.LT.150)GO TO 1380 WRITE(6,1371) 1371 FORMAT(3X,'DOES NOT CONVERGE -- 100 ITERATIONS') IERROR=1 RETURN 1380 A0=A B0=B IF(KATBIG.EQ.0)RETURN C C PRINT FINAL OUTPUT RESULTS (IF TOTAL RUNS > 400) C CALL OUTRSL(IKEY,A,B,X) RETURN END C--------------------------------------- SUBROUTINE LSTFIT(IKEY,A0,B0) C--------------------------------------- C C USE THE ESTIMATED PARAMETERS OBTAINED FROM "CUTOFF" C AS WELL AS THE A AND B PARAMETERS OBTAINED FROM "LABROC1" AS C INITIAL ESTIMATES IN ORDER TO GET FINAL ESTIMATES FROM C "MOSLOP" C INTEGER IARY(401) DIMENSION X(399),FPF(401),TPF(401) COMMON/E/IEND,IERROR COMMON/B1/MN,MS,KATBIG COMMON/B3/CATN4(400),CATS4(400),KAT4 character*80 ititle COMMON/B5/TRAPA,VTRAPA,ititle,ZZ(399) COMMON/F1/CATN(400),CATS(400),KAT c COMMON/F2/DIFFN(400),DIFFSN(400) C real aeff,beff C common/effi/aeff,beff common/err/iloss KAT=KAT4 NN=KAT+1 EMN=MN EMS=MS DO 1010 I=1,NN CATN(I)=CATN4(I) CATS(I)=CATS4(I) 1010 CONTINUE C C WRITE OUT HEADING, DATA DESCRIPTION, NUMBERS OF C CATEGORIES AND CASES, AND RESPONSE DATA FROM INPUT. C WRITE(6,1030)ititle(1:80) 1030 FORMAT(//,5X,'LABROC4 (Oct 1997 VERSION): ',A//) IF(IKEY.EQ.1)WRITE(6,1040) 1040 FORMAT(10X,'LARGER VALUES OF THE TEST RESULT REPRESENT STRONGER', + ' EVIDENCE THAT',/10x,' THE CASE IS ACTUALLY POSITIVE'/ + 10X,'(E.G., THAT THE PATIENT IS ACTUALLY ABNORMAL)') IF(IKEY.EQ.0)WRITE(6,1041) 1041 FORMAT(10X,'SMALLER VALUES OF THE TEST RESULT REPRESENT STRONGER', + ' EVIDENCE THAT',/10x,' THE CASE IS ACTUALLY POSITIVE'/ + 10X,'(E.G., THAT THE PATIENT IS ACTUALLY ABNORMAL)') WRITE(6,1050)EMN,EMS 1050 FORMAT(10X,'NO. OF ACTUALLY NEGATIVE CASES=',F6.0,/10X, 1 'NO. OF ACTUALLY POSITIVE CASES=',F6.0) cC C CHECK INPUT DATA. C SUMR=0. SUMRP=0. DO 1100 I=1,KAT SUMR=SUMR+CATN(I) SUMRP=SUMRP+CATS(I) 1100 CONTINUE C C WRITE OUT OBSERVED OPERATING POINTS C 1170 FPF(NN)=0. TPF(NN)=0. FPF(KAT)=CATN(KAT) TPF(KAT)=CATS(KAT) DO 1200 I=2,KAT FPF(KAT-I+1)=FPF(KAT-I+2)+CATN(KAT-I+1) TPF(KAT-I+1)=TPF(KAT-I+2)+CATS(KAT-I+1) 1200 CONTINUE DO 1210 I=1,NN FPF(I)=FPF(I)/EMN TPF(I)=TPF(I)/EMS 1210 CONTINUE WRITE(6,1220) 1220 FORMAT(/10X,'"OPERATING POINTS" CORRESPONDING TO THE ', + 'INPUT DATA:'/) INN=INT(NN/10) DO 1240 II=1,INN IB=NN-(II-1)*10 WRITE(6,1230) (FPF(IB-I+(IB-9)),I=IB-9,IB) WRITE(6,1235) (TPF(IB-I+(IB-9)),I=IB-9,IB) 1230 FORMAT(2X,' FPF: ',10(F6.4,1X)) 1235 FORMAT(2X,' TPF: ',10(F6.4,1X)/) 1240 CONTINUE IB=NN-INN*10 IF(IB.EQ.0)GO TO 1245 WRITE(6,1230) (FPF(IB-I+1),I=1,IB) WRITE(6,1235) (TPF(IB-I+1),I=1,IB) C C CHECK DEGENERACY OF DATA SET. IF IT IS A DEGENERATE DATA SET, C WRITE OUT MESSAGE AND ABORT THE EXECUTION OF THIS DATA SET. C 1245 CALL CHECK(ICLASS,KAT,FPF,TPF,ICON,IARY) IF (ICLASS.EQ.0) GO TO 1250 CALL DEGENE(ICLASS,ICON,IARY,FPF,TPF) RETURN 1250 A=A0 B=B0 C C CALCULATE INTIAL ESTIMATES FOR CUTOFFS C CALL CUTOFF(A,B,kat4,catn,cats,X) C C USE METHOD OF SCORING TO FIND FINAL ESTIMATES OF PARAMETERS C CALL MOSLOP(LL,A,B,X) c if (iloss.gt.0) write(*,*) 'L ',iloss IF(LL.LT.150)GO TO 1370 WRITE(6,1365) 1365 FORMAT(5X,'DOES NOT CONVERGE -- 100 ITERATIONS') RETURN 1370 IF(IERROR.EQ.0)GO TO 2190 WRITE(6,1375) 1375 FORMAT(5X,'UNDERFLOW OCCURS IN LABROC4') RETURN 2190 CALL OUTRSL(IKEY,A,B,X) RETURN END C------------------------------------- SUBROUTINE CUTOFF(A,B,kat,catn,cats,TJ) C------------------------------------- C C USE THE FINAL ESTIMATES OF A AND B FROM "LABROC1" C TO ESTIMATE CUTOFFS TJ(K) FOR THE MAIN PART OF THE C CURVE FITTING ALGORITHM. THE BASIC APPROACH IS TO C CALCULATE A MAXIMUM LIKELIHOOD ESTIMATE FOR EACH C TJ(K), GIVEN A AND B. C DIMENSION XJ(399),YJ(399),TJ(399) dimension catn(400),cats(400) COMMON/B3/CATN4(400),CATS4(400),KAT4 COMMON/D1/SUMN(400),SUMS(400),A0,B0 A0=A B0=B SUMN(1)=CATN(1) SUMS(1)=CATS(1) DO 1100 J=2,KAT SUMN(J)=SUMN(J-1)+CATN(J) SUMS(J)=SUMS(J-1)+CATS(J) 1100 CONTINUE KKK=KAT-2 EMN=SUMN(KAT) EMS=SUMS(KAT) C C FIND CUTOFFS TJ(J), J=2,KAT-2 C DO 2100 J=2,KKK P=SUMN(J)/EMN CALL ZDEV(P,X,DXJ,IE) XJ(J)=X P=SUMS(J)/EMS CALL ZDEV(P,Y,DYJ,IE) YJ(J)=(A0+Y)/B0 IF(XJ(J)-YJ(J))1200,1800,1300 1200 RLL=XJ(J) RUL=YJ(J) GO TO 1400 1300 RLL=YJ(J) RUL=XJ(J) C C IF PRODUCT OF H FUNCTION IS NOT LESS THAN 0 (DUE TO C ROUNDING ERROR), SET TJ EQUAL TO THE Average OF RLL AND RUL; C or numerical approximation,OTHERWISE USE C 'INTERVAL HAVING' FIND THE ROOT OF H FUNCTION C 1400 HFL=HF(J,RLL) HFU=HF(J,RUL) PHF=HFL*HFU IF(PHF)1600,1500,1500 1500 RLL=RLL-0.001 RUL=RUL+0.001 GO TO 1400 1600 CALL SETTJ(RLL,RUL,J,TJ(J)) GO TO 2100 1800 TJ(J)=XJ(J) 2100 CONTINUE C C FIND THE FIRST CUTOFF C RUL=TJ(2) TRIAL=TJ(2) 3000 TRIAL=TRIAL-0.1 if (RUL.eq.trial) goto 3200 IF(HF(1,TRIAL)*HF(1,RUL))3100,3200,3000 3100 RLL=TRIAL CALL SETTJ(RLL,RUL,1,TJ(1)) GO TO 4100 3200 TJ(1)=TRIAL C C FIND THE LAST CUTOFF C 4100 KK=KAT-1 RLL=TJ(KKK) TRIAL=TJ(KKK) 4200 TRIAL=TRIAL+0.1 if (RLL.eq.trial) goto 4400 jll= HF(KK,RLL) tll=HF(KK,TRIAL) if (tll.gt.1.0e+8) then tll = -1.0 endif IF(jll*tll)4300,4400,4200 4300 RUL=TRIAL CALL SETTJ(RLL,RUL,KK,TJ(KK)) RETURN 4400 TJ(KK)=TRIAL RETURN END C------------------------------------------- SUBROUTINE SETTJ(RLL,RUL,J,T) C------------------------------------------- C C USE THE 'INTERVAL HALVING' METHOD TO FIND ROOT OF C EQUATION HF=0 WITH GIVING LOWER AND UPPER BOUND (RLL,RUL) C C COMMON/D1/SUMN(400),SUMS(400),A0,B0 C COMMON/B1/MN,MS,KATBIG 1000 HR=(RLL+RUL)/2. IF(HF(J,RLL)*HF(J,HR))1100,1200,1300 1100 RUL=HR GO TO 1400 1200 T=HR RETURN 1300 RLL=HR 1400 VABSR=ABS(RUL-RLL) IF(VABSR.GE.0.001)GO TO 1000 T=(RUL+RLL)/2. RETURN END C------------------------------- FUNCTION HF(J,TJ) C------------------------------- C C HF FUNCTION ---- changed 1/16/91 C COMMON/D1/SUMN(400),SUMS(400),A0,B0 COMMON/B1/MN,MS,KATBIG DOUBLE Precision UF,VF,STP,BTA EMN=MN EMS=MS STP=2.506628275D0 BTA = B0*TJ-A0 RA=BTA CALL NDTR(TJ,PTJ,DJ) CALL NDTR(RA,QTJ,DJJ) IF (TJ.lt.-5) THEN UF = -(STP * EXP((TJ**2)*0.5)*SUMN(J)*TJ**3/(TJ**2-1.)+ + EMN - SUMN(J)) ELSE IF (TJ.GT.5) THEN UF = SUMN(J)- STP* EXP((TJ**2)*0.5)*TJ**3/(TJ**2-1.)* + (EMN - SUMN(J)) ELSE UF = SUMN(J)/PTJ - (EMN - SUMN(J))/(1. - PTJ) ENDIF IF (BTA.LT.-5) THEN VF = -B0*EXP((TJ**2)*0.5)*(STP*SUMS(J)*BTA**3./(BTA**2-1.)+ + (EMS - SUMS(J))*EXP(BTA**2*(-0.5))) ELSE IF (BTA.gt.5) THEN VF = B0*EXP((TJ**2)*0.5)*(SUMS(J)*EXP(BTA**2* (-0.5))- + STP*(EMS-SUMS(J)*BTA**3/(BTA**2-1.))) ELSE VF = B0*EXP(0.5*((1.-B0**2)*TJ**2+2.*A0*B0*TJ-A0**2))* + (SUMS(J)/QTJ-(EMS-SUMS(J))/(1.-QTJ)) ENDIF HF = UF + VF 1500 FORMAT('*** ATTEMPTING TO CORRECT FOR DIVIDE BY ZERO ERROR ***') RETURN END C------------------------------- C FUNCTION HF(J,TJ) C------------------------------- C C HF FUNCTION CC C COMMON/D1/SUMN(400),SUMS(400),A0,B0 C COMMON/B1/MN,MS,KATBIG C EMN=MN C EMS=MS C CALL NDTR(TJ,PTJ,DJ) C CALL NDTR(B0*TJ-A0,QTJ,DJJ) C HF=(SUMN(J)/PTJ-(EMN-SUMN(J))/(1.-PTJ))+ C + B0*EXP(0.5*((1.-B0**2)*TJ**2+2.*A0*B0*TJ-A0**2))* C + (SUMS(J)/QTJ-(EMS-SUMS(J))/(1.-QTJ)) C RETURN C END C------------------------------------------- SUBROUTINE OUTRSL(IKEY,A,B,X) C------------------------------------------- C C PRINT FINAL RESULTS C DIMENSION CX(401,401),FPVAL(26),TPVAL(26),BDLOW(26),BDUPP(26), + EFPF(8000),ETPF(8000),X(*) C real aeff,beff C common/effi/aeff,beff COMMON/XP/XP(401,401) COMMON/B1/MN,MS,KATBIG character*80 ititle COMMON/B5/TRAPA,VTRAPA,ititle,ZZ(399) COMMON/F1/CATN(400),CATS(400),KAT COMMON/D1/SUMN(400),SUMS(400),A0,B0 c COMMON/F2/DIFFN(400),DIFFSN(400) COMMON/DECIML/MAXDEC COMMON/B2/ACTN(4000),ACTSN(4000),V(2,8000) common/err/iloss character*3 plmi integer total,catsum,fr2 integer vt,vcount(8000) common/vc/vcount,vt real U(8000),ZU(8000) real vdiff,vsum,vsum2 DATA FPVAL/0.005,0.01,0.02,0.03,0.04,0.05,0.06, 2 0.07,0.08,0.09,0.10,0.11,0.12, 3 0.13,0.14,0.15,0.20,0.25,0.30, 4 0.40,0.50,0.60,0.70,0.80,0.90,0.95/ NN=KAT+1 KK=KAT-1 C C GET CORRELATION MATRIX ON FINAL ITERATION C 2190 DO 2195 I=1,NN DO 2195 J=1,NN CX(I,J)=XP(I,J)/SQRT(XP(I,I)*XP(J,J)) 2195 CONTINUE C C GET VALUE OF LOG L C CALL VLOGL(FLOGL) C C PRINT FINAL ESTIMATES C WRITE(6,2360)A,B 2360 FORMAT(//10X,'FINAL ESTIMATES OF THE BINORMAL ROC PARAMETERS ', 1 'A AND B:'/27X,'A =',F7.4,4X,'B =',F7.4) WRITE(6,2410) 2410 FORMAT(//10X,'ESTIMATED STD. ERROR AND CORRELATION OF ', 1 'THESE PARAMETER VALUES:') STDA=SQRT(XP(1,1)) STDB=SQRT(XP(2,2)) WRITE(6,2420)STDA,STDB,CX(1,2) 2420 FORMAT(17X,'STD. ERROR (A) =',F7.4/17X,'STD. ERROR (B) =',F7.4 1 /17X,'CORR(A,B) = ',F7.4) C C SUMMARY OF ROC CURVE C WORKV=1+B*B ADEV=A/SQRT(WORKV) CALL NDTR(ADEV,AREA,DENS) CA2=XP(1,1) CB2=XP(2,2) CAB=XP(1,2) VAREA=SQRT((DENS*DENS)*(CA2/WORKV+((A*B)**2)*CB2/(WORKV**3)- 1 2*A*B*CAB/(WORKV*WORKV))) WRITE(6,2490)AREA,VAREA,TRAPA,VTRAPA 2490 FORMAT(//10X,'AREA UNDER FITTED CURVE (A-SUB-Z) = ',F7.4, + /11x, 'ESTIMATED STD. ERROR OF A-SUB-Z = ',F7.4, + /10X,'TRAPEZOIDAL(WILCOXON) AREA = ',F7.4,5X, + /11x,'ESTIMATED STD. ERROR OF TRAPEZOIDAL AREA = ',F7.4) C if (katbig.ne.0) kat=-kat C WRITE(6,2495)ititle,kat C 2495 FORMAT(7X,A80,I4) C write(*,2491) Trapa,vtrapa C write(*,*) aeff,beff 2491 Format(F7.4,f7.4) c write(*,*) '#, Af, Bf, sdAz, sdZa, Mn, Ms' c write(*,2492) A,B,Area,Mn,Ms c write(*,2492) stdA,stdb,varea,Mn,Ms 2492 format(1x,f12.8,',',f12.8,',',f12.8, + ',',i4,',',i4) C C WITH GIVEN FALSE-POSITIVE FRACTIONS, COMPUTE THE TRUE- C POSITIVE FRACTIONS ON THE ESTIMATED ROC CURVE, WITH C LOWER AND UPPER BOUNDS OF THE ASYMMETRIC 95% CONFIDENCE C INTERVAL. C DO 3205 I=1,26 P=FPVAL(I) CALL ZDEV(P,Q,D,IER) TPDEV=A+B*Q SI=SQRT(Q*Q*CB2+CA2+2.*Q*CAB) BOUND1=TPDEV-1.96*SI BOUND2=TPDEV+1.96*SI CALL NDTR(TPDEV,P,D) TPVAL(I)=P CALL NDTR(BOUND1,P,D) BDLOW(I)=P CALL NDTR(BOUND2,P,D) BDUPP(I)=P 3205 CONTINUE WRITE(6,3210) (FPVAL(I),TPVAL(I),BDLOW(I),BDUPP(I),I=1,26) 3210 FORMAT(///10X,'ESTIMATED BINORMAL ROC CURVE, WITH LOWER AND' 1 /10X,'UPPER BOUNDS OF THE ASYMMETRIC 95% CONFIDENCE' 2 /10X,'INTERVAL FOR TRUE-POSITIVE FRACTION AT A VARIETY' 3 /10X,'OF FALSE-POSITIVE FRACTIONS:', 4 //5X,'*****************************************************', 5 '*****', 6 ///13X,'FPF',6X,'TPF',8X,'(LOWER BOUND, UPPER BOUND)'/ 7 26(/12X,F5.3,4X,F6.4,6X,'(',3X,F6.4,2X,',',3X,F6.4,3X,')')) C C COMPUTE THE TPF'S AND FPF'S FOR EACH EXPECTED OPERATING POINT C ESTIMATED ON THE FITTED CURVE. C total = 0 j = 1 i=1 do while (j.le.kat.and.i.le.vt) total = total + vcount(i) catsum = sumn(j) + sums(j) if (total.lt.catsum.and. + (cats(j).eq.0.or.catn(j).eq.0)) then fr2=catsum-total if (catn(j) .gt. 0) then U(i) = ccv(kat,j,fr2,0.0,1.0,catn(j),X) zu(i) = czv(i,vt,0.0,1.0,U) else U(i) = ccv(kat,j,fr2,a,b,cats(j),X) zu(i) = czv(i,vt,a,b,U) endif else U(i) = x(j) if (cats(j).ne.0.and.catn(j).ne.0) then zu(i)= (catn(j)*czv(i,vt,0.0,1.0,U) + + cats(j)*czv(i,vt,a,b,U))/(catn(j)+cats(j)) else if (catn(j).ne.0) then zu(i)= czv(i,vt,0.0,1.0,U) else zu(i)= czv(i,vt,a,b,U) endif j=j+1 endif Q=-zu(I) CALL NDTR(Q,P,D) EFPF(I)=P Q=a-B*zu(I) CALL NDTR(Q,P,D) ETPF(I)=P i=i+1 enddo c vt=vt-1 c DO 3305 I=1,vt c 3305 CONTINUE WRITE(6,3306) 3306 FORMAT(///,4X,'Estimated Relationship', 1 ' Between the Critical Test-Result Value' 1 /4X,'(which separates', 2 1H','POSITIVE',1H',' results from ',1H','NEGATIVE',1H', 3 ' results)',/,15X,'and its Corresponding Operating Point', 4 /10X,' Projected onto the Fitted Binormal ROC Curve:'/ 5/5X,'******************************************************', 6 '*****************') 3307 FORMAT(//,a5,7X,'CRITICAL TEST',5x,'Normal', 1 9X,'( FPF , TPF )',/12X,'RESULT VALUE',6x,'Deviate') c write(*,3308) 3308 format( 10x,'(or Run Boundary)',//) IF(katbig.eq.0) GO TO 3500 Write(6,3307) ' ' write(*,*) write(*,*) plmi=' ' DO 3320 I=1,kat-1 II=kat-I IF(MAXDEC.EQ.0) + WRITE(6,3310)plmi, zz(II), zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.1) + WRITE(6,3311)plmi, zz(II), zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.2) + WRITE(6,3312)plmi, zz(II), zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.3) + WRITE(6,3313)plmi, zz(II), zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.4) + WRITE(6,3314)plmi, zz(II), zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.5) + WRITE(6,3315)plmi, zz(II), zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.6) + WRITE(6,3316)plmi, zz(II), zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.7) + WRITE(6,3317)plmi, zz(II), zu(II),EFPF(II),ETPF(II) IF(MAXDEC.GE.8) + WRITE(6,3318)plmi, zz(II), zu(II),EFPF(II),ETPF(II) 3320 CONTINUE RETURN 3310 FORMAT(2x,a3,5X,F8.0,9x,F9.4,1X, + 8x,'(',F6.3,',',1X,F6.3,')') 3311 FORMAT(2x,a3,4X,F9.1,9x,F9.4,1X, + 8x,'(',F6.3,',',1X,F6.3,')') 3312 FORMAT(2x,a3,3X,F10.2,9x,F9.4,1X, + 8x,'(',F6.3,',',1X,F6.3,')') 3313 FORMAT(2x,a3,2X,F11.3,9x,F9.4,1X, + 8x,'(',F6.3,',',1X,F6.3,')') 3314 FORMAT(2x,a3,1X,F12.4,9x,F9.4,1X, + 8x,'(',F6.3,',',1X,F6.3,')') 3315 FORMAT(2x,a3,1X,F12.5,9x,F9.4,1X, + 8x,'(',F6.3,',',1X,F6.3,')') 3316 FORMAT(2x,a3,1X,F12.6,9x,F9.4,1X, + 8x,'(',F6.3,',',1X,F6.3,')') 3317 FORMAT(2x,a3,1X,F12.7,9x,F9.4,1X, + 8x,'(',F6.3,',',1X,F6.3,')') 3318 FORMAT(2x,a3,F13.8,9x,F9.4,1X, + 8x,'(',F6.3,',',1X,F6.3,')') 3500 i=i IB=vt j=kat X(kat)=99999.99 ivt=total+1 vsum=0 vsum2=0 write(*,3307) 'TRUTH' write(*,3308) DO 3520 I=1,vt II=vt-I+1 C REMEMBER TO ADD IN THESES LINES LATER OR SOMETHING SIMILAR C TO LIMIT THEMNUMBER OF POINTS PRONTED OUT c DELTA=SQRT((EFPF(II)-EFPF(IB))**2+(ETPF(II)-ETPF(IB))**2) c IF(DELTA.LT.0.02)GO TO 3520 ivt=ivt-vcount(II) if (v(2,ivt).ne.v(2,ivt+1).or.vcount(ii).gt.1) then insum=0 do inx=1,vcount(ii) insum=insum+v(2,ivt+inx-1) enddo if (j.gt.0.and. + ( vcount(ii)-abs(insum).gt.0.or. + v(2,ivt).ne.v(2,ivt+vcount(ii)) )) then if (j.lt.kat) + write(*,'(12x,a,f9.4)') 'run boundary = ', x(j) j=j-1 if (vcount(ii)-abs(insum).gt.0) v(2,ivt)=0 endif endif vtmp = v(2,ivt) if (vcount(ii).gt.1) then v(2,ivt)=vcount(ii) endif c vdiff=v(1,ivt)-zu(II) c vsum=vsum+vdiff c vsum2=vsum2+vdiff**2 if (ikey.eq.0) v(1,ivt)=-v(1,ivt) plmi='tie' if (v(2,ivt).eq.1) plmi='+' if (v(2,ivt).eq.-1) plmi='-' IF(MAXDEC.EQ.0)WRITE(6,3310) + plmi,v(1,ivt), zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.1)WRITE(6,3311) + plmi,v(1,ivt), zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.2)WRITE(6,3312) + plmi,v(1,ivt),zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.3)WRITE(6,3313) + plmi,v(1,ivt),zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.4)WRITE(6,3314) + plmi,v(1,ivt),zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.5)WRITE(6,3315) + plmi,v(1,ivt),zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.6)WRITE(6,3316) + plmi,v(1,ivt),zu(II),EFPF(II),ETPF(II) IF(MAXDEC.EQ.7)WRITE(6,3317) + plmi,v(1,ivt),zu(II),EFPF(II),ETPF(II) IF(MAXDEC.GE.8)WRITE(6,3318) + plmi,v(1,ivt),zu(II),EFPF(II),ETPF(II) IB=II v(2,ivt)=vtmp 3520 CONTINUE c vsum=vsum/vt c vsum2=sqrt((vsum2/(vt-1))-(vt*(vsum**2)/(vt-1))) c write(*,*) ' Diff = ', vsum, 'Sd = ', vsum2 c write(*,*) 't-stat = ',vt-1, c + vsum/(vsum2/sqrt(float(vt))) RETURN END C--------------------------------------- SUBROUTINE MOSLOP(LL,A,B,X) C--------------------------------------- C C USE METHOD OF SCORING TO GET THE BEST ESTIMATS C DIMENSION XXDUM(80700),ADJX(399),X(*) INTEGER SWT COMMON/XP/XP(401,401) COMMON/E/IEND,IERROR COMMON/B1/MN,MS,KATBIG COMMON/F1/CATN(400),CATS(400),KAT COMMON/F2/DIFFN(400),DIFFSN(400) COMMON/G/AAA,BBB,PXBA(401),DXBA(401),PX(401),DX(401),XX(401,401), + PARX(399) common/err/iloss real oadja,oadjb,ecor,esttest adja=0 adjb=0 alpha=1.0 NNN=401 NN=KAT+1 KK=KAT-1 EMN=MN EMS=MS C C INITIALIZE VARIABLES C DO 7010 J=1,NNN DO 7010 K=1,NNN XX(J,K)=0. XP(J,K)=0. 7010 CONTINUE SWT=0 LL=0 C C GET INTEGRALS AND DENSITIES FOR THE INITIAL ESTIMATES C PXBA(1)=0. PXBA(NN)=1. DXBA(1)=0. DXBA(NN)=0. PX(1)=0. PX(NN)=1. DX(1)=0. DX(NN)=0. DO 7020 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)) 7020 CONTINUE DO 7030 I=1,KAT DIFFSN(I)=PXBA(I+1)-PXBA(I) DIFFN(I)=PX(I+1)-PX(I) IF (DIFFSN(I).LT.5.0E-8) DIFFSN(I)=5.0E-8 IF (DIFFN(I).LT.5.0E-8) DIFFN(I)=5.0E-8 7030 CONTINUE CALL VLOGL(FUNLIK) C C ITERATION LOOPS START HERE C iloss=0 7050 LL=LL+1 CALL PDEV(A,B,X) oadja = adja oadjb = adjb C C INVERT PARTIAL DERIVATIVE MATRIX C DO 7135 I=1,NN DO 7135 J=1,NN 7135 XX(I,J)=-XX(I,J) CALL M2TOM1(NN,XX,XXDUM) CALL SINV(XXDUM,NN,.001,IER) C C IF PROBLEMS ARE ENCOUNTERED IN MATRIX INVERSION, REDUCE C ADJUSTMENTS TO SOLUTION VECTOR AND GO ON TO TEST CHANGE C IN LOG-LIKELIHOOD FUNCTION C IF(IER.GE.0) GO TO 7137 ALPHA=-0.5*ABS(ALPHA) IF(ABS(ALPHA).GT.1.0E-8)GO TO 7157 IERROR=1 RETURN cc 7137 if (ier.gt.0) iloss= iloss+1 7137 IF(IER.GT.0)WRITE(6,7140)IER 7140 FORMAT(' LOSS OF SIGNIFICANCE. STEP ',I5,'+1') CALL M1TOM2(NN,XXDUM,XP) IF(SWT.EQ.1)GO TO 7190 IF(KAT.LE.3)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 Sometimes Newtons method can get stuck in a loop C where it bounces back and forth b/w two values C Assume the true value is in b/w the extremes and C re-adjust accordingly. C (Only allow this to happen once per iteration. C The effect is like damping the damping) esttest=abs((oadja+adja+oadjb+adjb)/(adja+adjb)) if ((esttest.lt.0.15.and.abs(alpha).ge.0.5)) then ecor=ecor+1 else ecor=0 endif if (iter.gt.60) ecor=1 C C ITERATE C ALPHA=1.0 7157 A=A+ALPHA*ADJA*0.5**ecor B=B+ALPHA*ADJB*0.5**ecor DO 7160 I=1,KK 7160 X(I)=X(I)+ALPHA*ADJX(I)*0.5**ecor ecor=0 C C COMPUTE THE TERMS USED IN CALCULATION OF LOG-LIKELIHOOD FUNCTION C 7170 DO 7172 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)) 7172 CONTINUE DO 7173 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 7173 CONTINUE IF (ALPHA.NE.1.0 .OR. LL.EQ.1) GO TO 7175 C C CHECK FOR MAXIMIZATION OF LOG-LIKELIHOOD C TNET=ABS(ADJA)+ABS(ADJB) DO 7174 I=1,KK TNET=TNET+ABS(ADJX(I)) 7174 CONTINUE TNET=6.*TNET/FLOAT(1+KAT) IF (TNET.GT.0.001) GO TO 7175 SWT=1 GO TO 7050 C C CALCULATE LOG-LIKELIHOOD FUNCTION VALUE. C IF LOG-LIKELIHOOD FUNCTION VALUE FOR LL-TH STEP DECREASES, C HALVE ALPHA VALUE UNTIL ITS VALUE GOES UP AGAIN. C OTHERWISE, JUMP TO LL+1-TH STEP. C 7175 IF(ALPHA.EQ.1.)FUNLI1=FUNLIK CALL VLOGL(FUNLIK) IF (FUNLIK.GE.FUNLI1) GO TO 7185 FUNDIF=ABS((FUNLIK-FUNLI1)/tnet) IF(FUNDIF.LE.1.0E-4)GO TO 7185 ALPHA=-0.5*ABS(ALPHA) IF(ABS(ALPHA).GT.1.0E-8)GO TO 7157 C IERROR=1 C RETURN C C ALLOW UP TO 100 ITERATIONS C 7185 IF(LL.LE.150)GO TO 7050 7190 RETURN END C----------------------------------- SUBROUTINE PDEV(A,B,X) C----------------------------------- C C GET THE FIRST AND SECOND PARTIAL DERIVATIVES C COMMON/B1/MN,MS,KATBIG COMMON/F1/CATN(400),CATS(400),KAT COMMON/F2/DIFFN(400),DIFFSN(400) COMMON/G/AAA,BBB,PXBA(401),DXBA(401),PX(401),DX(401),XX(401,401), + PARX(399) DIMENSION X(*) EMN=MN EMS=MS C KK=KAT-1 C NN=KAT+1 C C GET FIRST PARTIAL DERIVATIVES 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)/DIFFSN(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 PARTIAL DERIVATIVES 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 RETURN END C----------------------------------- SUBROUTINE VLOGL(FLOGL) C----------------------------------- C C GET VALUE OF LOG L C COMMON/F1/CATN(400),CATS(400),KAT COMMON/F2/DIFFN(400),DIFFSN(400) NN=KAT+1 SUMN=0. SUMM=0. DO 7330 I=2,NN SUMM=SUMM+CATS(I-1)*ALOG(DIFFSN(I-1)) SUMN=SUMN+CATN(I-1)*ALOG(DIFFN(I-1)) 7330 CONTINUE FLOGL=SUMN+SUMM RETURN END C--------------------------------------------- SUBROUTINE M2TOM1(N,XX,VXX) C--------------------------------------------- C C THIS SUBROUTINE STORES THE UPPER TRIANGULAR PART OF C A SYMMETRIC (N BY N) MATRIX XX COLUMNWISE AS A VECTOR C VXX WITH LENGTH N*(N+1)/2. C DIMENSION XX(401,401),VXX(*) K=0 DO 100 J=1,N DO 100 I=1,J K=K+1 VXX(K)=XX(I,J) 100 CONTINUE RETURN END C--------------------------------------- SUBROUTINE M1TOM2(N,VXX,XX) C--------------------------------------- C C THIS SUBROUTINE CREATES A SYMETRIC (N BY N) MATRIX FROM C THE VECTOR VXX WITH LENGTH N*(N+1)/2, WHICH REPRESENTS C THE UPPER TRIANGULAR PART OF XX BY COLUMNS. C DIMENSION XX(401,401),VXX(*) K=0 DO 100 J=1,N DO 100 I=1,J K=K+1 XX(I,J)=VXX(K) IF(I.NE.J)XX(J,I)=XX(I,J) 100 CONTINUE RETURN END C------------------------------------------------------------------------------- SUBROUTINE CHECK(ICLASS,KAT,FPF,TPF,ICON,IARY) C------------------------------------------------------------------------------- C C PURPOSE C CHECK DEGENERACY OF DATA SET. C C DESCRIPTION OF PARAMETERS C ICLASS - OUTPUT CLASS NUMBER OF DEGENERACY. (USED BY C SUBROUTINE 'DEGENE' IF NECESSARY) C KAT - NUMBER OF CATEGORIES OF DATA SET. C FPF,TPF - OBSERVED OPERATING POINTS. C ICON,IARY - STORAGE SPACE FOR INFORMATION IN MESSAGE. (USED BY C SUBROUTINE 'DEGENE' IF NECESSARY) C REAL FPF(401),TPF(401) INTEGER IARY(401) ICLASS=0 do i=1,14 iary(i) = 0 enddo NN=KAT+1 ICON=0 DO 10 I=2,KAT T=FPF(I)*TPF(I)*(1.-FPF(I))*(1.-TPF(I)) IF (T.EQ.0.) GO TO 10 ICON=ICON+1 IARY(ICON)=I 10 CONTINUE IF (KAT.EQ.1) GO TO 20 IF (ICON.EQ.0) GO TO 40 I0=IARY(1)-1 I1=IARY(ICON)+1 IF (KAT.EQ.2.AND.ICON.NE.0) GO TO 30 IF (ICON.GE.2) GO TO 15 IF ((TPF(I0).NE.1.0.OR.TPF(I1).NE.0.0).AND. 1 (FPF(I0).EQ.1.0.AND.FPF(I1).EQ.0.0)) ICLASS=5 IF ((FPF(I0).NE.1.0.OR.FPF(I1).NE.0.0).AND. 1 (TPF(I0).EQ.1.0.AND.TPF(I1).EQ.0.0)) ICLASS=6 RETURN 15 IF ((TPF(IARY(1)).EQ.TPF(IARY(ICON))).AND. 1 (FPF(I0).EQ.1.AND.FPF(I1).EQ.0)) ICLASS=7 IF ((TPF(I0).EQ.1.AND.TPF(I1).EQ.0).AND. 1 (FPF(IARY(1)).EQ.FPF(IARY(ICON)))) ICLASS=8 RETURN 20 ICLASS=1 RETURN 30 ICLASS=2 RETURN 40 ICA=0 ICB=0 DO 50 I=1,NN IF (FPF(I).NE.0.AND.FPF(I).NE.1) ICA=ICA+1 IF (TPF(I).NE.0.AND.TPF(I).NE.1) ICB=ICB+1 50 CONTINUE I1=0 I2=0 DO 60 I=2,KAT IF (FPF(I).EQ.1.AND.TPF(I).NE.1) GO TO 55 IF (FPF(I).EQ.0.AND.TPF(I).NE.0) I2=I2+1 IF (I2.EQ.1) IARY(2)=I GO TO 60 55 I1=I1+1 IARY(1)=I 60 CONTINUE II1=0 II2=0 DO 70 I=2,KAT IF (FPF(I).NE.1.AND.TPF(I).EQ.1) GO TO 65 IF (FPF(I).NE.0.AND.TPF(I).EQ.0) II2=II2+1 IF (II2.EQ.1) IARY(4)=I GO TO 70 65 II1=II1+1 IARY(3)=I 70 CONTINUE IF ((ICA.NE.0.AND.ICB.NE.0).OR.((ICA.EQ.0.AND. 1(I1.EQ.0.OR.I2.EQ.0)).OR.(ICB.EQ.0.AND.(II1.EQ.0.OR.II2.EQ.0)))) 2 GO TO 80 IF (ICA.EQ.0.AND.I1.GE.1.AND.I2.GE.1) GO TO 75 IF (ICB.EQ.0.AND.II1.GE.1.AND.II2.GE.1) ICLASS=4 IARY(1)=IARY(3) IARY(2)=IARY(4) RETURN 75 ICLASS=3 RETURN 80 KK1=IARY(3)+1 K1=IARY(1)+1 KK2=IARY(4)-1 K2=IARY(2)-1 IF ((II1.NE.0.AND.FPF(KK1).EQ.0).OR.(I2.NE.0.AND. 1 TPF(K2).EQ.1)) ICLASS=9 IF ((I1.NE.0.AND.TPF(K1).EQ.0).OR.(II2.NE.0.AND. 1 FPF(KK2).EQ.1)) ICLASS=10 RETURN END C------------------------------------------------------------------------------- SUBROUTINE DEGENE(ICLASS,ICON,IARY,FPF,TPF) C------------------------------------------------------------------------------- C C PURPOSE C OUTPUT MESSAGE FOR EACH DEGENERACY CLASS. C C DESCRIPTION OF PARAMETERS C ICLASS - INPUT CLASS NUMBER OF DEGENERACY FROM SUBROUTINE C 'CHECK'. C ICON,IARY - INPUT INFORMATION FROM SUBROUTINE 'CHECK'. C FPF,TPF - OBSERVED OPERATING POINTS. C REAL FPF(*),TPF(*) INTEGER IARY(*) IF (ICLASS.EQ.1) GO TO 10 IF (ICLASS.EQ.2) GO TO 20 IF (ICLASS.EQ.3) GO TO 30 IF (ICLASS.EQ.4) GO TO 40 IF (ICLASS.EQ.5.OR.ICLASS.EQ.7) GO TO 50 IF (ICLASS.EQ.6.OR.ICLASS.EQ.8) GO TO 60 IF (ICLASS.EQ.9) GO TO 90 WRITE(6,5) 5 FORMAT(///,12X, 1 'DATA ARE DEGENERATE AND IMPLY PERFECT (BUT PERVERSE)' 2 ,/,12X,'DECISION PERFORMANCE.') RETURN 10 WRITE(6,15) 15 FORMAT(///,12X, 1 'DATA ARE DEGENERATE AND PROVIDE NO OPERATING POINTS',/ 2 ,12X,'OFF THE (0,0) AND (1,1) CORNERS.') RETURN 20 P=FPF(IARY(1)) CALL ZDEV(P,QF,D,IER) P=TPF(IARY(1)) CALL ZDEV(P,QT,D,IER) DD=QT-QF WRITE(6,25) DD 25 FORMAT(///,12X, 1 'DATA ARE DEGENERATE AND PROVIDE ONLY A SINGLE ',/, 2 12X,'OPERATING POINT OFF THE (0,0) AND (1,1) CORNERS.', 3 /,12X,'BINORMAL ROC CURVE CANNOT BE ESTIMATED UNIQUELY. ', 4 /,12X,'UNIT SLOPE ROC CURVE (B=1) WOULD HAVE A=D',2H'=, 5 F8.4,'.') RETURN 30 WRITE(6,35) TPF(IARY(2)),TPF(IARY(1)) 35 FORMAT(///,12X, 1 'DATA ARE DEGENERATE AND PRODUCE NO OPERATING POINTS', 2 /,12X,'OFF THE BORDERS OF THE UNIT SQUARE. IMPLIED', 3 /,12X,'EXACT-FIT BINORMAL ROC CURVE IS HORIZONTAL AT', 4 /,12X,'UNDETERMINED HEIGHT BETWEEN TPF=',F8.4,' AND', 5 /,12X,'TPF=',F8.4,' .') RETURN 40 WRITE(6,45) FPF(IARY(2)),FPF(IARY(1)) 45 FORMAT(///,12X, 1 'DATA ARE DEGENERATE AND PRODUCE NO OPERATING POINTS',/ 2 ,12X,'OFF THE BORDERS OF THE UNIT SQUARE. IMPLIED ', 3 /,12X,'EXACT-FIT BINORMAL ROC CURVE IS VERTICAL AT ', 4 /,12X,'UNDETERMINED POSITION BETWEEN FPF=',F8.4, 5 /,12X,'AND FPF=',F8.4,' .') RETURN 50 WRITE(6,55) TPF(IARY(1)) 55 FORMAT(//,12X,'DATA ARE DEGENERATE. IMPLIED EXACT-FIT', 1 /,12X,'BINORMAL ROC CURVE IS HORIZONTAL AT CONSTANT', 2 /,12X,'TPF=',F8.4,' .') RETURN 60 WRITE(6,65) FPF(IARY(1)) 65 FORMAT(//,12X,'DATA ARE DEGENERATE. IMPLIED EXACT-FIT', 1 /,12X,'BINORMAL ROC CURVE IS VERTICAL AT CONSTANT', 2 /,12X,'FPF=',F8.4,' .') RETURN 90 WRITE(6,95) 95 FORMAT(///,12X, 1 'DATA ARE DEGENERATE AND IMPLY PERFECT DECISION', 2 /,12X,'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 SINV(A,N,EPS,IER) C--------------------------------------------------- C C PURPOSE C INVERT A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX C C USAGE C CALL SINV(A,N,EPS,IER) C C DESCRIPTION OF PARAMETERS C A - UPPER TRIANGULAR PART OF THE GIVEN SYMMETRIC POSITIVE C DEFINITE N BY N COEFFICIENT MATRIX. ON RETURN A C CONTAINS THE RESULTANT UPPER TRIANGULAR MATRIX. C N - THE NUMBER OF ROW (COLUMNS) IN GIVEN MATRIX. C EPS - AN INPUT CONSTANT WHICH IS USED AS RELATIVE TOLERANCE C FOR TEST ON LOSS OF SIGNIFICANCE. C IER - RESULTING ERROR PARAMETER CODED AS FOLLOWS: C IER=0 - NO ERROR C IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAMETER N OR C BECAUSE SOME RADICAND IS NONPOSITIVE (MATRIX A C IS NOT POSITIVE DEFINITE, POSSIBLY DUE TO LOSS C OF SIGNIFICANCE) C IER=K - WARNING WHICH INDICATES LOSS OF SIGNIFICANCE. C THE RADICAND FORMED AT FACTORIZATION STEP K+1 C WAS STILL POSITIVE BUT NO LONGER GREATER THAN C ABS(EPS*A(K+1,K+1)). C C REMARKS C THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE C STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS. IN C THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGULAR MATRIX C IS STORED COLUMNWISE TOO. C THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL C CALCULATED RADICANDS ARE POSITIVE. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED. C MFSD C C METHOD C SOLUTION IS DONE USING THE FACTORIZATION BY SUBROUTINE MFSD. C DIMENSION A(*) DOUBLE PRECISION DIN,WORK C C FACTORIZE GIVEN MATRIX BY MEANS OF SUBROUTINE MFSD C C A=TRANSPOSE(T)*T C CALL MFSD(A,N,EPS,IER) IF(IER)9,1,1 C C INVERT UPPER TRIANGULAR MATRIX T C PREPARE INVERSION-LOOP C 1 IPIV=N*(N+1)/2 IND=IPIV C C INITIALIZE INVERSION-LOOP C DO 6 I=1,N DIN=1.D0/DBLE(A(IPIV)) A(IPIV)=DIN MIN=N KEND=I-1 LANF=N-KEND IF(KEND)5,5,2 2 J=IND C C INITIALIZE ROW-LOOP C DO 4 K=1,KEND WORK=0.D0 MIN=MIN-1 LHOR=IPIV LVER=J C C START INNER LOOP C DO 3 L=LANF,MIN LVER=LVER+1 LHOR=LHOR+L 3 WORK=WORK+DBLE(A(LVER)*A(LHOR)) C C END OF INNER LOOP C A(J)=-WORK*DIN 4 J=J-MIN C C END OF ROW-LOOP C 5 IPIV=IPIV-MIN 6 IND=IND-1 C C END OF INVERSION LOOP C C CALCULATE INVERSE(A) BY MEANS OF INVERSE(T) C INVERSE(A)=INVERSE(T)*TRANSPOSE(INVERSE(T)) C INITIALIZE MULTIPLICATION LOOP C DO 8 I=1,N IPIV=IPIV+I J=IPIV C C INITIALIZE ROW-LOOP C DO 8 K=I,N WORK=0.D0 LHOR=J C C START INNER LOOP C DO 7 L=K,N LVER=LHOR+K-I WORK=WORK+DBLE(A(LHOR)*A(LVER)) 7 LHOR=LHOR+L C C END OF INNER LOOP C A(J)=WORK 8 J=J+K C C END OF ROW- AND MULTIPLICATION-LOOP C 9 RETURN END C------------------------------------------- SUBROUTINE MFSD(A,N,EPS,IER) C------------------------------------------- C C PURPOSE C FACTOR A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX C C DESCRIPTION OF PARAMETERS C A - UPPER TRIANGULAR PART OF THE GIVEN SYMMETRIC POSITIVE C DEFINITE N BY N COEFFICIENT MATRIX. ON RETURN A C CONTAINS THE RESULTANT UPPER TRIANGULAR MATRIX. C N - THE NUMBER OF ROW (COLUMNS) IN GIVEN MATRIX. C EPS - AN INPUT CONSTANT WHICH IS USED AS RELATIVE TOLERANCE C FOR TEST ON LOSS OF SIGNIFICANCE. C IER - RESULTING ERROR PARAMETER CODED AS FOLLOWS: C IER=0 - NO ERROR C IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAMETER N OR C BECAUSE SOME RADICAND IS NONPOSITIVE (MATRIX A C IS NOT POSITIVE DEFINITE, POSSIBLY DUE TO LOSS C OF SIGNIFICANCE) C IER=K - WARNING WHICH INDICATES LOSS OF SIGNIFICANCE. C THE RADICAND FORMED AT FACTORIZATION STEP K+1 C WAS STILL POSITIVE BUT NO LONGER GREATER THAN C ABS(EPS*A(K+1,K+1)). C C REMARKS C THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE C STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS. IN C THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGULAR MATRIX C IS STORED COLUMNWISE TOO. C THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL C CALCULATED RADICANDS ARE POSITIVE. C THE PRODUCT OF RETURNED DIAGONAL TERMS IS EQUAL TO THE SQUARE C ROOT OF THE DETERMINANT OF THE GIVEN MATRIX. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C SOLUTION IS DONE USING THE SQUARE-ROOT METHOD OF CHOLESKY. C THE GIVEN MATRIX IS REPRESENTED AS THE PRODUCT OF 2 TRIANGULAR C MATRICES, WHERE THE LEFT HAND FACTOR IS THE TRANSPOSE OF THE C RETURNED RIGHT HAND FACTOR. C DIMENSION A(*) DOUBLE PRECISION DPIV,DSUM C C TEST ON WRONG INPUT PARAMETER N C IF(N-1)12,1,1 1 IER=0 C C INITIALIZE DIAGONAL-LOOP C KPIV=0 DO 11 K=1,N KPIV=KPIV+K IND=KPIV LEND=K-1 C C CALCULATE TOLERANCE C TOL=ABS(EPS*A(KPIV)) C C START FACTORIZATION-LOOP OVER K-TH ROW C DO 11 I=K,N DSUM=0.D0 IF(LEND)2,4,2 C C START INNER LOOP C 2 DO 3 L=1,LEND LANF=KPIV-L LIND=IND-L 3 DSUM=DSUM+DBLE(A(LANF)*A(LIND)) C C END OF INNER LOOP C C TRANSFORM ELEMENT A(IND) C 4 DSUM=DBLE(A(IND))-DSUM IF(I-K)10,5,10 C C TEST FOR NEGATIVE PIVOT ELEMENT AND FOR LOSS OF SIGNIFICANCE. C 5 IF(SNGL(DSUM)-TOL)6,6,9 6 IF(DSUM)12,12,7 7 IF(IER)8,8,9 8 IER=K-1 C C COMPUTE PIVOT ELEMENT C 9 DPIV=DSQRT(DSUM) A(KPIV)=DPIV DPIV=1.D0/DPIV GO TO 11 C C CALCULATE TERMS IN ROW C 10 A(IND)=DSUM*DPIV 11 IND=IND+I C C END OF DIAGONAL-LOOP C RETURN 12 IER=-1 RETURN END C--------------------------------------------------- SUBROUTINE MDCH(CS,DF,P,IER) C--------------------------------------------------- C C PURPOSE C CHI-SQUARED PROBABILITY DISTRIBUTION FUNCTION C C DESCRIPTION OF PARAMETERS C CS - INPUT VALUE FOR WHICH THE PROBABILITY IS C COMPUTED. CS MUST BE GREATER THAN OR EQUAL C TO ZERO. C DF - INPUT NUMBER OF DEGREES OF FREEDOM OF THE C CHI-SQUARED DISTRIBUTION. DF MUST BE GREATER THAN C OR EQUAL TO .5 AND LESS THAN OR EQUAL TO 200,000. C P - OUTPUT PROBABILITY THAT A RANDOM VARIABLE WHICH C FOLLOWS THE CHI-SQUARED DISTRIBUTION WITH DF C DEGREES OF FREEDOM IS LESS THAN OR EQUAL TO CS. C IER - OUTPUT ERROR CODE C =129 INDICATES THAT CS OR DF WAS SPECIFIED C INCORRECTLY. C =34 INDICATES THAT THE NORMAL PDF WOULD HAVE C PRODUCED AN UNDERFLOW C C REMARKS C FOR DEGREES OF FREEDOM GREATER THAN MAXDF=100, C THE NORMAL APPROXIMATION IS USED. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NDTR, MGAMAD C C METHOD C SEE EQUATION 26.4.14, 6.5.32, 6.5.29 , HANDBOOK OF C MATHEMATICAL FUNCTIONS, ABRAMOWITZ AND STEGUN, C DOVER PUBLICATIONS, INC., NEW YORK. C REAL CS,DF,P REAL PT2 DOUBLE PRECISION A,Z,DGAM,EPS,W,W1,B,Z1,HALF,ONE,THRTEN,THRD DOUBLE PRECISION MGAMAD c REAL X,C REAL X DATA EPS/1.0D-6/,HALF/5.D-1/,THRTEN/13.D0/,ONE/1.D0/ DATA PT2/.2222222E0/ DATA THRD/.3333333333333333D0/ FUNC(W,A,Z)=W*DEXP(A*DLOG(Z)-Z) C C FIRST EXECUTABLE STATEMENT C TEST FOR INVALID INPUT VALUES C IF (DF .GE. .5 .AND. DF .LE. 2.E5 .AND. CS .GE. 0.0) GO TO 5 IER=129 P=-1 GO TO 9000 5 IER=0 C C SET P=0. IF CS IS LESS THAN OR EQUAL TO 10.**(-12) C IF (CS .GT. 1.E-12) GO TO 15 10 P=0.0 GO TO 9005 15 IF (DF .LE. 100.) GO TO 20 C C USE NORMAL DISTRIBUTION APPROXIMATION FOR LARGE C DEGREES OF FREEDOM C IF (CS .LE. 2.0) GO TO 10 X=((CS/DF)**THRD-(ONE-PT2/DF))/SQRT(PT2/DF) IF (X .GT. 5.0) GO TO 50 IF (X .LT. -18.8055) GO TO 55 CALL NDTR(X,P,DEN) GO TO 9005 C C INITIALIZATION FOR CALCULATION USING INCOMPLETE C GAMMA FUNCTION C 20 IF (CS .GT. 200.) GO TO 50 A=HALF*DF Z=HALF*CS DGAM=MGAMAD(A) W=DMAX1(HALF*A,THRTEN) IF (Z .GE. W) GO TO 35 IF (DF .GT. 25. .AND. CS .LT. 2.) GO TO 10 C C CALCULATE USING EQUATION NO. 6.5.29 C W=ONE/(DGAM*A) W1=W DO 25 I=1,50 B=I W1=W1*Z/(A+B) IF (W1 .LE. EPS*W) GO TO 30 W=W+W1 25 CONTINUE 30 P=FUNC(W,A,Z) GO TO 9005 C C CALCULATE USING EQUATION NO. 6.5.32 C 35 Z1=ONE/Z B=A-ONE W1=B*Z1 W=ONE+W1 DO 40 I=2,50 B=B-ONE W1=W1*B*Z1 IF (W1. LE. EPS*W) GO TO 45 W=W+W1 40 CONTINUE 45 W=Z1*FUNC(W,A,Z) P=ONE-W/DGAM GO TO 9005 50 P=1.0 GO TO 9005 55 P=0.0 IER=34 9000 CONTINUE 9005 RETURN END C----------------------------------------------- DOUBLE PRECISION FUNCTION MGAMAD(X) C----------------------------------------------- C C PURPOSE C EVALUATE THE GAMMA FUNCTION OF A POSITIVE DOUBLE PRECISION C ARGUMENT C C DESCRIPTION OF PARAMETERS C X -- INPUT DOUBLE PRECISION ARGUMENT (X MUST BE A C POSITIVE REAL NUMBER) C C METHOD C SEE CODY, W. J. AND HILLSTROM, K. E. - CHEBYSHEV APPROXIMATIONS C FOR THE NATURAL LOGARITHM OF THE GAMMA FUNCTION, MATHEMATICS C OF COMPUTATION, 21(89) 1967, 198-208; ALSO SEE EQUATION 6.1.41, C HANDBOOK OF MATHEMATICAL FUNCTIONS, ABRAMOWITZ AND STEGUN, C DOVER PUBLICATIONS, INC., NEW YORK. C DOUBLE PRECISION P(6,3),Q(6,3),PSUM,QSUM,DLNGA,X DATA PI/3.14159265358979D0/ DATA P/-2.6094066054623D0,-4.1509018875434D01, 1 -9.5564117677317D01,-1.1811439967596D01, 2 2.8113744347038D01,3.5173589912443D0, 3 -2.16192292624703D02,-8.27790897809598D02, 4 1.82987822012009D02,7.06543700154966D02, 5 1.49903662709861D02,4.54827477723909D0, 6 -2.4273113085758D06,1.3860869828508D06, 7 1.8537773351564D06,-6.4279927530351D05, 8 -1.5515971577126D05,-4.2105209252847D03/ DATA Q/5.4587504274950D-01,1.6674969701154D01, 1 7.8556098036754D01,8.7289305773548D01, 2 2.3590762639739D01,1.0D0, 3 1.16412659461333D02,1.20459293663292D03, 4 1.85645035686087D03,7.05287069715149D02, 5 6.65573507467416D01,1.0D0, 6 -1.0542482321634D06,-2.4515705199457D06, 7 -8.6274186723037D05,-6.7771258633073D04, 8 -9.4136613234388D02,1.0D0/ I=1 C INITIALIZE FUNCTION VALUE MGAMAD=0.0D0 IF (X.LE.0) RETURN C C IF X>12, USE EQUATION 6.1.41 C IF (X.GT.12.0) GO TO 60 C C USE COEFFICIENTS P(J,I) AND Q(J,I) DEPENDING ON C VALUES OF X C IF ((X.GE.4.0).AND.(X.LE.12.0)) I=3 IF ((X.GE.1.5).AND.(X.LE.4.0)) I=2 IF (X.LE.0.5) X=X+1.0D0 C C CALCULATE LOG-GAMMA C PSUM=0.D0 QSUM=0.D0 DO 20 J=1,5 PSUM=(PSUM+P(6-J+1,I))*X QSUM=(QSUM+Q(6-J+1,I))*X 20 CONTINUE PSUM=PSUM+P(1,I) QSUM=QSUM+Q(1,I) IF (I.EQ.2) DLNGA=(X-2.0D0)*(PSUM/QSUM) IF (I.EQ.3) DLNGA=PSUM/QSUM IF ((I.EQ.1).AND.(X.LE.0.5)) DLNGA=-DLOG(X)+(PSUM/QSUM) IF ((I.EQ.1).AND.(X.GT.0.5)) DLNGA=(X-1.0D0)*PSUM/QSUM GO TO 70 C C IF X>12, CALCULATE LOG-GAMMA USING EQUATION NO. 6.1.41 C 60 DLNGA=(X-0.5D0)*DLOG(X)-X+0.5D0*DLOG(2.D0*PI)+1.D0/(12.D0*X)- 1 1.D0/(360.D0*X**3)+1.D0/(1260.D0*X**5) 70 MGAMAD=DEXP(DLNGA) RETURN END C ==================================================== real Function CCV (kat,j,fr2,ac,bc,cat,X) C ==================================================== c c CCV = Compute Critical Value c computes the value of the cutpoint Uj c st Uj's are equally space b/w c Xi and Xi-1 c real X(*),ac,bc,D integer fr1,fr2,kat,j if (j.lt.2) then ps2 = 0 else call ndtr((bc*x(j-1)-ac),ps2,D) endif if (j.gt.kat-1) then ps1 = 1 else call ndtr((bc*x(j)-ac),ps1,D) endif fr1 = cat-fr2 call zdev( ((fr1*ps1+fr2*ps2)/cat),temp,D,ier) c write(*,*) ' ccv:',fr1,fr2,cat CCV=(ac+temp)/bc return end C ==================================================== real Function CzV (k,m,ac,bc,U) C ==================================================== c c CzV = Compute 'z' Value c computes the value of the cutpoint Uk c st Uk's are equally space b/w c Xi and Xi-1 c real U(*),D1,D2,ps1,ps2 real temp1,temp2,zval,ac,bc integer k,m if (k.ge.m) then D1 = 0 temp1=1 else call ndtr((bc*U(k)-ac),temp1,D1) endif if (k.lt.2) then D2=0 temp2 = 0 else call ndtr((bc*U(k-1)-ac),temp2,D2) endif zval= (ac-(D1-D2)/(temp1-temp2))/bc CzV= zval return end