How to compute Monte Carlo Integration in R

Monte Carlo integration is a technique for estimating the value of a definite integral using random sampling and then applying a deterministic method to estimate the integral.

For example, suppose we want to estimate the value of the integral f(x)dx from a to b. We can do this by randomly generating points in the interval [a,b] and then estimating the integral as the average value of f(x) at those points.

In other words, we are approximating the integral by a Riemann sum. The advantage of Monte Carlo integration is that it can be used to estimate integrals that are difficult or impossible to compute using traditional methods.

Monte Carlo integration

Monte Carlo integration takes a number of random points, which are used to calculate the area under a curve or on a definite integral equation.

As the number of random points increases, the more accurate our integration becomes. This basically implies that more random variates = more accuracy. 

Example of Monte Carlo Integration

Approximate the integral KaTeX parse error: KaTeX doesn't work in quirks mode. using the Monte Carlo method.

The code for the above steps is shown below.

#Set the number of values
n <- 10000
#Generate n random numbers on the interval [0,1]
x  <- runif(n,min=0,max=1)
#Evaluate the function at the random values of s
fx <- x^2
#Estimate the integral
integral <- mean(fx)*1 # Multiply by the base length 1-0
integral

[1] 0.3336095

This output may vary due to the degree of randomness set in our R program. In order to stop this variation, we need to use the command set. seed( )

Applying the set.seed() function in our Monte Carlo integration.

set.seed(123)

#Set the number of values
n <- 10000
#Generate n random numbers on the interval [0,1]
x  <- runif(n,min=0,max=1)
#Evaluate the function at the random values
fx <- x^2
#Estimate the integral
integral <- mean(fx)*1 # Multiply by the base length 1-0
integral

After executing the above R code, our output should be the same.

[1] 0.3297405

Example 2

Compute the exact value of the integral KaTeX parse error: KaTeX doesn't work in quirks mode. and use Monte Carlo integration to obtain a value of the estimate.

set.seed(123)

#Set the number of values
n <- 10000
#Generate n random numbers on the interval [0,6]
x  <- runif(n,min=0,max=6)
#Evaluate the function at the random values
fx <- 4*x^2
#Estimate the integral
integral <- mean(fx)*6 # Multiply by the base length 6-0
integral

[1] 284.8957

The exact value is 2888.

Example 3

Compute the integral KaTeX parse error: KaTeX doesn't work in quirks mode. using monte carlo integration

set.seed(123)

#Set the number of values
n <- 10000
#Generate n random numbers on the interval 
t  <- runif(n,min=0,max=pi/3)
#Evaluate the function at the random values
ft <- sin(t)
#Estimate the integral
integral <- mean(ft)*pi/3 # Multiply by the base length Pi-0
integral

[1] 0.49808

Example 4: Estimation of Pi using Monte Carlo Integration

# Monte Carlo estimation of Pi
set.seed(123)  # for reproducibility
n <- 10000     # number of random points

x <- runif(n)
y <- runif(n)

inside_circle <- sum(x^2 + y^2 <= 1)
pi_estimate <- 4 * inside_circle / n

cat("Estimated Pi:", pi_estimate, "\n")

Estimated Pi: 3.18 

Example 4: Multidimensional Monte Carlo integration

Find the integral KaTeX parse error: KaTeX doesn't work in quirks mode.

To integrate a scalar function over a multidimensional rectangle, use the R function adaptIntegrate() from the library ‘cubature’

library(cubature)
f<- function(x){
  (2/3)*(x[1]+x[2]+x[3])
}
adaptIntegrate(f, lowerLimit = c(0,0,0), upperLimit = c(0.5,0.5,0.5))

> adaptIntegrate(f, lowerLimit = c(0,0,0), upperLimit = c(0.5,0.5,0.5))
$integral
[1] 0.0625

$error
[1] 1.387779e-17

$functionEvaluations
[1] 33

$returnCode
[1] 0

Example 5: Estimating a 2D region’s area

# Monte Carlo estimation of a 2D region's area
f <- function(x, y) x^2 + y^2 <= 1  # Function indicating the region

set.seed(123)
n <- 10000

x <- runif(n, -1, 1)
y <- runif(n, -1, 1)

inside_region <- sum(f(x, y))
area_estimate <- 4 * inside_region / n  # Multiplying by 4 to account for the entire square

cat("Estimated Area:", area_estimate, "\n")

Estimated Area: 3.1576 

Let’s visualize the above 2D region area:

# Monte Carlo estimation and plot of a 2D region's area
library(ggplot2)

# Function indicating the region
f <- function(x, y) x^2 + y^2 <= 1

set.seed(123)
n <- 1000

x <- runif(n, -1, 1)
y <- runif(n, -1, 1)

inside_region <- f(x, y)
outside_region <- !inside_region

# Create a data frame for the points
df <- data.frame(x = x, y = y, inside = inside_region)

# Create the plot
ggplot(df, aes(x, y, color = inside)) +
  geom_point() +
  scale_color_manual(values = c("red", "blue"), guide = FALSE) +
  theme_minimal() +
  labs(
    title = "Monte Carlo Estimation of 2D Region's Area",
    subtitle = paste("Estimated Area:", sum(inside_region) / n * 4),
    x = "X-Axis",
    y = "Y-Axis"
  )

Additional Resources

Why are Monte Carlo Simulations used?

How to plot a function curve in R?

How to plot ROC curve in Python?