Estimate Pi with Monte Carlo

One method to estimate the value of \(\pi\) is by applying the Monte Carlo method. Let’s consider a circle inscribed in a square. Then for a radius \(r\) we have:

Area of the Circle = \(\pi r^2\)

Are of the Square = \(4 r^2\)

Hence, we can say that \(\frac{AreaCircle}{AreaSquare}=Prob=\frac{\pi}{4}\). This implies that \(\pi = 4Prob\)

We know that the Equation of Circle with center (j,k) and radius (r) is:

\( (x_1-j)^2+(x_2-k)^2 = r^2\) where in our case the center is (0,0) and the radius is 1. This means, that in order to calculate the Area of the Circle we need to consider all the points from the uniform distribution which have a distance from the center (0,0) less than or equal to 1, i.e:

\(\sqrt{ x_1^2+x_2^2} \leq 1\)

Notice that since we chose the center to be (0,0) and the radius to be equal to 1, it means that the coordinates of the square will be from -1 to 1. For that reason, we simulate from the uniform distribution with min and max values the -1 and 1 respectively.

Example Using R

Below, we represent how we can apply the Monte Carlo Method in R to estimate the pi.

# set the seed for reproducility
set.seed(5)

# number of simulations
n=1000000

# generate the x1 and x2 co-ordinates from uniform 
# distribution taking values from -1 to 1
x1<-runif(n, min=-1, max=1)
x2<-runif(n, min=-1, max=1)

# Distrance of the points from the center (0,0)
z<-sqrt(x1^2+x2^2)

# the area within the circle is all the z
# which are smaller than the radius^2
# in our case radius=1
4*sum((z<=1))/length(z)

[1] 3.14204

Let’s have a look also at the simulated data points.

InOut<-as.factor(ifelse(z<=1, "In", "Out"))

plot(x1,x2, col=InOut, main="Estimate PI with Monte Carlo")

As we can see, with 1M simulations we estimated the PI to be equal to 3.14204. If want to get a more precise estimation we can increase the number of simulations.

Example Using Python

import pandas as pd
import numpy as np
import math
import seaborn as sns
%matplotlib inline

# number of simulations
n=1000000

x1=np.random.uniform(low=-1,high=1, size=n)
x2=np.random.uniform(low=-1,high=1, size=n)

z=np.sqrt((x1**2)+(x2**2))

print(4*sum(z<=1)/n)

3.143288

Let’s have a look also at the simulated data points.

inout=np.where(z<=1,'in','out')

sns.scatterplot(x1,x2,hue=inout,legend=False)