>  # quad2.r                October 2007, 2008
>  #
>  # simple one-dimensional integration
>  #
>  # define function to be integrated
>   pstar <- function(t) {
+      if ( t*(1-t) > 0 ) 
+         pstar <- (1-t)*(t**16)/( (-log(1-t))**10)
+         else
+         pstar <- 0 }
>  # log of function to be integrated
>   lpstar <- function(t) {
+      if ( t*(1-t) > 0 ) 
+         lpstar <- log(1-t) + 16*log(t) - 10*log( -log(1-t) )
+         else
+         lpstar <- -1e200 }     # close enough to -inf
>  # gauss quadrature code
>    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
>  #
>  # first find posterior mode
>     pmode <- optimize( lpstar,c(.05,.95),maximum=T)
>     pmode
$maximum
[1] 0.5271652

$objective
[1] -8.10283

>     xm <- pmode$maximum
>     lpstar0 <- pmode$objective
>     d2 <- (lpstar(xm+1e-5) - 2*lpstar0 + lpstar(xm-1e-5))/(1e-10)
>     d2
[1] -42.03599
>  # how many evaluations?
>  #     
>     for ( kint in c(2,3,5,10,20) ) {
+     cat("number of points",kint,"\n")
+  # first modified Hermite
+     qg <-qgaus(kint,4)
+     z <- qg[[1]]                        # standardized abscissas 
+     x <- xm + (qg[[1]])/sqrt(-d2)       # center and scale
+     w <- qg[[2]]
+     lpx <- sapply(x,lpstar) - lpstar0   # 'normalize'
+     px <- exp(lpx+z*z/2)*sqrt(6.2821853) # divide by normal
+     pt0 <- sum(px*w)
+     pt1 <- sum(x*px*w)
+     pt2 <- sum(x*x*px*w)
+     pt1 <- pt1/pt0               # posterior mean
+     pt2 <- pt2/pt0 - pt1*pt1       # posterior variance
+     pt0 <- pt0*exp(lpstar0)/sqrt(-d2) # return constants 
+     print("Gauss-Hermite")
+     cat("norm const",pt0,"mean",pt1,"var",pt2,"\n")
+                          }    # end loop
number of points 2 
[1] "Gauss-Hermite"
norm const 0.0001142720 mean 0.5273612 var 0.0237891 
number of points 3 
[1] "Gauss-Hermite"
norm const 0.0001084708 mean 0.5272444 var 0.02004276 
number of points 5 
[1] "Gauss-Hermite"
norm const 0.0001102353 mean 0.5272062 var 0.01908349 
number of points 10 
[1] "Gauss-Hermite"
norm const 0.0001097692 mean 0.5270053 var 0.01863429 
number of points 20 
[1] "Gauss-Hermite"
norm const 0.0001097960 mean 0.5270514 var 0.01871944 
>     for ( kint in c(2,3,5,10,20) ) {
+     cat("number of points",kint,"\n")
+  # how many evaluations?
+  #     
+  # next shifted Legendre so interval is (0,1)
+     qg <-qgaus(kint,2) 
+     x <- qg[[1]]       # no center and scale
+     w <- qg[[2]]
+     px <- sapply(x,pstar)
+     pt0 <- sum(px*w)
+     pt1 <- sum(x*px*w)
+     pt2 <- sum(x*x*px*w)
+     pt1 <- pt1/pt0               # posterior mean
+     pt2 <- pt2/pt0 - pt1*pt1       # posterior variance
+     print("shifted Legendre")
+     cat("norm const",pt0,"mean",pt1,"var",pt2,"\n")
+                          }    # end loop
number of points 2 
[1] "shifted Legendre"
norm const 3.973057e-05 mean 0.6292448 var 0.0666291 
number of points 3 
[1] "shifted Legendre"
norm const 0.0001346098 mean 0.5046054 var 0.002385241 
number of points 5 
[1] "shifted Legendre"
norm const 0.0001105693 mean 0.5253538 var 0.01631614 
number of points 10 
[1] "shifted Legendre"
norm const 0.000109807 mean 0.5270573 var 0.01871998 
number of points 20 
[1] "shifted Legendre"
norm const 0.000109807 mean 0.5270573 var 0.01871998 
>  # done
>   rm(list=ls())
>   q()