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

Community · Accepted Answer · 2017-05-23 12:41:02Z

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

mjolka · Accepted Answer · 2015-01-19 04:21:52Z

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

Stack Exchange Network

Generate random unit vectors around circle

2 Answers 2

You must log in to answer this question.

Hot Network Questions

Generate random unit vectors around circle

2 Answers 2

You must log in to answer this question.

Related

Hot Network Questions