{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Bayesian model comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will do an example of model comparison in the Bayesian framework. To do so, we'll simulate some data with a quadratic relation and fit both linear and quadratic functions to it to see if we can prefer one model to another. As you might expect, the ability to distinguish may depend on the uncertainties in the data points.\n",
    "<br>\n",
    "In the Bayesian framework, we compare two models with the posterior odds ratio:\n",
    "$${P(M_1 | D) \\over P(M_2|D)} = {P(D|M_1) P(M_1)\\over P(D|M_2) P(M_2)}$$\n",
    "In the case where there is no prior preference for one model over the other, $P(M_1) = P(M_2)$, this reduces to the Bayes factor:\n",
    "$$BF = {P(M_1 | D) \\over P(M_2|D)} = {P(D|M_1)\\over P(D|M_2)}$$\n",
    "For a model with parameters:\n",
    "$$P(\\theta|M,D) = {P(D|\\theta,M) P(\\theta|M)\\over P(D|M)}$$\n",
    "To compute the Bayes factor, we need to compute the evidence, $P(D|M)$ for each of the two models. This is given by:\n",
    "$$P(D|M) = \\int P(D|\\theta,M) P(\\theta|M) d\\theta$$\n",
    "Equivalently, this is the integral over the posterior:\n",
    "$$P(D|M) = \\int P(\\theta|M,D) d\\theta$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's simulate a data set with a second order polynomial, but with significant uncertainties, and see if we can distinguish between a linear and a quadratic fit. Use numpy.polyval() to generate the data, and add noise assuming a normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'y')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAVR0lEQVR4nO3df4zc913n8ecbx75sS2HjxgnJhp4TLrKAi5oNe2mCAZUE6rRXtSZqSRA6DERKEdeD6u6c2oAqQEBcIn6eECS0cLlTCQ7GcaIKSKO4BV0vdbqpk7glNU7T1MR27W0Tt0ANTtZv/pjvlvF6Zr3emc/Md+f7fEirmfnOdzxvfz1+7Xc+8/m8JzITSVJzfMOwC5AkDZbBL0kNY/BLUsMY/JLUMAa/JDXMecMuYDEuvPDCXLt27bDLkKRl5YknnvhSZq6Zv31ZBP/atWuZnp4edhmStKxExBc6bXeoR5IaxuCXpIYpFvwRsS4inmz7+WpEvCciVkfEIxFxoLq8oFQNkqQzFQv+zNyfmVdn5tXAdwFfAx4AtgCPZuaVwKPVbUnSgAxqqOdG4HOZ+QXg7cC91fZ7gY0DqkGSxOCC/1bgvur6xZl5BKC6vGhANUiSGEDwR8Qq4G3An53j426PiOmImJ6ZmSlTnCQ10CDO+N8MfCozj1a3j0bEJQDV5bFOD8rMezJzKjOn1qw5Y/2BJGmJBhH8P8K/DfMAPARsqq5vAh4cQA2StKzccvdj3HL3Y0X+7KLBHxGvAn4Q2Nm2eRvwgxFxoLpvW8kaJEmnK9qyITO/Brx23rYv05rlI0kaAlfuSlLDGPyS1DAGvyQ1jMEvSQ1j8EtSwxj8ktQwBr8kNYzBL0k1s2vvIfYePM6ez7/I+m272bX3UF//fINfkmpk195DbN25j5OzpwA4dPwEW3fu62v4G/ySVMBSe+3c9fB+Trw8e9q2Ey/PctfD+/tVmsEvSZ2UbJK2kMPHT5zT9qUw+CWpRi4dHzun7Uth8EtSjWzesI6xlStO2za2cgWbN6zr23MU7c4pSTo3GycnALhjx9OcnD3FxPgYmzes+/r2fjD4JalmNk5OcN/jBwHY/q7r+/7nO9QjSQ1j8EtSwxj8ktQwBr8k9Vnplgu9MvglaZ5egnsQLRd6ZfBLUpteg3sQLRd6ZfBLUpteg3sQLRd65Tx+SWrTa3BfOj7GoQ77nmvLhRLz9+d4xi9JbXrtlTOIlgu9MvglqU2vwb1xcoI7b76KVSta8ToxPsadN1/V15YLvXKoR5La9KNXTumWC70y+CVpnroHd68c6pGkhjH4JalhDH5JahiDX5IaxuCXpIYx+CWpYYoGf0SMR8SOiPhsRDwTEddHxOqIeCQiDlSXF5SsQZJ0utLz+H8H+KvMfEdErAJeBfwc8GhmbouILcAW4L2F65Ckgarz/P/IzDJ/cMQ3AU8BV2Tbk0TEfuCNmXkkIi4BPpaZC66Fnpqayunp6SJ1StKoiognMnNq/vaSQz1XADPAH0fE3oj4QES8Grg4M48AVJcXdSn49oiYjojpmZmZgmVKUrOUDP7zgGuA38/MSeCfaA3rLEpm3pOZU5k5tWbNmlI1SlLjlAz+F4AXMnNPdXsHrV8ER6shHqrLYwVrkCTNUyz4M/OLwN9HxNz4/Y3A3wIPAZuqbZuAB0vVIEk6U+lZPf8N+FA1o+c54Cdo/bK5PyJuAw4C7yxcgySpTdHgz8wngTM+UaZ19i9JGgJX7kpSwxj8kkbSLXc/xi13PzbsMmrJ4JekhjH4JalhDH5JahiDX5IaxuCXpIYx+CWpYQx+SWoYg1+SGsbgl1RLvSzA2rX3EHsPHmfP519k/bbd7Np7qM/VLW8Gv6SRsmvvIbbu3MfJ2VMAHDp+gq079xn+bQx+SSPlrof3c+Ll2dO2nXh5lrse3j+kiurH4Jc0Ug4fP3FO25vI4Jc0Ui4dHzun7U1k8EsaKZs3rGNs5YrTto2tXMHmDeu6PKJ5Sn8DlyQN1MbJCQDu2PE0J2dPMTE+xuYN676+XQa/pBG0cXKC+x4/CMD2d10/5Grqx6EeSWoYg19S7bgAqyyDX1KtuACrPINfUq24AKs8g19SrbgAqzyDX1KtuACrPINfUq24AKs85/FLqhUXYJVn8EuqnX4swHLhVncO9UhSwxj8ktQwBr8kNYzBL0kNY/BLUsMY/JLUMEWnc0bE88A/ALPAK5k5FRGrge3AWuB54Icz86WSdUiS/s0gzvi/PzOvzsyp6vYW4NHMvBJ4tLotSRqQYSzgejvwxur6vcDHgPcOoQ5JNeYCrHJKn/En8JGIeCIibq+2XZyZRwCqy4s6PTAibo+I6YiYnpmZKVymJDVH6TP+9Zl5OCIuAh6JiM8u9oGZeQ9wD8DU1FSWKlCSmqboGX9mHq4ujwEPANcCRyPiEoDq8ljJGiRJpysW/BHx6oh4zdx14E3Ap4GHgE3VbpuAB0vVIGl4brn7MW65+7Fhl6EOSg71XAw8EBFzz/MnmflXEfFJ4P6IuA04CLyzYA1SY82Frh+Sar5iwZ+ZzwGv77D9y8CNpZ5XkrQwV+5KUsMY/JLUMAa/JDWMwS9JDWPwS1LDGPyS1DAGv6S+27X3EHsPHmfP519k/bbd7Np7aNglqY3BL6mvdu09xNad+zg5ewqAQ8dPsHXnPsO/Rgx+aQQN84z7rof3c+Ll2dO2nXh5lrse3j+wGrQwg18aMcM+4z58/MQ5bdfgGfzSiBn2Gfel42PntF2Dd9bgj4h3R8QFgyhGUu+Gfca9ecM6xlauOG3b2MoVbN6wbiDPr7NbzBn/twCfjIj7I+KmqNptSqqnYZ9xb5yc4M6br2LVila8TIyPcefNV7FxcmIgz6+zO2vwZ+YvAFcCHwR+HDgQEb8WEd9WuDZJS9CvM+5e+ulvnJxg8nXjvOHy1Xx8yw2Gfs0saow/MxP4YvXzCnABsCMifr1gbZKWwDNunc1Z+/FHxM/Q+qasLwEfADZn5ssR8Q3AAeCOsiVKOlcbJye47/GDgF/EojMt5otYLgRuzswvtG/MzFMR8dYyZUmSSjlr8Gfm+xa475n+liNJKs15/JLUMAa/VFO9zKqRFmLwS1LDGPyS1DCLmdUjqWHmunuenD3F+m272bxh3TmvA3AaaX0Z/NKIWmrwduvuCbgIbEQ41CPpNMPu7qnyDH5Jpxl2d0+VZ/BLOs2wu3uqPINf6qKp8+jtpz/6DH6phob5nbl29xx9zuqRaqYOs2rs7jnaPOOXClnqUJGzalSawS/VjLNqVJrBL9WMs2pUWvHgj4gVEbE3Ij5c3b48IvZExIGI2B4Rq0rXIC0nzqpRaYM44/9ZoP0LW94P/FZmXgm8BNw2gBqkZcNZNSqtaPBHxGXAf6b1Xb1ERAA3ADuqXe4FNpasQVqONk5OMPm6cd5w+Wo+vuUGQ199VfqM/7dpfRn7qer2a4HjmflKdfsFoOMrOiJuj4jpiJiemZkpXKYkNUex4K++iP1YZj7RvrnDrtnp8Zl5T2ZOZebUmjVritQodTPMBVRSaSUXcK0H3hYRbwHOB76J1juA8Yg4rzrrvww4XLAG6ZzVYQFVHbhwa3QVO+PPzK2ZeVlmrgVuBXZn5o8CHwXeUe22CXiwVA0aruXa66YfC6h8x6A6G8Y8/vcC/z0inqU15v/BIdQgddXrAqpu7xgMf9XFQII/Mz+WmW+trj+Xmddm5n/IzHdm5r8MogZpsXpdQGXLBdWdK3dVW8MaKup1AZUtF1R3Br80T68LqGy5oLqzLbPUQS9tiTdvWMfWnftOG+5ZSssFZ9WoFINf6rO5dwZ37Hiak7OnmBgfY/OGdY2aCqp6M/ilAvwiE9WZY/yS1DAGvyQ1jMGvIuqwcnW5rhyWSjP41XeuXJXqzeBX39nrRqo3Z/Wo70r1uoHBdsd0No5GlWf86jt73Uj1ZvCr7+x1I9WbQz3qu15Xrl46PsahDiG/3HrdOFSkuvKMX0X08mXhvb5jkLQwz/hrbG4OetPOHO11I5Vl8KuWeu11Mzcd9OTsKdZv2+0vDqmNQz0aOS4gkxZm8GvkOB1UWpjBr5HjdFBpYQa/Ro5ffSgtzODXyHE6qLQwZ/WomGFNQ3U6qLQwg18jya8+lLpzqEeSGsYz/ppyAZJn6lIpnvHXkAuQJJVk8NdQXRYg+Z210mgy+GvIBUiSSjL4a8gFSJJKMvhryAVIkkpyVk9BS+2n7wIkSSUZ/DXVjwVITf0ilzlN/XtLZ1NsqCcizo+IxyPiqYj4TET8UrX98ojYExEHImJ7RKwqVYMk6Uwlx/j/BbghM18PXA3cFBHXAe8HfiszrwReAm4rWIMkaZ5iwZ8t/1jdXFn9JHADsKPafi+wsVQNkqQzFZ3VExErIuJJ4BjwCPA54HhmvlLt8gLQ8RPLiLg9IqYjYnpmZqZkmepgrmXEns+/yPptu101LI2QosGfmbOZeTVwGXAt8O2dduvy2Hsycyozp9asWVOyTM1jywhptA1kHn9mHgc+BlwHjEfE3Gyiy4DDpZ7XlgNLU5eWEZLKKDmrZ01EjFfXx4AfAJ4BPgq8o9ptE/BgqRq0NLaMkEZbyXn8lwD3RsQKWr9g7s/MD0fE3wJ/GhG/AuwFPliwhqHpR1vlYc1Dv3R8jEMdQt6WEdJoKBb8mfk0MNlh+3O0xvtHVrcxcmBZrL7dvGEdW3fuO224x5YR0uiwV08BdRgj72VWzsbJCe68+SpWrWi9PCbGx7jz5quWxS8tSWdny4YChj1G3o93HH5nrTS6PONfwFJnBQ27rXId3nFIqi+Dv4Bht1Ue9jsOSfVm8Bcw7DHyYb/jkFRvBn8hGycnmHzdOG+4fDUf33LDQD8YHfY7Dkn1NrLB3+ReM8N+xyGp3kZyVs9yn0ffD87KkdTNSJ7xO6tFkrobyeB3VoskdTeSQz32mukPh4ik0TSSZ/z9mNXS5A+HJY22kTzjn/sA944dT3Ny9hQT42Pn1B3TD4cljbKRDH7obVbLQh8OG/ySlruRDf5e9OvDYcfIJdXRSI7x98qWB5JGmcHfwai0PNj+rut91yHpDA71dNDrh8OSVGcGfxe2PJA0qhzqkaSGMfglqWFGeqjHIRpJOpNn/JLUMAa/JDWMwS9JDWPwS1LDGPyS1DAjPaunV84KkjSKPOOXpIYx+CWpYQx+SWoYg1+SGsbgl6SGKRb8EfGtEfHRiHgmIj4TET9bbV8dEY9ExIHq8oJSNUiSzlTyjP8V4H9k5rcD1wH/NSK+A9gCPJqZVwKPVrclSQNSLPgz80hmfqq6/g/AM8AE8Hbg3mq3e4GNpWqQJJ1pIGP8EbEWmAT2ABdn5hFo/XIALurymNsjYjoipmdmZgZRpiQ1QmRm2SeI+Ebgr4FfzcydEXE8M8fb7n8pMxcc54+IGeALSyzhQuBLS3zsIFhfb6yvN9bXm7rX9+8zc838jUVbNkTESuDPgQ9l5s5q89GIuCQzj0TEJcCxs/05nQo/hxqmM3NqqY8vzfp6Y329sb7e1L2+bkrO6gngg8AzmfmbbXc9BGyqrm8CHixVgyTpTCXP+NcD/wXYFxFPVtt+DtgG3B8RtwEHgXcWrEGSNE+x4M/M/wdEl7tvLPW8HdwzwOdaCuvrjfX1xvp6U/f6Oir+4a4kqV5s2SBJDWPwS1LDjEzwR8RNEbE/Ip6NiDPaQETEv4uI7dX9e6pFZYOqrWPfonn7vDEivhIRT1Y/7xtUfdXzPx8R+6rnnu5wf0TE71bH7+mIuGaAta1rOy5PRsRXI+I98/YZ6PGLiD+KiGMR8em2bYvqQxURm6p9DkTEpk77FKrvroj4bPXv90BEjHd57IKvhYL1/WJEHGr7N3xLl8cu+H+9YH3b22p7vm3SyvzHFj9+PcvMZf8DrAA+B1wBrAKeAr5j3j4/DfxBdf1WYPsA67sEuKa6/hrg7zrU90bgw0M8hs8DFy5w/1uAv6T1gf11wJ4h/lt/kdbClKEdP+D7gGuAT7dt+3VgS3V9C/D+Do9bDTxXXV5QXb9gQPW9CTivuv7+TvUt5rVQsL5fBP7nIv79F/y/Xqq+eff/BvC+YR2/Xn9G5Yz/WuDZzHwuM08Cf0qrJ1C79h5BO4Abq7UGxWX3vkXLyduB/5MtnwDGqwV4g3Yj8LnMXOpK7r7IzL8BXpy3eTF9qDYAj2Tmi5n5EvAIcNMg6svMj2TmK9XNTwCX9ft5F6vL8VuMxfxf79lC9VW58cPAff1+3kEZleCfAP6+7fYLnBmsX9+nevF/BXjtQKprM69v0XzXR8RTEfGXEfGdAy0MEvhIRDwREbd3uH8xx3gQbqX7f7hhHj9YXB+quhzHn6T1Dq6Ts70WSnp3NRT1R12Gyupw/L4XOJqZB7rcP8zjtyijEvydztznz1NdzD5FRatv0Z8D78nMr867+1O0hi9eD/wvYNcgawPWZ+Y1wJtptdD+vnn31+H4rQLeBvxZh7uHffwWqw7H8edptU3/UJddzvZaKOX3gW8DrgaO0BpOmW/oxw/4ERY+2x/W8Vu0UQn+F4Bvbbt9GXC42z4RcR7wzSztreaSROe+RV+XmV/NzH+srv8FsDIiLhxUfZl5uLo8BjxA6y11u8Uc49LeDHwqM4/Ov2PYx69ydG74K7r3oRrqcaw+TH4r8KNZDUjPt4jXQhGZeTQzZzPzFPCHXZ532MfvPOBmYHu3fYZ1/M7FqAT/J4ErI+Ly6qzwVlo9gdq19wh6B7C72wu/36oxwU59i9r3+Za5zxwi4lpa/zZfHlB9r46I18xdp/Uh4Kfn7fYQ8GPV7J7rgK/MDWsMUNczrWEevzaL6UP1MPCmiLigGsp4U7WtuIi4CXgv8LbM/FqXfRbzWihVX/tnRj/U5XkX83+9pB8APpuZL3S6c5jH75wM+9Plfv3QmnXyd7Q+8f/5atsv03qRA5xPa4jgWeBx4IoB1vY9tN6OPg08Wf28Bfgp4Keqfd4NfIbWLIVPAN89wPquqJ73qaqGuePXXl8Av1cd333A1ID/fV9FK8i/uW3b0I4frV9AR4CXaZ2F3kbrM6NHgQPV5epq3yngA22P/cnqdfgs8BMDrO9ZWuPjc6/BuVlulwJ/sdBrYUD1/d/qtfU0rTC/ZH591e0z/q8Por5q+/+ee8217Tvw49frjy0bJKlhRmWoR5K0SAa/JDWMwS9JDWPwS1LDGPyS1DAGvyQ1jMEvSQ1j8EtLEBH/qWomdn61WvMzEfEfh12XtBgu4JKWKCJ+hdaK8DHghcy8c8glSYti8EtLVPWK+STwz7RaRMwOuSRpURzqkZZuNfCNtL5V7fwh1yItmmf80hJFxEO0vgHqcloNxd495JKkRTlv2AVIy1FE/BjwSmb+SUSsAP5/RNyQmbuHXZt0Np7xS1LDOMYvSQ1j8EtSwxj8ktQwBr8kNYzBL0kNY/BLUsMY/JLUMP8KnP3iiq2o2t0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(42)\n",
    "pars=[0.1,0.5,25]  # true model parameters\n",
    "sig=3    # uncertainty on the data points\n",
    "xdata=np.arange(0,20)\n",
    "ydata = np.polyval(pars,xdata) + np.random.normal(0,sig,len(xdata))  # calculate the dependent variable with noise\n",
    "\n",
    "plt.errorbar(xdata, ydata, yerr=sig, marker='o',ls='none')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you tell by eye whether a linear or quadratic fit is preferred?\n",
    "<br>\n",
    "<B> ANSWER HERE: </B>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Linear model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with a linear model:\n",
    "$$y = \\alpha x + \\beta$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the normal Bayesian framework, let's supply routines for the log(prior), log(likelihood), and log(posterior). To keep things simple, let's just use a uniform prior in all parameters; however, we need to be careful here, because the prior must be normlized, so we need to set it to a value such that it integrates to unity over our trial grid. So we will first set up the trial grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "dalpha=0.05   # set sampling for alpha (slope parameter)\n",
    "alphas=np.arange(1,4,dalpha)  # set range for trial slopes\n",
    "dbeta=0.5     # set sampling for beta (intercept)\n",
    "betas=np.arange(15,35,dbeta)  # set range for trial intercepts\n",
    "\n",
    "def log_prior(theta):\n",
    "    \"\"\" Return log(prior)\n",
    "        For a constant uniform prior, need this to integrate to unity over the trial grid\n",
    "    \"\"\"\n",
    "    return np.log(1/len(alphas)/len(betas))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to supply routines for log(L) and log(posterior)\n",
    "As we've done many times, we will assume Gaussian uncertainties on the points, so we get:\n",
    "\n",
    "$$L(x|\\mu,\\sigma) = \\prod {1\\over \\sigma\\sqrt{2\\pi}} \\exp {-0.5(x_i-\\mu)^2\\over \\sigma^2}$$\n",
    "\n",
    "$$log(L) = -N \\log \\sigma -N \\log(2\\pi) + \\sum{-0.5(x_i-\\mu)^2\\over\\sigma^2} $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_likelihood(theta, x, y, sigma):\n",
    "    \"\"\" calculate log(likelihood) given polynomial parameters theta, and data points x,y,sigma\n",
    "    \"\"\"\n",
    "    y_model = np.polyval(theta, x)\n",
    "    n=len(y)\n",
    "    return np.sum(-0.5*(y-y_model)**2/sigma**2) - n*np.log(sigma) - n*np.log(2*np.pi) # expression for log(L)\n",
    "\n",
    "def log_posterior(theta, x, y, sigma):\n",
    "    \"\"\" calculate log_posterior from log_prior and log_likelihood\n",
    "    \"\"\"\n",
    "    return log_prior(theta) + log_likelihood(theta, x, y, sigma) # expression for log(posterior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the posterior, by calculating it over a grid in parameter space. You'll have to choose an appropriate range of parameters to sample and a spacing in each: you can probably get these from inspection of the data. Try some initial guesses and if your posterior is not roughly centered on a peak (in plots below), you can/should adjust them. Also, if the peak is not well-sampled, make your grid more finely spaced; since we will be integrating with the rectangle rule, we want the peak to be well-sampled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate log(posterior PDF)\n",
    "log_Ppdf=np.zeros([len(alphas),len(betas)])\n",
    "for i,alpha in enumerate(alphas) :\n",
    "    for j,beta in enumerate(betas) :\n",
    "        log_Ppdf[i,j] = log_posterior(np.array([alpha,beta]),xdata,ydata,sig)  #use your function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's display it using the plot2() routine from last time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine plots...\n",
    "def plot2(grid,xgrid,ygrid,vr=None,xt=None,yt=None) :\n",
    "    \"\"\"function to plot likelihood and marginal distributions\n",
    "        grid : input likelihood to display\n",
    "        xgrid : range of x values\n",
    "        ygrid : range of y values\n",
    "        vr : range for scaling log(L) image\n",
    "        xt, yt : titles for axes\n",
    "    \"\"\"\n",
    "    \n",
    "    # make a 2x2 plot grid\n",
    "    fig,ax=plt.subplots(2,2)\n",
    "    if vr is None : \n",
    "        # if no range is given, use min and max of data array\n",
    "        vr = [grid.min(),grid.max()]\n",
    "    # image of surface\n",
    "    ax[0,0].imshow(grid,vmin=vr[0],vmax=vr[1],interpolation='nearest',origin='lower',\n",
    "                   aspect='auto',extent=(xgrid.min(),xgrid.max(),ygrid.min(),ygrid.max()))\n",
    "    \n",
    "    # marginal distribution summed over y axis\n",
    "    ax[1,0].plot(xgrid,grid.sum(axis=0))\n",
    "    ax[1,0].set_xlabel(xt)\n",
    "    \n",
    "    # marginal distribution summed over x axis, flipped\n",
    "    ax[0,1].plot(grid.sum(axis=1),ygrid)\n",
    "    ax[0,1].set_ylabel(yt)\n",
    "    \n",
    "    # get rid of unused plot\n",
    "    ax[1,1].remove()\n",
    "    fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the posterior: if you don't have a well-centered, well-sampled, peak, adjust the grids above and rerun. We will need to integrate this, so you don't want significant probability to fall off of your grid!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3hcd33n8fdXo9FtZEmWJd9kO3JiNzfHcYJzI0mXpmEbCN3Q0gDZFkLL85juLhQe2N7Yp1tgt/eW0n26CzVNILQsISVZYFMuTSFZ4pA4sR07seNcnIvju+SLrPtt5rt/zBlbc5E0sjVzzkif1xM9M+fMTzNfn2jmM+d3fud3zN0RERGJmqqwCxARESlEASUiIpGkgBIRkUhSQImISCQpoEREJJIUUCIiEknVpXjSGqv1OhKleGqRshpmgFEfsXK8lpnFgG3AIXd/11Rt29ravLOzsxxliZTc9u3bj7t7e+76kgRUHQmus58vxVOLlNVW/1E5X+7jwF6gabqGnZ2dbNu2rfQViZSBme0vtF5dfCIRYGYrgNuBvw+7FpGoUECJRMMXgN8BUmEXIhIVCiiRkJnZu4Aud98+TbtNZrbNzLZ1d3dP2u6Jfcf51AO76Bsem+1SRcpKASUSvhuBf2dmbwD3A7eY2T/mNnL3ze6+0d03trfnHU8+49Xufh7ccZDRce2MSWVTQImEzN1/391XuHsn8H7gx+7+ayGXJRI6BZSIiERSSYaZi8i5cffHgMdCLkMkErQHJSIikaSAEhGRSFJAiYhIJCmgREQkkhRQIiISSQooERGJJAWUiIhEkgJKREQiSQElIiKRpIASEZFIUkCJiEgkKaBERCSSFFAiIhJJCigREYkkBZSIiESSAkpERCJJASUiIpGkgBIRkUhSQImISCQpoEREJJIUUCIiEknTBpSZ1ZnZ02a2y8z2mNlny1GYiIjMb9VFtBkBbnH3fjOLA1vM7Pvu/lSJaxMRkXls2j0oT+sPFuPBj5e0KpF5Rj0VIvmKOgZlZjEz2wl0AY+4+9bSliUy72R6Kq4ENgC3mdn1IdckEqqiAsrdk+6+AVgBXGtm63LbmNkmM9tmZtvGGJntOkXmNPVUiOSb0Sg+d+8BHgNuK/DYZnff6O4b49TOUnki88d0PRUTvwR2d3eHU6RIGRUziq/dzFqC+/XArcCLpS5MZL6Zrqdi4pfA9vb2cIoUKaNi9qCWAY+a2XPAM6S/2T1c2rJE5q+peipE5pNph5m7+3PAVWWoRWTeMrN2YMzdeyb0VPxZyGWJhKqY86BEpPSWAfeZWYx0z8YD6qmQ+U4BJRIB6qkQyae5+EREJJIUUCIiEkkKKBERiSQdgxIpATO7HbgcqMusc/fPhVeRSOXRHpTILDOzLwHvAz4GGHAncEGoRYlUIAWUyOx7q7t/EDjl7p8FbgBWlrsITeQnlU4BJTL7hoLbQTNbDowBq8v14mYGgCuhpMLpGJTI7Hs4mL/yL4AdpHdm/r5cL14TSwfUyHiyXC8pUhIKKJHZ9+fuPgI8aGYPkx4oMVyuF2+sjQPQOzQOC8v1qiKzT118IrPvycwddx9x99MT15XaxUsbAdh96HS5XlKkJLQHJTJLzGwp0AHUm9lVpEfwATQBDeWq48K2Rjpa6vmbH73C2y5pZ/GCuul/SSSCKiugzAqsK2In0FM5yzp6LCXxC8CHSF/P6fMT1vcCny5XEVVVxt994C3c+aUnee+XnuTut3Zyx4YOWhM15SpBZFaYl+DDusla/Tr7+Vl/XgWUlNtW/xG9frLAH97kzOw97v5gqWoC2Lhxo2/btm3KNk/sO86ffH8vuw/1Eo8Zt166hDs2dHDZsiY6FtYTq5rRP0ukZMxsu7tvzF1fWXtQBVgslr1cE89vVJUdYj46mtckb51CTM7dE2Z2D7Dc3d9hZpcBN7j7PeUs4sY1bTz8sZvZe6SXb20/yLefPcT3dx8FoCZWxQWLGriwPcGF7Y2sbktwUXuC1W2N2tOSyKisgCoiNKyuNn9d04Lpn7qvP2s5dbovv81YfrCJFPCV4Oe/BMsvA98EyhpQGZcua+IP3nUZv/eOS9h5oIdXu/p5/fgAr3YPsK+rnx+/2MVY8ux7q6Uhzuq2BMub62lpiNOaqKGloYaFDXEWNtTQEtwubKhhQV01VdoTkxKprIASqQxt7v6Amf0+gLuPm1noJyXFY1Vc09nKNZ2tWevHkykOnhriteP9vNY9wGvHB3itu5+9R3s5NTDK6aExUpN8N4xVGc318QnhFQRZYmKQxYP16XWJ2mrq4zF1Mcq0Kj6gPJnzvs9dBlKJ+qzloZX5e1Tmi7OW698sMET34NHsl+rL38tS16AAA2a2iGC2ITO7HojsmO/qWBWdbQk62xLcckn+46mU0zs8xqnBMU4NjtIzOMqpgcz9s7cnB0Y5eGqQ3YfS60bGU/lPNkFdvIpETTUNtbH0bU2MRG1wm7W+mkRt7MxtfXxCu5z2NbGqMzNpSOWr+IASiaBPAt8FLjSzJ4B24FfCLencVVUZLcHe0WoSRf/e0GiSU4OjWUF2emiMwZEkA6PjDI4mGRg5ezs0lr7t7htJPx60Gx6bOugmqq6y7OCaEGD1OcFXXxOjLh6jLl5FfTwdfHU1MeqqY9TXxCasq6IuuB+P6dTRcqr8gEpl7zGlBobymsT6BrKWk/XNeW1OXJ4z2OKqtrw27Tuzfy+x80Bem2T38axlHx/PayNz3gvA/wEGgT7g26SPQ80r9TUx6mvqWd5SP33jKSRTzmCBQBscTWYF2cT1g6PjDIwmGRxJ3x7rGz7bLridrNtyKtVVRn08Rm08Rn3NhGALfurjsYLBV18T/M6Z5ars38n8XnU6ELUnmFb5ASUSPV8jfe7THwfLdwH/QPqyGzJDsSpjQV2cBXUFRuieI3dnNJlieCzF8FiSodEkQ2PJ9P3M7WgqZznJ8Hh6fWZd5vGh0SQ9g6McHUs/NjSWZDhoP3EASrGqjDOhVT9hD7ChJt3VOfE2UROj/sy6YH1tjIagKzT39yppL3DOBVShkXap7hNZy4n9TXltTlyWvXfU/nOH89qcvC77m2DvD/InqF76r9ltUvsPFVWjzCkXu/uVE5YfNbNdoVUjecyM2uoYtdUxmutnL/gKGUueDbqRTIBNCMSzIZcfiIOZ29HMnmGSkwNDZ5dHxhkcS87o0HdNrGpCaE0MvBgNtdU0xHPu1559vL2xluUt9SxtrqMuHpv+xc7TnAuoQlJD2d1+sVcP5rVZ+lT2dDBvrG3Na/MPN2dPSP3TNWvz2nyl/Reyljsfyj+nJLXvjaxlHxnJayMV7Vkzu97dnwIws+uAJ0KuSUISj1URj1XN6h7gRO7OyHhqQvdmdqBl3Q8CbbBA2+P9owyeHMxaP9XeX1tjLctb6ljenO7GXd5SR0dLPcuC+22J2vM+BWFeBJRIOZjZ86RH7sWBD5rZm8HyBaSPS4nMOjM7czxr0Sw/9+h4KtiTG2dgJElX3zCHe4Y53DPE4Z4hDvUM8UpXH4++1JU3arMmVsWyljouXrKAL7x/Aw01M4+b+RFQOfu/ydO9eU3qnn09a3lFYk1em48tvCtrecfGb+a1ufiDX81a/p343XltVj+Q863ildfz2mivqiK961x/0cxWkj52tRRIAZvd/W9mqzCRmcjslWW6GfuGxzh6ephjvcMcPT3C0d5huvpGOBmcJ1fwlAJLf/SmzuPUm2kDSm8ckeK4+/7z+PVx4FPuvsPMFgDbzewRd9eel0wqs4eTGcV45njVWJLBkfT9obFkXhdf7nGt9O+efY7pRjkuStSwpKmOpc11rF/RwtKmOpY2155Zt7Spjub6+HmPRCxmD2ruvXEKJHry5Kms5cYt+/LaWDJ7r+rioQ/mtbnvmq9kLd9++9a8No+cvj5recVIgYEdb2QfJ9PAirnN3Y8AR4L7fWa2l/SlOyr3fTbPJVOeNeghMzJwKG9dcsK69ECJwaywyA6PzPD6odEk4zMcK58+yTk9MrAhfnaE38KGmjMDIfJHC8ZorI2zpCkdQIubaqmtLv0ACSgioPTGESkvM+sErgK25qzfBGwCWLVqVdnrmivGkymGgz2P3KHiZ4NjBqPtMkPQJ4TM0FiS0Wlm0phM5vypiQFRXxNjaVNdwfDIG2I+MWTiwZDz4ByrSps3cUbHoCZ748wJucepcvaoABI/eTFreVXfRXltPvDOj2a3eUv+MPP+9dlX/+45uDivzcLB7DbJY115bXQS8NxjZo3Ag8An3D3rYKm7bwY2Q/pyGyGUV1aZc5UmnojbPzLxxNz0gfvM7UBwUm7W+szJusFjQ6NJRpMzDw7LnJeUObm25uyJuI211bQ31uadqFufOVk3azn7RN3sdTFqq3WC7kRFB9RUb5zg8TPf7urKd/HQ0inUDdib/c+ufualvDZrj6/MWj6+tyOvTe0F2X+AQ/mTVrBgZfbK2MBAXpu8wR6aB7CimVmc9Hvs6+7+UNj1nCt3D87XGQ2mOhrj1MAofSPjWWGRGzCDo+PpAJowI8RMurASwbk7iZqz8/a11MfpaKkjUVN95qTVidManQmRvHXZUx1pZodwFBVQxbxxJn67a7JWfVKKzIClP/3uAfa6++ena18uk4XNyYH0pLEng4ljzz6eXp5uLyURzJM3ca68hYkaVixsODOH3pkJYs8ET3pdoQll6+OV130l0ytmFF8k3zhRkBoczFtnL2cPGW8/0ZPXZmHnkqzlkdb8a1ilqrOnI6lubMx/raHsbkANTa9oNwIfAJ43s53Buk+7+/dK/cKnB8fY193Pq939vNqVvj14amjasDHjzCU0WhtqWNnawJUrWmhJpJcXJtKX2GhNpC+30VQXJ1FbmcdCJBzF7EGF9sYRmS/cfQtQ0k/tnsFRdh08zb6u7DA63n92hGhNrIrVbQlWtjawfkUzCxM16bAJAqc1cfZihU31cV3TSUqqmFF8JX/jzCW5w8HHu47ntanKuXpvotAVf+uzp14qdHzJqrP/9/lYgUETqdCvkychcne2vn6Sbzz9Jt9//uiZvaHm+jhrFjdyyyWLWbO4kYva0z8rFtZTXUGTicrcNj9mkhCZh36w+wh//sOXeK17gAV11dx17UpuW7eMtUsaWZSo0UF/iTwFVKkV2INJ5YzIK3TsyGpyJpmNFTgxrsDVg0UAnn79JB/938+yZnEjf3nnldx+xTLqa8pzcqXIbFFARUCh85ny1lUV8eHi53ZioMwtyZTzsW/sYGVrA9/8yA0lv5yESKmos1lkjnmlq49jvSN87JY1CiepaNqDqhQa7CBFeuVYehDOZcvzL8wpUkm0ByUyxwyOpruHS3WBPJFyUUCJzDGjwVVQa6v19pbKpr9gkTnGg3PmNIhcKp0CSkREIkkBJSIikaSAEhGRSFJAiYhIJCmgREQkkhRQIiISSQooERGJJAWUiIhEkgJKREQiSQElIiKRpIASEZFIUkCJiEgkKaBERCSSFFAiIhJJCiiRCDCze82sy8x2h12LSFQooESi4avAbWEXIRIlCiiRCHD3nwAnw65DJEoUUCIVwsw2mdk2M9vW3d09abv2xlrecsFCqqv09pbKNu1fsPrGRaLB3Te7+0Z339je3j5pu3dcsYwH/8NbaW6Il7E6kdlXzFesr6K+cRERKbNpA0p94yIiEgZ1UotEgJl9A3gSuNjMDprZh8OuSSRs5u7TNzLrBB5293VTtNkEbAoW1wGVeMyqDTgedhHnQHWXzgXuPvkBn5CYWTewf4omlbBtp1LJ9Vdy7RBO/QXfZ7MWUDntt7n7xplWGDbVXV6VWnclqPRtW8n1V3LtEK361cUnIiKRVMwwc/WNi4hI2VVP18Dd7zqH5918Dr8TBaq7vCq17kpQ6du2kuuv5NohQvUXdQxKRESk3HQMSkREIum8AqrQNEhm9hkzO2RmO4Ofd55/mbPLzFaa2aNmttfM9pjZx4P1rWb2iJm9EtwuDLvWiaaoO9Lb3MzqzOxpM9sV1P3ZYP1qM9sabO9vmllN2LXOBWZ2m5m9ZGb7zOz3wq5nJip5arXJ3p+VYLL3aNjOq4vPzH4W6Ae+lhmCbmafAfrd/S9npcISMLNlwDJ332FmC4DtwLuBDwEn3f1Pgzf2Qnf/3RBLzTJF3e8lwtvczAxIuHu/mcWBLcDHgU8CD7n7/Wb2JWCXu38xzFornZnFgJeBtwMHgWeAu9z9hVALK1Khz5RKMdn7sxK2/WTvUXd/Ksy6zmsPqlKnQXL3I+6+I7jfB+wFOoA7gPuCZveR/vCPjCnqjjRP6w8W48GPA7cA3wrWR257V6hrgX3u/pq7jwL3k/67rgiV+pkClfv+hCnfo6Eq1TGoj5rZc8HueqS6yXIFJyFfBWwFlrj7EUj/sQGLw6tsajl1Q8S3uZnFzGwn0AU8ArwK9Lj7eNDkIBXyZo64DuDAhGVt1xAUeH9GXu571N1Dr70UAfVF4CJgA3AE+KsSvMasMLNG4EHgE+7eG3Y9xSpQd+S3ubsn3X0DsIL0t/xLCzUrb1VzkhVYp+1aRpX6uZL7HjWz0LtYZz2g3P1Y8A9NAV8m/WEUOUE/64PA1939oWD1saAfOdOf3BVWfZMpVHelbHMAd+8BHgOuB1rMLHMu3grgcFh1zSEHgZUTlrVdy2iSz5WKMuE9GvpllmY9oDIf8IFfIoKTxgYHBO8B9rr75yc89F3g7uD+3cB3yl3bVCarO+rb3MzazawluF8P3Eq6f/5R4FeCZpHb3hXqGWBtMEKyBng/6b9rKbEpPlcib5L36IvhVnX+o/i+AbyN9Oy3x4A/DJY3kO5WeAP4SOa4TlSY2U3A48DzQCpY/WnS/cUPAKuAN4E73T0yB2ynqPsuIrzNzWw96UEQMdJfih5w98+Z2YWkD+K3As8Cv+buI+FVOjcEpxl8gfT2vtfd/yjkkopW6DPF3e8JtagiTfb+dPfvhVdVcSZ7j4ZblWaSEBGRiNJMEiIiEkkKKBERiSQFlIiIRJICSkREIkkBJSIikaSAEhGRSFJAiYhIJCmgREQkkhRQIiISSQooERGJJAWUiIhEkgJKREQiqXr6JjPX1tbmnZ2dpXhqkbLavn37cXdvD7sOkfmoJAHV2dnJtm3bSvHUImVlZvvDrkFkvlIXn4iIRNK0AWVmdWb2tJntMrM9ZvbZchRWaXYfOs3Dzx1G19cSEZkdxXTxjQC3uHu/mcWBLWb2fXd/qsS1VQx35xPf3Mm+rn62XHOcz92xjppq7ZyKiJyPaT9FPa0/WIwHP9pNmODJV0+wr6ufGy5cxP3PHOAD92zl5MBo2GWJiFS0or7mm1nMzHYCXcAj7r61tGVVlvuefIPWRA1f+fVr+ML7NvDsgR7e/T+f4OVjfWGXJiJSsYoKKHdPuvsGYAVwrZmty21jZpvMbJuZbevu7p7tOiPrUM8Qj7xwjPdds5K6eIx3X9XB/ZuuZ3A0yXv/7kmGRpNhlygiUpFmdKDE3XuAx4DbCjy22d03uvvG9vb5c9rI159Kj0L+1etWnVl39aqF/MWd6+kZHOPpN06GVZqISEUrZhRfu5m1BPfrgVuBF0tdWCUYHkty/zMHuPXSJaxY2JD12PWrF1ETq+Lxl+fP3qSIyGwqZhTfMuA+M4uRDrQH3P3h0pZVGb73/BFODoxy91s78x6rr4mxsXMhW/YdL39hIiJzwLQB5e7PAVeVoZaKc9+T+7moPcFbL1pU8PGb17bzZz94ka7eYRY31ZW5OhGRyqaTdc7RzgM97DrQwwdv6MTMCra5eW0bgPaiRETOgQLqHH3tyTdI1MT45as7Jm1z2bImWhM1bHlFASUiMlMKqHPg7vz4xS7eccUyFtTFJ21XVWXcuKaNx/cd1xRIIiIzpIA6B2+cGKRncIyNFyyctu3Na9ro7hvhJZ20KyIyIwqoc7DzwCkANqxqmbbtTZnjUOrmExGZEQXUOXj2zR4SNTHWLl4wbdvlLfVc1J7gJwooEZEZUUCdg50Heli/ooVYVeHRe7luXtvO06+fYHhM0x6JiBRLATVDw2NJXjjcy1VFdO9l3LSmjeGxFDv2nyphZSIic4sCaob2HD7NeMrZsLL4gLr+okVUV5m6+UREZkABNUPPvtkDFDdAIqOxtpqrVy1kyz7NyyciUiwF1Aw9e6CHjpZ6Fi+Y2dRFN61tY8/hXk70j5SoMhGRuUUBNUM73+yZ0d5Txg0XLcL97B6YiIhMTQE1A119wxzqGeKqGRx/yrhsWRNmsPvw6RJUJiIy9yigZmBnsPczkxF8GYnaai5sS7D7UO9slyUiMicpoGbg2QM9xGPG5cubz+n313U0s0d7UCIiRVFAzcDON3u4dFkTdfHYOf3+uuXNHDk9rIESIiJFUEAVKZlynjvYM6Pzn3Jd3tEEwJ7D6uYTEZmOAqpIr3T1MTCaPKfjTxmZrkENlBARmZ4CqkhnTtBdOf0lNibTXB9nVWsDezRQQkRkWgqoIu18s4eWhjidixrO63nWdTRpD0pEpAgKqCLtPJA+/mRW3Azmk7l8eTP7TwxyemhslioTEZmbFFBFGBgZ5+WuvvMaIJGxriN9HOoFDZQQEZnStAFlZivN7FEz22tme8zs4+UoLEr2HunFPT1M/Hxdvjwzkk/dfCIiU6kuos048Cl332FmC4DtZvaIu79Q4toiIzMsPDNM/Hy0NdayrLmO3YcUUCIiU5l2D8rdj7j7juB+H7AX6Ch1YVGy5/BpFiVqWNo0sxnMJ3P58iZ2q4tPRGRKMzoGZWadwFXA1gKPbTKzbWa2rbt7bl33aM/hXi5b3nTeAyQyLl/ezGvd/QyOjs/K84mIzEVFB5SZNQIPAp9w97yv/+6+2d03uvvG9vb22awxVKPjKV4+1nfO8+8Vsq6jmZTD3iN9s/acIiJzTVEBZWZx0uH0dXd/qLQlRcvLx/oYS/qZwQ2zYV2HBkqIiEynmFF8BtwD7HX3z5e+pGjJDAefzYBa2lTHokSNBkqIiEyhmD2oG4EPALeY2c7g550lrisy9hw+TaImRueixKw9p5lxeUezrg0lIjKFaYeZu/sWYHZGB1SgPYd7uXRZE1VVs7sJ1i1vYvNPXmNkPElt9bldvkNEZC7TTBJTSKWcvUd6z8z+MJvWdTQznnJePto/688tIjIXKKCm8MaJAQZGk1w2i8efMjKzUjyv41AiIgUpoKawuwQDJDJWttbTXB9XQImITEIBNYU9h08TjxlrFy+Y9ec2M9avaOb5Qz2z/twiInOBAmoKLxzu5WeWLKCmujSbaV1HMy8d7WN4LFmS5xcRqWQKqEm4O3sO95akey9jfUczY0nnpaOaUUJEJJcCahJHe4c5OTA6q1Mc5bpiRfq5n9NxKBGRPAqoSewJTqJdNwuX2JhMR0s9rYkanj+o41AiIrkUUJPYffg0ZnDJ0tIFlJlxRUczzx3UHpSISC4F1CT2HO5ldVuCRG0x13Q8d+tXNPNKV78GSoiI5FBATeKFw70lPf6UcUVHM8mU88IRzcsnIjKRAqqAUwOjHOoZKukIvoz1K1oAeF7dfCIiWRRQBWRG1V1Rgjn4ci1pqqV9Qa2OQ4mI5FBAFbB9/ymqDK5c2VLy1zIz1ndoRgkRkVwKqAJ27D/FJUubaCzxAImMK1Y0s6+rn4GR8bK8nohIJVBA5UimnGffPMVbLlhYttdcv6KZlKOBEiIiEyigcrx4tJeB0WRZAypzvSkdhxIROUsBlWPH/lMAZQ2oxQvqWNZcpxklREQmUEDl2L7/FIsX1LJiYX1ZX3ddR7Pm5BMRmUABlWN7cPzJzMr6uus7mnmte4C+4bGyvq6ISFQpoCbo6h3mwMmhsnbvZWRmNt99SAMlRESgiIAys3vNrMvMdpejoDBtD44/XR1GQAUDJXQ+lIhIWjF7UF8FbitxHZGwff8paqqrWFeGOfhyLWqsZVVrA0+/fqrsry0iEkXTBpS7/wQ4WYZaQrf9zVNcuaK5ZJd4n86Na9p46rUTjCVToby+iEiU6BhUYHgsye5Dp0Pp3su4eW0b/SPj7Dqgbj4RkVkLKDPbZGbbzGxbd3f3bD1t2Tx/6DRjSectq8ILqLdetAgzePyV46HVICISFbMWUO6+2d03uvvG9vb22XrasglzgERGS0MN6zua2bJPASUioi6+wPb9p1jdlqCtsTbUOm5e287OAz306nwoEZnnihlm/g3gSeBiMztoZh8ufVnl5e7s2H+Kq0Ps3su4aW0byZTz1Ksnwi5FRCRU015Pwt3vKkchYdp/YpATA6OhnKCb6+pVC2moibFl33H+7eVLwy5HRCQ06uIDtr6e3luJQkDVVFdx3epWtmighIjMcwoo4HvPH2Vlaz0/s6Qx7FIAuGltO68dH+DgqcGwSxERCc28D6hTA6M8se84t1+xvOwTxE7m5rVtANqLEpF5bd4H1L+8cJTxlPOu9cvCLuWMtYsbWdJUy+Mabi4i89i8D6iHnzvCBYsauHx5U9ilnGFm3LimjZ/uO04q5WGXIyISinkdUCcHRvnpqye4/Yplkeney7h5bRunBsfYc1iX3xCR+WleB9QP9xwlmXJuj1D3XsaNa9LHoR7fV3nTRomIzIZ5HVD//NwRVrcluGxZdLr3MhYvqOOSpQv48d6usEsREQnFvA2oE/0j/PTV45Hs3st4z9Ur2Lb/FLsPnQ67FBGRspu3AfWDPUdJOZHs3st437UrSdTEuHfL62GXIiJSdvM2oP75uSNc2J7gkqULwi5lUk11ce7cuJL/+9xhunqHwy5HRKSs5mVAdfeN8NRrJ3hXhLv3Mn79xk7GU87XntwfdikiImU1LwPq+7uPBN17y8MuZVoXLErw9kuX8PWt+xkeS4ZdjohI2cy7gOofGedvf7yPK1c0R2buvel8+KbVnBoc46Edh8IuRUSkbOZdQP3Nv75Md/8In71jXeS79zKuXd3Kuo4m7n3iddw1s4SIzA/zKqBeOtrHvU+8wfuvWcmGlS1hl1M0M+PDN61mX1c//+9lnbgrIvPDvAkod+cPvrObBXXV/PYvXBJ2OTN2+xXLWbygli8//pr2okRkXpg3AfWdnYd5+vWT/O5tl9CaqAm7nBmrqf7HJ9kAAAgNSURBVK7iI//mIp7Yd4L/9dirYZcjIlJy017yfS7oHR7jv//zXq5c2cL7Nq4Mu5xz9hs3drL70Gn+4ocvsWJhPXds6Ai7JBGRkpnzATWWTPGZ7+zhxMAIX/nQNVRVVcbAiELMjD99zxUc7hnit//pOZY21XHdhYvCLktEpCTmdBdfV+8wv/rlrTz07CE+dstarljRHHZJ5622OsbmD2xkZWs9m/5hO69294ddkohISczZgNr62gne+T+28Pyh03zhfRv45Nt/JuySZk1zQ5yv/vq1xGPGB+95mh/tPaaBEyIy5xQVUGZ2m5m9ZGb7zOz3Sl3U+dh/YoDPP/Iy//7vt9JUV823/9ONvPuquXesZmVrA/fcfQ1VVfDh+7bxi3+7hX/Zc1RBJSJzhk33gWZmMeBl4O3AQeAZ4C53f2Gy39m4caNv27ZtNussaHQ8RXf/CAdODvKTl7v5173HePlYusvrnVcs5c/es54FdfGS1xGmsWSKbz97iL99dB/7TwxyydIF3Ly2jUuWNnHJsgWsWdxIbXUs7DIrlpltd/eNYdchMh8VM0jiWmCfu78GYGb3A3cAkwbUZMaTKT71T7uy1k3MRz+zznGHVHA7nnJGkylGx5OMjKcYGk3S3TfCiYHRs/+QKuPa1a28/5pV3HrpElYtaphpeRUpHqvizo0r+aWrOvjursPc9+R+7ntyP6PjKQBiVUZLfZzGumoaa9M/dfEYsSojVmVUVxlVZpD+78zsGrlDSSpk0o0Zu/GiNt57TeWO7BSZy4oJqA7gwITlg8B1uY3MbBOwCWDVqlUFn8iBXQd68tZPnHIoc6+qyjCgyoyqKqOmuora6ioaa6tpa6zl6gsWsmRBHUuaalnSVMfVqxbS3DC395amUh2r4pevXsEvX72C8WSKN04M8uLRXl462sfJgVH6R8bpHx6nb3icnqExkqkU40kn5U4y5ekvB8E3hNx96rncbbiqdX58kRGpRMUEVKHvznmfWO6+GdgM6S6+Qk8Uj1Xx2G//3IwKlJmrjlWxZnEjaxY38q71YVcjInJuihkkcRCY2AeyAjhcmnJERETSigmoZ4C1ZrbazGqA9wPfLW1ZIiIy303bxefu42b2UeCHQAy41933lLwyERGZ16YdZn5OT2rWDUx1jfI24Pisv/Dsq4Q6K6FGqNw6L3D39rCKEZnPShJQ076o2bZKOLekEuqshBpBdYrIzM3ZqY5ERKSyKaBERCSSwgqozSG97kxVQp2VUCOoThGZoVCOQYmIiExHXXwiIhJJJQ0oM7vXzLrMbPeEdZ8xs0NmtjP4eWcpayiGma00s0fNbK+Z7TGzjwfrW83sETN7JbhdGNE6I7VNzazOzJ42s11BnZ8N1q82s63B9vxmcOJ31Gr8qpm9PmFbbgirRpH5rqRdfGb2s0A/8DV3Xxes+wzQ7+5/WbIXniEzWwYsc/cdZrYA2A68G/gQcNLd/zS4DtZCd//dCNb5XiK0TS09+2/C3fvNLA5sAT4OfBJ4yN3vN7MvAbvc/YsRq/E3gYfd/Vth1CUiZ5V0D8rdfwKcLOVrzAZ3P+LuO4L7fcBe0rO43wHcFzS7j3QYhGaKOiPF0zLXoo8HPw7cAmQ++EPdnlPUKCIREdYxqI+a2XNBF2Co3Wa5zKwTuArYCixx9yOQDgdgcXiVZcupEyK2Tc0sZmY7gS7gEeBVoMfdx4MmBwk5XHNrdPfMtvyjYFv+tZnVhliiyLwWRkB9EbgI2AAcAf4qhBoKMrNG4EHgE+7eG3Y9kylQZ+S2qbsn3X0D6dnvrwUuLdSsvFXlvHhOjWa2Dvh94BLgGqAVCK1LV2S+K3tAufux4IMhBXyZ9IdX6ILjEA8CX3f3h4LVx4LjPpnjP11h1ZdRqM6oblMAd+8BHgOuB1rMLDNBcWQu2zKhxtuCblR39xHgK0RoW4rMN2UPqMwHfuCXgN2TtS2X4ID5PcBed//8hIe+C9wd3L8b+E65a5tosjqjtk3NrN3MWoL79cCtpI+XPQr8StAs1O05SY0vTvhCYqSPkYX+9ykyX5V6FN83gLeRniH6GPCHwfIG0t07bwAfyRznCYuZ3QQ8DjwPpILVnyZ9fOcBYBXwJnCnu4c26GOKOu8iQtvUzNaTHgQRI/0l6AF3/5yZXQjcT7rr7Fng14I9lSjV+GOgnfSVpHcCvzlhMIWIlJFmkhARkUjSTBIiIhJJCigREYkkBZSIiESSAkpERCJJASUiIpGkgIowM3vDzNrOt42ISCVSQImISCQpoCLCzL5tZtuDaxNtynms08xeNLP7gklMv2VmDROafMzMdpjZ82Z2SfA715rZT83s2eD24rL+g0REzpMCKjp+w93fAmwEfsvMFuU8fjGw2d3XA73Af5zw2HF3v5r0pLH/OVj3IvCz7n4V8F+BPy5p9SIis0wBFR2/ZWa7gKeAlcDanMcPuPsTwf1/BG6a8FhmctvtQGdwvxn4J0tfzfivgctLUbSISKkooCLAzN5GerLSG9z9StLz1NXlNMudk2ricmY+uySQmS38vwGPBlcy/sUCzyciEmkKqGhoBk65+2BwDOn6Am1WmdkNwf27SF+ifLrnPBTc/9CsVCkiUkYKqGj4AVBtZs+R3vN5qkCbvcDdQZtW0sebpvLnwJ+Y2ROkZ+wWEakoms28AgSXd3846K4TEZkXtAclIiKRpD0oERGJJO1BiYhIJCmgREQkkhRQIiISSQooERGJJAWUiIhEkgJKREQi6f8Dr/KZ84z08UQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot2(np.exp(log_Ppdf-log_Ppdf.max()),betas,alphas,xt='alpha',yt='beta')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the fit from the maximum a posteriori (MAP) parameters. We'll get the index of the maximum using numpy.argmax(), then use numpy.unravel_index() to convert this 1D index into 2D indices in the posterior array. Using those indices, we'll get the MAP values of the parameters and plot them with the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAP parameters:  2.1 21.5\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11ed83350>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXwU9f3H8deHABK5AnLIIeKJFwoY0Qoiioq3ERW8qdqitV61xYq2AlYFi60VtQqeeIGKiNbWKgQQtcghIHghooAEhHCEQwK5vr8/ZvAXYBM22Z2dze77+Xjkkd3Z2cyHYfPO7He/8xlzziEiIumjVtgFiIhIYin4RUTSjIJfRCTNKPhFRNKMgl9EJM3UDruAaDRr1sy1b98+7DJERGqUTz/9dK1zrvmuy2tE8Ldv3545c+aEXYaISI1iZssiLddQj4hImlHwi4ikmcCC38w6mNn8cl+bzOw2M2tqZpPMbLH/vUlQNYiIyO4CC37n3CLnXCfnXCfgWGAr8CZwJ5DrnDsEyPXvi4hIgiRqqKcXsMQ5twy4ABjjLx8D5CSoBhERIXHBfykw1r/d0jm3CsD/3iJBNYiICAkIfjOrC5wPvF7F5w0wszlmNic/Pz+Y4kRE0lAijvjPAuY651b791ebWSsA//uaSE9yzo12zmU757KbN9/t/AMREammRAT/Zfz/MA/A20B//3Z/4K0E1CAiUqP0GzWDfqNmBPKzAw1+M9sbOB2YUG7xcOB0M1vsPzY8yBpERGRngbZscM5tBfbZZdk6vFk+IiISAp25KyKSZhT8IiJpRsEvIpJmFPwiImlGwS8ikmYU/CIiaUbBLyKSZhT8IiJJZuK8POYtL2Dm9+vpNnwKE+flxfXnK/hFRAJQ3ZYLE+flMWjCQopKywDIKyhk0ISFcQ1/Bb+ISBIZ8d4iCotLd1pWWFzKiPcWxW0bCn4RkQiCbJJWmZUFhVVaXh0KfhGRJNI6K7NKy6tDwS8ikkQG9u5AZp2MnZZl1slgYO8OcdtGoN05RUSkanI6twHgjvELKCoto01WJgN7d/h5eTwo+EVEkkxO5zaMnbUcgFev/0Xcf76GekRE0oyCX0QkzSj4RUTSjIJfRCTOgm65ECsFv4jILmIJ7kS0XIiVgl9EpJxYgzsRLRdipeAXESkn1uBORMuFWGkev4hIObEGd+usTPIirFvVlgtBzN/fQUf8IiLlxNorJxEtF2Kl4BcRKSfW4M7p3IZhfTpSN8OL1zZZmQzr0zGuLRdipaEeEZFy4tErJ+iWC7FS8IuI7CLZgztWGuoREUkzCn4RkTSj4BcRSUbLZsCrV0Jx/Of/a4xfRCSZrPkaJg+Bb96Fhq1g3bewb8e4bkLBLyKSDDbmwbQHYP4rULcB9LoHjv8N1N077ptS8IuIhKmwAD7+B3zyBLgyL+xP+j3U3yewTQYa/GaWBTwNHAU44FpgEfAq0B5YCvR1zm0Isg4RkUTb4zTQku0w6yn48CEo3AAd+8Kpf4Im+wdeW9BH/I8A/3XOXWxmdYG9gbuAXOfccDO7E7gT+GPAdYiIVElg8/fLymDhazDlfti4HA46FU4bAq2OCWZ7EQQW/GbWCOgB/BLAOVcEFJnZBUBPf7UxwDQU/CKS6pyDb3O9D25XL/SC/vyRcNApCS8lyCP+A4F84DkzOwb4FLgVaOmcWwXgnFtlZi0iPdnMBgADANq1axdgmSIiAVs5DybdA99Ph6z94aJn4Mg+UCucGfVBBn9toAtws3Nuppk9gjesExXn3GhgNEB2drYLpkQRkQCt/w6m3AefvwF77wNnPgjZ10DtvUItK8jgXwGscM7N9O+Pxwv+1WbWyj/abwWsCbAGEZHE25IP00fAnGchow70GAgn3gL1GoVdGRBg8DvnfjSzH8ysg3NuEdAL+NL/6g8M97+/FVQNIiIJtX0LfPJP+PgR74zbLldBz0HQcN+wK9tJ0LN6bgZe9mf0fAdcg9cm4jUzuw5YDlwScA0iIsEqLYa5L8C04fDTGjjsXOg1GJofGnZlEQUa/M65+UB2hId6BbldEZGEcA6+ehty7/VaK+x3AvR7CdodH3ZlldKZuyKSkvqNmgEEOB9/6cfeTJ28OdCsA1w6FjqcBWbBbC+OFPwiIlWx+kvIHQrf/Ndronb+o3DM5ZBRc+K05lQqIhKmjStg6jD47BWo29Abwz/+hkCaqAVNwS8iUpnCDfDRwzBzlNdE7YQbvSZqezcNu7JqU/CLiERSvA1mPwXTH4JtG+HofnDKXQlpohY0Bb+ISHllpbDgNZh6P2z8AQ4+zWuiFueLoYRJwS8iAn4TtckwaTCs+QJadYILHoMDe4ZdWdwp+EVE8j71An/ph9CkPVz8LBxxYWhN1IKm4BeR9LVuiXfy1ZcTYe9mcNYIOPaXULtu2JUFSsEvIkkplhOwJs7LY97yAopKy+g2fAoDe3cgp3Ob/19hyxr44EH49HnIqAs97oATb06aJmpBU/CLSEqZOC+PQRMWUlRaBkBeQSGDJiwEIOeIRjDjcfjfo14TtWP7w8l3QsOWYZaccAp+EUkpI95bRGFx6U7LCotLGfH2HHIm/x5+yofDz/dOwGp2cEhVhkvBLyIpZWVBYeTlhbWg3aFeT539jktwVclFwS8iKaV1ViZ5EcK/df1a8Mt/14gmakFLzblKIpK2Bp5Qn0wr2WlZZp1aDDy3k0Lfp+AXkdRQ8AO8+RtypvVmWL0XyKrlHfW3ycpkWJ+jd57Vk+Y01CMiNVvhBvjwbzBztHf/xJvI6X47Y19cBATYj78GU/CLSM1UXAizRnuhv20THHOZ10Qta7+wK0t6Cn4RqVnKSuGzcTD1Adi0Ag45w5uaue9RYVdWYyj4RSTpRDzztlNrWPw+TB4Ca76E1l3gwifhgJPCLrfGUfCLSFKJeObtG/Phg7+SU/A8ND0QLnkejsjRLJ1qUvCLSFKJeOZtCYxYk01OzlFeE7WMOuEUlyIU/CKSVCo887asKXQ9J8HVpCbN4xeR5LF9M63rFUV8qHVWZoKLSV0KfhEJX0kRzHoKRnZmYMnoCGfeZjCwd4eQiks9GuoRkfA4B1+86V0MZcP3sH93ci67F/L35Y7xCygqLaNNVubu/fSjoBO3KqbgF5FwfD8dJt0DK+dBiyPg8tfhkNPBjJy2MHbWckABHgQFv4gk1o+fe3Pxv50EjdpCzhNwdD+olRF2ZWlDwS8iiVHwA0y93zvrtl5jOP0v0HUA1KkXdmVpR8EvIsHaut7rpzPrKe/+iTfDSbdDZpNw60pjCn4RCUZxIcx8Ej58GLZvgk6Xe03UGrcNu7K0p+AXkfgqK4X5r3hN1DavhEN6w2mDoeWRYVcmvkCD38yWApuBUqDEOZdtZk2BV4H2wFKgr3NuQ5B1iEgCOAff/BcmD4X8r6DNsdBntJqoJaFEnMB1inOuk3Mu279/J5DrnDsEyPXvi0hN9sNseO5sGHsplBbBJWPgV7kK/SQVxlDPBUBP//YYYBrwxxDqEJFYrV0MuUPhq39B/RZwzt+gS/+4NFHT/P3gBB38DnjfzBwwyjk3GmjpnFsF4JxbZWYtIj3RzAYAAwDatWsXcJkiUiWbf4Rpw2HuC1AnE3reBb/4LezVIOzKJApBB38359xKP9wnmdnX0T7R/yMxGiA7O9sFVaCIVMG2TfC/kTDjcW9IJ/taOPkOaBDx+E2SVKDB75xb6X9fY2ZvAl2B1WbWyj/abwWsCbIGEYmDkiL49Dn44EHYug6OvBBO/TPsc1DYlUk1BBb8ZlYfqOWc2+zfPgO4F3gb6A8M97+/FVQNIhKjsjL4YgJM+QtsWArtT4LTh3ozdvag36gZgMbqk1GQR/wtgTfNuzRabeAV59x/zWw28JqZXQcsBy4JsAaRtBVz8H43DSYNhlXzoeVRcMUbcHAvXe4wBQQW/M6574BjIixfB/QKarsiEqNVC7wmaktyofF+kPMkHN1XTdRSiM7cFRHPhmVeE7UFr3lN1M64D477tZqopSAFv0i627oepj8Es58CqwXdboXuv4PMrLArk4Ao+EXSVdFWr4naR/+Aos1eE7Wed0Hjql3pSmoeBb9Iuiktgc9eganDvCZqh57lNVFrcXjYlUmCKPhF0oVzsOhdr8VC/tfQ9ji46Glo3y3syiTBFPwi6eCHWd71bZfPgH0Ohr4vwuHnaWpmmlLwi6SwViU/wLjH4Ot3oEFLOPdh6HxVXJqoVWbivDzmLS+gqLSMbsOnMLB3B3I667ODZKHgF0lBE2d8weJly5lVVp/ZK3ozsGMPci6+EurWD37b8/IYNGEhRaVlAOQVFDJowkIAhX+SSEQ/fhFJlG2bmPjS4wx6axHryxrgMPLcPgz6an8mflGQkBJGvLeIwuLSnZYVFpcy4r1FCdm+7JmCXyQVlBTBJ0/CyE6M+Lw+hey108OJDN6VBYVVWi6Jp6EekZpsRxO13HuhYBkccDIr1zePuGqigrd1ViZ5EbbVOiszIduXPdMRv0hNtWQqPNUT3rgO6jWCKyfA1W9VGLCJCt6BvTuQWWfnvj6ZdTIY2LtDQrYve7bH4Dezm8ysSSKKEZEorFoAL14IL+ZA4Qbo8xQMmP5z58x4BW+/UTN+7vBZFTmd2zCsT0fqZnjx0iYrk2F9OuqD3SQSzVDPvsBsM5sLPAu855zTFbFEEm3DMphyHyx8DTKbQO8H4LhfQe2dx/N3BOwd4xdQVFpGm6zMhE+nzOnchrGzlgPqx5+M9hj8zrk/mdmf8S6kcg3wmJm9BjzjnFsSdIEiae+ndfDhQzD7abAM6H47dL/N66BZAQWvVCaqD3edc87MfgR+BEqAJsB4M5vknLsjyAJF0lbRVpj5hN9EbQt0vhJ6DoJGrcOuTGq4PQa/md2Cd4nEtcDTwEDnXLGZ1QIWAwp+kXgqLYH5L8O0YbB5FXQ4B3rdAy0OC7sySRHRHPE3A/o455aVX+icKzOzc4MpSyQNOQeL/gOTh8LaRbDf8XDxc7C/hmokvvY4q8c5d8+uoV/usa/iX5JIGlr+CTx7Joy7HHDQ72X6FQ+h33/DLkxSkU7gEglT/iLvCH/Rv6HBvnDeI9DpSsioDdOrPpVSJBoKfpEwbFrljeHPexHq1IdT/wQn3JiQJmoiCn6RRNq2ET5+BGb8E8pKoOv10GMg1N8n7MriTtNIk5eCXyQRSrbD7Gdg+ggoXA8dL4FT7oamB4RdWUTqp5/aFPwiQSorg8/Hw5S/QMFyOLAnnDYUWncKfNPVPeJWP/3UpyZtIkFZMgVGnwwTfg31suCqN+HqtxIS+rFQP/3UpyN+kXhbOR8mD4bvpkFWO6+J2lEXQ62acZylfvqpT8EvUoEdnSmjHjJZ/73XRO3z8ZDZFHoPg+Ou262JWrJTP/3Up+AXidVPa70PbWc/A7Vqw0m/h263VtpELZkN7N2BQRMW7jTco376qUXBL1JdRT950zI/fgSKf4LOV/lN1FrF/KPDnFWTDG2dJVgKfpGqKi2BeS/AtOGwZTUcdq7XRK35zkfEVR4q8iXDrBq1dU5tCn6RaDkHX7/jtVhYt9hrotb3BWh3Qlw3U9msGh11Szwo+EWisWwGTLoHVsyCZofCpa9Ah7PBLO6b0qwaCVrgwW9mGcAcIM85d66ZHQCMA5oCc4GrnHNFQdchUh1tipfB2JFeu+SGreC8kdDpCq+JWkA0q0aCloiJxbcC5ds3Pwg87Jw7BNgAXJeAGkSqZmMe1xc8zENrfwNLP4JT/ww3z4Vj+wca+kDcLpYuUpFAg9/M2gLn4F25CzMz4FRgvL/KGCAnyBpEqqSwACYPgUe70KMwl3frXwC3zIcef4C6eyekhJzObRjWpyN1M7xfzzZZmQzr01Hj+xI3QQ/1/APv0owN/fv7AAXOuRL//gpAr2YJX8l2mPWUd1Hzwg3QsS+3/Xg2+bX35ZwQOmdqVo0EKbAjfv+yjGucc5+WXxxhVVfB8weY2Rwzm5Ofnx9IjSKUlcFnr8Kj2fD+3dC6M1w/nYkHDuH9lfWY+f16ug2fwsR5eWFXKhI3QR7xdwPON7OzgXpAI7x3AFlmVts/6m8LrIz0ZOfcaGA0QHZ2dsQ/DiLV5hwsyYVJQ2D1Qtj3aDh/JBx0SlLMo08GeqeRugI74nfODXLOtXXOtQcuBaY4564ApgIX+6v1B94KqgYJV79RM34+iSmprJwHL5wPL10E2zfBRc/AgA/goFOA+HSn3HHmrd4xSDIKYx7/H4FxZnYfMA94JoQaJB2t/85vovYG7L0PnPkgZF+zWxO1WOfR6x2DJLuEBL9zbhowzb/9HdA1EdsVAbwmah/8FeY8Cxl1vEsdnngL1GsUcfVY59HrzFtJdjpzV5JWdXvd/KzoJ5jxOHw8Eoq3Qperoeed0HDfSp8Wa3dKnXkryU7BL6mntBjmvgAfPFiuidpgaH5oVE+PtTulzryVZKfgl9ThHHz1L8gdCuu+hXa/gH4vwX5VH1mMZR59vPrZa1aNBEXBL6lh2f/8JmqzoflhcNk4OPTMQJqo7Yn62UuyU/BLzbbmK69N8jfvek3Uzn8Ujrk88H46e6IzbyWZKfilZtqYB9MegPmvQN2G3hj+8TckrJ+OSE2m4JeapbAAPnoYZj4JrgxOuNG7xu3eTcOuTKTGUPBLIOJ+zdjibTD7KZj+EGzbCEf3hVPuhib7V/iUmKeDiqQoBb/EXTzPXDVXCvPHwtT7YeMPcFAvOG0ItDo6zlWLpI9EXIhF0kxcet3MXcEXy1Yzc2kB3cZtZ2Jpd7j6LbhqgkJfJEY64pe4i7nXzZQPGTRpHYWuDgB5NGfQxgtg4yEJvWqPhogkVemIX+KuojNU93jm6rol8PovGfH+kp9Df4fC4rIqvWMQkYop+CXuqnzN2C358O8/wONd4Zv3WEmziKup141IfGioR+Iu6jNXt2+BGY/B/x6F4kLvQuYn/5HWj3+ZEr1uNFQkyUrBL4Go9MzV0mKYOwamPQg/rYHDz4de90CzQwAY2Ls0Lr1uRCQyBX8SS7l56M7Bl29B7r2wfgm0OxEufQX2O26n1dTrRiRYCn5JjKUfeU3U8j6F5ofDZa/Cob0rbKIWa6+buJ9AJpJCFPwSqP2Kv4eXH4bF70GjNnDB43DMZVArY89PriZd+lCkcgp+CUbBD/ym4G/0KJwMmxvBaUPh+OuhTvAf0OrShyKVU/BLfBVugA//DjNH0a20jH/X78N5v30ooU3UdOlDkcop+CU+igth1mj48G+wbRMccym3rTyTtbVbcl6CO2fq0ocildMJXBKbslKY9zI8mu19eNu2K9zwIVz4JGtrtwylpCqfQCaSZnTEL9XjHCx+HyYPgTVfQuvOcOETcECPn1cJaxqqpoOKVE7BL1W3Yg5MGgzLPoImB8DFz8GRF4ZyfduK6NKHIhVT8Ev01i2B3KHeSVh7N4OzH4Iu/aF23bArE5EqUPAnqaQ6AWnLGpg23GuzkLEXnHwnnHgT7NUw0M3qSF0kGAr+JJQsJyD1fyKXc396g0u2T4TS7XDsL6HHHdAwnA9tRSQ+FPxJKPQTkEqKYO4YHsm/j6yyAjgix2uits9BwW9bRAKn4E9CoZ2A5Bx88abXRG3D9+TV7ciIhkO4v+81wW5XRBJKwZ+EQjkB6fvp3jz8lfOgxRFw+evcO6VBUs3UEZH40AlcAeo3asbPrZWrIqEnIP34Obx0MYw5z7sSVs4TcMNHcOgZCn2RFKUj/iSUkBOQCn6AqffDZ+OgXiM4/V7oOiAhTdREJFwK/iQVjxOQIl7IZet6r5/OrKe8+yfeDCfdDplNYqo3GWk6qEhkgQW/mdUDpgN7+dsZ75wbbGYHAOOApsBc4CrnXFFQdYivuBBmjoKP/u41Uet0OZxyFzRuG3ZlIpJgQR7xbwdOdc5tMbM6wEdm9i5wO/Cwc26cmT0JXAc8EWAdac1cKcx7CaY+AJvy4JDecNpgaHlk2KWJSEgCC37nnAO2+Hfr+F8OOBW43F8+BhiCgj/+nKPLtplcvvlZeGsZtDkW+oyG9t3DrkxEQhborB4zyzCz+cAaYBKwBChwzpX4q6wAIn5iaWYDzGyOmc3Jz88PsszUs2IOPH8Of9wwmAxXApeMgV/lVin0d7SMmPn9eroNn8LEeXkBFiwiiRRo8DvnSp1znYC2QFfg8EirVfDc0c65bOdcdvPmzYMsM3Ws/RZevQqe7gVrF/N0o5v4Q/NRcGROlaZmVtQyQuEvkhoSMo/fOVcATANOALLMbMcQU1tgZVDbre48+hpn82p453fweFdYMgV63gW3zGNS/XMptaqP5lXWMkJEar4gZ/U0B4qdcwVmlgmcBjwITAUuxpvZ0x94K6gaUt72zfDxSJjxGJQWQfa1cPId0KBFTD9W16wVSW1BzuppBYwxswy8dxavOefeMbMvgXFmdh8wD3gmwBpCE4+2yhXOQy8pgk+fhw8ehK1rvYugnPrnuDVR0zVrRVJbkLN6FgCdIyz/Dm+8P2UF1la5rAy+fBNy/wIbvof2J8HpQ70ZOxFqqO4fnoG9OzBowsKdhnt0zVqR1KEzdwMQSFvl7z7wmqitmg8tjoQrxsPBp0X80DbWPzy6Zq1IalPwVyJiy4MoxHWM/MeF3vVtl+RC4/0g50k4ui/UyqjwKfH4w6Nr1oqkLgV/AOIyRr5hmXe27YJXoV5jOOM+OO7XUKfeHp+qD2dFpDJqyxyAmNoqb10P790Nj2V7F0XpdgvcOt9rphZF6EPFf2D04ayIgI74A1GtMfKirTDzSfjoH1C0GY65HE4ZVK0mavpwVkQqk7LBH4/plLGIeoy8tAQ+ewWmDoPNK+HQM6HXYGh5REzbBn04KyKRpWTwBzadMp6cg0XvQu5QyP8a2mTDRU9D+25x+fH6cFZEKpKSY/xJ33Lgh1nw3Fkw7jIoK4G+L8CvJsct9EVEKpOSR/xJO6tl7WKYPAS+fgfqt4Bz/g5droaMOuHWJSJpJSWDP+laDmz+EaYNg7kvete0PeVuOOFG2KtBOPVESUNEIqkpJYM/WWa1ZJb9xPk/vQ4j34bSYjjuV9BjIDRQm2kRCU9KBn88ZrXENCuoZDvMeZaR+cNoVLYRjroITv0TND2wOv8cEZG4Ssngh9hmtVR7VlBZGXz+Bkz5CxQsY1ndTrzS8FqGXdy/ev8IEZEApGzwx6JavW6WTIXJg2HVZ9CyI1z5Bh0P6sWwKlz5SkQkERT8EVRpVtCqz7yZOkumQON2cOFo6HgJ1ErJmbIikgIU/BFENStow1KYcj8sfA0ym0DvByD7uqj76SSCZuWISCQK/ggqnRX00zr48CGY/TRYLej+O+h2G2RmhVixiEj0FPwRRJwV1Ks9OZvHwshHoGgLdLoCTrkLGrUOuVoRkapR8Fdgx6ygWq6Uscd9C1N/A1t+hA5ne03UWhwWdokiItWi4K+Ic2Rvm8Flm5+Dfy2Htl3hkudhf42bi0jNpuCPZPknMGkwAzd8Ql5GW+j3Ehx2bsTr24qI1DQpHfxVntWSvwgmD4VF/4YGLRnd+BamZvZm7OHdgylQRCQEKR38Udu0ymuiNu9FqFPfa69wwo3kPrcg7MpEROIuvYN/20b4+BGY8U+vL37X66HHH6B+s7ArExEJTHoGf8l2mP0MTB8BhevhqIv9JmoHhF2ZiEjg0iv4y8rg8/F+E7XlcGBPOG0otO4UdmUiIgmTPsH/ba7XRO3HhbBvR7hyAhzcq9KnqOWBiKSi1A/+lfO9wP9uGmS1gz5PeUM7aqImImkqtYP/X7fCp89DZlPoPQyOuw5q7xV2VSIioUrt4G9yAHS/HbrfBvUah12NiEhSSO3g735b2BWIiCQdDXSLiKSZwILfzPYzs6lm9pWZfWFmt/rLm5rZJDNb7H9vElQNIiKyuyCP+EuA3zvnDgdOAH5rZkcAdwK5zrlDgFz/voiIJEhgwe+cW+Wcm+vf3gx8BbQBLgDG+KuNAXKCqkFERHaXkDF+M2sPdAZmAi2dc6vA++MAtKjgOQPMbI6ZzcnPz09EmSIiaSHw4DezBsAbwG3OuU3RPs85N9o5l+2cy27evHlwBYqIpJlAg9/M6uCF/svOuQn+4tVm1sp/vBWwJsgaRERkZ0HO6jHgGeAr59zfyz30NtDfv90feCuoGkREZHfmnAvmB5t1Bz4EFgJl/uK78Mb5XwPaAcuBS5xz6/fws/KBZdUspRmwtprPTQTVFxvVFxvVF5tkr29/59xuY+WBBX+yMLM5zrnssOuoiOqLjeqLjeqLTbLXVxGduSsikmYU/CIiaSYdgn902AXsgeqLjeqLjeqLTbLXF1HKj/GLiMjO0uGIX0REylHwi4ikmZQJfjM708wWmdm3ZrZbx08z28vMXvUfn+n3D0pUbRFbVO+yTk8z22hm8/2vexJVn7/9pWa20N/2nAiPm5mN9PffAjPrksDaOpTbL/PNbJOZ3bbLOgndf2b2rJmtMbPPyy2LquW4mfX311lsZv0jrRNQfSPM7Gv//+9NM8uq4LmVvhYCrG+ImeWV+z88u4LnVvq7HmB9r5arbamZza/guYHvv5g552r8F5ABLAEOBOoCnwFH7LLOjcCT/u1LgVcTWF8roIt/uyHwTYT6egLvhLgPlwLNKnn8bOBdwPDabM8M8f/6R7wTU0Lbf0APoAvwebllfwXu9G/fCTwY4XlNge/87038200SVN8ZQG3/9oOR6ovmtRBgfUOAP0Tx/1/p73pQ9e3y+N+Ae8Laf7F+pcoRf1fgW+fcd865ImAcXvvn8sq3gx4P9PLbSgTOVdyiuia5AHjBeT4Bsnb0XEqwXsAS51x1z+SOC+fcdGDXM86jaTneG5jknFvvnNsATALOTER9zrn3nXMl/t1PgLbx3m60Kth/0Yjmdz1mldXn50ZfYGy8t5soqRL8bYAfyt1fwe7B+vM6/ot/I7BPQqorZ5cW1bv6hZl9ZmbvmtmRCS0MHPC+mX1qZgMiPB7NPk6ES6n4Fy7M/QfRtRxPlv14Ld47uEj29FoI0k3+UNSzFQyVJcP+OwlY7ZxbXMHjYe6/qKRK8Ec6ct91nmo06wRqDy2q5+INXxwDPApMTGRtQDfnXBfgLLyrpfXY5fFk2ARTLP8AAALfSURBVH91gfOB1yM8HPb+i1Yy7Me78a6Q93IFq+zptRCUJ4CDgE7AKrzhlF2Fvv+Ay6j8aD+s/Re1VAn+FcB+5e63BVZWtI6Z1QYaU723mtVikVtU/8w5t8k5t8W//R+gjpk1S1R9zrmV/vc1wJt4b6nLi2YfB+0sYK5zbvWuD4S9/3zRtBwPdT/6HyafC1zh/AHpXUXxWgiEc261c67UOVcGPFXBdsPef7WBPsCrFa0T1v6rilQJ/tnAIWZ2gH9UeCle++fyyreDvhiYUtELP978McFILarLr7Pvjs8czKwr3v/NugTVV9/MGu64jfch4Oe7rPY2cLU/u+cEYOOOYY0EqvBIK8z9V040LcffA84wsyb+UMYZ/rLAmdmZwB+B851zWytYJ5rXQlD1lf/M6MIKthvN73qQTgO+ds6tiPRgmPuvSsL+dDleX3izTr7B+8T/bn/ZvXgvcoB6eEME3wKzgAMTWFt3vLejC4D5/tfZwA3ADf46NwFf4M1S+AQ4MYH1Hehv9zO/hh37r3x9Bjzu79+FQHaC/3/3xgvyxuWWhbb/8P4ArQKK8Y5Cr8P7zCgXWOx/b+qvmw08Xe651/qvw2+BaxJY37d44+M7XoM7Zrm1Bv5T2WshQfW96L+2FuCFeatd6/Pv7/a7noj6/OXP73jNlVs34fsv1i+1bBARSTOpMtQjIiJRUvCLiKQZBb+ISJpR8IuIpBkFv4hImlHwi4ikGQW/iEiaUfCLVIOZHec3E6vnn635hZkdFXZdItHQCVwi1WRm9+GdEZ4JrHDODQu5JJGoKPhFqsnvFTMb2IbXIqI05JJEoqKhHpHqawo0wLuqWr2QaxGJmo74RarJzN7GuwLUAXgNxW4KuSSRqNQOuwCRmsjMrgZKnHOvmFkG8D8zO9U5NyXs2kT2REf8IiJpRmP8IiJpRsEvIpJmFPwiImlGwS8ikmYU/CIiaUbBLyKSZhT8IiJp5v8Agw+GepU3J8AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "imax = np.argmax(log_Ppdf)\n",
    "i=np.unravel_index(imax,log_Ppdf.shape)\n",
    "best_alpha = alphas[i[0]]\n",
    "best_beta = betas[i[1]]\n",
    "print('MAP parameters: ',best_alpha,best_beta)\n",
    "plt.errorbar(xdata, ydata, yerr=sig, marker='o',ls='none')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.plot(xdata,np.polyval([best_alpha,best_beta],xdata))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To do the model comparison, we need to integrate the posterior. But we may be dealing with very small numbers, so can't just take the exponential of the log. We have to be careful about just normalizing it, e.g. to the maximum, because the value matters. However, so long as we normalize all models that we consider with the same normalizing factor, we will be able to take the ratio of the integrals. So we'll get a normalization from this model, and use it for subsquent models.\n",
    "<br>\n",
    "Similarly, since we'll be comparing models, we need to do the integrals paying attention to the normalization; we will still use the rectangle rule (in 2D), but we will multiply by the bin widths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evidence from linear model:  0.5005677021697865\n"
     ]
    }
   ],
   "source": [
    "norm=log_Ppdf.max()    # we will use this normalization constant for both this and the quadratic model\n",
    "total_Ppdf = np.exp(log_Ppdf-norm).sum()*dalpha*dbeta   # integrate posterior using simple rectangle rule in 2D\n",
    "evidence2=total_Ppdf   #  the integral of the posterior is the evidence\n",
    "print('Evidence from linear model: ', total_Ppdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What do you think about the quality of your fit? From the Bayesian perspective, can you \"accept\" or \"reject\" it? Why or why not?\n",
    "<br>\n",
    "<B> ANSWER HERE: </B>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Quadratic model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, now we will do a quadratic model (3 parameters):\n",
    "$$y = \\alpha x^2 + \\beta x + \\gamma$$\n",
    "Again, let's write the expressions for prior, likelihood and posterior: note that if you're using numpy.polyval() to calculate the model, these routines will be identical: the order of the polynomial is just determined by the number of parameters passed to numpy.polyval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# again, choose trial grid limits and spacing, may need to iterate. Three parameters\n",
    "dalpha=0.005    # spacing in alpha (quadratic term)\n",
    "alphas=np.arange(-0.1,0.2,dalpha)  # trial grid for alpha\n",
    "dbeta=0.04      # spacing in beta (linear term)\n",
    "betas=np.arange(-2,2.5,dbeta)   # trial grid for beta\n",
    "dgamma=0.5      # spacing in gamma (constant term)\n",
    "gammas=np.arange(15,40,dgamma)    # trial grid for gamma\n",
    "\n",
    "def log_prior(theta):\n",
    "    \"\"\" Return log(prior)\n",
    "            For a constant uniform prior, need this to integrate to unity over the trial grid\n",
    "    \"\"\"\n",
    "    return np.log(1/len(alphas)/len(betas)/len(gammas))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_likelihood(theta, x, y, sigma):\n",
    "    \"\"\" calculate log(likelihood) given polynomial parameters theta, and data points x,y,sigma\n",
    "    \"\"\"\n",
    "    y_model = np.polyval(theta, x)\n",
    "    n=len(y)\n",
    "    return np.sum(-0.5*(y-y_model)**2/sigma**2) - n*np.log(sigma) - n*np.log(2*np.pi)\n",
    "\n",
    "def log_posterior(theta, x, y, sigma):\n",
    "    \"\"\" calculate log_posterior from log_prior and log_likelihood\n",
    "    \"\"\"\n",
    "    return log_prior(theta) + log_likelihood(theta, x, y, sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a 3D grid in likelihood to compute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate log(posterior PDF)\n",
    "log_Ppdf=np.zeros([len(alphas),len(betas),len(gammas)])\n",
    "for i,alpha in enumerate(alphas) :\n",
    "    for j,beta in enumerate(betas) :\n",
    "        for k,gamma in enumerate(gammas) :\n",
    "            log_Ppdf[i,j,k] = log_posterior(np.array([alpha,beta,gamma]),xdata,ydata,sig)   #use your function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's display it using the plot3() routine from last time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot3(grid,xgrid,ygrid,zgrid,xt=None,yt=None,zt=None) :\n",
    "     \"\"\" Plot projections of 3D PDF\n",
    "     \"\"\"\n",
    "     fig,ax=plt.subplots(3,3,figsize=(12,8))\n",
    "     # marginalize over z (background)\n",
    "     tmp=np.sum(grid,axis=0)\n",
    "     ax[1,0].imshow(tmp,interpolation='nearest',origin='lower',aspect='auto',\n",
    "                    extent=[xgrid[0],xgrid[-1],ygrid[0],ygrid[-1]])\n",
    "     ax[1,0].set_ylabel(yt)\n",
    "     ax[0,0].plot(xgrid,np.sum(tmp,axis=0))\n",
    "     ax[0,0].set_ylabel(xt)\n",
    "     ax[1,1].plot(ygrid,np.sum(tmp,axis=1))\n",
    "    \n",
    "     # marginalize over y (sigma)\n",
    "     tmp=np.sum(grid,axis=1)\n",
    "     ax[2,0].imshow(tmp,interpolation='nearest',origin='lower',aspect='auto',\n",
    "                    extent=[xgrid[0],xgrid[-1],zgrid[0],zgrid[-1]])\n",
    "     ax[2,0].set_xlabel(xt)\n",
    "     ax[2,0].set_ylabel(zt)\n",
    "     ax[2,2].plot(zgrid,np.sum(tmp,axis=1))\n",
    "     ax[2,2].set_xlabel(zt)\n",
    "    \n",
    "     # marginalize over x(scale)\n",
    "     tmp=np.sum(grid,axis=2)\n",
    "     ax[2,1].imshow(tmp,interpolation='nearest',origin='lower',aspect='auto',\n",
    "                    extent=[ygrid[0],ygrid[-1],zgrid[0],zgrid[-1]])\n",
    "     ax[2,1].set_xlabel(yt)\n",
    "     ax[1,2].remove()\n",
    "     ax[0,2].remove()\n",
    "     ax[0,1].remove()\n",
    "     fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the posterior: again, if you don't have a well-centered, well-sampled, peak, adjust the grids above and rerun. We will need to integrate this, so you don't want significant probability to fall off of your grid!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA14AAAI4CAYAAAB6AjgBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXxcZ3n3/+812i1buyzb8iKviZckduI4CVkIJIQk5GmAwvOEUkgpNLQNLX3gaQvtrz+6/OhDF6DQAm0gKaGlhBQopBDIvq92EideE8u7LNtabEu2Za1z/f6YM0HSGduyNTNnRvq8X695Seee+5xzaawoc81939dt7i4AAAAAQObEog4AAAAAACY6Ei8AAAAAyDASLwAAAADIMBIvAAAAAMgwEi8AAAAAyLDCqAPIprq6Om9qaoo6DABj8NJLL3W4e33UcQAAAKTDpEq8mpqatG7duqjDADAGZrY76hgAAADShamGAAAAAJBhJF4AAAAAkGEkXgAAAACQYSReAAAAAJBhJF7IS6+1HNHB7t6owwAAAADGhMQLeefZ5g695+vP6r1ff1ZtJF8AAADIAyReyCu7Oo7rd777subWTNHhnn79xr+u1dHegajDAgAAAE6JxAt5o7t3QB/7zjrFTLr7I2v09Q9eqNcPHtXv/PvL6h+MRx0eAAAAcFJZTbzMbI6ZPWZmW8xsk5l9Mmj/czPbZ2brg8eNw875rJk1m9nrZvbOYe3XB23NZvaZbP4cyL6huOv3v/eKdnUc19c/eJHm1k7R1edM1xfee56ebu7QZ374mtw96jABAACAlAqzfL9BSZ9295fNbJqkl8zsoeC5L7v73w/vbGbLJN0iabmkWZIeNrMlwdNfk/QOSS2S1prZfe6+OSs/BbLuCz/fosdfb9fn37NCly2sfbP9/avn6GB3r/7+wTc0s6pUf/jOcyOMEgAAAEgtqyNe7r7f3V8Ovj8qaYukxlOccrOke9y9z913SmqWtCZ4NLv7Dnfvl3RP0BcT0PM7OvXNp3bq1svm6YOXzAs9f/vbFul9F83W1x/frrajFNsAAABA7olsjZeZNUlaJemFoOkTZvaamd1lZtVBW6OkvcNOawnaTtae6j63mdk6M1vX3t6exp8A2fKz1/arrKhAn71xacrnzUwfu3K+3KVHtrRlOToAAADg9CJJvMxsqqQfSvoDd++W9A1JCyWtlLRf0heTXVOc7qdoDze63+Huq919dX19/bhjR3a5ux7eclBXLq5TaVHBSfud0zBNc2rK9NDmg1mMDgAAABibrCdeZlakRNL1XXf/kSS5+0F3H3L3uKRvKjGVUEqMZM0ZdvpsSa2naMcEs6m1W/u7evWOZQ2n7Gdmum7ZDD3d3KFjfYNZig4AAAAYm2xXNTRJd0ra4u5fGtY+c1i390jaGHx/n6RbzKzEzOZLWizpRUlrJS02s/lmVqxEAY77svEzILse3nJQZtLbz51+2r7vWNag/sG4nnyDKaUAAADILdmuani5pA9J2mBm64O2P5H0ATNbqcR0wV2SPi5J7r7JzO6VtFmJioi3u/uQJJnZJyQ9IKlA0l3uvimbPwiy4+EtB3XR3GrVTi05bd/V86pVPaVID20+qBvPm3na/gAAAEC2ZDXxcvenlXp91v2nOOfzkj6fov3+U52H/Le/64Q27uvWZ24YW4n4woKYrlnaoAc3HdDAUFxFBewPDgAAgNzAO1PkrIeDCoXXLj31+q7h3rGsQd29g1q781CmwgIAAADOGIkXctbDmw9qfl25FtaXj/mcqxbXq7QopgepbggAAIAcQuKFnHSsb1DPbe/UtUunK1GTZWzKigt05eJ6PbjpgNxT7jAAAAAAZB2JF3LSU2+0q38ofkbTDJPesaxBrV292tTanYHIAAAAgDNH4oWc9NCWg6qaUqSL5lWf8bnXnDtdMRPTDQEAAJAzSLyQcwaH4npsa5vefs50FZ5FZcLaqSVaPa9GD5F4AQAAIEeQeCHnvLzniA73DOjaZWc+zTDpuuUN2rK/W3sP9aQxMgAAAODskHgh5zy85aCKC2K6akn9WV8juTbs8dfb0hUWAAAAcNZIvJBznt/RqQvnVWlqydnv7z2vdorqphbrlb1H0hgZAAAAcHZIvJBTegeGtGV/t1bNPfOiGsOZmVbOqdJ6Ei8AAADkABIv5JRNrd0aGHKtnFM17mutnFOlHe3H1dUzkIbIAAAAgLNH4oWc8mowQpWexCsxavZqC6NeAAAAiBaJF3LK+r1HNLOyVA0VpeO+1vlzKmUmphsCAAAgciReyCnr9x5Jy2iXJFWUFmlR/VQSLwAAAESOxAs549Dxfu051KML0pR4SYkpi6/sOSx3T9s1AQAAgDNF4oWckc71XUkr51bpcM+A9rCRMgAAACJE4oWc8creI4qZdF5jZdqumUzimG4IAACAKJF4IWe8uveIljRMU/k4Nk4e7ZyGaSorKtAre0i8AAAAEJ2sJl5mNsfMHjOzLWa2ycw+GbTXmNlDZrYt+FodtJuZfdXMms3sNTO7cNi1bg36bzOzW7P5cyD93F2vtqSvsEZSYUFM5zVWMuIFAACASGV7xGtQ0qfdfamkSyXdbmbLJH1G0iPuvljSI8GxJN0gaXHwuE3SN6REoibpc5IukbRG0ueSyRry067OHh3pGUh74iUl1nltbu1W3+BQ2q8NAAAAjEVWEy933+/uLwffH5W0RVKjpJsl3R10u1vSu4Pvb5b0HU94XlKVmc2U9E5JD7n7IXc/LOkhSddn8UdBmiULa6SzomHSyjlV6h+Ka8v+o2m/NgAAADAWka3xMrMmSaskvSCpwd33S4nkTNL0oFujpL3DTmsJ2k7Wnuo+t5nZOjNb197ens4fAWm0fu8RTSku0JKGaWm/9psFNvYcTvu1AQAAgLGIJPEys6mSfijpD9y9+1RdU7T5KdrDje53uPtqd19dX19/5sEiK17Ze0TnNVaqIJbqn3Z8ZlaWqqGihHVeAAAAiEzWEy8zK1Ii6fquu/8oaD4YTCFU8LUtaG+RNGfY6bMltZ6iHXmob3BIW1q7M7K+S5LMTCvnVJF4AQAAIDLZrmpoku6UtMXdvzTsqfskJSsT3irpJ8PaPxxUN7xUUlcwFfEBSdeZWXVQVOO6oA15aMv+o+ofimcs8ZKklXOqtauzR4eP92fsHgAAAMDJpG/DpLG5XNKHJG0ws/VB259I+oKke83so5L2SHp/8Nz9km6U1CypR9JHJMndD5nZX0laG/T7S3c/lJ0fAemWLKyxcm4mE69gnVfLEb3tnOmn6Q0AAACkV1YTL3d/WqnXZ0nSNSn6u6TbT3KtuyTdlb7oEJX1e49o+rQSzagozdg9zp9dqZhJ6/eQeAEAACD7IqtqCCSt35vYODkxEzUzyksKtaRhGuu8AAAAEAkSL0Sqq2dAOzuOZ2T/rtHOa6zUptYuJQZSAQAAgOwh8UKkNrV2SUokRZm2fFaFOo71q+1oX8bvBQAAAAxH4oVIbWpNbOO2fFZFxu+1PEjuNu7ryvi9AAAAgOFIvBCpTa1dmlFRqtqpJRm/19KZFTL7ZbIHAAAAZAuJFyK1eX93Vka7JGlqSaGaasvfnN4IAAAAZAuJFyLTOzCk7e3Hs5Z4SdKyWRWMeAEAACDrSLwQma0Hjmoo7lo2K/OFNZKWz6pQy+ET6uoZyNo9AQAAABIvRCY55S+bI17LgyRv036mGwIAACB7SLwQmU2t3aosK9Ls6rKs3TOZ5G1muiEAAACyiMQLkdnU2q1lMytkZlm7Z93UEjVUlLDOCwAAAFlF4oVIDA7FtTWLFQ2HWz6rksqGAAAAyCoSL0RiR8dx9Q3GtbwxisSrQtvbj6t3YCjr9wYAAMDkROKFSPyysEb2KhomLZ9VoaG4a+uBo1m/NwAAACYnEi9EYtO+bpUUxrSgrjzr936zsiHTDQEAAJAlJF6IxKbWbp07s0KFBdn/FZxdXaaK0kIKbAAAACBrzvpdr5ldamZrzeyYmfWb2ZCZ8U4Wp+Xu2tTaFUlhDUkyMy2bVUHiBQAAgKwZz3DDP0n6gKRtksokfUzSP57uJDO7y8zazGzjsLY/N7N9ZrY+eNw47LnPmlmzmb1uZu8c1n590NZsZp8Zx8+BLGs5fELdvYORJV5SYrrh1v3dGhyKRxYDAAAAJo9xzfNy92ZJBe4+5O7/KultYzjt25KuT9H+ZXdfGTzulyQzWybpFknLg3O+bmYFZlYg6WuSbpC0TNIHgr7IA8mRpigKayQtn1WhvsG4dnQcjywGAAAATB6F4zi3x8yKJa03s7+VtF/SaSsluPuTZtY0xnvcLOked++TtNPMmiWtCZ5rdvcdkmRm9wR9N5/Zj4AobG7tUkHMdO6MaZHFMLzAxpKG6OIAAADA5DCeEa8PSSqQ9AlJxyXNkfSr47jeJ8zstWAqYnXQ1ihp77A+LUHbydqRBza1dmthfblKiwoii2FhfblKCmPatI91XgAAAMi8s0683H23u59w9253/wt3/1Qw9fBsfEPSQkkrlRg5+2LQbqlufYr2EDO7zczWmdm69vb2swwP6bSptTvSaYaSVFgQ07kzplFgAwAAAFkxnqqGN5nZK2Z2yMy6zezo2VY1dPeDwTqxuKRv6pfTCVuUGElLmi2p9RTtqa59h7uvdvfV9fX1ZxMe0qjzWJ8OdPdGWlgjadmsSm1q7ZJ7ypwdAAAASJvxTDX8B0m3Sqp19wp3n+buZ/Vu2sxmDjt8j6RkxcP7JN1iZiVmNl/SYkkvSlorabGZzQ/Wmd0S9EWOS44wLcuBxGtFY4W6ewfVcvhE1KEAAABgghtPcY29kjb6GQ4XmNn3JF0tqc7MWiR9TtLVZrZSiemCuyR9XJLcfZOZ3atE0YxBSbe7+1BwnU9IekCJdWZ3ufumcfwsyJI3KxrOjHaqoSStCKY7btjXpTk1UyKOBgAAABPZeBKvP5J0v5k9Iakv2ejuXzrVSe7+gRTNd56i/+clfT5F+/2S7h9ztMgJG/d1aW7NFFVOKYo6FJ0zY5oKY6YN+7p043kzT38CAAAAcJbGk3h9XtIxSaWSitMTDia61/Yd0fmNVVGHIUkqLSrQkoZp2rivK+pQAAAAMMGNJ/Gqcffr0hYJJrwjPf3ae+iEPnjJvKhDedP5syv1i00H5O4yS1UsEwAAABi/8RTXeNjMSLwwZhuCkaXzGqNf35W0orFSR3oGKLABAACAjBpP4nW7pF+Y2YnxlpPH5JBMvFZEvIfXcMkkcAPTDQEAAJBB49lAeZq7x9y9bLzl5DE5bGjp0rza3CiskTS8wAYAAACQKeNZ4yUzO19S0/DruPuPxhkTJqgN+7p0wZzcKKyRVFpUoHNmUGADAAAAmXXWiZeZ3SXpfEmbJMWDZpdE4oWQw8f71XL4hD50ae4U1kg6r5ECGwAAAMis8Yx4Xeruy9IWCSa0XCyskbSisVL3rN2rlsMn2EgZAAAAGTGe4hrPmRmJF8YkmXgtz8HEiwIbAAAAyLTxJF53K5F8vW5mr5nZBjN7LV2BYWLZ0NKlptopqizLncIaSRTYAAAAQKaNZ6rhXZI+JGmDfrnGC0hpw74uXTivOuowUqLABgAAADJtPCNee9z9Pnff6e67k4+0RYYJ49Dxfu07ckLnNebubgPnNVZqw74uuXvUoQAAAGACGk/itdXM/sPMPmBm700+0hYZJow3N07OwfVdSSsaK3WkZ0Ath09EHQoAAAAmoPFMNSyT1CfpumFtlJNHyIaWI5JyO/EaXmCDyoYAAABIt7NOvNz9I+kMBBPXhn1dml9XrorS3CuskXTuzGkqKkgU2LjxvJlRhwMAAIAJZjwbKJdK+qik5ZJKk+3u/ptpiAsTyIaWLq1uqok6jFMqKSzQkgYKbAAAACAzxrPG698kzZD0TklPSJot6Wg6gsLE0XGsT61dvTm5cfJoFNgAAABApown8Vrk7n8m6bi73y3pXZLOS09YmCiShTXOm537iRcFNgAAAJAp40m8BoKvR8xshaRKSU2nO8nM7jKzNjPbOKytxsweMrNtwdfqoN3M7Ktm1hxs0nzhsHNuDfpvM7Nbx/FzIIM2tiQSr+WzcreUfNLwAhsAAABAOo0n8bojSJD+H0n3Sdos6W/GcN63JV0/qu0zkh5x98WSHgmOJekGSYuDx22SviElEjVJn5N0iaQ1kj6XTNaQW17b16UF9eWalsOFNZKSBTZeayHxAgAAQHqNJ/GqlPQRSaslfU2JpGvQzFae6iR3f1LSoVHNN0u6O/j+bknvHtb+HU94XlKVmc1UYl3ZQ+5+yN0PS3pI4WQOEXN3vbr3SF6s75ISBTaWz6rUS7tH/3oCAAAA4zOexOsiSb8tqVHSLEm/JelqSd80sz86w2s1uPt+SQq+Tg/aGyXtHdavJWg7WXuImd1mZuvMbF17e/sZhoXx2HOoR21H+3K+ouFwa+bX6NW9XeodGIo6FAAAAEwg40m8aiVd6O6fdvdPKzHyVS/pKkm/kYbYJMlStPkp2sON7ne4+2p3X11fX5+msDAWL+5MjBytyaPEa/W8avUPxZluCAAAgLQaT+I1V1L/sOMBSfPc/YSkvjO81sFgCqGCr21Be4ukOcP6zZbUeop25JC1uw6psqxIi6dPjTqUMbs4SBLX7mK6IQAAANJnPInXf0h63sw+Z2afk/SMpO+ZWbkShTbOxH2SkpUJb5X0k2HtHw6qG14qqSuYiviApOvMrDooqnFd0IYcsnbXYV3cVK1YLNUAZW6qLi/W4ulTSbwAAACQVoVne6K7/5WZ3S/pCiWm/v22u68Lnv7gyc4zs+8psRaszsxalKhO+AVJ95rZRyXtkfT+oPv9km6U1CypR4liHnL3Q2b2V5LWBv3+0t15p5xD2o72amfHcX1gzZzTd84xq5tq9NNXWzUUdxXkUdIIAACA3HXWiZckuftLkl46w3M+cJKnrknR1yXdfpLr3CXprjO5N7Jn3a7Dkn45dS+frJlfre+9uEevHziqZXmw/xgAAABy33imGgIn9eLOQyorKtCKPCklPxzrvAAAAJBuJF7IiBd3HtKquVUqKsi/X7HGqjLNrCzViyReAAAASJP8e1eMnNfdO6AtB7rzcpqhJJmZLm6q0bpdh5SY7QoAAACMD4kX0u6l3YflntiMOF9dPL9GB7v7tPfQiahDAQAAwARA4oW0W7vzkApjplVzq6IO5axd3FQtSUw3BAAAQFqQeCHt1u46pBWNlZpSPK6imZFaMn2aKsuKtHYniRcAAADGj8QLadU7MKRX93bl9TRDSYrFTKvnVWvtbhIvAAAAjB+JF9Lq1b1H1D8Uz9vCGsOtbqrRjvbj6jjWF3UoAAAAyHMkXkir5N5Xq+dVRxzJ+K2Zn/gZ1rHOCwAAAONE4oW0enHXYS1pmKrq8uKoQxm3FY2VKimMae2uw1GHAgAAgDxH4oW0GYq7Xt59OO/XdyWVFBbogjlVb47iAQAAAGeLxAtps7m1W8f6BifE+q6kSxfUauO+LnWyzgsAAADjQOKFtHloy0HFTLp8UV3UoaTNdcsaFHfpka1tUYcCAACAPEbihbR5cNMBrZ5Xo7qpJVGHkjbLZ1VoVmWpHtx0MOpQAAAAkMdIvJAWuzuPa+uBo3rnihlRh5JWZqbrls/QU9va1dM/GHU4AAAAyFMkXkiLBzYdkJSYmjfRXLesQX2DcT35RkfUoQAAACBPkXghLX6x8YCWz6rQnJopUYeSdhfPr1FlWZEe3Hwg6lAAAACQp0i8MG5t3b16ec8RvXP5xJpmmFRUENM1507XI1vaNDgUjzocAAAA5KGcSrzMbJeZbTCz9Wa2LmirMbOHzGxb8LU6aDcz+6qZNZvZa2Z2YbTRT14Pbk4UnpioiZckXbd8hrpODOhF9vQCAADAWcipxCvwNndf6e6rg+PPSHrE3RdLeiQ4lqQbJC0OHrdJ+kbWI4WkxPquptopWtIwNepQMuaqJXUqKYxR3RAAAABnJRcTr9FulnR38P3dkt49rP07nvC8pCozmxlFgJNZ14kBPbe9U+9cPkNmFnU4GTOluFBXLq7XQ5sPyt2jDgcAAAB5JtcSL5f0oJm9ZGa3BW0N7r5fkoKv04P2Rkl7h53bErSNYGa3mdk6M1vX3t6ewdAnp8e2tmkw7rpuAk8zTLpueYP2HTmhTa3dUYcCAACAPJNridfl7n6hEtMIbzezq07RN9XwSmgowt3vcPfV7r66vr4+XXEi8MCmA5o+rUSr5lRFHUrGXXPudMUssVE0AAAAcCZyKvFy99bga5uk/5K0RtLB5BTC4Gtb0L1F0pxhp8+W1Jq9aNE7MKTHX2/XO5Y1KBabuNMMk2qnlmh1U82bxUQAAACAscqZxMvMys1sWvJ7SddJ2ijpPkm3Bt1ulfST4Pv7JH04qG54qaSu5JREZMeTb7TrxMDQhK5mONp1yxq09cBR7e48HnUoAAAAyCM5k3hJapD0tJm9KulFST9z919I+oKkd5jZNknvCI4l6X5JOyQ1S/qmpN/NfsiT209f269ppYW6dEFt1KFkzfUrZshMunfd3tN3BgAAAAKFUQeQ5O47JF2Qor1T0jUp2l3S7VkIDSm0HO7Rzzbs162XNam4MJfy98yaXT1F71w2Q//+/B797tWLVF6SM/8JAQAAIIdNnnfMSKs7n94pk/SxK+dHHUrW3fbWBeo6McCoFwAAAMaMxAtn7PDxft3z4l79yspZmlVVFnU4WXfh3Gpd3FStO5/eqcGheNThAAAAIA+QeOGM/dvzu3ViYEgfv2ph1KFE5reuXKCWwyd0/0ZKywMAAOD0SLxwRnoHhvTtZ3fp7edO1zkzpkUdTmSuXdqgBfXluuPJ7UosNwQAAABOjsQLZ+Q/1+3VoeP9+vhVC6IOJVKxmOm3rlygjfu69dyOzqjDAQAAQI4j8cKYDQ7F9c2ndmrV3CqtmV8TdTiRe8+qRtVNLdYdT+6IOhQAAADkOBIvjNnPNx7QnkM9+vhVC2VmUYcTudKiAt16WZMef71drx84GnU4AAAAyGEkXhiTeNz1z09s14L6cl23rCHqcHLGr186T2VFBfrHR7dFHQoAAAByGIkXxuRfntyhTa3d+sTbFikWY7Qrqbq8WB9/6wL99LX9+sn6fVGHAwAAgBxF4oXTWrfrkP7+wdf1rvNn6j2rGqMOJ+d84m2LtHpetf70vzZqd+fxqMMBAABADiLxwikdPt6v3//eK2qsKtP/fe95rO1KobAgpn+4ZaViJv3+915R/yCbKgMAAGAkEi+clLvrD3/wqtqP9emffm2VKkqLog4pZ82unqK/fd/5erWlS1986PWowwEAAECOIfHCSd31zC49vKVNn71hqc6fXRV1ODnv+hUz9cFL5upfntihJ99ojzocAAAA5BASL6T09LYOfeHnW3Tt0gZ95PKmqMPJG3920zItaZiqT927Xptau6IOBwAAADmCxAsjxOOurz/erA/f9YLm1Zbr799/Puu6zkBpUYG+/sELVVQQ03u//qx++FJL1CEBAAAgB5B44U1dPQO67d/W6W9/8bpuOG+mfnz75aqaUhx1WHln0fRp+u/fu0Kr5lbp0//5qv7sxxspuAEAADDJFUYdAHLDy3sO65P3vKIDXb368/+xTLe+pYmRrnGom1qif//oJfrbB17XHU/u0KbWLv3d+y/QwvqpUYcGAACACJB4TWIn+of036+16j9e2KP1e49oZmWpvv/xy3Th3OqoQ5sQCgti+pMbl+r82ZX6ox+8pmu++ISuWFSnX790nq5dOl2FBQw4AwAATBbm7lHHcNbM7HpJX5FUIOlb7v6FU/VfvXq1r1u3Liux5aJ43LWr87g2tXbrxZ2H9OP1+3S0d1AL6sv1a2vm6v2r56iyjJLxmdB2tFf3rt2r/3hhj1q7ejWjolQ3r5qlVXOqtKKxUo1VZYwwjmJmL7n76qjjAAAASIe8TbzMrEDSG5LeIalF0lpJH3D3zSc7ZyyJ1+BQXL0RrsdJ/nv4m8eJA5fLPdE+FHfF3TUUTzwG467egSGdGBhS78CQ+gbiOtzTr85j/eo43qfOY/3a09mjTa1dOt4/JEkqLozp+uUz9GuXzNUl82t4058lg0NxPbq1Tf/2/G49u71TQ/HEv3RNebGWz6rQzMpS1U4tUd3UEtVNLVZFWZHKigpUWlSg0qKYSgoLVBgzFQSPmJliJsXMZCaZTDJp9D/n6H/dKP+9SwpjKhrDaB+JFwAAmEjyearhGknN7r5DkszsHkk3Szpp4jUWT25r129+e+KMihUXxFQ3tVgzq8r0qxfN1opZlVo2q0JLGqapuJCpbtlWWBDTdctn6LrlM9Q7MKQt+7u1cV+XNuzr0pb9R/XGwaPqPNavwXh+fiAyFn/3vvP1/tVzog4DAAAgq/I58WqUtHfYcYukS0Z3MrPbJN0mSXPnzj3tRRfVT9Of3rg0TSGendBohZksaDcpMdIRMxUGIx6FBabSwsSoSElRTKVFBaqeUqzaqcWaVlLIaFaOKi0q0Kq51Vo1ak1dPO7qOjGgjmN9Oto3qN5gJLN3IK7egaE3RzwH46543BX3xEhpPBgRPd0odtSD3GzGDQAAJqN8TrxSZROht5TufoekO6TEVMPTXXRu7RT91lULxh8dcJZiMVN1ebGqyynlDwAAMFHk81yzFknD5yvNltQaUSwAAAAAcFL5nHitlbTYzOabWbGkWyTdF3FMAAAAABCSt1MN3X3QzD4h6QElysnf5e6bIg4LAAAAAELyNvGSJHe/X9L9UccBAAAAAKeSz1MNAQAAACAv5O0GymfDzNol7R5D1zpJHRkOJ52IN7OIN7NOFu88d6/PdjAAAACZMKkSr7Eys3XuvjrqOMaKeDOLeDMr3+IFAAA4G0w1BAAAAIAMI/ECAAAAgAwj8UrtjqgDOEPEm1nEm1n5Fi8AAMAZY40XAAAAAGQYI14AAAAAkGGTPvEys7vMrM3MNg5r+3Mz22dm64PHjVHGOJyZzTGzx8xsi5ltMrNPBu01ZvaQmW0LvlZHHat0ynhz8jU2s1Ize9HMXg3i/Yugfb6ZvRC8vt83s+KoY5VOGe+3zWznsNd3ZdSxJplZgZm9Ylj+rKwAACAASURBVGY/DY5z8rUFAABIp0mfeEn6tqTrU7R/2d1XBo/7sxzTqQxK+rS7L5V0qaTbzWyZpM9IesTdF0t6JDjOBSeLV8rN17hP0tvd/QJJKyVdb2aXSvobJeJdLOmwpI9GGONwJ4tXkv5w2Ou7ProQQz4pacuw41x9bQEAANJm0ide7v6kpENRxzFW7r7f3V8Ovj+qxBvYRkk3S7o76Ha3pHdHE+FIp4g3J3nCseCwKHi4pLdL+kHQnkuv78nizUlmNlvSuyR9Kzg25ehrCwAAkE6TPvE6hU+Y2WvBVMScmLY3mpk1SVol6QVJDe6+X0okO5KmRxdZaqPilXL0NQ6mwq2X1CbpIUnbJR1x98GgS4tyKHkcHa+7J1/fzwev75fNrCTCEIf7B0l/JCkeHNcqh19bAACAdCHxSu0bkhYqMXVrv6QvRhtOmJlNlfRDSX/g7t1Rx3M6KeLN2dfY3YfcfaWk2ZLWSFqaqlt2ozq50fGa2QpJn5V0rqSLJdVI+uMIQ5QkmdlNktrc/aXhzSm65sxrCwAAkC4kXim4+8HgzWxc0jeVePOdM8ysSIkk5rvu/qOg+aCZzQyen6nE6EdOSBVvrr/GkuTuRyQ9rsTatCozKwyemi2pNaq4TmZYvNcHUzzd3fsk/aty4/W9XNKvmNkuSfcoMcXwH5QHry0AAMB4kXilkExgAu+RtPFkfbMtWBNzp6Qt7v6lYU/dJ+nW4PtbJf0k27GlcrJ4c/U1NrN6M6sKvi+TdK0S69Iek/S+oFsuvb6p4t06LAk3JdZMRf76uvtn3X22uzdJukXSo+7+QeXoawsAAJBOk34DZTP7nqSrJdVJOijpc8HxSiWmPO2S9PHk+qmomdkVkp6StEG/XCfzJ0qsm7pX0lxJeyS9390jLxpying/oBx8jc3sfCUKPBQo8cHEve7+l2a2QIlRmhpJr0j69WA0KVKniPdRSfVKTOVbL+m3hxXhiJyZXS3p/7j7Tbn62gIAAKTTpE+8AAAAACDTmGoIAAAAABlG4gUAAAAAGUbiBQAAAAAZRuIFAAAAABlG4gUAAAAAGUbiBQAAAAAZRuIFAAAAABlG4gUAAAAAGUbiBQAAAAAZRuIFAAAAABlG4gUAAAAAGUbiBQAAAAAZRuIFAAAAABlG4gUAAAAAGUbiBQAAAAAZVhh1ANlUbCVeqvKowwAwBr06rn7vs6jjGIu6ujpvamqKOgwAY/TSSy91uHt91HEAmFwmVeJVqnJdYtdEHQaAMXjBH4k6hDFramrSunXrog4DwBiZ2e6oYwAw+TDVEAAAAAAyjMQLAAAAADKMxAsAAAAAMozECwAAAAAybFIV1wAA4Ey1H+3TnkM9Oto7oIubalRewv86AQBnjv97AACQQv9gXF97rFlfe6xZg3GXJFVPKdLHrlyg37pygYoLmTQCABg7Ei8AAEY51jeoD37zeb3a0qV3r5ylm1c1Kmamu5/dpb974HVtau3SV29ZpcICki8AwNiQeAEAMEw87vrU99drY2u3vvZrF+pd589887m3LqnXt57aof/vZ1s0pXiD/u5958ssL/b5BgBEjMQLAIBhvvroNj24+aD+35uWjUi6kj525QId7R3UVx7ZptXzqnXLmrkRRAkAyDfMkQAAILDt4FH946PNes+qRn3k8qaT9vuDaxdrzfwa/d+fb1XHsb7sBQgAyFskXgAABD5//xZNKS7Qn9207JRTCM1Mf/2eFerpH9Rf/2xLFiMEAOQrEi8AE56Z/W8z22RmG83se2ZWambzzewFM9tmZt83s+Kgb0lw3Bw83xRt9MiWJ95o1+Ovt+v3375YNeXFp+2/aPo0ffyqhfrRK/v0WsuRLEQIAMhneZt4mdkcM3vMzLYEb6g+GXVMAHKPmTVK+n1Jq919haQCSbdI+htJX3b3xZIOS/pocMpHJR1290WSvhz0wwTn7vrCz7dqbs0Uffgt88Z83sffukAVpYX6+mPbMxgdAGAiyNvES9KgpE+7+1JJl0q63cyWRRwTgNxUKKnMzAolTZG0X9LbJf0geP5uSe8Ovr85OFbw/DVG2boJ79ntndqyv1ufeNsilRQWjPm8aaVF+vBlTXpg8wE1tx3LYIQAgHyXt4mXu+9395eD749K2iKpMdqoAOQad98n6e8l7VEi4eqS9JKkI+4+GHRr0S//fjRK2hucOxj0r81mzMi+O5/eqbqpxfqVlbPO+NyPXN6kksKY/vkJRr0AACeXt4nXcMEajFWSXkjx3G1mts7M1g2IylPAZGNm1UqMYs2XNEtSuaQbUnT15CmneG74dd/829Le3p6ucBGB7e3H9OjWNv36pfNUWjT20a6k2qkluuXiufrxK/vUdrQ3AxECACaCvE+8zGyqpB9K+gN37x79vLvf4e6r3X11kUqyHyCAqF0raae7t7v7gKQfSXqLpKpg6qEkzZbUGnzfImmOJAXPV0o6NPqiw/+21NfXZ/pnQAZ9+5ldKi6M6dcvHfvartE+dNk8DcZdP3p5XxojAwBMJHmdeJlZkRJJ13fd/UdRxwMgJ+2RdKmZTQnWal0jabOkxyS9L+hzq6SfBN/fFxwreP5Rdw+NeGFi6B0Y0o9f2aebzpupuqln/+HcwvqpuripWveu3St+XQAAqeRt4hW8gbpT0hZ3/1LU8QDITe7+ghJFMl6WtEGJv3t3SPpjSZ8ys2Yl1nDdGZxyp6TaoP1Tkj6T9aCRNQ9uPqijfYN630Wzx32t/7l6jnZ0HNfaXYfTEBkAYKLJ28RL0uWSPiTp7Wa2PnjcGHVQAHKPu3/O3c919xXu/iF373P3He6+xt0Xufv73b0v6NsbHC8Knt8RdfzInB+93KJZlaW6dMH466e86/yZmlpSqO+v3ZuGyAAAE03eJl7u/rS7m7uf7+4rg8f9UccFAMgPbd29evKNdr3nwkbFYuPfMWBKcaH+xwUz9bMNrTreN3j6EwAAk0reJl4AAIzHj9fvU9yl9144/mmGSTevbFTvQFyPbm1L2zUBABMDiRcAYFL671f364LZlVpYPzVt17y4qUb100r0s9f2p+2aAICJgcRrMjA7uwcATFB7D/Vow74u3XDezLRetyBmunHFDD32epuOMd0QADAMiRdOjmQMwAT1wKYDkqQbVsxI+7Xfdf4s9Q3G9ciWg2m/NgAgf5F4TUYWG9sj5bkkYgDy3883HtDSmRWaV1ue9muvnlet6Uw3BACMQuIFAJhUDnb36qXdhzMy2iVJsZjpxvNm6ok32nWifygj9wAA5B8Sr4kmndMDxzoKBgB5JJPTDJOuXdqgvsG4nm7uyNg9AAD5hXfSk8Go5MlidtaPcDLGOjAA+eWhzQe1oK5cixumZewea+bXaFpJoR7ezDovAEACiRcSxrju6/SJGL9SAHLX8b5BvbDjkK5ZOj2j9ykujOmqc+r1yNY2xeOe0XsBAPID75IBAJPG080d6h+K623nZjbxkqRrl05Xx7E+vbavK+P3AgDkPhKvySjVKFXMQg8riIUeY5m2yPRDALnqsa1tmlZSqIubajJ+r6uXTFfMRFl5AIAkEi8ErKAg9FAsFnqcLhE7o9L0AJBF7q5Ht7bpqiX1KirI/N+l6vJirZ5Xo4e3tGX8XgCA3Mc74skoxehWqlEqKywMPVRUNOKRalSMUTAAuWhTa7fajvZlZZph0tXn1mvL/m61dfdm7Z4AgNxE4gUAmBQe29omM+nqc+qzds+3Lknc68ltlJUHgMmuMOoAkGaeonrWWAaXYily8IKC014qZa2uVDGkqurlbCwKIHue3Nau8xorVTe1JGv3XDqjQnVTS/TEG+1630Wzs3ZfAEDuIfGCJMlSTf1LtQZi1HqtlDldPJ6iLZx4ucdHN5w8QAAYh6O9A3p5zxH99lsXZPW+sZjpqiV1emxrm4biroIY06wBYLJiqiEAYMJ7dnunhuKuKxdnb5ph0luX1Otwz4A2UFYeACY1RryQkOpT2FRTDQtH/cqkGhUbPZIlSUPhaYU+FBvdcKoIAeCsPbWtXeXFBbpwbnXW733FojqZSU++0a6Vc6qyfn8AQG4g8cJJWap1X0WjfmUKw8lZyimDQ+FkzEa1haYenuxaAHCGntrWocsW1qq4MPsTPWqnlui8xko98Ua7fv+axVm/PwAgN5B4TQapEprRUu21lWLEy0cnXqkulapxcDDcNjAw8tqjR8ASjae9HwCcyu7O49rd2aPfvHx+ZDFcsahO//LkDh3tHdC00qLI4gAARIc1XgCACS1Zyv3KxXWRxXDFojoNxV0v7jwUWQwAgGgx4oWTS7V+a9SIl6eatpNiiqINpljj1T8wqk94VMw9xfgZ0w8BnIFnmzs0q7JU8+vKI4vhwnnVKimM6enmDl2ztCGyOAAA0SHxmoxS7qk1humICidaXhKeMuMppijGUhTXsL6+kQ2jph4mwmL/LwBnbyjuem5Hp65d2pB624wsKS0q0Jr5NXqmmY2UAWCyIvFCQqoEJ8WblNFJ1VBpil+hFCNSFi8NtcV6+0ee1tcf6pNyFGxsOSLwJjOrkvQtSSuU2Pf7NyW9Lun7kpok7ZL0P939sCXenX9F0o2SeiT9hru/HEHYSIPNrd060jOgyxfVRh2KLl9Upy/8fKvajvZq+rTw30QAwMTGGi8Ak8FXJP3C3c+VdIGkLZI+I+kRd18s6ZHgWJJukLQ4eNwm6RvZDxfp8sz2xAjTWxZGt74r6YpFiRiebe6MOBIAQBTydsTLzO6SdJOkNndfEXU8OW30CNRYh4xSraUqGDkKNlQSzt1TrfuyVANqJ0Z+4mu9veFr9YdHwUJ7grHmC6dgZhWSrpL0G5Lk7v2S+s3sZklXB93ulvS4pD+WdLOk77i7S3rezKrMbKa7789y6EiDZ5o7tHj6VDVURD/CtGxmhaqnFOmpbR1696rGqMMBAGRZ3iZekr4t6Z8kfSfiOPJOqnVTnmoN1himH8aLw0nWYFmKZMzCa8FifWUjjgt6womXnUiRjA2Mmn7Imi+c2gJJ7ZL+1cwukPSSpE9KakgmU+6+38ymB/0bJe0ddn5L0DYi8TKz25QYEdPcuXMz+gPg7PQNDmntrkO65eLc+PeJxUyXLazV8zs65e6RrjkDAGRf3iZe7v6kmTVFHceEkWrUaAwjSfGi8BuHwbJwWzzFRssFfSUjjmPd4U+k7XhxuG1UNUQ2XsZpFEq6UNLvufsLZvYV/XJaYSqp3g2HfqHc/Q5Jd0jS6tWr+YXLQa/sOaLegbjesjD69V1Jly2s0/0bDmjPoR7Nq42uyiIAIPsm/BovM7vNzNaZ2boB9Z3+BAATTYukFnd/ITj+gRKJ2EEzmylJwde2Yf3nDDt/tqTWLMWKNHp2e6diJl2yIIcSryCWZ7ezzgsAJpu8HfEaq+GfSldYDZ9KS6nXeMVTtKWYfjh6JMkLxjbipbJwU2xw5K9fUXe4U8HRklCbnTgxMoZw4UPgTe5+wMz2mtk57v66pGskbQ4et0r6QvD1J8Ep90n6hJndI+kSSV2s78pPz2/v1IrGSlWWhac6R2VhfbmmTyvRc9s79YE1uTEFEgCQHRM+8cLYeIqpeZ5i02MbPH1hjqGScOI1FM6fpPjIAdfirnCn0iPhZMyO94yKKUXJ+RRtmNR+T9J3zaxY0g5JH1FixP9eM/uopD2S3h/0vV+JUvLNSpST/0j2w8V4negf0it7D+s3L58fdSgjmJnesrBWTzezzgsAJhsSLwATnruvl7Q6xVPXpOjrkm7PeFDIqHW7D2lgyHVZDq3vSrpsYa1+vL5VzW3HtLhhWtThAACyJG8TLzP7nhKloOvMrEXS59z9zmijymOpKhjGU0w1HD3ileq0FLN6BqaG2zw28pPevqPhE4s7p4TaYl3HRjb0pVi7N4ZpkgAmrue2d6owZrq4qSbqUEKSe4o9u72TxAsAJpG8Tbzc/QNRx5C3UlYwHNsaLxvVlqrkfDxcwFCD5eF+o6cfnjgWnnJT1h6eflhyaNT+X8ePh/q4pagbQ9l5YNJ4bkenLphTpfKS3Pvf3JyaKWqsKtNz2zt161uaog4HAJAlufd/JEQi9d5e4WTMRu2hFRsY2yjSUFn4Wl448tzenvCvY29tilGwtpGjYNZ9LNTHUmy8PNZ9owHkt2N9g3qtpUu/89aFUYdyUpctrNXDWw4qHnfFYqzzAoDJYMKXkwcATC5rdx7SUDw313clXbagVkd6BrT1wNGoQwEAZAkjXkgY61TDUSNeBX3h82KpllcVhUfGCipGjkr19oc/BzjRHm6bUjlyqmFRZ3iTZfWkmO/Iui9gUnhuR6eKC2K6aF511KGcVDIpfG5Hp5bNqog4GgBANpB44aQ8VaLSPzDiMNYX7hMbCDVJFk5wppX3jjhO9bnvifpwcY2+mpHrvooOpCg5fzQ8/dAHUpSYZ90XMOE8t71TK+dWqbQoxQcwOWJWVZnm1U7Rc9s79dErcqvkPQAgM0i8kJBq5CfVuq+BkVlVQW84mSnoS1W8I7yGobJsZOJVXhJel9XaEC6ucaJ25JupKRXh5MwOh0fBLEX1Q9Z9ARNL14kBbWrt0u+9fXHUoZzWZQtq9bMN+zUUdxWwzgsAJjzWeAEAJowXdx5S3JXT67uSLltYq6O9g9rc2h11KACALGDECyeXajhocOQIl6UY8So8ER7xsoHwp7lTikaOcNWVhqcHdjSUh9pO1I/c92agqjTUp+RgeKRMPT3httHTKVnzBeS157Z3qqQwplVzq6IO5bQuW5Bc59Wh82ZXRhwNACDTSLxwUilLzI9aJ2V94emBRT3h82K9px9cXTilI9TWUR/eeXnb9JFtfTUpSs6Xp1j31R1eRRZa98WaLyCvPb+jUxfNq1ZJYe6u70qaXlGqBfXlem57p267KndL3wMA0oPECyc3hkqH1htOvAqPh5OXwhRVBvuGRv76NRR1hfqsrGoJtW2dMWPE8Yna8OhW+bTwuq9YiuqHo9d9seYLyF9Hevq15UC3/ve1S6IOZcwuW1Crn6xv1eBQXIUFzP4HgImMv/IAgAnh+R2H5HmyvivpsoW1OtY3qA37wh88AQAmFka8cHIp1jv50MghIe9PNeIVXvdVdCw82tTdO3Jt1pRYuOrgyvLdobanG0ZOyemsnxnqM1gVHgUrKR3Dui/2+gLy1nPbO1RWVKALZuf++q6kSxf8cj+vVXNzd98xAMD4kXjhzIyei9cf3rSr4Hi4rbg7nLx0HRuZeMVTDMAuKWoLtV1Yt3fE8X9Pbwj16asOr/sqKQsX4bCCkVMgSbGA/PXcjk6tbqpWcWH+TOaom1qiJQ1T9dz2Tv3u1YuiDgcAkEFZS7zM7F2Slkt6892vu/9ltu4PAJi4Oo716Y2Dx/SeVbOjDuWMvWVhnb6/dq/6B+N5lTQCAM5MVhIvM/tnSVMkvU3StyS9T9KL2bg30mt0pUMfDE8rjB3vDbUVd08LtQ10j5z61z4Y7lM3JTx6dvHUnSOO729YHupzoiZcXGPqtBSVDg+NGhnrDU93pNIhkPue39EpKb/WdyVduqBW3352l15tOaKLm2qiDgcAkCHZGvF6i7ufb2avuftfmNkXJf0oS/dGOo2eajgQToxSlZgv6U5R6fDwyF+/3SfqQn2KK8P7f60sGVnpcN70Q6E+rfXh/b8GK8JTDYtKRiZ/VnAi1MfjJF5Arnt2e6emlhRqxayKqEM5Y5cuqJGZ9GxzJ4kXAExg2Uq8ku9me8xslqROSfOzdG+k06hCE6OLbUiSpxg1KupKse7ryMjRpuaj4cTr6PTwqqt5hSOTsZXV4ZLzO2sbQ239lSn2+yodnXiFp/n4YDj5o+AGkFue396pS+bX5GVJ9qopxVo2s0LP7ejQJ7U46nAAABmSrf9D/dTMqiT9naSXJe2SdE+W7g0AmMAOdPVqR8fxvJxmmHT5ojq9vPuIegcYYQeAiSpbI15/6+59kn5oZj9VosBGeCEQ8o6nKr/eFx7xKjwWnn5Yemjkmqu9h8MloHcNVobaFhaNvOdF5btCfX44fVWorbc6xUbLU0ZNPywKj4qlqtzIui8gdzy7vUNSfq7vSrpsYa3ueHKH1u06rCsWh0f/AQD5L1uJ13OSLpSkIAHrM7OXk23IY6PXfEnygXDBDUtRcKP08Mi1GIfaw+uyNveGK5RdXTqyuMbi4oOhPnW1R0NtvTXhNV7xqSPbYsXhxMtOhKcapvixAUTkmeZO1ZQXa+mM/FvflbSmqUaFMdOz2ztIvABggspo4mVmMyQ1Siozs1WSku9gK5SocogJKGWlw55w4lVyZGS/4o5w0vPy0bmhts6KzSOOGwrC660WVHWG2l6trg+1DVSM3Ni5tDi80bNG7fUliY2WgRzh7np2e4cuW1CrWCzFesw8UV5SqJVzqvTs9vDfLgDAxJDpEa93SvoNSbMlfWlYe7ekP8nwvQEAE9zOjuPa39WrtyzK32mGSW9ZWKt/eqxZ3b0DqihNMe0ZAJDXMpp4ufvdku42s1919x9m8l6ISKpRnhSjQZ5i3VfRkZFtpR3h0aaNHTNDbVunj5ySuKzoeKjPovL2UNuL1UtCbf0VI/8TKC0LT0e0FCNejG3lHzMrkLRO0j53v8nM5itR5KdGiaI/H3L3fjMrkfQdSRcpUYH1f7n7rojCxmk8E4wQXb4w/6fnXbawTl99tFkv7jika5c1RB0OACDNsrXG6xkzu1PSLHe/wcyWSbrM3e/M0v2RRaM3WZYkT7G3V8HRkdMPyzrCa7wOtoXXbLw8v2nE8ZyCjaE+s4vDe3upKhxDX8XIghteFi7AYYUp/jOxFAVBKbiR6z4paYsSU50l6W8kfdnd7wk2ef+opG8EXw+7+yIzuyXo97+iCBin92xzhxqryjSvNv9nr184r0qlRTE93dxB4gUAE1C2Eq9/DR5/Ghy/Ien7kki8JqJUlSdSbbQ8at1X6aFw4lJ0MDzd5sUjI7eAu7hsR6jPtILwmrLyinBb/6hNlYfKw6NuBakKbqRYS0LBjdxlZrMlvUvS5yV9ysxM0tsl/VrQ5W5Jf65E4nVz8L0k/UDSP5mZubOIL9fE467ndnTq2qUNSvyT5reSwgKtmV+rp5s7og4FAJAB2drHq87d75UUlyR3H5TE8ACAbPkHSX+k4G+QpFpJR4K/RZLUokQhIAVf90pv/q3qCvojx2xs7dKRngFdsSj/pxkmXbGoVs1tx3Sgix1XAGCiydaI13Ezq1WwNMbMLlXizcy4mNn1kr4iqUDSt9z9C+O9JtIgxcCAD6UoOz9qxKv4SHgqYFlbeLRpU/uMkcc14ZLz5bHwmrLqKSdCbZ3TRu4dNlQevl9hqr29xlLpkAGSnGBmN0lqc/eXzOzqZHOKrj6G54Zf9zZJt0nS3Lnh6pvIvKe2JUaGLp9QiVe9pK16urlD77so/LcNAJC/spV4fUrSfZIWmNkzkuolvW88FwwWyn9N0juU+LR6rZnd5+6bT30mojCWjZYLU3zCW9oRXrfRdmDaiOOXZs0L9VlR3hpqm1ocTsYOTBv5fnqgPPyfRElpinVfFNzIJ5dL+hUzu1GJzdsrlBgBqzKzwmBUa7ak5C9Ni6Q5klrMrFBSpaTQokF3v0PSHZK0evVq/vkj8PS2Di2dWaH6aeH/RvPVuTOmqba8WE9vayfxAoAJJluJ12ZJ/yWpR9JRST9WYp3XeKyR1OzuOyTJzO5RYm0GiVcuSrXR8qj9vmLHwyNSZZ3hPcFKRq372ngoXPmwqih8rcJYOIahstGJV3j2rZeOcW+v0QU3KLaRE9z9s5I+K0nBiNf/cfcPmtl/KvEB0D2SbpX0k+CU+4Lj54LnH2V9V+7p6R/Uut2H9JHL55++cx6JxUyXL6rT082dcvcJsXYNAJCQrTVe35F0rqS/lvSPkhZL+rdxXvPNdRiB4Ws03mRmt5nZOjNbN6DwiAeASeuPlSi00azEGq5ksZ87JdUG7Z+S9JmI4sMpvLDzkAaGfEKt70q6YlGdOo716fWDR6MOBQCQRtka8TrH3S8YdvyYmb06zmuOaR3G8OlAFVbDp9Y5xAdGjmb5ifBUw+LDKdZ9tY8c8dp/sCrUZ9uU+lDbYDz8OUO8dOQo2MCUFFMIS1JUNUxRYn50pUOqHOYed39c0uPB9zuUGDkf3adX0vuzGhjO2NPbOlRcGNOa+TVRh5J2VyxOJJNPvtGuc2eEt9QAAOSnbCVer5jZpe7+vCSZ2SWSnhnnNZPrMJKGr9FArklZcGPUVLwUiVdhV4rph+1lI46PHQyv79hVHX4zNqU4XNJesZFxDYX3T9ZQWYqCGylKzIemGqaaIsSMNSAtntrWrjVNNSotSjHtN8/NqirTkoapevKNDt121cKowwEApElGpxqa2QYze03SJZKeNbNdZrZTibUTV43z8mslLTaz+WZWLOkWJdZmAAAmsP1dJ/TGwWO6cvHEm2aY9NYl9Xpx5yH19IfXuQIA8lOmR7xuytSF3X3QzD4h6QElysnf5e6bMnU/ZMCouXiji21IqQtulHaOHLkqbQuPeB2umxZqG6juCccwalBqqCQ8ShUvCX+ibilKzFvByM8xnPdLQEY88Xq7JOnqc6ZHHEnmXLWkXt98aqee39Gpt5/bEHU4AIA0yGji5e67M3z9+yXdn8l7IHtGr/mSwnt9SVLxkZFFUsraw1UHe9vDidHxglTzCEcmWvEU/0UMlaQYGC5K0XF0pcPRUw8lKh0CafD46+2aWVmqJQ1Tow4lYy5uqlFpUUxPvtFB4gUAE0S21ngBYaPWO41lry9JinWPHAUr6wzv9XWiI/yrfTxFWfjQ+FaKXClelKLEfIrEK9XeXgDSa2AormeaO3TTBTMndKn10qICXbagVk+80R51KACANMlWOXkAAMbtpd2HdbRvUG9dMnGnGSa9dUm9dnYc1+7O41GHAvz/7N15eBzndef772msuonsGgAAIABJREFU3AAQJLgBpLiIEiVuWihKjrxI1kavlBMrkSaJ5YwTJRN7Yt+M50aZ5NoeO77XuVnHceJEiWXLSWxZ3iLaVqzIWsaStRGSKW4iRYorwH0DCBDE1mf+qGqy0VUEm8RS3ejf53nq6e63q7oPSPSLfus9dV4RGQaa8ZLCEbfIck+0nHzudV9VMSXnq49EZ5/6xkfPM+QuoBxdkADSFTFn1fNINcwtLw8qMS8yVM9sPUx5yrjx0ilJhzLibl40jc/8YDNPbznEh8fYQtEiIqVIAy8paN4fMxg7PfC6r/ITMSXnj0fTCnsnRX/de/oGDo5SMQUxPO5SrfJoY25xjdhrvFRiXmRIntl6iBVzJzOpOmZJhzHmkikTWNAwgSc18BIRGRM08JLCkc9aX4CfHnjdV1zlw6qY676qJsUsjmznH3jFLdUdN/BK5S6qHDPjJSIXr+X4KbYcOMn/ePeipEMZNe9cNI2Hnt9NR3cfE6v0J1tEpJjpGi8RESkKT205BMCtV5ROlb93LppOT3+a57YdSToUEREZIp0+k8IWd91Xznpf3hUtOV/eHm2rPh5NTUqXD5wFi0srtLjrslIxO+Ze4xVT5dBVYl7koj2x+SDzGyYwv2HslpHPFaRVlvPUloOsWjIj6XBERGQINPCSwhZ3/VNu+mHMwCt1Mpp+WHkiutByf9XAtv7KaHqg9cekQMZlEeYOxsZwqWuR0XbydC8v7jjKb5TYtU4VZSnecVkDT205TDrtpJTCLCJStDTwkqLj6Zz1v/qiF2bZqejAq7w9uoByZfXAWam+CdFZqlRPdOBlcfUwcr8QxcyKqdKhyMV5dtsRevu9pNIMM269Yjo/XL+fdS0nuGbO5KTDERGRi6RrvEREpOD9ZPNB6sZXcM2cuqRDGXU3L5pGecp4fNOBpEMREZEh0IyXFJ+cKSLvjc54eXfM+l8d0ZTEipyS1ObRMvRxa3vFpkDmpBZaTKphbOF4lZgXGVRvf5ontxziliumUZ67bEMJqB1XwVsWTOHxjQe4f9Wi2L5FREQKnwZeUnxyBiWxJedjFl62mGvByjoGDrQ87ktd3DUVMdd95XVcXHENVFxDZDAvvHmUtq5e3rVkZtKhJGbVkhn80fc38sbBDi6fMSnpcERE5CKU3qlDEREpKv++cT8TKst428KpSYeSmNuunI4Z/Hij0g1FRIqVZryk+OVRch6iCy8D2KmBs2BlMQsje0W04Ib15jFLFTu7JSIXoq8/zX9sOsg7r5hOdcxnsVRMm1TNtXMm8/imA3z81oVJhyMiIhdBAy8pfvmUnCc+/TC3FL2VxwyyqmKu+0rHlCLM57oslYIWuSAv7zrG0c4e3q01rFi1ZAZ/8qPX2Xmkk3lTJyQdjoiIXCCdkpcxyfv7o1tfX2Sju3vAZl35bfT2Rbf+9MAtX5aKbiICBKl11RUp3nF5Q9KhJO49y2ZiBj94bV/SoYiIyEXQNzwRESlIff1pHtuwn3cumsb4SiVozKwdx3Vz61nz2j5clU9FRIqOBl5SOvr7I5v39g3cursjGz29kc16+yJb5PVjmFlkE5F4z20/wpGOHlZf1Zh0KAXjfctnsf1QB1sOnEw6FBERuUAaeMnY5B7ZPB2z5aYexgyyvLsnssXul04P2PCYTUTy9ui6fdRUl3OT0gzPePeSGZSlTOmGIiJFSAMvKR1xA6HcGbC468B6eyKb9/ZGNvr6Bm5pj24ikpdTPX08vukA71k2k6qYojelasrEKm68dCprXttHWn2KiEhR0cBLRMY0M5ttZk+b2etmtsnMPh6215vZE2a2LbydHLabmX3RzLab2XozuybZn6A0PbH5IKd6+nn/cqUZ5vrA1bNoOd7Fy7uOJR2KiIhcAA28pHTkk36Ye83XObbI7FZf9Bovd89ri2U2cJOh6AP+m7tfAdwAfNTMrgTuB55094XAk+FjgHcBC8PtPuDLox+yfO/VVmbWVnP9vPqkQyk4qxbPZFJVOY807006FBERuQAaeElpi0s/zNni0g9jB2T96QEb6ZhNRp2773f3V8P7J4HXgUZgNfBQuNtDwJ3h/dXA1z3wIlBnZjNHOeyS1nqii59uO8xd1zaR0tp3EeMqy3jv8pn8+4YDdHRHF4sXEZHCpIGXlLY8CnDEXqsVUyExr03XfSXKzOYCVwMvAdPdfT8EgzNgWrhbI5A9ldAStuW+1n1m1mxmzYcPHx7JsEvOI2uDf/67VsxOOJLC9cFrZ9PV28+P1qvIhohIsdDAS0RKgplNBL4LfMLd2wfbNaYtMkJ29wfcfYW7r2hoUNW94dKfdr7dvJe3XjqV2fXjkw6nYF0zp44FDRP41lqlG4qIFIuiHHiZ2V3hRfJpM1uRdDwyhsSlGvb3x2zpmC1nn3yv55IRZ2YVBIOuf3X374XNBzMphOHtobC9BcieamkCNK0wSp7ddph9bae5+7o5SYdS0MyMe1bO4dU9J9i0ry3pcEREJA9FOfACNgK/CPw06UBkjIkpwJHPdWB4Oj6NMHfT2l6jzoJVqr8CvO7uf5n11Brg3vD+vcCjWe0fCqsb3gC0ZVISZeT9y4u7mTKhktuunJ50KAXvrmtnU12R4p9f2J10KCIikoeiHHi5++vuvjXpOESkKNwI/DrwTjNbF27vBr4A3GZm24DbwscAjwE7gO3APwK/m0DMJWn30U6e3HKIX71+DpXlRfnnaVTVjq/gzqsa+bd1rbSd6k06HBEROY/ypAMQKUYeUxTDUpq9KkTu/hzx120B3BKzvwMfHdGgJNbXnt9Fecr4tRsuSTqUovHrb7mEh9fu5duv7OU33zY/6XBERGQQBXtK0cx+YmYbY7bVF/g6ZyqP9dI9UuHKWJZn+mFcRcR8NqUfisDJ0718u7mF9y6bxbSa6qTDKRqLZ9Vy3dzJfO35XfT1q+8QESlkBTvj5e63DtPrPAA8AFBj9apuIMMjtlDGMH7pUSEOKTHfWruXju4+fuPGuUmHUnTue/sCfuvrzfxg/T4+cHVT0uGIiMg5FOyMl4iIlIbTvf088NMd3DC/nmVNdUmHU3RuWTSNy6ZP5MvPvElaawOKiBSsohx4mdkHzKwFeAvwIzN7POmYRCLpiCKSl0ea93LoZDe/d8vCpEMpSqmU8V9uWsAbBzt4csuh8x8gIiKJKMqBl7t/392b3L3K3ae7+x1JxyQSEXdtmIgM0N3Xz5efeZPr5k7mLfOnJB1O0Xrfslk0TR7H3zy1TWsGiogUqKIceIkUrdhCHRqgSel6ZO1e9red5vduWUiw5JpcjPKyFB+/ZSHrW9p4bMOBpMMREZEYGniJiEgiTp7u5a9/so3r59Xz1kunJh1O0fvFa5q4bPpE/vw/ttKrCociIgVHAy8REUnEl595k6OdPfzxe67UbNcwKEsZ//2ORew80skjzXuTDkdERHJo4CUiIqOu5fgpvvLcTj5wdSNLm2qTDmfMuPWKaaycW8+fP76V4509SYcjIiJZNPASEZFR5e586tFNpMz45B2XJx3OmGJmfPbOxbSf7uNPf7wl6XBERCSLBl4iIjKqfrB+P09tOcQn77icxrpxSYcz5iyaUcNH3jqPh9fu5ZXdx5IOR0REQhp4iYjIqDnW2cP/XLOJ5U21fPgX5iYdzpj18VsW0lg3jk9+ez2d3X1JhyMiImjgJSIio8Td+eS3X+Pk6T6+8EvLKEupoMZImVBVzl/88nJ2He3kT360OelwREQEDbxERGSUfOW5nTy15RB/9J4ruGJmTdLhjHk3zJ/C77xjAd98eS+PbdifdDgiIiVPAy8RERlxL+44yhf+fQt3LJ7Oh95ySdLhlIz/69bLuHpOHf/tkdfYtK8t6XBEREqaBl4iIjKith86yX1fb2bu1An8/x9crjW7RlFleYp/+PVrqRtfwW891Myhk6eTDklEpGRp4CUiIiOm9UQXH/7qWirLU3z1w9dRO64i6ZBKzrRJ1fzjh1Zw/FQvv/5PL3NM63uJiCRCAy8RERkRe46e4pf//gXaunp58MPXMbt+fNIhlawljbX8070r2HW0k1/7p5c0+BIRSYAGXiIiMuw2tLTxy//wAp09fXzjN29gWVNd0iGVvBsvncoDH1rB9sMd/OLf/YydRzqTDklEpKRo4CUiIsPq337eygf//nnKUsY3f+sGljbVJh2ShN5xWQPf/K3raT/dxwf+7mc8veVQ0iGJiJQMDbxERGRYHO/s4b9+8+d84lvrWN5Ux6Mfu1Fl4wvQtZfU8/3f/QVm1FTzG19by6cf3ahFlkVERkF50gGIiEhx6+1P882X9/DXP9lGe1cvv3/bZfyXmxZQUaZze4XqkikT+LeP3sif/ngLX/3ZLh7fdJA/fPci3rdsFiktbC0iMiI08BIRkYvS2d3Ht5v38pWf7WTvsS5umF/Pp967mCtnaZarGFRXlPHp9y3mvctm8Zk1m/j4w+v4m6e28zvvWMB7l82kuqIs6RBFRMYUc/ekYxg1NVbv19stSYchInl4yZ+k3Y8Vxan3FStWeHNzc9JhjIqunn5e2HGEH67fz483HuBUTz/XzKnjozdfyjsXTdMaXUWqP+08tmE/X3pqO1sPnqSmupz3XzWLVYtncv38+jE3e2lmr7j7iqTjEJHSohkvEZEYZrYK+F9AGfBP7v6FhEMadem003qii9f3t7OxtY3m3cd5ZfdxuvvSTKou5/3LZ3HXitlce8nkpEOVISpLGe9bPov3LJ3JizuO8vDavXznlRb+5cU9jK8sY8XcelZcMpmljbVcPmMSM2urNcgWEblAGniJiOQwszLgb4HbgBZgrZmtcffNyUY2UCZjwR3S7jjhbfi4Px1svf1OXzpNX7/T3Zemu6+f0739nOrpp7O7n5One2nr6uVYZw9HO3o40H6afSe62H3sFD19aQBSBlfOquFXr7+Emxc1sHJePVXlSkUba1Ip4xcuncovXDqVrp5+nt12mOe2H+GFN4/yl28cPrPf+MoyZk8ez6y6aqbXVDN1YhWTJ1RSO66CSdXlTKwqZ1xlGdXlZVRVpKgsS1FRlqK8zChPGWUpI2XBrRmkLHhsgBka1InImKSBl4hI1Epgu7vvADCzh4HVwJAGXn/yw8189fldee+fmwruZ9qHEsW5laeM+gmVzKitZt7UCdy8aBpzp0xg0cxJXD59EhOq9CejlIyrLOP2xTO4ffEMAE6e7mXzvna2HergzcMdtBzvovV4Fxv3tXO0o5v0CF65kBmH5Q7Hcgdoz/3BzcysHTdygYiIDIH+ioqIRDUCe7MetwDXZ+9gZvcB9wHMmTMnrxe98dKpF1ywIPfEv2U9YVnPZ2YLUuEMQpmdnVEoLzMqwhmHijKjuqKMqvIUE6rKGV9ZRk11RbCNK9dMg5zTpOoKrp8/hevnT4k8l047J0/30dbVy8nuXjpO99HV209XTz89/Wm6+4IZ18zMa3/ag1nZzAxt2kk7OH7mxMKZcVxmZjfnPeNOQOjkgIgUMvVQIiJRcaOPAV/z3P0B4AEIimvk86I3L5rGzYumDT06kQKTShm14yuoHV+RdCgiIgVrbJUpEhEZHi3A7KzHTcC+hGIRERGRMUADLxGRqLXAQjObZ2aVwN3AmoRjEhERkSJWlAMvM/szM9tiZuvN7PtmVpd0TCIydrh7H/Ax4HHgdeARd9+UbFQiIiJSzIpy4AU8ASxx92XAG8AfJhyPiIwx7v6Yu1/m7gvc/fNJxyMiIiLFrSgHXu7+H+EZaYAXCa6/EBERERERKUiWu05MsTGzHwDfcvd/OcfzZ0o+A0uAjaMV2wiYChxJOoghKPb4ofh/hmKK/xJ3b0g6iHyY2WFgd567F9P/QRzFnyzFPzyKpn8RkbGjYAdeZvYTYEbMU3/k7o+G+/wRsAL4Rc/jBzGzZndfMbyRjh7Fn7xi/xmKPf6xoNj/DxR/shS/iEjxKth1vNz91sGeN7N7gfcCt+Qz6BIREREREUlKwQ68BmNmq4A/AN7h7qeSjkdERERERGQwRVlcA/gSMAl4wszWmdnf53ncAyMY02hQ/Mkr9p+h2OMfC4r9/0DxJ0vxi4gUqYK9xktERERERGSsKNYZLxERERERkaKhgZeIiIiIiMgIG5MDLzN70MwOmdnGrLbPmFlreE3YOjN7d5IxDsbMZpvZ02b2upltMrOPh+31ZvaEmW0LbycnHeu5DPIzFMX/g5lVm9nLZvZaGP//DNvnmdlL4f/Bt8ysMulY4wwS/9fMbGfWv/9VScdaiszsz8xsi5mtN7Pvm1ld0jFdCDO7K/y9SptZ0ZQGN7NVZrbVzLab2f1Jx3Mh4v6uFZNz/U0QESklY/IaLzN7O9ABfN3dl4RtnwE63P3Pk4wtH2Y2E5jp7q+a2STgFeBO4MPAMXf/QvilYbK7/0GCoZ7TID/DL1ME/w9mZsAEd+8wswrgOeDjwO8D33P3h8OiLq+5+5eTjDXOIPH/DvBDd/9OogGWODO7HXjK3fvM7E8BCvWzHMfMrgDSwD8An3T35oRDOi8zKwPeAG4DWoC1wD3uvjnRwPIU93etmJzrb0Kx/PuLiAyHMTnj5e4/BY4lHcfFcvf97v5qeP8k8DrQCKwGHgp3e4hgIFOQBvkZioIHOsKHFeHmwDuBzKClYP8PBolfCoC7/4e794UPXwSakoznQrn76+6+Nek4LtBKYLu773D3HuBhgj61KIzhv2siIiVjTA68BvGxMLXnwUJO08tmZnOBq4GXgOnuvh+CP2LAtOQiy1/OzwBF8v9gZmVmtg44BDwBvAmcyPrC3EIBf3HIjd/dM//+nw///f/KzKoSDFEC/xn496SDKAGNwN6sxwX9+R3LYv4miIiUhFIaeH0ZWABcBewH/iLZcM7PzCYC3wU+4e7tScdzMWJ+hqL5f3D3fne/imA2YiVwRdxuoxtV/nLjN7MlwB8Ci4DrgHqChchlBJjZT8xsY8y2OmufPwL6gH9NLtJ4+cRfZCymrWA/v2PVWPi7JiJyscqTDmC0uPvBzH0z+0fghwmGc17hdTnfBf7V3b8XNh80s5nuvj/Mlz+UXITnF/czFNv/A4C7nzCzZ4AbgDozKw9nvZqAfYkGl4es+FdlXVvXbWZfBT6ZXGRjm7vfOtjzZnYv8F7gFi/Ai23PF38RagFmZz0uis/vWHKOv2siIiWjZGa8woFKxgeAgq0MFRZG+Arwurv/ZdZTa4B7w/v3Ao+Odmz5OtfPUCz/D2bWkKk0Z2bjgFsJrkl4GvhguFvB/h+cI/4tmX//8P/nTgr033+sM7NVBLON73f3U0nHUyLWAgvDyqSVwN0EfaqMgkH+romIlIyxWtXwm8BNwFTgIPDp8PFVBKklu4DfzlwvVWjM7K3As8AGgsphAP+DIB/+EWAOsAe4y90L8mLrQX6GeyiC/wczW0ZQPKOM4ATFI+7+WTObT3BRfj3wc+DX3L07uUjjDRL/U0ADQdrVOuB3sopwyCgxs+1AFXA0bHrR3X8nwZAuiJl9APgbgt+lE8A6d78j2ajOz4LlK/6a4HPxoLt/PuGQ8hb3d83dv5JoUBfgXH8T3P2x5KISERldY3LgJSIiIiIiUkhKJtVQREREREQkKRp4iYiIiIiIjDANvEREREREREaYBl4iIiIiIiIjTAMvERERERGREaaBl4iIiIiIyAjTwEtERERERGSEaeAlIiIiIiIywjTwEhERERERGWEaeImIiIiIiIwwDbxERERERERGWKIDLzNbZWZbzWy7md0f8/zvm9lmM1tvZk+a2SVZz91rZtvC7d7RjVxECp36FxERESkk5u7JvLFZGfAGcBvQAqwF7nH3zVn73Ay85O6nzOy/ADe5+6+YWT3QDKwAHHgFuNbdj4/2zyEihUf9i4iIiBSaJGe8VgLb3X2Hu/cADwOrs3dw96fd/VT48EWgKbx/B/CEux8Lvww9AawapbhFpPCpfxEREZGCkuTAqxHYm/W4JWw7l48A/36Rx4pIaVH/IiIiIgWlPMH3tpi22LxHM/s1grSfd1zEsfcB9wGUUXbteGouPFIRGXWn6aTHu+M+6/kY8f5FfYtI8TrJ8SPu3pB0HPmYOnWqz507N+kwRCRPr7zyyjn7lyQHXi3A7KzHTcC+3J3M7Fbgj4B3uHt31rE35Rz7TNybuPsDwAMANVbv19stQ41bREbBS/7kUA4f8f5FfYtI8fqJf2d30jHka+7cuTQ3NycdhojkyczO2b8kmWq4FlhoZvPMrBK4G1iTvYOZXQ38A/B+dz+U9dTjwO1mNtnMJgO3h20iIqD+RURERApMYjNe7t5nZh8j+EJTBjzo7pvM7LNAs7uvAf4MmAh828wA9rj7+939mJl9juDLFcBn3f1YAj+GiBQg9S8iIiJSaJJMNcTdHwMey2n7VNb9Wwc59kHgwZGLTkSKmfoXERERKSSJLqAsIiIiIiKja++xU5zu7U86jJKjgZeIiIiISAnoTzt/8+Q23vFnT/Ohr7yswdco08BLRERERGSMO3yym3sffJm/eOINrp83hZd3HeP3H1lHOh272oqMAA28RERERETGsPUtJ3j3F59l7a5j/OkvLeUbv3U9f/yeK3hswwE+96PNuGvwNRoSLa4hIiIiIiIj6/97bAsGrPnYW7l8xiQAfvNt89l34jQP/mwns2rH8Vtvn59skCVAM14iIiIigzCz2Wb2tJm9bmabzOzjYXu9mT1hZtvC28lhu5nZF81su5mtN7Nrsl7r3nD/bWZ2b1I/k5SOHYc7eGHHUe79hblnBl0Zf/yeK3jP0pl8/rHXefNwR0IRlg4NvEREREQG1wf8N3e/ArgB+KiZXQncDzzp7guBJ8PHAO8CFobbfcCXIRioAZ8GrgdWAp/ODNZERsq31u6lPGXctaIp8lwqZXzqfVdiBo+u25dAdKVFAy8RERGRQbj7fnd/Nbx/EngdaARWAw+Fuz0E3BneXw183QMvAnVmNhO4A3jC3Y+5+3HgCWDVKP4oUmK6+/r59ist3HrFdKZNqo7dZ3pNNW+ZP4U161p1rdcI08BLREREJE9mNhe4GngJmO7u+yEYnAHTwt0agb1Zh7WEbedqFxkRT2w+yLHOHu65fs6g+915VSO7jp5ifUvbKEVWmjTwEhEREcmDmU0Evgt8wt3bB9s1ps0Hac99n/vMrNnMmg8fPnxxwYoA33hpD41143jbpVMH3e+OJTOoLEsp3XCEaeAlIiIich5mVkEw6PpXd/9e2HwwTCEkvD0UtrcAs7MObwL2DdI+gLs/4O4r3H1FQ0PD8P4gUjJ2Hunk+TePcs/K2aRScWP+s2rHVXDzogZ+sH4f/VrXa8Ro4CUiIiIyCDMz4CvA6+7+l1lPrQEylQnvBR7Nav9QWN3wBqAtTEV8HLjdzCaHRTVuD9tEht3Da/dQljLuWjH7/DsDq69q5PDJbl548+gIR1a6tI6XyPnY4GeJzu6Xcx7D0xf/nrq4VUbSuX6nM7/Dg/3u6ndTStONwK8DG8xsXdj2P4AvAI+Y2UeAPcBd4XOPAe8GtgOngN8AcPdjZvY5YG2432fd/djo/AhSSnr60nynuYVbFk1jek18UY1c71w0jYlV5Ty6rpW3Lhw8NVEujgZeIiIiIoNw9+eIvz4L4JaY/R346Dle60HgweGLTiSqefcxjnb28EvXRkvIn0t1RRmrlszgxxsP8Lk7l1BdUTaCEZYmpRpKaTMbuKXKIpuV5bmlbOAWsw+Wim4io8190JmrQX/Py8ux8vKzn5Hcz5CIiCTuxTePkjJ4y4IpF3Tc6qtmcbK7j6e3HDr/znLB9K1PRERERGQMeWHHUZY21lJTXXFBx71l/hSmTqzih+v3j1BkpU0DLxERERGRMaKrp591e09wwwXOdgGUl6V4+2VTeXHHUS2mPAISHXiZ2Soz22pm283s/pjn325mr5pZn5l9MOe5fjNbF25rRi9qKQq56U/5phFWlEe38pitsjK65exDWVlki6Qjpiw+/VDpW0Om/iUPmZTDM1saPI2n/cx25ncy8zsc/n6nqquCrSrYznwOKoItkoooIiKjonn3MXr7nbfMv/CBF8DKufUc7exhx5HOYY5MEiuuYWZlwN8CtxGsa7HWzNa4++as3fYAHwY+GfMSXe5+1YgHKiJFR/2LiIiUqhfePEp5yrhubv1FHb9yXnDcyzuPsaBh4nCGVvKSrGq4Etju7jsAzOxhYDVw5ouRu+8KnxtCXW4Z82LOpltZTCWemDYrH/gRsIqYj0S+BTD6+wc+7uuL7OKxixL2x7TJEKl/uRhn0krO/pN4+OtphJ+fsvDzFn6erDK4fsByPyfh58Ezt9mfh5znVKJeRGT4PP/mUZY11TKh6uK+5s+bOoGpEyt5eecx7lk5Z5ijK21Jpho2AnuzHreEbfmqNrNmM3vRzO48105mdl+4X3Mv3Rcbq4gUlxHvX9S3iIhIoeno7mNDa9sFVzPMZmasnFfPyzu1xNxwS3LgFZf0fyGnPee4+wrgPwF/bWYL4nZy9wfcfYW7r6ig6mLiFJHiM+L9i/oWEREpNGt3HqM/7bxl/tAWQL5ubj2tJ7poPdE1TJEJJJtq2ALMznrcBOzL92B33xfe7jCzZ4CrgTeHM0ApUDmphXFphVZZGdMWU1K1auAXZquI2SeuMEA6mp3m3TmzHnHpU7npiKBMw5Gh/mUosn93w1xD9/B33ssG7hN+1ixMaTnzGcqkImY+P1mfGe/pDe709gx47D3hY6UgiohclBd2HKWyLMW1l0we0utkrvNau/MYjVdfSMKIDCbJGa+1wEIzm2dmlcDdQF7Vw8xssplVhfenAjeSde2GiJQ89S8iIlJyXnjzKFfNqWNcZcy17hdg0YwaJlWV85LSDYdVYjNe7t5nZh8DHgfKgAfdfZOZfRZodvc1ZnYd8H1gMvA+M/uf7r4YuAL4h/Ci+BTwhZxqZTJWxBXOKB8APYSuAAAgAElEQVQ4K2XV0TQvq66Otk0YF2nzcQOPTcfNnsWdde/uie6XOwsWU1wj7ufROf3hp/5lBISfgzNFMjKzUrm344L9Mp9Brw5nn6vOzkJ7WKDD+sN9wxmvVFcwa+xdXeHt6eC2ty98794BsYiIyFltp3rZuK+N33vnwiG/VlnKWDF3Mmt3aeA1nJJMNcTdHwMey2n7VNb9tQQpQrnHPQ8sHfEARaRoqX8REZFS8tLOo7gzpMIa2VbOm8LTW7dwtKObKRN1LfNwSHQBZRERERERGboXdhylqjzF1XPqhuX1Vs4LrhNbu+v4sLyeJDzjJTJAHmmFAKlxA9MIbVJ0cT+vmRBp66uNphr2jxv4EfCYGMpOR1MGyzpiYu3N2e/0EEqMK5VKCl1O6mFmjbrM5yCTPmhx69lNGg9AX03w+U5XDvwMZz5zZZ1BSm9Z+6nguJMdwW0mBTEsaOPZRWv02RGREtW86zhXz6mjqnxo13dlLG2so6o8xcs7j7FqyYxhec1SpxkvEREREZEi1t3Xz5YD7SyfPTyzXQCV4ezZy7uODttrljrNeElycsvC5zG7BWC1NQMe90+tjezT0xCd3equjZ4B6q8cGEMqph5GVVv0/ESqJ6YGfFzZ+Ryus/EyVqXDsvOZ23AWysKZ31Q4SwVQFhbRSNVNAqC3PpgB664P+oBT04NbLws+x2U9wWe8sj14zaqjwWuljgczYLR1nHlt7+wMblWaXkRKyOv7T9Lb7yxvGr6BFwTXeX3pqW2cPN3LpOqYJXfkgmjGS0RERESkiK1vOQHAsqboyeihWDm3nrTDq3tODOvrlioNvEREREREithre9uYMqGSxrpoxs9QLJsdDOQ2tGjgNRyUaiijI4/CGbFphZOjZ276pg+cRj/VOD6yT+f06DmF3przpwJWtEfTkcp6oq9VmYp5rdxUpv6YdMR0TLqTp6NtIsUuJ/WwP+vzYGFRjFSYFljZHqQclnUEn/fTM4LP9KnpwZ+orobg89a3IPgspnorwuOCohzjD5/9DI07EKQhlh9uDxqOtwVxdAYFOtI9mbXAwmOUgigiY8D6lhMsa6rF8rjs4ULUVFcwf+oE1re0DevrlirNeImIiIiIFKmO7j62H+5g2TBf35WxpLGWja0aeA0HzXjJ8Iub3SqLFraw6oGL8cXObs2cHGk7OXfgDFf7JdHzB13To7NI6eroDFRZ58C4rD8au8dUZbW4mau+ga/vcTNemt2SUpU++3nw7nAWLCyAYZ1dAKTaglmqCceCvqD6SPAlonN28Jk/2RR8GE/NDD5/nfOCz9MxO/t5rGgP9h13IFhSYlLLtOA1W4L3KD8QpMv4ieBLRDpTmj6MRTNgIlJsNra24Q7LZw/v9V0Zy5pqWfPaPg6f7KZhkhZSHgrNeImIiIiIFKmzhTVGZsZraWMwoNOs19Bp4CUiIiIiUqTWt7TRWDeOqRNHZjZqcWMtZug6r2GgVEMZfhYdz1tVtDNI1Uwa8Di3aAZA+/xo4YwTlw18/e75pyP7TKnviLT19EVzBtsPDIzBD0U/EqneaOqRnY4u+OWZi/YzYlINPba4hlKbpESFv/veG6T59fcFnyHryqzTFfyRrzkQrN03viVIPe6YG6QRnlgQph5ecvbzWL0oOPM7fnnwmsc6gwpf+/YHx0zcGRTkqNsxI9hvd1Dgo/zQcQDSYQqihwVAQGuBiUhhW9/SNuxl5LNNrCpn/tQJbNCM15BpxktERERkEGb2oJkdMrONWW2fMbNWM1sXbu/Oeu4PzWy7mW01szuy2leFbdvN7P7R/jlk7Dne2cOeY6dGLM0wY1lTHRtaVVJ+qDTjJUOTis4ipSqjK5vbxAmRtnTDwE4icxY72/HLo+cGUovbBzy+uWlXZJ/6ys5I26a2mZG29kMTBzwu647sQnlndObKTsfsGJ61z/D+mEIaKq4hcm7nmgE7FZSCTx0LZqXqWoMzu5N21ANwcv7Zz/Hxy4JZsUOLgmIay+e0ALBowUEAOt4SzL6vPTwHgDe2TwWgdkvwmpO3Be9dvTvrC8bR4H29I+hXzpSkT8cU0JGx6mvAl4Cv57T/lbv/eXaDmV0J3A0sBmYBPzGzy8Kn/xa4DWgB1prZGnffPJKBy9i2PpyFWj6CM14QVDb8/s9bOdR+mmk10eV/JD+a8RIREREZhLv/FDiW5+6rgYfdvdvddwLbgZXhtt3dd7h7D/BwuK/IRVu/NzhJtGSEB16ZVEalGw6NBl4iIiIiF+djZrY+TEXMrH/SCOzN2qclbDtXu8hFe62ljfkNE6ipjmYbDacrZ9aQUoGNIUs01dDMVgH/CygD/sndv5Dz/NuBvwaWAXe7+3eynrsX+OPw4Z+4+0OjE3WJy1mjK3Z9rnHjosfVR8/EnGoamOZ34tLzpxUC3LXw5wMev33ilsg+rb3R9b/iUg1TOet4VbZFL56vaO+JtGXWHcqWzimuEb+Oly7OHy3qX8aATOphX9+A2/SZ4hvBmd7avWf7kpotQfphx6agz9m6aCEA65cG329XLXwdgD9e+CMAJlwWfL6ff0ew33d3LwegZWPDmdecvDm4X7ctLMTRGkx8nCnE0RX0ByrCUXK+DHwO8PD2L4D/DEQXhAz2iTvZHfvLYmb3AfcBzJkzZzhilTFqfcsJfmHBlBF/nwlV5SxomKgZryFKbMbLzMoIcp3fBVwJ3BPmRWfbA3wY+EbOsfXAp4HrCabuP511pklESpz6FxEZae5+0N373T0N/CNBfwHBTNbsrF2bgH2DtMe99gPuvsLdVzQ0NMTtIsKBttMcOtk94oU1MpY21bKhtQ3XyaWLluSM15lcZwAzy+Q6n7nI1N13hc/lViS4A3jC3Y+Fzz8BrAK+OfJhl7bcGS6rjpaJt5qJkbbu6ZMibe1zB/76dV0eLQt/94L1kbYP1b004PGkVPTk4ps90yNtrW3RWbeqIwPPPYw7Gi0TX9YWnd3ymOIa3ptzrAppJEn9y1gWFrVIh0U3yPo8WjgLNak1KNYzaXMwA9b182Ds/OyiawH40fIlALxrySYAPjL1pwB84uoNAGxfevbX4sGjbwXgB5uXAjD+tWAGYsrmYJZ7/M6wEMfhcCbsZLCcxZk+QUU4xiQzm+nu+8OHHwAyFQ/XAN8ws78kKK6xEHiZYCZsoZnNA1oJCnD8p9GNWsaSzMLJy2eP7PVdGcsaa/neq60cbO9mRq0KbFyMJAdecbnO1w/hWOVJi0iG+hcRGTZm9k3gJmCqmbUQzIrfZGZXEaQL7gJ+G8DdN5nZIwQnevqAj7p7f/g6HwMeJ0iBftDdN43yjyJjyMZ97aQMrphZMyrvtzQssLG+5QQzameMynuONUkOvM6VAz2sx2bnSVcTXYxXRMakEe9f1LeIlA53vyem+SuD7P954PMx7Y8Bjw1jaFLCNu9rY37DRMZXjs7X+Stn1pKyoLLh7Ys18LoYSQ688s51PsexN+Uc+0zcju7+APAAQI3VKyn1Qlj0+6eVD/yVsfHRQhr99dEzL52zKiNtJ+cNzPC6at7eyD531r0SaVtQMTCV8Y3e6Jpda9vnRdra90XTHafuH/grUX0omlZo7R2RtnR3TKphbjEN5UAnacT7F/UtBSQrlc+7g/v9PUHRDGsLCvSM238IgKbXgzO2vWuDi9HXXn41AD9efhUAS67aBcDHm54485p/PbMZgD+bEaQ5f//6IH3xiztuAWD7uiC1uX5D8Jp1b4RFOPYdDcLLFOHI6jdUiENEhmpjazs3zK8ftfcbV1nGZdMnqcDGECRZTn4tYa6zmVUS5DqvyfPYx4HbzWxyeNH77WGbiAiofxERkTHsSEc3B9pPs3jW6FzflbGksZYNLSqwcbESG3i5ex+QyXV+HXgkzIv+rJm9H8DMrgtzqe8C/sHMNoXHHiMo3bo23D6buRBeRET9i4iIjGWb9gWz+YsbR+f6rowls2o42tnDwfZo5o+cX6LreMXlOrv7p7LuryVI84k79kHgwRENsJTEpRXGrdFVOTBl0CZEr23paYimH3Y0Rsf4ZU0DU/jeWv9mZJ+55dE1tA71D2x7rGNxZJ9nd8+PtE3cGf11n9g68LXKjkXTCr3zVLStJxqXqhgWFvUvJS5nDbD+9nBNwI4gDbDscJAG2LA9SF2e+kqQJnh8wSUA/Ncrf/vMS9nVQVrNhy9/EYB7al4DYPXSbwHQfFnQV/7t24LUwxc2BGuCTX4teK36LcEXlKrdZ8fvfuw4AOlwTUDv6x0Qt4jIYDbtC/qlxTNHd8ZrcWPtmfdXZcMLl2SqoYiIiIiIXKBN+9ppmjyO2vEVo/q+V8yswSy4vkwuXKIzXlJALDoGzy2kAcC4gWc3+msnRHbpmho97nRD9CzunCknBjyeWh79EO/qixbleLVr7oDHD+9eEdnH34iuJVazK7qWTvX+nBmutosspAE6Uy1SDHLXAOsKZpzsWNAfTdodzNjXrDu7Znb3c8HF6w9fdjsAX15yKwALF7cCcE/jywD837N+DMCB6c8D8O2rrgPg6S2XAzBxw6wzr1m/JVgUd/zO8CL1Q8EMnIczcume3gHxiohk29TaxpJRvr4LYGJVOfOmTDgz4yYXRjNeIiIiIiJF4uTpXnYdPcXiWaN7fVfGlbNqzlxjJhdGAy8RERERkSKxORz0LGkc/RmvzPu2nujieGfM9e4yKKUalqK4QhqpmLbKaJpf7rpdfXVVkX1O10XH8301fZG2qrKBbTu7p0X22dQVrX3w0wMLBjw+sqkhss+UNyJNTNwdLZKROjrwjE26M7ommPdGY1daocgYkSnC0Rt8geg/EX6RyFq/r+JAsAbY9NeDFOZpLwaph53zgtTBv7jsgwB0XBEcu3RBCwArJ+8CYME1hwFonj/nzGuuWxYsMzf+9akATN4apDZO3HkSgPID4Rpg7cHjTEEfrf8lImcqGiY045V5383727nx0qmJxFCsNOMlIiIiIlIkNu5ro2FSFdNqkqkqmFk7bKMWUr5gmvGSQEzpeCqjlXJ83MAZrt4J0V+h/mg1+ViHOwcWwHimb2FknwMnomdz+nYOPK5ua/S1a3aejrSVH4p2EJmzyRkqEy8iwICiFucqxDFxV/ClZ9K64EtI76xgJuzw3LkAfGNBsKxF1/ygX5k+42xBoQVNwSzY0cnBkhwtc4PXmLijDoDaNyeF7xHMvJUdCMvPZ2bAwlg0AyZSejbva09stgugfkIls2qrdZ3XRdCMl4iIiIhIETjd28+2Qx2JVDTMduWsWlU2vAgaeImIiIiIFIGtB07Sn/ZEZ7wAljTWsONIJ53dMdfByzkp1bAUxa3ZFZNqaBXRVMP+qoFt6cpoUY44qVPR9zxyYGCncaQ3uk/1weivaN3ugSk1tTuj62xVth6PtPmJ6JkZPz3wWK3PJSLnlFuII7wlXHur7GCQPjj5jSB9cHJ9kDbYG56Zbp9zthDQvqag7+yaEaYyTwy+vHQsCNbvOl0f9H0dYdWymj3Bmonj92beK0hb9DD1MN11Nr3a+3oHxCsiY8fZwhrJzngtnlWLO2w50M61l9QnGksxyWvGy8wWmtl3zGyzme3IbCMdnIiIiIiIBDbua2NSdTmz6/O8oH6ELGkMTp5vbNV1Xhci3xmvrwKfBv4KuBn4DSC/qQ5JXk75+NjS8eUxvwqxbTlj9ZgTqmVd0bbqI9Exvh8bWK6+oiOyC+MPRgtbTGwZOEtVsf9EZB8/HlNIoytacOPMmeEzDTpDLCIXKCzEkT4dzph3B32UtQVfSMr3BUWJpmyZcOaQKfXB2eqeGUERjc4ZQX94alrQV/aGNYR6wpPaxy8N+uOuycGXnQkHglm1cfvDGbAjZ/u83FkwzYCJjB2bwsIaFrM00GiaUVNN/YRKXed1gfK9xmucuz8JmLvvdvfPAO8cubBERERERCSjtz/Nlv3tiRfWADAzFs+qUWXDC5TvwOu0maWAbWb2MTP7ABBd7VZERERERIbdtoMddPelWdqU/MALguu83jh4kp4+LbuTr3xTDT8BjAd+D/gcwWzXvSMVlIywuDW7YtIPKYsZl+ekqpSdjn7Yqk9EX6u8K9pW1jvwtaraYl7rcNx6XAPPrnhb9GzLmXV3svfrjam8o9QbERlumSIcfX0DbsnqlzJrgVW0BGmI9RODNMTJYSph75TgcfeUoKBR96Sg306H9Y26pgZ/vvurg1TF6tqzayxWHA6OLTse9I3psPhHppiQUg9FilNmweIljYUy8Kqht9954+DJgomp0OU18HL3teHdDoLru0REREREZJRsaG1jYlU586ZMOP/OoyAz2Nq0r00DrzzlNfAys8uA/w5ckn2Mu+s6r2IQUz4+r31izoZaz8BZo4qO6CxSqj96nMdMqJWfGli6vbw9OruVOhGtuOHtA9vynt1Kx5SKFxEZLVl9aqYkfeY23Rn0Y3bkGAAVLdUAVI4PKpdNrAmqbaRrg6IafRODqa/+yqDv7qvOymRoCGe8wuU/UieD17CT4cxXZgasJ5j50gyYSHHY0NrGlbNqSMVlKSXgkvrxTKwqZ2NrO79yXdLRFId8r/H6NvAq8McEA7DMNiRmtsrMtprZdjO7P+b5KjP7Vvj8S2Y2N2yfa2ZdZrYu3P5+qLGIyNii/kVERMaKvv40r+9vZ2kBzSylUsaSxhrWt6qyYb7yvcarz92/PJxvbGZlwN8CtwEtwFozW+Pum7N2+whw3N0vNbO7gT8FfiV87k13v2o4YxKRsUH9i4iIjCXbDoWFNQpo4AWwrKmOrz2/i56+NJW5Sw5JxKADLzPLLEX9AzP7XeD7wJlFlNz92BDeeyWw3d13hO/1MLAayP5itBr4THj/O8CXLOmFCwrdxf7zeExFmr5oap6dHrjuVdwvUFlHtHiH9ce8VlfPwMcdMSmDndG2dPfAdbyUVigx1L9IcQn7LA9v+zPpf2FaoB0PinFYZbDeVyYF0cYFKYk+vvrMS3ll0DN7+Ovs44PCG5YKvhSlKoIURO/qCm/D9b56wvTHTH+t1EORgrEhnFUqlIqGGUsba+npS6vARp7ONzR9BWgmqGD434Hnw7ZM+1A0AnuzHreEbbH7uHsf0AZMCZ+bZ2Y/N7P/bWZvO9ebmNl9ZtZsZs29dJ9rNxEZW0a8f1HfIlI6zOxBMztkZhuz2urN7Akz2xbeTg7bzcy+GKYxrzeza7KOuTfcf5uZqTq05G1jgRXWyFgWDgQ3KN0wL4POeLn7vBF877gzy7mn1861z35gjrsfNbNrgX8zs8XuHqkr7u4PAA8A1Fi9Tt+dSzqmIEZPb6Qtd0LAYvaJPUsaM3vmOTNX6a5ocY24GDx39kyzWxI14v2L+hYZUZl+1MOZsO7wNpyVIlOMI1wexCorzhxq1cEMVyqcHaM8/FOfWSIkvLXqcJYs069nlhrJvEdWX6tZML4GfAn4elbb/cCT7v6F8DrS+4E/AN4FLAy364EvA9eHWUSfBlYQ9DWvhCnQx0ftp5CiVWiFNTLm1I+nprqc9S1t3LMy6WgKX17JmGZWbWa/b2bfM7PvmtknzKz6/EcOqgWYnfW4Cdh3rn3MrByoBY65e7e7HwVw91eAN4HLhhiPiIwd6l9EZNi4+0+B3MsrVgMPhfcfAu7Mav+6B14E6sxsJnAH8IS7HwsHW08Aq0Y+eil2ff1pNu8rrMIaGWbGsqY6NrSeSDqUopDvVXBfBxYDf0NwxudK4J+H+N5rgYVmNs/MKoG7gTU5+6zh7ELNHwSecnc3s4bw4nnMbD7BWaUdQ4xHRMYO9S8iMtKmu/t+gPB2Wth+rlTnfFKgRSIKtbBGxtKmWrYeOMnpXmUgnU++VQ0vd/flWY+fNrPXhvLG7t5nZh8DHgfKgAfdfZOZfRZodvc1wFeAfzaz7QRnmu4OD3878Fkz6wP6gd8ZYqGP0hKXVtgXLVBhMdPZnpvWF1eLoD9aqCPu9entPe8+kbRCKOVUF8mT+hcZs3JTEDNFOfqy+tMwbftMGmImxTAsqmG5qYehM+0Z2f1vWMTobMphemA8ku1cacz5pEAHL2B2H3AfwJw5c4YvMilKmeunCrV4xbLGWnr7na0HTrJ8dl3S4RS0fAdePzezG8Ipc8zseuBnQ31zd38MeCyn7VNZ908Dd8Uc913gu0N9fxEZu9S/iMgIO2hmM919f5hKeChsP1eqcwtwU077M3EvnH0N6YoVKzS6LXEbW9uYUFnG/KmFVVgjI1NpcX1rmwZe55HvwOt64ENmtid8PAd43cw2AO7uy0YkOhkZMaXjY2eWeqJNWM6sVNzZznTc6+fxnnEl7XU2VUTk/LL7yshsWDjREhY0ysyEYWGRjczMV1huPi6TwXJmx/CwVH1pF93IpCt/Ibx9NKv9Y+EyFtcDbeHg7HHg/81UPwRuB/5wlGOWIrShtY3Fs2oLrrBGRmPdOOonVLKh5QRwSdLhFLR8B166+FNERERKkpl9k2C2aqqZtRBUJ/wC8IiZfQTYw9kZ9MeAdwPbgVPAb0Cw9qmZfY7gGlSAzyqNWc6nrz/N6/vb+U8rC3dAY2YsbaxlfYtKyp9Pvgson4x7Xh2GiIiIjHXufs85nrolZl8HPnqO13kQeHAYQ5MxbvvhDk73plnaVJN0KINa1lTL3z1zhK6efsZVliUdTsE634zXKwy8IDSTS2Dh/fkjFJdcrNh0j4EpfJ6OFrM0YtbZiinCEbNTfsfFpRFG9inJVBURkZGV07fmFjI6k4qYST1MDXw8qDP75PTx6s9FhsWGcBZpaWNhXzu1tLGW/rSzeX87114y+fwHlKi8F1AOZ78WAkNdv0tERERERM5j3d4TTKwqZ16BFtbIWNYUDAw3tJzQwGsQeV3jZWa/CXycoALPOuAG4HliptilCMTOUsWd2cxjPYZ8ZrJAZz9FRApVpER92B63XEhG7mxY5nHm4Myx6vtFhuTVPSe4anYdZQVaWCNjek0VDZOqWN+q67wGk+8Cyh8HrgN2u/vNwNXAkRGLSkRERESkhHV097H1QDvXFMEMkpmxrLH2TGqkxMt34HU6XPMGM6ty9y3A5SMXloiIiIhI6Xpt7wnSDtfMKezruzKWNtWy/XAHHd1959+5ROVbTr7FzOqAfwOeMLPjBIsBSjHIK9VDKYMiIpJlsP7ec1LRB0tLFJGL8uru4wBcPbvwZ7wArpkzGXf4+Z7jvG1hQ9LhFKS8Bl7u/oHw7mfM7GmgFvjxiEUlIiIiIlLCXt1znIXTJlI7viLpUPJy9Zw6UgbNuzTwOpd8Z7zOcPf/PRKBSMI0kyUiIhdLf0NEhlU67by65wSrFs9IOpS8TaquYNGMGl4JZ+okKt9rvEREREREZBTsONJJW1dv0ZVmXzF3Mq/uOU5ff56XsJQYDbxERERERArIq3uCWaNrLimOwhoZK+bWc6qnny0HTiYdSkHSwEtEREREpIC8uvs4NdXlzJ86MelQLsiKcIauedexhCMpTBp4iYiIiIgUkFf3HOeaSyaTKvCFk3PNqhvHrNpq1uo6r1gaeImIiIiIFIi2rl62HergmjnFdX1Xxoq59TTvOoar6E6EBl4iIiIiIgVi3d4TuFPEA6/JHGzvpuV4V9KhFJxEB15mtsrMtprZdjO7P+b5KjP7Vvj8S2Y2N+u5Pwzbt5rZHaMZt4gUPvUvIiJSjF7dfZyUwfLZtUmHclEylRhVVj4qsYGXmZUBfwu8C7gSuMfMrszZ7SPAcXe/FPgr4E/DY68E7gYWA6uAvwtfT0RE/YuIiBStV/cc57Lpk5hUXRwLJ+daNKOGiVXlNO9WgY1cSc54rQS2u/sOd+8BHgZW5+yzGngovP8d4BYzs7D9YXfvdvedwPbw9UREQP2LiIgUoXTaWbfnBNcU2fpd2cpSxtVz6mjepRmvXEkOvBqBvVmPW8K22H3cvQ9oA6bkeSwAZnafmTWbWXMv3cMUuogUuBHvX9S3iIjIcNu0r52T3X1cN7d4B14AKy6pZ+vBk7R19SYdSkFJcuAVVx8zt/zJufbJ59ig0f0Bd1/h7isqqLrAEEWkSI14/6K+RUREhtuz2w8DcOOlUxOOZGiumzsZd/j5Hs16ZUty4NUCzM563ATsO9c+ZlYO1ALH8jxWREqX+hcRESk6z207wqIZk5g2qTrpUIbkqjl1lKVM6YY5khx4rQUWmtk8M6skuJh9Tc4+a4B7w/sfBJ7yYFGANcDdYVWyecBC4OVRiltECp/6FxERKSpdPf007zrO2xYW92wXwPjKcpY21vKzN48kHUpBKU/qjd29z8w+BjwOlAEPuvsmM/ss0Ozua4CvAP9sZtsJzkTfHR67ycweATYDfcBH3b0/kR9ERAqO+hcRESk2L+08Sk9/mrcubEg6lGHxjssa+OJT2zje2cPkCZVJh1MQEht4Abj7Y8BjOW2fyrp/GrjrHMd+Hvj8iAYoIkVL/YuIiBSTZ7cdobI8xcq59UmHMixuuryB//XkNp7dfoT3L5+VdDgFIdEFlEVEREREJLi+67q5kxlXOTaWjlzWVMfk8RU8s/VQ0qEUDA28REREREQSdKj9NFsPnuRtYyTNEIL1vN62sIGfvnGYdDq2+HjJ0cBLRERE5CKZ2S4z22Bm68ysOWyrN7MnzGxbeDs5bDcz+6KZbTez9WZ2TbLRS6F4dltQhOKtRV5GPtdNlzdwpKOHzfvbkw6lIGjgJSIiIjI0N7v7Ve6+Inx8P/Ckuy8EngwfA7yLoFLqQuA+4MujHqkUpOe2H2HKhEqunFmTdCjD6u2XBTN4SjcMaOAlIiIiMrxWAw+F9x8C7sxq/7oHXgTqzGxmEgFK4XB3nt12hBsvnUoqZUmHM6ymTqxiaWMtz2w9nHQoBUEDLxEREZGL58B/mNkrZnZf2Dbd3fcDhLfTwvZGYG/WsS1h2wBmdp+ZNZtZ8+HD+sI61m05cJIjHd28dQys3xXnpssbeHXPcZTfxOEAABeASURBVNpO9SYdSuI08BIRERG5eDe6+zUEaYQfNbO3D7Jv3HRGpOqAuz/g7ivcfUVDw9gptiDxnguv7xoLCyfHuenyBtIepFOWOg28RERERC6Su+8Lbw8B3wdWAgczKYThbeYClxZgdtbhTcC+0YtWCtFPXj/IwmkTmVk7LulQRsTypjpqqst1nRcaeImIiIhcFDObYGaTMveB24GNwBrg3nC3e4FHw/trgA+F1Q1vANoyKYlSmva3dfHyrmO8d9nYXWC4vCzF2y5r4H+/cRj30i4rr4GXiIiIyMWZDjxnZq8BLwM/cvcfA18AbjOzbcBt4WOAx4AdwHbgH4HfHf2QpZD88LX9uMP7rxq7Ay+Amy5r4NDJbja2lnZZ+fKkAxAREREpRu6+A1ge034UuCWm3YGPjkJoUiQefa2V5U21zJs6IelQRtRtV06nsizF93/eytKm2qTDSYxmvERERERERtmbhzvY2NrO+5aP7dkugLrxldxyxTQeXddKb3866XASo4GXiIiIiMgoW7NuH2aUxMAL4IPXNnG0s6ek1/TSwEtEREREZBS5Oz94bR9vmT+F6TXVSYczKt5+WQNTJ1by3Vdakg4lMRp4iYiIiIiMoo2t7ew40snqMV5UI1tFWYrVVzXy5JaDHO/sSTqcRGjgJSIiIiIyih5d10pFmbFq8cykQxlVH7y2id5+Z81rpbl8XSIDLzOrN7MnzGxbeDv5HPvdG+6zzczuzWp/xsy2mtm6cJs2etGLSCFT/yIiIoWsP+38YP0+brp8GrXjK5IOZ1RdMbOGK2fW8N1XSzPdMKkZr/uBJ919IfBk+HgAM6sHPg1cT7AK/KdzvkD9qrtfFW5aCltEMtS/iIhIwXp222EOtnfz/hIpqpHrl65tYn1LG28cPJl0KKMuqYHXauCh8P5DwJ0x+9wBPOHux9z9OPAEsGqU4hOR4qX+RURECtbfPf0mM2uruWPxjKRDScTqq2ZRnrKSLLKR1MBrurvvBwhv41J5GoG9WY9bwraMr4ZpQP+Pmdm53sjM7jOzZjNr7qV7OGIXkcI2Kv2L+hYREblQL+04ysu7jvHbb59PZXlpllqYOrGKdy6axrea99J+ujfpcEbViP2Pm9lPzGxjzLY635eIafPw9lfdfSnwtnD79XO9iLs/4O4r3H1FBVUX9kOISEEqhP5FfYuIiFyoLz29nakTK7l75ZykQ0nUf33nQk6c6uXB53YmHcqoKh+pF3b3W8/1nJkdNLOZ7r7fzGYCcddQtAA3ZT1uAp4JX7s1vD1pZt8guEbj68MUuogUOPUvIiJSbF7be4Jntx3h/nctorqiLOlwErW0qZZVi2fwT8/u5N63zGXyhMqkQ/o/7d15nBxlncfxzzczuRNCQgIEEkgC4YiIHOEQEMMpuiCgURFkQVEUZFlF1oVl5YXXcrmrriIREIkKqESQwCohYIIoG44AuUiyDBEhEAIhEpIQcs1v/6hnoB1mJj0z3V3dM9/3a+rVXVf396maeqqerurqisjrHOdUoOkuYmcAd7YwzTTgWEmD05fejwWmSaqXNBRAUk/geGB+BTKbWW1w/WJmZlXnhzMaGNS3J586eOe8o1SFC47djbUbNjHpj8/kHaVi8mp4XQEcI+lp4JjUj6Txkm4AiIiVwDeBR1P3jTSsN9kB0lzgSeAF4PrKF8HMqpTrFzMzqyqLXnqd6U8t59OHjmJA77JdcFZTdttuICftsyOTH3qWl19/M+84FZHLmo+IV4GjWhj+GPDZgv4bgRubTbMW2L/cGc2sNrl+MTOzavOD+xvo36uOMw8ZlXeUqvKlo8dy15wXuWZGA18/ca+845Rd97ydipmZmZlZBfx+3jL+Z94yPnf4GLbu1z2+y1Ssnbfpz8fGj+SWR57j+ZVv5B2n7NzwMjMzMzMrg5dWvcnFd8xj7xGD+OIRu+Ydpyqdf9Su1PfowVdum8OmzY15xykrN7zMzMzMzEqssTH4lylzWL+xke99Yh961vmwuyXDB/Xl8o+8m0f+spKr712cd5yy8n+AmZmZmVmJ/fShZ3nw6RX8+/F7MmbYgLzjVLWT9t2R0w7aiR8/sIR7F7yUd5yyccPLzMzMzKyE5i1dxZX3LOLoPbfl1G7+Y8nFuvSEcew9YhBfuW0Of311bd5xysINLzMzMzOzEpm15FVOvX4WQ/v34oqP7o2kvCPVhN71dVxz6n70kPj8z2d3yVvMu+FlZmZmZlYC9y54iX+88RG2G9SHKeccwtABvfOOVFNGDunHNafux3Mr3+DEa/7M/BdW5R2ppNzwMjMzMzPrhIjg1kee4wu/mM244Vtx2+ffyw5b9807Vk06bOxQpnzhEARMnPQQv5u3LO9IJeOGl5mZmZlZB81buopTrpvFxbfP49Bdh3LzZw9icH//XldnjNthK+487zDGDd+Kc29+nEvumMeyVevyjtVp9XkHMDMzMzOrNQ0vr+ZHM57h9ideYEj/XnzzpL345AEjqfdt40ti2MDe3Hr2wVz+u0Xc/PBfue2xpZxy4EjOmbALwwfV5tlEN7zMzMzMKkjSccD3gTrghoi4IudIVoQNmxpZ9NLrTH9qOb+f/xINL6+hV30PvvD+XTj3iF3Yqk/PvCN2Ob3r67jsw+/is+8bzTUzGrjl4ee4+eHn2H/nwUzYfRhH7L4te2w/sGZuYOKGl5mZmVmFSKoDrgGOAZYCj0qaGhFP5Zus+4kINm4ONjU2sm7DZt5I3Zr1G3ll9QZeWbOeFavX8+yra1m0bDVLVqxh4+agh+Cg0dtw+sE788G9tmfbrfrkXZQub8Tgflz+kb05d8Ku3PrIc8xY/ApX3bOYq+5ZzMDe9YzZdgC7DO3PmGH9GTawN0P692ZI/14M6tuTvr3q6Nsz63rWiboeyq2h5oaXmZmZWeUcCDRExBIASb8ETgQ61fD6+KT/peGVNSWIV7siInt8qx8aIyD7ozEiddDYGGxqjFZfq4kEOwzqy+7bD+TIPbOzK4ftOpRtfLfCXIwc0o+vHrcHXz1uD5a//iYPLH6F+S+u4plX1vDQM69y+xMvFPU6dT2yBlidRA9BDwlSW0yAJFTQP2pof+4499BO53fDy8zMzKxydgSeL+hfChxUOIGks4GzAXbaqbgf3z18t6Hsvv3AEkWsXYUHy1l/dgAt0gF2D9EjHWz3rOtBr/oe1PfQW2dF+vWqp3/vOoYO6M22A7OzJv7OVnXabqs+fPyAkXyckW8NW7dhM6+uXc/KtRtYuXYDq9ZtZP3GRtZtzM5mbtrcyMbGYHNjI5s2R9Ygbww2NzXaU1u8eSO+VD8L4IaXmZmZWeW0dI3T3516iYjrgOsAxo8fv+XTMsB5R47tfDKzGte3Vx0jevVjxOB+eUdpkZvwZmZmZpWzFAo+oocRwIs5ZTGzCnLDy8zMzKxyHgXGShotqRdwCjA150xmVgG5NLwkDZE0XdLT6XFwK9PdI+k1SXc3Gz5a0sNp/l+lisvMzPWLmVW1iNgEnAdMAxYCv46IBfmmMrNKyOuM10XA/RExFrg/9bfkauD0FoZfCXw3zf834KyypDSzWuT6xcyqWkT8LiJ2i4hdIuLbeecxs8rIq+F1IjA5PZ8MnNTSRBFxP7C6cJiyG+8fCUzZ0vxm1i25fjEzM7Oqk9ddDbeLiGUAEbFM0rbtmHcb4LV0qh6yL6nu2NrEhbdkBdbfF1PmdyRwlRgKrMg7RCfUen6o/TLUUv6dOzhfReqXLla3tKSW/leK5TLVhkqUqaP1S8XNnj17haS/Fjl5tf0/OE/bqilPNWWB2s7Tav1StoaXpPuA7VsYdUlnX7qFYa3earXwlqySHouI8Z18/9w4f/5qvQy1nr9JNdQvXaluaYnLVBtcpq4vIoYVO221LTvnaVs15ammLNB185St4RURR7c2TtJyScPTp9HDgZfb8dIrgK0l1adPpX0bVrNuxvWLmZmZ1Zq8vuM1FTgjPT8DuLPYGSP7KekZwMSOzG9mXZ7rFzMzM6s6eTW8rgCOkfQ0cEzqR9J4STc0TSTpQeA24ChJSyV9II36V+ACSQ1k38n4SZHve12pCpAT589frZeh1vMXI4/6pSsuV5epNrhMVqjalp3ztK2a8lRTFuiieZR9wGtmZmZmZmblktcZLzMzMzMzs27DDS8zMzMzM7My65INL0k3SnpZ0vyCYZdJekHSk6n7UJ4Z2yJppKQZkhZKWiDpn9PwIZKmS3o6PQ7OO2tr2ihDTawHSX0kPSJpTsr/9TR8tKSH0zr4laReeWdtSRv5b5L0l4Llv0/eWbsKSVdLWiRprqQ7JG2dd6bOkvSx9P/TKKlqbuvbEZKOk7RYUoOki/LO01kt7edqWWv7DHunYo8FJN0j6TVJdzcbXtL9WDvynJGmeVrSGQXDZ6Zts2m/1J7fXmx6jTa3b0m9U1kbUtlHFYy7OA1fXPBd307paB5JoyStK1gWkyqU53BJj0vaJGlis3Etrrcc82wuWD5TK5TnAklPpf37/ZJ2LhjXvuUTEV2uAw4H9gPmFwy7DLgw72xF5h8O7JeeDwT+DxgHXAVclIZfBFyZd9YOlKEm1gPZ7zkNSM97Ag8DBwO/Bk5JwycB5+SdtZ35bwIm5p2vK3bAsUB9en5lNW+f7SjTnsDuwExgfN55OlGOOuAZYAzQC5gDjMs7VyfL9I79XC13re0z8s5VjV2xxwLAUcAJwN3Nhpd0P1ZMHmAIsCQ9Dk7PB6dxnapfitm+gXOBSen5KcCv0vNxafrewOj0OnWdXB6dyTOq1Nt0kXlGAXsDPys8RmhrveWRJ41bk8PyOQLol56fU7C+2r18uuQZr4j4I7Ay7xwdFRHLIuLx9Hw1sBDYETgRmJwmmwyclE/CLWujDDUhMmtSb8/UBXAkMCUNr9p10EZ+K5OIuDey3/4CmEX2G2A1LSIWRsTivHOUwIFAQ0QsiYgNwC/J6tOaVev7ueZqfZ9RYUUdC0TE/cDqwmGSROn3Y8Xk+QAwPSJWRsTfgOnAcZ183ybFbN+FGaeQ3c1WafgvI2J9RPwFaEivl1eecthinoh4NiLmAo3N5i3HeutMnnIoJs+MiHgj9Rbu39u9fLpkw6sN56XThDe2diq82qTTz/uSnbHYLiKWQbaTAtp9Oj4PzcoANbIeJNVJepLsB3ink30i8lrBwfVSqvjAoHn+iGha/t9Oy/+7knrnGLEr+wzw+7xD2Ft2BJ4v6K/qbbe7a2GfYX+vM8cC21D6/Vgxeba0Df40XTr2tQ40QIrZvt+aJpV9FdmyKEfd0Jk8AKMlPSHpAUnv62SWYvOUY95yvWYfSY9JmiWpFB9+tzfPWby9f293WbpTw+taYBdgH2AZ8J/5xtkySQOA3wBfiojX887TES2UoWbWQ0Rsjoh9yD7ZOJDssqt3TFbZVMVrnl/SXsDFwB7AAWSnxv81x4g1R9J9kua30J1YMM0lwCbg5vySFq+YMnUBLR3IVe222511hf1eKZRxu+zQtlCCPG2972kR8W7gfak7vcjXLOa1tzRNOeqGzuRZBuwUEfsCFwC3SNqqAnnKMW+5XnOniBgPnAp8T9Iulcoj6VPAeODq9s7bpL5d0WpYRCxvei7peuDuNibPnaSeZDufmyPi9jR4uaThEbFM0nCyMxlVq6Uy1Np6AIiI1yTNJPuO1NaS6tMnVCOAF3MNV4SC/MdFxHfS4PWSfgpcmF+y2hMRR7c1Pn2x9njgqEgXgFe7LZWpi1gKjCzor4ltt7tpZb/XLbW1XUrqzLHACjqwHytBnqXAhIL+EWTf7SIiXkiPqyXdQvZB58+KK85br72l7btpmqWS6oFBZJfqlqNu6HCetN9YDxARsyU9A+wGPFbmPG3NO6HZvDM7kaWzeYiIF9PjknRssy/ZFUllzSPpaOAS4P0Rsb5g3gnN5p3Z1pt1mzNeqTJocjJQtXeCSqfZfwIsjIj/Khg1FWi6Y8oZwJ2Vzlas1spQK+tB0jClu9JJ6gscTfadgxlA0x12qnYdtJJ/UdPyT+vnJKp0+dciSceRnUH8cMG14FYdHgXGKrubWy+yL7OX5G5YVhpt7PfsnTp8LJAO7Eu9HysmzzTgWEmD01cMjgWmSaqXNBTeangfT/v3S8Vs34UZJwJ/SMtiKnCKsrsMjgbGAo+08/1Lliftu+sAJI1JeZZUIE9rWlxveeVJOXqn50OBQ4Gnyp1H0r7Aj8n274UfLLR/+bR1541a7YBbyU7XbiRrjZ4F/ByYB8xNC3R43jnbyH8Y2anKucCTqfsQ2fW/9wNPp8cheWftQBlqYj2Q3U3niZRzPnBpGj6GrFJuAG4DeuedtZ35/5CW/3zgF6Q7H7oryTJvILvWu+n/fVLemUpQppNTHboeWA5MyztTJ8ryIbI75T0DXJJ3nhKU5x37ubwzdbI8Le4z8s5VjV1rxwJkl0DdUDDdg8ArwLr0P/KBNLyk+7F25PlMes8G4NNpWH9gdlrvC4Dv04G7Cra0fQPfIDtQBuiTytqQyj6mYN5L0nyLgQ+WaB11KA/w0bQc5gCPAydUKM8B6X9kLfAqsKCt9ZZXHuAQsmOYOemxJPVeEXnuI9sHNtVNUzu6fJRmMjMzMzMzszLpNpcampmZmZmZ5cUNLzMzMzMzszJzw8vMzMzMzKzM3PAyMzMzMzMrMze8zMzMzMzMyswNLzMz61YkjZJU9G8FSTpT0g7lzGRmXY+kZ5t+p6wz01jX4YaXmZlZ284E3PAyM7NOccPLciXpa5IWSZou6VZJF0r6nKRHJc2R9BtJ/dK0N0m6VtIMSUskvV/SjZIWSrqp4DXXSLpS0mxJ90k6UNLMNM+H0zSjJD0o6fHUHZLTIjCzfNRLmixprqQpkvpJ2l/SA6numCZpuKSJZD8Ee7OkJyX1lXRpqqPmS7pOkvIujJnlS9JvU92xQNLZzcaNSsc6f1fnFEzyT+lYZJ6kPdI8B0p6SNIT6XH3ihbIysINL8uNpPFkv9K+L/ARsoMbgNsj4oCIeA+wEDirYLbBwJHAl4G7gO8C7wLeLWmfNE1/YGZE7A+sBr4FHAOcTPZL5AAvA8dExH7AJ4D/Lkshzaxa7Q5cFxF7A68DXwR+AExMdceNwLcjYgrwGHBaROwTEeuAH6Y6ai+gL3B8PkUwsyrymVR3jAfOl7RNs/HN65xzC8atSMcj1wIXpmGLgMMjYl/gUuA/ypreKqI+7wDWrR0G3JkOZJB0Vxq+l6RvAVsDA4BpBfPcFREhaR6wPCLmpXkXAKOAJ4ENwD1p+nnA+ojYmOYZlYb3BH6YGmubgd3KU0Qzq1LPR8Sf0/NfAP8G7AVMTyew6oBlrcx7hKSvAv2AIcACsg+CzKz7Ol/Syen5SGBss/HN65zzge+k/tvT42yyD6IBBgGTJY0Fguy4xWqcG16Wp9Yuz7kJOCki5kg6E5hQMG59emwseN7U3/T/vDEiovl0EdEoqWmaLwPLgfeQnfl9s8OlMLNaFM36VwMLIuK9bc0kqQ/wI2B8RDwv6TKgT3kimlktkDQBOBp4b0S8IWkm76wXmtc5hf1NxzObeftY5pvAjIg4WdIoYGbpEltefKmh5elPwAmS+kgaAPxDGj4QWCapJ3Bamd57ELAsIhqB08k+3Taz7mMnSU2NrE8Cs4BhTcMk9ZT0rjR+NVm9BG8fTK1I9dbESgU2s6o1CPhbanTtARzcwjTN65w/FfGaL6TnZ5YkpeXODS/LTUQ8CkwF5pCdZn8MWAV8DXgYmE52jXM5/Ag4Q9IssssM15bpfcysOi0kqwPmkl0u+AOyRtSVkuaQXbbcdNOdm4BJkp4k+2T6erLLmH8LPFrh3GZWfe4hu2HPXLIzVbNamKZ5nXPtFl7zKuBySX/GHw53GXr7iiyzypM0ICLWpLv7/BE4OyIezzuXmZmZWSmkSwXvTjfksW7M3/GyvF0naRzZ5TuT3egyMzMzs67IZ7zMzMzMzMzKzN/xMjMzMzMzKzM3vMzMzMzMzMrMDS8zMzMzM7Myc8PLzMzMzMyszNzwMjMzMzMzK7P/B62tFz4gk3fqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot3(np.exp(log_Ppdf-norm),gammas,betas,alphas,zt='alpha',yt='beta',xt='gamma')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the fit from the maximum a posteriori (MAP) parameters. We'll get the index of the maximum using numpy.argmax(), then use numpy.unravel_index() to convert this 1D index into 2D indices in the posterior array. Using those indices, we'll get the MAP values of the parameters and plot them with the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.09500000000000017 0.280000000000002 27.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11f7ce150>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU5fnG8e+TECBsBmQRoghapa4solWx1oKKOxFRqRtVK7UurV3SShdr/dlCS/e6UpfS1gVFDGrdKGjdEI0GARdARJEEQwTCGkkyeX9/nBNNwiRMJnPmTGbuz3XlysyZM5mHIdy88573PMecc4iISObICrsAERFJLgW/iEiGUfCLiGQYBb+ISIZR8IuIZJgOYRcQi969e7tBgwaFXYaISLvyxhtvfOqc69N0e7sI/kGDBlFcXBx2GSIi7YqZfRRtu6Z6REQyjIJfRCTDBBb8ZjbEzBY3+NpiZteZWS8zm2dmK/3vPYOqQUREdhVY8DvnljvnhjnnhgFHADuAR4HrgfnOuQOA+f59ERFJkmRN9YwBVjnnPgLGATP97TOBgiTVICIiJC/4JwIP+Lf7OefWAfjf+yapBhERIQnBb2YdgbOAh1v5vMlmVmxmxRUVFcEUJyKSgZIx4j8VeNM5V+7fLzez/gD+9/XRnuScm+GcG+mcG9mnzy7nH4iISJySEfzf4ItpHoDHgEn+7UnA3CTUICLSrpx/50LOv3NhID870OA3sy7AScCcBpunASeZ2Ur/sWlB1iAiIo0F2rLBObcD2LPJtg14q3xERCQEOnNXRCTDKPhFRDKMgl9EJMMo+EVEMoyCX0Qkwyj4RUQyjIJfRCTDKPhFRFJMUUkpJWsqWbR6I6OmLaCopDShP1/BLyKSQopKSpkyZynVkToASiurmDJnaULDX8EvIhKAeHvtTH9mOVU1kUbbqmoiTH9meaJKU/CLiEQTZJO0lpRV7mhme1XCXkPBLyKSKl69gwF8GvWhAXm5CXsZBb+ISCpYeBs8/RMKB64gN6dxNOfmZFM4dkjCXkrBLyIStoW3wjNT4KAzKfj2r5g6/nA6ZnvxnJ+Xy9Txh1EwPD9hLxdoW2YREdmNV26BZ38GB50FE+6B7BwKhufzwGtrAJj17WMS/pIa8YuIhOXlv3qhf3DB56GfDAp+EZEwvPRnmPcLOGQ8nHN30kIfNNUjIpJ8L/4R5v8KDj0Hzp4B2cmNYo34RUQSrMWWCy/83gv9w84NJfRBI34RkV3UB3d1pI5R0xZQOHZIzKtqmmu5AFCw5X547mY4/HwouB2ysgP7M7REI34RkQba2iun2ZYLc1/zQ39iqKEPCn4RkUba2iunudYKZZ91gKEXQMFtoYY+aKpHRKSRZoM7xl45A/JyKY2y74COn8G4W2IO/SDW79fTiF9EpIHmeuLE2iuncOwQcnMah3tuVi2F444OfaRfT8EvItJA1OBuRa+cguH5TD37UHpmbceoI79TFVPPGUHBEfsEUW5cNNUjItJA/eqdH89eQnWkjvy83Fat6qEuQkHp7ynoeC/PdjmDk3/0L8hKrTG2gl9EpIm4e+VEaqDoO7D0YYq6nscD3S/l5BQLfVDwi4gkRs1nMPtSWP4kjPklD7wT3MHZtkq9/4pERNqbndvg/nO90D/t9/DVH4RdUYs04hcRaYuqTXDfuVD6Jpx9JwydGHZFu6XgFxGJ17b18K+z4dMVcN5MOOjMsCuKiYJfRCQelR/DvwpgSxlcMAv2Hx12RTFT8IuItNaGVTDzLNi5FS5+FAYeHXZFrRLowV0zyzOz2Wb2npm9a2bHmFkvM5tnZiv97z2DrEFEJKE+WQb3nAK1VfDNx9td6EPwI/6/AE875yaYWUegC/BTYL5zbpqZXQ9cD/wk4DpERNpubTH8ezzkdIVL/gN9Dmx21yB77bSVOeeC+cFmPYC3gP1cgxcxs+XACc65dWbWH3jeOdfiudAjR450xcXFgdQpIhKT1S/A/ROhW1+4ZC703DfsinbLzN5wzo1suj3IqZ79gArgXjMrMbO7zKwr0M85tw7A/963mYInm1mxmRVXVFQEWKaIyG4sfwr+PQHyBsJlT7eL0G9JkMHfARgB3O6cGw5sx5vWiYlzboZzbqRzbmSfPn2CqlFEpGVLZ8Osi6DfIXDpk9B9r7ArarMgg38tsNY5t8i/PxvvP4Jyf4oH//v6AGsQEYlf8b3wyLdgn6O96Z0uvcKuKCECC37n3CfAx2ZWP38/BngHeAyY5G+bBMwNqgYRkbg4By/+EZ64Dg44GS6aDZ17hF1VwgS9quda4D5/Rc8HwKV4/9k8ZGaXA2uAcwOuQUQkdnUReOon8Prf4dAJ3vVxO3QMu6qECjT4nXOLgV2OKOON/kVEUkvNZzDnCnj3MTj2WjjxppTrpZ8IOnNXRAS8ZmsPXABrFsLYqXDMVWFXFBgFv4ikpfPvXAjEeCLV5rXw73Ng4wcw4R44dHzA1YVLwS8ima38bW+NfvU2uOgRGHx82BUFTsEvIplr9Yvw4IXQsYt3Yla/Q8KuKCnS76iFiEgs3n7U67vTfS+4fF7GhD4o+EUkE716Bzx8KQwY4Y308/YJu6Kk0lSPiGSOujr47y/hlb96V8sa/3fIyQ27qqRT8ItIZqithrlXw9KH4Mhvwam/g6zssKsKhYJfRNLfZ1vgoYvhg+dhzA1w3A/ALOyqQqPgF5H0trUc7jsH1r/rtV8YdkHYFYVOwS8iKalVJ2A1UVRSSsmaSqojEUZNXUVhTm8KLpoFB5yY6DLbJa3qEZG0UlRSypQ5S6mO1AFGaV1PpkSuoGjbQWGXljIU/CKSVqY/s5yqmkijbVW13nbxKPhFJH3U1VFWuSPqQ2WVVUkuJnUp+EUkPdRUwSOXMYBPoz48IC/z1us3R8EvIu3f1k/gH6fD20UUjqgjN6dxtOXmZFM4dkgzT848WtUjIu3bJ0vh/oleP/2J91Pw5dPggFJ+PHsJ1ZE68vNyKRw7hILh+WFXmjIU/CLSfi1/CmZfDrl5Xs+d/ocDUDA8nwdeWwPEtxw03WmqR0TaH+fglb/BA9+APgfCFQs+D33ZPQW/iKSc+hOwFq3eyKhpCygqKf3iwdpqePy78OzP4eCz4JtPeq2VJWaa6hGRlNL4BCworaxiypylABQMyYWHLoEPX4TjC+GEn6blxdCDpuAXkZQS9QSsmgjTn1pGwYtTYPPHcPadMHRiSBW2fwp+EUkpzZ1oVbalBjpUwqTHYeDRSa4qvegzkoiklOZOtBqQvRm+NV+hnwAKfhFJKYVjh5Cb0/gCKblWQ+G4o6DX4JCqSi8KfhFJKQXD85l65n70sS0YdeR3+oypE0ZQcNQBYZeWNjTHLyKppfxtChZeyBmdPubeHt/hih/cHNeP0YlbzdOIX0RSx7JH4K4ToaaKG/f8Hf/tenrYFaUlBb+IhC9S652QNfsy2Otw+Pb/WNnx4LCrSlua6hGRcG3fALMvhdX/gyOvgLG/gQ4dgdVhV5a2FPwiEp6yxTDrYthWDuNug+EXhl1RRlDwi0g4Fj8AT1wHXXp7nTXzR4RdUcZQ8ItIckVq4JmfwmszYNBX4dx/QNfeYVeVUQINfjP7ENgKRIBa59xIM+sFzAIGAR8C5znnNgVZh4ikiK3l8PA3Yc0rcMw1cOKvIFvjz2RLxqqerzvnhjnnRvr3rwfmO+cOAOb790Uk3X38Osz4GpSVwDl3w9hfK/RDEsa7Pg44wb89E3ge+EkIdYhIshTfC08Wwh758K3/wl6H7vYpOgErOEGP+B3wrJm9YWaT/W39nHPrAPzvfaM90cwmm1mxmRVXVFQEXKaIBKJ2Jzx2rXcQd/DxcMVzMYW+BCvoEf8o51yZmfUF5pnZe7E+0Tk3A5gBMHLkSBdUgSISkM2l3kVTSovhqz+Er/8MsrJ3/zwJXKDB75wr87+vN7NHgaOAcjPr75xbZ2b9gfVB1iAiIVjxDDx6JUSq4bx/eZdIlJQR2FSPmXU1s+71t4GTgWXAY8Akf7dJwNygahCRJIvUeK0X7j+P1TU9uW6PPyv0U1CQI/5+wKNmVv869zvnnjaz14GHzOxyYA1wboA1iGSs8+9cCCTxIOmmj7xeO6XFcOS3+MXHZ1FjHZPz2tIqgQW/c+4DYGiU7RuAMUG9roiE4N0nYO5V4BycOxMOKaDG/49HUo8W0YpI/Gp3wrwbYNEdMGA4TLhXV8lqBxT8IhKfjR/Aw5fCusVw9FXeWbgdNLXTHij4RaT1ls2Bx77rLc+c+AB8+bSwK5JWUPCLSOxqqrwGa8X3wN5HwoR7IG9g2FVJKyn4RSQ2n670GqyVL4NR34PRv4DsnLCrkjgo+EVk996aBU98Hzp0ggsehgNPDrsiaQNdc1dEmle9HeZeDY9Ohv5D4cqXYgr9opJSStZUsmj1RkZNW0BRSWkSipVYacQvItGVv+NdC7diORxfCF+7PqY2ykUlpUyZs5TqSB0ApZVVTJmzFICC4fmBliyx0YhfJA21acRdVwev/M3rnb9jA1w8B0b/PObe+dOfWU5VTaTRtqqaCNOfWd6aP4IESCN+kTTTphF35Roougo+fBG+fAac+ZdWXxaxrLKqVdsl+TTiF0kzcY24nfMufn77KChbDONuhfP/Hde1cAfk5bZquyTfboPfzK4xs57JKEZE2q7VI+7tG7y++UVXQr9D4DsvwfCLwGuw2GqFY4eQm9O4735uTjaFY4fE9fMk8WKZ6tkLeN3M3gTuAZ5xzunCKCIpakBeLqVRQj7qiHvlPG/Vzo6NXsuFY69t88VS6qeTfjx7CdWROvLzcikcO0QHdlPIbkf8zrmfAwcAdwPfBFaa2W/MbP+AaxOROMQ04q7e7q3Lv28CdNkTJj8Hx13XKPTPv3Ph562dW6tgeD7DB+bxlcG9ePn60Qr9FBPTwV3nnDOzT4BPgFqgJzDbzOY5534cZIEi0jq7HXGvLYY5k70ma8deC1//OeR0DrFiSbbdBr+ZfRfvSlmfAncBhc65GjPLAlYCCn6RFFMwPJ8HXlsDNLgQS6QG/vc7ePEP0GMATHocBn81xColLLGM+HsD451zHzXc6JyrM7MzgilLRBKqYgXMucJroTz0G3Dqb6HzHmFXJSHZbfA7525o4bF3E1uOiCSSuTpYNAPm/QJyusB5/4SDx4VdloRMJ3CJpKlekQqurPwTPPUmfOkkGHcLdN8r7LIkBSj4RVJU3BdLr6uDN+7hDxW/IIsInP5HGHlZ3OvyJf3ozF2RdPLpSvjH6fCfH7Iq50AKe98BR16u0JdGNOIXSQeRGnj5L96qnZzOMO5Wbl60nwJfolLwi7R3ZSUw91ooX+oduD11OnTvB6/Fd/IVfNHdszpSx6hpC+I687bVU1SSNAp+kfaqegc8PxUW3gJd+3pN1Q468/OH4w1e9dNPf5rjF2mPVr8Atx8Lr/zVa6h29aJGod8W6qef/jTiF2lPqiq9Nflv/hN6DvbPvj0+oS+hfvrpT8Ev0l68+zj850ewfT0c+104YQp07JLwl2lVd09plzTVI9KMtnSnTKit5TDrYph1EXTtA1csgJP/L5DQB/XTzwQa8YukoC9W1UQY9ZuVFOZspuCkG7yRfnZOoK+tfvrpT8EvkmKKSkqZ8shbVEccYJS6PZlSdyV0G0pBwKFfL2p3T0kbmuoRCUhcU0VVlUwvepWq2sYXuauqdVpVIwmj4BdJBXV1UHIf3DKSsp2dou6iVTWSKJrqEQlb2WJ4shDWvgZ7H8mA6g6Ubo3ssptW1UiiBD7iN7NsMysxsyf8+4PNbJGZrTSzWWbWMegaRFLSjo3edW9nnACbVsO42+CyZyk87TCtqpFAJWOq53tAwwu2/Bb4k3PuAGATcHkSahBJHXURKL4X/nYEvDETvnIlXFMMwy+ErCwKhuczdfxhdMz2/nnm5+UydfxhWlUjCRPoVI+Z7Q2cDvwa+IGZGTAauMDfZSZwI3B7kHWIpIy1xfDkj7zGavuOglN/B3sdustuWlUjQQp6jv/PeBdj7+7f3xOodM7V+vfXAlGHMWY2GZgMMHDgwIDLFAnYtgqYfyOU/Bu67QXj74LDJqhtsoQisKke/0Ls651zbzTcHGVXF2UbzrkZzrmRzrmRffr0CaRGkebUn0C1aPVGRk1bQFFJaVw/J8tFvGve3nIEvPWgdwLWtcVw+LkKfQlNkCP+UcBZZnYa0BnogfcJIM/MOvij/r2BsgBrEGm1RLUlPnjnW3xzyx3w1GrY7wSvT36fAwOoOBiaYkpfgY34nXNTnHN7O+cGAROBBc65C4HngAn+bpOAuUHVIOFKmV43rdTmtsTl71B0y495tqwzJ2y7mVHZ91F06K3tKvQlvYWxjv8nwINmdjNQAtwdQg0izYq7LfHmUnj+NxQVr2ZKzbeowjsRq3Q7THl0GZhpZY6khKQEv3PueeB5//YHwFHJeF2ReLS6LfFnm+GlP8Ort4GrY3rWHZ+Hfr36TwwKfkkFatkgKSusqaKY2xLXVsOrt8NfhsFLf4SDzoJriinb2Tnqz1XLBUkVatkg0sRu2xLX1cE7j8L8m2DThzD4a3DSTTBgGAAD8lbpQiaS0hT8IlE0ewLV6hdg3g3eCVj9DoWLHoH9xzRamlk4dghT5ixtdIA4npYLWlUjQVHwi8Si/B347y9h5bPQY28ouB0OPx+ysnfZVRcykVSn4BdpQa9IBcy9GhbfDx27w4m/gq98G3JanrZRywVJZQp+kWi2refCLX/nlO2Pwwbg6Kvgqz+ELr3CrkykzRT8Ig1tLYeX/wLF93BG7U5eyj2B4yf/GXruG3ZlIgmj4JdAfHGx8DpGTVsQyhx3/VLQmKZatn7ircV/416IVMPh5/ODsjGs67A3xyv0Jc0o+CXhEtXrJim2lPmB/w+oq4WhE70pnT33Z107bDchEgsFvyRcS71uYg3+wD8xbF7rBf6bM8HVfRH4vfZL3GuIpCgFvyRc3L1ufIF+Yqj82DvLtuTfXuAPu8AL/J6DdtlVq3EkXSn4JeFa3eumiUR8YtjFpo/8wL/Puz/8Ijju+zpoKxlJwS8J19YzV9v6iaGhvrXr4LFrvXX4lgUjLvECP2+fVv8skXSh4JeEa+uZq239xADAuiVcXTmdY6ueh43ZcMSlcNx1sMfesf+MNtJUkaQqBb8Eoi1nrsb9iaGuzmupsPAW+PBFjrLOPNvlTE67cir0GNDqP4NIulLwp7BWrUNPI63+xFC9A5Y8CAtvgw0roUc+nHQT31lyMDuyunGaQl+kEQW/pKSYPjFsLYfX7/K+qjZC/6Ew/i44pICiJet5ee0SqiMbQzuBTCRVKfil/Sl/2xvdL30IIjUw5DQ45mrY91gwa18nkImEQMEv7YNzsGo+vHILfPAcdMj1VugcfRXsuX+jXQNZDiqSRhT8ktJyXDW8+U9YeCtUvAfd9oIxN3irdJrplJnI5aAi6UjBL6lp00ecv/UfjN7xNDxWCf0Og7PvhEPGQ4eOLT41IctBRdKYgl9SR6QGlj/pNUxb9RwFGCWdjuSI838Gg49vdHnDliTq0oci6UrBL4GJeRnqhlXedM7i+2B7hbcc84Trueadg9mQ3YdZ+7VuOasufSjSMgW/hKN2J7z7uNcdc/ULYNlw4ClwxDfhS2MgK5sNy+Nvi6xLH4o0T8EvyVWxwgv7xfd7a+/zBsLon8Owi6BH/7CrE8kICv4UlQpXsEqYmip4Zy68MRPWvAJZHeDLp3uj+8EnQFZW1KdppC4SDAV/CkqLE5Ccg3VveSP7JQ/CZ5u9i5ycdBMMvQC69Qm7QpGMpeBPQalyAlJcvYIqVsCyR2DZbNjwPmR3hIPO8kb3g46LeWWOiARHwZ+C2t0JSJVrYNkcL+w/WQqYF/LHXuuFfjMnWolIOBT8KahdnIC0bT28XeSF/ceLvG35I+GUaXBwgQ7UiqQwBX8KStkTkKo2wbtPeGG/+gXvmrV9D/FaKBwyHnoNDrc+EYmJgj9A8fbTT6UTkDrVfQZLZ3vz9u//FyLV0HOwd4HyQ8+BvgclvSYRaRsFf4pKxAlIcV/IZVsFvD+P72+6j2E7X4dHdkL3/nDkFXDYOTBgRLs4SKvloCLRBRb8ZtYZeAHo5L/ObOfcL81sMPAg0At4E7jYOVcdVB0SA+egfBmseBpWPANriwHHgVl78kLuiZx83lUw8BjIyg67UhFJgCBH/DuB0c65bWaWA7xkZk8BPwD+5Jx70MzuAC4Hbg+wDommeoc3T7/iae86tVtKve35R8DXfwoHjuU7c7eDGScP0shZJJ0EFvzOOQds8+/m+F8OGA1c4G+fCdyIgj85Nq/1RvQrnoHV/4Paz6BjN9j/617Yf+kk6N7vi/0t/l45IpK6Ap3jN7Ns4A3gS8CtwCqg0jlX6++yFoh6xNLMJgOTAQYOHBhkmWnLXIT9a1bCgvmw/Gko987+JW9f74SqA8fCvqOgQ6ddnptWLSNEpJFAg985FwGGmVke8CgQbQmIa+a5M4AZACNHjoy6jzRRF4FPlsCHL8GHL3NX+Yt0c9vgxSxvjv6km7wOmL0PbPHgbFq0jBCRZiVlVY9zrtLMngeOBvLMrIM/6t8bKAvqdeNe1dJeRGq9fjgfveSF/ZpXYecW77Fe+7Go83G83Wko3518ZavOnk2VlhEiEowgV/X0AWr80M8FTgR+CzwHTMBb2TMJmBtUDWknUgNlJV7If/SyF/TV/mGU3gd66+oHHedN3/Tozwz/P77vtrJlQrtrGSEirRLkiL8/MNOf588CHnLOPWFm7wAPmtnNQAlwd4A1hCYRc+SzLh/hBf0L0+HDl+Hj16Bmu/dgn4Ng6EQv5Pcd1figbBu1i5YRIhK3IFf1LAGGR9n+AXBUUK+bCuKaI6+pgvK3vaBft9ibwln/LtT5x8H7HQrDL4JBftB37R1Y/SnbMkJEEkJn7gZgt3Pk1du9Lpbr3oIyP+Qr3gPnP6fLntB/GIw62TtLdt9jW93hsi2fOFKpZYSIJJ6CPwDNz5HvgFuOgk9X8Plipq59YcAw74pU/Yd6t3vkt6klQiJW5eiatSLpS8HfgphXBdVFvJ70G1fBhlUM6NSb0p2dd9ltQNYm7ypUh473Qr7/sEDaF2tVjoi0RMEfK+dgW7l3VakN78OGVf7X+7Bptde10leYNYYpdglVLufzbbk5WRSOHw3DLw68VK3KEZGWKPjrVW+HrZ944e5//8aWN+kbKYc7K72Qr186CZDdyRu99z4AhpwCe37p86+Crn1gcVloc+RalSMiLUnv4HeOXLfDuw7stk9gazlsXdco3D//Xn/iUwNn0IEN2b2hyyGwz9F+sO/vfd9j7xa7VYY5R65VOSLSkrQN/qKSUj746CM+dd0Y9YdXKOwwi4IOr3gPdsj11r132wv6HQz7j/7ifvd+Xu/5bntx0b/ew1kWsy5uXwc3tSpHRFqSlsFfv6qlyvUAoJQ+TOEaGP1/FBx1AHTqEdOqGWdZQZcaGK3KEZHmtN9ka0HUVS21MP2VLdB5j3Zx9SgRkaCkZfBrVYuISPPScqpHq1oSQ1NEIukpLUf8hWOHkJvTeMVNa1e11Lc8WLR6I6OmLaCopDTRZYqIhCItR/xtXdWiC5GISDpLy+CHtq1qUcsDEUlnaRv8bZGog8OaIxeRVJSWc/xt1dxBYB0cFpF0oOCPIhEHh1PBrG8fo08dIrILTfVEoZYHIpLOFPzNUMsDEUlXmuoREckwCn4RkQyT1lM9mqIREdmVRvwiIhlGwS8ikmEU/CIiGUbBLyKSYRT8IiIZJq1X9bSVVgWJSDrSiF9EJMMo+EVEMoyCX0Qkwyj4RUQyjIJfRCTDBBb8ZraPmT1nZu+a2dtm9j1/ey8zm2dmK/3vPYOqQUREdhXkiL8W+KFz7iDgaOBqMzsYuB6Y75w7AJjv3xcRkSQJLPidc+ucc2/6t7cC7wL5wDhgpr/bTKAgqBpERGRXSZnjN7NBwHBgEdDPObcOvP8cgL7NPGeymRWbWXFFRUUyyhQRyQjmnAv2Bcy6Af8Dfu2cm2Nmlc65vAaPb3LOtTjPb2YVwEdxltAb+DTO5yaD6msb1dc2qq9tUr2+fZ1zfZpuDLRlg5nlAI8A9znn5viby82sv3NunZn1B9bv7udEK7wVNRQ750bG+/ygqb62UX1to/raJtXra06Qq3oMuBt41zn3xwYPPQZM8m9PAuYGVYOIiOwqyBH/KOBiYKmZLfa3/RSYBjxkZpcDa4BzA6xBRESaCCz4nXMvAdbMw2OCet0oZiTxteKh+tpG9bWN6mubVK8vqsAP7oqISGpRywYRkQyj4BcRyTBpE/xmdoqZLTez981slzYQZtbJzGb5jy/yTypLVm1R+xY12ecEM9tsZov9rxuSVZ//+h+a2VL/tYujPG5m9lf//VtiZiOSWNuQBu/LYjPbYmbXNdknqe+fmd1jZuvNbFmDbTH1oTKzSf4+K81sUrR9Aqpvupm95//9PWpmec08t8XfhQDru9HMShv8HZ7WzHNb/LceYH2zGtT2YYNFK02fG/j712bOuXb/BWQDq4D9gI7AW8DBTfa5CrjDvz0RmJXE+voDI/zb3YEVUeo7AXgixPfwQ6B3C4+fBjyFd8D+aGBRiH/Xn+CdmBLa+wccD4wAljXY9jvgev/29cBvozyvF/CB/72nf7tnkuo7Gejg3/5ttPpi+V0IsL4bgR/F8Pff4r/1oOpr8vgfgBvCev/a+pUuI/6jgPedcx8456qBB/F6AjXUsEfQbGCMf65B4FzzfYvak3HAP53nVSDPPwEv2cYAq5xz8Z7JnRDOuReAjU02x9KHaiwwzzm30Tm3CZgHnJKM+pxzzzrnav27rwJ7J/p1Y9XM+xeLWP6tt1lL9fm5cR7wQKJfN1nSJfjzgY8b3F/LrsH6+T7+L/9mYM+kVNdAk75FTR1jZm+Z2VNmdkhSCwMHPGtmb5jZ5CiPx/IeJ8NEmv8HF+b7B7H1oUqV9/EyvE9w0ezudyFI1/hTUfc0M1WWCu/fV4Fy59zKZh4P8/2LSboEf7SRe9N1qrHsEyjz+hY9AlznnNvS5OE38aYvhgJ/A4qSWRswyjk3AjgVr4X28U0eT4X3r5vgH+QAAALhSURBVCNwFvBwlIfDfv9ilQrv48/w2qbf18wuu/tdCMrtwP7AMGAd3nRKU6G/f8A3aHm0H9b7F7N0Cf61wD4N7u8NlDW3j5l1APYgvo+acbHofYs+55zb4pzb5t9+Esgxs97Jqs85V+Z/Xw88iveRuqFY3uOgnQq86Zwrb/pA2O+fr7x++sua70MV6vvoH0w+A7jQ+RPSTcXwuxAI51y5cy7inKsD/t7M64b9/nUAxgOzmtsnrPevNdIl+F8HDjCzwf6ocCJeT6CGGvYImgAsaO4XP9H8OcFofYsa7rNX/TEHMzsK7+9mQ5Lq62pm3etv4x0EXNZkt8eAS/zVPUcDm+unNZKo2ZFWmO9fA7H0oXoGONnMevpTGSf72wJnZqcAPwHOcs7taGafWH4Xgqqv4TGjs5t53Vj+rQfpROA959zaaA+G+f61SthHlxP1hbfqZAXeEf+f+dtuwvslB+iMN0XwPvAasF8SazsO7+PoEmCx/3UacCVwpb/PNcDbeKsUXgWOTWJ9+/mv+5ZfQ/3717A+A27139+lwMgk//12wQvyPRpsC+39w/sPaB1QgzcKvRzvmNF8YKX/vZe/70jgrgbPvcz/PXwfuDSJ9b2PNz9e/ztYv8ptAPBkS78LSarvX/7v1hK8MO/ftD7//i7/1pNRn7/9H/W/cw32Tfr719YvtWwQEckw6TLVIyIiMVLwi4hkGAW/iEiGUfCLiGQYBb+ISIZR8IuIZBgFv4hIhlHwi8TBzI70m4l19s/WfNvMDg27LpFY6AQukTiZ2c14Z4TnAmudc1NDLkkkJgp+kTj5vWJeBz7DaxERCbkkkZhoqkckfr2AbnhXVescci0iMdOIXyROZvYY3hWgBuM1FLsm5JJEYtIh7AJE2iMzuwSodc7db2bZwCtmNto5tyDs2kR2RyN+EZEMozl+EZEMo+AXEckwCn4RkQyj4BcRyTAKfhGRDKPgFxHJMAp+EZEM8/90GTMdXBNQSAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "i=np.unravel_index(np.argmax(log_Ppdf),log_Ppdf.shape)\n",
    "best_alpha = alphas[i[0]]\n",
    "best_beta = betas[i[1]]\n",
    "best_gamma = gammas[i[2]]\n",
    "print(best_alpha,best_beta,best_gamma)\n",
    "plt.errorbar(xdata, ydata, yerr=sig, marker='o',ls='none')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.plot(xdata,np.polyval([best_alpha,best_beta,best_gamma],xdata))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, now we will integrate over the posterior. We'll use the <b> same </b> normalizing factor as last time, so we can compare the models. Again, we have to incorporate the bin sizes in the integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evidence from quadratic model:  2.3918900606334224\n"
     ]
    }
   ],
   "source": [
    "total_Ppdf = np.exp(log_Ppdf-norm).sum()*dalpha*dbeta*dgamma\n",
    "evidence3 = total_Ppdf\n",
    "print('Evidence from quadratic model: ', total_Ppdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What do you think about the quality of your fit? From the Bayesian perspective, can you \"accept\" or \"reject\" it? Why or why not?\n",
    "<br>\n",
    "<B> ANSWER HERE: </B>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Bayes factor</h3>\n",
    "\n",
    "OK, now we can calculate the Bayes factor, which is the ratio of the evidences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bayes factor (linear/quadratic): 0.20927705265735574\n"
     ]
    }
   ],
   "source": [
    "bf=evidence2/evidence3\n",
    "print('Bayes factor (linear/quadratic):', bf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally, the Bayes factor needs to be >10 or <0.1 (depending on which is the numerator and denominator) to strongly prefer one model over another. Is either model strongly preferred? Why or why not?\n",
    "<br>\n",
    "<B> ANSWER HERE: </B>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now go back and reduce and/or increase your uncertainties, e.g., to 1 and to 10 and discuss how this affects the Bayes factor. Remember that you may need to adjust your trial grids.\n",
    "\n",
    "<br>\n",
    "<B> ANSWER HERE: </B>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}