11
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – whuber Mar 8 '11 at 16:04
  • $\begingroup$ @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? $\endgroup$ – Qiang Li Mar 8 '11 at 17:10
  • 2
    $\begingroup$ @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. $\endgroup$ – whuber Mar 8 '11 at 17:18
  • $\begingroup$ For me, it is not so obvious what "uniformly distributed" means... And probably a try to define it will automagically create a generating algorithm (= $\endgroup$ – user88 Mar 8 '11 at 19:33
  • $\begingroup$ @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$. $\endgroup$ – Qiang Li Mar 8 '11 at 21:12
14
$\begingroup$

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

$\endgroup$
  • $\begingroup$ +1. This is one of the techniques recommended by the comp.graphics.algorithms FAQ "Uniform random points on sphere". $\endgroup$ – David Cary Mar 17 '11 at 2:51
  • 1
    $\begingroup$ What if we want to do that for $n > 100$? $\endgroup$ – ares Oct 4 '18 at 9:35
  • 1
    $\begingroup$ 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. $\endgroup$ – Reb.Cabin Apr 4 at 13:24
13
$\begingroup$

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
$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – kingledion Feb 1 '17 at 22:37
  • 2
    $\begingroup$ The last bit of code can be adapted to higher dimension, say d. For this, replace all instances of 3 by d. $\endgroup$ – gui11aume Feb 1 '17 at 23:11
0
$\begingroup$

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à!

$\endgroup$
  • 2
    $\begingroup$ A uniformly distributed radius will not give a uniform point in the ball ... $\endgroup$ – kjetil b halvorsen Jun 26 '18 at 18:45
  • 1
    $\begingroup$ True. You'd need to adapt a little bit the distribution of the $U$ variable to take into account for lower / higher density regions. $\endgroup$ – Jean-Luc Bouchot Jun 26 '18 at 19:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.