Monte Carlo approximation of pi with a sphere

Question

I'm trying to estimate pi by uniform random sampling of points (x,y) inside a circle of radius 1 and then computing the corresponding z value in a sphere. Actually it is only a quarter-circle to simplify the computation. Then I compute the average z (which should be approximately 1/8 of the volume of the whole sphere) and then I compare it to the volume of a cube 1x1x1. It should be approximately 1/6 pi, but for some reason it isn't.

This is my matlab code:

r = rand(1000,1);
theta = rand(1000, 1) * pi/2;
x = zeros(1000);
y = zeros(1000);
z = zeros(1000);
for i = 1:1000
    x(i) = r(i)^(0.5) * sin(theta(i));
    y(i) = r(i)^(0.5) * cos(theta(i));
    z(i) = (1.0 - x(i) * x(i) - y(i) * y(i))^0.5;
end

mean(z)(1) * 6

It keeps saying that pi is approximately 4, which is nonsense - even if I increase the number of samples. Can you please explain me, where is the problem regardless of the fact, that I'm using pi to determine the angle when sampling random points inside a circle?

Answer 1

Aside from the fact that you need pi to call the radian-based trigonometric functions, your coordinates are incorrect as well.

In order to parametrize a sphere, you need two angles: theta and phi. You can also spare the loop using element-wise array operations:

N = 1000;
theta = 90*rand(N,1);
phi = 90*rand(N,1);
r = rand(N,1);

%hide the hiding of pi
%%hide pi:
%theta = deg2rad(theta);
%phi = deg2rad(phi);

%x = r.*sind(theta).*cosd(phi); %needless
%y = r.*sind(theta).*sind(phi); %needless
z = r.*cosd(theta);

pi_approx_maybe=mean(z(:))*6;

If z were an array (which it would if theta and phi came from a meshgrid rather than rand), then you'd need z(:) to get an array-wide mean, otherwise the result would be a vector. It's also clear that for the z component you don't even need the azimuthal angle (phi), only the polar angle (theta).