HotHands.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Hot Hot Hands"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "An exercise in Kruschke's [Doing Bayesian Data Analysis](http://amzn.com/0123814855) has you build a probability model for pairs of free throws, the second conditioned on the success of the first, and asks whether a \"hot hand\" phenomenon exists among NBA players. \n",
      "\n",
      "There's a [nice demo](https://www.wakari.io/sharing/bundle/wakari_demo/ipcluster_mcmc) running PyMC on IPCluster to solve this problem for many players at once, but why should they have all the fun? Here we work the original Larry Bird hot hands exercise, with some other bonuses.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Set up the environment for successive computational cells. \n",
      "# This stuff will not be repeated below, be sure to run this cell first.\n",
      "\n",
      "import pymc as mc\n",
      "import numpy as np\n",
      "from pymc import deterministic\n",
      "from IPython.display import Image\n",
      "\n",
      "class Struct:\n",
      "    # convenience 'bag of attributes' class\n",
      "    def __init__(self, **entries): self.__dict__.update(entries)\n",
      "        \n",
      "def show_graph(model):\n",
      "    \"\"\"display a graphical representation of the model\"\"\"\n",
      "    graph = mc.graph.graph(model)\n",
      "    Image(graph.create_png())\n",
      "\n",
      "def posterior_hist(trace, label=\"\", HDI=.95, alpha=1.0):   \n",
      "    \"\"\"Render a histogram of the posterior sampling distribution of a model variable.\"\"\"\n",
      "    t = sorted(trace)\n",
      "    h = hist(t, normed=True, alpha=alpha, label=\"%s:%.2f\" % (label, mean(t)))\n",
      "    ymax = max(h[0])\n",
      "    axvline(mean(t), c='red', linestyle='-.')\n",
      "    hdidx = int(len(t)*(1 - HDI)/2.) \n",
      "    annotate(\"\", xy=(t[hdidx],0), xytext=(t[-hdidx],0),\n",
      "             arrowprops=dict(arrowstyle=\"|-|\", connectionstyle=\"arc3\", color='k', lw=2))\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The text presents Larry Bird's free throw data as shot totals after an initial make or miss. We need to synthesize a list of throw events (success/fail) that match these proportions, to serve as the observed data for our Bernoulli likelihoods:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# free throw data for Larry Bird, 1980-1982. \n",
      "#\n",
      "throws = Struct(\n",
      "    hit_hit = [i < 251 for i in range(285)], # list of throw successes after an initial hit\n",
      "    hit_miss = [i < 48 for i in range(53)]   # list of throw successes after an initial miss\n",
      ")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "On to the model. We are estimating two parameters: the probabilities of making a second throw, after making or missing the first throw. We are given the make/miss counts, which can be treated as Bernoulli trial outcomes. To perform Bayesian estimation of the two probabilities, Kruschke suggests a Beta(30,10) prior for each, corresponding to the belief that NBA players make about 75% of their throws, $E[[Beta_{\\alpha,\\beta}]] = \\alpha/(\\alpha+\\beta)$.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def make_model():\n",
      "    \"\"\"\n",
      "    Return a dict populated with our model variables.\n",
      "    \"\"\"\n",
      "    # We need to feed a Model constructor a collection of named variables, and the \n",
      "    # prettiest way to do that is by exporting a dict of the local scope of a function, \n",
      "    # like we do here.\n",
      "\n",
      "    # Beta priors for our Bernoulli likelihoods. Assume NBA players make 75% of their shots.\n",
      "    phit_hit =  mc.Beta('phit_hit',  30, 10) # prob. of making second after making first\n",
      "    phit_miss = mc.Beta('phit_miss', 30, 10) # prob. of making second after missing first\n",
      "\n",
      "    # likelihood for observed second throws. Each throw is a Bernoulli trial.\n",
      "    hit_hit = mc.Container(\n",
      "                [mc.Bernoulli('hit_hit%d' % i, p=phit_hit, value=throw, observed=True) \n",
      "                for i, throw in enumerate(throws.hit_hit)])\n",
      "    hit_miss = mc.Container(\n",
      "                [mc.Bernoulli('hit_miss%d' % i, p=phit_miss, value=throw, observed=True)\n",
      "                for i, throw in enumerate(throws.hit_miss)])\n",
      "\n",
      "    # The statistic we care about is the difference between the hit/miss probabilities\n",
      "    # You could get a deterministic implicitly with 'diff = phit_hit - phit_miss' but it \n",
      "    # would default to having a strange name and trace=False. LinearCombination sets it \n",
      "    # up as a one-liner.\n",
      "    diff = mc.LinearCombination(\"diff\", [phit_hit, phit_miss], [1, -1])\n",
      "    \n",
      "    return locals() # give it up for the MC!\n",
      "\n",
      "# Instantiate and fit the model\n",
      "S = mc.MCMC(make_model(), db='ram')\n",
      "mc.MAP(S).fit() # load MAP estimates into our nodes as starting points\n",
      "S.sample(10000, burn=1000, thin=2)\n",
      "posterior_hist(S.phit_hit.trace(), \"hit_hit\", alpha=.5)\n",
      "posterior_hist(S.phit_miss.trace(), \"hit_miss\", alpha=.5)\n",
      "legend(loc='upper left')\n",
      "show()\n",
      "posterior_hist(S.diff.trace(), \"diff\", alpha=.5)\n",
      "legend(loc='upper left')\n",
      "show()\n",
      "S.summary()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " \r",
        "[****************100%******************]  10000 of 10000 complete"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEBCAYAAABysL6vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtwlNXdB/Dv5n4zCQaDxsQQICjehokWkEpQQU37ekHr\n66WI0o4ObW3ta5mhTDt13uo4RdAytb5aL2/r9c1oNaAgSdHUCaxcYk3xWkjTpVlMSDALuW6uu8/7\nR8yebLLJ7rN79jy7z/P9zDjm2ctzDj/O/vLw2/OcY9M0TQMREZlOgtEdICKi6GCCJyIyKSZ4IiKT\nYoInIjIpJngiIpNigiciMqmkYC/4/PPPceTIEXi9XixZsgSFhYXYsWMHPB4PAGDu3Lm46KKLot5R\nIiLSJ2iCd7lcuPnmmwEAu3btQmFhIdLS0nD11VdHvXNERBS+oAm+vLx80mMejwdVVVXQNA0lJSUo\nKyuLSueIiCh8QRP8mHfffReLFi0CAFRUVPger66ult8rIiKKWEhfstbV1aG0tBQzZ86c9FxycrL0\nThERUeSCXsHv3bsXs2bNwuzZs32PNTc3o7i4GAAwMDAw5Xt3796NxMTEyHtJRGQhubm5uOSSSyI+\nz7QJ/sSJE/jggw9w3nnn4fDhw+jp6cGaNWvQ3NyM+vp6AMDChQunfH9iYiLr80REOjU0NEg5z7QJ\nPj8/Hxs3bpz0eKAvXml6drsdl19+udHdiAmMhcBYCIyFfLzRiYjIpJjgFeGViWCVWKRt2hT0NbJj\nEUqbscoq40IlJngiIpMKeR68bL29vejq6oLNZjOqC0p1dXUhJyfH6G7EhKlioWkacnJykJWVZUCv\n5BsI8P3VRLLrzqG0GatYg5fPkATvcrkAAAUFBZZJ8AUFBUZ3IWZMFQtN03Dy5EkMDg4iLy9Pca+I\nzMeQEs3YB9gqyZ1CY7PZkJeXh8HBQaO7ogyvWAXGQj7W4ImITMqwGvxELvcQTrlHonb+GRlJyMtI\nidr5iSZK27QpaE1cdt05lDZjFWvw8sVMgj/lHsG+5q6onX9pcQ4TPCmVvnmz8mRrRJsUu1iiITIQ\nr1gFxkI+JvgJPvnkE9x0003Tvua+++7DsWPHJj3+5ptvoq2tTXp7stsERmesPProo3j88cfx61//\nGh0dHSG/t66uDps3b8bWrVvx4osvTnq+vb0dW7duxWOPPYbKykrdfTOL/g0bLNEmxa6YKdHEiosv\nvhhLliyZ9jUbN27EmWeeOenxoaEhDA8PS29PdpsAsHPnTlx55ZVYtGgRenp68MQTT+CXv/xlSO/d\ns2cPfvWrXwEAXn/9dfzrX//C3LlzAQDd3d2orKzE/fffb/mVRDkPXh/W4OVjgp/Ck08+CZvNBrfb\njRUrVqCsrAwDAwN47rnn0NDQgIceeghFRUW+1+/btw+7d+9GY2MjcnNzkZCQgHXr1iElJbS6f6D2\nAETc5ltvvYX169ejrq4OZ599tu/xw4cP4/rrrwcAnHbaabpiMzQ0BE3TYLPZ0NHRgfT0dN9zr732\nGsrLy/HYY48hPT0d99xzDzIyMnSdn4jkYIKfwre+9S3fVemjjz6KsrIypKWl4Sc/+UnAssPSpUvR\n3NyMZcuWobCwUEp7ACJu85JLLsH69esxa9Ysv8e9Xq/fcVpaWsh9vfLKK3Hrrbfi7LPPxsyZM/1u\nXGppacHx48fx4IMPoqurC3/4wx/ws5/9LORzWw2vWAXGQj4m+CmMJVsASEoKPUyapiltL1ibhYWF\n+OEPfzjp8ZER/ympod50NjQ0hLq6Ovz5z38GAFRWVuKTTz7BxRdfDABISUnB2rVrAYBLMxAZjF+y\nSpScnBxWPTyabba0tOCpp56a9JoLL7wQ+/fvBwD09PQEPMcbb7zh29hlzPDwsF/inj17tm/pCQC4\n5ppr8P777wMY/cUz8ReJlYSysqPdblfeZqySHQuKoSv4GRlJWFocvSu+GRmh/VEdDgcOHDiA1157\nDbfddhv27t2L/fv3w263Q9M0NDQ04NChQ7669w033ICSkhIAwFVXXYUnn3wSubm5sNlsyM3NxZo1\na8Jq74MPPsA3v/lN7N27N6I2P/zwQ/z2t7/FjTfe6FeDv/HGG7F582bs27cP3d3deOCBByb17Ykn\nnsCyZct8m60DQGZmJubPn4/HH38cSUlJSExMxI9//GPf85deein+9re/YcuWLRgYGMCdd94ZUtyJ\nSD6bFm5NIQS1tbUBt+xrbW3l4ls0JY4PsrqGhgasWLEi4vOwRENEZFIxU6Ixq1OnTuHll18O+EVo\nTk6O7wtJsibO/RYYC/mY4KNsxowZuP/++43uBpmEjEX5uPCedTDBE0VJKCs76r1iDbYo3+KX/gcH\n77pv2nPE6sJ7vHqXjzV4oihJ37xZeZtLXn5KeZsUu5jgiQzEud8CYyEfEzxRlBixsuOBNT9S3ibF\nrpipwfd2D6CvN3p7cWZmpSIrO/T1VogiFcrKjrLrzsHq77GMNXj5YibB9/UOoumLE1E7/7zz85ng\nichSWKKZQPWGH6Gaqs1oqaysxObNm/HII4/gww8/DPl9TqcTjzzyCJ544gls2rQp4Bo3b775Jq65\n5hqZ3Y1brDsLjIV8MXMFHytUb/gRqqnajIbjx4/D5XJhw9c15I0bN+Ib3/hGSO995ZVX8Itf/AI2\nmw1ffvklKisrcdddd/me7+jogMPhwMqVK6PSdyISmOCnoGrDj6effhq7du3C8uXL0d7ejvT0dLjd\nbvz0pz9FUVHRtG0CwNatW+H1epGamop///vf+O53v+tb/6eqqgpNTU1IT0+Hx+NBWloafvCDH/je\n++WXX+KKK67A1q1bfZt/AMChQ4dw5ZVX+o6LiorQ2dmJ3NzcoHFLTk5GT08PsrOz0dbWhuzsbL/n\nx9aH//3vfx/0XPEuGvPggwllHnysYg1ePib4Kaja8OO6665Deno61q5di4ceeggPPvggjh07hj17\n9mD16tXTtun1etHe3o6HH34YycnJOH78OFJTU33Pf/zxx3jggQeQm5sLr9eLI0eO+L0/Pz8f69ev\nx6WXXur3eGdnJy688EK/1/X09ISU4O+8807fL5nGxka/fu/evRvLli3jDk9EijDBT0Hlhh9jSTkz\nMxMAkJiYOGnHpUASEhJw33334dVXX8XIyAi6u7vx/e9/3/f8hg0b8Pbbb6Ovrw+9vb2TyiIpKSkB\nNwPJzs6Gy+Xy/Wuhq6sr5M07nn/+ebz++uvIyMhAY2MjXnjhBXzve98DMPovg9TUVBw6dAgHDhzA\nX/7yF1x77bUhnTceGbEna7xevQNciyYamOAlCmfDj0hWaz5y5AgyMzN9C5adOHECr732GtatW4fh\n4WF89NFHuOOOO3yvf/jhh/2uzIeHh/Hcc89h1apVfsvzXnrppaisrMTChQsBAM3NzZNKLfX19XA6\nnbjlllsm9WvsCr2goADd3d2+xzeMmxc+MDBg6uROFAtiJsFnZqVi3vn5UT1/KFRv+LFr1y40NDTg\nlltuQXd3N/bt24fZs2fj/fffx2233YaDBw9O2abT6YTdbseMGTOQkJCArq4u/OhHoze6DA4O4o9/\n/CP+/ve/w2azoa+vD9/+9rf92m5ra8Pjjz+OoqIivwQ/a9YsnHnmmdiyZQvcbje+853vTOr3tm3b\nsH///kkJ/qabbsJvfvMbZGZmore3128zkDHvvPMODhw4gJqaGlRUVIT092JWvGIVGAv5uOEHxRyO\nj6k1dbinXWwsFEuLczBvJr8HiWXc8IMoxhmxJ+vil/5H6vlU4jx4+WKmRGNW3PDDutI3bw7pi1aZ\nlrz8VFx/0UpyBU3wn3/+OY4cOQKv14slS5agsLAQbW1tqK2tRWJiIpYvX46zzjpLRV/jEjf8oOmw\n7iwwFvIFTfAulws333wzgNEvBAsLC1FfX4/Vq1cDALZv345Vq1ZFt5dEcYirSZLRgtbgy8vLJz02\n/kaVYHdqBpKamgqXyxXRFEEyH03T4HK5/G7WimehzoOXKZ7LM6zByxdyDf7dd9/FokWLAPjP3U5O\nTtbdaF5eHnp7e9Ha2gqbzab7/fFIz81CZjdVLDRNQ05ODrKysgzoFZH5hDSLpq6uDqWlpZg5cyYA\nwOPx+J4LlqDH/1a22+2+46ysLBw9ehQOhwMFBQUoKCiAw+Ew7fGCBQtiqj9GHi9YsCDg80ePHsWh\nQ4cCjhezHk/8rAR7/VGHw3fsdDrhdDrDPo6FP3+gGMRKf4w+liHoPPi9e/fijDPOwHnnned7bNu2\nbbjpppugaRp27tzpt1DVeFPNgyei8HAevDUomQd/4sQJfPDBBzh8+DC2b9+Ol19+GQCwePFiVFZW\norKyEosXL464E1bA+qJglVhwHrw+VhkXKk1bg8/Pz8fGAF8UFRQU+K1xQkREsceQpQqIKDws0VgD\nlyogIqJpMcErwvqiwFgIjIXAWMjHBE9EZFJM8IpwnQ3BKrEIZRZNNPZkjVdWGRcqMcETRUn65s3K\n21zy8lPK26TYxQSvCOuLAmMhMBYCYyEfEzxRlHA1STIa58ETxRHOg7cGzoMnIqJpccs+Rex2O2cJ\nfM2qsXC5h3DKPeL32FGHAyVz5oR8DvewJ/iL4pRVx0U0McETKXLKPTKpvOJs6cPxxNBLLvNZWiEd\nWKJRhFcmglViEco8+HPOOUdqm5wHT+MxwRMRmRQTvCKc4ytYJRah7Mk6foclGbgnK43HBE9EZFL8\nklUR1hcFxkKQXYMPxeCIF00d7rDfPyMjCXkZKRJ7NIrjQj4meCKL6Rn04KOWnrDfv7Q4JyoJnuRj\niUYR1hcFq8QilFk0smvw8TyLxirjQiUmeKIo4WqSZDQmeEVYXxQYC8GIGnys4riQjwmeKEq4miQZ\njQleEdYXBavEgvPg9bHKuFCJCZ6IyKSY4BVhfVFgLATW4AWOC/mY4ImITIoJXhHWFwWrxILz4PWx\nyrhQiQmeiMikmOAVYX1RsEosQplFI7sGH8+zaKwyLlRigiciMikmeEVYXxQYC0F2DT6ecVzIxwRP\nRGRSTPCKsL4oWCUW3JNVH6uMC5WY4ImihKtJktGY4BVhfVFgLATW4AWOC/mCJniv1wuPx6OiL0Sm\nwtUkyWjTbtlXU1MDh8OBVatWoaCgAACwY8cOX8KfO3cuLrroouj30gRYXxSsEgvOg9fHKuNCpWkT\nfEVFBb744gu/x9LS0nD11VdHtVNERBQ53TV4j8eDqqoqvPnmm2hoaIhGn0yJ9UWBsRBYgxc4LuSb\n9go+kIqKCt/P1dXVUjtDRETyRDSLJjk5WVY/TI/1RcEqseA8eH2sMi5U0p3gm5ubfT8PDAwEff34\nf3bZ7XYe89gyx06n0+/4qMPhV5JxOp26j12ujmmf7+rqiuj8eo6Njq/Zj2WwaZqmTfXke++9B4fD\ngbS0NMyePRvl5eXYs2cP2tvbAQALFy5EaWnplCevra1FWVmZ1A7HK7vdziuUr1k1Fk0dbuxr7vJ7\nzOl06rqKnz8zA40d7oj6Eek5lhbnYN7MjIj6EIhVx0UgDQ0NWLFiRcTnmbYGv3LlykmPlZeXR9wo\nERFFH+9kVYRXJgJjIXBPVoHjQj4meCIik2KCV0T2lyfxzCqx4J6s+lhlXKjEBE8UJVxNkozGBK8I\n64sCYyGwBi9wXMjHBE8UJVxNkozGBK8I64uCVWIRymqSsmvw8byapFXGhUpM8EREJsUErwjriwJj\nIbAGL3BcyMcET0RkUkzwirC+KFglFpwHr49VxoVKTPBERCbFBK8I64uCVWLBPVn1scq4UIkJnojI\npJjgFWF9UWAsBO7JKnBcyMcET0RkUkzwirC+KFglFtyTVR+rjAuVmOCJooSrSZLRmOAVYX1RYCwE\n1uAFjgv5mOCJooSrSZLRmOAVYX1RsEosOA9eH6uMC5WY4ImITIoJXhHWFwXGQmANXuC4kI8JnojI\npJjgFWF9UbBKLDgPXh+rjAuVmOCJiEyKCV4R1hcFq8SCe7LqY5VxoRITPBGRSTHBK8L6osBYCNyT\nVeC4kI8JnojIpJjgFWF9UbBKLLgnqz5WGRcqMcETRQlXkySjMcErwvqiwFgIrMELHBfyJRndASKz\nCrSaZF6iDRgaCfucCb0DyPN4pnz+k9U/mPr5lCS4PFrYbVP8YYJXxG638wrla1aJRcB58EMjaD7y\nle+wu6sL2Tk5IZ+zNyMZLvfwlM83L/pPYNz5xys+9wwgMTHktlSzyrhQiSUaIiKTCinBe71eeKb5\nZyEFxysTgbEQ9Fy9mx3HhXxBSzQ1NTVwOBxYtWoVCgoKAABtbW2ora1FYmIili9fjrPOOivqHSUi\nIn2CXsFXVFTgiiuu8Husvr4eq1evxu23346DBw9Gq2+mwjm+glViEco8+O6uLqltLqt+Uer5VLLK\nuFAprBp8RkaG7+eUlBRpnSEiInnCSvCaJqZaJScnS+uMmbG+KFglFqGsJim7Br/3W3dLPZ9KVhkX\nKoWV4Md/4Wqz2aR1hoiI5Akrwff39wMYvZIf+3kq4+tqdrvdssdjP8dKf4w8nhgTo/uj6viow4FT\nJ0/5jru7unC8tdXveHxNPtBxv9ut6/Xjj0+dPOW39o3T6YzoWHZ8nn766Zj6+zL6WAabNr7eEsB7\n770Hh8OBtLQ0zJ49G+Xl5WhtbUVdXR00TcPKlSuRn58f8L21tbUoKyuT2uF4ZbfzJo4xVo1FU4cb\nRxwnI7rRKS/IjU7TKT73DLgSEzF/ZgYaO9zB3zCFpcU5mDczI/gLdbLquAikoaEBK1asiPg8QadJ\nrly5ctJjBQUFuOOOOyJu3Eo4cAWrxCJt06agdXjZNfhl1S/GbR3eKuNCJd7JShQlRqwmWV7zkvI2\nKXYxwSsiu7YWzxgLQfY8+HjGcSEfEzxRlARaTTLa9lTcpbxNil1M8IqwvihYJRacB6+PVcaFSkzw\nREQmxQSvCOuLAmMhsAYvcFzIxw0/yFJ6uwfQ1ztoSNuDfcPI4I3fpBATvCKsLwpGxqKvdxBNX5xQ\n0tb8159D4633+o5PuoeRdUaW32s4D17gZ0Q+lmiIiEyKCV4R1hcFq8Ri/NX7VGTX4OP16h2wzrhQ\niQmeiMikWINXhPVFgbEQ4nFP1sERL5oiWKwMAGZkJCEvw3+zII4L+ZjgiUiXnkEPPmrpiegcS4tz\nJiV4ko8lGkVYXxSsEov5rz8X9DUq92Q9LcmGPI8HCb0DyPN4wv4vrPcnBp8fapVxoRKv4ImiZP4b\nz4f0RatM5TUvTflFq2fQg+YmF3ojWFMeQFjvLz73DCAxMew2KTy8gleE9UWBsRDisQYfLRwX8jHB\nE0VJ4y33KG+Tq0nSeEzwirC+KFglFpwHr49VxoVKTPBERCbFBK8I64sCYyGwBi9wXMjHBE9EZFJM\n8IqwvihYJRaxNg8+1lllXKjEBE9EZFJM8IqwvihYJRahzKLhnqyCVcaFSkzwREQmxQSvCOuLAmMh\ncE9WgeNCPiZ4IiKTYoJXhPVFIV5j0T/swUn3cMj/Fb36B7/jYa826ZzR2JM1XsXruIhlXE2SKET9\nw14c6xoM+fV3vfUn7LzqTt9xXkZyNLrlZ7rVJMl6eAWvCOuLAmMhsAYvcFzIxwRPFCVGrOzI1SRp\nPCZ4RVhfFKwSi1BKJZwHL1hlXKjEBE9EZFJM8IqwvigwFgJr8ALHhXycRUPKpaZkob3VmMQ2ODBi\nSLtERmCCV4T1ReGcwhI0fXHCkLZnnZ2trK1l1S8GrYlHYx58vNbh+RmRjyUaIiKTCivB79ixA9u3\nb8f27dvx6aefyu6TKbG+KDgcR43ughKhXElzT1aBnxH5wirRpKWl4eqrr5bdFyIikiisBO/xeFBV\nVQVN01BSUoKysjLZ/TId1heFOXOMq8HHGu7JKvAzIl9YCb6iosL3c3V1tbTOEBGRPBF/yZqcPP0C\nSuPrana73bLHYz/HSn+MPB6rwTudTjidTt/zKo5Pulxhv7/z1Cm/mnl3V9e0x4u2PeN33O92w93f\n7/f6462tIZ+vu6sL/W73tM8v2vbMlM+7+/t19V/28cR4ThwfTz/9dEyMz1g5lsGmadrkNUyDaG5u\nRnFxMQBg586duO666wK+rra2luWbr9ntdv4T9Gv77YfgHcg0pO1ZZ2ejvaU7rPeedA/rWk3ylz9d\ngUd+V+s7zstIRn5BNv7RJH7JdHd16SrT5GUkw+UeDrnN8RbMy8M/mlxBzxFpHwIpPvcMuBITfcdL\ni3Mwb2aG32v4GREaGhqwYsWKiM8TVommubkZ9fX1AICFCxdG3Akr4MAVWIMXWIMX+BmRL6wEX15e\nLrsfRKbD1STJaLzRSRHO8RU4D17gPHiBnxH5mOCJiEyKCV4R1heFOXNKjO5CzGANXuBnRD4meCIi\nk2KCV4T1RcEqNfhl1S8GfY3sGnwobcYqfkbkY4InIjIpJnhFWF8UrFKD556s+vAzIh8TPBGRSTHB\nK8L6omCVGnwouCerwM+IfNyyz6J6uwfQ1xv6uioyJSemwJiWiayFCV6RWKsv9vUOGrgvah7a+8Jb\n8CuecE9WfWLtM2IGLNEQRUl5zUuWaJNiFxO8IqwvCuPXZLc61uAFfkbkY4InihKuJklGY4JXhPVF\n4fS8PKO7oATnwevDz4h8TPBERCbFBK8I64sCa/ACa/ACPyPyMcETEZkU58ErwvqicHpeXtgbX4er\nf9iD/mEvsoY8OBnmhtPDXn3703MevHBakg0Y8fiO+1x9+PyU2+81M2adh8//2RHS+TKSE5CRkhj8\nhQAys1KRlZ0WemdNhAmeLKF/2ItjXYPIco/gWFd499HmZSRL7pV1eAY9aG4SpbnejGS4wvxFCwBF\nOak4PcS/j3nn51s2wbNEowjri4JVavDck1Uffh8hHxM8EZFJMcErwhq8YJV58KHgnqwCYyEfEzwR\nkUkxwSvCGrxglRo892TVhzV4+ZjgiaKEq0mS0ThNUpFANXgjN90YHBgxpF3AmHnwsYp1Z4GxkI8J\n3kDGbrqRbUi7VsLVJMloLNEowhq8YJUaPOfB68MavHxM8DHq5KkOaJq+W+MpPO6+LoyMhH9XJYVn\naGgA/e7eaV8z+nczpKhH5sMEr4ieefAff/oh7l73H6h6+5Uo9sg4sTQPfqC/F4/993fxf88/aEj7\nVq47P/+7/8LvfvM9eDyj3wdNjMXAQB8e+/VqvPqcMX83ZsAEH4OOt385+v+2Ywb3xPx6e05haGgA\nro4Wo7tiOa6OVvR2n8Tw0EDA5/t6OjE02M+/mwgwwSvCGrxglRo858Hrwxq8fEzwREQmxWmSinAt\nGkHvPPixtdwjoXctdxm4J6s+Vv4+IlosneCNvNEIMPZmo3gytpZ7JLiWu7l4vFrIG7d09A2jp8M9\n6fEZGUnIy0iR3bWYEnaCb2trQ21tLRITE7F8+XKcddZZMvulhMobjZxOJ8455xy/x6x6s9FoDZ4J\nFxitO/PKdZSeWAx6NLh6Qps+mXBqAK7uyb8MlhbnmD7Bh12Dr6+vx+rVq3H77bfj4MGDMvtEREQS\nhH0Fn5GR4fs5JSW834KaV0N3Z+ApUip4R9TVZSdevVuZVdai4Z6s+kTrXzIT94MdM3jKjfah6N/g\nZuSesGEn+PF3WSYnh/fPbY9XQ+NnxzE0NDn4KpwzN3ZuuAmkyXEYTY5/ROXcJ3uycPKr6e8ilME9\n5EHPhL/fszuy0dIaeoJPSUjAVxHs3wkAvWlJ6BoYQYqWg5YvxXS8zpPtAADXVy049OG7054jKzUJ\nvYOhf29SXPMSDs0s8Ht/dm4aWtvDj3uwPkxsc7wTziy0tvfq/nPo7cN0bY8/x+BAHwDg+JdNSEnL\nmPSek1+Nzn+P9C7jifvBjvHmpKJTwXczRu4JG3aC93jEh9Zms4V1DpsNyMvPgscT2QyJcCWE1+2w\nBKrBT+XTzz4CAPzr6BE8sHFtFHtF473xyia55wMAyeeMxTYj9b9Prp/2+e7OrxT1xHzCTvD9/f0A\nRq/kx36eKDc3Fw0NDeE2EXU9g0D66WraOvf0GQB6/B7r7u8J2H5RST6wb/Tn0tLS6HfOwoaGhtDc\n3AwASE1NNbg31jI4ODozas6cOUhMTJz0fEtLC9zu0dkvSy6YWE4dQuifjA4suSDQ4wNf/xddLW09\naGnT957c3Fwpbdu0MFe0am1tRV1dHTRNw8qVK5Gfny+lQwQ8++yzWLduHe699148++yzRnfH1Bob\nG3HuueeitLQUjY2NRnfHUnJyctDd3Y3Ozk7kBKi/NzU1obS0FHPnzkVTU5MBPYx/YV/BFxQU4I47\n7pDZFyIikohLFRARmRQTPBGRSTHBExGZVFg1eD3LFGzZssU3E6SjowP33HOP7nPEKhlx2LFjh2/K\n6dy5c3HRRRdFv+NRoCcWR48exUcffQSv14vLLrsMRUVFus8Ry2TEworj4sSJE/jrX/+K9PR0vz+v\nx+PBq6++aqlxMTEWc+bMARDGuNDC8NZbb/l+3rZt27Sv7e/v9/1cXV0d1jlilYw47N69e9Jrn3nm\nGQ2Adu+990ropRp6YrFr1y7fz++8805Y55DlyJEjGgCttLRU2jllxCLQuIhH08UiOztbA6B1dnZq\nmqZpNTU1vuf27t2r/fOf/9QAaGeeeeaU54gnesbFxFiM0TsuwrqC17NMQVpamu830viplDKWOjCa\njDh4PB5UVVVB0zSUlJSgrKwsOp2NMj2x8Hg88Hq90DTN745oM4wJQE4srDguUlNT0d/fj5SUFBw6\ndAhXXXUVACAhQVSSrTIuJsZi0aJFSElJ0T0uwkrwWhjLFHzxxRdYsmRJROeINTLiUFFR4fu5uroa\ngBjQ4wd2rNMTiwsvvBBbtmwBAKxZsyasc8gSjVjLiEWgcRGPpovFWMzH7oS//PLLUVNTg6GhIcyf\nP993I9T4O+XjNVcA+sbFxFj09fUhJSVF97gIK8GHs0zBwMCA7yo23HPEGhlxGG/sL/3aa6/F9ddf\n7/eBj3V6YvHxxx/j5z//ObxeL3bu3IkbbrhB9zlkmTNnDu6++25cdtll0s4pIxbjxXNSmy4WGzZs\nQHt7O7KM+WOCAAABTklEQVSzR5fNTkpKwnXXXQcA2L59OxYsWIC1a9ciMzNzynPEEz3jYmIsxsdg\nTCjjIqwEP9UyBZ999hkSEhJw/vnn+73e6/VOukIKZamDWCcjDs3NzSguLgYwmvwBoKioCG+//XY0\nuy6dnliM/fM0ISHBb5AaMSYSEhLwwgsvSD2njFgEGhfxaLpY3HjjjZM+IwDQ2dmJ3t5epKWl4U9/\n+hO2bdsW8BzxRm++AEQsxsaJ3nERVoJfvHgxKisrfcsUjNm/fz9sNtukjjqdTsybNy+kc8QTGXFo\nbm5GfX09AGDhwoXR73SU6IlFSUkJqqqqAAAXXHBB0HPEGxmxsOK4+PTTT3H48GFomoZbb7016Dni\njYxY6B0XYa9FQ0REsS1+vsUjIiJdmOCJiEyKCZ6IyKSY4ImITIoJnojIpJjgiYhMigmeiMikmOCJ\niEzq/wEJkvD2b31mNAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10ecfb890>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEBCAYAAACDu+UiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGltJREFUeJzt3W1sFNe9BvDH2Bh3IXgJjgEnBIeCcROiIJryJuQ0Ipek\nUpXSiN4blBClEmnURG0qtQmkVyJRP6S9bSPVTVoUJVJvw21dVcgl4SbkBUtYWVqwigW9GFFCF1jA\nNo5tXkyWF3u990Nj18b2eHfmzJ59Zp7fF7zr8czD+Pjv8X/PzilIp9NpiIhIoEywHUBERMxTcRcR\nCSAVdxGRAFJxFxEJIBV3EZEAGre49/f3I5VK5SKLiIgYUuT0yffeew/xeBxr1qxBRUUFAKC9vR0N\nDQ0oLCzEPffcg1mzZuUkqIiIZK5gvHnuhw8fRjQaHSzub7/9Nh588EEAwPbt27FmzRr/U4qISFay\n7rlHIpHBj4uLi42GERERM7Iu7kMv9CdOnGg0jIiImJF1cR/64mpBQYHRMCIiYobjC6qjuXz5MoB/\nXsEPfDyWDz74AIWFhe6SiYiEVDQaxRe/+EVP+3As7rt27UI8HkdJSQkqKytRU1ODpUuXoq6uDul0\nGvfdd5/jzgsLC7F48WJPAUVEwqa5udnzPhyL+2jFu6KiAuvWrfN8YAaxWAwrV660HcM15vzM2QHl\nt409vwl6h6qISACNO8/di4aGBrVlJK+U/OQnuLJpk+0YIo6am5uxatUqT/vQlbuISABlPVvGhHQ6\njY6Ojry/Z82FCxdQWlpqO4ZrXvIXFhaivLzc2nRXv3qmubpqZ+/5Kj8/K8W9o6MDN9xww7B3u+aj\ngVsusPKSP5lMoqOjAzNmzDCYSERyxUpbJpVK5X1hD7tIJGL1Lyv2qy7lt4s9vwnquYuIBJCVtsxo\nupLXcC7Z59v+p0WKMD2iG52x8KtnmqvZMuw9X+XnlzfF/VyyD38+ecG3/a+YU6riLvjcT3+qqZAS\nCnlT3PNRS0sLXnjhBWzbtm3wuaeffhqbNm3C7NmzR2z//vvv429/+xsikQgikQi++c1vDn7uH//4\nB37/+98jGo2iqKgIDz/8MKZNmzZuhrq6Opw6dQq9vb1YvXo1vvSlL2WUvbOzE7/85S8RjUaRSqXw\n7LPPDn6uvb0dv/3tb1FSUoJLly5hw4YNeffCKftVl/LbxZ7fBBV3B3fccQeWL18+7LlNmzZh5syZ\no25/4MABbNy4cdTPbd++HT/84Q+zupFaW1sburq68Nxzzw0eO9Pi/sYbb2Djxo2YPHky9u7di7fe\negtf+9rXAAA7duzAc889h4KCAvT29uLnP/85nn/++YxzMbv82bkUCTq9oHqdPXv24Ec/+hFefvll\nvPTSS/jkk08AAFeuXMErr7yCzZs3o729fdjXDHxu7969qK2tRW1tLU6cODH4+bq6OsRiMbz66quo\nra1FLBYb9vWnT5/GvHnzsGPHjmHPHzhwAPfee+/g49mzZ+P8+fMZ/T8KCwsxefJkAMCyZctw6NCh\nwc898cQTg/PXr1y5gpKSkoz2mUvXnyNTcjnPnZny89OV+xA9PT3YvXs3Nm/eDADo6+vD/fffDwAo\nKSnBd77zHdTV1Y34uoHPJZNJPPPMMyM+v27dOiQSiVE/BwDl5eX4/ve/j7vvvnvY8+fPn8fChQuH\nbdfT04NoNDru/6W/v39ExutdunQJv/jFL/C9731v3P2JCBcV9yGOHj06rO1RVFSEmpoa349bXFyM\nb3/72yOenzp1Krq6ugb7+9m847Svb/jMo+vfadrW1obXX38dP/jBD/LyPQeZ9kz9nmV1vUxnXbH3\nfJWfn4r7EFVVVaitrcXq1asBAL29vdi9ezdeeOEFX4/b29uL119/HWvWrBn2rtK7774bdXV1WLRo\nEQDg5MmTmDp16rCvbWpqQiKRwNq1a0fs9+LFi5g6dSr27t2LqqqqwecPHz48+BpAURH3EPB7ltX1\nNOtKWBS++OKLL/q18+PHj2PWrFkjnu/p6cENN9ww7LnuZC9OXbjqVxTMjpbgxojzmq+TJk1CQUEB\nfvOb3+DgwYPYtWsXpk6dijNnzuDTTz/F9u3bsWfPHhw/fhwHDx7EjBkzMG3aNHR2duKNN95ALBbD\n2bNn0dTUhAULFgy2Ql577TXs3r0bZ8+exb59+zBlyhSUl5cPHre1tRVPPvkk7rrrLixYsGDw+SlT\npuDEiRPYuXMndu3aha9//esjXsx99dVXUV9fj8cff3zY87fffjtqa2vR3NyMEydODJu5s3btWixa\ntAhNTU3Yt28fduzYMeod6Eb7PuVKLBbDrbfeOu522Y6bpW/+CmfuWuI6VybjCMg8f75Sfrva2tow\nd+5cT/uwcsvf1tbWEfc90ZuY8s9o36dcyfRNKMc6k1lduS9981fY99jTrnOtmFOKeWXjt7HY30Sj\n/HaZuOVv3vxNPj1SrOIrg/z6wfRS2LPBXFgA5Q8CTYUUEQkgFXfJS+zzlJXfLvb8JuRNW0aCw8Tr\nJ5eLozjWmRx3u2Rvfi/4ImKLtZWY0um0tVV+ZHwD3yM3jExPLJyOtgz2UZXBi5tDeX1BNVPsPV/l\n52elLVNaWoru7m4bh5YMdXd3Uy8xOJZlW39tO4JITli5cp8yZQquXr2K1tZWG4fPWJjXUJ00aRKm\nTJliOFHmEokE9Txl9ql4ys/PWs99+vTptg6dsXg8ji984Qu2Y7jGnt8Pe9c/ZTuCSE5YeROTBFu2\nbyzyoqosgqMZvPBqSqZvYhLxwsSbmDQVUkQkgFTcHbDPlWXOn0gkbEfwhPncA8ofBCruIiIBpOLu\ngP3Vdub8fs2UWfrmr3zZ7/WYzz2g/EGg4i4iEkAq7g7Y+3bM+f3quefqrpDM5x5Q/iBQcRcRCSAV\ndwfsfTvm/MzvTgW4zz2g/EGg4i4iEkAq7g7Y+3bM+f3quedqtgzzuQeUPwhU3CVUdFdICQtXNw47\nfvw49u/fj/7+fixfvhyzZ882nSsvsPftmPOr526X8vNzVdyPHDmCtWvXAgDefffdwBZ3CR7dFVLC\nwlVbJpVKob+/H6lUyvVqPQzY+3bM+TXP3S7l5+fqyn3hwoX42c9+BgBYv3690UAiIuKdqyv3gwcP\nYuPGjXj22Wfx17/+1XHbob9BY7EY1WPl9/Y4kUgMuwLP5vGtt96a0fZdXZ1Gjpfp4+PxeEb//5Ur\nV1r//nt5rPz2H3vlarGOnTt34itf+cqIj6+nxTrCSYt1iHhjbbGO2267DfX19aivr8fcuXM9Bchn\nJn+L2sCcX/Pc7VJ+fq567tXV1aiurjadRUREDNEaqmKc2jIi3mgNVRERGZWKuwP2vh1zfq2hapfy\n81NxFxEJIBV3B+z3p2DOrzVU7VJ+firuEiq6K6SEhYq7A/a+HXN+9dztUn5+Ku4SKrorpISF5rmL\ncZrnLuKN5rmLiMioVNwdsPftmPOr526X8vNTcRcRCSAVdwfsc2WZ82ueu13Kz0/FXUQkgFTcHbD3\n7Zjzaw1Vu5Sfn4q7iEgAqbg7YO/bMef3q+eeK8znHlD+IFBxFxEJIBV3B+x9O+b8WkPVLuXnp+Iu\noaK7QkpYqLg7YO/bMedXz90u5een4i6hortCSliouDtg79sx59c8d7uUn5+Ku4hIAKm4O2Dv2zHn\nV8/dLuXnp+IuIhJAKu4O2Pt2zPk1z90u5een4i4iEkBaQ1WM0xqqIt5oDVURERmVirsD9r4dc36t\noWqX8vNTcRcRCSAVdwfsc2WZ82sNVbuUn5+Ku4SK7gopYaHi7oC9b8ecXz13u5Sfn4q7hIruCilh\n4Xqeezwex/79+zFp0iQ88MADKC4uHrGN5rmHk+a5i3hjYp57kZsvOn/+PNra2vCNb3zD08FFRMQf\nrtoyBw4cwOTJk1FfX4+jR4+azpQ32Pt2zPnVc7dL+fm5unLv6OgAADz00EN45513MG/ePEyYoPa9\nBN/Vvn4cy6ANdLk4mtF245kWKcL0yMiWp8h4XBV3AKipqQEA3Hzzzeju7kZZWZmxUPmCfa4sc34/\n57l7WY2p52oK+8/0jL9h4XS0GXjdYcWcUivFnXnsAPz5TXB1uV1VVYV4PA4AOHfuHEpLS8fcduif\nR7FYTI9D9DiRSAxrr/jxuKurM6vtL1y4kNX2Xo/n9fHxz37OAPvfTz3O7WOvXM+W2bZtGwoLC3HT\nTTeN+VuSfbZMLBajvgKwld/EbJlEIpHR1XuuZ8tkerxM84/H1uwcjX27rM2WAYC1a9d6OrCIiPhH\nr4I6YP7ND3DnZ19DlT0/89gB+POboOIuIhJAKu4OTL64YQNzfvY1VDVP3y72/CaouEuo6K6QEhYq\n7g7Y+3bM+dl71uz5mccOwJ/fBBV3CRXdFVLCQsXdAXvfjjm/Xz1rL+9OzYZ67nax5zdBxV1EJIBU\n3B2w9+2Y87P3rNnzM48dgD+/CSruIiIBpOLugL1vx5xf89ztYh47AH9+E1TcRUQCSMXdAXvfjjm/\nXz3rXM2WUc/dLvb8Jqi4i4gEkIq7A/a+HXN+9p41e37msQPw5zdBxV1EJIBU3B2w9+2Y8/u5hmou\nqOduF3t+E1TcJVR0V0gJCxV3B+x9O+b87D1r9vzMYwfgz2+CiruEiu4KKWGh4u6AvW/HnF/z3O1i\nHjsAf34TVNxFRAJIxd0Be9+OOT97z5o9P/PYAfjzm6DiLiISQCruDtj7dsz5Nc/dLuaxA/DnN0HF\nXUQkgFTcHbD37Zjzaw1Vu5jHDsCf3wQVdxGRAFJxd8Det2POz96zZs/PPHYA/vwmqLiLiASQirsD\n9r4dc36toWoX89gB+POboOIuoaK7QkpYqLg7YO/bMedn71mz52ceOwB/fhNU3CVUdFdICYsi2wHy\nWSwWo74CGMjflbyGc8m+nB032ZvyvI9EIuHL1W8u57kzX70HZeyHmYp7CJxL9uHPJy/k7HhVZZGc\nHUtERqe2jAP23/zM+ZmvegH+/MxjB+DPb4Lr4n758mX8+Mc/Rmtrq8k8IiJigOvi3tjYiPvvv99k\nlrzDPleWOb/mudvFPHYA/vwmuCru586dQyQSQUlJiek8IiJigKviHpZXotn/j8z5tYaqXcxjB+DP\nb4Kr2TJnz57F22+/jY6ODlRWVqKiosJ0LhER8cBVcd+wYQMA4PDhw4hGo47bDr3KH+iDsTzesmUL\n7rzzzrzJ4zb/zOrFAP7VBx64qvTrcVVZtef9De1ZO20/vaAMQCRn/79Mj5dp/vEez0pNxryyhQBy\nO36G9qzzZTyHKb8JBel0Ou3mC7u7u7Ft2zZUVlZi9erVo27T0NCAxYsXewpoE3v7aSD/sc5kzue5\nH+1MetpHpm8CMnGsbGR6PFNvYloxpxTzLLxvIChjn1VzczNWrVrlaR+uZ8vceOON+Na3vjVmYQ8C\n5sEBcOfXGqp2MY8dgD+/CXoTk4SK7gopYaHi7oB9rixzfvZ54uz5mccOwJ/fBBV3CRXdFVLCQsXd\nAXvfjjm/5rnbxTx2AP78Jqi4i4gEkIq7A/a+HXN+9p41e37msQPw5zdBxV1EJIBU3B2w9+2Y82ue\nu13MYwfgz2+CiruISACpuDtg79sx5/erZ53LNVSZMY8dgD+/CSruIiIBpOLugL1vx5yfvWfNnp95\n7AD8+U1QcRcRCSAVdwfsfTvm/FpD1S7msQPw5zdBxV1CRXeFlLBwtRJTWLD37Zjzs/esTeW/2teP\nYzlcjGRapAjTI8XUYwfgHvumqLhLqLDdFbLnagr7z/Tk7Hgr5pRieqQ4Z8cT/6gt44C9b8ecX/Pc\n7WIeOwB/fhNU3EVEAkjF3QF73445v3rudjGPHYA/vwkq7iIiAaTi7oC9b8ecX/Pc7WIeOwB/fhNU\n3EVEAkjF3QF73445v9ZQtYt57AD8+U1QcRcRCSAVdwfsfTvm/Ow9a/b8zGMH4M9vgoq7iEgAqbg7\nYO/bMefXGqp2MY8dgD+/CSruEiq6K6SEhYq7A/a+HXN+9p41e37msQPw5zdBxV1Che2ukCJuqbg7\nYO/bMefXPHe7mMcOwJ/fBBV3EZEAUnF3wN63Y87P3rNmz888dgD+/CaouIuIBJCKuwP2vh1zfs1z\nt4t57AD8+U1QcRcRCSDXxb2lpQX19fXYtm0bTp8+bTJT3mDv2zHn1xqqdjGPHYA/vwlFbr+wq6sL\nDz30EADg3XffxS233GIslIiIeOP6yr2mpsZkjrzE3rdjzs/es2bPzzx2AP78JnjuuX/44YdYsmSJ\niSwiImKIp+Le2NiI+fPno6ysbMxthva+YrEY1eMtW7bkVR6v+ROJxLBecD4/Hvh4vO27ujqz2n/1\nKy95ypfp8TLNb/r/5/Xx8XgcwL/GUj6N52wes+c3oSCdTqfdfOFHH32Em266CdXV1WNu09DQgMWL\nF7sOZ1ssFqP+824g/7HOJP588kLOjltVFsHRzqSnfSQSiYxaG9ke65l/uwO1H7a4zpXp8TLNb+p4\npqyYU4p5ZZHAjH1Wzc3NWLVqlad9uHpBtaOjA3v27EF1dTWOHDmCnp4erF+/3lOQfMQ8OADu/Ow9\na/b8zGMH4M9vgqviXl5ejk2bNpnOIuI73RVSwkJvYnLAPleWOb/mudvFPHYA/vwmqLiLiASQirsD\n9r4dc372njV7fuaxA/DnN0HFXUQkgFTcHbD37Zjz+9WzztVdIdVzt4s9vwkq7iIiAaTi7oC9b8ec\nX2uo2sU8dgD+/CaouIuIBJCKuwP2vh1zfvaeNXt+5rED8Oc3QcVdRCSAVNwdsPftmPNrDVW7mMcO\nwJ/fBBV3CZVlW39tO4JITqi4O2Dv2zHnZ+9Zs+dnHjsAf34TVNwlVHRXSAkLFXcH7H075vya524X\n89gB+POboOIuIhJArhbrCAv2pbqY85taps4W1vxX+/pxrDOJ4/E4bps719djTYsUYXqk2Jd9M499\nU1TcLehKXsO5ZJ/vx7lcHMWxziSSvSnfjyXB0HM1hf1nepA48ynaCv1dd3fFnFLfiruouDvy6zf/\nuWRfbhasLpyOtpMXUFUW8f9Yhvk5zz0XfXfGq/ah2POH/aodUM9dRCSQVNwdsM+VZZ5rrTVU7WLP\nz/6za4KKu4hIAKm4O2Dv2zH3TZmzA8pvG/vPrgkq7iIiAaTi7oC9b8fcN9Uaqnax52f/2TVBxV1C\nRXeFlLBQcXfA3rdj7psyZweU3zb2n10TVNwlVHRXSAkLFXcH7H075r6p5rnbxZ6f/WfXBBX3ALvy\naQ96r121HUMC6mL3J0in07ZjyBhU3B0w9+16r13Ff//nE9iycYPtKK6w93yDnj9+qBmb/+PLeP9/\ntuQoUXaYf3ZNUXEPqGTPRSR7LqCzlfvPa8lPXW2nAQCdZzS+8pWKuwP17ezRPHe72PPrZ1e3/BUR\nSwYWBvHDwFoGQ/m5OEg+UnF3oL6dPVpD1a5c5B9YGMQXn61lMFTYFgdRW0ZEJIBcFff29nb87ne/\nwx/+8Ae0tbWZzpQ31Lezh73nq/x2sec3wVVxb2pqwiOPPIKHH34Y+/btM51JREQ8ctVzj0T+tSZn\ncXFwe1jqudujNVTtCmJ+P1/AHY3tF3BdFfeh70qbOHGisTC2nL10Dd3Jazk5VgGA/pwc6Z8udnfi\n48P/h1Pnr+TsmAUdJTk7XrbHWrn11zi17Ms5O55X+Xq8Ux+3AADS4HmHqq8v4I7C9gu4rop7KpUa\n/LigoMBYGFuu9fWjO9k34vkzp0/j5ltuMXqsAgDREv8nKaX6egc/fuaRB30/HouXAeDpf7cdIzD2\nN/wv1m/6L9sxRkgkEvR/fXjlqspcvnwZwD+v4Ac+Hk00GkVzc7O7ZDn2uVGemxcB0B03fqwrABb6\n/Duxs6hr8OP58+f7ezAJnY8//hgAUFZWhoUFLidVdPn3c7BwzkQA1+Xy8XijuZhoQ7PL13Wj0ajn\n4xekXdz5p7W1FY2NjUin07jvvvtQXl7uOYiY1dbWhoqKCsycOTPQM5rEjq1bt+Kxxx7Do48+iq1b\nt9qOI6NwdeVeUVGBdevWmc4iIiKG6E1MIiIBpOIuIhJAOS/u/f39w2bbsGHPL5Ip9rHOnt8rT3Py\n2tvb0dDQgMLCQtxzzz2YNWuW4/bvvfce4vE41qxZg4qKisHnd+zYMfhN+PznP48777zTS6yMmcqf\n7X5Myfa4Y22fy/OfTeaxtmU43/lwrjPNNJrxxvrevXtzEXkYE/lZzn9LSwv+/ve/o7+/H8uWLcMt\nn03Jzmrspz146623Bj/+05/+lNHXtLS0pM+cOTPsuQ8++MBLDNdM5XezHxOcjtva2poGkJ45c+a4\n2+fy/GdzrsbaNh/Pd6bb2hrr6XT2581prL/55ptpAOlHH33UbEgHJvKznP/GxsbBj9955x1X+/B0\n5W7qNgSpVAr19fVIp9O47bbbsHjxYi+xMmYqv63bMWR73LG2z+X5zybzWNsynO98ONeZZnK7j1wz\nkZ/l/NfU1Hjeh6finjZ0G4IHHnhg8OOdO3d6iZQVU/lN7cfkcSdMmDDsX6ftc3n+szlXY22bj+c7\n021tjXWnTG72Mdr48puJ/Gzn/8MPP8SSJUtc7SOj4t7Z2Yk//vGPw5776le/6sttCPz4YfU7v9+3\nY3CTv7y8HE8++SQWLFiQVU6/i2U252qsbW3d/sJE9qFyfV8mE+dtYB+rVq3C8uXL8fjjj5uIltWx\nATPf93w//42NjZg/fz7Kyspc7SOj4l5WVoannnpqxPP79+8HMPptCA4dOoQJEybg9ttvH3f/J0+e\nxJw5cwAAV66Yv0mS3/kzvR2DW27yt7S04Lvf/e6w/GPl9Pv8DzVWhtHO91jb+n2+x2Iiey7P9fWy\nyT/ePmbMmIHnn38e9957rz9hHY7tJT/L+f/oo48wY8YMVFZWZrSP0XhqyyxduhR1dXWDtyEY6i9/\n+QsKCgqGBd61axfi8ThKSkpQWVk52Fc6efIkmpqaAACLFi3yEslKfqf95FP+sbbP5fkfK0M2efPt\nfOfrub5eNvnzbaybys9w/js6OrBnzx5UV1fjyJEj6Onpwfr16x33MRpX95YREZH8pneoiogEkIq7\niEgAqbiLiASQiruISACpuIuIBJCKu4hIAKm4i4gEkIq7iEgA/T9yHzRdzzb4yAAAAABJRU5ErkJg\ngg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10d1b0f90>"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "diff:\n",
        " \n",
        "\tMean             SD               MC Error        95% HPD interval\n",
        "\t------------------------------------------------------------------\n",
        "\t0.022            0.042            0.001            [-0.054  0.106]\n",
        "\t\n",
        "\t\n",
        "\tPosterior quantiles:\n",
        "\t\n",
        "\t2.5             25              50              75             97.5\n",
        "\t |---------------|===============|===============|---------------|\n",
        "\t-0.051           -0.006          0.019          0.05          0.111\n",
        "\t\n",
        "\n",
        "phit_miss:\n",
        " \n",
        "\tMean             SD               MC Error        95% HPD interval\n",
        "\t------------------------------------------------------------------\n",
        "\t0.843            0.039            0.001            [ 0.765  0.911]\n",
        "\t\n",
        "\t\n",
        "\tPosterior quantiles:\n",
        "\t\n",
        "\t2.5             25              50              75             97.5\n",
        "\t |---------------|===============|===============|---------------|\n",
        "\t0.758            0.82            0.848          0.87          0.906\n",
        "\t\n",
        "\n",
        "phit_hit:\n",
        " \n",
        "\tMean             SD               MC Error        95% HPD interval\n",
        "\t------------------------------------------------------------------\n",
        "\t0.865            0.018            0.001            [ 0.829  0.9  ]\n",
        "\t\n",
        "\t\n",
        "\tPosterior quantiles:\n",
        "\t\n",
        "\t2.5             25              50              75             97.5\n",
        "\t |---------------|===============|===============|---------------|\n",
        "\t0.826            0.853           0.866          0.876         0.899\n",
        "\t\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice that a difference of 0 is within the 95% HDI, and thus it is credible to conclude no hot hand for Mr. Bird. But I'm bothered by using the same prior for hit/miss data having different numbers of observed hits and misses. To see this more clearly, consider a reparameterization of Beta in terms of its mean $\\mu$ and a \"confidence\" K:\n",
      "\n",
      "$Beta(a, b) = Beta(\\mu*K, (1-\\mu)*K)$\n",
      "\n",
      "K controls the prior mean's weight against new observations during updating. Since we use the same prior for both throw probabilities, we have built-in an initial bias towards these probabilities being the same. But because there are very different numbers of hit_hit vs hit_miss events, their posteriors feel the effect of an equivalent K differently. If 75% is the proper unconditional mean, then we will only be 1/3 as certain about prior hit_miss outcomes as hit_hit outcomes (because we only see 1/3 as many initial misses as initial hits): "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def make_model():\n",
      "    \"\"\"\n",
      "    Return a dict populated with our model variables.\n",
      "    \"\"\"\n",
      "    # Beta priors for our Bernoulli likelihoods. Assume NBA players make 75% of their shots.\n",
      "    # We reparameterize Beta in terms of mean (mu) and confidence (K):\n",
      "    # Beta(a, b) = Beta(mu*K, (1-mu)*K).\n",
      "    mu, K_hit, K_miss = (.75, 40, 13.3 )\n",
      "    phit_hit = mc.Beta('phit_hit',mu*K_hit, (1-mu)*K_hit) # prob. of making second after making first\n",
      "    phit_miss = mc.Beta('phit_miss', mu*K_miss, (1-mu)*K_miss) # prob. of making second after missing first\n",
      "\n",
      "    # likelihood for observed second throws. Each throw is a Bernoulli trial.\n",
      "    hit_hit = mc.Container(\n",
      "                [mc.Bernoulli('hit_hit%d' % i, p=phit_hit, value=throw, observed=True) \n",
      "                for i, throw in enumerate(throws.hit_hit)])\n",
      "    hit_miss = mc.Container(\n",
      "                [mc.Bernoulli('hit_miss%d' % i, p=phit_miss, value=throw, observed=True)\n",
      "                for i, throw in enumerate(throws.hit_miss)])\n",
      "\n",
      "    # The statistic we care about is the difference between the hit/miss probabilities\n",
      "    diff = mc.LinearCombination(\"diff\", [phit_hit, phit_miss], [1, -1])\n",
      "    \n",
      "    return locals()\n",
      "\n",
      "# Instantiate and fit the model\n",
      "S = mc.MCMC(make_model(), db='ram')\n",
      "mc.MAP(S).fit() # load MAP estimates into our nodes as starting points\n",
      "S.sample(10000, burn=1000, thin=2)\n",
      "posterior_hist(S.phit_hit.trace(), \"hit_hit\", alpha=.5)\n",
      "posterior_hist(S.phit_miss.trace(), \"hit_miss\", alpha=.5)\n",
      "legend()\n",
      "show()\n",
      "posterior_hist(S.diff.trace(), \"diff\", alpha=.5)\n",
      "legend()\n",
      "show()\n",
      "S.summary()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " \r",
        "[****************100%******************]  10000 of 10000 complete"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEBCAYAAABysL6vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1wVPW5B/DvZkmySZAEg0FCwou8XC8Fh0ENlEooBdrY\nQY1KdSzg1U5nqNI3x6ll7NSZ6jhVLGWwVm7HjlO1NlerAYGSFI0WCFCjpoiCgDR2FwJJzJL3zevu\nuX9g1izJ7p4957e75zz5fmYc9+yel9+z5+zDyXN+53ccmqZpICIicVKS3QAiIooPJngiIqGY4ImI\nhGKCJyISigmeiEgoJngiIqHGRJvh2LFjOHnyJAKBABYtWoSCggLs2rULfr8fADBjxgzMmzcv7g0l\nIqLYRE3wXq8Xt912GwBgz549KCgogMvlwsqVK+PeOCIiMi5qgi8uLh72nt/vR3l5OTRNw/Tp07Fg\nwYK4NI6IiIyLmuAHvfnmmygqKgIAlJSUBN+vqKhQ3yoiIjJN10XWffv2YdasWZgwYcKwz1JTU5U3\nioiIzIt6Bn/gwAFMnDgR06ZNC77ndrsxdepUAEBPT0/YZffu3Qun02m+lUREo0hOTg6uvfZa0+uJ\nmOCbmppw8OBBXH311Thx4gQ6Ojqwbt06uN1u1NTUAADmz58fdnmn08n6PBFRjGpra5WsJ2KCz8vL\nw8aNG4e9P9KF19GouroaN9xwQ7KbETeMz94kxyc5NpV4oxMRkVBM8CZIP4NgfNG5nnhC6XwqSd5/\nkmNTiQmeiEgo3f3gaTjpdUDGF13PCNeoLqVpGtz33gv/uXOmthWrtrY2ZGdnJ3SbiSIhNqfTiby8\nPDgcjrhtgwmeKM6amppw2WWXITMzM6Hbzc/PT+j2EklCbD6fD01NTZg4cWLctsESjQmSz24BxqeK\n3+9PeHIn68vMzAwO2hgvTPBEREKxRGMCa9T2piI+1xNP6KrDX8rr60OLb8DUtiMZnzkGuZlpcVs/\n2QMTPJEJGZs2GUrwLb4BHHK3xaFFFy2ems0ETyzRmCH57BZgfER2xwRPZEL3Qw8luwmGHT16FLfe\nemvEeTZs2IAzZ84Me//1119HQ0OD8u2p3iZwsZvqk08+ic2bN+NXv/oVmpubdS+7b98+bNq0CVu2\nbMELL7wQ8tnq1auxdetWbN26FZs3b8ZvfvObmNsWbyzRmMAatb0lqh+8VV1zzTVYtGhRxHk2btyI\nK6+8ctj7fX196O/vV7491dsEgN27d2PZsmUoKipCR0cHnn76afziF7/Qtez+/fvxy1/+EgDw6quv\n4t///jdmzJgBAHj++ecxbtw4AMA//vEPpKVZryTGBE80yj3zzDNwOBzw+XxYvnw5FixYgJ6eHjz3\n3HOora3Fo48+isLCwuD8hw4dwt69e3Hq1Cnk5OQgJSUF69ev153gRtoeANPbfOONN/Dggw9i3759\nmDx5cvD9EydO4KabbgIAXHbZZTF9N319fdA0DQ6HA83NzcjIyAh+NpjcAeDw4cMjDsyYbEzwJkg+\nuwUY32hx4403Bs9Kn3zySSxYsAAulws/+tGPUFZWNmz+xYsXw+12Y8mSJSgoKFCyPQCmt3nttdfi\nwQcfHHbjUCAQCJl2uVy627ps2TLccccdmDx5MiZMmDDiDVZ+vx8OhyOud6QaxQRPNMoNJlsAGDNG\nf0rQNC2h24u2zYKCAtx3333D3h8YCO2OqjcR9/X1Yd++ffjrX/8KACgrK8PRo0dxzTXXhMx36NAh\ny54s8CKrCdXV1cluQlwxvuiSMUqkFaSmphqqh8dzm/X19Xj22WeHzTN37lwcPnwYANDR0THiOl57\n7bXgQ4wG9ff3h4x3M23aNHi93mHLHjx4EIsXL44plkThGTxREozPHIPFU+M3WNb4zOg/7bq6Ovzz\nn//EK6+8gjvvvBMHDhzA4cOHUV1dDU3TUFtbiyNHjgTr3jfffDOmT58OAPjGN76BZ555Bjk5OXA4\nHMjJycG6desMbe/gwYP42te+hgMHDpja5nvvvYff/va3uOWWW0Jq8Lfccgs2bdqEQ4cOob29HQ88\n8MCwtj399NNYsmQJioqKgu9lZWVh9uzZ2Lx5M8aMGQOn04kf/vCHw5Z1uVxISbHmubJDM/p3lg5V\nVVV8ZB+NeufOnRMxOBapF+7YqK2txfLly02v35r/7BARkWks0ZjAfuL2Jj2+RGtpacFLL7004oXQ\n7Oxs3HPPPYlv1CjHBE9ESowfPx4//vGPk90MGoIlGhOkn/0xvuhGay8asgcmeCITMjZtSnYTiMJi\ngjeB/cTtTXp8REzwRCbYeTRJko8XWU1gjdreVMRndDTJzvYedHX2mt5+OFlj0zF2nP4xV0gmJnii\nJOjq7MXp401xW//MOXlM8MQSjRnSa7iMT7ZEP/BDr3DbjJeysjJs2rQJjz/+ON577z3dy3k8Hjz+\n+ON4+umn8cQTT4SMcfPTn/40+DCQrVu3jjjEQSLwDJ5olEr0Az/0CrfNeDh//jy8Xi8e+uJaysaN\nG3H99dfrWvbPf/4zHn74YTgcDpw9exZlZWW4++67AQCPPvpocLz4zs5OPP/88/EJIAomeBNYo7Y3\nVf3g7fxUJyBxD/zYtm0b9uzZg6VLl6KxsREZGRnw+Xz4yU9+gsLCwojbBIAtW7YgEAggPT0d//nP\nf/Dd7343ONZVeXk5Tp8+jYyMDPj9frhcLvzgBz8ILnv27Fl8/etfx5YtW4IP/wCAI0eOYNmyZcHp\nwsJCtLa2IicnJ+r3lpqaio6ODowbNw4NDQ0hDwAZ+vrvf/87vvWtb0VdXzwwwRONcol64MeqVauQ\nkZGBe+65B48++igeeeQRnDlzBvv378eaNWsibjMQCKCxsRGPPfYYUlNTcf78eaSnpwc///DDD/HA\nAw8gJycHgUAAJ0+eDFk+Ly8PDz74IK677rqQ91tbWzF37tyQ+To6OnQl+LVr1wb/kTl16tSI7QYu\nPlHq9ttvj7q+eGCCN0H6WCaMLzq7n70DiX3gx2BSzsrKAgA4nc5hT1waSUpKCjZs2ICXX34ZAwMD\naG9vx/e+973g5w899BB27tyJrq4udHZ2YsWKFSHLp6WljfgwkHHjxsHr9Qb/WmhrawsZAz6SP/7x\nj3j11VeRmZmJU6dO4U9/+hPuvffekHm6urowduxYXeuLByZ4IoqZkQd+mBmZ/OTJk8jKygoOWNbU\n1IRXXnkF69evR39/Pz744APcddddwfkfe+yxkDPz/v5+PPfccygtLQ0Znve6665DWVkZ5s+fDwBw\nu90h5RUAqKmpgcfjwerVq4e1KzMzEwCQn5+P9vb2YZ/v3bsX3/zmNw3HbRYTvAmSz24BxhdPWWPT\nMXNOXlzXH02iH/ixZ88e1NbWYvXq1Whvb8ehQ4cwbdo0vPPOO7jzzjvx7rvvht2mx+NBdXU1xo8f\nj5SUFLS1teH+++8HAPT29uL555/Hv/71LzgcDnR1deHb3/52yLYbGhqwefNmFBYWhiT4iRMn4sor\nr8RTTz0Fn883Yill+/btOHz48LAEf+utt+LXv/41srKy0NnZOWJPGbfbHbWnUjzxgR9EccYHflA4\nfOCHhUnvR834ouNokmRlLNEQmZCxaZOIC60q8IEf1hM1wR87dgwnT55EIBDAokWLUFBQgIaGBlRV\nVcHpdGLp0qWYNGlSItpqOaxR25v0+BKND/ywnqgJ3uv14rbbbgNw8SJJQUEBampqsGbNGgDAjh07\nUFpaGt9WElkUR5MkK4tagy8uLh723mDXIABR716TjDVqe1MRn57yjNPphM/nM70tksXn88HpdMZ1\nG7pr8G+++SaKiooAhPZnTU1NVd8qIovx+vrQ4huIebnxmWOQl5eHpqYmtLa2xqFl4cVy047dSIjN\n6XQiLy9+XWUBnQl+3759mDVrFiZMmAAA8Pv9wc8cDkfEZYfeLTh4xiRlmvHZezqW+Fp8A/i/Ax8B\nAKZMmQLg4miC0aYXTs7CyqK5mDhxYsLjq6urg9frtcz3rXI6Pz/fUu0xOv3pp59GPD7NitoP/sCB\nA7jiiitw9dVXB9/bvn07br31Vmiaht27d4cM3jMU+8GTFKebfTjkbot5ucVTszFzQmb0GYmGSEg/\n+KamJhw8eBAnTpzAjh078NJLLwEAFi5ciLKyMpSVlWHhwoWmG2FXrFHbm4r4Fr74e13zJaO/vOT9\nJzk2lSKWaPLy8rBxhItI+fn5IeM+EBGR9XCoAiIdWKKhROJQBUREFBETvAnS64CMz94kxyc5NpWY\n4ImIhGKCN0H6WCaMLzor96KRvP8kx6YSEzyRCYteelbXfBmbNsW5JUTDMcGbIL0OyPjsTXJ8kmNT\niQmeyIR/rrtf13wcdZKSgf3giXRgP3hKJPaDJyKiiJjgTZBeB2R89iY5PsmxqcQET0QkFBO8CdL7\n4jK+6NgPPjkkx6YSEzwRkVBM8CZIrwMyvujevXuDrvn0PLtVNcn7T3JsKjHBExEJxQRvgvQ6IOOz\nN8nxSY5NJSZ4IiKhmOBNkF4HZHzRWbkXjeT9Jzk2lZjgiUzgaJJkZUzwJkivAzI+e5Mcn+TYVGKC\nJzKBo0mSlTHBmyC9Dsj4oovWD753IIDTzT58/P0f43SzT/d/Xl+f6bZJ3n+SY1NpTLIbQCRZR68f\nH9R3xLzc4qnZyM1Mi0OLaDThGbwJ0uuAjM/eJMcnOTaVmOCJiIRigjdBeh2Q8UWntx+83vlUkrz/\nJMemEhM8EZFQfCYrkQ5Gn8k6e0ImTjX7Yl6Oz3Id3fhMViIiiogJ3gTpdUDGZ2+S45Mcm0pM8ERE\nQjHBmyC9Ly7ji87KvWgk7z/JsanEBE9kgt7RJPXOR6QSE7wJ0uuAUuPz+vpwutmHN2s+1j02jK/f\nn+xmx0zq/gNkx6ZS1LFoAoEANE2D0+lMRHuI4q7FN4BD7jZ46rtw3qmv6+PsMF0W9Y4mqXc+IpUi\n9oOvrKxEXV0dSktLkZ+fDwDYtWsX/P6LZzMzZszAvHnzwq6c/eDJioz0aTfan5394MkIVf3gI57B\nl5SU4Pjx4yHvuVwurFy50vSGiYgovmKuwfv9fpSXl+P1119HbW1tPNpkG9LrgNLj83g8yW5CXEne\nf5JjUynm8eBLSkqCrysqKpQ2hoiI1DHViyY1NVVVO2xJel9c6fFNmTLF9DrYDz45JMemUswJ3u12\nB1/39PREnX/on1LV1dWc5rRlpj0eT0iZxsh0W1tbxM+93mbD60/298Pp5E6rELEXzVtvvYW6ujq4\nXC5MmzYNxcXF2L9/PxobGwEA8+fPx6xZs8KuXHovmurqatFnElLjG+xF4/F4dJ/F27EXjdT9B8iO\nDUhQL5oVK1YMe6+4uNj0RomIKP54J6sJks8gAPnxqajBW5nk/Sc5NpWY4ImIhGKCN0H1BRGrkR6f\nin7wVu5FI3n/SY5NJSZ4IhM4miRZGRO8CdLrgNLjYw3eviTHphITPJEJHE2SrIwJ3gTpdUDp8amo\nwb979wal86kkef9Jjk0lJngiIqGY4E2QXgeUHh9r8PYlOTaVmOCJiIRigjdBeh1QenzsB29fkmNT\niQmeiEioiKNJmiV9NEmyJz6TlaxO1WiSPIMnIhKKCd4E6XVA6fHxmaz2JTk2lZjgiYiEYoI3QXpf\nXOnx8Zms9iU5NpWY4IlM4GiSZGVM8CZIrwNKj481ePuSHJtKTPBEJnA0SbIy9oOnUYf94Mnq2A+e\niIgiYoI3QXodUHp8rMHbl+TYVGKCJyISigneBOl9caXHx37w9iU5NpWY4ImIhGKCN0F6HVB6fHwm\nq31Jjk0lJngiIqGY4E2QXgeUHh+fyWpfkmNTiQmeiEgoJngTpNcBpcfHZ7Lal+TYVGKCJzKBo0mS\nlTHBmyC9Dig9Ptbg7UtybCoxwROZwNEkycqY4E2QXgeUHh/7wduX5NhUYoInIhJKV4IPBALw+/3x\nbovtSK8DSo+PNXj7khybSmOizVBZWYm6ujqUlpYiPz8fANDQ0ICqqio4nU4sXboUkyZNintDiYgo\nNlETfElJCY4fPx7yXk1NDdasWQMA2LFjB0pLS+PTOourrq4WfSYhPT6Px2P6LH7hi78Pqa/nOh1A\n30BwOqWzB7l+P6556VkcjeFCa2+LD419/SN+5nQ64PdHfxBbXd1nuOqq6bq3OVTW2HSMHecytGwi\nSD82VYma4EeSmfnlo8TS0tKUNYbI9voG4D75eXCyMzMVXl8/pnh9Ie9HE8hOR2tm6oifTZw8Do31\n7VHX4fF4EejJ0r3NoWbOybN0gid9DF1kHfoY19TUkQ/C0UD6GYT0+FTU4PX2jjlw4/+Y3lasJF9j\nkH5sqmIowQ+94OpwOJQ1hoiI1DGU4Lu7uwFcPJMffB3O0P6q1dXVoqa3bdtmqfYwvtim33///ZC+\n8B6Px/R0y4WW4HR7Wxu6fb6Q6fa2Nt3TZttjNr5k759I04OvrdKeeEyr4NCG1ltG8NZbb6Gurg4u\nlwvTpk1DcXExzp07h3379kHTNKxYsQJ5eXkjLltVVYUFCxYobbCVVFfLvtAjNb7TzT4ccrfFdJF1\n9oRMnGr2RZ0v1+8PqbXnflGDj1VhdjouN12DN34ReeacPEzMzza0bCJIPTYH1dbWYvny5abXE/Ui\n64oVK4a9l5+fj7vuusv0xu1O8gEGyI9P1TNZ9dThl1S8kPA6PGvwxDtZiUzQO0pkceWLcW4J0XBM\n8CaorpdZjfT4VIxFY2WS45N+bKrCBE9kgt5RIveX3B3nlhANxwRvgvQ6oPT42A/evqQfm6owwRMR\nCcUEb4L0OqD0+CTXqAHZ8Uk/NlUxNBYNkRV4fX1o8Q1En/ESvn4OfU2jAxO8CdLrgFaPr8U3gEPu\ntugzXmL2hIuD5bEfvH1Z/di0CpZoiIiEYoI3QXodUHp8iXwmazJ60bAGT0zwRERCMcGbIL0OKD0+\nyTVqQHZ80o9NVZjgiYiEYoI3QXodUHp8KmrUC1/8va75llS8ENN6/QENF3z9I/7X1ecf8f3uS7p/\nsgZPTPBEJsRrNMlev4Yzbb0j/nfBNzDi+939ASMhkGBM8CZIrwNKj09yjRqQHZ/0Y1MVJngiEzia\nJFkZE7wJ0uuA0uNjP3j7kn5sqsIET0QkFBO8CdLrgNLjk1yjBmTHJ/3YVIUJnohIKCZ4E6TXAaXH\nZ+V+8CqwBk9M8EREQjHBmyC9Dig9Pj6T1b6kH5uqMMETEQnFBG+C9Dqg9Pgk16gB2fFJPzZVYYIn\nIhKKz2Q1QXodUFp8uU4H0DeAlM4e5Pr9yJ08GfDrewD34DKXuualZ3F0yHAFmY6Rl+czWdWSdmzG\nCxM8jR59A3Cf/Bydmanw+vpjWjTcMmtf/l/sKvpOcPq/Z+aOuHxx5YtJudBKoxtLNCZIrwNKj6+9\nrS3ZTYgr1uCJCZ7IBL2jRHI0SUoGJngTpNcBpcc3Ljvb9Dr0ll3YD14t6cemKkzwRERCMcGbIL0O\nKD0+1uDtS/qxqQp70ZBldLb3oKuzV/f8vV39I3ZdDCdcF0YiqZjgTZBeB0x0fF2dvTh9vEn3/Bd8\n/TjTpv8fhEu7MKqowevt3263fvAD/X40novfXzhZY9MxdpzL8PLSf3uqMMET0TDdvn401nvjtv6Z\nc/JMJXjSx1ANfteuXdixYwd27NiBjz76SHWbbEN6HVB6fCpq8FbuRcMaPBk6g3e5XFi5cqXqthAR\nkUKGErzf70d5eTk0TcP06dOxYMEC1e2yBel1QOnxqajBWxn7wZOhBF9SUhJ8XVFRoawxRESkjul+\n8KmpqRE/H1orq66uFjW9bds2S7XH7vHV1X0WUjf2eDwRp1tbWkLq6O1tbRGnfd3dIdPnz52LafmR\npou2/yFk2tfdHTLd7fMBuNiLxsj6Y5lubWkJ+X7ef//9mL7PRE+bOV4GX1vp96J6WgWHpmlarAu5\n3W5MnToVALB7926sWrVqxPmqqqpEl2+qq6tF/6mY6Pgaz7XFvZvkJ6e9yP1iZMj2tjbdZZrcMKNJ\n/uIny/H41qph27h0uUvnM7q9kbYxqDA7HZdnfnnC5fF4DJdpJk4eh8b6dkPL6jFzTh4m5hsvkUn/\n7dXW1mL58uWm12OoRON2u1FTUwMAmD9/vulG2JXkAwxIXHxeXx9afAPo7erHhRiG8e0PxHxuEoI1\nePuS/ttTxVCCLy4uVt0OGsVafAM45G5Drt8f0xl5bmbk8mAicDRJsjKORWOC9L640uNjP3j7kn5s\nqsIET0QkFBO8CdLrgNLjYw3evqQfm6owwRMRCcUEb4L0OqD0+FTU4JdUvKB0PpVYgycmeCIioThc\nsAnS64DS45P2TFZ/QAu5j2DshEm67ivISE1BRqoznk1TTvqxqQoTPJEQvX4N3o6+mJcrzE63XYIn\nfViiMUF6HVB6fNKfySo5PunHpipM8EREQjHBmyC9Dig9PlXPZFU5n0qS+/lLPzZVYYInMqG48kWl\n8xGpxIusJkgfsnRofJ3tPejq1D8QWCx6u/qR6/cj0xGX1YcVy3DBdiQ5Pum/PVWY4EmXrs7emMZq\nj8XguO7/PTM3LuuPJ44mSVbGEo0J0s8gpMcnrR/8paSevQPyj01VmOCJiIRigjdBel9c6fFJ7icO\nyI5P+rGpChM8EZFQTPAmSK8DSo+P/eDtS/qxqQoTPBGRUEzwJkivA0qPT/ozWVmDJyZ4IiKhmOBN\nkF4HlB6f5Bo1IDs+6cemKkzwRERCMcGbIL0OKD0+6c9kZQ2emOCJTOBokmRlHGzMBOl1QCPxdff7\n0d0fiGmZ/oAW83ZUkFyjBmTHJ/23pwoTvADxHMp3UG/PgK75uvsDONMWW1tyM1ONNMkSOJokWRkT\nvAlWGZM6XkP5ejweTJkyBQAwcfI45etPNhXjpVu9H7zUs3ir/PasjjX4Uayl1Qt/wJ/sZgRpmoaO\nNm+ymyGGlb5Pf8CP1rYLyW7GqMMEb4KdzyAaGutxz/pV2Pbck2HnGTx7T5Sd5S/iyUfuwOkT7ydk\ne1LPbgcdO1KFJx+5A6eOv5vspuD3f/g1/mf9KjQ0nVOyPjv/9hKJCX6Uavz8PAJaAOcbzia7KUEN\n588AAC40n09yS2TwNtcDsMb3eb6xHoGAH59/nvy2jCaswZsgvQ44tAYvkYoa9ZKKF3TV1/XOp1Jf\nb19CtxeLgX4/Gs8Z76dfV/cZrrpqesR5ssamY+w4l+FtSMAET0QJ1+3rR2O98esDHo8XgZ6siPPM\nnJM36hM8SzQmSD57BxJfg0806c9kTUtPS/g2E0X6sakKz+CJRjl/QMMFX3/Ie2P7/MPeu1RGagoy\nUp3xbBqZZDjBNzQ0oKqqCk6nE0uXLsWkSZNUtssWWIO3N8n9xAH9NfhevwZvR+i8Y30DUW9YK8xO\nT1qCl35sqmI4wdfU1GDNmjUAgB07dqC0tFRZo6SJ952meu8y1WsgoOHzzj74UjJQ336x3emXDwRf\nh3NZGs/miKzEcILPzMwMvk5Lk1vri0Tv2Xu87jQdFI+7TL2+fvidGejtuvhn+hW9fjR3Rf6TPc3p\ngNPhUN6WeFH1TFar9qJhDZ4MJ3hN+3KAqNRUe44lEghoaGvxxXUbaWnWvsxxuu4ETtd9EvKePwDU\nX+jG0Htc07Rs1J+N3K2tLysVKQ6gsTPyPwSX6nSNQVvPANpaL97peOrEe0hLj977YWz6GHT26v/r\npckzFucaO2NeLtK2pla+iCMT8odt49LlLp3P6PZG2ka4Zf5z+kMAQGtLE+rPnAq7rcHvfyg9+3ug\nJRUXXPqO73/XndA1H6llOPv4/V/+/B02OmsbKhDQ0PJ5F/z+2EY/HHS2vh4FkydHnCcnN3JXrmRp\n+uKOQl93Fx7YeE9yG3OJEx8dxImPDia7Gbq8BgB/fkLZfPFQ/fYrqH77laRs+1KfNzcqWQ9r8PoY\nTvDd3d0ALp7JD76+VE5ODmpra41uwvKyJzjQ0Rv51uuOLz7OuDx+7Wjv7oh5/XkFGcHXs2bNUtwi\nY9xuN/r6Ll7sS09PT3Jr7K+39+I1k5ycHFxxxRVJbcunn34KALhicjoyLu8wdMwO9V+XjwfQEXGe\n+oYO1DcY30Yy5eTkqFmRZlB9fb32l7/8RXv55Ze1xsZGo6uhJKmqqtIAaMuWLUt2U4LWr1+vAdC2\nbduW7KaIsGHDBg2A9rvf/S7ZTdGWLl2qAdDeeeedZDdlVDF8Bp+fn4+77rpLzb8yRESkHO9kJSIS\nigmeiEgoJngiogQJBAIhPRDjzVANPpZhCp566qlgL43m5mZ8//vfj3kdiaQitl27dgV34owZMzBv\n3rz4N1ynwfg++eSTqPN+9tln+OCDDxAIBPDVr34VhYWFIeuw2r4DYmtbuPjssP/0xDfYu+3o0aOo\nq6vDVVddFfM6EimWdjU1NeHtt99GRkYG5s2bF4zNyvuusrISdXV1KC0tRX5+5Hsiwn0XMe87I1dm\n33jjjeDr7du3R5y3u7s7+LqiosLQOhJJRWx79+5V3zBFBuPT04tmz549wdd/+9vfhq1D09TuOxW9\naGJpW7j47LD/NC16fDfddFOwF82BAwcMrUMVPb1oYmlXZWVl8PXQ2Ky87zRN044dO6bV19dHnS/c\ndxHrvjN0Bh/LMAUu18U7EpuampCXl2doHYmkIja/34/y8nJomobp06djwYIF8WmsAUPji8bv9yMQ\nCEDTtJA7l62674DY2hYuPrvsv2jxOZ0XxwYKBAI4cuQIioqKkJaWZtn9F0u70tPT0d3djbS0tJDY\nrLzvYhHuu4h13xlK8JqBYQqOHz+ORYsWmVpHIqiIraSkJPi6oqJCXeMUGIwvJSUl5P8jmTt3Lp56\n6ikAwLp164atA1C77/S0KZpY2hYuPjvsPyB6fJO/uMv6448/xurVq9HV1YW0tLSk/PYG73aPtG9j\nadcNN9yAyspK9PX1Yfbs2cHYrLzvYhHuu4h13xlK8EaGKejp6Qme8RpdRyKoiG0oK/3jBXwZX1FR\nEZYsWYL77rsv7Lwffvghfv7znyMQCGD37t24+eabQ9YBqN13a9euxdmzZ3HjjTcaXkcsbQsX31BW\n3X9A9Pi0utHJAAABkklEQVTWrl0Lj8eDhx9+GEeOHEFWVlbM61Dl/vvvR25uLq6//vqw88TSrjFj\nxmDVqlUALo5mOxjbUFbbd7EI913Euu8MnSqFG6bg448/xvHjx4fNHwgEhv3LrWeog2RQEZvb7Q6+\n7unpiVNLjRmMKSMjAz/72c9w++23Axg5vsE/AVNSUkJ+LPHad4sXL8bOnTuDFzuNiGX/hYvPDvtP\nT3yLFi3Czp07kZOTg87OzmC8yfjtfec738Frr72GjIyMsPPE+tsDgNbW1pDYrLzvwhkpvnDfRaz7\nztAZ/MKFC1FWVgZN07BixYrg+4cPH4bD4cCcOXNC5vd4PJg5c6audSSbitjcbjdqamoAAPPnz49/\no2MQS3zTp09HeXk5AOArX/lK1HVYgYr4pOy/jz76CCdOnICmabjjjjuiriPZVMRm5X331ltvoa6u\nDi6XC9OmTUNxcTGAkeML913Euu8c2tCiDhERicEbnYiIhGKCJyISigmeiEgoJngiIqGY4ImIhGKC\nJyISigmeiEgoJngiIqH+H3IZKXuHQ5NOAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10f2f4d50>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEBCAYAAACDu+UiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGg1JREFUeJzt3WtsFNfZB/D/YuO4W4pNbAy4gA3B4BZQWrflJmSK4G1S\nKapQlKagyESVqpaGtvkQpUGtlEvVS0SVD9AWWqmv0oS2rioEJKS5gVPcLI1B2ApNjIhL12HjG46N\ncRYMju3d9wOvtzbePZ6dOTNznvH/9wXv7uzM3+Ph8fjZM2dCyWQyCSIiCpRpfgcgIiL9WNyJiAKI\nxZ2IKIBY3ImIAojFnYgogCYt7olEAiMjI15kISIiTXJVL7766quIRqPYsmULSktLAQBdXV2oq6tD\nTk4ONmzYgHnz5nkSlIiIrAtNNs793LlzKCwsTBX3F198EV/72tcAAEeOHMGWLVvcT0lERFnJuuce\nDodTX+fl5WkNQ0REemRd3Mee6E+fPl1rGCIi0iPr4j72w9VQKKQ1DBER6aH8QDWd69evA7h5Bj/6\ndSavv/46cnJy7CUjIpqiCgsL8YUvfMHROpTF/fjx44hGo8jPz0d5eTmqq6uxevVq1NbWIplMYvPm\nzcqV5+TkoKqqylFAIqKppqmpyfE6lMU9XfEuLS3Ftm3bHG9YgkgkgvXr1/sdwzbJ+SVnB5jfb9Lz\n68ArVImIAmjSce5O1NXVsS1DE+Q//TRu7NrldwwiYzU1NWHTpk2O1sEzdyKiAMp6tMxUIr1vZ2p+\nK2ftpma3yov8V69eRX9/vytDkvv7+1FQUKB9vV4xPX8ymURBQQFmzJjh2jZY3IkE6u3tBXBzgIMb\nxX10uhGpTM+fTCZx+fJlDA4OoqioyJVtsC2jIPnMEZCdX3J2wP38o0WBFxLKFAqFUFRUhMHBQde2\nweJORBRAbMsosO/rDiujZUzNbpUf+XsHPkbfwLBr658VzkVRmJMFSsHiPgW5XQRGZSoGn9i9m0Mh\nXdA3MIx/Xux3bf3rygpY3LM0PDyMUCjkyzQsLO4Kks8cgcz53S4Co5wUg6Du+6moubkZTzzxBA4e\nPJh6bufOndi1axcWLFgwYfnXXnsN//rXvxAOhxEOh/HNb34z9dp//vMf/PnPf0ZhYSFyc3OxdetW\nzJo1y1G+xsZGvPTSSwiHw1i4cCG+8Y1vWHpfMpnE7t27kZubi6tXr2Lnzp0oLi4GAFy4cAGHDx/G\nO++8g5/97Gdpv0+3sbiT567/8Id+RyAPLV++HGvXrh333K5duzB37ty0y7/99tt47LHH0r525MgR\n/OhHP9J6Jnz48GH89Kc/BQDs3bsXnZ2dlu4w99JLL2Hjxo1YtWoV4vE49u7dix//+McAgCVLluDR\nRx9FbW2ttpzZ4geqCpFIxO8Ijpia3+o4d8mk53fq5MmT+MlPfoJnnnkGP//5z/Hhhx8CAG7cuIFf\n/epXePzxx9HV1TXuPaOvNTQ0YM+ePdizZw/ef//91Ou1tbWIRCL49a9/jT179kzYx21tbViyZAmO\nHj1qOedHH3007qx6w4YNaGxstPTe8+fPY9WqVQCAT33qU5a36RWeuRORVvF4HCdOnMDjjz8O4Gbf\n+a677gIA5Ofn4/vf/37aM9rR1wYGBvDwww9PeH3btm2IxWJpXwOAkpISPPLII/jiF7847vnf/va3\naGhoGPfcjh07sGbNGvT392PmzJmp52fPno1z585Z+j4TicSE/CZhcVeQ3jeVnF9ydkB+fidaWlrw\npS99KfU4NzcX1dXVrm83Ly8P3/3udyc8v2PHDuzYsSPte2bOnJm6IAy4eWWr1R7+8PD4QQmmXXPA\ntgwRabV06VKcPn069XhoaAgnTpxwfbtDQ0PYt28fOjo6LL+noKAAbW1tqcdvvPEGPv/5z49bpr29\nHfv27cPQ0NC451esWIG33noLwM2/Vm593W85Tz755JNurby1tdXSBxOmikQiWLhwod8xbMuU//LA\nED7od+/KuFELCvNxe3jifXbzn34aw5Oc2QZ13+sSj8cn9Hnd/rlm+nne6rbbbkMoFMKzzz6Ls2fP\n4vjx45g5cyba29tx7do1HDlyBCdPnkRrayvOnj2LOXPmYNasWejp6cHvf/97RCIRXLp0CadPn8ay\nZctS7Y7f/e53OHHiBC5duoRTp05hxowZKCkpSW23o6MD3/nOd3DnnXdi2bJllr+vgoICPPvss2ho\naMD8+fNTffRRb7zxBp566incf//941o4lZWVqK2tRUNDA/7+97/je9/7Xipra2srnn/++dT3+c47\n7+DOO+9Ebu74Zkm6nyMAdHZ2YvHixZa/h3Q45a9CUC+kudAz4NlQyCXF4QnP8yIm5zo6OibMn8KL\nmORJ93ME9Ez5y567guTiApib38poGVOzW+VH/qJwHosvpbDnTkQUQCzuCtLHKkvOLzk7ID8/ycfi\nTkQUQCzuCuz7uiP/6acnXcbU7FZ5kf/Wi2hIFrd/fizu5LlP7N7tdwTxiouL0d7ezgIvVCKRQHt7\ne2qiMTdwtIwCh+P5R3J2wP38eXl5mDNnzoT5WXQx/R6kk5GQf86cOcjLc290E4s7eY6zQuqRl5fn\n2r1Co9EoPvOZz7iybi9Iz68DL2Kagvy+iImI1HRcxMSeOxFRALG4K0gfqyw5v+TsAPP7TXp+HVjc\niYgCiMVdQfJoDcDc/Bznbj7ml4/FnYgogFjcFaT37UzNz3uomo/55WNxJyIKII5zn4I4zp3IbBzn\nTkREabG4K0jv25ma38poGVOzW8X8/pKeXwcWd/IcZ4Ukcp+ticNaW1vR2NiIRCKBtWvXYsGCBbpz\nGUH6WFnJ+SVnB5jfb9Lz62CruJ8/fx733XcfAODll18ObHEnd3BWSCL32WrLjIyMIJFIYGRkBC4O\ntvGd9L6dqfk5zt18zC+frTP3FStW4Je//CUAoKamRmsgIiJyztaZ+9mzZ/HYY4/h0UcfxZkzZ5TL\njv0NGolERD0Oev5YLIZYLOba49Zo1Hb+9evX+77/nDxmfubX8f/XCVsXMb3yyiv46le/OuHrW/Ei\nJjPxIiYis+m4iMlWW2bRokU4dOgQAGD58uWOApgsEuF9PJ0YHE7gQs/AhOdL9z6Djh88onxvazSK\nRYsXW97WrHAuisLu3Y8yW37ve6eYXz5bxb2yshKVlZW6s1DAxAdH0Ngen/D86v5BnJrkL4dY+zV0\n5lj/62JdWYFRxZ3Ib5xbZgryqi2ztDiMljRn7m5gC4iChHPLEBFRWizuCjo/ufaD5PxjR91IJHnf\nA8wfBCzuREQBxOKuIP3TdlPzr37+N5Mus3DhQg+SuMfUfW8V88vH4k6eW3Ngn98RiAKPxV1Bet9O\ncn723P3F/PKxuJPnGmoe8jsCUeBxnPsUxHHuRGbjOHciIkqLxV1Bet9Ocn723P3F/PKxuBMRBRCL\nu4L0sbKm5uc4d/Mxv3ws7kREAcTiriC9b2dq/lPbd066DHvu/mJ++VjciYgCiMVdQXrfTnJ+9tz9\nxfzysbgTEQUQi7uC9L6dqfmtjJZhz91fzC8fizt5jrNCErmPxV1Bet9Ocn723P3F/PKxuJPnOCsk\nkftY3BWk9+1Mzc9x7uZjfvlY3ImIAojFXUF6305yfvbc/cX88rG4ExEFEIu7gvS+nan5Oc7dfMwv\nH4s7EVEA8R6qUxDvoUpkNt5DlYiI0mJxV5Det5Ocnz13fzG/fCzuREQBxOKuIH2srKn5eQ9V8zG/\nfCzu5DnOCknkvly/A5gsEomIPgOQnD8Wi2V19j44nMAFj0bmzArnoiicp1xG8r4HmD8IWNzJc27M\nChkfHEFje1z7etNZV1YwaXEn8pvt4h6NRtHY2IjbbrsNd999N/LygnewS//Nb2p+K7NCsufuL+aX\nz1Zxv3LlCjo7O/H1r39ddx4iItLA1geqb7/9Nj75yU/i0KFDaGlp0Z3JGNLHykrOz3Hu/mJ++WwV\n9+7ubly5cgX33nsv/v3vfyORSOjORUREDtgeClldXQ0A+PSnP43Lly9rC2QS6X07U/NznLv5mF8+\nW8V96dKliEajAIC+vj4UFBRkXHbsn0eRSISPDXoci8XGtT90P+7t7XF1/X4/9vvnx8fBfuyU7Vkh\nDx48iJycHMyePTvjb0nps0JGIrLHymbKL2FWyGzHuZs2A2VQjx0ppOfXMSuk7aGQ9913n6MNExGR\nezj9gILk3/yA7PzsufuL+eVjcSciCiAWdwWdH274wdT8vIeq+ZhfPhZ38hxnhSRyH4u7gvS+neT8\n7Ln7i/nlY3Enz7kxKyQRjccpfxW8HCvbO/Ax+gaGta6zNRrFosWLJzw/MDSidTvZsjIrZLbj3E0j\nfZw188vH4m6IvoFh7RcWxdqvoTNn4jqXTnIBDhHJx7aMgvTf/JLPfCVnB+QfO8wvH4s7EVEAsbgr\nSB8ra+pYcY5zNx/zy8fiTkQUQCzuCtL7dqb2rXkPVfMxv3ws7kREAcTiriC9bye5by05OyD/2GF+\n+VjciYgCiMVdQXrfztS+Ne+haj7ml4/FnTzHWSGJ3MfiriC9bye5by05OyD/2GF++VjcyXOcFZLI\nfSzuCtL7dqb2rTnO3XzMLx+LOxFRALG4K0jv20nuW0vODsg/dphfPhZ3IqIAYnFXkN63M7VvzXHu\n5mN++VjciYgCiMVdQXrfztS+tdV7qEom/dhhfvlY3ImIAojFXUF6305y31pydkD+scP88rG4ExEF\nEIu7gvS+nal9a95D1XzMLx+LO3mOs0ISuY/FXUF6305y31pydkD+scP88rG4k+c4KySR+1jcFaT3\n7UztW3Ocu/mYXz4WdyKiAGJxV5Det5Pct5acHZB/7DC/fLaL+/Xr1/GLX/wCHR0dOvMQEZEGtot7\nfX097rrrLp1ZjCO9b2dq35rj3M3H/PLZKu59fX0Ih8PIz8/XnYeIiDSwVdwjkciU6GlJ/x5N7Vvz\nHqrmY375cu286dKlS3jxxRfR3d2N8vJylJaW6s5FREQO2Cru3/rWtwAA586dQ2FhoXLZsWf5o30w\nKY/379+PlStXera90T7z6Fmr08dnzpxBSUnJhNeXFle6sr1bH/f29iAW67H1/rE9d7fy2X28rmwl\nAPXPc2zP15TjOZvHzO//Y6dCyWQyaeeNly9fxsGDB1FeXo6vfOUraZepq6tDVVWVo4B+8rL9dKFn\nAP+82K91nbFYLG17Y2lxGC09A1q3lY6T7WTK7sa2srWurABLisPKZaS3LpnfX01NTdi0aZOjddge\nLXP77bfj29/+dsbCHgSSDw7A3L4176FqPuaXjxcxkec4KySR+1jcFaSPlZU8VlxydkD+scP88rG4\nk+c4KySR+1jcFaT37UztW3Ocu/mYXz4WdyKiAGJxV5Det5Pct5acHZB/7DC/fCzuREQBxOKuIL1v\nZ2rfmuPczcf88rG4ExEFEIu7gvS+nal9a95D1XzMLx+LOxFRALG4K0jv20nuW0vODsg/dphfPhZ3\nIqIAYnFXkN63M7VvzXuomo/55WNxJ89xVkgi97G4K0jv20nuW0vODsg/dphfPhZ38hxnhSRyH4u7\ngvS+nal9a45zNx/zy8fiTkQUQCzuCtL7dpL71pKzA/KPHeaXj8WdiCiAWNwVpPftTO1bc5y7+Zhf\nPhZ3IqIAYnFXkN63M7VvzXuomo/55WNxJyIKIBZ3Bel9O8l9a8nZAfnHDvPLx+JORBRAuX4HMJn0\nvp2pfevVz/9m0r67qdkBYHA4gQs9A8pl5lZWTbqMFbPCuSgK5zleT7akH/vS8+vA4k6eW3Ngn6UP\nVU0VHxxBY3vck22tKyvwpbiTfGzLKEjv20nuW0vODsjPL/3Yl55fBxZ38hxnhSRyH4u7gvS+nal9\n66kwzl16funHvvT8OrC4ExEFEIu7gvS+neS+r+TsgPz80o996fl1YHEnIgogFncF6X07U/u+VmaF\nNDW7VdLzSz/2pefXgcWdiCiAbBf35uZmHDp0CAcPHkRbW5vOTMaQ3rczte87Fe6hKj2/9GNfen4d\nbF+h2tvbi3vvvRcA8PLLL2P+/PnaQhERkTO2z9yrq6t15jCS9L6d5L6v5OyA/PzSj33p+XVw3HM/\nduwYVq1apSMLERFp4qi419fXo6KiAsXFxRmXGdv7ikQioh7v37/f0+3FYrFxvVqnj8+cOaN1fdk+\n7u3tSfv66GgZ1ftHv/Yzv5PHuvK3RqOpx14e/6Nfm/T/cSrl1yGUTCaTdt745ptvYvbs2aisrMy4\nTF1dHaqqqmyH81skEvHsz7sLPQP458V+reuMxWJp2wNLi8No0TAd7WQybefh/1mOPceale/NlD3b\nbbnByrayzZ/JurICLCkOO15Ptrw89t0gPX9TUxM2bdrkaB22PlDt7u7GyZMnUVlZifPnzyMej6Om\npsZREBNJPjgA2X1fydkB+fmlH/vS8+tgq7iXlJRg165durMYp3fgY/QNDHuyrYGhEU+2YwLOCknk\nPt6sQ6Hp3RZ05hR5sq2lLvzpras1oJvVce4mZrdKen7pbQ3p+XXgFapERAHE4q6waPFivyM4IvnM\nUXJ2QH5+6We90vPrwOJORBRALO4KY8cYS2Tq/CZWZoU0NbtV0vNLn5tFen4dWNyJiAKIxV2BPXd3\n8B6q5pPes5aeXwcWdyKiAGJxV2DP3T+SswPy80vvWUvPrwOLOxFRALG4K7Dn7g7eQ9V80nvW0vPr\nwOJOnltzYJ/fEYgCj8VdgT13/0jODsjPL71nLT2/Dizu5DnOCknkPhZ3Bfbc3cFx7uaT3rOWnl8H\nTvlLZLDB4QQueHCHqVnhXBSF81zfDnmHxV2hNRoFPJrP3Q2S5xSXnB3Qlz8+OILG9riGRGrrygrG\nFXfp86FLz68D2zJERAHE4q7Anrs7OM7dfNLPeqXn14HFnYgogFjcFTjO3R1W76EqmfT80seJS8+v\nA4s7EVEAsbgrsOfuH8nZAfn5pfespefXgcWdiCiAWNwV2HN3B++haj7pPWvp+XVgcSfPcVZIIvex\nuCuw5+4fydkB+fml96yl59eBxZ08x1khidzH4q7Anrs7OM7dfNJ71tLz68Di7qJkMon+3g/9jkEB\nd+2jKxge+tjvGGQYFncFpz331/64H09s/TKi7zZpSpQdyX1fydkB7/JfvxbHUw9sxv8+8QOt65Xe\ns5aeXwcWdxf1dHxw89/OD3xOQkEV7+vFxzeuo6dDdhuI9BM3n3v/jWFcHRz2ZFvvX/wAwExPtuUG\nU+dEX/38bybtu5ua3Srp+aXPhy49vw7iintXfBCnYh+5vp1pIWBuTg4w4vqmiIi0E1fcvTR37lx0\neHAXHLeYeubIe6ia59bb+c2trHLt9n5e3NJvqp+1AyzuRATvbucHTLylH7nD1geqXV1d+NOf/oS/\n/OUv6Ozs1J3JGF1dXX5HcETyWGvJ2QHm9xvHudss7qdPn8YDDzyArVu34tSpU7ozERGRQ7baMuFw\nOPV1Xl5w/7xiz90dVkbLmJrdKubP7Nb+vhtGPzPwor9vKlvFPZlMpr6ePn26tjCB8//7qa3lHOaV\nLVEuGurOxwdXbniRyrNtZdrO+gP78MGaL3uyLTeYtK3utvcBAMPDQ57k0YH9fW/YKu4jI/8dHxgK\nhbSFsSI8PQd3FH3C9e2EAHR3dQKYYXsdZ+qOAgD+ceSP+MeRP+oJFgDPAMDO+/2OESh9lzq0rk/6\nOH3p+XWwVdyvX78O4OYZ/OjX6RQWFqKpSf+l9+6X9pvyARSG7J9hzJ49Gx9+eHNumYqKCk2piP6r\nra0t9X9wRcjB4IZeYMWY87QVZdMBuDRY4pZtuWE0/0exTjQJ/Gy4sLDQ8TpCybE9Fos6OjpQX1+P\nZDKJzZs3o6SkxHGQINq+fTsOHDiA5557Dtu3b/c7DgVQS0sLli1bhoqKCrS0tPgdhwxi68y9tLQU\n27Zt052FiIg04cRhREQBxOJORBRAnhf3RCIxbrSNNNLzE1kl/ViXnt8pR3PLdHV1oa6uDjk5Odiw\nYQPmzZunXP7VV19FNBrFli1bUFpamnr+6NGjqR/CHXfcgZUrVzqJZZmu/NmuR5dst5tpeS/3fzaZ\nMy0rYX+bsK+tZkrHtGM9222bVmuA7PI3NzfjvffeQyKRwJo1azB//vys14GkAy+88ELq68OHD1t6\nT3Nzc7K9vX3cc6+//rqTGLbpyp9pPTU1NUkAyeeee85h0vSyzZ9peS/3fzaZMy1r5+emg47suvf1\ne++9lwSQrKiomHTZbPdbNse6F3Tk96vWJJPZ5a+vr099/be//c3WOhydueuahmBkZASHDh1CMpnE\nokWLUFVV5SSWZbry+zUdQ7bbzbS8l/s/m8yZlpWwv03Y11Yzeb0Ou3RsW8r+r66udrwOR8U9qWka\ngrvvvjv19SuvvOIkUlZ05c+0nmnTpo37V7ds82da3sv9n03mTMvq+rllS0d23fs6m2NMx37za9/r\n2rZftQawl//YsWNYtWqVrXVYKu49PT3461//Ou65e+65x5VpCNw4YNzOn2k9Dz74IPr6+rB582bb\n6wb05beyvNv/YbPJnGlZv6a/0JF9LB37evHixXjwwQexdu3aSZfVsd/8nHpE97a9/uWUbf76+npU\nVFSguLjY1josFffi4mI89NBDE55vbGwEkH4agnfffRfTpk3DZz/72UnXf/HiRZSVlQEAbtzQPyGT\n2/kzTcewceNGbNy40Ul0APryZ8rp9v4fK1OGbPJanf5CNx3Zde/radOm4Q9/+IOlZbPJn+06vKAj\nv5fH+q2yyf/mm29izpw5KC8vt7SOdBy1ZVavXo3a2trUNARjvfXWWwiFQuMCHz9+HNFoFPn5+Sgv\nL0/1lS5evIjTp08DAD73uc85ieRLftV6TMqfaXkv93+mDNnkNW1/m7qvb5VNftOOdV35Jez/7u5u\nnDx5EpWVlTh//jzi8ThqamqU60jH1twyRERkNl6hSkQUQCzuREQBxOJORBRALO5ERAHE4k5EFEAs\n7kREAcTiTkQUQCzuREQB9H/71DixdY9howAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10ece4990>"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "phit_miss:\n",
        " \n",
        "\tMean             SD               MC Error        95% HPD interval\n",
        "\t------------------------------------------------------------------\n",
        "\t0.875            0.039            0.002            [ 0.8    0.941]\n",
        "\t\n",
        "\t\n",
        "\tPosterior quantiles:\n",
        "\t\n",
        "\t2.5             25              50              75             97.5\n",
        "\t |---------------|===============|===============|---------------|\n",
        "\t0.792            0.849           0.881          0.903         0.934\n",
        "\t\n",
        "\n",
        "phit_hit:\n",
        " \n",
        "\tMean             SD               MC Error        95% HPD interval\n",
        "\t------------------------------------------------------------------\n",
        "\t0.865            0.02             0.001            [ 0.826  0.904]\n",
        "\t\n",
        "\t\n",
        "\tPosterior quantiles:\n",
        "\t\n",
        "\t2.5             25              50              75             97.5\n",
        "\t |---------------|===============|===============|---------------|\n",
        "\t0.825            0.852           0.867          0.879         0.903\n",
        "\t\n",
        "\n",
        "diff:\n",
        " \n",
        "\tMean             SD               MC Error        95% HPD interval\n",
        "\t------------------------------------------------------------------\n",
        "\t-0.01            0.042            0.002            [-0.087  0.076]\n",
        "\t\n",
        "\t\n",
        "\tPosterior quantiles:\n",
        "\t\n",
        "\t2.5             25              50              75             97.5\n",
        "\t |---------------|===============|===============|---------------|\n",
        "\t-0.084           -0.037          -0.014         0.018         0.081\n",
        "\t\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Wha? The miss density is centered farther right than the hit density? Consider the original raw data: 48/53 is indeed greater than 251/285. The only reason our initial estimate showed a slightly hot hand is the unequal biasing towards the prior because of data size!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Continuing this example, recall that if you want to count successes in N Bernoulli(p) trials (like we do here for throws), that count will be distributed as Binomial(n,p). If we specify a Binomial likelihood for throw counts instead of a product of N Bernoulli likelihoods for individual throws, we should get the same estimate for the difference."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def make_model():\n",
      "    \"\"\"\n",
      "    Return a dict populated with our model variables.\n",
      "    \"\"\"\n",
      "    # Beta priors for our Bernoulli likelihoods. Assume NBA players make 75% of their shots.\n",
      "    mu, K_hit, K_miss = (.75, 40, 13.3 )\n",
      "    phit_hit = mc.Beta('phit_hit',mu*K_hit, (1-mu)*K_hit) # prob. of making second after making first\n",
      "    phit_miss = mc.Beta('phit_miss', mu*K_miss, (1-mu)*K_miss) # prob. of making second after missing first\n",
      "\n",
      "    # likelihood for observed number of second throw successs. Each is ~ Binomial \n",
      "    hit_hit = mc.Binomial('hit_hit', n = len(throws.hit_hit), p=phit_hit, value=sum(throws.hit_hit), observed=True)\n",
      "    hit_miss = mc.Binomial('hit_miss', n = len(throws.hit_miss), p=phit_miss, value=sum(throws.hit_miss), observed=True) \n",
      "\n",
      "    # The statistic we care about is the difference between the hit/miss probabilities\n",
      "    diff = mc.LinearCombination(\"diff\", [phit_hit, phit_miss], [1, -1])\n",
      "    \n",
      "    return locals()\n",
      "\n",
      "# Instantiate and fit the model\n",
      "S = mc.MCMC(make_model(), db='ram')\n",
      "mc.MAP(S).fit() # load MAP estimates into our nodes as starting points\n",
      "S.sample(10000, burn=1000, thin=2)\n",
      "posterior_hist(S.phit_hit.trace(), \"hit_hit\", alpha=.5)\n",
      "posterior_hist(S.phit_miss.trace(), \"hit_miss\", alpha=.5)\n",
      "legend(loc='upper left')\n",
      "show()\n",
      "posterior_hist(S.diff.trace(), \"diff\", alpha=.5)\n",
      "legend()\n",
      "show()\n",
      "S.summary()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " \r",
        "[****************100%******************]  10000 of 10000 complete"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEBCAYAAABysL6vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1wVOW9B/Dv5v1NEpoYJIAkvERLiWUiAlIMWtCmLZVo\nGd8Qr+1lxrZULdJrmd4ZZ4pjq6BldLzajh2vSm2K1YCCQIHUGwhQo8SIQAONi7saSNIsIW+b191z\n/0izZpPdkz27z+Y5z8n3M8OY3X32nC/Hw4/Db599jk3TNA1ERGQ5MbIDEBFRdLDAExFZFAs8EZFF\nscATEVkUCzwRkUWxwBMRWVTcaANOnTqFM2fOwOv1YtGiRZg6dSp27doFj8cDAJg5cyYKCgqiHpSI\niIwZtcC7XC7cfvvtAIA9e/Zg6tSpSEpKws033xz1cEREFL5RC3xRUdGI5zweD8rKyqBpGvLy8lBY\nWBiVcEREFL5RC/ygAwcOYMGCBQCA4uJi3/N79+4Vn4qIiCIW0oesFRUVmD17NrKyska8Fh8fLzwU\nERFFbtQr+MOHD2PSpEnIzc31PedwODB9+nQAQHd3d9D37t+/H7GxsZGnJCIaRzIyMnDttddGvB3d\nAt/U1IQjR47g6quvRm1tLdrb27FmzRo4HA5UVVUBAObNmxf0/bGxsezPExEZVF1dLWQ7ugU+Ozsb\nGzduHPF8oA9eraiyshJLliyRHSNszC8X88ujcnaR+EUnIiKLYoHXofoVAPPLZTR/0pNPRmVsuFQ+\n/ipnF4kFnojIokKeBy9aR0cHWltbYbPZZEUYVWtrK9LT02XHCJtq+TVNQ3p6OtLS0gCo30c1mr87\nwOddIsaGS+Xjr3J2kaQUeJfLBQDIyckxdYHPycmRHSEiquXXNA0XL15ET08PMjMzZcchUp6UFs3g\nH2AzF3caezabDZmZmejp6QGgfh+V+eVRObtI7METEVmUtB78cC53L1rc/VHb/sSUOGSmJERt+ySe\n6n1Uo/mTnnwy5N66kbHhUvn4q5xdJNMU+BZ3P446WqO2/cXT01ngydSSN28OuWgbGUvjF1s0ZFqq\nX4ExvzwqZxeJBX6YEydO4LbbbtMds27dOnz++ecjnn/rrbfQ0NAgfH+i9wkMzFh56qmn8Mwzz+BX\nv/oVmpubQ35vRUUFNm/ejK1bt+LVV18d8XpjYyO2bt2Kp59+GqWlpYazjVddjz4albE0frHAD3PN\nNddg0aJFumM2btyIK664YsTzvb296OvrE74/0fsEgN27d+Omm27Chg0b8Mgjj+D3v/99yO89dOgQ\nHn30Uaxfvx7Jycn49NNPfa+1tbWhtLQUDz30EH7+85/j7rvvNpxtUGVlZdjvNQOj+c04D15VKmcX\nyTQ9eLN5/vnnYbPZ4Ha7sWzZMhQWFqK7uxsvvfQSqqursWnTJkybNs03/ujRo9i/fz/Onj2LjIwM\nxMTE4IEHHkBCQmh9/0D7AxDxPt9++21s2LABFRUVmDJliu/52tpafO973wMAXHbZZYaOTW9vLzRN\ng81mQ3NzM5KTk32vbd++HUVFRXj66aeRnJyMtWvXIiUlxdD2iUgMFvggvv3tb2PmzJkAgKeeegqF\nhYVISkrCgw8+GLDtsHjxYjgcDtxwww2YOnWqkP0BiHif1157LTZs2IBJkyb5Pe/1ev0eJyUlhZz1\npptuwh133IEpU6YgKyvL7wtV9fX1uHDhAh577DG0trbid7/7HR555JGQtz2U6n1U5pdH5ewiscAH\nMVhsASAuLvTDpGnamO5vtH1OnToVP/7xj0c839/vPyU11C+d9fb2oqKiAn/5y18AAKWlpThx4gSu\nueYaAEBCQgLuv/9+AFBqmQQiK2IPXqD4+Piw+uHR3Gd9fT1eeOGFEWPmzp2LY8eOAQDa29sDbuPN\nN9/03dhlUF9fn1/hzs3N9S09AQC33HIL3nvvPQADf/EM/4vECNX7qEbzm201SZWPv8rZRTLNFfzE\nlDgsnh69K76JKaH9Vu12O/7+979j+/btuPPOO3H48GEcO3YMlZWV0DQN1dXVqKmp8fW9b731VuTl\n5QEAvvnNb+L5559HRkYGbDYbMjIysGbNmrD2d+TIEXzjG9/A4cOHI9rnBx98gN/+9rdYuXKlXw9+\n5cqV2Lx5M44ePYq2tjasX79+RLbnnnsON9xwg+9m6wCQmpqK/Px8PPPMM4iLi0NsbCx++tOf+l6f\nP38+PvzwQ2zZsgXd3d249957QzruRCSeTQu3pxCC8vLygLfsO3/+vHILYdHY4flB4111dTWWLVsW\n8XbYoiEisijTtGisqqWlBdu2bQv4QWh6errvA0kaSfX1RJhfHpWzi8QCH2UTJ07EQw89JDsGEY1D\nbNGQaal+BcZ7ssqjcnaRWOCJTCJ58+aojKXxiwWeTEv1uczML4/K2UVigScyCa4mSaKZZh58R1s3\nOjt6ohUFqWmJSJsQ+norJA/nwdN4J2oevGlm0XR29KDudFPUtj9rTjYLPBGNK2zRDDPWN/wIVbB9\nRktpaSk2b96MJ554Ah988EHI73M6nXjiiSfw3HPP4cknn/Rb4+ZnP/sZnn32Wd+voUscBKJ6H5X5\n5VE5u0imuYI3i7G+4Ueogu0zGi5cuACXy4VH/93n3bhxI6677rqQ3vvHP/4Rv/zlL2Gz2fDFF1+g\ntLQU9913HwBg06ZNmDBhAgCgo6MDL7/8cnR+A0QEgAU+qLG64ceLL76IPXv2YOnSpWhsbERycjLc\nbjcefvhhTJs2TXefALB161Z4vV4kJibis88+wz333OP73KOsrAx1dXVITk6Gx+NBUlISfvSjH/ne\n+8UXX+DGG2/E1q1bfTf/AICamhrcdNNNvsfTpk3DpUuXkJGRMepxi4+PR3t7OyZMmICGhgZfQQfg\n9/Nf//pXfOtb39LdlupzmcOZBx/qnZqMjA2Xysdf5ewiscAHMVY3/FixYgWSk5Nx//33Y9OmTXjs\nscfw+eef49ChQ1i9erXuPr1eLxobG/H4448jPj4eFy5cQGJiou/1jz/+GOvXr0dGRga8Xi/OnDnj\n9/7s7Gxs2LAB8+fP93v+0qVLmDt3rt+49vb2kAr8vffe6/tL5uzZs0HvyVpbW4vvf//7o26PiMLH\nAh/EWN7wY7Aop6amAgBiY2NH3HEpkJiYGKxbtw6vv/46+vv70dbWhh/+8Ie+1x999FG888476Ozs\nREdHB5YvX+73/oSEhIA3A5kwYQJcLpfvXwutra0h37zjD3/4A9544w2kpKTg7NmzeOWVV/CDH/zA\nb0xnZyfS0tJG3Zbq64kYzW/Ge7KqevxVzi4SC7xA4dzwI5JZqmfOnEFqaqpvwbKmpiZs374dDzzw\nAPr6+nD8+HG/m14//vjjflfmfX19eOmll1BSUuI3LXH+/PkoLS3FvHnzAAAOh8OvvQIAVVVVcDqd\nWLVq1Yhcg/dgzcnJQVtb24jX9+/fj1tuuSXs3zcRhcY0BT41LRGz5mRHdfuhGOsbfuzZswfV1dVY\ntWoV2tracPToUeTm5uK9997DnXfeiffffz/oPp1OJyorKzFx4kTExMSgtbUVP/nJTwAAPT09ePnl\nl/HRRx/BZrOhs7MT3/nOd/z23dDQgGeeeQbTpk3zK/CTJk3CFVdcgS1btsDtdgdspezYsQPHjh0b\nUeBvu+02/OY3v0Fqaio6OjoCzpRxOByjzlQC1O+jMr88KmcXyTRfdCIaxPODxjve8IMsT/W5zLwn\nqzwqZxfJNC0aq+INPyhUyZs3h/zhqZGxNH6NWuBPnTqFM2fOwOv1YtGiRZg6dSoaGhpQXl6O2NhY\nLF26FJMnTx6LrEriDT/Cp3oflfnlUTm7SKMWeJfLhdtvvx3AwAeCU6dORVVVFVavXg0A2LlzJ0pK\nSqKbkmgc4GqSJNqoPfiioqIRzw1OgwMw6jc1A0lMTITL5YpoiiBZj6ZpcLlcvu8FqN5HNZrfjPPg\nVaVydpFC7sEfOHAACxYsAOA/dzs+Pt7wTjMzM9HR0YHz58/DZrMZfv9YMfIFHzNSLb+maUhPTw/p\nS1BENLqQZtFUVFRg9uzZyMrKAgB4PB7fa6MV6KF/k1ZWVvoep6Wl4dy5c7Db7cjJyUFOTg7sdrup\nHrtcLlPlsXr+c+fOoaamJqTzR4XHw/MfOX4CB6pOoq7ZjbpmNw5UnRz18ZHjJ0yTX/bxNPJ4yZIl\npsoTzmMRRp0Hf/jwYVx++eW4+uqrfc/t2LEDt912GzRNw+7du/0Wqhoq2Dx4ovGortmNo45WQ+9Z\nPD0ds7JSRh9IljIm8+Cbmppw5MgR1NbWYufOndi2bRsAYOHChSgtLUVpaSkWLlwYcQizUr2Px/xy\nGc2/8LX/CXks58HrUzm7SLo9+OzsbGwM8GFOTk6O3xonRERkPlKWKiAaj9iioVBxqQIiItLFAq9D\n9T4e88vF/PKonF0kFngiIotigdeh+noWzC+X0fxmm0Wj8vFXObtIXE2SyCQWbXsB79+3zu+5nn4v\n6prdI8Zet3kzTq4duYjdxJQ4ZKYYXz6ErIlX8DpU7+Mxv1wi8rf3eHDU0TriF4CAz7e4+yPe5yCV\nj7/K2UVigScyib+v+UlUxtL4xXnwRGMknHnw+VkpOBugRRMM581bA+fBExGRLhZ4Har38ZhfLuaX\nR+XsIrHAExFZFAu8DtXn0jK/XNGcB29kbLhUPv4qZxeJBZ6IyKJY4HWo3sdjfrmM5h/+JSdRY8Ol\n8vFXObtILPBERBbFAq9D9T4e88vF/PKonF0kFngiIotigdeheh+P+eWK5j1Zx2IWjcrHX+XsIrHA\nE5nEom0vRGUsjV8s8DpU7+Mxv1zML4/K2UVigScyCa4mSaKxwOtQvY/H/HJxHrw8KmcXiQWeiMii\nWOB1qN7HY365mF8elbOLxAJPRGRRLPA6VO/jMb9cnAcvj8rZRWKBJyKyKN6TlWiM8J6sFCrek5WI\niHSxwOtQvY/H/HIxvzwqZxeJBZ6IyKJY4HWoPpeW+eXiPVnlUTm7SCzwRCbB1SRJNBZ4Har38Zhf\nLuaXR+XsIo1a4L1eLzwez1hkIRrXuJokiaY7D37fvn2w2+0oKSlBTk4OAGDXrl2+gj9z5kwUFBQE\n3TjnwRN9ifPgKVSi5sHH6b1YXFyM06dP+z2XlJSEm2++OeIdExFRdBnuwXs8HpSVleGtt95CdXV1\nNDKZhup9POaXi/nlUTm7SLpX8IEUFxf7ft67d6/QMEREJE5Es2ji4+NF5TAl1efSMr9cnAcvj8rZ\nRTJc4B0Oh+/n7u7uUccP/adSZWUlH/PxuH18zm6H0+n0PXY6nX6PW1tbR7zucjUHHR/o8Tm73TS/\nXz6O7LEIurNoDh48CLvdjqSkJOTm5qKoqAiHDh1CY2MjAGDevHmYPXt20I2rPoumsrJS6SsB5pdr\neH7VZtGofPxVzg6M0Sya5cuXj3iuqKgo4p0SEVH08ZusOlS+AgCYXzbml0fl7CKxwBMRWRQLvA7R\nH3iMNeaPLpe7F3XN7qC/DlSd9Hvs7tNf8sNss2jMfvz1qJxdJMPz4IloQIu7X/dDU2d9Jy7Efvl6\n/igffi7a9gLev29dSPs2MpbGL17B61C9j8f8cl155ZWyI0RE5eOvcnaRWOCJTIKrSZJoLPA6VO/j\nMb9cQ7+EFAojLZexaM+ofPxVzi4SCzwRkUWxwOtQvY/H/HKxBy+PytlFYoEnIrIoFngdqvfxmF8u\noz14zoMXR+XsIrHAExFZlO5qkpFSfTVJIj1GV4c0ujJkOO/hPVmtQdRqkryCJyKyKBZ4Har38Zhf\nLqM9eLNR+firnF0kFngiIotigdeh+lxa5pfL6Dx4s82iUfn4q5xdJBZ4IpNYtO2FqIyl8YsFXofq\nfTzml4s9eHlUzi4SCzyRSXA1SRKN8+CJwsR58BQtnAdPRES6WOB1qN7HY3652IOXR+XsIrHAExFZ\nFAu8DtXn0jK/XJwHL4/K2UVigScisigWeB2q9/GYXy7ek1UelbOLxAJPRGRRLPA6VO/jMb9cvCer\nPCpnF4kFnojIoljgdajex2N+uXhPVnlUzi4SCzyRSXA1SRKNBV6H6n085peLPXh5VM4uEgs8kUlw\nNUkSjQVeh+p9POaXi/Pg5VE5u0gs8EREFhVSgfd6vfB4PNHOYjqq9/GYXy724OVRObtIcaMN2Ldv\nH+x2O0pKSpCTkwMAaGhoQHl5OWJjY7F06VJMnjw56kGJiMiYUa/gi4uLceONN/o9V1VVhdWrV+Ou\nu+7C+++/H61s0qnex2N+uTgPXh6Vs4sUVg8+JeXLW4IlJCQIC0NEROKEVeCH3sY1Pj5eWBizUb2P\nx/xyGe3Bm20WjcrHX+XsIoVV4Id+4Gqz2YSFISIicUb9kDWQrq4uAANX8oM/B1NZWen723SwL6bK\n4xdffBEFBQWmycP85sp3zm6Hs77Td6U+2HMffPzhhx8iOzvb99jlaobT2Rx0fKDHmbYsACkhj58W\nnwVgsi8fAOTNmKH7uHBuPjJTEpQ7/nqPh/bgzZAnnMci2LSh/ZYADh48CLvdjqSkJOTm5qKoqAjn\nz59HRUUFNE3D8uXLkZ2dHfC95eXlKCwsFBZ2rA39y0lFzB9ddc1uHHW0Bn3d6XT6tWnys1Jwttlt\naB9G3xPOPhZPT8esrJQRz5v9+OtROTsAVFdXY9myZRFvZ9QCHwnVCzyRntEK/HCjFd+Fr/3PiN56\nsPcEGhvKPgIJVuBJHlEFPqwWDRFFLjPWBvT2+x4v2vYC6lb/yG9MTEc3MgN8yTDQWL3xSIiDyxO1\nazkyKRZ4Har/M4/55Rreohmhtx+OM//ye2r4446UeLjcfQHfPnys3vjpV10OxMaGkPpLKh9/lbOL\nxLVoiEziUPF9URlL4xd78ERhirQHn+nxBLwKHypT5wreyPi5X8tGe3/gP+q5E5OQlSru+yypaYlI\nm5AkbHvjEXvwRBQyT48HjjpXwNe86Ym4lCKuwM+ak80CbxJs0ehQfT0L5pfL6Fo0ZqNyftXPHVFY\n4ImILIoFXofqn8Izv1xG16K5Ye+rURkbLpXXs1f93BGFBZ6IyKJY4HWo3sdjfrmM9rAPf/s/ojI2\nXOzBq48FnojIoljgdajex2N+uVTuYQNq51f93BGFBZ6IyKJY4HWo3sdjfrmM9rDNNouGPXj1scAT\nmUTRvteiMpbGLxZ4Har38ZhfLpV72IDa+VU/d0RhgScyCa4mSaKxwOtQvY/H/HJxHrw8qp87orDA\nExFZFAu8DtX7eMwvl8o9bEDt/KqfO6KwwBMRWRQLvA7V+3jMLxfnwcuj+rkjCgs8EZFFscDrUL2P\nx/xyGe1hm20WDXvw6mOBJyKyKBZ4Har38ZhfLpV72IDa+VU/d0RhgScisigWeB2q9/GYXy7ek1Ue\n1c8dUVjgiUyCq0mSaCzwOlTv4zG/XCr3sAG186t+7ojCAk9kElxNkkRjgdeheh+P+eXiPHh5VD93\nRGGBJyKyKBZ4Har38Zg/dC53L+qa3YZ+ufs8uttUuYcNqJ1f9XNflDjZAYjMoMXdj6OOVkPvyc9K\niVIaIjF4Ba9D9T4e88vFefDyqH7uiMIreKIhMmNtQG9/SGNjOrqR6dFv0+iNT7EZjkdkSFgFfteu\nXfD8+0SdOXMmCgoKhIYyi8rKSqWvBJg/DL39cJz5V0hDO1Li4XL3BX29rbUVE9LTg47/6qxMv/Fm\nm0XjdDqVvYpX/dwXJawCn5SUhJtvvll0FiIiEiisAu/xeFBWVgZN05CXl4fCwkLRuUxB9SsA5pdr\n6NW7ilS9egfUP3dECavAFxcX+37eu3evsDBEpL7+Pg8azxubkRSu1LREpE1IGpN9qSjiD1nj4+N1\nXx/aCxucm6rK4xdffBEFBQWmycP80d2f0+lETPKXUx/bWgeK1OCV+PDHXW432lrdQV+/cP48UlNT\ng453d3X59ekX7Pg9Dn7zLr/tpSAFQPyI7d+w91W8u7hkRL5g4/UeIz3b9/sHvrxy//DDD5Gdne17\nPPz1YI8nTZmAxnpXyOMjeZybn4nrl8wD4P//c+g8eLOcz0Yfi2DTNE0z+iaHw4Hp06cDAHbv3o0V\nK1YEHFdeXq50+0b1D2qYP3R1zW4cdbQi0+MJ+UPWTIMfsg4f/9VZmfhHncv3+L8fXoYnni0PaR+B\nxuqNH76voaalJ+IrKSMv1ML9kHWgwLcZfl84Zs3JxqScka0w1c/96upqLFu2LOLthHUF73A4UFVV\nBQCYN29exCHMSuUTBGB+2diDl0f1c0eUsAp8UVGR6BxE4x5XkyTR+E1WHaqvZ8H8cvl63CEy4zx4\nVal+7ojCAk9EZFEs8DpU7+Mxv1zswcuj+rkjCgs8EZFFscDrUL2Px/xyGe3Bm201Sfbg1ccCT0Rk\nUSzwOlTv4zG/XEZ78GabRcMevPpY4ImILIoFXofqfTzml8toD95s2INXHws8EZFFscDrUL2Px/xy\nGe3Bm20WDXvw6mOBJzKJon2vRWUsjV8s8DpU7+Mxv1zswcuj+rkjCgs8kUlwNUkSLeI7OlmZ6n08\nK+TvaOtGZ0dP1PfV09mHTI8HKTZx2+Q8eHlUP/dFYYEnU+vs6EHd6aao7+eiuw+ft/bgq7Myo74v\ns/F4NVzUuTPVcMnxMUiOj41iIhKFBV6H6rf9skL+2TMKwnpvV58HXX3ekMf3eQ3fuXJUw2/ZZ1Y9\nHg2u9t4RzwfLPy090fQFXvVzXxQWeLKkrj4vPm8NvbWTGeCepESq44esOlS/AmB+uVSfB6/Cvz6C\nUf3cEYUFnojIoljgdag+l5b55VL9nqwqz+NX/dwRhQWeiMiiWOB1qN7HY365VO5hA2rnV/3cEYUF\nnojIoljgdajex2N+uVS/Jyt78OpjgScyCa4mSaKxwOtQvY/H/HKp3MMG1M6v+rkjCgs8kUlwNUkS\njQVeh+p9POaXi/Pg5VH93BGFBZ4i1tDQAE0Tv1hXNLW3upTLPBb6+/vQ2WG+wt5yyQWvN/TF42gA\nFxvToXofL1r5h67R/n8Vf8OP1/0nfrnxMay59wdC9zN7RgF6uvuFbhMATn5UgT+/sgkld23A/Ou/\nI3z7g1TsYf/5fzfh0zPH8V+/+rNp8tfZ/4H1G+/HqpL78B/3rAvpPar/2RWFBZ4MG7pGe/UHpwAA\nJ2pqUVcoft32SVMmCN+mq7keAHCx+bzwbavuYnM9+vp60N52ESmp4o99OBoaB/4/XWj4QnIS9bBF\no0P1Pp7q+QfvCerxaujzeEP+1R+Ftd3DwXnw8qh+7ovCK3gyvZaufjS0h762+2WJcbgs0dw3pFDZ\naHeASuv1+L3OO0DJwwKvQ/U+nur5B+8JqkGDx8BFudckH56qfk/WYPmD3QFqUJq73+9mKzLuAKX6\nuS8KC7xFjNXNqQFE5YNPIhIv7ALf0NCA8vJyxMbGYunSpZg8ebLIXKag0n0dA92c2ul0+q6CRYrG\nB5+BOJ1OTJoyd0z2FQ2q3JM1GBXy9/d50Hh+5GcFdvs5zJiRJ3x/qWmJSJuQJHy70RJ2ga+qqsLq\n1asBADt37kRJSYmwUEREoehy96Gx3jXieafTBW93qvD9zZqTrVSBD3sWTUpKiu/nhIQEIWHMRpWr\n92CicfU+llTPz3uyyqP6uSNK2FfwQ78FGB/PO9IHMp764p+eq0Wd/R/Ct3uxPQ2fOS/hX53BZ20M\ndyk+BmmJsWjsCPyetkvNAAD7Pz9CzQcHAABpiXHo6OlHkzMN5xs7QtrP4HtCNXz88H1N3/caarJy\nQtpHoLF64/V+X0Pf09Tg+Pd/P0N/f+APUjuS4tCqc74laOmo/+LLtkl/SzwuJoX/cV+dvRYA+E3W\nMIR91D0ej+9nm80mJEy0eT1eGJlfEWkPO1BfPFoC9cWj1YMf6njNMQDAP86cwPqN90d1X6J94ajF\nm44nZcfweRMA/hhaHiNjw7H9lcejtu1wHav6v5DHjsW5r4KwC3xXVxeAgSv5wZ+Hy8jIQHV1dbi7\nMIWWlosRvT/5K4KCjKKtq33Evq76ykQA7VHdV+6MyaiuGfh59uzZwvcVDQ0NDWhvHzguiYmJktOY\nS0/PwL84r7zyStMcm3/+858ABtrCyV/xP58DnfdA9M79+oZ21DcI3+wIGRkZYjakham+vl7705/+\npL3++utaY2NjuJshxT377LMaAO3BBx+UHSVkv/71rzUA2saNG2VHMZ05c+ZoALSTJ0/KjuLzxhtv\naAC0VatWyY6inLCv4HNycnD33XeL+VuGiIiE41o0REQWxQJPRGRRLPBERGPE6/X6zUCMNsM9eKNL\nFNjtdhw/fhyJiYkoLi5GQkICdu3a5ftNzpw5EwUFBeGlD4OR/Fu2bPHNDGlubsbatWsNb0M0EflV\nOf7nzp3D8ePH4fV6cf3112PatGmGtyGaiPyqHP+mpib87W9/Q3JyMgoKCjBjxgzD2xCpu7sbr7/+\nekTZZR77ffv2wW63o6SkBDk5I7/DMFSwY2z42Bv9VPbtt9/2/bxjxw7dsS0tLVplZeWI5/fv3290\nt8IYyd/V1eX7ee/evWFtQzQR+UUef6OzaIzk37Nnj+/nd999N6xtBBLJLBoR+c18/g+dRbNv3z7f\n84cPHw55G6INzqJZvHhxyPsNll3msdc0TTt16pRWX18/6rhgx9josTd8BW9kiYKamhpkZGSgrKwM\nc+fORX5+PoCBL0mVlZVB0zTk5eWhsLDQaIywGcmflDSw5kRTUxOys7PD2oZoIvKrcvw9Hg+8Xi80\nTfP75rQqxz9YflWOf2JiIrq6upCQkICamhosWLAACQkJ0o5/bOyXSw6Hm13msTci2DE2euwNF3jN\nwBIFTU0D3+K8/fbb8e6772LWrFmIiYlBcXGxb8zevXuNRoiIkfyDTp8+jUWLFkW0DVFE5Bd5/GNi\nYvz+Oxoj+efOnYstW7YAANasWRPWNgIxmnkoEfnNfP4PPTZLlizBvn370Nvbi/z8fHR2diIhIWHM\nz/9A/79G22+w7DKPvRHBjrHRY2+4wBtdoqCoqAgAMGXKFFy8eBFZWVl+r491gQxniYXu7m7f1XC4\n2xBFRP7si0jUAAAB5UlEQVShIj3+3/3ud7F//37cc889IY03kv/jjz/GL37xC3i9XuzevRu33nqr\n4W0EUlJSgsrKStxxxx2G3ysi/1BmO/8ffvhhVFRUID8/H7GxsVixYgWAgRVjU1NTQ9qGaDfeeCNW\nrlyJ6667LuT9xsXFBcw+lJnX0Ap2jI0ee8OXMMGWKDh58iROnz7tNzY/Px92ux0A0NLSgvR/r07n\ncDh8Y7q7u41GiIiR/MDAp97Dr/RCWaYhWkTkF3n88/Ly8M4772DBggUhjTeSf/CfoDExMX5/GCM9\n/ldddRV27dqFr3/964bfKyK/mc//tWvXYtu2bX7tkEuXLqGjo8P3+xnr8z8zMxM7d+7EnDlzdLMH\nMjy7zGMfTKD8wY6x0WNv+Ap+4cKFKC0thaZpWL58ue/5Y8eOwWaz+f4nAMC8efPw5ptv4pNPPsHl\nl1/uO8kdDgeqqqp8Y8aSkfzAwKJFs2bNCmkbY0FEflWOf15eHsrKygAAX/va10bdxlgQkV+V4//J\nJ5+gtrYWmqb5/WtH1vEXkV3msT948CDsdjuSkpKQm5vr624Eyh/s92r02Nu0oU0dIiKyDH7RiYjI\noljgiYgsigWeiMiiWOCJiCyKBZ6IyKJY4ImILIoFnojIoljgiYgs6v8BpAK8fSJFJiYAAAAASUVO\nRK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10f2faa50>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEBCAYAAABseY4cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHXhJREFUeJzt3XtsFOe5BvDH+IKzMdgExySUgKHGcQMo1EnAQZY5kUlC\nbxKKgJJGHKlq1FJoFOm0UVAq5dLTplWiSIekoalUKUdEDVWKwA1JoWATm6wJUOyaGCMHyDpsYhu7\nuwbHxODL7pw/kPd4zazxt/vN5R0/v3/w7I7nezyM3x2/++1MmmEYBoiIyDOmOB2AiIj0YmEnIvIY\nFnYiIo9hYSci8hgWdiIij2FhJyLymAzVb+ju7sahQ4dw0003YcmSJViwYIEVuYiIKEnKhf1f//oX\nNmzYAADw+/0s7ERELqPcipk6dSquXLmCSCSCpqYmDA4OWpGLiIiSlKb6ydPh4WHs378fg4ODyMnJ\nwX333YcZM2ZYlY+IiBQpF/bRqqqq8O1vfxtZWVk6MxERUQqUe+wjLl26hMuXL5sW9QMHDiA9PT2l\nYEREk01eXh7uueeelLejXNibm5vR2toKwzCwfv1603XS09NRWlqacjgiosmksbFRy3aU3zxdsmQJ\n1q1bh/Xr13u2BeP3+52OkBLmd5bk/JKzA/Lz68IPKBEReUxKb54mUlNTw1YMxcn+3e9wdetWp2MQ\nuVpjYyMqKytT3g7P2ImIPCbpWTFe5vf7UV5e7nSMpLkxv8rZuhvzq7A6v2EY6O7uRiQS0b7t3t5e\n5Obmat+uXSTkT09PR0FBAdLS0iwbg4WdSJju7m5MmzYNPp9P+7Znz56tfZt2kpC/v78f3d3dmDVr\nlmVjsBVjQvLZIsD8TrM6fyQSsaSokz18Pp8lf22NxjN2Qrh/EBf7hx0bf4YvAzN93pw6S+QEFnYT\nk63He7F/GEfO91qYCFi+43Uc+88tps+tmJcbV9gn2/5PldUvzHzhlYeFnWxR9tb2hIWdUmP1C/PY\nF166seHhYaSlpTl2aRUWdhOSzxYB5nea9Py6tLS04LnnnsOuXbtij23ZsgVbt27FHXfccd36//jH\nP/Dxxx/D5/PB5/Phhz/8Yey5Tz/9FG+//Tby8vKQkZGBDRs2pHxV2YaGBrz33nvw+XyYO3cuvv/9\n70/o+wzDwEsvvYSMjAxcvnwZW7ZsQX5+PgDg3Llz2LNnD5qbm/Gb3/zG9Oe0Aws72eLoxs1ORyCb\nLVq0CPfff3/cY1u3bsVtt91mun5TUxOefvpp0+eqqqrwzDPPaD0D3rNnD379618DAF599VV0dnbi\n9ttvv+H3vffee3jggQewbNky9PX14dVXX8Uvf/lLAEBRURGeeuop7Ny5U1vOZHBWjAnp15twY36V\nNowb86uQnj8V9fX1+NWvfoVXXnkFL774Iv79738DAK5evYrXXnsNzz77LC5cuBD3PSPPHT16FNu2\nbcO2bdvw2WefxZ7fuXMn/H4/fv/732Pbtm3X7d8vvvgCRUVF2Lt374Rzfvnll3Fn0ytXrkRDQ8OE\nvre1tRXLli0DAEybNm3CY9qJZ+xEpEVfXx9qa2vx7LPPArjWZ3744YcBANnZ2XjiiSdMz2RHnuvv\n78eTTz553fOPPvoogsGg6XMAUFBQgJ///Oe499574x5/4403cPTo0bjHNm3ahLKyMvT29mL69Omx\nx2+99VacPn16Qj9nNBq9Lr/bsLCbkN4jZX5nSc+frDNnzuC+++6LLWdkZKCiosLycbOysvDTn/70\nusc3bdqETZs2mX7P9OnTEQ6HY8u9vb0T7tkPD8fPQLLyE6TJYiuGiLQoLi7G8ePHY8tDQ0Oora21\nfNyhoSFs374dHR0dE/6e3NxcfPHFF7HlQ4cO4Zvf/GbcOu3t7di+fTuGhobiHl+8eDE++ugjANf+\nShn7vBukP//888/r3mhbW9uE3oRwK7/fj7lz5zodI2mq+Xv6h/B574CFia7NY2+/e5npc3fkZeMW\nX2ZsebLtf1V9fX1xvV2r///G/v8kMnXqVKSlpeHNN9/EyZMnUV1djenTp6O9vR1fffUVqqqqUF9f\nj7a2Npw8eRKzZs3CjBkzEAqF8Kc//Ql+vx9dXV04fvw47rzzzliL449//CNqa2vR1dWFY8eOIScn\nBwUFBbFxOzo68JOf/AR333037rzzzgn/XLm5uXjzzTdx9OhRzJkzJ9Y3H3Ho0CG88MILWL9+fVzb\npqSkBDt37sTRo0fxwQcf4Gc/+1ksa1tbG3bs2BH7OZubm3H33XcjIyO+OTL2/3BEZ2cnFixYMOGf\nIRFettfEZPuAzLlQv+MfUCrK//+PyE+2/a+qo6Mj7poo/ICSPGP/D0foumyvco+9ra0NDQ0NiEaj\nuP/++x2bp2klyUUFcGd+lVkxbsyvwu78M31ZLLwUR7mwt7a2Yu3atQCAv//9754s7EREkim/eRqJ\nRBCNRhGJRGBBF8cVpM9DZn5nSc9P8imfsS9evBgvv/wyAGDjxo3aAxERUWqUz9hPnjyJp59+Gk89\n9RROnDiRcL3RZy1+v1/U8mTL3xYIIBgMxpaDwaD25ZLXXkz4fFsgMKn3v+pyT0+PZ/9angwMw0BP\nT09s2ez4SZXyrJh9+/bhW9/61nVfjyZ9VsxkY8esmCcfXIRtB1tMnxs7K4bGd/nyZQwMDGDmzJlO\nR6EkhMNhTJ06FTk5Odc959ismPnz52P37t0Arl3kx4s43c5ZzD++nJwcDAwMKH0gZ6Ik3DN0PBLy\nJyrqOikX9pKSEpSUlFiRhTyMV3fUy6qz9UAggG984xuWbNsO0vPrwg8okS2tmPGwFUN0ja5WDK8V\nQ0TkMSzsJqTPQ2Z+Z0nOLzk7ID+/LizsREQew8JuQvKMDMCd+ZfveH3C67oxvwrJ+SVnB+Tn14WF\nnYjIY1jYTUjv07kxP+95KoPk7ID8/LqwsBMReQznsRPnsRO5BOexExGRKRZ2E9L7dG7MrzIrxo35\nVUjOLzk7ID+/LizsZIuyt7Y7HYFo0mBhNyF9LizzO0tyfsnZAfn5dWFhJ1vw6o5E9mFhNyG9T+fG\n/JzHLoPk7ID8/LqwsBMReQznsRPnsRO5hGO3xnv55ZexcOFCAEAoFMLjjz+ecggiItJHuRXzxBNP\nYM2aNVizZg3mzJljRSbHSe/TuTE/57HLIDk7ID+/Lspn7NnZ2QCA7u5uFBQUaA80GYX7B3Gxf1jb\n9q5k5eFcqH/C6/cPRbSNTUTOS7rHXltbi7KyslihH409djVO97iL8304o/BCoBt77ETXOH6tmKtX\nr5oW9RGj/yTy+/1cHme5LRBAMBiMLQeDQVuXw+GQo+O3BQKu+v/gMpedXNYhqTP2aDSK6upqPPTQ\nQ6bPSz9j9/v9tn6CTfcZezAYxNy5cye8vtvO2O3e/7pJzi85OyA/v6Nn7MFgEEVFRSkPTkRE+iVV\n2AsLC7FgwQLdWVxD8is+AKWzdbvwnqcySM4OyM+vCz95Srbg1R2J7MPCbkL3Gxl2G/3GpETS97/k\n/JKzA/Lz68LCTrbg1R2J7MNrxbgA57FzHjsR4IJ57ERE5E7KlxSYDKTPhVWdx+60geFo3CUQ2gIB\nzLdx1tUMXwZm+rK0bU/y8SM5OyA/vy4s7OS4voEIGtr7YsvB9q/QmW5fa2rFvFythZ3IaWzFmJD+\niu/Gs3WVeexuzK9C8vEjOTsgP78uLOxERB7Dwm5C+lxYN85jV7nnqRvzq5B8/EjODsjPrwsLOxGR\nx7Cwm5Dep5Peo5aeX/LxIzk7ID+/LizsREQew8JuQnqfzo09apVZMW7Mr0Ly8SM5OyA/vy4s7GQL\nXt2RyD4s7Cak9+mk96il55d8/EjODsjPrwsLO9mCV3cksk9ShT0QCOCvf/0r3n33XQwODurO5Djp\nfTo39qg5j10GydkB+fl1Ub5WzKVLl9DZ2Yl169ZZkYeIiFKkfMbe1NSEm2++Gbt378aZM2esyOQ4\n6X066T1q6fklHz+SswPy8+uiXNi7u7tx6dIlPPLIIzh79iyi0agVuYiIKElJ9dgrKioAAF/72tfQ\n09Njus7oXpff7xe1/Ic//MHW8doCgbi+cjAYTGn5xIkTSuuHwyGt45stl7z24oTHV82f6nJbICD6\n+NG5PPK1W/JMtvy6KN8ar6mpCTk5OSgqKsIHH3yA8vJyZGZmxq0j/dZ4fr+9F+vXfWs81Rtt2HFr\nvOU7Xk/4BurY8e2+UYjuW/PZffzoJDk7ID+/Y7fGW7p0KZqamrBnzx5kZmZeV9S9QPKBAbizR60y\nK8aN+VVIPn4kZwfk59clqTsorV27VncOIiLShB9QMiF9Lqz0eeDS80s+fiRnB+Tn14WFnYjIY1jY\nTUjv07mxR817nsogOTsgP78uLOxkC17dkcg+LOwmpPfppPeopeeXfPxIzg7Iz68LCzvZgld3JLIP\nC7sJ6X06N/aoOY9dBsnZAfn5dWFhJyLyGBZ2E9L7dNJ71NLzSz5+JGcH5OfXhYWdiMhjWNhNSO/T\nubFHzXnsMkjODsjPrwsLOxGRx7Cwm5Dep3Njj5r3PJVBcnZAfn5dWNiJiDyGhd2E9D6d9B619PyS\njx/J2QH5+XVhYSci8hilwr53715UVVWhqqoKzc3NVmVynPQ+nRt71CqzYtyYX4Xk40dydkB+fl2U\n7qCUnZ2NBx980Kos5GFlb21XegOViJKnVNgjkQh2794NwzAwf/580TesHo/0Pp30HrX0/JKPH8nZ\nAfn5dVEq7KtXr459vW/fPu1hyLt4dUci+yT95mlmZqbOHK4ivU/nxh4157HLIDk7ID+/LkqF/fz5\n87Gvr169Ou66o3ew3+8Xtdzc3GzreG2BQFwxCwaDKS13d3crrR8Oh7SOr7o8dnzV/KkutwUCoo8f\nLntrWYc0wzCMia58+PBhdHV1AQCWLl2KhQsXmq5XU1Pj2f67Fc6F+nHkfK9j4xfn+3Am1D9px18x\nLxdF+T7Hxica0djYiMrKypS3o9Rjr6ioSHlAIiKyFj+gZEJ6n86NPWrOY5dBcnZAfn5dWNiJiDyG\nhd2E9LmwbpwHznueyiA5OyA/vy4s7EREHsPCbkJ6n056j1p6fsnHj+TsgPz8urCwExF5DAu7Cel9\nOjf2qHnPUxkkZwfk59eFhZ1sUfbWdqcjEE0aLOwmpPfppPeopeeXfPxIzg7Iz68LCzvZgld3JLIP\nC7sJ6X06N/aoOY9dBsnZAfn5dWFhJyLyGBZ2E9L7dNJ71NLzSz5+JGcH5OfXhYWdiMhjWNhNSO/T\nubFHzXnsMkjODsjPrwsLOxGRx7Cwm5Dep3Njj5r3PJVBcnZAfn5dkirsV65cwW9/+1t0dHTozkNE\nRClKqrDX1dXh4Ycf1p3FNaT36aT3qKXnl3z8SM4OyM+vi3Jhv3jxInw+H7Kzs63IQ0REKVIu7H6/\n3/OvitL7dG7sUfOepzJIzg7Iz6+LcmHv6urCu+++C7/fj1OnTlmRiTyIV3ckso9yYX/88cexZs0a\nlJeXY/HixQnXG/3K6ff7RS3bnb8tEIg7Sw0Ggyktjzw20fXD4ZDW8c2Wx2Ybb3zV/KkutwUCoo8f\nncvl5eWuyjPZ8uuSZhiGofpNPT092LVrFwoLC/HQQw9d93xNTQ1KS0u1BJwMzoX6ceR8r2PjF+f7\ncCbUb+kYy3e8nnDKox3jj2fFvFwU5fscG59oRGNjIyorK1PeTlKzYm655Rb8+Mc/Ni3qXiC9T+fG\nHjXnscsgOTsgP78u/IASEZHHsLCbkD7rR/o8cOn5JR8/krMD8vPrwsJOROQxLOwmpPfp3Nij5jx2\nGSRnB+Tn14WFnYjIY1jYTUjv07mxR817nsogOTsgP78uLOxERB7Dwm5Cep9Oeo9aen7Jx4/k7ID8\n/LqwsBMReQwLuwnpfTo39qh5z1MZJGcH5OfXhYWdbMGrOxLZh4XdhPQ+nfQetfT8ko8fydkB+fl1\nYWEnWxzduNnpCESTBgu7Cel9Ojf2qDmPXQbJ2QH5+XVhYSci8hgWdhPS+3TSe9TS80s+fiRnB+Tn\n1yXD6QBEThsYjuKcxjs4XcnKU9reDF8GZvqytI1PpFzYjxw5gq6uLgDXeqH33HOP9lBOk96nc2OP\nerxb441ld/6+gQga2vv0bTB9JjoVbnW4Yl6uawq79GNfen5dlAv7ihUrYl8fOHBAaxgiIkpdUj32\nzz77DK+88gqKiop053EF6X06N/aoJ9M9TyXnl37sS8+vS1KFvbCwEFu2bME///lP3XmIiChFSc+K\nyc7OxrRp0xI+P/qV0+/3i1q2O39bIBB3lhcMBlNaHnlsouuHwyGt46sujx1fNb/u8e3e/22BgGuO\n//Lycsd//yZzfl3SDMMwVL6hp6cHt9xyCwDg/fffx3e+853r1qmpqUFpaamehJPAuVA/jii82aZb\ncb4PZzTOCuH4albMy0VRvs+x8ck9GhsbUVlZmfJ2lM/Y6+vrsWvXLuzatQuLFi1KOYAbSe/TubHH\nO5nueSo5v/RjX3p+XZRnxXzve9+zIgd5XNlb25XeQCWi5PGTpyakz4V14zx2FczvHOnHvvT8urCw\nky14dUci+7Cwm5Dep3Njj5fz2GWQfuxLz68LCzsRkcewsJuQ3qeT3OMFmN9J0o996fl1YWEnIvIY\nFnYT0vt0buzxch67DNKPfen5dWFhJyLyGBZ2E9L7dG7s8U6me55Kzi/92JeeXxcWdiIij2FhNyG9\nTye5xwswv5OkH/vS8+vCwk5E5DEs7Cak9+nc2ONVmRXjxvwqJOeXfuxLz68LCzvZouyt7U5HIJo0\nWNhNSO/TSe7xAszvJOnHvvT8urCwky14dUci+7Cwm5Dep3Njj5fz2GWQfuxLz6+L8h2UWlpa8Mkn\nnyAajaKsrAxz5syxIhcRESVJ+Yw9HA7jkUcewdq1a/Hxxx9bkclx0vt0knu8APM7SfqxLz2/LsqF\nvaKiwoocRESkSdI99oMHD2LZsmU6s7iG9D6dG3u8nMcug/RjX3p+XZIq7HV1dVi4cCHy8/MTrjP6\nTyK/38/lcZbbAoG4P9+DwaCty+FwyPLxent7HR3f6Z9/vOW2QMBVxyOXnV3WIc0wDEPlGz788EPc\neuutKCkpSbhOTU0NSktLUw7nFL/fb+sr/7lQP46c773xihMUDAaVzhqL8304E+rXNr6qseOr5tc9\nfqpU86+Yl4uifJ+28VNh97Gvm/T8jY2NqKysTHk7SrNiuru7UV9fj5KSErS2tqKvrw8bN25MOYTT\nwv2DuNg/HFu+kpWHczYWuv6hiG1jEZH3KRX2goICbN261aosjrnYPxx/xpw+E50az6BvpFjz2Zrk\nHi/A/E6SfLYLyM+vCz+gRETkMSzsJiTPQwbcmZ/3PJVB+jxw6fl1YWEnW/DqjkT2YWE3IblHCjC/\n0yTnl96jlp5fFxZ2sgWv7khkHxZ2E5J7pIA786tc3dGN+VVIzi+9Ry09vy4s7EREHsPCbkJyjxRg\nfqdJzi+9Ry09vy4s7EREHsPCbkJyjxRwZ37OY5dBeo9aen5dWNiJiDyGhd2E5B4p4M78vOepDNJ7\n1NLz68LCTkTkMSzsJiT3SAHmd5rk/NJ71NLz66J02V4i0m9gOGrr9f/HmuHLwExflmPjk34s7CYk\n90gBd+ZfvuP1CffZ3ZhfhWr+voEIGtr7LEpzYyvm5cYKu/QetfT8uiTViolGo4hEeNcfmjhe3ZHI\nPsqFff/+/XjjjTfQ1dVlRR5XkNwjBZjfaZLzS+9RS8+vi3IrZvXq1Th9+rQVWcjDeHVHIvtwVoyJ\nydbjtQPnscsgvUctPb8uLOxERB5jWWEf3evy+/2uXm4LBOL6oidOnIhbDgaDli6HwyGt21PNr3v8\nVH9+7n97939bIBD7fRj53XDT76fKsvT8uqQZhmGoftPp06eRl5eH2bNnmz5fU1OD0tLSlMPZ5Vyo\nH0fO98aWg8GgrX9OF+f7cEbjPGbV/LrHVzV2fO5/e62Yl4uifB+Aa0VGcjtDev7GxkZUVlamvB3l\nM/bq6mr4/X5UV1fj8OHDKQdwI8k9UsCd+VWu7ujG/Cok55dcFAH5+XVRnhWzatUqK3IQEZEmfPPU\nhOR5yIA78/OepzJInwcuPb8uLOxERB7Dwm5Cco8UYH6nSc4vvUctPb8uLOxERB7Dwm5Cco8UcGd+\n3vNUBuk9aun5dWFhJ1vw6o5E9uH12E1I7pECzO80aflH3+jjtpJS22/6ofNGH+yxX8PCTrbg1R3d\ny003+iA92IoxIblHCrgzP+exyyA5O8Ae+wgWdgXDQ4O43HvR6RhEKYtEhtF3Mex0DLIIC7uJRD3S\n//3v/8ILj63CV19esjmRGmk93rGY33rv/M8LeG7DAwhfaI97XEL28bDHfg0Lu4JQRxBDA1fxZU/I\n6ShEKQm1BxGNRnCxq8PpKGQBx988HRyOItBzBVH1qwdrkTFlCqakxT9m92VjdXNj/uU7Xp9wn92N\n+VVIzi85OyD/sr26OF7YI4aBkx19GIg4U9jzbsrAXQU3OzI2EZEV2IoxIfmMBXBnft7zVAbJ2QH2\n2EewsBMReYxyYb9w4QL+/Oc/4y9/+Qs6OzutyOQ46XN5md9ZkvNLzg5wHvsI5cJ+/PhxPPbYY9iw\nYQOOHTtmRSYiIkqB8punPp8v9nVWljc/Biy9z+jG/CqzYtyYX4Xk/E5kH32tmlQlc60bndeqcQvl\nwm6MmpaYmZmpNYzbDQ8NAQA+P9uC4cEBbdtN687G55euatueG8cvf2s7Pi/7D8fGH89kHP/zsy22\njjceXqtGP+XCHolEYl+npaWNs+bETElLw7wZNyHi0Dz27Izru1GJ5vKGOq71H99+6RnLc3nNKwCw\nZb3TMWiMi93x75NJn8cuPb8uyoX9ypUrAK6duY98PVZeXh4aGxsnvE2nXyu/BLB41GvU4nmZABK/\nMZyZmYnCwkKrYxFZ5uzZswCARXkGssKfxo7/Gx37lgjH//6lIpn8XwY70eiS94zz8vK0bCfNMNRO\nlTs6OlBXVwfDMLBq1SoUFBRoCSLBokWLcPr0aZw6dQqLFi1yOg5R0lauXInDhw+jtrYWK1eudDoO\naaZ8xj579mw8+uijVmQhIiIN+AElIiKPYWEnIvIY2wp7NBqNm1EjTTQadToCkW0k/75Kzq5LUld3\nvHDhAmpqapCeno6VK1fi9ttvH3f9/fv3IxAIYM2aNZg9e3bs8b1798b+A77+9a9jyZIlycRRlmz+\n4eHhlLaji+q4ida3c/+rZE60rlP7W3VsN+zvG2UaT6LfVwnHu9tqDaCWv6WlBZ988gmi0SjKysow\nZ84c5W0AAIwk/O1vf4t9vWfPngl9T0tLi9He3h732IEDB5IZPmXJ5i8uLjYAGKdOnUp6Ozqojpto\nfTv3v0rmROs6tb9Vx3bD/h5rbKaKigoDgFFbW2u6vtnvq5Tj3U21xjDU8tfV1cW+fv/995PahmEY\nRlJn7LouKxCJRLB7924YhoH58+ejtLQ06W2p0JXfqcsrqI6baH07979K5kTrOnk5Cx35nTrex8tk\n9zacGlfKvq+oqEh5G0CSrRhD02UFVq9eHft63759SW9HVbL5p0yZEvevrv2gSnXcROvbuf9VMida\n16n9rTq2G/b3jTKNPZaT2YZddIzrpn0/EQcPHsSyZcuS3sa4hT0UCuGdd96Je+y73/2u9ssKANYc\nKLrz/+hHP8LJkydRXFwMQP/lFcbSlX8i61v9i6qSOdG6Vu/vZDIlu67dL0xjM23evBn5+fm49957\nk96GXXSP6/S+v5G6ujosXLgQ+fn5SW9j3MKen5+PzZs3X/d4Q0MDAPPLCpw6dQpTpkzBXXfddcPB\nz58/j3nz5gEArl7VfxEk3fl/8IMf4Be/+EVseSKXV0iFrvyJclq9/0dLlEElr9X7ezw68tu5v8ca\nm2ndunVYt26d0u+rU/tfZd8n4qZ9P8Is/4cffohZs2Zdd8kS1X2fVCtm+fLl2LlzZ+yyAqN99NFH\nSEtLiwtbXV2NQCCA7OxsFBYWxvpI58+fx/HjxwEAS5cuTSZKUnTlH287bsqfaH0793+iDCp5ndrf\nuvI7dbyPl0ny8S6h1qjk7+7uRn19PUpKStDa2oq+vj5s3Lhx3G0konytGCIicjd+8pSIyGNY2ImI\nPIaFnYjIY1jYiYg8hoWdiMhjWNiJiDyGhZ2IyGNY2ImIPOb/AMM9NipeLEx2AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10f7563d0>"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "phit_hit:\n",
        " \n",
        "\tMean             SD               MC Error        95% HPD interval\n",
        "\t------------------------------------------------------------------\n",
        "\t0.863            0.019            0.001            [ 0.823  0.901]\n",
        "\t\n",
        "\t\n",
        "\tPosterior quantiles:\n",
        "\t\n",
        "\t2.5             25              50              75             97.5\n",
        "\t |---------------|===============|===============|---------------|\n",
        "\t0.823            0.852           0.863          0.876         0.901\n",
        "\t\n",
        "\n",
        "phit_miss:\n",
        " \n",
        "\tMean             SD               MC Error        95% HPD interval\n",
        "\t------------------------------------------------------------------\n",
        "\t0.875            0.042            0.002            [ 0.782  0.942]\n",
        "\t\n",
        "\t\n",
        "\tPosterior quantiles:\n",
        "\t\n",
        "\t2.5             25              50              75             97.5\n",
        "\t |---------------|===============|===============|---------------|\n",
        "\t0.776            0.85            0.879          0.905         0.94\n",
        "\t\n",
        "\n",
        "diff:\n",
        " \n",
        "\tMean             SD               MC Error        95% HPD interval\n",
        "\t------------------------------------------------------------------\n",
        "\t-0.011           0.046            0.002            [-0.102  0.083]\n",
        "\t\n",
        "\t\n",
        "\tPosterior quantiles:\n",
        "\t\n",
        "\t2.5             25              50              75             97.5\n",
        "\t |---------------|===============|===============|---------------|\n",
        "\t-0.095           -0.043          -0.016         0.016         0.097\n",
        "\t\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There is a very small difference between this and the Bernoulli likelihood results, which we attribute to not recalibrating K relative to our repackaging of all observations into two Binomial observations."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 25
    }
   ],
   "metadata": {}
  }
 ]
}