S03a_Simulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute before using this notebook if using google colab\n",
    "\n",
    "kernel = str(get_ipython())\n",
    "\n",
    "if 'google.colab' in kernel:    \n",
    "    !wget https://raw.githubusercontent.com/fredzett/rmqa/master/utils.py -P local_modules -nc \n",
    "    !npx degit fredzett/rmqa/data data\n",
    "    import sys\n",
    "    sys.path.append('local_modules')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from utils import plot_bar, plot_line, plot_hist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why is it so useful and important for statistics and research?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulation is a very useful tool for statistics / data analysis and research in general. Some reasons for this are:\n",
    "\n",
    "- **solve probability problems**: use numerical approaches rather than analytical approaches\n",
    "\n",
    "- **calculate features of probability densities**: calculating features of probability densities involves integration. Instead we could draw a randomly from that density and just calculate the feature from the sampled data. \n",
    "\n",
    "    > Example:\n",
    "    > - what is the variance of a random variable that is described by a $\\text{Beta}(\\alpha=3,\\beta=6)$ distribution?\n",
    "    > - if we don't know the formula for this we can simply draw from this distribution (e.g. 1.000 times)\n",
    "    > - we then calculate the variance of our randomly sampled data\n",
    "  \n",
    "\n",
    "- **understand the assumed data generation process**: simulate data from a proposed model to get a feel for the data\n",
    "   \n",
    "\n",
    "- **evaluating estimators**: estimator $\\approx$ function of the data. How do we know if model is well specified and appropriate for our model? We can simulate a world, randomly draw data from that world and apply our function to that data. We can then check if we get the right estimator. \n",
    "\n",
    "    > Example:\n",
    "    > - create world: $y = 5 + 17 x + \\epsilon$, where: $\\epsilon = N(0,\\sigma^2)$\n",
    "    > - randomly draw data from world: $x = \\text{1.000 random numbers}$ $\\implies y $\n",
    "    > - estimate linear regression function and check if parameters $\\beta_0$ and $\\beta_1$ can be recovered / are what we expect (i.e. if a linear regression model is well specified for our data)\n",
    "\n",
    "\n",
    "**Overall:** simulation is a very useful and powerful tool to solve complex problems without using complex math. This is especially true when dealing with probability given humans do not have a good intuition for probabilities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to simulation with `numpy`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In `python` there are numerous ways to simulate data.\n",
    "\n",
    "Let's look at some easy simulation examples using `numpy`. Specificall, we need to use:\n",
    "\n",
    "- `np.random`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Randomly draw one number between 1 and 6 \n",
    "np.random.randint(1,7) # note: upper bound is exclusive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Randomly draw 1.000 numbers between 1 and 6\n",
    "data = np.random.randint(1,7,size=1_000)\n",
    "data[:10] # show first 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at a typical highschool / undergradiuate probability problem and solve it via simulation:\n",
    "\n",
    "a. You roll one dice $8x$. What is the probability of $4$ sixes?\n",
    "\n",
    "b. You roll two dices once. What is the probability that the sum of eyes is 7?\n",
    "\n",
    "c. You make three draws (with replacement) from a Urn holding 20 red, 30 blue, 8 green balls. What is the probability of getting no blue ball?\n",
    "\n",
    "You can solve this problems mathematically (i.e. analytically). However, we will use simulation to answer the questions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a. What is the probability of 4 sixes out of 8 rolls?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "rolls = 8\n",
    "data = np.random.randint(1,7,size=(n,rolls)) # draw n samples of 8 rolls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sixes = np.sum(data == 6, axis=1) # determine the number of sixes in each sample (i.e. each 8 rolls)\n",
    "sixes[:10] # show first 10 samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exactly_four = sum(sixes == 4) / n  # determine number of sixes that are equal to 4 and divide by number of samples (to get probability)\n",
    "exactly_four"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With higher $n$ the probability will converge to the analytical correct probability which is:\n",
    "\n",
    "$$P(X = 4) = \\binom{8}{4}\\frac{1}{6}^4(1-\\frac{1}{6})^{8-4} \\approx 2.6048\\%$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**b. What is the probability of two dices summing to 7?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "dice = 2\n",
    "data = np.random.randint(1,7,size=(n,dice))\n",
    "data[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_eyes = np.sum(data,axis=1)\n",
    "sum_eyes[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(sum_eyes == 7) / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With higher $n$ the probability will converge to the analytical correct probability. Here two dices sum up to $7$ with the following combinations:\n",
    "\n",
    "- 1, 6\n",
    "- 6, 1\n",
    "- 2, 5\n",
    "- 5, 2\n",
    "- 3, 4\n",
    "- 4, 3\n",
    "\n",
    "with are $6$ combinations out of $36$ total possibilities. Therefore the analytical correct solution is $\\frac{6}{36} \\approx 16.7\\%$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What is probability of no blue ball with three draws?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we had assumed equally probabilities for each possible event. Now we need to introduce different probabilities of a specific event occuring\n",
    "\n",
    "- red = 20/58\n",
    "- blue = 30/58\n",
    "- green = 8/58\n",
    "\n",
    "We can model this easily using `numpy.random` by passing an additional argument `p` (for probability)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100000\n",
    "draws = 3\n",
    "probs = [20/58, 30/58, 8/58] # With replacement\n",
    "data = np.random.choice([\"red\",\"blue\",\"green\"],size=(n,draws),p=probs)\n",
    "data[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_blues = np.sum(data == \"blue\",axis=1)\n",
    "sum(num_blues == 0) / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With higher $n$ the probability will converge to the analytical correct probability. Here we have $58^3$ possible combinations to draw. Drawing no blue ball has possible combinations of $28^3$. Therefore the analytical solution is:\n",
    "\n",
    "$$P(\\text{no blue}) = \\frac{28^3}{58^3} \\approx  11.25\\%$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve probability problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at different and more relevant problems. \n",
    "\n",
    "Assume that you want to buy a used book from an online retailer. The book is offered from a resellers with the following reviews:\n",
    "\n",
    "- Seller A: 100% positive ratings; 10 total ratings\n",
    "\n",
    "- Seller B: 96% positive ratings; 50 total ratings\n",
    "\n",
    "The % of positive rating observed is likely not to be the real positive rating %. \n",
    "\n",
    "Assuming the real positive ratings are 95% for both sellers\n",
    "\n",
    "- how likely are the above rating observations?\n",
    "\n",
    "- plot a histogram with the distribution of positivity ratings (use: plot_hist() for this)\n",
    "\n",
    "\n",
    "(Note: this example is inspired by [John Cook](https://www.johndcook.com/blog/2011/09/27/bayesian-amazon/) and [3Blue1Brown](https://www.youtube.com/watch?v=8idr1WZ1A7Q&feature=youtu.be))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100000\n",
    "s = 0.95\n",
    "ks = [10,50]\n",
    "obs = [10,48]\n",
    "for k,o in zip(ks, obs):\n",
    "    data = np.random.choice([1,0],p=[s,1-s],size=(n,k))\n",
    "    sums = np.sum(data,axis=1)\n",
    "    print(\"Probability of observing\", o,\"/\",k,\" = \", sum(sums==o) / n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_hist(sums);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Monthe Hall Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem is loosely based on the American TV (game) show *Let's Make a Deal* and is named after its original host, Monty Hall. \n",
    "\n",
    "The problem became very famous as a quesiton from a reader's letter quoted in Marilyn vos Savant's \"Ask Marilyn\" column in a US magazine in 1990. \n",
    "\n",
    "The question was as follows:\n",
    "\n",
    "> Suppose you're on a game show, and you're given the choice of three doors: Behind one door is a car; behind the others, goats. You pick a door, say No. 1, and the host, who knows what's behind the doors, opens another door, say No. 3, which has a goat. He then says to you, \"Do you want to pick door No. 2?\" Is it to your advantage to switch your choice?\n",
    "\n",
    "### The Problem\n",
    "\n",
    "There are 3 doors in front of you, and there is a prize behind one of them (a brand new car!).\n",
    "Once you select a door, the moderator will open one of the two you had not selected which does not have a prize behind it.\n",
    "You will then have the opportunity to switch from the door you originally selected to an alternate door.\n",
    "\n",
    "Here are some examples\n",
    "\n",
    "| Door 1 | Door 2 | Door 3 | Monte Hall opens... | Result from Switch | Result from stay |\n",
    "|:------:|:------:|:------:|:------------------------:|:------------------:|:----------------:|\n",
    "|   Car  |  Goat  |  Goat  |        Door 2 or 3       |        Goat        |        Car       |\n",
    "|  Goat  |   Car  |  Goat  |          Door 3          |         Car        |       Goat       |\n",
    "|  Goat  |  Goat  |   Car  |          Door 2          |         Car        |       Goat       |\n",
    "\n",
    "**Question** Is it  better to switch or to stay with your initial choice? (Assuming you want to win the car!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem became so famous because Vos Savant's response was that the contestant should switch to the other door because switching would give you a probability of $\\frac{2}{3}$ of winning the car. Many famous mathematicians disputed the result but it turns out that under the above assumptions the answer is actually correct. See [here](https://en.wikipedia.org/wiki/Monty_Hall_problem) for more infos regarding the problem. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solution to Monte Hall Problem using simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 1000\n",
    "win_no_switch = 0\n",
    "win_switch = 0\n",
    "\n",
    "for i in range(n):\n",
    "    doors = [1,2,3]\n",
    "    win_door = np.random.choice(doors)\n",
    "    choice = np.random.choice(doors)\n",
    "    if win_door == choice: # 1/3\n",
    "        win_no_switch += 1\n",
    "    doors_left = doors[:]\n",
    "    doors_left.remove(choice)\n",
    "    if any(win_door == doors_left):\n",
    "        win_switch += 1\n",
    "print(\"Prob(Car | Switch) = \", win_switch/n)\n",
    "print(\"Prob(Car | No Switch) = \", win_no_switch/n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Birthday Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The birthday problem or birthday paradox concerns the probability that, in a set of $n$ randomly chosen people, some pair of them will have the same birthday. The problem is a good example that assessing probabilities is often difficult and results are often non-intuitive. \n",
    "\n",
    "**Question**: Given a room with 23 randomly selected people, what is the probability that at least two have the same birthday?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solution to Birthday Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_birthdays(n,group_size):\n",
    "    same_day = 0\n",
    "    alldays = np.arange(1,366,1)\n",
    "    for i in range(n):\n",
    "        bdays = np.random.choice(alldays,size=group_size,replace=True)\n",
    "        if len(bdays) != len(set(bdays)):\n",
    "            same_day += 1\n",
    "    return same_day / n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100000\n",
    "group_size = 23\n",
    "simulate_birthdays(n,group_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 1000\n",
    "ks = 75\n",
    "prob_by_group_size = [simulate_birthdays(n,s) for s in np.arange(2,ks+1,1)]\n",
    "fig, ax = plot_line(np.arange(2,ks+1),[prob_by_group_size], zero_origin=True, xlabel=\"Group size\",ylabel=\"Probability of at least two birthdays on same day\");"
   ]
  }
 ],
 "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.8.3"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}