COMMENT: Doug White's Hypergeometric Program, Interactive Version COMMENT: Max of 999 sample size NOTE: Check that the result is correct for table 1,2,3,4 (see "ADJUSTMENT") *3 way DIMENSION TITLE (20), FORMAT (20), RESLTS (56, 441) DIMENSION CMBTRL (12100), ICLASS (8), MSGDAT (441), ISUB (441) DIMENSION ACTUAL (50, 6), ILV (7), IUV (7) INTEGER DATATB (218, 50), PRINTQ, OUTFLE DOUBLE PRECISION NUMS (56, 2) REAL*4 x,y,z CHARACTER*14 FILE *2 way dimension f(150,999),d(2,2),r(2),c(2),h(1000) integer half,d,r,c,tables,rm,cm,sm real*8 p ! x,y C *3 way COMMON /ONE/ PRINTQ, FORMAT COMMON /TWO/ NVAR, NCASES, ITYPE, ICLASS C DATA MXNVAR, MXNCSS / 50, 218/ DATA INFILE, OUTFLE / 5, 6/ write (*,*) ' Fisher Exact Tests for 2 or 3 Way tables' write (*,*) ' (c) 1994, Douglas R. White' idev=0 call open (idev) 10 write (*,'(A\)') ' Enter 2 for 2-way Fisher exact, 3 for 3-way ex +act: ' read (*,*) nway if (nway.lt.2.or.nway.gt.3) goto 10 IF (NWAY.EQ.2) THEN write(*,*)' Enter 4 cell values for each of a series of tables:' if (idev.eq.1) write (1,'(a)')' Table A B C D p=' ELSE write(*,*)' Enter 8 cell values for each of a series of tables:' if (idev.eq.1) write (1,'(a)')' Table A1 B1 C1 D1 A2 +B2 C2 D2 p=' ENDIF ! NWAY IF (NWAY.EQ.2) THEN C READ MAXIMUM NUMBER OF CASES AND SET NUMBER OF TABLES write (*,'(a)') ' Maximum sample size, up to 999' read(*,*) max tables=250 half= max/2 write (*,'(a,i3,a)') ' Computing', half, ' N choose K tables ...' C COMPUTE COMBINATORIALS USING LOG 10 do 100 k=1,half write (*,'(a\)') '.' x=k+1 f(k,k+1)=alog10(x) kk=1 do 100 m=k+2,max kk=kk+1 x=m y=kk 100 f(k,m)=f(k,m-1)+alog10(x)-alog10(y) write (*,'(a)') '.' C READ DATA CARDS: D(1,1),(1,2),(2,1),(2,2) ARE DATA CELLS if (idev.eq.1) write (*,'(a)') * ' Enter all cells A, B, C, D = 0 to end program' do 1000 ii=1,tables BEGIN write (*,'(a)') ' Enter cells A, B, C, D (-1 0 0 0 to end)' read (*,*) (d(1,i),i=1,2),(d(2,j),j=1,2) do ix=1,2 do iy=1,2 if (d(ix,iy).lt.0) stop enddo enddo C COMPUTE MARGINALS r(1)=d(1,1)+d(1,2) r(2)=d(2,1)+d(2,2) c(1)=d(1,1)+d(2,1) c(2)=d(1,2)+d(2,2) n=r(1)+r(2) if (n.gt.max) then write (*,'(a,i4)') ' Re-enter: exceeds maximum of ', max goto 1000 endif if (n.eq.0) then close (1) stop endif C FIND SMALLEST CELL, L, LESS THAN EXPECTED l=n do 200 i=1,2 do 200 j=1,2 expect=r(i)*c(j)/n if (d(i,j).gt.l.or.d(i,j).gt.expect) go to 200 k=i kk=j l=d(i,j) 200 continue C COMPUTE MARGINALS OF SMALLEST CELL LESS THAN EXPECTED rm=r(k) cm=c(kk) C COMPUTE SMALLEST OF TWO MARGINALS sm=rm if(cm.lt.sm) sm=cm C COMPUTE HYPERGEOMETRIC DISTRIBUTION k=rm if (k.gt.n/2) k=n-k call hyper(f,cm,rm,sm,n,k,h) C COMPUTE PROBABILITY OF SAME OR MORE EXTREME DISTRIBUTION p=0.0 do 300 i=1, l+1 300 p=p+h(i) C WRITE RESULTS if (idev.eq.1) + write (1,311) ii, ':', (d(1,i),i=1,2),(d(2,j),j=1,2), p write (*,311) ii, ':', (d(1,i),i=1,2),(d(2,j),j=1,2), p 311 format (i5,a,4i5,f16.12) 1000 continue endif ! NWAY if (nway.eq.3) then ii=1 C BEGIN 3 way C N.B.: THIS VERSION OF THIS PROGRAM USES THE THIRD C VERSION OF THE EQUATION FOR THE PROBABILITY DISTRIBUTION, C AND THE EXACT METHOD OF COMBINING MULTIPLE SIGNIFICANCE TESTS. C C C THESE INTRODUCTORY COMMENTS DISCUSS: C C 1) DATA DECK SET-UP C 2) CHANGING PROGRAM LIMITS C 3) THE LINEAR ARRAY, CMBTRL, CONTAINING THE COMBINATORIALS C C C DATA DECK SET-UP C C FOR EACH DATA DECK: C THREE (3) CONTROL CARDS NECESSARY C (PLACED IN FRONT OF EACH DATA DECK): C CARD #1: C COL 1-3: NUMBER OF VARIABLES C CURRENT MAXIMUM: 50 C (IF MXNVAR IS INCREASED PAST 100, DATA C PRINTOUT PROCEDURE IN SUBROUTINE RDDATA C MUST BE ALTERED.) C COL 4-6: NUMBER OF CASES (SUBJECTS) C CURRENT MAXIMUM: 218 C COL 8: TYPE OF DATA (1 = TOTALS, ANY C OTHER CHARACTER = RAW) C COL 9: PRINT DATA? (1 = YES, ANY C OTHER CHARACTER = NO) C COL 10: A FOLLOWING DATA DECK? (1 = YES, C ANY OTHER CHARACTER = NO) C CARD #2: C COL 1-80: TITLE C CARD #3: C COL 1-80: DATA CARD FORMAT C THE DATA CARD(S) FOLLOW. IF THE DATA IS RAW, INDICATE C PRESENCE OF VARIABLE WITH "1", ABSENCE WITH "0" C AND NO DATA WITH "9". ANY OTHER CHARACTER WILL C CAUSE TERMINATION OF THE PROGRAM. C IF THE DATA IS ALREADY TOTALED FOR EACH OF THE C EIGHT POSSIBLE CATEGORIES (FOR ONLY THREE C VARIABLES), THEY MUST BE IN THE FOLLOWING ORDER: C VAR. 1: 0 0 0 0 1 1 1 1 C VAR. 2: 0 0 1 1 0 0 1 1 C VAR. 3: 0 1 0 1 0 1 0 1 C IT OF COURSE MAKES NO DIFFERENCE WHICH OF YOUR THREE C VARIABLES ARE DEFINED AS VAR. 1, 2 OR 3. C C C PROGRAM LIMITS C THE PRESENT VERSION OF THIS PROGRAM OPERATES WITHIN THE C FOLLOWING LIMITS: C MAXIMUM NUMBER OF VARIABLES ACCEPTABLE = 50 C (SYMBOLIC NAME: MXNVAR) C MAXIMUM NUMBER OF CASES ACCEPTABLE = 218 C (SYMBOLIC NAME: MXNCSS) C C THE FOLLOWING INDICATES ARRAY DIMENSIONING AND MUST BE C CONSIDERED WHEN MAKING ANY OF THE CHANGES DISCUSSED ABOVE: C CMBTRL ( (MXNCSS + 2) ** 2 / 4 ) C ACTUAL (MXNVAR, 6) C DATATB (MXNCSS, MXNVAR) C NUMS (INT (MXNCSS / 4 + 2), 2) C RESLTS (INT (MXNCSS / 4) + 2, (MXNVAR - 1)*(MXNVAR - 2) / 2) C MSGDAT ((MXNVAR - 1) * (MXNVAR - 2) / 2) C ISUB ((MXNVAR - 1) * (MXNVAR - 2) / 2) C ILV (LOG (MXNCSS / 4 + 2) / LOG 2 + 1) C ULV (LOG (MXNCSS / 4 + 2) / LOG 2 + 1) C C DIMENSIONING NEED ONLY BE CHANGED IN THE MAINLINE C DIMENSION AND SPECIFICATION STATEMENTS. C C FINALLY, THE FOLLOWING INPUT AND OUTPUT DEVICES ARE C SPECIFIED IN THE THIRD DATA STATEMENT: C READ STATEMENTS: SYMBOLIC NAME: INFILE = 5 C WRITE STATEMENTS: SYMBOLIC NAME: OUTFLE = 6 C C C STORAGE OF COMBINATORIALS: C IN ORDER TO SAVE STORAGE, COMBINATORIALS ARE STORED IN THE C LINEAR ARRAY, CMBTRL. THE SIZE OF THIS ARRAY IF GIVEN BY THE C FOLLOWING EQUATION: C SIZE = INT (MXNCSS + 2) ** 2 / 4) C C THE TRANSFORMATION EQUATION IS: C CMBTRL (K) = /X\ C \Y/ C C WHERE K = INT ((X + 1) ** 2 / 4 + Y + 1) C C C C DEL WRITE (*,'(A)') ' * or Output Filename?' DEL READ (*,'(A)') file DEL IF (file(1:1).eq.'*') then DEL open (unit=6, status='unknown') DEL else DEL open (unit=OUTFLE, file=file, status='unknown') DEL endif DEL WRITE (*,'(A)') ' Input Filename?' DEL READ (*,'(A)') file DEL open (unit=INFILE, file=file, status='unknown') DEL WRITE (OUTFLE, 100) DEL100 FORMAT (1H1, 12(1H:), 4X, 40HI N T E R A C T I O N A N A L Y S DEL 2I S, 4X, 12(1H:), //, 99X, 9HVERSION 5, 5X, 18H JULY 17, 1979) C FROM 2WAY 1001 write (*,'(a)') ' Enter cells A, B, C, D, A, B, C, D (-1 0 0 0 0 +0 0 0 to end)' read (*,*) (iclass(j),j=1,8) ncases=0 do j=1,8 if (iclass(j).lt.0) stop ncases=ncases+iclass(j) enddo NVAR=3 IDECKN = 1 ! EDITED BEGIN D IDECKN = 0 LIMIT1 = (MXNCSS + 2) ** 2 / 4 C FOR EACH DATA DECK: D110 IDECKN = IDECKN + 1 D READ (INFILE, 120, END = 210) NCASES, NVAR, ITYPE, PRINTQ, MORDAT, D 2 TITLE, FORMAT D120 FORMAT (2I3, 1X, 3I1, /, 20A4,/, 20A4) LIMIT2 = (NVAR - 1) * (NVAR - 2) / 2 LIMIT3 = NCASES / 4 + 2 ALIM3 = LIMIT3 LIMIT4 = ALOG (ALIM3) / ALOG (2.) + 1 D IF (NDECKS .NE. 1) WRITE (OUTFLE, 130) IDECKN D130 FORMAT (//, 13H DECK NUMBER , I1, 1H:) D WRITE (OUTFLE, 140) TITLE, NVAR, NCASES, FORMAT, LIMIT2 D140 FORMAT (26X, 20A4, ///, 29H NUMBER OF VARIABLES IS , I3, /, D 2 29H NUMBER OF CASES IS , I3, /, D 3 29H FORMAT FOR DATA IS , 20A4, / D 4 29H MAXIMUM NUMBER OF TESTS IS , I6, //) C C CHECK PARAMETERS. D IF (NVAR .LE. MXNVAR .AND. NCASES .LE. MXNCSS) GO TO 160 D WRITE (OUTFLE, 150) MXNVAR, MXNCSS D150 FORMAT (' YOU HAVE EXCEEDED THE MAXIMUM NUMBER OF VARIABLES (', I3, D 2 ') OR CASES (', I3, ') ALLOWED.') D GO TO 170 C C READ DATA, CHECK AND STORE IN DATATB. D160 CALL RDDATA (INFILE, OUTFLE, DATATB, IBDDTA, MXNCSS) D IF (IBDDTA .EQ. 1) GO TO 170 C C COMPUTE COMBINATORIALS UP TO MXNCSS TAKEN MXNCSS/2 AT A TIME, AND STORE C IN CMBTRL - BUT ONLY FOR THE FIRST DATA DECK. BEGIN STOP HERE TO CHECK COMBINATORIALS IF (IDECKN .EQ. 1) CALL CALCMB (MXNCSS, CMBTRL, LIMIT1) C C CALCULATE ACTUAL AND EXPECTED HYPERGEOMETRIC PROBABILITIES. iclass passed in common CALL PROB (MXNCSS, DATATB, CMBTRL, LIMIT1, ACTUAL, NUMS, 2 LIMIT3, OUTFLE, LIMIT2, RESLTS, MSGDAT, ISUB, LIMIT4, ILV, IUV, 3 CMBPRB) C C PRINT ACTUAL AND EXPECTED TABLES. if (idev.eq.1) + write (1,310) ii, ':', (iclass(j),j=1,8), cmbprb write (*,310) ii, ':', (iclass(j),j=1,8), cmbprb 310 format (i5,a,8i5,f12.8) need do loop here goto 1001 endif ! NWAY DEL CALL PRNTBS (OUTFLE, NVAR, ACTUAL, ITYPE) D170 IF (MORDAT .EQ. 1) GO TO 110 C C THAT'S ALL FOLKS! D WRITE (OUTFLE, 180) D180 FORMAT (////, 55X, 4(1H-), 2X, 10HEND OF RUN, 2X, 4(1H-)) STOP D190 WRITE (OUTFLE, 201) D201 FORMAT (' YOU DIDN''T PUT ANY CARDS IN!!') D210 WRITE (OUTFLE, 220) D220 FORMAT (' WHERE''S YOUR CONTROL CARDS?!!') D STOP END C C C D end subroutine hyper(f,cm,rm,sm,n,k,h) integer cm,rm,sm real*8 x,y dimension h(1000), f(150,999) h(1)=0.0 j=rm if (j.gt.(n-cm)/2) j=n-cm-j x=f(j,n-cm)-f(k,n) if (x.gt.-20.0) h(1)=10**x do 100 ii=2,sm+1 i=ii-1 j=rm-i ADJUSTMENT OF PROGRAM J>0 (J=0 for 1,2,3,4 etc.) if (j.gt.0) then x=f(i,cm)+f(j,n-cm)-f(k,n) h(ii)=0.0 if (x.gt.-20.0) h(ii)=10**x endif 100 continue h(sm+1)=0.0 x=0.0 y=0.0 if (sm.lt.cm) x=f(sm,cm) if (sm.lt.rm) y=f(rm-sm,n-cm) x=x+y-f(k,n) if (x.gt.-20.0) h(sm+1)=10**x return end subroutine open (idev) character*14 file character*1 Yes logical ex write (*,'(3(a))') ' Do you want to write results to a file?' 30 read (*,'(a1)') Yes if (Yes.eq.'Y'.or.Yes.eq.'y') goto 40 return 40 write (*,'(a\)') ' Filename? ' read (*,'(a)') File inquire (file=file, exist=ex) if (ex) goto 100 50 idev=1 open (idev, file=file) return 100 write (*,'(3(a))') ' File ', file, ' found - write over it?' read (*,'(a1)') Yes if (Yes.eq.'Y'.or.Yes.eq.'y') goto 50 write (*,'(3(a))') ' Write results to another file?' goto 30 end SUBROUTINE CALCMB (MXNCSS, CMBTRL, LIMIT1) DIMENSION CMBTRL (LIMIT1) C K = 1 DO 110 I = 1, MXNCSS + 1 X = I - 1 CMBTRL (K) = 0 K = K + 1 IF (X .LE. 1) GO TO 110 DO 100 Y = 1, X/2 Z = X - Y + 1 CMBTRL (K) = CMBTRL (K - 1) - ALOG10 (Y) + ALOG10 (Z) 100 K = K + 1 110 CONTINUE RETURN END C SUBROUTINE PROB (MXNCSS, DATATB, CMBTRL, LIMIT1, ACTUAL, NUMS, 2 LIMIT3, OUTFLE, LIMIT2, RESLTS, MSGDAT, ISUB, LIMIT4, ILV, IUV, 3 CMBPRB) COMMON /TWO/ NVAR, NCASES, ITYPE, ICLASS INTEGER OUTFLE, DATATB (NCASES, NVAR), CTNGCY (3, 2, 2) DIMENSION ICLASS (8), MRGNLS (2, 2, 2), MSGDAT (LIMIT2) DIMENSION IPOSTB (2, 2), CMBTRL (LIMIT1), ISUB (LIMIT2) DIMENSION RESLTS (LIMIT3, LIMIT2), ACTUAL (NVAR, 6) DIMENSION ILV (LIMIT4), IUV (LIMIT4) DOUBLE PRECISION NUMS (LIMIT3, 2), SMPSPC, DEN, DNUM, TOTSMP, 2 SMPLSP C C FOR EACH TRIPLE: D DO 190 IVAR = 1, NVAR D NUMTBS = 0 D IHLDX2 = 0 D TOTSMP = 0. D NMTSTS = 0 NMTSTS = 1 D DO 186 JVAR = 1, NVAR - 1 D IF (JVAR .EQ. IVAR) GO TO 186 D DO 185 KVAR = JVAR + 1, NVAR D IF (KVAR .EQ. IVAR) GO TO 185 D IF (ITYPE .EQ. 1) GO TO 125 D DO 110 I = 1, 8 D110 ICLASS (I) = 0 C C THERE ARE EIGHT PERMUTATIONS OF THREE DICHOTOMIZED VARIABLES. C THIS ALGORITHM CLASSIFIES AND COUNTS THE PARTICULAR COMBINATION OF C IVAR, JVAR AND KVAR FOR EACH CASE USING A BINARY-TO-DECIMAL C TRANSFORMATION. D DO 120 I = 1, NCASES D J = 4 * DATATB(I, IVAR) + 2 * DATATB(I,JVAR) + DATATB(I,KVAR) + 1 D120 IF (J .LE. 8) ICLASS (J) = ICLASS (J) + 1 C C INITIALIZE NECESSARY TABLES TO 0. 125 CTNGCY (3, 1, 1) = 0 CTNGCY (3, 1, 2) = 0 CTNGCY (3, 2, 1) = 0 CTNGCY (3, 2, 2) = 0 ITOTAL = 0 C C CONSTRUCT RAW AND CONTROLLED CONTINGENCY TABLES, IN CTNGCY, AND C CONTROLLED MARGINALS, IN MRGNLS, WITH 1ST SUBSCRIPT = 1 AS C ICNVAR PRESENT, = 2 AS ICNVAR ABSENT AND = 3 AS RAW, USING A DO 130 ICPORA = 1, 2 DO 130 IROW = 1, 2 DO 130 ICOL = 1, 2 K = 15 - 4 * ICPORA - 2 * IROW - ICOL CTNGCY (ICPORA, IROW, ICOL) = ICLASS (K) 130 CTNGCY (3, IROW, ICOL) = CTNGCY (3, IROW, ICOL) + 2 CTNGCY (ICPORA, IROW, ICOL) C C DETERMINE MARGINALS: MRGNLS (1, 1, 1) = CTNGCY (1, 1, 1) + CTNGCY (1, 1, 2) MRGNLS (1, 1, 2) = CTNGCY (1, 2, 1) + CTNGCY (1, 2, 2) MRGNLS (1, 2, 1) = CTNGCY (1, 1, 1) + CTNGCY (1, 2, 1) MRGNLS (1, 2, 2) = CTNGCY (1, 1, 2) + CTNGCY (1, 2, 2) MRGNLS (2, 1, 1) = CTNGCY (2, 1, 1) + CTNGCY (2, 1, 2) MRGNLS (2, 1, 2) = CTNGCY (2, 2, 1) + CTNGCY (2, 2, 2) MRGNLS (2, 2, 1) = CTNGCY (2, 1, 1) + CTNGCY (2, 2, 1) MRGNLS (2, 2, 2) = CTNGCY (2, 1, 2) + CTNGCY (2, 2, 2) C C FIND THE SMALLEST MARGINAL IN THE CONTROLLED TABLES. ISMMRG = MRGNLS (1, 1, 1) ICPORA = 1 IROW = 1 ICOL = 1 DO 140 I = 1, 2 DO 140 J = 1, 2 DO 140 K = 1, 2 IF (ISMMRG .LT. MRGNLS (I, J, K)) GO TO 140 ISMMRG = MRGNLS (I, J, K) ICPORA = I IROW = J ICOL = K 140 CONTINUE C C IF THE SMALLEST MARGINAL IS 0, DON'T CONSIDER TABLE. IF (ISMMRG .EQ. 0) GO TO 185 C ESTABLISH WHICH CELL IS TO BE USED AS THE BASE CELL, CHOOSING ONE C ASSOCIATED WITH THE SMALLEST MARGINAL. IRBASE = IROW ICBASE = ICOL IF (ICOL .EQ. 2) IRBASE = 3 - IRBASE C C EXAMINE ALL TABLES POSSIBLE GIVEN THE SMALLEST MARGINAL IN LIGHT C OF THE RAW CELL VALUES. MINVAL = -1 SMPSPC = 0 IHLDEX = 0 DO 150 ID = 1, ISMMRG + 1 I = ID - 1 C C FILL IN TABLE VALUES. IPOSTB (IRBASE, ICBASE) = I IF (I .GT. CTNGCY (3, IRBASE, ICBASE)) GO TO 150 IPOSTB (IRBASE, 3 - ICBASE) = MRGNLS (ICPORA, 1, IRBASE) - I IF (IPOSTB (IRBASE, 3 - ICBASE) .GT. CTNGCY (3, IRBASE, 3 - 2 ICBASE)) GO TO 150 IPOSTB (3 - IRBASE, ICBASE) = MRGNLS (ICPORA, 2, ICBASE) - I IF (IPOSTB (3 - IRBASE, ICBASE) .GT. CTNGCY (3, 3 - IRBASE, 2 ICBASE)) GO TO 150 IPOSTB (3 - IRBASE, 3 - ICBASE) = MRGNLS (ICPORA, 1, 3 - IRBASE) - 2 IPOSTB (3 - IRBASE, ICBASE) IF (IPOSTB (3 - IRBASE, 3 - ICBASE) .GT. CTNGCY (3, 3 - IRBASE, 2 3 - ICBASE)) GO TO 150 IF (MINVAL .EQ. -1) MINVAL = I IX = ID - MINVAL NUMS (IX, 1) = 0. NUMS (IX, 2) = IX DO 145 II = 1, 2 DO 145 JJ = 1, 2 ICELL = IPOSTB (II,JJ) IF (ICELL .GT. CTNGCY (3,II,JJ)/2) ICELL = CTNGCY (3,II,JJ) - 2 ICELL K = (CTNGCY (3, II, JJ) + 1) ** 2 / 4 + ICELL+ 1 145 NUMS (IX, 1) = NUMS (IX, 1) + CMBTRL (K) C C CALCULATE SAMPLE SPACE, USING OVERFOLW (10**35) CONTROL. IF ((NUMS (IX, 1) - IHLDEX) .LT. 35.) GO TO 146 IADD = NUMS (IX, 1) - IHLDEX - 34 SMPSPC = SMPSPC / 10. ** IADD IHLDEX = IHLDEX + IADD 146 DNUM = 10. ** (NUMS (IX, 1) - IHLDEX) IF ((10. ** 35. - SMPSPC) .GT. DNUM) GO TO 147 SMPSPC = SMPSPC / 10. DNUM = DNUM / 10. IHLDEX = IHLDEX + 1 GO TO 146 147 SMPSPC = SMPSPC + DNUM MAXVAL = I 150 CONTINUE C C DETERMINE RANGE OF TABLES POSSIBLE GIVEN ALL RESTRAINTS. IRANGE = MAXVAL - MINVAL + 1 C C IF IRANGE IS 1, DON'T CONSIDER TABLE. IF (IRANGE .EQ. 1) GO TO 185 DEL NMTSTS = NMTSTS + 1 MSGDAT (NMTSTS) = ITOTAL NUMTBS = NUMTBS + IRANGE C C CALCULATE TOTAL SAMPLE SPACE, USING OVERFLOW (10**35) CONTROL. SMPLSP = SMPSPC IF (IHLDX2 .GT. IHLDEX) GO TO 1501 IADD = IHLDEX - IHLDX2 TOTSMP = TOTSMP / 10. ** IADD IHLDX2 = IHLDEX GO TO 1502 1501 IADD = IHLDX2 - IHLDEX SMPLSP = SMPLSP / 10. ** IADD 1502 IF ((10. ** 35. - TOTSMP) .GT. SMPLSP) GO TO 1503 TOTSMP = TOTSMP / 10. SMPLSP = SMPLSP / 10. IHLDX2 = IHLDX2 + 1 GO TO 1502 1503 TOTSMP = TOTSMP + SMPLSP C C CALCULATE PROBABILITIES AND INCREMENT ACTUAL TABLE. ISTOP = CTNGCY (ICPORA, IRBASE, ICBASE) - MINVAL + 1 DEN = DLOG10 (SMPSPC) + IHLDEX C C SORT POINTS. IF (IRANGE .LE. 16) CALL MYSORT (NUMS, IRANGE, LIMIT3) IF (IRANGE .GT. 16) CALL QSORT (NUMS, IRANGE, LIMIT3, LIMIT4, ILV, 2 IUV) C C CALCULATE POINT PROBABILITIES. 158 IFLAG = 0 CMBPRB = 0. DO 180 I = 1, IRANGE PRBLOG = NUMS (I, 1) - DEN PTPROB = 0. debug IF (I+2.GT.LIMIT3.or.NMTSTS.GT.LIMIT2) THEN write (*,*) I, NMTSTS ENDIF RESLTS (I + 2, NMTSTS) = 10 ** PRBLOG IF (IFLAG .EQ. 0) CMBPRB = CMBPRB + 10 ** PRBLOG IF (ISTOP .EQ. NUMS (I, 2)) RESLTS (1, NMTSTS) = I IF (I .EQ. IRANGE) GO TO 180 IF (ISTOP .EQ. NUMS (I, 2) .AND. ISTOP .NE. NUMS (I + 1, 2)) 2 IFLAG = 1 IF (NUMS (I, 1) .NE. NUMS (I+1, 1)) GO TO 180 NUMS (I + 1, 2) = NUMS (I, 2) 180 CONTINUE RESLTS (2, NMTSTS) = IRANGE 185 CONTINUE 186 CONTINUE C C COMBINE TEST. RETURN DEL IF (ITYPE .NE. 1.AND.NMTSTS.GT.0) CALL CMBTST (RESLTS, NMTSTS, DEL 2 CMBPRB, LIMIT3, ISUB) AVMGDT = 0. DO 187 I = 1, NMTSTS 187 AVMGDT = AVMGDT + MSGDAT (I) IF (NMTSTS.GT.0) THEN AVMGDT = AVMGDT / NMTSTS ACTUAL (IVAR, 1) = CMBPRB ACTUAL (IVAR, 2) = NMTSTS ACTUAL (IVAR, 3) = NUMTBS ACTUAL (IVAR, 4) = AVMGDT ACTUAL (IVAR, 5) = TOTSMP ACTUAL (IVAR, 6) = IHLDX2 IF (Actual(IVAR,2).eq.0) Actual(iVAR,1)=1.0 WRITE (OUTFLE, 300) IVAR, (ACTUAL (IVAR, J), J = 1, 6) 300 FORMAT (4X, I3, 4X, F11.9, 3X, F4.0, 5X, F5.0, 4X, F4.0, 1PE16.7, 2 ' * 10 ** ', F8.0) prntbs ELSE ACTUAL (IVAR, 1) = 0 ACTUAL (IVAR, 2) = 0 ACTUAL (IVAR, 3) = 0 ACTUAL (IVAR, 4) = 0 ACTUAL (IVAR, 5) = 0 ACTUAL (IVAR, 6) = 0 ENDIF 190 CONTINUE RETURN END C C C SUBROUTINE MYSORT (NUMS, IRANGE, LIMIT3) DOUBLE PRECISION NUMS (LIMIT3, 2), HOLD1, HOLD2 C DO 300 I = 2, IRANGE HOLD1 = NUMS (I, 1) HOLD2 = NUMS (I, 2) J = I - 1 100 IF (HOLD1 .GE. NUMS (J, 1)) GO TO 200 NUMS (J + 1, 1) = NUMS (J, 1) NUMS (J + 1, 2) = NUMS (J, 2) J = J - 1 IF (J .GT. 0) GO TO 100 200 NUMS (J + 1, 1) = HOLD1 NUMS (J + 1, 2) = HOLD2 300 CONTINUE RETURN END C C C SUBROUTINE QSORT (NUMS, IRANGE, LIMIT3, LIMIT4, ILV, IUV) DOUBLE PRECISION NUMS (LIMIT3, 2), HOLD (2) DIMENSION ILV (LIMIT4), IUV (LIMIT4) C C INITIALIZE PUSHDOWN LIST. ILV (1) = 1 IUV (1) = IRANGE I = 1 100 IF (I .LT. 1) RETURN C C PARTITION ONE SEGMENT. 110 ILP = ILV (I) IUP = IUV (I) J = IUP - ILP + 1 IF (J .GE. 2) GO TO 120 I = I - 1 GO TO 100 120 IF (J .NE. 2) GO TO 150 IF (NUMS (ILP, 1) .LE. NUMS (IUP, 1)) GO TO 140 DO 130 II = 1, 2 HOLD (II) = NUMS (ILP, II) NUMS (ILP, II) = NUMS (IUP, II) 130 NUMS (IUP, II) = HOLD (II) 140 I = I - 1 GO TO 100 C C CHOOSE BOUND. 150 K = (ILP + IUP) / 2 IF (J .LT. 10) GO TO 190 IF (NUMS (ILP, 1) .GT. NUMS (K, 1)) GO TO 170 IF (NUMS (K, 1) .LE. NUMS (IUP, 1)) GO TO 190 IF (NUMS (ILP, 1) .GT. NUMS (IUP, 1)) GO TO 160 K = IUP GO TO 190 160 K = ILP GO TO 190 170 IF (NUMS (IUP, 1) .LE. NUMS (K, 1)) GO TO 190 IF (NUMS (ILP, 1) .GT. NUMS (IUP, 1)) GO TO 180 K = ILP GO TO 190 180 K = IUP 190 DO 200 II = 1, 2 HOLD (II) = NUMS (K, II) 200 NUMS (K, II) = NUMS (IUP, II) ILP = ILP - 1 C C MOVE LOWER POINTER. 210 IF ((IUP - ILP) .LT. 2) GO TO 240 ILP = ILP + 1 IF (NUMS (ILP, 1) .LE. HOLD (1)) GO TO 210 NUMS (IUP, 1) = NUMS (ILP, 1) NUMS (IUP, 2) = NUMS (ILP, 2) C C MOVE UPPER POINTER. 220 IF ((IUP - ILP) .LT. 2) GO TO 230 IUP = IUP - 1 IF (NUMS (IUP, 1) .GT. HOLD (1)) GO TO 220 NUMS (ILP, 1) = NUMS (IUP, 1) NUMS (ILP, 2) = NUMS (IUP, 2) GO TO 210 230 IUP = IUP - 1 C C TEST FOR EQUAL ITEMS. 240 IF (IUP .NE. IUV (I)) GO TO 300 ILP = ILV (I) - 1 250 IF ((IUP - ILP) .LT. 2) GO TO 280 ILP = ILP + 1 IF (NUMS (ILP, 1) .LT. HOLD (1)) GO TO 250 NUMS (IUP, 1) = NUMS (ILP, 1) NUMS (IUP, 2) = NUMS (ILP, 2) 260 IF ((IUP - ILP) .LT. 2) GO TO 270 IUP = IUP - 1 IF (NUMS (IUP, 1) .GE. HOLD (1)) GO TO 260 NUMS (ILP, 1) = NUMS (IUP, 1) NUMS (ILP, 2) = NUMS (IUP, 2) GO TO 250 270 IUP = IUP - 1 280 NUMS (IUP, 1) = HOLD (1) NUMS (IUP, 2) = HOLD (2) IF (IUP .NE. ILV (I)) GO TO 290 I = I - 1 GO TO 100 290 IUV (I) = IUP - 1 GO TO 110 C C IF UPPER SEGMENT IS NOT EMPTY, STORE BOTH SEGMENTS. 300 NUMS (IUP, 1) = HOLD (1) NUMS (IUP, 2) = HOLD (2) IF ((IUP - ILV (I)) .GE. (IUV (I) - IUP)) GO TO 310 ILV (I + 1) = ILV (I) IUV (I + 1) = IUP - 1 ILV (I) = IUP + 1 GO TO 320 310 ILV (I + 1) = IUP + 1 IUV (I + 1) = IUV (I) IUV (I) = IUP - 1 320 I = I + 1 GO TO 110 END C C C