Question: 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](http://mathworld.wolfram.com/SpherePointPicking.html) 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.

Answer: 
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)
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)