Trailing-Edge
-
PDP-10 Archives
-
decuslib10-09
-
43,50466/stp11.stp
There is 1 other file named stp11.stp 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