Poisson Process Simulation

Constructor

First off, dependencies. This class loads numpy, pandas, datetime, and scipy.stats for fairly obvious reasons. Next, our __init__ function with header: def __init__(self, numQueues, rate, start_hour, end_hour, appt_low, appt_high):. This initializes the data storage for the class and needs the basic info to set-up the simulation. Here’s the docstring for parameter descriptions:

"""
Arguments
------
numQueues: integer number of queues to handle arrivals
rate: float, 'scale' parameter of exponential distribution, time between arrivals
start_hour: integer 0:24, hour of day that arrivals can start
end_hour: integer 0:24, hour of day that arrivals stop (no arrivals after, last arrival may be before)
appt_low: integer, lower bound of number of minutes needed to process an arrival
appt_high: integer, upper bound of number of minutes needed to process an arrival

Description
------
This class will simulate a Poisson process where items arrive between 'start_hour' and 'end_hour',
with an item arriving every 'rate' minutes on average. Each arrival requires between 'appt_low'
and 'appt_high' minutes to be processed (uniformly random) and there are 'numQueues' available
to process the arrivals.

Example
-----
'simulation = QueueSimulator(3, 10, 9, 16, 5, 21)'
can be thought of as a clinic with 3 doctors, patients arriving every 10 minutes (exponentially dist.)
between 9AM and 4PM, and each patient requires an appointment between 5 and 20 minutes (uniformly dist.).
The appointment length is a property of each patient, and varies accordingly.

Calling 'simulation.run_simulation(1)' will simulate one day at the clinic. Results stored in the
dataframe 'simulation.results'
"""

To actually accomplish this set-up, there are a few details worth getting into. To work with datetime objects directly, I wanted to minimize work for the user so they just provide start_hour and end_hour as integers, the constructor then runs:

self.start = datetime.datetime.combine(datetime.date.today(), datetime.time(start_hour,0,0))
self.end = datetime.datetime.combine(datetime.date.today(), datetime.time(end_hour,0,0))

This creates a time object from the hours, then combines with today’s date to result in a full datetime object. This was necessary for later use with time arithmetic and timedelta objects which need a full datetime not just a time.

The next issue is that we don’t know how many patients will come each day so we create an expected_count attribute which gives an upper bound of how many patients to simulate arriving. It’s an upper bound because we want to have enough coming in that when we cut-off at the end_hour we haven’t already run out of patients. Here’s the main use of scipy.stats:

minutes_for_new_items = (end_hour-start_hour)*60 #minutes new patients seen
time_between_items = rate #exponential dist. time parameter
self.expected_count = int(np.ceil(stats.poisson.ppf(.9999, minutes_for_new_items/time_between_items)))

Lastly, we need the queues to handle patients (i.e. the doctors) which uses the datetime.datetime.combine tactic from the start and end time, as well as a list comprehension. This results in a list of when each doctor is available next. Initially this is the start_hour and will be updated as we run the simulation.

self.ques = [datetime.datetime.combine(datetime.datetime.today(), datetime.time(start_hour,0,0)) for i in range(0, self.numQueues)]

Everything else is just initializing class attributes.

Run Simulation

The other ‘public’ method is run_simulation which will populate a results attribute dataframe with one row per ‘day’ we simulate, the number of days is the only parameter (default is 1). The results dataframe has the info we need to answer the problems questions. This is primarily a driver function that loops over the number of simulations, resetting the doctors queues, generating a single day of results, aggregating the single run and merging into the results dataframe. The helper functions are: __single_sim_results(self): which calls __wait_time_update(self, item): on each item that arrives.

Single Simulation

The __single_sim_results function first generates all possible patients, arrival delays between patients (exponentially distributed), and then a list of actual arrival times.

itemID = np.arange(0, self.expected_count)
minutes_to_arrival = np.random.exponential(scale = self.rate, size = self.expected_count)
arrival_times = [self.start+datetime.timedelta(minutes = i) for i in minutes_to_arrival.cumsum()]

Next, we construct an arrivals dataframe with one row per patient. This includes generating the appointment length the patient will need, cutting off any who arrive after 4PM, and initializing some other columns:

arrivals = pd.DataFrame({
      'id':itemID,
      'min_btwn_arrival': minutes_to_arrival,
      'arrival_time': arrival_times,
      'appt_length': np.random.uniform(low=self.appt_low, high=self.appt_high, size=self.expected_count)
      })
arrivals = arrivals[arrivals['arrival_time']<=self.end]
arrivals['appt_length_minutes'] = arrivals.appt_length.apply(lambda x: datetime.timedelta(minutes=x))
arrivals['queue'] = np.nan
arrivals['wait_time'] = datetime.timedelta(minutes=0)
arrivals['appt_start_time'] = arrivals['arrival_time']
arrivals['appt_end_time'] = arrivals['arrival_time']+arrivals['appt_length_minutes']

You’ll notice we set the appt_end_time to the arrival time plus the appointment length. Now we have to update this if the patient has to wait.

Wait Time Update

The __wait_time_update routine adjusts the patients appt_end_time and the doctors queue of availability. First, we find the next available doctor and assign the patient to them:

first_que_avail = self.ques.index(min(self.ques))
item['queue'] = first_que_avail

If min(self.ques) (next time the doctor is available) is earlier or at the patients arrival_time, then we can just change the doctors availability to be at the already created appt_end_time since there is no wait. If there is a wait, we must adjust the appt_start_time and appt_end_time and then update the doctors availability tot he appointment end time.

We run all the __wait_time_update via a simple apply call in __single_sim_results:

arrivals = arrivals.apply(lambda x: self.__wait_time_update(x), axis=1)

The important thing about this is the axis=1 argument, which took a little debugging to fix. The __wait_time_update function must be run on each row, and usually the default (axis=0) accomplishes this when we hand a single column to apply (i.e. df.col_name.apply(lambda x: 2*x) would double the values in col_name). The pandas documentation has a nice example, showing that df.apply(func, axis=0) will apply func to each column’s values but df.apply(func, axis=1) applies func to each row’s values. I find this backwards from the usual “axis=1 is columns” standard.