! Last change: DOS 2 Aug 2000 9:17 pm ! *** copyright 2000 *** ! *** filename laplace2.f95 *** John F. Monahan ** ! ********************** ! ! *** SOME COMPILERS REQUIRE MODULE AT TOP *** MODULE PASSINFO ! KEY TRACKS COMBINATION OF MOMENTS INTEGER KEY ! 1=(0,0), 2=(1,0), 3=(0,1) END MODULE PASSINFO ! 4=(2,0), 5=(1,1), 6=(0,2) ! PROGRAM LAPLACE2 ! DEMONSTRATION OF LAPLACE APPROXIMATIONS FOR MEAN AND VARIANCE ! OF POSTERIOR DISTRIBUTIONS WITH EXTENDED EXAMPLE 9.4 OF CHAPTER 9 USE PASSINFO IMPLICIT NONE REAL, DIMENSION(6) :: FMAX,DETD2F REAL, DIMENSION(2) :: AMEAN REAL, DIMENSION(3) :: ACOV REAL T1,T2,F0,D1,D2,DD1,DD2,DD3,DET,S1,S2 REAL RLGLIK INTEGER J ! 21 FORMAT(/' ITERATION',I4,' THETA',2F12.6,' LOG-LIKE',F12.4) 22 FORMAT(' GRADIENT VECTOR',2F12.6) 23 FORMAT(/' LAPLACE MEAN APPROXIMATION ',2F12.6) 24 FORMAT(' HESSIAN MATRIX OF LOG-LIKE',3F14.6) 25 FORMAT(' UPDATE STEP ',2F12.6) 26 FORMAT(' COV MATRIX FROM MODAL APPRX',3F14.8) 27 FORMAT(' COV MATRIX FROM LAPLACE APPRX',3F14.8) 28 FORMAT(/' PROBLEM',I3,' NEWTON WITH ANALYTICAL DERIVATIVES') ! PUT OUTPUT HERE OPEN( UNIT=6, FILE='laplace2.out' ) ! KEY KEEPS TRACK OF WHICH FUNCTION DO KEY = 1,6 ! SECOND METHOD IS NEWTON WITH ANALYTIC 2'S WRITE(6,28) KEY ! COMMON STARTING VALUES T1 = .263 T2 = .074 J = 0 F0 = RLGLIK(T1,T2) !* WRITE(6,21) J,T1,T2,F0 ! DO SIX ITERATIONS DO J = 1,6 ! 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 (ITERATIONS) ! GET COVARIANCE MATRIX ESTIMATE FROM MODAL ! APPROXIMATION -- ONLY FOR KEY = 1 IF ( KEY .EQ. 1 ) THEN DD1 = -DD1 / DET DD3 = DD3 / DET DD2 = -DD2 / DET WRITE(6,26) DD2,DD3,DD1 END IF ! ( KEY .EQ. 1 ) ! SAVE RESULTS FOR EACH PROBLEM FMAX(KEY) = F0 DETD2F(KEY) = DET END DO ! LOOP ON KEY (MOMENTS) ! GET POSTERIOR MEAN VECTOR AMEAN(1) = SQRT(DETD2F(1)/DETD2F(2)) * EXP(FMAX(2)-FMAX(1)) AMEAN(2) = SQRT(DETD2F(1)/DETD2F(3)) * EXP(FMAX(3)-FMAX(1)) WRITE(6,23) AMEAN ! COV MATRIX ACOV(1) = SQRT(DETD2F(1)/DETD2F(4)) * EXP(FMAX(4)-FMAX(1)) ACOV(2) = SQRT(DETD2F(1)/DETD2F(5)) * EXP(FMAX(5)-FMAX(1)) ACOV(3) = SQRT(DETD2F(1)/DETD2F(6)) * EXP(FMAX(6)-FMAX(1)) ACOV(1) = ACOV(1) - AMEAN(1)*AMEAN(1) ACOV(2) = ACOV(2) - AMEAN(1)*AMEAN(2) ACOV(3) = ACOV(3) - AMEAN(2)*AMEAN(2) WRITE(6,27) ACOV(1),ACOV(2),ACOV(3) ! DONE STOP END PROGRAM LAPLACE2 REAL FUNCTION RLGLIK(T1,T2) ! COMPUTE LOG-LIKELIHOOD FOR EXTENDED EXAMPLE USE PASSINFO IMPLICIT NONE REAL, INTENT(IN) :: T1,T2 REAL P1,P2,P3,P4 INTEGER N1,N2,N3,N4 ! DATA N1,N2,N3,N4/1,10,4,9/ ! SHOULDN'T NEED TO CHECK FOR P'S OUTSIDE (0,1) 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) ! KEY = 1 IS POSTERIOR, 2,3 FOR FIRST MOMENTS ! 4,5,6 FOR SECOND MOMENTS IF( KEY .EQ. 2 ) RLGLIK = RLGLIK + LOG(T1) IF( KEY .EQ. 3 ) RLGLIK = RLGLIK + LOG(T2) IF( KEY .EQ. 4 ) RLGLIK = RLGLIK + 2.*LOG(T1) IF( KEY .EQ. 5 ) RLGLIK = RLGLIK + LOG(T1*T2) IF( KEY .EQ. 6 ) RLGLIK = RLGLIK + 2.*LOG(T2) ! RETURN END FUNCTION RLGLIK SUBROUTINE DELLL(T1,T2,D1,D2) ! COMPUTE THE GRADIENT FUNCTION OF LOG-LIKELIHOOD USE PASSINFO 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/1,10,4,9/ ! SHOULDN'T NEED TO CHECK IF P'S IN (0,1) 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)) IF( KEY .EQ. 2 ) D1 = D1 + 1./T1 IF( KEY .EQ. 4 ) D1 = D1 + 2./T1 IF( KEY .EQ. 5 ) D1 = D1 + 1./T1 ! D2 = (N1/P1)*(2.*T1) + (N2/P2)*(-2.*T1) + & & (N3/P3)*(2.*(1.-T1-T2)) + (N4/P4)*(-2.*(1.-T1-T2)) IF( KEY .EQ. 3 ) D2 = D2 + 1./T2 IF( KEY .EQ. 5 ) D2 = D2 + 1./T2 IF( KEY .EQ. 6 ) D2 = D2 + 2./T2 ! RETURN END SUBROUTINE DELLL SUBROUTINE DEL2L(T1,T2,DD1,DD2,DD3) ! COMPUTE HESSIAN MATRIX ANALYTICALLY USE PASSINFO 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/1,10,4,9/ ! SHOULDN'T NEED TO CHECK IF P'S IN (0,1) 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) IF( KEY .EQ. 2 ) DD1 = DD1 - 1./(T1*T1) IF( KEY .EQ. 4 ) DD1 = DD1 - 2./(T1*T1) IF( KEY .EQ. 5 ) DD1 = DD1 - 1./(T1*T1) ! 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) IF( KEY .EQ. 3 ) DD2 = DD2 - 1./(T2*T2) IF( KEY .EQ. 5 ) DD2 = DD2 - 1./(T2*T2) IF( KEY .EQ. 6 ) DD2 = DD2 - 2./(T2*T2) ! 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 ! *** end of filename laplace2.f95 *********************