{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting distributions to data with `paramnormal`.\n",
    "\n",
    "In addition to explicitly creating distributions from known parameters, `paramnormal.[dist].fit` provides a similar, interface to `scipy.stats` maximum-likelihood estimatation methods.\n",
    "\n",
    "Again, we'll demonstrate with a lognormal distribution and compare parameter estimatation with scipy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.simplefilter('ignore')\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn\n",
    "\n",
    "import paramnormal\n",
    "\n",
    "clean_bkgd = {'axes.facecolor':'none', 'figure.facecolor':'none'}\n",
    "seaborn.set(style='ticks', rc=clean_bkgd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by generating a reasonably-sized random dataset and plotting a histogram.\n",
    "\n",
    "The primary method of creating a distribution from named parameters is shown below.\n",
    "\n",
    "The call to `paramnormal.lognornal` translates the parameter to be compatible with scipy. We then chain a call to the `rvs` (random variates) method of the returned scipy distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "x = paramnormal.lognormal(mu=1.75, sigma=0.75).rvs(370)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a histogram to illustrate the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bins = np.logspace(-0.5, 1.75, num=25)\n",
    "fig, ax = plt.subplots()\n",
    "_ = ax.hist(x, bins=bins, density=True)\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlabel('$X$')\n",
    "ax.set_ylabel('Probability')\n",
    "seaborn.despine()\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pretending for a moment that we didn't generate this dataset with explicit distribution parameters, how would we go about estimating them?\n",
    "\n",
    "Scipy provides a maximum-likelihood estimation for estimating parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "print(stats.lognorm.fit(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately those parameters don't really make any sense based on what we know about our articifical dataset.\n",
    "\n",
    "That's where paramnormal comes in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "params = paramnormal.lognormal.fit(x)\n",
    "print(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This matches well with our understanding of the distribution.\n",
    "\n",
    "The returned `params` variable is a `namedtuple` that we can easily use to create a distribution via the `.from_params` methods.  From there, we can create a nice plot of the probability distribution function with our histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dist = paramnormal.lognormal.from_params(params)\n",
    "\n",
    "# theoretical PDF\n",
    "x_hat = np.logspace(-0.5, 1.75, num=100)\n",
    "y_hat = dist.pdf(x_hat)\n",
    "\n",
    "bins = np.logspace(-0.5, 1.75, num=25)\n",
    "fig, ax = plt.subplots()\n",
    "_ = ax.hist(x, bins=bins, density=True, alpha=0.375)\n",
    "ax.plot(x_hat, y_hat, zorder=2, color='g')\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlabel('$X$')\n",
    "ax.set_ylabel('Probability')\n",
    "seaborn.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recap\n",
    "#### Fitting data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "params = paramnormal.lognormal.fit(x)\n",
    "print(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Creating distributions \n",
    "The manual way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "paramnormal.lognormal(mu=1.75, sigma=0.75, offset=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From fit parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "paramnormal.lognormal.from_params(params)"
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}