0

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?

Rafael Korbas
  • 2,213
  • 3
  • 19
  • 30
  • The mean would be `1/6 pi` if you let `x` and `y` vary over the entire unit *square* `[0, 1] x [0, 1]`. As it is, you'd expect the mean to be `(1/6 pi) / (1/4 pi) = 2/3`. – Mark Dickinson Nov 22 '15 at 21:23
  • ok, thanks, I see. But I would like to avoid rejection sampling on a square, which in higher dimensions would be even more inefficient. Do you see any way to fix it? – Rafael Korbas Nov 22 '15 at 21:28
  • I don't really understand what you're trying to achieve. Why can't you do this using a circle instead of a sphere? (E.g., sample `x` uniformly on `[0, 1]` and look at the average value of `sqrt(1 - x*x)`.) – Mark Dickinson Nov 22 '15 at 21:31
  • Yes, I tried it with a circle and it worked, but I wanted to generalize it to higher dimensions. – Rafael Korbas Nov 22 '15 at 21:32
  • 1. `mean(z)(1)` is not valid matlab. Are you using octave? 2. I suggest avoiding magic numbers, set `N=1000` in the beginning and use `N` each time later 3. Use array operations instead of a loop: e.g. `z=(1-x.^2-y.^2).^(0.5)` (more efficient, same result). **4. In a sphere you need theta and phi, and `x=r.*sin(theta).*cos(phi); y=r.*sin(theta).*sin(phi); z=r.*cos(theta)`**. In your current version `x.^2+y.^2==r` for all points, so you're just averaging `sqrt(1-r)` which doesn't make sense. – Andras Deak -- Слава Україні Nov 22 '15 at 22:24
  • 1
    Also, you are using `pi` in a function that calculates `pi`, which is not very helpful! – David Nov 22 '15 at 22:28
  • @David's comment shows very well how your approach is inherently flawed. As Mark Dickinson suggested, you should generate points on a cube and count the fraction of points inside a unit sphere. Then use the formula of an n-dimensional sphere to approximate pi. – Andras Deak -- Слава Україні Nov 22 '15 at 22:38
  • @AndrasDeak Yes I know, that it isn't very helpful to use pi when approximating pi itselft, but it is given by the fact that the sin and cos get the parameters in radians, if they got them in degrees, it would not be an issue. – Rafael Korbas Nov 22 '15 at 22:43
  • Then use `dec2rad` and hide pi;) And try the boldface part of my earlier comment. – Andras Deak -- Слава Україні Nov 22 '15 at 22:52
  • @AndrasDeak maybe using approximation formulas, like Taylor series or there are other methods, like Bhaskara's formula. – Rafael Korbas Nov 22 '15 at 22:55
  • @AndrasDeak yes, thanks, I'm trying to figure out, how it works, because it keeps giving me weird results, when I'm trying the boldfaced part of your previous comment. – Rafael Korbas Nov 22 '15 at 23:00
  • I added an answer clarifying my point above. I'm highly uncertain about your formula from the `mean`, though. For instance, you're comparing a volume to a single component in your question. So have you performed the spherical integral necessary to compute `` on a sphere? Are you sure it's `1/6pi` times the size of a cube? – Andras Deak -- Слава Україні Nov 22 '15 at 23:12
  • @AndrasDeak Thanks for your answer. Probably I don't understand what you mean. Actually, I'm not very familiar with spherical integrals and in generally, integrals in more than 2 dimensions. But I expected, that it should be something like the average z coordinate if I pick uniformly randomly points from a circle at the x,y plane. And yes, it could be done with a cube and rejection sampling, but I wanted to try a different approach. – Rafael Korbas Nov 22 '15 at 23:22
  • Well, I found a new problem. Your points are not uniform on the unit sphere (another reason why you should use rejection sampling on a cube). And even if they were, you can't spare doing the integrals by hand first. The procedure would be this: compute on paper, approximate from Monte Carlo, then compare the two. – Andras Deak -- Слава Україні Nov 22 '15 at 23:27
  • Rafael, I keep running into problems. [This page](http://mathworld.wolfram.com/SpherePointPicking.html) discusses how to choose points uniformly *on a spherical surface*, but then you'd need `pi` to say anything about the result. I don't really see right now how we could salvage this approach... – Andras Deak -- Слава Україні Nov 22 '15 at 23:42
  • @Andras Deak yes, I see too many problems, too, but thank for your effort :) Also I finally figured out that in my original approach, if I want to get an approximation of the volume of the sphere, I need to consider as the base the quarter circle, over which I was trying to compute the average height, not the square, and of course, for that thing, I need pi :( – Rafael Korbas Nov 23 '15 at 08:36
  • Yes, exactly, that's what I realized when I tried to integrate over the sphere... A very traditional way to approximate pi from Monte Carlo (actually, one of the first Monte Caro methods) is [Buffon's needle problem](http://mathworld.wolfram.com/BuffonsNeedleProblem.html). – Andras Deak -- Слава Україні Nov 23 '15 at 11:10

2 Answers2

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).

Andras Deak -- Слава Україні
  • 33,737
  • 11
  • 83
  • 111
  • Using `deg2rad` still uses `pi`! – David Nov 22 '15 at 23:15
  • @David yes, yes, it does, but that's not *inherently* a flaw. I initially thought so too, but the Taylor series of sin/cos don't involve pi, so nothing stops you from computing those in degrees. It just happens that those functions are defined with radian inputs. – Andras Deak -- Слава Україні Nov 22 '15 at 23:16
  • 1
    @David I removed the radians from the picture, see my edit:P Although I'm still a bit uncertain. The fact that the Taylor series involves x^k/k! is probably related to radians. But you probably could formulate it without pi, just relating it to the number of `arbitrary degree`s in a full circle... – Andras Deak -- Слава Україні Nov 22 '15 at 23:18
  • :P :P :P :P I see what you did there!! – David Nov 22 '15 at 23:20
1

This is not at all how to compute pi without using it in the 1st place.

Actually this is not how MC works: it needs to fail sometimes, and the ratio of success/fail gives rise to a value. otherwise if you only have successes, then you must have 'baked' the searched value in the computation somehow...

Here's how to do it (bear with me, my Matlab is rusty) by sampling the unit cube's volume:

shoots = 100000;
hits = 0;
for i = 1:shoots
    % Sample the whole unit cube
    x = rand();
    y = rand(); 
    z = rand();
    % Are we inside the sphere?
    r = x^2 + y^2 + z^2;
    if (r < 1)
        hits = hits + 1; % We're inside! it's a hit
    end;
end
ratio = hits/shoots; % This ratio is statistically ~ Volume(quarter sphere) / Volume(cube)
myPi = ratio * 3;

Now if you really wanted to use the sphere's surface, then you need to sample it without knowing it's a sphere.

So for instance:

  • sample the unit cube's three outer faces (the three that do not go through the origin)
  • normalize() that point cloud to bring it onto the desired sphere surface
  • Then, triangulate that mesh (ouch)
  • Compute the triangulation's total surface (ouch again) ...and you have an approximation of the quarter sphere's surface (pi*r²), so that get you pi directly!
MrBrushy
  • 680
  • 3
  • 17