09_Summary.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    "This notebook attempts to provide a practical summary of everything covered previously, particularly with regard to choosing an appropriate workflow. The full range of methods available for model calibration, uncertainty estimation and comparison is huge, so this is by no means a comprehensive summary - it's simply an attempt to place the methods discussed earlier into context.\n",
    "\n",
    "## 1. Choosing an approach\n",
    "\n",
    "The flow chart below may be helpful when selecting an approach to take for your modelling assessment.\n",
    "\n",
    "<img src=\"https://github.com/JamesSample/enviro_mod_notes/blob/master/images/Bayesian_Workflow.png?raw=true\" alt=\"Bayesian workflow\" width=\"600\"> \n",
    "\n",
    "Regardless of whether you're performing single- or multi-model analysis, arguably the most difficult stages in either of the above workflows are:\n",
    "\n",
    "  1. Choosing an appropriate **likelihood function** and <br><br>\n",
    "  \n",
    "  2. Choosing an **efficient** method for estimating the **marginal posteriors** (or some summary statistics to represent them).\n",
    "\n",
    "## 2. Choosing a likelihood function\n",
    "\n",
    "Throughout this set of notes we have focused on **formal** statistical likelihoods, based on careful consideration of the expected error structure between the model and the observed data. If you don't know for certain what the error structure should be (and, in reality, you almost never do), it's usually best to start off with something simple and then progress in stages to more complex functions. At each stage, simple **diagnostic plots of the residuals** will show whether the assumptions have been satisfied, and may also suggest what modifications to try next.\n",
    "\n",
    "### 2.1.Simple independent Gaussian errors\n",
    "\n",
    "The simplest likelihood functions use independent and identically distributed (iid) Gaussian errors, such as those presented in notebooks [2](http://nbviewer.ipython.org/github/JamesSample/enviro_mod_notes/blob/master/notebooks/02_Calibration_Likelihood.ipynb#2.1.-The-likelihood-function) and [4](http://nbviewer.ipython.org/github/JamesSample/enviro_mod_notes/blob/master/notebooks/04_MCMC.ipynb#7.3.-Define-likelihood). These are rarely adequate for real-world environmental modelling, but it's a good place to start.\n",
    "\n",
    "When using these likelihoods with real data, residual plots will often show **[heteroscedasticity](https://en.wikipedia.org/wiki/Heteroscedasticity)**, which imples that the errors are *not* identically distributed (e.g. in hydrology it's common to find bigger errors at high flows). Autocorrelation plots may also show **serial dependence**, in which case the errors are not independent either. This latter problem is very common for models that simulate **time series** at fine-grained temporal resolution.\n",
    "\n",
    "### 2.2. Weighted errors\n",
    "\n",
    "One simple approach to deal with problems of heteroscedasticity is to use **weighted errors**. For example, the error variance may be made proportional to a parameter such as flow, exactly as we did with the hydrological model in [notebook 6](http://nbviewer.ipython.org/github/JamesSample/enviro_mod_notes/blob/master/notebooks/06_Beyond_Metropolis.ipynb#2.-Choosing-a-likelihood-function). An alternative way of achieving essentially the same thing is to consider **transforming** the data - for example by taking the **logarithm** or **square root** of the simulated and observed data series and performing the calibration using the transformed datasets instead.\n",
    "\n",
    "Simple weighted errors schemes often perform surprisingly well and are relatively easy to implement. \n",
    "\n",
    "### 2.3. Autoregressive errors\n",
    "\n",
    "If the residuals show autocorrelation, it may be worth considering an [autoregressive model](https://en.wikipedia.org/wiki/Autoregressive_model). An **AR(1) scheme** is often sufficient. We have not considered autoregressive models in these notebooks and there are some subtleties to be considered - see [Evin *et al.* (2013)](http://onlinelibrary.wiley.com/doi/10.1002/wrcr.20284/full) and [Evin *et al.* (2014)](http://onlinelibrary.wiley.com/doi/10.1002/2013WR014185/abstract) for further details. Such schemes are nevertheless well-established and not too difficult to code.\n",
    "\n",
    "## 3. Estimating marginal posteriors for the parameters\n",
    "\n",
    "Perhaps the most efficient approach for obtaining parameter estimates and confidence intervals for your model is to use the **Gaussian approximation** from [notebook 8](http://nbviewer.ipython.org/github/JamesSample/enviro_mod_notes/blob/master/notebooks/08_Gaussian_Approx.ipynb). This requires running an optimiser to find the [MAP](https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation) and then estimating the **covariance structure** in this vicinity using the **Hessian matrix**. This method has the additional advantage that it can be used to approximate the \"**evidence**\" in Bayes' equation, so it can also be used in a **model comparison** context. However, the technique is only valid *if* the posterior can be reasonably approximated by a multi-dimensional Gaussian, which is not always the case.\n",
    "\n",
    "* **Add something on how to test for this??**\n",
    "\n",
    "For more complicated posteriors, it may be necessary to use strategies based on **random sampling** or **random walks**. Of these **stochastic** approaches, the simplest are **Monte Carlo** techniques such as **Importance** or **Rejection** sampling (see [notebook 3](http://nbviewer.ipython.org/github/JamesSample/enviro_mod_notes/blob/master/notebooks/03_Monte_Carlo.ipynb#2.2.-Monte-Carlo-(MC)). Monte Carlo techniques are expected to be very inefficient in high dimensional parameter spaces, however, so for complex models **Markov chain Monte Carlo (MCMC)** is likely to be a much better option. The simplest MCMC approach uses the **Metropolis algorithm**, introduced in [notebook 4](http://nbviewer.ipython.org/github/JamesSample/enviro_mod_notes/blob/master/notebooks/04_MCMC.ipynb#3.-The-Metropolis-algorithm). More sophisticated MCMC methods include: **adaptive** algorithms (e.g. [DREAM](https://www.researchgate.net/profile/James_Hyman/publication/41035177_Accelerating_Markov_chain_Monte_Carlo_simulation_by_differential_evolution_with_self-adaptive_randomized_subspace_sampling/links/09e4150c0a96b6d034000000.pdf), which is very popular in hydrology); **affine invariant** algorithms, such as the one provided by [emcee](http://dan.iel.fm/emcee/current/) and used in [notebook 6](http://nbviewer.ipython.org/github/JamesSample/enviro_mod_notes/blob/master/notebooks/06_Beyond_Metropolis.ipynb); and algorithms which make use of gradient information, such as e.g. [Hamiltonian Monte Carlo](https://en.wikipedia.org/wiki/Hybrid_Monte_Carlo).\n",
    "\n",
    "All these approaches may be useful in different circumstances, but it is worth noting many of them - especially standard MC technqiues - may struggle to evaluate the evidence in complex models and are therefore **not necessarily well-suited to model comparison**. This is because the evidence, $P(D|M_i)$, is often a very computationally intensive integral, which is why we ignored it for the *parameter-inference-only* examples in previous notebooks\n",
    "\n",
    "$$P(D|M_i) = \\int_\\theta{P(D|\\theta, M_i)P(\\theta|M_i) d\\theta}$$\n",
    "\n",
    "A variety of methods exist for estimating this integral, but one of the most promising makes use of an analogy to statistical physics and is called **[thermodynamic integration](http://users.wpi.edu/~balnan/MarLik-BF-0.pdf)**. The details of this approach are beyond the scope of these notes, but it's worth pointing out that the emcee **[parallel-tempered ensemble sampler](http://dan.iel.fm/emcee/current/user/pt/)** includes a `thermodynamic_integration_log_evidence` method. In brief, parallel-tempering involves initialising multiple versions of the standrad affine invariant sampler, but with each one running at a different \"temperature\". Higher temperatures correspond to larger, more energetic steps in the Markov chains, in exactly the same way that molecules in a gas become more energetic at higher temperatures. Bigger steps mean chains running at high temperatures are more likely to move \"downhill\" into lower probability regions, whereas chains at lower temperatures focus their attemtion on exploring the \"peaks\" of the posterior distribution. At very low temperatures, the chains essentially become optimisers, in the sense that they will only ever move uphill.\n",
    "\n",
    "Running an ensemble of MCMC analyses at a range of temperatures provides a very robust exploration of the parameter space and can be very effective if the posterior is **multi-modal**. In addition, by computing the average of the log-likelihoods at various temperatures, it is possible to estimate the value of the **evidence integral** - emcee uses the method of [Goggans & Chi (2004)](http://scitation.aip.org/content/aip/proceeding/aipcp/10.1063/1.1751356).\n",
    "\n",
    "If your work involves comparing a number of complex models whose posteriors cannot be represented using the Gaussian approximation ([notebook 8](http://nbviewer.ipython.org/github/JamesSample/enviro_mod_notes/blob/master/notebooks/08_Gaussian_Approx.ipynb)), using the **parallel-tempered ensemble sampler** in emcee could be an excellent option."
   ]
  }
 ],
 "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
}