Trailing-Edge
-
PDP-10 Archives
-
decuslib10-09
-
43,50466/nonpar.f4
There are no other files named nonpar.f4 in the archive.
C WESTERN MICHIGAN UNIVERSITY
C NONPAR.F4 (FILE NAME ON LIBRARY DECTAPE)
C NONPAR, 1.1.1 (CALLING NAME, SUBLST NO.)
C NON-PARAMETRIC STATISTICS
C THIS IS A COMBINATION OF A PROGRAM GIVEN BY WAYNE STATE
C UNIVERSITY WITH SUBSTANTIAL REVISIONAL AND ADDITIONAL
C PROGRAMMING BY MR. SAM ANEMA. LATER IT WAS MODIFIED BY
C RUSS BARR.
C LIBRARY DECTAPE PROGS. USED: USAGE.MAC
C FORWMU PROGS. USED: TTYPTY, DEVCHG, EXISTS, PRINTS
C APLIB PROGS. USED: IO, GETFOR, CHIPRB, NORMCV
C INTERNAL SUBR. USED: KISQ, ONSMP, SNTST, WLTST, MDTST,
C MNWHU, WRTST, MSTST, COCHQ, FRANV, KWANV
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
DIMENSION VAR(20,200),NOB(20),IFMT(80)
COMMON VAR,NVAR,NOBS,NOB,NCODE,IND,NOUTD,IOUT
INP=2
IOUT=3
IND=-4
NOUTD=-1
WRITE(NOUTD,1)
1 FORMAT('1',20X,'WMU NONPARAMETRIC PROGRAM')
C CALL USAGE('NONPAR')
WRITE(NOUTD,2)
2 FORMAT(///)
C---------------TTYPTY RETURNS ZERO - TTY JOB, MINUS ONE - BATCH JOB
CALL TTYPTY(ICODER)
C---------------IDV1, IDV2 ARE RETURNED. OTHER ARGS. ARE INPUT.
C--------------- 1 MEANS OUTPUT? PRINTS, 0 MEANS INPUT? PRINTS
CALL IO(IOUT,IDV1,NOUTD,IND,1,ICODER)
34 CALL IO(INP,IDV2,NOUTD,IND,0,ICODER)
39 WRITE(NOUTD,4)
4 FORMAT(' THE NUMBER OF VARIABLES = ',$)
READ(IND,5)NVAR
5 FORMAT(I)
IF(NVAR.GT.20.OR.NVAR.LT.1)GO TO 37
C---------------IFMT, ISTD ARE RETURNED. OTHER ARGS. ARE INPUT.
C---------------80=NO. OF OBJ. TIME FORMAT WORDS (5 LINES). 2 MEANS
C--------------- F-TYPE ONLY FORMAT.
CALL GETFOR(IND,NOUTD,IFMT,ISTD,80,2)
IF(ISTD.EQ.1)IFMT(1)='(10F)'
WRITE(NOUTD,12)
12 FORMAT(' MISSING DATA SYMBOL = ',$)
READ(IND,7)FINAL
7 FORMAT(F)
IF(ISTD.EQ.1)GO TO 42
IF(IDV2.EQ.'TTY')WRITE(NOUTD,40)
40 FORMAT(' ENTER DATA'/)
GO TO 43
42 IF(IDV2.EQ.'TTY')WRITE(NOUTD,44)
44 FORMAT(' ENTER DATA (NO MORE THAN 10 PER LINE)'/)
43 NOBS = 0
6 NOBS=NOBS+1
READ(INP,IFMT,END=300)(VAR(I,NOBS),I=1,NVAR)
IF(NOBS.GT.200)GO TO 10
GO TO 6
300 NOBS=NOBS-1
WRITE(NOUTD,302)NOBS
302 FORMAT(1X,I4,' CASES WERE ENTERED')
DO 15 I=1,NVAR
JA=1
DO 16 J=1,NOBS
IF(VAR(I,J).EQ.FINAL)GO TO 16
VAR(I,JA)=VAR(I,J)
JA=JA+1
16 CONTINUE
NOB(I)=JA-1
15 CONTINUE
GO TO 36
10 WRITE(NOUTD,11)
11 FORMAT(' TOO MANY OBSERVATIONS'/)
GO TO 34
37 WRITE(NOUTD,38)
38 FORMAT(' ILLEGAL NO. OF VARIABLES'/)
GO TO 39
36 WRITE(NOUTD,22)
22 FORMAT(' ENTER CODE FOR TEST DESIRED.'/)
READ(IND,5)KOD
GO TO (23,24,25,26,27,28,29,30,31,32,33,34),KOD
WRITE(NOUTD,35)
35 FORMAT(' ILLEGAL CODE'/)
GO TO 36
23 CALL KISQ
GO TO 36
24 CALL ONSMP
GO TO 36
25 CALL SNTST
GO TO 36
26 CALL WLTST
GO TO 36
27 CALL MDTST
GO TO 36
28 CALL MNWHU
GO TO 36
29 CALL WRTST
GO TO 36
30 CALL MSTST
GO TO 36
31 CALL COCHQ
GO TO 36
32 CALL FRANV
GO TO 36
33 CALL KWANV
GO TO 36
END
C
C THIS SUBROUTINE WILL FOR A CONTINGENCY TABLE FROM
C RAW DATA AND CALCULATE CHI-SQUARE,CONTINGENCY COEFFICIENT,
C PHI-SQUARE, AND PHI-PRIME.
C
C
C---------------NVAR, NOBS, INP, NOUTD, IOUT ARE INPUT THRU COMMON
SUBROUTINE KISQ
DIMENSION NOB(20),SUMR(0/20),SUMC(0/20),VAR(20,200),CHI(0/20,0/20)
DIMENSION NCROSS(40),MM(2),NN(2)
COMMON VAR,NVAR,NOBS,NOB,NCODE,INP,NOUTD,IOUT
WRITE(IOUT,23)
23 FORMAT(10X,46('*'))
WRITE(IOUT,24)
24 FORMAT(/25X,'CHI-SQUARE'/)
WRITE (IOUT,23)
6 WRITE(NOUTD,1)
1 FORMAT(' ENTER PAIRS OF VARIABLES TO BE ANALYZED.'/)
READ(INP,2)NCROSS
2 FORMAT(40I)
DO 3 I=40,1,-1
IF(NCROSS(I).NE.0)GO TO 4
3 CONTINUE
GO TO 6
4 NPR=I/2
DO 5 I=1,NPR*2
5 IF(NCROSS(I).LT.1.OR.NCROSS(I).GT.NVAR)GO TO 90
DO 10 I=1,NPR
I1=I*2
I2=I1-1
K1=NCROSS(I1)
K2=NCROSS(I2)
DO 20 J=0,20
SUMR(J)=0.0
SUMC(J)=0.0
DO 21 K=0,20
21 CHI(J,K)=0.0
20 CONTINUE
DO 22 K=1,NOBS
J1=VAR(K1,K)
J2=VAR(K2,K)
IF(J1.GT.20.OR.J1.LT.0)GO TO 25
IF(J2.GT.20.OR.J2.LT.0)GO TO 25
CHI(J1,J2)=CHI(J1,J2)+1.0
SUMR(J1)=SUMR(J1)+1
SUMC(J2)=SUMC(J2)+1
22 CONTINUE
DO 30 J=20,0,-1
IF(SUMC(J).NE.0.0)GO TO 31
30 CONTINUE
31 NC=J
DO 32 J=20,0,-1
IF(SUMR(J).NE.0.0)GO TO 33
32 CONTINUE
33 NR=J
WRITE(IOUT,34)K2,K1
34 FORMAT(//5X,'ROW VARIABLE',I3,' VS COL VARIABLE',I3)
WRITE(IOUT,27)
27 FORMAT(4X,'ROW',13X,'CONTINGENCY TABLE'/' TOTAL'/)
WRITE(IOUT,38)(J,J=0,NC)
38 FORMAT(9X,15I4)
DO 35 J=0,NR
35 WRITE(IOUT,37)J,SUMR(J),(CHI(J,K),K=0,NC)
37 FORMAT(I3,F5.0,2X,15F4.0/(10X,15F4.0))
WRITE(IOUT,39)(SUMC(J),J=0,NC)
39 FORMAT(' COL'/3X,'TOTAL',2X,15F4.0/(10X,15F4.0))
C
C CALCULATE STATISTICS
C
CH=0.0
TOT=0.0
DO 50 J=0,NR
50 TOT=TOT+SUMR(J)
DO 51 J=0,NR
DO 52 K=0,NC
E=(SUMR(J)*SUMC(K))/TOT
IF(E.EQ.0.0)GO TO 52
CH=CH+(CHI(J,K)-E)**2/E
52 CONTINUE
51 CONTINUE
A=0.0
B=0.0
DO 53 J=0,20
IF(SUMC(J).NE.0.0)A=A+1.0
IF(SUMR(J).NE.0.0)B=B+1.0
53 CONTINUE
DF=(A-1.0)*(B-1.0)
C=SQRT(CH/(TOT+CH))
PHISQ=CH/TOT
BMAX=AMAX1(A,B)
IF(BMAX.EQ.1.0)GO TO 60
PHIPR=SQRT(PHISQ/(BMAX-1.0))
GO TO61
60 PHIPR=1.0
C---------------CH, INT ARE INPUT; PCH RETURNED
61 CALL CHIPRB(CH,INT(DF),PCH)
WRITE(IOUT,62)CH,PCH,DF,C,PHISQ,PHIPR
62 FORMAT(/' CHI-SQUARE = ',F12.4,/,' CHI SQUARE PROB. =',F8.5,
1/' DEGREES OF FREEDOM = ',F7.0/
1' CONTINGENCY COEFFICIENT = ',F12.4/' PHI-SQUARE = ',F12.4/
2' PHI-PRIME = ',F12.4)
IF(A.NE.2.0.OR.B.NE.2.0)GO TO 10
M1=0
M2=0
DO 56 J=0,20
IF(SUMC(J).NE.0.0)GO TO 54
57 IF(SUMR(J).NE.0.0)GO TO 55
GO TO 56
54 M1=M1+1
NN(M1)=J
GO TO 57
55 M2=M2+1
MM(M2)=J
56 CONTINUE
A=CHI(MM(1),NN(1))
B=CHI(MM(1),NN(2))
C=CHI(MM(2),NN(1))
D=CHI(MM(2),NN(2))
CH=TOT*(ABS(A*D-B*C)-TOT/2.0)**2/((A+B)*(C+D)*(A+C)*(B+D))
WRITE(IOUT,63)CH
63 FORMAT(' CORRECTED CHI-SQUARE = ',F12.4/)
10 CONTINUE
RETURN
25 WRITE(NOUTD,26)I1,I2
26 FORMAT(' A DATA ITEM GREATER THAN 20 WAS ENCOUNTERED WHILE
1TRYING TO CREATE THE CONTINGENCY TABLE BETWEEN VARIABLES',I4,
2' AND ',I4)
GO TO 10
90 WRITE(NOUTD,91)
91 FORMAT(' VARIABLE INDEX IS OUT OF RANGE'/)
GO TO 6
END
C---------------VAR, NVAR, NOB, IOUT ARE INPUT THRU COMMON
SUBROUTINE ONSMP
DIMENSION VARD(200),NOB(20)
DIMENSION VAR(20,200)
COMMON VAR,NVAR,NOBS,NOB,NCODE,INP,NOUTD,IOUT
C ONE SAMPLE RUNS TEST
WRITE (IOUT,9010)
9010 FORMAT(10X,46('*'))
WRITE (IOUT,9040)
9040 FORMAT(/23X,'ONE SAMPLE RUNS TEST'/)
WRITE (IOUT,9010)
C DO TEST FOR EACH VARIABLE SUPPLIED
DO 9990 I=1,NVAR
NO=NOB(I)
DO 9100 J=1,NO
9100 VARD(J)=VAR(I,J)
C FIND MEDIAN
9150 N1=NO-1
DO 9200 J=1,N1
J1=J+1
DO 9200 K=J1,NO
IF (VARD(J)-VARD(K))9200,9200,9170
9170 DUM=VARD(J)
VARD(J)=VARD(K)
VARD(K)=DUM
9200 CONTINUE
NOH=NO/2
IF (2*NOH-NO)9230,9220,9230
9220 XMED=(VARD(NOH)+VARD(NOH+1))/2.0
GO TO 9240
9230 XMED=VARD(NOH+1)
C FIND N1 AND N2
9240 NI=0
NS=0
NN=0
DO 9300 J=1,NO
IF (VAR(I,J)-XMED)9260,9270,9280
9260 NI=NI+1
VARD(J)=-1.0
GO TO 9300
9270 VARD(J)=0.0
NN=NN+1
GO TO 9300
9280 NS=NS+1
VARD(J)=1.0
9300 CONTINUE
C ADD THOSE EQUAL TO MEDIAN TO LESS FREQUENT COUNT
IF (NI-NS)9350,9350,9360
9350 DUM=-1.0
NI=NI+NN
GO TO 9400
9360 DUM=1.0
NS=NS+NN
9400 DO 9440 J=1,NO
IF (VARD(J))9440,9420,9440
9420 VARD(J)=DUM
9440 CONTINUE
C COUNT RUNS
R=1.0
DO 9490 J=1,N1
IF (VARD(J)-VARD(J+1))9450,9490,9450
9450 R=R+1.0
9490 CONTINUE
C FINAL RESULTS
WRITE (IOUT,9500)I,NI,NS,R
9500 FORMAT(//10X,'RESULTS FROM SAMPLE',I3/12X,'N1 = ',I3,
13X,'N2 = ',I3,3X,'NUMBER OF RUNS = ',F3.0)
IF (NI-20)9510,9540,9540
9510 IF (NS-20)9990,9540,9540
C RESULTS IF N1 OR N2 EXCEED 19
9540 A12=NI+NS
T12=NI*NS
FMEAN=((2.0*T12)/A12)+1.0
STD=(2.0*T12*(2.0*T12-A12))/(A12*A12*(A12-1.0))
STD=SQRT(STD)
Z=(R-FMEAN)/STD
WRITE (IOUT,9550)FMEAN,STD,Z
9550 FORMAT(12X,'MEAN = ',F8.5,4X,'ST. DEV. = ',F9.4,4X,
1'Z = ',F9.5)
9990 CONTINUE
400 RETURN
10000 RETURN
END
C---------------IOUT, NVAR, NOBS ARE INPUT THRU COMMON
SUBROUTINE SNTST
DIMENSION NOB(20),VAR(20,200)
COMMON VAR,NVAR,NOBS,NOB,NCODE,INP,NOUTD,IOUT
C ********* SIGN TEST *********
WRITE (IOUT,6010)
6010 FORMAT (1H0,10X,46H**********************************************)
WRITE (IOUT,6000)
6000 FORMAT (1H0,30X,9HSIGN TEST)
WRITE (IOUT,6010)
NVAR1=NVAR-1
DO 6500 I=1,NVAR1
I1=I+1
DO 6500 II=I1,NVAR
WRITE (IOUT,6050)I,II
6050 FORMAT (10H0 VAR(I2,10H) VS. VAR(,I2,1H))
C FIND SIGN OF DIFFERENCE IN EACH PAIR, STORE IN SIGNP OR SIGNN
SIGNP=0.0
SIGNN=0.0
DO 6150 J=1,NOBS
IF (VAR(I,J)-VAR(II,J))6120,6150,6130
6120 SIGNN=SIGNN+1.
GO TO 6150
6130 SIGNP=SIGNP+1.
6150 CONTINUE
C FIND LESS FREQUENT SIGN
IF (SIGNP-SIGNN)6220,6230,6230
6220 NLESS=SIGNP
GO TO 6250
6230 NLESS=SIGNN
6250 NTOT=SIGNN+SIGNP
C EVALUATE PROBABILITY
WRITE (IOUT ,6320)NTOT
IF (NTOT-25)6300,6300,6400
C FOR N LESS THAN OR EQUAL TO 25
6300 P=0.
DO 6310 K=0,NLESS
6310 P=P+BINO(NTOT,K)
P=P/2.0**NTOT
6320 FORMAT (45HL NUMBER OF PAIRS WITH SIGNED DIFFERENCE = ,I3)
WRITE (IOUT,6340)NLESS,P
6340 FORMAT (35H PROBABILITY OF OBSERVED FREQUENCY,I3,8X,2HIS,F12.7)
GO TO 6500
C FOR N GREATER THAN 25
6400 TOT=NTOT
FMEAN=TOT/2.0
SD=SQRT(TOT)/2.0
EESS=NLESS
IF (EESS-FMEAN)6410,6410,6420
6410 EESS=EESS+.5
GO TO 6430
6420 EESS=EESS-.5
6430 Z=(EESS-TOT/2.)/SD
WRITE (IOUT,6450)NLESS,FMEAN,SD,Z
6450 FORMAT(' ','OBSERVED FREQUENCY IS ',I3/'MEAN = ',F8.5,
1'ST. DEV. = ',F7.5,3X,'Z = ',F9.5)
6500 CONTINUE
400 RETURN
10000 RETURN
END
C---------------IOUT, NOB, NVAR, NOBS, VAR ARE INPUT THRU COMMON
SUBROUTINE WLTST
DIMENSION NOB(20),DIF(200)
DIMENSION RANK(200),VAR(20,200)
COMMON VAR,NVAR,NOBS,NOB,NCODE,INP,NOUTD,IOUT
C THE WILCOXON MATCHED-PAIRS SIGNED-RANKS TEST
C CALCULATES T, MEAN, SIGMA AND N FOR EACH PAIR OF VARIABLES
C N=NUMBER OF PAIRS WITH DIFFERENCE NOT EQUAL TO 0
C NVAR=TOTAL NUMBER OF VARIABLES
C NOBS=TOTAL NUMBER OF OBSERVATION
C PRINT TITLE
WRITE (IOUT,2000)
2000 FORMAT(10X,46('*'))
WRITE (IOUT,1050)
1050 FORMAT(/13X,'WILCOXON MATCHED-PAIRS SIGNED-RANKS TEST'/)
WRITE (IOUT,2000)
NOBS=NOB(1)
C BEGIN MAIN CALCULATIONS
NVA=NVAR-1
DO 1350 I=1,NVA
I1=I+1
DO 1350 J=I1,NVAR
C FIND DIFFERENCE OF MATCHED PAIRS - DIF(K)
DO 1150 K=1,NOBS
1150 DIF(K)=VAR(I,K)-VAR(J,K)
C ARRANGE DIFFERENCES IN RANK ORDER
N1=NOBS-1
DO 1170 K=1,N1
K1=K+1
DO 1170 K2=K1,NOBS
IF (ABS(DIF(K))-ABS(DIF(K2)))1170,1170,1160
1160 DUM=DIF(K)
DIF(K)=DIF(K2)
DIF(K2)=DUM
1170 CONTINUE
C DISCARD DIF=0
N=NOBS
LOW=1
DO 1180 K=1,NOBS
IF (DIF(K))1185,1175,1185
1175 N=N-1
1180 LOW=LOW+1
C CREATE RANK
1185 DO 1190 K=1,N
L=LOW+K-1
1190 RANK(L)=K
C ASSIGN RANKS AND EVALUATE TIES
IEQ=1
DO 1230 K=LOW,N1
IF (ABS(DIF(K))-ABS(DIF(K+1)))1220,1200,1220
1200 GO TO (1201,1210),IEQ
1201 RSUM=RANK(K)+RANK(K+1)
TIE=2.0
IEQ=2
GO TO 1230
1210 RSUM=RSUM+RANK(K+1)
TIE=TIE+1.0
GO TO 1230
1220 GO TO (1230,1221),IEQ
1221 IEQ=1
NTIE=TIE
RANKS=RSUM/TIE
DO 1225 K1=1,NTIE
L=K-K1+1
1225 RANK(L)=RANKS
1230 CONTINUE
GO TO (1234,1221),IEQ
C ASSIGN SIGN TO RANKS AND SUM
1234 SUMP=0.0
SUMN=0.0
DO 1250 K=LOW,NOBS
RANK(K)=(RANK(K)*DIF(K))/ABS(DIF(K))
IF (RANK(K))1235,1999,1240
1235 SUMN=SUMN+RANK(K)
GO TO 1250
1240 SUMP=SUMP+RANK(K)
1250 CONTINUE
C MAKE FINAL CALCULATIONS AND PRINT RESULTS
IF (SUMP-SUMN)1255,1255,1265
1255 T=SUMP
GO TO 1280
1265 T=ABS(SUMN)
1280 D=N
TMEAN=(D*D+D)/4.0
SIGMA=SQRT(TMEAN*(2.0*D+1.0)/6.0)
Z=(T-TMEAN)/SIGMA
C---------------Z INPUT, ZPRB RETURNED
CALL NORMCV(Z,ZPRB)
WRITE (IOUT,1300),I,J,T,TMEAN,SIGMA,Z,ZPRB,N
1300 FORMAT(//' ','LT(',I2,',',I2,') = ',F7.3,5X,'MEAN = ',F7.3,5X,
1'ST. DEV. = ',F8.4/3X,'Z = ',F8.4,5X,
2'Z-PROB =',F8.5,5X,'N = ',I3)
1350 CONTINUE
GO TO 400
1999 WRITE (IOUT,1950)
1950 FORMAT (20H INVALID CALCULATION)
400 RETURN
10000 RETURN
END
C---------------VAR RETURNED THRU COMMON; NVAR, NOBS, NOB,
C--------------- IOUT ARE INPUT THRU COMMON
SUBROUTINE MDTST
DIMENSION NOB(20),EXP(20)
DIMENSION VAR(20,200),TABLE(20,2)
COMMON VAR,NVAR,NOBS,NOB,NCODE,INP,NOUTD,IOUT
CALL OFILE(20,'DATA')
C MEDIAN AND EXTENSION OF MEDIAN TESTS
C DATA TO BE READ INTO AN ARRAY NAMED VAR
C NOBS=NUMBER OF OBSERVATIONS
C NVAR=NUMBER OF VARIABLES
C THE ARRY NAMED TABLE IS USED TO CALCULATE CHI
C NOB(I)=NUMBER OF OBSERVATIONS IN VARIABLE I
C NEXT = NEXT PROGRAM TO BE USED
C PRINT TITLE
WRITE (IOUT,2000)
2000 FORMAT(10X,46('*'))
WRITE (IOUT,2005)
2005 FORMAT(/15X,'MEDIAN AND EXTENSION OF MEDIAN TESTS'/)
WRITE (IOUT,2000)
DO12000I=1,NVAR
12000 WRITE(20)(VAR(I,J),J=1,NOBS)
C ORDER DATA
DO 2170 I=1,NVAR
NO=NOB(I)
DO 2170 J=1,NO
JJ=J
DO 2170 M=I,NVAR
NM=NOB(M)
DO 2165 N=JJ,NM
IF(VAR(I,J)-VAR(M,N))2165,2165,2160
2160 DUM=VAR(I,J)
VAR(I,J)=VAR(M,N)
VAR(M,N)=DUM
2165 CONTINUE
JJ=1
2170 CONTINUE
C FIND MEDIAN
MD=0
DO 2171 I=1,NVAR
2171 MD=MD+NOB(I)
M=0
DUM=MD
DUM=DUM/2.0
MD=MD/2
UD=MD
MUD=2.*UD+.0001
MUDD=2.*DUM+.0001
IF(MUD-MUDD)2172,2173,900
2172 MD=MD+1
M=1
C M=0 INDICATES MEDIAN REQUIRES AVERAGE OF 2 VARIABLES
2173 DO 2175 I=1,NVAR
MD=MD-NOB(I)
IF (MD)2178,2176,2175
2175 CONTINUE
C IF MEDIAN IS LAST VALUE FOR A VARIABLE - BRANCH HERE
2176 NO=NOB(I)
VMDN=VAR(I,NO)
IF (M)900,2177,2180
2177 VMDN=(VMDN+VAR(I+1,1))/2.0
GO TO 2180
C OTHERWISE - BRANCH HERE
2178 NO=NOB(I)+MD
VMDN=VAR(I,NO)
IF (M)900,2179,2180
2179 VMDN=(VMDN+VAR(I,NO+1))/2.0
C READ ORIGINAL DATA BACK
2180 CONTINUE
CALL IFILE(20,'DATA')
DO 12001 I=1,NVAR
12001 READ(20)(VAR(I,J),J=1,NOBS)
C CREATE TABLE INDICATING WHETHER ABOVE OR BELOW MEDIAN
C VAR=VMDN COUNTED AS BELOW MEDIAN
C ZERO OUT TABLE FIRST
DO 2210 J=1,2
DO 2210 I=1,NVAR
2210 TABLE(I,J)=0.0
DO 2240 I=1,NVAR
NO=NOB(I)
DO 2240 J=1,NO
IF (VAR(I,J)-VMDN)2220,2220,2230
2220 TABLE(I,1)=TABLE(I,1)+1.0
GO TO 2240
2230 TABLE(I,2)=TABLE(I,2)+1.0
2240 CONTINUE
C CREATE EXPECTED FOR CHI
DO 2260 I=1,NVAR
EX=NOB(I)
2260 EXP(I)=EX/2.0
C PRINT ACTUAL VS. EXPECTED
WRITE (IOUT,2280)
2280 FORMAT(//11X,'INTERMEDIATE RESULTS'//9X,'ABOVE MEDIAN',8X,
112HBELOW MEDIAN,/2X,41HVAR EXPECTED ACTUAL EXPECTED ACTUAL/)
DO 2290 I=1,NVAR
WRITE (IOUT,2285)I,EXP(I),TABLE(I,2),EXP(I),TABLE(I,1)
2285 FORMAT (2H ,I3,2F9.3,2X,2F9.3)
2290 CONTINUE
C CALCULATE CHI SQUARE
C IF NVAR=2, USE SPECIAL FORMULA
VARN=NOB(1)+NOB(2)
IF (NVAR-2)900,2300,2350
2300 CHISQ=(VARN*(ABS(TABLE(1,1)*TABLE(2,2)-TABLE(1,2)*TABLE(2,1))-VARN
1/2.0)**2)/((TABLE(1,1)+TABLE(1,2))*(TABLE(2,1)+TABLE(2,2))*(TABLE(
11,1)+TABLE(2,1))*(TABLE(1,2)+TABLE(2,2)))
MDF=1
GO TO 2370
C CHI SQUARE FOR NVAR GREATER THAN 2
2350 CHISQ=0.0
VARN=0.0
DO 2360 I=1,NVAR
V=NOB(I)
VARN=VARN+V
2360 CHISQ=CHISQ+((TABLE(I,1)-EXP(I))**2)/EXP(I)+((TABLE(I,2)-EXP(I))**
12)/EXP(I)
MDF=NVAR-1
C PRINT CHI SQUARE
2370 C=SQRT(CHISQ/(CHISQ+VARN))
PHISQ=CHISQ/VARN
EL=2.
PHIPR=SQRT(PHISQ/(EL-1.0))
5466 CALL CHIPRB(CHISQ,MDF,PCHI)
WRITE (IOUT,5470)CHISQ,PCHI,MDF,C,PHISQ,PHIPR
5470 FORMAT(' CHI SQUARE = ',F10.4/,' CHI SQUARE PROB. =',F8.5,
1/' DEGREES OF FREEDOM = ',
1I3/' CONTINGENCY COEFFICIENT = ',F8.3/' PHI SQUARE = ',
2F8.3/' PHI PRIME = ',F8.3)
400 RETURN
900 WRITE (IOUT,999)
999 FORMAT (11H IMPOSSIBLE)
10000 RETURN
END
C---------------IOUT, NVAR, NOB ARE INPUT THRU COMMON
SUBROUTINE MNWHU
DIMENSION NOB(20),VAR(20,200)
DIMENSION MVAR(2,200),VA(2,200),RVAR(2,200)
COMMON VAR,NVAR,NOBS,NOB,NCODE,INP,NOUTD,IOUT
C MANN-WHITNEY U FOR ALL PAIRS OF VARIABLES
C PRINT TITLE
WRITE (IOUT,2000)
2000 FORMAT(10X,46('*'))
WRITE(IOUT,3000)
3000 FORMAT(/23X,'MANN-WHITNEY U TEST'/)
WRITE (IOUT,2000)
NV=NVAR-1
DO 3460 KI=1,NV
KFP=KI+1
DO 3460 KS=KFP,NVAR
NI=NOB(KI)
NS=NOB(KS)
C CREARE DUMMY ARRAY MVAR WITH ORIGINAL PLACE VALUES IN IT
C CREATE ARRAY VA CONTAINING TWO VARIABLES IN USE
C STORE RANKS IN RVAR
MM=1
II=KI
NO=NI
RANK=1.0
DO 3080 I=1,2
DO 3070 J=1,NO
MVAR(I,J)=MM
VA(I,J)=VAR(II,J)
RVAR(I,J)=RANK
RANK=RANK+1.0
3070 MM=MM+1
II=KS
3080 NO=NS
C ORDER DATA AND DUMMY PLACE HOLDER
NO=NI
DO 3170 I=1,2
DO 3165 J=1,NO
DO 3155 M=I,2
GO TO (3081,3082),M
3081 NM=NI
GO TO 3083
3082 NM=NS
3083 IF (I-M)3090,3100,999
3090 K=1
GO TO 3120
3100 K=J
3120 DO 3155 N=K,NM
IF (VA(I,J)-VA(M,N))3155,3155,3150
3150 DUM=VA(I,J)
MDUM=MVAR(I,J)
VA(I,J)=VA(M,N)
MVAR(I,J)=MVAR(M,N)
VA(M,N)=DUM
MVAR(M,N)=MDUM
3155 CONTINUE
3165 CONTINUE
3170 NO=NS
C ALTER RANKS FOR TIES
IEQ=1
I=1
M=1
J=0
NO=NI
TI=1.0
3175 J=J+1
IF (J-NO)3190,3180,3300
3180 GO TO (3185,3240),I
3185 M=2
N=1
GO TO 3200
3310 NO=NS
I=2
J=1
3190 N=J+1
3200 IF (VA(I,J)-VA(M,N))3240,3210,999
C CHECK IF TIE OCCURRED
3210 GO TO (3220,3230),IEQ
3220 IEQ=2
AVE=RVAR(I,J)
II=I
JJ=J
C ADD TO TIE COUNT
3230 TI=TI+1.0
AVE=AVE+RVAR(M,N)
GO TO 3175
3240 GO TO (3175,3250),IEQ
3250 AVE=AVE/TI
ITT=TI
IT=0.0
TI=1.0
IEQ=1
3260 RVAR(II,JJ)=AVE
IT=IT+1
IF (IT-ITT)3265,3175,999
3265 JJ=JJ+1
IF (JJ-NI)3260,3260,3270
3270 GO TO (3280,3260),II
3280 II=2
JJ=1
GO TO 3260
3300 GO TO (3310,3350),I
C PUT PROPER RANKS FROM RVAR IN PLACE OF VALUES IN VA
3350 NO=NI
DO 3385 I=1,2
DO 3380 J=1,NO
JJ=MVAR(I,J)-NI
IF (JJ)3360,3360,3370
3360 M=1
N=MVAR(I,J)
GO TO 3380
3370 M=2
N=JJ
3380 VA(M,N)=RVAR(I,J)
3385 NO=NS
C CALCULATE U, MEAN, AND SD
3400 FNI=NI
FNS=NS
R1=0.0
R2=0.0
DO 3410 J=1,NS
3410 R2=R2+VA(2,J)
DO 3420 J=1,NI
3420 R1=R1+VA(1,J)
U=FNI*FNS+(FNI*(FNI+1.0))/2.0-R1
UP=FNI*FNS-U
WRITE (IOUT,3452)KI,KS
3452 FORMAT(//' ','VAR(',I2,') VS. VAR(',I2,')')
3430 FMN=(FNI*FNS)/2.0
SD=SQRT(FMN*(FNI+FNS+1.0)/6.0)
3450 Z=(U-FMN)/SD
C PRINT RESULTS
CALL NORMCV(Z,ZPRB)
WRITE(IOUT,3455)U,UP,FMN,SD,Z,ZPRB
3455 FORMAT(5X,'U = ',F8.3,10X,'U PRIME = ',F8.3,/5X,'MEAN = ',
1F8.3,10X,'ST. DEV. = ',F7.2/,5X,'Z = ',F8.3,
#10X,'Z-PROB. =',F8.5)
3460 CONTINUE
400 RETURN
999 WRITE (IOUT,900)
900 FORMAT (11H IMPOSSIBLE)
10000 RETURN
END
C---------------IOUT, NVAR, NOB, VAR ARE INPUT THRU COMMON
SUBROUTINE WRTST
DIMENSION NOB(20),PLACE(200),VA(200),VAR(20,200)
COMMON VAR,NVAR,NOBS,NOB,NCODE,INP,NOUTD,IOUT
C **** **** WALD-WOLFOWITZ RUNS TEST **** ****
C PRINT TITLE
WRITE (IOUT,7050)
7050 FORMAT(10X,46('*'))
WRITE (IOUT,7060)
7060 FORMAT(/21X,'WALD-WOLFOWITZ RUNS TEST'/)
WRITE (IOUT,7050)
C SET UP TO GO THROUGH EACH PAIR OF VARIABLES
NVAR1=NVAR-1
DO 7900 I=1,NVAR1
I1=I+1
DO 7900 II=I1,NVAR
C TRANSFER SIGNIFICANT VARIABLES TO VA
NO1=NOB(I)
DO 7120 J=1,NO1
PLACE(J)=1.0
7120 VA(J)=VAR(I,J)
NO2=NOB(II)
DO 7130 J=1,NO2
K=NO1+J
PLACE(K)=-1.0
7130 VA(K)=VAR(II,J)
C RANK VA AND REARRANGE PLACE
NNS=NO1+NO2
NNSS=NNS-1
DO 7240 K=1,NNSS
K1=K+1
DO 7240 KK=K1,NNS
IF (VA(K)-VA(KK))7240,7240,7210
7210 DUM=VA(K)
DUMPL=PLACE(K)
VA(K)=VA(KK)
PLACE(K)=PLACE(KK)
VA(KK)=DUM
PLACE(KK)=DUMPL
7240 CONTINUE
C COUNT RUNS
R=1.0
DO 7256 K=1,NNSS
IF (PLACE(K)-PLACE(K+1))7250,7256,7250
7250 R=R+1.0
7256 CONTINUE
C ALLOW FOR TIES
C IEQ=1 MEANS NO TIE IEQ=2 MEANS TIE OCCURRED
IEQ=1
DO 7380 K=1,NNSS
IF (VA(K)-VA(K+1))7300,7260,900
7260 GO TO (7270,7280),IEQ
7270 K1=K
C K1 IS THE FIRST OF A TIED SET
IEQ=2
RUN=1.0
IF (PLACE(K))7273,900,7275
7273 LS=1
LI=0
GO TO 7280
7275 LI=1
LS=0
7280 IF (PLACE(K)-PLACE(K+1))7290,7380,7290
7290 RUN=RUN+1.0
IF (PLACE(K+1))7293,900,7295
7293 LS=LS+1
GO TO 7380
7295 LI=LI+1
GO TO 7380
C CHECK IF TIE OCCURRED
7300 GO TO (7380,7310),IEQ
7310 K2=K+1
IEQ=1
IF (LI-LS)7320,7320,7330
7320 RU=2*LI+1
GO TO 7340
7330 RU=2*LS+1
7340 IF (K1-1)900,7360,7350
7350 IF (PLACE(K)-PLACE(K-1))7360,7362,7360
7360 RU=RU+1.0
7362 IF (K2-NNSS)7365,7365,7367
7365 IF (PLACE(K2)-PLACE(K2+1))7367,7370,7367
7367 RU=RU+1.0
7370 R=R-RUN+RU
7380 CONTINUE
L=R
O1=NO1
O2=NO2
WRITE (IOUT,7385)I,II,NO1,NO2
7385 FORMAT(//' ','VAR(',I2,') VS. VAR(',I2,')'/12X,'N1 = ',
1I3,3X,'N2 = ',I3)
C TEST SIZE OF SAMPLE
IF (NO1-20)7390,7390,7400
7390 IF (NO2-20)7500,7500,7400
C EVALUATION OF SAMPLE WITH N1 OR N2 GREATER THAN 20
7400 FMEAN=((2.0*O1*O2)/(O1+O2))+1.0
SD=SQRT((2.0*O1*O2)*(2.0*O1*O2-O1-O2)/((O1+O2)*(O1+O2)*(O1+O2-1.0)
1))
Z=(ABS(R-FMEAN)-.5)/SD
C PRINT VARIABLE AND N1,N2
CALL NORMCV(Z,ZPRB)
WRITE (IOUT,7450)L,FMEAN,SD,Z,ZPRB
7450 FORMAT(12X,'NUMBER OF RUNS = ',I3/12X,'MEAN = ',F8.3/
#12X,'STANDARD DEVIATION = ',F8.3/12X,'Z = ',F10.6,
#12X,'Z-PROB. =',F8.5)
GO TO 7900
7500 WRITE (IOUT,7550)L
7550 FORMAT(12X,'NUMBER OF RUNS = ',I3)
7900 CONTINUE
400 RETURN
900 WRITE (IOUT,999)
999 FORMAT (11H IMPOSSIBLE)
10000 RETURN
END
C---------------IOUT, NOUTD, INP, NOB, VAR, NVAR ARE INPUT THRU COMMON.
C---------------NCODE IS RETURNED THRU COMMON.
SUBROUTINE MSTST
DIMENSION NOB(20),VA(400)
DIMENSION MVAR(400),VAR(20,200)
COMMON VAR,NVAR,NOBS,NOB,NCODE,INP,NOUTD,IOUT
C THE MOSES TEST OF EXTREME REACTIONS
C PRINT TITLE
WRITE (IOUT,2000)
2000 FORMAT(10X,46('*'))
WRITE (IOUT,100)
100 FORMAT(/18X,'MOSES TEST FOR EXTREME REACTIONS'/)
WRITE (IOUT,2000)
C DO TEST FOR EACH PAIR OF VARIABLES
WRITE(NOUTD,103)
103 FORMAT(' ','WHICH VARIABLE IS THE CONTROL VARIABLE?'/)
READ(INP,102)KI
WRITE(NOUTD,101)
101 FORMAT(' ','ENTER NUMBER OF CONTROL SCORES TO BE',
1' SUBTRACTED FROM EACH EXTREME'/)
READ(INP,102)NCODE
102 FORMAT(I)
NI=NOB(KI)
NC=NI
DO 800 KS=1,NVAR
IF(KS-KI)104,800,104
104 NS=NOB(KS)
NE=NS
C READ VARIABLE VALUES INTO VA AND INDICATOR 1 OR 2 INTO MVAR
DO 150 J=1,NI
VA(J)=VAR(KI,J)
150 MVAR(J)=1
C
DO 160 J=1,NS
K=NI+J
VA(K)=VAR(KS,J)
160 MVAR(K)=2
C ORDER VA AND MVAR IN ASENDING ORDER
NIS=NI+NS
NISO=NIS-1
DO 200 K=1,NISO
J=K+1
DO 200 I=J,NIS
IF (VA(K)-VA(I))200,200,170
170 VDUM=VA(K)
MDUM=MVAR(K)
VA(K)=VA(I)
MVAR(K)=MVAR(I)
VA(I)=VDUM
MVAR(I)=MDUM
200 CONTINUE
C DETERMINE WHICH IS EXTREME
C NCODE IS THE VALUE OF H, COUNT H+1 OF KON FROM EACH END AND
C DENOTE THESE IC AND LC
230 K=0
DO 300 I=1,NIS
IF(MVAR(I)-1)300,250,300
250 K=K+1
IF (K-NCODE)300,300,310
300 CONTINUE
302 WRITE (IOUT,305)
305 FORMAT (57H0 INSUFFICIENT NUMBER OF CONTROL VARIABLES TO JUSTIFY
1 H)
GO TO 800
310 LI=I
K=0
DO 350 J=1,NIS
I=NIS-J+1
IF(MVAR(I)-1)350,320,350
320 K=K+1
IF (K-NCODE)350,350,360
350 CONTINUE
GO TO 302
360 LC=I
C CALCUCATE SH AND G
NSH=LC-LI+1
NG=NSH-NC+2*NCODE
C CALCULATE PROBABILITY USING SUBROUTINE BINO
NTWO=2*NCODE
J=NE+NTWO+1
K=NE
P=BINO(J,K)
DO 400 I=1,NG
J=I+NC-NTWO-2
A=BINO(J,I)
J=NE+NTWO+1-I
K=NE-I
B=BINO(J,K)
400 P=P+A*B
C=BINO(NIS,NC)
P=P/C
C PRINT RESULTS
WRITE (IOUT,410)KI,KS,NE,NC,NCODE,NSH,P
410 FORMAT(' ','CONTROL VAR(',I2,') VS. VAR(',
1I2,')'/' EXTREME VARIABLE HAS',
1I4,10H VARIABLES/21H CONTROL VARIABLE HAS,I4,10H VARIABLES/ 23H SP
1AN ALLOWING FOR H = ,I2,4H IS ,I3/15H PROBABILITY IS ,F15.6)
800 CONTINUE
10000 RETURN
END
C---------------IOUT, NOBS, NVAR, VAR ARE INPUT THRU COMMON
SUBROUTINE COCHQ
DIMENSION NOB(20),G(20)
DIMENSION VAR(20,200)
COMMON VAR,NVAR,NOBS,NOB,NCODE,INP,NOUTD,IOUT
C ********* COCHRAN Q TEST *********
WRITE (IOUT,8010)
8010 FORMAT(10X,46('*'))
WRITE (IOUT,8020)
8020 FORMAT(/26X,'COCHRAN Q TEST'/)
WRITE (IOUT,8010)
C SET SUM AREAS TO ZERO
GSQ=0.0
ELSQ=0.0
DO 8090 I=1,20
8090 G(I)=0.0
C FIND G, LSUM AND SUML-SQUARED
8100 DO 8190 J=1,NOBS
EL=0.0
DO 8150 I=1,NVAR
G(I)=G(I)+VAR(I,J)
8150 EL=EL+VAR(I,J)
WRITE (IOUT,8160)J,EL
8160 FORMAT(10X,2HL(,I2,4H) = ,F8.2)
8190 ELSQ =ELSQ+EL*EL
C FIND GSUM AND SUMG-SQUARED
GSUM=0.0
DO 8199 I=1,NVAR
GSUM=GSUM+G(I)
GSQ=GSQ+G(I)*G(I)
8199 WRITE (IOUT,8200)I,G(I)
8200 FORMAT(10X,2HG(,I2,4H) = ,F8.2)
C CALCULATE Q
X=NVAR
Q=((X-1.)*(X*GSQ-GSUM*GSUM))/(X*GSUM-ELSQ)
MDF=NVAR-1
WRITE (IOUT,8300)Q,MDF
8300 FORMAT(' ','Q = ',F10.4,5X,'D.F. = ',I3)
400 RETURN
10000 RETURN
END
C---------------NCODE IS RETURNED THRU COMMON. IOUT, NVAR,
C--------------- NOBS, VAR ARE INPUT THRU COMMON.
SUBROUTINE FRANV
DIMENSION NOB(20),CDUM(20)
DIMENSION SUM(20),VDUM(20)
DIMENSION KDUM(20),VAR(20,200)
COMMON VAR,NVAR,NOBS,NOB,NCODE,INP,NOUTD,IOUT
NCODE=0
C FRIEDMAN TWO-WAY ANALYSIS OF VARIANCE
C PRINT TITLE
WRITE (IOUT,4010)
4010 FORMAT(10X,46('*')/)
WRITE (IOUT,4020)
4020 FORMAT(14X,'FRIEDMAN TWO-WAY ANALYSIS OF VARIANCE'/)
WRITE (IOUT,4010)
4060 DO 4064 I=1,NVAR
4064 SUM(I)=0.0
C IF(NCODE) = 1 THE DATA IS ALREADY RANKED
IF (NCODE)900,4200,4100
C CHECK FOR NON-POSITIVE RANKING
4100 DO 4150 I=1,NVAR
DO 4150 J=1,NOBS
IF (VAR(I,J))4155,4155,4150
4150 CONTINUE
C CREATE SUM FOR PREVIOUSLY RANKED DATA
4400 DO 4450 I=1,NVAR
DO 4450 J=1,NOBS
4450 SUM(I)=SUM(I)+VAR(I,J)
GO TO 4500
4155 WRITE (IOUT,4160)I,J,VAR(I,J)
4160 FORMAT (11HJ VARIABLE,I3,15H IN OBSERVATION,I4,14H HAS A RANK OF,
1F6.5)
GO TO 400
4200 DO 4350 J=1,NOBS
DO 4210 I=1,NVAR
VDUM(I)=VAR(I,J)
KDUM(I)=I
4210 CDUM(I)=I
C RANK DATA
C RANK DUMMY LINE
N1=NVAR-1
DO 4240 I=1,N1
I1=I+1
DO 4240 II=I1,NVAR
IF (VDUM(I)-VDUM(II))4240,4240,4230
4230 DUM=VDUM(I)
KUM=KDUM(I)
VDUM(I)=VDUM(II)
KDUM(I)=KDUM(II)
VDUM(II)=DUM
KDUM(II)=KUM
4240 CONTINUE
C ALLOW FOR TIES
4250 IEQ=1
DO 4275 I=1,N1
IF (VDUM(I)-VDUM(I+1))4270,4255,900
4255 GO TO (4260,4265),IEQ
4260 T=2.0
AVE=CDUM(I)+CDUM(I+1)
IEQ=2
GO TO 4275
4265 T=T+1.0
AVE=AVE+CDUM(I+1)
GO TO 4275
4270 GO TO (4275,4272),IEQ
4272 K=T
AVE=AVE/T
KK=I-K+1
DO 4274 II=KK,I
4274 CDUM(II)=AVE
IEQ=1
4275 CONTINUE
I=NVAR
GO TO (4277,4272),IEQ
C PUT RANKS IN PROPER COLUMNS
4277 DO 4276 I=1,NVAR
K=KDUM(I)
4276 VDUM(K)=CDUM(I)
C ADD RANKS TO SUM
DO 4280 I=1,NVAR
4280 SUM(I)=SUM(I)+VDUM(I)
4350 CONTINUE
C WRITE OUT SUBTOTALS SUM
4500 DO 4560 I=1,NVAR
4560 WRITE (IOUT,4510)I,SUM(I)
4510 FORMAT (10X,22HSUM OF RANKS IN GROUP ,I2,4H IS ,F9.3)
C FINAL CALCULATIONS
SUMSQ=0.0
DO 4580 I=1,NVAR
4580 SUMSQ=SUMSQ+SUM(I)*SUM(I)
TNK3=NOBS*(NVAR+1)*3
TNKK=NOBS*(NVAR+1)*NVAR
XR=12.0/TNKK*SUMSQ-TNK3
IDF=NVAR-1
CALL CHIPRB(XR,IDF,PCHI)
WRITE(IOUT,4600)XR,PCHI,IDF
4600 FORMAT(15X,'CHI SQUARE IS ',F10.4,/,15X,'CHI SQUARE PROB. =',F8.5,
#/10X,'DEGREES OF FREEDOM = ',I3)
400 RETURN
900 WRITE (IOUT,999)
999 FORMAT (11H IMPOSSIBLE)
10000 RETURN
END
C---------------VAR, NVAR, NOBS, NOB, IOUT ARE INPUT THRU COMMON. VAR
C--------------- READ IN ST. 4500.
SUBROUTINE KWANV
DIMENSION NOB(20),VAR(20,200)
DIMENSION VARD(20,200)
COMMON VAR,NVAR,NOBS,NOB,NCODE,INP,NOUTD,IOUT
WRITE (IOUT,2000)
WRITE (IOUT,4000)
4000 FORMAT(/11X,'KRUSKAL-WALLIS ONE-WAY ANALYSIS OF VARIANCE'/)
WRITE (IOUT,2000)
2000 FORMAT(10X,46('*'))
C
C KRUSKAL-WALLIS
C CALCULATION SEGMENT OF PROGRAM
C DATA HAS BEEN SET
C
END FILE 21
END FILE 20
CALL OFILE(20,'DATA')
CALL OFILE(21,'DATE')
DO 4050 I=1,NVAR
4050 WRITE(20)(VAR(I,J),J=1,NOBS)
END FILE 20
C CREATE DUMMY ARRAY VARD WITH ORIGINAL PLACE VALUES IN IT
PLACE=1.0
DO 4070 I=1,NVAR
NO=NOB(I)
DO 4071 J=1,NO
VARD(I,J)=PLACE
4071 PLACE=PLACE+1.0
4070 CONTINUE
C ORDER DATA AND DUMMY PLACE HOLDER
DO 4170 I=1,NVAR
NO=NOB(I)
DO 4170 J=1,NO
DO 4170 M=I,NVAR
NM=NOB(M)
IF (M-I)900,4080,4090
4080 K=J
GO TO 4100
4090 K=1
4100 DO 4170 N=K,NM
IF (VAR(I,J)-VAR(M,N))4170,4170,4160
4160 DUM=VAR(I,J)
DUMD=VARD(I,J)
VAR(I,J)=VAR(M,N)
VARD(I,J)=VARD(M,N)
VAR(M,N)=DUM
VARD(M,N)=DUMD
4170 CONTINUE
C WRITE PLACE VALUES IN VARD ON TAPE
DO 4200 I=1,NVAR
NO=NOB(I)
DO 4200 J=1,NO
4200 WRITE(21)VARD(I,J)
C STORE RANKS IN VARD
PLACE=1.0
DO 4210 I=1,NVAR
NO=NOB(I)
DO 4210 J=1,NO
VARD(I,J)=PLACE
4210 PLACE=PLACE+1.0
C ALTER RANKS FOR TIES
T=0.0
ND=0
TI=1.0
I=0
4220 I=I+1
IF (I-NVAR)4221,4221,4300
4221 J=0
NO=NOB(I)
4222 J=J+1
IF (J-NO)4230,4225,4220
4225 IF (I-NVAR)4226,4260,900
4226 M=I+1
N=1
GO TO 4240
4230 M=I
N=J+1
4240 IF (VAR(I,J)-VAR(M,N))4270,4245,900
C ADD TO TIE COUNT
4245 IF (TI-1.0)900,4250,4255
4250 II=I
JJ=J
AVE=VARD(I,J)
4255 TI=TI+1.0
AVE=AVE+VARD(M,N)
GO TO 4222
4260 ND=1
C CHECK IF TIE OCCURRED
4270 IF (TI-1.0)900,4275,4271
4275 AVE=VARD(I,J)
IF (ND-1)4222,4300,900
4271 T=T+TI**3-TI
ITT=TI
AVE=AVE/TI
DO 4290 K=1,ITT
VARD(II,JJ)=AVE
JJ=JJ+1
IF (JJ-NOB(II))4290,4290,4280
4280 II=II+1
JJ=1
4290 CONTINUE
TI=1.0
GO TO 4222
C BRING PLACE INDICATORS BACK AND PUT RANKS IN PLACE
4300 CALL IFILE(21,'DATE')
DO 4400 M=1,NVAR
NO=NOB(M)
DO 4400 N=1,NO
READ(21)PLACE
I=0
4320 I=I+1
PJ=NOB(I)
PLACE=PLACE-PJ
IF (PLACE)4350,4350,4320
4350 J=PLACE
J=NOB(I)+J
VAR(I,J)=VARD(M,N)
4400 CONTINUE
END FILE 21
C FINAL CALCULATIONS
H=0.0
VN=0.0
DO 4410 I=1,NVAR
HH=0.0
NO=NOB(I)
VNO=NO
VN=VN+VNO
DO 4405 J=1,NO
4405 HH=HH+VAR(I,J)
WRITE (IOUT,4407)I,HH
4407 FORMAT(/10X,'SUM OF RANKS FOR VARIABLE ',I3,' IS ',F8.3)
4410 H=H+(HH*HH)/VNO
PLACE=VN+1.0
H=((12.0/(VN*PLACE))*H-(3.0*PLACE))/(1.0-(T/(VN*VN*VN-VN)))
MDF=NVAR-1
WRITE (IOUT,4420)H,MDF
4420 FORMAT(/10X,'H = ',F8.3,5X,'DEGREES OF FREEDOM = ',I2)
CALL IFILE(20,'DATA')
DO 4500 I=1,NVAR
4500 READ(20)(VAR(I,J),J=1,NOBS)
400 RETURN
900 WRITE (IOUT,999)
999 FORMAT (11H IMPOSSIBLE)
10000 RETURN
END
FUNCTION BINO(K,N)
KN=K-N
BINO=1.
DO 100 I=1,KN
YI=I
YK=K-I+1
Y=YK/YI
100 BINO=BINO*Y
RETURN
END