8

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?

Sameer Patel
  • 927
  • 3
  • 16
  • 20
  • How to generate uniformly distributed points at random on an N-sphere: http://en.wikipedia.org/wiki/N-sphere#Uniformly_at_random_from_the_.28n.C2.A0.E2.88.92.C2.A01.29-sphere – unutbu Apr 08 '13 at 13:38
  • 1
    wait, you want the points to be on a sphere, but uniformly distributed in x,y,z,(4th dimension)? that doesn't add up for me. I don't think that points uniformly distributed on a sphere would map to uniformly distributed in 4-space. – Nicu Stiurca Apr 08 '13 at 14:00
  • @SchighSchagh so you can't run monte carlo simulations in 4 dimensions? – Sameer Patel Apr 08 '13 at 14:08
  • @SameerPatel This doesn't have anything to do with Monte Carlo or any other sampling method. There are two different spaces here, (one is R^4, the other is the surface of the 4-sphere), and we need to know with respect to which you want to have a uniformly-at-random distribution. – Nicu Stiurca Apr 08 '13 at 17:41

4 Answers4

10

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

enter image description here

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

unutbu
  • 842,883
  • 184
  • 1,785
  • 1,677
  • Taking a uniform distribution would yield a sphere sampled with a density of points that is higher towards the corners of the cube it would occupy. How to prove that using a Gaussian distribution wouldn't cause this problem, too? It seems plausible if I picture it mentally, but is it really, and why? I assume numpy's Gaussian distributions to be independent for each coordinate when specifying a size. – Guillaume Chevalier Nov 15 '17 at 22:19
  • And what about divisions by zero? – Guillaume Chevalier Nov 15 '17 at 22:42
  • 1
    @GuillaumeChevalier: I confess I do not understand the details, but perhaps you will find [this proof outline useful](https://stats.stackexchange.com/questions/7977/how-to-generate-uniformly-distributed-points-on-the-surface-of-the-3-d-unit-sphe#comment13055_7984). On the issue of division by zero: I believe the probability of the divisor being 0 is itself 0. On a practical level, NumPy handles division by zero by issuing a *warning* and returning `NaN` (and Matplotlib skips points equal to NaN). – unutbu Nov 15 '17 at 22:58
  • That's a great proof! Thanks. – Guillaume Chevalier Nov 18 '17 at 11:08
  • how could we write `points` to a txt file? i used `np.savetxt('sphere.csv', points, delimiter=',')` but it saves them as three big rows instead of 3 columns – Snow May 04 '18 at 00:54
  • Greate solution, but how to change the center of the sphere and maximum radius? – user3379482 Sep 29 '18 at 07:53
  • @user3379482: You can shift `points` to a new location and with a new radius by using `points = new_radius * points + new_location`, where `new_radius` is a scalar, and `new_location` is a sequence (such as a tuple, list, or NumPy array). – unutbu Sep 29 '18 at 12:00
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)
Davide Fiocco
  • 5,350
  • 5
  • 35
  • 72
Manishearth
  • 14,882
  • 8
  • 59
  • 76
  • Sure, this will generate a random point on a 4-sphere, but is the distribution uniform? – Nicu Stiurca Apr 08 '13 at 17:40
  • @SchighSchagh: Uniformly distributed? Yes. – Manishearth Apr 08 '13 at 17:44
  • 1
    x,y,z,w are initially uniformly at random with respect to R^4, but then they undergo a non-linear transform, and it's still not clear to me if OP wants uniformly at random with respect to the surface of the sphere or with respect to R^4. EDIT: can you specify with respect to what you claim uniformly at random, and prove it? – Nicu Stiurca Apr 08 '13 at 17:54
  • @SchighSchagh: Oops, I forgot to normally distribute them. See http://enwp.org/N-sphere#Uniformly_at_random_from_the_.28n.C2.A0.E2.88.92.C2.A01.29-sphere. – Manishearth Apr 08 '13 at 17:59
  • What about divisions by zero? – Guillaume Chevalier Nov 15 '17 at 22:42
1

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.

Guillaume Chevalier
  • 9,613
  • 8
  • 51
  • 79
1

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 U1/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:

sphere

Rabbid76
  • 202,892
  • 27
  • 131
  • 174
Ron Klintskey
  • 9
  • 7
  • 1
    What are the variables `distr` and `sample_isotropic`? Could you modify the code example so that it works? – n1k31t4 Aug 13 '19 at 14:58