//NORMTEST JOB /*JOBPARM S=8,R=2048,L=4 //STEP1 EXEC FORTVCLG,RC=1024K,RL=1024K,RG=3072K, // PARM.FORT='LANGLVL(77),NOMAP,NOXREF,OPTIMIZE(1)', // PARM.LKED='LET' //FORT.SYSIN DD * IMPLICIT REAL*8(A-H,O-Z) DIMENSION ARRAY(1000), ABSDIF(1000), PSIF(1000), PSIDRV(1000), # U(1000) DIMENSION DX(2000), DY(2000) , DZ(2000) DIMENSION DT(2000), DL1(2000) ,INDEX(2000) DIMENSION DL2(2000), DL3(2000) ,DL4(2000) IPRINT = 6 IREAD = 5 IDISK = 22 BETA = 0.95 MODE = -1 CONST= 9.0D0 C C-------------------------------------------------- C READ TEST DATA C--------------------------------------------------- OPEN(IDISK,ERR=3000,STATUS='OLD',ACCESS='SEQUENTIAL', # IOSTAT=IERROR) C--------------------------------------------------- C READ TEST DATA 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 WRITE DATA C-------------------------------------------------- DO 3 M= 1,NDATA ARRAY(M) = DL1(M) C WRITE(6,*) VECTOR(M), M, NDATA 3 CONTINUE C----------------------------------------------------------- C PERFORM ROBUST TEST FOR DEPARTURE FROM NORMALITY C----------------------------------------------------------- CALL ROBNRM(ARRAY,NDATA,IPRINT,IREAD,CONST,MODE,BETA, # XMEDN,XMAD,PSIF,NPSI,PSIDRV,U,ABSDIF,STATI,STATIB,PERBET) DO 33 I=1,NDATA 33 WRITE(IPRINT,1000) I,PSIDRV(I) GO TO 2020 3000 WRITE(IPRINT,9040) IERROR GO TO 2020 C-------------------------------------------------- C CLOSE FILE ON DISK=10 C------------------------------------------------- 2010 WRITE(IPRINT,9030) IERROR GO TO 2020 2020 CLOSE(IDISK,ERR=4010,STATUS='KEEP',IOSTAT=IERROR) GO TO 2222 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') 4010 WRITE(IPRINT,9070) IERROR 2222 CONTINUE 142 FORMAT(I8,1X,3(F10.3,2X),F20.2,F20.2,1X,3(F12.2)) 1000 FORMAT(' ',5X,I5,F8.4,/) STOP END C SUBROUTINE ROBNRM(ARRAY,NDATA,IPRINT,IREAD,CONST,MODE,BETA, # XMEDN,XMAD,PSIF,NPSI,PSIDRV,U,ABSDIF,STATI,STATIB,PERBET) IMPLICIT REAL*8(A-H,O-Z) DIMENSION ARRAY(NDATA), PSIF(NDATA), PSIDRV(NDATA), U(NDATA) DIMENSION ABSDIF(NDATA) C---------------------------------------------------- C COMPUTE MEDIAN C----------------------------------------------------- CALL SORT8(ARRAY,NDATA) CALL XMEDIN(ARRAY,NDATA,XMEDN) C----------------------------------------------------- C COMPUTE MEDIAN ABSOLUTE DEVIATION (MAD) C----------------------------------------------------- CALL MAD(ARRAY,NDATA,XMEDN,ABSDIF,SUMSQR,XMAD) C C----------------------------------------------------- C COMPUTE VECTOR OF PSI-FUNCTION AND ITS DERIVATIVE C----------------------------------------------------- CALL PSIFUN(MODE,ARRAY,NDATA,XMEDN,CONST,XMAD, # PSIF,PSIDRV,U,NPSI) C C----------------------------------------------------- C COMPUTE A-ESTIMATE OF SCALE C----------------------------------------------------- CALL AESTIM(CONST,XMAD,PSIF,PSIDRV,NPSI,NDATA,ASCALE) C C------------------------------------------------------ C COMPUTE TEST STATISTIC I C------------------------------------------------------ CALL TESTI(SUMSQR,NDATA,ASCALE,STATI) C C------------------------------------------------------ C COMPUTE PERCENTILE IB C------------------------------------------------------ CALL TESTIB(NDATA,BETA,STATIB) PERBET = BETA * 100. C------------------------------------------------------ C WRITE COMPUTED DATA C------------------------------------------------------ WRITE(IPRINT,700) WRITE(IPRINT,99) WRITE(IPRINT,700) WRITE(IPRINT,100) XMEDN WRITE(IPRINT,200) XMAD IF(MODE.GT.0) THEN WRITE(IPRINT,300) ASCALE GO TO 3 ENDIF WRITE(IPRINT,400) ASCALE 3 WRITE(IPRINT,500) STATI WRITE(IPRINT,600) PERBET, STATIB IF(STATI.LE.STATIB) WRITE(IPRINT,800) PERBET IF(STATI.GT.STATIB) WRITE(IPRINT,900) PERBET WRITE(IPRINT,700) 99 FORMAT('1',6X,'ROBUST TEST FOR DEPARTURES FROM GAUSSIAN SHAPE') 100 FORMAT(' ',5X,'MEDIAN =',F10.3) 200 FORMAT(' ',5X,'MEDIAN ABSOLUTE DEVIATION =',F10.3) 300 FORMAT(' ',5X,'A-ESTIMATE OF SCALE =',F10.3,2X,'USING ANDREWS #SINE') 400 FORMAT(' ',5X,'A-ESTIMATE OF SCALE =',F10.3,2X,'USING BIWEIGHT') 500 FORMAT(' ',5X,'TEST STATISTIC =',F10.3) 600 FORMAT(' ',5X,'THE',1X,F5.1,'TH',1X,'PERCENTILE =',F10.3,//) 700 FORMAT(' ','------------------------------------------------------- #-------------------') 800 FORMAT(' ',5X,'TEST PASS AT',F5.1,'%',1X,'CONFIDENCE LEVEL') 900 FORMAT(' ',5X,'TEST FAILS AT',F5.1,'%',1X,'CONFIDENCE LEVEL') 1000 FORMAT(' ',5X,I5,F8.2,/) 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 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,SUMSQR,XMAD) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ARRAY(NDIM), ABSDIF(NDIM) SUMS = 0.0D0 DO 1 I=1,NDIM ABSDIF(I) = DABS(ARRAY(I) - XMEDN) 1 SUMS = SUMS + ABSDIF(I)**2 SUMSQR =SUMS CALL SORT8(ABSDIF,NDIM) CALL XMEDIN(ABSDIF,NDIM,XMAD) RETURN END C SUBROUTINE PSIFUN(MODE,ARRAY,NDIM,XMEDN,CONST,XMAD, # PSIF,PSIDRV,U,NU) C C--------------------------------------------------------------- C COMPUTES PSI-FUNCTION(PSIF) AND ITS DERIVATIVE (PSIDRV) C IF MODE----> NEGATIVE OR ZERO: BIWEIGHT FUNCTION (-1) C ----> POSITIVE : ANDREW'S SINE FUNCTION (+1) C AUTHOR : S.MERTIKAS 20 JUNE 1985 C---------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ARRAY(NDIM), PSIF(NDIM), PSIDRV(NDIM), U(NDIM) DATA PI, ICOUNT /3.141592654D0, 0/ DO 1 I=1,NDIM 1 U(I) =(ARRAY(I) - XMEDN)/(CONST * XMAD) IF(MODE) 2, 2, 4 2 DO 3 I=1,NDIM IF(DABS(U(I)) .GE. 1.) GO TO 3 ICOUNT = ICOUNT + 1 PSIF(I) = U(I) * (((1 - U(I)**2))**2) PSIDRV(I) = (1 - U(I)**2) * (1 -5*(U(I)**2)) 3 CONTINUE NU = ICOUNT RETURN 4 DO 5 I=1,NDIM IF(DABS(U(I)) .GE. 1.) GO TO 5 ICOUNT = ICOUNT + 1 PSIF(I) = DSIN( PI*U(I) ) PSIDRV(I) = PI * ( DCOS(PI*U(I)) ) 5 CONTINUE NU = ICOUNT RETURN END C SUBROUTINE AESTIM(CONST,XMAD,PSIF,PSIDRV,NPSI,NDATA,ASCALE) C---------------------------------------------------------- C COMPUTES A-ESTIMATES OF SCALE(VARIANCE) C AUTHOR :S. MERTIKAS 20 JUNE,1985. C---------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION PSIF(NPSI), PSIDRV(NPSI) SSPSI = 0.0D0 SUMDRV= 0.0D0 DO 1 I=1,NPSI SSPSI = SSPSI +(PSIF(I))**2 SUMDRV = SUMDRV + PSIDRV(I) 1 CONTINUE RSS = DSQRT(SSPSI) ABSUM = DABS(SUMDRV) XNUM =CONST * XMAD * (DSQRT(DFLOAT(NDATA))) * RSS ASCALE = XNUM/ABSUM RETURN END C SUBROUTINE TESTI(SUMSQR,NDATA,ASCALE,STATI) C-------------------------------------------------------- C COMPUTES THE OVERALL STATISTIC (STATI) C REFERENCE:CHAP.12 "UNDERSTANDING ROBUST & E.D.A." C AUTHOR: S.MERTIKAS 21 JUNE,1985. C-------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) VAR = ASCALE**2 STATI = SUMSQR/((NDATA - 1)*VAR) RETURN END C SUBROUTINE TESTIB(NDATA,BETA,STATIB) C--------------------------------------------------------- C COMPUTES THE 100B PERCENTAGE POINT (STATIB) C REFERENCE: CHAP.12 "UNDERSTANDING ROBUST & E.D.A." C AUTHOR: S.MERTIKAS 21 JUNE,1985. C--------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) A1 = 0.6376D0 A2 = 1.1535D0 A3 = 0.1266D0 B1 = 1.9065D0 B2 = 2.5465D0 B3 = 0.5652D0 C1 = 0.7824D0 C2 = 1.1021D0 C3 = 0.1021D0 ST = 0.982D0 XNDATA = DLOG10(DFLOAT(NDATA - 1)) IF(BETA.EQ.0.95) GO TO 1 RI90 = A1 - A2*XNDATA + A3*(XNDATA**2) STATIB = ST + 10**RI90 RETURN 1 IF(NDATA.GE.50) GO TO 2 RI95 = B1 - B2*XNDATA + B3*(XNDATA**2) STATIB = ST + 10**RI95 RETURN 2 RI95 = C1 - C2*XNDATA + C3*(XNDATA**2) STATIB = ST + 10**RI95 RETURN END //* //LKED.USERLIB DD DSN=A.M12129.SELIBOBJ,DISP=SHR //GO.FT22F001 DD DSN=A.M12128.MERTIKAS.FINAL, // UNIT=P3350,DISP=OLD //GO.SYSIN DD * //