Brian Quinlan wrote:
This is less a Python question and more a optimization/probability
question. Imaging that you have a list of objects and there frequency in
a population e.g.
lst = [(a, 0.01), (b, 0.05), (c, 0.50), (d, 0.30), (e, 0.04), (f, 0.10)]
and you want to drawn n items from that list (duplicates allowed), with
that probability distribution.
The fastest algorithm that I have been able to devise for doing so is:
O(n * log(len(lst))). Can anyone think or a solution with a better time
complexity? If not, is there an obviously better way to do this
(assuming n is big and the list size is small).
Any way I tried to slice and dice it, I could not get any faster.
draw2 and draw 3 generate code on the fly. draw4 sneakily tries to
trade memory and accuracy for speed but is even slower!
First the times, then the code:
$ ./timeit.py 'from probDistribution import draw as draw; draw(10000,
[0.01, 0.05, 0.50, 0.30, 0.04, 0.10])'
100 loops, best of 3: 13.4 msec per loop
$ ./timeit.py 'from probDistribution import draw2 as draw; draw(10000,
[0.01, 0.05, 0.50, 0.30, 0.04, 0.10])'
100 loops, best of 3: 15.2 msec per loop
$ ./timeit.py 'from probDistribution import draw3 as draw; draw(10000,
[0.01, 0.05, 0.50, 0.30, 0.04, 0.10])'
100 loops, best of 3: 16.2 msec per loop
$ ./timeit.py 'from probDistribution import draw4 as draw; draw(10000,
[0.01, 0.05, 0.50, 0.30, 0.04, 0.10])'
10 loops, best of 3: 30.5 msec per loop
=== CODE probDistribution.py ===
from random import random, randrange
from bisect import bisect
def draw(n, lst):
ps = []
last = 0
for p in lst:
ps.append(last + p)
last += p
# ps = [0.01, 0.06, 0.56, 0.86, 0.90, 1.00]
chosen = [0] * len(lst) # track frequency
for i in range(n):
r = random()
chosen[bisect(ps, r)] += 1 # binary search and increment
result = [] # rescale probability based on frequency
for c in chosen:
result.append(float(c) / n)
return result
def draw2(n, lst):
""" uses dynamicc code generation of this form:
chosen = [0] * 6
for i in xrange(10000):
r = random()
if r < 0.01: chosen[0]+=1
elif r < 0.06: chosen[1] +=1
...
elif r < 0.90: chosen[4] +=1
else chosen[5]+=1
"""
assert len(lst)>1, "Corner case NOT covered"
codestr = 'chosen = [0] * %i\n' % (len(lst),)
codestr += 'for i in xrange(%i):\n r = random()\n' % (n,)
last = 0.0
lstmax = len(lst)-1
for i,p in enumerate(lst):
last += p
if i==0: codestr += ' if r < %g: chosen[%i] +=1\n' %
(last, i)
elif i==lstmax: codestr += ' else: chosen[%i] +=1\n' % (i,)
else: codestr += ' elif r < %g: chosen[%i] +=1\n' %
(last, i)
exec codestr
result = [] # rescale probability based on frequency
for c in chosen:
result.append(float(c) / n)
return result
def draw3(n, lst):
""" uses dynamicc code generation of this form:
chosen = [0] * 6
for i in xrange(10000):
r = random()
chosen[-1+ (
((r<0.01) and 1)
or ((r<0.06) and 2)
...
or ((r<0.90) and 5)
or 6
)] +=1
"""
assert len(lst)>1, "Corner case NOT covered"
codestr = 'chosen = [0] * %i\n' % (len(lst),)
codestr += 'for i in xrange(%i):\n r = random()\n' % (n,)
codestr += ' chosen[-1+ (\n'
last = 0.0
lstmax = len(lst)-1
for i,p in enumerate(lst):
last += p
if i==0: codestr += ' ((r<%g) and 1)\n' % (last)
elif i==lstmax: codestr += ' or %i\n )] +=1\n' % (i+1,)
else: codestr += ' or ((r<%g) and %i)\n' % (last,
i+1)
#return codestr
exec codestr
result = [] # rescale probability based on frequency
for c in chosen:
result.append(float(c) / n)
return result
def draw4(n, lst, precision = 0.01, maxbins=10000):
""" Memory/speed tradeoff by coarse quantizing of frequency values.
"""
assert len(lst)>1, "Corner case NOT covered"
assert 0.0 < precision < 1.0 and (1.0/precision) < maxbins
binmax = int(1.0/precision)
chosenbin = [0] * binmax
for i in xrange(n):
chosenbin[randrange(binmax)] +=1
left, right = 0, 0 # extract bin range for summation
chosen = [0] * len(lst)
last = 0.0
for i,p in enumerate(lst):
last += p
right = int(last/precision)
chosen[i] = sum( chosenbin[left:right] )
left = right
result = [] # rescale probability based on frequency
for c in chosen:
result.append(float(c) / n)
return result
=== END CODE ===
