Trailing-Edge
-
PDP-10 Archives
-
decuslib10-09
-
43,50466/stp7.stp
There is 1 other file named stp7.stp in the archive. Click here to see a list.
C *** STAT PACK ***
C SUBROUTINE FOR 2 WAY ANALYSIS OF VARIANCE.
C CALLING SEQUENCE: CALL ANOV2(NV,NC,MV,MC,DATA,VMN,STD,NAMES)
C WHERE NV - IS THE NUMBER OF COLUMNS ACTUALLY FILLED (VARIABLES)
C NC - IS THE NUMBER OF ROWS ACTUALLY FILLED (CASES)
C MV - IS THE MAXIMUM NUMBER OF COLUMNS.
C MC - IS THE MAXIMUM NUMBER OF ROWS.
C DATA - DATA MATRIX, DIMENSIONED FOR MAXIMUM.
C VMN - IS A VECTOR CONTAINING THE MEANS.
C STD - IS A VECTOR CONTAINING THE STANDARD DEVIATIONS.
C NAMES - IS A VECTOR CONTAINING VARIABLE NAMES
C
C SUBROUTINE ALLOWS FOR EACH CELL TO BE AN INDIVIDUAL VARIABLE
C OR FOR A BREAKDOWN ON THE VALUE OF ANOTHER VARIABLE TO DETERMINE
C THE CELLS. THIS PARTICULAR ANOV IS USED SINCE IT DOES ALLOW
C EMPTY CELLS OR AN UNBALANCED DESIGN. ORIGINALLY WRITTEN
C AT W. M. U. BY EVA GAINES, UNDER DIRECTION OF MICHEAL STOLINE,
C IT WAS MODIFIED DURING THE SUMMER OF 71 BY SAM ANEMA,
C THE VERSION NOW IN STP IS A FURTHUR MODIFICATION OF THAT
C PROGRAM. ORIGINAL VERSION IS BASED ON BOOK BY BANCROFT
C "INTERMEDIATE STATISTICAL METHODS".
C
SUBROUTINE ANOV2(NV,NC,MV,MC,DATA,VMN,STD,NAMES)
DIMENSION DATA(MC,MV), VMN(1), STD(1), OPT(5), OPTD(5)
COMMON /DEV/ ICC, IDATA, IOUT, IDLG,IDSK
COMMON /PRNT/ LINPP,ICOPS,RUNPRG
COMMON/EXTRA/HEDR(70),NSZ
DIMENSION YS(20,20),Y1(20),Y2(20),S1(20),S2(20),B1(20)
DIMENSION T(20,20),C(20),AM(20),YT(20,20),IBK(2),BRNG(2,20,2)
DIMENSION NGRPS(2),IVA(120)
DIMENSION ICELLS(40),NAMES(1)
EQUIVALENCE (NGRPS(1),NGRPA),(NGRPS(2),NGRPB)
EQUIVALENCE (IBK(1),IBK1),(IBK(2),IBK2)
EQUIVALENCE (IVA,ICELLS),(IVA(41),BRNG)
ISQ=7
IF(IOUT.EQ.21) ISQ=15
ALL=0
5 IF(ICC.NE.2) WRITE (IDLG,1)
OPTD(1)=0
OPTD(2)=0
OPTD(3)=0
OPTD(4)=0
OPTD(5)=0
NGRPS(1)=20
NGRPS(2)=20
1 FORMAT ('0LIST OPTIONS SEPARATED BY COMMAS'/)
READ (ICC,3) OPT
2 FORMAT ('0THE 2 WAY ANALYSIS OF VARIANCE ASSUMES
1 CELLS TO '/' BE MADE UP OF INDIVIDUAL VARIABLES
2.'/' IT IS POSSIBLE HOWEVER TO BREAK
3 1 VARIABLE INTO'/' CELLS, BY MEANS OF
4 TWO OTHER VARIABLES.'/' IN THIS CASE, RANGES
5 FOR BREAK-DOWNS ARE ASSUMED'/' TO BE SUPPLIED
6 BY THE USER. WHEN THE PROGRAM'/' CALLS FOR
7 THE RANGES, THE USER MAY SPECIFY ANOTHER'/' OPTION
8 WHICH ALLOWS THE PROGRAM TO AUTOMATICALLY FORM'/' THE
9 BREAKDOWNS BY INDIVIDUAL VALUES OF THE')
632 FORMAT(' BREAKDOWN VARIABLE. IF MORE VALUES'/
2 ' EXIST THAN GROUPS, THE PROGRAM AUTOMATICALLY'/' FORMS
3 THE RANGES FOR BREAKDOWN'/ ' OPTIONS ARE:'/
5 ' "BREAK" - FORM GROUPINGS BASED ON VALUES OF BREAKDOWN
6 VARIABLES.'/' "GROUP" - SPECIFY MAXIMUM NUMBER OF GROUPS
7 (PRESET TO 20) MAXIMUM OF 20.'/' "AUTO" - AUTOMATIC BREAKDOWN
8 (SPECIFY THIS ONLY WHEN ASKED FOR RANGES)'/' "DISCR" - BOTH',
8' GROUPS ARE AUTO (SPECIFY NOW)'/' "HEADR" -',
9' OMIT SIZE, MEANS, AND STD. DEV. REPORT'/' "RANGE" - TYPE OUT',
1' RANGES FOR BREAKDOWNS IF AUTOMATICALLY FORMED'/
2'0IF NO OPTIONS ARE DESIRED TYPE A RETURN')
3 FORMAT (5(A5,1X))
IF(OPT(1).EQ.'!') RETURN
DO 4 I=1,5
IF (OPT(I).EQ.' ') GO TO 10
IF (OPT(I).NE.'HELP') GO TO 6
WRITE (IDLG,2)
WRITE(IDLG,632)
GO TO 5
6 IF (OPT(I).NE.'BREAK') GO TO 7
OPTD(1)=1
GO TO 4
7 IF (OPT(I).NE.'GROUP') GO TO 90
OPTD(2)=1
GO TO 4
90 IF(OPT(I).NE.'HEADR') GO TO 308
OPTD(3)=1
GO TO 4
308 IF(OPT(I).NE.'DISCR') GO TO 508
OPTD(4)=1
GO TO 4
508 IF(OPT(I).NE.'RANGE') GO TO 509
OPTD(5)=1
GO TO 4
509 IF(OPT(I).NE.'AUTO') GO TO 512
510 WRITE(IDLG,511)
511 FORMAT(' "AUTO" ONLY SPECIFIED WHEN RANGES ARE CALLED FOR')
GO TO 5
512 IF(OPT(I).EQ.'AUTO,') GO TO510
8 WRITE (IDLG,9) OPT(I)
9 FORMAT ('0OPTION ',A5,' DOES NOT EXIST TRY AGAIN')
GO TO 5
4 CONTINUE
10 IF (OPTD(1).EQ.1) GO TO 30
C
C CELLS INDIVIDUAL VARIABLES.
C
11 IF(ICC.NE.2) WRITE (IDLG,12)
12 FORMAT ('0HOW MANY CELLS IN FACTOR 1? ',$)
READ (ICC,13) NGRPA
13 FORMAT(I)
IF (NGRPA.LE.20) GO TO 15
WRITE (IDLG,14)
14 FORMAT ('0MAXIMUM OF 20')
GO TO 11
15 IF(ICC.NE.2) WRITE (IDLG,16)
16 FORMAT ('0HOW MANY CELLS IN FACTOR 2? ',$)
READ (ICC,13) NGRPB
IF (NGRPB.LE.20) GO TO 18
WRITE (IDLG,14)
GO TO 15
18 IF(ICC.NE.2) WRITE (IDLG,19)
19 FORMAT ('0TYPE IN EACH VARIABLE AFTER THE'/
1' CORRESPONDING LEVEL, "EMPTY"-INDICATES EMPTY CELL')
XT=NC
L=-1
DO 26 I=1,NGRPA
DO 26 J=1,NGRPB
22 IF(ICC.NE.2) WRITE (IDLG,21)I,J
21 FORMAT ('0CELL('I2,',',I2,')? ',$)
IRET=-99
CALL ALPHA(ICELL,1,K,IRET,IHELP,IERR,NAMES,NV)
IF(IRET.EQ.1) RETURN
IF(IHELP.EQ.1) GO TO 22
IF(IERR.EQ.1) GO TO 22
IF(ICELL.GE.0) GO TO 441
L=L-1
IF(L.GE.-NV) GO TO 28
445 WRITE(IDLG,27)
27 FORMAT(' TOO MANY CELLS FOR THE NUMBER OF VARIABLES')
GO TO 11
28 ICELL=L
GO TO 24
441 K=(I-1)*NGRPB+J-1
IF(K.LT.1) GO TO 24
DO 442 II=1,K
IF(IVA(II).LE.0) GO TO 442
IF(IVA(II).NE.ICELL) GO TO 442
WRITE(IDLG,443)
443 FORMAT(' THIS VARIABLE USED PREVIOUSLY')
GO TO 22
442 CONTINUE
24 K=(I-1)*NGRPB+J
26 IVA(K)=ICELL
NGPALL=NGRPA*NGRPB
K=0
DO 421 I=1,NGPALL
IF(IVA(I).EQ.0) GOTO 421
K=K+1
421 CONTINUE
IF(K.GT.NV) GO TO 445
GOTO 408
400 J=NGPALL
404 IF(IVA(J).GT.0) GO TO 405
IVA(J)=IVA(J)-1
IF(IVA(J).GE.-NV) GO TO 406
405 J=J-1
IF(J.GE.1) GO TO 404
RETURN
406 K=IVA(J)
IF(J.EQ.NGPALL) GO TO 408
DO 407 I=J+1,NGPALL
IF(IVA(I).GT.0) GO TO 407
K=K-1
IF(K.LT.-NV) GO TO 405
IVA(I)=K
407 CONTINUE
408 DO 403 I=1,NGPALL-1
DO 403 J=I+1,NGPALL
IF((IVA(I).EQ.0).OR.(IVA(J).EQ.0)) GO TO 403
IF(IVA(I).LT.0) GO TO 409
IF(IVA(J).GT.0) GOTO 410
GO TO 411
409 IF(IVA(J).LT.0) GO TO 410
411 IF(IVA(J).EQ.-IVA(I)) GO TO 400
GO TO 403
410 IF(IVA(J).EQ.IVA(I)) GO TO 400
403 CONTINUE
DO 88 I=1,20
AM(I)=0
Y1(I)=0
S1(I)=0
C(I)=0
B1(I)=0
S2(I)=0
Y2(I)=0
DO 88 J=1,20
T(I,J)=0
YT(I,J)=0
88 YS(I,J)=0
RAB=0
SST=0
SS=0
G=0
SSC=0
SSAB=0
AN=0
SSBA=0
TW=0
TC=0
DO 20 I=1,NGRPA
DO 25 J=1,NGRPB
K=(I-1)*NGRPB+J
ICELL=IVA(K)
IF(ICELL.EQ.0) GO TO 20
IF(ICELL.LT.0) ICELL=-ICELL
23 TC=TC+1.
TW=TW+XT-1.
YT(I,J)=VMN(ICELL)*XT
T(I,J)=XT
YS(I,J)=((STD(ICELL)**2)*(XT*(XT-1.))+YT(I,J)**2)/XT
Y1(I)=Y1(I)+YT(I,J)
S1(I)=S1(I)+T(I,J)
SSC=SSC+(YT(I,J)*YT(I,J))/T(I,J)
SST=SST+YS(I,J)
25 CONTINUE
IF(S1(I).EQ.0) GO TO 20
SSAB=SSAB+(Y1(I)**2)/S1(I)
SS=SS+S1(I)
20 G=G+Y1(I)
GO TO 70
C
C CELLS ON INDIV. VAR. BY BREAK-DOWN.
C
30 IF (OPTD(2).NE.1) GO TO 36
31 IF(ICC.NE.2) WRITE (IDLG,32)
32 FORMAT ('0MAXIMUM NUMBER OF GROUPS FOR LEVEL 1? ',$)
READ (ICC,13) NGRPA
IF ((NGRPA.GT.1).AND.(NGRPA.LE.20)) GO TO 34
WRITE (IDLG,33)
33 FORMAT ('0NO. OF GROUPS SPECIFIED WAS EITHER LESS THAN 1'/
1' OR GREATTER THAN THE MAXIMUM (20)')
GO TO 31
34 IF(ICC.NE.2) WRITE (IDLG,35)
35 FORMAT ('0MAXIMUM NUMBER OF GROUPS FOR LEVEL 2? ',$)
READ (ICC,13) NGRPB
IF ((NGRPB.GT.1).AND.(NGRPB.LE.20)) GO TO 36
WRITE (IDLG,33)
GO TO 34
36 IF(ICC.NE.2) WRITE (IDLG,38)
38 FORMAT ('0WHICH VARIABLES ARE TO BE ANALYSED? ',$)
CALL ALPHA(ICELLS,40,NZZ,IRET,IHELP,IERR,NAMES,NV)
IF(IRET.EQ.1) RETURN
IF(IHELP.EQ.1) GO TO 36
IF(IERR.EQ.1) GO TO 36
DO 302 I=1,NZZ
IF(ICELLS(I).NE.-1) GO TO 302
ALL=1
NZZ=NV
GO TO 39
302 CONTINUE
39 DO 40 I=1,2
42 IF(ICC.NE.2) WRITE (IDLG,41) I
41 FORMAT ('0WHICH VARIABLE IS BREAKDOWN VARIABLE',I2,'? ',$)
CALL ALPHA(ICELL,1,J,IRET,IHELP,IERR,NAMES,NV)
IF(IRET.EQ.1) RETURN
IF(IHELP.EQ.1) GO TO 42
IF(IERR.EQ.1) GO TO 42
IF(ICELL.GT.0)GO TO 40
WRITE(IDLG,450)
450 FORMAT(' *, ?, ALL MAY NOT BE USED HERE')
GO TO 42
40 IBK(I)=ICELL
DO 45 I=1,2
IF(OPTD(4).EQ.1) GO TO 50
44 IF(ICC.NE.2) WRITE (IDLG,46) NAMES(IBK(I))
46 FORMAT ('0LIST THE RANGES FOR BREAKDOWN VARIABLE: ',A5/)
DO 43 J=1,20
57 READ (ICC,47,END=37,ERR=37) HELP
47 FORMAT (A4,70X)
IF(HELP.EQ.' ') GO TO 37
IF(HELP.EQ.'!') RETURN
IF (HELP.NE.'HELP') GO TO 59
WRITE (IDLG,48)
48 FORMAT ('IF AUTOMATIC BREAK-DOWN BY GROUPS IS DESIRED TYPE "AUTO",
1'/' OTHERWISE ENTER THE RANGE FOR EACH BREAK-DOWN SMALLER FIRST,
2 1 PER LINE.'/'EXAMPLE:'/'17,35.2'/'CONTROL Z (^Z) WILL TERMINATE
3 ENTRY')
GO TO 44
59 IF (HELP.NE.'AUTO') GO TO 55
300 IF (J.EQ.1) GO TO 50
WRITE (IDLG,49)
49 FORMAT ('0AUTO CAN ONLY BE SPECIFIED FIRST TIME'/' CONFUSION
1 START OVER WITH SAME VARIABLE')
GO TO 44
55 REREAD 56, BRNG(I,J,1),BRNG(I,J,2)
56 FORMAT (2F)
IF (BRNG(I,J,1).LE.BRNG(I,J,2)) GO TO 43
WRITE (IDLG,58)
58 FORMAT (' SMALLER FIRST, PLEASE RETYPE'/)
GO TO 57
43 CONTINUE
J=J+1
37 NGRPS(I)=J-1
IF(NGRPS(I).LT.1) RETURN
GO TO 45
50 XMAX=-10.**11
XMIN=10.**11
N=0
L=IBK(I)
DO 51 K=1,NC
IF (XMAX.LT.DATA(K,L))XMAX=DATA(K,L)
IF (XMIN.GT.DATA(K,L))XMIN=DATA(K,L)
IF (N.GT.20) GO TO 51
IF (N.LT.1) GO TO 53
DO 52 M=1,N
IF (DATA(K,L).EQ.BRNG(I,M,1)) GO TO 51
52 CONTINUE
53 N=N+1
IF(N.GT.20) GO TO 51
BRNG (I,N,1)=DATA(K,L)
BRNG (I,N,2)=DATA(K,L)
51 CONTINUE
IF (N.LT.NGRPS(I)) NGRPS(I)=N
IF (N.EQ.NGRPS(I)) GO TO 89
XT=NGRPS(I)
RG=(XMAX-XMIN)/XT
BRNG(I,1,1)=XMIN
BRNG(I,1,2)=XMIN+RG
DO 54 K=2,NGRPS(I)
BRNG(I,K,1)=BRNG(I,K-1,2)
BRNG(I,K,2)=BRNG(I,K,1)+RG
54 CONTINUE
BRNG(I,NGRPS(I),2)=XMAX
89 IF(OPTD(5).NE.1) GO TO 420
IF(OPTD(4).EQ.1) WRITE(IDLG,301) NAMES(IBK(I))
301 FORMAT('0BREAKDOWN RANGES FOR VARIABLE: ',A5/)
C *******************************************************
C SORT BEFORE PUTTING OUT RANGES (BUBBLE SORT BECAUSE OF SIZE)
C
420 J=0
310 J=J+1
IF(J.EQ.NGRPS(I)) GO TO 312
TAB=BRNG(I,J+1,1)
TBA=BRNG(I,J+1,2)
IF(BRNG(I,J,1).LT.TAB) GO TO 310
K=J
311 BRNG(I,K+1,1)=BRNG(I,K,1)
BRNG(I,K+1,2)=BRNG(I,K,2)
K=K-1
IF(K.LT.1) GO TO 313
IF(TAB.LT.BRNG(I,K,1)) GO TO 311
313 BRNG(I,K+1,1)=TAB
BRNG(I,K+1,2)=TBA
GO TO 310
C ***********************************************************
312 IF(OPTD(5).NE.1) GO TO45
DO 71 J=1,NGRPS(I)
71 WRITE(IDLG,72)(BRNG(I,J,K),K=1,2)
72 FORMAT('+',G10.4,',',G10.4/)
GO TO 45
45 CONTINUE
C
C AT THIS POINT RANGES HAVE BEEN DET.
C BREAK DOWN VARIABLES
C
DO 2000 KCELL=1,NZZ
IF(ALL.EQ.1) GO TO 306
ICELL=ICELLS(KCELL)
GO TO 305
306 IF((KCELL.EQ.IBK1).OR.(KCELL.EQ.IBK2)) GO TO 2000
ICELL=KCELL
305 DO 388 I=1,20
AM(I)=0
Y1(I)=0
S1(I)=0
C(I)=0
B1(I)=0
S2(I)=0
Y2(I)=0
DO 388 J=1,20
T(I,J)=0
YT(I,J)=0
388 YS(I,J)=0
RAB=0
SST=0
SS=0
G=0
SSC=0
SSAB=0
AN=0
SSBA=0
TW=0
TC=0
DO 60 I=1,NC
C ----------------------------------------------------
C WHICH RANGE IN BREAKDOWN VARIABLE 1 DOES IT BELONG TO
DO 61 J=1,NGRPA
IF (DATA(I,IBK1).LT.BRNG(1,J,1)) GO TO 61
IF (DATA(I,IBK1).GT.BRNG(1,J,2)) GO TO 61
C -----------------------------------------------------------
C WHICH RANGE IN BREAKDOWN VARIABLE 2 DOES IT BELONG
62 DO 63 K=1,NGRPB
IF (DATA(I,IBK2).LT.BRNG(2,K,1)) GO TO 63
IF (DATA(I,IBK2).GT.BRNG(2,K,2)) GO TO 63
C --------------------------------------------------------
C ADD IT INTO THE CORRECT TOTALS
64 YT(J,K)=YT(J,K)+DATA(I,ICELL)
YS(J,K)=YS(J,K)+DATA(I,ICELL)**2
T(J,K)=T(J,K)+1.
GO TO 60
63 CONTINUE
GO TO 60
61 CONTINUE
60 CONTINUE
92 DO 65 I=1,NGRPA
DO 66 J=1,NGRPB
IF(T(I,J).EQ.0) GO TO 66
TC=TC+1
TW=TW+T(I,J)-1.
Y1(I)=Y1(I)+YT(I,J)
S1(I)=S1(I)+T(I,J)
SSC=SSC+YT(I,J)**2/T(I,J)
SST=SST+YS(I,J)
66 CONTINUE
IF(S1(I).EQ.0) GO TO 65
SSAB=SSAB+(Y1(I)*Y1(I))/S1(I)
SS=SS+S1(I)
65 G=G+Y1(I)
GO TO 70
C
C BEGIN ANALYSIS OF VARIANCE.
C
70 N1=NGRPA
N2=NGRPB
DO 93 J=1,N2
DO 93 I=1,N1
93 S2(J)=S2(J)+T(I,J)
SWTCH=0
DO 73 I=1,N1
IF(S1(I).NE.0) GO TO 73
IF(SWTCH.EQ.1) GO TO 75
WRITE(IDLG,74)
SWTCH=1
74 FORMAT('OTHE FOLLOWING LEVELS ARE EMPTY - WILL BE ELIMIN.')
75 WRITE(IDLG,76) I
76 FORMAT(' FACTOR 1 LEVEL',I3)
73 CONTINUE
DO 77 J=1,N2
IF(S2(J).NE.0) GO TO 77
IF(SWTCH.EQ.1) GO TO 78
WRITE(IDLG,74)
SWTCH=1
78 WRITE(IDLG,79) J
79 FORMAT(' FACTOR 2 LEVEL',I3)
77 CONTINUE
IF(SWTCH.EQ.0) GO TO 80
81 DO 82 I=1,NGRPA
IF(S1(I).NE.0) GO TO 82
NGRPA=NGRPA-1
IF(I.GT.NGRPA) GO TO 84
DO 83 K=I,NGRPA
DO 83 J=1,NGRPB
YT(K,J)=YT(K+1,J)
YS(K,J)=YS(K+1,J)
T(K,J)=T(K+1,J)
S1(K)=S1(K+1)
83 CONTINUE
GO TO 81
82 CONTINUE
84 DO 85 J=1,NGRPB
IF(S2(J).NE.0) GO TO 85
NGRPB=NGRPB-1
IF(J.GT.NGRPB) GO TO 87
DO 86 K=J,NGRPB
DO 86 I=1,NGRPA
YT(I,K)=YT(I,K+1)
YS(I,K)=YS(I,K+1)
T(I,K)=T(I,K+1)
S2(K)=S2(K+1)
86 CONTINUE
GO TO 84
85 CONTINUE
87 DO 91 I=1,20
Y1(I)=0
AM(I)=0
S1(I)=0
C(I)=0
B1(I)=0
S2(I)=0
91 Y2(I)=0
RAB=0
SST=0
SS=0
G=0
SSC=0
SSAB=0
AN=0
SSBA=0
TW=0
TC=0
GO TO 92
80 DO 100 J=1,N2
DO 101 I=1,N1
IF(T(I,J).EQ.0) GO TO 101
102 Y2(J)=Y2(J)+YT(I,J)
YS(I,J)=YS(I,J)-(YT(I,J)*YT(I,J))/T(I,J)
IF(T(I,J).EQ.1.) GO TO 101
106 YS(I,J)=YS(I,J)/(T(I,J)-1.)
IF(YS(I,J).LT.0) YS(I,J)=0.
YS(I,J)=SQRT(YS(I,J))
YT(I,J)=YT(I,J)/T(I,J)
101 CONTINUE
SSBA=SSBA+(Y2(J)*Y2(J))/S2(J)
100 CONTINUE
G1=(G*G)/SS
SST=SST-G1
SSC=SSC-G1
SSAB=SSAB-G1
SSBA=SSBA-G1
SSW=SST-SSC
NP=0
NP=NP+1
IF(TC.EQ.1.) GO TO 149
104 SMC=SSC/(TC-1.)
NP=NP+1
ITYPE=0
IF(TW.EQ.0) ITYPE=1
105 IF(ITYPE.EQ.1) GO TO 505
SMW=SSW/TW
NP=NP+1
IF(SMW.EQ.0) GO TO 149
103 F=SMC/SMW
505 IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(I),I=1,NSZ)
5566 FORMAT('1',70A1)
IF(IOUT.EQ.21) CALL PRNTHD
WRITE(IOUT,321)
321 FORMAT('0',15X,'*****2-WAY ANOVA*****')
LINES=4
IF(OPTD(1).NE.1) GO TO 323
WRITE(IOUT,320)NAMES(ICELL),NAMES(IBK1),NAMES(IBK2)
320 FORMAT(' ANALYSIS RUN ON VARIABLE: ',A5,' WITH CELLS DETERMINED'/
1' BY BREAKDOWNS ON VARIABLE: ',A5,' AND VARIABLE: ',A5)
LINES=LINES+2
323 IF(OPTD(3).EQ.1) GO TO 155
LINES=LINES+3+(N2+ISQ-1)/ISQ
324 WRITE (IOUT,107)
107 FORMAT('0FACTOR'/2X,'ONE',25X,'FACTOR TWO')
DO 322 K=1,N2,ISQ
NEND=K+ISQ-1
IF(NEND.GT.N2) NEND=N2
322 WRITE(IOUT,108)(I,I=K,NEND)
108 FORMAT(' LEVEL',5X,15(4X,I3,1X))
DO 116 I=1,N1
IF(IOUT.NE.21) GO TO 325
LINES=LINES+3*((N1+ISQ-1)/ISQ)+1
IF(LINES.LE.LINPP) GO TO 325
CALL PRNTHD
WRITE(IOUT,107)
DO 326 K=1,N2,ISQ
NEND=K+ISQ-1
IF(NEND.GT.N2) NEND=N2
326 WRITE(IOUT,108)(IK,IK=K,NEND)
LINES=5+((N2+ISQ-1)/ISQ)*4+1
325 DO 116 K=1,N2,ISQ
NEND=K+ISQ-1
IF(NEND.GT.N2) NEND=N2
IF(K.EQ.1) GO TO 112
WRITE(IOUT,110)(T(I,J),J=K,NEND)
110 FORMAT('0',5X,'N',4X,15F8.2)
GO TO 111
112 WRITE (IOUT,113)I,(T(I,J),J=K,NEND)
113 FORMAT('0',I3,2X,'N',4X,15F8.2)
111 WRITE (IOUT,114)(YT(I,J),J=K,NEND)
114 FORMAT(6X,'MEAN',1X,15F8.2)
115 WRITE (IOUT,154)(YS(I,J),J=K,NEND)
154 FORMAT(6X,'STDV',1X,15F8.2)
116 CONTINUE
IF(ITYPE.EQ.1) GO TO 335
155 IF(IOUT.NE.21) GO TO 327
LINES=LINES+9
IF(LINES.LE.LINPP) GO TO 327
CALL PRNTHD
LINES=11
327 WRITE (IOUT,117)
117 FORMAT(//24X,'PRELIMINARY ANOVA')
WRITE (IOUT,118)
118 FORMAT(6X,'SOURCE',11X,'DF',8X,'SS',8X,'MS',8X,'F',4X,'PROB')
KK=TC-1.
NDG=TW
PROB=FISHER(KK,NDG,F)
WRITE (IOUT,119)KK,SSC,SMC,F,PROB
119 FORMAT(6X'CELLS',9X,I5,3F10.2,2X,F6.4)
KK=N1-1
I=1
J=2
WRITE (IOUT,120)I,J,KK,SSAB
120 FORMAT(4X,I1,' IGNORING',I2,4X,I5,F10.2)
KK=N2-1
WRITE (IOUT,120)J,I,KK,SSBA
KK=TW
WRITE (IOUT,121)KK,SSW,SMW
121 FORMAT(6X,'WITHIN',8X,I5,2F10.2)
KK=SS-1.
WRITE (IOUT,122)KK,SST
122 FORMAT(6X,'TOTAL',9X,I5,F10.2)
335 DO 125 I=1,N2
DO 125 J=1,N2
YT(I,J)=0.
IF(I-J)124,123,124
123 YS(I,J)=1.
GO TO 125
124 YS(I,J)=0.
125 CONTINUE
DO 128 I=1,N2
DO 128 J=1,N2
DO 128 K=1,N1
IF(S1(K))128,126,128
126 TYPE 127,K,S1(K)
127 FORMAT(3HS1(,I3,2H)=,F6.0)
128 YT(I,J)=YT(I,J)+(T(K,I)*T(K,J))/S1(K)
DO 129 I=1,N2
129 YT(I,I)=YT(I,I)-S2(I)
DO 133 J=1,N2
130 DO 132 K=1,N1
IF(S1(K))132,131,132
131 TYPE 127,K,S1(K)
132 C(J)=C(J)+(T(K,J)*Y1(K))/S1(K)
YT(N2,J)=1.
133 C(J)=C(J)-Y2(J)
C(N2)=0.
C
C *************************************************************
C INVERSE BY LINEAR ROW TRANSFORMATION
C
DO 202 I=1,N2
DO 201 J=1,N2
201 YS(I,J)=0
202 YS(I,I)=1.
DO 210 I=1,N2
IF(((YT(I,I)+100.)-100.).NE.0.0) GO TO 220
IF(I.EQ.N2) GO TO 260
DO 211 J=I+1,N2
IF(((YT(J,I)+100.)-100.).NE.0.0) GO TO 212
211 CONTINUE
260 WRITE(IDLG,261)
261 FORMAT('0INVERSE DOES NOT EXIST')
GO TO 1000
212 DO 213 K=1,N2
YT(I,K)=YT(I,K)+YT(J,K)
213 YS(I,K)=YS(I,K)+YS(J,K)
220 GOB=YT(I,I)
DO 221 J=1,N2
YT(I,J)=YT(I,J)/GOB
221 YS(I,J)=YS(I,J)/GOB
DO 230 L=1,N2
IF(L.EQ.I) GO TO 230
GOB=YT(L,I)
DO 231 J=1,N2
YT(L,J)=YT(L,J)-GOB*YT(I,J)
231 YS(L,J)=YS(L,J)-GOB*YS(I,J)
230 CONTINUE
210 CONTINUE
C ***********************************************************
C
270 DO 134 J=1,N2
DO 134 K=1,N2
134 B1(J)=B1(J)+YS(J,K)*C(K)
DO 138 I=1,N1
DO 135 J=1,N2
135 AM(I)=AM(I)+T(I,J)*B1(J)
AM(I)=Y1(I)-AM(I)
IF(S1(I))137,136,137
136 TYPE 127,I,S1(I)
137 AM(I)=AM(I)/S1(I)
138 AN=AN+AM(I)
P1=N1
AN=AN/P1
DO 139 I=1,N1
A=AM(I)-AN
139 RAB=RAB+A*Y1(I)
DO 140 J=1,N2
140 RAB=RAB+B1(J)*Y2(J)
RAB=RAB+AN*G-G1
SI=SSC-RAB
P2=N2
TI=TC-P1-P2+1.
MP=0
MP=MP+1
IF(TI)141,151,141
141 SMI=SI/TI
TAB=RAB-SSBA
MP=MP+1
IF(P1-1.)142,151,142
142 TMAB=TAB/(P1-1.)
TBA=RAB-SSAB
MP=MP+1
IF(P2-1.)143,151,143
143 TMBA=TBA/(P2-1.)
MP=MP+1
NDG=TW
IF(ITYPE.EQ.1) NDG=TI
IF (ITYPE.EQ.1) SMW=SMI
IF(SMW)144,151,144
144 FI=SMI/SMW
FA=TMAB/SMW
FB=TMBA/SMW
IF(IOUT.NE.21) GO TO 328
LINES=LINES+10
IF(LINES.LE.LINPP) GO TO 328
CALL PRNTHD
LINES=10
328 WRITE (IOUT,145)
145 FORMAT(//29X11HFINAL ANOVA)
WRITE (IOUT,118)
KK=TC-1.
WRITE (IOUT,119)KK,SSC
KK=P1-1.
I=1
J=2
PROB=FISHER(KK,NDG,FA)
WRITE (IOUT,146)I,J,KK,TAB,TMAB,FA,PROB
146 FORMAT(3XI1,12H ELIMINATING,I2,2X,I5,3F10.2,2X,F6.4)
KK=P2-1.
PROB=FISHER(KK,NDG,FB)
WRITE (IOUT,147)J,I,KK,TBA,TMBA,FB,PROB
147 FORMAT(3XI1,12H ELIMINATING,I2,2X,I5,3F10.2,2X,F6.4)
KK=TI
IF(ITYPE.EQ.0) PROB=FISHER(KK,NDG,FI)
IF(ITYPE.EQ.0)WRITE (IOUT,148)I,J,KK,SI,SMI,FI,PROB
IF(ITYPE.EQ.1) WRITE(IOUT,148)I,J,KK,SI,SMI
148 FORMAT(6XI1,3H BY,I2,8X,I5,3F10.2,2X,F6.4)
KK=TW
IF(ITYPE.EQ.0)WRITE (IOUT,121)KK,SSW,SMW
KK=SS-1.
WRITE (IOUT,122)KK
GO TO 1000
C VALUES OF NP MAY BE 1,2,3
C 1 IMPLIES TC-1=0,2 IMPLIES TW=0,3 IMPLIES SMW=0
149 TYPE 150,NP
150 FORMAT(3HNP=,I2)
GO TO(104,105,103),NP
C VALUES OF MP MAY BE 1,2,3,4
C 1 IMPLIES TI=0,2 IMPLIES P1-1=0,3 IMPLIES P2-1=0,4 IMPLIES SMW=0
151 TYPE152,MP
152 FORMAT(3HMP=,I2)
GO TO(141,142,143,144),MP
1000 IF(OPTD(1).NE.1) GO TO 400
2000 CONTINUE
RETURN
END
C *** STAT PACK ***
C SUBROUTINE FOR FRIEDMAN TWO-WAY ANALYSES OF VARIANCE
C CALLING SEQUENCE: CALL FRIED(NV,NC,MV,MC,DATA,SUM,RANK,NAMES)
C WHERE: NV - NUMBER OF VARIABLES USED
C NC - NUMBER OF OBSERVATIONS USED
C MV - MAXIMUM NUMBER OF VARIABLES ALLOWED
C MC - MAXIMUM NUMBER OF OBSERVATIONS ALLOWED
C DATA - MATRIX CONTAINING DATA
C SUM - EXTRA VECTOR AT LEAST NV LONG
C RANK - EXTRA VECTOR AT LEAST NV LONG
C NAMES - VECTOR CONTAINING VARIABLE NAMES
C
C ROUTINE SUGGESTED TO BE INCLUDED IN STAT PACK BY LONNIE
C HANNAFORD (SPECIAL EDUCATION) AND ULDIS SMIDCHENS (EDUCATION).
C SOURCE IS NON PARAMETRIC STATISTICS BY SIEGEL PAGES 166-172.
C
SUBROUTINE FRIED(NV,NC,MV,MC,DATA,SUM,RANK,NAMES)
DIMENSION DATA(MC,MV),SUM(1),RANK(1),IV(20),IVB(20),NAMES(1)
COMMON /DEV/ ICC,IDATA,IOUT,IDLG,IDSK
COMMON /PRNT/ LINPP,ICOPS,RUNPRG
COMMON /EXTRA/ HEDR(70),NSZ
1 IF(ICC.NE.2) WRITE(IDLG,2)
2 FORMAT(' WHICH VARIABLES? '$)
CALL ALPHA(IV,20,N,IRET,IHELP,IERR,NAMES,NV)
IF(IRET.EQ.1) RETURN
IF(IERR.EQ.1) GO TO 1
IF(IHELP.NE.1) GO TO 4
WRITE(IDLG,3)
3 FORMAT('0ENTER UP TO 20 VARIABLES FOR THE FRIEDMAN TWO-WAY'/
1' ANALYSIS OF VARIANCE. VARIABLE NAMES OR NUMBERS'/
2' MAY BE USED TO INDICATED THE VARIABLES. RANGES OF'/
3' VARIABLES MAY BE ENTERED BY TYPING THE EXTREMES OF THE RANGE'/
4' SEPARATED BY A MINUS. THE ASTERIK MAY BE USED IN ONE'/
5' OR MORE POSITIONS TO INDICATED ALL POSSIBLE COMBINATIONS.')
GO TO 1
4 IF(N.LT.2) RETURN
DO 5 I=1,N-1
IF(IV(I).LE.0) GO TO 5
DO 6 J=I+1,N
IF(IV(I).NE.IV(J)) GO TO 6
WRITE(IDLG,7) NAMES(IV(J))
7 FORMAT(' SAME VARIABLE (',A5,') LISTED TWICE')
GO TO 1
6 CONTINUE
5 CONTINUE
C
C OK NOW INSERT FOR * AND ALL.
C
K=1
DO 8 I=1,N
IVB(I)=IV(I)
IF(IV(I).GT.0) GO TO 8
IVB(I)=K
K=K+1
8 CONTINUE
GO TO 18
C
C RETURN HERE TO PICK UP NEXT SET OF VARIABLES
C
10 J=N
14 IF(IV(J).GT.0) GO TO 15
IVB(J)=IVB(J)+1
IF(IVB(J).LE.NV) GO TO 16
15 J=J-1
IF(J.GE.1) GO TO 14
RETURN
16 K=IVB(J)
IF(J.EQ.N) GO TO 18
DO 17 I=J+1,N
IF(IV(I).GT.0) GO TO 17
K=K+1
IF(K.GT.NV) GO TO 15
IVB(I)=K
17 CONTINUE
18 DO 13 I=1,N-1
DO 13 K=I+1,N
IF(IVB(I).EQ.IVB(K)) GO TO 14
13 CONTINUE
C
C BEGIN THE FRIEDMAN ANALYSIS OF VARIANCE
C
DO 19 I=1,N
19 SUM(I)=0
DO 20 J=1,NC
DO 21 I=1,N
21 RANK(I)=DATA(J,IVB(I))
C
C RANKING ( ONLY 20 MAX, SO SORT RANK METHOD NOT USED )
C
DO 22 I=1,N
SAME=1.
ALESS=0
DO 23 K=1,N
IF(RANK(K).GT.RANK(I)) GO TO 23
IF(RANK(K).LT.RANK(I)) GO TO 24
SAME=SAME+1
GO TO 23
24 ALESS=ALESS+1.
23 CONTINUE
22 SUM(I)=SUM(I)+ALESS+SAME/2.
20 CONTINUE
C
C RANKING COMPLETE NOW CALCULATE THE CHI SQUARE
C
IF(IOUT.NE.21) WRITE(IOUT,5566) (HEDR(J),J=1,NSZ)
5566 FORMAT('1',70A1)
IF(IOUT.EQ.21) CALL PRNTHD
WRITE(IOUT,31)
31 FORMAT('0',8X,'**** FRIEDMAN TWO-WAY ANALYSIS OF VARIANCE ****'/
1'0VARIABLE SUM OF RANKS'/' ======== ============')
LINES=7
TOTSQ=0
DO 30 I=1,N
IF(IOUT.NE.21) GO TO 36
LINES=LINES+1
IF(LINES.LE.LINPP) GO TO 36
CALL PRNTHD
WRITE(IOUT,31)
LINES=8
36 WRITE(IOUT,32) NAMES(IVB(I)),SUM(I)
32 FORMAT(2X,A5,5X,F10.2)
30 TOTSQ=TOTSQ+SUM(I)**2
FRIED=(12./(NC*N*(N+1.)))*TOTSQ-3.*NC*(N+1.)
IDF=N-1
FD=FRIED/IDF
PROB=FISHER(IDF,1000,FD)
IF(IOUT.NE.21) GO TO 37
LINES=LINES+3
IF(LINES.LE.LINPP) GO TO 37
CALL PRNTHD
37 WRITE(IOUT,35) FRIED,IDF,PROB
35 FORMAT('0CHI SQUARE=',F12.3,', WITH ',I2,' DEGREES OF FREEDOM'/
1' HAVING A PROBABILITY OF ',F5.2)
GO TO 10
END
C *** STAT PACK ***
C SUBROUTINE FOR SIGN TEST.
C CALLING SEQUENCE: CALL SIGNT(NV,NC,MV,MC,DATA,NAMES)
C WHERE: NV - NUMBER OF VARIABLES USED.
C NC - NUMBER OF OBSERVATIONS USED.
C MV - MAXIMUM NUMBER OF VARIABLES ALLOWED
C MC - MAXIMUM NUMBER OF OBSERVATIONS ALLOWED
C DATA - MATRIX FOR DATA
C NAMES - VECTOR CONTAINING VARIABLE NAMES.
C
C ROUTINE SUGGESTED BY LONNIE HANNAFORD (SPECIAL EDUCATION) AND
C ULDIS SMIDCHENS (TEACHER EDUCATION). SOURCE IS NON-
C PARAMETRIC STATISTICS BY SIEGEL PAGES 68-75. THE
C BINOMIAL EXPANSION COULD BE USED FOR N>25, HOWEVER IT IS
C SLOW FOR LARGE N. WORKS BY ATTEMPTING TO NORMALIZE NUMBER
C AS IT IS BEING CALCULATED ABOUT THE VALUE 1.
C
SUBROUTINE SIGNT(NV,NC,MV,MC,DATA,NAMES)
DIMENSION DATA(MC,MV),NAMES(1)
COMMON /DEV/ ICC,IDATA,IOUT,IDLG,IDSK
COMMON /PRNT/ LINPP,ICOPS,RUNPRG
COMMON /EXTRA/ HEDR(70),NSZ
C EQUIVALENCE FOR BINARY EXPANSION USED TO FIND POWER OF 2
EQUIVALENCE (ISUM,SUM)
IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(I),I=1,NSZ)
5566 FORMAT('1',70A1)
IF(IOUT.EQ.21) CALL PRNTHD
WRITE(IOUT,1)
1 FORMAT('0',20X,'**** SIGN TEST ****'/'0',35X,'FREQUENCY OF'/
119X,'NUMBER OF',9X,'SIGN WHICH'/19X,'PAIRS WITH',6X,
2'OCCURED LEAST',7X,'Z',10X,'PROB'/' VAR',4X,'VAR',4X,
3'SIGNED DIFFERENCE',7X,'OFTEN',8X,'(N>25)',6X,'(N<=25)')
LINES=9
DO 2 I=1,NV-1
DO 2 J=I+1,NV
N=0
NMIN=0
DO 3 K=1,NC
IF(DATA(K,I)-DATA(K,J)) 4,3,5
4 NMIN=NMIN+1
5 N=N+1
3 CONTINUE
IF(NMIN.GT.(N-NMIN)) NMIN=N-NMIN
IF(N.GT.25) GO TO 20
C
C BINOMIAL EXPANSION (N<=25)
C
XSUM=0
IF(N.LT.127) XSUM=1./2.**N
IF(NMIN.LT.1) GO TO 13
DO 11 L=1,NMIN
SUM=1.
IN=N-L+1
XONE=1.
NTWOS=N
DO 12 K=IN,N
SUM=SUM*K/XONE
IF(NTWOS.LT.1) GO TO 12
IF(SUM.LT.2) GO TO 12
JJ=(ISUM.AND."377777777777)/2**27-129
NDIV=JJ
IF(JJ.GT.NTWOS) NDIV=NTWOS
SUM=SUM/2.**NDIV
NTWOS=NTWOS-NDIV
12 XONE=XONE+1.
IF(NTWOS.GT.127) GO TO 11
IF(NTWOS.GT.0) SUM=SUM/2.**NTWOS
XSUM=XSUM+SUM
11 CONTINUE
13 IF(IOUT.NE.21) GO TO 15
LINES=LINES+1
IF(LINES.LE.LINPP) GO TO 15
CALL PRNTHD
WRITE(IOUT,1)
LINES=10
15 WRITE(IOUT,14) NAMES(I),NAMES(J),N,NMIN,XSUM
14 FORMAT(1X,A5,2X,A5,8X,I5,13X,I5,20X,F7.4)
GO TO 2
C
C BINOMIAL DISTRIBUTION GIVE Z VALUE
C
20 A=N
XMN=A/2.
XSTD=SQRT(A)/2.
Z=N-XMN
IF(N.LT.XMN) Z=Z+.5
Z=Z/XSTD
IF(IOUT.NE.21) GO TO 22
LINES=LINES+1
IF(LINES.LE.LINPP) GO TO 22
CALL PRNTHD
WRITE(IOUT,1)
LINES=10
22 WRITE(IOUT,21) NAMES(I),NAMES(J),N,NMIN,Z
21 FORMAT(1X,A5,2X,A5,8X,I5,13X,I5,7X,G12.3)
2 CONTINUE
RETURN
END