MeanNormal_MCMC_Exercise9.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inferring the Parameters of a Normal with Markov Chain Monte Carlo (MCMC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are interested in learning about the trait values in a population (of salamanders!). We begin by making the assumption that the trait values in our population are normally distributed. We then sample some trait values from the population and try to infer the mean of the overall population based on the sample.\n",
    "\n",
    "First, we need to define how big our sample size will be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sampleSize <- 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to sample this many values from the population. In our case, we will know the true mean (and standard deviation) for our population, because we get to pick them! In `RevBayes`, functions that draw values from distributions all begin with __r__, because these are functions that draw **r**andom values. In our case, we want to draw values from a normal distribution, so we use `rnorm()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trueMean <- 10                                     # Define mean trait value in the population\n",
    "trueStandDev <- 1                                  # Define the standard deviation of trait values in the population\n",
    "data <- rnorm(sampleSize,trueMean,trueStandDev)    # Draw our sample from the population\n",
    "data                                               # Print our sample to the screen for inspection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have some data, we need to assign prior probability distributions to our model parameters and figure out how to calculate the likelihood. For our purposes, we will make the assumption that our true mean has uniform probability of being between 0 and 50, while the standard deviation has a uniform probability of being between 0 and 4. In this toy case, we will not specify how we came to these decisions, but for analyses of empirical data we would need to use outside information to set reasonable limits on the range of possible values for these parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Mean has a Uniform(0,50) prior\n",
    "# Standard deviation has a Uniform (0,4) prior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to give these variables starting values, to define a starting place for our exploration of parameter space. To do so, we'll draw random values from their prior distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m <- rUniform(1,0,50)  # Draw 1 value from a Uniform(0,50) distribution\n",
    "sd <- rUniform(1,0,4)  # Draw 1 value from a Uniform(0,4) distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print these values, just to see where we'll be starting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data\n",
    "m[1]  # Even though m has just one value, it's stored in a vector so we need to extract it\n",
    "sd[1] # The same for the standard deviation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need a way to calculate the likelihood, $P(D|H)$. In this case, our hypothesis ($H$) is a Normal distribution with a given mean and standard deviation. $D$ represents the data, the values sampled from our \"true\" Normal distribution. To do this, we'll create our own function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This function loops over values in our dataset and returns the log-likelihood of the data.\n",
    "# clear(calcLike) # Use this to remove previous function definition if defining again.\n",
    "function calcLike (Real[] dat, Real mean, Real stdev){\n",
    "    logLike <- 0.0\n",
    "    for (d in dat){\n",
    "        logLike += dnorm(d,mean,stdev)\n",
    "    }\n",
    "    return logLike\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate the likelihood for the starting values of our mean and standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "likes <- [calcLike(data,m[1],sd[1])]\n",
    "likes[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll begin sampling using the Metropolis-Hastings algorithm. For the proposal distribution of our sampler, we'll use a Normal distribution. Remember, this Normal is _not_ related to the distribution of trait values in the population, it is simply a function that gives us a way to propose moves through parameter space. \n",
    "\n",
    "Recall that the general outline of each Metropolis-Hastings step is:\n",
    "\n",
    "- Propose a new set of parameter values\n",
    "- Compare the posterior density of the proposed position to the current position\n",
    "- Decide whether to accept the proposal based on the comparison \n",
    "- Add the chosen values to the chain\n",
    "\n",
    "Because of the bug in the Jupyter notebook, I can't add comments in the `for` loop below. But you should copy the code to a text editor and add your own pseudocode to label the steps in the Metropolis-Hastings algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numGens = 3000\n",
    "for (gen in 1:numGens){\n",
    "    \n",
    "    if (gen % 500 == 0){\n",
    "        print(gen)\n",
    "    }\n",
    "    \n",
    "    currLike <- calcLike(data,m[gen],sd[gen])\n",
    "    \n",
    "    propMean <- rnorm(1,m[gen],1)\n",
    "    propSD <- rnorm(1,sd[gen],1)\n",
    "    propLike <- calcLike(data,propMean[1],propSD[1])\n",
    "    \n",
    "    LR <- (propLike-currLike)\n",
    "    \n",
    "    priorR <- (1/50)\n",
    "    if (propSD[1] < 0) { priorR <- 0.0 }\n",
    "    if (propSD[1] > 50) { priorR <- 0.0 }\n",
    "    if (propMean[1] < 0) { priorR <- 0.0 }\n",
    "    if (propMean[1] > 50) { priorR <- 0.0 }\n",
    "    \n",
    "    ranNumber <- runif(1,0,1)\n",
    "    \n",
    "    if (exp(LR)*priorR > 1){\n",
    "        likes.append(propLike)\n",
    "        m.append(propMean)\n",
    "        sd.append(propSD)\n",
    "    } else if (exp(LR)*priorR > ranNumber[1]){\n",
    "        likes.append(propLike)\n",
    "        m.append(propMean)\n",
    "        sd.append(propSD)\n",
    "    } else {\n",
    "        likes.append(currLike)\n",
    "        m.append(m[gen])\n",
    "        sd.append(sd[gen])\n",
    "    }\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "burnin = 1000\n",
    "\n",
    "posteriorMean = 0.0\n",
    "for (i in burnin:numGens){\n",
    "    posteriorMean += m[i]\n",
    "}\n",
    "posteriorMean = posteriorMean/(numGens-burnin)\n",
    "posteriorMean\n",
    "\n",
    "posteriorSD = 0.0\n",
    "for (i in burnin:numGens){\n",
    "    posteriorSD += sd[i]\n",
    "}\n",
    "posteriorSD = posteriorSD/(numGens-burnin)\n",
    "posteriorSD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, let's try plotting the values sampled by our chain in R. Select all the values sampled for m (click once at the beginning of the list and then hold down shift while clicking once at the end of the list), and then for sd. Use the R template provided to visualize the chain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "likes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explore the behavior of the sampler as you change:\n",
    "\n",
    "- The size of the proposal distributions\n",
    "- The number of generations\n",
    "- The amount of burn-in\n",
    "\n",
    "While the starting parameter values are chosen at random, how does the behavior of the chain vary with different starting values?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Use as needed to clear your RevBayes workspace and start again\n",
    "clear()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "RevBayes",
   "language": "bash",
   "name": "revbayes_kernel"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "r"
   },
   "file_extension": ".Rev",
   "help_links": [
    {
     "text": "RevBayes",
     "url": "https://revbayes.org"
    },
    {
     "text": "RevBayes Kernel",
     "url": "https://github.com/sdwfrost/revbayes_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-rsrc",
   "name": "RevBayes",
   "pygments_lexer": "R"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}