Monte Carlo exercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Play with the MCMC demonstration\n",
    "\n",
    "Import mcmc_demo.R (in the zip file) into R Studio (or whatever you use to run R) and run it. Experiment with changing the initialization (thetaInit), the step size (sampleSig), and the number of MCMC samples (Nsamp). Note that the resulting histogram of the samples is written out as a PDF file to your working directory.\n",
    "\n",
    "Explanation of the mcmc_demon.R code:\n",
    "\n",
    "The code samples the test function {\\tt testfunc}, or any other function you choose to provide.  I recommend you go through the code in chunks rather then sourcing the entire file.  With {\\tt plotev} set to {\\tt TRUE} the code will introduce a delay of {\\tt delay} seconds between each iteration.  (Actually it's just the plotting that is decelerated; all the samples are calculated in advance by {\\tt metrop}.)  At each iteration the code plots the current value of the chain in blue, and then the proposed candidate value in red and also prints this candidate value to the console. If the proposal is accepted it turns green. If it is not accepted it stays red and the current value turns green, because it is accepted again. The value to be used as the next current value is then printed to the console. The next iteration starts by turning the current point blue. The final plot can either overplot the individual samples or a histogram (density estimate) of the set of samples. The former is only suitable if {\\tt Nsamp} is less than a few hundred. I had to fiddle around a bit to get the histogram to appear right (see the comments in the code).  I suggest you experiment with different initial values, {\\tt thetaInt}, and different step sizes, {\\tt sampleSig}.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MCMC sampling of the signal/background estimation problem\n",
    "\n",
    "Write an R program to use the metrop function (in the file metropolis.R, in the zip file) to sample the 2D posterior from the signal/background estimation problem (in the R program notebook \"parameter estimation - multiple parameters\". It's also the R progam signal_background_estimation.R in the zip file). Use the data simulated at the beginning of signal_background_estimation.R. \n",
    "\n",
    "I recommend you follow the structure of the file signal_background_estimation.R. Copy it and delete everything which comes after the simulation of the data. Then use the syntax in the MCMC demontration (mcmc_demo.R) to work out how you can use the metropolis function to sample the posterior (logupost). Note that the function metrop requires the first argument of logupost to be a vector containing all the parameters we sample over. Therefore you will need to edit logupost to pass the two scalars a and b as a single vector (call it theta). You also need to edit logupost to return a two element vector (for compatibility with metrop). Just set the first value to 0. So if logupost is currently evaluating/returning the quantity which I'll now call \"lp\", you need to modify it to return \n",
    "\n",
    "c(0, lp)/log(10))\n",
    "\n",
    "(the factor 1/log(10) is also for compatibility, as metrop expects the log base 10, not the natural log, of the probability density).\n",
    "\n",
    "Think of sensible step sizes for the proposal distribution (i.e. values of sampleCov), as well as the initialization of the Markov Chain (thetaInit). \n",
    "\n",
    "The result of the Metropolis algorithm is a matrix with four columns: 0, log10(posterior density), a, b, with one row for each step of the Markov chain. Plot the resulting chains (a vs. row no., b vs. row no.). Do the chains look reasonable? Also plot a vs. b, which is the 2D posterior PDF. Use kernel density estimation (or a histogram) to plot the marginal (1D) posterior PDFs. Compute the mean, variances, and covariance of the 2D posterior. Note: once we've got a set of samples from the posterior, we only need these to compute everything we need!\n",
    "\n",
    "Experiment with changing the initialization of the Markov chain and its step sizes to see what impact this has on your results. You can also experiment with changing the data (by changing the true parameters of the model used to simulate the data) and investigate how well you can infer the parameters of the model.\n",
    "\n",
    "If you've managed all this, you can experiment with using a different prior. You could incorporate the prior into the logupost function. Remember that the prior does not need to be normalized."
   ]
  },
  {
   "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
}