Skip to content

/blog-resources

Permalink
Branch: master
Find file Copy path

blog-resources/mcmc-intro/mcmc_introduction.ipynb

Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
633 lines (633 sloc) 124 KB
Raw Blame History
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Markov chain Monte Carlo (MCMC) sampling, part 1: the basics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Markov chain Monte Carlo (MCMC)](https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) is a powerful class of methods to sample from probability distributions known only up to an (unknown) normalization constant. But before we dive into MCMC, let's consider why you might want to do sampling in the first place.\n",
"\n",
"The answer to that is: whenever you're either interested in the samples themselves (for example, inferring unknown parameters in Bayesian inference) or you need them to approximate expected values of functions w.r.t. to a probability distribution (for example, calculating thermodynamic quantities from the distribution of microstates in statistical physics). Sometimes, only the mode of a probability distribution is of primary interest. In this case, it's obtained by numerical optimization so full sampling is not necessary.\n",
"\n",
"It turns out that sampling from any but the most basic probability distributions is a difficult task. [Inverse transform sampling](https://en.wikipedia.org/wiki/Inverse_transform_sampling) is an elementary method to sample from probability distributions, but requires the cumulative distribution function, which in turn requires knowledge of the, generally unknown, normalization constant. Now in principle, you could just obtain the normalization constant by numerical integration, but this quickly gets infeasible with an increasing number of dimensions. [Rejection sampling](https://en.wikipedia.org/wiki/Rejection_sampling) does not require a normalized distribution, but efficiently implementing it requires a good deal of knowledge about the distribution of interest, and it suffers strongly from the curse of dimension, meaning that its efficiency decreases rapidly with an increasing number of variables. 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",
"\n",
"MCMC algorithms are a class of methods which do exactly that. These methods date back to a [seminal paper by Metropolis et al.](https://pdfs.semanticscholar.org/7b3d/c9438227f747e770a6fb6d7d7c01d98725d6.pdf), 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. In reality, they were looking for a general method to calculate expected values occurring in statistical physics.\n",
"\n",
"In this blog post, I introduce the basics of MCMC sampling; in subsequent posts I'll cover several important, increasingly complex and powerful MCMC algorithms, which all address different difficulties one frequently faces when using the Metropolis-Hastings algorithm. Along the way, you will gain a solid understanding of these challenges and how to address them. Also, this serves as a reference for MCMC methods in the context of the [monad-bayes](https://www.tweag.io/posts/2019-09-20-monad-bayes-1.html) series. Furthermore, I hope the provided notebooks will not only spark your interest in exploring the behavior of MCMC algorithms for various parameters/probability distributions, but also serve as a basis for implementing and understanding useful extensions of the basic versions of the algorithms I present."
]
},
{
"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",
"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.\n",
"In order to sample from a distribution $\\pi(x)$, a MCMC algorithm constructs and simulates a Markov chain whose stationary distribution is $\\pi(x)$, meaning that, after an initial \"burn-in\" phase, the states of that Markov chain are distributed according to $\\pi(x)$. We thus just have to store the states to obtain samples from $\\pi(x)$."
]
},
{
"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. 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 rain 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 that's lots of 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 continous, in which case $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 the transition kernel $T(x_{i+1}|x_i)$ into two parts:\n",
"a proposal step and an acceptance/rejection step. The proposal step features a proposal distribution $q(x_{i+1}|x_i)$, from which we can sample possible next states of the chain. In addition to being able to sample from it, we can choose this distribution arbitrarily. But, one should strive to design it such that samples from it are both as little correlated with the current state as possible and have a good chance of being accepted in the acceptance step. Said acceptance/rejection step is the second part of the transition kernel and corrects for the error introduced by proposal states drawn from $q \\neq \\pi$. It involves calculating an acceptance probability $p_\\mathrm{acc}(x_{i+1}|x_i)$ and accepting the proposal $x_{i+1}$ with that probability as the next state in the chain. Drawing the next state $x_{i+1}$ from $T(x_{i+1}|x_i)$ is then done as follows: first, a proposal state $x_{i+1}$ is drawn from $q(x_{i+1}|x_i)$. It is then accepted as the next state with probability \n",
"$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",
"\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, in the physics literature, *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",
"Now here's where the magic happens: we know $\\pi$ only up to a constant, but it doesn't matter, because that unknown constant cancels out in the expression for $p_\\mathrm{acc}$! It is this property of $p_\\mathrm{acc}$ which makes algorithms based on Metropolis-Hastings work for unnormalized distributions. 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",
"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",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Implementing the Metropolis-Hastings algorithm in Python"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All right, now that we know 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",
"Generally, 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": 12,
"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": 13,
"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": [
"## Testing our Metropolis-Hastings implementation and exploring its behavior"
]
},
{
"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": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Acceptance rate: 0.719\n",
"Last ten states of chain: 1.60194, 0.14617, -1.13052, 0.24360, -0.77219, -0.59475, 0.67935, 1.47293, 1.27589, 0.67263\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": [
"All right.\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": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAFlCAYAAADvZjI4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXiU1cH+8e95ZpKwhFUWkcWwrwkgq0IAURBBca9SeK1S1FqttlVbWluX1lrspvW1lrfV/uqOxV0RZJewKGsgAWQJhE0UEAhLCMk8z/P7YyAFIUAgkzPL/bmuXFlmJrnDXAl3zjnPOcb3fURERETimWM7gIiIiEikqfCIiIhI3FPhERERkbinwiMiIiJxT4VHRERE4p4Kj4iIiMS94Glu1zXrIiIiEitMWTdohEdERETingqPiIiIxD0VHhEREYl7KjwiIiIS91R4REREJO6p8IiIiEjcU+ERERGRuKfCIyIiInFPhUdERCTC9u7dy/PPP19pXy8tLY1du3YBcMkll5zyvk8++eQpbx86dCh79+4lPz+fTp06lSvH7NmzmT9/fun748eP5+WXXy7X56goxvdPuZmydloWERE5R/n5+Vx11VXk5uaecFsoFCIYPN3BB+WTlpbG4sWLqVev3mnvm5qayoEDB074uO/7+L6P44THRk71PZTlscceIzU1lQcffPDMw58b7bQsIiJiy9ixY8nLy6NLly489NBDzJ49m8zMTIYPH06HDh1OGD3505/+xGOPPQZAXl4eQ4YMoVu3bmRmZvLFF1+c8Pm/+eYbBg8eTMeOHRkzZgzHDmakpqYCsH37dvr160eXLl3o1KkTWVlZjB07lkOHDtGlSxdGjhxJfn4+bdu25dZbb6VTp05s2bLluNGiUCjEyJEjad++PTfeeCOFhYXA8SNKixcvZsCAAeTn5zN+/HiefvppunTpQlZWFo899hh/+tOfAMjOzqZ3795kZGRw3XXXsWfPHgAGDBjAz3/+c3r27EmbNm3IysqqkOegYiuliIhItJs8Fr7KqdjPeX46XDmuzJvHjRtHbm4u2dnZQHiqZ+nSpeTm5tK8eXPy8/PLfOydd97J+PHjad26NZ9//jk//OEPmTlz5nH3efzxx+nbty+PPPIIkyZN4sUXXzzh87z++utcccUVPPzww7iuS2FhIZmZmTz33HOlufLz81m3bh0vvfQSvXv3PuFzrFmzhhdffJE+ffowevRonn/++TJHb9LS0vjBD35w3AjPjBkzSm+/9dZb+d///V/69+/PI488wuOPP84zzzwDhIvVwoUL+fjjj3n88ceZPn16mf8+Z0qFR0RExIKePXvSvHnzU97nwIEDzJ8/n5tuuqn0Y4cPHz7hfnPmzOGdd94BYNiwYdSpU+eE+/To0YPRo0dTUlLCtddeS5cuXU76NS+88MKTlh2Apk2b0qdPHwBGjRrFs88+e1bTVQUFBezdu5f+/fsD8L3vfe+47/H6668HoFu3bqcsg+WhwiMilWPXesifA3s2gReCanWhcXdo2guSqthOJ4nkFCMxlal69eqlbweDQTzPK32/qKgIAM/zqF27dukIzLno168fc+bMYdKkSdx222389Kc/5dZbbz1lrm8zxpz0/WPzH81+LlJSUgAIBAKEQqFz/nygNTwiEkmeByvfg/GZ8Fw3+OgnsOBvsPhfMOM38PJw+Es7mPYoFO62nVYkYmrUqMH+/fvLvL1hw4bs2LGDb775hsOHD/PRRx8BULNmTZo3b87EiROB8ELi5cuXn/D4fv368frrrwMwefLk0vUwx9q0aRMNGzbkjjvuYMyYMSxduhSApKQkSkpKzuj72Lx5MwsWLADCU2R9+/YFwtNXS5YsAeDtt98+7fddq1Yt6tSpU7o+55VXXikd7YkUFR4RiYzdG8OFZuL3oOQQDHkK7lsGv9oBD2+HsZthxJvQvB/Mfxae6w4r/mM7tUhEnHfeefTp04dOnTrx0EMPnXB7UlISjzzyCD179mTQoEG0a9eu9LbXXnuNF198kc6dO9OxY0fef//9Ex7/6KOPMmfOHDp27Mg777xDs2bNTrjP7Nmz6dy5M127duXNN9/k/vvvB8JrhDIyMhg5cuRpv4+2bdvyt7/9jfbt27Nnzx7uvvvu0q9///330717dwKBQOn9r776at59993SRcvHeumll3jooYfIyMggOzubRx555LRf/1zosnQRqXh5s+A/3wN8GPwEdB0FTqDs+3+9Mjz6s+Vz6P59GDIOgsmVFldE4kaZl6Wr8IhIxVrxH3j3B1C/HYx4A+pceGaPc0Mw4/HwaE+bK+E7L0EwJbJZRSTeqPCISCXIfQfe/j5c2AdueR2q1Cz/51j0Akx6AFoPDn+OQFLF5xSReKWNB0UkwjYtgHfuDF919d03z67sAPQYA1c9A+umhqe5Tv1HmYjIGdFl6SJy7gq2wn/+B2o3DU9jJZd9WesZ6X477PsS5vwBGrSHi++pmJwikrA0wiMi5yZUDBNGQkkRjJgAVU/c8OysXPpLaHdV+JL1rUsq5nOKSMJS4RGRc/PpONieDdf9Heq3rbjPawxc8xzUaARv3Q6H9lbc5xaRhKPCIyJnb8simPs0dBkF7a+u+M9ftQ7c+K/wlNnUX1X85xepRMYYRo0aVfp+KBSifv36XHXVVaUfmzx5Mt27d6dDhw507dqVBx54AAifOm6MYf369aX3feaZZzDGsHjxYiB8DMVdd91Fy5Yt6datGwMGDODzzz+vpO/u5G677TbeeustqxmO0hoeETk7JUXw7l1QswkM+X3kvk7THnDJj2DeM5B+E7SI7G6skhjSxk6q0M+XP27Yae9TvXp1cnNzOXToEFWrVmXatGk0bty49Pbc3FzuvfdeJk2aRLt27XBdl3/84x+lt6enpzNhwgR+9atw+Z84cSIdO3YsvX3MmDE0b96cdevW4TgOGzduZNWqVRX4XcY2jfCIyNmZ/yzszoPhfy29Iitt7KQyX87JgLFQpzl8eH9412aRGDV06FAmTQr/PLzxxhuMGDGi9LY//OEPPPzww6W7LAcCgdKdjAGuvfba0l2W8/LyqFWrFvXq1St9//PPP+eJJ57AccL/tTdv3pxhw44vYq7rctttt9GpUyfS09N5+umnAfjnP/9Jjx496Ny5MzfccAOFhYVAeITm7rvvpnfv3rRo0YLZs2czevRo2rdvz2233Vb6eVNTU/nJT35Cx44dueyyy9i5c+cJ3/uSJUvo378/3bp144orrmD79u0APPvss3To0IGMjAxuueWWs//HPQ0VHhEpvz2bIOvP0OEaaDkw8l8vqSpc/Qzs2RieQhOJUbfccgsTJkygqKiIFStW0KtXr9LbcnNz6datW5mPrVmzJk2bNiU3N5cJEyZw8803l962cuVKunTpctyxDieTnZ3Ntm3byM3NJScnh9tvvx0In06+aNEili9fTvv27XnxxRdLH7Nnzx4WLFjA008/zfDhw/nJT37CypUrycnJKT3U9ODBg3Tv3p2VK1fSv39/Hn/88eO+bklJCT/60Y946623WLJkCaNHj+bhhx8GYNy4cSxbtowVK1Ywfvz4M/yXLD8VHhEpv09+CcaBK56svK/ZYgB0uBbmPRu+ZF0kBmVkZJCfn88bb7zB0KFDy/34o4Xpvffe47rrriv341u0aMGGDRv40Y9+xJQpU6hZMzw6m5ubS2ZmJunp6bz22musXLmy9DFXX301xhjS09Np2LAh6enpOI5Dx44dyc/PB8BxnNICNmrUKObOnXvc112zZg25ubkMGjSILl268MQTT7B169bSf5ORI0fy6quvEgxGbqWNCo+IlM+m+fDFR5D5ANRqUrlf+/LHwHdh5hOV+3VFKtDw4cN58MEHj5vOAujYsWPpieNlueqqq3jllVdo1qxZaVk5+tjly5fjuu4pH1+nTh2WL1/OgAEDGD9+PGPGjAHCU1fPPfccOTk5PProoxQVFZU+JiUlfMSL4zilbx99PxQKnfTrGHP8hse+79OxY0eys7PJzs4mJyeHqVOnAjBp0iTuueceli5dSo8ePcr8nOdKhUdEzpzvw/THw5eK9/5h5X/9us2h112Q/TpsX175X1+kAowePZpHH32U9PT04z7+0EMP8eSTT7J27VoAPM87YYqnWrVqPPXUU6XTQUe1bNmS7t278+ijj3L0yKj8/PzS9UJH7dq1C8/zuOGGG3jiiSdYunQpAPv376dRo0aUlJTw2muvlft78jyv9Gqs119/nb59+x53e9u2bdm5cycLFiwAwlNcK1euxPM8tmzZwqWXXspTTz1FQUEBBw4cKPfXPxO6SktEztzaT2DLZ3DV05BcrVwPPdXC5TO5wqVU5oOw7FWY+TsY+Z9yZRCJBk2aNOG+++474eMZGRk888wzjBgxgsLCQowxx12yflRZC3tfeOEFHnjgAVq1akXVqlWpV68ef/zjH4+7z7Zt27j99tvxPA+A3/8+fIXlb3/7W3r16kX9+vXp1asX+/fvL9f3VL16dRYuXMgTTzxBgwYNePPNN4+7PTk5mbfeeov77ruPgoICQqEQP/7xj2nTpg2jRo2ioKAA3/e57777qF27drm+9pnS4aEicmY8D/4vE0oK4Z6FJz3U82yvxjpV4TnZ5/xh4D1+lvQfuGMmNC57kaeIVI7U1NSIjcyUkw4PFZFztGYSfJ0LA35p/QTzl9wr2OOnwuxxVnOISOxQ4RGR0/N9yPpLeC+cjuW/MqSiHaQq/wwNC5+ornO2RKyLktGdU1LhEZHT2/gpfLkU+twPgehY+veSOzh89MTcv9iOIiIxQIVHRE4v6y+Q2hA6jzj9fSvJQapC9+/DF5PgmzzbcUQkykXHn2oiEr22LQmP8Az6DSRViciXOOujJ3reGT7i4rPnYdifKzaUiMQVjfCIyKnNexaq1IJut9tOcqIaDSHjO7DsNSjcbTuNiEQxFR4RKVvBNlj9IVx0a+kBoVHn4nshdAgWvXj6+4pIwlLhEZGyLfl/4HvhtTLRqkF7aDUIFv4fhA7bTiMiUUqFR0ROLnQYlvwb2gwJH+kQzXr/AA7uDI9GiYichAqPiJzcyvfCJaLXnbaTnF6LgeE9gjStJSJlUOERkZNb+H9wXmtoPsB2ktNzHOh+O2yeD1+vsp1GRKKQCo+InOjLZeHL0XveES4TsaDLKAgkw+J/2U4iIlEoRn6TiUilWvYqBKtAxs22k5y56udBh2th+QQ4HP3b3ItI5VLhEZHjlRyCFROh/XCoWtt2mvLp8X0o3g+5b9lOIiJRRoVHRI63+kM4XABdR9lOUn5Ne0H99uERKhGRY6jwiMjxlr0CtS+EtEzbScrPGOg6ErYugp1rbacRkSiiwiMi/7V7I2ycEx7diZXFyt+W/h0wAVj+uu0kIhJFYvQ3mohERPbrgIEu37Wd5OzVaAitB4UXL3uu7TQiEiVUeEQkzHMh+zVodRnUamI7zbnp8l3Yvx02zLKdRESihAqPiIRtnAP7tkGXkbaTnLs2Q6BqnSMjViIiKjwiclTOREipCW2vtJ3k3AVTIP0mWP0RHNprO42IRAEVHhEJ772z6gNofzUkVbWdpmJk3ALuYfjiI9tJRCQKqPCICKybGt6wL/1G20kqTuOLwgeK5mgTQhGBoO0AIhIFciZC9QaQ1u+0d00bO6kSAlUAY6DTDTD3L3BgB6Q2sJ1IRCxS4RFJdIf2wtqp4dPGA+FfCTFTak4n/UbI+hOsfA963Wk7jYhYpCktkUT3xUfhtS7pN9lOUvEatIcGHXW2loio8IgkvJyJ4bUujbvZThIZ6TfAls9hzybbSUTEIk1piSSy/V+F99/JfDC85iWOHJ2Wa2LqMjcFnvrTk/zdHQ5A/rhhNqOJiAUa4RFJZKs/BN+DTtfbThIxW/0GLPVaMTww33YUEbFIhUckka16H+q1Da91iWMfuJfQ3tlMS7PNdhQRsUSFRyRRHdwFm+ZBh+G2k0TcZLcnAEOcRZaTiIgtWsMjkqDGPvl7xiV5XDmtLqunxsll6GX4mros9VoxJLCQv7nX2o4jIhao8IgkqKHO5+R7DVntN7Md5ayVZ7+gyW5PHk56nSZmRwQTiUi00pSWSCI6tIeLnVVM8XoC8XV1VlmmeD0AuELTWiIJSYVHJBGtmUyScfn4yNqWRLDFb8hK70KuDKjwiCQiFR6RRLTqfbb69Vjht7CdpFJNcXvQ3Vkb3n9IRBKKCo9IoinaB3kzmeL2IFGms46a7B0Z0Vr9od0gIlLpVHhEEs3aT8AtLr1UO5Gs9xuT5zWC1R/YjiIilUyFRyTRrH4fUs9nqd/adhILTHiUJ38eHPzGdhgRqUQqPCKJpKQI1s+AdsPwE/THf4rbA3wX1nxsO4qIVKLE/I0nkqg2zoGSQmg71HYSa3L95lCrmdbxiCQYFR6RRLLmY0hOheaZtpNYZKDdUNj4KRQX2g4jIpVEhUckUfg+rJ0CLQdCMMV2GrvaDIFQEWyYbTuJiFQSFR6RRLE9G/ZvT+jprFIX9oGUmlrHI5JAVHhEEsWayWAcaD3YdhL7gsnQ6rLwJfqeZzuNiFQCFR6RRLHmY2jaC6qfZztJdGg7FA7ugC+X2k4iIpVAp6WLxLGjp4lfwC7mV8nhyZIR/KMcJ4zHtVaXgwmER76adLedRkQiTCM8IglgYGAZANO9bpaTRJFqdaFZ73DhEZG4p8IjkgAGOUvI8xqxwb/AdpTo0vZK2LES9myynUREIkyFRyTOpVLIxc5KpnsX2Y4SfdpcGX69dordHCIScSo8InEu08kh2bjMcFV4TlCvFZzXWtNaIglAhUckzl0eWMoeP5UlfhvbUaJT2yGQPxeK9tlOIiIRpMIjEsccPAY42czyuuASsB0nOrW5ErwSyJtpO4mIRJAKj0gcyzAbOM/sZ5bbxXaU6NW0F6TUgvXTbCcRkQhS4RGJY5cGsnF9wxwvw3aU6BUIQstLYd308HljIhKXVHhE4lh/J5tsvxUFpNqOEt1aD4YDX8FXObaTiEiEqPCIxKsDO+nibNB01plodXn49bqpdnOISMToaAmRGJdWxlER1zlZPJ0Ms73OlZwoBtVoCI06w/rp0O9B22lEJAI0wiMSpy4NZLPTr8VKP812lNjQahBsWQiH9thOIiIRoMIjEoccPPo5K/jU64yvH/Mz03ow+C7kzbKdREQiQL8JReJQF7Oe2uag1u+UR5PuUKV2eFpLROKOCo9IHBpw5HL0LK+T7SixwwlAq8tg3TTwPNtpRKSCqfCIxKFLnWyW+G3Yp8vRy6fVIDi4A75aYTuJiFQwFR6ROFOfvaQ7+cx2dXVWuR29PF27LovEHV2WLhJn+geWAzDb0/qdspR1KT9Afouu4Wmtfg9VYiIRiTSN8IjEmQHOcr72a7PKv9B2lNjUejBsXQSFu20nEZEKpMIjEkcCuGQ6K/jU7QwY23FiU6tB4Hs6PV0kzqjwiMSRrmYdtUwhszSddfYaXwRV6+rydJE4o8IjEkf6B1YQ8h3m6XL0s3f08vT103V5ukgcUeERiSOZzgqW+a3YR3XbUWJby4FwcCfsWGk7iYhUEBUekThRm/1kmI1kuRm2o8S+lgPDr9fPsJtDRCqMCo9InOjjrMQxPlleuu0osa/G+dCgoxYui8QRFR6ROJHprKDAr8YKv4XtKPGh5aWweQEUH7SdREQqgAqPSFzwyQzkMM/rhEvAdpj40OoycIth03zbSUSkAqjwiMSBluZLGptvNJ1VkZpdDMEqmtYSiRMqPCJxINPJAVDhqUhJVeHCS7RwWSROqPCIxIF+zgo2eOez1W9gO0p8aXkZ7FoDBVttJxGRc6TCIxLjkimht7OaOZ4uR69wRy9Pz5tlN4eInDMVHpEY181ZSzVzWNNZkdCgPdRoBHma1hKJdSo8IjEu08mhxA/wmdfBdpT4Y0x4lGfDbPBc22lE5Byo8IjEuExnBUv91hykqu0o8anlQDi0B77Mtp1ERM6BCo9ILDu4i3Qnnzk6TiJyWgwAjC5PF4lxQdsBROQcbJgN6HL0iKpeDxp1hryZpE0ue9owf9ywSgwlIuWlER6RWJY3kz1+Krl+c9tJ4lvLgbB1IakU2k4iImdJhUckVvk+5M1kntcJTz/KkdVyIHghLnZW2U4iImdJvyVFYtXOL2D/duZoOivymvaCpOqlO1qLSOxR4RGJVUeOPMjSguXICyZD80z6OStsJxGRs6TCIxKr8mZCvbZs5zzbSRJDy4GkOV/TzHxtO4mInAVdpSUSi0qKYNM86HY76JinCpU2dtJJP97COMxMCW/0+JrbsJJTici50giPSCzavABCRf8960kiboPfiK1+PU1ricQoFR6RWJQ3E5wkSOtjO0kCMcxx07nEWUmQkO0wIlJOmtISiQHfnmb5OPl99vqt+e4js+0ESlBZXgbfDc6is8ljid/WdhwRKQeN8IjEmPrspYOziSxPV2dVtvleR1zfkBnQ5ekisUaFRyTG9HFyAbT/jgUFpJLjt9B+PCIxSIVHJMZkBlawy6/JKv9C21ES0hwvnS5mPTU5aDuKiJSDCo9ITPHp5+Qw1+uErx9fK7LcDALG52Jnpe0oIlIO+o0pEkPamS3UNwXM1XSWNcv8Vhzwq2haSyTGqPCIxJDMI3vAZLkqPLaECLLA66DCIxJjVHhEYkimk8NarzFfU9d2lISW5aVzobNDx0yIxBAVHpEYkUIxPZ0vNJ0VBY5uCaBRHpHYocIjEiO6OWupYkrIUuGxbqN/Plv9eio8IjFEOy2LRImyDq08qp+TQ7Ef4HOvfSUlkrIZstx0hgU+I4CLS8B2IBE5DY3wiMSIvk4OS/02FFLFdhQhvI6npjlEZ5NnO4qInAEVHpEYUJd9dHLydXVWFJnndcLzjaa1RGKECo9IDDh6nMRcr5PlJHJUAams8JvrXC2RGKHCIxIDMp0c9vrVyfFb2I4ix8jyMuhi1lODQttRROQ0VHhEop5P30AO87yOePqRjSpZbjpB4+mYCZEYoN+eIlGupfmSC8xu7b8ThZb5rTnop2gdj0gMUOERiXJ9j6zf0f470adEx0yIxAwVHpEol+msIN9ryFa/ge0ochJZXgZpztewe6PtKCJyCio8IlEsSIjezmqN7kSx0udmwyy7QUTklFR4RKJYV7OeVFOk9TtRbIPfiG3+eZA303YUETkFFR6RKNY3kIPrGxZ4HWxHkTKFj5lgwxxwQ7bDiEgZVHhEolg/J4flfkv2Ud12FDmFLC8DDhfAl8tsRxGRMqjwiESpmhwgw+SF/zOVqDbP6wgYTWuJRDEVHpEodbGzioDxyXJ1nES020sNuKCLCo9IFFPhEYlSmU4O+/2qZPutbEeRM9FyIGxdBEUFtpOIyEmo8IhEqUwnh8+8DoQI2o4iZ6LlQPBdyJ9rO4mInIQKj0gUamq+5kJnB1k6HT12NOkJSdU1rSUSpVR4RKJQ5pHjJLT/TgwJJkNaXxUekSilwiMShfo6OWzzz2OD38h2FCmPlgNh9wbYk287iYh8ixYHiEQZB48+Ti6T3Z6AsR1HzlDa2Em0NA4zUuAXf3qWN9zLSm/LHzfMYjIRAY3wiESdDLOBWqZQ01kxKM+/gC/9ujo9XSQKqfCIRJm+Tg6eb45sZiexxZDlZtDHycXBsx1GRI6hwiMSZTIDOaz0L2QPNW1HkbOQ5aVTyxSSYTbYjiIix1DhEYki1TnERWadjpOIYfO8jni+IdNZYTuKiBxDhUckivRyVpNkXLK0fidm7aEmuX4amQGt4xGJJio8IlEk08nhkJ/MEq+N7ShyDrK8dLqa9aRSaDuKiByhwiMSRfo6uSz02lFMku0ocg6yvAySjEtvZ7XtKCJyhAqPSJQ4n29o7WxjjqazYt5SrzWFforW8YhEERUekSjRLxD+z1H778S+YpL4zGuv/XhEoogKj0iU6OesYLtflzV+U9tRpAJkeem0cL6iidlpO4qIoMIjEh3cEH2dXLLcdHScRHw4eqVdX43yiEQFFR6RaPDlUmqbg3zqdbadRCrIer8x2/26WscjEiVUeESiwfoZuL5hrtfJdhKpMIYsN50+zkrwXNthRBKeCo9INFg/neV+SwpItZ1EKlCWl05tcxC+zLYdRSThqfCI2Fa4G75cyhwdJxF35h0dscubaTeIiKjwiFi3YTb4HnNcFZ54s5ua5HhpKjwiUUCFR8S29TOgSm2W+y1tJ5EIyPIyYOtCOLzfdhSRhKbCI2KT70PeDGgxAJeA7TQSAVleOnghyJ9rO4pIQlPhEbFpx2rYvx1aXWY7iUTIEq8NBKtqWkvEMhUeEZvWTw+/bqnCE6+KSYK0vpA3y3YUkYSmwiNiU94MqN8eajW2nUQiqeVA+GYd7N1sO4lIwlLhEbGl+CBsmq/prETQ8tLwa43yiFijwiNiS/48cItVeBJB/XZQo1F4RE9ErFDhEbElb0Z4MWuzS2wnkUgzBlpdDnmzwS2xnUYkIanwiNiyfnp4MWtSFdtJpDK0HgSHC2DLQttJRBKSCo+IDXvy4Zv1ms5KJC0GgBOE9dNsJxFJSCo8IjasP7KWQ5ejJ44qtaDZxbBOhUfEBhUeERvyZkKtZlCvte0kUplaXQ5f58K+L20nEUk4KjwilS10OHxgaKvLwotZJXG0Hhx+rVEekUqnwiNS2TbNh+ID0OYK20mksjVoDzWbwLqptpOIJBwVHpHKtm4aBFKgeT/bSaSyGQOtL4cNn0Ko2HYakYSiwiNS2dZ9Er4cPbm67SRiQ+vBULwftnxmO4lIQgnaDiCSSPr/4kU+TVnPo1/14aWxk2zHERua9wcnKTytpVE+kUqjER6RSnSpkw3ALK+L5SRiTUoqXHgJrJtuO4lIQlHhEalEA51l5HmN2Ow3tB1FbGo9GHauhr1bbCcRSRgqPCKVpfggvZzVzPS62k4ith29PF27LotUGhUekcqy4VNSTEjTWRLecLJ2M+3HI1KJVHhEKsu6qRzwq7DIa2c7idhmTHiUZ8On4Y0oRSTiVHhEKoPvw7qpzPXSKdHFkQLhwlNyMLwRpYhEnAqPSGX4eiXs28ZMTWfJUWmZ4Q0oteuySKXQn5oiFSztJPvr3B34gDcyf7EAABiWSURBVJ8nwWxXhUeOSK4GzTNh7Scw5Pe204jEPY3wiFSCSwPLyPXS2EEd21EkmrQZArvzYNc620lE4p4Kj0iE1eQA3cxaTWfJidpeGX695mO7OUQSgAqPSIT1d1YQML6ms+REtZrA+emwZortJCJxT4VHJMIGBZaw069Jtt/KdhSJRm2Hhg8SPfiN7SQicU2LlkUiKIkQA5zlTHZ74unvCzmJq6fW4MMUj5/+7g+84x1/mGj+uGGWUonEH/0GFomgHs4X1DSFTPcush1FolSun8ZXfh0uCyy1HUUkrqnwiETQIGcJRX4SWV667SgSpXwcZrgX0d9ZQTIltuOIxC0VHpGI8RkUWEKWl04RKbbDSBSb7l1Eqimit7PKdhSRuKXCIxIh7cwWmphdTPe62Y4iUW6+15FCP4XLHE1riUSKCo9IhFzuLMHzDTNcrd+RUztMMlleOpcHlgK+7TgicUmFRyRCBgWWkO23ZBe1bEeRGDDdu4jG5hvam822o4jEJRUekQhoyG46OxuY7mo6S87MLLcrnm+43FliO4pIXFLhEYmAywLLAJiq9TtyhnZRi2y/5ZFpLRGpaCo8IhEwyFlMvteQ9X5j21Ekhkx3u9HZ2UAD9tiOIhJ3tNOySAWrRhGXOCt52R0MGNtxJIZM87rxM95kcGAxr7qDSBs7qcz7ahdmkfLRCI9IBevnrCDFhLR+R8ptnd+YPK8Rg53FtqOIxB2N8IhUsCGBhez2U1nst7EdRaLEqUZqjmf4xOvBHYFJ1OIABaRGNJdIItEIj0hFCh1moLOMqW53XAK200gMmuz2JMm4XK5NCEUqlAqPSEXaMJua5hBTvB62k0iMyvGbs80/jyGBRbajiMQVFR6RirTqA/b5VZnvdbKdRGKW4RO3B/2cFVSjyHYYkbihwiNSUdwSWDOJGd5FFJNkO43EsCluD1JMCZc62bajiMQNFR6RipI/Fw7tYYrb03YSiXGL/bbs9GsyJLDQdhSRuKHCI1JRVn8ASdX41MuwnURinIfDNLc7lzrZpFBsO45IXFDhEakIngurP4LWgygixXYaiQNTvB6kmiL6Ojm2o4jEBRUekYqwZSEc3AHth9tOInFigdeRAr8aQxxdrSVSEVR4RCrC6g8gkAytB9tOInGihCDTvYu4PLCUICHbcURingqPyLnyfVj9IbQcCFVq2k4jceQTtwd1zAF6OattRxGJeSo8Iufqy6VQsEXTWVLhPvU6c8CvwjDnM9tRRGKeCo/IuVr5HjhBaHul7SQSZw6TzDSvG1cGFmlaS+QcqfCInAvPg5XvhqezqtW1nUbi0Edub+qYA/R1cm1HEYlpKjwi52LrovB0VqcbbCeROJXlZVDgV+PqwALbUURimgqPyLnIfRsCKdB2qO0kEqeKSeITtweDnMXahFDkHARtBxCJWZ4bns5qM1hXZ0lEfehdzHeCn9LPWcE0rzsAaWMnlXn//HHDKiuaSMzQCI/I2cqfG95ssNONtpNInJvvdeQbv4amtUTOgQqPyNnKfRuSU7XZoEScS4Apbk8ud5ZSlSLbcURikgqPyNkIFcOq98Nrd5Kr2U4jCeBD72KqmcMMdLJtRxGJSVrDI3IWbn/kD/y/5L3cvuRCZi0qey2FSEVZ6LVjh1+bqwILmOT1th1HJOZohEfkLFwdWMBevzpzvXTbUSRBeDhMcnsx0MmmJgdtxxGJOSo8IuVVXMhgZzGT3Z6UaJBUKtE7biYppoShgc9tRxGJOSo8IuW15mNSTRHvuX1tJ5EEk+M3Z713AdcF5tqOIhJzVHhEymv5G2z167HQb2s7iSQcwztuX3o5X9DE7LQdRiSmqPCIlMe+7ZA3k3fcvvj68REL3nf7AHCto1EekfLQAgSRMpxsJ9s7Ah/xcJLHu26mhUQisI36fOa15/pAFs+51wLGdiSRmKDCIwntVNvzn8jnhkAWS71WbPQbRSyTyOm84/blD0n/pIvJI9tvZTuOSEzQmLzIGepgNtHO2cI7Gt0Ryya7vSjyk7gukGU7ikjMUOEROUPXB7Io9gN86F5sO4okuP1UY5rXjeGBBSQRsh1HJCao8IicgQAu1wTmMcO7iAJSbccR4R03kzrmAAN01ITIGVHhETkD/ZwV1Df7NJ0lUSPLS2enX4ubAp/ajiISE1R4RM7ATYFP2e2nMtvrYjuKCAAhgrzl9mOgs4z67LEdRyTqqfCInEY9ChjkLOFtt5+OkpCo8qY7gKDxuCkwx3YUkainwiNyGjcE5pBkXCa4l9qOInKcfL8Rn3nt+U5gNgbPdhyRqKbCI3JKPjcHZrHQa0ue39h2GJETTAhdSprzNb2d1bajiEQ1FR6RU+jtrKaF8xUTQhrdkeg02evJPr8aNwdm2Y4iEtVUeERO4ZbATPb51fjY62U7ishJHSaZd90+XOksohYHbMcRiVoqPCJlqMUBrnQW8Y7blyJSbMcRKdOb7qWkmBKuDcyzHUUkaqnwiJTh+kAWKaaECe5A21FETmmVn8YKrzm3BGYCvu04IlFJhUfkpHxGBGaS7bXkC7+Z7TAipzXBHUh7ZwsXmXW2o4hEJW0qInISPcwa2jjb+HnJHbajiJyR99w+jA2+wfeCU0kb26bM++WPG1aJqUSih0Z4RE7ie8GpFPjV+EAHhUqMKKQKE93+DHU+187LIiehwiPyLQ3ZzRBnIW+6l3KIKrbjiJyxV9zLSTIu3w3MtB1FJOqo8Ih8y8jgdBx8XnEvtx1FpFzy/UbMdjszMjiDJEK244hEFRUekWMkU8KIwExmeF3Z4je0HUek3F5yB9PA7OUKZ5HtKCJRRYVH5BhDnc+pb/bxknuF7SgiZ2W215lNXgNuDU61HUUkqqjwiBzjtuAn5HmNmOd1tB1F5Kz4OLzsDqKns4YOJt92HJGoocIjckRXs44uTh4vuYPx9aMhMWyi259DfjK3B6bYjiISNfRbXeSIO4Mfsdevzltuf9tRRM7JPlKZ6PbnmsA8GugSdRFAhUcEgDSznSucxbzqXk6hLkWXOPCCO5QAHrcFP7EdRSQqqPCIAGMCH1NCgJdCWqws8WGz35ApXg9GBaZTnUO244hYp8IjCe88CrgxMId33Ex2Utt2HJEK84/QVdQ0hdwSmGU7ioh1KjyS8G4NTqOKKeEFd6jtKCIVarnfis+9dowOTiaojQglwanwSEKrwmH+JzCVae5F5PmNbccRqXD/CA2jsfmGYc5ntqOIWKXCIwltRGAmdc0B/i90le0oIhEx0+vKWq8xdwc/xODZjiNijQqPJK6SIn4Q/JAFbgcW++1spxGJCB+H50LX0s7ZwhXOYttxRKxR4ZHEtfRlGpq9POteZzuJSER95F1MnteI+4PvgKdRHklMKjySmEKHYe7TLPTassDrYDuNSER5R0Z52jubYe1k23FErFDhkcS07FXY/yV/DV0PGNtpRCLuA+8S8r2G8OlT4Pu244hUOhUeSTyhYpj7NDTpyTyvk+00IpXCJcDf3Gtg+3JYp5PUJfEY/9RNX38GSExIGzupzNvyxw07/gOLXoBJD8DIt0l78XCEk4lEjyAh1jd6BKrUgjtmg6O/eSXulDlkr8IjceFUhedYVSni05SfstE/n5uLf42msyTR5I84AO/eCTf+CzrdYDuOSEUr85e66r0klNGBKTQwe3mq5BZUdiQhpd8EDTvBzCfALbGdRqTSqPBIwqjNfu4KfshUtxtL/Ta244jY4Thw2SOwewMsfcl2GpFKo8IjCeOe4PtUp4g/hG62HUXErtaDodkl8OkfoPig7TQilUKFRxLCBezi1sA03nb7sd5vYjuOiF3GwKDH4cDXsOB522lEKoUKjySEXyS9jg88HbrRdhSR6NC0J7S7KrxFw77tttOIRJwKj8S9nmY1Vwc+Y7x7Nds5z3Yckegx+AnwQjD9MdtJRCJOhUfimoPHY0kvs80/j/Ghq23HEYkudZvDJffCigmwZaHtNCIRpcIjcW1EYCYdnE08WTKSIlJsxxGJPn1/CjUaweSf6WBRiWsqPBK3anKAB4L/4TOvPZO8XrbjiESnlFQY9Bv4chlkv2Y7jUjEBG0HEImUscE3qMVBHi+5FW0yKHIK6Tex6K0/0uL9X3D5f4LsoeZxN59wPItIDNIIj8SlXmY13w3O4gV3KKv9C23HEYluxvDLkjHUoJBfJb1qO41IRKjwSNxJoZgnk15gs1dfl6GLnKF1fhP+7g7nhsBc+jnLbccRqXAqPBJ37g2+R0tnO78MjdFCZZFyeD50DXleI34X/BdVKbIdR6RCqfBIXGlnNvODwIe87WYy10u3HUckphwmmbEld9DU2ckDwYm244hUKC1alriRQjFPJ/2NvVTntyWjbMcRiUppYyed8vZFfjteDV3G6MAUpnvd+MzrUEnJRCJLIzwSNx4ITqS9s4WfldzFXmrYjiMSs34XGkm+35A/J/2dmuhwUYkPGuGRuHCxs5IxgY95JXQ5s7yutuOIxLRDVOEnJT/k7eTHeDzp36SNrV7mfXXJusQKjfBI7Du0hz8n/Z2N/vn8LjTSdhqRuLDcb8VfQ9dzXWAew535tuOInDMVHoltvg/v30t9CvhxyT26KkukAj3vXsMSrzVPJL1ImtGJ6hLbVHgkts37K3zxEb8PfZccv4XtNCJxxSXAfcX34hLg70nPUIXDtiOJnDUVHoldG7NgxuPQ4Vr+5Q6xnUYkLm2jPveX3ENbs5XfJb0I+LYjiZwVFR6JTfu+hLduh/NawTXPobOyRCJnjteZZ0I3cENgLiMDM2zHETkrKjwSe4oPwhsjoLgQvvMKpOgSdJFI+1/3Wma5nXks+BIXOyttxxEpNxUeiS2eC2/fAV+tgBv/BQ3a2U4kkhB8HO4r+REb/fMZn/Q0Lc0225FEykWFR2LL1F/DmkkwZBy01bodkcq0n2qMLvkZxSTx/5L+wHkU2I4kcsZUeCR2fPZ3+Oxv0Otu6HWX7TQiCWmrX58xxQ9Q3xTwYvKf4PAB25FEzogKj8SGpS/DlLHQ7iq44ne204gktOV+K+4ruZdOZiO8cQuUHLIdSeS0VHgk+q2YCB/cB60uD6/bcQK2E4kkvGled35acjfkz4U3R0FIe/RIdFPhkei28j149y5I6ws3vwpB7aQsEi0+8PrA1X+F9dPhrdEQKrYdSaRMOjxUoteyV+GDH0GTHjDiDdJ+PdN2IhH5tm7fC4/uTH4oPL118yuQXPZhoyK2GN8/5a6Z2lJTKlXa2EkAfD/wMb9OepU5bjp3lfyEQ1SxnExETuWmwGzGBf/JMr81o4sfYh/VdZK62FDmLrSa0pKoYvD4WXACv056lUluT8aUPKiyIxIDJroDuKfkfjJMHm8m/5ZGfGM7kshxNMIj0ePwAaY+MZzBgSW8HhrIr0Kj8dTJRWJKXyeH55Oe4TDJ1B/zFjTtaTuSJJYyR3hUeKTSHZ22OlZjdvJC8p9pY7bw29D/8G/3CnQ+lkhsamW28kLSn0lL2gNXPwtdRtiOJIlDU1oSvQY7i5iU8ksam13cVvJz/u0OQWVHJHat95twbfFvoGkveO8H8OH94bPvRCzSCI9UuqMjPCkUMzb4BrcHPyHHS+PekvvY5J9vOZ2IVJQgIR4ITuTu4Ies8Zpwb8l9rPObAGhBs0SKRngkumSYPN5P/jW3Bz/hhdCV3FD8uMqOSJwJEeSp0Aj+p3gsdc0+Pkx+mDsCHxHAtR1NEpBGeKRyFRfyj9+M4fuBj9lBHX5R8n1me11tpxKRCKvPXp5MepFBgSVkey3ocs+r0LCj7VgSf7RoWSzzfVj1Hkx9BAo281roMsaFRrCfaraTiUil8RnmfM7jSf+mXqAQet4J/X8GVescd6+TXdhwlKbC5DRUeMSibUvgk4dh8wJo2ImbN1/H535726lExJLa7OdnwQncHJjNPqrzdOgG3nAvo+QMNv9X4ZHTUOERC7Yuhk+fgnVToXp9GPhr6DqKtF9OsZ1MRKJAe7OJXwdf4ZLAKrb69fh7aDgT3f4Uk1TmY1R45DRUeKSS+D5smMWcfz9Cv0AOu/1UXggN42V3EAc0fSUiJ/Dp56zg/uA7dHPWsd2vywuhoUx0+7OPE8/kUuGR01DhkQgrKoDsN2DRC/DNOnb5NflnaBivuIMo1NEQInJaPn2cXO4Lvksv5wsK/RTec/vwsjuYL/xmpfdS4ZHTUOGRCAgVw/rpkDMR1kyG0CFo3B163kmbN1JOOSwtIlKWjmYjtwamcU1gHlVMCSu85rzv9uFD92IWjhtlO55ENxUeqSBFBZA3C9Z+wp7sD6hjDrDbT+Uj92Imuv3J8VvYTigicaI2+7k+MJdrA3PJcDbi+QaneV9oMwTaXAHntQKjXdnlOCo8ie6sL/MsKYIvl8Hm+bBhNmyaD14IqtTivYMded/tQ5aXTugMrq4QETlbLcyXXBOYx/0XrIEdq8IfrNsCWg2CtD7Q7BJIrW83pEQDFZ5Ed0aFx3Nh90b4akX4ZfNnsG0puIfDt9dvD20Gh/+6atKTtIc/qYTkIiL/lT9uGOzZFL76c+0nkD83PJ0O5HmNWOS1JddvzirvQlb7zTh0ZA2h1v4kDBWeeHG2IzXHPq4qRVxodnCh+Zpm5muam69o72ymrdlCNRMuNyV+gJV+Ggu9diz22rDIa8sealbcNyIichZO+D0XKobt2bBpPjM+eY+LnHXUMQcAcH3DRr8Ra/wmbPQbsdFrxEb/fDb657OHGhz9v1FlKK6UWXg0DxEnDB4c2gtFe8OvD+6C/dvhwFew/yvGJy2ngdlDE7OLBmbvcY/d46fyhdeMCd6lrPabscpLY53fWIuORST6BZOhaU9o2pPvf9Qa8LmAb+jo5IdfTD7tzWYGO0tICv73DK8DfhW+8uuy3a8L702CmhdAjUaQ2hCq1YWqdcM7QFerCwH9LowH8TvCc/T7Kv3+jn3/bG878v7pbvN98ErCa128UHiqyD32/W+9uMe+fRhKDn3rpRBCRVBSyHuL1lOVYqpRRC1zkFocpJY5SA0KCZgynq4qtVlbWJ0dfm2+9OuR7zdks9+QTX5DNvkN2Edq+f99RURiSJAQTcxOmpuvaGG209jsoqHZTSOzm4tqHwr/geh7J39wcg2oVgeq1ILkVEiqBsnVw28nH/N2UjUIpkAg+chL0jHvJx3z8SMvThCME1547QSOvO2AOeZtJxC+/aQfdzhuQOPoAu7EXsgdpVNaX+XCi4M5q1JR1m3xKJAMSVUhqRobCzyKSOEQyRT41Smg+nGv9x15vcuvxQ5qs9OvzWGSbX8HIiJRLYBLPQqoZwqoYw5QmwPUNgd4YvAFULgbDu0OX6VafDD8UlL437eLD/53rWPUOkkxOvbjZ/qx4z5+mo99W7fvwZDfn1Hac3B2hccYMwWoF4lEcaoesMt2CDmBnpfoo+ckOul5iT56Tspnl+/7Q052w+lGeKQcjDGLfd/vbjuHHE/PS/TRcxKd9LxEHz0nFcexHUBEREQk0lR4REREJO6p8FSsf9gOICel5yX66DmJTnpeoo+ekwqiNTwiIiIS9zTCIyIiInFPhSdCjDEPGGN8Y4wu67fMGPNHY8wXxpgVxph3jTG1bWdKZMaYIcaYNcaY9caYsbbzJDpjTFNjzCxjzCpjzEpjzP22M0mYMSZgjFlmjPnIdpZ4oMITAcaYpsBgYLPtLALANKCT7/sZwFrgF5bzJCxjTAD4G3Al0AEYYYzpYDdVwgsBD/i+3wHoDdyj5yRq3A+sth0iXqjwRMbTwM+I262fY4vv+1N93w8defczoInNPAmuJ7De9/0Nvu8XAxOAayxnSmi+72/3fX/pkbf3E/4PtrHdVGKMaQIMA16wnSVeqPBUMGPMNcA23/eX284iJzUamGw7RAJrDGw55v2t6D/XqGGMSQO6Ap/bTSLAM4T/cC7jgC8pL52WfhaMMdOB809y08PALwlPZ0klOtVz4vv++0fu8zDh4fvXKjObSCwwxqQCbwM/9n1/n+08icwYcxWww/f9JcaYAbbzxAsVnrPg+/7lJ/u4MSYdaA4sN+GD1JoAS40xPX3f/6oSIyacsp6To4wxtwFXAZf52ovBpm1A02Peb3LkY2KRMSaJcNl5zff9d2znEfoAw40xQ4EqQE1jzKu+74+ynCumaR+eCDLG5APdfd/XwW8WGWOGAH8B+vu+v9N2nkRmjAkSXjh+GeGiswj4ru/7K60GS2Am/NfZS8Bu3/d/bDuPHO/ICM+Dvu9fZTtLrNMaHkkEzwE1gGnGmGxjzHjbgRLVkcXj9wKfEF4c+x+VHev6AP8DDDzy85F9ZGRBJK5ohEdERETinkZ4REREJO6p8IiIiEjcU+ERERGRuKfCIyIiInFPhUdERETingqPiIiIxD0VHhEREYl7KjwiIiIS9/4/dHM3Gkot2KAAAAAASUVORK5CYII=\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 far away a proposal state can be from the current state of the chain. It is thus a parameter of the proposal distribution $q$ and controls how big the random steps are which the Markov chain takes. If the step size is too large, the proposal state will often be in the tails of the distribution, where probability is low.\n",
"The Metropolis-Hastings sampler rejects most of these moves, meaning that the acceptance rate decreases and convergence is much slower.\n",
"See for yourself:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Acceptance rate: 0.104\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAFlCAYAAADvZjI4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hUZcL+8e+ZmYQSCF1ABClKSwEMvavYaIKAwoKKrPrq2tsr7+qivOtPsezqa1tWZW2gKIgKRFQsiArSAwkqaCAICErvkMyc8/vjQBQhgYTMPDNn7s915cpMzpR7MleSO895znMsx3EQERER8TKf6QAiIiIi4abCIyIiIp6nwiMiIiKep8IjIiIinqfCIyIiIp6nwiMiIiKeFzjBdh2zLiIiIrHCKmqDRnhERETE81R4RERExPNUeERERMTzVHhERETE81R4RERExPNUeERERMTzVHhERETE81R4RERExPNUeERERMJs586dPP/88xF7voYNG7J161YAOnfuXOxtH3744WK39+7dm507d5KXl0dqamqJcsyZM4d58+YVXh8/fjyvvfZaiR6jrFiOU+xiylppWURE5BTl5eXRt29fcnJyjtkWDAYJBE504oOSadiwIYsXL6ZmzZonvG2lSpXYu3fvMV93HAfHcfD53LGR4l5DUR588EEqVarE3XffffLhT41WWhYRETFl9OjR5Obm0rp1a+655x7mzJlDt27d6N+/Py1btjxm9OSJJ57gwQcfBCA3N5eLL76YjIwMunXrxvfff3/M42/bto0LL7yQlJQUrr32Wn4/mFGpUiUANm3aRPfu3WndujWpqal8+eWXjB49mgMHDtC6dWuGDx9OXl4ezZo146qrriI1NZX169cfNVoUDAYZPnw4LVq0YPDgwezfvx84ekRp8eLF9OzZk7y8PMaPH8+TTz5J69at+fLLL3nwwQd54oknAMjKyqJjx46kp6czcOBAduzYAUDPnj259957ad++PU2bNuXLL78sk/egbCuliIhItJs1GjZnl+1j1kmDS8YVuXncuHHk5OSQlZUFuLt6li5dSk5ODo0aNSIvL6/I+15//fWMHz+es88+mwULFvCXv/yFzz777KjbjB07lq5duzJmzBgyMzOZMGHCMY/zxhtvcNFFF3HfffcRCoXYv38/3bp149lnny3MlZeXxw8//MCrr75Kx44dj3mMVatWMWHCBLp06cKoUaN4/vnnixy9adiwITfccMNRIzyffvpp4farrrqKZ555hh49ejBmzBjGjh3LU089BbjFauHChXzwwQeMHTuWTz75pMjvz8lS4RERETGgffv2NGrUqNjb7N27l3nz5jFkyJDCrx06dOiY282dO5dp06YB0KdPH6pVq3bMbdq1a8eoUaMoKChgwIABtG7d+rjPeeaZZx637ADUr1+fLl26ADBixAiefvrpUu2u2rVrFzt37qRHjx4AXH311Ue9xssuuwyAjIyMYstgSajwiIhIfClmJCaSkpKSCi8HAgFs2y68fvDgQQBs26Zq1aqFIzCnonv37sydO5fMzExGjhzJnXfeyVVXXVVsrj+yLOu413+f/0j2U1GuXDkA/H4/wWDwlB8PNIdHREqo4ejMYj9E5FiVK1dmz549RW6vXbs2v/76K9u2bePQoUPMnDkTgOTkZBo1asSUKVMAdyLx8uXLj7l/9+7deeONNwCYNWtW4XyY31u3bh21a9fmuuuu49prr2Xp0qUAJCQkUFBQcFKv46effmL+/PmAu4usa9eugLv7asmSJQC88847J3zdVapUoVq1aoXzc15//fXC0Z5wUeEREREJsxo1atClSxdSU1O55557jtmekJDAmDFjaN++PRdccAHNmzcv3DZp0iQmTJhAq1atSElJ4f333z/m/g888ABz584lJSWFadOm0aBBg2NuM2fOHFq1akWbNm146623uO222wB3jlB6ejrDhw8/4eto1qwZzz33HC1atGDHjh3ceOONhc9/22230bZtW/x+f+Ht+/Xrx7vvvls4afn3Xn31Ve655x7S09PJyspizJgxJ3z+U6HD0kWkRE40ipM3rk+EkoiIHEOHpYuIiEj8UuERERERz1PhEREREc9T4RERERHPU+ERERERz1PhEREREc9T4REREYkAy7IYMWJE4fVgMEitWrXo27dv4ddmzZpF27ZtadmyJW3atOGuu+4C3LOOW5bFjz/+WHjbp556CsuyWLx4MeCehuK//uu/aNKkCRkZGfTs2ZMFCxZE6NUd38iRI5k6darRDEfo1BIiIhJ3ynpV8JNZfyopKYmcnBwOHDhAhQoVmD17NvXq1SvcnpOTw80330xmZibNmzcnFArxwgsvFG5PS0tj8uTJ3H///QBMmTKFlJSUwu3XXnstjRo14ocffsDn87F27Vq+/fbbMnyVsU0jPCIiIhHSu3dvMjPdsvXmm28ybNiwwm2PPfYY9913X+Eqy36/v3AlY4ABAwYUrrKcm5tLlSpVqFmzZuH1BQsW8NBDD+HzuX/aGzVqRJ8+RxexUCjEyJEjSU1NJS0tjSeffBKAF198kXbt2tGqVSsGDRrE/v37AXeE5sYbb6Rjx440btyYOXPmMGrUKFq0aMHIkSMLH7dSpUrccccdpKSkcP7557Nly5ZjXvuSJUvo0aMHGRkZXHTRRWzatAmAp59+mpYtW5Kens7QoUNL/809ARUeERGRCBk6dCiTJ0/m4MGDrFixgg4dOhRuy8nJISMjo8j7JicnU79+fXJycpg8eTJXXHFF4baVK1fSunXro07rcDxZWVls3LiRnJwcsrOzueaaawD37OSLFi1i+fLltGjRggkTJhTeZ8eOHcyfP58nn3yS/v37c8cdd7By5Uqys7MLT2q6b98+2rZty8qVK+nRowdjx4496nkLCgq45ZZbmDp1KkuWLGHUqFHcd999AIwbN45ly5axYsUKxo8ff5LfyZJT4REREYmQ9PR08vLyePPNN+ndu3eJ73+kML333nsMHDiwxPdv3Lgxa9as4ZZbbuHDDz8kOTkZcMtWt27dSEtLY9KkSaxcubLwPv369cOyLNLS0qhduzZpaWn4fD5SUlLIy8sDwOfzFRawESNG8NVXXx31vKtWrSInJ4cLLriA1q1b89BDD7Fhw4bC78nw4cOZOHEigUD4Ztqo8IiIiERQ//79ufvuu4/anQWQkpJSeMbxovTt25fXX3+dBg0aFJaVI/ddvnw5oVCo2PtXq1aN5cuX07NnT8aPH8+1114LuLuunn32WbKzs3nggQc4ePBg4X3KlSsHuKXmyOUj14PB4HGfx7KOPqWV4zikpKSQlZVFVlYW2dnZfPzxxwBkZmZy0003sXTpUtq1a1fkY54qFR4REZEIGjVqFA888ABpaWlHff2ee+7h4YcfZvXq1QDYtn3MLp6KFSvy6KOPFu4OOqJJkya0bduWBx54gCMnBc/LyyucL3TE1q1bsW2bQYMG8dBDD7F06VIA9uzZQ926dSkoKGDSpEklfk22bRcejfXGG2/QtWvXo7Y3a9aMLVu2MH/+fMDdxbVy5Ups22b9+vWce+65PProo+zatYu9e/eW+PlPho7SEhERiaAzzjiDW2+99Zivp6en89RTTzFs2DD279+PZVlHHbJ+RFETe1966SXuuusuzjrrLCpUqEDNmjV5/PHHj7rNxo0bueaaa7BtG4BHHnkEgL///e906NCBWrVq0aFDB/bs2VOi15SUlMTChQt56KGHOO2003jrrbeO2p6YmMjUqVO59dZb2bVrF8FgkNtvv52mTZsyYsQIdu3aheM43HrrrVStWrVEz32yrCNNsAjFbhSR+HOiw3lP5vBcEfGWSpUqhW1kpoSsojZol5aIiIh4ngqPiIiInJIoGd0plgqPiIiIeJ4Kj4iIiHieCo+IiIh4ngqPiIiIeJ7W4RGRMlXcYes6ZF1ETNEIj4iIiHieCo+IiIh4ngqPiIiIeJ4Kj4iIiHieCo+IiIh4ngqPiIiIeJ4Kj4iIiHieCo+IiIh4ngqPiIiIeJ4Kj4iIiHieCo+IiIh4ngqPiIiIeJ4Kj4iIiHieCo+IiIh4ngqPiIiIeJ4Kj4iIiHieCo+IiIh4ngqPiIiIeJ4Kj4iIiHieCo+IiIh4XsB0ABER0xqOzixyW964PhFMIiLhohEeERER8TwVHhEREfE8FR4RERHxPBUeERER8TwVHhEREfE8FR4RERHxPB2WLiJRQ4eHi0i4aIRHREREPE+FR0RERDxPhUdEREQ8T4VHREREPE+FR0RERDxPhUdEREQ8T4VHREREPE+FR0RERDxPhUdEREQ8T4VHREREPE+FR0RERDxPhUdEREQ8T4VHREREPE9nSxeRiKjCXljzBez9BYKHoGINqNkUajQByzIdT0Q8ToVHRMKmHPkM8X/BAP/XtPWthteOc6NKtSF1EGRcE/F8IhI/VHhEpMxZ2Fzhn8NdgbepZe3mO7sB/ywYzJ2jhkOV+uBPhH1bYfMKyP0MFr4IC8bzREJXHi0Yyhaqmn4JIuIxKjwiUqZqs52nE5+lg+97FtjNuSn/NhY6LQC4s8l5v92w2plwRga0vQb2boF5/0e/r8fTq9wS7i8YxUy7k6FXICJepEnLIlJm2lvfMbPcX0m11nJPwfVckf+3wrJTrEq14MKH6J3/MGucujyb+Az3BSbiJxT+0CISF1R4RKRMXORbxMTEh9ntJHFp/t+ZEuoJlGwycq5Tj8vzx/By8CKuC3zA0wnPkEAwLHlFJL5ol5aInLJ+vnk8mfA8K5zGjMy/l90klfqxggQYG7yaDU4t/pYwkUSe4uaCW8swrYjEIxUeETklPX3LeDLheRY7zfhz/t3so0KZPO6EUG8OkcBDCS/zD8aDfSn4ih6Ubjg6s8hteeP6lEkmEYld2qUlIqWWYq3luYSn+c5pwKj8e8qs7BwxMXQB/6/gT/T1fwOfjCnTxxaR+KLCIyKlUosd/CfxcXZQmVH597Cf8mF5nhdDfXg1eAHMewaWvBqW5xAR71PhEZES8xPimcRnSWY/f86/my1UC+OzWYwNXg2Nz4UP7oFNy8P4XCLiVSo8IlJidwam0NH3HfcVjGKV0yDsz2fjg0EvuaejePtqOLAz7M8pIt6iwiMiJdLR9y03BaYzOdiTaXb3yD1xUk0Y8grsWg+z/jtyzysinqDCIyIn79BeHg/8m7V2bcYGr4r88zfoAN3uhhVvwXczI//8IhKzdFi6SBwq7hBuKOYw7tl/o561lSH5YzgQpknKJ9TtLliVCTNvhwadIKmGmRwiElM0wiMiJ2fdfFj8HyaELmGJ08xcjkAiDPy3O4/n4/vM5RCRmKLCIyInFgpC5l1QpT7/DA42nQZqp0DnW2D5m7Bunuk0IhIDVHhE5MQWvgC/roSLHzG3K+uPut8NVepD5t1uIRMRKYYKj4gUb89m+PxhOOsCaN7XdJrfJCbBxY+4RWzhC6bTiEiUU+ERkeJ9MhZCh+CSR8Eq2dnPw655XzirF3z+MDXYZTqNiEQxFR4RKdrmbHeeTIcboEYT02mOZVlw0SNQsI+bA++ZTiMiUUyFR0SK9slYKF8Fut1pOknRajWFNlcy3P8JDaxfTKcRkSilwiMix7d2Lvw42133pkI4z5VVBnr+DyH83B1423QSEYlSKjwicizHgdljIPkMaH+96TQnllyXl0K96e+fT6q1xnQaEYlCKjwicqxVH8DPy+Dcv0JClByGfgIvBPuy3anEXYEppqOISBRS4RGRP3Dgi0ehemNIv8J0mJO2h4q8FOzDuf7lpFu5puOISJTRubRE5Cjn+rJg03K49Dnwx9aviFdDF3JdIJNbA9O4tuCeMnnMUp93TESiikZ4ROR3HG4LTIOqDWJqdOeIfVRgQvASevmXaS6PiBxFhUdECnXzZdPal+semeVPMB2nVF4NXcQupyK3Bt41HUVEoogKj4gc5o7ubHRqQKs/mQ5TanuoyIRgby70L6GZ9ZPpOCISJVR4RASADtb3tPWtZnywHwQSTcc5Ja+FLmC/U47rAh+YjiIiUUKFR0QAuDaQyVYnmbdDPU1HOWU7qcxboZ70931NbbabjiMiUUCFR0RobP3MBf6lTAz14hCxPbpzxITQJfixuSbwkekoIhIFVHhEhD/7Z3HISeD14AWmo5SZDc5pzLI78Cf/J3Bwt+k4ImKYCo9InKvObgb55zIt1JVtVDEdp0z9O9iXZOsALH3VdBQRMUyFRyTOjfB/QnmrgJdCvU1HKXPZTmPmh1rCN/+CUIHpOCJikAqPSBwrRz5XBj7ms1Brcp16puOExYuh3rB7I3w3w3QUETFIhUckjvXzz6eWtduToztHzLFbQ7WGsPBF01FExCAVHpE4NsI/mx/sesyzU0xHCRsbH7S7Dn6aB5uzTccREUNUeETiVLqVS2vfGl4P9QIs03HCq81wCFSAhS+YTiIihqjwiMSpK/2z2eeUY1qom+ko4VehGqRfDiumwH4tRCgSj1R4ROJQVfbQzz+fd0Nd2UtF03Eio/11EDwAWZNMJxERA1R4ROLQEP8XlLcKeD3knYUGT6hOGjToDIteAjtkOo2IRJgKj0i8sW1G+D9hgd2cVU4D02kiq/11sCMPcj8znUREIixgOoCIRFjup5zp+5XH86+I+FM3HJ0Z8ec8SvO+ULEmLHkFzo6j0S0R0QiPSNxZ9BJbnCp8ZLcznSTyAonQ+k+w+kPYs9l0GhGJIBUekXiy+2f44WPeDvWgIF4HeM+5GuygJi+LxBkVHpF4kvUGODZvhc41ncScmmfBmV1h6Wtg26bTiEiExOm/eCKxr7j5MHnj+hz7RduGZa9Dw2789H3tMCaLARkjYdq1kDcXGvc0HEZEIkEjPCLxYt1X7hFK51xlOol5Lfq5ixEuecV0EhGJEBUekXix9HUoV8X9Yx/vEspD+lD4bibs22o6jYhEgAqPSDw4sAO+fd89vUJCBdNpokPG1WAXuPOaRMTzVHhE4kH2VAgdgnOuNJ0kepzWAup3cOc14ZhOIyJhpsIjEg+WvgZ10qFuK9NJokvr4bB1Na2sXNNJRCTMVHhEvO7nLNi8QpOVjydlAATKM9g/13QSEQkzFR4Rr1v2OgTKQ9pg00miT3l3End//zzKkW86jYiEkQqPiJcFD0H2lN8Ow5Zjtf4TVaz9nO9bajqJiISRCo+Il63+EA7uglbDTCeJXo168LNTXbu1RDxOhUfEy5ZPhkp1tJpwcXx+3g11pYdvObXYYTqNiISJTi0h4lX7tsEPH0PHG8HnN53GqOJOwwHQ2OrOTYHpDPB/zYuhvhFKJSKRpBEeEa/Kecc9K7h2Z53QGud0lthnH96tpTV5RLxIhUfEq1ZMhtppUDvFdJKY8E6oO818G0iz1pqOIiJhoMIj4kVbf4CNS6DVFaaTxIyZoY4cchIY7P/CdBQRCQMVHhEvWj4ZLB+kDTGdJGbsJomP7Qz6+ecTIGg6joiUMRUeEY+xsGHFW9D4XKhcx3ScmPJuqCvVrb10960wHUVEypgKj4jHtLdWwa71mqxcCnPtdLY7lRjo/8p0FBEpYyo8Ih4z0P8lJFaC5n1MR4k5QQJkhjrSy7eUJA6YjiMiZUjr8Ih4SDny6e1fAC0GQmJF03HK1InW0ikr74a6cmXgEy7yLWKa3T0izyki4acRHhEPucC3hGTrALQaajpKzFrqnM1Pdi0G+L82HUVEypAKj4iHXOb/kp+d6tCwm+koMczifbsLXXw51GKn6TAiUkZUeEQ8ojq76e5bwfuhLuDTj/apeC/UBb/l0M8/33QUESkj+q0o4hG9/QsIWLZbeOSU5Dr1yLYbMkBHa4l4hgqPiEf0889ntV2P7536pqN4wnuhLqT71tLY+tl0FBEpAyo8Ih5Ql2108H3P9FBnwDIdxxNmhDpjOxaXavKyiCeo8Ih4QN/Dc01m2J0MJ/GOX6nG13YKA3xfozOoi8Q+FR4RD7jUP48suzHrHJ1Koiy9b3fhTN+vnGP9YDqKiJwiFR6RGNfY+plUXx4zQp1NR/GcD0PtOOgkaLeWiAeo8IjEuP7+ediOxYyQdmeVtb1U5BM7g77+b3QGdZEYp8IjEtMc+vnms8Buwa9UMx3Gk6aHOlHD2kNn30rTUUTkFKjwiMSwFCuPJr5NTNdk5bD5wm7FbqcC/XxahFAklqnwiMSw/v55FDh+ZoXam47iWYdI5GO7HRf5F1GOfNNxRKSUVHhEYpSFTT//fOba6eyksuk4njY91Ilk6wA9fMtNRxGRUlLhEYlRba3VnG5tZ7omK4fdPDuFbU5lnVtLJIap8IjEqP7+eRxwEplttzUdxfOCBPgg1IFevqVU5KDpOCJSCio8IrEoVEBv/wI+tc9hP+VNp4kLM0KdqGDl08u31HQUESkFFR6RWLTmC2pYe7Q7K4IWOc3Y5FSnn3+e6SgiUgoqPCKxKGcqu52KzLFbm04SNxx8zAx1pIdvOcnsNR1HREooYDqAiJc1HJ1Z7Pa8cX1K/qAFB+C7mXwYakc+CaVMJqUxI9SJ6wIfcJF/MVNCPU3HEZES0AiPSKz54WPI38P7ts6dFWkrnMbk2bXp79NuLZFYo8IjEmuyp0LSacy3U0wniUMWM+xOdPatpCa7TIcRkRLQLi2RWHJwN6z+CDJGYs8N3/8rJ9oVF89mhDpxS+A9LvEv4PXQhabjiMhJ0giPSCz5PhNChyBtsOkkcWu1U5/v7fr019FaIjFFhUckluRMhaoN4Ix2ppPEtRmhTrTzreZ0tpqOIiInSYVHJFbs2wq5n0PqILAs02ni2ky7IwB9/N8YTiIiJ0uFRyRWfPseOCFI1e4s09Y5dciyG2u3lkgMUeERiRXZ70Ct5lBbR2dFgxmhTqT58mDrj6ajiMhJUOERiQW7NsJP89zRHe3OigozQ52wHQty3jEdRUROggqPSCxYOc39nHqZ2RxS6Beqs8hp5k4kdxzTcUTkBLQOj0gsyJ4Kp7eBGk1O6uZaRycyZoQ60WHry/DLSqiTajqOiBRDIzwi0W5bLmzK0mTlKPRBqANYfneUR0SimgqPSLTLeQewtDsrCm0nGRr3dN8j7dYSiWoqPCLRzHHc3Vlndobk002nkeNJHQQ7f4KNS0wnEZFiqPCIRLNfVsLWVe4fVYlOLfqCP9EtpiIStVR4RKJZzlR3jkjLAaaTSFHKV4GzL4SV74IdMp1GRIqgwiMSrRzHnRvS5FxIqmE6jRQndRDs3QzrvjadRESKoMIjEq02LHbnhujorOjX9GJIrATZU0wnEZEiqPCIRKucqeAvB837mE4iJ5JYEZr1hm+nQzDfdBoROQ4VHpFoZIfcOSFNL4TyyabTyMlIGwIHd0Lup6aTiMhxqPCIRKO8r2DvL9qdFUuanAsVqmu3lkiUUuERiUY5U905IU0vMp1ETpY/AVpeCqtmQf4+02lE5A9UeESiTTDfnQvSvA8kVDCdRkoibQgU7HdLj4hEFRUekWiT+5k7F0SLDcaeBp0guZ52a4lEIRUekWiTMxUqVIPG55pOIiXl80HKQPjxU9i/3XQaEfkdFR6RaJK/H77/wJ0LEkg0nUZKI20I2AXw3XTTSUTkd1R4RKLJ6g+hYJ92Z8Wyuq2gxlk6t5ZIlFHhEYkmOe9ApTpwZhfTSaS0LMtdTiDvK9j9s+k0InKYCo9ItDi4C36YDamXgc9vOo2cirTBgOMuHikiUSFgOoCIHPbdTAgdKtyd1XB0puFAUmo1z3Z3bWVPgU43mU4jImiERyR65LwDVc+Eehmmk0hZSB0MPy+Dbbmmk4gIKjwi0WHfVlgzxx3dsSzTaaQspA4CLLfIiohxKjwi0SBnGjihw3M/xBOq1IMzO7u7tRzHdBqRuKfCIxINst+G2mlQO8V0EilLqYNg62rYnG06iUjcU+ERMW1bLmxYBOlDTCeRstZyAPgC7urZImKUCo+IadlTgMNrt4i3JNWAJudB9jtg26bTiMQ1FR4RkxwHVrwFjbq5cz7Ee1IHw+4NsH6B6SQicU2FR8SkjUtg+xpIv8J0EgmX5r0hUF5nUBcxTIVHxKQVb7t/DFv0M51EwqVcZWjeB1ZOg2C+6TQicUuFR8SQAEF3jZamF0P5KqbjSDilD4UDO+CHj00nEYlbKjwihnT1ZcP+rdqdFQ+anAdJtWDFZNNJROKWCo+IIQP9X0OFanBWL9NRJNz8AUgbAqs/gv3bTacRiUsqPCIGJHGAC32LIeUyCCSajiORkH4FhPJ1BnURQ1R4RAy40LeYCla+dmfFk7qtoFZzdxkCEYk4FR4RAwb6v+InuxbUb286ikSKZUGroe56PDqDukjEqfCIRFgtdtDFl8N7dhedGT3epF0OWO5yBCISUSo8IhHW3z8fv+XwfqiL6SgSaVXquatqr5isM6iLRJgKj0iEDfB/xQq7EbmOTiURl1oNgx15OtWESISp8IhEUHPrJ9J8ebwT6m46ipjSoh8EKsByrckjEkkqPCIRNMg/l3zHz/uhzqajiCnlKkOLvodPNXHIdBqRuKHCIxIhAYIM8H/FJ3YGO6lsOo6Y1GooHNwFqz80nUQkbqjwiERIT99yalm7mardWdKoJ1SqDcu1Jo9IpKjwiETIYP9ctjhVmGunm44iph051cQPH8O+babTiMSFgOkAIvGgOrs537eUl0MXE/zdj13D0ZkGU4lRrYbB/Gchewp0vMF0GhHP0wiPSAT0988jwQrxTqib6SgSLeqkQt3WsOx1rckjEgEqPCIRMNg/lxV2I1Y5DUxHkWhyzpXwSw5syjKdRMTzVHhEwqyFtY5UXx5TQj1MR5FokzoYAuVh2UTTSUQ8T4VHJMwG++dyyAkwI9TJdBSJNhWqQov+sGIKFBwwnUbE01R4RMIooXDtnXO09o4cX5sRcGgXfDfTdBIRT9NRWiJhdK5vGTWsPUzV7izPKu5Iu7xxfU7iAbpB1TNh2WuQPqQMk4nI72mERySMhvo/Z7NTTWvvSNF8PneUZ+1c96SiIhIWKjwi4bJrAz18y5kS6kEIv+k0Es1aDQMsWDbJdBIRz9IuLZFwWTYRC3gr1NN0Eol2VetDk/Mg6w3oORp8Jy7Ip7wrTSTOaByT6ZoAABO8SURBVIRHJBzsECx9na/sVDY4p5lOI7GgzQjYvQHWzDGdRMSTVHhEwiH3M9i9gTdD55lOIrGieR+oUA2WvmY6iYgnqfCIhMOSV6BiTT6xM0wnkVgRKOfO5fk+E/b+ajqNiOeo8IiUtT2bYdUsaDOcAk2Tk5LIuAbsAq28LBIGKjwiZS1rEjghOOdq00kk1tRq6q7Ls+RlsG3TaUQ8Rf9+ipQl23bnYDTsBjWaAN+bTiQGFXckFRRxNFXba2DqKHce2Nm9wpRMJP5ohEekLK39wl08TqM7UlrN+0FSLVj8H9NJRDxFhUekLC2eABWqQ4t+ppNIrAokQpsrYfUs2LXBdBoRz1DhESkruza4R9iccxUklDedRmJZxtXgODpEXaQMqfCIlJXFL7uf244ym0NiX7WGcFYvWPIqhApMpxHxBBUekbIQPOSuvdP0Yqh2puk04gVtR8HezbD6Q9NJRDxBhUekLKx8D/ZvhfbXmU4iXnH2hZBcDxZNMJ1ExBNUeETKwqIXocZZ0Kin6STiFf6Ae4j6ms9hyyrTaURintbhETmBE66lcuvpsGERXDwOfPofQspQxjXwxeOwYDz0fdJ0GpGYpt/OIqdq4UuQkOSeB0mkLCXVhLQhsHwyHNhhOo1ITFPhETkFVdkDOVMh/XKoUNV0HPGijjdAwX4doi5yilR4RE7BUP/nEDyoycoSPnXS3FOVLHwRQkHTaURilgqPSCklEOTqwMfuH6PaKabjiJd1uAF2rYdVxc8nE5GiqfCIlFJv3zfUtbZD51tMRxGva3YJVD0TvhlvOolIzFLhESkVh+sCH/CjfTqcdYHpMOJ1Pj+0vx5+mgeblptOIxKTVHhESqGj7ztSfXlMCF2iQ9ElMtqMcI8G1CiPSKnoN7VIKVzrz2Srk8y0UDfTUSReVKgKbYZD9hTYvcl0GpGYo8IjUkKNrZ/p5V/GxFAvDpFoOo7Ek45/AScE3zxvOolIzFHhESmhP/tncchJ4PWg5u5IhFVvBCkDYfHLJLPPdBqRmKLCI1ICNdjFIP9cpoW6so0qpuNIPOpyG+TvYbj/U9NJRGKKCo9ICYwKzCKRIC+G+piOIvGqbitoch6jArMoR77pNCIxQ4VH5CQls48r/bP5wG7PGud003EknnW5nVrWLi7zf2k6iUjMUOEROUlX+meTbB3g+eClpqNIvGvUneV2Y673z8SHbTqNSExQ4RE5CRU4yJ8DH/BZqDXfOg1Nx5F4Z1mMD/ajke8XLvYtNJ1GJCYETAcQiQYNRxd/jqJh/s+pbu3l2eCACCUSKd5Hdjty7brcEniPWfntcfT/q0ix9BMicgKJFHB9YCbf2C1Y6jQ1HUcEABsfzwQH0sL3Exf6FpuOIxL1VHhETmCQfy51rB0a3ZGoM8PuRK5dl9sD07A0l0ekWCo8IsVIpICbAu+TZTfhKzvVdByRo4Tw/26UZ4npOCJRTYVHpBhX+D/nDGsr/wgOASzTcUSOMcPuxBq7DrdplEekWCo8IkUozyFuCbzHArs5X9pppuOIHNeRUZ6WvnUa5REpho7SEinCCP8nnGbt5Ob8W9DojoTDiY4OPFnT7c7cYr/LbYFpfJyfUSaPKeI1GuEROY4kDnBjYDpzQ2ksdFqYjiNSrBB+ng5eRkvfOvr6vjEdRyQqqfCIHMdI/0fUsPbwz+AQ01FETsp0uzPf2Q24O/A2BHWOLZE/UuER+YNk9nJ9YCazQ+eQ5ZxlOo7ISbHx8WhwKGf6foUlr5iOIxJ1NIdH5A9uCrxPZQ6c9OhOWc3DEDlVc+xWzA+1pNPcx6D1MChX2XQkkaihER6R3znD+pWR/o94J9SN75wzTccRKSGLccGhsG8LzH/OdBiRqKLCI/I79wYmY+PjieDlpqOIlMpy5yxoeSnMewb2/mo6jkjUUOEROayN9QP9/N/wQqgPv1DddByR0jtvDBQcgDmPmE4iEjVUeEQAcPhrwiS2OFV4IdjXdBiRU1PzLGh3rTt5eXO26TQiUUGFRwTo41tAO99q/hkczD4qmI4jcurO/R8oXxVmjQbHMZ1GxDgVHpFDe7k/YSI5dkPeCp1rOo1I2ahQDc67H9Z9BSvfNZ1GxDgVHpG5j1PX2s6YgpHY+pEQL8kYCbXT4OO/Qf5+02lEjNJvd4lvW3+A+c8xJdidpU5T02lEypbPD5c8Crs3wFdPmk4jYpQWHpS4cPzFAR1eSxhHa18C44LDIp5JJCIadoG0IfD1U+7nWir2Ep80wiNxq49vAd392fwzOJhtVDEdRyR8LnoYEirAzNvBtk2nETFChUfiUlX2MDbhFZbbjXk9dIHpOCLhVek0uODvsO5ryJpoOo2IESo8Epf+ljCRKuzj3oLrCeE3HUck/NpcCQ06uxOY924xnUYk4lR4JO708C1nkP9L/hXqx/dOA9NxRCLD54N+T0H+Ppj136bTiEScCo/ElSQO8P8SJvCjfTrPBgeajiMSWbWaQY97YeU0yJlmOo1IRKnwSFy5PzCR09nGfxdcTz4JpuOIRF7XO+D0cyDzTtjzi+k0IhGjwiNx40LfIoYFPmd8qJ/W3JH45Q/AwH+7JxedcatOOyFxQ4VH4kItdvJIwkvk2A15MjjYdBwRs2o1hfMfgNUfwjIdtSXxQYVHvM9xeCzh3yRxkNsKbqJA622KQIcboGE3mHUvbFltOo1I2Ok3v8SM46+W/Ju8cX2Ov2HBeM71L2dMwdXkOvXCkEwkBvl8cNkLML4rTBkJ133qLk4o4lEa4RFvW78QPr6f2aEMXgtdaDqNSHRJPh0GvgC/roQP/8d0GpGwUuER79q31f3PNbkedxXcAFimE4lEn7N7QZfbYcnLkD3VdBqRsFHhEW+yQzDtOrf0XP4au0kynUgkep13P9TvANNvhV9Wmk4jEhYqPOJNn/0dcj+D3o/B6a1NpxGJbv4EuPw1KJ8Mbw6FfdtMJxIpcyo84j1Zb8BXT0LGSDjnatNpRGJD5TpwxSR3McIpV0OowHQikTKlo7TEW9bNd4flG3WH3k+ApXk7En9KfUTjGRnQ7//gvRvgw9H6GRJPUeER79i+Ft4aDtXOdIfn/Tp1hEiJtR7mHrU17xmocoZ7KgoRD1DhEU+oxU54fSA4NvzpbahQzXQkkdjV639h98/wyYNQuS60Gmo6kcgpU+GRmJfMPl5LHAd7t8HV06FGE9ORRGKbzwcD/gX7tsD7N0FSTTirl+lUIqdEk5YlppXnEBMSH6eJtRGGToQz2pqOJOINgXLuJObTWsBbV0LeV6YTiZwSFR6JWRU4yH8SHifD+oHbCm6GJueZjiTiLeWTYcQ0qFIfJg2BdfNMJxIpNRUeiUlJHOCVxMfo4PuOuwpuYJbdwXQkEW+qdBpcPcOdwDxxsHskpEgMUuGRmFOZ/byWOI4MazW3F9zEu3Y305FEvK1ybbf0JNeFiYPgx09NJxIpMU1alphSm+28kvgYZ1kbubngVj6025uOJBJzSrVOT+U6MDLTHeV543IY+G9IGxymhCJlT4VHIqrUC6IBTa31vJL4KMns588FdzPXblXW8USkOJXrwDWZ8OYweOfP7lFcHW80nUrkpGiXlsSG3M+ZmjgWPzZX5I9R2RExpXwVdyJz877uaswzbodgvulUIiekwiPRzXHg6/+DiZex2anGZYfGstJpaDqVSHxLKO+uZt71DljyMrx2KezdYjqVSLFUeCR6HdrjnsRw9hho0Y+B+f/LRmqZTiUiAD4/9HoQBk2An5fCCz3hp28MhxIpmubwSHTasASmXQc71sIF/wudb2Xf0g+KvcuJ5geJyKk79uesAqnW/TyX/zT1JlxC4Py/Qtc73UIkEkU0wiPRJRSEOY/ChAsglO8eCtvlNp2xWSSK5TiN6ZP/MJl2R/jsIXcX144807FEjqLCI1EjxcqDCb1gzsOQOghu+AoadjUdS0ROwl4qclvBTXDpc/BzFjzfCb4ZD3bIdDQRQLu0JApU5CB3BKYyyj+LLRuTebDgVjIXdYRFX5uOJiIlYkGbEdCoB8y8Az68F3KmQu/H4fQ2psNJnNMIjxjjJ8RQ/2d8Xu5Orgt8wJuh8zj/0BPusLiIxK6q9WH4FHdxwu1r4YVz4b2/wO5NppNJHNMIjxjgcL5vKfcGJtPUt5HFdlNuyL+DZc7ZpoOJSFmxLGg1FJpdAl/+A775F6x8Dzrf7C5WWKGa6YQSZyzHcYrbXuxGkeMp6mgpHzYX+xZyY2A6ab48cu26PBYcykd2W0CTkkW8oMjV0revgdkPwHfTIbEydPgv6HQTVKwe2YDidUX+MVHhkTL3x8JTkYMM8H/Nn/0f0MS3iTV2HcaH+jEt1I2gBhlFPKW408MAsDkH5j4O374PCRWh9TBodx2c1jwyAcXrVHgkco4UnhbWOob7P+FS/zwqWwfIthvyfPBSPrLbYWv6mEjcOaoM/fodfP005LwDoUPQqDtkXOPuAkuoYC6kxDoVHomQbbk8/s9x9PXPp4VvPQedBDLtjkwKns9S52y060okfh139GffVlj6GiyaALs3QLlkaNEP0oa4JUgLGErJqPBImNg2bFoGP34K32fCpiwAFttNmRHqxHuhLuyikuGQIhINit3dZYdg7VzIngLfTof8PZBUC86+0P1och6UT45cWIlVKjxSRmwbtq6C9Qsg7yvI/Qz2bwMsqHcOpAyk8/Rkfqam6aQiEmVOOL/niIIDsGqW+0/Uj7Ph4C7wJUD9DtCwC5zZGc5oB4lJ4Q0ssUiFR0rBcWDXevhlpTvRcMNCWL8QDu50tyfVgibnw1m9oMm5kOSWHJ3TSkSO56QLz++Fgu4/WKs/hLVfwOZscGzwBaBuK6iXAXXS6TNlN6udMyg4zoEQpXpeiVVFFh4dIiNwaK973psda93P29fAL9/Cr9/Cod2/3a5Wc2jZH+p3dP/TqtFE57gSkfDyB9xRnYZd3OsHd7v/eP00D9bNh6w3IP8FMstBvuMn16lHrlOXtU5d1tp1WOvUhX3b3MPfi/h9Vdw/aSpL3qHCY1i4ftCajJ5OZfZTxdpHVfZSy9rFadZOTrN2UItdDE8pB3s2w64NsO/Xo+9cvgqc1hLSL4faKVA7FU5rAeUqlzqPiEiZKJ8MZ/dyP8Ddzb59DTf/8xVSfHk0s9bT0lrHxb5FBAK2e5vHH4BABUg+/fBHvd8uV6xOZ9+P7HQqscOpzHYqc4hEc69Pwsb7hcdx3A/3ym9fO9nrp3pfxwE7WORHupVLgBB+bAJW6PDlEAFsWFng7ssOHoCCg4c/H/4IHvztcsF+dx/37z5yy+897rfDdiy2URl2NoBKtd1CU70RVGsE1Rq6l0+wAqp2WYlIpJ34904nZtqdCq8FCFLf2kIjaxP/6VcDdm+E3T+7H+vmwZ6f3d/DwBt/6DcHnER2UIm9TgV48R/uXKHESoc/H7lcyT18PlAO/AngLwf+RPdy4Mjlwx+BRHcOks8Plh8s3+HLvt9dPomvH9lbUzhS9YfrGnEvltnCszkHJlxIycrESd42RkwvV8zGKcf7onX4h6z8b58TK0L5qlC9sfu5fBX+8eUv7CaJ3U5FdpHEFqcqvzpV2UYyQQLk3Vj06JEKjYjEuiABd7eWUxc6Hef3nW3D/q2wfxuXP5lJNWsP1ay9VGcPVa29VGMPSdZBmpavDPn7YP96yN/rXs7fBwX7Iv+iSux3BaioklTa66WRcTVc/Ejp73+Kip20bFnWh6DDbUqgJrDVdAg5ht6X6KP3JDrpfYk+ek9KZqvjOBcfb8OJjtKSErAsa7HjOG1N55Cj6X2JPnpPopPel+ij96TsaH1/ERER8TwVHhEREfE8FZ6y9YLpAHJcel+ij96T6KT3JfroPSkjmsMjIiIinqcRHhEREfE8FZ4wsSzrLsuyHMuydFi/YZZlPW5Z1veWZa2wLOtdy7Kqms4UzyzLutiyrFWWZf1oWdZo03ninWVZ9S3L+tyyrG8ty1ppWdZtpjOJy7Isv2VZyyzLmmk6ixeo8ISBZVn1gQuBn0xnEQBmA6mO46QDq4H/MZwnblmW5QeeAy4BWgLDLMtqaTZV3AsCdzmO0xLoCNyk9yRq3AZ8ZzqEV6jwhMeTwH8Ta8s+e5TjOB87jhM8fPUb4AyTeeJce+BHx3HWOI6TD0wGLjWcKa45jrPJcZylhy/vwf0DW89sKrEs6wygD/CS6SxeocJTxizLuhTY6DjOctNZ5LhGAbNMh4hj9YD1v7u+Af1xjRqWZTUE2gALzCYR4Cncf5xt00G8wvsnDw0Dy7I+AeocZ9N9wF9xd2dJBBX3njiO8/7h29yHO3w/KZLZRGKBZVmVgHeA2x3H2W06TzyzLKsv8KvjOEssy+ppOo9XqPCUguM4vY73dcuy0oBGwHLLPdHaGcBSy7LaO46zOYIR405R78kRlmWNBPoC5ztai8GkjUD9310/4/DXxCDLshJwy84kx3Gmmc4jdAH6W5bVGygPJFuWNdFxnBGGc8U0rcMTRpZl5QFtHcfRid8MsizrYuCfQA/HcbaYzhPPLMsK4E4cPx+36CwC/uQ4zkqjweKY5f539iqw3XGc203nkaMdHuG523GcvqazxDrN4ZF48CxQGZhtWVaWZVnjTQeKV4cnj98MfIQ7OfZtlR3jugBXAucd/vnIOjyyIOIpGuERERERz9MIj4iIiHieCo+IiIh4ngqPiIiIeJ4Kj4iIiHieCo+IiIh4ngqPiIiIeJ4Kj4iIiHieCo+IiIh43v8HGkSq7mmfsEMAAAAASUVORK5CYII=\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 as 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": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Acceptance rate: 0.985\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAFlCAYAAADvZjI4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5hU5cGG8fvMzC4iC4qAFRUkglQxIGJQwGCLINGgUQOiEuJniRqDGBSDYFCxRWMsxBJjFBTBLhawYMOogCBgRV0V1IAISGdn5nx/jGxEBSm7e2bP3L/rmmt29kx5Zofy7Hvec94gDEMkSZLiLBF1AEmSpMpm4ZEkSbFn4ZEkSbFn4ZEkSbFn4ZEkSbFn4ZEkSbGX+pHtHrMuSZKqi2B9GxzhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSSpki1evJibbrqpyl6vUaNGfPnllwD87Gc/2+B9L7vssg1uP+KII1i8eDGlpaW0atVqk3JMmjSJyZMnl98eOXIk//73vzfpOSpKEIYbPJmyZ1qWJGkLlZaW0qNHD2bNmvW9bel0mlTqxxY+2DSNGjViypQp1K9f/0fvW1JSwrJly773/TAMCcOQRCI3NrKh97A+Q4cOpaSkhPPOO2/jw28Zz7QsSVJUBg0axAcffEDbtm0ZOHAgkyZN4sADD6Rnz560aNHie6MnV199NUOHDgXggw8+4PDDD6ddu3YceOCBvPPOO997/oULF3LooYfSsmVL+vfvz7cHM0pKSgD4/PPP6dy5M23btqVVq1a8+OKLDBo0iJUrV9K2bVt69+5NaWkpzZo1o2/fvrRq1YpPP/10ndGidDpN7969ad68OccccwwrVqwA1h1RmjJlCl27dqW0tJSRI0dy7bXX0rZtW1588UWGDh3K1VdfDcD06dPp2LEjbdq04eijj2bRokUAdO3alT/96U906NCBpk2b8uKLL1bIZ1CxlVKSpHz3xCD4YmbFPueOreEXI9a7ecSIEcyaNYvp06cDuV0906ZNY9asWTRu3JjS0tL1PvbUU09l5MiR7Lnnnrz66qucccYZPPvss+vcZ9iwYRxwwAEMGTKE8ePHc/vtt3/veUaPHs1hhx3G4MGDyWQyrFixggMPPJAbbrihPFdpaSnvv/8+d955Jx07dvzec7z77rvcfvvtdOrUiX79+nHTTTetd/SmUaNGnHbaaeuM8DzzzDPl2/v27cvf//53unTpwpAhQxg2bBjXXXcdkCtWr732Go8//jjDhg3j6aefXu/PZ2NZeCRJikCHDh1o3LjxBu+zbNkyJk+ezLHHHlv+vdWrV3/vfi+88AIPPPAAAN27d6du3brfu8++++5Lv379KCsr46ijjqJt27Y/+Jq77777D5YdgF133ZVOnToB0KdPH66//vrN2l21ZMkSFi9eTJcuXQA46aST1nmPv/rVrwBo167dBsvgprDwSJIKywZGYqpSrVq1yr9OpVJks9ny26tWrQIgm82y7bbblo/AbInOnTvzwgsvMH78eE4++WT++Mc/0rdv3w3m+q4gCH7w9rfzr82+JWrUqAFAMpkknU5v8fOBc3gkKRYaDRq/zkX5pXbt2ixdunS923fYYQfmz5/PwoULWb16NY899hgAderUoXHjxowdOxbITSSeMWPG9x7fuXNnRo8eDcATTzxRPh/m2z7++GN22GEHfve739G/f3+mTZsGQFFREWVlZRv1Pj755BNeeeUVILeL7IADDgByu6+mTp0KwP333/+j73ubbbahbt265fNz7rrrrvLRnspi4ZEkqZLVq1ePTp060apVKwYOHPi97UVFRQwZMoQOHTpwyCGHsNdee5VvGzVqFLfffjt77703LVu25OGHH/7e4y+++GJeeOEFWrZsyQMPPMBuu+32vftMmjSJvffem3322YcxY8ZwzjnnALk5Qm3atKF3794/+j6aNWvGjTfeSPPmzVm0aBGnn356+eufc845tG/fnmQyWX7/I488kgcffLB80vK33XnnnQwcOJA2bdowffp0hgwZ8qOvvyU8LF2SYuC7ozqlI7pHlESKlIelS5KkwmXhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSSpCgRBQJ8+fcpvp9NpGjRoQI8ePcq/98QTT9C+fXtatGjBPvvsw4ABA4DcquNBEDBnzpzy+1533XUEQcCUKVOA3DIU//d//0eTJk1o164dXbt25dVXX62id/fDTj75ZMaNGxdphrVcWkKSVHAq+mzUG3Peo1q1ajFr1ixWrlxJzZo1mThxIrvsskv59lmzZvH73/+e8ePHs9dee5HJZLjlllvKt7du3Zp7772Xiy66CICxY8fSsmXL8u39+/encePGvP/++yQSCT766CPeeuutCnyX1ZsjPJIkVZEjjjiC8eNzZeuee+7hhBNOKN925ZVXMnjw4PKzLCeTyfIzGQMcddRR5WdZ/uCDD9hmm22oX79++e1XX32V4cOHk0jk/mtv3Lgx3buvW8QymQwnn3wyrVq1onXr1lx77bUA3Hrrrey7777svffe9OrVixUrVgC5EZrTTz+djh07ssceezBp0iT69etH8+bNOfnkk8uft6SkhHPPPZeWLVvSrVs3FixY8L33PnXqVLp06UK7du047LDD+PzzzwG4/vrradGiBW3atOH444/f/B/uj7DwSJJURY4//njuvfdeVq1axZtvvsl+++1Xvm3WrFm0a9duvY+tU6cOu+66K7NmzeLee+/luOOOK982e/Zs2rZtu86yDj9k+vTpzJs3j1mzZjFz5kxOOeUUILc6+euvv86MGTNo3rw5t99+e/ljFi1axCuvvMK1115Lz549Offcc5k9ezYzZ84sX9R0+fLltG/fntmzZ9OlSxeGDRu2zuuWlZVx1llnMW7cOKZOnUq/fv0YPHgwACNGjOCNN97gzTffZOTIkRv5k9x0Fh5JkqpImzZtKC0t5Z577uGII47Y5MevLUwPPfQQRx999CY/fo899uDDDz/krLPO4sknn6ROnTpArmwdeOCBtG7dmlGjRjF79uzyxxx55JEEQUDr1q3ZYYcdaN26NYlEgpYtW1JaWgpAIpEoL2B9+vThpZdeWud13333XWbNmsUhhxxC27ZtGT58OHPnzi3/mfTu3Zu7776bVKryZtpYeCRJqkI9e/bkvPPOW2d3FkDLli3LVxxfnx49enDXXXex2267lZeVtY+dMWMGmUxmg4+vW7cuM2bMoGvXrowcOZL+/fsDuV1XN9xwAzNnzuTiiy9m1apV5Y+pUaMGkCs1a79eezudTv/g6wTBuktahWFIy5YtmT59OtOnT2fmzJlMmDABgPHjx3PmmWcybdo09t133/U+55ay8EiSVIX69evHxRdfTOvWrdf5/sCBA7nssst47733AMhms9/bxbP11ltzxRVXlO8OWqtJkya0b9+eiy++mLWLgpeWlpbPF1rryy+/JJvN0qtXL4YPH860adMAWLp0KTvttBNlZWWMGjVqk99TNpstPxpr9OjRHHDAAetsb9asGQsWLOCVV14Bcru4Zs+eTTab5dNPP+Wggw7iiiuuYMmSJSxbtmyTX39jeJSWJElVqGHDhpx99tnf+36bNm247rrrOOGEE1ixYgVBEKxzyPpa65vYe9tttzFgwAB+8pOfULNmTerXr89VV121zn3mzZvHKaecQjabBeDyyy8H4C9/+Qv77bcfDRo0YL/99mPp0qWb9J5q1arFa6+9xvDhw9l+++0ZM2bMOtuLi4sZN24cZ599NkuWLCGdTvOHP/yBpk2b0qdPH5YsWUIYhpx99tlsu+22m/TaGytY2wTXY4MbJUn54buHWW/MYdJSRSkpKam0kZlNFKxvg7u0JElS7Fl4JEnSFsmT0Z0NsvBIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYs/BIkqTYS0UdQJK0rkaDxpd/XTqie4RJpPhwhEeSJMWehUeSJMWehUeSJMWehUeSJMWek5YlKQLfnpgMmzc5+bvPIWn9HOGRJEmx5wiPJFUjjupIm8cRHkmSFHsWHkmSFHsWHkmSFHvO4ZGkGHJ5CmldjvBIkqTYs/BIkqTYc5eWJOWB9R1u7mHoUsVwhEeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWehUeSJMWeS0tIKjiuJL5x/DkpThzhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsWfhkSRJsZeKOoAkKT80GjR+o7aVjuheFXGkCuUIjyRJij0LjyRJij13aUlSzH13V5W7pFSIHOGRJEmxZ+GRJEmxZ+GRJEmxZ+GRJEmxZ+GRJEmx51FakgqaJ9STCoMjPJIkKfYsPJIkKfbcpSUp9ja0RtTGPs7dXVL15giPJEmKPQuPJEmKPQuPJEmKPefwSFKBcW6SCpEjPJIkKfYsPJIkKfYsPJIkKfacwyOpaoQhLP8SVi2BbBqKa0HtHSFZFHUySQXAwiOpcmSz8MkrMGcifPQiLHgH1iz7zp0C2K4x7LofNO4Ce3WHrepEEldSvFl4JFWslYvh9Vth2l2w+GNIpGCX9rBPH6jbGLauB4kErF4KX38O/50F70+EGfdAsgY0PxI6nQ077R31O5EUIxYeSRWjbCVM/jtMvgFWL8mN2Pz8Imh2BNQo2fBjwxDmvg4zx8L0e2DWONjzMDh0ODRoWjX5JcWahUfSlnv/aXh8ACwqhb16QJfzN22EJghg1w65y0GDYcrt8NJ1cPP+sP+Z0PVCKNqq0uJLij+P0pK0+dKr4fHzYVQvSBbDSY/C8aO2bHdUzW3hwAFw1jTY+wR4+W9wS1f4/M0Kiy2p8Fh4JG2erz6E2w6G1/4B+50Op70EjTtX3POXNIBf3gC974eVX+Vea8aYint+SQXFXVqSNt2nr8E9x0M2A8ffA3sdUXmvtefBcPpkuO8kePBUmD8bul0MiWSFv9S3l1yQFC+O8EjaNG89DHceCTXqwO+erdyys1at+tD3IWj/29wurvt/C+k1lf+6kmLDER5JG2/GvfDgadCwPZxwb66IVJVkEfT4a+68PRMuyh0VduydTmaWtFEc4ZG0cdaWncYHQt9HqrbsfNvPzoIe18J7T8HoX+eKjyT9CAuPpB83c9z/ys4JY6B462jztO8HR90MH70AY0+BTFm0eSTlPQuPpA374Fl48P9g9075UXbWansCdL8a3nsCHj4zt5SFJK2Hc3gkrd/nb8KYvlC/GZwwOn/Kzlr79oeVi+DZ4VCrARx2adSJqh2PTFOhsPBI+mFL5sGoY2GrbaDPuNx1PjrwPFg2H165Aeo3hXYnRZ1IUh6y8Ej6vrJVMKYPrFkO/SdCnZ2jTrR+QQCHXQ4L58D4P8J2e+TmGknStziHR9K6wjBXHD6bBr/6B2zfPOpEPy6ZgmPuyJWd+06ERR9HnUhSnrHwSFrXa7fC9FHQ5U+wV/eo02y8mtvmzg2UzcDYk3PrfKlKNBo0vvwi5St3aUn6n09fg6cugKa/gC6Dok6z6eo1gV/eCPedyL+G9mFo+uSoE0nKE47wSMpZuRjG/Rbq7JLblZWopv88tOgJHc/k5NQEeiReiTqNpDxRTf9Fk1ShwhAePQeWfgbH/DN/j8jaWIcMY0q2KSOKbmX34Iuo00jKAxYeSTDt3/DWQ3DQ4Nw6WdVdsoiz1pxFhgTXFt1EkkzUiSRFzMIjFbov34cn/gR7dIVOf4g6TYX5nHr8uawfP03M4Yzkw1HHkRQxJy1LhSybya2RVbQVHF2N5+2sxyPZn/HzzDTOST3A89m9eTNsEnWkWNjQ0Vjf3lY6ohod5afYi9e/bpI2zeS/w7wpcMTVUHvHqNNUiiFlp/Bf6nJt0U1shYeqS4XKwiMVqvnvwHOXQvMjoVWvqNNUmq+pxXllp9Ek8TnnpsZFHUdSRNylJRWiTBoeOh2KS6D7X3PLM8TYK9mWjEp3o3/yccZnOka2a8sT80nRcYRHKkST/5ZbOqL7NVCyfdRpqsSI9AnMpy5XFN1CEemo40iqYhYeqdAs/AAmXQHNe0KrX0WdpsosZWsuKjuF5olPOT35SNRxJFUxd2lJhSQM4bFzIVUDjrgq6jQV7sd2GT2TbcfDmZ/x+9SDPJHtwPthwypKJilqjvBIhWTmWPjoeeg2JLZHZf2YYWV9WcrWjCi6lYBs1HEkVRELj1QoVnwFT10Iu7SD9v2iThOZr6jDZWW9aZd4n2OSL0QdR1IVsfBIheLpobnS0+M6SCSjThOp+7MH8lq2GRekRrMtS6OOI6kKWHikQvDJf2DandDxdNipTdRp8kDAn8tOoQ4rOD81JuowkqqAhUeKu0waHvsj1GkIXS+IOk3eeDfcjTsyh3N88jnaBnOijiOpkll4pLib8k+YPxsOvxxqlESdJq9cl+7FfLZleNE/STiBudI1GjS+/CJVNQuPFGfLF8Jzw6Fxl9wSElrHcmryl7ITaZUopU9y4gbv63/WUvVm4ZHi7LnhsHoZ/OKK2C8fsbnGZ/fjxUwrBqTGUpevo44jqZJYeKS4+vxNmHIHdPgdbN886jR5LOCSdF9qsYpzU/dHHUZSJbHwSHEUhvDEn2Dr7aDroKjT5L33w4aMynSjd/JpmgafRh1HUiWw8EhxNOt++GRy7ozKNetGnaZauDZ9DMuoyZ9TdwFh1HEkVTDX0pLiZs1ymDgEdmwD+5wYdZpqYzG1uS7di4uL7qJbZhrPZNut977fnbhcOqJ7ZceTtIUc4ZHi5uXr4et5ucVBC/yMypvqrswhzMnuzODUKIpIRx1HUgWy8EhxsvQLmHw9tDgKdusYdZpqJ02K4ek+7JH4gpOST0UdR1IFsvBIcTLpcsiUwcEXR52k2pqUbctzmb05O/UA9VgSdRxJFcQ5PFJczH8Hpv0bOpwK2+0RdZpqbXi6D08WD2JAaiwXpvtHHafa8iSNyieO8Ehx8fRQKC6BzudHnaTa+yDchbszB3Nc8jn2DOZGHUdSBbDwSHFQ+hK89wQccC7Uqhd1mli4Pn00y9mKP6XuiTqKpApg4ZGqu2wWJlyUWw294+lRp4mNRdThpvQvOTj5Bh0Tb0UdR9IWsvBI1d3sB+CzN+DnF0FRzajTxModmcOZF9bjwtQoAldTl6o1C49UnaVXwzOXwA6toc2vo04TO6sp5pqyY2mT+IgjE/+JOo6kLWDhkaqz12+DxR/DoZd4ksFK8mD2AGZnd2dgagzFlEUdR9JmsvBI1dXKRfD8ldDk57mLKkVIgsvSv2HXxAL6JidEHUfSZrLwSNXVi3+FVUvgkEuiThJ7L2db83ymDWelHmQblkUdR9JmsPBI1dHiT+DVf8DeJ8COraNOUxAuT/+G2qzkzNTDUUeRtBksPFJ19OxwCAL4+eCokxSMd8LdGJfpzEnJp2gYzI86jqRNZOGRqpvPpsObY3Ln3NmmYdRpCspf08eQJcF5qfuijiJpE1l4pOokDGHin6HmdrmzKqtKfUE9bs/8gqOSk2kdfBh1HEmbwMIjVSdznoaPXoAuf4Kttok6TUH6R/pIFoa1uTA1GgijjiNpI1l4pOoim4GJQ3IrobfvF3WagrWUrflb+lfsn3yLgxLTo44jaSNZeKTqYvpomP8WdLsYUsVRpyloozPd+DC7IxekRpMkE3UcSRvBwiNVB2uWw3OXQsN9ocUvo05T8NKkuDJ9PE0T8+iVfCHqOJI2goVHqg7+cxMs/RwO+UvucHRF7snsvkzL/oQBqbG5Qiopr1l4pHy3bAG89DfYqwfsvn/UaVQu4NKy3uwQLIZXboo6jKQfYeGR8t3zV0DZCjh4aNRJ9B1Tw2Y8mdkXXr4uV0wl5S0Lj5TPvpwDU++AdidD/T2jTqMfcGX6OChbmSumkvKWhUfKZ88MhdRW0HVQ1Em0Hh+GO+cK6dQ7cgVVUl6y8Ej56pP/wNuPQqdzoGT7qNNoQ7oOyhXTZ4ZGnUTSelh4pHwUhjDhz1CyI+x/ZtRp9GNKts8V07cfhU9ejTqNpB9g4ZHy0VsPw9zX4KALobhW1Gm0MfY/E0p2gAkX5QqrpLxi4ZHyTXoNPD0Utm8B+/SJOo02VnGtXEGd+1pupEdSXrHwSPlmyu2w6KPcSQYTyajTaFO07QMN9soV1kxZ1GkkfYuFR8onKxfnDm/eoyv8pFvUabSpkik4eBh89QFM/VfUaSR9i4VHyicvXpMrPS4hUX01PQx2PwAmjYBVX0edRtI3LDxSvlj0Mbw6Etr+BnZqE3Uaba4ggEMvgRVfwuTro04j6RupqANI+sYzl0CQhIMGR51EW2qXdtCqF0y+gQ4TGjGfulEnkgqeIzxSPpg3FWaNyx3avM0uUadRRfj5nyGb5tzUuKiTSMIRHil6a08yWKsBHPCHqNPkvUaDxkcdYeNs1xg6/I5fv3Iz/8z8gvfDhlEnkgqaIzxS1N59HD5+GbpeADVqR51Gm6HRoPHll3V0HshyajIodU80wSSVs/BIUcqUwcSLoX5T+OlJUadRRdt6O25M/5JuyTfYPzE76jRSQbPwSFGa+i9Y+D4ccknuHC6KnX9lDmNuWJ/BqVEkyEYdRypYFh4pKisXwXOXQaMDoenhUadRJVlNMVeUHU+rRCnHJJ+POo5UsCw8UlSevzJXeg67zJMMxtyj2f15PduUgakxlLAi6jhSQbLwSFFY8B68dgu0O8mTDBaEgEvK+tIg+Jrfpx6KOoxUkCw8UhSeuhCKtoaDLoo6iarIzHAPxqY70y/5BLsHX0QdRyo4zpKUqtp7E2DORDj0UihpEHWavFdtzruzEa5MH8cvkq8xODWKU8sGRB1HKiiO8EhVKb0mN7pT7yfQ4dSo06iKLaAuN6aP4tDkVDolZkYdRyooFh6pKr1+a+4w9MMug1Rx1GkUgX9mDueTbAOGpO4iSSbqOFLBsPBIVWX5lzDpCmjSDfY8NOo0ishqirk03YdmibmckHw26jhSwbDwSFXl2eGwZhkcfrmHoRe4p7LtmZxpwYDUWLZhWdRxpIJg4ZGqwrypubMqdzgVGjSLOo0iF3BJui+1WcHA1Jiow0gFwaO0pMqWzcD4AVCyPRx0QdRpqoU4HZm1Pu+Eu3Fn5jBOST7JfZmuvBk2iTqSFGuO8EiVbeod8NkbucPQt9om6jTKI9eme7GAbRhe9E/X2ZIqmYVHqkzLFsAzl+TWy2p9TNRplGeWsTWXlvWhTeIjfpN8Juo4UqxZeKTKNHEIrFkB3a9xorJ+0CPZ/Xk505KBqTHUY0nUcaTYcg6PVFk+fgVmjIYDznWisjYgYEj6ZJ4oHsQFRfdwXtlpUQeqEt+dp1U6ontESVQoHOGRKkOmDMb/EbbZFToPjDqN8twH4S7clunOMckX2Dd4J+o4Uiw5wiNVhv/cDPPfguPuhuJaUadRFdrcI8z+nj6KnsnJXFp0O93XXE6Z/zxLFcoRHqmiffUhPHcZNDsC9uoRdRpVEyvZiiFlJ9M0MY8zkg9HHUeKHQuPVJHCEB49B5JFTlTWJns2+1MeyezPmamH2DOYG3UcKVYsPFJFeuNu+OgFOHgo1Nk56jSqhoaV9WUZNbmi6BbPzSNVIAuPVFGW/hcmDIbdfgbtTok6jaqphWzDJWV9+WliDn2TE6KOI8WGhUeqKE8MhLJV0PN6SPhXS5vvoWwnnsvszcDUGBoGC6KOI8WC/ypLFeHtR+Gth6HL+VB/z6jTqNoLuKisHwCXpW4DwmjjSDFg4ZG21PIv4bFzYcfW0OmcqNMoJubRgCvTx9M5OZNjk89HHUeq9iw80pZYe1TWqiVw9C25o7OkCnJX5hBeze7FkNRd7tqStpCFR9oSb46Bdx6DgwbDDi2iTqOYyZJgQNlpBIRcXTSSwKO2pM1m4ZE215K58Pj5sGtH+NlZUadRTM0Nt2dYui8dE2/z2+QTUceRqi0Lj7Q5whAe/j1k03D0zZBIRp1IMTY204WJmXYMTI2hafBp1HGkasnCI22O126FD5+DQ/8C2+0RdRrFXsCgsv4sZWuuLbqJItJRB5KqHQuPtKm+mAkTLoI9D4X2/aJOowKxkG24oKw/LRMfMyB1X9RxpGrHwiNtijXLYewpULMuHHWza2WpSk3MtufudDdOSz1G18QbUceRqhULj7QpHj8fFs6BXrdCrfpRp1EB+kv6RN7O7sZfi25mRxZGHUeqNiw80sZ68z6Yfjd0HgiNO0edRgVqNcWcWXY2NSjjb8U3kiQTdSSpWrDwSBvjy/dzZ1PebX/o8qeo06jAfRjuzIVlv2W/xDv8IXV/1HGkasHCI/2Y1Uvh3t6Q2gp63QbJVNSJJB7OHsCYdFfOTD5M58SMqONIec/CI21IGMJDZ8DC9+HYO2CbhlEnkspdnD6J98KGXF90A7sG/406jpTX/FVV2pDJ18Pbj8Chw523o7yzihqcWlCwSJcAAAyESURBVPZHHim+iFuK/kqvNcNYwVZRx8prjQaNX+d26YjuESVRVXOER1qfDyfB00Oh5dGw/++jTiP9oE/CHTir7CyaBnO5qmgkEEYdabM0GjS+/CJVBkd4pB/y5ftw30lQvxn0vMHz7VQB/6PbfC9m2zAifQKDi0YzO/swN2WOijqSlHcsPNJ3LV8Io46FRAp+cy/UKIk6kfSjbs10p1WilPNSY/kw3Jknsx2ijiTlFXdpSd+WXg1j+sDXn8EJ90DdRlEnkjZSwPllp/JG+BOuK7qRnwbvRR1IyisWHmmtMIRHzoZPJudWQN/V35BVvaymmN+tGcAX4XbcWnwNuwdfRB1JyhsWHmmtZ4bBm/fCQYOhVa+o00ib5SvqcErZ+SQIuaPoSuryddSRpLxg4ZEAXr4eXro2t/p554FRp5G2yEfhTvRfM4BdgoXcUXwVtVgZdSQpchYeadpdMPHP0PJXcMTVHpGlWJgaNuP3ZWfRKviI24quoQZroo4kRcrCo8L21iPw6NnQpBsc/Q9IJKNOJFWYidn2/LHsdPZLvM3IomspIh11JCkyFh4VrrcfhXGnwC7t4bi7IFUcdSKpwj2S7cTgdD8OSs7guqIbXF1dBcvCo8I0+6HciQV3/in0GQfFtaJOJFWaezLdGF7Wm+7J1/h70d8d6VFB8sSDKjyz7of7f5c77Lz3WKhRO+pEUqW7LdOdLAmGFN1FDa7ljLJzWI2jmt8+w7frasWbIzwqLNPugvv7w24dofc4y44Kyj8zv+DCst/SLfkGtxddRU1WRR1JqjKO8KgwhCG8cDU8Nxya/ByOu9vdWCpIozPdWBUWcVXRPxhVfBn915zHV9SJOlalcp02gSM8KgTZDIwfkCs7e58Av7nPsqOC9kC2M2eUnUOL4GMeKL6YRsHnUUeSKp2FR/G2eimMORGm3A4HnAtH3QzJoqhTSZF7KtuB36wZTO1gBQ8UX0y74N2oI0mVysKj+Fr4Adx2CLz3JPziKjh4qCcVlL5lWtiUX60ZxuKwhNHFl9Er8ULUkaRKY+FRPM15Bm49CJZ9ASc+APudGnUiKS99HO7Ir9YMY2p2T64pHslfUv+kmLKoY0kVzknLipdsBl68BiZdDg2aw/GjYLvGUaeS8tpianNi2QUMDO/jtNSjtEqUcvqac/iCepHk8VBxVQZHeBQfS+bCnT3huUtz62L9doJlR9pIGZKMSJ/AaWv+wJ7BXB6vcQGHJV6LOpZUYSw8ioe3HoGbO8Fnb+QmJve6DWqURJ1KqnaezHag55rhzA0b8I/i67gy9Q9XW1csuEtL1duyBfDE+TD7AdipLRzzT6jXJOpU2gieGyV/fRjuTK81wzg79QBnJB9mv8TbDEr/jleyLaOOto4N7fryz5e+y8Kj6ikM4c0x8OQgWLMcDroIOp3jAqBSBSkjxTXpX/N8pg3XFI3knuJLGZfpzKVlv2FRhCcqtMhoc7lLS9XPFzPhziPhwf+D+k3htJegy0DLjlQJpoR7ceiaK7kx3ZNfJl7mmRrncWxyEgmyUUeTNomFR9XHsgXwyNkw8kD472w44mo45Qlo0CzqZFKsraaYq9LH033NZXwQ7sxVRbfwWPFgDkjMjDqatNHcpaX8t+IreOVGePUfkF4JHU+HLudDzbpRJ5MKynvhrvx6zRB6JP7D+akx3F18Oc9n2nBl+jhmhx4RqfwWhGG4oe0b3ChVqpWL4JWb4NWRsPpraPHL3FydBk2jTqYK4FyM6q2YMvomJ3BW6kG2CVbwbKYtN6SPYlpYff9+es6fWFjv6fQtPMo/X86BV2+G6aOhbAU0PxK6DIIdW0WdTBXIwhMPtVnBickJ9E89znbBMiZnWvCvzGE8k/0pGZJRx9skFp5YsPAoz2XK4P2JMO1OeO+p3AKfrY+FjmdYdGLKwhMvNVlF7+QznJJ6kl2ChcwN6zM63Y17MwfxVYRHdW0KC08sWHiUp+a/A9PvhhljYPl8qLU9tD8F2v8Wau8QdTpVIgtPPCXJcHBiGicmJ3BAcjZrwiSTsm15MHMAz2b3YTX5ezSlhScW1lt4nLSsqhWG8N9Z8PajubMjL3gbEiloejjs0wd+cnBudEdStZQhyVPZfXkquy9N0vM4PvkcPZOTOTQ5la/DmjyZ6cCEbHteyrZiFTWijqsC4giPKt+qJfDRi/DhJJjzNCz6CIIE7PYzaNEzt+5VSYOoU6qKOcJTOBJk2T8xm6OTL3NY4nVqBytZFRbxUrYVz2Z/ysvZlnwc7sAGfjmvEo7wxIK7tFSFls2HuVNg3hT46AWYNxXCLBTVgkadYK/u0Ky7JafAWXgKUxFpOiTe5uDENA5OTGPXxAIAPgu34z/ZFvwn25wp2WZ8FO5IGOGp4iw/1ZaFR5UgDOHrz3K7pea/DfOm5UrO4k9y2xOp3PpWTQ6CPQ6Chvt6NmSVs/AIQpoEn7F/4i06Jt5mv8RbNAi+BmBpWJPZYSPezO7BrGxjZoe783G4A+kqmolh4am2LDzaTGEIy7+ExR/DotLcZfHHsODd3ITj1Uv+d986DaFh+9xll/aw095QvHVUyZXnLDz6vpCfBPPYJzGH1sFHtE58RIvgY2oEZQCUhUk+Cbfng3BnPgh35sNwJz7J7sA86vHfcDvKKrAMWXiqLSct61uy2dyJ/FYtyV2Wz88t27Dsv7ndUcvn/+/rxZ9C2fJ1H791/dxyDq2Pge2b5y4NmkOtetG8H0kxETAnbMicTEPG0hWAFGmaBnNpFnxKk8RnNAlyl66J6RQHmfJHZsOA+WzLZ2E9PgvrsSDclq/C2nxFHRaGuctX1GZhWIevqUXWlZUKTmGN8Kx9r+XveT23N+Y+W/oc2cw3l/S3LptwO70qdylb+c31qtyyC9+7Xgmrvik3a0vO6qWs96Mt2hpKts8dHl6yPWy7G2y7O9TdPXe97W5Qo2SDP2ZpYzjCoy2RJMOuwXx2Cb5k52AhO7Mwd/3N7QbBYuoEK9f7+OVhDZZRk2VhzXWvqcmKsAYnHtAMUjUgtdX/rou2Wvd2shgSydzu+/JLEoLkurfXuU7lDtoIAiDYjOtvHhtEO8E7j+XpLq0vZsHth1LpRaNQBEkoqvnNX8xvXW+1Te5So843X9dZ93bJ9lCrAZTsYJlRlbHwqLIVU0ZdllIv+JrtgqVsx9fUC76mDiuoFayihBXUDlZSwkpKvrmuHaykJqupVyPM/TKZTUf9Nn7EBsrRBh+2oe1b8tgNPL7dSXD45T/y2C22eYUnCIIngfqVkSim6gNfRh1C3+Pnkn/8TPKTn0v+8TPZNF+GYXj4D234sREebYIgCKaEYdg+6hxal59L/vEzyU9+LvnHz6TiOGtLkiTFnoVHkiTFnoWnYt0SdQD9ID+X/ONnkp/8XPKPn0kFcQ6PJEmKPUd4JElS7Fl4KkkQBAOCIAiDIPCw/ogFQXBVEATvBEHwZhAEDwZBsG3UmQpZEASHB0HwbhAEc4IgGBR1nkIXBMGuQRA8FwTBW0EQzA6C4JyoMyknCIJkEARvBEHwWNRZ4sDCUwmCINgVOBT4JOosAmAi0CoMwzbAe8AFEecpWEEQJIEbgV8ALYATgiBoEW2qgpcGBoRh2ALoCJzpZ5I3zgHejjpEXFh4Kse1wPkU3Kme81MYhhPCMFx7utT/AA2jzFPgOgBzwjD8MAzDNcC9wC8jzlTQwjD8PAzDad98vZTcf7C7RJtKQRA0BLoDt0WdJS4sPBUsCIJfAvPCMJwRdRb9oH7AE1GHKGC7AJ9+6/Zc/M81bwRB0AjYB3g12iQCriP3i3M26iBx4WrpmyEIgqeBHX9g02DgQnK7s1SFNvSZhGH48Df3GUxu+H5UVWaTqoMgCEqA+4E/hGH4ddR5ClkQBD2A+WEYTg2CoGvUeeLCwrMZwjA8+Ie+HwRBa6AxMCPILa7WEJgWBEGHMAy/qMKIBWd9n8laQRCcDPQAuoWeiyFK84Bdv3W74TffU4SCICgiV3ZGhWH4QNR5RCegZxAERwBbAXWCILg7DMM+Eeeq1jwPTyUKgqAUaB+GoQu/RSgIgsOBvwJdwjBcEHWeQhYEQYrcxPFu5IrO68BvwjCcHWmwAhbkfju7E/gqDMM/RJ1H6/pmhOe8MAx7RJ2lunMOjwrBDUBtYGIQBNODIBgZdaBC9c3k8d8DT5GbHHufZSdynYATgZ9/8/dj+jcjC1KsOMIjSZJizxEeSZIUexYeSZIUexYeSZIUexYeSZIUexYeSZIUexYeSZIUexYeSZIUexYeSZIUe/8P/iirDt8pxxgAAAAASUVORK5CYII=\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 of the Markov chain (and thus samples) we'll 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": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Acceptance rate: 0.990\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAFlCAYAAADvZjI4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3yV5eH+8c99TgbLALJEVtgbwkaGIMgQEHFVLahIta1aV9WWb7Uqra1o7VdrrT/aitZaceFADSgBDFCWbAgzjAAie4cwcp7n+f1xIF8RCARycp/znOv9evHKeE5yrhAlV+7x3MbzPERERET8LGA7gIiIiEikqfCIiIiI76nwiIiIiO+p8IiIiIjvqfCIiIiI76nwiIiIiO8lnOO69qyLiIhIrDBnu6ARHhEREfE9FR4RERHxPRUeERER8T0VHhEREfE9FR4RERHxPRUeERER8T0VHhEREfE9FR4RERHxPRUeERGRCNu/fz+vvfZaiT1famoqu3fvBqBLly6FPvaPf/xjodcHDBjA/v37ycnJoUWLFkXKkZmZyezZswveHjNmDP/+97+L9DmKi/G8Qm+mrDsti4iIXKScnBwGDRpEVlbWaddCoRAJCec6+KBoUlNTWbBgAZUrVz7nY8uVK0dubu5p7/c8D8/zCATCYyOFfQ1n88wzz1CuXDkee+yx8w9/cXSnZREREVtGjhzJ+vXrSUtL4/HHHyczM5Pu3bszePBgmjVrdtroyYsvvsgzzzwDwPr16+nfvz/t2rWje/furF69+rTPv2fPHvr27Uvz5s25++67+f5gRrly5QDYtm0bV155JWlpabRo0YKZM2cycuRIjhw5QlpaGkOHDiUnJ4fGjRtzxx130KJFC7Zs2XLKaFEoFGLo0KE0bdqUm266iby8PODUEaUFCxbQs2dPcnJyGDNmDC+99BJpaWnMnDmTZ555hhdffBGAJUuW0LlzZ1q1asX111/Pvn37AOjZsye//vWv6dixI40aNWLmzJnF8j0o3kopIiIS7SaNhO3Li/dzXtYSrhl91sujR48mKyuLJUuWAOGpnkWLFpGVlUXdunXJyck568f+9Kc/ZcyYMTRs2JB58+Zx3333MW3atFMeM2rUKLp168ZTTz1Feno6Y8eOPe3zjBs3jn79+vHEE0/gOA55eXl0796dV199tSBXTk4O2dnZvPXWW3Tu3Pm0z7FmzRrGjh1L165dGTFiBK+99tpZR29SU1P5+c9/fsoIz9SpUwuu33HHHfz1r3+lR48ePPXUU4waNYqXX34ZCBerb775hokTJzJq1CimTJly1r+f86XCIyIiYkHHjh2pW7duoY/Jzc1l9uzZ3HzzzQXvO3bs2GmPmzFjBh9//DEAAwcOpGLFiqc9pkOHDowYMYL8/HyGDBlCWlraGZ+zTp06Zyw7ALVq1aJr164ADBs2jFdeeeWCpqsOHDjA/v376dGjBwB33nnnKV/jDTfcAEC7du0KLYNFocIjIiLxpZCRmJJUtmzZgtcTEhJwXbfg7aNHjwLgui4VKlQoGIG5GFdeeSUzZswgPT2d4cOH88tf/pI77rij0Fw/ZIw549vfz38y+8VITk4GIBgMEgqFLvrzgdbwiIhFqSPTC/0j4heXXHIJhw4dOuv1atWqsXPnTvbs2cOxY8f44osvAEhJSaFu3bp8+OGHQHgh8dKlS0/7+CuvvJJx48YBMGnSpIL1MN+3adMmqlWrxj333MPdd9/NokWLAEhMTCQ/P/+8vo7NmzczZ84cIDxF1q1bNyA8fbVw4UIAPvroo3N+3eXLl6dixYoF63PefvvtgtGeSFHhERERibBKlSrRtWtXWrRoweOPP37a9cTERJ566ik6duxInz59aNKkScG1d955h7Fjx9K6dWuaN2/OhAkTTvv4p59+mhkzZtC8eXM+/vhjateufdpjMjMzad26NW3atOH999/noYceAsJrhFq1asXQoUPP+XU0btyYv/3tbzRt2pR9+/Zx7733Fjz/Qw89RPv27QkGgwWPv/baa/nkk08KFi1/31tvvcXjjz9Oq1atWLJkCU899dQ5n/9iaFu6iFhzrlGcnNEDSyiJiPiEtqWLiIhI/FLhEREREd9T4RERERHfU+ERERER31PhEREREd9T4RERERHfU+EREREpAcYYhg0bVvB2KBSiSpUqDBo0qOB9kyZNon379jRr1ow2bdrw6KOPAuFTx40xrFu3ruCxL7/8MsYYFixYAISPofjZz35G/fr1adeuHT179mTevHkl9NWd2fDhwxk/frzVDCfpaAkREYk7xX0n7/O5Z1TZsmXJysriyJEjlC5dmoyMDGrUqFFwPSsri1/84hekp6fTpEkTHMfhH//4R8H1li1b8t577/Hkk08C8OGHH9K8efOC63fffTd169YlOzubQCDAxo0bWblyZTF+lbFNIzwiIiIlZMCAAaSnh8vWu+++y2233VZw7YUXXuCJJ54ouMtyMBgsuJMxwJAhQwrusrx+/XrKly9P5cqVC96eN28ezz77LIFA+Ed73bp1GTjw1CLmOA7Dhw+nRYsWtGzZkpdeegmAf/7zn3To0IHWrVtz4403kpeXB4RHaO699146d+5MvXr1yMzMZMSIETRt2pThw4cXfN5y5crxyCOP0Lx5c3r37s2uXbtO+9oXLlxIjx49aNeuHf369WPbtm0AvPLKKzRr1oxWrVpx6623Xvhf7jmo8IiIiJSQW2+9lffee4+jR4+ybNkyOnXqVHAtKyuLdu3anfVjU1JSqFWrFllZWbz33nvccsstBddWrFhBWlraKcc6nMmSJUvYunUrWVlZLF++nLvuugsIn04+f/58li5dStOmTRk7dmzBx+zbt485c+bw0ksvMXjwYB555BFWrFjB8uXLCw41PXz4MO3bt2fFihX06NGDUaNGnfK8+fn5PPDAA4wfP56FCxcyYsQInnjiCQBGjx7N4sWLWbZsGWPGjDnPv8mi05SWiETekf2wZiJsngMHvoX8o3BJNe4LJjLNbctq7/Rzf0BHT4j/tGrVipycHN59910GDBhQ5I8/WZi++uorpk6dyptvvlmkj69Xrx4bNmzggQceYODAgfTt2xcIl60nn3yS/fv3k5ubS79+/Qo+5tprr8UYQ8uWLalWrRotW7YEoHnz5uTk5JCWlkYgECgoYMOGDeOGG2445XnXrFlDVlYWffr0AcIjTdWrVy/4Oxk6dChDhgxhyJAhRf47OV8a4RGRyMndBV/8Ev7cBD69F1Z+BkcPgjGwbSm/SvyAL5NH8mHSM3QNLLedVqREDB48mMcee+yU6SwIF4iTJ46fzaBBg3j77bepXbs2KSkpp3zs0qVLcRyn0I+vWLEiS5cupWfPnowZM4a7774bCE9dvfrqqyxfvpynn36ao0ePFnxMcnIyAIFAoOD1k2+HQqEzPo8xpx5p5XkezZs3Z8mSJSxZsoTly5czefJkANLT07n//vtZtGgRHTp0OOvnvFgqPCISGUvfh7+2g0VvQaub4Z5p8KuNcM9UuGsiPLiYDkdf45n8O6hu9vJO0nO8lPg3UjhsO7lIRI0YMYKnn366YKTkpMcff5w//vGPrF27FgDXdU+b4ilTpgzPP/98wXTQSfXr16d9+/Y8/fTTnDwUPCcnp2C90Em7d+/GdV1uvPFGnn32WRYtWgTAoUOHqF69Ovn5+bzzzjtF/ppc1y3YjTVu3Di6det2yvXGjRuza9cu5syZA4SnuFasWIHrumzZsoWrrrqK559/ngMHDpCbm1vk5z8fmtISkeLlhGDiY7DwTajdBa79C1RpdMaH7qIC/3L6M87pzf0JE7gvOIFWSRv4af4vWe/VOOPHiMS6mjVr8uCDD572/latWvHyyy9z2223kZeXhzHmlC3rJ51tYe/rr7/Oo48+SoMGDShdujSVK1fmT3/60ymP2bp1K3fddReu6wLw3HPPAfD73/+eTp06UaVKFTp16sShQ4eK9DWVLVuWb775hmeffZaqVavy/vvvn3I9KSmJ8ePH8+CDD3LgwAFCoRAPP/wwjRo1YtiwYRw4cADP83jwwQepUKFCkZ77fJmTTfAsCr0oInKK/KMw/q7wep2uD1N/ajscCl9E+X0dzGr+X9LLBHEZdvx/WOHVLfTxWsMjEh3KlSsXsZGZIjJnu6ApLREpHk7oRNmZBANehD6jilR2AOZ7Tbj++CgOU4pxSX+gmcmJTFYRiTsqPCJy8TwvPI21ZiIM+BN0vOeCP9UWrxq3HPstuZTmjaQ/UZ09xRhURCIhSkZ3CqXCIyIXb86r4TU73X55UWXnpK1U4a7jv6IMR3kj6U8kc7wYQopIPFPhEZGLs+UbyHgamg6G3k8V26dd69XigfwHaRrYzNMJbxXb5xWR+KTCIyIXLm8vjB8BFWrBda+G769TjKa7rXktNJgfJ3zN4MDsYv3cIhJfVHhE5MJNfAwObYeb3oRS5SPyFH8O3cx8txF/SBzLZVrPIyIXSIVHRC7Mmi8h6yPo8Suo0TZiT+MQ5Jf595KAw+8T30R3yxCRC6HCIyJFd/QgfPEIVG0GXR+O+NNt8arx59DN9AkuYlBgbsSfT0T8R4VHRIru6z9A7nYY/CokJJXIU77p9GeJW4+nE9/iEvJK5DlFxD9UeESkaHatgW/+CW3vhJrtSuxpHYI8mT+CShzi/oQJJfa8IuIPOktLRIpm8pOQVBZ6PVniT53l1eNjtzt3BSfxjtPrnI9PHZle6HUdTSESPzTCIyLnb90UyJ4cXqhctrKVCC/k34JDkN8kjLPy/CISm1R4ROT8eB5MGQUV6kDHn1qLsZOKjAldyzXB+fDtQms5RCS2qPCIyPlZnQ7bl0HPkZCQbDXKWOca9nrlwounRUTOgwqPiJyb60Lmc3BpfWj5I9tpOExpxoSuhfVTYdMc23FEJAao8IjIua36DHZkQY9fQzA69jr82+kLZatqlEdEzosKj4gUznVh+vNQqSG0vMl2mgJHSYZuj0DOTNg8z3YcEYlyKjwiUri1X8LOleGdWYGg7TSnancnlK4Is/5iO4mIRDkVHhEp3Oy/Qvla0PwG20lOl1QWOtwDa9Jh11rbaUQkikXHZLyIRKdvF8Dm2dDvuahZu/N9qSPTqURdZiUn8ulfHmdkyN52eRGJbhrhEZGzm/1XSC4PbW+3neSs9lCeD5yeXB/8L1XZZzuOiESp6PuVTUSiw96N4d1ZXR6E5EvO+JBzHd1QUv7pDGBocArDE77ihdCttuOISBTSCI+InNm8MWCC0OlntpOc0xavGhlue24NTiOZ47bjiEgU0giPSBw72whNGY4yN/ktUloNgZTLSzjVhfm304f+wfkMCszlI/dK23FEJMqo8IjIaa4LziLFHOGGBc1ZND86pq3OZbbbnHXu5dyeMJmPjqvwiMipNKUlIj/gMSw4hVVubRZ5DW2HKQLDW05f0gIbaG3W2Q4jIlFGhUdETtHWZNM8sIm3nT6AsR2nSD5xupHrleKOhAzbUUQkyqjwiMgphiVM4ZBXmk+drrajFFkuZfjI6c6gwBwqcMh2HBGJIio8IlKgIgcZGJjLx0438ihlO84FGef0JtmEGBKcZTuKiEQRFR4RKTAkOItkE2Kc09t2lAu2xqvNUrcetwQzAc92HBGJEio8IlLg5uAMlrr1WOPVth3lonzg9KRpYDMtzEbbUUQkSqjwiAgAzU0OzQKb+NDpYTvKRfvcuYKjXiI/Ck63HUVEooTuwyMiANwUnM4xL5HPnCtsR7loBynLJLcjQ4Kz+ENoKMdIOuPjznU0Rs7ogZGIJyIWaIRHREginyHBWUx223GQcrbjFIsPnJ6kmDz6BebbjiIiUUCFR0ToHVhERZPLeB9MZ500123KZrfKicXLIhLvVHhEhJuD09nmXcpMt6XtKMXGI8CHTg+6BldQ0+yyHUdELFPhEYlzVdhHj8BSPna64frsn4RP3O4ADA7MtpxERGzz179uIlJkg4NzCBqPjxz/Hbj5rVeF+W4jrgvOQvfkEYlvKjwicW5wcDbL3Lps8C63HSUiJjhdaRz4liZmi+0oImKRCo9IHEs122gd2MBnThfbUSIm3elEvhfUURMicU6FRySODQ7MwfUMXzidbUeJmH2kMMNtxbXB2Rhc23FExBIVHpF45XlcF5zFPLcp26lkO01ETXC6UsPsoaNZYzuKiFiiwiMSr7YtpX5gG5+5sX9n5XPJcNty2Es+sXhZROKRCo9IvMoaz3EvyCSno+0kEXeEUkx22zMgOI8k8m3HERELVHhE4pHrwvKPmO62Zj+X2E5TIiY4XahgDtM1kGU7iohYoMIjEo82z4FD3/G5j3dn/dAstyUHvDIMDM6zHUVELFDhEYlHKydAQimmuG1tJykx+SSQ4banT2ABiYRsxxGREqbCIxJvXBdWfQYNriaPUrbTlKh0pxPlTR5dA8ttRxGREqbCIxJvvp0Ph7ZBsyG2k5S4WW4LDnplGBD4xnYUESlhKjwi8WblBAgmQaN+tpOUuOMkkuG2pW9wAQma1hKJKyo8IvHE88KFp35vKJViO40VE51OJ3ZrrbAdRURKkAqPSDzZuggOfgvNrrOdxJqZbisOeaW5JqDdWiLxRIVHJJ6s/BQCidC4v+0k1hwnkSluW/ppWkskrqjwiMSLk9NZ9XpA6Yq201g10elERZPLFYGVtqOISAlR4RGJF9uWwv5NcT2dddKME9NaAzStJRI3VHhE4sWqz8AEofFA20msO0YSX7tpXB1cSADXdhwRKQEqPCLxYuVnkNoNylaynSQqTHbaU8UcpI3Jth1FREqACo9IPNidDXuyoem1tpNEjUy3Nce9IH2CC21HEZESoMIjEg9Wp4dfNr7Gbo4okksZ5rjN6RtYAHi244hIhKnwiMSDNROhemsoX9N2kqiS4bajXmA79c13tqOISIQl2A4gIpGTOjKdyhzgm+Rv+EvoBv4yMt12pKiS4bTj2cQ36RdYwGtODdtxRCSCNMIj4nO9gosIGI8Mt53tKFFnB5eyxK1P3+AC21FEJMJUeER8rk9gEd96lVnp1bEdJSpNdtqRFlhPVfbZjiIiEaTCI+JjpThGt8ByMpx2gLEdJypNdtsDaLeWiM+p8Ij4WPfAckqb40xx29qOErXWeTXY6FajT0CFR8TPVHhEfKxPYCEHvTLMc5vajhLFDJPd9nQJZFGOPNthRCRCVHhE/Mp16BVczNduGiFtyCxUhtOOJOPQM7DUdhQRiRAVHhG/+nY+lc3BE+t3pDCLvEbs9lK0W0vEx1R4RPxqdTrHvSDT3da2k0Q9lwDTnDb0CCwliGM7johEgAqPiF+tmchctxmHKGM7SUyY6rahvMmjnVlrO4qIRIAKj4gf7c6GPet0s8Ei+K/bkuNekF7BxbajiEgEqPCI+NGaSQBM0fqd83aY0sxzm9I7oMIj4kcqPCJ+lD0ZqjZnG5VsJ4kp09w2NAxspZbZYTuKiBQzFR4Rvzl2CDbPhYZX204Sc6a5bQDoFVhiOYmIFDcVHhG/2TAd3Hxo0Md2kpizybuM9W51emlaS8R3VHhE/GZdBiRdArU7204Sk6a5begcWEkZjtqOIiLFSIVHxE88D7KnQL0eEEy0nSYmTXPbkGxCdA1k2Y4iIsVIhUfET3augoPfQsO+tpPErPluYw56pTWtJeIzKjwifrIuI/yygRYsX6gQCcxwW4Xvx+N5tuOISDFR4RHxk+wMqNocytewnSSmfe20oZrZD9t0mKiIX6jwiPiFtqMXm0y3Na5nwvczEhFfUOER8QttRy82eyjPUq8+rP3SdhQRKSYJtgOISDHJnqzt6MVoqtOGNls/pP3Iceym/BkfkzN6YAmnEpELpREeET/wPFg3Ber31Hb0YvL1ibsu9wzqrssifqDCI+IHO1fBwa2azipGK7w6bPcqanu6iE+o8Ij4gbajR4BhmpNG98ByEgjZDiMiF0mFR8QPtB09Iqa7aVxijtDOZNuOIiIXSYVHJNYdPQib50BDTWcVt1luc/K9ID2Cuh+PSKxT4RGJdRungxtS4YmAXMqw0GtEz4AKj0isU+ERiXXZGZCcArU62U7iS5lOa5oFNlGVfbajiMhFUOERiWUnt6PrdPSIyXTTADStJRLjVHhEYpm2o0fcaq8W272K9NC0lkhMU+ERiWUnz3rSdvQIMkx3WtM9sJwgju0wInKBVHhEYtm6KVCthbajR1im25ryJo822p4uErNUeERi1cnt6BrdibhZbgtCXoCeWscjErNUeERilbajl5iDlGWh10jreERimAqPSKzSdvQSNd1pRctADlXYbzuKiFwAFR6RWKTt6CVu+ont6VcGlllOIiIXIsF2ABG5ADtXwsGt/GrPAD5YnG47TVxY4dVhp1eBnsElfOReaTuOiBSRRnhEYlF2+HT06U5ry0HiiWG604rugeUEcG2HEZEiUuERiUUntqPv4FLbSeLKdLc1Fcxh0sw621FEpIhUeERijbajWzPTbYnjGR0zIRKDVHhEYo22o1tzgHIs9hrq9HSRGKTCIxJrsidrO7pFmU5rWgc2UIkDtqOISBGo8IjEEs+D7ClQr6e2o1sy3Q0vFO8eWG45iYgUhQqPSCzZuRIOfafpLIuyvFR2eSn0DC6xHUVEikCFRySWnNiOrgXL9ngEmOG2Ct+A0NXp6SKxQoVHJJacPB095XLbSeLadCeNS00ufKdRHpFYocIjEitObkfXdJZ1M90WuJ6BdRm2o4jIeVLhEYkVGzLD29EbqPDYto8Ulnr1/2+KUUSingqPSKxYd/J09I62kwjh7elsXQh5e21HEZHzoMIjEgu0HT3qZLqtAQ/WT7MdRUTOgwqPSCzQdvSos9yrB2UqaVpLJEao8IjEguzJ4Zfajh41XAJQv1d455yr09NFop0Kj0gsyJ4C1VpqO3q0adAH8nbDNm1PF4l2CbYDiMg5HD0IW+ZClwdsJ5EfaPuuy4Jkw0uv/Y2/Ojecdj1n9EALqUTkTDTCIxLttB09au0lhWVeXXoGdXq6SLRT4RGJdtqOHtWmu2mkmXWUJ9d2FBEphAqPSDTTdvSol+m0Jmi88NlaIhK1VHhEotmOFSe2o/e1nUTOYqlXn71eOU1riUQ5FR6RaLZOp6NHO5cAM91WXBlYikHb00WilQqPSDQr2I5e3XYSKUSm05oq5iDNTY7tKCJyFio8ItHq5Hb0hhrdiXYz3VYA9NA6HpGopcIjEq20HT1m7KY8y9y69AzqBoQi0UqFRyRaZU+G5PLajh4jMt3WtDXZpGh7ukhUUuERiUaeB+umQv2e2o4eIzKdNILGo3sgy3YUETkDFR6RaHRyO7qms2LGEq8B+72y9AxoWkskGqnwiEQjbUePOeHt6S3pEVym7ekiUUiFRyQaaTt6TMp00qhq9tPMbLYdRUR+QKeli0SbowcKTkdPHZluO40UwXS3NQA9AktY4aTaDSMip9AIj0i00Xb0mLWb8ix3U3XMhEgUUuERiTYF29E72U4iFyDTTTuxPf2w7Sgi8j0qPCLRxPMgOwMa9IKgZpxjUabTmgTj0i2w3HYUEfkeFR6RaLJ9GeTu0OnoMWyJ14ADXhl6BjStJRJNVHhEokn25PBLbUePWQ7BE9vTl4ZH7EQkKqjwiEST7Ay4vA2Uq2o7iVyE6W5rqpn9sF3TWiLRQoVHJFrk7YVv52s6ywcynfD2dNZNsRtERAqo8IhEi/XTwHNVeHxgFxVZ4dZR4RGJIio8ItEiOwPKVApPaUnMy3Rbw+a54RtJioh1Kjwi0cB1w+dnNbgaAkHbaaQYZDpp4DnhG0mKiHUqPCLR4LvFkLdH01k+sshrGL6BZHaG7SgiggqPSHTIngwmAPV72U4ixcQhCPV7wrqp2p4uEgVUeESiQfZkqNEeylxqO4kUpwZ94NB3sGOF7SQicU+FR8S23J3w3SJNZ/nRyRtIrtO0lohtKjwitq2bGn7ZUKej+05KdajWErK1PV3ENhUeEduyJ0O5anBZK9tJJBIa9IYt2p4uYpsKj4hNTgjWTw2v9Qjof0dfatQP3BCs/9p2EpG4pn9hRWzauiD8m7+ms/yrZkcoVQHWfmU7iUhcU+ERsSl7Mpgg1L/KdhKJlGBCuNBmTwbXsZ1GJG6p8IjYlD0Zal8BpcrbTiKR1Kg/5O2GrQttJxGJWwm2A4jErYPbYPtyuHqU7SQSIakj0wFIwWFRcoAxf3+VF0O3FFzPGT3QVjSRuKMRHhFbsk+s6dD9d3zvIOVY4DWmd2Cx7SgicUsjPCK2rPmSLW4Vur+0AdhoO41E2FSnDU8kjuNydvMdlW3HEYk7GuERseF4Hmz4miluW8DYTiMlYKrbFoBeQY3yiNigwiNiw8bpEDpa8ENQ/G+DV52NbjV6BxbZjiISl1R4RGxYMwmSLmGe29R2EikxhmluW7oEVlKao7bDiMQdFR6Rkua64ZvQNehFvpbRxZWpbhuSTT5dAzo9XaSkqfCIlLRtSyB3OzS6xnYSKWHz3SYc9ErTS9NaIiVOhUekpK39EkxA29HjUD4JzHBb0Tu4GPBsxxGJKyo8IiVtzaTw+UplK9lOIhZMc9pQzeynhdGtCERKkgqPSEk6sBW2L4PG/W0nEUsy3TRcz+gmhCIlTIVHpCSt/TL8Uut34tZeUljsNdD9eERKmAqPSEla+yVUTIUqjW0nEYumOm1pHdgAh7bbjiISN1R4RErK8cOwYXp4dMfo7srxbKrbJvxK9mS7QUTiiAqPSEnZkAnOMa3fEdZ4tfjWqxxewC4iJUKFR6SkrJkEySlQu4vtJGKdIcNpB+u/Dp+rJiIRp8IjUhIK7q58NSQk2U4jUWCy2x5CR2D9NNtRROKCCo9ISfhuMRzeCY21O0vC5ruNoVQFWJ1uO4pIXFDhESkJq78AEwyP8IgAIRKgUX9YOwmckO04Ir6nwiNSElZ/AandoMyltpNINGkyEI7sg82zbScR8T0VHpFI27UWdq+FptfaTiLRpkFvSCilaS2REqDCIxJpqz8Pv2wy0G4OiT5JZaF+r3Dh8XSYqEgkqfCIRNqqL6BGO0i53HYSiUZNBsKBLeEz1kQkYlR4RCLpwFb4bhE0GWQ7iUSrRp9cZ58AABZGSURBVP3BBMLFWEQiRoVHJJJOrs1Q4ZGzKVsZal+hdTwiEabCIxJJqz+Hyo2gSiPbSSSaNRkIO1fA3g22k4j4lgqPSKTk7YWcWRrdkXM7uaB99US7OUR8TIVHJFLWfgmeA01VeOQcKqZCtZaa1hKJIBUekUhZ9QWk1IDL29pOIrGgyUDYPAdyd9lOIuJLKjwikXD8MKyfGv4hZoztNBILmg4CvPBduUWk2KnwiETCuqkQOqr1O3L+qrWAS+vBygm2k4j4kgqPSCSs/BTKVII6XW0nkVhhDDS7DjbOgMN7bKcR8Z0E2wFEfCf/CKz5ElrdDEH9LyZnlzry1EXKzU1l0pMdfv3H53jfuYqc0TqORKS4aIRHpLitmwL5h6HZENtJJMas8FLZ5FZlYGCu7SgivqNfP0WK24pPwtNZqd1P+w1epHCGiW4n7gmmU4FDtsOI+IpGeESK08nprKbXajpLLki604kE49I3uMB2FBFfUeERKU7ZGeHprObX204iMSrLq8tmtwoDAt/YjiLiKyo8IsWpYHdWN9tJJGaFp7W6BrLCx5OISLFQ4REpLgXTWYM1nSUXZaLTiUTjwBqdrSVSXFR4RIpLwXSWdmfJxVnm1WOLWwVWfGo7iohvqPCIFJcVn0CZyprOkmJgmOh2hA2ZcGSf7TAivqDCI1IcjufB2q+0O0uKzUSnE7j54UNoReSiqfCIFId1ms6S4rXUqw8VUyFrvO0oIr6gwiNSHJaPh7JVNJ0lxchAy5vDZ2sd2m47jEjM09i7yMU6eiA8ndX+Lk1nSbHqnVGVqckuvxv9LG8415x2XWdtiZw/jfCIXKyVn4FzDFr+yHYS8Zn1Xg2y3FQGB2fZjiIS81R4RC7W8g+gYl2o0dZ2EvGhT52upAU2kGq22Y4iEtNUeEQuxsFtsHEmtPoRGGM7jfjQF05nXM8wODDHdhSRmKbCI3Ixsj4CvPDiUpEI2E4l5rlNuS44C/BsxxGJWSo8Ihdj+QdQPQ0qN7SdRHzsU7cr9QPbaGE22o4iErNUeEQu1K61sG1peDpLJIImOR045iVwXXC27SgiMUuFR+RCLf8QTABa3Gg7ifjcQcqR6aYxODibAK7tOCIxSYVH5EJ4Xrjw1L0SLrnMdhqJAxOcLlQz++kcWGk7ikhMUuERuRCb58K+jbr3jpSYqW5bDnqluTE403YUkZikwiNyIZb8BxLLQrPrbCeROHGMJL5wruCawDeU5YjtOCIxR4VHpKiOH4YVn0Lz6yG5nO00Ekc+dHpQxhxjYHCu7SgiMUeFR6SoVn4Gx3OhzVDbSSTOLPYasN6tzk3BGbajiMQcFR6RolryTvgoidpX2E4iccfwodODjoE1OmpCpIhUeESKYl8O5MyEtKE6SkKs+NjpjuMZjfKIFJEKj0hRLHkXMJB2m+0kEqd2UpEZbituCM4E17EdRyRmqPCInC/XhaXjoF5PKF/TdhqJYx86Pbjc7IUNmbajiMSMBNsBRGLGpv/C/s3Q66mCd6WOTLcYSOLVVLct+72yVFjyDjTobTuOSEzQCI/I+Vr0NiSXh6aDbCeROHeMJCY4XWDVF5C313YckZigwiNyPg7vgZUToPUtkFjadhoR3nV6g3MMlr5rO4pITFDhETkfS8eFf7i0u8t2EhEAVnu1oWZHWPBm+Gw3ESmUCo/IuXhe+IdKrc5QrZntNCL/p/0I2JMNOf+1nUQk6qnwiJzLxhmwd334h4tINGk+BEpVgAVv2E4iEvW0S0vkXBa8AaUr6qBQiTqpv53GbxOu4PasCVyxcBx7KH/K9ZzRAy0lE4k+GuERKUzuTlj9RfjOyomlbKcROc04pxdJxuHm4HTbUUSimgqPSGEWvw1uCNoNt51E5IzWezWY6zbltuA0DK7tOCJRS4VH5GxcFxb+C1K7Q+WGttOInNU7od7UCeyke2C57SgiUUuFR+RssieH76ysxcoS5b5yO7DbS+H2YIbtKCJRS4VH5GzmvgaXXA5Nr7WdRKRQx0nkHedqegcWk2q22Y4jEpVUeETOZMdK2DgdOt4DwUTbaUTO6T+hqwkR4M7gZNtRRKKSCo/ImcwbAwmltVhZYsYuKvC524UfBTNJ4bDtOCJRR4VH5Ify9sKy96HVj6DMpbbTiJy3N0L9KWuOcXMw03YUkaijwiPyQwv/BaGj0OnntpOIFMkKry7z3CbclfAVQRzbcUSiigqPyPc5+fDNP6FeT52bJTHpjVB/aprdXB1YaDuKSFRR4RH5vhWfwKHvoPN9tpOIXJAMtz1b3CqMSPjSdhSRqKLCI3KS58F/X4YqTaFBH9tpRC6IS4B/OX3pFFgN3y6wHUckaqjwiJyUnQE7V0C3hyGg/zUkdr3n9GK/Vxb++5LtKCJRQ/+qi5z035egfC1ocaPtJCIX5TClecvpGz74dtca23FEooLxPK+w64VeFPGNLd/A2D7QfzR0vrfg3akj0y2GErlwFTnI7OQHSXc781j+mXcc5oweWMKpRCLOnO2CRnhEILx2p3RFaHuH7SQixWIfKbznXMV1gVlczm7bcUSsU+ER2bkK1qRDx59BUlnbaUSKzT9D4RGcexI0UimiwiMy/XlIKgedfmY7iUix+o7KfOp05dbg11TigO04Ilap8Eh827ESVnwavquyjpEQH3rNuY4k8vl5wue2o4hYpcIj8e3k6M4V99tOIhIRG73qfOJ25/ZgBlXYZzuOiDUqPBK/dqyAlZ9CZ43uiL+9ErqeBBzuS/jMdhQRa1R4JH5ljobkFB0jIb632avGeOdKfhycymXssR1HxAoVHolP25fDqs+0dkfixqvO9Rg87k+YYDuKiBUqPBKfMp6GUhXgCo3uSHz41qvCB05Pbgl+TU2z03YckRKXYDuASIlb/zWsnwp9/xC+2aBInPhr6HpuDM7k8YQPeCj/F+e8k7juxCx+ohEeiS+uCxlPQfna0PEe22lEStQOLuWfzgCuC86mlVlvO45IiVLhkfiS9RFsXwa9noSEZNtpRErc30PXsstL4YnEd9BxiRJPVHgkfoSOwbTfwWUtoeXNttOIWHGY0vwldCOdAqu5OrDIdhyREqPCI/Fj7muwfzP0+R0E9J++xK/3nKtY71bnfxLGkUDIdhyREqF/9SU+HPwOpv8JGg+A+r1spxGxKkQCz4V+TP3ANu4IZtiOI1IitEtL4sPk34Ibgn5/POXd59qlIuJXU9y2ZDqteSRhPJ87ndmFdiyKv2mER/wv57+QNR66PQyX1rWdRiRKGJ4J3UES+fwmcZztMCIRp8Ij/uaEYOKvwtvQuz5sO41IVMnxqvN3ZxDXB2fRyayyHUckolR4xN/mjYGdK6DfHyCpjO00IlHntdB1bHGr8LvEN7WAWXxNhUf8a+9GmPYsNOoPTa+1nUYkKh0lmVGhO2gc+JafBrWmTfxLhUf8yfPg8wchkAAD/xeMsZ1IJGpNcdvxhdOJhxI+or7ZajuOSEQYzyv0Tpu6DafEpkX/hs8egEEvQfsRZ32YdmmJhFXmAJOTHyfHu4ybjj+Dex6/D+usLYlCZ/3tViM84j8Ht8FXT0KdbtB2uO00IjFhN+V5Jv9O2gbWcVdwku04IsVOhUf8xXVhwn3gHIPBr+iOyiJF8JnbhQynLY8nfEA9853tOCLFSj8NxF+++TusnxbelVWpvu00IjHG8ET+TzhCMn9JfJVE7doSH9GdlsU/tmdBxlPQ6Bpo/xNAa3REimonFfl1/j38I+klHkt4n+dCQ21HEikWGuERf8g/Ah/dDaUqwHWvaleWyEWY7Hbg7dDV/Cwhne6BZbbjiBQLFR7xhy9Hwq5VMOT/QdnKttOIxLxnQ8NY49bkfxP/H5U5YDuOyEVT4ZHYt+htWPgv6PYINLzadhoRXzhGEg/kP8Al5PFq0iu6C7PEPBUeiW1bF0H6o1CvJ/T6re00Ir6y1qvFyPx76BxYxW8SdMCoxDYtWpbYdXgPfHAHlKsKN74BgaDtRCK+86nbjZahjfwkYRLL3bp84na3HUnkgmiER2JT6Bh8cDvk7oAfvQVlK9lOJOJbz4VuY47TjOcSX6el2WA7jsgFUeGR2OO68Ol9sGlWeJFyjXa2E4n4WogE7s9/kF1eBd5IeoFaZoftSCJFpsIjsefrZyFrPPR+ClreZDuNSFzYSwrD839FIg7/SnyBihy0HUmkSFR4JLbMHwsz/wxt74Ruv7SdRiSurPdqcPfxR6lpdvN60p/D978SiREqPBI7loyD9F9Cw74w8M+6uaCIBQu8JjyUfz9tzDp4fxjkH7UdSeS8GM/zCrte6EWREpP1UfhOynWvhNveh8RSgI6OELHlR8GveSHxn9CwH9zyNiQk244kAnDW34S1LV2i38oJ8NE9UPsKuHVcQdkREXs+cK4iiMtz2WPJGHUN9+U/RP4PfqTkjB5oKZ3I6TSlJdFt0dvw4XCo2R5+/D4klbWdSEROeNfpzZP5d9EnuJB/JP6Z0mh6S6KXCo9Er1mvwGe/CN9F+fZPIPkS24lE5Af+4/RhZP7dXBlYxjtJf6QCh2xHEjkjFR6JPq4DXz0BGb+F5teH1+xoZEckar3n9OK+/IdpbjYxPmkUNdhlO5LIaVR4JLocPQDv3gpzXoUO98CNYyEhyXYqETmHr9wO3H58JFXNfiYk/5YOZrXtSCKnUOGR6LFnPbzeB9ZPg4H/CwNf1PlYIjHkG68p1x8fxQGvLOOS/gDzX4fCdwKLlBgVHokOy8fD33vA4Z3h9TodfmI7kYhcgPVeDa4//jtmui0h/dHwMTDHtK5H7NN9eMSu44dh4q9gyX+gVie48XWoULvgsu6zIxKbArhs6LccZvwJKtQJT0/X1Ll3EnFnvQ+PRnjEno0zYUw3WPIOdH8Mhk88peyISOxyCcBVv4Hh6eDkwxt9IXM0hI7bjiZxSoVHSt7RA/D5Q/DWoPD8/p2fQ+/fQlD3wRTxnTpd4N7/QrMhkPkc/L07bJ5rO5XEIU1pSclxHVj8H5j2e8jbA1fcT5NpbTmKbkkv4ken3Wl57VfhdT0HtkCb26HXk3DJZXbCiV9pSkss8jxY/3V4UfLnD8Kl9eDuqdD3WZUdkXjSqB/cNxeu+AUsfQ9eaQuZz4fX8olEmEZ4JHI8L7zFfPoLsGUulK8NfUaFbyZ44qRzLUoWiU91zHamp2XCqs+gTGXo8gto/xMolXLen+Nc/37oLK+4pMNDpQQ5+bA6PXzzwG/nQ0oNGPBieAhbB3+KCLDJu4zUxbfS1rTlIedjekx5hgMZL/Avpz9vh/qwm/IqLFKsNMIjxefQDlj0Fix4Aw5tgwp1+J9dV/ORcyXHSbSdTkSiWEuzgV8kfEq/4AKOe0G+cjvwn1Af5nlNKOSX9kKpMMWls/7HosIjF+dYLqz+ApZ9ABsywXOgfm/o+FNo2IfU33xpO6GIxJC6ZhtDg1O4KTiDCuYw693qfO5ewWdOFzZ4lxfrc6kQ+ZIKjxSjQzsgezKs/TK8Ric/L7w+p+VNkDYUKjcoeKjW6IjIhUjmONcG53BDYCadA6sIGI8sN5XJTnu+dtPI8lLxLnLfjQqPL6nwyEU4egC2fAObZodHcb5bFH5/Sk1o3B9a3BS+S3Lg9H98VHhE5GJVZR+DgnMZFJxDmllPwHjs8lKY7qYxx2nGfK8xm72qFHXqS4XHl1R45DyFjsHOVbAjC7Ythc1zYMcK8FwIJMDlbcNbSxv1h2rNC3ZbnY0Kj4gUp0s5SPfAMq4KLqFHYBkVTS4Au7zyzHcbs8htyEqvDqvc2uzj/Hd8nYkKUUxS4fGzIm/N9Dw4sg/2boC9G8Mv96wLF5vda8ANhR+XWDZ89k3tLuG7pdZsD0lli/TcIiKRYnBpYL6jQ2AN7QNraG/WUDuwq+D6dq8iK906rPNqsMmrxkbvMja51dhGpfDRFxdJhSgqqfBEs+IqDQmESCGPFHOYiuRS1eynitnPs70qQe728NqbQ9vhwObwNNX3pdSEas3gspbhP9VawqV1IRAskewiIsXhUg7SNLCJpmYzzQKbaGo2Uddsp5TJL3jMMS+BbV4ldlCRHV5FtnuXssOrwE6vIrspz36vXPgPZTlCMtolFlPitPB4XvhP+I3/e9/Jtwu7VvB2YdcK+TyuA24+uCGufH4KCTgEcUkkRBD3xNsOCSb8esHbOJTiOKVMPqU5RmmOU8ocp9TJ1zlOaXOMUhwnxeSRwmFSTB7lOUw5c/SMfw2OZ9hNeXZ6FdjhVWSbV4kc7zI2edV4/eEfQcU6kFj6Av6CVXhEJPoZXKqxj7qB7dQxO0g127nc7KGa2Uc19nGZ2XtKIfq+Y14i+ynLfq8chyhDnpfMEZLJI5kjXjKHKVXweh7JHCOJkBcknwReGdYRgkkQSIRgYvj14PdeDySG1z6aM/0JnnhpTr8WCP7gfRdWyHwqSgvP9iwY25diLyM+lOclc5REjpDMUS+JoyRxiDIc9MpwwCvLAcpy8MTLk2/v9Cqw06vAHsoXy/CtiIg/eaRwmMvMPiqZg5TnMBVMLhXIpYI5THlyqWBySSGPMuYYpTlGGY5S2hyjzIlfRgMmmn7+/OBn/mmFyJzftWL9WAPt7oT+z52WtphdWOExxnwJVI5EIp+qDOy2HUJOo+9L9NH3JDrp+xJ99D0pmt2e5/U/04VzjfBIERhjFnie1952DjmVvi/RR9+T6KTvS/TR96T4aJ5DREREfE+FR0RERHxPhad4/cN2ADkjfV+ij74n0Unfl+ij70kx0RoeERER8T2N8IiIiIjvqfBEiDHmUWOMZ4zRtn7LjDF/MsasNsYsM8Z8YoypYDtTPDPG9DfGrDHGrDPGjLSdJ94ZY2oZY742xqw0xqwwxjxkO5OEGWOCxpjFxpgvbGfxAxWeCDDG1AL6ApttZxEAMoAWnue1AtYC/2M5T9wyxgSBvwHXAM2A24wxzeyminsh4FHP85oBnYH79T2JGg8Bq2yH8AsVnsh4CfgVfr71cwzxPG+y53knTkRlLlDTZp441xFY53neBs/zjgPvAddZzhTXPM/b5nneohOvHyL8A7aG3VRijKkJDARet53FL1R4ipkx5jpgq+d5S21nkTMaAUyyHSKO1QC2fO/tb9EP16hhjEkF2gDz7CYR4GXCvzi7toP4RYLtALHIGDMFuOwMl54AfkN4OktKUGHfE8/zJpx4zBOEh+/fKclsIrHAGFMO+Ah42PO8g7bzxDNjzCBgp+d5C40xPW3n8QsVngvged7VZ3q/MaYlUBdYasKHptUEFhljOnqet70EI8ads31PTjLGDAcGAb093YvBpq1Are+9XfPE+8QiY0wi4bLzjud5H9vOI3QFBhtjBgClgBRjzH88zxtmOVdM0314IsgYkwO09zxPB79ZZIzpD/wv0MPzvF2288QzY0wC4YXjvQkXnfnAjz3PW2E1WBwz4d/O3gL2ep73sO08cqoTIzyPeZ43yHaWWKc1PBIPXgUuATKMMUuMMWNsB4pXJxaP/wL4ivDi2A9UdqzrCtwO9Drx/8eSEyMLIr6iER4RERHxPY3wiIiIiO+p8IiIiIjvqfCIiIiI76nwiIiIiO+p8IiIiIjvqfCIiIiI76nwiIiIiO+p8IiIiIjv/X+BmGHqmiUoSAAAAABJRU5ErkJggg==\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",
"\n",
"I highly encourage you to play around with the code in this notebook - this way, you can 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",
"\n",
"In my next post, I will discuss the Gibbs sampler - a special case of the Metropolis-Hastings algorithm that allows you to approximately sample from a multivariate distribution by sampling from the conditional distributions.\n",
"\n",
"Thanks for reading&mdash;go forward and sample!"
]
}
],
"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
}
You can’t perform that action at this time.