Generating solid spheres

Question

I'd like to generate random uniform samples from n-dimensional solid spheres.

My current method is this

def sample_sphere(d, npoints):
    points = np.zeros((npoints, d))
    for i in range(npoints):
        r = np.random.rand() #random radius
        v = np.random.uniform(low= -1, high=1, size=d) #random direction
        v = v / np.linalg.norm(v)
        points[i] = r * v
    return points

But unsurpringly, this method gives a higher concentration of points near 0, since the sampling density is not proportional to the volume:

How can I sample uniformly?

Answer 1

There are two basic interpretations of a spherical distribution:

1. Bounded solid sphere

The probability of getting a point at a given radius is given by the volume of a shell with thickness dr at that radius: p(r) ~ r^D up to a constant. The radial CDF is therefore (r / R)^(D+1), where R is the outer radius and D is the dimensionality. The sphere has to be bounded unless you want to select from all of space.

A standard way to generate samples given the CDF is to invert it:

random_radius = R * np.pow(np.random.uniform(0, 1), 1 / D)

By the way, a common way of picking properly uniform directions is to use the Gaussian distribution (see here):

random_direction = np.random.normal(size=D)
random_direction /= np.linalg.norm(random_direction)

You can, and should, vectorize your function. There is no need for looping:

def sample_sphere(r, d, n):
    radii = r * np.pow(1 - np.random.uniform(0, 1, size=(n, 1)), 1 / d)
    directions = np.random.normal(size=(n, d))
    directions *= radii / np.linalg.norm(directions, axis=1, keepdims=True)
    return directions

Notice the 1 - np.random.uniform(...). This is because the range of the function is [0, 1), while we want (0, 1]. The subtraction inverts the range without affecting the uniformity.

2. Unbounded decaying distribution

In this case you're looking for something like a Gaussian with spherical symmetry. The radial CDF of such a distribution scales with R^(D-1). For the Gaussian, this is called the Rayleigh distribution in 2D, Maxwell in 3D, and chi in the general case. These are available directly through the scipy.stats.rayleigh, scipy.stats.maxwell and scipy.stats.chi objects, respectively.

In this case, radius has less meaning, but it can be interpreted as the standards deviation of your distribution.

The unbounded version of your function (also vectorized), then becomes:

def sample_sphere(s, d, n):
    radii = scipy.stats.chi.rvs(df=d, scale=s, size=(n, 1))
    directions = np.random.normal(size=(n, d))
    directions *= radii / np.linalg.norm(directions, axis=1, keepdims=True)
    return directions