notebooks/3_Did_it_converge.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adding the data\n",
    "After having tried out various priors, we can now proceed and feed some data to our model.\n",
    "Let's first load the data and write down our model again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pymc3 as pm\n",
    "import arviz as az\n",
    "\n",
    "import sys\n",
    "sys.path.append('../src/')\n",
    "from utils import standardize_area, destandardize_area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use(\"fivethirthyeight\")\n",
    "plt.rcParams['figure.figsize'] = 9, 9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "berlin = pd.read_csv(\"../data/berlin.csv\", index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as mod:\n",
    "    alpha = pm.Normal(\"alpha\", mu=0, sigma=10)\n",
    "    beta = pm.Normal(\"beta\", mu=1, sigma=5)\n",
    "    \n",
    "    sigma = pm.HalfNormal(\"sigma\", sigma=5)\n",
    "    \n",
    "    mu = alpha + beta*berlin[\"livingSpace_s\"]\n",
    "    \n",
    "    rent = pm.Normal(\"rent\", mu=mu, sigma=sigma,\n",
    "                    observed = berlin[\"totalRent_s\"])\n",
    "    \n",
    "    priors = pm.sample_prior_predictive()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We already specified all necessary data in the model but we still need to tell it to start estimating the parameters (in ML speak: train the model).\n",
    "In Bayesian Modelling this process is also called sampling since we're trying to get a probability distributions, the posterior, by sampling from it.\n",
    "The command is thus also called `pm.sample()` in PyMC. The result is a sample from our target probability distribution and is often called the _trace_. This is because the successive sampling traces out the probability distribution. This [visualization](https://chi-feng.github.io/mcmc-demo/app.html#HamiltonianMC,standard) might help to give an intuition what happens in the background. \n",
    "\n",
    "With `draws`, we determine how big a sample we want to have and `tune` determines how much the algorithm should try to find and explore the target probability before starting to sample. Sampling can easily be parallelized by letting the sampler run multiple so called chains. It is recommended to have at least two chains but you can go up to as many cores you have (that's usually four).\n",
    "You will get `draws*chains` number of samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with mod:\n",
    "    trace = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We again feed all the different artificats into the ArviZ InferenceData object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm_data = az.from_pymc3(...)\n",
    "pm_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we first need to check, is if our model actually converged. For this, we will mostly look at two things:\n",
    "- trace plots\n",
    "- different summary statistics\n",
    "\n",
    "Let's start with the trace plots. Again, ArviZ provides a handy function for these:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the left, we see the estimated probability distribution for our parameters. On the right, we see how the samples trace through the distribution. The $y$-axes on the right are the same as the $x$-axes on the left. The $x$-axes on the right are the sample number.\n",
    "\n",
    "\n",
    "So what do we need to look out for when checking these trace plots? In short, any non-random looking patterns. These trace plots actually look all good, so it's probably best to show a few examples that fail. Sometimes the best way to understand something is trying to break it.\n",
    "\n",
    "\n",
    "Go play around with this model and try to make it fail!\n",
    "\n",
    "For most fails, you will also see warnings by PyMC. Another sign that something is wrong with your model is if the sampling process takes very long. This is a relatively simple model with not too much data, so it shouldn't take more than a few seconds."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some hints on what to try:\n",
    "- make the priors really really narrow\n",
    "- use priors with hard boundaries (e.g. Uniform) on a wrong range, e.g. a negative range for beta\n",
    "- use almost no tuning steps, e.g. `tune=10`\n",
    "- use a very small amount of samples, e.g. `draws=100`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as failed_mod:\n",
    "    \n",
    "    ...\n",
    "    \n",
    "    failed_trace = pm.sample(draws=100, tune=100)\n",
    "    \n",
    "failed_data = az.from_pymc3(trace = failed_trace, model = failed_mod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "az.plot_trace(failed_data)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next thing to check after the trace plots are the summary statistics.\n",
    "In particular, there are\n",
    "- ESS, effective sample size, the estimated number of samples\n",
    "- R_hat, a diagnostic that compares the different chains\n",
    "\n",
    "__ESS:__ As we're consecutively sampling from the probability space, sample $n$ is often slightly auto-correlated to sample $n-1$ and sample $n+1$. For this reason, the final sample cannot be treated as a completely independent sample from our target probability distribution, the posterior. ESS estimates how many _effective_ samples we have. If everything worked fine, this should be close to the actual number of samples. If it is much lower then there might be a problem with your model. This model is a very simple one and NUTS, the algorithm used by PyMC3 is very effective, so the effective number of samples is even higher than the actual number of samples! In general, you'll want the ESS to be at least larger than 25% of the actual number of samples. A ESS smaller than that often indicates that your model has some problems with the sampling. Also note that how many effective samples you need depends heavily on the goal of your analysis: If you only want to estimate the mean and a 50% interval, you will already get robust estimates with around 300 effective samples. If you need a high precision and want to estimate a 95% or even 99% interval, you will need much more effective samples to estimate the tails correctly!\n",
    "\n",
    "__R_hat:__ The R_hat diagnostic checks if the different chains converged to the same value. This diagnostic should be very close to 1, definitely smaller than 1.05, some argue even smaller than 1.005, everything larger indicates a problem.\n",
    "\n",
    "\n",
    "We can check these values with ArviZ:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this model, all summary and diagnostic statistics look good: the ESS is even larger than the actual number of samples and R_hat is very close to 1. Good!\n",
    "\n",
    "In the same summary, we can also see the point estimates for the three parameters in the column `mean`. The estimates for alpha and beta should be very close to the coefficients we would get from `sklearn.linear_model.LinearRegression()` (feel free to try this out).\n",
    "\n",
    "HPD is the highest probability density, so that 94% of probability mass is between hpd_3% and hpd_97%.\n",
    "\n",
    "Let's have a look at some visualizations of the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(20, 9))\n",
    "az.plot_forest(pm_data,\n",
    "              # combine the different chains\n",
    "              combined=True,\n",
    "              kind=\"ridgeplot\", ridgeplot_alpha=0.6, \n",
    "              linewidth=1,\n",
    "              ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can visualize the resulting model as before, this time not using the prior but the posterior, that is the trace object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import draw_models\n",
    "\n",
    "...\n",
    "\n",
    "plt.scatter(berlin[\"livingSpace\"], berlin[\"totalRent_s\"]*100, s=4)\n",
    "plt.title(\"Linear model\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of drawing a few sample lines from the model, we can also compute mu for each sample and use ArviZ to plot credibility intervals (hpd intervals):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "area_s = np.linspace(start=-2, stop=3.5, num=50)\n",
    "\n",
    "mu_pred = ...\n",
    "\n",
    "# destandardize area again\n",
    "area = destandardize_area(area_s)\n",
    "\n",
    "plt.plot(area, mu_pred.mean(1)*100, alpha=0.3, c=\"k\")\n",
    "plt.scatter(berlin[\"livingSpace\"], berlin[\"totalRent_s\"]*100, s=4, alpha=0.4)\n",
    "\n",
    "az.plot_hpd(area, mu_pred.T*100, credible_interval=0.83)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "plt.xlabel('Living Area [sqm]')\n",
    "plt.ylabel('Rent [€]')\n",
    "plt.title('Visualizing the uncertainty in our model')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our model is actually very certain about where the mu is. This makes sense, because we have quite a bit of data for this relatively simple model.\n",
    "\n",
    "But this is not all the uncertainty our model is aware of. We also have the estimate sigma that tells us how far the rent could be away from this mean mu. We can compute the predicted ranges for rent as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rent_pred = stats.norm.rvs(mu_pred, trace['sigma'])\n",
    "\n",
    "plt.plot(area, mu_pred.mean(1)*100, alpha=0.3, c=\"k\")\n",
    "plt.scatter(berlin[\"livingSpace\"], berlin[\"totalRent_s\"]*100, s=4, alpha=0.7)\n",
    "\n",
    "az.plot_hpd(area, rent_pred.T*100, credible_interval=0.83, \n",
    "            fill_kwargs={\"alpha\": 0.5})\n",
    "\n",
    "\n",
    "\n",
    "plt.xlabel('Living Area [sqm]')\n",
    "plt.ylabel('Rent [€]')\n",
    "plt.title('Visualizing the uncertainty in our model')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The red area should now contain around 83% of our data if the model is good.\n",
    "\n",
    "The same way as we computed the `mu_pred` and `rent_pred`, we could also predict the rent price for a new flat. \n",
    "\n",
    "Can you predict the rent price for example for your own flat?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "area_of_my_own_flat = []\n",
    "\n",
    "# don't forget to standardize\n",
    "area_of_my_own_flat_s = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is also a less manual way to get predictions on new data. Unfortunately, it is not as straight-forward as with scikit-learn. You will need to set up the predictor variables used inside the model as a data container:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as mod:\n",
    "    # create data containers for predictor variable\n",
    "    area = pm.Data(\"area\", berlin[\"livingSpace_s\"])\n",
    "    \n",
    "    # rest stays the same\n",
    "    alpha = pm.Normal(\"alpha\", mu=0, sigma=10)\n",
    "    beta = pm.Normal(\"beta\", mu=0, sigma=5)\n",
    "    \n",
    "    sigma = pm.HalfNormal(\"sigma\", sigma=5)\n",
    "    \n",
    "    mu = alpha + beta*area\n",
    "    \n",
    "    rent = pm.Normal(\"rent\", mu=mu, sigma=sigma,\n",
    "                    observed = berlin[\"totalRent_s\"])\n",
    "    \n",
    "    trace = pm.sample(draws=1000, tune=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then swap the content of this data container with new data and then compute the prediction for the new data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# unfortunately there is an issue that predicting a single obs\n",
    "# doesn't work correctly https://github.com/pymc-devs/pymc3/issues/3640#issuecomment-563897443\n",
    "# so just use two\n",
    "new_area = standardize_area([74, 97])\n",
    "with mod:\n",
    "    pm.set_data({\"area\": new_area})\n",
    "    \n",
    "    post_pred = pm.sample_posterior_predictive(trace, samples=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first, it might look unhandy and more difficult to work with samples as predictions instead of point estimates. While indeed this means we often need to handle multidimensional arrays, the advantage of this is that we can see the whole probability distribution over rents for one flat.\n",
    "\n",
    "Plot the histogram over rents for the first flat!\n",
    "(Don't forget to destandardize the rent price again!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sample also allows it to answer questions such as, how likely is it that the rent of this flat would be below 600€. \n",
    "\n",
    "For this, we check how many points in our samples are smaller than 600€ and then take the mean of the vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool, eh? (Though the result is maybe a bit depressing...)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PyLadies-Bayesian-Tutorial",
   "language": "python",
   "name": "pyladies-bayesian-tutorial"
  },
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}