Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap5_198111 - decus/20-0137/kolm/kolm.for
There is 1 other file named kolm.for in the archive. Click here to see a list.
C	WESTERN MICHIGAN UNIVERSITY
C	KOLM.F4 (FILE NAME ON LIBRARY DECTAPE)
C	KOLM, 1.7.1 (CALLING NAME, SUBLST NO.)
C	KOLMOGOROV-SMIRNOV TEST PROG.
C	ADAPTED FROM DECUS(NO.10-101, IBM-SSP) BY
C	 BERENICE GAN HOUCHARD.  MODIFICATION ON STATISTICAL DESIGN
C	 BY DR. MICHAEL STOLINE.
C	LIBRARY DECTAPE PROGS. USED:  USAGE.MAC
C	FORWMU PROGS. USED:  TTYPTY, DEVCHG, EXISS, PRINTS
C	APLIB PROGS. USED:  IO, GETFOR
C	INTERNAL SUBR. USED:  KOLMO, KOLM2, SMIRN, NDTR
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C
C
C	LIMITATIONS:
C
C	(1)  NUMBER OF DATA POINTS FOR EACH SAMPLE IS AT MOST 5001
C	(2)  ONLY 1 OBJECT TIME FORMAT LINE OR CARD
C
C***********************************************************************
C
      DIMENSION ID(16),NOTF(16),X(5001),Y(5001),TIT1(4,2),DIST(4,3)
      DOUBLE PRECISION TIT1
      DATA ((TIT1(I,J),J=1,2),I=1,4)/'NORMAL','      ','EXPONE',
     1 'NTIAL ','CAUCHY','      ','UNIFOR','M     '/
C
C***********************************************************************
C	DEVICES USED:
C
C	IDLG--DEVICE USED TO COMMUNICATE WITH USERS
C	      IT IS ALWAYS SET TO -1
C	ICC---DEVICE USED TO READ USER'S RESPONSES
C	      IT IS ALWAYS SET TO -4
C	INP---DEVICE USED TO READ INPUT DATA
C	      ITS LOGICAL NUMBER IS DETERMINED BY SUBROUTINE IO
C	IOUT--DEVICE USED TO WRITE THE REPORT
C	      ITS LOGICAL NUMBER IS DETERMINED BY SUBROUTINE IO
C***********************************************************************
C
      IDLG=-1
      ICC=-4
      INP=2
      IOUT=3
C
C**********************************************************************
C	CALL SUBROUTINE USAGE AND ADD 1 TO PROGRAM USAGE COUNTER
C***********************************************************************
C
C	CALL USAGE('KOLM')
C
C***********************************************************************
C	DETERMINE IF JOB IS ON TELETYPE OR PSEUDO-TELETYPE
C
C	IF ICODE = 0  JOB IS ON TELETYPE
C		 =-1  JOB IS ON PSEUDO-TELETYPE
C***********************************************************************
C
      CALL TTYPTY(ICODE)
C
C***********************************************************************
C	GATHER ALL I/O INFORMATION AND DETERMINE WHICH FORMAT TO USE
C***********************************************************************
C
C---------------IO3, IO2 ARE RETURNED.	OTHER ARGS. ARE INPUT.
C---------------1 MEANS OUTPUT? PRINTS, 0 - INPUT? PRINTS.
      CALL IO(IOUT,IO3,ICC,IDLG,1,ICODE)
100   CALL IO(INP,IO2,ICC,IDLG,0,ICODE)
      ITYPE=2
C---------------NOTF, ISTD ARE RETURNED.  OTHER ARGS. ARE INPUT.
C---------------16=NO. OF OBJ. TIME FORMAT WORDS (1 LINE).  ITYPE=2
C--------------- MEANS F - TYPE FORMAT ONLY.
      CALL GETFOR(IDLG,ICC,NOTF,ISTD,16,ITYPE)
C
C
      DO 9 I=1,4
      DO 9 J=1,3
9     DIST(I,J)=0
      WRITE(IDLG,10)
10    FORMAT(' ENTER HEADER'/)
      READ(ICC,11) ID
11    FORMAT(16A5)
12    WRITE(IDLG,13)
13    FORMAT(' 1-SAMPLE OR 2-SAMPLE TEST?--',$)
      READ(ICC,14) IS
14    FORMAT(I1)
      IF ((IS.LT.3).AND.(IS.GT.0)) GO TO 20
      WRITE(IDLG,15) IS
15    FORMAT(' ', I1, '-SAMPLE NOT POSSIBLE, TRY AGAIN'/)
      IF (ICODE) 16,12,12
16    CALL EXIT
20    IF (IS.EQ.2) GO TO 50
C
C***********************************************************************
C	GATHER INFORMATION FOR ONE-SAMPLE TEST
C***********************************************************************
C
C
C	DETERMINE WHICH PROBABILITY DENSITY FUNCTION TO USE
C
C	DIST(1,1)--F(X) IS THE NORMAL PDF
C	DIST(2,1)--F(X) IS THE EXPONENTIAL PDF
C	DIST(3,1)--F(X) IS THE CAUCHY PDF
C	DIST(4,1)--F(X) IS THE UNIFORM PDF
C***********************************************************************
C
21    WRITE(IDLG,22)
22    FORMAT(' ENTER PDF OPTION'/)
      READ(ICC,23) (DIST(I,1),I=1,4)
23    FORMAT(4F1.0)
      DO 24 I=1,4
      IF (DIST(I,1).NE.0) GO TO 30
24    CONTINUE
      WRITE(IDLG,25)
25    FORMAT(' ERROR IN PDF OPTION STATEMENT, TRY AGAIN'/)
      IF (ICODE) 16,21,21
30    WRITE(IDLG,300)
300   FORMAT('-PDF PARAMETERS, SEPARATE THE NUMBERS BY COMMA'/)
      DO 31 I=1,4
      IF (DIST(I,1).NE.1) GO TO 31
      IF ((I.EQ.1).OR.(I.EQ.2)) WRITE(IDLG,32) (TIT1(I,J),J=1,2)
32    FORMAT(' ENTER MEAN AND ST. DEV. OF THE ',2A6,'PDF'/)
      IF (I.EQ.3) WRITE(IDLG,33)(TIT1(I,J),J=1,2)
33    FORMAT(' ENTER 1ST QUARTILE AND MEDIAN OF THE ',2A6,'PDF'/)
      IF (I.EQ.4) WRITE(IDLG,34)(TIT1(I,J),J=1,2)
34    FORMAT(' ENTER LEFT AND RIGHT ENDPOINTS OF THE ',2A6,'PDF'/)
      READ(ICC,37)(DIST(I,J),J=2,3)
37    FORMAT(10F)
31    CONTINUE
C
C***********************************************************************
C	ASK AND CHECK ON THE SAMPLE SIZE
C***********************************************************************
C
380   WRITE(IDLG,38)
38    FORMAT(' SAMPLE SIZE =?  ',$)
      READ(ICC,39) N
39    FORMAT(2I)
      IF (N.GT.0) GO TO 40
      WRITE(IDLG,390)
390   FORMAT(' SAMPLE SIZE 0 IMPOSSIBLE, TRY AGAIN'/)
      IF (ICODE)16,380,380
40    IF (N.LE.5001) GO TO 60
43    WRITE(IDLG,44)
44    FORMAT(' SAMPLE SIZE TOO LARGE, MAXIMUM IS 5000. JOB TERMINATE'/)
      CALL EXIT
C
C***********************************************************************
C	GATHER INFORMATION FOR TWO-SAMPLES TEST
C***********************************************************************
C
50    WRITE(IDLG,51)
51    FORMAT(' ENTER SAMPLE SIZES,SEPARATE THEM BY COMMA'/)
      READ(ICC,39) N,M
      IF ((N.GT.0).AND.(M.GT.0)) GO TO 510
      WRITE(IDLG,390)
      IF (ICODE)16,50,50
510   IF ((N.GT.5000).OR.(M.GT.5000)) GO TO 43
52    WRITE(IDLG,53)
53    FORMAT(' ENTER INPUT DATA CODE:'/6X,'1--IF DATA FOR SAMPLE 2
     1 FOLLOW THOSE OF SAMPLE 1'/6X,'2--IF SAMPLE 1 AND 2 ARE ON THE
     2 SAME LINE'/)
      READ(ICC,39) LOR
      IF ((LOR.EQ.1).OR. (LOR.EQ.2)) GO TO 60
      WRITE(IDLG,54) LOR
54    FORMAT(' ERROR IN INPUT DATA CODE ',I2,' TRY AGAIN'/)
      IF (ICODE) 16, 52,52
C
C***********************************************************************
C	ADJUST FORMAT IF NECESSARY
C***********************************************************************
C
60    IF (ISTD.NE.1) GO TO 62
      NOTF(1)='(10F)'
      DO 61 I=2,16
61    NOTF(I)=' '
62    IF (IO2.NE.'TTY') GO TO 64
      WRITE(IDLG,63)
63    FORMAT(' ENTER DATA'/)
      GO TO 70
64    WRITE(IDLG,65)
65    FORMAT(' PLEASE WAIT, YOUR DATA IS BEING PROCESSED'/)
70    IF (IS-1) 71,71,90
C
C***********************************************************************
C	FOR ONE-SAMPLE TEST ONLY
C***********************************************************************
C
71    IF ((IO2.EQ.'TTY').AND. (ISTD.EQ.1)) WRITE(IDLG,72)
72    FORMAT('+AT MOST 10 NUMBERS PER LINE SEPARATED BY COMMAS'/)
      READ(INP,NOTF)(X(I),I=1,N)
      WRITE(IOUT,73) ID, N
73    FORMAT(1H1,16A5/'-KOLMOGOROV-SMIRNOV ONE-SAMPLE TEST'/' SAMPLE
     1 SIZE =', I5//)
      IF (N.LT.100) WRITE(IOUT,730)
730   FORMAT('-NOTE THE REMARKS CONCERNING ASYMPTOTIC RESULTS AND
     1 SAMPLE SIZE'/' IN WMU COMPUTER CENTER PROGRAM DOCUMENTATION
     2 #1.7.1 SECTION 1.0.'/)
      IES=0
      DO 74 I=1,4
      IF (DIST(I,1).LE.0) GO TO 74
      CALL KOLMO(X,N,Z,P,I,DIST(I,2),DIST(I,3),IER,D1)
      IES=IER+IES
      IF (IER.NE.0) GO TO 74
      WRITE(IOUT,75)(TIT1(I,J),J=1,2)
75    FORMAT('-THE HYPOTHESIS THAT THE SAMPLE IS FROM A(N) ', 2A6,
     1 'DISTRIBUTION')
      IF (I-3) 76,77,78
76    S2=DIST(I,3)*DIST(I,3)
      WRITE(IOUT, 760) DIST(I,2),S2
760   FORMAT(' WITH MEAN', F13.2, ' AND VARIANCE', F13.2)
      GO TO 79
77    S2=DIST(I,3)
      WRITE(IOUT, 770) DIST(I,2),S2
770   FORMAT(' WITH FIRST QUARTILE', F13.2, ' AND MEDIAN ', F13.2)
      GO TO 79
78    WRITE(IOUT, 780) DIST(I,2),DIST(I,3)
780   FORMAT(' IN THE INTERVAL', F13.2, ' TO', F13.2, ' INCLUSIVE')
79    WRITE(IOUT, 790) P,D1
790   FORMAT(' CAN BE REJECTED WITH PROBABILITY ',F6.3,' OF BEING
     1 INCORRECT.'/' THE STATISTIC D1 IS', E12.4,' FOR THIS SAMPLE.')
74    CONTINUE
C
C***********************************************************************
C	END OF	ONE-SAMPLE TEST, WRITE OUT SORTED DATA IF REQUESTED
C***********************************************************************
C
C
C***********************************************************************
C	ERROR IN PDF PARAMETER, WRITE OUT A MESSAGE
C***********************************************************************
C
81    IF (IES.EQ.0) GO TO 990
      WRITE(IOUT, 82)
82    FORMAT('-AT LEAST ONE (S) ENTRY PARAMTER FOR THE ONE-SAMPLE
     1 TEST WAS INCORRECT.'/' THE TEST FOR THE ASSOCIATED
     2 CONTINUOUS PDF WAS IGNORED.'/)
      GO TO 990
C
C***********************************************************************
C***********************************************************************
C	FOR TWO-SAMPLES TEST ONLY
C***********************************************************************
C
90    IF (LOR.EQ.2) GO TO 93
91    IF (IO2.NE.'TTY') GO TO 920
      WRITE(IDLG,92)
92    FORMAT('+SAMPLE 1 FIRST FOLLOWED BY SAMPLE 2'/)
      IF (ISTD.EQ.1) WRITE(IDLG,72)
920   READ(INP, NOTF)(X(I),I=1,N)
      READ(INP, NOTF)(Y(I),I=1,M)
      GO TO 96
93    IF (IO2.NE.'TTY') GO TO 941
      WRITE(IDLG,94)
94    FORMAT('+SAMPLE 1 AND 2 ON THE SAME LINE'/)
      IF (ISTD.EQ.1) WRITE(IDLG,940)
940   FORMAT('+SEPARATE THE TWO NUMBERS BY A COMMA'/)
941   LAST=N
      IF (M.GT.N) LAST=M
      DO 95 I=1,LAST
95    READ(INP,NOTF)X(I),Y(I)
96    WRITE(IOUT,960) ID,N,M
960   FORMAT(1H1,16A5/'-KOLMOGOROV-SMIRNOV TWO-SAMPLE TEST'/' SAMPLE
     1 SIZES =',I5,',',I5//)
      IF ((N.LT.100).OR.(M.LT.100)) WRITE(IOUT,730)
      CALL KOLM2(X,Y,N,M,Z,P,D2)
C
C***********************************************************************
C	END OF TWO-SAMPLES TESTS, WRITE OUT SORTED DATA IF REQUESTED
C***********************************************************************
C
98    WRITE(IOUT, 99) P,D2
99    FORMAT('-THE HYPOTHESIS THAT THE TWO SAMPLES ARE FROM THE SAME
     1 POPULATION CAN'/' BE REJECTED WITH (ASYMPTOTIC) PROBABILITY',
     2 F6.3, ' OF BEING INCORRECT.'/' THE STATISTIC D2 IS ',E12.4, ' FOR
     3 THESE SAMPLES.')
C
C***********************************************************************
C	END OF ONE JOB, BRANCH AND DETERMINE IF MORE IS TO BE DONE
C***********************************************************************
C
990   WRITE(IDLG,991)
991   FORMAT('-END OF TEST'/)
      GO TO 100
      END
C
C     ..................................................................
C
C	 SUBROUTINE KOLMO
C
C	 PURPOSE
C	    TESTS THE DIFFERENCE BETWEEN EMPIRICAL AND THEORETICAL
C	    DISTRIBUTIONS  USING THE KOLMOGOROV-SMIRNOV TEST
C
C	 USAGE
C	    CALL KOLMO(X,N,Z,PROB,IFCOD,U,S,IER,D1)
C
C	 DESCRIPTION OF PARAMETERS
C	    X	 - INPUT VECTOR OF N INDEPENDENT OBSERVATIONS.	ON
C		   RETURN FROM KOLMO, X HAS BEEN SORTED INTO A
C		   MONOTONIC NON-DECREASING SEQUENCE.
C	    N	 - NUMBER OF OBSERVATIONS IN X
C	    Z	 - OUTPUT VARIABLE CONTAINING THE GREATEST VALUE WITH
C		   RESPECT TO X OF  SQRT(N)*ABS(FN(X)-F(X)) WHERE
C		   F(X) IS A  THEORETICAL DISTRIBUTION FUNCTION AND
C		   FN(X) AN EMPIRICAL DISTRIBUTION FUNCTION.
C	    PROB - OUTPUT VARIABLE CONTAINING THE PROBABILITY OF
C		   THE STATISTIC BEING GREATER THAN OR EQUAL TO Z IF
C		   THE HYPOTHESIS THAT X IS FROM THE DENSITY UNDER
C		   CONSIDERATION IS TRUE.  E.G., PROB = 0.05 IMPLIES
C		   THAT ONE CAN REJECT THE NULL HYPOTHESIS THAT THE SET
C		   X IS FROM THE DENSITY UNDER CONSIDERATION WITH 5 PER
C		   CENT PROBABILITY OF BEING INCORRECT.  PROB = 1. -
C		   SMIRN(Z).
C	    IFCOD- A CODE DENOTING THE PARTICULAR THEORETICAL
C		   PROBABILITY DISTRIBUTION FUNCTION BEING CONSIDERED.
C		   = 1---F(X) IS THE NORMAL PDF.
C		   = 2---F(X) IS THE EXPONENTIAL PDF.
C		   = 3---F(X) IS THE CAUCHY PDF.
C		   = 4---F(X) IS THE UNIFORM PDF.
C		   = 5---F(X) IS USER SUPPLIED.
C	    U	 - WHEN IFCOD IS 1 OR 2, U IS THE MEAN OF THE DENSITY
C		   GIVEN ABOVE.
C		   WHEN IFCOD IS 3, U IS THE FIRST QUARTILE OF THE
C		   CAUCHY DENSITY.
C		   WHEN IFCOD IS 4, U IS THE LEFT ENDPOINT OF THE
C		   UNIFORM DENSITY.
C		   WHEN IFCOD IS 5, U IS USER SPECIFIED.
C	    S	 - WHEN IFCOD IS 1 OR 2, S IS THE STANDARD DEVIATION OF
C		   DENSITY GIVEN ABOVE, AND SHOULD BE POSITIVE.
C		   WHEN IFCOD IS 3, S SPECIFIES THE MEDIAN OF THE
C		   CAUCHY DISTRIBUTION.  S SHOULD BE NON-ZERO.
C		   IF IFCOD IS 4, S IS THE RIGHT ENDPOINT OF THE UNIFORM
C		   DENSITY.  S SHOULD BE GREATER THAN U.
C		   IF IFCOD IS 5, S IS USER SPECIFIED.
C	    IER  - ERROR INDICATOR WHICH IS NON-ZERO IF S VIOLATES ABOVE
C		   CONVENTIONS.  ON RETURN NO TEST HAS BEEN MADE, AND X
C		   AND Y HAVE BEEN SORTED INTO MONOTONIC NON-DECREASING
C		   SEQUENCES.  IER IS SET TO ZERO ON ENTRY TO KOLMO.
C		   IER IS CURRENTLY SET TO ONE IF THE USER-SUPPLIED PDF
C		   IS REQUESTED FOR TESTING.  THIS SHOULD BE CHANGED
C		   (SEE REMARKS) WHEN SOME PDF IS SUPPLIED BY THE USER.
C
C	 REMARKS
C	    N SHOULD BE GREATER THAN OR EQUAL TO 100.  (SEE THE
C	    MATHEMATICAL DESCRIPTION GIVEN FOR THE PROGRAM SMIRN,
C	    CONCERNING ASYMPTOTIC FORMULAE)  ALSO, PROBABILITY LEVELS
C	    DETERMINED BY THIS PROGRAM WILL NOT BE CORRECT IF THE
C	    SAME SAMPLES ARE USED TO ESTIMATE PARAMETERS FOR THE
C	    CONTINUOUS DISTRIBUTIONS WHICH ARE USED IN THIS TEST.
C	    (SEE THE MATHEMATICAL DESCRIPTION FOR THIS PROGRAM)
C	    F(X) SHOULD BE A CONTINUOUS FUNCTION.
C	    ANY USER SUPPLIED CUMULATIVE PROBABILITY DISTRIBUTION
C	    FUNCTION SHOULD BE CODED BEGINNING WITH STATEMENT 26 BELOW,
C	    AND SHOULD RETURN TO STATEMENT 27.
C
C	    DOUBLE PRECISION USAGE---IT IS DOUBTFUL THAT THE USER WILL
C	    WISH TO PERFORM THIS TEST USING DOUBLE PRECISION ACCURACY.
C	    IF ONE WISHES TO COMMUNICATE WITH KOLMO IN A DOUBLE
C	    PRECISION PROGRAM, HE SHOULD CALL THE FORTRAN SUPPLIED
C	    PROGRAM SNGL(X) PRIOR TO CALLING KOLMO, AND CALL THE
C	    FORTRAN SUPPLIED PROGRAM DBLE(X) AFTER EXITING FROM KOLMO.
C	    (NOTE THAT SUBROUTINE SMIRN DOES HAVE DOUBLE PRECISION
C	    CAPABILITY AS SUPPLIED BY THIS PACKAGE.)
C
C
C	 SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C	    SMIRN, NDTR, AND ANY USER SUPPLIED SUBROUTINES REQUIRED.
C
C	 METHOD
C	    FOR REFERENCE, SEE (1) W. FELLER--ON THE KOLMOGOROV-SMIRNOV
C	    LIMIT THEOREMS FOR EMPIRICAL DISTRIBUTIONS--
C	    ANNALS OF MATH. STAT., 19, 1948.  177-189,
C	    (2) N. SMIRNOV--TABLE FOR ESTIMATING THE GOODNESS OF FIT
C	    OF EMPIRICAL DISTRIBUTIONS--ANNALS OF MATH. STAT., 19,
C	    1948.  279-281.
C	    (3) R. VON MISES--MATHEMATICAL THEORY OF PROBABILITY AND
C	    STATISTICS--ACADEMIC PRESS, NEW YORK, 1964.  490-493,
C	    (4) B.V. GNEDENKO--THE THEORY OF PROBABILITY--CHELSEA
C	    PUBLISHING COMPANY, NEW YORK, 1962.  384-401.
C
C     ..................................................................
C
C---------------Z, PROB, IER, D1 ARE RETURNED.	OTHER ARGS. ARE INPUT.
      SUBROUTINE KOLMO(X,N,Z,PROB,IFCOD,U,S,IER,D1)
      DIMENSION X(1)
C
C	   NON DECREASING ORDERING OF X(I)'S  (DUBY METHOD)
C
      IER=0
      DO 5 I=2,N
      IF(X(I)-X(I-1))1,5,5
    1 TEMP=X(I)
      IM=I-1
      DO 3 J=1,IM
      L=I-J
      IF(TEMP-X(L))2,4,4
    2 X(L+1)=X(L)
    3 CONTINUE
      X(1)=TEMP
      GO TO 5
    4 X(L+1)=TEMP
    5 CONTINUE
C
C	    COMPUTES MAXIMUM DEVIATION DN IN ABSOLUTE VALUE BETWEEN
C	    EMPIRICAL AND THEORETICAL DISTRIBUTIONS
C
      NM1=N-1
      XN=N
      DN=0.0
      FS=0.0
      IL=1
    6 DO 7  I=IL,NM1
      J=I
      IF(X(J)-X(J+1))9,7,9
    7 CONTINUE
    8 J=N
    9 IL=J+1
      FI=FS
      FS=FLOAT(J)/XN
      IF(IFCOD-2)10,13,17
   10 IF(S)11,11,12
   11 IER=1
      GO TO 29
   12 Z =(X(J)-U)/S
      CALL NDTR(Z,Y,D)
      GO TO 27
   13 IF(S)11,11,14
   14 Z=(X(J)-U)/S+1.0
      IF(Z)15,15,16
   15 Y=0.0
      GO TO 27
   16 Y=1.-EXP(-Z)
      GO TO 27
   17 IF(IFCOD-4)18,20,26
   18 IF(S-U)19,11,19
   19 Y=ATAN((X(J)-S)/(S-U))*0.3183099+0.5
      GO TO 27
   20 IF(S-U)11,11,21
   21 IF(X(J)-U)22,22,23
   22 Y=0.0
      GO TO 27
   23 IF(X(J)-S)25,25,24
   24 Y=1.0
      GO TO 27
   25 Y=(X(J)-U)/(S-U)
      GO TO 27
   26 IER=1
      GO TO 29
   27 EI=ABS(Y-FI)
      ES=ABS(Y-FS)
      DN=AMAX1(DN,EI,ES)
      IF(IL-N)6,8,28
C
C	    COMPUTES Z=DN*SQRT(N)  AND	PROBABILITY
C
   28 Z=DN*SQRT(XN)
	D1=DN
      CALL SMIRN(Z,PROB)
      PROB=1.0-PROB
   29 RETURN
      END
C
C     ..................................................................
C
C	 SUBROUTINE KOLM2
C
C	 PURPOSE
C
C	    TESTS THE DIFFERENCE BETWEEN TWO SAMPLE DISTRIBUTION
C	    FUNCTIONS USING THE KOLMOGOROV-SMIRNOV TEST
C
C	 USAGE
C	    CALL KOLM2(X,Y,N,M,Z,PROB,D2)
C
C	 DESCRIPTION OF PARAMETERS
C	    X	 - INPUT VECTOR OF N INDEPENDENT OBSERVATIONS.	ON
C		   RETURN FROM KOLM2, X HAS BEEN SORTED INTO A
C		   MONOTONIC NON-DECREASING SEQUENCE.
C	    Y	 - INPUT VECTOR OF M INDEPENDENT OBSERVATIONS.	ON
C		   RETURN FROM KOLM2, Y HAS BEEN SORTED INTO A
C		   MONOTONIC NON-DECREASING SEQUENCE.
C	    N	 - NUMBER OF OBSERVATIONS IN X
C	    M	 - NUMBER OF OBSERVATIONS IN Y
C	    Z	 - OUTPUT VARIABLE CONTAINING THE GREATEST VALUE WITH
C		   RESPECT TO THE SPECTRUM OF X AND Y OF
C		   SQRT((M*N)/(M+N))*ABS(FN(X)-GM(Y)) WHERE
C		   FN(X) IS THE EMPIRICAL DISTRIBUTION FUNCTION OF THE
C		   SET (X) AND GM(Y) IS THE EMPIRICAL DISTRIBUTION
C		   FUNCTION OF THE SET (Y).
C	    PROB - OUTPUT VARIABLE CONTAINING THE PROBABILITY OF
C		   THE STATISTIC BEING GREATER THAN OR EQUAL TO Z IF
C		   THE HYPOTHESIS THAT X AND Y ARE FROM THE SAME PDF IS
C		   TRUE.  E.G., PROB= 0.05 IMPLIES THAT ONE CAN REJECT
C		   THE NULL HYPOTHESIS THAT THE SETS X AND Y ARE FROM
C		   THE SAME DENSITY WITH 5 PER CENT PROBABILITY OF BEING
C		   INCORRECT.  PROB = 1. - SMIRN(Z).
C
C	 REMARKS
C	    N AND M SHOULD BE GREATER THAN OR EQUAL TO 100.  (SEE THE
C	    MATHEMATICAL DESCRIPTION FOR THIS SUBROUTINE AND FOR THE
C	    SUBROUTINE SMIRN, CONCERNING ASYMPTOTIC FORMULAE).
C
C	    DOUBLE PRECISION USAGE---IT IS DOUBTFUL THAT THE USER WILL
C	    WISH TO PERFORM THIS TEST USING DOUBLE PRECISION ACCURACY.
C	    IF ONE WISHES TO COMMUNICATE WITH KOLM2 IN A DOUBLE
C	    PRECISION PROGRAM, HE SHOULD CALL THE FORTRAN SUPPLIED
C	    PROGRAM SNGL(X) PRIOR TO CALLING KOLM2, AND CALL THE
C	    FORTRAN SUPPLIED PROGRAM DBLE(X) AFTER EXITING FROM KOLM2.
C	    (NOTE THAT SUBROUTINE SMIRN DOES HAVE DOUBLE PRECISION
C	    CAPABILITY AS SUPPLIED BY THIS PACKAGE.)
C
C
C	 SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C	    SMIRN
C
C	 METHOD
C	    FOR REFERENCE, SEE (1) W. FELLER--ON THE KOLMOGOROV-SMIRNOV
C	    LIMIT THEOREMS FOR EMPIRICAL DISTRIBUTIONS--
C	    ANNALS OF MATH. STAT., 19, 1948.  177-189,
C	    (2) N. SMIRNOV--TABLE FOR ESTIMATING THE GOODNESS OF FIT
C	    OF EMPIRICAL DISTRIBUTIONS--ANNALS OF MATH. STAT., 19,
C	    1948.  279-281.
C	    (3) R. VON MISES--MATHEMATICAL THEORY OF PROBABILITY AND
C	    STATISTICS--ACADEMIC PRESS, NEW YORK, 1964.  490-493,
C	    (4) B.V. GNEDENKO--THE THEORY OF PROBABILITY--CHELSEA
C	    PUBLISHING COMPANY, NEW YORK, 1962.  384-401.
C
C     ..................................................................
C
C---------------Z, PROB, D2 ARE RETURNED.  OTHER ARGS. ARE INPUT.
      SUBROUTINE KOLM2(X,Y,N,M,Z,PROB,D2)
      DIMENSION X(1),Y(1)
C
C	 SORT X INTO ASCENDING SEQUENCE
C
      DO 5 I=2,N
      IF(X(I)-X(I-1))1,5,5
    1 TEMP=X(I)
      IM=I-1
      DO 3 J=1,IM
      L=I-J
      IF(TEMP-X(L))2,4,4
    2 X(L+1)=X(L)
    3 CONTINUE
      X(1)=TEMP
      GO TO 5
    4 X(L+1)=TEMP
    5 CONTINUE
C
C	 SORT Y INTO ASCENDING SEQUENCE
C
      DO 10 I=2,M
      IF(Y(I)-Y(I-1))6,10,10
    6 TEMP=Y(I)
      IM=I-1
      DO 8  J=1,IM
      L=I-J
      IF(TEMP-Y(L))7,9,9
    7 Y(L+1)=Y(L)
    8 CONTINUE
      Y(1)=TEMP
      GO TO 10
    9 Y(L+1)=TEMP
   10 CONTINUE
C
C	 CALCULATE D = ABS(FN-GM) OVER THE SPECTRUM OF X AND Y
C
      XN=FLOAT(N)
      XN1=1./XN
      XM=FLOAT(M)
      XM1=1./XM
      D=0.0
      I=0
      J=0
      K=0
      L=0
   11 IF(X(I+1)-Y(J+1))12,13,18
   12 K=1
      GO TO 14
   13 K=0
   14 I=I+1
      IF(I-N)15,21,21
   15 IF(X(I+1)-X(I))14,14,16
   16 IF(K)17,18,17
C
C	 CHOOSE THE MAXIMUM DIFFERENCE, D
C
   17 D=AMAX1(D,ABS(FLOAT(I)*XN1-FLOAT(J)*XM1))
      IF(L)22,11,22
   18 J=J+1
      IF(J-M)19,20,20
   19 IF(Y(J+1)-Y(J))18,18,17
   20 L=1
      GO TO 17
   21 L=1
      GO TO 16
C
C	 CALCULATE THE STATISTIC Z
C
   22 Z=D*SQRT((XN*XM)/(XN+XM))
	D2=D
C
C	 CALCULATE THE PROBABILITY ASSOCIATED WITH Z
C
      CALL SMIRN(Z,PROB)
      PROB=1.0-PROB
      RETURN
      END
C
C     ..................................................................
C
C	 SUBROUTINE SMIRN
C
C	 PURPOSE
C	    COMPUTES VALUES OF THE LIMITING DISTRIBUTION FUNCTION FOR
C	    THE KOLMOGOROV-SMIRNOV STATISTIC.
C
C	 USAGE
C	    CALL SMIRN(X,Y)
C
C	 DESCRIPTION OF PARAMETERS
C	    X	 - THE ARGUMENT OF THE SMIRN FUNCTION
C	    Y	 - THE RESULTANT SMIRN FUNCTION VALUE
C
C	 REMARKS
C	    Y IS SET TO ZERO IF X IS NOT GREATER THAN 0.27, AND IS SET
C	    TO ONE IF X IS NOT LESS THAN 3.1.  ACCURACY TESTS WERE MADE
C	    REFERRING TO THE TABLE GIVEN IN THE REFERENCE BELOW.
C	    TWO ARGUMENTS, X= 0.62, AND X = 1.87 GAVE RESULTS WHICH
C	    DIFFER FROM THE SMIRNOV TABLES BY 2.9 AND 1.9 IN THE 5TH
C	    DECIMAL PLACE.  ALL OTHER RESULTS SHOWED SMALLER ERRORS,
C	    AND ERROR SPECIFICATIONS ARE GIVEN IN THE ACCURACY TABLES
C	    IN THIS MANUAL.  IN DOUBLE PRECISION MODE, THESE SAME
C	    ARGUMENTS RESULTED IN DIFFERENCES FROM TABLED VALUES BY 3
C	    AND 2 IN THE 5TH  DECIMAL PLACE.  IT IS NOTED IN
C	    LINDGREN (REFERENCE BELOW) THAT FOR HIGH SIGNIFICANCE LEVELS
C	    (SAY, .01 AND .05) ASYMPTOTIC FORMULAS GIVE VALUES WHICH ARE
C	    TOO HIGH ( BY 1.5 PER CENT WHEN N = 80).  THAT IS, AT HIGH
C	    SIGNIFICANCE LEVELS, THE HYPOTHESIS OF NO DIFFERENCE WILL BE
C	    REJECTED TOO SELDOM USING ASYMPTOTIC FORMULAS.
C
C	 SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C	    NONE
C
C	 METHOD
C	    THE METHOD IS DESCRIBED BY W. FELLER-ON THE KOLMOGOROV-
C	    SMIRNOV LIMIT THEOREMS FOR EMPIRICAL DISTRIBUTIONS- ANNALS
C	    OF MATH. STAT., 19, 1948, 177-189, BY N. SMIRNOV--TABLE
C	    FOR ESTIMATING THE GOODNESS OF FIT OF EMPIRICAL
C	    DISTRIBUTIONS- ANNALS OF MATH. STAT., 19, 1948, 279-281,
C	    AND GIVEN IN LINDGREN, STATISTICAL THEORY, THE MACMILLAN
C	    COMPANY, N. Y., 1962.
C
C     ..................................................................
C
C---------------X IS INPUT.  Y IS OUTPUT.
      SUBROUTINE SMIRN(X,Y)
C     DOUBLE PRECISION X,Q1,Q2,Q4,Q8,Y
C
C	 IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE C
C	 IN COLUMN ONE OF THE DOUBLE PRECISION CARD ABOVE SHOULD BE
C	 REMOVED, AND THE C IN COLUMN ONE OF THE STATEMENTS NUMBERED
C	 C   3, C   5, AND C   8 SHOULD BE REMOVED, AND THESE CARDS
C	 SHOULD REPLACE THE STATEMENTS NUMBERED 3, 5, AND 8,
C	 RESPECTIVELY.	ALL ROUTINES CALLED BY THIS ROUTINE MUST ALSO
C	 PROVIDE DOUBLE PRECISION ARGUMENTS TO THIS ROUTINE.
C
C     ..................................................................
C
      IF(X-.27)1,1,2
    1 Y=0.0
      GO TO 9
    2 IF(X-1.0)3,6,6
    3 Q1=EXP(-1.233701/X**2)
C   3 Q1=DEXP(-1.233700550136170/X**2)
      Q2=Q1*Q1
      Q4=Q2*Q2
      Q8=Q4*Q4
      IF(Q8-1.0E-25)4,5,5
    4 Q8=0.0
    5 Y=(2.506628/X)*Q1*(1.0+Q8*(1.0+Q8*Q8))
C   5 Y=(2.506628274631001/X)*Q1*(1.0D0+Q8*(1.0D0+Q8*Q8))
      GO TO 9
    6 IF(X-3.1)8,7,7
    7 Y=1.0
      GO TO 9
    8 Q1=EXP(-2.0*X*X)
C   8 Q1=DEXP(-2.0D0*X*X)
      Q2=Q1*Q1
      Q4=Q2*Q2
      Q8=Q4*Q4
      Y=1.0-2.0*(Q1-Q4+Q8*(Q1-Q8))
    9 RETURN
      END
C
C.......................................................................
C
C	 SUBROUTINE NDTR
C
C	 PURPOSE
C	    COMPUTES Y = P(X) = PROBABILITY THAT THE RANDOM VARIABLE  U,
C	    DISTRIBUTED NORMALLY(0,1), IS LESS THAN OR EQUAL TO X.
C	    F(X), THE ORDINATE OF THE NORMAL DENSITY AT X, IS ALSO
C	    COMPUTED.
C
C	 USAGE
C	    CALL NDTR(X,P,D)
C
C	 DESCRIPTION OF PARAMETERS
C	    X--INPUT SCALAR FOR WHICH P(X) IS COMPUTED.
C	    P--OUTPUT PROBABILITY.
C	    D--OUTPUT DENSITY.
C
C	 REMARKS
C	    MAXIMUM ERROR IS 0.0000007.
C
C	 SUBROUTINES AND SUBPROGRAMS REQUIRED
C	    NONE
C
C	 METHOD
C	    BASED ON APPROXIMATIONS IN C. HASTINGS, APPROXIMATIONS FOR
C	    DIGITAL COMPUTERS, PRINCETON UNIV. PRESS, PRINCETON, N.J.,
C	    1955.  SEE EQUATION 26.2.17, HANDBOOK OF MATHEMATICAL
C	    FUNCTIONS, ABRAMOWITZ AND STEGUN, DOVER PUBLICATIONS, INC.,
C	    NEW YORK.
C
C.......................................................................
C
C---------------X IS INPUT.  P, D ARE RETURNED.
      SUBROUTINE NDTR(X,P,D)
C
      AX=ABS(X)
      T=1.0/(1.0+.2316419*AX)
      D=0.3989423*EXP(-X*X/2.0)
      P = 1.0 - D*T*((((1.330274*T - 1.821256)*T + 1.781478)*T -
     1	0.3565638)*T + 0.3193815)
      IF(X)1,2,2
    1 P=1.0-P
    2 RETURN
      END