Suppose we did not know the value of $\pi$ and we wanted to estimate its value using Monte Carlo methods. One practical approach is to draw a square of unit side, with its diagonal opposite corners extending from the coordinates origin $(0,0)$ to $(1,1)$. Now we try to simulate uniform random points from inside of this square by generating uniform random points along the $X$ and $Y$ axes, i.e., by generating two random uniform numbers (x,y) from the range $[0,1]$.

Now the generated random point $P$ has the coordinate $(x,y)$, so we can calculate its distance from the coordinate origin. Now suppose we also draw a quarter-circle inside of this square whose radius is unit and is centered at the origin $(0,0)$. The ratio of the area of this quarter-circle, $S_C$ to the area of the square enclosing it, $S_S$ is,

$\frac{S_C}{S_S} = \frac{\frac{1}{4}\pi r^2}{r^2} = \frac{1}{4}\pi$

This is because the area of the square of unit-length sides is just 1. Therefore, if we can somehow measure the area of the quarter $S_C$, then we can use the following equation, to get an estimate of $\pi$,

$\pi = 4S_C$

To obtain, $S_C$, we are going to throw random points in the square, just as described above, and then find the fraction of points, $f=n_C/n_{\rm total}$, that fall inside this quarter-circle. This fraction is related to the area of the circle and square by the following equation,

$f=\frac{n_C}{n_{\rm total}} = \frac{S_C}{S_S}$

Therefore, one can obtain an estimate of $\pi$ using this fraction,

$\pi \approx 4\frac{n_C}{n_{\rm total}}$

Now, write a script that takes in the number of points to be simulated, and then calculates an approximate value for $\pi$ based on the Monte Carlo algorithm described above. Write a second function that plot the estimate of $\pi$ versus the number of points stimulated, like the following,

Here is an example MATLAB script that estimates $\pi$: estimatePi.m,

nsample = 1000000;
X = rand(nsample,2);
InsideCircle = ( X(:,1).^2 + X(:,2).^2 ) <= 1;
ApproximatePI = zeros( size(InsideCircle) );
sumInsideCircle = 0.0;
for isample = 1:nsample
sumInsideCircle = sumInsideCircle + InsideCircle(isample) ;
ApproximatePI(isample) = 4.0 * sumInsideCircle / isample ;
end
disp( [ 'Approximate value of PI for '...
, num2str(nsample) ...
, ' samples: ' ...
, num2str(ApproximatePI(end)) ...
] ...
);
figure(); hold on; box on;
line( [ 1 nsample ] ...
, [ pi pi ] ...
, 'linestyle' , '--' ...
, 'linewidth' , 3 ...
, 'color' , 'black' ...
);
plot( ApproximatePI(2:end) ...
, 'linewidth' , 2 ...
, 'color' , 'red' ...
);
xlabel( 'Number of Sampled Points' ) ;
ylabel( 'Approximate Value of Pi' ) ;
set( gca , 'XScale' , 'log' );
saveas( gcf, 'approximatePI.png' );


Here is an example Python code that estimates $\pi$,

import numpy.random as rnd
import numpy.linalg as nlg
import numpy as np
import matplotlib.pyplot as plt

def estimatePi(n=100000):
x = np.zeros((n,4),dtype=np.dtype)
counter = 0
for i in range(n):
x[i,0:2] = rnd.uniform(0.0,1.0,2)  # generate two uniform random numbers
x[i,2] = nlg.norm(x[i,:], ord=2)
if (x[i,2]<=1.0): counter += 1
x[i,3] = 4.0*np.double(counter)/np.double(i+1)
return x # the running approximate of pi is returned as a vector in the fourth column of x

def plotPi(n=10000):
x = estimatePi(n)
plt.semilogx( list(range(1,n+1)) \
, x[:,3] \
) # plot with color red, as line
plt.xlabel('Number of simulated points')
plt.ylabel('Approximate value of Pi')
#plt.axis([1, n , 0.1, 4.0]) # [xmin, xmax, ymin, ymax]
plt.title('Estimating Pi by Monte Carlo simulation')
plt.savefig('approximatePi_{}.png'.format(n))
plt.show()

plotPi()