lectures/Lecture_6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Classical/Frequentist Statistical Inference: II\n",
    "\n",
    "*S. R. Taylor (2021)*\n",
    "\n",
    "Material in this lecture and notebook is based upon the \"Maximum Likelihood and Applications in Astronomy\" lectures of A. Connolly's & Ž. Ivezić's \"Astrostatistics & Machine Learning\" class at the University of Washington (ASTR 598, https://github.com/dirac-institute/uw-astr598-w18). Also the \"Inference\" lectures of G. Richards' \"Astrostatistics\" class at Drexel University (PHYS 440/540, https://github.com/gtrichards/PHYS_440_540), J. Bovy's mini-course on \"Statistics & Inference in Astrophysics\" at the University of Toronto (http://astro.utoronto.ca/~bovy/teaching.html). \n",
    "\n",
    "##### Reading:\n",
    "\n",
    "- [Textbook](http://press.princeton.edu/titles/10159.html) Chapter 4.\n",
    "- [David Hogg's \"Fitting A Model To Data\"](https://arxiv.org/abs/1008.4686)\n",
    "\n",
    "***Exercises required for class participation are in <font color='red'>red</font>.***\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "* [Fitting A Line To Data](#one)\n",
    "* [Goodness Of Fit & Model Comparison](#two)\n",
    "* [Confidence Estimating: Bootstrapping & Jackknifing](#three)\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting A Line To Data <a class=\"anchor\" id=\"one\"></a>\n",
    "\n",
    "\n",
    "Continuing the theme of last lecture, let's look at one of the most common tasks in statistical inference: **fitting a line to data**. We won't always fit a straight line, but they are prevalent in astronomy since we're potentially examining data over several orders of magnitude. Hence power-law relationships ($y\\propto x^\\alpha$) become linear relationships in log-log space ($\\ln y \\propto \\alpha\\ln x + \\mathrm{constant}$).\n",
    "\n",
    "Assume the noise/scatter in our measurements (the residuals) is generated by a Gaussian process, i.e.:\n",
    "\n",
    "$$ y_i = a x_i + b + r_i $$\n",
    "\n",
    "where $r_i$ is drawn from $N(0, \\sigma)$. Here, $\\sigma$ is the measurement uncertainty, which we take to be the same for all points. The data model includes a linear relationship with two parameters $a,b$: hence the model is written as $M(a,b)$.\n",
    "\n",
    "Let us compute the likelihood. First, we ask ourselves what is the probability $p(y_i|x_i, M(a, b), \\sigma)$ that a particular point $y_i$ would be measured. It is just the normal distribution:\n",
    "\n",
    "$$ p(y_i|x_i, M(a, b), \\sigma) = N(y_i - M(x)|\\sigma) = \\frac{1}{\\sqrt{2 \\pi \\sigma^2}} \\exp \\left( - \\frac{(y_i - M(x_i))^2}{2 \\sigma^2} \\right). $$\n",
    "\n",
    "The likelihood for all data points is given by the product over $N$ of these terms. Given our previous definitions we can then write down the $\\ln L$ as\n",
    "\n",
    "$$ \\ln L(a, b) = \\mathrm{constant} - \\frac{1}{2} \\sum_{i=1}^N \\frac{(y_i - M(x_i))^2}{\\sigma^2} = \\mathrm{constant} - \\frac{1}{2} \\chi^2$$\n",
    "\n",
    "This is the expression that we now ***maximize*** with respect to $a$ and $b$ to find ML estimators for those parameters. This is equivalent to ***minimizing*** the sum of the squares (the $\\chi^2$) in a *least-squares method*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'\n",
    "\n",
    "import scipy.stats\n",
    "from scipy.stats import norm\n",
    "from scipy.stats import uniform\n",
    "from scipy import optimize\n",
    "from astroML import stats as astroMLstats\n",
    "from astroML.datasets import fetch_hogg2010test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell will read in some data and peform a least-squares (or $\\chi^2$) fit of a linear relationship. It is minimizing the $\\chi^2$ (since this is Gaussian data this means it's also maximizing the likelihood) for both $a$ (the slope) and $b$ (the $y$-intercept). \n",
    "\n",
    "The data contains some **poorly modeled outliers** too, which have very different uncertainties from what we assume. We'll look at the case without outliers first to build intuition. \n",
    "\n",
    "***The outlier points are the first 4 in the dataset.***\n",
    "\n",
    "<font color='red'>Read through and understand the whole thing before you execute it.</font> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#------------------------------------------------------------\n",
    "# Get data from AstroML: this includes outliers\n",
    "data = fetch_hogg2010test()\n",
    "x = data['x'] # x data\n",
    "y = data['y'] # y data\n",
    "dy = data['sigma_y'] # uncertainties on y data\n",
    "\n",
    "# Define the standard squared-loss function.\n",
    "# This is just another name for chi^2\n",
    "def squared_loss(m, b, x, y, dy):\n",
    "    y_fit = m * x + b\n",
    "    return np.sum(((y - y_fit) / dy) ** 2, -1)\n",
    "\n",
    "# define a lambda function that defines the sum of squared errors.\n",
    "# these lambda functions are useful!\n",
    "# let's initially exclude the outliers by chopping off the first 4 points.\n",
    "f_squared = lambda beta: squared_loss(beta[0], beta[1], \n",
    "                                      x=x[4:], y=y[4:], \n",
    "                                      dy=dy[4:])\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# compute the maximum likelihood \n",
    "beta0 = (1, 30) # initial guess for a and b\n",
    "beta_squared = optimize.fmin(f_squared, beta0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#------------------------------------------------------------\n",
    "# Plot the results\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# plot the data\n",
    "ax.errorbar(x[4:], y[4:], dy[4:], \n",
    "            fmt='.k', lw=1, ecolor='gray')\n",
    "\n",
    "# plot the best fit linear relationship\n",
    "x_fit = np.linspace(0, 350, 10)\n",
    "ax.plot(x_fit, beta_squared[0] * x_fit + beta_squared[1], \n",
    "        ls='--', color='k',\n",
    "        label=\"squared loss:\\n $y=%.2fx + %.1f$\" % tuple(beta_squared))\n",
    "\n",
    "ax.set_xlim(0, 350)\n",
    "ax.set_ylim(100, 700)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "ax.legend(loc=4, prop=dict(size=14))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll now include the outlier data points to see how they contaminate the fit, i.e. we don't need to exclude any of the points from our fit. <font color='red'>Complete and execute the following cell.</font> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# complete to include the outliers\n",
    "f_squared_outlier = lambda beta_outlier: squared_loss(beta_outlier[0], \n",
    "                                                      beta_outlier[1], \n",
    "                                                      x=___, y=___, dy=___)\n",
    "\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# compute the maximum likelihood \n",
    "beta0 = (___, ___) # complete for initial guess\n",
    "beta_squared_outlier = optimize.fmin(f_squared_outlier, beta0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>Complete in order to plot the regular and outlier fits on the same figure.</font> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#------------------------------------------------------------\n",
    "# Plot the results\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# plot the data without outliers in gray\n",
    "ax.errorbar(x[___], y[___], dy[___], \n",
    "            fmt='.k', lw=1, ecolor='gray')\n",
    "# plot the outliers in red\n",
    "ax.errorbar(x[___], y[___], dy[___], \n",
    "            fmt='.k', lw=1, ecolor='red')\n",
    "\n",
    "x_fit = np.linspace(0, 350, 10)\n",
    "# plot the regular fit from before without outliers\n",
    "ax.plot(x_fit, beta_squared[0] * x_fit + beta_squared[1], \n",
    "        ls='--', color='k',\n",
    "        label=\"squared loss:\\n $y=%.2fx + %.1f$\" % tuple(beta_squared))\n",
    "# plot the fit that includes outliers\n",
    "ax.plot(x_fit, ___, \n",
    "        ls='--', color='red',\n",
    "        label=\"squared loss with outliers:\\n $y=%.2fx + %.1f$\" % tuple(___))\n",
    "\n",
    "ax.set_xlim(0, 350)\n",
    "ax.set_ylim(100, 700)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "ax.legend(loc=4, prop=dict(size=14))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How do we deal with outliers: modifying the likelihood\n",
    "\n",
    "- Summing the squares of the residuals ($\\chi^2=\\sum_{i=1}^N (y_i - M(x_i))^2/\\sigma^2$) is sensitive to outliers.\n",
    "\n",
    "- The $\\chi^2$ is an example of an $L_2$ norm (different L, not the likelihood), where the [$L_p$ norm](https://www.wikiwand.com/en/Lp_space) can be defined as $\\sum_{i=1}^N (y_i -M(x_i))^p/\\sigma^p$.\n",
    "\n",
    "- A number of approaches exist for correcting for outliers. See in particular [this](https://arxiv.org/abs/1008.4686) great paper. \n",
    "\n",
    "- These include \"sigma-clipping\", using interquartile ranges (which you've seen in previous lectures), taking the median of solutions of subsets of the data, and least trimmed squares (which searchs for the subset of points that minimizes $\\sum_i^N (y_i - \\theta_ix_i)^2$).\n",
    "\n",
    "\n",
    "We can also change the **likelihood** to reduce the weight of outliers. This is known as the **Huber loss function**\n",
    "\n",
    "$$ \\sum_{i=1}^N e(y_i|y),$$\n",
    "\n",
    "where\n",
    "\n",
    "$$ \n",
    "e(t) = \\left\\{\n",
    "            \\begin{array}{ll}\n",
    "            \\frac{1}{2} t^2 & \\mbox{if} \\; |t| \\leq c, \\\\\n",
    "            c|t| - \\frac{1}{2} c^2 & \\mbox{if} \\; |t| \\geq c,\n",
    "            \\end{array}\n",
    "        \\right ) \n",
    "$$\n",
    "\n",
    "and \n",
    "\n",
    "$$ t = \\left|\\frac{y-M(x)}{\\sigma}\\right|.$$\n",
    "\n",
    "This is continuous and differentiable, and transitions to an $L_1$ norm ($\\sum_{i=1}^N |y_i - M(x_i)|$) for large data excursions, which downweights the outlier points.\n",
    "\n",
    "In the image below, the blue curve is the usual $\\chi^2$, $L_2$ parabola shape that we attempt to minimize, whereas the green is the alternative Huber loss function that downweights outliers.\n",
    "\n",
    "![title](https://upload.wikimedia.org/wikipedia/commons/c/cc/Huber_loss.svg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the log-likelihood via the Huber loss function\n",
    "def huber_loss(m, b, x, y, dy, c=2):\n",
    "    y_fit = m * x + b # model\n",
    "    t = abs((y - y_fit) / dy) # define t\n",
    "    mask = t > c # define mask for large excursion points\n",
    "    \n",
    "    # apply masking to different points\n",
    "    # note '~' means 'not'\n",
    "    return np.sum((~mask) * (0.5 * t ** 2) - \\\n",
    "                  (mask) * c * (0.5 * c - t), -1)\n",
    "\n",
    "# lambda function for huber loss \n",
    "f_huber = lambda beta: huber_loss(beta[0], beta[1], \n",
    "                                  x=x, y=y, dy=dy, c=1)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# compute the maximum likelihood using the huber loss\n",
    "beta0 = (1, 30)\n",
    "beta_huber = optimize.fmin(f_huber, beta0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#------------------------------------------------------------\n",
    "# Plot the results\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# plot the data without outliers in gray\n",
    "ax.errorbar(x[4:], y[4:], dy[4:], \n",
    "            fmt='.k', lw=1, ecolor='gray')\n",
    "# plot the outliers in red\n",
    "ax.errorbar(x[:4], y[:4], dy[:4], \n",
    "            fmt='.k', lw=1, ecolor='red')\n",
    "\n",
    "x_fit = np.linspace(0, 350, 10)\n",
    "# plot the regular fit from before without outliers\n",
    "ax.plot(x_fit, beta_squared[0] * x_fit + beta_squared[1], \n",
    "        ls='--', color='k',\n",
    "        label=\"squared loss:\\n $y=%.2fx + %.1f$\" % tuple(beta_squared))\n",
    "# plot the fit that includes outliers\n",
    "ax.plot(x_fit, beta_squared_outlier[0] * x_fit + beta_squared_outlier[1], \n",
    "        ls='--', color='red',\n",
    "        label=\"squared loss with outliers:\\n $y=%.2fx + %.1f$\" % tuple(beta_squared_outlier))\n",
    "# plot the fit with the huber likelihood, downweighting outliers\n",
    "ax.plot(x_fit, beta_huber[0] * x_fit + beta_huber[1], \n",
    "        ls='-', color='blue',\n",
    "        label=\"huber loss:\\n $y=%.2fx + %.1f$\" % tuple(beta_huber))\n",
    "\n",
    "ax.set_xlim(0, 350)\n",
    "ax.set_ylim(100, 700)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "ax.legend(loc=4, prop=dict(size=14))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>Did this help to tame the outliers? Discuss with your breakout room colleagues about why.</font> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goodness Of Fit & Model Comparison <a class=\"anchor\" id=\"two\"></a>\n",
    "\n",
    "The MLE approach tells us what the \"best\" model parameters are, but not how good the fit actually is. If the model is wrong, \"best\" might not be particularly revealing! Remember this acronym from the dawn of the computer age: ***GIGO = Garbage In, Garbage Out***. \n",
    "\n",
    "If we have a crappy model then we shouldn't expect a good fit to the data. For example, if you have $N$ points drawn from a linear distribution, you can always fit the data perfectly with an $N-1$ order polynomial. But that won't help you predict future measurements.\n",
    "\n",
    "We can describe the **goodness of fit** in words as simply the followng\n",
    "\n",
    "> *The goodness of fit tells us whether or not it is likely to have obtained the maximum (log-)likelihood $\\ln L^0$ by randomly drawing from the data.* \n",
    "\n",
    "Using the best-fit parameters of a model, the maximum likelihood value $L^0$ should not be an unlikely occurence. If it is, then our model is not describing the data well. Thus we need to know the *distribution* of $\\ln L$ and not just its maximum. For the Gaussian case we have just described, we do a standard transform of variables and compute the so-called $z$ score for each data point (basically the number of standard deviations away from the mean that this point is), writing \n",
    "\n",
    "$$z_i = (x_i-\\mu)/\\sigma,$$ \n",
    "\n",
    "then\n",
    "\n",
    "$$\\ln L = {\\rm constant} - \\frac{1}{2}\\sum_{i=1}^N z_i^2 = {\\rm constant} - \\frac{1}{2}\\chi^2.$$\n",
    "\n",
    "Here, $\\chi^2$ is the thing whose distribution we discussed in early lectures.\n",
    "\n",
    "**So for Gaussian uncertainties, $\\ln L$ is distributed as $\\chi^2$.**  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run ./code/fig_chi2_distribution.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***The mean of the $\\chi^2$ distribution is $N-k$ and its standard deviation is $\\sqrt{2(N-k)}$.***\n",
    "\n",
    "We define the $\\chi^2$ per degree of freedom, $\\chi^2_\\mathrm{dof}$, as\n",
    "\n",
    "$$\\chi^2_\\mathrm{dof} = \\frac{1}{N-k}\\sum_{i=1}^N z^2_i.$$\n",
    "\n",
    "where $k$ is the number of model parameters determined from the data.\n",
    "\n",
    "- For a good fit, we would expect that $\\chi^2_\\mathrm{dof}\\approx 1$. \n",
    "- If $\\chi^2_\\mathrm{dof}$ is significantly larger than 1, or $(\\chi^2_\\mathrm{dof}-1)>> \\sqrt{2/(N-k)}$, then it is likely that we are not using the correct model.\n",
    "- If data uncertainties are **(over)under-estimated** then this can lead to improbably **(low)high $\\chi^2_\\mathrm{dof}$**, as seen below.  \n",
    "\n",
    "![Ivezic, Figure 4.1](http://www.astroml.org/_images/fig_chi2_eval_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read through and execute the following cell to compute $\\chi^2_\\mathrm{dof}-1$ for the line fitted without outliers and fitted with outliers. These will be compared to the standard deviation of the $\\chi^2_\\mathrm{dof}-1$ distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------\n",
    "# Get data from AstroML: this includes outliers\n",
    "data = fetch_hogg2010test()\n",
    "x = data['x'] # x data\n",
    "y = data['y'] # y data\n",
    "dy = data['sigma_y'] # uncertainties on y data\n",
    "\n",
    "# number of data points *without* outliers\n",
    "N = x[4:].shape[0]\n",
    "# number of data points *with* outliers\n",
    "N_outlier = x.shape[0]\n",
    "# number of model parameters (a,b)\n",
    "k = 2 \n",
    "\n",
    "# chi2 per dof *without* outliers\n",
    "chi2 = squared_loss(beta_squared[0], \n",
    "                    beta_squared[1], \n",
    "                    x=x[4:], y=y[4:], dy=dy[4:])\n",
    "chi2dof = chi2 / (N-k)\n",
    "\n",
    "# chi2 per dof *with* outliers\n",
    "chi2_outlier = squared_loss(beta_squared_outlier[0], \n",
    "                            beta_squared_outlier[1], \n",
    "                            x=x, y=y, dy=dy)\n",
    "chi2dof_outlier = chi2_outlier / (N_outlier - k)\n",
    "\n",
    "# without outliers\n",
    "print(chi2dof-1, np.sqrt(2/(N-k)))\n",
    "\n",
    "# with outliers\n",
    "print(chi2dof_outlier-1, np.sqrt(2/(N_outlier - k)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>For each case, divide the first value by the second value to get the sigma deviation from the desired value of zero. The lower the better. </font> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute your sigma values here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Comparison\n",
    "\n",
    "The distribution of $\\ln L$ can only be related to the $\\chi^2$ distribution whenever the likelihood is Gaussian. For non-Gaussian likelihoods we can still rank different models in terms of their respectively computed maximum likelihood values, $L^0$. This is only really fair if the models have the same number of parameters.\n",
    "\n",
    "Let's do that for the data above that contains outliers, both for a model based on the naive squared loss function ($\\chi^2$), and a Huber loss function model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute a quantity related to lnL0 for the \n",
    "# squared loss function\n",
    "lnL_sq = - squared_loss(beta_squared_outlier[0], \n",
    "                        beta_squared_outlier[1], \n",
    "                        x=x, y=y, dy=dy)\n",
    "\n",
    "# compute a quantity related to lnL0 for the \n",
    "# Huber loss function\n",
    "lnL_huber = - huber_loss(beta_huber[0], \n",
    "                         beta_huber[1], \n",
    "                         x=x, y=y, dy=dy)\n",
    "\n",
    "print(\"Maximum log likelihood for squared loss = %.1f\" % lnL_sq)\n",
    "print(\"Maximum log likelihood for Huber loss = %.1f\" % lnL_huber)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need a modified approach whenever the models we're comparing have different numbers of parameters. Such an approach should account for **model complexity** and **effectively penalize models with additional parameters that are not constrained by the data**. This is often called an ***Occam penalty***, because we're trying to incorporate [Occam's Razor](https://www.wikiwand.com/en/Occam%27s_razor#:~:text=Occam's%20razor%2C%20Ockham's%20razor%2C%20Ocham's,is%20usually%20the%20right%20one.).\n",
    "\n",
    "> *All else being equal (i.e., each model fits the data equally well), the less complex model is favored.*\n",
    "\n",
    "\n",
    "We'll meet this in extensive detail later, especially in a Bayesian context. But a popular general-purpose tool for model comparison is the **Akaike Information Criterion** (AIC):\n",
    "\n",
    "$$ \\mathrm{AIC}_M \\equiv -2\\ln[L^0(M)] + 2k + \\frac{2k(k+1)}{N-k-1}, $$\n",
    "\n",
    "where $k$ is the number of model parameters and $N$ is the number of data points.\n",
    "\n",
    "- For a Gaussian distribution, the first term is equal to $\\chi^2$.\n",
    "- **The model with lowest AIC is the most favored.**\n",
    "- If all models are equally successful at fitting the data (equal $L^0$ values) then the second and third terms penalize model complexity such that the model with fewest free parameters wins. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence Estimating: Bootstrapping & Jackknifing <a class=\"anchor\" id=\"three\"></a>\n",
    "\n",
    "Most ways of computing confidence limits and ucnertainties on measured model parameters assume that the distribution is Gaussian and our samples are large. But even if that is not the case, we can still compute good confidence intervals (e.g., $a<x<b$ with 95\\% confidence) using ***resampling*** strategies.\n",
    "\n",
    "Remember that we have a data set $\\{x_i\\}$ from which we have estimated the distribution as $f(x)$ for a true distribution $h(x)$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bootstrap Method\n",
    "\n",
    "In bootstrapping, we map the uncertainty of model parameters by re-sampling from our data (with replacement) $B$ times. This is relatively new (Efron, 1979) and was named in inspiration of \"*pulling oneself up by one's bootstraps*\".\n",
    "\n",
    "When sampling from observed data of length $N$ with replacement, there are $N!$ distinct combinations of new observed datasets, and there is only a $N!/N^N$ probability of getting back the original dataset (even for $N=10$ this is only $0.00036$). \n",
    "\n",
    "With $B$ new bootstrap datasets, we compute our statistics on each to obtain $B$ measures of our parameters. So, if we have $i=1,\\dots,N$ data points in $\\{x_i\\}$, we draw $N$ of them to make a new sample, where some values of $\\{x_i\\}$ will be used more than once (and this is ok).\n",
    "\n",
    "**EXAMPLE**\n",
    "- We have an original dataset of $1000$ points drawn from a Gaussian distribution. \n",
    "- We can only measure the standard deviation of the distribution once with this dataset. \n",
    "- Our previous strategies showed that we can use a *Fisher estimate* or an *analytic estimate* of the sample standard deviation uncertainty. \n",
    "- Instead, we use the bootstrap method to resample the data $10,000$ times, and compute the standard deviation and $\\sigma_G$ on each new dataset. \n",
    "- This will map out the uncertainty distribution of those statistics, allowing us to quote confidence limits in our actual measured value. \n",
    "\n",
    "<font color='red'>Read through the next cell carefully to understand what's happening when you execute it.</font> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ivezic v2, Figure 4.3, modified slightly by GTR and SRT\n",
    "# %load ../code/fig_bootstrap_gaussian.py\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "\n",
    "from astroML.resample import bootstrap\n",
    "from astroML.stats import sigmaG\n",
    "\n",
    "m = 1000  # number of points\n",
    "n = 10000  # number of bootstraps\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# sample values from a normal distribution\n",
    "np.random.seed(123)\n",
    "data = norm(0, 1).rvs(m)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Compute bootstrap resamplings of data\n",
    "mu1_bootstrap = bootstrap(data, n, np.std, \n",
    "                          kwargs=dict(axis=1, ddof=1))\n",
    "mu2_bootstrap = bootstrap(data, n, sigmaG, \n",
    "                          kwargs=dict(axis=1))\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Compute the theoretical expectations for the two distributions\n",
    "xgrid = np.linspace(0.8, 1.2, 1000)\n",
    "\n",
    "sigma1 = 1. / np.sqrt(2 * (m - 1))\n",
    "pdf1 = norm(1, sigma1).pdf(xgrid)\n",
    "\n",
    "sigma2 = 1.06 / np.sqrt(m)\n",
    "pdf2 = norm(1, sigma2).pdf(xgrid)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot the results\n",
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "\n",
    "ax.hist(mu1_bootstrap, bins=50, density=True, \n",
    "        histtype='step', color='blue', ls='dashed', \n",
    "        label=r'$\\sigma\\ {\\rm (std. dev.)}$')\n",
    "ax.plot(xgrid, pdf1, color='gray')\n",
    "\n",
    "ax.hist(mu2_bootstrap, bins=50, density=True, \n",
    "        histtype='step', color='red', \n",
    "        label=r'$\\sigma_G\\ {\\rm (quartile)}$')\n",
    "ax.plot(xgrid, pdf2, color='gray')\n",
    "\n",
    "ax.set_xlim(0.82, 1.18)\n",
    "ax.set_xlabel(r'$\\sigma$')\n",
    "ax.set_ylabel(r'$p(\\sigma|x,I)$')\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above shows the bootstrap uncertainty estimates for the sample standard deviation $\\sigma$ (blue dashed) and $\\sigma_G$ (red solid). The thin grey lines show Gaussians with the theoretical widths determined as $\\sigma_s = s\\, / \\sqrt{2(N - 1)}$ and $\\sigma_{\\sigma_G} = 1.06\\, s\\, / \\sqrt{N}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jackknife Method \n",
    "\n",
    "This is similar to bootsrapping except that we don't use a sample size of $N$. Rather, ***we leave off one or more of the observations*** from $\\{x_i\\}$. As with bootstrap, we do this multiple times, generating samples from which we can determine our uncertainties.\n",
    "\n",
    "- If we leave out just one datapoint each time, we can make $N$ such datasets that each contain $(N-1)$ data points.\n",
    "- We compute our statistic or find our model parameters in each of these jackknife datasets.\n",
    "\n",
    "For jackknifing by leaving out one datapoint, the bias-corrected jackknife estimate of a statistic $\\alpha$ is\n",
    "\n",
    "$$ \\alpha^J = \\alpha_N + \\Delta\\alpha,$$\n",
    "\n",
    "where $\\alpha_N$ is the statistic computed on the original dataset, and \n",
    "\n",
    "$$ \\Delta\\alpha = (N-1)\\left( \\alpha_N - \\frac{1}{N}\\sum_{i=1}^N \\alpha^*_i \\right)$$\n",
    "\n",
    "where $\\{\\alpha_i^*\\}$ are the statistics computed on the $N$ jackknife datasets. \n",
    "\n",
    "For asymptoticically normal estimators, the standard error on $\\alpha^J$ is\n",
    "\n",
    "$$ \\sigma_\\alpha = \\sqrt{\\frac{1}{N(N-1)}\\sum_{i=1}^N [N\\alpha_N - \\alpha^J - (N-1)\\alpha^*_i]^2}.$$\n",
    "\n",
    "Confidence limits on $\\alpha$ can be computed using the Student's $t$ distribution with $t=(\\alpha-\\alpha^J)/\\sigma_\\alpha$ and $(N-1)$ degrees of freedom.\n",
    "\n",
    "**NOTES**\n",
    "- The jackknife standard error is more reliable than the bias correction.\n",
    "- Standard jackknife (removing one point) does well for statistics like the mean and standard deviation, but is poor with rank-based statistics (e.g., median, quantiles, $\\sigma_G$). \n",
    "- This is because (as we've seen) rank-based statistics are not very sensitive to adding/removing single data points. So jackknifing can give many identical values of the statistic!\n",
    "- This can be resolved by modifying the jackknife to leave off more than one datapoint.\n",
    "\n",
    "\n",
    "**EXAMPLE**\n",
    "- We compute jackknife uncertainty estimates for the width of a Gaussian distribution, using the same data as in the bootstrap example.\n",
    "- The $N$ jackknife estimates of $\\sigma$ and $\\sigma_G$ will be shown, but the `astroML` jackknife method automatically incorporates the bias correction mentioned above.\n",
    "\n",
    "<font color='red'>Read through the next cell carefully to understand what's happening when you execute it.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ivezic v2, Figure 4.4, modified slightly by SRT\n",
    "# %load ./code/fig_jackknife_gaussian.py\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "\n",
    "from astroML.resample import jackknife\n",
    "from astroML.stats import sigmaG\n",
    "\n",
    "m = 1000  # number of points\n",
    "n = 10000  # number of bootstraps\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# sample values from a normal distribution\n",
    "np.random.seed(123)\n",
    "data = norm(0, 1).rvs(m)\n",
    "\n",
    "# mu1 is the mean of the standard-deviation-based width\n",
    "mu1, sigma_mu1, mu1_raw = jackknife(data, np.std,\n",
    "                                    kwargs=dict(axis=1, ddof=1),\n",
    "                                    return_raw_distribution=True)\n",
    "\n",
    "pdf1_theory = norm(1, 1. / np.sqrt(2 * (m - 1)))\n",
    "pdf1_jackknife = norm(mu1, sigma_mu1)\n",
    "print(mu1, sigma_mu1)\n",
    "\n",
    "# mu2 is the mean of the interquartile-based width\n",
    "#  WARNING: do not use the following in practice.  This example\n",
    "#           shows that jackknife fails for rank-based statistics.\n",
    "mu2, sigma_mu2, mu2_raw = jackknife(data, sigmaG,\n",
    "                                    kwargs=dict(axis=1),\n",
    "                                    return_raw_distribution=True)\n",
    "pdf2_theory = norm(data.std(), 1.06 / np.sqrt(m))\n",
    "pdf2_jackknife = norm(mu2, sigma_mu2)\n",
    "print(mu2, sigma_mu2)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# plot the results\n",
    "print(\"mu_1 mean: %.2f +- %.2f\" % (mu1, sigma_mu1))\n",
    "print(\"mu_2 mean: %.2f +- %.2f\" % (mu2, sigma_mu2))\n",
    "\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "fig.subplots_adjust(left=0.11, right=0.95, bottom=0.2, top=0.9,\n",
    "                    wspace=0.25)\n",
    "\n",
    "# plot all jacnknife results on histogram\n",
    "ax = fig.add_subplot(121)\n",
    "ax.hist(mu1_raw, np.linspace(0.996, 1.008, 100),\n",
    "        label=r'$\\sigma^*\\ {\\rm (std.\\ dev.)}$',\n",
    "        histtype='stepfilled', fc='white', \n",
    "        ec='black', density=False)\n",
    "ax.hist(mu2_raw, np.linspace(0.996, 1.008, 100),\n",
    "        label=r'$\\sigma_G^*\\ {\\rm (quartile)}$',\n",
    "        histtype='stepfilled', fc='gray', \n",
    "        density=False)\n",
    "ax.legend(loc='upper left', handlelength=2)\n",
    "\n",
    "ax.xaxis.set_major_locator(plt.MultipleLocator(0.004))\n",
    "ax.set_xlabel(r'$\\sigma^*$')\n",
    "ax.set_ylabel(r'$N(\\sigma^*)$')\n",
    "ax.set_xlim(0.998, 1.008)\n",
    "ax.set_ylim(0, 550)\n",
    "\n",
    "# plot the bias-corrected jackknife estimates\n",
    "ax = fig.add_subplot(122)\n",
    "x = np.linspace(0.45, 1.15, 1000)\n",
    "ax.plot(x, pdf1_jackknife.pdf(x),\n",
    "        color='blue', ls='dashed', \n",
    "        label=r'$\\sigma\\ {\\rm (std.\\ dev.)}$',\n",
    "        zorder=2)\n",
    "ax.plot(x, pdf1_theory.pdf(x), \n",
    "        color='gray', zorder=1)\n",
    "ax.plot(x, pdf2_jackknife.pdf(x),\n",
    "        color='red', \n",
    "        label=r'$\\sigma_G\\ {\\rm (quartile)}$', \n",
    "        zorder=2)\n",
    "ax.plot(x, pdf2_theory.pdf(x), \n",
    "        color='gray', zorder=1)\n",
    "plt.legend(loc='upper left', handlelength=2)\n",
    "\n",
    "\n",
    "ax.set_xlabel(r'$\\sigma$')\n",
    "ax.set_ylabel(r'$p(\\sigma|x,I)$')\n",
    "ax.set_xlim(0.45, 1.15)\n",
    "ax.set_ylim(0, 24)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **left panel** shows a histogram over all $N$ jackknife datasets of the widths determined using the sample standard deviation, and using $\\sigma_G$ (from the interquartile range). \n",
    "\n",
    "The **right panel** shows the bias-corrected jackknife estimates (see above) for the two methods. The gray lines show the theoretical results as before in the bootstrap example. The jackknife result for $\\sigma$ matches the theoretical result almost exactly, but note the failure of the jackknife to correctly estimate $\\sigma_G$. *Jackknifing does poorly for rank-based statistics!*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### When to bootstrap or jackknife?\n",
    "\n",
    "- Jackknife estimates are usually easier to calculate, easier to apply for complex sampling schemes, and automtically remove bias.\n",
    "- Bootstrap is better for computing confidence intervals because it maps out the full distribution of the statistic instead of assuming asymptotic normality.\n",
    "- Bootstrap is random resampling so gives different results each time. Whereas jackknifing gives repeatable results.\n",
    "\n",
    "It is generally a good idea to use both methods and compare the results. Use either/both with caution with $N$ is small!"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:astr8070] *",
   "language": "python",
   "name": "conda-env-astr8070-py"
  },
  "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.8.5"
  },
  "livereveal": {
   "scroll": true,
   "start_slideshow_at": "selected",
   "theme": "sky"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}