How to generate uniformly distributed points in the 3-d unit ball?

Question

I have posted a previous question, this is related but I think it is better to start another thread. This time, I am wondering how to generate uniformly distributed points inside the 3-d unit sphere and how to check the distribution visually and statistically too? I don't see the strategies posted there directly transferable to this situation.

Please help. Many thanks.

Answer 1

The easiest way is to sample points uniformly in the corresponding hypercube and discard those that do not lie within the sphere. In 3D, this should not happen that often, about 50% of the time. (Volume of the hypercube is 1, volume of the sphere is $\frac{4}{3}\pi r^3 = 0.523...$.)

Answer 2

You can also do this in spherical coordinates, in which case there is no rejection. First you generate the radius and the two angles at random, then you use the transition formula to recover $x$, $y$ and $z$ ($x = r \sin \theta \cos \phi$, $y = r \sin \theta \sin \phi$, $z = r \cos \theta$).

You generate $\phi$ unifomly between $0$ and $2\pi$. The radius $r$ and the inclination $\theta$ are not uniform though. The probability that a point is inside the ball of radius $r$ is $r^3$ so the probability density function of $r$ is $3 r^2$. You can easily check that the cubic root of a uniform variable has exactly the same distribution, so this is how you can generate $r$. The probability that a point lies within a spherical cone defined by inclination $\theta$ is $(1-\cos\theta)/2$ or $1 - (1-\cos (-\theta))/2$ if $\theta > \pi/2$. So the density $\theta$ is $sin(\theta)/2$. You can check that minus the arccosine of a uniform variable has the proper density.

Or more simply, we can simulate the cosine of $\theta$ uniformly beteen $-1$ and $1$.

In R this would look as shown below.

n <- 10000 # For example n = 10,000.
phi <- runif(n, max=2*pi)
r <- runif(n)^(1/3)
cos_theta <- runif(n, min=-1, max=1)
x <- r * sqrt(1-cos_theta^2) * cos(phi)
y <- r * sqrt(1-cos_theta^2) * sin(phi)
z <- r * cos_theta

In the course of writing and editing this answer, I realized that the solution is less trivial that I thought.

I think that the easiest and computationally most efficient method is to follow @whuber's method to generate $(x,y,z)$ on the unit sphere as shown on this post and scale them with $r$.

xyz <- matrix(rnorm(3*n), ncol=3)
lambda <- runif(n)^(1/3) / sqrt(rowSums(xyz^2))
xyz <- xyz*lambda