Picking random points in the volume of sphere with uniform probability

I have a sphere of radius $R_{s}$, and I would like to pick random points in its volume with uniform probability. How can I do so while preventing any sort of clustering around poles or the center of the sphere?

Since I'm unable to answer my own question, here's another solution:

Using the strategy suggested by Wolfram MathWorld for picking points on the surface of a sphere: Let $\theta$ be randomly distributed real numbers over the interval $[0,2\pi]$, let $\phi=\arccos(2v−1)$ where $v$ is a random real number over the interval $[0,1]$, and let $r=R_s (\mathrm{rand}(0,1))^\frac13$. Converting from spherical coordinates, a random point in $(x,y,z)$ inside the sphere would therefore be: $((r\cos(\theta)\sin(\phi)),(r\sin(\theta)\sin(\phi)),(r\cos(\phi)))$.

A quick test with a few thousand points in the unit sphere appears to show no clustering. However, I'd appreciate any feedback if someone sees a problem with this approach.

• Summary: Nate proposed a transformation method, while Kevin proposed a rejection method. – J. M. is a poor mathematician Dec 1 '11 at 3:08
• I think the solution you took from MathWorld is sound, if a bit slower than the approaches proposed in the answer (evaluating transcendental functions is expensive!). – J. M. is a poor mathematician Dec 1 '11 at 3:58
• I don't know much but, my approach would be something along the lines of choosing a random pole from the surface to the center and then choosing a point on that line where the probability is greater the closer to the surface, to account for the expansion of the sphere. – Sebastian Garrido Dec 10 '13 at 20:17

Let's say your sphere is centered at the origin $(0,0,0)$.

For the distance $D$ from the origin of your random pointpoint, note that you want $P(D \le r) = \left(\frac{r}{R_s}\right)^3$. Thus if $U$ is uniformly distributed between 0 and 1, taking $D = R_s U^{1/3}$ will do the trick.

For the direction, a useful fact is that if $X_1, X_2, X_3$ are independent normal random variables with mean 0 and variance 1, then $$\frac{1}{\sqrt{X_1^2 + X_2^2 + X_3^2}} (X_1, X_2, X_3)$$ is uniformly distributed on (the surface of) the unit sphere. You can generate normal random variables from uniform ones in various ways; the Box-Muller algorithm is a nice simple approach.

So if you choose $U$ uniformly distributed between 0 and 1, and $X_1, X_2, X_3$ iid standard normal and independent of $U$, then $$\frac{R_s U^{1/3}}{\sqrt{X_1^2 + X_2^2 + X_3^2}} (X_1, X_2, X_3)$$ would produce a uniformly distributed point inside the ball of radius $R_s$.

• This also generalizes to the $n$-sphere (with $D = R_s U^{1/n}$). – Michael Lugo Dec 1 '11 at 4:58
• @nullgeppetto: The pdf is equal to a constant on the ball $B(R_s)$, and 0 outside the ball. That's what "uniformly distributed" means. The constant will be 1 over the volume of the ball (so that the pdf integrates to 1). You can find a formula for this volume (with several derivations) on Wikipedia. Yes, some ways of writing the formula involve the gamma function. – Nate Eldredge Jan 30 '15 at 15:53
• @NateEldredge, thank you for your answer. Please clarify this for me: I have found that the pdf is equal to the inverse of the volume of the $n$-ball of radius $r$, that is, $f(\mathbf{x})=\frac{\Gamma\Big(\frac{n}{2}+1\Big)}{\pi^{\frac{n}{2}}}r^{-n}$, when $\mathbf{x}$ is in the $n$-ball, and zero everywhere else. Is that true? Finally, in terms of computational efficiency (e.g. in a C program), what would be your suggestion about the formula of the pdf (should I prefer the gamma function or not?). Thank you very much! – nullgeppetto Jan 31 '15 at 19:01
• @j13r: I don't understand your objection. My claim is that if we generate independent random variables $(X_1, X_2, X_3, U)$ with $X_i \sim N(0,1)$ and $U \sim U(0,1)$, then the random vector $\frac{R_s U^{1/3}}{\sqrt{X_1^2 + X_2^2 + X_3^2}} (X_1, X_2, X_3)$ is uniformly distributed on the ball. Of course there are many other ways to do this, and you are welcome to post another answer describing alternative approaches. I understand that you don't need to generate 4 random numbers for this vector - but you can. – Nate Eldredge May 10 '15 at 14:34
• @NateEldredge Cant we say that since the volume of the sphere is $\frac{4}{3}\pi {R_s}^3$, pdf can be defined by $P(r,\theta,\phi)= \frac{r^2 \sin(\theta)}{\frac{4}{3}\pi {R_s}^3}$. Integrating $\int\int Pd\phi d\theta = \frac{3r^2}{{R_s}^3}=P(r)$, where $\int_0^r P(x)dx$ would give the same CDF. This can be still inverted, to sample from. I don't really get the full detail of what you wrote. More precisely why $D=R_SU^{\frac{1}{3}}$ – Alexander Cska Apr 7 '16 at 20:00

An alternative method in $3$ dimensions:

Step 1: Take $x, y,$ and $z$ each uniform on $[-r_s, r_s]$.

Step 2: If $x^2+y^2+z^2\leq r_s^2$, stop. If not, throw them away and return to step $1$.

Your success probability each time is given by the volume of the sphere over the volume of the cube, which is about $0.52$. So you'll require slightly more than $2$ samples on average.

If you're in higher dimensions, this is not a very efficient process at all, because in a large number of dimensions a random point from the cube is probably not in the sphere (so you'll have to take many points before you get a success). In that case a modified version of Nate's algorithm would be the way to go.

• Could you please explain why this technique works? – Marc Ourens May 29 '16 at 22:17
• The idea is to pick points uniformly from the cube (which is easy, since you can just pick each coordinate separately), then toss out any point not in the sphere. As long as your points are uniform in the larger set, they'll stay uniform when you restrict to the smaller set. – Kevin P. Costello May 30 '16 at 1:06
• A 3-cube has 8 corners; a 16-cube has 32,768 corners. Won't this method spend most of its time (in 16 dimensions) generating points that are in the corners but not inside the ball? – Reb.Cabin Feb 15 at 14:04
• This relates to what I said in the last paragraph about the process being inefficient in higher dimensions. It turns out the corners themselves aren't the problem -- in high dimensions there are many corners, but very little of the volume of a cube is near a corner. But even if you're away from a corner, you're still very likely to not be in the sphere -- if I pick $n$ random numbers from $[-1,1]$, the expected sum of squares is $\frac{n}{3}$, and I'm throwing out all the points whose square sum is less than $1$. – Kevin P. Costello Feb 15 at 21:05
• An alternative way of seeing just how inefficient this process is in high dimensions: The cube $[-1,1]^n$ has volume $2^n$, but the unit sphere has volume going to $0$ faster than any exponential (see, for example, the argument due to Fedya (listed as Greg Kuperberg's due to an editing quirk) at mathoverflow.net/questions/8258/… – Kevin P. Costello Feb 15 at 23:31

Nate and Kevin already answered the two I knew. Recalling this and this, I think that another way to generate a uniform distribution over the sphere surface would be to generate a uniform distribution over the vertical cylinder enclosing the sphere, and then project horizontally.

That is, generate $z \sim U[-R,R]$, $\theta \sim U[0,2\pi]$, and then $x=\sqrt{R^2-z^2} \cos(\theta)$, $y=\sqrt{R^2-z^2} \sin(\theta)$. This (if I'm not mistaken) gives a uniform distribution over the sphere surface. Then, apply Nate's recipe to get a uniform distribution over the sphere volume.

This method is a little simpler (and more efficient) than the accepted answer, though it's not generalizable to other dimensions.

• In case it isn't immediately clear to people, this method works because of Archimedes' theorem on slicing cylinders and spheres: mathworld.wolfram.com/ArchimedesHat-BoxTheorem.html – Michael Lugo Dec 1 '11 at 5:00
• I believe you need a factor of $\sqrt{R^2 - z^2}$ in the choice of $x$ and $y$, in order to be on the sphere's surface. – Erik P. Feb 21 '13 at 15:57
• @ErikP. Of course, you're right, I forgot to do the projection, I'd left the points over the cylinder surface. Fixed, thanks! – leonbloy Feb 21 '13 at 18:23
• I think you need $z \sim U[-R,R]$? Or a coin flip for the sign. – copper.hat Dec 22 '15 at 21:16
• In pure mathematics, this method works fine. If your goal is to generate a pseudorandom point in a computer program, however, the uniform distribution of points on the cylinder will resemble a square grid of lattice points wrapped around the cylinder. Near where the sphere is tangent to the cylinder, this will map down to a nice rectangular grid of points on the sphere, but farther away it creates two "poles" on the sphere around which you have discrete rings of closely-spaced points. For that reason I'd be inclined to use one of the slightly less efficient methods. – David K Feb 4 '18 at 17:43

I just want to add a small derivation to leonbloy's answer, which uses calculus instead of geometrical intuition.

Changing from cartesian $(x,y,z)$ to spherical $(r,\theta,\phi)$ coordinates, we have for the volume element $$dx dy dz =r^2 \sin \theta ~ dr d\theta d\phi$$ The coordinates $(r,\theta,\phi)$ don't work for a uniform distribution because we still have a non-constant factor in front of $dr d\theta d\phi$ (see "EDIT" at the bottom, if you do not see why they don't work). Therefore we introduce $$u=-\cos \theta \Rightarrow du= \sin \theta d\theta$$ $$\lambda=r^3/R^3 \Rightarrow d \lambda=\frac{3}{R^3}r^2dr$$ with which we obtain an expression with a constant pre-factor $$dx dy dz= \frac{R^3}{3} d\lambda du d\phi$$ The range of our variables is $\lambda \in [0,1], ~u \in [-1,1], \phi \in [0, 2\pi)$. Choosing those numbers uniformly we get cartesian coordinates

\begin{align} x&=r \sin(\theta) \cos (\phi) =&R \lambda^{1/3} \sqrt{1-u^2}\cos(\phi)\\ y&=r \sin(\theta) \sin (\phi) =&R \lambda^{1/3} \sqrt{1-u^2}\sin(\phi) \\ z&=r \cos (\theta)=&R \lambda^{1/3} u \end{align}

EDIT: I want to add an argument why we want a constant prefactor in front of $d\lambda du d\phi$.

Consider the one dimensional case (uniform distribution of points on the line $[0,L]$). For $0<x<L$ the probability to find a point in $(x,x+dx)$ is $P(x)dx$. Since we assume a uniform probability, $P(x)$ has to be $P(x)=1/L$, and hence the probability $P(x)dx=dx/L$ is directly proportional to the volume element $dV=dx$.

Now consider we have a variable $y$, for which we do not know the probability density $Q(y)$ but we know that the volume element is $dV=dx=c dy$ with some constant $c$. Furthermore, we know that $Q(y)dy$ has to be $P(x)dx$ (by definition of probability density). Hence $Q(y)=P(x)dx/dy=c$.

In summary we have shown:

"Variable $y$ is uniformly distributed" $\Leftrightarrow$ "The volume element is $dV=c dy$ for some constant $c$. (For the correct normalization of the probability density the value of $c$ is not arbitrary)

• Why do you have $r^3/R^3$. I don't get this. I also think that $u\in[-1,1]$, since $u$ is a cosine and it goes form $-1$ to $1$. I agree only for the $\pi \in [0,2\pi)$ since if $2\pi$ was included, we would have double counting of the positive $x$ axis. – Alexander Cska Apr 7 '16 at 17:24
• You are correct about $u \in [-1,1]$, thanks for pointing that out! Regarding $\lambda=r^3/R^3$: This is a choice to guarantee that for the new variables $(\lambda,u,\phi)$ the volument element is $dV=const. d\lambda du d\phi$ instead of $dV=f(\lambda,u,\phi) d\lambda du d\phi$ – thomasfermi Apr 9 '16 at 9:32
• Actually by setting $\lambda=r^3$ without the $\frac{1}{R^3}$ you would still get a constant factor. I think that things are more complicated than a simple substitution. Could you please have a look at my comment under the first answer to this question. I proposed some explanation, why you need the $\frac{r^3}{R^3}$ . Regarding your method, of reducing the factor before the differential element to unity, could you give me some reference where I can read more about it. – Alexander Cska Apr 9 '16 at 10:27
• The choice $\lambda'=r^3$ is fine as well. Instead of $\lambda \in [0,1]$, we get $\lambda' \in [0,R^3]$ and $dV=d\lambda' du d\phi/3$ and $R\lambda^{1/3}=\lambda'^{1/3}$. Regarding the reference: I have none, but made this up by myself. I edited my answer to add an explanation at the end. – thomasfermi Apr 9 '16 at 12:38

Much simpler way would be to pair surfaces at different distances from the center.

The surface of a sphere of radius R has area 4 π R^2

So just pair up various surfaces that add up to the same area.

For example:

Surface of sphere with r = 0 paired with surface of sphere with r = R to have a total area of 4π R^2

Surface of sphere with r = R/2 paired with surface of sphere with r = sqrt(3)R/2 to have a total area of 4π R^2

To get the random point, do the following steps:

1) For choosing theta (it can be uniformly chosen over its entire range, since all values of theta are equally likely):

Choose a random "theta" in range [0, 2π]

2) For choosing phi (Use same strategy as for areas, but this time it will apply by pairing circumference of circle for a phi value 2π R cos(phi) with its complement to produce total length of 2π R):

Choose a random "a" in range [0, π/2]

Choose a random "b" in range [0, 1]

Choose a random "c" to be either 1 or -1 (up or down)

if b <= cos(a) then
phi = c a
else
phi = c inv_cos(1 - cos(a))


3) For choosing r (Pair up every possible radius value with its complement, in such a way that the sum of the surface areas of the spheres will be equal to 4π R^2):

Choose a random "d" in range [0, 1]

Choose a random "e" in range [0, 1]

One r value is "d R" the complementary r value is "sqrt(1 - d^2) R".

Total area = 4π (d R)^2 + 4π (sqrt(1 - d^2) R)^2

Total area = 4π (d^2 R^2 + (1 - d^2) R^2)

Total area = 4π (d^2 R^2 + R^2 - d^2 R^2)

Total area = 4π R^2 (as we wanted)

if e <= d^2 then
r = d R
else
r = sqrt(1 - d^2) R


Your random point uniformly chosen inside the sphere in radial coordinates is:

(r, theta, phi)