MATLAB: Volume of n-sphere? Please help me with the version:

MATLABmonte carlosimulationspherevolume

Hello Experts,
I have to write a function for 10-dimensional sphere using Monte – Carlo method.
What I did is this:
N = 1000000;
g = 1;
S = 0;
RandVars = -1 + 2*rand(1,10,N);
for i = 1:N
    if (sum(RandVars(1,10,:).^2) < 1)
        S = S+1;
    end
end
Please revise my code and assist me with computing the volume using the Monte Carlo.
Many thanks, Steve

Best Answer

My guess is this homework was chosen to show you how inefficient Monte Carlo can be when used in a high number of dimensions. THe curse of dimensionality strikes here.
The known volume of a unit 10-sphere is simple enough to compute. Google is your friend here, but I'll just use a tool I've posted on the FEX.
format long g
spheresegmentvolume([],10)
ans =
                2.55016403987735
Use of Monte Carlo to compute volume is a dartboard approach. Thus throw a series darts at a hypercube of edge length 2 that just contains the unit 10-sphere. Count the number of times the dart falls inside the sphere, and divide by the total number of darts thrown. Multiply by the known 10-cube volume, and you get an approximation of the 10-sphere volume.
The trick is to do this efficiently. Do NOT use a loop.
% the dimension of the problem
dim = 10
% the volume of a 10-cube of edge length 2
cubevol = 2^dim;
% sample size
N = 10e6;
% draw the samples
X = 2*rand(N,dim) - 1;
% Distance from the origin for each sample
% Note that I could have skipped the square root, since this is a unit sphere.
R = sqrt(sum(X.^2,2));
% count the number of hits
C = sum(R <= 1)
C =
                   24925
% compute the estimated volume
V = C/N*cubevol
V =
                   2.55232
As you can see, the computed volume was not terribly inaccurate compared to the known volume of 2.5502… 3 significant digits is as much as we can expect.
To be honest, if I were your instructor, I'd have specified a 20 or even 50 dimensional sphere, to make that hit rate vanishingly small. For the 10-sphere, the fraction of darts that actually hit the sphere is only about 0.25%. But that is a large enough fraction that we got almost 3 digits of accuracy from our pretty large sample size.
hitrate = spheresegmentvolume([],10)/cubevol
hitrate =
       0.00249039457019272
A nice thing about Monte Carlo is is allows you to compute a measure of your uncertainty in the final estimate.
A binomial random variable with p = 0.00249039… has variance n*p*(1-p). So if we wish to put +/- 2*sigma confidence limits around the number of hits we should have expected:
C + [-2 2]*sqrt(N*hitrate*(1-hitrate))
ans =
          24609.7735731207          25240.2264268793
And of course, the final volume will have confidence bounds of:
V + cubevol/N*[-2 2]*sqrt(N*hitrate*(1-hitrate))
ans =
          2.52004081388756          2.58459918611244
Again, learn how to avoid loops when you can do so.
Related Question