[R] multivariate integral with ADAPT when the parameter is close to boundary

[Muhtar Osman]

Thanks so much for the comments.
Of course, integrating
"dbeta(theta[1],0.005,0.005)*dbeta(theta[2],0.005,0.005)" over the
whole range is not something I am really interested in. I am just
concerned about the results I got for other integrands with extreme
priors using ADAPT function. It is true that some of the functions
that I want to integrate are products. But then there should be a
similar problem with 1-D function INTEGRATE for the extreme integrands
in my setting, that is, lower and upper limits have to be set very
very close to the boundaries. However, the closer I set those limits
to the boundaries, the longer it takes to run. I was just wondering if
anybody had encountered the similar problem before and tried some way
to get around such as transforming the variables.
Again, thank you for the comments.

On Sun, Oct 19, 2008 at 1:51 AM, Prof Brian Ripley
<ripley at stats.ox.ac.uk> wrote:
> 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
>