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
First, here are two methods that do not give uniform distributions:
Random angles 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


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.


Random in cube, then rerandom if out
in = randN(1, 1) // random vector in cube
while (length(in) > 1) { // if outside sphere, rerandom
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%
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
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 x^{2}. The inverse of that integral
function is x^{1/2} so we would
use
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
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 x^{2} x^{1/2} sqrt(rand(0,1)) square 3x^{2} x^{3} x^{1/3} pow(rand(0,1), 1/3) cube 4x^{3} x^{4} x^{1/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.
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
To help understand why this works, we can show samples at regular
intervals between min and max, in place of


Random in 3sphere
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
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 3x^{2} 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


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.


The same is true for the volumes of the elements within the sphere after adjusting the radius with the cube root function.
Random in 4sphere
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(1u) // 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
A 3D rotation can be represented by a unit quaternion, or 4vector, 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.
Random in Nsphere
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
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 BoxMuller transform:
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 Nvector as above in steps 2 and 3.
This polarcoordinate 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
BoxMuller transform
Rayleigh distribution
Back to Karl Sims homepage
© 2024, Karl Sims, All rights reserved.