lectures/Lecture_8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Bayesian Statistical Inference: I\n",
    "\n",
    "*S. R. Taylor (2021)*\n",
    "\n",
    "Material in this lecture and notebook is based upon the \"Introduction To Bayesian Inference\" 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\" and \"Inference2\" 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 5.\n",
    "- [David Hogg's \"Fitting A Model To Data\"](https://arxiv.org/abs/1008.4686)\n",
    "- [Jake VanderPlas's workshop \"Bayesian Astronomy\"](https://github.com/jakevdp/BayesianAstronomy)\n",
    "- [Jake VanderPlas's blog \"Frequentism and Bayesianism: A Practical Introduction\"](http://jakevdp.github.io/blog/2014/03/11/frequentism-and-bayesianism-a-practical-intro/)\n",
    "\n",
    "\n",
    "##### Highly recommended supplemental background reading:\n",
    "\n",
    "- [Jake VanderPlas: \"Frequentism and Bayesianism: A Python-driven Primer\"](https://arxiv.org/abs/1411.5018)\n",
    "- [Hogg, Bovy and Lang: \"Data analysis recipes: Fitting a model to data\"](https://arxiv.org/abs/1008.4686)\n",
    "\n",
    "\n",
    "##### For those who want to dive deep:\n",
    "\n",
    "- [D. Sivia and J. Skilling: \"Data Analysis: A Bayesian Tutorial\"](https://www.amazon.com/Data-Analysis-Bayesian-Devinderjit-Sivia/dp/0198568320)\n",
    "- [E.T. Jaynes: \"Probability Theory: The Logic of Science\"](http://bayes.wustl.edu/etj/prob/book.pdf)\n",
    "- [E.T. Jaynes: \"Confidence Intervals vs. Bayesian intervals\"](http://bayes.wustl.edu/etj/articles/confidence.pdf)\n",
    "- [This great explanation of confidence levels versus credible regions on Stackexchange](https://stats.stackexchange.com/questions/2272/whats-the-difference-between-a-confidence-interval-and-a-credible-interval/2287#2287)\n",
    "\n",
    "---\n",
    "\n",
    "## Contents\n",
    "* [From Bayes Rule To Bayesian Inference](#one)\n",
    "* [Bayesian Priors: What Are They & How Do I Choose Them?](#two)\n",
    "* [Bayesian Credible Regions](#three)\n",
    "* [Simple Parameter Estimation Examples](#four)\n",
    "\n",
    "---\n",
    "\n",
    "\n",
    "***Exercises required for class participation are in <font color='red'>red</font>.***\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## From \"Bayes Rule\" To \"Bayesian Inference\" <a class=\"anchor\" id=\"one\"></a>\n",
    "\n",
    "We had that \n",
    "\n",
    "$$p(x,y) = p(x|y)p(y) = p(y|x)p(x)$$\n",
    "\n",
    "We also used the notation for intersection ($p$ that both $A$ **and** $B$ will happen) \n",
    "\n",
    "$$p(A \\cap B) \\equiv p(A,B) = p(A|B)p(B) = p(B|A)p(A)$$\n",
    "\n",
    "We can define the **marginal probability** as\n",
    "\n",
    "$$p(x) = \\int p(x,y)dy,$$\n",
    "\n",
    "where **marginal means essentially projecting on to one axis**, and **conditional means taking a slice at a fixed value of one axis**.\n",
    "\n",
    "We can re-write this as\n",
    "\n",
    "$$p(x) = \\int p(x|y)p(y) dy$$\n",
    "\n",
    "Since $$p(x|y)p(y) = p(y|x)p(x)$$ we can write that\n",
    "\n",
    "$$p(y|x) = \\frac{p(x|y)p(y)}{p(x)} = \\frac{p(x|y)p(y)}{\\int p(x|y)p(y) dy}$$\n",
    "\n",
    "which in words says that\n",
    "\n",
    "> the (conditional) probability of $y$ given $x$ is just the (conditional) probability of $x$ given $y$ times the (marginal) probability of $y$ divided by the (marginal) probability of $x$, where the latter is just the integral of the numerator.\n",
    "\n",
    "This is **Bayes' rule**, which itself is not at all controversial, though its application can be as we'll discuss in detail. \n",
    "\n",
    "<font color='red'>Btw, what are the units of the various terms in the above expression? Discuss this in your breakout rooms.</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Recall Maximum Likelihood Estimation (MLE) applied to a Heteroscedastic Gaussian\n",
    "\n",
    "Assume $N$ measurements, $\\{x_i\\}$, where the uncertainty for each measurement is Gaussian with\n",
    "a known $\\sigma_i$. The likelihood of one measurement is \n",
    "\n",
    "$$L \\equiv p(x_i|\\mu,\\sigma_i) = \\frac{1}{\\sigma_i\\sqrt{2\\pi}} \\exp\\left(\\frac{-(x_i-\\mu)^2}{2\\sigma_i^2}\\right).$$\n",
    "\n",
    "\n",
    "and therefore the likelihood of all N measurements is *(little clue here for HW3 if you haven't done it yet...)*\n",
    "\n",
    "$$L \\equiv p(\\{x_i\\}|\\mu,\\sigma_i) = \\prod_{i=1}^N \\frac{1}{\\sigma_i\\sqrt{2\\pi}} \\exp\\left(\\frac{-(x_i-\\mu)^2}{2\\sigma_i^2}\\right).$$"
   ]
  },
  {
   "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 numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.stats import norm\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's draw a homoscedastic sample of $\\{x_i\\}$ from a Gaussian and see what happens with $L$. First generate a sample of $N$ points drawn from $\\mathcal{N}(\\mu,\\sigma)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleSize = 3\n",
    "mu = 1.0\n",
    "sigma = 0.2 \n",
    "sample = norm(mu, sigma).rvs(sampleSize) \n",
    "print(sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we did in the MLE lectures, let's now compute probabilities for each point centered at the measured value across a grid, and multiply the probabilities together to find the likelihood for the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "muGrid = np.linspace(0,2,1000)\n",
    "\n",
    "L1 = norm(sample[0], sigma).pdf(muGrid) \n",
    "L2 = norm(sample[1], sigma).pdf(muGrid) \n",
    "L3 = norm(sample[2], sigma).pdf(muGrid) \n",
    "L = L1 * L2 * L3\n",
    "\n",
    "# plot\n",
    "fig, ax = plt.subplots(figsize=(7, 5))\n",
    "plt.plot(muGrid, L1, ls='-', c='green', \n",
    "         label=r'$L(x_1)$')\n",
    "plt.plot(muGrid, L2, ls='-', c='red', \n",
    "         label=r'$L(x_2)$')\n",
    "plt.plot(muGrid, L3, ls='-', c='blue', \n",
    "         label=r'$L(x_3)$')\n",
    "plt.plot(muGrid, L, ls='-', c='black', \n",
    "         label=r'$L(\\{x\\})$')\n",
    "\n",
    "plt.xlim(0.2, 1.8)\n",
    "plt.ylim(0, 8.0)\n",
    "plt.xlabel('$\\mu$')\n",
    "plt.ylabel(r'$p(x_i|\\mu,\\sigma)$')\n",
    "plt.title('MLE for Gaussian Distribution')\n",
    "plt.legend()\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>But what if I told you that mu>0.9? That's prior information! Complete and execute the following. </font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "muMin = 0.9\n",
    "L1[muGrid < ___] = ___\n",
    "L2[muGrid < ___] = ___\n",
    "L3[muGrid < ___] = ___\n",
    "L = ___ * ___ * ___\n",
    "\n",
    "# plot\n",
    "fig, ax = plt.subplots(figsize=(7, 5))\n",
    "plt.plot(muGrid, L1, ls='-', c='green', label=r'$L(x_1)$')\n",
    "plt.plot(muGrid, L2, ls='-', c='red', label=r'$L(x_2)$')\n",
    "plt.plot(muGrid, L3, ls='-', c='blue', label=r'$L(x_3)$')\n",
    "plt.plot(muGrid, L, ls='-', c='black', label=r'$L(\\{x\\})$')\n",
    "\n",
    "plt.xlim(0.2, 1.8)\n",
    "plt.ylim(0, 8.0)\n",
    "plt.xlabel('$\\mu$')\n",
    "plt.ylabel(r'$p(x_i|\\mu,\\sigma)$')\n",
    "plt.title('MLE for Gaussian Distribution')\n",
    "plt.legend()\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly the *maximum is unchanged in this trivial example*, but the distribution is truncated leading to very different uncertainty estimates than what one would naively get by assuming that all $\\mu$ values are equally likely."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The Essence of the Bayesian Method \n",
    "\n",
    "- The basic premise of the Bayesian method is that probability statements are not limited to data,  but can be made for model parameters and models themselves. \n",
    "- Inferences are made by producing  probability density functions (pdfs); most notably, **model parameters are treated as random variables**.\n",
    "- These **pdfs represent our \"belief spread\" in what the model parameters are**. They have nothing to do with outcomes of repeated experiments (although the shape of resulting distributions can often coincide).\n",
    "\n",
    "\n",
    "### Brief History \n",
    "\n",
    "- The **Reverend Thomas Bayes (1702–1761)** was an English amateur mathematician who wrote a manuscript \n",
    "on how to combine an initial belief with new data to arrive at an improved belief. \n",
    "- The manuscript \n",
    "was published posthumously in 1763 and gave rise to the name Bayesian statistics. \n",
    "- **Laplace** rediscovered the Bayesian approach a decade after it was originally published, and greatly clarified some principles.  \n",
    "- Howevever, Bayesian statistics did not find its ways into mainstream science until well into the 20th century, and widespread usage has been hindered until the 1990s with the advent of cheap computing that can map out the Bayesian probability distributions.\n",
    "\n",
    "![](https://upload.wikimedia.org/wikipedia/commons/d/d4/Thomas_Bayes.gif?1613667187659)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bayesian Statistical Inference\n",
    "\n",
    "Up to now we have been computing the **likelihood** $p(D\\,|\\,M)$.  In Bayesian inference, we instead evaluate the **posterior probability** taking into account **prior** information.\n",
    "\n",
    "Recall that Bayes' Rule is:\n",
    "\n",
    "$$p(M\\,|\\,D) = \\frac{p(D\\,|\\,M)\\,p(M)}{p(D)},$$\n",
    "\n",
    "where $D$ is for data and $M$ is for model. Or in words,\n",
    "\n",
    "$${\\rm Posterior \\,\\, Probability} = \\frac{{\\rm Likelihood}\\times{\\rm Prior}}{{\\rm Evidence}}.$$\n",
    "\n",
    "If we explicitly recognize prior information, $I$, and the model parameters, $\\theta$, then we can write:\n",
    "\n",
    "$$p(M,\\theta \\,|\\,D,I) = \\frac{p(D\\,|\\,M,\\theta,I)\\,p(M,\\theta\\,|\\,I)}{p(D\\,|\\,I)},$$\n",
    "\n",
    "where we will omit the explict dependence on $\\theta$ by writing $M$ instead of $M,\\theta$ where appropriate.  However, as the prior can be expanded to \n",
    "\n",
    "$$p(M,\\theta\\,|\\,I) = p(\\theta\\,|\\,M,I)\\,p(M\\,|\\,I),$$\n",
    "\n",
    "it will still appear in the term $p(\\theta\\,|\\,M,I)$.\n",
    "\n",
    "**NOTE** \n",
    "\n",
    "We don't often care about **the evidence $p(D\\,|\\,I)$** because it does not depend on model parameters. We usually set it to $1$ for parameter estimation. **BUT** it's at the heart of Bayesian model selection (which we'll look at in the future) since it gives us a way of ranking different model descriptions of the data.  \n",
    "\n",
    "**The Bayesian Statistical Inference process** is then\n",
    "1. formulate the likelihood, $p(D\\,|\\,M,\\theta,I)$\n",
    "2. chose a prior$^1$, $p(M,\\theta\\,|\\,I)$, which incorporates *other information beyond the data in $D$*\n",
    "3. determine the posterior pdf, $p(M,\\theta \\,|\\,D,I)$\n",
    "4. search for the model parameters that maximize $p(M,\\theta \\,|\\,D,I)$ \n",
    "5. quantify the uncertainty of the model parameter estimates\n",
    "6. perform model selection to find the most apt description of the data\n",
    " \n",
    "$^1$: Note that $p(M,\\theta\\,|\\,I) = p(\\theta\\,|\\,M, I)\\, p(M\\,|\\,I)$.  \n",
    "\n",
    "Before applying this expression, we need to discuss how to choose priors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayesian priors: What Are They & How Do I Choose Them? <a class=\"anchor\" id=\"two\"></a>\n",
    "\n",
    "Priors can be **informative** or **uninformative**.  As it sounds, informative priors are based on existing information (including previously obtained data, but not the data considered right now) that might be available.  Uniformative priors can be thought of as \"default\" priors, i.e., what your prior is when you never used\n",
    "any data, e.g, a \"flat\" prior like $p(\\theta|M,I) \\propto {\\rm C}$.\n",
    "\n",
    "Detailed discussion can be found in Section 5.2 in the textbook. In general, we want our inferences to be ***data dominated*** rather than prior dominated, so we try to use ***weakly-informative priors***. There are three\n",
    "main principles used to choose a prior: \n",
    "\n",
    "\n",
    "### (i) The Principle of Indifference\n",
    "\n",
    "Essentially this means adopting a uniform prior, though you have to be a bit careful.  Saying that an asteroid is equally likely to hit anywhere on the Earth is not the same as saying that all latitudes of impact are equally likely.  \n",
    "\n",
    "Assuming $1/6$ for a six-side die, or $1/2$ for heads and tails of a fair coin, would be an example of indifference.\n",
    "\n",
    "### (ii) The Principle of Invariance (or Consistency)\n",
    "\n",
    "This applies to location and scale invariance.  \n",
    "\n",
    "**Location invariance** suggests a uniform prior, within the accepted bounds: $p(\\theta|I) \\propto 1/(\\theta_{max}-\\theta_{min})$ for $\\theta_{min} \\le \\theta \\le \\theta_{max}$. \n",
    "\n",
    "**Scale invariance** gives us priors that look like $p(\\theta|I) \\propto 1/\\theta$, which implies a uniform\n",
    "prior for ln($\\theta$), i.e. a prior that gives equal weight over many orders of magnitude. \n",
    "\n",
    "### (iii) The Principle of Maximum Entropy\n",
    "\n",
    "We will not discuss it here - for more details, see Section 5.2.2 in the textbook.\n",
    " \n",
    "It is often true that Bayesian analysis and traditional MLE are essentially equivalent.  \n",
    "However, in some cases, considering the priors can have significant consequences, as\n",
    "we will see later. \n",
    "\n",
    "We will skip examples of very steep priors and their consequences called in astronomy\n",
    "literature **Eddington-Malmquist** and **Lutz-Kelker** biases (see Chapter 5 in the textbook\n",
    "if you are interested). \n",
    "\n",
    "### Conjugate Priors\n",
    "\n",
    "In special combinations of priors and likelihood functions, the resulting posterior probability distribution is from the same function family as the prior. These priors are called **conjugate priors** and give a convenient way for generalizing computations. There are exhaustive tables [here](https://www.wikiwand.com/en/Conjugate_prior#/Table_of_conjugate_distributions). \n",
    "\n",
    "**EXAMPLE**\n",
    "\n",
    "If the likelihood is Gaussian and the prior function is Gaussian, then so too is the posterior distribution! So the conjugate prior for a Gaussian likelihood is a Gaussian.\n",
    "\n",
    "For data drawn from a Gaussian likelihood equal to $\\mathcal{N}(\\bar{x},s)$ (where $\\bar{x}$ is the sample mean and $s$ is the sample standard deviation), with a prior on the underlying parameters $\\mathcal{N}(\\mu_p,\\sigma_p)$, the posterior is $\\mathcal{N}(\\mu^0,\\sigma^0)$, where\n",
    "\n",
    "$$\\mu^0 = \\frac{\\mu_p/\\sigma_p^2 + \\bar{x}/s^2}{1/\\sigma_p^2 + 1/s^2},\\quad \\sigma^0 = \\left( 1/\\sigma_p^2 + 1/s^2 \\right)^{-1/2} $$\n",
    "\n",
    "\n",
    "### Hierarchical Bayes\n",
    "\n",
    "You may hear of ***hierarchical Bayesian modeling*** a great deal these days. It's become sort of a buzz-term for people wanting to sound fancy, and [I am not immune to this](https://arxiv.org/abs/1806.08365). \n",
    "\n",
    "But there's nothing terribly fancy about it. Essentially, we will look at employing prior distributions today that have fixed shapes (e.g. Gaussian distributions centered around fixed values with fixed widths). \n",
    "\n",
    "But in hierarchical Bayesian modeling, the parameters of the prior distribution (called ***hyperparameters***) become part of the search! The data informs not only properties of individual events but also the shape of the prior. Those prior parameters then get their own priors, called ***hyperpriors***. \n",
    "\n",
    "**The whole analysis is then hierarchical, corresponding to multiple layers of inference.** \n",
    "\n",
    "For example, we have lots of [exoplanet discoveries](https://exoplanets.nasa.gov/discovery/discoveries-dashboard/). \n",
    "- Each of those discoveries started with some lightcurve data, where we fit a likelihood model to that data to deduce parameters of the system. \n",
    "- The prior on those parameters was likely weakly informative. \n",
    "- ***BUT*** the prior is really describing the underlying distribution of orbital periods, eccentricies, etc. So we can use our data not only to inform the properties of each system, but to map out the demographic distribution of exoplanet periods and more! Pretty cool! (e.g. https://arxiv.org/abs/1406.3020)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian credible regions <a class=\"anchor\" id=\"three\"></a>\n",
    "\n",
    "In the **frequentist paradigm**, the meaning of the *confidence interval* $\\mu_0 \\pm \\sigma_{\\mu}$ is \n",
    "the interval that would contain the true $\\mu$ (from which the data were drawn) in $68\\%$ (or $X\\%$) cases\n",
    "of a large number of *imaginary repeated experiments* (each with a different N values of $\\{x_i\\}$). \n",
    "\n",
    "However, the meaning of the so-called ***Bayesian credible region*** is *fundamentally different*: it is the interval that contains the true $\\mu$ with a probability of $68%\\$ (or $X\\%$), given the given dataset (our dear one and only dataset - there are no imaginary experiments in Bayesian paradigm). This credible region is the \n",
    "relevant quantity in the context of scientific measurements. \n",
    "\n",
    "There are several important features of a Bayesian posterior distribution:\n",
    "- They represent how well we believe a parameter is constrained within a certain range\n",
    "- We often quote the posterior maximum, or the **Maximum A Posteriori (MAP)**\n",
    "- We also often the posterior marginalized mean, $\\bar{\\theta} = \\int \\theta\\, p(\\theta|D)d\\theta$\n",
    "- ***The credible regions are not unique***. We can compute them in two different ways (visualized below)\n",
    "    1. We can integrate downwards from the MAP to enclose $X\\%$ (\"highest probability density interval\"), or\n",
    "    2. We can integrate inwards from each tail by $X/2\\%$ (\"equal-tailed interval\")\n",
    "\n",
    "<img src=\"figures/fig_credibleregion.png\" alt=\"\" style=\"width: 500px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple parameter estimation examples <a class=\"anchor\" id=\"four\"></a>\n",
    "\n",
    "### The first exercise in all Bayesian lectures: the coin flip \n",
    "\n",
    "Imagine you met a gambler who wanted you to **bet on whether the coin would come up heads or tails**. You watch the gambler work as a number of other people guess (as the coin is flipped). \n",
    "\n",
    "<font color='red'>Calculate and plot the Bayesian posterior of the probabilty of drawing a head as a function of watching [5, 50, 500] coin flips. Talk and work in your breakout rooms to do this. Share screen and co-code if necessary.</font>\n",
    "\n",
    "Remember the probability of getting $h$ heads in $n$ coin flips is given by the binomial probability distribution\n",
    ">$P(h|\\theta) = \\theta^h (1-\\theta)^{(n-h)}$\n",
    "\n",
    "with $\\theta$ the probability of a head.\n",
    "\n",
    "You will need to\n",
    "\n",
    "* formulate the likelihood, $p(D\\,|\\,M,\\theta,I)$\n",
    "* choose a prior$^1$, \n",
    "* plot the posterior pdf, $p(M,\\theta \\,|\\,D,I)$\n",
    " \n",
    " \n",
    " $^1$: Start with just a uniform prior on $\\theta$ that is constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell for some data to use\n",
    "data = np.array([1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1,\n",
    "       1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1,\n",
    "       1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1,\n",
    "       1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1,\n",
    "       1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0,\n",
    "       1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0,\n",
    "       1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1,\n",
    "       1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,\n",
    "       1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1,\n",
    "       1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0,\n",
    "       1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0,\n",
    "       0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1,\n",
    "       1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0,\n",
    "       1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1,\n",
    "       1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0,\n",
    "       0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0,\n",
    "       1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1,\n",
    "       0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1,\n",
    "       1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
    "       0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1,\n",
    "       1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0,\n",
    "       0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Modify and execute this cell.\n",
    "\n",
    "def p_nheads(n, h, theta):\n",
    "    # code the binomial distribution here\n",
    "    return ___\n",
    "\n",
    "# likelihood\n",
    "def L(data, ntrials, theta):\n",
    "    nheads = data[:ntrials].sum()\n",
    "    return p_nheads(ntrials, nheads, theta)\n",
    "\n",
    "# prior (ignore loc and scale until the beta prior)\n",
    "def prior(theta, loc=3.0, scale=3.0):\n",
    "    # put a prior distribution here\n",
    "    return ___\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "\n",
    "probs = []\n",
    "# define the number of trials\n",
    "n_trials = [5, 50, 500]\n",
    "theta = np.linspace(0, 1, 100)\n",
    "for nt in n_trials:\n",
    "    prob = L(data, nt, theta) * prior(theta)\n",
    "    # posterior plot\n",
    "    ax.plot(theta, prob, \n",
    "            label=r'ntrial = '+str(nt))\n",
    "    # likelihood plot\n",
    "    ax.plot(theta, L(data, nt, theta))\n",
    "    \n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if I used a $\\beta$ function for the prior?\n",
    "\n",
    "$P(\\alpha, \\beta, \\theta) = \\theta^{\\alpha-1} (1-\\theta)^{(\\beta-1)}$\n",
    "\n",
    "Modify your code above to use a $\\beta$ function prior on $\\theta$ with $a=3$ and $b=3$. You can code this directly or grab from `scipy.stats`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--> Paste and modify here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Nuisance parameters and marginalization\n",
    "\n",
    "#### Heteroscedastic Gaussian as an example\n",
    "\n",
    "\n",
    "Consider the case of measuring a rod.  We want to know the posterior pdf for the length of the rod, $p(M,\\theta|D,I) \\equiv p(\\mu|\\{x_i\\},\\{\\sigma_i\\},I)$.\n",
    "\n",
    "For the likelihood we have\n",
    "\n",
    "$$L = p(\\{x_i\\}|\\mu,I) = \\prod_{i=1}^N \\frac{1}{\\sigma_i\\sqrt{2\\pi}} \\exp\\left(\\frac{-(x_i-\\mu)^2}\n",
    "{2\\sigma_i^2}\\right).$$\n",
    "\n",
    "---\n",
    "\n",
    "**In the Bayesian case, we also need a prior.**  We'll adopt a *flat uniform distribution* given by\n",
    "\n",
    "$$p(\\mu|I) = C, \\; {\\rm for} \\; \\mu_{\\rm min} < \\mu < \\mu_{\\rm max},$$\n",
    "\n",
    "where $C = \\frac{1}{\\mu_{\\rm max} - \\mu_{\\rm min}}$ between the min and max and is $0$ otherwise.\n",
    "\n",
    "The log of the posterior pdf is then\n",
    "\n",
    "$$\\ln L = {\\rm constant} - \\sum_{i=1}^N \\frac{(x_i - \\mu)^2}{2\\sigma_i^2}.$$\n",
    "\n",
    "This is exactly the same as we saw before, except that the value of the constant is different.  Since the constant doesn't come into play, we get the same result as before:\n",
    " \n",
    "$$\\mu^0 = \\frac{\\sum_i^N (x_i/\\sigma_i^2)}{\\sum_i^N (1/\\sigma_i^2)},$$\n",
    "\n",
    "with uncertainty\n",
    "\n",
    "$$\\sigma_{\\mu} = \\left( \\sum_{i=1}^N \\frac{1}{\\sigma_i^2}\\right)^{-1/2}.$$\n",
    " \n",
    "We get the same result because we used a flat prior. If the case were homoscedastic instead of heteroscedastic, we obviously would get the result from our first example.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now let's consider the case where **$\\sigma$ is not known**, but rather it needs to be determined from the data, too.\n",
    "\n",
    "In this case, the posterior pdf that we seek is not $p(\\mu|\\{x_i\\},\\{\\sigma_i\\},I)$, but rather the joint $2$D pdf $p(\\mu,\\sigma|\\{x_i\\},I)$.\n",
    "\n",
    "As before we have\n",
    "\n",
    "$$L = p(\\{x_i\\}|\\mu,\\sigma,I) = \\prod_{i=1}^N \\frac{1}{\\sigma\\sqrt{2\\pi}} \\exp\\left(\\frac{-(x_i-\\mu)^2}{2\\sigma^2}\\right),$$\n",
    "\n",
    "except that now $\\sigma$ is uknown.\n",
    "\n",
    "---\n",
    "\n",
    "Our Bayesian prior is *now 2D instead of 1D* and we'll adopt \n",
    "\n",
    "$$p(\\mu,\\sigma|I) \\propto \\frac{1}{\\sigma},\\; {\\rm for} \\; \\mu_{\\rm min} < \\mu < \\mu_{\\rm max} \\; {\\rm and} \\; \\sigma_{\\rm min} < \\sigma < \\sigma_{\\rm max}.$$\n",
    "\n",
    "With proper normalization, we have\n",
    "\n",
    "$$p(\\{x_i\\}|\\mu,\\sigma,I)p(\\mu,\\sigma|I) = C\\frac{1}{\\sigma^{(N+1)}}\\prod_{i=1}^N \\exp\\left( \\frac{-(x_i-\\mu)^2}{2\\sigma^2}  \\right),$$\n",
    "\n",
    "where\n",
    "\n",
    "$$C = (2\\pi)^{-N/2}(\\mu_{\\rm max}-\\mu_{\\rm min})^{-1} \\left[\\ln \\left( \\frac{\\sigma_{\\rm max}}{\\sigma_{\\rm min}}\\right) \\right]^{-1}.$$\n",
    "\n",
    "The log of the posterior pdf is\n",
    "\n",
    "$$\\ln[p(\\mu,\\sigma|\\{x_i\\},I)] = {\\rm constant} - (N+1)\\ln\\sigma - \\sum_{i=1}^N \\frac{(x_i - \\mu)^2}{2\\sigma^2}.$$\n",
    "\n",
    "Since we are assuming a Gaussian distribution, the mean, $\\overline{x}$, and the variance, $V (=s^2)$, completely characterize the distribution. So we can write this expression in terms of those variables instead of $x_i$.  Skipping over the math details (see textbook $\\S$5.6.1), we find\n",
    "\n",
    "$$\\ln[p(\\mu,\\sigma|\\{x_i\\},I)] = {\\rm constant} - (N+1)\\ln\\sigma - \\frac{N}{2\\sigma^2}\\left( (\\overline{x}-\\mu)^2 + V  \\right).$$\n",
    "\n",
    "Note that this expression only contains the 2 parameters that we are trying to determine: $(\\mu,\\sigma)$ and 3 values that we can determine directly from the data: $(N,\\overline{x},V)$. A side note: these three data-based values fully encapsulate our dataset and are called *sufficient statistics*.\n",
    "\n",
    "<font color='red'>Execute the next cell to visualize the posterior pdf for the case of $(N,\\overline{x},V)=(10,1,4)$. Try playing around with and changing the values of $(N,\\overline{x},V)$.</font> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load code/fig_likelihood_gaussian.py\n",
    "\"\"\"\n",
    "Log-likelihood for Gaussian Distribution\n",
    "----------------------------------------\n",
    "Figure5.4\n",
    "An illustration of the logarithm of the posterior probability density\n",
    "function for :math:`\\mu` and :math:`\\sigma`, :math:`L_p(\\mu,\\sigma)`\n",
    "(see eq. 5.58) for data drawn from a Gaussian distribution and N = 10, x = 1,\n",
    "and V = 4. The maximum of :math:`L_p` is renormalized to 0, and color coded as\n",
    "shown in the legend. The maximum value of :math:`L_p` is at :math:`\\mu_0 = 1.0`\n",
    "and :math:`\\sigma_0 = 1.8`. The contours enclose the regions that contain\n",
    "0.683, 0.955, and 0.997 of the cumulative (integrated) posterior probability.\n",
    "\"\"\"\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "from astroML.plotting.mcmc import convert_to_stdev\n",
    "\n",
    "def gauss_logL(xbar, V, n, sigma, mu):\n",
    "    \"\"\"Equation 5.57: gaussian likelihood\"\"\"\n",
    "    return (-(n + 1) * np.log(sigma)\n",
    "            - 0.5 * n * ((xbar - mu) ** 2 + V) / sigma ** 2)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Define the grid and compute logL\n",
    "sigma = np.linspace(1, 5, 70)\n",
    "mu = np.linspace(-3, 5, 70)\n",
    "xbar = 1\n",
    "V = 4\n",
    "n = 10\n",
    "\n",
    "logL = gauss_logL(xbar, V, n, sigma[:, np.newaxis], mu)\n",
    "logL -= logL.max()\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot the results\n",
    "fig = plt.figure(figsize=(10, 7.5))\n",
    "plt.imshow(logL, origin='lower',\n",
    "           extent=(mu[0], mu[-1], sigma[0], sigma[-1]),\n",
    "           cmap=plt.cm.binary,\n",
    "           aspect='auto')\n",
    "plt.colorbar().set_label(r'$\\ln(L)$', fontsize=15)\n",
    "plt.clim(-5, 0)\n",
    "\n",
    "plt.contour(mu, sigma, convert_to_stdev(logL),\n",
    "            levels=(0.683, 0.955, 0.997),\n",
    "            colors='k')\n",
    "\n",
    "plt.text(0.5, 0.93, r'$L(\\mu,\\sigma)\\ \\mathrm{for}\\ \\bar{x}=1,\\ V=4,\\ n=10$',\n",
    "         bbox=dict(ec='k', fc='w', alpha=0.9),\n",
    "         ha='center', va='center', \n",
    "         fontsize=15, transform=plt.gca().transAxes)\n",
    "\n",
    "plt.xlabel(r'$\\mu$', fontsize=15)\n",
    "plt.ylabel(r'$\\sigma$', fontsize=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The plot from the previous cell is described by \n",
    "\n",
    "$$\\ln[p(\\mu,\\sigma|\\{x_i\\},I)] = {\\rm constant} - (N+1)\\ln\\sigma - \\frac{N}{2\\sigma^2}\\left( (\\overline{x}-\\mu)^2 + V  \\right).$$\n",
    "\n",
    "\n",
    "**What if we don't care that much about $\\sigma$?** e.g. it may simply represent how good our measuring apparatus is - what we might get a Nobel prize for is $\\mu$. \n",
    "- Therefore, what we want is $p(\\mu|\\{x_i\\},I)$, rather than $p(\\mu,\\sigma|\\{x_i\\},I)$. \n",
    "- We can get the former from the latter by the **marginalization over $\\sigma$**, which means integration of $p(\\mu,\\sigma|\\{x_i\\},I)$ over $\\sigma$: \n",
    "\n",
    "$$ p(\\mu\\,|\\,\\{x_i\\},I) = \\int_0^\\infty p(\\mu,\\sigma|\\{x_i\\},I) d\\sigma$$\n",
    "and thus (using the substitution $t$ = 1/$\\sigma$ and integration by parts)\n",
    "\n",
    "$$ p(\\mu\\,|\\,\\{x_i\\},I) \\propto \\left(1 + \\frac{(\\overline{x}-\\mu)^2}{V} \\right)^{-N/2}. $$\n",
    "\n",
    "It is easy to show that this result corresponds to *Student’s $t$ distribution* (Google it!) with \n",
    "$k = N-1$ degrees of freedom for the variable $t = (x-\\mu)/(s/\\sqrt{N})$, where $s$ is the sample \n",
    "standard deviation. As we've seen the Student’s $t$ distribution is symmetric and bell shaped, but with heavier \n",
    "tails than  a Gaussian distribution.  \n",
    "\n",
    "**Hold on!** This is not a Gaussian distribution promised by the Central Limit Theorem!!! \n",
    "\n",
    "That's fine! The CLT promised a Gaussian only for large $N$. Indeed, the above expression\n",
    "(Student's $t$ distribution) morphs into a Gaussian for large $N$ - let's check!\n",
    "\n",
    "<font color='red'>Let's see how Student's t distribution morphs into Gaussian distribution. let's first choose mean x and variance as in the previous example. </font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### let's see how Student's t distribution morphs into Gaussian distribution\n",
    "# let's first choose mean x and variance as in the previous example\n",
    "mu = np.linspace(-3, 5, 1000)\n",
    "xbar = 1.0\n",
    "V = 4.0\n",
    "\n",
    "# and now let's generate p(mu) for N=3, 10 and 20\n",
    "def compute_pmu(mu, xbar, V, N):\n",
    "    p = (1 + (xbar - mu) ** 2 / V) ** (-0.5 * N)\n",
    "    return p / np.sum(p)\n",
    "#\n",
    "N1 = 3\n",
    "pmu1 = compute_pmu(mu,xbar,V,N1)\n",
    "G1 = norm(xbar, np.sqrt(V/N1)).pdf(mu) \n",
    "G1norm = G1 / np.sum(G1)\n",
    "\n",
    "N2 = 10\n",
    "pmu2 = compute_pmu(mu,xbar,V,N2)\n",
    "G2 = norm(xbar, np.sqrt(V/N2)).pdf(mu) \n",
    "G2norm = G2 / np.sum(G2)\n",
    "\n",
    "N3 = 30\n",
    "pmu3 = compute_pmu(mu,xbar,V,N3)\n",
    "G3 = norm(xbar, np.sqrt(V/N3)).pdf(mu) \n",
    "G3norm = G3 / np.sum(G3)\n",
    " \n",
    "## plot\n",
    "fig, ax = plt.subplots(figsize=(8, 5.75))\n",
    "plt.plot(mu, pmu1, ls='-', c='blue', label=r'$N=3$')\n",
    "plt.plot(mu, G1norm, ls='--', c='blue', label=r'Gauss')\n",
    "plt.plot(mu, pmu2, ls='-', c='red', label=r'$N=10$')\n",
    "plt.plot(mu, G2norm, ls='--', c='red', label=r'Gauss')\n",
    "plt.plot(mu, pmu3, ls='-', c='black', label=r'$N=30$')\n",
    "plt.plot(mu, G3norm, ls='--', c='black', label=r'Gauss')\n",
    "\n",
    "plt.xlim(-3, 5)\n",
    "plt.ylim(0, 0.01)\n",
    "plt.xlabel('$\\mu$', fontsize=15)\n",
    "plt.ylabel(r'$p(\\mu | \\{x_i\\},\\sigma)$', fontsize=15)\n",
    "plt.legend()\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The top panels of the figure below (textbook figure 5.5, for the code see \n",
    "[astroML website](http://www.astroml.org/book_figures/chapter5/fig_posterior_gaussian.html)) show marginal distributions $p(\\mu)$ and $p(\\sigma)$.  The solid line is \n",
    "analytic Bayesian result with uninformative prior and the dotted line is the result for a uniform prior (note \n",
    "that there is not that much difference).  The dashed line is the MLE result, which is quite different.  The bottom panels show the cumulative distributions.\n",
    "\n",
    "![Ivezic, Figure 5.5](http://www.astroml.org/_images/fig_posterior_gaussian_1.png)\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- The main result here is that **for smallish N ($<$10 or so), $p(\\mu)$ is not Gaussian!** \n",
    "- The code above can be used to compute $p(\\mu)$ for arbitrary values of N, $\\overline{x}$ and V.\n",
    "- For large N, Gaussian is a good approximation of $p(\\mu)$. "
   ]
  }
 ],
 "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
}