* * $Id: dmaxlk.F,v 1.1.1.1 1996/04/01 15:02:20 mclareni Exp $ * * $Log: dmaxlk.F,v $ * Revision 1.1.1.1 1996/04/01 15:02:20 mclareni * Mathlib gen * * #include "gen/pilot.h" SUBROUTINE DMAXLK(SUB,K,M,N,NX,X,A,AL,AU,MODE,EPS,MAXIT,IPRT, + MFR,IAFR,PHI,DPHI,W,NERROR) #include "gen/imp64.inc" DIMENSION AL(*),AU(*),A(*),DPHI(*),IAFR(*),X(*),W(*) EXTERNAL SUB *********************************************************************** * LEAMAX, VERSION: 15.03.1993 *********************************************************************** * * DMAXLK IS THE STEERING ROUTINE FOR MAXIMUM LIKELIHOOD ESTIMATION. * * SUBROUTINE CALLED: D501L2 * * * THE CONSTANTS, VARIABLES AND ARRAYS HAVE THE FOLLOWING MEANING. * * SUB NAME OF USER-SUPPLIED SUBROUTINE SUBPROGRAM, DECLARED * EXTERNAL IN THE CALLING PROGRAM. THIS SUBPROGRAM MUST PROVIDE * THE VALUES OF THE FUNCTION AND, IF MODE=1, THE VALUES OF THE * DERIVATIVES (SEE EXAMPLE) . * K (INTEGER) DIMENSION OF A SINGLE DATA POINT (OBSERVATION) X . * M (INTEGER) NUMBER OF DATA POINTS (OBSERVATIONS). * N (INTEGER) NUMBER OF UNKNOWN PARAMETERS A. * NX (INTEGER) DECLARED FIRST DIMENSION OF ARRAY X IN THE * CALLING PROGRAM, WITH NX .GE. K . * X (DOUBLE PRECISION) TWO-DIMENSIONAL ARRAY OF DIMENSION (NX,M). * ON ENTRY, X MUST CONTAIN THE DATA SET (X(I)) (THE I-TH * COLUMN OF X BELONGS TO THE DATA POINT X(I), I=1,...,M). * A (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH N . * ON ENTRY, A MUST CONTAIN THE STARTING VALUES OF THE UNKNOWN * PARAMETERS FOR THE LEVENBERG-MARQUARDT ALGORITHM. * ON EXIT, A CONTAINS AN APPROXIMATION OF THE MINIMUM POINT. * AL (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH N . * ON ENTRY, AL MUST CONTAIN THE LOWER BOUNDS OF A . * AU (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH N . * ON ENTRY, AU MUST CONTAIN THE UPPER BOUNDS OF A . * MODE (INTEGER) * = 0: THE DERIVATIVE IS COMPUTED NUMERICALLY. * = 1: THE DERIVATIVE HAS TO BE EVALUATED IN SUBPROGRAM SUB . * EPS (DOUBLE PRECISION) USER-SUPPLIED TOLERANCE USED TO CONTROL * THE TERMINATION CRITERION. EPS SHOULD BE CHOSEN ACCORDING * TO THE ACCURACY REQUIRED BY THE UNDERLYING PROBLEM AND TO * THE MACHINE ACCURACY ALSO (RECOMMENDED VALUE ON ENTRY: * 1D-6 ... 1D-12 ). * MAXIT (INTEGER) MAXIMUM PERMITTED NUMBER OF ITERATIONS. * IPRT (INTEGER) PRINTING CONTROL. G * = 0 : NO PRINTING OF INTERMEDIATE RESULTS * = +/- L : PRINTING OF INTERMEDIATE RESULTS AT EVERY ABS(L)-TH * ITERATION; IF IPRT < 0, PRINTING OF ALL INPUT * PARAMETERS OF DSUMSQ IN ADDITION. * MFR (INTEGER) ON EXIT, MFR CONTAINS THE NUMBER OF FREE VARIABLES * AT THE SOLUTION POINT. * IAFR (INTEGER) ONE-DIMENSIONAL ARRAY OF LENGTH >= N , USED AS * WORKING SPACE. ON EXIT, THE FIRST MFR ELEMENTS OF IAFR * CONTAIN THE INDICES OF THE FREE VARIABLES AT THE SOLUTION * POINT. * PHI (DOUBLE PRECISION) ON EXIT, PHI CONTAINS THE VALUE OF THE * OBJECTIVE FUNCTION AT THE MINIMUM POINT. * DPHI (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH N . * ON EXIT, DPHI CONTAINS THE DERIVATIVES OF THE OBJECTIVE * FUNCTION WITH RESPECT TO A (THE GRADIENT) AT THE LAST * ITERATION POINT. * W (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH * 7*N+2*N*N , USED AS WORKING SPACE. * NERROR (INTEGER) ERROR INDICATOR. ON EXIT: * = 0: NO ERROR OR WARNING DETECTED. * = 1: AT LEAST ONE OF THE CONSTANTS K, M, N, NX, MAXIT IS * ILLEGAL OR AT LEAST FOR ONE J THE RELATION * AL(J) .LE. AU(J) IS NOT TRUE. * = 2: THE MAXIMUM NUMBER MAXIT OF ITERATIONS HAS BEEN * REACHED. * = 3: THE OBJECTIVE FUNCTION PHI OR ITS DERIVATIVE IS NOT * DEFINED FOR THE CURRENT VALUES OF THE UNKNOWN * PARAMETER VECTOR A. * = 4: THE ROUTINE DSINV (F012) WAS UNABLE TO SOLVE THE * NORMAL EQUATIONS. * ************************************************************************* * * THE FOLLOWING SUBROUTINE IS A SIMPLE EXAMPLE FOR SUB. * * SUBROUTINE SUB(K,X,N,A,F,DF,MODE,NERROR) * * IMPLICIT DOUBLE PRECISION (A-H,O-Z) * DIMENSION A(*),X(*),DF(*) * PARAMETER (PIR = 0.56418 95835 47756 287D0) * NERROR=1 * IF(A(1) .EQ. 0) RETURN * T=0.5D0*((X(1)-1)/A(1))**2) * F=PIR*EXP(-T)/A(1) * IF(F .EQ. 0) RETURN * NERROR=0 * IF(MODE .NE. 0) DF(1)=-F*(1-2*T)/A(1)**2 * RETURN * END * ************************************************************************* CALL D501L2(K,M,X,NX,MODE,EPS,MAXIT,IPRT,N,A,AL,AU, 1 PHI,DPHI,IAFR,MFR,W(1),W(N+1),W(2*N+1),W(3*N+1), 2 W(4*N+1),W(5*N+1),W(6*N+1),W(7*N+1),W(7*N+N*N+1), 3 SUB,NERROR) RETURN END