Double Integration

[Brian Beckage <bb2@duke.edu>]

Hi,

I need to integrate the following function in two dimensions:

fct<-function(x, y){ (a^-b)*b*exp(s*x-(((x^2 + y^2)^0.5)/a)^b)*(x^2 +
y^2)^(0.5*(b-1))}

over x (0, Inf) and y (-Inf, Inf).  a,b, and s are fixed for any particular
call to this function.  For example if 

a<-1.1;  b<-2.1;  s<-1.8;

then the value for this integral equals 20.6541 (answer from Mathematica).
I have tried several different functions from previous S-news queries or
from StatLib (see below), but they either do not allow Inf or -Inf as
limits or give answers that are different from Mathematica. (Mathematica
appears to be giving correct answer based on comparisons across functions
using finite limits.)  In addition, using large but finite limits does not
appear to provide consistent results (i.e. the resulting answer depends on
the limits chosen, rather than asymptotically approaching some answer). 

Any suggestions on how to perform this integration within Splus would be
greatly appreciated.  I am using Splus 2000 R3 (Windows).

Brian

#############################################################################
############### Function 1 (Gives different answer than Mathematica; see
below)
# Sorry but I don't recall who contributed this particular function.

dfn <- function(IL, IU, OL, OU, f11)
{
 IL <<- IL      
 IU <<- IU      
 OL <<- OL      # y lower limit
 OU <<- OU      # y upper limit
 fi1 <<- f11
 fun <- function(y)
 unlist(lapply(y, function(y)
 {
  integrate(f = fi1, lower = IL, upper = IU, y = y, subdivisions=100)$
   integral
 }
 ))
 integrate(f = fun, lower = OL, upper = OU,subdivisions=100)$integral
}
########### Here is the call to the function.
a<-1.1;b<-2.1;s<-1.8;
dfn(0, Inf, -Inf, Inf, function(x, y){ (a^-b)*b*exp(s*x-(((x^2 +
y^2)^0.5)/a)^b)*(x^2 + y^2)^    (0.5*(b-1))})

> 31.23751 (compared to Mathematica's 20.6541).  
 
############################################################################
###
################################## Function 2 (Can't use Inf or -Inf as
limits)
# From Diego Kuonen 
MC.integrate.2D <- function(fct, low, upp, npoints, exact.value=NULL)
  {
  points.x<-runif(n=npoints, min=low[1], max=upp[1])
  points.y<-runif(n=npoints, min=low[2], max=upp[2])
  approx.tmp <- fct(points.x, points.y)
                # approx.tmp-points.x^2*points.y   yields vector of 0's
  V.tmp <- diff(c(low[1], upp[1])) * diff(c(low[2], upp[2]))
  approximate <- mean(approx.tmp) * V.tmp
  varapproximate<- var(approx.tmp) * V.tmp/approximate
  if(!is.null(exact.value)) {
    cat("To achieve an error of 0.0001 you need at least",
         floor(abs(varapprox - exact.value^2)/(0.0001^2)),
         "points.\n")
    }
  return(approximate,varapproximate)
  }

############################################################################
###
################################## Function 3 (Can't use Inf or -Inf as
limits)
# from Brian Ripley's Integrat library

f<-function(x){ 
        (a^-b)*b*exp(s*x[1]-(((x[1]^2 + x[2]^2)^0.5)/a)^b)*(x[1]^2 +
x[2]^2)^(0.5*(b-1))}
dcuhre(f,low=c(0,-10),hi=c(10,10))

# which seems to agree with Mathematica for these finite limits but can't
be used 
with Inf or -Inf limits.

##############################################