{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# A Markov Chain Monte Carlo Example\n", "\n", "Here is a short example of using a Markov Chain Monte Carlo (MCMC) algorithm to fit a model to some data. \n", "\n", "**Warning** This is not a very well written MCMC, but hopefully it's quite clear how they work from this. Let me know if you find problems or bugs! (hannah.middleton@unimelb.edu.au)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bayes Theorem\n", "\n", "When we want to do inference, this normally falls into one of two categories\n", " - Model selection\n", " - Parameter estimation\n", "Here we will just look at parameter estimation. Where the model is already choosen and we would like to use the data to figure out what values might be a good fit for the model parameters. \n", "\n", "Bayes theorem\n", "\n", "$$P(\\theta|d,M) = \\frac{P(d|\\theta,M) P(\\theta|M)}{P(d|M)}$$\n", "\n", "where \n", " - $P(\\theta|d,M)$ the posterior distribution of our parameters\n", " - $P(d|\\theta,M)$ the likelihood of the data given some parameter selection\n", " - $P(\\theta|M)$ the prior on our parameters\n", " - $P(d|M)$ the evidence\n", "\n", "Given some data, we would like to work out the probability distributions on the parameters of our model. " ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [], "source": [ "import matplotlib\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "import corner" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example using a quadratic\n", "\n", "Here we will use a simple example and apply Bayes Theorem to it. We make some fake data, which is of the form\n", "$$y = A x^2 + Bx + C$$\n", "where we have three parameters $A$, $B$ and $C$. \n", "\n", "As a first step, we make our fake data and add some Gaussian noise to it. \n" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XlcVPX+x/HXl10RQREFBQUVNxTUcF/S3Cr3tFxazEp/lWbdut02W+5tue3dyqzcrUwtNbUst8pyyQUFTVxxR0RxB2Xn+/vjgKKCDDAzZ4DP8/HwMTBz5pzPmPHmfFeltUYIIYQoipPZBQghhCgbJDCEEEJYRAJDCCGERSQwhBBCWEQCQwghhEUkMIQQQlhEAkMIIYRFJDCEEEJYRAJDCCGERVzMLsCaatSooYODg80uQwghypStW7ee1lr7FXVcuQqM4OBgoqKizC5DCCHKFKXUEUuOkyYpIYQQFpHAEEIIYREJDCGEEBaRwBBCCGERCQwhhBAWkcAQQghhEQkMIYQQFpHAEEIIYREJDCGEmNnX+CNuSgJDCCGERSQwhBBCWEQCA+R2VAghLCCBIYQQZZ2dfumVwBBCCGERCQwhhBAWKVf7YZRYZhpcPg3ZmeDsanY1QgjhkOQOAyAzBc4fhmObzK5ECCEclgQGgIcPoGD/SrMrEUIIhyWBAeDkAu5VYf9qsysRQgiH5RCBoZTyUUotUErtUUrtVkp1UEpVV0qtUkrtz32sZssaMjx84VQsXIi35WWEEKLMcojAAD4GlmutmwARwG7geeBXrXUo8Gvu9zaxPDmEEadHG9/sX2WrywghRJlmemAopaoCXYHpAFrrDK31eWAgMDv3sNnAIFvV0NLjFPt0HU47+aElMIQQokCmBwZQH0gCZiqlopVS05RSnkAtrfUJgNzHmrYqwH/sd7xyV1tWZISTFfc7ZGXY6lJCCFFmOUJguACtgc+11q2ASxSj+UkpNVYpFaWUikpKSipxEUNvCeR0QFdcsy+TuPP3Ep9HCCHsbc75ZhzL8LL5dRwhMOKBeK113iSIBRgBclIpFQCQ+3iqoDdrradorSO11pF+fn4lLkIpxfB77icDF7asmk9Oji7xuYQQwl42HTzDxJO3MuN8uM2vZXpgaK0TgWNKqca5T/UAdgFLgVG5z40Clti6llp+vpyrEUmT5I3M/uuwrS8nhBClcik9i38u2E6Q60X+WcP2E49ND4xcTwBzlFI7gJbAW8DbQC+l1H6gV+73NlezdX9CnY7zzfI/OXz6kj0uKYQQJfLWz7uJP5fK+/6/4emUZfPrOURgaK1jcpuVwrXWg7TW57TWZ7TWPbTWobmPZ+1RiwrtDUA35x38a8EOaZoSQjikP/clMWfTUR7uFELbyifsck2HCAyHUiMUfOrxcM39bD58llkbDptdkRBCXONCaibPLdxBAz9P/tmncdFvsBIJjOspBaG9CDi3hd6NfHh3xR5pmhJCOJTXf9rFyYtpfHBPSzxcne12XQmMgoT2RmVe5u3IZFydnaRpSgjhMFbtOsmCrfE83q0hLYN87HptCYyCBHcBZ3eqJ/zBq/3DpGlKCOEQzl3K4IVFf9PE34sJPUKvvnD2IKTavptXAqMgbpUhuDPsX8mQ1nW4rUlN3l2xh0PSNCWEMNHLS3ZyITWDD+9piZtL7o/vxJ2QnAAZl21+fQmMwoT2hjNxqHOHeGtwi9ymqe3SNCWEMMVPOxL4accJJtwWSrPaVa++sOkLUE5QpZbNa5DAKExoL+Nx/2r8vT14tX8YWw6fk6YpIcobreH8UchIMbuSQiUlp/Py4p1EBHrzWLcGV1+4dAb+/h48a9ple2kJjML4NoDq9a/swidNU0KUU+ePwoWjcGo3pJ4zu5obaK158Ye/uZSRzQf3RODinO/H9rZZkJUGXgF2qUUC42ZCe8PhtZCZilJKmqaEKI9OxBiP2emwdIJxx+FAfog+zqpdJ3m2d2Ma1sy3wGB2JmyeBiG3gpunXWqRwLiZhr2M9D68DuCapqmZ0jQlRPmQEA0o8K4Lu5fC1llmV3TFiQupvLo0lsh61Xioc8i1L+5eanR2t3/MbvVIYNxMcCdwqXSlWQquNk29J01TQpQPCdHGyEjvIKjfHZY/bzRPmUxrzXML/yYrW/P+3RE4O6lrD9j0JVQLgdA+dqtJAuNmXCtBSFcjMHJvU5VS/PeuFrg5O/Hs99vJlqYpIcourSEhBtyqGKs8DP7S+HrBQ5CZamppczcf4899SbxwZxOCa1zX5HR8GxzbBG3HgpP9foxLYBQltBecOwxnDlx5qlZVo2kq6oiMmhKiTDt3GNLOg1tu34BXLSM0Tu2CFS9Zfp6ZfY0/JXXd+4+dvcyby3bRqaEv97Wrd+Pxm3KDrdW9Jb9mCUhgFKVhT+MxX7MUwF3SNCVE2ZcQbTy6Vbn6XGhP6DAeoqbD7h/tXlJOjuaf329HKcW7QyNwur4pKvkk7FwILe8FD2+71iaBUZTqIeAbCnGrrnlamqaEKAdOxICzm9GHkV+PVyGgJSwZD+eP2bWk2X8dZtOhs7zSrxl1fCrdeEDUDMjJNJqj7EwCwxKhvY2RUhnX3knkb5qauf6QScUJIUosIRpqhRkzpfNzcYOhMyAnCxaNhWzbb04EcDAphXeW76F7Yz/ujgy88YCsdCMwQntDjYZ2qSk/CQxLhPaC7Aw4tPaGl+5qXYceTWry3oq9xJ1KNqE4IUSJaA0J26F2q4Jf920AfT+Aoxvgz/dsXk62Vjzz/XbcXZx5e0g4SqkbD4r9AS6dgnaP2ryegkhgWKJeR3D1vKEfA642TVVxd+HxOdtIzcg2oUAhRLGdPQjpF4ymp8JEDIfw4fDnu3B4vU3LmXK2JdFHz/OfgWHUqupx4wFaw8bPoUYjaHCbTWspjASGJVzcof6tsH9VgbNAay4cwkfVFrL/VAqvLt1pQoFCiGLL6/Au7A4jT9/3oVowLBoDl22zhPje9Op8dKYtdzT3Z0BE7YIPOrbZ6HNp93/GEGATSGBYKrSXsd7M6X0FvtzV8xjjujXku6h4Fm6Nt3NxQohiS4gGZ3eo2fTmx7l7Gf0ZKaeMTnArLx2SnpXN0yd64OWUzhuDmhfcFAWw6XNw9zbueEwigWGphnmr197YLJXnqZ6htA2pzsTFO6U/QwhHd2I7+De3bJXX2q2g52uwdxlsmWbVMt5ctpvYdD/+W2sNvlXcCz7ownHYtRRa3w/uVQo+xg4kMCzlEwR+TY1mqUK4ODvx6YhWVHZzlv4MIRxZTo4xw7uo5qj82j9uzMta8ZKxaZEVLN2ewFd/HWFMtRh6ex0u/MAt0wBtylDa/CQwiiO0JxzZAOmF3z3UqurBR8NaSn+GEI7s7AHISC5eYDg5waAvoJKPsXRIKXe4izuVwvMLdxBZrxr/8ttY+IGZqcaCiI3vhGoFzPq2IwmM4gjtbUyYOfjHTQ/r2shP+jOEcGQJuUua32yEVEGq+BlLh5zeZyxSWEKXM7J4fM5WPFyd+XRkK1xVTuEH//29sV+3SUNp85PAKI6g9saaM3GFN0vlkf4MIRxYQjS4eIBfk+K/t0F36PQkbJttzIsoJq01E3/Yyf5TKXw8vCUB3gXM5r56MGz8Amo1h+DOhR83epnxx8YkMIrDxe2mw2uvOVT6M4RwXAnR4N8CnF1K9v7bJkKdW2Dpk3DuSLHeOm/LMRZFH+fJHqF0CfW7+cGH18GpWFOH0uYngVFcob3h4nGL1suX/gwhHFBONiTuKF7/xfWcXWHIdEDDwkcsHmq78/gFXl0aS5fQGjxxW2jRb9j0BVSqDi3uLnmtViSBUVyFrF5bmPz9GYu2SX+GEKY7EwcZKaULDDAWJu33EcRvNvYFL8LFtEzGfbuN6pXd+N+wljduiHS9c4dhzzKIHG3szeMAJDCKy7uO0Z54k+G118vrz3jpB+nPEMJ0eTO8i9vhXZAWQ6HVfXDxGFw+XehhWmue/X47x8+l8tm9rQqfb5Hf5qnGooiRD5e+TiuRwCiJ0F5wbCOkXbDo8Pz9GePmREt/hhBmSogB18rGmkzWcMe74OYJSXtg9b8LXNl2+rpDrIg9yfN3NOGWetWLPmd6Cmz7GpoNNH5JdRASGCXRsJex7PHBNRa/Ja8/Y9+pZF5bGmu72oQQN5cQDf7hJe/wvp6bJ9QKhyq1YN2H8NVASE688vLWI2d5+5c99AmrxcOdQyw75455xsKI7R+zTo1WIoFREkFtjTVdLOzHyJPXnzE/6pj0Zwhhhisd3lZojsrPydnYaG3QF5CwDb7oAof+5ExKOuPmRFOnWiXeuzui8HWirqkxx9iCtXYrCGxj3TpLSQKjJJxdjbHY+1cXeyEy6c8QwkSn90Hm5dJ3eBem5QgY8xtU8kF/NZDfpzzLuctpTL63NVU9LFizCuDgb0ad7R5ziKG0+UlglFRoL0hJhMS/i/U26c8QwkSWLmleGjWbwpjf2VujF0MvzmaN/2eEeRdjx75NX4JnTQgbZLsaS0gCo6TyhtdaMOv7etKfIYRJEmKMzdB8r9ve1MozpdceTeWO+Af4PuAZ/M9tgS+7wNFNRb/xdJzR1N3mYWMfHgcjgVFSXv5Gx1kxhtfmJ/0ZQpggIRoCIow+Bxs5cSGVJ+fFEFrTi76jX0Q9vMpoxp51J2yYdPNm7M1TwMkVIh+yWX2lIYFRGqG94dimEm8QL/0ZQthRdpbRhGzD5qjM7BzGfxtNemY2n993C5XdXIwO9rF/QKPbYeVLMP8+SD1/45tzsiBmDjQfAlVq2qzG0nCYwFBKOSulopVSP+V+H6KU2qSU2q+Umq+UcjO7xhuE9gadA2nnSvT26/szLmeULHiEEBY4vReyUq0/Qiqfd5fvYeuRc7w9JJwGfvk2OqrkA8O+gT5vwb7l8GXXq/0peVJOGjPQ25u/Km1hHCYwgCeB/As0vQN8pLUOBc4BjjPdMU9gJHj4QGrJAgPyrzeVzFPzYsjOse72j0KIXDbu8F6eHMLUtYd4oEM9+he0L7dS0GEcjP7FuJuY3hu2TDeaqLSG5BPGiti27JAvJYcIDKVUINAXmJb7vQJuAxbkHjIbcMAhA87QsIcRGKXY57drIz9e7teMlbtO8vYvRS9qKIQogYRoY3uC6g2sfuojGVV5NvE2IgK9ealvEXuEB7WF/1sLIV1h2dOwaIyxrEhWmrEqrQNziMAA/gf8C8jbRcQXOK+1zmujiQccZ358fnmbKmVcKtVpRncKYVSHekxde4g5m4q3XLIQwgIJMbkd3tb9sXc5I4vHEvrghOaze1vj7mJBh7qnL4z83lgmfedCo7nM2Q2a9rdqbdZmemAopfoBp7TWW/M/XcChBf4Kr5Qaq5SKUkpFJSUl2aTGm2rQw3hMPVvqU73crxndG/vxypJY/thnwmcRorzKzszt8LZu/0VWdg5PfBvNnnRf/hewmsBqlS1/s5MTdH0WHlhibObkXdcYTeXATA8MoBMwQCl1GJiH0RT1P8BHKZW32EsgkFDQm7XWU7TWkVrrSD+/IjYjsYUqfuBWpcQd3/m5ODvx6cjWNKrlxbg529ibKCOnhLCKU7shO92q/QNaa/794y5+3XOKf9dcS/cqRS9xXqCQrlAn0hiq7+BMDwyt9Qta60CtdTAwHPhNa30v8DswNPewUcASk0osWqVqkJ4Mlwpf3thSVdxdmPFgJJXdnHlo1hZOJadZoUAhKrgTuXt4WzEwpq49yNcbj/B/t9bn/moVYwKu6YFxE88BTyul4jD6NKabXE/hKtcwHqNmWOV0Ad6VmD6qDWcvZTBmdpQsHyJEaSVEGwuGVrNwtdgi/LQjgbd+3kO/8ACe61OCfcHLKIcKDK31Gq11v9yvD2qt22qtG2qt79Zap5tdX6HcPI1tFDdONtaxt4IWgd58PLwlO45f4B/zY8iR4bZClFxCNNS2Tof3lsNneXr+dtoEV+P9uyNwKmrnvHLEoQKjTPMONIbXbptttVP2DvNnYt9mLI9N5J3le6x2XiEqlKwMOBlrlR32DiSlMOarKAKrVWLK/ZF4uNpuiRFHZKUdRATuVSG4C2z4FNo8YrWFwx7qFMzh05f48s+DBNfwZETbulY5rxAVxqldkJ1R6v6L0ynpPDhzM85KMWt0W6p55lt8wooLFzoyucOwpi7PGLM1t88t3vtm9jX+FEApxav9m9GtsR8TF+9k7X4ZbitEsVhhhndqRjYPz44iKTmd6Q+2oa5vMYbPliMSGNZUvxvUbg3r/lfiBQkLkrfmVGjNKjz+zTb2nZThtkJY7ESMsYRPteASvT07RzNhXjQ74s/zyfBWtAzysW59ZYgEhjUpZdxlnDsEuxZb9dReHq5Mf7ANHm7OjJ65haRkxx0DIIRDSYg2JuyVYPc6rTWv/7SLVbtO8lr/MHqHOf5cCVuSwLC2xneCXxNY+4GxN68V1fGpxPRRkZy5lM6Yr6JIy5ThtkLcVFY6nNxV4uao6esOMWvDYR7pHMKojsHWra0MksCwNicn6Py00dG2b7nVTx8e6MPHw1uxPf48T38nw22FuKmTscZabyUYIfXL3yd48+fd3NHcnxfvLGJBwQpCAsMWmg8Bn3qw9v1SrWJbmD5h/rx4R1N+/juRd1fstfr5hSg3StjhvfXIWZ6aH0OrIB8+GtayQs21uBkJDFtwdoHOT8HxrXDoT5tc4pEuIYxsV5cv/jjAvM0lXMNGiPIuIdpYusfH8uHoh05f4pHZUQR4ezBtVJsKN9fiZiQwbCViJFTxN/oybEApxX8GhNG1kTHcdn1c6dexEqLcORFj3F1Y2OF9JiWd0TM3o3LnWlT3dLyNPs0kgWErrh7QcTwc+gPio2xyCRdnJz4b2YoGLqcZN+N3YhMu2OQ6QpRJmanGKrUWNkelZWbzyFdRnLiQxtQHIgmu4WnjAsseCQxbumW0Mf577Yc2u4SXhyvfVp/KetfxfDF1MrsSLtrsWkKUKSdjja1QLejwzs7RPDUvhphj5/l4eEtuqVfNDgWWPRIYtuReBdo/BnuXGUP7bCEhBt8LO6ms0nhTf8q/pi4xJzRuMltdCFNY2OGttebNZbtZHpvIxL7NuL15gB2KK5skMGyt7Vhw9YR1H1n/3JdOw/z7wMkV5R+Bp7sLH/ABD039Q+40hEiIMbYe8A4s9BCtNe8s38uM9Yd4sGMwD3e2zvLn5ZUEhq1Vrg5tHoKdC+DsIeudNzsTvn8QLiXBwyvg0T9xHjKVxvoQL6qZ3DttI7tPSGiICqyIGd55dxZf/HGAe9vV5ZV+zexcYNkjgWEPHcaDkwus/9h651z5MhxeC/0/vnrL3agPdH2WATm/MlT9zr3TNrEnUUJDVEAZlyFpT6HNUXnbq05bZ9xZvDGoucy1sIAEhj14+UOr+yBmDlwscGvy4omZC5s+h3aPQcTwa1/r9gLU78aLTKe5OszIqRIaogI6uRN0doGBkZOjeXnJzitLfrzavxmqBOtMVUQSGPbScQLkZMNfn5XuPMe3wY9PGntv9H79xtednGHIdJRnDaZX/hRfp8sSGqLiyevwvm6EVE6O5sUf/uabjUd59NYGvNS3qYRFMUhg2Ev1EGgx1Nj3+/LZkp0j5ZTRyV2lJtw9C5xdCz7OswbcPRvXlAQW1/kaNyfNyKmb2Jsoy6KLCiIhBjxrQtXaV57KztE8u2AH87Yc44nbGvLc7Y0lLIpJAsOeOj8NmZdh0xfFf29eJ/flszB8jhEKNxPUBvq8hefhVfx8yzZcnRUjp26U0BAVQ0L0NTO8s7JzeOa7GBZui+cfPRvxTG8Ji5KQwLCnmk2gST8jMNKL+YN7xYtwZD0M+BQCIix7T9sx0Hwo1Te9ww+3Z+LsZISGbMAkyrWMS3B6rzFCCsjMzuGp+TEsjkng2T6NebJnqMkFll0SGPbW5WlIu2A0TVkq+hvYPMUYbRV+t+XvU8oYRVWjEbVXj+f7kXVxdlKMmHJdaMikO1GeJP4NOgdqtyIjK4cJc6P5accJXryzCeO6NzS7ujJNAsPe6twC9bvDhkmQmVb08fFR8NM/IORW6Pnv4l/PvQrc8zVkpVHv13HMe7j1lTuN/XKnIcqj3A7v9JoteHzONn7ZmcjL/ZoxtmsDkwsr+yQwzNDlGbh0CmK+uflxySeNTm4v/9xObpeSXc+vEQycBPGbqR/9DnPHtsdJKUZIaIjyKCEGXcWfRxcnsHr3Sf4zMExmcFuJBIY1jF5m/LFUcGcIbGtM5MvOLPiYrAz47gFIPQ/DvzVmjJdG2GBoPw42fUGDxOXMHdsepRQjpm5if7ostCbKj5yEbURnBfP73iTeGtyCBzoEm11SuSGBYQaljLuM80dh58KCj1n+HBzbCIM+A/8W1rlur39DUHtYOoEGOp65Y9qjFIw4NpA4CQ1RDlxOPgen97MmJZB3h4Yzsp3lGyeZqri/dJpEAsMsjfpArebG0ufXb+O6dbbRKd7pSWO7V2txdjWattwqw3f309BbG6GBZvixgew8LvtpiLLrUnoWb89agBOa9p1u457IILNLKnckMMyiFHT+hzH8L/XM1eePbYaf/wkNboMer1r/ulUDYOhMOBMHS5+goZ8nc4OW4KpyuPuLv1gRm2j5uS6dhu3zYeEYSNhmzDERwgTJaZmMmrEZt5PbAejYuafJFZVPEhhmChsM1evDhXjjLuPiCZh/vzE7dch0Y5kPWwjpYoRR7A+w6Qsaup9nSb0FNPL34v++3srkNXHo6+96wFja5Ogm+O1NmNIN3msIP4yFA79BVprxOYSws9Mp6dz/9lfEHDnN6JBzULUOeNUyu6xyqYTDboRVODlDp6fgxwmQetbo5E5PhvsXlb6TuyidnoT4LbByIvg1o6YHzH+wPf/8fjvvLt9L3KkU/ntXC9xTkyDuV4hbBQd+h7TzoJyMTvvuL0FoT/CPgI+aQXKi8cfL37a1C5Fr5/ELjP0qijPpNfis9krqXN5r0Q57omQkMMwWMcJogkraA2i4ezbUCrP9dZWCQZONO4XTeyCgJR6uznx6T3Nudd9HUvQHJOzbSUjWQeP4Kv7GLPXQnlC/G1S6rpPcqzYkn4DNU6HHy7avX1R4P25P4NkF26lW2Y0FdX+ghdsJOLYfwoeZXVq5JYFhNhc34xb63CFjramwQfa7toe3Manvyy5wajfMvw918A/uTr9IjqsLWzJDWeF6H7cPeoDgsLaFbkQDgGslqOQLUdON2exunvb7HKJCycnRvL9yL5PXHCCyXjU+v+8W/BZ8CGkpxgFFbMkqSq7IPgyl1GqllIWLF4kS8aoNtVrAbRPtf23/5lC9AWQkG0unhw2GYd/g9NxBPB5Zzgw1mH7fn+f3vUlFn6tqbUg9B9vn2r5uUSFdTMtkzFdRTF5zgBFtg/h2THv8vNyNFzPyAkOapGzFkk7vfwEfKaVmKqVkd3RbUMr4bd9WndxFqVILAtvBP2JhwCfQtD94eBMR5MOS8Z2o51uZh2dvYdragwV3hudxr2osffLXZMjJsV/9okI4dPoSgz9bz5p9SfxnYBhvDW6Bm0u+H2HpKeAdVPRKzqLEigwMrfU2rfVtwE/AcqXUq0qpSrYvTdiVs2uBTU4B3pX4/tEO9G7mzxvLdvPiD3+TkVVIGCgFHcbB2QOwb7mNCxYVyR/7khg4aR1nL2XwzcPteKBD8I3Lk2ekyN2FjVk0rFYZ/2X2Ap8DTwD7lVL327Iw4Tgqu7kw+d7WjOvegLmbj/HAjE2cu5RR8MFNBxq/5ZV2Z0EhMPbenvrnQUbP3Extn0osHd+ZDg18bzwwJ8sY2i0jpGzKkj6MdcBx4COgDvAg0A1oq5SaYsvihONwclI826cJH94TwbYj5xk8eT1xp1JuPNDZBdo9CkfWXd0mU4gSSMvM5pnvtvPmz7vp3cyfhY91JKh65YIPTpcOb3uw5A7jUaCO1rqX1vplrfVPWus4rfUTQJfSFqCUClJK/a6U2q2UilVKPZn7fHWl1Cql1P7cR1nsyAHc1TqQuWPbkZyWxeDJ61m7v4DO8Nb3g5uX3GWIEku8kMawL/9iUfRxnu7ViMn3tsbT/SaDOjMkMOzBkj6Mnbrwnk5r7LqTBTyjtW4KtAfGKaWaAc8Dv2qtQ4Ffc78XDuCWetVZPK4Ttb0r8eDMLXz91+FrD/DwhtYPGDPJZfa3KKZtR8/Rf9I64k6l8OX9tzChRyhOTkVsp5qRAs7utp/wWsGVamkQrfXB0hagtT6htd6W+3UysBuj6WsgMDv3sNmAHScoiKIEVa/Mwsc7cmsjP15eEsurJzuTqfP9c2r/qLHr2aYvzStSlDnfRx1j+Jcb8XB1YtHjnegTZuGqARkpxmZhwqYcai0ppVQw0ArYBNTSWp8AI1SAmoW8Z6xSKkopFZWUZMFcAWE1VdxdmPpAJGO6hDD7fDiDjgxh94mLxos+daHZQGPl3eLuXy4qnKzsHP79YyzPLthBZHA1lo7rTGN/L8vevM3YURL3qrYtUjhOYCilqgALgae01hctfZ/WeorWOlJrHenn52e7AkWBnJ0UL/Vtxhe1f+FklicDJq1j0m/7ycrOMfYgT78A0XPMLtP2ZF/0EjtxIZVRMzczc/1hHuwYzFcPtaWap5tlb972NSx9Ajx8wEumidmaQwSGUsoVIyzmaK0X5T59Mm+iYO7jKbPqE0W73esQK4Pn0SfMn/dX7uOuzzew37WxsWHTxsnGSrdC5JOTo5mz6Qi9PvyTrUfO8e6QcF4bEIaLs4U/lvLCosFtULOZsSimsCnT/4Zz53hMB3ZrrT/M99JSYFTu16OAJfaurcKwxm5fo5dRfcxCJo1szWcjWxN/LpW+n6xjhfcQOH8E9vxknVpFuXDo9CVGTN3ISz/sJDzQm5VP3co9bYqx4VH0N1fDYvi3EhZ24giLD3YC7gf+VkrF5D73IvA28J1S6mHgKHC3SfWJYuobHkC7+tWZ+MNOHovKYkNlf7z//IRKzQaaXZowWVZ2DjPWH+KDlftwc3bi7btaMKxN0I2ztm8m+htYMv5qWLh62K5gcQ3TA0NrvQ4o7F9LD3vWIqynRhV3Pr+vNUu3+zNr8Z08nziDJT8tpt+dA3EuaoiksL+8/hcfx6OWAAAZ10lEQVQb7iu9+8RFnlu4gx3xF+jZtBZvDGqOv7dH8a5/JSy6S1iYwPTAEOWXUoqBLevQMXAilyZ/h/OmyQyPr8l7QyMIriHLn1cU6VnZfPb7ASb/Hod3JVcmjWxF3xYBxburAGPwhISFqaThT9icXw1fKnd4hL7OW7iYeIDbP/6TWesPkZNzk5VvRbmw7eg5+n2yjk9+3U//iNqsfvpW+oXXLmFYjDM27xr+rbH/irA7CQxhF6rtWJSTE4ta7aB9fV9e+3EXI6Zu5OiZy2aXJmzgckYWr/+0iyGfbyAlPYuZD7bho2EtLR8um1/+sBgxV8LCRBIYwj6860DYXXjGfsvM4Y14d0g4uxIucvvHf/L1xiNyt1GOrI87TZ///cn0dYe4r109Vv6jK92bFDjvtmgx30pYOBAJDGE/HcZBRgpq21fc0yaI5f/oyi31qvHy4p3cP2MTexNlRnhZdiE1k+cX7uDeaZtwcXJi/tj2vD6oOV4eriU7Ycy3sPhxCQsHIoEh7Kd2SwjuYqwvlZ1JHZ9KfPVQW94a3ILtxy5w+8d/8sTcaOJOSXCUJVprlu9MpNeHf/D91ngevbUBvzzZhXb1C9i3wlJXwuJWCQsHIqOkhH11GAdzh8OuJdBiKEopRraryx3N/Zm27iAz1x/mpx0JDIyozYQeodT3kwXlHJXWmtW7TzHpt/1sj79A04CqTB/VhhaB3qU78TVhMU/CwoFIYAj7Cu0Dvg3hr0nQfMiVbWGrebrxbJ8mPNQphClrD/LVhiMs3Z7A4FaBTOjRkHq+MgzXUWTnaH7ZeYJJv8WxJzGZoOqVeGtwC+6ODMTV0mU9ChMz92pYDC/GnYUN54+IqyQwhH05OUH7x2HZ03D0L6jX0Xg+d+KW7+hlvHBHUx7pXJ8v/zjA1xuPsDjmOENbBzL+toaF77gmbC4zO4elMQl8tiaOg0mXaODnyUfDIugfXhuXr/rDLkr3gzvlJCx+DEK6GmHhJv+tHY0EhrC/iBHw2xvGjnx5gXEdPy93JvZrxtiu9Zm85gDfbj7Kwm3x3NMmiHHdG1LHR5op7CU9K5sFW+P5fM0B4s+l0jSgKpPvbU2fMP/Sz9rPyjCCIvkEnD0AIbnNUBIWDkkCQ9ifW2Vo8zD8+T6cOQC+DQo9tGZVD14bEMajtzZg8po45m0+xoKoeIblBseVpSWE1aVmZDN381Gm/HmQxItpRAT58Fr/MHo0rVn0xLuMS5CcmBsGN3lMPXv1PR7eEhYOTgJDmKPNGFj/MWz8HPq+X+Th/t4e/Gdgc/7v1gZ89nscczcfZX7UMUa2rcvj3RpQs6oEh7Ukp2XyzcajTFt7kDOXMmgbUp33746gU0PfgoMiK924O8i4BJ+0NsIgb4/t/JxcoUot8KoF1etD3Q7g5W88t+kLY08LCQuHJoEhzOFVC1rcAzFzoPuLFr+tjo/RwfrYrQ2Y9FscX288wtzNR7mvfT1GtqtLAxlVVWIXst2ZuXofM9cf5kJqJl0b+TG+e0PahhSxT/aa/xpNSu5e4N8CvHrlBoN/vkd/Y7/twu5Mdnxn/Q8krE4CQ5inw+MQ8w1snVnstwZVr8w7/evzdMMTRK9fQdqm3bywvgfJ/u3oHxFA//Da0kFuAa01u05cZOmpDsy5EEZK3H56Nq3FE7c1JCLIp+gTHNts3Cl61oIaoXDPbNsXLUwjgSHMUysM6neHTVOMJoqiNsE5fwyObbr6J3EntXQ2t6PIcXFmkPMGfknry7PLB/Pu8r1EBPnQPzyAvuEBBHhLJ3menBzNtqPnWBGbyPLYRI6dTcWJCO7wOsi40aNpVtvCvbEzLsMPj0LVOsYfUe5JYDiCijyGvON4+GYIuFeFKvnWG8rOgpN/w9F8AXHxuPGaqycE3gJdnoGgdhAYidO3w+H8Ee5I+YVeNaJZWf95Jh+vyhvLdvPGst20Da5Ov4gA7mgegJ+Xu3U/Q042pF0EJ2dITwF3E5rFithPIjM7h78OnGFFbCIrd50kKTkdV2dF54Y1GNetIT23P00Nl1So/aTl11z9mtF3MepHWPNO6T+DcHgSGMJcDXqAXxPj7sHZBX593QiH41shM3cl26qBULe9EQ5B7aBWc+PY/B5ebjzGR+GyZBx37pjAneHDODLoZZbuS+PHHQm8siSW15bG0qGBL/3Da3N7c398Kpdg9dQ8OdmwazH88S4k7TGe+28d8KwJ1YKheghUC8l9DDa+rlKz8HZ8K0vNyObP/Ums2JnI6t0nuZiWRWU3Z7o3rknvsFp0b1KTqnnrPMWmFu/kB/+AzV9Cu0eNeRMhXa3/AYTDUVqXn1VCIyMjdVRUlNlliOLa9pWxPzOAcgb/5hDUHurmBoR3YPHOl5UOaz8w/nj4wJ3vQdhg9p5M4acdCfy4PYHDZy7j4qToElqDfuG16RVW6+oPz6JcHxR+TSAnxwiCiGFw9hCcO2z8uRAP5Pt/zLXy1fDIHyRr/gsuHiW/28y9w7gwfDG/7znF8p2JrNl3irTMHHwqu9KzaS36hPnTJbQGHq7Ohb7fouunXYDPO4GzGzy6zjojm+yw458onFJqq9Y6sqjj5A5DmC9ipNFx6uwBD68ofZOOi7sx8qrpAGNp7AWjYedCGt/5Po17N+bpXo2ITbjIjzsS+Gn7CZ75fjuuixSNannR2N+LJv5eNPGvShN/L/y83K8OJS0oKIbOhGaDYHZ/45guz1xbS1Y6nD9qhMfZQ3DukPF49iAc+A2y8v1m79uoWB9Ta82ZSxkcTLrEnnPNWZ0SzF9vrCIzW1Orqjv3RAbRJ8yftiHVS79kR34rXjSaBx9aKcNgKxgJDGE+Zxdj2CVYt/3fvzk88its/Ax+fws+awd93kC1up/mdbxpXseb529vQvSx86yITWRXwkXW7T/Nom3Hr5yiWmVXmtSqzBD3KHqdno13ykFyajTGKS8onIr4QeziboweqhF642taG5PXzh2GOUPg/BFj5rPLtc1kGVk5HDlziQNJlzh4OoUDp/IeU7iYlpV7VFeCXc/zUOcQ+oT50zLQBydb7J2+d7mxr3bnpyGojfXPLxyaBIYo35xdoNOT0KSf0ey19AnYuRD6fwzVglFK0bpuNVr/ep+x2P9Lyzh7KYM9iRfZd+I8bnuX0iVhBkHZx9iXU4cXsybwy/G2BP1Shcbbthl3IwFVaZRejcpOmThdSMPJCZyVwtlJ4eSkrn595dHY7xyloGoAVA1AD52J+vYeDv06lY3V+nMwKcUIiKQUjp1LJTvfBlO1qrpTv0YVBrSsTf0aVajv50mDPyYQ6JKMuuNe2/1dXj4LP04w+pC6PW+76wiHJYEhKgbfBjDqJ2POx6pXYXIH6PEKtB1rjG7Kp3olZzpeXkPHmHfh9F7wa0JO15m4B/RiwMlLNE5MZk/iRfYkJrN690mMn+UjjDf/91eLynFSXBMiWmu+VQ3x2/Ahr6QH4uTiRkgNT8Jqe9M/ojYN/IxgCKnhWfCGRBvssIfIsmeM0LhvoXHnJCocCQxRcTg5GWtYNeoDPz4Fy5+HnYtg4CTjda3h7wVGH0VuUOT1UTg5OVEPqFfDiz5h/ldOmZaZTdypFOIWvEpajgvZHZ8gJ0eTo41lwHO0JjtHk601OTma7Byufn3lOY0GTmU/Tavox9l0RwI+XR61TZNSSe1cCLGL4LaJxmxua5PO7jJBAkNUPN6BcO/38Pf38Mtz8EVnqOwLaedh4fprO7OL6KPwcHU2+kOq7jeeaFev5HXppnB6NtW3fgodR4OTg6yPlZxo3F3UuQU6/cPsaoSJZItWUTEpBeH3wLjN0KSvsRYSygiKx/6C5ncV3aFti5q6545A2vaVfa9dGK1h6QTITIVBX9w4/0VUKBIYomKr4gd3z4I6bSCglTlBkV/IrVC3I6z70Pghbbbob2D/CujxKvgVb9ivKH8kMIQAoxPXTjOwbyrvLiP5BGydZW4t54/C8hegXmdjRreo8OT+UjgG6fS8KqQLBHeBdR9B61HmTI7LyTH21kbDoM/MvesSDkP+FQjhiLq/aGxEFDXDnOtvmQqH10Kft4zlS4RA7jCEMDjaHU69jlC/m3GXETka3Dztd+3TccZclYa9oPUD9ruucHhyhyGEo+r2Ilw+DVum2e+a2Vmw+FGjT2fAp47RryMchgSGEI6qbjtj+ff1Hxv7bNjDho8hfgv0/cBYtkSIfCQwhHBk3V+Ey2dg8xTbXyvjEvz+X2g2EJoPsf31RJkjgSGEIwuMhNDesOETY1c/W9E5cHofVPKBvh9JU5QokASGEI6u2wuQeg42fWm7a5w/BpmXoP8n4Olru+uIMk0CQwhHV6c1NLoD/vrU2O3O2qLnwMVjxtayTe60/vlFueHwgaGUul0ptVcpFaeUkkX4RcXU/QUjLDZ+bt3zrv8YljxubGVbvYF1zy3KHYcODKWUM/AZcAfQDBihlGpmblVCmCAgwtgE6q/JRvNUaWkNKyfCqlcg7C6o2eyGfUGEuJ5DBwbQFojTWh/UWmcA84CBJtckxI1GL7P95L9uL0D6BSM0SiM7y9jrfMOn0OYRGDINlKP/KBCOwNH/ldQBjuX7Pj73OSEqHv/mxpDXjZ8bO9+VRGYqzL8PYuYYAXTn+3JnISzm6IFR0Ng+fc0BSo1VSkUppaKSkpLsVJYQJrn1echIgb8mFf+9qefh67tg33IjKLo9L8NnRbE4emDEA0H5vg8EEvIfoLWeorWO1FpH+vn52bU4IeyuVjMIG2wMsb10xvL3JSfCrL7GLO6h06HtGNvVKMotRw+MLUCoUipEKeUGDAeWmlyTEObq9rwxK3vDJ5Ydf/YgzOgDZw/ByPkyi1uUmEMHhtY6CxgPrAB2A99prWPNrUoIk/k1hhZDjeVCUopohj2xA6b3MWaJj/oRGvawT42iXHLowADQWv+stW6ktW6gtX7T7HqEcAi3PgdZacZigYU5vM5ohnJ2hYeWQ+At9qtPlEsOHxhCiALUCIXwYbB5GiSfvPH1PcuMDm4vf3h4pXFXIkQpSWAIUVZ1fRayM2D9/659PvobY+isf3N4aAV4B5pTnyh3JDCEKKt8G0DECGMb16x047l1/zMm5dXvBg8shcrVzaxQlDMSGEKUZV3/CTlZcDEezh2C1a8aS32MmA/uVcyuTpQzEhhClGXVQ6DlSEg+ARePQ5sxMGQ6uLiZXZkohyQwhCjruj4Lzm7gXRfufA+c5H9rYRsuZhcghCgln7pQp42xzIcs9SFsSH4VEaI8kKAQdiCBIYQQwiISGEIIISwigSGEEMIiEhhCCCEsIoEhhBDCIuV+WG1mZibx8fGkpaWZXUqZ4uHhQWBgIK6urmaXIoRwEOU+MOLj4/Hy8iI4OBglQw8torXmzJkzxMfHExISYnY5QggHUe6bpNLS0vD19ZWwKAalFL6+vnJXJoS4RrkPDEDCogTk70wIcb0KERhm69ixY5HHPPLII+zatavY546JieHnn3++8v3SpUt5++23i30eIYQoSrnvw3AEGzZsKPKYadOmlejcMTExREVFceeddwIwYMAABgwYUKJzCSHEzcgdhh1UqWLsS7BmzRq6devG0KFDadKkCffeey9aawC6detGVFQUACtXrqRDhw60bt2au+++m5SUFAC2bNlCx44diYiIoG3btly4cIFXXnmF+fPn07JlS+bPn8+sWbMYP348AEeOHKFHjx6Eh4fTo0cPjh49CsCDDz7IhAkT6NixI/Xr12fBggX2/isRQpRBFeoO498/xrIr4aJVz9msdlVe7R9m8fHR0dHExsZSu3ZtOnXqxPr16+ncufOV10+fPs0bb7zB6tWr8fT05J133uHDDz/k+eefZ9iwYcyfP582bdpw8eJFKleuzH/+8x+ioqKYNGkSALNmzbpyrvHjx/PAAw8watQoZsyYwYQJE1i8eDEAJ06cYN26dezZs4cBAwYwdOhQ6/yFCCHKrQoVGI6gbdu2BAYaeyy3bNmSw4cPXxMYGzduZNeuXXTq1AmAjIwMOnTowN69ewkICKBNmzYAVK1atchr/fXXXyxatAiA+++/n3/9619XXhs0aBBOTk40a9aMkydPWu3zCSHKrwoVGMW5E7AVd3f3K187OzuTlZV1zetaa3r16sXcuXOveX7Hjh2lHrmU//3568hrFhNCiJuRPgwH0759e9avX09cXBwAly9fZt++fTRp0oSEhAS2bNkCQHJyMllZWXh5eZGcnFzguTp27Mi8efMAmDNnzjV3MkIIUVwSGA5EKYWfnx+zZs1ixIgRhIeH0759e/bs2YObmxvz58/niSeeICIigl69epGWlkb37t3ZtWvXlU7v/D755BNmzpxJeHg4X3/9NR9//LFJn0wIUR6o8tQcERkZqfNGGuXZvXs3TZs2Nakiy7Vo0YKlS5c61FIcZeXvTgAz+xqPo5eZ835RpimltmqtI4s6Tu4wHECvXr1o0aKFQ4WFEEJcr0J1ejuqVatWmV2CqOjkzkJYQO4whBBCWEQCQwghhEUkMIQQQlhEAqMgM/teHTUihBACkMAwxWuvvcb7779f6OuLFy8u0VLnQghhSxIYDkgCQwjhiCQw7OTNN9+kcePG9OzZk7179wIwdepU2rRpQ0REBEOGDOHy5cts2LCBpUuX8uyzz9KyZUsOHDhQ4HFCCGFvFWsexi/PQ+LfRR+XuMN4tKQfw78F3HHzHe62bt3KvHnziI6OJisri9atW3PLLbdw1113MWbMGAAmTpzI9OnTeeKJJxgwYAD9+vW7suS4j49PgccJIYQ9VazAMMnatWsZPHgwlStXBriyI97OnTuZOHEi58+fJyUlhT59+hT4fkuPE0IIWzI1MJRS7wH9gQzgADBaa30+97UXgIeBbGCC1npFqS9YxJ3AFTZYV6egpckffPBBFi9eTEREBLNmzWLNmjUFvtfS44QQwpbM7sNYBTTXWocD+4AXAJRSzYDhQBhwOzBZKeVsWpWl1LVrV3744QdSU1NJTk7mxx9/BIwlygMCAsjMzGTOnDlXjr9+yfLCjhNCCHsyNTC01iu11nk7CG0EAnO/HgjM01qna60PAXFAWzNqtIbWrVszbNgwWrZsyZAhQ+jSpQsAr7/+Ou3ataNXr140adLkyvHDhw/nvffeo1WrVhw4cKDQ44QQwp4cZnlzpdSPwHyt9TdKqUnARq31N7mvTQd+0VovKOB9Y4GxAHXr1r3lyJEj17xeoiW6ZalnQJY3L1Pk36woBUuXN7d5H4ZSajXgX8BLL2mtl+Qe8xKQBeS1txS0F2mByaa1ngJMAWM/jFIXDPI/nSh75N+ssAObB4bWuufNXldKjQL6AT301dudeCAo32GBQIJtKhRCCGEJU/swlFK3A88BA7TW+WejLQWGK6XclVIhQCiw2YwahRBCGMyehzEJcAdW5Q473ai1flRrHauU+g7YhdFUNU5rnV3Si2itCxzWKgrnKH1bQgjHYWpgaK0b3uS1N4E3S3sNDw8Pzpw5g6+vr4SGhbTWnDlzBg8PD7NLEUI4ELPvMGwuMDCQ+Ph4kpKSzC6lTPHw8CAwMLDoA4UQFUa5DwxXV1dCQkLMLkMIIco8s2d6CyGEKCMkMIQQQlhEAkMIIYRFHGZpEGtQSiUBR4o80DHVAE6bXYTJKvrfgXx++fxmff56Wmu/og4qV4FRlimloixZy6U8q+h/B/L55fM7+ueXJikhhBAWkcAQQghhEQkMxzHF7AIcQEX/O5DPX7E5/OeXPgwhhBAWkTsMIYQQFpHAcEBKqX8qpbRSqobZtdiTUuo9pdQepdQOpdQPSikfs2uyB6XU7UqpvUqpOKXU82bXY09KqSCl1O9Kqd1KqVil1JNm12QGpZSzUipaKfWT2bXcjASGg1FKBQG9gKNm12KCVUBzrXU4sA94weR6bE4p5Qx8BtwBNANGKKWamVuVXWUBz2itmwLtgXEV7PPneRLYbXYRRZHAcDwfAf+ikC1pyzOt9UqtdVbutxsxdlos79oCcVrrg1rrDGAeMNDkmuxGa31Ca70t9+tkjB+adcytyr6UUoFAX2Ca2bUURQLDgSilBgDHtdbbza7FATwE/GJ2EXZQBziW7/t4KtgPzDxKqWCgFbDJ3Ers7n8YvyTmmF1IUcr98uaORim1GvAv4KWXgBeB3vatyL5u9vm11ktyj3kJo6lijj1rM0lBu3pVuLtLpVQVYCHwlNb6otn12ItSqh9wSmu9VSnVzex6iiKBYWda654FPa+UagGEANtzdwYMBLYppdpqrRPtWKJNFfb58yilRgH9gB66Yoz5jgeC8n0fCCSYVIsplFKuGGExR2u9yOx67KwTMEApdSfgAVRVSn2jtb7P5LoKJPMwHJRS6jAQqbWuMIuxKaVuBz4EbtVaV4gtEpVSLhgd/D2A48AWYKTWOtbUwuxEGb8dzQbOaq2fMrseM+XeYfxTa93P7FoKI30YwpFMAryAVUqpGKXUF2YXZGu5nfzjgRUYHb7fVZSwyNUJuB+4Lfe/eUzub9vCAckdhhBCCIvIHYYQQgiLSGAIIYSwiASGEEIIi0hgCCGEsIgEhhBCCItIYAghhLCIBIYQQgiLSGAIYUO5ez30yv36DaXUJ2bXJERJyVpSQtjWq8B/lFI1MVZiHWByPUKUmMz0FsLGlFJ/AFWAbrl7PghRJkmTlBA2lLsKcQCQLmEhyjoJDCFsRCkVgLGnx0DgklKqj8klCVEqEhhC2IBSqjKwCGO/6t3A68BrphYlRClJH4YQQgiLyB2GEEIIi0hgCCGEsIgEhhBCCItIYAghhLCIBIYQQgiLSGAIIYSwiASGEEIIi0hgCCGEsMj/AxKPZHr7JBVJAAAAAElFTkSuQmCC\n", "text/plain": [ "