# chex103.s                         J F Monahan, June 2001
 # integrate variance components posterior using Korobov
 rlgpst <- function(t) {
   p <- 5                              # number of groups
   ni <- c(2, 1, 7, 3, 2)              # group sample sizes
   yb <- c( 6.8, 10.3, 5.9, 6.6, 8.5)  # sample mean vector
   w <- 13.10                          # within ss
   a1 <- 1; b1 <- 1;                   # prior on error, invg(a1,b1)
   a2 <- 4; b2 <- 4;                   # prior on group, invg(a2,b2)
   b0 <- 8; phi0 <- 16                 # prior on mean, N(b0,phi0)
   nbig <- sum(ni)
   pst <- 1000000
   if( min( t[1], t[2] ) > 0 ) {
   pst <- (a2+p/2+1)*log(t[2]) + (a1+nbig/2+1)*log(t[1])
   sg <- sum( log( ni/t[1] + 1/t[2] ) )
   sf <- sum( (t[3]-yb)*(t[3]-yb)/(t[2]+t[1]/ni) )
    pst <- pst + b2/t[2] + (b1+w/2)/t[1] + 
                (t[3]-b0)*(t[3]-b0)/(2*phi0) + sg/2 + sf/2
                        }          # end of if block
                         }         # end of function rlgpst
 # check mode
 pstmx <- rlgpst( c( 1.200429, .905818, 7.327349 ) )
 pstmx
 # Korobov rule
 np <- 10007                      # many, many points
 g <- c( 1, 544, 5733 )
#  np <- 101                      # debug with fewer
#  g <- c( 1, 40, 85)
  scale <- c( 5, 4, 6)
  shft <- runif(3)                # Cranley-Patterson randomization
 # initialize sums
  pnc <- 0
  pmn <- rep( 0, 3)
  pcv <- matrix( 0, 3, 3)
 # loop
 for( j in 1:np ) {
    k <- ( (j*g) %% np )
    t <- c( 0, 0, 4) + ((k+shft)/np)*scale  # locate & scale
    rll <- rlgpst(t) - pstmx
    if( rll < 30 ) {
      f <- exp(-rll)
      pnc <- pnc + f
      pmn <- pmn + f*t
      pcv <- pcv + f*outer(t,t)
                   }     # end of if block to avoid underflow
                  }      # end of loop on j
 # fix posteriors
    pmn <- pmn / pnc
    pcv <- pcv / pnc - outer(pmn,pmn)
    pnc <- pnc * prod(scale)*exp(-pstmx)/np
    "posterior mean"
    pmn
    "posterior covariance matrix"
    pcv
    "normalization constant"
    pnc
 rm(list=ls())