!     Last change:  DOS  30 Jul 2000    1:56 pm
!                                       *** copyright 2000 ***
! *** filename chex94rp.f95             *** John F. Monahan **
!                                       **********************
      PROGRAM CHEX94RP
!  EXTENDED EXAMPLE 9.4 OF CHAPTER 9  -- REPARAMETERIZED
!  DEMONSTRATION OF ONE METHOD FOR COMPUTING MAXIMUM LIKELIHOOD
!   ESTIMATES -- NEWTON'S METHOD WITH NUMERICAL DERIVATIVES
      IMPLICIT NONE
      REAL, EXTERNAL :: RLGLKA
      REAL A1,A2,F0,D1,D2,DD1,DD2,DD3,DET,S1,S2         
      INTEGER KREP
!
  21  FORMAT(/' ITERATION',I4,' ALPHA',2F12.6,'  LOG-LIKE',F12.4)
  22  FORMAT('   GRADIENT VECTOR',2F12.6)
  24  FORMAT('   HESSIAN MATRIX OF LOG-LIKE',3F14.6)
  25  FORMAT('   UPDATE STEP ',2F12.6)
  26  FORMAT('  COVARIANCE MATRIX FOR ESTIMATES',3F14.8)
  29  FORMAT(/' NEWTON WITH NUMERICAL DERIVATIVES')
!
      OPEN( UNIT=6, FILE='chex94rp.out' )
!                      METHOD IS NEWTON WITH NUMERICAL 2'S
      WRITE(6,29)
!                      COMMON STARTING VALUES
      A1 = 0.0
      A2 = 0.0                                                         
      F0 = RLGLKA(A1,A2)
      WRITE(6,21) 0,A1,A2,F0
!                      DO FIVE ITERATIONS
      DO KREP = 1,5
!                      GET GRADIENT VECTOR
      CALL DELLA(A1,A2,D1,D2)
      WRITE(6,22) D1,D2                                                 
!                      HESSIAN MATRIX NUMERICALLY
      CALL DEL2LN(RLGLKA,A1,A2,DD1,DD2,DD3)
      WRITE(6,24) DD1,DD2,DD3                                        
      DET = DD1*DD2 - DD3*DD3                                       
!                      COMPUTE THE STEP IN S1,S2                 
      S1 = (DD2*D1 - DD3*D2) / DET
      S2 = (-DD3*D1 + DD1*D2) / DET                                   
      WRITE(6,25) S1,S2                                               
!                      GET NEW POINT
      A1 = A1 - S1
      A2 = A2 - S2                                                      
      F0 = RLGLKA(A1,A2)                                                
      WRITE(6,21) KREP,A1,A2,F0
      END DO  ! LOOP ON KREP
!                      GET COVARIANCE MATRIX ESTIMATE
      DD1 = -DD1 / DET
      DD3 = DD3 / DET 
      DD2 = -DD2 / DET
      WRITE(6,26) DD2,DD3,DD1
!                      DONE                           
      STOP
      END PROGRAM CHEX94RP
      REAL FUNCTION RLGLIK(T1,T2)
!  COMPUTE LOG-LIKELIHOOD FOR EXTENDED EXAMPLE
      IMPLICIT NONE
      REAL, INTENT(IN) :: T1,T2
      REAL P1,P2,P3,P4
      INTEGER N1,N2,N3,N4                                               
      DATA N1,N2,N3,N4/17,182,60,176/                                   
      P1 = 2.*T1*T2                                                     
      P2 = T1*(2. - T1 - 2.*T2)                                         
      P3 = T2*(2. - T2 - 2.*T1)                                         
      P4 = (1.-T1-T2)*(1.-T1-T2)                                        
!
      RLGLIK = N1*LOG(P1) + N2*LOG(P2) + N3*LOG(P3) + N4*LOG(P4)
!
      RETURN                                                            
      END FUNCTION RLGLIK
      SUBROUTINE TALPH(T1,T2,A1,A2)
!  CONVERT A'S TO T'S -- REPARAMETERIZATION
      IMPLICIT NONE
      REAL, INTENT(OUT) :: T1,T2
      REAL, INTENT(IN) :: A1,A2
      REAL S
      T1 = EXP(A1)
      T2 = EXP(A2)
      S = T1 + T2 + 1.
      T1 = T1 / S
      T2 = T2 / S
      RETURN
      END SUBROUTINE TALPH
      REAL FUNCTION RLGLKA(A1,A2)
!  LOG-LIKELIHOOD FUNCTION OF ALPHA'S
!  CALL RLGLIK WITH T'S
      IMPLICIT NONE
      REAL, INTENT(IN) :: A1,A2
      REAL T1,T2,RLGLIK
      CALL TALPH(T1,T2,A1,A2)
      RLGLKA = RLGLIK(T1,T2)
      RETURN
      END FUNCTION RLGLKA
      SUBROUTINE DELLA(A1,A2,D1,D2)
!  COMPUTE GRADIENT FUNCTION OF LOG-LIKELIHOOD IN TERMS OF ALPHA
      IMPLICIT NONE
      REAL, INTENT(IN) :: A1,A2
      REAL, INTENT(OUT) :: D1,D2
      INTEGER N1,N2,N3,N4,N
      REAL EA1,EA2
      DATA N1,N2,N3,N4,N/17, 182, 60, 176, 435 /
!  FIRST TRANSFORM
      EA1 = EXP(A1)
      EA2 = EXP(A2)
!  GRADIENT IS ONLY MODERATELY COMPLICATED
      D1 = N1 + N2*(2.*(1.+EA1)/(2.+EA1)) - 2.*N*(EA1/(1.+EA1+EA2))
      D2 = N1 + N3*(2.*(1.+EA2)/(2.+EA2)) - 2.*N*(EA2/(1.+EA1+EA2))
      RETURN
      END SUBROUTINE DELLA
      SUBROUTINE DEL2LN(FN,T1,T2,DD1,DD2,DD3)
!  COMPUTE HESSIAN MATRIX NUMERICALLY WITH CENTRAL DIFFERENCES
      IMPLICIT NONE
      REAL FN
      REAL, INTENT(IN) :: T1,T2
      REAL, INTENT(OUT) :: DD1,DD2,DD3
      REAL, DIMENSION(3,3) :: T
      REAL EP1,EP2,TT1,TT2
      INTEGER I1,I2
!  BE CAREFUL THAT THE DIFFERENCE ISN'T ZERO
      EP1 = MAX( 0.01, 0.01*ABS(T1) )
      EP2 = MAX( 0.01, 0.01*ABS(T2) )
!                    FULL CENTRAL DIFFERENCES
      DO I1 = 1,3
      DO I2 = 1,3
!                    USE .01 RELATIVE FOR DIFFERENCE H
      TT1 = EP1*REAL(I1-2) + T1
      TT2 = EP2*REAL(I2-2) + T2
!                    FN HERE IS FOR LIKELIHOOD FUNCTION
      T(I1,I2) = FN(TT1,TT2)
      END DO  ! LOOP ON I2
      END DO  ! LOOP ON I1
!                    11 ELEMENT
      DD1 = ( (T(1,1)-T(2,1)) - (T(2,1)-T(3,1))                       &
     &    +  (T(1,2)-T(2,2)) - (T(2,2)-T(3,2))                        &
     &    +  (T(1,3)-T(2,3)) - (T(2,3)-T(3,3)) )/(EP1*EP1*3.)
!                    22 ELEMENT
      DD2 = ( (T(1,1)-T(1,2)) - (T(1,2)-T(1,3))                       &
     &    +  (T(2,1)-T(2,2)) - (T(2,2)-T(2,3))                        &
     &    +  (T(3,1)-T(3,2)) - (T(3,2)-T(3,3)) )/(EP2*EP2*3.)
!                    12 ELEMENT USES CORNER POINTS
      DD3 = ( (T(1,1) - T(1,3)) - (T(3,1) - T(3,3)) )/(EP1*EP2*4.)
      RETURN                                                            
      END SUBROUTINE DEL2LN
! *** end of filename chex94rp.f95      *********************