Re: Double integration(summary)

[Diego Kuonen <Diego@Kuonen.com>]

"Dr. Takashi Kikuchi" wrote:
> 
> Thanks to Steven, Carlisle and Alec. This is a summary of responses. I hope
> that someone will find new information about double integration and add it
> to this summary.

First, let us start with the proposed Monte Carlo integration.  The analytical 
solution is 20*log(5)/6 = 5.364793, hence using the following function we see 
that we would need M=2823687451 to get an accuracy of 0.0001!

MC.integrate.2D <- function(fct, low=c(-1,1), upp=c(1,1),
                            npoints=100, 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 <- mean(approx.tmp)
  varapprox <- var(approx.tmp)
  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")
    }
  approx
  }

> test.fct <- function(x,y) 20/(y*x^2)
> MC.integrate.2D(test.fct, c(2,1), c(3,5), 100000, 20*log(5)/6)
To achieve an error of 0.0001 you need at least 2823687451 points.
[1] 1.337588

So, you may want to try

> MC.integrate.2D(test.fct, c(2,1), c(3,5), 2823687451, 20*log(5)/6)


Second, the real answer to your question lies in the use of adaptive quadrature
methods. Right now, there are two Fortran routines available,  DCUHRE 
and ADAPT, which can be dynamically loaded into S-Plus. The Fortran routine 
DCUHRE is implemented in the function dcuhre, which is contained
in the integrate2 library (lib.stat.cmu.edu/S/integrate2), and
ADAPT comes with the function adapt included in the library adapt 
(lib.stat.cmu.edu/S/adapt). For details I refer to the included 
references.

For your example a simple use would result in

> test.fct <- function(z) as.single(20/(z[2]*z[1]^2))
> adapt(ndim=2, lower=c(2,1), upper=c(3,5),
        minpts=100, maxpts=500,
        functn=test.fct,
        eps=1.0e-4)
Allocating a work space of size 73
$finest:
[1] 5.364787

$relerr:
[1] 2.773141e-06

$minpts:
[1] 125

$ifail:
[1] 0


> test.fct <- function(z) 20/(z[2]*z[1]^2)
> dcuhre(f=test.fct,
         low=c(2,1),
         hi=c(3,5),
       name="z",
       epsabs=0.0,
       epsrel=1.0e-4,
       mincall=0,
       maxcall=10000,
       key=0)

$result:
[1] 5.364793

$abserr:
[1] 8.438938e-05

$neval:
[1] 195


Note that adapt took, in a single run, 0.03 seconds and dchure 0.19 seconds 
in CPU time on a Sun SPARC Ultra 60 workstation with 1Gb RAM using 
S-Plus 3.4 Release 1.

Hope that's what you wanted.

Unfortunately, traditional adaptive quadrature methods have been almost
forgotten in the recent rush to MCMC methods. 

Greets

  Diego Kuonen
 
-- 
Diego DOT Kuonen AT epfl DOT ch      diego AT kuonen DOT com
http://stat.kuonen.com                 http://www.Statoo.com 
`If you can imagine it, you can achieve it; if you can dream 
it, you can become it.'              Powered by Open Source!