## uniform on the sphere [or not]

**W**hile looking at X validated questions, I came upon this comment that simulating a uniform distribution on a d-dimensional unit sphere does not proceed from generating angles at random on (0,2π) and computing spherical coordinates… Which I must confess would have been my initial suggestion! This is obvious, nonetheless, when computing the Jacobian of the spherical coordinate transform, which involves powers of the sines of the angles, in a decreasing sequence from d-1 to zero. This means that the angles should be simulated according to their respective sine-power densities. However, except for the d=3 case, where simulating from the density sin(φ) is straightforward by inverse cdf, i.e. φ=acos(1-2u), the cdfs for the higher powers are combinations of sines and cosines, and as such are not easily inverted. Take for instance the eighth power:

F⁸(φ)=(840 φ – 672 sin(2 φ) + 168 sin(4 φ) – 32 sin(6 φ) + 3 sin(8 φ))/3072

While the densities are bounded by sin(φ), up to a constant, and hence an accept-reject can be easily derived, the efficiency decreases with the dimension according to the respective ratio of the Wallis’ integrals, unsurprisingly. A quick check for d=4 shows that the Normal simulation+projection-by-division-by-its-norm is faster.

Puzzling a bit further about this while running, I wondered at the simultaneous simulations from sin(φ), sin(φ)², sin(φ)³, &tc., but cannot see a faster way to recycle simulations from sin(φ). Points (φ,u) located in-between two adjacent power curves are acceptable simulations from the corresponding upper curve but they need be augmented by points (φ,u) under the lower curve to constitute a representative sample. In the end, this amounts to multiplying simulations from the highest power density as many times as there are powers. No gain in sight… Sigh!

However, a few days later, while enjoying the sunset over Mont Blanc(!), I figured out that there exists a direct and efficient way to simulate from these powers of the sine function. Indeed, when looking at the density of cos(φ), it happens to be the signed root of a Beta(½,(d-1)/2), which avoids the accept-reject step. Presumably this is well-known, but I have not seen this proposal associated with the uniform distribution on the sphere.

March 13, 2018 at 5:54 pm

I realize your real question is how to generate randoms from a pdf which is powers of sines.

But as for the proximal question of how to generate randoms on a hypersphere, here’s how I did it (unsure of how this will render in comments… the lines in the original were up to 95 char wide):

nextVector <- function(N) { # Uniform on 1st 2^N-ant of sphere

## Method of Muller & Marsaglia for generating uniform randoms on the hypersphere:

##

## Muller, M. E. "A Note on a Method for Generating Points Uniformly on N-Dimensional

## Spheres", Comm ACM 2 (1959), 19-20.

## Marsaglia, G. "Choosing a Point from the Surface of a Sphere", Ann Math Stat 43 (1972),

## 645-646. Marsaglia's email is: George Marsaglia .

## See also: Knuth, _Seminumerical Algorithms_, 2nd. ed., pp. 130-131.

## See also: http://mathworld.wolfram.com/HyperspherePointPicking.html

dotProduct <- function(x, y) { drop(crossprod(x, y)) } # Dot product of 2 vectors

x 1 as N -> +oo, and you waste a lot of work!

} #

March 13, 2018 at 5:56 pm

Ok, looks like the code was totally mangled by WordPress. Here’s the gist of it, for an N-dimensional hypersphere:

x <- rnorm(N)

abs(x / sqrt(drop(crossprod(x, x))))

March 13, 2018 at 6:44 pm

Thank you. I am aware of this proposal and so is the person who asked the initial question on X validated.

March 9, 2018 at 2:19 pm

Hello,

You may be interested in these papers:

Marsaglia G (1972) Choosing a point from the surface of a sphere.

The Annals of Mathematical Statistics, 43, 645-646

http://chalkdustmagazine.com/blog/sylvesters-convex-hull-problem/

Using Marsaglia’s method, this is a function in R to generate points in the d-dimensional unit ball.

gen.dBall<- function(npts, d)

{

### generate npts points within the d-dimensional unit ball

matd<- matrix( rnorm( npts*d), ncol=d)

radii<- runif( npts)^(1/d)

d0<- apply(matd, 1, function(x){sqrt(sum(x^2))})

matd<- radii * matd/d0

invisible( matd)

}

Best wishes,

Mario

March 9, 2018 at 6:24 pm

Dear Mario, thanks for the comment. I was aware of George Marsaglia’s paper and method. As stressed in other replies, I am not claiming any novelty in my post or X validated answer, just a resolution to a particular puzzle that I found amusing!

March 8, 2018 at 4:43 pm

You can generate points (x,y,z) uniformly from the surface of a sphere as illustrated below:

n = 500

a = runif(n,0,2*pi)

z = runif(n,-1,1)

r = sqrt(1-z^2)

x = r*cos(a)

y = r*sin(a)

I’ll leave the proof, and possible generalization to higher dimensions, as an exercise for the reader.

March 9, 2018 at 6:26 pm

Thank you, Radford. The property that one coordinate is marginally uniform is quite helpful to build a recursive algorithm, recursive in the dimension that is.

March 8, 2018 at 3:23 pm

Pourquoi s’embeter? Si sont iid N(0,1) alors est uniforme sur la sphere

March 8, 2018 at 3:56 pm

Je sais, et c’est sûrement plus efficace!, mais bon l’auteur de la question voulait une façon de simuler directement en dimension n-1…

March 8, 2018 at 1:28 am

The reverse transformation — from polynomials to powers of trig functions — is part of the process of expressing the Beta function in terms of Gamma functions.

March 8, 2018 at 7:10 am

Thanks! I deem the simulation method to be less efficient that the x/|x| approach.