#  chex102s.r                 Oct 2007, Nov 2008, Oct 2009
 #
 # Bayesian analysis of Extended Example ** new data
 #  integration using Simpson's rule
 #
  # 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
  #
  #
  # Simpson's rule, 4 se's in each dimension
  #
    for (kint in c(3,5,10,20)) {          # begin loop
     t1 <- tm[1] + se[1]*(4/kint)*(-kint:kint)    # first component
     t2 <- tm[2] + se[2]*(4/kint)*(-kint:kint)    # second component
     np <- 2*kint + 1                             # number of points
     w <- c(1,rep(c(4,2),kint-1),4,1)/(6*kint)    # Simpson weights
  #               first row    second row  
     xx <- matrix(c(rep(t1,np),rep(t2,each=np)),  # lattice of points
  #            rows 2, np*np columns
                    2, np*np, byrow=TRUE)         #  to vector of pts
     ww <- as.vector(w %o% w)                     # weights to match
     px <- exp(-apply(xx,2,ull) - lpstar0)        # eval & normalize
     p0 <- sum(px*ww)                             # add for norm const
     print("number of intervals, norm const, mean, cov")
     print(c(kint,p0*(2*4)*(2*4)*prod(se)))  
     pmean <- apply(t(xx)*ww*px,2,sum)/p0            # post mean vector
     print(pmean) 
     pcov <- (xx %*% (t(xx)*(ww*px/p0))) - outer(pmean,pmean)
     print(pcov)                                  # post cov matrix
                                  }          # end loop  
  #
  # Simpson's rule, use covariance to reorient
  #
    for (kint in c(3,5,10,20,100)) {          # begin loop
     z <- (4/kint)*(-kint:kint)                   # (-4,4) in z space
     np <- 2*kint + 1                             # number of points
     w <- c(1,rep(c(4,2),kint-1),4,1)/(6*kint)    # Simpson weights
  #               first z1    second z2
     zz <- matrix(c(rep(z,np),rep(z,each=np)),
       2, np*np, byrow=TRUE)                      # col vector of pts
     ww <- as.vector(w %o% w)                     # weights to match
     xx <- matrix(tm,2,np**2) + backsolve(t(L),zz)# Simpson absiccas
     px <- exp(-apply(xx,2,ull) - lpstar0)        # eval & normalize
     p0 <- sum(px*ww)                             # add for norm const
     print("number of intervals, norm const, mean, cov")
     print(c(kint,p0*(2*4)*(2*4)/prod(diag(L))))
     pmean <- apply(t(xx)*ww*px,2,sum)/p0         # post mean vector
     print(pmean) 
     pcov <- xx%*%(t(xx)*(ww*px/p0))-outer(pmean,pmean)
     print(pcov)                                  # post cov matrix
                                  }          # end loop  
  #  done
    rm(list=ls())
  q()