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

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

Answer 1

[John D'Errico]

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.

[Steve]

Thanks a lot for such a detailed answer, I have a question: Why do I need to multiply by the known 10-cube volume to get an approximation of the 10-sphere volume?

[Walter Roberson]

The number you get out is the fraction of occupancy. To convert that to a volume, you need to multiply by the volume of the full cubeoid.

You are not using a cubeoid that is 1 unit per side: you are using a cubeoid that is two units per side, -1 to +1

[Steve]

Dear Walter, thanks a lot got it!

[John D'Errico]

Think of it like this: the fraction C/N has no units. It is simply a ratio, a unit-less number. To compute the sphere volume, you need to see that the fraction must be applied to the total volume of the hypercube. What fraction of the hypercube does the sphere occupy? So even if we look at it purely from a units perspective, you can see this multiplication makes some sense.

Another way of viewing the Monte Carlo integration is to see it as a variation of rectangle rule, but in many dimensions. Again, we would need to scale things by the cube volume.

[Marc]

Does this mean John is out of retirement?? Or are you simply taking pity on poor Steve? Either way it is good to see your name in this user forum.

[John D'Errico]

@Marc - I can never truly retire because there are constantly things I need/want to maintain on the tools I've posted here. A few neat things are on the horizon too. (Major VPI upgrade! SLM gui tool one day...) So I look in occasionally, and if I see something on which I can add something meaningful, I might chime in.

Answer 2

[Walter Roberson]

To calculate volume, you need to divide the successes, S, by the number of attempts, to get the fill fraction. Then multiply the fill fraction by the area of a hypercube each side of which is the same as the "2" you used in "2 * rand(1,100,N)"

Note, the leading 1 in your rand() is not serving any purpose.

In each iteration of your loop, you are summing the same thing. Where are you indexing by "i" ?

[Steve]

Can you please explain me with more details? What should I do within the for loop?

[Walter Roberson]

Please look again:

In each iteration of your loop, you are summing the same thing. Where are you indexing by "i" ?

That's a hint.

[Walter Roberson]

Another hint:

r = reshape(randperm(21), 3, 7);   %some data
disp(r)
disp(sum(r.^2))     %think about this

[Steve]

I need to sum up vectors of 10 squared elements: x1^2 + ... + x10^2 and check for each vector if the rule x1^2 + ... + x10^2 < 1.

I have generated the matrix RandVar = -1+2*rand(1,10,N) and want to check for each level 1<=i<=N the rule above. How can I do this?

[Steve]

Here is my new version, please have a look what is wrong:

N = 100000;
RandVars = -1 + 2*rand(10,N);
V = zeros(N,1);

for i = 1:N
    if (sum(RandVars(:,i).^2) < 1)
        V(i) = 1;
    end
end

% Mean values using the Monte-Carlo method:
M = (1/N)*sum(V,1);
% Error estimation:
C = V - M*ones(N,1);
Error = sqrt((1/(N-1))*sum(C.^2,2));