//RESTON JOB /*JOBPARM S=7,R=2048,L=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*8(A-H,O-Z) INTEGER STEP CHARACTER LABEL*5(20) , OPTION*2(3) C REAL * 4 XBOX(2000) INTEGER KBOX, NIBOX(50), MAXL, IERBOX C 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 TABLE(2000,4) C DIMENSION ARRAY(2000), Y(2000) DIMENSION CLEAN(2000), VECT(2000) DIMENSION YSORT(2000), ABSDIF(2000) DIMENSION SD(200), ST(200), SAD(200), SMAD(200), SF(200) DIMENSION SG(200), SW(200), SMW(200) DIMENSION SB7(200) , SB8(200), SB9(200), SB10(200) DIMENSION SH14(200), SH17(200) , SH20(200) DIMENSION SHL12(200), SHL15(200), SHL22(200), SHL25(200) DIMENSION SGN(200) C DIMENSION SCALE(20,200) , VARLOG(20), EFF(20) DIMENSION WORK(200), VEC(2000) C OPTION(1) = ' A' OPTION(2) = ' B' OPTION(3) = ' C' C LABEL(1) = 'SD ' LABEL(2) = 'ST ' LABEL(3) = 'SAD ' LABEL(4) = 'SMAD ' LABEL(5) = 'SF ' LABEL(6) = 'SG ' LABEL(7) = 'SH14 ' LABEL(8) = 'SH17 ' LABEL(9) = 'SH20 ' LABEL(10)= 'SB7 ' LABEL(11)= 'SB8 ' LABEL(12)= 'SB9 ' LABEL(13)= 'SB10 ' LABEL(14)= 'SW ' LABEL(15)= 'SMW ' LABEL(16)= 'SHL12' LABEL(17)= 'SHL15' LABEL(18)= 'SHL22' LABEL(19)= 'SHL25' LABEL(20)= 'SGN ' DATA PI /3.141592654D0/ C------------------------------------------------------------------- C IOPT = 1 SITUATION A (RAW DATA ) C IOPT = 2 SITUATION B (OUTLIERS HAVE BEEN REMOVED) C IOPT = 3 SITUATION C (STABLE DATA , N =600) C-------------------------------------------------------------------- IOPT = 2 ISAMPL = 1000 RANGE = 1117.0D0 IGEN = 200 ISTART = 10 ISCAL = 20 C MAXL = 129 C IPRINT = 6 ITEMP =15 IDISK =22 IREAD = 5 ALFAT = 0.20D0 ALFAS = 0.10D0 C READ DATA C--------------------------------------------------- OPEN(IDISK,ERR=3000,STATUS='OLD',ACCESS='SEQUENTIAL', # IOSTAT=IERROR) C 9 DO 1 K=1,2000 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 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 R =DSQRT( TABLE(I,1)**2 +TABLE(I,2)**2+TABLE(I,3)**2) TABLE(I,4) = ( R ) 3 CONTINUE L = 1 CALL SET(ARRAY, NDATA,1,NDATA, TABLE(1,L), NDATA) C IF(IOPT .EQ. 2) THEN CALL YINFO8(ARRAY, NDATA, XMED, HL, HH, ADJL,ADJH, IADJL,IADJH, # XSTEP,IOUTL,IOUTH,IERR) C--------------------------------------------------------------------- C REMOVE FAR OUT VALUES FROM ARRAY AND CREATE CLEAN VECTOR C--------------------------------------------------------------------- MAXOUT = IADJH MINOUT = IADJL XMAX = ARRAY(MAXOUT) XMIN = ARRAY(MINOUT) CALL CLEAR(TABLE(1,L),NDATA, XMAX, XMIN, CLEAN,ICLEAN) NDATA = ICLEAN IF(L .EQ .1) THEN RANGE = 1060.D0 IGEN = 200 ENDIF IF(L .EQ. 2) THEN RANGE = NDATA IGEN = 200 ENDIF IF(L .EQ. 3) THEN RANGE = NDATA IGEN = 200 ENDIF CALL SET(ARRAY, NDATA,1,NDATA, CLEAN, NDATA) ENDIF IF(IOPT .EQ. 3 ) THEN IGEN = 200 NDATA = 600 RANGE = 600.0D0 ENDIF NOBS = NDATA C DO 9999 K= 1, IGEN CALL GENRDM(ARRAY, NDATA, K, Y, ISAMPL, RANGE) NDATA = ISAMPL CALL SET(YSORT, ISAMPL,1,ISAMPL, Y, ISAMPL) C------------------------------------------------------------------ C COMPUTE STANDARD DEVIATION C------------------------------------------------------------------ CALL STRDEV(YSORT, NDATA, SD(K)) C----------------------------------------------------------------- C TRIMMED STANDARD DEVIATION C----------------------------------------------------------------- CALL SORT8(YSORT,NDATA) CALL TRIM2A(YSORT, NDATA, ALFAT, XTRIM) CALL TRIMSD(YSORT, NDATA, ALFAS, XTRIM, ST(K)) C---------------------------------------------------------------------- C MEAN ABSOLUTE DEVIATION FROM MEDIAN C---------------------------------------------------------------------- CALL MEANAD(YSORT, NDATA, SAD(K)) C---------------------------------------------------------------------- C MEDIAN ABSOLUTE DEVIATION (MAD) C---------------------------------------------------------------------- CALL XMEDIN(YSORT,NDATA,XMEDN) CALL MAD(YSORT,NDATA,XMEDN,ABSDIF,SMAD(K)) C-------------------------------------------------------------------- C FOURTH-SPREAD C---------------------------------------------------------------------- CALL ESIGMA(YSORT, NDATA, SF(K)) C------------------------------------------------------------------- C GAUSSIAN SKIP C------------------------------------------------------------------- C = 6.0D0 CALL TBIW(YSORT, NDATA, XMEDN, C, SMAD(K), XBIW) CALL GSKIP(YSORT, NDATA, XBIW, 7.D0, SG(K)) C----------------------------------------------------------------- C HUBER'S M-ESTIMATE OF SCALE C----------------------------------------------------------------- B = 1.4D0 CALL HUBER(YSORT, NDATA, XMEDN, SMAD(K), B, SH14(K)) B = 1.7D0 CALL HUBER(YSORT, NDATA, XMEDN, SMAD(K), B, SH17(K)) B = 2.0D0 CALL HUBER(YSORT, NDATA, XMEDN, SMAD(K), B, SH20(K)) C----------------------------------------------------------------- C BIWEIGHT A-ESTIMATE C------------------------------------------------------------------ CONST = 7.0D0 MODE = -1 CALL AESTIM(MODE,YSORT,NDATA,XMEDN,CONST,SMAD(K),SB7(K)) CONST = 8.0D0 MODE = -1 CALL AESTIM(MODE,YSORT,NDATA,XMEDN,CONST,SMAD(K),SB8(K)) CONST = 9.0D0 MODE = -1 CALL AESTIM(MODE,YSORT,NDATA,XMEDN,CONST,SMAD(K),SB9(K)) CONST =10.0D0 MODE = -1 CALL AESTIM(MODE,YSORT,NDATA,XMEDN,CONST,SMAD(K),SB10(K)) C----------------------------------------------------------------- C ANDREWS' SINE A-ESTIMATOR (WAVE) C----------------------------------------------------------------- MODE = 1 CWAVE= PI * 2.1D0 CALL AESTIM(MODE,YSORT,NDATA,XMEDN,CWAVE,SMAD(K),SW(K)) C WRITE(ITEMP,*) SB7(K) , SD(K), K C WRITE(6 ,*) SB7(K) , SD(K), K C----------------------------------------------------------------- C MODIFIED ANDREWS' SINE A-ESTIMATOR (M-WAVE) C----------------------------------------------------------------- MODE = 0 CALL AESTIM(MODE,YSORT,NDATA,XMEDN,CWAVE,SMAD(K),SMW(K)) C------------------------------------------------------------------- C HAMPEL'S A-ESTIMATE(THREE-PART REDESCENDING ESTIMATOR) C------------------------------------------------------------------- A = 1.71D0 B = 2.85D0 C = 11.4D0 CALL HAMPEL(YSORT, NDATA, A, B, C, XMEDN, SMAD(K), SHL12(K)) A = 1.71D0 B = 2.28D0 C = 5.70D0 CALL HAMPEL(YSORT, NDATA, A, B, C, XMEDN, SMAD(K), SHL15(K)) A = 2.25D0 B = 3.75D0 C =15.00D0 CALL HAMPEL(YSORT, NDATA, A, B, C, XMEDN, SMAD(K), SHL22(K)) A = 2.25D0 B = 3.00D0 C = 7.50D0 CALL HAMPEL(YSORT, NDATA, A, B, C, XMEDN, SMAD(K), SHL25(K)) C-------------------------------------------------------------------- C GINI'S ROBUST SCALE ESTIMATOR C-------------------------------------------------------------------- CALL GINI(YSORT, NDATA, SGN(K)) NDATA = NOBS 9999 CONTINUE C C-------------------------------------------------------------------- C COMPUTE RELATIVE EFFICIENCIES OF VAR{ LOG(S) } C--------------------------------------------------------------------- WRITE(6,67) NDATA, IGEN, L, IOPT, ISAMPL 67 FORMAT(2X,'NDATA=',I7,2X,'IGEN=',I4,2X,'L=',I4,2X,'IOPT=',I3,2X, #'ISAMPLE=',I5,//) CALL GENMTX(SD, ST, SAD, SMAD, SF, SG, # SH14, SH17, SH20, # SB7, SB8, SB9, SB10, SW, SMW, # SHL12, SHL15, SHL22, SHL25, # SGN, SCALE , ISCAL, IGEN) DO 60 J=1,IGEN C WRITE(6,61) J,(SCALE(I,J),I=1,ISCAL) 60 CONTINUE 61 FORMAT(1X,I3,2X,12(F9.3,1X),/) CALL VRNLOG(SCALE, ISCAL, IGEN, VARLOG, EFF) C--------------------------------------------------------------------- C BOXPLOT DISPLAYS C--------------------------------------------------------------------- CALL GENBOX(SCALE, ISCAL, IGEN, XBOX, KBOX, NIBOX) CALL USBOX(XBOX,ISCAL, NIBOX, MAXL, IERBOX) C--------------------------------------------------------------------- C LETTER-VALUE DISPLAYS C--------------------------------------------------------------------- WRITE(6,16) OPTION(IOPT), ISAMPL WRITE(6,111) DO 64 K=1, ISCAL DO 65 I=1,IGEN WORK(I) = SCALE(K,I) 65 CONTINUE CALL YINFO8(WORK,IGEN, XMED, HL, HH, ADJL,ADJH, IADJL,IADJH, # XSTEP,IOUTL,IOUTH,IERR) WRITE(6,15) K,LABEL(K),WORK(1),ADJL,HL,XMED,HH,ADJH,WORK(IGEN) 64 CONTINUE 15 FORMAT(2X,I2,2X,A5,2X,7(F10.3),//) 16 FORMAT('1','NUMBER',1X,'SCALE',4X,'MIN',6X,'ADJ-L',6X,'FOURTH-L', # 1X,'MEDIAN',4X,'FOURTH-H',2X,'ADJ-H',6X,'MAX',2X,'SITUATION=', # A2,2X,'SAMPLE SIZE=',I5,/) C-------------------------------------------------------------------- C WRITE FOR PLOTING C------------------------------------------------------------------- 90 WRITE(6, 91) DO 87 I=1,ISCAL WRITE(6,89) I, LABEL(I), VARLOG(I), EFF(I) 87 CONTINUE 89 FORMAT(12X,I2,2X,A5,1X,F7.5,2X,F6.2,2X,//) 91 FORMAT('1',7X,'NUMBER',2X,'SCALE',2X,'VAR{LN(S)}',1X,'EFF',//) C DO 40 I=1,ISCAL C WRITE(ITEMP,*) I, EFF(I) 40 CONTINUE 50 FORMAT(1X, 9(F5.1,1X), I2) 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--------------------------------------------------------------------- 77 FORMAT('1',107X,'VAR(LOG(S))',2X,'EFFICIENCY') 111 FORMAT(2X,'-------------------------------------------------------- #------------------------------------------------------------') 200 FORMAT(2X,'SG=',F10.3,2X,'SH=',F10.3,2X,'K=',I2) 139 FORMAT(2X,I8,4(F20.2,2X)) 142 FORMAT(I8,1X,3(F10.3,2X),F20.4,F20.4,1X,3(F14.2)) 300 FORMAT(2X,'SCALE #',I3,2X,'VARLOG=',F15.8,2X,'EFF=',F10.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 SUBROUTINE GINI(ARRAY, N, SGINI) IMPLICIT REAL * 8 (A-H, O-Z) DIMENSION ARRAY( N ) C-------------------------------------------------------------- C COMPUTES GINI'S ROBUST SCALE ESTIMATE C ARRAY = MUST BE ORDERED C AUTHOR : S. MERTIKAS C DATE : OCTOBER,1986 C-------------------------------------------------------------- EN = DFLOAT( N ) DATA PI /3.141592654D0/ SUM = 0.D0 DO 10 I =1, N VALUE = DFLOAT(I) - (EN + 1.D0) / 2.D0 SUM = SUM + VALUE * ARRAY(I) 10 CONTINUE SGINI = 2.D0 * DSQRT( PI ) * SUM SGINI = SGINI / ( EN * (EN - 1.D0) ) RETURN END C FUNCTION HAM(U, A, B, C) IMPLICIT REAL * 8 (A-H, O-Z) C----------------------------------------------------------------- C COMPUTES HAMPEL'S PSI-FUNCTION C----------------------------------------------------------------- Z = DABS( U ) IF(Z .LT. A) THEN HAM = U RETURN ENDIF IF(Z .GE. A .AND. Z .LT. B) THEN HAM = A RETURN ENDIF IF(Z .GE. B .AND. Z .LT. C) THEN SIGN= DSIGN( 1.D0, U ) HAM = SIGN * A * (C - Z) / (C - B) RETURN ENDIF HAM = 0.D0 RETURN END C FUNCTION DRVHAM(U, A, B, C) IMPLICIT REAL* 8(A-H, O-Z) C-------------------------------------------------------------- C COMPUTES PSI-DERIVATIVE FOR HAMPEL'S PSI-FUNCTION C (THREE-PART REDESCENDING FUNCTION) C------------------------------------------------------------- Z = DABS( U) IF(Z .LT. A) THEN DRVHAM = 1.D0 RETURN ENDIF IF(Z .GE. B .AND. Z .LT. C) THEN DRVHAM = - A / (C - B) RETURN ENDIF DRVHAM = 0.D0 RETURN END C SUBROUTINE GENRDM(ARRAY, NDATA, IGEN, VEC, ISAMPL, RANGE) IMPLICIT REAL * 8(A-H, O-Z) DIMENSION ARRAY( NDATA ), VEC(ISAMPL) C----------------------------------------------------------------- C GENERATES A VECTOR OF SIZE ISAMPL OF RANDOMLY C SELECTED VALUES FROM ARRAY(NDATA) C---------------------------------------------------------------- IY = IGEN DO 10 I=1, ISAMPL RANDOM = URAND(IY) IR = RANDOM * RANGE VEC(I) = ARRAY( IR ) C WRITE(6,*) IR, ,VEC(I), I 10 CONTINUE RETURN END C SUBROUTINE HAMPEL(ARRAY, N, A, B, C, XMEDN, SMAD, HMPL) IMPLICIT REAL * 8 (A-H, O-Z) DIMENSION ARRAY(N) C------------------------------------------------------------------- C COMPUTES HAMPEL'S A-ESTIMATE OF SCALE C AUTHOR : S. MERTIKAS C DATE : OCTOBER,1986 C------------------------------------------------------------------- SUM1 = 0.D0 SUM2 = 0.D0 EN = DFLOAT(N) DO 10 I = 1, N U = (ARRAY(I) - XMEDN) / SMAD SUM1 = SUM1 + (HAM(U, A, B, C)**2) SUM2 = SUM2 + DRVHAM(U, A, B, C) 10 CONTINUE SCALE = EN * ( SMAD ** 2) * SUM1 HMPL = SCALE / ( SUM2 ** 2) HMPL = DSQRT ( HMPL ) RETURN END C SUBROUTINE GENBOX(SCALE, ISCAL, IGEN, XBOX, NBOX, NIBOX) IMPLICIT REAL*8(A-H, O-Z) DIMENSION SCALE(ISCAL, IGEN) REAL*4 XBOX(2000) INTEGER NBOX, NIBOX(ISCAL) C DO 10 I=1, ISCAL K = I - 1 NIBOX(I) = IGEN WRITE(6,*) NIBOX(I),I DO 20 J=1, IGEN JJ = IGEN * K + J XBOX(JJ) = SNGL( SCALE(I,J) ) 20 CONTINUE 10 CONTINUE NBOX = JJ RETURN END C SUBROUTINE CLEAR(ARRAY, N, XMAX, XMIN, CLEAN,J) IMPLICIT REAL *8(A-H,O-Z) DIMENSION ARRAY(N), CLEAN(N) 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----------------------------------------------------------------- J = 0 DO 1 I = 1, N IF( ARRAY(I). GT. XMIN .AND. ARRAY(I) .LT. XMAX) THEN J = J + 1 CLEAN( J ) = ARRAY(I) C WRITE(6,*) J, CLEAN(J), XMAX, XMIN ENDIF 1 CONTINUE RETURN END C SUBROUTINE YINFO8(Y, N, XMED, HL, HH, ADJL, ADJH, IADJL, IADJH, # XSTEP,IOUTL,IOUTH,IERR) C----------------------------------------------------------------- C GET GENERAL INFORMATION ABOUT Y(). USEFUL FOR PLOT SCALING C SORTS Y() AND RETURNS IT SORTED. IT ALSO RETURNS: C XMED = 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---------------------------------------------------------------- IMPLICIT REAL *8 (A-H,O-Z) DIMENSION Y(N) C--------------------------------- C LOCAL VARIABLES C--------------------------------- REAL HFENCE, LFENCE, OUTFL, OUTFH INTEGER J, K, TEMP1, TEMP2 C CALL SORT8(Y, N, IERR) IF(IERR .NE. 0) GO TO 999 K = N J = (K/2) + 1 C TEMP1 = N + 1 - J XMED = (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 XSTEP = (HH - HL) * 1.5 HFENCE = HH + XSTEP LFENCE = HL - XSTEP OUTFH = HH + 2 * XSTEP OUTFL = HL - 2 * XSTEP 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 SUBROUTINE VRNLOG(SCALE, ISCAL, IGEN, VARLOG, EFF) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION SCALE(ISCAL, IGEN), VARLOG(ISCAL), SLOG(200) DIMENSION YY(200), EFF(ISCAL) C------------------------------------------------------------------- C COMPUTES THE VARIANCE OF THE LOGARITHM OF SCALE ESTIMATES C AND CORRESPONDING EFFICIENCIES C AUTHOR : S. MERTIKAS C DATE : AUGUST,1986 C------------------------------------------------------------------- DO 10 I=1, ISCAL DO 20 J=1, IGEN SLOG(J) = DLOG( SCALE(I,J) ) 20 CONTINUE CALL STRDEV( SLOG, IGEN, SIGMA ) VARLOG(I) = SIGMA ** 2 10 CONTINUE CALL SET(YY, ISCAL,1,ISCAL, VARLOG, ISCAL) CALL SORT8(YY, ISCAL) VMIN = YY(1) DO 30 I=1, ISCAL EFF(I) = (VMIN*100.D0) / VARLOG(I) 30 CONTINUE RETURN END C SUBROUTINE GENMTX(SD, ST, SAD, SMAD, SF, SG, # SH14, SH17, SH20, # SB7, SB8, SB9, SB10, SW, SMW, # SHL12, SHL15, SHL22, SHL25, # SGN, SCALE , ISCAL, IGEN) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION SD(IGEN), ST(IGEN), SAD(IGEN), SMAD(IGEN), SF(IGEN) DIMENSION SG(IGEN), SW(IGEN), SMW(IGEN) DIMENSION SB7(IGEN), SB8(IGEN), SB9(IGEN), SB10(IGEN) DIMENSION SH14(IGEN), SH17(IGEN), SH20(IGEN) DIMENSION SHL12(IGEN), SHL15(IGEN), SHL22(IGEN), SHL25(IGEN) DIMENSION SGN(IGEN) DIMENSION SCALE(ISCAL, IGEN) C--------------------------------------------------------------------- C GENERATES MATRIX OF ALL SCALE ESTIMATES FOR COMPUTING EFFICIENCY C--------------------------------------------------------------------- DO 10 K=1, IGEN SCALE(1,K) = SD(K) SCALE(2,K) = ST(K) SCALE(3,K) = SAD(K) SCALE(4,K) = SMAD(K)/ 0.6745D0 SCALE(5,K) = SF(K) / (2. * 0.6745D0) SCALE(6,K) = SG(K) SCALE(7,K) = SH14(K) SCALE(8,K) = SH17(K) SCALE(9,K) = SH20(K) SCALE(10,K)= SB7(K) SCALE(11,K)= SB8(K) SCALE(12,K)= SB9(K) SCALE(13,K)= SB10(K) SCALE(14,K)= SW(K) SCALE(15,K)= SMW(K) SCALE(16,K)= SHL12(K) SCALE(17,K)= SHL15(K) SCALE(18,K)= SHL22(K) SCALE(19,K)= SHL25(K) SCALE(20,K)= SGN(K) 10 CONTINUE RETURN END C SUBROUTINE GENERT(Y, NDIM, ISTART, STEP, ARRAY, NDATA) IMPLICIT REAL * 8(A-H,O-Z) INTEGER STEP DIMENSION Y(NDIM), ARRAY(NDATA) C--------GENERATES SAMPLES FOR SCALE ESTIMATION------------------ DO 10 I=1, NDIM J = ISTART + (I - 1) * STEP IF(J .GT. NDATA) WRITE(6,1000) Y(I) = ARRAY(J) C WRITE(6,*) I, ARRAY(J), J 10 CONTINUE 1000 FORMAT(2X,'ERROR ENCOUNTERED IN SUBROUTINE GENERT (J>NDATA)') RETURN END C FUNCTION HPSI( U, B ) IMPLICIT REAL * 8(A-H,O-Z) C---------------- COMPUTES HUBER'S PSI FUNCTION------------ IF( DABS(U) .LE. B) THEN HPSI = U GO TO 999 ENDIF IF( U .GT. B ) THEN HPSI = B GO TO 999 ENDIF HPSI = - B 999 RETURN END C SUBROUTINE HUBER(ARRAY, NDATA, XMEDN, SMAD, B, SHUBER) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION ARRAY(NDATA) C----------------------------------------------------------------- C COMPUTES HUBER'S M-ESTIMATE OF SCALE C INPUT : ARRAY = VECTOR OF SAMPLE VALUES C NDATA = NUMBER OF DATA C XMEDN = MEDIAN OF SAMPLE DATA C SMAD = MEDIAN ABSOLUTE DEVIATION OF SAMPLE DATA C B = HUBER'S CONSTANT C OUTPUT: SHUBER= HUBER'S M-ESTIMATE OF SCALE C C AUTHOR : S. MERTIKAS C DATE : AUGUST,1986 C---------------------------------------------------------------- ITER = 40 EN = DFLOAT( NDATA ) T = XMEDN S = SMAD DO 10 K = 1,ITER SUMA = 0.D0 SUMB = 0.D0 DO 20 I=1, NDATA U = (ARRAY(I) - T ) / S PSI = HPSI(U, B) SUMA = SUMA + PSI ** 2 IF( DABS(U) .LE. B ) SUMB = SUMB + U * PSI 20 CONTINUE AN = SUMA - (EN -1.D0) BD = SUMB * 2.D0 SOLD = S S = SOLD * (1.D0 + (AN / BD) ) DIFF = DABS(S - SOLD) IF( DIFF .LT. 0.5D0) GO TO 40 10 CONTINUE 40 IF( DIFF .GT. 0.5D0) WRITE(6,1000) ITER SHUBER = S IF(SHUBER .LT. 0D0) WRITE(6,2000) SHUBER, B 1000 FORMAT('1',2X,'SUBROUTINE HUBER COULD NOT CONVERGE AFTER',2X,I3, #'ITERATIONS') 2000 FORMAT(2X,'HUBER SCALE CONVERGES TO NEGATIVE NUMBER',F10.2,2X, #'WITH B=',F8.2) RETURN END C SUBROUTINE GSKIP(ARRAY, NDATA, XBIW, C, SKIP) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION ARRAY(NDATA), U(2000), ABSDIF(2000) C----------------------------------------------------------------- C COMPUTES THE GAUSSIAN SKIP AS A SCALE ESTIMATE C INPUT : ARRAY = VECTOR OF SAMPLE VALUES C NDATA = NUMBER OF DATA C XBIW = BIWEIGHT LOCATION ESTIMATE C C = TUNING CONSTANT C OUTPUT: SKIP = GAUSSIAN SKIP SCALE ESIMATE C C AUTHOR : S. MERTIKAS C DATE : AUGUST,1986 C---------------------------------------------------------------- ALFA = 1.D0 / 3.D0 DO 10 I=1, NDATA PROB = (I - ALFA) / (DFLOAT(NDATA)+ ALFA) CALL MDNRIS( PROB, Z, IER ) V = ( ARRAY(I) - XBIW) / Z U(I) = V ** 2 10 CONTINUE C CALL SORT8( U, NDATA) CALL XMEDIN(U,NDATA, XMEDN) CALL MAD(U,NDATA,XMEDN,ABSDIF,SMAD) CALL TBIW(U, NDATA, XMEDN, C, SMAD, XBIW) C SKIP = DSQRT( XBIW ) RETURN END C SUBROUTINE TBIW(ARRAY, NDATA, XMED, C, SMAD, XBIW) IMPLICIT REAL *8 (A-H,O-Z) DIMENSION ARRAY(NDATA) C----------------------------------------------------------------- C COMPUTES THE BIWEIGHT ESTIMATE OF LOCATION C INPUT : ARRAY = VECTOR OF SAMPLE VALUES C NDATA = NUMBER OF DATA C XMED = MEDIAN OF DATA C C = TUNING CONSTANT (C=6,7,8,9) C SMAD = MEDIAN ABSOLUTE DEVIATION OF DATA C OUTPUT: XBIW = BIWEIGHT LOCATION ESTIMATE C C AUTHOR : S. MERTIKAS C DATE : AUGUST,1986 C---------------------------------------------------------------- SUMW = 0.0D0 WXI = 0.0D0 DO 10 I= 1, NDATA U = ( ARRAY(I) - XMED ) / ( C * SMAD ) IF( DABS(U) .GE. 1.0D0 ) GO TO 10 W = (1.D0 - U ** 2 )**2 WXI = WXI + W * ARRAY(I) SUMW = SUMW + W 10 CONTINUE XBIW = WXI / SUMW RETURN END C SUBROUTINE ESIGMA(ARRAY, NDATA, DF) IMPLICIT REAL *8 (A-H, O-Z) DIMENSION ARRAY( NDATA ) C------------------------------------------------------------ C COMPUTES FOURTH-SPREAD (E-PSEUDOSIGMA) C INPUT : ARRAY = VECTOR OF ORDERED SAMPLE VALUES C NDATA = NUMBER OF DATA C OUTPUT: DF = FOURTH-SPREAD C C AUTHOR: S. MERTIKAS C DATE : AUGUST, 1986 C----------------------------------------------------------- NMED = (DFLOAT(NDATA) + 1.0D0 ) / 2.0 ILOW = (DFLOAT(NMED) + 1.0D0 )/ 2.0 IHI = NDATA - ILOW C DF = ARRAY(IHI) - ARRAY(ILOW) RETURN END C SUBROUTINE MEANAD(ARRAY, NDATA, XMENAD) IMPLICIT REAL * 8(A-H, O-Z) DIMENSION ARRAY(NDATA) C------------------------------------------------------------ C COMPUTES THE MEAN ABSOLUTE DEVIATION FROM MEDIAN C INPUT : ARRAY = VECTOR OF ORDERED SAMPLE VALUES C NDATA = VNUMBER OF DATA C OUTPUT: XMENAD= NEAN ABSOLUTE DEVIATION FROM THE MEDIAN C C AUTHOR: S. MERTIKAS C DATE : AUGUST, 1986 C----------------------------------------------------------- CALL XMEDIN(ARRAY,NDATA,XMEDN) SUM = 0.0D0 DO 10 I = 1, NDATA SUM = SUM + DABS( ARRAY(I) - XMEDN ) 10 CONTINUE XMENAD = SUM / DFLOAT( NDATA ) RETURN END C SUBROUTINE TRIMSD(ARRAY, NDATA, ALFA, XTRIM, STRIM) IMPLICIT REAL * 8 (A-H, O-Z) INTEGER G DIMENSION ARRAY( NDATA ), U(2000) C------------------------------------------------------------------- C COMPUTES THE ONE-SIDED TRIMMED STANDARD DEVIATION C INPUT : ARRAY = VECTOR OF SAMPLE VALUES C NDATA = NUMBER OF DATA C ALFA = TRIMMING PROPORTION C XTRIM = TWO-SIDED TRIMMED MEAN C OUTPUT: STRIM = ONE-SIDED STANDARD DEVIATION C C AUTHOR: S. MERTIKAS C DATE : AUGUST,1986 C----------------------------------------------------------------- EN = DFLOAT(NDATA) DO 10 I=1, NDATA U(I) = (ARRAY(I) - XTRIM)**2 10 CONTINUE CALL SORT8(U, NDATA) G = ALFA * EN + 0.5 C SUM = 0.0D0 MAX = NDATA - G DO 20 I = 1, MAX SUM = SUM + U(I) 20 CONTINUE DEN = DFLOAT(MAX) -1.0D0 S2 = SUM / DEN STRIM = DSQRT( S2 ) RETURN END C SUBROUTINE TRIM2A(ARRAY, NDATA, ALFA, XTRIM) IMPLICIT REAL * 8(A-H, O-Z) INTEGER G DIMENSION ARRAY(NDATA) C----------------------------------------------------------------- C COMPUTES AN ALFA% TWO-SIDED TRIMMED MEAN C INPUT : ARRAY = ORDERED VECTOR OF SAMPLE VALUES C NDATA = NUMBER OF DATA C ALFA = TRIMMING PROPORTION C OUTPUT : XTRIM = TWO-SIDED TRIMMED MEAN C C AUTHOR : S. MERTIKAS C DATE : AUGUST,1986 C---------------------------------------------------------------- EN = DFLOAT( NDATA ) G = ALFA * EN+ 0.5 R = ALFA * EN - G SUM = 0.0D0 MIN = G + 2 MAX = NDATA - G - 1 C DO 10 I= MIN, MAX SUM = SUM + ARRAY(I) 10 CONTINUE SUM = SUM + (1.D0 - R) * ( ARRAY(G+1) + ARRAY(NDATA-G) ) DEN = EN - 2.D0* G - 1.D0 XTRIM = SUM / DEN RETURN END C SUBROUTINE STRDEV(ARRAY, NDATA, SD) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION ARRAY( NDATA ) C----------------------------------------------------------------- C COMPUTES STANDARD DEVIATION OF A DATA SET ARRAY C C INPUT : ARRAY = VECTOR OF SAMPLE DATA C NDATA = SAMPLE SIZE C OUTPUT: SD = STANDARD DEVIATION C C AUTHOR: S. MERTIKAS C DATE : AUGUST,1986 C--------------------------------------------------------------- CALL XMEAN(ARRAY,NDATA,XBAR) SUM = 0.0D0 DO 10 I =1, NDATA SUM = SUM + (ARRAY(I) - XBAR) ** 2 10 CONTINUE SD = SUM / DFLOAT(NDATA - 1) SD = DSQRT( SD ) RETURN END C SUBROUTINE SET(ARRAY, N, IFIRST, ILAST, ARRAY1, NDATA) IMPLICIT REAL*8(A-H,O-Z) 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 AESTIM(MODE,ARRAY,NDATA,XMEDN,C,SMAD,ASCALE) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ARRAY(NDATA) C--------------------------------------------------------------- C COMPUTES A - ESTIMATES OF SCALE C INPUT : MODE = - 1 USES THE BIWEIGHT FUNCTION C = + 1 USES THE ANDREW'S SINE FUNCTION C = 0 USES THE MODIFIED ANDREW'S SINE FUNCTION C ARRAY = VECTOR OF SAMPLE DATA C NDATA = NUMBER OF DATA C XMEDN = MEDIAN OF SAMPLE DATA C C = TUNING CONSTANT C SMAD = MEDIAN ABSOLUTE DEVIATION OF DATA C OUTPUT: ASCALE= A-ESTIMATE OF SCALE C C AUTHOR : S.MERTIKAS C DATE : AUGUST,1986 C---------------------------------------------------------------- DATA PI /3.141592654D0/ EN = DFLOAT( NDATA) SQN = DSQRT( EN ) FACTOR = EN / DSQRT( EN - 1.0D0) IF(MODE) 2, 6, 4 2 SUMA = 0.D0 SUMB = 0.D0 DO 3 I = 1, NDATA U = (ARRAY(I) - XMEDN)/(C * SMAD) IF(DABS(U) .GE. 1.D0) GO TO 3 QA = (ARRAY(I) - XMEDN) ** 2 QB = (1.D0 - U **2) ** 4 PSI = QA * QB SUMA = SUMA + PSI PSID = (1.D0 - U**2 ) * (1.D0 - 5.D0 * (U ** 2)) SUMB = SUMB + PSID 3 CONTINUE ANUM = ( FACTOR ) * DSQRT( SUMA ) BDEN = DABS( SUMB ) ASCALE = ANUM / BDEN RETURN 4 SUMA = 0.D0 SUMB = 0.D0 DO 5 I=1,NDATA U = (ARRAY(I) - XMEDN)/(C * SMAD) IF(DABS(U) .GE. 1.D0) GO TO 5 PSI = DSIN( PI * U ) SUMA = SUMA + PSI ** 2 PSID = DCOS( PI * U ) SUMB = SUMB + PSID 5 CONTINUE ANUM = (C * SMAD ) * FACTOR * DSQRT( SUMA ) BDEN = PI * DABS(SUMB ) ASCALE = ANUM / BDEN RETURN 6 SUMA = 0.D0 SUMB = 0.D0 DO 7 I=1,NDATA U = (ARRAY(I) - XMEDN)/(C * SMAD) IF(DABS(U) .GE. 1.D0) GO TO 7 PSI = DSIN( PI * U ) SUMA = SUMA + PSI ** 2 PSID = DCOS( PI * U ) SUMB = SUMB + PSID 7 CONTINUE ANUM = DSQRT( SUMA ) BDEN = PI * DABS(SUMB ) QA = DATAN( ANUM / BDEN) ASCALE = (C * SMAD ) * FACTOR * QA RETURN END C SUBROUTINE XMEDIN(ARRAY,NDIM,XMEDN) IMPLICIT REAL*8(A-H,O-Z) DIMENSION ARRAY(NDIM) M = (NDIM +1) / 2 N = (NDIM +2) / 2 XMEDN = .5D0*(ARRAY(M) + ARRAY(N)) RETURN END C SUBROUTINE MAD(ARRAY,NDIM,XMEDN,ABSDIF,XMAD) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ARRAY(NDIM), ABSDIF(NDIM) C------------------------------------------------------------------ C COMPUTES THE MEDIAN ABSOLUTE DEVIATION C INPUT : ARRAY = VECTOR OF ORDERED SAMPLE VALUES C NDIM = NUMBER OF DATA C XMEDN = MEDIAN OF THE SAMPLE (I.E.,ARRAY) C OUTPUT: XMAD = MEDIAN ABSOLUTE DEVIATION C C AUTHOR :S. MERTIKAS C DATE :JUNE,1985 C------------------------------------------------------------------ DO 1 I=1,NDIM ABSDIF(I) = DABS(ARRAY(I) - XMEDN) 1 CONTINUE CALL SORT8(ABSDIF,NDIM) CALL XMEDIN(ABSDIF,NDIM,XMAD) RETURN END C SUBROUTINE XMEAN(ARRAY,NDATA,XBAR) IMPLICIT REAL*8(A-H,O-Z) DIMENSION ARRAY(NDATA) XBAR = 0.D0 DO 1 I=1,NDATA 1 XBAR=XBAR+ARRAY(I) XBAR=XBAR/DFLOAT(NDATA) RETURN END C SUBROUTINE SORT8(V1,N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V1(N) I = 1 1 I = I + I IF(I.LE.N) GO TO 1 M = I - 1 2 M = M / 2 IF(M.EQ.0) RETURN K = N - M DO 4 J=1,K L = J 5 IF(L .LT. 1) GO TO 4 IF(V1(L+M) .GE. V1(L)) GO TO 4 X = V1(L+M) V1(L+M) = V1(L) V1(L) = X L = L - M GO TO 5 4 CONTINUE GO TO 2 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 //* //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 Y1 Y2 X; *-----------------------------------------------------------; * PLOT OF QUANTILE-QUANTILE; *-----------------------------------------------------------; PROC GPLOT DATA=PVA; PLOT Y1*X Y2*X /OVERLAY HAXIS=0 TO 200 BY 20 VAXIS=0 TO 15 BY 1; LABEL Y1=SCALE ESTIMATES X=NUMBER OF SAMPLE OF SIZE 100; SYMBOL1 V=SQUARE C=BLACK I=NONE L=1; SYMBOL2 V=STAR C=BLACK I=JOIN L=2; TITLE1 .C=BLACK .F=TITALIC STANDARD DEV & BIWEIGHT VS SAMPLE NUMBER; TITLE2 .C=BLACK .H=2 .F=TITALIC X-AXIS (N=100, SITUATION:B); //