1

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

Lucian
  • 351
  • 2
  • 10
  • 1
    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. – IRTFM Jun 05 '21 at 18:53

2 Answers2

3

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.

IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • I understand that result is much closer to real one, but why should I widen it? What's the intuition behind it? I thought that it's good to close this curve in rectangle [1, 5], becuase the minimal and maximal value of x^2 + 1 on [-2, 2] is 1 and 5. – Lucian Jun 05 '21 at 19:12
  • 1
    @Lucian When we talk about the area under a curve we generally mean the area between the curve and the x axis. What you have calculated is the area between the curves y=x^2+1 and y=1 on the interval x=[-2, 2]. – AkselA Jun 05 '21 at 19:18
  • @AkselA: see if the added plots are helpful in understanding the problem. – IRTFM Jun 05 '21 at 21:16
  • @AkselA. My apologies. I addressed the comment incorrectly. Should have been to Lucian. It was clear to me that you understood and clearly described the issues. – IRTFM Jun 06 '21 at 21:20
0

We can try the Monte Carlo simulation like this

> n <- 1e6

> x <- runif(n, -2, 2)

> y <- runif(n, 0, 5)

> mean(x^2 + 1 - y >= 0) * 4 * 5
[1] 9.33014

where the area can be calculated as the mean number of points laying under the curve x^2 + 1 -y >=0

ThomasIsCoding
  • 96,636
  • 9
  • 24
  • 81