C*********************************************************************** C THIS SUBROUTINE PROVIDES THE FRAMEWORK FOR MAXIMUM LIKELIHOOD C ESTIMATION FOR GROUPED OR UNGROUPED PERPENDICULAR DISTANCES. IT C IS BASED ON NEWTONS METHOD WITH THE MODIFICATION OF USING A MOD- C IFIED MARQUARDT APPROACH TO ALLOW THE PARAMETERS TO BE CONSTRAINED. C IF THE PARAMETERS ATTEMPT TO GO OUT OF THE FEASIBLE REGION THEN C SUBROUTINE MARQUA IS CALLED TO FIND AN ACCEPTABLE STEP. MARQUA C IS ALSO CALLED IF NEWTONS METHOD DOES NOT FIND A STEP WHICH WILL C INCREASE THE VALUE OF THE LOG-LIKELIHOOD FUNCTION. SUBROUTINE C START IS USED TO GET STARTING VALUES FOR THE PROCEDURE; HOWEVER C USER SUPPLIED STARTING VALUES CAN ALSO BE USED(STORED IN PSTRT). C FOR ESTIMATORS WHICH HAVE A VARIABLE NUMBER OF PARAMETERS(I.E., C FOURIER SERIES) THE LOGICAL VARIABLE LOOP IS SET TO TRUE IN C WHICH CASE ADDITIONAL PARAMETERS ARE ADDED UNTIL THE LIKELIHOOD C RATIO TEST IS NOT SIGNIFICANT AT THE .05 LEVEL. IF THE C NUMBER OF PARAMETERS HAS BEEN SET (MSET) THEN NO LOOPING IS DONE. C C SUBROUTINES CALLED: C FOR MAXIMIZATION- START,DOMAIN,UPDATE,MARQUA,INVERT C FOR LIKELIHOOD EQUATIONS C GROUPED DATA - CELLS C UNGROUPED DATA - UMLE C FOR F(0) ESTIMATION - CALC C FOR OUTPUT - HEADER,NAME C*********************************************************************** SUBROUTINE MLE (IND) C*********************************************************************** C DECLARATIONS C*********************************************************************** INCLUDE 'PARMTR.INC' LOGICAL GRP, POOL, PEST, SEST, DESC, DEF, CUTP, TRUNC, WARN, LOOP, 1 OK, TRAN REAL LOGL DOUBLE PRECISION RF0, PRF0A DOUBLE PRECISION PCELL, CELL, SVCMAT(MAXPAR,MAXPAR), SPAR(MAXPAR), 1 SXLL, PRF0B DOUBLE PRECISION SCELL(MAXCEL) DOUBLE PRECISION OXLL, W(MAXPR3) INTEGER CNT, ERR(4), NERR CHARACTER*1 DASH CHARACTER*8 ERRSTR LOGICAL STRT, MSET, SET, MSTRT, HELP C*********************************************************************** C COMMON STATEMENTS C*********************************************************************** COMMON /KEEP/ SUM, SSQ COMMON /INTEG/ RF0, PRF0A, PRF0B COMMON /LPASS/ MSET, MSTRT COMMON /STORE/ D(13), DL(13), DU(13), STD(13), LOGL(10), PB(5,10) COMMON /PDOPT/ STRT, SET, PSTRT(10,MAXPAR), NPSET(10) DOUBLE PRECISION PAR, VCMAT, G, XLL COMMON /DPAR/ PAR(MAXPR2), VCMAT(MAXPAR,MAXPAR), G(MAXPAR), 1 XLL, NPAR, INDEX COMMON /INTER/ K, CUT(MAXCEL), F(MAXCEL), NCUT, NKC(5), 1 NK(5), RFREQ(5,MAXCEL), CCUT(5,MAXCEL) COMMON /PART/ PCELL(MAXCEL,MAXPAR), CELL(MAXCEL) COMMON /SOLN/ F0, VARF0, FMAX, FMIN, FZ(100) COMMON /OPTION/ GRP, POOL, PEST, SEST, DESC, DEF, CUTP, TRUNC, 1 HELP COMMON /NUM/ XL(MAXLIN), WIDTH, N, CNT, CONV(3), VARN, IDF, WARN COMMON /ERROR/ IER(4), LER DOUBLE PRECISION BNDL, BNDU COMMON /BOUNDY/ BNDL(MAXPAR), BNDU(MAXPAR) LOGICAL LER CHARACTER*60 ERRMSG(129:131) EXTERNAL XLIKE C*********************************************************************** C DATA STATEMENTS C*********************************************************************** DATA DASH /'-'/ DATA ERRMSG/ 1 'Initial Hessian matrix is not positive definite.', 2 'Rounding errors forced termination of iterations.', 3 'Number of function calls exceeds program limit.'/ C*********************************************************************** C INITIALIZE VARIABLES - THE VARIABLE SET IS TRUE IF STRT (STARTING C VALUES) IS TRUE. IF SET IS TRUE AND STRT IS FALSE THEN ONLY THE C NUMBER OF PARAMETERS HAS BEEN SET. SINCE SET AND STRT NEED NOT C BE TRUE FOR ALL ESTIMATORS MSET AND MSTRT ARE USED FOR THE PAR- C TICULAR ESTIMATOR BEING EXAMINED. IF THE NUMBER OF PARAMETERS C (NPSET(IND)) IS ZERO THEN FOR THIS ESTIMATOR NEITHER STARTING C VALUES NOR THE NUMBER OF PARAMETERS HAVE BEEN SET (I.E., BOTH C MSET AMD MSTRT ARE FALSE). C*********************************************************************** NP=MAXPAR DO 10 I=1,NP PAR(I)=0.0 DO 10 J=1,NP 10 VCMAT(I,J)=0.0 NPAR=1 CALL HEADER (1) MSET=SET MSTRT=STRT IF (NPSET(IND).EQ.0) THEN MSET=.FALSE. MSTRT=.FALSE. ENDIF 20 LER=.FALSE. LOOP=.FALSE. IF (INDEX.EQ.1) LOOP=.TRUE. CALL NAME (INDEX) WRITE (6,350) C C LOOP FOR NUMERICAL SOLUTION STARTS NOW C 30 CALL START (IND) C C SET STARTING VALUES TO ADMISSIBLE RANGE C TRAN=.FALSE. CALL DOMAIN (INDEX,PAR,PAR(MAXPAR+1),NPAR,TRAN) WRITE (6,410) (DASH,I=1,MIN(NPAR*14+24,79)) DO 60 I=1,4 60 IER(I)=0 C*********************************************************************** C CALCULATE MAXIMUM LIKELIHOOD ESTIMATES FOR EITHER GROUPED C OR UNGROUPED DATA C*********************************************************************** NSIG=5 IOPT=3 MAXFN=5000 CALL OPTMIZ(XLIKE,NPAR,NSIG,MAXFN,IOPT, 1 PAR(MAXPAR+1),VCMAT,G,OXLL,W,IER1) ITER=INT(W(2)) NSIG=INT(W(3)) WRITE (6,410) (DASH,I=1,MIN(NPAR*14+24,79)) IF (IER1.EQ.0) THEN WRITE (6,400) 'Convergence achieved',ITER ELSE WRITE (6,400) 'No convergence',ITER WRITE (6,'(1X,A)') ERRMSG(IER1) ENDIF WRITE (6,'(1X,A,I4)') 1 'Probable number of significant digits for estimates is ',NSIG C*********************************************************************** C CALCULATE LOG-LIKELIHOOD EQUATIONS AND THEIR PARTIALS AND THE LOG- C LIKELIHOOD VALUE FOR GROUPED AND UNGROUPED DATA. C THEN INVERT THE VAR-COV MATRIX FOR PARAMETER VARIANCES C*********************************************************************** TRAN=.TRUE. CALL DOMAIN(INDEX,PAR,PAR(MAXPAR+1),NPAR,TRAN) IF (GRP) THEN CALL CELLS (1) ELSE CALL UMLE (1) ENDIF CALL INVERT (VCMAT,NPAR,NP,ACC) IF (ACC.NE.1.0) IER(2)=-1 C C IF THE NUMBER OF PARAMETERS IN THE MODEL IS VARIABLE THEN THE C LIKELIHOOD VALUES,PARAMETER ESTIMATES AND THE LIKE MUST BE STORED C WHILE ANOTHER MODEL IS BEING EXAMINED WITH MORE PARAMETERS. ALSO C A LIKELIHOOD RATIO TEST MUST BE PERFORMED AND CHECKS MUST BE MADE C SO THAT THE NUMBER OF PARAMETERS DOES NOT EXCEED THE MAXIMUM NO. C POSSIBLE NOR THE NUMBER OF INTERVALS IF THE DATA ARE GROUPED. C 230 IF (MSET.OR.(.NOT.LOOP)) GO TO 330 IF (NPAR.NE.1) GO TO 240 WRITE (6,420) GO TO 250 240 RATIO=-2.*(SXLL-XLL) CALL ONE (RATIO,P) P=1.-P WRITE (6,430) RATIO,P IF(RATIO.GE.0) GO TO 245 WRITE(6,435) GO TO 305 245 IF (P.GT.0.05) GO TO 300 NP1=NPAR+1 WRITE (6,440) NP1 250 SXLL=XLL DO 260 I=1,K 260 SCELL(I)=CELL(I) DO 270 I=1,NPAR SPAR(I)=PAR(I) DO 270 J=1,NPAR 270 SVCMAT(I,J)=VCMAT(I,J) IF (.NOT.GRP) GO TO 280 IF (NPAR+1.LT.K) GO TO 280 WRITE (6,450) GO TO 330 280 IF (NPAR+1.LE.NP) GO TO 290 WRITE (6,460) NP GO TO 330 290 NPAR=NPAR+1 GO TO 30 300 WRITE (6,470) 305 XLL=SXLL DO 310 I=1,K 310 CELL(I)=SCELL(I) NPAR=NPAR-1 DO 320 I=1,NPAR PAR(I)=SPAR(I) DO 320 J=1,NPAR 320 VCMAT(I,J)=SVCMAT(I,J) C C CALL CALC TO GET AN ESTIMATE OF F(0) AND ITS VARIANCE. C 330 IF (LER) WRITE (6,490) CALL HEADER (1) CALL NAME (INDEX) CALL CALC LOGL(IND)=SNGL(XLL) RETURN C C FORMAT STATEMENTS C 350 FORMAT ('0Some results during the iteration for maximum', 1' likelihood estimates'/ 2' are given below:'/ 3' ITER = the iteration count'/ 2' XLL = the log-likelihood value at the parameter value(s)', 3' for that iteration'/ 4' ERROR = error conditions encountered during that iteration.'/ 5' The numbers given correspond to conditions explained below.'/ 6' Parameter values are printed in order, i.e., A(1),...,A(NPAR)'/ 7'0ITER XLL ERROR Parameter Values') 400 FORMAT (1X,A,', number of function evaluations was ',I4) 410 FORMAT (1X,79A1) 420 FORMAT (//'0Now a model with 2 parameters will be fit.') 430 FORMAT ('0Likelihood ratio test value =',g13.6/ 1' Probability of a greater value =',f7.4) 435 FORMAT('0The log-likelihood value is less than the previous', 1' estimator as a'/ 2' result of numerical problems. Therefore, the previous', 3' estimator will be used.') 440 FORMAT (//'0This ratio is significant at the .05 level, so an', 1' attempt to fit'/ 2' a model with ',I2,' parameters will be made.') 450 FORMAT ('0Further model fits must be discontinued, otherwise'/ 1' the number of parameters would equal the number of cells.') 460 FORMAT ('0Further model fits must be discontinued, otherwise', 1' the number'/ 2' of parameters would exceed the maximum of ',I2,' allowed.') 470 FORMAT ('0This ratio is not significant at the .05 level, so'/ 1' the previous model will be used to represent the data.') 490 FORMAT (//'0The following are error conditions which may be', 1' encountered during iteration:'/ 2' 1 - Numerical integration achieved less than 6 significant', 3' digits of accuracy,'/ 4' 2 - Parameter has approached the boundary of the search', 5' space,'/ 6' 3 - For grouped data, some excessively large or small', 7' probabilities for an'/5X,'interval of F(X) were encountered,'/ 8' 4 - For grouped data, the sum of the probabilities for the', 9' intervals of F(X)'/5X,'was not 1, so they were normalized.'// A' These error conditions may suggest that the estimator is', B' having problems'/ C' fitting the data; however, convergence can be achieved even', D' if problems'/ E' are encountered. If convergence is not achieved then try', F' the final parameter '/ G' estimates as starting values for another run. If no', H' acceptable step was '/ I' found, try various starting values.') END