Monte-Carlo method for definite integral in R

user_10798020 0
Viewed 610

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.

montecarlo r
2 Answers

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

Here's an approach that I've found helpful in estimating integrals. Consider the following:

So if we have a distribution G over the support A, then we can estimate the integral of f(x) over A with samples from a distribution G. The function f need not be strictly positive.

In your case, let X be uniformly distributed, then

since g(x)=1 for all x in the support. So you can estimate the integral with

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

You could also let X be a beta random variable and estimate the integral with

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

Here's another example for estimating a function over the positive real line.

f = function(x)
    abs(sin(x)) / (x+1)^2
x = rgamma(N, 2, 1/4)
fx = f(x)
gx = dgamma(x, 2, 1/4)
mean(fx / gx)

I don't remember all the theory (there are likely to be some limitations on when you can use this approach), but generally if you can sample from any distribution on the same support as the integral, you can estimate the integral as above. And if your choice of G has a density that is "closer" to the function f, the smaller your Monte Carlo error would be.

user_105829030
Related