Generate random unit vectors around circle

Question

I'm trying to generate a bunch of uniformly distributed unit vectors around the unit circle. Here's my code, which is working, but unfortunately I have a for-loop. How do I get rid of this for-loop?

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
def gen_rand_vecs(dims, number):
 vecs = np.random.uniform(low=-1, size=(number,dims))
 mags = np.sqrt((vecs*vecs).sum(axis=-1))
 # How to distribute the magnitude to the vectors
 for i in range(number):
 vecs[i,:] = vecs[i, :] / mags[i]
 return vecs
theta = np.linspace(0, 2*np.pi, 100)
circle = np.array([np.cos(theta), np.sin(theta)])
fig = plt.figure()
plt.plot(circle[0], circle[1])
rand_vecs = gen_rand_vecs(2, 100)
for e in rand_vecs:
 plt.plot([0,e[0]], [0,e[1]], 'r')
plt.show()

Answer 1

To answer your question, you need to add a new dimension to the ndarray:

vecs /= mags[..., np.newaxis]

However

uniformly distributed unit vectors around the unit circle

No it's not, at least not in \$\theta\$. You're generating uniformly distributed points on the unit n-sphere and modifying it to the unit circle; effectively reducing it to an angle. These angles will not be uniformly distributed, and this is easiest to show in 2D:

The area of an angular slice is not invariant.

Notice how the corner piece is larger than the side piece, despite both being \30ドル^{\circ}\$.

Instead, generate the original points with a normal distribution.

This is because the standard normal distribution has a probability density function of

$$ \phi(x) = \frac{1}{\sqrt{2\pi}}\exp\left(-x^2/2\right) $$

so a 2D one has probability density

$$ \phi(x, y) = \phi(x) \phi(y) = \frac{1}{2\pi}\exp\left(-(x^2 + y^2)/2\right) $$

which can be expressed solely in terms of the distance from the origin, \$r = \sqrt{x^2 + y^2}\$:

$$ \phi(r) = \frac{1}{2\pi}\exp\left(-r^2/2\right) $$

This formula applies for any number of dimensions (including \1ドル\$). This means that the distribution is independent of rotation (in any axis) and thus must be evenly distributed along the surface of an n-sphere.

Thanks to Michael Hardy for that explanation.

This is as simple as using instead

np.random.normal(size=(number,dims))

Generating mags can also be just

np.linalg.norm(vecs, axis=-1)

A few minor changes then gets one to

import numpy as np
import matplotlib.pyplot as plt
def gen_rand_vecs(dims, number):
 vecs = np.random.normal(size=(number,dims))
 mags = np.linalg.norm(vecs, axis=-1)
 return vecs / mags[..., np.newaxis]
def main():
 fig = plt.figure()
 for e in gen_rand_vecs(2, 1000):
 plt.plot([0, e[0]], [0, e[1]], 'r')
 plt.show()
main()

The plotting is quite painful, so here's a cleaner, faster version (credit HYRY):

import numpy
from numpy import linalg, newaxis, random
from matplotlib import collections, pyplot
def gen_rand_vecs(dims, number):
 vecs = random.normal(size=(number,dims))
 mags = linalg.norm(vecs, axis=-1)
 return vecs / mags[..., newaxis]
def main():
 ends = gen_rand_vecs(2, 1000)
 # Add 0 vector to start
 vectors = numpy.insert(ends[:, newaxis], 0, 0, axis=1)
 figure, axis = pyplot.subplots()
 axis.add_collection(collections.LineCollection(vectors))
 axis.axis((-1, 1, -1, 1))
 pyplot.show()
main()

Answer 2

If you just want to do this for the unit circle, you can use numpy.random.uniform to choose the angle \$\theta\$ defining the unit vector.

If you want to deal with the more general case of the \$n\$-sphere, there is a method mentioned on Wikipedia.

To generate uniformly distributed random points on the \$(n − 1)\$-sphere (i.e., the surface of the \$n\$-ball), Marsaglia (1972) gives the following algorithm.

Generate an \$n\$-dimensional vector of normal deviates (it suffices to use \$N(0, 1)\,ドル although in fact the choice of the variance is arbitrary), \$\mathbf{x}=(x_1,x_2,\ldots,x_n)\$.

Now calculate the "radius" of this point, \$r=\sqrt{x_1^2+x_2^2+\cdots+x_n^2}\$.

The vector \$\frac{1}{r} \mathbf{x}\$ is uniformly distributed over the surface of the unit n-ball.

My attempt at translating that into Python would be

def get_rand_vec(dims):
 x = np.random.standard_normal(dims)
 r = np.sqrt((x*x).sum())
 return x / r
def gen_rand_vecs(dims, number):
 return [get_rand_vec(dims) for _ in xrange(number)]

---

@Veedrac's answer alerted me to the existence of np.linalg.norm, which gives us

def get_rand_vec(dims):
 x = np.random.standard_normal(dims)
 return x / np.linalg.norm(x)

Or

def gen_rand_vecs(dims, number):
 return map(lambda x: x / np.linalg.norm(x),
 [np.random.standard_normal(dims) for _ in xrange(number)])