Bayesian model of TC occurrence with MCMC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pymc\n",
    "import seaborn as sns\n",
    "sns.set_context('poster')\n",
    "sns.set_style('whitegrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"C:/WorkSpace/data/TC/Allstorms.ibtracs_wmo.v03r05.csv\", delim_whitespace=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnt = df[df['Basin'].isin([' SP' ])].groupby('Season')['Serial_Num'].nunique()\n",
    "years = cnt.index\n",
    "y0, y1 = years[0], years[-1]\n",
    "arr = cnt.values\n",
    "plt.figure(figsize=(12,6));\n",
    "plt.plot(years, arr, '-ok');\n",
    "plt.xlim(y0, y1);\n",
    "plt.xlabel(\"Year\");\n",
    "plt.ylabel(\"Number of storms\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "switchpoint = pymc.DiscreteUniform('switchpoint',\n",
    "                                   lower=0, upper=len(arr))\n",
    "early_mean = pymc.Exponential('early_mean', beta=1)\n",
    "late_mean = pymc.Exponential('late_mean', beta=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@pymc.deterministic(plot=False)\n",
    "def rate(s=switchpoint, e=early_mean, l=late_mean):\n",
    "    out = np.empty(len(arr))\n",
    "    out[:s] = e\n",
    "    out[s:] = l\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "storms = pymc.Poisson('storms', mu=rate, value=arr, observed=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = pymc.Model([switchpoint, early_mean, late_mean, rate, storms])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "mcmc = pymc.MCMC(model)\n",
    "mcmc.sample(iter=10000, burn=1000, thin=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,12))\n",
    "plt.subplot(311);\n",
    "plt.plot(mcmc.trace('switchpoint')[:]+y0);\n",
    "plt.ylabel(\"Switch point\"); \n",
    "plt.subplot(312);\n",
    "plt.plot(mcmc.trace('early_mean')[:]);\n",
    "plt.ylabel(\"Early mean\");\n",
    "plt.subplot(313);\n",
    "plt.plot(mcmc.trace('late_mean')[:]);\n",
    "plt.xlabel(\"Iteration\");\n",
    "plt.ylabel(\"Late mean\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,3))\n",
    "plt.subplot(131);\n",
    "plt.hist(mcmc.trace('switchpoint')[:] + y0, 15);\n",
    "plt.xlabel(\"Switch point\")\n",
    "plt.ylabel(\"Distribution\")\n",
    "plt.subplot(132);\n",
    "plt.hist(mcmc.trace('early_mean')[:], 15);\n",
    "plt.xlabel(\"Early mean\");\n",
    "plt.subplot(133);\n",
    "plt.hist(mcmc.trace('late_mean')[:], 15);\n",
    "plt.xlabel(\"Late mean\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "yp = y0 + mcmc.trace('switchpoint')[:].mean()\n",
    "em = mcmc.trace('early_mean')[:].mean()\n",
    "lm = mcmc.trace('late_mean')[:].mean()\n",
    "print((yp, em, lm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6));\n",
    "plt.plot(years, arr, '-ok');\n",
    "plt.axvline(yp, color='k', ls='--');\n",
    "plt.plot([y0, yp], [em, em], '-b', lw=3);\n",
    "plt.plot([yp, y1], [lm, lm], '-r', lw=3);\n",
    "plt.xlim(y0, y1);\n",
    "plt.xlabel(\"Year\");\n",
    "plt.ylabel(\"Number of storms\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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
}