Trailing-Edge
-
PDP-10 Archives
-
decuslib10-09
-
43,50466/kolm.f4
There are no other files named kolm.f4 in the archive.
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