为什么我的 Monte Carlo 集成错误了 2 倍？

Question

我正在尝试使用 Monte Carlo 集成来集成以下功能。我要积分的区间是x <- seq(0, 1, by = 0.01)和y <- seq(0, 1, by = 0.01).

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

我使用 cubature 包计算了积分。

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))

结果是3.134606。但是当我使用我的 Monte Carlo 集成代码时，见下文，我的结果大约是 1.396652。我的代码错误超过 2 倍！

我做了什么：

因为我需要一个体积来进行Monte Carlo积分，所以我计算了上述区间的函数值。这将给我一个函数的最大值和最小值的估计值。

# 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()

所以，我们需要的体积就是函数值的最大值，因为 1 * 1 * 4.559753 就是 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")

所以，结果是完全错误的。但是如果我稍微改变一下功能，突然就可以了。

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")

谁能解释一下为什么结果突然变了？对我来说，这两个函数似乎做的是完全相同的事情。

Answer 1

在您 monte 的（第一个）代码中，这一行是错误的：

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

根据您对 my.f 的定义，它肯定是

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

或者...，鉴于您不应该不必要地重复自己：

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

为什么我的 Monte Carlo 集成错误了 2 倍？

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

r

montecarlo