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.