Расчет объема n-шара с использованием метода Монте-Карло

Question

Я пытаюсь создать функцию в MATLAB, которая вычисляет объем n-шара в R^n. Для этого я использую метод Монте-Карло случайной генерации точек в n-кубе, а затем используя отношение точек внутри n-сферы ко всем полученным точкам, умноженным на объем n-куба. Вот код, который я подготовил до сих пор:

function [ approximate_volume ] = MonteCarloHypersphereVolume(radius, dimension, number_of_generations) 
%MonteCarloHypersphereVolume Computes the volume of a 
%'dimension'-dimensional hypersphere by the Monte Carlo method 

number_within_sphere = 0; 
parfor i = 1 : number_of_generations 
    randoms = zeros(1, dimension); 
    for j = 1 : dimension 
     randoms(j) = randi(radius * 2) - radius; 
    end 

    if sum(randoms .^ 2) <= radius^2 
     number_within_sphere = number_within_sphere + 1; 
    end 
end 

approximate_volume = (number_within_sphere/number_of_generations) * (2*radius)^dimension; 

end 

Однако это представляется крайне неточным; согласно Википедии, объем единицы 10-шара должен быть: V_10 = pi^5/5! = 2.5502, однако при запуске функции с 1000000 итерациями она вернулась 11.0067, но ее запуск несколько раз всегда возвращает значение около 11, что намного выше, чем должно быть?

Кроме того, есть ли способ использовать программирование GPGPU для повышения производительности этой функции? Казалось бы, он легко распараллеливается, за исключением зависимости данных от number_within_sphere?

Answer 1

Вы должны использовать rand, а не randi, чтобы пробовать каждое измерение с непрерывным равномерным распределением. То есть, заменить линии

randoms(j) = randi(radius * 2) - radius; 

от

randoms(j) = rand*radius*2 - radius; 

Answer 2

этот подход не масштабируется, так как отношение объема п-размерных единичного шара к объему п мерного куба [-1, 1]^п стремится к нулю экспоненциально быстро (и, следовательно, почти каждая случайная точка внутри единичного куба будет находиться за пределами единичного шара, например, при п = 30 объем куба примерно в 5 * 10 ~ 13 раз больше чем объем шара).

Вместо следует использовать здесь алгоритм Монте-Карло полиномиальной сложности для нахождения объема выпуклого тела, как описано, например, в

http://www.cs.berkeley.edu/~sinclair/cs294/n16.pdf

В форме оно записано, однако, один уже принимает формулу для объема n мерный шар (нам нужно знать объем B_0 в тексте). Вместо последовательности увеличивающихся концентрических шаров, как и в тексте, можно взять последовательность увеличивающихся кубов со схожими свойствами (первый куб - тот, который вписан в единичный шар, последний - [-1,1]^n, а отношение между сторонами последовательных кубов (самое большее) 1 + 1/n), выпуклое тело K является единичным шаром, а затем тот же алгоритм может быть использован для нахождения объема единичного шара.