# 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. 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. 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
So, what is the probability density function (pdf) of such a distribution (say we are in $n$ dimensions). Somehow the gamma function is introduced. Could you help me a bit? Thanks! – nullgeppetto Jan 30 '15 at 9:59
@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
@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 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 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 Costello May 30 at 1:06

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[0,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

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 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 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 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 at 12:38