basic monte carlo question: area of circle

Question: Use monte carlo method to find area of circle radius 5 inside 7x7 square.

So this is the code I've put together so far to determine how many points land inside or on the circle (hits). My output for hits is 0 but I can't for the life of me figure out why, to me everything seems fine but clearly isn't.

clear;
N= 1000; % number of points generated
a = -7;
b = 7;
hits= 0;
for k = 1:N
     x = a + (b-a).*rand(N,1);
     y = a + (b-a).*rand(N,1);
     % Count it if it is in/on the circle
     radii = sqrt(x.^2+y.^2);
         if radii <=5;
            hits = hits + 1;
         end
end
disp(hits)

Answer 1

[pfb]

You are mixing things up. Your code is part vectorized, part scalar.

You have a loop over N, but then at each iteration you generate N points (so you're generating N^2 points overall). The variables x, y, and radii are vectors.

Note that

 [6 7 4 2 8 5] <=5

gives

 [0 0 1 1 0 1]

Therefore condition radii<=5 has a very low probability to be met ( all of the N radii should be less than 5).

As a result hits is (almost) never accumulated.

Your code is structured for a scalar radius. It would work if you substitute the lines

 x = a + (b-a).*rand(N,1);
 y = a + (b-a).*rand(N,1);

with

 x = a + (b-a).*rand;
 y = a + (b-a).*rand;

In this case hits will be a number between 0 and N, and hits/N is going to be proportional to the ratio of the areas of the circle and square, pi*(5/14)^2 (note that you can use simple power operator ^ instead of elementwise .^)

You can avoid the loop altogether using the vectorized variables.

 N= 1000; % number of points generated
 a = -7;
 b = 7;
 x = a + (b-a).*rand(N,1);
 y = a + (b-a).*rand(N,1);
 radii = sqrt(x.^2+y.^2);
 hits = sum(radii<=5);

Answer 2

[Eric]

Thanks both really helped a bunch.

How would I go about plotting the graph so that the hits are one colour and the rest are another?

I tried the following but I just get a linear line.

x = linspace(-7,7);
y = linspace(-7,7);
figure(1);
plot(x,y);
xlabel('x');
ylabel('y');

axis square
grid on;
title('Area of Circle');

[pfb]

Yes, your code above clearly gives you a straight line bisecting the 1st and 3rd quadrant: x and y are the same vector, so you're plotting the graph of the function y=f(x)=x.

So, if you want to plot green dots for the hits and red dots for the misses, you first build your random points

 N= 1000; % number of points generated
 a = -7;
 b = 7;
 x = a + (b-a).*rand(N,1);
 y = a + (b-a).*rand(N,1);

then build a (logical) vector

 i = radii<=5

that is 1 for the hits and 0 for the misses (this is what you used above for counting the hits). Then, exploiting logical indexing

 plot(x(i),y(i),'.g');
 plot(x(~i),y(~i),'.r');

Where ~i (not i) is a vector that is 1 for the misses and 0 for the hits.

You can also do this in a single plot, using scatter

 scatter(x,y,10,i,'filled')

but then you have to hack the colors and possibly other things.

All in all, I'd go with the two separate plots.

[Eric]

Thanks for helping me out pfb.

I've tried incoroprating your hints into my code.

I get the error " Subscript indices must either be real positive integers or logicals."

I'm not sure if im building the logical vector right. Instead of i, im trying to use hits.

clear;
N= 10000; % number of points generated
a = -7;
b = 7;
r=5;
hits= 0;
misses = 0;
for k = 1:N 
x = a + (b-a).*rand(N,1);
y = a + (b-a).*rand(N,1);
% Count it if it is in/on the circle
radii = sqrt(x.^2+y.^2);
if radii <=r;
    hits = hits + 1;
elseif radii > r
    misses = misses + 1;
end
end
disp('hits')
disp(hits)
disp('misses')
disp(misses)

x = linspace(-7,7);
y = linspace(-7,7);
figure(1);
plot(x(hits),y(hits),'.g');
plot(x(misses), y(misses), '.b');
xlabel('x');
ylabel('y');

[pfb]

No wait. You are mixing things again.

How many points do you need? I guess N=10000, not N^2.

If this is the case, why do you insist with the loop??? Forget the loop.

What you do in the last part also makes no sense. If you write

 x = linspace(-7,7);
 y = linspace(-7,7);

then your old x and y are canceled and replaced by two 1x100 vectors. You entirely loose the information you used to build hits.

Furthermore hits and misses are just numbers, they do not contain the information of which points hitted or missed, right?

For the plot, you need the data you used to build hits.

Finally, a suggestion: when you report an error, please include the relevant line (which matlab always tells you). Otherwise it can be hard to pinpoint the problem. In your case it is easy: x an y have 100 elements each (check the default usage of linspace), while hits and misses are quite likely larger than that. So x(hits) is going to give you the error.

The following is going to work.

N= 10000; % number of points generated (no loops!!)
a = -7;
b = 7;
r=5;
% both x and y have already 10000 elements. No need to loop!
x = a + (b-a).*rand(N,1);
y = a + (b-a).*rand(N,1);

% now radii is a 1x10000 vector of radii.
radii = sqrt(x.^2+y.^2);
% now i is a 1x10000 logical vector of "hits"
i = radii <= r;
% if you want to count the hits, simply do
hits = sum(i);
% misses are simply
misses = N-hits;

% plot the result
plot(x(i),y(i),'.g');
hold;
plot(x(~i),y(~i),'.r');

% you can put the info in the title
ttl = sprintf('%d trials, %d hits, %d misses',N,hits,misses);
title(ttl);

% or perhaps, even better, in the labels 
plot(x(i),y(i),'.g','DisplayName',sprintf('%d hits',hits));
hold;
plot(x(~i),y(~i),'.r','DisplayName',sprintf('%d misses',misses));
legend('show');
% and some info in the title
ttl = sprintf('%d trials, Area of the circle: %1.3f',N,hits/N*14^2);
title(ttl);

[Eric]

Sorry I should have pointed out, I used loops because we unfortunately aren't allowed to use the sum function as we haven't been taught it yet.

Ive been trying to use a loops to add up the hits.

[pfb]

ok, then you can write the function sum with a loop. After creating the vector i (which you need for the plotting)

 hits = 0;
 for j = 1:N
  hits=hits+i(j);
 end
 misses= N-hits;

If you're not allowed to use logical indexing in the plot you can pretend you do not know it and do this

 % vectors for the hits
 xh = zeros(1,hits);
 yh = xh;
 % vector for the misses
 xm = zeros(1,misses);
 ym = xm;

 % counters
 mc=1;
 hc=1;

 for j=1:N
  if(i(j))
    % put the point in the hits vector
    xh(hc)=x(j);
    yh(hc)=y(j);
    % increment the hit counter
    hc=hc+1;
  else 
    % put the point in the misses vector
    xm(mc)=x(j);
    ym(mc)=y(j);
    % increment the miss counter
    mc=mc+1;
 end
 end

 plot(xh,yh,'g.');
 hold
 plot(xm,ym,'r.');

etcetera

Answer 3

[Eric]

Thanks a lot, appreciate it.

[pfb]

you're welcome. Note that there was an error in my code above (in one instance I wrote hits instead of misses). Now it should be ok.

Answer 4

[munesh pal]

Calculate the area of the circle using Monte Carlo Simulation

clc
clear all
close all 
%%input circle radius and coordinate 
r=2;
c_x=7;
c_y=7;
%%position of the outer rectangle
p_x=c_x-r;
p_y=c_y-r;
pos=[p_x,p_y,r^2,r^2] 
%%random number generation within the sqaure
N=1000;
a=p_x;
b=p_x+(2*r);
x = a + (b-a).*rand(N,1);
y = a + (b-a).*rand(N,1);
radii = sqrt((x-c_x).^2+(y-c_y).^2);
i = radii <= r;
hits = sum(i);
misses = N-hits;
plot(x(i),y(i),'.g');
hold;
plot(x(~i),y(~i),'.r');
rectangle('Position',pos,'Curvature',[1 1])
rectangle('Position',pos,'EdgeColor','r')
actual_a=22/7*r^2;
ttl = sprintf('Estimate Actual Area of circle: %1.3f, Area of the circle: %1.3f',actual_a,hits/N*(2*r)^2);
title(ttl);
axis equal

[hosein bashi]

can you explain please why we use random numbers? why we don't assume a grid and put our numbers in the middle of the cells? it would be more uniform in this way.