SUBROUTINE GAMPLT(X,N,GAMMA) C C PURPOSE--THIS SUBROUTINE GENERATES A GAMMA C PROBABILITY PLOT C (WITH TAIL LENGTH PARAMETER VALUE = GAMMA). C THE PROTOTYPE GAMMA DISTRIBUTION USED C HEREIN HAS MEAN = GAMMA C AND STANDARD DEVIATION = SQRT(GAMMA). C THIS DISTRIBUTION IS DEFINED FOR ALL POSITIVE X, C AND HAS THE PROBABILITY DENSITY FUNCTION C F(X) = (1/CONSTANT) * (X**(GAMMA-1)) * EXP(-X) C WHERE THE CONSTANT = THE GAMMA FUNCTION EVALUATED C AT THE VALUE 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 GAMMA 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 GAMMA 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 GAMMA 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, ABS, EXP, DEXP, DLOG. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION AND DOUBLE PRECISION C LANGUAGE--ANSI FORTRAN. C REFERENCES--WILK, GNANADESIKAN, AND HUYETT, 'PROBABILITY C PLOTS FOR THE GAMMA DISTRIBUTION', C TECHNOMETRICS, 1962, PAGES 1-15. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 257, FORMULA 6.1.41. C --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 166-206. 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--NOVEMBER 1974. C UPDATED --SEPTEMBER 1975. C UPDATED --NOVEMBER 1975. C UPDATED --FEBRUARY 1976. C C--------------------------------------------------------------------- C DOUBLE PRECISION Z,Z2,Z3,Z4,Z5,DEN,A,B,C,D DOUBLE PRECISION DEXP,DLOG DIMENSION D(10) DIMENSION X(1) DIMENSION Y(7500),W(7500) COMMON /BLOCK2/ WS(15000) EQUIVALENCE (Y(1),WS(1)),(W(1),WS(7501)) DATA C/ .918938533204672741D0/ DATA D(1),D(2),D(3),D(4),D(5) 1 /+.833333333333333333D-1,-.277777777777777778D-2, 1+.793650793650793651D-3,-.595238095238095238D-3,+.8417508417508417 151D-3/ DATA D(6),D(7),D(8),D(9),D(10) 1 /-.191752691752691753D-2,+.641025641025641025D-2,-.2955065359 147712418D-1,+.179644372368830573D0,-.139243221690590111D1/ 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 GAMPLT SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6 1H *****) 17 FORMAT(1H , 98H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 GAMPLT SUBROUTINE IS OUTSIDE THE ALLOWABLE (1,,I6,16H) INTERVAL * 1****) 18 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUME 1NT TO THE GAMPLT SUBROUTINE HAS THE VALUE 1 *****) 25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE 1 GAMPLT 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 DGAMMA=GAMMA C C COMPUTE THE GAMMA FUNCTION USING THE ALGORITHM IN THE C NBS APPLIED MATHEMATICS SERIES REFERENCE. C THIS GAMMA FUNCTION NEED BE CALCULATED ONLY ONCE. C IT IS USED IN THE CALCULATION OF THE CDF BASED ON C THE TENTATIVE VALUE OF THE PPF IN THE ITERATION. C Z=DGAMMA DEN=1.0D0 150 IF(Z.GE.10.0D0)GOTO160 DEN=DEN*Z Z=Z+1.0D0 GOTO150 160 Z2=Z*Z Z3=Z*Z2 Z4=Z2*Z2 Z5=Z2*Z3 A=(Z-0.5D0)*DLOG(Z)-Z+C B=D(1)/Z+D(2)/Z3+D(3)/Z5+D(4)/(Z2*Z5)+D(5)/(Z4*Z5)+ 1D(6)/(Z*Z5*Z5)+D(7)/(Z3*Z5*Z5)+D(8)/(Z5*Z5*Z5)+D(9)/(Z2*Z5*Z5*Z5) G=DEXP(A+B)/DEN 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 GENERATE GAMMA DISTRIBUTION ORDER STATISTIC MEDIANS C C DETERMINE LOWER AND UPPER BOUNDS ON THE DESIRED I-TH GAMMA C ORDER STATISTIC MEDIAN. C FOR EACH I, A LOWER BOUND IS GIVEN BY C (Y(I)*GAMMA*THE GAMMA FUNCTION OF GAMMA)**(1.0/GAMMA) C WHERE Y(I) IS THE CORRESPONDING UNIFORM (0,1) ORDER STATISIC C MEDIAN. C FOR EACH I EXCEPT I = N, AN UPPER BOUND IS GIVEN BY THE C (I+1)-ST GAMMA ORDER STATISTIC MEDIAN (ASSUMEDLY ALREADY C CALCULTATED). C FOR I = N, AN UPPER BOUND IS DETERMINED BY COMPUTING C MULTIPLES OF THE LOWER BOUND FOR I = N UNTIL A LARGER C VALUE IS OBTAINED. C DUE TO THE ABOVE CONSIDERATIONS, THE GAMMA ORDER STATISTIC C MEDIANS WILL BE CALCULATED LARGEST TO SMALLEST, THAT IS, C IN THE FOLLOWING SEQUENCE: W(N), W(N-1), ..., W(2), W(1). C NOTE ALSO THAT 1) THE CODE IS COMPLICATED SLIGHTLY BY THE C FACT THAT PERCENT POINT VALUES INVOLVED IN THE CALCULATION OF C THE TAIL LENGTH MEASURE TAU (SEE LABEL 605) ARE GOING ON C 'SIMULATNEOUSLY'. AND 2) THE VECTOR W WILL AT VARIOUS TIMES C IN THE PROGRAM HAVE UNIFORM ORDER STATISTIC MEDIANS AND C THEN LATER GRADUALLY FILL UP WITH GAMMA ORDER STATISTIC C MEDIANS. C I=N ITAIL=0 310 IF(ITAIL.EQ.0)U=W(I) DP=U XMIN0=(U*GAMMA*G)**(1.0/GAMMA) XMIN=XMIN0 IF(I.EQ.N.OR.ITAIL.GE.1)GOTO320 IP1=I+1 XMAX=W(IP1) GOTO370 320 ILOOP=1 ICOUNT=1 350 ACOUNT=ICOUNT XMAX=ACOUNT*XMIN0 DX=XMAX GOTO1000 360 IF(PCALC.GE.DP)GOTO370 XMIN=XMAX ICOUNT=ICOUNT+1 IF(ICOUNT.LE.30000)GOTO350 370 XMID=(XMIN+XMAX)/2.0 C C AT THIS STAGE WE NOW HAVE LOWER AND UPPER LIMITS ON C THE DESIRED I-TH GAMMA ORDER STATISITC MEDIAN W(I). C NOW ITERATE BY BISECTION UNTIL THE DESIRED ACCURACY IS ACHIEVED C FOR THE I-TH GAMMA ORDER STATISITIC MEDIAN. C ILOOP=2 XLOWER=XMIN XUPPER=XMAX ICOUNT=0 550 DX=XMID GOTO1000 560 IF(PCALC.EQ.DP)GOTO570 IF(PCALC.GT.DP)GOTO580 XLOWER=XMID XMID=(XMID+XUPPER)/2.0 GOTO590 580 XUPPER=XMID XMID=(XMID+XLOWER)/2.0 590 XDEL=ABS(XMID-XLOWER) ICOUNT=ICOUNT+1 IF(XDEL.LT.0.0000001.OR.ICOUNT.GT.100)GOTO570 GOTO550 570 IF(ITAIL.GE.1)GOTO605 W(I)=XMID IF(I.LE.1)GOTO595 I=I-1 GOTO310 595 CONTINUE C C AT THIS POINT, THE GAMMA ORDER STATISTIC MEDIANS ARE ALL COMPUTED. C NOW PLOT OUT THE GAMMA PROBABILITY PLOT C CALL PLOT(Y,W,N) C 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 605 IF(ITAIL.EQ.0)GOTO600 IF(ITAIL.EQ.1)GOTO610 IF(ITAIL.EQ.2)GOTO620 IF(ITAIL.EQ.3)GOTO630 GOTO640 600 U=.9975 ITAIL=1 GOTO310 610 PP9975=XMID U=.0025 ITAIL=2 GOTO310 620 PP0025=XMID U=.975 ITAIL=3 GOTO310 630 PP975=XMID U=.025 ITAIL=4 GOTO310 640 PP025=XMID TAU=(PP9975-PP0025)/(PP975-PP025) WRITE(IPR,655)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 DO660I=1,N SUM1=SUM1+Y(I) SUM2=SUM2+W(I) 660 CONTINUE YBAR=SUM1/AN WBAR=SUM2/AN SUM1=0.0 SUM2=0.0 SUM3=0.0 DO670I=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) 670 CONTINUE CC=SUM2/SQRT(SUM3*SUM1) YSLOPE=SUM2/SUM3 YINT=YBAR-YSLOPE*WBAR WRITE(IPR,675)CC,YINT,YSLOPE C 655 FORMAT(1H ,46HGAMMA PROBABILITY PLOT WITH SHAPE PARAMETER = , 1E17.10,1X,7H(TAU = ,E15.8,1H),16X,16HSAMPLE SIZE N = ,I7) 675 FORMAT(1H ,43HPROBABILITY PLOT CORRELATION COEFFICIENT = ,F8.5,5X, 122HESTIMATED INTERCEPT = ,E15.8,3X,18HESTIMATED SLOPE = ,E15.8) C RETURN C C******************************************************************** C THIS SECTION BELOW IS LOGICALLY SEPARATE FROM THE ABOVE. C THIS SECTION COMPUTES A CDF VALUE FOR ANY GIVEN TENTATIVE C PERCENT POINT X VALUE AS DEFINED IN EITHER OF THE 2 C ITERATION LOOPS IN THE ABOVE CODE. C C COMPUTE T-SUB-Q AS DEFINED ON PAGE 4 OF THE WILK, GNANADESIKAN, C AND HUYETT REFERENCE C 1000 SUM=1.0/DGAMMA TERM=1.0/DGAMMA CUT1=DX-DGAMMA CUT2=DX*10000000.0 DO700J=1,1000 AJ=J TERM=DX*TERM/(DGAMMA+AJ) SUM=SUM+TERM CUTOFF=CUT1+(CUT2*TERM/SUM) IF(AJ.GT.CUTOFF)GOTO750 700 CONTINUE WRITE(IPR,705) WRITE(IPR,707)GAMMA 705 FORMAT(1H ,48H*****ERROR IN INTERNAL OPERATIONS IN THE GAMPLT , 1 53HSUBROUTINE--THE NUMBER OF CDF ITERATIONS EXCEEDS 1000) 707 FORMAT(1H ,33H THE INPUT VALUE OF GAMMA IS ,E15.8) 750 T=SUM PCALC=(DX**DGAMMA)*(EXP(-DX))*T/G IF(ILOOP.EQ.1)GOTO360 GOTO560 C END