<jupyter_start><jupyter_text>The Waiting Time Paradox, or, Why Is My Bus Always Late? *This notebook originally appeared as a [post](http://jakevdp.github.io/blog/2018/09/13/waiting-time-paradox/) on the blog [Pythonic Perambulations](http://jakevdp.github.io).* ![Rapid Ride Bus](http://jakevdp.github.io/images/rapid-ride.jpg)*Image Source: [Wikipedia](https://en.wikipedia.org/wiki/RapidRide/media/File:King_County_Metro_Rapid_Ride_New_Flyer_DE60LFR_6060.JPG) License CC-BY-SA 3.0*If you, like me, frequently commute via public transit, you may be familiar with the following situation:> *You arrive at the bus stop, ready to catch your bus: a line that advertises arrivals every 10 minutes. You glance at your watch and note the time... and when the bus finally comes 11 minutes later, you wonder why you always seem to be so unlucky.*Naïvely, you might expect that if buses are coming every 10 minutes and you arrive at a random time, your average wait would be something like 5 minutes.In reality, though, buses do not arrive exactly on schedule, and so you might wait longer.It turns out that under some reasonable assumptions, you can reach a startling conclusion:**When waiting for a bus that comes on average every 10 minutes, your average waiting time will be 10 minutes.**This is what is sometimes known as the *waiting time paradox*.I've encountered this idea before, and always wondered whether it is actually true... how well do those "reasonable assumptions" match reality?This post will explore the waiting time paradox from the standpoint of both simulation and probabilistic arguments, and then take a look at some real bus arrival time data from the city of Seattle to (hopefully) settle the paradox once and for all.  The Inspection Paradox If buses arrive exactly every ten minutes, it's true that your average wait time will be half that interval: 5 minutes.Qualitatively speaking, it's easy to convince yourself that adding some variation to those arrivals will make the average wait time somewhat longer, as we'll see here.The waiting time paradox turns out to be a particular instance of a more general phenomenon, the *inspection paradox*, which is discussed at length in this enlightening post by <NAME>: [The Inspection Paradox Is Everywhere](http://allendowney.blogspot.com/2015/08/the-inspection-paradox-is-everywhere.html).Briefly, the inspection paradox arises whenever the probability of observing a quantity is related to the quantity being observed.Allen gives one example of surveying university students about the average size of their classes. Though the school may truthfully advertise an average of 30 students per class, the average class size *as experienced by students* can be (and generally will be) much larger. The reason is that there are (of course) more students in the larger classes, and so you oversample large classes when computing the average experience of students.In the case of a nominally 10-minute bus line, sometimes the span between arrivals will be longer than 10 minutes, and sometimes shorter, and if you arrive at a random time, you have more opportunities to encounter a longer interval than to encounter a shorter interval. And so it makes sense that the average span of time *experienced by riders* will be longer than the average span of time between buses, because the longer spans are over-sampled.But the waiting time paradox makes a stronger claim than this: when the average span between arrivals is $N$ minutes, the average span *experienced by riders* is $2N$ minutes.Could this possibly be true?  Simulating Wait Times To convince ourselves that the waiting time paradox is making a reasonable claim, let's start by simulating a stream of buses that arrive at an average of 10 minutes.For the sake of numerical accuracy, we will simulate a large number of bus arrivals: one million buses (or approximately 19 years of round-the-clock 10-minute headways):<jupyter_code>import numpy as np

N = 1000000  # number of buses
tau = 10  # average minutes between arrivals

rand = np.random.RandomState(42)  # universal random seed
bus_arrival_times = N * tau * np.sort(rand.rand(N))<jupyter_output><empty_output><jupyter_text>Just to confirm we've done things correctly, let's check that the mean interval is close to $\tau = 10$:<jupyter_code>intervals = np.diff(bus_arrival_times)
intervals.mean()<jupyter_output><empty_output><jupyter_text>With these bus arrivals simulated, we can now simulate the arrival of a large number of passengers to the bus stop during this span, and compute the wait time that each of them experiences.Let's encapsulate this in a function for later use:<jupyter_code>def simulate_wait_times(arrival_times,
                        rseed=8675309,  # Jenny's random seed
                        n_passengers=1000000):
    rand = np.random.RandomState(rseed)
    
    arrival_times = np.asarray(arrival_times)
    passenger_times = arrival_times.max() * rand.rand(n_passengers)

    # find the index of the next bus for each simulated passenger
    i = np.searchsorted(arrival_times, passenger_times, side='right')

    return arrival_times[i] - passenger_times<jupyter_output><empty_output><jupyter_text>We can then simulate some wait times and compute the average:<jupyter_code>wait_times = simulate_wait_times(bus_arrival_times)
wait_times.mean()<jupyter_output><empty_output><jupyter_text>The average wait time is also close to 10 minutes, just as the waiting time paradox predicted.  Digging Deeper: Probabilities & Poisson Processes How can we understand what's going on here?Fundamentally, this is an instance of the inspection paradox, in which the probability of observing a value is related to the value itself.Let's denote by $p(T)$ the distribution of intervals $T$ between buses as they arrive at a bus stop. In this notation, the expectation value of the arrival times is$$E[T] = \int_0^\infty T~p(T)~dT$$In the above simulation, we had chosen $E[T] = \tau = 10$ minutes. When a rider arrives at a bus stop at a random time, the probability of the time interval they experience will be affected by $p(T)$, but also by $T$ itself: the longer the interval, the larger the probability is that a passenger will experience it.So we can write the distribution of arrival times experienced by passengers:$$p_{exp}(T) \propto T~p(T)$$The constant of proportionality comes from normalizing the distribution:$$p_{exp}(T) = \frac{T~p(T)}{\int_0^\infty T~p(T)~dT}$$Comparing to above we see this simplifies to$$p_{exp}(T) = \frac{T~p(T)}{E[T]} $$The expected wait time $E[W]$ will then be half of the expected interval experienced by passengers, so we can write$$E[W] = \frac{1}{2}E_{exp}[T] = \frac{1}{2}\int_0^\infty T~p_{exp}(T)~dT$$which can be rewritten in a more suggestive way:$$E[W] = \frac{E[T^2]}{2E[T]}$$and now all that remains is for us to choose a form for $p(T)$ and compute the integrals.  Choosing  *p(T)*With this formalism worked out, what is a reasonable distribution to use for $p(T)$?We can get a picture of the $p(T)$ distribution within our simulated arrivals by plotting a histogram of the intervals between arrivals:<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn')

plt.hist(intervals, bins=np.arange(80), density=True)
plt.axvline(intervals.mean(), color='black', linestyle='dotted')
plt.xlabel('Interval between arrivals (minutes)')
plt.ylabel('Probability density');<jupyter_output><empty_output><jupyter_text>The vertical dotted line here shows the mean interval of about 10 minutes.This looks very much like an exponential distribution, and that is no accident: our simulation of bus arrival times as uniform random numbers very closely approximates a [Poisson process](https://en.wikipedia.org/wiki/Poisson_point_process), and for such a process it can be shown that the distribution of intervals between arrivals is exponential.(Side note: In our case this is only approximately exponential; in reality the intervals $T$ between $N$ uniformly sampled points within a timespan $N\tau$ follow the [Beta distribution](https://en.wikipedia.org/wiki/Beta_distribution): $T/(N\tau) \sim \mathrm{Beta}[1, N]$, which in the large $N$ limit approaches $T \sim \mathrm{Exp}[1/\tau]$.See, e.g. [this StackExchange post](https://math.stackexchange.com/questions/68749/difference-of-order-statistics-in-a-sample-of-uniform-random-variables), or  [this twitter thread](https://twitter.com/jakevdp/status/1038086415190351872) for more details).An exponential distribution of intervals implies that the arrival times follow a Poisson process.To double-check this reasoning, we can confirm that it matches another property of a Poisson process: the number of arrivals within a fixed span of time will be Poisson-distributed.Let's check this by binning our simulated arrivals into hourly blocks:<jupyter_code>from scipy.stats import poisson

# count the number of arrivals in 1-hour bins
binsize = 60
binned_arrivals = np.bincount((bus_arrival_times // binsize).astype(int))
x = np.arange(20)

# plot the results
plt.hist(binned_arrivals, bins=x - 0.5, density=True, alpha=0.5, label='simulation')
plt.plot(x, poisson(binsize / tau).pmf(x), 'ok', label='Poisson prediction')
plt.xlabel('Number of arrivals per hour')
plt.ylabel('frequency')
plt.legend();<jupyter_output><empty_output><jupyter_text>The close match between the empirical and theoretical values gives confidence that our interpretation is correct: for large $N$, the arrival times we simulated above are well-described by a Poisson process, which implies exponentially-distributed arrival intervals. That means we can write the probability distribution:$$p(T) = \frac{1}{\tau}e^{-T/\tau}$$Plugging this into the above results, we find that the average waiting time experienced by a person is$$E[W] = \frac{\int_0^\infty T^2~e^{-T/\tau}}{2\int_0^\infty T~e^{-T/\tau}} = \frac{2\tau^3}{2(\tau^2)} = \tau$$For bus arrivals consistent with a Poisson process, the expected wait time for a passenger is identical to the average interval between arrivals. A complementary way to reason about this is this: a Poisson process is a *memoryless* process, meaning that the history of events has no bearing on the expected time to the next event.So when you arrive at the bus stop, the average waiting time until the next bus is always the same: in our case, it is 10 minutes, and this is regardless of how long it has been since the previous bus!Along the same lines, it does not matter how long you have been waiting already: the expected time to the next arrival is always exactly 10 minutes: for a Poisson process, you get no "credit" for time spent waiting.  Waiting Times In RealityThe above is well and good if real-world bus arrivals are actually described by a Poisson process, but are they?[![Seattle Transit Map](http://jakevdp.github.io/images/seattle-transit-map.png)](https://seattletransitmap.com/)*Image Source: https://seattletransitmap.com/*To determine whether the waiting time paradox describes reality, we can dig into some data, available for download here: [arrival_times.csv](https://gist.githubusercontent.com/jakevdp/82409002fcc5142a2add0168c274a869/raw/1bbabf78333306dbc45b9f33662500957b2b6dc3/arrival_times.csv) (3MB CSV file).The dataset contains scheduled and actual arrival times for Seattle's [Rapid Ride](https://kingcounty.gov/depts/transportation/metro/travel-options/bus/rapidride.aspx) lines C, D, and E at the 3rd & Pike bus stop in downtown Seattle, recorded during the second quarter of 2016 (huge thanks to <NAME> of the Washington State Transportation Center for providing this data!).<jupyter_code>import pandas as pd
df = pd.read_csv('arrival_times.csv')
df = df.dropna(axis=0, how='any')
df.head()<jupyter_output><empty_output><jupyter_text>The reason I sought data from Rapid Ride routes in particular is that for much of the day, the buses are scheduled at regular intervals of between 10 and 15 minutes — not to mention the fact that I'm a frequent passenger on the C-line.  Data CleanupTo start with, let's do a little bit of data cleanup to get it into a form that's easier to work with:<jupyter_code># combine date and time into a single timestamp
df['scheduled'] = pd.to_datetime(df['OPD_DATE'] + ' ' + df['SCH_STOP_TM'])
df['actual'] = pd.to_datetime(df['OPD_DATE'] + ' ' + df['ACT_STOP_TM'])

# if scheduled & actual span midnight, then the actual day needs to be adjusted
minute = np.timedelta64(1, 'm')
hour = 60 * minute
diff_hrs = (df['actual'] - df['scheduled']) / hour
df.loc[diff_hrs > 20, 'actual'] -= 24 * hour
df.loc[diff_hrs < -20, 'actual'] += 24 * hour
df['minutes_late'] = (df['actual'] - df['scheduled']) / minute

# map internal route codes to external route letters
df['route'] = df['RTE'].replace({673: 'C', 674: 'D', 675: 'E'}).astype('category')
df['direction'] = df['DIR'].replace({'N': 'northbound', 'S': 'southbound'}).astype('category')

# extract useful columns
df = df[['route', 'direction', 'scheduled', 'actual', 'minutes_late']].copy()

df.head()<jupyter_output><empty_output><jupyter_text>How Late Are Buses? There are essentially six different datasets within this table: the northbound and southbound directions for each of the C, D, and E lines.To get a feeling for their characteristics, let's plot a histogram of the actual minus scheduled arrival times for each of these six:<jupyter_code>import seaborn as sns
g = sns.FacetGrid(df, row="direction", col="route")
g.map(plt.hist, "minutes_late", bins=np.arange(-10, 20))
g.set_titles('{col_name} {row_name}')
g.set_axis_labels('minutes late', 'number of buses');<jupyter_output><empty_output><jupyter_text>You might expect that the buses stick closer to their schedule near the beginning of each one-way trip and show more spread near the end, and this is borne out in the data: the southbound C-line and northbound D and E lines are near the beginning of their respective routes, and in the opposite direction they are near the end.  Scheduled and Observed Arrival Intervals Let's next take a look at the observed and scheduled intervals between arrivals for these six routes.We'll start by using Pandas `groupby` functionality to compute these intervals:<jupyter_code>def compute_headway(scheduled):
    minute = np.timedelta64(1, 'm')
    return scheduled.sort_values().diff() / minute

grouped = df.groupby(['route', 'direction'])
df['actual_interval'] = grouped['actual'].transform(compute_headway)
df['scheduled_interval'] = grouped['scheduled'].transform(compute_headway)
g = sns.FacetGrid(df.dropna(), row="direction", col="route")
g.map(plt.hist, "actual_interval", bins=np.arange(50) + 0.5)
g.set_titles('{col_name} {row_name}')
g.set_axis_labels('actual interval (minutes)', 'number of buses');<jupyter_output><empty_output><jupyter_text>It's already clear that these don't look much like the exponential distribution of our model, but that is not telling us much yet: the distributions may be affected by non-constant scheduled arrival intervals.Let's repeat the above chart, examining the scheduled rather than observed arrival intervals:<jupyter_code>g = sns.FacetGrid(df.dropna(), row="direction", col="route")
g.map(plt.hist, "scheduled_interval", bins=np.arange(20) - 0.5)
g.set_titles('{col_name} {row_name}')
g.set_axis_labels('scheduled interval (minutes)', 'frequency');<jupyter_output><empty_output><jupyter_text>This shows that the buses come at variety of arrival intervals thorughout the week, so we cannot evaluate the accuracy of the waiting time paradox from the distributions of raw arrival times.  Constructing Uniform Schedules Even though the scheduled arrival intervals are not uniform, there are a few particular intervals that have a large number of arrivals: for example, there are nearly 2000 northbound E-line buses with a scheduled interval of 10 minutes.In order to explore whether the waiting time paradox applies, let's group the data by line, direction, and scheduled interval, then re-stack these similar arrivals together as if they had happened in sequence.This should maintain all the relevant characteristics of the raw data, while making it easier to directly compare with the predictions of the waiting time paradox.<jupyter_code>def stack_sequence(data):
    # first, sort by scheduled time
    data = data.sort_values('scheduled')
    
    # re-stack data & recompute relevant quantities
    data['scheduled'] = data['scheduled_interval'].cumsum()
    data['actual'] = data['scheduled'] + data['minutes_late']
    data['actual_interval'] = data['actual'].sort_values().diff()
    return data

subset = df[df.scheduled_interval.isin([10, 12, 15])]
grouped = subset.groupby(['route', 'direction', 'scheduled_interval'])
sequenced = grouped.apply(stack_sequence).reset_index(drop=True)
sequenced.head()<jupyter_output><empty_output><jupyter_text>Using this cleaned data, we can plot the distribution of "actual" arrival intervals for each route, direction, and arrival frequency:<jupyter_code>for route in ['C', 'D', 'E']:
    g = sns.FacetGrid(sequenced.query(f"route == '{route}'"),
                      row="direction", col="scheduled_interval")
    g.map(plt.hist, "actual_interval", bins=np.arange(40) + 0.5)
    g.set_titles('{row_name} ({col_name:.0f} min)')
    g.set_axis_labels('actual interval (min)', 'count')
    g.fig.set_size_inches(8, 4)
    g.fig.suptitle(f'{route} line', y=1.05, fontsize=14)<jupyter_output><empty_output><jupyter_text>We see that for each line and schedule, the distribution of observed arrival intervals is nearly Gaussian, is peaked near the scheduled arrival interval, and has a standard deviation that is smaller near the beginning of the route (southbound for C, northbound for D/E) and larger near the end.Even without a statistical test, it's clear by eye that the actual arrival intervals are definitely **not** exponentially distributed, which is the basic assumption on which the waiting time paradox rests.We can make use of the wait time simulation function we used above in order to find the average wait time for each bus line, direction, and schedule:<jupyter_code>grouped = sequenced.groupby(['route', 'direction', 'scheduled_interval'])
sims = grouped['actual'].apply(simulate_wait_times)
sims.apply(lambda times: "{0:.1f} +/- {1:.1f}".format(times.mean(), times.std()))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09199411412066053, 'fraction_numerical': 0.014387705052046434, 'mean_word_length': 5.436338126972992, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 61, 'https://': 9, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '477452', 'license': 'MIT', 'n_tokens_mistral': 5423, 'n_tokens_neox': 5025, 'n_words': 2192, 'path': 'content/downloads/notebooks/WaitingTimeParadox.ipynb', 'repo_name': 'Jerska/jakevdp.github.io-source', 'size': '18364'}