R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
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 an HTML browser interface to help.
Type 'q()' to quit R.

>  # 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
[1] "number of intervals"
[1] 2
[1] "trapezoid"
norm const 0.0001490066 mean 0.5 var 0 trapezoid 
[1] "midpoint"
norm const 7.027747e-05 mean 0.5900083 var 0.0543985 midpoint 
[1] "Simpson"
norm const 9.65205e-05 mean 0.5436906 var 0.02842901 Simpson 
[1] "number of intervals"
[1] 3
[1] "trapezoid"
norm const 0.0001090315 mean 0.5352534 var 0.02653498 trapezoid 
[1] "midpoint"
norm const 0.0001105749 mean 0.5192481 var 0.01092115 midpoint 
[1] "Simpson"
norm const 0.0001100604 mean 0.5245333 var 0.01613376 Simpson 
[1] "number of intervals"
[1] 5
[1] "trapezoid"
norm const 0.0001097788 mean 0.5274963 var 0.01877648 trapezoid 
[1] "midpoint"
norm const 0.0001098368 mean 0.5266331 var 0.01866792 midpoint 
[1] "Simpson"
norm const 0.0001098175 mean 0.5269208 var 0.01870426 Simpson 
[1] "number of intervals"
[1] 10
[1] "trapezoid"
norm const 0.0001098078 mean 0.5270646 var 0.01872237 trapezoid 
[1] "midpoint"
norm const 0.0001098063 mean 0.5270506 var 0.01871787 midpoint 
[1] "Simpson"
norm const 0.0001098068 mean 0.5270553 var 0.01871937 Simpson 
[1] "number of intervals"
[1] 20
[1] "trapezoid"
norm const 0.0001098071 mean 0.5270576 var 0.01872012 trapezoid 
[1] "midpoint"
norm const 0.0001098069 mean 0.527057 var 0.01871986 midpoint 
[1] "Simpson"
norm const 0.0001098070 mean 0.5270572 var 0.01871995 Simpson 
>  # done
>   rm(list=ls())
>   q()