Re: [R] How to double integrate a function in R
Tiago V. Pereira tiago.pereira at mbe.bio.br writes:
I am trying to double integrate the following expression:
# expression
(1/(2*pi))*exp(-y2/2)*sqrt((y1/(y2-y1)))
for y2y10.
I am trying the following approach
# first attempt
library(cubature)
fun - function(x) { (1/(2*pi))*exp(-x[2]/2)*sqrt((x[1]/(x[2]-x[1])))}
adaptIntegrate(fun, lower = c(0,0), upper =c(5, 6), tol=1e-8)
However, I don't know how to constrain the integration so that y2y10.
Any ideas?
Tiago
You could use integral2() in package 'pracma'. It implements the
TwoD algorithm and has the following properties:
(1) The boundaries of the second variable y can be functions of the first
variable x;
(2) it can handle singularities on the boundaries (to a certain extent).
library(pracma)
fun - function(y1, y2) (1/(2*pi))*exp(-y2/2)*sqrt((y1/(y2-y1)))
integral2(fun, 0, 5, function(x) x, 6, singular=TRUE)
$Q
[1] 0.7706771
$error
[1] 7.890093e-11
The relative error is a bit optimistic, the absolute error here is 0.5e-6.
The computation time is 0.025 seconds.
Hans Werner
__
R-help@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.
Re: [R] How to double integrate a function in R
On Jul 26, 2013, at 8:44 AM, Tiago V. Pereira wrote:
I am trying to double integrate the following expression:
# expression
(1/(2*pi))*exp(-y2/2)*sqrt((y1/(y2-y1)))
for y2y10.
I am trying the following approach
# first attempt
library(cubature)
fun - function(x) { (1/(2*pi))*exp(-x[2]/2)*sqrt((x[1]/(x[2]-x[1])))}
adaptIntegrate(fun, lower = c(0,0), upper =c(5, 6), tol=1e-8)
However, I don't know how to constrain the integration so that y2y10.
Generally incorporating boundaries is accomplished by multiplying the integrand
with logical vectors that encapsulate what are effectively two conditions:
Perhaps:
fun - function(x) { (x[1]x[2])*(x[1]0)* (1/(2*pi))*exp(-x[2]/2)*
sqrt((x[1]/(x[2]-x[1])))}
That was taking quite a long time and I interrupted it. There were quite a few
warnings of the sort
1: In sqrt((x[1]/(x[2] - x[1]))) : NaNs produced
2: In sqrt((x[1]/(x[2] - x[1]))) : NaNs produced
Thinking the NaNs might sabotage the integration process, I added a conditional
to the section of that expression that was generating the NaNs. I don't really
know whether NaN's are excluded from the summation process in adaptIntegrate:
fun - function(x) { (x[1]x[2])*(x[1]0)* (1/(2*pi))*exp(-x[2]/2)*
if(x[1]x[2]){ 0 }else{ sqrt((x[1]/(x[2]-x[1]))
)} }
adaptIntegrate(fun, lower = c(0,0), upper =c(5, 6) )
I still didn't have the patience to wait for an answer, but I did plot the
function:
fun2 - function(x,y) { (xy)*(x0)* (1/(2*pi))*exp(-y/2)* sqrt((x/(y-x)))}
persp(outer(0:5, 0:6, fun2) )
So at least the function is finite over most of its domain.
--
David Winsemius
Alameda, CA, USA
__
R-help@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.
Re: [R] How to double integrate a function in R
On Jul 26, 2013, at 9:33 AM, David Winsemius wrote:
fun2 - function(x,y) { (xy)*(x0)* (1/(2*pi))*exp(-y/2)* sqrt((x/(y-x)))}
persp(outer(0:5, 0:6, fun2) )
There does seem to be some potential pathology at the edges of the range,
Restricting it to x = 0.03 removes most of that concern.
fun2 - function(x,y) { (xy)*(x0)* (1/(2*pi))*exp(-y/2)* sqrt((x/(y-x)))}
persp(outer(seq(0.01,5,by=.01), seq(.02,6,by=.01), fun2) ,ticktype=detailed)
fun - function(x) { (x[1]x[2])*(x[1]0)*
(1/(2*pi))*exp(-x[2]/2)*if(x[1]x[2]){0}else{ sqrt((x[1]/(x[2]-x[1])) )}}
adaptIntegrate(fun, lower = c(0.03,0.03), upper =c(5, 6), tol=1e-2 )
$integral
[1] 0.7605703
$error
[1] 0.00760384
$functionEvaluations
[1] 190859
$returnCode
[1] 0
I tried decreasing the tolerance to 1e-3 but the wait exceeds the patience I
have allocated to the problem.
--
David Winsemius
Alameda, CA, USA
__
R-help@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.
