Dear All
Thanks to Steven, Carlisle and Alec. This is a summary of responses. I hope
that someone will find new information about double integration and add it
to this summary.
Question:
How to do a double integration of f(x,y) = 20/(y*x^2) over {x,2,3} and
{y,1,5}
Solutions:
1) Steven introduces MC method;
#2-dimensional Monte Carlo integration
M<-100000 # points
rangex<-c(2,3)
rangey<-c(1,5)
#f<-function(x,y){20/(x*y^2)}
f<-function(x,y){x*y^2}
fint<-mean(f(runif(M,rangex[1],rangex[2]),runif(M,rangey[1],rangey[2])))*
diff(rangex)*diff(rangey)
#2-d Riemann sums
N=100
xx<-rep(seq(from=2,to=3,length=N),N)
yy<-rep(seq(from=1,to=5,length=N),N)
int2<-0
for(yz in yy){
int2<-int2+sum(f(xx,rep(yz,10) ) )
#print(sum(f(xx,rep(yz,10) ) ) )
}
int2.x<-int2*diff(rangex)*diff(rangey)/N^2
2)Carlisle shows the approximation.
First make vectors appropriate for dx=dy=0.1
> test.x <- (21:30)/10 - 0.05
> test.y <- (11:50)/10 - 0.05
The integral is a sum:
> sum(20/outer(test.y,test.x^2))*(0.1)^2
[1] 5.361106
Compare with the analytical solution:
> 20*log(5)/6
[1] 5.364793
3)Alec forwards a past discussion by Andrzej Galecki.
In response to the recent message on double integration: the subject was
discussed in March 99 (see S-Archives).
However the functions given there should be used cautiously. The basic
code proposed was of the form
double.integral.demo.fn <- function() {
fun <- function(y)
unlist(lapply(y, function(y)
{
integrate(f = function(x, y)
exp(-0.5 * (x^2 + y^2))/(2 * pi), lower = IL, upper =
IU, y = y)$integral
}))
integrate(f = fun, lower = OL, upper = OU)$integral
}
where IL,IU,OL and OU represent limits of integration (inner lower, inner
upper etc..) and could in general be given as function arguments.
the function given within the code is the joint density of two independent
standard normal rvs
consider those square regions of integration where L=IL=OL, U=IU=OU
on R2 (L=-Inf,U=Inf)
> double.integral.demo.fn()
[1] 0.6853893
on R2+ (L=0,U=Inf)
> double.integral.demo.fn()
[1] 0.1713473
on R2- (L=-Inf,U=0)
> double.integral.demo.fn()
[1] 0.1713473
All clearly incorrect. No warnings or errors were printed.
I've tested the function over finite squared regions of integration -
similar errors occur in some circumstances, sometimes with warnings and
sometimes without. In some cases the algorithm seems to get stuck since it
continues to evaluate the inner integral (in x) at the same y value
(y=0) and hence increasing the subdivisions has no affect on the
output. The only way around these problems appears to be to lessen the
tolerance levels.
e.g try, with default tolerance levels (abs.tol = 10^-4 say)
(L=-c,U=c): works for c<=2.4, not for c>=2.5
(L=-c,U=5): works for c<=0.6, not for c>=0.7
Other functions put on the list based on this code (e.g the
integrate.double function from Karen Johnson [original code was from Bill
Venables]) yield the same conclusions.
|