Hands-on: Bayesian credible interval

In [7]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate

We study an exponential decay

$$f(t | \tau) = \frac{1}{\tau} e^{-t/\tau} $$

and observe one event after $t = 1\,$h.

Construct a 68\% Bayesian credible interval for the mean life time $\tau$ using the Jeffreys' prior $\pi(\tau) = 1/\tau$.

In [61]:
# posterior distribution for prior 1/tau
def post(tau):
    # tau in hours
    return 1./tau**2 * np.exp(-1/tau)
In [37]:
# numerical integration of the pdf from 0 to tau_max
def cdf(tau_max):
    result, error = integrate.quad(post, 0, tau_max)
    return result
In [64]:
import scipy.optimize

# return inverse of the cumulative distribution function
def inv_cdf(quantile):
    return scipy.optimize.bisect(lambda tau_max: cdf(tau_max) - quantile, 0.01, 20.)

a) Determine interval limits [$\tau_\mathrm{min}$, $\tau_\mathrm{max}$]

In [ ]:
# you code here

b) Plot the posterior distribution and indicate the credible interval as a shaded area

In [1]:
# plot posterior distribution and indicate credible interval

# your code here
# hint: shaded area under curve: plt.fill_between(x, f(x))
In [ ]: