SUBROUTINE DPMLAD(Y,N, 1XTEMP,MAXNXT,DALPHA,DBETA,DH,ITEMP, 1AKML,ALOCML,SCALEML, 1ICAPSW,ICAPTY, 1ISUBRO,IBUGA3,IERROR) C C PURPOSE--THIS ROUTINE COMPUTES THE MAXIMUM LIKELIHOOD C ESTIMATES FOR THE ASYMMETRIC DOUBLE EXPONENTIAL C DISTRIBUTION. THE ALGORITHM IS: C C 1) SORT THE DATA C C 2) COMPUTE C C h(x(j)) = 2*LOG[alpha(theta)) + SQRT(beta(theta))] + C SQRT(alpha(theta)*SQRT(beta(theta)) C C WHERE C C alpha(theta) = (1/N)*SUM[j=1 to N][(x(j) - theta)+] C beta(theta) = (1/N)*SUM[j=1 to N][(x(j) - theta)-] C C WHERE C C (x(j) - theta)+ = x(j) - theta x(j) >= theta C = 0 x(j) < theta C (x(j) - theta)- = theta - x(j) x(j) <= theta C = 0 x(j) > theta C C 3) SET R EQUAL TO THE VALUE OF J WHERE H(x(j)) HAS C IT'S MINIMUM VALUE. C C 4) IF R=1 OR R=N, THE MAXIMUM LIKELIHOOD ESTIMATES C DO NOT EXIST. HOWEVER, THESE CASES SUGGEST C POSITIVE AND NEGATIVE EXPONENTIAL DISTRIBUTIONS, C RESPECTIVELY. C C 5) OTHERWISE, THE MAXIMUM LIKELIHOOD ESTIMATES ARE: C C THETAHAT = X(R) C C KHAT = (BETA(THETAHAT))**(1/4)/ C (ALPHA(THETAHAT))**(1/4) C C SIGMAHAT = SQRT(2)*(BETA(THETAHAT))**(1/4)* C (ALPHA(THETAHAT))**(1/4)* C (SQRT(ALPHA(THETAHAT)) + C SQRT(BETA(THETAHAT))) C C EXAMPLE--ASYMMETRIC DOUBLE EXPONENTIAL MAXIMUM LIKELIHOOD Y C REFERENCES--KOTZ, KOZUBOWSKI, AND PODGORSKI, "THE LAPLACE C DISTRIBUTION AND GENERALIZATIONS: A REVISIT WITH C APPLICATIONS TO COMMUNICATIONS, ECONOMICS, C ENGINEERING, AND FINANCE", BIRKHAUSR, 2001, C PP. 133-178. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING DIVISION C INFORMATION TECHNOLOGY LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899-8980 C PHONE--301-975-2855 C NOTE--DATAPLOT IS A REGISTERED TRADEMARK C OF THE NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY. C THIS SUBROUTINE MAY NOT BE COPIED, EXTRACTED, C MODIFIED, OR OTHERWISE USED IN A CONTEXT C OUTSIDE OF THE DATAPLOT LANGUAGE/SYSTEM. C LANGUAGE--ANSI FORTRAN (1977) C VERSION NUMBER--2004/8 C ORIGINAL VERSION--AUGUST 2004. C C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES------------------- C DOUBLE PRECISION DSUM C CHARACTER*4 ICAPSW CHARACTER*4 ICAPTY CHARACTER*4 ISUBRO CHARACTER*4 IBUGA3 CHARACTER*4 IERROR C CHARACTER*4 IWRITE CHARACTER*1 IBASLC C CHARACTER*4 ISUBN1 CHARACTER*4 ISUBN2 CHARACTER*4 ISTEPN C C--------------------------------------------------------------------- C DIMENSION Y(*) DIMENSION XTEMP(*) DIMENSION ITEMP(*) C DOUBLE PRECISION DALPHA(*) DOUBLE PRECISION DBETA(*) DOUBLE PRECISION DH(*) C DOUBLE PRECISION DHMIN DOUBLE PRECISION DSUM1 DOUBLE PRECISION DSUM2 DOUBLE PRECISION DTERM1 DOUBLE PRECISION DTERM2 DOUBLE PRECISION DTERM3 DOUBLE PRECISION DTERM4 C DIMENSION FISH(3,3) DIMENSION COV(3,3) C INCLUDE 'DPCOF2.INC' C CHARACTER*4 ISUBN0 C CHARACTER*80 IFILE1 CHARACTER*12 ISTAT1 CHARACTER*12 IFORM1 CHARACTER*12 IACCE1 CHARACTER*12 IPROT1 CHARACTER*12 ICURS1 CHARACTER*4 IERRF1 CHARACTER*4 IENDF1 CHARACTER*4 IREWI1 C C--------------------------------------------------------------------- C CHARACTER*4 IFEEDB CHARACTER*4 IPRINT CHARACTER*240 ICOUT C COMMON /MACH/IRD,IPR,CPUMIN,CPUMAX,NUMBPC,NUMCPW,NUMBPW COMMON /PRINT/IFEEDB,IPRINT COMMON /TEXTOU/ICOUT,NCOUT,ILOUT C C-----START POINT----------------------------------------------------- C ISUBN1='DPML' ISUBN2='DE ' C IERROR='NO' C IF(IBUGA3.EQ.'ON'.OR.ISUBRO.EQ.'MLAD')THEN WRITE(ICOUT,999) 999 FORMAT(1X) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,51) 51 FORMAT('**** AT THE BEGINNING OF DPMLAD--') CALL DPWRST('XXX','WRIT') WRITE(ICOUT,52)IBUGA3 52 FORMAT('IBUGA3 = ',A4) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,55)N 55 FORMAT('N = ',I8) CALL DPWRST('XXX','WRIT') DO56I=1,MIN(N,100) WRITE(ICOUT,57)I,Y(I) 57 FORMAT('I,Y(I) = ',I8,E15.7) CALL DPWRST('XXX','WRIT') 56 CONTINUE ENDIF C C ******************************************** C ** STEP 11-- ** C ** CHECK THE INPUT ARGUMENTS FOR ERRORS ** C ******************************************** C ISTEPN='11' IF(IBUGA3.EQ.'ON'.OR.ISUBRO.EQ.'MLAD') 1CALL TRACE2(ISTEPN,ISUBN1,ISUBN2) C IF(N.LE.1)THEN WRITE(ICOUT,999) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,1111) 1111 FORMAT('***** ERROR IN ASYMMETRIC DOUBLE EXPONENTIAL MAXIMUM ', 1 'LIKELIHOOD--') CALL DPWRST('XXX','WRIT') WRITE(ICOUT,1113) 1113 FORMAT(' THE NUMBER OF OBSERVATIONS ', 1 'FOR VARIABLE 1 IS NON-POSITIVE') CALL DPWRST('XXX','WRIT') WRITE(ICOUT,1115)N 1115 FORMAT('SAMPLE SIZE = ',I8) CALL DPWRST('XXX','WRIT') IERROR='YES' GOTO9000 ENDIF C HOLD=Y(1) DO1135I=2,N IF(Y(I).NE.HOLD)GOTO1139 1135 CONTINUE 1130 CONTINUE WRITE(ICOUT,999) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,1131) 1131 FORMAT('***** ERROR FROM ASYMMETRIC DOUBLE EXPONENTIAL MAXIMUM ', 1 'LIKELIHOOD--') CALL DPWRST('XXX','WRIT') WRITE(ICOUT,1133)HOLD 1133 FORMAT(' VARIABLE 1 HAS ALL ELEMENTS = ',E15.7) CALL DPWRST('XXX','WRIT') GOTO9000 1139 CONTINUE C C ********************************** C ** STEP 21-- ** C ** CARRY OUT CALCULATIONS ** C ** FOR ASYMMETRIC DOUBLE EXPONENTIAL MLE ** C ** ESTIMATE (FULL SAMPLE CASE) ** C ********************************** C 2100 CONTINUE C ISTEPN='21' IF(IBUGA3.EQ.'ON'.OR.ISUBRO.EQ.'MLAD') 1CALL TRACE2(ISTEPN,ISUBN1,ISUBN2) C IERROR='NO' IWRITE='OFF' C CALL SORT(Y,N,XTEMP) DO2105I=1,N Y(I)=XTEMP(I) 2105 CONTINUE XMIN=Y(1) XMAX=Y(N) C CALL MEDIAN(Y,N,IWRITE,XTEMP,MAXNXT,XMED,IBUGA3,IERROR) CALL MEAN(Y,N,IWRITE,XMEAN,IBUGA3,IERROR) CALL SD(Y,N,IWRITE,XSD,IBUGA3,IERROR) C DHMIN=DBLE(CPUMAX) DO2110I=1,N THETA=Y(I) DSUM1=0.0D0 DSUM2=0.0D0 DO2120J=1,N IF(Y(J).GE.THETA)DSUM1=DSUM1 + DBLE(Y(J) - THETA) IF(Y(J).LE.THETA)DSUM2=DSUM2 + DBLE(THETA - Y(J)) 2120 CONTINUE DALPHA(I)=DSUM1/DBLE(N) DBETA(I)=DSUM2/DBLE(N) DH(I)=2.0D0*DLOG(DSQRT(DALPHA(I)) + DSQRT(DBETA(I))) 1 + DSQRT(DALPHA(I))*DSQRT(DBETA(I)) IF(DH(I).LT.DHMIN)THEN DHMIN=DH(I) IR=I ENDIF 2110 CONTINUE C IF(IR.EQ.1)THEN WRITE(ICOUT,999) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2131) 2131 FORMAT('***** ERROR FROM ASYMMETRIC DOUBLE EXPONENTIAL ', 1 'MAXIMUM LIKELIHOOD--') CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2133) 2133 FORMAT(' ESTIMATE OF LOCATION PARAMTER EQUALS DATA ', 1 'MINIMUM. THE MAXIMUM') CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2135) 2135 FORMAT(' LIKELIHOOD ESTIMATES DO NOT EXIST. HOWEVER, ') CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2137) 2137 FORMAT(' THIS IMPLIES THAT AN EXPONENTIAL MODEL IS ', 1 'APPROPRIATE.') CALL DPWRST('XXX','WRIT') IERROR='YES' AKML=-999.0 SCALEML=1.0 ALOCML=0.0 GOTO9000 ELSEIF(IR.EQ.N)THEN WRITE(ICOUT,999) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2131) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2143) 2143 FORMAT(' ESTIMATE OF LOCATiON PARAMTER EQUALS DATA ', 1 'MAXIMUM. THE MAXIMUM') CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2145) 2145 FORMAT(' LIKELIHOOD ESTIMATES DO NOT EXIST. HOWEVER, ') CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2147) 2147 FORMAT(' THIS IMPLIES THAT A NEGATIVE EXPONENTIAL ', 1 'MODEL IS APPROPRIATE.') CALL DPWRST('XXX','WRIT') IERROR='YES' AKML=-999.0 SCALEML=1.0 ALOCML=0.0 GOTO9000 ELSE ALOCML=Y(IR) ENDIF C DTERM1=DBETA(IR)**(1.0D0/4.0D0) DTERM2=DALPHA(IR)**(1.0D0/4.0D0) DTERM3=DSQRT(DBETA(IR)) DTERM4=DSQRT(DALPHA(IR)) IF(DTERM2.LE.0.0D0)THEN WRITE(ICOUT,999) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2131) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2153) 2153 FORMAT(' INFINITE VALUE FOR THE SHAPE PARAMETER.') CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2155) 2155 FORMAT(' THE CAUSE OF THIS IS TIES FOR THE ', 1 'MAXIMUM DATA VALUE AND') CALL DPWRST('XXX','WRIT') WRITE(ICOUT,2157) 2157 FORMAT(' THE ESTIMATE OF THE LOCATION PARAMETER ', 1 'OCCURS AT THE DATA MAXIMUM.') CALL DPWRST('XXX','WRIT') IERROR='YES' GOTO9000 ENDIF AKML=DTERM1/DTERM2 SCALEML=DSQRT(2.0D0)*DTERM2*DTERM1*(DTERM3 + DTERM4) C C NOW COMPUTE THE FISHER INFORMATION MATRIX, THEN INVERT TO C OBTAIN THE ASYMPTOTIC VARIANCE-COVARIANCE MATRIX. THEN WRITE C TO DPST1F.DAT. C FISH(1,1)=2.0/(SCALEML**2) FISH(1,2)=-SQRT(2.0)*2.0/(SCALEML*(1.0+AKML**2)) FISH(3,1)=0.0 FISH(2,1)=FISH(1,2) FISH(2,2)=(1.0/(AKML**2)) + 4.0/((1.0+AKML**2)**2) FISH(2,3)=-(1.0-AKML**2)/(SCALEML*AKML*(1.0+AKML**2)) FISH(1,3)=FISH(3,1) FISH(3,2)=FISH(2,3) FISH(3,3)=1.0/(SCALEML**2) CCCCC print *,'fish(1,1) = ',fish(1,1) CCCCC print *,'fish(2,2) = ',fish(2,2) CCCCC print *,'fish(3,3) = ',fish(3,3) CCCCC print *,'fish(1,2) = ',fish(1,2) CCCCC print *,'fish(2,3) = ',fish(2,3) CCCCC print *,'fish(3,1) = ',fish(3,1) CALL SGECO(FISH,3,3,ITEMP,RCOND,XTEMP) IJOB=1 CALL SGEDI(FISH,3,3,ITEMP,XTEMP,XTEMP(MAXNXT/2),IJOB) DO2810J=1,3 DO2815I=1,3 COV(I,J)=FISH(I,J) 2815 CONTINUE 2810 CONTINUE C IOUNI1=IST1NU IFILE1=IST1NA ISTAT1=IST1ST IFORM1=IST1FO IACCE1=IST1AC IPROT1=IST1PR ICURS1=IST1CS ISUBN0='MLAD' IERRF1='NO' C IREWI1='ON' CALL DPOPFI(IOUNI1,IFILE1,ISTAT1,IFORM1,IACCE1,IPROT1,ICURS1, 1IREWI1,ISUBN0,IERRF1,IBUGA3,ISUBRO,IERROR) IF(IERRF1.EQ.'YES')GOTO9000 C WRITE(IOUNI1,2905) 2905 FORMAT(' THETA K SIGMA') DO2910I=1,3 WRITE(IOUNI1,'(3(E15.7,2X))')COV(I,1),COV(I,2),COV(I,3) 2910 CONTINUE C IENDF1='OFF' IREWI1='ON' CALL DPCLFI(IOUNI1,IFILE1,ISTAT1,IFORM1,IACCE1,IPROT1,ICURS1, 1IENDF1,IREWI1,ISUBN0,IERRF1,IBUGA3,ISUBRO,IERROR) C IF(IPRINT.EQ.'ON')THEN IF(ICAPSW.EQ.'ON' .AND. ICAPTY.EQ.'HTML')THEN C C STEP 2: START TABLE AND DEFINE A CAPTION C 5001 FORMAT('') 5011 FORMAT('') 5099 FORMAT('
')
        WRITE(ICOUT,5091)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5093)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5099)
        CALL DPWRST('XXX','WRIT')
C
      ELSEIF(ICAPSW.EQ.'ON' .AND. ICAPTY.EQ.'LATE')THEN
C
C  STEP 1: END ASIS FORMAT, START TABLE ENVIRONMENT, WRITE A HEADER, AND
C          WRITE A TABLE CAPTION
C
 8001 FORMAT(A1,'end{verbatim}')
 8003 FORMAT(A1,'begin{table}')
 8007 FORMAT(5X,'$ ',A1,1X,A1,' $ ',A1,A1,' ')
 8009 FORMAT(A1,'begin{center}')
 8011 FORMAT(5X,'{',A1,'bf Asymmetric Exponential Maximum ',
     1       'Likelihood Estimate}')
 8013 FORMAT(A1,'end{center}')
 8015 FORMAT(5X,'} ',A1,A1)
C
        CALL DPCONA(92,IBASLC)
C
        WRITE(ICOUT,8001)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8003)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8009)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8011)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8007)IBASLC,IBASLC,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8007)IBASLC,IBASLC,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8013)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
C
C STEP 2: START TABULAR ENVIRONMENT, WRITE SUBSESQUENT ROWS, END
C         TABULAR ENVIRONMENT
C
 8020 FORMAT(5X,A1,'begin{tabular} {lr}')
 8021 FORMAT(5X,'Number of Observations: & ',I8,2X,A1,A1)
 8022 FORMAT(5X,'Sample Mean: & ',G15.7,2X,A1,A1)
 8023 FORMAT(5X,'Sample Median: & ',G15.7,2X,A1,A1)
 8024 FORMAT(5X,'Sample Standard Deviation: & ',G15.7,2X,A1,A1)
 8025 FORMAT(5X,'Estimate of Location Parameter ($',A1,'theta$): & ',
     1       G15.7,2X,A1,A1)
 8026 FORMAT(5X,'Estimate of Scale Parameter ($,',A1,'sigma$): & ',
     1       G15.7,2X,A1,A1)
 8027 FORMAT(5X,'Estimate of Shape Parameter ($,',A1,'kappa$): & ',
     1       G15.7,2X,A1,A1)
 8028 FORMAT(5X,'Sample Minimum: & ',G15.7,2X,A1,A1)
 8029 FORMAT(5X,'Sample Maximum: & ',G15.7,2X,A1,A1)
 8030 FORMAT(5X,' & ',2X,A1,A1)
 8040 FORMAT(5X,A1,'hline')
 8049 FORMAT(A1,'end{tabular}')
        WRITE(ICOUT,8009)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8020)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8021)N,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8022)XMEAN,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8023)XMED,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8024)XSD,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8028)XMIN,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8029)XMAX,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8030)IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8025)IBASLC,ALOCML,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8026)IBASLC,SCALEML,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8027)IBASLC,AKML,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8049)IBASLC
        CALL DPWRST('XXX','WRIT')
C
C STEP 3: END TABLE ENVIRONMENT, RESET ASIS MODE
C
 8091 FORMAT(A1,'end{center}')
 8093 FORMAT(A1,'end{table}')
 8099 FORMAT(A1,'begin{verbatim}')
        WRITE(ICOUT,8091)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8093)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8099)IBASLC
        CALL DPWRST('XXX','WRIT')
      ELSEIF(ICAPSW.EQ.'ON' .AND. ICAPTY.EQ.'RTF')THEN
      ELSE
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4105)
 4105   FORMAT(6X,'ASYMMETRIC DOUBLE EXPONENTIAL MAXIMUM LIKELIHOOD ',
     1         'ESTIMATION')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4109)N
 4109   FORMAT('THE NUMBER OF OBSERVATIONS                 = ',I8)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4111)XMEAN
 4111   FORMAT('THE SAMPLE MEAN                            = ',G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4112)XMED  
 4112   FORMAT('THE SAMPLE MEDIAN                          = ',G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4113)XSD
 4113   FORMAT('THE SAMPLE STANDARD DEVIATION              = ',G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4115)XMIN
 4115   FORMAT('THE SAMPLE MINIMUM                         = ',G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4117)XMAX
 4117   FORMAT('THE SAMPLE MAXIMUM                         = ',G15.7)
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,4121)ALOCML
 4121   FORMAT('ESTIMATE OF THE LOCATION PARAMETER (THETA) = ',
     1         G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4123)SCALEML
 4123   FORMAT('ESTIMATE OF THE SCALE PARAMETER (SIGMA)    = ',
     1         G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4125)AKML
 4125   FORMAT('ESTIMATE OF THE SHAPE PARAMETER (AK)       = ',
     1         G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,4131)
 4131   FORMAT('THE LOCATION, SCALE, AND SHAPE PARAMETERS WILL ',
     1         'BE SAVED AS')
        CALL DPWRST('XXX','BUG ')
        WRITE(ICOUT,4133)
 4133   FORMAT('THE INTERNAL PARAMETERS LOCML, SCALEML, AND ',
     1         'KML.')
        CALL DPWRST('XXX','BUG ')
        WRITE(ICOUT,4135)
 4135   FORMAT('THE ASYMPTOTIC PARAMETER VARIANCE-COVARIANCE ')
        CALL DPWRST('XXX','BUG ')
        WRITE(ICOUT,4137)
 4137   FORMAT('WRITTEN TO THE FILE dpst1f.dat.')
        CALL DPWRST('XXX','BUG ')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','BUG ')
      ENDIF
      ENDIF
C
C               *****************
C               **  STEP 90--  **
C               **  EXIT       **
C               *****************
C
 9000 CONTINUE
      IF(IBUGA3.EQ.'OFF'.AND.ISUBRO.NE.'MLAD')GOTO9090
      WRITE(ICOUT,999)
      CALL DPWRST('XXX','WRIT')
      WRITE(ICOUT,9011)
 9011 FORMAT('***** AT THE END       OF DPMLAD--')
      CALL DPWRST('XXX','WRIT')
      WRITE(ICOUT,9012)N,IBUGA3,IERROR
 9012 FORMAT('N,IBUGA3,IERROR = ',I8,2X,A4,2X,A4)
      CALL DPWRST('XXX','WRIT')
      WRITE(ICOUT,9015)N
 9015 FORMAT('N = ',I8)
      CALL DPWRST('XXX','WRIT')
 9090 CONTINUE
C
      RETURN
      END
      SUBROUTINE DPMLAE(Y,X,N,NVAR,
     1TEMP1,TEMP2,TEMP3,DTEMP1,
     1THETMO,PMOM,THETFR,PFR,THETAF2,PF2,THETML,PML,
     1ICAPSW,ICAPTY,MAXNXT,
     1ISUBRO,IBUGA3,IERROR)
C
C     PURPOSE--THIS ROUTINE COMPUTES THE MAXIMUM LIKELIHOOD
C              ESTIMATES FOR THE POLYA-AEPPLI DISTRIBUTION.
C
C              THE MOMENT ESTIMATORS ARE:
C
C                  THETAHAT = 2*XBAR**2/(s**2+XBAR)
C                  PHAT     = (S**2 - XBAR)/(S**2 + XBAR)
C
C              THE MEAN AND ZERO FREQUENCY ESTIMATORS ARE:
C
C                  THETAHAT = LOG(f0/N)
C                  PHAT     = 1 - THETHAT/XBAR
C
C              THE FIRST TWO FREQUENCIES ESTIMATORS ARE:
C
C                  THETAHAT = -LOG(f0/N)
C                  PHAT     = -f1/{f0*LOG(f0/N)}
C
C              THE MAXIMUM LIKELIHOOD ESTIMATES ARE THE SOLUTION
C              THE FOLLOWING EQUATIONS:
C
C                  XBAR - THETAHAT/(1-PHAT) = 0
C                  XBAR - SUM[J=1 to N][fj*(J-1)*P(J-1)/(N*P(J))} = 0
C
C              WHERE P(J) = THE POLYA-AEPPLI PDF USING THE
C              ESTIMATED PARAMETERS.
C
C              THERE ARE TWO CASES:
C
C              1) ONE VARIABLE CASE: Y IS RAW DATA
C              2) TWO VARIABLE CASE: Y IS FREQUENCY, X IS CLASS
C                 MID-POINT.
C
C     EXAMPLE--POLYA-AEPPLI MAXIMUM LIKELIHOOD Y
C            --POLYA-AEPPLI MAXIMUM LIKELIHOOD Y X
C     REFERENCES--JOHNSON, KOTZ, AND KEMP (1992).  "UNIVARIATE
C                 DISCRETE DISTRIBUTIONS", SECOND EDITION, 
C                 WILEY, PP. 378-382.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBUG, MD 20899-8980
C                 PHONE--301-975-2855
C     NOTE--DATAPLOT IS A REGISTERED TRADEMARK
C           OF THE NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY.
C           THIS SUBROUTINE MAY NOT BE COPIED, EXTRACTED,
C           MODIFIED, OR OTHERWISE USED IN A CONTEXT
C           OUTSIDE OF THE DATAPLOT LANGUAGE/SYSTEM.
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--2006/6
C     ORIGINAL VERSION--JUNE      2006.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-----------------
C
      CHARACTER*4 ICAPSW
      CHARACTER*4 ICAPTY
      CHARACTER*4 ISUBRO
      CHARACTER*4 IBUGA3
      CHARACTER*4 IERROR
C
      CHARACTER*4 IWRITE
C
      CHARACTER*4 ISUBN1
      CHARACTER*4 ISUBN2
      CHARACTER*4 ISTEPN
C
      CHARACTER*1 IBASLC
C
      CHARACTER*4 ISUBN0
      CHARACTER*4 IRELAT
      CHARACTER*4 IRHSTG
C
C-------------------------------------------------------------------
C
      DIMENSION Y(*)
      DIMENSION X(*)
      DIMENSION TEMP1(*)
      DIMENSION TEMP2(*)
      DIMENSION TEMP3(*)
      DOUBLE PRECISION DTEMP1(*)
C
      DOUBLE PRECISION TOL
      DOUBLE PRECISION XPAR(2)
      DOUBLE PRECISION FVEC(2)
C
      EXTERNAL PAPFUN
      DOUBLE PRECISION XBAR
      COMMON/PAPCOM/MAXROW,NTOT2,XBAR
C
C-------------------------------------------------------------------
C
      CHARACTER*4 IFEEDB
      CHARACTER*4 IPRINT
      CHARACTER*240 ICOUT
C
      COMMON /MACH/IRD,IPR,CPUMIN,CPUMAX,NUMBPC,NUMCPW,NUMBPW
      COMMON /PRINT/IFEEDB,IPRINT
      COMMON /TEXTOU/ICOUT,NCOUT,ILOUT
C
C-----START POINT---------------------------------------------------
C
      ISUBN1='DPML'
      ISUBN2='AE  '
C
      IERROR='NO'
C
      IF(IBUGA3.EQ.'ON'.OR.ISUBRO.EQ.'MLAE')THEN
        WRITE(ICOUT,999)
  999   FORMAT(1X)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,51)
   51   FORMAT('**** AT THE BEGINNING OF DPMLAE--')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,52)IBUGA3
   52   FORMAT('IBUGA3 = ',A4)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,55)N,NVAR
   55   FORMAT('N,NVAR = ',2I8)
        CALL DPWRST('XXX','WRIT')
        DO56I=1,MIN(N,100)
          WRITE(ICOUT,57)I,Y(I),X(I)
   57     FORMAT('I,Y(I),X(I) = ',I8,2G15.7)
          CALL DPWRST('XXX','WRIT')
   56   CONTINUE
      ENDIF
C
C               ********************************************
C               **  STEP 11--                             **
C               **  CHECK THE INPUT ARGUMENTS FOR ERRORS  **
C               ********************************************
C
      ISTEPN='11'
      IF(IBUGA3.EQ.'ON' .OR. ISUBRO.EQ.'MLAE')
     1CALL TRACE2(ISTEPN,ISUBN1,ISUBN2)
C
      IF(N.LE.2)THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,1111)
 1111   FORMAT('***** ERROR IN POLYA-AEPPLI ',
     1         'MAXIMUM LIKELIHOOD--')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,1113)
 1113   FORMAT('      THE NUMBER OF OBSERVATIONS ',
     1         'FOR VARIABLE 1 IS LESS THAN 3')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,1115)N
 1115   FORMAT('      SAMPLE SIZE = ',I8)
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
      ENDIF
C
      IF(NVAR.EQ.1)THEN
        HOLD=Y(1)
        DO1135I=2,N
          IF(Y(I).NE.HOLD)GOTO1139
 1135   CONTINUE
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,1131)
 1131   FORMAT('***** ERROR IN POLYA-AEPPLI ',
     1         'MAXIMUM LIKELIHOOD--')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,1133)HOLD
 1133   FORMAT('      HAS ALL ELEMENTS = ',E15.7)
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
 1139   CONTINUE
C
        DO1145I=1,N
          IF(Y(I).LT.0.0)THEN
            WRITE(ICOUT,999)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,1147)
 1147       FORMAT('***** ERROR IN POLYA-AEPPLI MAXIMUM ',
     1             'LIKELIHOOD ESTIMATION')
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,1148)I,Y(I)
 1148       FORMAT('      ROW ',I8,' IS NON-POSITIVE (VALUE = ',
     1             G15.7,')')
            CALL DPWRST('XXX','WRIT')
            IERROR='YES'
            GOTO9000
          ENDIF
 1145   CONTINUE
C
        CALL SORT(Y,N,TEMP2)
        DO1160I=1,N
          Y(I)=TEMP2(I)
 1160   CONTINUE
        XMIN=Y(1)
        XMAX=Y(N)
C
        IRELAT='OFF'
        IRHSTG='OFF'
        XSTART=XMIN-0.5
        XSTOP=XMAX+0.5
        CLWID=1.0
        CALL DPBINI(Y,N,IRELAT,CLWID,XSTART,XSTOP,IRHSTG,
     1              TEMP2,TEMP1,N2,IBUGA3,IERROR)
        ICNT=0
        DO101I=1,N2
          IF(TEMP2(I).GT.0)THEN
            ICNT=ICNT+1
            TEMP2(ICNT)=TEMP2(I)
            TEMP1(ICNT)=TEMP1(I)
          ENDIF
 101    CONTINUE
        N2=ICNT
        IF(IERROR.EQ.'YES')GOTO9000
        F1=TEMP2(1)
        F2=TEMP2(2)
C
      ELSEIF(NVAR.EQ.2)THEN
        CALL SORTC(X,Y,N,TEMP1,TEMP2)
        NTOT=0
        DO1210I=1,N
          X(I)=TEMP1(I)
          Y(I)=TEMP2(I)
          NTOT=NTOT + Y(I)
 1210   CONTINUE
        F1=Y(1)
        F2=Y(2)
C
        DO1220I=1,N
          IF(Y(I).LT.0.0)THEN
            WRITE(ICOUT,999)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,1221)
 1221       FORMAT('***** ERROR IN POLYA-AEPPLI ',
     1             'MAXIMUM LIKELIHOOD--')
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,1223)
 1223       FORMAT('      A NEGATIVE FREQUENCY WAS SPECIFIED.')
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,1225)I,Y(I)
 1225       FORMAT('      ROW ',I8,' (AFTER SORTING) HAS FREQUENCY ',
     1             G15.7)
            CALL DPWRST('XXX','WRIT')
          ENDIF
 1220   CONTINUE
      ENDIF
C
      IF(IBUGA3.EQ.'ON'.OR.ISUBRO.EQ.'MLAE')THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,1301)N
 1301   FORMAT('AFTER SORT, N = ',I8)
        CALL DPWRST('XXX','WRIT')
        DO1310I=1,MAX(N,100)
          WRITE(ICOUT,1311)I,X(I),Y(I)
 1311     FORMAT('I,X(I),Y(I) = ',I8,2G15.7)
          CALL DPWRST('XXX','WRIT')
 1310   CONTINUE
      ENDIF
C
C               *********************************************
C               **  STEP 21--                              **
C               **  CARRY OUT CALCULATIONS                 **
C               **  FOR POLYA-AEPPLI MLE ESTIMATION        **
C               *********************************************
C
 2100 CONTINUE
C
      ISTEPN='21'
      IF(IBUGA3.EQ.'ON' .OR. ISUBRO.EQ.'MLAE')
     1CALL TRACE2(ISTEPN,ISUBN1,ISUBN2)
C
      IERROR='NO'
      IWRITE='OFF'
C
      IF(NVAR.EQ.1)THEN
        NTOT=N
        DO2105I=1,N
          ITEMP=INT(Y(I)+0.5)
          Y(I)=REAL(ITEMP)
 2105   CONTINUE
        CALL MEAN(Y,N,IWRITE,AMEAN,IBUGA3,IERROR)
        CALL SD(Y,N,IWRITE,ASD,IBUGA3,IERROR)
        AVAR=ASD**2
        AMIN=Y(1)
        AMAX=Y(N)
C
        IRELAT='OFF'
        IRHSTG='OFF'
        XSTART=AMIN-0.5
        XSTOP=AMAX+0.5
        CLWID=1.0
        CALL DPBINI(Y,N,IRELAT,CLWID,XSTART,XSTOP,IRHSTG,
     1              TEMP1,TEMP2,N2,IBUGA3,IERROR)
        IINDX=MAXNXT/2
        IF(N2.LE.IINDX)THEN
          IML=0
          DO2110I=1,N2
            TEMP3(I)=TEMP1(I)
            TEMP3(IINDX+I)=TEMP2(I)
 2110     CONTINUE
          IK=N2
        ELSE
          IML=1
        ENDIF
C
      ELSE
        AMIN=X(1)
        AMAX=X(N)
        CALL WEMEAN(X,Y,N,IWRITE,AMEAN,IBUGA3,IERROR)
        CALL WESD(X,Y,N,IWRITE,ASD,IBUGA3,IERROR)
        AVAR=ASD**2
        IINDX=MAXNXT/2
        IF(N.LE.IINDX)THEN
          IML=0
          DO2210I=1,N
            NTOT=NTOT+Y(I)
            TEMP3(I)=Y(I)
            TEMP3(IINDX+I)=X(I)
 2210     CONTINUE
          IK=N
        ELSE
          IML=1
        ENDIF
      ENDIF
C
      THETMO=2.0*AMEAN**2/(AVAR+AMEAN)
      PMOM=(AVAR - AMEAN)/(AVAR + AMEAN)
      THETFR=-LOG(F1/REAL(NTOT))
      PFR=1.0 - THETFR/AMEAN
      THETF2=-LOG(F1/REAL(NTOT))
      PF2=-F2/(F1*LOG(F1/REAL(NTOT)))
      THETML=THETMO
      PML=PMOM
C
      IF(IML.EQ.0)THEN
        IOPT=2
        TOL=1.0D-5
        NPAR=2
        NPRINT=-1
        INFO=0
        LWA=MAXNXT
        MAXROW=MAXNXT
        NTOT2=NTOT
C
        XBAR=DBLE(AMEAN)
        XPAR(1)=DBLE(THETML)
        XPAR(2)=DBLE(PML)
        CALL DNSQE(PAPFUN,JAC,IOPT,NPAR,XPAR,FVEC,TOL,NPRINT,INFO,
     1             DTEMP1,LWA,TEMP3,IK)
C
        THETML=REAL(XPAR(1))
        PML=REAL(XPAR(2))
      ENDIF
C
C               ***********************************************
C               **   STEP 42--                               **
C               **   WRITE OUT EVERYTHING                    **
C               **   FOR POLYA-AEPPLI MLE ESTIMATION     **
C               ***********************************************
C
      ISTEPN='42'
      IF(IBUGA3.EQ.'ON' .OR. ISUBRO.EQ.'MLNM')
     1CALL TRACE2(ISTEPN,ISUBN1,ISUBN2)
C
      IF(IPRINT.EQ.'ON')THEN
      IF(ICAPSW.EQ.'ON' .AND. ICAPTY.EQ.'HTML')THEN
C
C  STEP 1: END ASIS MODE AND WRITE A HEADER
C
 5001   FORMAT('
') 5011 FORMAT('') 5099 FORMAT('
')
        WRITE(ICOUT,5091)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5093)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5099)
        CALL DPWRST('XXX','WRIT')
C
      ELSEIF(ICAPSW.EQ.'ON' .AND. ICAPTY.EQ.'LATE')THEN
C
C  STEP 1: END ASIS FORMAT, START TABLE ENVIRONMENT, WRITE A HEADER,
C          AND WRITE A TABLE CAPTION
C
 8001 FORMAT(A1,'end{verbatim}')
 8003 FORMAT(A1,'begin{table}')
 8007 FORMAT(5X,'$ ',A1,1X,A1,' $ ',A1,A1,' ')
 8009 FORMAT(A1,'begin{center}')
 8011 FORMAT(5X,'{',A1,'bf Polya-Aeppli ',
     1       'Parameter Estimation}')
 8013 FORMAT(A1,'end{center}')
 8015 FORMAT(5X,'} ',A1,A1)
C
        CALL DPCONA(92,IBASLC)
C
        WRITE(ICOUT,8001)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8003)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8009)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8011)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8007)IBASLC,IBASLC,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8007)IBASLC,IBASLC,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8013)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
C
C STEP 2: START TABULAR ENVIRONMENT, WRITE SUBSESQUENT ROWS, END
C         TABULAR ENVIRONMENT
C
 8020 FORMAT(5X,A1,'begin{tabular} {lr}')
 8021 FORMAT(5X,'Number of Observations: & ',I8,2X,A1,A1)
 8022 FORMAT(5X,'Sample Mean: & ',G15.7,2X,A1,A1)
 8023 FORMAT(5X,'Sample Standard Deviation: & ',G15.7,2X,A1,A1)
 8024 FORMAT(5X,'Sample Minimum: & ',G15.7,2X,A1,A1)
 8025 FORMAT(5X,'Sample Maximum: & ',G15.7,2X,A1,A1)
 8026 FORMAT(5X,'First Frequency: & ',G15.7,2X,A1,A1)
 8027 FORMAT(5X,'Second Frequency: & ',G15.7,2X,A1,A1)
 8029 FORMAT(5X,'Estimate of Theta: & ',
     1       G15.7,2X,A1,A1)
 8030 FORMAT(5X,'Estimate of Lambda: & ',
     1       G15.7,2X,A1,A1)
 8031 FORMAT(5X,'Method of Moments: & ',2X,A1,A1)
 8032 FORMAT(5X,'Method of Zero Frequency and Mean: & ',
     1       2X,A1,A1)
 8033 FORMAT(5X,'Method of First Two Frequencies: & ',
     1       2X,A1,A1)
 8034 FORMAT(5X,'Method of Maximum Likelihood: & ',
     1       2X,A1,A1)
 8039 FORMAT(5X,' & ',2X,A1,A1)
 8040 FORMAT(5X,A1,'hline')
 8049 FORMAT(A1,'end{tabular}')
        WRITE(ICOUT,8009)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8020)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8021)N,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8022)AMEAN,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8023)ASD,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8024)AMIN,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8025)AMAX,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8026)F1,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8027)F2,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8039)IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,8031)IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8029)THETMO,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8030)PMOM,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,8039)IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8032)IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8029)THETFR,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8030)PFR,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,8039)IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8033)IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8029)THETF2,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8030)PF2,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,8039)IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8034)IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8029)THETML,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8030)PFML,IBASLC,IBASLC
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,8049)IBASLC
        CALL DPWRST('XXX','WRIT')
C
C STEP 3: END TABLE ENVIRONMENT, RESET ASIS MODE
C
 8091 FORMAT(A1,'end{table}')
 8093 FORMAT(A1,'end{center}')
 8099 FORMAT(A1,'begin{verbatim}')
        WRITE(ICOUT,8093)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8091)IBASLC
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,8099)IBASLC
        CALL DPWRST('XXX','WRIT')
      ELSEIF(ICAPSW.EQ.'ON' .AND. ICAPTY.EQ.'RTF')THEN
C
C
      ELSE
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4211)
 4211   FORMAT(10X,
     1  'POLYA-AEPPLI MAXIMUM LIKELIHOOD ESTIMATION:')
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4221)N
 4221   FORMAT('NUMBER OF OBSERVATIONS                   = ',I8)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4222)AMEAN
 4222   FORMAT('SAMPLE MEAN                              = ',G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4223)ASD  
 4223   FORMAT('SAMPLE STANDARD DEVIATION                = ',G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4225)AMIN  
 4225   FORMAT('SAMPLE MINIMUM                           = ',G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4227)AMAX  
 4227   FORMAT('SAMPLE MAXIMUM                           = ',G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4228)F1
 4228   FORMAT('FIRST FREQUENCY                          = ',G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4229)F2
 4229   FORMAT('SECOND FREQUENCY                         = ',G15.7)
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4231)
 4231   FORMAT('METHOD OF MOMENTS:')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4235)THETMO
 4235   FORMAT('ESTIMATE OF THETA:                       = ',G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4237)PMOM
 4237   FORMAT('ESTIMATE OF P                            = ',G15.7)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,4241)
 4241   FORMAT('METHOD OF ZERO FREQUENCY AND MEAN:')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4235)THETFR
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4237)PFR
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,4251)
 4251   FORMAT('METHOD OF FIRST TWO FREQUENCIES:')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4235)THETF2
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4237)PF2
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,4261)
 4261   FORMAT('METHOD OF MAXIMUM LIKELIHOOD:')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4235)THETML
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,4237)PML
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,4291)
 4291   FORMAT('ESTIMATES ARE SAVED IN THE INTERNAL PARAMETERS')
        CALL DPWRST('XXX','BUG ')
        WRITE(ICOUT,4292)
 4292   FORMAT('THETAMOM, PMOM, THETAFR, PFR, THETAF2, PF2 ',
     1         'THETAML, AND PML.')
        CALL DPWRST('XXX','BUG ')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','BUG ')
C
      ENDIF
      ENDIF
      GOTO9000
C
C               *****************
C               **  STEP 90--  **
C               **  EXIT       **
C               *****************
C
 9000 CONTINUE
      IF(IBUGA3.EQ.'ON'.OR.ISUBRO.EQ.'MLAE')THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,9011)
 9011   FORMAT('***** AT THE END       OF DPMLAE--')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,9012)N,IBUGA3,IERROR
 9012   FORMAT('N,IBUGA3,IERROR = ',I8,2X,A4,2X,A4)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,9015)N
 9015   FORMAT('N = ',I8)
        CALL DPWRST('XXX','WRIT')
      ENDIF
C
      RETURN
      END
      SUBROUTINE DPMLBE(Y,N,
     1XTEMP,MAXNXT,AUSER,BUSER,
     1A,B,ALPHA2,BETA2,ALPHA1,BETA1,
     1ALPHSE,BETASE,COVSE,
     1QP,XQPHAT,XQPLCL,XQPUCL,XQPSE,NPERC,
     1IOUNI1,IOUNI2,ALPHAP,
     1ICAPSW,ICAPTY,DTEMP1,
     1ISUBRO,IBUGA3,IERROR)
C
C     PURPOSE--THIS ROUTINE COMPUTES THE METHOD OF MOMENT AND
C              MAXIMUM LIKELIHOOD ESTIMATES FOR THE LOWER AND UPPER
C              LIMITS OF THE BETA DISTRIBUTION
C     EXAMPLE--BETA MOMENTS Y
C     REFERENCE--EVANS, HASTINGS, AND PEACOCK.  "STATISTICAL
C                DISTRIBUTIONS", THIRD EDITION, WILEY, 2000,
C                PP. 34-42.
C                JOHNSON, KOTZ, AND BALAKRISHNAN.  "CONTINUOUS
C                UNIVARIATE DISTRIBUTIONS, VOLUME II", SECOND
C                EDITION, WILEY, 1994.
C              --KARL BURY, "STATISTICAL DISTRIBUTIONS IN
C                ENGINEERING", CAMBRIDGE UNIVERSITY PRESS,
C                1999, CHAPTER 14.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2855
C     NOTE--DATAPLOT IS A REGISTERED TRADEMARK
C           OF THE NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY.
C           THIS SUBROUTINE MAY NOT BE COPIED, EXTRACTED,
C           MODIFIED, OR OTHERWISE USED IN A CONTEXT
C           OUTSIDE OF THE DATAPLOT LANGUAGE/SYSTEM.
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--2003/11
C     ORIGINAL VERSION--NOVEMBER  2003.
C     UPDATED         --DECEMBER  2004. CONFIDENCE INTERVALS FOR
C                                       SHAPE PARAMETERS
C     UPDATED         --JULY      2005. SOME COSMETIC CHANGES TO THE
C                                       OUTPUT
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
      CHARACTER*4 ICAPSW
      CHARACTER*4 ICAPTY
      CHARACTER*4 ISUBRO
      CHARACTER*4 IBUGA3
      CHARACTER*4 IERROR
C
      CHARACTER*4 IWRITE
C
      CHARACTER*4 ISUBN1
      CHARACTER*4 ISUBN2
      CHARACTER*4 ISTEPN
C
      CHARACTER*1 IBASLC
C
C---------------------------------------------------------------------
C
      DIMENSION Y(*)
      DIMENSION XTEMP(*)
      DOUBLE PRECISION DTEMP1(*)
C
      PARAMETER (NUMALP=6)
      DIMENSION ALPHA(NUMALP)
      DIMENSION ALOWAL(NUMALP)
      DIMENSION AUPPAL(NUMALP)
      DIMENSION ALOWBE(NUMALP)
      DIMENSION AUPPBE(NUMALP)
      DIMENSION ALOWA2(NUMALP)
      DIMENSION AUPPA2(NUMALP)
      DIMENSION ALOWB2(NUMALP)
      DIMENSION AUPPB2(NUMALP)
C
      DIMENSION QP(*)
      DIMENSION XQPHAT(*)
      DIMENSION XQPSE(*)
      DIMENSION XQPLCL(*)
      DIMENSION XQPUCL(*)
C
      DOUBLE PRECISION TOL
      DOUBLE PRECISION XPAR(2)
      DOUBLE PRECISION FVEC(2)
      DOUBLE PRECISION DAE
      DOUBLE PRECISION DRE
      DOUBLE PRECISION DXSTRT
      DOUBLE PRECISION DXLOW
      DOUBLE PRECISION DXUP
C
      EXTERNAL BETFUN
      DOUBLE PRECISION BETFU2
      EXTERNAL BETFU2
      DOUBLE PRECISION BETFU5
      EXTERNAL BETFU5
      EXTERNAL BETFU7
      EXTERNAL BETFU8
C
      DOUBLE PRECISION DANS(10)
      DOUBLE PRECISION DA
      DOUBLE PRECISION DB
      DOUBLE PRECISION DC
      DOUBLE PRECISION DALPHA
      DOUBLE PRECISION DBETA
      DOUBLE PRECISION DALPBE
      DOUBLE PRECISION DTERM1
      DOUBLE PRECISION DTERM2
      DOUBLE PRECISION DTERM3
      DOUBLE PRECISION DTERM4
C
C---------------------------------------------------------------------
C
      COMMON /BETAML/ BETALL, BETAUL
C
      INTEGER N2
      DOUBLE PRECISION DSUM3
      DOUBLE PRECISION DSUM4
      DOUBLE PRECISION DLLAB
      DOUBLE PRECISION DK
      COMMON/BETCO9/DSUM3,DSUM4,DLLAB,DK,N2
C
      DOUBLE PRECISION DBETA2
      COMMON/BETCO2/DBETA2
C
      DOUBLE PRECISION DALPH2
      COMMON/BETCO5/DALPH2
C
      DOUBLE PRECISION DBETA3
      COMMON/BETCO3/DBETA3
C
      DOUBLE PRECISION DALPH3
      COMMON/BETCO4/DALPH3
C
      COMMON/BETCO7/P7,BETA3
      COMMON/BETCO8/P8,ALPHA3
C
      DOUBLE PRECISION DN
C
      CHARACTER*4 IFEEDB
      CHARACTER*4 IPRINT
      CHARACTER*240 ICOUT
C
      COMMON /MACH/IRD,IPR,CPUMIN,CPUMAX,NUMBPC,NUMCPW,NUMBPW
      COMMON /PRINT/IFEEDB,IPRINT
      COMMON /TEXTOU/ICOUT,NCOUT,ILOUT
C
      DATA ALPHA /0.50, 0.25, 0.10, 0.05, 0.01, 0.001/
C
C-----START POINT-----------------------------------------------------
C
      ISUBN1='DPML'
      ISUBN2='BE  '
C
      IERROR='NO'
C
      IF(IBUGA3.EQ.'OFF'.AND.ISUBRO.NE.'MLBE')GOTO90
      WRITE(ICOUT,999)
  999 FORMAT(1X)
      CALL DPWRST('XXX','WRIT')
      WRITE(ICOUT,51)
   51 FORMAT('**** AT THE BEGINNING OF DPMLBE--')
      CALL DPWRST('XXX','WRIT')
      WRITE(ICOUT,52)IBUGA3
   52 FORMAT('IBUGA3 = ',A4)
      CALL DPWRST('XXX','WRIT')
      WRITE(ICOUT,55)N
   55 FORMAT('N = ',I8)
      CALL DPWRST('XXX','WRIT')
      DO56I=1,MIN(N,100)
      WRITE(ICOUT,57)I,Y(I)
   57 FORMAT('I,Y(I) = ',I8,E15.7)
      CALL DPWRST('XXX','WRIT')
   56 CONTINUE
   66 CONTINUE
   90 CONTINUE
C
C               ********************************************
C               **  STEP 11--                             **
C               **  CHECK THE INPUT ARGUMENTS FOR ERRORS  **
C               ********************************************
C
      ISTEPN='11'
      IF(IBUGA3.EQ.'ON')CALL TRACE2(ISTEPN,ISUBN1,ISUBN2)
C
      IF(N.LE.2)THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,1111)
 1111   FORMAT('***** ERROR IN DPMLBE--THE NUMBER OF OBSERVATIONS ',
     1        'FOR VARIABLE 1 IS LESS THAN 3')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,1112)N
 1112   FORMAT('SAMPLE SIZE = ',I8)
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
      ENDIF
C
      HOLD=Y(1)
      DO1135I=2,N
      IF(Y(I).NE.HOLD)GOTO1139
 1135 CONTINUE
 1130 CONTINUE
      WRITE(ICOUT,999)
      CALL DPWRST('XXX','WRIT')
      WRITE(ICOUT,1131)HOLD
 1131 FORMAT('***** NOTE FROM DPMLBE--VARIABLE 1 ',
     1'HAS ALL ELEMENTS = ',E15.7)
      CALL DPWRST('XXX','WRIT')
      GOTO9000
 1139 CONTINUE
C
      IF(NPERC.GT.0)THEN
        DO1145I=1,NPERC
          IF(QP(I).LE.0.0 .OR. QP(I).GE.100.0)THEN
            WRITE(ICOUT,999)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,1141)
 1141       FORMAT('***** WARNING IN BETA MAXIMUM LIKELIHOOD--')
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,1143)QP(I)
 1143       FORMAT('      REQUESTED PERCENTILE (',G15.7,') IS ',
     1             'OUTSIDE THE (0,100) INTERVAL')
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,1144)
 1144       FORMAT('      NO PERCENTILE CONFIDENCE LIMITS WILL BE ',
     1             'COMPUTED.')
            CALL DPWRST('XXX','WRIT')
            NPERC=0
          ENDIF
 1145   CONTINUE
      ENDIF
C
 1290 CONTINUE
C
C               *************************************
C               **  STEP 31--                      **
C               **  CARRY OUT CALCULATIONS         **
C               **  FOR BETA MOMENT/MLE ESTIMATION **
C               *************************************
C
 3100 CONTINUE
C
      ISTEPN='41'
      IF(IBUGA3.EQ.'ON')CALL TRACE2(ISTEPN,ISUBN1,ISUBN2)
C
      IERROR='NO'
      IWRITE='OFF'
      CALL MINIM(Y,N,IWRITE,XMIN,IBUGA3,IERROR)
      CALL MAXIM(Y,N,IWRITE,XMAX,IBUGA3,IERROR)
      CALL MEAN(Y,N,IWRITE,XMEAN,IBUGA3,IERROR)
      CALL VAR(Y,N,IWRITE,XVAR,IBUGA3,IERROR)
      XSD=SQRT(XVAR)
C
CCCCC ALLOW FOR USER SPECIFIED LOWER AND UPPER LIMITS
C
      IF((AUSER.EQ.CPUMIN .OR. BUSER.EQ.CPUMIN) .OR.
     1   (AUSER.GE.XMIN .OR. BUSER.LE.XMAX))THEN
        IF((XMIN.GE.0.0 .AND. XMIN.LE.1.0) .AND.
     1     (XMAX.GE.0.0 .AND. XMAX.LE.1.0))THEN
          A=0.0
          B=1.0
        ELSE
          A=XMIN - 1.0E-12
          B=XMAX + 1.0E+12
        ENDIF
        BETALL=A
        BETAUL=B
      ELSE
        BETALL=AUSER
        BETAUL=BUSER
        A=AUSER
        B=BUSER
      ENDIF
C
      XMEAN1=(XMEAN-A)/(B-A)
      VAR1=XVAR/((B-A)**2)
      ALPHA1=XMEAN1*(XMEAN1*(1.0-XMEAN1)/VAR1 - 1.0)
      BETA1=(1.0-XMEAN1)*(XMEAN1*(1.0-XMEAN1)/VAR1 - 1.0)
C
      XPAR(1)=DBLE(ALPHA1)
      XPAR(2)=DBLE(BETA1)
      DPROD1=1.0D0
      DPROD2=1.0D0
      DN=DBLE(N)
C
      DO3101I=1,N
        DTERM1=DBLE((B-Y(I))/(B-A))**(1.0D0/DN)
        DTERM2=DBLE( (Y(I)-A)/(B-A))**(1.0D0/DN)
        IF(DTERM1.NE.0.0D0)DPROD1=DPROD1*DTERM1
        IF(DTERM2.NE.0.0D0)DPROD2=DPROD2*DTERM2
 3101 CONTINUE
      XPAR(1)=0.5D0*(1.0D0 - DPROD1)/(1.0D0 - DPROD2 - DPROD1)
      DO3103I=1,N
        DTERM1=DBLE((Y(I)-A)/(B-A))**(1.0D0/DN)
        DTERM2=DBLE( (B-Y(I))/(B-A))**(1.0D0/DN)
        IF(DTERM1.NE.0.0D0)DPROD1=DPROD1*DTERM1
        IF(DTERM2.NE.0.0D0)DPROD2=DPROD2*DTERM2
 3103 CONTINUE
      XPAR(2)=0.5D0*(1.0D0 - DPROD1)/(1.0D0 - DPROD1 - DPROD2)
C
      IOPT=2
      TOL=1.0D-6
      NVAR=2
      NPRINT=-1
      INFO=0
      LWA=MAXNXT
      CALL DNSQE(BETFUN,JAC,IOPT,NVAR,XPAR,FVEC,TOL,NPRINT,INFO,
     1           DTEMP1,MAXNXT,Y,N)
C
      ALPHA2=REAL(XPAR(1))
      BETA2=REAL(XPAR(2))
C
CCCCC CONFIDENCE INTERVALS FOR SHAPE PARAMETERS
C
      DN=DBLE(N)
      DALPHA=DBLE(ALPHA2)
      DBETA=DBLE(BETA2)
      DALPBE=DBLE(ALPHA2 + BETA2)
C
      KODE=1
      NTEMP=1
      M=1
      NZ=0
C
      CALL DPSIFN(DALPHA,NTEMP,KODE,M,DANS,NZ,IERR)
      DA=DANS(1)
      IF(IERR.EQ.1)THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3201)
 3201   FORMAT('***** ERROR FROM BETA MAXIMUM LIKELIHOOD--')
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3203)
 3203   FORMAT('      UNABLE TO COMPUTE TRIGAMMA FUNCTION.')
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
      ELSEIF(IERR.EQ.2)THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3201)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3205)
 3205   FORMAT('      OVERFLOW IN COMPUTING THE TRIGAMMA FUNCTION.')
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
      ELSEIF(IERR.EQ.3)THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3201)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3207)
 3207   FORMAT('      OVERFLOW IN COMPUTING THE TRIGAMMA FUNCTION.')
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
      ENDIF
C
      CALL DPSIFN(DBETA,NTEMP,KODE,M,DANS,NZ,IERR)
      DB=DANS(1)
      IF(IERR.EQ.1)THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3201)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3203)
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
      ELSEIF(IERR.EQ.2)THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3201)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3205)
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
      ELSEIF(IERR.EQ.3)THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3201)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3207)
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
      ENDIF
C
      CALL DPSIFN(DALPBE,NTEMP,KODE,M,DANS,NZ,IERR)
      DC=DANS(1)
      IF(IERR.EQ.1)THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3201)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3203)
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
      ELSEIF(IERR.EQ.2)THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3201)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3205)
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
      ELSEIF(IERR.EQ.3)THEN
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3201)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,3207)
        CALL DPWRST('XXX','WRIT')
        IERROR='YES'
        GOTO9000
      ENDIF
C
      DTERM1=1.0D0/(DN*(DA*DB - DC*(DA+DB)))
      DTERM2=DTERM1*(DB-DC)
      ALPHSE=REAL(DSQRT(DTERM2))
      DTERM2=DTERM1*(DA-DC)
      BETASE=REAL(DSQRT(DTERM2))
      DTERM2=DTERM1*DC
      COVSE=REAL(DSQRT(DTERM2))
C
      DO3310I=1,NUMALP
        ALP=ALPHA(I)
        P=1.0-(ALP/2.0)
        CALL NORPPF(P,PPF)
        ALOWAL(I)=ALPHA2 - PPF*ALPHSE
        AUPPAL(I)=ALPHA2 + PPF*ALPHSE
        ALOWBE(I)=BETA2 - PPF*BETASE
        AUPPBE(I)=BETA2 + PPF*BETASE
 3310 CONTINUE
C
      N2=N
      DA=DBLE(A)
      DB=DBLE(B)
      DALPH2=DBLE(ALPHA2)
      DALPH3=DBLE(ALPHA2)
      DBETA2=DBLE(BETA2)
      DBETA3=DBLE(BETA2)
      DSUM3=0.0D0
      DSUM4=0.0D0
      DO3320I=1,N
        DTEMP1(I)=DBLE(Y(I))
        DSUM3=DSUM3 + DLOG(DBLE(Y(I)) - DA)
        DSUM4=DSUM4 + DLOG(DB - DBLE(Y(I)))
 3320 CONTINUE
      DSUM3=DSUM3/(DN*(DB - DA))
      DSUM4=DSUM4/(DN*(DB - DA))
C
      DTERM1=-DN*DLBETA(DALPH2,DBETA2)
      DTERM2=DN*(DALPH2-1.0D0)*DSUM3
      DTERM3=DN*(DBETA2-1.0D0)*DSUM4
      DLLAB=DTERM1 + DTERM2 + DTERM3
C
      DAE=1.D-7
      DRE=1.D-7
      NUTEMP=1
C
      DO3410I=1,NUMALP
        ALP=ALPHA(I)
        CALL CHSPPF(1.0-ALP,NUTEMP,APPF)
        DK=DBLE(APPF)
C
        DXSTRT=DBLE(ALOWAL(I))
        DXLOW=DXSTRT/5.0D0
        DXUP=DBLE(ALPHA2)
        CALL DFZER2(BETFU2,DXLOW,DXUP,DXSTRT,DRE,DAE,IFLAG,DTEMP1)
        ALOWA2(I)=REAL(DXLOW)
C
        DXSTRT=DBLE(AUPPAL(I))
        DXUP=DXSTRT*5.0D0
        DXLOW=DBLE(ALPHA2)
        CALL DFZER2(BETFU2,DXLOW,DXUP,DXSTRT,DRE,DAE,IFLAG,DTEMP1)
        AUPPA2(I)=REAL(DXLOW)
C
        DXSTRT=DBLE(ALOWBE(I))
        DXLOW=DXSTRT/5.0D0
        DXUP=DBLE(BETA2)
        CALL DFZER2(BETFU5,DXLOW,DXUP,DXSTRT,DRE,DAE,IFLAG,DTEMP1)
        ALOWB2(I)=REAL(DXLOW)
C
        DXSTRT=DBLE(AUPPBE(I))
        DXUP=DXSTRT*5.0D0
        DXLOW=DBLE(BETA2)
        CALL DFZER2(BETFU5,DXLOW,DXUP,DXSTRT,DRE,DAE,IFLAG,DTEMP1)
        AUPPB2(I)=REAL(DXLOW)
C
 3410 CONTINUE
C
C  CONFIDENCE INTERVALS FOR SELECTED PERCENTILES.
C
C  1. STANDARD ERROR USES TECHNIQUE DEMONSTRATED IN EXAMPLE 14.3
C     (PP. 256-257) OF BURY.  THIS IS BASED ON PROPOGATION OF ERROR.
C 
C  2. CONFIDENCE INTERVAL IS THEN GENERATED USING NORMAL
C     APPROXIMATION (EXAMPLE 14.3 OF BURY).
C
      IF(NPERC.GE.1)THEN
C
        ALPHL=ALPHAP/2.0
        ALPHU=1.0 - ALPHAP/2.0
        CALL NORPPF(ALPHU,Z95)
C
        ALPHA3=ALPHA2
        BETA3=BETA2
        IORD=1
        EPS=0.001
        ACCUR=0.0
C
        WRITE(IOUNI1,3531)
        WRITE(IOUNI1,3532)
        DO3529I=1,NPERC
          QPTEMP=QP(I)/100.0
          CALL BETPPF(QPTEMP,ALPHA2,BETA2,APPF)
          XQPHAT(I)=APPF
C
          P7=QPTEMP
          P8=QPTEMP
C
          IFAIL=0
C
          ALPHAT = ALPHA2
          ALPHMN = 0.0001
          ALPHMX = ALPHA2 + 10.0
          CALL DIFF(IORD,ALPHAT,ALPHMN,ALPHMX,BETFU7,EPS,ACCUR,
     1              D1,ERROR,IFAIL)
C
          IF(IFAIL.EQ.1)THEN
            WRITE(ICOUT,999)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3501)
 3501       FORMAT('***** WARNING IN NUMERICAL DERIVATIVE FOR BETA ',
     1             'MAXIMUM LIKELIHOOD PERCENTILES.')
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3503)
 3503       FORMAT('      THE ESTIMATED ERROR IN THE RESULT ',
     1             'EXCEEDS THE')
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3505)
 3505       FORMAT('      REQUESTED ERROR, BUT THE MOST ACCURATE ',
     1             'RESULT')
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3507)
 3507       FORMAT('      POSSIBLE HAS BEEN RETURNED.')
            CALL DPWRST('XXX','BUG ')
          ELSEIF(IFAIL.EQ.2)THEN
            WRITE(ICOUT,999)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3511)
 3511       FORMAT('***** ERROR IN NUMERICAL DERIVATIVE FOR BETA ',
     1             'MAXIMUM LIKELIHOOD PERCENTILES.')
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3513)
 3513       FORMAT('      ERROR IN THE INPUT TO THE DIFF ROUTINE.')
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3515)
 3515       FORMAT('      NO PERCENTILES WILL BE GENERATED.')
            CALL DPWRST('XXX','BUG ')
            NPERC=0
          ELSEIF(IFAIL.EQ.3)THEN
            WRITE(ICOUT,999)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3511)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3523)
 3523       FORMAT('      THE INTERVAL FOR DIFFERENTIATION, (',G15.7,
     1             ',',G15.7,')')
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3525)
 3525       FORMAT('      IS TOO SMALL.')
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3515)
            CALL DPWRST('XXX','BUG ')
            D1=0.0
            NPERC=0
          ENDIF
C
          BETAT = BETA2
          BETAMN = 0.0001
          BETAMX = BETA2 + 10.0
          CALL DIFF(IORD,BETAT,BETAMN,BETAMX,BETFU8,EPS,ACCUR,
     1              D2,ERROR,IFAIL)
C
          IF(IFAIL.EQ.1)THEN
            WRITE(ICOUT,999)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3501)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3503)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3505)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3507)
            CALL DPWRST('XXX','BUG ')
          ELSEIF(IFAIL.EQ.2)THEN
            WRITE(ICOUT,999)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3511)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3513)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3515)
            CALL DPWRST('XXX','BUG ')
            NPERC=0
          ELSEIF(IFAIL.EQ.3)THEN
            WRITE(ICOUT,999)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3511)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3523)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3525)
            CALL DPWRST('XXX','BUG ')
            WRITE(ICOUT,3515)
            CALL DPWRST('XXX','BUG ')
            D2=0.0
            NPERC=0
          ENDIF
          V11=ALPHSE**2
          V22=BETASE**2
          V21=COVSE
          V12=V21
          TERM11=(D1*ALPHSE)**2
          TERM22=(D2*BETASE)**2
          TERM12=2.0*D2*D1*COVSE**2
          SEXQP=TERM11+TERM12+TERM22
          IF(SEXQP.GE.0.0)THEN
            SEXQP=SQRT(SEXQP)
          ELSE
            SEXQP=0.0
          ENDIF
          XQPSE(I)=SEXQP
          XQPLCL(I)=XQPHAT(I) - Z95*SEXQP
          XQPUCL(I)=XQPHAT(I) + Z95*SEXQP
          WRITE(IOUNI1,'(5E15.7)')
     1         QP(I),XQPHAT(I),XQPSE(I),XQPLCL(I),XQPUCL(I)
 3529   CONTINUE
 3531   FORMAT(15X,'       POINT     ','   STANDARD   ',
     1         '     LOWER     ',
     1         '     UPPER')
 3532   FORMAT('    PERCENTILE ','     ESTIMATE   ',
     1         '     ERRROR     ',
     1         'CONFIDENCE LIMIT ','CONFIDENCE LIMIT')
      ENDIF
C               *********************************
C               **   STEP 42--                 **
C               **   WRITE OUT EVERYTHING      **
C               **   FOR BETA MLE ESTIMATION   **
C               **********************************
C
      ISTEPN='42'
      IF(IBUGA3.EQ.'ON')CALL TRACE2(ISTEPN,ISUBN1,ISUBN2)
C
      IF(IPRINT.EQ.'ON')THEN
      IF(ICAPSW.EQ.'ON' .AND. ICAPTY.EQ.'HTML')THEN
C
C  STEP 1: END ASIS MODE AND WRITE A HEADER
C
 5001   FORMAT('
') 5002 FORMAT('

BETA DISTRIBUTION PARAMETER ESTIMATION

') 5004 FORMAT('

') WRITE(ICOUT,5001) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,999) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,5002) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,5004) CALL DPWRST('XXX','WRIT') WRITE(ICOUT,999) CALL DPWRST('XXX','WRIT') C C STEP 2: START TABLE AND DEFINE A CAPTION C 5011 FORMAT('') 5099 FORMAT('
')
        WRITE(ICOUT,5091)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5093)
        CALL DPWRST('XXX','WRIT')
C
C  STEP 2: START TABLE AND DEFINE A CAPTION
C
 5111   FORMAT('')
        WRITE(ICOUT,5191)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5193)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
C
C  STEP 2: START TABLE AND DEFINE A CAPTION
C
 5217   FORMAT('      Confidence Limits for the Beta ',
     1         'Parameter
')
        WRITE(ICOUT,5111)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5113)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5115)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5217)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5119)
        CALL DPWRST('XXX','WRIT')
C
C  STEP 3: DEFINE HEADER ROW
C
        WRITE(ICOUT,5121)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5123)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5134)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5127)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5136)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5137)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5127)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5136)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5138)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5127)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5139)
        CALL DPWRST('XXX','WRIT')
C
        WRITE(ICOUT,5121)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5123)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5125)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5127)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5129)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5131)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5127)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5129)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5133)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5127)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5129)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5131)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5127)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5129)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5133)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5127)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5139)
        CALL DPWRST('XXX','WRIT')
C
C  FOLLOWING ADDS A RULE LINE BETWEEN HEADER AND DATA LINES
C
        WRITE(ICOUT,5121)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5161)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5162)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5147)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5139)
        CALL DPWRST('XXX','WRIT')
C
C  STEP 4: DEFINE DATA ROW
C
        DO5240I=1,NUMALP
          ATEMP=100.0*(1.0 - ALPHA(I))
          WRITE(ICOUT,5141)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5149)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5151)ATEMP
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5147)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5149)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5151)ALOWBE(I)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5147)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5149)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5151)AUPPBE(I)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5147)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5149)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5151)ALOWB2(I)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5147)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5149)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5151)AUPPB2(I)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5147)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5159)
          CALL DPWRST('XXX','WRIT')
 5240   CONTINUE
C
C  STEP 4: END THE TABLE AND RESET ASIS MODE
C
        WRITE(ICOUT,5191)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,5193)
        CALL DPWRST('XXX','WRIT')
        WRITE(ICOUT,999)
        CALL DPWRST('XXX','WRIT')
C
C  START THE TABLE FOR PERCENTILE CONFIDENCE INTERVALS
C
        IF(NPERC.GT.0)THEN
          WRITE(ICOUT,5801)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,999)
          WRITE(ICOUT,5811)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5813)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5815)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5817)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5816)ALPHAP
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5819)
          CALL DPWRST('XXX','WRIT')
C
C  STEP 4: DEFINE DATA ROW
C
          WRITE(ICOUT,5841)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5849)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5861)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5848)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5849)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5862)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5848)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5849)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,55862)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5848)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5849)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5863)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5848)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5849)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5864)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5848)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5859)
          CALL DPWRST('XXX','WRIT')
C
          WRITE(ICOUT,5841)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5870)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5872)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5847)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5859)
          CALL DPWRST('XXX','WRIT')
C
          DO55880I=1,NPERC
            WRITE(ICOUT,5841)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5843)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5851)QP(I)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5847)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5843)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5851)XQPHAT(I)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5847)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5843)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5851)XQPSE(I)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5847)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5843)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5851)XQPLCL(I)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5847)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5843)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5851)XQPUCL(I)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5847)
            CALL DPWRST('XXX','WRIT')
            WRITE(ICOUT,5859)
            CALL DPWRST('XXX','WRIT')
55880     CONTINUE
C
C  END THE TABLE AND RESET ASIS MODE
C
          WRITE(ICOUT,5191)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,5193)
          CALL DPWRST('XXX','WRIT')
          WRITE(ICOUT,999)
          CALL DPWRST('XXX','WRIT')
        ENDIF
C
 5801     FORMAT('
') 5811 FORMAT('') C 5801 FORMAT('') 5811 FORMAT('