{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Markov chain Monte Carlo\n", "\n", "A major takeaway from the previous section is the inherent difficulty in calculating or estimating posterior distributions for use in Bayesian inference. The two alternative strategies to obtaining posteriors for moderate to large models involve either analytic **approximations** or stochastic **sampling**. Approximations are usually valid conditional on assumptions regarding the true posterior distribution, which are typically impossible to validate. Direct sampling strategies rely on our ability to sample from the posterior distribution, and this is frequently not possible. Indirect sampling methods, such as rejection sampling, can be plagued with sampling efficiency issues.\n", "\n", "The sampling approaches we have introduced so far have each attempted to obtain *independent* samples from the posterior distribution. It turns out, however, that it is possible to generate samples from the posterior distribution using a *dependent* sampling algorithm, and despite the dependence of the samples, one may extract valid inference from them. A class of algorithms called **Markov chain Monte Carlo** yields a Markovian sample (explained below) which, provided that certain conditions are satisfied, is guaranteed to be indistinguishable from a sample drawn from the true posterior itself." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Markov Chains\n", "\n", "A Markov chain is a special type of *stochastic process*. The standard definition of a stochastic process is an ordered collection of random variables:\n", "\n", "$$\\{X_t: t \\in T\\}$$\n", "\n", "where $t$ is frequently (but not necessarily) a time index. If we think of $X_t$ as a state $X$ at time $t$, and invoke the following dependence condition on each state:\n", "\n", "$$Pr(X_{t+1}=x_{t+1} | X_t=x_t, X_{t-1}=x_{t-1},\\ldots,X_0=x_0) = Pr(X_{t+1}=x_{t+1} | X_t=x_t)$$\n", "\n", "then the stochastic process is known as a Markov chain. This conditioning specifies that the future depends on the current state, but not past states. Thus, the Markov chain wanders about the state space,\n", "remembering only where it has just been in the last time step. The collection of transition probabilities is sometimes called a *transition matrix* when dealing with discrete states, or more generally, a\n", "*transition kernel*.\n", "\n", "In the context of Markov chain Monte Carlo, it is useful to think of the Markovian property as “mild non-independence”. MCMC allows us to indirectly generate independent samples from a particular posterior distribution.\n", "\n", "### Markov Chain Jargon\n", "\n", "Before we move on, it is important to define some general properties of Markov chains. They are frequently encountered in the MCMC literature, and some will help us decide whether MCMC is producing a useful sample from the posterior.\n", "\n", "---\n", "\n", "**Homogeneity**: A Markov chain is homogeneous at step $t$ if the transition probabilities are independent of time $t$. \n", " \n", "**Irreducibility** A Markov chain is irreducible if every state is accessible in one or more steps from any other state. That is, the chain contains no absorbing states. This implies that there is a non-zero probability of eventually reaching state $k$ from any other state in the chain. \n", " \n", "**Recurrence** States which are visited repeatedly are *recurrent*. If the expected time to return to a particular state is bounded, this is known as *positive recurrence*, otherwise the recurrent state is *null recurrent*. Further, a chain is *Harris recurrent* when it visits all states $X \\in S$ infinitely often in the limit as $t \\to \\infty$; this is an important characteristic when dealing with unbounded, continuous state spaces. Whenever a chain ends up in a closed, irreducible set of Harris recurrent states, it stays there forever and visits every state with probability one. \n", " \n", "**Stationarity** A stationary Markov chain produces the same marginal distribution when multiplied by the transition kernel. Thus, if $P$ is some $n \\times n$ transition matrix: \n", " \n", "$$\\mathbf{\\pi} P = \\mathbf{\\pi}$$\n", " \n", "for Markov chain $\\pi$. Thus, $\\pi$ is no longer subscripted, and is referred to as the *limiting distribution* of the chain. In MCMC, the chain explores the state space according to its limiting marginal distribution. \n", " \n", "**Ergodicity**: Ergodicity is an emergent property of Markov chains which are irreducible, positive Harris recurrent and aperiodic. Ergodicity is defined as: \n", " \n", "$$\\lim_{n \\to \\infty} Pr^{(n)}(\\theta_i \\rightarrow \\theta_j) = \\pi(\\theta) \\quad \\forall \\theta_i, \\theta_j \\in \\Theta$$ \n", " \n", "or in words, after many steps the marginal distribution of the chain is the \n", "same at one step as at all other steps. This implies that our Markov chain, \n", "which we recall is dependent, can generate samples that are independent if \n", "we wait long enough between samples. If it means anything to you, \n", "ergodicity is the analogue of the strong law of large numbers for Markov \n", "chains. For example, take values $\\theta_{i+1},\\ldots,\\theta_{i+n}$\n", "from a chain that has reached an ergodic state. A statistic of interest can \n", "then be estimated by:\n", " \n", "$$\\frac{1}{n}\\sum_{j=i+1}^{i+n} h(\\theta_j) \\approx \\int f(\\theta) h(\\theta) d\\theta$$\n", "\n", "---\n", "\n", "## Why MCMC Works: Reversible Markov Chains\n", "\n", "Markov chain Monte Carlo simulates a Markov chain for which some function of interest\n", "(*e.g.* the joint distribution of the parameters of some model) is the unique, invariant limiting distribution. An invariant distribution with respect to some Markov chain with transition kernel $Pr(y \\mid x)$ implies that:\n", "\n", "$$\\int_x Pr(y \\mid x) \\pi(x) dx = \\pi(y).$$\n", "\n", "Invariance is guaranteed for any *reversible* Markov chain. Consider a Markov chain in reverse sequence:\n", "$\\{\\theta^{(n)},\\theta^{(n-1)},...,\\theta^{(0)}\\}$. This sequence is still Markovian, because:\n", "\n", "$$Pr(\\theta^{(k)}=y \\mid \\theta^{(k+1)}=x,\\theta^{(k+2)}=x_1,\\ldots ) = Pr(\\theta^{(k)}=y \\mid \\theta^{(k+1)}=x$$\n", "\n", "Forward and reverse transition probabilities may be related through Bayes theorem:\n", "\n", "$$\\frac{Pr(\\theta^{(k+1)}=x \\mid \\theta^{(k)}=y) \\pi^{(k)}(y)}{\\pi^{(k+1)}(x)}$$\n", "\n", "Though not homogeneous in general, $\\pi$ becomes homogeneous if:\n", "\n", "- $n \\rightarrow \\infty$\n", "\n", "- $\\pi^{(i)}=\\pi$ for some $i < k$\n", "\n", "If this chain is homogeneous it is called reversible, because it satisfies the ***detailed balance equation***:\n", "\n", "$$\\pi(x)Pr(y \\mid x) = \\pi(y) Pr(x \\mid y)$$\n", "\n", "Reversibility is important because it has the effect of balancing movement through the entire state space. When a Markov chain is reversible, $\\pi$ is the unique, invariant, stationary distribution of that chain. Hence, if $\\pi$ is of interest, we need only find the reversible Markov chain for which $\\pi$ is the limiting distribution.\n", "This is what MCMC does!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gibbs Sampling\n", "\n", "The Gibbs sampler is the simplest and most prevalent MCMC algorithm. If a posterior has $k$ parameters to be estimated, we may condition each parameter on current values of the other $k-1$ parameters, and sample from the resultant distributional form (usually easier), and repeat this operation on the other parameters in turn. This procedure generates samples from the posterior distribution. Note that we have now combined Markov chains (conditional independence) and Monte Carlo techniques (estimation by simulation) to yield Markov chain Monte Carlo.\n", "\n", "Here is a stereotypical Gibbs sampling algorithm:\n", "\n", "1. Choose starting values for states (parameters):\n", " ${\\bf \\theta} = [\\theta_1^{(0)},\\theta_2^{(0)},\\ldots,\\theta_k^{(0)}]$\n", "\n", "2. Initialize counter $j=1$\n", "\n", "3. Draw the following values from each of the $k$ conditional\n", " distributions:\n", "\n", "\\begin{aligned}\n", "\\theta_1^{(j)} &\\sim& \\pi(\\theta_1 | \\theta_2^{(j-1)},\\theta_3^{(j-1)},\\ldots,\\theta_{k-1}^{(j-1)},\\theta_k^{(j-1)}) \\\\\n", "\\theta_2^{(j)} &\\sim& \\pi(\\theta_2 | \\theta_1^{(j)},\\theta_3^{(j-1)},\\ldots,\\theta_{k-1}^{(j-1)},\\theta_k^{(j-1)}) \\\\\n", "\\theta_3^{(j)} &\\sim& \\pi(\\theta_3 | \\theta_1^{(j)},\\theta_2^{(j)},\\ldots,\\theta_{k-1}^{(j-1)},\\theta_k^{(j-1)}) \\\\\n", "\\vdots \\\\\n", "\\theta_{k-1}^{(j)} &\\sim& \\pi(\\theta_{k-1} | \\theta_1^{(j)},\\theta_2^{(j)},\\ldots,\\theta_{k-2}^{(j)},\\theta_k^{(j-1)}) \\\\\n", "\\theta_k^{(j)} &\\sim& \\pi(\\theta_k | \\theta_1^{(j)},\\theta_2^{(j)},\\theta_4^{(j)},\\ldots,\\theta_{k-2}^{(j)},\\theta_{k-1}^{(j)})\n", "\\end{aligned}\n", "\n", "4. Increment $j$ and repeat until convergence occurs.\n", "\n", "As we can see from the algorithm, each distribution is conditioned on the last iteration of its chain values, constituting a Markov chain as advertised. The Gibbs sampler has all of the important properties outlined in the previous section: it is aperiodic, homogeneous and ergodic. Once the sampler converges, all subsequent samples are from the target distribution. This convergence occurs at a geometric rate." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example: Inferring patterns in UK coal mining disasters\n", "\n", "Let's try to model a more interesting example, a time series of recorded coal mining \n", "disasters in the UK from 1851 to 1962.\n", "\n", "Occurrences of disasters in the time series is thought to be derived from a \n", "Poisson process with a large rate parameter in the early part of the time \n", "series, and from one with a smaller rate in the later part. We are interested \n", "in locating the change point in the series, which perhaps is related to changes \n", "in mining safety regulations." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "disasters_array = np.array([4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,\n", " 3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,\n", " 2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,\n", " 1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,\n", " 0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,\n", " 3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,\n", " 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAu8AAAEBCAYAAADIJUs+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHv5JREFUeJzt3XuYJVV97vHvOwMIAma4KwwwgFwMJPFyRkUZUYwaL6gxHo13jSZBT2LIhagkJ6CJxghGiKgEhUiEoELQxEuUeFACiooRRFG5ySADchlgRDiCOPzyR1XLZk93T++Z7uqu7u/neeqZ3bVW1Vp7r+6ed69eVTtVhSRJkqS5b9Fsd0CSJEnS1BjeJUmSpJ4wvEuSJEk9YXiXJEmSesLwLkmSJPWE4V2SJEnqCcO7pAUtyauTXDBN53pZknOmu+7GSvLhJH/TPl6R5PIu2pUkTT/Du6QZk6SSPHxo39FJTmsfPznJqoGyzZKcneTLSR7SdX83VlWdXlVPn+6606mqzq+qfWfq/INvFLqW5IAkn0+yOsk6H2KSZFmSzya5PcmNSU5IsslAeSW5K8md7fahgbKnJPlikh8nWTmFvkxaP8kTknw9yU+SXJrkoKHyHZL8S5I1bX9PHyg7NsmV7bHfT/LKqb9KkvrO8C5pTkjyIOBsYAnw9Kq6Y5a7pFkwGKY3wL3Ax4HXTlD+fuBm4GHAI4GDgTcM1fm1qtqq3V43sP8u4BTgiCn2ZcL6SbYF/h04hub7/V3Ap5JsM1DtbOBGYHdgR+DYoXMfCvwS8Crg+CRPmGK/JPWc4V3SrEvyYOBTwKbAs6vqrgnqbZHk3UmubWc0L0iyRVv23CSXtTOVX0ryiIHj3pzk6nam8rtJfnOK/VrWzsa+Jsl17QzoYUmWt7Ola5KcMFD/AUtw2mMPa2dJb0/yviTZgLqL2+e9Osk1Sf6grT9u0E3yqCTfbJ/vx4DNB8qG/9rxpiTXt3UvT/LUdv9jk1zYPscftbPUm7VlSfKeJDe343BpO+v9e8DLgD9vZ64/1dbfOcm/Jrml7f8bB9o/OslZSU5Lcgfw6rbtbyS5I8lNSf5+KuNVVZdX1cnAZRNU2QP4eFXdXVU3Ap8D9p/iub9eVR8BfjAN9Z8A3FRVZ1bV2qo6DbgFeAFAkqcDuwJHVNWPq+reqrp44NxHVdX3q+q+qvoacD5w4FT6Jan/DO+SZtuDgP8A7gaeW1U/naTuscBjaMLPtsCfA/cl2Qc4Azgc2AH4LM1M5mbtcVcDK2hmKt8KnJbkYSP08XHA3sCLgeOAvwB+nSb4vSjJwZMc+xxgOfBrwIuAZ2xA3d8FnkkzW/xo4PkTnaB9zp8EPkLzGp0J/NYEdfcF/gBYXlVbt+2tbIvXAn8MbE8TDJ/K/bPUTweeBOxDM3P8YuDWqjoJOB14VztzfWiSRTRvzL4F7NKe5/Akg6/D84Cz2nOdDhwPHF9VDwH2oplNnw7HA7+d5MFJdqF5TT83VOe/0iypOTvJsmlqd1jabXjfAe3jxwOXA6cmuTXJRRN9j7VvXpcz8RsWSfOM4V3SbNuaJhyeWlX3TFSpDYG/A/xRVV3fzlh+pT3mxcBnquo/q+pempC/BU3Ip53hvKGdqfwYcCXw2BH6+NftbO05NEsWzqiqm6vqeppZz0dNcuw7q2pNVf0Q+CJNAB+17otowuyqqrodeOck53g8zV8wjmtnbM8CLpqg7lqaN0+/nGTTqlpZVVcDVNV/V9VXq+rnVbUS+EeaZSbQLE/ZGtgPSFV9r6p+NEEby4EdquptVfWzqvoB8EHgtwfqXFhVn2zH56ft+R+eZPuqurOqvjrJ8x3FeTRvuO4AVgHfoHmjM+ZgYFn7vG4APr2Ry3gm8hVg5yQvSbJpklfRvEl5cFu+lOYN0heBhwLvBv4tyfbjnOtEmjdGn5+BfkqagwzvkmbSWpogOWhTmnA2ZjVNkDt1aDZ22PY0yz+uHqdsZ+DasS+q6j7gOpqZXpK8Mskl7RKQNTQznOMFoYncNPD4p+N8vdUkx9448Pj/b2DdnWmez5jBx8N2Bq6vqsELNq8dr2JVXUXz14qjgZuTfDTJzgBJ9kny6XYW+g7gHbSvWVWdC5wAvA+4KclJmfgC491pguqagdf/SGCnSZ7Pa2lm9b/fzjo/Z5LnOyXtm7/P06wl37J9LtsAfzdWp6r+q32DsQb4I5plNo8Y53TD5z4y91/keuL66lfVrTR/bfgTmu+l3wC+QPOGAprvqZVVdXL7BuyjNK/RE4faPYbme/lFQ+MtaR4zvEuaST+kmckctAdDYbKqzqZZGnJWkqdMcK7VNEtr9hqn7AaakAg0a7Jp1gxfn2R3mpnePwC2q6olwHdYd9nCXPYjmtnYMbuup+4uY+vlW7tNVLmq/qWqDqJ5/Yr7w+wHgO8De7fLV45k4DWrqn+oqsfQzGTvw/0XZg6HyOuAa6pqycC2dVU9a7AbQ326sqpeQnOh5t/RfF9sOclznoptaV63E6rqnjZA/xPwrEmOKabwfVJV7xi4yPWwqXSmqs6rquVVtS3wCmBf4Ott8aWs+zo+QJK30iz78eJuaYExvEuaSR8D/jLJ0iSLkvw6zV0yzhquWFVn0ATsf0vyxHHK76O5e8fftxdALk5yYJq71HwceHaSpybZFPhT4B6a5Qlb0gShWwCSvIb71xb3xceBP0qyS5IlwJsmqXsh8HPgjUk2SfICJlgilGTfJIe0r+HdNDO+a9virWmWl9yZZD/g9QPHLU/yuPa1vqs9duy4m4A9B5r5OnBHe2HsFu24HZBk+URPIMnLk+zQjvmadvfatmxlkldPcFySbA6MXVi7efvcqKrVwDXA69vXZQnNnVq+1dbdP8kj2/5tRbNU5Xrge235ovbcm7ZNbT5wTcV4fZm0fpqLijdt/2JxLLCqqsaWvnwC2CbJq9r+vJDmr0hfbo99C/BS4GntmxBJC4jhXdJMehtNgL4AuJ3mlngvq6rvjFe5qk6lCd6fSTJe4Pwz4Ns0a7hvo5mVXVRVlwMvB95LM0N/KHBouwTiuzRB7EKaYPkrtCGoRz4InEMzI3sxzQW5P+f+wPwLVfUzmruWvJrmNX8xzVKR8TyIZv38apolOzvSzLBD81q/FPhJ2/7HBo57SLvvdpq/otzK/bcyPJlmDf2aJJ+sqrU04/FImvC8GvgQzcXDE/kN4LIkd9JeZFpVd7fhdztgojXwu9O8ARm7ePOnNBd+jnlBe+5bgKtoXsM/bst2ap/jHTR3iFkGPKe9hgKaC3R/SvPa79Y+nuxDttZX/89pXovraG5d+Ys7IFXVbcBzacbgx8Cbgee1b0CgWcK0G3DlwHKdI5G0IMRlcpLUL0meCZxYVbuvt/I8kuaDjP5Pu6RGkhYkw7skzXHt7QCfQjNzuxPwr8BXq+rwWe2YJKlzhndJmuPSfIjVeTS3MPwp8BmaW2Z6oaIkLTCGd0mSJKknvGBVkiRJ6omZ+OS4adXe5ms5zb2L17mzgiRJktRDi2nuNnXRZJ8wPmzOh3ea4H7+bHdCkiRJmgEraG6pPCV9CO8/Ajj//PNZunTp+upKkiRJc96qVatYsWIFtFl3qvoQ3tcCLF26lGXLls1yVyRJkqRpNdKycC9YlSRJknrC8C5JkiT1RGfhPcnmST6Q5Mok305yUldtS5IkSfNBl2ve3wXcDexTVZVkpw7bliRJknqvk/CeZCvglcDSaj/Stapu6qJtSZIkab7oauZ9L+BW4KgkTwHuBP6yqh5wT8skS4AlQ8d6f0hJkiSJ7sL7JsCewMVVdUSSxwGfSvLwqrpjoN7hwFEd9WnWLT/m3HX2XXTEIbPQk/vNxT5JkiSp0dUFq9cCPwfOAKiqrwGrgX2G6h0H7DG0reioj5IkSdKc1snMe1WtTvJF4GnAOUn2AXYErhqqtwZYM7gvSRddlCRJkua8Lu82cxhwSpJ3A/cCr2jDuiRJkqQp6Cy8V9UPgCd31Z4kSZI03/gJq5IkSVJPGN4lSZKknjC8S5IkST1heJckSZJ6wvAuSZIk9YThXZIkSeoJw7skSZLUE4Z3SZIkqScM75IkSVJPGN4lSZKknjC8S5IkST1heJckSZJ6wvAuSZIk9YThXZIkSeoJw7skSZLUE4Z3SZIkqScM75IkSVJPGN4lSZKknjC8S5IkST1heJckSZJ6wvAuSZIk9cQmXTWUZCVwd7sBvKmqPt9V+5IkSVLfdRbeWy+squ903KYkSZI0L7hsRpIkSeqJrmfeT08S4ALgyKpaM1iYZAmwZOiYpV11TpIkSZrLugzvK6rquiQPAo4DTgBePlTncOCoyU6y/Jhz19l30RGHTFcfNY02ZKzm6vhO1K8u+jtXXxNJktS9zpbNVNV17b/3AO8HnjhOteOAPYa2FV31UZIkSZrLOpl5T7IlsElV/bhdNvPbwCXD9dplNMNLabrooiRJkjTndbVsZifgX5MsBhYD3wXe0FHbkiRJ0rzQSXivqh8Aj+qiLUmSJGm+8laRkiRJUk8Y3iVJkqSeMLxLkiRJPWF4lyRJknrC8C5JkiT1hOFdkiRJ6gnDuyRJktQThndJkiSpJwzvkiRJUk8Y3iVJkqSeMLxLkiRJPWF4lyRJknrC8C5JkiT1hOFdkiRJ6gnDuyRJktQThndJkiSpJwzvkiRJUk8Y3iVJkqSeMLxLkiRJPWF4lyRJknrC8C5JkiT1hOFdkiRJ6onOw3uSo5JUkgO6bluSJEnqs07De5JHA48Hfthlu5IkSdJ80Fl4T/Ig4H3AG4Dqql1JkiRpvtikw7beBpxWVdckGbdCkiXAkqHdS2e6Y5IkSVIfdBLekxwILAfevJ6qhwNHjVfw3H/8Cpf+7bJJD15+zLnr7LvoiEOm1MeNPdd0tr0hZrv9UXXR34namM3Xaq6O01ztlyTNF/7fo+ky5WUzSZ6SZI/28cOSnJrklCQPncLhBwP7AdckWUkzm/75JE8fqnccsMfQtmKqfZQkSZLms1HWvL8fWNs+fjewKc3a9ZPWd2BVvbOqdq6qZVW1DFgFPKOqzhmqt6aqVg5ubV1JkiRpwRtl2cwuVfXDJJsAzwB2B34G3DAjPZMkSZL0AKOE9zuS7AQcAHy3qu5MshnNDPxI2tl3SZIkSSMYJby/F7gI2IzmwlKAJwLfn+5OSZIkSVrXKOH9GOATwNqqurrddz3wumnvlSRJkqR1TCm8J1kM3Aksqap7xvZX1RUz1TFJkiRJDzSlu81U1VrgCmC7me2OJEmSpImMsmzmdODTSY6nuX1jjRVU1bp3/5ckSZI0rUYJ769v/z16aH8Be05LbyRJkiRNaMrhvar2mMmOSJIkSZrcKJ+wSpJNk6xI8uL26y2TbDkzXZMkSZI0aMrhPcmv0Fy0+kHg5Hb3wcApM9AvSZIkSUNGmXn/APBXVbUfcG+77zzgoGnvlSRJkqR1jBLe9wdOax8XQFXdBWwx3Z2SJEmStK5RwvtK4DGDO5I8FrhqOjskSZIkaXyj3Cry/wKfSXIisFmStwCHAb87Iz2TJEmS9ABTnnmvqk8DzwR2oFnrvjvwgqo6Z4b6JkmSJGnAlGfek/zvqjoTeMPQ/hdW1VnT3jNJkiRJDzDKmveTJ9h/0nR0RJIkSdLk1jvznmTP9uGiJHsAGSjeE7h7JjomSZIk6YGmsmzmKppbQwa4eqjsRuDoae6TJEmSpHGsN7xX1SKAJOdV1cEz3yVJkiRJ4xnlbjMPCO5J9kyy+/R3SZIkSdJ4phzek5yR5Ant49cAlwHfTfLameqcJEmSpPuNcreZpwLfaB//CfDrwGOBN0/l4CSfTPKtJBcnOT/JI0frqiRJkrSwjfIJq5tV1c+S7AJsW1VfBkiy0xSPf1VV/bg95nnAKcCjR+qtJEmStICNEt4vSfIWmk9W/QxAG+TvmMrBY8G99UvAfSO0LUmSJC14o4T31wJ/DdwLHNHuOxA4faonSPIh4Ok0t538jXHKlwBLhnYvHaGPkiRJ0rw15fBeVVcDLx3adxZw1gjneB1AklcAxwDPGqpyOHDUVM83G5Yfc+46+y464pBO2uii7cnMZr9m+7mPaq7217FaV9/6uyGm6zkuhNdKM2uu/v8m9ckoM+9j69sfC2zPwCetVtUpo5ynqj6S5KQk21XVrQNFxwEfHqq+FDh/lPNLkiRJ89GUw3uS5wOnAVcC+9PcKvIA4AKai08nO3YrYJuquq79+lDgtnb7hapaA6wZOnaqXZQkSZLmtVFm3v8GeE1VnZnk9qp6VHu/9/2ncOyWwJlJtgTW0oT2Q6uqRu+yJEmStDCNEt53q6ozh/adCtwI/NlkB1bVTcDjR+ybJEmSpAGjfEjTzQP3dF+Z5EBgL2Dx9HdLkiRJ0rBRwvsHgYPax+8Bvgh8C3j/dHdKkiRJ0rpGuVXk3w08/uckXwK2rKrvzUTHJEmSJD3QKDPvw/YCdpiujkiSJEma3JTDe5Lzkjyxffwm4KPAGUmOnKnOSZIkSbrfKDPvBwBfbR//LvBkmjvIHDbNfZIkSZI0jlFuFbkIqCR7ARlb655kmxnpmSRJkqQHGCW8XwCcADwM+ARAG+RXz0C/JEmSJA0ZZdnMq4E1wKXA0e2+/YDjp7dLkiRJksYzyq0ibwWOHNr3mWnvkSRJkqRxTRrek/xFVb29ffy2iepV1V9Nd8ckSZIkPdD6Zt6XDjzedSY7IkmSJGlyk4b3qnr9wJfHACuAbYHbgAuq6rIZ7JskSZKkAetd854kwMnAK4HrgRuAXYCdk3wE+J2qqhntpSRJkqQp3W3m92g+kOnAqtq9qg6sqt2AA2lm4n9/BvsnSZIkqTWV8P4K4I1VddHgzvbrw9tySZIkSTNsKuH9l4HzJig7ry2XJEmSNMOmEt4XV9VPxito94/yQU+SJEmSNtBUPqRp0yRPAbIR55AkSZK0kaYSvG8GTllPuSRJkqQZtt7wXlXLOuiHJEmSpPXoZL16ku2SfDbJ5UkuTXJ2kh26aFuSJEmaL7q62LSAd1XVvlX1q8DVwDs7aluSJEmaFzoJ71V1W1V9aWDXV4Hdu2hbkiRJmi86v1NMkkXA64F/H6dsCbBkaPfSLvolSZIkzXWzcZvH9wJ3AieMU3Y4cFS33VmYlh9z7jr7LjrikFnoicaMOiaT1Z+u8Z3O75Mu+qvZM9EYOrZTN59eq9l8Ln1rez6Nu7rRaXhPciywN3BoVd03TpXjgA8P7VsKnD/DXZMkSZLmvM7Ce5K3A48Bnl1V94xXp6rWAGuGjuugd5IkSdLc10l4T7I/cCRwBfCVNpBfU1W/2UX7kiRJ0nzQSXivqssAp9AlSZKkjdDVfd4lSZIkbSTDuyRJktQThndJkiSpJwzvkiRJUk8Y3iVJkqSeMLxLkiRJPWF4lyRJknrC8C5JkiT1hOFdkiRJ6gnDuyRJktQThndJkiSpJwzvkiRJUk8Y3iVJkqSeMLxLkiRJPWF4lyRJknrC8C5JkiT1hOFdkiRJ6gnDuyRJktQThndJkiSpJwzvkiRJUk8Y3iVJkqSe6CS8Jzk2yTVJKskBXbQpSZIkzTddzbx/EngScG1H7UmSJEnzziZdNFJVFwAk6aI5SZIkaV7qJLxPVZIlwJKh3Utnoy+SJEnSXDOnwjtwOHBUV40tP+bcdfZddMQhXTWvaTCdY+j3w2i6eL0mamPU/ZOdazqN2kbf+tuF8foEGz7uo7SzoWPVhbn489aF6fwZ2ZBzzcXXZDJzdaxm02zmhJl6TeZaeD8O+PDQvqXA+d13RZIkSZpb5lR4r6o1wJrBfa6TlyRJkhpd3SryH5KsoplF/0KSy7poV5IkSZpPurrbzBuBN3bRliRJkjRf+QmrkiRJUk8Y3iVJkqSeMLxLkiRJPWF4lyRJknrC8C5JkiT1hOFdkiRJ6gnDuyRJktQThndJkiSpJwzvkiRJUk8Y3iVJkqSeMLxLkiRJPWF4lyRJknrC8C5JkiT1hOFdkiRJ6gnDuyRJktQThndJkiSpJwzvkiRJUk8Y3iVJkqSeMLxLkiRJPWF4lyRJknrC8C5JkiT1RGfhPck+SS5MckX7795dtS1JkiTNB13OvJ8IvK+q9gHeB/xjh21LkiRJvddJeE+yI/Bo4Ix21xnAo5Ps0EX7kiRJ0nywSUft7ApcX1VrAapqbZIb2v23jFVKsgRYMnTs7gA/W3MLK1eu5J7bblzn5CtXrgSYtGw8G3KuUfdPZjrbmK7n3tVrMl3nsr/2d6ba2BBd/A6azn7N5u/MUdpY37mm8/fvdNXfEF19P4za/mz+/+bvoNF00caobc+2ufo7HmDVqlVjDxeP0o9U1Sj1N0iSxwD/XFX7D+z7LvDyqvrmwL6jgaNmvEOSJEnS3LCiqi6YauWuwvuOwBXAdu2s+2LgVmDvqlrfzPtmwJ7AlcDaGe9s95YC5wMrgFXrqav5w3FfmBz3hccxX5gc94Vp1HFfDDwMuKiq7plqI50sm6mqm5NcArwEOK399+LB4N7WWwOsGecUV8x8L2dHkrGHq6pq5Sx2RR1y3Bcmx33hccwXJsd9YdrAcb961Ha6WvMOcBhwapK/Am4HXtlh25IkSVLvdRbeq+r7wOO6ak+SJEmab/yEVUmSJKknDO+zbw3wVsZf66/5y3FfmBz3hccxX5gc94Wpk3Hv5G4zkiRJkjaeM++SJElSTxjeJUmSpJ4wvE+zJMcmuSZJJTlgYP9zklyc5JIklyZ5wUDZ5kk+kOTKJN9OctJA2T5JLkxyRfvv3l0/J63fBo77ZGWOew9MMu7PTvLN9uf5vCR7DJRNOLaOez+MOu5Jtkvy2SSXtz/rZyfZYeC4xyf5Vjvu57QfbKg5ZEN+1gfqHDXOcY55D2zg7/iZz3RV5TaNG3AQsCuwEjig3Reae9uPff2rwE+ARe3X/wC8h/uvQdhp4HznAi9vH78cOHe2n6Pbxo/7FL4nHPcebBOM+zbAamCfgfH73MAxE46t496PbdRxB7YFnjxw/DHAye3jAFcBB7Vf/yVwymw/R7eNG/OB4x4N/Adw7dD/DY55D7YN/B0/45lu1l+Y+bqNE+JuBZ7Yfv0k4Ir28VY0VyVvNc45dmzLFrdfL26/3mG2n5/bRo/7ZGWOe8+2oXFfDlw2ULYtUMD2k42t496/barjPs5xvwV8YeC47wyUbQ/cOdvPzW3jxxx4EHAhsMc4xznmPdpG+B3fSaZz2UwHqhmlFwH/luRa4JPAq9rivWhC3FFJvpHkS0kOast2Ba6vqrXtedYCN7T7NcdNNu7r+Z5w3PvtCuChSZa3X7+s/Xc3Jh9bx73fJhv3X0iyCHg98O8D5deOlVfVamBRkm1ntruaBusb87cBp1XVNUPHOeb9Ntm4d5LpDO8dSLIJ8BbgeVW1O3Ao8LEkW9F8yu2ewMVV9b+ANwFnJ3nIrHVY02KycV/P94R6rKp+DLwYeE+Sb3D/bMu9s9oxzagRxv29wJ3ACd32UNNtsjFPciDNDO37Z7GLmgHr+VnvJNNtMp0n04QeCexcVV8GqKovJ7kLeARwDfBz4Iy27GtJVgP7AD8EdkmyuKrWJlkM7AxcNxtPQiObbNxrkrJrcdx7raq+AHwBIMlOwBHAD4AtmXhsM0mZemCScafddyywN3BoVd3X7v4hsPtAne2bU9VtXfVbG26SMf9DYD/gmiQAS4HPJ3kNjnnvTTLuW9BBpnPmvRurgKVJ9gVI8gjgocDV7Z/Lvgg8rS3bh+Zd3FVVdTNwCfCS9jwvoXk3d0vH/deGmXDcJytz3PsvyUPbfxcB7wBOrKq7Jhtbx73/Jhr3dt/bgccAz6+qewYO+29gi4E/rR8GfLy7XmtjTPKz/s6q2rmqllXVMprf+c+oqnNwzHtvknHvJtPN9kUA822jucp4Fc07rxtpL2qgWRP1beBb7fb8gWP2BL7Uln8TeOZA2X7A12jWWH0N2He2n6PbtI37ZGWOew+2Scb9Q8D3aN6ofQDYfCpj67j3Yxt13IH9af7adjnNf96XAJ8YON8T2t8FVwL/ycDdKdzmxrYhP+tDx6+kveDRMe/PtoG/42c8043dxkaSJEnSHOeyGUmSJKknDO+SJElSTxjeJUmSpJ4wvEuSJEk9YXiXJEmSesLwLkmSJPWE4V2S5qkkpyc5ZWjfwUluTfKw2eqXJGnDGd4laf56I/CsJGOf9rc58EHgT6vqR9PVSPsx35KkDhjeJWmeqqpbgT8ETkqyJXAUcHVVfTjJoiRHJrk6yeokH02yDTQf+Z3krCQ3JlmT5EtJHjF23iSnJXlfks8luQtYMStPUJIWIMO7JM1jVXUm8N/AGcDvAb/fFv0J8GzgScBS4C6ajwIf82lgb+ChwHeAjwyd+qXAW4GtgQtnqPuSpCGpqtnugyRpBiXZCbga+IuqOr7ddyXwuqo6r/16V+AqYIuqum/o+O2BW4CtququJKcBP6uq3+nyeUiSYJPZ7oAkaWZV1U1JVgOXDezeDfhUksGgXsCOSW4B/hZ4IbA9MFZne5oZeoDrZrbXkqTxuGxGkhamVcDTqmrJwLZ5Vd0IvBJ4FnAI8EvAw9tjMnC8f7aVpFlgeJekhelE4B1JdgNIsmOS57ZlWwP3ALcCDwbePjtdlCQNM7xL0sL098DngP+X5CfAV4Dlbdk/ATe022VtmSRpDvCCVUmSJKknnHmXJEmSesLwLkmSJPWE4V2SJEnqCcO7JEmS1BOGd0mSJKknDO+SJElSTxjeJUmSpJ4wvEuSJEk9YXiXJEmSeuJ/AGMBYhzWjf+qAAAAAElFTkSuQmCC\n", "text/plain": [ "