Python Uniform distribution of points on 4 dimensional sphere

Question

I need a uniform distribution of points on a 4 dimensional sphere. I know this is not as trivial as picking 3 angles and using polar coordinates.

In 3 dimensions I use

from random import random

u=random()
costheta = 2*u -1 #for distribution between -1 and 1
theta = acos(costheta)
phi = 2*pi*random

x=costheta
y=sin(theta)*cos(phi)
x=sin(theta)*sin(phi)

This gives a uniform distribution of x, y and z.

How can I obtain a similar distribution for 4 dimensions?

Answer 1

A standard way, though, perhaps not the fastest, is to use Muller's method to generate uniformly distributed points on an N-sphere:

import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as axes3d

N = 600
dim = 3

norm = np.random.normal
normal_deviates = norm(size=(dim, N))

radius = np.sqrt((normal_deviates**2).sum(axis=0))
points = normal_deviates/radius

fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
ax.scatter(*points)
ax.set_aspect('equal')
plt.show()

Simply change dim = 3 to dim = 4 to generate points on a 4-sphere.

Answer 2

Take a point in 4D space whose coordinates are distributed normally, and calculate its unit vector. This will be on the unit 4-sphere.

from random import random
import math
x=random.normalvariate(0,1)
y=random.normalvariate(0,1)
z=random.normalvariate(0,1)
w=random.normalvariate(0,1)
r=math.sqrt(x*x + y*y + z*z + w*w)
x/=r
y/=r
z/=r
w/=r
print (x,y,z,w)

Answer 3

I like @unutbu's answer if the gaussian sampling really creates an evenly spaced spherical distribution (unlike sampling from a cube), but to avoid sampling on a Gaussian distribution and to have to prove that, there is a simple solution: to sample on a uniform distribution on a sphere (not on a cube).

1. Generate points on a uniform distribution.
2. Compute the squared radius of each point (avoid the square root).
3. Discard points:
   • Discard points for which the squared radius is greater than 1 (thus, for which the unsquared radius is greater than 1).
   • Discard points too close to a radius of zero to avoid numerical instabilities related to the division in the next step.
4. For each sampled point kept, divide the sampled point by the norm so as to renormalize it the unit radius.
5. Wash and repeat for more points because of discarded samples.

This obviously works in an n-dimensional space, since the radius is always the L2-norm in higher dimensions.

It is fast so as avoiding a square-root and sampling on a Gaussian distribution, but it's not a vectorized algorithm.

Answer 4

I found a good solution for sampling from N-dim sphere. The main idea is:

If Y is drawn from the uncorrelated multivariate normal distribution, then S = Y / ||Y|| has the uniform distribution on the unit d-sphere. Multiplying S by U^1/d, where U has the uniform distribution on the unit interval (0,1), creates the uniform distribution in the unit d-dimensional ball.

Here is the python code to do this:

Y = np.random.multivariate_normal(mean=[0], cov=np.eye(1,1), size=(n_dims, n_samples))
Y = np.squeeze(Y, -1)
Y /= np.sqrt(np.sum(Y * sample_isotropic, axis=0))
U = np.random.uniform(low=0, high=1, size=(n_samples)) ** (1/n_dims)
Y *= distr * radius # in my case radius is one

This is what I get for the sphere: