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

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:

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

StackExchange Code Review Q#77927, answer score: 6