tail(3f) - [M_datapac:ANALYSIS] performs a symmetric distribution
tail length analysis
SUBROUTINE TAIL(X,N)
REAL(kind=wp) :: X(:)
INTEGER :: N
TAIL(3f) performs a symmetric distribution tail length analysis on
the data in the input vector X.
The analysis consists of the following--
1. Various test statistics to test
the specific hypothesis of normality;
2. A uniform probability plot
(a short-tailed distribution);
3. A normal probability plot
(a moderate-tailed distribution);
4. A tukey lambda = -0.5 probability plot
(a moderate-long-tailed distribution);
5. A cauchy probability plot
(a long-tailed distribution);
6. A determination of the best-fit
symmetric distribution
to the data set from an
admissible set consisting of
43 symmetric distributions.
The admissible set of symmetric distributions considered includes
the uniform, normal, logistic, double exponential, cauchy, and 37
distributions drawn from the the tukey lambda distributional family.
The goodness of fit criterion is the maximum probability plot
correlation coefficient criterion.
X The vector of unsorted or sorted) observations.
N The integer number of observations in the vector X.
The maximum allowable value of N for this subroutine is 3000.
6 pages of automatic printout--
1. various test statistics for normality;
2. a uniform probability plot;
3. a normal probability plot;
4. a tukey lambda = -0.5 probability plot;
5. a cauchy probability plot;
6. a determination of the best-fit symmetric distribution to the
data set.
Sample program:
program demo_tail
use M_datapac, only : tail, label
implicit none
real,allocatable :: x(:)
integer :: i
call label('tail')
x=[(real(i)/10.0,i=1,2000)]
x=x**3.78-6*x**2.52+9*x**1.26
call tail(x,size(x))
end program demo_tail
Results:
The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.
John Urban, 2022.05.31
CC0-1.0
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp) | :: | X(:) | ||||
| integer | :: | N |
| Type | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|
| real | :: | WS(15000) |
SUBROUTINE TAIL(X,N) REAL(kind=wp) :: X(:) INTEGER :: N REAL(kind=wp) :: a2, a3, a4, aa, ai, al, alamba, am2, am3, am4, an, arg, asub1, asubn, b1, b2, bb, bs, cc, coef REAL(kind=wp) :: coefi, constn, corr, corrmx, cox1xn, dd, del, eandev, eb1, eb2, ecc, ee, egeary, ei, er, ers, ersq, erssq, es, essq REAL(kind=wp) :: ewilks, ex1, ex1xn, exn, exnsq, g, gamma, geary, hold, P, p1, pi, picons, pn, ppfnor, PTEnth, q, rp1, rpn, rs REAL(kind=wp) :: s, sdb1, sdb2, sdcc, sdgear, sdrs, sdwilk, sfp1, sfpn, sum, sum1, sum2, sum3, sum4, varrs, varxn, wilksh, WS REAL(kind=wp) :: xbar REAL(kind=wp) :: xline, Y, YM, Z, zb1, zb2, zcc, zgeary, zrs, zwilks INTEGER :: i, icount, idis, idis2, idismx, ievodd, imax, imin, irev, iupper, mx, nhalf, nhalfp, nm1, numdis CHARACTER(len=4) :: iflag1 CHARACTER(len=4) :: iflag2 CHARACTER(len=4) :: iflag3 CHARACTER(len=4) :: iline1 CHARACTER(len=4) :: iline2 ! CHARACTER(len=4) :: alpham CHARACTER(len=4) :: alphaa CHARACTER(len=4) :: blank CHARACTER(len=4) :: hyphen CHARACTER(len=4) :: alphai CHARACTER(len=4) :: alphax character(len=256) :: message integer :: ios ! DIMENSION Y(3000) , Z(3000) , YM(3000) DIMENSION P(3000) , PTEnth(3000) DIMENSION corr(50) , iflag1(50) , iflag2(50) , iflag3(50) DIMENSION iline1(130) , iline2(130) DIMENSION xline(13) COMMON /BLOCK2_real32/ WS(15000) EQUIVALENCE (Y(1),WS(1)) EQUIVALENCE (Z(1),WS(3001)) EQUIVALENCE (YM(1),WS(6001)) EQUIVALENCE (P(1),WS(9001)) EQUIVALENCE (PTEnth(1),WS(12001)) ! DATA alpham , alphaa/'M' , 'A'/ DATA blank , hyphen , alphai , alphax/' ' , '-' , 'I' , 'X'/ DATA picons/3.14159265358979_wp/ DATA constn/.3989422804_wp/ ! iupper = 3000 ! ! 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 TAIL(3f) 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 TAIL(3f) 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 TAIL(3f) has all elements = ', & & E15.8,' *****') RETURN ENDIF ! !-----START POINT----------------------------------------------------- ! 50 an = N ! ! COMPUTE THE SAMPLE MEAN ! sum = 0.0_wp DO i = 1 , N sum = sum + X(i) ENDDO xbar = sum/an ! ! COMPUTE S, BIASED S, B1, AND B2 ! sum2 = 0.0_wp sum3 = 0.0_wp sum4 = 0.0_wp DO i = 1 , N del = X(i) - xbar a2 = del*del a3 = del*a2 a4 = a2*a2 sum2 = sum2 + a2 sum3 = sum3 + a3 sum4 = sum4 + a4 ENDDO am2 = sum2/an am3 = sum3/an am4 = sum4/an s = SQRT(sum2/(an-1.0_wp)) bs = SQRT(am2) b1 = am3/(bs**3) b2 = am4/(bs**4) ! ! COMPUTE THE EXPECTED VALUE AND STANDARD DEVIATION OF B1 AND B2 ! UNDER THE NORMALITY ASSUMPTION ! REFERENCE--CRAMER, page 386 ! eb1 = 0.0_wp sdb1 = 6.0_wp*(an-2.0_wp)/((an+1.0_wp)*(an+3.0_wp)) sdb1 = SQRT(sdb1) zb1 = (b1-eb1)/sdb1 eb2 = 3.0_wp - 6.0_wp/(an+1.0_wp) sdb2 = 24.0_wp*an*(an-2.0_wp)*(an-3.0_wp) / ((an+1.0_wp)*(an+1.0_wp)*(an+3.0_wp)*(an+5.0_wp)) zb2 = (b2-eb2)/sdb2 ! ! COMPUTE GEARY'S STATISTIC ! sum = 0.0_wp DO i = 1 , N sum = sum + ABS(X(i)-xbar) ENDDO eandev = sum/an geary = eandev/bs ! ! COMPUTE THE EXPECTED VALUE AND STANDARD DEVIATION ! OF GEARY'S STATISTIC UNDER THE NORMALITY ASSUMPTION ! REFERENCE--BIOMETRIKA, 1936, page 296 ! aa = SQRT(2.0_wp/picons) bb = SQRT(2.0_wp/(an-1.0)) IF ( N>=100 ) cc = SQRT(an/2.0_wp) * (1.0_wp-(1.0_wp/(8.0_wp*an/2.0_wp))+(1.0_wp/(128.0_wp*an*an/ 4.0_wp))) IF ( N<100 ) THEN coef = 1.0 imax = N - 1 ievodd = N - 2*(N/2) imin = 2 IF ( ievodd==0 ) imin = 3 IF ( imin<=imax ) THEN DO i = imin , imax , 2 ai = i coef = ((ai-1.0_wp)/ai)*coef ENDDO ENDIF coef = coef*(an-1.0_wp) IF ( ievodd==0 ) THEN coef = coef/SQRT(picons) ELSE coef = coef*(SQRT(picons)/2.0_wp) ENDIF cc = coef ENDIF egeary = aa/(bb*cc) dd = (2.0_wp/picons)*SQRT(an*(an-2.0_wp)) arg = 1.0_wp/(an-1.0_wp) arg = arg/SQRT(1.0_wp-arg*arg) ee = ATAN(arg) sdgear = (1.0_wp/an)*(1.0_wp+dd+ee) sdgear = sdgear - egeary*egeary sdgear = SQRT(sdgear) zgeary = (geary-egeary)/sdgear ! ! SORT THE DATA, ! THEN COMPUTE RANGE/S. ! CALL SORT(X,N,Y) rs = (Y(N)-Y(1))/s ! ! COMPUTE THE EXPECTED VALUE AND STANDARD DEVIATION OF THE RANGE/S ! UNDER THE NORMALITY ASSUMPTION ! REFERENCE--BIOMETRIKA, 1954, page 483 ! g = .33000598_wp + ((an-2.0_wp)**.16_wp)/41.785_wp pn = (an-g)/(an-2.0_wp*g+1.0_wp) p1 = 1.0_wp - pn CALL NORPPF(pn,rpn) CALL NORPPF(p1,rp1) exn = rpn ex1 = rp1 er = exn - ex1 CALL NORPPF(p1,ppfnor) sfp1 = 1.0_wp/(constn*EXP(-(ppfnor*ppfnor)/2.0_wp)) CALL NORPPF(pn,ppfnor) sfpn = 1.0_wp/(constn*EXP(-(ppfnor*ppfnor)/2.0_wp)) varxn = pn*(1.0_wp-pn)*sfpn*sfpn/(an+2.0_wp) exnsq = varxn + exn*exn cox1xn = p1*p1*sfp1*sfpn/(an+2.0_wp) ex1xn = cox1xn + ex1*exn ersq = 2.0_wp*(exnsq-ex1xn) es = bb*cc essq = 1.0_wp ers = er/es erssq = ersq/essq varrs = erssq - ers*ers sdrs = SQRT(varrs) zrs = (rs-ers)/sdrs ! ! COMPUTE THE WILK-SHAPIRO STATISTIC ! al = LOG10(an) gamma = .327511_wp + .058212_wp*al - .009776_wp*al*al sum = 0.0_wp IF ( N<=20 ) arg = N IF ( N>20 ) arg = N + 1 asubn = SQRT((1.0_wp+(1.0_wp/(4.0_wp*arg)))/SQRT(arg)) asub1 = -asubn sum = sum + asub1*Y(1) + asubn*Y(N) IF ( N>2 ) THEN nm1 = N - 1 DO i = 2 , nm1 ai = i pi = (ai-gamma)/(an-2.0_wp*gamma+1.0_wp) CALL NORPPF(pi,ei) coefi = 2.0_wp*ei/SQRT(-2.722_wp+4.083_wp*an) sum = sum + coefi*Y(i) ENDDO ENDIF wilksh = sum*sum/(an*bs*bs) ! ! COMPUTE THE EXPECTED VALUE AND STANDARD DEVIATION OF THE WILK-SHAPIRO ! STATISTIC UNDER THE NORMALITY ASSUMPTION ! REFERENCE--JJF APPROXIMATION TO MOMENTS ON page 601 OF BIOMETRIKA (1965) ! IF ( N==3 ) ewilks = .9135_wp IF ( N==4 ) ewilks = .9012_wp IF ( N>=5 ) ewilks = .9026_wp + (an-5.0_wp)/(44.608_wp+13.593_wp*SQRT(an)+10.267_wp*an) IF ( N==3 ) sdwilk = .0755_wp IF ( N==4 ) sdwilk = .0719_wp IF ( N>=5 ) sdwilk = .0670_wp + (an-5.0_wp)/(-42.368_wp-5.026_wp*SQRT(an)-14.925_wp*an) zwilks = (wilksh-ewilks)/sdwilk ! ! COMPUTE THE CORRELATION COEFFICIENT BETWEEN THE ORDERED OBSERVATIONS ! AND THE ORDER STATISIC MEDIANS FROM 44 DIFFERENT SYMMETRIC DISTRIBUTIONS ! numdis = 44 nhalf = N/2 nhalfp = nhalf + 1 ievodd = N - 2*(N/2) CALL UNIMED(N,Z) DO i = 1 , N PTEnth(i) = Z(i)**(0.1_wp) ENDDO DO idis = 1 , numdis IF ( idis==20 ) THEN DO i = 1 , nhalf irev = N - i + 1 CALL NORPPF(Z(i),YM(i)) YM(irev) = -YM(i) ENDDO IF ( ievodd==1 ) YM(nhalfp) = 0.0_wp ELSEIF ( idis==22 ) THEN DO i = 1 , nhalf irev = N - i + 1 YM(i) = LOG(Z(i)/(1.0_wp-Z(i))) YM(irev) = -YM(i) ENDDO IF ( ievodd==1 ) YM(nhalfp) = 0.0_wp ELSEIF ( idis==23 ) THEN DO i = 1 , nhalf irev = N - i + 1 IF ( Z(i)<=0.5_wp ) YM(i) = LOG(2.0_wp*Z(i)) IF ( Z(i)>0.5_wp ) YM(i) = -LOG(2.0_wp*(1.0_wp-Z(i))) YM(irev) = -YM(i) ENDDO IF ( ievodd==1 ) YM(nhalfp) = 0.0_wp ELSEIF ( idis==33 ) THEN DO i = 1 , nhalf irev = N - i + 1 arg = picons*Z(i) YM(i) = -COS(arg)/SIN(arg) YM(irev) = -YM(i) ENDDO IF ( ievodd==1 ) YM(nhalfp) = 0.0_wp ELSE IF ( idis<20 ) idis2 = idis IF ( idis==21 ) idis2 = idis - 1 IF ( 23<idis .AND. idis<33 ) idis2 = idis - 2 IF ( 33<idis ) idis2 = idis - 3 alamba = -(0.1_wp)*FLOAT(idis2) + 2.1_wp IF ( idis==1 ) THEN DO i = 1 , nhalf irev = N - i + 1 P(i) = Z(i)*Z(i) P(irev) = Z(irev)*Z(irev) YM(i) = (P(i)-P(irev))/alamba YM(irev) = -YM(i) ENDDO IF ( ievodd==1 ) YM(nhalfp) = 0.0_wp ELSEIF ( idis==11 ) THEN DO i = 1 , nhalf irev = N - i + 1 P(i) = Z(i) P(irev) = Z(irev) YM(i) = (P(i)-P(irev))/alamba YM(irev) = -YM(i) ENDDO IF ( ievodd==1 ) YM(nhalfp) = 0.0_wp ELSEIF ( idis==24 ) THEN DO i = 1 , nhalf irev = N - i + 1 P(i) = Z(i)**(-0.1_wp) P(irev) = Z(irev)**(-0.1_wp) YM(i) = (P(i)-P(irev))/alamba YM(irev) = -YM(i) ENDDO IF ( ievodd==1 ) YM(nhalfp) = 0.0_wp ELSEIF ( idis==34 ) THEN DO i = 1 , nhalf P(irev) = 1.0_wp/Z(irev) P(i) = 1.0_wp/Z(i) irev = N - i + 1 YM(i) = (P(i)-P(irev))/alamba YM(irev) = -YM(i) ENDDO IF ( ievodd==1 ) YM(nhalfp) = 0.0_wp ELSEIF ( idis==44 ) THEN DO i = 1 , nhalf irev = N - i + 1 P(i) = 1.0_wp/(Z(i)*Z(i)) P(irev) = 1.0_wp/(Z(irev)*Z(irev)) YM(i) = (P(i)-P(irev))/alamba YM(irev) = -YM(i) ENDDO IF ( ievodd==1 ) YM(nhalfp) = 0.0_wp ELSE DO i = 1 , nhalf irev = N - i + 1 P(i) = P(i)/PTEnth(i) P(irev) = P(irev)/PTEnth(irev) YM(i) = (P(i)-P(irev))/alamba YM(irev) = -YM(i) ENDDO IF ( ievodd==1 ) YM(nhalfp) = 0.0_wp ENDIF ENDIF sum1 = 0.0_wp sum2 = 0.0_wp DO i = 1 , N sum1 = sum1 + Y(i)*YM(i) sum2 = sum2 + YM(i)*YM(i) ENDDO sum2 = SQRT(sum2) sum3 = s*SQRT(an-1.0_wp) corr(idis) = sum1/(sum2*sum3) ENDDO ! ! DETERMINE THAT DISTRIBUTION WITH THE MAXIMUM 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 cc = corr(20) ! ! COMPUTE THE EXPECTED VALUE AND STANDARD DEVIATION OF THE PROBABILITY PLOT ! CORRELATION COEFFICIENT UNDER THE NORMALITY ASSUMPTION ! REFERENCE--JJF UNPUBLISHED MANUSCRIPT ! IF ( N==2 ) ecc = 1.0_wp IF ( N==3 ) ecc = .95492958_wp IF ( N>=4 ) ecc = .94947355_wp + (an-4.0_wp)/(196.815_wp-2.9418_wp*SQRT(an)+19.7916_wp*an) IF ( N==2 ) sdcc = 99999999.9999_wp IF ( N==3 ) sdcc = .04007697_wp IF ( N>=4 ) sdcc = .039492_wp + (an-4.0_wp)/(-127.0_wp-25.3_wp*an) zcc = (cc-ecc)/sdcc ! ! WRITE OUT THE NORMAL TAIL LENGTH STATISTICS page ! WRITE (G_IO,99044) WRITE (G_IO,99005) ! 99005 FORMAT (' ',48X,'TAIL LENGTH ANALYSIS') WRITE (G_IO,99045) WRITE (G_IO,99006) N 99006 FORMAT (' ',46X,'(THE SAMPLE SIZE N = ',I0,')') WRITE (G_IO,99007) xbar 99007 FORMAT (' ',40X,'(THE SAMPLE MEAN = ',E15.8,')') WRITE (G_IO,99008) s 99008 FORMAT (' ',35X,'(THE SAMPLE STANDARD DEVIATION = ',E15.8,')') WRITE (G_IO,99045) WRITE (G_IO,99009) 99009 FORMAT (' ',35X,'REFERENCE--SHAPIRO, WILK, AND CHEN, JASA, 1968, pages 1343-1372') WRITE (G_IO,99010) 99010 FORMAT (' ',35X,'REFERENCE--CRAMER, pages 386-387') WRITE (G_IO,99011) 99011 FORMAT (' ',35X,'REFERENCE--GEARY, BIOMETRIKA, 1947, pages 209-242') WRITE (G_IO,99012) 99012 FORMAT (' ',35X,'REFERENCE--BIOMETRIKA TABLES, VOLUME 1, pages 67-69, 59-60, 207-208, AND 200') WRITE (G_IO,99013) 99013 FORMAT (' ',35X,'REFERENCE--SHAPIRO AND WILK, BIOMETRIKA, 1965, pages 591-611') DO i = 1 , 6 WRITE (G_IO,99045) ENDDO WRITE (G_IO,99014) 99014 FORMAT (' ',49X,'TAIL LENGTH STATISTICS') WRITE (G_IO,99015) 99015 FORMAT (' ',5X,& & 'THE EXPECTED VALUE AND STANDARD DEVIATION OF STATISTICS ON THIS page ARE BASED ON THE NORMALITY ASSUMPTION') WRITE (G_IO,99045) WRITE (G_IO,99045) WRITE (G_IO,99016) 99016 FORMAT (' ',& & 'FORM OF STATISTIC VALUE OF STAT EXP(STAT) SD(STAT) (STAT-EXP(STAT))/SD(STAT) TABLE REFERENCE') WRITE (G_IO,99045) WRITE (G_IO,99017) b1 , eb1 , sdb1 , zb1 99017 FORMAT (' ','STANDARDIZED THIRD CENTRAL MOMENT ',F10.5,6X,F10.5,2X,F10.5,9X,F10.5,13X,'BIOMETRIKA TABLES') WRITE (G_IO,99018) 99018 FORMAT (' ',111X,'VOL. 1, page 207') WRITE (G_IO,99019) b2 , eb2 , sdb2 , zb2 99019 FORMAT (' ','STANDARDIZED FOURTH CENTRAL MOMENT ',F10.5,6X,F10.5,2X,F10.5,9X,F10.5,13X,'BIOMETRIKA TABLES') WRITE (G_IO,99020) 99020 FORMAT (' ',111X,'VOL. 1, page 208') WRITE (G_IO,99021) geary , egeary , sdgear , zgeary 99021 FORMAT (' ','GEARY STATISTIC (MEAN DEVIATION/S) ',F10.5,6X,F10.5,2X,F10.5,9X,F10.5,13X,'BIOMETRIKA TABLES') WRITE (G_IO,99022) 99022 FORMAT (' ',111X,'VOL. 1, page 207') WRITE (G_IO,99023) rs , ers , sdrs , zrs 99023 FORMAT (' ','RANGE/S ',F10.5,6X,F10.5,2X,F10.5,9X,F10.5,13X,'BIOMETRIKA TABLES') WRITE (G_IO,99024) 99024 FORMAT (' ',111X,'VOL. 1, page 200') WRITE (G_IO,99025) wilksh , ewilks , sdwilk , zwilks 99025 FORMAT (' ','WILK-SHAPIRO STATISTIC (BLUE FOR SCALE/S)',F10.5,6X,F10.5,2X,F10.5,9X,F10.5,13X,'BIOMETRIKA (1965)') WRITE (G_IO,99026) 99026 FORMAT (' ',111X,'page 605') WRITE (G_IO,99027) cc , ecc , sdcc , zcc 99027 FORMAT (' ','PROBABILITY PLOT CORRELATION COEFFICIENT ',F10.5, & & 6X,F10.5,2X,F10.5,9X,F10.5,13X,'UNPUBLISHED JJF') WRITE (G_IO,99028) 99028 FORMAT (' ',111X,'MANUSCRIPT') ! ! COMPUTE THE LINE PLOT WHICH SHOWS THE DISTRIBUTION OF THE OBSERVED ! VALUES IN TERMS OF MULTIPLES OF SAMPLE STANDARD DEVIATIONS AWAY FROM ! THE SAMPLE MEAN ! DO i = 1 , 130 iline1(i) = blank iline2(i) = blank ENDDO icount = 0 DO i = 1 , N mx = 10.0_wp*(((X(i)-xbar)/s)+6.0_wp) + 0.5_wp mx = mx + 7 IF ( mx<7 .OR. mx>127 ) icount = icount + 1 IF ( mx>=7 .AND. mx<=127 ) iline1(mx) = alphax ENDDO DO i = 7 , 127 iline2(i) = hyphen ENDDO DO i = 7 , 127 , 10 iline2(i) = alphai ENDDO xline(7) = xbar DO i = 1 , 6 irev = 13 - i + 1 ai = i xline(i) = xbar - (7.0_wp-ai)*s xline(irev) = xbar + (7.0_wp-ai)*s ENDDO ! ! WRITE OUT THE LINE PLOT SHOWING THE DEVIATIONS OF THE OBSERVATIONS ! ABOUT THE SAMPLE MEAN IN TERMS OF MULTIPLES OF THE SAMPLE STANDARD ! DEVIATION ! DO i = 1 , 8 WRITE (G_IO,99045) ENDDO WRITE (G_IO,99029) 99029 FORMAT (& & ' LINE PLOT SHOWING THE DISTRIBUTION OF THE OBSERVATIONS ABOUT THE SAMPLE MEAN ',& & 'IN TERMS OF MULTIPLES OF THE SAMPLE STANDARD DEVIATION') WRITE (G_IO,99045) WRITE (G_IO,99045) WRITE (G_IO,99042) (iline1(i),i=1,130) WRITE (G_IO,99042) (iline2(i),i=1,130) WRITE (G_IO,99030) 99030 FORMAT (' ', & &' -6 -5 -4 -3 -2 -1 & & 0 1 2 3 4 5 6') WRITE (G_IO,99031,iostat=ios,iomsg=message) (xline(i),i=1,13) 99031 FORMAT (' ',13F10.4) if(ios.ne.0)write(G_IO,'(A)')message(:len_trim(message)) WRITE (G_IO,99045) WRITE (G_IO,99032) icount 99032 FORMAT (' ',10X,I0, & &' OBSERVATIONS WERE IN EXCESS OF 6 SAMPLE STANDARD DEVIATIONS FROM& & THE SAMPLE MEAN AND SO WERE NOT PLOTTED') ! ! GENERATE UNIFORM, NORMAL, LAMBDA = -0.5, AND CAUCHY PROBABILITY PLOTS ! nhalf = (N/2) + 1 CALL PLOT(Y,Z,N) WRITE (G_IO,99033) N 99033 FORMAT (' ',35X,'UNIFORM PROBABILITY PLOT (THE SAMPLE SIZE N = ',I0,')') WRITE (G_IO,99043) corr(11) DO i = 1 , nhalf irev = N - i + 1 CALL NORPPF(Z(i),YM(i)) YM(irev) = -YM(i) ENDDO CALL PLOT(Y,YM,N) WRITE (G_IO,99034) N 99034 FORMAT (' ',35X,& & 'NORMAL PROBABILITY PLOT (THE SAMPLE SIZE N = ',I0, & & ')') WRITE (G_IO,99043) corr(20) alamba = -0.5_wp DO i = 1 , nhalf irev = N - i + 1 q = Z(i) YM(i) = (q**alamba-(1.0_wp-q)**alamba)/alamba YM(irev) = -YM(i) ENDDO CALL PLOT(Y,YM,N) WRITE (G_IO,99035) alamba , N 99035 FORMAT (' ',35X,'LAMBDA = ',F4.1,' PROBABILITY PLOT (THE SAMPLE SIZE N = ',I0,')') WRITE (G_IO,99043) corr(28) DO i = 1 , nhalf irev = N - i + 1 arg = picons*Z(i) YM(i) = -COS(arg)/SIN(arg) YM(irev) = -YM(i) ENDDO CALL PLOT(Y,YM,N) WRITE (G_IO,99036) N 99036 FORMAT (' ',35X,'CAUCHY PROBABILITY PLOT (THE SAMPLE SIZE N = ',I0,')') WRITE (G_IO,99043) corr(33) ! ! WRITE OUT THE PROBABILITY PLOT CORRELATION COEFFICIENT page ! WRITE (G_IO,99044) DO idis = 1 , numdis IF ( idis==20 ) THEN WRITE (G_IO,99037) N , corr(idis) , iflag1(idis) , iflag2(idis) , iflag3(idis) 99037 FORMAT (' ','THE CORRELATION BETWEEN THE ',I0, & &' ORDERED OBS. AND THE ORDER STAT. MEDIANS FROM THE NORMAL DISTRIBUTION IS ',F8.5,1X,3A1) ELSEIF ( idis==22 ) THEN WRITE (G_IO,99038) N , corr(idis) , iflag1(idis) , iflag2(idis) , iflag3(idis) 99038 FORMAT (' ','THE CORRELATION BETWEEN THE ',I0, & &' ORDERED OBS. AND THE ORDER STAT. MEDIANS FROM THE LOGISTIC DIST. IS ',F8.5,1X,3A1) ELSEIF ( idis==23 ) THEN WRITE (G_IO,99039) N , corr(idis) , iflag1(idis) , iflag2(idis) , iflag3(idis) 99039 FORMAT (' THE CORRELATION BETWEEN THE ',I0, & & ' ORDERED OBS. AND THE ORDER STAT. MEDIANS FROM THE DOUBLE EXP. DIST. IS ',F8.5,1X,3A1) ELSEIF ( idis==33 ) THEN WRITE (G_IO,99040) N , corr(idis) , iflag1(idis) , iflag2(idis) , iflag3(idis) 99040 FORMAT (' ','THE CORRELATION BETWEEN THE ',I0, & &' ORDERED OBS. AND THE ORDER STAT. MEDIANS FROM THE CAUCHY DISTRIBUTION IS ',F8.5,1X,3A1) ELSE IF ( idis<20 ) idis2 = idis IF ( idis==21 ) idis2 = idis - 1 IF ( 23<idis .AND. idis<33 ) idis2 = idis - 2 IF ( 33<idis ) idis2 = idis - 3 alamba = -(0.1)*FLOAT(idis2) + 2.1 WRITE (G_IO,99041) N, alamba, corr(idis), iflag1(idis), iflag2(idis), iflag3(idis) 99041 FORMAT (' ','THE CORRELATION BETWEEN THE ',I0, & & ' ORDERED OBS. AND THE ORDER STAT. MEDIANS FROM THE LAMBDA = '& & ,F4.1,' DIST. IS ',F8.5,1X,3A1) ENDIF ENDDO ENDIF 99042 FORMAT (' ',130A1) 99043 FORMAT (' ',34X,'(PROBABILITY PLOT CORRELATION COEFFICIENT = ', & & F8.5,')') 99044 FORMAT ('1') 99045 FORMAT (' ') ! END SUBROUTINE TAIL