SUBROUTINE EV2PLT(X,N,GAMMA) C C PURPOSE--THIS SUBROUTINE GENERATES A EXTREME VALUE TYPE 2 C PROBABILITY PLOT C (WITH TAIL LENGTH PARAMETER VALUE = GAMMA). C THE PROTOTYPE EXTREME VALUE TYPE 2 DISTRIBUTION USED N C HEREIN IS DEFINED FOR ALL NON-NEGATIVE X, C AND HAS THE PROBABILITY DENSITY FUNCTION C F(X) = GAMMA * (X**(-GAMMA-1)) * EXP(-(X**(-GAMMA))). C AS USED HEREIN, A PROBABILITY PLOT FOR A DISTRIBUTION C IS A PLOT OF THE ORDERED OBSERVATIONS VERSUS C THE ORDER STATISTIC MEDIANS FOR THAT DISTRIBUTION. C THE EXTREME VALUE TYPE 2 PROBABILITY PLOT IS USEFUL IN C GRAPHICALLY TESTING THE COMPOSITE (THAT IS, C LOCATION AND SCALE PARAMETERS NEED NOT BE SPECIFIED) C HYPOTHESIS THAT THE UNDERLYING DISTRIBUTION C FROM WHICH THE DATA HAVE BEEN RANDOMLY DRAWN C IS THE EXTREME VALUE TYPE 2 DISTRIBUTION C WITH TAIL LENGTH PARAMETER VALUE = GAMMA. C IF THE HYPOTHESIS IS TRUE, THE PROBABILITY PLOT C SHOULD BE NEAR-LINEAR. C A MEASURE OF SUCH LINEARITY IS GIVEN BY THE C CALCULATED PROBABILITY PLOT CORRELATION COEFFICIENT. 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 --GAMMA = THE SINGLE PRECISION VALUE OF THE C TAIL LENGTH PARAMETER. C GAMMA SHOULD BE POSITIVE. C OUTPUT--A ONE-PAGE EXTREME VALUE TYPE 2 PROBABILITY PLOT. C PRINTING--YES. C RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N C FOR THIS SUBROUTINE IS 7500. C --GAMMA SHOULD BE POSITIVE. C OTHER DATAPAC SUBROUTINES NEEDED--SORT, UNIMED, PLOT. C FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, ALOG. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C REFERENCES--FILLIBEN, 'TECHNIQUES FOR TAIL LENGTH ANALYSIS', C PROCEEDINGS OF THE EIGHTEENTH CONFERENCE C ON THE DESIGN OF EXPERIMENTS IN ARMY RESEARCH C DEVELOPMENT AND TESTING (ABERDEEN, MARYLAND, C OCTOBER, 1972), PAGES 425-450. C --HAHN AND SHAPIRO, STATISTICAL METHODS IN ENGINEERING, C 1967, PAGES 260-308. C --JOHNSON AND KOTZ, 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--DECEMBER 1972. C UPDATED --SEPTEMBER 1975. C UPDATED --NOVEMBER 1975. C UPDATED --FEBRUARY 1976. C C--------------------------------------------------------------------- C DIMENSION X(1) DIMENSION Y(7500),W(7500) COMMON /BLOCK2/ WS(15000) EQUIVALENCE (Y(1),WS(1)),(W(1),WS(7501)) C IPR=6 IUPPER=7500 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(N.LT.1.OR.N.GT.IUPPER)GOTO50 IF(N.EQ.1)GOTO55 IF(GAMMA.LE.0.0)GOTO60 HOLD=X(1) DO65I=2,N IF(X(I).NE.HOLD)GOTO90 65 CONTINUE WRITE(IPR, 9)HOLD RETURN 50 WRITE(IPR,17)IUPPER WRITE(IPR,47)N RETURN 55 WRITE(IPR,18) RETURN 60 WRITE(IPR,25) WRITE(IPR,46)GAMMA RETURN 90 CONTINUE 9 FORMAT(1H ,109H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT (A VECTOR) TO THE EV2PLT SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6 1H *****) 17 FORMAT(1H , 98H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 EV2PLT SUBROUTINE IS OUTSIDE THE ALLOWABLE (1,,I6,16H) INTERVAL * 1****) 18 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUME 1NT TO THE EV2PLT SUBROUTINE HAS THE VALUE 1 *****) 25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE 1 EV2PLT SUBROUTINE IS NON-POSITIVE *****) 46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****) 47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8 ,6H *****) C C-----START POINT----------------------------------------------------- C AN=N C C SORT THE DATA C CALL SORT(X,N,Y) C C GENERATE UNIFORM ORDER STATISTIC MEDIANS C CALL UNIMED(N,W) C C COMPUTE EXREME VALUE TYPE 2 DISTRIBUTION ORDER STATISTIC MEDIANS C DO100I=1,N W(I)=(-ALOG(W(I)))**(-1.0/GAMMA) 100 CONTINUE C C PLOT THE ORDERED OBSERVATIONS VERSUS ORDER STATISTICS MEDIANS. C COMPUTE THE TAIL LENGTH MEASURE OF THE DISTRIBUTION. C WRITE OUT THE TAIL LENGTH MEASURE OF THE DISTRIBUTION C AND THE SAMPLE SIZE. C CALL PLOT(Y,W,N) Q=.9975 PP9975=(-ALOG(Q))**(-1.0/GAMMA) Q=.0025 PP0025=(-ALOG(Q))**(-1.0/GAMMA) Q=.975 PP975 =(-ALOG(Q))**(-1.0/GAMMA) Q=.025 PP025 =(-ALOG(Q))**(-1.0/GAMMA) TAU=(PP9975-PP0025)/(PP975-PP025) WRITE(IPR,105)GAMMA,TAU,N C C COMPUTE THE PROBABILITY PLOT CORRELATION COEFFICIENT. C COMPUTE LOCATION AND SCALE ESTIMATES C FROM THE INTERCEPT AND SLOPE OF THE PROBABILITY PLOT. C THEN WRITE THEM OUT. C 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 SUM1=SUM1+(Y(I)-YBAR)*(Y(I)-YBAR) SUM2=SUM2+(Y(I)-YBAR)*(W(I)-WBAR) SUM3=SUM3+(W(I)-WBAR)*(W(I)-WBAR) 300 CONTINUE CC=SUM2/SQRT(SUM3*SUM1) YSLOPE=SUM2/SUM3 YINT=YBAR-YSLOPE*WBAR WRITE(IPR,305)CC,YINT,YSLOPE C 105 FORMAT(1H ,63HEXTREME VALUE TYPE 2 (CAUCHY TYPE) PROB. PLOT WITH E 1XP. PAR. = ,E17.10,1X,7H(TAU = ,E15.8,1H),1X,16HSAMPLE SIZE N = ,I 17) 305 FORMAT(1H ,43HPROBABILITY PLOT CORRELATION COEFFICIENT = ,F8.5,5X, 122HESTIMATED INTERCEPT = ,E15.8,3X,18HESTIMATED SLOPE = ,E15.8) C RETURN END