!     Last change:  DOS  30 Jul 2000    2:06 pm
!                                       *** copyright 2000 ***
! *** filename conclk.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 conclk
!  Example of concentrated likelihood -- given one parameter,
!   problem becomes regression through the origin
!  From Bates and Watts, p. 41, Example: BOD 1
!
      use passinfo
      implicit none
      real, external :: sse
      real xl,xr
      integer mit
!
  20  format(f3.0,f6.1)
!
      open( unit=6, file='conclk.out' )
!                    Golden Section to find max of -sse
!                    starting values
      xl = 0.
      xr = 1.
      mit = 25
      call golden(sse,mit,xl,xr)
      stop
      end program conclk
      real function sse(t2)
      use passinfo
      implicit none
      real t1,t2,z,yz,zz
      integer i
!
  21  format(/' in sse, estimates',2f8.4,'  -sse',f9.4,'  var est',f9.4)
!                        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 -- ML divides by n
      z = sse / real(nobs)
!                        flip sign to maximize -sse
      sse = - sse
!                        write out results
      write(6,21) t1,t2,sse,z
      return
      end function sse
      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 conclk.f95        *********************