elephants_soln.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Think Bayes\n",
    "\n",
    "Copyright 2018 Allen B. Downey\n",
    "\n",
    "MIT License: https://opensource.org/licenses/MIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# import classes from thinkbayes2\n",
    "from thinkbayes2 import Pmf, Cdf, Suite, Joint\n",
    "\n",
    "from thinkbayes2 import MakePoissonPmf, EvalBinomialPmf, MakeMixture\n",
    "\n",
    "import thinkplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cats and rats and elephants\n",
    "\n",
    "Suppose there are six species that might be in a zoo: lions and tigers and bears, and cats and rats and elephants.  Every zoo has a subset of these species, and every subset is equally likely.\n",
    "\n",
    "One day we visit a zoo and see 3 lions, 2 tigers, and one bear.  Assuming that every animal in the zoo has an equal chance to be seen, what is the probability that the next animal we see is an elephant?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "I'll start by enumerating all possible zoos with `itertools`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import combinations\n",
    "\n",
    "def power_set(s):\n",
    "    n = len(s)\n",
    "    for r in range(1, n+1):\n",
    "        for combo in combinations(s, r):\n",
    "            yield ''.join(combo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can enumerate only the zoos that are possible, given a set of animals known to be present."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def enumerate_zoos(all_species, present):\n",
    "    \"\"\"Enumerate all zoos that contain `present`.\n",
    "    \n",
    "    all_species: sequence of all species\n",
    "    present: sequence of species present\n",
    "    \n",
    "    yields: possible zoos\n",
    "    \"\"\"\n",
    "    present = set(present)\n",
    "    for combo in power_set(species):\n",
    "        intersect = set(combo) & present\n",
    "        if len(intersect) == len(present):\n",
    "            yield len(combo), combo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the possible zoos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3 LTB\n",
      "4 LTBC\n",
      "4 LTBR\n",
      "4 LTBE\n",
      "5 LTBCR\n",
      "5 LTBCE\n",
      "5 LTBRE\n",
      "6 LTBCRE\n"
     ]
    }
   ],
   "source": [
    "species = 'LTBCRE'\n",
    "present = 'LTB'\n",
    "\n",
    "for n, zoo in enumerate_zoos(species, present):\n",
    "    print(n, zoo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To represent the prior and posterior distributions I'll use a hierarchical model with one Dirichlet object for each possible zoo.\n",
    "\n",
    "At the bottom of the hierarchy, it is easy to update each Dirichlet object just by adding the observed frequencies to the parameters.\n",
    "\n",
    "In order to update the top of the hierarchy, we need the total probability of the data for each hypothetical zoo.  When we do an update using grid algorithms, we get the probability of the data free, since it is the normalizing constant.\n",
    "\n",
    "But when we do an update using a conjugate distribution, we don't get the total probability of the data, and for a Dirichlet distribution it is not easy to compute.\n",
    "\n",
    "However, we can estimate it by drawing samples from the Dirichlet distribution, and then computing the probability of the data for each sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Dirichlet(object):\n",
    "    \"\"\"Represents a Dirichlet distribution.\n",
    "\n",
    "    See http://en.wikipedia.org/wiki/Dirichlet_distribution\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, n, conc=1, label=None):\n",
    "        \"\"\"Initializes a Dirichlet distribution.\n",
    "\n",
    "        n: number of dimensions\n",
    "        conc: concentration parameter (smaller yields more concentration)\n",
    "        label: string label\n",
    "        \"\"\"\n",
    "        if n < 2:\n",
    "            raise ValueError('A Dirichlet distribution with '\n",
    "                             'n<2 makes no sense')\n",
    "\n",
    "        self.n = n\n",
    "        self.params = np.ones(n, dtype=np.float) * conc\n",
    "        self.label = label if label is not None else '_nolegend_'\n",
    "\n",
    "    def update(self, data):\n",
    "        \"\"\"Updates a Dirichlet distribution.\n",
    "\n",
    "        data: sequence of observations, in order corresponding to params\n",
    "        \"\"\"\n",
    "        m = len(data)\n",
    "        self.params[:m] += data\n",
    "\n",
    "    def random(self):\n",
    "        \"\"\"Generates a random variate from this distribution.\n",
    "\n",
    "        Returns: normalized vector of fractions\n",
    "        \"\"\"\n",
    "        p = np.random.gamma(self.params)\n",
    "        return p / p.sum()\n",
    "\n",
    "    def mean(self):\n",
    "        \"\"\"Array of means.\"\"\"\n",
    "        return self.params / self.params.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an example that represents a zoo with 4 animals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<__main__.Dirichlet at 0x7f7ca0ca5320>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d4 = Dirichlet(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a sample from it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.22243685, 0.13145446, 0.59255112, 0.05355757])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = d4.random()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compute the probability of the data, given these prevalences, using the multinomial distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.006761575978171484"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import multinomial\n",
    "\n",
    "data = [3, 2, 1, 0]\n",
    "m = sum(data)\n",
    "multinomial(m, p).pmf(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since I only observed 3 species, and my hypothetical zoo has 4, I had to zero-pad the data.  Here's a function that makes that easier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def zero_pad(a, n):\n",
    "    \"\"\"Why does np.pad have to be so complicated?\n",
    "    \"\"\"\n",
    "    res = np.zeros(n)\n",
    "    res[:len(a)] = a\n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3., 2., 1., 0.])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = [3, 2, 1]\n",
    "zero_pad(data, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's pull all that together.  Here's a function that estimates the total probability of the data by sampling from the dirichlet distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_likelihood(dirichlet, data, iters=1000):\n",
    "    \"\"\"Estimate the total probability of the data.\n",
    "    \n",
    "    dirichlet: Dirichlet object\n",
    "    data: array of observed frequencies\n",
    "    iters: number of samples to draw\n",
    "    \"\"\"\n",
    "    data = zero_pad(data, dirichlet.n)\n",
    "    m = np.sum(data)\n",
    "    likes = [multinomial(m, dirichlet.random()).pmf(data) \n",
    "             for i in range(iters)]\n",
    "    return np.mean(likes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.012789548125351505"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_likelihood(d4, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready to solve the problem.\n",
    "\n",
    "Here's a Suite that represents the set of possible zoos.  The likelihood of any zoo is just the total probability of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Zoo(Suite):\n",
    "    \n",
    "    def Likelihood(self, data, hypo):\n",
    "        \"\"\"\n",
    "        data: sequence of counts\n",
    "        hypo: Dirichlet object\n",
    "        \"\"\"\n",
    "        return sample_likelihood(hypo, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can construct the prior by enumerating the possible zoos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "suite = Zoo([Dirichlet(n, label=''.join(zoo))\n",
    "             for n, zoo in enumerate_zoos(species, present)]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.125 LTB\n",
      "0.125 LTBC\n",
      "0.125 LTBR\n",
      "0.125 LTBE\n",
      "0.125 LTBCR\n",
      "0.125 LTBCE\n",
      "0.125 LTBRE\n",
      "0.125 LTBCRE\n"
     ]
    }
   ],
   "source": [
    "def print_zoos(suite):\n",
    "    for d, p in suite.Items():\n",
    "        print(p, d.label)\n",
    "        \n",
    "print_zoos(suite)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can update the top level of the hierarchy by calling `Update`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.010737491081704907"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "suite.Update(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have to update the bottom level explicitly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "for hypo in suite:\n",
    "    hypo.update(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the posterior for the top level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.39926034376784497 LTB\n",
      "0.1364576975358686 LTBC\n",
      "0.1365266963559979 LTBR\n",
      "0.1369789254977338 LTBE\n",
      "0.057347250996821975 LTBCR\n",
      "0.06027653795703476 LTBCE\n",
      "0.04856388605939697 LTBRE\n",
      "0.02458866182930092 LTBCRE\n"
     ]
    }
   ],
   "source": [
    "print_zoos(suite)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we can get the posterior distribution of `n`, the number of species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_n = Pmf()\n",
    "for d, p in suite.Items():\n",
    "    pmf_n[d.n] += p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.81610465490401\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAG55JREFUeJzt3X+0XWV95/H3xyDaij+oZGpNAokaOwVtwUlBxyl1WsSobcJMcQwdLXbsUDtE7bK/oGWhUmmrndpqoVVaaKmKkWJtb20s0iqdoQ6ai1I1YIYQUG5jazQgWigY/M4fZ4c5HM/9Re7OfZK8X2vdlb2f/Tx7f++Gez9377PPc1JVSJLUmkcsdgGSJI1jQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJnSRfT/KUnvb9vCRTQ+tbkzxvgfb9X5N8eGi9kjxtIfbd7a+38zLDMb8tyV8m+WqSP92fx1Y7DCjtd0luT3Jv94vvn5P8UZIj9mF/K7tfyoftS11VdURV7diXfczjWMdV1bUz9Znr91VV76mqUxeiriTXJvmpkf3vt/My5HTgO4EnVtVL9vOx1QgDSovlR6vqCOBZwPcD5y1WIfsabPs6/kA9ds+OAf5vVe1Z7EK0eAwoLaqq+kfgQ8AzAJI8OclEkt1Jtif573v7JjkxyWSSu7srr7d2m/5X9+9d3VXZc7r+/y3JzUnuTHJ1kmOG9lVJzk5yC3DLUNvTuuXHJ/mTJLuSfD7JeUke0W17RZK/T/LbSXYDbxj9vrpbVH/cHfsmBiE8vP32JKfM9/sad+yu7bqREl6UZEeSLyf5zaHa35Dk3UN1PHiVluRC4AeAi7rjXfQwzst1Sf5n933fluSF0/23T/I93RXbXd0tz3Vd+xuB84GXdnW8cszYNyS5sqvla934NdMdSwemg/WvLx0gkqwAXgT8Wdf0XmAr8GTg3wLXJNlRVX8LvA14W1W9q7sl+IxuzMnAbcAT9v7FneQ04JeBH2UQQOd0+/73Q4c/DTgJuHdMab8LPB54CvBE4MPAF4FLu+0nAZuAfwM8csz41wNP7b4ewyCEpzOf7+u7xxz7pWP2+Z+ANcARwN8A24A/nKEGqupXkjwXeHdVTdd3LuflcuAo4Czg0iTLamROtSSPBP4SuAw4FfgPwF8kWVNVr09SwNOq6mUzlLwO+M/ATwJvAi4Cnj3T96gDi1dQWix/nuQu4Drg74Bf68LqPwC/VFX/WlU3Mvil+vJuzDeApyU5qqq+XlXXz7D/nwZ+vapu7n65/xpw/PBVVLd9d1U9JKCSLGHwS//cqvpaVd0O/NZQHQA7q+p3q2rP6PjOfwEu7PZ/B/D2GWqdz/c1l2MDvLk79heA3wHOmGWfs5rjefl8Vf1BVT3AIKi+i8FrSaOezSA8f6Oq7q+qjwAfnGed11XV5u5Y7wK+b97flJpmQGmxnFZVT6iqY6rqf3S/aJ8M7K6qrw31+zywrFt+JfB04HNJtiT5kRn2fwzwtu720V3AbiBD+wK4Y5qxRwGHd8ceV8dMY/d68kifz0/Xkfl9X3M59mifz3f17Ku5nJd/2rtQVfd0i+MegHkycEdVfXOGfc3mn4aW7wEefRC/JndIMqDUkp3AdyR57FDb0cA/AlTVLVV1BoNbW28GrkryGGDclPx3AD/dheDer2+rqo8N9ZluKv8vM7iqGb7aerCOWcbu9UVgxcj4seb5fc3l2Iw59s5u+V+Abx/a9qR57Hsu52WudgIr9r5+tY/70kHKgFIzulthHwN+Pcmjk3wvg6uL9wAkeVmSpd1f3Xd1wx4AdgHfZPC6yF7vAM5Nclw39vFJ5vS4cnfL6ErgwiSP7W4Lvg5498wjH+LK7vhHJlkOvHq6jvP8vubqF7pjrwBeC7yva78RODnJ0UkeD5w7Mu6fpzveAp2XvT7OICx/MckjM3hP2I8yeG1NAgwotecMYCWDv7A/ALy+qq7ptq0Ftib5OoMHCzZ0r1XdA1wI/H13S+/ZVfUBBlcjm5LcDXwWmPaJsjFezeAX6A4Gr5NdweAF/bl6I4NbVrcxeJDgXTP0nfP3NY/j/wVwA4NA+iu6hxi6c/k+4NPd9g+OjHsbcHr3FN6418329bzQ1XE/g4ccXsjgyuz3gJ+oqs/Nd186eMUPLJQktcgrKElSkwwoSVKTDChJUpMMKElSkw6aN7UdddRRtXLlysUuQ5I0ixtuuOHLVbV0tn4HTUCtXLmSycnJxS5DkjSLJDPNrPIgb/FJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKadNDMJKE2bXzTexe7hAPKReedsdglSM3wCkqS1KReAyrJ2iTbkmxPcs4M/U5PUknWDLWd243bluQFfdYpSWpPb7f4kiwBLgaeD0wBW5JMVNVNI/0eC7wG+PhQ27HABuA44MnA3yR5elU90Fe9kqS29HkFdSKwvap2VNX9wCZg/Zh+vwq8BfjXobb1wKaquq+qbgO2d/uTJB0i+nxIYhlwx9D6FHDScIckJwArquqDSX5+ZOz1I2OXjR4gyVnAWQBHH330PhfsC/rz4wv6kvrU5xVUxrTVgxuTRwC/DfzcfMc+2FB1SVWtqao1S5fO+tlXkqQDSJ9XUFPAiqH15cDOofXHAs8Ark0C8CRgIsm6OYyVJB3k+ryC2gKsTrIqyeEMHnqY2Luxqr5aVUdV1cqqWsnglt66qprs+m1I8qgkq4DVwCd6rFWS1JjerqCqak+SjcDVwBLgsqramuQCYLKqJmYYuzXJlcBNwB7gbJ/gk6RDS68zSVTVZmDzSNv50/R93sj6hcCFvRUnSWqaM0lIkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkprUa0AlWZtkW5LtSc4Zs/1VST6T5MYk1yU5tmtfmeTerv3GJO/os05JUnt6+0TdJEuAi4HnA1PAliQTVXXTULcrquodXf91wFuBtd22W6vq+L7qkyS1rc8rqBOB7VW1o6ruBzYB64c7VNXdQ6uPAarHeiRJB5A+A2oZcMfQ+lTX9hBJzk5yK/AW4DVDm1Yl+VSSv0vyA+MOkOSsJJNJJnft2rWQtUuSFlmfAZUxbd9yhVRVF1fVU4FfAs7rmr8IHF1VJwCvA65I8rgxYy+pqjVVtWbp0qULWLokabH1GVBTwIqh9eXAzhn6bwJOA6iq+6rqK93yDcCtwNN7qlOS1KA+A2oLsDrJqiSHAxuAieEOSVYPrb4YuKVrX9o9ZEGSpwCrgR091ipJakxvT/FV1Z4kG4GrgSXAZVW1NckFwGRVTQAbk5wCfAO4EzizG34ycEGSPcADwKuqandftUqS2tNbQAFU1WZg80jb+UPLr51m3PuB9/dZmySpbc4kIUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWpSrwGVZG2SbUm2JzlnzPZXJflMkhuTXJfk2KFt53bjtiV5QZ91SpLa01tAJVkCXAy8EDgWOGM4gDpXVNUzq+p44C3AW7uxxwIbgOOAtcDvdfuTJB0i+ryCOhHYXlU7qup+YBOwfrhDVd09tPoYoLrl9cCmqrqvqm4Dtnf7kyQdIg7rcd/LgDuG1qeAk0Y7JTkbeB1wOPBDQ2OvHxm7bMzYs4CzAI4++ugFKVqS1IY+r6Aypq2+paHq4qp6KvBLwHnzHHtJVa2pqjVLly7dp2IlSW3pM6CmgBVD68uBnTP03wSc9jDHSpIOMn0G1BZgdZJVSQ5n8NDDxHCHJKuHVl8M3NItTwAbkjwqySpgNfCJHmuVJDWmt9egqmpPko3A1cAS4LKq2prkAmCyqiaAjUlOAb4B3Amc2Y3dmuRK4CZgD3B2VT3QV62SpPb0+ZAEVbUZ2DzSdv7Q8mtnGHshcGF/1UmSWuZMEpKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJvUaUEnWJtmWZHuSc8Zsf12Sm5J8OsnfJjlmaNsDSW7sviZGx0qSDm69faJukiXAxcDzgSlgS5KJqrppqNungDVVdU+SnwHeAry023ZvVR3fV32SpLb1eQV1IrC9qnZU1f3AJmD9cIeq+mhV3dOtXg8s77EeSdIBpM+AWgbcMbQ+1bVN55XAh4bWH51kMsn1SU4bNyDJWV2fyV27du17xZKkZvR2iw/ImLYa2zF5GbAG+MGh5qOrameSpwAfSfKZqrr1ITurugS4BGDNmjVj9y1JOjD1eQU1BawYWl8O7BztlOQU4FeAdVV13972qtrZ/bsDuBY4ocdaJUmN6TOgtgCrk6xKcjiwAXjI03hJTgDeySCcvjTUfmSSR3XLRwHPBYYfrpAkHeR6u8VXVXuSbASuBpYAl1XV1iQXAJNVNQH8JnAE8KdJAL5QVeuA7wHemeSbDEL0N0ae/pMkHeT6fA2KqtoMbB5pO39o+ZRpxn0MeGaftUmS2uZMEpKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmzRhQSf54aPnM3quRJKkz2xXU9w0tv7bPQiRJGjZbQDm/nSRpUcz2Rt3lSd7OYOLXvcsPqqrX9FaZJOmQNltA/cLQ8mSfhUiSNGzGgKqqy/dXIZIkDZsxoJJMzLS9m9hVkqQFN9stvucw+FTc9wIfZ/yHEEqStOBmC6gnAc8HzgB+HPgr4L1VtbXvwiRJh7YZHzOvqgeq6q+r6kzg2cB24Nokr94v1UmSDlmzfh5U98m2L2ZwFbUSeDvwZ/2WJUk61M32kMTlwDOADwFvrKrP7peqJEmHvNlmkng58HQG0xz9nyR3d19fS3L3bDtPsjbJtiTbk5wzZvvrktyU5NNJ/jbJMUPbzkxyS/flPICSdIiZ7X1QD3u28yRLgIsZPGQxBWxJMlFVNw11+xSwpqruSfIzwFuAlyb5DuD1wBoG0y3d0I298+HWI0k6sMw2m/mjk/xskouSnJVk1teshpwIbK+qHVV1P7AJWD/coao+WlX3dKvXA8u75RcA11TV7i6UrgHWzuPYkqQD3GxXSJczuIr5DPAi4Lfmse9lDN5DtddU1zadVzJ4rWvOY7vQnEwyuWvXrnmUJklq3WxXRMdW1TMBklwKfGIe+x73pt6xs6MneRmDIPzB+YytqkuASwDWrFnjzOuSdBCZ7QrqG3sXqmrPPPc9BawYWl8O7BztlOQU4FeAdVV133zGSpIOXrN+YOHwk3vA987jKb4twOokq5IcDmwAHjK3X5ITgHcyCKcvDW26Gjg1yZFJjgRO7dokSYeI2Z7iW/Jwd1xVe5JsZBAsS4DLqmprkguAyaqaAH4TOAL40yQAX6iqdVW1O8mvMgg5gAuqavfDrUWSdOCZz1N581ZVm4HNI23nDy2fMsPYy4DL+qtOktSyh/0+J0mS+mRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKa1OtcfJIW3sY3vXexSzigXHTeGYtdgh4mr6AkSU0yoCRJTTKgJElNMqAkSU3qNaCSrE2yLcn2JOeM2X5ykk8m2ZPk9JFtDyS5sfuaGB0rSTq49fYUX5IlwMXA84EpYEuSiaq6aajbF4BXAD8/Zhf3VtXxfdUnSWpbn4+Znwhsr6odAEk2AeuBBwOqqm7vtn2zxzokSQegPm/xLQPuGFqf6trm6tFJJpNcn+S0cR2SnNX1mdy1a9e+1CpJakyfAZUxbTWP8UdX1Rrgx4HfSfLUb9lZ1SVVtaaq1ixduvTh1ilJalCfATUFrBhaXw7snOvgqtrZ/bsDuBY4YSGLkyS1rc+A2gKsTrIqyeHABmBOT+MlOTLJo7rlo4DnMvTalSTp4NdbQFXVHmAjcDVwM3BlVW1NckGSdQBJvj/JFPAS4J1JtnbDvweYTPIPwEeB3xh5+k+SdJDrdbLYqtoMbB5pO39oeQuDW3+j4z4GPLPP2iRJbXMmCUlSkwwoSVKTDChJUpMMKElSkwwoSVKTDChJUpMMKElSkwwoSVKTDChJUpMMKElSkwwoSVKTDChJUpMMKElSkwwoSVKTDChJUpMMKElSkwwoSVKTeg2oJGuTbEuyPck5Y7afnOSTSfYkOX1k25lJbum+zuyzTklSe3oLqCRLgIuBFwLHAmckOXak2xeAVwBXjIz9DuD1wEnAicDrkxzZV62SpPb0eQV1IrC9qnZU1f3AJmD9cIequr2qPg18c2TsC4Brqmp3Vd0JXAOs7bFWSVJj+gyoZcAdQ+tTXduCjU1yVpLJJJO7du162IVKktrTZ0BlTFst5NiquqSq1lTVmqVLl86rOElS2/oMqClgxdD6cmDnfhgrSToI9BlQW4DVSVYlORzYAEzMcezVwKlJjuwejji1a5MkHSJ6C6iq2gNsZBAsNwNXVtXWJBckWQeQ5PuTTAEvAd6ZZGs3djfwqwxCbgtwQdcmSTpEHNbnzqtqM7B5pO38oeUtDG7fjRt7GXBZn/VJktrlTBKSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCb1GlBJ1ibZlmR7knPGbH9Ukvd12z+eZGXXvjLJvUlu7L7e0WedkqT29PaJukmWABcDzwemgC1JJqrqpqFurwTurKqnJdkAvBl4abft1qo6vq/6JElt6/MK6kRge1XtqKr7gU3A+pE+64HLu+WrgB9Okh5rkiQdIPoMqGXAHUPrU13b2D5VtQf4KvDEbtuqJJ9K8ndJfmDcAZKclWQyyeSuXbsWtnpJ0qLqM6DGXQnVHPt8ETi6qk4AXgdckeRx39Kx6pKqWlNVa5YuXbrPBUuS2tFnQE0BK4bWlwM7p+uT5DDg8cDuqrqvqr4CUFU3ALcCT++xVklSY/oMqC3A6iSrkhwObAAmRvpMAGd2y6cDH6mqSrK0e8iCJE8BVgM7eqxVktSY3p7iq6o9STYCVwNLgMuqamuSC4DJqpoALgXelWQ7sJtBiAGcDFyQZA/wAPCqqtrdV62SpPb0FlAAVbUZ2DzSdv7Q8r8CLxkz7v3A+/usTZLUNmeSkCQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1qdf3QUlS6za+6b2LXcIB56Lzztgvx/EKSpLUJANKktQkA0qS1CQDSpLUJANKktQkA0qS1CQDSpLUJANKktQkA0qS1KReAyrJ2iTbkmxPcs6Y7Y9K8r5u+8eTrBzadm7Xvi3JC/qsU5LUnt4CKskS4GLghcCxwBlJjh3p9krgzqp6GvDbwJu7sccCG4DjgLXA73X7kyQdIvq8gjoR2F5VO6rqfmATsH6kz3rg8m75KuCHk6Rr31RV91XVbcD2bn+SpENEqqqfHSenA2ur6qe69ZcDJ1XVxqE+n+36THXrtwInAW8Arq+qd3ftlwIfqqqrRo5xFnBWt/rdwLZ5lnkU8OV5jlkMB0qdcODUap0L70Cp1ToX3nxrPaaqls7Wqc/ZzDOmbTQNp+szl7FU1SXAJfMvrTt4MllVax7u+P3lQKkTDpxarXPhHSi1WufC66vWPm/xTQErhtaXAzun65PkMODxwO45jpUkHcT6DKgtwOokq5IczuChh4mRPhPAmd3y6cBHanDPcQLY0D3ltwpYDXyix1olSY3p7RZfVe1JshG4GlgCXFZVW5NcAExW1QRwKfCuJNsZXDlt6MZuTXIlcBOwBzi7qh7oocyHfXtwPztQ6oQDp1brXHgHSq3WufB6qbW3hyQkSdoXziQhSWqSASVJatJBHVBJViT5aJKbk2xN8toxfZ6X5KtJbuy+zl+kWh+d5BNJ/qGr9Y1j+kw7NVRjdb4iya6hc/pT+7vOoVqWJPlUkg+O2bbo53Oknplqbemc3p7kM10dk2O2J8nbu/P66STParTOVn72n5DkqiSf635XPWdkexPnc461Lug57fN9UC3YA/xcVX0yyWOBG5JcU1U3jfT731X1I4tQ37D7gB+qqq8neSRwXZIPVdX1Q30enBoqyQYGU0O9tME6Ad43/KbsRfRa4GbgcWO2tXA+h81UK7RzTgH+Y1VN98bMFzJ48nY1gzfe/37372KYqU5o42f/bcBfV9Xp3RPP3z6yvaXzOVutsIDn9KC+gqqqL1bVJ7vlrzH44V+2uFWNVwNf71Yf2X2NPsEy3dRQ+80c62xCkuXAi4E/nKbLop/PveZQ64FkPfAn3f8r1wNPSPJdi11Ui5I8DjiZwRPNVNX9VXXXSLcmzucca11QB3VADetu35wAfHzM5ud0t6w+lOS4/VrYkO4Wz43Al4Brqmq01mXAHTB4jB/4KvDE/VvlnOoE+LHudsRVSVaM2b4//A7wi8A3p9nexPnszFYrtHFOYfAHyYeT3JDBdGOjHjyvnSkW5w/D2eqExf/ZfwqwC/ij7vbuHyZ5zEifVs7nXGqFBTynh0RAJTkCeD/ws1V198jmTzKYF+r7gN8F/nx/17dXVT1QVcczmDnjxCTPGOkypymg+jaHOv8SWFlV3wv8Df//KmW/SfIjwJeq6oaZuo1p2+/nc461Lvo5HfLcqnoWg1tPZyc5eWR7E+eV2ets4Wf/MOBZwO9X1QnAvwCjH03UyvmcS60Lek4P+oDqXid5P/Ceqvqz0e1VdffeW1ZVtRl4ZJKj9nOZozXdBVzL4KNGhk03NdSimK7OqvpKVd3Xrf4B8O/2c2kAzwXWJbmdwUz6P5Tk3SN9Wjmfs9bayDndW8vO7t8vAR/gWz9poImpymars5Gf/SlgauguxFUMQmC0z6KfT+ZQ60Kf04M6oLrXEy4Fbq6qt07T50l7X3dIciKDc/KV/Vflg3UsTfKEbvnbgFOAz410m25qqP1mLnWO3B9fx+C1v/2qqs6tquVVtZLBDCUfqaqXjXRb9PMJc6u1hXPa1fGY7oEjuts7pwKfHek2AfxE9/TZs4GvVtUXW6uzhZ/9qvon4I4k3901/TCDGXSGLfr5hLnVutDn9GB/iu+5wMuBz3SvmQD8MnA0QFW9g8Evpp9Jsge4F9iwGL+kgO8CLs/ggxkfAVxZVR/MHKaGarDO1yRZx+Apyt3AKxahzrEaPJ/TavScfifwge530GHAFVX110leBQ/+TG0GXsTgc9zuAX6y0Tpb+dl/NfCe7qm4HcBPNng+95qt1gU9p051JElq0kF9i0+SdOAyoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoKQGJVmZweft/EEGn7v14W7mDumQYUBJ7VoNXFxVxwF3AT+2yPVI+5UBJbXrtqraO0XXDcDKRaxF2u8MKKld9w0tP8DBP3em9BAGlCSpSQaUJKlJzmYuSWqSV1CSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCb9PxX783vNs6OkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "thinkplot.Hist(pmf_n)\n",
    "print(pmf_n.Mean())\n",
    "thinkplot.decorate(xlabel='n', \n",
    "                   ylabel='PMF', \n",
    "                   title='Posterior distribution of n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, to answer the question, we have to compute the posterior distribution of the prevalence of elephants.  Here's a function that computes it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def enumerate_posterior(suite):\n",
    "    for d, p in suite.Items():\n",
    "        mean = d.mean()\n",
    "        index = d.label.find('E')\n",
    "        p_elephant = 0 if index == -1 else mean[index]\n",
    "        yield d, p, p_elephant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the possible zoos, the posterior probability of each, and the conditional prevalence of elephants for each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LTB 0.39926034376784497 0\n",
      "LTBC 0.1364576975358686 0\n",
      "LTBR 0.1365266963559979 0\n",
      "LTBE 0.1369789254977338 0.1\n",
      "LTBCR 0.057347250996821975 0\n",
      "LTBCE 0.06027653795703476 0.09090909090909091\n",
      "LTBRE 0.04856388605939697 0.09090909090909091\n",
      "LTBCRE 0.02458866182930092 0.08333333333333333\n"
     ]
    }
   ],
   "source": [
    "for d, p, p_elephant in enumerate_posterior(suite):\n",
    "    print(d.label, p, p_elephant)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can use the law of total probability to compute the probability of seeing an elephant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02564153170370892"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total = np.sum(p * p_elephant \n",
    "               for d, p, p_elephant in enumerate_posterior(suite))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}