R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

>  #  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
[1] 0.2657811 0.1109792
>    lpstar0 <- -mll$minimum    # reverse sign
>    lpstar0
[1] -28.02714
>   #
>   # now get inverse of -Hessian for std errors
>     cov <-backsolve( t(L), forwardsolve(L,diag(1,nrow(L))) ) 
>     cov
              [,1]          [,2]
[1,]  0.0048188423 -0.0006078584
[2,] -0.0006078584  0.0021986700
>     se <- sqrt(diag(cov))
>     se
[1] 0.06941788 0.04688998
>   #
>    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  
[1] "number of intervals, norm const, mean, cov"
[1] 7.00000000 0.01937902
[1] 0.2720017 0.1244538
             [,1]         [,2]
[1,]  0.004485033 -0.000662668
[2,] -0.000662668  0.002293747
[1] "number of intervals, norm const, mean, cov"
[1] 11.00000000  0.01928117
[1] 0.2717955 0.1252678
             [,1]         [,2]
[1,]  0.004472892 -0.000639844
[2,] -0.000639844  0.002222759
[1] "number of intervals, norm const, mean, cov"
[1] 21.00000000  0.01929588
[1] 0.2718347 0.1251295
              [,1]          [,2]
[1,]  0.0044754230 -0.0006456731
[2,] -0.0006456731  0.0022423093
>   #  done
>     rm(list=ls())
>   q()