Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-05 - decus/20-0137/stp/stp11.for
There is 1 other file named stp11.for in the archive. Click here to see a list.
C                                            *** STAT PACK ***
C     SUBROUTINE FOR CHI SQUARES
C     CALLING SEQUENCE: CALL CHI(NV,NC,MV,MC,DATA,ISUB,LEVEL,NAMES)
C     WHERE NV - NUMBER OF VARIABLES ACTUALLY USED
C           NC - NUMBER OF OBSERVATION ACTUALLY USED
C           MV - MAX. NUMBER OF VARIABLES POSS. (SPEC. IN MAIN)
C           MC - MAX. NUMBER OF OBSERVATIONS POSS. (SPEC. IN MAIN)
C           DATA - STORAGE FOR DATA DIMENSIONED FOR MAXIMUM
C           ISUB IS A VECTOR DIMENSIONED AT LEAST FOR NC
C           LEVEL IS A VECTOR DIMENSIONED AT LEAST FOR NC
C           NAMES - IS A VECTOR CONTAINING VARIABLE NAMES
C
C      CHI SQUARE OPTIONS ARE TO CALAPSE THE CONTINGENCY TABLE, AND
C     TO DETERMINE CERTAIN GROUPINGS OF CELLS PREVIOUS TO COLAPSEING
C     THE ROUTINE MAKES USE OF CHITAB, CHISOR, AND FISHER.
C
      SUBROUTINE CHI(NV,NC,MV,MC,DATA,ISUB,LEVEL,NAMES)
      DIMENSION DATA(MC,MV),ISUB(1),LEVEL(1),IV(2),ANS(7),BRK(2,3,15)
      DIMENSION LVL(2),MCC(3,3),IVB(2)
      DIMENSION NAMES(1),STRG(50),B(5),LTBM(20)
      COMMON/DEV/ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON/EXTRA/HEDR(70),NSZ
      ILMAP="007777700000
      IRMAP="000000077777
      IJUST="000000100000
1     BREAK=0
      CONTG=0
      COLPS=0
      FISHR=0
      LVL(2)=0
      IF(ICC.NE.2) WRITE(IDLG,3)
3     FORMAT(' LIST OPTIONS SEPARATED BY COMMAS'/)
      READ(ICC,4)ANS
4     FORMAT(7(A5,1X))
      IF (ANS(1).EQ.'!') RETURN
      DO 5 J=1,7
      IF (ANS(J).EQ.' ') GO TO 12
      IF (ANS(J).NE.'HELP') GO TO 7
      WRITE(IDLG,6)
6     FORMAT('0THE CHI SQUARE COMMAND OPERATES FROM RAW DATA NOT'/
     1 ' PRE CALCULATED TABLES, OPTIONS AVAILABLE ARE:'/
     2 ' "GROUP" - USED TO ESTABLISH GROUPINGS PRIOR TO COLPS'/
     3' "CONTG" - USED TO ELIMINATE CONTINGENCY TABLE FROM FINAL OUTPUT'
     4/ ' "COLPS" - COLLAPSE CONTINGENCY TABLE (INCLUDES OMITTING)'/
     5' "FISHR" - FISHER''S EXACT TEST PROBABILITY (2X2 ONLY)'/
     6'0IF NO OPTIONS ARE DESIRED TYPE A RETURN')
      GO TO 1
7     IF (ANS(J).NE.'GROUP') GO TO 8
      BREAK=1
      GO TO 5
8     IF (ANS(J).NE.'CONTG') GO TO 9
      CONTG=1
      GO TO 5
9     IF (ANS(J).NE.'COLPS') GO TO 10
      COLPS=1
      GO TO 5
10    IF(ANS(J).NE.'FISHR') GO TO 19
      FISHR=1
      GO TO 5
19    WRITE(IDLG,11) ANS(J)
11    FORMAT(' OPTION "',A5,'" DOES NOT EXIST')
      GO TO 1
5     CONTINUE
12    IF(ICC.NE.2) WRITE(IDLG,13)
13    FORMAT(' WHICH VARIABLE OR VARIABLES? ',$)
      IRET=-98
      CALL ALPHA(IV,2,NN,IRET,IHELP,IERR,NAMES,NV)
      IF (IRET.EQ.1) RETURN
      IF (IERR.EQ.1) GO TO 12
      IF (IHELP.EQ.1) GO TO 12
      IVB(1)=IV(1)
      IVB(2)=IV(2)
      IF (IV(1).LT.0) IVB(1)=1
      IF (NN.LT.2) GO TO 18
      IF (IV(2).LT.0) IVB(2)=1
      IF ((IV(1).LT.0).AND.(IV(2).LT.0)) IVB(2)=2
      GO TO 18
15    IF (IV(2).LT.0) GO TO 16
      IF (IV(1).LT.0) GO TO 17
      RETURN
16    IVB(2)=IVB(2)+1
      IF (IVB(2).LE.NV) GO TO 18
      IF(IV(1).GE.0) RETURN
      IVB(1)=IVB(1)+1
      IVB(2)=IVB(1)+1
      IF(IVB(2).LE.NV) GO TO 18
      RETURN
17    IVB(1)=IVB(1)+1
      IF (IVB(1).LE.NV) GO TO 18
      RETURN
18    DO 2 I=1,2
2     LVL(I)=0
      IF (NN.LT.2) GO TO 20
      IF (IVB(1).EQ.IVB(2)) GO TO 15
C
C      HERE BEGINS ACTUAL CHI SQUARE
C
20    IF (BREAK.NE.1) GO TO 35
      DO 33 I=1,NN
      IF(ICC.NE.2) WRITE(IDLG,22) NAMES(IVB(I))
22    FORMAT(' ENTER RANGES, 1 PER LINE FOR VARIABLE: ', A5/)
      DO 23 J=1,15
27    IF(ICC.NE.2) WRITE(IDLG,24)
24    FORMAT('+?',$)
      READ(ICC,25,END=21) ANS(1)
25    FORMAT(A5)
      IF (ANS(1).EQ.'!') RETURN
      IF (ANS(1).EQ.'STOP') GO TO 21
      IF (ANS(1).EQ.' ') GO TO 21
      IF(ANS(1).EQ.'AUTO') GO TO 32
      IF (ANS(1).NE.'HELP') GO TO 28
      WRITE(IDLG,26)
26    FORMAT(' ENTER EACH RANGE AFTER THE ? , SMALLER FIRST SEPARATED
     1 BY A COMMA, ONE PER LINE'/
     2 ' EXAMPLE:'/' ?13,2'/' TO SIGNAL THE END OF RANGES TYPE
     3 "STOP", ^Z, OR <CR>'/)
      GO TO 27
32    LVL(I)=0
      GO TO 33
28    REREAD 29, (BRK(I,K,J),K=1,2)
29    FORMAT(2F)
      IF (BRK(I,1,J).LE.BRK(I,2,J)) GO TO 23
      WRITE(IDLG,30)
30    FORMAT(' SMALLER FIRST'//)
      GO TO 27
23    BRK(I,3,J)=0
      WRITE(IDLG,31)
31    FORMAT(' 15 RANGES MAXIMUM')
      J=16
21    LVL(I)=J-1
33    CONTINUE
C
C     SET LEVELS FOR EACH SET OF VARIABLES
C
35    DO 34 I=1,NC
      ISUB(I)=I
34    LEVEL(I)=0
      NK=NC
      DO 36 I=1,NN
      KONS=1
      IF(I.EQ.2) KONS=IJUST
      IF (LVL(I).EQ.0) GO TO 39
      DO 37 J=1,NK
      DO 38 K=1,LVL(I)
      IF (DATA(ISUB(J),IVB(I)).LT.BRK(I,1,K)) GO TO 38
      IF (DATA(ISUB(J),IVB(I)).GT.BRK(I,2,K)) GO TO 38
      BRK(I,3,K)=1
      LEVEL(J)=K*KONS+LEVEL(J)
      GO TO 37
38    CONTINUE
      ISUB(J)=0
37    CONTINUE
      DO 130 J=1,LVL(I)
      IF(BRK(I,3,J).NE.0) GO TO 130
      WRITE(IDLG,131) BRK(I,1,J),BRK(I,2,J)
131   FORMAT(' RANGE ',G10.4,',',G10.4,' HAS NO VALUES - LEVEL DROPPED')
      DO 132 K=1,NK
      IF(LEVEL(K).GT.J) LEVEL(K)=LEVEL(K)-1*KONS
132   CONTINUE
130   CONTINUE
      GO TO 41
39    LL=I+2
      CALL CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,LL)
      L=1
      VALUE=DATA(ISUB(1),IVB(I))
      DO 40 J=1,NK
      IF (DATA(ISUB(J),IVB(I)).EQ.VALUE) GO TO 40
      VALUE=DATA(ISUB(J),IVB(I))
      L=L+1
40    LEVEL(J)=L*KONS+LEVEL(J)
      LVL(I)=L
41    J=0
      K=1
42    IF (ISUB(K).EQ.0) GO TO 43
      J=J+1
      ISUB(J)=ISUB(K)
      LEVEL(J)=LEVEL(K)
43    K=K+1
      IF (K.LE.NK) GO TO 42
      NK=J
36    CONTINUE
      IF(NK.LT.1) RETURN
C
C     LEVELS DETERMINED NOW
C
      IF (COLPS.NE.1) GO TO 100
C
C     COLLAPSING PORTION
C
      CALL CHITAB(NK,LEVEL,ISUB,MV,MC,DATA,NAMES,IVB,1,NN,LINES)
      IF(ICC.NE.2) WRITE(IDLG,50)
50    FORMAT('0COLLAPSING PORTION INSERT 1 INSTRUCTION PER LINE AFTER
     1 THE ?')
51    IF(ICC.NE.2) WRITE(IDLG,52)
52    FORMAT(' ? ',$)
      READ(ICC,53,END=95) STRG
53    FORMAT(50A1)
      IF (STRG(1).EQ.' ') GO TO 95
      IF (STRG(1).EQ.'!') RETURN
      IF ((STRG(1).EQ.'S').AND.(STRG(2).EQ.'T').AND.(STRG(3).EQ.'O')
     1.AND.(STRG(4).EQ.'P')) GO TO 95
      IF ((STRG(1).NE.'H').OR.(STRG(2).NE.'E').OR.(STRG(3).NE.'L')
     1.OR.(STRG(4).NE.'P')) GO TO 55
      WRITE(IDLG,54)
54    FORMAT('0TWO OPERATIONS ARE AVAILABLE TO COLLAPSE THE CHI-SQUARE:'/
     1'"C"-COMBINE AND "D"-COLLAPSE.  THEY MUST BE THE FIRST CHARACTER
     2 OF'/' THE LINE.  THIS MUST BE FOLLOWED BY THE VARIABLE NAME
     3 OR NUMBER'/' ENCLOSED IN PARENTHESIS.  THIS IS FOLLOWED BY
     4 A STRING OF NUMBERS'/' REFERENCING LEVELS OF THE VARIABLES.
     5 RANGES MAY BE SPECIFIED'/' BY SEPARATEING THE EXTREMES OF',
     5' THE RANGE WITH A "-".'/' A <CR>, ^Z, OR "STOP"',
     6' WILL EXIT FROM THE'/' COLLAPSING PORTION AND CONTINUE THE
     7 CHI SQUARE'/' EXAMPLE:'/' C(AGE)1-5'/' COMBINE THE LEVELS
     8 1 THRU 5 FOR THE VARIABLE AGE.'/' D(3)1,9,3'/' DELETE
     9 LEVELS 1,9,AND 3 FOR VARIABLE 3')
      GO TO 51
55    IF ((STRG(1).EQ.'C').OR.(STRG(1).EQ.'D')) GO TO 57
      WRITE(IDLG,56)
56    FORMAT(' INSTRUCTION MUST BE "C" OR "D"')
      GO TO 51
57    IF (STRG(2).EQ.'(') GO TO 60
58    WRITE(IDLG,59)
59    FORMAT(' VARIABLE NOT SPECIFIED, OR NOT IN PARENTHESIS')
      GO TO 51
C
C     DETERMINE VARIABLE TO BE USED IN COLLAPSING
C
60    DO 61 I=1,5
61    B(I)=' '
      NUM=0
      IF  ((STRG(3).LE.'9').AND.(STRG(3).GE.'0')) NUM=1
      K=1
      J=3
62    IF (STRG(J).EQ.')') GO TO 65
      IF (NUM.EQ.0) GO TO 63
      IF ((STRG(J).LE.'9').AND.(STRG(J).GE.'0')) GO TO 63
      GO TO 58
63    IF (K.GT.5) GO TO 64
      B(K)=STRG(J)
      K=K+1
64    J=J+1
      IF (J.LE.50) GO TO 62
65    IF (NUM.EQ.0) GO TO 165
      IF (B(5).NE.' ') GO TO 165
      DO 166 I=4,1,-1
166   B(I+1)=B(I)
      B(1)=' '
      GO TO 65
165   ENCODE(5,66,NNS) B
66    FORMAT(5A1)
      IF (NUM.EQ.0) GO TO 70
      DECODE(5,67,NNS) NS
67    FORMAT(I5)
72    IF ((IVB(1).EQ.NS).OR.(IVB(2).EQ.NS)) GO TO 75
68    WRITE(IDLG,69)
69    FORMAT(' VARIABLE SPECIFIED IS NOT BEING ANALYZED')
      GO TO 51
70    DO 71 NS=1,NV
      IF (NAMES(NS).EQ.NNS) GO TO 72
71    CONTINUE
      WRITE(IDLG,73) NNS
73    FORMAT(' VARIABLE "',A5,'" DOES NOT EXIST')
      GO TO 51
C
C     DETERMINE THE LEVELS TO BE ACTED UPON
C
75    K=1
      IF (NS.EQ.IVB(2)) K=2
      THRU=0
      M=1
175   J=J+1
      L=1
      DO 76 I=1,5
76    B(I)=' '
78    IF (STRG(J).EQ.' ') GO TO 79
      IF (STRG(J).EQ.',') GO TO 79
      IF (STRG(J).EQ.'-') GO TO 79
      IF ((STRG(J).LT.'0').OR.(STRG(J).GT.'9')) GO TO 181
      IF (L.GT.5) GO TO 77
      B(L)=STRG(J)
      L=L+1
77    J=J+1
      IF (J.LE.50) GO TO 78
79    IF (L.EQ.1) GO TO 83
80    IF (B(5).NE.' ') GO TO 82
      DO 81 I=4,1,-1
81    B(I+1)=B(I)
      B(1)=' '
      GO TO 80
181   WRITE(IDLG,182)
182   FORMAT(' LEVELS MAY ONLY BE SPECIFIED BY LEVEL NUMBERS')
      GO TO 51
82    ENCODE(5,66,NNS) B
      DECODE(5,67,NNS)NS
      IF ((NS.LE.0).OR.(NS.GT.LVL(K))) GO TO 183
      LTBM(M)=NS
      IF (THRU.EQ.1) LTBM(M)=-NS
      THRU=0
      M=M+1
      IF (STRG(J).EQ.'-') THRU=1
      IF (STRG(J).EQ.' ') GO TO 83
      GO TO 175
183   WRITE(IDLG,184) LVL(K)
184   FORMAT(' LEVEL SPECIFIED DOES NOT EXIST, MAXIMUM LEVEL IS',I4)
      GO TO 51
83    IF (M.GT.1) GO TO 85
      WRITE(IDLG,84)
84    FORMAT(' NO LEVELS TO BE ACTED UPON')
      GO TO 51
85    M=M-1
      KMIN=LTBM(1)
      DO 86 I=2,M
      L=LTBM(I)
      IF (L.LT.0) L=-L
      IF (KMIN.LE.L) GO TO 87
      KMIN=L
87    IF (LTBM(I).GT.0) GO TO 86
      L=-LTBM(I)
      IF (LTBM(I-1).GT.0) GO TO 89
      WRITE(IDLG,88)
88    FORMAT(' TWO - IN A ROW')
      GO TO 51
89    IF (LTBM(I-1).LT.L) GO TO 86
      LTBM(I)=-LTBM(I-1)
      LTBM(I-1)=L
86    CONTINUE
C
C     ACTUALLY GO THRU AND RECODE THOSE VARIABLES USED IN COLAPSE
C     PORTION
C
      J=1
190   L=LEVEL(J).AND.IRMAP
      IF(K.EQ.2) L=(LEVEL(J).AND.ILMAP)/IJUST
      DO 91 I=1,M
      IF (I.GE.M) GO TO 92
      IF (LTBM(I+1).GT.0) GO TO 92
      IF ((L.LT.LTBM(I)).OR.(L.GT.-LTBM(I+1))) GO TO 91
      GO TO 93
92    IF (L.NE.LTBM(I)) GO TO 91
93    IF(STRG(1).EQ.'D') GO TO 193
      IF(K.EQ.2) LEVEL(J)=(IRMAP.AND.LEVEL(J))+KMIN*IJUST
      IF(K.EQ.1) LEVEL(J)=(ILMAP.AND.LEVEL(J))+KMIN
      GO TO 90
193   IF(J.EQ.NK) GO TO 194
      DO 94 II=J,NK-1
      ISUB(II)=ISUB(II+1)
94    LEVEL(II)=LEVEL(II+1)
194   NK=NK-1
      IF(J.GT.NK) GO TO 51
      GO TO 190
91    CONTINUE
90    J=J+1
      IF(J.LE.NK) GO TO 190
      GO TO 51
C
C     RECODE LEVELS AFTER COLLAPSING
C
95    DO 96 I=1,NN
      L=I
      CALL CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,L)
      L=1
      KVAL=LEVEL(1).AND.IRMAP
      IF(I.EQ.2) KVAL=(LEVEL(1).AND.ILMAP)/IJUST
      DO 97 J=1,NK
      K=LEVEL(J).AND.IRMAP
      IF(I.EQ.2) K=(LEVEL(J).AND.ILMAP)/IJUST
      IF (K.EQ.KVAL) GO TO 98
      L=L+1
      KVAL=K
98    IF(I.EQ.2) LEVEL(J)=(IRMAP.AND.LEVEL(J)) +L*IJUST
      IF(I.EQ.1) LEVEL(J)=(ILMAP.AND.LEVEL(J))+L
97    CONTINUE
96    LVL(I)=L
      IF(NK.LT.1) RETURN
C
C     BEGIN OUTPUT PHASE
C
100   IF(IOUT.NE.21)WRITE(IOUT,5566) (HEDR(I),I=1,NSZ)
5566  FORMAT('1',70A1)
      IF(IOUT.EQ.21) CALL PRNTHD
      WRITE(IOUT,103)
103   FORMAT('0',20X,'***** CHI SQUARE *****')
      LINES=4
      IF (CONTG.EQ.1) GO TO 101
      CALL CHITAB(NK,LEVEL,ISUB,MV,MC,DATA,NAMES,IVB,2,NN,LINES)
101   IF (NN.GT.1) GO TO 105
110   A=LVL(1)
      IF(A.EQ.0) GO TO 106
      EXP=NK/A
      CALL CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,5)
      CHI=0
      KVAL=LEVEL(1)
      O=1
      DO 102 J=2,NK
      IF (KVAL.EQ.LEVEL(J)) GO TO 102
      CHI=CHI+(O-EXP)**2/EXP
      O=0
      KVAL=LEVEL(J)
102   O=O+1
      CHI=CHI+(O-EXP)**2/EXP
      NF=LVL(1)-1
      IF (NF.LT.1) GO TO 106
      IF(IOUT.NE.21) GO TO 111
      LINES=LINES+2
      IF(LINES.LE.(LINPP-1)) GO TO 111
      CALL PRNTHD
      LINES=4
111   WRITE(IOUT,104) CHI,NF
104   FORMAT('0CHI SQUARE =',G,' WITH',I4,' DEGREES OF FREEDOM')
      GO TO 123
C
C     CHECK TO SEE IF IT MAY BE ANALY. AS A ONE SAMPLE CHI SQ.
C
106   IF(IOUT.NE.21) GO TO 112
      LINES=LINES+2
      IF(LINES.LE.LINPP) GO TO 112
      CALL PRNTHD
      LINES=4
112   WRITE(IOUT, 107)
107   FORMAT('0CONTINGENCY TABLE DEGENERATED TO 1 CELL - NO CHI
     1 SQUARE PERFORMED')
      GO TO 15
108   IF (LVL(2).GT.LVL(1)) LVL(1)=LVL(2)
      LVL(2)=0
      WRITE(IDLG,109)
109   FORMAT(' CHI SQUARE ANALYSIS PERFORMED FOR 1 SAMPLE CASE')
      GO TO 110
105   IF(LVL(1).LT.2) GO TO 108
      IF(LVL(2).LT.2) GO TO 108
C
C     TWO SAMPLE CHI SQUARE
C
      NF=(LVL(1)-1)*(LVL(2)-1)
      CALL CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,5)
C
C     BREAK IT INTO TWO GROUPS FOR CALCULATION OF CHI SQUARE
C
      DO 115 I=1,NK
      ISUB(I)=(LEVEL(I).AND.ILMAP)/IJUST
115   LEVEL(I)=LEVEL(I).AND.IRMAP
C
C     CALCULATION FOR CHI SQUARE
C
      CHI=0
      K=1
      DO 116 I=1,LVL(2)
      RO=0
117   IF (ISUB(K).NE.I) GO TO 120
      RO=RO+1
      K=K+1
      IF (K.LE.NK) GO TO 117
120   DO 118 J=1,LVL(1)
      CO=0
      OBS=0
      DO 119 L=1,NK
      IF (LEVEL(L).NE.J) GO TO 119
      CO=CO+1
      IF (ISUB(L).EQ.I) OBS=OBS+1
119   CONTINUE
      EXP=(CO*RO)/NK
118   CHI=CHI+(OBS-EXP)**2/EXP
116   CONTINUE
      IF(IOUT.NE.21) GO TO 113
      LINES=LINES+2
      IF(LINES.LE.(LINPP-1)) GO TO 113
      CALL PRNTHD
      LINES=4
113   WRITE(IOUT,104) CHI,NF
      IF (NF.NE.1) GO TO 123
C
C     YATES CORRECTION FACTOR AUTOMATICALLY PERFORMED IF DEGREES OF
C     FREEDOM IS EQUAL TO ONE
C
      DO 125 I=1,2
      DO 125 J=1,2
125   MCC(I,J)=0
      DO 121 I=1,NK
121   MCC(ISUB(I),LEVEL(I))=MCC(ISUB(I),LEVEL(I))+1.
      CHI=MCC(1,1)*MCC(2,2)-MCC(1,2)*MCC(2,1)
      IF (CHI.LT.0) CHI=-CHI
      CHI=(CHI-NK/2.)**2
       CHI=(CHI*NK)/((MCC(1,1)+MCC(1,2))*(MCC(1,1)+MCC(2,1))*(MCC(2,2)
     1+MCC(1,2))*(MCC(2,2)+MCC(2,1)))
      IF(IOUT.NE.21) GO TO 128
      LINES=LINES+2
      IF(LINES.LE.(LINPP-1)) GO TO 128
      CALL PRNTHD
      LINES=4
128   WRITE(IOUT,122) CHI
122   FORMAT('0CORRECTED CHI SQUARE =',G)
123   CHI=CHI/NF
      PROB=FISHER(NF,1000,CHI)
      IF(IOUT.NE.21) GO TO 129
      LINES=LINES+1
      IF(LINES.LE.LINPP) GO TO 129
      CALL PRNTHD
      LINES=3
129   WRITE(IOUT,124) PROB
124   FORMAT(' HAVING A PROBABILITY OF ',F5.2)
      IF((NF.NE.1).OR.(FISHR.EQ.0).OR.(LVL(2).EQ.0)) GO TO 137
C
C     SECTION TO CALCULATE FISHERS EXACT PROBABILITY
C
      DO 138 I=1,2
      MCC(I,3)=0
138   MCC(3,I)=0
      DO 141 I=1,2
      MCC(I,3)=MCC(I,1)+MCC(I,2)
141   MCC(3,I)=MCC(1,I)+MCC(2,I)
      MCC(3,3)=MCC(1,3)+MCC(2,3)
      CALL CONP(MCC,3,3,PT,PS,PC)
      IF(IOUT.NE.21) GO TO 142
      LINES=LINES+4
      IF(LINES.LE.LINPP) GO TO 142
      CALL PRNTHD
      LINES=6
142   WRITE(IOUT,127) PT,PS
127   FORMAT('0FISHER''S EXACT PROBABILITY FOR OBTAINING THE ',
     1'GIVEN TABLE =',F8.5/23X,'FOR OBTAINING A TABLE AS ',
     2'PROBABLE, OR'/23X,'LESS PROBABLE THAN THE GIVEN',
     3' TABLE =',F8.5)
C
C
137   IF(NN.EQ.1) GO TO 15
      CONCOF=SQRT(CHI/(NK+CHI))
      IF(IOUT.NE.21) GO TO 143
      LINES=LINES+2
      IF(LINES.LE.LINPP) GO TO 143
      CALL PRNTHD
      LINES=4
143   WRITE(IOUT,126) CONCOF
126   FORMAT('0CONTINGENCY COEFFICIENT =',G)
      GO TO 15
      END
C                                           *** STAT PACK ***
C     SUBROUTINE CALLED BY CHI TO WRITE OUT CONTINGENCY TABLE
C     CALLING SEQUENCE: CALL CHITAB(NL,LEVEL,ISUB,MV,MC,DATA,NAMES,IX,ITM,NUMV,LINES)
C     WHERE NL - IS THE NUMBER OF OBS. IN CONTG TABLE
C           LEVEL - IS A VECTOR DIMENSIONED AT LEAST NL
C           ISUB - IS A VECTOR DIMENSIONED AT LEAST NL
C           MV - MAX. NUMBER OF VARIABLES
C           MC - MAX. NUMBER OF OBSERVATIONS
C           DATA - MATRIX FOR DATA, DIMENSIONED FOR MAXIMUM
C           NAMES - A VECTOR CONTAINING VARIABLE NAMES
C           ITM - SWITCH TO SHOW HOW THINGS SHOULD BE OUTPUT
C           NUMV - ONE OR TWO WAY CHISQ?
C           LINES - NUMBER OF LINES PRESENTLY USED ON THE LPT PAGE
C
C      SIMPLY TAKES DATA REFERENCED BY ISUB, LEVELS DETERMINED
C      BY LEVEL AND CREATES A CONTINGENCY TABLE
C
      SUBROUTINE CHITAB (NL,LEVEL,ISUB,MV,MC,DATA,NAMES,IX,ITM,NUMV,
     1LINES)
      COMMON/DEV/ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      DIMENSION LEVEL(1), ISUB(1), DATA(MC,MV), NAMES(1), XOUT(20),
     1IX(2), ICT(20), ITOT(20),IHDR(20)
      ILMAP="007777700000
      IRMAP="000000077777
      IJUST="000000100000
      ILL=10
      IF(ITM.EQ.1) GO TO 8
      IF(IOUT.EQ.21) ILL=20
8     DOTS='.....'
      ITAB0=IOUT
      IF(ITM.EQ.1) ITAB0=IDLG
      DO 10  KKK=1,NUMV
      KKKS=KKK
      CALL CHISOR(NL,LEVEL,ISUB,IX,DATA,MV,MC,KKKS)
      IF(ITAB0.NE.21) GO TO 40
      LINES=LINES+2
      IF(LINES.LE.(LINPP-10)) GO TO 40
      CALL PRNTHD
      LINES=4
40    IF((KKK.EQ.2).AND.(NUMV.NE.1)) WRITE(ITAB0,1) NAMES(IX(KKK))
1     FORMAT ('0LEVELS FOR VERTICAL VARIABLE: ',A5)
      IF((KKK.EQ.1).AND.(NUMV.NE.1)) WRITE(ITAB0,7) NAMES(IX(KKK))
7     FORMAT ('0LEVELS FOR HORIZONTAL VARIABLE: ',A5)
      IF(NUMV.EQ.1) WRITE(ITAB0,30)NAMES(IX(1))
30    FORMAT('0LEVELS FOR VARIABLE: ',A5)
      WRITE(ITAB0,2)
2     FORMAT ('0LEVEL',2X,'VALUES COMPRISING LEVEL')
      LEVL=1
      VALUE=DATA(ISUB(1),IX(KKK))
      IFRQ=1
      K=1
      XOUT(K)=VALUE
      K=K+1
      ISW=0
      DO 3 I=2,NL
      LH=LEVEL(I).AND.IRMAP
      IF(KKK.EQ.2) LH=(LEVEL(I).AND.ILMAP)/IJUST
      IF(LEVL.NE.LH) GO TO 4
      IF(VALUE.EQ.DATA(ISUB(I),IX(KKK))) GO TO 3
      XOUT(K)=','
      K=K+1
      IF(K.GT.ILL) GO TO 4
      VALUE=DATA(ISUB(I),IX(KKK)) 
      XOUT(K)=VALUE
      K=K+1
      GO TO 3
4     IF(ITAB0.NE.21) GO TO 41
      LINES=LINES+1
      IF(ISW.EQ.0) LINES=LINES+1
      IF(LINES.LE.(LINPP-1)) GO TO 41
      CALL PRNTHD
      LINES=3
      IF(ISW.EQ.0) LINES=LINES+1
41    IF(ISW.EQ.0) WRITE(ITAB0,5) LEVL,(XOUT(L),L=1,K-1)
5     FORMAT ('0',I3,4X,10(G10.4,A2))
      IF(ISW.EQ.1) WRITE (ITAB0,6)(XOUT(L),L=1,K-1)
6     FORMAT (8X,10(G10.4,A2))
      ISW=1
      K=1
      VALUE=DATA(ISUB(I),IX(KKK))
      XOUT(K)=VALUE
      K=K+1
      IF(LEVL.EQ.LH) GO TO 3
      IF(NUMV.EQ.1) WRITE(ITAB0,31) IFRQ
31    FORMAT(9X,'FREQUENCY =',I4)
      IFRQ=0
      ISW=0
      LEVL=LEVL+1
      GO TO 3
3     IFRQ=IFRQ+1
      IF(ITAB0.NE.21) GO TO 43
      LINES=LINES+1
      IF(ISW.EQ.0) LINES=LINES+1
      IF(LINES.LE.(LINPP-1)) GO TO 43
      CALL PRNTHD
      LINES=3
      IF(ISW.EQ.0) LINES=LINES+1
43    IF(ISW.EQ.0) WRITE(ITAB0,5) LEVL,(XOUT(L),L=1,K-1)
      IF(ISW.EQ.1) WRITE(ITAB0,6)(XOUT(L),L=1,K-1)
      IF(NUMV.EQ.1) WRITE(ITAB0,31)IFRQ
      IF(KKK.EQ.1) MAX=LEVL+1
10    CONTINUE
      IF(NUMV.EQ.1) RETURN
      LL=MAX/ILL
      IF((MAX-ILL*LL).NE.0) LL=LL+1
      DO 11  I=1,LL
      IADD=(I-1)*ILL
      IEND=ILL
      IF((IEND+IADD).GT.MAX) IEND=MAX-IADD
      DO 12  J=1,IEND
      L=IADD+J
      ENCODE(5,20,IHDR(J)) L
20    FORMAT (I5)
      ITOT(J)=0
12    CONTINUE
      IF((IADD+IEND).EQ.MAX) IHDR(IEND)='TOTAL'
      IF(ITAB0.NE.21) GO TO 45
      LINES=LINES+4
      IF(LINES.LE.(LINPP-10)) GO TO 45
      CALL PRNTHD
      LINES=6
45    WRITE (ITAB0,13) NAMES(IX(1)), (IHDR(J),J=1,IEND)
      WRITE (ITAB0,14) NAMES(IX(2)),(DOTS,DOTS,J=1,IEND)
13    FORMAT ('0',30X,A5/' LEVEL',2X,20(1X,A5))
14    FORMAT (1X,A5,1X,'.',40A3)
      LEVL=1
      DO 15  J=1,IEND
15    ICT(J)=0
      K=1
23    LH=(LEVEL(K).AND.ILMAP)/IJUST
      IRH=LEVEL(K).AND.IRMAP
      IF(LEVL.NE.LH) GO TO 17
      IF((IEND+IADD).EQ.MAX) ICT(IEND)=ICT(IEND)+1
      IF((IRH.LT.IADD+1).OR.(IRH.GT.IADD+IEND)) GO TO 16
      ICT(IRH-IADD)=ICT(IRH-IADD)+1
      GO TO 16
17    IF(ITAB0.NE.21) GO TO 46
      LINES=LINES+1
      IF(LINES.LE.LINPP) GO TO 46
      CALL PRNTHD
      WRITE (ITAB0,13) NAMES(IX(1)), (IHDR(J),J=1,IEND)
      WRITE (ITAB0,14) NAMES(IX(2)),(DOTS,DOTS,J=1,IEND)
      LINES=7
46    WRITE (ITAB0,18) LEVL,(ICT(J),J=1,IEND)
18    FORMAT (1X,I3,3X,'.',20(1X,I5))
      DO 19  J=1,IEND
      ITOT(J)=ITOT(J)+ICT(J)
19    ICT(J)=0
      LEVL=LEVL+1
      GO TO 23
16    K=K+1
      IF(K.LE.NL) GO TO 23
      DO 22 J=1,IEND
22    ITOT(J)=ITOT(J)+ICT(J)
      IF(ITAB0.NE.21) GO TO 47
      LINES=LINES+2
      IF(LINES.LE.LINPP) GO TO 47
      CALL PRNTHD
      WRITE (ITAB0,13) NAMES(IX(1)), (IHDR(J),J=1,IEND)
      WRITE (ITAB0,14) NAMES(IX(2)),(DOTS,DOTS,J=1,IEND)
      LINES=8
47    WRITE(ITAB0,18) LEVL,(ICT(J),J=1,IEND)
      WRITE (ITAB0,21)(ITOT(J),J=1,IEND)
21    FORMAT (1X,'TOTAL',1X,'.',20(1X,I5))
11    CONTINUE
      RETURN
      END
C                                                      *** STAT PACK ***
C      SUBROUTINE PART OF CHI USED TO SORT
C      CALLING SEQUENCE: CALL CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,ITYP)
C     WHRER NK - IS THE NUMBER OF OBSERVATIONS TO BE SORTED
C           LEVEL IS A VECTOR DIMENSIONED AT LEAST NK
C           ISUB IS A VECTOR DIMENSIONED AT LEAST NK
C           IBV - IS A VECTOR CONTAINING THE VARIABLE NUMBERS TO  BE USED
C           DATA - IS MATRIX CONTAINING DATA
C           MV - IS MAXIMUM NUMBER OF VARIABLES
C           MC - IS MAXIMUM NUMBER OF OBSERVATIONS
C           ITYP - TELLS IN WHICH WAY THEY ARE TO BE SORTED
C
C     SORTS NUMBERS BY SEVERAL DIFFERENT SEQUENCES
C
      SUBROUTINE CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,ITYP)
      DIMENSION DATA(MC,MV),LEVEL(1),ISUB(1),IVB(2)
      DIMENSION IU(16),IL(16)
      ILMAP="007777700000
      IRMAP="000000077777
      IBOTH="007777777777
      IF((ITYP.NE.3).AND.(ITYP.NE.4)) GO TO 1
      IFIN=0
      IS=ITYP-2
      GO TO 2
1     IF(ITYP.EQ.1) IFIN=IRMAP
      IF(ITYP.EQ.2) IFIN=ILMAP
      IF(ITYP.EQ.5) IFIN=IBOTH
      IS=ITYP
      IF(IS.GT.2) IS=1
2     IS=IVB(IS)
      M=1
      II=1
      J=NK
71    IF(II.GE.J) GO TO 78
72    K=II
      IJ=(J+II)/2
      T=DATA(ISUB(IJ),IS)
      IT=IFIN.AND.LEVEL(IJ)
      KT=LEVEL(II).AND.IFIN
      IF(KT.LT.IT) GO TO 73
      IF((KT.EQ.IT).AND.(DATA(ISUB(II),IS).LE.T)) GO TO 73
      ISAV=LEVEL(IJ)
      LEVEL(IJ)=LEVEL(II)
      LEVEL(II)=ISAV
      ISAV=ISUB(IJ)
      ISUB(IJ)=ISUB(II)
      ISUB(II)=ISAV
      IT=KT
      T=DATA(ISUB(IJ),IS)
73    LL=J
      KT=LEVEL(J).AND.IFIN
      IF(KT.GT.IT) GO TO 75
      IF((KT.EQ.IT).AND.(DATA(ISUB(J),IS).GE.T)) GO TO 75
      ISAV=LEVEL(IJ)
      LEVEL(IJ)=LEVEL(J)
      LEVEL(J)=ISAV
      ISAV=ISUB(IJ)
      ISUB(IJ)=ISUB(J)
      ISUB(J)=ISAV
      IT=KT
      T=DATA(ISUB(IJ),IS)
      KT=LEVEL(II).AND.IFIN
      IF(KT.LT.IT) GO TO 75
      IF((KT.EQ.IT).AND.(DATA(ISUB(II),IS).LE.T)) GO TO 75
      ISAV=LEVEL(IJ)
      LEVEL(IJ)=LEVEL(II)
      LEVEL(II)=ISAV
      ISAV=ISUB(IJ)
      ISUB(IJ)=ISUB(II)
      ISUB(II)=ISAV
      IT=KT
      T=DATA(ISUB(IJ),IS)
      GO TO 75
74    ISAV=LEVEL(LL)
      LEVEL(LL)=LEVEL(K)
      LEVEL(K)=ISAV
      ISAV=ISUB(LL)
      ISUB(LL)=ISUB(K)
      ISUB(K)=ISAV
75    LL=LL-1
      KT=LEVEL(LL).AND.IFIN
      IF(KT.GT.IT) GO TO 75
      IF((KT.EQ.IT).AND.(DATA(ISUB(LL),IS).GT.T)) GO TO 75
76    K=K+1
      KT=LEVEL(K).AND.IFIN
      IF(KT.LT.IT) GO TO 76
      IF((KT.EQ.IT).AND.(DATA(ISUB(K),IS).LT.T)) GO TO 76
      IF(K.LE.LL) GO TO 74
      IF((LL-II).LE.(J-K)) GO TO 77
      IL(M)=II
      IU(M)=LL
      II=K
      M=M+1
      GO TO 79
77    IL(M)=K
      IU(M)=J
      J=LL
      M=M+1
      GO TO 79
78    M=M-1
      IF(M.EQ.0) GO TO 90
      II=IL(M)
      J=IU(M)
79    IF((J-II).GE.11) GO TO 72
      IF(II.EQ.1) GO TO 71
      II=II-1
80    II=II+1
      IF(II.EQ.J) GO TO 78
      T=DATA(ISUB(II+1),IS)
      IT=LEVEL(II+1).AND.IFIN
      KT=LEVEL(II).AND.IFIN
      IF(KT.LT.IT) GO TO 80
      IF((KT.EQ.IT).AND.(DATA(ISUB(II),IS).LE.T)) GO TO 80
      K=II
      ISAV=LEVEL(II+1)
      ISAV2=ISUB(II+1)
81    LEVEL(K+1)=LEVEL(K)
      ISUB(K+1)=ISUB(K)
      K=K-1
      KT=LEVEL(K).AND.IFIN
      IF(IT.LT.KT) GO TO 81
      IF((IT.EQ.KT).AND.(T.LT.DATA(ISUB(K),IS))) GO TO 81
      LEVEL(K+1)=ISAV
      ISUB(K+1)=ISAV2
      GO TO 80
90    RETURN
      END
C                                  *** STAT PACK ***
C     SUBROUTINE IS PART OF STAT PACK CHISQUARE USED TO CALCULATE
C     THE EXACT PROBABILITY OF A 2X2 TABLE OR THE PROBABILITY OF
C     HAVING THAT TABLE OR A TABLE LESS PROBABLE THAN IT.
C     ROUTINE FROM COMMUNICATIONS OF ACM NOVEMBER 1972
C      WHERE MATRIX - IS A 3X3 TABLE FOR WHICH THE PROB IS FOUND
C            NR - IS THE NUMBER OF ROWS (HERE ALWAYS 2)
C            NC - IS THE NUMBER OF COLUMNS (HERE ALWAYS 2)
C            PT - PROBABILITY OF GIVEN TABLE
C            PS - PROBABILITY OF TABLE AS PROBABLE AS GIVEN TABLE
C            PC - NOT USED IN STP, BUT IT IS THE PROBABILITY OF
C                 OBTAINING SOME OF THE TABLES POSSIBLE WITHIN
C                 THE CONSTRAINTS OF THE MARGINAL TOTALS
C
C
      SUBROUTINE CONP(MATRIX,NR,NC,PT,PS,PC)
      DIMENSION MATRIX(NR,NC)
      COMMON/IFLAG1/IFLAG
      INTEGER R,C,TEMP
      IFLAG=0
      R=NR-1
      C=NC-1
      QXLOG=-FACLOG(MATRIX(NR,NC))
      DO 10 I=1,R
10    QXLOG=QXLOG+FACLOG(MATRIX(I,NC))
      DO 20 J=1,C
20    QXLOG=QXLOG+FACLOG(MATRIX(NR,J))
      RXLOG=0.0
      DO 50 I=1,R
      DO 50 J=1,C
50    RXLOG=RXLOG+FACLOG(MATRIX(I,J))
      PT=10.0**(QXLOG-RXLOG)
      PS=0
      PC=0
      DO 100 I=2,R
      DO 100 J=2,C
100   MATRIX(I,J)=MIN0(MATRIX(I,NC),MATRIX(NR,J))
      GO TO 300
200   DO 220 I=2,R
      DO 220 J=2,C
      MATRIX(I,J)=MATRIX(I,J)-1
      IF(MATRIX(I,J).GE.0) GO TO 300
220   MATRIX(I,J)=MIN0(MATRIX(I,NC),MATRIX(NR,J))
      RETURN
300   DO 320 I=2,R
      TEMP=MATRIX(I,NC)
      DO 310 J=2,C
310    TEMP=TEMP-MATRIX(I,J)
      IF(TEMP.LT.0) GO TO 200
320   MATRIX(I,1)=TEMP
      DO 340 J=1,C
      TEMP=MATRIX(NR,J)
      DO 330 I=2,R
330   TEMP=TEMP-MATRIX(I,J)
      IF(TEMP.LT.0) GO TO 200
340   MATRIX(1,J)=TEMP
      RXLOG=0.0
      DO 350 I=1,R
      DO 350 J=1,C
350   RXLOG=RXLOG+FACLOG(MATRIX(I,J))
      PX=10.0**(QXLOG-RXLOG)
      PC=PC+PX
      IF((PT/PX).GT.0.9999) PS=PS+PX
      GO TO 200
      END
C                                            *** STAT PACK ***
C     FUNCTION IS PART OF STP TAKEN FOR COMMUNICATIONS OF ACM
C     NOVEMBER 1972. USED TOR FINDING LOG BASE 10 OF N FACTORIAL.
C     USES STIRLINGS APPROX. IF N.GT.100
C
C     FIRST TIME THRU IT CREATES A TABLE IT USES.  IFLAG IS THE 
C     INDICATOR USED TO ESTABLISH RATHER A TABLE NEED BE CREATED
C
      FUNCTION FACLOG(N)
      DIMENSION TABLE (101)
      COMMON /IFLAG1/IFLAG
      TPILOG=.3990899342
      ELOG=.4342944819
      IF(N.GT.100) GO TO 50
      IF(IFLAG.EQ.0) GOTO 100
10    FACLOG=TABLE(N+1)
      RETURN
50    X=FLOAT(N)
      FACLOG=(X+.5)*ALOG10(X)-X*ELOG+TPILOG+ELOG/(12.0*X)-ELOG
     1/(360.0*X**3)
      RETURN
C
C     CREATE TABLE TO BE USED FOR REST OF TIME
C
100   TABLE(1)=0.0
      DO 120 I=2,101
      X=FLOAT(I-1)
120   TABLE(I)=TABLE(I-1)+ALOG10(X)
      IFLAG=1
      GOTO 10
      END