68
$\begingroup$

I am wondering how to generate uniformly distributed points on the surface of the 3-d unit sphere? Also after generating those points, what is the best way to visualize and check whether they are truly uniform on the surface $x^2+y^2+z^2=1$?

$\endgroup$
  • $\begingroup$ If by uniform you mean "regular", there is no way to do so outside of $n$=2, 4, 6, 8, 12, 20. $\endgroup$ – Marcos Feb 6 '16 at 21:10
  • 1
    $\begingroup$ what's wrong with sample from a MultiVariateGaussian and that vector just normalize it: v = MultivariateNormal(torch.zeros(10000), torch.eye(10000)) and then v = v/v.norm(10000) $\endgroup$ – Pinocchio Apr 1 '18 at 1:24
72
$\begingroup$

A standard method is to generate three standard normals and construct a unit vector from them. That is, when $X_i \sim N(0,1)$ and $\lambda^2 = X_1^2 + X_2^2 + X_3^2$, then $(X_1/\lambda, X_2/\lambda, X_3/\lambda)$ is uniformly distributed on the sphere. This method works well for $d$-dimensional spheres, too.

In 3D you can use rejection sampling: draw $X_i$ from a uniform$[-1,1]$ distribution until the length of $(X_1, X_2, X_3)$ is less than or equal to 1, then--just as with the preceding method--normalize the vector to unit length. The expected number of trials per spherical point equals $2^3/(4 \pi / 3)$ = 1.91. In higher dimensions the expected number of trials gets so large this rapidly becomes impracticable.

There are many ways to check uniformity. A neat way, although somewhat computationally intensive, is with Ripley's K function. The expected number of points within (3D Euclidean) distance $\rho$ of any location on the sphere is proportional to the area of the sphere within distance $\rho$, which equals $\pi\rho^2$. By computing all interpoint distances you can compare the data to this ideal.

General principles of constructing statistical graphics suggest a good way to make the comparison is to plot variance-stabilized residuals $e_i(d_{[i]} - e_i)$ against $i = 1, 2, \ldots, n(n-1)/2=m$ where $d_{[i]}$ is the $i^\text{th}$ smallest of the mutual distances and $e_i = 2\sqrt{i/m}$. The plot should be close to zero. (This approach is unconventional.)

Here is a picture of 100 independent draws from a uniform spherical distribution obtained with the first method:

100 uniform spherical points

Here is the diagnostic plot of the distances:

Diagnostic plot

The y scale suggests these values are all close to zero.

Here is the accumulation of 100 such plots to suggest what size deviations might actually be significant indicators of non-uniformity:

Simulated values

(These plots look an awful lot like Brownian bridges...there may be some interesting theoretical discoveries lurking here.)

Finally, here is the diagnostic plot for a set of 100 uniform random points plus another 41 points uniformly distributed in the upper hemisphere only:

Simulated non-uniform values

Relative to the uniform distribution, it shows a significant decrease in average interpoint distances out to a range of one hemisphere. That in itself is meaningless, but the useful information here is that something is non-uniform on the scale of one hemisphere. In effect, this plot readily detects that one hemisphere has a different density than the other. (A simpler chi-square test would do this with more power if you knew in advance which hemisphere to test out of the infinitely many possible ones.)

$\endgroup$
  • 1
    $\begingroup$ @whuber: very nice! thanks a lot for your post! "$(X_1/\lambda, X_2/\lambda, X_3/\lambda)$ is uniformly distributed on the sphere." where can I find reference about its proof, or is it simply provable? $\endgroup$ – Qiang Li Mar 8 '11 at 1:12
  • 25
    $\begingroup$ @Qiang, Here's the essence of the proof: Let $X \sim \mathcal{N}(0, I_n)$ where $I_n$ denotes the $n \times n$ identity matrix. Then for any orthogonal matrix $Q$, $Q X \sim N(0, I_n)$. Hence the distribution of $X$ is invariant under rotations. Let $Y = X / \|X\|_2$ and note that $Y_Q = Q X / \|Q X\|_2 = Q X / \|X\|_2$ for any orthogonal $Q$. Since $X$ is invariant to rotations, so is $Y$, and since $\|Y\|_2 = 1$ almost surely, then it must be uniformly distributed on the sphere. $\endgroup$ – cardinal Mar 8 '11 at 1:53
  • 3
    $\begingroup$ @Mike No, because a uniform distribution of the latitude $\phi$ does not give a uniform distribution on the sphere. (Most of the sphere's surface lies at lower latitudes near the Equator away from the Poles. You need a uniform distribution of $\cos(\phi)$ instead.) $\endgroup$ – whuber Jul 18 '13 at 20:38
  • 1
    $\begingroup$ @Ahsan Because the orthogonal matrices form a transitive group of area-preserving transformations of the sphere, the distribution is uniform over the subset of the sphere of the form $X/||X||_2$: but this is the entire sphere. $\endgroup$ – whuber Jun 15 '15 at 18:51
  • 1
    $\begingroup$ @Cesar "Uniform distribution" (on the sphere). $\endgroup$ – whuber Sep 19 '16 at 12:29
19
$\begingroup$

Here is some rather simple R code

n     <- 100000                  # large enough for meaningful tests
z     <- 2*runif(n) - 1          # uniform on [-1, 1]
theta <- 2*pi*runif(n) - pi      # uniform on [-pi, pi]
x     <- sin(theta)*sqrt(1-z^2)  # based on angle
y     <- cos(theta)*sqrt(1-z^2)     

It is very simple to see from the construction that $x^2+y^2 = 1- z^2$ and so $x^2+y^2+z^2=1$ but if it needs to be tested then

mean(x^2+y^2+z^2)  # should be 1
var(x^2+y^2+z^2)   # should be 0

and easy to test that each of $x$ and $y$ are uniformly distributed on $[-1,1]$ ($z$ obviously is) with

plot.ecdf(x)  # should be uniform on [-1, 1]
plot.ecdf(y)
plot.ecdf(z)

Clearly, given a value of $z$, $x$ and $y$ are uniformly distributed around a circle of radius $\sqrt{1-z^2}$ and this can be tested by looking at the distribution of the arctangent of their ratio. But since $z$ has the same marginal distribution as $x$ and as $y$, a similar statement is true for any pair, and this too can be tested.

plot.ecdf(atan2(x,y)) # should be uniform on [-pi, pi]
plot.ecdf(atan2(y,z))
plot.ecdf(atan2(z,x))

If still unconvinced, the next steps would be to look at some arbitrary 3-D rotation or how many points fell within a given solid angle, but that starts to get more complicated, and I think is unnecessary.

$\endgroup$
  • $\begingroup$ I am just wondering if your method of generating points (x,y,z) is essentially the same as whuber's method? $\endgroup$ – Qiang Li Mar 8 '11 at 1:33
  • 3
    $\begingroup$ No it isn't: whuber uses three random numbers while I use two. Mine is a special case of "generate a point on $[-1,1]$ with suitable density [proportional to $(1-z^2)^{n/2-1}$] and then step down a dimension". Here conveniently $n=2$ as this is formally a 2-sphere. $\endgroup$ – Henry Mar 8 '11 at 1:46
  • 3
    $\begingroup$ Or, more generally, generate uniform points on the map using any equal-area projection (yours is a cylindrical equal-area one), then project back. (+1) $\endgroup$ – whuber Mar 8 '11 at 2:14
  • $\begingroup$ @whuber: Indeed. Offtopic, but for anyone interested I have an interactive selection of world map projections here, some of which are equal area $\endgroup$ – Henry Mar 8 '11 at 7:10
  • 2
    $\begingroup$ This is pretty much the standard approach used in computer graphics, based Archimedes' Hat-Box Theorem: mathworld.wolfram.com/ArchimedesHat-BoxTheorem.html $\endgroup$ – Edward KMETT Sep 11 '16 at 16:40
10
$\begingroup$

If you want to sample points uniformly distributed on the 3D sphere (i.e., the surface of a 3D ball), use a simple rejection, or the method of Marsaglia (Ann. Math. Statist., 43 (1972), pp. 645–646). For low dimensions, the rejection ratio is quite low.

If you want to generate random points from higher-dimensional spheres and balls, then it depends on the purpose and scale of the simulation. If you do not want to perform large simulations, use the method of Muller (Commun. ACM, 2 (1959), pp. 19–20) or its "ball" version (see the paper of Harman & Lacko cited above). That is:

to get a sample uniformly distributed on an n-sphere (surface) 1) generate X from n-dimensional standard normal distribution 2) divide each component of X by the Euclidean norm of X

to get a sample uniformly distributed on an n-ball (interior) 1) generate X from (n+2)-dimensional standard normal distribution 2) divide each component of X by the Euclidean norm of X and take only first n components

If you want to perform large simulations, then you should investigate more specialized methods. Upon request, I can send you the paper of Harman and Lacko on conditional distribution methods, which provides the classification and generalizations of some algorithms mentioned in this discussion. The contact is available at my website (http://www.iam.fmph.uniba.sk/ospm/Lacko)

If you want to check, whether yout points are truly uniform on the surface or interior of a ball, look at the marginals (all should by the same, because of the rotational invariance, the squared norm of a projected sample is beta distributed).

$\endgroup$
  • $\begingroup$ what's wrong with sample from a MultiVariateGaussian and that vector just normalize it: v = MultivariateNormal(torch.zeros(10000), torch.eye(10000)) and then v = v/v.norm(10000) $\endgroup$ – Pinocchio Apr 1 '18 at 1:27
8
$\begingroup$

I had a similar problem (n-sphere) during my PhD and one of the local 'experts' suggested rejection sampling from a n-cube! This, of course, would have taken the age of the universe as I was looking at n in the order of hunderds.

The algorithm I ended up using is very simple and published in:

W.P. Petersen and A. Bernasconic Uniform sampling from an n-sphere: Isotropic method Technical Report, TR-97-06, Swiss Centre for Scientific Computing

I also have this paper in my bibliography that I havent looked at. You may find it useful.

Harman, R. & Lacko, V. On decompositional algorithms for uniform sampling from $n$-spheres and $n$-balls Journal of Multivariate Analysis, 2010

$\endgroup$
  • $\begingroup$ is it possible to post the links where I can find the full text of these references? thanks. $\endgroup$ – Qiang Li Mar 8 '11 at 0:30
  • $\begingroup$ I dont have the paper on me, but this page seems to describe the algorithm (and several others) mlahanas.de/Math/nsphere.htm $\endgroup$ – emakalic Mar 8 '11 at 5:38
  • 3
    $\begingroup$ As I understand, (from the paper of Petersen and Bernasconic) for a d-dimensional ball, one can generate the radius by raising a U(0,1) variate to the (1/d) power, and the last angle as a U(0,2$\pi$) variate. The intermediate angles can be obtained as $C. asin(\sqrt[k]{u})$, where ${C^{-1}}$ is $$\frac{\sqrt{\pi} \Gamma(\frac{k}{2} + 0.5)}{\Gamma(\frac{k}{2} + 1)}$$. To me this sounds rather simple. What I'm wondering is this: if I use a quasi random sequence for my uniforms, will I get the niceness in the ball as well? $\endgroup$ – Mohit Feb 28 '12 at 6:29
3
$\begingroup$

I have had this problem before, and here is an alternative I found,

As for the distribution itself, the formula I found that works decently is to use polar coordinates (I actually use a variation of poler coordinates that developed), then convert to Cartesian coordinates.

The radius is of course the radius of the sphere on which you are plotting . Then you have the second value for angle on the flat plane, followed by the third value which is the angle above or below that plane.

To get a decent distribution, assume that U is a uniformly distributed random number, r is radius, a is the second polar coordinate, and b is the third polar coordinate,

a=U*360 b=U+U-1 then convert to cartesian via x=r*sin(b)sin(a) z=rsin(b)cos(a) y=rsin(b)

I recently found the following which is better mathematically speaking, a=2(pi)*U b=cos^-1(2U-1)

Not much different from my original formula actually, though mine is degrees vs radians.

This recent version supposedly can be used for hyperspheres, though no mention was made on how to achieve it.

Though I check the uniformity visually by the rather cheap method of making maps for Homeworld 2 and then "playing" those maps. In fact, because the maps are made with lua scripts, you can build your formula right into the map and thus check multiple samplings without ever leaving the game. Not scientific perhaps, but is a good method for visually seeing the the results.

$\endgroup$
2
$\begingroup$

Here is the pseudocode:

  1. $v \sim MultiVariateGaussian(\mu,\sigma I)$
  2. $v = \frac{v}{ \| v \| }$

In pytorch:

v = MultivariateNormal(torch.zeros(10000), torch.eye(10000))
v = v/v.norm(2)

I don't understand this well enough but I've been told by whuber that:

v = torch.normal(torch.zeros(10000), torch.eye(10000))
v = v/v.norm(2)

is also correct i.e. sampling from a univariate normal for each coordinate.

$\endgroup$
0
$\begingroup$

My best guess would be to first generate a set of uniformly distributed points in 2 dimensional space and to then project those points onto the surface of a sphere using some sort of projection.

You will probably have to mix and match the way you generate the points with the way that you map them. In terms of the 2D point generationgeneration, I think that scrambled low-discrepancy sequences would be a good place to start (i.e. a scrambled Sobol sequence) since it usually produces points that are not "clumped together". I'm not as sure about which type of mapping to use, but Woflram popped up the Gnonomic projection... so maybe that could work?

MATLAB has a decent implementation of low discrepancy sequences which you can generate using q = sobolset(2) and scramble using q = scramble(q). There is also a mapping toolbox in MATLAB with a bunch of different projection functions you could use in case you did not want to code the mapping and graphics yourself.

$\endgroup$
  • 1
    $\begingroup$ can any of these projections still preserve the uniformity of randomness? Again, how can I check whether the final distribution of these points are truly uniformly distributed on the surface of the sphere? Thanks. $\endgroup$ – Qiang Li Mar 7 '11 at 23:35
  • $\begingroup$ Sorry I was just speaking hypothetically... I think the mapping functions on MATLAB would allow you to check that since they have some visualizations embedded in them. If not, I also found a nice website that talks about how to generate uniformly distributed points on a sphere in 3D using things like randomized angles etc. They have some C code on there too. Take a look $\endgroup$ – Berk U. Mar 8 '11 at 0:01
  • 3
    $\begingroup$ Uniform random points on a gnomonic projection will not be uniform on the sphere, because the gnomonic is not equal area. The projection proposed by Henry, $(\lambda, \phi)$ --> $(\lambda, \sin(\phi))$ (from longitude-latitude to a rectangle in $\mathbb{R}^2$), is equal-area. $\endgroup$ – whuber Mar 8 '11 at 6:30

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.