How to approximate the solution to a double integral in R using hit and miss and Halton?

Question

I have the following question:

what I do is:

set.seed(1)
N = 10000
f = function(x,y) x^y * y^x

then I don't know how to proceed. Can somebody please explain to me how to do this? thanks!

Answer 1

The Monte-Carlo approximation of the integral is obtained as follows:

set.seed(1)
f <- function(x,y) x^y * y^x

N <- 10000
mean(f(runif(N), runif(N)))
# 0.4293375

This is probably what is called the "uniform method" in the question.

To compare with a numerical approximation (better):

library(cubature)
f <- function(xy){
  x <- xy[1]; y <- xy[2]
  x^y * y^x
} 
integral <- pcubature(f, lowerLimit = c(0,0), upperLimit = c(1,1))
integral$integral
# 0.4280186

Wolfram|Alpha gives 0.42802.

Now, the approximation obtained by sampling a Halton sequence:

library(randtoolbox)
setSeed(1)
xy <- halton(n = N, dim = 2)
f <- function(x,y) x^y * y^x
mean(f(xy[,1], xy[,2]))
# 0.4277797

Note that I'm not sure this is the desired answer, because I don't know yet what is the "Hit-or-miss" method.

EDIT

I took a look at the Hit-or-miss method. It consists in sampling points in the 3D space, and taking the proportion of these points that fall below the surface defined by f. Here the function f takes its values in [0,1], so the approximation is obtained as follows:

library(randtoolbox)
setSeed(1)
xyz <- halton(n = N, dim = 3)
f <- function(x,y) x^y * y^x
mean(f(xyz[,1],xyz[,2]) > xyz[,3])
# 0.4287