Приближение Монте-Карло пи сферой

Я пытаюсь оценить пи путем равномерного случайного отбора точек (x, y) внутри круга радиуса 1, а затем вычислить соответствующее значение z в сфере. На самом деле это всего лишь четверть круга для упрощения вычислений. Затем я вычисляю среднее значение z (которое должно составлять примерно 1/8 объема всей сферы), а затем сравниваю его с объемом куба 1x1x1. Оно должно быть примерно 1/6 пи, но по какой-то причине это не так.

Это мой код 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

Он все время повторяет, что пи примерно равно 4, что бессмысленно, даже если я увеличу количество отсчетов. Не могли бы вы объяснить мне, в чем проблема, несмотря на то, что я использую число Пи для определения угла при выборке случайных точек внутри круга?


matlab approximation montecarlo pi random-sample
person Rafael K.    schedule 22.11.2015    source источник
comment
Среднее значение будет 1/6 pi, если вы позволите x и y изменяться по всей единице квадрат [0, 1] x [0, 1]. Как бы то ни было, можно ожидать, что среднее значение будет (1/6 pi) / (1/4 pi) = 2/3.   -  person Mark Dickinson    schedule 23.11.2015
comment
хорошо, спасибо, понятно. Но я бы хотел избежать выборки отбраковки на квадрате, которая в более высоких измерениях была бы еще более неэффективной. Вы видите способ исправить это?   -  person Rafael K.    schedule 23.11.2015
comment
Я действительно не понимаю, чего вы пытаетесь достичь. Почему вы не можете сделать это, используя круг вместо сферы? (Например, выберите x равномерно на [0, 1] и посмотрите на среднее значение sqrt(1 - x*x).)   -  person Mark Dickinson    schedule 23.11.2015
comment
Да, я попробовал это с кругом, и это сработало, но я хотел обобщить его на более высокие измерения.   -  person Rafael K.    schedule 23.11.2015
comment
1. mean(z)(1) не валидный матлаб. Вы используете октаву? 2. Я предлагаю избегать магических чисел, устанавливать N=1000 в начале и использовать N каждый раз позже 3. Используйте операции с массивами вместо цикла: например, z=(1-x.^2-y.^2).^(0.5) (более эффективно, результат тот же). 4. В сфере вам нужны тэта и фи и x=r.*sin(theta).*cos(phi); y=r.*sin(theta).*sin(phi); z=r.*cos(theta). В вашей текущей версии x.^2+y.^2==r для всех точек, поэтому вы просто усредняете sqrt(1-r), что не имеет смысла.   -  person Andras Deak    schedule 23.11.2015
comment
Кроме того, вы используете pi в функции, которая вычисляет pi, что не очень помогает!   -  person David    schedule 23.11.2015
comment
Комментарий @David очень хорошо показывает, насколько ваш подход изначально ошибочен. Как предложил Марк Дикинсон, вы должны создать точки на кубе и подсчитать долю точек внутри единичной сферы. Затем используйте формулу n-мерной сферы для аппроксимации числа пи.   -  person Andras Deak    schedule 23.11.2015
comment
@AndrasDeak Да, я знаю, что использовать пи при аппроксимации пи само по себе не очень полезно, но это обусловлено тем фактом, что sin и cos получают параметры в радианах, если бы они получили их в градусах, это не было бы проблема.   -  person Rafael K.    schedule 23.11.2015
comment
Затем используйте dec2rad и скройте пи;) И попробуйте выделенную жирным шрифтом часть моего предыдущего комментария.   -  person Andras Deak    schedule 23.11.2015
comment
@AndrasDeak может использовать формулы аппроксимации, такие как ряд Тейлора, или есть другие методы, такие как формула Бхаскары.   -  person Rafael K.    schedule 23.11.2015
comment
@AndrasDeak да, спасибо, я пытаюсь понять, как это работает, потому что это продолжает давать мне странные результаты, когда я пробую выделенную жирным шрифтом часть вашего предыдущего комментария.   -  person Rafael K.    schedule 23.11.2015
comment
Я добавил ответ, разъясняющий мою точку зрения выше. Однако я очень не уверен в вашей формуле из mean. Например, вы сравниваете объем с отдельным компонентом в своем вопросе. Итак, выполнили ли вы сферический интеграл, необходимый для вычисления <z> на сфере? Вы уверены, что он в 1/6pi раз больше куба?   -  person Andras Deak    schedule 23.11.2015
comment
@AndrasDeak Спасибо за ответ. Наверное, я не понимаю, о чем вы. На самом деле, я не очень хорошо знаком со сферическими интегралами и вообще интегралами в более чем двух измерениях. Но я ожидал, что это должно быть что-то вроде средней координаты z, если я выбираю равномерно случайным образом точки из круга в плоскости x, y. И да, это можно сделать с помощью куба и отбраковки, но я хотел попробовать другой подход.   -  person Rafael K.    schedule 23.11.2015
comment
Что ж, я обнаружил новую проблему. Ваши точки неоднородны на единичной сфере (еще одна причина, почему вы должны использовать выборку отклонения на кубе). И даже если бы они были, вы не можете не делать интегралы сначала вручную. Процедура будет такая: вычислить ‹z› на бумаге, приблизительное ‹z› из Монте-Карло, а затем сравнить их.   -  person Andras Deak    schedule 23.11.2015
comment
Рафаэль, у меня все время возникают проблемы. На этой странице обсуждается, как равномерно выбирать точки на сферической поверхности, но тогда вам нужно pi, чтобы сказать что-нибудь о результате. Я действительно не понимаю прямо сейчас, как мы могли бы спасти этот подход ...   -  person Andras Deak    schedule 23.11.2015
comment
@Andras Deak: да, я тоже вижу слишком много проблем, но спасибо за ваши усилия :) Также я наконец понял, что в моем первоначальном подходе, если я хочу получить приближение объема сферы, мне нужно рассматривать как основание четверть круга, по которому я пытался вычислить среднюю высоту, а не квадрат, и, конечно же, для этого мне нужно число пи :(   -  person Rafael K.    schedule 23.11.2015
comment
Да, именно это я понял, когда попытался интегрировать по сфере ... Очень традиционный способ аппроксимации числа Пи из Монте-Карло (фактически, один из первых методов Монте-Каро) - это Проблема с иглой Буффона.   -  person Andras Deak    schedule 23.11.2015

Ответы (2)


arrow_upward
1
arrow_downward

Помимо того факта, что вам нужно 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), тогда вам нужно было бы z(:), чтобы получить среднее значение по всему массиву, иначе результат был бы вектором. Также ясно, что для компонента z вам даже не нужен азимутальный угол (phi), только полярный угол (theta).

person Andras Deak    schedule 22.11.2015
comment
При использовании deg2rad по-прежнему используется pi! - person David; 23.11.2015
comment
@ Дэвид: да, да, но это не по своей сути недостаток. Сначала я тоже так думал, но ряд Тейлора sin / cos не включает пи, поэтому ничто не мешает вам вычислить их в градусах. Просто так получилось, что эти функции определены с входными данными в радианах. - person Andras Deak; 23.11.2015
comment
@David Я удалил радианы с изображения, см. Мое редактирование: P Хотя я все еще немного не уверен. Тот факт, что ряд Тейлора включает x ^ k / k! вероятно связано с радианами. Но вы, вероятно, могли бы сформулировать это без числа Пи, просто связав его с количеством arbitrary degree в полном круге ... - person Andras Deak; 23.11.2015
comment
: P: P: P: P Я вижу, что ты там делал !! - person David; 23.11.2015

arrow_upward
1
arrow_downward

Это совсем не то, как вычислить Пи, не используя его в первую очередь.

На самом деле MC работает не так: иногда он должен терпеть неудачу, и соотношение успех / неудача дает значение. в противном случае, если у вас есть только успехи, вы должны каким-то образом `` запечь '' искомое значение в вычислении ...

Вот как это сделать (терпите, мой 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!
person MrBrushy    schedule 11.10.2016