#  chex102g.r                 Oct 2007, Oct, Nov 2008, Oct 2009
 #
 # Bayesian analysis of Extended Example ** new data
 #  integration using Gauss-Hermite quadrature
 #
  # 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
  #
   qgaus <-function(n,type=3,alph=0) {
 #            default=Hermite, alph parameter in Laguerre
        nm1 <- 1:(n-1)        # vector one smaller
 #            main branch
        if( type == 1 ) {     # Legendre
          alpha <- nm1/sqrt(4*nm1*nm1-1)
          beta <- rep(0,n)    }
        if( type == 2 ) {     # shifted Legendre
          alpha <- (nm1/2)/sqrt(4*nm1*nm1-1)
          beta <- rep(1/2,n)  }
        if( type == 3 ) {     # Hermite
          alpha <- sqrt(nm1/2)
          beta <- rep(0,n)    }
        if( type == 4 ) {     # modified Hermite
          alpha <- sqrt(nm1)
          beta <- rep(0,n)    }
        if( type == 5 ) {     # Laguerre, alph = parameter
          alpha <- sqrt(nm1*(nm1+alph))
          beta <- 2*(1:n)+alph-1  }
 # make matrix
   A <- rbind(cbind(rep(0,n-1),diag(alpha,n-1)),rep(0,n))
   A <- A + t(A) + diag(beta,n)
   e <- eigen(A,symmetric=T)  # eigenvalues and vectors
   abscissas <- e$values
   weights <- (e$vectors[1,])*(e$vectors[1,])
   qgaus <- list(abscissas,weights) 
                       }      # end of function qgaus
  #
  # Gauss-Hermite rule, use covariance to reorient
  #
    for (np in c(7,11,21)) {            # begin loop
     zw <- qgaus(np,type=4)
     z <- as.vector(zw[[1]])            # abscissas
     w <- as.vector(zw[[2]])            # modified Gauss-Hermite weights
  #               row1 is z1, row2 is z2 
     zz <- matrix(c(rep(z,np),rep(z,each=np)),
       2, np*np, byrow=TRUE)                      # vector of pts
     zz2 <- apply(zz*zz,2,sum)/2                  # z'z/2
     ww <- as.vector(w %o% w)                     # weights to match
     xx <- matrix(tm,2,np**2) + backsolve(t(L),zz)# G-H absiccas
     px <- exp(-apply(xx,2,ull) - lpstar0 + zz2)  # eval & normalize
     p0 <- sum(px*ww)                             # add for norm const
     print("number of intervals, norm const, mean, cov")
     print(c(np,p0*(2*pi)/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()