Generate random points inside a sphere

For reasons I don’t understand, it’s more common, at least in my experience, to see references on how to generate points on a sphere than in a sphere. But the latter is easily obtained from the former.

The customary way to generate points on a sphere in n dimensions is to generate points from an n-dimensional normal (Gaussian) distribution and normalize by dividing by the length. And the way to generate samples from an n-dimensional normal is simply to generate a list of n one-dimensional normals [1].

The way to generate points in a sphere is to generate a point S on the sphere then randomly adjust its radius r. If you generate a uniform random point in [0, 1] and set r = u then you’ll over-sample points near the origin. Instead, the thing to do is to set

r = u^1/n

That’s because volume is proportional to rⁿ. Taking the nth root of u mades the volume inside a sphere of radius r uniformly distributed.

Demonstration

To illustrate this, I’ll use the same technique as the post on sampling from a tetrahedron: generate random points in a sphere and see whether the right proportion fall inside a (hyper)cube inside the sphere.

We’ll work in n = 4 dimensions and use a sphere of radius R = 5. Our test cube will be the cube with all vertex coordinates equal to 1 or 2. We’ll have to use a lot of sample points. (This is why naive Monte Carlo breaks down in high dimensions: it may be that few if any points land in the region of interest.)

Here’s the code to do the sampling.

import numpy as np
from scipy.stats import norm

def incube(v):
    retval = True
    for x in v:
        retval = retval and 1 <= x <= 2
    return retval

def sample(R, n):
    x = norm.rvs(size=n)
    x /= np.linalg.norm(x)
    u = np.random.random()
    x *= R*u**(1/n)
    return x

R = 5
n = 4
N = 100_000_000
V = R**n * np.pi**(n/2) / gamma(n/2 + 1)
print(V)

count = 0

for _ in range(N):
    v = sample(R, n)
    if incube(v):
        count += 1        

print(count/N, 1/V)

This prints

0.0003252 0.000324

The two results agree to four decimal places, which is about what we’d expect given 1/√N = 0.0001.

[1] If you’re porting the sampler above to an environment where you don’t have a function like norm.rvs for generating normal random samples, you can roll your own using the Box-Muller transformation. For example, here is a C++ implementation from scratch.