Реализуйте метод симуляции Монте-Карло для оценки интеграла в R

Question

Я пытаюсь реализовать метод симуляции Монте-Карло для оценки интеграла в R. Однако я все еще получаю неправильный ответ. Мой код выглядит следующим образом:

f <- function(x){
  ((cos(x))/x)*exp(log(x)-3)^3
}
t <- integrate(f,0,1)

n <- 10000 #Assume we conduct 10000 simulations
int_gral <- Monte_Car(n)
int_gral

Answer 1

Вы не делаете Монте-Карло здесь. Монте-Карло - это метод моделирования, который помогает аппроксимировать интегралы с использованием сумм / средних на основе случайных величин.

Вы должны сделать что-то в этом духе (возможно, вам придется убедиться, что правильно сказать, что среднее значение вывода f может приблизительно соответствовать вашему интегралу:

f <- function(n){
  x <- runif(n)
  return(
    ((cos(x))/x)*exp(log(x)-3)^3
  )
}
int_gral <- mean(f(10000))

То, что делает ваш код, таково: взяв число n и вернув ((cos(n))/n)*exp(log(n)-3)^3, в этом нет случайности

Обновление

Теперь, чтобы получить более точные оценки, вам нужно повторить этот шаг K раз Вместо использования loop вы можете использовать функцию replicate:

K <- 100
dist <- data.frame(
  int = replicate(K, mean(f(10000)))
)

Вы получаете распределение оценок для вашего интеграла:

library(ggplot2)
ggplot(dist) + geom_histogram(aes(x = int, y = ..density..))

и вы можете использовать mean для получения числового значения:

mean(dist$int)
# [1] 2.95036e-05

Вы можете оценить точность своих оценок с помощью

sd(dist$int)
# [1] 2.296033e-07

Здесь он маленький, потому что N уже большой, что дает вам хорошую точность первого шага.

Answer 2

Мне удалось изменить коды следующим образом. Пожалуйста, подтвердите мне, что я поступаю правильно.

С уважением.

f <- function(x){
      ((cos(x))/x)*exp(log(x)-3)^3
    }
    set.seed(234)
    n<-10000
    for (i in 1:10000) {
      x<-runif(n)
      I<-sum(f(x))/n
    }
    I