Hi, thanks for the pointer to cubature (which i had probably dismissed too quickly). Your tests with "f" should not work: the domain of f(x,.) is restricted to positive reals, but this domain of integration is then transformed in mixedrule() to map the semi-infinite range to a more reasonable domain (for example [-4,4]).
Thanks, baptiste On 21 September 2010 14:16, David Winsemius <dwinsem...@comcast.net> wrote: > Baptiste; > > You should see if this meets your requirements: > > help(adaptIntegrate, package=cubature) > > (I got errors when I ran the code and NaN's when I looked at the output of > test function, "f".) > >> vAverage(mixedrule, -4, 4, 0.0, 1, 20, f) - c(pi, pi/2, 2*pi) > Error: object 'mixedrule' not found > >> f(-4,0) > [1] NaN NaN NaN > Warning message: > In sqrt(x) : NaNs produced >> f(-3.9,0) > [1] NaN NaN NaN > Warning message: > In sqrt(x) : NaNs produced >> f(-3.9,0.1) > [1] NaN NaN NaN > Warning message: > In sqrt(x) : NaNs produced > -- > David > On Sep 21, 2010, at 4:11 AM, baptiste auguie wrote: > >> Dear list, >> >> I'm seeking some advice regarding a particular numerical integration I >> wish to perform. >> >> The integrand f takes two real arguments x and y and returns a vector >> of constant length N. The range of integration is [0, infty) for x and >> [a,b] (finite) for y. Since the integrand has values in R^N I did not >> find a built-in function to perform numerical quadrature, so I wrote >> my own after some inspiration from a post in R-help, >> >> library(statmod) >> >> ## performs 2D numerical integration >> ## using Gauss-Legendre quadrature >> ## with N points for x and y >> >> vAverage <- function(f, a1,b1, a2,b2, N=5, ...){ >> >> GL <- gauss.quad(N) >> nodes <- GL$nodes >> weights <- GL$weights >> >> C2 <- (b2 - a2) / 2 >> D2 <- (b2 + a2) / 2 >> y <- nodes*C2 + D2 >> >> C1 <- (b1 - a1) / 2 >> D1 <- (b1 + a1) / 2 >> x <- nodes*C1 + D1 >> >> value <- 0 >> for (ii in seq_along(x)){ >> tmp <- 0 >> for (jj in seq_along(y)){ >> tmp <- tmp + C1 * weights[jj] * f(x[jj], y[ii], ...) >> } >> value <- value + C2 * weights[ii] * tmp >> } >> value >> } >> >> ## test function, the result is pi for y=1 >> f <- function(x, y) { >> res <- 1 / (sqrt(x)*(1+x)) >> c(res, res/2, 2*res) >> } >> >> ## Transformation rule from Numerical Recipes >> ## to deal with the [0, infty) range of x >> >> mixedrule <- function(x, y, f, ...) >> { >> t <- exp(pi*sinh(x)) >> dtdx <- t*(pi*cosh(x)) >> f(t, y, ...)*dtdx >> } >> >> >> vAverage(mixedrule, -4, 4, 0.0, 1, 20, f) - c(pi, pi/2, 2*pi) >> ## -3.535056e-06 -1.767528e-06 -7.070112e-06 >> >> >> So it seems to work. I wonder though if I may have missed an easier >> (and more reliable) way to perform such integration using base >> functions or an add-on package that I may have overlooked. >> >> Best regards, >> >> baptiste >> >> ______________________________________________ >> R-help@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide >> http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. > > David Winsemius, MD > West Hartford, CT > > ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.