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)