<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Observer Bias
The Red Line problem
In Massachusetts, the Red Line is a subway that connects Cambridge and
Boston. When I was working in Cambridge I took the Red Line from Kendall
Square to South Station and caught the commuter rail to Needham. During
rush hour Red Line trains run every 7–8 minutes, on average.
When I arrived at the station, I could estimate the time until the next
train based on the number of passengers on the platform. If there were
only a few people, I inferred that I just missed a train and expected to
wait about 7 minutes. If there were more passengers, I expected the
train to arrive sooner. But if there were a large number of passengers,
I suspected that trains were not running on schedule, so I would go back
to the street level and get a taxi.
While I was waiting for trains, I thought about how Bayesian estimation
could help predict my wait time and decide when I should give up and
take a taxi. This chapter presents the analysis I came up with.
This chapter is based on a project by Brendan Ritter and Kai Austin, who
took a class with me at Olin College. The code in this chapter is
available from http
Step1:  pmf is the actual distribution; new_pmf is the biased
distribution. Inside the loop, we multiply the probability of each
value, x, by the likelihood it will be observed, which is
proportional to x. Then we normalize the result.
Figure [fig.redline0] shows the actual distribution of gaps, labeled
z, and the distribution of gaps seen by passengers, labeled
zb for “z biased”.
Wait times
[fig.redline2]
Wait time, which I call y, is the time between the arrival
of a passenger and the next arrival of a train. Elapsed time, which I
call x, is the time between the arrival of the previous
train and the arrival of a passenger. I chose these definitions so that
zb = x + y.
Given the distribution of zb, we can compute the
distribution of y. I’ll start with a simple case and then
generalize. Suppose, as in the previous example, that zb is
either 5 minutes with probability 1/3, or 10 minutes with probability
2/3.
If we arrive at a random time during a 5 minute gap, y is
uniform from 0 to 5 minutes. If we arrive during a 10 minute gap,
y is uniform from 0 to 10. So the overall distribution is a
mixture of uniform distributions weighted according to the probability
of each gap.
The following function takes the distribution of zb and
computes the distribution of y
Step2:  PmfOfWaitTime makes a meta-Pmf that maps from each uniform
distribution to its probability. Then it uses MakeMixture,
which we saw in Section [mixture], to compute the mixture.
PmfOfWaitTime also uses MakeUniformPmf,
defined here
Step3:  low and high are the range of the uniform
distribution, (both ends included). Finally, MakeUniformPmf
uses MakeRange, defined here
Step4:  MakeRange defines a set of possible values for wait time
(expressed in seconds). By default it divides the range into 10 second
intervals.
To encapsulate the process of computing these distributions, I created a
class called WaitTimeCalculator
Step5:  The parameter, pmf_z, is the unbiased distribution of z.
pmf_zb is the biased distribution of gap time, as seen by passengers.
pmf_y is the distribution of wait time. pmf_x is the distribution of
elapsed time, which is the same as the distribution of wait time. To see
why, remember that for a particular value of zp, the
distribution of y is uniform from 0 to zp.
Also
Step6:  So the distribution of x is also uniform from 0 to
zp.
Figure [fig.redline2] shows the distribution of z,
zb, and y based on the data I collected from
the Red Line web site.
To present these distributions, I am switching from Pmfs to Cdfs. Most
people are more familiar with Pmfs, but I think Cdfs are easier to
interpret, once you get used to them. And if you want to plot several
distributions on the same axes, Cdfs are the way to go.
The mean of z is 7.8 minutes. The mean of zb
is 8.8 minutes, about 13% higher. The mean of y is 4.4,
half the mean of zb.
As an aside, the Red Line schedule reports that trains run every 9
minutes during peak times. This is close to the average of
zb, but higher than the average of z. I
exchanged email with a representative of the MBTA, who confirmed that
the reported time between trains is deliberately conservative in order
to account for variability.
Predicting wait times
[fig.redline3]
Let’s get back to the motivating question
Step7:  pmf_z is the given distribution of gap times.
The next step is to make an ElapsedTimeEstimator (defined
below), which encapsulates the posterior distribution of x
and the predictive distribution of y.
Step8:  The parameters are the WaitTimeCalculator, the passenger
arrival rate, lam (expressed in passengers per second), and
the observed number of passengers, let’s say 15.
Here is the definition of ElapsedTimeEstimator
Step9:  prior_x and posterior_x are the prior and posterior distributions of
elapsed time. pmf_y is the predictive distribution of wait time.
ElapsedTimeEstimator uses Elapsed and
PredictWaitTime, defined below.
Elapsed is a Suite that represents the hypothetical
distribution of x. The prior distribution of x
comes straight from the WaitTimeCalculator. Then we use the
data, which consists of the arrival rate, lam, and the
number of passengers on the platform, to compute the posterior
distribution.
Here’s the definition of Elapsed
Step10:  As always, Likelihood takes a hypothesis and data, and
computes the likelihood of the data under the hypothesis. In this case
hypo is the elapsed time since the last train and
data is a tuple of lam and the number of
passengers.
The likelihood of the data is the probability of getting k
arrivals in x time, given arrival rate lam. We
compute that using the PMF of the Poisson distribution.
Finally, here’s the definition of PredictWaitTime
Step11:  pmf_zb is the distribution of gaps between trains; pmf_x is the
distribution of elapsed time, based on the observed number of
passengers. Since y = zb - x, we can compute
Step12:  The subtraction operator invokes Pmf.__sub__, which enumerates all
pairs of zb and x, computes the differences,
and adds the results to pmf_y.
The resulting Pmf includes some negative values, which we know are
impossible. For example, if you arrive during a gap of 5 minutes, you
can’t wait more than 5 minutes. RemoveNegatives removes the
impossible values from the distribution and renormalizes.
Step13:  Figure [fig.redline3] shows the results. The prior distribution of
x is the same as the distribution of y in
Figure [fig.redline2]. The posterior distribution of x
shows that, after seeing 15 passengers on the platform, we believe that
the time since the last train is probably 5-10 minutes. The predictive
distribution of y indicates that we expect the next train
in less than 5 minutes, with about 80% confidence.
Estimating the arrival rate
[fig.redline1]
The analysis so far has been based on the assumption that we know (1)
the distribution of gaps and (2) the passenger arrival rate. Now we are
ready to relax the second assumption.
Suppose that you just moved to Boston, so you don’t know much about the
passenger arrival rate on the Red Line. After a few days of commuting,
you could make a guess, at least qualitatively. With a little more
effort, you could estimate $\lambda$ quantitatively.
Each day when you arrive at the platform, you should note the time and
the number of passengers waiting (if the platform is too big, you could
choose a sample area). Then you should record your wait time and the
number of new arrivals while you are waiting.
After five days, you might have data like this
Step14:  This Likelihood might look familiar; it is almost identical
to Elapsed.Likelihood in Section [elapsed]. The difference
is that in Elapsed.Likelihood the hypothesis is
x, the elapsed time; in ArrivalRate.Likelihood
the hypothesis is lam, the arrival rate. But in both cases
the likelihood is the probability of seeing k arrivals in
some period of time, given lam.
ArrivalRateEstimator encapsulates the process of estimating
$\lambda$. The parameter, passenger_data, is a list of k1, y,
k2 tuples, as in the table above.
Step15:  __init__ builds hypos, which is a sequence of
hypothetical values for lam, then builds the prior
distribution, prior_lam. The for loop updates the prior
with data, yielding the posterior distribution, post_lam.
Figure [fig.redline1] shows the prior and posterior distributions. As
expected, the mean and median of the posterior are near the observed
rate, 2 passengers per minute. But the spread of the posterior
distribution captures our uncertainty about $\lambda$ based on a small
sample.
Incorporating uncertainty
[fig.redline4]
Whenever there is uncertainty about one of the inputs to an analysis, we
can take it into account by a process like this
Step16:  Next I bias pmf_z to get the distribution of zb, draw a
sample, and then add in delays of 30, 40, and 50 minutes (expressed in
seconds)
Step17:  Cdf.Sample is more efficient than Pmf.Sample,
so it is usually faster to convert a Pmf to a Cdf before sampling.
Next I use the sample of zb to estimate a Pdf using KDE,
and then convert the Pdf to a Pmf
Step18:  Finally I unbias the distribution of zb to get the
distribution of z, which I use to create the
WaitTimeCalculator
Step19:  This process is complicated, but all of the steps are operations we have
seen before. Now we are ready to compute the probability of a long wait.
<ASSISTANT_TASK:>
Python Code:
def BiasPmf(pmf):
    new_pmf = pmf.Copy()
    for x, p in pmf.Items():
        new_pmf.Mult(x, x)
    new_pmf.Normalize()
    return new_pmf
Explanation: Observer Bias
The Red Line problem
In Massachusetts, the Red Line is a subway that connects Cambridge and
Boston. When I was working in Cambridge I took the Red Line from Kendall
Square to South Station and caught the commuter rail to Needham. During
rush hour Red Line trains run every 7–8 minutes, on average.
When I arrived at the station, I could estimate the time until the next
train based on the number of passengers on the platform. If there were
only a few people, I inferred that I just missed a train and expected to
wait about 7 minutes. If there were more passengers, I expected the
train to arrive sooner. But if there were a large number of passengers,
I suspected that trains were not running on schedule, so I would go back
to the street level and get a taxi.
While I was waiting for trains, I thought about how Bayesian estimation
could help predict my wait time and decide when I should give up and
take a taxi. This chapter presents the analysis I came up with.
This chapter is based on a project by Brendan Ritter and Kai Austin, who
took a class with me at Olin College. The code in this chapter is
available from http://thinkbayes.com/redline.py. The code I used to
collect data is in http://thinkbayes.com/redline_data.py. For more
information see Section [download].
The model
[fig.redline0]
Before we get to the analysis, we have to make some modeling decisions.
First, I will treat passenger arrivals as a Poisson process, which means
I assume that passengers are equally likely to arrive at any time, and
that they arrive at an unknown rate, $\lambda$, measured in passengers
per minute. Since I observe passengers during a short period of time,
and at the same time every day, I assume that $\lambda$ is constant.
On the other hand, the arrival process for trains is not Poisson. Trains
to Boston are supposed to leave from the end of the line (Alewife
station) every 7–8 minutes during peak times, but by the time they get
to Kendall Square, the time between trains varies between 3 and 12
minutes.
To gather data on the time between trains, I wrote a script that
downloads real-time data from
http://www.mbta.com/rider_tools/developers/, selects south-bound
trains arriving at Kendall square, and records their arrival times in a
database. I ran the script from 4pm to 6pm every weekday for 5 days, and
recorded about 15 arrivals per day. Then I computed the time between
consecutive arrivals; the distribution of these gaps is shown in
Figure [fig.redline0], labeled z.
If you stood on the platform from 4pm to 6pm and recorded the time
between trains, this is the distribution you would see. But if you
arrive at some random time (without regard to the train schedule) you
would see a different distribution. The average time between trains, as
seen by a random passenger, is substantially higher than the true
average.
Why? Because a passenger is more like to arrive during a large interval
than a small one. Consider a simple example: suppose that the time
between trains is either 5 minutes or 10 minutes with equal probability.
In that case the average time between trains is 7.5 minutes.
But a passenger is more likely to arrive during a 10 minute gap than a 5
minute gap; in fact, twice as likely. If we surveyed arriving
passengers, we would find that 2/3 of them arrived during a 10 minute
gap, and only 1/3 during a 5 minute gap. So the average time between
trains, as seen by an arriving passenger, is 8.33 minutes.
This kind of observer bias appears in many contexts.
Students think that classes are bigger than they are because more of
them are in the big classes. Airline passengers think that planes are
fuller than they are because more of them are on full flights.
In each case, values from the actual distribution are oversampled in
proportion to their value. In the Red Line example, a gap that is twice
as big is twice as likely to be observed.
So given the actual distribution of gaps, we can compute the
distribution of gaps as seen by passengers. BiasPmf does
this computation:
End of explanation
def PmfOfWaitTime(pmf_zb):
    metapmf = thinkbayes.Pmf()
    for gap, prob in pmf_zb.Items():
        uniform = MakeUniformPmf(0, gap)
        metapmf.Set(uniform, prob)
    pmf_y = thinkbayes.MakeMixture(metapmf)
    return pmf_y
Explanation: pmf is the actual distribution; new_pmf is the biased
distribution. Inside the loop, we multiply the probability of each
value, x, by the likelihood it will be observed, which is
proportional to x. Then we normalize the result.
Figure [fig.redline0] shows the actual distribution of gaps, labeled
z, and the distribution of gaps seen by passengers, labeled
zb for “z biased”.
Wait times
[fig.redline2]
Wait time, which I call y, is the time between the arrival
of a passenger and the next arrival of a train. Elapsed time, which I
call x, is the time between the arrival of the previous
train and the arrival of a passenger. I chose these definitions so that
zb = x + y.
Given the distribution of zb, we can compute the
distribution of y. I’ll start with a simple case and then
generalize. Suppose, as in the previous example, that zb is
either 5 minutes with probability 1/3, or 10 minutes with probability
2/3.
If we arrive at a random time during a 5 minute gap, y is
uniform from 0 to 5 minutes. If we arrive during a 10 minute gap,
y is uniform from 0 to 10. So the overall distribution is a
mixture of uniform distributions weighted according to the probability
of each gap.
The following function takes the distribution of zb and
computes the distribution of y:
End of explanation
def MakeUniformPmf(low, high):
    pmf = thinkbayes.Pmf()
    for x in MakeRange(low=low, high=high):
        pmf.Set(x, 1)
    pmf.Normalize()
    return pmf
Explanation: PmfOfWaitTime makes a meta-Pmf that maps from each uniform
distribution to its probability. Then it uses MakeMixture,
which we saw in Section [mixture], to compute the mixture.
PmfOfWaitTime also uses MakeUniformPmf,
defined here:
End of explanation
def MakeRange(low, high, skip=10):
    return range(low, high+skip, skip)
Explanation: low and high are the range of the uniform
distribution, (both ends included). Finally, MakeUniformPmf
uses MakeRange, defined here:
End of explanation
class WaitTimeCalculator(object):
    def __init__(self, pmf_z):
        self.pmf_z = pmf_z
        self.pmf_zb = BiasPmf(pmf)
        self.pmf_y = self.PmfOfWaitTime(self.pmf_zb)
        self.pmf_x = self.pmf_y
Explanation: MakeRange defines a set of possible values for wait time
(expressed in seconds). By default it divides the range into 10 second
intervals.
To encapsulate the process of computing these distributions, I created a
class called WaitTimeCalculator:
End of explanation
x = zp - y
Explanation: The parameter, pmf_z, is the unbiased distribution of z.
pmf_zb is the biased distribution of gap time, as seen by passengers.
pmf_y is the distribution of wait time. pmf_x is the distribution of
elapsed time, which is the same as the distribution of wait time. To see
why, remember that for a particular value of zp, the
distribution of y is uniform from 0 to zp.
Also
End of explanation
wtc = WaitTimeCalculator(pmf_z)
Explanation: So the distribution of x is also uniform from 0 to
zp.
Figure [fig.redline2] shows the distribution of z,
zb, and y based on the data I collected from
the Red Line web site.
To present these distributions, I am switching from Pmfs to Cdfs. Most
people are more familiar with Pmfs, but I think Cdfs are easier to
interpret, once you get used to them. And if you want to plot several
distributions on the same axes, Cdfs are the way to go.
The mean of z is 7.8 minutes. The mean of zb
is 8.8 minutes, about 13% higher. The mean of y is 4.4,
half the mean of zb.
As an aside, the Red Line schedule reports that trains run every 9
minutes during peak times. This is close to the average of
zb, but higher than the average of z. I
exchanged email with a representative of the MBTA, who confirmed that
the reported time between trains is deliberately conservative in order
to account for variability.
Predicting wait times
[fig.redline3]
Let’s get back to the motivating question: suppose that when I arrive at
the platform I see 10 people waiting. How long should I expect to wait
until the next train arrives?
As always, let’s start with the easiest version of the problem and work
our way up. Suppose we are given the actual distribution of
z, and we know that the passenger arrival rate, $\lambda$,
is 2 passengers per minute.
In that case we can:
Use the distribution of z to compute the prior
    distribution of zp, the time between trains as seen by
    a passenger.
Then we can use the number of passengers to estimate the
    distribution of x, the elapsed time since the last
    train.
Finally, we use the relation y = zp - x to get the
    distribution of y.
The first step is to create a WaitTimeCalculator that
encapsulates the distributions of zp, x, and
y, prior to taking into account the number of passengers.
End of explanation
ete = ElapsedTimeEstimator(wtc,
                           lam=2.0/60,
                           num_passengers=15)
Explanation: pmf_z is the given distribution of gap times.
The next step is to make an ElapsedTimeEstimator (defined
below), which encapsulates the posterior distribution of x
and the predictive distribution of y.
End of explanation
class ElapsedTimeEstimator(object):
    def __init__(self, wtc, lam, num_passengers):
        self.prior_x = Elapsed(wtc.pmf_x)
        self.post_x = self.prior_x.Copy()
        self.post_x.Update((lam, num_passengers))
        self.pmf_y = PredictWaitTime(wtc.pmf_zb, self.post_x)
Explanation: The parameters are the WaitTimeCalculator, the passenger
arrival rate, lam (expressed in passengers per second), and
the observed number of passengers, let’s say 15.
Here is the definition of ElapsedTimeEstimator:
End of explanation
class Elapsed(thinkbayes.Suite):
    def Likelihood(self, data, hypo):
        x = hypo
        lam, k = data
        like = thinkbayes.EvalPoissonPmf(k, lam * x)
        return like
Explanation: prior_x and posterior_x are the prior and posterior distributions of
elapsed time. pmf_y is the predictive distribution of wait time.
ElapsedTimeEstimator uses Elapsed and
PredictWaitTime, defined below.
Elapsed is a Suite that represents the hypothetical
distribution of x. The prior distribution of x
comes straight from the WaitTimeCalculator. Then we use the
data, which consists of the arrival rate, lam, and the
number of passengers on the platform, to compute the posterior
distribution.
Here’s the definition of Elapsed:
End of explanation
def PredictWaitTime(pmf_zb, pmf_x):
    pmf_y = pmf_zb - pmf_x
    RemoveNegatives(pmf_y)
    return pmf_y
Explanation: As always, Likelihood takes a hypothesis and data, and
computes the likelihood of the data under the hypothesis. In this case
hypo is the elapsed time since the last train and
data is a tuple of lam and the number of
passengers.
The likelihood of the data is the probability of getting k
arrivals in x time, given arrival rate lam. We
compute that using the PMF of the Poisson distribution.
Finally, here’s the definition of PredictWaitTime:
End of explanation
pmf_y = pmf_zb - pmf_x
Explanation: pmf_zb is the distribution of gaps between trains; pmf_x is the
distribution of elapsed time, based on the observed number of
passengers. Since y = zb - x, we can compute
End of explanation
def RemoveNegatives(pmf):
    for val in pmf.Values():
        if val < 0:
            pmf.Remove(val)
    pmf.Normalize()
Explanation: The subtraction operator invokes Pmf.__sub__, which enumerates all
pairs of zb and x, computes the differences,
and adds the results to pmf_y.
The resulting Pmf includes some negative values, which we know are
impossible. For example, if you arrive during a gap of 5 minutes, you
can’t wait more than 5 minutes. RemoveNegatives removes the
impossible values from the distribution and renormalizes.
End of explanation
class ArrivalRate(thinkbayes.Suite):
    def Likelihood(self, data, hypo):
        lam = hypo
        y, k = data
        like = thinkbayes.EvalPoissonPmf(k, lam * y)
        return like
Explanation: Figure [fig.redline3] shows the results. The prior distribution of
x is the same as the distribution of y in
Figure [fig.redline2]. The posterior distribution of x
shows that, after seeing 15 passengers on the platform, we believe that
the time since the last train is probably 5-10 minutes. The predictive
distribution of y indicates that we expect the next train
in less than 5 minutes, with about 80% confidence.
Estimating the arrival rate
[fig.redline1]
The analysis so far has been based on the assumption that we know (1)
the distribution of gaps and (2) the passenger arrival rate. Now we are
ready to relax the second assumption.
Suppose that you just moved to Boston, so you don’t know much about the
passenger arrival rate on the Red Line. After a few days of commuting,
you could make a guess, at least qualitatively. With a little more
effort, you could estimate $\lambda$ quantitatively.
Each day when you arrive at the platform, you should note the time and
the number of passengers waiting (if the platform is too big, you could
choose a sample area). Then you should record your wait time and the
number of new arrivals while you are waiting.
After five days, you might have data like this:
k1      y     k2
--     ---    --
17     4.6     9
22     1.0     0
23     1.4     4
18     5.4    12
4      5.8    11
where k1 is the number of passengers waiting when you
arrive, y is your wait time in minutes, and k2
is the number of passengers who arrive while you are waiting.
Over the course of one week, you waited 18 minutes and saw 36 passengers
arrive, so you would estimate that the arrival rate is 2 passengers per
minute. For practical purposes that estimate is good enough, but for the
sake of completeness I will compute a posterior distribution for
$\lambda$ and show how to use that distribution in the rest of the
analysis.
ArrivalRate is a Suite that represents
hypotheses about $\lambda$. As always, Likelihood takes a
hypothesis and data, and computes the likelihood of the data under the
hypothesis.
In this case the hypothesis is a value of $\lambda$. The data is a pair,
y, k, where y is a wait time and
k is the number of passengers that arrived.
End of explanation
class ArrivalRateEstimator(object):
    def __init__(self, passenger_data):
        low, high = 0, 5
        n = 51
        hypos = numpy.linspace(low, high, n) / 60
        self.prior_lam = ArrivalRate(hypos)
        self.post_lam = self.prior_lam.Copy()
        for k1, y, k2 in passenger_data:
            self.post_lam.Update((y, k2))
Explanation: This Likelihood might look familiar; it is almost identical
to Elapsed.Likelihood in Section [elapsed]. The difference
is that in Elapsed.Likelihood the hypothesis is
x, the elapsed time; in ArrivalRate.Likelihood
the hypothesis is lam, the arrival rate. But in both cases
the likelihood is the probability of seeing k arrivals in
some period of time, given lam.
ArrivalRateEstimator encapsulates the process of estimating
$\lambda$. The parameter, passenger_data, is a list of k1, y,
k2 tuples, as in the table above.
End of explanation
import thinkbayes
n = 220
cdf_z = thinkbayes.MakeCdfFromList(gap_times)
sample_z = cdf_z.Sample(n)
pmf_z = thinkbayes.MakePmfFromList(sample_z)
Explanation: __init__ builds hypos, which is a sequence of
hypothetical values for lam, then builds the prior
distribution, prior_lam. The for loop updates the prior
with data, yielding the posterior distribution, post_lam.
Figure [fig.redline1] shows the prior and posterior distributions. As
expected, the mean and median of the posterior are near the observed
rate, 2 passengers per minute. But the spread of the posterior
distribution captures our uncertainty about $\lambda$ based on a small
sample.
Incorporating uncertainty
[fig.redline4]
Whenever there is uncertainty about one of the inputs to an analysis, we
can take it into account by a process like this:
Implement the analysis based on a deterministic value of the
    uncertain parameter (in this case $\lambda$).
Compute the distribution of the uncertain parameter.
Run the analysis for each value of the parameter, and generate a set
    of predictive distributions.
Compute a mixture of the predictive distributions, using the weights
    from the distribution of the parameter.
We have already done steps (1) and (2). I wrote a class called
WaitMixtureEstimator to handle steps (3) and (4).
class WaitMixtureEstimator(object):
def __init__(self, wtc, are, num_passengers=15):
    self.metapmf = thinkbayes.Pmf()
    for lam, prob in sorted(are.post_lam.Items()):
        ete = ElapsedTimeEstimator(wtc, lam, num_passengers)
        self.metapmf.Set(ete.pmf_y, prob)
    self.mixture = thinkbayes.MakeMixture(self.metapmf)
wtc is the WaitTimeCalculator that contains
the distribution of zb. are is the
ArrivalTimeEstimator that contains the distribution of
lam.
The first line makes a meta-Pmf that maps from each possible
distribution of y to its probability. For each value of
lam, we use ElapsedTimeEstimator to compute
the corresponding distribution of y and store it in the
Meta-Pmf. Then we use MakeMixture to compute the mixture.
Figure [fig.redline4] shows the results. The shaded lines in the
background are the distributions of y for each value of
lam, with line thickness that represents likelihood. The
dark line is the mixture of these distributions.
In this case we could get a very similar result using a single point
estimate of lam. So it was not necessary, for practical
purposes, to include the uncertainty of the estimate.
In general, it is important to include variability if the system
response is non-linear; that is, if small changes in the input can cause
big changes in the output. In this case, posterior variability in
lam is small and the system response is approximately
linear for small perturbations.
Decision analysis
[fig.redline5]
At this point we can use the number of passengers on the platform to
predict the distribution of wait times. Now let’s get to the second part
of the question: when should I stop waiting for the train and go catch a
taxi?
Remember that in the original scenario, I am trying to get to South
Station to catch the commuter rail. Suppose I leave the office with
enough time that I can wait 15 minutes and still make my connection at
South Station.
In that case I would like to know the probability that y
exceeds 15 minutes as a function of num_passengers. It is easy enough
to use the analysis from Section [elapsed] and run it for a range of
num_passengers.
But there’s a problem. The analysis is sensitive to the frequency of
long delays, and because long delays are rare, it is hard estimate their
frequency.
I only have data from one week, and the longest delay I observed was 15
minutes. So I can’t estimate the frequency of longer delays accurately.
However, I can use previous observations to make at least a coarse
estimate. When I commuted by Red Line for a year, I saw three long
delays caused by a signaling problem, a power outage, and “police
activity” at another stop. So I estimate that there are about 3 major
delays per year.
But remember that my observations are biased. I am more likely to
observe long delays because they affect a large number of passengers. So
we should treat my observations as a sample of zb rather
than z. Here’s how we can do that.
During my year of commuting, I took the Red Line home about 220 times.
So I take the observed gap times, gap_times, generate a sample of 220
gaps, and compute their Pmf:
End of explanation
cdf_zp = BiasPmf(pmf_z).MakeCdf()
sample_zb = cdf_zp.Sample(n) + [1800, 2400, 3000]
Explanation: Next I bias pmf_z to get the distribution of zb, draw a
sample, and then add in delays of 30, 40, and 50 minutes (expressed in
seconds):
End of explanation
pdf_zb = thinkbayes.EstimatedPdf(sample_zb)
xs = MakeRange(low=60)
pmf_zb = pdf_zb.MakePmf(xs)
Explanation: Cdf.Sample is more efficient than Pmf.Sample,
so it is usually faster to convert a Pmf to a Cdf before sampling.
Next I use the sample of zb to estimate a Pdf using KDE,
and then convert the Pdf to a Pmf:
End of explanation
pmf_z = UnbiasPmf(pmf_zb)
wtc = WaitTimeCalculator(pmf_z)
Explanation: Finally I unbias the distribution of zb to get the
distribution of z, which I use to create the
WaitTimeCalculator:
End of explanation
def ProbLongWait(num_passengers, minutes):
    ete = ElapsedTimeEstimator(wtc, lam, num_passengers)
    cdf_y = ete.pmf_y.MakeCdf()
    prob = 1 - cdf_y.Prob(minutes * 60)
Explanation: This process is complicated, but all of the steps are operations we have
seen before. Now we are ready to compute the probability of a long wait.
End of explanation
<END_TASK>