bacterial-colony.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How Long Will The Bacterial Colony Last?\n",
    "\n",
    "**2020-06-12**: https://fivethirtyeight.com/features/how-long-will-the-bacterial-colony-last/\n",
    "\n",
    "## Riddler Express\n",
    "\n",
    "From Dave Mattingly comes a high-rolling question of craps:\n",
    "\n",
    "There’s a technique for rolling dice called “bowling,” in which you place your index finger and thumb on two opposite sides of the die and roll it along the table. If done correctly, the die will never land on the faces on which you were holding the die, leaving you with a 25 percent chance of landing on each of the remaining four faces.\n",
    "\n",
    "You’d like to apply this technique to improve your chances of winning a simplified game of craps, in which your goal is simply to roll a 7 or 11 using two dice. With a standard rolling technique, your chances of rolling a 7 or 11 are 2/9, or about 22.2 percent.\n",
    "\n",
    "Now suppose you’re using your bowling technique, and you roll the dice one at a time (i.e., you know the outcome of the first die before rolling the second). If you play to maximize your chances of rolling a 7 or 11, what will be your chances of winning?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution\n",
    "\n",
    "There are only a small number of outcomes, so let's just enumerate them and add the probabilities. \n",
    "\n",
    "First, let's consider the possible finger configurations for \"bowling\" style rolls. There are only 3 different configurations:\n",
    "\n",
    "* Placing your fingers on the 1 and 6\n",
    "* Placing your fingers on the 2 and 5\n",
    "* Placing your fingers on the 3 and 6\n",
    "\n",
    "This leads to the following 4 sides being possible rolls in each scenario.\n",
    "\n",
    "```\n",
    "{2, 3, 4, 5}\n",
    "{1, 3, 4, 6}\n",
    "{1, 2, 5, 6}\n",
    "```\n",
    "\n",
    "Based on the outcome of the first roll, you'll simply select a configuration that allows you to reach both the 7 and 11 total, or a 7 if it's not possible to reach both. \n",
    "\n",
    "```\n",
    "1: {1, 3, 4, 6} or {1, 2, 5, 6} P(Win) = 25%\n",
    "2: {2, 3, 4, 5} or {1, 2, 5, 6} P(Win) = 25%\n",
    "3: {2, 3, 4, 5} or {1, 3, 4, 6} P(Win) = 25%\n",
    "4: {2, 3, 4, 5} or {1, 3, 4, 6} P(Win) = 25%\n",
    "5: {1, 2, 5, 6}                 P(Win) = 50%\n",
    "6: {1, 2, 5, 6}                 P(Win) = 50%\n",
    "```\n",
    "\n",
    "\n",
    "\n",
    "If your first die is a 5 or 6 it becomes possible to reach both a 7 and 11, giving you a 50% probability of winning with the correct \"bowling\" throw. As a result, for our first throw let's choose the `{1, 2, 5, 6}` configuration to give ourselves a 50% chance of getting one of those numbers. Now we have a 50% chance of entering into a state with a 25% chance to win and a 50% chance of entering into a state with a 50% chance of winning.\n",
    "\n",
    "\\begin{equation*}\n",
    "  \\frac{1}{2} \\times (\\frac{1}{2} + \\frac{1}{4} ) = \\frac{3}{8} = 37.5\\%\n",
    "\\end{equation*}\n",
    "\n",
    "Let's confirm the result through simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.375073"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "n_sims = 10**6\n",
    "\n",
    "# The bowling dictionary will take the result of the first roll as a key.\n",
    "# The value is the possible outcomes of the bowling configuration we choose as a result of the first roll.\n",
    "bowling = dict()\n",
    "bowling[1] = [1, 3, 4, 6]\n",
    "bowling[2] = [1, 2, 5, 6]\n",
    "bowling[3] = [2, 3, 4, 5]\n",
    "bowling[4] = [2, 3, 4, 5]\n",
    "bowling[5] = [1, 2, 5, 6]\n",
    "bowling[6] = [1, 2, 5, 6]\n",
    "\n",
    "first_rolls = np.random.choice([1, 2, 5, 6], size=n_sims)\n",
    "\n",
    "second_rolls = np.zeros(n_sims)\n",
    "\n",
    "for i, roll in enumerate(first_rolls):\n",
    "    second_rolls[i] = np.random.choice(bowling[roll])\n",
    "\n",
    "totals = first_rolls + second_rolls\n",
    "\n",
    "win_percentage = np.mean((totals == 7) | (totals == 11))\n",
    "win_percentage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This confirms our expected win percentage of `3/8`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Riddler Classic\n",
    "\n",
    "From Adam Wagner comes a curious case of colonies:\n",
    "\n",
    "You are studying a new strain of bacteria, Riddlerium classicum (or R. classicum, as the researchers call it). Each R. classicum bacterium will do one of two things: split into two copies of itself or die. There is an 80 percent chance of the former and a 20 percent chance of the latter.\n",
    "\n",
    "If you start with a single R. classicum bacterium, what is the probability that it will lead to an everlasting colony (i.e., the colony will theoretically persist for an infinite amount of time)?\n",
    "\n",
    "Extra credit: Suppose that, instead of 80 percent, each bacterium divides with probability `p`. Now what’s the probability that a single bacterium will lead to an everlasting colony?\n",
    "\n",
    "### Solution\n",
    "\n",
    "Unfortunately I've seen this type of question before, so I know it's very similar to [Bobo the Amoeba](https://www.quora.com/Bobo-the-amoeba-has-a-25-25-and-50-chance-of-producing-0-1-or-2-offspring-respectively-Each-of-Bobos-descendants-also-have-the-same-probabilities-What-is-the-probability-that-Bobos-lineage-dies-out) or the [Cave Escape](https://github.com/Booleans/math-through-simulation/blob/master/cave-escape.py). We simply set up an equation:\n",
    "\n",
    "\\begin{equation*}\n",
    "  D = .20 + .80D^2\n",
    "\\end{equation*}\n",
    "\n",
    "Where `D` represents the probability of a single bacterium dying off. It has a 20% chance of dying right away, but an 80% chance of transitioning to a state where 2 bacteria would have to die off. Expanding and solving the quadratic for `D` yields a 25% chance the bacterium dies out, meaning there's a 75% chance the colony lives on in perpetuity. Let's verify through simulation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.250062"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_sims = 10**6\n",
    "success = 0\n",
    "failure = 0\n",
    "\n",
    "for _ in range(n_sims):\n",
    "    num_amoebas = 1\n",
    "    while True:\n",
    "        num_amoebas = np.sum(np.random.choice(np.array([0, 2]), p=[.2, .8], size=(num_amoebas)))\n",
    "        # With such a high probability of reproduction I am choosing a threshold of 20 bacteria alive at once.\n",
    "        # Above this threshold I assume the colony lives on in perpetuity. \n",
    "        if num_amoebas > 20:\n",
    "            success += 1\n",
    "            break\n",
    "        elif num_amoebas == 0:\n",
    "            failure += 1\n",
    "            break\n",
    "            \n",
    "failure/n_sims"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This confirms that the colony dies out 25% of the time. "
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}