I'm trying to estimate pi using the monte carlo method and the area of a circle (see code for reasoning and approach).
I'm getting outputs that seem to converge to pi (see command window screenshot)
But I'm not quite satisfied with my visuals. I think they're illustrating what's happening, but I'd like to make the points that fall inside of the circle a different color than the points that fall outside of the circle. I think the problem may be in how I'm identifying the points that fall inside the circle in my 'for' loops, but I'm not sure how to go about the solution with an 'if else' statement.
Any and all help would be appreciated!
(Just to note: I have the barebones, student version of Matlab. If your suggestion includes using a certain toolbox, please assume I don't have it)
%% Pi Estimation using Monte Carlo (Area of a Circle)
% Reasoning and Approach
% Imagine a circle that exists inside of a square, such that each side of
% the square touches tangent to an edge of the circle. Now, the circle is
% going to have some nonzero radius, denoted as 'r'. Because the circle is
% perfectly enclosed within this square, then by extension, each side of the
% square is going to be the same length as the diameter of the circle, or
% '2r'. Now, the area of a circle is given by pi*r^2. Then it stands to
% reason that the area of the square that encloses the circle is given by
% 2r*2r, or 4r^2 (area is given by length*width). My goal is to implement
% the Monte Carlo Algorithm to generate a set of random points within the
% square, and use the number of points that end up falling within the
% circle, divided by the total number of points generated, and multiplied
% by 4 will give me an approximation of the value of pi.
pi_ref=3.141592; %this is pi to 6 decimal places. I'm going to use this as
%my 'accepted' value of pi so that I can make error calculations with what ke error calculations with what
%the algorithm spits out and see how it changes as the number of randomly
%generated points increases.
%using up to 100 points to estimate
numpoint1=100; %this is the total number of points we want to generate.
%these points are going to fall either inside of the circle or inside of
%the square. we'll identify the ones that fell inside of the circle in the
%for loop
pi_est1=zeros(); %pre-loaded with zeros for speed
circ_stdev1=zeros(); %standard deviation array. pre-loaded with zeros for speed
for i=1:numpoint1 %start by using 1 point and work up to the designated amount of points to be generated
x1=rand(i,1); %'i' row, 1 column vector of random x values will be generated with each iteration of the for loop
y1=rand(i,1); %same thing for y
check_in1=(x1-0.5).^2 + (y1-0.5).^2 < 0.5^2; %this is how we'll isolate the points that fell inside of the circle
pi_est1(i)=sum(check_in1)*4/i; %calculation of pi. sum up the number of points within the circle,
%multiply by 4, and then divide by the number of points generated
%during THAT particular iteration
circ_stdev1(i)=std(pi_est1); %calculate the standard deviation each time pi is calculated
endcirc_avg1=mean(pi_est1);circ_error1=abs((circ_avg1 - pi_ref)/pi_ref)*100; %calculate the percentage
%error in the calculation of pi using 'n' points
%using up to 1000 points to estimate
numpoint2=1000; %this is the total number of points we want to generate.%these points are going to fall either inside of the circle or inside of%the square. we'll identify the ones that fell inside of the circle in the%for looppi_est2=zeros(); %pre-loaded with zeros for speedcirc_stdev2=zeros(); %standard deviation array. pre-loaded with zeros for speedfor i=1:numpoint2 %start by using 1 point and work up to the designated amount of points to be generated x2=rand(i,1); %'i' row, 1 column vector of random x values will be generated with each iteration of the for loop y2=rand(i,1); %same thing for y check_in2=(x2-0.5).^2 + (y2-0.5).^2 < 0.5^2; %this is how we'll isolate the points that fell inside of the circle pi_est2(i)=sum(check_in2)*4/i; %calculation of pi. sum up the number of points within the circle, %multiply by 4, and then divide by the number of points generated %during THAT particular iteration circ_stdev2(i)=std(pi_est2); %calculate the standard deviation each time pi is calculatedendcirc_avg2=mean(pi_est2);circ_error2=abs((circ_avg2 - pi_ref)/pi_ref)*100; %calculate the percentage %error in the calculation of pi using 'n' points%using up to 10000 points to estimate
numpoint3=10000; %this is the total number of points we want to generate.%these points are going to fall either inside of the circle or inside of%the square. we'll identify the ones that fell inside of the circle in the%for looppi_est3=zeros(); %pre-loaded with zeros for speedcirc_stdev3=zeros(); %standard deviation array. pre-loaded with zeros for speedfor i=1:numpoint3 %start by using 1 point and work up to the designated amount of points to be generated x3=rand(i,1); %'i' row, 1 column vector of random x values will be generated with each iteration of the for loop y3=rand(i,1); %same thing for y check_in3=(x3-0.5).^2 + (y3-0.5).^2 < 0.5^2; %this is how we'll isolate the points that fell inside of the circle pi_est3(i)=sum(check_in3)*4/i; %calculation of pi. sum up the number of points within the circle, %multiply by 4, and then divide by the number of points generated %during THAT particular iteration circ_stdev3(i)=std(pi_est3); %calculate the standard deviation each time pi is calculatedendcirc_avg3=mean(pi_est3);circ_error3=abs((circ_avg3 - pi_ref)/pi_ref)*100; %calculate the percentage %error in the calculation of pi using 'n' points%display the results in the command window
fprintf('Monte Carlo (Circle Area) Results:\n ');fprintf('Using %d points: %f \n Percent Error: %3.2f (compared to accepted value of pi) \n ',numpoint1,circ_avg1,circ_error1);fprintf('\n Using %d points: %f \n Percent Error: %3.2f (compared to accepted value of pi) \n ',numpoint2,circ_avg2,circ_error2);fprintf('\n Using %d points: %f \n Percent Error: %3.2f (compared to accepted value of pi)\n \n \n',numpoint3,circ_avg3,circ_error3);%illustrate the results
figure(4); subplot(3,1,1); plot(pi_est1,'.b'); title('Fig. 4: Estimating Pi with Monte Carlo (Area of a Circle)'); grid on;subplot(3,1,2); plot(pi_est2,'.r'); ylabel('Expected Value'); grid on;subplot(3,1,3); plot(pi_est3,'.k'); xlabel('Number of points'); grid on;figure(5);subplot(3,1,1); plot(circ_stdev1,'b','linewidth',1.5); title('Fig. 5: Standard Deviation (corresponds to Fig. 4)'); grid on;subplot(3,1,2); plot(circ_stdev2,'r','linewidth',1.5); ylabel('Standard Deviation'); grid on;subplot(3,1,3); plot(circ_stdev3,'k','linewidth',1.5); xlabel('Number of points'); grid on;%Below is the code to generate 'snapshots' of what is happening, to aid
%visualization of the circle inside of the square and the population of
%random points.
theta=0:0.01:2*pi; %set of angles for the circle
circlex=0.5+0.5*cos(theta); %function to create the x values of the circle.
%circle is centered at (0.5,0.5) and has a radius of 0.5
circley=0.5+0.5*sin(theta); %function to create y values of circle
%using 100 points to estimate
popnum1=100; %number of random points we want to generate
xvec1=rand(popnum1,1); yvec1=rand(popnum1,1); %vectors containing 'n' randomly
%generated numbers in a single column.
figure(6); subplot(2,2,1); plot(xvec1,yvec1,'.'); hold on; %populate the square with random points
plot(circlex,circley,'r','linewidth',1.5); hold off; %plot the circle inside of the
%square. each side of the square should be tangential to the edge of the
%circle
title('Fig. 6: Visualization: Using Circle Area to Estimate Pi, number of points = 100');%using 500 points to estimate
popnum2=500; %number of random points we want to generatexvec2=rand(popnum2,1); yvec2=rand(popnum2,1); %vectors containing 'n' randomly %generated numbers in a single column.subplot(2,2,2); plot(xvec2,yvec2,'.'); hold on; %populate the square with random pointsplot(circlex,circley,'r','linewidth',1.5); hold off; %plot the circle inside of the %square. each side of the square should be tangential to the edge of the%circletitle('Number of points = 500');%using 1000 points to estimate
popnum3=1000; %number of random points we want to generatexvec3=rand(popnum3,1); yvec3=rand(popnum3,1); %vectors containing 'n' randomly %generated numbers in a single column.subplot(2,2,3); plot(xvec3,yvec3,'.'); hold on; %populate the square with random pointsplot(circlex,circley,'r','linewidth',1.5); hold off; %plot the circle inside of the %square. each side of the square should be tangential to the edge of the%circletitle('Number of points = 1,000');%using 10000 points to estimate
popnum4=10000; %number of random points we want to generatexvec4=rand(popnum4,1); yvec4=rand(popnum4,1); %vectors containing 'n' randomly %generated numbers in a single column.subplot(2,2,4); plot(xvec4,yvec4,'.'); hold on; %populate the square with random pointsplot(circlex,circley,'r','linewidth',1.5); hold off; %plot the circle inside of the %square. each side of the square should be tangential to the edge of the%circletitle('Number of points = 10,000');
Best Answer