{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Rock Hyrax Problem\n",
    "----------------------\n",
    "\n",
    "Allen Downey\n",
    "\n",
    "This notebook contains a solution to a problem I posed in my Bayesian statistics class:\n",
    "\n",
    "> Suppose I capture and tag 10 rock hyraxes.  Some time later, I capture another 10 hyraxes and find that two of them are\n",
    "> already tagged.  How many hyraxes are there in this environment?\n",
    "\n",
    "This is an example of a mark and recapture experiment, which you can [read about on Wikipedia](http://en.wikipedia.org/wiki/Mark_and_recapture).  The Wikipedia page also includes the photo of a tagged hyrax shown above.\n",
    "\n",
    "As always with problems like this, we have to make some modeling assumptions.\n",
    "\n",
    "1) For simplicity, you can assume that the environment is reasonably isolated, so the number of hyraxes does not change between observations.\n",
    "\n",
    "2) And you can assume that each hyrax is equally likely to be captured during each phase of the experiment, regardless of whether it has been tagged.  In reality, it is possible that tagged animals would avoid traps in the future, or possible that the same behavior that got them caught the first time makes them more likely to be caught again.  But let's start simple.\n",
    "\n",
    "My solution uses the ThinkBayes2 framework, which is described in [Think Bayes](http://thinkbayes.com), and summarized in [this notebook](http://nbviewer.ipython.org/github/AllenDowney/ThinkBayes2/blob/master/code/framework.ipynb).\n",
    "\n",
    "I'll start by defining terms:\n",
    "\n",
    "$N$: total population of hyraxes\n",
    "$K$: number of hyraxes tagged in the first round\n",
    "$n$: number of hyraxes caught in the second round\n",
    "$k$: number of hyraxes in the second round that had been tagged\n",
    "\n",
    "So $N$ is the hypothesis and $(K, n, k)$ make up the data.  The probability of the data, given the hypothesis, is the probability of finding $k$ tagged hyraxes out of $n$ if (in the population) $K$ out of $N$ are tagged.  There are two ways we can compute this:\n",
    "\n",
    "1) If you are familiar with the [hypergeometric distribution](http://en.wikipedia.org/wiki/Hypergeometric_distribution), you might recognize this problem and use an implementation of the hypergeometric PMF, evaluated at $k$.\n",
    "\n",
    "2) Otherwise, you can figure it out using combinatorics.\n",
    "\n",
    "I'll do the second one first.  Out of a population of $N$ hyraxes, we captured $n$; the total number of combinations is $N \\choose n$.\n",
    "\n",
    "$k$ of the ones we caught are tagged, so $n-k$ are not.  The total number of combinations is ${K \\choose k}{N-K \\choose n-k}$.  So the probability of the data is\n",
    "\n",
    "${K \\choose k}{N-K \\choose n-k}/{N \\choose n}$\n",
    "\n",
    "`scipy.special` provides `binom(x, y)`, which computes the binomial coefficient, $x \\choose y$.\n",
    "\n",
    "So let's see how that looks in code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# first a little house-keeping\n",
    "from __future__ import print_function, division\n",
    "% matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import thinkbayes2\n",
    "from scipy.special import binom\n",
    "\n",
    "class Hyrax(thinkbayes2.Suite):\n",
    "    \"\"\"Represents hypotheses about how many hyraxes 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": "markdown",
   "metadata": {},
   "source": [
    "Again $N$ is the hypothesis and $(K, n, k)$ is the data.  If we've tagged $K$ hyraxes and then caught another $n-k$, the total number of unique hyraxes we're seen is $K + (n - k)$.  For any smaller value of N, the likelihood is 0.\n",
    "\n",
    "Notice that I didn't bother to compute $K \\choose k$; because it does not depend on $N$, it's the same for all hypotheses, so it gets cancelled out when we normalize the suite.\n",
    "\n",
    "Next I construct the prior and update it with the data.  I use a uniform prior from 0 to 999."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0010248876709531896"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hypos = range(1, 1000)\n",
    "suite = Hyrax(hypos)\n",
    "\n",
    "data = 10, 10, 2\n",
    "suite.Update(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the posterior distribution looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEPCAYAAACOU4kjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XHW9//HXZCYzWdok3Zc0bZou0NKFspRy2cKi1qpF\nfqLAxYsXvRdcql5X9HK9ptfrRb2iCChWWfUKeHHBohRkC+gtlgKlLXRN2nRJ6ULatGmWmSST3x/f\nMzNnpsnMJJnJzOS8n4/HPPI9Z77fM98c6HzyXQ+IiIiIiIiIiIiIiIiIiIiIiIiIiIiIZJ2lwDZg\nJ3BLH3nutN7fCCxKouxi4BVgA7AeODe1VRYRkWzmBuqASiAfeAOYE5NnGfCklT4P+FsSZWuB91jp\n9wIvpLriIiIyOHlpvPZiTIBoADqBR4ErY/IsBx6y0uuAMmBigrJvA6VWugxoTEflRURk4DxpvHY5\nsM92vB/TOkmUpxyYHKfs14C/At/HBMfzU1dlERFJhXS2XHqSzOfq53XvAz4HTAW+ANzfz/IiIpJm\n6Wy5NAIVtuMKTAskXp4pVp78OGUXA1dY6d8A9/b24TNmzOipr68fUMVFRBysHpiZ6UrE48FUshLw\nknhAfwmRAf14ZV8HLrHSl2NmjPWmR4xvfvObma5C1tC9iNC9iNC9iCD5Xqe40tly6QJWAE9jZn/d\nB2wFbrbeX4UJLMswg/etwI0JygLcBPwY8AHt1rGIiGSRdAYXgDXWy25VzPGKfpQFeJVTJwaIiEgW\nSeeAvmSJ6urqTFcha+heROheROhepF5/Z2rlEqv7UEREkuVyuSAFsUEtFxERSTkFFxERSTkFFxER\nSTkFFxERSTkFFxERSTkFFxERSTkFFxERSTkFFxERSTkFFxERSTkFFxERSTkFFxERSbl074qcczo6\nAjz8u7/QsO8I73/3OSxelNXPzBERyUrauDLGt25/jKee3wCAKy+Pu/7rEyyaPz3VdRMRyUrauDIN\nGg8eDQcWgJ5gkHv/59kM1khEJDcpuNj85W9bTzn3xpu72VF/IAO1ERHJXQouNpu37On1fO3at4a4\nJiIiuU3BxdLT08PmrZHg8vG/vzycfunlLZmokohIzkp3cFkKbAN2Arf0kedO6/2NwKIkyj4KbLBe\nu62fg9Z0rIWmoy0AFBX5uPaqC/Dkm8l0u/ccoulYSyo+RkTEEdIZXNzA3ZggMRe4DpgTk2cZMBOY\nBdwE3JNE2WsxQWgR8FvrNWh7978TTldWjKe4qIAzZk8Jn3tj8+5UfIyIiCOkM7gsBuqABqAT0+K4\nMibPcuAhK70OKAMmJlnWBXwEeCQVld3/dlM4PWXyGADOtE1BfuOthlR8jIiII6QzuJQD+2zH+61z\nyeSZnETZi4BDQH0qKrv/gC24TLKCyzxbcHmzIRUfIyLiCOkMLsmuYBzoYp3rgIcHWPYUexsj3WIV\n5Sa4zJ8zlTy3uUW79hyita0jVR8nIjKspXP7l0agwnZcgWmBxMszxcqTn6CsB7gKOCteBWpqasLp\n6upqqqur+8x7+MjxcHrShFEAFBZ4mT51PPW7D0JPD9vrDnDWgqp4HykiklNqa2upra1N+XXTGVxe\nxQzUVwIHgGswrQ271cAKzJjKEqAZ09XVlKDsFcBW670+2YNLIu8cPRFOjxtTGk7PmTXFBBdge72C\ni4gML7F/eK9cuTIl101ncOnCBI6nMbO/7sMEhJut91cBT2JmjNUBrcCNCcqGXEOKBvIBuruDHG1u\nDR+PGTUinJ49Y3I4vW1nY6o+UkRkWEv3rshrrJfdqpjjFf0oG3JjH+cH5GjzSXqCQQDKSovJz4/c\nljmzIvMIttUpuIiIJEMr9CG8eBJg7JiSqPdmVE7E7XEDsL/xHU62alBfRCQRBReix1vGjBoZ9Z7P\nl8+Mygnh4+1qvYiIJKTgArxjb7mMHnnK+7OrIuMudQ0Hh6ROIiK5TMEFOH4iMphfVlp8yvtVtpbL\nroZDQ1InEZFcpuACHG9pD6dLRxad8v6MyonhdN1utVxERBJRcAGOH4+0XEpLegku0yItl917D9Pd\nHRySeomI5CoFF+B4S1s43VtwGVU2gjHWWIzfH4ja5FJERE6l4AKcsHWLlfTSLQZQZesa271H4y4i\nIvEouBDTcukjuMzUuIuISNIUXIATJyLBpaSXbjEgaq2LpiOLiMTn+ODS3R3kxMn4s8UgesZYvaYj\ni4jE5fjgcrK1A3rMo2eKiwtwu3u/JZUV48PPdjnwdhNt7f4hq6OISK5RcLE9AGxkcWGf+bxeD1PL\nx4aPd2lQX0SkT44PLm1tkRZIUZEvbt6Z0yeF0+oaExHpm4KLrXurOEFwqZo2PpxWy0VEpG+ODy6t\n9pZLYfzgYh/U36UZYyIifXJ8cOlPy2W6bRuY+oZD9FgTAUREJJrjg4v94V+Jgsuk8WUUFHgBs5Py\n0eaTaa2biEiucnxwiR7QL4ibNy8vjypb60XjLiIivUt3cFkKbAN2Arf0kedO6/2NwKIky34W2Aq8\nCXx3MBVstU1FLk4w5gLRK/U1Y0xEpHeeNF7bDdwNXAE0AuuB1ZigELIMmAnMAs4D7gGWJCh7KbAc\nWAB0AuMGU8nW9uSnIkP0BpYa1BcR6V06Wy6LgTqgARMEHgWujMmzHHjISq8DyoCJCcp+CrjNOg9w\nZDCVbG8PhNOJxlwAdYuJiCQhncGlHNhnO95vnUsmz+Q4ZWcBFwN/A2qBcwZTyagB/WS6xaKCy2GC\nQT04TEQkVjq7xZKdp+vq53U9wChM99m5wP8CVb1lrKmpCaerq6uprq4+JU90t1j8AX0wDw4bPWok\nR4+14PcHaDx4jIrJY/r1C4iIZIva2lpqa2tTft10BpdGoMJ2XIFpgcTLM8XKkx+n7H7gd1Z6PRAE\nxgCnPB7SHlz6Yp8tlky3GJiV+kePtQDmwWEKLiKSq2L/8F65cmVKrpvObrFXMV1YlYAXuAYzKG+3\nGrjBSi8BmoFDCco+DlxmpWdb7w/4ucNRs8WSDS5R2+9rUF9EJFY6Wy5dwArgaczsr/sws71utt5f\nBTyJmTFWB7QCNyYoC3C/9doMBIgEpwGxr9BPtP1LyIxpmo4sIhJPOoMLwBrrZbcq5nhFP8qCmSX2\nD4OsV1h/9hYLqYoKLmq5iIjEcvQK/Z6enqjgkmy32PRpE8Bl5iHsf7sJv78zQQkREWdxdHBp7wiE\nn0Lp83nxeNxJlSss8DJ54mgAgt1B9uwf1FIbEZFhx9HBxT7eUliQ36+y9m1gtJhSRCSao4NLh607\nq9Da7ThZWqkvItI3RwcX+1iJz9fPlotmjImI9MnRwWUwLZeop1Lu0YwxERE7ZweXjsimlQW+/gWX\nivKxePLNTO7DR45zoqUtpXUTEcllzg4ug+gWc7vzqKyI7Pa/e+/hlNVLRCTXOTq42MdcCvoZXCC6\na0zjLiIiEY4OLh2DDS5aqS8i0isFF8tAgst0TUcWEemVo4OL3x8Z0Pf1c7YYnLqQsqcn2UfYiIgM\nb44OLu1Rs8X633IZP7aU4mLzgLGTJ9s5/M7xlNVNRCSXOTq4RM0W8/Z/g2iXyxWz3kUzxkREQMEl\nnPb1c51LSNXU8eH0Lg3qi4gADg8ug52KDNFPpdSgvoiI4ejgMpjtX0Lsg/r1Ci4iIoDjg8vgBvQh\nenfk3XsP09XVPeh6iYjkOkcHl8HsihxSMrKICePKAOjq7KJhnwb1RUTSHVyWAtuAncAtfeS503p/\nI7AoibI1wH5gg/VaOtDK+QNd4fRAWy4As2dODqe31x0Y8HVERIaLdAYXN3A35st/LnAdMCcmzzJg\nJjALuAm4J4myPcAPMIFoEfDUQCsY3S02sDEXgNNm2IJLvYKLiEg6g8tioA5oADqBR4ErY/IsBx6y\n0uuAMmBiEmVdqahgR8fgZ4sBzLYFl5273h5UnUREhoN0BpdyYJ/teL91Lpk8kxOU/SymG+0+TEAa\nkFSMuQCcZusW21F/gGAwOOBriYgMB/1flp68ZDfa6m8r5B7gP6z0t4DbgU/0lrGmpiacrq6uprq6\nOur9wW5cGTJ2dAljRo+k6WgLHR0B9h1oYtqUcYkLiohkWG1tLbW1tSm/bjqDSyNQYTuuwLRA4uWZ\nYuXJj1PWPh3rXuCJvipgDy69iQouA1znEjKrajJNR7cDpvWi4CIiuSD2D++VK1em5Lrp7BZ7FTNQ\nXwl4gWuA1TF5VgM3WOklQDNwKEHZSbbyVwGbB1K5YDAYvSvyAPYWsztNM8ZERMLS2XLpAlYAT2Nm\nf90HbAVutt5fBTyJmTFWB7QCNyYoC/Bd4ExMt9tu2/X6xT4N2efzkpc3uDhrnzG2QzPGRMTh0hlc\nANZYL7tVMccr+lEWIi2dQenoSF2rBaJnjG2vP0BPTw8uV0omtYmI5BzHrtAPdEZaLt4UBJeJ48so\nGVkEmGe7HDh0bNDXFBHJVQ4OLpE9wLzegc8UC3G5XFHjLjvVNSYiDubY4NJpG3PJ97hTcs3YrjER\nEadybHCxd4v5fKkZerIHl607Ymddi4g4h4ILkO9JTXCZM3tKOL1lx36t1BcRx1JwAbz5qQkukyeM\noqy0GIDW1g72NTal5LoiIrnGscGl0zagn+9NzZiLy+Vi7mmRjQXe2r4vTm4RkeHLwcEl9S0XgDPs\nwWXb3pRdV0Qklzg2uKRjzAXgjNMjwWXLTg3qi4gzOTe42Ld/ScEiypA5s6aAtTK/bvfBqJ0ARESc\nwrnBxd5yyU/NmAvAiOKC8I7Iwe4g2+oaU3ZtEZFcoeAC5KdwzAXgjNNsU5K3q2tMRJzHscHFvkI/\nlQP6gGaMiYjjOTe4dNn3Fkt1y8U2qL9DwUVEnMexwSV6tljqxlwAZlROxOczT7Y8fOQ4h440p/T6\nIiLZLl5wedCW/lia6zHkomeLDX5XZDu3O4/5cyKtl41vNaT0+iIi2S5ecFloS/9Luisy1OwtF08K\nZ4uFLJhbGU5vfLMh5dcXEclmju0WS9cK/ZCF8yrD6TfUchERh4n3rToFuBNwAeW2NJjn138uvVVL\nr0AaZ4sBzDutArfHTXdXNw17D3P8RBulJUUp/xwRkWwUr+XyFeA16/VV4HXb8WtJXn8psA3YCdzS\nR547rfc3Aov6UfZLQBAYnWRdoqRzthhAQYGX02eWh483bWlI+WeIiGSreN+qDw7y2m7gbuAKoBFY\nD6wGttryLANmArOA84B7gCVJlK0A3gXsGWjl/IHOcDrVs8VCFp5RGd688o03G7hoydy0fI6ISLaJ\nF1yewHR/uXp5rwdYnuDai4E6oME6fhS4kujgshx4yEqvA8qAicD0BGV/gGlN/SFBHfpk33Lfm+LZ\nYiELz5jGw781ac0YExEniRdclgD7gUcwX/wQPeaSSDlgX0G4H9M6SZSnHJgcp+yV1vGmJOrQJ/uY\nSyr3FrObP2daOL29/gBt7X6KCn1p+SwRkWwSL7hMwnQ9XWe9/oQJNG8lee1kAhD03jLqSyHwr1a9\nEpavqakJp6urq6murg4fd3ald0AfoLSkiKrKiexqOEiwO8ib2/axeNHMtHyWiMhA1NbWUltbm/Lr\nxvtW7QLWWC8fJsC8CNRgxkMSacSMjYRUYFoc8fJMsfLk91F2BlCJGfwP5X8N0wV3OLYC9uASKx2P\nOe7NmfMq2dVwEIDXNtYruIhIVon9w3vlypUpuW6idS4FwIeA/wE+A/wI+H2S134VM1BfCXiBazCD\n8nargRus9BKgGTgUp+ybwATMmMx0TMA5i14CSyJRU5HTMFss5NwzI8Hk1Tfq0/Y5IiLZJN636i+B\nM4Angf8ANvfz2l3ACuBpzOyv+zAD8jdb76+yrr0MM3jfCtyYoGysZLveThGwDeinest9u0Xzp+PK\ny6MnaJ7tovUuIuIE8b5VrwfaMC2Iz8e81wOUJHH9ULea3aqY4xX9KBurKok69KprCMZcAEaOKGTO\nrHK2bN8HPT1s2LyL6gvmpe3zRESyQbxusTxgBDCyl1cygSWr+f3pny0Wcq5tnGX9hrq0fpaISDaI\nF1wKgS9gBu9vJn4rJ+fYZ4ulelfkWOcsnBFOr9+ocRcRGf7iBZeHgLMxg+jLgNuHpEZDJGrMJU0r\n9EPmnT41/HyXxgNNHDh4NK2fJyKSafGCyxzgo8BPMTPGLh6SGg2BYDBIV+fQdYt5vR4WzossqHxt\n0660fp6ISKbFCy5dfaRznn3rF0++h7y89D954NyFkXGXV17fmfbPExHJpHjfqguAFttrvi19Iv1V\nS5+oHZHT3GoJWXyWLbhsqKPLVgcRkeEmXnBxEz1DzMMwmS0WvSPy0MxTmFE5kQnjygA4ebKdTVsG\nvKGziEjWc+STKO3dYj7f0AQXl8vF+eeeFj5eu377kHyuiEgmODK42PcV8wxRywXggsWnh9MKLiIy\nnDkyuHRGbVo5NGMuAGcvqApPSd6z7zD7324ass8WERlKjgwuQ7VpZSyfL5+zF0Z2rFHrRUSGK0cG\nl+jZYkO78UBU19gr24b0s0VEhoojg0smZouFnH/O7HD69c27OdnaMaSfLyIyFBwZXOyzxYayWwxg\nwrgyZs+cDEB3Vzd/XdfbkwRERHKbQ4PL0G390ptLbVvuv/B/bw7554uIpJsjg8tQPeK4L/bnuax7\nfSetbeoaE5HhxZnBJZDZ4DK1fCwzqyYB0Bno0qwxERl2nBlcOjMzFdkuqmvsr+oaE5HhxfHBJT8D\nLReASy+MBJeXX91BW7s/I/UQEUmHdAeXpcA2YCdwSx957rTe3wgsSqLst6y8bwDPARX9rVRXZ+bW\nuYRMmzKOqsqJAAQCneoaE5FhJZ3BxY15RPJSYC5wHeYBZHbLgJnALOAm4J4kyn4PWAicCTwOfLO/\nFbOvc8lUtxjAZbbWy5rnNmSsHiIiqZbO4LIYqAMagE7gUeDKmDzLMY9TBlgHlAETE5RtsZUfAbzT\n34plcoW+3bsvPTOcfmXDTpqOtcTJLSKSO9IZXMqBfbbj/da5ZPJMTlD228Be4GPAd/pbMftssXzP\n0K9zCSmfOJqFZ1QCEOwO8syLmzJWFxGRVErnn+09SeZzDeDat1qvrwE/BG7sLVNNTU04XV1dTXV1\nNRAzoJ/BbjGApZcvYuNbDQA89fwGrv3gBRmtj4g4S21tLbW1tSm/bjq/WRuJHmyvwLRA4uWZYuXJ\nT6IswMPAk31VwB5c7DozvIjS7rIL5/ODnz5BZ6CLnfUHqG84yAxroF9EJN3sf3gDrFy5MiXXTWe3\n2KuYgfpKwAtcA6yOybMauMFKLwGagUMJys6ylb8S6PdIeGcWzBYLGVFcwMVL5oaPNbAvIsNBOoNL\nF7ACeBrYAvwa2ArcbL3AtDp2YQbvVwGfTlAW4DZgM2YqcjXwpf5WzJ+h57n0ZellkRnYa557PWpM\nSEQkF6X7m3WN9bJbFXO8oh9lAa4ebKU6u7KnWwzgvLNmMX5cKYePHKf5eCsvvvwW77pkYaarJSIy\nYFqhnwXBxe3OY/l7zg0fP/7kKxmsjYjI4DkyuHQGImMumdhyvzfvf/c55LnNf4433tzN7j2HMlwj\nEZGBc2RwyfSW+70ZN6YkamD/92vUehGR3OXI4JJtYy4hH1y2OJx+6vkN2sxSRHKWI4OL35+dweXs\nBVVMKR8LQGtrB3985rUM10hEZGAcGVzse4tleoW+XV5eHh9Z/nfh48f+sJbu7mAGayQiMjDODC5R\nYy7ZMaAf8r4rzqK0pBiAAweP8uLatzJcIxGR/nNkcInauDKLusUACgq8UWMvj/z+r/T0JLtNm4hI\ndnBmcLG1XHxZ1C0W8qH3L8FjBb0t2/eFN7YUEckVjg8u+Z7sCy5jRo3kvZdFnvXy4KMvZLA2IiL9\n57jgEgwG6e7KvkWUsa6/+pLwosr1G+rYtGVPhmskIpI8xwUX+47I+V4PLtdAHieTfhWTx/Bu2/5i\n9z/8XAZrIyLSP44LLv5AZzidTWtcevOxay/FlafWi4jkHscFF/sal2wPLlPLx/Ke6kjr5We/eEYz\nx0QkJzguuNinIWd7cAHTegmNvWzYvIu167dnuEYiIok5L7jYZop5snQw325q+ViWv/uc8PFPHnhK\nq/ZFJOs5Lrh0Rq1xyc9gTZL38esvp7DQB0DD3sPac0xEsp7jgkvU6nxP9rdcwKx7uf5DF4WP7/vV\ns7S2dWSwRiIi8TkvuNj3FcvC1fl9ufaqCxk7pgSApqMt3PsrTU0Wkew1FMFlKbAN2Anc0keeO633\nNwKLkij738BWK//vgNJkK5NLs8XsCgu8fObjS8PHj61+mR31BzJYIxGRvqU7uLiBuzFBYi5wHTAn\nJs8yYCYwC7gJuCeJsn8GzgAWAjuArydboWzetDKRd12ykLMXzgCgJxjk9ntWEwxqcF9Esk+6g8ti\noA5oADqBR4ErY/IsBx6y0uuAMmBigrLPAEFbmSnJVsg+oJ+tW7/0xeVy8aVPLcdtjRW9uXWvBvdF\nJCulO7iUA/tsx/utc8nkmZxEWYCPA08mW6FADs4Ws5tWMS5qcP/u+9Zw6EhzBmskInKqdPcLJbuc\nfKAbfN0KBICHe3uzpqYmnK6urqa6ujonZ4vF+thHqnn2pc0ceLuJ1tYObvvR7/jht27M2n3SRCR7\n1dbWUltbm/Lrpju4NAIVtuMKTAskXp4pVp78BGX/ETNec3lfH24PLiEB28aVuTRbzK6gwMs3vng1\nn/rqz6Cnh/Ub6nh8zStctey8TFdNRHJM6A/vkJUrV6bkuunuFnsVM1BfCXiBa4DVMXlWAzdY6SVA\nM3AoQdmlwFcwYzD9WvAR9SyXHBvQt1swdxrXXXVh+Pju+9aw70BTBmskIhKR7uDSBawAnga2AL/G\nTCG+2XqBGS/ZhRm8XwV8OkFZgLuAEZiB/Q3AT5KuUGdu7S0Wzz9/9AqmVYwHoKMjwL9/9xH8/s4E\npURE0m8ovl3XWC+7VTHHK/pRFkyLZkD8gdydLRbL58vnm1/+MDd9eRVdnV3sqDvAnfc+yVc+Ezsh\nT0RkaDluhb59EaUvR8dc7E6bWc5nP/He8PHjT67j2Zc2ZbBGIiIODC5RYy6e3A8uAB96/xIuvXB+\n+Pg7d/6e+oaDGayRiDid84JLIDf3FovH5XLx9c9fRfnkMQC0t/v56spfcrT5ZIZrJiJO5bjg0jlM\nZovFKi4q4Du3fpSiIrM1/8HDx/jXb/8qKpiKiAwVxwWXqF2Rc3xAP1ZV5QRqvnINWIspN2/Zw7fv\n+K32HxORIee84JJjjznurwsWn84K2wD/sy9u5I6f/YmenmQ3SxARGTzHBZeoLfeHyZhLrGs/eAEf\ntK3W/+0TL3P/w89nsEYi4jSOCy6BQGSRoWeYzBaL5XK5+OInP8DlFy8In7v/4ed49PH/y2CtRMRJ\nHBdcOjuH1zqXvrjdeXzji1ez+OzIetO7fv4nHv7tXzJYKxFxCscFl0AOP8+lv/LzPfzXv17PgjMq\nw+d+fP8aHnz0hcxVSkQcwXHBxW/rFsvF57n0V2GBl9tXfoxF86vC537+y2f4yQNPaRaZiKSN44JL\nIKpbbPgHF4CiQh/fr7mBcxfNDJ/71W9eYuX3H9M6GBFJC+cFF9uuwcO9W8yuoMDLd7/xD1xw3pzw\nuWdf3MgX//1BTrS0ZbBmIjIcOS+42B9z7HNGyyXE58vntluvj5qmvGHzLj7xhXu0F5mIpJSjg8tw\nXESZiNudx5c/vZxP37g0fO7A20388xd/yjMvbsxgzURkOHFccOnwO2tAvzcul4vrr76Y//z631NQ\n4AXA7w9Q871fc8fP/qhxGBEZNEcFl+7uIN22FfpOGnPpzaUXzuPnt3+Kiinjwuce+8Na/umLP2FX\nw6EM1kxEcp2jgot9R2SvNx+XtcGjk1VVTuDeH3yKi86fGz5Xv/sgH//Cj/nf1Ws1XVlEBmQogstS\nYBuwE7iljzx3Wu9vBBYlUfbDwFtAN3BWshXxD8NnuaTCiOICbrv1er7wyQ/gtboKOwNd/GjVH/nM\n1+6lYe/hDNdQRHJNuoOLG7gbEyTmAtcBc2LyLANmArOAm4B7kii7GbgKeKk/lYmaKabgEsXlcnH1\nB87n/js+zcyqSeHzm95q4IbP3sV9v3pOYzEikrR0B5fFQB3QAHQCjwJXxuRZDjxkpdcBZcDEBGW3\nATv6Wxmnrc4fiOnTJvDz2z/FDddU4/aYManurm7uf/g5bvjsXfx13VZt3y8iCaU7uJQD+2zH+61z\nyeSZnETZfgn4h+dTKFPN6/Vw8w3v5oEffYa5p1WEz+/bf4Rb/uOXfP7W+9lRfyCDNRSRbJfu4JLs\nn7hDMrIevYBSwSWRGZUT+el/38wXPvkBiosLwudf21jPjZ//Md/+4W9oPHg0gzUUkWyV7m/YRqDC\ndlyBaYHEyzPFypOfRNm4ampqwunq6mpGj58ePnbiAsqBcLvzuPoD53P5RfO5/5HneXzNKwS7g9DT\nw5PPvs5TL7zBssvP4oZrqimfODrT1RWRfqqtraW2tjbl1013i8EDbAcuBw4Ar2AG5rfa8iwDVlg/\nlwB3WD+TKfsC8GXgtV4+uyd2bGD9G3X8y633A3DWgiruuu2fBvXLOVHD3sPcdd8a/vbq9qjzbo+b\n91Qv5LqrLqKqckKGaicig2Ut0Rh0bEj3n+9dmMDxNGb2132Y4HCz9f4q4ElMYKkDWoEbE5QFM1Ps\nTmAs8CdgAxB5cHwf/P7odS7Sf5VTx3P7yo+xYfNu7n/4OV7ftAswg/5PPvs6Tz77OovPnsW1H7yQ\nxYtmai2RiEMN53/5p7RcXvjrm/zbbQ8DcMnfncF/3Xp9Juo1rLy+aRcPPPJ8OMjYTZ82gSuXLuY9\nly6kZGRRBmonIv2VKy2XrGKfiqxFlKlx1oIqzlpQxaYte/j14//Hiy9vocda1b97zyHuWPUEP3ng\nKS69cB7L33MOC8+oVGtGxAEc9Q0bvYhS3WKptGDuNBbMnUbjwaM8tnotf/zza7S3+wEIBDp5+vkN\nPP38BqaUj+VdFy/gXZcsZFrFuARXFZFcNZz/hDylW+yx1S9zx6onALjqfUv48qeXZ6JejnCytYNn\nXtzI6qd5fvtuAAAPpklEQVTXs6Ou9zUxs2dO5oqLFnDZRfOZNGHUENdQRHqjbrEB8PsD4XSBwx4U\nNtRGFBdw1bLzuGrZeWzb2cgTT6/nzy9upK3NH86zo+4AO+oO8JMHnmJm1SQuXHw6F543h9NmTiYv\nz1F7qooMO44KLm0dkeBSVOjNYE2c5fRZ5Zw+q5zP/fP7WLt+O8+8uJG1r26n07ZXWd2ut6nb9TYP\nPvoCY8eUcOHi01lyzmwWza9ihG0Bp4jkBkcFl3ZbcCks9GWwJs7k8+Vz6YXzuPTCeZxs7eCll7fw\n7EubeG3TLrps42HvNJ3g8TWv8PiaV3Dl5TF39hTOOXMG5yycwbzTp2oyhkgOcNS/0tAAM0BhgVou\nmTSiuIBlV5zFsivOorWtg/Ub6vjL37aydv12TrS0hfP1BIO8tW0vb23by0OPvoDP52Xe6RXMmzOV\n+XOmMX/OVLVsRLKQs4JLR2QqsrrFskdxUQHVF8yj+oJ5dHcH2bx1D2vXb+e1jbvYVtcItokZfn+A\n1zbW89rGenPC5aJq2gQWzJnK/LnTOH1mORXlY3G7NWYjkkmOCi5ttpZLkbrFspLbnceZ86Zz5jyz\nD9zxE21s2LyL9RvqWL+xnsYDTdEFenrY1XCQXQ0HeXzNKwD4fF5mz5jE7BmTOW3GZE6fVU5lxXgF\nHJEh5Kjg0t5umy2mbrGcUFpSFG7VABw60szmrXvZtGUPm7fuYeeug+FFmyF+f4DNW/awecue8Dmv\nN5/KqeOomjqB6VPHUzl1AtOnjWfS+DLNTBNJA2cFF/tsMQWXnDRhXBkTxpVxxcULANMafWv7PjZv\n2cPWnY1sr2uk6WjLKeUCgc7w1Ge7ggIvlRXjqZo2nqlTxjFl0himTDYvjcuJDJyjgkubBvSHnaJC\nH+eeOZNzz5wZPvfO0RNsrzvAjvoDbKtrZGf92xw60txr+Y6OANt27mfbzlOf5jBm9EgqJo8NB5sp\nk8YwacIoJo4vo7SkSNvYiMThqODS4ddUZCcYO7qEsYtLuGDx6eFzzcdb2bXnELv3Hmb33sM07D3M\nrj2HOH6itc/rNB1toeloC2+8ufuU93w+LxPGlTJhfBmTxpvW1ETr54RxpYwbU6KnnYqjOer//tY2\ntVycqqy0OLzJpt3R5pPstoLO/reb2NfYROPBoxw4eJTuru4+r+f3B9i7/wh79x/pM09pSTFjRo9k\n7JiRjB01krFjShgzaiRjR48050ebY63bkeHIMf9X9/T0RI+5aCqyAKPLRjC6bARnL5wRdb67O8jB\nI83sa3yH/QdMwNnX+A6Hjhzn4OFjUdvY9OX4iVaOn2hlV8PBuPlGjiykrHQEo0qLKS0poqykiLLS\nEZSVFlvpYspKiikrLWZUaTE+bV0kOcAxwSUQ6DKP58U8NVFdFhKP251H+cTR5tHNZ5/6fsvJdg4e\nbubQkWYOHm6OSh860szR5tZTZrH1paWlnZaWdvbFaQXZFRR4KSspZuTIQkYUFTByZCElIwoZUVxA\nychCRhQXhn+OHFHAyBGFjBxh8ng87v7cBpEBc8w3rFotkkqhL+xZVZN6fb+7O8jR5pM0HW3hnaMn\neMcav2k61sI7R1t4p+kETcdaaDp2MukgFNLREeBgR4CDh4/1u94FBV5GFBdQVOijuMhHUaGPoiIf\nxdZP+/nY94uLCigs9FJc5KOwwKsp3BKXY4JLS2tHOF1cqO1CJL3c7jzGjSlh3JgSoLzPfN3dQY63\ntHH8eCvHjrdy/EQbx4630nz8JM0n2jh+oo3m4yetc600n2iLOxaUSEdHgA7bH1qD4fN5KfDlU+DL\nx+fLp7DAOrZ+ho59tmN7vtBxgS8fn9dcw5vvwev14POadH6+W7PycpRjgsux5pPh9Kiy4gzWRCTC\n7c4Lj/tMTyJ/T08PrW1+jre0me60VtOlduJkOydPmp+trR2cONlOy8l2Wk52hN9rae3odyspHr8/\ngN8f4HjKrtg7rzffCjgeK/jk4/NZaevYa73n83rIz/eEA1WoTL7XQ77HTb7Hjcf6mZ/vMen8yHG+\nx43H7cYTPucm3+PB48kL/1SLLTnpDi5LgTsAN3Av8N1e8twJvBdoA/4R2JCg7Gjg18A0oAH4CND7\nIgYb+5TTslIFF8lNLpeLEcUFZrPOif0rGwwGaWsP0NrWQVubn9Z2P23tAdraOmi1jtvb/Sbd5qet\n3XqFjjtM3pOt/qhnI6VbINBJINDJycRZh4TbCkD5tgDk8XjwuPPCwSg/3wpSnjzcbjdudx4edx5u\n6+Wxn/O4o86581x4bOfyrHyhc+68PHPdPNu50Gfl9Z4vL1Q+z3x+Xp7LXDvPFc4bOp8q6QwubuBu\n4AqgEVgPrAa22vIsA2YCs4DzgHuAJQnKfg14BvgecIt1/LVElTnaHAkuoxwWXGpra6murs50NbKC\nk+9FXl5eJDAxuHsRDAbxB7ro6AjQ3hHA7++kvSNAh7+TDn9n9HEoT6Dz1Ly2MoHOLjoDXfgDJu0P\ndA2qC7A/jjftoXTMtKTydnd1093VjT/xhEFHS2dwWQzUYVoXAI8CVxIdXJYDD1npdUAZ5u+x6XHK\nLgcusc4/BNSSRHCxD36OHVPSr18k1zn5CzWW7kXEYO5FXl4ehQVeCgu8pPMB1d3dQTqtQBPotF7+\nSPDxBzoJBLqi8vj9nXRaef2BLgJWvtC1Oru66ezspqvb+tnVzdoXNzP3tArrPZOny5ans6vbOteV\nuNICpDe4lAP7bMf7Ma2TRHnKgclxyk4ADlnpQ9Zxr76y8hfh9NpXtoXTlRXjk6i+iGSa6fLxpn2j\n2ZqafdTUfCphvp6eHrqsQNMZE3y6QkGpKxKMuruD1suku7qDdHVF0uanOQ5a6a6uIN3BU/MFY851\ndXVH8gWt8vayQXMcSnd3BQn2BCN1CgYJBnvCZYLWNVIlncGlJ3EWAJLp5HP1cb2eeJ9jDyh2s/uY\nPioiEo/L5TID//keKMx0bVKvp6eHvLzbMl2NhJYAT9mOv44ZI7H7KXCt7XgbpiUSr+w2IkOZk6zj\n3tQRCT566aWXXnol96ojy3mAeqAS8AJvAHNi8iwDnrTSS4C/JVE2NJAPZqzlOymvuYiIZLX3Atsx\nkfDr1rmbrVfI3db7G4GzEpQFMxX5WWAH8GfMJAAREREREZHcshQzDrOTU8d4hqMK4AXgLeBN4HPW\n+dGY9UC9tfC+jrk/24B3D1lNh44bsxj3CevYqfeiDPgNZgr/FsyMS6fei69j/o1sBh4GfDjnXtyP\nmVm72XZuIL/72dY1dgI/SmN9s5Ib041WCeTT+zjPcDMRONNKj8B0Jc7BjE191Tp/C5GxqbmY+5KP\nuU91wHDbz+KLwK8wC2/BuffiIeDjVtoDlOLMe1EJ7MIEFDA7fHwM59yLi4BFRAeX/vzuoRm9r2DW\nL4IZK1+athpnofOJnmWW1Or9YeZxzM4GoZl3YAJQaFZd7Ky9pzCTKYaLKZgxuUuJtFyceC9KMV+o\nsZx4L0Zj/ugahQmyTwDvwln3opLo4NLf330S0Qvgr8XM9u1TLkfj3vS1KNMpKjF/oayj78WmkzH3\nJWS43aMfAl8B7KvBnHgvpgNHgAeA14GfA8U4814cBW4H9gIHMHsRPoMz70VIf3/32PONJLgnwy24\n9GS6Ahk0Avgt8HmgJea90Pz1vgyX+/Z+4DBmvKWvxblOuRcezOzLn1g/Wzm1Fe+UezED+BfMH1+T\nMf9WPhqTxyn3ojeJfvcBGW7BpREzwB1SQXS0Ha7yMYHll5huMTB/jdgXmx620rH3aIp1bjj4O8ze\nc7uBR4DLMPfEifdiv/Vabx3/BhNkDuK8e3EOsBZoArqA32G60J14L0L6829iv3V+Ssz54XZP4kpm\n4eZw4wJ+gekOsutrsWlowM6L6TqpJ7kteHLNJUTGXJx6L14CZlvpGsx9cOK9WIiZSVmI+Z0eAj6D\ns+5FJacO6Pf3d1+HmXHowoED+tD34svh6kLM+MIbmO6gDZj/6PEWm/4r5v5sA94zlJUdQpcQmS3m\n1HuxENNy2Yj5a70U596LrxKZivwQprXvlHvxCGasKYAZk76Rgf3uoanIdZjncImIiIiIiIiIiIiI\niIiIiIiIiIiIiIiISEQQ+L7t+MvAN1N07QeBD6XoWvF8GLMd/nMx56uJLAIVyTnDbfsXcZYAcBUw\nxjpO5f5Ig7mWpx95PwH8E3D5ID6vN/q3LRml/wEll3UCPwO+0Mt7DxLd8jhp/awGXsTswVaP2fbi\nHzDPqtgEVNnKXIFZ4b4deJ91zg38t5V/I3CT7bp/Af6AWQke6zrr+puJbLXx78AFmIc5fS8mfw9m\ng8XHMFud/491/jLg97Z878Ksvg/9jt/H7NZwPvANq56bgVVWHo917hLr+DbgP630RzFbfGzAbKee\nZ/2+D1rX2ITZAFJEZFhrAUZiNqosAb5EpFvsAaKDS2in6GrgGGaLcS9m870a673PEdmj7UHM/kkA\nMzHbZvgwweRW67wPE3wqreueBKb1Us/JwB5MC8uN6QK70nrvBcyGkrGqMVvDT8bs5bQWszEnmGAT\naq09TCTwBYGrbdcYZUv/ArNrNJj9o7ZggufrmIAzB7NdjtvK82NM0D0Lsz1ISGkvdRU5hVoukuta\nMF+cn0uU0WY9ZlfYAGafpKet829iAgWYlsP/Wuk6zIO3Tsc89vUGzF/3f8Ps0TTTyvcKJojEOhcT\nRJqAbsxTMi+2vd/XpoivYPaE6sG0RqZb53+J+eIvwzzIaY11vhuzO3bIZVYdN1npM6zzWzAtoScw\n+0x1YbrlzgZetX63y63P24Vpzd2J2WfqRB91FYnSn75hkWx1B+Yv8Ads57qI/PGUh2mlhPht6aDt\nOEj8fxOhcZgVmIdN2VVjnpnSVzl7AHERPabT1/iOvZ7dtro9gAkMHZgAGHowWoftWgWY1sfZmNbZ\nN61zIfOJtOBCHsJsWhhrAWYz1E8CH8GME4nEpZaLDAfHMF+ynyDy5dqA+WIF84yX/H5e04WZyeXC\nPGyqCrNL7NPAp4l80c8GihJcaz1mjCPULXYtZtxnoN7GtGj+jeiAahcKJE2YsZsPE7k3/w/T6rkE\nuAvT1fUcpkttnJVnNDDVqrMHM67zDXrvwhM5hVouksvsf/HfjmlRhPwcM7j+BuY54Cf7KBd7vR5b\nei+ma6oEuBnTjXYvpuvsdUzgOYyZsRbvaX5vY56Z8YJV5o8knmbc2/Xsxw8DYzGTDXp7vxlzD97E\nPBRrnXV+DGYQ/zJMi+Zu4EfAP2KC1Z8xf3R2YoJoByaAhf4QjX2apYiIDCN3Y8ZLREREUuI1oJb+\nd/WJiIiIiIiIiIiIiIiIiIiIiIiIiIiIiGSz/w8G8ibXdDINzAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa7e16cfa50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import thinkplot\n",
    "thinkplot.Pdf(suite)\n",
    "thinkplot.Config(xlabel='Number of hyraxes', ylabel='PMF', legend=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here are some summaries of the posterior distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior mean 185.570957948\n",
      "Maximum a posteriori estimate 50\n",
      "90% credible interval (36, 618)\n"
     ]
    }
   ],
   "source": [
    "print('Posterior mean', suite.mean())\n",
    "print('Maximum a posteriori estimate', suite.MaximumLikelihood())\n",
    "print('90% credible interval', suite.CredibleInterval(90))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The combinatorial expression we computed is the PMF of the hypergeometric distribution, so we can also compute it using `thinkbayes2.EvalHypergeomPmf`, which uses `scipy.stats.hypergeom.pmf`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import thinkbayes2\n",
    "\n",
    "class Hyrax2(thinkbayes2.Suite):\n",
    "    \"\"\"Represents hypotheses about how many hyraxes 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 = thinkbayes2.EvalHypergeomPmf(k, N, K, n)\n",
    "        return like\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the result is the same:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEPCAYAAACOU4kjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XOV97/HPaEajxbIk77Is2bKxDcZsZjGGsIgljXFb\n09wshNtmoc2Fe1OXNM1CaJIip71JSZM0obTUNyGJk5StWQgpEEiIRUITDAZjG++bbEm25U2ydo1m\nuX88Z2bOjKSZkTSjmdH5vl+veeksz3P06IDPb571gIiIiIiIiIiIiIiIiIiIiIiIiIiIiEjOWQ3s\nAfYD942Q5iHr/DZgRQp5VwKvAVuB14Gr0ltkERHJZW7gAFAHFAJvAcvi0qwBnrO2rwZeTSFvI/Au\na/s2YFO6Cy4iIuNTkMFrr8QEiCZgEHgCuD0uzVpgo7W9GagEqpLkPQ5UWNuVQGsmCi8iImPnyeC1\n5wHNtv0WTO0kWZp5QHWCvJ8FXgG+igmO16SvyCIikg6ZrLmEUkznGuV1HwXuBeYDnwC+M8r8IiKS\nYZmsubQCtbb9WkwNJFGaGitNYYK8K4Fbre0fAd8e7pefd955oYMHD46p4CIiDnYQWJztQiTiwRSy\nDvCSvEN/FdEO/UR53wRutLZvwYwYG05IjAceeCDbRcgZuhdRuhdRuhdRpN7qlFAmay5+YB3wAmb0\n16PAbuAe6/wGTGBZg+m87wHuSpIX4G7gX4EioM/aFxGRHJLJ4ALwvPWx2xC3v24UeQG2MHRggIiI\n5JBMduhLjqivr892EXKG7kWU7kWU7kX6jXakVj6xmg9FRCRVLpcL0hAbVHMREZG0U3AREZG0U3AR\nEZG0U3AREZG0U3AREZG0U3AREZG0U3AREZG0U3AREZG0U3AREZG0U3AREZG0U3AREZG0y/SqyHln\nYMDP08/tpeVYJ7feuJBLL6rKdpFERPKOFq6M880Nm3n5v48AZgG3L95fz/ILZqW7bCIiOUkLV2bA\niZPdkcACEAqFePzHb2exRCIi+UnBxeb1N48NObZr7ykONbVnoTQiIvlLwcVmz/7Twx5/dUvrBJdE\nRCS/KbhYQqEQe/adiey//0+WR7Y3v6HgIiIyGpkOLquBPcB+4L4R0jxknd8GrEgh7xPAVutz2Po5\nbu3n+mk/1wdASXEha29bisdjbk9z6znaO/rT8WtERBwhk8HFDTyMCRIXAncCy+LSrAEWA0uAu4FH\nUsj7AUwQWgH82PqM27HjXZHtmupySksKWbJoRuTYrr2n0vFrREQcIZPBZSVwAGgCBjE1jtvj0qwF\nNlrbm4FKoCrFvC7g/cDj6Sjs8bbuyPbcOWUAMUOQFVxERFKXyeAyD2i27bdYx1JJU51C3uuBNuBg\nOgprDy5VkeAyM3Js5x4FFxGRVGUyuKQ6g3Gsk3XuBB4bY94h7M1i1VUmuJy/eCYFBaZ4R1s66e0b\nTNevExGZ1DK5/EsrUGvbr8XUQBKlqbHSFCbJ6wHeDVyeqAANDQ2R7fr6eurr60dMe/psb2R79qwp\nABQXe6idV8GR5g4gxKGmdi5aNjvRrxQRySuNjY00Njam/bqZDC5bMB31dcAx4A5MbcPuGWAdpk9l\nFdCBaeo6kyTvrcBu69yI7MElmfb26GiwGdNKItuLF06zggsKLiIy6cR/8V6/fn1arpvJZjE/JnC8\nAOwCnsQEhHusD8BzwCFM5/0G4GNJ8obdQZo68gECgSAdndHgUllZHNleVDctsn3gsGbqi4ikItOr\nIj9vfew2xO2vG0XesLvGU6h45zoHCC9yWT61iEKPO3Ju8cJocDmo4CIikhLN0AfOdvRFtqfbmsQA\nFtRW4nab23S8rYueXt+Elk1EJB8puBDb3zLN1iQG4PW6WVBTEdk/1NQxYeUSEclXCi7E1lymVZQM\nOb9wQWVkO9y5LyIiI1NwAbq6ByLb5eVFQ84vqI3WXI40n5uQMomI5DMFF6CrO9qPMrXMO+T8fFuz\nWJOCi4hIUgouQGeXreZSlrjm0tzaSSAQnJByiYjkKwUXoDtJzaWivDjSF+Pz+TlxsntIGhERiVJw\nIbZZrGyY4AIw31Z7OdrSmfEyiYjkMwUX4vtchjaLAdTFdOprxJiISCIKLsSOFhuuWQxiay5NR9Wp\nLyKSiOODSyAQpLsnupT+SMHFPpHySIuCi4hIIo4PLuYdLWZdsdKSwshSL/Fq5pVH3u3SdrKbvn69\n20VEZCSODy49vdEgMaV0+FoLgLfQzby5UyP76tQXERmZ44NLX58/sl1SkniR6AW19mVg1DQmIjIS\nBRdb81ZpSWHCtPNryiPbzep3EREZkYJL/yhqLjW2mouCi4jIiBRc+qI1l5LixDWXWlvN5UjzucgL\nxkREJJbjg4u9Qz9Zs9jsmVMoKjK1m67uATrODSRMLyLiVI4PLqNpFisocDF/nm0ZmFY1jYmIDCfT\nwWU1sAfYD9w3QpqHrPPbgBUp5v0rYDfwNvDgeAo4mmYxiF0h+ahGjImIDCvxV/XxcQMPA7cCrcDr\nwDOYoBC2BlgMLAGuBh4BViXJexOwFrgEGARmjaeQvX32ZrHkt2O+XhwmIpJUJmsuK4EDQBMmCDwB\n3B6XZi2w0dreDFQCVUny/h/gy9ZxgFPjKWR/TLNY8prL/HnRTn01i4mIDC+TwWUe0Gzbb7GOpZKm\nOkHeJcANwKtAI3DleAoZW3MZZbNYSyfBoEaMiYjEy2SzWKpPXdcor+sBpmGaz64CngIWDZewoaEh\nsl1fX099ff2QNL32GfrFyW9HRXkxlRXFdJzrx+fz03aym7lVU5PmExHJRY2NjTQ2Nqb9upkMLq1A\nrW2/FlMDSZSmxkpTmCBvC/ATa/t1IAjMAM7EF8AeXEYS06GfQs0FYP68CjrO9QNwtLVTwUVE8lb8\nF+/169en5bqZbBbbgmnCqgO8wB2YTnm7Z4APWdurgA6gLUnep4Gbre2l1vkhgSVVscu/pBZr1akv\nIpJYJmsufmAd8AJm9NejmNFe91jnNwDPYUaMHQB6gLuS5AX4jvXZAfiIBqcxiWkWS7HmEvNuFwUX\nEZEhMhlcAJ63PnYb4vbXjSIvmFFiHxxnuSJi57mkWHOxLQNzVGuMiYgM4egZ+qFQKK5DP7WaS+28\nCsLjEI63deHzBTJRPBGRvOXo4NI/4Cc8qM3r9eDxpHY7ios9zJk9BYBgMETLMb04TETEztHBxb6u\nWHGRe1R57f0umkwpIhLL0cHFNxBtziouGl330/ya2MmUIiIS5ejgMuCL1lyKRltz0QKWIiIjcnZw\nsdVcisZRc9FbKUVEYjk7uNhGeRV5R1dzqa4qiwwAOHO2l+5uX1rLJiKSzxweXGzNYt7R1Vzc7gJq\nqrVCsojIcJwdXGKaxUZXc4G4EWNqGhMRiXB0cLFPfvSOslkMtMaYiMhIHB1cxtMsBurUFxEZibOD\nSxqbxZpbOgmF9OIwERFweHAxy78YY6m5zJheEnl7ZU+vjzNn+9JWNhGRfObo4DIQ0+cy+lvhcrli\n+l3UqS8iYjg6uPjG2ecC5q2UYep3ERExHB1c7H0uYxktBnHLwCi4iIgATg8uvrEvXBmmNcZERIZy\ndnAZGPvClWH2ZrHmY534/cFxl0tEJN85O7jErC02tppLWZmXmTNKAfD7g3pxmIgImQ8uq4E9wH7g\nvhHSPGSd3wasSCFvA9ACbLU+q8daON/g+Oa5hC1aMC2yfaipfczXERGZLDIZXNzAw5iH/4XAncCy\nuDRrgMXAEuBu4JEU8oaAr2MC0QrgF2MtoC8NHfoAi+psweWIgouISCaDy0rgANAEDAJPALfHpVkL\nbLS2NwOVQFUKeV3pKGA6msUAFtVVRrYPHekYV5lERCaDTAaXeUCzbb/FOpZKmuokef8K04z2KCYg\njcl43kRpZ6+5HD7SQTCoZWBExNnG/nU9uVSfsKOthTwCfNHa/nvga8BfDJewoaEhsl1fX099fX3M\n+Zi1xcbRLDa9soRpFSW0n+tjYMDP8bYu5s0tT55RRCTLGhsbaWxsTPt1MxlcWoFa234tpgaSKE2N\nlaYwQd6TtuPfBn4+UgHswWU46WoWA1i4oJL27WZtsUNNHQouIpIX4r94r1+/Pi3XzWSz2BZMR30d\n4AXuAJ6JS/MM8CFrexXQAbQlyTvXlv/dwI6xFC4YDMUs/zKeDn2I69TXiDERcbhM1lz8wDrgBczo\nr0eB3cA91vkNwHOYEWMHgB7griR5AR4ELsM0ux22XW9UYl8U5qGgYHxjBGI79RVcRMTZMhlcAJ63\nPnYb4vbXjSIvRGs642LvzPcWjr8CF1tz6SAUCuFypWVQm4hI3nHsDP3BwegyLYWF42sSA5g1o5Sy\nKV7AvNul7VTPuK8pIpKvHBxc0jOBMszlcg0Zkiwi4lTODS62BSYLPem5DerUFxExHBtc7OuKedPQ\nLAawaEG0U3//obNpuaaISD5ybHCx97l40lRzWbxoemT7wKGzmqkvIo7l4OASrbkUpmG0GMCcWVMo\nn1oEQG/fIMdOdKXluiIi+caxwcX+Uq90jBYD06m/9LwZkf39B9U0JiLO5NjgkokOfYAl50WbxvYe\nOJO264qI5BPnBhdbs5jHk56aC8DSxdGaywF16ouIQzk4uERrLl5v+m7D4oXTCC/03NTcwcCAP3EG\nEZFJyLHBxT4UuTCNNZcppV5qqqcCZnHMg5rvIiIO5NjgEtPnkqbRYmH2fpd96tQXEQdybnDJUM0F\niBsxpk59EXEexwYXfyZrLotUcxERZ3NscMnEDP2wBbUVeK03W54528vpM71pvb6ISK5L9FT9nm37\nwxkux4RL96rIdm53ARfYhiTv2nsqrdcXEcl1iYLLpbbtv850QSbaoN/e55L+CtwFS2dGtnftPZ32\n64uI5DI1i5G+5V/sLjw/Glx271PNRUScJdFrjmuAhzAzAufZtsG8v/7ezBYts2LnuaQ/xp6/eAZu\ndwGBQJDm1k66ugeYWlaU9t8jIpKLEj1VPw28YX0+A7xp238jxeuvBvYA+4H7RkjzkHV+G7BiFHk/\nCQSB6cOcSyoTC1faFRV5OM/28rDdahoTEQdJVHP53jiv7QYeBm4FWoHXgWeA3bY0a4DFwBLgauAR\nYFUKeWuBdwJHxlo4n8++tlhmWgcvPH8m+6x5Lrv2nmblFfMy8ntERHJNouDyc0zzl2uYcyFgbZJr\nrwQOAE3W/hPA7cQGl7XARmt7M1AJVAELk+T9OqY29bMkZRiRP2BfWyz9NReAZUtn8fRzewHYpX4X\nEXGQRMFlFdACPI558ENsn0sy84Bm234LpnaSLM08oDpB3tut/e0plGFEE1FzuWBpdDjyoaZ2+voH\nKSkuzMjvEhHJJYmCy1xM09Od1udZTKDZmeK1U33H73A1o5GUAH9rlStp/oaGhsh2fX099fX1kX1/\nht7nYje1rIj5NRUcbTlHMBhi34EzXHpRVUZ+l4jIWDQ2NtLY2Jj26yYKLn7geetThAkwLwMNmP6Q\nZFoxfSNhtZgaR6I0NVaawhHyngfUYTr/w+nfwDTBnYwvgD24xItZcj8DHfphyy+YxdGWcwBs33VS\nwUVEckr8F+/169en5brJvrIXA+8Bfgj8JfBN4KcpXnsLpqO+DvACd2A65e2eAT5kba8COoC2BHnf\nBuZg+mQWYgLO5QwTWJKJmUSZ5rXF7C5ZPieyvX3nqIspIpKXEtVcfgAsB54DvgjsGOW1/cA64AXM\n6K9HMR3y91jnN1jXXoPpvO8B7kqSN16qTW9D+DI8iTJs+QWzcLlchEIhDh5u13wXEXGERMHlT4Fe\nTA3i43HnQkB5CtcPN6vZbYjbXzeKvPEWpVCGYU1EnwtA2RQvixdOZ/+hM0CIt3ef4pqrajL2+0RE\nckGip2oBUAZMHeaTSmDJabGjxTJXcwG49KJo09i2t9sy+rtERHJBouBSAnwC03l/D4lrOXnHXnPx\nejO7xNoly2dHtrfvUnARkckv0VN1I3AFphN9DfC1CSnRBIldFTmzNZeli2dE3u9yoq2btlM9Gf19\nIiLZlii4LAP+DPh3zIixGyakRBMgGAzF1FwyNYkyzFvojlkleYdqLyIyySV6qvpH2M578YGloGA0\n8zjHxj4k+a0dCi4iMrklCi6XAF22z8W27c7MFy1zJrJJLOyyuE59e4ATEZlsEgUXN7EjxDxMktFi\nPt/ENYmFLaitYOaMUgB6en3s2a8l+EVk8nLkmyj9tppLplZEjudyubji0rmR/TfeOj4hv1dEJBsc\nGVwGJ7Az3+7KFdWR7Te2KbiIyOTlzOBie8WxN4PrisW7eNnsyJDklmOdnGjrnrDfLSIykRwaXOxL\nv0xMsxiYJrhLLoxOqNyi2ouITFLODC72dcUmsOYCcMVl0X6XLVuPTejvFhGZKI4MLr7BzL+FciT2\nTv2de07R0+ub0N8vIjIRHBlcYtcVm7hmMYCZM0pZtGAaAIFAkNdVexGRSciRwcXeoe9xT/wtsC+5\n//vX4l/OKSKS/xwaXCbmRWEjuWZlNLi8taON3r7BCS+DiEgmOTS4TMwrjkdSXTWVuvmVpiz+gCZU\nisik48zgYu9zyULNBeKaxl5X05iITC7ODC5ZrrkAXGtrGntz23H6+tU0JiKTR6afrKuBPcB+4L4R\n0jxknd8GrEgh799bad8CXgJqR1uomHkuEziJ0m7e3HLm11QAZmi0msZEZDLJZHBxY16RvBq4ELgT\n8wIyuzXAYmAJcDfwSAp5vwJcClwGPA08MNqC+Xz2mkt2ggvAtSujcbHxlSNZK4eISLplMrisBA4A\nTcAg8ARwe1yatZjXKQNsBiqBqiR5u2z5y4BRr13vz+IMfbsbrp0f2X7r7RO0d/RnrSwiIumUySfr\nPKDZtt9iHUslTXWSvP8XOAp8GPjH0RbMPhR5omfo21XNLmPZ0lmAefXyK68ezVpZRETSyZPBa4dS\nTDeWdwx/zvp8Fvhn4K7hEjU0NES26+vrqa+vB+LeRJnFZjGA+usWsHvfKQAaX2nij1cvzWp5RMRZ\nGhsbaWxsTPt1MxlcWontbK/F1EASpamx0hSmkBfgMeC5kQpgDy52sasiZ3fA3LUra/j297cy6A9w\n+GgHR5o7WFBbmdUyiYhz2L94A6xfvz4t183kk3ULpqO+DvACdwDPxKV5BviQtb0K6ADakuRdYst/\nO7B1tAWz11yyNc8lbEqpl5VXRF8ipo59EZkMMhlc/MA64AVgF/AksBu4x/qAqXUcwnTebwA+liQv\nwJeBHZihyPXAJ0dbMJ8vNzr0w+qvq4tsb3qlKWbVZhGRfJTJZjGA562P3Ya4/XWjyAvw3vEWyp9D\nfS4Al100hxnTSzlztpfOrgE2b2nl+mvmJ88oIpKjsv+1PQtiJ1Fm/xa43QW8s35RZP+FXx/MYmlE\nRMYv+0/WLBjM4svCRnLrjQspKDAD53btPcXRlnNZLpGIyNjlxpN1gmV7yf3hTJ9WwtVXRKfyvLhJ\ntRcRyV+ODC65MkM/3h/cfF5ku/GVI1rMUkTyVu48WSeQfTRWtoci2128bDZz50wFoLdvkJd+05Td\nAomIjJEjg0uudeiHFRS4+KN3RafxPPvCPgKBYIIcIiK5KXeerBMo9n0uuVNzAbj5+jqmlhUB0Haq\nh1e3tGa5RCIio+fQ4JIbC1cOp6jIw7tsfS/PPL+PUCjVZdpERHJDbj1ZJ0jM8i/e3Kq5ANx26+JI\n0Nt/6Ay79436rQIiIlnlzOCSwzUXgGmVxdS/oy6y/9TTu7JXGBGRMci9J2uGBYOhmE7yXOrQt3v3\nH54fmVS5fWebai8ikldy88maQTFzXDxuXK6xvE4m8+ZWTY1ZX+zJn+7MYmlEREbHccHF57OPFMvt\nP/99t18YCX6qvYhIPsntp2sGxMxxybFhyPGqq6Zyw7XR2stjP9qhkWMikhccGFzyp+YCpvYS7nvZ\nuecUb2w7nuUSiYgkl/tP1zTLxRWRE6mumsqtN0aX4//Bk9s1a19Ecl7uP13TzD4MOZfWFUvkjncv\np7jYvNetubVTa46JSM5zYHDJr5oLmHkvf7Lmgsj+kz/ZSW+fVkwWkdyVH0/XNLJ36OdLzQVg7eql\nTKssAaD9XB9P/ERDk0Ukd01EcFkN7AH2A/eNkOYh6/w2YEUKef8J2G2l/wlQkWphBnP0XS7JFBd7\n+PAHLonsP/vifg41tWexRCIiI8v009UNPIwJEhcCdwLL4tKsARYDS4C7gUdSyPsisBy4FNgH3J9q\ngWJWRPbkT80F4Ppr5nPxhbMBCIVCfOv7bxIMamiyiOSeTAeXlcABoAkYBJ4Abo9LsxbYaG1vBiqB\nqiR5fwkEbXlqUi1Qrq8rlojL5eJ/fehy3G5T7r0HzvDSbw5nuVQiIkNl+uk6D2i27bdYx1JJU51C\nXoA/B55LtUD2mksuroicTE11OX+y5vzI/sbHt3H6TG8WSyQiMpQnw9dPtc1mrAt8fQ7wAY8Nd7Kh\noSGyXV9fT319fV7XXMLeu3YZr2xupu1kN719gzz87dd54DM35Ow6aSKSuxobG2lsbEz7dTMdXFqB\nWtt+LaYGkihNjZWmMEnej2D6a24Z6Zfbg0tYvs3QH05RkYd7717J5/5hExBi+842Xtx0KOYlYyIi\nqQh/8Q5bv359Wq6b6afrFkxHfR3gBe4AnolL8wzwIWt7FdABtCXJuxr4NKYPpn80BYqZRJlnHfp2\ny5bO5Pbblkb2v/f4No6f6MpiiUREojIdXPzAOuAFYBfwJGYI8T3WB0x/ySFM5/0G4GNJ8gL8C1CG\n6djfCvxbqgXK16HIw7nzPRdRU10OwMCAn6/966sxqz6LiGTLZG6kDw23gvATP9nJU0+bCYjvXXsh\n//O9F010udLqYFM793/xpch7at5183nc85ErslwqEclXVt/tuGNDfn91HwP7y8K83vz/88+rm8ZH\n7rw0sv/Crw/yyqtHs1giEREHBhd7h74nj/tc7G67dTHXXBUd+/Bv39nCkeaOLJZIRJzOecElZlXk\nyfHnu1wu/vKjV1I1pwyA/n4/X/r6f9NxblRjHURE0mZyPF1HIXZV5MlRcwEoLSnksx9/ByXFhQCc\nOtPDVx76Hb5BdfCLyMRzYHCZPKPF4s2vqeATH7uacF/cnv2nefhbr2v9MRGZcJPr6ZoCe59LPi25\nn6orL6vmI3dGV09+5dWjfOc/tjLcyDkRkUxxYHCx1VzydPmXZP549dKY2frP/fIAT/50VxZLJCJO\nMzmfrgnYJxnm69piybhcLj76wRW84+roCLKnnt7Jz3+xL4ulEhEnmZxP1wRi57lMvmaxMLe7gHvv\nWcllF1dFjn33sbd4+rm9WSyViDiF44LLZFgVOVWFHjef+atruWDpzMix7z+xjf/8mZrIRCSzJvfT\ndRj2obmTsUM/XnGxhy986nqWXzArcuzxH7/N95/crlFkIpIxjgsu+f6ysLEoKS7k85+8nkuWz4kc\ne/rZPXzj3zdrHoyIZITjgov9YVrogJpLWFGRh7/9xHVcuaI6cuyVV4/y9//0G7q7fVksmYhMRo4L\nLrHLvzgnuICpqd1377Uxw5R37jnFpxt+pbXIRCStHBhc8v9NlOPhdhdw94cv54N3RCdatp3s5r71\nv+a3v9dqyiKSHo57ug74nNfnEs/lcvHuP7yAT627hqIi86Zrn8/PPz/yKo/+cKv6YURk3BwVXAKB\nIIHA5J+hn6prV9by4AO3UF01NXLs2Rf3c1/DSxxtOZfFkolIvnPU0zVmAmWhO/zGNUebX1PBgw23\nsPLyeZFjR5o7+PTf/Yr/enG/hiuLyJhMRHBZDewB9gP3jZDmIev8NmBFCnnfB+wEAsDlqRbEqSPF\nkplS6uW+j1/LRz+4IjLIYdAf4Ds/3Mrnv7SJ5tbOLJdQRPJNpoOLG3gYEyQuBO4ElsWlWQMsBpYA\ndwOPpJB3B/Bu4DejKYyTR4ol43K5WPPOJXxl/a3Uza+MHN+z7zR/8/kXeeInO9UXIyIpy3RwWQkc\nAJqAQeAJ4Pa4NGuBjdb2ZqASqEqSdw8w6lUYfTGd+Y5qEUzZ/JoK/vGBW3jPHy/D7Tb3KBAI8tTT\nO/mbz73I628e0/L9IpJUpp+w84Bm236LdSyVNNUp5B0VNYulxlvo5k/fdzFf/eKtLFk0I3L82Iku\nvvyNV2h48GUONbVnsYQikusyHVxS/Yo7IT3rTltXbLwW1FbypS/cxEc/uILSksLI8R27TvKpv/sV\n//Kt1zhxsjuLJRSRXOXJ8PVbgVrbfi2mBpIoTY2VpjCFvAk1NDREtuvr65lTfVFk34kTKMfC7S5g\nzTuX8I6ra3nq6V288OuD1giyEJt+28TL/32Em66r4z1rl1E1uyzbxRWRUWpsbKSxsTHt1810jcED\n7AVuAY4Br2E65nfb0qwB1lk/VwHfsH6mkncT8CngjWF+dyi+b2DbzjbWP/gyABctm80X768fz9/m\nSM2tnWx8fBtvbj8ec9ztLuCGa+dz+23nM7+mIkulE5HxsqZojDs2ZLrm4scEjhcwo78exQSHe6zz\nG4DnMIHlANAD3JUkL5iRYg8BM4Fnga3AbckKE9Ohr2axMamdV87nP3U9O/ec4smf7uTt3ScB0+m/\n6bdNbPptE5ddXMXa1Uu59KI5mksk4lCT+V/+kJrL715r5qsP/x6AVVfW8Jl7r81GuSaVt3ef5Kmn\nd0WCjF3tvAr+4KZF3HjtAsrKvFkonYiMVr7UXHKKz+GLVmbCRctmc9Gy2ezed5r/emEfr25pjQxV\nbm49x6M/3MoPntzONVfV8M6bFrFs6UzVZkQcwFHBRZMoM2fZ0pksWzqTEye7efbF/bz0m8P09/sB\nE9Rf/t0RXv7dEebOmcr119Ry3ar51FSXZ7nUIpIpk/kr5JBmsWdf3M+jP9wKwOpbFnP3h1NeOUZG\nqafXxyuvNvPLTYc4dGT4OTGLFkzjHatqecfVtcyeOWWCSygiw1Gz2BjYl9svcuhy+xNlSqmXd918\nHu+6+TwOHD7LrxoP89vfH6WvfzCS5tCRdg4daecHT26nbn4lV62o5qoV1Syqm0ZBwWT+3iMy+Tkq\nuISbaQAXCvE6AAAO5ElEQVSKix31p2fV4oXTWbxwOn/+p5ex5a1jvPJqM2+8dZxBfzTYNx3toOlo\nB//5s11MqyzhqhXVXH5JFcuXzWJKqQYDiOQbRz1h+weiwaVEwWXCeb1url1Zy7Ura+np9fHaG8d4\nZfNRduw6GfM6hPaOPl7cdJAXNx3E5XKxZNF0Llk+h0uWz2bp4hnqLxPJA456wtprLuE3MEp2TCn1\nctP1ddx0fR29fYNse7uN195s5Y23jtPd44ukC4VC7Dt4hn0Hz/CjZ3bh9Xo4f/F0zl88kwuWzOD8\nJTNUsxHJQY56wqrmkptKSwq55qoarrmqhkAgyN4DZ9jy1nF27DrJwcPt2Jeo8/n87Nh1kh27wvNq\nXMyvKeeCJTO5YOkMzqubTnVVWWRFZxHJDkc9Ye01l5LiwgQpJVvc7gIuPH8WF54/C4Cu7gHe3n2K\nbW+3sX1XGyfa4hfKDHG05RxHW87x4qaDAHi9HhYtqGRRXSWLFkzjvIXTqKkuV8ARmUCOCi72kUpF\nRWq3zwdTy4oitRqA02d62bP/tPnsO8Phox1D3i/j8/kjacK8hW5qqsuprSmndl4FtfPKqZ1XzuyZ\nUzQyTSQDHBVc+geio5PULJafZs4o5boZ87lu1XzAfGHYd/Ase/ad5uDhdg4ebqf9XN+QfL7BQGTo\ns11RkYea6nLm15Qzr2oqVXPKmDtnKnPnlGlEocg4OOpfjzr0J5+S4kIuXT6HS5fPiRw729HHoaZ2\nDjV1cPDwWQ4f7eD0md5h8w8M+Dl4+CwHD58dcm5aRQlzq8qYO8d8quaUMXvWFGbPnMLUMq+WsRFJ\nwFFPWHXoO8P0yhKmX1bClZdVR451dg1wtOUcza2d1uccR1s66eoeGPE67ef6aD/Xx669p4ac83o9\nzJpRyswZpcyeaX7OmlnKrBlTmDmjlOnTiin0qOlVnMtRT9i+vmifi5o8nKV8alFkkU27jnP9kWBz\nvK2b4ye6OXGym7ZTPQQCwRGuZvp1Wo930nq8c8Q0U8uKmFZZzPRpJeZnpfk5zfo5fVoJlRXFmrcj\nk5JjnrChUCim5qLgIgCVFcVUVhRz8YWxQScQCHLqTC/HT3RzvK2LEyd7OH6ii1Nnejl1ujdmcMhI\nuroH6Oo2NaZEyqZ4KZ9aRHl5EeVlRWZ7ahEV5UVMneqlYmqx9bOIivJivFq6SPKAY56wg4NB6/W8\nZrirmiwkEbe7gKrZZVTNLmMFVUPOd/f4OHW6l9Nnejl1podTp3sjgef0mV46OvuHjGIbSXePj+4e\nH8dOdKWUvqjIQ3lZEWVTCiktLaRsipeyKV6mlHopKys0P6d4mVJayJQpXsqmRI95PBqOLRPDMcFF\n/S2STuEH+sIFlcOeDwSCnOsc4GxHH+3t/eZnR7/px+no52y7+dlxLvUgFDYw4OfUgJ9TZ0Zf7qIi\nD6UlhZSWFFJS4qGk2ENJsdkuLSmkuNhDqbVfUlJozpcUUlocu19c5NEQbknIMU/ZHtuSIppAKZnm\ndhcwfVoJ06eVwMKR0wUCQbq6fXR2DcR8znWaJrXOLh+dnQOc6xqgs3OAzu6BhH1ByQwM+BkY8NPe\nMXS49mh5vR6KvG6KitwUeT2Rn8XFsfuR40VuvEXmp9n34C1yU+R14y001yksdFPoceP1FljbBRqV\nl6ccE1w6OqOjgioqirJYEpEot7sg0u+TilAoRG/fIF3dpimtp3eQbvt2j49e66f9WI+1PdpaUiI+\nnx+fz09X/KIJaeYtNEHH/CzA641u24/HpzFBKnqu0FOAx12A2+Oi0OPG4ymIfMLn7Mfs59zhc+4C\n1dhSlOngshr4BuAGvg08OEyah4DbgF7gI8DWJHmnA08CC4Am4P1AR7KCdHXZgstUBRfJTy6Xiyml\n3jEt1hkMhujrH6Svz09v3yB9/X76+wfp7fPT1zdIb9x+X5/fpO836fv7/fT1+enpG8Tn8yf/hWni\nGwzgGwzQM2G/MTG3u8Dqt40LRO74gOQyQcldQIHbFQlM4ePuAuu8x0VBQQEet/kZc94TTmcdc7tw\n29PY9z3Ra7oLXAl/p8vlwuNx4XK5ItcoKHClNXBmMri4gYeBW4FW4HXgGWC3Lc0aYDGwBLgaeARY\nlSTvZ4FfAl8B7rP2P5usMPaaS3m5s4JLY2Mj9fX12S5GTnDyvSgoiA1M47kXwWAIny/AgM9P/0Ag\n0tw24AuYz4CfgYHo+f5+P77B6PH+AT8+n/kZTjc4GGTQH8DnCzA4GMQ3GBhXE+BonDqxk1lVy1NK\nGwgECQSC+HzJ0zpZJoPLSuAApnYB8ARwO7HBZS2w0dreDFQCVZhW6pHyrgVutI5vBBpJIbicOh39\n3jO9smRUf0i+c/IDNZ7uRdR47kVBgYviYg/FxR4q0lusGIFAEL/fBJpwwPH5Agz6g/h8Aeu4dc4X\niAlOZjvIoFXz8VvXivlYx37d8ixLFt2A32/SDVrnA/4Qg/7YvJKaTAaXeUCzbb8FUztJlmYeUJ0g\n7xygzdpus/aH9aWvvxLZ3vLWsch2TXV5CsUXkWwLN0FlerkmV38jDQ23JE0XCoXiglTIBCS/LSDZ\nglMwGCIQDBHwB83PgAlowUAIfyBEMBg0PwO28/4QwWAoki4QPhc0wS4QDBIIH7edN9cMxuT1h9NY\neYMhky4QNOlMWlOO8H66ZPK/WKqlTKWRzzXC9UKJfo89oNiNNHxURCQRl8sMBpis8+RCoRAFP7gj\n28VIahXwC9v+/Zg+Ert/Bz5g29+DqYkkyrsHIrPa5lr7wzlANPjoo48++uiT2ucAOc4DHATqAC/w\nFrAsLs0a4DlrexXwagp5wx35YPpa/jHtJRcRkZx2G7AXEwnvt47dY33CHrbObwMuT5IXzFDkXwH7\ngBcxgwBERERERETyy2pMP8x+hvbxTEa1wCZgJ/A2cK91fDpmPtBwNbz7MfdnD/AHE1bSiePGTMb9\nubXv1HtRCfwIM4R/F2bEpVPvxf2YfyM7gMeAIpxzL76DGVm7w3ZsLH/7FdY19gPfzGB5c5Ib04xW\nBxQyfD/PZFMFXGZtl2GaEpdh+qY+Yx2/j2jf1IWY+1KIuU8HgMm2VO7fAP+BmXgLzr0XG4E/t7Y9\nQAXOvBd1wCFMQAGzwseHcc69uB5YQWxwGc3fHh7R+xpm/iKYvvLVGStxDrqG2FFmKc3en2Sexqxs\nEB55ByYAhUfVxY/a+wVmMMVkUYPpk7uJaM3FifeiAvNAjefEezEd86VrGibI/hx4J866F3XEBpfR\n/u1ziZ0A/wHMaN8R5XM0Hs5IkzKdog7zDWUzI082rcbcl7DJdo/+Gfg0YJ9K7cR7sRA4BXwXeBP4\nFjAFZ96Ls8DXgKPAMcxahL/EmfcibLR/e/zxVpLck8kWXELZLkAWlQE/Bj4OxL91Kjx+fSST5b79\nEXAS098y0uRcp9wLD2b05b9ZP3sYWot3yr04D/hrzJevasy/lT+LS+OUezGcZH/7mEy24NKK6eAO\nqyU22k5WhZjA8gNMsxiYbyP2yaYnre34e1RjHZsMrsWsPXcYeBy4GXNPnHgvWqzP69b+jzBB5gTO\nuxdXAr8DzgB+4CeYJnQn3ouw0fybaLGO18Qdn2z3JKFUJm5ONi7g+5jmILuRJpuGO+y8mKaTg6S2\nBE++uZFon4tT78VvgKXWdgPmPjjxXlyKGUlZgvmbNgJ/ibPuRR1DO/RH+7dvxow4dOHADn0YefLl\nZHUdpn/hLUxz0FbMf/REk03/FnN/9gDvmsjCTqAbiY4Wc+q9uBRTc9mG+bZegXPvxWeIDkXeiKnt\nO+VePI7pa/Jh+qTvYmx/e3go8gHMe7hERERERERERERERERERERERERERERERCQqCHzVtv8p4IE0\nXft7wHvSdK1E3odZDv+luOP1RCeBiuSdybb8iziLD3g3MMPaT+f6SOO5lmcUaf8C+Chwyzh+33D0\nb1uySv8DSj4bBP4f8Ilhzn2P2JpHt/WzHngZswbbQcyyFx/EvKtiO7DIludWzAz3vcAfWsfcwD9Z\n6bcBd9uu+1vgZ5iZ4PHutK6/g+hSG38HvAPzMqevxKUPYRZY/E/MUuc/tI7fDPzUlu6dmNn34b/x\nq5jVGq4BvmCVcwewwUrjsY7daO1/GfgHa/vPMEt8bMUsp15g/b3fs66xHbMApIjIpNYFTMUsVFkO\nfJJos9h3iQ0u4ZWi64F2zBLjXsziew3WuXuJrtH2Pcz6SQCLMctmFGGCyees40WY4FNnXbcbWDBM\nOauBI5galhvTBHa7dW4TZkHJePWYpeGrMWs5/Q6zMCeYYBOurT1GNPAFgffarjHNtv19zKrRYNaP\n2oUJnm9iAs4yzHI5bivNv2KC7uWY5UHCKoYpq8gQqrlIvuvCPDjvTZbQ5nXMqrA+zDpJL1jH38YE\nCjA1h6es7QOYF29dgHnt64cw3+5fxazRtNhK9xomiMS7ChNEzgABzFsyb7CdH2lRxNcwa0KFMLWR\nhdbxH2Ae/JWYFzk9bx0PYFbHDrvZKuN2a3u5dXwXpib0c8w6U35Ms9wVwBbrb7vF+n2HMLW5hzDr\nTHWOUFaRGKNpGxbJVd/AfAP/ru2Yn+iXpwJMLSVswLYdtO0HSfxvItwPsw7zsim7esw7U0bKZw8g\nLmL7dEbq37GXM2Ar23cxgaEfEwDDL0brt12rGFP7uAJTO3vAOhZ2MdEaXNhGzKKF8S7BLIb6v4H3\nY/qJRBJSzUUmg3bMQ/YviD5cmzAPVjDveCkc5TVdmJFcLszLphZhVol9AfgY0Qf9UqA0ybVex/Rx\nhJvFPoDp9xmr45gazeeJDah24UByBtN38z6i9+Z/YGo9NwL/gmnqegnTpDbLSjMdmG+V2YPp1/kC\nwzfhiQyhmovkM/s3/q9hahRh38J0rr+FeQ949wj54q8Xsm0fxTRNlQP3YJrRvo1pOnsTE3hOYkas\nJXqb33HMOzM2WXn+i+TDjIe7nn3/MWAmZrDBcOc7MPfgbcxLsTZbx2dgOvFvxtRoHga+CXwEE6xe\nxHzpHMQE0X5MAAt/EY1/m6WIiEwiD2P6S0RERNLiDaCR0Tf1iYiIiIiIiIiIiIiIiIiIiIiIiIiI\niIjksv8PtfEoZNyLxMwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa7d51b61d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hypos = range(1, 1000)\n",
    "suite = Hyrax2(hypos)\n",
    "\n",
    "data = 10, 10, 2\n",
    "suite.Update(data)\n",
    "\n",
    "thinkplot.Pdf(suite)\n",
    "thinkplot.Config(xlabel='Number of hyraxes', ylabel='PMF', legend=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior mean 185.570957948\n",
      "Maximum a posteriori estimate 50\n",
      "90% credible interval (36, 618)\n"
     ]
    }
   ],
   "source": [
    "print('Posterior mean', suite.mean())\n",
    "print('Maximum a posteriori estimate', suite.MaximumLikelihood())\n",
    "print('90% credible interval', suite.CredibleInterval(90))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we run the analysis again with a different prior (running from 0 to 1999), the MAP is the same, but the posterior mean and credible interval are substantially different:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior mean 233.983126212\n",
      "Maximum a posteriori estimate 50\n",
      "90% credible interval (36, 890)\n"
     ]
    }
   ],
   "source": [
    "hypos = range(1, 2000)\n",
    "suite = Hyrax2(hypos)\n",
    "data = 10, 10, 2\n",
    "suite.Update(data)\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": "markdown",
   "metadata": {},
   "source": [
    "This difference indicates that we don't have enough data to swamp the priors, so a more definitive answer would require either more data or a prior based on more background information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}