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.

[Previously saved workspace restored]

>  #  chex102z.r                 October 2007
>  # Bayesian analysis of Extended Example ** new data
>  #  integration using importance sampling from normal approximation
>  #  coded using loops
>  #
>   # 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
+   #                           initialize
+      p0 <- 0
+      pmean <- rep(0,2)
+      pcov <- matrix(0,2,2)
+   #                           loop this time
+      for ( i in (1:np)) {     # begin loop
+      z <- rnorm(2)            # abscissas
+      x <- tm + backsolve(t(L),z)
+      px <- exp(-ull(x) - lpstar0 + as.real(crossprod(z))/2)    # normalize
+      p0 <- p0 + px
+      pmean <- pmean + px*x
+      pcov <- pcov + outer(px*x,x)
+                         }     # end loop
+      print("number of points, norm const, mean, cov")
+      print(c(np,(p0*2*pi)/(np*prod(diag(L)))))
+      pmean <- pmean/p0
+      print(pmean) 
+      pcov <- pcov/as.real(p0) - outer(pmean,pmean)
+      print(pcov) 
+                                   }          # end loop  
[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] "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] "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()