SUBROUTINE EXTREM(X,N) C C PURPOSE--THIS SUBROUTINE PERFORMS AN EXTREME VALUE ANALYSIS C ON THE DATA IN THE INPUT VECTOR X. C THIS ANALYSIS CONSISTS OF DETERMINING THAT PARTICULAR C EXTREME VALUE TYPE 1 OR EXTREME VALUE TYPE 2 DISTRIBUTION C WHICH BEST FITS THE DATA SET. C THE GOODNESS OF FIT CRITERION IS THE MAXIMUM PROBABILITY C PLOT CORRELATION COEFFICIENT CRITERION. C AFTER THE BEST-FIT DISTRIBUTION IS DETERMINED, C ESTIMATES ARE COMPUTED AND PRINTED OUT FOR THE C LOCATION AND SCALE PARAMETERS. C TWO PROBABILITY PLOTS ARE ALSO PRINTED OUT-- C THE BEST-FIT TYPE 2 PROBABILITY PLOT C (IF THE BEST FIT WAS IN FACT A TYPE 2), C AND THE TYPE 1 PROBABILITY PLOT. C PREDICTED EXTREMES FOR VARIOUS RETURN PERIODS ARE C ALSO COMPUTED AND PRINTED OUT. C INPUT ARGUMENTS--X = THE SINGLE PRECISION VECTOR OF C (UNSORTED OR SORTED) OBSERVATIONS. C N = THE INTEGER NUMBER OF OBSERVATIONS C IN THE VECTOR X. C OUTPUT--6 PAGES OF AUTOMATIC PRINTOUT. C PRINTING--YES. C RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N C FOR THIS SUBROUTINE IS 7500. C OTHER DATAPAC SUBROUTINES NEEDED--SORT, UNIMED, EV1PLT, C EV2PLT, PLOT. C FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, ALOG. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C REFERENCES--FILLIBEN (1972), 'TECHNIQUES FOR TAIL LENGTH C ANALYSIS', PROCEEDINGS OF THE EIGHTEENTH C CONFERENCE ON THE DESIGN OF EXPERIMENTS IN C ARMY RESEARCH AND TESTING, PAGES 425-450. C --FILLIBEN, 'THE PERCENT POINT FUNCTION', C UNPUBLISHED MANUSCRIPT. C --JOHNSON AND KOTZ (1970), CONTINUOUS UNIVARIATE C DISTRIBUTIONS-1, 1970, PAGES 272-295. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE--301-921-2315 C ORIGINAL VERSION--JUNE 1972. C UPDATED --DECEMBER 1974. C UPDATED --NOVEMBER 1975. C UPDATED --MAY 1976. C C--------------------------------------------------------------------- C CHARACTER*4 BLANK,ALPHAM,ALPHAA,ALPHAX CHARACTER*4 ALPHAI,ALPHAN,ALPHAF,ALPHAT,ALPHAY CHARACTER*4 ALPHAG,EQUAL C CHARACTER*4 IFLAG1 CHARACTER*4 IFLAG2 CHARACTER*4 IFLAG3 C DIMENSION W(3000) DIMENSION X(1) DIMENSION Y(7500),Z(7500) DIMENSION GAMTAB(50),CORR(50) DIMENSION YI(50),YS(50),T(50) DIMENSION IFLAG1(50),IFLAG2(50),IFLAG3(50) CCCCC DIMENSION C(10) DIMENSION AM(50) DIMENSION SCRAT(50) C DIMENSION AINDEX(50) CCCCC DIMENSION P0(10) DIMENSION H(60,2) COMMON /BLOCK2/ WS(15000) EQUIVALENCE (Y(1),WS(1)),(Z(1),WS(7501)) DATA BLANK,ALPHAM,ALPHAA,ALPHAX/' ','M','A','X'/ DATA ALPHAI,ALPHAN,ALPHAF,ALPHAT,ALPHAY/'I','N','F','T','Y'/ DATA ALPHAG,EQUAL/'G','='/ DATA GAMTAB(1),GAMTAB(2),GAMTAB(3),GAMTAB(4),GAMTAB(5), 1GAMTAB(6),GAMTAB(7),GAMTAB(8),GAMTAB(9),GAMTAB(10), 1GAMTAB(11),GAMTAB(12),GAMTAB(13),GAMTAB(14),GAMTAB(15), 1GAMTAB(16),GAMTAB(17),GAMTAB(18),GAMTAB(19),GAMTAB(20), 1GAMTAB(21),GAMTAB(22),GAMTAB(23),GAMTAB(24),GAMTAB(25) 1/1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12., 113.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25./ DATA GAMTAB(26),GAMTAB(27),GAMTAB(28),GAMTAB(29),GAMTAB(30), 1GAMTAB(31),GAMTAB(32),GAMTAB(33),GAMTAB(34),GAMTAB(35), 1GAMTAB(36),GAMTAB(37),GAMTAB(38),GAMTAB(39),GAMTAB(40), 1GAMTAB(41),GAMTAB(42) 1/30.,35.,40.,45.,50.,60.,70.,80.,90.,100.,150.,200.,250., 1350.,500.,750.,1000./ CCCCC DATA C(1),C(2),C(3),C(4),C(5),C(6),C(7),C(8),C(9),C(10) CCCCC1/60.,75.,100.,150.,250.,500.,1000.,10000.,100000.,1000000./ CCCCC DATA P0(1),P0(2),P0(3),P0(4),P0(5),P0(6),P0(7),P0(8),P0(9),P0(10) CCCCC1/.0,.5,.75,.9,.95,.975,.99,.999,.9999,.99999/ DATA T(1),T(2),T(3),T(4),T(5),T(6),T(7),T(8),T(9),T(10), 1T(11),T(12),T(13),T(14),T(15),T(16),T(17),T(18),T(19),T(20), 1T(21),T(22),T(23),T(24),T(25) 1/10.18011,3.39672,2.47043,2.14609,1.98712,1.89429,1.83394, 11.79175,1.76069,1.73691,1.71814,1.70297,1.69045,1.67996, 11.67103,1.66335,1.65667,1.65082,1.64564,1.64102,1.63689, 11.63316,1.62979,1.62672,1.62391/ DATA T(26),T(27),T(28),T(29),T(30), 1T(31),T(32),T(33),T(34),T(35),T(36),T(37),T(38),T(39),T(40), 1T(41),T(42),T(43) 1/1.61287,1.60516,1.59947,1.59510,1.59164,1.58651,1.58289, 11.58019,1.57811,1.57645,1.57152,1.56908,1.56763,1.56666, 11.56546,1.56377,1.56330,1.56187/ DATA AINDEX(1),AINDEX(2),AINDEX(3),AINDEX(4),AINDEX(5), 1AINDEX(6),AINDEX(7),AINDEX(8),AINDEX(9),AINDEX(10), 1AINDEX(11),AINDEX(12),AINDEX(13),AINDEX(14),AINDEX(15), 1AINDEX(16),AINDEX(17),AINDEX(18),AINDEX(19),AINDEX(20), 1AINDEX(21),AINDEX(22),AINDEX(23),AINDEX(24),AINDEX(25) 1/1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12., 113.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25./ DATA AINDEX(26),AINDEX(27),AINDEX(28),AINDEX(29),AINDEX(30), 1AINDEX(31),AINDEX(32),AINDEX(33),AINDEX(34),AINDEX(35), 1AINDEX(36),AINDEX(37),AINDEX(38),AINDEX(39),AINDEX(40), 1AINDEX(41),AINDEX(42),AINDEX(43),AINDEX(44),AINDEX(45), 1AINDEX(46),AINDEX(47),AINDEX(48),AINDEX(49),AINDEX(50) 1/26.,27.,28.,29.,30.,31.,32.,33.,34.,35.,36.,37.,38., 139.,40.,41.,42.,43.,44.,45.,46.,47.,48.,49.,50./ C IPR=6 IUPPER=7500 NUMDIS=43 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(N.LT.1.OR.N.GT.IUPPER)GOTO50 IF(N.EQ.1)GOTO55 HOLD=X(1) DO60I=2,N IF(X(I).NE.HOLD)GOTO90 60 CONTINUE WRITE(IPR, 9)HOLD RETURN 50 WRITE(IPR,17)IUPPER WRITE(IPR,47)N RETURN 55 WRITE(IPR,18) RETURN 90 CONTINUE 9 FORMAT(1H ,109H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT (A VECTOR) TO THE EXTREM SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6 1H *****) 17 FORMAT(1H , 98H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 EXTREM SUBROUTINE IS OUTSIDE THE ALLOWABLE (1,,I6,16H) INTERVAL * 1****) 18 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUME 1NT TO THE EXTREM SUBROUTINE HAS THE VALUE 1 *****) 47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8 ,6H *****) C C-----START POINT----------------------------------------------------- C AN=N C C COMPUTE THE SAMPLE MINIMUM AND SAMPLE MAXIMUM C XMIN=X(1) XMAX=X(1) DO140I=2,N IF(X(I).LT.XMIN)XMIN=X(I) IF(X(I).GT.XMAX)XMAX=X(I) 140 CONTINUE C C COMPUTE THE PROB PLOT CORRELATION COEFFICIENTS FOR THE VARIOUS VALUES C OF GAMMA C CALL SORT(X,N,Y) CALL UNIMED(N,Z) C DO100IDIS=1,NUMDIS IF(IDIS.EQ.NUMDIS)GOTO150 A=GAMTAB(IDIS) DO110I=1,N W(I)=(-ALOG(Z(I)))**(-1.0/A) 110 CONTINUE GOTO170 150 DO160I=1,N W(I)=-ALOG(ALOG(1.0/Z(I))) 160 CONTINUE C 170 SUM1=0.0 SUM2=0.0 DO200I=1,N SUM1=SUM1+Y(I) SUM2=SUM2+W(I) 200 CONTINUE YBAR=SUM1/AN WBAR=SUM2/AN SUM1=0.0 SUM2=0.0 SUM3=0.0 DO300I=1,N SUM2=SUM2+(Y(I)-YBAR)*(W(I)-WBAR) SUM1=SUM1+(Y(I)-YBAR)*(Y(I)-YBAR) SUM3=SUM3+(W(I)-WBAR)*(W(I)-WBAR) 300 CONTINUE SY=SQRT(SUM1/(AN-1.0)) CC=SUM2/SQRT(SUM3*SUM1) YSLOPE=SUM2/SUM3 YINT=YBAR-YSLOPE*WBAR CORR(IDIS)=CC YI(IDIS)=YINT YS(IDIS)=YSLOPE 100 CONTINUE C C DETERMINE THAT DISTRIBUTION WITH THE MAX PROB PLOT CORR COEFFICIENT C IDISMX=1 CORRMX=CORR(1) DO400IDIS=1,NUMDIS IF(CORR(IDIS).GT.CORRMX)IDISMX=IDIS IF(CORR(IDIS).GT.CORRMX)CORRMX=CORR(IDIS) 400 CONTINUE DO500IDIS=1,NUMDIS IFLAG1(IDIS)=BLANK IFLAG2(IDIS)=BLANK IFLAG3(IDIS)=BLANK IF(IDIS.EQ.IDISMX)GOTO550 GOTO500 550 IFLAG1(IDIS)=ALPHAM IFLAG2(IDIS)=ALPHAA IFLAG3(IDIS)=ALPHAX 500 CONTINUE C C WRITE OUT THE TABLE OF PROB PLOT CORR COEFFICIENTS FOR VARIOUS GAMMA C WRITE(IPR,998) WRITE(IPR,305) WRITE(IPR,999) WRITE(IPR,310)N WRITE(IPR,311)YBAR WRITE(IPR,312)SY WRITE(IPR,313)XMIN WRITE(IPR,314)XMAX WRITE(IPR,999) WRITE(IPR,323) WRITE(IPR,324) WRITE(IPR,325) WRITE(IPR,999) C NUMDM1=NUMDIS-1 IF(NUMDM1.LT.1)GOTO850 DO800I=1,NUMDM1 WRITE(IPR,805)GAMTAB(I),CORR(I),IFLAG1(I),IFLAG2(I),IFLAG3(I), 1YI(I),YS(I),T(I) 800 CONTINUE 850 I=NUMDIS WRITE(IPR,806)ALPHAI,ALPHAN,ALPHAF,ALPHAI,ALPHAN,ALPHAI, 1ALPHAT,ALPHAY,CORR(I),IFLAG1(I),IFLAG2(I),IFLAG3(I), 1YI(I),YS(I),T(I) C C PLOT THE PROB PLOT CORR COEFFICIENT VERSUS GAMMA VALUE INDEX C CALL PLOT(CORR,AINDEX,NUMDIS) WRITE(IPR,810)ALPHAG,ALPHAA,ALPHAM,ALPHAM,ALPHAA,EQUAL, 1GAMTAB(1),GAMTAB(12),GAMTAB(23),GAMTAB(34), 1ALPHAI,ALPHAN,ALPHAF,ALPHAI,ALPHAN,ALPHAI,ALPHAT,ALPHAY WRITE(IPR,999) WRITE(IPR,812) WRITE(IPR,813) C C IF THE OPTIMAL GAMMA IS FINITE, PLOT OUT THE EXTREME VALUE C TYPE 2 PROBABILITY PLOT FOR THE OPTIMAL VALUE C OF GAMMA. C IF(IDISMX.LT.NUMDIS)CALL EV2PLT(X,N,GAMTAB(IDISMX)) C C PLOT OUT AN EXTREME VALUE TYPE 1 PROBABILITY PLOT C CALL EV1PLT(X,N) C C FORM THE VARIOUS RETURN PERIOD VALUES C 1650 K=0 DO2100I=1,4 DO2200J=1,9 K=K+1 AM(K)=J*(10**(I-1)) 2200 CONTINUE 2100 CONTINUE K=K+1 AM(K)=10000. K=K+1 AM(K)=50000. K=K+1 AM(K)=100000. K=K+1 AM(K)=500000. K=K+1 AM(K)=1000000. K=K+1 AM(K)=N NUMAM=K CALL SORT(AM,NUMAM,SCRAT) DO2300I=1,NUMAM AM(I)=SCRAT(I) 2300 CONTINUE C C IF THE OPTIMAL GAMMA IS FINITE, COMPUTE THE C PREDICTED EXTREME (= F(1-(1/M)) FOR VARIOUS RETURN PERIODS M C FOR THE OPTIMAL EXTREME VALUE TYPE 2 DISTRIBUTION. C IF(IDISMX.EQ.NUMDIS)GOTO2450 A=GAMTAB(IDISMX) YINT=YI(IDISMX) YSLOPE=YS(IDISMX) DO2400I=2,NUMAM R=1.0/AM(I) P=1.0-R ARG=-ALOG(P) IF(ARG.LE.0.0)GOTO2400 H(I,1)=YINT+YSLOPE*(ARG**(-1.0/A)) 2400 CONTINUE C C COMPUTE THE PREDICTED EXTREME (= F(1-(1/M)) FOR VARIOUS RETURN C PERIODS M FOR THE EXTREME VALUE TYPE 1 DISTRIBUTION. C 2450 YINT=YI(NUMDIS) YSLOPE=YS(NUMDIS) DO2500I=2,NUMAM R=1.0/AM(I) P=1.0-R ARG=-ALOG(P) IF(ARG.LE.0.0)GOTO2500 H(I,2)=YINT+YSLOPE*(-ALOG(ARG)) 2500 CONTINUE C C WRITE OUT THE PAGE WITH THE RETURN PERIODS AND THE PREDICTED EXTREMES C FOR THE 2 DISTRIBUTIONS--OPTIMAL EXTREME VALUE TYPE 2, AND EXTREME C VALUE TYPE 1. C WRITE(IPR,998) IF(IDISMX.EQ.NUMDIS)GOTO2750 WRITE(IPR,2602) WRITE(IPR,2604) WRITE(IPR,2606) WRITE(IPR,2608) WRITE(IPR,2610)GAMTAB(IDISMX) WRITE(IPR,999) DO2700I=2,NUMAM WRITE(IPR,2705)AM(I),H(I,1),H(I,2) J=I-1 JSKIP=J-5*(J/5) IF(JSKIP.EQ.0)WRITE(IPR,999) 2700 CONTINUE RETURN C 2750 WRITE(IPR,2802) WRITE(IPR,2804) WRITE(IPR,2806) WRITE(IPR,2808) WRITE(IPR,999) DO2900I=2,NUMAM WRITE(IPR,2705)AM(I),H(I,2) J=I-1 JSKIP=J-5*(J/5) IF(JSKIP.EQ.0)WRITE(IPR,999) 2900 CONTINUE C 998 FORMAT(1H1) 999 FORMAT(1H ) 305 FORMAT(1H ,40X,22HEXTREME VALUE ANALYSIS) 310 FORMAT(1H ,37X,20HTHE SAMPLE SIZE N = ,I7) 311 FORMAT(1H ,34X,18HTHE SAMPLE MEAN = ,F14.7) 312 FORMAT(1H ,28X,32HTHE SAMPLE STANDARD DEVIATION = ,F14.7) 313 FORMAT(1H ,32X,21HTHE SAMPLE MINIMUM = ,F14.7) 314 FORMAT(1H ,32X,21HTHE SAMPLE MAXIMUM = ,F14.7) 323 FORMAT(1H ,85H EXTREME VALUE PROBABILITY PLOT LOCATIO 1N SCALE TAIL LENGTH) 324 FORMAT(1H ,83H TYPE 2 TAIL LENGTH CORRELATION ESTIMAT 1E ESTIMATE MEASURE) 325 FORMAT(1H ,37H PARAMETER (GAMMA) COEFFICIENT) 805 FORMAT(1H ,3X,F10.2,13X,F8.5,1X,3A1,2X,F14.7,2X,F14.7,3X,F10.5) 806 FORMAT(1H ,5X,8A1,13X,F8.5,1X,3A1,2X,F14.7,2X,F14.7,3X,F10.5) 810 FORMAT(1H ,12X,5A1,1X,A1,F14.7,11X,F14.7,11X,F14.7,11X,F14.7, 115X,8A1) 812 FORMAT(1H ,96HTHE ABOVE IS A PLOT OF THE 46 PROBABILITY PLOT CORRE 1LATION COEFFICIENTS (FROM THE PREVIOUS PAGE)) 813 FORMAT(1H ,16X,41HVERSUS THE 46 EXTREME VALUE DISTRIBUTIONS) 2602 FORMAT(1H ,43H RETURN PERIOD PREDICTED EXTREME WIND, 1 27H PREDICTED EXTREME WIND) 2604 FORMAT(1H ,43H (IN YEARS) BASED ON OPTIMAL , 1 20H BASED ON) 2606 FORMAT(1H ,42H EXTREME VALUE TYPE 2, 1 27H EXTREME VALUE TYPE 1) 2608 FORMAT(1H ,43H DISTRIBUTION , 1 22H DISTRIBUTION) 2610 FORMAT(1H ,30H (GAMMA = ,F12.5,1H)) 2705 FORMAT(1H ,2X,F9.1,13X,F10.2,17X,F10.2) 2802 FORMAT(1H ,43H RETURN PERIOD PREDICTED EXTREME WIND) 2804 FORMAT(1H ,36H (IN YEARS) BASED ON) 2806 FORMAT(1H ,42H EXTREME VALUE TYPE 1) 2808 FORMAT(1H ,38H DISTRIBUTION) C RETURN END