{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "## Problem Type\n",
    "\n",
    "Poisson regression is the kind of regression we do when we want to estimate the effect that our explanatory variables have on the dependent variable, which is of type \"count data\". If we're trying to find a linear combination of the explanatory variables, then our Poisson regression is a subset of generalized linear models.\n",
    "\n",
    "It's \"Poisson\" mainly because we use the Poisson distribution to model the likelihood of the dependent variable.\n",
    "\n",
    "What we get out of this type of model is the relative contribution of each explanatory variable to the value of the dependent variable.\n",
    "\n",
    "## Data structure\n",
    "\n",
    "To use it with this model, the data should be structured as such:\n",
    "\n",
    "- Each row is one measurement.\n",
    "- The columns should be:\n",
    "    - One column per explanatory variable.\n",
    "        - Use ordinal data where possible; otherwise, strictly categorical data should be binarized.\n",
    "    - One column for the dependent variable.\n",
    "\n",
    "## Extensions to the model\n",
    "\n",
    "None.\n",
    "\n",
    "## Reporting summarized findings\n",
    "\n",
    "Here are examples of how to summarize the findings.\n",
    "\n",
    "> For every increase in $X_i$, we expect to see an increase in Y by `mean` (95% HPD: [`lower`, `upper`].\n",
    "\n",
    "## Other notes\n",
    "\n",
    "None."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploratory Data Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"../datasets/ship-damage.txt\")\n",
    "# Log10 transform months\n",
    "df[\"months\"] = df[\"months\"].apply(lambda x: np.log10(x))\n",
    "df[\"period_op\"] = df[\"period_op\"] - 1\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x=df[\"months\"], y=df[\"n_damages\"])\n",
    "plt.xlabel(\"months\")\n",
    "plt.ylabel(\"n_damages\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import theano.tensor as tt\n",
    "\n",
    "with pm.Model() as model:\n",
    "    xs = pm.floatX(df[[\"yr_construction\", \"period_op\", \"months\"]].values)\n",
    "    betas = pm.Normal(\"betas\", mu=0, sd=100 ** 2, shape=(3, 1))\n",
    "    n_damages = tt.dot(xs, betas)\n",
    "    n_damages_like = pm.Poisson(\n",
    "        \"likelihood\", mu=np.exp(n_damages), observed=df[\"n_damages\"]\n",
    "    )\n",
    "    trace = pm.sample(draws=2000)  # , start=pm.find_MAP(), step=pm.Metropolis())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.traceplot(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.forestplot(trace, ylabels=[\"yr_construction\", \"period_op\", \"months\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best interpretation of this is that the log10 number of months that a boat has been used is the strongest positive contributor to the number of damages that a ship takes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.summary(trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Posterior Predictive Checks\n",
    "\n",
    "Let's see what the PPC looks like. We will sample 10,000 predicted values for each row in the dataframe, and plot the 95% HPD of the values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:\n",
    "    ppc = pm.sample_ppc(trace, samples=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppc[\"likelihood\"].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the posterior distribution of predictions vs. actual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_true = df[\"n_damages\"].values\n",
    "lower, med, upper = np.percentile(ppc[\"likelihood\"], [2.5, 50, 97.5], axis=0)\n",
    "\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.errorbar(y_true, med, yerr=[lower, upper], fmt=\"o\")\n",
    "\n",
    "\n",
    "def x_eq_y(ax):\n",
    "    xmin, xmax = min(ax.get_xlim()), max(ax.get_xlim())\n",
    "    ymin, ymax = min(ax.get_ylim()), max(ax.get_ylim())\n",
    "\n",
    "    ax.plot([xmin, xmax], [ymin, ymax])\n",
    "    return ax\n",
    "\n",
    "\n",
    "ax = x_eq_y(ax)\n",
    "ax.set_xlabel(\"true values\")\n",
    "ax.set_ylabel(\"predicted values\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bayesian",
   "language": "python",
   "name": "bayesian"
  },
  "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.6"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "156px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}