R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), 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 a HTML browser interface to help.
Type 'q()' to quit R.

>  #  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
[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.0021986699
>     se <- sqrt(diag(cov))
>     se
[1] 0.06941788 0.04688998
>   #
>   # 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  
[1] 0.9816849
[1] 0.1664055
[1] "number of points, norm const, mean, cov"
[1] 200.00000000   0.01972402
[1] 0.2714114 0.1239343
              [,1]          [,2]
[1,]  0.0043253176 -0.0005712711
[2,] -0.0005712711  0.0020557529
[1] 0.9526774
[1] 0.1280419
[1] "number of points, norm const, mean, cov"
[1] 1.000000e+03 1.914120e-02
[1] 0.2731846 0.1239405
              [,1]          [,2]
[1,]  0.0044049401 -0.0004379248
[2,] -0.0004379248  0.0018661081
[1] 0.962362
[1] 0.2555061
[1] "number of points, norm const, mean, cov"
[1] 4.000000e+04 1.933578e-02
[1] 0.2718387 0.1251584
              [,1]          [,2]
[1,]  0.0044810984 -0.0006673045
[2,] -0.0006673045  0.0022172711
>   #  done
>     rm(list=ls())
>   q()