fonnesbeck/Mark-Recapture.ipynb

## Mark-Recapture.ipynb
{
  "cells": [
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Original BUGS model from Link and Barker 2009\n\n```\n# Data \nlist(T=2, Cells=4, f=c(NA,22,363,174), \nomegas=structure(.Data=c(0,0, \n1,0, \n0,1, \n1,1),.Dim=c(4,2))) \n\nmodel \n{ \n    for (t in 1:T) \n    { \n        p[t] ~ dunif(0,1) \n    }     \n    for (j in 1:Cells) \n    { \n        for (t in 1:T) \n        { \n            c[j,t] <- pow(p[t],omegas[j,t]) * pow(1-p[t],1-omegas[j,t]) \n        } \n        pi[j] <- prod(c[j,1:T]) \n    } \n    for (j in 2:Cells) \n    { \n        pi.obs[j] <- pi[j]/(1-pi[1]) \n    } \n    n <- sum(f[2:Cells]) \n    f[2:Cells] ~ dmulti(pi.obs[2:Cells],n) \n    pn <- 1-pi[1] \n    n ~ dbin(pn,N) \n    cN ~ dunif(n,10000) \n    N <- round(cN) \n} \n\n```"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "%matplotlib inline\nimport numpy as np\nimport pymc3 as pm\nimport theano.tensor as tt",
      "execution_count": 21,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# Number of observers\nT = 2\n# Observations by neither, both, first and second observer\nf = np.array([np.nan, 22, 363, 174])\n# Total observations\nn = np.nansum(f)\n# Observer indicators\nomega = np.reshape((0,0, \n                    1,0, \n                    0,1, \n                    1,1), (4, 2))",
      "execution_count": 14,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "with pm.Model() as Mt2:\n    \n    # Individual detection probabilities\n    p = pm.Uniform('p', 0, 1, shape=T)\n    \n    # Joint probabilities\n    c = p**omega * (1 - p)**(1-omega)\n    π = pm.Deterministic('π', tt.prod(c, 1))\n    \n    # Conditional (on observation) probabilities\n    π_obs = π[1:] / (1 - π[0])\n    \n    f_obs = pm.Potential('f_obs', pm.Multinomial.dist(n=n, p=π_obs).logp(f[1:]))\n    \n    # Latent population size\n    N = pm.Uniform('N', 0, 10000)\n    \n    # Probability of being observed by one or more observers\n    n_obs = pm.Potential('n_obs', pm.Binomial.dist(n=tt.floor(N), p=1-π[0]).logp(n))",
      "execution_count": 26,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "with Mt2:\n    trace = pm.sample(1000, tune=1000)",
      "execution_count": 30,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\n100%|██████████| 2000/2000 [04:26<00:00,  7.56it/s]/Users/fonnescj/Repos/pymc3/pymc3/step_methods/hmc/nuts.py:459: UserWarning: Chain 0 reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.\n  'reparameterize.' % self._chain_id)\n/Users/fonnescj/Repos/pymc3/pymc3/step_methods/hmc/nuts.py:451: UserWarning: The acceptance probability in chain 0 does not match the target. It is 0.715818227235, but should be close to 0.8. Try to increase the number of tuning steps.\n  % (self._chain_id, mean_accept, target_accept))\n/Users/fonnescj/Repos/pymc3/pymc3/step_methods/hmc/nuts.py:467: UserWarning: Chain 0 contains 16 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize.\n  % (self._chain_id, n_diverging))\n\n",
          "name": "stderr"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "pm.plot_posterior(trace, varnames=['N'])",
      "execution_count": 31,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 31,
          "data": {
            "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x1295292b0>"
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<matplotlib.figure.Figure at 0x129179ef0>",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcgAAAEYCAYAAADYn2bFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FNXCx/HvSS8kEJqUgGAUwdDFCwpeehNFUARBBUR8\nr/reqyhYuFIsXAURyysIKoogIMgVFORSBEVFsSFKU7zSpUNCCISElHn/mM26m0zIBlN24fd5nn2S\nPXNm9pzMZH87M2dmjWVZiIiIiLegsm6AiIiIP1JAioiIOFBAioiIOFBAioiIOFBAioiIOFBAioiI\nOFBAioiIOFBAigQwY8wuY8whY0y0R9lQY8yaMmyWyHlBASkS+EKAB8q6ESLnGwWkSOCbCIwwxlQo\n64aInE8UkCKB73tgDTCijNshcl5RQIqcH8YA/zDGVCnrhoicLxSQIucBy7I2Ax8Bj5V1W0TOFwpI\nkfPHWOBuoGZZN0TkfKCAFDlPWJb1GzAfuL+s2yJyPlBAipxfngKiC60lIoUy+sJkERGR/LQHKSIi\n4kABKSIi4kABKSIi4kABKSIi4iCkrBtQBBpNJCIixcH4Ukl7kCIiIg4UkCIiIg4UkCIiIg4UkCIi\nIg4UkCIiIg4UkCIiIg4UkCIiIg4UkCIiIg4UkCIiIg4UkCIiIg4UkHLeSs/MLtH6InJ+C6QvTA6Y\nhor/qPPYUp/r7hrfowRbIiJ+RPdiFREROVcKSBEREQcKSBEREQcKSBEREQcKSBEXjXoVEU8hZd0A\nEX8RERqsUa8i4qY9SBEREQcKSBEREQcKSBEREQcKSBEREQcKSBEREQcKSBE/pktPRMqOLvMQ8WO6\n9ESk7GgPUkRExIECUkRExIECUkRExIECUkRExIECUkRExIECUkRExIECUkRExIECUkRExIECUkRE\nxIECUqQEff3113Tr1o0KFSoQHR1No0aNmDdvnled9PR0Hn74YapXr05kZCRXX301n3/+eb5lnfh2\nEYf//SS/T76D3ROu5/jaOT63wxhT4GP8+PGO82RmZtKoUSOMMUyfPj3f9L1799KnTx/Kly9PbGws\nN910E3v27PG5TSL+TreaEykhS5cupXfv3gwYMIC5c+cSFhbG1q1bSU9P96p31113sXTpUiZOnMgl\nl1zClClT6Nq1K+vWraNp06bueqk/rSAoPIrIy1px8sdlRWrLunXr8pVNmTKF2bNnc8MNNzjO8/zz\nz3P06FHHaWlpaXTo0IHw8HBmzpyJMYZRo0bRvn17Nm7cSHR0dJHaJ+KPFJAiJSA1NZU777yT++67\nj5deesld3qlTJ696P/30E3PnzuWtt97izjvvBKBt27YkJiYyZswYFi9e7K5bY+irGBOElZNd5IBs\n1apVvrLbb7+dFi1akJiYmG/ajh07GDduHK+//jq33357vulvvPEGO3bsYNu2bVx66aUANG7cmMsu\nu4zXXnuNhx56qEjtE/FHOsQqpeqJJ57AGMMvv/xC165diY6Opnbt2syYMQOAd955h/r161OuXDna\nt2/P9u3bveZ/4403aNKkCREREVSuXJm77rqLpKQkrzqTJ0/m6quvpmLFiux5qR8HZg0nbft3XnWy\nUg6xe8L1pP64jONfzOb3yXdQoUIFDv/7SbJOOO81FcWCBQs4cuQIw4cPP2u9xYsXExoaSr9+/dxl\nISEh3HrrraxYsYKMjAx3uTHF9++6du1atm/fzqBBgxyn33vvvdx66620bt26wHa3atXKHY4AdevW\npXXr1nz44YfF1k6RsqSAlDJxyy230KNHDz744AOuvPJKhgwZwj//+U+mTp3K+PHjmTFjBtu2bWPA\ngAHueR577DHuu+8+OnXqxOLFi5k4cSLLly+ne/fuZGf/8TVPu3btYujQoSxYsIAqPR8hvPqlHPn3\nk5ze/n2+dqSsW0Bm8gEqdX+Al19+mYz92zj60fNedSwrBysnO98jKyvL62FZlnuetWvXUrFiRTZt\n2kSjRo0ICQmhVq1aPPnkk15t3bJlC3Xr1iUqKsrrNRMTEzlz5gy//fbbn/5bO5k5cyZhYWH0798/\n37Q5c+bw/fffM2HChALn37JlCw0bNsxXnpiYyNatW4u1rSJlRYdYpUw8/PDDDBw4EIAWLVqwZMkS\nXnvtNXbu3ElsbCwABw4c4IEHHmD37t1YlsXEiRMZO3YsY8aMcS+nXr16tGnThiVLltCrVy/APneW\nK2JlGhF1mpKZtJ/UH5cRmdDCqx0h5atSpefDAAwa1INhb3/O8TVvkZV6jJCYSgAc+8/LnNq8Ol8f\nQid6P58xYwaDBw8GYP/+/aSlpTFgwABGjx7NlVdeyapVq3j66ac5fvw4L774IgBJSUnExcXlW3bF\nihXd04tbeno6CxYsoEePHlSqVMlrWnJyMg899BATJkygcuXKnDx50nEZZ2t3cnJysbdZpCwoIKVM\ndO/e3f17XFwcVatWpVmzZu5wBKhfvz5gj5b8+eefycnJ4bbbbiMrK8tdp2XLlsTGxvL555+7A3L9\n+vWMHTuW7777jsOHjwD2nl1Ixfh87cgbmGFVLgYg+8QRd0BWaDOAmObX55t3yT/aeD2vW7eu+/ec\nnBzS09P517/+5T4f165dO44dO8aUKVN44oknKF++PJZlYYzJt2zPvdHi9sEHH5CSkuIOc08PP/ww\nCQkJ3HXXXYUup7TbLVLaFJBSJvLufYSFhTmWgb3Hc/jwYQCvc16ejh07Bthh2rFjR6644gpeeeUV\nHliyF4KCOf7FbDKP7c03X1BEjNdzExIKgJV9xl0WHFuF4JjK+eb1HGEKEBwc7P49d8+sc+fOXnW6\ndOnCtGnT2LJlC9dcc419ntTh0ojcvTB7T/KEY5/P1axZs6hSpYrXhxSAb775hrfffpvVq1eTkpIC\nwIkT9mufPn2a48ePU758eYwxxMXFOe7dJicnO+5ZigQiBaQEhNzAWblypeMbcO705cuXk5KSwnvv\nvUd8fDyP/LAUACsrI988vjqXQ6y5I0Pz7mXl7mEFBQW56y1atIi0tDSv85Bbt24lLCzM9YFg1zm3\nPa+DBw+ycuVK/v73vxMaGuo17eeffyY7O5t27drlm+/+++/n/vvvJzk5mQoVKpCYmMiWLVvy1du6\ndStXXHFFsbVXpCwpICUgdO7cmaCgIPbs2ZNvr8xTWloagNebf2bSPjJ+3+q4F+iLcznE2qtXL0aP\nHs3y5cu9BrOsWLGCiIgId1nPnj0ZO3YsCxYscI8ozcrKYv78+XTp0oXw8PBzanNBZs+eTXZ2tuPo\n1W7duvHpp596lR08eJD+/fszYsQIevToQbly5dztHjFiBDt27OCSSy4B7MFRX375ZYE3HhAJNApI\nCQgJCQk8+uij/P3vf2fbtm20bduWiIgI9u7dy8cff8zQoUNp3749nTp1IiQkhIEDBzJ8+HBOblrN\n8bVzCImtcs7nx0LKX0RI+Yvylbdo0cKhtq1hw4YMHjyYMWPGkJOTQ/PmzVm1ahXTp09n9OjR7qBp\n2rQp/fr1Y9iwYWRmZlK3bl2mTp3Kzp07mTPH+045GQf+S1bKIXD1I/PoXk79shawz6UGhUYA9o0H\nZs6c6XWuNtesWbNo1KgRzZo1yzetWrVqVKtWzats165dAFx++eVee5Z33303kydP5sYbb2TcuHEY\nYxg9ejS1atXib3/7W4F/F5FAooCUgPHMM8/QoEEDpkyZwpQpUzDGUKtWLTp27Mhll10G2Ics58yZ\nw5gxY+jZsyc55S4iru0gTu/8gfQ9m0q1va+99ho1a9bklVde4dChQ9SpU4cXXniBBx54wKvejBkz\nePzxxxk1ahTHjx+nSZMmLF++nObNm3vVS/3hI69DvWnb1pK2zQ7Imve8SVB5OyCzs7O9LiXJtWHD\nBjZt2uQ1yvdcRUdH88knn/Dggw9yxx13YFkWHTt25KWXXnKHv0igMwE06ixgGir+o85jS32uu2t8\njyLXLw3+2CaRAJd/CLYD3ShARETEgQJSRETEgQJSRETEgQJSRETEgQJSRETEgQJSAkJ6Zv7LFspa\nUdvkj30QkYLpOkgJCBGhwUW63AFK/pKHorZJl2CIBBbtQYqIiDhQQIqUEh1iFQksOsQqUkr88TCx\niBRMe5AiIiIOFJAiIiIOFJAiIiIOFJAiIiIOFJAiIiIOFJAiIiIOFJAiIiIOFJAiIiIOFJAi5xHd\nQF2k+OhOOiLnEd1AXaT4aA9SRETEgQJSRETEgQJSRETEgQJSRETEgQJSRETEgQJSRETEgQJSRETE\ngQJSRETEgQJSRETEgQJSRETEgQJSRETEgQJSRETEgQJSRETEgQJSRETEgQJS5AKm748UKZi+D1Lk\nAqbvjxQpmPYgRUREHCggRUREHCggRUREHCggRUREHCggRUREHCggRUREHCggRUREHCggRUREHCgg\nRUREHCggRUREHCggRUREHCggRUREHCggRUREHCggRUREHCggRaRE6TsnJVDp+yBFpETpOyclUGkP\nUkRExIECUkRExIECUkRExIECUoqFBmKIyPlGg3SkWGgghoicb7QHKSIi4kABKSIi4kABKSIi4kAB\nKWVCg3RExN9pkI6UCQ3qERF/pz1IERERBwpIERERBwpIERERBwpIERERBwpIERERBwpIERERBwpI\nERERBwpIERERBwpIERERBwpIERERBwpIERERBwpIERERBwpIERERBwpIERERBwpIERERBwpIERER\nBwpIERERBwpIERERBwpIERERBwpIERERBwpIERERBwpIERERBwpIySc9M7usmyAiUuZCyroB4n8i\nQoOp89jSIs2za3yPEmqNiEjZ0B6kiIiIAwWkiPhMh9/lQqJDrCLiMx1+lwuJ9iBFREQcKCBFREQc\nKCBFREQcKCBFREQcKCBFREQcKCBFREQcKCBFREQcKCBFREQcKCBFREQcKCBFREQcKCBFREQcKCBF\nREQcKCBFREQcKCBFREQcKCBFREQcKCBFREQcKCBFJKClZ2aXaH25cIWUdQNERP6MiNBg6jy21Of6\nu8b3KMHWyPlEe5AiIiIOFJAiIiIOFJAiIiIOFJAiIiIOFJAiIiIOFJAiIiIOFJABSNd9iYiUPF0H\nGYB03ZeISMnTHqSIiIgDBaSIiIgDBaSIiIgDBaSI+BUNKhN/oUE6IuJXNAhN/IX2IEVERBwoIEVE\nRBwoIEVERBwoIEVERBxokI4EpPTdGzn+xWzOHPoNExJGZMJVxLUfQnB0nLvOrl272D3hesf5az0w\nj6CIcgDkZKaT9PFrVHz9DlKtcCr8dSDRDf7qVT/lm39zassaqg9+GRMUXGj7Ds59DHJyqHb7c/mm\npf60gqTlr1DznjcJKX8RAEeXvsipzavddYIiYwmtVIvl7bxfy6s/QcEEhUcTWimeiDrNiGnavdB2\niYjvFJAScNL3bubQe6OJrNucKr1Gkn06leNfvMOheY9TfdDLmJBQr/qxrW4h6tKWXmUmLNL9+4mv\n/036rh+ZO/ttBj6/gKMfTSLsogRCK9YEIOvEUVK+mk/VW57yKRzPVVBUeareNBqA7FPJnPhuEddd\ndx1V+j5NZJ2m7nrRDTsR07QblmWRk36CjP3bSP3hI1LXL+Gr/nVLrH0iFxoFpASclC/fJSS2KlVu\nGuUOrNBK8Ryc9RAnN64kprn3sP+QCtUIr1m/wOWd3rGemCt70LNnTyp8FcypLWtI3/2TOyCTV79O\ndP1riYhvUHKdAkxQiFc7Iy5uwokZd5O6frFXQAbHVPKqF3VpS2KvvIGDcx7lpptuImzAFILCIkq0\nrRea9MxsIkJ9/3BU1PrinxSQEnAy9m8jOrG9195cePV6BEXGkvbrunwBWRgrOxMTEu5+HhQajpV1\nBrDDM33vZmrcPa14Gl8EQeFR1KtXjx93HCi0bnB0HHHth3Bo4Tgq/vw5MU26lEILLxy6NvPCpEE6\nEnhMECY4/2c7ExxC5tHd+cqPfzaT3c/1ZM+LfTn8/lOcObLLa3p4jcs5tXk1Bw4c4PSO9Zw5vJPw\nGpdjZWWStGoaFdoOIjgy9pyaauVkux9ZWVlYOdlgWT7Pu3fvXoLCo32qH1GnGSEhIWTs23pObRUR\nb9qDlIATWqkmGfu3eZVlpRwm+2QyBHvsVYaHU65pNyLrNCcoKpbMY7+T8vUCDs5+mOp3vEBo5VoA\nlG/dn8MLnqBGjRoAxP7lJsJrNuD42rkER5anXONz2xvL2LeVPRNv/KPdEwufx8qxb7OWfTKZlK/m\ncfLgQeI63ODT6wWFhlO5cmWOn0w+p/aKiDcFpAScmCt7cuyjSSR//g6xV95ATnoqx5ZPBmMw5o+D\nItWrV6dS17+7n0fUakjkJVey/837SFk3n8o3jAAgJKYy1e98hU/+pwHt/u9bgiNjyTx+kBPfLaLa\nbROwsjJI+mQ6ab9+jQkNJ/aqXsReWXhohVatS6Vu97ufL/lHG254ZS2n//s1Kevm56ufffKYV6Ca\nsEieeuoppp9smq9uQSzLAuNzdRE5CwWkBJxyie3JOvY7J75bxIl18wFDVINriUxoQeaR/IdYPYXE\nViE8/goyDv7Xq9wYQ0JCAsGRvwCQ/PE0yjXuQljVS0j+fBZnDv5GjbumkJ16jINzHyW0Ui3g7OeZ\ngkIjCa9+mft5ixYtCK9+iDOHdzjXj6pA1T5jwRiCI2MIjqnM6NE9edPHc185mRkcPXqUiKqNfap/\nodIAGvGVAlICUoW/3kFsq1vISjlIcFR5gqPj2PfGPYTHX1H4zBacbTcr7devOHN4B5V7PgJA+o71\nRDfqZL9OVHki6zQjfecPxdMRDyYo2CtQiyp95w9kZ2cT4cvf4AJW1AE3oEE3FyoN0pGAFRQWQViV\nOgRHx3F6x3qykn6nXNPrzjpP1onDZOzbSniNeo7TczLTSVr9BnEd7iYoPMpdbmWm/1HnzGmfB9qU\nluxTx0leM4Pq1asTlecmByJybrQHeQE43w4pnTm0ndM71hN2UQIA6b9v5cS37xPb8mavaxWHDx9O\n0ufbCa9Rn+Co8mQm2YN0MEGUb9XXcdkpX80jtGI80Q2udZdF1GlK6vqPCK0YT/bJJNJ3/0TsX3qX\nbCfPIjv1GBn7fnHdKCCVjP3bOPnTCsBi0bIV9Ft0tMzaJnI+UUBeAM67a7iCQji9/XtSvnkfsjMJ\nrVSLSl3+l3KNO3tVS0xMJGPuEk5tWkXOmdMERcUSUbsJFVr3J7RSfL7F/vLLL6T+sJTqg17yKi9/\nza1kn0rh2LKXMSFhVGg7mMi6zUu0i2dzavMqTm1e9cet5irGE3Pl9cQ07U7Lli1hUdEOH4qIMwWk\nBJywKhc73uM0ryFDhvDUrxf5vNz69etT+8EF+cqDwiKp3GNYkdpYbcD4AqfFNOlKTJOuXmWVezzo\n03IvfvSjIrVDRM6dzkGKiIg4UECKiIg4UECKiBSz9MzsEq0vpUPnIEVEitl5NzDuAqU9SBEREQcK\nSBEREQcKSBEREQcKSBEREQcKSBEREQcKSBEREQcKSBEREQcKSBEREQcKSBEREQcKSD+g20yJiPgf\n3WrOD+i2VCIi/kd7kCIiIg4UkCIiIg4KDUhjTDtjjOXwOO5R5+0C6ljGmF/yLK+gek19aXBycjLD\nhg2jdu3ahIeHEx8fz+DBg73qZGdn8+KLL9KwYUOio6OpXr06vXv3ZuPGjV716tSpgzHG8XHPPff4\n0hwRkSI5vf07Ds55lD0v9GHPi7dwYOYwPvnkEwBWr17N7bffTkJCApGRkSQkJHDvvfdy+PDhfMsp\n6L3rxx9/LLQNM2fO5Oabb+biiy/GGJPvPdTT1KlTqV+/PuHh4dSuXZvRo0eTmZnpWHfWrFlcddVV\nREVFERcXR5s2bdi0aZNvfxg/VJRzkPcD33k8z/L4/WlgWp76dYB3gcUOy3obeC1P2a+FNSA5OZk2\nbdpgjGHcuHHUqVOH/fv38+WXX3rVGz16NBMmTGDkyJF06NCBo0ePMm7cONq3b89PP/1EfHw8AIsW\nLSIjI8Nr3oULFzJx4kR69uwJ2ANoIkKDC2uaW1Hri8iFI/XHZSR9PI2Y5tdT/ppbwcrhzOGdpKWl\nATBt2jROnjzJqFGjuOSSS/jvf//L2LFjWbFiBRs3bqRcuXJeyxs8eDB/+9vfvMrq1atXaDtmz57N\nkSNH6Ny5MwsWLCiw3rPPPsvjjz/Ogw8+SLdu3fjxxx8ZO3YsBw4cYPr06V51//nPf/LSSy/xyCOP\n8Nxzz5GWlsa3337r7lsgKkpA/mxZ1tdOEyzL2g5s9ywzxnR2/TrTYZZ9BS3rbEaOHMnJkyfZtGkT\nsbGx7vJbb73Vq97bb79Nv379GDdunLuscePGNGjQgKVLl7o3qGbNmuV7jccff5xq1arRtWtXQANo\nRKR4ZKUcInn1G8S1G0LsVTe6yyMvuZLrr7ffN1599VWqVKninta2bVvq1atH27Ztee+99xgyZIjX\nMmvWrEmrVq2K3JYVK1YQFGQfQFy+fLljnfT0dJ555hkGDhzIpEmTAOjcuTPGGB555BEefPBBEhMT\nAVi3bh3jx49n4cKF9OrVy72MHj0C+/2wJM9BDgTWW5a1pTgWdurUKWbNmsXQoUO9wtHJmTNn8tWp\nUKECADk5OQXOt2fPHj799FNuu+02goO1Fygixefkxo/BGGKadS+wjmc45rrqqqsA2LdvX7G1JTcc\nz2bz5s2cPHmS7t2929utWzcsy+KDDz5wl02dOpW6det6heP5oCgBOccYk22MOWaMmWuMqV1QRWNM\na+BSnPceAe41xmQYY9KMMZ8YY64t7MXXr1/P6dOnueiii+jTpw+RkZGUK1eOXr16sXPnTq+69913\nH7Nnz+bDDz/kxIkT7Nixg/vuu4/4+Hj69etX4Gu88847WJbFoEGDCmuOiEiRpP++ldCK8Zz6+XP2\nvTaU3c/1ZN9rd5P6w0dnne+zzz4DoEGDBvmmTZ06lfDwcKKioujQoQNffPFFsbU3dychLCzMqzw8\nPBywAzTX2rVradKkCc899xw1a9YkJCSEhg0bnvXwbSDwJSBTgEnAUKAD9vnGTsA6Y0zVAuYZCGRi\nn4PMazZwn2sZ/wNUAj4xxrQ7WyP2798PwIgRIwgODmbx4sW8/vrrbNiwgXbt2pGamuqu+9RTTzFy\n5EhuuukmypcvT0JCAlu2bGHNmjVUrFixwNd45513aNasGY0aNTpbU0REiiz7ZBKZyftJ/vQtYlv2\noWq/p4mo05Skj6fx8ssvO86TmprKsGHDaNCgQb69s9tvv51XX32VVatW8frrr3Ps2DE6dOjAmjVr\niqW9l112GUFBQXz9tffZsHXr1gGQlJTkLtu/fz+rVq1i6tSpTJw4kWXLltGgQQP69u3Lhx9+WCzt\nKQuFnoO0LGsDsMGj6DNjzOfAt9gDd0Z51jfGhAN9gY8syzrqsLw7PJ5+YYz5ENgMjAPaFNSO3EOj\ndevWZd68eRhjAEhISKBVq1bMnj2be++9F7A/VY0bN45Ro0bRvn17jh49yvjx4+nSpQtffPEFNWrU\nyLf8r7/+mm3bthW4oYqI/ClWDtaZ01Tu9SBRl18DQOTFTchKOcyzzz7L/fff735fA8jKyqJ///7s\n27ePL7/8kpAQ77frd955x/37tddey4033kjDhg0ZNWoUa9eu/dPNLVeuHEOGDGHy5Mk0a9aMbt26\nsWHDBkaOHElwcLDXYdqcnBxSU1NZs2YNzZs3B6Bjx440btyYZ555hhtvvLGgl/Fr53QO0rKsH7BH\nnV7lMPlGoAIFH17Nu6xUYGkBy3KrVKkSAJ06dfLaiFq2bElsbCwbNtgZnpSUxIMPPsiIESN48skn\nadeuHX369GHlypUcOXKEiRMnOi5/1qxZhIaG0r9/f1+aLSJSJEGRMQBE1PG+oi2ybjMOHTrEgQMH\n3GU5OTkMGjSIVatW8cEHH9C4ceNClx8TE0OPHj347rvvCq3rq0mTJtG1a1cGDBhAXFwc1113HcOG\nDSMuLo7q1au761WqVImKFSu6wxHs85wdO3b06bITf/VnBukYwHIoHwQcBf5TDMtyyx0t5RmOnnI/\nzfz6669kZGS4T2znqlixIgkJCfz888/55s3IyGD+/Plcd911jifJRUT+rLDKFztPsOy3Ps89snvu\nuYf58+czb948Onbs6PNrWJZV4HvkuYiNjWXhwoUcOnSIjRs3cvjwYQYOHMjRo0dp0+aPA36JiYmO\nr1vc7Slt5xSQxpgWQD3gmzzlFwFdgLmWZTlfSZp/WbFAj7zLyis+Pp4WLVqwcuVKLOuPLF23bh0n\nTpxwB2K1atUA+Pbbb73mT0pK4rfffqNmzZr5lr1kyRKSkpKKZXCObjwuIk4i610NwOmdP3iVn975\nA/Hx8e73ruHDhzN9+nRmzJhRpFGhJ06cYOnSpbRs2dJxelHfmzzrV6lShUaNGhETE8OLL75I5cqV\nueWWW9zTe/fuzbFjx/j+++/dZTk5OaxatSrfzkogKfQcpDFmDrAT+AE4DjQDRgL7gFfyVL/NtUzH\nw6vGmBHA5cCnwH7gYmAEUM01r2fd34DdlmW5Pz6NHz+erl270qdPH4YOHcqRI0d4/PHHqV+/PgMG\nDADsu+Ncf/31TJw4kaCgINq2bcuxY8d47rnnyMjIcJ+n9DRr1iwqVapULNfsFPW6SdC1kyIXgshL\nWhBeuzFJK6aQc/oEIRWqkbbtS9J3beDpGTMAmDBhAi+88AJDhgzhsssu8xogU6VKFRISEgB4/vnn\n2bZtG+3bt6dGjRrs3r2b559/noMHDzJnzhyv17300ku5+OKLWb16tfu96czRPWQe3QNA0rETzPv0\nB5b2GglARO1GBEeVZ9f4HsyfP5+kpCQuv/xykpOTWbRoEfPnz+f9998nJibG/Rp33XUXU6ZM4eab\nb2bcuHFUrlyZ119/nW3btrFy5cqS+6OWMF9uFLAZ6A/8A4gCDgILgbEOg3AGAZtd5yidbAN6ux7l\ngRPAl8BdlmV9m6duCOB1MWLHjh1ZsmQJY8aMoXfv3kRHR9OjRw8mTpxIZGSku978+fOZNGkS7777\nLpMmTSKt9J+GAAALdElEQVQ2NpbmzZuzdu1aWrRo4fUiR44cYdmyZdxzzz35hjOLiBQXYwxVbxpF\n8mczOb52LjnpJwmtFE/lG0a4b/W2bNkyAN566y3eeustr/kHDRrE22+/DcDll1/OokWLWLRoESkp\nKcTGxtK6dWvefPNN/vKXv3jNl5WVRXa2995j2i9fkPLlHxcZZOzZRMYe+5ZwF/V/huDajd1tfvXV\nV9m+fTshISG0atWKNWvW0Lp1a6/lRUREsHr1ah5++GEeeOAB0tLSaNasGcuWLSvSIWJ/48so1meB\nZ31ZmGVZTQqZvgRY4uOy6jiVd+/ePd+Fq3lFRUUxevRoRo8enW9a3lvBValSpcD7CoqIFKeg8Cgq\ndbmXSl28j2Tlvi/5eolG527XccMNN/hUd9euXfnKKrS5jQptbstfOU+b+vbtS9++fX16nbjKVZk9\ne7ZPdXOX7++35bzgvg9St44TEX9zLu9LJX0qp6TbFAjvrfq6KxEREQcKSBEREQcKSBEREQcKSBER\nKXXncs14aV9nfsEN0hERkbIXCNeMaw9SRETEgQJSRETEgQ6xnsXuCddjJhRtnqLWP5d5/K1+abzG\n+dAm9SEw6xfna1z86Nm/HFn8i/YgRUREHCggRUREHCggRUREHOgc5Flc/OhH53R/wXMZulzS9zy8\n0Prgj21SHwKzfmm9hvgf7UGKiIg4MJZllXUbfPLkk0/+hv2FzeejGthfIH0+Ut8Ck/oWmNQ33xwd\nO3Zst0JrWZYVEI8nnnjCKus2qG/qm/oW+A/1LTAfZdE3HWIVERFxEEgB+WRZN6AEqW+BSX0LTOpb\nYCr1vgXMOUgREZHSFEh7kCIiIqVGASkiIuJAASkiIuKg1APSGNPOGGM5PI571Hm7gDqWMeaXPMuL\nMMZMNMYcMMacNsasM8b8tbT75WpLoX1z1Us0xiw0xuw3xpwyxmwxxgw3xoTkqRdkjBlpjNlljEk3\nxvxkjLm5dHvlbouvfWtqjFlujDlpjDlhjFlsjLnUYXl+s9482nSdMeZzj7Z/b4zp4DE9zhgz3Rhz\n1LXeVhljGjksJ6D6ZoyJMcY8b4xZ45pmGWPaFbAcv9kmPdp0tr51NMbMNsZsd62L7caYqcaYqg7L\nCbT1dqXrf22fa10cNMb8xxhztcNyfNp2S1Nh/2956r7m2i5nO0wrufVW2teVAO0AC/gH0Mrj0cKj\nTkKeaa2AW13zPZdneXOwbyBwN9ARWAicBpr6ad9qAEeAH4G+QAfgX0AOMCHP8v4FZAAjgPbAa656\n1/lp3y4DUoF1wI3AzcBG7It7q/rrenO1529AJvAi0BnoCjwKXO+aboAvgN+B/kA34DPgKBAf4H2r\nAyQBq4D3Xeu5XQHL8ptt0se+LQCWAXcCbYGhwD5gB1AuwNdbR2Ay9ntjW+AW4CvgDPAXj+X4vO36\nS9/y1L0GOAmkALMdppfYeiuLP0w71z9gpyLON9o1X6JHWRNX2Z0eZSHANmCxP/YN+B9XnXp5yucB\nBzyeV3W9ET2Zp95qYKOf9m26a0Ot4FEWD6Tj8cHGD9dbHdc/1LCz1LnR1eb2HmXlsYPl/wK8b8bj\n904UEJB+uE360rcqDmV/dfVxSCCvtwLmi3Gto1c8ynzadv2xb0AosBkYCewiT0CW9HoLpHOQA4H1\nlmVt8Sjrif0pZH5ugWVZWdhh09UYE166TfRJmOvniTzlx/E+5N3VVTfvIYXZQCNjTN2Sad6f0gpY\nZ1mW+7CrZVm/Y2/gvT3q+dt6G4K9FzTtLHV6Avsty/o0t8CyrBRgCfYbkGe9gOqb5XpX8YG/bZO+\n9O2IQ/F3rp81PcoCbr0V4BR2QGZ6lPm67ZaWovTtYSAYmFTA9BJdb2UZkHOMMdnGmGPGmLnGmNoF\nVTTGtAYuBWbmmZQI7LQsKy1P+Rbsf+R8575Kydn6tgD70MZkY0xdY0ysMaY3cAfeG0Ei9ob+W55l\n535AuKKkGl+Is/UtG/vwTl4ZQIIxJsL13N/WWxvgF+BW1zmqLGPMb8aY//Wok4gd9HltAWobY8p5\n1Au0vvnK37bJc+1bW9fPnz3KAna9uc4Lh7r+Fye7iqd7VPF12y0tPvXNGJMAjALusyzL6X0FSni9\nlcXXXaVgB8Fn2HtRzYB/AuuMMc0syzrsMM9A7E8J7+YprwgkO9RP8phemgrtm2VZh1wn0T/EPg8C\n9iGCJyzLes5jWRWB4w6f7v22b9iHNa4xxoRalpUJ9gAQ7I3YAHHAAfxvvdVwPSZi92k79vmcycaY\nEMuyXna1aZfDvLltjsM+TxKIffOVv22TRe6ba3t8CTscP/CYFMjr7T3s8/0Ah7HPB2/1mO7rtlta\nfO3bNGCh556vgxJdb6UekJZlbQA2eBR9Zoz5HPgWuB/7E4Obaxe5L/CRZVlH8yzOYIdLXqb4Wuw7\nX/pmjKmCfRL5FNAHOIY9UGeUMSbDsqwJrnkDrm/Ay9gb+jRjzBjs7WsSkPsJNcf106/6hn0kJQYY\nbFnWQlfZJ8aYOsBIY8z/4XubA65vRTjEGtB9M/Yo8XexD622dh2Kc08mcPv2CDABqAX8L/CRMaaT\nZVnfu6YHXN+A24CrgPqFLKtE++YX5yAty/oB+BX7D5LXjUAF8h9eBftTgtMnhDiP6WXKoW+PYJ+k\n7mpZ1vuWZa2xLGsM9qepp40xlV31koA4Y0zeFe23fbMs60vsf9A+2CPmdvHHujvDH232t/V2zPXz\n4zzlK4GLgOoU3ubcT7GB2Ddf+ds26XPfjDFB2NthJ6CXZVkb88wTsOvNsqwdlmV95wqb7th7keM8\n5vF12y0thfWtFvACduinG2MqGGMqYOdVqOt5qGueEl1vfhGQLgV9EhiEfc7uPw7TtgB1jTFRecqv\nwH5DznuupKx49q0R8JtlWXk3ym+xR2zlHjPfAoRjX/LiKfc8z1b8g9d6syzrVezRjg2B2pZldcI+\nnPJN7mFX/G+9bSmgPDcIclx1Eh3qXAHssSwr9xBVIPatKMvyp22yKH2bBvQDbrUsa3UBywr49eY6\nV7cR73Nvvm67paWwvtUAqgDPYId37qMW9tHEZKCHx7JKbL35RUAaY1oA9YBv8pRfBHQB5nq8uXpa\njB0qt3jME4L9j7DSsqyMEmu0jxz6dhC41BgTl6dqS9fPfa6fy7FX8G156t0ObLYsa2cJNLdIClpv\nlmVlWJa1xbKsvca+GLkTMNWjir+tt0Wun13zlHcFfrcs6yB2m2saY3IHeGCMiQVucE3LFYh985W/\nbZM+9c0YMwn7+sc7Lcv6AGfnxXpzBUUL7PN6uXzddkvLWfuGfY14e4fHIexrddsDa13zlOx6+7PX\niRT1gX1R5zjgJuxzb8Ox9xD3AJXz1H0Ie++k+VmWNw/7E8VQ7ItE/4193V2B85Rl37AvhcjEHmre\n19Xmp7HfeBbmWd54V18ewr4OcSr2J8cb/LRv8dgXkvfADsVHsQf3zPfz9WaAT7AP/dyD/aHsdde2\nN9hVJwj7Iuy92BdmdwXWYB/CqRXIfXPV6459aPxJ17Sxrufd/Xib9GW9Pep6/ib5bz6SEMjrDfsm\nDc9iD9Bpiz2Y8VvskcbXeizL523XX/pWwHy7cL5RQImtt1L9w7g6MxL7EEAKdlDsdf1xqjvU/QnY\nVMjyIrGPVx90/VG+oYC7gPhL31z/nP/BHtF5CvswwSggMk+9YFf5btdGvxHo4699wz5/sAo7ODOw\nD7kNB0L8eb252hMLTMH+lJp7mGpAnjoVgbdcbyxp2BfINzlP+rbL9QaV97HLX7dJX/qGHQRO/bKA\ntwN5vWFfT/gVdtCkY+81zgUaOSzLp23XX/pWwDy7cA7IEltv+j5IERERB35xDlJERMTfKCBFREQc\nKCBFREQcKCBFREQcKCBFREQcKCBFREQcKCBFREQcKCBFREQc/D9gqtPXGYy7jgAAAABJRU5ErkJg\ngg==\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "pm.summary(trace, varnames=['p'])",
      "execution_count": 35,
      "outputs": [
        {
          "output_type": "stream",
          "text": "\np:\n\n  Mean             SD               MC Error         95% HPD interval\n  -------------------------------------------------------------------\n  \n  0.325            0.020            0.002            [0.290, 0.367]\n  0.894            0.022            0.002            [0.853, 0.936]\n\n  Posterior quantiles:\n  2.5            25             50             75             97.5\n  |--------------|==============|==============|--------------|\n  \n  0.287          0.311          0.323          0.337          0.365\n  0.853          0.879          0.893          0.910          0.937\n\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "pm.forestplot(trace, varnames=['π'])",
      "execution_count": 33,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 33,
          "data": {
            "text/plain": "<matplotlib.gridspec.GridSpec at 0x127ae55f8>"
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<matplotlib.figure.Figure at 0x129a50d68>",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEiCAYAAACobJZDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF6pJREFUeJzt3X2w5FV95/HPZ2Yw+MBCVpItFMw4o91khlJc0cUilTTG\nSrHagzGhElNmV1IJPeviItqY3Y3EdeMaE2k0y+ImEnXNJiTiQ2rlN0m0TAKmzAZQzPAwg93KLBVY\nXCM+gARGefjuH/2bS+d679y+9/7uPX1Ov19VXbcfzj2/76G5/Z3v75zza0eEAACYVVtSBwAAwNGQ\nqAAAM41EBQCYaSQqAMBMI1EBAGYaiQoAMNNIVMiG7TfYvt32AdsXTzz/Ntv/1/b++vby+vmzbN9q\n+3O2n1M/d4LtT9n2Msc4xvZv2P5SfaybbP/LdcS83fbt9f0zbF8xEfMlR2u/imM8OEWbi20/ZTX9\nroXtD9k+b6OPg/lCokIWbJ8m6QJJL5b0fEld28+daPKeiDi9vv1p/Vxf0k9L+hVJr6uf+1VJvx7L\nbyB8u6STJJ0WEadJ2iPpuCXi2braMUTE5yPiotX+XkMulrSqRLWWMQIbgUSFXPywpBsi4qGIeFTS\nZyS9aoXfeUTSkzX+gH7E9k5Jz4yIzyzVuK44LpD07yLiO5IUEV+NiI/Urz9o+9ds3yjpJbZfaPsz\ntm+uq7ST6nYvtH2L7b+RdOFE/x3b+yYO+Xzbf1lXbxcsEc9W25fVFeGttvcebbB1/9fb/pjtL9q+\n2mMXSXqGpOtsX1e3/Qnbf2P7C7Y/avtp9fN32X6r7c9K+mXbN030v932rfX9t9Zx3W77qqUq1Loy\nPVjHPjha7MDRkKiQi9sl/ajtp9cJ5eWSTpl4/fX1B+IHbX9//dw7JV2lcTVxpaR3aFxRLec5kv4u\nIh5Y5vWnSro9Iv6FpBsl/TdJ50XECyV9sO5fkv6HpIsi4iUrjOl5kl4h6SWS3mr7GYte/0VJ90fE\niyS9SNIFtp+9Qp8v0Hi8uyTtkHRWRFwh6V5JZ0fE2bZPlHSppJdFxD+X9HlJb5ro43BE/EhEvFPS\nk2zvqJ//WUkfqe9fGREvqqvOJ0vqTgZh+59q/A+J3RHxPEn/ZYW4gWWRqJCFiLhD0m9K+rSkT0q6\nRdKj9cu/LWmnpNMlfUXS5fXv7I+IMyPibI0/tO+VZNvX2P4D2/9slWE8Junj9f22pNMkfdr2fo0/\n+E+2fbykEyaqtt8/Sn+fiIiHI+I+SddpfFpz0k9I+td1/zdKerqk5+roboqIeyLicUn7JW1fos2Z\nGieyv677fq2kH5p4/ZqJ+x+R9DP1/Z+deO1s2zfavk3SSyXtXnSMByQdlvR+2z8l6aEV4gaWtS11\nAMC0IuIDkj4gSbZ/XdI99fNfPdLG9u9Kmjy9pvq01KUaf9BeKek/afwBfpGkt0w0/bKkZ9k+LiK+\nvUQIhyPisSPdSjqwuGqyfYKkaS+gubjd4sfW+DTkp6bsT5K+M3H/MS39N25Jn46In1umj3+YuH+N\npI/a/mNJERFfsn2spP8u6YyIuNv22yQdO9lBRDxq+8WSflzSqyW9XuOEBqwaFRWyYfsH65/PkvRT\nkv6ofnzSRLNXaXyacNJrJf1JRHxT4/mqx+vbP1pcEBEPaZwIr7D9pCN92/75JcIZSvoB2y+p2x1j\ne3dEfEvS/bZ/pG73mqMM6ZW2j7X9dEkdSZ9b9PqnJL3O9jH1MVq2n3qU/o7m23piUcgNks7yEysh\nn2K7tdQvRcSdGie8X9UT1dSRpHRfPbf1Pav86uePrxe2XKxxtQusCRUVcvLx+kP9EUkX1olHkt5l\n+3SNK5K7JC0sOqjns16r8Wk0SXq3xqfvvitpqYriUo3nUw7aPqxxdfHWxY0i4rv1Muwr6tN92yT9\nlqQDkn5B0gdtP6RxslnOTZL+RNKzJL09Iu61vX3i9fdrXPl9oa4KvybpJ4/S39FcJenPbH+lnqc6\nX9If2f6++vVLJY2W+d1rJF0m6dmSFBHfqivX2zT+7704wUrjpPiJuvqypDeuMW5A5ms+AACzjFN/\nAICZRqICAMw0EhUAYKaRqAAAM41EBQCYaTO7PP2cc86J++67T5J0+PBhSdKxxx57tF8BAGTk5ptv\n/lREnLNSu1lenr4QWKfTkSRdf/31iUIBAGyAJb9uZzFO/QEAZloWFRUAoEhUVACA/GWRqHq9nnq9\nXuowAAAJzOyqv0mj0XLXygQAlC6LRNXv91OHAABIhMUUAIBUyllMUVWVqqpKHQYAIIEsKio2/AJA\nkaaqqLKYo2q1lvyWbADAHMiiogIAFKmcOSoAwPzKIlF1Op2FeSoAwHzJIlEBAOZXFnNUw+FQktRu\nt5MFAwBo3FRzVFkkKgBAkcpZTDEYDDQYDFKHAQBIIIuKig2/AFCkcjb8drvd1CEAABLJoqICABSp\nnDmq4XC4sPIPADBfsqiomKMCgCKVU1EBAOZXFhUVAKBIVFQAgPxlkah6vZ56vV7qMAAACTSSqGw/\nZnu/7WfUj99h+27bDy5q90bbf2f7ytX0PxqNNBqNmggVAJCZpjb8PhwRp088riRdKelLk40i4j22\nvynpjNV03u/31x/hBjl0SNqzRxoOpXZbqippx47UUQFAOY66mML2FZIulHSbpK2S2pLeEhGXLWr3\nYEQ8bYnf/57nbZ8v6YyIeP0Ksc3cYgpPNe33vWZ3vQoAJLX+xRQRcZGeqJbOlnTv4iS1GaqqUlVV\nG3oMe+XbRva9nv4BoGRZXOvv8ssvlyTt2bNnw46x1qpn927pi1+UHn9c2rJFOvVU6cCBZmMDgHmW\nxaq/VqulVquVOowlVdU4OW3dOv65wYUfAMydFTf8Hplnsn2CpIMR8Yzl2kzzfM5zVACARjW+4fd+\nSV+zfdWKR7bfZfseSU+xfY/tt63iOAAALGjkEkrLVVTLtD1fq6youCgtABRpUy+h9MDkht/l2H6j\npP8o6YGGjgsAKFwWF6U98l1U7XY7WTAAgMZNVVFlkagAAEUq5+rpg8FAg8EgdRgAgASyqKhYTAEA\nRZqqosriyhTdbjd1CACARLKoqAAARSpnjmo4HC6s/AMAzJcsKirmqACgSOVUVACA+ZVFRQUAKBIV\nFQAgf1kkql6vp16vlzoMAEACWeyjGo1GqUMAACSSRaLq9/upQwAAJMJiCgBAKuUspqiqSlVVpQ4D\nAJBAFhUVG34BoEjlXJS21WqlDgEAkEgWFRUAoEjlzFEBAOZXFomq0+kszFMBAOZLFokKADC/spij\nOvJdVO12O1kwAIDGTTVHlUWiAgAUqZzFFIPBQIPBIHUYAIAEsqio2PALAEUqZ8Nvt9tNHQIAIJEs\nKioAQJHKmaMaDocLK/8AAPMli4qKOSoAKFI5FRUAYH5lUVEBAIpERQUAyF8WiarX66nX66UOAwCQ\nQBb7qEajUeoQAACJZJGo+v1+6hAAAImwmAIAkEo5iymqqlJVVanDAAAkkEVFxYZfAChSORelbbVa\nqUMAACSSRUUFAChSOXNUAID5lUWi6nQ6C/NUAID5kkWiAgDMryzmqI58F1W73U4WDACgcVPNUWWR\nqAAARSpnMcVgMNBgMEgdBgAggSwqKjb8AkCRytnw2+12U4cAAEgki4oKAFCkcuaohsPhwso/AMB8\naSRR2d5u+2Hb+22fYvs623fYPmD7DRPtLrP9/2xfspr+9+7dq7179zYRKgAgM01WVHdGxOmSHpXU\nj4gflnSmpAtt75KkiHizpN9p8JirduiQtHu3tG3b+OehQymjAQCspPHFFBHxFUlfqe9/2/Ydkp4p\n6eBa+2xitZ+XOBN68KC0c6c0u9N0AIANXfVne7ukF0i6cSOPs/zx196O5AUAs2HDFlPYfpqkj0u6\nOCIeWE9fvV5PvV5v1b8X8b23XbukLfWot2wZP16qHQBgNmxIorJ9jMZJ6uqI+OP19jcajTQajdYf\nmKSqkk49Vdq6dfyTb7gHgNnW+Kk/25b0AUl3RMS7m+iz3+830Y0kaccO6cCBxroDAGywjZijOkvS\nv5J0m+399XO/EhF/utYO9+zZ00hgAID8bMSqv89qyt3G06rq83MkLACYP41cQsn2KZL+t6Sv13up\nlmt3maRXSbo8In57hW65KC0AlG3zLkobEXdLOmWKdm+W9ObV9t9qtdYSFgCgAFyUFgCQSjkXpQUA\nzK8sElWn01mYpwIAzJcsEhUAYH5lMUd15Luo2u12smAAAI2bao4qi0QFAChSOYspBoOBBoNB6jAA\nAAlkUVGx4RcAirR5G343WrfbTR0CACCRLCoqAECRypmjGg6HCyv/AADzJYuKijkqAChSORUVAGB+\nZVFRAQCKREUFAMhfFomq1+up1+ulDgMAkEAW+6hGo1HqEAAAiWSRqPr9fuoQAACJsJgCAJBKOYsp\nqqpSVVWpwwAAJJBFRcWGXwAoUjkXpW21WqlDAAAkkkVFBQAoUjlzVACA+ZVFoup0OgvzVACA+ZJF\nogIAzK8s5qiOfBdVu91OFgwAoHFTzVFlkagAAEUqZzHFYDDQYDBIHQYAIIEsKio2/AJAkcrZ8Nvt\ndlOHAABIJIuKCgBQpHLmqIbD4cLKPwDAfMmiomKOCgCKVE5FBQCYX1lUVACAIlFRAQDyl0Wi6vV6\n6vV6qcMAACSQxT6q0WiUOgQAQCJZJKp+v586BABAIiymAACkUs5iiqqqVFVV6jAAAAlkUVGx4RcA\nilTORWlbrVbqEAAAiWRRUQEAilTOHBUAYH5lkag6nc7CPBUAYL5kkagAAPMrizmqI99F1W63kwUD\nAGjcVHNUWSQqAECRNm8xhe3tth+2vd/2sbZvsn2L7QO2//NEu6ttf8P2eavpfzAYaDAYNBEqACAz\nTc5R3RkRp0v6jqSXRsTzJZ0u6RzbZ0pSRLxG0rWr7Xjfvn3at29fg6ECQNkOHZJ275a2bRv/PHQo\ndURr1/iG3xifS3ywfnhMfVvXabxut7vesACgOJ7qxJl08KC0c+fyr8/uDNDYhlyZwvZWSTdLeo6k\n90bEjevp75JLLmkkLgDIwbQJKOXxNjO5bcjy9Ih4rD4NeLKkF9s+bT39DYfDhZV/AFC6iPXfdu2S\nttSf8Fu2jB830e+R22ba0H1UEfEtSddLOmc9/ezdu1d79+5tJCYAmAdVJZ16qrR16/hnzl9A0fip\nP9s/IOmRiPiW7SdLepmk32z6OACA5e3YIR04kDqKZmzEHNVJkn6vnqfaIukjEbGuJXt8vQcAzK+N\nWPV3q6QXNN0vAGA+NTVH9Zik423vP1oj21dL+jFJh1fTea/XU6/XW0d4AIBcNVJRRcTdkk6Zot1r\n1tL/aDRay68BAAqQxTf89vv91CEAABLhorQAgFTK+YbfqqpU5bwJAACwZllUVEe+3Zdl6gBQlKkq\nqizmqFqtVuoQAACJZFFRAQCKVM4cFQBgfmWRqDqdzsI8FQBgvmSRqAAA8yuLOaoj30XVbreTBQMA\naNxUc1RZJCoAQJHKWUwxGAw0GAxShwEASCCLiooNvwBQpHI2/Ha73dQhAAASyaKiAgAUqZw5quFw\nuLDyDwAwX7KoqJijAoAilVNRAQDmVxYVFQCgSFRUAID8ZZGoer2eer1e6jAAAAlksY9qNBqlDgEA\nkEgWiarf76cOAQCQCIspAACplLOYoqoqVVWVOgwAQAJZVFRs+AWAIpVzUdpWq5U6BABAIllUVACA\nIpUzRwUAmF9ZJKpOp7MwTwUAmC9ZJCoAwPzKYo7qyHdRtdvtZMEAABo31RxVFokKAFCkchZTDAYD\nDQaD1GEAABLIoqJiwy8AFKmcDb/dbjd1CACARLKoqAAARSpnjmo4HC6s/AMAzJcsKirmqACgSOVU\nVACA+ZVFRQUAKBIVFQAgf1kkql6vp16vlzoMAEACWeyjGo1GqUMAACSSRaLq9/upQwAAJMJiCgBA\nKuUspqiqSlVVpQ4DAJBAFhUVG34BoEjlXJS21WqlDgEAkEgWFRUAoEibN0dle7vth23vn3huq+2/\ntb1v4rmrbX/D9nlNHBcAUL4mF1PcGRGnTzx+g6Q7JhtExGskXbvajjudzsI8FbDZDh2Sdu+Wtm0b\n/zx0KHVEwHzZkDkq2ydLeoWkd0h600YcA2iSpzoBIR08KO3cufzrs3smHcjXRi2m+C1JvyzpuCY6\ne9/73tdEN5gj0yaeWT0uCQ94QuP7qGx3Jf19RNzcVJ/tdlvtdrup7jAHIpq77dolban/UrZsGT9u\nsv+lbgCesBEbfs+SdK7tuyR9WNJLbf/BejocDAYaDAZNxAasWlVJp54qbd06/snec2BzNbI83fZ2\nSfsi4rRFz3ckXRIR3YnnPlS3/dgK3bLhFwDKVs6G3263u3IjAECRNrSiWqbth7TKigoAUKRNvSjt\nY5KOn9zwuxTbV0v6MUmHV9P5cDjUcDhcR3gAgFxlcQkl5qgAoEjlfM0HAGB+ZVFRAQCKREUFAMhf\nFomq1+up1+ulDgMAkEAW+6hGo1HqEAAAiWSRqPr9fuoQAACJsJgCAJBKOYspqqpSxZVAAWAuZVFR\nseEXAIpUzkVpW61W6hAAAIlkUVEBAIpUzhwVAGB+ZZGoOp3OwjwVAGC+ZJGoAADzK4s5qiPfRdVu\nt5MFAwBo3FRzVFkkKgBAkcpZTDEYDDQYDFKHAQBIIIuKig2/AFCkcjb8drvd1CEAABLJoqICABSp\nnDmq4XC4sPIPADBfsqiomKMCgCKVU1EBAObXzFZUtj8p6cSJp06UdF+icDZDyeMreWxS2eMreWxS\n2ePLYWz3RcQ5KzWa2US1mO3PR8QZqePYKCWPr+SxSWWPr+SxSWWPr6SxceoPADDTSFQAgJmWU6K6\nKnUAG6zk8ZU8Nqns8ZU8Nqns8RUztmzmqAAA8ymnigoAMIdmLlHZPsf20PaXbf+HJV7/PtvX1K/f\naHv75ke5dlOM70dtf8H2o7bPSxHjWk0xtjfZPmj7Vtt/YfuHUsS5FlOM7d/Yvs32ftuftb0rRZxr\ntdL4JtqdZztsZ7OabIr37nzbX6vfu/22fylFnGs1zXtn+2fqv70Dtv9ws2Nct4iYmZukrZLulLRD\n0pMk3SJp16I2/1bS79T3Xy3pmtRxNzy+7ZKeJ+l/SjovdcwNj+1sSU+p778ul/duyrH9k4n750r6\nZOq4mxxf3e44SX8l6QZJZ6SOu8H37nxJV6aOdQPH91xJfyvp++vHP5g67tXeZq2ierGkL0fEoYj4\nrqQPS3rlojavlPR79f2PSfpx21NdhmMGrDi+iLgrIm6V9HiKANdhmrFdFxEP1Q9vkHTyJse4VtOM\n7YGJh09VXhdVnubvTpLeLuldkg5vZnDrNO3YcjXN+C6Q9N6I+KYkRcTfb3KM6zZrieqZku6eeHxP\n/dySbSLiUUn3S3r6pkS3ftOML1erHdsvSvqzDY2oOVONzfaFtu/U+MP8ok2KrQkrjs/2CySdEhH7\nNjOwBkz7/+VP16ekP2b7lM0JrRHTjK8lqWX7r23fYHvFK0HMmllLVEtVRov/ZTpNm1mVc+wrmXps\ntn9e0hmSLtvQiJoz1dgi4r0RsVPSv5d06YZH1Zyjjs/2FknvkdTftIiaM817V0naHhHPk/TneuKM\nTQ6mGd82jU//dST9nKT32z5hg+Nq1KwlqnskTf5r5mRJ9y7XxvY2ScdL+samRLd+04wvV1ONzfbL\nJL1F0rkR8Z1Nim29Vvu+fVjST25oRM1aaXzHSTpN0vW275J0pqRrM1lQseJ7FxFfn/h/8XclvXCT\nYmvCtJ+Zn4iIRyLi/0gaapy4sjFriepzkp5r+9m2n6TxYolrF7W5VtJr6/vnSfrLqGcIMzDN+HK1\n4tjq00fv0zhJ5XSefJqxTf7hv0LSlzYxvvU66vgi4v6IODEitkfEdo3nF8+NiM+nCXdVpnnvTpp4\neK6kOzYxvvWa5jPlf2m8kEm2T9T4VOChTY1yvVKv5lhiFcvLJY00Xsnylvq5X9P4D0OSjpX0UUlf\nlnSTpB2pY254fC/S+F9A/yDp65IOpI65wbH9uaSvStpf365NHXODY/uvkg7U47pO0u7UMTc5vkVt\nr1cmq/6mfO/eWb93t9Tv3ampY254fJb0bkkHJd0m6dWpY17tjStTAABm2qyd+gMA4B8hUQEAZhqJ\nCgAw00hUAICZRqICAMw0EhUAYKaRqAAAM41EBQCYaf8fdIp/VRCA0YgAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "name": "python",
      "version": "3.6.1",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "gist": {
      "id": "",
      "data": {
        "description": "Mt2.ipynb",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}