Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap5_198111 - decus/20-0137/nonpar/nonpar.for
There is 1 other file named nonpar.for in the archive. Click here to see a list.
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