>Does anybody know whether S-Plus has a function for the Gauss-Hermite
>quadrature formula which is an approximation of an integral?
Here are a number of functions.
--------------------------------------------------
# Tim Hesterberg, last modified 3 Nov 03.
#In this file
#Integrate: 10-point Gaussian integration
#Integrate: also allows integration with respect to Gaussian density
#Integrate2: front end, call Integrate for each range defined by cuts
#Integrate2g: like Integrate2, for Gaussian density
#Integrate6: Like first Integrate, but 4point Gaussian
#CumIntegrate: 10-point Gaussian integration, cumulative
Integrate <- function(f, a = 0, b = 1, nintervals = 5, ...)
{
# Integrate f from a to b
# using 10-point Gaussian on each of nintervals subintervals
# Assume that f can take a vector argument.
ends <- seq(a, b, len = nintervals + 1)
h <- (b - a)/(2 * nintervals) #half-width of a single subinterval
mids <- (ends[-1] + ends[1:nintervals])/2
x <- c(0.14887433898163099, 0.43339539412924699, 0.67940956829902399,
0.86506336668898498, 0.97390652851717197)
wt <- c(0.29552422471475298, 0.26926671930999602, 0.21908636251598201,
0.149451349150581, 0.066671344308687999)
xvalues <- outer(h * c( - x, x), mids, "+")
h * sum(wt * f(c(xvalues), ...))
}
Integrate <- function(f, a = 0, b = 1, nintervals = 5, ..., normal=F)
{
# Integrate \int_a^b f(x) dx
# If normal==T, Integrate \int_qnorm(a)^qnorm(b) f(x) \phi(x) dx
# 10-point Gaussian on each of nintervals subintervals
# Assume that f can take a vector argument
ends <- seq(a, b, len = nintervals + 1)
h <- (b - a)/(2 * nintervals) #half-width of a single subinterval
mids <- (ends[-1] + ends[1:nintervals])/2
x <- c(0.14887433898163099, 0.43339539412924699, 0.67940956829902399,
0.86506336668898498, 0.97390652851717197)
wt <- c(0.29552422471475298, 0.26926671930999602, 0.21908636251598201,
0.149451349150581, 0.066671344308687999)
xvalues <- outer(h * c( - x, x), mids, "+")
if(normal) xvalues <- qnorm( xvalues )
h * sum(rep(wt,length=length(xvalues)) * f(c(xvalues), ...))
}
Integrate2 <- function(f, cuts, nintervals = rep(5,length(cuts)-1), ...,
normal=F)
{
# This is a front end for Integrate(), call for each range defined by cuts.
Integrate2sum <- 0
for(i in 1:(length(cuts)-1)) Integrate2sum <- Integrate2sum +
Integrate( f, cuts[i], cuts[i+1], nintervals[i], ..., normal=normal)
Integrate2sum
}
#Experimental, based on first version of Integrate, return cumulative integral
CumIntegrate <- function(f, a = 0, b = 1, nintervals = 5, ...){
# 10-point Gaussian on each of nintervals subintervals
# Assume that f can take a vector argument
ends <- seq(a, b, len = nintervals + 1)
h <- (b - a)/(2 * nintervals) #half-width of a single subinterval
mids <- (ends[-1] + ends[1:nintervals])/2
x <- c(0.14887433898163099, 0.43339539412924699, 0.67940956829902399,
0.86506336668898498, 0.97390652851717197)
wt <- c(0.29552422471475298, 0.26926671930999602, 0.21908636251598201,
0.149451349150581, 0.066671344308687999)
xvalues <- outer(h * c( - x, x), mids, "+")
list( x = ends,
y = c(0,h*cumsum(colSums( matrix( wt*f(c(xvalues), ...), 10)))))
}
Integrate2g <- function(f, xcuts, nintervals = rep(5,length(ucuts)-1),
ucuts=pnorm(xcuts),...){
# Front end: call Integrate for each subinterval defined by cuts.
# Do integrals of the form \int f(x) phi(x) dx
# Supply either xcuts or ucuts (u = pnorm(x))
Integrate2sum <- 0
for(i in 1:(length(ucuts)-1)) Integrate2sum <- Integrate2sum +
Integrate( f, ucuts[i], ucuts[i+1], nintervals[i], ..., normal=T)
Integrate2sum
}
Integrate6 <- function(f, a = 0, b = 1, nintervals = 5, ...){
# 4-point Gaussian on each of nintervals subintervals
# Assume that f can take a vector argument
ends <- seq(a, b, len = nintervals + 1)
h <- (b - a)/(2 * nintervals) #half-width of a single subinterval
mids <- (ends[-1] + ends[1:nintervals])/2
x <- c(0.339981043584856, 0.861136311594053)
wt <- c(0.652145154862546, 0.347854845137454)
xvalues <- outer(h * c( - x, x), mids, "+")
h * sum(wt * f(c(xvalues), ...))
}
|