You have to use rand
, not randi
, to sample each dimension with a continuous uniform distribution. That is, replace line
randoms(j) = randi(radius * 2) - radius;
by
randoms(j) = rand*radius*2 - radius;
سؤال
I am trying to create a function in MATLAB which computes the volume of an n-ball in R^n. To do this I am using the Monte-Carlo method of randomly generating points in an n-cube and then using the ratio of points inside the n-sphere to all points generated multiplied by the volume of the n-cube. Here is the code I have produced so far:
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
However, this appears to be extremely inaccurate; according to Wikipedia, the volume of a unit 10-ball should be: V_10 = pi^5/5! = 2.5502, however when running the function with 1000000 iterations it returned 11.0067, indeed running it multiple times always returns a value around 11, which is much higher than it should be?
Moreover, is there a way in which to use GPGPU programming to improve the performance of this function? It would appear to be easily parallelizable except for the data-dependency on number_within_sphere
?
المحلول
You have to use rand
, not randi
, to sample each dimension with a continuous uniform distribution. That is, replace line
randoms(j) = randi(radius * 2) - radius;
by
randoms(j) = rand*radius*2 - radius;
نصائح أخرى
this approach does not scale, since the ratio of the volume of n-dimesional unit ball to the volume of the n dimensional cube [-1,1]^n tends to zero exponentially fast (and thus almost every random point inside the unit cube will be outside of the unit ball; for example, for n=30, the volume of the cube is about 5*10^13 times bigger than the volume of the ball).
Instead one should use here the polynomial complexity Monte Carlo algorithm for finding the volume of a convex body as described, for example, in
http://www.cs.berkeley.edu/~sinclair/cs294/n16.pdf
In the form it is written down, however, one already assumes the formula for the volume of the n dimensional ball (we need to know the volume of B_0 in the text). One may take, however, instead of the sequence of increasing concentric balls as in the text, a sequence of increasing cubes with similar properties (the first cube is the one inscribed in the unit ball, the last is [-1,1]^n, and the ratio between the sides of consecutive cubes is (at most) 1+1/n), the convex body K being the unit ball, and then the same algorithm can be used to find the volume of the unit ball.