Monte Carlo integration (hit or miss) to find the area of a circle of radius R

Hi everyone. Trying to solve an old exam topic regarding Monte Carlo integration, I wrote the following code for which I based it on a code from a professor in the C language. For R = 1 it worked so I assume I did it correctly. Then I tried R > 1 and the results were wrong compared to the formula πR^2 so I tried a few modifications by trial and error method. I fixed the problem and it seems to works correctly but there is a modification that I don't understand. If you could help me understand this it would be great. Thank you in advance.

clear, clc, format short

%---------- Computation for R = 1 ----------
R = 1;
N = 1e+7;
insideCircle = 0;

for i = 1:N
    
    x = rand(1);
    y = rand(1);
    
    if (x^2 + y^2) < R
        insideCircle = insideCircle + 1;
    end
end

Area_MC = 4*insideCircle/N;
Area = pi*R^2;
error = 100*abs(Area - Area_MC)/Area;

fprintf('\n\tExact area\t\t%f\n\tMonte Carlo area\t%f\n\tPercentage error\t%.5f%%\n',Area,Area_MC,error);

	Exact area		3.141593
	Monte Carlo area	3.142050
	Percentage error	0.01456%

%---------- Computation for R > 1 (e.g. R = 4 but you can try any real value) ----------
R = 4; %input('Enter radius R: R = ','s'); R = str2double(R);

% while R <= 0 || isreal(R) == 0 || isnan(R) == 1 || isnumeric(R) == 0
%     disp('Invalid input')
%     R = input('Enter radius R: R = ','s');
%     R = str2double(R);
% end

N = 1e+7;
insideCircle = 0;

for i = 1:N
    
    x = rand(1)*R;
    y = rand(1)*R;
    
    if (x^2 + y^2) < R
        insideCircle = insideCircle + 1;
    end
end

Area_MC = 4*R^3*insideCircle/N; % Modification needed that I don't understand (I had to multiply with R^3)
Area = pi*R^2;
error = 100*abs(Area - Area_MC)/Area;

fprintf('\n\tExact area\t\t%f\n\tMonte Carlo area\t%f\n\tPercentage error\t%.5f%%\n',Area,Area_MC,error);

	Exact area		50.265482
	Monte Carlo area	50.285594
	Percentage error	0.04001%

Answer 1

[Torsten]

You compare the area of the square with side length R and corner points (0,0), (R,0), (R,R) and (0,R) with the area of the quarter circle of radius R. This square has area R^2, and x and y are generated to cover this square uniformly.

N = 1e+7;
insideCircle = 0;

for i = 1:N
    x = R*rand();
    y = R*rand();
    
    if (x^2 + y^2) < R^2
        insideCircle = insideCircle + 1;
    end
end

Area_MC = 4*R^2*insideCircle/N; % Modification needed that I don't understand (I had to multiply with R^3)
Area = pi*R^2;
error = 100*abs(Area - Area_MC)/Area;

fprintf('\n\tExact area\t\t%f\n\tMonte Carlo area\t%f\n\tPercentage error\t%.5f%%\n',Area,Area_MC,error);

[Left Terry]

How do I plot this in order to check if I have uniform distribution of points?

[Torsten]

R = 4;
N = 5e+3;
insideCircle = 0;

for i = 1:N
    x = R*rand();
    y = R*rand();    
    if (x^2 + y^2) < R^2
        insideCircle = insideCircle + 1;
        xinsideCircle(insideCircle) = x;
        yinsideCircle(insideCircle) = y;
    end
end
Area_MC = 4*R^2*insideCircle/N; % Modification needed that I don't understand (I had to multiply with R^3)
Area = pi*R^2;
error = 100*abs(Area - Area_MC)/Area;
fprintf('\n\tExact area\t\t%f\n\tMonte Carlo area\t%f\n\tPercentage error\t%.5f%%\n',Area,Area_MC,error);

	Exact area		50.265482
	Monte Carlo area	50.496000
	Percentage error	0.45860%

plot(xinsideCircle,yinsideCircle,'x')

●

[Left Terry]

So, the mistakes i made were:

1) I compared x^2 + y^2 with R instead of R^2

2) The x and y range was [0,R] instead of [-R,R] and because of that I didn't have uniform distribution in [-R,R] but I had uniform distribution in [0,R].

Am I correct ?

[Torsten]

Thus the only problem with your code was that you had to use

if (x^2 + y^2) < R^2    

instead of

if (x^2 + y^2) < R

[Left Terry]

By the way I don't see any difference between rand() and rand(1) in the plots

clear, clc, format short

R = 2;
N = 2e+3;

%---------- Using rand(1) ----------
insideCircle = 0;

for i = 1:N
    
    x = (-1 + 2*rand(1))*R;
    y = (-1 + 2*rand(1))*R;
    
    if (x^2 + y^2) < R^2
         
        insideCircle = insideCircle + 1;
        subplot(1,2,1)
        plot(x,y,'ro','MarkerFaceColor','r'), hold on
       
    else
        subplot(1,2,1)
        plot(x,y,'bo','MarkerFaceColor','b'), hold on
        title('rand(1)')
    end
end

Area_MC = 4*R^2*insideCircle/N;
Area = pi*R^2;
error = 100*abs(Area - Area_MC)/Area;
fprintf('\n\tUsing rand(1)\n\tExact area\t\t\t%f\n\tMonte Carlo area\t%f\n\tPercentage error\t%.5f%%\n',Area,Area_MC,error);

	Using rand(1)
	Exact area			12.566371
	Monte Carlo area	12.568000
	Percentage error	0.01297%

%---------- Using rand() ----------
insideCircle = 0;
for i = 1:N
        
    X = (-1 + 2*rand())*R;
    Y = (-1 + 2*rand())*R;
    
    if (X^2 + Y^2) < R^2
         
        insideCircle = insideCircle + 1;
        subplot(1,2,2)
        plot(X,Y,'go','MarkerFaceColor','g'), hold on
    else
        
        subplot(1,2,2)
        plot(X,Y,'bo','MarkerFaceColor','b'), hold on
        title('rand()')
    end
end

Area_MC = 4*R^2*insideCircle/N;
Area = pi*R^2;
error = 100*abs(Area - Area_MC)/Area;
fprintf('\n\tUsing rand()\n\tExact area\t\t\t%f\n\tMonte Carlo area\t%f\n\tPercentage error\t%.5f%%\n',Area,Area_MC,error);

	Using rand()
	Exact area			12.566371
	Monte Carlo area	12.608000
	Percentage error	0.33128%

[Torsten]

I didn't realize at first that you are working in a quarter circle instead of a full circle.

For this approach, the only thing that was wrong in your original code is that you used "if (x^2 + y^2) < R" instead of "if (x^2 + y^2) < R^2".

Of course, "Area_MC = 4*R^3*insideCircle/N;" has to be reset to "Area_MC = 4*R^2*insideCircle/N;"

[Walter Roberson]

By the way I don't see any difference between rand() and rand(1) in the plots

X = rand returns a random scalar drawn from the uniform distribution in the interval (0,1).

X = rand(n) returns an n-by-n matrix of uniformly distributed random numbers.

... so when n = 1, rand(1) returns a 1 x 1 matrix of uniformly distributed random numbers, which is the same as returning a random scalar. rand() and rand(1) mean the same thing.