MHB Computing Triple Integral in 'R'

[WMDhamnekar]

I would like to compute the triple integral of a function of three variables $f(x,y,z)$in R. I am using the package Cubature, Base, SimplicialCubature and the function adaptIntegrate(), Integrate and adaptIntegrateSimplex(). The integrand is equal to 1 only in certain domain(x<y<z, 0 otherwise).
Followings are the different ways to answer this question but i didn't understand how to compute the answer 0.166666 manually. If any member knows the answer, may reply. Using Cubature package,

library(cubature)
lower <- rep(0,3)
upper <- rep(1,3)

fxyz <- function(w) {
x <- w[1]
y <- w[2]
z <- w[3]
as.numeric(x <= y)*as.numeric(y <= z)
}
adaptIntegrate(f=fxyz,lowerLimit=lower,upperLimit= upper,doChecking=TRUE, maxEval=2000000,absError=10e-5,tol=1e-5)
\$integral
[1]0.1664146

\$error
[1]0.0001851699

\$functionEvaluations
[1]2000031

\$returnCode
[1]0

Using base package in 'R'

f.xyz <- function(x, y, z) ifelse(x < y & y < z, 1, 0)
f.yz <- Vectorize(function(y, z) integrate(f.xyz, 0, 1, y=y, z=z)\$value,
vectorize.args="y")

f.z <- Vectorize(function(z) integrate(f.yz, 0, 1, z=z)\$ value,
vectorize.args="z")
integrate(f.z, 0, 1)
>0.1666632 with absolute error < 9.7e-05

Using SimpicialCubature package in 'R'

library(SimplicialCubature)
f <- function(x) 1
S <- CanonicalSimplex(3)
> adaptIntegrateSimplex(function(x) 1, S)
\$integral
[1] 0.1666667

\$estAbsError
[1] 1.666667e-13

\$functionEvaluations
[1] 55

$returnCode
[1] 0

\$message
[1] "OK"

Note that integrating the constant function f(x)=1 over the simplex simply gives the volume of the simplex, which is 1/6. The integration is useless for this example

SimplexVolume(S)
[1] 0.1666667

Can any member explain me what is this actual question and how we can solve this question manually?

 

[HallsofIvy]

I don't know "Cubature" but I will reply anyway!

First, your domain, x< y< z, is not well defined. The simple "x< y< z" is infinite and the integral of "1" on that domain is not finite. Is there not some other condition bounding the domain? In your program you have "lower<- rep(0,3)" and "upper<- rep(1,3)". Could that possibly mean that x lies between 0,3 and 1,3 (I would have said 0.3< z< 1.3)?

 

[WMDhamnekar]

I don't know "Cubature" but I will reply anyway!

First, your domain, x< y< z, is not well defined. The simple "x< y< z" is infinite and the integral of "1" on that domain is not finite. Is there not some other condition bounding the domain? In your program you have "lower<- rep(0,3)" and "upper<- rep(1,3)". Could that possibly mean that x lies between 0,3 and 1,3 (I would have said 0.3< z< 1.3)?

Hello,
lower<-rep(0,3) means in adaptIntegrate function lower limit 0,repeats 3 times and upper<-(1,3) means upper limit 1, repeats 3 times.
e.g.$\displaystyle\int_0^1\displaystyle\int_0^1\displaystyle\int_0^1$

 

[HallsofIvy]

Hello,
lower<-rep(0,3) means in adaptIntegrate function lower limit 0,repeats 3 times and upper<-(1,3) means upper limit 1, repeats 3 times.
e.g.$\displaystyle\int_0^1\displaystyle\int_0^1\displaystyle\int_0^1$

That's strange. What would happen if you entered different number of repetitions for the the lower limits and upper limits? Also \int_0^1\int_0^1\int_0^1 doesn't match your "x< y< z".

 

[WMDhamnekar]

That's strange. What would happen if you entered different number of repetitions for the the lower limits and upper limits? Also \int_0^1\int_0^1\int_0^1 doesn't match your "x< y< z".

Hello,
Because of the following programming command while using 'Cubature' package in 'R', adaptIntegrate function gives answer 0.1664146.

fxyz <- function(w) {
x <- w[1]
y <- w[2]
z <- w[3]
as.numeric(x <= y)*as.numeric(y <= z)
}

This programming command seems to be difficult to understand.