R help - Oct 2008 - multivariate integral with ADAPT when the parameter is close to boundary

Dear All,

There is one problem I encountered when I used ADAPT to compute some
2-D integral w.r.t beta density.

For example, when I try to run the following comments:

fun2<-function(theta){return(dbeta(theta[1],0.005,0.005)*dbeta(theta[2],0.005,0.005))}
int.fun2<-adapt(ndim=2,lo = c(0,0), up = c(1,1),functn = fun2,eps = 1e-4)

It seems it will take very long time to run. Acturally, I stopped the
program after it was running for like 20 minutes.

I thought this might be due to the inclusion of the lower and upper in
to the integral computation, so I tried to change the lower and upper
limits:

fun2<-function(theta){return(dbeta(theta[1],0.005,0.005)*dbeta(theta[2],0.005,0.005))}
int.fun2<-adapt(ndim=2,lo = c(0.0001,0.0001), up c(0.9999,0.9999),functn =
fun2,eps = 1e-4)

It only took few seconds to run, but it gave me the wrong result:
int.fun2= 0.00202210665273673, whereas the correct result should be int.fun2=1.

I guess the reason for this is beta(0.005,0.005) has very high density
close to the boundary (theta=0).
So even letting  "lo = c(0.0001,0.0001)" will cause some loss of
probability mass in the integral computation.

I was wondering if anybody has encountered the similar problem before.
Any comments are appreciated.
Thanks.

Muhtar Osman
Dept.of Stats
NCSU

1) That integrand is a product, so you can do this a product of integrals, 
and do those analytically.

2) Do you have any idea how extreme beta(0.005, 0.005) is?  See the 
comment in the help for integrate:

      Like all numerical integration routines, these evaluate the
      function on a finite set of points.  If the function is
      approximately constant (in particular, zero) over nearly all its
      range it is possible that the result and error estimate may be
      seriously wrong.

delta <- 1e-4
x <- seq(delta, 1-delta, delta)
plot(x, dbeta(x, 0.005, 0.005), type="l")
pbeta(0.9999, 0.005, 0.005) - pbeta(0.0001, 0.005, 0.005)

so 95% of the mass is outside the limits you set.

On Sun, 19 Oct 2008, Muhtar Osman wrote:
> Dear All,
>
> There is one problem I encountered when I used ADAPT to compute some
> 2-D integral w.r.t beta density.
>
> For example, when I try to run the following comments:
>
>
fun2<-function(theta){return(dbeta(theta[1],0.005,0.005)*dbeta(theta[2],0.005,0.005))}
> int.fun2<-adapt(ndim=2,lo = c(0,0), up = c(1,1),functn = fun2,eps =
1e-4)
>
> It seems it will take very long time to run. Acturally, I stopped the
> program after it was running for like 20 minutes.
>
> I thought this might be due to the inclusion of the lower and upper in
> to the integral computation, so I tried to change the lower and upper
> limits:
>
>
fun2<-function(theta){return(dbeta(theta[1],0.005,0.005)*dbeta(theta[2],0.005,0.005))}
> int.fun2<-adapt(ndim=2,lo = c(0.0001,0.0001), up >
c(0.9999,0.9999),functn = fun2,eps = 1e-4)
>
> It only took few seconds to run, but it gave me the wrong result:
> int.fun2= 0.00202210665273673, whereas the correct result should be
int.fun2=1.
No, that's the correct answer for the problem you set.
>
> I guess the reason for this is beta(0.005,0.005) has very high density
> close to the boundary (theta=0).
> So even letting  "lo = c(0.0001,0.0001)" will cause some loss of
> probability mass in the integral computation.
>
> I was wondering if anybody has encountered the similar problem before.
> Any comments are appreciated.
> Thanks.
>
> Muhtar Osman
> Dept.of Stats
> NCSU
>
> ______________________________________________
> R-help at 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.
>
-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595