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