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 School of Mathematics, Statistics and Actuarial Science, University of Kent ################################################### -- View this message in context: http://r.789695.n4.nabble.com/R-numerical-integration-tp4500095p4500095.html Sent from the R help mailing list archive at Nabble.com.
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
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.
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
School of Mathematics, Statistics and Actuarial Science, University of Kent
###################################################
--
View this message in context:
http://r.789695.n4.nabble.com/R-numerical-integration-tp4500095p4503766.html
Sent from the R help mailing list archive at Nabble.com.