Evaluate the area of a circle of radius $1= \pi$ using Monte Carlo method . Hence we can generate pairs of random numbers $(x_i,y_i) \in [-1,1]$.
Thus :
$$ \pi= \frac {Number Of Samples Inside The Circle}{Total Number Of Samples} X 4 $$
Consider a unit circle inscribed in a square, each of the small circles drawn on this figure represents a random point that was generated in the square, the red and blue circles represent points inside and outside the unit circle respectively (I got the figure from google) . If we choose a point uniformaly at random within the square, then the probability that the point lies in the unit circle is
$$ \pi= \frac {Number Of Samples Inside The Circle}{Total Number Of Samples} $$
We know that the area of the circumscribed square is $4$, if we knew $p$, then we could compute the area of the unit circle:
Area of the unit circle = $p$ x area of the circumscribed square = $4p$
$\endgroup$ 3How to Estimate $\pi$ using the Monte Carlo Method in MATLAB, then what's the error a quarter of a circle ?
2 Answers
$\begingroup$The following code is vectorized, and thus typically faster than using loops:
N = 1e6; % Number of samples
r = 1; % Circle radius
x = 2*rand(1,N)-1; % N samples between -1 and 1
y = 2*rand(1,N)-1; % N samples between -1 and 1
r2 = x.^2 + y.^2; % compute squared distance to origin
area = mean(r2<=r^2)*4; % the area is 4 (area of square) times the proportion % of points in the circle. Note that the threshold % should be r^2 (1^2)
error = pi - area % error in the estimation of piN_all = round(10.^(1:.5:7)); % Numbers of samples
r = 1; % Circle radius
errors = NaN(size(N_all)); % Preallocate result
for k = 1:numel(N_all) % Loop over size of N_all N = N_all(k); % Current N x = 2*rand(1,N)-1; % N samples between -1 and 1 y = 2*rand(1,N)-1; % N samples between -1 and 1 r2 = x.^2 + y.^2; % compute squared distance to origin area = mean(r2<=r^2)*4; % the area is 4 (area of square) times the % proportion of points in the circle error_all(k) = pi-area; % error in the estimation of pi
end % End loop
semilogx(N_all, error_all, 'o-'); % Plot with logarithmic x axis
grid % Use gridSome code like this may work:
clear;
N=1000; # the experiment event number
r=1; #the circle radius
n=0; # sucessful event number
for i=1:N x=-r+2*r*rand(); y=-r+2*r*rand(); if ((x^2+y^2)<=r^2) n = n+1; end
end
pi_sim=4*n/N
