# How to generate uniformly random points on n-spheres and in n-balls

For many Monte Carlo methods, such as those in graphical computing, it is critical to uniformly sample from $d$-dimensional spheres and balls. This post describes over twenty different methods to uniformly random sample from the (surface of) a $d$-dimensional sphere or the (interior of) a $d$-dimensional ball.

First published: 10th March, 2019 Last updated:    16 April, 2020

# Introduction.

This post gives a comprehensive list of the twenty most frequent and useful methods to uniformly sample from a the surface of a  $d$-sphere, and the interior of the $d$-ball.

Many of them are simple and intuitive; others utilize well-known math methods in sophisticated ways; whilst others can only be described as truly remarkable in their structure. The choice of which one suits your needs will most likely be a balance between code complexity and time efficiency. However, other aspects such as variance reduction may also be of considerable importance.

It is of course, always possible to use the most generalised versions of these formulae as outlined in methods 19-22 for all $d$, but I specifically list some of the specialised methods for smaller $d$ because they may be more efficient, relevant to your application,  and/or amenable to understanding than their  more general expressions. Furthermore, many of the methods applicable to the lower dimensions do not have simple or natural generalization.

### Definition of the $d-$sphere and $d$-ball.

Let us recap the standard definitions of spheres and balls in $d$-dimensions.

A unit $d-$dimensional sphere is defined such that:

$$S^d = \left\{ x \in \Bbb{R}^{d+1} : \left \lvert x \right \rvert = 1 \right\}$$

And a unit $d-$dimensional ball is defined such that:

$$B^d = \left\{ x \in \Bbb{R}^{d} : \left \lvert x \right \rvert \leq 1 \right\}$$

So the perimeter of a circle is a 1-sphere, and the interior of the circle (a disk) is a 2-ball.

And if the earth was perfectly spherical, than the surface of earth would be a 2-sphere, and the interior of  the earth would be a 3-ball.

The way I remember these definitions, is that for both the $d$-sphere and the $d$-ball, the $d$ signifies how many degrees of freedom it has.

### Uniform Random Sampling

In this post, we adopt the convention that ‘uniform random sampling‘, ‘(continuous) simple random sampling‘, and even ‘uniform sampling‘ are all fully synonymous. And that for discrete probabilities, this means that all possible elements of $S$ have an equal probability of being selected. For continuous probabilities, this means that the likelihood of an element falling in any subinterval is directly proportional to the length of the subinterval.

Virtually all computer languages have a built-in function for this concept. For example , $\text{rand}(\cdot)$, that produces a random number that is (pseudo)-randomly drawn from uniform distribution $[0,1)$.

### Note 1.

In contrast to many of my other posts on this blog, we are not talking about quasirandom sampling (aka ‘evenly distributed’), we are just talking about ‘plain old simple vanilla’ uniform random sampling.

### Note 2.

For the coding examples, I generally use readable python code (built on numpy), except for acceptance/rejection methods where I used pseudo-code. The purpose of the code is aimed for clarity and so does not necessarily include language-specific conventions or optimizations. It uses the notation that $x**y$ is equal to $x^y$.

Although each of them have their strengths and weaknesses, an asterisk (*) indicates methods that I believe are an excellent default option for that section.

# Uniformly sampling a 1-sphere (a circle).

That is, uniformly sampling points on the circumferende of a circle. ### *Method 1. Polar

This should be very clear to everyone, and is the foundation for many of the other methods.

theta = 2*PI * random()
x = cos(theta)
y = sin(theta)


### Method 2. Rejection Method

A lesser-known method (Cook, 1957) that does not directly require trigonometric functions.

You might notice that the expressions have almost identical structure to Pythagorean triplets!

(This method is better known when generalised to the 2-sphere.)

Select u,v ~ U(-1,1)
d2 = u^2+v^2
If d2 >1 then
reject and go back to step 1
Else
x = (u^2-v^2)/d2
y = (2*u*v)/d2
Return (x,y)


### Method 3. Muller , Marsaglia (‘Normalised Gaussians’)

This method was first proposed by Muller and then Marsaglia and the popularised by Knuth, it is very elegant. It naturally generalizes and remains efficient for higher dimensions.

Presumably, it is less well-known due to the less obvious (magical?!) relationship between a circle and the Normal Distribution.

Note that in this version (and its generalizations), the standard deviation does not have to be 1, it just has to be identical for all the random variates. However, in nearly all implementations, a standard error of 1 is picked due to its convenience.

Also, in the absence of a native function, a simple and common method to draw from a normal distribution is via the #

# Uniformly sampling a 2-ball (disk).

That is, uniformly sampling points inside a circle.

### *Method 4. Rejection Method

This is probably the most intuitive method, and is quite fast for a 2-ball.

However, (as discussed later) it generally gets a bad rap, because if you generalize this algorithm to higher dimensions, it can become catastrophically very inefficient. The rejection method for $B^2$ has a ~781. Select x,y ~ U(-1,1) 2. If (x^2+y^2 >1) then reject and go to step 1. 3. Return (x,y) Method 5. Polar + Radial CDFFor balls of any dimension, a very common method of uniformly picking from a $d$-ball, is to first select a random directional unit-vector from a $(d-1)$-sphere, and then multiply this vector by a scalar radial multiplicative factor.The crucial aspect of this algorithm is the multiplicative factor, $\sqrt{\cdot}$ function. This is intuitively required because area, which corresponds to the cumulative distribution function (CDF), grows according to $r^2$, so we need to apply the inverse of this. See here for a more visual explanation of this.This method is notably asymmetrical. Geometrically, one variable is used for radial distance and the other for the angular displacement.Figure 5. Polar mapping. Horizontal lines map to concentric circles.u = random() v = random() r = u**0.5 # sqrt function theta = 2* PI *v x = r*cos(theta) y = r*sin(theta) *Method 6. Concentric Map (Shirley 1997)This variance reduction method  (Shirley, 1997) is fascinating and motivated by the idea that ideal mappings from the square-based random variate $uv$-space to the coordinate $xy$-space (disks) should be area-preserving, two-way continuous, and with low distortion. By comparing the figures for method 5 and method 6, it should be clear why this one generally produces better results in (finite) MC methods. Additional visual explanation and analysis of this method is given here (which also include some other variations including ‘squircle’, elliptical and approximate equal areas). Figure 6. Concentric Mapping maps concentric squares to concentric circles.u= random() v = random() if (u==0 && v==0) return (0,0) theta=0 r = 1 a=2*u-1 b=2*v-1 if(a*a>b*b) {r = a; phi = PI/4*b/a} else {r = b; PI/2- PI/4*a/b} x= r*cos(theta) y= r*sin(theta)  Method 7. Muller / Marsaglia (‘Normalized Gaussians’)This method is a hybrid of methods 3 and 5, because it selects a random direction via method 3 and a random radius via method 5.u = np.random.normal(0,1) v = np.random.normal(0,1) norm = (u*u+v*v)**0.5 r = pow( np.random(0,1),1/2.) (x,y) = r*(u,v)/norm Method 8. Exponential DistributionAnother stunningly elegant method (Barthe, 2005) that is valid for balls of any dimension, is to draw $d$ random variates from the normal distribution, and a single random variable from the exponential distribution (with $\lambda=1/2$).(Note that I believe there is a typo in the original paper, as the original paper says $\lambda = 1$, which has then also been copied into Wolfram’s mathworld. )Recall, that in the absence of a native function, a simple and common method to draw from an exponential distribution with parameter $\lambda$, is via the negative logarithm of a uniform random variate divided by $\lambda$.u = np.random.normal(0,1) v = np.random.normal(0,1) e = np.random.exponential(0.5) denom = (u*u + v*v + e)**0.5 (x,y) = (u,v)/denom  Method 9. Dropped CoordinatesA somewhat magical method, first noted by (Harman and Lacko 2010)  and then proven ( Voelker 2017) that is not well-known is that if $(x_1,x_2,x_3,x_4)$ is a random vector uniformly distributed on the $3$-sphere, the random vector $(x_1,x_2)$, (that is, with the last two coordinates dropped) is uniformly distributed in the $2$-ball.See method 20 for a more detailed commentary of this method.s = np.random.normal(0,1) t = np.random.normal(0,1) u = np.random.normal(0,1) v = np.random.normal(0,1) norm = (s*s + t*t + u*u + v*v)**(0.5) (x,y) = (u,v)/norm Results for n=3000, for methods 8 (radial CDF), 9 (exponential) and 10 (dropped coordinates). It is truly remarkable that such different techniques can produce the same result!  Uniformly sampling a 2-sphere (surface of a sphere).That is, uniformly sampling points on the surface of a sphere. Method 10. PolarThis is the natural generalization of method 1. Note that the inverse cosine is required for the same reason that the sqrt is required in method 5.u = random() v = random() theta = 2 * PI * u phi = np.arccos(2v-1) x = np.sin(theta) np.cos(phi) y = np.sin(theta) np.sin(phi) z = np.cos (theta) Using $u=\cos \theta$, an equivalent version of this that does not directly involve inverse trig functions isu = 2*random()-1 phi = 2* PI * random() x = cos(phi) * (1-z^2)**0.5 y = sin(phi)* (1-z^2)**0.5 z = u*Method 11. MullerThis is a direct generalization of method 3.u = np.random.normal(0,1) v = np.random.normal(0,1) w = np.random.normal(0,1) norm = (u*u + v*v + w*w)**(0.5) (x,y,z) = (u,v,w)/norm  Method 12. Rejection Method #1This method (Marsaglia, 1972) is a direct generalization of Method 2 .Select u,v ~ U(-1,1) d2 = u*u + v*v If d2>1 then reject and go to step 1. Else     x = 2*u* sqrt (1-d2) y = 2*v* sqrt (1-d2) z = 1-2*d2 Return (x,y,z) Method 13. Rejection Method #2This method (Cook, 1957) is quite different to any of the previously discussed methods.For those variables that are accepted, the rules for quaternion transformations, imply that the following points are uniformly distributed.Select s,t,u,v ~ U(-1,1) d2 = s^2 + t^2 + u^2 + v^2 If d2>1 then reject and go to step 1. Else     x = 2*(s*u+t*v)/d2 y = 2*(u*v-s*t)/d2 z = (s*s-t*t+u*u-v*v)/d2 Return (x,y,z)  Uniformly sampling a 3-ball (interior of a sphere).That is, uniformly sampling points on the inside of a sphere. *Method 14. Rejection MethodThis method is a direct generalization of Method 2. Again, it is quite efficient as the acceptance rate is still relatively high and the calculations can be done blazingly fast. However, it carries a bad reputation as it is catastrophically inefficient for higher dimensions.1. Select x,y,z ~ U(-1,1) 2. If (x^2+y^2 +z^2>1) then reject and go to step 1. 3. Return (x,y,z)  The rejection method for $B^3$ has a ~52*Method 15. MullerThis is combination of method 7.u = np.random.normal(0,1) v = np.random.normal(0,1) w = np.random.normal(0,1) r = random()**(1./3) norm= (u*u+v*v+w*w)**(0.5) (x,y,z) = r*(u,v,w)/norm  *Method 16. PolarThis is the natural generalization of method 5.u = 2*random()-1 phi = 2* PI * random() r = random()**(1/3.) x = r* cos(phi) * (1-u^2)**0.5 y = r* sin(phi)* (1-u^2)**0.5 z = r*u Method 17. Exponential DistributionThis is a direct generalisation of method 8. Again, note that I believe there is a typo in the original paper, as the original paper says $\lambda = 1$, which has then also been copied into Wolfram’s mathworld.u = np.random.normal(0,1) v = np.random.normal(0,1) w = np.random.normal(0,1) e = np.random.exponential(0.5) denom = (e + u*u + v*v + w*w)**0.5 (x,y,z) = (u,v,w)/denom  Method 18. Dropped CoordinatesA natural extension of method X.Based on the fact that if $(x_1,x_2,x_3,x_4,x_5)$ is a random vector uniformly distributed on the $4$-sphere, the random vector $(x_1,x_2,x_3)$, (that is, with the last two coordinates dropped) is uniformly distributed in the $3$-ball.See method 20 for a more detailed commentary of this method.s = np.random.normal(0,1) t = np.random.normal(0,1) u = np.random.normal(0,1) v = np.random.normal(0,1) w = np.random.normal(0,1) norm = (s*s + t*t + u*u + v*v + w*w)**(0.5) (x,y,z) = (u,v,w)/norm  Uniformly sampling a $d$-sphere. *Method 19. MullerFor higher dimensions where $d>3$, most of the above methods do not neatly simplify. Therefore, the one that is almost always used for a $d$-sphere is the generalised version of methods 3 and 9.u = np.random.normal(0,1,d)  d=np.sum(u**2) **(0.5) (x,y,z) = (u,v,w)/d Uniformly sampling a $d$-ball. Although the rejection methods 4 and 13 are extremely intuitive and trivially generalise to higher dimensions, one quickly encounters ‘the curse of dimensionality’. As mentioned, for $d=1$ the acceptance rate is ~71By comparing the hypervolume of a $d$-ball with the enclosing cube, one can deduce that for $d=10$, the acceptance rate is less than 0.25*Method 20. MullerIn practice, the almost exclusive method to calculate points inside a $d-$ball is the generalized version of method 14. Namely,u = np.random.normal(0,1,d)  # an array of d normally distributed random variables norm=np.sum(u**2) **(0.5) r = random()**(1.0/d) x= r*u/norm  Method 21. Exponential DistributionThis is a direct generalisation of method 8 and 16.x = np.random.normal(0,1,d) e = np.random.exponential(0.5) denom = (e + np.sum(u**2))**0.5 (x,y,z) = (u,v,w)/denom  *Method 22. Dropped CoordinatesA very elegant method, first noted by (Harman and Lacko 2010)  and then proven ( Voelker 2017) that is not well-known is that if $(x_1,x_2,x_3,…,x_d,x_{d+1},x_{d+2})$ is a random vector uniformly distributed on the $(d+1)$-sphere, the random vector $(x_1,x_2,x_3,…,x_d)$, (that is, with the last two coordinates dropped) is uniformly distributed in the $d$-ball.This is non-trivial result where even developing an intuition is difficult, (even the base case of $d=1$ is not trivial!)However, it may be interesting/useful to note that an $\exp(1/2)$ variable $Y$ in method 19, has the same distribution as $(Y_1^2+Y_2^2)$  where $Y_1,Y_2$ are standard normal – thus connecting this methods 19 and 20.Thus,u = np.random.normal(0,1,d+2)  # an array of (d+2) normally distributed random variables norm=np.sum(u**2) **(0.5) u = u/norm x = u[0:d] #take the first d coordinates  ConclusionThis post gives a comprehensive list of the most frequent and useful methods to uniformly sample from a the surface of a  $d$-sphere, an the interior of the $d$-ball.Many of them are simple and intuitive, whilst others are truly remarkable in their structure.It is of course, possible to always use the most generalised versions of these formulae as outlined in methods 19-22, but sometimes specialised methods for smaller $d$ may be more efficient and/or amenable to understanding. ***My name is Dr Martin Roberts, and I’m a freelance Principal Data Science consultant, who loves working at the intersection of maths and computing. “I transform and modernize organizations through innovative data strategies solutions.”You can contact me through any of these channels.LinkedIn: https://www.linkedin.com/in/martinroberts/Twitter: @TechSparx  https://twitter.com/TechSparxemail: Martin (at) RobertsAnalytics (dot) comMore details about me can be found here.Other Posts you may likeA new algorithm for multi-armed banditsHow to evenly distribute points on a sphere more effectively than the canonical Fibonacci lattice

## 6 Replies to “How to generate uniformly random points on n-spheres and in n-balls”

1. Dimitris Kechrakos says:

Excellent summary and presentation of algorithms. Thank you for this contribution.

1. Martin Roberts says:

my pleasure.