Last active
July 7, 2020 15:06
-
-
Save simeoncarstens/ab1a0bc6f00a4403783b0bfc860573d3 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| from numpy import sin, cos, arctan, sqrt, exp, random, pi, linspace | |
| import matplotlib.pyplot as plt | |
| def draw_sample(xold, sigma): | |
| t = 3.0 | |
| vold = random.normal() | |
| phi = arctan(-vold / xold * sigma) | |
| A = vold * sigma * sqrt(xold ** 2 / sigma ** 2 / vold ** 2 + 1) | |
| xnew = A * cos(t / sigma + phi) | |
| vnew = -A / sigma * sin(t / sigma + phi) | |
| E = lambda x: 0.5 * x ** 2 / sigma ** 2 | |
| K = lambda v: 0.5 * v ** 2 | |
| H = lambda x, v: E(x) + K(v) | |
| p_acc = min(1, exp(-(H(xnew, vnew) - H(xold, vold)))) | |
| if random.random() < p_acc: | |
| return xnew, True | |
| else: | |
| return xold, False | |
| sigma = 2.0 | |
| samples = [2.0] | |
| accepted = 0 | |
| n_samples = 100000 | |
| for _ in range(n_samples): | |
| new_state, acc = draw_sample(samples[-1], sigma) | |
| samples.append(new_state) | |
| accepted += acc | |
| fig, ax = plt.subplots() | |
| ax.hist(samples, bins=40, density=True) | |
| gaussian = lambda x: exp(-0.5 * x ** 2 / sigma ** 2) / sqrt(2 * pi * sigma ** 2) | |
| xspace = linspace(-5, 5, 300) | |
| ax.plot(xspace, list(map(gaussian, xspace))) | |
| plt.show() | |
| print("Acceptante rate:", accepted / n_samples) |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/bin/bash | |
| img_list=$(ls -v output*.png) | |
| b=$(<$2) | |
| while read strA <&3 && read strB <&4; do | |
| rstring="..\/..\/img\/posts\/${strB}" | |
| echo $rstring | |
| sed -i "s/${strA}/${rstring}/g" $1 | |
| mv $strA $strB | |
| # cp $strB ~/projects/tweag/www/app/assets/img/posts/ | |
| done 3<<<"$img_list" 4<<<"$b" | |
| # cp $1 ~/projects/tweag/www/app/views/posts/ |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "# Markov Chain Monte Carlo (MCMC) sampling, part I: the basics" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Markov Chain Monte Carlo (MCMC) is a powerful class of methods to sample from probability distributions known only up to an (unknown) normalization constant.\n", | |
| "But before we dive into MCMC, let's consider why you might want to do sampling in the first place. \n", | |
| "The answer to that is: whenever you need to fully characterize a distribution because you're either interested in the samples themselves or because you need them to approximate expectations of functions w.r.t. to a probability distribution.\n", | |
| "Sometimes, only the mode of a probability distribution is of primary interest. In this case, it can be obtained by numerical optimization and full sampling is not necessary. \n", | |
| "It turns out that sampling from any but the most basic probability distributions is a difficult task, because for those, you usually don't know the normalization constant.\n", | |
| "You then can't easily get random numbers distributed according to this distribution: [inverse transform sampling](https://en.wikipedia.org/wiki/Inverse_transform_sampling) requires proper, normalized distributions.\n", | |
| "Now in principle you could just obtain the normalization constant by numerical integration, but this quickly gets infeasible with increasing number of dimensions.\n", | |
| "[Rejection sampling](https://en.wikipedia.org/wiki/Rejection_sampling) does not require normalized distributions, but implementing it efficiently can be very hard and requires a good deal of knowledge about the distribution of interest.\n", | |
| "That's when you need a smart way to obtain representative samples from your distribution which doesn't require knowledge of the normalization constant. \n", | |
| "MCMC algorithms are a class of methods which do exactly that.\n", | |
| "These methods date back to a seminal paper by Metropolis et al., who developed the first MCMC algorithm, correspondingly called [Metropolis algorithm](https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm), to calculate the equation of state of a two-dimensional system of hard spheres, but really were looking for a general method to calculate expectation values occurring in statistical physics. \n", | |
| "\n", | |
| "In this series of blog posts, I will introduce several important, increasingly complex MCMC algorithms.\n", | |
| "Along the way, I hope you will acquire a solid understanding of challenges you may face when sampling from non-standard probability distributions and how to address them with the methods presented here. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Markov chains" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Now that we know why we want to sample, let's get to the heart of MCMC — Markov chains.\n", | |
| "What is a Markov chain?\n", | |
| "Really informally and without all the technical details, a Markov chain is a random sequence of states in some state space in which the probability of picking a certain state next depends only on the current state in the chain and not on the previous history: it is memory-less.\n", | |
| "Under certain conditions, a Markov chain has a unique stationary distribution of states to which it will converge after a certain number of states.\n", | |
| "From that number on, states in the Markov chain will be distributed according to the invariant distribution.\n", | |
| "MCMC algorithms work by constructing a Markov chain with the probability distribution you want to sample from as the stationary distribution." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "For didactic purposes, let's for now consider both a discrete state space and discrete \"time\".\n", | |
| "The key quantity characterizing a Markov chain is the transition operator $T(x_{i+1}|x_i)$ which gives you the probability of being in state $x_{i+1}$ at time $i+1$ given that the chain is in state $x_i$ at time $i$." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Now just for fun (and for illustration), let's quickly whip up a Markov chain which has a unique stationary distribution.\n", | |
| "We'll start with some imports and settings for the plots:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "%matplotlib notebook\n", | |
| "%matplotlib inline\n", | |
| "import numpy as np\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "plt.rcParams['figure.figsize'] = [10, 6]\n", | |
| "np.random.seed(42)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "The Markov chain will hop around on a discrete state space which is made up from three weather states:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "state_space = (\"sunny\", \"cloudy\", \"rainy\")" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "In a discrete state space, the transition operator is just a matrix.\n", | |
| "Columns and rows correspond, in our case, to sunny, cloudy, and rainy weather.\n", | |
| "We pick more or less sensible values for all transition probabilities:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "transition_matrix = np.array(((0.6, 0.3, 0.1),\n", | |
| " (0.3, 0.4, 0.3),\n", | |
| " (0.2, 0.3, 0.5)))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "The rows indicate the states the chain might currently be in and the columns the states the chains might transition to.\n", | |
| "If we take one \"time\" step of the Markov chain as one hour, then, if it's sunny, there's a 60% chance it stays sunny in the next hour, a 30% chance that in the next hour we will have cloudy weather and only a 10% chance of it raining immediately after it had been sunny before.\n", | |
| "This also means that each row has to sum up to one." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Let's run our Markov chain for a while:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "n_steps = 20000\n", | |
| "states = [0]\n", | |
| "for i in range(n_steps):\n", | |
| " states.append(np.random.choice((0, 1, 2), p=transition_matrix[states[-1]]))\n", | |
| "states = np.array(states)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "We can monitor the convergence of our Markov chain to its stationary distribution by calculating the empirical probability for each of the states as a function of chain length:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 5, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAFzCAYAAAB2A95GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZxcVZ3//9e5t/be093ZdyALZIEQFmVJgBEQGBVEVEAUQUe+iqMjoP5EDQruztcfDi7AT0RnQFHGGWQTFZDFBcIWQkjISval9632Or8/TnV1p5PuFKQr1Qnv5+NRj6pbd/vUrap7P/ecc8811lpERERE5MDyyh2AiIiIyFuRkjARERGRMlASJiIiIlIGSsJEREREykBJmIiIiEgZKAkTERERKYNAuQN4o84++2z78MMPlzsMERERkWKYwUYcdCVhTU1N5Q5BREREZL8ddEmYiIiIyKFASZiIiIhIGSgJExERESkDJWEiIiIiZaAkTERERKQMlISJiIiIlIGSMBEREZEyUBImIiIiUgZKwkRERETKoGRJmDHmZ8aYncaY5YOMN8aYm40xa4wxy4wxC0oVi4iIiBx6uru7Offcc5k/fz5z5szh17/+NVOnTi3cXWfp0qUsXrwYgCVLlvDRj36UxYsXM336dG6++WYANmzYwOzZs/nYxz7GUUcdxZlnnkk8Hmft2rUsWNCXmqxevXq34eFQyntH/hz4D+AXg4x/J3BE/nEC8OP8s4iIiBxEbvj9K6zY2jGsyzxyfDVf/eejhpzm4YcfZvz48TzwwAMAtLe38/nPf37Q6VeuXMljjz1GZ2cnM2fO5KqrrgJcgnX33Xdz2223cdFFF3Hvvfdy6aWXUlNTw4svvsjRRx/NHXfcweWXXz58H5ASloRZa58AWoaY5N3AL6zzd6DWGDOuVPEUa1dnksdW7aQzkS53KCIiIjKEuXPn8sc//pHPf/7zPPnkk9TU1Aw5/bnnnks4HKahoYHRo0ezY8cOAKZNm8bRRx8NwLHHHsuGDRsAuPLKK7njjjvIZrP8+te/5uKLLx7W+EtZErYvE4BN/YY359/bNnBCY8zHgY8DTJ48uaRBPb+xlX/55XM88OmTOWr80F+miIiIsM8Sq1KZMWMGzz//PA8++CDXX389Z5xxBoFAgFwuB0Aikdht+nA4XHjt+z6ZTGav78fjcQDe+973csMNN3D66adz7LHHUl9fP6zxHxQN8621t1prF1prFzY2NpY7HBERERkBtm7dSiwW49JLL+Xaa6/l+eefZ+rUqTz33HMA3Hvvvfu1/EgkwllnncVVV1017FWRUN6SsC3ApH7DE/PviYiIiOzTyy+/zLXXXovneQSDQX784x8Tj8e54oor+PKXv1xolL8/LrnkEn73u99x5pln7n/AA5QzCbsP+JQx5le4Bvnt1to9qiJFRERE9uass87irLPO2uP91157bY/3lixZstvw8uXL9/r6mmuu2W26p556issvvxzf9/cz2j2VLAkzxtwNLAYajDGbga8CQQBr7U+AB4FzgDVADzD85Xz7wdpyRyAiIiLldP7557N27VoeffTRkiy/ZEmYtfaD+xhvgU+Wav1vlil3ACIiIjIi/O53vyvp8g+KhvkiIiIihxolYSIiIiJloCRMREREpAyUhA1gjFqFiYiISOkpCRMREZFDxpIlS/je9743LMtavHgxS5cuHZZl7Y2SMBEREZEyUBImIiIiB61f/OIXzJs3j/nz5/OhD31ot3EvvvgiJ554IvPmzeP888+ntbUV2L2Eq6mpialTpwIQj8f5wAc+wOzZszn//PML95D82c9+xmc+85nCcm+77TY++9nP7nfs5ewxf0RTZ60iIiJFeugLsP3l4V3m2Lnwzm8NOckrr7zCjTfeyF//+lcaGhpoaWnh5ptvLoy/7LLL+OEPf8iiRYv4yle+wg033MAPfvCDQZf34x//mFgsxquvvsqyZctYsGABABdddBE33XQT3/3udwkGg9xxxx389Kc/3e+PqJKwAdQsX0RE5ODw6KOP8r73vY+GhgYARo0aVRjX3t5OW1sbixYtAuDDH/4wTzzxxJDLe+KJJ7j00ksBmDdvHvPmzQOgsrKS008/nfvvv5+VK1eSTqeZO3fufsevkjARERHZP/sosRppAoEAuVwOgEQiUdQ8V155Jd/4xjeYNWsWl18+PHdaVEmYiIiIHJROP/10fvOb39Dc3AxAS0tLYVxNTQ11dXU8+eSTAPzyl78slIpNnTqV5557DoDf/va3hXlOPfVU7rrrLsDd1HvZsmWFcSeccAKbNm3irrvu4oMfHPLOjEVTSdggLGoUJiIiMpIdddRRfOlLX2LRokX4vs8xxxxTaGQPcOedd/KJT3yCnp4epk+fzh133AHANddcw0UXXcStt97KueeeW5j+qquu4vLLL2f27NnMnj2bY489drf1XXTRRbz44ovU1dUNS/zGHmQt0BcuXGhL2WfHn1/dwRV3LuW+T53EvIm1JVuPiIiIHFzOO+88PvvZz3LGGWe8kdkGbW6u6kgRERGRIbS1tTFjxgyi0egbTcCGpOpIERERkSHU1tby2muvDftyVRImIiIiUgZKwgZxkDWVExERkYOMkrABjHprFRERkQNASZiIiIhIGSgJExERkUPalVdeyYoVK8odxh50deQg1CRMRETk4GGtxVqL5+1ZvnT77beXIaJ9U0nYAEa38BYRETkobNiwgZkzZ3LZZZcxZ84crrjiChYuXMhRRx3FV7/61cJ0ixcvprej98rKSr70pS8xf/58TjzxRHbs2EFnZyfTpk0jnU4D0NHRsdtwqagkTERERPbLt5/5NitbVg7rMmeNmsXnj//8PqdbvXo1d955JyeeeCItLS2MGjWKbDbLGWecwbJly5g3b95u03d3d3PiiSdy0003cd1113Hbbbdx/fXXs3jxYh544AHe85738Ktf/YoLLriAYDA4rJ9pIJWEiYiIyEFrypQpnHjiiQDcc889LFiwgGOOOYZXXnllr+3AQqEQ5513HgDHHnssGzZsAFy7sd57S95xxx1cfvnlJY9dJWEiIiKyX4opsSqViooKANavX8/3vvc9nn32Werq6vjIRz5CIpHYY/pgMIjJ90fl+z6ZTAaAk046iQ0bNvD444+TzWaZM2dOyWNXSdggDrYbm4uIiLyVdXR0UFFRQU1NDTt27OChhx56w8u47LLLuPjiiw9IKRgoCduT2uWLiIgcdObPn88xxxzDrFmzuPjiiznppJPe8DIuueQSWltb+eAHP1iCCPek6kgRERE5KE2dOpXly5cXhn/+85/vdbrHH3+88Lqrq6vw+sILL+TCCy8sDD/11FNceOGF1NbWDnuse6MkTERERN7yrr76ah566CEefPDBA7ZOJWGDUIswERGRt44f/vCHB3ydahM2gJqEiYiIyIGgJExERESkDJSEiYiIiJSBkjARERGRMlASNgj11SoiIiKlpCRsgN5bGYiIiIiUkpIwERERkTJQEiYiIiJSBkrCBqVGYSIiIlI6SsIGUIswERERORCUhImIiIiUgZIwERERkTJQEiYiIiJSBkrCBqHOWkVERKSUlIQNoL5aRURE5EBQEiYiIiJSBkrCRERERMpASdgg1CRMRERESklJ2ABG3bWKiIjIAaAkTERERKQMSpqEGWPONsasMsasMcZ8YS/jJxtjHjPGvGCMWWaMOaeU8YiIiIiMFCVLwowxPnAL8E7gSOCDxpgjB0x2PXCPtfYY4APAj0oVj4iIiMhIUsqSsOOBNdbaddbaFPAr4N0DprFAdf51DbC1hPG8IeqsVUREREopUMJlTwA29RveDJwwYJolwCPGmKuBCuCfShhPUdRZq4iIiBwI5W6Y/0Hg59baicA5wC+NMXvEZIz5uDFmqTFm6a5duw54kCIiIiLDrZRJ2BZgUr/hifn3+rsCuAfAWvs3IAI0DFyQtfZWa+1Ca+3CxsbGEoUrIiIicuCUMgl7FjjCGDPNGBPCNby/b8A0G4EzAIwxs3FJ2Igo6rJqFCYiIiIlVLIkzFqbAT4F/AF4FXcV5CvGmK8ZY96Vn+xzwMeMMS8BdwMfsWXOftQkTERERA6EUjbMx1r7IPDggPe+0u/1CuCkUsYgIiIiMhKVu2G+iIiIyFuSkjARERGRMlASNgg1yxcREZFSUhI2kFrmi4iIyAGgJExERESkDJSEiYiIiJSBkrBBqK9WERERKSUlYQMYNQoTERGRA0BJmIiIiEgZKAkTERERKQMlYSIiIiJloCRsEFbdtYqIiEgJKQkbwKhdvoiIiBwASsJEREREykBJmIiIiEgZKAkbjJqEiYiISAkpCRtATcJERETkQFASJiIiIlIGSsJEREREykBJmIiIiEgZKAkbhNrli4iISCkpCRvAqLdWEREROQCUhImIiIiUgZIwERERkTJQEiYiIiJSBkrCBmHVMl9ERERKSEnYAGqXLyIiIgeCkjARERGRMlASJiIiIlIGSsIGYdVdq4iIiJSQkrAB1CRMREREDgQlYSIiIiJloCRMREREpAyUhImIiIiUwT6TMGPM1caYugMRzEiizlpFRESklIopCRsDPGuMuccYc7Yxh3Z3pof2pxMREZGRYp9JmLX2euAI4P8DPgKsNsZ8wxhzWIljExERETlkFdUmzFprge35RwaoA35rjPlOCWMTEREROWQF9jWBMeZfgcuAJuB24FprbdoY4wGrgetKG2J5qEmYiIiIlNI+kzBgFHCBtfb1/m9aa3PGmPNKE1Y5qVGYiIiIlF4x1ZHTByZgxphfAlhrXy1JVCPAh3/2TLlDEBERkUNYMUnYUf0HjDE+cGxpwhERERF5axg0CTPGfNEY0wnMM8Z05B+dwE7gfw9YhCIiIiKHoEGTMGvtN621VcB3rbXV+UeVtbbeWvvFAxijiIiIyCFn0Ib5xphZ1tqVwG+MMQsGjrfWPl/SyMqkf2et1loO8b5pRUREpEyGujryc8DHgO/vZZwFTi9JRCNIzoKvHExERERKYNAkzFr7sfzzaQcunJElm7P4nrIwERERGX5DVUdeMNSM1tr/Hv5wRpZsTl22ioiISGkMVR35z0OMs8AhmYT1L/fKWiVhIiIiUhpDVUdefiADGYlUEiYiIiKlMlR15KXW2v80xvzb3sZba/+9dGGNDErCREREpFSG6jG/Iv9cNchjn4wxZxtjVhlj1hhjvjDINBcZY1YYY14xxtz1BmIvOSVhIiIiUipDVUf+NP98w5tZcP72RrcA7wA2A88aY+6z1q7oN80RwBeBk6y1rcaY0W9mXaWSU5swERERKZF93jvSGDPdGPN7Y8wuY8xOY8z/GmOmF7Hs44E11tp11toU8Cvg3QOm+Rhwi7W2FcBau/ONfoDh1r9z1oxKwkRERKREirmB913APcA4YDzwG+DuIuabAGzqN7w5/15/M4AZxpinjTF/N8acvbcFGWM+boxZaoxZumvXriJWPTxySsJERESkRIpJwmLW2l9aazP5x38CkWFafwA4AlgMfBC4zRhTO3Aia+2t1tqF1tqFjY2Nw7TqfVNJmIiIiJTKUFdHjsq/fCjfqP5XuP7B3g88WMSytwCT+g1PzL/X32bgH9baNLDeGPMaLil7trjwS0sN80VERKRUhuqs9Tlc0tXbSOpf+o2zuAb1Q3kWOMIYMw2XfH0AuHjANP+DKwG7wxjTgKueXFdc6KXRv7NWNcwXERGRUhnq6shp+7Nga23GGPMp4A+AD/zMWvuKMeZrwFJr7X35cWcaY1YAWeBaa23z/qx3OGWySsJERESkNIYqCSswxswBjqRfWzBr7S/2NZ+19kEGVF1aa7/S77UF/i3/GHFUEiYiIiKlss8kzBjzVVzD+SNxCdU7gaeAfSZhBzu1CRMREZFSKebqyAuBM4Dt+ftJzgdqShrVCKGrI0VERKRUiknC4tbaHJAxxlQDO9n9qsdDSr++WlUdKSIiIiVTTJuwpfm+u27DXTHZBfytpFGNEKqOFBERkVLZZxJmrf0/+Zc/McY8DFRba5eVNqyRQUmYiIiIlEqxV0deAJyM6x/sKUBJmIiIiMh+KOYG3j8CPgG8DCwH/sUYc0upAysX06+71qzahImIiEiJFFMSdjowO9+nF8aYO4FXShrVCJFVZ60iIiJSIsVcHbkGmNxveFL+vUOeSsJERESkVIa6gffvcW3AqoBXjTHP5EcdDzwz2HyHkpzahImIiEiJDFUd+b0DFsUIpc5aRUREpFSGuoH3X3pfG2PGAMflB5+x1u4sdWDlos5aRURE5EAo5urIi3DVj+8DLgL+YYy5sNSBjQQZNcwXERGREinm6sgvAcf1ln4ZYxqBPwG/LWVgI4Ea5ouIiEipFHN1pDeg+rG5yPkOemqYLyIiIqVSTEnYw8aYPwB354ffDzxYupBGDjXMFxERkVIp5t6R1/a7bRHArdba35U2rJHh9ebucocgIiIih6ghkzBjjA/8yVp7GvDfByakkeO2J9fzhXfOxvfMvicWEREReQOGbNtlrc0COWNMzQGKZ8S59/nN5Q5BREREDkHFtAnrAl42xvwRKNTPWWs/XbKoRpDmrlS5QxAREZFDUDFJ2H/zFqqKNANqHsfXRsoTiIiIiBzSimmYf6cxJgTMwt1LcpW19i1TPJRIZ8sdgoiIiByC9pmEGWPOAX4KrAUMMM0Y8y/W2odKHdxIEE8pCRMREZHhV0ynq/8OnGatXWytXQScBvzf0oZVPtt7tlA1+wv4sdUAJDK5MkckIiIih6JikrBOa+2afsPrgM4SxVN2q9uXARCseR5QdaSIiIiURjEN85caYx4E7sG1CXsf8Gy+A1estYdUo33P7J6XJtIqCRMREZHhV0wSFgF2AIvyw7uAKPDPuKTsEEvC8i/yzyoJExERkVIo5urIyw9EICOF7/WWhLn7RiYzSsJERERk+BXTJuwtxR/QUZiqI0VERKQUlIQN4HkDkzCVhImIiMjwUxI2QP+CsLkTapSEiYiISEkM2ibMGPNvQ81orf334Q+n/AL92oRFgp6qI0VERKQkhmqYX3XAohhB+ndREQn6dCUzZYxGREREDlWDJmHW2hsOZCAjhdevPjIc8GjqUkmYiIiIDL9i7h0ZAa4AjsL1GQaAtfajJYyrbPpfHRkO+uqiQkREREqimIb5vwTGAmcBfwEmcgjftqjv6khLJOCTVJswERERKYFikrDDrbVfBrqttXcC5wInlDas8imUhBmIBD22tMX5/iOryhuUiIiIHHKKScLS+ec2Y8wcoAYYXbqQyqt/m7BI0Afgh4+uGWxyERERkTelmHtH3mqMqQOuB+4DKoGvlDSqMjKmrzqyMlzM5hERERF544q5d+Tt+ZdPANNLG84I0K+z1smjYuWLQ0RERA5p+6yONMZ8wxhT22+4zhhzY2nDGhmqo8FyhyAiIiKHqGLahL3TWtvWO2CtbQXOKV1I5WX6FYUFBtxHUkRERGS4FJOE+caYcO+AMSYKhIeY/pAwvTHG+NpoYbitJ1XGaERERORQU0zL8/8C/myMuSM/fDlwZ+lCKq/ekrCZY6qYObbvzk1NXSlqY6FyhSUiIiKHmH2WhFlrvw3cBMzOP75urf1OqQMrF2N2r4K88NiJAHTrHpIiIiIyjIqpjsRa+5C19pr84w+lDmoksFgA3nP0BADefcvT5QxHREREDjGDVkcaY56y1p5sjOmEfEaSHwVYa211yaMrg/4N8wGiIb9MkYiIiMihbNAkzFp7cv65arBpDmXWurzzzV4gmcnm2NaeYNIgfY3lcrbffSpFRETkrWbI6khjjG+MWXmgghkJekvCeqsjQ4G+TZTL2b3Oszffe+Q1TvnOY+zsSOwx7o8rdjD9/3mQV7a272e0IiIicrAaMgmz1maBVcaYyQconvIbUDh11Pga5k9yfdV+/YEVRS/mL6/tAuD3y7bt9n57T5qP/WIpAOfe/BQPL9++H8GKiIjIwaqYhvl1wCvGmD8bY+7rfRSzcGPM2caYVcaYNcaYLwwx3XuNMdYYs7DYwA+k9+WvkLzj6Q2FasrNrT2s2t651+nbe9K8uq0DgK/f35e47exIMP9rj+w27Sf+8zne/R9P8ZE7nim8Z63l6TVNb6jkTURERA4uxfQT9uU3s2BjjA/cArwD2Aw8a4y5z1q7YsB0VcC/Av94M+spFdvvWoRIsK9x/s+e3sCssVVccrsLd8O3zt1j3uc3te42nMtZUtkcx3/jz4X3Hv3cIk7//l8AeGmzq5b844odvOPIMXz74VX85C9rAVh149mEA7o44M3qTKSJhQL4an834rX3pGmLp0hmciTTOXLWUl8ZoqEyTDjg7dF9zBthraUnlXVdzRiojgSJBH2yOcuG5m7iqSwNlWFqY8FB15XJ5vCMUVtOERk2xdzA+y9vctnHA2ustesAjDG/At4NDKzT+zrwbeDaN7meYVW4OrJfIVT/fe6anZ38dU1TYXhvDewvv+NZAI6fOopnNrQw88sP8c454wrjl17/TzRUhnnyutM45TuPFd7vrabsb+b1DwNwy8ULOHfeuD3GH8q2tMUJeobqaJCHlm/jgWXbOHvOON67YALGGKy1bGzpYVJdDM8z5HKWv61r5rYn15HK5NjSFuf15h6qwgHOnjOWT59xxG4XSqSzOZ5e08SWtjiHNVZy4vT6wjhrLTmLkrdBWOtOLDJZSyzkk8zkCknNqu2dvLCplWQ6RzydJZezNHenaOlOEQv5NHWlmFofo64iRFtPij+v3Im1sLGlh+wQpb/RoM9R46s5e85Y6itD1EZDJNJZtrTFmVAbJRL0Wb6lndaeNKMqguQsbM3/BpZvbacz0dfXn2dgVEWIlu4UA1fZUBkmEvQI+h6tPSm6kxnSWYtnoCIcYGJdjIqQz47OBAaDMXB4YyVT6iuYUh9jV2cS3zNEQz6xkE9jZZhoyCcc8JlcH2N8TUS/LREBSttFxQRgU7/hzcAJA9axAJhkrX3AGDNoEmaM+TjwcYDJk0vbPG1gFxUAQb+v1nbSqBibW+OF4V1dScZURwrDiXS28PqL58zi/B/9lXTWct9LWwH4wfuPpqEyXFjW2m+cw/aOBP/vn17jnqWbC/N+7JRp3Pbk+sLwJ+96npljT+WXf3udRTMbueLOpfzk0mM566ixbGuPs2JrB8dOqTsoe/V/7vVWvn7/Cpq7k2Sylm3tCUZXhdnZmdxj2j+9upNrfvMScyfU8PKWfV/YUF8RoiYa5DfPbebe5zezaEYjhzVWsmxLO8+sb9lt2nE1EY6fNoqaaJD/+sdGGivDHDGmkoBnmDm2mqMn1bC+qYdkJsv42ihbWuOsb+pme0eC1u4UcyfWcOS4aiJBn0Q6y2GjK1l0RCPpXA5rdy9RHaitJ8XK7Z2kszlmja2mLuaSiP4XhoBLfnpLaXoTxa5kBmstq3d2sXJbBxYI+R7pbI62njSJTJZM1nLc1FHEwj67OpMs39JOIp0jk7N0JzNsaO4mlclRGQ64xMlCVzLNzDHVjK4OY4CW7hTt8TRb2uJsbYuTzrrdQjToE09nCQc8ctYW3u8vEvQYFQvRk85SXxHiL6/tLEx34nS3zc88cgwzxlQRDnqEAz4G2NTaQyKfzLX3pFj6eis3PvDqkN95OOCRzOQAqIsFGVcT5dQZjUyrr6CuIoRv3B0wtrTFGV0VZtKoGLWxIDs6krR0J9nWniCbs3TE0xw7pY7GqjC+MWSte29rW5yeVJZ5E2oxBqyFtbu6eHJNE6lMDs+wR2K3t/jG1USYPa6aSXVRKsKutNYYg2dwJW7GdR5dEfKpjYUYXRUm4HvEQm57J9JZWrpTbG1LEAv59KTcd9CZyBDwDbPHVVEZDmKM62y6KhLE9wxjqyNURgJEAh4Bv6iuIgFXEtgbo7UWa8HzDNmcZVt7nHW7umnpTpHJWYK+IeB5BH1D0PcI5X8bPSkXdyKdJZPfSPFUlngqSyKTZWtbwiWwQZ9Uxn3vPakskaBHdTRINOgTDrjlVUWCJDNZtrbFCQd8JtVFmVAXK/wmjYFpDRWMrY6QyuboSmYIeKZQ0tqTztCdzBAO+ETy62vqStKZyNAeT7OrM4kxUBMNEgp4RIIe9RVhAvnP5nvuewp4bpznmXzCnmNnR5LWnjRdyTRb2xIk0lnGVEeoigQKJywh39BYFaGxKoS1kMzk8tsmRyzk05XMsL09wZiaCNWRAOGATyjgEQ265D4W8mmPpzHG4HuGgOe2te8Zgr6hIhzY7dglI5PpbeM07As25kLgbGvtlfnhDwEnWGs/lR/2gEeBj1hrNxhjHgeusdbuWRzUz8KFC+3SpUNOsl8e2/gYn37s0yyauIj/OOM/AOhIpFly3yv89/Nb9jpPLORz71VvZ/a4an7/0lauvvsFPnnaYVxz5kzec8vTherGaNDn1a+fPei61+zs4vXmbh5evp3vXDiP9niaq/7zeaoiAR5ZsWOv8yyYXMvzGwv3V+fr75nDh06cUhj+44odfOwXSzlyXDUPfPrkPapZcjnLb5/fzNT6CmIhnzkTaorbUPtgreWepZu4f9k2KkIB6iqCgOG8eeM46fAGcjnLtb9dxjMbmtnUEqci5NOTztL7c3zX/PEYA/cv20Yk4HHRcZP4l1MP44GXt/H4qp2kMjleb+5hQl2UnZ0JslnLuNooY6rDfPOCefSkMoR8j/rKMNZaVm7v5N7nNvOLv79OKn+AHlsd4dITJzNrbDUbmrv5zsOrSGVzhc8wtjpCMGDY1BLfyyd0wgGPMdURGqvCrN3VRVtPeq/TBTzDqTMaqa8IMaoiRDKTozORYV1TF0HP45kNLXud75jJtVgLzd1JupNZctYyY0wVrzd30x5Pk0jn9jrfQAMTg4BnCPjuADq6KuK2Y0eCMdURKsPuQNHUlaS5O8WuziSjKkLUxYLUxUKMrg5TEw3SWBkmHPTZ2ZGgNhaiqStJRTjA7HFVLJhcR0U4QMAzGAw1seBu8aSzOXpSWaJBf49EcyjWWtbu6iKVscTTWUK+x8S6KGt2deEZOLyxiuqoSyQ9Y4ZMfIdbVzJDVyJDfWUI3xgSmSxdyQxNnSkSGZdorN3VxaaWHgtCYnAAACAASURBVCJBn00tPazY1sG29gTxlPtuD3Qz0KpIAIMrlfPM7klgKp90BT2XgO7oSBaSOPfby2JxZ+WZYQjc9wxjqtwJak8+qY8GXYLUW5WczORIZrK7Jfo1UZeMDfZfCPle4USoHCpC7jM0dx/4+w/XxYLUVYTc/y2ZxfcMsZBPVSRILOTT0p3C9wyNVWEaq8JkspaAb4inslTnk9x42p3EpbI5Qr4rIa6rCDFjTGU+WcwypjpMzrqTSWtdiXFlOEBDZSh/zLHEQi4pDPkelZEAdbFgUU0MMtkciUyOcMBz+5P9aJZQRoMGXVQSli+xOhlXIvaUtfaFIuZ5G7DEWntWfviLANbab+aHa4C1QFd+lrFAC/CuoRKxUidhj296nKsfvXq3JKzXzOsfKpxh782Gb53LMV97hNaeNHd85DhOmzUagN8s3cS1v13G8hvOojJcTDO8PU39wgNFTzt3Qg05a3lla8ce49537ERuOn8uAc/w93XNXHz7nk3xrj79cLa2JdjQ3M3EuigvbWpjQ3NPYdkbmrqZUBflylOmc+qMBkZXRdjenmDZ5jaeWd/CP9a3DFlKNX9SLdWRAE+udtW6Hz1pGp87cwbG9J7B9v1eszlLIp2l4k1ut4E6Eml2diRJZXLMHle1xx96Z2ciXwphGFXhShW3tyfY3NpDS3eKURUhJtRF2djcw5T6CsbW9JWC9vYNt6G52y2rI8mqHZ35g2yG5Vs66EykSWZyhAIe8XSWoOdxxJhKTjq8gcNHVxLyPdbu6uIf61uIhXzaetL4nqEqEmBcTYSuZJbt7XGiIbcTsxamNlRQHQkwvjbKvIk1hHyPtnia2lgQzxiqI0EslmfXt9Lak2JaQwVHjKlUW8MRqreEM5uz5KylM5FhV2eSjkSaeDpLKpMjGvSJhnwqw+53kUjniAQ9MjlLRShAOpfjte2dhekrwgE6ExkMsL0jQVcyQ08qS0fcnTS4BNCt11pLKmOJBD2yOVtIsMbXRMjkLJtb44yqCBVKPj1jmNpQweRRMcbWRAh6LunJZC3prCttTeZLV6siAaL5pMQ3rio3EvSJBn0CXvHt7bI5S2ciTdD3qAgHsNbS1JViU2uPSyAzOTK5HGt3dbO5tYdwwKc2GiSXL0muiQYLpUmJtEvsgr7nSs5jQSpCAcbWRLAW2uIpMlm3H2rqSpHNWfewllx++/RWu7sSTRhdFaGhMkws7FMVDmCMIZ3N0Z3MEMyXUqfyJWZNXcnCCUMk6Er5elJu3zCxLkpzd5L2uNtvpPKlZT2pLN2pLBUhH88YLK4EOpO1ZHI50llLVyLDrq4Erd1pl9CGXHOBnlSWtnianmSGxqow2ZxlR0eCXZ1JggGPVMaVxHUmMu67Cvj4niEU8NyyM5btHQna43s/4SxWJOgVfit1sRC+Z3aLP5PNkc7Zwklzr95Sv4DntllVJEAmZwvNHWpjocLv0jPuBCHgGepiISzuBDSddfvg2miI84+ZwD8dOWa/PksR3nwSZoz5CvA+4L/zb70H+I219sZ9zBcAXgPOALYAzwIXW2tfGWT6xxkBJWG9SdipE0/lljNu2W3cvCV/oCPfrmRCbZQtbbuXkNx0/hxu/vNqsjnL0uvfMaxxdSZcsXNlOEAqfxD/0eNr+M7Dq1jyz0fykZOm0d6T3uPqS2PgZx8+jpsefJU1O7sGWTpMrY+xsaVnWM/Ej5lcy48uWYDBsKXNJXE/e3oDz6xvYVdnko+eNI1/PeOIPUpJDmW9O+2gb0ikc7ojg4gcdHqT3ngqS7RQouYuePE9Q3cyS0ciTUu3S1o9zxBPZUhlXVLV1pNie3uCbD5ZautJY60rhQv4HkHPPQd8QyTgTjjSGZfQZ3LuOZt1yW9nIoNnoL4yTFciQ2cyjcElprmcOwams5b2eAqDIZNzx89UJkd7PM2Vp0zng8eXvBeuQZOwYooXLgHmW2sTAMaYbwEvAkMmYdbajDHmU8AfAB/4mbX2FWPM14Cl1tqiurk40Aqdte4lOf30GUdw4wOvEgp4/OXaxVx8+z9Ytb2ToO/R1JXka79fQTqb45OnHT7scVVF+hKV3uqb/7P4cP7P4r511cSCPPOlM7j3uS20dCe55IQpTG2oAOC0WaNZsbWDc25+crflvvvo8fzg/UcXSoS2tsX57h9Wcdqs0YytjnDXP14n4Ht864K5BHyPlza1kbOWlza10dyd4oePriks63PvmMGZR43lpU1tvHPu2N1i7i0xOnbKKMAllf3Hv1V4niGUP9tXAiYiByNjXBVmr/6vAeorD3REB69ikrCtQATo7fo9jCvZ2idr7YPAgwPe+8og0y4uZpmlNlR9c3XUJQ2nHtFAwPe451/eBriE7bib/lS40qp/Q/0DbXRVhKsWH7bXcUeOr2bDt84ll7Osb+5mQ1M3Z8zevRh2fG2U//v+owvDx08btdv43o5rj5lcB8Dnzpy5x3pmjt33na7eigmYiIhIf8W0iG3Hddb6c2PMHcByoM0Yc7Mx5ubShjeyVEdczjqwkMwYwyUnTClU5c0YM7Jvt+l5hsMaK/dIwEREROTAKaYk7Hf5R6/HSxPKyPLklif3eG+oy31n9Sv9mVgXLUlMIiIicugoprPWOw9EIAeDuRNrmDGmkk/spbqvt6PPoG/KWh0pIiIiB4d9JmHGmPNwvdpPyU9fbGeth5zRVREe+eyivY6rqwjt9RZGIiIiIntTTHXkD4ALgJdtqXp2HUHS2f3r+0RERESkGMU0zN8ELH8rJGAANeHh6TFeREREZCjFlIRdBzxojPkLULiZn7X230sWVRkFPLdJxleML3MkIiIicigrJgm7CXdroQhw8N0d+k2yvCUK/kRERKRMiknCxltr55Q8khEmZ4u7MbKIiIjIm1FMm7AHjTFnljySEUYlYSIiIlJKxSRhVwEPG2MSxpgOY0ynMaaj1IGVnXIwERERKaFiOmsd2ffgKZEcqo4UERGR0imms1YDXAJMs9Z+3RgzCRhnrX2m5NGV0VukRw6Rg0/nDgjFIHwQnh/mcuDlKyCshc7tkE1C104wHkTr8tNl3XM2CTYHfggyScim3HMm4Z79EIya5rZFIAKRGjCm3/qykO6BVA/YrBvu3gVdO6BjKyQ7oLdvRGPc62wK0gk3fSACXgByaTcu1Q3BqHuvcjRUjoGKRqid4r6TWL0bb61bbyYff7LDLbM37kDYxZ5NuuVWjnHLC8b6YilG7356sOmthWQneD6EKor/nkQOkGIa5v8IyAGn43rO7wJuAY4rYVxlt0ebsFsXQy4Dn3iqLPEMm84d4AchNqrckRyakl3Qvhm6tsO2Ze7AWjka2l53B6FYvTtYBaLuQJhog5pJ0DADKhogUgurH4FEuzuYZdPQOAMmHg/hSuhpcQfsZEf+YJj/Lm0OKkb3HeD3pqfFHfRrJoAfhmAEMinYtRJ2rYKeJuhucjF5QXfA9IMQqoTm1ZDNuINnNuX+C9E6Ny7e6ubtackfpGNuXO/8kVqIVLvXbZugerz7fJmk+0yB/L1WA2E3TbLTbcdoHaS63LoySTfPzlch1emmrxwDmL7EI9nlEpJA2E2fy+YTgwaoP8zF6vluHX7IxWd8t5xonUsaQpUu+cimXeKRSblxfsht486t7nWswR3UsymI1rpxiXZIdLj1t29229IPQPsWF/fOFe59P+S2RybptvVw8oJue+TSkI67+N4oP+QeXsBto2w6/14wv21TbhulugZZgOHNt+cw7j/jB91/JToKAiH33Wbi+UQun4TmMi6O6Cj3vWdS/abJP+fSfdsgVOmSSj/kvhdMXyJvTF/C6AXc+r2gmz7V5T6P8dx4a926vYB7Bve77/+5IzVueuPl/8cp9zoQgXC1+++lut3nClW49fVOb4yL3Q+631423Ze8ZpP5OK37LqJ1bvmBcH65+QTZ8138oZhbZjbt3sO4WHqTbc9368rmE/xEu9uOvcO9r8EtJ5N0nzmXzf8vE249uYwb73nuRKPwefP/aZvrm6/3v0U+cfZ8991Eqt1nCITz02fd9ow1uO0ZCPd9j+Eq9/BDbv7oKLdNje8+f+9+KtXt9ifxFvd/zCTd9g5E8tsoAGPnuX1smRSThJ1grV1gjHkBwFrbaow5ZLuq6E2+disJW7LvDlyttSzdsZSFYxZiij2LO5A2Pwev3gdP/8ANR2ph/gfh5M9C1Rj3XjrhDkS7VsGUt7mz57//CJbfC7POg+px7g+w9UX3g08n3A6jZZ3bUR33MRhzlHuvZnJ+R4dLGsJVfWfIW5+H6ol9691f25a5JGD6aUOfEe9cAX/+mksW4q3QOBOmL3Z/WC/g/rTtm2Hj313iVD0B6qa4g3TbJjcu0QYNR7hl9DS7UgUv4BKLnSvh9RIl6cZ3B/ue5sGnqT8cxs51O9K2jW4HlehwO1Ev4BKw/gfHUFVfQlNYj5c/iGfddkl2uh1t7RQ3rx9ySZMxsOMVt5OraHAHzNrJLgFLdbsYIjV923RHm3uv/nBY/6TbEYYqoHVDfoedT6LA7YyDFbDlORdDqMI9Bytg/gegdpLboTevBYxLSMEto21j/oCIm6d7l/telv1qeL6HN8MLwKjpMHEhzLnAHYQSbS4Rbpzl/heV+f9C146+JADyB2ffbaNA2M0TCLnt3JtItqxznzkdd/+DeJsb15vsB/MPPwgY931VjoWqsX0JJtb9R/xg8aVQqW4Xb+cOd5KRSbj/Rbonn0iF8kmOccl2MOYOfoGwm9dm3efxA9CxzcWe7MofgNPQ0+oOnpmE+y/2/g56l+EF3YG0c5v7f1ZE3H6pdz3BqBsfq3e/566dLrZctq9UMNXlPrfN5kvhrFt/NpNPwpNQNc59nt5EDOOWm033JRjh6r7x1rrvt3e56Xi/0sF4/oQs/164yp2wZdN9CZ7Nufizqb7fSaDfw88nI/FW9933JtvJjr6k8A0z7ncXre37jfV/BhdX72/J+G4bBMJuO3p+X+wYN18um0/4UvnEMNCX+HjBvlXnMvkTr07oXue2uefntzXQ9bT7nt7MCUUx3vG1EZ+EpY0xPvm9tzGmEQ79BlOFkrAiqyWf3PIkn/zzJ7lizhV85tjPlDCyQexc2XfmCvDk92HLUndQ2ptEG/zjx/DCf8Ixl8DLvxn6AL9zxb5jWP9E3+tgzMVSNR46Nrvh+sNg+8t900xY6P7QE4+DEz6RT3BmDH0Q6G6Gx25y60p2uJ1V2+tuXLjG/Zl6E5HmtS4RSHW7bdHLeO4gtObPsPL+PdfRWyWS7tn9/d7kI93jdkyR6r4d0mt/cEnI/IvhsNPdAbB2iquqad8MlY1uW6Q6XWzxVpeQBKKuJKpto9v+8RYXf8NMN1/VGJeMbF8Oreth9JHuIBod5XZWxnM78US72yabl7qDzujZ7rebTefPIkMunppJLjHp3O52fhUN7v3xR7tYI7V7lqZl0+63NRJPLoqV6nHJaC6bT9pMX0lj+6b8QTHmdva9pQie735fiba+KrtIbb50KH/Q8wLu4NFb4hGrd99L9QSXWGVTLqH3/NJ9tknHl27ZQwlVuORy1HR30ibll81XG2Pd7zDd4/YDXiCfINl84ht0j2w6n9C/geS7XGz+RCHd3Ze09ZbKdze71zZfQpdO5E/gKt3vNFbv9m+9paq9+4Jc1o0rI7Ovtk/GmEuA9wMLgDuBC4HrrbW/KX14e1q4cKFdunTpvid8k17Y+QKXPXQZNeEanvrAU+7LvKlfic2xl8M//2CP+T796Kd5bNNjbhkfeqHQ8/6wy2Xhvqvh9b+6g/IbcdEv4Yh3uAPL6391B/+//Yc7m+o1dq47IDWvccMf+h93xvz0zTDnve6A0rLOHWDGzXcJjh90B6r2zW6ZK+5z49f8ESa/zZ3tb1nqfvzpHldi1bUDVj7QV9TdK1QFE45xZ/O1k2HDk+6sdOwcV63TsdlNVzvZJRc9Le4gMOl418Zl49/ddglVuiqK9s3ubHPSCTD+GPf5pp6ST6bi8NLd7s+5+VmY8U43TTDi1pnscEXYLeugbqpLEG3WHYRDlQPa3uT2TF5EREQKda97GVFMA3RjzCzgjPyC/mytfXX4YntjDlQSVhWq4q8f/Ks7CH/3sN2rbpa07zHf3DvnFl7ff/79TKmeMryBWQur/whP/Tts/NvQ046ZC+/+oUsachmX0NQNEk8mCa/+3iVnkX7Vrh1b88Xw+3F2ZO3Q82czripiy3Pw8m+habU7e2nd4EoV2ja5Up7Jb3dtPFrWuRKfd/0HLPjQ3peZy7qSsdqpSopERGQkGPRAWFRxjbV2JbBy2MIZwQpJaW9u2tv49OxvkPz91fx7XR1X9eyiNtZYmGdVy6rdlrG2be3wJGGpbvjGIPewPPyfXInOmKNg2qmuhKn3yqqBIkO0aQuEYe6Fe75fPQz3ztxXAtfbZmzCse4xHDzflYyJiIiMcCoqGEShTVihsXANT805l7tqqvj200tcg+dHb4J0nBXNrr3U5KrJAKxrX7eXJRazUgt3vgu+PdW1GXrouj2nidXDp5bCpffCyZ/JVy+GB0/AREREZEQqUcOlg18hCestCQtVEhk7H9a8zM4Nj0HXda49UdMquqe4qsj/Oue/OOXXp3D7y7dz5dwrB1nwgH5tXnvEVbltfQGW/qxvum9Pdc9+CBZ9Ho750PBdTSgiIiJlpyRsEDmbg6Y18I+fujdClWQbjoA1sNaD+LJf8bapk8jGX4CVLwBQE3bVft3pbh5e/zBnTzvbzdvbOeF9n4bn73TvVU/MX9LdtOfKa6f0XfH32RXuyjoRERE5pCgJG8qti/qVhMVI5CoBaA74XDJ+DNkBbZ5MOs4NR17JV1fcztNbn+bsSafBjaO5ua6G22pruGfLNmb3Ttx7lV9//3QDLLjMteFaeT8ccZa7Uk9EREQOOUrCBtits9b+PUIHonSnuwuDq0O791d7aXsHfGMcFwC3TxzH/6z5H97VnWCy73NbrSshu2jCOF6uOQVe/K++GT9wF0w5yTXCr5nQ9/6R7x72zyYiIiIjh5KwQexx26KGw+naOXjXEFe39nVbkcmXkH1028MwecJu0zWfcT317/nRnguI1r75YEVEROSgo6sjB5Gze94UoCPlbo/imb7N9v1F3+flc39HrLfB/Tu/y9Ute94P7oYTvwLA33Y9X4JoRURE5GCjJGwQFutuPdLPT176CQAnTzi58F5dpM7djgdcr+onfJzzzruNK9va8a3lozVH8cT7n+A9M94LwBef/OJu96VsS7Rx499v5AP3f4C5d85l3p3zWPLXJXskga80v8KyXctK8ElFRESkHFQdOcBunbWOPwY6tuwxzZK3LeH035wOwIy6/L0O+/Wib2afx7+Oncu/vnwPnHJNoTuKw2sPZ03bGub9Yh6fWfAZfvD8nrc/sljuXX0v966+F4APH/lhTp14Klc8ckVhmktnX8rnj/88y3Yt41crf8UNb7+BYP6Gv9979nvctfIuvnnKN1k0cRGRwN4b9mdzWXb07GB85Z6dsmZyGXzjj8wbkYuIiBwiirpt0UhS6tsWLd2+lMv/cDm+8XnRnwGr/+BGLGnn7HvP5pjRx/DNU77JL1f8krVta1ny9iVFL3tT5ybO+e9zBh1/3XHXcfvLt9OSaHnDcX9mwWdoT7Vzx/I7dnv/c8d+jvMOO49RkVG0J9sJekEe2vAQX/vb1wrT1IZr+fSCT3P6pNO5ddmt3LPqHjI2w/iK8Zw2+TQ+Ouej1Efq8YzHI68/QmO0kaZ4E5FAhO8v/T7r2tdx4rgTec/h76El0cKatjWMjY3l3OnnMrl68m7x7OjeQUeqg+VNy1k4diGTqia94c8qIiJyENm/e0eOJAcqCTMYlplpsO5xN2JJOyfedSIXHHEB1x23l57si/TUlqe46k9XFYZ//57fM7FqIs9se4a3jX8bxhistaxoXkHQD/Le+1w15mVHXsY1C6/hr1v/yif+9Ikh13HjSTfyzWe+udvVnOUyvWY6iyYuYnnzcnb27OT1jtd3G7944mLePuHtTK6aTCaXYWPnRtqT7SwYvQAMbOnaQnWomjMmn1G6m6KLiIiUzv7dO/KtyGIhnXAD7/wuXakuutPd1IX37/ZAJ084mccvepzrnriOf1v4b0ytmQrA2ye8vTCNMYajGo4C4OUPv7zb/CdNOInHLnqMzz72WW446Qam10zn8U2Pc/WjVwPwv+/+X6bXTufdh7+bVDbFt5/5Nve8ds8ecdx97t0cWX8kO7p3EM/GueWFW3jk9UcAeOL9T1AXqeOZbc9w98q7aUu28dyO57BYDq89nFmjZvHs9meZ0zCHd0x5B+dMO4edPTvZ3rOdF3e+yNGjj2Zl80ruee0eXmt9rXAbp9HR0UyumszRo4/mlImncP/a+/nH9n/w+ObHi9p2teFaPOPx9vFv57ixx1ETqqEh1sCR9UcS9ILs7NmJtZYxFcN7Z4GtXVsxGOoidYNW75ZaIpOgK91FQ7ShLOsXERluqWyKnM2xrn0dGzs2kswmCfkhQl6IoB8k5IcwGCyWeDpOxmYwGDpTnWzu2kxTvIkd3TvI2AxhP0xPuofWRCtZmwUg7IfdIxAuLDeRSZCxGZLZJJ2pTj585Id5b77NdjmoJGyA3pIwgJfXb4TDzmDXe39aaAP2rsPexU0n31Sy9ZdCOpcmnU0TC8ZIZBIEvSC+5x+QdW/u3My69nVMrppcSDj7y+QybOvaxt+3/53NnZuZWTeToB+kO93NqpZVzBo1i6pQFb9f+3v+tPFPAAS8AJlcprCMaCBKPBMvDM+pn0NVqIrqcDVTq6cytmIsa9vW8uLOF1nevJyZdTMZVzmOhmgDEysnMqFyAn/e+Gf+9PqfiAaidGe6OXXiqXh4LN2xdLerYsN+mHgmjmc8jqo/isnVk/GNT2WwknQuzbSaaQS8AKNjown7YZZuX0prspVoIEp7sp2Xdr1EU7yJTC7DmNgY5o+ez45u1zavJlxDwARY37GeMbExZHIZmuJNrGlbw5Yu1zaxIlhByAsxtmIsvvEJ+SGiwSjJTJKmeBPRQJTacC3GGIJekIZoA5OqJmGxdKW6qAxVks1lGV85nkggwuSqyUyqmkQsGCvhr0BEhls6m6Y73U06lyaZTWKtpT3VTnuyndZkK52pTjqSHXSmOunJ9BT2F737qngmTjKbdMvKpWlLtlEdqiYWiBUSl4gfKSQyAS+A7/kEvADtyXZWt64mmU2SsznimTgGQ33UNVtJZpNs795Oc6K5cOyBvjbXzfFmmhPNhWTpzfCMR32kntGx0QS8AKlsimggyqjIKIwx+MYnmU0WHqlsimQ2ScSPuATPC1Edrua86eexeNLi/f4+9kHVkcV6dvuzfPQPHwXgbxs28fapkxhfOaFwEDwYk7BDRc7m8IxHa6KV53Y8x5auLaSyKXb07ODJzU8yrWYadZE6dsV30ZZo4/WO10lkE4X5D689nO50N9FAlM5UJ7viu3ZbfjQQ5YSxJ7C6bTXbu7cT8kNMqJzAcWOPoy5SRyaXoTvdTVeqC894rGhewarWVYX5q4JVdKY794g75IVI5VKMjo7msNrDmFYzjbAfZmXLSta2raUh1sDGjo1kbZaczRV2jBXBCsJ+mHmN85hQOYH6SD0bOjYAsKNnB+3JdmKBGC2JFuoidTREG4hn4rQmWulKdxFPx9kV31Xo8673jHIgg2FOwxzqo/WMr3DJYMgPkc6mWdG8gqZ4E1NqpjCpahKN0UYao42Mjo0mZ3N0pbsYFRlFfbSemnANBvOWrDbO5rJs7d5KT7qHqlAVjbFGgl6QZDaJb/y9bpOczdGZ6mR793ZiwRg5myOTyxSS/aZ4EyHfnblv795OZ7qTTC5DOpumIdpQWEc6l2Zb9zZ849OT6aEn3YPBsKFjAxZLLBAjFnQH1vZkO/FM3B1oA+7gGvFd6W4mlyESiBD0gmRtlmggSkeqA2stFcEKxsTGEM/G6Ux1EvSCtCRaaEm0FBKASCDiTgr8KMYYooEo0UCUWDBGR7KDVDblSiP8EJFAhHEV40hmk2RymcJ2slgMhnQuTUeqg+5UN9FglEwuQ2eqk/poPb7xCwfvdDYNwPr29SSyCTpTnRgMQT9I0AuyunU18UycdC5NzuZIZBJEA1GCfrDQ1VAykySdS+MZj0ggQsSPEA1EiQTcs2c8etI95GyOVC5FwAsQ9IJ9JSqZJIlsAt/4hP3wbp+793UkEKEp3kRPuofOVCfxTJyczdGSaKEj1UE6m8YY998JmAAhP+Re5383qWyKVC5FxI+QzqVpTbTSk+kp6rfZG0MmlykkX57xiAVihZIm3/OpCdfQmeokkUmQzCZJZBJ73V/0qg3XEgvEMMZQEawgZ3M0xZsK239MbAyN0UYCXsAlaf0u9KoL19EYayTkhfA9nynVU5hWM41oIEo6myaVS7nPnE0VfhPRQLSwPWLBGGMrxhL0gkVtgxFASVix+idhe/Po+x6lMaZ7OR4MetI9bOzcSHO8mWNGH7NHac+unl1kchlWtKxgavVUDqs9rDAunU2TI0fYDw+5jpzNkc1l8YyHZzzWd6wvFIn3Vt+OrxyPtfYNXW2ayWWGLZlpT7aTyqZoiDYUzlw3dGzAWsvmrs2salnFczueozXZSlO8ic5UXyI5oXICo2Oj2dixkdZk6177z+svYAJMqp7EEbVHUBmqJJPL0JJoKSQD1aFqdzDMJoin41SFqpg5aibTaqYVqnxrw7XUhmv///buPkiyq7zv+PeZl57ZmZ3dndWuCJG0rLAUHPFiSWwpJAVOxSayoBzJsYkjjGNskagIJg5WxbESVVy2U5UY48SJK8QYB8XYls2LYsImgQjZ4BcoS0KgNZIQQoteIuSFFUi7q9XszGz3PPnjnhn19k7Pzkpz587L91PVNbdP33v7nHu7p3997ul7ac+18iGr7QAAE15JREFUGR4c5vjscbYMb2G6Pc2RmSM8fuxxHnvmMR479hhBMNGa4IGnHmBoYIiJ4YlTQkSSbGttY1trG9tHtrN1eCsznRmennmaIzNHGB4YZnx4nPGhcaY70/zl8b/ksWOPMRADzDHHUAyRJMMDw2QmnexUh71JZjozCz2jh6cOn9IbOzwwvNAjMDwwzEU7LmLnlp20O20OHjnITGeGqfbUGbfnCzHRmqA10GKqPbVQt5HBEUaHRheCwws1HxigCgrtbJ9hiXoMxiBbhrawrbWtqkv5EN+7fS+TI5MMDwwTEYwMjlThaa5NJzskSWugxcjgCEky3Z7mROdE9bdd/Z1/Lc2Hr/n30HyPyvDAMGNDYwuHuE6cPMGJ9qm36c40kyOTjA+Ps7W1lfHhcQajCj47RnbQGmwt/C+ZX09nrrPwxWxkcOS5UD8wyOTIZNXb39rG8MAwI0MjZCaTo5NsH9nOztGdTLQmqvfD4KlBZbm/fs9MTs6dZLozzUx7hk52aM+1ac+1GR8eZ9eWXf6CfvkMYct1phDWO0ZL2mjmeyZmOjNMjkwu/KOdD1RPTj3J4anDHJs9xu6x3RyfPc5T00/x9PTTzHRmeOToIzz49IML36R3jOxguj1Na7C18GOR1mD1wXf85HEOTx1+XvXcObqT9lyb6fY0e7btWRgX2M42W4e3Lnzjf2b2mYXexW7DA8O059qnfNsfGxpjz7bqRyKtwdbC4ZOTc1VPxcnOSXaP7V7oRZjtzDI6NMruLbu5ePJiJloTHJ05yqNHq16o+d6Frz79VZ6ZrXqy9m7byzlbzmHL0Ba2j2xn15ZdCz1g8wFytjPL1lZ1rdrRwVF2bdnF5OjkQjA/Mn2Ew1OH6WSH1mCL3VuqL4bjw+MLvWoTrYmF3p75Q0bzh5Xmy2Y7swuBY9vINqbb08x2ZokIptvTbB+pLrk21Z6qeocHWuwYrQLy1uGtp30hnX/tzPc6TbWnePbks0y0JtgytGXhkNBUe4pDxw8thJuRwRE62SEI5nKOoYGhhe0z25mt9s3wWDU0IFkIyaODo3Syw4vGXrQQBtei+V58bVoOzJe0PPOHQcaHx08rP3fsXM4dO5eX8/IVe77DU4c59OyhakxLGc9yZPoIz7afZWRwhNZAi052GBsaY/vIdi6YuIA92/Yw0ZpYdg/jTGeGYzPHODZ7jMEYrHrbRneQmZxon2CqPcXwwDDbWtvWxbf7XVt2cdHkRcuefyAGTtufC4ffhkbZQXXZtN555k0yyXlbz1v0sW7dh9AmWhNLzvvyc87+NbRef5hiAFM/hjBJjZoPds/HcgPTyOAIu8d2n9ZzExGMDY/5wwRJjTCen4XPXvvZpqsgSZI2CEPYWZgfHyFJkvRCGcIkSZIaYAiTJElqgAPze/SesuPC0V381GtuWhe/mJIkSeuHIewM9r/senjJ65uuhiRJ2mA8HLmEf3zkKHSdBVuSJGmlGMKWEABzz/8Co5IkSf0Ywnp0X8IkElhnl3WSJEnrgyHsTGq8wK4kSdq8DGFLuHe0Ba/8B01XQ5IkbUCGsCU8Pr4Txs9puhqSJGkDMoQtYa413nQVJEnSBlVrCIuIqyLiwYg4GBE3LvL4DRHx5Yj4UkT8UUS8pM76LEf3wHxJkqS61BbCImIQeC/wBuAS4M0RcUnPbPcA+zLzVcCtwC/XVZ/nw7PkS5KkutTZE3YFcDAzH87MWeBDwDXdM2TmZzJzqty9Azi/xvqctTlPTyFJkmpSZwg7D3i86/7XS1k/bwM+WWN9zloHT08hSZLqsSauHRkRPwrsA/52n8evB64H2LNnz6rVa/vwxKo9lyRJ2lzq7Al7Arig6/75pewUEfF64Cbg6sycWWxFmfn+zNyXmft2795dS2W7nmth+r9e8W9qfS5JkrR51RnCPg9cHBEXRkQLuBbY3z1DRFwG/AZVADtcY13O2i88+W3+yti5TVdDkiRtULWFsMxsA+8EbgMeAD6SmfdHxC9GxNVltvcAW4GPRsSBiNjfZ3Wrbu/JNoSnUZMkSfWodUxYZn4C+ERP2c91Tb++zud/Xjrt56YNYZIkqSamjB7Znn7ujiFMkiTVxJRxmmpgfpCGMEmSVBtTxlIMYZIkqSamjNN0nSXfECZJkmpiyliKIUySJNXElNGr+3KRI1sbq4YkSdrYDGG9cn5gPjA00mhVJEnSxmUIkyRJaoAh7DR55lkkSZJeIENYjzSESZKkVWAI62UGkyRJq8AQdhpTmCRJqp8hrJ/v+pGmayBJkjYwQ9hpSk/Y+K5mqyFJkjY0Q1iPTA9HSpKk+hnC+ojqdK2SJEm1MISdpvSEhSFMkiTVxxAmSZLUAENYD8eESZKk1WAI6yPCTSNJkupj0jiNPWGSJKl+hrDTzIcwB+ZLkqT6GML6MYNJkqQaGcJ6ODBfkiStBkNYH+GmkSRJNTJpnMaeMEmSVD9DWD+OCZMkSTUyhPVyTJgkSVoFhrDTVCHMMWGSJKlOJo1+PBwpSZJqZAjr5eFISZK0CgxhfdkVJkmS6mMI65GeokKSJK0CQ1ivcjgywp4wSZJUH0NYX4YwSZJUH0OYJElSAwxhPRbGhHk4UpIk1cgQ1mthXL4hTJIk1ccQdhp/HSlJkupnCOvHw5GSJKlGhrBenjFfkiStAkNYDyOYJElaDYaw05STtTowX5Ik1cgQ1mv+cKQZTJIk1cgQ1pcpTJIk1ccQ1sMLeEuSpNVgCOsjwk0jSZLqY9Lo5SkqJEnSKjCESZIkNcAQ1pcD8yVJUn0MYT0y55qugiRJ2gRqDWERcVVEPBgRByPixkUeH4mID5fH74yIvXXW52yE146UJEk1qi2ERcQg8F7gDcAlwJsj4pKe2d4GPJ2ZFwG/Cry7rvosnwPzJUlS/ersCbsCOJiZD2fmLPAh4Jqeea4BPlimbwW+N5rugvrWQ2XCnjBJklSfOkPYecDjXfe/XsoWnScz28BR4JzeFUXE9RFxd0Tc/eSTT9ZU3co5576CvxVbGf+rr671eSRJ0uY21HQFliMz3w+8H2Dfvn21Hi+89JVv4Tde+ZY6n0KSJKnWnrAngAu67p9fyhadJyKGgO3At2uskyRJ0ppQZwj7PHBxRFwYES3gWmB/zzz7gbeW6TcBn870lPWSJGnjq+1wZGa2I+KdwG3AIHBzZt4fEb8I3J2Z+4EPAL8TEQeBp6iCmiRJ0oYX663jad++fXn33Xc3XQ1JkqTl6Hu6Bc+YL0mS1ABDmCRJUgMMYZIkSQ0whEmSJDXAECZJktQAQ5gkSVIDDGGSJEkNMIRJkiQ1wBAmSZLUgHV3xvyIeBJ4rOan2QV8q+bnWMs2c/s3c9thc7d/M7cdNnf7bfvmtRrt/1ZmXrXYA+suhK2GiLg7M/c1XY+mbOb2b+a2w+Zu/2ZuO2zu9tv2zdl2aL79Ho6UJElqgCFMkiSpAYawxb2/6Qo0bDO3fzO3HTZ3+zdz22Fzt9+2b16Ntt8xYZIkSQ2wJ0ySJKkBhrAeEXFVRDwYEQcj4sam67MSIuKCiPhMRHw5Iu6PiH9eyn8+Ip6IiAPl9sauZf5V2QYPRsT3dZWvy+0TEY9GxL2lnXeXsp0RcXtEPFT+TpbyiIhfK238UkRc3rWet5b5H4qItzbVnuWKiJd17d8DEXEsIt61kfd9RNwcEYcj4r6ushXb1xHx6vJaOliWjdVtYX992v6eiPhKad/HImJHKd8bESe6XgPv61pm0Tb2245rQZ+2r9jrPCIujIg7S/mHI6K1eq07sz7t/3BX2x+NiAOlfKPt+36fcWv/fZ+Z3soNGAS+BrwUaAF/AVzSdL1WoF0vBi4v0xPAV4FLgJ8H/sUi819S2j4CXFi2yeB63j7Ao8CunrJfBm4s0zcC7y7TbwQ+CQTwGuDOUr4TeLj8nSzTk0237Sy2wSDwDeAlG3nfA98NXA7cV8e+Bu4q80ZZ9g1Nt/kMbb8SGCrT7+5q+97u+XrWs2gb+23HtXDr0/YVe50DHwGuLdPvA/5p020+U/t7Hv8PwM9t0H3f7zNuzb/v7Qk71RXAwcx8ODNngQ8B1zRcpxcsMw9l5hfL9DPAA8B5SyxyDfChzJzJzEeAg1TbZqNtn2uAD5bpDwI/0FX+21m5A9gRES8Gvg+4PTOfysyngduBRU/At0Z9L/C1zFzqZMfrft9n5p8CT/UUr8i+Lo9ty8w7svrP/Ntd62rcYm3PzE9lZrvcvQM4f6l1nKGN/bZj4/rs937O6nVeej2+B7i1LL+m2g5Lt7/U/4eB319qHet43/f7jFvz73tD2KnOAx7vuv91lg4r605E7AUuA+4sRe8s3bE3d3Uv99sO63n7JPCpiPhCRFxfyl6UmYfK9DeAF5Xpjdh+gGs59Z/wZtn3sHL7+rwy3Vu+XlxH9S1+3oURcU9E/ElEvK6ULdXGfttxLVuJ1/k5wJGuMLve9vvrgG9m5kNdZRty3/d8xq35970hbBOJiK3A/wDelZnHgF8HvgO4FDhE1V29Ub02My8H3gD8ZER8d/eD5dvNhv2pcBm/cjXw0VK0mfb9KTb6vu4nIm4C2sAtpegQsCczLwNuAH4vIrYtd33rZDtu2td5jzdz6hewDbnvF/mMW7BW62wIO9UTwAVd988vZeteRAxTvThvycw/AMjMb2ZmJzPngN+k6oqH/tth3W6fzHyi/D0MfIyqrd8s3czz3fCHy+wbrv1U4fOLmflN2Fz7vlipff0Epx7OWxfbISJ+HPh+4C3lw4hyKO7bZfoLVGOh/hpLt7HfdlyTVvB1/m2qQ1ZDPeVrXqnzDwIfni/biPt+sc841sH73hB2qs8DF5dfwbSoDt/sb7hOL1gZD/AB4IHM/I9d5S/umu3vA/O/qtkPXBsRIxFxIXAx1aDEdbl9ImI8Iibmp6kGKt9HVff5X7+8Ffh4md4P/Fj5Bc1rgKOlS/s24MqImCyHNa4sZevBKd+EN8u+77Ii+7o8diwiXlPeVz/Wta41KSKuAv4lcHVmTnWV746IwTL9Uqp9/fAZ2thvO65JK/U6L8H1M8CbyvJrvu1dXg98JTMXDqdttH3f7zOO9fC+P5tR/JvhRvWria9SfTO4qen6rFCbXkvVDfsl4EC5vRH4HeDeUr4feHHXMjeVbfAgXb8CWY/bh+qXTn9RbvfP15tqnMcfAQ8BfwjsLOUBvLe08V5gX9e6rqMaxHsQ+Imm27bM9o9TfZPf3lW2Yfc9Vdg8BJykGrvxtpXc18A+qg/zrwH/hXLS67Vw69P2g1TjXObf++8r8/5QeT8cAL4I/L0ztbHfdlwLtz5tX7HXefk/clfZnh8FRppu85naX8p/C3h7z7wbbd/3+4xb8+97z5gvSZLUAA9HSpIkNcAQJkmS1ABDmCRJUgMMYZIkSQ0whEmSJDXAECZpXYuIP46IfavwPD8VEQ9ExC3LmHdHRLyj7jpJWt8MYZI2ra4zoC/HO4C/m5lvWca8O8r8ktSXIUxS7SJib+lF+s2IuD8iPhURW8pjCz1ZEbErIh4t0z8eEf8zIm6PiEcj4p0RcUO56PAdEbGz6yn+UUQciIj7IuKKsvx4VBdtvqssc03XevdHxKepTuTYW9cbynrui4h3lbL3UZ2s85MR8dM987+8PMeBqC4UfTHwS8B3lLL3lPl+JiI+X+b5ha7t8pWIuKVsn1sjYqw89ksR8eUy/6+s2M6QtGaczbdASXohLgbenJn/JCI+QnXW7t89wzKvAC4DRqnOYP2zmXlZRPwq1aVD/lOZbywzL43qwuw3l+VuAj6dmddFxA7groj4wzL/5cCrMvOp7ieLiFcDPwH8Daqzat8ZEX+SmW8vl//5O5n5rZ46vh34z5l5S7nUzSBwI/CKzLy0rPfK0v4rynr3l7r+P+BlVGc3/1xE3Ay8IyL+O9Vldr4zM7PUX9IGY0+YpNXySGYeKNNfAPYuY5nPZOYzmfkkcBT4X6X83p7lfx8gM/8U2FZCy5XAjRFxAPhjqiC3p8x/e28AK14LfCwzn83M48AfAK87Qx3/HPjXEfGzwEsy88Qi81xZbvdQXSbmO6lCGcDjmfm5Mv27pQ5HgWngAxHxg8AUkjYcQ5ik1TLTNd3huZ74Ns/9LxpdYpm5rvtznNqT33v9taTqcfqhzLy03PZk5gPl8WefR/0XlZm/B1wNnAA+ERHfs8hsAfz7rrpclJkf6Ff3zGxT9ZrdCnw/8H9Xqr6S1g5DmKSmPQq8uky/6Xmu4x8CRMRrgaOZeRS4DfhnERHlscuWsZ4/A34gIsYiYpzqkOCfLbVARLwUeDgzfw34OPAq4Blgomu224DrImJrWea8iDi3PLYnIv5mmf4R4LNlvu2Z+Qngp4HvWkbdJa0zjgmT1LRfAT4SEdcD/+d5rmM6Iu4BhoHrStm/pRoz9qWIGAAeoepV6iszvxgRvwXcVYr+W2bec4bn/mGqHwacBL4B/LvMfCoiPhcR9wGfzMyfiYi/Dvx5yYTHgR+l6hF8EPjJMh7sy8CvA9uBj0fEKFUv2g3L3RCS1o/I7O0JlySthojYC/zvzHxFw1WR1AAPR0qSJDXAnjBJkqQG2BMmSZLUAEOYJElSAwxhkiRJDTCESZIkNcAQJkmS1ABDmCRJUgP+P7r39Ap3yEpIAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 720x432 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "def despine(ax, spines=('top', 'left', 'right')):\n", | |
| " for spine in spines:\n", | |
| " ax.spines[spine].set_visible(False)\n", | |
| "\n", | |
| "fig, ax = plt.subplots()\n", | |
| "width = 1000\n", | |
| "offsets = range(1, n_steps, 5)\n", | |
| "for i, label in enumerate(state_space):\n", | |
| " ax.plot(offsets, [np.sum(states[:offset] == i) / offset \n", | |
| " for offset in offsets], label=label)\n", | |
| "ax.set_xlabel(\"number of steps\")\n", | |
| "ax.set_ylabel(\"empirical probability\")\n", | |
| "ax.legend(frameon=False)\n", | |
| "despine(ax, ('top', 'right'))\n", | |
| "plt.show()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## The mother of all MCMC algorithms: Metropolis-Hastings" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "So all that is fun, but back to sampling an arbitrary probability distribution $\\pi$. \n", | |
| "It could either be discrete, in which case we would keep talking about a transition matrix $T$, or it could be continous, so that $T$ would be a transition *kernel*.\n", | |
| "From now on, we're considering continuous distributions, but all concepts presented here transfer to the discrete case. \n", | |
| "If we could design the transition kernel in such a way that the next state is already drawn from $\\pi$, we would be done, as our Markov Chain would... well... immediately sample from $\\pi$.\n", | |
| "Unfortunately, to do this, we need to be able to sample from $\\pi$, which we can't.\n", | |
| "Otherwise you wouldn't be reading this, right? \n", | |
| "A way around this is to split up the transition kernel $T(x_{i+1}|x_i)$ into two parts:\n", | |
| "a proposal distribution $q(x_{i+1}|x_i)$, from which we can sample possible next states, and an acceptance probability $p_\\mathrm{acc}(x_{i+1}|x_i)$ which determines if a proposal $x_{i+1}$ is accepted as the next state in the chain or not.\n", | |
| "This acceptance / rejection step corrects for the error introduced by proposal states from $q \\neq \\pi$.\n", | |
| "We thus have \n", | |
| "$$\n", | |
| "T(x_{i+1}|x_i)=q(x_{i+1} | x_i) \\times p_\\mathrm{acc}(x_{i+1}|x_i) \\ \\mbox .\n", | |
| "$$\n", | |
| "A sufficient condition for a Markov chain to have $\\pi$ as its stationary distribution is the transition kernel obeying *detailed balance* or, more intuitively, *microscopic reversibility*:\n", | |
| "$$\n", | |
| "\\pi(x_i) T(x_{i+1}|x_i) = \\pi(x_{i+1}) T(x_i|x_{i+1})\n", | |
| "$$\n", | |
| "This means that the probability of being in a state $x_i$ and transitioning to $x_{i+1}$ must be equal to the probability of the reverse process, namely, being in state $x_{i+1}$ and transitioning to $x_i$.\n", | |
| "Transition kernels of most MCMC algorithms satisfy this condition. \n", | |
| "For the two-part transition kernel to obey detailed balance, we need to choose $p_\\mathrm{acc}$ correctly, meaning that is has to correct for any asymmetries in probability flow from / to $x_{i+1}$ or $x_i$.\n", | |
| "One possibility is the Metropolis acceptance criterion: \n", | |
| "$$\n", | |
| "p_\\mathrm{acc}(x_{i+1}|x_i) = \\mathrm{min} \\left\\{1, \\frac{\\pi(x_{i+1}) \\times q(x_i|x_{i+1})}{\\pi(x_i) \\times q(x_{i+1}|x_i)} \\right\\} \\ \\mbox .\n", | |
| "$$\n", | |
| "Often, symmetric proposal distributions with $q(x_i|x_{i+1})=q(x_{i+1}|x_i)$ are used, in which case the Metropolis-Hastings algorithm reduces to the original, but less general Metropolis algorithm developed in 1953 and for which\n", | |
| "$$\n", | |
| "p_\\mathrm{acc}(x_{i+1}|x_i) = \\mathrm{min} \\left\\{1, \\frac{\\pi(x_{i+1})}{\\pi(x_i)} \\right\\} \\ \\mbox .\n", | |
| "$$\n", | |
| "Formally, we can then write the complete Metropolis-Hastings transition kernel as\n", | |
| "$$\n", | |
| "T(x_{i+1}|x_i) = \\begin{cases}\n", | |
| " q(x_{i+1}|x_i) \\times p_\\mathrm{acc}(x_{i+1}|x_i) &: x_{i+1} \\neq x_i \\mbox ; \\\\\n", | |
| " 1 - \\int \\mathrm{d}x_{i+1} \\ q(x_{i+1}|x_i) \\times p_\\mathrm{acc}(x_{i+1}|x_i) &: x_{i+1} = x_i\\mbox .\n", | |
| " \\end{cases} \n", | |
| "$$\n", | |
| "A step is then performed as follows: a proposal state $x_{i+1}$ is drawn from $q(x_{i+1}|x_i)$.\n", | |
| "It is then accepted as the next state with probability $p_\\mathrm{acc}(x_{i+1}|x_i)$ or rejected with probability $1-p_\\mathrm{acc}(x_{i+1}|x_i)$, in which case the current state is copied as the next state. \n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Alright, now that we now how Metropolis-Hastings works, let's go ahead and implement it.\n", | |
| "First, we set the log-probability of the distribution we want to sample from - without normalization constants, as we pretend we don't know them.\n", | |
| "Let's work for now with a standard normal distribution:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 6, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def log_prob(x):\n", | |
| " return -0.5 * np.sum(x ** 2)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Next, we choose a symmetric proposal distribution.\n", | |
| "In general, including information you have about the distribution you want to sample from in the proposal distribution will lead to better performance of the Metropolis-Hastings algorithm. \n", | |
| "A naive approach is to just take the current state $x$ and pick a proposal from $\\mathcal{U}(x-\\frac{\\Delta}{2}, x+\\frac{\\Delta}{2})$, that is, we set some step size $\\Delta$ and move left or right a maximum of $\\frac{\\Delta}{2}$ from our current state:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 7, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def proposal(x, stepsize):\n", | |
| " return np.random.uniform(low=x - 0.5 * stepsize, \n", | |
| " high=x + 0.5 * stepsize, \n", | |
| " size=x.shape)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Finally, we calculate our acceptance probability:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 8, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def p_acc_MH(x_new, x_old, log_prob):\n", | |
| " return min(1, np.exp(log_prob(x_new) - log_prob(x_old)))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Now we piece all this together into our really brief implementation of a Metropolis-Hastings sampling step:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 9, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def sample_MH(x_old, log_prob, stepsize):\n", | |
| " x_new = proposal(x_old, stepsize)\n", | |
| " # here we determine whether we accept the new state or not:\n", | |
| " # we draw a random number uniformly from [0,1] and compare\n", | |
| " # it with the acceptance probability\n", | |
| " accept = np.random.random() < p_acc_MH(x_new, x_old, log_prob)\n", | |
| " if accept:\n", | |
| " return accept, x_new\n", | |
| " else:\n", | |
| " return accept, x_old" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Apart from the next state in the Markov chain, `x_new` or `x_old`, we also return whether the MCMC move was accepted or not.\n", | |
| "This will allow us to keep track of the acceptance rate.\n", | |
| "Let's complete our implementation by writing a function that iteratively calls `sample_MH` and thus builds up the Markov chain:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 10, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def build_MH_chain(init, stepsize, n_total, log_prob):\n", | |
| "\n", | |
| " n_accepted = 0\n", | |
| " chain = [init]\n", | |
| "\n", | |
| " for _ in range(n_total):\n", | |
| " accept, state = sample_MH(chain[-1], log_prob, stepsize)\n", | |
| " chain.append(state)\n", | |
| " n_accepted += accept\n", | |
| " \n", | |
| " acceptance_rate = n_accepted / float(n_total)\n", | |
| " \n", | |
| " return chain, acceptance_rate" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Now you're probably excited to see this in action.\n", | |
| "Here we go, taking some informed guesses at the `stepsize` and `n_total` arguments:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 11, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Acceptance rate: 0.720\n", | |
| "Last ten states of chain: 1.05847, 1.59966, 0.14389, -1.13281, 0.24131, -0.77448, -0.59703, 0.67707, 1.47065, 1.27361\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "chain, acceptance_rate = build_MH_chain(np.array([2.0]), 3.0, 10000, log_prob)\n", | |
| "chain = [state for state, in chain]\n", | |
| "print(\"Acceptance rate: {:.3f}\".format(acceptance_rate))\n", | |
| "last_states = \", \".join(\"{:.5f}\".format(state) \n", | |
| " for state in chain[-10:])\n", | |
| "print(\"Last ten states of chain: \" + last_states)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Alright.\n", | |
| "So did this work?\n", | |
| "We achieved an acceptance rate of around 71% and we have a chain of states.\n", | |
| "We should throw away the first few states during which the chain won't have converged to its stationary distribution yet.\n", | |
| "Let's check whether the states we drew are actually normally distributed:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 12, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAFlCAYAAADvZjI4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXiU1cH+8e95JgtLCCCbCGog7JAAsiqETUEEwf1VCnXhZ7V20bZqS+tbl1Yttrb6Wmup1bbuKIorgoJsAZE9ZBERAlFAVJYQ1sDM8zy/PwZSkDWQyZnl/lxXriwzk9xJSLhzznnOMb7vIyIiIhLPHNsBRERERCJNhUdERETingqPiIiIxD0VHhEREYl7KjwiIiIS91R4REREJO4lneB2XbMuIiIiscIc6waN8IiIiEjcU+ERERGRuKfCIyIiInFPhUdERETingqPiIiIxD0VHhEREYl7KjwiIiIS91R4REREJO6p8IiIiETY9u3beeqpp6rt42VkZLBlyxYALrjgguPe9+GHHz7u7cOGDWP79u2UlJTQqVOnSuWYPXs2H3/8ccXrEyZM4Pnnn6/U+6gqxvePu5mydloWERE5TSUlJVx66aUUFhYecVsoFCIp6UQHH1RORkYGS5YsoWHDhie8b1paGrt27Tri7b7v4/s+jhMeGzne53As999/P2lpadx1110nH/70aKdlERERW8aNG0dxcTFdunTh7rvvZvbs2eTk5DBy5Eg6dOhwxOjJo48+yv333w9AcXExQ4cOpVu3buTk5PDZZ58d8f63bt3KkCFD6NixIzfffDOHDmakpaUBsGnTJvr160eXLl3o1KkTubm5jBs3jr1799KlSxdGjx5NSUkJbdu25frrr6dTp06sX7/+sNGiUCjE6NGjad++PVdffTV79uwBDh9RWrJkCQMGDKCkpIQJEybw2GOP0aVLF3Jzc7n//vt59NFHAcjLy6N3795kZ2dzxRVXUFpaCsCAAQP41a9+Rc+ePWnTpg25ublV8j2o2kopIiIS7aaOg68LqvZ9npkFl4w/5s3jx4+nsLCQvLw8IDzVs2zZMgoLC2nRogUlJSXHfOwtt9zChAkTaN26NQsXLuRHP/oRM2fOPOw+DzzwAH379uXee+9lypQpPPvss0e8n5dffpmLL76Ye+65B9d12bNnDzk5OTz55JMVuUpKSli9ejXPPfccvXv3PuJ9rFq1imeffZY+ffowduxYnnrqqWOO3mRkZPDDH/7wsBGejz76qOL266+/nr/+9a/079+fe++9lwceeIDHH38cCBerRYsW8f777/PAAw8wY8aMY359TpYKj4iIiAU9e/akRYsWx73Prl27+Pjjj7nmmmsq3rZv374j7jd37lwmT54MwPDhw6lfv/4R9+nRowdjx44lGAxy+eWX06VLl6N+zHPPPfeoZQfg7LPPpk+fPgCMGTOGJ5544pSmq8rKyti+fTv9+/cH4IYbbjjsc7zyyisB6Nat23HLYGWo8IhI9diyBkrmQukX4IWg1hnQrDuc3QuSa9hOJ4nkOCMx1al27doVLyclJeF5XsXr5eXlAHieR7169SpGYE5Hv379mDt3LlOmTOHGG2/kF7/4Bddff/1xc32XMeaorx+a/2D205GamgpAIBAgFAqd9vsDreERkUjyPCh6CybkwJPd4L2fw4K/wZJ/wUe/g+dHwl/awfT7YM8222lFIqZOnTrs3LnzmLc3adKEb7/9lq1bt7Jv3z7ee+89ANLT02nRogWTJk0CwguJV6xYccTj+/Xrx8svvwzA1KlTK9bDHOqLL76gSZMm/OAHP+Dmm29m2bJlACQnJxMMBk/q8/jyyy9ZsGABEJ4i69u3LxCevlq6dCkAb7zxxgk/77p161K/fv2K9TkvvPBCxWhPpKjwiEhkbFsXLjSTboDgXhj6CNy+HP73W7hnE4z7Eka9Ci36wcdPwJPdIf8126lFIqJBgwb06dOHTp06cffddx9xe3JyMvfeey89e/Zk8ODBtGvXruK2l156iWeffZbOnTvTsWNH3n777SMef9999zF37lw6duzI5MmTOeecc464z+zZs+ncuTNdu3bl1Vdf5Y477gDCa4Sys7MZPXr0CT+Ptm3b8re//Y327dtTWlrKbbfdVvHx77jjDrp3704gEKi4/4gRI3jzzTcrFi0f6rnnnuPuu+8mOzubvLw87r333hN+/NOhy9JFpOoVz4LXbgB8GPIgdB0DTuDY9/+mKDz6s34hdP9/MHQ8JKVUW1wRiRvHvCxdhUdEqlb+a/DmD6FROxj1CtQ/9+Qe54bgowfCoz1tLoH/eQ6SUiObVUTijQqPiFSDwsnwxv+Dc/vAdS9DjfTKv4/Fz8CUO6H1kPD7CCRXfU4RiVfaeFBEIuyLBTD5lvBVV9979dTKDkCPm+HSx2H1h+FpruP/USYiclJ0WbqInL6yDfDa96He2eFprJRjX9Z6UrrfBDu+grl/hMbt4fwfV01OEUlYGuERkdMT2g8TR0OwHEZNhJpHbnh2Sgb+BtpdGr5kfcPSqnmfIpKwVHhE5PTMGQ+b8uCKv0OjtlX3fo2By56EOk3h9Ztg7/aqe98iknBUeETk1K1fDPMegy5joP2Iqn//NevD1f8KT5l9+L9V//5FqpExhjFjxlS8HgqFaNSoEZdeemnF26ZOnUr37t3p0KEDXbt25c477wTCp44bY1izZk3FfR9//HGMMSxZsgQIH0Nx6623kpmZSbdu3RgwYAALFy6sps/u6G688UZef/11qxkO0hoeETk1wXJ481ZIbw5D/xC5j3N2D7jgpzD/cci6BlpGdjdWSQwZ46ZU6fsrGT/8hPepXbs2hYWF7N27l5o1azJ9+nSaNWtWcXthYSE/+clPmDJlCu3atcN1XZ5++umK27Oyspg4cSL/+7/h8j9p0iQ6duxYcfvNN99MixYtWL16NY7jsG7dOj799NMq/Cxjm0Z4ROTUfPwEbCuGkf936ldknawB46B+C3j3jvCuzSIxatiwYUyZEi5br7zyCqNGjaq47Y9//CP33HNPxS7LgUCgYidjgMsvv7xil+Xi4mLq1q1Lw4YNK15fuHAhDz74II4T/q+9RYsWDB9+eBFzXZcbb7yRTp06kZWVxWOPPQbAP//5T3r06EHnzp256qqr2LNnDxAeobntttvo3bs3LVu2ZPbs2YwdO5b27dtz4403VrzftLQ0fv7zn9OxY0cuvPBCNm/efMTnvnTpUvr370+3bt24+OKL2bRpEwBPPPEEHTp0IDs7m+uuu+7Uv7gnoMIjIpVX+gXk/hk6XAaZgyL/8ZJrwojHoXRdeApNJEZdd911TJw4kfLycvLz8+nVq1fFbYWFhXTr1u2Yj01PT+fss8+msLCQiRMncu2111bcVlRURJcuXQ471uFo8vLy2LhxI4WFhRQUFHDTTTcB4dPJFy9ezIoVK2jfvj3PPvtsxWNKS0tZsGABjz32GCNHjuTnP/85RUVFFBQUVBxqunv3brp3705RURH9+/fngQceOOzjBoNBfvrTn/L666+zdOlSxo4dyz333APA+PHjWb58Ofn5+UyYMOEkv5KVp8IjIpX3wW/AOHDxw9X3MVsOgA6Xw/wnwpesi8Sg7OxsSkpKeOWVVxg2bFilH3+wML311ltcccUVlX58y5YtWbt2LT/96U+ZNm0a6enh0dnCwkJycnLIysripZdeoqioqOIxI0aMwBhDVlYWTZo0ISsrC8dx6NixIyUlJQA4jlNRwMaMGcO8efMO+7irVq2isLCQwYMH06VLFx588EE2bNhQ8TUZPXo0L774IklJkVtpozU8IlI5X3wMn70Hg34LdZuf8O7HWytxMuseDnPR/bDqfZj5IFz+VOUeKxIlRo4cyV133cXs2bPZunVrxds7duzI0qVL6dy58zEfe+mll3L33XfTvXv3irJy8LErVqzAdd3jjvLUr1+fFStW8MEHHzBhwgRee+01/vWvf3HjjTfy1ltv0blzZ/7zn/8we/bsisekpoaPeHEcp+Llg6+HQqGjfhxjDt/w2Pd9OnbsWHHS+qGmTJnC3Llzeffdd3nooYcoKCiISPHRCI+InDzfhxkPhC8V7/2j0353GeOmHPPpqM5oAb1uhbyXYdOK0/74IjaMHTuW++67j6ysrMPefvfdd/Pwww/z+eefA+B53hFTPLVq1eKRRx6pmA46KDMzk+7du3Pfffdx8MiokpKSivVCB23ZsgXP87jqqqt48MEHWbZsGQA7d+6kadOmBINBXnrppUp/Tp7nVVyN9fLLL9O3b9/Dbm/bti2bN2+uKDzBYJCioiI8z2P9+vUMHDiQRx55hLKyMnbt2lXpj38yNMIjIifv8w9g/Sdw6WOQUstOhpy7YPmLMPMhGP2anQwip6F58+bcfvvtR7w9Ozubxx9/nFGjRrFnzx6MMYddsn7QsRb2PvPMM9x55520atWKmjVr0rBhQ/70pz8ddp+NGzdy00034XkeAH/4Q/gKy9///vf06tWLRo0a0atXL3bu3Fmpz6l27dosWrSIBx98kMaNG/Pqq68edntKSgqvv/46t99+O2VlZYRCIX72s5/Rpk0bxowZQ1lZGb7vc/vtt1OvXr1KfeyTpcNDReTkeB78IweCe+DHi076UM9Tvfz3uNNdcx+Fmb+HH8yEZsde5Cki1SMtLS1iIzOVdMzDQzXCIyInZ9UU+KYQrnzmiLJT1XuanFCvW2HBkzB7PIyeVL0fW0RiktbwiMiJ+T7k/iW8F07Hyl8ZUuVS64Q3I1z9oc7ZEokCUTK6c1wqPCJyYuvmwFfLoM8dEIiSgeGet4SPnpj3F9tJRCQGqPCIyInl/gXSmkDnUSe+b3VJrQPd/x98NgW2FttOIyJRToVHRI5v49LwCM/5P4bkGrbTHK7nLeH1RJ9oTx4ROT4VHhE5vvlPQI260O0m20mOVKcJZP8PLH8J9myznUZEoliUTMaLSFQq2wgr34XzfxT5A0Ir4dCrwlqbbKanvsijD/2SJ90rKr97s4gkBI3wiMixLf03+F54rUyUWu03Z5bbmRuSPiCFoO04IhKlVHhE5OhC+2Dpf6DN0PCRDlHs3+5QGpkdDHUW244iIlFKhUdEjq7oLdi9GXrdYjvJCeV6WZR4TRidNMN2FBGJUio8InJ0i/4BDVpDiwG2k5yQj8PL7iB6OZ/BN5/ajiMiUUiFR0SO9NXy8OXoPX8ATmz8mpjk9mefnwRL/mU7iohEodj4TSYi1Wv5i5BUA7KvtZ3kpJWSzvteL1gxEfZF/zb3IlK9VHhE5HDBvZA/CdqPhJr1bKeplBdDF8H+nVD4uu0oIhJlVHhE5HAr34V9ZdB1jO0klbbUbwON2odHqEREDqHCIyKHW/4C1DsXMnJsJzkFBrqOhg2LYfPntsOISBRR4RGR/9q2DtbNDY/uxMhi5SNk/Q+YAKx42XYSEYkiMfobTUQiIu9lwECX79lOcurqNIHWg8OLlz3XdhoRiRIqPCIS5rmQ9xK0uhDqNred5vR0+R7s3ARrZ9lOIiJRQoVHRMLWzYUdG6HLaNtJTl+boVCz/oERKxERFR4ROahgEqSmQ9tLbCc5fUmpkHUNrHwP9m63nUZEokCS7QAiUr0yxk054m2p7GdJ6mSmuT24+7czj/q4kvHDIx2tamVfB4uehs/ei8lL7EWkaqnwiAgDnTzqmL287fWxHaXC0YpZpTQ7D+q3gILXVXhERFNaIgKXBeaz2a/LAq+D7ShVxxjodBWsmwO7vrWdRkQsU+ERSXDp7GaQk8d7bm9cArbjVK2sq8H3oOgt20lExDIVHpEEd3FgMakmyNtu9ExnVZnG7aFxR52tJSIqPCKJbqTzMSVeE/L8TNtRIiPrKli/EEq/sJ1ERCzSomWRBNaIUi5winjSvRwwx73vaS8itqXTVfDR76DwDcj5he00ImKJRnhEEtjQwGICxuc993zbUSKnfgY07xEuPCKSsFR4RBLYMGcRq71mrPZj/CiJE+l0NXxTqBPURRKYCo9IgjqDHfR0VjLV62E7SuR1GBl+vvJtuzlExBoVHpEENSSwhIDxmer2sh0l8tLPCk9rrXzXdhIRsUSLlkUS1DBnISVeE1b659iOUqWOtbj6B4HW3JP8MpSWhNf1iEhC0QiPSAJKZxfnO58yzevJia7OihfTDk7drXzPbhARsUKFRyQBDXaWkWxc3nd72o5Sbdb7TeDMLFj5ju0oImKBCo9IAroksJANfkPy/Za2o1Sv9iPDmxDu/Np2EhGpZio8IgkmjT3kOAVMc3uQKNNZFdofvFpLi5dFEo0Kj0iCGeQsJ9WEmJpA01kVGrWFBq01rSWSgFR4RBLMJYFFfOPXY5nf2naU6mdMeE+ekvmwe6vtNCJSjVR4RBJJsJz+Tj4fut3xE/XHv/0I8F1Y9b7tJCJSjRL0N55Iglo3l1pmHzO8braT2NO0C9Q9R+t4RBKMCo9IIln1Prv8GizwOthOYo8x0G4YrJsD+/fYTiMi1USFRyRR+D58Po25Xjb7Sbadxq42QyFUDmtn204iItVEhUckUWzKg52bmOGeZzuJfef2gdR0reMRSSAqPCKJYtVUMA6zvC62k9iXlAKtLoTPPwDPs51GRKqBCo9Iolj1Ppzdi1LSbSeJDm2Hwe5v4atltpOISDVQ4RFJBNvXw9cF0PYS20miR6uLwATCI18iEvdUeEQSwefTws/bDrObI5rUOgPO6a3CI5IgVHhEEsGqqdCgFTRMwN2Vj6ftJfBtEZR+YTuJiERYku0AIlL1MsZNqXg5jT0sS53Nv92h/OGQtwvQ5hL48H/DI2C9brWdRkQiSCM8InEuxykgxbh8pMvRj9SwVfgwUU1ricQ9FR6ROHdRYBmlfhpL/Ta2o0SntkOhZB6U77CdREQiSIVHJI45eAxw8pjldcElYDtOdGpzCXhBKJ5pO4mIRJAKj0gcyzZraWB2MsvVZoPHdHYvSK0La6bbTiIiEaTCIxLHBgbycH3DXC/bdpToFUiCzIGwekb4vDERiUsqPCJxrL+TR57fijLSbEeJbq2HwK6vw5szikhcUuERiVMNKKOLs1bTWSej1UXh56s/tJtDRCJGhUckTvVz8gGY7XW2nCQG1GkCTTvDmhm2k4hIhKjwiMSpgYE8Nvt1KfIzbEeJDa0Gw/pFsLfUdhIRiQAVHpE45ODRz8lnjtcZXz/mJ6f1EPBdKJ5lO4mIRICOlhCJURnHOSbiPLOGema31u98x/G+ZiUPD4Ua9cLTWp2urMZUIlIdVHhE4tCAA5ej53qdbEeJGRm/mcYTye05f/l79PzkksNGxkrGD7eYTESqgsa6ReLQQCePpX4bduhy9EqZ7XamkSmjg9Hp6SLxRoVHJM40YjtZTgmzXV2dVVlzDlzRNsBZYTmJiFQ1FR6RONM/EP7Peran9TuVtZW6rPBaMjCQZzuKiFQxFR6RODPAWcE3fj0+9c+1HSUmzfa60NWspi67bEcRkSqkwiMSRwK45Dj5zHE7A8Z2nJg02+1MwPgVGzeKSHxQ4RGJI13NauqaPczSdNYpW+Fnss1PY0BA63hE4okKj0gc6R/IJ+Q7zNfl6KfMw2Gul00/ZwUGz3YcEakiKjwicSTHyWe534od1LYdJablutk0MjtoZ9bbjiIiVUSFRyRO1GMn2WYduW627Sgxb66XBaB1PCJxRIVHJE70cYpwjE/ugf+s5dRtpj4rvbPJUeERiRsqPCJxIsfJp8yvRb7f0naUuJDrZdPDWUVNym1HEZEqoMIjEhd8cgIFzPc64RKwHSYuzPWySTUhejmf2Y4iIlVAhUckDmSar2hmtmo6qwot9tpS7ieT4xTYjiIiVUCFRyQOHPxPWYWn6uwjhUVeOy1cFokTKjwicaCfk89a70w2+I1tR4krc7xsWjsboWyD7SgicppUeERiXApBejsrmevpcvSqlnvwa1o8y24QETltSbYDiMjp6eZ8Ti2zT9NZEfC535yv/fosefNFfvLaGUe9T8n44dWcSkROhUZ4RGJcjlNA0A/widfBdpQ4ZMh1s+jjFOLomAmRmKbCIxLjcpx8lvmt2U1N21HiUq6XTX2ziyyz1nYUETkNKjwiMewMdpDllDBXx0lEzDyvE55vdHm6SIxT4RGJYX2dQkCXo0fSNtIp9DPICajwiMQyFR6RGJbj5FPqp1Hot7AdJa7lelmcZ1aTxh7bUUTkFKnwiMSs/x4n4elHOaJyvWySjcv5zqe2o4jIKdJvSZEY1dps5ExTylxNZ0XcUq8Nu/1UreMRiWEqPCIxqp+zAoBcLViOuCBJLPA66JgJkRimwiMSo/o5Baz2mrGJBrajJIRcL5sM5xvOMd/YjiIip0CFRyQWBcvp5azU1VnV6ODXWtNaIrFJhUckFn25gBomqPU71Wit35QNfkNNa4nEKBUekVhUPJP9foCFXnvbSRKIYa6bxQVOEUmEbIcRkUpS4RGJRcWzWOK1ZS81bCdJKLleNnXMXjqbYttRRKSSVHhEYs3Ob+CbAnI9XZ1V3T72OuL6Rrsui8QgFR6RWLN2NoDW71hQRhoFfkstXBaJQSo8IrGmeCbUasin/rm2kySkuV4WXcwa0tltO4qIVIIKj0gs8f1w4ckciK8fXyty3WwCxud8p8h2FBGpBP3GFIkl3xTB7m+h5UDbSRLWcr8Vu/wamtYSiTEqPCKxpHhm+HmmCo8toQPHTKjwiMQWFR6RWLJ2FjRqB+ln2U6S0HK9LM51vtUxEyIxRIVHJFYEy+GLjzWdFQUObgmgUR6R2JFkO4CIHFvGuCkVL1/gFPJySjk35aYxa86U4zxKIm2dfyYb/IYqPCIxRCM8IjGin1Og4ySihiHXzeICpxBcHTMhEgtUeERiRF+ngGV+G/boOImokOtlkW72wsaltqOIyElQ4RGJAWewg05OCbmudleOFvO9Tni+CS8kF5Gop8IjEgP6OIUAzPM6WU4iB5WRRr7f4r9bBYhIVFPhEYkBOU4B2/3aFPgtbUeRQ+R62bBhCZSX2Y4iIiegwiMS9Xz6BgqY73XE049sVMl1s8B3YV2u7SgicgL67SkS5TLNV5xltjFPp6NHneV+a0iurWktkRigwiMS5foeWL+Tq8ITdYIkQYscFR6RGKDCIxLlcpx8SrwmbPAb244iR5M5CErXwbZ1tpOIyHGo8IhEsSRC9HZWanQnmmUOCj/X5ekiUU2FRySKdTVrSDPlWr8TzRq0gvTmmtYSiXIqPCJRrG+gANc3LPA62I4ix2IMZA6EtXN1zIRIFFPhEYli/ZwCVviZ7KC27ShyPJmDYF8ZfLXcdhIROQYVHpFotbeUbFMc3txOolvLAYDRtJZIFFPhEYlW63IJGJ9cV8dJRL1aZ8BZXVR4RKKYCo9ItFo7i51+TfL8VraTyMnIHAQbFuuYCZEopcIjEq2KZ/KJ14EQSbaTyMnIHBQ+ZqJknu0kInIUKjwi0WjbOigtIVeno8eO5j11zIRIFFPhEYlGBzax0/47MSQpBTL6qvCIRCkVHpFoVDwL0puz1m9qO4lURuYg2LYWSktsJxGR79DiABHLMsZNOex1B4/lqTOY6vYEjJ1QcmoyB4afF8+C7jfZzSIih9EIj0iUyTZrqWv2aDorFjVsA+nNNK0lEoVUeESiTF+nAM83zPc62o4ilXXwmIl1c8BzbacRkUOo8IhEmZxAAUX+uZSSbjuKnIrMQeG9eHTMhEhUUeERiSK12ct5ZrWOk4hlLQagYyZEoo8Kj0gU6eWsJNm45Gr9Tuyq3QCadlbhEYkyKjwiUSTHKWCvn8JSr43tKHI6Ko6Z2GE7iYgcoMIjEkX6OoUs8tqxn2TbUeR0ZA4CL6RjJkSiiPbhEYkSZ7KV1s5GJgYH2I4ilfDdfZQAUgiSl5pKreKZ0G6YhVQi8l0a4RGJEv0C+YCOk4gH+0nmE6+91vGIRBEVHpEo0c/JZ5N/Bqv8s21HkSqQ62XBtmIo/cJ2FBFBhUckKgRw6esUkutmoeMk4kPFlXYHDoIVEbtUeESiQGdTTD2zmzleZ9tRpIqs8ZtBnbM0rSUSJVR4RKJAv0A+rm+Y53WyHUWqjAlfrbVWx0yIRAMVHpEo0N/JZ4WfSRlptqNIVcocCOXb4as820lEEp4Kj4hlddlFtilmro6TiD8tB4Sfa1pLxDoVHhHL+jqFBIzPXFeFJ+7UbqhjJkSihAqPiGX9nRVs92uzws+0HUUiIXMQbFgE+3baTiKS0FR4RGzyffoF8pnndcIlYDuNRELLgTpmQiQKqPCI2PTtSs40pVq/E8/O6Q1JNTWtJWKZCo+ITWtmAGj9TjxLSoWMvlCsDQhFbFLhEbGp+CNWec35mga2k0gkZQ6Crath+5e2k4gkLJ2WLlINjnaidk3KyUudxxzvYguJpFplDgw/L54F3W6wm0UkQWmER8SSXs5KUk1I63cSQaN2UKcpFH9kO4lIwlLhEbGkv5PPXj+FxV5b21Ek0oyBVhdB8Wxwg7bTiCQkFR4RS/o5+XzitWcfKbajSHVoPRj2lcH6RbaTiCQkFR4RC5qbb8l0Nmk6K5G0HABOEqyZbjuJSEJS4RGxoL+TD6DCk0hq1IVzzofVKjwiNqjwiFjQz8lng9+QYv8s21GkOrW6CL4phB1f2U4iknBUeESqWQpB+jiFzHE7A8Z2HKlOrYeEn2uUR6TaqfCIVLMezmekmXJmel1sR5Hq1rg9pDeH1R/aTiKScFR4RKrZQCePfX4yH3sdbUeR6mYMtL4I1s6B0H7baUQSigqPSDUb6OTxideevdSwHUVsaD0E9u+E9Z/YTiKSUFR4RKrRueZrMp1NzPS62o4itrToD06yprVEqpkKj0g1GujkATBL63cSV2oanHsBrJ5hO4lIQlHhEalGg5zlFHtN+dJvYjuK2NR6CGxeCdvX204ikjBUeESqSU3K6eWs1HSW/PfydO26LFJtVHhEqkkfp4hUE9J0lkDD1lDvHO3HI1KNkmwHEEkUA508dvk1WOy1sx1FbDMGWg9h96IX6DruLfaTfMRdSsYPtxBMJH5phEekWvgMDCxnnpdFUH9nCEDrIdQ2++jhfGY7iUhCUOERqQbtzHrOMvVxabkAABh6SURBVNu0u7L8V0YO+/xkBh24ck9EIkuFR6QaHLwcfbarwiMHpNRigdeBgc5y20lEEoLG1kWqwcDAcgq9DL6lvu0oUo0yxk057u1jAufxYPK/aWm+Yq1/VjWlEklMGuERibS9pXQzn2s6S47wkXseABc5Sy0nEYl/GuERqSLH+mt+hPMxf03xNZ0lR9hEA4q8c7kwsJyn3RG244jENY3wiETY4MBSNvvp5PmtbEeRKDTD60Z3s4r67LAdRSSuqfCIRFAyIQY4K5jpnoenHzc5ihnueQSMX7GwXUQiQ7+BRSKoh/MZ6WYPM7zzbEeRKFXoZ/C1X58LA8tsRxGJayo8IhE02FlKuZ9MrpdlO4pEKR+Hj9zz6O/kk0LQdhyRuKXCIxIxPoMDS8n1sign1XYYiWIzvPNIM+X0dj61HUUkbqnwiERIO7Oe5mYLM7xutqNIlPvY68geP5ULHU1riUSKCo9IhFzkLMXzTcVeKyLHso8Ucr0sLgosA3zbcUTikgqPSIQMDiwlz89kC3VtR5EYMMM7j2ZmK+3Nl7ajiMQlFR6RCGjCNjo7a5nhajpLTs4styueb7TrskiEqPCIRMCFgfCBkB9q/Y6cpC3UJc/PPDCtJSJVTYVHJAIGO0so8Zqwxm9mO4rEkBluNzo7a2lMqe0oInFHhUekitWinAucIqZ73QBjO47EkOkHRgSHBJZYTiISf1R4RKpYPyefVBPS+h2ptNV+M4q9pgxxVHhEqpoKj0gVGxpYxDY/jSV+G9tRJOYYPvB6cL7zKezZZjuMSFxR4RGpQikEGeQs50O3Oy4B23EkBk11e5JsXPh8mu0oInFFhUekCvVxCkk3e5nm9bAdRWJUgd+CjX4DWPmu7SgicUWFR6QKXeIsYodfk4+9TrajSMwyfOD2gDUfwb5dtsOIxA0VHpEqkkSIwYGlfOSdx36SbceRGDbN7QHuPlj9oe0oInFDhUekivRyVlLf7GKa29N2FIlxS/y2ULuRprVEqpAKj0gVucRZxB4/lTletu0oEuM8HGg3PDzCEyy3HUckLqjwiFQFz+XiwBJmeZ0pJ9V2GokH7UfA/l2wdpbtJCJxQYVHpCqsX0QjU8YHms6SqpLRD2rU1bSWSBVR4RGpCivfYZ+fxEyvi+0kEi+SUqDNJbDqfXCDttOIxDwVHpHT5fuw8l1yvSx2Uct2Gokn7UfA3lIomWc7iUjMU+EROV1fLYOy9UzzNJ0lVazVhZCSBkVv2k4iEvNUeEROV9Fb4CQxXYeFSlVLrglth8HKdzStJXKaVHhETofnhf/6zhxEGWm200g86nRleFpr7WzbSURimgqPyOnYsBjK1kOnq2wnkXiVOSh8tVbhG7aTiMQ0FR6R01H4BgRSw9MOIpGQlArtRsBnU7QJochpUOEROVWeG57OajMEaqTbTiPxrNOVsG8HFH9kO4lIzEqyHUAklmSMm1Lx8vlOEa+kfMttK1oydfmU4zxK5DS16A+1GoRHFNsNP+zf4XeVjB9ejcFEYodGeERO0QjnY3b5NZilzQYl0gJJ0H4krJoK+3fbTiMSk1R4RE5BMiGGBRYx3eums7OkenS6CoJ74PMPbCcRiUkqPCKnoK9TQD2zm3fcC2xHkURx7gWQ1gSKJttOIhKTtIZH5BSMCCxgu1+beV6W7SgSp462Tue+pC58b+c00hnJDmpbSCUSuzTCI1JJNdjHEGcJU92eBPU3g1SjyW4OqSbIsMBC21FEYo4Kj0glDXGWkmbKecvtazuKJJgCvwVrvLO4IqDDREUqS4VHpJKuDOSywW/IIr+t7SiScAyT3b70cj6judlsO4xITFHhEamExpSS4+Qz2e2Lrx8fseBttw8Alzsa5RGpDP3GFqmEywLzCRifN90c21EkQW2kEZ947bkykAv4tuOIxAwVHpGT5ftcFchlmdeKdX5T22kkgU12+9LS+Zoupth2FJGYocIjcrK+LqCds57JGt0Ry6a6vSj3k7kikGs7ikjMUOEROVkrJrLfD/Cue77tJJLgdlKL6V43RgYWkEzIdhyRmKDCI3Iy3BAUvMZH3nmUkWY7jQiT3Rzqm10McPJsRxGJCdo1TeQ7jrbD7UBnOf9O2cxk9/sWEokcKdfLYrNfl2sCc5judbcdRyTqaYRH5CRcE5jDNj+N2ToZXaJEiCRed/sxyFlOI0ptxxGJeio8IifQkDIGO0t5w+2noyQkqrzqDiDJeFwTmGs7ikjU029vkRO4KjCXZOMy0R1oO4rIYUr8pnzited/ArP5uzsCH+eoU7IV9x8/vBrTiUQXjfCIHJfPtYFZLPLaUuw3sx1G5AgTQwPJcL6ht7PSdhSRqKbCI3IcvZ2VtHS+ZmJIozsSnaZ6Pdnh1+LawCzbUUSimgqPyHFcF5jJDr8W73u9bEcROap9pPCm24dLnMXUZZftOCJRS4VH5BjqsotLnMVMdvtSTqrtOCLH9Ko7kFQT5PLAfNtRRKKWCo/IMVwZyCXVBJnoDrIdReS4PvUzyPdacF1gJjpQVOToVHhEjspnVGAmeV4mn/nn2A4jckIT3UG0d9ZznlltO4pIVFLhETmKHmYVbZyNvKLRHYkRb7l92OHX4oakD21HEYlKKjwiR3FD0oeU+bV4RweFSozYQw0muf0Z5izUzssiR6HCI/IdTdjGUGcRr7oD2UsN23FETtoL7kUkG5fvBWbajiISdbTTsiSk4+1G+4ukGTj4vOBeVI2JRE5fid+U2W5nRid9xFPuZToKReQQGuEROUQKQUYFZvKR15X1fhPbcUQq7Tl3CI3Ndi52FtuOIhJVVHhEDjHMWUgjs4Pn3IttRxE5JbO9znzhNeZ6LV4WOYwKj8ghbkz6gGKvKfO9jrajiJwSH4fn3cH0dFbRwZTYjiMSNVR4RA7oalbTxSnmOXcIvn40JIZNcvuz10/hpsA021FEooZ+q4sccEvSe2z3a/O62992FJHTsoM0Jrn9uSwwn8a6RF0EUOERASDDbOJiZwkvuhexR5eiSxx4xh1GAI8bkz6wHUUkKqjwiAA3B94nSIDnQlqsLPHhS78J07wejAnMoDZ7bccRsU6FRxJeA8q4OjCXyW4Om6lnO45IlXk6dCnpZg/XBWbZjiJinQqPJLzrk6ZTwwR5xh1mO4pIlVrht2Kh146xSVNJImQ7johVKjyS0Gqwj+8HPmS6ex7FfjPbcUSq3NOh4TQzWxnufGI7iohVKjyS0EYFZnKG2cU/QpfajiISETO9rnzuNeO2pHfB82zHEbFGhUcSVir7+WHSuyxwO7DEb2c7jkhE+Dg8Gbqcds56+Oxd23FErFHhkYR1bWAWTcx2nnCvsB1FJKLe886n2GsKc/6oUR5JWCo8kpBSCHJb0rss8tqywOtgO45IRHkHRnn4phA+n2o7jogVSbYDiERKxrgpx7xtdGAOTc027greCpjqCyViyTveBTxWfxrMeQTaDgOjf/eSWDTCIwknmRC3Jb3DUq81871OtuOIVAuXAPS7CzatgNU6SV0SjwqPJJxrA7NobrbwROhKNLojCSX7WqifAbMe0loeSTgqPJJQalLO7UlvstBrxxwv23YckeoVSIYBvwmP8nz6pu00ItVKhUcSytjANBqb7TwSvA6N7khCyroGmnSCmQ+CG7SdRqTaqPBIwqjHTm5NepcP3W4s89vYjiNih+PAhffCtrWw7DnbaUSqjQqPJIwfJ71Nbcr5Y+ha21FE7Go9BM65ILwvz/7dttOIVAsVHkkIZ7GF6wPTecPtxxq/ue04InYZA4MfgF3fwIKnbKcRqRbah0cSwq+TX8YHHgtdbTuKiDXf3ZtqQnIPcmb+iYFTm7Jo/BhLqUSqh0Z4JO71NCsZEfiECe4INtHAdhyRqPFQ6Hsk4fGr5FdsRxGJOBUeiWsOHvcnP89GvwETQiNsxxGJKuv9JvzTHcZVgXmwfpHtOCIRpcIjcW1UYCYdnC94ODiaclJtxxGJOk+FLuNrvz5M/aU2I5S4psIjcSudXdyZ9BqfeO2Z4vWyHUckKu2hBn8IjoKvlkPeS7bjiESMCo/ErXFJr1CX3TwQvB5tMihybG97feDs3jDjPti91XYckYgwvu8f7/bj3igStUrmwX+G84/QcP4QGm07jUjUa202MCXl17zrnc+dwR8ddlvJ+OGWUolU2jH/utUIj8SfYDm8ewdfeo10GbrISVrtN+fv7kiuCsyjn7PCdhyRKqfCI/Fn7p9g6xp+E7pZC5VFKuGp0GUUe015KOlf1KTcdhyRKqWNByWmfXcjtXbmS95NeYx3vBzmeVmWUonEpn2kMC74Ayal/o47kybxYOj7tiOJVBmN8EjcSGU/jyX/je3U5vdB7RorcioW++14MXQhYwPT6O18ajuOSJVR4ZG4cWfSJNo76/ll8Fa2U8d2HJGY9VBoNCV+E/6c/HfS0eGiEh9UeCQunO8UcXPgfV4IXcQsr6vtOCIxbS81+HnwRzShlAeS/2M7jkiVUOGRmJfOLv6c/HfW+WfykC5BF6kSK/xW/F/oSq4IzIeC123HETltKjwS43z+lPw0jSjjZ8Ef66oskSr0lHsZS73W8N7PYWux7Tgip0VXaUlMuzXwHhcHlvC74Pcp8FvajiMSV1wC3L7/J7xn7uHr/7uMK/b/7og/KrQpocQK7bQssWtdLu5/RjDV68lPgrej4yNEIqOfs4L/JP+RN70+3Bm8jZP9WVMZEgu007LEmR1fwes3sc5vyq+Ct6CyIxI5c73OPB66iqsC8xgd+Mh2HJFTosIjsWf/bnhlFOzfww+DP2M3NW0nEol7f3UvZ5bbmfuTnuN8p8h2HJFKU+GR2OK58MYP4Ot8uPpfrPGb204kkhB8HG4P/pR1/plMSH6MTLPRdiSRSlHhkdjy4W9h1RQYOh7aDrWdRiSh7KQWY4O/ZD/J/Dv5jzSgzHYkkZOmwiOx45O/wyd/g163Qa9bbacRSUgb/EbcvP9OGpkynk15lFo6ZFRihAqPxIZlz8O0cdDuUrj4IdtpRBLaCr8Vtwd/QiezjmeTHyWV/bYjiZyQCo9Ev/xJ8M7t0OoiuPpf4ARsJxJJeNO97vwieBu9nJX8I/kxUgjajiRyXCo8Et2K3oI3b4WMvnDti5CknZRFosU7Xh9+HbqZAYEV/DX5ryQTsh1J5JhUeCR6LX8RXr8JmneHUa9Asi4/F4k2r7oDuTd4AxcHlvBM8qPU1JoeiVIqPBKdFvwN3v4xtBwA338TUuvYTiQix/C8ezF3B2+hr1PACynjSWe37UgiR1DhkejieTDjfvjgN9DhMhg1EVJq204lIicwyR3Aj4N3kG2KeTXl9zRlq+1IIofRWVoSPfbtgsm3hPfZ6XYjDP8LOAEyxk2xnUxETlJfp4Cnkh9nHyk0uvl1OLun7UiSWI55zpAKj1Sb4xWXknFZ4eMivv0ULv5DeJ8dY074OBGJPq3MBp5J/jMZyaUw4gnoMsp2JEkcOjxUotcQZzFMyIHt62H069D7hxVlR0Rizxq/OZfv/x2c3Qve+iG8ewfs32M7liS4JNsBJHGlsp9xSa9wU9IHFOzJ4CfB2/nimXJAIzoisW47dcIXHMx8EOY/Dl8uhGv+DY3b244mCUojPGJFtinm7ZTfclPSBzwTuoSr9j/AF/6ZtmOJSFUKJMPgB2DMZNizBZ4eAPOfAFf79Uj1U+GRalWDffwm6SXeTLmXumY3N+6/mwdD32c/ybajiUiktLoQfjgfMgfB9N/Cs4PhmyLbqSTBaEpLqofvM8z5hN8kv0xzs4WXQhcyPjSKndSynUxEqkOdJnDdy1D0Jrx/N/yjH/S8Bfr/EmrWP+yux73AYfzwSCeVOKWrtCTyNi6FD+6BLxew0juH+4M3sNDXPL5IvDtmOdmzjZcfHsu1gdnsoDaPha7iFfdCgifxN7gKj5yALksXCzYsgTmPwOoPoXYjfrX9cia5/fE0kyoiQHvzBb9NeoELAp+ywW/I30MjmeT2P+4UtwqPnIAKj1SPjHHv0dcp5JbAe/QLFLDNT+OZ0HCedwezS9NXInIEn35OPnckTaabs5pN/hk8ExrGJLc/Ozhyl3UVHjkBFR6JsPIyyHuF4vcfJ9PZxBY/nX+GhvOCO5g91LCdTkSink8fp5Dbk96kl/MZe/xU3nL78Lw7hM/8cyrupcIjJ6DCIxEQ2g9rZkDBJFg1FUJ7We614rnQEN73eunKKxE5JR3NOq4PTOeywHxqmCD5Xgvedvvwrns+i8aPsR1PopsKj1SR8jJu+92fuTCwnAudZdQ3u9jmp/Geez6T3P4U+C1tJxSROFGPnVwZmMflgXlkO+vwfIPToi+0GQptLoYGrbQru3yXCo9UzsHLQlPZT5ZZS09nFX2cQno6n5FsXMr8WszyuvC224dcL4uQdjgQkQhqab7issB87jhrVfjMPYAzWkKrwZDRB865ANIa2Q0p0UCFJ55V2Z4Vngvb1sHX+Tw18U26O6vobNaSaoIArPKaM8vrykduV5b5rXEJnG50EZFKKRk/HEq/CF/9+fkHUDIPQnvDNzZoDef0hqadw09NOkLKkQufJa6p8MSzyp4mXpNyzjXfMu365uGCs60Yvi4M/9UUDB/wF/QDFPkZLPLascRrw2KvLaWkRyK+iMgpSyZElllLD2cVv+6wDTYsgr2l4RuNE572atw+/LxBKzgjM/y81hkV02Ha6DCuqPDEs4xxUzB41GEv6WY3ddlNQ7ODxqaUxmyniSmlsdlOY1NKc7OFxmb7YY8v9dP4zDuHlX746VMvg9V+My06FpGYUjJ+OPg+lG2Ar/NhUz5sWgFbPoftX4B3yBleKWmQfhakn8Wkzz02cQbf+Gew2a9LqV+HUtIo89MoJe2YU/YqQ1EpAQvPwc+r4vM79PVTve3A6ye6zffBC4Z/uLxQeKrIPfT17zy5h768D4J7v/O0B0Ll4ecH37Z/F+zdDuXb2b5tM3XYQ8Ac/du13a/Nt349vvXr8ZXfkBK/CV/6TfjCb8IXfmN2kHYKX2ARkdhR8tAQ2P4lbF0Tftq+HnZ+BTu+4qv1xTSh9Ji/Q3f6Ndnup7GDWuymBnv9VHZTg2HntYKUWuFps5Q0SK4FSakQSDnwlHzI68mHvP3Ak5MUHoUyBpzAgZcdMIe87ATCtx/17Q6H/f9+cAF3Yi/kjtLC83UhPDuEUyoVx7otHgVSILlm+Ifp0Oc16kHNejy3fDtl1KbMr82OA8+3+HX5lnps9uuxjxTbn4GISFQL4NKQMhqaMuqbXdRjF/XMLuqzM/y62UU6u6nFPmqZcmqxjzb1Hdi/O/zk7rP9KZzAUYrRoW8/2bcd9vYTvO27ut0AQ/9wUmlPw6kVHmPMNKBhJBLFqYbAFtsh5Aj6vkQffU+ik74v0Uffk8rZ4vv+0KPdcKIRHqkEY8wS3/e7284hh9P3JfroexKd9H2JPvqeVB2d4igiIiJxT4VHRERE4p4KT9V62nYAOSp9X6KPvifRSd+X6KPvSRXRGh4RERGJexrhERERkbinwhMhxpg7jTG+MUaX9VtmjPmTMeYzY0y+MeZNY0w925kSmTFmqDFmlTFmjTFmnO08ic4Yc7YxZpYx5lNjTJEx5g7bmSTMGBMwxiw3xrxnO0s8UOGJAGPM2cAQ4EvbWQSA6UAn3/ezgc+BX1vOk7CMMQHgb8AlQAdglDGmg91UCS8E3On7fgegN/BjfU+ixh3AStsh4oUKT2Q8BvySuN36Obb4vv+h7/sHD9H5BGhuM0+C6wms8X1/re/7+4GJwGWWMyU03/c3+b6/7MDLOwn/B9vMbioxxjQHhgPP2M4SL1R4qpgx5jJgo+/7K2xnkaMaC0y1HSKBNQPWH/L6BvSfa9QwxmQAXYGFdpMI8DjhP5w920HixdGPgJXjMsbMAM48yk33AL8hPJ0l1eh43xPf998+cJ97CA/fv1Sd2URigTEmDXgD+Jnv+zts50lkxphLgW99319qjBlgO0+8UOE5Bb7vX3S0txtjsoAWwAoTPkitObDMGNPT9/2vqzFiwjnW9+QgY8yNwKXAhb72YrBpI3D2Ia83P/A2scgYk0y47Lzk+/5k23mEPsBIY8wwoAaQbox50ff9MZZzxTTtwxNBxpgSoLvv+zr4zSJjzFDgL0B/3/c3286TyIwxSYQXjl9IuOgsBr7n+36R1WAJzIT/OnsO2Ob7/s9s55HDHRjhucv3/UttZ4l1WsMjieBJoA4w3RiTZ4yZYDtQojqwePwnwAeEF8e+prJjXR/g+8CgAz8feQdGFkTiikZ4REREJO5phEdERETingqPiIiIxD0VHhEREYl7KjwiIiIS91R4REREJO6p8IiIiEjcU+ERERGRuKfCIyIiInHv/wMJJTaRBaQ/zwAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<Figure size 720x432 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "def plot_samples(chain, log_prob, ax, orientation='vertical', normalize=True,\n", | |
| " xlims=(-5, 5), legend=True):\n", | |
| " from scipy.integrate import quad\n", | |
| " \n", | |
| " ax.hist(chain, bins=50, density=True, label=\"MCMC samples\",\n", | |
| " orientation=orientation)\n", | |
| " # we numerically calculate the normalization constant of our PDF\n", | |
| " if normalize:\n", | |
| " Z, _ = quad(lambda x: np.exp(log_prob(x)), -np.inf, np.inf)\n", | |
| " else:\n", | |
| " Z = 1.0\n", | |
| " xses = np.linspace(xlims[0], xlims[1], 1000)\n", | |
| " yses = [np.exp(log_prob(x)) / Z for x in xses]\n", | |
| " if orientation == 'horizontal':\n", | |
| " (yses, xses) = (xses, yses)\n", | |
| " ax.plot(xses, yses, label=\"true distribution\")\n", | |
| " if legend:\n", | |
| " ax.legend(frameon=False)\n", | |
| " \n", | |
| "fig, ax = plt.subplots()\n", | |
| "plot_samples(chain[500:], log_prob, ax)\n", | |
| "despine(ax)\n", | |
| "ax.set_yticks(())\n", | |
| "plt.show()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Looks great!" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Now, what's up with the parameters `stepsize` and `n_total`?\n", | |
| "We'll discuss the step size first: it determines how big the random steps are the Markov chain takes.\n", | |
| "If the step size is too large, we will jump around a lot in the tails of the distribution, where probability is low.\n", | |
| "The Metropolis-Hastings sampler will reject most of these moves, meaning that the acceptance rate decreases and convergence is much slower.\n", | |
| "See for yourself:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 13, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Acceptance rate: 0.104\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAFlCAYAAADvZjI4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xUZaLG8d+ZmSSU0KRLL9JCIEjvKCpIbxYWVGSx91Wv7NULunrtu3pddbGwu9hABRQkoGJBQHoJEJBOBAtK7ySZOef+cSArQgIJybwzZ57v55NPMjnJzBOGJE/e877vsRzHQURERMTLfKYDiIiIiBQ1FR4RERHxPBUeERER8TwVHhEREfE8FR4RERHxPBUeERER8bzAWY5rzbqIiIhECyu3AxrhEREREc9T4RERERHPU+ERERERz1PhEREREc9T4RERERHPU+ERERERz1PhEREREc9T4RERERHPU+EREREpYvv37+fVV18N2+PVrl2b3bt3A9ChQ4c8P/bJJ5/M83ivXr3Yv38/GRkZNG3aNF855syZw4IFC3Jujxs3jrfeeitf91FYLMfJczNl7bQsIiJynjIyMujTpw/p6emnHQsGgwQCZ7vwQf7Url2bZcuWUaFChbN+bGJiIocPHz7t/Y7j4DgOPp87NpLX15CbRx99lMTERB544IFzD39+tNOyiIiIKaNHj2bLli2kpKTw4IMPMmfOHDp37ky/fv1o0qTJaaMnzz//PI8++igAW7ZsoWfPnrRs2ZLOnTuzfv360+5/z549XHHFFSQlJTFq1Ch+O5iRmJgIwM8//0yXLl1ISUmhadOmzJs3j9GjR3Ps2DFSUlIYNmwYGRkZNGzYkOuvv56mTZuyY8eOU0aLgsEgw4YNo3HjxgwZMoSjR48Cp44oLVu2jG7dupGRkcG4ceN44YUXSElJYd68eTz66KM8//zzAKSlpdGuXTuaNWvGwIED2bdvHwDdunXjoYceok2bNjRo0IB58+YVynNQuJVSREQk0s0aDTvXFO59VkmGK5/O9fDTTz9Neno6aWlpgHuqZ8WKFaSnp1OnTh0yMjJy/dybb76ZcePGcdFFF7F48WJuv/12vvrqq1M+5rHHHqNTp06MGTOG1NRUxo8ff9r9vPfee/To0YOHH36YUCjE0aNH6dy5My+//HJOroyMDDZt2sSECRNo167dafexYcMGxo8fT8eOHRk5ciSvvvpqrqM3tWvX5tZbbz1lhOfLL7/MOX799dfz97//na5duzJmzBgee+wxXnzxRcAtVkuWLGHmzJk89thjfPHFF7n++5wrFR4RERED2rRpQ506dfL8mMOHD7NgwQKuuuqqnPdlZmae9nFz585l6tSpAPTu3Zty5cqd9jGtW7dm5MiRZGdnM2DAAFJSUs74mLVq1Tpj2QGoUaMGHTt2BGD48OG89NJLBTpddeDAAfbv30/Xrl0BuOGGG075GgcNGgRAy5Yt8yyD+aHCIyIisSWPkZhwKlmyZM7bgUAA27Zzbh8/fhwA27YpW7ZszgjM+ejSpQtz584lNTWVESNG8Kc//Ynrr78+z1y/Z1nWGW//Nv/J7OcjISEBAL/fTzAYPO/7A83hEZF8qj06Nc8XETldqVKlOHToUK7HK1euzK+//sqePXvIzMxkxowZAJQuXZo6derw4YcfAu5E4lWrVp32+V26dOG9994DYNasWTnzYX7r+++/p3Llytx0002MGjWKFStWABAXF0d2dvY5fR3bt29n4cKFgHuKrFOnToB7+mr58uUATJky5axfd5kyZShXrlzO/Jy33347Z7SnqKjwiIiIFLHy5cvTsWNHmjZtyoMPPnja8bi4OMaMGUObNm24/PLLadSoUc6xd999l/Hjx9O8eXOSkpKYNm3aaZ8/duxY5s6dS1JSElOnTqVmzZqnfcycOXNo3rw5LVq04P333+eee+4B3DlCzZo1Y9iwYWf9Oho2bMgrr7xC48aN2bdvH7fddlvO499zzz20atUKv9+f8/F9+/blo48+ypm0/FsTJkzgwQcfpFmzZqSlpTFmzJizPv750LJ0EcmXs43iZDzdO0xJREROo2XpIiIiErtUeERERMTzVHhERETE81R4RERExPNUeERERMTzVHhERETE81R4REREwsCyLIYPH55zOxgMUrFiRfr06ZPzvlmzZtGqVSuaNGlCixYtuP/++wH3quOWZbF58+acj33xxRexLItly5YB7mUobrnlFurVq0fLli3p1q0bixcvDtNXd2YjRoxg8uTJRjOcpEtLiIhIzCnsXcHPZf+pkiVLkp6ezrFjxyhevDizZ8+mWrVqOcfT09O58847SU1NpVGjRoRCIV5//fWc48nJyUyaNIlHHnkEgA8//JCkpKSc46NGjaJOnTps2rQJn8/Htm3bWLduXSF+ldFNIzwiIiJh0qtXL1JT3bI1ceJEhg4dmnPs2Wef5eGHH87ZZdnv9+fsZAwwYMCAnF2Wt2zZQpkyZahQoULO7cWLF/PEE0/g87m/2uvUqUPv3qcWsVAoxIgRI2jatCnJycm88MILALzxxhu0bt2a5s2bM3jwYI4ePQq4IzS33XYb7dq1o27dusyZM4eRI0fSuHFjRowYkXO/iYmJ3HfffSQlJdG9e3d27dp12te+fPlyunbtSsuWLenRowc///wzAC+99BJNmjShWbNmXHvttQX/xz0LFR4REZEwufbaa5k0aRLHjx9n9erVtG3bNudYeno6LVu2zPVzS5cuTY0aNUhPT2fSpElcc801OcfWrl1LSkrKKZd1OJO0tDR+/PFH0tPTWbNmDTfeeCPgXp186dKlrFq1isaNGzN+/Picz9m3bx8LFy7khRdeoF+/ftx3332sXbuWNWvW5FzU9MiRI7Rq1Yq1a9fStWtXHnvssVMeNzs7m7vuuovJkyezfPlyRo4cycMPPwzA008/zcqVK1m9ejXjxo07x3/J/FPhERERCZNmzZqRkZHBxIkT6dWrV74//2Rh+vjjjxk4cGC+P79u3bps3bqVu+66i08//ZTSpUsDbtnq3LkzycnJvPvuu6xduzbnc/r27YtlWSQnJ1O5cmWSk5Px+XwkJSWRkZEBgM/nyylgw4cPZ/78+ac87oYNG0hPT+fyyy8nJSWFJ554gh9++CHn32TYsGG88847BAJFN9NGhUdERCSM+vXrxwMPPHDK6SyApKSknCuO56ZPnz68/fbb1KxZM6esnPzcVatWEQqF8vz8cuXKsWrVKrp168a4ceMYNWoU4J66evnll1mzZg1jx47l+PHjOZ+TkJAAuKXm5NsnbweDwTM+jmWdekkrx3FISkoiLS2NtLQ01qxZw+effw5Aamoqd9xxBytWrKB169a53uf5UuEREREJo5EjRzJ27FiSk5NPef+DDz7Ik08+ycaNGwGwbfu0UzwlSpTgmWeeyTkddFK9evVo1aoVY8eO5eRFwTMyMnLmC520e/dubNtm8ODBPPHEE6xYsQKAQ4cOUbVqVbKzs3n33Xfz/TXZtp2zGuu9996jU6dOpxxv2LAhu3btYuHChYB7imvt2rXYts2OHTu45JJLeOaZZzhw4ACHDx/O9+OfC63SEhERCaPq1atz9913n/b+Zs2a8eKLLzJ06FCOHj2KZVmnLFk/KbeJvW+++Sb3338/9evXp3jx4lSoUIHnnnvulI/58ccfufHGG7FtG4CnnnoKgMcff5y2bdtSsWJF2rZty6FDh/L1NZUsWZIlS5bwxBNPUKlSJd5///1TjsfHxzN58mTuvvtuDhw4QDAY5N5776VBgwYMHz6cAwcO4DgOd999N2XLls3XY58r62QTzEWeB0Uk9pxtOe+5LM8VEW9JTEwsspGZfLJyO6BTWiIiIuJ5KjwiIiJyXiJkdCdPKjwiIiLieSo8IiIi4nkqPCIiIuJ5KjwiIiLieSo8IiIi4nnaeFBEClVe+/Rojx4RMUUjPCIiIuJ5KjwiIiLieSo8IiIi4nkqPCIiIuJ5KjwiIiLieSo8IiIi4nkqPCIiIuJ5KjwiIiLieSo8IiIi4nkqPCIiIuJ5KjwiIiLieSo8IiIi4nkqPCIiIuJ5KjwiIiLieSo8IiIi4nkqPCIiIuJ5KjwiIiLieSo8IiIi4nkqPCIiIuJ5AdMBRERMqz06NddjGU/3DmMSESkqGuERERERz1PhEREREc9T4RERERHPU+ERERERz1PhEREREc9T4RERERHP07J0EYkYWh4uIkVFIzwiIiLieSo8IiIi4nkqPCIiIuJ5KjwiIiLieSo8IiIi4nkqPCIiIuJ5KjwiIiLieSo8IiIi4nkqPCIiIuJ5KjwiIiLieSo8IiIi4nkqPCIiIuJ5KjwiIiLieSo8IhJeoSBkHQHHMZ1ERGJIwHQAEfG2iuyjj38Rl/jS4Nm74Oge90CgGFS4CGp1hKaDoXprs0FFxNNUeESkSFRhD/cEpjLYP5d4K8QGuzo07AVla4I/Ho7sgp1rYPkEWDwOqjSjm68Hc+wUwDIdX0Q8RoVHRAqZwwj/ZzwQ+IA4QkwMXcqEUA+2OheS0b/36R+eeQjSp8D8F/h3/HN8FmrFI9kj2UXZ8EcXEc/SHB4RKTSJHOUfcS/yaNxbLLUb0j3rOcYGb2Src2Hun5RQClqOgDuW8mT2ULr5VjEzYTStrfVhyy0i3qfCIyKFojJ7+Sh+LJf7lvN49jBuzP4vfnAqnfsdBOJ5PdSX3ln/y0GnJO/F/y/9ffOLLrCIxBQVHhE5b9WtXXwQ/xeqWHu5LvvPjA/1pqDzcDY71RmQ9ThL7Ea8EPcPrvF/XbhhRSQmqfCIyHmpxD4mxT9OGesIw7L+m4V20nnf5yFKMDL7Qb6xm/FM3Btc5Z9z/kFFJKZp0rKIFFhJjvHv+Gcpy2Guyfof1jp1Cu2+M4nnluw/8SbP81TgTdjcA+pfluvH1x6dmuuxjKfPMFlaRGKKRnhEpEB82Lwa9380sHZwe/Y9hVp2Tsoijtuz72GjUwM+uAF2phf6Y4hIbFDhEZECuTcwma7+1TwSHMlcu3mRPc5hSnBj1oPuaq4ProPjB4rssUTEu1R4RCTfuvnSuDvwMe8HuzEpdGmRP94vXABD/gX7vodpd+qyFCKSbyo8IpIvFdnH3+Je5Tu7JmOCI8L3wLXaw2Vj4bvpsPTN8D2uiHiCCo+InDvH4am4NylBJndk300m8eF9/A53uxOXZ4+BPVvC+9giEtW0SkskBuW1ognyWNWU9h6X+Vfyl+zr8t49uahYFvR9CV5tDx/fDjfOBJ8//DlEJOpohEdEzs2hnfDpaBbbjfhXqIe5HGWqwZXPwI5FsOR1czlEJKqo8IjIufn8EQhm8l/ZN+OY/tHR/Fr31NZX/wsHfzabRUSiggqPiJzdtnmw5kPodC/fO1VMp3FPbV35LIQyYfb/mE4jIlFAhUdE8hbKhpkPQNla0Ok+02n+o3w96HivW8S2zTWdRkQinAqPiORt8TjYtd4dUYkrbjrNqTr/CcrWhNQH8BMynUZEIpgKj4jk7uhemPsc1L8cGvY0neZ0ccWhx1OwewNX6wKjIpIHFR4Ryd28v8Lxg3D5Y6aT5K5Rb6jRlvsCUyjOcdNpRCRCqfCIyJnt3+4u+075A1ROMp0md5YFlz9OJWs/f/TPMp1GRCKUCo+InNnXT4Llg0v+23SSs6vZls9CrbglMIMLOGg6jYhEIBUeETndr9/BqknQ9hYoU910mnPybPAaipPJHYFppqOISARS4RGR0819DuJLusu+o8QWpxof250Y5v+Ciuw3HUdEIowKj4icop71I6RPhTY3QYkLTMfJl5eD/YkjyM2BGaajiEiE0cVDReQUdwSmucu9299pOkq+ZThV+djuyHD/F7wW7Mtuypz3fRb4QqsiElE0wiMiOWpZO+nv+xZa/xFKVjAdp0BeDg4knmxGBfIuKiISW1R4RCTHHf5pZBOADnebjlJg25yqTLc7cL1/NuW0YktETlDhEREAqlu7GOSfx3uh7pBYyXSc8/JqsD8lrEyu839hOoqIRAgVHhEBYKR/Fg4Wrwejf07KJqc6X4VSuD7wOQlkmY4jIhFAhUdEKM1hrvF/zSd2e3ZS3nScQvF6qA8VrIMM9s8zHUVEIoAKj4jwB/9XlLQyeTPYy3SUQrPIbswquy6j/Klg60rqIrFOhUckxsURZETgM+aHkljn1DYdpxBZvB7sQ13fTtgw03QYETFMhUckxvXxLaSKtY83Q9E/d+f3PrVbs92uCN++ZDqKiBimwiMS0xxuCsxko12NOXZz02EKXQg/40O94IclsGOp6TgiYpAKj0gMa+9bRxPf97wZ6gVYpuMUicmhLhBfCpa+YTqKiBikwiMSw67zz2afk8i0UEfTUYrMEYpDyh/c64Md/tV0HBExRIVHJEZVYQ9X+JbxfqgbmcSbjlO02twEdjYsn2A6iYgYosIjEqOGBr7Gh8O7oe6moxS9ChdBvUth2T8hlG06jYgYoMIjEoPiCDLU/xVz7ObscCqbjhMebW6GQz/Bel1UVCQWqfCIxKArfMuoZO3nrdDlpqOEz0VXQNlasESTl0VikQqPSAy6PvA52+2KzPXgUvRc+fzQehR8Px9+WWc6jYiEmQqPSKz5ZR1tfet5J3QZdqz9CEgZBv54WKHJyyKxJmA6gIiE2bLxZDpxfBjqGvaHrj3a8PyZkuWhcT9YNREuexTiipvNIyJhE2N/3onEuKwjsOp9Ztht2Udp02nMaHkDHD8A66aZTiIiYaTCIxJL1k2DrENMCl5qOok5tTvDBXW1J49IjNEpLZEoldfpoYync7kQ6Iq3oXx9lv7YsIhSRQHLgotvgC/Gwq4NUDGG/y1EYohGeERixe5NsH0BtLgOr14365ylDANfHKx4y3QSEQkTFR6RWLHybbD80Hyo6STmJVaERr0g7T3IPm46jYiEgQqPSCwIZbu/3BteCaViZGfls2k5Ao7thfUzTCcRkTBQ4RGJBRs/gyO7TpzOEgDqdHN3XtaePCIxQYVHJBaseAtKVYX6l5lOEjl8Pncuz7a5VGOX6TQiUsRUeES87uBPsHk2pPwB/FqYeYrm1wIwyD/PcBARKWoqPCJel/YuODa0GG46SeQpVwtqd2awfx7gmE4jIkVIhUfEyxwHVr77n8325HQpw6jt+4VW1gbTSUSkCKnwiHjZjsWwb5s7V0XOrEk/DjvFGOKfazqJiBQhFR4RL1s1CeJKQOO+ppNErviSzAq1obd/McXRnjwiXqUZjCJelX0c1k51y05Couk0Rp3tKu1tra5cFZhLD98yPrY7hSmViISTRnhEvGrjp+5VwU+sRJLcLXEast2uyBD/N6ajiEgRUeER8arV77t779TpajpJxHPwMSXUhQ6+dVzIbtNxRKQIqPCIeNGR3bDpc0i+Cnx+02miwhS7Mz7LYaB/vukoIlIEVHhEvCh9KthBnc7Khx+cSiyyGzPYPxftySPiPSo8Il60aiJUSYbKSaaTRJWpoU7U9e2kubXFdBQRKWQqPCIeU8/6EX5aAc00upNfn4bakOkEGOD/1nQUESlkKjwiHjPQPx8snzt/R/LlICX50r6YPv6F+AmZjiMihUj78Ih4iIXtFp563aFUZdNxCtXZ9tIpLNNCHenlX0InXzrf2M3D8pgiUvQ0wiPiIe1831HN2qPJyufhazuFA04J+uu0loinqPCIeMhA33wOOcWhYS/TUaJWFnGkhtrSw7dUl5oQ8RAVHhGPSCCLK/1LmBVqA/ElTMeJatNCnShpZXK5b7npKCJSSFR4RDziUt9KSlnHmGZ3MB0l6i1xGvKjU16bEIp4iAqPiEf08y9gl1OGhbb23jlfDj6mhzrQ2beG8hwwHUdECoEKj4gHlOIol/rSmBFqh61v60LxUagTAcumt3+R6SgiUgj0k1HEA67wLSPBymZ6SKezCstGpwbf2TUZqNVaIp6gwiPiAf3937LdrshKp77pKJ7ycagjLXybqWXtNB1FRM6TCo9IlCvPATr41vKJ3R6wTMfxlOmhDtiOxQCfRnlEop0Kj0iU6+VfTMCydTqrCPxMeRbbjU9sQqgrqItEMxUekSjXz7+ADXZ1Njg1TUfxpGl2B+r6dpJkZZiOIiLnQYVHJIpdyG5a+zYyTaM7RWZWqA3Zjp9+/gWmo4jIeVDhEYliff0LAU7M35GicIBE5trN6ONfhIVtOo6IFJAKj0gU6+dfwEq7Pjscb10ZPdJMD7WnmrWHltZG01FEpIBUeESiVD3rR5J83zM9pNGdovaF3ZLjTlzOiJqIRB8VHpEo1c+/gJBjMSPUznQUzztCcb6wL6aXfzF+QqbjiEgBqPCIRCPHoZ9vAYvsJuyinOk0MeGTUHsqWgdp71tnOoqIFIAKj0g0+mkldXy/6MroYTTHTuGQU5x+Pq3WEolGKjwi0Sh9ClmOn09DrU0niRmZxPO53Yqe/qXEk206jojkU8B0ABEvqz06Nc/jGU/3zv+d2jakT+UbO4WDJBYwmRTE9FAHBvvn0dW3itl2K9NxRCQfNMIjEm22L4BDP2l1lgHf2knsdRK1WkskCqnwiESbNZMhrgRf2BebThJzggSYGWrLZb4VFOe46Tgikg86pSUSTYJZsO5jaNiLY8uKFdnDnO1UXCz7JNSB4YEvucy3gk80aVwkamiERySabP0aju2D5KtMJ4lZS5yG7HTK0U+ntUSiigqPSDRZMxmKlYV6l5pOErMcfMwItaOrL43SHDYdR0TOkQqPSLTIOgrrU6FJfwjEm04T06aHOhBvhejhX2Y6ioicIxUekWix8VPIPgLJQ0wniXmrnbpk2JXp69NpLZFoocIjEi3Sp0BiFajV0XQSweITuz0dfelw+FfTYUTkHKjwiESDY/th0+fQdBD4/KbTCO61tfyWA+ummY4iIudAhUckGqyfAaEsaKrTWZFio1OD9XYNdyK5iEQ87cMjEg3WTIZydaDauW02qH10wuOTUHsa7fgA9u+AsjVMxxGRPGiERyTSHd4F276BpoPBskynkd/4xD5xeY+1U80GEZGzUuERiXTrPgbH1uqsCLTdqQwXXuxOKBeRiKbCIxLp1kyGSk2gUmPTSeRMkofAz6tg92bTSUQkDyo8IpFs/w7Yscg9nSWRKWkgYGmURyTCqfCIRLKTv0RVeCJX6QvdvZHSJ4PjmE4jIrlQ4RGJZOlToForuKCO6SSSl6aDYPdG2LnGdBIRyYUKj0ik2r0Jdq7W6E40aDIAfAF3lEdEIpIKj0ikWjMZsE7MEZGIVrK8ewX7NVPAtk2nEZEzUOERiUSO444W1O4EpauaTiPnIvkqOPgD7FhsOomInIEKj0gk2rka9mzW3jvRpGEvCBSHNR+aTiIiZ6DCIxKJ1kx254Q07mc6iZyrhERoeKW7UWQo23QaEfkdFR6RSGPbkD4V6nWHEheYTiP5kTwEju6BrXNMJxGR31HhEYk0Oxa7c0F0Oiv61L8MipXRFdRFIpAKj0ikSZ8CgWLu6RGJLoEE9zTk+hmQddR0GhH5DRUekUgSCrpzQBr0hIRSptNIQSQPgazDsOkz00lE5DdUeEQiybZv4Mgunc6KZrU7Q2JlndYSiTAqPCKRJH0KJJSG+pebTiIF5fND0iDY9Dkc2286jYicoMIjEimCmfDdJ9CoD8QVM51GzkfyVRDKcufyiEhECJgOICInbJoNmQch2b12Vu3RqYYDSYFVuxjK1XE3IWwx3HQaEUEjPCKRI30ylKgAdbqZTiLny7LceVjb5sKhX0ynERFUeEQiQ+Zh2PApJA0AvwZePaHpEHBsWPuR6SQiggqPSGRYnwrBY9B0sOkkUlgqNYLKye7InYgYp8IjEglWvw9la0KNdqaTSGFKHgw/LIW920wnEYl5Kjwiph36BbZ+DclXg0/fkp5ycsQufYrZHCKiwiNiXPoUd65Hs6tNJ5HCdnLUTpsQihinwiNi2poPoGoKVGxoOokUheQhsOs7+GWt6SQiMU2FR8SkXRvhp5Ua3fGyJgPA8rt78oiIMSo8Iiat+QAsn1ZneVliRah3Kaz+EGzbdBqRmKXCI2KM467OqtsNSlUxHUaKUvNr4eAP8P1800lEYpYKj4ghLa2NsH87NLvGdBQpag17QXwpWPW+6SQiMUuFR8SQgf75EFfCvVioeFt8CWjSH9Z9DFlHTacRiUkqPCIGxBGkj38RNOoNCYmm40g4NL8Gsg7Dhpmmk4jEJBUeEQO6+dIoax3R6axYUqsTlK4OqyaZTiISk1R4RAwY4J/Pbqc01L3EdBQJF5/P3X5gy5e6grqIASo8ImFWiqNc5lvJJ6H2ujJ6rGl+rburti4oKhJ2KjwiYdbTv4QEK5uPQp1MR5Fwq9gQLmwBqyaaTiISc1R4RMJskG8+W+yqrHbqmo4iJjQfCjvX6FITImGmwiMSRtWtX2nvX3didMcyHUdMaDoYfAFNXhYJMxUekTAa7JuH7VhMDXU2HUVMKVkB6l/uXlvLDplOIxIzVHhEwsTCZrB/Lt/aSfxEBdNxxKTm18Chn2HbN6aTiMQMFR6RMGnrW09N3y4+DHU1HUVMa3AlJJTRpSZEwkiFRyRMhvjnctApzud2K9NRxLS4YpA0AL6bDpmHTacRiQnaBEQkDEpwnCt9i5kW6sBxEnLeX3t0qsFUYlTzobBiAqybBi2GmU4j4nka4REJg97+RZS0Mpms01lyUs12UL4+rHzHdBKRmKDCIxIGQ/xz2WJXZYVzkekoEiksC1oMh+0LYPdm02lEPE+FR6SI1bR+oa1vPVNCXdDeO3KK5kPB8sPKt00nEfE8FR6RIjbYP/fE3ju6lIT8TqkqcNEV7qUmQkHTaUQ8TYVHpAi5e+/MY77dlJ2UNx1HIlGL4XD4F9g823QSEU9T4REpQu1966hu7dbeO5K7Bj2gZCVNXhYpYlqWLlKErvV/zQGnhPbe8bC8thbIeLr32e/AHwfNr4VFr8LhXyGxUiGmE5GTNMIjUlSO7KGHbylTQ53JJN50GolkLa4DO6gLiooUIY3wiBSVVRNJsIJMDF1qOolEuooNoEZbd7VWh7vcJetncd4jSyIxRiM8IkXBcWDFBJbbF7HRqWE6jUSDFtfB7tEcUCkAABOeSURBVI2wY4npJCKepMIjUhS2L4TdG5kUusR0EokWSQMgriSsfMt0EhFPUuERKQrLJ0BCaWaE2plOItEioRQ0HQTpU+H4AdNpRDxHhUeksB3bB+s+huSrOEYx02kkmrQaCdlHYfUHppOIeI4Kj0hhW/0BBI9DyxtMJ5FoU+1iqJoCy/7pzgMTkUKjVVoihclx3NNZF7aAqs2BH0wnEoPyWkkFuaymav1HmH4XbF8EtdoXUTKR2KMRHpHC9MMy+HUtXKzRHSmgpoMhobQ7yiMihUaFR6QwLRsP8YmQPMR0EolW8SXdnZfXfQxH9phOI+IZKjwiheXIbneFTfOh7oobkYJqNRJCWZCm62uJFBYVHpHCsuItCGVC61Gmk0i0q9QYanaAZf8C2zadRsQTVHhECkMo6M65qNMFKjUynUa8oNVI2LcNts0xnUTEE1R4RArDxk/hwA5ofZPpJOIVTfpBifKwdLzpJCKeoMIjUhiWvgGlq0HDXqaTiFcEEtzra22YCfu3m04jEvW0D4/IWZx1L5X768PWOXDp/4Bf31JSiNrcBAv+DkvegCseN51GJKpphEfkfC19E/zx2ntHCl+Z6tC4L6yYAFlHTKcRiWoqPCLnoSTHIG0iJA2ExIqm44gXtbvdvZjoqommk4hENRUekfMwyD8Psg5psrIUnRpt3EuVLH5NS9RFzoMKj0gB+bAZ6Z8FF14M1VuZjiNeZVnQ9jbYvRG2fmU6jUjUUuERKaDLfMup4/sFOtzp/lISKSpJAyGxMiz6h+kkIlFLhUekgP4YmMUPTgVo3N90FPG6QLy7g/fmL2DXRtNpRKKSCo9IATSzttDWt55/BXtqKbqER8sb3dWAS14znUQkKqnwiBTATYFUDjrFeT/UzXQUiRWJFSH5akh7D47uNZ1GJOqo8IjkUzV2caVvCZNCl3KYEqbjSCzpcCdkH4Ulr5tOIhJ1VHhE8mlE4DMA/h3sYTiJxJxKjaHBlbD4NYpz3HQakaiiwiOSD6U5wrX+r0m12/ETFUzHkVjU6V44tper/d+YTiISVVR4RPLhOv9sSlnHeD3Yx3QUiVU120GNdtwUSCVA0HQakaihwiNyjopznD8GZvJVKIW1Tm3TcSSWdbqX6tZuevsWmU4iEjVUeETO0VD/11xgHeaVoPbdEcMu6sFGuxq3BmYAjuk0IlFBhUfkHMSTzc2BGSyyG7PcaWg6jsQ6n4/Xgn1p7NtON1+a6TQiUUE7pokAtUen5nl8kH8eVax9PJB9a5gSieRtut2B+5zJ3B34iDlZKYAubyKSF43wiJyFnxC3+aeTZtdlvt3UdBwRALIJ8GqwPxf7NtPFt9p0HJGIp8IjchZ9fAup5fuVV4ID0F/REkk+DHXlR6c89wamoLk8InlT4RHJg58Qdwc+4ju7Bl/YF5uOI3KKbAK8EhzAxb7NdPatMR1HJKKp8IjkYaB/PvV8P/NCcAiOvl0kAmmUR+Tc6Ce4SC7iCHKPfyqr7Tp8brcyHUfkjE7O5Wnp26RRHpE8aJWWSC6u9s+hhm8Xj2SNRHN3pCicbXXgufow1JXbA9O4NzCFeVnJhXKfIl6jER6RM0ggi7sCH7HUbsA3djPTcUTylEUcrwQH0NK3ie6+FabjiEQkFR6RMxjm/5Iq1j7+Grwaje5INPgg1JWtdhUeCkwCO2Q6jkjEUeER+Z0SHOe2wDTmh5JYZDcxHUfknAQJ8FzwGhr4foRVE03HEYk4msMj8ju3BGZQ0TrIzcGrz+njC2sehsj5mmW3Ic2uR8rXT0LTwRBX3HQkkYihER6R36jMXm72z2BGqB0rnYtMxxHJJ4unsv8AB3+EJa+bDiMSUVR4RH7jT4HJ+AnxTPAa01FECmSx0xguugLm/RWO7jUdRyRiqPCInNDI2s5V/m+YEOrBDqey6TgiBdd9LBw/6JYeEQFUeEROcPjvwLscpAR/Dw4wHUbk/FRpCi2GweJxsHuT6TQiEUGFRwS41LeSLv41/D04kIMkmo4jcv66j4W4EvDpn00nEYkIKjwi2ccZG3iLzfaFvBW6wnQakcKRWAm6PgSbZ8PGz0ynETFOhUfk2/+jlu9XxgRHkK2dGsRL2twMFRrAp6MhmGk6jYhRKjwS2/ZlwPy/MSPUjgV2U9NpRApXIB56PgV7t8KiV02nETFKf85KTMhtc8DX4/5KR5/DE9nDwpxIJEzqXwaN+sCcZyBpIJSrbTqRiBEa4ZGY1d23nCv8y/l7cCA7KW86jkjRufIZ8Plhxp/AcUynETFChUdiUimO8kTcv1hv12B8qJfpOCJFq0x16D4GtnwJ6VNMpxExQoVHYtJDgYlUYh8PZd+kicoSG1qPgmot3QnM2oFZYpAKj8ScttZ3DA98yT9DV7LKqW86jkh4+PzQ9//csvP5/5hOIxJ2KjwSUxLI4qm4N/jersTfgkNMxxEJryrJ0OleSHsHNswynUYkrFR4JKY8FJhEXd9O/hwcxTGKmY4jEn5dH4LKTWH63XBkj+k0ImGjwiMxo7NvNSMDn/KvYA/tuSOxK5AAA1+DY/sgVau2JHao8EhMKMdBno8bx0a7Gk8Hh5qOI2JWlaZwyZ9h3cdatSUxQ4VHvM9xeDJuPOU4xL3Zd5BJvOlEIuZ1uAeqt4EZ98HebabTiBQ5rceVqJHbbsknZTzd+8wHlv+bK/1LeTJ7KOuc2oUfTCQa+QMw+E14rTNMvhFGfuae7hLxKI3wiLf9lAazHmJuKJk3Q7kUIpFYVa4W9H8VfloJs8eaTiNSpFR4xLuO7YMProeSFbgn+w5s/XcXOV3jPtD2Vlj8D/huhuk0IkVGvwHEm2wbProNDv4EV01gH6VNJxKJXJf/BS5sAR/fBrs2mk4jUiRUeMSb5j4HG2dBj/+FGq1NpxGJbIEEuPpt9/WkoXBsv+lEIoVOhUe8J30qzHkSml0LbW42nUYkOpSt4Zaefd/DlD+CHTKdSKRQaZWWeMsPy91h+RrtoN9LYFmmE4mEXYFXNNZqD72egxn3wuwx7gipiEeo8Ih37N/hDscnVoZr39USW5GCaHUj/LoOFr4MZWpAu1tNJxIpFCo84gnlOAjvDILsY3D9NChZwXQkkejV82l3wv+noyGxEjQdZDqRyHnTHB6Jeokc5d/xz8L+7fCH96FSY9ORRKKbz+9uSlizHXx0C2ybZzqRyHlT4ZGolkAWr8f9jSQrA66aALU6mI4k4g1xxWHoRLigHkwcCjuWmE4kcl5UeCRqJZDFa3Ev0MG/jvuzb4WGPU1HEvGW4uXguo/c01pvD4IdS00nEikwFR6JSsXI5M245+niW81D2Tcxze5kOpKIN5WuCiNmuPPi3hnkroQUiUIqPBJ1SnCcf8c/SwffWh7IvoX3Q5eYjiTibaUvdEtPiQvgrf6a0yNRSau0JKpU4ABvxj9HUyuD+7LvYLqtOTsi+VWgfXrKVIcRM91RnncGuZOam/QvooQihU+FR8KqwBuiAXWtn/h33DNUsA5yS/Z9fGm3LOx4IpKXMtXgxlnw3jXwwQ3Q+6/Q+o+mU4mcE53SkuiQ8S1T4h+luJXJtVmPqOyImFLiAnevq4uugNQ/wcz/glC26VQiZ6XCI5HNcWDRP2BCX/Y6pRiU9RirnXqmU4nEtvgScO170O4OWPKae4rr6F7TqUTypMIjkSvrCEy9yd3ttUFPBmQ9zg6nsulUIgLgD0DPJ2HAP2D7Ini9m1ZwSUTTHB6JTD+luWVn9ya49BHodD+HVs3K81PONj9IRM7f6d9nZUixHuHl4EtUfuNy4i4fAx3uBp/+npbIov+RElnsEMz7K7zZHTIPuZuedXlQPzxFIliaU59emU8y224JX4x1T3Ht32E6lsgp9FtEIkZDazv8swd8+Rdo1AduWwD1tMeOSDQ4SCK3Z98DfV6AHYvh1Xaw9E2wbdPRRACd0pIIUIxM7g1MZZQ/lT07SvKX7NuZtqIjrFhoOpqI5IsFrUZCvUvhk3sg9X5YMxmufBaqNjMdTmKcRnjEGB82Q/zf8FXC/dwa+ITJoS50z3z+xGUiLNPxRKSgytWG6z6G/q/Arg3wWheYfhcc/tV0MolhGuERAxy6+VbxUGAijX07SLPrck/WnSx1GpkOJiKFxbKgxXBo1BvmPg+Lx0H6VHdCc9tboHhZ0wklxliO4+R1PM+DImeS22opC5srfMu5LTCdFN8WMuzKPBu8hpl2WzSiI+INue6WvnszzB4DG1IhoQy0uxXa3upuZChSeHL9ZaLCI4Xu94WnOMfp51/ITf5U6vt+4nu7Eq+F+vJhqCvZGmQU8ZS8Lg8DwM+r4JtnYf0MiE+ElGHQ5iaocFF4AorXqfBI+JwsPA2t7fzB/yUD/fMpbR1jnV2LfwT7MtNuSwi/4ZQiEm6nlKFf1sK3/wdrP4JQFtS9BFrdCA16QiDBXEiJdio8EiZ7t/HMX5+mj38hSb7vyXTimGm34d1gd5Y5DdGpK5HYdcbRn8O/wooJsPSfcOgn93RXUn9IvhpqddQeXJJfKjxSRGwbdq6CzV/A+pnw0woAVtr1+STUno9CHdlHacMhRSQS5Hm6KxSEbd/A6g/c011ZhyGxMlx0uTvqU7cbJJQKV1SJXio8UkhsG/ZscjcWy5gPm7+Eo7vdYxe2gCYD6JRalh+cimZzikjEOev8npOyjsKGmbA+1f0Zk3kA/PFQoy3U7gQ120P11u5FTEVOpcIjBeA4cOAH91z7L2tgx1L4YQkc2+ceL1HB3WCs/mXu60S35OiaViJyJudceH4rlO1enHTjp+4I0M50wAFfAKqmQLWWULUZvT44wCan+hkXQhTocSVa5Vp4tERGIPMw7P8e9m6DfRmwdyv8ug5+Wef+ZXVShYbuJR9qtnP/0ipf391rQ0SkqPjjoE5n9wXg2H7YsQS2L4DvF8LKd2DJEWYmQJbjZ7NTna1OVbY6VdhmV2WbUxWO7oXi5XL9eZXXH2kqS96hwmNYUX2j1Rs9nVIcpYx1hLIcpqJ1gErWfipZ+6jIAYYlJcChne4IzpHf7X6aUAYqNYbkIVA5CSo3dW8X01wcETGseFlocIX7Au4Fh/du5a6//Zsk3/c0tLaTZG2jp28JgcCJ63g9OwbiSkDpC0+8VHNfl6oKJcrT3reZ/U4i+5xE9lGKTOLNfX1SZLxfeBzHfXFv/Od953r7fD/XccAO5vrSzNpCgBB+bAJW6MTbIQLYsDYbso9B8BhkHz/x+sRL8Ph/3s4+CscPnPKypdjhXP9J9jilYF9NKFUZKjeBcnXggjrudvDl6uT5lxDolJWIhN/Zf+504BO7Q86tOILUsH6ljvUz4/uWh4M/wcEf3dcZ893XTgiAib/rN0edBPaRyBGnGLzxN4gv6e4ZFF/yNy+JEFfcXULvjwN/gjvPyB934n3x/3kJxIMvDnx+sPxg+U687fvN2+fw/pNna3J+Pv/utkbc82S28OxMh/FXkL8ycY4fGyWm57XdxIe5vD+uBASKnfhmK+ZO3CtWFi6o674uVoa/zfuFg5TgoFOCA5Rkl1OWX52y7KYMQQJk3J776JEKjYhEu2wCbHUuZKtzIbQ/w887OwRHdsPRPVzz4kzKWoe4wDpEOQ5T1jpMOQ5R0jpOg2KlIOuIe1os67D7dtYRyD4S/i8q335TgHIrSQW9XRAtb4CeTxX8889TnpOWLcv6FKgQvjhRrwKw23QIOY2el8ij5yQy6XmJPHpO8me34zg9z3TgbKu0JB8sy1rmOE4r0znkVHpeIo+ek8ik5yXy6DkpPNrCUkRERDxPhUdEREQ8T4WncL1uOoCckZ6XyKPnJDLpeYk8ek4KiebwiIiIiOdphEdEREQ8T4WniFiWdb9lWY5lWVrWb5hlWc9ZlrXesqzVlmV9ZFlWWdOZYpllWT0ty9pgWdZmy7JGm84T6yzLqmFZ1teWZa2zLGutZVn3mM4kLsuy/JZlrbQsa4bpLF6gwlMELMuqAVwBbDedRQCYDTR1HKcZsBH4s+E8McuyLD/wCnAl0AQYallWE7OpYl4QuN9xnCZAO+AOPScR4x7gO9MhvEKFp2i8APwX0bbts0c5jvO54zjBEzcXAdVN5olxbYDNjuNsdRwnC5gE9DecKaY5jvOz4zgrTrx9CPcXbDWzqcSyrOpAb+BN01m8QoWnkFmW1R/40XGcVaazyBmNBGaZDhHDqgE7fnP7B/TLNWJYllUbaAEsNptEgBdx/3C2TQfxCu9fPLQIWJb1BVDlDIceBv4b93SWhFFez4njONNOfMzDuMP374Yzm0g0sCwrEZgC3Os4zkHTeWKZZVl9gF8dx1luWVY303m8QoWnABzHuexM77csKxmoA6yy3AutVQdWWJbVxnGcnWGMGHNye05OsixrBNAH6O5oLwaTfgRq/OZ29RPvE4Msy4rDLTvvOo4z1XQeoSPQz7KsXkAxoLRlWe84jjPccK6opn14ipBlWRlAK8dxdOE3gyzL6gn8DejqOM4u03limWVZAdyJ491xi85S4A+O46w1GiyGWe5fZxOAvY7j3Gs6j5zqxAjPA47j9DGdJdppDo/EgpeBUsBsy7LSLMsaZzpQrDoxefxO4DPcybEfqOwY1xG4Drj0xPdH2omRBRFP0QiPiIiIeJ5GeERERMTzVHhERETE81R4RERExPNUeERERMTzVHhERETE81R4RERExPNUeERERMTzVHhERETE8/4fSQmjCgXQoRkAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<Figure size 720x432 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "def sample_and_display(init_state, stepsize, n_total, n_burnin, log_prob):\n", | |
| " chain, acceptance_rate = build_MH_chain(init_state, stepsize, n_total, log_prob)\n", | |
| " print(\"Acceptance rate: {:.3f}\".format(acceptance_rate))\n", | |
| " fig, ax = plt.subplots()\n", | |
| " plot_samples([state for state, in chain[n_burnin:]], log_prob, ax)\n", | |
| " despine(ax)\n", | |
| " ax.set_yticks(())\n", | |
| " plt.show()\n", | |
| " \n", | |
| "sample_and_display(np.array([2.0]), 30, 10000, 500, log_prob)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Not so cool, right?\n", | |
| "Now you could think the best thing to do is do choose a tiny step size.\n", | |
| "Turns out that this is not too smart, either, because then the Markov chain will explore the probability distribution only very slowly and thus again won't converge as rapidly as with a well-adjusted step size:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 14, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Acceptance rate: 0.985\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAFlCAYAAADvZjI4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5hU5eH28e+Z2V1EFhABKyiIglQxIGKQotgiSDTozwKiIjGWqDGIwWBQDCq2aIwaYkleo6AIdhEVC5ZgVEAQsKKuCmooAtLZmTnvHyMbCIKU3T2zZ76f65prdvZMuWcHlpvnPOc8QRiGSJIkxVki6gCSJEkVzcIjSZJiz8IjSZJiz8IjSZJiz8IjSZJiz8IjSZJir+BHtnvMuiRJqiqCTW1whEeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSpAq2ZMkS7rzzzkp7vUaNGrFw4UIAfvrTn272vtdee+1mtx977LEsWbKEkpISWrVqtVU5Jk2axOTJk8tujxw5kn/+859b9RzlJQjDzZ5M2TMtS5K0nUpKSujZsyezZs3aaFsqlaKg4McWPtg6jRo1YsqUKdSrV+9H71tcXMzy5cs3+n4YhoRhSCKRHRvZ3HvYlKuuuori4mIuvfTSLQ+/fTzTsiRJURk8eDCffPIJbdu2ZdCgQUyaNInOnTvTq1cvWrRosdHoyU033cRVV10FwCeffMIxxxxDu3bt6Ny5Mx988MFGz79o0SKOOuooWrZsyYABA1h/MKO4uBiAr7/+mi5dutC2bVtatWrFa6+9xuDBg1m1ahVt27alT58+lJSU0KxZM/r160erVq348ssvNxgtSqVS9OnTh+bNm3PiiSeycuVKYMMRpSlTptCtWzdKSkoYOXIkt9xyC23btuW1117jqquu4qabbgJg+vTpdOzYkTZt2nDCCSewePFiALp168bvfvc7OnToQNOmTXnttdfK5TMo30opSVKumzAYvplZvs+5W2v42YhNbh4xYgSzZs1i+vTpQHZXz7Rp05g1axaNGzempKRkk48955xzGDlyJPvttx9vvvkm559/Pi+99NIG9xk2bBiHHnooQ4cOZfz48dx7770bPc/o0aM5+uijGTJkCOl0mpUrV9K5c2duv/32slwlJSV8/PHH3HfffXTs2HGj5/jwww+599576dSpE/379+fOO+/c5OhNo0aNOPfcczcY4XnxxRfLtvfr14+//OUvdO3alaFDhzJs2DBuvfVWIFus3nrrLZ555hmGDRvGCy+8sMmfz5ay8EiSFIEOHTrQuHHjzd5n+fLlTJ48mZNOOqnse2vWrNnofq+++iqPPvooAD169KBOnTob3eeggw6if//+lJaWcvzxx9O2bdsffM299977B8sOQMOGDenUqRMAffv25bbbbtum3VVLly5lyZIldO3aFYAzzjhjg/f4i1/8AoB27dpttgxuDQuPJCm/bGYkpjLVqFGj7OuCggIymUzZ7dWrVwOQyWTYaaedykZgtkeXLl149dVXGT9+PGeeeSa//e1v6dev32Zz/a8gCH7w9vr512XfHtWqVQMgmUySSqW2+/nAOTySFCuNBo8vuyh31KxZk2XLlm1y+6677sr8+fNZtGgRa9as4emnnwagVq1aNG7cmLFjxwLZicQzZszY6PFdunRh9OjRAEyYMKFsPsz6Pv/8c3bddVd++ctfMmDAAKZNmwZAYWEhpaWlW/Q+vvjiC9544w0gu4vs0EMPBbK7r6ZOnQrAI4888qPvu3bt2tSpU6dsfs79999fNtpTUSw8kiRVsLp169KpUydatWrFoEGDNtpeWFjI0KFD6dChA0ceeST7779/2bZRo0Zx7733csABB9CyZUueeOKJjR5/5ZVX8uqrr9KyZUseffRR9tprr43uM2nSJA444AAOPPBAxowZw8UXXwxk5wi1adOGPn36/Oj7aNasGXfccQfNmzdn8eLFnHfeeWWvf/HFF9O+fXuSyWTZ/Y877jgee+yxsknL67vvvvsYNGgQbdq0Yfr06QwdOvRHX397eFi6JMXI+iM7JSN6RJhEioSHpUuSpPxl4ZEkSbFn4ZEkSbFn4ZEkSbFn4ZEkSbFn4ZEkSbFn4ZEkqRIEQUDfvn3LbqdSKerXr0/Pnj3LvjdhwgTat29PixYtOPDAAxk4cCCQXXU8CALmzJlTdt9bb72VIAiYMmUKkF2G4le/+hVNmjShXbt2dOvWjTfffLOS3t0PO/PMMxk3blykGdZxaQlJUt4p7zNRb8k5j2rUqMGsWbNYtWoV1atXZ+LEiey5555l22fNmsWvf/1rxo8fz/777086neauu+4q2966dWseeughrrjiCgDGjh1Ly5Yty7YPGDCAxo0b8/HHH5NIJPjss8947733yvFdVm2O8EiSVEmOPfZYxo/Plq0HH3yQU089tWzbDTfcwJAhQ8rOspxMJsvOZAxw/PHHl51l+ZNPPqF27drUq1ev7Pabb77J8OHDSSSy/7Q3btyYHj02LGLpdJozzzyTVq1a0bp1a2655RYA7r77bg466CAOOOAAevfuzcqVK4HsCM15551Hx44d2WeffZg0aRL9+/enefPmnHnmmWXPW1xczCWXXELLli3p3r07CxYs2Oi9T506la5du9KuXTuOPvpovv76awBuu+02WrRoQZs2bTjllFO2/Yf7Iyw8kiRVklNOOYWHHnqI1atX8+6773LwwQeXbZs1axbt2rXb5GNr1apFw4YNmTVrFg899BAnn3xy2bbZs2fTtm3bDZZ1+CHTp09n3rx5zJo1i5kzZ3LWWWcB2dXJ3377bWbMmEHz5s259957yx6zePFi3njjDW655RZ69erFJZdcwuzZs5k5c2bZoqYrVqygffv2zJ49m65duzJs2LANXre0tJQLL7yQcePGMXXqVPr378+QIUMAGDFiBO+88w7vvvsuI0eO3MKf5Naz8EiSVEnatGlDSUkJDz74IMcee+xWP35dYXr88cc54YQTtvrx++yzD59++ikXXnghzz77LLVq1QKyZatz5860bt2aUaNGMXv27LLHHHfccQRBQOvWrdl1111p3bo1iUSCli1bUlJSAkAikSgrYH379uX111/f4HU//PBDZs2axZFHHknbtm0ZPnw4c+fOLfuZ9OnThwceeICCgoqbaWPhkSSpEvXq1YtLL710g91ZAC1btixbcXxTevbsyf33389ee+1VVlbWPXbGjBmk0+nNPr5OnTrMmDGDbt26MXLkSAYMGABkd13dfvvtzJw5kyuvvJLVq1eXPaZatWpAttSs+3rd7VQq9YOvEwQbLmkVhiEtW7Zk+vTpTJ8+nZkzZ/L8888DMH78eC644AKmTZvGQQcdtMnn3F4WHkmSKlH//v258sorad269QbfHzRoENdeey0fffQRAJlMZqNdPDvuuCPXX3992e6gdZo0aUL79u258sorWbcoeElJSdl8oXUWLlxIJpOhd+/eDB8+nGnTpgGwbNkydt99d0pLSxk1atRWv6dMJlN2NNbo0aM59NBDN9jerFkzFixYwBtvvAFkd3HNnj2bTCbDl19+yWGHHcb111/P0qVLWb58+Va//pbwKC1JkipRgwYNuOiiizb6fps2bbj11ls59dRTWblyJUEQbHDI+jqbmth7zz33MHDgQPbdd1+qV69OvXr1uPHGGze4z7x58zjrrLPIZDIAXHfddQD88Y9/5OCDD6Z+/focfPDBLFu2bKveU40aNXjrrbcYPnw4u+yyC2PGjNlge1FREePGjeOiiy5i6dKlpFIpfvOb39C0aVP69u3L0qVLCcOQiy66iJ122mmrXntLBeua4CZsdqMkKbesf7j1lhwqLZWH4uLiChuZ2UrBpja4S0uSJMWehUeSJG2XHBnd2SwLjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJij0LjyRJir2CqANIkjbWaPD4DW6XjOgRURIpHhzhkSRJsWfhkSRJsWfhkSRJsWfhkSRJseekZUmKSHlMTP7f55D0wxzhkSRJsecIjyRVMY7qSFvPER5JkhR7Fh5JkhR7Fh5JkhR7zuGRpJhaf66PS1Mo3znCI0mSYs/CI0mSYs9dWpKUIzZ3uLmHokvbxxEeSZIUexYeSZIUexYeSZIUexYeSZIUexYeSZIUexYeSZIUexYeSZIUexYeSZIUexYeSZIUexYeSZIUey4tISkvuZL4lvHnpLhwhEeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMVeQdQBJEm5o9Hg8Vu8rWREj4qOI5UbR3gkSVLsWXgkSVLsuUtLkvKAu6OU7xzhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsedRWpLynkcwSfHnCI8kSYo9C48kSYo9d2lJyhubWydqSx7jri6p6nKER5IkxZ6FR5IkxZ6FR5IkxZ5zeCQpDzk3SfnGER5JkhR7Fh5JkhR7Fh5JkhR7zuGRVDnCEFYshNVLIZOCohpQczdIFkadTFIesPBIqhiZDHzxBsyZCJ+9Bgs+gLXL/+dOAezcGBoeDI27wv49YIdakcSVFG8WHknla9USePtumHY/LPkcEgWwZ3s4sC/UaQw71oVEAtYsg+++hv/Mgo8nwowHIVkNmh8HnS6C3Q+I+p1IihELj6TyUboKJv8FJt8Oa5ZmR2wOvwKaHQvVijf/2DCEuW/DzLEw/UGYNQ72OxqOGg71m1ZOfkmxZuGRtP0+fgGeGQiLS2D/ntD1sq0boQkCaNghezlsCEy5F16/Ff56CBxyAXT7PRTuUGHxJcWfR2lJ2napNfDMZTCqNySL4Iyn4JRR27c7qvpO0HkgXDgNDjgV/vVnuKsbfP1uucWWlH8sPJK2zbefwj1HwFt/g4PPg3Nfh8Zdyu/5i+vDz2+HPo/Aqm+zrzVjTPk9v6S84i4tSVvvy7fgwVMgk4ZTHoT9j62419rvCDhvMjx8Bjx2DsyfDd2vhESywl5y/WUXJMWDIzySts57T8B9x0G1WvDLlyq27KxTox70exzan53dxfXI2ZBaW/GvKyk2HOGRtOVmPASPnQsN2sOpD2WLSGVJFkLPP2XP2/P8Fdmjwk66z8nMkraIIzyStsy6stO4M/R7snLLzvp+eiH0vAU+eg5G/1+2+EjSj7DwSPpxM8f9t+ycOgaKdow2T/v+cPxf4bNXYexZkC6NNo+knGfhkbR5n7wEj/0K9u6UG2VnnbanQo+b4KMJ8MQF2aUsJGkTnMMjadO+fhfG9IN6zeDU0blTdtY5aACsWgwvDYca9eHoa6JOVCV5VJrygYVH0g9bOg9GnQQ71Ia+47LXuajzpbB8PrxxO9RrCu3OiDqRpBxk4ZG0sdLVMKYvrF0BAyZCrT2iTrRpQQBHXweL5sD438LO+2TnGknSepzDI2lDYZgtDl9Ng1/8DXZpHnWiH5csgBP/kS07D58Oiz+POpGkHGPhkbSht+6G6aOg6+9g/x5Rp9ly1XfKnhsok4axZ2bX+VKlaTR4fNlFykXu0pL0X1++Bc9dDk1/Bl0HR51m69VtAj+/IzvK8/wVNHq1W9SJJOUIR3gkZa1aAuPOhlp7ZndlJaror4cWvaDjBfDWXfRMvBF1Gkk5oor+RpNUrsIQnroYln0FJ/49d4/I2lJHDoOGBzOi8G72Dr6JOo2kHGDhkQTT/gnvPQ6HDcmuk1XVJQvhxL+TJsEthXeSJB11IkkRs/BI+W7hxzDhd7BPN+j0m6jTlJ/aDfhDaX9+kpjD+cknok4jKWJOWpbyWSadXSOrcAc4oQrP29mEJzM/5fD0NC4ueJRXMgfwbtgk6kixsqkjsv73+yUjqtDRfoqteP12k7R1Jv8F5k2BY2+CmrtFnaZCDC09i/9Qh1sK72QHPFRdylcWHilfzf8AXr4Gmh8HrXpHnabCfEcNLi09lyaJr7mkYFzUcSRFxF1aUj5Kp+Dx86CoGHr8Kbs8Q4y9kWnJqFR3BiSfYXy6Y6S7tjwxnxQNR3ikfDT5z9mlI3rcDMW7RJ2mUoxIncp86nB94V0Ukoo6jqRKZuGR8s2iT2DS9dC8F7T6RdRpKs0yduSK0rNonviS85JPRh1HUiVzl5aUT8IQnr4ECqrBsTdGnaZCbG6X0YuZdjyR/im/LniMCZkOfBw2qMRkkqLkCI+UT2aOhc9ege5DY3tU1o8ZVtqPZezIiMK7CchEHUdSJbHwSPli5bfw3O9hz3bQvn/UaSLzLbW4trQP7RIfc2Ly1ajjSKokFh4pX7xwVbb09LwVEsmo00TqkUxn3so04/KC0ezEsqjjSKoEFh4pH3zxb5h2H3Q8D3ZvE3WaHBDwh9KzqMVKLisYE3UYSZXAwiPFXToFT/8WajWAbpdHnSZnfBjuxT/Sx3BK8mXaBnOijiOpgll4pLib8neYPxuOuQ6qFUedJqfcmurNfHZieOHfSTiBuVI0Gjy+7CJVJg9Ll+JsxSJ4eTg07ppdQkIbWEF1/lh6OncU3Ubf5ET+mT56s/d3UUyp6nKER4qzl4fDmuXws+tjv3zEthqfOZjX0q0YWDCWOnwXdRxJFcTCI8XV1+/ClH9Ah1/CLs2jTpPDAq5O9aMGq7mk4JGow0iqIBYeKY7CECb8DnbcGboNjjpNzvs4bMCodHf6JF+gafBl1HEkVQALjxRHsx6BLyZnz6hcvU7UaaqEW1Inspzq/KHgfiCMOo6kcmbhkeJm7QqYOBR2awMHnh51mipjCTW5NdWbzslZdE9M26LHeMSRVHVYeKS4+ddt8N287OKgeX5G5a11f/pI5mT2YEjBKApJRR1HUjmy8EhxsuwbmHwbtDge9uoYdZoqJ0UBw1N92SfxDWckn4s6jqRyZOGR4mTSdZAuhSOujDpJlTUp05aX0wdwUcGj1GVp1HEklRNPPCjFxfwPYNo/ocM5sPM+Uaep0oan+vJs0WAGFozl96kBUcep8pzjpFzgCI8UFy9cBUXF0OWyqJNUeZ+Ee/JA+ghOTr7MfsHcqONIKgcWHikOSl6HjybAoZdAjbpRp4mF21InsIId+F3Bg1FHkVQOLDxSVZfJwPNXZFdD73he1GliYzG1uDP1c45IvkPHxHtRx5G0nSw8UlU3+1H46h04/AoorB51mlj5R/oY5oV1+X3BKAJXU5eqNAuPVJWl1sCLV8OuraHN/0WdJnbWUMTNpSfRJvEZxyX+HXUcSdvBwiNVZW/fA0s+h6Ou9iSDFeSxzKHMzuzNoIIxFFEadRxJ28jCI1VVqxbDKzdAk8OzF1WIkATXpk6jYWIB/ZLPRx1H0jay8EhV1Wt/gtVL4ciro04Se//KtOaVdBsuLHiM2iyPOo6kbWDhkaqiJV/Am3+DA06F3VpHnSYvXJc6jZqs4oKCJ6KOImkbWHikquil4RAEcPiQqJPkjQ/CvRiX7sIZyedoEMyPOo6krWThkaqar6bDu2Oy59yp3SDqNHnlT6kTyZDg0oKHo44iaStZeKSqJAxh4h+g+s7ZsyqrUn1DXe5N/4zjk5NpHXwadRxJW8HCI1Ulc16Az16Frr+DHWpHnSYv/S11HIvCmvy+YDQQRh1H0hay8EhVRSYNE4dmV0Jv3z/qNHlrGTvy59QvOCT5HoclpkcdR9IWsvBIVcX00TD/Peh+JRQURZ0mr41Od+fTzG5cXjCaJOmo40jaAhYeqSpYuwJevgYaHAQtfh51mryXooAbUqfQNDGP3slXo44jaQtYeKSq4N93wrKv4cg/Zg9HV+SezRzEtMy+DCwYS3VWRx1H0o+w8Ei5bvkCeP3PsH9P2PuQqNOoTMA1pX3YNVjC2ckJUYeR9CMsPFKue+V6KF0JR1wVdRL9j6lhM55NH8S5BU9li6mknGXhkXLZwjkw9R/Q7kyot1/UafQDbkidzA6szRZTSTnLwiPlshevgoIdoNvgqJNoEz4N9+DB9OHZYrpwTtRxJG2ChUfKVV/8G95/CjpdDMW7RJ1Gm/HnVO9sMX3xqqijSNoEC4+Ui8IQnv8DFO8Gh1wQdRr9iIXUzhbT95+CL96MOo6kH2DhkXLRe0/A3LfgsN9DUY2o02hLHHIBFO8Kz1+RLayScoqFR8o1qbXwwlWwSws4sG/UabSlimpkC+rct7IjPZJyioVHyjVT7oXFn2VPMphIRp1GW6NtX6i/f7awpkujTiNpPRYeKZesWpI9vHmfbrBv96jTaGslC+CIYfDtJzD1/0WdRtJ6LDxSLnnt5mzpcQmJqqvp0bD3oTBpBKz+Luo0kr5n4ZFyxeLP4c2R0PY02L1N1Gm0rYIAjroaVi6EybdFnUbS9wqiDiDpey9eDUESDhsSdRJtrz3bQaveMPl2aH82ja6dFnUiKe85wiPlgnlTYda47KHNtfeMOo3Kw+F/gEwKJl0bdRJJOMIjRW/dSQZr1IdDfxN1miqh0eDxZV+XjOgRYZLN2LkxdPglvDmS/YJWfBw2iDqRlNeCcPMnyPLsWVJF+2A8PHQa9PgTHHR21GmqhPULTy4qK2Erv4U/t+XFlftwdumgaEPloJwtq6rKNnm0h7u0pCilS2HilVCvKfzkjKjTqLztuDN0/i3dk+9wSGJ21GmkvGbhkaI09f/Boo/hyKuz53BR/Bz8K+aG9RhSMIoEmajTSHnLwiNFZdViePlaaNQZmh4TdRpVlMLqXF96Cq0SJZyYfCXqNFLesvBIUXnlhmzpOfpaTzIYc09lDuHtTFMGFYyhmJVRx5HykoVHisKCj+Ctu6DdGZ5kMC8EXF3aj/rBd/y64PGow0h5ycIjReG530PhjnDYFVEnUSWZGe7D2FQX+icnsHfwTdRxpLzjLEmpsn30PMyZCEddA8X1o05TZeT6oehb4obUyfws+RZDCkZxTunAqONIecURHqkypdZmR3fq7gsdzok6jSrZAupwR+p4jkpOpVNiZtRxpLxi4ZEq09t3Zw9DP/paKCiKOo0i8Pf0MXyRqc/QgvtJko46jpQ3LDxSZVmxECZdD026w35HRZ1GEVlDEdek+tIsMZdTky9FHUfKGxYeqbK8NBzWLodjrvMw9Dz3XKY9k9MtGFgwltosjzqOlBcsPFJlmDc1e1blDudA/WZRp1HkAq5O9aMmKxlUMCbqMFJe8CgtqaJl0jB+IBTvAoddHnWaKiMOR2VtzgfhXtyXPpqzks/ycLob74ZNoo4kxZojPFJFm/oP+Oqd7GHoO9SOOo1yyC2p3iygNsML/+46W1IFs/BIFWn5Anjx6ux6Wa1PjDqNcsxyduSa0r60SXzGackXo44jxZqFR6pIE4fC2pXQ42YnKusHPZk5hH+lWzKoYAx1WRp1HCm2nMMjVZTP34AZo+HQS5yorM0IGJo6kwlFg7m88EEuLT036kCV5n/naZWM6BFREuUDR3ikipAuhfG/hdoNocugqNMox30S7sk96R6cmHyVg4IPoo4jxZIjPFJF+PdfYf57cPIDUFQj6jSqZNtyhNlfUsfTKzmZawrvpcfa6yj117NUrhzhkcrbt5/Cy9dCs2Nh/55Rp1EVsYodGFp6Jk0T8zg/+UTUcaTYsfBI5SkM4amLIVnoRGVttZcyP+HJ9CFcUPA4+wVzo44jxYqFRypP7zwAn70KR1wFtfaIOo2qoGGl/VhOda4vvMtz80jlyMIjlZdl/4Hnh8BeP4V2Z0WdRlXUImpzdWk/fpKYQ7/k81HHkWLDwiOVlwmDoHQ19LoNEv7V0rZ7PNOJl9MHMKhgDA2CBVHHkWLB38pSeXj/KXjvCeh6GdTbL+o0qvICrijtD8C1BfcAYbRxpBiw8Ejba8VCePoS2K01dLo46jSKiXnU54bUKXRJzuSk5CtRx5GqPAuPtD3WHZW1eimccFf26CypnNyfPpI3M/sztOB+d21J28nCI22Pd8fAB0/DYUNg1xZRp1HMZEgwsPRcAkJuKhxJ4FFb0jaz8EjbaulceOYyaNgRfnph1GkUU3PDXRiW6kfHxPucnZwQdRypyrLwSNsiDOGJX0MmBSf8FRLJqBMpxsamuzIx3Y5BBWNoGnwZdRypSrLwSNvirbvh05fhqD/CzvtEnUaxFzC4dADL2JFbCu+kkFTUgaQqx8Ijba1vZsLzV8B+R0H7/lGnUZ5YRG0uLx1Ay8TnDCx4OOo4UpVj4ZG2xtoVMPYsqF4Hjv+ra2WpUk3MtOeBVHfOLXiabol3oo4jVSkWHmlrPHMZLJoDve+GGvWiTqM89MfU6byf2Ys/Ff6V3VgUdRypyrDwSFvq3Ydh+gPQZRA07hJ1GuWpNRRxQelFVKOUPxfdQZJ01JGkKsHCI22JhR9nz6a81yHQ9XdRp1Ge+zTcg9+Xns3BiQ/4TcEjUceRqgQLj/Rj1iyDh/pAwQ7Q+x5IFkSdSOKJzKGMSXXjguQTdEnMiDqOlPMsPNLmhCE8fj4s+hhO+gfUbhB1IqnMlakz+ChswG2Ft9Mw+E/UcaScZuGRNmfybQpt/KUAAAyfSURBVPD+k3Dk1c7bUc5ZTTXOKf0tAHcV/okdWR1xoqqj0eDxZRflBwuPtCmfToIXroKWJ8Ahv446jfSDvgh35cLSC2kazOXGwpFAGHWkbWYJUUVyMoL0QxZ+DA+fAfWaQa/bPd9OJfEfum3zWqYNI1KnMqRwNLMzT3Bn+vioI0k5x8Ij/a8Vi2DUSZAogNMegmrFUSeSftTd6R60SpRwacFYPg334NlMh6gjSTnFXVrS+lJrYExf+O4rOPVBqNMo6kTSFgq4rPQc3gn35dbCO/hJ8FHUgaScYuGR1glDePIi+GJydgX0hv4PWVXLGor45dqBfBPuzN1FN7N38E3UkaScYeGR1nlxGLz7EBw2BFr1jjqNtE2+pRZnlV5GgpB/FN5AHb6LOpKUEyw8EsC/boPXb8muft5lUNRppO3yWbg7A9YOZM9gEf8oupEarIo6khQ5C4807X6Y+Ado+Qs49iaPyFIsTA2b8evSC2kVfMY9hTdTjbVRR5IiZeFRfnvvSXjqImjSHU74GySSUSeSys3ETHt+W3oeByfeZ2ThLRSSijqSFBkLj/LX+0/BuLNgz/Zw8v1QUBR1IqncPZnpxJBUfw5LzuDWwttdXV15y8Kj/DT78eyJBff4CfQdB0U1ok4kVZgH090ZXtqHHsm3+EvhXxzpUV7yxIPKP7MegUd+mT3svM9YqFYz6kRShbsn3YMMCYYW3k81buH80otZg6OasOEZvktG9IgwiSqSIzzKL9Puh0cGwF4doc84y47yyt/TP+P3pWfTPfkO9xbeSHUXG1UecYRH+SEM4dWb4OXh0ORwOPkBd2MpL41Od2d1WMiNhX9jVNG1DFh7Kd9SK+pYFc512uQIj+Ivk4bxA7Nl54BT4bSHLTvKa49munB+6cW0CD7n0aIraRR8HXUkqcJZeBRva5bBmNNhyr1w6CVw/F8hWRh1Kilyz2U6cNraIdQMVvJo0ZW0Cz6MOpJUoSw8iq9Fn8A9R8JHz8LPboQjrvKkgtJ6poVN+cXaYSwJixlddC29E69GHUmqMBYexdOcF+Huw2D5N3D6o3DwOVEnknLS5+Fu/GLtMKZm9uPmopH8seDvFFEadSyp3DlpWfGSScNrN8Ok66B+czhlFOzcOOpUUk5bQk1OL72cQeHDnFvwFK0SJZy39mK+oW5kmTxUXOXNER7Fx9K5cF8vePma7LpYZz9v2ZG2UJokI1Kncu7a37BfMJdnql3O0Ym3oo4llRsLj+LhvSfhr53gq3eyE5N73wPViqNOJVU5z2Y60GvtcOaG9flb0a3cUPA3V1tXLLhLS1Xb8gUw4TKY/Sjs3hZO/DvUbRJ1Km0Fz4+Sez4N96D32mFcVPAo5yef4ODE+wxO/ZI3Mi2jjraR//3zs/7uL/9saX0WHlVNYQjvjoFnB8PaFXDYFdDpYhcAlcpJKQXcnPo/Xkm34ebCkTxYdA3j0l24pvQ0Fkd8okKLjLaFu7RU9XwzE+47Dh77FdRrCue+Dl0HWXakCjAl3J+j1t7AHale/DzxL16sdiknJSeRIBN1NGmrWHhUdSxfAE9eBCM7w39mw7E3wVkToH6zqJNJsbaGIm5MnUKPtdfySbgHNxbexdNFQzg0MTPqaNIWc5eWct/Kb+GNO+DNv0FqFXQ8D7peBtXrRJ1MyisfhQ35v7VD6Zn4N5cVjOGBout4Jd2GG1InMzv0iEjltiAMw81t3+xGqUKtWgxv3AlvjoQ130GLn2fn6tRvGnUylSPnY1RNRZTSL/k8FxY8Ru1gJS+l23J76nimhVX776fn/KnyNnk6fQuPcs/COfDmX2H6aChdCc2Pg66DYbdWUSdTBbDwVG01WcnpyecZUPAMOwfLmZxuwf9LH82LmZ+QJhl1vK1m4anyLDzKcelS+HgiTLsPPnouu8Bn65Og4/kWnZiz8MRDdVbTJ/kiZxU8y57BIuaG9Rid6s5D6cP4NuKjuraGhafKs/AoR83/AKY/ADPGwIr5UGMXaH8WtD8bau4adTpVAgtPvCRJc0RiGqcnn+fQ5GzWhkkmZdryWPpQXsocyBpy+2hKC0+Vt8nC46RlVa4whP/Mgvefyp4decH7kCiApsfAgX1h3yOyozuSqqQ0SZ7LHMRzmYNokprHKcmX6ZWczFHJqXwXVufZdAeez7Tn9UwrVlMt6rjKI47wqOKtXgqfvQafToI5L8DizyBIwF4/hRa9suteFdePOqUi4ghP/CXIcEhiNick/8XRibepGaxidVjI65lWvJT5Cf/KtOTzcFc285/zSuMIT5XnLi1VouXzYe4UmDcFPnsV5k2FMAOFNaBRJ9i/BzTrYckRYOHJN4Wk6JB4nyMS0zgiMY2GiQUAfBXuzL8zLfh3pjlTMs34LNyNMAdOFWcBqnIsPKoAYQjffZXdLTX/fZg3LVtylnyR3Z4oyK5v1eQw2OcwaHCQZ0PWRiw8+SykSfAVhyTeo2PifQ5OvEf94DsAloXVmR024t3MPszKNGZ2uDefh7uSquSZGBaeKsfCo20UhrBiISz5HBaXZC9LPocFH2YnHK9Z+t/71moADdpnL3u2h90PgKIdo0quKsLCo/8K2TeYx4GJObQOPqN14jNaBJ9TLSgFoDRM8kW4C5+Ee/BJuAefhrvzRWZX5lGX/4Q7U1oBZcjCU+U4aVnryWSyJ/JbvTR7WTE/u2zD8v9kd0etmP/fr5d8CaUrNnz8jvWyyzm0PhF2aZ691G8ONepG834kxUTAnLABc9INGEs3AApI0TSYS7PgS5okvqJJkL10S0ynKEiXPTITBsxnJ74K6/JVWJcF4U58G9bkW2qxKMxevqUmi8JafEcNMjmwu0yVK78Kz7rRrLJRrU3c3pL7bO9zZNLfX1LrXbbidmp19lK66vvr1dllFza6XgWrvy8360rOmmVscvCucEco3iV7eHjdfaHJ4bDT3lBn7+z1TntBteLN/pglqbykKOC9sBHvhY1Yf73SJGkaBvPZM1jIHsEi9mBR9jpYSPPgC7ok3qVWsGqTz7sirMZyqrM8rL7hNdVZGVZjDYWsoRAmzYKCalCwAxTukL1edztZBIlkdvd92SUJQXLD2xtcF2QP2ggCINiG6+8fG0Q/wbuqiXaX1jez4N6jqPCikS+CJBRW//4v5nrXO9TOXqrV+v7rWhveLt4FatSH4l0tM6p07tJSRSmilDoso27wHTsHy9iZ76gbfEctVlIjWE0xK6kZrKKYVRR/f10zWEV11lCNUqpRSuF6o0i5aTPlaLMP29z27XnsZh7f7gw45rofeex227Y5PEEQPAvUq4hEMVUPWBh1CG3EzyX3+JnkJj+X3ONnsnUWhmF4zA9t+LERHm2FIAimhGHYPuoc2pCfS+7xM8lNfi65x8+k/DhrS5IkxZ6FR5IkxZ6Fp3zdFXUA/SA/l9zjZ5Kb/Fxyj59JOXEOjyRJij1HeCRJUuxZeCpIEAQDgyAIgyDwsP6IBUFwYxAEHwRB8G4QBI8FQbBT1JnyWRAExwRB8GEQBHOCIBgcdZ58FwRBwyAIXg6C4L0gCGYHQXBx1JmUFQRBMgiCd4IgeDrqLHFg4akAQRA0BI4Cvog6iwCYCLQKw7AN8BFwecR58lYQBEngDuBnQAvg1CAIWkSbKu+lgIFhGLYAOgIX+JnkjIuB96MOERcWnopxC3AZeXeq59wUhuHzYRimvr/5b6BBlHnyXAdgThiGn4ZhuBZ4CPh5xJnyWhiGX4dhOO37r5eR/Qd2z2hTKQiCBkAP4J6os8SFhaecBUHwc2BeGIYzos6iH9QfmBB1iDy2J/Dlerfn4j+uOSMIgkbAgcCb0SYRcCvZ/zhnfuyO2jL5tXhoOQmC4AVgtx/YNAT4PdndWapEm/tMwjB84vv7DCE7fD+qMrNJVUEQBMXAI8BvwjD8Luo8+SwIgp7A/DAMpwZB0C3qPHFh4dkGYRge8UPfD4KgNdAYmBFkF1drAEwLgqBDGIbfVGLEvLOpz2SdIAjOBHoC3UPPxRCleUDD9W43+P57ilAQBIVky86oMAwfjTqP6AT0CoLgWGAHoFYQBA+EYdg34lxVmufhqUBBEJQA7cMwdOG3CAVBcAzwJ6BrGIYLos6Tz4IgKCA7cbw72aLzNnBaGIazIw2Wx4Ls/87uA74Nw/A3UefRhr4f4bk0DMOeUWep6pzDo3xwO1ATmBgEwfQgCEZGHShffT95/NfAc2Qnxz5s2YlcJ+B04PDv/35M/35kQYoVR3gkSVLsOcIjSZJiz8IjSZJiz8IjSZJiz8IjSZJiz8IjSZJiz8IjSZJiz8IjSZJiz8IjSZJi7/8D7JaqJzlQTikAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<Figure size 720x432 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "sample_and_display(np.array([2.0]), 0.1, 10000, 500, log_prob)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "No matter how you choose the step size parameter, the Markov chain will eventually converge to its stationary distribution.\n", | |
| "But it may take a looooong time.\n", | |
| "The time we simulate the Markov chain for is set by the `n_total` parameter - it simply determines how many states we will end up with.\n", | |
| "If the chain converges slowly, we need to increase `n_total` in order to give the Markov chain enough time to forget it's initial state.\n", | |
| "So let's keep the tiny step size and increase the number of samples by increasing `n_total`:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 15, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Acceptance rate: 0.990\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAFlCAYAAADvZjI4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3yV5cH/8c91zkkChC0KKELYI0DCEBCEIFuWuLXyWOTRttq6anm0PyxqtYrjqaPWB1eptSpuURE1jAACIkOQMBMwDBmyEkYgyTn3/fvjQApCQgInuc74vl8vXiG5T3K+IZB8ucZ9Gdd1EREREYlmHtsBRERERCqaCo+IiIhEPRUeERERiXoqPCIiIhL1VHhEREQk6qnwiIiISNTznea69qyLiIhIpDAlXdAIj4iIiEQ9FR4RERGJeio8IiIiEvVUeERERCTqqfCIiIhI1FPhERERkainwiMiIiJRT4VHREREop4Kj4iISAXLzc3lxRdfrLTnS0pKYvfu3QD07Nmz1Mc+9thjpV4fOnQoubm55OTk0L59+3LlyMjIYMGCBcWvT5o0iX/961/l+hihYly31Jsp607LIiIiZyknJ4fhw4eTmZl50jW/34/Pd7qDD8onKSmJJUuWUK9evdM+tnr16hw8ePCkt7uui+u6eDzBsZHSPoeSPPTQQ1SvXp0//OEPZQ9/dnSnZREREVvuv/9+NmzYQGpqKuPGjSMjI4PevXszcuRI2rVrd9LoydNPP81DDz0EwIYNGxgyZAhdunShd+/erF279qSPv2fPHgYNGkRycjK33HILxw9mVK9eHYDt27fTp08fUlNTad++PfPmzeP+++/n8OHDpKamcuONN5KTk0Pr1q256aabaN++PVu2bDlhtMjv93PjjTfStm1brr76avLz84ETR5SWLFlC3759ycnJYdKkSTzzzDOkpqYyb948HnroIZ5++mkAli9fTo8ePejYsSNXXHEF+/btA6Bv377cd999dOvWjVatWjFv3ryQfA1CWylFRETC3fT7YcfK0H7MBh3gsoklXp44cSKZmZksX74cCE71LFu2jMzMTJo2bUpOTk6J7/urX/2KSZMm0bJlSxYtWsTtt9/OrFmzTnjMww8/zCWXXMKECROYNm0ar7322kkf56233mLw4MGMHz+eQCBAfn4+vXv35oUXXijOlZOTQ1ZWFq+//jo9evQ46WOsW7eO1157jV69ejF27FhefPHFEkdvkpKS+M1vfnPCCM/MmTOLr99000387W9/Iy0tjQkTJvDwww/z7LPPAsFi9e233/L555/z8MMPM2PGjBL/fMpKhUdERMSCbt260bRp01Ifc/DgQRYsWMA111xT/LaCgoKTHjd37lw+/PBDAIYNG0adOnVOesxFF13E2LFjKSoqYtSoUaSmpp7yOZs0aXLKsgNw4YUX0qtXLwBGjx7N888/f0bTVXl5eeTm5pKWlgbAL3/5yxM+xyuvvBKALl26lFoGy0OFR0REYkspIzGVKTExsfj3Pp8Px3GKXz9y5AgAjuNQu3bt4hGYs9GnTx/mzp3LtGnTGDNmDL///e+56aabSs31c8aYU75+fP5j2c9GQkICAF6vF7/ff9YfD7SGR0QsSrp/Wqm/RKJFjRo1OHDgQInX69evz08//cSePXsoKCjgs88+A6BmzZo0bdqU9957DwguJF6xYsVJ79+nTx/eeustAKZPn168HuZ4mzZton79+tx6663ccsstLFu2DIC4uDiKiorK9Hls3ryZhQsXAsEpsksuuQQITl8tXboUgA8++OC0n3etWrWoU6dO8fqcN954o3i0p6Ko8IiIiFSwc845h169etG+fXvGjRt30vW4uDgmTJhAt27dGDhwIG3atCm+9uabb/Laa6+RkpJCcnIyU6dOPen9H3zwQebOnUtycjIffvghjRs3PukxGRkZpKSk0KlTJ9555x3uuusuILhGqGPHjtx4442n/Txat27N3//+d9q2bcu+ffu47bbbip//rrvuomvXrni93uLHjxgxgo8++qh40fLxXn/9dcaNG0fHjh1Zvnw5EyZMOO3znw1tSxcRa043ipMzcVglJRGRKKFt6SIiIhK7VHhEREQk6qnwiIiISNRT4REREZGop8IjIiIiUU+FR0RERKKeCo+IiEglMMYwevTo4tf9fj/nnnsuw4cPL37b9OnT6dq1K+3ataNTp07ce++9QPDUcWMM2dnZxY999tlnMcawZMkSIHgMxa9//WuaN29Oly5d6Nu3L4sWLaqkz+7UxowZw/vvv281wzE6WkJERGJOqO/kXZZ7RiUmJpKZmcnhw4epWrUq6enpXHDBBcXXMzMz+d3vfse0adNo06YNgUCAl19+ufh6hw4dmDJlCg888AAA7733HsnJycXXb7nlFpo2bUpWVhYej4cffviB1atXh/CzjGwa4REREakkQ4cOZdq0YNl6++23ueGGG4qvPfnkk4wfP774Lster7f4TsYAo0aNKr7L8oYNG6hVqxb16tUrfn3RokU8+uijeDzBH+1NmzZl2LATi1ggEGDMmDG0b9+eDh068MwzzwDwyiuvcNFFF5GSksJVV11Ffn4+EByhue222+jRowfNmjUjIyODsWPH0rZtW8aMGVP8catXr84999xDcnIy/fv3Z9euXSd97kuXLiUtLY0uXbowePBgtm/fDsDzzz9Pu3bt6NixI9dff/2Z/+GehgqPiFQs14XNi2DaH2DSJTCxCTxaH57tyEtxf+UabwbVOPvDBkUiwfXXX8+UKVM4cuQI33//Pd27dy++lpmZSZcuXUp835o1a3LhhReSmZnJlClTuO6664qvrVq1itTU1BOOdTiV5cuX8+OPP5KZmcnKlSu5+eabgeDp5IsXL2bFihW0bduW1157rfh99u3bx8KFC3nmmWcYOXIk99xzD6tWrWLlypXFh5oeOnSIrl27smrVKtLS0nj44YdPeN6ioiLuuOMO3n//fZYuXcrYsWMZP348ABMnTuS7777j+++/Z9KkSWX8kyw/TWmJSIW58o/PMCHuDVI9G8h3E1jitGKj241C4qhfsI8unvUM9i7h//ne4mX/cF4JDMV/3LclHT0h0aZjx47k5OTw9ttvM3To0HK//7HC9OWXXzJz5kwmT55crvdv1qwZGzdu5I477mDYsGEMGjQICJatBx54gNzcXA4ePMjgwYOL32fEiBEYY+jQoQP169enQ4cOACQnJ5OTk0Nqaioej6e4gI0ePZorr7zyhOddt24dmZmZDBw4EAiONDVs2LD4z+TGG29k1KhRjBo1qtx/JmWlwiMioecvhPQJfJjwf2x36zK+aCwfBS4hnyo/e6BLZ5PFb31TuS9uCiO9C/hd0R1scC845YcViQYjR47kD3/4AxkZGezZs6f47cnJySxdupSUlJQS33f48OGMGzeOrl27UrNmzRPed8WKFQQCgVJHeerUqcOKFSv48ssvmTRpEu+++y7/+Mc/GDNmDB9//DEpKSn885//JCMjo/h9EhISAPB4PMW/P/a63+8/5fMYc+KRVq7rkpycXHzS+vGmTZvG3Llz+fTTT/nLX/7CypUr8flCX080pSUioZW/F14fDov+j8n+wfQveJo3AwNOUXYADMvcVvx30ThuLfw99UwuH8dP4FLPd5UeW6SyjB07lgcffLB4pOSYcePG8dhjj7F+/XoAHMc5aYqnWrVqPPHEE8XTQcc0b96crl278uCDD3LsUPCcnJzi9ULH7N69G8dxuOqqq3j00UdZtmwZAAcOHKBhw4YUFRXx5ptvlvtzchyneDfWW2+9xSWXXHLC9datW7Nr167iwlNUVMSqVatwHIctW7Zw6aWX8sQTT5CXl8fBgwfL/fxloREeEQmdAzvgX6Ng70a4ejIP/zvh9O9zVLrTlcyCprwc/7+8Eve/3Fn0Oz53elRgWBE7GjVqxJ133nnS2zt27Mizzz7LDTfcQH5+PsaYE7asH1PSwt5XX32Ve++9lxYtWlC1alXq1avHU089dcJjfvzxR26++WYcxwHg8ccfB+CRRx6he/funHvuuXTv3p0DBw6U63NKTEzk22+/5dFHH+W8887jnXfeOeF6fHw877//PnfeeSd5eXn4/X7uvvtuWrVqxejRo8nLy8N1Xe68805q165drucuK3OsCZag1IsiIsWO5MHkobD3B/jFFGja54y2/iZymMnxT9LZZHFb0d2kO11LfKzW8IiEh+rVq1fYyEw5mZIuaEpLRM6evxDeGQ271sJ1b0DTPmf8oQ5RlTGF95HpNuX5uBfoaDaEMKiIxCoVHhE5e1/cBz/MhZEvQIv+Z/3h8qnCLYV/YJdbi1fj/5d65IUgpIhUlDAZ3SmVCo+InJ3MD2HJP6DnnZB6w+kfX0a7qcWviu6lJof4a9yLGJyQfWwRiT0qPCJy5vZuhE/uhEbdoP+EkH/4tW5jHvbfRB/vSm7zfhLyjy8isUOFR0TOjOPAx78Fjweufg28cRXyNG8H+vFpoAf3+D6gjdlcIc8hItFPhUdEzsyyf8LmBTD4MajduAKfyPCnopvJI5GJcS/j0dSWiJwBFR4RKb/92yD9weBurNQbK/zpcqnBw0U3kerZyM3eLyr8+UQk+qjwiEj5fTkeAoUw4jkwJd72IqQ+dS5mRqATv/e9x3nsq5TnFJHoocIjIuWz+RtY9SH0ugvqNqvEJzY87L8JHwHG+d45/cNFRI6jwiMiZec48MUfoUbDYOGpZFvc+kwODOEq7zySzQ+V/vwiErlUeESk7DLfh23LglvQ4xOtRPi7fxR7qcGEuDeg9KNxRESK6fBQESmbQBHM/gs06AAdT314YWU4QDWe9V/Fo3GTYcNMaDGgxMee7iwvncUlEjs0wiMiZbNiCuzLgb7/L3jvHYveCVzKVrcezPqLRnlEpExUeETk9AJFMPcpaJgKrS+znYYifDzvvyI4vbZe29RF5PRUeETk9Ja/BbmboO8fK20b+ul8GOgNdZKC02wa5RGR01DhEZHSBYpg3tNwfmdoNdh2mmJ+fNBnHOxYCdkzbccRkTCnwiMipcv8EHI3Q9p9YTO6U6zDtcEt8gues51ERMKcdmmJSMlcFxb+Deq1hpaDbKc5SdID6dzqvZTxB95ixB//xkq3Mm+EKCKRRCM8IlKyH+YEp4x6/s76zqySvB3ox363Kr/2fWY7ioiEsfD8DiYi4WHB3yDx3ODUUZg6SDXeDAzgMs8iGpudtuOISJjSlJaInNrO1ZA9Ay59AOKqnPIhp7uxX2X5h38IY73TGeudzkP+MbbjiEgY0giPiJzaNy+Crypc9N+2k5zWLuowzenBVd55JHLYdhwRCUMa4RGJYSWN0NTkIIsS3uHjQC/++OeFlZzqzLzhH8iVCV9zhfdr/h0YaDuOiIQZjfCIyEmu8s6jqimMqOLwnduClU4S/+VNB3QjQhE5kQqPiPyMy2jvDL5zWrDKTbIdphwM/woMorVnK93NWtthRCTMqPCIyAku9qymuWc7b/hLPoU8XH0auJhcN5H/8n1lO4qIhBkVHhE5wWhvOvvc6kxzetiOUm5HSODdQF+GeBZzLrm244hIGFHhEZFi57KPQZ6lvBdIo4B423HOyJTApfiMwxXeebajiEgYUeERkWJXer8mzgR4O9DPdpQzttE9n8VOK671zkGLl0XkGBUeETnK5WrvXJY4rfjBbWg7zFl5N9CXFp5tdDZZtqOISJhQ4RERAFLNBlp6fuS9QJrtKGft80B3DrkJXOOdYzuKiIQJFR4RAeAa7xwOu/FMC3S3HeWsHaIq0wI9GOFdSFWO2I4jImFAd1oWERIoZIR3IdOdbhykmu04IfFuII1rfXMY6vmWD5w+p3zM6c4Cy5k4rCKiiYgFGuEREQZ7llDT5EfFdNYxS9zWbHQacK0vw3YUEQkDKjwiwtXeOWx16/GN09Z2lBAyvBfoS3fPWhqbnbbDiIhlKjwiMa4Be7jEk8kHgT64UfYt4eNALwBGehZYTiIitkXXdzcRKbeR3gV4jMuHgUtsRwm57ZzDIqcNo7zz0T15RGKbCo9IjBvpXchypzmb3Aa2o1SIjwO9aOHZRrLZZDuKiFikwiMSw5qbH2nvyeGTQE/bUSrM54HuFLpeRnrn244iIhap8IjEsJHehTiu4bNA5B0UWlZ5VGeOk8pI70I8OLbjiIglKjwiscp1GeFZwEKnHT9Rx3aaCjU10JOGZi/dPGttRxERS1R4RGLV9uU08+zgEyd6p7OOmeF05qBbhcs9mtYSiVUqPCKxauX7FLpepgcusp2kwh0hgS+drgz1LiKeIttxRMQCFR6RWOQ4kPkhc5xU9lPddppK8WmgJ7VMPpd4VtqOIiIWqPCIxKLNC+HANj4JXGw7SaWZ77Qnz63GMO8i21FExAIVHpFYtPpj8FVhptPZdpJKU4SPrwJdGehZShx+23FEpJKp8IjEGseBNZ9CiwHkU8V2mkr1udOdmiafXp5M21FEpJKp8IjEmq2L4cB2aHe57SSVbr7Tnv1uVYZ6NK0lEmtUeERizZpPwBsPrQbbTlLpCokj3enCIO8SfJrWEokpKjwiscR1YfUn0LwfVKllO40V0wPdqW0O0dOzynYUEalEKjwisWTbd5C3GdqOtJ3EmnlOBw64VbnM863tKCJSiVR4RGLJ6qng8UHry2wnsaaAeGY6nRjsXYyXgO04IlJJVHhEYoXrBtfvNO0D1eraTmPV9EB36pqD9PCsth1FRCqJCo9IrNiZCXs3xuTurJ/LcFI45CYwTLu1RGKGCo9IrFg9FYwH2gy3ncS6AuKZ7aQywLsMg2M7johUAhUekVix5lNo0gsS69lOEha+CnTlPJNLqtlgO4qIVAKf7QAiUnGS7p8GQBOzgzkJa3loWzf+efRtsS7DSaXI9TLQu5Tv/C1txxGRCqYRHpEYMMCzFIAZThfLScLHfhJZ5LRhkGeJ7SgiUglUeERiwCDvUtY4jdnqnms7SlhJd7rSwrONpma77SgiUsFUeESiXG0O0NWsIz2GTkYvq/RAcMRroEZ5RKKeCo9IlOvn+Q6vcUkPdLUdJexsox6ZThIDvUttRxGRCqbCIxLlBniXscOtQ6abZDtKWEoPdKGLyeIc8mxHEZEKpMIjEsUSKCTNs4IZgc64+ud+SulOFzzGpb93me0oIlKB9B1QJIpd7FlNoinQ7qxSrHabsNWtx0CPprVEopkKj0gUG+hZykG3CgucZNtRwpghPdCF3p6VVOWI7TAiUkFUeESilePQ37uMuU5HComznSasfeV0pYopordnpe0oIlJBVHhEotX25TQw+4q3XkvJFjutyXOrMUi7tUSilgqPSLRa9zkB1zDbSbWdJOz58THL6URfz3IdJioSpVR4RKLVuukscVuTSw3bSSLC7EAn6pn9pJiNtqOISAVQ4RGJRvtyYGcmX2k6q8zmOB0JuIZ+2p4uEpVUeESi0fovAZip4yTKLI/qLHVb0c+z3HYUEakAKjwi0SgrHeo2J8dtaDtJRJkV6ER7Tw7nsc92FBEJMRUekWhTdBhy5kHLgbaTRJxZTicALvVqlEck2qjwiESbnPngPwItVHjKa73biK1uPfp7tI5HJNqo8IhEm+x08FWBpF62k0Qgw+xAKr08mcRTZDuMiISQCo9ItMlKh6TeEFfVdpKINMvpRKIpoLtnje0oIhJCKjwi0WTvRti7Qet3zsICJ5nDbjz9PN/ZjiIiIaTCIxJNsmYEX7YYYDdHBCsgngVOcrDwuK7tOCISIio8ItEkewbUbQbnNLedJKLNdlJp4vkJdmfZjiIiIaLCIxItio7AD3M1uhMCswLB7elkfWk3iIiEjAqPSLTYNB/8h7UdPQS2UY81zoXFd6wWkcinwiMSLbJngDcBki6xnSQqzHY6weaFcCTPdhQRCQGf7QAiEiJZ6cGyE1/NdpKoMCuQyu3OJ9z+56f53OlxysfkTBxWyalE5ExphEckGuzLgT1Z2o4eQt+5Lcl1E+mnYyZEooIKj0g0yEoPvtT6nZAJ4GWOk0KaZzkGx3YcETlLKjwi0SB7BtRJ0nb0EMsIpHCu2U+yybEdRUTOkgqPSKQr3o4+EIyxnSaqzHU6ApDm+d5yEhE5Wyo8IpFu80Ioytf9dyrAHmrxvdOUvlrHIxLxVHhEIl32DPDGQ9PetpNEpQwnhc4mi5octB1FRM6CCo9IpMtKhya9ID7RdpKoNCeQgte4XOLJtB1FRM6CCo9IJMvdDLvXaTt6BVrutiDPraZ1PCIRToVHJJJpO3qFC+BlntOBNO8KQKeni0QqFR6RSJY9A2o3hnotbSeJanOcFBqYfbQxW2xHEZEzpMIjEqn8BbBxjrajV4I5gRQA0jwrLCcRkTOlwiMSqTYvhKJDWr9TCX6iDqudJvRV4RGJWCo8IpEqKz24HT1J29ErwxynI10960jksO0oInIGVHhEIlX2TGh8MSRUt50kJsxxUogzAXppe7pIRFLhEYlEeVth1xpNZ1WipU4rDrhVtT1dJEL5bAcQkTNwdDv6gE/jyf5kmuUwsaEIHwuc5OD2dL8LaKG4SCTRCI9IJMqeAbUuJNu9wHaSmJLhpNDI7Ka52WY7ioiUkwqPSKTxF8LGjKOHhWqUoTId256u3VoikUeFRyTSbPkGCg9q/Y4F26jHeucC3Y9HJAKp8IhEmqx08MRB0z62k8SkOU4K3T1rqMoR21FEpBxUeEQiTfYMaHIxJNSwnSQmZTgpJBg/PTxrbEcRkXJQ4RGJJHlb4afVOizUosVOG/LdBE1riUQYFR6RSJI9M/iyxQC7OWJYIXEscNpp4bJIhFHhEYkk2elQ8wI4r63tJDFtjpNCkmcn7NlgO4qIlJEKj0ikCBQdPR19gE5HtyzDCW5PJ3uG3SAiUmYqPCKRYssiKNiv7ehhYItbn41OAxUekQiiwiMSKbLSweODpmm2kwjBaS1+mAdF2p4uEglUeEQiRfaM4OnoVWraTiIcLTz+w7Bpvu0oIlIGKjwikWD/NtiZqd1ZYeQbpy34qmhaSyRCqPCIRIJjP1S1fidsHCEBmvRS4RGJECo8IpEgKx1qnA/ntbOdRI7XciDsXg/7NtlOIiKn4bMdQEROI1AUPB293eXajh5m+k31MSsBxj/9LG8GTp5uzJk4zEIqETkVjfCIhLuti7UdPUxtdBuy2TlXd10WiQAqPCLh7th29GZ9bSeRkxjmOCn09GQST5HtMCJSChUekXCXnQ4XdocqtWwnkVPIcFJINAV08ay3HUVESqHCIxLODuyAHSu1HT2MLXSSKXB99PUstx1FREqhwiMSzrQdPezlU4XFTmut4xEJcyo8IuEsKx1qNIT67W0nkVJkOKm09mylIXtsRxGREqjwiISrgB82zoYW/bUdPczNOXp6eppXozwi4UqFRyRcbV0MR/Kghaazwl2WewE/uudoWkskjKnwiISr7HQwXm1HjwiGOYEUenky8eG3HUZETkGFRyRcZc+AC7tB1dq2k0gZzHFSqGEO08Vk2Y4iIqegwiMSjg7shO0rtB09gsx3kilyvVrHIxKmdJaWSDjaMBOAYdOrsurzaZbDSFkcpBpLjm5Pf5LrbccRkZ/RCI9IOMpKh+r1WeU2sZ1EymGO05F2nk2cxz7bUUTkZ1R4RMJNwA8bZh2dztJ29EiS4aQC2p4uEo5UeETCzY9L4Uiu1u9EoLXuhexw65Cm7ekiYUeFRyTcHNuO3vxS20mk3ILb03t7VuIlYDuMiBxHhUck3GR9dXQ7eh3bSeQMZDgp1DL5pJps21FE5DgqPCLh5Nh2dB0WGrHmO+3xux76ah2PSFhR4REJJ8Wnow+ym0PO2H4SWea2pK9nue0oInIcFR6RcJL1lU5HjwJzAil08OTAwZ9sRxGRo1R4RMJFwA8bZgd3Z+l09IiWcfT0dLJn2g0iIsVUeETCxdZvoSBP63eiwGq3CbvcWv+ZohQR61R4RMJF1lfg8el09Cjg4mGu0zF4RIij7eki4UCFRyRcZKVD44uhSi3bSSQEMgIpcHgf/LjMdhQRQYVHJDzs3wY7M3V35Sgyz+kAxqNpLZEwocIjEg6y0oMvtR09auRSAy7oErxztohYp8IjEg6y06FmIzivre0kEkotBgantA7ttp1EJOap8IjY5i+EDRnQUtvRo06LAYAbvN2AiFilwiNi25ZvoPCAprOi0fmdoNo5mtYSCQMqPCK2ZaWDJw6aptlOIqHm8UDz/sEbEDqO7TQiMU2FR8S2rHRo0hMSqttOIhWhxQDI3w3bvrOdRCSmqfCI2JS7BXat0XRWNGsxILg9PetL20lEYpoKj4hNx9Z26DiJ6JV4DjTqBuu/sJ1EJKap8IjYlJUOtRtDvVa2k0hFajUYtq8I3mBSRKxQ4RGxxV8AG+cEp7O0HT26tRoSfJn1ld0cIjFMhUfElk3zoehQ8OZ0Et3Oawu1GsN6reMRscW4rlva9VIvishZmH4fR755jdSClzlCgu00UsH+7JvM1d65dCp4iQLiAciZOMxyKpGoU+JwuUZ4RGxwXVg3na+d9io7MWKW04lqpoAenjW2o4jEJBUeERt2rYXcTcx0OttOIpVkodOOfDeBfp5ltqOIxCQVHhEb1k0HYFagk+UgUlkKiGe+057+3u/QagGRyqfCI2LD+i+hYQo7qWs7iVSimU4nGpndtDJbbUcRiTkqPCKV7dAe2PottLrMdhKpZLMDqQD09+iYCZHKpsIjUtmyvgLXCd6MTmLKTuqy0kmin1freEQqmwqPSGVbPx2qN4CGqbaTiAWznE50NlnUYb/tKCIxRYVHpDL5CyF7VnB0x6N/frFoVqATXuOS5vnedhSRmKLvuCKVadN8KDwArbV+J1Z97zZjl1uTfl6t4xGpTCo8IpVp/RfgqwJN02wnEUtcPMwOdCLNswICRbbjiMQMFR6RynL07so0TYP4arbTiEUznU7UMvmwZZHtKCIxQ4VHpLLsWge5m6D1ENtJxLKvnQ4Uut7iG1CKSMVT4RGpLOuP/nBrpcIT6w5RlW+cdsHCU/oBziISIio8IpVl3RfQoCPUPN92EgkDXzldYe8G2L3edhSRmKDCI1IZjt1dWaINU2YAABcXSURBVLuz5KgZgaMHx679zG4QkRihwiNSGYrvrqzpLAnawTlwfmdYO812FJGYoMIjUhnWfgY1zofzdTq6HKfNMPhxKezfZjuJSNRT4RGpaIX5kD0z+MPNGNtpJJy0GR58ue5zuzlEYoAKj0hF2zAL/Ieh7XDbSSTcnNsa6jbXtJZIJVDhEaloaz+DKrWhSS/bSSTcGBMc+fthHhzJs51GJKqp8IhUpIA/eK+VVkPAG2c7jYSjNsPAKYKsdNtJRKKaCo9IRdo0H47kajpLStboIkg8V9NaIhVMhUekIq2dBr6q0Ly/7SQSrjze4P2ZstLBX2A7jUjUUuERqSiuGyw8zfvpsFApXZvhUHgguJZHRCqECo9IRdn2HezfquksOb2maRCXqLsui1QgFR6RirJ2Ghiv7q4spxdXBVoOCN6Px3FspxGJSio8IhVl7WfQpCdUq2s7iUSCNsPh4E74cYntJCJRSYVHpCLszoZda6HtCNtJJFK0HASeOFg91XYSkaikwiNSEdZ8EnzZeqjdHBI5qtaG5pfC6k+CC95FJKSMW/o/LP2rEzkTL/UJ/m/91pkk3a/7q0jZXO2dw9NxLzGy4BG+d5uTM3GY7UgikabEAws1wiMSant/gO0rIHmU7SQSYdIDXShyvQz1LrIdRSTqqPCIhNrqj4Mv211uN4dEnDyqM99pz1DPIjTALhJaKjwiobbqY7igC9RubDuJRKDPnW409uyivfnBdhSRqKLCIxJKe3+A7cuhnaaz5Mx8FeiK3/Uw1Put7SgiUUWFRySUjk1naf2OnKFcarDASQ5Oa2m3lkjIqPCIhJKmsyQEPne6k+TZCTtW2o4iEjVUeERCRdNZEiLHprWKRwxF5Kyp8IiEyrE75Gp3lpylvdTkG6dtcMRQ01oiIaHCIxIqqz+G8ztDnSa2k0gUmO50h70bYOcq21FEooIKj0go7MuBbd9psbKEzBeBi8B4YdWHtqOIRAUVHpFQWPVR8KWmsyRE9lALmvWFle9pWkskBFR4REJh5fvQ6CKok2Q7iUSTDtdA7mbYuth2EpGI57MdQCTi7VwNOzPhsqdsJ5Eo036KjyUJcUx56Ske8o856boOFxUpO43wiJytzPeDay20fkdC7CDVmOF0Zrj3G7wEbMcRiWgqPCJnw3WDayya9YXq59lOI1Hok0BP6pn99PJk2o4iEtFUeETOxpZvg2ssOl5rO4lEqQwnlf1uNS73LrAdRSSiqfCInI2V74KvCrTRWgqpGIXE8XmgG4M9i0mg0HYckYilwiNypgJFwe3orS+DhBq200gUm+r0oro5Qn/PMttRRCKWCo/ImdqYAfl7oIOms6RiLXLastOtrWktkbOgwiNypr5/F6rUhhYDbCeRKOfg4dPAxfT1LKcmB23HEYlIKjwiZ6LgIKydFryzsi/edhqJAVMDvUgwfi7z6iaEImdChUfkTKz5BIoOQcoNtpNIjFjpNiXbOZ+rvXNsRxGJSCo8Imdi+VtQtxk07mE7icQMw/uBPlzkWU+S2W47jEjEUeERKa99OZAzD1J/AcbYTiMx5MNAbwKu4SrvPNtRRCKOCo9Iea2YAhjoeL3tJBJjfqIOc52OXOWdiwfHdhyRiKLCI1IejhOczmqWBrUvtJ1GYtD7gTTON3vp6VllO4pIRFHhESmPTfMhdxOkjradRGLUDKczuW6iFi+LlJMKj0h5LH8LEmrqKAmxpoB4Pgn0ZIhnMRzOtR1HJGL4bAcQiRgFB2H1VOhwNcRXK35z0v3TLIaSWPR+oA83+dKDR5t0vdl2HJGIoBEekbJaPTV4751Oms4Su753m7HOaQTL37QdRSRiqPCIlNWy1+GcltDoIttJJOYZ3gukwdbF8NMa22FEIoIKj0hZ7FwFWxZBlzG6946EhQ8CvcEbD0sm244iEhFUeETKYslk8CYEbzYoEgb2UTN4ltuKKVCYbzuOSNhT4RE5ncJD8P07kDwKqtW1nUbkP7rcDAV5sOpD20lEwp4Kj8jpZH4IBfuDP1xEwkmTnlCvtaa1RMrAuK5b2vVSL4rEhJcvhaJ8uP2bU67f0bZ0sWmsdzoT4t7gsoLHWeM2OeFazkTdL0piTomLLDXCI1Ka7Stg27Lg6I4WK0sY+iDQmyNuHL/wzrQdRSSsqfCIlGbJZPBVgZTrbCcROaU8qjPN6cEo73yqccR2HJGwpcIjUpKCA7DyPUi+EqrWsZ1GpERv+vtTwxzmcu9821FEwpYKj0hJlr8FhQfhov+2nUSkVMvclqx2mvBL71do6aXIqanwiJyK48Cil+CCrtCoq+00IqdhmBwYTBvPFi72rLYdRiQsqfCInEp2OuzdAD1us51EpEw+CfRkt1uTsd7ptqOIhCUVHpFT+eb/oEbD4J1sRSJAAfG8GehPf893NDY7bccRCTsqPCI/99Na2Dg7uHbHG2c7jUiZ/ds/AD8exni/tB1FJOyo8Ij83LcvBc/N0p2VJcLsog6fORdzjXcO1dH5WiLHU+EROd7hfcHDGDtcA4n1bKcRKbfJ/iHUMIe5xjvHdhSRsKLCI3K8xa8Gj5HQYmWJUCvdZixxWgWntZyA7TgiYUOFR+SYosPwzSRoMQAatLedRuSMvea/jCaen2DNp7ajiIQNFR6RY5a/Cfm7odfdtpOInJUvnYvY6DSAr5+B0g+IFokZPtsBRMJCwA8L/gYXdIGkS064pNPQJdI4eHgpMIIntr8S3HHYvJ/tSCLWqfCIAKyZCvtyYOAjOhVdosJHgUu4x/c+G/75ADcWjT/lY3ImDqvkVCL2aEpLxHVh/nNwTgtoox8AEh0KieNV/1B6eVeRYrJtxxGxToVHZMNM2L4Cet4JHq/tNCIh83agH7luIr/xafGyiAqPxDbXhYyJULMRpNxgO41ISB2iKq8HBjHYs4SWZqvtOCJWqfBIbNswE7Yuhj73gi/edhqRkJvsH8IhqnCX7wPbUUSsUuGR2HVsdKfWhZA62nYakQqRSw0mBwYz3LuINmaz7Tgi1qjwSOzKPjq601ujOxLdXvUPZb9blbs1yiMxTIVHYpPrQsbjR0d3brSdRqRC7ac6r/mHMsS7mGSTYzuOiBUqPBKbsr6CH5dA799rdEdiwj8Cl5HnVtMoj8QsFR6JPU4AZjwEdZpq7Y7EjANU42X/cAZ6l9LJZNmOI1LpjFv6OSs6hEWiz3dvwtTb4erJ0P7K0z5cR0tItKjGEeYk3MMPbgOuLZwAlH5Xcd2JWSJQiX+pNcIjsaXoMMz+C5zfGZKvsJ1GpFLlU4Vn/FfTzbOOQZ4ltuOIVCoVHokti16C/T/CwD/rzCyJSe8E+pLlXMB9vin48NuOI1JpVHgkduTvha//Ci0HQdPettOIWBHAy0T/9TT3bOd672zbcUQqjQqPxI6Mx6HgAAx4yHYSEatmOp1Z5LThLt8HVCffdhyRSuGzHUCkUuzIhMWvQtexUD/5hEtalCyxx/Bo0Wimxv+JO30f8Zhf96KS6KcRHol+rgvT/weq1IZLx9tOIxIWVrrNeCfQl5u9X+hgUYkJKjwS/TI/gE3zof+foFpd22lEwsaT/us4RBX+7PsnuguJRDsVHoluBQfhqwegYQp0/qXtNCJhZR81ecp/HRd7VzPCs9B2HJEKpcIj0W32Y3BgO1z2FHi8ttOIhJ23A/1Y6SQxPu5NEjlsO45IhVHhkei1dSks+r/gQuXG3W2nEQlLDh7+VDSW88jlf3xTbMcRqTAqPBKd/IXwye+gegMY8LDtNCJhbbnbgsmBIfzSl043s8Z2HJEKobO0JDrNeTJ4hMQNU6D1ZaU+VNvSRaAqR/gi/n4cDJcVTuQICad9H521JWFIZ2lJDPlpLcx9CpKvPG3ZEZGgw1Thfv+tNPXs5B7f+7bjiIScCo9EF38hfHgrJNSAy560nUYkoix0knnT359bvJ/T2ay3HUckpFR4JLrMfhR2fA8jX4Dq59pOIxJxHvffwDa3Hs/F/Z0aOnZCooiOlpDo8cNcmP88dBkDbYYWv1lrdETK7iDVuKvot7wb/2ceifsHdxf9llKWRYhEDI3wSHQ4vA8++g2c0xwGP2Y7jUhEW+a24jn/lYzyLuAKz9e244iEhAqPRD7HgY9ug4M74apXIT7RdiKRiPf3wCgWOW14JG4ySWa77TgiZ02FRyLf13+F9dNh8ONwfifbaUSigoOHewpvpwgfk+KepRpHbEcSOSsqPBLZsmfCrEehw7XQ7VbbaUSiyjbqcUfRHbQ0W3ky7mV0azaJZCo8ErlyN8MHt8B57WDEs2C0sFIk1L52OvCU/zqGe7/hV97PbMcROWMqPBKZjuTBW9eBE4Dr3tC6HZEKNCkwgmmBbtznm0IfzwrbcUTOiAqPRJ5AEbz7S9i9Hq77V3BnlohUIMO4ot+wzm3Mi3HP0c7k2A4kUm4qPBJZXBc+uwc2zoYRz0GzvrYTicSEfKpwc+E48khkcvyTnM9u25FEykWFRyLLnCfguzegzzjoNNp2GpGYspO6jCm8j6oU8s/4J4L3vxKJECo8EjnmPw8Zj0PKL+DS8bbTiMSkLLcRvyr6PU3MTvj31XBkv+1IImViXLfUbYbagyjhYdHLMH0cJF8BV74K3v+ciqKjI0Qq3wDPUl6t8hxc0AVGfxA8sFfEvhK366rwSPhb8o/gup3Ww+Da18Ebd8JlFR4RO4Z4vuWFuOdZ5rZkTOF95FPlhOs5E4dZSiYxrMTCoyktCW9fPxMsOy0HwTWTTyo7ImLPF0437i76LV3Mev4d/xi1OGg7kkiJVHgkPLkupE+AGQ9B+6vgujfBl2A7lYj8zGfOxdxedBfJJof34h+mAXtsRxI5JRUeCT/+Avj4dpj/HHT9b7jyFfDF204lIiX40unGmKL7aGj28kHCQ7QwW21HEjmJCo+El4O74PURsOItSLsfhv0veLy2U4nIaSx0krm+8E/E4+ej+Afp71lqO5LICVR4JHzsWAmvXArbV8DV/4BL/6jzsUQiyCo3iZEFj5Lj1ueVuL/C3KeD09MiYUCFR+xzXVj8KrzSHxw/3Dw9uG5HRCLOds7h6sKHmOr0hFmPwJQbIX+v7Vgi2pYulh3eB5/cAWs+hRYDYNQkqH7uCQ/RtnORSOSSc/nm4OaDxHpwxUvQLM12KIl+2pYuYWjddHjx4uDLgY/AL947qeyISKQycPHtcOtMiK8O/7ocvhwPhYdsB5MYpcIjle/Qbnh/LLx9PVStC/+dDr3uBI/+OopEnYYp8Os50PVmWPhC8D852TNtp5IYpCktqTz+wuBanTkToTAf0v4Het1N0gPptpOJSAU46U7LOfPh07tgTxYkXwkDHoI6TWxEk+hV4pSWr6QLIiHjusFpq/Q/wZ5saN4PBj8O57WxnUxEKlNSL/jN18E7qM9/DtZ+Bt1/A73vhaq1baeTKKcRHqk4rgtrp8GcJ2DH93BOCxj8WPCYiOO2m2tRskjsacAevum+AFa8DQm1oPuvocdtUK1umT/G6b536CyvmKQRHqlERUdg9cew4AXYuRLqNoPLX4SO1+osLBEBYAfnkLRoBMmmA78LfMRlc5/k4Jzn+XdgAK/7B7Odc1RYJKQ0wiOhk7sFlk6Gpa9D/m6o1wp630vzt6sRQHdLFpGStTJbuN03lRGehQDMdDrzRmAgXzvtcc9wf40KU0wqcYRHhUfOzuFcWPMJfP8u5HwdnKpqPRS63QpN08AYTVmJSJldwC5u8M3ieu9s6pn9bHHO5VPnYqYGerLObRzS51IhikoqPBJCeVth/ZfBXxszIFAQXJ/T4VpIvQFqn/hNSYVHRMorniKGeBZzpXcel3hW4jMO65xGfOV0ZXYgleVuC5yzvLOKCk9UUuGRs5C/FzZ/A5sXBAvOjpXBt9dJCo7mdLgazu9c4rlXKjwicjbOIY/LvN8ywruQLmY9PuOwz63OHKcjC51kljit2OCeTyk/605JhScqqfBIGRUdhp9WB0vN9hXBovPT6uA1bzw0ughaDYZWQ4JrdMpwuKcKj4iESk0O0sezkku9y0nzrKCe2Q/AHrcGS51WLHVasdptwhqnCbupdVbPpUIUkVR4olm5t2a6LuTvgb0b//NrTzbsXAW714PrBB+XUBMadYUmPaFxT7igC8RVKffzi4hUDJdmZjtdPeu4yKyjq2cdTT07i6/ucmux2mlCtnsBP7gN2OTWJ8etz4/uuWc9HQYqRGFKhSechaow+PBTk3xqmUPU4QDnmlzqm338uV89OLDjP7/ytkDB/uPe00DtC+G8ZGjQ4T+/ajcp03EPKjwiEi5qc4C2ns20M5to69lMW7OJpmYH1UxB8WMKXS873LrsoC4/uXXY4dZhhxv8/S5qkecmkutWJ5fq5JNAeafKjlEhsiJGC4/rBn8FX/nP2469Xtq14tdLu1bKx3EC4BSB46fPEzOIw48XBx8BvDj/ed0Ejr7dj+/o9SoUUsUUUpUCqlJIgimkKsHXq1BIVVNIFQqoYQ5Ti0PUNIeoxSESj/sHfTy/62EXtfnJrc1Pbh22uXXJcRuQc/R/PLMeuQl8CWfwBxykwiMi4c3lPHJJMjtI8uygqdlBQ7OHBmYf57GPBmbfCYXoeAWujzyqk+smcoBq5LsJHCaBfBLIdxPIpwr5JHDYDb4sOPrdvdD18fyN3YJLAbzxwXuQeeOO+308eOKC/6k0p/rlPfrSnHzN4/3Z286skEWpMC08OzLhtUGEvIxEoeA/sniOEM8RN/jyANXIcxPZ71Yjj0Ty3ETySGS/m0guiexy67DTrcNeaoRk+FZEJDq51OAw9c1e6pn91OIgtc0hanOQ2uZg8es1OUQ1U0BVCqhGwQm/95hw+vnzs5/5JxUiU7ZrIX1fA11+CUMePyltiJ1Z4THGfAHUq4hEUaoesNt2CDmJvi7hR1+T8KSvS/jR16R8druuO+RUF043wiPlYIxZ4rpuV9s55ET6uoQffU3Ck74u4Udfk9DRPIeIiIhEPRUeERERiXoqPKH1su0Ackr6uoQffU3Ck74u4UdfkxDRGh4RERGJehrhERERkainwlNBjDH3GmNcY4y29VtmjHnKGLPWGPO9MeYjY0xt25limTFmiDFmnTEm2xhzv+08sc4Yc6ExZrYxZrUxZpUx5i7bmSTIGOM1xnxnjPnMdpZooMJTAYwxFwKDgM22swgA6UB713U7AuuBP1rOE7OMMV7g78BlQDvgBmNMO7upYp4fuNd13XZAD+C3+pqEjbuANbZDRAsVnorxDPA/RPOtnyOI67pfua7rP/rqN0Ajm3liXDcg23Xdja7rFgJTgMstZ4pprutud1132dHfHyD4A/YCu6nEGNMIGAa8ajtLtFDhCTFjzOXAj67rrrCdRU5pLDDddogYdgGw5bjXt6IfrmHDGJMEdAIW2U0iwLME/+Ps2A4SLXy2A0QiY8wMoMEpLo0H/h/B6SypRKV9TVzXnXr0MeMJDt+/WZnZRCKBMaY68AFwt+u6+23niWXGmOHAT67rLjXG9LWdJ1qo8JwB13UHnOrtxpgOQFNghQkemtYIWGaM6ea67o5KjBhzSvqaHGOMGQMMB/q7uheDTT8CFx73eqOjbxOLjDFxBMvOm67rfmg7j9ALGGmMGQpUAWoaY/7tuu5oy7kimu7DU4GMMTlAV9d1dfCbRcaYIcBfgTTXdXfZzhPLjDE+ggvH+xMsOouBX7iuu8pqsBhmgv87ex3Y67ru3bbzyImOjvD8wXXd4bazRDqt4ZFY8AJQA0g3xiw3xkyyHShWHV08/jvgS4KLY99V2bGuF/BfQL+j/z6WHx1ZEIkqGuERERGRqKcRHhEREYl6KjwiIiIS9VR4REREJOqp8IiIiEjUU+ERERGRqKfCIyIiIlFPhUdERESingqPiIiIRL3/D6G6c88QmL3TAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 720x432 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "sample_and_display(np.array([2.0]), 0.1, 500000, 25000, log_prob)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Sloooowly getting there..." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Conclusions" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "With these considerations, I conclude the first blog post of this series.\n", | |
| "I hope you now understand the intuition behind the Metropolis-Hastings algorithm, its parameters and why it is an extremely useful tool to sample from non-standard probability distributions you might encounter out there in the wild. \n", | |
| "I highly encourage you to play around with the code in this notebook - this way you will learn how the algorithm behaves in various circumstances and deepen your understanding of it.\n", | |
| "Go ahead and try out a non-symmetric proposal distribution!\n", | |
| "What happens if you don't adjust the acceptance criterion accordingly?\n", | |
| "What happens if you try to sample from a bimodal distribution?\n", | |
| "Can you think of a way to automatically tune the stepsize?\n", | |
| "What are pitfalls here?\n", | |
| "Try out and discover yourself! \n", | |
| "In the next post, I will discuss the Gibbs sampler - a special case of the Metropolis-Hastings algorithm which allows you to approximately sample from a multivariate distribution by sampling from the conditional distributions." | |
| ] | |
| } | |
| ], | |
| "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.3" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 2 | |
| } |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/usr/bin/env python3 | |
| import sys | |
| from itertools import cycle | |
| import re | |
| with open(sys.argv[1]) as ipf: | |
| lines = ipf.readlines() | |
| # ## replace \{ with \\{ and \} with \\} | |
| lines = [l.replace('\\{', r'\\{') for l in lines] | |
| lines = [l.replace('\\}', r'\\}') for l in lines] | |
| ## replace \\ with \\\\ | |
| lines = [l.replace(r' \\', r' \\\\') for l in lines] | |
| ## replace ^* with ^\* | |
| lines = [l.replace(r'^*', r'^\*') for l in lines] | |
| ## alternatingly replace $ with \\( and \\) | |
| ## if it's not part of $$ | |
| lines2 = [] | |
| for line in lines: | |
| if '$$' in line: | |
| lines2.append(line) | |
| continue | |
| else: | |
| cycler = cycle((True, False)) | |
| matches = re.finditer('\$', line) | |
| offset = 0 | |
| for match in matches: | |
| replacement = '\\\(' if next(cycler) else '\\\)' | |
| line = line[:match.start()+offset] + replacement + line[match.start()+1+offset:] | |
| offset += 2 | |
| lines2.append(line) | |
| with open(sys.argv[2]) as ipf: | |
| header = ipf.readlines() | |
| with open(sys.argv[3], 'w') as opf: | |
| for line in header + lines2[2:]: | |
| opf.write(line) |
I don't think the method you use for the mixture of Gaussians can be called Gibbs sampling. As you say, in Gibbs sampling you'd have to draw from p(k|x), but it does not equal p(k) because k and x are not independent. You'd have to use the Bayes theorem to calculate p(k|x).
Your method is to draw hierarchically: first draw k from the marginal distribution, then draw x from the conditional distribution, whereas in Gibbs sampling you draw from both conditional distributions, each time using the value from the previous step.
Your method works, but it works for a different reason than the reason Gibbs sampling works. I don't think it has a name.
Edit: moreover, I think your method works well precisely because you are not using GIbbs sampling. If you did, you'd have the same problem of being stuck in one of the mixture components, because p(k|x) would very likely give you the same k as the one that generated that x.
Finished reading — very nice! I corrected a few typos here: https://gist.github.com/feuerbach/5ccaeee166b45be13a8e375f980f405a.
Maybe one day I'll make a similar post except everything is done in Stan :)
Thanks for your feedback, @MMesch and @feuerbach! Highly appreciated!
For now I mostly worked on the MH part and put it into a separate notebook, where I addressed some of @MMesch's points.
@MMesch: I fear the \pi is indeed somewhat standard - as p often denotes other probabilities such as p_acc.
@feuerbach: This is very interesting - and you're right. The background to this example is that, a few years ago, I did something like that to sample from a different kind of mixture model and back then we used to call it Gibbs sampling. Looking this up, it seems like a very common method to sample from Gaussian (and other) mixture models (where you can easily marginalize out x to obtain p(k)). It seems to be called "Collapsed Gibbs sampling". I will rewrite this part to feature an actual Gibbs sampler. I just need to think of a nice, concise example where you can sample from the conditional distributions without resorting to MCMC. And as for being stuck in one of the mixture components: highly correlated samples often are a problem when using Gibbs sampling and I should discuss this.
I am also thinking of discussing how to assess convergence for MCMC methods, but this is a very difficult topic. It would make the series much more complete, but is also a really ugly rock to look under...
@feuerbach: Just FYI, I implemented the actual Gibbs sampler for this problem and it turns out that
you'd have the same problem of being stuck in one of the mixture components, because p(k|x) would very likely give you the same k as the one that generated that x.
is a slight understatement. If the sampler starts in the mode on the right, there's (on average) a ~1/1e6 chance for it to jump to the mode on the left. That was fun and very instructive and might even serve in a next version of that blog post / notebook as a negative example.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Reads nicely @simeoncarstens !
Here are a few thoughts that I had while reading the beginning (will continue another day). Some of them are resolved by reading it several times, but maybe it is a good hint for where a reader might hang a bit:
what is y in the equation?
Calling the stationary distribution eigenvector
\piconfuses me because I associate it so much with a number. But maybe it is the standard. In that case never mind ...I feel this is a nice paragraph but it could be easier to understand without prior knowledge: You don't really introduce what "design our transition kernel" means. Also, how can I "draw from \pi" I thought I only draw a new state based on the row in T of the previous step?
the
^*is difficult to see and read. Probably better if it has different support. Otherwise, can you use a different symbol?how can we know when it is set correctly? Is this difficult to grasp intuitively?