How to compute volume of 10-Dimentional sphere with Monte-Carlo-Method in Python?

Question

I am trying to compute the volume of a 10 dimensional sphere with python, but my computation doesn't work.

Here is my code:

def nSphereVolume(dim,iter):
    i = 0
    j = 0
    r = 0
    total0 = 0
    total1 = 0

    while (i < iter):
        total0 = 0;
        while (j < dim):
            r = 2.0*np.random.uniform(0,1,iter) - 1.0

            total0 += r*r
            j += 1

        if (total0 < 1):
            total1 += 1
        i += 1

    return np.pow(2.0,dim)*(total1/iter)

nSphereVolume(10,1000)

Here the error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-253-44104dc9a692> in <module>()
     20     return np.pow(2.0,dim)*(total1/iter)
     21 
---> 22 nSphereVolume(10,1000)

<ipython-input-253-44104dc9a692> in nSphereVolume(dim, iter)
     14             j += 1
     15 
---> 16         if (total0 < 1):
     17             total1 += 1
     18         i += 1

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Could someone who knows this algorithm try to run it and tell me what would be corret to implement, to get the volume of this 10-dim sphere? Thanks!

Answer 1

You have multiple problems in your routine.

The error message you get comes from your line

r = 2.0*np.random.uniform(0,1,iter) - 1.0

The function call np.random.uniform(0,1,iter) does not create a single random number. Instead, like most of numpy's functions, it returns an array--in this case, a vector of the length you state (in this case iter). So r is an array as well, since it uses this array, and total0 is an array for the same reason.

Then later you try to evaluate total0 < 1. But the left is an array and the right is a scalar, so numpy does not like the comparison. I could go into more detail on the meaning of the error message, but that is the basic idea.

The solution is to treat r as a vector--in fact, as the random point in the sphere with side-length 2 that you want. You made a typo and used iter rather than dim as the size of the random vector, but if you make that change you will get the point that you want. You can use numpy to quickly get its length and see if the point lies within the sphere of radius one centered at the origin.

Here is some corrected code. I removed one loop--the one where you tried to build up a vector of the correct size, but we now have numpy build it all at once. I also replaced your variable names with more descriptive names and made some other changes.

import numpy as np

def nSphereVolume(dim, iterations):
    count_in_sphere = 0

    for count_loops in range(iterations):
        point = np.random.uniform(-1.0, 1.0, dim)
        distance = np.linalg.norm(point)
        if distance < 1.0:
            count_in_sphere += 1

    return np.power(2.0, dim) * (count_in_sphere / iterations)

print(nSphereVolume(10, 100000))

Note that a mere 1000 iterations of Monte-Carlo gives very bad precision. You will need more iterations to get a useful answer, so I changed the number of repetitions to 100,000. The routine is now slower but gives more consistent answers of around 2.5. This agrees well with the theoretical answer of

np.pi**(10 // 2) / math.factorial(10 // 2)

which evaluates to 2.550164039877345.

---

(This is a response to a comment, to explain the returned value np.power(2.0, dim) * (count_in_sphere / iterations.)

This routine generates random points, where each coordinate is in the interval [-1.0, 1.0). That means these points are uniformly distributed in the hypercube of dimension dim. The volume of a hypercube is the product of its sides. Each side has length 2.0 so we can calculate the hypercube's volume with 2.0 power dim, or np.power(2.0, dim).

We generated iterations points, with count_in_sphere of them in the hypersphere of radius 1 centered at the origin. So the fraction or proportion of points in the hypercube that were also in the hypersphere is count_in_sphere / iterations. Those points were uniformly distributed through the hypercube, so we can estimate that the fraction of the volume of the hypersphere to the volume of the hypercube is the same as the fraction of the random points in those sets. By high-school proportions we therefore get

volume_of_hypersphere / volume_of_hypercube = points_in_hypersphere / points_in_hypercube

realizing that the equation is merely approximate. Multiplying both sides of that equation by volume_of_hypercube, we get

volume_of_hypersphere = volume_of_hypercube * points_in_hypersphere / points_in_hypercube

Substituting, we get

volume_of_hypersphere = np.power(2.0, dim) * count_in_sphere / iterations

which is the value returned by the function. Again, it is merely approximate, but probability theory tells us that the more random points we generate, the better the approximation.