<jupyter_start><jupyter_text>Exercise 1If buses arrive at a bus stop every 20 minutes, and you arrive at the bus stop at a random time, your wait time until the bus arrives is uniformly distributed from 0 to 20 minutes.But in reality, there is variability in the time between buses. Suppose you are waiting for a bus, and you know **the historical distribution of time between buses.**Compute your distribution of wait times.Hint: Suppose that the time between buses is either 5 or 10 minutes with equal probability. What is the probability that you arrive during one of the 10 minute intervals?<jupyter_code># supposed the historical distribution of time between buses is pmf_zb
def PmfOfWaitTime(pmf_zb: Pmf):
    """
    To estimate waiting time distribution based on distribution of time between buses.
    Arguments:
        1, pmf_zb: the historical distribution of time between buses
    Return:
        distribution of waiting time of buses
    """
    metapmf = Pmf()
    for gap, prob in pmf_zb.Items():
        uniform = MakeUniformPmf(0, gap)
        metapmf.Set(uniform, prob)
    pmf_y = MakeMixture(metapmf)
    return pmf_y<jupyter_output><empty_output><jupyter_text>Exercise 2Suppose that passengers arriving at the bus stop are well-modeled by a **Poisson process with parameter λ**. If you arrive at the stop and find 3 people waiting, what is your **posterior distribution for your waiting time since the last bus arrived**? Analysis* Poisson distribution equation    $$ P(x) = { \lambda^k e^{-\lambda} \over k!} $$<jupyter_code>def BiasPmf(pmf):
    """
    Compute Bias pmf based on observation of customer, as customer believe that what they saw ~ likelihood of happening
    """
    new_pmf = pmf.Copy()
    for x, p in pmf.Items():
        # actually doing this: p[x] * likelihood(~x)
        new_pmf.Mult(x, x)
    new_pmf.Normalize()
    return new_pmf
class WaitTimeCalculator(object):
    """
    calculator for waiting time.
    """
    def __init__(self, pmf_z):
        # prior of z for timeslot between trains distribution
        self.pmf_z = pmf_z
        self.pmf_zb = BiasPmf(pmf=pmf_z)
        
        self.pmf_y = self.PmfOfWaitTime(self.pmf_zb)
        self.pmf_x = self.pmf_y
        
    def PmfOfWaitTime(self, pmf_zb: Pmf):
        """
        To estimate waiting time distribution based on distribution of time between buses.
        Arguments:
            1, pmf_zb: the historical distribution of time between buses
        Return:
            distribution of waiting time of buses
        """
        metapmf = Pmf()
        for gap, prob in pmf_zb.Items():
            uniform = MakeUniformPmf(0, gap)
            metapmf.Set(uniform, prob)
        pmf_y = MakeMixture(metapmf)
        return pmf_y
# group by 'date' then calculate the gap between
def readSublineGaps(kendall_redline_path:str = 'redline.csv.gz'):
    kendall_redline_df = pd.read_csv("redline.csv.gz", compression='gzip', parse_dates=['arrival_time', 'departure_time'])
    subline_gaps = []
    for date_ in kendall_redline_df['date'].unique():
        arrival_times = kendall_redline_df[kendall_redline_df['date'] == date_][['arrival_time']]
        for i in range(len(arrival_times) - 1):
            difference = arrival_times.iloc[i+1].arrival_time - arrival_times.iloc[i].arrival_time
            subline_gaps.append(difference.seconds / 60)
    return subline_gaps
# for calculation of true pmf
pmf_z = Pmf()
subline_gaps = readSublineGaps()
for gap in subline_gaps:
    pmf_z.Incr(gap, 1)
pmf_z.Normalize()
pmf_zb = BiasPmf(pmf_z)

# thinkplot for visualization
thinkplot.Clf()
thinkplot.PrePlot(num=2)
thinkplot.Plot(pmf_z, label='real z')
thinkplot.Plot(pmf_zb, label='observed z')
plt.legend()
plt.show()
class EstimatedPdf(Pdf):
    """
    Kernel density estimation from sample
    """
    def __init__(self, sample):
        self.kde = scipy.stats.gaussian_kde(sample)
        
    def Density(self, x):
        return self.kde.evaluate(x)
    
def kde_sample(sample):
    """
    sampling based on KDE
    """
    pdf = EstimatedPdf(sample)
    # reference: internal implementation of pdf.MakePmf
    low, high = min(sample), max(sample)
    n = 101
    xs = np.linspace(low, high, n)
    pmf = pdf.MakePmf(xs=xs)
    return pmf
subline_gaps = readSublineGaps()
pmf_zsmooth = kde_sample(subline_gaps)
pmf_zbsmooth = BiasPmf(pmf_zsmooth)

# thinkplot for visualization
thinkplot.Clf()
thinkplot.PrePlot(num=2)
thinkplot.Plot(pmf_zsmooth, label='real z(kde smooth)')
thinkplot.Plot(pmf_zbsmooth, label='observed z(kde smooth)')
plt.legend()
plt.show()
def RemoveNegatives(pmf):
    leftover = pmf.Copy()
    count = 0
    for val in leftover.Values():
        if val < 0:
            pmf.Remove(val)
            count += 1
    print("original size: ", len(leftover.Values()), "current size: ", len(pmf.Values()))
    pmf.Normalize()
    
def PredictWaitTime(pmf_zb, pmf_x):
    pmf_y = pmf_zb - pmf_x
    RemoveNegatives(pmf_y)
    return pmf_y

class Elapsed(Suite):
    """
    Elapsed time distribution
    """
    def Likelihood(self, data, hypo):
        x = hypo
        lam, k = data
        like = EvalPoissonPmf(lam * x, k)
        return like
class ElapsedTimeEstimator(object):
    def __init__(self, wtc, lam, num_passengers):
        self.prior_x = Elapsed(wtc.pmf_x)
        self.post_x = self.prior_x.Copy()
        # update with (lam, num_passengers) for waiting time prediction
        self.post_x.Update((lam, num_passengers))
        self.pmf_y = PredictWaitTime(wtc.pmf_zb, self.post_x)
wtc = WaitTimeCalculator(pmf_zsmooth)
ete = ElapsedTimeEstimator(wtc, lam=2.0/60, num_passengers=3)

# thinkplot for visualization
thinkplot.Clf()
thinkplot.PrePlot(num=2)
thinkplot.cdf(ete.prior_x.MakeCdf(label='x prior'), color='C1')
thinkplot.cdf(ete.post_x.MakeCdf(label='x post'), color='C2')
thinkplot.cdf(ete.pmf_y.MakeCdf(label='predict y'), color='C3')
plt.legend()
plt.show()<jupyter_output>original size:  71506 current size:  36481<jupyter_text>Exercise 3Suppose that you are an ecologist sampling the insect population in a new environment. You deploy 100 traps in a test area and come back the next day to check on them. You find that 37 traps have been triggered, **trapping an insect inside**. Once a trap triggers, it cannot trap another insect until it has been reset.If you reset the traps and come back in two days, how many traps do you expect to find triggered? Compute a posterior predictive distribution for the number of traps.<jupyter_code># Question: how to define prior of distribution, then update with posterior?
# Prior for traps triggering with uniform probability
class Trap(Suite):
    """
    Trap suite for estimation of triggered trap.
    """
    def __init__(self, label=None, trap_num=100):
        """
        Initialize trap with uniform probablity.
        """
        self.trap_num = trap_num
        pmf = MakeUniformPmf(low=0, high=trap_num)
        super().__init__(pmf, label=label)
        
    def Likelihood(self, data, hypo):
        """
        average of hypo as PoissonPmf for lam.
        """
        lam = hypo
        k = data
        like = EvalPoissonPmf(lam, k)
        return like
def MakeTrapPmf(suite: Trap):
    """
    Makes the distribution of traps scored to restrict posterior with triggered trap <= high.
    """
    mix = Pmf()

    for lam1, prob1 in suite.Items():
        for lam2, prob2 in suite.Items():
            if lam1 + lam2 <= suite.trap_num:
                # impossible larger than suite.trap_num
                mix.Incr(lam1 + lam2, prob1 * prob2)         
    mix.Normalize()     
    return mix
trap_test = Trap(trap_num=100)
trap_test.UpdateSet([37])
trap_in_2day_pmf = MakeTrapPmf(trap_test)
maxProb, maxLam = (0, 0)
for lam, prob in trap_in_2day_pmf.Items():
    if prob > maxProb:
        maxProb = prob
        maxLam = lam
print("MLE trap number: ", maxLam, ", with probablity: ", maxProb)

# thinkplot for visualization
thinkplot.Clf()
thinkplot.PrePlot(num=2)
thinkplot.pmf(trap_in_2day_pmf, label="set traps in 2 days")
plt.legend()
plt.show()<jupyter_output>MLE trap number:  73.0 , with probablity:  0.046400520393770635<jupyter_text>Exercise 4Suppose you are the manager of an apartment building with 100 light bulbs in common areas. It is your responsibility to replace light bulbs when they break.On January 1, all 100 bulbs are working. When you inspect them on February 1, you find 3 light bulbs out. If you come back on April 1, how many light bulbs do you expect to find broken?In the previous exercise, you could reasonably assume that an event is equally likely at any time. For light bulbs, **the likelihood of failure depends on the age of the bulb**. Specifically, old bulbs have an increasing failure rate due to evaporation of the filament.This problem is more open-ended than some; you will have to make modeling decisions. You might want to read about the Weibull distribution (http://en.wikipedia.org/wiki/Weibull_distribution). Or you might want to look around for information about light bulb survival curves.  AnalysisMixture distribution for **Weibull distribution** of natural failure depends on the age of the bulb, not equal probablity. And failure happend in **Poisson distribution** way.Refer to http://blog.csdn.net/mage_ee/article/details/75320838.$$ f(x;\lambda, k) = \begin{cases}{\frac k\lambda}({\frac x\lambda})^{k-1} e^{-({x/\lambda})^{k}},  & \text{$x$ >= 0} \\0, & \text{$x$ <0}\end{cases} $$* **Pre-requisiste**  1. Make $\lambda = 1$? So need to estimate valid $k$ value  2. Simplify equation as follow$$ f(x;1, k) = \begin{cases}kx^{k-1} e^{-x^{k}},  & \text{$x$ >= 0} \\0, & \text{$x$ <0}\end{cases} $$     $\Rightarrow$ $\lambda$ = 1  3. k > 1, then make aging function for failure  $\Rightarrow$ k= 2  * **Sample**  1. $ {\frac k\lambda}({\frac x\lambda})^{k-1} e^{-({x/\lambda})^{k}} $, $ x=1 $, we have $ f(x=1;\lambda, k) = 0.03, x $ is the first month passed  2. Need to estimate Pmf of $ f(x=3;\lambda, k) $  3. High $\approx$ 100/3 = 33.4 without aging function speed-up $\Rightarrow$ $ f(x=33.4;\lambda, k) $ = 1e-5  4. Starting point $\Rightarrow$ $ f(x=0;\lambda, k) $ = 0<jupyter_code># import Weibull distribution library
from thinkbayes2 import EvalWeibullPdf, EvalWeibullCdf, MakeWeibullPmf
# My gut feeling is that we need to use CDF instead of PMF, as bulb broken case will be cumulatived
min_difference = math.inf
candidates = []
high = 33.4
estimate_time = 4
last_time = 33.4
for lam in np.arange(1, 100, 0.5):
    for k in np.arange(1, 100, 0.5):
        pmf = MakeWeibullPmf(lam, k, high)
        cdf = pmf.MakeCdf()
        for x, prob in cdf.Items():
#         for x, prob in pmf.Items():
            if abs(x - 1) < 1e-2:
                # x close to 1
                probOf4 = EvalWeibullCdf(estimate_time, lam, k)
                probOflast = EvalWeibullCdf(last_time, lam, k)
#                 probOf4 = EvalWeibullPdf(estimate_time, lam, k)
#                 probOflast = EvalWeibullPdf(last_time, lam, k)
                # suppose probablity of x=4 is increasing by reality estimate
                if probOf4 > 0.03 and abs(prob - 0.03) < 1e-2 and abs(probOflast - 1.0) < 1e-2:
                    # make probOf4 not close to 1, as only 33.4 will lead to close to 1
                    if abs(probOf4 - 1) > 1e-1:
                        candidates.append((lam, k, probOf4))
assert(len(candidates) > 0)
candidates = sorted(candidates, key=lambda x:x[2])
# select the minimal probablity of bulb
bulb_pmf = MakeWeibullPmf(candidates[0][0], candidates[0][1], high)
probfor4 = candidates[0][2]
print('lambda={}, k={}, we estimate in April 1st, having broken buld = {}'.format( \
    candidates[0][0], candidates[0][1], probfor4 * 100 ))
# thinkplot for visualization
thinkplot.Clf()
thinkplot.PrePlot(num=2)
thinkplot.pmf(bulb_pmf, label="buld life curve")
plt.legend()
plt.xlabel("Elapse month")
plt.ylabel("Probablity of broken")
plt.show()
bulb_cdf=bulb_pmf.MakeCdf()
thinkplot.Clf()
thinkplot.PrePlot(num=2)
thinkplot.cdf(bulb_cdf, label="buld life curve")
plt.legend()
plt.xlabel("Elapse month")
plt.ylabel("Probablity of cumulative broken")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10876876384217865, 'fraction_numerical': 0.019194487736855056, 'mean_word_length': 3.3372465314834576, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '431095', 'license': 'Apache-2.0', 'n_tokens_mistral': 4242, 'n_tokens_neox': 3944, 'n_words': 1261, 'path': 'chapter7/chapter7_exercises.ipynb', 'repo_name': 'llv22/thinkbayes', 'size': '12191'}