parameter estimation - line fitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory\n",
    "\n",
    "Suppose we have a two-dimensional set of $N$ points $\\{x_i,y_i\\}$.\n",
    "The model $M$ predicts the values of $y$ as being \n",
    "\\begin{alignat}{2}\n",
    "y \\,&=\\, f(x) + \\epsilon \\ \\ \\ \\ \\ {\\rm where} \\\\\n",
    "f(x) \\,&=\\, b_0 + b_1x\n",
    "\\end{alignat}\n",
    "which is a straight line with parameters $b_0$ (intercept) and $b_1$ (gradient). \n",
    "$f(x)$ is the generative model: it gives the noise-free predictions of the data given the parameters.\n",
    "The residuals\n",
    "$\\epsilon = y - f(x)$ are modelled as a zero-mean Gaussian random variable with standard deviation $\\sigma$, i.e.\\\n",
    "$\\epsilon \\sim {\\cal N}(0, \\sigma)$.\n",
    "This is the \n",
    "{\\em noise model}.\n",
    "Assuming the $\\{x\\}$ are noise free, this tells us that the likelihood is\n",
    "\\begin{equation}\n",
    "P(y_i  |  x_i, \\theta, M) \\,=\\, \\frac{1}{\\sigma\\sqrt{2 \\pi}} \\exp{ \\left[ -\\frac{[y_i - f(x_i; b_0, b_1)]^2}{2\\sigma^2} \\right] } \n",
    "\\end{equation}\n",
    "where $\\theta = (b_0, b_1, \\sigma)$ are the parameters of the model.\n",
    "It may come as a surprise that we will try to infer the uncertainty in the data points from the data. Yet $\\sigma$ is as a model parameter just like the others.\n",
    "Although the $x$ values are supplied with the data, we assume them to be fixed: they are not described by a measurement model.\n",
    "Thus the data are $D=\\{y_i\\}$. Assuming that the various $y$ measurements are independent, the log likelihood for all $N$ data points is\n",
    "\\begin{equation}\n",
    "\\ln P(\\{y_i\\}  |  \\{x_i\\}, \\theta, M) \\,=\\, \\sum_{i=1}^N \\ln P(y_i  |  x_i, \\theta, M) \\ .\n",
    "\\end{equation}\n",
    "In general none of the three parameters are known in advance, so we want to infer their posterior PDF from the data, which is given as usual by\n",
    "\\begin{equation}\n",
    "P(\\theta  |  D) \\,\\propto\\, P(D  | \\theta) P(\\theta) \\ .\n",
    "\\end{equation}\n",
    "As we will be sampling from the posterior we do not need to compute the normalization constant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Procedure\n",
    "\n",
    "Given a set of data, the procedure to compute the posterior is as follows:\n",
    "* define the prior PDF over the parameters. I will use plausible yet convenient priors, and I will make use of a variable transformation;\n",
    "* define the covariance matrix of the proposal distribution. I will use a diagonal, multivariate Gaussian;\n",
    "* define the starting point (initialization) of the MCMC;\n",
    "* define the number of burn-in iterations and the number of sampling iterations.\n",
    "\n",
    "Once we have run the MCMC we perform the following analyses:\n",
    "* thin the chains;\n",
    "* plot the chains and the one-dimensional marginal posterior PDFs over the parameters. I do the latter via kernel density estimation;\n",
    "* plot the two-dimensional posterior distributions of all three pairs of parameters, simply by plotting the samples (we could do two-dimensional kernel density estimation instead). I do this to look for correlations between the parameters;\n",
    "* calculate the maximum a posteriori (MAP) values of the model parameters from the MCMC chains, calculate and plot the resulting model, and compare to the original data;\n",
    "* calculate the predictive posterior distribution over $y$ at a new data point \n",
    "Note that because we have samples drawn from the posterior, we don't need the actual values of the posterior density in order to plot the posteriors. We likewise don't have to do any integration to get the one-dimensional marginal distributions.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code\n",
    "\n",
    "The analysis described above is all done with the R script linearmodel_posterior.R, with explanations provided as comments in the code. The code looks long, but a reasonable chunk of it is actually concerned with plotting and analysing the results.  The code initially sources two other files. The first file, linearmodel_functions.R, defines functions that compute the (logarithm of the) prior, likelihood, and posterior. The second file is the Metropolis algorithm metrop in metropolis.R. The rest of the code is then executed. Part of the code at the end of the file is concerned with doing prediction, but we are not covering that in this course (so ignore or delete that part of the program)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "To better appreciate how MCMC works in this example, I experiment with the code. Change in particular the following (one by one):\n",
    "* the initialization of the parameters;\n",
    "* the parameter step sizes;\n",
    "* the number of iterations and the length of the burn-in;\n",
    "* the value of the standard deviation of the prior on $b_0$;\n",
    "* the prior on $b_0$ from a Gaussian to an improper uniform prior. Do this by setting {\\tt b0Prior} to unity in the function {\\tt logprior.linearmodel};\n",
    "* the amount of data. Try both more data points (e.g.\\ 100) and fewer, including just 3, 2, and 1 data points.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Things to think about and to investigate\n",
    "\n",
    "* Why is the inferred straight line sometimes quite different from the true straight line?\n",
    "* How can we infer the errors bars as well as the line?\n",
    "* What happens if we have much more data?\n",
    "* How are the results affected by the choice of priors?\n",
    "* What happens if we reduce the size of the data set to two points, or even just one?\n",
    "* How, logically thinking, can we estimate three parameters at all given only one or two data points? (Isn't the the solution is somehow \"underdetermined''?)\n",
    "* What would the result be if we had no data? (Don't try to run my code with no data, however, as it's not robust to this)\n",
    "* How would you extend the code to fit higher order polynomials, or even different types of functions?\n",
    "\n",
    "If your current work involve fitting curves to data, think about taking this Bayesian approach of determining the posterior PDF over the parameters. It has several advantages over simply using linear least squares: we get full posteriors on the parameters, rather than simplified symmetric error bars; we can take into account prior information; it can easily be extended to include uncertainty on both axes, upper/lower limits, outliers, etc. You don't have to implement all of this by hand. Packages exist in R, and other programming languages, whi\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.3.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}