One method to estimate the value of \(\pi\) is by applying the **Monte Carlo method**. Let’s consider a circle inscribed in square. Then for a radious \(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 co-ordinates 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)