Stochastic simulation event timing

Question

So I'm learning bits and pieces of python. Here's some code that just doesn't seem pythonic, and I'd like to improve it. I'm doing stochastic simulations with a Gillespie style approach, but if that doesn't mean anything to you, no worries. I'm trying to avoid some iterations and replace them with something like a list comprehension. The code will work, only I'm looking for a better way to do the same thing.

First I calculate a stopping time (maxTime). Then I calculate the time of an event (and store it if it's less than maxTime). Then the time of the next event (and store again). I repeat until I finally get an event happening after maxTime.

import random
maxTime = random.expovariate(1)

L = []
eventTime = random.expovariate(10)

while eventTime<maxTime:
    L.append(eventTime)
    eventTime += random.expovariate(10)

Any cleaner way to write this code?

Answer 1

If you're writing a simulation, it's probably worthwhile to add some abstractions to your code so that you can free your mind to think at a more abstract level. Ideally, I would like to see

L = list(until_total(stop_time, poisson_process(10)))

(Consider not calling list() if you just need an iterable and not a list.)

Here is one way to get there:

from itertools import takewhile
from random import expovariate

def poisson_process(rate):
    while True:
        yield expovariate(rate)

def until_total(limit, iterable):
    total = [0]                 # http://stackoverflow.com/q/2009402
    def under_total(i):
        total[0] += i
        return total[0] < limit
    return takewhile(under_total, iterable)

stop_time = next(poisson_process(1))
L = until_total(stop_time, poisson_process(10))

Also, consider using more meaningful identifiers:

customer_arrivals = poisson_process(10)
cashier_yawns = poisson_process(1)
customer_interarrival_times = until_total(next(cashier_yawns), customer_arrivals)

Answer 2

This isn't a large change, but you could decouple the creation of the list from the generation of event times. This example creates a simple generator which allows for a list comprehension.

import random

def event_time_gen(inc_lambd, stop_lambd):
    max_time = random.expovariate(stop_lambd)
    event_time = random.expovariate(inc_lambd)
    while event_time < max_time:
        yield event_time
        event_time += random.expovariate(inc_lambd)

L = [rand for rand in event_time_gen(10, 1)]