Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/stp/stp15.for
There is 1 other file named stp15.for in the archive. Click here to see a list.
C **** STAT PACK ****
C ROUTINE FOR 1 AND 2 SAMPLE KOLMOGOROV-SMIRNOV TESTS.
C CALLING SEQUENCE: CALL KOLMG(NV,NC,MV,MC,DATA,VMN,STD,SAMP1,SAMP2,NAMES)
C WHERE NV NUMBER OF VARIABLES.
C NC - NUMBER OF OBSERVATIONS.
C MV - MAXIMUM NUMBER OF VARIABLES POSSIBLE.
C MC - MAXIMUM NUMBER OF CASES POSSIBLE.
C DATA - MATRIX CONTAINING DATA (MV X MC).
C VMN - VECTOR CONTAINING VARIABLE MEANS.
C STD - VECTOR CONTAINING STANDARD DEVIATIONS.
C SAMP1 - EXTRA VECTOR AT LEAST NC LONG.
C SAMP2 - EXTRA VECTOR AT LEAST NC LONG.
C NAMES - VECTOR CONTAINING VARIABLE NAMES.
C
C KOLM ORIGINALLY REQUESTED BY DAVE DICKERSON GEOGRAPHY DEPARTMENT.
C MAIN STATISTICAL ROUTINES ARE ADAPTED FROM IBM SCIENTIFIC SUBROUTINE
C PACKAGE WITH SOME MODIFICATION AS SUGESTED BY MIKE STOLINE
C COMPUTER CENTER STATISTICAL CONSULTANT.
C
SUBROUTINE KOLMG(NV,NC,MV,MC,DATA,VMN,STD,SAMP1,SAMP2,NAMES)
DIMENSION DATA(MC,MV),VMN(1),STD(1),SAMP1(1),SAMP2(1)
DIMENSION NAMES(1),ISPDF(4),INPUT(80),VALUE(2,4)
DIMENSION IVAR(40),RNGE(20,3),IPLACE(15),EXTRA(3)
EQUIVALENCE(ISPDF(1),ISNORM),(ISPDF(2),ISEXPN)
EQUIVALENCE (ISPDF(3),ISCAUC),(ISPDF(4),ISUNIF)
COMMON /DEV/ ICC,IDATA,IOUT,IDLG,IDSK
COMMON /PRNT/ LINPP,ICOPS,RUNPRG
COMMON /EXTRA/ HEDR(70),NSZ
1 IS2SMP=0
ISSUPL=0
ISBREK=0
ISDISC=0
ISRNGE=0
ISNORM=0
ISEXPN=0
ISCAUC=0
ISUNIF=0
IF(ICC.NE.2) WRITE(IDLG,2)
2 FORMAT(' ENTER OPTIONS SEPARATED BY COMMAS'/)
READ(ICC,3,END=150)(INPUT(I),I=1,16)
3 FORMAT(16(A5,1X))
IF(INPUT(1).EQ.'!') RETURN
IF(INPUT(1).EQ.' ') GO TO 35
I=1
4 IF(INPUT(I).EQ.' ') GO TO 25
IF(INPUT(I).NE.'NORML') GO TO 8
IF(ISNORM.EQ.0) GO TO 7
5 WRITE (IDLG,6) INPUT(I)
6 FORMAT(' SWITCH "',A5,'" LISTED TWICE')
GO TO 1
7 ISNORM=1
GO TO 22
8 IF(INPUT(I).NE.'EXPON') GO TO 9
IF(ISEXPN.NE.0) GO TO 5
ISEXPN=1
GO TO 22
9 IF(INPUT(I).NE.'CAUCH') GO TO 10
IF(ISCAUC.NE.0) GO TO 5
ISCAUC=1
GO TO 22
10 IF(INPUT(I).NE.'UNIFM') GO TO 11
IF(ISUNIF.NE.0) GO TO 5
ISUNIF=1
GO TO 22
11 IF(INPUT(I).NE.'SUPPL') GO TO 12
IF(ISSUPL.NE.0) GO TO 5
ISSUPL=1
GO TO 22
12 IF(INPUT(I).NE.'2SAMP') GO TO 13
IF(IS2SMP.NE.0) GO TO 5
IS2SMP=1
GO TO 22
13 IF(INPUT(I).NE.'BREAK') GO TO 14
IF(ISBREK.NE.0)GO TO 5
ISBREK=1
GO TO 22
14 IF(INPUT(I).NE.'DISCR') GO TO 15
IF(ISDISC.NE.0) GO TO 5
ISDISC=1
GO TO 22
15 IF(INPUT(I).NE.'RANGE') GO TO 16
IF(ISRNGE.NE.0) GO TO 5
ISRNGE=1
GO TO 22
16 IF(INPUT(I).NE.'HELP') GO TO 19
WRITE(IDLG,17)
17 FORMAT(
1'0THE KOLMOGOROV-SMIRNOV TEST ASSUMES ONE SAMPLE TESTS UNLESS'/
2' THE TWO SAMPLE OPTION IS SPECIFIED. OPTIONS INCLUDE:'/
3' NORML - TEST AGAINST A NORMAL DISTRIBUTION (WILL BE ASSUMED'/
4' IF NO OPTION ARE USED)'/
5' EXPON - TEST AGAINST AN EXPONENTIAL DISTRIBUTION'/
6' CAUCH - TEST AGAINST A CAUCHY DISTRIBUTION'/
7' UNIFM - TEST AGAINST A UNIFORM DISTRIBUTION'/
7' TOTAL - TEST AGAINST ALL ABOVE DISTRIBUTIONS'/
8' SUPPL - USER SUPPLIES DETAILED INFORMATION ABOUT'/
9' DISTRIBUTION TO BE TESTED AGAINST'/
1'0------TWO SAMPLE OPTIONS '/
2' 2SAMP - 2 SAMPLE KOLMOGOROV-SMIRNOV TEST'/
3' IF A TWO SAMPLE TEST IS SPECIFIED, IT IS ASSUMED THAT THE'/
4' SAMPLES ARE TWO DIFFERENT VARIABLES. IF THE SAMPLES ARE TO BE'/
5' SELECTED FROM A SINGLE VARIABLE THE BREAK OPTION SHOULD BE'/
6' USED'/
7' BREAK - SELECT THE 2 SAMPLES FROM A SINGLE VARIABLE BASED ON'/
8' THE VALUE OF A SECOND VARIABLE (BREAKDOWN VARIABLE)')
WRITE(IDLG,18)
18 FORMAT(
1' DISCR - RATHER THAN USING USER SPECIFIED RANGES FOR THE'/
2' BREAKDOWN VARIABLE, USE THE INDIVIDUAL VALUES OF THE'/
3' BREAKDOWN VARIABLE. IF THERE ARE MORE THAN 20'/
4' INDIVIDUAL VALUES, THE USER WILL BE REQUIRED'/
3' TO ENTER RANGES.(ONLY AVAILABLE IF BREAK IS USED)'/
4' RANGE - LIST THE RANGES USED. (ONLY AVAILABLE IF DISCR IS'/
5' USED)')
GO TO 1
19 IF(INPUT(I).NE.'TOTAL') GO TO 20
ISNORM=1
ISEXPN=1
ISCAUC=1
ISUNIF=1
GO TO 22
20 WRITE(IDLG,21) INPUT(I)
21 FORMAT(' OPTION "',A5,'" DOES NOT EXIST')
GO TO 1
22 I=I+1
GO TO 4
I=1
25 IF(IS2SMP.NE.1) GO TO 33
IF((ISNORM.NE.1).AND.(ISEXPN.NE.1).AND.(ISCAUC.NE.1).AND.
1(ISUNIF.NE.1)) GO TO 27
WRITE(IDLG,26)
26 FORMAT(
1' THE TYPE OF PROBABILITY DISTRIBUTION FUNCTION CANNOT BE'/
2' SPECIFIED IN A 2 SAMPLE KOLMOGOROV-SMIRNOV TEST')
GO TO 1
27 IF(ISBREK.EQ.1)GO TO 31
IF(ISDISC.NE.1) GO TO 29
WRITE(IDLG,28)
28 FORMAT(' DISCR MAY NOT BE USED WITHOUT BREAK')
GO TO 1
29 IF(ISRNGE.NE.1) GO TO 200
32 WRITE(IDLG,30)
30 FORMAT(' RANGE MAY NOT BE USED WITHOUT DISCR')
GO TO 1
31 IF((ISDISC.EQ.0).AND.(ISRNGE.EQ.1)) GO TO 32
GO TO 200
33 IF((ISBREK.NE.1).AND.(ISDISC.NE.1).AND.(ISRNGE.NE.1)) GO TO 40
WRITE(IDLG,34)
34 FORMAT(' BREAK, DISCR, AND RANGE ARE ONLY AVAILABLE ON 2 SAMPLE ',
1' TESTS')
GO TO 1
35 ISNORM=1
C
C ONE SAMPLE TESTS (FIRST PICK UP SPECIAL DIST FUNCT) SUPPL
C
40 IF(ISSUPL.EQ.0) GO TO 70
IF(ISNORM.NE.1) GO TO 50
43 IF(ICC.NE.2) WRITE(IDLG,41)
41 FORMAT(' ENTER MEAN, STDEV FOR NORMAL DISTRIBUTION? ',$)
CALL NUMBR(INPUT,VALUE(1,1),IERR,IHELP,IRET)
IF(IRET.EQ.1) RETURN
IF(IERR.EQ.1) GO TO 43
IF(IHELP.NE.1) GO TO 50
WRITE(IDLG,44)
44 FORMAT(' ENTER THE MEAN AND STANDARD DEVIATION OF THE NORMAL'/
1' DISTRIBUTION TO BE TESTED AGAINST')
GO TO 43
50 IF(ISEXPN.NE.1) GO TO 60
51 IF(ICC.NE.2) WRITE(IDLG,52)
52 FORMAT(' ENTER MEAN, STDEV FOR EXPONENTIAL DISTRIBUTION? ',$)
CALL NUMBR(INPUT,VALUE(1,2),IERR,IHELP,IRET)
IF(IRET.EQ.1) RETURN
IF(IERR.EQ.1) GO TO 51
IF(IHELP.NE.1) GO TO 60
WRITE(IDLG,53)
53 FORMAT(' ENTER THE MEAN AND STANDARD DEVIATION FOR THE'/
1' EXPONENTIAL DISTRIBUTION TO BE TESTED AGAINST, SEPARATED'/
2' BY A COMMA')
GO TO 51
60 IF(ISCAUC.NE.1) GO TO 65
61 IF(ICC.NE.2) WRITE(IDLG,62)
62 FORMAT(' ENTER FIRST QUARTILE AND MEDIAN FOR CAUCHY DISTRIBUTION?',
1' ',$)
CALL NUMBR(INPUT,VALUE(1,3),IERR,IHELP,IRET)
IF(IRET.EQ.1) RETURN
IF(IHELP.NE.1) GO TO 65
WRITE(IDLG,63)
63 FORMAT(' ENTER FIRST QUARTILE AND MEDIAN SEPARATED BY A COMMA'/
1' FOR THE CAUCHY DISTRIBUTION TO BE TESTED AGAINST')
GO TO 61
65 IF(ISUNIF.NE.1) GO TO 70
66 IF(ICC.NE.2) WRITE(IDLG,67)
67 FORMAT(' ENTER END POINTS FORM UNIFORM DISTRIBUTION? ',$)
CALL NUMBR(INPUT,VALUE(1,4),IERR,IHELP,IRET)
IF(IRET.EQ.1) RETURN
IF(IERR.EQ.1) GO TO 61
IF(IHELP.NE.1) GO TO 70
WRITE(IDLG,68)
68 FORMAT(' ENTER END POINTS SEPARATED BY COMMAS FOR THE UNIFROM'/
1' DISTRIBUTION TO BE TESTED AGAINST')
GO TO 66
70 IF(ICC.NE.2) WRITE(IDLG,71)
71 FORMAT(' ENTER VARIABLES TO BE TESTED SEPARATED BY COMMAS'/)
CALL ALPHA(IVAR,40,NN,IRET,IHELP,IERR,NAMES,NV)
IF(IRET.EQ.1) RETURN
IF(IERR.EQ.1) GO TO 70
IF(IHELP.NE.1) GO TO 80
WRITE(IDLG,72)
72 FORMAT(' ENTER VARIABLES TO BE TESTED AGAINST THE VARIOUS'/
1' DISTRIBUTIONS. VARIABLE NAMES MAY BE USED OR VARIABLE'/
2' NUMBERS. RANGES MAY BE SPECIFIED BY TYPING THE EXTREMES'/
3' OF THE RANGE SEPARATED BY A -')
GO TO 70
80 ALL=0
DO 81 I=1,NN
IF(IVAR(I).GT.0) GO TO 81
ALL=1
NN=NV
GO TO 82
81 CONTINUE
82 IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(J),J=1,NSZ)
5566 FORMAT('1',70A1)
IF(IOUT.EQ.21) CALL PRNTHD
WRITE(IOUT,83)
83 FORMAT('0',10X,'***** KOLMOGOROV-SMIRNOV TEST *****')
WRITE(IOUT,422)
422 FORMAT('0VAR.',3X,'DISTRIBUTION TESTED AGAINST',21X,'D1',10X,
1'PROB')
LINES=6
EXCPN=0
DO 90 I=1,NN
IVB=I
IF(ALL.EQ.1) GO TO 91
IVB=IVAR(I)
91 DO 92 J=1,NC
92 SAMP1(J)=DATA(J,IVB)
IF((ISCAUC.EQ.1).OR.(ISUNIF.EQ.1)) CALL SRTDAT(SAMP1,NC)
ISPACE=NAMES(IVB)
IF(ISNORM.NE.1) GO TO 100
VAL1=VMN(IVB)
VAL2=STD(IVB)
IF(ISSUPL.EQ.0) GO TO 94
VAL1=VALUE(1,1)
VAL2=VALUE(2,1)
94 EXCP=' '
IF(NC.GE.100) GO TO 96
EXCP='*'
EXCPN=1
96 CALL KOLMO(SAMP1,NC,Z,P,1,VAL1,VAL2,IERR,D1)
IF(IOUT.NE.21) GO TO 400
LINES=LINES+1
IF(LINES.LE.LINPP) GO TO 400
CALL PRNTHD
WRITE(IOUT,422)
LINES=7
400 IF(IERR.EQ.0) WRITE(IOUT,93) ISPACE,VAL1,VAL2,D1,P,EXCP
93 FORMAT(1X,A5,2X,'NORMAL MEAN=',G13.6,' STDEV=',G13.6,2X
1,G10.4,1X,F6.3,A1)
IF(IERR.EQ.1) WRITE(IOUT,95) ISPACE
95 FORMAT(1X,A5,2X,'NORMAL HAD INCORRECT PARAMETER')
ISPACE=' '
100 IF(ISEXPN.NE.1) GO TO 105
VAL1=VMN(IVB)
VAL2=STD(IVB)
IF(ISSUPL.EQ.0) GO TO 101
VAL1=VALUE(1,2)
VAL2=VALUE(2,2)
101 CALL KOLMO(SAMP1,NC,Z,P,2,VAL1,VAL2,IERR,D1)
IF(IOUT.NE.21) GO TO 401
LINES=LINES+1
IF(LINES.LE.LINPP) GO TO 401
CALL PRNTHD
WRITE(IOUT,422)
ISPACE=NAMES(IVB)
LINES=7
401 IF(IERR.EQ.0) WRITE(IOUT,102) ISPACE,VAL1,VAL2,D1,P,EXCP
102 FORMAT(1X,A5,2X,'EXPONENTIAL MEAN=',G10.4,' STDEV=',G10.4,
13X,G10.4,1X,F6.3,A1)
IF(IERR.EQ.1) WRITE(IOUT,103)
103 FORMAT(1X,A5,2X,'EXPONENTIAL HAD INCORRECT PARAMETER')
ISPACE=' '
105 IF(ISCAUC.NE.1) GO TO 115
X=(.25*NC)+.5
K=X
IF((K.LE.0).OR.(K.GE.NC)) GO TO 108
Y=K
YY=X-Y
VAL1=((SAMP1(K+1)-SAMP1(K))*YY)+SAMP1(K)
X=.5*(NC+1.)
K=X
IF((K.LE.0).OR.(K.GE.NC)) GO TO 108
Y=K
YY=X-Y
VAL2=((SAMP1(K+1)-SAMP1(K))*YY)+SAMP1(K)
IF(ISSUPL.EQ.0) GO TO 106
VAL1=VALUE(1,3)
VAL2=VALUE(2,3)
106 CALL KOLMO(SAMP1,NC,Z,P,3,VAL1,VAL2,IERR,D1)
IF(IOUT.NE.21) GO TO 402
LINES=LINES+1
IF(LINES.LE.LINPP) GO TO 402
CALL PRNTHD
WRITE(IOUT,422)
ISPACE=NAMES(IVB)
LINES=7
402 IF(IERR.EQ.0) WRITE(IOUT,107) ISPACE,VAL1,VAL2,D1,P,EXCP
107 FORMAT(1X,A5,2X,'CAUCHY 1ST QUART.=',G10.4,' MEDIAN=',G10.4,
11X,G10.4,1X,F6.3,A1)
IF(IERR.EQ.0) GO TO 110
108 WRITE(IDLG,109)
109 FORMAT(1X,A5,2X,'CAUCHY HAD INCORECT PARAMETER')
110 ISAMP=' '
115 IF(ISUNIF.NE.1) GO TO 90
VAL1=SAMP1(1)
VAL2=SAMP1(NC)
IF(ISSUPL.EQ.0) GO TO 116
VAL1=VALUE(1,4)
VAL2=VALUE(2,4)
116 CALL KOLMO(SAMP1,NC,Z,P,4,VAL1,VAL2,P,D1)
IF(IOUT.NE.21) GO TO 403
LINES=LINES+1
IF(LINES.LE.LINPP) GO TO 403
CALL PRNTHD
WRITE(IOUT,422)
ISPACE=NAMES(IVB)
LINES=7
403 IF(IERR.EQ.0) WRITE(IOUT,117) ISPACE,VAL1,VAL2,D1,P,EXCP
117 FORMAT(1X,A5,2X,'UNIFORM FROM',G14.7,' TO ',G14.7,3X,
1G10.4,1X,F6.3,A1)
IF(IERR.EQ.1) WRITE(IOUT,118) ISPACE
118 FORMAT(1X,A5,2X,'UNIFORM HAD INCORRECT PARAMETER')
90 CONTINUE
IF(IOUT.NE.21) GO TO 404
LINES=LINES+4
IF(LINES.LE.LINPP) GO TO 404
CALL PRNTHD
LINES=6
404 WRITE(IOUT,120)
120 FORMAT(
1'0',1X,'THE HYPOTHESIS THAT THE SAMPLE IS FROM THE SPECIFIED'/
22X,'DISTRIBUTION CAN BE REJECTED WITH THE INDICATED PROBABILITY'/
32X,'(PROB) OF BEING INCORRECT.')
IF(EXCPN.NE.1) RETURN
IF(IOUT.NE.21) GO TO 405
LINES=LINES+4
IF(LINES.LE.LINPP) GO TO 405
CALL PRNTHD
LINES=6
405 WRITE(IOUT,121)
121 FORMAT('0*THE PROBABILITY MAY NOT BE CORRECT FOR THE SAMPLE'/
1' SIZES BEING USED. FOR MORE ACCURATE PROBABILITIES'/
2' CHECK THE TABLES IN NON-PARAMETRIC STATISTICS BY SIEGEL')
150 RETURN
200 IF(ISBREK.EQ.1) GO TO 250
201 IF(ICC.NE.2) WRITE(IDLG,202)
202 FORMAT(' ENTER VARIABLES SEPARATED BY COMMAS'/)
CALL ALPHA(IVAR,40,NN,IRET,IHELP,IERR,NAMES,NV)
IF(IERR.EQ.1) GO TO 1
IF(IRET.EQ.1) RETURN
IF(IHELP.NE.1) GO TO 204
WRITE(IDLG,203)
203 FORMAT(' ENTER VARIABLES TO BE TESTED SEPARATED BY COMMAS.'/
1' UP TO 40 VARIABLES MAY BE ENTERED. ALL POSSIBLE PAIRS'/
2' WILL BE SELECTED FROM THE LIST AND A KOLMOGOROV-SMIRNOV'/
3' TEST WILL BE RUN FOR EACH PAIR')
GO TO 201
204 ALL=0
DO 205 I=1,NN
IF(IVAR(I).GT.0) GO TO 205
ALL=1
NN=NV
GO TO 206
205 CONTINUE
206 IF(NN.GE.2) GO TO 208
WRITE(IDLG,207)
207 FORMAT(' IN ORDER TO RUN A TWO SAMPLE TEST, AT LEAST 2 VARIABLES'/
1' MUST BE ENTERED.')
GO TO 201
208 IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(J),J=1,NSZ)
IF(IOUT.EQ.21) CALL PRNTHD
WRITE(IOUT,212)
212 FORMAT('0',10X,'***** KOLMOGOROV-SMIRNOV 2 SAMPLE TEST *****')
WRITE(IOUT,421)
421 FORMAT('0VAR.',2X,'SIZE',5X,'VAR.',2X,'SIZE',8X,'D2',15X,'PROB')
LINES=6
EXCPN=0
DO 209 I=1,NN-1
IVAR1=I
IF(ALL.NE.1) IVAR2=IVAR(I)
DO 210 J=I+1,NN
IVAR2=J
IF(ALL.NE.1)IVAR2=IVAR(J)
DO 211 K=1,NC
SAMP1(K)=DATA(K,IVAR1)
SAMP2(K)=DATA(K,IVAR2)
211 CONTINUE
EXCP=' '
IF(NC.LT.100) EXCP='*'
IF(NC.LT.100) EXCPN=1
CALL KOLM2(SAMP1,SAMP2,NC,NC,Z,P,D2)
IF(IOUT.NE.21) GO TO 406
LINES=LINES+1
IF(LINES.LE.LINPP) GO TO 406
CALL PRNTHD
WRITE(IOUT,421)
LINES=7
406 WRITE(IOUT,213) NAMES(IVAR1),NC,NAMES(IVAR2),NC,D2,P,EXCP
213 FORMAT(1X,A5,1X,I4,5X,A5,1X,I4,5X,G15.7,4X,F6.3,A1)
210 CONTINUE
209 CONTINUE
IF(IOUT.NE.21) GO TO 407
LINES=LINES+4
IF(LINES.LE.LINPP) GO TO 407
CALL PRNTHD
LINES=6
407 WRITE(IOUT,214)
214 FORMAT(
1'0THE HYPOTHESIS THAT THE 2 SAMPLES ARE DRAWN FROM THE SAME'/
2' POPULATION CAN BE REJECTED WITH THE INDICATED PROBABILITY'/
3' (PROB) OF BEING INCORRECT.')
IF(EXCPN.NE.1) RETURN
IF(IOUT.NE.21) GO TO 408
LINES=LINES+4
IF(LINES.LE.LINPP) GO TO 408
CALL PRNTHD
LINES=6
408 WRITE(IOUT,121)
RETURN
250 IF(ICC.NE.2) WRITE(IDLG,251)
251 FORMAT(' ENTER VARIABLES TO BE ANALYZED, SEPARATED BY COMMAS'/)
CALL ALPHA(IVAR,40,NN,IRET,IHELP,IERR,NAMES,NV)
IF(IERR.EQ.1) GO TO 250
IF(IRET.EQ.1) RETURN
IF(IHELP.NE.1) GO TO 266
WRITE(IDLG,252)
252 FORMAT(' ENTER VARIABLES WHICH ARE TO BE ANALYZED'/
1' VARIABLES MAY BE SPECIFIED BY VARIABLE NAMES OR VARIABLE'/
2' NUMBERS. LIST VARIABLES SEPARATED BY COMMAS OR IF RANGES ARE'/
3' TO BE USED ENTER THE EXTREMES OF THE RANGE SEPARATED BY A -')
GO TO 250
266 ALL=0
DO 267 I=1,NN
IF(IVAR(I).GE.1) GO TO 267
ALL=1
NN=NV
GO TO 253
267 CONTINUE
253 IF(ICC.NE.2) WRITE(IDLG,254)
254 FORMAT(' WHICH IS THE BREAKDOWN VARIABLE? ',$)
CALL ALPHA(IBRK,1,I,IRET,IHELP,IERR,NAMES,NV)
IF(IERR.EQ.1) GO TO 253
IF(IRET.EQ.1) RETURN
IF(IHELP.NE.1) GO TO 256
WRITE(IDLG,255)
255 FORMAT(' ENTER THE VARIABLE NUMBER OR NAME OF THE VARIABLE'/
1' TO BE USED AS THE BREAKDOWN VARIABLE. IF DISCR HAS BEEN USED'/
2' RANGES WILL BE AUTOMATICALLY CREATED IF LESS THAN 21'/
3' DISTINCT VALUES EXIST FOR THE BREAKDOWN VARIABLE OTHERWISE'/
4' THE USER WILL BE EXPECTED TO ENTER RANGES.')
GO TO 253
256 IF(ISDISC.EQ.1) GO TO 310
IF(ICC.NE.2) WRITE(IDLG,257) NAMES(IBRK)
257 FORMAT(' ENTER RANGES FOR BREAKDOWN VARIABLE ',A5/)
I=1
258 IF(ICC.NE.2) WRITE(IDLG,259)
259 FORMAT('+? ',$)
READ(ICC,260,END=301) INPUT
260 FORMAT(80A1)
IF(INPUT(1).EQ.'!') RETURN
IF((INPUT(1).NE.'H').OR.(INPUT(2).NE.'E').OR.(INPUT(3).NE.'L').
1OR.(INPUT(4).NE.'P')) GO TO 262
WRITE(IDLG,261)
261 FORMAT(' ENTER RANGES TO BE USED WITH THE BREAKDOWN VARIABLE.'/
1' RANGES ARE LISTED ON SEPARATE LINES, EACH RANGE BEING'/
2' COMPOSED OF A MINIMUM FOR THE RANGE, A COMMA, AND THE'/
3' MAXIMUM OF THE RANGE. WHEN FINISHED ENTERING RANGES TYPE A'/
4' CARRIAGE RETURN OR A ^Z(CONTROL Z).'/)
GO TO 258
262 IF((INPUT(1).EQ.'A').AND.(INPUT(2).EQ.'U').AND.(INPUT(3).EQ.'T')
1.AND.(INPUT(4).EQ.'O')) GO TO 310
IF(INPUT(1).EQ.' ') GO TO 301
263 K=1
L=1
265 DO 270 J=1,15
270 IPLACE(J)=' '
DCLPT=0
EXPN=0
J=1
271 IF(INPUT(L).EQ.',') GO TO 287
IF(INPUT(L).EQ.' ') GO TO 285
IF(INPUT(L).EQ.'E') GO TO 278
IF(INPUT(L).EQ.'-') GO TO 276
IF(INPUT(L).EQ.'.') GO TO 273
IF((INPUT(L).LE.'9').AND.(INPUT(L).GE.'0')) GO TO 281
WRITE(IDLG,272)
272 FORMAT(' NON NUMERIC CHARACTER IN RANGE'/)
GO TO 258
273 IF(DCLPT.EQ.0) GO TO 275
WRITE(IDLG,274)
274 FORMAT(' ONLY 1 DECIMLE POINT PER NUMBER IN RANGE'/)
GO TO 258
275 DCLPT=1
GO TO 281
276 IF(J.EQ.1) GO TO 281
WRITE(IDLG,277)
277 FORMAT(' MINUS MUST BE FIRST CHARACTER OF NEGATIVE NUMBER'/)
GO TO 258
278 IF(EXPN.EQ.0) GO TO 280
WRITE(IDLG,279)
279 FORMAT(' ONLY 1 E(EXPONENT) PER NUMBER'/)
GO TO 258
280 EXPN=1
281 IF(J.GT.15) GO TO 282
IPLACE(J)=INPUT(L)
J=J+1
282 L=L+1
IF(L.LE.80) GO TO 271
285 IF(K.GT.1) GO TO 287
WRITE(IDLG,286)
286 FORMAT(' RANGES MUST BE SEPARATED BY A COMMA'/)
GO TO 258
287 IF(J.GT.1) GO TO 289
WRITE(IDLG,288)
288 FORMAT(' NO SPACES WHEN TYPING RANGE'/)
GO TO 258
289 IF(IPLACE(15).NE.' ') GO TO 291
DO 290 J=15,2,-1
290 IPLACE(J)=IPLACE(J-1)
IPLACE(1)=' '
GO TO 289
291 ENCODE(15,260,EXTRA) IPLACE
DECODE(15,303,EXTRA) RNGE(I,K)
303 FORMAT(F15.0)
K=K+1
IF(K.NE.2) GO TO 293
L=L+1
IF(L.LE.80) GO TO 265
WRITE(IDLG,292)
292 FORMAT(' NO SECOND NUMBER'/)
GO TO 258
293 IF(RNGE(I,1).LE.RNGE(I,2)) GO TO 299
SAVE=RNGE(I,1)
RNGE(I,1)=RNGE(I,2)
RNGE(I,2)=SAVE
299 IF(I.LE.1) GO TO 502
DO 500 J=1,I-1
IF(RNGE(J,2).LT.RNGE(I,1)) GO TO 500
IF(RNGE(J,1).GT.RNGE(I,2)) GO TO 500
WRITE(IDLG,501) J
501 FORMAT(' THE RANGE JUST ENTERED IS IN CONFLICT WITH RANGE NUMBER'
1,I2/' SAMPLES MUST BE INDEPENDENT; REENTER LAST RANGE'/)
GO TO 258
500 CONTINUE
502 IF(INPUT(L).EQ.',') GO TO 295
ENCODE(5,294,RNGE(I,3)) I
294 FORMAT(I3,2X)
GO TO 300
295 L=L+1
DO 296 J=1,5
296 IPLACE(J)=' '
297 IF(INPUT(L).EQ.' ') GO TO 298
IF(J.GT.5) GO TO 298
IPLACE(J)=INPUT(L)
J=J+1
L=L+1
IF(L.LE.80) GO TO 297
298 ENCODE(5,260,RNGE(I,3)) (IPLACE(J),J=1,5)
300 I=I+1
IF(I.LE.20) GO TO 258
301 NRNG=I-1
IF(NRNG.GE.2) GO TO 350
WRITE(IDLG,302)
302 FORMAT(' MUST BE AT LEAST 2 RANGES'/)
RETURN
C
C DISCR USED
C
310 DO 311 J=1,NC
311 SAMP1(J)=DATA(J,IBRK)
CALL SRTDAT(SAMP1,NC)
I=1
RNGE(1,1)=SAMP1(1)
RNGE(1,2)=SAMP1(1)
ENCODE(5,294,RNGE(I,3)) I
DO 312 J=2,NC
IF(SAMP1(J).EQ.RNGE(I,1)) GO TO 312
I=I+1
IF(I.GT.20) GO TO 330
RNGE(I,1)=SAMP1(J)
RNGE(I,2)=SAMP1(J)
ENCODE(5,294,RNGE(I,3)) I
312 CONTINUE
GO TO 340
C MORE THAN 20 INDIVIDUAL VALUES
330 WRITE(IDLG,331)
331 FORMAT(' MORE THAN 20 INDIVIDUAL VALUES EXIST, IT WILL BE'/
1' NECESSARY FOR YOU TO ENTER RANGES')
ISDISC=0
GO TO 256
340 NRNG=I
IF(NRNG.GE.2) GO TO 343
WRITE(IDLG,302)
RETURN
343 IF(ISRNGE.EQ.0) GO TO 350
DO 341 I=1,NRNG
WRITE(IDLG,342)(RNGE(I,J),J=1,2)
342 FORMAT(1X,G10.4,',',G10.4)
341 CONTINUE
350 IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(J),J=1,NSZ)
IF(IOUT.EQ.21) CALL PRNTHD
WRITE(IOUT,351)
351 FORMAT('0',10X,'***** KOLMOGOROV-SMIRNOV 2 SAMPLE TEST *****')
WRITE(IOUT,420)
420 FORMAT(' VAR.',4X,'SMP A SIZE',4X,'SMP B SIZE',5X,'D2',13X,'PROB')
LINES=6
EXCPN=0
DO 352 I=1,NN
IV=I
IF(ALL.EQ.0) IV=IVAR(I)
IF((ALL.EQ.1).AND.(IV.EQ.IBRK)) GO TO 352
ISPACE=NAMES(IV)
IF((IV.EQ.IBRK).AND.(ALL.EQ.1)) GO TO 352
DO 353 J=1,NRNG-1
DO 353 K=J+1,NRNG
J1=0
K1=0
DO 354 M=1,NC
IF(DATA(M,IBRK).LT.RNGE(J,1)) GO TO 355
IF(DATA(M,IBRK).GT.RNGE(J,2)) GO TO 355
J1=J1+1
SAMP1(J1)=DATA(M,IV)
355 IF(DATA(M,IBRK).LT.RNGE(K,1)) GO TO 354
IF(DATA(M,IBRK).GT.RNGE(K,2)) GO TO 354
K1=K1+1
SAMP2(K1)=DATA(M,IV)
354 CONTINUE
EXCP=' '
IF((K1.GE.100).AND.(J1.GE.100)) GO TO 356
EXCP='*'
EXCPN=1
356 CALL KOLM2(SAMP1,SAMP2,J1,K1,Z,P,D2)
IF(IOUT.NE.21) GO TO 409
LINES=LINES+1
IF(LINES.LE.LINPP)GO TO 409
CALL PRNTHD
WRITE(IOUT,420)
LINES=7
409 WRITE(IOUT,357) ISPACE,RNGE(J,3),J1,RNGE(K,3),K1,D2,P,EXCP
ISPACE=' '
357 FORMAT(1X,A5,3X,A5,1X,I4,4X,A5,1X,I4,2X,G15.7,1X,F6.3,A1)
353 CONTINUE
352 CONTINUE
IF(IOUT.NE.21) GO TO 410
LINES=LINES+4
IF(LINES.LE.LINPP) GO TO 410
CALL PRNTHD
LINES=6
410 WRITE(IOUT,214)
IF(EXCPN.NE.1) RETURN
IF(IOUT.NE.21) GO TO 411
LINES=LINES+4
IF(LINES.LE.LINPP) GO TO 411
CALL PRNTHD
LINES=6
411 WRITE(IOUT,121)
RETURN
END
C **** STAT PACK ****
C ROUTINE IS PART OF KOLMOGOROV-SMIRNOV TEST. USED TO SORT DATA.
C CALLING SEQUENCE: CALL SRTDAT(SAMP1,NC)
C WHERE SAMP1 - IS A VECTOR CONTAINING DATA TO BE SORTED.
C NC - NUMBER OF OBSERVATIONS TO BE SORTED.
C
C STANDARD SINGLETON SORT AS USED ELSE WHERE IN THE PACKAGE IS
C ALSO USED HERE. TO INCREAS NUMBER OF OBSERVATIONS WHICH MAY
C BE SORTED INCREASE THE SIZE OF THE DIMENSIONS FOR
C IL AND IU.
C
SUBROUTINE SRTDAT(SAMP1,NC)
DIMENSION SAMP1(1),IL(16),IU(16)
M=1
II=1
J=NC
71 IF(II.GE.J) GO TO 78
72 K=II
IJ=(J+II)/2
T=SAMP1(IJ)
IF(SAMP1(II).LE.T) GO TO 73
SAMP1(IJ)=SAMP1(II)
SAMP1(II)=T
T=SAMP1(IJ)
73 LL=J
IF(SAMP1(J).GE.T) GO TO 75
SAMP1(IJ)=SAMP1(J)
SAMP1(J)=T
T=SAMP1(IJ)
IF(SAMP1(II).LE.T) GO TO 75
SAMP1(IJ)=SAMP1(II)
SAMP1(II)=T
T=SAMP1(IJ)
GO TO 75
74 SAMP1(LL)=SAMP1(K)
SAMP1(K)=TT
75 LL=LL-1
IF(SAMP1(LL).GT.T) GO TO 75
TT=SAMP1(LL)
76 K=K+1
IF(SAMP1(K).LT.T) GO TO 76
IF(K.LE.LL) GO TO 74
IF((LL-II).LE.(J-K)) GO TO 77
IL(M)=II
IU(M)=LL
II=K
M=M+1
GO TO 79
77 IL(M)=K
IU(M)=J
J=LL
M=M+1
GO TO 79
78 M=M-1
IF(M.EQ.0) RETURN
II=IL(M)
J=IU(M)
79 IF((J-II).GE.11) GO TO 72
IF(II.EQ.1) GO TO 71
II=II-1
80 II=II+1
IF(II.EQ.J) GO TO 78
T=SAMP1(II+1)
IF(SAMP1(II).LE.T) GO TO 80
K=II
81 SAMP1(K+1)=SAMP1(K)
K=K-1
IF(T.LT.SAMP1(K)) GO TO 81
SAMP1(K+1)=T
GO TO 80
END
C **** STAT PACK ****
C ROUTINE IS PART OF KOLMOGOROV-SMIRNOV TEST. USED TO INPUT PARAMETERS
C FOR USER SPECIFIED ONE SAMPLE TESTS.
C CALLING SEQUENCE: CALL NUMBR(INPUT,VALUE,IERR,IHELP,IRET)
C WHERE INPUT - VECTOR FOR INPUT OF DATA AT LEAST 80 LONG.
C VALUE - VECTOR OF 2 WORDS USED TO RETURN PARAMETERS.
C IERR - ERROR RETURN: 0 IF NO ERROR, 1 IF ERROR.
C IHELP - HELP RETURN: 0 IF NO HELP REQUESTED, 1 IF HELP REQUESTED.
C IRET - RETURN TO MAIN: 0 IF NO RETURN REQUESTED, 1 IF RETURN
C REQUESTED.
C
C ROUTINE ACCEPTS DATA AND RETURNS VALUES TO MAIN LINE FOR USE AS
C PARAMETERS FOR THE 1 SAMPLE TESTS.
C
SUBROUTINE NUMBR(INPUT,VALUE,IERR,IHELP,IRET)
DIMENSION INPUT(1),VALUE(1),IVL(10),IXTRA(2)
COMMON /DEV/ICC,IDATA,IOUT,IDLG,IDSK
COMMON /PRNT/ LINPP,ICOPS,RUNPRG
IERR=0
IRET=0
IHELP=0
READ(ICC,1,END=27) (INPUT(I),I=1,80)
1 FORMAT(80A1)
IF(INPUT(1).EQ.'!') GO TO 27
IF((INPUT(1).EQ.'H').AND.(INPUT(2).EQ.'E').AND.(INPUT(3).EQ.'L').
1AND.(INPUT(4).EQ.'P')) GO TO 28
I=1
N=1
25 DO 2 J=1,10
2 IVL(J)=' '
J=1
IDECL=0
IEXPN=0
3 IF(INPUT(I).EQ.'-') GO TO 5
IF(INPUT(I).EQ.'.') GO TO 7
IF(INPUT(I).EQ.'E') GO TO 10
IF(INPUT(I).EQ.',') GO TO 15
IF(INPUT(I).EQ.' ') GO TO 15
IF((INPUT(I).LE.'9').AND.(INPUT(I).GE.'0')) GO TO 13
WRITE(IDLG,4)
4 FORMAT(' ILLEGAL CHARACTER "',A1,'" IN PARAMETER')
GO TO 26
5 IF(J.EQ.1) GO TO 13
WRITE(IDLG,6)
6 FORMAT(' MINUS MUST BE FIRST CHARACTER OF NUMBER')
GO TO 26
7 IF(IDECL.EQ.0) GO TO 9
WRITE(IDLG,8)
8 FORMAT(' ONLY 1 DECIMLE PER NUMBER')
GO TO 26
9 IDECL=1
GO TO 13
10 IF(IEXPN.EQ.0) GO TO 12
WRITE(IDLG,11)
11 FORMAT(' ONLY ONE EXPONENT PER NUMBER')
GO TO 26
12 IEXPN=1
13 IF(J.GT.10) GO TO 14
IVL(J)=INPUT(I)
J=J+1
14 I=I+1
IF(I.LE.80) GO TO 3
15 IF(J.GT.1) GO TO 17
WRITE(IDLG,16)
16 FORMAT(' NO SPACES BEFOR THE NUMBERS OR BETWEEN THE NUMBERS')
GO TO 26
17 IF(IVL(10).NE.' ') GO TO 20
DO 18 J=9,1,-1
18 IVL(J+1)=IVL(J)
IVL(1)=' '
GO TO 17
20 ENCODE(10,1,IXTRA) IVL
DECODE(10,22,IXTRA) VALUE(N)
22 FORMAT(F10.0)
N=N+1
IF(N.EQ.3) RETURN
IF(INPUT(I).EQ.',') GO TO 24
WRITE(IDLG,23)
23 FORMAT(' THERE MUST BE A COMMA BETWEEN THE TWO NUMBERS')
GO TO 26
24 I=I+1
GO TO 25
26 IERR=1
RETURN
27 IRET=1
RETURN
28 IHELP=1
RETURN
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
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
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
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
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