Метод Монте-Карло для определенного интеграла в R - PullRequest
Новая
       13

Метод Монте-Карло для определенного интеграла в R

0 голосов
/ 16 декабря 2018

Я начал Msc, где мы изучаем пакет R, но у меня проблемы с упражнением.Упражнение таково:

"Предположим, что мы хотим оценить ∫x2 (с определенными значениями от 1 до 0), используя базовый метод Монте-Карло. По сути, мы бросаем дротики на кривую и подсчитываем количество выпадающих дротиковниже кривой Алгоритм метода состоит из:

1) Инициализация: хиты = 0

2) для (i в 1: N): генерировать два случайных числа, U1, U2между 0 и 1. Если U2 <(U1) ^ 2, то хиты = хиты + 1 конец для. </p>

3) Оценка площади = хиты / N

Предоставить R-код с использованием цикловдля реализации этого алгоритма Монте-Карло. Сколько времени занимает ваша функция? Предоставьте более эффективный код, избегая предыдущих циклов. Проиллюстрируйте выигрыш в эффективности, который можно получить путем векторизации вашего кода. "

У меня есть эти коды,но я думаю, что делаю это неправильно. 

montecarlo <- function(N){
  hits=0
  for (i in 1:N){
    U1 <- runif(1, 0, 1)
    U2 <- runif(1, 0, 1)
    if (U2 < (U1)^2){
      hits = hits+1}}
  return(hits/N)
}

montecarlo2 <- function(N){
  hits=0
  U1 <- runif (1:N, 0, 1)
  U2 <- runif (1:N, 0, 1)
  hits= hits+1 [U2<(U1)^2]
  return(hits/N)
}

Для первого метода с помощью циклов я получил (например): 

> montecarlo(23)
[1] 0.3478261
> montecarlo(852)
[1] 0.3591549
> montecarlo(8563255)
[1] 0.3332472

Можете ли вы мне помочь?Большое спасибо: S.

Ответы [ 2 ]

0 голосов
/ 17 декабря 2018

Вот подход, который я нашел полезным при оценке интегралов.Учтите следующее:

Таким образом, если у нас есть распределение G по носителю A, то мы можем оценить интеграл от f (x) по A с выборками из распределения G. Функция f не обязательно должна быть строго положительной.

В вашем случае, пусть X будет равномерно распределен, тогда

, поскольку g (x) = 1 для всех x в поддержке.Таким образом, вы можете оценить интеграл с помощью 

N = 100000
mean(runif(N)^2)

Вы также можете позволить X быть бета-случайной переменной и оценить интеграл с помощью 

x = rbeta(N, 2, 1)
fx = x^2
gx = dbeta(x, 2, 1)
mean(fx / gx)

Вот еще один пример оценки функции по положительномуРеальная линиякак интеграл, вы можете оценить интеграл, как указано выше.И если ваш выбор G имеет плотность, которая «ближе» к функции f, тем меньше будет ошибка Монте-Карло.

0 голосов
/ 16 декабря 2018

Один из способов: 

montecarlo_for <- function(N) {
  hits <- 0
  for (i in 1:N) {
    U1 <- runif(1)
    U2 <- runif(1)
    if (U2 < (U1) ^ 2) hits <- hits + 1
  }
  return(hits / N)
}

Векторизация 

montecarlo_vec <- function(N) {
  sum(runif(N) < runif(N)^2) / N
}

Сравните скорость, например, используя пакет microbenchmark: 

microbenchmark::microbenchmark(times = 50,
  montecarlo_for(1e5),
  montecarlo_vec(1e5)
)

Сравнение скоростина моей машине видно, что векторизованный метод примерно в 100 раз быстрее (среднее и среднее время показано ниже): 

Unit: milliseconds
expr                  mean       median
montecarlo_for(1e+05) 509.927001 497.238904
montecarlo_vec(1e+05) 5.214527   4.922007

Просто для удовольствия, если вы хотите посмотреть, насколько быстро алгоритм сходится к результату(1/3) с растущим размером выборки: 

plot(sapply(1:1000, montecarlo_vec), type = "line")
...