{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Chapter 12: Markov Chain Monte Carlo\n",
" \n",
"This Jupyter notebook is the Python equivalent of the R code in section 12.4 R, pp. 512 - 517, [Introduction to Probability, 1st Edition](https://www.crcpress.com/Introduction-to-Probability/Blitzstein-Hwang/p/book/9781466575578), Blitzstein & Hwang.\n",
"\n",
"----"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Metropolis-Hastings\n",
"\n",
"Here's how to implement the Metropolis-Hastings algorithm for Example 12.1.8, the Normal-Normal model. First, we choose our observed value of $Y$ and decide on values for the constants $\\sigma$, $\\mu$, and $\\tau$:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"y = 3\n",
"sigma = 1\n",
"mu = 0\n",
"tau = 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We also need to choose the standard deviation of the proposals for step 1 of the algorithm, as explained in Example 12.1.8; for this problem, we let $d = 1$. We set the number of iterations to run, and we allocate a NumPy array `theta` of length 104 which we will fill with our simulated draws:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"d = 1\n",
"niter = 10**4\n",
"theta = np.zeros(niter)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now for the main loop. We initialize $\\theta$ to the observed value $y$, then run the algorithm described in Example 12.1.8:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"theta[0] = y\n",
"\n",
"np.random.seed(1134903170)\n",
"\n",
"from scipy.stats import binom, norm\n",
"\n",
"for i in range(1, niter):\n",
" theta_p = theta[i-1] + norm.rvs(loc=mu, scale=2, size=1)[0]\n",
" numer = norm.pdf(y, loc=theta_p, scale=sigma) * norm.pdf(theta_p, loc=mu, scale=tau)\n",
" denom = norm.pdf(y, loc=theta[i-1], scale=sigma) * norm.pdf(theta[i-1], loc=mu, scale=tau)\n",
" r = numer / denom\n",
" flip = binom.rvs(1, np.min([r, 1]), size=1)[0]\n",
" theta[i] = theta_p if flip==1 else theta[i-1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's step through each line inside the loop. The proposed value of $\\theta$ is `theta_p`, which equals the previous value of $\\theta$ plus a Normal random variable with mean 0 and standard deviation $d$ (recall that [`scipy.stats.norm.rvs`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html) function takes as parameter `scale` the standard deviation and _not_ the variance). The ratio `r` is\n",
"\n",
"\\begin{align}\n",
" \\frac{f_{\\theta|Y}(x^{\\prime}|y)}{f_{\\theta|Y}(x|y)} &= \\frac{e^{-\\frac{1}{2 \\, \\sigma^2}(y-x^{\\prime})^2} \\,\\, e^{-\\frac{1}{2 \\, \\tau^2}(x^{\\prime}-\\mu)^2}}{e^{-\\frac{1}{2 \\, \\sigma^2}(y-x)^2} \\,\\, e^{-\\frac{1}{2 \\, \\tau^2}(x-\\mu)^2}}\n",
"\\end{align}\n",
"\n",
"where `theta_p` is playing the role of $x^{\\prime}$ and `theta[i-1]` is playing the role of $x$. The coin flip to determine whether to accept or reject the proposal is `flip`, which is a coin flip with probability `numpy.min([r, 1])` of Heads (encoding Heads as 1 and Tails as 0). Finally, we set `theta[i]` equal to the proposed value if the coin flip lands Heads, or keep it at the previous value otherwise.\n",
"\n",
"The array `theta` now contains all of our simulation draws. We typically discard some of the initial draws to give the chain some time to approach the stationary distribution. The following line of code discards the first half of the draws:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"theta = theta[-int(niter/2):]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To see what the remaining draws look like, we can create a histogram using [`matplotlib.axes.Axes.hist`](https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.hist.html) function. We can also compute summary statistics such as `numpy.mean(theta)` and `numpy.var(theta)`, which give us the sample mean and sample variance."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAFRCAYAAACVJc8eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHuRJREFUeJzt3XuYZHV95/H3JwwIiAjIaBBGxwtLNCYCGQ1Zk2gAo4CC+sgGYyK6JOyuJN42q6PrRjZrNvg8iRfW9YKiAeIFRaIoGoMiMeRRYECiXFQQRmcYhMkCIhcdwe/+cU5r0VPdXd1d1dXV5/16nnr61Klz+f5OndOfOpc6lapCkiR1wy+MuwBJkrR0DH5JkjrE4JckqUMMfkmSOsTglySpQwx+SZI6xOCXJKlDDH5JkjrE4BcASTYmObzn+dVJnjHGkuZkjZMhyQFJvpbkh0leMe56BjXq9y7J3yZ587DnNX1a07ftYU9/KSxmHUpycpKTR1TaRDL4h6DdsLYl2Xta/yuTVJK1A05jaBvnYlXVL1fVRfMZp18bkrw0ycWLrafftBdS4zAN0t5Ba1xu7/+QvRa4qKoeUlWnLsUMh7E8l3L9GuZ6Msy6l9F2N+c6lOSgJP+S5J4klyZ5VJ9h9mw/PDx9Wv+zkpybJIspMsnfJbk5yZ1Jvp3kjxYzvVEx+IfnRuBFU0+S/Aqwy7AmnmTVsKalblkG686jgavHXMPAFru8xrW8l8H7PEqzrkNJ9gM+C7wFeBhwA/DG6cNV1e3AacCre8b9H8ATgT+oxd/D/q+AtVW1O3A08OYkv7bIaQ6dwT88ZwEv6Xl+PHBm7wBJHpnkE0m2Jrlx6pBVkrOARwGfTnJXkte2/TcmeV2SrwN3J1mV5AlJLkpyR3vI7eie6W9M8vok1yS5PckHk+zc8/qM407X+0m/reGm9pPyt5IcttCFlGR9ku+007omyfOnvb7dvOZYPof3dP9Zkq8n+UGSs6e1/eCeQ4Ufb19/80zzXGj7+rR3+imU+bRvtvd6xvb0zHf6ujPjsm+H/2/t8rs7yelJHpHkc+3wX0iy5wxtnK3OC4HfAd7Ztu3fzbCM+q63c62z81yefbe/WZZX7/o1Vx3bjd+nnQcluaKt9Wxg52njz7q99WvXXHX3eMoMy7eSPL6njt7TD4Nsd3P9P5pxm+yzfPpOa5B1CPgb4H1VdV5V3Qt8FHjKDLN6K/CsJI9LcixwIvDcqrpnptoGVVVXV9WPp562j8ctdrpDV1U+FvkANgKHA98CngDsAGyi+ZRawFqaD1mXA38O7AQ8luZT6bN6p9FnulcCa2iOHuwIXA+8oZ3GocAPgQN6hr+qHX4v4F+AN7evDTLu4X3adEDblke2/dcCj5ttOUzr91Lg4p7nxwKPbJfH7wF3A/u0r804r1mWz+E93Ze2094LuBb4z+1rOwHfBV7ZLocXANuAN8/VPuBdwLsW0d7eGgdu32zv12ztmWndGWDZbwS+CjwC2Be4FbgCOAh4EHAh8KY+y2DW9aod5iLgj+bYfrZbb+ea9jyX5yDb3/TltZFmGxikjduNP62NU+/Zq9vpvRD4CT/fPqfmNdf6OL1dM9Y91/JtXyvg8T3D/i3br0d9t7sB3p+NzLBNznc9YpZ1CNgd+BHwqJ5+LwC+0nafDJw8bZwPAP8AbAUOnmG6nwHumOHxmVnW53cB97TL9gpgt5mGHdfDPf7hmtrrfybwTeCmnteeAqyuqr+oqm1VdQPwPuC4OaZ5alVtquZT7CHAbsAp7TQupFk5X9Qz/Dvb4W8D/rLntUHG7ed+mn/8T0yyY1VtrKrvzDL8J9tP7HckuYNmI/iZqvp4VW2pqp9W1dnAdcBTFziv6U5tp30b8GngwJ62r2pf/0lVnUvzD2nOeVbVy6vq5Qtt7zTzad9s79ds7Zm+PKbWnbmWPcD/qapbquom4J+BS6rqa9Xswfw9zYeA+dQ5H/3W27mmPZ/lOcj294DltYA2zjT+1DR2BN7evmfnAJf1GW4h28Bs850y0/+FxRhkucy0TS5kWjM5jGbZfr1nO/wQzQetmbwVeBZwUlVd0W+AqnpOVe0xw+M5M024/X/xEOC3gHOBH8807LgY/MN1FvD7NHt9Z0577dHAI6eFxBto9rBms6mn+5HApqr6aU+/79LsofUb/rvtOIOOu52quh54Fc2n5luTfLQ9ZPri9rDbXUk+1zPK83o3EOABoZnkJWkuepxaBk8C9p5tXrPVN833e7rvoflHMtX2m6qq9/zdpiHNc9b29prnvGZ7v2ZszzQP6Dfbsm/d0tN9b5/nu7G9Ba1Xc9Q6td7OOu15Ls9Btr9+y5C56hhg/KlpTH/PtgumBa6Ps8233zC9/xcWY5DlMtM2uZBpzWQtcN607fBLNHv0M9mJJpDPHWD681ZV91fVxcB+wH8ZxTwWw+Afoqr6Ls1Ffkey/Qq1Cbhx2qfGh1TVkVOjzzTZnu4twJokve/bo3jgkYU1017bMo9xZ2rXh6vqN/n5qYu3VNWHqmq39nHEXNMASPJomr2sPwEe1m6gVwE/u5K237ymXhpkHjO4Gdg3ecAVuz9bTrPMc+jm0b7Z3q9Z29M7u6mOQZb9Ai14vZqm33o757TnsTzn2v76jTNl0DbOto72e8+2u+oc5lwf+81jkG1jpv8L9wC79rz2i/OY9rDe+8VO60E07QAgyWOAdcB5s4zzZOCqqrpvpgHSXN9y1wyPz8003jSrWIbn+A3+4TsBOLSq7p7W/1LgzvZCnF2S7JDkSUmmLkC5hea842wuoTkv+9okO6b5Lu1zaS5kmXJSkv2S7EWzR3P2PMbdTprvzx6a5EE059HupTkcuRAPpvlHsrWd9sto9joHmdcgy2cmX2mn8yftxU/H0B7iHnL7ZjXP9s32fs3YnlnMuuwXYUHrVR/91ttZpz3P5TnX9jfqNn4FuA94RfuevYA+79kA6+NCt4OZ/i9cCfx+uzyeDTx92nizzW9Y7/1ip3UZ8PT2SOQa4MPAf29PL8zkQJq2z6iqjujZuZn+2G5nJ8nDkxyXZLd2eT6L5lTFhQO0YUkZ/ENWVd+pqg19+t9PsyIfSHNU4N+A9wMPbQf5K+CN7WHIP5th2ttoviJyRDv+u4CXVNU3ewb7MPCPNBcu3UBzkdSg4/bzIOCUdpzvAw+n+ccxb1V1Dc3Vt1+h+YfyKzQXGg0yrzmXzyzz3UZzsc8JNBfm/AHN+cMfzzFPkrwnyXvm1dCZDdy+2d6vOdrT1wDLfkEWsV5Nt916O8C057M859r+RtrGnvfspcDtNBdX9jvMPNf2ttDtoO//BZoLRJ9Lsx69GPjktPFmnN8Q3/vFTutCmusHvg1cDJxVVe+bY5wnM0fwL0DRHNbfTPMe/zXwqqr61JDns2h54CknTbIkG2mufP3CuGtZ7pJcArynqj447lqGYZLb43qrUUp7176qOnm8lSwf7vGrE5I8PckvtodZjwd+ldkv/lnWVlp7JC2dlXynJ6nXAcDHaK4q/g7wwqq6ebwlLcpKa480KheNu4DlxkP9kiR1iIf6JUnqEINfkqQOWZHn+Pfee+9au3btuMuQJGnJXH755f9WVavnGm5FBv/atWvZsGG7r9JLkrRiJZnt9wl+xkP9kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdciKvFe/pO2tXX/+ksxn4ylHLcl8JC2Me/ySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKH+CM90jKwVD+gI0nu8UuS1CEGvyRJHWLwS5LUISML/iQfSHJrkqt6+u2V5IIk17V/92z7J8mpSa5P8vUkB/eMc3w7/HVJjh9VvZIkdcEo9/j/Fnj2tH7rgS9W1f7AF9vnAEcA+7ePE4F3Q/NBAXgT8OvAU4E3TX1YkCRJ8zey4K+qLwO3Tet9DHBG230G8Lye/mdW46vAHkn2AZ4FXFBVt1XV7cAFbP9hQpIkDWipz/E/oqpuBmj/Prztvy+wqWe4zW2/mfpLkqQFWC4X96VPv5ql//YTSE5MsiHJhq1btw61OEmSVoqlDv5b2kP4tH9vbftvBtb0DLcfsGWW/tupqtOqal1VrVu9evXQC5ckaSVY6uA/D5i6Mv944FM9/V/SXt1/CPCD9lTA54HfTbJne1Hf77b9JEnSAozslr1JPgI8A9g7yWaaq/NPAT6W5ATge8Cx7eCfBY4ErgfuAV4GUFW3JflfwGXtcH9RVdMvGJQkSQMaWfBX1YtmeOmwPsMWcNIM0/kA8IEhliZJUmctl4v7JEnSEjD4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqkLEEf5JXJ7k6yVVJPpJk5ySPSXJJkuuSnJ1kp3bYB7XPr29fXzuOmiVJWgmWPPiT7Au8AlhXVU8CdgCOA94CvK2q9gduB05oRzkBuL2qHg+8rR1OkiQtwLgO9a8CdkmyCtgVuBk4FDinff0M4Hlt9zHtc9rXD0uSJaxVkqQVY8mDv6puAv4a+B5N4P8AuBy4o6ruawfbDOzbdu8LbGrHva8d/mFLWbMkSSvFOA7170mzF/8Y4JHAg4Ej+gxaU6PM8lrvdE9MsiHJhq1btw6rXEmSVpRxHOo/HLixqrZW1U+Ac4F/D+zRHvoH2A/Y0nZvBtYAtK8/FLht+kSr6rSqWldV61avXj3qNkiSNJHGEfzfAw5Jsmt7rv4w4BrgS8AL22GOBz7Vdp/XPqd9/cKq2m6PX5IkzW0c5/gvoblI7wrgG20NpwGvA16T5Hqac/int6OcDjys7f8aYP1S1yxJ0kqxau5Bhq+q3gS8aVrvG4Cn9hn2R8CxS1GXJEkrnXfukySpQwx+SZI6xOCXJKlDDH5JkjrE4JckqUPGclW/pJVr7frzRz6PjaccNfJ5SCuVwS/NYilCTJKWkof6JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDhko+JM8adSFSJKk0Rt0j/89SS5N8vIke4y0IkmSNDIDBX9V/SbwYmANsCHJh5M8c6SVSZKkoRv4HH9VXQe8EXgd8HTg1CTfTPKCURUnSZKGa9Bz/L+a5G3AtcChwHOr6glt99tGWJ8kSRqiVQMO907gfcAbqureqZ5VtSXJG0dSmSRJGrpBg/9I4N6quh8gyS8AO1fVPVV11siqkyRJQzXoOf4vALv0PN+17SdJkibIoMG/c1XdNfWk7d51NCVJkqRRGTT4705y8NSTJL8G3DvL8JIkaRka9Bz/q4CPJ9nSPt8H+L3RlCRJkkZloOCvqsuS/BJwABDgm1X1k5FWJkmShm7QPX6ApwBr23EOSkJVnTmSqiRJ0kgMFPxJzgIeB1wJ3N/2LsDglyRpggy6x78OeGJV1TBm2v7Qz/uBJ9F8gPiPwLeAs2mOKmwE/kNV3Z4kwDto7iVwD/DSqrpiGHVIktQ1g17VfxXwi0Oc7zuAf6iqXwKeTHMr4PXAF6tqf+CL7XOAI4D928eJwLuHWIckSZ0y6B7/3sA1SS4FfjzVs6qOnu8Mk+wO/Dbw0nYa24BtSY4BntEOdgZwEc0PAh0DnNkebfhqkj2S7FNVN8933pIkdd2gwX/yEOf5WGAr8MEkTwYuB14JPGIqzKvq5iQPb4ffF9jUM/7mtp/BL0nSPA10qL+q/onmvPuObfdlwELPs68CDgbeXVUHAXfz88P6/aRfSdsNlJyYZEOSDVu3bl1gaZIkrWyD/izvHwPnAO9te+0LfHKB89wMbK6qS9rn59B8ELglyT7t/PYBbu0Zfk3P+PsBW5imqk6rqnVVtW716tULLE2SpJVt0Iv7TgKeBtwJUFXXAQ+fdYwZVNX3gU1JDmh7HQZcA5wHHN/2Ox74VNt9HvCSNA4BfuD5fUmSFmbQc/w/rqptzTfrIMkq+hxun4c/BT6UZCfgBuBlNB9CPpbkBOB7wLHtsJ+l+Srf9TRf53vZIuYrSVKnDRr8/5TkDcAuSZ4JvBz49EJnWlVX0twbYLrD+gxbNEccJEnSIg16qH89zZX43wD+E81e+BtHVZQkSRqNQX+k56fA+9qHJEmaUIPeq/9G+pzTr6rHDr0iSZI0MvO5V/+UnWkuvNtr+OVIkqRRGvQGPv+v53FTVb0dOHTEtUmSpCEb9FD/wT1Pf4HmCMBDRlKRJEkamUEP9f9NT/d9tD+bO/RqJEnSSA16Vf/vjLoQSZI0eoMe6n/NbK9X1VuHU44kSRql+VzV/xSa++YDPBf4Mg/8uVxJkrTMDRr8ewMHV9UPAZKcDHy8qv5oVIVJkqThG/SWvY8CtvU83wasHXo1kiRppAbd4z8LuDTJ39Pcwe/5wJkjq0qSJI3EoFf1/2WSzwG/1fZ6WVV9bXRlSZKkURj0UD/ArsCdVfUOYHOSx4yoJkmSNCIDBX+SNwGvA17f9toR+LtRFSVJkkZj0D3+5wNHA3cDVNUWvGWvJEkTZ9Dg31ZVRfvTvEkePLqSJEnSqAwa/B9L8l5gjyR/DHwBeN/oypIkSaMw6FX9f53kmcCdwAHAn1fVBSOtTJIkDd2cwZ9kB+DzVXU4YNhLkjTB5jzUX1X3A/ckeegS1CNJkkZo0Dv3/Qj4RpILaK/sB6iqV4ykKkmSNBKDBv/57UOSJE2wWYM/yaOq6ntVdcZSFSRJkkZnrnP8n5zqSPKJEdciSZJGbK7gT0/3Y0dZiCRJGr25gr9m6JYkSRNorov7npzkTpo9/13abtrnVVW7j7Q6SZI0VLMGf1XtsFSFSNKg1q4f/ZeMNp5y1MjnIY3DoPfqlyRJK4DBL0lShxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1yNiCP8kOSb6W5DPt88ckuSTJdUnOTrJT2/9B7fPr29fXjqtmSZIm3Vz36h+lVwLXAlP3+38L8Laq+miS9wAnAO9u/95eVY9Pclw73O+No2AtL0tx21ZJWmnGssefZD/gKOD97fMAhwLntIOcATyv7T6mfU77+mHt8JIkaZ7Gdaj/7cBrgZ+2zx8G3FFV97XPNwP7tt37ApsA2td/0A4vSZLmacmDP8lzgFur6vLe3n0GrQFe653uiUk2JNmwdevWIVQqSdLKM449/qcBRyfZCHyU5hD/24E9kkxdc7AfsKXt3gysAWhffyhw2/SJVtVpVbWuqtatXr16tC2QJGlCLXnwV9Xrq2q/qloLHAdcWFUvBr4EvLAd7HjgU233ee1z2tcvrKrt9vglSdLcltP3+F8HvCbJ9TTn8E9v+58OPKzt/xpg/ZjqkyRp4o3z63xU1UXARW33DcBT+wzzI+DYJS1MkqQVajnt8UuSpBEz+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6pAlD/4ka5J8Kcm1Sa5O8sq2/15JLkhyXft3z7Z/kpya5PokX09y8FLXLEnSSjGOPf77gP9aVU8ADgFOSvJEYD3wxaraH/hi+xzgCGD/9nEi8O6lL1mSpJVh1VLPsKpuBm5uu3+Y5FpgX+AY4BntYGcAFwGva/ufWVUFfDXJHkn2aacjSSOxdv35SzKfjacctSTzkaaM9Rx/krXAQcAlwCOmwrz9+/B2sH2BTT2jbW77SZKkeRpb8CfZDfgE8KqqunO2Qfv0qz7TOzHJhiQbtm7dOqwyJUlaUcYS/El2pAn9D1XVuW3vW5Ls076+D3Br238zsKZn9P2ALdOnWVWnVdW6qlq3evXq0RUvSdIEG8dV/QFOB66tqrf2vHQecHzbfTzwqZ7+L2mv7j8E+IHn9yVJWpglv7gPeBrwh8A3klzZ9nsDcArwsSQnAN8Djm1f+yxwJHA9cA/wsqUtV5KklWMcV/VfTP/z9gCH9Rm+gJNGWpQkSR3hnfskSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDhnHvfrVAWvXnz/uEiRJfbjHL0lShxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKHrBp3AZLUZWvXnz/yeWw85aiRz0OTwz1+SZI6xOCXJKlDDH5JkjrE4JckqUMMfkmSOsSr+jtmKa4gliQtX+7xS5LUIQa/JEkdYvBLktQhBr8kSR3ixX2StMJ5W2D1co9fkqQOmZjgT/LsJN9Kcn2S9eOuR5KkSTQRwZ9kB+D/AkcATwRelOSJ461KkqTJMynn+J8KXF9VNwAk+ShwDHDNUhXgOTJJ0kowKcG/L7Cp5/lm4NfHVMvIeFc9SZNqqf5/LcUO0krf0ZuU4E+ffvWAAZITgRPbp3cl+dbIq9re3sC/jWG+o2BblqeV1BZYWe2xLUsgb5n3KMuyLQtoB8zdlkcPMpFJCf7NwJqe5/sBW3oHqKrTgNOWsqjpkmyoqnXjrGFYbMvytJLaAiurPbZlebIt25uIi/uAy4D9kzwmyU7AccB5Y65JkqSJMxF7/FV1X5I/AT4P7AB8oKquHnNZkiRNnIkIfoCq+izw2XHXMYexnmoYMtuyPK2ktsDKao9tWZ5syzSpqrmHkiRJK8KknOOXJElDYPAPWZJjk1yd5KdJJvJK0pVye+QkH0hya5Krxl3LYiVZk+RLSa5t169XjrumhUqyc5JLk/xr25b/Oe6aFivJDkm+luQz465lMZJsTPKNJFcm2TDuehYjyR5JzknyzXa7+Y1x17QQSQ5o34+px51JXrWoaXqof7iSPAH4KfBe4M+qaqI2nvb2yN8GnknzNcrLgBdV1ZLdJXFYkvw2cBdwZlU9adz1LEaSfYB9quqKJA8BLgeeN6HvS4AHV9VdSXYELgZeWVVfHXNpC5bkNcA6YPeqes6461moJBuBdVW17L73Pl9JzgD+uare334bbNequmPcdS1G+//5JuDXq+q7C52Oe/xDVlXXVtU4bh40LD+7PXJVbQOmbo88carqy8Bt465jGKrq5qq6ou3+IXAtzR0tJ0417mqf7tg+JnYPJMl+wFHA+8ddixpJdgd+GzgdoKq2TXrotw4DvrOY0AeDX9vrd3vkiQyYlSrJWuAg4JLxVrJw7aHxK4FbgQuqamLbArwdeC3Nkb5JV8A/Jrm8vRvqpHossBX4YHsK5v1JHjzuoobgOOAji52Iwb8ASb6Q5Ko+j4ncM55mztsja3yS7AZ8AnhVVd057noWqqrur6oDae7C+dQkE3kqJslzgFur6vJx1zIkT6uqg2l+CfWk9nTZJFoFHAy8u6oOAu4GJvZ6JYD2dMXRwMcXO62J+R7/clJVh4+7hhGa8/bIGo/2fPgngA9V1bnjrmcYquqOJBcBzwYm8SLMpwFHJzkS2BnYPcnfVdUfjLmuBamqLe3fW5P8Pc2pvy+Pt6oF2Qxs7jmSdA4THvw0H8auqKpbFjsh9/g1nbdHXobaC+JOB66tqreOu57FSLI6yR5t9y7A4cA3x1vVwlTV66tqv6paS7OtXDipoZ/kwe2Fo7SHxX+XyfwwRlV9H9iU5IC212Es4c+4j8iLGMJhfjD4hy7J85NsBn4DOD/J58dd03xU1X3A1O2RrwU+Nqm3R07yEeArwAFJNic5Ydw1LcLTgD8EDu35Ws+R4y5qgfYBvpTk6zQfNC+oqon+GtwK8Qjg4iT/ClwKnF9V/zDmmhbjT4EPtevZgcD/HnM9C5ZkV5pvWg3lSJ9f55MkqUPc45ckqUMMfkmSOsTglySpQwx+SZI6xOCXJKlDDH5JkjrE4JckqUMMfklD0f7wzjuSXN3+pvtjx12TpO0Z/JKG5fXADVX1y8CpwMvHXI+kPvyRHkmL1t7b/flV9WttrxtpfqNe0jJj8EsahsOBNUmubJ/vBXxhjPVImoGH+iUNw4HAn1fVgVV1IPCPwJVzjCNpDAx+ScOwJ3APQJJVND/p+umxViSpL4Nf0jB8Gzik7X41zU+63jjGeiTNwJ/llbRoSfYEPgfsDXwFOLGq7h1vVZL6MfglSeoQD/VLktQhBr8kSR1i8EuS1CEGvyRJHWLwS5LUIQa/JEkdYvBLktQhBr8kSR3y/wFw7OOB5yBZ7wAAAABJRU5ErkJggg==\n",
"text/plain": [
"