HiveBrain v1.2.0
Get Started
← Back to all entries
patternpythonMinor

Stochastic simulation event timing

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
timingstochasticsimulationevent

Problem

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?

Solution

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)

Code Snippets

L = list(until_total(stop_time, poisson_process(10)))
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))
customer_arrivals = poisson_process(10)
cashier_yawns = poisson_process(1)
customer_interarrival_times = until_total(next(cashier_yawns), customer_arrivals)

Context

StackExchange Code Review Q#73698, answer score: 2

Revisions (0)

No revisions yet.