Search code examples
rsimulationmontecarlointegral

Monte carlo simulation - what's wrong?

I want to calculate area under the curve using monte carlo simulation of function

enter image description here. I want to calculate it on interval [-2, 2]

My work so far

# Define function f
f <- function(x) x^2 + 1

# I want to close my function in rectangle (-2, 2) - x axis and (1, 5) y -axis
n <- 10^6
# Randomize from x axis of rectangle
x_n <- runif(n, min = -2, max = 2)
# Randomize from y axis of rectangle
y_n <- runif(n, min = 1, max = 5)
# Calculate function values of randomized points
values <- f(x_n)

# Formula for are under the curve for monte carlo simulation is 
# Area of rectangle * (Points below curve) / (Total number of points)

So my result is:

> sum(y_n < values) / n * (4 * 4)
[1] 5.329888

which is bad result (correct result is 9.33333). What I'm doing wrong? For sure algorithm should have been much closer to 9.3333 after milion sampling

Solution

  • Here's a plot to show what you were working with. I'm hoping it will help you understand what I wrote in my comment better:

    enter image description here

    You seem to be ignoring the rectangle below y=1. It's area (=4) is the missing quantity. So the code is correct for calculating the non-offset expression x^2. Change to y_n <- runif(n, min = 0, max = 5) and re-run the calculations

    The comment was half-the answer, i.e. that you hadn't simulated the point that were between 0 and 1 for y_n's. Those need to be in the Monte Carlo model of integration of an area. The other modification is to add in the correct total area of [-2 < x <2]x[0<y<5] = 4*5 to the calculation of the total "area" under consideration

    f <- function(x) x^2 + 1

# I want to close my function in rectangle (-2, 2) - x axis and (1, 5) y -axis
n <- 10^6
# Randomize from x axis of rectangle
x_n <- runif(n, min = -2, max = 2)
# Randomize from y axis of rectangle
y_n <- runif(n, min = 0, max = 5)
# Calculate function values of randomized points
values <- f(x_n)

# Formula for are under the curve for monte carlo simulation is 
# Area of rectangle * (Points below curve) / (Total number of points)
 sum(y_n < values) / n * (5 * 4)
#[1] 9.3429  inaccurate to 1 or 2 parts in 933

    Plot of 100 points to show second situation:

    enter image description here

    The other mod you might consider is using set.seed to make your calculations reproducible.