# How to generate uniformly distributed points on the surface of the 3-d unit sphere?

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$?

-

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:

Here is the diagnostic plot of the distances:

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:

(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:

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

-
@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? –  Qiang Li Mar 8 '11 at 1:12
@whuber what software did you use to generate the sphere plot? –  Berk U. Mar 8 '11 at 1:18
I bet it is mathematica... just tell from the color... :) –  Qiang Li Mar 8 '11 at 1:26
@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. –  cardinal Mar 8 '11 at 1:53
@whuber, (+1), nice answer. –  cardinal Mar 8 '11 at 2:00

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.

-
I am just wondering if your method of generating points (x,y,z) is essentially the same as whuber's method? –  Qiang Li Mar 8 '11 at 1:33
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. –  Henry Mar 8 '11 at 1:46
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) –  whuber Mar 8 '11 at 2:14
@whuber: Indeed. Offtopic, but for anyone interested I have an interactive selection of world map projections here, some of which are equal area –  Henry Mar 8 '11 at 7:10

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

-
is it possible to post the links where I can find the full text of these references? thanks. –  Qiang Li Mar 8 '11 at 0:30
I dont have the paper on me, but this page seems to describe the algorithm (and several others) mlahanas.de/Math/nsphere.htm –  emakalic Mar 8 '11 at 5:38
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? –  Mohit Feb 28 '12 at 6:29

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

-

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.

-
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. –  Qiang Li Mar 7 '11 at 23:35
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 –  Berk U. Mar 8 '11 at 0:01
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. –  whuber Mar 8 '11 at 6:30