!     Last change:  DOS  31 Jul 2000    7:57 pm
!                                       *** copyright 2000 ***
! *** filename chex103.f95              *** John F. Monahan **
!                                       **********************
      program chex103
!  Integration of posterior from
!    Chapter 10's Variance Components Example 10.3
!             using Korobov/good lattice point scheme
!             hopefully located & scaled to capture most posterior mass
!
      implicit none
      real t1,t2,t3,rll,pstmx,f,fnp,rlgpst
      real, dimension(10) :: ss
      real, dimension(3) :: sc
      integer k1,k2,k3,j,np,g2,g3
!
!*  21  format(4f16.8)
  28  format(/'Korobov rule',i8,' points -- Normalization constant',  &
     &       e14.6/8x,'Posterior mean',10x,'Covariance Matrix')
  29  format('  theta(1)',f12.6,6x,f12.6/                             &
     &       '  theta(2)',f12.6,6x,2f12.6/                            &
     &       '  theta(3)',f12.6,6x,3f12.6)
!                   constants of good lattice point schemes
!   due to Hua & Wang (1981) and given also in K-T Fang & Y Wang (1994)
!    Number-theoretic Methods in Statistics, Chapman & Hall, New York
!
!*      data np/ 101/, g2/ 40/, g3/ 85/           ! smaller rules
!*      data np/ 5037/, g2/ 580/, g3/ 1997/
      data np/ 10007/, g2/ 544/, g3/ 5733/
!                   scale factors for components
      data sc/ 5., 4., 6. /
!
!                   output unit
      open( unit=6, file='chex103.out' )
!                   max of log posterior
      pstmx = 18.466364
!                   initialize sums to zero
      ss = 0.
      fnp = real(np)
      k1 = 0
      k2 = 0
      k3 = 0
!                   main loop through points
      do j = 1,np
!                   g1 = 1 -- Riemann sum
      k1 = k1 + 1
      k2 = k2 + g2
      if( k2 .gt. np ) k2 = k2 - np
      k3 = k3 + g3
      if( k3 .gt. np ) k3 = k3 - np
!                   avoid zero -- right-handed sum
      t1 = sc(1)*real(k1+1)/fnp
      t2 = sc(2)*real(k2+1)/fnp
      t3 = 4. + sc(3)*real(k3+1)/fnp
!                   actually getting computing -log(posterior)
      rll = rlgpst(t1,t2,t3) - pstmx
!*      WRITE(6,21) t1,t2,t3,rll
      f = 0.
!                   avoid underflow if too small
      if( rll .lt. 30. ) f = exp(-rll)
!                   equally weighted rule
      ss(1) = ss(1) + f
      ss(2) = ss(2) + f*t1
      ss(3) = ss(3) + f*t2
      ss(4) = ss(4) + f*t3
      ss(5) = ss(5) + f*t1*t1
      ss(6) = ss(6) + f*t1*t2
      ss(7) = ss(7) + f*t2*t2
      ss(8) = ss(8) + f*t1*t3
      ss(9) = ss(9) + f*t2*t3
      ss(10) = ss(10) + f*t3*t3
      end do  ! loop on j
!                   posterior moments 
      ss(2) = ss(2) / ss(1)
      ss(3) = ss(3) / ss(1)
      ss(4) = ss(4) / ss(1)
      ss(5) = ss(5) / ss(1) - ss(2)*ss(2)
      ss(6) = ss(6) / ss(1) - ss(2)*ss(3)
      ss(7) = ss(7) / ss(1) - ss(3)*ss(3)
      ss(8) = ss(8) / ss(1) - ss(2)*ss(4)
      ss(9) = ss(9) / ss(1) - ss(3)*ss(4)
      ss(10) = ss(10) / ss(1) - ss(4)*ss(4)
!                   fix normalization constant
      ss(1) = ss(1) * sc(1) * sc(2) * sc(3) * exp(-pstmx) / fnp
!                   write out results
      write(6,28) np,ss(1)
      write(6,29) ss(2),ss(5),ss(3),ss(6),ss(7),ss(4),ss(8),ss(9),ss(10)
      stop
      end program chex103
      real function rlgpst(t1,t2,t3)
!  - log( posterior ) for exchange variance components problem
!      data are 100*(rate-1.6) for USD/GBP rate in May 1997
!
!  t1  error/within variance component
!  t2  trader/tmt   variance component
!  t3  mean parameter
!
      implicit none
      real, intent(in) :: t1,t2,t3
      integer, parameter :: p = 5       ! number of groups
      integer, dimension(p) :: ni
      real, dimension(p) :: yb
      real w,a1,b1,a2,b2,b0,phi0,sf,sg
!
      integer i,nbig
!                   data for problem
!  p         number of groups
!  ni        number of obs in group i
!  nbig      sum of ni
!  w         within/error sum of squares
!  a1,b1     prior for error component is inversegamma( a1, b1 )
!  a2,b2     prior for treatment component is inversegamma( a2, b2 )
!  b0,phi0   prior for mean in Normal( b0, phi0 )
!  yb        mean for group i
      data ni/ 2, 1, 7, 3, 2 /
      data yb/ 6.8, 10.3, 5.9, 6.6, 8.5/
!
      w = 13.10     ! within ss
      nbig = 15
      a1 = 1.
      b1 = 1.
      a2 = 4.
      b2 = 4.
      b0 = 8.
      phi0 = 16.
!                   assign outlanding posterior (tiny) if out of range
      rlgpst = 1.e6
      if( t1 .le. 0.d0 ) return
      if( t2 .le. 0.d0 ) return
!                   first part of posterior
      rlgpst = (a2 + p/2. + 1)*log( t2 ) + (a1 + nbig/2. + 1)*log( t1 )
      sg = 0.
      sf = 0.
      do i = 1,p
      sg = sg + log( ni(i)/t1 + 1./t2 )
      sf = sf + (t3-yb(i))*(t3-yb(i))/(t2+t1/ni(i))
      end do  ! loop on i
      rlgpst = rlgpst + b2/t2 + (b1+w/2.)/t1                &
     &            + (t3-b0)*(t3-b0)/(2.*phi0)               &
     &            + sg/2. + sf/2.
!
!                  using - log( posterior)  *** reminder
      rlgpst = rlgpst
      return
      end function rlgpst
! *** end of filename chex103.f95       *********************