! Last change: DOS 29 Jul 2000 4:50 pm ! *** copyright 2000 *** ! *** filename finder.f95 *** John F. Monahan ** ! ********************** ! ! SOME COMPILERS WANT MODULE AVAILABLE FOR ALL PROGRAMS, SO THIS ! MODULE IS AT THE TOP -- USED BY PROGRAM FINDER AND ALL SUBPROGRAMS MODULE PASSINFO ! USE THIS MODULE TO PASS INFORMATION REAL :: SUM ! FROM MAIN PROGRAM (FINDER) TO THE INTEGER :: N ! CALCULATION OF F WITHOUT GOING THROUGH END MODULE PASSINFO ! THE BLACK BOX ROOT-FINDING ROUTINES ! PROGRAM FINDER ! DEMONSTRATION PROGRAM OF SEARCH ROUTINES ! PROBLEM IS TO FIND MAXIMUM OF LIKELIHOOD FUNCTION FOR THE ! LOGARITHMIC SERIES DISTRIBUTION ! ! METHODS HERE ARE GOLDEN SECTION FOR MAXIMIZATION, AND ! BISECTION, REGULA FALSI(FALSE POSITION), ILLINOIS, ! SECANT AND NEWTON FOR ROOT-FINDING ! ! USE INFORMATION HERE USE PASSINFO ! FIRST EXTERNAL SINCE PASSING SUPROGRAM NAMES REAL, EXTERNAL :: F,G,GP ! REAL YL,YR ! 21 FORMAT(' MAXIMUM LIKELIHOOD ESTIMATE=',F12.6,6X, & & 'MAX OF LIKELIHOOD=',F12.6//) 22 FORMAT(' MAXIMUM LIKELIHOOD ESTIMATE=',F12.6,6X, & & 'DERIVATIVE OF LIKE=',F12.6/) 25 FORMAT(/' GOLDEN SECTION FOR MAXIMIZATION') 26 FORMAT(/' BISECTION FOR ROOT OF LIKELIHOOD EQUATION') 27 FORMAT(/' REGULA FALSI FOR ROOT OF LIKELIHOOD EQUATION') 28 FORMAT(/' ILLINOIS METHOD FOR ROOT OF LIKELIHOOD EQUATION') 29 FORMAT(/' SECANT FOR ROOT OF LIKELIHOOD EQUATION') 30 FORMAT(/' NEWTON FOR ROOT OF LIKELIHOOD EQUATION') ! OPEN( UNIT=6, FILE='finder.out' ) ! ! GIVE NECESSARY STATISTICS FOR PROBLEM N = 10 SUM = 15 ! ! DO 15 ITERATIONS OF GOLDEN SECTION SEARCH ON (YMIN,YMAX) ! NEED STARTING VALUES, BEWARE OF ENDPOINTS YL = .02 YR = .98 ! ASSUMING F IS UNIMODAL IN INTERVAL, DOESN'T CHECK! CALL GOLDEN(F,15,YL,YR) WRITE(6,25) WRITE(6,21) YL,YR ! ! DO 15 ITERATIONS OF BISECTION SEARCH ON (YMIN,YMAX) ! AGAIN NEED STARTING VALUES, BEWARE OF ENDPOINTS YL = .02 YR = .98 ! ASSUMING G HAS SIGN CHANGE IN INTERVAL, DOESN'T CHECK! CALL BISECT(G,YL,YR,15,1.E-6) WRITE(6,26) WRITE(6,22) YL,YR YR = F(YL) WRITE(6,21) YL,YR ! ! DO 15 ITERATIONS OF REGULA FALSI SEARCH ON (YMIN,YMAX) ! AGAIN NEED STARTING VALUES, BEWARE OF ENDPOINTS YL = .02 YR = .98 ! ASSUMING G HAS SIGN CHANGE IN INTERVAL, DOESN'T CHECK! CALL REGULA(G,YL,YR,15,1.E-6,1.E-6) WRITE(6,27) WRITE(6,22) YL,YR YR = F(YL) WRITE(6,21) YL,YR ! ! DO 15 ITERATIONS OF ILLINOIS METHOD SEARCH ON (YMIN,YMAX) ! AGAIN NEED STARTING VALUES, BEWARE OF ENDPOINTS YL = .02 YR = .98 ! ASSUMING G HAS SIGN CHANGE IN INTERVAL, DOESN'T CHECK! CALL ILLINI(G,YL,YR,15,1.E-6,1.E-6) WRITE(6,28) WRITE(6,22) YL,YR YR = F(YL) WRITE(6,21) YL,YR ! ! DO 15 ITERATIONS OF SECANT SEARCH ! AGAIN NEED STARTING VALUES, TAKE 'EM IN THE MIDDLE YL = .40 YR = .60 ! NOT MAKING MANY ASSUMPTIONS HERE CALL SECANT(G,YL,YR,15,1.E-6,1.E-6) WRITE(6,29) WRITE(6,22) YL,YR YR = F(YL) WRITE(6,21) YL,YR ! ! DO 15 ITERATIONS OF NEWTON SEARCH YL = .50 ! ONLY NEED ONE STARTING VALUE HERE FOR NEWTON CALL NEWTON(G,GP,YL,15,1.E-6,1.E-6) WRITE(6,30) YR = G(YL) WRITE(6,22) YL,YR YR = F(YL) WRITE(6,21) YL,YR STOP END PROGRAM FINDER REAL FUNCTION F(T) ! FUNCTION TO BE MAXIMIZED USE PASSINFO REAL, INTENT(IN) :: T F = SUM*LOG(T) - N*LOG( - LOG(1.-T) ) RETURN END FUNCTION F REAL FUNCTION G(T) ! FUNCTION FOR WHICH ROOT IS TO BE FOUND USE PASSINFO REAL, INTENT(IN) :: T G = SUM/T + N/( (1.-T)*LOG(1.-T) ) RETURN END FUNCTION G REAL FUNCTION GP(T) ! derivative of FUNCTION FOR WHICH ROOT IS TO BE FOUND USE PASSINFO REAL, INTENT(IN) :: T GP = - SUM/(T*T) + N/( (1.-T)*(1.-T)*LOG(1.-T) ) & & + N/( (1.-T)*(1.-T)*ALOG(1.-T)*LOG(1.-T) ) RETURN END FUNCTION GP SUBROUTINE GOLDEN(F,MIT,XL,XR) ! GOLDEN SECTION SEARCH FOR THE MAX OF F ON THE INTERVAL (XL,XR) ! F MUST BE UNIMODAL; ON RETURN XL=MODE, XR=VALUE OF F THERE IMPLICIT NONE REAL F REAL, INTENT(IN OUT) :: XL,XR INTEGER, INTENT(IN) :: MIT REAL, PARAMETER :: TAU = .618034 ! GOLDEN RATIO REAL, PARAMETER :: TAUSQ = .381966 REAL XLM,XRM,FXLM,FXRM INTEGER I,IT ! KEEP TRACK OF BETTER INSIDE POINT ! I=1 LEFT, I=-1 RIGHT I=0 DO IT = 1,MIT ! INTERVAL OF UNCERTAINTY IS (XL,XR) XL LE XLM LE XRM LE XR XLM = XL*TAU + XR*TAUSQ XRM = XL*TAUSQ + XR*TAU ! ARE ANY OF THESE ALREADY EVALUATED? IF(I.NE.-1) FXLM=F(XLM) IF(I.NE.1) FXRM=F(XRM) write(6,21) IT,XL,XR,FXLM,FXRM 21 format(' GOLD: IT=',I2,' INT IS (',F8.5,' , ',F8.5,' ), VALS ',& &'INSIDE',2F12.5) ! MAIN COMPARISON IF( FXLM .GT. FXRM ) THEN ! F(XLM) .GT. F(XRM) XR = XRM FXRM = FXLM I = 1 ELSE ! F(XLM) .LT. F(XRM) XL = XLM FXLM = FXRM I = -1 ENDIF END DO ! LOOP ON IT ! ALL DONE, TAKE THE BEST POINT SO FAR IF( I .EQ. -1 ) THEN XL = XRM XR = FXRM ELSE XL = XLM XR = FXLM ENDIF RETURN END SUBROUTINE GOLDEN SUBROUTINE BISECT(F,XL,XR,MIT,EPS) ! BISECTION SEARCH TO FIND THE ROOT OF F ON (XL,XR) ! ON RETURN, XL IS THE ROOT AND XR IS VALUE OF F THERE IMPLICIT NONE REAL F REAL, INTENT(IN OUT) :: XL,XR REAL, INTENT(IN) :: EPS INTEGER, INTENT(IN) :: MIT REAL ESP,FL,FR,XN,FN,S INTEGER IT ! DON'T ALLOW EPSILON TOO SMALL ESP=AMAX1(.000001,EPS) FL=F(XL) ! S=1 IF F INCREASES, S=-1 IF F DECREASES S=1. IF(FL.GT.0.) S=-1. FR=F(XR) DO IT = 1,MIT ! THE MIDPOINT OF THE CURRENT INTERVAL IS THE NEW POINT XN = ( XL + XR ) / 2. FN = F(XN) write(6,21) IT,XN,FN,XL,XR 21 format(' BISECT IT=',I2,' XNEW=',F8.5,' F(XN)=',F9.5,' INT IS ('& &,F8.5,' , ',F8.5,' )') IF( ABS(FN) .LT. ESP) EXIT IF(FN*S .LT. 0.) THEN XL=XN FL=FN ELSE XR=XN FR=FN ENDIF END DO ! LOOP ON IT XL=XN XR=FN RETURN END SUBROUTINE BISECT SUBROUTINE SECANT(F,X1,X2,MIT,EPSX,EPSF) ! THE SECANT METHOD FOR FINDING ROOTS TO F ! RETURNS ROOT IN X1, VALUE OF F THERE IN X2 ! EPSF: STOP IF ABS(F(X)) .LE. EPSF ! EPSX: STOP IF CHANGE IN X IS .LE. EPSX IMPLICIT NONE REAL F REAL, INTENT(IN OUT) :: X1,X2 REAL, INTENT(IN) :: EPSX,EPSF INTEGER, INTENT(IN) :: MIT REAL ESPX,ESPF,FX1,FX2,XN,FXN INTEGER IT ! DON'T LET EPSILONS BE TOO SMALL ESPX=MAX(.000001,EPSX) ESPF=MAX(.000001,EPSF) ! TWO STARTING POINTS FX1=F(X1) FX2=F(X2) DO IT = 1,MIT ! FIND NEW POINT XN=X2-FX2*(X2-X1)/(FX2-FX1) FXN=F(XN) write(6,21) IT,XN,FXN 21 format(' SECANT IT=',I3,' CURRENT X=',F12.6,' F(X)=',F12.6) ! TEST FOR CONVERGENCE IF( ( ABS(FXN) .LT. ESPF) .OR. (ABS(XN-X2) .LT. ESPX) ) EXIT ! UPDATE X1=X2 FX1=FX2 X2=XN FX2=FXN END DO ! LOOP ON IT X1=XN X2=FXN RETURN END SUBROUTINE SECANT SUBROUTINE NEWTON(F,FP,X,MIT,EPSX,EPSF) ! NEWTON'S METHOD FOR FINDING ROOT OF F RETURNED IN X ! EPSF: STOP IF ABS(F(X)) .LE. EPSF ! EPSX: STOP IF CHANGE IN X IS .LE. EPSX IMPLICIT NONE REAL F,FP REAL, INTENT(IN OUT) :: X REAL, INTENT(IN) :: EPSX,EPSF INTEGER, INTENT(IN) :: MIT REAL ESPX,ESPF,FX,FPX,XN INTEGER IT ! DON'T LET EPSILONS GET TOO SMALL ESPX = MAX(.000001,EPSX) ESPF = MAX(.000001,EPSF) DO IT = 1,MIT FX=F(X) write(6,21) IT,X,FX 21 format(' NEWTON IT=',I3,' CURRENT X=',F12.6,' F(X)=',F12.6) ! STOP IF ROOT IS FOUND IF( ABS(FX) .LT. ESPF ) RETURN FPX = FP(X) ! GET DERIVATIVE ! FIND NEW POINT TO GO TO XN = X - FX/FPX ! CONVERGENCE TEST: STOP IF NOT MOVING IF( ABS(XN-X) .LT. ESPX ) RETURN X=XN END DO ! LOOP ON IT RETURN END SUBROUTINE NEWTON SUBROUTINE REGULA(F,XL,XR,MIT,EPSX,EPSF) ! USES REGULA FALSI, SECANT FORMULA BUT USE LAST 2 PTS STRADDLING ROOT ! RETURNS ROOT IN XL, VALUE OF F THERE IN XR ! EPSF: STOP IF ABS(F(X)) .LE. EPSF ! EPSX: STOP IF CHANGE IN X IS .LE. EPSX IMPLICIT NONE REAL F REAL, INTENT(IN OUT) :: XL,XR REAL, INTENT(IN) :: EPSX,EPSF INTEGER, INTENT(IN) :: MIT REAL ESPX,ESPF,FL,S,FR,XN,FN INTEGER IT ! SET LOWER LIMIT ON EPSILONS ESPX = MAX(.000001,EPSX) ESPF = MAX(.000001,EPSF) ! GET DIRECTION FL=F(XL) S=1. ! INCREASING (S=1) OR DECREASING (S=-1) ? IF( FL .GT. 0. ) S = -1. FR=F(XR) ! START LOOP DO IT = 1,MIT ! TEST FOR CONVERGENCE IF(XR-XL .LT. ESPX) EXIT ! FIND NEW POINT XN = XL - FL*(XL-XR)/(FL-FR) ! EVALUATE AT NEW POINT FN = F(XN) write(6,21) IT,XN,FN,XL,XR 21 format(' REGULA IT=',I2,' XNEW=',F8.5,' F(XN)=',F9.6,' INT IS ('& &,F8.5,' , ',F8.5,' )') ! STOP IF ROOT IS CLOSE ENOUGH IF( ABS(FN) .LT. ESPF ) EXIT ! ! BRANCH IF( FN*S .LT. 0. ) THEN XL=XN FL=FN ELSE XR=XN FR=FN ENDIF END DO ! LOOP ON IT ! MOVE ROOT TO XL, FUNCTION VAL TO XR XL=XN XR=FN RETURN END SUBROUTINE REGULA SUBROUTINE ILLINI(F,XL,XR,MIT,EPSX,EPSF) ! USES 'ILLINOIS METHOD' OF ROOT-FINDING, A VARIANT OF REGULA FALSI ! RETURNS ROOT IN XL, VALUE OF F THERE IN XR ! EPSF: STOP IF ABS(F(X)) .LE. EPSF ! EPSX: STOP IF CHANGE IN X IS .LE. EPSX IMPLICIT NONE REAL F REAL, INTENT(IN OUT) :: XL,XR REAL, INTENT(IN) :: EPSX,EPSF INTEGER, INTENT(IN) :: MIT REAL ESPX,ESPF,FL,S,FR,XN,FN INTEGER IT,LAST ! SET LOWER LIMIT ON EPSILONS ESPX = MAX(.000001,EPSX) ESPF = MAX(.000001,EPSF) ! GET DIRECTION FL=F(XL) S=1. ! INCREASING (S=1) OR DECREASING (S=-1) ? IF( FL .GT. 0. ) S = -1. FR=F(XR) LAST = 2 ! FIRST TIME THROUGH ! START LOOP DO IT = 1,MIT ! TEST FOR CONVERGENCE IF(XR-XL .LT. ESPX) EXIT ! FIND NEW POINT XN = XL - FL*(XL-XR)/(FL-FR) ! EVALUATE AT NEW POINT FN = F(XN) write(6,21) IT,XN,FN,XL,XR 21 format(' ILLINI IT=',I2,' XNEW=',F8.5,' F(XN)=',F9.6,' INT IS ('& &,F8.5,' , ',F8.5,' )') ! STOP IF ROOT IS CLOSE ENOUGH IF( ABS(FN) .LT. ESPF ) EXIT ! ! BRANCH IF( FN*S .LT. 0. ) THEN XL=XN FL=FN ! IF MOVED THIS POINT LAST TIME, DECREASE OTHER POINT BY HALF IF( LAST .EQ. 1 ) FR = FR/2. LAST = 1 ! LAST = 1 MOVED LEFT ENDPOINT LAST TIME ELSE XR=XN FR=FN ! IF MOVED THIS POINT LAST TIME, DECREASE OTHER POINT BY HALF IF(LAST.EQ.0) FL=FL/2. LAST=0 ! LAST = 0 MOVED RIGHT ENDPOINT LAST TIME ENDIF END DO ! LOOP ON IT ! MOVE ROOT TO XL, FUNCTION VAL TO XR XL=XN XR=FN RETURN END SUBROUTINE ILLINI ! *** end of filename finder.f95 *********************