Why is my Monte Carlo Integration wrong by a factor of 2?

I am trying to integrate the following function using a Monte Carlo Integration. The interval I want to integrate is x <- seq(0, 1, by = 0.01) and y <- seq(0, 1, by = 0.01).

my.f <- function(x, y){
  result = x^2 + sin(x) + exp(cos(y))
  return(result)
}

I calculated the integral using the cubature package.

library(cubature)
library(plotly)

# Rewriting the function, so it can be integrated
cub.function <- function(x){
  result = x[1]^2 + sin(x[1]) + exp(cos(x[2]))
  return(result)
}
cub.integral <- adaptIntegrate(f = cub.function, lowerLimit = c(0,0), upperLimit = c(1,1))

The result is 3.134606. But when I use my Monte Carlo Integration Code, see below, my result is about 1.396652. My code is wrong by more than a factor of 2!

What I did:

Since I need a volume to conduct a Monte Carlo Integration, I calculated the function values on the mentioned interval. This will give me an estimation of the maximum and minimum of the function.

# My data range
x <- seq(0, 1, by = 0.01)
y <- seq(0, 1, by = 0.01)

# The matrix, where I save the results
my.f.values <- matrix(0, nrow = length(x), ncol = length(y))

# Calculation of the function values
for(i in 1:length(x)){
  for(j in 1:length(y)){
    my.f.values[i,j] <- my.f(x = x[i], y = y[j])
  }
}

# The maximum and minimum of the function values
max(my.f.values)
min(my.f.values)

# Plotting the surface, but this is not necessary
plot_ly(y = x, x = y, z = my.f.values) %>% add_surface()

So, the volume that we need is simply the maximum of the function values, since 1 * 1 * 4.559753 is simply 4.559753.

# Now, the Monte Carlo Integration
# I found the code online and modified it a bit. 

monte = function(x){
  tests = rep(0,x)
  hits = 0
  for(i in 1:x){
    y = c(runif(2, min = 0, max = 1),     # y[1] is y; y[2] is y
    runif(1, min = 0, max = max(my.f.values))) # y[3] is z
    if(y[3] < y[1]**2+sin(y[1])*exp(cos(y[2]))){
      hits = hits + 1
    }
    prop = hits / i
    est = prop * max(my.f.values)
    tests[i] = est
  }
  return(tests)
}
size = 10000
res = monte(size)
plot(res, type = "l")
lines(x = 1:size, y = rep(cub.integral$integral, size), col = "red")

So, the result is completely wrong. But if I change the function a bit, suddenly is works.

monte = function(x){
  tests = rep(0,x)
  hits = 0
  for(i in 1:x){
    x = runif(1)
    y = runif(1)
    z = runif(1, min = 0, max = max(my.f.values))
    if(z < my.f(x = x, y = y)){
      hits = hits + 1
    }
    prop = hits / i
    est = prop * max(my.f.values)
    tests[i] = est
  }
  return(tests)
}
size = 10000
res = monte(size)
plot(res, type = "l")
lines(x = 1:size, y = rep(cub.integral$integral, size), col = "red")

Can somebody explain why the result suddenly changes? To me, both functions seem to do the exact same thing.

Solution

In your (first) code for monte, this line is in error:

y[3] < y[1]**2+sin(y[1])*exp(cos(y[2]))

Given your definition of my.f, it should surely be

y[3] < y[1]**2 + sin(y[1]) + exp(cos(y[2]))

Or..., given that you shouldn't be repeating yourself unnecessarily:

y[3] < my.f(y[1], y[2])