!     Last change:  DOS  30 Jul 2000    2:08 pm
!                                       *** copyright 2000 ***
! *** filename concll.f95               *** John F. Monahan **
!                                       **********************
!
!  SOME COMPILERS REQUIRE MODULE TO PREVIOUSLY EXIST OR BE
!   AT THE BEGINNING
      module passinfo
      real, dimension(6) :: y,x
      integer, parameter :: nobs = 6
      data x/ 1., 2., 3., 4., 5., 7./
      data y/ 8.3, 10.3, 19.0, 16.0, 15.6, 19.8/
      end module passinfo
!
      program concll
!  Example of concentrated likelihood -- given one parameter,
!   problem becomes regression through the origin
!  Recoded so that covariances are easier
!
!  From Bates and Watts, p. 41, Example: BOD 1
!
      implicit none
      real, external :: rlc
      real xl,xr
      integer mit
!
  20  format(f3.0,f6.1)
!
      open( unit=6, file='concll.out' )
!                    Golden Section to find max of rlc
!                    starting values
      xl = 0.
      xr = 1.
      mit = 25
      call golden(rlc,mit,xl,xr)
      stop
      end program concll
      real function rlc(t2)
!  concentrated likelihood function -- just of t2
      implicit none
      real t1,t2,t3,rll
!                         get MLE's for t1 and t3 from t2
      call t1t3(t1,t2,t3)
!                         get log-likelihood
      rlc = rll(t1,t2,t3)
      return
      end
      real function rll(t1,t2,t3)
!  compute the log likelihood for given values of t1,t2,t3
      use passinfo
      implicit none
      real, intent(in) :: t1,t2,t3
      real s,z
      integer i
!
  21  format(/' in rll * t1,t2,t3 are',3f9.4,' rll = ',f12.6)
!
      s = 0.
      do i = 1,nobs
      z = 1. - exp( -t2*x(i) )
      z = y(i) - t1*z 
      s = s + z*z
      end do  ! loop on i
      rll = -( nobs*log(t3) + s/t3 )/2.
!  write out results
      write(6,21) t1,t2,t3,rll
      return
      end function rll
      subroutine t1t3(t1,t2,t3)
!  computes mle's for t1 and t3 for given t2
      use passinfo
      implicit none
      real, intent(in) :: t2
      real, intent(out) :: t1,t3
      real z,yz,zz,sse
      integer i
!                        first do regression to get t1
      yz = 0.
      zz = 0.
      do i = 1,nobs
      z = 1. - exp( -t2*x(i) )
      zz = zz + z*z
      yz = yz + y(i)*z
      end do  ! loop on i
!                        estimate
      t1 = yz/zz
!                        compute sse
      sse = 0.
      do i = 1,nobs
      z = y(i) - t1*(1. - exp( -t2*x(i) ))
      sse = sse + z*z
      end do  ! loop on i
!                        variance estimate
      t3 = sse / real(nobs)
      return
      end subroutine t1t3
      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
! *** end of filename concll.f95        *********************