Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-09 - 43,50466/stp15.stp
There is 1 other file named stp15.stp 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