prob.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sets, Combinatorics, & Probability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probability versus Statistics\n",
    "\n",
    "In probability, we assume we know what the world is like, and predict what we would then see.\n",
    "\n",
    "In statistics, we make observations and use them to learn something about the world.\n",
    "That generally requires a mathematical model for how the obervations were generated, to relate the observations to properties of the world.\n",
    "Sometimes, that model arises naturally from knowledge of the process that generated the data,\n",
    "for instance, a randomized experiment or a physical measurement where the\n",
    "underlying physics is understood.\n",
    "In other cases, the model is basically an assumption.\n",
    "When the model does not have a firm basis in actual knowledge, inferences are suspect.\n",
    "As Robert L. Parker says, \"the more you assume, the less you know.\"\n",
    "\n",
    "Probability is like a _forward problem_ in science or engineering: the the physics, including any parameters,\n",
    "is known. The issue is to predict what will be observed, up to instrumental error and other noise.\n",
    "\n",
    "Statistics, at least inferential statistics, is like an _inverse problem_. The physics is known\n",
    "but there are unknown parameters. The issue is to use the (noisy) data to learn about parameters\n",
    "of the system.\n",
    "\n",
    "Probability is a combination of mathematics and philosophy.\n",
    "The philosophy connects the mathematics to the world.\n",
    "\n",
    "We will start with the mathematics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sets\n",
    "\n",
    "The mathematics of probability is expressed most naturally using set theory.\n",
    "We will review the basic terminology and reviews naive set theory: how to define and\n",
    "manipulate sets, operations on sets that yield other sets, special\n",
    "relationships among sets, and so on.\n",
    "\n",
    "A _set_ is a collection of objects, called _members_ or _elements_ of the set, without\n",
    "regard for their order.\n",
    "$a \\in A$, pronounced \"$a$ is an element of $A$,\"  \"$a$ is in $A$,\" or \"$a$ is a member of\n",
    "$A$\" means that $a$ is an element of the set $A$.\n",
    "This is the same as writing $A \\ni a$, which is pronounced \"$A$ contains $a$.\"\n",
    "If $a$ is not an element of $A$, we write $a \\not \\in A$.\n",
    "Sets may be described explicitly by listing their contents, or implicitly\n",
    "by specifying a property that all elements of the set share, or a condition\n",
    "that they satisfy.\n",
    "The contents of sets are enclosed in curly braces: $\\{ \\}$.\n",
    "\n",
    "Here are some examples:\n",
    "+ $A = \\{a, b, c, d \\}$: the set containing the four elements $a$, $b$, $c$, and $d$.\n",
    "+ $\\emptyset = \\{ \\}$: the empty set, the set that contains no elements.\n",
    "+ ${\\mathbf Z}  \\equiv \\{ \\ldots, -2, -1, 0, 1, 2, \\ldots \\}$: the integers.\n",
    "+ ${\\mathbf N}  \\equiv \\{1, 2, 3, \\ldots \\}$: the natural (counting) numbers.\n",
    "+ $\\Re \\equiv (-\\infty, \\infty)$: the real numbers.\n",
    "+ $\\Re^+   \\equiv [-\\infty, \\infty]$: the extended real numbers.\n",
    "+ ${\\mathbf C}  \\equiv \\{ a + bi: a, b \\in \\Re \\}$, where $i = \\sqrt{-1}$: the complex numbers.\n",
    "+ ${\\mathbf Q}  \\equiv \\{ a/b: a, b \\in {\\mathbf Z}, b \\ne 0\\}$: the rational numbers.\n",
    "\n",
    "\n",
    "$B$ is a _subset_ of $A$, written $B \\subset A$ or $A \\supset B$, if every element\n",
    "of the set $B$ is also an element of the set $A$.\n",
    "Thus ${\\mathbf N}  \\subset {\\mathbf Z}  \\subset {\\mathbf Q}  \\subset \\Re \\subset {\\mathbf C}$.\n",
    "The empty set $\\emptyset$ is a subset of every set.\n",
    "If $A \\subset B$ and $B \\subset A$, $A$ and $B$ are the same set, and we write $A = B$.\n",
    "If $B$ is not a subset of $A$, we write $B \\not \\subset A$ or $A \\not \\supset B$.\n",
    "$B$ is a _proper subset_ of $A$ if $B \\subset A$ but\n",
    "$A \\not \\subset B$.\n",
    "\n",
    "The _complement_ of $A$ (with respect to the universe ${\\mathcal X}$), written $A^c$ or $A'$,\n",
    "is the set of all objects under consideration (${\\mathcal X}$) that are not elements of $A$.\n",
    "That is, $A^c \\equiv \\{ a \\in {\\mathcal X} : a \\not \\in A\\}$.\n",
    "\n",
    "The _intersection_ of $A$ and $B$, written $A \\cap B$ or $AB$, is the set of all objects that\n",
    "are elements of both $A$ and $B$:\n",
    "$$\n",
    "A \\cap B \\equiv \\{a: a \\in A \\mbox{ and } a \\in B \\}.\n",
    "$$\n",
    "If $A \\cap B = \\emptyset$, we say $A$ and $B$ are _disjoint_ or _mutually\n",
    "exclusive_.\n",
    "\n",
    "The _union_ of $A$ and $B$, written $A \\cup B$, is the set of all objects that\n",
    "are elements of $A$ or of $B$ (or both):\n",
    "$$\n",
    "A \\cup B \\equiv \\{a: a \\in A \\mbox{ or } a \\in B \\mbox{ or both} \\}.\n",
    "$$\n",
    "\n",
    "The _difference_ or _set difference_ of $A$ and $B$, $A \\setminus B$, pronounced \"$A$ minus $B$,\" is\n",
    "the set of all elements of $A$ that are not elements of $B$:\n",
    "$$\n",
    "A \\setminus B \\equiv \\{a \\in A : a \\not \\in B \\} = A \\cap B^c.\n",
    "$$\n",
    "\n",
    "_Intervals_ are special subsets of ${\\mathbf R}$:\n",
    "$$\n",
    "[a, b] \\equiv \\{x \\in \\Re : a \\le x \\le b\\}\n",
    "$$\n",
    "$$\n",
    "(a, b] \\equiv \\{x \\in \\Re : a < x \\le b\\}\n",
    "$$\n",
    "$$\n",
    "[a, b) \\equiv \\{x \\in \\Re : a \\le x < b\\}\n",
    "$$\n",
    "$$\n",
    "(a, b) \\equiv \\{x \\in \\Re : a < x < b\\}.\n",
    "$$\n",
    "Sometimes we have a collection of sets, _indexed_ by elements of another\n",
    "set: $\\{A_\\beta : \\beta \\in B \\}$.\n",
    "Then $B$ is called an _index set_.\n",
    "If $B$ is a subset of the integers ${\\mathbf Z}$, usually we write $A_i$ or $A_j$, etc.,\n",
    "rather than $A_\\beta$.\n",
    "If $B = {\\mathbf N}$, we usually write $\\{A_j\\}_{j=1}^\\infty$ rather than\n",
    "$\\{A_\\beta : \\beta \\in {\\mathbf N}  \\}$.\n",
    "$$\n",
    "\\bigcap_{\\beta \\in B} A_\\beta \\equiv \\{a: a \\in A_\\beta \\;\\;\\forall \\beta \\in B \\}.\n",
    "$$\n",
    "($\\forall$ means \"for all.\")\n",
    "If $B = \\{1, 2, \\ldots, n\\}$, we usually write $\\bigcap_{j=1}^n A_j$ rather than\n",
    "$\\bigcap_{j \\in \\{1, 2, \\ldots, n\\}} A_j$.\n",
    "The notation $\\bigcup_{\\beta \\in B} A_\\beta$ and $\\bigcup_{j=1}^n A_j$ are defined analogously.\n",
    "\n",
    "A collection of sets $\\{A_\\beta : \\beta \\in B \\}$ is _pairwise disjoint_\n",
    "if $A_\\beta \\cap A_{\\beta'} = \\emptyset$ whenever $\\beta \\ne \\beta'$.\n",
    "The collection $\\{A_\\beta : \\beta \\in B\\}$ _exhausts_ or _covers_\n",
    "the set $A$ if $A \\subset \\bigcup_{\\beta \\in B} A_\\beta$.\n",
    "The collection $\\{A_\\beta : \\beta \\in B \\}$ is a _partition_ of the\n",
    "set $A$ if $A = \\cup_{\\beta \\in B} A_\\beta$ and the \n",
    "sets $\\{A_\\beta : \\beta \\in B \\}$\n",
    "are pairwise disjoint.\n",
    "If $\\{A_\\beta : \\beta \\in B \\}$ are pairwise disjoint and exhaust $A$, then\n",
    "$\\{A_\\beta \\cap A : \\beta \\in B \\}$ is a partition of $A$.\n",
    "\n",
    "A set is _countable_ if its elements can be put in one-to-one correspondence with\n",
    "a subset of ${\\mathbf N}$.\n",
    "A set is _finite_ if its elements can be put in one-to-one correspondence with\n",
    "$\\{1, 2, \\ldots, n\\}$ for some $n \\in {\\mathbf N}$.\n",
    "If a set is not finite, it is infinite.\n",
    "${\\mathbf N}$, ${\\mathbf Z}$, and ${\\mathbf Q}$ are infinite and countable\n",
    "(_countably infinite_);\n",
    "${\\mathbf R}$ is infinite and uncountable.\n",
    "\n",
    "The notation $\\# A$, pronounced \"the cardinality of $A$\" is the size of the set $A$, suitably defined.\n",
    "If $A$ is finite, $\\# A$ is the number of elements in $A$.\n",
    "If $A$ is not finite but $A$ is countable (if its elements can be put in one-to-one\n",
    "correspondence with the elements of ${\\mathbf N}$), then $\\# A = \\aleph_0$ (aleph-null).\n",
    "If the elements of $A$ can be put in one-to-one correspondence with the elements of $\\mathbf{R}$, then\n",
    "$\\#A = c$, _the cardinality of the continuum_.\n",
    "\n",
    "The _power set_ of a set $A$, denoted $2^A$, is the set of all subsets of the set $A$.\n",
    "For example, the power set of $\\{a, b, c\\}$ is\n",
    "$$\n",
    "\\{ \\emptyset, \\{a\\}, \\{b\\}, \\{c\\}, \\{a, b\\}, \\{a, c\\}, \\{b, c\\}, \\{a, b, c\\} \\}.\n",
    "$$\n",
    "If $A$ is a finite set, the cardinality of the power set of $A$ is $2^{\\# A}$.\n",
    "This can be seen as follows: suppose $\\# A = n$ is finite.\n",
    "Consider the elements of $A$ to be written in some canonical order.\n",
    "We can specify an element of the power set by an $n$-digit binary number.\n",
    "The first digit is 1 if the first element of $A$ is in the subset, and 0 otherwise.\n",
    "The second digit is 1 if the second element of $A$ is in the subset, and 0 otherwise,\n",
    "etc.\n",
    "There are $2^n$ $n$-digit binary numbers, so there are $2^n$ subsets.\n",
    "The cardinality of the power set of ${\\mathbf N}$ is not $\\aleph_0$.\n",
    "\n",
    "If $A$ is a finite set, $B$ is a countable set\n",
    "and $\\{A_j : \\beta \\in B \\}$ is a partition of $A$, then\n",
    "$$\n",
    "\\# A = \\sum_{\\beta \\in B} \\# A_\\beta.\n",
    "$$\n",
    "\n",
    "Suppose $A = A_1 \\cup A_2$, but that possibly $A_1 A_2 \\ne \\emptyset$,\n",
    "so $\\{A_1, A_2\\}$ might not be a partition of $A$, because\n",
    "$A_1$ and $A_2$ are not necessarily disjoint.\n",
    "Then still \n",
    "$$\\#A = \\#A_1 + \\#A_2 - \\#A_1A_2.$$\n",
    "\n",
    "This is seen most easily using a Venn diagram, and can be proved\n",
    "by constructing a partition of $A$,\n",
    "$A = A_1A_2^c \\cup A_1^cA_2 \\cup A_1A_2$, and noting that\n",
    "$\\#A_1 + \\#A_2 = \\#A_1A_2^c + \\#A_1^cA_2 + 2\\#A_1A_2$.\n",
    "\n",
    "If \n",
    "$A = A_1 \\cup A_2 \\cup A_3$ but $\\{A_1, A_2, A_3\\}$ are not necessarily disjoint,\n",
    "then still\n",
    "\n",
    "$$\\#A = \\#A_1 + \\#A_2 + \\#A_3 - \\#A_1A_2 - \\#A_1A_3 - \\#A_2A_3 + \\#A_1A_2A_3.$$\n",
    "\n",
    "More generally, if $A = \\cup_{i=1}^n A_i$, then the _Inclusion-Exclusion Principle_ holds:\n",
    "\n",
    "$$ \\#A = \\sum_{i=1}^n \\#A_i - \\sum_{1 \\le i_1 < i_2 \\le n} \\#(A_{i_1}A_{i_2}) +\n",
    "\\sum_{1 \\le i_1 < i_2 < i_3 \\le n} \\#(A_{i_1}A_{i_2}A_{i_3}) - \\cdots \n",
    "+(-1)^{k-1} \\sum_{1 \\le i_1 < i_2 < \\cdots < i_k \\le n} \\# (A_{i_1}A_{i_2} \\cdots A_{i_k}) + \\cdots\n",
    "$$\n",
    "\n",
    "### Some useful set relations\n",
    "\n",
    "+ If $A\\subset B$ and $B\\subset A$ then $A = B$. (If two sets are subsets of each other, they are the same set.)\n",
    "\n",
    "+ $\\emptyset \\subset A$. (The empty set is a subset of every set.)\n",
    "\n",
    "+ $\\emptyset \\cup A = A$.\n",
    "(The union of the empty set and any other set is that set.)\n",
    "\n",
    "+ $\\emptyset \\cap A = \\emptyset$.\n",
    "(The intersection of the empty set and any other set is empty.)\n",
    "\n",
    "+ If $A \\subset B$ and $B \\subset C$ then $A \\subset C$. (The subset relation is transitive.)\n",
    "\n",
    "+ If $A \\subset B$ then $B^c \\subset A^c$.\n",
    "(Complementation reverses the subset relation.)\n",
    "\n",
    "+ $A \\cap B \\subset  A$.  Moreover, $A\\cap B = A$ if and only if $A \\subset B$.\n",
    "\n",
    "+ $A \\subset  A\\cup B$.   Moreover, $A = A\\cup B$ if and only if $B \\subset A$.\n",
    "\n",
    "+ $(A\\cap B)^c = A^c \\cup B^c$. (de Morgan)\n",
    "\n",
    "+ $(A\\cup B)^c = A^c\\cap B^c$. (de Morgan)\n",
    "\n",
    "+ $A \\cap B = B \\cap A$. (Intersection is commutative.)\n",
    "\n",
    "+ $A\\cup B = B\\cup A$. (Union is commutative.)\n",
    "\n",
    "+ $A\\cap (B\\cap C) = (A\\cap B)\\cap C$. (Intersection is associative.)\n",
    "\n",
    "+ $A\\cup (B\\cup C) = (A\\cup B)\\cup C$. (Union is associative.)\n",
    "\n",
    "+ $A\\cap (B\\cup C) = (A\\cap B)\\cup (A\\cap C)$. (Distribution of intersection over union.)\n",
    "\n",
    "+ $A\\cup (B\\cap C) = (A\\cup B)\\cap (A\\cup C)$. (Distribution of union over intersection.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "1. Prove that if $\\{A_\\beta : \\beta \\in B \\}$ are pairwise disjoint and exhaust $A$, then $\\{A_\\beta \\cap A : \\beta \\in B \\}$ is a partition of $A$.\n",
    "\n",
    "2. Prove de Morgan's identity $(A\\cup B)^c = A^c\\cap B^c$.\n",
    "\n",
    "3. Prove that $\\{A_1A_2^c, A_1^cA_2, A_1A_2\\}$ is a partition of $A_1 \\cup A_2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# boilerplate\n",
    "%matplotlib inline\n",
    "from __future__ import division\n",
    "import math\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy import special\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A:  [1, 2, 3, 3]\n",
      "A as a set (duplicates removed): set([1, 2, 3])\n",
      "B: set(['a', 3, 'b', 'green'])\n",
      "C: set([1, 2])\n",
      "1 is in A: True\n",
      "\"green\" is in A: False\n",
      "\"green\" is in B: True\n",
      "A is a subset of C: False\n",
      "C is a subset of A: True\n",
      "A intersect B: set([3])\n"
     ]
    }
   ],
   "source": [
    "# Sets in Python\n",
    "# Make three sets, A = {1, 2, 3}, B = {\"a\", \"b\", \"green\", 3}, C = {1, 2}\n",
    "A = [1, 2, 3, 3] # this is a list but not a set, because 3 is listed twice\n",
    "print 'A: ', A\n",
    "A = set(A)  # this is a set containing the unique elements of A\n",
    "print 'A as a set (duplicates removed):',A\n",
    "B = set([\"a\", \"b\", \"green\", 3])\n",
    "print 'B:', B\n",
    "C = set(range(1,3)) # a different construction\n",
    "print 'C:', C\n",
    "#\n",
    "# Set membership\n",
    "print '1 is in A:', 1 in A  # is 1 an element of A?\n",
    "print '\"green\" is in A:',\"green\" in A\n",
    "print '\"green\" is in B:',\"green\" in B\n",
    "print 'A is a subset of C:', A.issubset(C)\n",
    "print 'C is a subset of A:', C.issubset(A)\n",
    "print 'A intersect B:', A.intersection(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cartesian Products\n",
    "The notation $(s_1, s_2, \\ldots, s_n) \\equiv (s_j)_{j=1}^n$ denotes an\n",
    "_ordered $n$-tuple_ consisting of\n",
    "$s_1$ in the first position, $s_2$ in the second position, etc.\n",
    "The parentheses are used instead of curly braces to distinguish\n",
    "$n$-tuples from sets: $(s_j)_{j=1}^n \\ne \\{s_j\\}_{j=1}^n$.\n",
    "The $k$th\n",
    "_component_ of the $n$-tuple $s = (s_j)_{j=1}^n$, is $s_k$,\n",
    "$k = 1, 2, \\ldots, n$.\n",
    "Two $n$-tuples are equal if their components are equal.\n",
    "That is, $(s_j)_{j=1}^n = (t_j)_{j=1}^n$ means that\n",
    "$s_j = t_j$ for $j = 1, \\ldots, n$.\n",
    "In particular, $(s, t) \\ne (t, s)$ unless $s=t$.\n",
    "In contrast, $\\{s, t \\} = \\{ t, s \\}$ always.\n",
    "\n",
    "The _Cartesian product of $S$ and $T$_ is\n",
    "$S \\bigotimes T \\equiv \\{(s, t): s \\in S \\mbox{ and } t \\in T\\}$.\n",
    "Unless $S = T$, $S \\bigotimes T \\ne T \\bigotimes S$.\n",
    "${\\mathbf R}^n$ is the Cartesian product of ${\\mathbf R}$ with itself, $n$ times; its elements\n",
    "are $n$-tuples of real numbers.\n",
    "If $s$ is the $n$-tuple $(s_1, s_2, \\ldots, s_n) = (s_j)_{j=1}^n$\n",
    "and $t$ is the $n$-tuple $(t_1, t_2, \\ldots, t_n) = (t_j)_{j=1}^n$,\n",
    "then $s = t$ iff $s_j = t_j$ for all $j=1, \\ldots, n$.\n",
    "\n",
    "The Python library <tt>itertools</tt> has utilities for generating Cartesian products, power sets, permutations, and the like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('a', 0), ('a', 1), ('a', 2), ('b', 0), ('b', 1), ('b', 2), ('c', 0), ('c', 1), ('c', 2), ('d', 0), ('d', 1), ('d', 2)]\n"
     ]
    }
   ],
   "source": [
    "# Cartesian products using itertools\n",
    "import itertools\n",
    "A = ['a','b','c','d']\n",
    "B = range(3)\n",
    "print [i for i in itertools.product(A,B)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combinatorics\n",
    "Let $A$ be a finite set with $\\# A = n$.\n",
    "A _permutation_ of (the elements of) $A$ is an element $s$ of\n",
    "$\\bigotimes_{j=1}^n A = A^n$\n",
    "whose components are distinct elements of $A$.\n",
    "That is, $s = (s_j)_{j=1}^n \\in A^n$ is a permutation of $A$ if\n",
    "$\\# \\{s_j\\}_{j=1}^n = n$.\n",
    "There are $n! = n(n-1)\\cdots 1$ permutations of a set with $n$ elements:\n",
    "there are $n$ choices for the first component of the permutation, $n-1$ choices for\n",
    "the second (whatever the\n",
    "first might be), $n-2$ for the third (whatever the first two might be), etc.\n",
    "\n",
    "This is an illustration of the _Fundamental Rule of Counting_:\n",
    "If specifying each (unique) element of a set can be expressed as a series $n$ of choices, and the number of options at step $i$ of the series is $n_i$ regardless of the choices made before step $i$, then the total number of elements of the set is $\\prod_{i=1}^n n_i$. \n",
    "\n",
    "The number of permutations of $n$ things taken $k$ at a time, ${}_nP_k$,\n",
    "is the number of ways there are of selecting $k$ of $n$ things, then permuting\n",
    "those $k$ things.\n",
    "There are $n$ choices for the object that will be in the first place in the permutation,\n",
    "$n-1$ for the second place (regardless of which is first), etc., and $n-k+1$ choices\n",
    "for the item that will be in the $k$th place.\n",
    "By the fundamental rule of counting, it follows that\n",
    "${}_nP_k = n(n-1)\\cdots(n-k+1) = n!/(n-k)!$.\n",
    "\n",
    "The number of subsets of size $k$ that can be formed from $n$ objects is\n",
    "$$\n",
    "{}_nC_k = {{n}\\choose{k}} =\n",
    "{}_nP_k/k! = n(n-1)\\cdots(n-k+1)/k! = \\frac{n!}{k!(n-k)!},\n",
    "$$\n",
    "since each subset of size $k$ can be ordered into $k!$ permutations.\n",
    "\n",
    "Here are some useful identities:\n",
    "\n",
    "1. $ {{n}\\choose{k}} = {{n}\\choose{n-k}} $.\n",
    "1. $ {{n}\\choose{0}} = {{n}\\choose{n}} = 1$.\n",
    "1. $ {{n}\\choose{1}} = {{n}\\choose{n-1}} = n $.\n",
    "1. $ {{n} \\choose {k}} = \\prod_{m=0}^{k-1} \\frac{n-m}{m+1} $.\n",
    "1. $ {{n} \\choose {k}} = {{n-1}\\choose{k-1}} + {{n-1} \\choose {k}}$.\n",
    "\n",
    "Because the power set of a set with $n$ elements can be partitioned as\n",
    "\n",
    "$$\n",
    "\\cup_{k=0}^n \\left \\{ \\mbox{all distinct subsets of size $k$}\n",
    "\\right \\},\n",
    "$$\n",
    "\n",
    "and since the power set has $2^n$ elements, it follows that\n",
    "$$\n",
    "\\sum_{j=0}^n {{n} \\choose {j}} = 2^n.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr />\n",
    "### Exercise\n",
    "Prove the five identities about binomial coefficients given above.\n",
    "<hr />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n!: 120\n",
      "nCk: 10.0\n",
      "nPk: 60\n"
     ]
    }
   ],
   "source": [
    "# Combinatorics in Python\n",
    "\n",
    "n = 5\n",
    "k = 3\n",
    "print 'n!:',math.factorial(n)\n",
    "\n",
    "# number of combinations of n objects taken k at a time\n",
    "print 'nCk:',sp.special.binom(n,k)\n",
    "\n",
    "# number of permutations of n objects taken k at a time\n",
    "\n",
    "def permu(n, k):\n",
    "    permu = 1\n",
    "    s = 0\n",
    "    while (s < k):\n",
    "        permu = permu*(n-s)\n",
    "        s = s+1\n",
    "    return(permu)\n",
    "\n",
    "print 'nPk:',permu(n, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A: ['a', 'b', 'c', 'd', 'e']\n",
      "combinations of  2 elements of A: [('a', 'b'), ('a', 'c'), ('a', 'd'), ('a', 'e'), ('b', 'c'), ('b', 'd'), ('b', 'e'), ('c', 'd'), ('c', 'e'), ('d', 'e')]\n",
      "permutations of  2 elements of A: [('a', 'b'), ('a', 'c'), ('a', 'd'), ('a', 'e'), ('b', 'a'), ('b', 'c'), ('b', 'd'), ('b', 'e'), ('c', 'a'), ('c', 'b'), ('c', 'd'), ('c', 'e'), ('d', 'a'), ('d', 'b'), ('d', 'c'), ('d', 'e'), ('e', 'a'), ('e', 'b'), ('e', 'c'), ('e', 'd')]\n"
     ]
    }
   ],
   "source": [
    "# enumerate combinations and permutations using itertools\n",
    "A = ['a','b','c','d','e']\n",
    "print 'A:', A\n",
    "k = 2\n",
    "print 'all combinations of', k, 'elements of A:', [i for i in itertools.combinations(A,k)]\n",
    "print 'all permutations of', k, 'elements of A:', [i for i in itertools.permutations(A,k)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Strategies for counting\n",
    "\n",
    "There are several approaches to counting the elements of a set that can be helpful in practice, aside from simply enumerating the elements.\n",
    "\n",
    "1. Divide and conquer: partition the set.\n",
    "1. Use the Fundamental Rule of Counting. (This is how we found the number of permutations of $k$ of $n$ things.)\n",
    "1. Overcount: count each item $k$ times, then divide by $k$. (This is how we derived the number of combinations of $k$ of $n$ things from the number of permutations of $k$ of $n$ things: each combination of $k$ things occurs $k!$ times among the set of permutations.)\n",
    "1. Overinclude, then adjust: count a set that is larger than the set in question, then subtract the extras.\n",
    "1. Use the Inclusion-Exclusion Principle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computational exercise\n",
    "\n",
    "Write a Python function that constructs all permutations of a finite list.\n",
    "\n",
    "The function should take as input a list s, and produce as output all permutations of the elements of the list. Items should be considered unique based on their position in the list.\n",
    "\n",
    "For instance,\n",
    "\n",
    "    s = ['a','b','c']\n",
    "    permute(s)\n",
    "\n",
    "should return the 6 permutations of a, b, and c:\n",
    "\n",
    "    a b c\n",
    "    a c b\n",
    "    b a c\n",
    "    b c a\n",
    "    c a b\n",
    "    c b a\n",
    "\n",
    "Write your function from scratch; i.e., do not use any libraries that are not part of the core of Python (e.g., do not use the <span class=\"code\">itertools</span> library).\n",
    "\n",
    "<hr />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical stability and order of operations\n",
    "\n",
    "As a matter of mathematics, ${{n}\\choose{k}} = \\frac{n!}{k!(n-k)!}$, but this is a bad way to compute ${{n}\\choose{k}}$.\n",
    "When $n$ is large, the numerator $n!$ is enormous, and can overflow finite-precision calculations even for values of $n$ and $k$ for which ${{n}\\choose{k}}$ is small.\n",
    "Moreover, relying on cancellation in finite-precision arithmetic can lead to large errors, whether the cancellation is through division or subtraction.\n",
    "\n",
    "So, rather than compute ${{n}\\choose{k}}$ by dividing factorials, it's better to build it by multiplying ratios."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000!: 402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000\n",
      "\n",
      "998!: 402790050127220994538240674597601587306681545756471103647447357787726238637266286878923131618587992793273261872069265323955622495490298857759082912582527118115540044131204964883707335062250983503282788739735011132006982444941985587005283378024520811868262149587473961298417598644470253901751728741217850740576532267700213398722681144219777186300562980454804151705133780356968636433830499319610818197341194914502752560687555393768328059805942027406941465687273867068997087966263572003396240643925156715326363340141498803019187935545221092440752778256846166934103235684110346477890399179387387649332483510852680658363147783651821986351375529220618900164975188281042287183543472177292257232652561904125692525097177999332518635447000616452999984030739715318219169707323799647375797687367013258203364129482891089991376819307292252205524626349705261864003453853589870620758596211518646408335184218571196396412300835983314926628732700876798309217005024417595709904449706930796337798861753941902125964936412501007284147114260935633196107341423863071231385166055949914432695939611227990169338248027939843597628903525815803809004448863145157344706452445088044626373001304259830129153477630812429640105937974761667785045203987508259776060285826091261745049275419393680613675366264232715305430889216384611069135662432391043725998805881663054913091981633842006354699525518784828195856033032645477338126512662942408363494651203239333321502114252811411713148843370594801145777575035630312885989779863888320759224882127141544366251503974910100721650673810303577074640154112833393047276025799811224571534249672518380758145683914398263952929391318702517417558325636082722982882372594816582486826728614633199726211273072775131325222240100140952842572490801822994224069971613534603487874996852498623584383106014533830650022411053668508165547838962087111297947300444414551980512439088964301520461155436870989509667681805149977993044444138428582065142787356455528681114392680950815418208072393532616122339434437034424287842119316058881129887474239992336556764337968538036861949918847009763612475872782742568849805927378373244946190707168428807837146267156243185213724364546701100557714520462335084082176431173346929330394071476071813598759588818954312394234331327700224455015871775476100371615031940945098788894828812648426365776746774528000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000\n",
      "\n",
      "1000C2: 499500.0\n",
      "1000C2 (2nd way): 499500.0\n"
     ]
    }
   ],
   "source": [
    "# illustrating issues calculating nCk as a ratio of factorials for large n\n",
    "n = 1000\n",
    "k = 2\n",
    "print '1000!:', math.factorial(n) # Really Big Number: \n",
    "                                  # would overflow in many languages, \n",
    "                                  # including R, MATLAB, Fortran\n",
    "# however, 1000 choose 2 is only (1000*999)/2 = 499500\n",
    "print '\\n998!:', math.factorial(n-k)\n",
    "print '\\n1000C2:', sp.special.binom(n, k)\n",
    "\n",
    "# let's code this from scratch in a stable way\n",
    "def myChoose(n, k):\n",
    "    ''' \n",
    "    no error checking. If this were for production, we'd check that the arguments are integers\n",
    "    and that n >= k >= 0.\n",
    "    '''\n",
    "    m = 0\n",
    "    myc = 1\n",
    "    while (m < k):\n",
    "        myc = myc * (n-m)/(m+1)\n",
    "        m = m+1\n",
    "    return(myc)\n",
    "\n",
    "print '1000C2 (2nd way):', myChoose(n, k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A related same issue comes up in subtracting large numbers.\n",
    "\n",
    "For instance, algebraically, $x^2 - (x-1)^2 = 2x-1$. \n",
    "But if the mantissa of $x$ is large, squaring\n",
    "$x$ and $x-1$ can overflow the precision of the calculation, and their difference of the computed squares will be numerically wrong."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Connecting Probability to Set Theory\n",
    "\n",
    "A _random experiment_ or _random trial_\n",
    "is basically any\n",
    "situation whose outcome is not perfectly predictable, but for which we\n",
    "can specify all possible outcomes, and that shows long-term regularities.\n",
    "For example, when we toss a coin, we do not know how it will land,\n",
    "but it certainly must land heads, tails, on its edge, or not land at all.\n",
    "There is no other possibility.\n",
    "The set of all possible outcomes of a random experiment is called the\n",
    "_outcome space_.\n",
    "The letter $\\mathcal{S}$ will denote outcome space.\n",
    "We are free to choose the outcome space to correspond to what we deem\n",
    "relevant for the experiment, as long as it is essentially inevitable that the\n",
    "random experiment will result in some outcome in the outcome space.\n",
    "For example, the outcome space we just described was {heads, tails, edge, doesn't land}.\n",
    "It might be adequate for our purposes for the outcome space to be {heads, not heads}.\n",
    "\n",
    "Often we shall tailor outcome spaces for specific problems.\n",
    "Here is an example: Imagine a box containing tickets that are indistinguishable\n",
    "except that each has written upon it a unique number between 1 and the number\n",
    "of tickets, $n$.\n",
    "We shake the box, draw a ticket from the box without looking, and record\n",
    "the number written on the ticket we happened to draw.\n",
    "The natural outcome space of this experiment is the set of numbers\n",
    "$\\{1, 2, \\ldots, n\\}$.\n",
    "However, suppose we are interested only in whether the number on the\n",
    "ticket we draw is even.\n",
    "The outcome space then could be reduced to {even number on ticket, odd number on ticket},\n",
    "or coded even more abstractly as $\\{0, 1\\}$, \n",
    "where the outcome is the number of even-numbered tickets drawn.\n",
    "\n",
    "An _event_ is a subset of outcome space: a collection of outcomes in the\n",
    "outcome space.\n",
    "$A$ is an event if $A \\subset \\mathcal{S}$.\n",
    "For example, in the experiment of drawing a numbered ticket from the box,\n",
    "suppose there are 10 tickets in all, and that we choose the outcome space\n",
    "$\\mathcal{S}$ to be the numbers\n",
    "$\\{1, 2, 3, \\ldots, 9, 10\\}$.\n",
    "Then &quot;we draw the number 1&quot; is the event $\\{1\\}$, and\n",
    "&quot;we draw an even number&quot; is the event $\\{2, 4, 6, 8, 10\\}$,\n",
    "both of which are subsets of the set of possible outcomes, the outcome space\n",
    "$\\mathcal{S}$.\n",
    "\n",
    "The treatment is quite informal, omitting important technical concepts such as sigma-algebras,\n",
    "measurability, and the like.\n",
    "\n",
    "Modern probability theory starts with Kolmogorov's Axioms; here is an informal startement of the axioms.\n",
    "For more (but still a very informal treatment), see: \n",
    "[Probability: Philosophy and Mathematical Background](http://www.stat.berkeley.edu/~stark/SticiGui/Text/probabilityPhilosophy.htm),\n",
    "[Set theory](http://www.stat.berkeley.edu/~stark/SticiGui/Text/sets.htm),\n",
    "and\n",
    "[Probability: Axioms and Fundaments](http://www.stat.berkeley.edu/~stark/SticiGui/Text/probabilityAxioms.htm).\n",
    "\n",
    "Let $\\mathcal{S}$ denote the _outcome space_, the set of all possible outcomes of a random experiment, \n",
    "and let $\\{A_i\\}_{i=1}^\\infty$ be a countable collection of \n",
    "subsets of $\\mathcal{S}$.\n",
    "Then any probability function ${\\mathbb P}$ must satisfy these axioms:\n",
    "\n",
    "1. For every $A \\subset \\mathcal{S}$, ${\\mathbb P}(A) \\ge 0$ (probabilities are nonnegative)\n",
    "2. ${\\mathbb P}(\\mathcal{S}) = 1$ (the chance that _something_ happens is 100%)\n",
    "3. If $A_i \\cap A_j = \\emptyset$ for $i \\ne j$, then ${\\mathbb P} \\cup_{i=1}^\\infty A_i = \\sum_{i=1}^\\infty {\\mathbb P}(A_i)$\n",
    "(If a countable collection of events is _pairwise disjoint_, then the chance that any of the\n",
    "events occurs is the sum of the chances that they occur individually.)\n",
    "\n",
    "These axioms have many useful consequences, among them that ${\\mathbb P}(\\emptyset) = 0$, ${\\mathbb P}(A^c) = 1 - {\\mathbb P}(A)$,\n",
    "and ${\\mathbb P}(A \\cup B) = {\\mathbb P}(A) + {\\mathbb P}(B) - {\\mathbb P}(AB)$.\n",
    "More generally, the Inclusion-Exclusion Principle gives us\n",
    "\n",
    "$$\n",
    "\\mathbb{P} \\cup_{i=1}^n A_i  = \\sum_{i=1}^n \\mathbb{P} A_i \n",
    "- \\sum_{1 \\le i_1 < i_2 \\le n} \\mathbb{P}(A_{i_1}A_{i_2}) +\n",
    "\\sum_{1 \\le i_1 < i_2 < i_3 \\le n} \\mathbb{P}(A_{i_1}A_{i_2}A_{i_3}) - \\cdots \n",
    "+(-1)^{k-1} \\sum_{1 \\le i_1 < i_2 < \\cdots < i_k \\le n} \\mathbb{P} (A_{i_1}A_{i_2} \\cdots A_{i_k}) + \\cdots\n",
    "$$\n",
    "$$\n",
    "= \\sum_{i=1}^n \\left((-1)^{i-1}\\sum_{\\mathcal{I} \\subset\\{1,\\ldots,n\\}: \\#\\mathcal I =i} \\mathbb{P}(\\cap_{j \\in \\mathcal{I}} A_j \\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr />\n",
    "### Definitions\n",
    "\n",
    "Let $A$ and $B$ be subsets of outcome space $\\mathcal{S}$.\n",
    "\n",
    "+ If $AB = \\emptyset$, then $A$ and $B$ are _mutually exclusive_.\n",
    "+ If $A \\subset B$ then $A$ _implies_ $B$:  if $A$ occurs, $B$ must\n",
    "also occur (if the outcome that occurs is in $A$, the outcome\n",
    "that occurs is also in $B$, because every element of $A$\n",
    "is an element of $B$).\n",
    "+ If ${\\mathbb P}(AB) = {\\mathbb P}(A){\\mathbb P}(B)$, then $A$ and $B$ are _independent_.\n",
    "+ If ${\\mathbb P}(B) > 0$, then the _conditional probability of $A$ given $B$_ is\n",
    "${\\mathbb P}(A | B) \\equiv {\\mathbb P}(AB)/{\\mathbb P}(B)$.\n",
    "<hr />\n",
    "\n",
    "Independence is an extremely specific relationship.\n",
    "At one extreme, $AB = \\emptyset$; then ${\\mathbb P}(AB) = 0 \\le {\\mathbb P}(A){\\mathbb P}(B)$. \n",
    "At another extreme, either $A$ is a subset of $B$ or vice versa; then ${\\mathbb P}(AB) = \\min({\\mathbb P}(A),{\\mathbb P}(B)) \\ge {\\mathbb P}(A){\\mathbb P}(B)$.\n",
    "Independence lies at a precise point in between."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises.\n",
    "\n",
    "1. Show that if $A$ and $B$ are independent and $\\mathbb{P}(B) >0$, then\n",
    "$\\mathbb{P}(A|B) = \\mathbb{P}(A)$.\n",
    "1. Show that conditional probability behaves like \"ordinary\" probability: suppose $\\mathbb{P}(B) > 0$ and $\\{A_i\\}$ are such that $\\{A_i \\cap B\\}_{i=1}^\\infty$ are pairwise disjoint. Then \n",
    "   - $\\mathbb{P}(A | B) \\ge 0$\n",
    "   - $\\mathbb{P}(B | B) = 1$\n",
    "   - $\\mathbb{P}(\\cup_{i=1}^\\infty A_i | B) = \\sum_{i=1}^\\infty \\mathbb{P} (A_i | B)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayes' Rule\n",
    "\n",
    "Suppose $\\mathbb{P}(A)$, $\\mathbb{P}(B) \\ne 0$.\n",
    "Then\n",
    "\n",
    "$$\n",
    "\\mathbb{P}(A|B) = \\frac{\\mathbb{P}(B|A)\\mathbb{P}(A)}{\\mathbb{P}(B)}.\n",
    "$$\n",
    "\n",
    "Bayes' rule lets you \"invert\" conditional probabilities to find $\\mathbb{P}(A|B)$ in terms of $\\mathbb{P}(B|A)$.\n",
    "\n",
    "\n",
    "## The Law of Total Probability\n",
    "\n",
    "Let $A$ be an event and suppose $\\{ A_i \\}$ (a finite or countable collection of sets) is a partition of outcome space $\\mathcal{S}$ such that $\\mathbb{P}(A_i) > 0$ $\\forall i$.\n",
    "Then\n",
    "\n",
    "$$\n",
    "\\mathbb{P}(A) = \\sum_i \\mathbb{P}(A|A_i) \\mathbb{P}(A_i).\n",
    "$$\n",
    "\n",
    "The Law of Total Probability is basically just an application of the definition of conditional probability and Kolmogorov's third axiom, but it's extremely useful."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examples.\n",
    "\n",
    "Here is a classic application of Bayes' Rule. Suppose that a screening test for a medical condition\n",
    "is 95% accurate, in the sense that if one has the condition, there is a 95% chance that the test will be positive, and if one does not have the condition, there is a 95% chance that the test will be negative.\n",
    "Suppose that 1% of the population actually has the condition.\n",
    "We select someone at random from the population and test him or her.\n",
    "The test result is positive.\n",
    "What's the chance that the person actually has the condition?\n",
    "\n",
    "$$ \\mathbb{P}(\\mbox{has condition}|\\mbox{tests positive})\n",
    "= \\frac{\\mathbb{P}(\\mbox{tests positive }|\\mbox{ has condition}) \\mathbb{P}(\\mbox{has condition })}{\\mathbb{P}(\\mbox{tests positive})}\n",
    "=\n",
    "\\frac{\\mathbb{P}(\\mbox{tests positive }|\\mbox{ has condition}) \n",
    "\\mathbb{P}(\\mbox{has condition})}\n",
    "{\\mathbb{P}(\\mbox{tests positive } | \\mbox{ has condition})\\mathbb{P}(\\mbox{has condition}) + \\mathbb{P}(\\mbox{tests positive } | \\mbox{ doesn't have condition})\\mathbb{P}(\\mbox{doesn't have condition})}\n",
    "=\n",
    "\\frac{0.95}{0.95 \\times 0.01 + 0.05 \\times 0.99}\n",
    "= 16.1\\%.\n",
    "$$\n",
    "(The denominator involves applying the Law of Total Probability.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises.\n",
    "[To do]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Previous: [Jupyter & Python](./jupyter.ipynb) Next: [Theories of Probability](./probTheory.ipynb)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}