a rapid method for two dimensional integration in R?

TheEcologist

TheEcologist

Global Moderator
#1
#1
Does anybody know a rapid method for two dimensional integration in R?

library(adapt)
?adapt

is painfully slow.

Example:

I need something that can (numerically) integrate this function (it's a 2d isotropic version of the weibull):

fx=function(x,y,k=c(2000,20))
{(1/(2*pi)) * (k[2]/(k[1]^k[2])) * sqrt(x^2+y^2)^(k[2]-2) *
exp(-(sqrt(x^2+y^2)/k[1])^k[2])}

I don’t need an analytic solution as only the example I give will be simple enough for that... and not the actual models I work with (which is not the above).

Anybody have a way of rapid 2d numerical integration?


Thanks,
 
TheEcologist

TheEcologist

Global Moderator
#3
#3
Martingale said:
Do you have access to MATLAB or Maple or Mathmatica?

Edit: I don't know of a faster way in R
Click to expand...
I do, however the integration needs to take place as part of a ML fit procedure.
I have written the entire fit procedure in R and it would be a major pain to translate it into another code!

Nevertheless that is an option, however why do you think it would be faster in e.g. Mathlab?
 
Dragan

Dragan

Super Moderator
#5
#5
TheEcologist said:
Does anybody know a rapid method for two dimensional integration in R?

library(adapt)
?adapt

is painfully slow.

Example:

I need something that can (numerically) integrate this function (it's a 2d isotropic version of the weibull):

fx=function(x,y,k=c(2000,20))
{(1/(2*pi)) * (k[2]/(k[1]^k[2])) * sqrt(x^2+y^2)^(k[2]-2) *
exp(-(sqrt(x^2+y^2)/k[1])^k[2])}


Anybody have a way of rapid 2d numerical integration?


Thanks,
Click to expand...

I'm using Mathematica 6.0 and I can't get an answer using various integration techniques e.g. AdapativeMonteCarlo.

Can you more explicitly write fx (with the upper and lower limits that your using in R) so that I can double check my Mathematica notation.

I'm guessing - but I think I might be making a mistake by not expressing fx correctly in NIntegrate.
 
T

Tart

New Member
#6
#6
Hey TheEcologist

How accurate results should be? If you don't need super accurate results, then you can simply use summation to approximate two dimensional integration. I had similar problem once (only in MATLAB), and default integration function was super slow. I gave up on default integration and simply used summation over grid. (recall we define integrals as limits of sums over grids). If you make grid small enough and have more or less efficient code :) your approximation will work much faster than default functions and results will be very reasonable.

Open source alternative to mathematica and maple is SAGE -http://www.sagemath.org/ I recently installed it, don't ask me any questions about. Looks cool and powerfull.

One other solution is to find C, C++, or Fortran library for numerical integration, I never did it before, but within R you can send computationally intense procedures to such libraries. I know that many R packages do exactly that. So it is theoretically possible, I just don't know how.
 
Last edited:
F

fasteddie

New Member
#7
#7
If you do take the advice of Tart (which makes sense IMO) make sure to multiply your answer by 'dx' and 'dy', the grid units ... easily forgotten :)
 
TheEcologist

TheEcologist

Global Moderator
#8
#8
Tart said:
Hey TheEcologist

How accurate results should be? If you don't need super accurate results, then you can simply use summation to approximate two dimensional integration. I had similar problem once (only in MATLAB), and default integration function was super slow. I gave up on default integration and simply used summation over grid. (recall we define integrals as limits of sums over grids). If you make grid small enough and have more or less efficient code :) your approximation will work much faster than default functions and results will be very reasonable.

Open source alternative to mathematica and maple is SAGE -http://www.sagemath.org/ I recently installed it, don't ask me any questions about. Looks cool and powerfull.

One other solution is to find C, C++, or Fortran library for numerical integration, I never did it before, but within R you can send computationally intense procedures to such libraries. I know that many R packages do exactly that. So it is theoretically possible, I just don't know how.
Click to expand...
Yes, that is exactly what I have been using.

x=seq(xmin,xmax,dx)
y=seq(ymin,ymax,dy)
outer(x,y,fx [dx.dy])

Only problems is that then I do have limits (xmin,xmax,ymin,ymax) and these limits are again limited to the limits of my limited laptop memory. Right now I'm starting to write my own (more flexible) alogoritm, using the grid method.
It is used in a model fit procedure it needs to be flexible. E.g. let the grid size increase with increasing distance from (0,0).

I'll update this with code if I gets something decent working.
Anyway thanks for the comments :D
 
You must log in or register to reply here.