negative_binomial_dist.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "import pymc as mc"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Parameterizing the Negative Binomial Distribution\n",
      "=======================================================\n",
      "\n",
      "The PyMC docstring for pymc.negative_binomial_like says:\n",
      "\n",
      "Negative binomial log-likelihood. The negative binomial\n",
      "distribution describes a Poisson random variable whose rate\n",
      "parameter is gamma distributed. PyMC's chosen parameterization is\n",
      "based on this mixture interpretation.\n",
      "\n",
      "$$\n",
      "    f(x \\mid \\mu, \\alpha) = \\frac{\\Gamma(x+\\alpha)}{x! \\Gamma(\\alpha)} (\\alpha/(\\mu+\\alpha))^\\alpha (\\mu/(\\mu+\\alpha))^x\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I like calling the first parameter $\\mu$ for mean, but I like calling the second parameter $\\delta$ for dispersion.  And let's call the value $k$ instead of $x$, since it is a count:\n",
      "\n",
      "$$\n",
      "    f(k \\mid \\mu, \\delta) = \\frac{\\Gamma(k+\\delta)}{k! \\Gamma(\\delta)} (\\delta/(\\mu+\\delta))^\\delta (\\mu/(\\mu+\\delta))^k\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# plot the log-likelihood of this distribution for a few values of delta\n",
      "\n",
      "x = arange(0,30,.5)\n",
      "mu = 10.\n",
      "for delta in [5., 10., 100.]:\n",
      "    y = [mc.negative_binomial_like(x_i, mu, delta) for x_i in x]\n",
      "    plot(x, y, label='$\\\\delta=%.2f$'%delta, drawstyle='steps-post')\n",
      "legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "<matplotlib.legend.Legend at 0x3bf6b90>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD9CAYAAACyYrxEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGOFJREFUeJzt3XuMVOXBx/HfYbktBeosNxt2A1pB2LvcaUEHcIuwYKFq\nIIoJbf8gEZVC2tKXpmG3KddKKUqbaEuxgX8KlVUbrqtkiyndghBSEWwj1MqQLq17QWChoDzvH7zM\ny7IX5pw5Zy7PfD/JJDtnzpzznDkzv3n2Oc88j2OMMQIAWKtTsgsAAAgWQQ8AliPoAcByBD0AWI6g\nBwDLEfQAYDnfg37JkiXKz89Xfn6+ZsyYofr6er93AQBwwfegnzlzpo4fP64TJ06osLBQP/nJT/ze\nBQDABd+DftKkSerU6cZmv/rVr+rs2bN+7wIA4ELnIDf+yiuvaO7cuW0+5jhOkLsGACt5GczAU42+\nrKxMRUVFrW5/+MMfouusWLFCXbt21VNPPdVhgW28LV++POll4Pg4Po7PvptXnmr01dXVHT7+29/+\nVjt37tT+/fs9FQpJkJMjNTbGvn5lZWzrhUJSQ4O3MgHwhe9t9Hv27NHatWv15ptvqnv37n5vHm7k\n5EiOE9tNkoyJ7bZ8eezruvnyABAI39von3vuOV29elVlZWWSpPHjx+uXv/yl37tJaeFwONlFuKGx\n8UbY+szV8YVC//9FEsu6KVD7T5nzFxCOL/M4Jp6Gn3h27DhxtTkhBo4TSNAHJt3KCySY19wMtNcN\ngPSUk5OjRprdkiYUCqnBx/9uqdGnGxcXTRuzHeUs9f81DnUPqWFpAE0sbi8Ip0hTj434fCZXe6+/\n1/NC0Kcbx5FTEduqQQVyzpocNV6JLZAD+1KQaOoJEJ/P5CLoM53jyFEAtfSAKseBfikQ9IHh85lc\nfgc9bfQpoNP/5Mh0jy0MjYLJtps9MWPh5kvBTXA7lfxaGggCQZ8CTPdGmeUxpndFMGHopjYf1JeC\na2nYdRNIBoI+AG6vKaoiqJIEw01euhnSKNQ95KpWH1oaUkOsF5sZWwkZjKAPgNvfKTX0UOw19VDI\nU5mSxV2lu0HGzZcITT2I03333adIJKK77rpLP/3pT/X000+3ud727dv14Ycfqnv37urbt290vfaW\npxouxgbA+UGOlB17ld5UiIuKcn9t1al0Ym/y4sKtK+ny+dy/f79OnjypqqoqrV69WqNGjXL1/Jdf\nfllTp05Vbm6uOnduu957/vx5TZ48WUeOHJF049f+b7zxhrp169bm8v79+8d3UOJibHrIdtHmLgXW\n7p5u3NT+Jcn5QexNPakfWfDi4MGDWrhwoRYuXOjp+d27d9fgwYM7XKempkb5+fnR+6Wlpdq3b596\n9erV5vJ58+Z5KkuQCPoYuWp3rwiyJPZye63UcRpir6TzZWqlRx55RHPmzNGuXbuiNfIPP/xQmzZt\navc5Y8eO1axZsyRJtbW1unr1qi5cuKD77rtPjz76aKv1T506pb59+0bv9+nTR6dOnVLv3r3bXJ6K\nCPoYuWl3d2IcwReJ05At5dBDx1d+Xd/22kL00Ucf6cCBA8rKytLatWu1bNkySTfa3VetWhXTNqZO\nnRoN/VGjRunBBx/UXXfd1WKd5ubmFs06nTt31uXLl5WVldXm8lTk+zDFQKLcbOqJ5dZnqRha2Wex\nvpx3unnx6aefatGiRVq0aJGWLVum3bt3e9pOeXl59O9evXqppqam1To9evRocf+///2vevTo0e7y\nVESNHmnLVTdP2vOtsm3bNs2YMUNZWVnKyclRv379oo/dqelmzJgxmj17tl599VXt2LFDb775piTp\n0qVLbV6QHTJkiE6ePBm9f+7cOU2ZMkW9e/dusbyurk4PP/ywH4fnO3rdxMhNpw1XvUHcbhyeuHqJ\nOR8p//n8+c9/ruHDh2vq1KnauHGjBgwYoCeeeMLVNt555x199tlnmjRpki5fvqz8/Hy9//776tGj\nhz744AMNHTpUnTp10qeffqrJkyfr3XfflSSVlJSourpa3bp105QpU1otT8VeNxkd9G4usDo/iH2Y\nAsZsST1uzrUR5yMVPp8daWpq0oYNG3T//ffr1KlT+uEPf+hpO5s3b1Zzc7MikYhmzZqlsWPHSpKG\nDh2qrVu3asyYMZKkLVu26J///Kccx9HAgQM1f/78DpfHi6D3tQwB1dIZbje98cWbEp/PTEY/+nQQ\n0BR+SJxY2/MDHYYZ8AlBD9ym0QnJVMT2H1lDdqO0NOACAXEi6IHbhK7HXkOPuW8+kET0owcAy1Gj\nB+IUa3u+cyWk66toz0fiEfRAnGLtjcWwykgW65pubs5+FMstzYZ2BwBPrKvRu+nZmLMmR05l7D+C\nAoB0ZF3Qu9F4xeW48UCcGEATyZDRQQ/EzcVsKQ3ZjnJiHBe/8XJIEkkPfxD0sXIzrAGN/5nDRbU7\n5MQ+jAYXbhPn2LFj2rJli9atWxddFutcsGkzl6wJyAsvvGAcxzH19fVtPh7Urt1sVhVuVg7spUKm\ncPEecvXeDECA0eCrt99+22zcuNFMmTLFHD582PXz161bZ2bPnm3mz58fXdbU1GRGjBgRvT9u3Dhz\n7ty5Vs9tb71Yn9+R9l5/r+clkF43Z86cUXV1tQYNGhTE5oGMEGvvMce58Q9nJjp48KCefPJJvfXW\nW64nBpekJUuW6Otf/3qLZe3NEXu79taL9fmJFEjTzZIlS7R27dpWLyCA2LkZFy9TR2KId85YSa1G\ng2xvjtjbpdNcsr4H/RtvvKHc3FwVFxffcd2Kioro3+FwWOFw2O/iAAiIX9cRvPZ882POWOnG0L+3\nam+O2NslYi7ZmpqaNqc3dMtT0JeVlamurq7V8hUrVmjVqlUt/k25/dvyVrcGPYD/F+oe+9SHkqSl\nie+lk8yuyTfnjN2xY4dGjx6tZcuWRYPerdszqkePHjp//nz0fntzwba3XqzPj8XtFeDKykpP2/EU\n9NXV1W0uP378uP7xj3+opKREkhSJRDRy5EgdOnTIl+m1gEzhdox7p9LJqD76fswZe9PtNfrb54ht\nby7YtJpL1tMl3BgNHjyYXjfATQG+h9y8l2MpRsDRELf169ebPXv2GGOMeemll8y2bds8b2vz5s0t\net2cP3/ejBw5Mnq/uLg42mvm5MmT5vPPP+9wvaampnafH6v2Xn+v5yXQqQTvvfdevfvuu8ppo0tA\nUFOVBTY9INPLIV5uf4vhotrt5r0cy1s51acS9GvO2I0bN2rbtm06c+aM5s+fr8WLF6t3797tzgWb\nqLlkM3LO2KA+HwQ9UpbL91umBb3tMnLOWAYqAzrm5uKt84OQHCfNG+nhSloEvRsMVIZM5ObirVPp\nxFCjj7NASCnWjUcPAGiJoAcAyxH0AGA5gh4ALGfdxVhXGGMeGcj18ApIe2nRj54fQSHjJPm9ST/6\n5PK7Hz1NN0AqujlFIQPSwweZ3XQDpCq3o47R8R0dIOgBZDQ3c8am7VyynoZC84GbXTMiJXAHPr+X\nkxgNriRyzthEziXb3uvv9bzQRg8gbSVyzth0nkuWphsAaSuRc8bGOhdsKs4lS9AD8MavC8Aeu3Em\nes7YWOeCTcRcsm7RdAPAmxtXBuK/eXBzzthFixZp2bJl2r17dxyH0XrO2FvdOhdsW8tvF+/zg0CN\nHkDaSeScsW7ngk3JuWQ9XcL1gZtd0+sGuIMM63WTjDljO5oL1u+5ZNt7/b2el/QYAuEHOVJ27LNG\nxTwJA0MgwBY+z7eZ6kMgJGvO2ETNJZuRc8a6Gr/GXSEIemSeGN73qR70tmOsGwCAKwQ9AFguqUEf\n68B8AADvkhr0AXazBQD8H5puAMBy/GAKQCuhUKjVD4mQOCGfpy4l6AG00uB24hOXctbkqPFKAL+N\ncSlTeljbF/RM+A2kPDfBzUTm8Qukjf6ll15SSUmJioqK9L3vfS+IXbSvsTH2q7wB11oApDY3U/Om\n87S8vtfod+7cqb179+rIkSPq3Lmz6uvr/d4FAPjCTV0vnS9Z+F6j//Wvf62lS5dGx13u06eP37sA\nEI9MqcYiyvca/QcffKC9e/fqueeeU3Z2ttatW6evfOUrba5bUVER/TscDiscDvtdHAC3y5RqrAVq\nampUU1MT93Y8BX1ZWZnq6upaLV+xYoWuX7+uCxcu6NixYzp8+LAee+yx6Ghtt7s16AEALd1eAa6s\nrPS0HU9BX11d3e5jL774or7xjW9IkkaPHq2uXbvq3Llzuvvuuz0VEEBmC3UPuep5E1R3zJstXrGu\nm0p9PXxvuikvL9f+/fv10EMP6e9//7uam5vVv39/v3cDIEO4De2gumOmc4uX70H/7LPP6lvf+pYK\nCwslSa+++qo6dWKkBQBIFt+DvkuXLtqyZYvfmwUAeERVGwAsR9ADgM9S7acK9o11AwBJlmoXbgl6\nAFZx0x0zyJExUwlBD8AqjIzZGm30ANqXao3N8IQaPYD2pVpjs4U6+sWtX5OiEPQAkESJGCqBphsA\nsBxBDwCWo+kGQMbKlK6YSQ16Ny8wAPgtU7piJjXozXKfLikDANpFGz0AWI6gBwDLEfQAYDmCHgAs\nR9ADgOXoRw/AHx0N2tLWuon47b+P0rnPvWOMX8PmuNyx4yjmXefkSI2Nsa2bhm8gIOM4jn8jdqUg\np9IJpPu4q9y8RXrU6BsbrX5TAECQaKMHAMsR9ABgOYIeACxH0AOA5Qh6ALBcevS6AYA0kmp97tOj\nH73lfW6BjMNnOspNn3uv/eh9b7r505/+pNLSUhUWFqqkpEQHDx70excAABd8D/qlS5dqzZo1On78\nuFavXq2lS5f6vQsAgAu+B31eXp7Onz8vSWpqatKgQYP83gUAwAXfL8auXr1aEyZM0He/+11dv35d\nf/7zn9tdt6KiIvp3OBxWOBz2uzgAUpHlA6D5paamRjU1NXFvx9PF2LKyMtXV1bVavmLFCr344ota\nuHChZs+ere3bt+uVV15RdXV16x1zMRZALCz//CfiYqzvvW569uypixcvSpKMMerVq1f0fosdE/QA\nYmH55z9nTY4ar7Q9Ou/tXwApM3rloEGD9Mc//lEPPfSQ9u/fr3vuucfvXQCANRIxbr3vQf+rX/1K\nzzzzjK5du6Zu3bpp06ZNfu8CAOACP5gCkNr4/EelzA+mAACphaAHAMsR9ABgOYIeACxH0AOA5Qh6\nALAcQQ8AliPoAcByBD0AWI6gB5Dabg5pHOstJyfZJU45DIEAwC4W5wVDIAAA2kTQA4DlCHoAsBxB\nDwCWI+gBwHIEPQBYjqAHAMsR9ABgOYIeACxH0AOA5Qh6ALAcQQ8Aluuc1L07TmzrhULBlgOAPW6O\ndhnrug0NwZYnBaTH6JUAEIQ0G+mS0SsBAG0i6AHAcgQ9AFjOc9Bv375dBQUFysrK0tGjR1s8tmrV\nKuXn56uoqEj79u2Lu5AAAO8897opKipSVVWVFixY0GL5kSNHtGPHDr333nuqq6vThAkT9Le//U1d\nu3aNu7AAAPc81+iHDRumoUOHtlq+c+dOzZ07V1lZWRo4cKAKCgp06NChuAoJAPDO9370Z8+e1eTJ\nk6P3c3NzFYlE2ly3oqIi+nc4HFY4HPa7OACQtmpqalRTUxP3djoM+rKyMtXV1bVavnLlSs2cOTPu\nnd8a9ACAlm6vAFdWVnraTodBX11d7XqDubm5OnPmTPR+JBJRXl6e+5IBAHzhS/fKW3+pNX36dP3u\nd7/TZ599pkgkouPHj2vMmDF+7AYA4IHnoK+qqlJeXp5qa2tVXl6uadOmSZJGjhyp2bNnq7i4WI88\n8ohefvlldenSxbcCAwDcYawbAJmLsW4AADYg6AHAcgQ9gMx1c+z6WG45OckurWe00QNALFKgPZ82\negBAmwh6ALAcQQ8AliPoAcByBD0AWI6gBwDLEfQAYDmCHgAsR9ADgOUIegCwHEEPAJYj6AHAcgQ9\nAFiOoAcAyxH0ABCLNB67nvHoAcBvAY1dz3j0AIA2EfQAYDmCHgAsR9ADgOUIegCwHEEPAJYj6AHA\ncgQ9AFiOoAcAy3kO+u3bt6ugoEBZWVk6cuRIdPm+ffs0YsQIFRcXq6ioSHv37vWloAAAbzp7fWJR\nUZGqqqq0YMECOY4TXX733Xdr79696tevn95//31NmTJF//rXv1qsAwBIHM9BP2zYsDaXFxcXR/8u\nKCjQ9evXdeXKFWVnZ3vdFQAgDp6DPha///3vVVJS0m7IV1RURP8Oh8MKh8NBFgcA0kpNTY1qamri\n3k6Ho1eWlZWprq6u1fKVK1dq5syZkqRJkyZp3bp1GjFiRIt1Tpw4oUcffVTV1dW65557Wu+Y0SsB\n2CrFRq/ssEZfXV3tqTCRSESzZ8/Wli1b2gx5ALDazbHrY123oSHQ4vjSdHPrN0xTU5PKy8u1evVq\njR8/3o/NA0B6cRPcCeio4rl7ZVVVlfLy8lRbW6vy8nJNmzZNkrRx40adOnVKP/7xj/XAAw/ogQce\n0CeffOJbgQEA7jDDFAAkk4v2fGaYAgC0iaAHAMsR9ABgOYIeACxH0AOA5Qh6ALAcQQ8AliPoAcBy\nBD0AWI6gBwDLEfQAYDmCHgAsR9ADgOUCnUoQAHAHHU1S4tMIvwxTDABpgmGKAQBtIugBwHIEPQBY\njqAHAMsR9ABgOYIeACxH0AOA5Qh6ALAcQQ8AliPoAcByBD0AWI6gBwDLEfQBqKmpSXYRAsXxpTeO\nL/N4Dvrt27eroKBAWVlZOnr0aKvHP/74Y/Xs2VPr1q2Lq4DpyPY3GseX3ji+zOM56IuKilRVVaUH\nH3ywzceXLFmi8vJyzwUDAPjD88Qjw4YNa/ex119/Xffee6++8IUveN08AMAvJk7hcNgcOXIkev/C\nhQtm/Pjx5tKlS6aiosK88MILbT5PEjdu3Lhxc3nzosMafVlZmerq6lotX7lypWbOnNnmcyoqKrR4\n8WL16NGjw5lQOnoMAOCfDoO+urra9QYPHTqk1157Td///vfV1NSkTp06KTs7W88884znQgIAvPNl\ncvBba+cHDhyI/l1ZWalevXoR8gCQRJ573VRVVSkvL0+1tbUqLy/XtGnT/CwXAMAvnlr247R7925T\nWFhohg8fblavXp2MIgRq0KBBpqioyJSWlprRo0cnuzhx+eY3v2n69+9vCgsLo8vq6+vNww8/bIqK\niszXvvY109jYmMQSxqet41u+fLkZOHCgKS0tNaWlpWb37t1JLGF8Pv74YzNx4kRTWFhohg4datas\nWWOMsecctnd8NpzDy5cvm1GjRpnS0lIzZMgQ853vfMcYY8zp06fNuHHjTGFhoZkzZ465evXqHbeV\n8KC/cuWKGTx4sIlEIubatWtm1KhR5ujRo4kuRqAGDx5s6uvrk10MXxw4cMAcPXq0RRA+++yzZv36\n9cYYY9avX2+ef/75ZBUvbm0dX0VFhVm3bl0SS+Wfuro689577xljbvSIGzJkiDl27Jg157C947Pl\nHDY3NxtjjLl27ZoZO3as2b9/v5kxY4apqqoyxhizaNEi87Of/eyO20n4EAh/+ctfVFBQoIEDB6pz\n586aM2eOdu7cmehiBM5Y0qto4sSJCoVCLZbt2rVLTz/9tCRp3rx5aX3+2jo+yZ7zN2DAABUWFkqS\nevbsqeLiYp09e9aac9je8Ul2nMPs7GxJ0tWrV/X555+rf//+qq2t1axZsyTFfu4SHvSRSER5eXnR\n+7m5uYpEIokuRqAcx1FZWZmKi4u1cePGZBfHd//5z3/Up08fSVLfvn3173//O8kl8t8vfvELDR8+\nXPPmzVNDQ0Oyi+OLjz76SIcPH9aECROsPIc3j2/ixImS7DiH169fV2lpqQYMGKBJkyYpFAqpb9++\n0ccHDhwYU34mPOgdx0n0LhOutrZWR48e1dtvv63NmzfrrbfeSnaR4MLChQt16tQpnThxQl/+8pf1\n/PPPJ7tIcbt48aIef/xxbdiwQb179052cXx38eJFPfHEE9qwYYN69eplzTns1KmTjh07pkgkogMH\nDngexyfhQZ+bm6szZ85E7585c6ZFDd8G/fv3lyT169dPjz/+uA4fPpzkEvmrX79++uSTTyTdqN3f\nPF5b9O3bV47jyHEcLViwIO3P37Vr1/TYY4/pqaeeiv7Lb9M5vHl8Tz75ZPT4bDuHX/ziF1VeXq7T\np09Hz5t0o4UkNzf3js9PeNCPHj1ax48f19mzZ3Xt2jVt27bNqq6Zzc3Nam5uliRdunRJe/bsUUFB\nQZJL5a/p06dr69atkqStW7dq+vTpSS6Rv25txnjttdfS+vwZY/Ttb39b+fn5Wrx4cXS5LeewveOz\n4RzW19frwoULkqTLly+rurpapaWlGjdunF5//XVJLs5dUFeLO7Jr1y5TUFBghg8fblauXJmMIgTm\n9OnTpri42JSUlJghQ4aYH/3oR8kuUlzmzp1rvvSlL5kuXbqY3Nxc85vf/KZF17yysrK07ZpnTOvj\n27Rpk5k3b54pLi42w4YNM1OnTjWRSCTZxfTsnXfeMY7jmJKSkhZdDW05h20d365du6w4h3/9619N\naWmpKSkpMffff7+prKw0xnjrXukYY8GlaQBAu5hhCgAsR9ADgOUIegCwHEEPAJYj6AHAcgQ9AFju\nfwHDB6YUp5ReYAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3bf6c10>"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I am sure that there is a representation of the negative binomial distribution as a hierarchical model:\n",
      "$$\n",
      "Y \\sim \\text{Poisson}(\\lambda)\n",
      "$$\n",
      "$$\n",
      "\\lambda \\sim \\text{Gamma}(\\alpha, \\beta)\n",
      "$$\n",
      "\n",
      "The PyMC docstring for pymc.gamma_like says the distribution is\n",
      "$$\n",
      "    f(x \\mid \\alpha, \\beta) = \\frac{\\beta^{\\alpha}x^{\\alpha-1}e^{-\\beta x}}{\\Gamma(\\alpha)}\n",
      "$$\n",
      "\n",
      "But since I want to represent the negative binomial as a Poisson with a Gamma-distributed rate, I need to write the distribution of Gamma in terms of $\\mu$ and $\\delta$.\n",
      "\n",
      "\n",
      "I believe that $\\alpha = \\delta$ and $\\beta = \\delta/\\mu$ is the appropriate transformation.  Let's check."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "lmb = mc.Gamma('lmb', alpha=delta, beta=delta/mu)\n",
      "Y = mc.Poisson('Y', mu=lmb)\n",
      "\n",
      "mc.MCMC([lmb, Y]).sample(10000, progress_bar=False)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# compare hierarchical model to closed form distribution\n",
      "hist(Y.trace(), range(30), histtype='step', label='Observed', normed=True)\n",
      "\n",
      "x = arange(0,30,.5)\n",
      "y = [exp(mc.negative_binomial_like(x_i, mu, delta)) for x_i in x]\n",
      "plot(x, y, drawstyle='steps-post', label='Predicted')\n",
      "\n",
      "legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "<matplotlib.legend.Legend at 0x3e6a610>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD9CAYAAABdoNd6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1clGW+P/DPPYBhOtDwuMiwaC0p4PCkZPwUo1RILFvM\nXTwkeTqVxuraupun9pQHxn0oa63sUBzPHltNsuSFJscFZzN9zQHWRCwJyfK3aSszgJKAzzA8zPX7\nw21+jsAwgwwDc3/er9e8ZOa+uO/v5S0fLq+557olIYQAERG5NYWrCyAiIudj2BMRyQDDnohIBhj2\nREQywLAnIpIBhj0RkQwMGPY6nQ4ajQZRUVHYsGFDr+3l5eVISEiAl5cXdu3a1Wv7pUuXoFar8fOf\n/3xoKiYiIofZDHuTyYScnBzodDrU1taiuLgYx44ds2oTHh6Obdu2ISsrq899rFu3Dvfdd9/QVUxE\nRA6zGfZVVVWIjo5GaGgoPD09kZmZidLSUqs24eHh0Gg0UCh67+qzzz5Dc3MzUlNTh7ZqIiJyiKet\njUajEWFhYZbnarUaer3erh2bzWY899xzeP/997F///4+20iSZH+lRERk4ejiBzZH9rcSxu+88w7S\n09MxYcIEm0UJIdz2kZub6/Ia2D/2T259k0P/BsPmyF6tVsNgMFieGwwGq5H+zW785XD48GFUVFTg\nnXfewZUrV9DZ2QmlUonf//73gyqUiIgGz2bYJyYmoq6uDg0NDQgKCkJRURE2b97cZ9ubf+MUFhZa\nvt62bRuOHj3KoCcichGb0zje3t4oKChAWloaYmNjsWjRIiQkJCA3Nxd79+4FAFRXVyMsLAzFxcVY\nsWIFNBpNn/uS4/x8SkqKq0twKvZv9HLnvgHu37/BkMRgJ4CG4uCSNOj5JyIiuRpMdtqcxiEi+fHz\n80NbW5uryyAAKpUKra2tQ7IvjuyJyAp/LkeO/s7FYM4R18YhIpIBhj0RkQww7ImIZIBhT0Sy98//\n/M9Yt24dAKCiogJTpkwZluMqFAqcPn16eI41LEchIhoCEydOxO233w6lUgmVSoWsrCxcvHjxlvcr\nSZLls0DJycn4+uuvB/yerVu3Ijk5+ZaPPVwY9kQ0akiShD//+c+4fPkyvvzySxw/fhy5ublWbbq7\nuwe1b3e/AolhT0Sj0oQJEzB//nx89dVXUCgUeOeddzB58mTLFMzOnTsxZcoU+Pj4ICEhAdXV1Zbv\n/fTTTxEZGQlfX18sWbIEHR0dlm16vd5qDbBTp04hPT0dvr6+8Pf3R05ODr7++ms888wz+PTTT6FU\nKuHn5wcAaG9vR05ODoKCgqBSqbBs2TK0t7db9qXVauHn54cf/vCHePfdd539V2SFYU9Eo8r3I3CD\nwYCysjLExMQAuH5Xvc8++wwnTpxAZWUlfv7zn6OoqAiXLl3Cc889h0ceeQQmkwkdHR3IyMjAs88+\ni4sXL2Lp0qXYvXt3n0u6dHV14cEHH0R8fDzOnz+Ps2fP4oknnsCUKVOwefNmJCUl4fLly5YPPj37\n7LNobm7GqVOn0NjYiEuXLuHXv/41AOCjjz7Cf/3Xf+Ho0aP45ptvUF5ePkx/Y9fxQ1VEZGWgn0tJ\nOzTrXIlcx3/2J06ciJaWFnh6esLb2xsPPfQQ3nzzTSiVSlRUVGDmzJkAgCeeeAJhYWFYv3695Xun\nTJmCt956CwDw9NNP48yZM5Zt999/P5KTk7F+/Xro9XpkZ2fDYDDg4MGDWLp0KRobG3vVsnXrVmzZ\nsgUVFRUAgM7OTvj6+uLLL7/EnXfeCeD6/yAeffRRNDY2IisrCxEREdBqtQCAM2fOYNKkSfjmm28s\n7W82lB+q4nIJROSQwYT0UJEkCSUlJXjggQd6bQsJCbF8bTQaUVRUhP/4j/+wvNbV1YWWlhaYzWaE\nhoZafa9are7zeE1NTZg4caJdtX333XcwmUyYNm2a5TUhhOU9hO+++w7333+/ZdvNNTgbp3GIyO2E\nhIQgLy8PbW1tlseVK1fwT//0TwgKCkJDQ4NV+xvv23GjCRMmWP0P4EY3T/v4+/vDy8sLf/vb3yzH\nvHDhAq5cuQIACAoKgtFotLS/8evhwLAnIrfz1FNPoaCgAMeOHQMAdHR04OOPP8aVK1eQnJwMk8lk\nuTdHaWkpDh8+3Od+kpOTMW7cOKxbtw6dnZ3o7OxEVVUVgOsLxjU1NaGrqwvA9SXhs7Oz8atf/QoX\nLlwAAJw9exYHDhwAACxevBhbtmzB6dOnYTKZrKaYhgPDnohGvZtH2bNnz8Zrr72GZcuWQalUIjw8\n3BLu3t7e2LVrF9588034+vpi+/btyMjI6HN/np6e2LdvH6qrqxEQEICQkBBs374dADB37lzceeed\n8Pf3R1BQEAAgPz8fKpUKkZGR8PHxwX333Ye6ujoAQEZGBp588klMmzYNERERSE5OHtb7fPANWiKy\nwp/LkYOrXhIRkUMY9kREMsCwJyKSAYY9EZEMMOyJiGSAYU9EJAMMeyIiGRgw7HU6HTQaDaKiorBh\nw4Ze28vLy5GQkAAvLy/s2rXL8vqxY8cwY8YMxMTEIDIyEu+9997QVk5ERHazGfYmkwk5OTnQ6XSo\nra1FcXGx5ePH3wsPD8e2bduQlZVl9bpSqURRURFqa2tx8OBBrF27Fi0tLUPfAyKiG+Tl5SE7O9vV\nZThkOO56ZXPVy6qqKkRHR1tWZ8vMzERpaSni4+MtbcLDwwFcv5fijX70ox9Zvg4JCUFYWBiam5vh\n7+9v1S4vL8/ydUpKClJSUgbVERoafhv80NbRZldblbcKrc+3OrkicjU/P6DNvn8Sg6JSAa0O/jPa\nunUrXn/9dXz77bcYP348MjIy8PLLL8PX13dYlyAYLnq9Hnq9/pb2YTPsjUaj1R1b1Gr1oA545MgR\nXLt2DZGRkb223Rj25HptHW12L2E7VOua08jW1gY4c/UER7N548aNeP311/Hhhx8iOTkZ586dw8qV\nKzFv3jz89a9/HZalHrq7u+HpOXwrxN88EP5+TXxH2JzGGYrfkE1NTXj88cexdevWW94XjSxjoYKk\nlex6eL3kh4MHXV0xjXaXLl1CXl4e/vjHP1qmPYKDg/H++++jqakJhYWFkCQJHR0dWLJkCXx9fREV\nFYUjR45Y9pGXl4fg4GAolUpERERYVqU0m81Yt24dQkND4evri4ULF+L8+fMAgL///e9QKBR49913\nMWnSJMydOxfp6el4++23reqLjY3Fnj17AAA1NTVITk6Gj48PwsPDrd63PHfuHObOnQulUol7770X\np06dcurfGzDAyF6tVlut82wwGKxG+je7+ZfDpUuX8NBDD+F3v/sd7rnnnlsslUaa+w63YvZswJ5T\nO7dSQk0N0Mc9J4jsdujQIZjNZsyfP9/q9dtuuw3p6enYv38/Jk+ejP/5n//B7t278eGHHyI/Px8Z\nGRn49ttvcfLkSbz77rv44osv8IMf/ACNjY3o7OwEALzyyis4ePAgampq4Ofnh9WrV+Opp56yhDdw\nfWr75MmTEEKgqKgImzdvxsqVKwEAJ06cQH19PRYsWIALFy4gLS0Nr7/+Oh577DGcOHECs2fPxtSp\nU5GQkIDly5djwoQJaG1txZkzZ/DAAw/YfZOUQRM2tLe3i/DwcGE0GkVnZ6eYPn26+Oyzz/psu2zZ\nMlFcXGx5bjKZxAMPPCDefPPNfvc/wOHJBZBn/zl58EEhysrs3+/GjYMsiobVzT+Xzv4xdWT/27dv\nF2FhYX1ue+mll8S8efNEXl6emD17ttU2tVotPvnkE/HNN9+IoKAgceDAAdHZ2WnVZtKkSeLAgQOW\n542NjcLDw0O0t7eLb7/9VkiSJIxGo2X7pUuXxLhx40R9fb0QQoh/+7d/E08++aQQQoitW7eK5ORk\nq/0vX75c/PrXvxbXrl0Tnp6e4vTp05ZteXl5YtasWb361F9GDiY7bU7jeHt7o6CgAGlpaYiNjcWi\nRYuQkJCA3Nxc7N27FwBQXV2NsLAwFBcXY8WKFdBoNACAoqIiVFRUYOvWrYiPj0d8fDxqa2ud+5uL\niNxaQEAAWlpa+pyXP3fuHAICAgD0vuWfWq3GuXPncNddd2Hjxo1Yt24dgoODsXjxYssdowwGAzIy\nMqBSqaBSqRAVFYUxY8ZYXUV4460PlUolFixYgA8++AAA8OGHH+Kxxx4DcP39zqqqKsu+VCoVduzY\ngba2NrS2tqKnp8fqVojDcotCh389DCEXH576gDyIwEBh18PLSwi93v79cmQ/Otz8czmSRvYXLlwQ\n48aNE3/+85+tXm9vbxdqtVps2bKlz5F9WFiY+OSTT6xeu3z5ssjKyhJLliyxtKmuru7zuN+P7Ht6\neqxe/+ijj0RcXJw4dOiQmDBhguX1LVu2iPT09D731dfIPjc317Uje5KnsjKgrm7gR2MjcN99rq6W\n5MTX1xe5ublYsWIFKioqAFwf0S9duhTBwcFYunQphBA4fPgwysrKAADvvPMOenp6kJycjL/97W+o\nqKhAd3c3xowZg9tuu81y2fjy5cvx4osvoqmpCQDQ1taGffv22awnPT0dZ86cQW5uLpYsWWJ5PSMj\nAzU1NSguLkZPTw/MZjOOHTuGkydPYuzYsUhPT4dWq0VnZydOnTqFP/3pT06/ZJRhT70EBABBQQM/\n/vE/ZnJzKtX1yyOd9VCpHKtn7dq1WL9+PX72s59BqVQiNjYWgYGBOHDgAMaMGQNJkvDII4/gvffe\nwx133IH8/Hzs3r0bY8aMQUdHB9asWQOVSoWAgAA0NjZaVgZ48cUXMWvWLMyYMQM+Pj5ISEhAeXm5\n5bh9hfGYMWOwaNEiHDhwwOqDpSqVCjqdDv/5n/8JPz8/+Pv7Y82aNejo6AAAbN68GQaDAf7+/njs\nscewbNmyQZwZx/C2hGRF0kr4dpnAUF8YIGklbFQK/PKXQ7tfGnr8uRw5eFtCIiJyCMOeiEgGGPZE\nRDIwfIs7ENGooFKp3HIxsdFI5ei71zYw7GnYbNoEFBXZ1/aJJ4AVK5xbD/Wt1dElKGlUYNjTsNm5\n0752JSXA0aMMe6KhxLB3c46sTw8ACtPQ/bfxZvfea1+748eBfyw2SERDhGHv5hxZnx7A9evrlzut\nHCJyEV6NQ0QkAxzZ07BQeavsvrPV7ZIKWeCbhERDiWEvAxcvAiaTfW17epxTgyP3quXtDomGHsNe\nBoKCAF9f+9redhvg4+Pceoho+DHsZUChAJqbXV0FEbkS36AlIpIBhj0RkQww7ImIZIBhT0QkAwx7\nIiIZYNgTEcnAgGGv0+mg0WgQFRVluTHvjcrLy5GQkAAvLy/s2rXLatu2bdsQHR2N6OhovPfee0NX\nNREROcTmdfYmkwk5OTmorKxEcHAwkpKSkJqaivj4eEub8PBwbNu2DX/4wx+svrepqQm/+c1vUFNT\nAwCIi4tDWloagoODndANIiKyxebIvqqqCtHR0QgNDYWnpycyMzNRWlpq1SY8PBwajQYKhfWu9u/f\nj/nz52P8+PEYP348HnzwQezfv3/oe0BERAOyObI3Go0ICwuzPFer1dDr9XbtuKGhAWq12up7jUZj\nr3Z5eXmWr1NSUpCSkmLX/omI5EKv19udvf2xGfbDcR/KG8OeiIh6u3kgrNVqHd6HzWkctVoNg8Fg\neW4wGKxG+je78ZeDo99LRETOYzPsExMTUVdXh4aGBnR1daGoqAjz58/vs60QAkL8/zsizZ07Fzqd\nDpcvX8bly5eh0+kwd+7coa2eiIjsYjPsvb29UVBQgLS0NMTGxmLRokVISEhAbm4u9u7dCwCorq5G\nWFgYiouLsWLFCmg0GgBASEgIXnzxRcyYMQMzZszAv//7v/NKHCIiF5HEjcPx4T64JMGFh5cFSSvB\n+xWB9nZXV2I/SSvhKaPAH//o6kqIRqbBZCc/QUtEJAMMeyIiGeCdqmhE+ugj4B8fvh7QAw8Afazk\nQUQ3YNjTiLRvn33tamqAbducWwuRO2DY04iUmGhfO5PJuXUQuQvO2RMRyQDDnohIBhj2REQywDl7\nGnFU3ipIWvsW4VN6qhCDVidXRDT6MexHIb8NfmjraLOrrcpbhVH04VkAQOvz9oe3vb8UiOSOYT8K\ntXW0QeTa91Hp9nbAL8+59RDRyMc5eyIiGeDIfpQ6eBDo7h64Ha9DJyKAYT9qPfIIkJRkX9uf/tS5\ntRDRyMewH6UiIoCPP3Z1FUQ0WnDOnohIBhj2REQywLAnIpIBhj0RkQww7ImIZIBhT0QkAwx7IiIZ\nYNgTEcnAgGGv0+mg0WgQFRWFDX3c1dlkMiEzMxMajQYzZ87EmTNnAABXr17F4sWLERkZiYiICOTl\n5Q158UREZB+bYW8ymZCTkwOdTofa2loUFxfj2LFjVm3y8/MREhKC48ePY+3atVi9ejUA4IMPPoCX\nlxe++uor1NbWYtu2bTh9+rTzekJERP2yGfZVVVWIjo5GaGgoPD09kZmZidLSUqs2ZWVlyM7OBgAs\nXLgQhw4dghACYWFhuHr1Knp6enD16lWMGTMGfn5+zusJERH1y+baOEajEWFhYZbnarUaer2+3zYK\nhQL+/v5obm5GWloatm/fjpCQEFy7dg1vvvkm7rjjjl7HuHF6JyUlBSkpKYPvDRGRG9Lr9b2y11E2\nw16SBn8XoMLCQrS3t6OpqQmtra1ITk7GnDlzMGnSJKt2nMsnIrLt5oGwVqt1eB82p3HUajUMBoPl\nucFgsBrpf9+mvr4eAGA2m9HS0oKAgABUVlYiIyMDHh4eCAwMxMyZM3HkyBGHCyQioltnM+wTExNR\nV1eHhoYGdHV1oaioCPPnz7dqk56ejsLCQgBASUkJkpKS4OHhgbvuugsHDx4EcP3KnE8//RQ/+tGP\nnNQNIiKyxWbYe3t7o6CgAGlpaYiNjcWiRYuQkJCA3Nxc7N27FwCwatUqNDY2QqPR4LXXXsNbb70F\nAFi5ciWuXLmCyZMnIy4uDllZWZg2bZrze0RERL1IQgj77lztjINLElx4+FFL0kqILxH4/HNXV+J6\nklbCzP0ClZWuroRo+AwmO/kJWiIiGWDYExHJAMOeiEgGGPZERDLAsCcikgGGPRGRDDDsiYhkgGFP\nRCQDNhdCIxrplJ4q/HWeBMmOdaFU3iq0Pt/q/KKIRiCGPY1qZfe14oUXYNcnaCXt4FdxJRrtOI1D\nRCQDDHsiIhngNM4IIQRw8eL1P8kx334LPPecHQ2VwKZNwLPPOr0kohGHq16OEHv2AD/9KTBu3MBt\nL/xCQsYXArt3O7+uke7yZeC//xvo6Rm47dqrEvzeEmhpcX5dRM40mOxk2I8Q778PlJVd/3MgklaC\nyOXfm6MkLcOe3AOXOCYioj4x7ImIZIBv0I4QK0774erdbdhh54eDiIgcwbAfIa6a25D1f4Vdc/ZE\nRI7iNA4RkQww7ImIZIBhT0QkAwx7IiIZGDDsdTodNBoNoqKisGHDhl7bTSYTMjMzodFoMHPmTJw5\nc8ayrba2FsnJyYiPj4dGo4HJZBra6omIyC42w95kMiEnJwc6nQ61tbUoLi7GsWPHrNrk5+cjJCQE\nx48fx9q1a7F69WoAQEdHB7KysrBlyxYcO3YMlZWV8PLycl5PiIioXzbDvqqqCtHR0QgNDYWnpycy\nMzNRWlpq1aasrAzZ2dkAgIULF+LQoUMwm83Q6XS45557cPfddwMAfH19oVBw1oiIyBVsXmdvNBoR\nFhZmea5Wq6HX6/tto1Ao4O/vj+bmZpw8eRKdnZ1ISUlBS0sLMjMz8dJLL/U6Rl5enuXrlJQUpKSk\nDL43RERuSK/X98peR9kMe0ka/J19enp6cOjQIRw9ehRjx47FnDlzMG3aNMyfP9+q3Y1hT0REvd08\nENZq7fio/U1szquo1WoYDAbLc4PBYDXS/75NfX09AMBsNqOlpQVBQUH44Q9/iNmzZ8PPzw9jx45F\neno6ampqHC6QiIhunc2wT0xMRF1dHRoaGtDV1YWioqJeI/P09HQUFhYCAEpKSpCUlASFQoE5c+bg\niy++QHt7O7q7u/G///u/iIyMdF5PiIioXzancby9vVFQUIC0tDSYzWZkZ2cjISEBubm5mD59Oh5+\n+GGsWrUK2dnZ0Gg0UCqV2LFjBwAgJCQEzz33HBITE9HV1YX09HT8+Mc/HpZOERGRNd68ZISQtBIX\nQnMy3ryE3AVvXkJERH1i2BMRyQDDnohIBhj2REQywLAnIpIB3paQZOXiRSAqyr62t98OlJdf/5No\ntGPYk2yovFVoWyfhKzvbSx0qXLnSyrAnt8CwJ9lofb7VofaSdvBrQxGNNJyzJyKSAYY9EZEMMOyJ\niGSAYU9EJAMMeyIiGWDYExHJAMOeiEgGGPZERDLAsCcikgGGPRGRDDDsiYhkgGFPRCQDDHsiIhlg\n2BMRyQDDnohIBgYMe51OB41Gg6ioKGzYsKHXdpPJhMzMTGg0GsycORNnzpyx2l5fX4/x48dj48aN\nQ1c1ERE5xObNS0wmE3JyclBZWYng4GAkJSUhNTUV8fHxljb5+fkICQnBzp07sWfPHqxevRolJSWW\n7b/85S+xYMEC5/VgBPPb4Ie2jja72o5TqJxcDRHJmc2wr6qqQnR0NEJDQwEAmZmZKC0ttQr7srIy\nvPrqqwCAhQsX4umnn4YQApIkYc+ePbjzzjsxbtw4J3Zh5GrraMMmlbCr7dGjQI+T6yEi+bIZ9kaj\nEWFhYZbnarUaer2+3zYKhQL+/v5obm7GuHHj8Oqrr+KTTz7Ba6+91u8x8vLyLF+npKQgJSXF8V6M\nYDt2APfcM3C7O+4AUlOdXw8RjT56vb5X9jrKZthL0uDuwSmEQF5eHtasWYPbb78dQvQ/ur0x7N3R\nQw8BL73k6iqIaDS7eSCs1Wod3ofNsFer1TAYDJbnBoPBaqT/fZv6+noEBQXBbDajpaUFgYGBOHLk\nCHbt2oV//dd/xYULF6BQKDB27Fj87Gc/c7hIIiK6NTbDPjExEXV1dWhoaEBQUBCKioqwefNmqzbp\n6ekoLCzE9OnTUVJSgqSkJHh4eKC8vNzSRqvVQqlUMuiJiFzEZth7e3ujoKAAaWlpMJvNyM7ORkJC\nAnJzczF9+nQ8/PDDWLVqFbKzs6HRaKBUKrFjx47hqp2IiOwkCVsT6s4+uCTZnM8f7SSthN94CM7Z\nj1KSVsK5HIGgIFdXQmRtMNlpc2RPJGdShwrBBfZdpKDyVqH1+VYnV0Q0eAx7on4EbGlFXR3sGtlL\n2sFduUY0XLg2DhGRDHBkT2RDdTXg52df254ewMPDufUQDRbDnqgfqanAb39rZ+MHgcOHgZkznVoS\n0aAx7In6UVhof1tJC3R3O68WolvFOXsiIhlg2BMRyQDDnohIBhj2REQywLAnIpIBhj0RkQww7ImI\nZIBhT0QkAwx7IiIZYNgTEckAw56ISAYY9kREMsCwJyKSAYY9EZEMcIljB124ADQ2uroKGmk8u1RI\n0UuAfuC2Pl4qnF/bCi8vp5dFZMGwd9AzzwCVlYCPjx2NM4HJk51eEo0Aq6624i9/sa/tV5kStm8H\n/uVfnFsT0Y0Y9g7q6gI2bQIefXTgtpIW+MlPnF8Tud4bb9jfVtICnZ3Oq4WoLwOGvU6nw9q1a9HT\n04Nly5bh+eeft9puMpnw+OOP48SJE/Dx8cGOHTsQHh6Ojz/+GC+88AK6u7shhMAf/vAHpKWlOa0j\nw2VvlB9217UBdQO3VXmrnF8QEZEdbIa9yWRCTk4OKisrERwcjKSkJKSmpiI+Pt7SJj8/HyEhIdi5\ncyf27NmD1atXo6SkBD/4wQ/wl7/8BYGBgfjyyy8xZ84cNDU1QZIkp3fKmbo821A8Vdg1siciGils\nXo1TVVWF6OhohIaGwtPTE5mZmSgtLbVqU1ZWhuzsbADAwoULcejQIQghEBMTg8DAQABAdHQ0zGYz\nOjo6nNQNIiKyxebI3mg0IiwszPJcrVZDr9f320ahUMDf3x/Nzc0IDg62tCkuLkZsbCzGjh3b6xh5\neXmWr1NSUpCSkjKIbhARuS+9Xt8rex1lM+yHYsrlxIkTeOGFF7B///4+t98Y9kRE1NvNA2GtVuvw\nPmxO46jVahgMBstzg8FgNdL/vk19fT0AwGw2o6WlxTJ9YzQakZGRge3bt2PSpEkOF0dEREPDZtgn\nJiairq4ODQ0N6OrqQlFREebPn2/VJj09HYWFhQCAkpISJCUlQaFQ4MKFC1iwYAFeeeUVJCUlOa8H\nREQ0IJth7+3tjYKCAqSlpSE2NhaLFi1CQkICcnNzsXfvXgDAqlWr0NjYCI1Gg9deew1vvfUWgOtX\n6Zw6dQrr169HfHw84uPjcf78eef3iIiIepGEEMJlB5ckuPDwgyJpJV56SbdE0kooCBZ45hlXV0Kj\n1WCykwuhERHJAJdLIHKBixeBs2fta+vvDy6aRreMYU80zMZChRc6JLyw2Y7G7Sq8ILXi5ZedXha5\nOYY90TC7lttqd1tJK6G9zYnFkGxwzp6ISAYY9kREMsCwJyKSAYY9EZEMMOyJiGSAYU9EJAO89JJo\nhKuqAn77W/vaJiYCbnD3T3IChj0Avw1+aOuw72Jmr27eV5aG15w5gD03eTt1CqioYNhT3xj2ANo6\n2jB2g4DCjkktSQICSwduRzQUVN4q/A723URofLQK/6fC/g9skbww7P/h6FEgKsrVVRBZa33esU/b\nEvWHb9ASEckAw56ISAYY9kREMsCwJyKSAb5BS+Qmxnuo8PH/kSBpB26r8lY59OYvjX4MeyI3UZzY\niuXLYdf9kd+AhLffBlaudH5dNDIw7IncxIwZwJo1gNlsR+PLwIYNDHs5cduw93rJD91edt7ip10F\nT7f9myC5uOMO4Be/sK/tr7SAENcf9pJ4Gf+o5rYR1+3VhuoFwq5/oGPGABERQ1+DXq9HSkrK0O94\nhGD/Rq/xDeNhfEqCYr2d39CuQoe2Fbfd5tSyhow7n7vBGjDsdTod1q5di56eHixbtgzPP/+81XaT\nyYTHH38cJ06cgI+PD3bs2IHw8HAAwMsvv4zt27fDw8MDGzduRGpqqnN60Y9p01w7GnH3f3Ds3+j1\nqwm/Ql4KsFfWAAAGVUlEQVRunt3tJa2EhgZg3LiB23p6Av7+g69tKLjzuRssm2FvMpmQk5ODyspK\nBAcHIykpCampqYiPj7e0yc/PR0hICHbu3Ik9e/Zg9erVKCkpwWeffYbdu3fj+PHjOHv2LGbNmoWT\nJ09izJgxTu8UEQ0tj04V7tpu58ipXYW67FZERzu3JnKMzbCvqqpCdHQ0QkNDAQCZmZkoLS21Cvuy\nsjK8+uqrAICFCxfi6aefhtlsRmlpKZYsWQIPDw+EhoYiOjoaR44cwaxZswZVqOoVP1ww2TkHDwDt\nXJ2SaKh0/87+yzQ9X/TD1GIJKB7aGsYpVPj6iVao1UO7X9kQNrz//vvimWeesTz/4IMPxIoVK6za\n3H333eLcuXOW55MnTxZNTU1i+fLl4sMPP7S8vmLFCvHBBx9YfS8APvjggw8+BvFwlM2RveTkCW/h\nyKUAREQ0aDaXS1Cr1TAYDJbnBoMBYWFhvdrU19cDAMxmM1paWhAYGNjre41GY6/vJSKi4WEz7BMT\nE1FXV4eGhgZ0dXWhqKgI8+fPt2qTnp6OwsJCAEBJSQmSkpLg4eGB9PR07Ny5E93d3TAajairq8M9\n99zjvJ4QEVG/bE7jeHt7o6CgAGlpaTCbzcjOzkZCQgJyc3Mxffp0PPzww1i1ahWys7Oh0WigVCqx\nY8cOAMC0adOQkZGBmJgYKBQKbN68GV5eXsPSKSIiuonDs/xDZN++fWLq1KkiMjJSvPLKK64qw2nC\nw8OFRqMRcXFxIjEx0dXl3LInnnhCBAUFialTp1pea2lpEXPnzhUajUakpqaKtrY2F1Z4a/rqX25u\nrggNDRVxcXEiLi5O7Nu3z4UVDl59fb1ITk4WU6dOFXfffbfYsGGDEMJ9zl9//XOX89fe3i6mT58u\n4uLiREREhPjFL34hhBDi9OnT4t577xVTp04VmZmZorOz0+Z+XBL2HR0dYuLEicJoNIquri4xffp0\n8fnnn7uiFKeZOHGiaGlpcXUZQ6a8vFx8/vnnVmG4atUq8cYbbwghhHjjjTfE6tWrXVXeLeurf3l5\neWLjxo0urGponD17Vhw/flwIIcTly5dFRESEqKmpcZvz11//3OX8CSHEtWvXhBBCdHV1iRkzZoiD\nBw+Khx56SHz00UdCCCGeffZZ8frrr9vch0vWs7/x+n1PT0/L9fvuRrjR1UbJyclQqaw/u1BWVobs\n7GwAwNKlS0f1Oeyrf4B7nMPg4GBMnToVADB+/HjExMSgoaHBbc5ff/0D3OP8AcDYsWMBAJ2dnejp\n6UFQUBAOHz6MH//4xwDsO38uCfubr8xRq9UwGo2uKMVpJEnCvHnzEBMTg/z8fFeX4xTfffcd/P/x\nufiAgAA0Nze7uKKh9/bbbyMyMhJLly5Fa+voX//973//O6qrqzFr1iy3PH/f9y85ORmA+5w/s9mM\nuLg4BAcH4/7774dKpUJAQIBle2ho6IAZ6pKwd/b1+yPB4cOH8fnnn+PAgQP405/+hE8++cTVJZGD\nVq5ciVOnTuHEiRO46667sHr1aleXdEuuXLmCxYsXY9OmTfDx8XF1OUPuypUr+MlPfoJNmzZBqVS6\n1flTKBSoqamB0WhEeXk59Hq94/sY+rIGZs/1+6NdUFAQACAwMBCLFy9GdXW1iysaeoGBgTh//jyA\n66P87/vsLgICAiBJEiRJwooVK0b1Oezq6sKjjz6Kxx57zPJff3c6f9/3Lysry9I/dzp/3/P19cWC\nBQtw+vRpy7kDrs+WqAdYR8IlYW/P9fuj2bVr13Dt2jUAwNWrV6HT6RDthqtC3fgZi8LCQqSnp7u4\noqF147TGrl27Ru05FELgySefRFRUFNasWWN53V3OX3/9c5fz19LSgsuXLwMA2tvbsX//fsTFxeHe\ne+/Fnj17ANh5/pz5DrItZWVlIjo6WkRGRorf//73rirDKU6fPi1iYmJEbGysiIiIEOvWrXN1Sbds\nyZIlIiQkRHh5eQm1Wi3effddq0v35s2bN2ov3ROid/+2bNkili5dKmJiYsSUKVNEWlqaMBqNri5z\nUCoqKoQkSSI2NtbqMkR3OX999a+srMxtzl9tba2Ii4sTsbGxYvLkyUKr1QohHL/0UhLCTd6uJiKi\nfrlkGoeIiIYXw56ISAYY9kREMsCwJyKSAYY9EZEMMOyJiGTg/wEHy84ylL98IwAAAABJRU5ErkJg\ngg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3bf6950>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    }
   ],
   "metadata": {}
  }
 ]
}