Extreme value introduction.ipynb

{
 "metadata": {
  "name": "",
  "signature": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Introduction to extreme value statistics in Python"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The concept fo a return period can be confusing and misleading. Often I've heard people say about a 1 in 100 year event that \"We've just had two events in consecutive years! How can they be 1 in 100 year events?\". My response is that \"the event has a small chance of happening each year. It happens on average once every 100 years, but it could occur in any year\". I prefer to use the term \"annual exceedance probability\" to describe extreme values. This way, the focus is on the probability of an event, not the (expected) time period between occurences."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We first generate some 'observations' to work with. In this case, we simulate them from a random distribution to look like the observations of annual maximum wind speed from a weather station. The generalised extreme value distribution (``scipy.stats.genextreme``) represents the distribution of annual maxima of wind speeds well. If we took the annual maximum wind speed observations from each year from 1901-2000, then it might look something like this."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "from matplotlib import pyplot\n",
      "import numpy as np\n",
      "np.random.seed(1234)\n",
      "from scipy.stats import genextreme\n",
      "\n",
      "def simulate(loc=25., shp=0.27, scale=12.5553):\n",
      "    years = np.arange(1901, 2001)\n",
      "    ranks = np.arange(1, 101)\n",
      "    obs = -1*np.ones(100)\n",
      "    while obs.min() < 5.:\n",
      "        obs = genextreme.rvs(shp, loc, scale, size=100)\n",
      "\n",
      "    pyplot.plot(years, obs, '.')\n",
      "    pyplot.xlabel('Year',fontsize='small')\n",
      "    pyplot.ylabel('Wind speed (m/s)',fontsize='small')\n",
      "    pyplot.ylim((0,100))\n",
      "    pyplot.grid(True)\n",
      "    return obs\n",
      "\n",
      "shp = 0.27\n",
      "loc = 25.\n",
      "scale = 12.5553\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Import widgets for interactive notebook\n",
      "from IPython.html.widgets import interact, fixed\n",
      "from IPython.html import widgets\n",
      "\n",
      "interact(simulate,\n",
      "         loc=widgets.FloatSliderWidget(min=0.,max=50.,step=0.1,value=25.,description=\"loc\"),\n",
      "         shp=widgets.FloatSliderWidget(min=0.,max=1.,step=0.01,value=0.27,description=\"shp\"),\n",
      "         scale=widgets.FloatSliderWidget(min=0.,max=20.,step=0.05,value=12.5,description=\"scale\"))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "See how there are some clusters of our 'observations' where there are very high values in (possibly) consecutive years? \n",
      "If we sort the observations into ascending order, then we see the following pattern: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sortObs = obs.copy()\n",
      "sortObs.sort()\n",
      "\n",
      "pyplot.plot(ranks, sortObs,'.')\n",
      "pyplot.ylabel('Wind speed (m/s)', fontsize='small')\n",
      "pyplot.xlabel('Rank',fontsize='small')\n",
      "pyplot.xticks(np.arange(0,101,10.) )\n",
      "pyplot.ylim((0, 100))\n",
      "pyplot.grid(True)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The histogram of the observations:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n, b, p = pyplot.hist(sortObs, np.arange(0, 100, 5))\n",
      "pyplot.xlabel('Wind speed (m/s)', fontsize='small')\n",
      "pyplot.ylabel('Count', fontsize='small')\n",
      "pyplot.grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pdf = genextreme.pdf(obs, shp, loc, scale)\n",
      "\n",
      "simWind = np.arange(1,150, 0.5)\n",
      "simsf = genextreme.sf(simWind, shp, loc, scale)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The cumulative distribution function describes the probability that an observation $v$ will have a value less than or equal to a given wind speed $V$\n",
      "\n",
      "$$F_V(v) = P(V \\leq v) = \\int_{0}^{v} f_V(t)dt$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cdf = genextreme.cdf(obs, shp, loc, scale)\n",
      "\n",
      "pyplot.plot(obs, cdf,'.')\n",
      "pyplot.xlabel('Wind speed (m/s)',fontsize='small')\n",
      "pyplot.ylabel('Cumulative probability',fontsize='small')\n",
      "pyplot.grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The probability of an observation being greater than a given threshold is given by the survival function. The survival function is also 1 - cdf:\n",
      "\n",
      "$$S_V(v) = P(\\{V > v\\}) = \\int_v^{\\infty} f(t)dt = 1 - F_V(v)$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sf = genextreme.sf(obs, shp, loc, scale)\n",
      "\n",
      "pyplot.plot(obs, sf, '.')\n",
      "pyplot.ylabel('Exceedance probability', fontsize='small')\n",
      "pyplot.xlabel('Wind speed (m/s)',fontsize='small')\n",
      "pyplot.yticks(np.arange(0, 1.01, 0.1))\n",
      "pyplot.grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rp = 1./sf\n",
      "simrp = 1./simsf\n",
      "pyplot.figure()\n",
      "pyplot.semilogx(rp, obs, '.')\n",
      "pyplot.semilogx(simrp, simWind, 'r-')\n",
      "pyplot.xlabel('Return period (years)')\n",
      "pyplot.ylabel('Wind speed (m/s)')\n",
      "pyplot.xlim(1, 1000)\n",
      "pyplot.grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pyplot.plot( obs, sf, '.')\n",
      "pyplot.plot( simWind, simsf, 'r-')\n",
      "pyplot.ylabel('Exceedance probability')\n",
      "pyplot.xlabel('Wind speed (m/s)')\n",
      "pyplot.xlim(0,100)\n",
      "pyplot.grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pyplot.semilogy(obs, sf, '.')\n",
      "pyplot.semilogy(simWind, simsf, 'r-')\n",
      "pyplot.ylabel('Exceedance probability')\n",
      "pyplot.xlabel('Wind speed (m/s)')\n",
      "pyplot.xlim(0,100)\n",
      "pyplot.ylim(0.001,1)\n",
      "pyplot.grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}