WEIB Subroutine

public subroutine WEIB(X, N)

NAME

weib(3f) - [M_datapac:ANALYSIS] perform a Weibull distribution
analysis (Weibull PPCC analysis)

SYNOPSIS

   SUBROUTINE WEIB(X,N)

DESCRIPTION

WEIB(3f) performs a Weibull distribution analysis on the data in the
input vector X.

This analysis consists of determining that particular Weibull
distribution which best fits the data set.

The goodness of fit criterion is the maximum probability plot
correlation coefficient criterion.

After the best-fit distribution is determined, estimates are computed
and printed out for the location and scale parameters.

Two probability plots are also printed out-- the best-fit Weibull
probability plot and an extreme value type 1 probability plot (this
is due to the fact that as the Weibull parameter gamma approaches
infinity, the Weibull distribution approaches the extreme value type
1 distribution).

OPTIONS

 X   description of parameter
 Y   description of parameter

EXAMPLES

Sample program:

program demo_weib
use M_datapac, only : weib
implicit none
! call weib(x,y)
end program demo_weib

Results:

AUTHOR

The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.

MAINTAINER

John Urban, 2022.05.31

LICENSE

CC0-1.0

REFERENCE

  • Filliben (1972), ‘Techniques for Tail Length Analysis’, Proceedings of the Eighteenth Conference on the Design of Experiments in Army Research and Testing, pages 425-450.
  • Filliben, ‘The Percent Point Function’, UNpublished Manuscript.
  • Johnson and Kotz (1970), Continuous Univariate Distributions-1, pages 250-271.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), dimension(:) :: X
integer :: N

Common Blocks

CAUPLT (subroutine)
CHSCDF (subroutine)
CODE (subroutine)
DECOMP (subroutine)
DEMOD (subroutine)
DEXPLT (subroutine)
EV1PLT (subroutine)
EV2PLT (subroutine)
EXPPLT (subroutine)
EXTREM (subroutine)
EXTREM (subroutine)
FREQ (subroutine)
GAMPLT (subroutine)
GEOPLT (subroutine)
HFNPLT (subroutine)
INVXWX (subroutine)
LAMPLT (subroutine)
LGNPLT (subroutine)
LOGPLT (subroutine)
MEDIAN (subroutine)
norout (subroutine)
NORPLT (subroutine)
PARPLT (subroutine)
PLOTU (subroutine)
POIPLT (subroutine)
RUNS (subroutine)
SAMPP (subroutine)
SPCORR (subroutine)
TAIL (subroutine)
TPLT (subroutine)
TRIM (subroutine)
UNIPLT (subroutine)
WEIPLT (subroutine)
WIND (subroutine)
"> common /BLOCK2_real64/

Type Attributes Name Initial
real :: WS(15000)

Source Code

SUBROUTINE WEIB(X,N)
REAL(kind=wp) :: a , aindex , an , cc , corr , corrmx , gamtab , hold , sum1 ,&
     &     sum2 , sum3 , sy , t , w , wbar , WS , X , xmax , xmin , Y
REAL(kind=wp) :: ybar , yi , yint , ys , yslope , Z
INTEGER i , idis , idismx , iupper , N , numdis , numdm1
!
!     INPUT ARGUMENTS--X      = THE  VECTOR OF
!                               (UNSORTED OR SORTED) OBSERVATIONS.
!                      N      = THE INTEGER NUMBER OF OBSERVATIONS
!                               IN THE VECTOR X.
!     OUTPUT--4 pages OF AUTOMATIC PRINTOUT.
!     PRINTING--YES.
!     RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N
!                   FOR THIS SUBROUTINE IS 7500.
!
!---------------------------------------------------------------------
!
CHARACTER(len=4) :: iflag1
CHARACTER(len=4) :: iflag2
CHARACTER(len=4) :: iflag3
!
CHARACTER(len=4) :: blank
CHARACTER(len=4) :: alpham
CHARACTER(len=4) :: alphaa
CHARACTER(len=4) :: alphax
CHARACTER(len=4) :: alphai
CHARACTER(len=4) :: alphan
CHARACTER(len=4) :: alphaf
CHARACTER(len=4) :: alphat
CHARACTER(len=4) :: alphay
CHARACTER(len=4) :: alphag
CHARACTER(len=4) :: equal
!
DIMENSION w(3000)
DIMENSION X(:)
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_real64/ WS(15000)
EQUIVALENCE (Y(1),WS(1))
EQUIVALENCE (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) ,  &
     &     gamtab(6) , gamtab(7) , gamtab(8) , gamtab(9) , gamtab(10) , &
     &     gamtab(11) , gamtab(12) , gamtab(13) , gamtab(14) ,          &
     &     gamtab(15) , gamtab(16) , gamtab(17) , gamtab(18) ,          &
     &     gamtab(19) , gamtab(20) , gamtab(21) , gamtab(22) ,          &
     &     gamtab(23), gamtab(24), gamtab(25)/1.0_wp, 2.0_wp, 3.0_wp, 4.0_wp, 5.0_wp,&
     &     6.0_wp, 7.0_wp, 8.0_wp, 9.0_wp, 10.0_wp, 11.0_wp, 12.0_wp, 13.0_wp, 14.0_wp, 15.0_wp, 16.0_wp,&
     &     17.0_wp, 18.0_wp, 19.0_wp, 20.0_wp, 21.0_wp, 22.0_wp, 23.0_wp, 24.0_wp, 25.0_wp/
DATA gamtab(26) , gamtab(27) , gamtab(28) , gamtab(29) ,          &
     &     gamtab(30) , gamtab(31) , gamtab(32) , gamtab(33) ,          &
     &     gamtab(34) , gamtab(35) , gamtab(36) , gamtab(37) ,          &
     &     gamtab(38) , gamtab(39) , gamtab(40) , gamtab(41) ,          &
     &     gamtab(42)/30.0_wp, 35.0_wp, 40.0_wp, 45.0_wp, 50.0_wp, 60.0_wp, 70.0_wp, 80.0_wp,   &
     &     90.0_wp, 100.0_wp, 150.0_wp, 200.0_wp, 250.0_wp, 350.0_wp, 500.0_wp, 750.0_wp, 1000.0_wp/
DATA t(1) , t(2) , t(3) , t(4) , t(5) , t(6) , t(7) , t(8) ,      &
     &     t(9) , t(10) , t(11) , t(12) , t(13) , t(14) , t(15) ,       &
     &     t(16) , t(17) , t(18) , t(19) , t(20)/1.63474 , 1.36116 ,    &
     &     1.34278_wp , 1.35854_wp , 1.37836_wp , 1.39657_wp , 1.41225_wp , 1.42557_wp ,  &
     &     1.43690_wp , 1.44660_wp , 1.45496_wp , 1.46223_wp , 1.46860_wp , 1.47422_wp ,  &
     &     1.47921_wp , 1.48368_wp , 1.48769_wp , 1.49132_wp , 1.49461_wp , 1.49761_wp/
DATA t(21) , t(22) , t(23) , t(24) , t(25) , t(26) , t(27) ,      &
     &     t(28) , t(29) , t(30) , t(31) , t(32) , t(33) , t(34) ,      &
     &     t(35) , t(36) , t(37) , t(38) , t(39) , t(40) , t(41) ,      &
     &     t(42) , t(43)/1.50036_wp , 1.50288_wp , 1.50521_wp , 1.50736_wp ,        &
     &     1.50935_wp , 1.51748_wp , 1.52344_wp , 1.52798_wp , 1.53157_wp , 1.53447_wp ,  &
     &     1.53888_wp , 1.54206_wp , 1.54447_wp , 1.54636_wp , 1.54788_wp , 1.55248_wp ,  &
     &     1.55480_wp , 1.55620_wp , 1.55781_wp , 1.55902_wp , 1.55997_wp , 1.56044_wp ,  &
     &     1.62391_wp/
DATA aindex(1) , aindex(2) , aindex(3) , aindex(4) , aindex(5) ,  &
     &     aindex(6) , aindex(7) , aindex(8) , aindex(9) , aindex(10) , &
     &     aindex(11) , aindex(12) , aindex(13) , aindex(14) ,          &
     &     aindex(15) , aindex(16) , aindex(17) , aindex(18) ,          &
     &     aindex(19) , aindex(20) , aindex(21) , aindex(22) ,          &
     &     aindex(23), aindex(24), aindex(25)/1.0_wp, 2.0_wp, 3.0_wp, 4.0_wp, 5.0_wp,&
     &     6.0_wp, 7.0_wp, 8.0_wp, 9.0_wp, 10.0_wp, 11.0_wp, 12.0_wp, 13.0_wp, 14.0_wp, 15.0_wp, 16.0_wp,&
     &     17.0_wp, 18.0_wp, 19.0_wp, 20.0_wp, 21.0_wp, 22.0_wp, 23.0_wp, 24.0_wp, 25.0_wp/
DATA aindex(26) , aindex(27) , aindex(28) , aindex(29) ,          &
     &     aindex(30) , aindex(31) , aindex(32) , aindex(33) ,          &
     &     aindex(34) , aindex(35) , aindex(36) , aindex(37) ,          &
     &     aindex(38) , aindex(39) , aindex(40) , aindex(41) ,          &
     &     aindex(42) , aindex(43) , aindex(44) , aindex(45) ,          &
     &     aindex(46) , aindex(47) , aindex(48) , aindex(49) ,          &
     &     aindex(50)/26.0_wp, 27.0_wp, 28.0_wp, 29.0_wp, 30.0_wp, 31.0_wp, 32.0_wp, 33.0_wp,   &
     &     34.0_wp, 35.0_wp, 36.0_wp, 37.0_wp, 38.0_wp, 39.0_wp, 40.0_wp, 41.0_wp, 42.0_wp, 43.0_wp,  &
     &     44.0_wp, 45.0_wp, 46.0_wp, 47.0_wp, 48.0_wp, 49.0_wp, 50.0_wp/
!
      iupper = 7500
      numdis = 43
!
!     CHECK THE INPUT ARGUMENTS FOR ERRORS
!
      IF ( N<1 .OR. N>iupper ) THEN
         WRITE (G_IO,99001) iupper
99001    FORMAT (' ',                                                   &
     &'***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE WEIB   SUBROU&
     &TINE IS OUTSIDE THE ALLOWABLE (1,',I0,') INTERVAL *****')
         WRITE (G_IO,99002) N
99002    FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,' *****')
         RETURN
      ELSE
         IF ( N==1 ) THEN
            WRITE (G_IO,99003)
99003       FORMAT (' ',                                                &
     &'***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUMENT TO THE WEIB&
     &   SUBROUTINE HAS THE VALUE 1 *****')
            RETURN
         ELSE
            hold = X(1)
            DO i = 2 , N
               IF ( X(i)/=hold ) GOTO 50
            ENDDO
            WRITE (G_IO,99004) hold
99004       FORMAT (' ',                                                &
     &'***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUMENT (A VECTOR) &
     &TO THE WEIB   SUBROUTINE HAS ALL ELEMENTS = ',E15.8,' *****')
            RETURN
         ENDIF
!
!-----START POINT-----------------------------------------------------
!
 50      an = N
!
!     COMPUTE THE SAMPLE MINIMUM AND SAMPLE MAXIMUM
!
         xmin = X(1)
         xmax = X(1)
         DO i = 2 , N
            IF ( X(i)<xmin ) xmin = X(i)
            IF ( X(i)>xmax ) xmax = X(i)
         ENDDO
!
!     COMPUTE THE PROB PLOT CORRELATION COEFFICIENTS FOR THE VARIOUS VALUES
!     OF GAMMA
!
         CALL SORT(X,N,Y)
         CALL UNIMED(N,Z)
!
         DO idis = 1 , numdis
            IF ( idis==numdis ) THEN
               DO i = 1 , N
                  w(i) = -LOG(LOG(1.0_wp/Z(i)))
               ENDDO
            ELSE
               a = gamtab(idis)
               DO i = 1 , N
                  w(i) = (-LOG(1.0_wp-Z(i)))**(1.0_wp/a)
               ENDDO
            ENDIF
!
            sum1 = 0.0_wp
            sum2 = 0.0_wp
            DO i = 1 , N
               sum1 = sum1 + Y(i)
               sum2 = sum2 + w(i)
            ENDDO
            ybar = sum1/an
            wbar = sum2/an
            sum1 = 0.0_wp
            sum2 = 0.0_wp
            sum3 = 0.0_wp
            DO i = 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)
            ENDDO
            sy = SQRT(sum1/(an-1.0_wp))
            cc = sum2/SQRT(sum3*sum1)
            yslope = sum2/sum3
            yint = ybar - yslope*wbar
            corr(idis) = cc
            yi(idis) = yint
            ys(idis) = yslope
         ENDDO
!
!     DETERMINE THAT DISTRIBUTION WITH THE MAX PROB PLOT CORR COEFFICIENT
!
         idismx = 1
         corrmx = corr(1)
         DO idis = 1 , numdis
            IF ( corr(idis)>corrmx ) idismx = idis
            IF ( corr(idis)>corrmx ) corrmx = corr(idis)
         ENDDO
         DO idis = 1 , numdis
            iflag1(idis) = blank
            iflag2(idis) = blank
            iflag3(idis) = blank
            IF ( idis==idismx ) THEN
               iflag1(idis) = alpham
               iflag2(idis) = alphaa
               iflag3(idis) = alphax
            ENDIF
         ENDDO
!
!     WRITE OUT THE TABLE OF PROB PLOT CORR COEFFICIENTS FOR VARIOUS GAMMA
!
         WRITE (G_IO,99005)
!
99005    FORMAT ('1')
         WRITE (G_IO,99006)
99006    FORMAT (' ',40X,'WEIBULL ANALYSIS')
         WRITE (G_IO,99020)
         WRITE (G_IO,99007) N
99007    FORMAT (' ',37X,'THE SAMPLE SIZE N = ',I0)
         WRITE (G_IO,99008) ybar
99008    FORMAT (' ',34X,'THE SAMPLE MEAN = ',F14.7)
         WRITE (G_IO,99009) sy
99009    FORMAT (' ',28X,'THE SAMPLE STANDARD DEVIATION = ',F14.7)
         WRITE (G_IO,99010) xmin
99010    FORMAT (' ',32X,'THE SAMPLE MINIMUM = ',F14.7)
         WRITE (G_IO,99011) xmax
99011    FORMAT (' ',32X,'THE SAMPLE MAXIMUM = ',F14.7)
         WRITE (G_IO,99020)
         WRITE (G_IO,99012)
99012    FORMAT (' ',                                                   &
     &'       WEIBULL          PROBABILITY PLOT     LOCATION         SCA&
     &LE       TAIL LENGTH')
         WRITE (G_IO,99013)
99013    FORMAT (' ',                                                   &
     &'     TAIL LENGTH          CORRELATION        ESTIMATE        ESTI&
     &MATE       MEASURE')
         WRITE (G_IO,99014)
99014    FORMAT (' ','   PARAMETER (GAMMA)      COEFFICIENT')
         WRITE (G_IO,99020)
!
         numdm1 = numdis - 1
         IF ( numdm1>=1 ) THEN
            DO i = 1 , numdm1
               WRITE (G_IO,99015) gamtab(i) , corr(i) , iflag1(i) ,      &
     &                           iflag2(i) , iflag3(i) , yi(i) , ys(i) ,&
     &                           t(i)
99015          FORMAT (' ',3X,F10.2,13X,F8.5,1X,3A1,2X,F14.7,2X,F14.7,  &
     &                 3X,F10.5)
            ENDDO
         ENDIF
         i = numdis
         WRITE (G_IO,99016) alphai , alphan , alphaf , alphai , alphan , &
     &                     alphai , alphat , alphay , corr(i) ,         &
     &                     iflag1(i) , iflag2(i) , iflag3(i) , yi(i) ,  &
     &                     ys(i) , t(i)
99016    FORMAT (' ',5X,8A1,13X,F8.5,1X,3A1,2X,F14.7,2X,F14.7,3X,F10.5)
!
!     PLOT THE PROB PLOT CORR COEFFICIENT VERSUS GAMMA VALUE INDEX
!
         CALL PLOT(corr,aindex,numdis)
         WRITE (G_IO,99017) alphag , alphaa , alpham , alpham , alphaa , &
     &                     equal , gamtab(1) , gamtab(12) , gamtab(23) ,&
     &                     gamtab(34) , alphai , alphan , alphaf ,      &
     &                     alphai , alphan , alphai , alphat , alphay
99017    FORMAT (' ',12X,5A1,1X,A1,F14.7,11X,F14.7,11X,F14.7,11X,F14.7, &
     &           15X,8A1)
         WRITE (G_IO,99020)
         WRITE (G_IO,99018)
99018    FORMAT (' ',                                                   &
     &'THE ABOVE IS A PLOT OF THE 46 PROBABILITY PLOT CORRELATION COEFFI&
     &CIENTS (FROM THE PREVIOUS page)')
         WRITE (G_IO,99019)
99019    FORMAT (' ',16X,'VERSUS THE 46 WEIBULL DISTRIBUTIONS')
!
!     IF THE OPTIMAL GAMMA IS FINITE, PLOT OUT THE WEIBULL
!     PROBABILITY PLOT FOR THE OPTIMAL VALUE
!     OF GAMMA.
!
         IF ( idismx<numdis ) CALL WEIPLT(X,N,gamtab(idismx))
!
!     PLOT OUT AN EXTREM VALUE TYPE 1 PROBABILITY PLOT
!     (WHICH IS IDENTICALLY A WEIBULL PROBABILITY
!     WITH GAMMA = INFINITY)
!
         CALL EV1PLT(X,N)
      ENDIF
99020 FORMAT (' ')
!
END SUBROUTINE WEIB