# 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
    xm <- pmode$maximum
    lpstar0 <- pmode$objective
    d2 <- (lpstar(xm+1e-5) - 2*lpstar0 + lpstar(xm-1e-5))/(1e-10)
    d2
 # 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
    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
 # done
  rm(list=ls())
  q()