Attribution of a single event to climate change.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can we confidently attribute observed extreme events to climate change? A statitical formulation (from Hansen _et al._, 2014) provides a way to construct confidence intervals for so-called attributable risk. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "from itertools import product\n",
    "\n",
    "import seaborn; seaborn.set_context(\"poster\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 0.5\n",
    "beta = 1.5\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The length of records are $T_1$ and $T_2$, the ratio of the two is $\\tau$. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "T1 = 20\n",
    "T2 = 10\n",
    "tau = T1/T2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The conditional probability of observing $x$ and $y$ is \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathrm{prob}(X=x, Y=y) = \\frac{[\\mu(T_1 + \\beta T_2)]^{x+y}\\exp{[-\\mu (T_1 + \\beta T_2)]}}{(x+y)!} \\times \\frac{(x+y)!}{x!y!}[\\tau / (\\tau+\\beta)]^x [\\beta/(\\tau + \\beta)]^y]$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def jointProb(x, y, mu, beta, T1, T2):\n",
    "    tau = T1/T2\n",
    "    xyfac = np.misc.factorial(x + y)\n",
    "    xfac = np.misc.factorial(x)\n",
    "    yfac = np.misc.factorial(y)\n",
    "    totalobs = np.power(mu*(T1 + beta * T2), x + y) * np.exp(-mu * (T1 + beta * T2)) / xyfac\n",
    "    condprob = (xyfac / (xfac * yfac) ) * np.power(tau/(tau + beta), x) * \\\n",
    "                np.power(beta/(tau + beta), y)\n",
    "    pr = totalobs * condprob\n",
    "    \n",
    "    return pr\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The attributable risk "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def attributable(x, y, tau):\n",
    "    return 1. - x / (tau * y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the confidence intervals $p_L$ and $p_U$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def confidence(x, y, alpha, tau):\n",
    "    \n",
    "    \n",
    "    if y == 0:\n",
    "        return None\n",
    "    \n",
    "    n = x + y \n",
    "    q = x/n\n",
    "    if x == 0:\n",
    "        p_U = 1 - np.power(alpha, 1./n)\n",
    "        p_L = 1. - (1. - np.power(alpha, 1./n)) / (tau * np.power(alpha, 1./n))\n",
    "        return p_L, p_U\n",
    "    \n",
    "    z = norm.ppf(alpha/2.)\n",
    "    \n",
    "    det = np.sqrt(q * (1 - q) + z*z/(4 * n))\n",
    "    a = (x + z*z/2) / (n + z*z) \n",
    "    b = z*np.sqrt(n) / (n + z*z)\n",
    "    L = a - b * det\n",
    "    U = a + b * det\n",
    "    \n",
    "    beta_U = tau * (1. - U) / U\n",
    "    beta_L = tau * (1. - L) / L\n",
    "    \n",
    "    p_L = 1. - 1./beta_L\n",
    "    p_U = 1. - 1./beta_U\n",
    "    \n",
    "    return p_L, p_U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_results(x, y):\n",
    "    tau = np.array([0.5, 1., 2., 3., 4., 5.])\n",
    "    alpha = 0.05\n",
    "    vatt = np.vectorize(attributable)\n",
    "    vconf = np.vectorize(confidence)\n",
    "    \n",
    "    p = vatt(x, y, tau)\n",
    "    pL, pU = vconf(x, y, alpha, tau)\n",
    "    \n",
    "    fig, ax = plt.subplots(1, figsize=(8,4))\n",
    "    \n",
    "    ax.plot(tau, p, c='red', label='Probability')\n",
    "    if x == 0.:\n",
    "        ax.fill_between(tau, pL, 1, color='red', alpha=0.25)\n",
    "    else:\n",
    "        ax.fill_between(tau, pL, pU, color='red', alpha=0.25)\n",
    "    ax.set_xlabel(r\"Ratio of time periods ($\\tau$)\")\n",
    "    ax.set_ylabel(\"Probability\")\n",
    "    ax.set_xlim((tau.min(), tau.max()))\n",
    "    ax.set_ylim((0, 1))\n",
    "    plt.show()\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.html.widgets import interact, FloatSlider\n",
    "x = FloatSlider(min=0,max=100, step=1, value=10, description=\"Pre-climate change events\")\n",
    "y = FloatSlider(min=1,max=200, step=1, value=10, description=\"Post-climate change events\")\n",
    "interact(plot_results, x=x, y=y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}