How to estimate \(\pi\) using Monte Carlo?

Introduce a square \(\mathcal{S\subseteq}\mathbb{R}^{2}\), the sides being of length 2, with inscribed disk \(\mathcal{D}\) of radius 1. We are interested in computing the area \(I\) of \(\mathcal{D}\) using Monte Carlo. We have \[\begin{aligned} I & =\pi=\int\int_{\mathcal{D}}dx_{1}dx_{2}\\ & =\int\int_{\mathcal{S}}\mathbb{I}_{\mathcal{D}}\left( x_{1},x_{2}\right) dx_{1}dx_{2}\text{ as }\mathcal{D\subset S}\\ & =4\int\int_{\mathbb{R}^{2}}\mathbb{I}_{\mathcal{D}}\left( x_{1}% ,x_{2}\right) \mu\left( x_{1},x_{2}\right) dx_{1}dx_{2}% \end{aligned} \] where \(\mathcal{S}:=\left[ -1,1\right] \times\left[ -1,1\right]\) and \[ \mu\left( x_{1},x_{2}\right) =\frac{1}{4}\mathbb{I}_{\mathcal{S}}\left( x_{1},x_{2}\right) \] is the uniform density on the square \(\mathcal{S}\). In this case, we have \[ \widehat{I}_{n}=4\frac{n_{\mathcal{D}}}{n}% \] where \(n_{\mathcal{D}}\) is the number of samples which fell within the disk.

# Asign basic parameters of the experiment.
numPoints <- 1000  # number of samples
r <- 1  # radius of disk
C <- c(1, 1)  # center of disk and square
l <- 2  # size of square

# Draw 2-by-numPoints matrix of independent U[0, l] random variables.
U <- matrix(data = l * runif(n = 2 * numPoints), nrow = numPoints, ncol = 2)

# Compute index mask of all points which are inside the circle.
inside <- (U[, 1] - C[1])^2 + (U[, 2] - C[2])^2 < r^2

# Count the number points in U that fell inside the disk.
M <- sum(inside)
cat("Monte Carlo estimate of pi: ", (M/numPoints) * (l^2)/(r^2))
## Monte Carlo estimate of pi:  3.08
# A simple graph of the points: green for the points inside the circle and
# black for the others.
plot(U[inside, 1], U[inside, 2], type = "p", xlab = "X", ylab = "Y", xlim = c(0, 
    l), ylim = c(0, l), col = "green")
points(U[!inside, 1], U[!inside, 2], col = "black")

# Graph using ggplot2.
library(ggplot2)
theme_set(theme_bw())  # customize the graph style

# Plot the points.
graph_sample <- qplot(x = U[, 1], y = U[, 2], colour = inside, geom = "point")  # plot points
graph_sample <- graph_sample + xlab("X") + ylab("Y")  # change axis labels

# Draw the circle.
graph_sample <- graph_sample + annotate("path", x = C[1] + r * cos(seq(0, 2 * 
    pi, length.out = 100)), y = C[2] + r * sin(seq(0, 2 * pi, length.out = 100)), 
    colour = "blue", size = 2)

# Draw the rectangle.
graph_sample <- graph_sample + geom_rect(aes(xmin = 0, ymin = 0, xmax = l, ymax = l), 
    colour = "red", size = 2, fill = NA)

# Show the graph.
print(graph_sample)

# Repeat the above experiment 100 times.
numRepeats <- 100

# Store the results per experiment and per number of points used.
onlineErrors <- matrix(nrow = numPoints, ncol = numRepeats)
# For each experiment,
for (i in 1:numRepeats) {
    # draw a uniform sample,
    U <- matrix(data = l * runif(n = 2 * numPoints), nrow = numPoints, ncol = 2)
    # count points that fell inside the disk for each sample size up to
    # numPoints, and
    M <- cumsum((U[, 1] - C[1])^2 + (U[, 2] - C[2])^2 < r^2)
    # compute the relative error = (estimate - true value)/ true value.
    onlineEstimate <- M/(1:numPoints) * (l^2)/(r^2)
    onlineErrors[, i] <- (onlineEstimate - pi)/pi
}
# To plot the results nicely using ggplot2, first reshape the results into a
# data frame.
library(reshape2)
onlineErrors <- melt(onlineErrors)
names(onlineErrors) <- c("samples", "run", "error")
# Plot error vs number of samples for each experiment.
graph_error <- ggplot(onlineErrors, aes(x = samples, y = error, colour = factor(run), 
    group = factor(run)))
graph_error <- graph_error + geom_line(alpha = 1)
graph_error <- graph_error + theme(legend.position = "none") + ylim(-0.3, 0.3)
print(graph_error)