//EDA JOB /*JOBPARM S=5,R=2048,L=4 /*SERVICE -4 //STEP1 EXEC FORTVCLG,RC=2048K,RL=2048K,RG=2048K, // PARM.FORT='LANGLVL(77),NOMAP,NOXREF,OPTIMIZE(1)', // PARM.LKED='LET' //FORT.SYSIN DD * IMPLICIT REAL*4(A-H,O-Z) EXTERNAL CDF DIMENSION DX(2000), DY(2000) , DZ(2000) DIMENSION DT(2000), INDEX(2000) DIMENSION DL1(2000), DL2(2000), DL3(2000) ,DL4(2000) DIMENSION GDOP(2000), XDOP(2000), YDOP(2000) ,ZDOP(2000) DIMENSION INDX(2000), B(2000), C(2000) DIMENSION QUANT(2000), PROB(2000) DIMENSION PSIGMA(1000), SPRD(1000), SMIDS(1000) DIMENSION ZSQR(1000), YY(1500) DIMENSION GP(15), RES(15) DIMENSION CELLS(30), COMP(30) DIMENSION SMOOTH(1000), ROUGH(1000) C REAL ARRAY(2000), CLEAN(2000, 4) REAL TABLE(2000, 4) REAL ARRAY1(5000) REAL SORTY(5000), Y(5000), D(15), YLV(15,2) REAL IEPSI, MAXINT REAL NICNOS(2000), FRACT, UNIT REAL ON(2000), WITH(2000), Z(2000) REAL X, LO, HI, HL, HH, NPW, MED, ADJL, ADJH, STEP C INTEGER GSUB(5000), ICLEAN(4) INTEGER IOUTL, IOUTH INTEGER VERSN INTEGER IA, IQ, ND INTEGER POSN, CHAR INTEGER IADJL, IADJH INTEGER IOUNIT, IPMIN, IPMAX, IMAXIN INTEGER P(130), PMAX, PMIN, OUTPTR, MAXPTR, OUNIT INTEGER N, NLV, ERR, NG INTEGER NDATA, I, W INTEGER NN, MAXP, PTOTL C LOGICAL LINE3, NOTCH LOGICAL MZERO COMMON/CHRBUF/P, PMAX, PMIN, OUTPTR, MAXPTR, OUNIT COMMON/NUMBRS/ EPSI, MAXINT C COMMON/CHARIO/ CHARS, CMAX, 1 CHA, CHB, CHC, CHD, CHE, CHF, CHG, CHH, CHI, CHJ, CHK, 2 CHL, CHM, CHN, CHO, CHP, CHQ, CHR, CHS, CHT, CHU, CHV, 3 CHW, CHX, CHY, CHZ, CH0, CH1, CH2, CH3, CH4, CH5, CH6, 4 CH7, CH8, CH9, CHBL, CHEQ, CHPLUS, CHMIN, CHSTAR, CHSLSH, 5 CHLPAR, CHRPAR, CHCOMA, CHPT IPRINT = 6 ITEMP =15 IDISK =22 IREAD = 5 IOPT = 0 C------------------------------------------------------------------- C IOPT = 0 WITH OUTLIERS IN THE ANALYSIS C IOPT = 1 WITHOUT OUTLIERS IN THE ANALYSIS(REMOVED) C-------------------------------------------------------------------- IOUNIT = 6 IPMIN = 1 IPMAX = 130 IEPSI = 1.0E-06 IMAXIN = 2**15 C CALL CINIT(IOUNIT, IPMIN, IPMAX, IEPSI, IMAXIN, ERR) WRITE(6,1001) ERR C-------------------------------------------------- C READ DATA C--------------------------------------------------- OPEN(IDISK,ERR=3000,STATUS='OLD',ACCESS='SEQUENTIAL', # IOSTAT=IERROR) C 9 DO 1 K=1,2000 READ(17,139,ERR=2010,IOSTAT=IERROR) L, GDOP(K), XDOP(K), # YDOP(K), ZDOP(K) READ(IDISK,142,ERR=2010,IOSTAT=IERROR,END=2) INDEX(K), DX(K), # DY(K), DZ(K), DT(K), DL1(K), DL2(K), DL3(K), DL4(K) 1 CONTINUE 2 NDATA = K - 1 NDATA = 600 C--------------------------------------------------- C FORM DATA TABLE C-------------------------------------------------- DO 3 I= 1,NDATA TABLE(I,1) = DX(I) - 1.877 TABLE(I,2) = DY(I) - 2.280 TABLE(I,3) = DZ(I) - 0.635 TABLE(I,4) = SQRT( DX(I)**2 + DY(I)**2 + DZ(I)**2 ) 3 CONTINUE C C----------------------------------------------------------- C DERIVE LETTER-VALUES C----------------------------------------------------------- NG = 1 LINE3 = .TRUE. NOTCH = .FALSE. IFIRST = 1 ILAST = 0 C DO 9999 K= 1, 3 CALL SET(ARRAY, NDATA,1,NDATA, TABLE(1,K), NDATA) C-------------------------------------------------------------------- C CHECK FOR 'OUTSIDE' AND 'FAR OUT' VALUES C------------------------------------------------------------------- NPL = 2 CALL SORTS(ARRAY, NDATA, 1, NDATA, -1, INDX, NPL, # B, C, ARRAY, ARRAY, ARRAY, ARRAY, ARRAY, ARRAY, ARRAY, ARRAY) C IF(IOPT .EQ. 0) THEN CALL SET(CLEAN(1,K), NDATA,1,NDATA, ARRAY, NDATA) ICLEAN(K) = NDATA GO TO 7 ENDIF CALL YINFO(YY, NDATA, MED, HL, HH, ADJL, ADJH, IADJL, IADJH, # STEP,IOUTL,IOUTH,ERR) C WRITE(6,1006) MED C WRITE(6,1007) HL, HH, STEP C WRITE(6,1008) ADJL, IADJL, XDOP(INDX(IADJL)), C # ADJH, IADJH, XDOP(INDX(IADJH)) C WRITE(6,1009) ARRAY(IOUTL), IOUTL, XDOP(INDX(IOUTL)), C # ARRAY(IOUTH), IOUTH, XDOP(INDX(IOUTH)) C------------------------------------------------------------------- C PRINT OUTSIDE AND FAR OUT VALUES C------------------------------------------------------------------- C CALL OUT(YY, NDATA, IOUTL, IOUTH, IADJL, IADJH, GDOP, INDX) C--------------------------------------------------------------------- C REMOVE FAR OUT VALUES FROM ARRAY AND CREATE CLEAN VECTOR C RERUN LETTER-VALUE DISPLAYS AND BOX PLOTS C--------------------------------------------------------------------- MAXOUT = IADJH MINOUT = IADJL CALL CLEAR(YY, NDATA, MAXOUT, MINOUT, CLEAN(1,K)) ICLEAN (K) = (MAXOUT - 1) - (MINOUT + 1) + 1 C 7 CALL YINFO(CLEAN(1,K), ICLEAN(K), MED, HL, HH, ADJL, ADJH, # IADJL, IADJH, STEP,IOUTL,IOUTH,ERR) DO 4 II = 1, ICLEAN(K) CLEAN(II , K) = CLEAN(II , K) - MED 4 CONTINUE C CALL LVALS(CLEAN(1,K),ICLEAN(K), D, YLV, NLV, SORTY, ERR) WRITE(6,1002) ERR CALL LVPRNT(NLV, D, YLV, ERR) WRITE(6,1003) ERR C ILAST = ILAST + ICLEAN(K) CALL SET(Y, 5000, IFIRST, ILAST,CLEAN(1,K), NDATA) CALL SETCON(GSUB,5000, IFIRST, ILAST, K) IFIRST = IFIRST + ICLEAN(K) 9999 CONTINUE NG = 3 WRITE(6,1010) CALL BOXES(Y, ILAST, GSUB, NG, LINE3, NOTCH, SORTY, ERR) C--------------------------------------------------------------- C PERFORM QUANTILE-QUANTILE PLOTS C--------------------------------------------------------------- L = 1 CALL QQTEST( CLEAN(1,L),ICLEAN(L),6, QUANT, PROB,IERR) C DO 5 I=1, ICLEAN(L) C WRITE(ITEMP,*) CLEAN(I,L), QUANT(I) C WRITE(6,*) CLEAN(I,L), QUANT(I) C5 CONTINUE C-------------------------------------------------------------- C FULL NORMAL PROBABILITY PLOT/ADRERSON-DARLING TEST C CHI-SQUARED TEST C-------------------------------------------------------------- C CALL FUNOP(CLEAN(1,L), ICLEAN(L), QUANT, Z) C DO 8 I=1,ICLEAN(L) C WRITE(ITEMP,*) Z(I), I C WRITE(6,*) Z(I), I C8 CONTINUE C CALL ANDRSN(CLEAN(1,L),ICLEAN(L), AVALUE, ACRT) C CALL AGOSTN(CLEAN(1,L),ICLEAN(L), DA, YA) C ALFA = 0.95 C CALL ZTEST(CLEAN(1,L), ICLEAN(L), ALFA, R, ZVALUE, ZCRIT) C------------------------------------------------------------------ C CHI-SQUARED TEST C----------------------------------------------------------------- C IDFR = 2 C KCELL = 30 C CALL GFIT(CDF,KCELL,CLEAN(1,L),ICLEAN(L),CELLS,COMP,CS,IDFR, C # Q, IER) C CHI = CHISQ(0.95, 27) C-------------------------------------------------- C WRITE(6,77) ZVALUE, ZCRIT, R C77 FORMAT('1','Z-VALUE =',F10.4,2X,'CR=',F10.3,2X,'R =',F8.5,//) C------------------------------------------------------------- C PLOTS FOR SKEWNESS AND ELONGATION DIAGNOSIS C------------------------------------------------------------- M = (ICLEAN(L) + 1.0) / 2.0 - 1.0 NOBS = ICLEAN(L) DO 555 I=1, M WRITE(ITEMP,*) ABS(CLEAN(I,L)), CLEAN(NOBS+1-I,L) WRITE(6,*) ABS(CLEAN(I,L)), CLEAN(NOBS+1-I,L) 555 CONTINUE C C CALL LVALS(CLEAN(1,L),ICLEAN(L), D, YLV, NLV, SORTY, ERR) C CALL SKEW(YLV, D, NLV, SPRD, SMIDS, PSIGMA, ZSQR) C CALL PUSH(YLV, NLV, PSIGMA,ZSQR, IPRINT,ITEMP, SMED, RES) CALL GSKEW(YLV, NLV, D, GP) C CALL RSM(GP, NLV, SMOOTH, ROUGH, VERSN, ERR) DO 6 I=1, NLV C WRITE(ITEMP,*) SMOOTH(I), ZSQR(I) C WRITE(6,*) GP(I),SMOOTH(I), ZSQR(I), I 6 CONTINUE CALL TAILW(CLEAN(1,L),ICLEAN(L), WTAIL, WGAUSS) C-------------------------------------------------- C CLOSE FILE ON DISK C------------------------------------------------- GO TO 2020 3000 WRITE(IPRINT,9040) IERROR GO TO 2020 2010 WRITE(IPRINT,9030) IERROR GO TO 2020 2020 CLOSE(IDISK,ERR=4010,STATUS='KEEP',IOSTAT=IERROR) GO TO 2222 C--------------------------------------------------------------------- 139 FORMAT(2X,I8,4(F20.2,2X)) 142 FORMAT(I8,1X,3(F10.3,2X),F20.4,F20.4,1X,3(F14.2)) 1001 FORMAT(2X,'ERROR FLAG FROM SUBROUTINE CINIT=',I5) 1002 FORMAT(2X,'ERROR FLAG FROM SUBROUTINE LVALS=',I5) 1003 FORMAT(2X,'ERROR FLAG FROM SUBROUTINE LVPRTN=',I5) 1004 FORMAT(2X,'ERROR FLAG FROM SUBROUTINE BOXES=',I5) 1005 FORMAT('1',2X,' BOX PLOT # ',//) 1006 FORMAT(34X,'MEDIAN=',F10.3,/) 1007 FORMAT(4X,'LOW HINGE=',F10.3,24X,'HIGH HINGE=',F10.3,4X,'STEP=', # F10.3,/) 1008 FORMAT(4X,'LOW ADJACENT VALUE=',F10.3,2X,'INDEX=',I5,2X, #'ITS GDOP=',F10.3,2X, #'HIGH ADJACENT VALUE=',F10.3,2X,'INDEX=',I5,2X,'ITS GDOP=',F10.3/) 1009 FORMAT(4X,'OUTSIDE LOW VALUE=',F10.3,2X,'INDEX=',I5,2X, #'ITS GDOP=',F10.3,2X, #'OUTSIDE HIGH VALUE=',F10.3,2X,'INDEX=',I5,2X,'ITS GDOP=',F10.3/) 1010 FORMAT('1',30X,'SIMULTANEOUS BOXPLOT DISPLAYS',////) 9030 FORMAT(' ','ERROR',I4,3X,'WHILE READING FILE OF RESULTS') 9040 FORMAT(' ','ERROR',I4,3X,'WHILE OPENING FILE OF RESULTS') 9070 FORMAT(' ','ERROR',I4,3X,'WHILE CLOSING FILE OF RESULTS') C------------------------------------------------------ 4010 WRITE(IPRINT,9070) IERROR 2222 CONTINUE STOP END C BLOCK DATA C------------------------------------------------------------------- C CHARS = SYMBOLS OF THE STANDARD FORTRAN CHARACTER SET C CHA-CHPT ARE THE CORRESPONDING INDICES INTO CHARS C PUTCHR= PRIMARY USER OF THIS TRANSLATION VECTOR C------------------------------------------------------------------ COMMON/CHARIO/ CHARS, CMAX, 1 CHA, CHB, CHC, CHD, CHE, CHF, CHG, CHH, CHI, CHJ, CHK, 2 CHL, CHM, CHN, CHO, CHP, CHQ, CHR, CHS, CHT, CHU, CHV, 3 CHW, CHX, CHY, CHZ, CH0, CH1, CH2, CH3, CH4, CH5, CH6, 4 CH7, CH8, CH9, CHBL, CHEQ, CHPLUS, CHMIN, CHSTAR, CHSLSH, 5 CHLPAR, CHRPAR, CHCOMA, CHPT C INTEGER CHARS(46), CMAX INTEGER CHA, CHB, CHC, CHD, CHE, CHF, CHG, CHH, CHI INTEGER CHJ, CHK, CHL, CHM, CHN, CHO, CHP, CHQ, CHR INTEGER CHS, CHT, CHU, CHV, CHW, CHX, CHY, CHZ INTEGER CH0, CH1, CH2, CH3, CH4, CH5, CH6, CH7, CH8, CH9 INTEGER CHBL, CHEQ, CHPLUS, CHMIN, CHSTAR, CHSLSH INTEGER CHLPAR, CHRPAR, CHCOMA, CHPT C DATA CHARS( 1),CHARS( 2),CHARS( 3),CHARS( 4) /1HA,1HB,1HC,1HD/ DATA CHARS( 5),CHARS( 6),CHARS( 7),CHARS( 8) /1HE,1HF,1HG,1HH/ DATA CHARS( 9),CHARS(10),CHARS(11),CHARS(12) /1HI,1HJ,1HK,1HL/ DATA CHARS(13),CHARS(14),CHARS(15),CHARS(16) /1HM,1HN,1HO,1HP/ DATA CHARS(17),CHARS(18),CHARS(19),CHARS(20) /1HQ,1HR,1HS,1HT/ DATA CHARS(21),CHARS(22),CHARS(23),CHARS(24) /1HU,1HV,1HW,1HX/ DATA CHARS(25),CHARS(26),CHARS(27),CHARS(28) /1HY,1HZ,1H0,1H1/ DATA CHARS(29),CHARS(30),CHARS(31),CHARS(32) /1H2,1H3,1H4,1H5/ DATA CHARS(33),CHARS(34),CHARS(35),CHARS(36) /1H6,1H7,1H8,1H9/ DATA CHARS(37),CHARS(38),CHARS(39),CHARS(40) /1H ,1H=,1H+,1H-/ DATA CHARS(41),CHARS(42),CHARS(43),CHARS(44) /1H*,1H/,1H(,1H)/ DATA CHARS(45),CHARS(46) /1H,,1H./ C DATA CMAX /46/ DATA CHA,CHB,CHC,CHD,CHE,CHF / 1, 2, 3, 4, 5, 6/ DATA CHG,CHH,CHI,CHJ,CHK,CHL / 7, 8, 9,10,11,12/ DATA CHM,CHN,CHO,CHP,CHQ,CHR /13,14,15,16,17,18/ DATA CHS,CHT,CHU,CHV,CHW,CHX /19,20,21,22,23,24/ DATA CHY,CHZ,CH0,CH1,CH2,CH3 /25,26,27,28,29,30/ DATA CH4,CH5,CH6,CH7,CH8,CH9 /31,32,33,34,35,36/ DATA CHBL,CHEQ,CHPLUS,CHMIN /37,38,39,40/ DATA CHSTAR,CHSLSH,CHLPAR,CHRPAR /41,42,43,44/ DATA CHCOMA,CHPT /45,46/ END C SUBROUTINE PUSH(YLV, NLV, PSIGMA,ZSQR, IPRINT,ITEMP, SMED, RES) IMPLICIT REAL * 4(A-H, O-Z) REAL MEDIAN DIMENSION YLV(15,2), PSIGMA(NLV), RES(NLV), ZSQR(NLV) DIMENSION WORK(15) CHARACTER * 1 TAG(15) C------------------------------------------------------------ C PUSHBACK ANALYSIS FROM THE QUANTILE-QUANTILE PLOT C AUTHOR : S.MERTIKAS C DATE : JULY,1986 C------------------------------------------------------------ DATA TAG( 1), TAG( 2), TAG( 3), TAG( 4) /'M','F','E','D'/ DATA TAG( 5), TAG( 6), TAG( 7), TAG( 8) /'C','B','A','Z'/ DATA TAG( 9), TAG(10), TAG(11) /'Y','X','1'/ C WRITE(IPRINT,1000) DO 30 I= 2, NLV J = I - 1 WORK(J) = PSIGMA(I) 30 CONTINUE CALL SORT(WORK, NLV-1, IERROR) SMED = MEDIAN(WORK, NLV-1) J = NLV DO 10 I = 1, NLV Z = - SQRT( ZSQR(J) ) DOT = SMED * Z RES(J) = YLV(J,1) - DOT WRITE(ITEMP, *) RES( J ), Z IF(J .EQ. NLV) TAG(J) = TAG(11) WRITE(IPRINT,2000) TAG(J), YLV(J,1), Z, DOT, RES(J) J = J - 1 10 CONTINUE C DO 20 I=1,NLV Z = SQRT( ZSQR(I) ) DOT = SMED * Z RES(I) = YLV(I,2) - DOT WRITE(ITEMP, *) RES( I ), Z WRITE(IPRINT,2000) TAG(I), YLV(I,2), Z, DOT, RES(I) 20 CONTINUE C 1000 FORMAT('1',3X,'TAG',7X,'LETTER',8X,'Z',10X,' SMED*Z ', #5X,'RESIDUAL',//) 2000 FORMAT(' ',4X,A1,3X,F10.4,3X,F10.4,3X,F10.4,3X,F10.4) RETURN END C SUBROUTINE TAILW(ARRAY, N, WTAIL, WGAUSS) IMPLICIT REAL * 4 (A-H, O-Z) DIMENSION ARRAY( N ) C------------------------------------------------------------------ C COMPUTES THE SAMPLE ESTIMATE OF TAIL-WEIGHT INDEX AND C THE CORRESPONDING THEORETICAL ONE (GAUSSIAN). C C AUTHOR : S.MERTIKAS C DATE : JULY, 1986 C RE : UNDERSTANDING ROBUST & EXPLORATORY DATA ANALYSIS C----------------------------------------------------------------- C LOCAL DIMENSIONS C DIMENSION P( 4 ), Q( 4 ), G( 4 ) DATA P(1), P(2), P(3), P(4) /0.90, 0.75, 0.25, 0.10 / ALFA = (1.0 / 3.0) DO 10 I =1,4 D = (N + ALFA) * P(I) + ALFA ID = D R = D - ID Q(I) = (1.0 - R) * ARRAY(ID) + R * ARRAY(ID + 1) CALL MDNRIS( P(I), G(I), IER) 10 CONTINUE WTAIL = ( Q(1) - Q(4) ) / ( Q(2) - Q(3) ) WGAUSS= ( G(1) - G(4) ) / ( G(2) - G(3) ) RATIO = WTAIL / WGAUSS WRITE(6, 1000) WTAIL, WGAUSS, RATIO 1000 FORMAT('1',4X,'SAMPLE INDEX OF TAIL-WEIGHT =',F10.5,3X, #'GAUSSIAN VALUE=',F10.5,3X,'RATIO =',F10.5) RETURN END C SUBROUTINE GSKEW(YLV, NLV, D, GP) IMPLICIT REAL * 4 (A-H, O-Z) DIMENSION YLV(15,2), D(15), GP(NLV) CHARACTER*1 TAG(15) C------------------------------------------------------------------- C COMPUTES THE SKEWNESS PARAMETER GP FROM A SET OF LETTER-VALUES C C YLV = VECTOR OF LETTER-VALUES AS IN LVALS SUBROUTINE(IN) C NLV = NUMBER OF LETTER-VALUES AS IN LVALS SUBROUTINE(IN) C D = VECTOR OF CORRESPONDING DEPTHS(IN) C GP = VECTOR OF SKEWNESS PARAMETERS(OUT) C C AUTHOR : S. MERTIKAS C DATE : JULY, 1986 C------------------------------------------------------------------- DATA TAG( 1), TAG( 2), TAG( 3), TAG( 4) /'M','F','E','D'/ DATA TAG( 5), TAG( 6), TAG( 7), TAG( 8) /'C','B','A','Z'/ DATA TAG( 9), TAG(10), TAG(11) /'Y','X','1'/ C NDATA = 2.0 * D(1) - 1.0 EN = FLOAT(NDATA) ALFA = (1.0 / 3.0) WRITE(6, 100) C DO 10 I = 2, NLV XUPPER = YLV(I,2) - YLV(1,1) XLOWER = YLV(1,1) - YLV(I,1) RATIO = ALOG( XUPPER / XLOWER ) TAIL = ( D(I) - ALFA ) / (EN + ALFA) CALL MDNRIS( TAIL, Z, IER) Z = ABS( Z ) GP(I) = RATIO / Z T = TAIL * 100. WRITE(6,1000) TAG(I),T, YLV(I,1), YLV(I,2), RATIO, Z, GP(I), I 10 CONTINUE 100 FORMAT('1',3X,'TAG',7X,'TAIL',6X, #'LOWER',6X,'UPPER',6X,'RATIO',6X,'Z',7X,'GP') 1000 FORMAT(3X,A1,7X,F6.3,'%',1X,F10.3,1X,F10.3,1X, #F10.5,2X,F7.4,2X,F6.4) RETURN END C SUBROUTINE SKEW(YLV, D, NLVS, SPRD, SMIDS, PSIGMA, ZSQR) IMPLICIT REAL*4(A-H,O-Z) DIMENSION SPRD(NLVS), SMIDS(NLVS), PSIGMA(NLVS), ZSQR(NLVS) DIMENSION D(15), YLV(15,2) C------------------------------------------------------------------- C COMPUTES MIDSUMMARIES, SPREADS, AND Z**2 FOR SKEWNESS DIAGNOSIS C YLV = VECTOR OF LETTER-VALUES AS IN LVALS SUBROUTINE(IN) C D = VECTOR OF CORRESPONDING DEPTHS(IN) C NLVS = NUMBER OF LETTER-VALUES AS IN LVALS SUBROUTINE(IN) C SPRD = VECTOR OF LETTER-VALUE SPREADS(OUT) C SMIDS = VECTOR OF LETTER-VALUE MIDSUMMARIES(OUT) C PSIGMA = VECTOR OF LETTER-VALUE PSEUDOSIGMA WHEN GAUSSIAN(OUT) C ZSQR = VECTOR OF CORRESPONDING GAUSSIAN QUANTILES SQUARED(OUT) C C AUTHOR : S. MERTIKAS C DATE : JULY, 1986 C------------------------------------------------------------------- NDATA = 2.0 * D(1) - 1.0 ALFA = 1.0 / 3.0 DO 10 I=1, NLVS SPRD(I) = YLV(I, 2) - YLV(I, 1) SMIDS(I) =(YLV(I, 2) + YLV(I, 1)) / 2.0 TAIL = (D(I) - ALFA ) / ( NDATA + ALFA) CALL MDNRIS( TAIL, QNRML, IER) WRITE(6,*) YLV(I, 1) , YLV(I,2), I ZSQR(I) = QNRML ** 2 PSIGMA(I)= SPRD(I) / (2 * ABS(QNRML)) 10 CONTINUE RETURN END C SUBROUTINE CDF(E,P) IMPLICIT REAL*4(A-H,O-Z) U = -0.2342 S = 10.1223 T = (E - U) / S P = .5 * ERFC(-T*.7071068) RETURN END C SUBROUTINE RSM(Y, N, SMOOTH, ROUGH, VERSN, ERR) INTEGER N, VERSN, ERR REAL Y(N), SMOOTH(N), ROUGH(N) INTEGER I C IF (N.GT.6) GO TO 10 ERR = 61 GO TO 999 10 DO 20 I=1,N SMOOTH(I) = Y(I) 20 CONTINUE IF (VERSN .EQ. 1) CALL S3RSSH(SMOOTH, N, ERR) IF (VERSN .EQ. 2) CALL S4253H(SMOOTH, N, ERR) IF (ERR .NE. 0) GO TO 999 C DO 30 I=1,N ROUGH(I) = Y(I) - SMOOTH(I) 30 CONTINUE C IF (VERSN .EQ. 1) CALL S3RSSH(ROUGH, N, ERR) IF (VERSN .EQ. 2) CALL S4253H(ROUGH, N, ERR) IF (ERR .NE. 0) GO TO 999 DO 40 I=1, N SMOOTH(I) = SMOOTH(I) + ROUGH(I) ROUGH(I) = Y(I) - SMOOTH(I) 40 CONTINUE 999 RETURN END C SUBROUTINE S3RSSH(Y, N, ERR) INTEGER N, ERR REAL Y(N) LOGICAL CHANGE CALL S3R(Y, N) CHANGE = .FALSE. CALL SPLIT(Y, N, CHANGE) IF( .NOT. CHANGE) GO TO 10 CALL S3R(Y, N) CHANGE = .FALSE. CALL SPLIT(Y, N, CHANGE) IF(CHANGE) CALL S3R(Y, N) 10 CALL HANN(Y, N) 999 RETURN END C SUBROUTINE S4253H(Y, N, ERR) INTEGER N, ERR REAL Y(N) REAL ENDSAV, WORK(5), SAVE(5) INTEGER NW LOGICAL CHANGE DATA NW /5/ CHANGE = .FALSE. CALL S4(Y, N, ENDSAV, WORK, SAVE, NW, ERR) IF(ERR .EQ. 0) CALL S2(Y, N, ENDSAV) IF(ERR .EQ. 0) CALL S5(Y, N, WORK, SAVE, NW, ERR) IF(ERR .EQ. 0) CALL S3(Y, N, CHANGE) IF(ERR .EQ. 0) CALL ENDPTS(Y, N) IF(ERR .EQ. 0) CALL HANN(Y, N) 999 RETURN END C SUBROUTINE S4(Y, N, ENDSAV, WORK, SAVE, NW, ERR) INTEGER N, NW, ERR REAL Y(N), ENDSAV, WORK(NW), SAVE(NW) REAL ENDM1, TWO DATA TWO /2.0/ C ENDSAV = Y(N) ENDM1 = Y(N-1) CALL RUNMED(Y, N, 4, WORK, SAVE, NW, ERR) Y(2) = (Y(1) + Y(2)) / TWO Y(N) = (ENDM1 + ENDSAV) / TWO 999 RETURN END C SUBROUTINE S2(Y, N, ENDSAV) INTEGER N REAL Y(N), ENDSAV INTEGER NM1, I REAL TWO DATA TWO /2.0/ NM1 = N - 1 DO 10 I=2, NM1 Y(I) = (Y(I+1) + Y(I)) / TWO 10 CONTINUE Y(N) = ENDSAV 999 RETURN END C SUBROUTINE S5(Y, N, WORK, SAVE, NW, ERR) INTEGER N, NW, ERR REAL Y(N), WORK(NW), SAVE(NW) LOGICAL CHANGE REAL YMED1, YMED2 CHANGE = .FALSE. CALL MEDOF3(Y(1), Y(2), Y(3), YMED1, CHANGE) CALL MEDOF3(Y(N), Y(N-1), Y(N-2), YMED2, CHANGE) CALL RUNMED(Y, N, 5, WORK, SAVE, NW, ERR) Y(2) = YMED1 Y(N-1) = YMED2 999 RETURN END C SUBROUTINE HANN(Y, N) INTEGER N REAL Y(N) INTEGER I, NM1 REAL Y1, Y2, Y3 NM1 = N -1 Y2 = Y(1) Y3 = Y(2) DO 10 I=2,NM1 Y1 = Y2 Y2 = Y3 Y3 = Y(I + 1) Y(I) = (Y1 + Y2 + Y2 + Y3) / 4.0 10 CONTINUE 999 RETURN END C SUBROUTINE S3(Y, N, CHANGE) INTEGER N REAL Y(N) LOGICAL CHANGE REAL Y1, Y2, Y3 INTEGER NM1 Y2 = Y(1) Y3 = Y(2) NM1 = N -1 DO 10 I=2,NM1 Y1 = Y2 Y2 = Y3 Y3 = Y(I+1) CALL MEDOF3(Y1, Y2, Y3, Y(I), CHANGE) 10 CONTINUE 999 RETURN END C SUBROUTINE S3R(Y, N) INTEGER N REAL Y(N) LOGICAL CHANGE 10 CHANGE = .FALSE. CALL S3(Y, N, CHANGE) IF(CHANGE) GO TO 10 CALL ENDPTS(Y, N) 999 RETURN END C SUBROUTINE MEDOF3(X1, X2, X3, XMED, CHANGE) REAL X1, X2, X3, XMED LOGICAL CHANGE REAL Y1, Y2, Y3 Y1 = X1 Y2 = X2 Y3 = X3 XMED = Y2 IF( (Y2 - Y1) * (Y3 - Y2) .GE. 0.0 ) GO TO 999 CHANGE = .TRUE. XMED = Y1 IF( (Y3 - Y1) * (Y3 - Y2) .GT. 0.0 ) GO TO 999 XMED = Y3 999 RETURN END C SUBROUTINE ENDPTS(Y, N) INTEGER N REAL Y(N) REAL Y0, YMED LOGICAL CHANGE CHANGE = .FALSE. Y0 = 3.0 * Y(2) - 2.0 * Y(3) CALL MEDOF3(Y0, Y(1), Y(2), YMED, CHANGE) Y(1) = YMED C Y0 = 3.0 * Y(N-1) - 2.0 * Y(N-2) CALL MEDOF3(Y0, Y(N), Y(N-1), YMED, CHANGE) Y(N) = YMED 999 RETURN END C SUBROUTINE SPLIT(Y, N, CHANGE) INTEGER N REAL Y(N) LOGICAL CHANGE REAL W(6), Y1 INTEGER I1, I, NM2 C NM2 = N - 2 DO 10 I=1,4 W(I + 2) = Y(I) 10 CONTINUE W(2) = Y(3) I1 = 1 20 IF( W(3) .NE. W(4) ) GO TO 40 IF( (W(3) - W(2)) * (W(5) - W(4)) .GE. 0.0 ) GO TO 40 IF( I1 .LT. 3) GO TO 30 C Y1 = 3.0 * W(2) - 2.0 * W(1) CALL MEDOF3(Y1, W(3), W(2), Y(I1), CHANGE) 30 IF(I1 .GE. NM2) GO TO 40 Y1 = 3.0 * W(5) - 2.0 * W(6) CALL MEDOF3(Y1, W(4), W(5), Y(I1+1), CHANGE) 40 DO 50 I=1,5 W(I) = W(I+1) 50 CONTINUE I1 = I1 + 1 IF(I1 .GE. NM2) GO TO 60 W(6) = Y(I1 + 3) GO TO 20 60 W(6) = W(3) IF( I1 .LT. N) GO TO 20 999 RETURN END C SUBROUTINE RUNMED(Y, N, LEN, WORK, SAVE, NW, ERR) INTEGER N, LEN, NW, ERR REAL Y(N) , WORK(NW), SAVE(NW) REAL MEDIAN REAL TEMP, TWO INTEGER SAVEPT, SMOPT, LENP1, I, J DATA TWO / 2.0/ IF(LEN .LE. NW) GO TO 5 ERR = 62 GO TO 999 5 DO 10 I=1, LEN WORK(I) = Y(I) SAVE(I) = Y(I) 10 CONTINUE SAVEPT = 1 SMOPT = INT( (FLOAT(LEN) + TWO) / TWO) LENP1 = LEN + 1 DO 50 I= LENP1, N CALL SORT(WORK, LEN, ERR) IF(ERR .NE. 0) GO TO 999 Y(SMOPT) = MEDIAN(WORK, LEN) TEMP = SAVE(SAVEPT) DO 20 J=1, LEN IF(WORK(J) .EQ. TEMP) GO TO 30 20 CONTINUE ERR = 63 GO TO 999 30 WORK(J) = Y(I) SAVE(SAVEPT) = Y(I) SAVEPT = MOD(SAVEPT, LEN ) + 1 SMOPT = SMOPT + 1 50 CONTINUE CALL SORT(WORK, LEN, ERR) IF(ERR .NE. 0) GO TO 999 Y(SMOPT) = MEDIAN(WORK, LEN) 999 RETURN END C REAL FUNCTION MEDIAN(Y, N) C FIND THE MEDIAN OF THE SORTED VALUES Y(1), Y(2),...,Y(N) INTEGER N REAL Y(N) INTEGER MPTR, MPT2 C MPTR =(N/2) +1 MPT2 = N - MPTR + 1 MEDIAN = ( Y(MPTR) + Y(MPT2) ) / 2.0 RETURN END C SUBROUTINE SORT(Y, N, ERR) INTEGER N, ERR REAL Y(N) INTEGER I, J, J1, GAP, NMG REAL TEMP IF(N .GE. 1) GO TO 10 ERR = 1 GO TO 999 10 IF(N .EQ. 1) GO TO 999 GAP = N 20 GAP = GAP / 2 NMG = N - GAP DO 40 J1 = 1, NMG I = J1 + GAP J = J1 30 IF( Y(J) .LE. Y(I)) GO TO 40 TEMP = Y(I) Y(I) = Y(J) Y(J) = TEMP I = J J = J - GAP IF ( J .GE. 1) GO TO 30 40 CONTINUE IF(GAP .GT. 1) GO TO 20 999 RETURN END SUBROUTINE CORREL(V1,V2,NDATA,COEFF) C------------------------------------------------------------ C COMPUTES CORRELATION COEFFICIENT OF TWO VECTORS V1 AND V2 C------------------------------------------------------------ IMPLICIT REAL*4(A-H,O-Z) DIMENSION V1(2000),V2(2000) SUMNUM = 0.00 SUMD1 = 0.00 SUMD2 = 0.00 CALL XMEAN(V1,NDATA,XBAR1) CALL XMEAN(V2,NDATA,XBAR2) DO 1 I=1,NDATA SUMNUM = SUMNUM + (V1(I) - XBAR1)*(V2(I) - XBAR2) SUMD1 = SUMD1 + (V1(I) - XBAR1)**2 SUMD2 = SUMD2 + (V2(I) - XBAR2)**2 1 CONTINUE COEFF = SUMNUM/(SQRT(SUMD1*SUMD2)) RETURN END C SUBROUTINE ZTEST(ARRAY, N, ALFA, R, ZVALUE, ZCRIT) IMPLICIT REAL * 4(A-H,O-Z) DIMENSION ARRAY(N), Y(2000) C-------------------------------------------------------------------- C PERFORMS AN ENTROPY Z-TEST OF A DATA SET C ARRAY = VECTOR OF DATA (INPUT) C N = NUMBER OF DATA (INPUT) C ALFA = PROBABLITY VALUE FOR COMPUTING THE CRITICAL Z (INPUT) C R = CORRELATION COEFFICIENT C ZVALUE= COMPUTED TEST STATISTIC ( MUST FOLLOW N(0,N/3) ) C ZCRIT = COMPUTED CRITICAL VALUE AT ALFA PROBABILITY (PERCENTILE) C C AUTHOR : S. MERTIKAS C DATE : JULY, 1986 C RE : BIOMETRIKA, VOL.67, NO.2, PP.455-461 (1980) C-------------------------------------------------------------------- C COMPUTE Y'S C----------------------------------------------- EN = FLOAT( N ) DO 10 I = 1, N SUM1 = 0.0 SUM2 = 0.0 DO 20 IA = 1, N IF(I .EQ. IA) GO TO 20 SUM1 = SUM1 + ARRAY(IA) ** 2 SUM2 = SUM2 + ARRAY(IA) 20 CONTINUE A = SUM1 B = (SUM2 ** 2) / (EN - 1.0) Y(I) = (A - B) ** (1.0 / 3.0) Y(I) = Y(I) / EN 10 CONTINUE C---------------------------------------------- C COMPUTE CORRELATION COEFFICIENT C---------------------------------------------- CALL CORREL(ARRAY,Y,N,R) C---------------------------------------------- C COMPUTE TEST STATISTIC C---------------------------------------------- ZVALUE = 0.5 * ALOG10( (1.0 + R)/(1.0 - R) ) C--------------------------------------------- C COMPUTE THE UPPER 100*ALFA PERCENTILE C FOR THE NULL DISTRIBUTION OF Z STATISTIC C--------------------------------------------- SIGMA = 3.0 / (EN) - 7.324 / (EN ** 2) + 53.005 /(EN ** 3) SIGMA = SQRT(SIGMA) C GAMMA = -11.70 /( EN ) + 55.06 /( EN ** 2 ) C CALL MDNRIS(ALFA, U, IER) ZCRIT = SIGMA * ( (U + ( U ** 3 - 3.0 * U ) * GAMMA) / 24.) WRITE(6,*) ALFA, U, SIGMA, GAMMA , EN RETURN END C SUBROUTINE AGOSTN(ARRAY, N, DA, YA) IMPLICIT REAL * 4 (A-H,O-Z) C------------------------------------------------------------------- C COMPUTES D'AGOSTINO STAISTICAL TEST FOR NORMALITY DEPARTURES C IT IS AN EXTENSION OF THE SHAPIRO-WILK TEST C ARRAY = VECTOR OF ORDERED OBSERVATIONS (INPUT) C N = NUMBER OF ORDERED OBSERVATIONS (INPUT) C DA = D'AGOSTINO TEST STATISTIC (OUTPUT) C YA = APPROXIMATE STANDARDIZED VERSION OF D (OUTPUT) C C DATA : JULY, 1986 C AUTHOR : S.MERTIKAS C RE : BIOMETRIKA, VOL.58, PP.341-348,(1971) C------------------------------------------------------------------ DIMENSION ARRAY(N) PI = 4.0 * ATAN(1.0) SUM1 = 0.0 SUM2 = 0.0 DO 10 I=1,N B = FLOAT(I) - ((FLOAT(N) + 1.0)/2.0) SUM1 = SUM1 + B * ARRAY(I) 10 CONTINUE CALL XMEAN(ARRAY, N, XBAR) C DO 20 I=1,N SUM2 = SUM2 + (ARRAY(I) - XBAR) ** 2 20 CONTINUE S = SQRT(SUM2) C-------------------------------------------------- C TEST STATISTIC (D) C-------------------------------------------------- F = FLOAT(N) ** (3.0/2.0) DA = SUM1 / ( F * S ) C------------------------------------------------- C STANDARDIZED STATISTIC C------------------------------------------------- SQRN = SQRT(FLOAT(N)) DMEAN = 1.0 / ( 2.0 * SQRT(PI) ) DNUM = (SQRN * DA) - (SQRN * DMEAN) YA = DNUM / 0.02998598 WRITE(6,1000) WRITE(6,2000) DA, YA 1000 FORMAT(20X,'D-AGOSTINO STATISTICAL TEST FOR NORMALITY #DEPARTURE',//) 2000 FORMAT(5X,'CALCULATED TEST STATISTIC D=',F10.4, #2X,'STANDARDIZED VALUE=',F10.4) RETURN END C SUBROUTINE ANDRSN(ARRAY, N, AVALUE, ACRT) C---------------------------------------------------------------- C COMPUTES THE ANDERSON-DARLING A**2 VALUE FOR NORMALITY TESTING C ARRAY(N) = VECTOR OF N ELEMENTS (ORDERED SAMPLE) C AVALUE = COMPUTED A**2 VALUE C ACRT = CRITICAL VALUE WITH 95% CONFIDENCE. C C AUTHOR : S.MERTIKAS C DATE : JULY, 1986 C RE : APPLIED STATISTICS, VOL.26, NO.2, PP.156,1977 C EXTERN : MDNOR FROM IMSL LIBRARY C--------------------------------------------------------------------- IMPLICIT REAL * 4 (A-H,O-Z) DIMENSION ARRAY(N) C---------------------------------------------------- C COMPUTE MEAN AND STANDARD DEVIATION C---------------------------------------------------- SUM1 = 0.0 SUM2 = 0.0 CALL XMEAN(ARRAY, N, XBAR) DO 10 I=1, N SUM2 = SUM2 +( ARRAY(I) -XBAR) **2 10 CONTINUE SCALE = SQRT( SUM2/ (FLOAT(N) - 1.0) ) C----------------------------------------------- C COMPUTE TEST STATISTIC (ANDERSON-DARLING) C----------------------------------------------- DO 20 I=1,N J = N + 1 - I Z1 = (ARRAY(I) -XBAR)/ SCALE Z2 = (ARRAY(J) -XBAR)/ SCALE CALL MDNOR( Z1, PROB1) CALL MDNOR( Z2, PROB2) ZLN = ALOG( PROB1) + ALOG(1.0 - PROB2) COFF = 2.0 * FLOAT(I) - 1.0 SUM1 = SUM1 + COFF * ZLN 20 CONTINUE C AVALUE = - FLOAT(N) - (SUM1 / FLOAT(N) ) C----------------------------------------------------------------- C COMPUTE CRITICAL VALUE FOR P=95.0% PR{ AVALUE < ACRT}= 95%. C------------------------------------------------------------------ A00 = 0.7514 B1 =-0.89 B0 =-0.795 ACRT = A00 * ( 1.0 + B0 / (FLOAT(N)) + B1/(FLOAT(N)**2) ) WRITE(6,1000) WRITE(6,2000) AVALUE , ACRT 1000 FORMAT(20X,'ANDERSON DARLING STATISTICAL TEST FOR NORMALITY # DEPARTURES',//) 2000 FORMAT(4X,'COMPUTED TEST STATISTIC A**2=',F10.4,2X,'CRITICAL #VALUE =',F10.4,'AT 95% C.L.') RETURN END C SUBROUTINE XMEAN(ARRAY,NDATA,XBAR) IMPLICIT REAL*4(A-H,O-Z) DIMENSION ARRAY(NDATA) XBAR = 0.0 DO 1 I=1,NDATA 1 XBAR=XBAR+ARRAY(I) XBAR=XBAR/FLOAT(NDATA) RETURN END C SUBROUTINE PSEUDO(ARRAY, N, SPRD, SMIDS, PSIGMA, ZSQR,NSIGMA) IMPLICIT REAL * 4 (A-H,O-Z) DIMENSION ARRAY(N), SPRD(NSIGMA), SMIDS(NSIGMA), PSIGMA(NSIGMA) DIMENSION ZSQR(NSIGMA) C------------------------------------------------------------------- C COMPUTES PSEUDOSIGMA VALUES FROM ORDERED DATA C ARRAY = VECTOR OF ORDERED DATA C N = NUMBER OF ORDERED DATA C SPRD = VECTOR OF SPREAD VALUES FROM DATA C SMIDS = VECTOR OF MIDSUMMARIES C PSIGMA= VECTOR OF PSEUDOSIGMA VALUES, GAUSSIAN REFFERENCE C ZSQR = SQUARED STANDARD GAUSSIAN DEVIATE C NSIGMA= NUMBER OF PSEUDOSIGMA VALUES, USUALLY N/2. C C EXTERNAL : MDNRIS (FROM IMSL LIBRARY) C AUTHOR : S.MERTIKAS (JULY,1986) C------------------------------------------------------------------- DO 10 I=1, NSIGMA SPRD( I ) = ARRAY(N + 1 - I) - ARRAY(I) SMIDS(I) =(ARRAY(N + 1 - I) + ARRAY(I)) / 2.0 PROB = (I - (1.0/3.0) ) / (FLOAT(N) + (1.0/3.0) ) CALL MDNRIS(PROB,QNRMAL,IER) ZSQR(I) = QNRMAL ** 2 PSIGMA(I) = SPRD(I) / (2.0 * ABS(QNRMAL)) 10 CONTINUE RETURN END C SUBROUTINE FUNOP(ARRAY, N, QUANT, Z) IMPLICIT REAL * 4 (A-H, O-Z) DIMENSION ARRAY(N), QUANT(N), Z(N) C----------------------------------------------------------------- C ARRAY = VECTOR OF ORDERED DATA WITH MEDIAN = 0. C N = NUMBER OF DATA C QUANT = STANDARD GAUSSIAN DEVIATE CORRESPONDING TO A SELECTED C PROBABILITY (E.G., (3*I - 1)/(3*N +1),I=1,2,.,N ). C Z = VECTOR OF ESTIMATED SLOPE (SIGMAS) C----------------------------------------------------------------- MIN = N / 3.0 MAX = 2.0 * (N / 3.0) DO 10 I =1, MIN Z(I) = ARRAY(I) / QUANT(I) 10 CONTINUE C DO 20 I= MAX, N Z(I) = ARRAY(I) / QUANT(I) 20 CONTINUE RETURN END C SUBROUTINE SETCON(ARRAY, N, IFIRST, ILAST, CONST) C------------------------------------------------------------ C FUNCTION : SETS ARRAY TO A CONSTANT C C HISTORY : VERSION DATE PROGRAMMER C 1 JUNE 1984 AL. KLEUSBERG C 2 FEB 1986 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IARR I/O VECTOR C N I DIMENSION OF VECTOR C CONST I CONSTANT C C EXTERNALS : NONE C------------------------------------------------------------- IMPLICIT REAL*4(A-H,O-Z) DIMENSION ARRAY(N) DO 10 I=IFIRST, ILAST ARRAY(I) = CONST 10 CONTINUE RETURN END C SUBROUTINE SET(ARRAY, N, IFIRST, ILAST, ARRAY1, NDATA) INTEGER NDATA, I, IFIRST, ILAST, N DIMENSION ARRAY( N ), ARRAY1(NDATA) C J = 0 DO 1 I = IFIRST, ILAST J = J + 1 ARRAY( I ) = ARRAY1( J ) 1 CONTINUE RETURN END C SUBROUTINE CLEAR(ARRAY, N, IOUTH, IOUTL, CLEAN) C----------------------------------------------------------------- C CREATES A NEW VECTOR CLEAN(.) WHICH HAS THE FAR OUT VALUES C REMOVED FROM IT. C AUTHOR : S. MERTIKAS C DATE : JUNE 1986 C----------------------------------------------------------------- INTEGER IOUTH, IOUTL INTEGER I, J, N, MIN, MAX REAL ARRAY(N), CLEAN(N) J = 0 MIN = IOUTL + 1 MAX = IOUTH - 1 DO 1 I = MIN, MAX J = J + 1 CLEAN( J ) = ARRAY(I) 1 CONTINUE RETURN END C SUBROUTINE OUT(ARRAY,NDATA,IOUTL,IOUTH, IADJL, IADJH, GDOP, INDX) C--------------------------------------------------------------------- C PRINTS OUT FAR OUT AND OUTSIDE VALUES C--------------------------------------------------------------------- INTEGER INDX(NDATA) INTEGER I, NDATA, IOUTL, IOUTH, IADJL, IADJH REAL ARRAY(NDATA), GDOP(NDATA) C WRITE(6,1012) DO 5 I = 1, IOUTL WRITE(6,100) ARRAY(I), I, GDOP(INDX(I)) 5 CONTINUE C WRITE(6,1010) DO 4 I = IOUTL+ 1, IADJL WRITE(6,100) ARRAY(I), I, GDOP(INDX(I)) 4 CONTINUE C WRITE(6,1013) DO 7 I=NDATA ,IOUTH, -1 WRITE(6,100) ARRAY(I), I, GDOP(INDX(I)) 7 CONTINUE C WRITE(6,1011) DO 6 I=IOUTH- 1 ,IADJH, -1 WRITE(6,100) ARRAY(I), I, GDOP(INDX(I)) 6 CONTINUE 100 FORMAT(' ',4X,F10.3,2X,I7,2X,F10.3) 1012 FORMAT('1',4X,'----------FAR OUT VALUES (LOW)--------------',//) 1010 FORMAT(4X,'---------LOW VALUES BETWEEN INNER AND OUTER FENCE #-------------------------------------',//) 1013 FORMAT('1',4X,'##########FAR OUT VALUES (HIGH)#############',//) 1011 FORMAT(4X,'##########HIGH VALUES BETWEEN INNER AND OUTER # FENCE ##############################',//) RETURN END C SUBROUTINE PSORT(ON, WITH, N, ERR) INTEGER N, ERR REAL ON(N), WITH(N) C-------------------------------------------------------------------- C PAIR SHELL SORT N VALUES IN ON(N) FROM SMALLEST TO LARGEST C CARRYING ALONG THE VALUES IN WITH(N) C-------------------------------------------------------------------- INTEGER I, J, J1, GAP, NMG REAL TON, TWITH C IF(N .GE. 1) GO TO 10 ERR = 1 GO TO 999 10 IF(N .EQ. 1) GO TO 999 C------------------------------------------- C ONE ELEMENT IS ALWAYS SORTED C------------------------------------------- GAP = N 20 GAP = GAP / 2 NMG = N - GAP DO 40 J1= 1, NMG I = J1 + GAP J = J1 30 IF(ON(J) .LE. ON(I)) GO TO 40 C------------------------------------- C SWAP OUT OF ORDER PAIR C------------------------------------- TON = ON(I) ON(I) = ON(J) ON(J) = TON TWITH = WITH(I) WITH(I) = WITH(J) WITH(J) = TWITH C KEEP OLD POINTER FOR NEXT TIME THROUGH I = J J = J - GAP IF(J .GE. 1) GO TO 30 40 CONTINUE IF(GAP .GT. 1 ) GO TO 20 999 RETURN END C SUBROUTINE NPOSW(HI, LO, NICNOS, NN, MAXP, MZERO, PTOTL, FRACT, # UNIT, NPW, ERR) C-------------------------------------------------------------------- C FIND A NICE(SIMPLE) DATA-UNITS VALUE TO ASSIGN TO ONE PLOT C POSITION IN ONE DIMENSION OF A PLOT. A PLOTPOSITION IS C TYPICALLY ONE CHARACTER POSITION HORIZONTALLY, OR ONE LINE C VERTICALLY. C C INPUT: C HI, LO = THE HIGH AND LOW EDGES OF THE DATA RANGE TO BE PLOTTED C NICNOS = A VECTOR LENGTH NN CONTAINING NICE MANTISSAS FOR PLOT C MAXP = MAXIMUM NUMBER OF PLOT POSITIONS ALLOWED IN THIS DIMENSION C MZERO = .TRUE. IF A POSITION LABELED C = -0 AS ALLOWED IN THIS DIMENSION C = .FALSE. OTHERWISE C C OUTPUT: C PTOTL = HOLDS THE TOTAL NUMBER OF PLOT POSITIONS TO BE USED IN THIS C DIMENSION (MUST BE .LE. MAXP) C FRACT = MANTISSA OF THE NICE POSITION WIDTH, SELECTED FROM THE C NUMBERS IN NICNOS. C UNIT = INTEGER POWER OF 10 SUCH THAT NPW= FRACT * UNIT. C NPW = NICE POSITION WIDTH. ONE PLOT POSITION WIDTH WILL REPRESENT C A DATA-SPACE DISTANCE OF NPW. C-------------------------------------------------------------------- C INTEGER NN, MAXP, PTOTL, ERR REAL HI, LO, NICNOS(NN), FRACT, UNIT, NPW LOGICAL MZERO C------------------------------ C FUNCTIONS C------------------------------ INTEGER FLOOR, INTFN C------------------------------ C LOCAL VARIABLES C------------------------------ INTEGER I REAL APRXW IF(MAXP .GT. 0) GO TO 5 ERR = 8 GO TO 999 5 APRXW = (HI - LO) / FLOAT(MAXP) IF(APRXW .GT. 0.0) GO TO 10 C------------------------------ C IF(HI .LE. LO) IS AN ERROR C------------------------------ ERR = 9 GO TO 999 10 UNIT = 10.0 ** FLOOR(ALOG10(APRXW)) FRACT = APRXW / UNIT DO 20 I=1, NN IF( FRACT .LE. NICNOS(I) ) GO TO 30 20 CONTINUE 30 FRACT = NICNOS(I) NPW = FRACT * UNIT PTOTL = INTFN(HI / NPW, ERR) - INTFN(LO / NPW, ERR) + 1 IF( ERR .NE. 0) GO TO 999 C-------------------------------- C IF MINUS ZERO POSITION POSSIBLE AND( SGN(HI) .NE. SGN(LO) C ALLOW IT C-------------------------------- IF( MZERO .AND. (HI*LO .LT. 0.0 .OR. HI .EQ. 0.0)) PTOTL=PTOTL+1 C-------------------------------- C PLOT POSITION REQUIRED WITH THIS WIDTH C-------------------------------- IF( PTOTL .LE. MAXP) GO TO 999 C-------------------------------- C TOO MANY POSITIONS NEEDED, SO BUMP NPW UP ONE NICE NUMBER C-------------------------------- I = I + 1 IF( I .LE. NN) GO TO 30 I = 1 UNIT = UNIT * 10.0 GO TO 30 999 RETURN END C INTEGER FUNCTION FLOOR(Y) REAL Y C----------------------------------------------- C FIND FLOOR(Y), THE LARGEST INTEGER NOT EXCEEDING Y C----------------------------------------------- FLOOR = INT(Y) IF(Y .LT. 0.0 .AND. Y .NE. FLOAT(FLOOR)) FLOOR= FLOOR - 1 RETURN END C SUBROUTINE BOXTOP(LO, HI, HL, HH, NPW, ERR) REAL LO, HI, HL, HH, NPW INTEGER ERR C------------------------------------------------- C PRINT THE TOP OR BOTTOM OF A BOXPLOT DISPLAY C HI AND LO ARE EDGES OF THE PLOTTING REGION C USED BY THE PLTPOS FUNCTION C HL AND HH ARE THE LOW AND HIGH HINGES C NPW= NICE POSITION WIDTH SET BY THE PLOT SCALING ROUTINES C C LOCAL VARIABLES C-------------------------------------------------- INTEGER I, IFROM, ITO, CHMIN C--------------------------------------------------- C FUNCTION C--------------------------------------------------- INTEGER PLTPOS DATA CHMIN/40/ IFROM = PLTPOS(HL, LO, HI, NPW, ERR) ITO = PLTPOS(HH, LO, HI, NPW, ERR) IF(IFROM .GT. ITO) GO TO 11 DO 10 I= IFROM, ITO CALL PUTCHR(I, CHMIN, ERR) 10 CONTINUE 11 CONTINUE IF( ERR .EQ. 0) CALL PRINT 999 RETURN END C INTEGER FUNCTION PLTPOS(X, LO, HI, NPW, ERR) C-------------------------------------------------- C FIND THE POSITION CORRESPONDING TO X ON PLOT C BOUNDED BETWEEN LO AND HI AND SCALED ACCORDING TO NPW C--------------------------------------------------- REAL X, LO, HI, NPW INTEGER ERR INTEGER INTFN C----------------------- C COMMON VARIABLES C----------------------- COMMON /CHRBUF/P, PMAX, PMIN, OUTPTR, MAXPTR, OUNIT INTEGER P(130), PMAX, PMIN, OUTPTR, MAXPTR, OUNIT C PLTPOS = INTFN((X -LO)/NPW, ERR) + PMIN IF(PLTPOS .LT. PMIN) PLTPOS = PMIN IF(PLTPOS .GT. PMAX) PLTPOS = PMAX RETURN END C SUBROUTINE PUTCHR(POSN, CHAR, ERR) INTEGER POSN, CHAR, ERR C----------------------------------------------------------------- C PLACE THE CHARACTER CHAR AT POSITION POSN IN THE OUTPUT C LINE P. IF POSN=0 PLACE CHAR IN THE NEXT AVAILABLE POSITION C IN P. MAXPTR IS TO BE INITIALIZED TO PMIN AND PRINT MUST RESET IT C----------------------------------------------------------------- COMMON /CHARIO/ CHARS, CMAX, 1 CHA, CHB, CHC, CHD, CHE, CHF, CHG, CHH, CHI, CHJ, CHK, 2 CHL, CHM, CHN, CHO, CHP, CHQ, CHR, CHS, CHT, CHU, CHV, 3 CHW, CHX, CHY, CHZ, CH0, CH1, CH2, CH3, CH4, CH5, CH6, 4 CH7, CH8, CH9, CHBL, CHEQ, CHPLUS, CHMIN, CHSTAR, CHSLSH, 5 CHLPAR, CHRPAR, CHCOMA, CHPT C COMMON /CHRBUF/P, PMAX, PMIN, OUTPTR, MAXPTR,OUNIT C INTEGER CHARS(46), CMAX INTEGER CHA, CHB, CHC, CHD, CHE, CHF, CHG, CHH, CHI INTEGER CHJ, CHK, CHL, CHM, CHN, CHO, CHP, CHQ, CHR INTEGER CHS, CHT, CHU, CHV, CHW, CHX, CHY, CHZ INTEGER CH0, CH1, CH2, CH3, CH4, CH5, CH6, CH7, CH8, CH9 INTEGER CHBL, CHEQ, CHPLUS, CHMIN, CHSTAR, CHSLSH INTEGER CHLPAR, CHRPAR, CHCOMA, PHPT INTEGER P(130), PMAX, PMIN, OUTPTR, MAXPTR, OUNIT C IF(CHAR .GT. 0 .AND. CHAR .LE. CMAX) GO TO 10 ERR = 4 RETURN 10 IF(POSN .NE. 0) OUTPTR = MAX0(PMIN, POSN) OUTPTR = MIN0(OUTPTR, PMAX) P(OUTPTR) = CHARS(CHAR) MAXPTR = MAX0(MAXPTR, OUTPTR) OUTPTR = OUTPTR + 1 RETURN END C SUBROUTINE PRINT C--------------------------------------------------------------- C PRINT THE OUTPUT LINE P ON UNIT OUNIT C MAXPTR INDICATES THE RIGHTMOST POSITION WHICH HAS BEEN USED IN C LINE. THEN, RESET P TO SPACES AND MAXPTR AND OUTPTR TO PMIN. C--------------------------------------------------------------- COMMON /CHRBUF/ P, PMAX, PMIN, OUTPTR, MAXPTR, OUNIT INTEGER P(130), PMAX, PMIN, OUTPTR, MAXPTR, OUNIT C--------------------------------- C LOCAL VARIABLES C--------------------------------- INTEGER BLANK, I DATA BLANK /1H / C WRITE(OUNIT, 10) (P(I), I=1, MAXPTR) 10 FORMAT(1X, 130A1) C DO 20 I=1, MAXPTR P(I) = BLANK 20 CONTINUE C OUTPTR = PMIN MAXPTR = PMIN RETURN END C SUBROUTINE YINFO(Y, N, MED, HL, HH, ADJL, ADJH, IADJL, IADJH, # STEP,IOUTL,IOUTH,ERR) C----------------------------------------------------------------- C GET GENERAL INFORMATION ABOUT Y(). USEFUL FOR PLOT SCALING C SORTS Y() AND RETURNS IT SORTED. IT ALSO RETURNS: C MED= MEDIAN C HL = LOW HINGE HH = HIGH HINGE C ADJL= LOW ADJACENT VALUE, ADJH = HIGH ADJACENT VALUE C IADJL = ITS INDEX(LOCATION) IADJH = ITS INDEX C---------------------------------------------------------------- INTEGER N, IADJL, IADJH, IOUTL, IOUTH, ERR REAL Y(N), MED, HL, HH, ADJL, ADJH, STEP C--------------------------------- C LOCAL VARIABLES C--------------------------------- REAL HFENCE, LFENCE, OUTFL, OUTFH INTEGER J, K, TEMP1, TEMP2 C CALL SORT(Y, N, ERR) IF(ERR .NE. 0) GO TO 999 K = N J = (K/2) + 1 C TEMP1 = N + 1 - J MED = (Y(J) + Y(TEMP1)) / 2.0 C K = (K + 1)/2 J = (K /2) + 1 TEMP1 = K + 1 - J HL = (Y(J) + Y(TEMP1))/ 2.0 TEMP1= N - K + J TEMP2 = N + 1 -J HH = (Y(TEMP1) + Y(TEMP2)) / 2.0 C STEP = (HH - HL) * 1.5 HFENCE = HH + STEP LFENCE = HL - STEP OUTFH = HH + 2 * STEP OUTFL = HL - 2 * STEP C-------------------------------------- C FIND ADJACENT VALUES AND WRITE OUTSIDE VALUES C-------------------------------------- IADJL = 0 20 IADJL = IADJL + 1 IF( Y(IADJL) .LE. OUTFL) IOUTL = IADJL IF( Y(IADJL) .LE. LFENCE) GO TO 20 ADJL = Y(IADJL) C IADJH = N + 1 30 IADJH = IADJH - 1 IF( Y(IADJH) .GE. OUTFH) IOUTH = IADJH IF( Y(IADJH) .GE. HFENCE) GO TO 30 ADJH = Y(IADJH) 999 RETURN END C INTEGER FUNCTION INTFN(X, ERR) C-------------------------------------------------------------------- C FIND THE INTEGER EQUAL TO OR NEXT CLOSER TO ZERO THAN X C CHECKS TO SEE THAT X IS NOT TOO LARGE TO FIT IN AN INTEGER VARIABL C-------------------------------------------------------------------- REAL X INTEGER ERR C COMMON/NUMBRS/ EPSI, MAXINT REAL EPSI, MAXINT IF( ABS(X) .LE. MAXINT) GO TO 10 C--------------------------------------------------------- C X IS TOO LARGE IN MAGNITUDE TO FIT IN AN INTEGER C RETURN THE LARGEST LEGAL INTEGER AND SET THE ERROR FLAG C--------------------------------------------------------- ERR = 3 INTFN = IFIX( SIGN(MAXINT, X) ) GO TO 999 10 INTFN = INT ( (1.0 + EPSI )* X) 999 RETURN END C SUBROUTINE BOXES(Y, N, GSUB, NG, LINE3, NOTCH, SORTY, ERR) C--------------------------------------------------------------- C PRINTS ADJACENT BOXPLOTS ON A SINGLE SCALE FOR ALL VARIABLES IN Y() C--------------------------------------------------------------- INTEGER N, NG, ERR INTEGER GSUB(N) REAL Y(N), SORTY(N) LOGICAL LINE3, NOTCH C--------------------------------------------------------------- C INPUT : C Y() = VECTOR THAT CONTAINS ALL DATA C GSUB()= CONTAINS INTEGERS BETWEEN 1 TO NG, IDENTIFYING THE C DATA SET EACH ELEMENT OF Y() BELONGS TO. C E.G., (GSUB(I)=1,I=1,IFIRST), (GSUB(I)=2,I=1,SECOND).. C THE USE OF THE GSUB() VECTOR IS TO SUGGEST BOXPLOTS C OF EITHER THE ROWS OR THE COLUMNS OF A MATRIX. C NG = NUMBER OF GROUPS. C IF LINE3 IS .TRUE. A 3-LINE BOXPLOT (FULL BOXES) IS PRINTED. C IF NOT, THE SIMPLE ONE-LINE BOX PLOT IS PRINTED C BOTH CONVEY THE SAME INFORMATION, BUT THE 3-LINE VERSION C MAY LOOK NICER. C IF NOTCH IS .TRUE. A CONFIDENCE INTERVAL AROUND THE MEDIAN C IS INDICATED WITH PARENTHESES. C---------------------------------------------------------------- COMMON/CHRBUF/P, PMAX, PMIN, OUTPTR, MAXPTR, OUNIT INTEGER P(130), PMAX, PMIN, OUTPTR, MAXPTR, OUNIT C------------------------------------------------------ C LOCAL VARIABLES C----------------------------------------------------- INTEGER NN, NPMAX, NPOS, LPMIN, SPMIN INTEGER CHRPAR, LBLW, OPOS, I, J, K REAL NICNOS(9), FRACT, UNIT, NPW, LO, HI C------------------------------------------------------- C FUNCTIONS C------------------------------------------------------- INTEGER WDTHOF C------------------------------------------------------- C CALL SUBROUTINES BOXP, NPOSW, PUTCHR, PUTNUM C------------------------------------------------------- DATA NN, NICNOS(1),NICNOS(2),NICNOS(3)/9,1.0,1.5,2.0/ DATA NICNOS(4),NICNOS(5),NICNOS(6)/2.5,3.0,4.0/ DATA NICNOS(7),NICNOS(8),NICNOS(9)/5.0,7.0,10.0/ DATA CHRPAR/44/ C------------------------------------------------------------ C CHECK FOR AT LEAST 2 DATA VALUES. OTHERWISE HIGHEST AND C LOWEST WILL BE EQUAL AND PLOT SCALING WILL FAIL C----------------------------------------------------------- IF(N .GT. 1) GO TO 5 ERR = 31 GO TO 999 5 LPMIN = PMIN + 7 LO = Y(1) HI = Y(N) DO 10 I=1,N IF(LO .GT. Y(I)) LO=Y(I) IF(HI .LT. Y(I)) HI=Y(I) 10 CONTINUE C--------------------------------------------------------- C SCALE THE EXTREMES C--------------------------------------------------------- NPMAX = PMAX - LPMIN + 1 CALL NPOSW(HI, LO, NICNOS, NN, NPMAX,.FALSE., NPOS, FRACT, # UNIT,NPW, ERR) IF(ERR .NE. 0) GO TO 999 C---------------------------------------------------------- C NOW PRINT ALL THE BOXES C DATA SETS ARE IDENTIFIED BY THEIR CODES IN GSUB() C----------------------------------------------------------- IF(NG .GT. 1) GO TO 17 DO 15 K=1, N SORTY(K) = Y(K) 15 CONTINUE CALL BOXP(SORTY, N, LINE3, NOTCH, LO, HI, NPW, ERR) GO TO 999 17 SPMIN = PMIN DO 30 I=1, NG K = 0 WRITE(6,100) 100 FORMAT(///) DO 20 J= 1, N IF(GSUB(J) .NE. I) GO TO 20 K = K + 1 SORTY(K) = Y(J) 20 CONTINUE PMIN = SPMIN LBLW = WDTHOF(I) OPOS = PMIN + 5 - LBLW CALL PUTNUM(OPOS, I, LBLW, ERR) OPOS = PMIN + 6 CALL PUTCHR(OPOS,CHRPAR, ERR) IF(ERR .NE. 0) GO TO 999 PMIN = LPMIN CALL BOXP(SORTY, K, LINE3, NOTCH, LO, HI, NPW, ERR) IF(ERR .NE. 0) GO TO 999 30 CONTINUE PMIN = SPMIN 999 RETURN END C INTEGER FUNCTION WDTHOF(I) INTEGER I INTEGER IA, IQ, ND C IA = IABS(I) ND = 1 IF(I .LT. 0) ND = 2 10 IQ = IA/10 IF(IQ .EQ. 0) GO TO 20 IA = IQ ND = ND + 1 GO TO 10 20 WDTHOF = ND RETURN END C SUBROUTINE PUTNUM(POSN, N, W, ERR) INTEGER POSN, N, W, ERR C INTEGER CHD, CH0, CHMIN, DSTK(20), INUM, IP, IQ, IW, ND C COMMON/CHRBUF/P, PMAX, PMIN, OUTPTR, MAXPTR, OUNIT INTEGER P(130), PMAX, PMIN, OUTPTR, MAXPTR, OUNIT C DATA CH0, CHMIN /27, 40/ C IW = W IF(N .LT. 0) IW = IW - 1 INUM = IABS(N) C ND = 1 10 IQ = INUM / 10 DSTK(ND) = INUM - IQ * 10 IF(ND .LE. 20 .AND. ND .LE. IW) GO TO 20 ERR = 2 GO TO 999 20 IF(IQ .EQ. 0) GO TO 30 INUM = IQ ND = ND + 1 GO TO 10 C 30 IP = POSN IF(IP .EQ. 0) IP = OUTPTR IP = IP + IW -ND IF( N .GE. 0) GO TO 40 CALL PUTCHR(IP, CHMIN, ERR) IP = IP + 1 40 CHD = CH0 + DSTK(ND) CALL PUTCHR(IP, CHD, ERR) IF(ND .EQ. 1) GO TO 50 ND = ND -1 IP = IP + 1 GO TO 40 50 CONTINUE 999 RETURN END C SUBROUTINE BOXP(SORTY, N, LINE3, NOTCH, LO, HI, NPW, ERR) INTEGER N, ERR REAL SORTY(N), LO, HI, NPW LOGICAL LINE3, NOTCH C COMMON/CHRBUF/P, PMAX, PMIN, OUTPTR, MAXPTR, OUNIT INTEGER P(130), PMAX, PMIN, OUTPTR, MAXPTR, OUNIT C INTEGER PLTPOS C INTEGER I, IADJL, IADJH, IFROM, ITO, LPMAX, LPMIN INTEGER OPOS, CHI, CHO, CHSTAR, CHMIN, CHPLUS, CHRPAR, CHLPAR REAL MED, HL, HH, ADJL, ADJH, STEP REAL FLOATN, NSTEP, LNOTCH, HNOTCH, OFENCL, OFENCH C DATA CHI, CHO, CHPLUS, CHMIN, CHSTAR /9, 15, 39, 40, 41/ DATA CHLPAR, CHRPAR /43, 44/ C LPMAX = PMAX LPMIN = PMIN CALL YINFO(SORTY, N, MED, HL, HH, ADJL, ADJH, IADJL, IADJH, # STEP, IOUTL, IOUTH, ERR) IF(ERR .NE. 0) GO TO 999 FLOATN = FLOAT(N) NSTEP = 1.7 * (1.25 *(HH - HL)/(1.35 * SQRT(FLOATN))) LNOTCH = MED - NSTEP HNOTCH = MED + NSTEP C------------------------------------------------------------ C PRINT TOP OF BOX, IF 3-LINE VERSION C------------------------------------------------------------ IF(LINE3) CALL BOXTOP(LO, HI, HL, HH, NPW, ERR) IF(ERR .NE. 0) GO TO 999 C------------------------------------------------------------ C FILL CENTER LINE OF DISPLAY C NOTE CAREFUL HIERARCHY OF OVERPRINTING C LAST PLACED CHARACTER INS ONLY ONE TO APPEAR C C MARK WHISKERS C------------------------------------------------------------ IFROM = PLTPOS(ADJL, LO, HI, NPW, ERR) ITO = PLTPOS(HL, LO, HI, NPW, ERR) - 1 IF( IFROM .GT. ITO) GO TO 21 DO 20 I= IFROM, ITO CALL PUTCHR(I, CHMIN, ERR) 20 CONTINUE 21 CONTINUE IFROM = PLTPOS(HH, LO, HI, NPW, ERR) ITO = PLTPOS(ADJH, LO, HI, NPW, ERR) IF(IFROM .GT. ITO) GO TO 31 DO 30 I=IFROM, ITO CALL PUTCHR(I, CHMIN, ERR) 30 CONTINUE 31 CONTINUE C---------------------------------------------------------- C MARK LOW OUTLIERS IF ANY C---------------------------------------------------------- IF(IADJL .EQ. 1) GO TO 41 OFENCL = HL - 2.0 * STEP ITO = IADJL - 1 DO 40 I= 1, ITO OPOS = PLTPOS(SORTY(I), LO, HI, NPW, ERR) IF(SORTY(I) .LT. OFENCL) CALL PUTCHR(OPOS, CHO, ERR) IF(SORTY(I) .GE. OFENCL) CALL PUTCHR(OPOS, CHSTAR, ERR) 40 CONTINUE 41 CONTINUE C---------------------------------------------------------- C MARK HIGH OUTLIERS, IF ANY C---------------------------------------------------------- IF(IADJH .EQ. N) GO TO 51 OFENCH = HH + 2.0 * STEP IFROM = IADJH + 1 DO 50 I= IFROM, N OPOS = PLTPOS(SORTY(I), LO, HI, NPW, ERR) IF(SORTY(I) .GT. OFENCH) CALL PUTCHR(OPOS, CHO, ERR) IF(SORTY(I) .LE. OFENCH) CALL PUTCHR(OPOS, CHSTAR, ERR) 50 CONTINUE 51 CONTINUE C----------------------------------------------------------- C MARK HINGES, NOTCHES AND MEDIAN C----------------------------------------------------------- OPOS = PLTPOS(HL, LO, HI, NPW, ERR) CALL PUTCHR(OPOS, CHI, ERR) OPOS = PLTPOS(HH, LO, HI, NPW, ERR) CALL PUTCHR(OPOS, CHI, ERR) OPOS = PLTPOS(LNOTCH, LO, HI, NPW, ERR) IF(NOTCH) CALL PUTCHR(OPOS, CHLPAR, ERR) OPOS = PLTPOS(HNOTCH, LO, HI, NPW, ERR) IF(NOTCH) CALL PUTCHR(OPOS, CHRPAR, ERR) OPOS= PLTPOS(MED, LO, HI, NPW, ERR) CALL PUTCHR(OPOS, CHPLUS, ERR) C--------------------------------------------------------- C PRINT THE POX PLOT C-------------------------------------------------------- IF(ERR. NE. 0) GO TO 999 CALL PRINT C--------------------------------------------------------- C PRINT THE BOTTOM OF THE BOX C-------------------------------------------------------- IF(LINE3) CALL BOXTOP(LO, HI, HL, HH, NPW, ERR) 999 RETURN END C SUBROUTINE CINIT(IOUNIT, IPMIN, IPMAX, IEPSI, IMAXIN, ERR) INTEGER IOUNIT, IPMIN, IPMAX, IMAXIN, ERR REAL IEPSI C------------------------------------------------------------------- C INITIALIZATION. TO BE CALLED AT START OF ANY MAIN PROGRAM C WHICH CALLS ONE OF THE EDA SUBROUTINES C EITHER INDIRECTLY OR INDIRECTLY C C IOUNIT = OUTPUT LOGICAL UNIT C IPMIN = LEFT MARGIN C IPMAX = RIGHT MARGIN C IPEPSI = MACHINE-RELATED EPSILON C IMAXIN = MAXIMUM PERMITTED INTEGER VALUE C C ERR = ERROR FLAG (ROUTINE EXECUTED SUCCESSFULLY?) C---------------------------------------------------------------- COMMON/CHRBUF/P, PMAX, PMIN, OUTPTR, MAXPTR, OUNIT COMMON/NUMBRS/EPSI, MAXINT C INTEGER P(130), PMAX, PMIN, OUTPTR, MAXPTR, OUNIT REAL EPSI, MAXINT C------------------------------------------------ C LOCAL VARIABLES C------------------------------------------------ INTEGER BLANK, I DATA BLANK /1H / C ERR = 6 IF(IPMIN .LT. 1) GO TO 999 IF(IPMAX .GT. 130) GO TO 999 IF(IPMAX .LE. IPMIN) GO TO 999 ERR = 7 IF( (1.0 +IEPSI) .LE. 1.0) GO TO 999 ERR = 0 OUNIT = IOUNIT PMIN = IPMIN OUTPTR = IPMIN MAXPTR = IPMIN PMAX = IPMAX EPSI = IEPSI MAXINT = FLOAT(IMAXIN) C DO 50 I=1,130 P(I) = BLANK 50 CONTINUE 999 RETURN END C SUBROUTINE LVALS(Y, N, D, YLV, NLV, SORTY, ERR) INTEGER N, NLV, ERR REAL Y(N), D(15), YLV(15,2), SORTY(N) C------------------------------------------------------------ C FOR A BATCH OF VALUES IN Y C FIND THE SELECTED QUANTILES (LETTER VALUES) C YLV (OUT) LETTER VALUES C D (OUT) THE CORRESPONDING DEPTHS C NLV (OUT) THE NUMBER OF PAIRS OF LETTER VALUES C YLV(1,1) AND YLV(1,2) ARE SET EQUAL TO THE MEDIAN WHOSE C DEPTH D(1) IS (N+1)/2 C THE REST OF LETTER VALUES ARE COME IN PAIRS AND ARE STORED IN YLV C IN ORDER FROM THE HINGES OUT TO THE EXTREMES C THUS YLV(2,1) @ YLV(2,2) ARE THE LOWER AND UPPER HINGES AND C YLV(NLV,1) @ YLV(NLV,2) ARE THE LOWER (MIN) AND UPPER(MAX) EXTREMES C--------------------------------------------------------------------- C C LOCAL VARIABLES C------------------------------------------------ INTEGER I, J, K, PT1, PT2 IF((N .GT. 3) .AND. (N .LE. 24576)) GO TO 10 NLV = 0 ERR = 21 GO TO 999 C----------------------------------------------- C SORT Y INTO SORTY C----------------------------------------------- 10 DO 15 I=1,N SORTY(I) = Y(I) 15 CONTINUE CALL SORT(SORTY, N, ERR) IF(ERR .NE. 0) GO TO 999 C----------------------------------------------- C HANDLE MEDIAN SEPERATELY BECAUSE IT IS NOT C A PAIR OF LETTER VALUES. C----------------------------------------------- D(1) = FLOAT(N + 1) / 2.0 J = (N / 2) + 1 PT2 = N + 1 - J YLV(1,1) = (SORTY(J) + SORTY(PT2)) / 2.0 YLV(1,2) = YLV(1,1) C K = N I = 2 C 20 K = (K +1) / 2 J = (K / 2) + 1 D(I) = FLOAT(K + 1) / 2.0 C PT2 = K + 1 - J YLV(I,1) =(SORTY(J) + SORTY(PT2)) / 2.0 PT1 = N - K + J PT2 = N + 1 - J YLV(I,2) = (SORTY(PT1) + SORTY(PT2)) /2.0 C I = I+ 1 IF( D(I-1) .GT. 2.0) GO TO 20 NLV = I D(I) = 1.0 YLV(I,1) = SORTY(1) YLV(I,2) = SORTY(N) 999 RETURN END C SUBROUTINE LVPRNT(NLV, D, YLV, ERR) INTEGER NLV, ERR REAL D(15), YLV(15,2) REAL G(15) C----------------------------------------------------------------- C PRINTS A LETTER-VALUE DISPLAY C THE NLV PAIRS OF LETTER VALUES ARE IN YLV C YLV(I,1) = LOWER LETTER VALUE IN THE PAIR C YLV(I,2) = UPPER LETTER VALUE IN THE PAIR C YLV(1,1) = YLV(1,2) = MEDIAN (EXEPTION !!) C D = CORRESPONDING DEPTHS C--------------------------------------------------------------- COMMON /CHRBUF/ P, PMAX, PMIN, OUTPTR, MAXPTR, OUNIT INTEGER P(130), PMAX, PMIN, OUTPTR, MAXPTR, OUNIT C------------------------------------------------ C LOCAL VARIABLES C------------------------------------------------ INTEGER I, N, TAGS(14) REAL MID, SPR, SIGMA C DATA TAGS( 1), TAGS( 2), TAGS( 3), TAGS( 4) /1HM, 1HF, 1HE, 1HD/ DATA TAGS( 5), TAGS( 6), TAGS( 7), TAGS( 8) /1HC, 1HB, 1HA, 1HZ/ DATA TAGS( 9), TAGS(10), TAGS(11), TAGS(12) /1HY, 1HX, 1HW, 1HV/ DATA TAGS(13), TAGS(14) /1HU, 1HT/ C DATA G( 2), G( 3), G( 4), G( 5) /1.349, 2.301, 3.068, 3.725/ DATA G( 6), G( 7), G( 8), G( 9) /4.308, 4.835, 5.320, 5.771/ DATA G(10), G(11), G(12), G(13) /6.195, 6.594, 0.0 , 0.0 / C C PMAX = 100 C OUNIT= 6 IF((NLV .GE. 3) .AND. (NLV .LE.15)) GO TO 10 ERR = 22 GO TO 999 10 IF(PMAX .GE. 64) GO TO 20 ERR = 23 GO TO 999 20 WRITE(OUNIT, 1001) 1001 FORMAT(5X,5HDEPTH,7X,5HLOWER,8X,5HUPPER,11X, #8HMIDRANGE,5X,6HSPREAD,2X,11HPSEUDOSIGMA) C------------------------------------------------- C RECOVER N FROM D(1)=DEPTH OF THE MEDIAN C------------------------------------------------- N = INT(2.0 * D(1) ) - 1 WRITE(OUNIT, 1002) N 1002 FORMAT(1X,2HN=,I5) C------------------------------------------------ C WRITE LINE CONTAINING MEDIAN (AND FIRST MID) C------------------------------------------------ WRITE(OUNIT, 1003) D(1), YLV(1,1), YLV(1,1) 1003 FORMAT(1X,1HM,1X,F7.1,8X,F10.3,13X,F10.3) C N = NLV - 1 DO 30 I = 2, N MID = (YLV(I,1) + YLV(I,2) )/ 2.0 SPR = YLV(I,2) - YLV(I,1) SIGMA = SPR / G(I) WRITE(OUNIT, 1004) TAGS(I), D(I), YLV(I,1), # YLV(I,2), MID, SPR, SIGMA 1004 FORMAT(1X,A1,1X,F7.1,3X,F10.3,3X,F10.3,5X,F10.3,3X,F10.3, # 3X,F10.3) 30 CONTINUE MID = (YLV(NLV,1) + YLV(NLV,2)) / 2.0 SPR = YLV(NLV,2) - YLV(NLV,1) WRITE(OUNIT, 1005) YLV(NLV,1), YLV(NLV,2), MID, SPR 1005 FORMAT(7X,1H1,5X,F10.3,3X,F10.3,5X,F10.3,3X,F10.3/) 999 RETURN END C SUBROUTINE QQTEST(ARRAY,NDATA,IPRINT,QUANT,PROB,IERR) C-------------------------------------------------------------------- C COMPUTES Q-Q PLOTS C EXTERNALS: MDNRIS IS A SUBROUTINE FROM IMSL C INPUT : ARRAY (NDATA) C : NDATA NUMBER OF DATA C : IPRINT = LOGICAL UNIT OF THE PRINTER C OUTPUT : QUANT = ORDERED GAUSSIAN QUANTILES C : PROB = ORDERED PROBABILITY VALUES C : IERR = ERROR CODE C-------------------------------------------------------------------- IMPLICIT REAL*4(A-H,O-Z) DIMENSION ARRAY(NDATA), PROB(NDATA), QUANT(NDATA) C----------------------------------------------------- C COMPUTE CORRESPONDING PROBABILITIES C----------------------------------------------------- DO 1 I=1,NDATA 1 PROB(I) = (I - (1.00/3.00) ) / (NDATA + (1.0/ 3.0) ) C C---------------------------------------------------------- C CALCULATE QUANTILES OF GAUSSIAN DIASTRIBUTION C--------------------------------------------------------- DO 2 I=1,NDATA CALL MDNRIS(PROB(I),QUANT(I),IER) 2 CONTINUE C C--------------------------------------------------------- C PRINT COMPUTED DATA C--------------------------------------------------------- RETURN END //* //LKED.USERLIB DD DSN=A.M12129.SELIBOBJ,DISP=SHR //GO.FT15F001 DD DSN=&&TEMPPVA,DISP=(NEW,PASS), // UNIT=DASD,SPACE=(TRK,(1,1),RLSE), // DCB=(RECFM=FB,LRECL=80,BLKSIZE=8000) //GO.FT22F001 DD DSN=A.M12128.MERTIKAS.FINAL, // UNIT=P3350,DISP=OLD //GO.FT17F001 DD DSN=A.M12128.MERTIKAS.GEDOP, // UNIT=P3350,DISP=OLD //GO.SYSIN DD * 148. 154. 158. 160. 161. 162. 166. 170. 182. 195. 236. //* //STEP2 EXEC SAS,REGION=4096K //PVAFILE DD DSN=&&TEMPPVA,DISP=(OLD,DELETE) //SYSIN DD * GOPTIONS DEVICE=ZETA HSIZE=11.0 VSIZE=8.0 BORDER; * *------------------------------------------------------------ * INPUT DATA *------------------------------------------------------------ *; DATA PVA; INFILE PVAFILE; INPUT Y X; * PROC UNIVARIATE; * VAR Y; *-----------------------------------------------------------; * PLOT OF QUANTILE-QUANTILE; *-----------------------------------------------------------; PROC GPLOT DATA=PVA; PLOT Y*X; LABEL Y=DISTANCE TO PTNS ABOVE MEDIAN X=DISTANCE TO PNTS BELOW MEDIAN; SYMBOL V=X C=BLACK I=NONE; TITLE1 .C=BLACK .F=TITALIC PLOT FOR SYMMETRY; TITLE2 .C=BLACK .H=2 .F=TITALIC X-AXIS (RAW DATA,N=600); //