Generate random unit vectors around circle
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()