Monte-Carlo method for definite integral in R

Question

I started a Msc where we are learning R package, but I am having problems with an exercise. The exercise is this:

"Suppose we want to estimate ∫x2 (definite between 1 and 0) using a basic Monte-Carlo method. Essentially, we throw darts at the curve and count the number of darts that fall below the curve. The algorithm of the method consist on:

1) Initialise: hits=0

2) for (i in 1:N): Generate two random numbers, U1,U2 between 0 and 1. If U2<(U1)^2, then hits=hits+1 end for.

3) Area estimate = hits/N

Provide a R-code using loops to implementing this Monte-Carlo algorithm. How much time takes your function? Provide a more efficient code avoiding the previous loops. Illustrate the efficiency gains that can be made by vectorising your code."

I have these codes, but I think I am doing it wrong.

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

For the first method, with loops, I have obtained (for example):

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

Can you help me? Thank you so much: S.

Answer 1

One of the ways:

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

Vectorized

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

Compare speed for example using the microbenchmark package:

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

The speed comparison on my machine shows that the vectorized method is roughly a 100 times faster (mean and median times shown below):

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

Just for the fun, if you want to look how quickly the algorithm converges to the result (1/3) with a growing sample size:

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