SUBROUTINE OUTLY(X,N,K,EK,ESD,KUR,RST,INDEX,DET,IFAULT) C C THIS ROUTINE CALCULATES FOUR OUTLIER DETECTION STATISTICS. C C C REFERENCE: SIMONOFF, J.S. (1984) "THE CALCULATION OF OUTLIER C DETECTION STATISTICS," COMMUNICATIONS IN STATISTICS - C ------------------------------ C SIMULATION AND COMPUTATION, 13, 275-285. C -------------------------- C C C THE STRUCTURE OF THE PARAMETERS IS AS FOLLOWS: C X REAL ARRAY (N) INPUT: THE SAMPLE (X1,...,XN) C N INTEGER INPUT: THE SAMPLE SIZE C K INTEGER INPUT: THE MAXIMUM NUMBER OF C OUTLIERS TO BE TESTED C EK REAL ARRAY (K) OUTPUT: THE SAMPLE VALUES OF THE C EK STATISTIC; EK(I) IS THE C VALUE WHEN CHECKING FOR I C OUTLIERS C ESD REAL ARRAY (K) OUTPUT: THE SAMPLE VALUES OF THE C ESD STATISTIC C KUR REAL ARRAY (K) OUTPUT: THE SAMPLE VALUES OF THE C KUR STATISTIC C RST REAL ARRAY (K) OUTPUT: THE SAMPLE VALUES OF THE C RST STATISTIC C INDEX INTEGER ARRAY OUTPUT: INDEX(1,I) IS THE LOCATION C (2,K) IN X OF THE ITH OUTLIER C WHEN USING EK,ESD OR KUR. C INDEX(2,I) IS THE LOCATION C IN X OF THE ITH OUTLIER C WHEN USING RST C DET LOGICAL ARRAY SCRATCH FILE C (2,N) C IFAULT INTEGER OUTPUT: FAULT INDICATOR C 0 = NO FAULT C 1 = MAXIMUM NUMBER OF C OUTLIERS TOO SMALL (K<1) C OR TOO LARGE (K>N-1) ; NO C OUTLIERS TESTED C 2 = SAMPLE X NOT ORDERED ; C RST NOT CALCULATED C 3 = N < 2K+2 ; RST NOT C CALCULATED C 4 = TRIMMED S.D. B=0 ; RST C NOT CALCULATED C C THE FOLLOWING RESTRICTIONS ARE ASSUMED: C (1) THE SAMPLE X IS ORDERED (EITHER ASCENDING OR C DESCENDING ORDER). THIS IS REQUIRED FOR C CALCULATION OF RST ONLY ; THE OTHER STATISTICS C ARE CALCULATED AS NORMAL IF X IS NOT ORDERED. C (2) THE SAMPLE SIZE IS AT LEAST 2K+2. THIS IS REQUIRED C SO THAT THE TRIMMED SAMPLE IS AT LEAST OF SIZE 2, C AND IS REQUIRED FOR CALCULATION OF RST ONLY. C (3) THE MAXIMUM NUMBER OF OUTLIERS TESTED MUST BE AT C LEAST 1 AND AT MOST N-1. IF THIS RESTRICTION IS C VIOLATED, NO STATISTICS WILL BE CALCULATED. C REAL X(N),EK(K),ESD(K),KUR(K),RST(K) INTEGER INDEX(2,K) LOGICAL DET(2,N) DATA XMN,VAR,SK,XKUR,TRMN,TRSD,CHXO/7*0./,IMAXT/1/ IFAULT=0 C C TEST FOR VALID PARAMETERS C IF (K.LT.1) GOTO 85 IF (K.LT.N) GOTO 90 85 IFAULT=1 RETURN 90 IF (N.GT.(2*K+1)) GOTO 95 IFAULT=3 95 XN=FLOAT(N) KP1=K+1 NMK=N-K NMKP1=NMK+1 XFIO=X(1) C C FORM MEAN, SUM OF SQUARES, SUM OF CUBES,SUM OF QUARTICS, C TRIMMED MEAN AND TRIMMED SD C DO 100 I=1,K C C CHECK THAT X IS ORDERED C XFI=X(I) IF (IFAULT.GT.0) GOTO 99 CHX=XFI-XFIO IF ((CHX*CHXO).LT.0.) IFAULT=2 CHXO=CHX XFIO=XFI 99 INDEX(1,I)=0 INDEX(2,I)=0 DET(1,I)=.FALSE. DET(2,I)=.FALSE. XI=FLOAT(I) XIM1=XI-1. DX=(XFI-XMN)/XI DX2=DX*DX XKUR=XKUR-DX*(4.*SK-DX*(6.*VAR+XIM1*(1.+XIM1**3)*DX2)) SK=SK-DX*(3.*VAR-XI*XIM1*(XI-2.)*DX2) VAR=VAR+XI*XIM1*DX2 XMN=XMN+DX 100 CONTINUE IF (IFAULT.NE.3) GOTO 101 NMKP1=K+1 TRSD=1. GOTO 106 101 DO 105 I=KP1,NMK XFI=X(I) IF (IFAULT.GT.0) GOTO 104 CHX=XFI-XFIO IF ((CHX*CHXO).LT.0.) IFAULT=2 CHXO=CHX XFIO=XFI 104 DET(1,I)=.FALSE. DET(2,I)=.FALSE. XI=FLOAT(I) XIM1=XI-1. XIMK=XI-FLOAT(K) XIMKM1=XIMK-1. DX=(XFI-XMN)/XI DX2=DX*DX XKUR=XKUR-DX*(4.*SK-DX*(6.*VAR+XIM1*(1.+XIM1**3)*DX2)) SK=SK-DX*(3.*VAR-XI*XIM1*(XI-2.)*DX2) VAR=VAR+XI*XIM1*DX2 XMN=XMN+DX DX=(XFI-TRMN)/XIMK DX2=DX*DX TRSD=TRSD+XIMK*XIMKM1*DX2 TRMN=TRMN+DX 105 CONTINUE TRSD=SQRT(TRSD/(XN-2.*FLOAT(K)-1.)) 106 DO 110 I=NMKP1,N XFI=X(I) IF (IFAULT.GT.0) GOTO 109 CHX=XFI-XFIO IF ((CHX*CHXO).LT.0.) IFAULT=2 CHXO=CHX XFIO=XFI 109 DET(1,I)=.FALSE. DET(2,I)=.FALSE. XI=FLOAT(I) XIM1=XI-1. DX=(XFI-XMN)/XI DX2=DX*DX XKUR=XKUR-DX*(4.*SK-DX*(6.*VAR+XIM1*(1.+XIM1**3)*DX2)) SK=SK-DX*(3.*VAR-XI*XIM1*(XI-2.)*DX2) VAR=VAR+XI*XIM1*DX2 XMN=XMN+DX 110 CONTINUE IF (VAR.EQ.0.) GOTO 155 IF (IFAULT.EQ.2) GOTO 114 IF (TRSD.EQ.0.) IFAULT=4 C C FIND MOST OUTLYING POINT IN SUBSAMPLE C 114 TVAR=VAR DO 145 KC=1,K IF (VAR.NE.0.) GOTO 120 IKC=KC DO 115 II=IKC,K ESD(II)=0. RST(II)=0. KUR(II)=0. EK(II)=1. 115 CONTINUE RETURN 120 XMAX=0. XMAXT=0. DO 130 I=1,N IF (DET(1,I)) GOTO 125 T=ABS(X(I)-XMN) IF (T.LE.XMAX) GOTO 125 IMAX=I XMAX=T 125 IF (TRSD.EQ.0.) GOTO 130 IF (DET(2,I)) GOTO 130 T=ABS(X(I)-TRMN) IF (T.LE.XMAXT) GOTO 130 IMAXT=I XMAXT=T 130 CONTINUE DET(1,IMAX)=.TRUE. DET(2,IMAXT)=.TRUE. TARG=X(IMAX) XNMKC=XN-FLOAT(KC) XNMKC1=XNMKC+1. C C FORM DETECTION STATISTICS AND UPDATE MEAN,SUM OF SQUARES, C SUM OF CUBES, SUM OF QUARTICS C ESD(KC)=ABS(TARG-XMN)/SQRT(VAR/XNMKC) IF (IFAULT.EQ.0) GOTO 135 RST(KC)=0. GOTO 140 135 RST(KC)=ABS(X(IMAXT)-TRMN)/TRSD INDEX(2,KC)=IMAXT 140 KUR(KC)=(XKUR/(VAR*VAR))*XNMKC1 XMN=(XNMKC1*XMN-TARG)/XNMKC DX=(TARG-XMN)/XNMKC1 DX2=DX*DX VAR=VAR-XNMKC*XNMKC1*DX2 EK(KC)=VAR/TVAR SK=SK+DX*(3.*VAR-XNMKC1*XNMKC*(XNMKC-1.)*DX2) XKUR=XKUR+DX*(4.*SK-DX*(6.*VAR+XNMKC*(1.+XNMKC**3)*DX2)) INDEX(1,KC)=IMAX 145 CONTINUE RETURN 155 DO 160 I=1,K ESD(I)=0. RST(I)=0. KUR(I)=0. EK(I)=1. 160 CONTINUE RETURN END