Copyright 2017 Allen Downey

MIT License: https://opensource.org/licenses/MIT

```
In [1]:
```from __future__ import print_function, division
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import pandas as pd
import thinkbayes2
import thinkplot

Suppose that you are an auto insurance company interested in the time between collisions for a particular driver. If the probability of a collision is roughly constant over time, the time between collisions will follow an exponential distribution.

Here's an example with parameter $\lambda = 0.5$ collisions / year.

```
In [2]:
```from thinkbayes2 import MakeExponentialPmf
pmf = MakeExponentialPmf(lam=0.5, high=30)
thinkplot.Pdf(pmf)
thinkplot.Config(xlabel='Lifetime', ylabel='PMF')

```
```

For the exponential distribution, the mean and standard deviation are $1/\lambda$.

In this case they are only approximate because we truncated the distribution.

```
In [3]:
```pmf.Mean(), pmf.Std()

```
Out[3]:
```

From the PMF, we can compute the CDF.

```
In [4]:
```cdf = pmf.MakeCdf()
thinkplot.Cdf(cdf)
thinkplot.Config(xlabel='Lifetime', ylabel='CDF')

```
```

And from the CDF, we can compute the survival function, which is the complement of the CDF.

$SF(x) = Prob\{X > x\} = 1 - Prob\{X \le x\} = 1 - CDF(x)$

```
In [5]:
```from survival import MakeSurvivalFromCdf
sf = MakeSurvivalFromCdf(cdf)
thinkplot.Plot(sf)
thinkplot.Config(xlabel='Lifetime', ylabel='Survival function')

```
```

```
In [6]:
```hf = sf.MakeHazardFunction()
thinkplot.Plot(hf)
thinkplot.Config(xlabel='Lifetime', ylabel='Hazard function')

```
```

If the distribution is truly exponential, the hazard function is constant for all $x$.

In this case it goes to 1 at the end, again because we truncated the distribution.

**Exercise:** Go back and increase the value of `high`

, and confirm that the hazard function is a constant until we approach the point where we cut off the distribution.

```
In [7]:
```def RemainingLifetime(sf):
"""Computes remaining lifetime as a function of age.
sf: survival function
returns: Series that maps from age to remaining lifetime
"""
pmf = sf.MakePmf()
d = {}
for t in sorted(pmf.Values()):
pmf[t] = 0
if pmf.Total():
pmf.Normalize()
d[t] = pmf.Mean() - t
return pd.Series(d)

And here's what it looks like for the exponential survival function.

```
In [8]:
```mean_rem_life = RemainingLifetime(sf)
thinkplot.Plot(mean_rem_life)
thinkplot.Config(xlabel='Lifetime', ylabel='Survival function')

```
```

The Weibull distribution is a generalization of the exponential distribution that takes an additional "shape" parameter, `k`

.

When `k=1`

, the Weibull is an exponential distribution. Other values of `k`

yield survival curves with different shapes, and hazard functions that increase, decrease, or both. So the Weibull family can capture a wide range of survival patterns.

```
In [9]:
```from thinkbayes2 import MakeWeibullPmf
pmf = MakeWeibullPmf(lam=2.0, k=1.5, high=30)
thinkplot.Pdf(pmf)
thinkplot.Config(xlabel='Lifetime', ylabel='PMF')

```
```

**Exercise**: In the previous section, replace the exponential distribution with a Weibull distribituion and run the analysis again. What can you infer about the values of the parameters and the behavior of the hazard function and remaining lifetime?

Suppose you are the manager of a large building with many light fixtures. To figure out how often you will need to replace lightbulbs, you install 10 bulbs and measure the time until they fail.

To generate some fake data, I'll choose a Weibull distribution and generate a random sample (let's suppose it's in years):

```
In [10]:
```def SampleWeibull(lam, k, n=1):
return np.random.weibull(k, size=n) * lam
data = SampleWeibull(lam=2, k=1.5, n=10)
data

```
Out[10]:
```

**Exercise:** Write a class called `LightBulb`

that inherits from `Suite`

and provides a `Likelihood`

function that takes an observed lifespan as data and a tuple, `(lam, k)`

, as a hypothesis. It should return a likelihood proportional to the probability of the observed lifespan in a Weibull distribution with the given parameters.

Test your method by creating a `LightBulb`

object with an appropriate prior and update it with the data above.

Plot the posterior distributions of `lam`

and `k`

. As the sample size increases, does the posterior distribution converge on the values of `lam`

and `k`

used to generate the sample?

```
In [11]:
```# Hint
from thinkbayes2 import Suite, Joint, EvalWeibullPdf
class LightBulb(Suite, Joint):
def Likelihood(self, data, hypo):
lam, k = hypo
x = data
like = 1
return like

```
In [12]:
```# Solution
class LightBulb(Suite, Joint):
def Likelihood(self, data, hypo):
lam, k = hypo
x = data
like = EvalWeibullPdf(x, lam, k)
return like

```
In [13]:
```from itertools import product
lams = np.linspace(0.001, 6, 101)
ks = np.linspace(0.001, 8, 101)
suite = LightBulb(product(lams, ks))

```
In [14]:
```suite.UpdateSet(data)

```
Out[14]:
```

```
In [15]:
```thinkplot.Contour(suite)

```
```

```
In [16]:
```pmf_lam = suite.Marginal(0)
thinkplot.Pdf(pmf_lam)
pmf_lam.Mean()

```
Out[16]:
```

```
In [17]:
```pmf_k = suite.Marginal(1)
thinkplot.Pdf(pmf_k)
pmf_k.Mean()

```
Out[17]:
```

**Exercise:** Go back and run this analysis again with `n=20`

and see if the posterior distributions seem to be converging on the actual parameters.

**Exercise:** Now suppose that instead of observing a complete lifespan, you observe a lightbulb that has operated for 1 year and is still working. Write another version of `LightBulb`

that takes data in this form and performs an update.

```
In [18]:
```# Hint
from thinkbayes2 import EvalWeibullCdf
class LightBulb2(Suite, Joint):
def Likelihood(self, data, hypo):
lam, k = hypo
x = data
like = 1
return like

```
In [19]:
```# Solution
class LightBulb2(Suite, Joint):
def Likelihood(self, data, hypo):
lam, k = hypo
x = data
like = 1 - EvalWeibullCdf(x, lam, k)
return like

```
In [20]:
```from itertools import product
lams = np.linspace(0.001, 10, 101)
ks = np.linspace(0.001, 10, 101)
suite = LightBulb2(product(lams, ks))

```
In [21]:
```suite.Update(1)

```
Out[21]:
```

```
In [22]:
```thinkplot.Contour(suite)

```
```

```
In [23]:
```pmf_lam = suite.Marginal(0)
thinkplot.Pdf(pmf_lam)
pmf_lam.Mean()

```
Out[23]:
```

```
In [24]:
```pmf_k = suite.Marginal(1)
thinkplot.Pdf(pmf_k)
pmf_k.Mean()

```
Out[24]:
```

Note: based on this data alone, we can rule out some small values of `lam`

and `k`

, but we can't rule out large values. Without more data or a more informative prior, the results are not useful.

To see why, try increasing the upper bounds in the prior distribition.

```
In [25]:
```# Hint
class LightBulb3(Suite, Joint):
def Likelihood(self, data, hypo):
lam, k = hypo
x = data
like = 1
return like

```
In [26]:
```# Solution
class LightBulb3(Suite, Joint):
def Likelihood(self, data, hypo):
lam, k = hypo
x = data
like = EvalWeibullCdf(x, lam, k)
return like

```
In [27]:
```from itertools import product
lams = np.linspace(0.001, 20, 101)
ks = np.linspace(0.001, 20, 101)
suite = LightBulb3(product(lams, ks))

```
In [28]:
```suite.Update(1)

```
Out[28]:
```

```
In [29]:
```thinkplot.Contour(suite)

```
```

```
In [30]:
```pmf_lam = suite.Marginal(0)
thinkplot.Pdf(pmf_lam)
pmf_lam.Mean()

```
Out[30]:
```

```
In [31]:
```pmf_k = suite.Marginal(1)
thinkplot.Pdf(pmf_k)
pmf_k.Mean()

```
Out[31]:
```

**Exercise:** Suppose you have 15 lightbulbs installed at different times over a 10 year period. When you observe them, some have died and some are still working. Write a version of `LightBulb`

that takes data in the form of a `(flag, x)`

tuple, where:

- If
`flag`

is`eq`

, it means that`x`

is the actual lifespan of a bulb that has died. - If
`flag`

is`gt`

, it means that`x`

is the current age of a bulb that is still working, so it is a lower bound on the lifespan. - If
`flag`

is`lt`

, it means that`x`

is the elapsed time between installation and the first time the bulb is seen broken, so it is an upper bound on the lifespan.

To help you test, I will generate some fake data.

First, I'll generate a Pandas DataFrame with random start times and lifespans. The columns are:

`start`

: time when the bulb was installed`lifespan`

: lifespan of the bulb in years`end`

: time when bulb died or will die`age_t`

: age of the bulb at t=10

```
In [32]:
```import pandas as pd
lam = 2
k = 1.5
n = 15
t_end = 10
starts = np.random.uniform(0, t_end, n)
lifespans = SampleWeibull(lam, k, n)
df = pd.DataFrame({'start': starts, 'lifespan': lifespans})
df['end'] = df.start + df.lifespan
df['age_t'] = t_end - df.start
df.head()

```
Out[32]:
```

Now I'll process the DataFrame to generate data in the form we want for the update.

```
In [33]:
```data = []
for i, row in df.iterrows():
if row.end < t_end:
data.append(('eq', row.lifespan))
else:
data.append(('gt', row.age_t))
for pair in data:
print(pair)

```
```

```
In [34]:
```# Hint
class LightBulb4(Suite, Joint):
def Likelihood(self, data, hypo):
lam, k = hypo
flag, x = data
like = 1
return like

```
In [35]:
```# Solution
class LightBulb4(Suite, Joint):
def Likelihood(self, data, hypo):
lam, k = hypo
flag, x = data
if flag == 'eq':
like = EvalWeibullPdf(x, lam, k)
elif flag == 'gt':
like = 1 - EvalWeibullCdf(x, lam, k)
elif flag == 'lt':
like = EvalWeibullCdf(x, lam, k)
else:
raise ValueError('Invalid data')
return like

```
In [36]:
```from itertools import product
lams = np.linspace(0.001, 10, 101)
ks = np.linspace(0.001, 10, 101)
suite = LightBulb4(product(lams, ks))

```
In [37]:
```suite.UpdateSet(data)

```
Out[37]:
```

```
In [38]:
```thinkplot.Contour(suite)

```
```

```
In [39]:
```pmf_lam = suite.Marginal(0)
thinkplot.Pdf(pmf_lam)
pmf_lam.Mean()

```
Out[39]:
```

```
In [40]:
```pmf_k = suite.Marginal(1)
thinkplot.Pdf(pmf_k)
pmf_k.Mean()

```
Out[40]:
```

Suppose we know that, for a particular kind of lightbulb in a particular location, the distribution of lifespans is well modeled by a Weibull distribution with `lam=2`

and `k=1.5`

. If we install `n=100`

lightbulbs and come back one year later, what is the distribution of `c`

, the number of lightbulbs that have burned out?

The probability that any given bulb has burned out comes from the CDF of the distribution.

```
In [41]:
```lam = 2
k = 1.5
p = EvalWeibullCdf(1, lam, k)
p

```
Out[41]:
```

The number of bulbs that have burned out is distributed Binom(n, p).

```
In [42]:
```from thinkbayes2 import MakeBinomialPmf
n = 100
pmf_c = MakeBinomialPmf(n, p)
thinkplot.Pdf(pmf_c)

```
```

Or we can approximate the distribution with a random sample.

```
In [43]:
```n = 100
sample = np.random.binomial(n, p, 1000)
pdf_c = thinkbayes2.EstimatedPdf(sample)
thinkplot.Pdf(pdf_c)
np.mean(sample), np.std(sample)

```
Out[43]:
```

**Exercise:** Now suppose that `lam`

and `k`

are not known precisely, but we have a `LightBulb`

object that represents the joint posterior distribution of the parameters after seeing some data. Compute the posterior predictive distribution for `c`

, the number of bulbs burned out after one year.

```
In [44]:
```# Solution
n = 100
ns = []
params = suite.Sample(1000)
for lam, k in params:
p = EvalWeibullCdf(1, lam, k)
sample = np.random.binomial(n, p, 10)
ns.extend(sample)
pdf_c_mix = thinkbayes2.EstimatedPdf(ns)
thinkplot.Pdf(pdf_c_mix)
np.mean(ns), np.std(ns)

```
Out[44]:
```

```
In [45]:
```# Solution
# Alternatively, you can get a more precise solution by making
# a mixture of binomials.
from thinkbayes2 import Pmf
n = 100
t_return = 1
metapmf = Pmf()
for (lam, k), prob in suite.Items():
p = EvalWeibullCdf(t_return, lam, k)
pmf = MakeBinomialPmf(n, p)
metapmf[pmf] = prob

```
In [46]:
```# Solution
from thinkbayes2 import MakeMixture
mix = MakeMixture(metapmf)
thinkplot.Pdf(mix)

```
```

```
In [47]:
```# Solution
mix.Mean(), mix.Std()

```
Out[47]:
```

```
In [ ]:
```