Calculating volume of n-ball using Monte-Carlo method

Question

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?

Answer 1

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;

Answer 2

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.