SUBROUTINE WEIB(X,N) C C PURPOSE--THIS SUBROUTINE PERFOMS A WEIBULL DISTRIBUTION ANALYSIS C ON THE DATA IN THE INPUT VECTOR X. C THIS ANALYSIS CONSISTS OF DETERMINING THAT PARTICULAR C WEIBULL 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 WEIBULL PROBABILITY PLOT C AND AN EXTREME VALUE TYPE 1 PROBABILITY PLOT C (THIS IS DUE TO THE FACT THAT AS THE WEIBULL PARAMETER C GAMMA APPROACHES INFINITY, THE WEIBULL DISTRIBUTION C APPROACHES THE EXTREME VALUE TYPE 1 DISTRIBUTION). 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--4 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, WEIPLT, C EV1PLT, PLOT. C FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT AND ALOG. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C REFERENCE--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, PAGES 250-271. 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 --AUGUST 1975. C UPDATED --NOVEMBER 1975. C UPDATED --MAY 1976. C C--------------------------------------------------------------------- C CHARACTER*4 IFLAG1 CHARACTER*4 IFLAG2 CHARACTER*4 IFLAG3 C CHARACTER*4 BLANK CHARACTER*4 ALPHAM CHARACTER*4 ALPHAA CHARACTER*4 ALPHAX CHARACTER*4 ALPHAI CHARACTER*4 ALPHAN CHARACTER*4 ALPHAF CHARACTER*4 ALPHAT CHARACTER*4 ALPHAY CHARACTER*4 ALPHAG CHARACTER*4 EQUAL 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) DIMENSION AINDEX(50) 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./ 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) 1/1.63474,1.36116,1.34278,1.35854,1.37836,1.39657,1.41225,1.42557, 1 1.43690,1.44660, 1 1.45496,1.46223,1.46860,1.47422,1.47921,1.48368,1.48769,1.49132, 1 1.49461,1.49761/ DATA T(21),T(22),T(23),T(24),T(25),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.50036,1.50288,1.50521,1.50736,1.50935,1.51748,1.52344,1.52798, 1 1.53157,1.53447, 1 1.53888,1.54206,1.54447,1.54636,1.54788,1.55248,1.55480,1.55620, 1 1.55781,1.55902, 1 1.55997,1.56044,1.62391/ 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 WEIB SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6 1H *****) 17 FORMAT(1H , 98H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 WEIB SUBROUTINE IS OUTSIDE THE ALLOWABLE (1,,I6,16H) INTERVAL * 1****) 18 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUME 1NT TO THE WEIB 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(1.0-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 WEIBULL C PROBABILITY PLOT FOR THE OPTIMAL VALUE C OF GAMMA. C IF(IDISMX.LT.NUMDIS)CALL WEIPLT(X,N,GAMTAB(IDISMX)) C C PLOT OUT AN EXTREM VALUE TYPE 1 PROBABILITY PLOT C (WHICH IS IDENTICALLY A WEIBULL PROBABILITY C WITH GAMMA = INFINITY) C CALL EV1PLT(X,N) C 998 FORMAT(1H1) 999 FORMAT(1H ) 305 FORMAT(1H ,40X,16HWEIBULL 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 WEIBULL PROBABILITY PLOT LOCATIO 1N SCALE TAIL LENGTH) 324 FORMAT(1H ,83H 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,35HVERSUS THE 46 WEIBULL DISTRIBUTIONS) C RETURN END