!     Last change:  DOS  31 Jul 2000    7:50 pm
!                                       *** copyright 2000 ***
! *** filename chex102s.f95             *** John F. Monahan **
!                                       **********************
      program chex102s
!  ext102 -- integration of posterior from Chapter 9's Extended Example
!             using Simpson's rule centered at posterior mode
!             and scaled to go out 4 stdev in each parameter direction
!
      implicit none
      real, dimension(6) :: ss
      real t1,t2,tm1,tm2,se1,se2,rll,pstmx,f,w1,w2,wgt,fint
      real rlglik
      integer kint,i,j,np
!
  28  format(/'Simpson''s rule',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='chex102s.out' )
!                   mode of posterior
      tm1 = .265782
      tm2 = .110980
!                   standard deviations of normal approximation
      se1 = .0694
      se2 = .0469
!                   max of log posterior
      pstmx = -28.0271
!                   number of intervals on a side
      do kint = 2,10
!                   initialize to zero
      ss = 0.
!                   how many points for Simpson
      np = 2*kint+1
      fint = real(kint)
      do i = 1,np
!                   set the weight for Simpson as 1 4 2 4 2 4 2 4 1
      w1 = 2.*real(2 - mod(i,2))
      if( (i .eq. 1) .or. (i .eq. np) ) w1 = 1.
!                   go out 4 stdev to make sure to cover (& hope)
      t1 = tm1 + 4.*se1*real(i-1-kint)/fint
      do j = 1,np
!                   set the weight for Simpson as 1 4 2 4 2 4 2 4 1
      w2 = 2.*real(2 - mod(j,2))
      if( (j .eq. 1) .or. (j .eq. np) ) w2 = 1.
!                   go out 4 stdev to make sure to cover (& hope)
      t2 = tm2 + 4.*se2*real(j-1-kint)/fint
      wgt = w1*w2
      rll = rlglik(t1,t2) - pstmx
      f = 0.
      if( rll .gt. -30. ) f = exp(rll)*wgt
!                   Simpson rule -- use all points
      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 -- Simpson
      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) = ss(1) * (8.*se1)* (8.*se2) / (36.*real(kint*kint))
      write(6,28) np*np,ss(1)
      write(6,29) ss(2),ss(4),ss(3),ss(5),ss(6)
      end do  ! loop on kint
      stop
      end program chex102s
      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 chex102s.f95      *********************