sampling integers uniformly efficiently in python using numpy/scipy

Question

I have a problem where depending on the result of a random coin flip, I have to sample a random starting position from a string. If the sampling of this random position is uniform over the string, I thought of two approaches to do it: one using multinomial from numpy.random, the other using the simple randint function of Python standard lib. I tested this as follows:

from numpy import *
from numpy.random import multinomial
from random import randint
import time

def use_multinomial(length, num_points):
    probs = ones(length)/float(length)
    for n in range(num_points):
    result = multinomial(1, probs)

def use_rand(length, num_points):
    for n in range(num_points):
    rand(1, length)

def main():
    length = 1700
    num_points = 50000

    t1 = time.time()
    use_multinomial(length, num_points)
    t2 = time.time()
    print "Multinomial took: %s seconds" %(t2 - t1)

    t1 = time.time()
    use_rand(length, num_points)
    t2 = time.time()
    print "Rand took: %s seconds" %(t2 - t1)    

if __name__ == '__main__':
    main()

The output is:

Multinomial took: 6.58072400093 seconds Rand took: 2.35189199448 seconds

it seems like randint is faster, but it still seems very slow to me. Is there a vectorized way to get this to be much faster, using numpy or scipy?

thanks.

Answer 1

I changed your code to actually return values (and used randint instead of rand - isn't that what you meant?) like this...

def use_multinomial(length, num_points):
    probs = ones(length)/float(length)
    return multinomial(1, probs, num_points)

def use_rand(length, num_points):
    return [randint(1,length) for _ in range(num_points)]

Then I tried my own version, using numpy.random.randint to generate a numpy array of random points on the string:

def use_np_randint(length, num_point):
    return nprandint(1, length, num_points)

The results:

Multinomial took: 13.6279997826 seconds
Rand took: 0.185000181198 seconds
NP randint took: 0.00100016593933 seconds

Multinomial is obviously really slow comparitively, but is that even what you want? I thought you said you wanted a uniform distribution? Using numpy's randint is clearly the fastest of the bunch.