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: 
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](http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform) 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$.
Nate and Kevin already answered the two I knew. Recalling [this](https://math.stackexchange.com/a/73358/312) and [this](http://mathcentral.uregina.ca/QQ/database/QQ.09.99/wilkie1.html), 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.