# quad1.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 }
 # how many intervals?
    for ( kint in c(2,3,5,10,20) ) {
    print("number of intervals")
    print(kint)
 # trapeziod rule
    x <- (0:kint)/kint
    w <- c(1/2,rep(1,kint-1),1/2)
    px <- sapply(x,pstar)
    pt0 <- sum(px*w)
    pt1 <- sum(x*px*w)
    pt2 <- sum(x*x*px*w)
 # midpoint rule
    x <- (2*(1:kint)-1)/(2*kint)
    px <- sapply(x,pstar)
    pm0 <- sum(px)
    pm1 <- sum(x*px)
    pm2 <- sum(x*x*px)
 # simpson's rule
    ps0 <- (2*pm0 + pt0)/3
    ps1 <- (2*pm1 + pt1)/3
    ps2 <- (2*pm2 + pt2)/3
 # convert to get posterior means and variances
 # first trapeziod
    pt1 <- pt1/pt0               # posterior mean
    pt2 <- pt2/pt0 - pt1*pt1       # posterior variance
    pt0 <- pt0/kint
    print("trapezoid")
    cat("norm const",pt0,"mean",pt1,"var",pt2,"trapezoid","\n")
 # second midpoint
    pm1 <- pm1/pm0               # posterior mean
    pm2 <- pm2/pm0 - pm1*pm1       # posterior variance
    pm0 <- pm0/kint
    print("midpoint")
    cat("norm const",pm0,"mean",pm1,"var",pm2,"midpoint","\n")
 # third simpson rule
    ps1 <- ps1/ps0               # posterior mean
    ps2 <- ps2/ps0 - ps1*ps1       # posterior variance
    ps0 <- ps0/kint
    print("Simpson")
    cat("norm const",ps0,"mean",ps1,"var",ps2,"Simpson","\n")
                         }    # end loop
 # done
  rm(list=ls())
  q()