!     Last change:  DOS  31 Jul 2000    7:45 pm
!                                       *** copyright 2000 ***
! *** filename chex102g.f95             *** John F. Monahan **
!                                       **********************
      program chex102g
!  ext102 -- integration of posterior from Chapter 9's Extended Example
!             using Gauss-Hermite quadrature centered at posterior mode
!             and scaled posterior normal approximation
!
      implicit none
      real, dimension(6) :: ss
      real, dimension(10) :: absc,wgt
      real t1,t2,tm1,tm2,h11,h12,h22,rll,pstmx,f
      real rlglik
      integer i,j,np
!
  22  format(4x,2f16.8)
  28  format(/'Gauss-Hermite',i6,' points -- Normalization constant', &
     &       f12.8/10x,'Posterior mean',10x,'Covariance Matrix')
  29  format('  theta(1)',f12.6,6x,f12.6/                             &
     &       '  theta(2)',f12.6,6x,2f12.6)
!                   output unit
      open( unit=6, file='chex102g.out' )
!                   file with Gauss-Hermite abscissas & weights
      open( unit=8, file='qgaushm.tab' )
!
!                   mode of posterior
      tm1 = .265782
      tm2 = .110980
!                   hessian matrix from normal approximation
      h11 = 215.1117
      h12 =  59.4624
      h22 = 472.0310
!                   max of log posterior
      pstmx = -28.0271
!                   find Cholesky factor of hessian
      h11 = sqrt(h11)
      h12 = h12 / h11
      h22 = sqrt( h22 - h12*h12 )
!*      write(6,22) h11,h12,h22
!                   number of points on a side
      do np = 3,10
!                   initialize to zero
      ss = 0.
!                   read table
      do i = 1,np
      read(8,22) absc(i),wgt(i)
!*      write(6,22) absc(i),wgt(i)
      end do  ! loop on i
!
!                   use i for first component, j for second
      do i = 1,np
      do j = 1,np
!                   first L-inv-transpose * z
      t2 = absc(j) / h22
      t1 = (absc(i) - h12*t2) / h11
!                   now center at mode
      t1 = t1 + tm1
      t2 = t2 + tm2
      rll = rlglik(t1,t2) - pstmx + (absc(i)*absc(i)+absc(j)*absc(j))/2.
      f = 0.
      if( rll .gt. -30. ) f = exp(rll)*wgt(i)*wgt(j)
!                   Gauss-Hermite -- remember to divide
      ss(1) = ss(1) + f
      ss(2) = ss(2) + f*t1
      ss(3) = ss(3) + f*t2
      ss(4) = ss(4) + f*t1*t1
      ss(5) = ss(5) + f*t1*t2
      ss(6) = ss(6) + f*t2*t2
      end do  ! loop on j
      end do  ! loop on i
!                   posterior moments -- Gauss-Hermite
      ss(2) = ss(2) / ss(1)
      ss(3) = ss(3) / ss(1)
      ss(4) = ss(4) / ss(1) - ss(2)*ss(2)
      ss(5) = ss(5) / ss(1) - ss(2)*ss(3)
      ss(6) = ss(6) / ss(1) - ss(3)*ss(3)
!                   fix normalization constant
      ss(1) = 6.28318 * ss(1) / (h11 * h22)
      write(6,28) np*np,ss(1)
      write(6,29) ss(2),ss(4),ss(3),ss(5),ss(6)
      end do  ! loop on np
      stop
      end program chex102g
      REAL FUNCTION RLGLIK(T1,T2)
!  COMPUTE LOG-LIKELIHOOD FOR EXTENDED EXAMPLE
!   *** modified from Example 9.4 by dropping sample size ***
      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/                                   
!                   **************** numbers are changed from Chapter 9
      P1 = 2.*T1*T2                                                     
      P2 = T1*(2. - T1 - 2.*T2)                                         
      P3 = T2*(2. - T2 - 2.*T1)                                         
      P4 = (1.-T1-T2)*(1.-T1-T2)                                        
!                   CHECK IF OUT OF RANGE
!                   ASSIGN TINY LIKELIHOOD IF OUTSIDE TRIANGLE
      RLGLIK = -1.E30
      IF( P1 .LE. 0. ) RETURN
      IF( P2 .LE. 0. ) RETURN
      IF( P3 .LE. 0. ) RETURN
      IF( P4 .LE. 0. ) RETURN
!
      RLGLIK = N1*LOG(P1) + N2*LOG(P2) + N3*LOG(P3) + N4*LOG(P4)
!
      RETURN                                                            
      END FUNCTION RLGLIK
! *** end of filename chex102g.f95      *********************