F0E-00-Probability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FLIP(00):  Advanced Data Science\n",
    "**(Module 03: A Concrete Introduction to Probability)**\n",
    "\n",
    "---\n",
    "- Materials in this module include resources collected from various open-source online repositories.\n",
    "- You are free to use, but NOT allowed to change or distribute this package.\n",
    "\n",
    "Prepared by and for \n",
    "**Student Members** |\n",
    "2006-2023 [TULIP Lab](http://www.tulip.org.au)\n",
    "\n",
    "---\n",
    "\n",
    "\n",
    "# Session 00 - A Concrete Introduction to Probability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "In 1814, Pierre-Simon Laplace [wrote](https://en.wikipedia.org/wiki/Classical_definition_of_probability):\n",
    "\n",
    ">*Probability theory is nothing but common sense reduced to calculation. ... [Probability] is thus simply a fraction whose numerator is the number of favorable cases and whose denominator is the number of all the cases possible ... when nothing leads us to expect that any one of these cases should occur more than any other.*\n",
    "\n",
    "![Laplace](https://upload.wikimedia.org/wikipedia/commons/thumb/3/30/AduC_197_Laplace_%28P.S.%2C_marquis_de%2C_1749-1827%29.JPG/180px-AduC_197_Laplace_%28P.S.%2C_marquis_de%2C_1749-1827%29.JPG)\n",
    "<center><a href=\"https://en.wikipedia.org/wiki/Pierre-Simon_Laplace\">Pierre-Simon Laplace</a><br>1814</center>\n",
    "\n",
    "\n",
    "Laplace nailed it. To untangle a probability problem, all you have to do is define exactly what the cases are, and careful count the favorable and total cases. Let's be clear on our vocabulary words:\n",
    "\n",
    "\n",
    "- **[Trial](https://en.wikipedia.org/wiki/Experiment_(probability_theory%29):**\n",
    "  A single occurrence with an outcome that is uncertain until we observe it. \n",
    "  <br>*For example, rolling a single die.*\n",
    "- **[Outcome](https://en.wikipedia.org/wiki/Outcome_(probability%29):**\n",
    "  A possible result of a trial; one particular state of the world. What Laplace calls a **case.**\n",
    "  <br>*For example:* `4`.\n",
    "- **[Sample Space](https://en.wikipedia.org/wiki/Sample_space):**\n",
    "  The set of all possible outcomes for the trial. \n",
    "  <br>*For example,* `{1, 2, 3, 4, 5, 6}`.\n",
    "- **[Event](https://en.wikipedia.org/wiki/Event_(probability_theory%29):**\n",
    "  A subset of the sample space, a set of outcomes that together have some property we are interested in.\n",
    "  <br>*For example, the event \"even die roll\" is the set of outcomes* `{2, 4, 6}`. \n",
    "- **[Probability](https://en.wikipedia.org/wiki/Probability_theory):**\n",
    "  As Laplace said, the probability of an event with respect to a sample space is the \"number of favorable cases\" (outcomes from the sample space that are in the event) divided by the \"number of all the cases\" in the sample space (assuming \"nothing leads us to expect that any one of these cases should occur more than any other\"). Since this is a proper fraction, probability will always be a number between 0 (representing an impossible event) and 1 (representing a certain event).\n",
    "<br>*For example, the probability of an even die roll is 3/6 = 1/2.*\n",
    "\n",
    "This notebook will explore these concepts in a concrete way using Python code. The code is meant to be succint and explicit, and fast enough to handle sample spaces with millions of outcomes. If you need to handle trillions, you'll want a more efficient implementation. I also have  [another notebook](http://nbviewer.jupyter.org/url/norvig.com/ipython/ProbabilityParadox.ipynb) that covers  paradoxes in Probability Theory. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# `P` is for Probability\n",
    "\n",
    "The code below implements Laplace's quote directly: *Probability is thus simply a fraction whose numerator is the number of favorable cases and whose denominator is the number of all the cases possible.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "from fractions import Fraction\n",
    "\n",
    "def P(event, space): \n",
    "    \"The probability of an event, given a sample space.\"\n",
    "    return Fraction(cases(favorable(event, space)), \n",
    "                    cases(space))\n",
    "\n",
    "favorable = set.intersection # Outcomes that are in the event and in the sample space\n",
    "cases     = len              # The number of cases is the length, or size, of a set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    " \n",
    "# Warm-up Problem: Die Roll"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "What's the probability of rolling an even number with a single six-sided fair die? Mathematicians traditionally use a single capital letter to denote a sample space; I'll use `D` for the die:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "D     = {1, 2, 3, 4, 5, 6} # a sample space\n",
    "even  = {   2,    4,    6} # an event\n",
    "\n",
    "P(even, D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Good to confirm what we already knew. We can explore some other events:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prime = {2, 3, 5, 7, 11, 13}\n",
    "odd   = {1, 3, 5, 7, 9, 11, 13}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P(odd, D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P((even | prime), D) # The probability of an even or prime die roll"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P((odd & prime), D) # The probability of an odd prime die roll"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# Card Problems\n",
    "\n",
    "Consider dealing a hand of five playing cards. An individual card has a rank and  suit, like `'J♥'` for the Jack of Hearts, and a `deck` has 52 cards:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "suits = u'♥♠♦♣'\n",
    "ranks = u'AKQJT98765432'\n",
    "deck  = [r + s for r in ranks for s in suits]\n",
    "len(deck)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I want to define `Hands` as the sample space of all 5-card combinations from `deck`. The function `itertools.combinations` does most of the work; we than concatenate each combination into a space-separated string:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "\n",
    "def combos(items, n):\n",
    "    \"All combinations of n items; each combo as a space-separated str.\"\n",
    "    return set(map(' '.join, itertools.combinations(items, n)))\n",
    "\n",
    "Hands = combos(deck, 5)\n",
    "len(Hands)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are too many hands to look at them all, but we can sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "random.sample(Hands, 7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.sample(deck, 7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Now we can answer questions like the probability of being dealt a flush (5 cards of the same suit):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "flush = {hand for hand in Hands if any(hand.count(suit) == 5 for suit in suits)}\n",
    "\n",
    "P(flush, Hands)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Or the probability of four of a kind:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "four_kind = {hand for hand in Hands if any(hand.count(rank) == 4 for rank in ranks)}\n",
    "\n",
    "P(four_kind, Hands)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "\n",
    "\n",
    "# Urn Problems\n",
    "\n",
    "Around 1700, Jacob Bernoulli wrote about removing colored balls from an urn in his landmark treatise *[Ars Conjectandi](https://en.wikipedia.org/wiki/Ars_Conjectandi)*, and ever since then, explanations of probability have relied on [urn problems](https://www.google.com/search?q=probability+ball+urn). (You'd think the urns would be empty by now.) \n",
    "\n",
    "![Jacob Bernoulli](http://www2.stetson.edu/~efriedma/periodictable/jpg/Bernoulli-Jacob.jpg)\n",
    "<center><a href=\"https://en.wikipedia.org/wiki/Jacob_Bernoulli\">Jacob Bernoulli</a><br>1700</center>\n",
    "\n",
    "For example, here is a three-part problem [adapted](http://mathforum.org/library/drmath/view/69151.html)  from mathforum.org:\n",
    "\n",
    "> *An urn contains 6 blue, 9 red, and 8 white balls.  We select six balls at random. What is the probability of each of these  outcomes:*\n",
    "\n",
    "> - *All balls are red*.\n",
    "- *3 are blue, and 1 is red, and 2 are white, *.\n",
    "- *Exactly 4 balls are white*.\n",
    "\n",
    "We'll start by defining the contents of the urn. A `set` can't contain multiple objects that are equal to each other, so I'll call the blue balls `'B1'` through `'B6'`, rather than trying to have 6 balls all called `'B'`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def balls(color, n):\n",
    "    \"A set of n numbered balls of the given color.\"\n",
    "    return {color + str(i)\n",
    "            for i in range(1, n + 1)}\n",
    "\n",
    "urn = balls('B', 6) | balls('R', 9) | balls('W', 8)\n",
    "urn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Now we can define the sample space, `U6`, as the set of all 6-ball combinations:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "U6 = combos(urn, 6)\n",
    "\n",
    "random.sample(U6, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define  `select` such that `select('R', 6)` is the event of picking 6 red balls from the urn:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def select(color, n, space=U6):\n",
    "    \"The subset of the sample space with exactly `n` balls of given `color`.\"\n",
    "    return {s for s in space if s.count(color) == n}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I can answer the three questions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P(select('R', 6), U6) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P(select('B', 3)  & select('R', 1) & select('W', 2), U6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P(select('W', 4), U6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## Urn problems via arithmetic\n",
    "\n",
    "Let's verify these calculations using basic arithmetic, rather than exhaustive counting. First, how many ways can I choose 6 out of 9 red balls? It could be any of the 9 for the first ball, any of 8 remaining for the second, and  so on down to any of the remaining 4 for the sixth and final ball. But we don't care about the *order* of the six balls, so divide that product by the number of permutations of 6 things, which is 6!, giving us \n",
    "9 &times; 8 &times; 7 &times; 6 &times; 5 &times; 4 / 6! = 84. In general, the number of ways of choosing *c* out of *n* items is (*n* choose *c*) = *n*! / ((*n* - *c*)! &times; c!).\n",
    "We can translate that to code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "from math import factorial\n",
    "\n",
    "def choose(n, c):\n",
    "    \"Number of ways to choose c items from a list of n items.\"\n",
    "    return factorial(n) // (factorial(n - c) * factorial(c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "choose(9, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can verify the answers to the three problems. (Since `P` computes a ratio and `choose` computes a count,\n",
    "I multiply the left-hand-side by `N`, the length of the sample space, to make both sides be counts.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = len(U6)\n",
    "\n",
    "N * P(select('R', 6), U6) == choose(9, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N * P(select('B', 3) & select('W', 2) & select('R', 1), U6) == choose(6, 3) * choose(8, 2) * choose(9, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N * P(select('W', 4), U6) == choose(8, 4) * choose(6 + 9, 2)  # (6 + 9 non-white balls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can solve all these problems just by counting; all you ever needed to know about probability problems you learned from Sesame Street:\n",
    "\n",
    "![The Count](http://img2.oncoloring.com/count-dracula-number-thir_518b77b54ba6c-p.gif)\n",
    "<center><a href=\"https://en.wikipedia.org/wiki/Count_von_Count\">The Count</a><br>1972&mdash;</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# Non-Equiprobable Outcomes\n",
    "\n",
    "So far, we have accepted Laplace's assumption that *nothing leads us to expect that any one of these cases should occur more than any other*.\n",
    "In real life, we often get outcomes that are not equiprobable--for example, a loaded die favors one side over the others.  We will introduce three more vocabulary items:\n",
    "\n",
    "* [Frequency](https://en.wikipedia.org/wiki/Frequency_%28statistics%29): a non-negative number describing how often an outcome occurs. Can be a count like 5, or a ratio like 1/6.\n",
    "\n",
    "* [Distribution](http://mathworld.wolfram.com/StatisticalDistribution.html): A mapping from outcome to frequency of that outcome. We will allow sample spaces to be distributions. \n",
    "\n",
    "* [Probability Distribution](https://en.wikipedia.org/wiki/Probability_distribution): A probability distribution\n",
    "is a distribution whose frequencies sum to 1. \n",
    "\n",
    "\n",
    "I could implement distributions with `Dist = dict`, but instead I'll make `Dist` a subclass `collections.Counter`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "        \n",
    "class Dist(Counter): \n",
    "    \"A Distribution of {outcome: frequency} pairs.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because a `Dist` is a `Counter`, we can initialize it in any of the following ways:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A set of equiprobable outcomes:\n",
    "Dist({1, 2, 3, 4, 5, 6})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A collection of outcomes, with repetition indicating frequency:\n",
    "Dist('THHHTTHHT')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A mapping of {outcome: frequency} pairs:\n",
    "Dist({'H': 5, 'T': 4})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Keyword arguments:\n",
    "Dist(H=5, T=4) == Dist({'H': 5}, T=4) == Dist('TTTT', H=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I will modify the code to handle distributions.\n",
    "Here's my plan:\n",
    "\n",
    "- Sample spaces and events can both be specified as either a `set` or a `Dist`.\n",
    "- The sample space can be a non-probability distribution like `Dist(H=50, T=50)`; the results\n",
    "will be the same as if the sample space had been a true probability distribution like `Dist(H=1/2, T=1/2)`.\n",
    "- The function `cases` now sums the frequencies in a distribution (it previously counted the length).\n",
    "- The function `favorable` now returns a `Dist` of favorable outcomes and their frequencies (not a `set`).\n",
    "- I will redefine `Fraction` to use `\"/\"`, not `fractions.Fraction`, because frequencies might be floats.\n",
    "- `P` is unchanged.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cases(outcomes): \n",
    "    \"The total frequency of all the outcomes.\"\n",
    "    return sum(Dist(outcomes).values())\n",
    "\n",
    "def favorable(event, space):\n",
    "    \"A distribution of outcomes from the sample space that are in the event.\"\n",
    "    space = Dist(space)\n",
    "    return Dist({x: space[x] \n",
    "                 for x in space if x in event})\n",
    "\n",
    "def Fraction(n, d): return n / d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, here's the probability of rolling an even number with a crooked die that is loaded to prefer 6:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Crooked = Dist({1: 0.1, 2: 0.1, 3: 0.1, 4: 0.1, 5: 0.1, 6: 0.5})\n",
    "\n",
    "P(even, Crooked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As another example, an [article](http://people.kzoo.edu/barth/math105/moreboys.pdf) gives the following counts for two-child families in Denmark, where `GB` means a family where the first child is a girl and the second a boy (I'm aware that not all births can be classified as the binary \"boy\" or \"girl,\" but the data was reported that way):\n",
    "\n",
    "    GG: 121801    GB: 126840\n",
    "    BG: 127123    BB: 135138"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "DK = Dist(GG=121801, GB=126840,\n",
    "          BG=127123, BB=135138)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_girl  = {'GG', 'GB'}\n",
    "P(first_girl, DK)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "second_girl = {'GG', 'BG'}\n",
    "P(second_girl, DK)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "This says that the probability of a girl is somewhere between 48% and 49%. The probability of a girl is very slightly higher for the second child. \n",
    "\n",
    "Given the first child, are you more likely to have a second child of the same sex?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "same = {'GG', 'BB'}\n",
    "P(same, DK)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Yes, but only by about 0.3%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# Predicates as events\n",
    "\n",
    "To calculate the probability of an even die roll, I originally said\n",
    "\n",
    "    even = {2, 4, 6}\n",
    "    \n",
    "But that's inelegant&mdash;I had to explicitly enumerate all the even numbers from one to six. If I ever wanted to deal with a twelve or twenty-sided die, I would have to go back and redefine `even`.  I would prefer to define `even` once and for all like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def even(n): return n % 2 == 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Now in order to make `P(even, D)` work, I'll allow an `Event` to be either a collection of outcomes or a `callable` predicate (that is, a function that returns true for outcomes that are part of the event). I don't need to modify `P`, but  `favorable` will have  to convert a callable `event` to a `set`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def favorable(event, space):\n",
    "    \"A distribution of outcomes from the sample space that are in the event.\"\n",
    "    if callable(event):\n",
    "        event = {x for x in space if event(x)}\n",
    "    space = Dist(space)\n",
    "    return Dist({x: space[x] \n",
    "                 for x in space if x in event})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "favorable(even, D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P(even, D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll define `die` to make a sample space for an *n*-sided die:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def die(n): return set(range(1, n + 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "favorable(even, die(12))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P(even, die(12))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P(even, die(2000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P(even, die(2001))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We can define more interesting events using predicates; for example we can determine the probability that the sum of  rolling *d* 6-sided dice is prime:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def sum_dice(d): return Dist(sum(dice) for dice in itertools.product(D, repeat=d))\n",
    "\n",
    "def is_prime(n): return (n > 1 and not any(n % i == 0 for i in range(2, n)))\n",
    "\n",
    "for d in range(1, 9):\n",
    "    p = P(is_prime, sum_dice(d))\n",
    "    print(\"P(is_prime, sum_dice({})) = {}\".format(d, round(p, 3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fermat and Pascal: The Unfinished Game\n",
    "\n",
    "<table>\n",
    "<tr><td><img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/9/98/Pierre_de_Fermat2.png/140px-Pierre_de_Fermat2.png\"><center><a href=\"https://en.wikipedia.org/wiki/Pierre_de_Fermat\">Pierre de Fermat</a><br>1654\n",
    "<td><img src=\"https://www.umass.edu/wsp/images/pascal.jpg\"><center><a href=\"https://en.wikipedia.org/wiki/Blaise_Pascal\">Blaise Pascal]</a><br>1654\n",
    "</table>\n",
    "\n",
    "Consider a gambling game consisting of tossing a coin repeatedly. Player H wins the game as soon as a total of 10 heads come up, and T wins if a total of 10 tails come up before H wins. If the game is interrupted when H has 8 heads and T has 7 tails, how should the pot of money (which happens to be 100 Francs) be split?  Here are some proposals, and arguments against them:\n",
    "- It is uncertain, so just split the pot 50-50. \n",
    "<br>*No, because surely H is more likely to win.*\n",
    "- In proportion to each player's current score, so H gets a 8/(8+7) share. \n",
    "<br>*No, because if the score was 0 heads to 1 tail, H should get more than 0/1.*\n",
    "- In proportion to how many tosses the opponent needs to win, so H gets 3/(3+2). \n",
    "<br>*This seems better, but no, if H is 9 away and T is only 1 away from winning, then it seems that giving H a 1/10 share is too much.*\n",
    "\n",
    "In 1654, Blaise Pascal and Pierre de Fermat corresponded on this problem, with Fermat [writing](http://mathforum.org/isaac/problems/prob1.html):\n",
    "\n",
    ">Dearest Blaise,\n",
    "\n",
    ">As to the problem of how to divide the 100 Francs, I think I have found a solution that you will find to be fair. Seeing as I needed only two points to win the game, and you needed 3, I think we can establish that after four more tosses of the coin, the game would have been over. For, in those four tosses, if you did not get the necessary 3 points for your victory, this would imply that I had in fact gained the necessary 2 points for my victory. In a similar manner, if I had not achieved the necessary 2 points for my victory, this would imply that you had in fact achieved at least 3 points and had therefore won the game. Thus, I believe the following list of possible endings to the game is exhaustive. I have denoted 'heads' by an 'h', and tails by a 't.' I have starred the outcomes that indicate a win for myself.\n",
    "\n",
    ">       h h h h *       h h h t *       h h t h *       h h t t *\n",
    ">       h t h h *       h t h t *       h t t h *       h t t t\n",
    ">       t h h h *       t h h t *       t h t h *       t h t t\n",
    ">       t t h h *       t t h t         t t t h         t t t t\n",
    "\n",
    ">I think you will agree that all of these outcomes are equally likely. Thus I believe that we should divide the stakes by the ration 11:5 in my favor, that is, I should receive (11/16)&times;100 = 68.75 Francs, while you should receive 31.25 Francs.\n",
    "\n",
    "\n",
    ">I hope all is well in Paris,\n",
    "\n",
    ">Your friend and colleague,\n",
    "\n",
    ">Pierre\n",
    "\n",
    "Pascal agreed with this solution, and [replied](http://mathforum.org/isaac/problems/prob2.html) with a generalization that made use of his previous invention, Pascal's Triangle. There's even [a book](https://smile.amazon.com/Unfinished-Game-Pascal-Fermat-Seventeenth-Century/dp/0465018963?sa-no-redirect=1) about it.\n",
    "\n",
    "We can solve the problem with the tools we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def win_unfinished_game(h, t):\n",
    "    \"The probability that H will win the unfinished game, given the number of points needed by H and T to win.\"\n",
    "    return P(at_least(h, 'h'), finishes(h, t))\n",
    "\n",
    "def at_least(n, item):\n",
    "    \"The event of getting at least n instances of item in an outcome.\"\n",
    "    return lambda outcome: outcome.count(item) >= n\n",
    "    \n",
    "def finishes(h, t):\n",
    "    \"All finishes of a game where player H needs h points to win and T needs t.\"\n",
    "    tosses = ['ht'] * (h + t - 1)\n",
    "    return set(itertools.product(*tosses))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can generate the 16 equiprobable finished that Pierre wrote about:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "finishes(2, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can find the 11 of them that are favorable to player `H`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "favorable(at_least(2, 'h'), finishes(2, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can answer the question:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "100 * win_unfinished_game(2, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We agree with Pascal and Fermat; we're in good company!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Newton's Answer to a Problem by Pepys\n",
    "\n",
    "<table>\n",
    "<tr><td><img src=\"http://scienceworld.wolfram.com/biography/pics/Newton.jpg\"><center><a href=\"https://en.wikipedia.org/wiki/Isaac_Newton\">Isaac Newton</a><br>1693</center>\n",
    "<td><img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/f/f8/Samuel_Pepys_portrait.jpg/148px-Samuel_Pepys_portrait.jpg\"><center><a href=\"https://en.wikipedia.org/wiki/Samuel_Pepys\">Samuel Pepys</a><br>1693</center>\n",
    "</table>\n",
    "\n",
    "Let's jump ahead from 1654 all the way to 1693, [when](http://fermatslibrary.com/s/isaac-newton-as-a-probabilist) Samuel Pepys wrote to Isaac Newton posing the problem:\n",
    "\n",
    "> Which of the following three propositions has the greatest chance of success? \n",
    "  1. Six fair dice are tossed independently and at least one “6” appears. \n",
    "  2. Twelve fair dice are tossed independently and at least two “6”s appear. \n",
    "  3. Eighteen fair dice are tossed independently and at least three “6”s appear.\n",
    "  \n",
    "Newton was able to answer the question correctly (although his reasoning was not quite right); let's see how we can do. Since we're only interested in whether a die comes up as \"6\" or not, we can define a single die like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "die6 = Dist({6: 1/6, '-': 5/6})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can define the joint distribution formed by combining two independent distribution like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def joint(A, B, combine='{}{}'.format):\n",
    "    \"\"\"The joint distribution of two independent distributions. \n",
    "    Result is all entries of the form {'ab': frequency(a) * frequency(b)}\"\"\"\n",
    "    return Dist({combine(a, b): A[a] * B[b]\n",
    "                 for a in A for b in B})\n",
    "\n",
    "joint(die6, die6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the joint distribution from rolling *n* dice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dice(n, die):\n",
    "    \"Joint probability distribution from rolling `n` dice.\"\n",
    "    if n == 1:\n",
    "        return die\n",
    "    else:\n",
    "        return joint(die, dice(n - 1, die))\n",
    "    \n",
    "dice(4, die6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to determine which proposition is more likely to have the required number of sixes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P(at_least(1, '6'), dice(6, die6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P(at_least(2, '6'), dice(12, die6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P(at_least(3, '6'), dice(18, die6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We reach the same conclusion Newton did, that the best chance is rolling six dice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# More Urn Problems: M&Ms and Bayes\n",
    "\n",
    "Here's another urn problem (actually a \"bag\" problem) [from](http://allendowney.blogspot.com/2011/10/my-favorite-bayess-theorem-problems.html) prolific Python/Probability pundit [Allen Downey ](http://allendowney.blogspot.com/):\n",
    "\n",
    "> The blue M&M was introduced in 1995.  Before then, the color mix in a bag of plain M&Ms was (30% Brown, 20% Yellow, 20% Red, 10% Green, 10% Orange, 10% Tan).  Afterward it was (24% Blue , 20% Green, 16% Orange, 14% Yellow, 13% Red, 13% Brown). \n",
    "A friend of mine has two bags of M&Ms, and he tells me that one is from 1994 and one from 1996.  He won't tell me which is which, but he gives me one M&M from each bag.  One is yellow and one is green.  What is the probability that the yellow M&M came from the 1994 bag?\n",
    "\n",
    "To solve this problem, we'll first create distributions for each bag: `bag94` and `bag96`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "bag94 = Dist(brown=30, yellow=20, red=20, green=10, orange=10, tan=10)\n",
    "bag96 = Dist(blue=24, green=20, orange=16, yellow=14, red=13, brown=13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Next, define `MM` as the joint distribution&mdash;the sample space for picking one M&M from each bag. The outcome `'94:yellow 96:green'` means that a yellow M&M was selected from the 1994 bag and a green one from the 1996 bag. In this problem we don't get to see the actual outcome; we just see some evidence about the outcome, that it contains a yellow and a green."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "MM = joint(bag94, bag96, '94:{} 96:{}'.format)\n",
    "MM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We observe that \"One is yellow and one is green\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def yellow_and_green(outcome): return 'yellow' in outcome and 'green' in outcome\n",
    "\n",
    "favorable(yellow_and_green, MM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given this observation, we want to know \"What is the probability that the yellow M&M came from the 1994 bag?\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def yellow94(outcome): return '94:yellow' in outcome\n",
    "\n",
    "P(yellow94, favorable(yellow_and_green, MM))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "So there is a 74% chance that the yellow comes from the 1994 bag.\n",
    "\n",
    "Answering this question was straightforward: just like all the other probability problems, we simply create a sample space, and use `P` to pick out the probability of the event in question, given what we know about the outcome.\n",
    "But in a sense it is curious that we were able to solve this problem with the same methodology as the others: this problem comes from a section titled **My favorite Bayes's Theorem Problems**, so one would expect that we'd need to invoke Bayes Theorem to solve it.  The computation above shows that that is not necessary. \n",
    "\n",
    "![Bayes](http://img1.ph.126.net/xKZAzeOv_mI8a4Lwq7PHmw==/2547911489202312541.jpg)\n",
    "<center><a href=\"https://en.wikipedia.org/wiki/Thomas_Bayes\">Rev. Thomas Bayes</a><br>1701-1761\n",
    "</center>\n",
    "\n",
    "Of course, we *could* solve it using Bayes Theorem. Why is Bayes Theorem recommended? Because we are asked about the probability of an outcome given the evidence&mdash;the probability the yellow came from the 94 bag, given that there is a yellow and a green. But the problem statement doesn't directly tell us the probability of that outcome given the evidence; it just tells us the probability of the evidence given the outcome. \n",
    "\n",
    "Before we see the colors of the M&Ms, there are two hypotheses, `A` and `B`, both with equal probability:\n",
    "\n",
    "    A: first M&M from 94 bag, second from 96 bag\n",
    "    B: first M&M from 96 bag, second from 94 bag\n",
    "    P(A) = P(B) = 0.5\n",
    "    \n",
    "Then we get some evidence:\n",
    "    \n",
    "    E: first M&M yellow, second green\n",
    "    \n",
    "We want to know the probability of hypothesis `A`, given the evidence:\n",
    "    \n",
    "    P(A | E)\n",
    "    \n",
    "That's not easy to calculate (except by enumerating the sample space, which our `P` function does). But Bayes Theorem says:\n",
    "    \n",
    "    P(A | E) = P(E | A) * P(A) / P(E)\n",
    "    \n",
    "The quantities on the right-hand-side are easier to calculate:\n",
    "    \n",
    "    P(E | A) = 0.20 * 0.20 = 0.04\n",
    "    P(E | B) = 0.10 * 0.14 = 0.014\n",
    "    P(A)     = 0.5\n",
    "    P(B)     = 0.5\n",
    "    P(E)     = P(E | A) * P(A) + P(E | B) * P(B) \n",
    "             = 0.04     * 0.5  + 0.014    * 0.5   =   0.027\n",
    "    \n",
    "And we can get a final answer:\n",
    "    \n",
    "    P(A | E) = P(E | A) * P(A) / P(E) \n",
    "             = 0.04     * 0.5  / 0.027 \n",
    "             = 0.7407407407\n",
    "             \n",
    "You have a choice: Bayes Theorem allows you to do less calculation at the cost of more algebra; that is a great trade-off if you are working with pencil and paper. Enumerating the sample space allows you to do less algebra at the cost of more calculation; usually a good trade-off if you have a computer. But regardless of the approach you use, it is important to understand Bayes theorem and how it works.\n",
    "\n",
    "There is one important question that Allen Downey does not address: *would you  eat twenty-year-old M&Ms*?\n",
    "&#128552;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "<hr>\n",
    "\n",
    "# Simulation\n",
    "\n",
    "Sometimes it is inconvenient, difficult, or even impossible to explicitly enumerate a sample space. Perhaps the sample space is infinite, or perhaps it is just very large and complicated (perhaps with a bunch of low-probability outcomes that don't seem very important). In that case, we might feel more confident in writing a program to *simulate* a random outcome. *Random sampling* from such a simulation\n",
    "can give an accurate estimate of probability.\n",
    "\n",
    "# Simulating Monopoly\n",
    "\n",
    "![Mr. Monopoly](http://buckwolf.org/a.abcnews.com/images/Entertainment/ho_hop_go_050111_t.jpg)<center>[Mr. Monopoly](https://en.wikipedia.org/wiki/Rich_Uncle_Pennybags)<br>1940&mdash;\n",
    "\n",
    "Consider [problem 84](https://projecteuler.net/problem=84) from the excellent [Project Euler](https://projecteuler.net), which asks for the probability that a player in the game Monopoly ends a roll on each of the squares on the board.  To answer this we need to take into account die rolls, chance and community chest cards, and going to jail (from the \"go to jail\" space, from a card, or from rolling doubles three times in a row). We do not need to take into account anything about acquiring properties or exchanging money or winning or losing the game, because these events don't change a player's location. \n",
    "\n",
    "A game of Monopoly can go on forever, so the sample space is infinite. Even if we limit the sample space to say, 1000 rolls, there are $21^{1000}$ such sequences of rolls, and even more possibilities when we consider drawing cards. So it is infeasible to explicitly represent the sample space. There are techniques for representing the problem as\n",
    "a Markov decision problem (MDP) and solving it, but the math is complex (a [paper](https://faculty.math.illinois.edu/~bishop/monopoly.pdf) on the subject runs 15 pages).\n",
    "\n",
    "The simplest approach is to implement a simulation and run it for, say, a million rolls.  Here is the code for a simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import deque as Deck # a Deck of community chest or chance cards\n",
    "\n",
    "# The Monopoly board, as specified by https://projecteuler.net/problem=84\n",
    "(GO,   A1, CC1, A2,  T1, R1, B1,  CH1, B2, B3,\n",
    " JAIL, C1, U1,  C2,  C3, R2, D1,  CC2, D2, D3, \n",
    " FP,   E1, CH2, E2,  E3, R3, F1,  F2,  U2, F3, \n",
    " G2J,  G1, G2,  CC3, G3, R4, CH3, H1,  T2, H2) = board = range(40)\n",
    "\n",
    "# A card is either a square, a set of squares meaning advance to the nearest, \n",
    "# a -3 to go back 3 spaces, or None meaning no change to location.\n",
    "CC_deck = Deck([GO, JAIL] +  14 * [None])\n",
    "CH_deck = Deck([GO, JAIL, C1, E3, H2, R1, -3, {U1, U2}] \n",
    "               + 2 * [{R1, R2, R3, R4}] + 6 * [None])\n",
    "\n",
    "def monopoly(rolls):\n",
    "    \"\"\"Simulate given number of dice rolls of a Monopoly game, \n",
    "    and return the counts of how often each square is visited.\"\"\"\n",
    "    counts = [0] * len(board)\n",
    "    doubles = 0 # Number of consecutive doubles rolled\n",
    "    random.shuffle(CC_deck)\n",
    "    random.shuffle(CH_deck)\n",
    "    goto(GO)\n",
    "    for _ in range(rolls):\n",
    "        d1, d2 = random.randint(1, 6), random.randint(1, 6)\n",
    "        doubles = (doubles + 1 if d1 == d2 else 0)\n",
    "        goto(here + d1 + d2)\n",
    "        if here == G2J or doubles == 3:\n",
    "            goto(JAIL)\n",
    "            doubles = 0\n",
    "        elif here in (CC1, CC2, CC3):\n",
    "            do_card(CC_deck)\n",
    "        elif here in (CH1, CH2, CH3):\n",
    "            do_card(CH_deck)\n",
    "        counts[here] += 1\n",
    "    return counts\n",
    "\n",
    "def goto(square):\n",
    "    \"Update 'here' to be this square (and handle passing GO).\"\n",
    "    global here\n",
    "    here = square % len(board)        \n",
    "\n",
    "def do_card(deck):\n",
    "    \"Take the top card from deck and do what it says.\"\n",
    "    card = deck.popleft()       # The top card\n",
    "    deck.append(card)           # Move top card to bottom of deck\n",
    "    if card == None:            # Don't move\n",
    "        pass\n",
    "    elif card == -3:            # Go back 3 spaces\n",
    "        goto(here - 3)\n",
    "    elif isinstance(card, set): # Advance to next railroad or utility\n",
    "        next1 = min({place for place in card if place > here} or card)\n",
    "        goto(next1)\n",
    "    else:                       # Go to destination named on card\n",
    "        goto(card)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let's run the simulation for a million dice rolls:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "counts = monopoly(10**6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "And print a table of square names and their percentages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "property_names = \"\"\"\n",
    " GO,   A1, CC1, A2,  T1, R1, B1,  CH1, B2, B3,\n",
    " JAIL, C1, U1,  C2,  C3, R2, D1,  CC2, D2, D3, \n",
    " FP,   E1, CH2, E2,  E3, R3, F1,  F2,  U2, F3, \n",
    " G2J,  G1, G2,  CC3, G3, R4, CH3, H1,  T2, H2\"\"\".replace(',', ' ').split()\n",
    "\n",
    "for (c, n) in sorted(zip(counts, property_names), reverse=True):\n",
    "    print('{:4} {:.2%}'.format(n, c / sum(counts)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "There is one square far above average: `JAIL`, at a little over 6%. There are four squares far below average: the three chance squares, `CH1`, `CH2`, and `CH3`, at around 1% (because 10 of the 16 chance cards send the player away from the square), and the \"Go to Jail\" square, which has a frequency of 0 because you can't end a turn there. The other squares are around 2% to 3% each, which you would expect, because 100% / 40 = 2.5%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Central Limit Theorem \n",
    "\n",
    "We have covered the concept of *distributions* of outcomes. You may have heard of the *normal distribution*, the *bell-shaped curve.* In Python it is called `random.normalvariate` (also `random.gauss`). We can plot it with the help of the `repeated_hist` function defined below, which samples a distribution `n` times and displays a histogram of the results. (*Note:* in this section I am using \"distribution\" to mean a function that, each time it is called, returns a random sample from a distribution. I am not using it to mean a mapping of type `Dist`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline \n",
    "import matplotlib.pyplot as plt\n",
    "from statistics import mean\n",
    "from random import normalvariate, triangular, choice, vonmisesvariate, uniform\n",
    "\n",
    "def normal(mu=0, sigma=1): return random.normalvariate(mu, sigma)\n",
    "\n",
    "def repeated_hist(dist, n=10**6, bins=100):\n",
    "    \"Sample the distribution n times and make a histogram of the results.\"\n",
    "    samples = [dist() for _ in range(n)]\n",
    "    plt.hist(samples, bins=bins, density=True)\n",
    "    plt.title('{} (μ = {:.1f})'.format(dist.__name__, mean(samples)))\n",
    "    plt.grid(axis='x')\n",
    "    plt.yticks([], '')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normal distribution\n",
    "repeated_hist(normal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is this distribution called *normal*? The **Central Limit Theorem** says  that it is the ultimate limit of other distributions, as follows (informally):\n",
    "- Gather *k* independent distributions. They need not be normal-shaped.\n",
    "- Define a new distribution to be the result of sampling one number from each of the *k* independent distributions and adding them up.\n",
    "- As long as *k* is not too small, and the component distributions are not super-pathological, then the new distribution will tend towards a normal distribution.\n",
    "\n",
    "Here's a simple example: summing ten independent die rolls:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sum10dice(): return sum(random.randint(1, 6) for _ in range(10))\n",
    "\n",
    "repeated_hist(sum10dice, bins=range(10, 61))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As another example, let's take just *k* = 5 component distributions representing the per-game scores of 5 basketball players, and then sum them together to form the new distribution, the team score. I'll be creative in defining the distributions for each player, but [historically accurate](https://www.basketball-reference.com/teams/GSW/2016.html) in the mean for each distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SC(): return max(0, normal(12.1, 3) + 3 * triangular(1, 13, 4))                     # 30.1\n",
    "def KT(): return max(0, triangular(8, 22, 15.3) + choice((0, 3 * triangular(1, 9, 4)))) # 22.1\n",
    "def DG(): return max(0, vonmisesvariate(30, 2) * 3.08)                                  # 14.0\n",
    "def HB(): return max(0, choice((normal(6.7, 1.5), normal(16.7, 2.5))))                  # 11.7\n",
    "def BE(): return max(0, normal(17, 3) + uniform(0, 40))                                 # 37.0\n",
    "\n",
    "team = (SC, KT, DG, HB, BE)\n",
    "\n",
    "def Team(team=team): return sum(player() for player in team)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for player in team: \n",
    "    repeated_hist(player, bins=range(70))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that none of the players have a distribution that looks like a normal distribution: `SC` is skewed to one side (the mean is 5 points to the right of the peak); the three next players have bimodal distributions; and `BE` is too flat on top. \n",
    "\n",
    "Now we define the team score to be the sum of the *k* = 5 players, and display this new distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "repeated_hist(Team, bins=range(50, 180))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sure enough, this looks very much like a normal distribution. The **Central Limit Theorem** appears to hold in this case. But I have to say: \"Central Limit\" is not a very evocative name, so I propose we re-name this as the **Strength in Numbers Theorem**, to indicate the fact that if you have a lot of numbers, you tend to get the expected result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# Conclusion\n",
    "\n",
    "We've had an interesting tour and met some giants of the field: Laplace, Bernoulli, Fermat, Pascal, Bayes, Newton, ... even Mr. Monopoly and The Count.\n",
    "\n",
    "The conclusion is: be methodical in defining the sample space and the event(s) of interest, and be careful in counting the number of outcomes in the numerator and denominator. and you can't go wrong. Easy as 1-2-3. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "# Appendix: Continuous Sample Spaces\n",
    "\n",
    "Everything up to here has been about discrete, finite sample spaces, where we can *enumerate* all the possible outcomes. \n",
    "\n",
    "But a reader asked about *continuous* sample spaces, such as the space of real numbers. The principles are the same: probability is still the ratio of the favorable cases to all the cases, but now instead of *counting* cases, we have to (in general) compute integrals to compare the sizes of cases. \n",
    "Here we will cover a simple example, which we first solve approximately by simulation, and then exactly by calculation.\n",
    "\n",
    "## The Hot New Game Show Problem: Simulation\n",
    "\n",
    "Oliver Roeder posed [this problem](http://fivethirtyeight.com/features/can-you-win-this-hot-new-game-show/) in the 538 *Riddler* blog:\n",
    "\n",
    ">Two players go on a hot new game show called *Higher Number Wins.* The two go into separate booths, and each presses a button, and a random number between zero and one appears on a screen. (At this point, neither knows the other’s number, but they do know the numbers are chosen from a standard uniform distribution.) They can choose to keep that first number, or to press the button again to discard the first number and get a second random number, which they must keep. Then, they come out of their booths and see the final number for each player on the wall. The lavish grand prize — a case full of gold bullion — is awarded to the player who kept the higher number. Which number is the optimal cutoff for players to discard their first number and choose another? Put another way, within which range should they choose to keep the first number, and within which range should they reject it and try their luck with a second number?\n",
    "\n",
    "We'll use this notation:\n",
    "- **A**, **B**: the two players.\n",
    "- *A*, *B*: the cutoff values they choose: the lower bound of the range of first numbers they will accept.\n",
    "- *a*, *b*: the actual random numbers that appear on the screen.\n",
    "\n",
    "For example, if player **A** chooses a cutoff of *A* = 0.6, that means that **A** would accept any first number greater than 0.6, and reject any number below that cutoff. The question is: What cutoff, *A*, should player **A** choose to maximize the chance of winning, that is, maximize P(*a* > *b*)?\n",
    "\n",
    "First, simulate the number that a player with a given cutoff gets (note that `random.random()` returns a float sampled uniformly from the interval [0..1]):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number= random.random\n",
    "\n",
    "def strategy(cutoff):\n",
    "    \"Play the game with given cutoff, returning the first or second random number.\"\n",
    "    first = number()\n",
    "    return first if first > cutoff else number()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "strategy(.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compare the numbers returned with a cutoff of *A* versus a cutoff of *B*, and repeat for a large number of trials; this gives us an estimate of the probability that cutoff *A* is better than cutoff *B*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Pwin(A, B, trials=20000):\n",
    "    \"The probability that cutoff A wins against cutoff B.\"\n",
    "    return mean(strategy(A) > strategy(B) \n",
    "                for _ in range(trials))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pwin(0.6, 0.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define a function, `top`, that considers a collection of possible cutoffs, estimate the probability for each cutoff playing against each other cutoff, and returns a list with the `N` top cutoffs (the ones that defeated the most number of opponent cutoffs), and the number of opponents they defeat: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def top(N, cutoffs):\n",
    "    \"Return the N best cutoffs and the number of opponent cutoffs they beat.\"\n",
    "    winners = Counter(A if Pwin(A, B) > 0.5 else B\n",
    "                      for (A, B) in itertools.combinations(cutoffs, 2))\n",
    "    return winners.most_common(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import arange\n",
    "\n",
    "top(10, arange(0.5, 1.0, 0.01))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get a good idea of the top cutoffs, but they are close to each other, so we can't quite be sure which is best, only that the best is somewhere around 0.60. We could get a better estimate by increasing the number of trials, but that would consume more time.\n",
    "\n",
    "## The Hot New Game Show Problem: Exact Calculation\n",
    "\n",
    "More promising is the possibility of making `Pwin(A, B)` an exact calculation. But before we get to `Pwin(A, B)`, let's solve a simpler problem: assume that both players **A** and **B** have chosen a cutoff, and have each received a number above the cutoff. What is the probability that **A** gets the higher number? We'll call this `Phigher(A, B)`. We can think of this as a two-dimensional sample space of points in the (*a*, *b*) plane, where *a* ranges from the cutoff *A* to 1 and *b* ranges from the cutoff B to 1. Here is a diagram of that two-dimensional sample space, with the cutoffs *A*=0.5 and *B*=0.6:\n",
    "\n",
    "<img src=\"http://norvig.com/ipython/probability2da.jpg\" width=413>\n",
    "\n",
    "The total area of the sample space is 0.5 &times; 0.4 = 0.20, and in general it is (1 - *A*) &middot; (1 - *B*). What about the favorable cases, where **A** beats **B**? That corresponds to the shaded triangle below:\n",
    "\n",
    "<img src=\"http://norvig.com/ipython/probability2d.jpg\" width=413>\n",
    "\n",
    "The area of a triangle is 1/2 the base times the height, or in this case, 0.4<sup>2</sup> / 2 = 0.08, and in general, (1 - *B*)<sup>2</sup> / 2. So in general we have:\n",
    "\n",
    "    Phigher(A, B) = favorable / total\n",
    "    favorable     = ((1 - B) ** 2) / 2 \n",
    "    total         = (1 - A) * (1 - B)\n",
    "    Phigher(A, B) = (((1 - B) ** 2) / 2) / ((1 - A) * (1 - B))\n",
    "    Phigher(A, B) = (1 - B) / (2 * (1 - A))\n",
    "    \n",
    "And in this specific case we have:\n",
    "\n",
    "     A = 0.5; B = 0.6\n",
    "     favorable         = 0.4 ** 2 / 2 = 0.08\n",
    "     total             = 0.5 * 0.4    = 0.20\n",
    "     Phigher(0.5, 0.6) = 0.08 / 0.20  = 0.4\n",
    "\n",
    "But note that this only works when the cutoff *A* &le; *B*; when *A* > *B*, we need to reverse things. That gives us the code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Phigher(A, B):\n",
    "    \"Probability that a sample from [A..1] is higher than one from [B..1].\"\n",
    "    if A <= B:\n",
    "        return (1 - B) / (2 * (1 - A))\n",
    "    else:\n",
    "        return 1 - Phigher(B, A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Phigher(0.5, 0.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're now ready to tackle the full game. There are four cases to consider, depending on whether **A** and **B** gets a first number that is above or below their cutoff choices:\n",
    "\n",
    "| first *a* | first *b* | P(*a*, *b*) | P(A wins &vert; *a*, *b*) |  Comment |\n",
    "|:-----:|:-----:| ----------- | ------------- | ------------ |\n",
    "| *a* > *A* | *b* > *B* | (1 - *A*) &middot; (1 - *B*) | Phigher(*A*, *B*) | Both above cutoff; both keep first numbers |\n",
    "| *a* < *A* | *b* < *B* | *A* &middot; *B* | Phigher(0, 0)  | Both below cutoff, both get new numbers from [0..1] |\n",
    "| *a* > *A* | *b* < *B* | (1 - *A*) &middot;  *B* | Phigher(*A*, 0) | **A** keeps number; **B** gets new number from [0..1] |\n",
    "| *a* < *A* | *b* > *B* | *A* &middot; (1 - *B*) | Phigher(0, *B*) | **A** gets new number from [0..1]; **B** keeps number |\n",
    "\n",
    "For example, the first row of this table says that the event of both first numbers being above their respective cutoffs has probability (1 - *A*) &middot; (1 - *B*), and if this does occur, then the probability of **A** winning is Phigher(*A*, *B*).\n",
    "We're ready to replace the old simulation-based `Pwin` with a new calculation-based version:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Pwin(A, B):\n",
    "    \"With what probability does cutoff A win against cutoff B?\"\n",
    "    return ((1-A) * (1-B) * Phigher(A, B)  # both above cutoff\n",
    "            + A * B       * Phigher(0, 0)  # both below cutoff\n",
    "            + (1-A) * B   * Phigher(A, 0)  # A above, B below\n",
    "            + A * (1-B)   * Phigher(0, B)) # A below, B above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pwin(0.5, 0.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Pwin` relies on a lot of algebra. Let's define a few tests to check for obvious errors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test():\n",
    "    assert Phigher(0.5, 0.5) == Phigher(0.75, 0.75) == Phigher(0, 0) == 0.5\n",
    "    assert Pwin(0.5, 0.5) == Pwin(0.75, 0.75) == 0.5\n",
    "    assert Phigher(.6, .5) == 0.6\n",
    "    assert Phigher(.5, .6) == 0.4\n",
    "    return 'ok'\n",
    "\n",
    "test()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's repeat the calculation with our new, exact `Pwin`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top(10, arange(0.5, 1.0, 0.01))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is good to see that the simulation and the exact calculation are in rough agreement; that gives me more confidence in both of them. We see here that 0.62 defeats all the other cutoffs, and 0.61 defeats all cutoffs except 0.62. The great thing about the exact calculation code is that it runs fast, regardless of how much accuracy we want. We can zero in on the range around 0.6:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top(10, arange(0.5, 0.7, 0.001))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This says 0.618 is best, better than 0.620. We can get even more accuracy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top(10, arange(0.617, 0.619, 0.000001))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So 0.618034 is best. Does that number [look familiar](https://en.wikipedia.org/wiki/Golden_ratio)? Can we prove that it is what I think it is?\n",
    "\n",
    "To understand the strategic possibilities, it is helpful to draw a 3D plot of `Pwin(A, B)` for values of *A* and *B* between 0 and 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from mpl_toolkits.mplot3d.axes3d import Axes3D\n",
    "\n",
    "def map2(fn, A, B):\n",
    "    \"Map fn to corresponding elements of 2D arrays A and B.\"\n",
    "    return [list(map(fn, Arow, Brow))\n",
    "            for (Arow, Brow) in zip(A, B)]\n",
    "\n",
    "cutoffs = arange(0.00, 1.00, 0.02)\n",
    "A, B = np.meshgrid(cutoffs, cutoffs)\n",
    "\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "ax.set_xlabel('A')\n",
    "ax.set_ylabel('B')\n",
    "ax.set_zlabel('Pwin(A, B)')\n",
    "A= np.array(A)\n",
    "B = np.array(B)\n",
    "\n",
    "Z = np.array(map2(Pwin, A, B));\n",
    "ax.plot_surface(A, B, Z);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does this [Pringle of Probability](http://fivethirtyeight.com/features/should-you-shoot-free-throws-underhand/) show us? The highest win percentage for **A**, the peak of the surface, occurs when *A* is around 0.5 and *B* is 0 or 1. We can confirm that, finding the maximum `Pwin(A, B)` for many different cutoff values of `A` and `B`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutoffs = (set(arange(0.00,    1.00,    0.01)) | \n",
    "           set(arange(0.500,   0.700,   0.001)) | \n",
    "           set(arange(0.61803, 0.61804, 0.000001)))\n",
    "\n",
    "def Pwin_summary(A, B): return [Pwin(A, B), 'A:', A, 'B:', B]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(Pwin_summary(A, B) for A in cutoffs for B in cutoffs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So **A** could win 62.5% of the time if only **B** would chose a cutoff of 0. But, unfortunately for **A**, a rational player **B** is not going to do that. We can ask what happens if the game is changed so that player **A** has to declare a cutoff first, and then player **B** gets to respond with a cutoff, with full knowledge of **A**'s choice. In other words, what cutoff should **A** choose to maximize `Pwin(A, B)`, given that **B** is going to take that knowledge and pick a cutoff that minimizes `Pwin(A, B)`? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(min(Pwin_summary(A, B) for B in cutoffs)\n",
    "    for A in cutoffs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And what if we run it the other way around, where **B** chooses a cutoff first, and then **A** responds?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min(max(Pwin_summary(A, B) for A in cutoffs)\n",
    "    for B in cutoffs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In both cases, the rational choice for both players in a cutoff of 0.618034, which corresponds to the  \"saddle point\" in the middle of the plot. This is a *stable equilibrium*; consider fixing *B* = 0.618034, and notice that if *A* changes to any other value, we slip off the saddle to the right or left, resulting in a worse win probability for **A**. Similarly, if we fix *A* = 0.618034, then if *B* changes to another value, we ride up the saddle to a higher win percentage for **A**, which is worse for **B**. So neither player will want to move from the saddle point.\n",
    "\n",
    "The moral for continuous spaces is the same as for discrete spaces: be careful about defining your sample space; measure carefully, and let your code take care of the rest."
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}