IntroToProbability_Exercise2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2 (Frequentist Probability)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are going to change our _sample space_, so that the outcomes have a more interesting biological meaning. In this case, we are going to imagine that we are sampling individual organisms from some population (say, a population of salamanders). Individuals can have four possible phenotypes:\n",
    "\n",
    "1. large and striped\n",
    "2. large and spotted\n",
    "3. small and striped\n",
    "4. small and spotted\n",
    "\n",
    "We will say that the `population` vector holds all individuals in our population. While, in this case, we could simply tally the phenotypes of every individual in this vector, we will not do that here in order to maintain our thought experiment. In almost any real biological population, we would not know precisely how many individuals it contained and we could not survey all of their phenotypes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "population = [\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"large_striped\",\"small_striped\",\"small_striped\",\"small_striped\",\"small_striped\",\"small_striped\",\"small_striped\",\"small_striped\",\"small_striped\",\"small_striped\",\"small_striped\",\"small_striped\",\"large_spotted\",\"large_spotted\",\"large_spotted\",\"large_spotted\",\"large_spotted\",\"large_spotted\",\"large_spotted\",\"large_spotted\",\"large_spotted\",\"large_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\",\"small_spotted\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate what the frequencies of the four phenotypes might be in our population, we will draw random samples. To do this, we can use the same RevBayes syntax as above, with one new skill. As we draw samples, we will want to tally the number of times we've sampled different phenotypes. To do that, we will need a way to \"ask\" RevBayes which phenotype we have in hand. \n",
    "\n",
    "The standard way to test a condition, which is used in nearly every programming language, is an `if...else` statement. In these statements, we include an argument to `if` that we use to test the condition we are interested in. Often, we are evaluating whether two values are equal to one another. To do that we use the equality operater, `==`. Here are some examples of the operator on its own."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"a\" == \"a\"\n",
    "\"a\" == \"b\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the values produced by `==` are either `TRUE` or `FALSE` - and __only__ one of those two. What `if` is actually doing is returning a _Boolean_ variable. These are variables used to control the flow of code by evaluating logical statements. In the case of an `if...else` statement, the code controlled by `if` is executed when the _Boolean_ is `TRUE`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if (TRUE){\n",
    "    print(\"Code for if statement.\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that above, we used a fixed value of `TRUE` for our `if` statement. In most cases, we'll want to replace that with a logical test involving some variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "myVar = 1\n",
    "if (myVar == 1){\n",
    "    print(\"myVar has value 1.\")\n",
    "} else {\n",
    "    print(\"myVar does not have value 1.\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, now we have all the tools we need to sample from our population and keep track of the results. The first thing we'll need to do is to set up a variable to hold the number of times we've sampled each phenotype. In this case, we'll use a vector (`phenoCounts`) that will have four integers. Each value stores the number of times we've sampled one of the four phenotypes listed above. Then, we'll build a `for` loop to sample from our population and store the counts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for (i in 1:4){\n",
    "    phenoCounts[i] = 0\n",
    "}\n",
    "\n",
    "numSamples = 20\n",
    "for (s in 1:numSamples){\n",
    "    sample = population[runifInt(1,1,population.size())[1]]\n",
    "    if (sample == \"large_striped\"){\n",
    "        phenoCounts[1] += 1\n",
    "    } else if (sample == \"small_striped\"){\n",
    "        phenoCounts[2] += 1\n",
    "    } else if (sample == \"large_spotted\"){\n",
    "        phenoCounts[3] += 1\n",
    "    } else if (sample == \"small_spotted\"){\n",
    "        phenoCounts[4] += 1\n",
    "    }\n",
    "}\n",
    "print(phenoCounts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can transform the counts to frequencies by dividing the entire vector by the number of samples we've drawn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "phenoCounts/numSamples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, now we'll run this whole sampling procedure with 20 samples per replicate 10 times so that we can see how consistent our estimates are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for (rep in 1:10){\n",
    "    for (i in 1:4){\n",
    "        phenoCounts[i] = 0\n",
    "    }\n",
    "\n",
    "    numSamples = 20\n",
    "    for (s in 1:numSamples){\n",
    "        sample = population[runifInt(1,1,population.size())[1]]\n",
    "        if (sample == \"large_striped\"){\n",
    "            phenoCounts[1] += 1\n",
    "        } else if (sample == \"small_striped\"){\n",
    "            phenoCounts[2] += 1\n",
    "        } else if (sample == \"large_spotted\"){\n",
    "            phenoCounts[3] += 1\n",
    "        } else if (sample == \"small_spotted\"){\n",
    "            phenoCounts[4] += 1\n",
    "        }\n",
    "    }\n",
    "    print(phenoCounts/numSamples)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How variable are each of the frequencies across replicates? Now we'll run this same procedure again, but drawing 100 samples per replicate. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for (rep in 1:10){\n",
    "    for (i in 1:4){\n",
    "        phenoCounts[i] = 0\n",
    "    }\n",
    "\n",
    "    numSamples = 100\n",
    "    for (s in 1:numSamples){\n",
    "        sample = population[runifInt(1,1,population.size())[1]]\n",
    "        if (sample == \"large_striped\"){\n",
    "            phenoCounts[1] += 1\n",
    "        } else if (sample == \"small_striped\"){\n",
    "            phenoCounts[2] += 1\n",
    "        } else if (sample == \"large_spotted\"){\n",
    "            phenoCounts[3] += 1\n",
    "        } else if (sample == \"small_spotted\"){\n",
    "            phenoCounts[4] += 1\n",
    "        }\n",
    "    }\n",
    "    print(phenoCounts/numSamples)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the __frequentist__ interpretation of probability, the probability of an event (like sampling a large, striped individual) is defined by its long-run frequency. The long-run frequency is the frequency of an event over a very large number of trials. So, the frequencies above (from 100 samples) are estimates of each event's true probability. However, because we have only drawn a finite number of samples (e.g., 100), our estimates have some error. The fact that our estimates become more precise as we draw more samples is known as the [Law of Large Numbers](https://en.wikipedia.org/wiki/Law_of_large_numbers). Let's run this one more time with 1,000 samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for (i in 1:4){\n",
    "    phenoCounts[i] = 0\n",
    "}\n",
    "\n",
    "numSamples = 1000\n",
    "for (s in 1:numSamples){\n",
    "    sample = population[runifInt(1,1,population.size())[1]]\n",
    "    if (sample == \"large_striped\"){\n",
    "        phenoCounts[1] += 1\n",
    "    } else if (sample == \"small_striped\"){\n",
    "        phenoCounts[2] += 1\n",
    "    } else if (sample == \"large_spotted\"){\n",
    "        phenoCounts[3] += 1\n",
    "    } else if (sample == \"small_spotted\"){\n",
    "        phenoCounts[4] += 1\n",
    "    }\n",
    "}\n",
    "print(phenoCounts/numSamples)"
   ]
  }
 ],
 "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
}