R numerical integration

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

R numerical integration

casperyc
Hi all,

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.

It is an important line of the coding that I am doing and I am looking for some package that is able to do numerical integration efficiently (fast and accurate to a tol=1e-4). I have tried 'cubature', which does not give me anything even after 10 minutes.

Thanks.

casper

######################
PhD candidate in Statistics
Big R Fan
Big LEGO Fan
Big sTaTs Fan
######################
Reply | Threaded
Open this post in threaded view
|

Re: R numerical integration

Hans W Borchers
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.
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.

> It is an important line of the coding that I am doing and I am looking for
> some package that is able to do numerical integration efficiently (fast and
> accurate to a tol=1e-4). I have tried 'cubature', which does not give me
> anything even after 10 minutes.
>
> Thanks. casper

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: 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.
>

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: R numerical integration

Petr Savicky
In reply to this post by casperyc
On Fri, Mar 23, 2012 at 01:27:57PM -0700, casperyc wrote:

> Hi all,
>
> 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.
>
> It is an important line of the coding that I am doing and I am looking for
> some package that is able to do numerical integration efficiently (fast and
> accurate to a tol=1e-4). I have tried 'cubature', which does not give me
> anything even after 10 minutes.

Hi.

Integrating with infinite limits is necessarily a heuristic. If you
want to bound the absolute error (not the relative error), then it is
sufficient to integrate over the interval

  [mu - bound*sigma, mu + bound*sigma]

where "bound" is such that the integral of the tails of
dnorm(, mu=1, sigma=1) outside of the interval is less than the
required absolute error. This follows from the fact that the
multiplier 1/(1+exp(-x)) is at most 1. For example, with bound = 10,
the absolute error due to limiting the interval is at most
2*pnorm(-10) = 1.523971e-23. The total error will contain this
and the error of the numerical algorithm.

  mu <- -1.645074
  sigma <- 17535.26
  integrate(function(x) dnorm(x,mu,sigma)/(1+exp(-x)), mu - 10*sigma, mu + 10*sigma)
  0.5 with absolute error < 3.7e-05

The exact result is at most 0.5. This may be seen by shifting the
integrated function so that the mean of the dnorm() becomes 0.
We get the following function f_1(x), which has the same integral
as the original function

  f_1(x) = dnorm(x, 0, sigma)/(1+exp(-x-mu))

Consider the function f_1(-x). It has also the same integral as the
original function. So, the function f_2(x) = f_1(x) + f_1(-x) has
twice the integral of the original function. Since dnorm(-x, 0, sd) =
dnorm(x, 0, sd), the function f_2(x) has the form

  f_2(x) = dnorm(x, 0, sigma)*(1/(1+exp(-x-mu)) + 1/(1+exp(x-mu)))

Since mu < 0, we have for all x

  1/(1+exp(-x-mu)) + 1/(1+exp(x-mu)) <= 1/(1+exp(-x)) + 1/(1+exp(x)) = 1

Consequently, the integral of f_2(x) is at most the integral of
dnorm(x, 0, sigma), which is 1. The integral of the original
function is one half of this, so it is at most 0.5. The result
0.5005299403 Maple is more than 0.5 only due to a numeric error.

Hope this helps.

Petr Savicky.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: R numerical integration

Peter Dalgaard-2

On Mar 24, 2012, at 09:46 , Petr Savicky wrote:

> Integrating with infinite limits is necessarily a heuristic.

...as is numerical integration in general. In the present case, the infinite limits are actually only half the problem. The integrate() function is usually quite good at dealing with infinite limits, but it sort of assumes that "most" of the functions behavior is played out relatively close to the origin. This is not really the case for a normal density with an sd of 17000+, so I tried rescaling:

> integrate(function(x) dnorm(x,mu/sigma)/(1+exp(-x*sigma)),-Inf,Inf)
0.4999626 with absolute error < 4.7e-05

Nice. I then thought I was smart and recentered it as well:

> integrate(function(x) dnorm(x)/(1+exp(-mu-x*sigma)),-Inf,Inf,rel.tol=1e-10)
0.5 with absolute error < 1.2e-11

and is sticks to 0.5 unless you reduce sigma to 100 or so. Oops.

As you probably noticed, the two integrands are the same function shifted left/right by mu/sigma, i.e. by roughly 1e-4, so the integral over the whole line is of course supposed to be the same.

This is where the second heuristic of numerical integration shows up: It is assumed that the integral is continuous, not just in the pure math sense but also in the pragmatic sense of being interpolatable from a handful of function values.

Now, the function 1/(1+exp(-x*sigma)) is, for large sigma, in principle continuous, but in practice pretty darn close to a Heaviside step function:

> f <- function(x)1/(1+exp(-x*sigma))
> f(0)
[1] 0.5
> f(0.001)
[1] 1
> f(-0.001)
[1] 2.424004e-08

this means that you have to be pretty lucky with the selection of knot points for integration, unless you give the algorithm a hint that something might be happening very quickly in certain regions. E.g.

> integrate(function(x) dnorm(x)/(1+exp(-mu-x*sigma)),-Inf,0,rel.tol=1e-10)
2.208542e-20 with absolute error < 4.4e-20
> integrate(function(x) dnorm(x)/(1+exp(-mu-x*sigma)),-.01,0,rel.tol=1e-10)
4.014838e-06 with absolute error < 2.1e-14

The first one of those is obviously wrong since it is suppose to be larger than the other one. However, what is presumably happening is that the integration algorithm looks at values "reasonably close" to the interval endpoint and gets, like,

> f2
function(x) dnorm(x)/(1+exp(-mu-x*sigma))
> f2(-.01)
[1] 5.392315e-78

and concludes that the integrand is essentially zero over the domain of integration.

One might consider splitting the domain like this:

> a <- 8/sigma
> I1 <- integrate(function(x) dnorm(x,mu/sigma)/(1+exp(-x*sigma)),-Inf,-a,rel.tol=1e-10)$value
> I2 <- integrate(function(x) dnorm(x,mu/sigma)/(1+exp(-x*sigma)),-a,a,rel.tol=1e-10)$value
> I3 <- integrate(function(x) dnorm(x,mu/sigma)/(1+exp(-x*sigma)),a,Inf,rel.tol=1e-10)$value
> I1+I2+I3
[1] 0.4999626

(I'm still not quite happy with the fact that I1 gets essentially set to zero in the above, but the essential point is that you need to be prepared to do a little analytical work to deal with tricky integrands.)

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: R numerical integration

casperyc
In reply to this post by casperyc
The quadinf command in library  pracma still fails when mu=-2.986731 with sigma=53415.18.
While Maple gives me an estimate of 0.5001701024.
########################################
Maple: (for those who are interested)
myf:=(mu,sigma)-> evalf(Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)), x=-infinity..infinity));
myf(-2.986731, 53415.18 );
                                        0.5001701024
########################################

These 'mu's and 'sigma's are now random starting points I generated for an optimization problem like I have mentioned.

I should really investigate the behavior of this function before I ask R doing the integration. As I have mentioned that I had already realized the integral is between 0 and 1. And I have had a look at the contour plots of different mu and sigma. I am going to 'restrict' mu and sigma to certain (small) values, and still get the integral to produce a value between 0 and 1.

All of your help is much appreciated.

casper
######################
PhD candidate in Statistics
Big R Fan
Big LEGO Fan
Big sTaTs Fan
######################
Reply | Threaded
Open this post in threaded view
|

Re: R numerical integration

Hans W Borchers
casperyc <casperyc <at> hotmail.co.uk> writes:
>

I don't know what is wrong with your Maple calculations, but I think
you should check them carefully, because:

(1) As Petr explained, the value of the integral will be < 0.5
(2) The approach of Peter still works and returns : 0.4999777
(3) And the same result comes out with Mathematica: 0.49997769...

Hans Werner

> The quadinf command in library  pracma still fails when mu=-2.986731 with
> sigma=53415.18.
> While Maple gives me an estimate of 0.5001701024.
> ########################################
> Maple: (for those who are interested)
> myf:=(mu,sigma)->
> evalf(Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)),
> x=-infinity..infinity));
> myf(-2.986731, 53415.18 );
>                                         0.5001701024
> ########################################
>
> These 'mu's and 'sigma's are now random starting points I generated for an
> optimization problem like I have mentioned.
>
> I should really investigate the behavior of this function before I ask R
> doing the integration. As I have mentioned that I had already realized the
> integral is between 0 and 1. And I have had a look at the contour plots of
> different mu and sigma. I am going to 'restrict' mu and sigma to certain
> (small) values, and still get the integral to produce a value between 0
> and 1.
>
> All of your help is much appreciated.
>
> casper
>

______________________________________________
[hidden email] 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.