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 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