!     Last change:  DOS   2 Aug 2000    9:16 pm
!                                       *** copyright 2000 ***
! *** filename laplace1.f95             *** John F. Monahan **
!                                       **********************
!
! *** SOME COMPILERS REQUIRE MODULE AT TOP ***
!
      MODULE PASSINFO
      REAL :: SUMY
      INTEGER :: N
      END MODULE PASSINFO
!
      PROGRAM LAPLACE1
!  DEMONSTRATION PROGRAM OF LAPLACE'S METHOD
!   PROBLEM IS TO FIND MAXIMUM OF LIKELIHOOD FUNCTION FOR THE
!   LOGARITHMIC SERIES DISTRIBUTION
!
!                   USE INFORMATION HERE
      USE PASSINFO
      IMPLICIT NONE
!                   FIRST EXTERNAL SINCE PASSING SUPROGRAM NAMES
      REAL, EXTERNAL :: F,G,GP
!
      REAL YL,YR,FMAX0,FMAX1,FMAX2,D2F0,D2F1,D2F2,AMEAN,ASCND,AVAR
!
  21  FORMAT(' SECOND DERIV AT MAX=',F12.6,6X,              &
     &        'MAX OF LIKELIHOOD=',F12.6//)
  22  FORMAT(' MAXIMUM LIKELIHOOD ESTIMATE=',F12.6,6X,      &
     &        'DERIVATIVE OF LIKE=',F12.6/)
  24  FORMAT(' LAPLACE APPROXIMATIONS',/' MEAN',F12.6,' SECOND MOMENT'&
     &       ,F12.6,'  VARIANCE',F12.6)
  28  FORMAT(/' ILLINOIS METHOD FOR ROOT OF LIKELIHOOD EQUATION')
!
      OPEN( UNIT=6, FILE='laplace1.out' )
!
!  DO 15 ITERATIONS OF ILLINOIS METHOD SEARCH ON (YMIN,YMAX)            
!              AGAIN NEED STARTING VALUES, BEWARE OF ENDPOINTS
      YL = .02
      YR = .98   
!              FIND MAX OF DENOMINATOR
      SUMY = 15.
      N = 10
!              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                                                 
      FMAX0 = F(YL)
      D2F0 = GP(YL)                                                      
      WRITE(6,21) D2F0,FMAX0                                                 
!
!              FIND MAX OF NUMERATOR FOR MEAN
      SUMY = 16.
      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                                                 
      FMAX1 = F(YL)
      D2F1 = GP(YL)                                                      
      WRITE(6,21) D2F1,FMAX1                                              
!              FIND MAX OF NUMERATOR FOR SECOND MOMENT
      SUMY = 17.
!
      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                                                 
      FMAX2 = F(YL)
      D2F2 = GP(YL)                                                      
      WRITE(6,21) D2F2,FMAX2
!              FIND LAPLACE APPROXIMATIONS FOR MEAN AND VARIANCE
      AMEAN = SQRT(D2F0/D2F1) * EXP( FMAX1 - FMAX0 )
      ASCND = SQRT(D2F0/D2F2) * EXP( FMAX2 - FMAX0 )
      AVAR = ASCND - AMEAN*AMEAN
      WRITE(6,24) AMEAN,ASCND,AVAR
      STOP                                                              
      END PROGRAM LAPLACE1
      REAL FUNCTION F(T)
!  FUNCTION TO BE MAXIMIZED
      USE PASSINFO
      REAL, INTENT(IN) :: T
      F = LOG(1.-T) + (SUMY+1.)*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 = -1./(1.-T) + (SUMY+1.)/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 = -1./((1.-T)*(1.-T))                                   &
     &        - (SUMY+1.)/(T*T) + N/( (1.-T)*(1.-T)*LOG(1.-T) )  &
     &        + N/( (1.-T)*(1.-T)*LOG(1.-T)*LOG(1.-T) )
      RETURN                                                            
      END FUNCTION GP
      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 laplace1.f95      *********************