Fitting a distribution to radii of maximum winds v2.0.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting a distribution to radius to maximum wind observations\n",
    "\n",
    "The historical record has only sparse observations of radius to maximum winds $R_{max}$ in the Australian region (2002 onwards). As in Vickery _et al._ (2000), we assume that $R_{max}$ fits a log-normal distribution. Powell _et al._ (2005) provide a functional form for the distribution (their Eq. 7), and we will use this as a first estimate. The resulting model is intended for application in setting $R_{max}$ values for stochastically generated storms in TCRM.\n",
    "\n",
    "Note that this model describes the log normal distribution of $R_{max}$ in kilometres -- Powell _et al._ define their model in nautical miles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from __future__ import division, print_function\n",
    "import os\n",
    "from os.path import join as pjoin\n",
    "from matplotlib import pyplot as plt\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "from Utilities.metutils import convert\n",
    "\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "\n",
    "import pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.nonparametric.api as smnp\n",
    "from six import string_types\n",
    "    \n",
    "from statsmodels.tools.tools import ECDF\n",
    "\n",
    "from lmfit import Model, Minimizer, fit_report, conf_interval, printfuncs, report_fit\n",
    "import corner\n",
    "\n",
    "import seaborn as sns\n",
    "from seaborn.utils import _kde_support\n",
    "sns.set_style(\"darkgrid\")\n",
    "sns.set_context(\"poster\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First a short function to convert the formatted latitude/longitude values to actual numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def convertLatLon(strval):\n",
    "    \"\"\"\n",
    "    Convert a string representing lat/lon values from '140S to -14.0, etc.\n",
    "    \n",
    "    :param str strval: string containing the latitude or longitude.\n",
    "    \n",
    "    :returns: Latitude/longitude as a float value.\n",
    "    \n",
    "    \"\"\"\n",
    "    hemi = strval[-1].upper()\n",
    "    fval = float(strval[:-1]) / 10.\n",
    "    if (hemi == 'S') | (hemi == 'W'): \n",
    "        fval *= -1\n",
    "    if (hemi == 'E') | (hemi == 'W'):\n",
    "        fval = fval % 360\n",
    "    return fval\n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the data structure and a small function to load a file. This uses the JTWC data format, described [here](http://www.usno.navy.mil/NOOC/nmfc-ph/RSS/jtwc/best_tracks/shindex.php). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "COLNAMES = ['BASIN','Number', 'Datetime','TECHNUM', 'TECH','TAU', 'Latitude', 'Longitude', 'Windspeed','Pressure',\n",
    "            'Status', 'RAD', 'WINDCODE','RAD1', 'RAD2','RAD3', 'RAD4','Poci', 'Roci','rMax', 'GUSTS','EYE',\n",
    "            'SUBREGION','MAXSEAS', 'INITIALS','DIR', 'SPEED','STORMNAME', 'DEPTH','SEAS',\n",
    "            'SEASCODE','SEAS1', 'SEAS2','SEAS3', 'SEAS4'] \n",
    "\n",
    "COLTYPES = ['|S2', 'i', datetime, 'i', '|S4', 'i', 'f', 'f', 'f', 'f', \n",
    "            '|S4', 'f', '|S3', 'f', 'f', 'f', 'f', 'f', 'f', 'f', 'f', 'f',\n",
    "            '|S1', 'f', '|S3', 'f', 'f', '|S10', '|S1', 'f', \n",
    "            '|S3', 'f', 'f', 'f', 'f']\n",
    "COLUNITS = ['', '', '', '', '', '', '', '', 'kts', 'hPa', \n",
    "            '', 'nm', '', 'nm', 'nm', 'nm', 'nm', 'hPa', 'nm', 'nm', 'kts', 'nm',\n",
    "            '', '', '', 'degrees', 'kts', '', '', '',\n",
    "            '', '', '', '', '']\n",
    "DATEFORMAT = \"%Y%m%d%H\"\n",
    "dtype = np.dtype({'names':COLNAMES, 'formats':COLTYPES})\n",
    "converters = {\n",
    "    1: lambda s: s.strip(' ,'),\n",
    "    2: lambda s: datetime.strptime(s.strip(' ,'), DATEFORMAT),\n",
    "    6: lambda s: float(convertLatLon(s.strip(' ,'))),\n",
    "    7: lambda s: float(convertLatLon(s.strip(' ,'))),\n",
    "    8: lambda s: s.strip(' ,'),\n",
    "    9: lambda s: s.strip(' ,'),\n",
    "    10: lambda s: s.strip(' ,'),\n",
    "    11: lambda s: convert(float(s.strip(' ,') or 0), COLUNITS[11], 'km'),\n",
    "    12: lambda s: s.strip(' ,'),\n",
    "    13: lambda s: convert(float(s.strip(' ,') or 0), COLUNITS[13], 'km'),\n",
    "    14: lambda s: convert(float(s.strip(' ,') or 0), COLUNITS[14], 'km'),\n",
    "    15: lambda s: convert(float(s.strip(' ,') or 0), COLUNITS[15], 'km'),\n",
    "    16: lambda s: convert(float(s.strip(' ,') or 0), COLUNITS[16], 'km'),\n",
    "    17: lambda s: float(s.strip(' ,')),\n",
    "    18: lambda s: convert(float(s.strip(' ,') or 0), COLUNITS[18], 'km'),\n",
    "#    19: lambda s: float(s.strip(' ,'))\n",
    "    19: lambda s: convert(float(s.strip(' ,') or 0), COLUNITS[19], 'km')\n",
    "}\n",
    "delimiter = (3,4,12,4,6,5,7,7,5,6,4,5,5,6,6,6,6,6,6,5,5,5,5)\n",
    "skip_header = 0\n",
    "usecols = tuple(range(23))\n",
    "missing_value = \"\"\n",
    "filling_values = 0\n",
    "\n",
    "def loadData(filename):\n",
    "    try:\n",
    "        data = np.genfromtxt(filename, dtype, delimiter=delimiter, skip_header=skip_header, \n",
    "                             converters=converters, missing_values=missing_value, \n",
    "                             filling_values=filling_values, usecols=usecols, autostrip=True, invalid_raise=False)\n",
    "    except IndexError:\n",
    "        try:\n",
    "            data = np.genfromtxt(filename, dtype, delimiter=delimiter, skip_header=skip_header, \n",
    "                             converters=converters, missing_values=missing_value, \n",
    "                             filling_values=filling_values, usecols=tuple(range(18)), autostrip=True, invalid_raise=False)\n",
    "        except IndexError:\n",
    "            data = np.genfromtxt(filename, dtype, delimiter=[3,4,12,4,6,5,7,7,5], skip_header=skip_header, \n",
    "                             converters=converters, missing_values=missing_value, \n",
    "                             filling_values=filling_values, usecols=tuple(range(9)), autostrip=True, invalid_raise=False)\n",
    "    return data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often the b-deck files contain multiple records with the same time stamp. This is to record information on different wind speed radii (e.g. the radius to 34-knot winds, 48-knot winds, etc.). We can quickly filter out this extra information using [`numpy.unique()`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.unique.html). Additional filtering restricts to a known domain and only those storms that are of Tropical Storm or Typhoon strength, those that have radius to maximum wind records, and those that include a pressure for the outermost closed isobar ($p_{oci}$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def filterData(data):\n",
    "    datetimes, idx = np.unique(data['Datetime'], True)\n",
    "    filter1 = (data['Status'][idx] == 'TS') | (data['Status'][idx] == 'TY')\n",
    "    filter2 = (data['Longitude'][idx] >= 60.) & (data['Longitude'][idx] <= 180.)\n",
    "    filter3 = (data['rMax'][idx] >= 0.1)\n",
    "    filter4 = (data['Poci'][idx] > 0.1)\n",
    "    subsidx = np.nonzero(filter1 & filter2 & filter3 & filter4)\n",
    "    return data[subsidx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now churn through the best-track files (unmodified) and pull out $R_{max}$, $p_c$, $p_{oci}$ and latitude values. This assumes you have the JTWC best track files somewhere locally - no download performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def processfiles(path, basin):\n",
    "    rmax = np.array([])\n",
    "    prs = np.array([])\n",
    "    lat = np.array([])\n",
    "    poci = np.array([])\n",
    "    for root, dirs, files in os.walk(path):\n",
    "        if root.endswith(basin):\n",
    "            for file in files:\n",
    "                data = loadData(pjoin(root, file))\n",
    "                if 'Status' in data.dtype.names:\n",
    "                    data = filterData(data)\n",
    "                    if 'rMax' in data.dtype.names:\n",
    "                        rmax = np.append(rmax, data['rMax'])\n",
    "                        prs = np.append(prs, data['Pressure'])\n",
    "                        poci = np.append(poci, data['Poci'])\n",
    "                        lat = np.append(lat, data['Latitude'])\n",
    "    return rmax, prs, poci, lat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "inputPath = \"C:\\\\WorkSpace\\\\data\\\\Raw\\\\best_tracks\"\n",
    "rmax, prs, poci, lat = processfiles(inputPath, 'sh')\n",
    "#outputFile = pjoin(inputPath, \"rmax-sh.csv\")\n",
    "#np.savetxt(outputFile, np.column_stack((rmax, prs, poci, lat)), delimiter=',', fmt='%6.1f')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we test the first hypothesis - that the distribution of $R_{max}$ is represented by a log-normal distribution. Plot up the probability distribution function, with a kernel density estimate and a fitted log-normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"Parameter estimates:       Shape; Location (fixed);    Scale;    Mean\")\n",
    "fig, ax = plt.subplots(1,1)\n",
    "sns.distplot(rmax, bins=np.arange(0, 101, 10),\n",
    "             kde_kws={'clip':(0, 100), 'label':\"KDE\"}, ax=ax)\n",
    "\n",
    "shape, loc, scale = stats.lognorm.fit(rmax, scale=np.mean(rmax), floc=0)\n",
    "print(\"Southern hemisphere basin: \", shape, loc, scale, np.mean(rmax))\n",
    "x = np.arange(1, 201)\n",
    "v = stats.lognorm.pdf(x, shape, loc=loc, scale=scale)\n",
    "fcdf = stats.lognorm.cdf(np.sort(rmax), shape, loc=loc, scale=scale)\n",
    "\n",
    "ax.plot(x, v, label=\"Lognormal fit\")\n",
    "ax.legend(loc=0)\n",
    "ax.set_xlabel(r'$R_{max}$ (km)')\n",
    "ax.set_ylabel('Probability')\n",
    "ax.set_xlim((0, 100))\n",
    "ax.set_title(\"Southern hemisphere (2002-2014)\")\n",
    "\n",
    "\n",
    "fig.tight_layout()\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the empirical CDF to the fitted CDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ecdf = ECDF(rmax, side='left')\n",
    "\n",
    "plt.plot(np.sort(rmax), ecdf.y[1:])\n",
    "plt.plot(np.sort(rmax), fcdf, 'r' )\n",
    "rsq = stats.pearsonr(np.sort(rmax), fcdf)[0]**2\n",
    "plt.text( 10, 0.9, r\"$R^{2}$ = %f\"%rsq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting to multiple parameters\n",
    "\n",
    "In this approach, the natural logarithm of $R_{max}$ is modelled as follows:\n",
    "\n",
    "$\\ln R_{max} = \\alpha + \\beta \\Delta p + \\gamma \\exp^{(-\\delta \\Delta p^2)} + \\zeta |\\lambda| + \\varepsilon$\n",
    "\n",
    "The constants are determined by a generalised linear model, $\\Delta p$ is the central pressure deficit (hPa), $\\lambda$ the latitude and $\\varepsilon$ is a normal random variable with zero mean. We choose an exponential decay function to ensure physically realistic behaviour at large $\\Delta p$. Other models examined were quadratic in $\\Delta p$, which produced unrealistic behaviour when $\\Delta p$ increased beyond 100 hPa.\n",
    "\n",
    "Additional filtering is needed here to remove records where the pressure of the outermost closed isobar ($p_{oci}$) is not known."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def filterPoci(field, poci):\n",
    "    filter1 = (poci >= 0.1)\n",
    "    subsidx = np.nonzero(filter1)\n",
    "    return field[subsidx]\n",
    "\n",
    "rmax = filterPoci(rmax, poci)\n",
    "dp = filterPoci(poci, poci) - filterPoci(prs, poci)\n",
    "dp = np.extract(np.nonzero(rmax), dp)\n",
    "dpsq = dp*dp\n",
    "expdp = np.exp(-dp)\n",
    "expdpsq = np.exp(-dpsq)\n",
    "lat = filterPoci(lat, poci)\n",
    "lat = np.extract(np.nonzero(rmax), lat)\n",
    "rmax = np.extract(np.nonzero(rmax), rmax)\n",
    "\n",
    "latsq = lat*lat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now fit a model, based on the functional form given above. Start by defining the independent variables ($\\Delta p$, $\\lambda$) and dependent variable (in this case $\\ln(R_{max})$). This is followed by defining individual components of the hypothesised model, namely a linear component in $\\Delta p$, a linear component in $|\\lambda|$ and an exponential decay term for $\\Delta p^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X = np.column_stack((dp, lat))\n",
    "y = np.log(rmax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def exp_dpsq(x, gamma, delta):\n",
    "    dp = x[:,0]\n",
    "    return gamma*np.exp(-delta*dp*dp)\n",
    "\n",
    "def lin_dp(x, alpha, beta):\n",
    "    dp = x[:,0]\n",
    "    return alpha + beta*dp\n",
    "\n",
    "def lin_lat(x, zeta):\n",
    "    lat = np.abs(x[:,1])\n",
    "    return zeta*lat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct the model using the `lmfit` `Model` objects, define the model parameters (i.e. the coefficients), and build a fuction that returns the residual of the model fit for a given solution, which will be used in a minimizer fuction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rmod = Model(lin_dp) + Model(exp_dpsq) + Model(lin_lat)\n",
    "params = rmod.make_params(alpha=1., beta=-0.001, gamma=.1, delta=.001, zeta=.001)\n",
    "def resid(p):\n",
    "    return p['alpha'] + p['beta']*X[:,0] + p['gamma']*np.exp(-p['delta']*X[:,0]*X[:,0]) + p['zeta']*np.abs(X[:,1]) - y\n",
    "\n",
    "mini = Minimizer(resid, params)\n",
    "result = mini.minimize()\n",
    "print(fit_report(result.params))\n",
    "ci = conf_interval(mini, result)\n",
    "printfuncs.report_ci(ci)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can examine the range of solutions that are supported by the observations, using Markov Chain Monte Carlo sampling of the posterior probability distribution. We can use this to obtain the credible intervals for the parameters (the Bayesian equivalent of confidence intervals)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rr = mini.emcee(burn=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ll = [r'$\\{0}$'.format(v) for v in rr.var_names]\n",
    "with sns.plotting_context(\"notebook\"):\n",
    "    corner.corner(rr.flatchain, labels=ll, truths=list(rr.params.valuesdict().values()),\n",
    "              no_fill_contours=True, fill_contours=False, plot_density=False,\n",
    "              quantiles=[0.05, 0.5, 0.95],\n",
    "                 data_kwargs=dict(color='r', alpha=0.01),\n",
    "             contour_kwargs=dict(color='g'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(report_fit(rr.params))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(rr.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(r'alpha = {0}'.format(rr.params['alpha'].value))\n",
    "print(r'beta = {0}'.format(rr.params['beta'].value))\n",
    "print(r'gamma = {0}'.format(rr.params['gamma'].value))\n",
    "print(r'delta = {0}'.format(rr.params['delta'].value))\n",
    "print(r'zeta = {0}'.format(rr.params['zeta'].value))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "result = rmod.fit(y, x=X, params=params)\n",
    "print(result.fit_report())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "result.params = rr.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(X[:,0], y,         'bo')\n",
    "#plt.plot(X[:,0], result.init_fit, 'k--')\n",
    "plt.plot(X[:,0], result.best_fit, 'r.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can inspect the individual model components of the model and their contribution to the solution. Using the `eval_components` method from the `result` object allows us to inspect each model component separately. If one wanted, you could plot them to visualise the relationships."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "comps = result.eval_components()\n",
    "print(comps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now examine the residuals to formulate a method of describing the random innovations required in a stochastic modelling framework. We assume the residuals are normally distributed and test for normality to confirm our hypothesis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, (ax0, ax1) = plt.subplots(1, 2)\n",
    "\n",
    "ax = sns.distplot(result.residual, kde_kws={'label':'Residuals', 'linestyle':'--'}, ax=ax0, norm_hist=True)\n",
    "pp = sm.ProbPlot(result.residual, stats.norm, fit=True)\n",
    "pp.qqplot('Normal', 'Residuals', line='45', ax=ax1, color='gray',alpha=0.5)\n",
    "fig.tight_layout()\n",
    "x = np.linspace(-2, 2, 1000)\n",
    "\n",
    "ax0.legend(loc=0)\n",
    "\n",
    "fp = stats.norm.fit(result.residual)\n",
    "ax0.plot(x, stats.norm.pdf(x, fp[0], fp[1]), label='Normal', color='r')\n",
    "print(fp)\n",
    "print(stats.mstats.normaltest(result.residual))\n",
    "ax0.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The p-value indicates the normal distribution is a good description of the residuals of the model. The quantile-quantile plot (right) shows the model residuals are slightly fat-tailed for higher values, but are close to normal within 2 standard deviations of the mean. We can then use normally distributed values, with variancce of 0.33, to inform our stochastic model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functional form\n",
    "The basic functional form of the model can readily be plotted, using the `eval` method on the `result` object. All that is needed are input values representing the pressure deficit and latitude values we want to explore. Here, we calculate the mean $R_{max}$ value for $\\Delta p$ values between 0 and 100 hPa, and latitude ranging from 30$^{\\circ}$S to 2$^{\\circ}$S. This allows us to demonstrate the influence of $\\Delta p$ and $\\lambda$ on $R_{max}$, without the random variations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "deltap = np.linspace(0, 100, 100)\n",
    "lats = np.arange(-30, -1, 4)\n",
    "#lats = np.arange(2, 31, 4)\n",
    "fig, ax = plt.subplots(1,1)\n",
    "sns.set_palette(\"RdBu\", 10)\n",
    "for l in lats:\n",
    "    xx = np.column_stack((deltap, l*np.ones(len(deltap))))\n",
    "    yy = result.eval(x=xx)\n",
    "    ax.plot(deltap, np.exp(yy), label=\"%d\"%l)\n",
    "ax.set_ylabel(r\"$R_{max}$ (km)\")\n",
    "ax.set_xlabel(r\"$\\Delta p$ (hPa)\")\n",
    "ax.legend(loc=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the mean, $R_{max}$ starts out around 80 km for weak storms and rapidly decreases as intensity increases to $\\Delta p$ = 40 hPa. At this point, $R_{max}$ is close to 30 km, and then steadily reduces for increasing intensity. Lower latitude storms are, in the mean, smaller than at high latitudes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Complete model\n",
    "We construct the complete model by combining the functional form with a random innovation which is based on the residuals from the fitted model. In this case, the random innovations are sampled from a normal distribution with zero mean and $\\sigma$ = 0.33.\n",
    "\n",
    "Using observed values of $\\Delta p$ and $\\lambda$, we can then generate a plot of observed and simulated $R_{max}$ values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nx = len(dp)\n",
    "ind = np.random.choice(np.arange(nx), nx, replace=True)\n",
    "dp0 = dp[ind]\n",
    "l0 = lat[ind]\n",
    "\n",
    "xx = np.column_stack((dp0, l0))\n",
    "yy = result.eval(x=xx) + np.random.normal(scale=0.33, size=nx)\n",
    "\n",
    "\n",
    "rm = np.exp(yy)\n",
    "fig, ax = plt.subplots(1, 2, sharey=True)\n",
    "ax[0].scatter(dp0, rm, c=np.abs(l0), cmap=sns.light_palette('blue', as_cmap=True), s=40, label='Model', alpha=0.5)\n",
    "ax[0].scatter(dp, rmax, c='w', edgecolor='r', s=50, marker='+', label='Observations')\n",
    "ax[0].set_xlim(0, 100)\n",
    "ax[0].set_xlabel(r\"$\\Delta p$ (hPa)\")\n",
    "ax[0].set_ylabel(r\"$R_{max}$ (km)\")\n",
    "ax[0].set_ylim(0, 200)\n",
    "ax[0].legend(loc=1)\n",
    "ax[0].grid(True)\n",
    "\n",
    "ax[1].scatter(l0, rm, c=np.abs(l0), cmap=sns.light_palette('blue', as_cmap=True), s=40, label='Model', alpha=0.5)\n",
    "ax[1].scatter(lat, rmax, c='w', edgecolor='r', s=50, marker='+', label='Observations')\n",
    "ax[1].set_xlim(-30, 0)\n",
    "ax[1].set_xlabel(r\"Latitude\")\n",
    "ax[1].set_ylim(0, 200)\n",
    "ax[1].legend(loc=1)\n",
    "ax[1].grid(True)\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def bivariate_kde(x, y, bw='scott', gridsize=100, cut=3, clip=None):\n",
    "    if isinstance(bw, string_types):\n",
    "        bw_func = getattr(smnp.bandwidths, \"bw_\" + bw)\n",
    "        x_bw = bw_func(x)\n",
    "        y_bw = bw_func(y)\n",
    "        bw = [x_bw, y_bw]\n",
    "    elif np.isscalar(bw):\n",
    "        bw = [bw, bw]\n",
    "\n",
    "    if isinstance(x, pd.Series):\n",
    "        x = x.values\n",
    "    if isinstance(y, pd.Series):\n",
    "        y = y.values\n",
    "\n",
    "    kde = smnp.KDEMultivariate([x, y], \"cc\", bw)\n",
    "    x_support = _kde_support(x, kde.bw[0], gridsize, cut, [x.min(), x.max()])# clip[0])\n",
    "    y_support = _kde_support(y, kde.bw[1], gridsize, cut, [y.min(), y.max()])#clip[1])\n",
    "    xx, yy = np.meshgrid(x_support, y_support)\n",
    "    z = kde.pdf([xx.ravel(), yy.ravel()]).reshape(xx.shape)\n",
    "    return xx, yy, z\n",
    "\n",
    "def l2score(obs, model):\n",
    "    return np.linalg.norm(obs - model, np.inf)\n",
    "\n",
    "xx, yy, odp_rmax = bivariate_kde(dp,  rmax, bw='scott')\n",
    "xx, yy, mdp_rmax = bivariate_kde(dp0, rm, bw='scott')\n",
    "\n",
    "xx, yy, olat_rmax = bivariate_kde(lat,  rmax, bw='scott')\n",
    "xx, yy, mlat_rmax = bivariate_kde(l0, rm, bw='scott')\n",
    "\n",
    "\n",
    "l2rmdp = l2score(odp_rmax, mdp_rmax)\n",
    "l2rmlat = l2score(olat_rmax, mlat_rmax)\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "levs = np.arange(0.01, 0.11, 0.01)\n",
    "ax = sns.kdeplot(dp, rmax, cmap='Reds', kwargs={'levels':levs}, shade=True, shade_lowest=False)\n",
    "ax = sns.kdeplot(dp0, rm, cmap='Blues', kwargs={'levels':levs})\n",
    "ax.set_xlim(0, 100)\n",
    "ax.set_xlabel(r\"$\\Delta p$ (hPa)\")\n",
    "ax.set_ylabel(r\"$R_{max}$ (nm)\")\n",
    "ax.set_ylim(0, 100)\n",
    "ax.grid(True)\n",
    "\n",
    "red = sns.color_palette(\"Reds\")[-2]\n",
    "blue = sns.color_palette(\"Blues\")[-2]\n",
    "ax.text(80, 90, \"Observed\", color=red)\n",
    "ax.text(80, 80, \"Model\", color=blue)\n",
    "ax.text(80, 70, r\"$l_2=${0:.3f}\".format(l2rmdp))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the model reproduces the observations reasonably well on visual inspection. Modelled values of $R_{max}$ are generally less than 50 km for the most intense storms, while for weak storms, $R_{max}$ values tend to be higher, with maximum values occuring for those storms with $\\Delta p < 20$ hPa. \n",
    "\n",
    "The overall distribution of $R_{max}$ is also well reproduced. Here, we present the distribution from the $R_{max}$ model with the fitted log-normal distribution for the observations. There's a slight over-representation of smaller storms ($R_{max} < 25$ km), but above this the distributions match well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "ax = sns.kdeplot(lat, rmax, cmap='Reds', kwargs={'levels':levs}, shade=True, shade_lowest=False)\n",
    "ax = sns.kdeplot(l0, rm, cmap='Blues', kwargs={'levels':levs})\n",
    "ax.set_xlim(-30, 0)\n",
    "ax.set_xlabel(\"Latitude\")\n",
    "ax.set_ylabel(r\"$R_{max}$ (nm)\")\n",
    "ax.set_ylim(0, 100)\n",
    "ax.grid(True)\n",
    "\n",
    "\n",
    "\n",
    "ax.text(-5, 90, \"Observed\", color=red)\n",
    "ax.text(-5, 80, \"Model\", color=blue)\n",
    "ax.text(-5, 70, r\"$l_2=${0:.3f}\".format(l2rmlat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.arange(1, 101)\n",
    "v = stats.lognorm.pdf(x, shape, loc=loc, scale=scale)\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "sns.distplot(rm, bins=np.arange(0, 101, 5),\n",
    "             kde_kws={'clip':(0, 100), 'label':\"Model data (KDE)\"},)\n",
    "ax.plot(x, v, label=\"Lognormal fit from observations\", color='r')\n",
    "ax.legend(loc=0)\n",
    "ax.set_xlabel(r'$R_{max}$ (km)')\n",
    "ax.set_xlim((0, 100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1)\n",
    "sns.distplot(rmax, bins=np.arange(0, 151, 10),\n",
    "             kde_kws={'clip':(0, 150), 'label':\"Observations\"}, ax=ax, \n",
    "             hist_kws={ \"linewidth\":3})\n",
    "sns.distplot(rm, bins=np.arange(0, 151, 10),\n",
    "             kde_kws={'clip':(0, 150), 'label':\"Model\"}, ax=ax, color='r',\n",
    "             hist_kws={ \"linewidth\":3})\n",
    "ax.set_ylabel(\"Probability\")\n",
    "ax.set_xlabel(r\"$R_{max}$ (km)\")\n",
    "ax.set_xlim((0, 120))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "A cursory examination of the AIC scores from these three models indicates the Powell _et al._ (2005) model provides the best balance between magnitude of residuals and degrees of freedom. \n",
    "\n",
    "Heuristically, in the limit of high intensity the simplest model will lead to $\\lim_{\\Delta p\\to\\infty}R_{max} = 0$, which is physically unrealistic. The Wang and Rosowsky model is comparatively less reliable than the Powell _et al._ model across the observed range. And in fact, the OLS model performs better than the corresponding GLS model, using the same form as Powell _et al._ This latter form will have increasing $R_{max}$ for large $\\Delta p$\n",
    "\n",
    "So our final model (based on southern hemisphere observations) is:\n",
    "\n",
    "$\\ln R_{max} = 4.4726 -0.04057 \\Delta p + 0.000313182 \\Delta p^2 + 0.0001455 \\lambda^2 + \\varepsilon$\n",
    "\n",
    "where $\\varepsilon = \\mathcal{N}(0, 0.353)$.\n",
    "\n",
    "<a id='references'></a>\n",
    "## References\n",
    "\n",
    "1. Powell, M., Soukup, G., Cocke, S., Gulati, S., Morisseau-Leroy, N., Hamid, S., Dorst, N. and Axe, L. (2005): State of Florida hurricane loss projection model: Atmospheric science component. _Journal of Wind Engineering and Industrial Aerodynamics_, __93__, pp 651-674.\n",
    "2. Vickery, P. J., Skerlj, P. F. and Twisdale, L. A. (2000): Simulation of Hurricane Risk in the U.S. Using Empirical Track Model. _Journal of Structural Engineering_, __126__, pp 1222-1237.\n",
    "3. Wang, Y. and Rosowsky, D. V. (2012): Joint distribution model for prediction of hurricane wind speed and size. _Structural Safety_, __35__, pp 40-51."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dparray = np.arange(10, 51, 5)\n",
    "latarray = np.arange(-23, -5, 2)\n",
    "testinput = np.column_stack((dparray, latarray))\n",
    "np.random.seed(10)\n",
    "print(result.params)\n",
    "yy = result.eval(x=testinput) #+ np.random.normal(scale=0.335)\n",
    "print(np.exp(yy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "eps = 0#np.random.normal(0, scale=0.335)\n",
    "rmw = np.exp(result.eval(x=np.column_stack((np.array([25]), np.array([-15]))))+eps)\n",
    "print(rmw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "eps = 0#np.random.normal(0, scale=0.335)\n",
    "newparams = rmod.make_params(alpha=3.5, beta=-0.004, gamma=.7, delta=.002, zeta=.001)\n",
    "yy = result.eval(params=newparams, x=testinput)+eps\n",
    "rmw = np.exp(yy)\n",
    "print(rmw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.scatter(dp, np.abs(lat), rmax, c=rmax, cmap=sns.light_palette('blue', as_cmap=True))\n",
    "ax.set_xlabel(r\"$\\Delta p$ (hPa)\")\n",
    "ax.set_ylabel(r\"$|\\lambda|$ (degrees)\")\n",
    "ax.set_zlabel(r\"$R_{max}$ (km)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xbins = np.arange(0, 101, 1)\n",
    "ybins = np.arange(0, 30, 1)\n",
    "(count, xedge, yedge, img) = plt.hist2d(dp, np.abs(lat), bins=[xbins, ybins], weights=rmax, normed=True)\n",
    "plt.xlabel(r\"$\\Delta p$ (hPa)\")\n",
    "plt.ylabel(r\"$|\\lambda|$ (degrees)\")\n",
    "plt.colorbar(label=r\"$R_{max}$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from thinkbayes import Pmf, Suite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#rm = filterPoci(rmax, poci)\n",
    "print(len(rmax))\n",
    "pmf = Pmf()\n",
    "for r in rmax:\n",
    "    pmf.Incr(r, 1)\n",
    "    \n",
    "pmf.Normalize()\n",
    "for x,y in pmf.Items():\n",
    "    print( x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def log_prior(theta):\n",
    "    alpha, beta, gamma, zeta, sigma = theta\n",
    "    if sigma < 0:\n",
    "        return -np.inf  # log(0)\n",
    "    else:\n",
    "        return -1.5 * np.log(1 + beta ** 2 + gamma ** 2 + zeta ** 2) - np.log(sigma)\n",
    "\n",
    "def log_likelihood(theta, x1, x2, x3, y):\n",
    "    alpha, beta, gamma, zeta, sigma = theta\n",
    "    y_model = alpha + beta * x1 + gamma * x2 + zeta * x3\n",
    "    return -0.5 * np.sum(np.log(2 * np.pi * sigma ** 2) + (y - y_model) ** 2 / sigma ** 2)\n",
    "\n",
    "def log_posterior(theta, x1, x2, x3, y):\n",
    "    return log_prior(theta) + log_likelihood(theta, x1, x2, x3, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ndim=5\n",
    "nwalkers=50\n",
    "nburn=1000\n",
    "nsteps=2000\n",
    "np.random.seed(0)\n",
    "starting_guesses=np.random.random((nwalkers, ndim))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import emcee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ydata = rm\n",
    "xdata1 = dp\n",
    "xdata2 = dp * dp\n",
    "xdata3 = np.abs(lat)**2\n",
    "sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=[xdata1, xdata2, xdata3, ydata])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sampler.run_mcmc(starting_guesses, nsteps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "emcee_trace = sampler.chain[:, nburn:, :].reshape(-1, ndim).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_sigma_level(trace1, trace2, nbins=20):\n",
    "    \"\"\"From a set of traces, bin by number of standard deviations\"\"\"\n",
    "    L, xbins, ybins = np.histogram2d(trace1, trace2, nbins)\n",
    "    L[L == 0] = 1E-16\n",
    "    logL = np.log(L)\n",
    "\n",
    "    shape = L.shape\n",
    "    L = L.ravel()\n",
    "\n",
    "    # obtain the indices to sort and unsort the flattened array\n",
    "    i_sort = np.argsort(L)[::-1]\n",
    "    i_unsort = np.argsort(i_sort)\n",
    "\n",
    "    L_cumsum = L[i_sort].cumsum()\n",
    "    L_cumsum /= L_cumsum[-1]\n",
    "    \n",
    "    xbins = 0.5 * (xbins[1:] + xbins[:-1])\n",
    "    ybins = 0.5 * (ybins[1:] + ybins[:-1])\n",
    "\n",
    "    return xbins, ybins, L_cumsum[i_unsort].reshape(shape)\n",
    "\n",
    "\n",
    "def plot_MCMC_trace(ax, xdata1, xdata2, xdata3, ydata, trace, scatter=False, **kwargs):\n",
    "    \"\"\"Plot traces and contours\"\"\"\n",
    "    xbins, ybins, sigma = compute_sigma_level(trace[0], trace[1])\n",
    "    ax.contour(xbins, ybins, sigma.T, levels=[0.683, 0.955], **kwargs)\n",
    "    if scatter:\n",
    "        ax.plot(trace[0], trace[1], ',k', alpha=0.1)\n",
    "    ax.set_xlabel(r'$\\alpha$')\n",
    "    ax.set_ylabel(r'$\\beta$')\n",
    "    \n",
    "    \n",
    "def plot_MCMC_model(ax, xdata1, xdata2, xdata3, ydata, trace):\n",
    "    \"\"\"Plot the linear model and 2sigma contours\"\"\"\n",
    "    ax.plot(xdata1, ydata, 'ok')\n",
    "\n",
    "    alpha, beta, gamma, zeta = trace[:4]\n",
    "    xfit = np.linspace(-20, 120, 10)\n",
    "    yfit = alpha[:, None] + beta[:, None] * xfit\n",
    "    mu = yfit.mean(0)\n",
    "    sig = 2 * yfit.std(0)\n",
    "\n",
    "    ax.plot(xfit, mu, '-k')\n",
    "    ax.fill_between(xfit, mu - sig, mu + sig, color='lightgray')\n",
    "\n",
    "    ax.set_xlabel(r'$\\Delta p$ (hPa)')\n",
    "    ax.set_ylabel(r'$R_{max}$ (km)')\n",
    "\n",
    "def plot_MCMC_results(xdata1, xdata2, xdata3, ydata, trace, colors='k'):\n",
    "    \"\"\"Plot both the trace and the model together\"\"\"\n",
    "    fig, ax = plt.subplots(1, 2, figsize=(10, 4))\n",
    "    plot_MCMC_trace(ax[0], xdata1, xdata2, xdata3, ydata, trace, True, colors=colors)\n",
    "    plot_MCMC_model(ax[1], xdata1, xdata2, xdata3, ydata, trace)\n",
    "    fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_MCMC_results(xdata1, xdata2, xdata3, ydata, emcee_trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import corner\n",
    "fig = corner.corner(emcee_trace.T, labels=[r\"$\\alpha$\", r\"$\\beta$\", r\"$\\gamma$\", r\"$\\zeta$\", r\"$\\ln\\,f$\"], \n",
    "                    plot_density=False, no_fill_contours=True, data_kwargs={'color':'r'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pymc\n",
    "print(pymc.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "alpha = pymc.Normal('alpha', 0, 1)\n",
    "\n",
    "@pymc.stochastic(observed=False)\n",
    "def beta(value=0):\n",
    "    return -1.5*np.log(1+value**2)\n",
    "\n",
    "@pymc.stochastic(observed=False)\n",
    "def gamma(value=0):\n",
    "    return -1.5*np.log(1+value**2)\n",
    "\n",
    "@pymc.stochastic(observed=False)\n",
    "def zeta(value=0):\n",
    "    return -1.5*np.log(1+value**2)\n",
    "\n",
    "@pymc.stochastic(observed=False)\n",
    "def sigma(value=1):\n",
    "    return -np.log(abs(value))\n",
    "\n",
    "@pymc.deterministic\n",
    "def y_model(x1=xdata1, x2=xdata2, x3=xdata3, alpha=alpha, beta=beta, gamma=gamma, zeta=zeta):\n",
    "    return alpha + beta * x1 + gamma * x2 + zeta * x3\n",
    "\n",
    "y = pymc.Normal('y', mu=y_model, tau=1. / sigma ** 2, observed=True, value=ydata)\n",
    "\n",
    "model1 = dict(alpha=alpha, beta=beta, gamma=gamma, zeta=zeta, sigma=sigma, y_model=y_model, y=y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "S = pymc.MCMC(model1)\n",
    "S.sample(iter=10000,burn=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pymc_trace = [S.trace('alpha')[:],\n",
    "              S.trace('beta')[:],\n",
    "              S.trace('gamma')[:],\n",
    "              S.trace('zeta')[:],\n",
    "              S.trace('sigma')[:]]\n",
    "plot_MCMC_results(xdata1, xdata2, xdata3, ydata, pymc_trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = corner.corner(np.array(pymc_trace).T, labels=[r\"$\\alpha$\", r\"$\\beta$\", r\"$\\gamma$\", r\"$\\zeta$\", r\"$\\ln\\,f$\"], \n",
    "                    plot_density=False, no_fill_contours=True, data_kwargs={'color':'r'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "def func(x, a, b, c, d, f):\n",
    "    dp = x[:,0]\n",
    "    lat = x[:,1]\n",
    "    return a + b*dp + c*np.exp(-d*dp*dp) + f*np.abs(lat)\n",
    "\n",
    "xx = np.column_stack((dp, lat))\n",
    "\n",
    "popt, pcov = curve_fit(func, xx, np.log(rmax))\n",
    "perr = np.sqrt(np.diag(pcov))\n",
    "print(popt)\n",
    "print(perr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nx = len(dp)\n",
    "ind = np.random.choice(np.arange(nx), nx, replace=True)\n",
    "dp0 = dp[ind]\n",
    "l0 = lat[ind]\n",
    "xx = np.column_stack((dp0, l0))\n",
    "yy = func(xx, *popt) + np.random.normal(scale=0.3, size=nx)\n",
    "rm = np.exp(yy)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.scatter(dp0, rm, c=np.abs(l0), cmap=sns.light_palette('blue', as_cmap=True), s=40, label='Model', alpha=0.5)\n",
    "ax.scatter(dp, rmax, c='w', edgecolor='r', s=50, marker='+', label='Observations')\n",
    "ax.set_xlim(0, 100)\n",
    "ax.set_xlabel(r\"$\\Delta p$ (hPa)\")\n",
    "ax.set_ylabel(r\"$R_{max}$ (km)\")\n",
    "ax.set_ylim(0, 100)\n",
    "ax.legend(loc=1)\n",
    "ax.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "deltap=np.linspace(0, 200, 200)\n",
    "lats = np.arange(-30, -1, 4)\n",
    "\n",
    "fig, ax = plt.subplots(1,1)\n",
    "sns.set_palette(sns.color_palette(\"coolwarm\", 8))\n",
    "for l in lats:\n",
    "    xx = np.column_stack((deltap, l*np.ones(len(deltap))))\n",
    "    yy = func(xx, *popt)\n",
    "    ax.plot(deltap, np.exp(yy), label=\"%d\"%l)\n",
    "    \n",
    "ax.set_ylabel(r\"$R_{max}$ (km)\")\n",
    "ax.set_xlabel(r\"$\\Delta p$ (hPa)\")\n",
    "ax.legend(loc=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xx = np.column_stack((dp, lat))\n",
    "yy = func(xx, *popt)\n",
    "#rm = filterPoci(rmax, poci)\n",
    "\n",
    "resid = np.log(rmax) - yy\n",
    "print(yy)\n",
    "print(np.log(rmax))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, (ax0, ax1) = plt.subplots(1, 2)\n",
    "\n",
    "ax = sns.distplot(resid, kde_kws={'label':'Residuals', 'linestyle':'--'}, ax=ax0, norm_hist=True)\n",
    "pp = sm.ProbPlot(resid, stats.norm, fit=True)\n",
    "pp.qqplot('Normal', 'Residuals', line='45', ax=ax1, color='gray',alpha=0.5)\n",
    "fig.tight_layout()\n",
    "\n",
    "#ppfit = pp.fit_params\n",
    "\n",
    "x = np.linspace(-2, 2, 1000)\n",
    "\n",
    "ax0.legend(loc=0)\n",
    "\n",
    "fp = stats.norm.fit(resid)\n",
    "ax0.plot(x, stats.norm.pdf(x, fp[0], fp[1]), label='Normal', color='r')\n",
    "print(fp)\n",
    "print(stats.mstats.normaltest(resid))\n",
    "ax0.legend()\n",
    "#p = list(results.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}