08_Gaussian_Approx.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaussian approximations and model comparison\n",
    "\n",
    "The previous notebooks have introduced a variety of methods for performing **parameter inference** using complex models. All these methods depend on **Bayes' Theorem** to estimate something proportional to the **posterior probability density**.\n",
    "\n",
    "We've also covered a variety of ways for searching complex posteriors, ranging from simple Monte Carlo approaches like importance sampling and GLUE, to basic MCMC (e.g. Metropolis) and some more sophisticated MCMC algorithms such as the affine invariant sampler in [emcee](http://dan.iel.fm/emcee/current/).\n",
    "\n",
    "So far, we've only focussed on inferring parameter distributions for a **single model**, but what if we have several different models that we'd like to compare? This notebook introduces some of the basic ideas behind **Bayesian model comparison**. At the same time, we'll take a step back from the earlier considerations of sampling schemes (MC, MCMC etc.) to take a look at the **Gaussian approximation** approach to model calibration.\n",
    "\n",
    "The Gaussian approximation method is both simpler and more computationally efficient than the methods discussed in earlier notebooks. It also provides a very complete approach to the inference process, allowing both parameter & uncertainty estimation *and* model comparison (as we will see, some of the more complex methods introduced earlier can struggle when it comes to model comparison). \n",
    "\n",
    "The Gaussian approximation method is only valid in particular circumstances, but you may be surprised by just how often you can get away with using it. Because of it's power and simplicity, **this approach is well worth considering** for an initial attempt at the inference process - especially if you're planning to compare multiple model structures.\n",
    "\n",
    "\n",
    "## Work in progress\n",
    "\n",
    "### Plan\n",
    "\n",
    "* Ability to approximate the posterior in the vicinity of the MAP by a multi-dimensional Gaussian <br><br>\n",
    "* Justification for why this often works OK for large $N$ (link to **Central Limit Theorem**?) <br><br>\n",
    "* Estimates of uncertainty from parameter covariance matrix (from **Hessian**?) <br><br>\n",
    "* Introduce model comparison and the idea that the **normalising constant** from a parameter inference step becomes the **model likelihood** in the model comparison step <br><br>\n",
    "* Describe how the Bayesian method automatically penalises \"**over-fitting**\" <br><br>\n",
    "* Estimate the \"**evidence**\" for a model from **determinant** of **Hessian** <br><br>\n",
    "* Remember to stress that for model comparison we need **proper** priors - see page 354 of [Information Theory, Inference and Learning Algorithms](http://www.inference.phy.cam.ac.uk/itila/book.html) <br><br>\n",
    "* As long as the Gaussian approximation is appropriate, provides very efficient and comprehensive method of finding:\n",
    "  1. The optimum parameter set for each model under consideration <br><br>\n",
    "  2. The uncertainty associated with the parameter optimum for each model <br><br>\n",
    "  3. The relative likelihoods of each model <br><br>  \n",
    "* Only requires a single optimisation run; no very intensive numerical sampling (unless required for optimisation) <br><br>\n",
    "* Good place to start! <br><br>\n",
    "* How do we know if the Gaussian approximation is valid/appropriate?\n",
    "\n",
    "Notes so far are based on Chapter 28 of David MacKay's superb book **[Information Theory, Inference and Learning Algorithms](http://www.inference.phy.cam.ac.uk/itila/book.html)**. \n",
    "\n",
    "Still need to produce some working code\n",
    "\n",
    "## 1. A simple Bayesian workflow\n",
    "\n",
    "Suppose we have several different models for a particular process or system of interest. We have a single observed dataset that we will use to calibrate each of the models. Some of the models are more complex (e.g. have more poorly-constrained parameters) than others. We are interested in:\n",
    "\n",
    "  1. **Calibrating each model**. This involves identifying plausible parameter sets and associated **credible intervals**. <br><br>\n",
    "  \n",
    "  2. **Comparing the models**. Are the more complex models *significantly* better at simulating the data than the simple ones? Is the additional complexity \"worth it\"? <br><br>\n",
    "  \n",
    "  3. **Making predictions**. What is our best estimate for the next observation? Should we select the \"best\" model and use it, or can we integrate output from all the models, weighted according to their relative probability?\n",
    "\n",
    "We have already considered the problem of **model calibration** (i.e. inferring posterior parameter distributions) in some detail. However, in the discussions so far, we have generally ignored the normalising constant in Bayes' equation by writing \n",
    "\n",
    "$$posterior \\propto likelihood \\times prior$$\n",
    "\n",
    "instead of the full version\n",
    "\n",
    "$$posterior = \\frac{likelihood \\times prior}{normalising \\; constant}$$\n",
    "\n",
    "Recall from [notebook 3](http://nbviewer.ipython.org/github/JamesSample/enviro_mod_notes/blob/master/notebooks/03_Monte_Carlo.ipynb#1.1.-A-Bayesian-model-calibration-framework) that the normalising constant is sometimes referred to as the \"**probability of the data**\" or the \"**evidence**\".\n",
    "\n",
    "If we're just interested in parameter inference for a single model, it is OK to ignore the normalising constant. However, if we want to compare several competing models, this constant becomes important. Below, we'll consider a slightly more general case of model calibration where more than one model is available.\n",
    "\n",
    "### 1.1. Parameter inference\n",
    "\n",
    "If we assume that a particular model, $M_i$, is correct, we can write Bayes' equation as\n",
    "\n",
    "$$P(\\theta|D, M_i) = \\frac{P(D|\\theta, M_i)P(\\theta|M_i)}{P(D|M_i)}$$\n",
    "\n",
    "where $\\theta$ is the vector of parameters for model $M_i$ and $D$ is the data (i.e. the observations) used for calibration.\n",
    "\n",
    "If all of the model parameters, $\\theta$, are continuous variables, we can re-write the denominator in this equation as\n",
    "\n",
    "$$P(D|M_i) = \\int_\\theta{P(D|\\theta, M_i)P(\\theta|M_i) d\\theta}$$\n",
    "\n",
    "(See the end of [notebook 1](http://nbviewer.ipython.org/github/JamesSample/enviro_mod_notes/blob/master/notebooks/01_Distributions.ipynb#1.4.-Marginal-and-conditional-distributions) for an explanation of where these equations come from).\n",
    "\n",
    "To calibrate a single model to a dataset, we are primarily interested in finding two things:\n",
    "\n",
    "  1. Some kind of **central estimate** for an appropriate parameter set (e.g. the median, mean or \"best\" parameter set), and <br><br>\n",
    "  \n",
    "  2. An indication of the **uncertainty** associated with the estimate for each parameter.\n",
    "  \n",
    "There are many ways of achieving these goals, depending on the complexity of the problem. In the previous notebooks we have taken a very comprehensive approach, by using MC or MCMC techniques to estimate **marginal posterior distributions** for each parameter. From these distributions we can calculate a whole range of different statistics for the central estimate (e.g. the MAP) and the uncertainty bounds (e.g. a 95% credible interval). The truly Bayesian approach, however, is simply to report the distribution (i.e. a histogram) in its entirety.\n",
    "\n",
    "As we have seen, though, estimating marginal posteriors is not always easy and can be computationally expensive. In some cases, simpler approaches may provide adequate answers much more quickly.\n",
    "\n",
    "#### 1.1.1. Gaussian approximation\n",
    "\n",
    "For the simplest possible example of a one parameter model with an approximately Gaussian posterior distribution, we could simply report the **mean** and the **variance** to achieve the two aims listed above. Most real problems are more complex, but for some multi-parameter models it may still be reasonable to assume an approximately Gaussian posterior in the region around the most promising parameter set. In this case, the two aims listed above can be achieved by:\n",
    "\n",
    "  1. Finding the **[Maximum a posteriori (MAP)](https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation)** estimate. This is just the maximum of the posterior distribution, so finding it is an optimisation problem rather than a sampling problem, so it is likely to be faster (see examples in notebooks 2 and 6). <br><br>\n",
    "  \n",
    "  2. Obtain **confidence intervals** for each parameter by considering the **local curvature** of the posterior distribution in the region surrounding the MAP. This can be done by evaluating the **[Hessian](https://en.wikipedia.org/wiki/Hessian_matrix)**, which can provide an estimate of the parameter **covariance matrix**.\n",
    "  \n",
    "As long as the posterior can be reasonably approximated by a multi-dimensional Gaussian in the vicinity of the MAP, reporting the **location of the MAP** and the **parameter covariance matrix** is simply the multi-dimensional equivalent of reporting the mean and variance in the 1D example mentioned above. This approach therefore provides a very efficient (though less generally applicable) method of model auto-calibration and uncertainty estimation.\n",
    "  \n",
    "Note also that, as the amount of data collected increases, the Gaussian approximation can often be surprisingly effective. **(Link to the Central Limit Theorem?)**\n",
    "\n",
    "### 1.2. Model comparison\n",
    "\n",
    "For a discrete set of models, $M_i$, we can write another version of Bayes' equation\n",
    "\n",
    "$$P(M_i|D) = \\frac{P(D|M_i)P(M_i)}{P(D)}$$\n",
    "\n",
    "In this case, because the set of models is discrete, we can re-write the denominator as a sum, rather than as an integral like we did above\n",
    "\n",
    "$$P(D) = \\sum_iP(D|M_i)P(M_i)$$\n",
    "\n",
    "As we have done previosuly, in many cases it is usual to ignore the \"evidence\" in this equation and just write\n",
    "\n",
    "$$P(M_i|D) \\propto P(D|M_i)P(M_i)$$\n",
    "\n",
    "The prior probability for each model, $P(M_i)$, represents our initial belief about how likely model $M_i$ is to be correct. If we have no reason to prefer some models over others, $P(M_i)$ will be the same for all models (i.e. a uniform prior). \n",
    "\n",
    "The likelihood term, $P(D|M_i)$, encapsulates what the data can tell us about the probability of $M_i$ being correct. Note from section 1.1 above that this term is exactly the same as the **evidence** (i.e. the denominator) in the form of Bayes' equation used for parameter inference. This is generally the case: **the \"evidence\" from the parameter inference stage becomes a crucial component of Bayes' equation at the model comparison stage**. \n",
    "\n",
    "In previous notebooks, we have generally ignored the normalising constant in Bayes' equation when performing **parameter inference**. This is fine as long as you're only interested in evaluating parameter-related uncertainty for a **single model**. However, as you can see, it is necessary to include it if you're interested in comparing multiple models.\n",
    "\n",
    "#### 1.2.1 Occam's razor\n",
    "\n",
    "**[Occam's razor](https://en.wikipedia.org/wiki/Occam%27s_razor)** embodies the intuitive principle that, if two models provide equally good explanations of the data, the simplest model should be preferred. This principle is automatically incorporated into Bayesian model comparison in a quantitative way, which is extremely useful. This is because simple models have less flexibility and therefore make a narrower range of predictions than complex models. Because the total volume under the probabiltiy distribution $P(D|M_i)$ must equal $1$, the probabiltiy density for more complex models tends to be more \"spread out\" (and therefore lower) than for simple models. ***If*** a simple model manages to match the data just as well, it is therefore likely to have a higher probability density in the region of credible parameter sets than a more complex model with a wider range of output. The simpler model will therefore \"win\" during model comparison.\n",
    "\n",
    "#### 1.2.2. Gaussian approximation\n",
    "\n",
    "If the posterior can be well approximated by a Gaussian in the vicinity of the MAP (which, perhaps surprisingly, is often the case given enough data), then the **evidence** for model $M_i$, $P(D|M_i)$, can be estimated using the **best-fit likelihood** and **prior** and the **determinant** of the corresponding covariance matrix (which is obtained from the **Hessian**). See [equation 28.10](http://www.inference.phy.cam.ac.uk/itprnn/book.pdf) in David MacKay's book for further details.\n",
    "\n",
    "For the case where the Gaussian approximation is valid, a full Bayesian model evaluation can therefore be performed in a surprisingly efficient way, without requiring any complicated and computationally intensive numerical sampling (e.g. MCMC). First, for each model:\n",
    "\n",
    "  1. Find the MAP. <br><br>\n",
    "  2. Use the Hessian evaluated at the MAP to find the parameter covariance matrix.\n",
    " \n",
    "These two steps provide, for each model, an estimate of the \"best\" parameter set and an indication of the parameter-related uncertainty. Then, for model comparison:\n",
    "\n",
    "  1. Estimate the **evidence** for each model using the best-fit likelihood, $P(D|\\theta_{MAP}, M_i)$, the best-fit prior, $P(\\theta_{MAP}, M_i)$ and the **determinant** of the Hessian (from above). <br><br>\n",
    "  2. Multiply the evidence by the prior for each model and use the result to rank the models according their probability. Predictions can be made by calculating a weighted sum of model predictions, with the weights based on the posetrior probabiltiy of each model.\n",
    " \n",
    "Note that this workflow uses only a single **optimisation** step - no complex sampling algorithms are required. As long as the Gaussian approximation is valid, it is therefore an extremely efficient way to perform model calibration, uncertainty analysis and comparison."
   ]
  }
 ],
 "metadata": {
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}