Lectures/lecture_4_notes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture #4: Bayesian versus Frequentist Inference\n",
    "## AM 207: Advanced Scientific Computing\n",
    "### Stochastic Methods for Data Analysis, Inference and Optimization\n",
    "### Fall, 2021"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"fig/logos.jpg\" style=\"height:150px;\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "### Import basic libraries\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import pandas as pd\n",
    "import sklearn as sk\n",
    "from sklearn.linear_model import LinearRegression\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Outline\n",
    "1. Review of Bayesian Modeling\n",
    "2. Examples of Conjugate and Non-Conjugate Models\n",
    "2. Connections to Frequentist Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Review of Bayesian Modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Bayesian Modeling Process\n",
    "\n",
    "In order to make statements about $Y$, the outcome, and $\\theta$, parameters of the distribution generating the data, we form the joint distribution over both variables and use the various marginals/conditional distributions to reason about $Y$ and $\\theta$.\n",
    "\n",
    "1. we form the ***joint distribution*** over both variables $p(Y, \\theta) = p(Y | \\theta) p(\\theta)$.\n",
    "\n",
    "2. we can condition on the observed outcome to make inferences about $\\theta$,\n",
    "$$\n",
    "p(\\theta | Y) = \\frac{p(Y, \\theta)}{p(Y)}\n",
    "$$\n",
    "where $p(\\theta | Y)$ is called the ***posterior distribution*** and $p(Y)$ is called the ***evidence***.\n",
    "3. before any data is observed, we can simulate data by using our prior\n",
    "$$\n",
    "p(Y^*) = \\int_\\Theta p(Y^*, \\theta) d\\theta = \\int_\\Theta p(Y^* | \\theta) p(\\theta) d\\theta\n",
    "$$\n",
    "where $Y^*$ represents new data and $p(Y^*)$ is called the ***prior predictive***.\n",
    "4. after observing data, we can simulate new data simliar to the observed data by using our posterior\n",
    "$$\n",
    "p(Y^*|Y) = \\int_\\Theta p(Y^*, \\theta|Y) d\\theta = \\int_\\Theta p(Y^* | \\theta) p(\\theta | Y) d\\theta\n",
    "$$\n",
    "where $Y^*$ represents new data and $p(Y^*|Y)$ is called the ***posterior predictive***."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayesian Update\n",
    "\n",
    "<img src=\"fig/bayes.jpg\" style=\"height:300px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Model Evaluation\n",
    "\n",
    "How do we know that our Bayesian model is a good fit for the data?\n",
    "\n",
    "1. **(Log-likelihood)** We can compute the marginal log-likelihood of the data under our posterior. That is, give a set of test data $\\{y^*_1, \\ldots, y^*_M\\}$, compute<br><br>\n",
    "\\begin{align}\n",
    "\\log \\prod_{m=1}^M p(\\mathbf{y}^*_m | \\text{Data}) &= \\sum_{m=1}^M \\log p(\\mathbf{y}^*_m | \\text{Data}) = \\sum_{m=1}^M \\log \\int_\\Theta p(\\mathbf{y}^*_m | \\theta) p(\\theta| \\text{Data}) d\\theta\n",
    "\\end{align}\n",
    "<br>This is simply called the **log-likelihood** of the data under the Bayesian model.<br><br>\n",
    "\n",
    "2. **(Posterior Predictive Check)** We can also compare the synthetic data generated from our posterior predictive:\n",
    "  1. sample from the posterior $\\theta_s \\sim p(\\theta | Y)$\n",
    "  2. plug the posterior samples into the likelihood, and sample synthetic data from the likelihood $Y_s \\sim p(Y | \\theta_s)$.\n",
    "  \n",
    "  We can then compare the synthetic data from the posterior predictive to the observed data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Components of Bayesian Inference\n",
    "\n",
    "We see that in order to evaluate Bayesian models we need to be able to perform two tasks: \n",
    "1. integration over the posterior (required by log-likelihood)\n",
    "2. sampling from the posterior (required by the posterior predictive check). \n",
    "\n",
    "Both requirements becomes easier if know the closed form expression for the posterior, e.g. if the prior is conjugate to the likelihood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Examples of Conjugate and Non-Conjugate Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayesian Model for (Univariate) Gaussian Likelihood with Known Variance\n",
    "### The Bayesian Model\n",
    "\n",
    "Let $Y \\sim \\mathcal{N}(\\mu, \\sigma^2)$, with $\\sigma^2$ known. We place a normal prior on $\\mu$, $\\mu\\sim\\mathcal{N}(m, s^2)$. \n",
    "\n",
    "**Question:** is our choice of prior appropriate?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayesian Model for (Univariate) Gaussian Likelihood with Known Variance\n",
    "**Inference:** The posterior $p(\\mu|Y)$ is then:\n",
    "\\begin{aligned}\n",
    "p(\\mu | Y) = \\frac{p(Y| \\mu)p(\\mu)}{p(Y)} = \\frac{\\overbrace{\\frac{1}{\\sqrt{2\\pi \\sigma^2}} \\mathrm{exp} \\left\\{-\\frac{(Y - \\mu)^2}{2\\sigma^2}\\right\\}}^{\\text{likelihood}} \\overbrace{\\frac{1}{\\sqrt{2\\pi s^2}} \\mathrm{exp} \\left\\{-\\frac{(m - \\mu)^2}{2s^2}\\right\\}}^{\\text{prior}}}{p(Y)}\n",
    "\\end{aligned}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "We can simplify the posterior as:\n",
    "\\begin{aligned}\n",
    "p(\\mu | Y) &= const *\\frac{\\mathrm{exp} \\left\\{ -\\frac{s^2(Y - \\mu)^2 + \\sigma^2(m - \\mu)^2}{2s^2\\sigma^2}\\right\\}}{p(Y)} \\\\\n",
    "&= const *\\mathrm{exp} \\left\\{ \\frac{s^2Y^2 + \\sigma^2m^2}{\\sigma^2 s^2}\\right\\}\\mathrm{exp} \\left\\{ -\\frac{(s^2 + \\sigma^2)\\mu^2 - 2(s^2Y + \\sigma^2m)\\mu}{2s^2\\sigma^2}\\right\\}\\\\\n",
    "&= const* \\mathrm{exp} \\left\\{ -\\frac{\\left(\\mu - \\frac{s^2Y + \\sigma^2m}{s^2 + \\sigma^2} \\right)^2}{2s^2\\sigma^2}\\right\\}\\quad \\text{(Completing the square)}\n",
    "\\end{aligned}\n",
    "Thus, we see that the posterior is a normal distribution, $\\mathcal{N}\\left(\\frac{s^2Y + \\sigma^2m}{s^2 + \\sigma^2}, s^2\\sigma^2\\right)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayesian Model for (Univariate) Gaussian Likelihood with Known Mean\n",
    "### The Bayesian Model\n",
    "\n",
    "Let $Y \\sim \\mathcal{N}(\\mu, \\sigma^2)$, with $\\mu$ known. We place an inverse-gamma prior on $\\sigma^2$, $\\sigma^2\\sim IG(\\alpha, \\beta)$. \n",
    "\n",
    "**Question:** is our choice of prior appropriate?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayesian Model for (Univariate) Gaussian Likelihood with Known Mean\n",
    "**Inference:** The posterior $p(\\sigma^2|Y)$ is then:\n",
    "\\begin{aligned}\n",
    "p(\\sigma^2 | Y) = \\frac{p(Y| \\sigma^2)p(\\sigma^2)}{p(Y)} = \\frac{\\overbrace{\\frac{1}{\\sqrt{2\\pi \\sigma^2}} \\mathrm{exp} \\left\\{-\\frac{(Y - \\mu)^2}{2\\sigma^2}\\right\\}}^{\\text{likelihood}} \\overbrace{\\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} \\left( \\sigma^2\\right)^{-\\alpha -1}\\mathrm{exp} \\left\\{-\\frac{\\beta}{\\sigma^2}\\right\\}}^{\\text{prior}}}{p(Y)}\n",
    "\\end{aligned}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can simplify the posterior as:\n",
    "\n",
    "\\begin{aligned}\n",
    "p(\\sigma^2 | Y) &= const * \\left( \\sigma^2\\right)^{-(\\alpha + 0.5) -1}\\mathrm{exp} \\left\\{-\\frac{\\frac{(Y-\\mu)^2}{2} + \\beta}{\\sigma^2}\\right\\}\n",
    "\\end{aligned}\n",
    "\n",
    "Thus, we see that the posterior is an inverse gamma distribution, $IG\\left(\\alpha + 0.5, \\frac{(Y-\\mu)^2}{2} + \\beta\\right)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayesian Model for (Univariate) Gaussian Likelihood with Unknown Mean and Variance\n",
    "\n",
    "Let $Y \\sim \\mathcal{N}(\\mu, \\sigma^2)$, with both parameters unknown. We place a normal prior on $\\mu$, $\\mu\\sim\\mathcal{N}(m, s^2)$, and an inverse-gamma prior on $\\sigma^2$, $\\sigma^2\\sim IG(\\alpha, \\beta)$.\n",
    "\n",
    "The posterior $p(\\sigma^2|Y)$ is then:\n",
    "\n",
    "\\begin{aligned}\n",
    "p(\\mu, \\sigma^2 | Y)  = \\frac{\\overbrace{\\frac{1}{\\sqrt{2\\pi \\sigma^2}} \\mathrm{exp} \\left\\{-\\frac{(Y - \\mu)^2}{2\\sigma^2}\\right\\}}^{\\text{likelihood}} \\overbrace{\\frac{1}{\\sqrt{2\\pi s^2}} \\mathrm{exp} \\left\\{-\\frac{(m - \\mu)^2}{2s^2}\\right\\}}^{\\text{prior on $\\mu$}}\\overbrace{\\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} \\left( \\sigma^2\\right)^{-\\alpha -1}\\mathrm{exp} \\left\\{-\\frac{\\beta}{\\sigma^2}\\right\\}}^{\\text{prior on $\\sigma^2$}}}{p(Y)}\n",
    "\\end{aligned}\n",
    "\n",
    "Can the posterior be simplified so that we recognize the form of the distribution?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Bayesian Model for Poisson Likelihood\n",
    "### The Bayesian Model\n",
    "\n",
    "With the kidney cancer dataset in mind, let $Y\\sim Poi(N\\theta)$, where $N$ is the total population and $\\theta$ is the underlying cancer rate. We place a gamma prior on $\\theta$, $\\theta \\sim Ga(\\alpha, \\beta)$. \n",
    "\n",
    "**Question:** is our choice of prior appropriate?\n",
    "\n",
    "### Inference\n",
    "The posterior $p(\\theta|Y)$ is then:\n",
    "\\begin{aligned}\n",
    "p(\\theta | Y) = \\frac{p(Y| \\theta)p(\\theta)}{p(Y)} = \\frac{\\overbrace{\\frac{(N\\theta)^Y}{Y!} \\mathrm{exp} \\left\\{-N\\theta\\right\\}}^{\\text{likelihood}} \\overbrace{\\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} \\left( \\theta\\right)^{\\alpha -1}\\mathrm{exp} \\left\\{-\\beta\\theta\\right\\}}^{\\text{prior}}}{p(Y)}\n",
    "\\end{aligned}\n",
    "\n",
    "We can simplify the posterior as:\n",
    "\n",
    "\\begin{aligned}\n",
    "p(\\theta | Y) &= const * \\theta^{\\alpha + Y - 1}\\mathrm{exp}\\left\\{-(N+ \\beta)\\theta \\right\\}\n",
    "\\end{aligned}\n",
    "\n",
    "Thus, we see that the posterior is a gamma distribution, $Ga\\left(\\alpha + Y, N + \\beta\\right)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Non-Conjugate Models\n",
    "\n",
    "We know that conjugate priors yield closed-form expressions for the posterior. In all our examples, these posteriors have been both easy to sample from and easy to integrate over. That is, we can evaluate our Bayesian models: can compute the log-likelihood of the data and perform posterior predictive checks.\n",
    "\n",
    "So why would we ever consider non-conjugate priors?\n",
    "\n",
    "Suppose that $Y\\sim \\mathcal{N}(\\mu, 2)$ represent the height of a person randomly selected from a population. Would the conjugate prior $\\mu\\sim \\mathcal{N}(5.7, 1)$ be appropriate for this application?\n",
    "\n",
    "If we wanted to use the prior $\\mu \\sim Ga(5.7, 1)$, would we be able to derive a closed form expression for the posterior?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Connections with Frequentist Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Point Estimates from the Posterior\n",
    "If you absolutely wanted to derive a point estimate for the parameters $\\theta$ in the likelihood from your Bayesian model, there are two common ways to do it:\n",
    "\n",
    "1. ***the posterior mean***\n",
    "$$\n",
    "\\theta_{\\text{post mean}} = \\mathbb{E}_{\\theta\\sim p(\\theta|Y)}\\left[ \\theta|Y \\right] = \\int \\theta p(\\theta|Y) d\\theta\n",
    "$$\n",
    "2. ***the posterior mode*** or ***maximum a posterior (MAP)*** estimate\n",
    "$$\n",
    "\\theta_{\\text{MAP}} = \\mathrm{argmax}_{\\theta} p(\\theta|Y)\n",
    "$$\n",
    "\n",
    "**Question:** is it better to summarize the entire posterior using a point estimate? I.e. why should we keep the posterior distribution around?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Point Estimates Can Be Misleading\n",
    "\n",
    "The posterior mode can be an atypical point: \n",
    "<img src=\"fig/map.jpg\" style=\"height:250px; width:450px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Point Estimates Can Be Misleading\n",
    "\n",
    "The posterior mean can be an unlikely point: \n",
    "\n",
    "<img src=\"fig/mean.jpg\" style=\"height:250px; width:450px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Comparison of Posterior Point Estimates and MLE\n",
    "\n",
    "**Beta-Binomial Model for Coin Flips**\n",
    " - Likelihood: $Bin(N, \\theta)$ \n",
    " - Prior: $Beta(\\alpha, \\beta)$\n",
    " - MLE: $\\frac{Y}{N}$\n",
    " - MAP: $\\frac{Y + \\alpha - 1}{N + \\alpha + \\beta - 2}$\n",
    " - Posterior Mean: $\\frac{Y + \\alpha}{N + \\alpha + \\beta}$ \n",
    "\n",
    "\n",
    "**Question:** What is the effect of the prior on the posterior point estimates? Imagine if $Y=10$, $N=11$, $\\alpha=100$, $\\beta=300$. What if $Y=1,000$, $N=11,000$, $\\alpha=1$, $\\beta=3$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Coin Toss Example: Revisited Yet Again\n",
    "Recall that one way to prevent the MLE from overfitting is to add ***regularization terms***:\n",
    "$$\n",
    "\\theta_{\\text{MLE Reg}} = \\frac{Y + \\alpha}{N + \\beta}.\n",
    "$$\n",
    "This is very similar to the MAP and posterior mean estimates:\n",
    " - MAP: $\\frac{Y + \\alpha - 1}{N + \\alpha + \\beta - 2}$\n",
    " - Posterior Mean: $\\frac{Y + \\alpha}{N + \\alpha + \\beta}$ \n",
    "\n",
    "In fact, we have see that one effect of adding a prior is that it **regularizes** our inference about $\\theta$.\n",
    "\n",
    "**Question:** What happens to the MAP and posterior mean estimates as $N$ (and hence $Y$) becomes very large?\n",
    "\n",
    "$$\n",
    "\\lim_{N\\to \\infty} \\frac{Y + \\alpha - 1}{N + \\alpha + \\beta - 2} = ?\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Law of Large Numbers for Bayesian Inference\n",
    "\n",
    "In general, in Bayesian inference we are **less interested asymptotic behavior**. But the properties of the asymptotic distribution of the posterior can be useful.\n",
    "\n",
    "**Theorem: (Berstein-von Mises**)\n",
    "\n",
    "\"Under some conditions, as $N\\to \\infty$ the posterior distribution converges to a Gaussian distribution centred at the MLE with covariance matrix given by a function of the Fisher information matrix at the true population parameter value.\"\n",
    "\n",
    "**Consequences**\n",
    "1. The posterior point estimates approach the MLE, with large samples sizes. <br><br>\n",
    "2. It may be valid to approximate the posterior with a Gaussian, with large samples sizes. This will become a very important idea during the second half of the course!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Computational Comparisons\n",
    "\n",
    "1. **Computation of the MLE is an optimization problem.** Although difficult, there are many established methods for performing optimization (even when the objective function is not convex -- i.e. many local optima). \n",
    "\n",
    "  More importantly, there are algorithms to perform general, automatic optimization (e.g. gradient descent) on a large class of functions -- that is, we do not need to artisanally solve an optimization problem for each statistical model.<br><br>\n",
    "  \n",
    "2. **Computation of the posterior (thus far) is an process of choosing the right priors and noting that the posterior distribution is of the same type as the prior.** The derivation is simple so long as we use conjugate priors. But many intuitively appropriate priors (like the inverse gamma and normal priors for a univariate gaussian) are not conjugate. In those cases, it becomes intractable to\n",
    " - compute posterior mean\n",
    " - simulate samples from the posterior (and hence simulate samples from the posterior predictive)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}