[R] triple integral: adapt package question

[Ravi Varadhan]

Here is a simple-minded approach using indicator functions:

fxyz <- function(w) {
x <- w[1]
y <- w[2]
z <- w[3]
x*y*z*(y < x) * (z < x)

}


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

# exact answer is 1/24 = 0.041666

adapt(ndim=3, lower=lower, upper=upper, functn = fxyz, eps=1.e-03)

This approach can be very inefficient depending upon how big the 3-dim
hyper-rectangle is compared to the actual region over which the integrand is
non-zero.


Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 

----------------------------------------------------------------------------
--------

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Davood Tofighi
Sent: Friday, March 07, 2008 3:01 PM
To: r-help at r-project.org
Subject: [R] triple integral: adapt package question

Dear All,

I have a function f(x,y,z)=exp(x^3+y^4+x^2*y+x*z^2+y/z) over D, where is D={
(x,y,z)| 0 <z<Inf, 0<y<c1*z, 0<x<c2*/y}. x,y,z are all vectors and c1 and c2
are constants. I tried the "adapt" package and I get some error. This is the
error message:

"Error in function (z, y, x)  : argument "x" is missing, with no default"

I included my R code. Can anyone please let me know how I can calculate the
numerical integral of such a function where some of the boundaries are
functions of the other variables?

require(adapt)
fn<-function(z,y,x) {exp(x^3+y^4+x^2*y+x*z^2+y/z)}
x<-runif(200);y<-runif(200);z<-runif(200);
c1<-.5;c2<-5;M<-100;        #M to represent infinity
i1<-adapt(3,lo=c(.0001,0,0),up=c(M,c1*z,c2/y),functn=fn)$value
print(i1);
 Thanks,
-- 
Davood Tofighi

	[[alternative HTML version deleted]]

______________________________________________
R-help at 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.