{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "b8051b2d", "metadata": {}, "outputs": [], "source": [ "import pymc as pm\n", "import numpy as np\n", "import arviz as az\n", "import pandas as pd" ] }, { "cell_type": "markdown", "id": "ea9efb0d", "metadata": {}, "source": [ "# Predicting Using Censored Data*" ] }, { "cell_type": "markdown", "id": "7e9007bd", "metadata": {}, "source": [ "## Icelandic volcano eruptions\n", "\n", "Adapted from [Unit 10: katla.odc](https://raw.githubusercontent.com/areding/6420-pymc/main/original_examples/Codes4Unit10/katla.odc).\n", "\n", "Data can be found [here](https://raw.githubusercontent.com/areding/6420-pymc/main/data/r.txt).\n", "\n", "```{epigraph}\n", "In April 2010, the volcano Eyjafjallajökull in Iceland erupted. The resulting ash cloud was blown towards Western Europe and caused severe disruption to air travel for the following few weeks. A report into the eruption and its impact ({cite:t}`UCLIcelandVolcano`) reviewed how well the risk had been managed. One question was whether potentially more devastating eruptions from the larger neighbouring volcano Katla can be predicted from a recent eruption of Eyjafjallajökull. The report provides the dates of all 18 eruptions of Katla since the year 1177, with a corresponding indicator of whether Eyjafjallajökull had also erupted within the previous year.\n", "\n", "-- {cite:t}`Lunn2012BugsBook` p. 254\n", "```\n", "\n", "Thanks to William Naramore, this example has been updated to work with PyMC v5!" ] }, { "cell_type": "code", "execution_count": 2, "id": "456933fb", "metadata": {}, "outputs": [], "source": [ "# fmt: off\n", "D = np.array(\n", " (1177, 1262, 1311, 1357, 1416, 1440, 1450, 1500, \n", " 1550, 1580, 1612, 1625, 1660, 1721, 1755, 1823, \n", " 1860, 1918, np.inf) # \n", ")\n", "# fmt: on\n", "\n", "# probabilities\n", "ps = [1, 5, 10, 50]\n", "\n", "# time between eruptions\n", "t = np.diff(D)\n", "t[t > 100] = 100" ] }, { "cell_type": "code", "execution_count": 6, "id": "81c33012", "metadata": { "tags": [ "hide-output" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag_grad...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [α, σ]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [16000/16000 00:01<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 1 seconds.\n" ] } ], "source": [ "with pm.Model() as m:\n", " α = pm.TruncatedNormal(\"α\", mu=0, sigma=5, lower=0) # v in BUGS model\n", "\n", " σ = pm.Gamma(\"σ\", 0.1, 0.1)\n", " λ = 1 / σ**α\n", " β = λ ** (-1 / α)\n", "\n", " _t = pm.Weibull.dist(α, β)\n", " pm.Censored(\"likelihood\", _t, lower=None, upper=100, observed=t)\n", "\n", " median = pm.Deterministic(\"median tte\", σ * np.log(2) ** (1 / α))\n", "\n", " for p in ps:\n", " pm.Deterministic(\n", " f\"p_erupt_{p}\", 1 - pm.math.exp((100 / σ) ** α - ((100 + p) / σ) ** α)\n", " )\n", "\n", " trace = pm.sample(3000, init=\"jitter+adapt_diag_grad\", target_accept=.95)" ] }, { "cell_type": "code", "execution_count": 7, "id": "18e87fe7-453d-4425-89ef-6cd4fbee68a5", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meansdhdi_3%hdi_97%mcse_meanmcse_sdess_bulkess_tailr_hat
α1.9370.3801.2082.6270.0050.0044801.05240.01.0
σ50.2776.44838.13262.4640.0860.0615584.05754.01.0
median tte41.4066.06429.79752.4940.0820.0585416.05553.01.0
p_erupt_10.0750.0310.0260.1330.0000.0005439.05645.01.0
p_erupt_50.3220.1110.1300.5270.0010.0015423.05693.01.0
p_erupt_100.5350.1470.2680.8070.0020.0015408.05657.01.0
p_erupt_500.9660.0530.8681.0000.0010.0015281.05581.01.0
\n", "
" ], "text/plain": [ " mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk \\\n", "α 1.937 0.380 1.208 2.627 0.005 0.004 4801.0 \n", "σ 50.277 6.448 38.132 62.464 0.086 0.061 5584.0 \n", "median tte 41.406 6.064 29.797 52.494 0.082 0.058 5416.0 \n", "p_erupt_1 0.075 0.031 0.026 0.133 0.000 0.000 5439.0 \n", "p_erupt_5 0.322 0.111 0.130 0.527 0.001 0.001 5423.0 \n", "p_erupt_10 0.535 0.147 0.268 0.807 0.002 0.001 5408.0 \n", "p_erupt_50 0.966 0.053 0.868 1.000 0.001 0.001 5281.0 \n", "\n", " ess_tail r_hat \n", "α 5240.0 1.0 \n", "σ 5754.0 1.0 \n", "median tte 5553.0 1.0 \n", "p_erupt_1 5645.0 1.0 \n", "p_erupt_5 5693.0 1.0 \n", "p_erupt_10 5657.0 1.0 \n", "p_erupt_50 5581.0 1.0 " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "az.summary(trace)" ] }, { "cell_type": "code", "execution_count": 5, "id": "633613a1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Last updated: Mon Nov 27 2023\n", "\n", "Python implementation: CPython\n", "Python version : 3.11.5\n", "IPython version : 8.15.0\n", "\n", "pytensor: 2.17.1\n", "\n", "pandas: 2.1.0\n", "arviz : 0.16.1\n", "pymc : 5.9.2\n", "numpy : 1.25.2\n", "\n" ] } ], "source": [ "%load_ext watermark\n", "%watermark -n -u -v -iv -p pytensor" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.5" } }, "nbformat": 4, "nbformat_minor": 5 }