3
\$\begingroup\$

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()
\$\endgroup\$
6
\$\begingroup\$

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()
\$\endgroup\$
5
\$\begingroup\$

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)])
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.