{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Think Bayes\n",
    "\n",
    "This notebook presents example code and exercise solutions for 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 classes from thinkbayes2\n",
    "from thinkbayes2 import Hist, Pmf, Suite\n",
    "import thinkplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Grizzly Bear Problem\n",
    "\n",
    "In 1996 and 1997 Mowat and Strobeck deployed bear traps in locations in British Columbia and Alberta, in an effort to estimate the population of grizzly bears.  They describe the experiment in \"Estimating Population Size of Grizzly Bears Using Hair Capture, DNA Profiling, and Mark-Recapture Analysis\"\n",
    "\n",
    "The \"trap\" consists of  a lure and several strands of barbed wire intended to capture samples of hair from bears that visit the lure.  Using the hair samples, the researchers use DNA analysis to identify individual bears.\n",
    "\n",
    "During the first session, on June 29, 1996, the researchers deployed traps at 76 sites.  Returning 10 days later, they obtained 1043 hair samples and identified 23 different bears.  During a second 10-day session they obtained 1191 samples from 19 different bears, where 4 of the 19 were from bears they had identified in the first batch.\n",
    "\n",
    "To estimate the population of bears from this data, we need a model for the probability that each bear will be observed during each session.  As a starting place, we'll make the simplest assumption, that every bear in the population has the same (unknown) probability of being sampled during each round.\n",
    "\n",
    "We also need a prior distribution for the population.  As a starting place, let's suppose that, prior to this study, an expert in this domain would have estimated that the population is between 100 and 500, and equally likely to be any value in that range.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution:\n",
    "\n",
    "Define:\n",
    "\n",
    "* N: population size\n",
    "\n",
    "* K: number of bears that have ever been identified\n",
    "\n",
    "* n: number of bears observed in the second second\n",
    "\n",
    "* k: the number of bears in the second session that had previously been identified\n",
    "\n",
    "\n",
    "For given values of N, K, and n, the distribution of k is the hypergeometric distribution:\n",
    "\n",
    "$PMF(k) = {K \\choose k}{N-K \\choose n-k}/{N \\choose n}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "from scipy.special import binom\n",
    "\n",
    "class Grizzly(Suite):\n",
    "    \"\"\"Represents hypotheses about how many bears there are.\"\"\"\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        \"\"\"Computes the likelihood of the data under the hypothesis.\n",
    "\n",
    "        hypo: total population (N)\n",
    "        data: # tagged (K), # caught (n), # of caught who were tagged (k)\n",
    "        \"\"\"\n",
    "        N = hypo\n",
    "        K, n, k = data\n",
    "\n",
    "        if hypo < K + (n - k):\n",
    "            return 0\n",
    "\n",
    "        like = binom(N-K, n-k) / binom(N, n)\n",
    "        return like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.05801258299152e-06"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "hypos = range(100, 501)\n",
    "suite = Grizzly(hypos)\n",
    "\n",
    "data = 23, 19, 4\n",
    "suite.Update(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt8VdWZ//HPQ+4h4R65a7ipBFTAiChVqyigVVBLW2yd6tT5OW21rbUzrba/6a/jvDpT27HaTnVaq52qbUVEHVOr4gWslyoQQJSLSOQikVuAQAgJCUme3x97cziEXElOzknyfb9eeWXvtdfe59mbJA9rr73XMndHRETkRPWIdwAiItK5KZGIiEibKJGIiEibKJGIiEibKJGIiEibKJGIiEibKJGIiEibKJGIiEibKJGIiEibJMc7gI4wYMAAz83NjXcYIiKdxvLly3e7e05L6naLRJKbm0thYWG8wxAR6TTMbEtL6+rWloiItIkSiYiItIkSiYiItIkSiYiItIkSiYiItIkSiYiItEm3ePy3Pbk7xTv3UXmomsyMVIae1Aczi3dYIiJxo0TSQvsPVPLkwuW8sXwD5RVVkfJeWRmce2Yu1142iZP6ZccxQhGR+FAiaYG3393I/Y+/RuWh6uO2lZVX8vLf1rFoyXquuXQin59xNklJumMoIt2H/uI1Y+Gba7jnf146Jon0yspg1PAcemVlRMpqa+tYsHA5/+9Xf2b/gcp4hCoiEhdqkTRh6fubefDJNyLrJ/XL5qY5n+LsvJMxM9yd1Ru28cfnlrBhyy4A1m3czr/88ll++PUrGdA3K16hi4h0GLVIGvHJrn384rFXI+ujhufwk9uvJX/cKZHOdTPjjFOH8u+3Xc11n5mMRe37w/8qYN+BijhELiLSsZRIGlB5qJqfPrSQQ1WHAcjpm82/fO0z9M7OaLB+jx49mDN9Erf//WWR/pGde8r4t/9+norK4/tVRES6EiWSBrjD4JzeAKQkJ/G9f5hBds/0Zvc7f8Io/unvp0daJps/2c39f1qMu8cwWhGR+FIiaUBmRirf+4cZzL3iHL429yJGDBvQ4n0nn5HL1667KLL+znubKFj8XizCFBFJCDFNJGY208zWm1mRmd3RwPY0M3si3L7EzHKjtt0Zlq83sxlR5d82szVmttrMHjez5psKJxY7n5txNhedc2qr9502ZSyXXzA+sv6HgndYU7StPcMTEUkYMUskZpYE3A9cDuQB15lZXr1qNwGl7j4auBe4O9w3D5gLjANmAg+YWZKZDQW+CeS7+3ggKayXcG68+jzGnHISAHXu/OKxVzlYWdXMXiIinU8sWySTgSJ33+ju1cA8YHa9OrOBR8LlBcA0Cx6Jmg3Mc/cqd98EFIXHg+CR5QwzSwYygYT8r35ychL/9PfTycpMA2DPvoP8/pm34xyViEj7i2UiGQpsjVovDssarOPuNcB+oH9j+7r7J8B/Ah8D24H97v5STKJvBwP6ZnHz5y+MrC9a8gGFa1o8e6WISKcQy0TS0EiG9R9faqxOg+Vm1pegtTICGAL0NLPrG/xws5vNrNDMCktKSloRdvuaOnEU508cFVn/9by/6haXiHQpsUwkxcDwqPVhHH8bKlInvFXVG9jbxL6XApvcvcTdDwNPA+c39OHu/qC757t7fk5OTjuczom7+XMXRN5BKS2rYP4Ly+Maj4hIe4plIlkGjDGzEWaWStApXlCvTgFwQ7g8B1jkwUsXBcDc8KmuEcAYYCnBLa0pZpYZ9qVMA9bF8BzaRXbPdL5yzdTI+vOvv8+WbXviGJGISPuJWSIJ+zxuBRYS/LGf7+5rzOwuM5sVVnsY6G9mRcDtwB3hvmuA+cBa4EXgFnevdfclBJ3yK4D3w/gfjNU5tKepk0YxfswQIHiK6zfz39CLiiLSJVh3+GOWn5/vhYWF8Q6DrTtK+c5Pn6S2tg6AW794MRefe1qcoxIROZ6ZLXf3/JbU1ZvtHWj4oL7M+vSZkfVHC95Rx7uIdHpKJB1szoyz6d+nJxBMivXMyyvjHJGISNsokXSw9LQUvjzrvMj6n//6PiV7D8QxIhGRtlEiiYOpk0YxanjwSHJNTS1/+svSOEckInLilEjiwMy44eqjrZLXCzfw0cfxe2lSRKQtlEjiZNzoIUw+Izey/sizb+txYBHplJRI4uj6WVPo0SP4J1hTtI2V67Y2s4eISOJRIomjoSf14bLzxkbWH39+mVolItLpKJHE2ZwZk0hJTgJg49YSlry3Kc4RiYi0jhJJnPXr3fOY2RTnPb+Murq6OEYkItI6SiQJ4JpLJ5CWmgIEw6i8uaIozhGJiLScEkkC6JWVwVUXHx065YkXCqmpqY1jRCIiLadEkiBmXXwmPTOCaXl37C5j8dL1cY5IRKRllEgSRM+MNGZPOyuy/uTC5VQfroljRCIiLaNEkkA+c+EZ9MoKZlLcs+8gL/8t4efsEhGJbSIxs5lmtt7Miszsjga2p5nZE+H2JWaWG7XtzrB8vZnNCMtOM7N3o77KzOy2WJ5DR0pPS+Gzl02MrD/98kq1SkQk4cUskZhZEnA/cDmQB1xnZnn1qt0ElLr7aOBe4O5w3zyCqXnHATOBB8wsyd3Xu/sEd58AnA1UAM/E6hziYfrUPPr2ygRg34EKtUpEJOHFskUyGShy943uXg3MA2bXqzMbeCRcXgBMC+dinw3Mc/cqd98EFIXHizYN+Mjdt8TsDOIgNSWZay492ip55hW1SkQkscUykQwFogePKg7LGqwTzvG+H+jfwn3nAo+3Y7wJ47Lzx9InO2iVlJZV8MrbapWISOKKZSKxBsrqDyTVWJ0m9zWzVGAW8GSjH252s5kVmllhSUnnGqI9aJVMiKw/88q7HD6s90pEJDHFMpEUA8Oj1ocB2xqrY2bJQG9gbwv2vRxY4e47G/twd3/Q3fPdPT8nJ+eETyJepk/Ni7RK9u4/yKvvfBDniEREGhbLRLIMGGNmI8IWxFygoF6dAuCGcHkOsMiD4W8LgLnhU10jgDFA9DSC19FFb2sdkZqSzNXTjrZKnn5lhVolIpKQYpZIwj6PW4GFwDpgvruvMbO7zGxWWO1hoL+ZFQG3A3eE+64B5gNrgReBW9y9FsDMMoHLgKdjFXuimD51LL2zj75XsmiJWiUiknisO8x/kZ+f74WFhfEO44Q8u2gVjz77NgD9+/TkgX/5IsnhsPMiIrFiZsvdPb8ldfVme4KbMTXvmLfdFy3RGFwikliUSBJceloKsy85OgbXUy+v0MjAIpJQlEg6gZmfGkd2z3QAdpeWa2RgEUkoSiSdwHGtkpdWqlUiIglDiaSTuPyC8ZFWSUnpAf5a+GGcIxIRCSiRdBLpaSnMuvhoq2TBQvWViEhiUCLpRC6/YBxZmcEsirv2qlUiIolBiaQTyUhPZdYlapWISGJRIulkrrhgvFolIpJQlEg6GbVKRCTRKJF0QmqViEgiUSLphDLSU5l9ydGRgZ98Ua0SEYkfJZJOKvoJLr1XIiLxpETSSalVIiKJQomkE7viwvHHtEpeW6ZWiYh0PCWSTiwYg+toq0RPcIlIPMQ0kZjZTDNbb2ZFZnZHA9vTzOyJcPsSM8uN2nZnWL7ezGZElfcxswVm9oGZrTOz82J5DonuiguPHYNLIwOLSEeLWSIxsyTgfuByIA+4zszy6lW7CSh199HAvcDd4b55BHO8jwNmAg+ExwP4BfCiu58OnEUwjW+3pZGBRSTeYtkimQwUuftGd68G5gGz69WZDTwSLi8AppmZheXz3L3K3TcBRcBkM+sFXEgw1zvuXu3u+2J4Dp1C/ZGBNYuiiHSkWCaSocDWqPXisKzBOu5eA+wH+jex70igBPgfM1tpZg+ZWc+GPtzMbjazQjMrLCkpaY/zSVjpaSlcPe1oX4lmURSRjhTLRGINlHkL6zRWngxMAv7b3ScCB4Hj+l4A3P1Bd8939/ycnJyWR91J1Z9FUa0SEekosUwkxcDwqPVhwLbG6phZMtAb2NvEvsVAsbsvCcsXECSWbk+tEhGJl1gmkmXAGDMbYWapBJ3nBfXqFAA3hMtzgEXu7mH53PCprhHAGGCpu+8AtprZaeE+04C1MTyHTmXmp8bRKysDCFolr77zQZwjEpHuIGaJJOzzuBVYSPBk1Xx3X2Nmd5nZrLDaw0B/MysCbie8TeXua4D5BEniReAWdz/y3+tvAH80s/eACcC/x+ocOpv6rZIFL62g+nBNHCMSke7AggZA15afn++FhYXxDqNDHKo6zC3/9jj7DlQAcMPV5x0zRa+ISEuY2XJ3z29JXb3Z3sWkp6UwZ8bRbqOnXlpBRWV1HCMSka5OiaQLuuy8seT0zQagvKKKP7/2XpwjEpGuTImkC0pOTmLuFUdbpAWLV1FWXhnHiESkK1Mi6aIuzB/DsIF9gaDf5JlX3o1zRCLSVSmRdFE9evTgus+cE1l//o3V7C4tj2NEItJVKZF0YeeeOYJRw4O3+mtqalnw0vI4RyQiXZESSRdmZnzpqnMj66++s57tJfvjGJGIdEVKJF3cmacOZfyYIQDU1dXx+PPL4hyRiHQ1SiRdnJnxxc9Mjqy/taKIjVu79mjIItKxlEi6gdNGDOKc8bmR9UeefZvuMKKBiHQMJZJu4vpZ59LDgtH5V2/Yxoq1H8c5IhHpKpRIuolhA/ty6fljI+uPFbxDbW1dHCMSka5CiaQb+cLl+aSlpgCwdUcpi5dq8isRaTslkm6kT3YmV087OhLwvOeXcajqcBwjEpGuQImkm5l18Vn07ZUJQGlZBQWLV8U5IhHp7GKaSMxsppmtN7MiMztubvVwBsQnwu1LzCw3atudYfl6M5sRVb7ZzN43s3fNrHtMMtKO0tNSjhk65X9fXUVpWUUcIxKRzi5micTMkoD7gcuBPOA6M8urV+0moNTdRwP3AneH++YRTM07DpgJPBAe74iL3X1CSyddkWNdPPk0hg/uB0BV9WEe/8vSOEckIp1ZLFskk4Eid9/o7tXAPGB2vTqzgUfC5QXANDOzsHyeu1e5+yagKDyetIMePXrw5VlTIuuL3vlALymKyAmLZSIZCmyNWi8OyxqsE87xvh/o38y+DrxkZsvN7ObGPtzMbjazQjMrLCnRH8n6JuWdzKS8k4Hggj701Ft6SVFETkgsE4k1UFb/L1VjdZrad6q7TyK4ZXaLmV3Y0Ie7+4Punu/u+Tk5OS2NuVu58Zrz6dEj+BFYv2kHb634KM4RiUhnFMtEUgwMj1ofBmxrrI6ZJQO9gb1N7evuR77vAp5Bt7xO2NCT+nDlRWdE1h8teFuPA4tIq8UykSwDxpjZCDNLJeg8L6hXpwC4IVyeAyzy4P5KATA3fKprBDAGWGpmPc0sG8DMegLTgdUxPIcub86MSfTKygBgz76DPPOqZlIUkdaJWSIJ+zxuBRYC64D57r7GzO4ys1lhtYeB/mZWBNwO3BHuuwaYD6wFXgRucfdaYCDwppmtApYCf3H3F2N1Dt1Bz4w0rr/qaKPu2VffZdfeA3GMSEQ6G+sOHaz5+fleWKhXThrj7nz3nqcjT25NOWsk//yV6XGOSkTiycyWt/QVC73ZLpgZN107NbL+zqqNrFpfHMeIRKQzUSIRAE4fOYgL88dE1n/75BtUH66JY0Qi0lkokUjEl2efR2Z6KgDbS/bzzCvqeBeR5jWZSMzs91HLNzRRVbqAvr0y+dKV50bWn35lJdt27YtjRCLSGTTXIjkravlbsQxEEsP0qWMZNTx4gbOmppbfPvmm3ngXkSY1l0j0F6Sb6dGjB1/9woWRoQXe+7CYt1bqjXcRaVxyM9uHmdkvCYYsObIc4e7fjFlkEjcjh+dw+YXjef714F3P3z/zNyaOHU7PjLQ4RyYiiai5Fsk/A8uBwqjl6C/pouZecc4xE2D96TkNNS8iDWuyReLujzS1Xbqunhlp3HjN+dz7yCsAvPjmGqZOGk3eqMFxjkxEEk2TicTM6o+NdQx3n9XUduncpk4cxevLNrB87RYAHnj8NX7+vc+RmtLcHVER6U6a+4twHsG8II8DS2h4eHfposyMmz9/Abf9ZDuVh6rZXrKfec8v48uzz4t3aCKSQJrrIxkEfB8YD/wCuAzY7e5/dfe/xjo4ib8BfbO4YfbR2RQLFq2iaMuuOEYkIommyUTi7rXu/qK73wBMIZjy9jUz+0aHRCcJ4dLzxnLGqcEElQ786vHXqKmpjW9QIpIwmh0iJZwT5FrgD8AtwC+Bp2MdmCQOM+Nrcy+K9I1s3b6XBS+viHNUIpIomhsi5RHgb8Ak4F/d/Rx3/zd3/6RDopOEMbB/L66/6ujwKU+9tJJNxbvjGJGIJIrmWiR/B5xKMDzK22ZWFn4dMLOy5g5uZjPNbL2ZFZnZHQ1sTzOzJ8LtS8wsN2rbnWH5ejObUW+/JDNbaWbPteQkpX1cceF4ThsxCIC6ujrue/RVjRAsIs32kfRw9+yor17hV7a792pqXzNLAu4HLgfygOvMLK9etZuAUncfDdwL3B3um0cwNe84YCbwQHi8I75FMOuidCAz49Yvfjpyi6t4Zyl/+POSOEclIvHW3K2tdDO7zcx+ZWY3m1lrXiCYDBS5+0Z3rwbmAbPr1ZkNHHnpcQEwzcwsLJ/n7lXuvomgk39yGNMw4DPAQ62IRdrJkJP68PfXnB9Z/8tf39ckWCLdXHO3th4B8oH3gSuAe1px7KEE76AcURyWNVgnnON9P9C/mX3vA74L1DX14WHiKzSzwpKSklaELc257PyxnJ13SmT9v/6wiAMHD8UxIhGJp+YSSZ67X+/uvwHmABe04tgNvbxYfzThxuo0WG5mVwK73L3Zcb7c/UF3z3f3/JycnOajlRYzM77+xYvolZUBBGNx/Wb+GxpuXqSbai6RHD6yELYYWqMYGB61PgzY1lid8LZZb2BvE/tOBWaZ2WaCW2WXmNkfWhmXtIM+2Zl8be5FkfW33/2I1ws3xDEiEYmXZie2in5SCzizFU9tLQPGmNkIM0sl6DyvP3ZXAXBk5sU5wCIP/ltbAMwNn+oaAYwBlrr7ne4+zN1zw+MtcvfrW3y20q4mn5HLtCmnR9Z/M/8Nzago0g0199RWUr0ntZJb+tRW2IK5FVhI8ITVfHdfY2Z3mdmRwR4fBvqbWRFwO3BHuO8aYD6wFngRuMXd9Sp1AvrKtVMZNCD4UaiqPsw9v39FjwSLdDPWHe5r5+fne2FhYbzD6LI2Fe/mez9/mtra4PmHGVPHcfPnW9OdJiKJxsyWu3t+S+o2O0SKSHNGDBvAV66ZGllf+NYaTc8r0o0okUi7mPGpPKacNTKy/sDjr7G9ZH8cIxKRjqJEIu3CzPj6dRcxsH/QX3Ko6jD3/P5l9ZeIdANKJNJuemak8Z0bLyMpKfix2lS8m98++abeLxHp4pRIpF2NOjmHG68+OoPioiUfsPDNtXGMSERiTYlE2t3lF4znwvwxkfWHn36LtR9tj2NEIhJLSiTS7o5MhDVi2AAgGHL+Z797id2l5XGOTERiQYlEYiI1JZnv3TQjMh5XWXklP314oTrfRbogJRKJmZx+2XznxkvpYcEYnB9tLdHgjiJdkBKJxNT4MUO5MWr+kteWruepl1fGMSIRaW9KJBJzV1w4novPPS2y/vhflvLm8qI4RiQi7UmJRGLOzPjq5y9k/JghkbL/+tNiPti4I45RiUh7USKRDpGcnMR3b5rBsIF9AaipqeU/fvuChlER6QKUSKTD9MxI4/v/eHnkSa7yiip+/JvnNU2vSCenRCIdamD/Xtz5f2aSkpwEwPaS/fz4N89zqOpwM3uKSKKKaSIxs5lmtt7Miszsjga2p5nZE+H2JWaWG7XtzrB8vZnNCMvSzWypma0yszVm9q+xjF9i49TcgXzr76Zh4fqGLbu4+6GFHD6suctEOqOYJRIzSwLuBy4H8oDrzCyvXrWbgFJ3Hw3cC9wd7ptHMJXuOGAm8EB4vCrgEnc/C5gAzDSzKbE6B4md8yaM5B/mHJ386r0Pi7nvsVepq6uLY1QiciJi2SKZDBS5+0Z3rwbmAbPr1ZkNPBIuLwCmmZmF5fPcvcrdNwFFwGQPHBlnIyX80tttndTMC8Yx94pzIuvvrNrIr594XS8sinQysUwkQ4GtUevFYVmDdcI53vcD/Zva18ySzOxdYBfwsrsviUn00iHmTJ/ElRedGVl/9Z0PeKzgHSUTkU4klonEGiir/9ehsTqN7uvute4+ARgGTDaz8Q1+uNnNZlZoZoUlJSWtCFs6kplx4zXn8enJR19YfHbRKh7/yzIlE5FOIpaJpBgYHrU+DNjWWB0zSwZ6A3tbsq+77wNeI+hDOY67P+ju+e6en5OTc+JnITFnZnx97kWcMz43UvbUyyuY97ySiUhnEMtEsgwYY2YjzCyVoPO8oF6dAuCGcHkOsMiDvxwFwNzwqa4RwBhgqZnlmFkfADPLAC4FPojhOUgHSUrqwXduvIyz806JlC14aQXzXihUMhFJcDFLJGGfx63AQmAdMN/d15jZXWY2K6z2MNDfzIqA24E7wn3XAPOBtcCLwC3uXgsMBhab2XsEiepld38uVucgHSslJYl//sp0JuWdHClbsHA5814ojGNUItIc6w7/28vPz/fCQv0x6iyqD9fw04cXsnLd0ectrr10Il+8cjJmDXWfiUh7M7Pl7p7fkrp6s10STmpKMt+9aQYTxx7tJnv6lZX89sk3dZtLJAEpkUhCOpJMovtMFr61hvsee5WaGr0BL5JIlEgkYQXJZDpTJ42OlL25vIifPvySpuwVSSBKJJLQkpOTuO3vLmH61KOj6yxfu4W7HvgL5RVVcYxMRI5QIpGE16NHD27+3AV89rJJkbJ1G7fz/XufYeeesjhGJiKgRCKdhJnxxSsn8+XZ50XKPtm1jzt+/gwbtuyMY2QiokQincrsS87i2zdcSnI4n0lZeSX/8ssC3lm1Mc6RiXRfSiTS6Xxq0mh+9PUrycpMA+BwTS3/+buXeHbRKj0eLBIHSiTSKY0dNZj/+PY1DBrQCwhG9Hz02be577FXqarWbIsiHUmJRDqtISf14T++fQ2njRgUKXtzeRHfv+9Zdu09EMfIRLoXJRLp1HplZXDXrVdx2fljI2WbP9nNd//zKVZv+CSOkYl0H0ok0uklJyfx1S9cxD9+/kKSkoIf6QMHD/Gv9z/H0y+vVL+JSIwpkUiXMX1qHnfdOove2RkA1Lnzx+eW8OPfPM/+A5Vxjk6k61IikS7l9JGD+Nk/ffaYfpOV67bynZ8+ydqPtscxMpGuS4lEupz+fbK469aruGbahEhZaVkFP/zls8x/sZDa2ro4RifS9cQ0kZjZTDNbb2ZFZnZHA9vTzOyJcPsSM8uN2nZnWL7ezGaEZcPNbLGZrTOzNWb2rVjGL51XcnIS18+awvdvvjzyvokDT7xQyA9+8b9s27UvvgGKdCExSyRmlgTcD1wO5AHXmVlevWo3AaXuPhq4F7g73DePYGrecQRzsj8QHq8G+I67jwWmALc0cEyRiLPHncI93/0cY0cOjpRt2LKL7/x0AS+9tVYd8SLtIJYtkslAkbtvdPdqYB4wu16d2cAj4fICYJoFU+DNBua5e5W7bwKKgMnuvt3dVwC4+wGCKXyHxvAcpAsY0DeLu75xFV+68tzIU13Vh2v4zfzX+fcHX2B3aXmcIxTp3GKZSIYCW6PWizn+j36kTjjH+36gf0v2DW+DTQSWtGPM0kX16NGDay+byN23X8uwgX0j5SvWfsy3/uMJXnxjjVonIicolomkocm16/+mNlanyX3NLAt4CrjN3RscR9zMbjazQjMrLCkpaWHI0tWNGDaAn/3zZ7nq02dGyg5VHea3C97gB794luKdpXGMTqRzimUiKQaGR60PA7Y1VsfMkoHewN6m9jWzFIIk8kd3f7qxD3f3B909393zc3Jy2ngq0pWkpiRz4zXn82/fnM2QnN6R8vWbdnD73U8y/8VCTecr0gqxTCTLgDFmNsLMUgk6zwvq1SkAbgiX5wCLPLi/UADMDZ/qGgGMAZaG/ScPA+vc/ecxjF26gbxRg7nne5/js5dNokeP4FehtraOJ14o5Ns/mc/KdVubOYKIAFgs7wub2RXAfUAS8Dt3/7GZ3QUUunuBmaUDjxH0dewF5rr7xnDfHwBfIXhS6zZ3f8HMPgW8AbwPHHkZ4Pvu/nxTceTn53thYWEMzlC6ii3b9nD/n17jo63H3gY9Z3wuN15zfmSUYZHuwsyWu3t+i+p2hw5GJRJpibq6Op5/fTXzXiik8lB1pDw5OYmrp03g2ksnkJaaEscIRTqOEkk9SiTSGqVlFfzhz0t4ben6Y8r79e7J3CvyuXjyaZFbYSJdlRJJPUokciI+3LyThxa8edztrmED+3L9rHPJH3cKQbedSNejRFKPEomcKHdn0ZIP+ONzS48bQfj0kYP4u6umcPrIQY3sLdJ5KZHUo0QibXWo6jAFi1fxv6+uOm4q34ljh/P5mfmcmjswTtGJtD8lknqUSKS97D9QyYKXlrPwrbXHjSI84fThfH7m2ccMYS/SWSmR1KNEIu1te8l+5r2wjLeWFx03XMOZpw7js9MnMm70EPWhSKelRFKPEonESvHOUhYsXMGbyzccl1BGDc9h1iVncd5ZIyODRYp0Fkok9SiRSKx9smsfCxYu543C4xNKTt9srrr4TKZNOZ30NL2HIp2DEkk9SiTSUbbt2kfB4lUsXvrhceN19cxI49LzTmf61HF6U14SnhJJPUok0tH2H6jk+TdW8+IbqymvqDpmmwETxg5n5gXjmTR2uF5ulISkRFKPEonES1X1YRYv+ZA/v7aKHbuPn/Egp28206fmccmU0+iTnRmHCEUapkRSjxKJxFtdXR0r1m1l4ZtrWLn24+P6UXqYMSnvZC4+9zTyx51CcnJSXOIUOUKJpB4lEkkkO3aX8fLf1vLK2+uOu+0FkN0znQvzx3DJuaeRO3RAHCIUUSI5jhKJJKLqwzX8beVHvPz2Oj7YuKPBOicP7sfUSaOZOnEUg6Mm4RKJNSWSepRIJNFtL9nP4iXreW3ZevbsO9hgnVHDc5g6aTTnTxhJTr/sDo5QupuESSRmNhPryj84AAARkklEQVT4BcHEVg+5+0/qbU8DHgXOBvYAX3D3zeG2O4GbgFrgm+6+MCz/HXAlsMvdx7ckDiUS6Szq6up4f8M2Xn3nA5a+t4nDjUz5e9qIQUw5awTnjM9VS0ViIiESiZklAR8ClxHMwb4MuM7d10bV+Tpwprt/1czmAte4+xfMLA94HJgMDAFeAU5191ozuxAoBx5VIpGurPJQNYWrt/DmiiJWfrD1uLG9jhg+qC+TzxjB5DNyGXVyjoZlkXbRmkSSHMM4JgNFUVPnzgNmA2uj6swGfhQuLwB+Fc7LPhuY5+5VwCYzKwqP97a7v25muTGMWyQhZKSnckH+GC7IH0N5RRVL3tvIWys+4v0PP6Eu6j+AW3eUsnVHKU+9vIK+vTI554xcJo49mTPGDCEjPTWOZyDdRSwTyVBga9R6MXBuY3XcvcbM9gP9w/J36u07NHahiiS2rMw0pk0Zy7QpY9l/oJJlqzez9L3NrPqw+Jg36EvLKnjprbW89NZakpJ6cFruQCaOPZmJY4eTO7S/WisSE7FMJA39xNa/j9ZYnZbs2/SHm90M3Axw8sknt2ZXkYTWOzuDS88by6XnjeVQ1WHe/WArS9/fzPI1W455nLi2to61H21n7Ufb+eNzS+iVlcGE04dx1mnDGD9mKAP6ZsXxLKQriWUiKQaGR60PA7Y1UqfYzJKB3sDeFu7bJHd/EHgQgj6SVkUu0kmkp6Uw5ayRTDlrJLW1dXywaQcr1n7MynVb2bJtzzF1y8oreb1wA68XbgBgYP9e5I0ezPjRQxg3eoieBJMTFsvO9mSCzvZpwCcEne1fdPc1UXVuAc6I6my/1t0/b2bjgD9xtLP9VWCMu9eG++UCz6mzXaRxe/cf5L31xaxYt5X31hdz4OChJuuf1C+bvNFDGDtyEKeNGMSwgX10K6wbS4intsJArgDuI3j893fu/mMzuwsodPcCM0sHHgMmErRE5kZ1zv8A+ApQA9zm7i+E5Y8DnwYGADuB/+fuDzcVhxKJdHfuzkcfl/Du+mLWbNjGuo3bG320+IjM9FROzR3IqbkDOW3EQMacchI9M9I6KGKJt4RJJIlCiUTkWDU1tRR9XMLqom2s2bCNDzbtoPpwTZP7GDBsUF9Gn3ISo4bnMHLYAHKH9ictVXOsdEVKJPUokYg0raamlo+2lrCmaDsfbt7J+s07KSuvbHa/I8llxLABQXIZnsOIof312HEXoERSjxKJSOu4Ozv3HODDzTtYvylILFs+2XPM+yuNMWBQTm9OGdyPYYP7cfLgfpwypB+DB/TWlMOdiBJJPUokIm13qOowH20tYePW3WwsDr5/srO0xc/lJyX1YNjAvpwypB/DB/Xj5CH9GD6oLzl9szS5VwJKlDfbRaQLSU9LYVz4qPARh6oOs/mTPUGCKd7Nxq0lFO8obbDlUltbx5Zte457LDk5OYnBOb0ZmtObISf1Cb+C5eye6TE/L2k7JRIROWHpaSmcPnIQp48cFCmrPlxD8Y5SPt6+ly3b9rJ1x14+3r630VGNa2pq2bp9L1u37z1uW1ZmGkMH9mVwTm8GDejFoP69OKl/NgMH9KJ3VoYeT04QSiQi0q5SU5IZGXa8RyuvqGLr9r3HJJjinfua7NQvr6hi/aYdrN90/HwtqSnJDOyfzUn9ejFwwJHvvcKybHX4dyAlEhHpEFmZaYwdNZixowYfU36wsoptu/axbdd+tu3axye79vPJzlK2l+xv8l2X6sM1kQErG5KRnsqAPj0Z0DeL/n2y6N+nJwP6ZNG/bxYD+mYxoE9PPbrcTpRIRCSuemakMeaUgYw5ZeAx5e7O7tJytpUECWbXngPs2lPGjvB7xaHqJo9beaiarTuqG000ECS3fr2DZNMnO5O+vTLp2zszstynVyb9emeSmqI/lU3R1RGRhGRm5PTLJqdfNmedNuy47eUVVezcXcbOvWXs2nOAnXuOfi8pLT9mVOTGlFdUUV5RxccN9M9Ey0xPpU92Bn16ZdK3d0/6ZmfSp1cGfbIzyc5Kp3dWOtk9M+idlU56Wkq367tRIhGRTikrM42sk3MYdXLOcdvcnbLyQ+zZV87ufQfZs6+cPaXllOwrZ0/pQXaXlrO37GCjk4XVV3GomopD1Wwr2d9s3eTkJHr1TKdXVga9eqZHEs2R9V5ZGfTKSie7Z3pwDplpnb7F07mjFxFpgJnROzuD3tkZx3X6H+Hu7DtQyZ7ScvbsP8i+sgr2llWw/0AFpfsrKC2rYN+BCkrLKqmra1nCgeAptL37D7J3f8NPqTUkJTmJrMw0emak0TMzjayMNHpmpoaJ5mjCObotXM9IJSU5Ke4tICUSEemWzCzoE+mVyegm6rk7Bw4eorSskn0HKthXFiSZ0v0VlB2spKz8EPvLKzlw8BD7D1Q2OxhmQw7X1AbHLKto9b5JST3ITE+lZ0YqGempZB75ygjKRgwdwCVTTm/1cVtDiUREpAlmFt6OyuAU+jVb/1DVYcoOHqLsQCVlBw8FCaa8MrJeVn4oKC+vpLyymoOVVS2+xdaQ2to6DoSf05DJZ+QqkYiIdCbpaSmkp6VwUgsnCnN3qqprKK+o4mBlVfi9mvKDhyivrOJg+EBA9PLBiqoWJ6HMDhj6X4lERCSOzCySfE5k+uPqwzVUHKrmYGU1lZXVR5fD70MH9olB1MeKaSIxs5nALwgmtnrI3X9Sb3sa8ChwNrAH+IK7bw633QncBNQC33T3hS05pohId5KakkxqSjJ9sjPjFkPMhtw0syTgfuByIA+4zszy6lW7CSh199HAvcDd4b55wFxgHDATeMDMklp4TBER6UCxHLt5MlDk7hvdvRqYB8yuV2c28Ei4vACYZsFzbLOBee5e5e6bgKLweC05poiIdKBYJpKhwNao9eKwrME67l4D7Af6N7FvS44pIiIdKJaJpKE3ZOpPUtBYndaWH//hZjebWaGZFZaUlDQZqIiInLhYJpJiYHjU+jBgW2N1zCwZ6A3sbWLflhwTAHd/0N3z3T0/J6fhN1tFRKTtYplIlgFjzGyEmaUSdJ4X1KtTANwQLs8BFnkw928BMNfM0sxsBDAGWNrCY4qISAeK2eO/7l5jZrcCCwke1f2du68xs7uAQncvAB4GHjOzIoKWyNxw3zVmNh9YC9QAt7h7LUBDx4zVOYiISPPMG5hbuasxsxJgywnuPgDY3Y7htBfF1TqKq3UUV+t0xbhOcfcW9Qt0i0TSFmZW6O758Y6jPsXVOoqrdRRX63T3uGLZRyIiIt2AEomIiLSJEknzHox3AI1QXK2juFpHcbVOt45LfSQiItImapGIiEibdOtEYma/M7NdZrY6qqyfmb1sZhvC733DcjOzX5pZkZm9Z2aTOjiuH5nZJ2b2bvh1RdS2O8O41pvZjBjGNdzMFpvZOjNbY2bfCsvjes2aiCuu18zM0s1sqZmtCuP617B8hJktCa/XE+HLtYQv4D4RxrXEzHI7OK7fm9mmqOs1ISzvsJ/98POSzGylmT0Xrsf1ejURV6Jcr81m9n4YQ2FY1rG/k+7ebb+AC4FJwOqosp8Cd4TLdwB3h8tXAC8QjPc1BVjSwXH9CPinBurmAauANGAE8BGQFKO4BgOTwuVs4MPw8+N6zZqIK67XLDzvrHA5BVgSXof5wNyw/NfA18LlrwO/DpfnAk/E6Ho1FtfvgTkN1O+wn/3w824H/gQ8F67H9Xo1EVeiXK/NwIB6ZR36O9mtWyTu/jrBG/XRooe2fwS4Oqr8UQ+8A/Qxs8EdGFdjGhtyPxZxbXf3FeHyAWAdwejLcb1mTcTVmA65ZuF5l4erKeGXA5cQTJsAx1+vhqZV6Ki4GtNhP/tmNgz4DPBQuG7E+Xo1FFczOux6NRNDh/1OdutE0oiB7r4dgj9QwElheSIMYX9r2Bz93ZGmarziCm8jTCT432zCXLN6cUGcr1l4O+RdYBfwMkHrZ58H0ybU/+zGplWIeVzufuR6/Ti8XvdaMIPpMXE1EHN7uw/4LnBkIvL+JMD1aiCuI+J9vSD4T8BLZrbczG4Oyzr0d1KJpOVaPIR9jPw3MAqYAGwH7gnLOzwuM8sCngJuc/eypqo2UBaz2BqIK+7XzN1r3X0CwUjVk4GxTXx23OIys/HAncDpwDlAP+B7HRmXmV0J7HL35dHFTXx2POOCOF+vKFPdfRLBzLG3mNmFTdSNSWxKJMfbeaSpF37fFZa3eAj7WHD3neEvfx3wW47eiunQuMwsheCP9R/d/emwOO7XrKG4EuWahbHsA14juC/dx4JpE+p/dmPTKnREXDPDW4Tu7lXA/9Dx12sqMMvMNhPMfnoJQUsg3tfruLjM7A8JcL0AcPdt4fddwDNhHB36O6lEcrzooe1vAJ6NKv9y+NTDFGD/kaZjR6h3H/Ma4MgTXY0NuR+LGIxgxOZ17v7zqE1xvWaNxRXva2ZmOWbWJ1zOAC4l6L9ZTDBtAhx/vRqaVqEj4vog6g+PEdxTj75eMf93dPc73X2Yu+cSdJ4vcvcvEefr1Uhc18f7eoWf3dPMso8sA9PDODr2d7I9euw76xfwOMEtj8MEmfomgnusrwIbwu/9wroG3E9wj/t9IL+D43os/Nz3wh+GwVH1fxDGtR64PIZxfYqgGfwe8G74dUW8r1kTccX1mgFnAivDz18N/DAsH0mQuIqAJ4G0sDw9XC8Kt4/s4LgWhddrNfAHjj7Z1WE/+1ExfpqjT0fF9Xo1EVfcr1d4bVaFX2uAH4TlHfo7qTfbRUSkTXRrS0RE2kSJRERE2kSJRERE2kSJRERE2kSJRERE2kSJRLosM3Mzuydq/Z/M7EftdOzfm9mc5mu2+XM+Z8GoxovrlX/awlFoReJNiUS6sirgWjMbEO9AoplZUiuq3wR83d0vjlU80aLeIBdpMSUS6cpqCKYa/Xb9DfVbFGZWHn7/tJn91czmm9mHZvYTM/uSBfN3vG9mo6IOc6mZvRHWuzLcP8nMfmZmy8LB/P4x6riLzexPBC+C1Y/nuvD4q83s7rDshwQvW/7azH7WwPn1MrNnzGytmf3azHqE+003s7fNbIWZPRmOQYaZ/TCMa7WZPRi+kY2ZvWZm/25mfwW+FbaCVlswX8nrrb/s0t3ofx/S1d0PvGdmP23FPmcRDK64F9gIPOTuky2YMOsbwG1hvVzgIoKBIReb2WjgywTDTpxjwWiwb5nZS2H9ycB4D4aujzCzIcDdwNlAKcFIrle7+11mdgnBnCqFDcQ5mWBulS3AiwStr9eA/wtc6u4Hzex7BPNo3AX8yt3vCj/zMeBK4M/hsfq4+0XhtveBGe7+yZGhVESaohaJdGkejAL8KPDNVuy2zIMB+aoIhpI4kgjeJ0geR8x39zp330CQcE4nGOvoyxYM0b6EYKiKMWH9pfWTSOgc4DV3L/FgOPQ/Ekxu1pyl7r7R3WsJhtX5FMGgkHkECexdgnGWTgnrX2zBTILvEwyIOC7qWE9ELb8F/N7M/g/Qmttw0k2pRSLdwX3ACoIRWo+oIfyPVHiLJzVqW1XUcl3Ueh3H/s7UH1/ICcYy+oa7L4zeYGafBg42Et+JTsbU2Oe/7O7X1fv8dOABgrGVtoYPHaRHVYnE5u5fNbNzCSZyetfMJrj7nhOMUboBtUiky3P3vQTTtd4UVbyZ4FYSBLPGpZzAoT9nZj3CfpORBANALgS+ZsGw9pjZqeGorE1ZAlxkZgPCjvjrgL+24PMnWzCfeQ/gC8CbwDvA1PA2G2aWaWancjRp7A77TBp94szMRrn7Enf/IbCbY4cdFzmOWiTSXdwD3Bq1/lvgWTNbSjA6amOthaasJ/iDPxD4qrsfMrOHCG5/rQhbOiUcnea0Qe6+3czuJBgu3YDn3f3ZpvYJvQ38BDgDeB14xt3rzOxG4HE7OmPf/3X3D83stwS35zYDy5o47s/MbEwYy6sEI8uKNEqj/4qISJvo1paIiLSJEomIiLSJEomIiLSJEomIiLSJEomIiLSJEomIiLSJEomIiLSJEomIiLTJ/we+p+PUHbGlBgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "thinkplot.Pdf(suite)\n",
    "thinkplot.Config(xlabel='Number of bears', ylabel='PMF', legend=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior mean 193.93845891363907\n",
      "Maximum a posteriori estimate 109\n",
      "90% credible interval (105, 379)\n"
     ]
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "print('Posterior mean', suite.Mean())\n",
    "print('Maximum a posteriori estimate', suite.MaximumLikelihood())\n",
    "print('90% credible interval', suite.CredibleInterval(90))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# Alternatively, we can take advantage of the `hypergeom`\n",
    "# object in scipy.stats.\n",
    "\n",
    "from scipy import stats\n",
    "\n",
    "class Grizzly2(Suite):\n",
    "    \"\"\"Represents hypotheses about how many bears there are.\"\"\"\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        \"\"\"Computes the likelihood of the data under the hypothesis.\n",
    "\n",
    "        hypo: total population (N)\n",
    "        data: # tagged (K), # caught (n), # of caught who were tagged (k)\n",
    "        \"\"\"\n",
    "        N = hypo\n",
    "        K, n, k = data\n",
    "\n",
    "        if hypo < K + (n - k):\n",
    "            return 0\n",
    "\n",
    "        like = stats.hypergeom.pmf(k, N, K, n)\n",
    "        return like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.07135370142238903"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "hypos = range(100, 501)\n",
    "suite = Grizzly2(hypos)\n",
    "\n",
    "data = 23, 19, 4\n",
    "suite.Update(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior mean 193.9384589136376\n",
      "Maximum a posteriori estimate 109\n",
      "90% credible interval (105, 379)\n"
     ]
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "print('Posterior mean', suite.Mean())\n",
    "print('Maximum a posteriori estimate', suite.MaximumLikelihood())\n",
    "print('90% credible interval', suite.CredibleInterval(90))"
   ]
  },
  {
   "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": 2
}