! Last change: DOS 30 Jul 2000 1:54 pm ! *** copyright 2000 *** ! *** filename chex94.f95 *** John F. Monahan ** ! ********************** PROGRAM CHEX94 ! EXTENDED EXAMPLE 9.4 OF CHAPTER 9 ! DEMONSTRATION OF THREE METHODS FOR COMPUTING MAXIMUM LIKELIHOOD ! ESTIMATES -- SCORING, AND NEWTON'S METHOD WITH ANALYTICAL AND ! NUMERICAL DERIVATIVES IMPLICIT NONE REAL, EXTERNAL :: RLGLIK REAL T1,T2,F0,D1,D2,DD1,DD2,DD3,DET,J1,J2,J3,S1,S2 INTEGER J,N ! 21 FORMAT(/' ITERATION',I4,' THETA',2F12.6,' LOG-LIKE',F12.4) 22 FORMAT(' GRADIENT VECTOR',2F12.6) 23 FORMAT(' FISHER''S INFORMATION (11, 12, 22)',3F14.6) 24 FORMAT(' HESSIAN MATRIX OF LOG-LIKE',3F14.6) 25 FORMAT(' UPDATE STEP ',2F12.6) 26 FORMAT(' COVARIANCE MATRIX FOR ESTIMATES',3F14.8) 27 FORMAT(/' FIRST METHOD -- SCORING ') 28 FORMAT(/' SECOND METHOD -- NEWTON WITH ANALYTICAL DERIVATIVES') 29 FORMAT(/' THIRD METHOD -- NEWTON WITH NUMERICAL DERIVATIVES') ! OPEN( UNIT=6, FILE='chex94.out' ) ! NUMBER OF OBSERVATIONS N = 435 ! COMMON STARTING VALUES T1 = .263 T2 = .074 J = 0 ! FIRST METHOD IS SCORING WRITE(6,27) ! EVALUATE FUNCTION AT START F0 = RLGLIK(T1,T2) WRITE(6,21) J,T1,T2,F0 ! DO FOUR ITERATIONS DO J = 1,4 ! GET GRADIENT VECTOR CALL DELLL(T1,T2,D1,D2) WRITE(6,22) D1,D2 ! GET FISHER'S INFORMATION CALL INFO(T1,T2,J1,J2,J3) WRITE(6,23) J1,J2,J3 DET = J1*J2 - J3*J3 ! COMPUTE THE STEP IN S1,S2 S1 = (J2*D1 - J3*D2) / (DET*N) S2 = (-J3*D1 + J1*D2) / (DET*N) WRITE(6,25) S1,S2 ! GET NEW POINT T1 = T1 + S1 T2 = T2 + S2 F0 = RLGLIK(T1,T2) WRITE(6,21) J,T1,T2,F0 END DO ! LOOP ON J ! GET COVARIANCE MATRIX ESTIMATE J2 = J2 / (DET*N) J1 = J1 / (DET*N) J3 = - J3 / (DET*N) WRITE(6,26) J2,J3,J1 ! SECOND METHOD IS NEWTON WITH ANALYTIC 2'S WRITE(6,28) ! COMMON STARTING VALUES T1 = .263 T2 = .074 J = 0 F0 = RLGLIK(T1,T2) WRITE(6,21) J,T1,T2,F0 ! DO FOUR ITERATIONS DO J = 1,4 ! GET GRADIENT VECTOR CALL DELLL(T1,T2,D1,D2) WRITE(6,22) D1,D2 ! GET HESSIAN MATRIX CALL DEL2L(T1,T2,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 T1 = T1 - S1 T2 = T2 - S2 F0 = RLGLIK(T1,T2) WRITE(6,21) J,T1,T2,F0 END DO ! LOOP ON J ! GET COVARIANCE MATRIX ESTIMATE DD1 = -DD1 / DET DD3 = DD3 / DET DD2 = -DD2 / DET WRITE(6,26) DD2,DD3,DD1 ! THIRD METHOD IS NEWTON WITH NUMERICAL 2'S WRITE(6,29) ! COMMON STARTING VALUES T1 = .263 T2 = .074 J = 0 F0 = RLGLIK(T1,T2) WRITE(6,21) J,T1,T2,F0 ! DO FOUR ITERATIONS DO J = 1,4 ! GET GRADIENT VECTOR CALL DELLL(T1,T2,D1,D2) WRITE(6,22) D1,D2 ! HESSIAN MATRIX NUMERICALLY CALL DEL2LN(RLGLIK,T1,T2,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 T1 = T1 - S1 T2 = T2 - S2 F0 = RLGLIK(T1,T2) WRITE(6,21) J,T1,T2,F0 END DO ! LOOP ON J ! GET COVARIANCE MATRIX ESTIMATE DD1 = -DD1 / DET DD3 = DD3 / DET DD2 = -DD2 / DET WRITE(6,26) DD2,DD3,DD1 ! DONE STOP END PROGRAM CHEX94 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 DELLL(T1,T2,D1,D2) ! COMPUTE THE GRADIENT FUNCTION OF LOG-LIKELIHOOD IMPLICIT NONE REAL, INTENT(IN) :: T1,T2 REAL, INTENT(OUT) :: D1,D2 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) ! D1 = (N1/P1)*(2.*T2) + (N2/P2)*(2.*(1.-T1-T2)) + & & (N3/P3)*(-2.*T2) + (N4/P4)*(-2.*(1.-T1-T2)) ! D2 = (N1/P1)*(2.*T1) + (N2/P2)*(-2.*T1) + & & (N3/P3)*(2.*(1.-T1-T2)) + (N4/P4)*(-2.*(1.-T1-T2)) ! RETURN END SUBROUTINE DELLL SUBROUTINE INFO(T1,T2,J1,J2,J3) ! COMPUTE FISHER'S INFORMATION MATRIX FOR EXTENDED EXAMPLE IMPLICIT NONE REAL, INTENT(IN) :: T1,T2 REAL, INTENT(OUT) :: J1,J2,J3 REAL P1,P2,P3,P4 P1 = 2.*T1*T2 P2 = T1*(2. - T1 - 2.*T2) P3 = T2*(2. - T2 - 2.*T1) P4 = (1.-T1-T2)*(1.-T1-T2) ! 11 ELEMENT J1 = 4.*T2*T2/P1 + 4.*(1.-T1-T2)*(1.-T1-T2)/P2 + & & 4.*T2*T2/P3 + 4.*(1.-T1-T2)*(1.-T1-T2)/P4 ! 22 ELEMENT J2 = 4.*T1*T1/P1 + 4.*T1*T2/P2 + & & 4.*(1.-T1-T2)*(1.-T1-T2)/P3 + 4.*(1.-T1-T2)*(1.-T1-T2)/P4 ! 12 ELEMENT J3 = 4.*T1*T2/P1 - 4.*(1.-T1-T2)*T1/P2 - & & 4.*T2*(1.-T1-T2)/P3 + 4.*(1.-T1-T2)*(1.-T1-T2)/P4 RETURN END SUBROUTINE INFO SUBROUTINE DEL2L(T1,T2,DD1,DD2,DD3) ! COMPUTE HESSIAN MATRIX ANALYTICALLY IMPLICIT NONE REAL, INTENT(IN) :: T1,T2 REAL, INTENT(OUT) :: DD1,DD2,DD3 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) ! 11 ELEMENT DD1 = N1*(-(4.*T2*T2)/(P1*P1)) & & + N2*(-2./P2-4.*((1.-T1-T2)/P2)**2) & & + N3*(-4.*T2*T2/(P3*P3)) & & + N4*(2/P4-4.*((1.-T1-T2)/P4)**2) ! 22 ELEMENT DD2 = N1*(-(4.*T1*T1)/(P1*P1)) & & + N3*(-2./P3-4.*((1.-T1-T2)/P3)**2) & & + N2*(-4.*T1*T1/(P2*P2)) & & + N4*(2/P4-4.*((1.-T1-T2)/P4)**2) ! 12 ELEMENT DD3 = N1*(2/P1-(4.*T1*T2)/(P1*P1)) & & + N2*(-2.+4.*T1*(1.-T1-T2)/P2)/P2 & & + N3*(-2.+4.*T2*(1.-T1-T2)/P3)/P3 & & + N4*(2/P4-4.*((1.-T1-T2)/P4)**2) ! RETURN END SUBROUTINE DEL2L 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 TT1,TT2 REAL, DIMENSION(3,3) :: T INTEGER I1,I2 ! FULL CENTRAL DIFFERENCES DO I1 = 1,3 DO I2 = 1,3 ! USE .01 RELATIVE FOR DIFFERENCE H TT1 = (1. + .01*(I1-2) )*T1 TT2 = (1. + .01*(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)) )/(T1*T1*3.E-4) ! 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)) )/(T2*T2*3.E-4) ! 12 ELEMENT USES CORNER POINTS DD3 = ( (T(1,1) - T(1,3)) - (T(3,1) - T(3,3)) )/(T1*T2*4.E-4) RETURN END SUBROUTINE DEL2LN ! *** end of filename chex94.f95 *********************