Evaluate multivariate integral (of dependent variables) using monte carlo in R

Question

I would like to evaluate

$\int_0^1\int_0^1 \exp\left(-a\sqrt{(x-b)^2+(y-c)^2}\right)dxdy$

using a Monte Carlo approach ($a$ is a scalar and $b,c$ are in $(0,1)$).

Given that $a,b,c$ are parameters, I think the most sensible thing to do would be to set up a grid of $a,b,c$ and optimize over these values, but I still need to calculate the integral given a fixed $a,b,c$. Can I do something like the following:

x<-rnorm(100)
y<-x^2
func_to_int<-function(a,b,c,xval,yval){
  exp(-a*(sqrt((xval-b)^2+(yval-c)^2)))
}

means_vec<-c()
for(j in 1:100){
  set.seed(j)
  u_y<-runif(10000)
  u_x<-runif(10000)
  sim_vec<-c()
  for(i in 1:length(u_x)){
    sim_vec[i]<-func_to_int(a=1,b=.5,c=.5,xval=u_x,yval=u_y)
  }
  means_vec[j]<-mean(sim_vec)
  if(j%%10==0){
    print(j)
  }

}

mean(means_vec)

This doesn't really seem correct to me as $x$ and $y$ are not independent. The correct (analytical) solution given $a=b=c=1$ is 0.484999, but what I want to know if this is a reasonable approach given that I know $x$ and $y$ are dependent.

Please note that I only mention R as that is the language I primarily use, but I am more just looking for a general method (if you happen to know a way in R, then great)!

Thanks in advance for your time.

Edit: to

Answer 1

Fixing $a, b, c$, you should be evaluating points at u_x and u_y rather than x and y.

Formula for u_x and u_y are unif(10000).

Edit:

func_to_int<-function(a,b,c,xval,yval){
  exp(-a*(sqrt((xval-b)^2+(yval-c)^2)))
}

means_vec<-c()
for(j in 1:5){
  set.seed(j)
  u_y<-runif(10000)
  u_x<-runif(10000)
  sim_vec<-c()
  for(i in 1:length(u_x)){
    sim_vec[i]<-func_to_int(a=1,b=1,c=1,xval=u_x[i],yval=u_y[i])
  }
  means_vec[j]<-mean(sim_vec)
}
print(means_vec)