Почему моя интеграция Монте-Карло неверна в 2 раза? - PullRequest
Новая
       7

Почему моя интеграция Монте-Карло неверна в 2 раза?

0 голосов
/ 24 мая 2019

Я пытаюсь интегрировать следующую функцию, используя интеграцию Монте-Карло. Интервал, который я хочу интегрировать, составляет 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. Но когда я использую код интеграции Монте-Карло, см. Ниже, мой результат составляет около 1,396652. Мой код неверен более чем в 2 раза!

Что я сделал:

Так как мне нужен объем для проведения интеграции Монте-Карло, я рассчитал значения функций в указанном интервале. Это даст мне оценку максимума и минимума функции. 

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

Может кто-нибудь объяснить, почему результат внезапно меняется? Мне кажется, что обе функции делают одно и то же.

1 Ответ

4 голосов
/ 24 мая 2019

В вашем (первом) коде для 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])
...