[R] R numerical integration

[Hans W Borchers]

Hans W Borchers <hwborchers <at> googlemail.com> writes:

> 
> casperyc <casperyc <at> hotmail.co.uk> writes:
> 
> > Is there any other packages to do numerical integration other than the
> > default 'integrate'?
> > Basically, I am integrating:
> >
> >  integrate(function(x) dnorm(x,mu,sigma)/(1+exp(-x)),-Inf,Inf)$value
> >
> > The integration is ok provided sigma is >0.
> > However, when mu=-1.645074 and sigma=17535.26 It stopped working.
> > On the other hand, Maple gives me a value of 0.5005299403.
> 
> Using `integrate()` to integrate from -1e-8 to 1e-8 will give quite a correct
> result, while integrating from -1e-10 to 1e-10 will return 0.

Saturday morning... Well, of course i meant integrating from -1e8 to 1e8 and
from -1e10 to 1e10. The first one returns almost the correct result, while the
other returns 0. The same happens for `adaptIntegrate` in package cubature.

It shows that one cannot automatically set the limits very high. Therefore, 
transforming to a finite intervall is to be preferred. There are several way to
do that, depending also on the convergence rate of your function at infinity.

Hans Werner

> It seems more appropriate to transform the infinite into a finite interval.
> Function `quadinf` in package pracma does this automatically, applying the
> `integrate` routine to the finite interval [-1, 1].
> 
>     library(pracma)
>     quadinf(fun, -Inf, Inf, tol=1e-10)
>     # [1] 0.4999626
> 
> I am astonished to see the result from Maple as this does not appear to be
> correct --- Mathematica, for instance, returns 0.499963.
>