83 lines
2.1 KiB
R
83 lines
2.1 KiB
R
# Copyright (C) 1997-2005 The R Core Team
|
|
|
|
## Adaptive integration: Venables and Ripley pp. 105-110
|
|
## This is the basic integrator.
|
|
|
|
area <- function(f, a, b, ..., fa = f(a, ...), fb = f(b, ...), limit
|
|
= 10, eps = 1.e-5)
|
|
{
|
|
h <- b - a
|
|
d <- (a + b)/2
|
|
fd <- f(d, ...)
|
|
a1 <- ((fa + fb) * h)/2
|
|
a2 <- ((fa + 4 * fd + fb) * h)/6
|
|
if(abs(a1 - a2) < eps)
|
|
return(a2)
|
|
if(limit == 0) {
|
|
warning(paste("iteration limit reached near x = ",
|
|
d))
|
|
return(a2)
|
|
}
|
|
area(f, a, d, ..., fa = fa, fb = fd, limit = limit - 1,
|
|
eps = eps) + area(f, d, b, ..., fa = fd, fb =
|
|
fb, limit = limit - 1, eps = eps)
|
|
}
|
|
|
|
|
|
## The function to be integrated
|
|
|
|
fbeta <- function(x, alpha, beta)
|
|
{
|
|
x^(alpha - 1) * (1 - x)^(beta - 1)
|
|
}
|
|
|
|
|
|
## Compute the approximate integral, the exact integral and the error
|
|
|
|
b0 <- area(fbeta, 0, 1, alpha=3.5, beta=1.5)
|
|
b1 <- exp(lgamma(3.5) + lgamma(1.5) - lgamma(5))
|
|
c(b0, b1, b0-b1)
|
|
|
|
|
|
## Modify the function so that it records where it was evaluated
|
|
|
|
fbeta.tmp <- function (x, alpha, beta)
|
|
{
|
|
val <<- c(val, x)
|
|
x^(alpha - 1) * (1 - x)^(beta - 1)
|
|
}
|
|
|
|
|
|
## Recompute and plot the evaluation points.
|
|
|
|
val <- NULL
|
|
b0 <- area(fbeta.tmp, 0, 1, alpha=3.5, beta=1.5)
|
|
plot(val, fbeta(val, 3.5, 1.5), pch=0)
|
|
|
|
|
|
## Better programming style -- renaming the function will have no effect.
|
|
## The use of "Recall" as in V+R is VERY black magic. You can get the
|
|
## same effect transparently by supplying a wrapper function.
|
|
## This is the approved Abelson+Sussman method.
|
|
|
|
area <- function(f, a, b, ..., limit=10, eps=1e-5) {
|
|
area2 <- function(f, a, b, ..., fa = f(a, ...), fb = f(b, ...),
|
|
limit = limit, eps = eps) {
|
|
h <- b - a
|
|
d <- (a + b)/2
|
|
fd <- f(d, ...)
|
|
a1 <- ((fa + fb) * h)/2
|
|
a2 <- ((fa + 4 * fd + fb) * h)/6
|
|
if(abs(a1 - a2) < eps)
|
|
return(a2)
|
|
if(limit == 0) {
|
|
warning(paste("iteration limit reached near x =", d))
|
|
return(a2)
|
|
}
|
|
area2(f, a, d, ..., fa = fa, fb = fd, limit = limit - 1,
|
|
eps = eps) + area2(f, d, b, ..., fa = fd, fb =
|
|
fb, limit = limit - 1, eps = eps)
|
|
}
|
|
area2(f, a, b, ..., limit=limit, eps=eps)
|
|
}
|