# Generate random unit vectors around circle

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

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:

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.axis((-1, 1, -1, 1))

pyplot.show()

main()

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