#  chex102n.r                 October 2007, 2008
 #   *** revised version ***
 #
 # Bayesian analysis of Extended Example ** new data
 #  integration using importance sampling from normal approximation
 #
  # define shortened function
    ull <- function(t) { # un-log-likelihood for extended example
        t1 <- t[1]         #  without gradient and Hessian
        t2 <- t[2]
        n <- c(1,10,4,9)   # note new data
        p <- c( 2*t1*t2, t1*(2-t1-2*t2), t2*(2-t2-2*t1), (1-t1-t2)**2 )
        if( (min(p)>0)&(max(p)<1) ) ull <- -crossprod( n, log(p) )
              else ull <- 1e100 }
  #
 # starting value
    t <- c(.263,.074)
 # use nlm this time instead of Newton
   mll <- nlm( ull, t, hessian=T, typsize=c(1,1), print.level=0 )
  # have posterior mode, now do two dimensional Simpson
   L <- t(chol(mll$hessian))  # hessian now pos def, not neg
   tm <- mll$estimate
   tm
   lpstar0 <- -mll$minimum    # reverse sign
   lpstar0
  #
  # now get inverse of -Hessian for std errors
    cov <-backsolve( t(L), forwardsolve(L,diag(1,nrow(L))) ) 
    cov
    se <- sqrt(diag(cov))
    se
  #
  # set seed
    set.seed(5151917)
  #
    for (np in c(200,1000,40000)) {          # begin loop
     zz <- matrix(rnorm(2*np),2,np)               # random abscissas
     zz2 <- apply(zz*zz,2,sum)/2                  # z'z/2
     xx <- matrix(tm,2,np) + backsolve(t(L),zz)   # tfm to t's
     px <- exp(-apply(xx,2,ull) - lpstar0 + zz2)  # eval & normalize
     p0 <- sum(px)                                # add for norm const
     print(mean(px))
     print(var(px))
     print("number of points, norm const, mean, cov")
     print(c(np,p0*(2*pi)/(np*prod(diag(L)))))    
     pmean <- apply(t(xx)*px,2,sum)/p0            # post mean vector
     print(pmean) 
     pcov <- xx%*%(t(xx)*(px/p0))-outer(pmean,pmean)
     print(pcov)                                  # post cov matrix
                                  }          # end loop  
  #  done
    rm(list=ls())
  q()