How to generate random points in N-dimensional spheres

Karl Sims


Here is a collection of methods for generating random vectors on, and in, spheres of various dimensions.

We want the results uniformly distributed on the surface of a sphere, or within the volume of a sphere. Each random point will be independent. All example spheres will have radius 1 and center at 0.

Assume we have a rand(min, max) function that can generate uniformly distributed random values from min to max, and randN(min, max) functions that give N-dimensional vectors where each component is independently random from min to max. The normalize(v) function generates a unit vector from any input vector.

First, here are two methods that do not give uniform distributions:


Random angles and radius

A random point within a circle can be generated from a random angle and radius: 
a = rand(0, 2*pi)      // angle
on = [cos(a), sin(a)]  // xy on circle
in = on * rand(0, 1)   // xy in circle with random radius
But note the higher density of random points in the center of the circle.

Similarly a random point in a 3D sphere can be generated from random spherical coordinates such as latitude, longitude and radius: 

lat = rand(-pi/2, pi/2)      // latitude
long = rand(-pi, pi)         // longitude
on = [cos(lat) * sin(long),  // xyz on sphere surface
      sin(lat),
      cos(lat) * cos(long)]
in = on * rand(0, 1)   // xyz in sphere with random radius
nonuniform random in circle

in circle, nonuniform

nonuniform random on sphere

on sphere, nonuniform
 nonuniform random in sphere

in sphere, nonuniform

Here the density is higher around the poles and near the center of the sphere. These methods will be improved below by adjusting the random radius and latitude values.


Random in cube, then warp to sphere 

v = randN(-1, 1)      // random vector in cube
on = normalize(v)     // warp to sphere surface
in = on * rand(0, 1)  // put in sphere with random radius
This is a little better, but note the higher density near the corners and edges of the cube, and again near the center of the sphere.

nonuniform random on sphere

on sphere, nonuniform
 nonuniform random in sphere

in sphere, nonuniform



Random in cube, then re-random if out 

in = randN(-1, 1)         // random vector in cube
while (length(in) > 1) {  // if outside sphere, re-random
   in = randN(-1, 1)   
}
on = normalize(in)  // warp to sphere surface, if needed

This method is simple and gives correct results with a uniform distribution in or on the sphere. It may be practical in some cases, but it can have performance issues. The number of attempts varys for each random point and is unbounded, which may be problematic for parallel GPU computation where the slowest affects the total performance. Also, for higher dimensions, the chances of a random point falling within the sphere decreases sharply, which can cause an unreasonable number of iterations.

sphere volume    cube volume    chances in
2D     π 4 79%
3D 4/3 π 8 52%
4D 1/2 π2 16 31%
5D 8/15 π2 32 16%
6D 1/6 π3 64 8%
random in square

in square, uniform


How to generate random numbers with nonuniform distributions

In general, to make random numbers with higher chances of some values than others, we can use this procedure:

1. Define a desired probability distribution curve, or density function.
2. Find the integral function of that desired distribution.
3. Find the inverse function of that integral.
4. Apply that inverse integral function to the rand(min, max) values, with min and max matching the integral function's range.

This will transform uniform random values into the desired distribution.

For example, let's say we want random numbers from 0 to 1 but with a linear distribution that gives more higher values. The probability distribution curve could be x, but for convenience we can scale that so the area under the curve sums to 1, which instead gives 2x with an integral of x2. The inverse of that integral function is x1/2 so we would use sqrt(rand(0,1)) to generate a linear distribution of random numbers.

Note that we can scale the desired distribution curve and/or offset the integral by any constants without affecting the results, as long as the rand range matches the adjusted output range of the integral function (which is the same as the input domain of the inverse integral function). If the integral function's range is from 0 to 1, it may be called a "cumulative distribution function".

This table shows various probability distributions and the corresponding inverse integral functions that can generate them.

desired distribution integral    inverse of integral    applied to rand
uniform 1 x x rand(0,1)
linear 2x x2 x1/2 sqrt(rand(0,1))
square 3x2 x3 x1/3 pow(rand(0,1), 1/3)
cube 4x3 x4 x1/4 pow(rand(0,1), 1/4)
sin hump      sin(x) -cos(x) acos(-x) acos(rand(-1,1))
cos hump cos(x) sin(x) asin(x) asin(rand(-1,1))
Gaussian 2/√π e-x2      erf(x) erfinv(x) erfinv(rand(-1,1))
Rayleigh x e-x2/2 -e-x2/2 (-2 log(-x))1/2 sqrt(-2*log(rand(0,1)))

Some of these formulas will be used below.


Random in 2-sphere

For a circle we can use a random angle and radius if we adjust the random radius to avoid a higher density in the center of the circle. The relative probability of a given radius should be proportional to the circumference at that radius, which is a simple linear relationship. So as described in the example above, we square root a random number to give a radius with a linear probability distribution, and a final uniform distribution within the circle. 

a = rand(0, 2*pi)          // random angle
on = [cos(a), sin(a)]      // random xy vector on circle
in = on * sqrt(rand(0,1))  // radius is sqrt of random
uniform random points in circle

in circle, uniform

To help understand why this works, we can show samples at regular intervals between min and max, in place of rand(min, max). For uniform results, the "rectangles" between neighboring samples should have equal areas (in the limit with more samples and smaller areas). The densities of the regular samples along single dimensions don't need to equal each other, as long as the final multidimensional densities are equal. These images show how this is achieved in a circle by adjusting the radius with the square root function.

nonuniform areas in circle

nonuniform areas
 uniform areas in circle

uniform areas



Random in 3-sphere

For a 3D sphere we can also use random angles and radius if we adjust the probability distribution of both the latitude and the radius. We can also benefit from the property that uniformly distributed points on a sphere's surface are also uniformly distributed on any axis through the sphere. So we can start by setting y to a uniform random value from -1 to 1. If the y axis is aligned with the poles, the latitude with a correctly adjusted distribution would be asin(y) but we don't need to actually calculate that. The radius on the xz plane at y equals cos(latitude) which is also equal to sqrt(1 - y*y). These values, along with a random longitude, can generate an evenly distributed point on the sphere surface.

There are other spherical coordinate systems that can be used with similar results. Some methods use arccos to calculate a polar angle, or arcsin to calculate the latitude, but avoiding the inverse trig functions should be more efficient.

For a random point within the sphere, we need a final random radius with a probability proportional to the spherical surface area at that radius. The scaled distribution curve for that is 3x2 and the integral inverse of that is the cube root function. 

y = rand(-1, 1)       // random y coordinate
r = sqrt(1 - y*y)     // radius on xz plane at y
long = rand(-pi, pi)  // random longitude
on = [r * sin(long),  // xyz on sphere surface
      y,
      r * cos(long)]
in = on * pow(rand(0,1), 1/3) // radius is cube root of random


uniform random points on sphere

on sphere, uniform
 uniform random points in sphere

in sphere, uniform


Again we can use samples at regular intervals to show why this works. For uniform results on the sphere surface, the latitudes have been adjusted so the rectangles between the samples have equal areas.

nonuniform areas on sphere

nonuniform areas
 uniform areas on sphere

uniform areas

The same is true for the volumes of the elements within the sphere after adjusting the radius with the cube root function.



Random in 4-sphere

For a 4D sphere we can use two random angles to generate points on two circles, interpolate those with square roots of another random value, and finally scale by a random radius using the 4th root function. 
a1 = rand(0, 2*pi)     // xy angle
a2 = rand(0, 2*pi)     // zw angle
u = rand(0, 1)         // interp
r1 = sqrt(u)           // xy radius
r2 = sqrt(1-u)         // zw radius
on = [cos(a1) * r1,    // xyzw vector on 4D sphere
      sin(a1) * r1, 
      cos(a2) * r2, 
      sin(a2) * r2)]
in = on * pow(rand(0,1), 1/4) // radius is 4th root


It can be harder to visualize why this 4D method works, but note that the areas between regularly spaced samples within the xy and zw circles are all equal because their radii are square rooted, as explained for the individual circle above. Therefore the volume elements between regular samples will also have equal volumes. Also note that r12+r22=1 and sin2+cos2=1 so the surface vector will have magnitude 1 as needed. The 4th root of the final radius follows the expected pattern, and causes equal 4D volumes between the regularly spaced samples within the sphere.

A 3D rotation can be represented by a unit quaternion, or 4-vector, so a random point on a 4D sphere can be used to perform a random rotation. The example image shows a small colored line on the sphere that has been rotated by many different random unit quaternions.

uniform random rotations

uniform random rotations


Random in N-sphere

Finally, here is a general case solution:

1. Generate a vector with each component a random number with a Gaussian, or normal, distribution centered on zero. Two options for calculating this are given below. The Gaussian function is separable, so this distribution takes on a spherical blob shape with no orientation.

2. Normalize that vector to give a uniform distribution on the sphere surface.

3. Scale by the Nth root of another uniform random number to get a random point uniformly distributed within the sphere.

The integral of the Gaussian function is the "error function" or erf, and its inverse is simply the "inverse error function" or erfinv, which can transform uniform random numbers into a Gaussian distribution. 

v = erfinv(randN(-1, 1))       // random Gaussian vector
on = normalize(v)              // put on sphere surface
in = on * pow(rand(0,1), 1/N)  // radius is Nth root

An efficient alternative to this without the inverse error function, is to generate pairs of Gaussian distributed values using the Box-Muller transform:

Gaussian distribution
Gaussian distribution
a = rand(0, 2*pi)             // random angle
r = sqrt(-2 * log(rand(0,1))) // random radius with a Rayleigh distribution
x = r * cos(a)                // xy pair, each with a Gaussian distribution
y = r * sin(a)
Repeat that for each of the N/2 pairs of vector values needed. If N is odd, just discard one value from the last pair. Then normalize and scale the whole N-vector as above in steps 2 and 3.

This polar-coordinate transform method is convenient because when 2D coordinates have standard normal distributions, their magnitudes (radii) have a Rayleigh distribution: x e-x2/2 which has a more favorable integral: -e-x2/2. The inverse of that: (-2 log(-x))1/2 can transform uniform random numbers into a Rayleigh distribution as needed for the radii. Then combining the radius with a random angle generates a pair of values, each with an independent Gaussian distribution.



References

Inverse transform sampling
Uniform Random Rotations, Ken Shoemake, 1992
Quaternions and spatial rotation
Error function
Box-Muller transform
Rayleigh distribution



Back to Karl Sims homepage

© 2024, Karl Sims, All rights reserved.