11

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.

Community
  • 1
Qiang Li
  • 1,145
  • 2
  • 9
  • 10
  • 4
    The techniques in the previous question apply directly once you observe that the number of points within distance $r$ of the origin must be proportional to $r^3$. Thus if you generate an independent uniform variate $u$ in $[0,1]$ along with a point $w$ on the surface of the sphere, scaling $w$ by $u^{1/3}$ does the trick. – whuber Mar 08 '11 at 16:04
  • @whuber: maybe I just did not get the essence of the previous techniques. Let me try what you described. Additionally, what are the ways to check the uniformity here again? – Qiang Li Mar 08 '11 at 17:10
  • 2
    @Qiang Ripley's K function and chi-squared tests. You could also check the uniformity of the radial projection of the points on the sphere's surface, the uniformity of the cube of the lengths of the points, and the independence of those two. – whuber Mar 08 '11 at 17:18
  • For me, it is not so obvious what "uniformly distributed" means... And probably a try to define it will automagically create a generating algorithm (= –  Mar 08 '11 at 19:33
  • @mbq, I think to define the term, we need to have a p.d.f. of $f_{R, \Theta, \Phi}(r,\theta, \phi)=r^2$. – Qiang Li Mar 08 '11 at 21:12
  • @mbq: correction, I think the p.d.f is $f_{R, \Theta, \Phi}(r,\theta, \phi)=\frac{3}{4\pi}r^2sin\theta$ – Qiang Li Mar 08 '11 at 23:51
  • @whuber: is the formula above correct? – Qiang Li Mar 09 '11 at 00:25
  • @whuber another superfast technique is to generate spherically symmetric multinormals in five dimensions, normalize the vectors (so they're uniform on the surface of the 5-sphere), then drop any two coordinates from all the vectors (the same two from all vectors, but any two will do). This method generalizes to large dimensions, wherein avoiding powers of 1/d is a good idea. It works because of Archimedes' Hat-Box theorem (http://mathworld.wolfram.com/ArchimedesHat-BoxTheorem.html), generalized to any number of dimensions. – Reb.Cabin Apr 04 '19 at 13:21
  • @Reb.Cabin Very nice! – whuber Apr 04 '19 at 13:25

3 Answers3

15

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 than 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
Community
  • 1
gui11aume
  • 13,383
  • 2
  • 44
  • 89
  • 3
    This is a much better answer due to the lack of rejection. In high dimensional spaces, rejection sampling can be very costly due to the low probability of acceptance. – kingledion Feb 01 '17 at 22:37
  • 2
    The last bit of code can be adapted to higher dimension, say `d`. For this, replace all instances of `3` by `d`. – gui11aume Feb 01 '17 at 23:11
15

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...$.)

bayerj
  • 12,735
  • 3
  • 35
  • 56
  • +1. This is one of the techniques recommended by the comp.graphics.algorithms FAQ ["Uniform random points on sphere".](http://www.cgafaq.info/wiki/Uniform_random_points_on_sphere) – David Cary Mar 17 '11 at 02:51
  • 1
    What if we want to do that for $n > 100$? – ares Oct 04 '18 at 09:35
  • 3
    This is called the "rejection method." While working well in three dimensions, by twenty-seven dimensions, only one in a trillion points lies in the 27-ball and not in the rest of the 27-cube, so the rejection method doesn't generalize well. I mention this because I currently need samples uniformly in a ball of 2,440 dimensions. – Reb.Cabin Apr 04 '19 at 13:24
1

In my opinion, the easiest option which also generalizes to higher dimensional balls (which is not the case of spherical coordinates and even less the case of rejection sampling) is to generate random points $P$ that are products of two random variables $P = N/||N|| * U^{1/n}$ where $N$ is a Gaussian random variable (i.e. isotropic, i.e. pointing in any direction uniformly) normalized so that it lies on the sphere and $U$ which is a uniform random variable in $[0,1]$ to the power $1/n$, $n$ being the dimensionality of the data, taking care of the radius.

Et voilà!

Jean-Luc Bouchot
  • 141
  • 2