>  #  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
[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
>   #
>   #
>   # 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  
[1] "number of intervals, norm const, mean, cov"
[1] 3.00000000 0.02136211
[1] 0.2723867 0.1233853
              [,1]          [,2]
[1,]  0.0033561413 -0.0003944579
[2,] -0.0003944579  0.0018037765
[1] "number of intervals, norm const, mean, cov"
[1] 5.00000000 0.01908056
[1] 0.2719069 0.1257032
              [,1]          [,2]
[1,]  0.0044369488 -0.0006068924
[2,] -0.0006068924  0.0021192348
[1] "number of intervals, norm const, mean, cov"
[1] 10.00000000  0.01925860
[1] 0.2718876 0.1248082
              [,1]          [,2]
[1,]  0.0044611742 -0.0006256599
[2,] -0.0006256599  0.0021794746
[1] "number of intervals, norm const, mean, cov"
[1] 20.0000000  0.0192575
[1] 0.2718878 0.1248090
              [,1]          [,2]
[1,]  0.0044612220 -0.0006259058
[2,] -0.0006259058  0.0021803585
>   #
>   # 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  
[1] "number of intervals, norm const, mean, cov"
[1] 3.00000000 0.02140289
[1] 0.2712248 0.1239636
              [,1]          [,2]
[1,]  0.0034342207 -0.0005493519
[2,] -0.0005493519  0.0018505786
[1] "number of intervals, norm const, mean, cov"
[1] 5.00000000 0.01909049
[1] 0.2718919 0.1257385
              [,1]          [,2]
[1,]  0.0044156112 -0.0006036965
[2,] -0.0006036965  0.0021216190
[1] "number of intervals, norm const, mean, cov"
[1] 10.00000000  0.01925864
[1] 0.2718888 0.1248057
              [,1]          [,2]
[1,]  0.0044614818 -0.0006263064
[2,] -0.0006263064  0.0021795751
[1] "number of intervals, norm const, mean, cov"
[1] 20.00000000  0.01925755
[1] 0.2718888 0.1248066
              [,1]          [,2]
[1,]  0.0044615928 -0.0006265664
[2,] -0.0006265664  0.0021804627
[1] "number of intervals, norm const, mean, cov"
[1] 100.00000000   0.01925749
[1] 0.2718887 0.1248070
             [,1]         [,2]
[1,]  0.004461590 -0.000626556
[2,] -0.000626556  0.002180425
>   #  done
>     rm(list=ls())
>   q()