{
 "metadata": {
  "name": "",
  "signature": "sha256:4573d6b3f33faaaac0570e250160a84b36eafe64824a39b91ed926f2ad33ff45"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Afternoon Breakout: Do-It-Yourself MCMC\n",
      "\n",
      "The standard Bayesian approach to model fitting involves sampling the posterior, usually via a variant of Markov Chain Monte Carlo (MCMC). Though there are many very sophisticated MCMC samplers out there, the most simple algorithm (Metropolis-Hastings) is rather straightforward to code.\n",
      "\n",
      "Here we'll walk through creating our own Metropolis-Hastings sampler from scratch.\n",
      "\n",
      "(The following is based on material put together by Adrian Price-Whelan)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Preliminaries\n",
      "\n",
      "As usual, we start with some imports:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "from scipy import stats\n",
      "\n",
      "# use seaborn plotting defaults\n",
      "# If this causes an error, you can comment it out.\n",
      "import seaborn as sns\n",
      "sns.set()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And load some data:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from fig_code import linear_data_sample\n",
      "x, y, dy = linear_data_sample()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Exercise\n",
      "\n",
      "Walk through all the following steps, filling-in the code along the way."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First plot the data to see what we're looking at (Use a ``plt.errorbar()`` plot with the provided data)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.errorbar(x, y, dy, fmt='o');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFVCAYAAADc5IdQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHJtJREFUeJzt3X+QnVd93/H33bW0ws61uVO2VimMNC3DGf2BEmFaB8nj\nHwMBTNJCUdtohpJAoIkrJTIZSU5QjCmBpKRaua1AgtZA5EyY7gT/oCUdXM8E12q1rTHpEsNIPbbD\nSBlaXLaaxbsWWv3Yvf1j711d7d69v/Y+99wf79cMw3Of++P57tF6P/c8z3nOyRWLRSRJUjpDqQuQ\nJGnQGcaSJCVmGEuSlJhhLElSYoaxJEmJGcaSJCV2Xa0nQwjDwEPAG4EicA+wHvhT4PnSyz4fY/yT\nLIuUJKmf1Qxj4BeAhRjjbSGEO4DfA74OHI4xPph5dZIkDYBcvUk/QgjDMcb5EMIvA3cBPwECi0H+\nAvDRGOMrmVcqSVKfqnvNuBTEx4F/A3wF+BawP8Z4B/B94BOZVihJUp+rd5oagBjjB0MINwPPANtj\njP+n9NTXgCP13l8sFou5XK71KiVJ6i1NhV69AVwfAF4XY/wXwAVgAXgshPAbMcZngbcB365bUS7H\n1NRsM3WpSaOjeds4Y7ZxZ9jO2bONszc6mm/q9fV6xo8Ax0MITwPrgHuBvwKOhhAuAz8EfrWFOiVJ\nUknNMI4xXgB+scpTt2VTjiRJg8dJPyRJSswwliQpMcNYkqTEDGNJkhIzjCVJSswwliQpMcNYkqTE\nDGNJkhIzjCVJSswwliQpMcNYkqTEDGNJkhIzjCVJSswwliQpMcNYkqTEDGNJkhIzjCVJSswwliQp\nMcNYkqTEDGNJkhIzjCVJA+PAsQkOHJtIXcYKhrEkSYkZxpIkJWYYS5KUmGEsSVJihrEkSYkZxpIk\nJWYYS5KUmGEsSVJihrEkSYkZxpIkJXZdrSdDCMPAQ8AbgSJwD3AROA4sAN8D9sQYi9mWKUnS2oyN\nT3JuZm5pe/+ubYkruqpez/gXgIUY423A/cDvA4eBgzHG24Ec8J5sS5QkaW3Gxic5dWZ66fGpM9Ps\nO3qSsy/Nrnhtivmra4ZxjPE/AL9WergZmAZuiTGeKO37BvD2zKqTJKkNTlcEcdn07EWOPPpcgmpW\nqnmaGiDGOB9COA68F/hHwM9VPP0KcFMjBxodzbdSn5pgG2fPNu4M2zl7A9fGORYvti4zNJRb0RbD\nwzmgs21UN4wBYowfDCHcDHwL2FDxVB74cSOfMTW18lSA2md0NG8bZ8w27gzbOXuD2MZbNhWuOU0N\nUMiP8Ov/4E0r2mJ+fjG1q7VR+fT1od3bax6v2SCveZo6hPCBEMLHSg8vAPPAt0MId5T23Q2cqPpm\nSZK6xP5d2yjkR5YeF/IjHN6zg00bu+MMQb0BXI8APxNCeBp4ArgX+HXgkyGECRZ71o9kW6IkSWu3\nd+fWqtvdoOZp6hjjBeAXqzx1ZybVSJKUkU0b8wzlrm53Eyf9kCQpMcNYkqTEDGNJkhIzjCVJSsww\nliQpMcNYkqSS8mIS52bmGBuf7NhxDWNJkqi/mESWQW0YS5IGwtj4JAtFWChSNUxrLSbRzKpPrTCM\nJUl9b61hmvWqT4axJKnvNRKmWzYXVrymkB/pyNSZhrEkSdReTCLroDaMJUl9r9Ew3btzK0M5GMpd\nu5hE1qs+GcaSpL7XaJhu2pinkN9AIb9hxXOrBXU7GMaSpI44cGyCA8cmkh1/rWFaK6jXquYSipIk\n9YtymJa3u4k9Y0mSEjOMJUlKzDCWJCkxw1iSpMQMY0mSEjOMJamPpb6dSI0xjCVJSswwliR1lUHs\nzTvphyT1iXKAHdq9PXEl3auRtqn1mqza1p6xJKlv9Gqv2jCWJGVubHySczNznJuZY2x8MnU5Xccw\nliRlamx8klNnppcenzozzb6jJzn70mzCqrqLYSxJytTpiiAum569yJFHn0tQTXcyjCVJSswwlqQe\n1gsDlrZsLqzYV8iPtLSmcL+qeWtTCGEd8GVgEzACfBr4AfCnwPOll30+xvgnWRYpSepd+3dtY9/R\nk0zPXgQWg/jwnh2Jq+ou9XrG7wemYoy3A+8CjgJvBg7HGO8q/c8gliTVtHfnVoZyMJTDHnEV9cL4\nq8ADFa+9DNwC/HwI4ekQwhdDCD+VZYGSpN63aWOeQn4DhfwGNm3Mr/q6Qb0FqmYYxxjPxxhfCSHk\nWQzm3wG+BeyPMd4BfB/4RPZlSpJq6YcQW+stUL3cBnWnwwwhvB54DDgaYxwPIdwUY3y59PTXgCON\nHGh0dPVvQmoP2zh7tnFn2M6NGx7OMT17kXMzc0v7Tp2Z5sDnJ7j/V25leDgHrGzTZtv4w59+EoAv\n3f+ONdVa69inz1a/Bepzj3+X4w+8s+Znf/wLEyuCvNwGb3jdq1uuuVPqDeC6GXgS2B1jfKq0+4kQ\nwt4Y47PA24BvN3KgqSlv7s7S6GjeNs6YbdwZtnNz5ueLXL6ysGL/uZfnOHDkBFfmiwD81mdPsH/X\nNqC1Np4vfc5a/m3qfkax+u6FhWLd4/7FC1Mr9p17eY7f/eL/SDJYrNkvO/WuGR8EbgIeCCE8FUJ4\nCvgo8K9K229lcYS1JKnLlIMYemPWq0G+BapmzzjGeC9wb5WnbsumHElSs9YND3F5fmXveLnyrFd/\n9KbXdqCq5q3lFqgtmwvXnKYuv79XgtxJPySpx914w3oK+ZGlx4X8CLmE9axFq7dA7d+1bUUbHN6z\no+bI7W5iGEtSH1geYr16yrfRW6Cq6eV7mQ1jSeoDy0Os13uKrVhLkKdmGEtSn+rlnuKgqXufsSSp\nN5V7iuXt1A7t3p66hK5lGEtSjyrPOFXeVu8yjCWpB1WbOnIoB/nr1yesSq3ymrEk9aDTZ1ZOHblQ\nhNmfXGr7sXp5zudeYRhLkla11sUb1BjDWJJ6ULX7iLM4TV2tB16eyUvt4zVjSepB1aaOHMr16rxb\n1xrEUdeGsST1qL07t/Kph59d2s7i9qVem/O5V4Pc09SS1KM6MePUIM7klYJhLEmqyZm8sudpaklS\nTd02k1c/smcsSVJihrEkSYkZxpIkJeY1Y0nqY716q8+gsWcsSVJihrEkSYkZxpIkJWYYS5KUmGEs\nSVJihrEkSYnlisViJ45TnJpyIeosjY7msY2zZRt3hu2cPds4e6Oj+abWs7RnLKlvHTg2wYFjE6nL\n6IhGftZBao9eYxhL6gkGifqZYSxJUmKGsSRJidWcmzqEsA74MrAJGAE+DZwGjgMLwPeAPTHGjowC\nkySpH9XrGb8fmIox3g68CzgKHAYOlvblgPdkW6IkSf2tXhh/FXig4rWXgTfHGE+U9n0DeHtGtUmS\nNBBqnqaOMZ4HCCHkWQzm+4Gxipe8AtyUWXWSJA2AuusZhxBeDzwGHI0x/vsQwr+seDoP/LiRA42O\n5lurUA2zjbNnG3dGtXYeHs6t+txqWnlPr2rkZy2/pt7r1Hn1BnDdDDwJ7I4xPlXaPRlCuCPG+DRw\nN/BnjRzI2V6y5Yw62bONO2O1dp6fXxwn2sy/QSvv6VX1ftax8Ul+NH0BgI9/YYK9O9/UsdoGUbNf\ndupdMz7I4mnoB0IIT4UQnmLxVPUnQwgTLIb5I60UKkmNGhuf5NzMHOdm5hgbn0xdTs8ZG5/k1Jnp\npcffeWGKfUdPcval/v+S0ivqXTO+F7i3ylN3ZlKNJC2zPEhOnZlm39GT7N25lU0bPdXaiNMV7Vc2\nPXuRI48+x+E9OxJUpOWc9ENSV6sVJFK/MIwlqc9t2VxYsa+QH2Hvzq0JqlE1hrGktvrIH3yTj/zB\nN9v2eQZJffWuqe/ftY1CfmTp8V+7aQOH9+zwNH8XMYwldZUPf/rJa1ZnWh4khfyIQVJhtWvqywdn\n7d25laEcDOXg/l+5tdNlqg7DWFLXqwySRnvEgzICu9Fr6ps25inkN1DIb+ANr3t1p8pTg+pO+iFJ\nqVX2guv1iA8cm2Dm/CUuzy8s7XMEtrqdPWNJfacyiMv6dQS219T7gz1jSV2t8vqxVtq/axv7jp5k\nevYicPWaunqLPWNJfWfd8Mo/bf3cW2zlmrq6iz1jSX3nxhvWs1AsDkxvsTw4q7yt3mPPWFJfsreo\nXmLPWOpC5eukh3ZvT1xJ77K3qF5iz1iSVnHg2ETLA8jW8l4NHnvGknpCuZer1nmmpXvZM5YkKTHD\nWFLXmzl/aSCmttTgMowldbXVprZcvhCC1MsMY0ldrZmpLQdlcQj1H8NYUtuMjU+yUISFIh0Pw2pL\nCU7PznGlSphL3cYwltQWja6r2+xnVlNtastqSwkuFGH2J5daPr7UKYaxpLZodF3dWsbGJ/nR9AXO\nzczxG//6xDXhXpa/fh2H9+xwIg/1FcNYUldY3rM+P3elqfdXW0pwKAf569evuTYpa4axpLZY67q6\n1XrW1VxXZUUmWFxKsJAfuebYhfyGVV9fT68NBju0e7uTevQww1hSz6i36EO7FofI4vq3VIthLKkt\n1nrNuFrPOpe7uj2UW5wSs9a14vLiEPVeV087rn9LzTCMJXWFaqeZH/jlv7PU0/Xar/qZYSypLdZ6\nzRhWnmau7Om2cu231euoa/1Zeu16s9IzjKUu06t/yKv1bJu9BWnTxjyvefWr1nyaea3W8rN4vVmt\nMIylLtLrf8gre45rGUDVDVodDOb1ZrXCMJa6xIFjE1UnueilP+SbNuaXAqzXJ+Vo12AwqRHXpS5A\nkhrRK/fQbtlcWPGlqtlr5xo8DfWMQwi3hhCeKm1vCyH8IITwVOl//zjbEqXBsa7KICX/kPeWdlw7\n1+Cp2zMOIdwH/BPgldKuW4AHY4wPZlmYNIhuvGE9C8Ui07MXgat/yLN04NgEUL/n2ejrtHiN+VMP\nP7u0LdXTSM/4ReB9QPn2+1uAnw8hPB1C+GII4acyq04aQO2aRUrpeL1ZzarbM44xPhZC2Fyx6xng\n38UYJ0MIB4FPAAfqfc7oqL+QWbONs5dlGw8PL37ffcubXstrXv2qpe2slY9b72dr9HXlabPW0lbD\nw7ml9zd83Ir3rvX47fisdtaRhW6ta1C1MoDr8Rjjy6XtrwFHGnnT1FRv3JrRq0ZH87ZxxrJu4/n5\nIrD430rldtYaPVbDNRXXXvv8fHHp/c22RTvbbi2f1cl/w2b59yJ7zX7ZaSWMnwgh7I0xPgu8Dfh2\nC58hqUd1+7Xjbq1LqqWZMC6W/v8e4GgI4TLwQ+BX216VJJUYrhoEDYVxjPEMsL20/RfAbRnWJEld\nwS8C6hQn/ZDUVoX8htQlSD3H6TAlSUrMnrG0Rt0+oKkdyitJweIsYTfe4NrCUjsZxpJqWr6S1OX5\nBaZn5zj70mwmE1p86f53eNuNBo6nqSXVVG1JwIUiPbOSlNQLDGNJkhLzNLWkmqotCei82fX18xgC\ntZ89Y0k1LV8ScCiHCyBIbWYYS1py4NjE0ujwSpUrSeWvdyS11G6eppZUV3lJwJnzl3j5/CVgcZT1\n/l3bVrzW07NS8wxjqQtU3sc7Nj7ZlYE2c/4Sl+cXlh6fOjPNvqMn2btzq6espTXyNLWU2PL7eMsh\nd/al7O+1LX8JODczx9j4ZM3XVgZx2fTsRW9xktrAMJbWoJkwW021+3g7EXLVvgRMz85xpUroSsqW\nYSy1KGWPth1Wm8xj9ieXqr5+3fDKPxeF/Ii3OEltYBhLLWpXj3bL5sKKfd0YcjfesJ6h3NXHhfwI\nh/fs8Hqx1AaGsZTY8vt4OxVy1b4E1Lt1KX/9+qVbnLrty4LUywxjqUXt7NFW3sfbqZCr9iWgkN/A\ndVVOR5ddNzxEIb/BST+kNjOMpRa1s0dbvo+30yHXzJeAQ7u3d+UtV1I/MIylNUjRo22nVF8CJF3L\nST+kNSiHWXlbklphz1gqWW1eZknKmmEsCWjPBCaSWmMYS2Lm/KWensBE6nWGsSTnnZYSM4ylDvGa\ntKTVGMaSnHdaSsxbmyRx4w3rWSgWmZ69CFydwERSZ9gzlgT0/gQmUi+zZyytUb9MEdnoBCb98vNK\n3cSesSRJiTXUMw4h3Ap8JsZ4VwjhDcBxYAH4HrAnxljMrkRJkvpb3Z5xCOE+4CGgvDzNg8DBGOPt\nQA54T3blSZLU/xo5Tf0i8D4WgxfgzTHGE6XtbwBvz6IwqZOcClJSSnXDOMb4GHClYleuYvsV4KZ2\nFyV10tj4pFNBSkqqldHUlfPm5YEfN/Km0VGXl8uabdya02enV+ybnr3I5x7/LscfeOc1+9fSxsPD\nuZqfUe/5rFQeN1UNy6U+/iCwjbtLK2E8GUK4I8b4NHA38GeNvGlqyl5GlkZH87Zxq1YZfriwULym\nTdfSxmPjk/xo+gIAv/XZE+zftW3Faz7za28FOv/fSuVx5+eLSWqo5O9y9mzj7DX7ZaeZW5vKf7L2\nAZ8MIUywGOaPNHVEqcts2VxYsa+dU0F6GlxSPQ31jGOMZ4Dtpe0XgDuzK0m9oLzgQT9MALF/1zb2\nHT2Z2VSQp89UPw1+5NHnnHJSEuCkHxLgVJCS0jKMJa5OBVnIb6g5FWQrsj4NLqn3GcZSmy1ft3j/\nrm0U8iNLj8unwdsd+pJ6lwtFSB2wd+dWPvXws0vb3aofxgBIvcgwljqg0RWRJA0mT1NLkpSYYSxJ\nUmKGsSRJiRnGkiQlZhhLkpSYYSxJUmKGsSRJiRnGkiQl5qQfUomzT0lKxTCWOsSwl7QaT1NLkpSY\nYSxJUmKGsSRJiRnGkiQlZhiraWPjk5ybmePczBxj45Opy+kqto2kVhjGasrY+CSnzkwvPT51Zpp9\nR09y9qXZhFV1B9tGUqsMYzXldEXYlE3PXuTIo88lqKa72DaSWmUYS5KUmGGspmzZXFixr5AfYe/O\nrQmq6S62jaRWGcZqyv5d2yjkR5YeF/IjHN6zg00b8wmr6g62jaRWGcZq2t6dWxnKwVAOe33L2DaS\nWuHc1Grapo15CvkNS9u6yraR1Ap7xpIkJWYYS5KUmGEsSVJiLV8zDiH8T+Dl0sPvxxg/3J6SJEka\nLC2FcQhhA0CM8a72ltP9DhybAFwoXpLUPq32jH8auD6E8J9Ln3EwxvhM+8qSJGlwtHrN+DxwKMb4\nTuAe4CshBK8/S5LUglZ7xs8DLwLEGF8IIZwD/gbwv1d7w+hof9xzOTycA7rz5+lkTd3cDllq5Ocd\n1LZpJ9sue7Zxd2k1jD8EbAX2hBBeC9wI/LDWG6am+mMZufn5ItB9P8/oaL6jNXVrO2Sp0TYexLZp\np07/Lg8i2zh7zX7ZaTWMvwT8YQjhROnxh2KMCy1+liRJA62lMI4xXgE+0OZaJEkaSA66kiQpMcNY\nkqTEDGNJkhJzCcUmjI1Pcm5mbml7/65tiStSN3J2NknNsmfcoLHxSU6dmV56fOrMNPuOnuTsS94e\nIElaG8O4QacrgrhsevYiRx59LkE1kqR+YhhLkpSY14wbtGVz4ZrT1ACF/Ah7d25NVFFaXheVpPax\nZ9yg/bu2UciPLD0u5Ec4vGcHmzY6v6skaW0M4ybs3bmVoRwM5RjYHrEkqf08Td2ETRvzFPIblrYl\nSWoHe8aSJCVmGEuSlJhhLElSYoaxJEmJGcaSJCVmGEuSlJhhLElSYoZxxg4cm+DAsYnUZUiSuphh\nLElSYs7A1SQXSJAktZs9Y0mSEjOMJUlKzDCWJCmxgQ5jRzpLkrpB8jA2ECVJgy55GEuSNOgMY0mS\nEjOMJUlKzDCWJCkxw1iSpMRamg4zhDAEHAO2AheBj8QY/7KdhUmSNCha7Rm/F1gfY9wO/DZweC1F\nTM/OZXJ7k7dNSZJ6QathvAN4AiDG+AzwlrZVJEnSgGl11aYbgZmKx/MhhKEY48JqbxgdzVfdPzyc\ng1yO4eHcqq9p1fBwbtVjf/wLE5ybmQPgyKPf5VP3ZLMaU60a2q0Txxh0tnFn2M7Zs427S6thPANU\n/kvWDGKAqanZqvvn54tQLDI/X1z1Na2any9WPfbY+CSnzkwvPf7OC1P80j9/gr07t7JpY3t/QVer\nod1GR/OZH2PQ2cadYTtnzzbOXrNfdlo9TX0SeDdACOFngeda/JwkTlcEcdn07EWOPNreH2NsfJJz\nM3Ocm5ljbHyyrZ8tSeofrYbx48BcCOEki4O3frN9JfWH5b3vU2em2Xf0JGdf8tuoJOlaLZ2mjjEW\ngX/W5lo6ZsvmwjVBCVDIj7B359a2HaNW7/vwnh1tO44kqfcN5KQf+3dto5AfWXpcyI9weM+Otl8v\nliSpEQMZxgB7d25lKAdDOdraIy7bsrmwYl+7e9+SpP4wsGG8aWOeQn4DhfyGTHrE9r4lSY1KGsbl\n0cYLRZg5fyllKZnIuvctSeoPycJ4+Wjjy/MLbR1t3A23FWXd+5Yk9YdkYZzlvb7eViRJ6iV9ec24\nU5N6SJLUDsnC2NHGkiQtShbGy0cbD+Vo22hjg16S1EuSnqauDMeFIm0baOVtRZKkXpI0jL/6X168\n5nE7B1p5W5EkqVckDeMsB1p5W5EkqVe0up5xXzi0e3vqEiRJStszdqCVJEmJw3j/rm0M5a4+dqCV\nJGkQJZ/0I3/9esCBVpKkwZU8jK8bHmIohwOtJEkDK3kYS5I06AxjSZISM4wlSUqsr+8z9j5iSVIv\n6Osw7gZ+IZAk1eNpakmSEjOMJUlKzDCWJCkxw1iSpMSSh/Gh3dsp5DekLkOSpGSSh7EkSYPOMJYk\nKTHDWJKkxJqe9COEkAN+ADxf2vXfY4wH21qVJEkDpJUZuP428Ocxxr/friKcpUqSNMhaCeNbgL8Z\nQvgmcAH4zRjj83XeI0mSVlEzjEMIHwY+umz3buD3Y4yPhhB2AH8M/N2M6pMkqe/lisViU28IIbwK\nuBJjvFx6/IMY4+uyKE6SpEHQymjqByj1lkMIPw38VVsrkiRpwLRyzfgzwB+HEN4NXAE+2NaKJEka\nME2fppYkSe3lpB+SJCVmGEuSlJhhLElSYoaxJEmJtTKauiEhhCHgGLAVuAh8JMb4l1kdb1CFENYB\nXwY2ASPAp2OMX09bVX8KIfx14M+BtznrXPuFED4G/D1gHfC5GOPDiUvqK6W/yV8E3ggsAP80xhjT\nVtU/Qgi3Ap+JMd4VQngDcJzFdv4esCfGWHO0dJY94/cC62OM24HfBg5neKxB9n5gKsZ4O/Au4HOJ\n6+lLpS89/xY4n7qWfhRCuBN4a+nvxZ3A30paUH96B3BDjPE24HeB30tcT98IIdwHPMRihwjgQeBg\n6e9yDnhPvc/IMox3AE8AxBifAd6S4bEG2VdZnIgFFv89rySspZ8dAj4P/DB1IX3qHcB3QwhfA74O\n/MfE9fSjC8BNpZX3bgIuJa6nn7wIvI/F4AV4c4zxRGn7G8Db631AlmF8IzBT8Xi+dJpEbRRjPB9j\nfCWEkGcxmH8ndU39JoTwQRbPPjxZ2pWr8XK1ZpTFRWj+IXAP8JW05fSlk8AG4H+xeJbns2nL6R8x\nxse4tiNU+TfiFRa//NSUZTjOAPnKY8UYFzI83sAKIbwe+CbwRzHG8dT19KEPAT8XQngK+Bng4RDC\nzYlr6jf/D3gyxnildD1+LoTwmtRF9Zn7gJMxxsDV3+P1iWvqV5VZlwd+XO8NWYbxSeDdACGEnwWe\ny/BYA6sUCk8C98UYjycupy/FGO+IMd4ZY7wL+A7wSzHG/5u6rj7z31gc80AI4bXADcC5pBX1nxu4\nerZymsWBcsPpyulrkyGEO0rbdwMnar0YMhxNDTzOYm/iZOnxhzI81iA7yOIpkAdCCOVrx3fHGOcS\n1iQ1Jcb4n0IIt4cQvsViJ2F3vdGnatoh4A9DCP+VxSD+WIzxQuKa+k35d3Yf8FDpzMMp4JF6b3Ru\nakmSEnNAlSRJiRnGkiQlZhhLkpSYYSxJUmKGsSRJiRnGkiQlZhhLkpTY/wdV0O/5C0q7rQAAAABJ\nRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10a033c50>"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We're going to fit a line to the data, as we've done through the lecture:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def model(theta, x):\n",
      "    # the `theta` argument is a list of parameter values, e.g., theta = [m, b] for a line\n",
      "    return theta[0] + theta[1] * x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---\n",
      "\n",
      "We'll start with the assumption that the data are independent and identically distributed so that the likelihood is simply a product of Gaussians (one big Gaussian). We'll also assume that the uncertainties reported are correct, and that there are no uncertainties on the `x` data. We need to define a function that will evaluate the (ln)likelihood of the data, given a particular choice of your model parameters. A good way to structure this function is as follows:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def ln_likelihood(theta, x, y, dy):\n",
      "    # we will pass the parameters (theta) to the model function\n",
      "    # the other arguments are the data\n",
      "    return -0.5 * np.sum(np.log(2 * np.pi * dy ** 2) + ((y - model(theta, x)) / dy) ** 2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "What about priors? Remember your prior only depends on the model parameters, but be careful about what kind of prior you are specifying for each parameter. Do we need to properly normalize the probabilities?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def ln_prior(theta):\n",
      "    # flat prior: log(1) = 0\n",
      "    return 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can define a function that evaluates the (ln)posterior probability, which is just the sum of the ln prior and ln likelihood:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def ln_posterior(theta, x, y, dy):\n",
      "    return ln_prior(theta) + ln_likelihood(theta, x, y, dy)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now write a function to actually run a Metropolis-Hastings MCMC sampler. Ford (2005) includes a great step-by-step walkthrough of the Metropolis-Hastings algorithm, and we'll base our code on that. Fill-in the steps mentioned in the comments below:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def run_mcmc(ln_posterior, nsteps, ndim, theta0, stepsize, args=()):\n",
      "    \"\"\"\n",
      "    Run a Markov Chain Monte Carlo\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    ln_posterior: callable\n",
      "        our function to compute the posterior\n",
      "    nsteps: int\n",
      "        the number of steps in the chain\n",
      "    theta0: list\n",
      "        the starting guess for theta\n",
      "    stepsize: float\n",
      "        a parameter controlling the size of the random step\n",
      "        e.g. it could be the width of the Gaussian distribution\n",
      "    args: tuple (optional)\n",
      "        additional arguments passed to ln_posterior\n",
      "    \"\"\"\n",
      "    # Create the array of size (nsteps, ndims) to hold the chain\n",
      "    # Initialize the first row of this with theta0\n",
      "    chain = np.zeros((nsteps, ndim))\n",
      "    chain[0] = theta0\n",
      "    \n",
      "    # Create the array of size nsteps to hold the log-likelihoods for each point\n",
      "    # Initialize the first entry of this with the log likelihood at theta0\n",
      "    log_likes = np.zeros(nsteps)\n",
      "    log_likes[0] = ln_posterior(chain[0], *args)\n",
      "    \n",
      "    # Loop for nsteps\n",
      "    for i in range(1, nsteps):\n",
      "        # Randomly draw a new theta from the proposal distribution.\n",
      "        # for example, you can do a normally-distributed step by utilizing\n",
      "        # the np.random.randn() function\n",
      "        theta_new = chain[i - 1] + stepsize * np.random.randn(ndim)\n",
      "        \n",
      "        # Calculate the probability for the new state\n",
      "        log_like_new = ln_likelihood(theta_new, *args)\n",
      "        \n",
      "        # Compare it to the probability of the old state\n",
      "        # Using the acceptance probability function\n",
      "        # (remember that you've computed the log probability, not the probability!)\n",
      "        log_p_accept = log_like_new - log_likes[i - 1]\n",
      "        \n",
      "        # Chose a random number r between 0 and 1 to compare with p_accept\n",
      "        r = np.random.rand()\n",
      "        \n",
      "        # If p_accept>1 or p_accept>r, accept the step\n",
      "        # Else, do not accept the step\n",
      "        if log_p_accept > np.log(r):\n",
      "            chain[i] = theta_new\n",
      "            log_likes[i] = log_like_new\n",
      "        else:\n",
      "            chain[i] = chain[i - 1]\n",
      "            log_likes[i] = log_likes[i - 1]\n",
      "            \n",
      "    return chain"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now run the MCMC code on the data provided."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "chain = run_mcmc(ln_posterior, 10000, 2, [0, 1], 0.1, (x, y, dy))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plot the position of the walker as a function of step number for each of the parameters. Are the chains converged? "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(chain[:, 0])\n",
      "plt.plot(chain[:, 1]);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAFVCAYAAADPM8ekAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd409b6x7+ynb1DHHYgQGL2KrTsAt2le962t7/bPW73\nHnTQXSjde+9Je4G2dNCWsPcIG4dAICF7T8fx0O8PW7IkS7I84hh4P8/Tp44t6xxk6bzn3QzLsiAI\ngiAIInLRdfUECIIgCIJQh4Q1QRAEQUQ4JKwJgiAIIsIhYU0QBEEQEQ4Ja4IgCIKIcEhYEwRBEESE\nYwjmyyaTKRPAFgCnmM3mgtBMiSAIgiAIIQFr1iaTKQrA+wBaQzcdgiAIgiCkBGMGfwnAuwDKQzQX\ngiAIgiBkCEhYm0ymawBUm83mpe63mJDNiCAIgiAIEUwg5UZNJtMKAKz7v9EAzADON5vNlXLHsyzL\nMgzJc4IgCOK4IaRCLyBhLcRkMuUBuNlHgBlbXd0c1DiEOkZjEugady50jcMDXefOh65x52M0JoVU\nWFPqFkEQBEFEOEGlbgGA2WyeEYqJEARBEAQhD2nWBEEQBBHhkLAmCIIgiAiHhDVBEARBRDgkrAmC\nIAgiwiFhTRAEQRARDglrgiAIgohwSFgTBEEQRIRDwpogCIIgIhwS1gRBEAQR4ZCwJgiCIIgIh4Q1\nQRAEQUQ4JKwJgiAIIsIJupEHQRBEKLDYLWCgA5DU1VMhiIiDhDVBEBHB/SufhI7R4bvL3u7qqRBE\nxEFmcOKYwMk6u3oKRAg4Fn5Hh9OBg42HwLJsV0/lmGJ37T40WBu7ehpeHGkuQ5utDU7WCYfT0Wnj\nkGZNaKLWUo/DzSUYmzmyq6fixTf7fsKWyu14dvIjYFkgPioubGOzLAuGYYI+j81px4GGIuSmDYSO\nOf720O32dtHfTtYJB+tElO7oW6LuXP4IAOCc7DNwVvYpXTybo5sOhw01llroGR3e2f4JEgzxmDdt\nTqePW9/egFhDDOIM6mtJcfMRzN30BkYbhwNgUNJcitkn3YsYfXTI53T0PQl+svLIOmTEpWNoN1NX\nT6XTaLQ2w9bchijEd9oYL2x6FRZ7O7qfeA8y441ei+imim1IjUlGTtrATpuDHFVtNVhTtgGAy4wK\nALeOvBbDM4aEZfzH176Anondcduo64M6z+LC35B3ZDUA4MbhV2N05ohQTO+ooaD+AP/a6XTijryH\nAQBvz5zXVVPyi0ZrM2rb6zAgpR//3q9Ff+LXoj/x1oy5AW3oGq3NqGqrCvsz1dVUtFYiOToJ8VHx\neGf7x9jfcBA9E7oDAFrtbfxx++sPIiMuHWmxqSGfw2Nrn4eB0eP1GS8oHuNknZi76Q0AQH71Lv79\ne1c8hrvH3AyjcXRI53RMC+sGayO+L1gI4Oh56APh0TXPAPDv37irZi8W7P8Z94y9BakxKT6Pt7g1\nn1e3vgeL3YKpvSdCx+gwpddJiI+Kw2d7vgUAvDH9Beh1+gD+FR7Wlm3E1/t+RLfYdFyccw7iDHHI\nVViw/jj0j9d7Gyq2hE1Y11sbUG9tCPo8u2v38a8/3PXlMX2/SilsKML7Oz/n/3awHlNifXtDpyzG\noebxtc/DwTpgjOvm9dlP+3/BJbnn+TwHy7L4X+GvWFayCgAQZ4iFxd6OZyc9elRcAw6H0wE76whI\nuyxqLMb8LW95vV/eWsm/bu5ogY7R4bVt7wEA3pzxYkisUQcbD6G4qRTT+kwEANhZdZN2dVuN4mev\nbXsfk3LfDXpOQo5pe5twt3484I+/790dn6LGUovZa57zy7dmsVsAAKtK12HFkTV4buMreHf7p/zn\nwp1voHy970cAQG17HT7Y+QVe3/Y+CuoLZY8tbDjo9d7Wqh1o6WgNeh5K1Fhqsa5sE4oai0N2ziqL\n+ME/3FQS1PlYlsXSw3nYV7c/qPN0JiXNZbht2YN4dat4UdtStpN//dja58M2nzabRSQUtNLc0cJv\nMKottV6f5x1ZjcKGItiddtnv17XXY3ftPpGgBjwb5OaOFr/nJEeNpRZtNtfz+VvBMuSVrA6Jj7Ww\noQjvbP8EbTYLFhX+hjuXP4J7VzzGz18re2sLZAW1lN8P/QOrw8r/fUfew7ApXFsl1pRuwG3LHsR9\nK55AVVsNDjQcwstb3sGC/YtxoKGIP+5g4yHFc2yr3qn4WWdwzGjWRY2HYXPaEKWLxhv5H+D+E25D\nSnRyV08rrDhZp+IOk2VZlLdW4nBTCaSieWvVDpzQfZTiedVuWAA40lLGv35k9TN4edoziDXEaJ22\niFWl62Tfr2itQm7aIK/3GzuaZY9/aPVTnaKdOlknnlw3V9OxHY4O7KzZg5HG4Yq+12XFK2X/DfM2\nvxnU/Bs7mrD4wO8AusaqpHYvlraU453tnygGC72y9kPR37PXPIcbhl+N7JQsAIDNYYNBZwhJrICQ\nZzfMR2NHM16aOgfxUdpcSnanHQ+vftrnca9ufRfZyf1w/7jb+PccTgfe3fEp9tYVqH537uY3FLXH\nNlsbHlg1BxN6jsPVQy5TPIfD6cCT6+ZiQEp/3D3mZny2bQEA4Mf9P2Nougm3jQ7MjcOyLL/Z+nLv\nD9hRs5v/LK9kFYakmxCtj0LvxJ78+/XtDahoq0L/5CzEGWL599/a/pGmMUuaj3htqu5e/igAbZa9\n5o4WfGP+CQDQ7mjHU+vFz8dr297nX++o3oPNlflIi0nFaf2m8++zLItfDv6pab6hIuI0ayfrxObK\nfL93SvO3vI3Xt32A+VveQoejA89vfBWVbVX853trC/w+ZzDUtzd0amSgHPM2v6m4C//54B94buMr\n+GrfAny9b4Hosx8KFnkdv658M/5X+CtYlsXLW97xax7c4sOyLBYV/ob99QfRYG1EZWsVihqLUd+u\nbDb+zrxQ8zjNHS2K2gqATrn+c9b5FnxVbdU42HgI/ytcgk92f4M/Dy2TPa6g/gB+KvwVfxevkP38\nrrxHNM2pzdaGz/d8hwr3AtZgbfQK2AonDdZG3JH3MH4VLGYsy2Jd+WY0WBvxy8E//YrqbbA2Yv6W\nt9But8Jit+DuFbPx8e6vNX3X4XSgwdqIpo5m3irEsixWla5DraUOAFDZWoVfDvzBb5qaO1rgZJ2K\nc3SyTny77yfsrSvAXW4hoYWipsOivw82HvYpqDnWlG30es9it+Aftya+vnyz7Pc2V2zDu9s/4TX+\ng42HcP/KJ0TH7Kkzo0lh06vEkeYy3LviMT62AIBIUAPAkqK/MH/LW3h+46ui9x9b+zzeyv8I9698\nAqUt5X6NC7iu2zvbP5H97K/i5ai11OOvw8sxb9ObONjoueYsy6KwoQi7BC4nX/xVvBwrjqzFogO/\n4au9rnVzVek63J73kN/zDpaI06zXlG3Ed+b/YUTGUNw84j+ads/Vbd6mJwD4XiCE3tr+EWb0maLJ\ndxQsR5rL8MKm1zAmcyRuGP5v1WMrWqsQo48OyCd1sPGQSJCWtpTj8z3f4fbRN4iO40yiSrTYWlHS\nXIbk6ES02S3omdAdX+39AYDLNOgvvxX9hTGZI7C5Mh9/FS/HX8XLvY5JiU7Gc5Nna9aOvi9YhDGZ\nI5EUnQgAqLXU4Yl1L6p+p8PZgTidtsjw78wLsap0HU7qcQL+b+jl/PvSaO/a9jrZ7zucDn5H/9T6\nl0SfbanMxzkDTgfgWujXlG3AiIyheF2wg5fDzjpQ2lIOg86AWkudYpDkX8UrsLFiKzZWbMXgtBzs\nqxebvkMVsS495xvbPkBBwwGcO+AMpMem4cQeYwG4An8Al7nynAFnAHC5NtaVb5Kdn1buW/k4rz1u\nq9qBitZK9HAHHnFsrNiKxKgELCxcgjP6zcCqsvUoFJg1n5zwIKotNfjOvBAJUfGYN3UOnt4wX3SO\nN/I/FAnqpyc+gm5xafzf7+74FHtqzVjtDmz0hwUFi7H8yBoAwPjuYzR/T/rrlbaUewlBKctL1mDB\n/sUAIBJQHU6b17Ef7fwK955wq+x5WJbFjpo9yE0byGvCL2x6TfPcAZfZeWKv8V7C+aOdX+LJiQ/6\ndS41fjn4p0jjfXnL2zgtazr+Kl6OnNQB2C/jNtPKuvJN2FW710shCva8Wok4YV3k3gntrNmD2/Me\nwn9HXYdh3QbLHttma0OcIQ5z1mszS+YdWY2x3UeJIjZDTbvdyt/I26p2wGJvF5l6pDzjXig4U2VZ\nSwX0jA6Z8UbVBbbWUier8e6tK8C+uv3ITRvI73ovGHi2z3n/XbwcmyvzAQCZcRn8++vKN/n8rpSy\n1go0Wpv4oDM5GjuacHveQ7h15LXIiOuGZze8DNbLQC9mQcFi9ErsidOyTsaHu770OY92u9Vn6kWH\nowNRuije/L6hYgsvrOdvfhvlrRWYN3UO9Do96izKFoE7lz+iaG7mfkcn68QrW95BUVOxZguCcEF+\naepTXmlpFrtFtBGTE4S35z0UclP4gv0/o6DBFRPCLY5Du5mQGJUgCsxpsDZi6eHl/H0UqKDm+NK9\niQSAZza8jNHG4Tg7+zT0TuyJ3bX78Pme7/jPP5W5/4Qmz1ZbG9bJaKRSjfqJdS/gvAFn4uQ+kxFr\niMGeWnPA8+cENQBsqtym+Xvfmv+HKb0noK69HqkxKbKCut3ejlj3WvPrwT/xu0zwpRIHGotQ114P\nm9OOzLgM0dqzvWY3Ptz5BVJjUvDc5NkBWSi/Mf/Em56FVFlqRBvdzoBTFkIhUOUsl7ePvgGvbX0P\nRU2hi2GRQz9nzpxOHcDNnLa2DtUDmjqaUd5aiR8KFover2qrwdjuI8GAEflsNlZsxUtb3sJvh/72\nayLryjdhfPfRSIhKUD0uv2onDjcdQZ+kXorHOFknNlVuQ3psGqraqvHixtdQY6lFicCHu7Zso8jX\nIeW3or8AALOyT8OWyny8kf8BVpSuRY+E7siIS0eb3SIbVfn1vgWoEJj5hWys2Io2u4UPUtKyQJa1\nVvCvQxEk9k/JSk3HFTYUoaD+gNcCeWKPsbA7HWi1eQLFylsrUVBfiIrWKpgVAs6ELCtZhVP6TsU7\n2z9BckwSH6nb3NECm9MGi92K2Wuew69FYt9T78Se+HjXVyhpKYWddeD3Q//A4XTg9Y3qPrUz+83E\nt/t+Ev3+AJCbNghjjCOwr24/lspYGYQ8MeEBrDiyVvYzHaODKd3jt2/uaMEfh5ZpWiT21poxvsdY\n6EMQNcuyLN7d4W2GbLA2YpRxGOZuep1/b1nJqoCC5f5lulCTubKirQqrStejrKUiIB+i1HSrhLm+\nEEsP56FfUh9+Uyvkucmzcd7As/jshIsHnaPZxK2VrKTeeGnzW/hdYc1LiIpHdnI/3J73kMiaoJW8\nktVYcWQtWu0WXkFqsbXy6UntDiuGpueCBbDcnWYYCqb1mYRoXbRoLX9q4sM4K/sUXDhoFr9Gqj0b\nXcm5A87ApF4nYlb2aRiXOQorSl1zvHT4OU+FcpyIEdb3r3ySz5cV0tjRhKWHl2ND+RbMzJrKv+/L\nBKTG0HQTMuONXu+zLIvP9nwLB+vAx7u/xvaa3ZiVfRoAV3Tu70V/Y0h6LhgwONRUgsfXPo/t1btQ\n116P7wsWwero8FqoO5w2nNJ3Kh+pXdpSgUNNxeiRkIktlfl8ROGZ/WbieYFpaVv1Tvx5eBn+Ll6B\nVlsbhqTn4kDjISRHJ0HH6LC6dD1q2+sV/41aFsj5057C0sPLfR4nZUTGEMzsOxVFTcV4fvJspMak\nYE9dYNpGu8Mq6xt89MR7cHKfSfyDKkRuk/LilCdw3oAzvbSJvw4vR017HW8eXVz4G74vWIS/i1fg\nn5KVshH0W6u2o9km3kEfaPS9+P1+6B+v3x9w5Y2WtVZCz+h8bpwuy70AZ/c/VXYTeqCxCJN6juct\nNY+ufkazttBgbUSbrc0rpY1lWSwsXAKDziAy9arx/MZXva4P4NrwZSX1xpaq7ZrOo4QxPh3/GXqF\n7G+vhNLGNdTICWoAuDjnXOgZHRraG2DQGXC56UIcbipBu70dFwyaJUrNU+LpiY/Axtpx/fCrRBq4\nr7E5esRnoltcOlaWahNoD427U3bNPdxUglnZp6G5o8UreG5t+Sakx6aGdCMyNnMkPt71Fb8OvDDl\ncaTFpiLaraQMTOmPyb1OQlZSb8V74pKc8wJeg+Tg1jctcDICABKjE3BW/1PQJ6kXBhr7hlRYR4QZ\n3Mk6faYd1VsbsLvWjNzUAfhKEiDlLztq9sjm4VZbarG5Ml/0UDicDqwp28D7v1eXbcDkXieKAj62\nVu1QHe8+SUAHANwz9lZ8svsb/m+u6pEcK46sgY5hkFeyGgwYjO8xhjdBBoMvE7GU0/vNwNLDebhw\n4Cx0T8jE1N6ufMTpfSdjbPeRKGkuVQz8CJQRGUOws2av6jHjuo/mfdlShKZ1uUC6cJFfvRP5GlM9\nGIbBlYMvxjf7vM2GL256HXOnuoq/tAvSV7SwsnQdzux/KlJiXI0y9tXth7m+EP+UrMQ/JSt5UzmX\nyic0hVa2VoFhdEiNSRFZYaQI86W1MiJjKKyODiRFJWBL1Xakx7s2Ddz9djRx1ZBL+dfC2BGle69f\ncl+kRCdjUGo2usWl4QrTRQCAl6bOwQOr5vg19vIjazCp14majp2SNR5ZyX1Uj1FKY5P7TeZMeMjL\nHXmF6SJ8a/6fz7m8tFmcrpUcLW7kMjg9h399Rr+Z+POwd8DmjL5TsOLIGtnUOSWidAYvk/4TEx7g\n3QDCNDol5FyMOkbnrmgWWgIS1iaTSQ/gQwC5AFgAt5jNZm32JBna7NqCmN7Z/jFy0wYp5txyJETF\no9XWhmuHXoGUmBQ+eZ5jTdkGXDn4Yq/vyeXsyglRaWRmIPWMpbmlvsgrcZmdWLDYWLHV6/MBaVkY\nmT4ciw785vdctPDI+LvRJ6kXzh94luznydFJGJpuQqw+xm8hosZ1w/6Ne1bMVvx8aDcTLsu9gP87\nwRCPVnsb0mPTUKdieZBjYs/xmnz0z09+HOWtFShvrcSP+3/2aww5ZvadimUlq0SL+4Qe42SFdYut\nFbctexD/N+Ryr8+08OiaZzCt90RcknMe3swXp0g5nA7oGB0f6fr2zHluk/envGbYN6l3QOMqISz4\n0dTRDIPOgKtPuACwAOcNOBOTep6oOSalqxhtDLza3KU55yFbJoYmPioer09/Hr8eXCoboKmEFovj\nqVkn46aJ/0J1tXIE+O9F/2B12XrZz1ps3jUM0mNTMfvEexFniEWb3YLk6CQkRSdiWclqUVZOsMi5\nBIemuwIvn5zwoM8o7QEp/fgI8StMF+OLvd/j4pxz0SuhBwalZsPgR3nbW0ZegxEZQ/2YfXAE6sA6\nB4DTbDZPAfAYgOeCmURJc6nmY9UE9aiMYbhu2JWYN3UO3poxF+N6jEFO2gDFADUpXDGOo5FJWeM0\n7apnn3gv7h5zC56d5Eo5GZiSLXtcVpJn182ZdXzBMAxfFjAYLsnxROxH66PQK6GH4rG3jboeCYKc\n2Mcn3I/HT7oPFw6a5fe4M/tO9X0QgJSYJAxOz8GMvlNw3bCrFI87s796XeispN6474T/4uKcc/H2\nzHkYkp7Lf6bX6XHH6BsVv/vF3u+93jupxwkaZu/SsOU2fHcuf0S02N227EF8ufcHkQnXn2dVC8Is\niOToJPzf0MuRmegKcGQYBsb4brh15LUhGevCQbMwqac2zRNwBVo+MeEB1WMuGHg2/j3kEp/nmtJ7\nguz7coKaw6AzYFx375KVsfoYvDzNd263kDdnvIhnJj2CU7Kmie7LcwecKXv8r0Xa0usGpPTH0xMf\ngV6nR6/EHkiLTUXvxJ68peuuMTfJfu+qwb6vmRxR+ij+9fxpT+OSnPN4xYthGNwz9lbcOPxq/pgs\nweayf3IW7h37X/7vERlD8cb0FzCz71QMTs/xEtQGxhX09sj4u3mTvBCpBaCzCUizNpvNi00m06/u\nP/sD8E+FkfBWvrZkeDXOG3Amzug/k/9baMK7deS1WFi4RBT01OGwufJRHe3ISuoT8m4/b82Yy6e2\ndDYPjbsTo7NzUVxRLXr/phH/wQdus+SVpouREpOMXoliwZca4ykcM6PvFIwxjsQ/JStxee6FfBlT\nf3abqbGpgMDXc8/YW72sCNcOuxKfClwAHHePuQUx+mgv89ywboNFptdYfSzaHe1IlAkSTIpORFJ0\nInokdMcaP9KELhh4NtJjfftsOTMlxwndR8EY1w1zN7/hdezZ/U9FjC4aiw/+7vXZYyfd53NjIzT/\nXZxzLn7a/4vq8XLzv2rwJRhtHO5lUj0i41uXY0PFFtXPn5s8G43WJszb/KbqcSnRSeiT1Fsk+DPj\nM1S+4WF4xhAMSs0WBU29NWOuaGMRrY9Gh6MD1w//N9JiUpESk4TH13rqOs/oMwWnZp0MALg093wY\ndHpRjjDgWtjjDHG4eshl2FyZj8HpOegeb8RrJz+Hu2WsO7lpg1SDR4VcYboIl+WcjyfXzUW9tQGD\n03Iwve9kn9+T2ySPzhyBWEMsTu4zGStkfNuAS6DPn/Y0H+GtY3RIj03DRYPOER13Zv+ZGNd9NJ70\nkQYpx/8NuRxjMkciWiBApaTEJCMrqTeKJZu8MZkjXPnqGszkQganeZ6JOEMsZvSdIvp8UKpL+Tgh\ncxS2VG3H/SfcjoL6A8hJG8CvY7OyT4PNaffZ8OfRk+7F4aYS9EnqhacmPoRHVrvWw1tGXoMGayP6\nJff1a+7BErDP2mw2O0wm02cALgTg9zap1daGaH00HCEqVKL20DAMgwsGnY3GjibeHy00rb4x/QUU\nNx8JyTwAl7bKMAzuGHOj14IQauZMeAjG+G7Q6/RIMMRjaLoJQ7rlemmJidGJsn76xGiPwDt3wJmI\n0UdjYGp/0TEdDu+8TCXOG3AG7E47ctMGwpQ2SFS5aO6UJ1FtqUX/5L6ywjonbYDsOYWLAVepqdHa\n5HMuFwyahRc15oNK75+h3UyyKTpyGlJWch+vBemc7DOg1+lxev8ZvLDOSR2A3ok9cUnOeZrznh87\n6T6Y6wpxcp9JisL6v6OuQ7QuCsb4DOyq3YuLBp2DHwoWoby1EmmxqdDLbLbkApj8ZUKPcUiNSUFq\nTAoSoxJE5tG3Z87DoaZivLT5LZzV/xQ+3/q2ZQ+653w9sv1Y7G4ffSPa7e0ob63EwJT+YBgGD4y7\nHe/t+AxzJjyEemsDNlfmY7RxOJ81wlX9arO18SlNgOd+mn3ivdhTZ8bCwiUAgKm9J/LWKeH9EKWP\nwtsz5/FzB1wL/ln9T/XnckGv0+PZydqLqHD0Tuwpyk/+l3vDeFnu+bg05zysKdvgJfRenPIEGIZB\nt7h0n+fP0HDMWzPmitwhAHBST22WnIfG34XvzYugYxg0WJuwr24/YvQxSJLRTKf3Ud/A9ErsgblT\nn0ScXjkdFgCuG34VroPL6jWkW67os7MFAWFqdI83ors7EDk5OglvTH8BDMN0WVe8oALMzGbzNSaT\n6SEAG0wm0xCz2azofDYaPT+M0+nEbQseRPdEI6b18zZL9UzKRHmzNj/HM6fcDz2jR/duvptRPGC8\nCdcuvM+r0IdacJdWXjpjNjITMmB32pEU4wl2unTYLCzYvUT1u3Nm3Is5ea8ENO7Qfv3515mZyZhz\n2t2yxyUmRYt+A44hzQP5dIgexhQY9N63RFl7mex35TAiCcP63SF6771zX0C0PgqJMQnIhkuzn3f6\nbDy41OM9mTPjHsUxJjKjscQdBdpkb4TRmAQjfM/HaDThDt21aLQ24Yt8b/+v+FjX+T698GU0d7Qi\nPS4Vi/b+iR8lv53SHC8beQ7mr3kffVN6Yf4Zj4mE8cVDz0ZyTCLOyp3hc85y8xrZ37vMqpDpg8fz\nr1/u8xgAYHjWAOyu2o+JfceqVnkLhhsnXM7f62+d8wy+2bEISw+sxLDMXNdvZByGHwaKrSpT+52I\n7RV7MCV3DAwqubVK13kAegqOGYbxA10FaPoiQ+U6yZ/LaEzCqOwclFpKsbE0H6OycmFMU76vrhp5\nIb7e4cqP75vRHZmZ4Sln/NzpD+Ca/93L/92ru9iCMkw3EBDsK2884Ur06uFbAGt9pgHX2nLHpP/g\nll8eQW63AXj6lPv8Elq3G12maWGBnqS0Ufi+IAmNVpfv/OtL3hCZuRXnreHZPxYJNMDsagB9zGbz\nCwAsAJzu/xQRBjN0OFxpXJUt1SJB9sj4u/Hj/p9xzbArYHPYkR6biiVFf2FQajbe3v6x6HzXDr0C\nte31SGczARaqwRJCUqNT/KrKdePwq/kCHCnRyWjscGl0aTGpuGLwRXz0c7wtBS0NLg20HZ65TMqY\nhAVQF9ZGpgdePflZfG9ehGl9JqqaFIXzuTz3Qv7fbTQmqV6Dthab7OfDEofh7OzTYIzrhvo68XU5\nq/+p+P3Q3xiWOlTz9ZVHDwucsAiuSwJScEnOedhduw/XDbsK8Uyc4hhp8KTZGWOMfs1lcMIQIAE4\ncYZrU+hgHWizW/DshpfR6m5q0Dept+icesSi0dKOGd1PxpjU0Wi0NmF16XrcOunfimNnxwzEy9Oe\nRqwhFjU14rSmmT2mA9B+j/qL/HkZ5MTl8nPJTu7nVe4yGB476T60N7Gie/38fufg/H7nqMwJuGzA\nRbio//mor1XO5fd1L4eaKwddijP7noYEe6rquJMyJuJruIR1Cpse1jm+PO1pvL7tfczoO9VrXOHz\nMaPPFIxOGe1zbtJrfHHOuUiOcrmPpNXJLsu9wH2sHvOnPYVYfSxqa0LTKOf5yY9je/VuxBli0VDX\nDqDryuSGGn82Q1oIVLP+EcBnJpNpBYAoAHeZzWbNIcBKO/0+Sb1w99hbRO+dN9AVAPHmjBcxe81z\n6JXQA7eMujbgpvRqaScA8OSEB0TlIoV9hYd0y+Vr8D496WE+TUCtIlq0PgojM4ZhR81unNFvpque\nNWvnA3y4Eo3R+mhcPVRciH9W9mm8Rgm4Ug1GZ47A/GlP40BDkaYe3Q+NuxMbK7byEZNSdIxOlCco\nHX+UcTj6CEzZoWRG3ylePicl7h37Xyw68JvfpkcObjdvYAxIjk7CvKlzsL16NxYULMa9Y+XLLALg\nzbz9kvsuVs2OAAAgAElEQVQixhANQPk2j1WpVBdKTGmDkBiVgA5nB2Zln67pO/ePuw1f7vkBabEp\nODXrZGyqzEeP+EwkRMVj4YEl2FNrxoCU/qhrr5cNLIrSRcHmLlOZEp0UcCChjtEhWh9ZLQmi9FG8\nudMXqTEpaLA2IjXWtyUvlMQaYvHQ+LsUP39ywoPYVbvXpxlZCaHbbN7UOWiwNiItJhV76swYmzmS\n/8zfdE8tjDIOC/k5j0UYf9ojBgHL7eIcToes2blvUm88rHIzAi7fqUGnD8pnIPQ7ycHlmq4t2wQd\nw2BCz3F4deu7KGwowoWDZvH+LWF50Iy4dNloQQ6bw4by1kqfeY0cz254GeWtlbh5xH9Q217Ppwc9\nOeEB2WIuQPi1keORrr7Gq0rX4TvzQk3Baf7CrQPFzUd4y060Lgq5aYOQkzYAp/SdBqvDigUFP+OU\nrGlegYqhpKuvsxqttjbUtdeHPIUt3ETyNT5WMBqTQlqQP+zC+qFVT3nl6d015ibZ9oedwZKDSxVL\nlPZO7IlHT7zH6/02mwXba3ZjfPfRaLVZYLG3eTUQCCWN1ibsqSvAhB4naA5Eooev84mEa6zWejKU\nNFqb+eIp4SYSrvOxDl3jzifUwjqs9qj69gbZhPrkMPadVst9zU0dKPt+fFQcJvYcB4POgJSYpE4V\n1IAr3WFiz3Eh75REHP2EKxK1qwQ1QRDyhE1Y/128Qral3PkDz9LsLwoFep0eY9wVh/4lyJkdYxyB\n8wf57k5FEARBEOEmLLXBLbZ23tcr5KrBl2iuZRtKbhjhqXAzKDUbmyu2YdaA07ssf44gCIIg1AiL\nsFaqDjYmM/CauqGiZ0J3nDtQvuQeQRAEQUQCYVElHTLC+vR+MzolDYAgCIIgjjXCI6ydDq/3lLo3\nEQRBEAQhpkuEtVwnGYIgCIIg5AmLz9rOeoT1/SfcjuyUrHAMSxAEQRDHBGHRrJ1uzTondQAJaoIg\nCILwk7AIa7tbWIe6RCJBEARBHA+ENRpczyi3xCMIgiAIQp6wBpjpdFR0hCAIgiD8JaxmcNKsCYIg\nCMJ/wmQGJ2FNEARBEIESVjO4nmpvEwRBEITfhEV6djg6AJBmTRAEQRCBEBZh3dDeBABod1jDMRxB\nEARBHFOERVgbdK5CaakxyeEYjiAIgiCOKcLqs47Rx4RjOIIgCII4pghr6pZBRz5rgiAIgvAXSt0i\nCIIgiAgnTJq1HQCgJ82aIAiCIPwmvGZwJiwdOQmCIAjimCK85UapNjhBEARB+E1YpKe55gAAgAnP\ncARBEARxTBGQXdpkMkUB+ARAPwAxAJ41m82/KB2fFJ0AAEiPTQ1kOIIgCII4rglU1b0KQLXZbJ4G\n4EwAb6kdbLG7KpclRMUHOBxBEARBHL8EGvG1AMCP7tc6AHa1g7eV7wIAROujAxyOIAiCII5fAhLW\nZrO5FQBMJlMSXIJ7tpbv6ajrFkEQBEH4TcC5VCaTqS+A/wF422w2f6d2bEJ0POKj4mA0JgU6HKEB\nur6dD13j8EDXufOha3x0EWiAWXcASwH812w25/k6nmVZRDPRqK5uDmQ4QgNGYxJd306GrnF4oOvc\n+dA17nxCvRkKVLN+FEAKgCdMJtMT7vfOMpvN7XIHO5wOKjVKEARBEAESqM/6LgB3aT3e7rRTEw+C\nIAiCCJAwNfJwkmZNEARBEAEStvBsg47qghMEQRBEIIRNWOspbYsgCIIgAiJsErTNbgnXUARBEARx\nTBE2YT0gpX+4hiIIgiCIY4qwCetoXVS4hiIIgiCIY4rwCWuqC04QBEEQARE2YR2lJ82aIAiCIAIh\nbMJaByZcQxEEQRDEMUXYhHVqTEq4hiIIgiCIY4qwCevkGOrwQhAEQRCBEMZocAowIwiCIIhACJuw\njqFocIIgCIIICErdIgiCIIgIh4Q1QRAEQUQ44cuzpq5bBEEQBBEQ4cuzpq5bBEEQBBEQYZGgud0G\nhGMYgiAIgjgmCYuw1jFUvYwgCIIgAiVMwppM4ARBEAQRKCSsCYIgCCLCIWFNEARBEBEO+awJgiAI\nIsIhzZogCIIgIhwS1gRBEAQR4YRHWOtIWBMEQRBEoAQtRU0m00kmkylPdRDSrAmCIAgiYIIq2G0y\nmR4E8G8ALWrHkbAmCIIgiMAJVooWArgIgGq4NwlrgiAIggicoKSo2Wz+HwC770EodYsgCIIgAiUs\nfSuzUnvDaEwKx1DHNXSNOx+6xuGBrnPnQ9f46CIswvq8waehuro5HEMdtxiNSXSNOxm6xuGBrnPn\nQ9e48wn1ZihUzmQ2ROchCIIgCEJC0Jq12Ww+BGBS8FMhCIIgCEIOCtMmCIIgiAiHhDVBEARBRDgk\nrAmCIAgiwiFhTRAEQRARDglrgiAIgohwSFgTBEEQRIRDwpogCIIgIhwS1gRBEAQR4ZCwJgiCIIgI\nh4Q1QRAEQUQ4JKwJgiAIIsIJi7C22R149ovN+G39YVg7HOEYkiAIgiCOGcIirB95ew0OljXhx+UH\ncOsrK8IxJEEQBEEcM4RFWJuL68MxDEEQBEEck5DPmiAIgiAiHBLWBEEQBBHhkLAmCIIgiAiHhDVB\nEARBRDgkrAmCIAgiwiFhTRAEQRARTpcI6+a2jq4YliAIgiCOSrpEWFttVMWMIAiCILTSJcLa6WS7\nYliCIAiCOCoJq7DunZEAAHCQsCYIgiAIzYRVWJuyUgGQsCYIgiAIfzCEa6DTxvUFC5eQdjhIWBME\nQRCEVsKmWV86YyAMOtdwpFkTBEEQhHYC0qxNJpMOwDsARgKwArjBbDYfUPuOTsdAp2MAUIAZQRAE\nQfhDoJr1BQCizWbzJAAPA3jZ1xcYAHaHEwDw0ZI9AQ5LEARBEMcfgQrryQD+AACz2bwBwDhfX2AY\nBvuPNAIAquotvOAmCIIgCEKdQAPMkgE0Cf52mEwmndlsVpTARmMSehkTUVTu+tqfm4/gmnOGBTg8\nIYfRmNTVUzjmoWscHug6dz50jY8uAhXWTQCEv7SqoAaA6upmWK02/u/dB2pQXd0c4PCEFKMxia5n\nJ0PXODzQde586Bp3PqHeDAVqBl8D4GwAMJlMEwDs8PcE/XskBzg0QRAEQRxfBKpZLwRwmslkWuP+\n+1otX6IYcIIgCILwn4CEtdlsZgHcqvX4meP6ul4IpTUTyMgEQRAEcfwRlqIo91wxFgDQv6fHhk+5\n1gRBEAShjbDWBj+V07ABlFS1hHNogiAIgjhqCauwjonS869TE2PCOTRBEARBHLWEvZ/1+VOyAQAO\nJxVFIQiCIAgthF1Yn3JCHwDAnkP14R6aIAiCII5Kwi6s9e5mHi0WG1ZtLwv38ARBEARx1BF2YR0X\n48kW+/T3feEeniAIgiCOOsIurAmCIAiC8A8S1gRBEAQR4ZCwJgiCIIgIp0uE9RWn5gAAUhKju2J4\ngiAIgjiq6BJhfdq4vtAxDBLjorpieIIgCII4qugyM7iTZVFa3Yr6ZmtXTYEgCIIgjgq63Gd9uIIa\noBMEQRCEGl0urKlVJkEQBEGo0+XCmmWpVSZBEARBqNHlwpphSLUmCIIgCDW6XFjHCtpmEgRBEATh\nTZcLaweZwQmCIAhClS4T1hefPAAA4HBQX2uCIAiCUKPLhLVe5xr6tQU7umoKsFjtmP3hevy+/nCX\nzYEgCIIgfNGFwrrrA8tKa1pRXtuGBcsPdPVUCIIgCEKRLhPWNoH5u6vSt4TbhfvfWdMlcyAIgiAI\nX3SZsO6wOfjXpdWtXTIH4RahrsmKzfuqumQeBEEQBKFGRGjWBkPXTEMa3GanYDeCIAgiAuk6YW3z\nCMauEpJOp9j8Xt9CTUUIgiCIyCNoYW0ymS40mUxf+/u9pARPL2uHI7w+67Z2G+Z9sxW/bSgWvb8g\njwLNCIIgiMgjKGFtMpleB/A8AmjHMXNsb/61Lcya9fbCWuwrbsDuorqwjksQBEEQgRCsZr0GwK0I\nQFgnxEbhnEn9AAC7DtYGOQ3/OERtOQmCIIijCIOWg0wm0/UA7pa8fY3ZbP7BZDJN13IOozHJ672q\nBpeP+Oc1h3DjRaO0nCYk/LW5RPb9aWN6y87zaOFonvvRAl3j8EDXufOha3x0oUlYm83mjwF8HMxA\n1dXe2my71ab6ebiJ1jMRMY9AMBqTjtq5Hy3QNQ4PdJ07H7rGnU+oN0Nd2sgjSt/lfUQAAOMGZwLw\njg4nCIIgiEggFNKShbi+iGb0+q4vOQoA6UkxAAAHCWuCIAgiAglaWJvN5hVms/nKQL6bkhAT7PAh\noaaxHQCwIr8M7R32Lp4NQRAEQYjpUjv0eVP686+7snrYgbJG/vWeQ/VdNg+CIAiCkKNLhXVCbBT/\nemeY0rda221e73VPi+dfk9+aIAiCiDQiI8ILwPuLd4dlHDnNeXh2eljGJgiCIIhA0JS6FQ467OEx\ngy9eXcS/vvLUHPTNTESvjAT8b+XBsIxPEARBEP7S5Zp1WlJ4g8xy+6QAAP59ei5OHdcXpqw0GAQp\nZJFmBP8hrxB3vr6KOoIRBEEcx3S5sL7p3KEAgLG5xrCMp3cL5kG9Uzzv6SIjhUyOPzYUo8ViQ0Mz\ndQQjCII4XulyYd2jWwIAQBcmgckFkAkFtHBslo003dqFM0LnRRAEQXQ+XS6sDe7CKI4wmXm5wid6\ngek7kjVrDgpSJwiCOH6JAGHtmkK42mRymrVQm2YYoWYdlmn4zZ8bi30fRBCdRFu7Da8t2C6qSUAQ\nRPjocmHN1Qe3hykanNesGXlt2smycDpZ2OyOsMxHKyvyy7p6CkQE4GRZtLWHv8reP1uOYMeBWrzy\nfX7YxyYIIgKEtU7HQMcwsDvCo9Jyvl8lH3l5bRue/GQjbp6/QvBeK+qa2sMyPzWKK6lLzvHO6wt2\n4PbXVqLF4l3cxxcsy+JwRXNAmQWWDtfm1RGm55QgCDERkWet1zNwOMOkWbsXKiU/tUHHoLSmVfTe\n7A83AAA+eXhm507OB0eqW5DVnXrQHs9wlf4q69uQGJfi42gxG/dW4f2fd2NwVioGZ6Xh3Mn9RS4g\nKcvzS/HFH2Z0T4vD8AHdAABRhi7f3xPEcUlEPHkGPRO2HbtDxmctJFy+8+LKZqza4Z9p+5c1hzpn\nMsRRR1NLh9/f2XvYVb1vX3EDFq0uUqyDv/dwPV5bsB1f/GEGAFTWW2Bzu6kMJKwJokuIiCdPr9PB\nHqZwZ7nULSHtHR5ftVwaV02jBe8s3ImmNv8XSyFzPt2ET3/bh8ZW7echAyTB8cVSs9/fkVqvLFZ5\n3/cHP+/GjgPiWv0rt7s2ltEkrAmiS4iIJ0+vZ8KXusX6ENaCBUxOOL7+4w5sNlfj59VFYFkWc7/e\nil/WFMkc6cHJsvh9w2FUNVi8PlNaMOWYOaa35mOJY5vGADRrab/2Mom7hz+3ygYyyqD3e1yCIIIn\nIoS1Qec7wKy9w44f8gpR2xhcoJdc6paQhhar6Fipds1VEmPAoMPmhLmkAQtXqQvr/P01WJB3AM9/\nsdnrs50HtHcbS4jzdCmL1OItROQifcYWrVa/b+Uoq2lFUXlTqKZEEIRGIkJYu8zg6pr10o0l+GND\nMR54d21QwWhKwnpMTgYAYLfAj8eyrFfedWsAaTOrd5QDAJravCN4pdqOEKlA5o7N21aK6+fmobKu\nze+5EMcnH/y8G5v3Vfk8Toul55nPvTedBEF0LpEhrPWMzz7Slg7PIvLHhsALhNidLBgG0EmiYK+f\nNcTr2BX5ZYplPp0sq6kEaEFJA/ILaxQ//yGv0HuODifKa1u9Umy4a/Tlny5/5SYNiy9BAMD6PZWa\njvtt/eFOngnRWYTC2lZU3oTbXl2JJz/ZiCNVLSGYFREqIkNY63xHgze1erTSqnpv369WnE4Wep33\nP1vOF/fN3/sVNV+7w6np4agIQPu96aXlmP3hBnz9137R+9K5kCk8OBparFiy7lDEFcBRojnIoEZf\n1DRY+NQw4ujC7nDi+rl5uO7FZUGd58Nf9sBitaOkqgU/LPdWJIiuI0KEtU7VHAwA63ZX8K9X7SgP\n2HftcLKywWUGPQM5L/baneWy57HZnZrqdQvPmSjwOWuBi8Dl+HXdIZTXeoKCSFT7h83uwBs/7sAu\nt0B6b9Eu/LTiIJZuKunimWmjMyuXddgcePC9dSiuJG3qaOTxjzbwr60dgW8+he7BQNIDic4jMoR1\nAEVRPl6yJ6CxHA5WNriMYRjZgg9fLi2QPU+b1e7TdO86sedlj27xXh+nJ2vv550UF8UXaAEit455\npLLnUD3yC2vwyg/bAQAFR1x1rmuCDFoMF8J5Th/dC6U1rbDaAl+Yhb3kP/19X1BzI7qWSoG1cbWC\ngqEF4dJY29TOx9sQXU9ECGudjvGpWUvZV9wQ0FhOVl6zBnxXZ9pzqI5/veNArchn3douX/5RqK9H\nCTp9cQ1MUhKifc45IyUWADC0f7rofTKD+0d8rKdgX4dAyOXvV44piCS+X+Zxi1TUteHxjzbgqz/9\nz7fmSI733HsbFHzaU0f2lH0/kJKlRHiwBdFnQZgx0Npuxye/7cXhCipzHAlEhLCub7KCZYGFKw/K\nft7kR+EQXyiZwQHf1Zmku0yhrGyWifQGAGEcGydc87Ye4Rc7LZXbuHae0mNJVvuH0HUi3GiFq5d6\nsByp9rhAuM3qml0VSof7hNsgq236oqPk86r9qQ9AdD4DeyXzr0urA3dlJMsoD/WCdFai64gIYV3r\nbpLxy9pDsp83hPBmcTqdiouzr0IT0h2rMHVKSdMQCmuuSpvQtO7QIHG5nt/SUqjHk6y2WO0ibTgQ\nPvjF4zoRWnKyMhODOm+4GJyVGtLzcRuWgyp501xdAQC49YLh/Gt/LWHHChar3a+qg6Fi5fYyPPXp\nJkULnrACZDAbuGqZwk3UvCUyiAhh7Qu9PnTTVNOsfdEhEdbzvt3Gv1625Yjsd4RmcIeDxYe/iH3t\nWvze3DmkAWfHkxn8tldX4paXV6BNYbHyF6Gw2e5HYZquJFDXjxIt7ujyhmZl4aPXe+7f8YMzMXFY\nDwDH7wL+8vf5ePi9dWEf97Pf9+FwZTP2KtRzD9XvIYxj4Phny9ERgHmsc3QI6xCaKYMR1mopPss1\n9Juub24XRbVz8/GF0m76eOTnEDUzCbew0ZKTHy4+fHA6AFeRno+X7MHbC3eKPr/45AGIiXaZvxmG\nwfWzhuDx/4wD4BHe4eqSF2kcLGuC1ebost+zqELeChKq+UTJKEbCTaKTZfHln2asCSKITQrLshQD\noQG/hbXJZEoxmUy/mEym5SaTaa3JZJrQGRPrLBpbOgL2UQai2QijdeWeJ6Fm7aqY5n2Q0oPIdVE6\n1hE+yKFyidz39pqQnEcLb/60A7cI+qMHSkpiNDLT4oI+j7DOwJqd4s1jckI0Jg3vCb3bf2PtcGDy\niJ7I7unyiRp0nLCOnM1HV6ApE6QT+H29fEGoUPVWiIvx7po8amA3/nXhkUbkbSvFx0v2hmQ8AHhn\n0S7c9NLyoFLOjgcC0azvAfCX2WyeDuAaAG+HckJyKD0Yf270r5IZl6da3xy+gIlVgqA0OV8XF5g2\n9+utuH5uHkqrvZsr9Ez3TvkCXLv84wGhr9pXDXkl/HEZWKx2r57mwbBtf43mIjpq2O1OWc3HF8LN\n3vM3qe+tX7tjCtKSYnDjuUMRH2PA5acMEn3OCfpAf4djha7crMgF3AY7n7KaVjz6wXrsLfZWALq5\ns1EA8RoWKm1+i7kaQGhjk45FAhHWrwL4wP06CkDg5cQ0wi1yMyRdp5ZvK/XrPB1uM/aQfmmyn186\nY2AAs1OmocXqs+mB1eZAdYMF5hKX1r5uj3dwSL8eSUiI9d7xHgs4nE4UVzarCjJh8Eyg5rJmi7Ir\nQRq49f7Pu/H4RxtQIxNsEwzBLm12J8un/HH4SjcEXEIeAIZnp6OHe+PHacpCzjwpi389alAG3rpn\nGrqnSTaKbqNUsA11jkaE7qiu0qwByAZaBlsr4MflB1BR1yar3S7b6llnhc9puzV4TVi48WCOkqyM\nrkJVAphMpusB3C15+xqz2bzFZDL1APAlgLu0DGQ0Jil+9vB/xuPFzzchPtYge1y9xaURJ8SL0woY\nHaN6XilOvcsPl5YSJ/u9hHjtBUrkkJ7TJlsTzRuHwCwZE+Nd5SwxIQbGtHi0ygh+4Zj+XItI4eOf\nd2HRigO454qxmDmur+wx+ibPQuTvbw64FphaGY2BP79BLzon18tZHxPlNVYw1zgjIymo+AuHw4lY\nmU2brzm1uDcqCfHR/LF3XzEWd72yXHRcQ2uHz3PFxrruz+TkWJ/HOpwsNu6uwPih3b02Gb6IxHvZ\nIagemJaegKR43zUSQjKuZGOQ3i0RRoG1TS7oUsv1Ex6j1r8AABra7RjQKwXvLd7NvxcTHw1jtwSf\n46gx+8O/PfOR/LsIMarC2mw2fwzgY+n7JpNpBIBvAdxnNptXaRmoulo5sT63ZxJ6ZySgtKYVVVVN\nYCRNNhYvdxWDaGxux5knZuEPt/nbbneqnldKpfths9vsst9bsVU+olsLGSmxXuesq/ekdg3qk4JC\nd8UsjtGDMpBfWINawSJQLCOQrVY7Dilo6NyYRmOSX9ciUli04gAAYOOuMozoJ5+aVCcQ1i2tHX7/\nOzfvq8I7i3Ypfm61eu4H4cLX0NCGA6wTSzeW4Lf1h7H4pfNQW6ucw2rtcCAqSufVJIajqqrJb6HF\n4QrCYcHK1If3dT047cXh8DwviVHec4w16HyeKznOtWTUN7TJHrvFXIU1Oytw3awh2LS3El8uLcC0\nUT1xzVmuRjkWq13WLyokUu/lKsFzWlXVjHYNBY1CgTSnvba2BTqHR6sVPh857nXG1/WTXuPoKB06\nbMpWq3tfW4nURPG/t6a2BfogAw2F5ZNr61rAOI4dv3WoN5yBBJgNBbAAwBVms/nPUE2E8xFu3Ovd\nSYrLbx6dk4HLZnp8aP6afrgIYKUFU/hQZHX3L/dWzn8jfCtGprhEb6NrV9ohiDKX86dLTY6mvqHN\nt+0qiis9i4Wa+0tocowLwB2wzUeFMqHmYhUsWNsLa3D3G6v5TlRb9il3rrLZnbj1lRV46ZttiscE\ngpNlkV9Ywz8fUsVci9uQcx1ECdKwpBtiQOybVELvI8Ds7YW7kF9Yg32H63nXDheYebiiGbe9ulKx\n+FGkI/w3h9Nn3S4xTUtN8E2CBi86hgEL/9M6hYczDPDBA9MxqHeK6JgGSR2KULsCjvOYRZ8Ess1/\nHkA0gDdMJlOeyWRaGMoJFZY2er3HBVL16y7eqfjrx213+3qUqjIJ75UnrxmP3L6pSEn03j13l4nI\nrWtyCVlrhwO7DtZ6pSNIo3jH5GTwmwbhjraq3rtLV3FVM8aZjPzfJ4/phcy0OK+d7tGGsIyh2toi\nLBwjXUC0INNkTYSomplAhm12B75wrFWpk8y1cOUElBwlAbQc3Li3Em/8uAPz3JsAqRldy4LJFdPx\npdVPHdnL57m4TArpuHaHE8sElikWno03t0ne7ja1KhU/inSE6X4lVeHT/KU+amncBufm6JYc4/l9\n/BbWwqwU171y5Wk5qt8JtXCVWo0IMX4La7PZfIHZbB5gNptnuP+7MJQTklt8uDKgUh+Rvzdkq/um\nVup+JVwHGYbB/f8ajRdumiASlAD4HFQpdocTSzcV45UftmP1znLRw60T+K/n3jIRd1w8kl94hZp1\nq0xnpdgoPW4531M9imVdi3ZXBrmEgoWrhBqW8r9FnN7m/zi+/MSHK5p5641DJZjt703K2QcVtfKt\nUIVzbwggC2HLPteGgVuQB/QSb1a0XA4uwExaTveDB6aL/pYrNSmFS+mSapbf/L0fXyk0veEqFIb6\nbmVZNqxlT4X/5tcW7AjbuFLNWkmrnza6N7+G+arGKEX4XHGFUXxlHoRi/ekmaGQkt/YRHiKuKAqn\nyTqdLN78aQfW7a6A08kiIdbgFfnqb/pInjt6vL5JftE840RXNCwXFWvQ6xAbbcAtgjKLAJCSIB+I\n9uLXW/HXZpd2cbCsSbTYW90CWa9jYEyN488PaIisZVwazfABrkYevTMSXD3Aj3JhLTSr7VHJGXcG\nGQ2u86VaA/jHXYFOOJY/uat/bFDIfxWcry0AwbKlQKzdCyuKaYV7TqSLr0Gv09RIRginuUnzbKWZ\nGeEomvL1XwW47dWVIr9nZyK89+QqfXUW0s5q3n3tXf/XMZ4qi+UKm0clhIoPd5/0NibypY45uiV7\nXCWhENaJcZ77L1y/49FKxAhrLn2GE2Tlta3Ytr8GH/6yB+W1bbIpKja7E0vWHdJ0frvDyUf5bi7w\n9osDwLRRvfDqHVNw2QxxbqlO0D6TYYCdB+XLUx4sa+I1IJ2O4R/u3D4ebUjoFzxY5jL5a63KdcdF\nI/DMDSchq3tSQJ3KIhk1TSBYX6GWAOyGFitYlsUTn2zk36uV2dQpLSgKMWW86RdwdcoKFr2OwYs3\ni3Olffkn1czgL/13kn/jCxZvte5OchXiQp2Yw6UUqbkeQsmLX2/lX/fKCC4K2h+q6sUphNJrywlN\nHcNg1KAM1zF+bpaEt5DwN54uSZetbZJvhhMowuc5NYwboEDYYq5GVYjTOf0hYoT1BHfNYW7hqZTc\noNLgBo6fVhzUpG0Ju2IpResCyi0ruYWJ20wIkVsEWSfL5wcPzU7HOZP6Y2DvZNwwayh/jNQnqgRX\nGzzKoEdv9yLR1NrhZR47VhEuCjsP1qKqwYJ9Gqu3FZU3ibSM8yb3lz3O4WTBwtvkKEVuU1HX1M4H\nsUkDE4VR6EvWHQ66CYRep0OmJPfZ1waGC1pkZQzR/kanC58doRtDakKXzsnhdAYlrZvbOrBmZ7n8\nsx6GPas0kyMtMXyChbvk3OWTCmLu8WAYxhMAGETRGqEFRm2tDIWwFm40g2nt2dnUNFjw9sKdeEdS\nmjecRIywlkaZqi1Aj/3fONHfN720nDdjKiGs6x3MTVFVb/Fa8K841TsQo7XdzpvbDXodMlPjMPvq\ncSkzXEYAACAASURBVBgk0LJnTeynOlZSvMu3LrfIcpuXUFT9KSxtRFkIK3aFGuHapGcYPPzeOsz7\ndpuXeVCOZz7fLCrLet7kbP61MHbB6ZQv9SpF2swFEFtGfCk0cz7dqH6AD+TM4L4WZm7BZUKs23J9\nwEuqWryqakldCGU1bWi1BO6T/Paf/fh4yV7ZvtvhsC89/9UW0d+rQ1gb2xecZaRHN9cmzS5ZGzmh\nqWM8m69grG7C2Aa10swrNfRD8IVwnnlb/StyFU44f3pxZeDtR4MlYoQ1t6A4nSysHQ78seGw4rED\nenlXX/r6L/ngFg6bYEELJKJYSGy0OAo9R+Z80QYdv7AqBcH4SsHiWzeqPHdywsMfnE4Wz3+5BY99\ntCGo83QmQt+YcA3SIqylCBcfYcW6iro2TcFrchu9Tfs8bhVfAr+xpSOozmEGmcVT2jpVCqeJpScH\nrw0KLU+ua8biw192ex0nFShPfrIRf20OvHvTJndkuawrIQgNz2pz4OnPNuG1BdsDPkdnw91zse7A\nVulGiLvnhJq1P7Ed0ntWeI+padarVLIjlPhjQzHWuVt4trXbRL/nzoO1EVtCWU5hCjeRI6zdN0hB\nSSMWLC9EUXloUyPsgkV2+hjfKSpqSNdLuUWww+7kHxivko1u1EyQl0xXL306LNsVbBbsD9jcFv7e\nvEoo1WwXmv2Erzt8mKx9Cc4kQXBLQlyUNmEtU7QhXlDkQ4tGc9/bawMOwJLTdHwFwnGadyi610nz\ns39dewgWmbKToe5q5uQFkvdn0pEOVTRpthQtXl2EQxXN2HGg1q/cZH8D84LBzgtr130mvLaf/7EP\nS9a5FBudjvGZBy+H9NgCgclfQ2ymX/yQV4gPf3W1CZZL043UDoOR0DQvYoQ1F2VY02jRlI86YVh3\nv86/VRBVO3Jghn+TE5CZFocESepXrExFpk37qgRFWOQXSTWt2JSVygeLjMk1en1eWu26RnU+0oFY\nlpWtJcyhppU1tnbg8z/2iSokaaG+2Yrvl+1Ho4qJXs58rzRP4Voi1LK3+SiRKF2E5lw7XvS3weD5\nXfQ6RtNiLddoRSgutPjxrDYHCkq8FyotcI00hBs9tawIm93JL45Km8NuybHorrHMo/QaRRn0sgLU\nl7bvL56IZ+/BpKbJpz/brNlSJIzi90cblaut7ovCI4349Le9fmc0dEg1a8F9vSK/DIfc9QoYxuMm\nabHYVJ97IWrzUdOs/cXr3pG5H0M5XighYS1A6MvV0iptcJZ8Mw4l1u7ybpARCHddMhLjTJmiwihy\nN1hSfBT/ECgtkn2MnohSafGVgb1ScMoJffDM9SeKGixwcD7rFT78RvO/y8ctL6+Q7dQDiB986THf\nL9uPFfll+PwPs+oYUhavPog/N5Zg8eoixWPWyfweSs+pUyEa3NeiJwwqBIAsd1Gd+/41GhOGdceQ\nfmm46rRc13kdTr5ojhqcFqOENHJXKWr4pW+1VTob2l98n3OL8et3TuEFRpPEOtLc1oECd4S00F+v\npFnPvXUinrvxJE3zkSK1Kl19uut6ylUp46r4KW1etSBXeW3l9uB9p4C4ep0vAkkhfP6rLVi1o5zP\nStEKZwbnyrTa3VaZed9sFR2nYxh+rfl+WSHueF1TJWivrIcTBMpBKIWntJwpp/QYUz0ZMoG2L+5s\nIqEffcQIa479Rxo1+WGnjOyJGWM9aQU9u6lrBmPdN6C0Cpq/6HUMYqL1eOHmiarHTR7Rk9d4lHJj\nE2I9GvpN5wz1+pxhGPQ2Jqo+MD1k/t1OQV9sbrE+pNC0Xs1c2eIWdv6YpliWxcrtLl+WUqoSy7JY\nsPyAzLzlzykU0EJTud3HfSLnSwWAYf3TcdO5w6DX6XgLjcPJ8i1U/Uf8+wj/3fGxBjAAzp3UP6Az\n90wXC3tO4MbFGNC/h+telua7frxkL178equ7RKnvPG8dw2helKVWJb2OEW2y1FKauHsyGBO5r1kG\n04ZUSROV2wwE2v0tkO96+6xd/0aujCsHw4hjKmxuV5yvgNrPfhfnzAvTS0Op6LZ3iO8/TgCOyRFu\nDkI3XigRuprMKk2BOpOIE9aAtnxUHcNg8vCe/N9yO26O2sZ29Eh3acKnnyjf2UkrwnHSkmIUqz45\nnSzvl1TSrIWV0KIMegxza1FayqhecYorAr2HxHzZYrHhhrl5mP9dvliwKSyQDqeyCVfNT6jEYUG9\nb6We0DsP1vmci5C65nbZY9bLRAYLkS5mcnBCyuFkQ5a3/vJ3Hq2ZdbLQ6Rh0T/dO+dOyW3dKPLLC\nhZ7LFpAWlOE0N6mrIRSRrH2MiaKKfja7E9UNnt9Hr+Lk5IQGi8ALavi6F4XXp6CkAT8sK9SsFSkF\nLH72+z7R3zqGQWW9Bc98vpmvq9CZcC4F3mftlHdt6RhGVMIXAB56bx1unr9c9fyVdWJrUJTGaHB/\nET5fLMvycQXC3zRSNWthf4G532wLS9EfKREprLUiDPJQ0v7qmtrxwLtr8aW7FKKvEnpKcJq5cMyX\nbp2EV26bLHv8ofImvpqZXASvlCiDDndfNgov3ToJb949zefxXHqFdNHjolr3Hq4XbXqUdvOi4C2V\nYgtaEabnSM3QG/dWoriy2Sv4J8ftAlGqDbytQN437W+VJjk4k6zDyYasfGttk5XPp3ayLmF94hDv\nGItWHwu9k2W90lmE+b1cIxulxhhL1h7yatAQCoQV/KQWErUKa8KrG+jGyFcevNAqN/+7fPyxsRgH\nZAKZZL8rYwaXE8YGA4P6Zlev+js1mpqDweYWzMJo8G//2e91HMMwXtdfKWhTiLDcMSB2l/hbtlQN\nu8Sd9dGvLo0+f38Nb3nS8gweqW5BtY/iJCzLitJ1g2XpJnEmg7/VM0PBUS2shf4ypZtKWkxFrhKa\nFm67cDg+eGC6qAmITsco7gQLjjTyPmClxiFCDAYd9Dqdps5HgOeBkt7cSqkPSqYwoYB2sCwWry7C\nS99uA8uyvGbqzxovffA5rDYH3lu8G3M+3eS1cejHmXMVNKBRg7r5MQMPwpgAJThN0OFwqu6W1e8b\n73l/8LPLBO90in2JQuRcAUL2yhR+Ed5vvoTvvuIGfPGnJ95g+ujeKkdrR6ixS4WB1C+uNMdANZPf\nNxSr+qiFFeO4+2xBnuc6t1hs/DFSk7mcZv2ROziP4/aLRsAQ6hBpH3ABZJzFze5kYZaxGul0UAwU\nVHMP2CSbFOG93r+n2G04ZUTPgBsICU3JwjWg3eZQXM/keOLjjXjovXWyn1U3WOBkWbz+4w7cPH8F\n2tptQblGlAh1toMWjmphrWb6DuV3uO+ppVqNlYnY5pDe8HL4u4ng/hlqJr5KTZq15/vWDgcWry7C\n3sP1Io2s4Ij2yGWlOufCzYJ0LkJTtOwcA3wwTH09wVmcyVgKJ/z2FTeo+vYy3ZXr5Jq4pCV5b7A4\nq4bDbQaXY6MPM75cUKDw/hW2XeWyA6QIFz9/274qovII6SXPiLCWtJBgXA6f/b5P8b5vkSm8UtPo\n0cLmf7sNr/+4A3sP13v93nIZCtJgsLG5RlGA3Jknegd/+sLfNSjeLaS5cpwOB+sVyAi4IpZPHCyf\nJaOmCUo/Ea5zQ/un86/Tk2Nw3awhXg2VtCLMpBBee4NA6XEEIVh3H6rDQ++tw3NfbOF/t9tfW4V3\nVXrZB4r9eDeDyy0mwrragSB9qIWRh6Hk9otG4O17puH6WUMAeBb1wVmpqn48Dn/N874EHACRVqVk\nPhTudp8U1MUONPpRbiGqrGvDakEBhWaBafGFmyfwu+pnPt8se85AHwxh+hD3u6ghbUbBMWNsbzxw\nxRhk9UiSdWlER3n/dtz1Y1mWD5qRdnsLpKCNUHMV3lePf+z67aSR4cLfMdCNqhS1s6RL6jvLlecF\ngu8HLdWu/3YXXJHT2BtaOrB+tyv7oNidFrp6RxlueXmF6DjOLOsLYYUvrcUyhMGLcpkQajicLAx6\nnccK5HTKPp9rdpYjyiD/6/jjWxcKa+H9xg15pUzFRi0IM3KEmwdhMRdfmrVacN4a9xpTVC62Lm42\nV+PFr7fi3UW7FJ9xX0gLWB33mrU0krRbciymuyO+sxW002dvUE85kf74MRpM0oESF2PApOGuGudc\n+pmWICfXvPz7KfibW6NQ/X5Zoez7SotmoP5bufk88sF6/JDnGZ9rFRkbrUf3tHiR5nlERkNU0grk\nhKT4e+4Ul1smasqtV9IYrj7dhOSEaOgYRjZiXe5acZH0nM8aAP51yiCv4/xFJxLW4oXZZnfi7jdW\ni97jNgQXTs1GqMhIkRfAgOs31RKgFMxil5oY7WXB+eZvlw9X6V754Jc9orz/dbu9rRrc/eLrmRKa\nwe12bf8OYfW2fB/1AeTmpdd7/NEOJytriatuaFdUDPwxBQstB0LBPXWkK6DXlJWGXhkJiq2GlcgU\npLsK/ckWq92jWftYd3YVyQenAuoBpwUlDdi0r0qkwPiDNLX2uA8wkz7Wj18zDuNMmbhgSjZukElt\nAsQNBLjdsxDvIg6d+092mcu1azC3XjAcF07NRnysfze+zg8fD0eBTHciaVlIjkA1a2mQmNwiwWl4\nF04b4P7b85lcVSOlCl0dNqdqEAm3+EpNs0pwU73zkpF46MoxiDLoRIVUdDpG9rrYHaxXEB636DS0\ndPCBS/5qtnJahHAcRvLPkqvaxW0a+xhDZAKHJ9BJDoZhRCZOhpGv9uXvZlD4TGX3TMYwgXmWY+fB\nWtVqblosGfXNVtwwNw+/rPGuETAmx7XhE2rWchXt5MjbKu5d4E9ZTYeThUHH8FYdh4PFtFHeVRhr\nm9oVN0dq15vreMghFPjC637+lGzBMQxsDqemNLSqBguWbysVaafC8s+XCroc+rIA+NO2Vom2dpti\npooShyRR9tIMgXAQUcI6TlIJLDk+Gga9DudNyUbPbvLBQkLz8Se/eV9A6U3qb5ehQNAqHABg/OBM\nnDvZf62HW7T9aeTx+o/e9Y+VNJxAN45Sn5OczOeqyXGaofAYYTAQf06VhWZ7oXKBiY3uetJaN0+e\nhggMTFlpeP/+6XwhFW6+cotee4dDvgwmy8JitfOBS/5aoX+WKSojXIylG4QWlXz4QPpgK3HKuD6a\nLVQ6hsHlMhYFpbx/OViWFWnM2/bXYJ5MUZkFeYWqxY98XYFzJvXj298uXOV97f/lTpcUukK0NgVq\nkmRGPPuFvMtHjtLqVrS223khanc6kaSi1Y4a6B2Qqbb5ln4kvK2Ea5lws1le2wZrhwO3v7bS1/Tx\n8Hvr8MWfZvy+UVAtTnDdJgztzjdG+VRmDRfSqlILQUtAKQA8+sF6PP7RBr/qRyRLrG57DoU/1zqi\nhHUg/X6FAT9yUYrStVUuQCjUaEnVChbuuRGWS/SlDcuZyA5Xyi+aja2BdfOSCjM1c5GcKd9itXtZ\nANR270pNUo4IStZq3aBxu/bkBOVgNDlhfaS6RXZDIf095FLg1Hb4cv20hcJ6cD9xdbNfVCrGaYmb\n0EpCbBTuvXyUpmN1DCMbj7HajyYQWvPDj1S38ulssvh4LKMNepQojHXKCX14/7tIsw6ikc4PCq4p\nIcL7jdtw7S6qk3UXPXTlGADyfX/UDBl2p1PkUhHep0ppm9wzKZfupoSwzaiwQE+UQYcy2TK+3qhZ\nCMYPztR0jia+2JP2IkjSgi5dQUQJ62B3Kx12J6w2B8pqWnHXG6uw51CdaMHM7ZsaltqzoWiY4AuT\nu9yq3cGist61yfFlWpQzMf26Vr585ldLPV3MlCJ65eCEFicgbSo+PU6AsJLn/TtJDqmaf1PJjP+E\nIFjOl2bNmTe5cpNK94jLZy0eb9UO5TQibiEf4haqcuetqPXPHCe8tYb1T8eEoZ7o39E5yhkJob4n\n5f4tXK91YYlUhgG2yuTJbz9Qi+teXMZrsmo89dmmIGbqwVeL0MLSRlFRH8AT4Hr5TI91QLjB0ppv\nK5eN8MfGYp8d2IRCjfsNiytbsNlc7XVsb7er4+wJ/bw+U7NO2e2sOKgshFaYqnp5BYxraXyCyQiG\nYXDSUG29HtQ27n4HLfrh6vOV3x8OIkpYjxgQWD4tp1E3tXbgvy+vwG/rD6O5zeZK8RD8gL18lCQN\nFcL0EWlt51AhDO74aYWrKIavph7+ICwkUetHIw/uenOxAWo+ZU7Zk+Zmnzxa7I+TCuSRAjOfr5Kj\ngG+tkluouHko+f30egYsK/bDC8120ohzLpCIuxbDstPRLTkW15w1mD/m7YXyaSVKAUHxEleRsMyu\nWsBdqGM15K7RXZeOBCAO4jRlpalGTCsFPmpFuglRazv7zd/qbXR3HKgVxUzUN1vBwrXhEAozYe+C\n/MIaxbr7QgYoNP4QVsaSQ+ij9bXh4n7j3L6p/AaRQ6ngEODSrA16hn+upF0CH/33CXjk32MVv68W\nvPaOQtoUlwrJCcGzJ7o2GL60Y2m5WyGhqkCo9dzhrhceUcJ6WICC7Y6LR/KvpZevM1JXfJEhSA+7\n+GT1VpehgNu9fve3d1UjKVpdDdLFXWs9Y+56R/PCWvl7nHlzmY+m83ygmHuxEgqDeA2lWX2VMORM\n6dyiq7Qo8j3XFR7SySN64u5LPeZhrpqbsJ73S/+dJBscJEVp4UlJFKdGCeM8hNaQE0xiLTvUsRry\nzWtcm+bROZ7I+1NP6IOpI5X/vf4ssNPHeBd1kW5Q1O5TX4JRyn1v/397Zx4mR1Xu/2/1MvuezGQm\ny0wyWU72hewLWYEEQTYhBEQiBhCJeEGUqygXFES4/kD00Xv1KhdBrl71et0XUNkuiCyBsIRQgRAC\nCYFMQjJJhsls3b8/qk71qVN7d1VPz8z7eR4eMt3VVdWnT533nPe87/d9AqlU2tIf5kw0Zxb4SQdy\nyh/2ijnh7TNr/DBLTI+M+MzKgjpuhmVvWwcS8RguPX0q/vnCOYZIEWfC6GpMHO08CXL7DWUVQ86f\n9f1rvhfM6yR4PavmVDJpyy1ghkEQW2vXfmEpHvqloIx1c5ZFNsSUAMD8I4hbpvmqvlYjRL/mI6CN\nG71tbzqnNXD8rAIA636O30G1T1pZu0XgOj1cota0dhyP6s78gB9Z0Wq6jsxZAVKVeLQyv1cnNzgP\nSNmz3+q65lsF4qqfbzvYGX+vYBi7SY5d4NCJNoZ/02lTLK5QP3rzQbAbVPn3XChIq8ZiCsY0uESi\nB0kpsrmmXEv7aGePMVH0y/Ufm+v4np2ojbxqddIUF+ErW1k34o//sN+G4vDJR6UebOuG2G/XzB1t\nes/JWHMvQXtHNypKk8b2mhejhP671WUS5JUytv+wtnhI+pjcy+/L++VB9RiCmFq7rzGkjbUcMOMX\nuRPzDnKg/bips+SrVqr4cOdSDtALOaDOT7CLHMXLZ9HDPWROD/t0sfMOzFcMTgFggBZ9C1hTe+QB\njE8AuBGYNq7OWMU5Gfxf20TzOsH3Bd/S9yudZvc79QCZb//yRct740dZ3Zy8RKjd+ZYIRWj22+gc\nv2ojNTrVJl3JLiq7rCRhmSD4lbH1i9134s+hPIFyM55Bhjs/E9/9hzoDu/zdBIlSKWtanpvil+N5\n9A9dvG6y6XWvffTd+kT88Zf2BSpysX7VBFy8lhm50U52jAsHyd4CL8StHCdXN+A9F+Nt69dYi54T\neQIStFaAH2PLr2F3bJtbMGMEFJSxBoDF0xoDf0Z+2ETXaH+4weMOCkBhs1F68Of5jIYU6etLo7Q4\n4ag0xfntE2/6O5/eqXngjZgeIZcx5S7d1Se4a1bzc31kxXjcdMl8nDizyWhXL/e8OKg4wd2jvKt4\nTersvBNuLla7PiD2UbsgIzuPxPhR/tT8NLUrs3hK2H3frVsrioIbPz4fd+hFbkRtfIthC2Ct5bxy\nQNufnd5qnsSI0doNHv0acHe99qWtxlrGScrWdB79uZCrr3nlfour1iCLjWQihpVzRhmT2lvue9ZW\nh4I/P0EnOHLshBPiYslOPZJ/o6QRkOrupXjnQMYgyyVtvbYJLPfm8f7ftuzBpbc/jH0HO2w9BE4S\nv1FRcMZ6bRYlLOWHTfzRxJma34jDXImbVtbRNbE4KKVSaVTrD+a1G2abjvvKJxagWXdFptJp7Hj7\nsLGK7NPTNrzqLftV7OGR3dwdzF1VVWVJnL86I1P4yTOm+TofkMnLTiZiaB5RadJpd4oG5/g1cCJO\ng3eTHu1st10jB+WI2AW4iSsZO++A3V5860j7ICVZSSoRUyKfMMrnvGHjPNPfLY2VqBWkR/k9coU/\nTq5KUPXVJfjs+tmmiWoiFsPV583CVR+Z4SsIyM1Y9/SkrNHR0in9RISn0mkoitXgek02eezB2MZK\n1228uQ61CcSu9x+/e8VyPR6gGXSc8iuP3CW4qm+/Yomln2TuU5Mc9VpZVwgTo2/9j9nDFYZgish/\n/UWLAfm/F/fZ9qMtNhH5UVJwxprPBL3csm784cmMG5XPaDesmeg42IXNcTHdIkJjLT74famUsWdT\nV1mM6y7Qci7XLhiDMQ0VmDpOW32k08Bt//UcbrrnGf1zacTjiucWhLyP7AQPpOHuWe4imzauzuRK\nEidOsnrbshn2LmJxQOHbC16DXUMWWvBOg/clp2sTDNHQjtDjJT53QWaCNEmKSLY7n1NgGEc0bjdd\nMh/XrHfOa5avF4/HzLKkEWzFyEZnnEO0c+Z4/f8xxYg3AOxzyZ2wcxnX6+0/TgiKSiRimDl+GOZM\nrLfsL8tcu2G262Rm/+FOy+8nD9uilK4TaT1Qzc7D4Tah4CveuazesV9+6eK52HzODNv35N/p8m88\ngqe3Z2Q5+WQ3sLH2uRIvLTZv01jaWjH3U6/JtxjVLueaexn6BVPMnke/e85/fuotIzZipZCp8syr\n+319PiwCWxLGWDlj7DeMsUcZY39hjHmHtgagtrIY/3zhHFynJ/jnCpeF8+u2CQOxSlWUbnCxROiz\n29/Do1vfMa45uaUW//mF1cZqls+gX5GC0Pr60r5EXPxql7fpkemy+EwyEcf01jqMa6rElWdNN723\nfNZIrF0wBp/SX39ScNe9LeS9ikaHT4K8IkCTCW8RnEaprKDTb8ZlNkWXWBqabrCocCSnD9mdTxx4\n5dxeIDPJrK4oQvOISte0RjkuIh5XEDcN0hEY66D9mt9POo3TFo/N6pqKYt7qWTClAafr5xKNh9ge\nXm7mqS21nt9FLiHJf99aqWiJG302e9/Gey59mGvMx2MxR9W4mnLn+7C75vd+s834d2ZlHez3TNgY\na01pLiW9Zj5Grv8gEo/FPJ9nt4mNl7GWvWxBA8TGNlZaYg7ySTbLvksBPKOq6goA9wO4Ltxb0nIz\n3YoFZINX0YcwEWfzURprsd/eKghH2M2S/6qLEPzvY2+YXtfc4N5t45bfKLJNF7aRH7pUWhNeuGHj\nfMveejIRw/mrJxqrVDHy3ElNydBJdnCjVpQmLXvkTsgFNpx+Mn55cf/ebrJzVKp8ZWcMvPYf1be0\ndpT35eyQf++U7i3huAX5ZUvQPXB+dC7xs2Mbq0zBauevnmi0bZGpXGgmWt9LUlJRFM+J/AdS+7U0\nVuKOzUtx/UXOUeQyqXQailNRE5etAF5Pm68iZV32K8+a7ho86HRNDvfGBV1Z2/XfO3/+Ai7/xiMm\ng8rzprmuuNt14jHFc1vEzVi7TczWr5pguXbQvGze19brWuZhZ1h4Xj/oB1RV/RaAW/U/WwDkXyQ1\nC4KWoMwFMU8xysIh3Q4pI0Fc733CwG4XjLV8luaSdsuztKNFKnfqRxzGTv7PabDh39Fpv1DLj/XX\nDjVS7rKTIeJGry+Vxp+e2o0D7Z2mqlqcR7Y6q5rZcfI8a5wGF7p51Me55Nvt7OoLvvLNAV+rsszC\nGkBGGlNOMZLhJT+bGyo0V7DwZcWVrdOqU/aaDKuyrkKDVo/i1x5WXYKaiiJfQWzvH+lyXlk7GI1n\nXt1vTG4O6dkYJ80zt5dnTIbHnj2vGhbUWJcUxY19cr4g2aZXxBKFivj4N32cNZMBMPt84nHF04C6\nfR2xbrnI/7tyCdYuGGNZOAWVEOW/37qFzSgtTljGjahx/YUYY5sYYy9J/81VVTXFGPsbgM0Afp2f\nW82NqKttiTyzPbOXEWWAmVP1oyof0amcvr6M6MPyWSMxXtrXbx1ZrR8XLHijSBo8d+71LtwgX9sN\ntz3rnt4+fNDVa6ufbIc8iDrpx4uT/l88vBO33LcFh452ee7ne4lmVJTmNkOXJxesuSYv+vSc6y50\nVrfiCFIWAICSIu07uy3QO473GCU/uWdHXuVypjkYA/m5v+XSRZ73GoTDx7qx/3CnpY64zLHOHkc3\nrZPr99+FlCj+jJ59Yis2rmPG616R6F7Gjz8/Qd3giqJg8zkzMKWlFn0ps/tbvOaf9NoFolynuCIV\nAykTMcXbDe7wfXa/e9QiwHLCpHp895rlqKsqgaJYJwKqz/LFHPGRKimKY+8B+yjxqHAdJVRVvRvA\n3Q7vrWGMMQB/AOBZqLe+PjvBE79cdOpk3O9Stmz48IrI74EjumqivKbTuRsa/Bm9+vpKpNJpFBcn\njHMtmN6EnUL5vmF6lHNxaTLQd6mqMq82Esm4r8+XlyZRX1NqHFspFFYQP39IVwcrKrbe1zOvvGv7\nGSc6pQFihEP7bd9jLt8ppnCJ19lwMsN//yVTN7cvlXa9j54UUFVTZrs6LC32breyUm2/vK6qGPfe\nuA6A1fUddj+MF2eMxOLZ7qtjAJg+fjgef+EdjB1Vg/r6ShzT94FLbH4/QLvfY+9k2rtU739JoY3E\nz4mx0C2NlcZ7F6ydgj8/9ZaxIhvZVI1vX7sSn7njEcs5nFg4rdHzuDf3d2DNfPu4Aj6gl5UkbM+T\njnv/xiUlmXZqHpVJX2pqdF9Zl5XZr/66e/pQX1+Jt9/XVqPVVaVZ9ZFKfb+8qjrjwaiuKbMo7XX2\nZZ6BH/3LWnzj/i3YufcwNp05w3i9qCiBru5e1/soFoJRT5rfbBy7+ZvW6l+11SVoHp3x6I2UssEl\nqQAAHiRJREFUsjiqq0sCfeciYZzkhj9WlMRwH56VMAg8pWeMfRHAHlVVfwygA4AvX0JbmzWIJkxW\nzxrpaqw7jh2P/B444qoyX9cU8XvNtraj6O1LIZ1KG59ZNasJzcPLcPtPtBKEx3UVrrfeafc8rzhJ\niUl7T9WlSV/3FVeAru5e49hDesDahxa1mD5/9Ig2yBxq77Scd58esDWuqdLXNdsl95nTZ9xWKeJn\nYlJlktOXtNie87PrZ+HOn7+A3z+xC49t3Yu7rlpmOWbOxHrP79Ddpf1GfX2Z31HOVw27H6ZSaTTU\nlmL+5AZf5z572TiMqCnBQqZ9n3ZduWrbGwctn6+v1363Q0IRiFRvCm1tRw19/1MXNjte99yV403v\nffniebj5Xq0k5cGDx1CRjGFKSy3qqoqN42ZPGI6trx+AAuu+enEi5nitlsZK7H73KLqO9zgew1fU\nY+rL0dZ2FOUlCdOWzzV3PYq7rlqGKkEcSN677evtM85fmYxBAbDqhFGebd/hUD3vHy/vw5TR1bj/\nj9qe+HsHjmXXR/S+vndfZmL13v6j2LHrgOl5mTq62nT+K86YavzbeD2dRndPyvU+OoQJ8hFhTLeL\ny+jo6Dada0JjBTafPR273zuK3/99N44d7cL9f9iGHW8fxqfPmeEZh/F+e+Z6rU1V2Pr6Aby8Y7+j\nVyfsCXI2Ptq7AVzIGHsYwE8AXBLqHeXAKfOdc7TzuWddobt28uF6l/MrZelVN9LptMkNDmhue9Zc\na7jSeUWvx17wLmnIXVRV5UWY3GLe4142s8nuIxZ6+lIm/XI+AaiT9hp529q5mJ9/Tct/9NoP5fj1\nGCcc9sB5wAlH3kd0ymEvEQKbZKEVvrf6sbUMXvBBRnTJJRNxx22SMIjFFNz2ycW+te9rK4tx2uKx\nRspahR4977ZlIxqsF3Zq4iCnLmzGNetn4ezlrU4fszzrdh6Lz18wB5tOyxiMK86chhs/Ph93f2G1\n5Vi3IFGu8+4WGJVxNWv3ZbfPLBfLkcuCils1w2tK8b3PrcBHT57keE2O0/Syu6cP6XQar+neIj/V\nz+zgbSvGz/T2pXDzvc/ia/dtMV7zo+Efj8U83fbiguDp7e6pU3KAaTwWw1zWgGm6EmBfKo2fPfQ6\nnn/tAO76hVWVUOaQoOLIU12jCN50IpsAs/2qqp6qquoqVVWXq6r6ZBQ3lg1iGTuZfO5Zl+odM6iE\nXzacoUdZcm69zP+eXCqt1UOyG4w+f8EcXH3eTMtD6Ho+/UFrbqhASVECm8/O5H76DeLp7OpDOp15\n+PnDK9+jGMErPzBbdBGV0fUumtQCfpWhFs2wn3DI+ZvjmqpMedFO6WNyHAAfiA4c7jTyj50Cp0T4\niqxWKmV68TpvQ99f8H1Ltz4lBg/yMToRj2FG6zDbWJDzVo7HlJZajGsyr2hGDi/HR0+ehK98YoHj\ntYqScUsBC45bsJ5Rk93FyPCVNR+DNq6bjFUnjMKiaRmtAbkLyn3yTEnrPpmI+4zIt7+vomTcMNSA\nc1UwL7gXsUsaJ7LZytUCzIJFg/MgsVV6kZdLhCBZpyqEvN3ElEk+Wfnh71/BJ257CPtsSteKrc1T\nWXe+0245LioKThQlF9w6bz5X1hvXTsbSGY24YM1E74NzRA4McRpYNtjcS/sxbTVnJ0Qwqr4CM8cP\nN5Wse/Jlq1yhiLGPY1THyr7NuQFOSefkiKlkYqpOXyozUPidIPiNnHaa8Nn1Oz/R711S/i433kGr\nQ61dMAbrFjRj89nm/PVJegS/XbR5fxOPaRInbuk2zwVUiDp1UQs+f8Ec28nRmrmj3QuKuOC2subv\nudaL5itrvf/UVhbjY6cwTG2xd58CVuPtd+Ip42Q0e/vSOChoW5+7yjPsyJbMytq5wIZf/ASYyd/n\nt4+/CSCzcm8cVmZMBJ3sAf/NXtxp9Sb8XR/juK6/iOgd4IV72n0WRgqDQWWsC4Vh1SXYdNpUS5BF\nFPiNNpdXf0Amd9Ot0IK4srNLrRLhblhjtZFDpGSXPog7rawT8ZghxdgjrM7+tiXjFveb4+43h9wJ\n+wpU3r+LnM7So6fRBM3/LClKYP3qCRZtgrqqEnz/cytwwUnRTxqDoigKkomYbYT0Uy/vwxa1zSij\nCGSnex8WbpO5mA9jzQd0+buK6l7yo5KtwZNxegTv/cMrpudTLqbjFz4hF1fW3/vtNqfDXdHyrNOu\nEdayB+NhfRtMXCh4Tb6d3hdX0zv3WlfMrDmztcdrKYi1zaNmyBjrupCrDhUKfo21XU4gL+Eo56KK\niLNTr7xEeWVdUZrdAAAA3fpD4LSyBjRFIcDsSv3rs28b//abb16cjOP2KxYD8FeUQSbbGhnyd+Lf\n4813vdPc/OJHwa2/SCZipokWoG0B3HLP0/jur14yvX7Z6VPRX7htkwRxg1eWmZ+HxmEZNa9X3zLL\nVXjVuc6V94+EUzGKu8GfeDkT0/Le+8GqX3EMVUKXtpQXAF09fXh510E8oKeIOcWViDj9nl/6wVPG\nv+1kcMVYA77VmU89gyFjrPNVHjPfiG5wL31mWR6R5z66GXyxM3pV3uJ7jHwlkYsW+y7dYGVW1tZ7\n5FsbogjDnImZgLsg6nH1NaW4YeM83HzpwsD3Glbf4sZ6sPZVmY7jvdjb1oHtQjnQbwnlR8VtjHzG\nnIyoKzNde8ce53xcww3u4r7l8Rfy8yeuZnmxGk5YVdLSLppxYaQI8zHnH9ve8zjSGy7O5GSsO7t6\n8fiL1kDXO3/2gvEtYzHFUEtb5FC4KVsDu3xmRlk74eN3D5shY6wHK6KhPXGWe8S17F7ixtptIAwy\nZnBh+9f35B50cc8ftTQ8Y2VtcyN8D1BcnYmuxaBCD+Oaqkwa334Ja3Jt1PEeGrba4Bs/fd74txh/\nwEunynvxUfO1yxaa0uiOu0T88kmk22qQxybIwYLihCCd1laNv318F/Yd7PAM5vSLm0H2KfDnilxE\nxg5egc8Lvip2MoCPbHUXFwK0ydPqE0bj+59b4XhvQZ5XMejQVBzHx+8eNmSsBzjiABD3sKxcjYzD\n3dpuK2u7Gf5LbxzEq7utKrNth+3l/oLA06B45a2Ma916LL9vcS+Qu/YBf/vGQbFbrTutgm66ZD6+\nusk5Alnmoee0wWiI2WoT1RXWyZJclS1qYoq87+m9Z+2uWa1Nit0CLjuO92L77kP49eO78KUfPGUY\n6wmjq3ObrLjYEu4Jm+XTmGbL4mmN3gch82x1Ocgov7nPh06Dfg637Z8gK2unCX/GCxBuWU43yFgP\ncMSO59UJT1/SYvqbF4lwW4Ha2aFv/vwF/KuwGgKAHW8fxt/0YiG5wKOoeY5wKu3DDS6sQrjBA6LZ\nT7JrD6ciMc0jKrOK4h0qbnA7KoXVplNwYb45dWGz43u+osF91IxuO9Rpkojl3qJVs0dhLss+uE52\ng5cUafn3rSOrjYqEfle+2bLf5ySeb5919zrUPPDRD8LuK/xZrJEmkW7bH0c6ui0xF6HcS+hn7GcK\nObc0aryMkxwtvP+Q9hC5pbXJZ1Tfsq/b4mSoN6yZiEtPn+J6XyJcl/u4R541kHGDOxXziIJ10sB9\nw8Z5Oa3gr/qItQ4x9w5ccea0rM870OBRtWJbGr99BPW4/fDZ9bNw9XmzsHi688qQd0u7ALP3jxzH\nV+55Brv2afEXbt8jlU4bqoFAxgjYlaIMgjxBWDqjCYoCvCHkB+c6qZ3c7O4Kn8fqXd/njNZT65ye\nZz8lSf18lyCu69f2tON7167Av35qiel1RVEQs9EbB4DfPLELWwKmHfph0Bnrijy7zAoJr1llRWkS\n3/rMMiOVhyuFuQWzyF1RHFA4x7t7TYXYxfs4Zf4YLJnuT70MAEq4yIKPaHDuEbCLUv+sIEoSJmvm\nmnOWc53Jz5lYb1LeO/pBN17Q8z8neFVUGkR83ybdh+f195enYXrrMM9Vp+EGtxm0f/34Lux+76hR\nhc1PpDKnJ8s60zKr5oxC84gKfH7DbNx6+SKcv3oCOrvMK9dc23fKWOd8cQBoGuZcw1okE7Rl71rm\nbeJWmtRPBkjQeICiZNzWK+Ik4uJVVztbBp2xHl5jTdG62Idk42CAz+DdqCwrMiY0PLDMTZayqqwI\nwx3S3njd5ivvNIvoy/KbQeAr69f3tuNYZ48hTmBnFHnayFOvWCNRo0qpqCxLYvaEjDJdLrnknIVC\n1Cr3dgD5jYDub7a+rgnBiIMfd4u+m2UqUD7gngC7fiD3WbuVtZOd5AYlVzGnuqoS3HTJAkwZW4fG\nujJbo5Prs/KKXhrTDj8KfJyER+rWX3XvnRJzlqn1M3nOVmDG7lofdFld9lHt2gy60WBsYxWuPX82\nrlk/C3OZViJtpS5FN9iRZ8xO8EGDCy94ibecKUmacm6571lbAYNc3Jb84a6rLPZ01XGFLllXGciU\nYAybmKLgM+fONP4OM/1lXFOVadafTZ3lgUw6ncaeNqvMY9QBULngJooir6Tttks+e/5s2/PyrZAo\nS+xycjXWbuVBZX1uN+IuZW9FYoqCT55pH3Tnx0sgt+kZS8fiLIcxTpQvlTne3Yf33v/AUkf78LFo\nVM0GnbEGtNq2M1qHYfPZM4zCAUOBehuvgh189skDObyMq9PD3Hb4uFGvViSXFaGiKChKxnCg/bip\n09vNmDOeAetvLOtDR0WQlYPreYriSKXSpkE/rFzbQuULHzXXwXZSyMt3NHgQ3ERR5OfKzqXttG33\n4DOauE+ue9Z2yFWicu1lblsFQbaJxKCtIx9048cPqmi3EYeJxxTH7YFsxp6lM5psx8Dls0bixFkj\nbT5h5rW3zamqlRFNsgelsR6q1Pusq8pn+LxYu1cncHvc/ueRnZbXcl3VxhQF7R3dRrQq4PzQj6gt\ntV3VRG3obr50IS5exzByuL/9OC/ierCKXXsOVuQ82O/80rvyUaHBu9kHNrnYw6Tto2ziG3Lds7ZD\nlhfO1Tu0bKazQZsdoJgR32/uTaXwzZ+9gIef24trvvOE5biYopjKiYpkY6yLknHTgmTpjEbUVhZj\no89gZfGznV29eO618IPLgCzqWROFx+VnTMXT29t86yfLrmtZ8lEmqOHLpYAHYD+oOa3uP+jqNSYd\n2/S9s2wLNgRh1PByjArJUAPa9+tLpbCnrXD3Z6NmRwhiOvmGu1Qff3EfNq5jeHX3YdRUFGFUfYXN\nnrX1ufB6tKJwg8uu4sYArmo7nALCNqyegFUn+CtTC2QmJr19aVNFLBklpphUC3MlGVdMWxQXncyM\n2BlfnxcmCD/43Su+tyODQivrQcCiqY24+YolvqM629rlPV73zwU1frm6hu0Mc5HDObmh3rm3HXf8\nbCuA/s/LzYZ4THHVlx6sbDrNf1pfISLuyXZ29eGOn23FDXc/DcBmH9vm5/WaCEdRLXDOpMxqt7go\n7ilT7MXI4eWY3FyDc6Qa46y5NtBKN6NgZjbEstZ4TLGPszlxpv+sE5GS4oRpzEgk/I0fvCynqLOw\n3SG1NQzIWA9B5K7oZeNHDi83SsL5IVc3ODfAIqUeM929BzKBSQPSWMeVvOaLFwpLHWqEc/qz2pYf\nFEXBYr0utSxLKk++jvdYXeVNw8owbWwtFk8bgfNXT7DkLEexsp4lZDOMDsE7lIjHcN2FJ+D0JWNN\nr/cGVPfi+8bf/dXLptc7pdTMmKIYaaic6y+ai0s+lN3EL6Yopj1rv7oJPEtGFEaJcuQhYz0Eycbt\nNaLO3344gEAuJL94DVpiqbp8VsIJi5ii2Ea1DwXkbANxsnXlWfnVBc8G7kn61f+9YXr9l4+a4w/s\ndOcT8Riu3TAHl314GtYuaMa5K81pj1EEmJULQW1RTg/rKoNVOnR6xne9cwRvCW5xRe8fYhUzJWAz\nnS15AbIpyGHssffx0rYpo95CFNCe9RCES3py/Ji2Uxe14JU3/bl4onDdebnTjgt6wgNyZS3d84Y1\nhVeDOipkz04+iyOEATcyT0qVp2RPSfMI7wwFOfgrGZF625xJ9Xh+RzSBUADwuQ2zfSmOiTg9tz9+\ncIfpb7uRIKiwixyx7SQZ7EZC0gd/fseBwOcIAq2shyDZKBZNG1tnqkTkRjLHADM7nB7k0xZreue/\neuwNz2MLGTn4KJ8FAvqbgZ6i9t4hb+3r2/R66V7Ik1K/NdmDwlq09C2uVRA22QR5+tVnsOsvQbuQ\nLGKTjfKlrA8uu+vDhlbWQxC5sys+jZsf9/KGNROzKjPphdOAbqe+1hWR3F+UyBOoiRENooXIAJxb\nmfCaHJ61bBwafKZVyq7gKLxUAHDOqgkoTSihxwTMaB2Gl944mJW+hd+9YruhQAm4W8wzOWa0ajni\nI+q0rcHxI/0H2yWEVLN3DnQYZX2jgow14SnEz/EzqIo612HhVmLPLvI8jHra+UYe8IeSLrg8UVk1\ndzQe3rIHdVXB3Kj9hZ3bfve7mT3WVSf4V1CU96ijir8oLU74EvwIytXnzURfKp1VYJxfj5jdxD1o\nO7HmWtywcZ6hkzByeDlu/Pj8QLE54sr6yz98KtD1s4GM9RBn2YymADNaGwWm0iRmjh+Gv+tFF8Lm\nO1cvR1mJczcNS0Gsv7GTbR0qyP3qvDWTcO7y1gGjPnjO8la89MZB02tiWcjKAJ6mqPao84WiOKuL\neeHX4NrFr2SzkyKnrLU0Zqd6eN8DalafCwrtWQ9RKnUh/CCFKOz2ui88eSLOWDo2rNsy8fXLF7ka\nagB495BVRMSt/nCh4qavPNiRx+hEPDZgDDVgP8i/e9Cqce6HfGiBFyp+V9ZlBdI3HtAlYfPF0O0Z\nQxyey/x8AGk829lrOroShn4GroMWgZdwKmHlm7AqAQ1E5JiJgRggKFemezbLesZRpGoNFPysrOuq\nim2P6w9BoRWzw99GcGPo9gwCgP9KXYD9w1RRloxsX81PVPnaBdZV9ANP53fGGwZDeUUlT/ZyqdrW\nX1x0yiTT3502WuF+ENvCb2GewULcx6T//SPWwh5AONXvgpJNBHkuFIY/geg3lkx3Dt6SEZ+l265Y\njO1vvo/p44bhsE1lnFy49vzZePPdI76iyu3c5BNHD7zgrKFUu1rm/aNm70giHsNA2xSQU6wO2Hh8\ngjIAHUQ5kcukvz9y8ysdampfdc4MoxZ7mGRtrBljkwH8A0CDqqoD7dka8nzz00vx1Pb9OHmef6F9\ncdbfUFOKhtlalGtFaRINtaWhuXKnjauzlPFzoihhDTD7cER76FEylFfWDz+31/T3QHSDJyK458ES\nPOmXEbXWSOymYWXYd9C7uE1/uMEnjanB6PoK7Gk7ZrxWU1GEOZPqI7leViMEY6wKwB0AhqY+4iCg\nuqIYp8wfE0iQwunYRDyGr1++CJvPzr80pN2AFlVuapTIXouhhLwqikoIJErCnGydt3I8AGB6q78J\n62ChuqLYpPddVV7kaYT56rayPP81zxVFweq55rS8qWOj+80Cr6wZYwqA7wP4IoDfhH5HxICkv1So\n7NzgA3GVyu95dH25bwGNwcK5K8fjp399zfh7QK6sHfrc8lnBK0GdNG80GmrLMG1crffBgwwxze2M\npWPxwNNvmd6X68d/ddNC7G07hhG1uZX5zJaElPb6sbX+amBndS23NxljmwBcLb28G8B/q6r6ImMM\niLbQCFFgFBfFMSGAyk+/MAB75MnzRmPbroM4f/XQ0QTnNNWZB9qBuLJ2Coo7aV5wkaBkIo65LBpX\n6kCgOBlHV08fSorikFV35Rzu6vIiVJf3nweiu9ccoBvl1oWrsVZV9W4Ad4uvMcZeA7BJN+SNAB4A\nsNLrQvX12SWcE/7JRxv//GunQVEKU895WuswbHvjIBbOHBXZ/UXVxvX1lfjRjdnV4x3o1Era2vGY\nMuDGi26HGWJDfSXqCzQtr1Db+I5/Wo5d+45gyYwmPPT8O6ZqdBPG1BbUfbeMyqg/rl3UEum9BXaD\nq6pqTP0ZY7sAnOLnc21tR70PIrKmvr5yyLfxtetnAQAOHDjmcWR2UBtHw4iqYsyb3IBnX91vvDbQ\n2vnIYftiHgcPHkNRpIUos6OQ+3JZQsG0MdVoP/wBLjttCm665xl06VX11swZWVD3nRZW1sumN5ru\nLWzDnau/qfB6IUEQA4pEPDYg6la7UVlehNLihGW/failX4XNiLoyXPbhqcbfdtkf/UmdUAY0iowA\nkZzyrFVVbfU+iiAIYnBTnIzjrquWIpUGPnXHo8brVeXhV6AbaogToKgEmLIllsd7I1EUYsBz5rJx\n6Dje09+3QeTI9RfNHdDiMEmbVV9Faf5TigYbYvxJobWneG9RTyPIWBMDnjOXjevvWyBCYMIAVJ4j\noqemIuOdKOjJXMTWuoC/OUEQBDHUKWQDXSKkakWtOEcra4IgiAiopv3qUGisK8OaE0ZjxvjCU3Qr\nLorjMx+Zic6u3kB1y7OBjDVBEEQEXKOnEhK5oSgKPipVNSskZk8cnpfrkLEmCIIIkTs/vRQH24+j\neUThiHcQAx8y1gRBECFSU1GMmopi7wMJIgCFu3NPEARBEAQAMtYEQRAEUfCQsSYIgiCIAoeMNUEQ\nBEEUOGSsCYIgCKLAIWNNEARBEAUOGWuCIAiCKHDIWBMEQRBEgUPGmiAIgiAKHDLWBEEQBFHgkLEm\nCIIgiAKHjDVBEARBFDhkrAmCIAiiwCFjTRAEQRAFDhlrgiAIgihwyFgTBEEQRIFDxpogCIIgChwy\n1gRBEARR4JCxJgiCIIgCJxH0A4wxBcAeADv0l55UVfX6UO+KIAiCIAiDwMYawHgAW1RVPSPsmyEI\ngiAIwko2xnougFGMsYcAdAK4RlXVHR6fIQiCIAgiS1yNNWNsE4CrpZevBHCrqqq/ZIwtBXA/gAUR\n3R9BEARBDHmUdDod6AOMsVIAvaqq9uh/71FVdXQUN0cQBEEQRHbR4P8CfbXNGJsF4K1Q74ggCIIg\nCBPZ7FnfBuB+xtiHAPQC+Hiod0QQBEEQhInAbnCCIAiCIPILiaIQBEEQRIFDxpogCIIgChwy1gRB\nEARR4JCxJgiCIIgCJ5tocF8wxmIA/g3ATABdAC5VVXVnVNcbjDDGkgD+E0ALgGIAtwDYDuBHAFIA\nXgawWVXVNGPsMgCXQ4vQv0VV1T/oOfH3A6gHcBTARlVVD+T9iwwAGGMNALYAWAOtbX8EauNQYYx9\nEcCHASQBfAfAE6B2Dg19zP0hgEnQ2vQyAH2gNg4FxthCALepqrqKMTYBObYrY2wRgLv0Yx9UVfWr\nbtePcmV9FoAiVVWXAPgCgDsivNZg5aMA2lRVXQ5gHYDvQmvH6/XXFABnMsYaAVwFYAmAtQC+zhgr\nAvApAC/ox94H4Mv98B0KHn1S9H0AHdDa9E5QG4cKY2wlgMX6eLASQCuoL4fNKQDKVVVdBuCrAG4F\ntXEoMMauA/ADaIsmIJwx4nsALtB/r4WMsdlu9xClsV4K4M8AoKrqUwDmRXitwcovoInQANpv1QPg\nBFVVH9Nf+xOAkwDMB/CEqqo9qqoeAfA6NI+G8Rvo/z8pXzc+wPgGgH8HsE//m9o4fE4B8BJj7NcA\nfgfgtwDmUjuHSieAar0yYjWAblAbh8XrAM6BZpiBHMcIxlgltMXsLv31B+DR3lEa6yoAR4S/+3Q3\nDeETVVU7VFU9pv+wv4A2IxPb8Ci0h7IKQLvD60ek1wgBxtjHoXkvHtRfUpB5IAFq47Coh1YE6FwA\nVwD4Caidw+YJACUAXoXmKfo2qI1DQVXV/4Xmrubk2q6yffRs7yiN5xEAleK1VFVNRXi9QQljbAyA\nhwDcp6rqT6HtkXCqAByGta0rbV7nrxFmLgFwMmPsYQCzAdwLzbBwqI3D4QC0fblevUrfcZgHJ2rn\n3LkO2sqOQevL90GLD+BQG4dHruOwfCw/hyNRGusnAHwIAPSN9BcjvNaghDE2AsCDAK5TVfVH+svP\nM8ZW6P8+FcBjAJ4GcCJjrJgxVg1gCrSgB+M3EI4lBFRVXaGq6kpVVVcB2ArgYgB/pjYOncehxV2A\nMTYSQBmAv1E7h0o5Mqu1Q9ACiGm8iIac2lVV1aMAuhljrfq2xSnwaO/I5Eb1G+DR4ABwCdW9DgZj\n7FsAzgOgCi//EzT3VhGAVwBcpkchXgotCjEG4Guqqv5Kj0K8F0ATtIj8C1VV3Z/P7zCQ0FfXnwSQ\nhhZMQm0cIoyx2wGsgtZ+XwTwJqidQ4MxVgPgHgDDoa2o74KW4UBtHAKMsbEAfqKq6hLG2ETk2K56\ndPldAOIAHlBV9Qa365M2OEEQBEEUOBTwRRAEQRAFDhlrgiAIgihwyFgTBEEQRIFDxpogCIIgChwy\n1gRBEARR4JCxJgiCIIgCh4w1QRAEQRQ4/x/0gNJ6qYKwKQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x109e27750>"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Make histograms of the samples for each parameter. Should you include all of the samples? "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.hist(chain[1000:, 0], alpha=0.5)\n",
      "plt.hist(chain[1000:, 1], alpha=0.5);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAFVCAYAAADPM8ekAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFudJREFUeJzt3X1w5Hd92PG3tDrpJO1KNzrJ0BYmzDTlE3dyTgNNSHGK\n7akbYo87NoxLJ9AM0AaHh6HQpsM0F5fOdMxDS6FTt44nxUlNcEsndg1tuMFPmGBzEx6TkqN0vmDS\n0JpgczqtHleWdLvbP3bPCPl80q4eft9dvV8zDKuffj/t53uy9N7d32p3oNlsIkmS8jVY9ACSJOnS\njLUkSZkz1pIkZc5YS5KUOWMtSVLmjLUkSZkbutQnI+II8DvAjwEjwG3A/wbuBhrAN4B3pJSaEfEW\n4BbgPHBbSulURIwC9wAzwBLwxpTS7D6tRZKkvrTdPes3AGdTSq8CfhG4A/gwcLK9bQC4MSJeCLwT\neCXwauADETEMvA34envf3wVu3Z9lSJLUv7aL9b3AezftuwG8LKX0WHvbZ4BrgZ8BTqeUNlJKi8AT\nwBXAlcAD7X0faO8rSZI6cMmHwVNKKwARUaEV7luBf7NplyVgEpgAFp5n++KWbZIkqQOXjDVARLwY\nuB+4I6X0iYj415s+PQHM0wpyZdP2ykW2X9h2Sc1mszkwMLCz6SVJ6n3bRm+7J5i9AHgIeHtK6XPt\nzX8cEVellD4PXAd8Fvgy8L6IGAGOApfTevLZaeB64CvtfR9jGwMDA5w9u7Tdblmbman0/BrAdeSk\nH9YA/bGOflgDuI6czMxUtt1nu3vWJ2k9dP3eiLhw7vpdwO3tJ5B9E7iv/Wzw24HHaZ3bPplSWouI\nO4GPRcTjwBrw+u6WIknS4bXdOet30YrzVldfZN+7gLu2bFsFXreL+SRJOvR8URRJkjJnrCVJyty2\nzwaXJKmX1Ot1FhZ++MdHk5PHKJVKBU60e8ZaktRXFhbm+dSZTzM+UWZlcZmbTtzA1NTxosfaFWMt\nSeo74xNlypPb/0lUr/CctSRJmTPWkiRlzlhLkpQ5Yy1JUuZ8gpkkqedt/nOtarVKo9kseKK9Zawl\nST1v859r/eB7T1OZmqD1xpD9wYfBJUl94cKfa41XxoseZc8Za0mSMmesJUnKnLGWJClzxlqSpMwZ\na0mSMuefbmlPbH1Luk70w9vXSdJ+MtbaEwsL89z3yBnGyp39XWNteZGbrz3R829fJ0n7yVhrz4yV\nJyhXjhU9hiT1Hc9ZS5KUOWMtSVLmjLUkSZkz1pIkZc5YS5KUOWMtSVLmjLUkSZkz1pIkZc5YS5KU\nOWMtSVLmjLUkSZkz1pIkZc5YS5KUOWMtSVLmjLUkSZkz1pIkZc5YS5KUOWMtSVLmjLUkSZkz1pIk\nZc5YS5KUOWMtSVLmhooeQJKk/dKoN6hWq89+PDl5jFKpVOBE3THWkqS+tbpS48G5Rzk+c5yVxWVu\nOnEDU1PHix6rY8ZaktTXxspjlCcrRY+xK56zliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTM\nGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJypyxliQp\nc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJytzQTnaKiFcAH0wpXRMRPw38PvDt\n9qd/M6V0b0S8BbgFOA/cllI6FRGjwD3ADLAEvDGlNLvnq5AkqY9tG+uIeA/w94Hl9qaXAx9JKX1k\n0z4vBN7Z/two8IWIeBh4G/D1lNK/jIi/B9wKvHtvlyBJUn/byT3rJ4DXAh9vf/xy4KURcSOte9fv\nBn4WOJ1S2gA2IuIJ4ArgSuBftY97APjnezi7JEmHwrbnrFNK99N6aPuCLwH/NKV0FfCnwL8AKsDC\npn2WgElgAljcsk2SJHVgR+est/hkSulCmD8J/HvgMVrBvqACzNMKdWXLtm3NzFS23ylzvbiGer3O\n/PwPv0Xnzp1jcIdPQRwc3ODo0SOMj490dJ2N+jDT0xWOH9/ff69e/H5s1Q9rgP5YRz+sAfprHYOD\n64yeHWZsfISjo0cYGh56zuX6xtqB/L7ZD93E+oGI+Ecppa8A1wJfBb4MvC8iRoCjwOXAN4DTwPXA\nV4DraEV9W2fPLnUxVj5mZio9uYa5uXPc98gZxsoTAIyNDVOrre/o2NmnnqQ8Oc3QkfGOrrNWW2d2\ndolGY7jjeXeqV78fm/XDGqA/1tEPa4D+W8fc3BKrtXVKR9Z4ZnWDwfMNais/enn1AH7fdGMnN5o6\niXWz/f9vBe6IiA3g+8AtKaXliLgdeJzWQ+snU0prEXEn8LGIeBxYA17fyQJ08MbKE5QrxwAYHx9h\nsLS2o+NWlhe230mS1JUdxTql9GfAK9uXvw78/EX2uQu4a8u2VeB1u55SkqRDzBdFkSQpc8ZakqTM\nGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJypyxliQp\nc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJ\nypyxliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc0NFD6DDrdGoU61Wuz5+cvIYpVJpDyeSpPwYaxVq\ntbbMqdPnmJq+rONja8uL3HztCaamju/DZJKUD2Otwo2NT1CuHCt6DEnKluesJUnKnLGWJClzxlqS\npMwZa0mSMmesJUnKnLGWJClzxlqSpMwZa0mSMmesJUnKnLGWJClzxlqSpMwZa0mSMmesJUnKnLGW\nJClzvkWmJKkn1et1zp07x9zcEtVqlUazWfRI+8ZYS5J60sLCPA9+5yEGh4b5wfeepjI1AUwUPda+\nMNaSpJ41PlGmdGSElcXlokfZV56zliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJn\nrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJytxQ0QNob9XrdRYW\n5rs6tlqt0mw293giSdJuGes+s7Awz32PnGGsPNHxsbNPPUl5cppK54dKkvaRse5DY+UJypVjHR+3\nsrywD9NIknbLc9aSJGVuR/esI+IVwAdTStdExI8DdwMN4BvAO1JKzYh4C3ALcB64LaV0KiJGgXuA\nGWAJeGNKaXYf1iFJUt/a9p51RLwH+Cgw0t70EeBkSulVwABwY0S8EHgn8Erg1cAHImIYeBvw9fa+\nvwvcuvdLkCRpe416g2q1ytzcOebmzlGv14seacd28jD4E8BraYUZ4GUppcfalz8DXAv8DHA6pbSR\nUlpsH3MFcCXwQHvfB9r7SpJ04FZXajz4nUd5+Lt/wKfOfLrrv5wpwrYPg6eU7o+Il2zaNLDp8hIw\nCUwAC8+zfXHLtm3NzFR2slvWilrD4OA6Y2PDjI+PbL/zFqOjw5SGjvzIsTv9Ohc7dj+PA2jUh5me\nrnD8+Pb/1v43lY9+WEc/rAF6fx2Dg+twFsbGRzg6eoSh4aFtL5cnR5n5C9Msje3890cOunk2eGPT\n5QlgnlaQN6+4cpHtF7Zt6+zZpS7GysfMTKWwNczNLVGrrTNYWuv42NXVdUqlAVZWWseOj488e7nT\nY7u9zk7UauvMzi7RaAxfcr8ivx97pR/WAP2xjn5YA/THOubmWvPXVtZ4ZnWDwfONHV9e3eHvj4Ow\nkxtN3Twb/I8j4qr25euAx4AvA38zIkYiYhK4nNaTz04D12/ZV5IkdaCTe9YXXtrq14CPtp9A9k3g\nvvazwW8HHqd1A+BkSmktIu4EPhYRjwNrwOv3cHYdco1GnWq1uu1+g4Prz94C32xy8hilUmk/RpOk\nPbWjWKeU/ozWM71JKX0buPoi+9wF3LVl2yrwut0OKV3Mam2ZU6fPMTV92SX3GxsbplZb/5FtteVF\nbr72BFNTx/dzREnaE76CmXra2Pj2r9Y2Pj7S1Tl8ScqFr2AmSVLmjLUkSZkz1pIkZc5YS5KUOWMt\nSVLmjLUkSZkz1pIkZc5YS5KUOWMtSVLmjLUkSZkz1pIkZc7XBs9QvV5nYWFHb/39HNVqlWazuf2O\nkqSeYawztLAwz32PnGGsPNHxsbNPPUl5cppK54dKkjJlrDM1Vt7+3aQuZmV5YR+mkSQVyXPWkiRl\nzlhLkpQ5Yy1JUuaMtSRJmTPWkiRlzlhLkpQ5Yy1JUuaMtSRJmTPWkiRlzlhLkpQ5Yy1JUuaMtSRJ\nmTPWkiRlzlhLkpQ5Yy1JUuaMtSRJmTPWkiRlzlhLkpQ5Yy1JUuaMtSRJmTPWkiRlzlhLkpQ5Yy1J\nUuaMtSRJmTPWkiRlzlhLkpQ5Yy1JUuaMtSRJmTPWkiRlzlhLkpQ5Yy1JUuaMtSRJmTPWkiRlbqjo\nAaQiNBp1qtVq18dPTh6jVCrt4USS9PyMtQ6l1doyp06fY2r6so6PrS0vcvO1J5iaOr4Pk0nScxlr\nHVpj4xOUK8eKHkOStuU5a0mSMuc9a0lSz6jX6ywszANQrVZpNhoFT3QwjLUkqWcsLMzzqTOfZnyi\nzA++9zQveNE0QyOjRY+173wYXJLUU8YnypQnK4xXxose5cAYa0mSMmesJUnKnOes98nmJ0F0qlqt\n0mw293giSVKvMtb7ZGFhnvseOcNYeaLjY2efepLy5DSVzg+VJPUhY72PxsrdvejGyvLCPkwjSepV\nnrOWJClzxlqSpMwZa0mSMmesJUnKnLGWJClzxlqSpMwZa0mSMmesJUnKXNcvihIRfwRcePWOPwU+\nANwNNIBvAO9IKTUj4i3ALcB54LaU0qldTSxJ0iHTVawj4ihASumaTdv+B3AypfRYRNwJ3BgRXwTe\nCbwcGAW+EBEPp5TWdz+6JEmHQ7f3rH8KGIuIB9tf4zeAl6WUHmt//jPALwB14HRKaQPYiIgngCuA\nr+5ubEmSDo9uz1mvAB9KKb0aeCvwn7d8fgmYBCb44UPlm7dLkqQd6vae9beAJwBSSt+OiHPAT2/6\n/AQwDywClU3bK0B1uy8+M1PZbpfsTU9XGBsbZnx8pONjR0eHKQ0dyeLYnX6dbq/3oNa6dZ/dXG+j\nPsz0dIXjxw/2v9N++LmA/lhHP6wBenMdg4PrjJ4dZmx8hKOjRwCevTw0PLTjy/WNtUJ+jrvVbazf\nTOvh7HdExF+kFeGHIuKqlNLngeuAzwJfBt4XESPAUeByWk8+u6SzZ5e6HCsPMzMVZmeXqNXWGSyt\ndXz86uo6pdIAKyvFHjs+PrLjr9Pt9R7EWi+2jt1cb622zuzsEo3GcMfHdmtmptLzPxfQH+vohzVA\n765jbm6J1do6pSNrPLO6QXl4iNpK6/Lg+caOL68W8HP8fHZyo6nbWP828J8i4sI56jcD54CPRsQw\n8E3gvvazwW8HHqf1kPtJn1wmSVJnuop1Suk88MsX+dTVF9n3LuCubq5HkiT5oiiSJGXPWEuSlDlj\nLUlS5oy1JEmZM9aSJGWu6zfykCSpVzXqDarVH75G1+TkMUqlUoETXZqxliQdOqsrNR6ce5TjM8dZ\nWVzmphM3MDV1vOixnpexliQdSmPlMcqTvfFyo56zliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8Za\nkqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJypyx\nliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJn\nrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjJnrCVJypyxliQpc8ZakqTM\nGWtJkjJnrCVJypyxliQpc8ZakqTMGWtJkjI3VPQAufvDr/5P/nx2taNjyuURnvr+WZ5eavLSyrF9\nmkySdFgY622sbzQpjb+go2MGx0YojUN94Qf7NJUk6TDxYXBJkjJnrCVJypyxliQpc56zliQdao16\ng2q1CsDk5DFKpVLBEz2XsZY61GjUn/3B7lSuvwikw2x1pcaDc49ydGSEm07cwNTU8aJHeg5jLXVo\ntbbMqdPnmJq+rKPjasuL3HztiSx/EUi5+8IffYE1NlhaXGR+YJ7yZGVPv/5YeYzR0dE9/Zp7yVhL\nXRgbn6Ds39BLB2ZuY57SzDBrA+d5ZvaZosc5cD7BTJKkzBlrSZIyZ6wlScqcsZYkKXPGWpKkzBlr\nSZIyZ6wlScqcsZYkKXPGWpKkzO37K5hFxCDwm8AVwBrwKyml7+z39UqS1C8O4p71TcBwSumVwD8D\nPnwA1ylJUkcuvPvW3Nw55ubOUa/Xix7pWQcR6yuBBwBSSl8C/voBXKckSR1ZXanx4Hce5eHv/gGf\nOvNpFhbmix7pWQfxRh4TwOKmj+sRMZhSahzAde9e8zy16p93dsz6MKvzszyzsszyUuff7NXlJQaH\n1lheOlrosY36MLXa+r5e70Gs9WLrKOLfuLa8uP1Oki7q/Mo6tfOrrCwus7paY3lhiZWlFQYGGjSb\nJVaWVigNDz27vZvL9fXzlIbzfH+rgWazua9XEBEfBr6YUrq3/fH/Sym9eF+vVJKkPnIQD4OfBq4H\niIifA/7kAK5TkqS+cRD39z8J/O2ION3++M0HcJ2SJPWNfX8YXJIk7Y4viiJJUuaMtSRJmTPWkiRl\nzlhLkpS57P76OyIGgCeBb7U3/WFK6WSBI3UtIn4C+CJwWUppZ68ukpGIGAf+C3AMWAfemFLq8BVi\nihURk8A9QAUYBv5JSumLxU61OxHxGuDmlNIbip5lJ/rt/QEi4hXAB1NK1xQ9Szci4gjwO8CPASPA\nbSml3y92qs5ERAn4KPBSoAm8NaX0v4qdqnsRcRnwNeBvpZS+dbF9crxn/ZeBr6WUrmn/r1dDPUHr\nddCfKXqWXfgV4CsppatoBe89Bc/TjX8MPJxSuhp4E3BHodPsUkT8O+D9wEDRs3Sgb94fICLeQysS\nI0XPsgtvAM6mlF4F/CLwHwqepxs3AI2U0s8DtwLvK3ierrVvPP0WsHKp/XKM9cuBvxQRj0bEqYh4\nadEDdar96MBvAb8OrBY8TtdSShfCAK1b4dUCx+nWvwX+Y/vyEXr4+9F2GngbvRXrfnp/gCeA19Jb\n//5b3Qu8t315EDhf4CxdSSn9d+BX2x++hN783XTBh4A7ge9faqdCHwaPiH8IvHvL5rcD708p/beI\nuJLWPbqfPfDhduh51vBd4L+mlP4kIqAHfrCfZx1vSil9LSI+C/wk8AsHP9nObbOGFwIfB9518JN1\n7hJr+b2IuLqAkXajt98fYJOU0v0R8ZKi59iNlNIKQERUaIX7N4qdqDsppXpE3A28Bri54HG6EhFv\novUox0MR8etcohXZvShKRIwC51NKG+2Pn0wpvajgsToSEd+mdd4d4OeAL7Ufhu1Z0brVcSql9ONF\nz9KpiDgBfAL4tZTSg0XPs1vtWP9qSumXip5lJ/rt/QHasf5ESulvFD1LtyLixcD9wB0ppbsLHmdX\nIuIFwJeAy1NKPfXIWUR8ntY59ybw14AE3JhSenrrvtk9wYzWwzNzwIci4qeA/1vwPB1LKf2VC5cj\n4v+Q+T3S59O+pfdkSunjtM6n9NzDZRHxV2nde/i7KaUzRc9zSJ0G/g5wr+8PULx23B4C3p5S+lzR\n83QjIn4ZeFFK6QO0Tm012v/rKe3nAwEQEZ+jdSP8OaGGPGP9QeCeiLieVhzeVOw4u5bXQxed+W3g\nYxHxD4ASvfm67u+n9Szw29unJOZTSq8pdqRdu3BLvFf04/sD9NK//1YngUngvRFx4dz1dSmlXnoy\n7H3A3e17pkeAd6WU1gqeaV9l9zC4JEn6UTk+G1ySJG1irCVJypyxliQpc8ZakqTMGWtJkjJnrCVJ\nypyxliQpc/8f2AqYwdPZAicAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10a2c4510>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It's also sometimes useful to view a two-dimensional histogram:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.hist2d(chain[1000:, 0], chain[1000:, 1], bins=20)\n",
      "plt.xlabel('intercept')\n",
      "plt.ylabel('slope')\n",
      "plt.grid(False);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAFgCAYAAABEyiulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF1JJREFUeJzt3X+0XfVZ5/H3vTeXkJCQDtAYAo6AkgemtYbSxRodpNXa\nTqtxTVNwdVoWP1qLoFZpnU47WNRxBkSXjgyoozS2E7GDLEtdg51a7I9UaFW0hQYQ67NEWXUVyo8R\nhJBccnN/+Mc+kcPlJjmkee6935P3ay3WvefsfZ7zZWff8znfs/fZz8js7CySJKkto4s9AEmS9OIZ\n4JIkNcgAlySpQQa4JEkNMsAlSWrQssUewKBmi06XHxkZqSgrSdKhMm9QOQOXJKlBBrgkSQ0ywCVJ\napABLklSgwxwSZIaZIBLktQgA1ySpAYZ4JIkNcgAlySpQQa4JEkNMsAlSWqQAS5JUoMMcEmSGmSA\nS5LUoGbaiVaZnp4uqz02NlZWW5J0eHMGLklSgwxwSZIaZIBLktQgA1ySpAYZ4JIkNcgAlySpQQa4\nJEkNMsAlSWqQAS5JUoMMcEmSGmSAS5LUIANckqQGGeCSJDXIAJckqUEGuCRJDWqmH/js7GxJ3T17\n9pTUrTQ6Wve+a2RkpKx2har9AtrbFpIOL87AJUlqkAEuSVKDDHBJkhpkgEuS1CADXJKkBhngkiQ1\nqOxrZBExBmwBNgCzwGWZef88630Q+MfMvKJqLJIkDZvKGfgmYCYzzwauBK6eu0JEXAq8nC7gJUnS\ngMoCPDNvBS7t3TwJeLJ/eUR8F3AWcAPgFTMkSXoRSo+BZ+Z0RGwFrgdu2nt/RBwP/CzwLgxvSZJe\ntPKT2DLzYrrj4FsiYkXv7vOA44A/At4PvC0iLqweiyRJw6LyJLYLgBMz8xpgApihd6w7M38N+LXe\nehcBp2XmjVVjkSRp2FTOwG8BNkbE7cBtwOXA5oi4ZJ51PYlNkqQXYaSym9OhNDMzUzLQycnJirIA\njI+Pl9S1G9lz7EYm6TAw74tRMwE+NTVVMtDKF+mJiYmSupVjXrlyZUndFsOw6m+jxW0haVHN+6Lh\nldgkSWqQAS5JUoMMcEmSGmSAS5LUIANckqQGGeCSJDXIAJckqUEGuCRJDTLAJUlqkAEuSVKDDHBJ\nkhpkgEuS1CADXJKkBhngkiQ1yACXJKlBzfQDn56eLhno448/XlEWgB07dpTUPfHEE0vqAkxNTZXU\nPfLII0vqLlu2rKQu2Ldb0pJhP3BJkoaFAS5JUoMMcEmSGmSAS5LUIANckqQGGeCSJDXIAJckqUEG\nuCRJDTLAJUlqkAEuSVKDDHBJkhpkgEuS1CADXJKkBhngkiQ1qJl2ojMzMyUDnZmZqSgLwKOPPlpS\n95FHHimpC7B+/fqSuitWrCipu2rVqpK6ULdvVLZAHR31Pbk0hGwnKknSsDDAJUlqkAEuSVKDDHBJ\nkhpkgEuS1CADXJKkBhngkiQ1yACXJKlBBrgkSQ0ywCVJapABLklSgwxwSZIaZIBLktSgsrZIETEG\nbAE2ALPAZZl5f9/yc4H395b9n8y8vmoskiQNm7q+hrAJmMnMsyPi1cDVwJvgX8L9GuBMYCfw1xHx\nkcx8Yl/Fqlo7VrYTXb58eUndxx57rKRupbGxsZK6lf9+p556akndypafVe2Bq/79JB28sleSzLwV\nuLR38yTgyb5l08BpmbkDeCkwBkxWjUWSpGFTegw8M6cjYitwPXDTnGUzEfFm4MvA54BdlWORJGmY\nlJ/ElpkX0x0H3xIRK+Ys+wPgBGA5cGH1WCRJGhaVJ7FdAJyYmdcAE8AM3QlrRMTRwMeB12XmZETs\nBKarxiJJ0rCpnIHfAmyMiNuB24DLgc0RcUlmPg18BLgjIj5PF+4fKRyLJElDZaTqrNVDbWpqqmSg\nlWcxP/300yV1v/jFL5bUBVi7dm1JXc9Cf86KFSsOvNJBGhkZKanrWejSopr3D9sLuUiS1CADXJKk\nBhngkiQ1yACXJKlBBrgkSQ0ywCVJapABLklSgwxwSZIaZIBLktSgyn7gh9SyZTVDrbySV9UVt844\n44ySugDPPvtsSd277767pO7OnTtL6gKsXLmypG7V1dIA1q9fX1J3fHy8pO4RRxxRUhdq+65LS4F7\nuCRJDTLAJUlqkAEuSVKDDHBJkhpkgEuS1CADXJKkBhngkiQ1yACXJKlBBrgkSQ0ywCVJapABLklS\ngwxwSZIaZIBLktQgA1ySpAaNzM7OLvYYBtXMQPeampoqqfvQQw+V1AXYvn17Sd0HH3ywpO7JJ59c\nUhdg9erVJXUzs6QuwFlnnVVSt6qd7ymnnFJSF+Coo44qqWubUi2CeXsQuydKktQgA1ySpAYZ4JIk\nNcgAlySpQQa4JEkNMsAlSWqQAS5JUoMMcEmSGmSAS5LUIANckqQGGeCSJDXIAJckqUEGuCRJDTLA\nJUlqkAEuSVKD7AdeaGZmpqRuVZ9xgEcffbSk7sTEREnde+65p6QuwM6dO0vqTk5OltStVNUPfNOm\nTSV1AV7ykpeU1K18zVy+fHlZbTXNfuCSJA0LA1ySpAYZ4JIkNeiAB7YiYjnwXiCAnwQuB67JzP0e\nyIuIMWALsIHu+PVlmXl/3/K39mpNAfcBP5aZzR3nliRpMQwyA/8NYBVwJl3YfhvwoQEetwmYycyz\ngSuBq/cuiIgVwH8HXtNbvqa3viRJGsAgAX5mZl4BTGbmM8CFwCsP9KDMvBW4tHfzJODJvsXPAt+Z\nmc/2bi8Dak5TliRpCA3y3ZCZiDii7/ZxwEDfj8rM6YjYCmwGzuu7fxZ4HCAifgI4KjM/M+igJUk6\n3A0yA78O+AywLiKuA+4C/uegT5CZF9MdB9/S++gcgIgYjYhfAV4LnPtiBi1J0uHugDPwzLwxIu4C\nXgOMAZsy894DPS4iLgBOzMxr6D4en+H5F2O5ge6j9M2evCZJ0oszyFnoRwCvB76X7iS2iYi4b4DQ\nvQXYGhG3A+N0Z5xvjohVwJeAdwB3ANsiAuC6zPy/B/1/IknSYWSQY+C/DRwJfJBuBn4h8HK6QN6n\nzJwA3rKfVcYGHKMkSZpjkAA/Czh974w7Iv4QuH//D5EkSZUGOYnta8ApfbfXAg/XDEeSJA1i0BZD\n90TEZ+iOgX8P8FBEfBKYzczvLxudJEma1yABflXv596T1n69b5lnj+/H6GjNpeZHRubtLHdIHH/8\n8SV1d+3aVVJ3zZo1JXUBHnjggZK6N910U0ldgGOOOaakblVr1ZUrV5bUBdi4cWNJ3RNOOKGkLsDY\nWM2pQVWvRVV1NZhBvkb2JxHx/XTf114GbOtdZU2SJC2SA759ioj3AT8HfBV4EPhARHygemCSJGnf\nBvkI/QLgrN7XwoiIDwJ309ecRJIkLaxBDmCM0F0xba9ngT01w5EkSYMYZAa+Dbil15RkBLiod58k\nSVokgwT4u4HL6K7ANkoX3jdUDkqSJO3fPgM8Iv51381P9P7baz3wD1WDkiRJ+7e/GfgdvPB73v23\nT0GSJC2KfZ7ElpknZebJwFvpLt5yOvB3wNHA+xdmeJIkaT6DnIV+HV37zzcDu4BXYoBLkrSoBgnw\n0cy8HfgB4GOZ+Q/YClSSpEU1SIDvioj30l1K9f9FxOXAjtphSZKk/RkkwM8HVgJvzswngHXA20pH\nJUmS9muQZiZfA/5b3+0rSkckSZIOaNB+4FpCxsfHy2rPzMyU1K1qG1m5LVavXl1S97zzziupC7Bt\nW81FEq+99tqSum984xtL6gJEREndhx9+uKQuwLp160rqrlq1qqTu9PR0SV2oa606TGzmKklSgwxw\nSZIaZIBLktQgA1ySpAYZ4JIkNcgAlySpQQa4JEkNMsAlSWqQAS5JUoMMcEmSGmSAS5LUIANckqQG\nGeCSJDXIAJckqUEGuCRJDRqZnZ1d7DEMqpmBauHs3r27rHZVP+KvfOUrJXUBtm/fXlL3U5/6VEnd\nkZGRkroAb3jDG0rqPvbYYyV1Ac4999ySuitXriype/TRR5fUBRgdrZlfNtpnfN4/FGfgkiQ1yACX\nJKlBBrgkSQ0ywCVJapABLklSgwxwSZIaZIBLktQgA1ySpAYZ4JIkNcgAlySpQcuqnyAixoAtwAa6\ny6Felpn3z1lnJfBp4B2ZmdVjkiSpdQsxA98EzGTm2cCVwNX9CyPiVcAdwMl4vXNJkgZSHuCZeStw\nae/mScCTc1Y5AngT4MxbkqQBlX+EDpCZ0xGxFdgMnDdn2Z8BRMRCDEWSpKGwIAEOkJkXR8T7gb+I\niNMzc2KhnluLr6pt7fLly0vqQt2YN2zYUFIX6tpGVtV98MEHS+oC3HvvvSV1n3rqqZK6AHfeeWdJ\n3XXr1pXUPeOMM0rqAixbVhNPjbYTnVf5R+gRcUFEXNG7OQHM4LFuSZK+IQtxEtstwMaIuB24Dbgc\n2BwRlyzAc0uSNJTKP0LvfVT+lgHW+57qsUiSNCy8kIskSQ0ywCVJapABLklSgwxwSZIaZIBLktQg\nA1ySpAYZ4JIkNcgAlySpQQa4JEkNMsAlSWqQAS5JUoMMcEmSGjRS1fO4QDMDlfan8m9u165dJXWf\neOKJkro7duwoqQuwffv2krpV2xjghBNOKKm7fv36krqV+/Jpp51WUreyH/j4+HhV6ZH57nQGLklS\ngwxwSZIaZIBLktQgA1ySpAYZ4JIkNcgAlySpQQa4JEkNMsAlSWqQAS5JUoMMcEmSGmSAS5LUIANc\nkqQGGeCSJDXIAJckqUG2E5WGSEN/zwBMTEyU1Z6cnCypW9VaFeChhx4qqbt79+6SupX72yte8YqS\nuscdd1xJXShtVWo7UUmShoUBLklSgwxwSZIaZIBLktQgA1ySpAYZ4JIkNcgAlySpQQa4JEkNMsAl\nSWqQAS5JUoMMcEmSGmSAS5LUIANckqQGGeCSJDVo2WIPQNKhMzIyb9fBJWv58uVltataO65du7ak\nLsD09HRJ3UceeaSkbqVnnnmmpG7lPrdmzZqSuvv6u3YGLklSg8pm4BExBmwBNgCzwGWZeX/f8h8E\nfgaYAj6cmb9dNRZJkoZN5Qx8EzCTmWcDVwJX710QEePArwKvA14N/EhE1H0uJUnSkCkL8My8Fbi0\nd/Mk4Mm+xacDD2TmU5m5B/gCcE7VWCRJGjalJ7Fl5nREbAU2A+f1LToaeKrv9g6g5ui/JElDqPwk\ntsy8mO44+JaIWNG7+ylgdd9qq3n+DF2SJO1H5UlsFwAnZuY1wAQwQ3cyG8DfAKdGxL8CdtJ9fP7L\nVWORJGnYVM7AbwE2RsTtwG3A5cDmiLikd9z7p4A/Bv4M+FBmfr1wLJIkDZWR2dnZA6+1NDQzUEmD\nqbpwCcDU1FRJ3T179pTUBfj612vmMS1eyGX9+vUldY899tiSulB6IZd5r+TihVwkSWqQAS5JUoMM\ncEmSGmSAS5LUIANckqQGGeCSJDXIr5FJGkozMzMldStfMycnJ0vq7ty5s6Turl27SupC3VffNm7c\nWFIXYNmymmujjY6O+jUySZKGhQEuSVKDDHBJkhpkgEuS1CADXJKkBhngkiQ1yACXJKlBBrgkSQ0y\nwCVJapABLklSgwxwSZIaZIBLktQgA1ySpAYZ4JIkNch2opKGUtVr28jIvJ0dD4mqMe/evbuk7tjY\nWEldgKmpqZK6lf9+Rx55ZFVp24lKkjQsDHBJkhpkgEuS1CADXJKkBhngkiQ1yACXJKlBBrgkSQ0y\nwCVJapABLklSgwxwSZIaZIBLktQgA1ySpAYZ4JIkNcgAlySpQQa4JEkNsh+4JA25qtf56enpkrpQ\n17e7sh/46GjZnNh+4JIkDQsDXJKkBhngkiQ1yACXJKlBBrgkSQ0ywCVJatCyqsIRMQ58GPgWYDlw\nVWZ+vG/5W4H/DDwLfDQzr60aiyRJw6ZyBn4+8HhmngO8Afj1vQsi4ljgF4DvBf4d8B8i4ozCsUiS\nNFQqA/yjwM/2Pc9U37JvBe7JzH/KzFngTuCcwrFIkjRUyj5Cz8ydABGxmi7MP9C3+G+Bl0XEWuAZ\n4LXAH1SNRZKkYVN6EltEfDOwDbgxM2/ee39mPgm8B/gYcBNwN/D/K8ciSdIwKQvwiPgm4FPA+zJz\n65xly4BXZeZ3A28BvgP4bNVYJEkaNmXNTCLiOuCHgOy7ewtwVGZuiYifAd4ETAO/lZkfPkBJm5lI\n0kGwmUl9XVj4ZiZ2I5OkIWeA19cFA1ySJA3AK7FJktQgA1ySpAYZ4JIkNcgAlySpQQa4JEkNMsAl\nSWpQ2bXQDwcRcRpdI5a1mTk5Z9klwI/QNXG5KjM/sQhDXHQRcRTd5XJfAkwCF2Xmw3PWuY6uK90O\nuu/7vykzn17osS6mAbeT+xQQEWuAjwCrgSOAn8rMO+esc9jvUzDwtnK/6hMRm4HzMvP8eZYtqf3K\nAD9IEXE08D/o+pnPXbYO+AngTGAF8IWI+PTckD9MvBP4YmZeFREXAe8D3j1nnVcCr8/MJxZ8dEvH\nfreT+9TzvAf4dGZeHxEbgN+j2y793Kc6+91W7lfP1wvo1wNf3scqS2q/8iP0gxARI8ANwBXAxDyr\nnAX8aWbu6b07ewB4xQIOccnIzOvoer8DfAvwZP/yiBgFTgW2RMQXIuLtCzzEJeFA2wn3qX7XAh/s\n/T7OnL9B96nn2e+2wv1qrj8FfpR5rny2FPcrZ+AHEBE/zAtnjF8Fbs7MeyMCXviPvRp4qu/2DmBN\n2SCXiH1sq4sz866I+Czwcrp3t/1WAtcDv0q3P34uIr6UmfeVD3iRHOR2cp96zt5ttQ74XeDyOcsP\nu30KDnpbuV895+LM/P2IeM0+Hrbk9isD/AAy80PAh/rvi4i/BX64txOsA/4YeE3fKk/T/WHstZoX\nzqiGznzbqm/Za6N7t/MJ4Nv6Fu0Crs/MZwEiYhtdd7qhfbE9yO3kPtUnIr6d7uPg/5SZn5+z+LDb\np+Cgt5X71eCW3H5lgB+EzDx17+8R8SAvnC39JXB1RCwHjgROB/5q4Ua4dETEFcDXMvN3gZ10J8o8\nbxXg9yLilcAYcDawdUEHuQQMsJ3cp3oi4t8AHwV+aB+zH/epngG2lfvV4JbcfmWAf+P+pRtMRLwH\neCAzPx4R1wOfpzvP4KcP15NC6N7l/k5EvINup387vGBb3Qj8ObAH2JqZX1m00S6eQbaT+1TnF+jO\nqL6+dwjrnzJzs/vUvAbZVu5XzzfLvl/Xl9R+ZTcySZIa5FnokiQ1yACXJKlBBrgkSQ0ywCVJapAB\nLklSgwxwSZIaZIBLQyYizoyILftZ/oO977YuuIj4+Yg4ezGeWxo2XshFGjKZeRdwyX5WOZO+C1Us\nsHOAbYv03NJQ8UIu0pDpNWP4r3Qh/ZfAdwMvpWsb+VXgc71l/wX4GPAbwMvorgD3S5l5c0RcDFwE\nHAv8IbAF+N+9OruAd2bmfRFxIV2DjFHgLuDHM3N3RHwN+Cywka5BxvnAq3vP9XVgc2beX7kdpGHn\nR+jS8Om/FOR4Zn4XXV/oq3qXfvxN4Dcz83eAK4EvZear6AL2AxFxcu+xJwAbM/NK4H8BH83Mb6d7\nc3Bl7zrb7wS+MzPPAB4H3tt77Hrgk5n5HcDNdE0gbgS+RBf+hrf0DfIjdGm43db7eT9wTO/3/va3\n3wes6F2DHbqWiS+jewNwd2bO9O4/B3gLQGZ+EvhkRLyLrj/yX/Sus30E3Swc4OnMvLn3+43ANX3P\n+YJey5JePANcGm67ez9neX5w7p2hjwLnZ+Z2gF7f6H8E3gZM9K2/p//xvdn3KPD7mXl5775VPPea\n0t9NbXTObY/bSYeAH6FLw+dAM9w9wHjv923AjwFExPHAl4FvnqfGHcB/7K33OuAG4E+AzRHx0ogY\nofto/id76x8TEf++9/vbgT/q/T7V99ySvgEGuDR8Zvt+zs5z/x3A+RHx48DP032Efh/dSWfvy8y/\nn+ex7wLOjYgvAz8HXJKZ9/Yev43nekj/Yu/nHuCCiLgHeB3w7t79twG/FRH/9pD8n0qHMc9Cl3TI\nRcREZq5Y7HFIw8wZuKQKzgykYs7AJUlqkDNwSZIaZIBLktQgA1ySpAYZ4JIkNcgAlySpQQa4JEkN\n+megWzD7SPzrcQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10a25a650>"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Report to us your constraints on the model parameters -- you have some freedom in interpreting what this means..."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "theta_best = np.mean(chain[1000:], 0)\n",
      "print(\"slope, intercept =\", theta_best)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "slope, intercept = [-2.74153022  3.1508756 ]\n"
       ]
      }
     ],
     "prompt_number": 13
    }
   ],
   "metadata": {}
  }
 ]
}