interval estimates from map

iefm.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!date"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Tue Dec 10 16:35:57 PST 2013\r\n"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# An exploration of a bootstrap approach to getting Bayesian uncertainty estimates with MAP point estimates\n",
      "\n",
      "This notebook shows what goes wrong with a simple bootstrap approach to finding interval estimates in a minimal model inspired by the `dismod_pde` model.  It also shows a potential solution, based on incorporating additional uncertainty with a parametric bootstrap on priors."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pymc as pm, pandas as pd"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 57
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The model in mathematics is this:\n",
      "\n",
      "$ y_i \\sim \\text{Normal}(\\beta_0, \\sigma^2) $\n",
      "\n",
      "$ \\beta_0 \\sim \\text{Uninformative} $\n",
      "\n",
      "$ \\beta_1 \\sim \\text{Normal}(\\beta_0, 1^2) $\n",
      "\n",
      "$ \\sigma \\sim \\text{Uniform}(0, 10) $"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# the data is 100 draws from a standard normal\n",
      "\n",
      "y = randn(100)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 58
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# here is the model above, implemented in PyMC3\n",
      "\n",
      "with pm.Model() as m:\n",
      "    beta_0 = pm.Flat('beta_0')\n",
      "    beta_1 = pm.Normal('beta_1', beta_0, 1.)\n",
      "    sigma = pm.Uniform('sigma', 0., 10.)\n",
      "    likelihood = pm.Normal('likelihood', beta_0, sigma**-2, observed=y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 59
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# this finds the maximum a posteriori value, which can serve as a point estimate\n",
      "\n",
      "with m: X = pm.find_MAP()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 60
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As expected, the MAP value for `beta_0` is very close to the mean of the data.  And the MAP value for `beta_1` is very close to the MAP value for `beta_0`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'mean(y) =', y.mean()\n",
      "print 'MAP values:', X"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "mean(y) = 0.0103592766465\n",
        "MAP values: {'beta_1': array(0.01035925643482589), 'beta_0': array(0.010359276450654365), 'sigma': array(0.975474761072395)}\n"
       ]
      }
     ],
     "prompt_number": 61
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "MCMC computation provides a way to draw samples from the model posterior distribution:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%time with m: trace = pm.sample(10000, pm.HamiltonianMC(), start=X)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [                  1%                  ] 116 of 10000 complete in 0.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [                  2%                  ] 238 of 10000 complete in 1.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-                 3%                  ] 384 of 10000 complete in 1.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [--                6%                  ] 600 of 10000 complete in 2.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [---               8%                  ] 811 of 10000 complete in 2.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [---              10%                  ] 1013 of 10000 complete in 3.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [----             12%                  ] 1205 of 10000 complete in 3.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----            14%                  ] 1403 of 10000 complete in 4.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [------           16%                  ] 1611 of 10000 complete in 4.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [------           18%                  ] 1816 of 10000 complete in 5.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-------          20%                  ] 2023 of 10000 complete in 5.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [--------         22%                  ] 2221 of 10000 complete in 6.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [---------        24%                  ] 2425 of 10000 complete in 6.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [---------        26%                  ] 2628 of 10000 complete in 7.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [----------       28%                  ] 2831 of 10000 complete in 7.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------      29%                  ] 2992 of 10000 complete in 8.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------      31%                  ] 3120 of 10000 complete in 8.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [------------     33%                  ] 3328 of 10000 complete in 9.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-------------    35%                  ] 3544 of 10000 complete in 9.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [--------------   37%                  ] 3731 of 10000 complete in 10.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [--------------   38%                  ] 3857 of 10000 complete in 10.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [---------------  40%                  ] 4046 of 10000 complete in 11.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [---------------- 42%                  ] 4242 of 10000 complete in 11.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [---------------- 44%                  ] 4437 of 10000 complete in 12.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------46%                  ] 4628 of 10000 complete in 12.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------48%                  ] 4829 of 10000 complete in 13.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------50%                  ] 5029 of 10000 complete in 13.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------52%                  ] 5223 of 10000 complete in 14.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------54%                  ] 5416 of 10000 complete in 14.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------56%-                 ] 5608 of 10000 complete in 15.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------58%--                ] 5805 of 10000 complete in 15.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------60%--                ] 6001 of 10000 complete in 16.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------61%---               ] 6199 of 10000 complete in 16.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------64%----              ] 6401 of 10000 complete in 17.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------66%-----             ] 6601 of 10000 complete in 17.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------67%-----             ] 6797 of 10000 complete in 18.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------69%------            ] 6982 of 10000 complete in 18.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------71%-------           ] 7143 of 10000 complete in 19.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------72%-------           ] 7297 of 10000 complete in 19.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------74%--------          ] 7487 of 10000 complete in 20.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------76%---------         ] 7683 of 10000 complete in 20.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------78%---------         ] 7856 of 10000 complete in 21.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------79%----------        ] 7986 of 10000 complete in 21.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------81%-----------       ] 8192 of 10000 complete in 22.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------83%-----------       ] 8397 of 10000 complete in 22.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------86%------------      ] 8601 of 10000 complete in 23.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------88%-------------     ] 8802 of 10000 complete in 23.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------90%--------------    ] 9003 of 10000 complete in 24.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------92%---------------   ] 9212 of 10000 complete in 24.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------93%---------------   ] 9387 of 10000 complete in 25.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------95%----------------  ] 9522 of 10000 complete in 25.6 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------97%----------------  ] 9720 of 10000 complete in 26.1 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------99%----------------- ] 9915 of 10000 complete in 26.6 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------100%-----------------] 10000 of 10000 complete in 26.8 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 34 s, sys: 1.34 s, total: 35.3 s\n",
        "Wall time: 35.7 s\n"
       ]
      }
     ],
     "prompt_number": 62
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It is important to check convergence somehow.  For this I will do so by eye:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(arange(len(trace)), trace['beta_0']);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD9CAYAAAC4EtBTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl0FUX2x78PgsAoKBAYkRcVWcbsEIk/EJEIhpCwKIjE\nBRiVQUFABcVtRhN0gEFHUFHQ0fGARkcxisAEwogYwFFZgyHqqCBgEsUECBIUQyD1+6PpvO5+ve8v\nuZ9zcvLe6+qq29VVdWu5dSvAGGMgCIIgiDM081oAgiAIwl+QYiAIgiBEkGIgCIIgRJBiIAiCIESQ\nYiAIgiBEkGIgCIIgRFhWDIWFhUhMTERcXBzmz5+vGO7dd99Fs2bNsHPnTqtJEgRBEA5iSTHU1tZi\nypQpKCwsRElJCfLz81FcXBwWrqamBs8++yz69u1rJTmCIAjCBSwphi1btiA+Ph5dunRBVFQUsrOz\nUVBQEBbu0UcfxUMPPYSWLVuC9tMRBEH4G0uKoby8HDExMQ3fg8EgysvLRWF27tyJiooKZGVlAQAC\ngYCVJAmCIAiHibJys1YjX19fj5kzZ2LZsmUNvymNGEhhEARBmMPumRhLI4ZgMIiysrKG72VlZaIR\nRE1NDb744gukpaWha9eu+OyzzzBy5EjFBWjGGP0xhpycHM9l8Msf5QXlBeWF+p8TWFIMqampKC0t\nRUVFBerq6rB8+XJkZmY2XD/33HNRVVWFffv2Yd++fejbty9Wr16NlJQUy4ITBEEQzmBJMbRq1QpL\nlixBRkYGkpOTMXr0aKSkpCAnJwerV6+2S0aCIAjCRQLMqbGIQQKBgGPDokijqKgIaWlpXovhCygv\nQlBehKC8COFE20mKgSAIIoJxou0klxgEQRCECFIMBEEQhAhSDARBEIQIUgwEQRCECFIMBEEQhAhS\nDARBEIQIUgwEQRCECFIMBEEQhAhSDARBEIQIUgyNmEAA+O47r6UgrPLNN15LQDQ1SDE0cqqqvJaA\nsMof/gB8/rnXUhBNCVIMBBEB/Pab1xIQTQlSDARBEIQIUgwEQRCECFIMRKPl228B8uROEMYhxUA0\nWnr2BLZt81oKgog8SDEQjZoTJ7yWgCDso6wM+OEH59MhxUAQRKOlshKIjvZaCvvo0QNISXE+HcuK\nobCwEImJiYiLi8P8+fPDri9evBjJyclISkpCnz59sGPHDqtJEgRB6GLvXuDwYa+lsI/aWuDnn51P\nx5JiqK2txZQpU1BYWIiSkhLk5+ejuLhYFGbChAn4/PPPUVJSgpycHMycOdOSwARBEISzWFIMW7Zs\nQXx8PLp06YKoqChkZ2ejoKBAFOacc85p+Hz8+HF07tzZSpIEYYhAwGsJ7KGxPAdhHTcs7aKs3Fxe\nXo6YmJiG78FgEEVFRWHhFi9ejAULFuCXX37BJ598ohhfbm5uw+e0tDSkpaVZEc9zXnqJs4y5+mqv\nJSEIorFQX1+E3NwiR9OwpBgCOrsxd911F+666y7861//wu23346PPvpINpxQMTQGJk/mFopoWYUg\njPHZZ0BSEvC733ktif9o1iwNublpDd9nz55tfxpWbg4GgygrK2v4XlZWJhpBSMnOzsa2JmZYThus\nCMI4/foBr7/utRRNF0uKITU1FaWlpaioqEBdXR2WL1+OzMxMUZj9+/c3fC4oKEBsbKyVJCMOrxWD\n1+l7Dc3NRy7Hj3stQdPF0lRSq1atsGTJEmRkZKC+vh7jx49HSkoKcnJy0KdPH4wYMQJPP/00Nm7c\niPr6enTo0AGvvfaaXbJHBJHUMDdvDhw9CrRp47UkBGEP1DEwhyXFAACZmZlhowThnNeiRYusJkFY\nwEjFqK8HjhxpXIohkhQzEZmcPs11qhoTtPPZYahh8pZ33vFaAsIrjhzhFrHleP99e9yllJQAUZa7\n18aorXU+DVIMDuO1YvA6fa/Zu9drCQizWJ0GWrAAmDFD/tqoUUB+vrX4AeDHH63H4UciWjFs3Ag8\n+KDXUqgTKQ3zwYNeS2AfVVXAl196LQXhNvfdJ/ama2bx+sABOkYViHDF8NxzwJNPWosjM5M8cALA\nNdd4LUGI3bu59Q6z3HQTEB9vnzxEZLBgAfDKK9biyMwEevWyR55IJqIVgx0UFjrrxtbrEYPe4fgX\nXzgrhxGSkoAVK8zf35gU/bFjXksQWQjrm5m699VXxsJr1S+r9aq62hvLqohWDNXV6tcrK92RQw2v\nFYPX6Zvlt9+8lsB7Nm4Ezj3Xaym8ozGYmiYkGFc2Qo4etU8WI0S0Yti6Vf36738P/PqrO7Io4ceG\nec8e5Wt+lLep4oeOTaThZPlt1y58JKtHeZ086Yw8ThLRikEPp097LYH/6NGDOwlKD19/7awsTkDK\njXCCo0cBFR+gjYqIVgyRMNT0upFSyiOlXoww/J49wKWX2i+TXgIB7elCPXG4TW0tcNFF7qdLWF9j\nIDgiWjHY5UslEhSMWaxUDuFGmi++8Mak1epUoBfv9uhR4Pvv3U+XcB895SsS2xdfKQYj1iTl5erX\n339fO466Ou6/kz0Lrbh//pnbPWkn+/fbv3ibkABMnWpvnIR+IrFx8QIjdVnmJGJHiMSRi68UgxHf\n67/8on59/HjtOP77X/3pOcXMmUBysr1xdu0KOHG0hRdn57pVqTZtAtascScts5SVqRsOEGK0yo4Z\nU9KmoqBd9vLhL06dcj4NrcLplFkmPzdvpSAbvXfPHs5FwIAB5tM0m7YQYZ7rjWfUKM63jp97dzfe\nyO252bfPa0n8i9vvz2lF4ZUi8tWIwW3cKEReNzTC9F96CWjRQn94KVqFNDsbuOoq/bI1VpyqzF99\nxU0TarFuHXDDDc7IEEn4pXfvFzmM4DvFcOSI8Xuqq4EXXgB++slcmk6+OC3FsH27vempra1s3Roa\nJZlRWF4oOatpeq2YvaCw0B4HcUbgjQROnADS0+2JU65efvihen112yopEht9PfhOMWitHciRnw9M\nmwacfz7nZyeS+OYbe+MbNYr7L1cpzOStEeTS9MKlg9DVst6K25gUiNao0G4qKoCzz+Y+//ADsH69\nc2nt3Olc3HpQKk8bNwJvv209fsY4R37C717gO8WgF6UXxJuw6slQPoyXVklGGTMGePxx/eGF+WR1\nPcNo76iuzrxLB7vyTSizV+4F7EJvnpx1lnMynDoVsubjIX9OwO23c2tAVtm4Ebj4YuvxWCViFYMQ\nP/f2jMhWXw/84x/qYd5919gh6UbzRtiQWh0mW/GQaifPPssp1Hbt3ElPT77t2uVc+k6OGLKygJQU\n6/E8/TTXCBJiamq4/++8o113b72VMyN3gkahGOSQ9mq8wkjDfPgwcOedxuJ/9135/R9y6fKFTi9+\nULjr1xszY5bj3nu5fFLD7LMOHWrOf3/v3s4ddq9HMWzapM8q7+RJzo05z5YtQGmpedl47r+fM6kO\nBJRdZZvpmETKzufKSqBTJ+XrY8dql49ly5zzimxZMRQWFiIxMRFxcXGYL7Nj5KmnnkJ8fDwSEhJw\n1VVXYZ9LtnZGHFf5ZQGpmYm3MWYM8N57+sJu2GA8fq/ZssXfbrTXrQP+/W9z90obLsbseVY9inTg\nQGDVKs6XmNpRkVVVwFtvhf+enQ3062deRiGTJtkTjxvY1Vbs3cvlLc/y5ZxnAWGZOH06Qs1Va2tr\nMWXKFBQWFqKkpAT5+fkoLi4Whenbty927tyJ0tJS3HzzzZg5c6YlgeXwc8/AiGxOFAKjcZrZAyB3\nr5n75e7Vyr8dO4Crr9aOh2fVKuAvfzEvlxxmR6fSZ3v+efVGXW9ZUuuJCjl9GrjrLvUptkOH5H9f\nv175PGWn2LJFuxctzKMXXrBfBukmT6t1ljcIyc7mpteELmDc2GelhCXFsGXLFsTHx6NLly6IiopC\ndnY2CgoKRGEGDBiAli1bAgD69++PiooKK0k2oORDx8iL8ts+Br0jBukzyqWhtQhv1rQX4E65+vln\n8/frQe80S2EhUFSkP95584A5c8J/t1IW6uq4RtbqKXh27Wo28iwlJeqjFCtnCRhh6VLguuvUw/Tt\nC/ztb66Io4jdbca8eeLpIGG/2UvP0JZ2PpeXlyMmJqbhezAYRJFKLX3ppZdw7bXXqsSYi4ULgbZt\ngbS0NKSlpSmG7N3buLxe4IRi0BOn3PBfyHvvAf3760uPh1dIhYWcO+7LLxdf5xdUT5wAWrc2FrcU\nvftZ/DANWFfHWXx9+KGxyiw3lWSVvXuB//3PXPpyWMnfXbuAmBigQwftsKtWAStXaofT6kVH2s7n\nH34ITaVt2MB5DuBRftYi5OYWWUtYA0uKIWAgV9544w3s3LkTG1VNEXIxYwZXmIxgtjDI3de6Nfey\n3LJgESJdfnniCc7c8+67xb/b0bO0ajGklud33KFuObVnD9C9u7X0ebScKbqB0JPq0KHAm296J0tq\nqjFX5fx7rKvj9gDZYXEEcKbBvXtzi6h22Pfz6Knrhw9zi+tOIG3yPvrIepy8ua90j4ayYkhDbm4a\nZs/mv89WCmgaS1NJwWAQZYITX8rKykQjCJ7169djzpw5WLVqFVpomEx4vV7w22/OngGtxqBB3H8+\nDx57DHj0Ue37zOSZk4qBny1U6jf06KE9Inj4YX1yLFmift2NEYVw34CVg1y03qP0+uDB4c9vdr1j\n6VLgssvM3SsnNz89pWYEIpxV1luG9eTRM88Ao0fri4+/x+weF35PkRVXMkr3ejmVZEkxpKamorS0\nFBUVFairq8Py5cuRmZkpClNcXIzJkydj9erViI6OtiTsd9/Z613SjUbDSKPNN5ZuKEelQie3eW74\ncGNx65Ff78Ka8B2tXWttbUQNO9dMlMrVvHnWDx4SsmEDZ+9uFsZC039enLH97bfupynHypXOzhCY\nrc+nT3unHCwphlatWmHJkiXIyMhAcnIyRo8ejZSUFOTk5ODfZ2z4HnjgAfzyyy8YM2YMevfujes0\nVpjUGuvkZKBnT3WZnFx83rPHuqWOHoS9eaeUhNKIIS8v9Jl/Vok9AQBjO8uNXhOyd2/oc1YWZ7mh\nxl136YvXThgLPY/QKCIxUTw9+Mgj6q69zeSX1ZGfmqkq4I81HCOYqS/CeX2zGMmnkyf1W8YpTUuu\nWqU/PTNYdrudmZkZNkqYHZr8wgcffGA1iQZ++827qaZffzXn8M7MHLgbz+jkVJKeiqb3GaV+d7Ts\n/LWmlpxAafRTWgps28adjyGEnx7UWnweOzb0Wel9+WF3uR3KQ6488PEeOqT/kKg33uAUsp+Q5s/e\nveI9RWp1Qcn6UtWGxwYiauezngJopFH98ENxvGprC0OHineAOokblV1piGrHXO/XX3P/1d7XtGnA\n5s360rITJ/JWrfcml89//avxNAIBedm9XpMTEgiERiA7doR+E/LVV8bXQbZt4zaA6cVJR5p8fs+d\na8/GvOpqZZNgL5W+rxXDCy/os2c3Wzn43iV//8CBymH5xk6K0iKqlX18RqeSzDy/kqdVvvc7apT6\nkaNWG6R331V2haCG1d7p1q3q16urucVLI0gbf6GMZs1X5RoFubpQX8/ZwfPnNNg9laonjNDVO7/Y\nPGuWfNi4OO5cED188w33PFlZ+sJr0akTt0nODNOnc///+U/u/wsvmCu/gPgdqRlhXHqpd2eH+1ox\nTJsG2DgTpYmaEpKrcL/8omyjvXCheTnc6AX+61/q199/H1i9WvybMA+MyGjH87jVM54+HZgxw774\npNNMwp6v9JmE+w+ki9TCxfHCwtDn+nrOkZrKlh/D/PqrMQd3jzxiLH65Tonc+7XbxXZVVXgHT68i\nff558+laMThQ2nnuNL5WDIDYNznPokXK4c30mA4e1A4rt/nMqS3rbg8h9fZqvXRQxo/ujI4YjIY3\nc1BUba1yfqxbJ74mnHb605/EYdVs4jMyQp+FS3p83GasVwRLgSIZFy/WVjRK+TR5Mvff7MhOOEqV\nq/t28OKL8r875b5GOhNhxIOxV24xfK8YhL03PhOlG76sNlJ6dgDLvUBpuitX6i8wf/mL8ghFWMn5\nNAIB5akdxqyZW5pRRHryXM42/M9/Vr9Hye6dX3y04kZcD2b3ASjNE7/1lvKOdr3OD+Xg5dR6d127\nii3NTp8ONerCEYqwzAk/8/mnp/MEhKZq+PtWrNB3H4/QF5Gd+4mEvXal0bySwvASUgwOMXq0PSv4\ndvvYnzNHeb5TaKIpRKkH9fLLwHnncSfY2QHf6Hz5pfh3o1NJy5aF/yadoujRQ1z4tXYNL1okv5fB\njIfTV14BPv5Y/JvW6WOffy4/76tkPeIUEydy/3nFoOR0cP9+8Ujkuefkpz+Fa2JyjhSVyuSnn+oW\n2VO++047jBVTZ6dGGxG5j8EtCgrMD03Xr1e2GvHa86kSRj2c/ve/3H+7Nn899BD3X26Olze/NbvG\nINwlHAhwe0OEtvR6ekhyJ4aNGCEfVq3HOmkScN992ukJ6dULGDYs/Het/LjhBmPpaMXHz5ULR2U1\nNdr7Eux2IaLH/cexY0C3bqHvXuyN8JP1lhAjecEvfLuB7xSD3AvkeytKmSjwymE5LSXkKhR/f02N\nuDev92xlu81v7ULJnfKGDaFRjhm5Tp2Sn0f/5z+5aZXYWG9dDetFbvpG6CFTrnedn2+vDLx1lTCt\n88+3/4AqOxrxigpxj12uA2O2nG/aZP3AI7VnLCsDOne2FocdtG/vnP8nOXynGOTQ8tRp1hWvUmEs\nLTXWQLVtKz6nVct9sBZKctmtJKqrzfuJMWPqWFkpH+6ee4Drr+fmvN0+eU/pOQIBbsFXeD6C2gl4\nvDkjYG6fglUCAXtc0QPKRgZGzthSS3PBAmPyqDFwoPZueCF//GP4b2pl+csv9a2vmJ3y4c2MleBl\nq65Wrj9OEBGK4dln9Ye1Q3PfeKPYRtloQ212BGMULZt8gDtlS2kPRvv2nBmeVZcAL78sP4XB5w//\nv1Ura+m4zX/+E9ppXVLCdQCEKB0XakcZNDr1pBejezR49u8HLrlE/7O5OV2kx3hCrfFXa9SVnkP6\nuxELKuG9WopBiFNHwcrhO8Ug9yKcOtdUrfAKXS+sW2csXr0NoPDMBKE/ImmP7bzzuM/SgvF//6cv\nnUsvVb4mtfDSi1DGO+7Qd5qXNL/1WHr5BenJXR9/zB2r6hRGDh8yglYjumKF/BoX72TPzPvRc8yu\nm+9dmpZwtCfFa19RwvSlhhJOYtlXkpvofUl6NauwgPBb+OUweg7vmQPrZKmvB5o35z6//HLod6EH\nU6li4E1Rb7nFmBx6MVIp+ekKI+aqr7wCXHRRuCIz2xg89BDnD8fqiWlWZFCbfnPawZkcdjVgkydz\nU3tSrDTcjz2mHcbNjax6mTEjVFftxGtlo4eIUgx6XQMLF3+bNdM31OzTB9DrFXzLFu70MqXKonYS\nm9EKZsfh8FoIXRpoMWEC91/LARzAuZkGuMYmJoYz9bSD997j/tSUuV70OEYMBIDbbrOell7sdM0N\ncIuW1dXmXEvz78yKYtAzRy+ciuR9mKmxYQN3HgXATQUZMRfW2zA/84zYik5PHIsXc2eQx8bql0cN\nvac62o3vppIALtOlNvRmsbOR5unblzuQxanTw/TIZWevQ2qXr3dheeHCcM+hauhR0A8+qD8+NxE6\nZgsEnJ36aN/e3vi+/dZYnMJn82Ih/T//0Q4j7GTMmcOZEduJ1r4HpbI8dSrw5JPhvwvrq56pNa/x\npWIAlBdMzSK3+UxauYV+SfRsVTdjXeFHE1WzO583bNC/eFZWpm9E5sWBMXIb8aSY9RPlJWYsZZTK\n5x13qF+3E7mGVYp0c6jWoT9G84Lfd6H0vMLT5/TC1zMjzve8Kmu+VQxm7dnr68Xn1vLxVFVZl0nK\na6/J/y59mVqbjrxm6VLrcZgpwFYbGbumXG691Z54/IZZKxY1x216Xa84vR/F6DnS/fqZS8dMGVUz\ngQaM5Y2bC85CfKsYzNqz19Xp97Gi1pjNmMHNdSrJwc+f6+HOO/WH1ZLLLHb39FautB6HVYXE7/h2\nA6mljt3vyOr7sfP9GjFpVcoHoR8mI/htNObECCkSFp99pxj4guFUj+OKK0Kuj7UK4XvvKS8+GTFh\ndcpLpJe8+KL/KrERnDo72iiffGJPPHJuQvQ6vhPidb64MZVotdwaHa1EIr61SjJbQLS08aefcjup\nO3fWDsvbrhcX60tbq8BNnqzvrNdIbHD1WCn5CbscDlqlf3978krODbba/hW70NP7NdJDHjnSvCx6\nMZLfclaBN95oLf4mMWIoLCxEYmIi4uLiMH/+/LDrmzZtQkpKClq0aIF3lbaKymDWjbRapvMHnGzY\nAFx1lfYL5HsGTzxhThapTC+9pM/iwutG1W5zSbeQs7+3C6FZq5uuCZoaXnkTdQOhC32/Y0kx1NbW\nYsqUKSgsLERJSQny8/NRLOleX3TRRVi2bBluvvlmQ3ErTeFYQXjAiRMovXCh4zivG309GNmmzyN9\ndi+e87nn3EnH7mebO9fe+PwGn19eLaT6herqkOdiv2NJMWzZsgXx8fHo0qULoqKikJ2djQKhbwdw\niiExMRHNdO7U4BuYe++1Ipk+7K7g/EK1EWdjctTWAmPHWpfHLHqn8YQmg4MGia/ZtQ/Fj9jd49M6\nvMgucnPdSUcKn19Wjru1Ezc6LYFA+OFNH30EzJ8fGSMGS2sM5eXliImJafgeDAZRZMnJS66g8KSd\n+TOGkUy3+whNfhhsdarh6FHgnXesy+M0amaNehz8RTKRMPKTIjzK0034OqnmnbYxwRubOKfwi878\nOYclxRCwXfXlYsYM//Qs7IRXFnoaFKEPJcJ/OLEnpingF39ITit1/pRC6SFRfBugtRlPmzSIO832\na3xLU0nBYBBlAh/TZWVlohGEFPsViVwajidhit//Xn9YfpGc8C+ROGLwCj2ed93Eq3fHn1tu5hha\nt7GkGFJTU1FaWoqKigrU1dVh+fLlyFRY4WWMgbnwRryssGYPDCKISCUSFaTVNcCmgCXF0KpVKyxZ\nsgQZGRlITk7G6NGjkZKSgpycHKxevRoAsG3bNsTExCA/Px933nknEhMTbRGc8DeR2GAQxrHDwy3h\nPyxvcMvMzAwbJcwWrHKlpqaKppu0sHq0Y//+8r/7ZYqpqTSYkeBB0gpN5T1qYcTdNRE5+M4lhlNn\n/vqlIsvtUG2M3HWX1xI4i93efyMVv9Qrwl58pxgMbI6OSCJlgwuhzsMPey2B8+ixDhSeU0E0Hnyn\nGJxyouWXqSSCiBQa+14UQhnfKQan/bgTBKEPmiZquvhOMRg554AgCOegTlrTxXeKwSn27PFaAoKI\nLBr7eh+hTIC5setMB9yuaF+IQhAEEUEEbN883GRGDARBEIQ+SDEQBEEQIkgxEARBECJIMRAEQRAi\nSDEQBEEQIkgxEARBECJIMRAEQRAiSDEQBEEQIkgxEARBECJIMRAEQRAiSDEQBEEQIkgxEARBECIs\nK4bCwkIkJiYiLi4O8+fPD7teW1uL7OxsJCYmon///jhw4IDVJAmCIAgHsaQYamtrMWXKFBQWFqKk\npAT5+fkoLi4WhXn++efRuXNn7N69G7NmzcLdd99tSWCCIAjCWSwphi1btiA+Ph5dunRBVFQUsrOz\nUVBQIAqzZs0ajB8/HgAwcuRIfPLJJ7a7iCUIgiDsI8rKzeXl5YiJiWn4HgwGUVRUpBimWbNm6NCh\nAyorK/H73/9eJsZcwee0M38EQRBEiKIzf85hSTFwh+vYSa7N8REEQTQ20iDuNM+2PQVLU0nBYBBl\nZWUN38vKykQjCD7M999/DwCor6/H4cOH0bFjRyvJEgRBEA5iSTGkpqaitLQUFRUVqKurw/Lly5GZ\nmSkKk5WVhby8PADAypUr0a9fPzRrRlayBEEQfsXSVFKrVq2wZMkSZGRkoL6+HuPHj0dKSgpycnLQ\np08fjBgxAtOmTcP48eORmJiINm3a4M0337RLdoIgCMIBAswnJkLceoUvRCEIgoggArZbetKcDkEQ\nBCGCFANBEAQhghQDQRAEIYIUA0EQBCGCFMMZ3n/fawkIgiD8ASkGgiAIQgQpBoIgCEIEKQaCIAhC\nBCkGgiAIQgQphjPY4Sh22jTtMAkJ1tNpTMTGei1B4+Xii72WgIhUmoxiyMhwPg3ZIyYktG/vvByR\nxH//67UEjZehQ72WgIhUmoxiGDUq9Dkuzpk0yJu4cX73O68laLzYflwK0WRoMorBDS6/XDvMbPvP\n1IhoqPFyDvJuT5iFio6L/Por0L2711L4C1IMztG5s9cSiNm0yWsJCL00ScXwxhvhv7nVQFFDGHm8\n+qrXEphj2DCvJRAzYIDXEhB6aTKKgW+Qr7sO6NUr/Lo/TqXg0LOI3ViIBEV5+rTXEhBSrr2W+19d\nDdx6q/PpnXWW82n4iSajGHhWrAh9jo+3N26tRi4Q0NcQHjgA/PWv9sjkd0gxuMell3otgX08+ij3\n/7zzgChL51Dqo6bG+TS0SE52L60moxjkGqDU1NBnpRHD3/6mPw27FvtatgT+/Gd74tLDoEH2x3n7\n7frCRYJiOHXKawkinwsvtDc+N5SBkObN3U2P58YbQ587dXIv3SatGPRMHz3wgP40jCqGxmaqOWaM\n8XsiQTF4MWL4v/+zHoefpkcnTXIubj89p9306eNNuqYVw5EjR5Ceno6kpCRkZGTg6NGjsuGGDh2K\ndu3aYcSIEaaFdAOlRtpIw9W6tXYYYXxr1+qPOxJo0cL4PYEAsGyZ/bJocdtt+sN6oRjatbMeR4cO\n4u9299qNYLfprNsdCjfSu/nm8N9uuMH5dOUw/bpycnIwbNgwlJSUIDMzEzk5ObLhHnjgAbz++uum\nBZRy9tm2RSUiJcV6HN26qZvkSdcYzjkn9JlfTGtMXHaZvnDChb22bZ2RxQotWzoT7/Dhytes9oJv\nvlmsCAoLgZkzrcVpBa+mYuzGSeUqZy3p1V4U08muWbMG48ePBwCMGzcOBQUFsuEGDRqEc4QtoEXM\nLgKZnUoyipxJXlKSdnrXXGO/LHoJBpWvWVnwGjfO/L1+YtIkYNs2++N95hnla/X11uKWGi+0bett\n4xzpikFELjf/AAAd80lEQVRrxODUOpSwzXBzlGRaMVRVVaHDmbFqdHQ0KisrbRAnV/BXJBvCzswR\nZrpao/3889bS+fzz0Gdefr+MEG67TX3+l5/rvuoq7bik7+bee/UpX+F9VpT1vHn6wxpJp0WLkHK3\nEzUZrCqGrl2t3Q/YazbtpmK4+mrn4lZqf5o3d9OktQjittJ+VBVDeno6EhMTw/5WrVrliDDih01z\nKA3jTJ1qLLySlY9ec1W36dFD+Zob8vbvb0886en2xOMWaopBaoHSrZv19Iy+y4MH5ff8mCEz0554\neNSexYkyqydO/n0OH859lnamLrrILmnS4Kli+OCDD7B79+6wv5EjR6Jjx444dOgQAG700EnFlirg\ng9ZQa2HYbcsGxsSFzQdZJIuRfDGbh8FgqLfv1nswavXjxPtRe9aXXxZ/V1jCM4SZZ3jtNevp5uXZ\nv4ciGOT2MLiNdCQnXD6VXsvPD33++mtgzx7n5LIb01NJWVlZyMvLAwDk5eUhKytLMSwzUNuLi81K\npM7IkeGLb3JiSS05zKD0uIzZY/1kBiWLIb+Y+vHPb0UeI3k4eTJgy+ynBdSmi5wysogE9IxS2rXj\ndj0D4WXGybqkthjMyyFXx3v2tL73IiLWGGbPno2CggIkJSVh7dq1ePzxxwEAO3bswCTBpPWAAQMw\nduxYfPjhh4iJicEHH3ygGq8Tc7kAZwH09NPK1/mXeugQMGuWtbTUnJedcw7nnjsx0Z2zGfhprSlT\nrMXjRKEU+vLh899NfzpOWRvphX/m6Gjuv5qNhl8U+IwZzsXNlzGr1lNe5RW/G1vLcum668zF76Zi\nMK3D2rdvL9vIX3bZZXhZMA7evHmzoXj9UgEA87IoVXD+xf7vf1wY4WKVUy990SJrrj+cfB9yzzxv\nHmdaaVd8arRty7lIEZ7VYVfceuDztmNHrkPidNm34xm6dDF+j9Hn8rO7cDXZcnPFbvWV8vvcc63L\nERcHfPWVc2XGx6/APHo3Wum1SjKKVlzt2+uzYJDuAzBqaWKnvyW5Z7r+eu0wUpSWovhKpLfxGjoU\nKCvTF1YNu3af692zIQffv2LMfv9dQvTkrZL1EO8BQPiO27RRjkfNoEELo2sHU6YADz0U+i73nLt3\nm5dHSCAAnDjBfVZq4Pk8UqoPRhS0UhxffOHsvhRfKobcXGv3b9gg/q5U0LxYcFaCL2xqGK1sjzxi\nLLxR7JwH5/Pm/PPDr331VfhvLVqo77+wG63KvGWL8Tj5Z+7SBbjlFmD8eOBPfzIej1msbtbSsgxS\n26chZPJk8Xejc/F9+mibKvNnrcuVLyMEAkCrVtznu+/21mjEyfbLl4ph4kT19QC/sGiRfXHJvWRp\noTNTCLUKj9BJF49cT0Sa9pklJdUwWsiFF9rOP/ss9//CC/VVAjd2TW/YID9VaObdCBef8/KAF1+0\nvn9BDV7Gt9+2HocVHnxQ/J0fpRgdNSph9n6ja3533mkuHSPw5f7TTyNk8dlJgkH5xok32XvzTW6V\nn0fFIMpRjG4AUnuxcg2f9De1+U25DXN6ClJGRng6QqWs1CDfdJO9BVUuLjVLJbnwdtj6a3HllSF5\nrDYMcs+lVzFMmGA8PT7PrrxSOX2ld2p0E6JSPDk5Ic+7vIKVulFzsgFUk90uL8NO9OTtcKpoBF8q\nBi1uukncSJqtoPyQEDC33qA0X26mYOhpEPjhsBzvv288TatIn9NKhfCT0YEcdjVW338f+mxEMUjT\n58ueXo+2rVuH98j59NXWCsyilF+5uaFOHZ9+Roa+e62mbcd9dozijaDmEsNRs1znojaHWq9YaUOY\nmUZl1y6xBYGZOAYOVL8eE2M8TiHCZ+zaVf5sCOEi9rp14dcvvthYOm5iR7p33SX+ztu22yGDVqNr\npswIy4QVxWDk1LLFi9X3Bhw7pj8uIxjNH/4Z9Xgp1hOPHGoy6TFEcGJTqt7d5WPGAGlpwEcf2ZOu\nGr5SDG3bmtvqbqaCJiebXzzVe4hOx46hz9u3y9vNS3tPSrRsKa80R41S99JpdYFYTS6rvfy4uNBn\ntfcu11Dwv0krs527YYWLoGY6Im++qX7dylQS38vn90CoIZ1i4++VS/+KK8TflToWK1cC776rnKaV\nRtOM+3Y7UDp75f77Q5+1nqtHD+0pKWm+q71DYXq33cYphbQ09fjtwFeKQcrEifrCSTNabcpFjvnz\ntSu7sGeemBj6/Omn2vIA2uaMehuEL78EPvkk9P2tt4DVq/XdawW7RxVt2wJz5oS+q62xSBXiggXA\nkiX2yHHZZcac7wUCXMfg4Ye1PV9qWU3JvXO95SAqijPXXbBAX3ggJCOv8KR5PnUq8OGH4t+k00z8\nPSNH6vcIrFS3amvV5XQb4dSykKeeCn3Wms755hvxxk09XHKJ8jW1fSNNaipJyCuviH3EKGWEtDIZ\nnTM1ckqblL59la9Z9TMkN2yNjXX37F4jhS82Vv16dXUovrPOMr+RKS7OnoNsAM4SRWgDr4dHHgHm\nzg3/Xdrw8gwdKh+P3OhmyBCgXz9tGZo35xSPlWkXac88Kor7s2MvhZ5yI3VVbZdVkh0byNxk0SLl\nKVC3jzDl8bViAPQVErXpko4d9Vkb2LnZzcz9fl985VHq5WdmAn/5i/J9w4apT/P4yYmg3P4PPed5\n8Gs8o0frS6dbN+D4cfFvl18uHhEqYcaNtfQZpO5H+OcZNUo8DaoWh5r7GK13etll2vPr0k2UQLi5\nq5R588L3MvFI39mXX3Kjz3PPNVcG7bjnrLO8cQiohq8Uw+LF6teFGSqsGFKrBiGffqo8F+qmRYES\nfAXUCiuU1clexIoVoc+DBoXkuvVWYOzYUKMplVdrBKBlwaQnr5zcFSxEqPzUUHoP0n0CfM9cqgQA\n/WtARkyXlZD2yJXK/xNP6HcwqOb+hZdZac0gOTnkNPPtt7lpUalcy5eH33fHHcryBAKcwhWas6sR\nG8ttsFM4mTgMq5tv7eTqq+105S3GV4rhllvCfxMWEmEhFM7ZBwLi81L19r6F4YwcMqcVv5GprNWr\ngfJy9amkHTsA4QF5bdrIb/G3Y9QhlF3Y8/3jH61tjDKDtOEqLQ0P4+SmNq2dwfPmAZ99xn1Wy3u+\n8bdzp7jaiEGp4ZTmp5HyYsX7rR6T2rFjudHBunVA796h33kFqPdccF4h8JvVxo4VX7daR6xaGmph\n5CTHrCxg/35n5PCVYtBC6CFUWsiF56VGRRnfrHLuufZN58yZI+/GQY527bgFJrW0U1LCrUO0CqjS\nVIAWRixChOdka+Wdnf5jhGEfeMA+f0dS/vxn4MgR5evnncdtPNq7l2vA5J5t40b59Qij76dPH/F3\n6UIpb/F2/vnASy+pxyUn58CB4RvNhFjJY7mR1bJl8tZ9Q4bIj4aE5VKtrPH38msvRkZWcmVQy9uv\n2VkHleNrZHF7qtn3ikGY8XoaLb6C2LWL0Qxnn218gVjuxd9wAzB4sDkZzBwzuGyZuCL06KFeIGfM\nCF03agkmNRgwW/BbttQ/bWCU5s3DF7nl5FSyKmGMO8VLzhzx4Ye10xeeHCg8J4SxcNNnPQ2UdD5f\n+CxFRcq91XvvDcmiNOqQuofRkmfCBHVrHLNIFYH0u55yxne6WrQANm0SX5M+l9nj7I0c2lNQ4L7H\nWd8rBiF6Cr9T8+967MWlWLVKuv9+YP16/XGo2dqrue/gG53+/cVTBkOGmEsbUF73+cc/uP9OHZ7u\nJXbuswHEZ41rlf1AgGvQ1M4zkVow6ZVr4UJzZ4d4YVAhzScjDeoFF4jjUHPTwpOVxZmoyqGmNIxM\nN3vh8sf3isGr7fFSduywZtaqhdFKZOSQmR9/VJedV6Z2zoErze3yIzmpDbsdrojVwrVubd2bqJ2Y\n7Wmq8fXX7uxpEaLHt5Ld6ajFrTViEN4rPSPc6Cg7EOD+lLwe84rYzOjda3yvGHjmzbO/8Wjblhsm\n67n3/PONLzwZaeyNetVs1Uq/Zc/554ePpOSmuqz4zDHaEJw8aT4tKStWaLu97twZOHDAWjrx8eqL\nvnrf9//+x+1itZvWrfU1QlrWQkZwa1TAmD433lJFIH1f69YBO3faJ5cepJsx1fLML2bbvlcMTmZU\n8+bcMFkvf/qTc35K3B52S0ccVVXyIwY7FovlsNPs7+KLOft/pykttWYmyvOHP9g/5WnmfbRs6Vy5\nc6Le3nOPdhjpme3S99W1q9jqScq2baGDk+Sw+7n0rGF5gUf76pyDz1gnCmarVpyfEt7eWq8sdoc1\nG9/8+aHNQSkpwOefh66ZWUPRy+zZoUVtvkEU7mh3A6uO2fTgl95eY0RP3h44EG6+3L+/8hSpXJxS\n6y8n4NPt2RNhpzI2CsVw5MgRZGdn46effkLnzp3x9ttv4zzJFr7i4mJMnjwZJ06cQF1dHR5++GFM\nMOBMns9EJywYzOLEy3OjQAinisyc3WuWxx4LfW7dWn6jlxxqPTsjvPNOuHM4PRht6L2q1EOHeucG\nQqkjFgioL4SbTUf6WYjcGlL79qEDn/7xD+s7jIXPadVd+e7d/u1MWJpKysnJwbBhw1BSUoLMzEzk\nyHQD27Rpg+XLl6OkpAQbNmzArFmzcPjwYcNpGV2Zt1JJJ07UZ05opwx2NypaFhVmFnCVGDJE3cOr\nFL2L3MJDcYQYrUw9eoQsTjIzjd0bCaxdq38Uy/PYY+J9QVaRTuUFApyptdc9YOFU0qRJ5vfLSH9j\nzLoBwVlneedJVgtLI4Y1a9Zg69atAIBx48ahb9++eJZXz2fo3r17w+fOnTsjJiYGlZWV6CCdDFSA\nL1hmN9iY0chxcfKbkpyEb7icxM5KKoxL7hwIJ7HyHE76pLGys9hthGeR2EFqKtdQ6h0NGsVMPV67\n1thOYqfkiEQsKYaqqqqGBj46OhqVGg5Wtm7dil9//RWxCm44cwUrkmlpaUhLS8Pll3MbvfQs+t1x\nhzOnUUmRFo577gkNV4UYaRjuvJPrzb73nnm5lE6kk0NvAY/EinDNNcb2fyjhp4bdCfNWM+g1k3bC\nTbtRlLzamuXss835JrK7HBUVFaGoqMjeSCVoKob09HQcPHgw7Pc5ej2NneHHH3/EhAkT8NprrymG\nyZUxVUlPD7c3VkLoCsDJSj16NFBYGPp+++3yisEIZ50FdO9uba/EVVfpD+uVO183EC40P/ec+EAg\nJ1Eqc61bc+ayZund297Fc7N1Y+fOcEeGbijPr7925jxvo50np0ZDSmkqTZ/znWae2XYP/6BDMXyg\nYrvVsWNHHDp0CNHR0aiqqkInBQcgx44dw/DhwzFnzhxc7oZdIZzdQn7WWeoeXXnc7nFq9dL4Befq\nauCf/9QXp596zXoRyjx9uviamyMg4fkMP/xgLo6yMncsqvSgZQwwbFjI0aKd+eyUyxO/I/RF5jaW\nms+srCzk5eUBAPLy8pAlo+JOnjyJUaNGYcKECbhezrm6Q0TiFIgdvPqq8rURI4CaGm6e3Skrr717\nzd135ZXurLM4iRNrDMFguG2+XxEugNu5sC3FrdGu0oluTuKXjpglxTB79mwUFBQgKSkJa9euxeOP\nPw4A2L59OyZNmgQAWL58OTZv3oylS5eid+/e6N27N0pKSqxLHgF48ZLVdtQGAqG56lGjlI9WtIJZ\nhdy/P1BRYa8sbuOXSu0mo0Zxa4BCLrmE2wzqBJ9/7twZBEL27AmdFWEn/fure7H1C5Z0b/v27WWn\nmvr06YM+Z3aKjBs3DuPGjbOSjCz33cftUvziC/nrblbSSB2dGHGfEElEosxuYmf+9Oolf5iOU9i5\nN0INO9c0hJvYPv5Y/31elmPfu8RQIiEBmDnTaykIJ9i3z2sJGheRqii//trZ+N3o0J0+bewIADc3\nnqoRsYqB8Cd2VDar3ij90hD6RQ7CO4wawTzzTGjPlpczEaQYHMTLhoEaJXm8tPTwCmkDo+Yh1iof\nf+y+6+/GRLt2gAnHELbTaBUDf350pM7/+4XGpmDuu0//M1kdufil7Amft7bWWWub/v3t2zfil/xz\nG/65jexLspuI3uakVnCkZyQ7Sffu8uZ5ja1R5Vm8WNnc1Y7KbDXf5s0Lt5RxWwa/EomHxjRV9G7s\ndYKIVgx+oXVrrrH0E072Cp20UbeDpCT3rFcigaba89biyiu9lsC/kGJwEK96ncuX2+cszi+7bt1g\n8GDgww/ti8+J99+YXZm4jVeuyiMBKmaNEKvTKEKefBKYOlV/+PPPt56+mlJ78UXn5l7tPPPaKcx4\nlInEabF77/XXGd1NjYhWDH4fIkdihZTSti2QmKg/fMuW1jc8nX22ct7deae1uNXw+/t67TX5s7ob\nI0aO3DWL39sPL4loxaAHevlEY2H8eHP3UR2ILJw0J9ZLozVX9QN+74ESYuyskEuWcIfX+AEqh5FF\nVJT37yyiFYPeQ0MIQg9DhwKbNtkT1+TJkbFmQRByRLRiuO46QOW4CM/xWusTxmjeHBgwwGsp7Iem\nkgijRLRiaN1a+0zXpmRuSRByUAeFMEqjXnw+ftzb4TxVSILwLwpHzxNo5IqB5ngJI1hV5H49aY2m\nksKhTps6ET2V5HeGDAEuu8xrKQi3SEnhztP2G9QIEkZp1CMGr1myxGsJCCNYbUADAftckRCEl5ge\nMRw5cgTp6elISkpCRkYGjh49GhbmwIEDSElJQe/evdGzZ088+eSTloQlCIIgnMe0YsjJycGwYcNQ\nUlKCzMxM5OTkhIW54IILsHXrVhQXF6O4uBiLFy/G/v37rchLEI5BUy4EwWFaMaxZswbjz+zRHzdu\nHAoKCsLCtGjRAlFn3EGeOHECLVq0QNu2bc0mSRAEQbiAacVQVVWFDmfMMKKjo1FZWSkbrry8HElJ\nSbjwwgsxY8YMtG/f3mySBEEQhAuoLj6np6fj4MGDYb/PmTNHdwLBYBAlJSX48ccfMXDgQAwZMgTd\nu3eXDZubm9vwOS0tDWlpabrTIQir0FQSEQkUFRWhqKjI0TQCjJmrDt26dcOWLVsQHR2Nqqoq9OvX\nD3v27FG9Z+LEiUhPT8eNN94YLkggAJOiEIRlAgHOimzyZK8lsZdAACgsBDIyvJaEcAon2k7TU0lZ\nWVnIy8sDAOTl5SErKysszI8//oja2loAQHV1NTZv3oyEhASzSRKEo1C/hCA4TCuG2bNno6CgAElJ\nSVi7di0ef/xxAMCOHTswadIkAMDu3buRmpqK5ORk9O/fH/fddx8pBoIgCJ9jeirJbmgqifCSQAB4\n5RVg4kSvJbEXmkpq/DjRdtLOZ4IAsH07kJTktRQE4Q9IMRAEyKcVQQghJ3oEQRCECFIMBEEQhAhS\nDARBEIQIUgwEQRCECFIMBEEQhAhSDARBEIQIUgwEQRCECFIMBEEQhAhSDARBEIQIUgwE0ciJjvZa\nAiLSIJcYBNGIqakBzjnHaymISIO8qxIEQUQwvjqohyAIgmickGIgCIIgRJBiIAiCIESQYiAIgiBE\nkGIgCIIgRJhWDEeOHEF6ejqSkpKQkZGBo0ePKoY9duwYgsEgpk+fbja5JkVRUZHXIvgGyosQlBch\nKC+cxbRiyMnJwbBhw1BSUoLMzEzk5OQohn300UcxcOBAs0k1OajQh6C8CEF5EYLywllMK4Y1a9Zg\n/PjxAIBx48ahoKBANtyOHTtQWVmJIUOGmE2KIAiCcBHTiqGqqgodOnQAAERHR6OysjIsTH19Pe6/\n/348/fTT5iUkCIIgXEV153N6ejoOHjwY9vucOXMwbtw4HDt2rOG3tm3bir4DwPPPP48TJ05g1qxZ\nWLp0KXbs2IFFixbJCxIImH0GgiCIJo3dO59VfSV98MEHitc6duyIQ4cOITo6GlVVVejUqVNYmM8+\n+wybN2/G4sWLcfz4cZw8eRJt2rTB3Llzw8KSOwyCIAh/YNpX0vTp09GtWzfce++9WLhwIfbt24fn\nnntOMfyyZcuwfft2xREDQRAE4Q9MrzHMnj0bBQUFSEpKwtq1a/H4448D4BabJ02aJHsPTRcRBEFE\nAMxj1q5dyxISElhsbCz729/+5rU4jvD999+zAQMGsISEBNazZ082f/58xhhjhw8fZtdccw1LTExk\nQ4YMYdXV1Q33TJ8+ncXFxbHevXuznTt3Nvy+dOlSFhcXx+Li4tiyZctcfxa7OHXqFOvVqxcbPnw4\nY4yx7777jvXt25clJCSw7OxsdvLkScYYY7/99hsbO3YsS0hIYFdccQXbv39/Qxxz585lsbGxLCEh\nga1bt86T57BKdXU1GzNmDEtKSmKXXnop+/TTT5tsuXjsscdYjx492B/+8Ad2/fXXs19++aXJlIvb\nbruNderUiSUkJDT8Zmc52L59O+vVqxeLi4tjd999t6Y8niqG3377jV188cWsvLyc1dXVsT59+oge\nsrFw8OBBtnv3bsYYYzU1NaxHjx5s165dbNq0aWzhwoWMMcYWLlzY8MLy8/PZtddeyxhjbOfOnSw5\nOZkxxtgPP/zAunXrxmpqalhNTQ3r1q0bO3jwoAdPZJ2nn36a3XzzzWzEiBGMMcaGDx/OVqxYwRhj\n7J577mELFixgjDH297//nd1zzz2MMcZWrFjBRo4cyRjjCnqfPn3YqVOnWHl5Obv44otZbW2tB09i\njTFjxrA333yTMcbY6dOn2c8//9wky8W3337Lunbt2vAOx44dy1555ZUmUy42bdrEdu7cKVIMdpSD\nn376iTHGWGJiYkPbeu2117L33ntPVR5PFcPGjRvZsGHDGr4/9dRT7IknnvBQIne4/vrrWUFBAbvk\nkkvYoUOHGGOMVVVVsW7dujHGuN5Dfn5+Q/j4+HhWVlbGli1bxqZNm9bw+9SpU9nrr7/urvA2UFZW\nxgYPHsw2bNjAhg8fzk6dOsWio6Mbrm/bto0NHjyYMcbYoEGD2Pbt2xljXMMZHR3NTp8+zWbPns3+\n/ve/N9wzbNgwtnnzZncfxCKHDh1i3bt3D/u9KZaLw4cPs549e7IjR46wuro6Nnz4cPaf//ynSZWL\nffv2iRSDXeXgwIEDLD4+vuH3d955h02cOFFVFk99JZWXlyMmJqbhezAYRHl5uYcSOc/+/fuxbds2\nXHnllYp7QSoqKmTzpaKiAsFgMOz3SGPGjBl46qmn0KwZV/wqKysRLTh/skuXLg3PJSwjzZo1Q4cO\nHVBZWdko8uLbb79Fx44dMXbsWCQkJGDChAmoqalpkuWiffv2uO+++3DhhRfiggsuwHnnnYeEhIQm\nWS547CoH0vDCfFTCU8XQ1Bajjx8/jjFjxuDZZ59F27ZtVcOyRmq+++9//xudOnVC7969G56xsT6r\nFvX19di2bRtmzZqF0tJStG/fHk888YTqPY01r/bu3YtnnnkG+/fvxw8//IDjx4+rmss3dZwuB54q\nhmAwiLKysobvZWVlIs3WmKirq8P111+PW265Bddddx2A0F4QAKK9INJ84XtHjSG/PvnkE6xatQpd\nu3bFTTfdhA0bNuDBBx9syAeAe16+5xMMBvH9998D4BrSw4cPo2PHjop5FEnExMSgS5cuSE1NBQCM\nGTMGu3btQqdOnZpcudi6dSuuuOIKdOjQAVFRURg9ejQ2bdrUJMsFj13tg1x44chCDk8VQ2pqKkpL\nS1FRUYG6ujosX74cmZmZXorkCIwxTJw4EXFxcZgxY0bD71lZWcjLywMA5OXlISsrq+H3N954AwCw\nc+dONG/eHF26dMHgwYNRWFiImpoa1NTUoLCwENdcc437D2SBuXPnoqysDPv27cNbb72FQYMG4fXX\nX0ffvn3x/vvvAwjPCz6PVq5ciX79+qF58+bIysrC22+/jVOnTqG8vBylpaW4/PLLPXsuM8TExCA6\nOhrffPMNAGD9+vWIjY1FZmZmkysX3bt3x2effYYTJ06AMYb169fj0ksvbZLlgseu9iEmJgbNmjVD\ncXExAOCNN95oiEsRa8sl1lmzZg2Lj49nsbGxbO7cuV6L4wibN29mgUCAJScns169erFevXqxtWvX\niszR0tPTReZoU6dObTBH27FjR8Pvr776KouNjWWxsbFs6dKlXjyObRQVFTVYJamZJd5www0sISGB\n9evXj+3bt6/h/jlz5rDY2FgWHx/PCgsLvXgEy+zatYv16dOHxcXFsczMTHbkyJEmWy5ycnJY9+7d\nWc+ePVl2djY7ceJEkykXN954I+vcuTNr0aIFCwaD7NVXX7W1HAjNVadPn64pj+mdzwRBEETjhE5w\nIwiCIESQYiAIgiBEkGIgCIIgRJBiIAiCIESQYiAIgiBEkGIgCIIgRPw/rJTRDS8aRRcAAAAASUVO\nRK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x9108ac90>"
       ]
      }
     ],
     "prompt_number": 63
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot([0,1], [trace['beta_0'].mean(), trace['beta_1'].mean()], 'ks-')\n",
      "t = pm.stats.hpd(trace)\n",
      "fill_between([0,1], [t['beta_0'][0], t['beta_1'][0]], [t['beta_0'][1], t['beta_1'][1]], color='grey')\n",
      "axis(ymin=-2, ymax=2, xmin=-.1, xmax=1.1);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD9CAYAAABHnDf0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3WlwVGX6NvDrZCGIsrgQ1HQLighkARoCEoE/rSUgwaDj\nAjKKio5DIYJbWaNVM2OoGnGkxqlxqUGcchsZHXFBsUBqnNHGQgdy+pzu7AkJpJPuhKyQfevleT84\n6RdIQprel+v3Kd39JOc5LBeHu+8+tySEECAioqgSF+oNEBGR/zHciYiiEMOdiCgKMdyJiKIQw52I\nKAox3ImIopBP4W61WvF///d/yMjIwPTp07Fjx44h123duhVpaWmYO3cuTCaTL4ckIiIPJPjyzaNG\njcJf//pXpKeno7OzE3PnzsWKFSswe/Zs95rPP/8cNTU1KC4uhslkwoYNG2A2m33eOBERDc+nK/dJ\nkyYhPT0dAHDJJZdg1qxZqKurO2vNgQMHsH79egCATqeDw+GAzWbz5bBERDQCv9XcLRYLZFnG4sWL\nz3reZrNBq9W6H2s0GoY7EVGA+VSWGdDZ2Yl7770Xr732GsaOHTvo9XPvcCBJ0qA1Qz1HREQjG+ou\nMj6Hu91ux913341f/vKXuPPOOwe9rtFoYLVaceONNwL4+Upeo9F4vMFAyc3NRW5ubtCOF0zRfG4A\nzy/Shfv5CSHQ0NAAVVVRUFAAIQT6+/s9+t7ExER0dHTg5ZdfDvAu/7/hLox9CnchBB599FGkpqbi\n6aefHnJNdnY2du/ejXvuuQeqqiI+Ph4pKSm+HJaIyO86OjpQUFAAo9GIrq4uOByOoF5w+ptP4f7j\njz9i9+7dmDVrFnQ6HQBg+/btqKmpAQBs3LgRd999N77//nukpaUhKSkJ7733nu+7JiLyA7vdjrKy\nMsiyjLq6OkiSBIfDEept+YVP4b548WK4XK4R17355pu+HCYg9Hp9qLcQMNF8bgDPL9KF+vyEEKiu\nroaiKCgrK0NcXJzHZRdPnNtUEipSuNzPXZKkiP4vEBGFt5aWFpjNZqiqCofD4ddAH5CYmIjNmzdj\n/Pjxfv/ZwxkuO/3SLUNEFI56enpQVFQEWZZx+vRpCCHgdDpDva2gYLgTUVRxOp2orKyE0WhEVVUV\n4uLiYLfbQ72toGO4E1HEE0Lg5MmTUFUVhYWFAOAuu8TKlfq5GO5EFLHa29vd7Yvd3d1wOp0eNXnE\nAoY7EUWU/v5+lJWVIS8vD/X19VHVvuhPDHciCntCCFgsFiiKgvLycr+3L0YjhjsRha3m5maYTCaY\nTCY4nU4G+gVguBNRWOnu7na3L7a2tsLlcrGO7gWGOxGFnNPpREVFBWRZRnV1dcy2L/oTw52IQkII\ngbq6OqiqiqKiIgBsX/QnhjsRBVVbWxvy8/OhKAp6enoi/u6L4YrhTkQB19/fj9LSUuTl5aGxsREA\n2L4YYAx3IgoIl8sFi8UCo9GIiooKti8GGcOdiPyqqakJJpMJZrOZ7YshxHAnIp91d3ejsLAQsiyj\nra2N7YthgOFORF5xOBzu9sWamhq2L4YZhjsReUwIgdraWiiKguLiYkiSxPbFMMVwJ6IRtba2utsX\ne3t72b4YAXwO90ceeQT79+9HcnKy+z7KZzIYDLjjjjtw3XXXAQDuvvtu/Pa3v/X1sEQUYH19fSgp\nKYEsy2hqaoqpKUbRwOdw37BhA7Zs2YIHH3xw2DVLly7Fvn37fD0UEQWYy+VCVVUVjEYjKisr2b4Y\nwXwO9yVLlsBisZx3Df/7RhTeGhsb3e2LLpeLgR4FAl5zlyQJ//3vf5GRkYHk5GT8+c9/xuzZswN9\nWCIaQVdXl7t9saOjg1OMokzAw33evHmw2WwYPXo0/vWvf+HOO+9EVVXVkGtzc3PdX+v1euj1+kBv\njyimOBwOHDt2DHl5ebDZbJxiFIEMBgMMBsOI6yThh5qJxWJBTk7OkG+onmv69Ok4dOgQrrzyyrM3\nIkks3xAFgBACNpsNiqKgpKTkrPZF8q/ExERs3rwZ48ePD9oxh8vOgF+5Nzc344orrgAAKIqCrq4u\nJCcnB/qwRDHv9OnTMJvNUBQFdrsddrudF1AxxOdwX7duHQ4dOoTm5mZotVps27bN/Sm1jRs34uOP\nP8bbb78NABg1ahQ++ugjxMXF+XpYIhpCb2+vu32xubmZ7YsxzC9lGX9gWYbIOy6XCydOnIDRaMTx\n48chSRJvAxAiMVWWIaLAaGhogKqqyM/PhxCCdXQ6C8OdKIJ0dnaioKAARqMRnZ2dbF+kYTHcicKc\n3W5HeXk5ZFlGbW0t2xfJIwx3ojAkhIDVanW3L/I2AHShGO5EYeTUqVMwm81QVZXti+QThjtRiPX2\n9qK4uBh5eXk4deoU2xfJLxjuRCHgdDpx/PhxGI1GnDhxglOMyO8Y7kRBIoRwty8WFBSc1b7IK3Xy\nN4Y7UYB1dHS42xe7uro4xYiCguFOFAB2ux1lZWWQZRl1dXVsX6SgY7gT+YkQAtXV1VAUBWVlZWxf\npJBiuBP5qKWlxd2+6HA4GOgUFhjuRF7o6elBUVERZFnG6dOn2b5IYYfhTuQhp9OJyspKGI1GVFVV\nsX2RwhrDneg8hBA4efIkVFV1Txpj+yJFAoY70RDa29vd7Yvd3d1sX6SIw3An+p/+/n6UlZUhLy8P\n9fX1bF+kiMZwp5gmhIDFYoGiKCgvL2f7IkUNhjvFpObmZphMJphMJjidTgY6RR2fw/2RRx7B/v37\nkZyc7H7D6Vxbt27Ff/7zHyQlJeGdd96BTqfz9bBEF6y7u9vdvtja2gqXy8UpRhS1fA73DRs2YMuW\nLXjwwQeHfP3zzz9HTU0NiouLYTKZsGHDBpjNZl8PS+QRp9OJiooKyLKM6upqti9SzPA53JcsWQKL\nxTLs6wcOHMD69esBADqdDg6HAzabDRqNxtdDEw1JCIG6ujqoqoqioiIAbF+k2BPwmrvNZoNWq3U/\n1mg0w4Z7bm6u+2u9Xg+9Xh/o7VEUaWtrQ35+PhRFQU9PD9sXKSoZDAYYDIYR1wXlDdVz/4JJkjTk\nujPDncgT/f39KC0tRV5eHhobGwGA7YsU1c698N22bduQ6wIe7hqNBlarFTfeeCMAsCRDPnO5XLBY\nLDAajaioqGD7ItEQAh7u2dnZ2L17N+655x6oqor4+HikpKQE+rAUhZqammAymWA2m9m+SDQCn8N9\n3bp1OHToEJqbm6HVarFt2zZ3N8LGjRtx99134/vvv0daWhqSkpLw3nvv+bxpih1dXV3u9sW2tja2\nLxJ5SBJh8o6TJEl884sA/FwzH2hfrKmpYfsiRYzExERs3rwZ48ePD9oxh8tOfkKVwoIQArW1tVAU\nBcXFxZAkie2LRD5guFNItba2utsXe3t72b5I5CcMdwq6vr4+lJSUQJZlNDU1cYoRUQAw3CkoXC4X\nqqqqYDQaUVlZyfZFogBjuFNANTY2QlVV5Ofnw+VyMdCJgoThTn7X1dWFwsJCyLKMjo4OOJ1Oti8S\nBRnDnfzC4XDg2LFjyMvLg81mY/siUYgx3MlrQgjYbDYoioKSkhK2LxKFEYY7XbDTp0/DbDZDURTY\n7XbY7Xa2LxKFGYY7eaS3t9fdvtjc3Mz2RaIwx3CnYblcLpw4ccLdvhgfH89uF6IIwXCnQRoaGtzt\ni0II1tGJIhDDnQAAnZ2dKCgogNFoRGdnJ9sXiSIcwz2G2e12lJeXQ5Zl1NbWQpIkTjEiihIM9xgj\nhIDVaoXRaERpaSlvA0AUpRjuMeLUqVMwm81QVZXti0QxgOEexXp7e1FcXIy8vDycOnWK7YtEMYTh\nHmWcTieOHz8Oo9GIEydO8DYARDGK4R4FhBBoaGiAoigoKCgAALYvEsW4OF9/wMGDB5GRkYHU1FS8\n8sorg15///33MXHiROh0Ouh0Orz77ru+HpL+p6OjAz/++CNef/11vPvuu1AUBf39/XyDlIh8u3Lv\n6+vDpk2bcPjwYUyaNAlZWVlYvnw5dDqde40kSVi3bh1ef/11nzdLP7cvlpWVQZZl1NXVsX2RiIbk\nU7gfPXoUaWlpSElJAQCsXbsW+/fvPyvchRDsyvCREALV1dVQFAVlZWVsXySiEfkU7jabDVqt1v1Y\no9HAYDCctUaSJHzxxRf47rvvMHXqVLzxxhuYPHnykD8vNzfX/bVer4der/dlexGvpaXF3b7ocDgY\n6EQEg8EwKGeH4lO4S5I04prVq1fj/vvvR0JCAt555x3cf//9OHz48JBrzwz3WNXT04OioiLIsozT\np0+zfZGIznLuhe+2bduGXOdTuGs0GlitVvdjq9V61pU8AFx66aXurx999FE89dRTvhwyKjmdTlRW\nVkKWZVgsFrYvEpHPfAr3+fPno6ioCLW1tUhOTsaePXuwa9eus9Y0NTVh4sSJAICvv/4a06ZN8+WQ\nUUMIgZMnT0JVVRQWFgJg+yIR+Y9P4T569Gjs3LkTK1asgMvlwvr16zF37ly8+OKLyMzMRE5ODl59\n9VUcOHAATqcTl156KT788EN/7T0itbe3u+++2N3dDYfDwTecicjvJBEmySJJUtSGXH9/P8rKypCX\nl4f6+nq2LxJFqcTERGzevBnjx48P2jGHy05+QjVAhBCwWCxQFAXl5eVsXySioGK4+1lzczNMJhNM\nJhOcTicDnYhCguHuB93d3e72xdbWVrhcLk4xIqKQYrh7yel0oqKiArIso7q6mu2LRBRWGO4XQAiB\nuro6qKqKoqIiAGxfJKLwxHD3QFtbG/Lz86EoCnp6eti+SERhj+E+jP7+fpSWliIvLw8NDQ1sXySi\niMJwP4PL5YLFYoHRaERFRQXbF4koYjHc8fMtEkwmE8xmM9sXiSgqxGy4d3V1udsX29ra2L5IRFEl\nJsP9q6++QmFhIdsXiShqxWS4l5eXw+l0sn2RiKKWzwOyiYgo/DDciYiiEMOdiCgKMdyJiKIQw52I\nKAox3ImIolBMtUI+/PDDsFgssFqtZ31gacKECbjzzjtDuDMiIv/yOdwPHjyI5557Dk6nEw899BB+\n85vfnPV6X18fHnzwQZSUlGDcuHH46KOPMHnyZF8P6xWLxYJDhw4Nen7KlCnB3wwRRY0vv/wSra2t\nkCQJ33zzDeLj4wH8nC3vv/9+SPbkU7j39fVh06ZNOHz4MCZNmoSsrCwsX74cOp3OvebNN9/EVVdd\nhU8++QRffvkltm7diq+++srnjfuTEAIOhwOSJI241pM1gVhHROGrtbUVFosFAFBVVRXazfyPT+F+\n9OhRpKWlISUlBQCwdu1a7N+//6xwP3DgAHbs2AEAWL16NR577DEIIcIq1Kqrq/Hyyy+f9Zyn92v3\n97rh+PMfi1D9AxWKdeG8N3+vC+e9ebounPd2vnWnTp3y6PuDyadwt9ls0Gq17scajQYGg2HYNXFx\ncbj88svR2NiISZMmDfp5ubm57q/1ej30er0v2/PYlClT8PDDDwflWCMJxT8qofqHjOfq3bpw3pu/\n14Xz3s5ct2/fPpw8edKj7/GVwWAYlLND8Snc/X31fWa4xyqWc4giT1JSUtCOde6F77Zt24Zc51O4\nazQaWK1W92Or1XrWlfzAmpqaGiQnJ8PlcqGlpQUTJ0705bBeG3jjdKhuGSKiaOJTuM+fPx9FRUWo\nra1FcnIy9uzZg127dp21Jjs7G7t370ZmZia++uorZGVlIS4uNO31A+9a79ixAz09PSHZAxFFnwkT\nJmDKlCmQJAkpKSlndcuEik/hPnr0aOzcuRMrVqyAy+XC+vXrMXfuXLz44ovIzMxETk4OnnjiCaxf\nvx4ZGRkYO3YsPvroI3/tnYgoLAx8TiYxMRGbN2/G+PHjQ7wjQBK+tnH4iSRJPneUeIpX7kQUCKEI\n9+Gyk7cfICKKQgx3IqIoxHAnIopCDHcioigUk+H+i1/8AhqNBvHx8UhIiKkbYxJRjIjJZJs2bRqm\nTZuG9vZ2FBQUwGg0oru7Gw6HI2gdO0REgRST4T5g3LhxWLx4MRYtWoSGhgaoqoqCggIIIdDf3x/q\n7REReS0m+9zPx+Vy4cSJEzAajTh+/DgkSYLdbg/1togoAoRTn3tMX7kPJS4uDtdffz2uv/569PX1\noaSkBLIso6mpCQDgcDhCvEMiopEx3M8jKSkJOp0OOp0ObW1tyM/Ph6Io6OnpYX2eiMIayzIXSAiB\nkydPQlEUFBUVAQDr80QEgGWZiCZJEq6++mpcffXVyM7ORmVlJYxGI6qqqhAXF8f6PBGFBYa7D+Lj\n4zF9+nRMnz4dvb29KC4uhizLaGlpgRACTqcz1FskohjFcPeT0aNHY968eZg3bx5Onz4Ns9kMRVFg\nt9tht9sjouRERNGDNfcAEkLAZrNBURSUlJRAkiTW54miGGvuMUKSJGi1Wmi1Wtx+++2oqKiALMuo\nqalhfZ6IAorhHiQJCQmYOXMmZs6cie7ubhQVFUGWZbS2trI+T0R+x3APgTFjxmDBggVYsGABWlpa\nYDKZYDKZ4HA4WJ8nIr9gzT1MCCFQU1MDRVFQWlqKuLg41ueJIkxU1NxPnTqFtWvXoqGhAVdddRU+\n+eQTTJgwYdC6+Ph4zJo1CwAwefJkfPnll94eMqpJkoTJkydj8uTJsNvtOHbsGGRZhs1mgyRJvO0B\nEV0Qr6/ct2zZgqlTp+Kpp57CX/7yF1RVVeG1114btG7s2LHo6OgYeSMxfuU+nK6uLhQWFkKWZbS3\nt8PlcsHlcoV6W0Q0hHC6cvc63KdOnYq8vDxcfvnlaG5uxsKFC1FZWTloHcPdf5qammAymWA2m+F0\nOlm2IQoz4RTuXpdlmpqacPnllwMArrjiCjQ2Ng65rre3F5mZmXC5XHj++eexZs2aYX9mbm6u+2u9\nXg+9Xu/t9qLSxIkTsXz5ctx6662wWCxQFAXHjh1jfZ4ohhgMBhgMhhHXnffKfdmyZaivrx/0/Esv\nvYQHHngA7e3t7ufGjRt31uMBjY2NSE5ORlVVFW655RYcPHgQ06dPH7wRXrl7pb+/H2VlZZBlGSdP\nnmR9niiEIubK/dtvvx32tYkTJ6K5uRlXXHEFmpqakJycPOS6geevvfZaLF++HKqqDhnu5J1Ro0Zh\n1qxZmDVrFjo6OtxjA7u6unhbYqIY5vWA7OzsbOzevRsAsHv3bmRnZw9a09bW5v4UZktLCw4dOoS0\ntDRvD0kjGDt2LBYtWoQnn3wSjzzyCObPn4+kpCSMGjUq1FsjoiDz+g3VM1shr7zySuzZswcTJkyA\noih466238Le//Q0//fQTNm7ciLi4OPT19WHr1q14/PHHh94IyzIBce7YQNbniQInnMoy/BBTDOnr\n60NpaSlkWXa/Ac76PJH/hFO48/YDMSQpKQlz5szBnDlzODaQKMrxyj3GDYwNVFUVhYWFADg2kMhb\nvHKnsHHm2MCVK1fi+PHjMBqNOHHiBG9LTBTBGO7kFh8fjxtuuAE33HADxwYSRTiGOw3p3LGBA/X5\nvr4+1ueJIgBr7uQxIQRqa2uhKAqKi4s5NpDoHKy5U0SSJAkajQYajQarVq1CRUUFjEYjqqurWZ8n\nCjMMd/LKuWMDB+rzp0+fZn2eKAww3MlnY8aMwfz58zF//ny0tLTAbDZDVVWODSQKIdbcKSAGxgaq\nqoqSkhLe9oBiQjjV3BnuFHAOhwPl5eXusYGsz1O0CqdwZ1mGAi4hIQFpaWlIS0tzjw00Go1oa2vj\n2ECiAGG4U1BdfPHFWLhwIRYuXMixgUQBxLIMhZwQwj02sLy8nPV5ilgsyxCdQZIkXHvttbj22mth\nt9vdYwPr6uo4NpDISwx3CiuJiYnIyMhARkaGe2ygoijo7OyE0+lkfZ7IQwx3ClsDYwMXLVqEhoYG\nqKqKgoICuFwulm2IRsCaO0UUl8uFqqoqGI1GVFZWsj5PYYU1dyIvxcXFYerUqZg6depZYwMbGhoA\ngLc9IPqfOG+/8dNPP0VaWhri4+Ohquqw6w4ePIiMjAykpqbilVde8fZwRIMMjA187LHHsHXrVixd\nuhTjxo1DYmIiJEkK9faIQsrrK/eMjAzs3bsXGzduHHZNX18fNm3ahMOHD2PSpEnIysrC8uXLodPp\nvD0s0ZDGjRuHJUuWYPHixTh58iRMJhMKCgoAcGwgxSavw33GjBkjrjl69CjS0tKQkpICAFi7di32\n79/PcKeAOXNs4G233caxgRSzAlpzt9ls0Gq17scajQYGg2HY9bm5ue6v9Xo99Hp94DZHUe/csYEl\nJSWQZRnNzc28LTFFLIPBcN4cHXDecF+2bBnq6+sHPb99+3bk5OSM+MMvtO55ZrgT+dPo0aMxd+5c\nzJ07F62tre7bEvf29nJsIEWUcy98t23bNuS684b7t99+69MmNBoNrFar+7HVaj3rSp4oFCZMmAC9\nXo+lS5eitrYWqqqiuLgYAOvzFD38UpYZ7qpn/vz5KCoqQm1tLZKTk7Fnzx7s2rXLH4ck8tlwYwMt\nFgvr8xTxvG6F3Lt3L7RaLY4cOYJVq1Zh5cqVAIC6ujqsWrUKwM//Fd65cydWrFiB2bNn46677sLc\nuXP9s3MiP4qPj8eMGTPwwAMP4Nlnn8WyZcswceJEJCQkID4+PtTbI7pg/IQq0Xm0tLQgPz8fqqrC\nbrdzbCCdVzh9QpXhTuQBIQSsVisURUFpaSkkSWJ9ngZhuA+B4U6RwuFw4NixY5BlGVarlbclJrdw\nCnfeW4boAiUkJCA1NRWpqano6upCUVERZFnm2EAKKwx3Ih9cfPHFuPHGG3HjjTeiubkZJpMJJpOJ\nYwMp5FiWIfIzIQSqq6thNBo5NjDGsCxDFMUkScKUKVMwZcoUjg2kkGG4EwXQuWMDCwsLYTQaOTaQ\nAo7hThQkY8eOxU033YSbbroJDQ0NMJlMyM/P59hACgjW3IlCaGBsoKIoqKioYH0+wrHmTkQAzh4b\n2N/f7x4bWF9fz/o8+YThThQmRo0ahdmzZ2P27Nlob29HQUEBjEYjuru7eVtiumAMd6IwNG7cOCxe\nvBiLFi1CfX09VFVFYWEhhBAs25BHWHMnihAul8s9NvD48eO8LXEYYs2diC5YXFwcpk2bhmnTpqG3\ntxelpaXIy8vj2EAaEsOdKAKNHj0aOp0OOp0Ora2tyM/Ph6IoHBtIbizLEEUJIQTq6uqgKArHBoYI\nyzJE5HeSJCElJQUpKSkcG0gMd6JoNDA2cMaMGejp6UFxcTFkWcapU6dYn48RXof7p59+itzcXPdN\nkYabjTplyhSMGzcO8fHxSExMRF5entebJaILd9FFFyEzMxOZmZk4deoUzGYzxwbGAK/DPSMjA3v3\n7sXGjRvPu06SJBgMBlx22WXeHoqI/OSyyy7DLbfcgptvvhlWqxWqqqKkpIRjA6OQ1+E+Y8YMj9fy\nyoAovEiShGuuuQbXXHMNbr/9do4NjEIBr7lLkoRly5bB4XDg17/+NZ544olAH5KILsBQYwONRiNa\nW1s5NjCCnTfcly1bhvr6+kHPb9++HTk5OR4d4MiRI0hOTkZTUxNuu+02zJgxA7feeuuQa3Nzc91f\n6/V66PV6j45BRP5x7tjAgfo8xwaGD4PBAIPBMOI6n/vcb775Zrz66qvDvqF6ppdffhkA8MILLwze\nCPvcicLSwNhARVFQVlbG2xKfR9T1uQ8Xyt3d3QCAMWPGoKurCwcPHsSzzz7rj0MSUZCcOzawvLwc\nsiyjtraW9fkwFuftN+7duxdarRZHjhzBqlWrsHLlSgBAXV0dVq1aBQCor69HVlYW5syZA51Oh6VL\nl2L16tX+2TkRBV1iYiLS09OxYcMGPPXUU7jllltw6aWXIjExEXFxXscJBQBvP0BEPmtsbITJZILZ\nbI7psYHhVJZhuBOR37hcLlgsFhiNxpgcGxhO4c7bDxCR38TFxeG6667Dddddx7GBIcZwJ6KA4NjA\n0GK4E1HAnTk2sKGhAaqqoqCggGMDA4g1dyIKiYGxgYqioLKyMipuS8yaOxHFvDPHBvb19aGkpASy\nLKOxsREAeFtiHzHciSjkkpKSODbQz1iWIaKwNDA2UFVVFBUVAQj/sYEsyxARjeDMsYHZ2dmorKyE\n0WhEVVVVVNTnA43hTkRhLz4+HtOnT8f06dPR09ODkpIS5OXlcWzgeTDciSiiXHTRRZg3bx7mzZvH\nsYHnwZo7EUU8IQRsNhsURQnp2EDW3ImI/EiSJGi1Wmi1Wtx+++2oqKiALMuoqamJ2fo8w52IokpC\nQgJmzpyJmTNnoru7G0VFRZBlGa2trTFVn2e4E1HUGjNmDBYsWIAFCxa4xwaaTCY4HI6wb6v0FWvu\nRBRTAjk2kDV3IqIQiZWxgQx3IopZA2MD09PT0dnZicLCQhiNRrS3t8PlcsHlcoV6i15juBMRAbjk\nkkuQlZWFrKysqBgb6PVE22eeeQapqalITU3F7bffjpaWliHXHTx4EBkZGUhNTcUrr7zi9Ub9zWAw\nhHoLARPN5wbw/CJdJJxfcnIyVqxYgeeeew5r167FzJkzkZCQgFGjRo34vYcPHw7CDkfmdbjn5OSg\nqKgIJSUlSE9Pxx/+8IdBa/r6+rBp0yYcPHgQBQUF+Oyzz2AymXzasL9Ewh8wb0XzuQE8v0gXSec3\nMDZwzZo1eO6557Bq1SpoNBrEx8cjIWHowke4hLvXZZmbb77Z/fWiRYvw4YcfDlpz9OhRpKWlISUl\nBQCwdu1a7N+/HzqdztvDEhGFxKhRozBr1izMmjULHR0d7rGBXV1dYXlbYr/U3N9++23cd999g563\n2WzQarXuxxqNJqL+1SYiGsrYsWOxaNEiLFq0CPX19e6xgX19faHemtt5+9yXLVuG+vr6Qc9v374d\nOTk5AICXXnoJqqri888/H7Tu448/xg8//ICdO3cCAP75z3/CYDDgrbfeGrwRSfL6JIiIYtkF97l/\n++235/2BH3zwAfbv34/vvvtuyNc1Gg2sVqv7sdVqPetKfqTNERGRd7x+Q/XgwYPYsWMH9u3bh9Gj\nRw+5Zv78+SgqKkJtbS3sdjv27NmDlStXer1ZIiLyjNfhvmXLFnR2dmLZsmXQ6XR4/PHHAQB1dXVY\ntWoVAGAtsJsFAAAE4ElEQVT06NHYuXMnVqxYgdmzZ+Ouu+7C3Llz/bNzIiIanohi33zzjUhPTxcz\nZ84Uf/zjHwe93tvbK9asWSPS09PFTTfdJCwWSwh26b2Rzm/Hjh0iNTVVpKWliSVLlogTJ06EYJfe\nG+n8Bnz22WdCkiShKEoQd+c7T87vk08+EXPmzBEZGRli3bp1Qd6hb0Y6v9LSUrFgwQKRlpYmZs6c\nKb788ssQ7NI7GzZsEMnJySI9PX3YNVu2bBGpqalCp9MJVVWDuLufRW249/b2iilTpgibzSbsdrvI\nzMwc9Av8pz/9STz55JNCCCH27t0rVq9eHYqtesWT8/vhhx9Eb2+vEEKInTt3ijvvvDMUW/WKJ+cn\nhBDt7e1iyZIlIisrK6LC3ZPzM5vNYsGCBaKzs1MIIURLS0sotuoVT87v/vvvF2+99ZYQQoiSkhKh\n0WhCsVWv/PDDD0JV1WHD/bPPPhN33HGHEEIIVVXF7Nmzg7k9IYQQXpdlwt2ZPfYJCQnuHvszHThw\nAOvXrwcArF69Gj/99FPEvLHryfktWbIESUlJAH7+LEJtbW0otuoVT84PAH73u9/h+eefR1JSUsT8\n3gGend97772HJ554AhdffDEA4LLLLgvFVr3iyflptVq0tbUBAFpbWzF58uRQbNUrS5YswaWXXjrs\n62dmi06ng8PhgM1mC9b2APhQcw93Q/XYn/uLe+aauLg4XH755WhsbAzqPr3lyfmdadeuXbjjjjuC\nsTW/8OT8VFVFbW0tsrOzAURWO60n51deXg6z2YzMzEzMmzcP+/btC/Y2vebJ+b3wwgv44IMPoNVq\nsWrVKrzxxhvB3mbAXOjfz0CI2huHRdJfdG9cyPn94x//gKqqOHToUAB35F8jnZ/L5cIzzzyDDz74\nwP1cJF25e/L753K5YLFYcPToUVitVtx0001YvHhxRFzBe3J+zzzzDH71q1/h6aefxpEjR/DAAw+g\nuLg4CLsLjnP/PAY7k6L2yt2THnuNRoOamhoAP/9FamlpwcSJE4O6T295+hmCf//733jppZewb98+\nJCYmBnOLPhnp/Do6OlBcXAy9Xo9rr70WR44cwerVq6Gqaii2e8E8+f3TarXIyclBfHw8pkyZgtTU\nVBw7dizYW/WKJ+d3+PBhrFmzBgCwcOFC9Pb2Rsz/nEdy7vnbbDZoNJrgbiLoVf4g6enpEZMnTxY2\nm0309/eLzMzMQW+4nfmG6hdffCFycnJCsVWveHJ+qqqKqVOnisrKyhDt0nuenN+Z9Hp9RL2h6sn5\nffHFF+Khhx4SQgjR1NQkrr76atHY2BiC3V44T84vOztbvP/++0KIn99QnTRpknA4HKHYrleqqqrO\n+4bqQAODoihi1qxZwdyaECKKu2WEEOLAgQPuNqvt27cLIYT4/e9/L/bt2yeE+Pkd/XvvvVekp6eL\nrKwsUVVVFcLdXrjhzu/rr78WQghx6623iiuvvFLMmTNHzJkzx/3ufaQY6ffvTJEW7kJ4dn7PPPOM\nSE1NFdOnTxd///vfQ7VVr4x0fmVlZWLhwoUiNTVVzJw50/3nNhLcd9994qqrrhKJiYlCo9GId955\nR7z11lvu7h8hhNi8ebO7FTIUfzbDZoYqERH5T9TW3ImIYhnDnYgoCjHciYiiEMOdiCgKMdyJiKIQ\nw52IKAr9P5rvttIasy+wAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x88e24250>"
       ]
      }
     ],
     "prompt_number": 64
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Interval estimates without MCMC\n",
      "\n",
      "It is appealing sometimes to avoid using MCMC. How can we get an interval estimate on `beta_0` and `beta_1` with a MAP estimate machine?\n",
      "\n",
      "One idea is to use bootstrap to represent the uncertainty in the data.  This bootstrap is not going to get the uncertainty on `beta_1` right, however, because it does not capture the uncertainty in the prior."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%time\n",
      "\n",
      "X_rep = []\n",
      "\n",
      "for rep in range(100):\n",
      "    y_resampled = y[randint(low=0, high=100, size=100)]\n",
      "    \n",
      "    with pm.Model() as m_rep:\n",
      "        beta_0 = pm.Flat('beta_0')\n",
      "        beta_1 = pm.Normal('beta_1', beta_0, 1.)\n",
      "        sigma = pm.Uniform('sigma', 0., 10.)\n",
      "        likelihood = pm.Normal('likelihood', beta_0, sigma**-2, observed=y_resampled)\n",
      "        \n",
      "        X_rep.append(pm.find_MAP())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 15min 20s, sys: 1min 26s, total: 16min 46s\n",
        "Wall time: 17min 17s\n"
       ]
      }
     ],
     "prompt_number": 65
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X_rep = pd.DataFrame(X_rep)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 66
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot([0,1], [X_rep['beta_0'].mean(), X_rep['beta_1'].mean()], 'ks-')\n",
      "t={}\n",
      "t[0] = pm.stats.hpd(array(X_rep['beta_0']))\n",
      "t[1] = pm.stats.hpd(array(X_rep['beta_1']))\n",
      "fill_between([0,1], [t[0][0], t[1][0]], [t[0][1], t[1][1]], color='grey')\n",
      "axis(ymin=-2, ymax=2, xmin=-.1, xmax=1.1);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD9CAYAAABHnDf0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFd1JREFUeJzt3X9s1dX9x/HXbYuQbSiItLjeG2rIBtzbFm4thDqq183C\naNdKJMIYVkWXEHSwSWImf2y0ibBBZjInW9GFARvgRH6JoTZj0VvCGGzhgg5wW8zacG8JA2pQQekK\nPd8//HJDvfe2t/fTe0sPz0dCcj+f+773vI8tLz+ce+69LmOMEQDAKlkD3QAAoP8R7gBgIcIdACxE\nuAOAhQh3ALAQ4Q4AFnIU7uFwWPfee6+Kioo0fvx4rVmzJm7d0qVL5fP5VFJSoqNHjzoZEgCQhBwn\nD77lllv0m9/8RoWFhbp48aJKSko0c+ZMTZo0KVqzY8cOnTp1SidOnNDRo0e1cOFCHTt2zHHjAIDE\nHF255+XlqbCwUJL0la98RcXFxTp9+nS3msbGRtXW1kqS/H6/rly5okgk4mRYAEAv+m3NvbW1VX//\n+981ffr0bucjkYg8Hk/02O12E+4AkGaOlmWuuXjxoh5++GG9+OKLGj58eMz9X/yEA5fLFVMT7xwA\noHfxPkXG8ZV7Z2en5syZo+9973uaPXt2zP1ut1vhcDh6HIlE5Ha7EzaYqT8rVqzI6HjMjfkxP+aX\njj+JOAp3Y4yefPJJeb1ePfPMM3FrKisrtWXLFklSKBRSdna28vPznQwLAOiFo2WZv/zlL9q8ebOK\ni4vl9/slSatWrdKpU6ckSYsWLdKcOXP0zjvvyOfzaejQodqwYYPzrgEAPXIU7tOnT1dXV1evdWvX\nrnUyTFoEAoGBbiFtbJ6bxPwGO+aXGS7T06JNBrlcrh7XjwAAsRJlJx8/AAAWItwBwEKEOwBYiHAH\nAAsR7gBgIcIdACxEuAOAhQh3ALAQ4Q4AFiLcAcBChDsAWIhwBwALEe4AYCHCHQAsRLgDgIUIdwCw\nEOEOABYi3AHAQo7D/YknnlBeXp6Kiori3h8MBnXbbbfJ7/fL7/fr+eefdzokAKAXjr4gW5IWLlyo\nJUuW6NFHH01Yc99992nPnj1OhwIAJMnxlXt5eblGjhzZYw1ffA0AmZX2NXeXy6W//vWvKioq0re+\n9S29++676R4SAG56jpdlenP33XcrEolo2LBh+tOf/qTZs2erpaUlbm1dXV30diAQUCAQSHd7ADCo\nBINBBYPBXutcph/WTFpbW1VdXa1//OMfvdaOHz9ezc3NGjNmTPdGXC6WbwCgjxJlZ9qXZc6fPx+9\nfeTIEV26dEm5ubnpHhYAbmqOl2Xmz5+v5uZmnT9/Xh6PR/X19ers7JQkLVq0SK+++qpeeeUVSdIt\nt9yirVu3KiuL7fUAkE79sizTH1iWAYC+G7BlGQBA5hHuAGAhwh0ALES4A4CFCHcAsBDhDgAWItwB\nwEKEOwBYiHAHAAsR7gBgIcIdACxEuAOAhQh3ALAQ4Q4AFiLcAcBChDsAWIhwBwALEe4AYCHCHQAs\n5Djcn3jiCeXl5amoqChhzdKlS+Xz+VRSUqKjR486HRIA0AvH4b5w4UI1NTUlvH/Hjh06deqUTpw4\nofXr12vhwoVOhwQA9MJxuJeXl2vkyJEJ729sbFRtba0kye/368qVK4pEIk6HBQD0ICfdA0QiEXk8\nnuix2+1WJBKR2+2Oqa2rq4veDgQCCgQC6W4PAAaVYDCoYDDYa13aw12SjDHdjl0uV9y668MdABDr\nixe+9fX1cevSvlvG7XYrHA5HjxNdtQMA+k/aw72yslJbtmyRJIVCIWVnZys/Pz/dwwLATc3xssz8\n+fPV3Nys8+fPy+PxqL6+Xp2dnZKkRYsWac6cOXrnnXfk8/k0dOhQbdiwwXHTAICeucwXF8QHiMvl\nilmbBwD0LFF28g5VALAQ4Q4AFiLcAcBChDsAWIhwBwALEe4AYCHCHQAsRLgDgIUIdwCwEOEOABYi\n3AHAQoQ7AFiIcAcACxHuAGAhwh0ALES4A4CFCHcAsBDhDgAWItwBwEKOw72pqUlFRUXyer1avXp1\nzP0bN27U6NGj5ff75ff79bvf/c7pkACAXuQ4eXBHR4cWL16sAwcOKC8vT2VlZZoxY4b8fn+0xuVy\naf78+frVr37luFkAQHIcXbkfPnxYPp9P+fn5ysnJ0bx587R3795uNcaYuN/MDQBIH0dX7pFIRB6P\nJ3rsdrsVDAa71bhcLu3cuVNvv/22xo0bp5deekljx46N+3x1dXXR24FAQIFAwEl7AGCdYDAYk7Px\nOAp3l8vVa01NTY0WLFignJwcrV+/XgsWLNCBAwfi1l4f7gCAWF+88K2vr49b52hZxu12KxwOR4/D\n4XC3K3lJGjlypHJyPv9/yJNPPql3333XyZAAgCQ4CvcpU6bo+PHjamtrU2dnp7Zt26ZZs2Z1qzl3\n7lz09ptvvqmvfe1rToYEACTB0bLMsGHD1NDQoJkzZ6qrq0u1tbUqKSnRihUrVFpaqurqar3wwgtq\nbGzU1atXNXLkSP3hD3/or94BAAm4zA2ylcXlcrGrBgD6KFF28g5VALAQ4Q4AFiLcAcBChDsAWIhw\nBwALEe4AYCHCHQAsRLgDgIUIdwCwEOEOABYi3AHAQoQ7AFjI0adCDla//vWv1d7ePtBtALBMVlaW\nli1bpi996UsD3crNGe6XLl3iEygB9LusrCx1dnYOdBuSWJYBACsR7gBgIcIdACxEuAOAhQh3ALCQ\n4+9QbWpq0rPPPqurV6/qscce049//ONu93d0dOjRRx/VyZMndeutt2rr1q0aO3ZsbCMZ+A7Vxx9/\nXK2trQqHw+rq6oqeHzFihGbPnp3WsQHYa/fu3bpw4YJcLpfy8/OVnZ0tSSooKNDGjRvTOnai7HS0\nFbKjo0OLFy/WgQMHlJeXp7KyMs2YMUN+vz9as3btWt1555167bXXtHv3bi1dulRvvPGGk2FT1tra\nqubm5pjzBQUFmW8GgDUuXLig1tZWSVJLS8vANvP/HC3LHD58WD6fT/n5+crJydG8efO0d+/ebjWN\njY2qra2VJNXU1OjgwYPsMQeANHN05R6JROTxeKLHbrdbwWAwYU1WVpZGjRqls2fPKi8vL+b56urq\norcDgYACgYCT9pLW2trabWwAuFEFg8GYnI3HUbi7XC4nD48xUAFbUFCgxx9/fEDGBjD4bdy4Mbos\nk25fvPCtr6+PW+doWcbtdiscDkePw+Fwtyv5azWnTp2SJHV1dam9vV2jR492MiwAoBeOrtynTJmi\n48ePq62tTbm5udq2bZtefvnlbjWVlZXavHmzSktL9cYbb6isrExZWQOzA/PaC6fxdssAQKpGjBih\ngoKCuLtlBorjrZBvvfWWnn32WXV1dam2tlbLly/XihUrVFpaqurqanV0dKi2tlbvv/++hg8frq1b\nt8adcCa2Ql6zZs0affbZZxkZC8DNY8iQIXr66ad12223ZWzMRNnpONz7C+EOYLC7kcKdd6gCgIUI\ndwCwEOEOABYi3AHAQoQ7AFjopgz3G+HLawHY5+rVqxoyZMhAtyHpJt0KCQC2YCskANxECHcAsBDh\nDgAWItwBwEKEOwBYiHAHAAsR7gBgIcIdACxEuAOAhQh3ALAQ4Q4AFiLcAcBCKYf7hx9+qIqKChUX\nF2vmzJm6cOFC3Lrs7Gz5/X75/X7Nnj075UYBAMlL+VMhlyxZonHjxulHP/qRfvnLX6qlpUUvvvhi\nTN3w4cP1ySef9N4InwoJAH2WKDtTDvdx48bpb3/7m0aNGqXz589r2rRp+uCDD2LqCHcASJ9E2ZmT\n6hOeO3dOo0aNkiTdcccdOnv2bNy6y5cvq7S0VF1dXXruuec0d+7chM9ZV1cXvR0IBBQIBFJtDwCs\nFAwGFQwGe63r8cq9oqJCZ86ciTm/cuVKPfLII/r444+j52699dZux9ecPXtWubm5amlp0Te/+U01\nNTVp/PjxsY1w5Q4AfZbSlfu+ffsS3jd69GidP39ed9xxh86dO6fc3Ny4ddfO33XXXZoxY4ZCoVDc\ncAcA9J+Ud8tUVlZq8+bNkqTNmzersrIypuajjz5SZ2enJKm9vV3Nzc3y+XypDgkASFLKL6h++OGH\nmjdvnv773/9qzJgx2rZtm0aMGKEjR45o3bp1+u1vf6uDBw9q0aJFysrKUkdHh5YuXaqnnnoqfiMs\nywBAn/X7bpn+RrgDQN/xBdkAcBMh3AHAQoQ7AFiIcAcACxHuAGAhwh0ALES4A4CFCHcAsBDhDgAW\nItwBwEKEOwBYiHAHAAsR7gBgIcIdACxEuAOAhQh3ALAQ4Q4AFiLcAcBChDsAWCjlcH/99dfl8/mU\nnZ2tUCiUsK6pqUlFRUXyer1avXp1qsMBAPog5XAvKirSrl27dO+99yas6ejo0OLFi9XU1KT33ntP\n27dv19GjR1MdEgCQpJTDfcKECfr617/eY83hw4fl8/mUn5+vnJwczZs3T3v37k11SABAknLS+eSR\nSEQejyd67Ha7FQwGE9bX1dVFbwcCAQUCgfQ1BwCDUDAY7DFHr+kx3CsqKnTmzJmY86tWrVJ1dXWv\nT+5yuXqtud714Q4AiPXFC9/6+vq4dT2G+759+xw14Xa7FQ6Ho8fhcLjblTwAID36ZSukMSbu+SlT\npuj48eNqa2tTZ2entm3bplmzZvXHkACAHqQc7rt27ZLH49GhQ4dUVVUVDe3Tp0+rqqpKkjRs2DA1\nNDRo5syZmjRpkh566CGVlJT0T+cAgIRcJtFld4a5XK6E/wIAAMSXKDt5hyoAWIhwBwALEe4AYCHC\nHQAsRLgDgIUIdwCwEOEOABYi3AHAQoQ7AFiIcAcACxHuAGAhwh0ALES4A4CFCHcAsBDhDgAWItwB\nwEKEOwBYiHAHAAulHO6vv/66fD6fsrOzFQqFEtYVFBSouLhYfr9fU6dOTXU4AEAf5KT6wKKiIu3a\ntUuLFi3qsc7lcikYDOr2229PdSgAQB+lHO4TJkxIupYvvgaAzEr7mrvL5VJFRYWKi4u1du3adA8H\nAFAvV+4VFRU6c+ZMzPlVq1apuro6qQEOHTqk3NxcnTt3Tt/+9rc1YcIEPfDAA3Fr6+rqorcDgYAC\ngUBSYwDAzSIYDCoYDPZa5zIO10zuv/9+vfDCCyopKem19mc/+5kkafny5bGNuFws3wBAHyXKzn5Z\nlkkUyp9++qk+/fRTSdKlS5fU1NQkn8/XH0MCAHqQcrjv2rVLHo9Hhw4dUlVVlWbNmiVJOn36tKqq\nqiRJZ86cUVlZmSZPniy/36/77rtPNTU1/dM5ACAhx8sy/YVlGQDou7QuywAAbiyEOwBYiHAHAAsR\n7gBgIcIdACxEuAOAhQh3ALAQ4Q4AFiLcAcBChDsAWIhwBwALEe4AYCHCHQAsRLgDgIUIdwCwEOEO\nABYi3AHAQoQ7AFiIcAcAC6Uc7suWLZPX65XX69V3vvMdtbe3x61rampSUVGRvF6vVq9enXKj/S0Y\nDA50C2lj89wk5jfYMb/MSDncq6urdfz4cZ08eVKFhYV6/vnnY2o6Ojq0ePFiNTU16b333tP27dt1\n9OhRRw33lxvlB5AONs9NYn6DHfPLjJTD/f7771dW1ucP/8Y3vqG2traYmsOHD8vn8yk/P185OTma\nN2+e9u7dm3q3AICk9Mua+yuvvKIHH3ww5nwkEpHH44keu91uRSKR/hgSANADlzHGJLqzoqJCZ86c\niTm/atUqVVdXS5JWrlypUCikHTt2xNS9+uqr2r9/vxoaGiRJf/zjHxUMBrVu3brYRlyulCcBADez\neDGe09MD9u3b1+MTbtq0SXv37tXbb78d9363261wOBw9DofD3a7ke2sOAJCalJdlmpqatGbNGu3Z\ns0fDhg2LWzNlyhQdP35cbW1t6uzs1LZt2zRr1qyUmwUAJCflcF+yZIkuXryoiooK+f1+PfXUU5Kk\n06dPq6qqSpI0bNgwNTQ0aObMmZo0aZIeeughlZSU9E/nAIDEjMXeeustU1hYaCZOnGh+/vOfx9x/\n+fJlM3fuXFNYWGjuuece09raOgBdpq63+a1Zs8Z4vV7j8/lMeXm5+c9//jMAXaaut/lds337duNy\nucyRI0cy2J1zyczvtddeM5MnTzZFRUVm/vz5Ge7Qmd7m9/7775upU6can89nJk6caHbv3j0AXaZm\n4cKFJjc31xQWFiasWbJkifF6vcbv95tQKJTB7j5nbbhfvnzZFBQUmEgkYjo7O01paWnMf+Bf/OIX\n5oc//KExxphdu3aZmpqagWg1JcnMb//+/eby5cvGGGMaGhrM7NmzB6LVlCQzP2OM+fjjj015ebkp\nKysbVOGezPyOHTtmpk6dai5evGiMMaa9vX0gWk1JMvNbsGCBWbdunTHGmJMnTxq32z0QraZk//79\nJhQKJQz37du3mwcffNAYY0woFDKTJk3KZHvGGGOs/fiBZPbYNzY2qra2VpJUU1OjgwcPDpoXdpOZ\nX3l5uYYOHSop8XsRblTJvkfiJz/5iZ577jkNHTp00PzspOTmt2HDBv3gBz/Ql7/8ZUnS7bffPhCt\npiSZ+Xk8Hn300UeSpAsXLmjs2LED0WpKysvLNXLkyIT3X58tfr9fV65cyfg2cGvDPZk99tfXZGVl\nadSoUTp79mxG+0xVX99D8PLLL8d9L8KNKpn5hUIhtbW1qbKyUtLg2k6bzPz+9a9/6dixYyotLdXd\nd9+tPXv2ZLrNlCUzv+XLl2vTpk3yeDyqqqrSSy+9lOk20+ZGeI9Pj1shB7PB9Bc9FX2Z35YtWxQK\nhdTc3JzGjvpXb/Pr6urSsmXLtGnTpui5wXTlnszPr6urS62trTp8+LDC4bDuueceTZ8+fVBcwScz\nv2XLlun73/++nnnmGR06dEiPPPKITpw4kYHuMuOLv4+ZziRrr9yT2WPvdrt16tQpSZ//RWpvb9fo\n0aMz2meqkn0PwZ///GetXLlSe/bs0ZAhQzLZoiO9ze+TTz7RiRMnFAgEdNddd+nQoUOqqalRKBQa\niHb7LJmfn8fjUXV1tbKzs1VQUCCv16t///vfmW41JcnM78CBA5o7d64kadq0abp8+fKg+Zdzb744\n/0gkIrfbndkmMr7KnyGfffaZGTt2rIlEIuZ///ufKS0tjXnB7foXVHfu3Gmqq6sHotWUJDO/UChk\nxo0bZz744IMB6jJ1yczveoFAYFC9oJrM/Hbu3Gkee+wxY4wx586dM1/96lfN2bNnB6DbvktmfpWV\nlWbjxo3GmM9fUM3LyzNXrlwZiHZT0tLS0uMLqtc2MBw5csQUFxdnsjVjjMW7ZYwxprGxMbrNatWq\nVcYYY37605+aPXv2GGM+f0X/4YcfNoWFhaasrMy0tLQMYLd9l2h+b775pjHGmAceeMCMGTPGTJ48\n2UyePDn66v1g0dvP73qDLdyNSW5+y5YtM16v14wfP978/ve/H6hWU9Lb/P75z3+aadOmGa/XayZO\nnBj9vR0Mvvvd75o777zTDBkyxLjdbrN+/Xqzbt266O4fY4x5+umno1shB+J3s8fPlgEADE7WrrkD\nwM2McAcACxHuAGAhwh0ALES4A4CFCHcAsND/AY6gxGO8yXlIAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xe56e1590>"
       ]
      }
     ],
     "prompt_number": 67
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# (Semi)Parametric bootstrap\n",
      "\n",
      "There is a different bootstrap, where I think it will get things right.  In addition to including a resampling of the data, I will include a parametric offset of the conditional prior `beta_1 | beta_0`.\n",
      "\n",
      "In mathematics, I will change the prior on $\\beta_1$ for each bootstrap replicate, to be\n",
      "\n",
      "$ \\beta_1 \\sim \\text{Normal}(\\beta_0 + \\text{random noise}, 1^2) .$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%time\n",
      "\n",
      "X_rep_2 = []\n",
      "\n",
      "for rep in range(100):\n",
      "    y_resampled = y[randint(low=0, high=100, size=100)]\n",
      "    \n",
      "    with pm.Model() as m_rep_2:\n",
      "        beta_0 = pm.Flat('beta_0')\n",
      "        beta_1 = pm.Normal('beta_1', beta_0 + randn(), 1.)  # note the randn() added in here\n",
      "        sigma = pm.Uniform('sigma', 0., 10.)\n",
      "        likelihood = pm.Normal('likelihood', beta_0, sigma**-2, observed=y_resampled)\n",
      "        \n",
      "        X_rep_2.append(pm.find_MAP())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 26min 20s, sys: 1min 38s, total: 27min 59s\n",
        "Wall time: 30min 21s\n"
       ]
      }
     ],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X_rep_2 = pd.DataFrame(X_rep_2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot([0,1], [X_rep_2['beta_0'].mean(), X_rep_2['beta_1'].mean()], 'ks-')\n",
      "t={}\n",
      "t[0] = pm.stats.hpd(array(X_rep_2['beta_0']))\n",
      "t[1] = pm.stats.hpd(array(X_rep_2['beta_1']))\n",
      "fill_between([0,1], [t[0][0], t[1][0]], [t[0][1], t[1][1]], color='grey')\n",
      "axis(ymin=-2, ymax=2, xmin=-.1, xmax=1.1);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD9CAYAAABHnDf0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtwVOX9BvDn7CaEyl0kAbMZwjWQC5AQkKDIegnUYICR\nKqWKFGnLIEIVh1Fn2hpmKq1M7dTqFLDjiC1qxQuIA2Zqq4tDFXQIZC+538huQsgdSEI2e3l/f+ie\nX65k2fuefT4zzGR3X/a8xySPL+f7nu9KQggBIiJSFFWwJ0BERL7HcCciUiCGOxGRAjHciYgUiOFO\nRKRADHciIgXyKtzNZjPuvvtupKWlISkpCfv27Rt03M6dO5GSkoKMjAycP3/em0MSEZEborz5yyNG\njMDf/vY3pKamoqOjAxkZGVi5ciXmz58vj/noo49QW1sLk8mE8+fPY/Pmzbhw4YLXEycioqF5tXKP\ni4tDamoqAGD06NGYN28e6uvr+4w5efIkNm7cCABIT0+H3W6HxWLx5rBERDQMn11zr6mpwXfffYe7\n7rqrz/MWiwUJCQnyY41Gw3AnIvIzry7LuHR0dODhhx/Gq6++ijFjxgx4vX+HA0mSBowZ7DkiIhre\nYF1kvF6522w2rFu3Dj/72c+wdu3aAa9rNBqYzWb5scVigUajGXKCgfrz4osvBvR4PDeeH8+P5+eP\nP0PxKtyFENiyZQuSk5PxzDPPDDomJycH77zzDgCgoKAAarUa8fHx3hyWiIiG4dVlmf/97384fPgw\n5s2bh/T0dADA3r17UVtbCwDYunUr1q1bhy+//BIpKSmIiYnBW2+95f2siYjohrwK97vuugtOp3PY\nca+//ro3h/ELrVYb7Cn4jZLPDeD5hTueX2BI4kYXbQJIkqQbXj8iIqKBhspOth8gIlIghjsRkQIx\n3ImIFIjhTkSkQAx3IiIFYrgTESkQw52ISIEY7kRECsRwJyJSIIY7EZECMdyJiBSI4U5EpEAMdyIi\nBWK4ExF5SQiBxsZGFBQUuNUGPRB88hmqRESRxmazoaamBsXFxSgtLYXNZoPdbseMGTMwbty4YE+P\n4U5E5K6rV6+irKwMRqMRFosFUVFRsFqt8uvR0dFBnF1fDHcioiE4nU7U1dWhtLQUJpMJ165dg0ql\ngs1mAwA4HI4gz3BoDHciol6uX7+OyspKmEwmVFZWQpIk2Gw2+dOOQjnQe/M63J944gmcOHECsbGx\nMBgMA17X6XRYs2YNpk+fDgBYt24dfvOb33h7WCIinxBCoKmpSV6dNzc3Q61Wo6enJ9hT84rX4b55\n82bs2LEDjz/++JBjli9fjuPHj3t7KCIinxisGOp0OuVVebiszm/E63BftmwZampqbjiGH3xNRMF2\n5coVlJeXD1kMVRq/X3OXJAnffPMN0tLSEBsbiz//+c+YP3++vw9LRBHO6XTCYrGgtLQURUVF6Ojo\nkK+fA8pYnd+I38N94cKFsFgsGDlyJP79739j7dq1qK6uHnRsXl6e/LVWq4VWq/X39IhIQa5fv46K\nigoUFRUNWgxVAp1OB51ON+w4SfjgrGtqapCbmztoQbW/pKQknDp1CpMnT+47EUlS1DeAiPwv1Iqh\n0dHR2L59e0BvYhoqO/2+cm9ubsZtt90GADh37hw6OzsRGxvr78MSkUK5iqFFRUUoLS2F3W5XXDHU\nF7wO9w0bNuDUqVNobm5GQkIC9uzZI1/T2rp1K9577z288cYbAIARI0bg3XffhUrFljZE5D5XMdRg\nMKCurk7xxVBf8MllGV/gZRkicnEVQ0tKSlBUVITOzs4+xdBQFVGXZYiI3OEqhppMJlRVVSmyGBpI\nDHciCgpXm1xXI66WlhZF3BkaKhjuRBQwNpsN1dXV8p2hLIb6D8OdiPzqypUr8uqcxdDAYbgTkU8N\nVwzl6jwwGO5E5DUWQ0MPw52IbpqrGOq6M5TF0NDDcCcit7iKoUVFRSgrK2MxNMQx3IloSO3t7XKb\nXBZDwwvDnYhkTqcTZrNZbpPLYmj4YrgTRbiuri65TS6LocrBcCeKMCyGRgaGO1EEYDE08jDciRTK\nVQw1GAyor69nMTTCMNyJFMJVDC0pKUFxcTGLoRGO4U4UxlzFUJPJhOrqahZDScZwJwojLIaSuxju\nRCHOZrOhqqoKxcXFLIaS2xjuRCGovb1dbpPLYih5wutwf+KJJ3DixAnExsbCYDAMOmbnzp3473//\ni5iYGLz55ptIT0/39rBEisJiKPma1+G+efNm7NixA48//vigr3/00Ueora2FyWTC+fPnsXnzZly4\ncMHbwxKFvf7FUJVKhZ6eHhZDySe8Dvdly5ahpqZmyNdPnjyJjRs3AgDS09Nht9thsVig0Wi8PTRR\nWBFC4PLly/LlltbWVhZDyW/8fs3dYrEgISFBfqzRaIYM97y8PPlrrVYLrVbr7+kR+VVPT4/8maEs\nhpIv6HQ66HS6YccFpKDa/5+ZkiQNOq53uBOFq7a2NrlNLouh5Gv9F7579uwZdJzfw12j0cBsNuOO\nO+4AAF6SIcVxOBx9iqFdXV0AALvdLr9OFGh+D/ecnBwcPnwYP/nJT1BQUAC1Wo34+Hh/H5bIr7q6\nulBeXg6TyYSamhreGUohx+tw37BhA06dOoXm5mYkJCRgz5498vatrVu3Yt26dfjyyy+RkpKCmJgY\nvPXWW15PmijQXMVQ152hLIZSqJNEiCw1JEniqodCiqsYWlRUhPLy8gHFUKL+oqOjsX37dowbNy5g\nxxwqO3mHKlEvLIaSUjDcKaKxGEpKxXCniNPZ2SnfGcpiKCkVw50Uj8VQikQMd1Kk3sXQsrIyOBwO\n3hlKEYXhTorR1tYm9225dOkSi6EU0RjuFLZYDCUaGsOdwkrvYmh1dbV87ZzFUKK+GO4U0oQQaGho\nkIuhbW1tfYqhrlU6EfXFcKeQ09PT0+czQ12FUBZDidzHcKeQwGIokW8x3CkoXMXQ4uJiFBcX4/r1\n6wBYDCXyFYY7BUxnZ2efNrkqlYp3hhL5CcOd/Ga4YigR+Q/DnXzKVQx1tcllMZQoOBju5LXW1laU\nlZXBZDKxGEoUIhjudNMcDgdqa2vlO0NZDCUKPQx3ckv/YijvDCUKbQx3GhTvDCUKb16He35+Pnbv\n3g2Hw4FNmzbhueee6/P6oUOHsHv3bmg0GgDAjh078MQTT3h7WPIDq9Uq3xnKYihRePMq3K1WK7Zt\n24bTp08jLi4OWVlZWLFiBdLT0+UxkiRhw4YN+Otf/+r1ZMn3XMVQo9GIhoYGblUkUgivwv3s2bNI\nSUlBfHw8AGD9+vU4ceJEn3AXQvC6bAhhMZQoMngV7haLBQkJCfJjjUYDnU7XZ4wkSfj444/xxRdf\nYMaMGXjttdcwderUQd8vLy9P/lqr1UKr1XozPfpBR0dHn88MZTGUKHzpdLoBOTsYr8JdkqRhx6xe\nvRqPPvoooqKi8Oabb+LRRx/F6dOnBx3bO9zJc0IIXLp0SS6Gtre3sxhKpBD9F7579uwZdJxX4a7R\naGA2m+XHZrO5z0oeACZMmCB/vWXLFjz99NPeHJKG4CqGuu4MFUKwGEoUwbwK90WLFsFoNKKurg6x\nsbE4cuQIDh482GdMU1MTJk2aBAD49NNPMWvWLG8OSb30L4byzlAicvEq3EeOHIn9+/dj5cqVcDqd\n2LhxIzIyMvDiiy8iMzMTubm5eOWVV3Dy5Ek4HA5MmDAB//znP30194jjKoYWFxejpKSExVAiGpIk\nQqSqJkkSC3yDcBVDjUYjLl68yDa5RCEsOjoa27dvx7hx4wJ2zKGyk3eohpjBiqGuQCcichfDPQQM\nVgy12+1wOp0AeLmFiG4ewz1IWlpa5Da5LIYSka8x3APE4XDg4sWL8p2h3d3dAFgMJSL/YLj7UUdH\nh9wml8VQIgokhrsPsRhKRKGC4e4lFkOJKBQx3D3gKoYajUZcvnyZbXKJKOQw3N3gKoa67gxlMZSI\nQh3DfQj9i6Fsk0tE4YTh/gMhBOrr6+W95/2LoWyTS0ThJKLD3Wq1orKyUv7MUBZDiUgpIjLcjUYj\nvvnmGxZDiUixIjLcT548KbfL5eqciJRIFewJEBGR7zHciYgUiOFORKRADHciIgWKyIIqEZEvHTt2\nDO3t7ZAkCZ999hnUajUAIDExEYcOHQrKnLwO9/z8fOzevRsOhwObNm3Cc8891+d1q9WKxx9/HEVF\nRRg7dizeffddTJ061dvDeuTnP/85ampqYDab5b3sADB+/HisXbs2KHMiovDX3t6OmpoaAEB1dXVw\nJ/MDr8LdarVi27ZtOH36NOLi4pCVlYUVK1YgPT1dHvP6669jypQpeP/993Hs2DHs3LkTn3zyidcT\n90RNTQ1OnTo14PnExMTAT4YoRPRuqdG/vYYnr/niPfz9/r5+j1C8V8arcD979ixSUlIQHx8PAFi/\nfj1OnDjRJ9xPnjyJffv2AQBWr16NX/7ylxBCQJIkbw7tU62trfjss8/kx+H2Q8Y5hs57hOMce+v9\neznU1/54LVTew9P3b2trQ6jxKtwtFgsSEhLkxxqNBjqdbsgxKpUKEydORGNjI+Li4ga8X15envy1\nVquFVqv1Znpui4qKwoQJEwBEzg9xpJwn53hz70GeOXTokHxZxt90Ot2AnB2MV+Hu6x+M3uEeSGPH\njsWSJUuCcmwiopvRf+G7Z8+eQcd5Fe4ajQZms1l+bDab+6zkXWNqa2sRGxsLp9OJlpYWTJo0yZvD\nEhGFlPHjxyMxMRGSJCE+Pr7Pbplg8SrcFy1aBKPRiLq6OsTGxuLIkSM4ePBgnzE5OTk4fPgwMjMz\n8cknnyArKwsqVXC217v+Qw+2W4aIyFOu3XbR0dHYvn07xo0bF+QZeRnuI0eOxP79+7Fy5Uo4nU5s\n3LgRGRkZePHFF5GZmYnc3Fw89dRT2LhxI9LS0jBmzBi8++67vpr7TXPtN923b5/cOIyISIkkESIf\nLSRJUsA+5YjhTkT+EIyV+1DZyfYDREQKxHAnIlIghjsRkQIx3ImIFIjhTkSkQBEZ7rzlmoiULiLD\nfd26dZg1axaioqIQExMT7OkQEflcRH5Yx/Tp0zF9+nTY7XZUVVVBr9ejvLwcAGCz2QK2356IyF8i\nMtxdoqKiMHv2bMyePRtOpxMXL16E0WhEUVERnE4ng56IwlZE3qE6HCEE6uvrYTKZYDAYYLVa4XQ6\n4XA4gj01IgphoXSHakSv3Ifi6uwWHx+P7OxsNDU1oaioCHq9HteuXQMA2O32IM+SiGhoDPdhSJKE\n2NhYxMbGQqvVoq2tDcXFxSgsLERLSwtUKhVsNluwp0lE1Acvy3iho6MDJSUlKCwsxKVLl6BWq0Py\nsxSJKDB4WUYhRo8ejczMTGRmZuL69esoKyuDXq/HxYsXERUVBavVGuwpElGEYrj7yI9+9CPMnz8f\n8+fPR09PDyorK6HX61FRUQG1Ws2gJ6KAYrj7wYgRIzB37lzMnTsXDocD1dXVMBgMKCkpAcC99ETk\nfwx3P1Or1Zg5cyZmzpwJp9MJi8UCo9EIk8kEu90Ou93e5yP/iIh8gQXVIBFCoKGhQd5i2dXVBSEE\n99IThTEWVAmSJGHKlCmYMmUK7rvvPjQ3N8tbLK9cuQKAe+mJyHMer9xbW1uxfv16XL58GVOmTMH7\n77+P8ePHDxinVqsxb948AMDUqVNx7NixwScSYSv3G7ly5Yq8xbKxsZF76YnCRCit3D0O9x07dmDG\njBl4+umn8Ze//AXV1dV49dVXB4wbM2aMfFenJxOMdJ2dnSgtLUVhYSHq6uq4l54ohCki3GfMmIFv\nv/0WEydORHNzM5YsWYKKiooB4xjuvmO1WlFeXg69Xo/q6mpusSQKMaEU7h5fc29qasLEiRMBALfd\ndhsaGxsHHdfd3Y3MzEw4nU48//zzeOSRR4Z8z7y8PPlrrVYLrVbr6fQUKSYmBqmpqUhNTYXNZkNV\nVRUMBgPKysqgUqnQ09PD/0ESKZxOp4NOpxt23A1X7tnZ2WhoaBjw/EsvvYTHHnsMV69elZ8bO3Zs\nn8cujY2NiI2NRXV1Ne69917k5+cjKSlp4ES4cveYq12xwWBAcXEx2xUTBUnYrNw///zzIV+bNGkS\nmpubcdttt6GpqQmxsbGDjnM9P23aNKxYsQIFBQWDhjt5TqVSYdq0aZg2bRpyc3NRV1cHk8kEo9HI\ndsVEEcrjj9nLycnB4cOHAQCHDx9GTk7OgDFXrlyRd3m0tLTg1KlTSElJ8fSQ5AZJkqDRaLBy5Urs\n2rULW7ZswZ133onx48cjKioKUVHc/UoUCXyyFXLy5Mk4cuQIxo8fj3PnzuHAgQP4+9//jq+//hpb\nt26FSqWC1WrFzp078eSTTw4+EV6W8bvW1lZ5L31rayu3WBL5WChdluEdqhHq2rVr8l76hoYGbrEk\n8oFQCnf+Gz1CjRkzBosWLcKiRYvkdsWFhYWora1l0BMpAMOdBrQrrqiogF6vR2VlJffSE4Uphjv1\nMWLECCQnJyM5OZntionCGMOdhtS/XbHZbJbbFTscDrYrJgphLKjSTRNC4NKlSygqKoLBYMD169e5\nl54ILKhSmJMkCbfffjtuv/123H///WhqakJxcTH0ej3bFROFCIY7eW3SpEmYNGkS7r77brS3t8tb\nLJuamriXnihIGO7kU+PHj8eSJUuwZMkSuV3xhQsXUF9fzy2WRAHEcCe/GTVqFDIyMpCRkYHu7u4+\n7YqjoqK4xZLIjxjuFBAjR45EWloa0tLS5HbFer0e5eXlkCSJWyyJfIzhTgEXHR2NpKQkJCUlweFw\n4OLFizAajWxXTORD3ApJIUMIgbq6OhiNRhiNRvT09MDhcHAvPYUNboUkGoSrXbGrZXFjYyNMJhMM\nBgM6OjoAcIslkbsY7hSSJElCXFwc4uLicO+996K1tRVFRUXQ6/VsV0zkBl6WobBz9epVeYvl5cuX\nucWSQgYvyxB5YezYsXK74q6uLrldsdlsZtAT/YDhTmHtlltuwYIFC7BgwQK2KybqheFOitG7XbHd\nbpfbFZeWlgJgu2KKLAx3UqSoqCjMmjULs2bNgtPpRG1tLUwmE9sVU8RQefoXP/jgA6SkpECtVqOg\noGDIcfn5+UhLS0NycjJefvllTw9H5DGVSoXExESsWrUKu3fvxqZNm7BkyRKMGTMG0dHRUKvVwZ4i\nkc95vHJPS0vD0aNHsXXr1iHHWK1WbNu2DadPn0ZcXByysrKwYsUKpKene3pYIq/0blecnZ2NpqYm\neYvl1atXAXAvPSmDx+E+Z86cYcecPXsWKSkpiI+PBwCsX78eJ06cYLhTyJg0aRKWL1+O5cuXo729\nHcXFxSgsLERzczP30lNY8+s1d4vFgoSEBPmxRqOBTqcbcnxeXp78tVarhVar9d/kiPoZP348srKy\nkJWVhY6ODpSWlqKwsJDtiimk6HS6G+aoyw3DPTs7Gw0NDQOe37t3L3Jzc4d9c0mShh3TW+9wJwqm\n0aNHY+HChVi4cKHcrriwsBA1NTUMegqq/gvfPXv2DDruhuH++eefezUJjUYDs9ksPzabzX1W8kTh\noH+74srKSuj1elRUVECSJAY9hSSfXJYZau/wokWLYDQaUVdXh9jYWBw5cgQHDx70xSGJgiI6Ohpz\n5szBnDlz4HA4UFNTI7crFkJwLz2FDI+3Qh49ehQJCQk4c+YMVq1ahQceeAAAUF9fj1WrVgH4fsWz\nf/9+rFy5EvPnz8dDDz2EjIwM38ycKMjUajVmzJiBNWvW4LnnnsNjjz2GRYsWYdSoUYiOjoZK5fGv\nF5HX2DiMyMeEELh8+bK8xbKzsxMAt1hGAjYOI1IwSZIwefJkTJ48Gffeey9aWlrkLZZtbW3cYkkB\nwZU7UQBdvXoVJSUlKCwsxOXLlxn0CsOVO1GEGjt2LBYvXozFixfL7YovXLgAi8XCLZbkUwx3oiDp\n3a7YarXK7YqrqqqgUqkY9OQVhjtRCIiJiUFKSgpSUlJgt9tRVVUFg8GAsrIyAGxXTDeP4U4UYqKi\nojB79mzMnj1bbldsNBpRVFTEdsXkNhZUicKEEAL19fUwmUwwGo3o7u6G0+mEw+EI9tToByyoEtFN\nkyQJ8fHxiI+PR3Z2Npqbm2EymaDX63Ht2jUA3EtP/48rdyIFYLvi0MCVOxH51GDtii9cuIBLly5x\ni2WEYrgTKUz/dsVlZWXQ6/VsVxxhGO5ECjZy5EjMmzcP8+bNg81mQ0VFBQwGAyoqKqBSqWC1WoM9\nRfIThjtRhIiOjsbcuXMxd+5cuV2xwWBASUkJ2xUrEMOdKAK52hXPmDEDTqcTFosFRqMRJpMJNpsN\nDoeDe+nDHHfLEJHM1a7YZDLBYDCgs7MTQgjupXcTd8sQUUjq3a74vvvuQ0tLi9yXvq2tDZIkcS99\nmODKnYjccvXqVXkvfWNjI/fSD4IrdyIKO2PHjsUdd9yBO+64A11dXSgtLUVhYSHbFYcoj8P9gw8+\nQF5eHkpKSvDdd98N+dmoiYmJGDt2LNRqNaKjo/Htt996PFkiCg233HIL0tPTkZ6eDqvVivLycuj1\nelRXV0OtVnOLZQjwONzT0tJw9OhRbN269YbjJEmCTqfDrbfe6umhiCiExcTEIDU1FampqXK7Yr1e\nj/LycgBsVxwsHof7nDlz3B7LbyxRZOjfrvjixYtyu2Kn08mgDyC/X3OXJAnZ2dmw2+341a9+haee\nesrfhySiEKBSqTBt2jRMmzYNDz74oNyu2GAwwGq1sl2xn90w3LOzs9HQ0DDg+b179yI3N9etA5w5\ncwaxsbFoamrCj3/8Y8yZMwf333//oGPz8vLkr7VaLbRarVvHIKLQ1r9dcVNTk7zFku2Kb45Op4NO\npxt2nNdbIe+55x688sorQxZUe/vDH/4AAHjhhRcGToRbIYkiUltbm7zFsqWlJay3WCpuK+RQodzV\n1QXg+8p6Z2cn8vPz8eyzz/rikESkEBMmTMDSpUuxdOlSdHR0oKSkBIWFhWxX7CWVp3/x6NGjSEhI\nwJkzZ7Bq1So88MADAID6+nqsWrUKANDQ0ICsrCwsWLAA6enpWL58OVavXu2bmROR4owePRqZmZnY\nsmULnn32WeTk5GD69OlQq9UYMWJEsKcXVniHKhGFvJ6eHlRWVkKv16OyshKSJIXkil5xl2WIiPxp\nxIgRfdoVV1dXy+2KAe6lHwzDnYjCilqtxsyZMzFz5ky2K74BXpYhIkUQQqChoUHeYtnV1RXwdsW8\nLENE5GOSJGHKlCmYMmUK7rvvPjQ3N8tbLNvb2yOuXTFX7kSkeFeuXJG3WPqzXTFX7kREATRu3Di5\nXXFnZ6fcrriurk6xe+kZ7kQUUUaNGoWMjAxkZGT0aVdcVVWlqKBnuBNRxOrdrthms6GqqgoGgwFl\nZWWQJCmst1gy3ImI8P318qSkJCQlJcntig0GA4qLi8OyXTELqkRENyCEQF1dHUwmE4xG4w3bFbOg\nSkQUJiRJgkajgUajwYoVK9DY2Cjvpe/o6Aj4Xnp3MdyJiNwkSRLi4uIQFxeHe+65B62trSguLoZe\nr0dLS0tI7aPnZRkiIh+4du0aamtrkZycDEmSAnbcobKT4U5EFMaGyk6P+7kTEVHoYrgTESkQw52I\nSIEY7kRECsRwJyJSII/DfdeuXUhOTkZycjIefPBBtLS0DDouPz8faWlpSE5Oxssvv+zxRH1Np9MF\newp+o+RzA3h+4Y7nFxgeh3tubi6MRiOKioqQmpqK3//+9wPGWK1WbNu2Dfn5+dDr9fjwww9x/vx5\nrybsK6HyDfAHJZ8bwPMLdzy/wPA43O+55x6oVN//9TvvvBN1dXUDxpw9exYpKSmIj49HVFQU1q9f\njxMnTng+WyIicotPrrm/8cYbWLNmzYDnLRYLEhIS5McajQYWi8UXhyQiohu44R2q2dnZaGhoGPD8\n3r17kZubCwB46aWXUFBQgI8++mjAuPfeew9fffUV9u/fDwD417/+BZ1OhwMHDgycSABv1yUiUpKb\n7gr5+eef3/AN3377bZw4cQJffPHFoK9rNBqYzWb5sdls7rOSH25yRETkGY8vy+Tn52Pfvn04fvw4\nRo4cOeiYRYsWwWg0oq6uDjabDUeOHMEDDzzg8WSJiMg9Hof7jh070NHRgezsbKSnp+PJJ58EANTX\n12PVqlUAgJEjR2L//v1YuXIl5s+fj4ceeggZGRm+mTkREQ1NKNhnn30mUlNTxdy5c8Uf//jHAa93\nd3eLRx55RKSmpoqlS5eKmpqaIMzSc8Od3759+0RycrJISUkRy5YtE1VVVUGYpeeGOz+XDz/8UEiS\nJM6dOxfA2XnPnfN7//33xYIFC0RaWprYsGFDgGfoneHOr7i4WCxevFikpKSIuXPnimPHjgVhlp7Z\nvHmziI2NFampqUOO2bFjh0hOThbp6emioKAggLP7nmLDvbu7WyQmJgqLxSJsNpvIzMwc8B/4T3/6\nk/j1r38thBDi6NGjYvXq1cGYqkfcOb+vvvpKdHd3CyGE2L9/v1i7dm0wpuoRd85PCCGuXr0qli1b\nJrKyssIq3N05vwsXLojFixeLjo4OIYQQLS0twZiqR9w5v0cffVQcOHBACCFEUVGR0Gg0wZiqR776\n6itRUFAwZLh/+OGHYs2aNUIIIQoKCsT8+fMDOT0hhBCKbT/gzh77kydPYuPGjQCA1atX4+uvvw6b\nwq4757ds2TLExMQAGPpehFDl7j0Sv/3tb/H8888jJiYmbL53gHvn99Zbb+Gpp57CqFGjAAC33npr\nMKbqEXfOLyEhAVeuXAEAtLe3Y+rUqcGYqkeWLVuGCRMmDPl672xJT0+H3W4P+DZwxYa7O3vse49R\nqVSYOHEiGhsbAzpPT93sPQQHDx4c9F6EUOXO+RUUFKCurg45OTkAwms7rTvnV1paigsXLiAzMxML\nFy7E8ePHAz1Nj7lzfi+88ALefvttJCQkYNWqVXjttdcCPU2/CYV7fBT7Garh9IvuiZs5v3feeQcF\nBQU4deqUH2fkW8Odn9PpxK5du/D222/Lz4XTyt2d75/T6URNTQ3Onj0Ls9mMpUuX4q677gqLFbw7\n57dr1y6qpyEFAAACGElEQVT84he/wDPPPIMzZ87gscceg8lkCsDsAqP/z2OgM0mxK3d39thrNBrU\n1tYC+P4XqaWlBZMmTQroPD3l7j0E//nPf/DSSy/h+PHjiI6ODuQUvTLc+V27dg0mkwlarRbTpk3D\nmTNnsHr1ahQUFARjujfNne9fQkICcnNzoVarkZiYiOTkZJSVlQV6qh5x5/xOnz6NRx55BACwZMkS\ndHd3h82/nIfT//wtFgs0Gk1gJxHwq/wBcv36dTF16lRhsVhET0+PyMzMHFBw611Q/fjjj0Vubm4w\npuoRd86voKBAzJgxQ1RUVARplp5z5/x602q1YVVQdef8Pv74Y7Fp0yYhhBBNTU3i9ttvF42NjUGY\n7c1z5/xycnLEoUOHhBDfF1Tj4uKE3W4PxnQ9Ul1dfcOCqmsDw7lz58S8efMCOTUhhIJ3ywghxMmT\nJ+VtVnv37hVCCPG73/1OHD9+XAjxfUX/4YcfFqmpqSIrK0tUV1cHcbY3b6jz+/TTT4UQQtx///1i\n8uTJYsGCBWLBggVy9T5cDPf96y3cwl0I985v165dIjk5WSQlJYl//OMfwZqqR4Y7v5KSErFkyRKR\nnJws5s6dK//choOf/vSnYsqUKSI6OlpoNBrx5ptvigMHDsi7f4QQYvv27fJWyGD8bN6wtwwREYUn\nxV5zJyKKZAx3IiIFYrgTESkQw52ISIEY7kRECsRwJyJSoP8DEp+wRefmjHsAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0xeb872d10>"
       ]
      }
     ],
     "prompt_number": 70
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 70
    }
   ],
   "metadata": {}
  }
 ]
}