Monte Carlo 圆周率的近似值

Question

我正在尝试通过对半径为 1 的圆内的点 (x,y) 进行均匀随机抽样来估计 pi，然后计算球体中相应的 z 值。实际上只是一个四分之一圆来简化计算。然后我计算平均 z（它应该大约是整个球体体积的 1/8），然后我将它与 1x1x1 立方体的体积进行比较。它应该大约是 1/6 pi，但由于某些原因它不是。

这是我的matlab代码：

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

它一直说pi大约是4，这是胡说八道-即使我增加了样本数。你能解释一下吗，不管事实如何，我在对圆内的随机点进行采样时使用 pi 来确定角度，问题出在哪里？

Answer 1

除了您需要 pi 来调用基于弧度的三角函数之外，您的坐标也不正确。

为了参数化球体，您需要两个角度：theta 和 phi。您还可以使用逐元素数组操作来节省循环：

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;

如果 z 是一个数组（如果 theta 和 phi 来自 meshgrid 而不是 rand），那么你' d 需要 z(:) 来获得数组范围的平均值，否则结果将是一个向量。同样清楚的是，对于 z 组件，您甚至不需要方位角 (phi)，只需要极角 (theta)。

Answer 2

这根本不是如何计算 pi 而不首先使用它。

其实这不是MC的工作方式：它有时需要失败，success/fail的比率会产生一个值。否则，如果你只有成功，那么你必须以某种方式'baked'计算中的搜索值...

这是通过采样单位立方体的体积:

来实现的方法（请耐心等待，我的 Matlab 已经生锈了）

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;

现在如果你真的想使用球体的表面，那么你需要在不知道它是球体的情况下对其进行采样。

例如：

采样单位立方体的三个外面（不经过原点的三个）

normalize() 将点云带到所需的球体表面

然后，对该网格进行三角剖分（哎哟）

计算三角剖分的总面积（再次哎哟） ...并且您有四分之一球体表面的近似值 (pi*r²)，因此可以直接得到 pi！

Monte Carlo 圆周率的近似值

Monte Carlo approximation of pi with a sphere

matlab

random-sample

montecarlo

pi

approximation