{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Survival Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Think Bayes, Second Edition\n",
    "\n",
    "Copyright 2020 Allen B. Downey\n",
    "\n",
    "License: [Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:34.889427Z",
     "iopub.status.busy": "2021-04-16T19:37:34.888019Z",
     "iopub.status.idle": "2021-04-16T19:37:34.891992Z",
     "shell.execute_reply": "2021-04-16T19:37:34.891386Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# If we're running on Colab, install empiricaldist\n",
    "# https://pypi.org/project/empiricaldist/\n",
    "\n",
    "import sys\n",
    "IN_COLAB = 'google.colab' in sys.modules\n",
    "\n",
    "if IN_COLAB:\n",
    "    !pip install empiricaldist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:34.895783Z",
     "iopub.status.busy": "2021-04-16T19:37:34.895192Z",
     "iopub.status.idle": "2021-04-16T19:37:34.896815Z",
     "shell.execute_reply": "2021-04-16T19:37:34.897202Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Get utils.py\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://github.com/AllenDowney/ThinkBayes2/raw/master/soln/utils.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:34.900465Z",
     "iopub.status.busy": "2021-04-16T19:37:34.899738Z",
     "iopub.status.idle": "2021-04-16T19:37:35.562654Z",
     "shell.execute_reply": "2021-04-16T19:37:35.562138Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from utils import set_pyplot_params\n",
    "set_pyplot_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This chapter introduces \"survival analysis\", which is a set of statistical methods used to answer questions about the time until an event.\n",
    "In the context of medicine it is literally about survival, but it can be applied to the time until any kind of event, or instead of time it can be about space or other dimensions.\n",
    "\n",
    "Survival analysis is challenging because the data we have are often incomplete.  But as we'll see, Bayesian methods are particularly good at working with incomplete data.\n",
    "\n",
    "As examples, we'll consider two applications that are a little less serious than life and death: the time until light bulbs fail and the time until dogs in a shelter are adopted.\n",
    "To describe these \"survival times\", we'll use the Weibull distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Weibull Distribution\n",
    "\n",
    "The [Weibull distribution](https://en.wikipedia.org/wiki/Weibull_distribution) is often used in survival analysis because it is a good model for the distribution of lifetimes for manufactured products, at least over some parts of the range.\n",
    "\n",
    "SciPy provides several versions of the Weibull distribution; the one we'll use is called `weibull_min`.\n",
    "To make the interface consistent with our notation, I'll wrap it in a function that takes as parameters $\\lambda$, which mostly affects the location or \"central tendency\" of the distribution, and $k$, which affects the shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.566415Z",
     "iopub.status.busy": "2021-04-16T19:37:35.565864Z",
     "iopub.status.idle": "2021-04-16T19:37:35.569596Z",
     "shell.execute_reply": "2021-04-16T19:37:35.569017Z"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.stats import weibull_min\n",
    "\n",
    "def weibull_dist(lam, k):\n",
    "    return weibull_min(k, scale=lam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, here's a Weibull distribution with parameters $\\lambda=3$ and $k=0.8$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.575101Z",
     "iopub.status.busy": "2021-04-16T19:37:35.574198Z",
     "iopub.status.idle": "2021-04-16T19:37:35.576663Z",
     "shell.execute_reply": "2021-04-16T19:37:35.577224Z"
    }
   },
   "outputs": [],
   "source": [
    "lam = 3\n",
    "k = 0.8\n",
    "actual_dist = weibull_dist(lam, k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is an object that represents the distribution.\n",
    "Here's what the Weibull CDF looks like with those parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.584820Z",
     "iopub.status.busy": "2021-04-16T19:37:35.583794Z",
     "iopub.status.idle": "2021-04-16T19:37:35.757065Z",
     "shell.execute_reply": "2021-04-16T19:37:35.756515Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from empiricaldist import Cdf\n",
    "from utils import decorate\n",
    "\n",
    "qs = np.linspace(0, 12, 101)\n",
    "ps = actual_dist.cdf(qs)\n",
    "cdf = Cdf(ps, qs)\n",
    "cdf.plot()\n",
    "\n",
    "decorate(xlabel='Duration in time', \n",
    "         ylabel='CDF',\n",
    "         title='CDF of a Weibull distribution')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`actual_dist` provides `rvs`, which we can use to generate a random sample from this distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.760484Z",
     "iopub.status.busy": "2021-04-16T19:37:35.760006Z",
     "iopub.status.idle": "2021-04-16T19:37:35.762412Z",
     "shell.execute_reply": "2021-04-16T19:37:35.761932Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "np.random.seed(17)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.766596Z",
     "iopub.status.busy": "2021-04-16T19:37:35.765998Z",
     "iopub.status.idle": "2021-04-16T19:37:35.768601Z",
     "shell.execute_reply": "2021-04-16T19:37:35.768230Z"
    }
   },
   "outputs": [],
   "source": [
    "data = actual_dist.rvs(10)\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, given the parameters of the distribution, we can generate a sample.\n",
    "Now let's see if we can go the other way: given the sample, we'll estimate the parameters.\n",
    "\n",
    "Here's a uniform prior distribution for $\\lambda$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.773021Z",
     "iopub.status.busy": "2021-04-16T19:37:35.772572Z",
     "iopub.status.idle": "2021-04-16T19:37:35.778527Z",
     "shell.execute_reply": "2021-04-16T19:37:35.778919Z"
    }
   },
   "outputs": [],
   "source": [
    "from utils import make_uniform\n",
    "\n",
    "lams = np.linspace(0.1, 10.1, num=101)\n",
    "prior_lam = make_uniform(lams, name='lambda')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a uniform prior for $k$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.783338Z",
     "iopub.status.busy": "2021-04-16T19:37:35.782851Z",
     "iopub.status.idle": "2021-04-16T19:37:35.784420Z",
     "shell.execute_reply": "2021-04-16T19:37:35.784771Z"
    }
   },
   "outputs": [],
   "source": [
    "ks = np.linspace(0.1, 5.1, num=101)\n",
    "prior_k = make_uniform(ks, name='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll use `make_joint` to make a joint prior distribution for the two parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.788549Z",
     "iopub.status.busy": "2021-04-16T19:37:35.787996Z",
     "iopub.status.idle": "2021-04-16T19:37:35.789987Z",
     "shell.execute_reply": "2021-04-16T19:37:35.790430Z"
    }
   },
   "outputs": [],
   "source": [
    "from utils import make_joint\n",
    "\n",
    "prior = make_joint(prior_lam, prior_k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a `DataFrame` that represents the joint prior, with possible values of $\\lambda$ across the columns and values of $k$ down the rows.\n",
    "\n",
    "Now I'll use `meshgrid` to make a 3-D mesh with $\\lambda$ on the first axis (`axis=0`), $k$ on the second axis (`axis=1`), and the data on the third axis (`axis=2`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.794031Z",
     "iopub.status.busy": "2021-04-16T19:37:35.793470Z",
     "iopub.status.idle": "2021-04-16T19:37:35.798668Z",
     "shell.execute_reply": "2021-04-16T19:37:35.798094Z"
    }
   },
   "outputs": [],
   "source": [
    "lam_mesh, k_mesh, data_mesh = np.meshgrid(\n",
    "    prior.columns, prior.index, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use `weibull_dist` to compute the PDF of the Weibull distribution for each pair of parameters and each data point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.803606Z",
     "iopub.status.busy": "2021-04-16T19:37:35.802941Z",
     "iopub.status.idle": "2021-04-16T19:37:35.829401Z",
     "shell.execute_reply": "2021-04-16T19:37:35.828903Z"
    }
   },
   "outputs": [],
   "source": [
    "densities = weibull_dist(lam_mesh, k_mesh).pdf(data_mesh)\n",
    "densities.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The likelihood of the data is the product of the probability densities along `axis=2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.832478Z",
     "iopub.status.busy": "2021-04-16T19:37:35.831851Z",
     "iopub.status.idle": "2021-04-16T19:37:35.835448Z",
     "shell.execute_reply": "2021-04-16T19:37:35.835040Z"
    }
   },
   "outputs": [],
   "source": [
    "likelihood = densities.prod(axis=2)\n",
    "likelihood.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compute the posterior distribution in the usual way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.839158Z",
     "iopub.status.busy": "2021-04-16T19:37:35.838606Z",
     "iopub.status.idle": "2021-04-16T19:37:35.843194Z",
     "shell.execute_reply": "2021-04-16T19:37:35.842756Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from utils import normalize\n",
    "\n",
    "posterior = prior * likelihood\n",
    "normalize(posterior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function encapsulates these steps.\n",
    "It takes a joint prior distribution and the data, and returns a joint posterior distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.847328Z",
     "iopub.status.busy": "2021-04-16T19:37:35.846908Z",
     "iopub.status.idle": "2021-04-16T19:37:35.849172Z",
     "shell.execute_reply": "2021-04-16T19:37:35.848696Z"
    }
   },
   "outputs": [],
   "source": [
    "def update_weibull(prior, data):\n",
    "    \"\"\"Update the prior based on data.\"\"\"\n",
    "    lam_mesh, k_mesh, data_mesh = np.meshgrid(\n",
    "        prior.columns, prior.index, data)\n",
    "    \n",
    "    densities = weibull_dist(lam_mesh, k_mesh).pdf(data_mesh)\n",
    "    likelihood = densities.prod(axis=2)\n",
    "\n",
    "    posterior = prior * likelihood\n",
    "    normalize(posterior)\n",
    "\n",
    "    return posterior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we use it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.852532Z",
     "iopub.status.busy": "2021-04-16T19:37:35.852026Z",
     "iopub.status.idle": "2021-04-16T19:37:35.875653Z",
     "shell.execute_reply": "2021-04-16T19:37:35.875092Z"
    }
   },
   "outputs": [],
   "source": [
    "posterior = update_weibull(prior, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's a contour plot of the joint posterior distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:35.889429Z",
     "iopub.status.busy": "2021-04-16T19:37:35.888348Z",
     "iopub.status.idle": "2021-04-16T19:37:36.051329Z",
     "shell.execute_reply": "2021-04-16T19:37:36.050848Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from utils import plot_contour\n",
    "\n",
    "plot_contour(posterior)\n",
    "decorate(title='Posterior joint distribution of Weibull parameters')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like the range of likely values for $\\lambda$ is about 1 to 4, which contains the actual value we used to generate the data, 3.\n",
    "And the range for $k$ is about 0.5 to 1.5, which contains the actual value, 0.8."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Marginal Distributions\n",
    "\n",
    "To be more precise about these ranges, we can extract the marginal distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.055939Z",
     "iopub.status.busy": "2021-04-16T19:37:36.055450Z",
     "iopub.status.idle": "2021-04-16T19:37:36.057535Z",
     "shell.execute_reply": "2021-04-16T19:37:36.057916Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from utils import marginal\n",
    "\n",
    "posterior_lam = marginal(posterior, 0)\n",
    "posterior_k = marginal(posterior, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "And compute the posterior means and 90% credible intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.083939Z",
     "iopub.status.busy": "2021-04-16T19:37:36.077060Z",
     "iopub.status.idle": "2021-04-16T19:37:36.242353Z",
     "shell.execute_reply": "2021-04-16T19:37:36.241898Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.axvline(3, color='C5')\n",
    "posterior_lam.plot(color='C4', label='lambda')\n",
    "decorate(xlabel='lam',\n",
    "         ylabel='PDF', \n",
    "         title='Posterior marginal distribution of lam')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The vertical gray line show the actual value of $\\lambda$.\n",
    "\n",
    "Here's the marginal posterior distribution for $k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.258259Z",
     "iopub.status.busy": "2021-04-16T19:37:36.255943Z",
     "iopub.status.idle": "2021-04-16T19:37:36.436794Z",
     "shell.execute_reply": "2021-04-16T19:37:36.436423Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.axvline(0.8, color='C5')\n",
    "posterior_k.plot(color='C12', label='k')\n",
    "decorate(xlabel='k',\n",
    "         ylabel='PDF', \n",
    "         title='Posterior marginal distribution of k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The posterior distributions are wide, which means that with only 10 data points we can't estimated the parameters precisely.\n",
    "But for both parameters, the actual value falls in the credible interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.440827Z",
     "iopub.status.busy": "2021-04-16T19:37:36.440375Z",
     "iopub.status.idle": "2021-04-16T19:37:36.443293Z",
     "shell.execute_reply": "2021-04-16T19:37:36.442768Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(lam, posterior_lam.credible_interval(0.9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.448649Z",
     "iopub.status.busy": "2021-04-16T19:37:36.448148Z",
     "iopub.status.idle": "2021-04-16T19:37:36.450658Z",
     "shell.execute_reply": "2021-04-16T19:37:36.450297Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(k, posterior_k.credible_interval(0.9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Incomplete Data\n",
    "\n",
    "In the previous example we were given 10 random values from a Weibull distribution, and we used them to estimate the parameters (which we pretended we didn't know).\n",
    "\n",
    "But in many real-world scenarios, we don't have complete data; in particular, when we observe a system at a point in time, we generally have information about the past, but not the future.\n",
    "\n",
    "As an example, suppose you work at a dog shelter and you are interested in the time between the arrival of a new dog and when it is adopted.\n",
    "Some dogs might be snapped up immediately; others might have to wait longer.\n",
    "The people who operate the shelter might want to make inferences about the distribution of these residence times.\n",
    "\n",
    "Suppose you monitor arrivals and departures over 8 weeks and 10 dogs arrive during that interval.\n",
    "I'll assume that their arrival times are distributed uniformly, so I'll generate random values like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.454085Z",
     "iopub.status.busy": "2021-04-16T19:37:36.453519Z",
     "iopub.status.idle": "2021-04-16T19:37:36.455347Z",
     "shell.execute_reply": "2021-04-16T19:37:36.455719Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "np.random.seed(19)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.459627Z",
     "iopub.status.busy": "2021-04-16T19:37:36.459000Z",
     "iopub.status.idle": "2021-04-16T19:37:36.461880Z",
     "shell.execute_reply": "2021-04-16T19:37:36.462341Z"
    }
   },
   "outputs": [],
   "source": [
    "start = np.random.uniform(0, 8, size=10)\n",
    "start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's suppose that the residence times follow the Weibull distribution we used in the previous example.\n",
    "We can generate a sample from that distribution like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.466403Z",
     "iopub.status.busy": "2021-04-16T19:37:36.465743Z",
     "iopub.status.idle": "2021-04-16T19:37:36.467654Z",
     "shell.execute_reply": "2021-04-16T19:37:36.468115Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "np.random.seed(17)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.473358Z",
     "iopub.status.busy": "2021-04-16T19:37:36.472568Z",
     "iopub.status.idle": "2021-04-16T19:37:36.476129Z",
     "shell.execute_reply": "2021-04-16T19:37:36.475500Z"
    }
   },
   "outputs": [],
   "source": [
    "duration = actual_dist.rvs(10)\n",
    "duration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll use these values to construct a `DataFrame` that contains the arrival and departure times for each dog, called `start` and `end`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.481746Z",
     "iopub.status.busy": "2021-04-16T19:37:36.480901Z",
     "iopub.status.idle": "2021-04-16T19:37:36.483162Z",
     "shell.execute_reply": "2021-04-16T19:37:36.483712Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "d = dict(start=start, end=start+duration)\n",
    "obs = pd.DataFrame(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For display purposes, I'll sort the rows of the `DataFrame` by arrival time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.489558Z",
     "iopub.status.busy": "2021-04-16T19:37:36.488816Z",
     "iopub.status.idle": "2021-04-16T19:37:36.497799Z",
     "shell.execute_reply": "2021-04-16T19:37:36.497392Z"
    }
   },
   "outputs": [],
   "source": [
    "obs = obs.sort_values(by='start', ignore_index=True)\n",
    "obs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that several of the lifelines extend past the observation window of 8 weeks.\n",
    "So if we observed this system at the beginning of Week 8, we would have incomplete information.\n",
    "Specifically, we would not know the future adoption times for Dogs 6, 7, and 8.\n",
    "\n",
    "I'll simulate this incomplete data by identifying the lifelines that extend past the observation window:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.501084Z",
     "iopub.status.busy": "2021-04-16T19:37:36.500671Z",
     "iopub.status.idle": "2021-04-16T19:37:36.503464Z",
     "shell.execute_reply": "2021-04-16T19:37:36.503887Z"
    }
   },
   "outputs": [],
   "source": [
    "censored = obs['end'] > 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`censored` is a  Boolean Series that is `True` for lifelines that extend past Week 8.\n",
    "\n",
    "Data that is not available is sometimes called \"censored\" in the sense that it is hidden from us.\n",
    "But in this case it is hidden because we don't know the future, not because someone is censoring it.\n",
    "\n",
    "For the lifelines that are censored, I'll modify `end` to indicate when they are last observed and `status` to indicate that the observation is incomplete. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.509515Z",
     "iopub.status.busy": "2021-04-16T19:37:36.508879Z",
     "iopub.status.idle": "2021-04-16T19:37:36.510902Z",
     "shell.execute_reply": "2021-04-16T19:37:36.511348Z"
    }
   },
   "outputs": [],
   "source": [
    "obs.loc[censored, 'end'] = 8\n",
    "obs.loc[censored, 'status'] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot a \"lifeline\" for each dog, showing the arrival and departure times on a time line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.517051Z",
     "iopub.status.busy": "2021-04-16T19:37:36.516471Z",
     "iopub.status.idle": "2021-04-16T19:37:36.519153Z",
     "shell.execute_reply": "2021-04-16T19:37:36.518641Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plot_lifelines(obs):\n",
    "    \"\"\"Plot a line for each observation.\n",
    "    \n",
    "    obs: DataFrame\n",
    "    \"\"\"\n",
    "    for y, row in obs.iterrows():\n",
    "        start = row['start']\n",
    "        end = row['end']\n",
    "        status = row['status']\n",
    "        \n",
    "        if status == 0:\n",
    "            # ongoing\n",
    "            plt.hlines(y, start, end, color='C0')\n",
    "        else:\n",
    "            # complete\n",
    "            plt.hlines(y, start, end, color='C1')\n",
    "            plt.plot(end, y, marker='o', color='C1')\n",
    "            \n",
    "    decorate(xlabel='Time (weeks)',\n",
    "             ylabel='Dog index',\n",
    "             title='Lifelines showing censored and uncensored observations')\n",
    "\n",
    "    plt.gca().invert_yaxis()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.541732Z",
     "iopub.status.busy": "2021-04-16T19:37:36.537157Z",
     "iopub.status.idle": "2021-04-16T19:37:36.697499Z",
     "shell.execute_reply": "2021-04-16T19:37:36.697069Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_lifelines(obs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And I'll add one more column to the table, which contains the duration of the observed parts of the lifelines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.701567Z",
     "iopub.status.busy": "2021-04-16T19:37:36.701132Z",
     "iopub.status.idle": "2021-04-16T19:37:36.703232Z",
     "shell.execute_reply": "2021-04-16T19:37:36.702836Z"
    }
   },
   "outputs": [],
   "source": [
    "obs['T'] = obs['end'] - obs['start']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we have simulated is the data that would be available at the beginning of Week 8."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using Incomplete Data\n",
    "\n",
    "Now, let's see how we can use both kinds of data, complete and incomplete, to infer the parameters of the distribution of residence times.\n",
    "\n",
    "First I'll split the data into two sets: `data1` contains residence times for dogs whose arrival and departure times are known; `data2` contains incomplete residence times for dogs who were not adopted during the observation interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.707495Z",
     "iopub.status.busy": "2021-04-16T19:37:36.706904Z",
     "iopub.status.idle": "2021-04-16T19:37:36.709424Z",
     "shell.execute_reply": "2021-04-16T19:37:36.708865Z"
    }
   },
   "outputs": [],
   "source": [
    "data1 = obs.loc[~censored, 'T']\n",
    "data2 = obs.loc[censored, 'T']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.714502Z",
     "iopub.status.busy": "2021-04-16T19:37:36.713805Z",
     "iopub.status.idle": "2021-04-16T19:37:36.717518Z",
     "shell.execute_reply": "2021-04-16T19:37:36.716979Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "data1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.722147Z",
     "iopub.status.busy": "2021-04-16T19:37:36.721569Z",
     "iopub.status.idle": "2021-04-16T19:37:36.724025Z",
     "shell.execute_reply": "2021-04-16T19:37:36.724393Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "data2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the complete data, we can use `update_weibull`, which uses the PDF of the Weibull distribution to compute the likelihood of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.727904Z",
     "iopub.status.busy": "2021-04-16T19:37:36.727229Z",
     "iopub.status.idle": "2021-04-16T19:37:36.745191Z",
     "shell.execute_reply": "2021-04-16T19:37:36.744717Z"
    }
   },
   "outputs": [],
   "source": [
    "posterior1 = update_weibull(prior, data1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the incomplete data, we have to think a little harder.\n",
    "At the end of the observation interval, we don't know what the residence time will be, but we can put a lower bound on it; that is, we can say that the residence time will be greater than `T`.\n",
    "\n",
    "And that means that we can compute the likelihood of the data using the survival function, which is the probability that a value from the distribution exceeds `T`.\n",
    "\n",
    "The following function is identical to `update_weibull` except that it uses `sf`, which computes the survival function, rather than `pdf`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.748842Z",
     "iopub.status.busy": "2021-04-16T19:37:36.748416Z",
     "iopub.status.idle": "2021-04-16T19:37:36.750080Z",
     "shell.execute_reply": "2021-04-16T19:37:36.750423Z"
    }
   },
   "outputs": [],
   "source": [
    "def update_weibull_incomplete(prior, data):\n",
    "    \"\"\"Update the prior using incomplete data.\"\"\"\n",
    "    lam_mesh, k_mesh, data_mesh = np.meshgrid(\n",
    "        prior.columns, prior.index, data)\n",
    "    \n",
    "    # evaluate the survival function\n",
    "    probs = weibull_dist(lam_mesh, k_mesh).sf(data_mesh)\n",
    "    likelihood = probs.prod(axis=2)\n",
    "\n",
    "    posterior = prior * likelihood\n",
    "    normalize(posterior)\n",
    "\n",
    "    return posterior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the update with the incomplete data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.753329Z",
     "iopub.status.busy": "2021-04-16T19:37:36.752847Z",
     "iopub.status.idle": "2021-04-16T19:37:36.760899Z",
     "shell.execute_reply": "2021-04-16T19:37:36.760377Z"
    }
   },
   "outputs": [],
   "source": [
    "posterior2 = update_weibull_incomplete(posterior1, data2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's what the joint posterior distribution looks like after both updates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.778276Z",
     "iopub.status.busy": "2021-04-16T19:37:36.777834Z",
     "iopub.status.idle": "2021-04-16T19:37:36.913674Z",
     "shell.execute_reply": "2021-04-16T19:37:36.913302Z"
    }
   },
   "outputs": [],
   "source": [
    "plot_contour(posterior2)\n",
    "decorate(title='Posterior joint distribution, incomplete data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to the previous contour plot, it looks like the range of likely values for $\\lambda$ is substantially wider.\n",
    "We can see that more clearly by looking at the marginal distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.917569Z",
     "iopub.status.busy": "2021-04-16T19:37:36.916847Z",
     "iopub.status.idle": "2021-04-16T19:37:36.920372Z",
     "shell.execute_reply": "2021-04-16T19:37:36.919902Z"
    }
   },
   "outputs": [],
   "source": [
    "posterior_lam2 = marginal(posterior2, 0)\n",
    "posterior_k2 = marginal(posterior2, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the posterior marginal distribution for $\\lambda$ compared to the distribution we got using all complete data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:36.960290Z",
     "iopub.status.busy": "2021-04-16T19:37:36.959500Z",
     "iopub.status.idle": "2021-04-16T19:37:37.098684Z",
     "shell.execute_reply": "2021-04-16T19:37:37.099061Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "posterior_lam.plot(color='C5', label='All complete',\n",
    "                   linestyle='dashed')\n",
    "posterior_lam2.plot(color='C2', label='Some censored')\n",
    "\n",
    "decorate(xlabel='lambda',\n",
    "         ylabel='PDF', \n",
    "         title='Marginal posterior distribution of lambda')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distribution with some incomplete data is substantially wider.\n",
    "\n",
    "As an aside, notice that the posterior distribution does not come all the way to 0 on the right side.\n",
    "That suggests that the range of the prior distribution is not wide enough to cover the most likely values for this parameter.\n",
    "If I were concerned about making this distribution more accurate, I would go back and run the update again with a wider prior.\n",
    "\n",
    "Here's the posterior marginal distribution for $k$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.139823Z",
     "iopub.status.busy": "2021-04-16T19:37:37.139316Z",
     "iopub.status.idle": "2021-04-16T19:37:37.272916Z",
     "shell.execute_reply": "2021-04-16T19:37:37.272485Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "posterior_k.plot(color='C5', label='All complete',\n",
    "                   linestyle='dashed')\n",
    "posterior_k2.plot(color='C12', label='Some censored')\n",
    "\n",
    "decorate(xlabel='k',\n",
    "         ylabel='PDF', \n",
    "         title='Posterior marginal distribution of k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, the marginal distribution is shifted to the left when we have incomplete data, but it is not substantially wider.\n",
    "\n",
    "In summary, we have seen how to combine complete and incomplete data to estimate the parameters of a Weibull distribution, which is useful in many real-world scenarios where some of the data are censored.\n",
    "\n",
    "In general, the posterior distributions are wider when we have incomplete data, because less information leads to more uncertainty.\n",
    "\n",
    "This example is based on data I generated; in the next section we'll do a similar analysis with real data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Light Bulbs\n",
    "\n",
    "In 2007 [researchers ran an experiment](https://www.researchgate.net/publication/225450325_Renewal_Rate_of_Filament_Lamps_Theory_and_Experiment) to characterize the distribution of lifetimes for light bulbs.\n",
    "Here is their description of the experiment:\n",
    "\n",
    "> An assembly of 50 new Philips (India) lamps with the rating 40 W, 220 V (AC) was taken and installed in the horizontal orientation and uniformly distributed over a lab area 11 m x 7 m.\n",
    ">\n",
    "> The assembly was monitored at regular intervals of 12 h to look for failures. The instants of recorded failures were [recorded] and a total of 32 data points were obtained such that even the last bulb failed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.276799Z",
     "iopub.status.busy": "2021-04-16T19:37:37.276346Z",
     "iopub.status.idle": "2021-04-16T19:37:37.278661Z",
     "shell.execute_reply": "2021-04-16T19:37:37.278236Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "download('https://gist.github.com/epogrebnyak/7933e16c0ad215742c4c104be4fbdeb1/raw/c932bc5b6aa6317770c4cbf43eb591511fec08f9/lamps.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can load the data into a `DataFrame` like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.282359Z",
     "iopub.status.busy": "2021-04-16T19:37:37.281849Z",
     "iopub.status.idle": "2021-04-16T19:37:37.291574Z",
     "shell.execute_reply": "2021-04-16T19:37:37.292000Z"
    }
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('lamps.csv', index_col=0)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Column `h` contains the times when bulbs failed in hours; Column `f` contains the number of bulbs that failed at each time.\n",
    "We can represent these values and frequencies using a `Pmf`, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.296777Z",
     "iopub.status.busy": "2021-04-16T19:37:37.296262Z",
     "iopub.status.idle": "2021-04-16T19:37:37.298609Z",
     "shell.execute_reply": "2021-04-16T19:37:37.298966Z"
    }
   },
   "outputs": [],
   "source": [
    "from empiricaldist import Pmf\n",
    "\n",
    "pmf_bulb = Pmf(df['f'].to_numpy(), df['h'])\n",
    "pmf_bulb.normalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because of the design of this experiment, we can consider the data to be a representative sample from the distribution of lifetimes, at least for light bulbs that are lit continuously."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The average lifetime is about 1400 h."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.302433Z",
     "iopub.status.busy": "2021-04-16T19:37:37.301719Z",
     "iopub.status.idle": "2021-04-16T19:37:37.305072Z",
     "shell.execute_reply": "2021-04-16T19:37:37.304694Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "pmf_bulb.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming that these data are well modeled by a Weibull distribution, let's estimate the parameters that fit the data.\n",
    "Again, I'll start with uniform priors for $\\lambda$ and $k$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.309226Z",
     "iopub.status.busy": "2021-04-16T19:37:37.308759Z",
     "iopub.status.idle": "2021-04-16T19:37:37.310609Z",
     "shell.execute_reply": "2021-04-16T19:37:37.311018Z"
    }
   },
   "outputs": [],
   "source": [
    "lams = np.linspace(1000, 2000, num=51)\n",
    "prior_lam = make_uniform(lams, name='lambda')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.315203Z",
     "iopub.status.busy": "2021-04-16T19:37:37.314719Z",
     "iopub.status.idle": "2021-04-16T19:37:37.316795Z",
     "shell.execute_reply": "2021-04-16T19:37:37.317307Z"
    }
   },
   "outputs": [],
   "source": [
    "ks = np.linspace(1, 10, num=51)\n",
    "prior_k = make_uniform(ks, name='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, there are 51 values in the prior distribution, rather than the usual 101.  That's because we are going to use the posterior distributions to do some computationally-intensive calculations.\n",
    "They will run faster with fewer values, but the results will be less precise.\n",
    "\n",
    "As usual, we can use `make_joint` to make the prior joint distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.321446Z",
     "iopub.status.busy": "2021-04-16T19:37:37.320851Z",
     "iopub.status.idle": "2021-04-16T19:37:37.322797Z",
     "shell.execute_reply": "2021-04-16T19:37:37.323250Z"
    }
   },
   "outputs": [],
   "source": [
    "prior_bulb = make_joint(prior_lam, prior_k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although we have data for 50 light bulbs, there are only 32 unique lifetimes in the dataset.  For the update, it is convenient to express the data in the form of 50 lifetimes, with each lifetime repeated the given number of times.\n",
    "We can use `np.repeat` to transform the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.327777Z",
     "iopub.status.busy": "2021-04-16T19:37:37.327199Z",
     "iopub.status.idle": "2021-04-16T19:37:37.330359Z",
     "shell.execute_reply": "2021-04-16T19:37:37.330844Z"
    }
   },
   "outputs": [],
   "source": [
    "data_bulb = np.repeat(df['h'], df['f'])\n",
    "len(data_bulb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use `update_weibull` to do the update."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.334786Z",
     "iopub.status.busy": "2021-04-16T19:37:37.334174Z",
     "iopub.status.idle": "2021-04-16T19:37:37.363731Z",
     "shell.execute_reply": "2021-04-16T19:37:37.364153Z"
    }
   },
   "outputs": [],
   "source": [
    "posterior_bulb = update_weibull(prior_bulb, data_bulb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the posterior joint distribution looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.376746Z",
     "iopub.status.busy": "2021-04-16T19:37:37.372073Z",
     "iopub.status.idle": "2021-04-16T19:37:37.531370Z",
     "shell.execute_reply": "2021-04-16T19:37:37.530980Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_contour(posterior_bulb)\n",
    "decorate(title='Joint posterior distribution, light bulbs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To summarize this joint posterior distribution, we'll compute the posterior mean lifetime."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior Means\n",
    "\n",
    "To compute the posterior mean of a joint distribution, we'll make a mesh that contains the values of $\\lambda$ and $k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.534652Z",
     "iopub.status.busy": "2021-04-16T19:37:37.534197Z",
     "iopub.status.idle": "2021-04-16T19:37:37.536544Z",
     "shell.execute_reply": "2021-04-16T19:37:37.536146Z"
    }
   },
   "outputs": [],
   "source": [
    "lam_mesh, k_mesh = np.meshgrid(\n",
    "    prior_bulb.columns, prior_bulb.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for each pair of parameters we'll use `weibull_dist` to compute the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.540341Z",
     "iopub.status.busy": "2021-04-16T19:37:37.539900Z",
     "iopub.status.idle": "2021-04-16T19:37:37.543970Z",
     "shell.execute_reply": "2021-04-16T19:37:37.544466Z"
    }
   },
   "outputs": [],
   "source": [
    "means = weibull_dist(lam_mesh, k_mesh).mean()\n",
    "means.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is an array with the same dimensions as the joint distribution.\n",
    "\n",
    "Now we need to weight each mean with the corresponding probability from the joint posterior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.549477Z",
     "iopub.status.busy": "2021-04-16T19:37:37.548505Z",
     "iopub.status.idle": "2021-04-16T19:37:37.551293Z",
     "shell.execute_reply": "2021-04-16T19:37:37.550670Z"
    }
   },
   "outputs": [],
   "source": [
    "prod = means * posterior_bulb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we compute the sum of the weighted means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.556060Z",
     "iopub.status.busy": "2021-04-16T19:37:37.555354Z",
     "iopub.status.idle": "2021-04-16T19:37:37.559241Z",
     "shell.execute_reply": "2021-04-16T19:37:37.558666Z"
    }
   },
   "outputs": [],
   "source": [
    "prod.to_numpy().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the posterior distribution, we think the mean lifetime is about 1413 hours.\n",
    "\n",
    "The following function encapsulates these steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.564528Z",
     "iopub.status.busy": "2021-04-16T19:37:37.563826Z",
     "iopub.status.idle": "2021-04-16T19:37:37.566786Z",
     "shell.execute_reply": "2021-04-16T19:37:37.566099Z"
    }
   },
   "outputs": [],
   "source": [
    "def joint_weibull_mean(joint):\n",
    "    \"\"\"Compute the mean of a joint distribution of Weibulls.\"\"\"\n",
    "    lam_mesh, k_mesh = np.meshgrid(\n",
    "        joint.columns, joint.index)\n",
    "    means = weibull_dist(lam_mesh, k_mesh).mean()\n",
    "    prod = means * joint\n",
    "    return prod.to_numpy().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Incomplete Information\n",
    "\n",
    "The previous update was not quite right, because it assumed each light bulb died at the instant we observed it.  \n",
    "According to the report, the researchers only checked the bulbs every 12 hours.  So if they see that a bulb has died, they know only that it died during the 12 hours since the last check.\n",
    "\n",
    "It is more strictly correct to use the following update function, which uses the CDF of the Weibull distribution to compute the probability that a bulb dies during a given 12 hour interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.572709Z",
     "iopub.status.busy": "2021-04-16T19:37:37.572014Z",
     "iopub.status.idle": "2021-04-16T19:37:37.574370Z",
     "shell.execute_reply": "2021-04-16T19:37:37.574922Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def update_weibull_between(prior, data, dt=12):\n",
    "    \"\"\"Update the prior based on data.\"\"\"\n",
    "    lam_mesh, k_mesh, data_mesh = np.meshgrid(\n",
    "        prior.columns, prior.index, data)\n",
    "    dist = weibull_dist(lam_mesh, k_mesh)\n",
    "    cdf1 = dist.cdf(data_mesh)\n",
    "    cdf2 = dist.cdf(data_mesh-12)\n",
    "    likelihood = (cdf1 - cdf2).prod(axis=2)\n",
    "\n",
    "    posterior = prior * likelihood\n",
    "    normalize(posterior)\n",
    "\n",
    "    return posterior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The probability that a value falls in an interval is the difference between the CDF at the beginning and end of the interval.\n",
    "\n",
    "Here's how we run the update."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.580024Z",
     "iopub.status.busy": "2021-04-16T19:37:37.579225Z",
     "iopub.status.idle": "2021-04-16T19:37:37.617913Z",
     "shell.execute_reply": "2021-04-16T19:37:37.617496Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "posterior_bulb2 = update_weibull_between(prior_bulb, data_bulb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "And here are the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.637793Z",
     "iopub.status.busy": "2021-04-16T19:37:37.634750Z",
     "iopub.status.idle": "2021-04-16T19:37:37.784827Z",
     "shell.execute_reply": "2021-04-16T19:37:37.785224Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_contour(posterior_bulb2)\n",
    "decorate(title='Joint posterior distribution, light bulbs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Visually this result is almost identical to what we got using the PDF.\n",
    "And that's good news, because it suggests that using the PDF can be a good approximation even if it's not strictly correct.\n",
    "\n",
    "To see whether it makes any difference at all, let's check the posterior means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.789388Z",
     "iopub.status.busy": "2021-04-16T19:37:37.788691Z",
     "iopub.status.idle": "2021-04-16T19:37:37.793773Z",
     "shell.execute_reply": "2021-04-16T19:37:37.794199Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "joint_weibull_mean(posterior_bulb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.799345Z",
     "iopub.status.busy": "2021-04-16T19:37:37.798512Z",
     "iopub.status.idle": "2021-04-16T19:37:37.803590Z",
     "shell.execute_reply": "2021-04-16T19:37:37.804088Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "joint_weibull_mean(posterior_bulb2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "When we take into account the 12-hour interval between observations, the posterior mean is about 6 hours less.\n",
    "And that makes sense: if we assume that a bulb is equally likely to expire at any point in the interval, the average would be the midpoint of the interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior Predictive Distribution\n",
    "\n",
    "Suppose you install 100 light bulbs of the kind in the previous section, and you come back to check on them after 1000 hours.  Based on the posterior distribution we just computed, what is the distribution of the number of bulbs you find dead?\n",
    "\n",
    "If we knew the parameters of the Weibull distribution for sure, the answer would be a binomial distribution.\n",
    "\n",
    "For example, if we know that $\\lambda=1550$ and $k=4.25$, we can use `weibull_dist` to compute the probability that a bulb dies before you return:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.810920Z",
     "iopub.status.busy": "2021-04-16T19:37:37.810150Z",
     "iopub.status.idle": "2021-04-16T19:37:37.814082Z",
     "shell.execute_reply": "2021-04-16T19:37:37.813392Z"
    }
   },
   "outputs": [],
   "source": [
    "lam = 1550\n",
    "k = 4.25\n",
    "t = 1000\n",
    "\n",
    "prob_dead = weibull_dist(lam, k).cdf(t)\n",
    "prob_dead"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If there are 100 bulbs and each has this probability of dying, the number of dead bulbs follows a binomial distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.818684Z",
     "iopub.status.busy": "2021-04-16T19:37:37.818078Z",
     "iopub.status.idle": "2021-04-16T19:37:37.820328Z",
     "shell.execute_reply": "2021-04-16T19:37:37.819814Z"
    }
   },
   "outputs": [],
   "source": [
    "from utils import make_binomial\n",
    "\n",
    "n = 100\n",
    "p = prob_dead\n",
    "dist_num_dead = make_binomial(n, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "And here's what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:37.902574Z",
     "iopub.status.busy": "2021-04-16T19:37:37.892613Z",
     "iopub.status.idle": "2021-04-16T19:37:38.033693Z",
     "shell.execute_reply": "2021-04-16T19:37:38.034038Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "dist_num_dead.plot(label='known parameters')\n",
    "\n",
    "decorate(xlabel='Number of dead bulbs',\n",
    "         ylabel='PMF',\n",
    "         title='Predictive distribution with known parameters')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But that's based on the assumption that we know $\\lambda$ and $k$, and we don't.\n",
    "Instead, we have a posterior distribution that contains possible values of these parameters and their probabilities.\n",
    "\n",
    "So the posterior predictive distribution is not a single binomial; instead it is a mixture of binomials, weighted with the posterior probabilities.\n",
    "\n",
    "We can use `make_mixture` to compute the posterior predictive distribution.  \n",
    "It doesn't work with joint distributions, but we can convert the `DataFrame` that represents a joint distribution to a `Series`, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:38.037138Z",
     "iopub.status.busy": "2021-04-16T19:37:38.036652Z",
     "iopub.status.idle": "2021-04-16T19:37:38.042709Z",
     "shell.execute_reply": "2021-04-16T19:37:38.042213Z"
    }
   },
   "outputs": [],
   "source": [
    "posterior_series = posterior_bulb.stack()\n",
    "posterior_series.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a `Series` with a `MultiIndex` that contains two \"levels\": the first level contains the values of `k`; the second contains the values of `lam`.\n",
    "\n",
    "With the posterior in this form, we can iterate through the possible parameters and compute a predictive distribution for each pair."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:38.059067Z",
     "iopub.status.busy": "2021-04-16T19:37:38.056752Z",
     "iopub.status.idle": "2021-04-16T19:37:40.599814Z",
     "shell.execute_reply": "2021-04-16T19:37:40.600201Z"
    }
   },
   "outputs": [],
   "source": [
    "pmf_seq = []\n",
    "for (k, lam) in posterior_series.index:\n",
    "    prob_dead = weibull_dist(lam, k).cdf(t)\n",
    "    pmf = make_binomial(n, prob_dead)\n",
    "    pmf_seq.append(pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use `make_mixture`, passing as parameters the posterior probabilities in `posterior_series` and the sequence of binomial distributions in `pmf_seq`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:40.673531Z",
     "iopub.status.busy": "2021-04-16T19:37:40.638022Z",
     "iopub.status.idle": "2021-04-16T19:37:40.898812Z",
     "shell.execute_reply": "2021-04-16T19:37:40.899237Z"
    }
   },
   "outputs": [],
   "source": [
    "from utils import make_mixture\n",
    "\n",
    "post_pred = make_mixture(posterior_series, pmf_seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the posterior predictive distribution looks like, compared to the binomial distribution we computed with known parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:40.938328Z",
     "iopub.status.busy": "2021-04-16T19:37:40.915391Z",
     "iopub.status.idle": "2021-04-16T19:37:41.080448Z",
     "shell.execute_reply": "2021-04-16T19:37:41.080783Z"
    },
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "dist_num_dead.plot(label='known parameters')\n",
    "post_pred.plot(label='unknown parameters')\n",
    "decorate(xlabel='Number of dead bulbs',\n",
    "         ylabel='PMF',\n",
    "         title='Posterior predictive distribution')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior predictive distribution is wider because it represents our uncertainty about the parameters as well as our uncertainty about the number of dead bulbs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "This chapter introduces survival analysis, which is used to answer questions about the time until an event, and the Weibull distribution, which is a good model for \"lifetimes\" (broadly interpreted) in a number of domains.\n",
    "\n",
    "We used joint distributions to represent prior probabilities for the parameters of the Weibull distribution, and we updated them three ways: knowing the exact duration of a lifetime, knowing a lower bound, and knowing that a lifetime fell in a given interval.\n",
    "\n",
    "These examples demonstrate a feature of Bayesian methods: they can be adapted to handle incomplete, or \"censored\", data with only small changes.  As an exercise, you'll have a chance to work with one more type of censored data, when we are given an upper bound on a lifetime.\n",
    "\n",
    "The methods in this chapter work with any distribution with two parameters.\n",
    "In the exercises, you'll have a chance to estimate the parameters of a two-parameter gamma distribution, which is used to describe a variety of natural phenomena.\n",
    "\n",
    "And in the next chapter we'll move on to models with three parameters!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Exercise:** Using data about the lifetimes of light bulbs, we computed the posterior distribution from the parameters of a Weibull distribution, $\\lambda$ and $k$, and the posterior predictive distribution for the number of dead bulbs, out of 100, after 1000 hours.\n",
    "\n",
    "Now suppose you do the experiment:  You install 100 light bulbs, come back after 1000 hours, and find 20 dead light bulbs.\n",
    "Update the posterior distribution based on this data.\n",
    "How much does it change the posterior mean?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "tags": []
   },
   "source": [
    "Suggestions:\n",
    "\n",
    "1. Use a mesh grid to compute the probability of finding a bulb dead after 1000 hours for each pair of parameters.\n",
    "\n",
    "2. For each of those probabilities, compute the likelihood of finding 20 dead bulbs out of 100.\n",
    "\n",
    "3. Use those likelihoods to update the posterior distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.085218Z",
     "iopub.status.busy": "2021-04-16T19:37:41.084757Z",
     "iopub.status.idle": "2021-04-16T19:37:41.088304Z",
     "shell.execute_reply": "2021-04-16T19:37:41.088622Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.093414Z",
     "iopub.status.busy": "2021-04-16T19:37:41.092096Z",
     "iopub.status.idle": "2021-04-16T19:37:41.095793Z",
     "shell.execute_reply": "2021-04-16T19:37:41.096159Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.117695Z",
     "iopub.status.busy": "2021-04-16T19:37:41.112432Z",
     "iopub.status.idle": "2021-04-16T19:37:41.236742Z",
     "shell.execute_reply": "2021-04-16T19:37:41.237307Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.241318Z",
     "iopub.status.busy": "2021-04-16T19:37:41.240598Z",
     "iopub.status.idle": "2021-04-16T19:37:41.246161Z",
     "shell.execute_reply": "2021-04-16T19:37:41.246507Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** In this exercise, we'll use one month of data to estimate the parameters of a distribution that describes daily rainfall in Seattle.\n",
    "Then we'll compute the posterior predictive distribution for daily rainfall and use it to estimate the probability of a rare event, like more than 1.5 inches of rain in a day.\n",
    "\n",
    "According to hydrologists, the distribution of total daily rainfall (for days with rain) is well modeled by a two-parameter\n",
    "gamma distribution.\n",
    "\n",
    "When we worked with the one-parameter gamma distribution in <<_TheGammaDistribution>>, we used the Greek letter $\\alpha$ for the parameter.\n",
    "\n",
    "For the two-parameter gamma distribution, we will use $k$ for the \"shape parameter\", which determines the shape of the distribution, and the Greek letter $\\theta$ or `theta` for the \"scale parameter\". "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The following function takes these parameters and returns a `gamma` object from SciPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.250070Z",
     "iopub.status.busy": "2021-04-16T19:37:41.249652Z",
     "iopub.status.idle": "2021-04-16T19:37:41.252217Z",
     "shell.execute_reply": "2021-04-16T19:37:41.252548Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import scipy.stats\n",
    "\n",
    "def gamma_dist(k, theta):\n",
    "    \"\"\"Makes a gamma object.\n",
    "    \n",
    "    k: shape parameter\n",
    "    theta: scale parameter\n",
    "    \n",
    "    returns: gamma object\n",
    "    \"\"\"\n",
    "    return scipy.stats.gamma(k, scale=theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Now we need some data.\n",
    "The following cell downloads data I collected from the National Oceanic and Atmospheric Administration ([NOAA](http://www.ncdc.noaa.gov/cdo-web/search)) for Seattle, Washington in May 2020."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.256571Z",
     "iopub.status.busy": "2021-04-16T19:37:41.256055Z",
     "iopub.status.idle": "2021-04-16T19:37:41.258114Z",
     "shell.execute_reply": "2021-04-16T19:37:41.257746Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Load the data file\n",
    "\n",
    "download('https://github.com/AllenDowney/ThinkBayes2/raw/master/data/2203951.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Now we can load it into a `DataFrame`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.261365Z",
     "iopub.status.busy": "2021-04-16T19:37:41.260827Z",
     "iopub.status.idle": "2021-04-16T19:37:41.274723Z",
     "shell.execute_reply": "2021-04-16T19:37:41.274264Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "weather = pd.read_csv('2203951.csv')\n",
    "weather.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "I'll make a Boolean Series to indicate which days it rained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.278994Z",
     "iopub.status.busy": "2021-04-16T19:37:41.278485Z",
     "iopub.status.idle": "2021-04-16T19:37:41.280829Z",
     "shell.execute_reply": "2021-04-16T19:37:41.281256Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "rained = weather['PRCP'] > 0\n",
    "rained.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "And select the total rainfall on the days it rained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.288529Z",
     "iopub.status.busy": "2021-04-16T19:37:41.287802Z",
     "iopub.status.idle": "2021-04-16T19:37:41.291569Z",
     "shell.execute_reply": "2021-04-16T19:37:41.291072Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "prcp = weather.loc[rained, 'PRCP']\n",
    "prcp.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Here's what the CDF of the data looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.327363Z",
     "iopub.status.busy": "2021-04-16T19:37:41.311942Z",
     "iopub.status.idle": "2021-04-16T19:37:41.416224Z",
     "shell.execute_reply": "2021-04-16T19:37:41.416543Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "cdf_data = Cdf.from_seq(prcp)\n",
    "cdf_data.plot()\n",
    "decorate(xlabel='Total rainfall (in)',\n",
    "         ylabel='CDF',\n",
    "         title='Distribution of rainfall on days it rained')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The maximum is 1.14 inches of rain is one day.\n",
    "To estimate the probability of more than 1.5 inches, we need to extrapolate from the data we have, so our estimate will depend on whether the gamma distribution is really a good model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I suggest you proceed in the following steps:\n",
    "\n",
    "1. Construct a prior distribution for the parameters of the gamma distribution.  Note that $k$ and $\\theta$ must be greater than 0.\n",
    "\n",
    "2. Use the observed rainfalls to update the distribution of parameters.\n",
    "\n",
    "3. Compute the posterior predictive distribution of rainfall, and use it to estimate the probability of getting more than 1.5 inches of rain in one day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.420568Z",
     "iopub.status.busy": "2021-04-16T19:37:41.420152Z",
     "iopub.status.idle": "2021-04-16T19:37:41.424495Z",
     "shell.execute_reply": "2021-04-16T19:37:41.424149Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.428911Z",
     "iopub.status.busy": "2021-04-16T19:37:41.428490Z",
     "iopub.status.idle": "2021-04-16T19:37:41.431090Z",
     "shell.execute_reply": "2021-04-16T19:37:41.431568Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.435824Z",
     "iopub.status.busy": "2021-04-16T19:37:41.435397Z",
     "iopub.status.idle": "2021-04-16T19:37:41.438583Z",
     "shell.execute_reply": "2021-04-16T19:37:41.438240Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.441827Z",
     "iopub.status.busy": "2021-04-16T19:37:41.441245Z",
     "iopub.status.idle": "2021-04-16T19:37:41.444095Z",
     "shell.execute_reply": "2021-04-16T19:37:41.444570Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.448288Z",
     "iopub.status.busy": "2021-04-16T19:37:41.447637Z",
     "iopub.status.idle": "2021-04-16T19:37:41.450859Z",
     "shell.execute_reply": "2021-04-16T19:37:41.450497Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.455536Z",
     "iopub.status.busy": "2021-04-16T19:37:41.454731Z",
     "iopub.status.idle": "2021-04-16T19:37:41.462470Z",
     "shell.execute_reply": "2021-04-16T19:37:41.462038Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.466277Z",
     "iopub.status.busy": "2021-04-16T19:37:41.465631Z",
     "iopub.status.idle": "2021-04-16T19:37:41.468004Z",
     "shell.execute_reply": "2021-04-16T19:37:41.468364Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.472830Z",
     "iopub.status.busy": "2021-04-16T19:37:41.472258Z",
     "iopub.status.idle": "2021-04-16T19:37:41.475159Z",
     "shell.execute_reply": "2021-04-16T19:37:41.474721Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.492554Z",
     "iopub.status.busy": "2021-04-16T19:37:41.488649Z",
     "iopub.status.idle": "2021-04-16T19:37:41.647173Z",
     "shell.execute_reply": "2021-04-16T19:37:41.646742Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.651740Z",
     "iopub.status.busy": "2021-04-16T19:37:41.651267Z",
     "iopub.status.idle": "2021-04-16T19:37:41.653711Z",
     "shell.execute_reply": "2021-04-16T19:37:41.653349Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.689643Z",
     "iopub.status.busy": "2021-04-16T19:37:41.675794Z",
     "iopub.status.idle": "2021-04-16T19:37:41.824132Z",
     "shell.execute_reply": "2021-04-16T19:37:41.824491Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.831274Z",
     "iopub.status.busy": "2021-04-16T19:37:41.830519Z",
     "iopub.status.idle": "2021-04-16T19:37:41.834487Z",
     "shell.execute_reply": "2021-04-16T19:37:41.835136Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:41.855628Z",
     "iopub.status.busy": "2021-04-16T19:37:41.854320Z",
     "iopub.status.idle": "2021-04-16T19:37:42.012073Z",
     "shell.execute_reply": "2021-04-16T19:37:42.011701Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:42.016693Z",
     "iopub.status.busy": "2021-04-16T19:37:42.016116Z",
     "iopub.status.idle": "2021-04-16T19:37:42.018530Z",
     "shell.execute_reply": "2021-04-16T19:37:42.018903Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:42.022546Z",
     "iopub.status.busy": "2021-04-16T19:37:42.021767Z",
     "iopub.status.idle": "2021-04-16T19:37:42.027862Z",
     "shell.execute_reply": "2021-04-16T19:37:42.027491Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:42.030840Z",
     "iopub.status.busy": "2021-04-16T19:37:42.030357Z",
     "iopub.status.idle": "2021-04-16T19:37:42.033846Z",
     "shell.execute_reply": "2021-04-16T19:37:42.033353Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:42.064699Z",
     "iopub.status.busy": "2021-04-16T19:37:42.052712Z",
     "iopub.status.idle": "2021-04-16T19:37:45.138382Z",
     "shell.execute_reply": "2021-04-16T19:37:45.137916Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:45.155754Z",
     "iopub.status.busy": "2021-04-16T19:37:45.150818Z",
     "iopub.status.idle": "2021-04-16T19:37:45.515419Z",
     "shell.execute_reply": "2021-04-16T19:37:45.514852Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:45.537004Z",
     "iopub.status.busy": "2021-04-16T19:37:45.534520Z",
     "iopub.status.idle": "2021-04-16T19:37:45.659852Z",
     "shell.execute_reply": "2021-04-16T19:37:45.659355Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:45.664168Z",
     "iopub.status.busy": "2021-04-16T19:37:45.663749Z",
     "iopub.status.idle": "2021-04-16T19:37:45.667926Z",
     "shell.execute_reply": "2021-04-16T19:37:45.668412Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:45.671553Z",
     "iopub.status.busy": "2021-04-16T19:37:45.671143Z",
     "iopub.status.idle": "2021-04-16T19:37:45.675804Z",
     "shell.execute_reply": "2021-04-16T19:37:45.675337Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}