{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2016-03-02T07:49:55.087088", "start_time": "2016-03-02T07:49:55.073850" } }, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "from scipy.stats import norm,uniform\n", "import matplotlib.pyplot as plt\n", "import seaborn as sbn\n", "from scipy.optimize import minimize\n", "import pandas as pd\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", " \n", "sbn.set_style('white')\n", "sbn.set_context('talk')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this workbook, I will briefly go over the theory of MCMC. We'll cover questions like why does MCMC work (why does it allow us to draw from the posterior)? and what conditions are necessary to get this result.\n", "\n", "## Asymptotics\n", "\n", "A decent starting point is to go back to our simple example using conjugate priors and analytical Bayes. Our data is distributed normal with mean $\\mu$ (unknown) and standard deviation equal to $\\sigma$ (known). We have priors describing our beliefs about the distribution of $\\mu$: distributed normal with mean $\\mu_0$ and standard deviation $\\sigma^2_0$. In this setting, our posterior can be solved for analytically as: \n", "\n", "$$\n", "N \\left (\\frac{\\sigma_0^2}{\\frac{\\sigma^2}{n} + \\sigma^2_0}\\frac{\\sum y_i}{n} + \\frac{\\frac{\\sigma^2}{n}}{\\frac{\\sigma^2}{n} + \\sigma^2_0} \\mu_0,\\left ( {\\frac{1}{\\sigma_0^2} + \\frac{n}{\\sigma^2}} \\right )^{-1} \\right ) \n", "$$\n", "\n", "How might we expect the posterior to behave as sample size get increasingly large? First examine the mean of the analytical posterior: \n", "\n", "$$\n", "\\lim_{n\\to \\infty} \\frac{\\sigma_0^2}{\\frac{\\sigma^2}{n} + \\sigma^2_0}\\frac{\\sum y_i}{n} + \\frac{\\frac{\\sigma^2}{n}}{\\frac{\\sigma^2}{n} + \\sigma^2_0} \\mu_0 \\to \\frac{\\sum y_i}{n}\n", "$$\n", "which means that the mean of the analytical posterior will converge on the Maximum Likelihood Estimate of the mean: the data completely dominates the prior.\n", "\n", "What about the variance? Doing a similar excercise shows that\n", "\n", "$$\n", "\\lim_{n\\to \\infty} \\left ( {\\frac{1}{\\sigma_0^2} + \\frac{n}{\\sigma^2}} \\right )^{-1} = \\lim_{n\\to \\infty} \\left ( \\frac{\\sigma^2_0 \\sigma^2}{(\\sigma^2 + n \\sigma^2_0)}\\right ) \\to 0\n", "$$\n", "\n", "So, the asymptotic properties of our Bayesian estimator, exactly parallels the frequentist worldview: our posterior collapses to the Maximum Likelihood Mean with zero variance for super huge sample sizes. \n", "\n", "We can also see this using numerical methods and revisit our example from 2 worksheets ago:\n", "\n", "1. Mean of data is unknown, variance is known ($\\sigma=3$)\n", "2. Priors $\\mu_0 = 8$, $\\sigma_0 = 2$\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2016-03-02T07:49:59.646410", "start_time": "2016-03-02T07:49:59.641277" } }, "outputs": [], "source": [ "# plot the MLE mean, and MAP estimate at different sample \n", "# sizes of data given below\n", "sample_size = [1,10,10**2,10**3,10**6]\n", "\n", "sigma = 3\n", "mu_0 = 8\n", "sigma_0 = 2" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "ExecuteTime": { "end_time": "2016-03-02T07:50:00.287112", "start_time": "2016-03-02T07:50:00.280967" } }, "outputs": [], "source": [ "def log_posterior_numerical(mu,data,sigma,mu_0,sigma_0):\n", " log_prior = norm(mu_0,sigma_0).logpdf(mu)\n", " log_like = norm(mu,sigma).logpdf(data).sum()\n", " return -1*(log_prior + log_like)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "ExecuteTime": { "end_time": "2016-03-02T07:50:08.910751", "start_time": "2016-03-02T07:50:00.945060" } }, "outputs": [], "source": [ "# for each sample size, generate a dataset and calculate\n", "# 1. MAP\n", "# 2. MLE (just the mean of the data)\n", "# and store results\n", "\n", "results = np.zeros((len(sample_size),4))\n", "\n", "count=0\n", "for i in sample_size:\n", " data = norm(10,3).rvs(i)\n", " # for this dataset, solve for map\n", " res = minimize(log_posterior_numerical,np.mean(data), \n", " args=(data,sigma,mu_0,sigma_0),options={'disp': False},\n", " method='BFGS')\n", " \n", " results[count,:] = [i,mu_0,res.x,np.mean(data)]\n", " count+=1" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "ExecuteTime": { "end_time": "2016-03-02T07:50:08.918755", "start_time": "2016-03-02T07:50:08.912329" } }, "outputs": [], "source": [ "# convert to a pandas dataframe\n", "results=pd.DataFrame(results,columns=['Sample Size','Prior','MAP','MLE'])" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "ExecuteTime": { "end_time": "2016-03-02T07:50:08.970270", "start_time": "2016-03-02T07:50:08.921221" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Sample SizePriorMAPMLE
01.08.09.01941311.313092
110.08.08.0713628.087419
2100.08.09.99967910.044672
31000.08.010.07535910.080028
41000000.08.09.9990809.999085
\n", "
" ], "text/plain": [ " Sample Size Prior MAP MLE\n", "0 1.0 8.0 9.019413 11.313092\n", "1 10.0 8.0 8.071362 8.087419\n", "2 100.0 8.0 9.999679 10.044672\n", "3 1000.0 8.0 10.075359 10.080028\n", "4 1000000.0 8.0 9.999080 9.999085" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results.head()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "ExecuteTime": { "end_time": "2016-03-02T07:50:09.371755", "start_time": "2016-03-02T07:50:08.973372" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAE/CAYAAABW/Dj8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAG4xJREFUeJzt3XmYXFWd//F3J+kkLMoiSwgIOAhf1mgQREFFIoIGNRJhFBhFRbYRZFzBhRFcQBGYEUQHERAkLLJKiAsKI0KcARHkh4BfFgGBIYJGlrAkIenfH7fSXV3dnaRSla6uW+/X8/QDfe+tW99cwqdOnXvOPV09PT1IktrfqFYXIElqDgNdkkrCQJekkjDQJakkDHRJKokxrXrjiBgH7Ag8DixqVR2S1EZGAxsAv8vM+bU7WxboFGF+YwvfX5La1ZuBm2o3tjLQHweYMWMGEyZMaGEZktQe5syZwwEHHACV/KzVykBfBDBhwgQ22mijFpYhSW1n0G5qb4pKUkkY6JJUEga6JJWEgS5JJWGgS1JJGOiSVBIGuiSVRCvHoY9o25233XIfe+eBd67ESiRp+Rjoklpq02NmDcv7PPSNvYblfVrJQJekQUQE48ePZ9Soome6q6uLyZMnc/TRR7PFFlsM+prJkydz2WWXsdlmmw1nqb0MdEkawqWXXtob3gsXLuTUU0/l4IMP5vrrr2f06NEDjr/99tuHu8R+vCkqScuhu7ub6dOnM2fOHJ5++mmuuOIK9t9/f/bdd1922mknHn74YSKCe++9F4DZs2czffp0tt9+e6ZNm8YNN9zQe66I4Pjjj2fHHXfkzDPPbFqNttAlaTk8/fTT/OhHP2KLLbZg7bXXBuC2227j3HPPZdttt+VlL3tZ77H33Xcfhx9+OCeffDJTpkxh9uzZHHXUUVxyySVEBADz589n9uzZLFiwoGk1GuiSNIQPfOADvX3oY8eOZdKkSZx22mm9+9ddd13e+MY3DnjdrFmz2Hnnndljjz0A2HXXXZkyZQozZ87sDfS99tqLsWPHMnbs2KbVa6BL0hAuvvjiIW+AQhHog5k7dy4TJ07st23ixInMmTOn9/d11lmnOUVWsQ9dkppsgw024LHHHuu37dFHH+0X4l1dXU1/XwNdkpps6tSp3HzzzVx77bUsWrSIG264geuvv56pU6eu1Pe1y0VSS5Vxws8mm2zCGWecwcknn8zRRx/NhhtuyCmnnMKkSZNW6vsa6JI0iMxc6v7p06czffr0IV+zyy67sMsuu6zQuVeUXS6SVBIGuiSVhIEuSSVhoEtSSRjoklQSBroklYSBLkklYaBLUknUNbEoIl4PXJWZEyu/rwWcA0wBngaOz8yzm16lpPI6bo1hep+nh+d9Wmi5WugR0RURHwWuBaqf9XgWMA9YH9gHOCkiVu7cVkkaBhHBa17zGubNm9dv+8KFC9lpp52YMmXKgNfsv//+vOENb2D+/Pn9tp9++ulsvfXWTJ48ufdnypQpnHHGGU2teXm7XL4AHAV8fcmGiFgdeC/w5cx8MTNvAS4EDm5qhZLUIuPHj+e6667rt+3GG29k4cKFA4594IEHmDNnDltttRUzZ84csH/33Xfn9ttv7/0566yzmDFjBhdffHHT6l3eLpdzgBOAXau2bQ4szMw/V21LoP/DDaQS2e687Zb72DsPvHMlVrIS1Nn1sd2rNl7uY9vuWlTsueeezJo1i2nTpvVumzlzJnvssQe33HJLv2MvueQS3va2tzFp0iTOOecc9tlnn6Wee7PNNmOHHXboXbKuGZYr0DPzcaB3pY2K1YAXag59Hli1KZVJw6He/ts6Qkxt7Lg1gAlMffJ7HHrTWvzj82uy1rge5i3s4nf/vQ7Hvu4Zbnnq5b1/f7Z95StZ+OOFjDlgDBc9cRELH1zINl/bhlGvLDpBXvrDS/Q82cN2523HnQfeyaJFi7jjjju4+eab+cpXvtK0sht52uLzwCo121al6FOXpLa39rjF7LjeAq59ZDzvf/UL/PKRcbx14nzG1nRWL87FdK3Vxaj1ih2jXjOKxb9f3BvoAD339bDglAXscPoO9PT0MGHCBA455BD23HPPptXbSKDfB3RHxMaZ+ZfKtgDubrwsSRoZ3rXJC1z+51V4/6tfYObDq3DYNvN4bmH/RF98+2J6nuxhwbcrCz4vAhZAz7weulYvVibq2ryL7vd1c+uBt660Wld4HHpmPgv8BDgxIlaNiB2B/YEZzSpOklrt7RvN549zu7lr7hj+Mm80O67b/4bog8+MpuexHroP6u77OaSbrg26WHTbomGttdGJRQcD3cCjwOXAZzPz5oarkqQRYrXuHt46cT6f+581mLrxi9QuBfrjB1al65+66Fq7i67V+35GTRpVtNwX9QxbrXV1uWTmr4F1qn6fC/xzk2uS1EnaYMLPuzd9kcN/sxbf3vSpftsXLIKrHhzPqHcMbBuP2moUi365iMV/WjxcZboEnSQNJveb0/vvUzac3+/33Tacz24bPgnA/0x/ctAhnF3juxj7ubEDtq9MPstFkkrCQJekkjDQJakkDHRJKgkDXZJKwkCXpJIw0CWpJAx0SSoJJxZJaql6njHfiHqfyR4XTWD86B5m7/0Eq3f3Td9fuBjedOV6rNbdw/XveZJH541mwQkL6P5MN11juwacZ8EZC+A5oLJr8n9OBmDrrbdmxozmPvrKQJekIYwf3cN1j45j2qte7N124+PjWFjnbP4xe49h1OZFh8jtB97ezBL7sctFkoaw58YvMusv4/ttm/nQePZ45fwhXtFandVCr2d1GlemkTre1I1f5NAb1uIf87v6Vix6cmyxYtETw/ucluVhC12ShlC9YhEw5IpFy/LSVS+x4JQFxYpFO+zADjvswAUXXND0ejurhS5JdVqeFYuWZcx7+/rQR+SKRZLUCZa1YtFIYgtdkpZiWSsW9XoWeqqGN9INXasMdfDKYaBL0jIMtWJRtYVn9m+5j9pmFGOmFRH70pUvDRiHDnDrrbcyevToptVpoEtqqXon/AyX5V2xaKPVFzH2C0OPeBn78f77HIcuSVomA12SSsJAl6SSMNAlqSQMdEkqCQNdkkrCQJekkjDQJakkGp5YFBE7A6cBWwCPA8dn5oWNnleSVJ+GWugRMRq4CvhGZr4c+BhwXkRs2oTaJEl1aLSFviawLjAmIrqAxcACYFGjhWklq2Oxj+3qWOxjpE7jljpBQy30zPw78F3gImAhcCNwRGY+0oTaJEl1aKiFHhGjgOeBfYGrgbcDF0bEbZl5RxPqW6pNj5lV1/EPjV/2Me3Ka9GnnmtR5usAXotqnXAtGh3lMh3YKTMvy8wFmTkLmAV8qPHSJEn1aDTQNwbG1WxbWPmRJA2jRm+K/hI4MSI+AvwQeAuwNzClwfNKkurU6E3RO4F9gKOAp4EzgAMzc+WtgipJGlTDE4sycyYwswm1SJIa4NR/SSoJA12SSsJAl6SSMNAlqSQMdEkqCQNdkkrCQJekkjDQJakkDHRJKgkDXZJKwkCXpJIw0CWpJAx0SSoJA12SSsJAl6SSMNAlqSQMdEkqCQNdkkrCQJekkjDQJakkDHRJKgkDXZJKwkCXpJIw0CWpJAx0SSoJA12SSmJMoyeIiI2A/wLeAjwDnJSZpzV6XklSfRpqoUdEF3AVcA/wCmBP4LiI2LkJtUmS6tBoC30nYCJwTGYuAu6KiDcCf2u4MklSXRrtQ98euAs4KSLmRMS9wBsy8++NlyZJqkejgb42sBtFi3xj4MPA6RHx5gbPK0mqU6NdLvOBuZl5YuX330bE5cA04MYGzy1JqkOjLfQEVouI6g+G0UBXg+eVJNWp0Rb6L4F/AN+IiGOA1wN7A29vtDBJUn0aaqFn5gvAW4FtgSeAC4FPZOb/Nl6aJKkeDU8sysz7gXc0oRZJUgOc+i9JJWGgS1JJGOiSVBIGuiSVhIEuSSVhoEtSSRjoklQSBroklYSBLkklYaBLUkkY6JJUEga6JJWEgS5JJWGgS1JJGOiSVBIGuiSVhIEuSSVhoEtSSRjoklQSBroklYSBLkklYaBLUkkY6JJUEga6JJWEgS5JJWGgS1JJGOiSVBJNC/SIWD8inoiIdzXrnJKk5dfMFvrZwCuaeD5JUh2aEugRcRjwHPBIM84nSapfw4EeEZsDnwYOb7wcSdKKaijQI2IMcAFwVGbObU5JkqQV0WgL/VjgD5n502YUI0lacY0G+vuBD0TEUxHxFLAxcHFEHNN4aZKkeoxp5MWZuWX17xHxEHBEZl7TyHklSfVzYpEklURDLfRamblpM88nSVp+ttAlqSQMdEkqCQNdkkrCQJekkjDQJakkDHRJKgkDXZJKwkCXpJIw0CWpJAx0SSoJA12SSsJAl6SSMNAlqSQMdEkqCQNdkkrCQJekkjDQJakkDHRJKgkDXZJKwkCXpJIw0CWpJAx0SSoJA12SSsJAl6SSMNAlqSQMdEkqCQNdkkpiTKMniIg3AacAWwJ/A07KzDMbPa8kqT4NtdAjYi3gauA0YC1gX+DEiNi9CbVJkurQaAt9E2BWZs6o/H5bRPw3sDPwqwbPLUmqQ0OBnpl/AD645PdKi/3NwPkN1iVJqlPTbopGxBrATOD3lX9KkoZRUwI9Il4F/BaYC0zPzMXNOK8kafk1HOgRsT1wM/AL4L2Z+ULDVUmS6tZQH3pErA/8HDglM7/ZnJIkSSui0VEuBwHrAsdGxLFV27+dmV9s8NySpDo0OsrlBOCEJtUiSWqAU/8lqSQMdEkqCQNdkkrCQJekkjDQJakkDHRJKgkDXZJKwkCXpJIw0CWpJAx0SSoJA12SSsJAl6SSMNAlqSQMdEkqCQNdkkrCQJekkjDQJakkDHRJKgkDXZJKwkCXpJIw0CWpJAx0SSoJA12SSsJAl6SSMNAlqSQMdEkqiTGNniAiJgNnAtsA9wGHZeb/NnpeSVJ9GmqhR8R4YCZwLrAmcBpwRUSMbUJtkqQ6NNrlshuwODO/l5kLM/Mc4O/AuxsvTZJUj0a7XLYE7q7ZlhTdL5cv47WjAebMmbPi7/7c3LoOf/Sl0ct9bM9TPct/3kcfrauOlcJr0aeOa1HPdQCvRTWvRZ/huhZVeTlogV09PctfSK2I+BKwfWZOr9p2PvB/mXnMMl77JuDGFX5zSepcb87Mm2o3NtpCfx5YpWbbqsC85Xjt74A3A48DixqsQ5I6wWhgA4r8HKDRQL8HOKJmWwAXLuuFmTkfGPAJI0laqgeG2tFooF8PjIuII4H/Aj4IrA/8osHzSpLq1NAol0or+53AfsBc4EjgPZn5XBNqkyTVoaGbopKkkcOp/5JUEga6JJWEgS5JJWGgS1JJGOiSVBIGuiSVhIEuScvQLo8E78hx6BGxDsWTIpc8d+aezPxHa6tqjcoCJdvQdy3uzsw/tLaq4RcRYyge+9zvWgC/yMwXWlnbcIuIVwAfYuC1uDQzH2tlbcMlIrqAfwUOBbaiaPy+RPG4k/OB/8jMEReeHRXoEbEu8D1gGvAUxcPFVgXWAn4KfDQz/9a6CodPRGwOXAT8E8WzIZZci1cDjwDTM/P+1lU4fCJiJ+AKimf5J33XIoD1gHdn5qAPQyqbiNiL4llMNzDwWuwKfCAzf9a6CodHRPwH8CbgRPpfhy2BY4AbM/NTratwcA0vQddmzgb+Cqyfmb0PR6602E+kWHmpUxbn+AFwNfC1zFy8ZGNEjAKOpbhWu7aotuH2PeBLmXlu7Y6I+CjFc4peN+xVtcYpwPsz8+e1OyLincCpQOkDneIbylaZ+UTN9rsi4ibgj8CIC/RO60N/K3BEdZgDVFrlR1I8zrdTbA+cUB3mAJXfTwBe25KqWmMLiq/Rgzmf4ltLp5gI/HKIfddRPLq1Eyztkd5jgAXDVUg9Oq2F/iSwHXDrIPu2p2i9d4qHgD0ouppq7QU8OKzVtNZdwEHA9wfZdyhFa6xT3AIcHxFfrTx8D4CI6Aa+AnTKAvBnA9dGxH9S3D9YsvZDAJ+t7B9xOi3QvwT8KiKuYuB/pPcBB7ewtuH2SeDyiLid/tdiS4ruhWktrG24HQZcHRFfYODfi3EUTxTtFAcBlwKfioiH6LsWm1Jcm+lDvbBMMvPzEfEYxd+NLYHVKK7FPcBZwHdbWN6QOuqmKEBETAL2Z+B/pIsy845W1jbcImI9iv9Ba6/FlZnZSd9WiIhxFIue116L6zNzRH69XpkiYktqrkVm/qm1VWlZOi7QJalREbFfZl7U6jpqddpN0aWKiBH5NaoVImJWq2sYKSLizlbXMFJExDOtrmGE+GKrCxiMgd5fV6sLGEFc77XPia0uYASZ2uoCRoLM3LbVNQzGLhdJQ4qIsZ14D2GJdptJbaB3sIjYlmJkT+0U7xmdMjNyicoMyUMYeC3Oz8wrWlnbcGrXKe/N1q4zqTtq2GJlGFL30o7JzPWGqZyWqsyAPIlimveV9J/i/dOI+ExmntfCEodNZbjix4DvUIxFr74WJ0fE5pn5zRaWOJxOpZjyfhyDT3nfiBE4Q3IlaMuZ1B0V6MDewC+ArzH45KJO8u/AOwdriUfEjyjGIndEoFOZJTxIi2tWRMwEfgN0SqC35ZT3lWB74G2DzaSOiBMYodegowI9M2+JiE8DR2bmKa2up8XWBIbqC/wj8LJhrKXVxgJzhtj3Nzpr8EBbTnlfCR6iDWdSd1SgA2TmORExISLW77TJMzWuB34QEf+emQ8v2RgRG1GM6vhVyyobflcBV0TE1yn6iqtnih4HXN660oZdW055Xwnacia1N0U7VESsSdFPOI2i1fUCxV/YbuAnwCGd8oz4yuIFJwAHAOtX7XoCuIDiSYwvtqK2VoiII4B/YeCs2QuA79Z2Q5RVRKxP0U3bNjOpDfQOFxGrUTxtcFWKv7D3Zea81lbVOpUPutWA5zvlA01Da7fFcAx0iaUO4bwgMzvqBrpDOAcshvMPim+wI34xHAO9QzmEs0/NEM7aVXr2Azp1CGfttTgCOLMThnBGxNUUj9M+eojFcCZk5ohbDKfjboqql0M4+ziEs49DOAtvBdatfiY8FIvhRMSRDD0qqqUM9A7lEM5+HMLZxyGchbZcDMdA72AO4ezlEM4+DuEstOViOPahq+M5hLOPQzj7tONiOAa6VOEQzv4cwtn7sLINqYz2ycz/a3FJS2WgS+rHIZwQEasAXwU+CqxRtetpim8qnxuJ31TsQ1fHcwhnn2U8hfNnHfQUzjOBVwC7M/Cpk8dTPJXzQy2rbggGuuQQzmoO4Sy8B3hlZj5bte054PcRsT/Fw7tGHANdHc8hnP04hLMwD1gPeHaQfRtRdL2MOAa6hEM4qziEs3Ay8OuIOJuBwxYPZYROrvKmqKReDuHsExHvYoinTmbmNa2sbSgGuqQBHMLZX7sslm2gS1KNqsWyDwG2pk0Wy7YPXVIvh3D2WrJY9vG00WLZBrqkag7hLLTlYtkGuqReDuHs1ZaLZRvokvpxCCfQpotle1NUkgbRjotlG+iSVBJ2uUjSINpxsWxb6JJUo10XyzbQJalGRDzO4ItlExGbA7/JzA2Gv7Kl65QFXyWpHm25WLZ96JI0UFsulj0iP2UkqcUOB/4fcDFFS/0Z4K/AZRTPix9xs0TBPnRJWqp2WizbQJekQbTjYtl2uUhSjcpi2b8GuigWyz6LYoGPURSLZR/YuuqG5k1RSRqoLRfLtoUuSQO15WLZBrokDbRksexNqjdWFsv+PiN0sWy7XCRpoI9SLJZ9f0QMulh2C2sbkqNcJGkI7bZYtoEuSSVhl4sk1WjXxbINdEkaqC0Xy7bLRZIGUZlcdGRmTm51LcvLYYuSNIjMPAe4NCLWb3Uty8sWuiSVhC10SSoJA12SSsJRLmqZiBgDfA74MLAJ8BTFyIJjM/PhFtX0a+DWzPzMCr5+N+DLwA6VTXcAJ2XmTyr7NwUeBLbLzD82XLBUxRa6WukE4CDgkxRLe00DNgBuiIhVW1nYioiI1wI/q/y8rvJzDcVSZu+qHPYIxZ/xTy0pUqVmC12tdBDwb5k5q/L7QxGxL/AEMJViua928kFgdmZ+s2rbiRGxPXAYcE1mLmLoxYelhhjoaqXFwFsi4sJK0JGZT1VWinkcertlvg7sR9GyfRL4YWZ+obL/h8Bcisedvp/iw+BwYCLwFYqlw87NzE9VHb8AWAN4D/Aw8PnMvHKwAiPiI8AXgA2BuyrHDvWkvcXAlhGxQWY+XrX9CCqzDqu7XCi6Zc4d5Dw/zMyPREQ3xbeYDwHjgJuAT2Tmn4d4f3U4u1zUSqcCHwMejojvR8R+EbF2Zt6bmc9WjjkG+Gdgf2Bz4KvA5yPiLVXn+TjFyuyTgNsoFvbdD9gT+CzwyYjYper4D1ME/2TgfIqxxq+tLS4ipgKnAJ+vnPtHwDURsd0Qf54fUDwn+8GIuCYiPhUR22bmXzPz0UGOv4TiQ2rJz2HAi8Dplf1fB3YH3ge8keJD7vqIWGWI91eHM9DVMpl5IkVY3wt8BLgQeDwiToyIrsphfwQ+nJk3ZeZDmfk9ii6LbfqfKr+VmQ8AZ1O0vj+dmXdl5g8ownvrquMfoGjp/ikzTwB+S7F2ZK1jgG9l5mWZeX9mngZcARw1xJ8nKfrNZwA7UXwY3BkRN9Y+V7ty/AuZOScz5wAvB04CjsjM2yqh/QngXyt/9nuAQym+Vb9v6KuqTmaXi1oqMy+laCG/DNiN4jnUxwCPAd/JzKsiYreI+BbFjdPXAhOA0VWnqe6CeL7yz4eqtr1A0WWxxOzMrJ5R9zv6RqVU2xp4fUR8sWrbWODmpfx57gMOiohRFN8A9qYI5suAHQd7TUSsTrFu5WWZeXZl82aVmn8VEdW1rkJxHaQBDHS1RERMAg7OzCMBKl0sVwNXR8Q1wDuA70TEccCRwDkUXRT/BtxQc7qXBnmLxUt5+9rjRwOLBjluDEV3y8ya7fMHO2nlQ+fSzLwlMxcDvwd+HxF/oPjQWmeIes6tnPPjNe8N8HaKbxjVnhriPOpwBrpaZTRwRERckpk31ex7hr7Q/Djwmcw8FyAi1gTWp1iNfUXVPmzp9Qze6r4H2CQz71+yISJOoOjyOW2Q4/cA1gJuqdn+DEXf+LPA6tU7IuJoin7y12Xmi1W77qf44FkvM39bObYbuAj4LsUSaVI/BrpaIjNvj4grgcsqXRq/puhHfifF6JM3VQ79OzA1In4DvIJi1Ec3/btQ6vX6iDiWosW/H7A9xY3SWicBF0XE3cB1wLuBo4G9hjjvl4HLI+IFilb308BrgG8C387M+RF9vSURsTvF41kPAOZFxIQl+zJzTkR8F/h2RLxI0e//JWAKxTcWaQBviqqV9qNobX4KuBP4DbAr8LbMXLLi+oeBV1PcHP0xRTfG5QxsZdfj2srr76CYzPSOzLy39qDKUMZPUIyUuZvipuQHM/Png500M6+iCPutKT4A7qYI7DMohj7W+heKRtUlwF8pRrEs+YFiFu0VwHkUK9C/CtijZkik1MunLaqjVMahr56Z+7S6FqnZbKFLUkkY6JJUEna5SFJJ2EKXpJIw0CWpJAx0SSoJA12SSsJAl6SSMNAlqST+PzfR5qUzABF+AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "results.set_index('Sample Size').plot(kind='bar')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Properties of MCMC\n", "\n", "This section draws **very heavily** from Patrick Lam's excellent course notes [here](http://patricklam.org/teaching/mcmc_print.pdf)\n", "\n", "The last notebook asserted that by using the Metropolis Hasting algorithm to construct a Monte-Carlo Markov Chain, is akin to drawing randomly from the posterior distribution of the parameters. Why is this the case and what conditions do we need for it to hold true? \n", "\n", "There are proofs showing that with appropriate choices of the proposal distribution $q(\\theta^P_{t+1},\\theta_t)$, the sequence of values we construct using MH is a markov chain approximately equal to taking random draws from the posterior: $p(\\theta|\\mathbf{y,x}$) [see page 124 \"Contemporary Bayesian Econometrics and Statistics\", John Geweke].\n", "\n", "### Markov Chains\n", "\n", "> is a **stochastic process** in which future states are independent of past states given the present state\n", "\n", "* **stochastic process**: a consecutive set of random quantities defined on some known state space $\\Theta$\n", " * $\\Theta$: set of all possible parameter values\n", " * consecutive implies a time component defined by $t$\n", "\n", "#### The Strong Law of Large Numbers (SLLN)\n", "\n", "Let the sequence $\\theta_1,\\theta_2,\\ldots,\\theta_N$ be iid distributed random variables where $E(x_i)=\\mu$. Then we know that\n", "\n", "$$\n", "\\frac{\\theta_1+\\theta_2+\\ldots+\\theta_N}{N} \\to \\theta_\\mu \\text{ as } N\\to\\infty\n", "$$\n", "**Note: N is chain length, not size of your dataset**\n", "\n", "So think of these $\\theta_i$ as being draws of *parameter values* from the posterior. Remembering back to last time, we see that there is likely to be autocorrelation, so that our sequence isn't totally iid. MH samples (if done correctly) are nearly independent but not completely independent from one value in the sequence to another. \n", "\n", "#### Markov Chain Dependence and the Ergodic Theorem\n", "\n", "Since the sequence of parameter values generated by MH are not independent, we can invoke the **Ergodic Theorem** to be confident that our sequence can be used for inference. \n", "\n", "> Let $\\theta_1,\\theta_2,\\ldots,\\theta_N$ be N values from a Markov chain that is *aperiodic*, *irreducible*, and *positive recurrent* then the chain is ergodic and $E[\\theta] < \\infty$, then we know that\n", "\n", "$$\n", "\\frac{\\theta_1+\\theta_2+\\ldots+\\theta_N}{N} \\to \\int_\\Theta \\theta \\pi(\\theta) d\\theta \\to \\theta_\\mu \\text{ as } N\\to\\infty\n", "$$\n", "where $\\pi(\\theta)$ is the probability of transitioning from to another state given a current state (where state is defined by the current value of our parameters). In a MH context, this transition probability is governed by the MH accept reject criteria. This also generalizes to other moments, like standard deviations, etc." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Aperiodic\n", "\n", "A Markov Chain is said to be **aperiodic** if the only length of time for which the sequence of values repeats a cycle is the trivial case of 1 (a single value in the sequence is repeated). The following shows a **periodic** sequence of values, where $\\theta$ always cycles in order $[.5,0,1]$ \n", "\n", "![](../site_pics/aperiodicity.png)\n", "\n", "So if we see a trace plot that looks like this, we need to be very worried about our MH chain failing due to periodicity:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "ExecuteTime": { "end_time": "2016-03-02T07:50:16.900847", "start_time": "2016-03-02T07:50:16.586370" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAENCAYAAADzFzkJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXmcHFd17789u2bRMtpXS7KkI0vehLHxgvFuI4MNj6yOn5MYYiC850ASyGIgJARISIAHJMFOwE4gdpxAMDhGluUdjHeDZVuLj/ZdMxqNttGMZp/3x63bU9PqnqmZ3qrl8/189LGnqrr6Vvfte+4993fOSQwMDGAYhmEYUSgrdgMMwzCM0sGMhmEYhhEZMxqGYRhGZMxoGIZhGJExo2EYhmFExoyGYRiGEZmKYjfAiC8iUgnsASYC81S1eYTrU/Xb/cBx4BfAX6jqz7Nszw7gP1X1z7K4x78BS1X1QhGZD2wHVqrqI2O41+XAU8Gfi1R1a5prngCuBH5fVe/Kot1/CXxUVWeM9R5jfN9y4CPAbwNLgQFgPfBNVf1+6LqngSZV/c0s3usvKcIzGqPDVhrGcNwIjAMOAr8X8TV/CswM/s0B3gm0A2tEZF6W7Tkf+EKW9wizG9fOJ7O8Tw/wG6kHRWQm8K4s7100RKQaeAL4c+DbwIXAJcAa4D9E5K9z/JZfAc7K8T2NHGMrDWM4fg83k24CPiwif6uqfSO85piqNoX+3i8iHwH2Av8L+MZYG6OqLWN9bYb79eGeLVsexRmNL6Uc/03g58DlOXiPYvB54G3AWaq6M3R8g4j0A58XkftU9c1cvJmqHsetTI0YY0bDSIuIzAWuBW4HFPgw8F7gwTHcrjf4b2dw7wTwx8DvA7OArcCdwLdUdSDkNvo08H9xLpGLgZ8Sck+JyEW4ge384P6PAZ/0A1zgXvsizrVSD/wHUBN6Rv8+SfeUiPwqbma9DLfC+i7wuRGM5f3AvSKyNGUAvQn4F1KMhoj8LvAHOHdPGc7d8xlVXR2cfxrYBiwBzgb+IvUNReRm4N+A21X1rsCNdDvue1qAM4bfBr6sqn0i8hjQoarvC91jLrADuEFVH065fyVu0nB3isHwfANnEHeEjtWKyJ3ArwPVOGP6+96tKSIXA3+F+75qca7PO1X174Pzf0nIPRW4Oz+Cm2xcBrQAPwT+RFV7MYqCuaeMTHwQN1j/EHgaNwh9bLQ3EZE5wD8CbYAfmL4IfAL4JLAc+Gvgc5w8U78NuA74QOrAJSIXBO3ahXOBrcQZoGdEZGJw2TeAW3HG6e1AF2ncSKF7vh/4PvBj4BzcgPX7uIFuOH4BbMatLPy9FuEG/P9OeY8bcYP5nTjDdDFwAGd0xoUu/V3gX4ELgP9KucfNwD3A74X2Sb6KM6B/j3PxfA7nKvx6cP4eYKWITAnd6hbc97omzTMtBBqB59I9sKoeV9Wfqmpn6PANOFfkO3AD/RXA3wVtnokzIho805nAD4C/E5F3pHuPgK/gvpNzgH8A/hD4rWGuN/KMrTSMkxCRMtxg+2RolvifwMdFZLGqbh7m5V8Xka8E/18JVOFm0r+qqrtFpA73w79NVX8UXLdNRBqBr4lIeIC+U1XfyPA+f4wbqG9T1f6gjb+CWzncKiLfAT4E/LF/HxH5A4bfY/gU8KCqel/9psC1FmVj9n6cQfrL4O/fAlar6hERCV93KGjzvwV/7xCRr+EG1DnBMwG8qap3+xf5ewQG4zvALX4jWkTG4wz6Z1T1X4OXbBGRScBXReTzwAPAPwVt/KfgmluA72ZYRTUG/z0c4dk9r6vqJ0Pv/wOcQQe3wvtr4Kt+lSAinwX+BGdcX8xwz/tCz/QVEfkQbl/le6Nol5FDbKVhpOMa4DTcQOi5F0jgZt7D8UXg3ODfUmCCqp6pqo8G55fhBpB/EZHj/h/wteD4gtC9hjNOZwM/9wYDINhLUdysVHAG66XQ+QEyzJwDzgFeCB9Q1f9W1X8c5jWe+4GlInJ28PdNwH2pFwUKsudE5NMi8j0ReRa3sgEoD12a7tkn41YfFTj3lecMnIH+Wcr1Pw3uuVxVu3Duuf8NyZXa0uB+6fD7R5MznE/HppS/D+GEFKjqdtwK66Mi8s8i8jjOPQVDnzuV1P2So7jv1SgSZjSMdHil1LdFpFdEehmcCf5uihsllRZV3RL8266qx1LO+z53C4PG5VycS2Uxbn/Dc2IMbS/DuaE8iZTz3cO8drhzwxLsZawFfkNEzsW5yn6Sep2I/AZu5bUM59b6LM7ApJLu2QeA9+OMwXcDddNw+M/afx53AxcGrrNbgGeGWTVuA5qBi9KdFJEGEXlSRN4dOpxuxZIIrj8DZ9B/DbcP8k3c9z4SXWmOpX6nRgExo2EMIfB534jzo5+b8u8zwCTSD3JReRMnUZ0fMi5bcH7wLxB9QHgdeGfgSvNtn4nbPF4PbMRtvF+a8roLhrnnhtTzIvIJEVkbsU334wbFm4Afpvj7PZ/GuVxuVtVvqOqTwPzg3EjPfijYsL4NtxL0brSNuM801fV2GU6EsAlAVV/FGbabgnbek+mNghXcd4APBhvmqdyO27NIt0mejt/H7Wtdrqp/o6r/A/j9FTMCJYTtaRip/A5u+f9lVV0XPiEiW3Eb2B9jmAFnOFT1aKCw+ZyIHMZtZp8HfAv4L1XtStkDyMRXceqdbwd7Ag3BsVbgXlVtF5FvAJ8VkX3AK7jN/fOD/0/Hl4CfiMif4PYAzsCtBKK4pwD+E/hb4KPABzJcsws3278A5wK6msHBf6SVA+BcPSLyadwe0I9V9TkR+RbwGRE5gPtcLsbtr3xHVVtDL78b50Isw21ED8cXg/Y9F+w//BynQrsJJ2L4jKpujNJm3HPPBG4UkVdxAoivBeciPbcRD2ylYaTyIdxewaupJ1T1BG4APS8Y9MbKH+EG10/jVh5/j3NX/J+oN1DVl3B7L0twRmAVLljvYlU9FFx2B05983fAazj313eGuefDOHnubxNEPQP/j5HVU/71u3B7Ju0MRoqn8n9x7pnHce6p38EppU7gVltR+Qfc/su/iUgtThjwZZxqagNOpvu3wfuFuQ83KfgvVW0f4XlO4FYTd+K+s1/iAiEvAn5NVb84yvZ+D/f5b8R95/+M24cZzXMbRSZhlfsMwzCMqNhKwzAMw4hMzvY0AnfFj1V1Vobzq3CJ25IKC1WtD86twC1Vl+Okhh9V1RfS3ccwDMMoHlkbjSAlxK24Ta3hQvvPBS5V1SGbkCJSAzyE23T7Dk4K+ICIzFfVMUsgDcMwjNyTi5XGHbhcM1/EpS04CRGZBkwD1qU5fQXQr6p3Bn/fIyJ/iEtJ8MOR3jzQqp8P7Ce9TtwwDMM4mXKcou3lIPgzErkwGvfgpIqXDXPNCpxG+ycicg5ON/5JVX0eF5W6IeV6xbmqRjQaOIPxzGgbbRiGYQAulilyrZusjYaq7ofB3DgZqAGex+WZ2YLTy68WkaVAHdCRcn0HLgtmFPYD3HfffcyYYbVbDMMwotDU1MTNN98MwRgalYIE96nqgwxNqX2niHwM55rqIMhPE6KW6Hn1+wBmzJjBnDlzsm2qYRjGW41RufULIrkVkV8VkV9POVyDS/OwEZdcbshLONllZRiGYRSZQqURqQf+VkTW4SS1n8CtLh7FKa6qReR24C6cemo66XP8G4ZhGEUkbysNEblLRO4CCGoHfAN4BDiCS4i3UlXbg137lbh8NodwidBuHCnFgWEYhlF4Sj6NiC/Z+cQTT9iehmEYRkT27NnDVVddBbBAVXdEfZ2lETEMwzAiY0bDMAzDiIwZDcMwDCMyZjQMwzCMyJjRMAzDMCJjRsMwDMOIjBkNwzAMIzJmNAzDMIzImNEwDMMwImNGwzAMw4iMGQ3DMAwjMmY0DMMwjMiY0TAMwzAiY0bDMAzDiIwZDcMwDCMyZjQMwzCMyJjRMAzDMCKTsxrhInIB8GNVnZXh/G3An+DqfyvwR6r6THDuU8AXge7QS1b684ZhGEY8yNpoiEgCuBX4GtCb4ZorgC8B1wCvA7cAD4nI6araCpwL3KGqX8m2PYZhGEb+yIV76g7g47iVQibmAH+vqmtVtV9Vvwv0AcuD8yuAtTloi2EYhpFHcuGeuge3irgs0wWq+u/hv0XkEqAB2CAitcAS4OMici9wGGdg7slB2wzDMIwckvVKQ1X3q+pA1OtFZBnwQ+AvVPUgbo/jWeBOYB7wYeBrIrIy27YZhmEYuSVnG+FREJFrgf8CvqqqfwugqtsZukp5RkT+HXg/sLqQ7TMMwzCGp2CSWxG5Ffhv4GOq+oXQ8beJyJ+lXF4DdBaqbYZhGEY0CrLSEJGrgG8B16aR0R4HPiciW4AHgCuA32SYPRLDMAyjOOTNaIjIXQCq+lHgT4EqYLWIhC/7VVV9RER+Hae++i6wB7hVVX+Zr7YZhmEYYyNnRkNVnwamhP7+aOj/rx3htQ8BD+WqLYZhGEZ+sDQihmEYRmTMaBiGYRiRMaNhGIZhRMaMhmEYhhEZMxqGYRhGZMxoGIZhGJExo2EYhmFExoyGYRiGERkzGoZhGEZkzGgYhmEYkTGjYRiGYUTGjIZhGIYRGTMahmEYRmTMaBiGYRiRMaNhGIZhRMaMhmEYhhEZMxqGYRhGZHJWuU9ELgB+rKqzMpy/CVfSdRrwNPAhVW0Ozl0NfB1YAPwyOLcpV20zDMMwckPWKw0RSYjIB4FHcXXA011zNnAXcBMwFWgCvhWcmw48APw5MAl4HLg/23YZhmEYuScX7qk7gI/jVhGZuBl4UFVfVNUTwJ8C7xORacAHgLWq+pCqdgNfABaKyHk5aFtJ8dL2Qzy+obnYzcgpAwMD/OjVPWzcf6zYTckph9u7+fcXdnKss6fYTckpG/Yd48G1exkYGCh2U3LK4xuaeXnHoWI3I6ec6O7je8/vYM/hjoK+by6Mxj3AucDLw1yzFNjg/1DVVuBIcDz1XB+wFVieg7aVDMc6e/ide17i9773Cr/cdbjYzckZT2w8wB/+12vccvdL9PT1F7s5OeMrjyqf/fE6/vLB9cVuSs4YGBjgtu+9wsf/cy0Pvb6/2M3JGdrUxu997xVu/s6LtLR1Fbs5OeO7z+/gLx5czxd+srGg75u10VDV/ao60rSkDkg1hx1A7Qjn3jJsbm7jRE8fAA+9tq/Irckda3cfAeDg8S6e39pa5NbkDv9cj6xvojP43kqdluNd7D1yAoD/Wbu3yK3JHa8F31V3bz9r1jcVuTW5wz/XuKrygr5vodRTHcC4lGO1wPERzr1l2Nw8+LiPrGuiv//UcA9sPtCW/P/V606N2Wtf/wBbDrjvq6O7j59uailyi3LDllAf/Nmmg7SdIq63U7EPAmwO+uCiafUFfd9CGY2NgPg/RGQK0BgcTz1XDiwi5LJ6K+A7AMD+o52s3XOkiK3JHeHnWrO+md5TwEW19/AJunoHn2P1G6fGQBT+rrr7+nnyzQNFbE3uCD/XC9sO0Xq89F1U3b397DjYDsDiU9Ro3A/8ioi8U0RqgL8BVgd7Gz8C3i4iHxCRKuAzwB7g1QK1LRaEOzbAw6eAT7mrt4+drYOex0Pt3by4vfQ3I8MzV4DHNx44JVxUqc+16hTogzB0Fd/XP8Cjp4DYZGdrO72BN2Lx9IaCvnfejIaI3CUidwGo6lrgNtym+QFgFnBrcK4JeB/wOaAVuBr4QIR9klOKLc3uBzulvhqA1euaSl7BsuNgB31Bx/bP9fApMCv3Bn5ynVOYH+/q5eebDxazSTnBD67+u3p6UwvHu3qL2aSsae/qTe7TnIp9sKqijHmNhd3+zZnRUNWnVXVK6O+PqupHQ39/X1WXqOp4VX2Pqh4InXtKVc9R1QZVvfStFtjX1tnDvqOdAHz0soUA7D1ygtf3HC1ms7LGz1xrKsv47YtOA2DN+qakISlV/OD69vmTePtpkwB4+BTwlft9mg+9cwHlZQm6e/t5qsRdVFtbBlcZ/rf13NZWDrd3F6tJOcH3wdOn1lNelijoe1sakRiwJeSauvGcWZw22c0cSn1GFO7Y7z17JgAHj3fzUom7qLYExnDxtAZWnuWe67ENzXT1lq6LqvV4F63BQHrBgklcfPpk4NTpgw01FfzmBfOoqiijr3+Ax0rcRbU52QcLu58BZjRigV9qThhXydSGalae6Qaih9ftL2kXVbhjL5xaz9IZzvdaygqW/v6B5Pe1eHo9K8+cAUBbZy/PbSldSXF4T23R1IZkH3xKD9DRXbouqk2hPlhfXcFlS6YCpb8y9MbQjMZbFL/SWDytnkQiwXuC2evuQydYv690I6mTHTvYqLs+eK7VJSwp3nf0BB3dbkWxaFo9syaOY8W8iUBpz8q90ZjWUM2E2kquWz6dsgR09vTztJaupHhLcnB1fdD/tp7dcpCjHaUpKe7t62fbwcGJS6ExoxEDNgeb4L4DnDl7PHMmudCVVSU6EPX09bM9kAR6Hfn1Z7lZeUtbF6/sLM2odz+4JhLO7QZwfTArf3RDc8lGvW9J6YOT66u5cKFzUZVqHwSGrAoBrjxjGlXlZfT0DfDYxtJ0Ue081EFPn5t0LZpWWOUUmNGIBYNBOq4DJBKJwVn5G6XpohoiCQyMxqJpDcn/L9VZuZ+5zmuspabSReK+O3BRHT3Rw3MlGvWeHFxDg5Dfr3nqzQOc6C69/ZoT3X3sDvIy+YnL+JpKLl3s9DqlGl/jV/CV5Ynk/mchMaNRZDq6e9lz2EkCl4SWmt5o7GjtYOP+trSvjTO+Y6dKAv1zlWrU++bQJrhnbmMt58yZAJTwQHTgZHfHu5fPIJEo3aj3rS3H8fOtJaFYBt8Hn9l8sCQTTnohxoIpdVSWF34IN6NRZLYeaE/+f3ggOmfOBGZNqAFKc+PYD0ILp9RREerY/gfbdKyTV3eXXtR7usEVBmfla9Y3lVzU+5GO7mQiv3AfnNpQzQXzG4HSXBn6vcK6qnJmBr8lgKvPmE5leYLuvn6eKEEX1WAfLLxrCsxoFJ1NgS+5obqC6eOrk8cTiURyIFpVgi6qTB17yfR6Fk6tA0pvIBoYGAhtrA41Gn5f43BHDy9sKy1JcVjyfdJzBX3wiY3NJRf17leFi6Y3kEgMxjJMqK3kkkXORfXwG6WXwLCYyikwo1F0kvsZ0+uHdGwY/MFua2lnUygVQimQ3NxP6dhhdVip7dc0HeukLYiQXpyyATlvci1nzh4PlJ6cc1MyEryKSXVD66i9+0znomrv7uOZEot63zTM4Op/Wz8tsaj3vv6BZMBiah8sFGY0isyWYYJ0VsydyIzxblldSrPy3r5+trVkTqbmYwD2He3ktRKKeg/nMDp9Wt1J5/1zrVlXWlHvyRl5mu9q+viawaj3EuqDMFTKnsq1y6ZTEUS9l5KLavehjmSyzGLIbcGMRtFJp1rxlJUlksqcUvrB7jrUQXdf5o59xswG5pdg1Lv/ruZMGkdtVcVJ5/3stbW9mxe3l46KasswfRAGjeHjJRT13tnTx87WYOKSpg9OrK3ioiDqfXUJuah8HywvSzB/8skTl0JgRqOIdPb0seuQkwRmmjW8J0i/sfnA8aTLJ+74ju0kgSd37LCk+OESclH5VeGSDBuQC6bUccZM56IqqYEoWEEtydAHVwbxNW1dvTy7pTRcVNsPtuMXe5mMoXeTPqUHaC8RF5VfFc6fXEtVRXGGbzMaRSQsCcykhDhv3iSmNQxmvi0F/Mx1OEmgNxp7Dp9g3d7SiHqPsgF5fbAyfKREEjMe6+yh6ZhLlpkpUGzmhHG8LYh6X/V6afRBP3EZV1nO7ImpNd4c1y6fQXlZgq7efp7S0kjMuCVp4IuznwFmNIqKH4TqqsqT8tpUStFFNbgJnrljL581nrmNpRP1PjAwQJRKaV7x1tLWxSs74q+iGqKcGsZHfn0yMWMT3b3xlxT7CPdF0+opy5AFtrGuigsXOklxqawMNw+zT1MozGgUkfAGZKpyKoz/wb7Z1DYk1XNc8aqV4QbXIVHvJZCYsaWti6MnXCDYcPr4RdPqkek+MWP8ByJv4CfVVibrg6TDG8Njnb08tzX+LqrhlFNhfB98sgSi3vtDZYYX2Urjrcnm5OA6fAc4f34jU+rdD/qRmA9EQySBI6g7fGzDztYONuyPt4tqSBbYEQYivwewet3+2Ee9bw4l9Btu4jJ74jjOmVs6iRkHYzSG/66uXTaDsgSc6Onj6Zi7qPYeOcGJIFbGVhpvUbZkiC5OpbwswXXL3UAU9xKcew6HJIEjGMOz50xI+pvjPhD5GfmsCTXUV5+snArjZ6/Nx7r45a54J2YMxwmNhN+viXtixu7efnYEZYZH6oNTG6q5YEEQ9R7zCZk3hGUJt19YLLI2GiKyQkReEpF2EVkrIhemuWa1iBwP/esQkQERuTg4/08i0pVyzbxs2xZnunr72BFIAjOpVsL4gWjD/mPJgvJxxM9cy8sSI3bsRCKRrEfx8BvxLm87mtQNi6fVc3oy6j3eA9FwsQyp+D54pKOHF7bFV1K8o7U9KUIYzXM9GfOod//bOm1yXTJZZjHIymiISA3wEPCvwETgm8ADIjLEOaqqK1W13v8DfgD8h6o+F1xyLnBz+BpV3ZVN2+JOFElgmHcsaKQx8DnH2VfuB9eoksDrA0nx9oPtaIwlxaPZgBwS9R5jF9XxUP3sKGqcuY21nDXbJWaM88rQD67VFWXMjVA/2ydmbI95YsYoQoxCkO1K4wqgX1XvVNUeVb0HaAVuyPQCEXk/cCXw0eDvMuBsYG2WbSkp/EZdTWVZRklgmIryMq5bPh2IdwLDdFlgh+PcOROTyeTiPCuP6kr0+I3j/Uc7eW1PPBMzbh0m51Qm/H7NmvXNsU3M6Ptg1PrZ08bXcP5pXkUV59/W8PE0hSJbo7EU2JByTIHl6S4WkQrga8AnVdVPKxcDtcBXRKRFRF4Vkfdm2a7YE0USmIqPzH19z1F2B0GBcWO0g2spSIpbj3dxKKifHbXozdIZDUn3XFyfyw9C42sqmNpQPcLVDi9eONQe31rvmTIRD4c3ho9vPBDLqHeXLHN0E7J8ka3RqANSR68OnBFIx28AnTj3lGcS8DTwd8As4K+A74vIWVm2LdYMlz4kExedPpmJtZVAPFcb/f0DIUVY9B+sd+VsiWnUezhZZNTncpLieO/XDFaMHF45FWb+lDqWBVHvcY2vyZQsczj8hOx4Vy/PbIqfpHjf0U7aQ2WGi0m2RqMDSPWt1AKZggluBf5FVZPrWlV9QVWvUtXnAhfXj4EngFN6tTEW/2RleRnXLnMuqji6coZKAqMbw7eFot7j+Fw+fcj08dVMGFcZ+XV+INp75ARv7I1fYsaxBopdn3RRxS/qfWiZ4eh9cMaEGs7ziRljOCHzhjBcZrhYZGs0NgKSckw42WWFiDQAlwHfTzl+lYh8JOXyGtyK5JSku7c/qYAabToA7ytfu/tIchMzLnjXVFmCZM2MKJSVhVVUMfzBJn3Jo/uuls8an6xaGMdZ+XDZbYfD98GDx7t5OWZR7ztbB+tnjzYLrO+Dj21ojl3Uu/9tzZ1Uy7iq4imnIHuj8SRQLSK3i0iliHwQmA6sSXPt24F9qrov5Xgf8FURuVREykXkJuAdpBiXU4l09bOjcsnpUxhf4+IE4rZp5wehsUgCvexRm9uGpLaIA2NxuUFqrfd4uaiGlhkenTE8fWo9S2e418TNyPtVYVV5GadFUE6F8cawrTN+iRmLXXgpTFZGQ1W7gJXATcAh4HbgRlVtD2Iz7ghdPh84qYep6tPAx4F7gGPAp4AbVHVvNm2LM5tC9bOjSALDVFWUcc0yH3EcL1fOWAdXgLfPb2RKvXNRPRIz98BY9p883pWz61AH6/fFJ+p9W0t7KFnm6L8v73pbHbNa774PLpw6tMxwFGZPHMe5MY16T6oSi5g+xDN8aGsEVPV14OI0x1em/P2vuHiOdPe4G7g727aUCqOVBKZy/Vkz+OEv9/CLnYfZf/QEMyeMLNktBNkkUysvS/DuM6dz7wu7WPVGE//3ysW5bt6YONzezcHjQf3sMQyuZ82ewJxJ49hz+AQPv7GfM4M4h2Lj+2B9dUWy0NdouP6sGfy/xzfR0tbFL3Yd5vyglnixyTaW4fqzZrB29xEe3dDMl/r6M2ZpLiThZJklv9Iwxka2HeCdi6fQEKSyiEsuqoGBgVHLbVPxrpyN+48lNzOLzZCcU2PYgIxr7ZBwUsmoyqkwi6c3JPtvnFLbbMpSlupXUEdP9PDc1nhEvTcf66KtMygzXOQYDTCjURS2ZOmfrK4o5+pARRWXlM77j3Ymay2P9Qd7wfzGZKbVuEiK/Yx8Sn31SfWzo+I3WHe0dvBmUzwkxbnwkfs9gEdi4qLq7etn28HM1fqiMLexlrPnuNVgXPYMfR+E4iunwIxGwXEdO3oeo0z4gejlnYc4cKz4QjM/I89GElhRXsa1y+Olohqpql0Uzp07MVkvJS7PlaxNn8Vz+f2apmOdvLq7+FHvuw+fSKqesjKGvtb7+qZYRL37Pjh74jjqRkiWWQjMaBSYnYfGLgkM864lU6mrKmdgwFWJKzZeR56tJNAH+q3be4xdrcWPeh9NQr9MJBKJ5Kx8VQxcVEPKDGcRXSzTG5LS6jgYQ98HK8rSlxmOijeGhzt6eGFb8SXFY4lwzydmNAqMnzVUlidGLQkMU1NZzlVn+EC/4v9gczG4ArxjYSOTYhT1PliXITvVih+ItrW0D9knKQbbWkLJMrMYiBKJRDKtyOoYGEP/uS6YUpdV/ezTJtexfJaLeo9DoF9yVRiDTXAwo1Fw/Gxo4ZT6UUsCU/ED0UvbD9HS1pV127JhNHUZhsNFvQcuqiJv8h890UPzsUA5leUPdsXcSUwf76PeizsQeUNYW1XOrCyVdz5n076jnby2p7hR79kKMcJ48cKadcWNeh8YGBisQhgDuS2Y0Sg4uRpcAS5bMo1xleX0Dzj/a7EYGBiIVBc8Kn4gem33EfYcLp6LassYssBmwkW9D6qoismWkCw1arLMTCybOZ7TJrsVc7GfazDCPQd9MNgzbG3v5sXtxVNRHTzePVhm2FYab01yqbceV1XOlWdMA4rryjnQ1sUxLwmCIDgEAAAgAElEQVTMwXNdsmhKMsdTMSXF3hA21lUxuT5aFtjh8LPXTc3Hky6HYrCpeWzpQ9IRF0lxuH52LvrgwlDUezEViuEEnsVOVOgxo1FAwvWzR5u6IRPep/zCtkO0Hi+Oi2rzGLLADkdleRnXLCv+fk2uA6rOO21SMgV5UQeiLCLc0+H74J7DJ1i3tzhR73sOn6CzJ1BO5WjD2BvDR4qYmNF/VzMn1NBQEz1ZZj4xo1FAdh/qyIkkMMwVS6dSU1lGX/8Aj25ozsk9R4t3C+RSEuj3a3656wj7ipSYMdeqlfKyBO/2td6LZAy7evvYmayfnZvnOnP2eOZMcnsjxXqufNTP9n2wpa2LV4qUmHGsSSXziRmNAuIHoWwlgWFqqyq4QpyLqliz8nxIAi9ZNIWGmuJGveej6I2fvb7Z1Ma2lsKrqHYc7Bisn52j7yu1vG0xXFSDZYbrqK7ITRbYRdMakvE5xcrzNhiEGY9NcDCjUVC8L3l+lpLAVHwMwHNbWzkcVJgrJNlGuKejuqKca84oXnnbts4e9h11QZO5fK4LFoSj3gs/EPmZa01lGXMmjV3ynYrvgztbO9iwv/AuquTgmuNYhsHEjMWp9b7lQPbBpbnGjEYByeVGXZgrl06jqsK5qB4rsItqYGCATaOsCx4VPxC9svMwzQWOet/aMpj7KhdKN095WYLrilg7xA+uY02WmYlz5kxI1rovxnNtyVMf9CvD5mNd/HLX4ZzeeyRaj3fRGkwC4xLYB2Y0CsrmPAXp1FdXcPmSqUDhg5Fa27s50uEkgbkcXAEuXTyF+uoKF/Ve4Fm5V61MGFfJ1Bwop8J4V876fcfY2VrYxIz5mrgkEuFCWoWtHTIkC2yO++CS6fWcnox6L2wf3DIkWaa5p95yDJEE5iFIx8+Int1ykKPBIF4INoUkgbkeiFzUe3H2azaH3AJjyQI7HO9YEI56L+xAtKk5f3UZ/Mpw+8F2tIC13vceOUFHnupnDymkVWAX1aagD05rqGZCbTyUU2BGo2DsPZJ7SWCYK8+YRlV5GT19Azy2sXAuqi15lgR6n/JLOw5xoK1wLqrNzbkLFEuloryM64qQmHFo/ezc98EVcycma3MUclaei2SZw+H74P6jnby2p3CJGZNCjBi5psCMRsHIhyQwzPiaSt61ZApQ2JTO2VTri8LlMpXaIDHjmvWFM4b5LnrjZ6+v7znK7kOFiXrPpsxwFMrKEslo/kIaQy/EmNdYO+oyw1E4Y2ZD8jdbyOfKdTxNrsjaaIjIChF5SUTaRWStiFyY4br1ItIhIseDf+tD564WkXXBPZ4RkSXZtitu+PwxuZQEpuJnRM9sPsixzsK4qDbnaQPSU1NZzpVLg6j3Av1gw/Wz8zXLu+j0yQWPet8cKjM8L4tkmcPhjeGWA8eHRDPnk3z3wWLt12RbhTBfZGU0RKQGeAhXxnUi8E3gARGpSrluHCDAaapaH/xbHpybDjwA/DkwCXgcuD+bdsWRfM/IAa5eNp3K8gTdff08ufFA3t4nTC6TxGXCD0QvbGstSNT71gODm9P5GohcYsYg6r1A4gU/CC2cMvr62VE5b94kpjX4xIwFMoYF7IN7j5zgjb35T8x4pKM7mYQ0V9kjckW2PecKoF9V71TVHlW9B2gFbki57iygSVVb0tzjA8BaVX1IVbuBLwALReS8LNs2Ir19/bQVaEaei6I3IzFhXCXvXORcVIWIzD3U3s3B44EkMI/G8AoJJ2bMv4vKz1wbqiuSmWnzwfVnu4Ho1QJFvW/OoxDD4xIzFs5FNTAwkJc4oVSWzxqfXJ0V4reVy2SZuSZbo7EU2JByTIHlKcdWAD0i8ryItIjIoyJyRrp7qGofsDXNPXLOn/z366z4/GP8Ymd+UwQMLQyf31mDV7D8dFNLsvxqvtg8RDmVv+caV1XOFUsDSXEBfrDJ+tl5UE6FueT0waj3QjzXYCbi/A5Cvg9qc9uQwS8fNB3rpC3LMsNRcIW0nDFcXQAXle+DU+qrxlxmOF9kazTqgNRdvA4gncP0ZeAmYB7wCvBw4LYazT1yStOxTnr7B7jvhV15fZ99RzuTksB8KyGuXTadirIE3b39PPlmfl1UmwsoCfT7Nc9vy3/Uu18VLsmzga+qGEzMmO99jd6+frYFAYv5Nhrnz29kSr0b6B7Js+stnCzz9Gm5F5iE8YkZdx3qYP2+/Ea9xzHnlCdbo9EBpFZxqQWGTC9U9Z9V9ddVdYeqngA+DTQC50a9Rz64OkhT8diGZrp6+/L2Pn6Gly9JYJiJtVVcdPpkAB5+Pb8/2ELsZ3iuXDqN6gqfmDG/A2why2v6geiVnYdpOpo/SfGuQx109+VP8h2mvCyRlBSvyvO+hv+u5kwaR21Vfutnn13AqPctMVVOQfZGYyNugzuMkOKyEpEPi8jVoUPlQCXQmXoPESkHFqXeIx/45WZbVy/PbjmYt/fZnGdJYCo+4vgpPUB7Hl1U+VathKmrruBy8S6q/A1E4frZhZjlXbrERb1DfmflfnCtLM9dsszh8H1w4/5jydiQfFDIUqgu0G9wvyafLqp85dLKBdkajSeBahG5XUQqReSDwHRgTcp1s4BviMjcwCX1NeBN4DXgR8DbReQDgerqM8Ae4NUs2zYiMyeM423zJgL5HYjylT4kE9cun0F5WYKu3n6eTqs9yA2FUISFKUTU+9aW4wwk62fn3xhWV5RztY96z6OLys9cF0ypozJPyqkwFyxopDGZmDGPxrA5t/VpRsLv1+xo7eDNpvxIio919tB0zCfLPMVWGqraBazE7VUcAm4HblTVdhFZLSJ3BJd+EWdIXgIOAKcD71fVflVtAt4HfA6nvLoa+ICqFkQM7QeiR9c3JWtd5JpBvXVhOkBjXRUXLmwE8reMPtrRw4G23NTPjopPzNibRxeVH1zrqsqZNaEmL++Rih+IXt5xiAN5SsyYy3K8UXBR7/ktpBUWmBRq4nLunInMnOCj3vPzXEOUUzFcaWTtBFTV14GL0xxfGfr/HuCPgn/p7vEUcE62bRkLK8+ayRdWbeRYZy/PbT3I5UFtilxRKElgKtefNZNnt7Ty5JsHONHdx7iq3LrFtrSElFMFmuU11FTyrsVTeXxjM6vXNfFrb5+b8/cIr57yqZwKc9mSqdRVldPe3cea9U3cctH8nL9HMQLFrj9rJve/tJt1e4+xq7WDeZNzq21pOd41WD+7QH3Q13q/59ntrHpjP390zZKc9xM/XkyqrUym0Y8Tb/k0IrMnjuOcuc5FlY8SnM3HupKSwEIG6Vy7bAZlCTjR08dPN+VeRRWWBDYWsGN7n/Izm1vyEvWez4R+maipLOfKM/ysPPd9sG9IsszCGY0LF05mYjIxY+5n5bkuMxwV3we3tbQnjXEu2RRaFRZq4jIa3vJGA+D6IBhpzYYmevpy66Ly+xmQf0lgmKkN1VywwLmo8qFgKfR+hueqM1zUe0/fAI/noXZIvlKHj4Tvgy9ub+VgjqPe9xzuoCtZZrhwxnBI1HseXDne5TZrQk1STFAI3jZvUjLoc1UeFIrJVWEMXVNgRgMY3Nc40tHDC9tac3pvPyMvhCQwFa9geXJjM509uZUUF1I5FWbCuEouXZwfFVVXbx87gvoWhfYlXz4k6j23z+UNfHlZgvlT8h7+NAT/23ptz1H2HM5tYsbBwbWwfdC7qCA/K6hiTVyiYkYDmNtYy1mzJwC5H4gKKQlM5brlM0gkoL27j59tyq2KqhjuDo9PU/GzzS05TQOz/WA7vlxCoY3huKpwYsYcG41k/ezavCXLzMTFp09hfJ5qvSdrnhTht+X74Kbm48nfeC443tXL3iClTByVU2BGI4mP2VizvoneHLqoBvXWhe8A08bXcP5puVdRtXX2sD8IRCtGxOo1eYp6999VTWVZMoirkPg++Py2Vg7lMOq9WKtC8FHvPtAvt7PyYk5c3j6/kSn1uU/MuDW0RxKnuuBhzGgE+MjcQ+3dvLQ9N7moiiEJTMVv2j2+8UDOot6HJlMr/EA0sbaKS4LEjLk0huHvqiyH9bOjcoVMo6YyiHrPoYuqmIMrwHvOdn0wl4kZW493JQ1roaTsYcrzlJjR98HxNRVMbchfssxsMKMRMH9KHctmjgdyl6o6LAksVnrjdwfG8HhXLz/fnJuo980hSaDPMVRovDF8WltyFvXuN1bznXMqE3XVFVy+JLeBfv39A0UTLXguWTSFhurcuqg2FUk5FcavDN9samNbS25UVJtD6r04KqfAjMYQ/ED0yLpm+nJQC3hLDDr2jAk1nHfaJCB37oGwu6NYHfuaZYNR77lyUcVBteIHoue2HORIR/Yuqr1HTnAiEEEUy0deXVHO1YGKKlcbx34fYfr46mQxq0JzwfzGZBxFrmq957tiZC4woxHCKz0OHu/i5R3Zu6h8Byi0JDAV/1yPbWjOSdR7HAbXxroqLg4SM+ZiIOru7WfHQZ8FtngbkFedMT0U9Z69pNi7psoSsHBq4STfqfg++MrOwzTnIOo9DqVQK8rLuC7HLqo4Z7f1mNEIsXBqPUtnuE6Yi9KiPkin0JLAVN4ddOy2zl6e3Zq9i2pzESLc0+Flj0+92UJHd3YuqnzXz45KfXUFly1xkuJc9EE/CBUqWWYmLl08hbqg1nsuXFTFdrl5/F7o+n3H2NmaXWLGoWWG46mcAjMaJzGov26iP0sXVVyWmrMnjuPcIOo923Tp7TGSBF67fHoy6j3bxIz+u6qqKGNunupnR8W7SX++5WByT2ysFFO9F6amspyrzshdoF9Sblvk57pwYSOTklHv2RnDbS3tyWSZcVVOgRmNk/A/2ANtXfxi1+Gs7hWnIB0f6Pfohuasot63hjb8ip1MbUp9NRcuDGqHZDkQ+cH19Kn1lBdBORXmqjOmU1VelpOo97hMXGDQRfXSjkMcaBu7i+pwe3cyar7YfdAlZsyNi8qvCuurK5gxvjDJMseCGY0UFk9vSP7AsukEYUlgsWd5MOiiOnqih+e3jj3q3Q+u42sqmBYDSaDPEPvkmweyinr3P9g4zPDG11Ry6WInKc5mv2ZgoDg5pzJxuUylNnBRZVPrfUto4rIoz0XNouD74Ot7jrL70Nij3ouRLHMsmNFIg+8Ej2ThogonMiu23xVc1PvZc3zU+9gHok0H4iUJvG75dBIJ6OjOzkUVl30aj++DP9t0cMyJGfcf7UzWiS+2KxGci+qKZNR7Fn0w2CucUl8di/rZF58+Oangyma/ZlPM+mAmzGikwbuo9h/t5NXdR8Z0D280iikJTMW7B7KJei9GmvfhmNZQwwXzXdT7WGflvX39bDtY2JonI3FNkJixu6+fJzeOTVLs+2AhygxHxbtJX9g29sSMcTPw4cSM2cjakymHYrAqHA4zGmmQ6Q1JeeJYZ0SFLnoTBR/BerijhxfHGPVe7Aj3dHhj+MTGsbmodh7qoKcvUE7F5Ac7obYy66h33wfnTBqX83oqY+VymUpNZRn9A/DoGF1UcXK5eXwfXLv7SFIoMhrCZYbjNGakw4xGGhKJRFJKt3pd05hqAcdFEhjmtMl1LJ/lot7HMiM60d3H7iBTaRz2aTzvPtMlZjze1cszY4h6999VZXmC04qsnArj++DTm1qSbqbRsCWZ0C8+31VtVQVXBIXOxroy3Hyg8DVPRuKSRVNoyCIx47aWwWSZcRoz0mFGIwM+MnfvkRO8tufoqF+/OYazIQi5qNY1jTrqfUj97Bh17Onja3h7EPU+lpWhdwssnFJPRQHqZ0fl2uXZJWaMQxBmOnwffG7r6BMzHj3RQ/OxwpYZjoJLzDh2SbE3hLVV5UVJljkasv6FiMgKEXlJRNpFZK2IXJjhus+IyC4ROSIiT4vImaFz/yQiXSJyPPRvXrZty4ZlM8czPyhPOdqBKCwJLLaOPBXvomodQ2JGP3Otr65I1kmOCz6+5rGNzaNOzBhXAz+xtoqLfNT7KPvgwMBALF2kAFcsnUZ1hUvM+Ngoa70PTZYZr+/Lrwx/sfMwTUdHJyneUuRkmaMhK6MhIjXAQ8C/AhOBbwIPiEhVynW/C/w2cDkwBXgcWCUi/v3PBW5W1frQv13ZtC1bEolEUsHy8Lr9o3JRxU0SGCYc9T7aGVE4xUEclFNh/MqwrbOXZ7eMzkU1uLEar8EVBmflT+mBUUW9t7R1cazTK6fi1QfDUe+jTSvuV4WNdVVMri++5DvMOxdPSaYLemSUrrc4urMzke1K4wqgX1XvVNUeVb0HaAVuSLluCvBFVd2mqr3AN4B5wJzAcJwNrM2yLTnHzxx2HzrBur3HIr9uc7J+djwkgal4Bcsj60fnooqzJHDmhHG8bV4Q9T6KgaivfyAZsBi3lQa4QlrlZQk6e/p56s3okuJwFtjTY/h9veds1wefHWVixk0xHlxrKsu5+owgS/EojeGmItY8GS3ZGo2lwIaUYwosH3JA9Suq+t3QoRtxxmUPsBioBb4iIi0i8qqIvDfLduWEM2ePZ26j8y+OJl36YGH4+HVsGIwBaGnr4hc7o0e9x1G1EsbPyh9d3xQ5MePuQ+H62fF7rsa6Ki5cGBTSGkUf9KvC2RPHFTVZZiauXDqNqnKXmPGxUUS9xynCPR3+t/XyzkMciJiYsau3j52tXjkVz+cKk63RqANSQyA7cEYgLSLyLuAu4A9UtR+YBDwN/B0wC/gr4PsiclaWbcuaISqqN6K7qOI+uC6aVp+MfI7qours6UsmZIvrbMhHvR/r7OW5iIkZ/SBUUZbgtMnFywI7HH6/5smNBzjRHW2/Jo7S6DANNZW8a4mPeo8+K98S8wnZZUumJhMzRq31vuNgR3LFH9cxI0y2RqMDSN3qrwXSViQRkVuAVcDtqvofAKr6gqpeparPBS6uHwNPALFYbfiZw47WDjbuj1YLeHMR64JHxc/KV6/bHynqPVw/O64D0ZxJtZwTJGaMWmfbf1fzp9RRVREf5VSY65bPSCZm/OmmaCoqH4QZh7QomfB98JnNLZGi3ts6e9gXbDDHTWDiqaks58pkYsbR9cGayjLmTIqP5DsT2f5KNgKSckw42WWFiHwW+DrwPlX9t9Dxq0TkIymX1wDZJ93PAefMmZCUwEWZlQ+RBMa0Y8PgD7b5WBev7h7ZReVnrnGXBF4frDbWbGiKlJixFAbXqQ3VXLDA13ofeSAaGBgoCR/5VUHUe9TEjFtbBlOPx01GHMb3wRe3R4t6j1OyzChkazSeBKpF5HYRqRSRDwLTgTXhi0TkVuAPgUtU9cmUe/QBXxWRS0WkXERuAt4BfD/LtuWERGJoLeCRXFRxlgSGWTytntODqPdVr488EHm3QNwlgd4YHuno4YVtIydmHHTjxHdwhXDUe/OIUe+t7d0c6XAz9zgPrhPGVfLOZNT7yH3QS4gnjKtkasyUU2Eul2mMqyynP6KLKk7ZsKOQldFQ1S5gJXATcAi4HbhRVdtFZLWI3BFc+udAA/BKSizGGar6NPBx4B7gGPAp4AZV3ZtN23KJd1FtO9g+RJWSjjhLAsMkEomkiiqKiyrOqpUwcxtrOWu2T8w4/A+2v3+gZH6w717uot7bu/v42abhVVReiAHx/768MfzZ5hbaRnBRhTfB4yb5DjOuqpwrk4kZRzYam5rjF+E+HFnLKlT1deDiNMdXhv5/yQj3uBu4O9u25IsVcycyY3wNTcc6WfXGfmRG5i+3VAZXcMbwm09uYf/RTl7bc4QV8yZlvHZzCbg7PCvPmsEbe4/y6Pom/vp9yzNGeQ+pnx3jGTnAtPE1nH9aIy/tOMTqdU1cG9RwSIc3hDPG1zC+Jh7JMjNxzbKhUe/vO3d2xmuTwYox/67A9cFVb+zn+W0u6r0xg/S+p6+f7UGZ4VIYM8DSiESirCyRDB4bKTI37pLAMEtnNLBginNRDbdf093bz44SkgR6xdtIUe/eEJYlSH4Occb3wcc3DB/1PlitL/7f1cTaKi6OmJixVFyJAFfIYNT7o8O4qOJSZng0mNGIiF9Gbz5wPDnjSUfcJYFhEolEMg38w29kTsy4o7W9pCSB86fUsWymS8w4XGyDH1znT66juiIeWWCHw0tv27p6+fkwiRlLaVUI8J6gDz6tLbRnSMwYrp8dZ9GCp656MDHjw8NIin0frCovY16MkmUOhxmNiJw3b1KyUl0mX3kpSAJT8QPR3iMneGNv+sSMvmOXiiQQBmuiPLKuOWPUe1xzTmVixoQazgsSMw63XxP3OKFUrlnmot67hknMuPXAoHKqVIyhXxk+N0zUu++DC6fWxSpZ5nCURitjQFnZoIoqU0rnUpEEhlk+a3xyhpMpXbqfuZaKJBAGxQsHj3fx8o70LqpBV2JpDEIwmHDysQ3po94PtXdz8HhQZrgEVrvgRCMXBbXeM/22fB9sqK5g+vj4CkzCXHXGdKoqXNT7oxkkxYMTl9Lpg2Y0RoEfiN5sakvmKwpTKpLAMM5F5aPe07uoSmmfxnN6KDFjun2ogYGBQVdiiRh4GOyDxzp7eTZN1HtY8l0qG6sQSsz4ZkvaxIzhNO9xVk6FCSdmzLQXurmE3NkeMxqj4Pz5jUwJjEG6TlAqksBUvCtn16EO1u87OTHj5hKTBHpWhgpppUqK9x3tpD1IyVFKg+vsieM4Nxn1fnIf9PLNqQ3VTKyNX7LMTFy7fHoy6j1drfdSHFxh8Lf18y0HOXpiqKS4t6+fbS0+NU/pPJcZjVFQXpbg3WdmThFQSpLAMGfNnsCcSemj3ktREujxP9gDbV38YtfQqHf/XcWpfnZU/HM9uqH5pKj3Uok7SWVKfTXvWOBcVOlUVKXoSoTARVVeljbqfdehDrqD76+UxgwzGqPEyzk37D/GjoPtQ86VkiQwTNhFlRr1vrM1VD+7xAaixdMbkm1OHYj84DqvsZaayvgrp8L4FdSRjh6e3zo06r0U8p5l4vogXfqTbw6t9R6un10qe4We8TWVXLrYJ2Yc2gf9eFFZHt9kmekwozFKLljQmAzUCcs5w5LAUvzB+g3WHa0dvNk0KCn2Ee6lJAkM4/cAHklxUW2OcW2QkZjbWMvZc1zU+0kDUXPpbax6rls+nUQCOrqHuqjCZYZLRZUYxvfBn206OCQxo5+4LJhSR2WJKKfAjMaoqSgv47rlzkUVThEQlgSWYsc+d+5EZgUlXMOzcj8IlZIkMIx35ew/2smru48kjydn5CX4XcHgamPN+mZ6AxfH0Y4eDrTFr352VKY11HD+fJeYMWwM/eBaV1We7KOlxDVBYsbuvn6e3DgoKY5rOd6RKL1RIAZ4V84be4+yO1g2l6IkMEy4vO2qkIuqFCWBYWR6AwuDxIx+43hgYKAkFWFhvDE81N7Ni0HU+5aWwRViqX5f70kmZhx0UYVLoZaSwMQzobaSS9JEvce95kkmzGiMgQsXTmZircvp42dEpSgJTMUPRNta2pPPE/cqhCMxpJDWOicpbj7WRVuyfnZpDq6nTa5j+awg6j0YiHzes8l1VRlzHcUdX0jreCjqfVMyw3JpflcwuBf69KYWjnf10hdOllli+zRmNMZAZXkZ1y5zLqpVgYuqVCWBYVbMncSM8W75v+r1/U4SeLD0JIGp+MjcvUdO8Nqeo8lVIcDp00pnAzIVv+JdE9R631xCyTIzMX18DW9PRr07Y1iqg2uYa5ZNpzyUmHHP4XCZ4dIyhmY0xoj/wb62+wh7DneUrCQwTFlZIjnTW71uP7sPn0hGHZfyD3bZzPHMn+w28Ve/sT85uM6ZNI7aqvjVz46KFy8cPO4SMw7u05TudwWDv63HNjbT1tnDjtbSn7hMqqvi4tODqPdQHywvSzB/SmkJTMxojJGLT5/C+Bo34Dy4dl/JSgJT8T/YTc3Hk663UpMEphLer3l43f6SlqWGWRiOel+3PzkjL0UhRhg/cWnr7OXeF3YlywyX+nMlo971AK/vcaKM+ZNrSyJZZhgzGmOkqqKMa5a5zn33z7eXtCQwzHmnTWJqkJjxO89sB0pPEpgO71PefegEa9a7IKtS/65gcCB66LV97A+SZZayewpg1sRxrJjnot6//cw2wCXLjHOZ4ShcG7ioOnv6uffFXUBpeiZKeyQoMu85e1DBAqUrCQxTHkrM6J+rFDt2KmfOHs/cRjfo+Ocq9cEVBo3G4Y5B/f+p8H15FVX4u4pzmeEoTK6v5sKFTlKc/G2VoGciJ0ZDRFaIyEsi0i4ia0XkwgzXfUJE9orIMRG5V0TqQuduEpFtQRnYn4jI9Fy0LZ9csmgKDdWDPvFSlQSm4mMAPKfC4BpWUXlKVZYaZtG0+iH1JSbWVjKlvjSVU2G8i8pzKhhCODV+W1kbDRGpAR4C/hWYCHwTeEBEqlKuey+u/vcVwFygEfh8cO5s4C5crfGpQBPwrWzblm+qK8q5etmgbStlSWCYCxY0Dhl4SnE2lA6/r+EpxR9sOsIDUakly8zEnEm1nBNEvcOp811dt3wG4QVTKRrDXKw0rgD6VfVOVe1R1XuAVuCGlOtuAe5W1U2qehT4LPAhESkHbgYeVNUXVfUE8KfA+0SC0lcx5vrQQHSqDK7lZQmuC9WgLsWOnY5z5kxI+sVnTaihvrp0lVNh3nP2YB88VSYukPLbOkWMxtSGai5Y4FxUZQmSgaelRC6MxlJgQ8oxBZaPcJ0CE4DZqedUtRU4EhyPNZcunpIMpHrbvElFbk3ueN+5swHn7ig1SWAmEokE7zt3FgArTjt1vqvF0+o5Iyhv+7ZgA/lU4PqzZlJVXkZleYIzZ08Y+QUlwvuD39byWRNKLlkmQC6mWnVAR8qxDiB1pEm9zv9/7SjuETtqKsu5/7YL2X2oIzmDOBW4YEEj3/3gBUyuqyo5SeBwfOLqJSyeXp9M63AqkEgk+PZvn8cvdh7mhrNnFbs5OWNuYy3/cds76OsfYFaJK6fC/Mb5c2moqeTM2eOL3ZQxkQuj0QGkfqO1QGppu9TrvEE4PvMAJPMAAAxFSURBVIp7xBKZ0YDMOHXcAh5fdexUoqqijP+1Yk6xm5Fz5kyqLZn67aPh7fNPnYmYJ5FIDHEplhq5cE9tBCTlmHCyyyr1OgGOAvtSz4nIFNxG+cYctM8wDMPIEblYaTwJVIvI7TgF1C3AdGBNynX3AneJyA+B3Tjl1H2q2i8i9wM/FZF7gFeAvwFWB3sbhmEYRkzIeqWhql3ASpxc9hBwO3CjqraLyGoRuSO47iHgy8AqYBduo/tTwbm1wG3APcABYBZwa7ZtMwzDMHJLIlzasxQRkfnA9ieeeII5c049X7VhGEY+2LNnD1dddRXAAlXdEfV1lkbEMAzDiIwZDcMwDCMyZjQMwzCMyJjRMAzDMCJjRsMwDMOIjBkNwzAMIzJmNAzDMIzImNEwDMMwImNGwzAMw4iMGQ3DMAwjMmY0DMMwjMiY0TAMwzAiY0bDMAzDiIwZDcMwDCMyZjQMwzCMyJjRMAzDMCJjRsMwDMOITNY1wkXkE7iyrQ3A/wAfUdX2NNfNAf4RuBToAX4AfFJVu0QkARwDEqGXPKOqK7Ntn2EYhpE7sjIaIvJenMG4AmgG7gc+D/xxmsvvBdYBs4GJwI+BzwKfARYF1zSoamnXnzUMwziFydY9dQtwt6puUtWjOCPwIREpD18kIlVAO/AFVe1U1SbgPuDi4JIVwGtmMAzDMOLNiCsNEakA6tOc6geWAj8KHVNgAm41sSt5ULUbeE/K628AXgv+fwUwQUTWArOAnwEfV9W90R7DMAzDKARR3FOXA4+lOb4T6AU6Qsf8/9dmulmwf/ENnMH538HhLuB53EqlMzj/Q+DCCO0zDMMwCsSIRkNVH2foBnUSEXkdGBc65I3F8QzXjwP+HTgLuExVDwTv8Zcp1/0xcFBEZqrq/pHaaBiGYRSGbPc0NgIS+luAo8C+1AtFpBH4KdAIXKSq20Pn/kxE3ha6vCb4b2eW7TMMwzBySLaS23uBu0Tkh8BunHLqPlXtD18UuKQeAJqAX1HVnpT7LAWuE5Ffxbm8vg48qKqHs2yfYRiGkUOyMhqq+pCILABW4WS0q3ASXERkHrABWAbMAS7DrRwOiyQXJ79U1XcBf4Dbx3gTqAru8+Fs2mYYhmHknqyD+1T1m8A30xzfxaDqahcZ9kWCa48Bt2bbFsMwDCO/WBoRwzAMIzJmNAzDMIzImNEwDMMwImNGwzAMw4iMGQ3DMAwjMmY0DMMwjMiY0TAMwzAiY0bDMAzDiIwZDcMwDCMyZjQMwzCMyJjRMAzDMCJjRsMwDMOIjBkNwzAMIzJmNAzDMIzImNEwDMMwImNGwzAMw4iMGQ3DMAwjMllX7hORT+BKvDYA/wN8RFXb01x3PvACcCJ0+Euq+qWghviXgN8L2vQ94I9UtS/b9hmGYRi5IyujISLvxRmMK4Bm4H7g88Afp7n8XGC1qr43zbn/A7wHOBsYAH4CfAz4h2zaZxiGYeSWbN1TtwB3q+omVT0KfBb4kIiUp7l2BbB2mPt8XVX3q2oT8DfAh7Nsm2EYhpFjRlxpiEgFUJ/mVD+wFPhR6JgCE4DZwK6U61cAnSKyHSgHvg98WlW7gvtsSLnPMhFJqOrACE0sB2hqahrpUQzDMIyA0JiZbpKfkSjuqcuBx9Ic3wn0Ah2hY/7/a9Nc3wI8DfwzMB34AfBXwJ8BdWnuUwZUA50jtG8mwM033zzCZYZhGEYaZgJbo148otFQ1ceBRLpzIvI6MC50yBuL42nuc2Poz20i8iXc5vef4YxE6n16VXUkgwHwMnApsB+wjXPDMIxolOMMxsujeVG26qmNgIT+FuAosC98kYhMAu4APq+qbcHhGgZXEf4+L4buszFKAwL31s/H0njDMIy3OJFXGJ5sjca9wF0i8kNgN045dZ+q9qdcdxT4AFAmIn8GnAZ8GviX0H0+JSJPAj3AnwP/nmXbDMMwjByTlXpKVR8Cvgyswm18H8FJcBGReSJyXETmBUbkBuAc4CBuZfAD4BvBrb4FPAi8hNsQfxb4WjZtMwzDMHJPYmBgJHGSYRiGYTgsjYhhGIYRGTMahmEYRmTMaBiGYRiRMaNhGIZhRCbrLLelgIiswEWiLwc2Ax9V1RfSXBcpY2++EZF3Al/FpVc5CPydqv5zmutWAVcSCmpU1XQpX/KOiHwK+CLQHTq8UlWfSbnupuC6abgMAR9S1eZCtTPUjptxfSJMLfAdVf1wyrVF/5xF5ALgx6o6K/h7EnBP0K6jwF+p6t0ZXns18HVgAfBL3Ge+qQhtngP8Iy4YtwenoPxkEGuV+tpI/akAbc6YnTvNayONM/lqr4jMY2g6JnBZNbar6pI0rx3TZ3zKGw0RqQEewn0438ElR3xAROaranfoutFk7M1neyfhDNbtQRvOBR4Xka1BdH6Yc4FLVfWVQrYxA+cCd6jqVzJdICJnA3cB1wKv47IYfwv4lYK0MISq3gfcF2rbVbjYoM+nubxon3NQNuBWnAS9N3Tq27jMC9Nx2aFXi8jLqvp6yuunAw8ANwNrcDFQ9wPnFaHN9wLrcLnpJgI/xiU5/Uya24zYn3LJMG0eLjt3+PWRxpl8tldVdxHKEygiM4BfAH+Q4TZj+ozfCu6pK4B+Vb1TVXtU9R6gFRc3EmY0GXvzyWnAKlW9T1X7VfWXwFPAxeGLRGQabra+rsDty8RwWYw9NwMPquqLqnoC+FPgfcGzFA0RqQe+C3xMVfeknCv253wH8HHcYOTbVA+8H/icqnaq6kvAfwC3pXn9B4C1qvpQMHh9AVgoInkzGhnaXAW0A18I2tyEM9oXp79FpP6US05q8yjbEXWcyRWZ2hvmLuAHqvpIhvNj+ozfCkYjNYMuuCy6y0e4Lpyxt2Co6lpVvcX/Haw8LgVeS7l0BdAG/EREWkTkWRG5qIBNTSIitcAS4OMi0iQiG0Xkg2kuHfIZq2orLiB0aWFampE/Ad5Q1R+nOVfsz/ke3IwwnB9oMdCjqttCx9L1aTj5M+/DpY5Id22uOKnNqtqtqu8JjIXnBk7u16PpT3ltc8AK4BIR2S4iu0TkKyJSneb1UceZXJGpvQCIyJXAJaRfxWX1Gb8VjEZqBl2Cv1Mz8abLtEua6wqGiEzALXl/Efw3TA3wPG62MQe39F8dLEkLzXRcFP+dwDxcLZSvicjKlOuifhcFI5i1347LuJyOon7OQY2Z1AjcOob62CHz51jwzzxDm5OISEJEvokbaP8mzSVR+1POGKbNLbjf3pm4jN9XkL6vFPRzHukzxrkhv6qqJyWPDRjzZ3zK72lwcgZdcF9k6oeZLtMuaa4rCCKyAFfBcCvwG6n5vFT1QVzqFc+dIvIxXKe+v2ANdW3ZDlwWOvSMiPw7zoWyOnQ86ndRSN4P7My0YRmnzznEaD7HWH3mIjIOt3d0FnCZqh5IvWYU/SnvjJCdO0xsPmcRmYv7/H4r0zXZfMZvhZVGaiZegr9Tl5KRMvYWAhF5Gy7j7xrg/YH/P/WaXxWRX085HM4cXDBE5G1BIsqR2jLkMxaRKUAjETMa54kbcAXB0hKnzznEZqAyUMt40vVpOPkzLwcWZbg2r4hII/BT3Hd+UTBwpbsuan/KKyIySUT+XkQaIrQj6jhTCN4LPK2qLZkuyOYzfiusNJ4EqkXkdtzG0C24pdmalOuiZuzNK4Ha5RHc0vLLw1xaD/ytiKzDDSKfwM10Hs1/K0/iOPA5EdmCU+pcAfwmQ2cy4GbmPxWRe4BXcK6J1cHeRrG4ENcvMhGnzxkAVW0TkQeBvxGR23B+898Crk9z+Y+AL4vIB3Ar1z8H9gCvFqq9kFT7PAA0Ab+iqj3DXB61P+WbkbJzh4k6zhSCC3Eu1eEY82d8yq80Ag34SuAm4BDOf32jqraLyGoRuSO4LmPG3gLzIWAq8NkgS7D/90URuUtE7gra+2+4LMGPBG29EaexLnhcSaD5/3XgL3Cbxt8CblXVX6a0eS1O4XMPcACYhZMNFoVg1j0HV8ArfDyWn3MKtwGVOAPwQ+BTqvoigIjcISKrAYKN5/cBn8Opea4GPhChjHKuuQg3IF0DHA71658Fbb5ZRNYHbc7YnwrZ4JGyc6e0OeM4U8g2B8wnpU9D7j5jy3JrGIZhROaUX2kYhmEYucOMhmEYhhEZMxqGYRhGZMxoGIZhGJExo2EYhmFExoyGYRiGERkzGoZhGEZkzGgYhmEYkfn/vZ2+xbcy/xkAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "theta = [.5,0,1]*6\n", "time = np.arange(3*6)\n", "plt.plot(time,theta)\n", "plt.title('A Periodic Markov Chain')\n", "plt.ylim(-.5,1.5)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Remember, failing the aperiodicity assumption means that we **cannot use** our markov chain for any statistical inference, since we aren't guaranteed that it is capturing the true underlying posterior.\n", "\n", "#### Irreducibility\n", "\n", "A Markov Chain is irreducible if it is possible to go from any state to any other state (doesn't have to be in 1 step). \n", "\n", "This chain is reducible, since once we transition away from $\\theta=0.5$, we never go back:\n", "![](../site_pics/irreducible.png)\n", "\n", "The problem is again, if our chain is reducible, we can't rely on the ergodic theorem to allow us to use our constructed chain for inference. This one is much more difficult to visualize since something that looks reducible might re-occur if we persisted with a sufficiently long chain length.\n", "\n", "#### Positive Recurrence\n", "\n", "For positive recurrence, we need the time it takes for the chain to return to some value to be finite. \n", "\n", "Example. Suppose we are modeling the mean of a dataset as we have been doing. A candidate value that is included in our chain is -10. This value is very unlikely given our data:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "ExecuteTime": { "end_time": "2016-02-23T07:27:13.413459", "start_time": "2016-02-23T07:27:13.255741" } }, "outputs": [ { "data": { "text/plain": [ "array([9.89591756e-14, 3.16951052e-08, 2.63872403e-13, ...,\n", " 6.54483182e-13, 1.83032511e-13, 2.84562968e-12])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "norm(-10,sigma).pdf(data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We do however, need for MH to have some positive probability of returning back to this value even if it is very small. And MH does this because of the accept-reject criteria it employs- inferior values of our parameter vector are sometimes included. So while it may take a very long time for the chain to return to -100 it will eventually (probably an extremely long time)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Chain Pathologies\n", "\n", "Ergodicity ensures that we explore the entire parameter space and never get \"stuck\" in any given place. Sometimes however, we have a breakdown in the MCMC process where our chain is limited and can't fully explore $\\theta$. The following example shows a case where our chain is never allowed to go below -2. It is easy to visualize the problem simply by examining the trace plot." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "ExecuteTime": { "end_time": "2016-03-02T07:56:13.423137", "start_time": "2016-03-02T07:56:12.990807" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEfCAYAAABbIFHdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXmYHEXZwH+zu9kkm4NckHCGQyjuS4iA3KACKoeKiIoioIJcKvIpYpQbBEFABJH7FpGbKIGEKxBIOMKVo5IQEgi5E5JssvfOfH9Uz25PT09Pd0/39Mzu+3uePJmd6a56u7u63nqPqkplMhkEQRAEwUlN0gIIgiAIlYkoCEEQBMEVURCCIAiCK6IgBEEQBFdEQQiCIAiuiIIQBEEQXKlLWoDeiFKqD7AQGAJsobVe6uOcYcBvgOOA0cBqYApwhdb6TdtxGeAMrfU/SpDvJWCJ1vp7JZSxI7Ct1vrJCMtMAT8AXtRafxa2nN5GFG0iDpRSdwPba633KaGMg4EXgR201rMiEi2MHJsDVwGHAAOA94ALtdaTkpIpCsSCSIajgf7ACuC0Ygcrpb6AaXCHAP8H7AwcAzQDryqljoxYvm8BPy+xjHHAvhHIYudA4D5gUMTlCkJolFL9gfHAFpgB3BjgA+B5pdTOScpWKmJBJMNpmFHPEuBnSqmrtNadHsc/ACwCDtZat1rfzVNK/QAYDtyilNpWa90ehXBa61URFJOKoIxylCkIpXIYsAOwudZ6IYBS6izga8BJwG8TlK0kREGUGcsU/SpwNqCBnwHfAJ4scPwXMSOSY23KAQCtdUYpdQYwGOiw/bS1UmocxuJYB9wL/DarhJRSJwPnANtjrMjpwB+01v+zfn8Jyx1kM+GPxJjQCvgIuE5rfWcBmedj3GC/VUp9T2u9pfVTg1LqFuC7QF/gOYzrY6l13n7AxcDeQAPGDXeL1voamxwAM5VSF2utL3Kp+yXgbYyZfzyQAV5x1FNv1fNDYBjmOfxZa/2w9fts4H9a63Nt5e4PTAJ20lrPUErtC1wJfAlYBUwA/s9Wx0vAPGA7YFfgj1rr613kHQxcjRl5NmBGnmO11hOt3x8FDrDqXa6UGghMA+YCRwF3WectwAw8MsATwK+01mvyHo4p8yjgD8AuQAvwFHB+dmBgPb/HgUMxo+Kfaq3/o5T6OnCRdd5i67w/aK0brfP2Av4CfNGS41XgN1rrGW5yWNQqpa62ZK8F/g382lam57NyubZazLv1M2ArzCDsNuDPlkxLgBu01pdbx38HeAQ4UWv9L+u73wKnaq2385DbzjvAUVnlAF3vJpbMVYu4mMrPKZiG+ijwEqbB/sLj+L2s/ye7/ai1nqe1fldrbV8z5SzgMWAnzMt1HvAjAKXU0ZgX5hZgR2A/YBlwv2UqF+KvGPfWFzEd1D+VUlsVOHZvTOd+k/U5yzeB9ZhO9TiMArvakmtjjMLQGIW4M+bFvVop9SXr+r9tlXMApiMqxFlWPfth7vehmM48y90YN98pmM77VuBOpVTWrXYXcKJSyj6A+hEwxVIOuwIvYBTG7ta1bAxMVko12M452SprDJDXoVkxlf9Z1/otzL19EnjW5jb8GdAO3Gz9fTNGIfzI9syPwXTaB2OU7yGYe5eHUuo44BmMst0L+B7GFTjRcb1nAr+3ynxeKXUEpk3dZ8l7MrA/MF4plVJK1QBPYxTXHph73wejaLzYCzPoOIjuNvGQ7fe78X5WTq4FLgGuse7JnzAj+Ou11mnr2r9mO/5rmPfxMNt33/Ahdxda60XZwVUWpdQJwDZWfVWLWBBlxHqJfgK8YBtp/gs413IRzXE5LTsC+TxAVbdpre+wPv9dKXUupqO+CzPa/anW+m7r9/lKqeswnfNmgJsMABdprZ+zZP4NZkS3D/Cx80BrpNsJrNdaL7f99L7W+jfW57lKqUcwnQxAP+BS4FqtdYdVz1iMUtpVaz1FKZV1fa3QWq/zuP55WuvzsuIopR7CKIlsPOdE4ACt9avWMR8ppbYGfofpgO6xZPkaME4p1RdjjWRdBecDr2mtx2YrtEaiK6zj7rG+nmV7Dm4ciulIN7eNPv+slNrTuu7/aa1XWhbfeKXU7cD3gcMc93U9cELWYlBKnWnJvYvW+gNHnb8DntVaX2i7PycC72IUeLZjnKC1Hme7vguB+7TWN1pfzVVK/RCYgYkNvQ+MtO7BAq11uyX3F5RSNVbn7MYK4Pta6/VWPb+wrnUHjGIs9qy6sKyxX2Csmrtscg4FrlVKXYJRwP9RSg3WWq/FWPOP0d0+hmMU5vkF5C2KUuoQ4A7g8WySRrUiCqK8fAXjernY9t39wC+BM4Bfu5yT7QiGA0WznSxmO/5ehQmKo7V+VSm1zHrhFWaUs7t1XK1HmfYMkazrot6nPH7k+lgpdRtwulJqF0uuXX3IVUxWMPJmZd3D+n+8ld2TpQ7oq5Tqr7VepJR6FqMEx2FG6PXAv6xj9wS2U0o5lVQdxirLUkjZZtkzK6/ljshSj21AoLV+Xil1E8Z18met9cuOct5yuJNes/7fDeOysrMrue0PrfV7SqnV1vFZBeGUfU9gjFLKLQttR631y0qpK4ELgF9YLrZngQc8lAPAtKxysJhikzN7XsFn5ShrB4zV8orj+5cxbWgnzECoHThUKaUxlt9lwDSl1GiMsltmkyMHl2d+pD1TyVKatwMTMcqtqhEFUV6yGUu3WZ2hnZOVUhdqrZsd32ddS/tifMs5KKUOxIw2z9ZaZ0fzbgHvlHX8CRil9G9gKsaEH0iBGIiNVpfvggaNveTaAfNiz8B0LOOAt4Aw6axesmbdqodhRq+Fzr0DeEApNQgTaHzUGnFmy3gE+KPL+attn53P0kmNVd/uLr913SvLr74nJs50mFKqjyMhwZmckFWoXokPTmrJvW9O2WuAG4B/upy7HEBrfaFS6mZMbORQ4ArgPKXUfh6p3E4Zs8+nju64WrFnVYxsma1a62al1PPAEZjB2hta63etuMuh1vdPOFy2dpzPqqt9KqUuxrSJO4DTs5ZwNSMKokwopUZgfKkPY0Ysdo6xvjsRyAn8Wj7v14ALlFL/1Vq32cqswfiJdwY+9SnKhZhR3cm2cs6xPkaZJRR0HfkzgEZMplbGkitrQWTlimJt+uyIelOt9RvZL5VSv8bM2zjD+uoZS54fY1xNX3WUsRPGlZW2zt8Qo2yvwcSW/MrSF9hAa/22TZarMJ3fn6yvLsR0TIdiAsOXYlwsWXZXStXb2saXrf/fcqnzfcwo+SpbfXtgUoenF5F1B631XNt5O1nlXGC5dy4ALtVa34YZBG2DiUkchBmQuLGbUqrWlsV3oPX/NNsxBZ+Vw/KaiVGWBwJv2L4/CKNsshbsk8BYYBQmuQDr/6Mwz/k7hW6C/frtWM/st5hkhEsLnV9tiIIoHz/GuA7+rLX+0P6DUuojjJvpFzgUhMVPMZ3OS0qpyzAvwmaYiXMHA18PMFr5BNhHKTUGM/I7HNPhgOmsoqIR2FYptanPSW2fYMz9o5VS0zAd8HUOuRqt/3dXSi0ulKXjhaVwnwJusjqXaZh7cBUmQyd7XLtS6j7gcozytbt1/oLJ0Lndit/UW7JuR75Lx4vxmAyYB5VSZ2M60+9jLMKTAKwA/VjgXK31JKXUr4A7lFLPaq1fssrZxPruCkzW0U3AIwViWlcBjymlLsdkt20G/A2jOJ71kPVK67wrMTGW4cA/MKNzbf1/LLCF5b5chwmwt2GyygoxEmOpXWFdx02YEfwM6/qLPqssWuu1lgXzB6XUMswz2s869nat9Urr0KcxiRobYyVJYBTEg8Ba/Ct4LBm/glEON2EU4yjbz81h2mmlIFlM5eNU4FWt9TTnD5Zb6Sbgi1bH7fx9JibIPA3zMs/A+MM7gH2zKZE+OQuYj3kh3sYorpMxLoUvBSinGNdiRrzvO7JjCvE3TId1O0YBXoMJQr5ik+sDjI88G0QOy/cwc0tusOo6D5MGfIXjuDswKcR32V0OWuupmHjSNhg33QRMnOMQWydUFGvU/FVMNtT9mBH8d4GTtNYPKJPS+gDmHtxinXM3piO/zwq+glEyKzB+83sxGXI/KlDn45hssCMxSuFBq/5D7dapx3lfwQS0H7f+P1xr3a5NCvbXMIPOF4EPMSP5r2utP/K4Df8DVmI68wcxFtIPbL/7fVZZzsOktP4J8578EaNQzrJdy3LgdYyVNtX6eiLGUh2ng88n+rH1/1mY9F/7v1sLnVQNpGRHOUGoXlQEy1UIQiHEghAEQRBcEQUhCIIguCIuJkEQBMGVHpHFZM103RsTFAqS+y0IgtCbqcVkc72pHWu9QQ9REBjlUNXrrguCICTIAZhMshx6ioJYDPDAAw8watSoYscKgiAIwJIlS/jBD34AVh/qpKcoiE6AUaNGsdlmmyUtiyAIQrXh6pqXLCZBEATBFVEQgiAIgiuiIARBEARXREEIgiAIroiCEARBEFwRBSEIgiC4IgpCEARBcEUUhCAIguCKKAhBEATBlYqaSa2U+i5wMbA5sAC4UGv9RLJSCYIg9E4qxoJQSm0H3AWcqrUeCJwLPKyUGpGsZIIgCL2TilEQWuvZwEit9WRrL96NMZvUF9wnVxAEQYiPinIxaa3XKaW2AuZiNhA/Q2u9NmGxBEEQeiUVY0HY+BToBxwOXKuUOjRheQRBEHolFWVBAGitO6yPLyilHgWOBV5IUCRBEIReScVYEEqpo5RSExxf1wOrk5BHEASht1NJFsQ7wF5KqZOAB4AjgKOALyUqlSAIQi+lYiwIrfUS4JuY9NbVwCXAsVrrWYkKJgiC0EupJAsCrfUkYK+k5RAEQRAqyIIQBEEQKgtREIIgCIIroiAEoUJYsa418DmvzF7OrCUyl1SIB1EQglAB3PTCHPa6bAJ3vvqx73PenL+KH905lSOun0RnOhOjdEJvRRSEIFQAf3luNgCXPDPD9zmvzF7e9bmtI+3rnMaWdlaGsFSE3klFZTEJQm/jzlc/Zn1rR/EDXUgFPL6tI80XL51AW2eaqRcexkaD+oWqV+g99HoLoqW9k/ten8/cZY1JiyL0Mj5d1cQlz8zg2udnhzo/qFPps9XNtHUaS2P8h0tC1enG/BXrOeL6V7jrNf/uMaE66PUK4q/Pz2bsk9M5/LpXkhZF6GV83hTdSvYpH+ZEUIvDL798+F1mLWnk4qf9u8eE6qDXK4gn3v0saREEoWQyCcaol61tSa5yIVZ6vYJI8sUSejeltr24LAJByNLrFYRkBxanM51hTVN70mL0OKJsepmApUmzF/zQ6xVERkyIovzw9insdslz6CWNTJ67gv3//AJPimuuZKJse36K8hOnCFV3PMUKFUCvT3NNi4IoyuvzVgJwzXjNhJlLATj3X+9yzO6bJimWYCNoK5ZmL/ih1ysIeU96Hn9/cS6frW7msmN2pqam8jz1mUyGu16bz6LVzaUVZDMJ/FgjKYId7xdRNj2XXq8g0hKE6FGsWt/GNeM1AF/eZgRf33XjhCXKZ/z0JYFmTPvBTyuOy8Uk9FwkBiH6ITTL1rawrLGyUhxb2ju7Pi+t0PTLF2Yti7zMJNtx0AC5UD30egUhMYjwjLliImMun8iaZslwCsJnpbqWhKpm7rJ1fLqqKWkxfCEKQvRDyby/cHXSIlQV7Z0xNLqARUaaYivvkG8+W93M4de9zAFXvxh6Da5y0usVhJjH3jw/Y6ntL/d7JR1E8kg7rg7e+Ghl1+dPqsCK6PUKohotiOdnLOW65zQdnf6WeA7L/BXr+em9b8VaR1AymQxzl61LPLmgqa2yRn9BlXSUSl2C3/6psfW41eDe7vUKohonyv303re48YW53PfGgljrmbdina/jynkH//bCXA6/7mUufOKDMtaay+2T5rHzn8bzxLRwkwWj6k/t5STZiqvwFUqMmpzU5AQF8UlFKQil1P5KqSlKqTVKqY+UUj+Pu844B6Ir17Vy+6R5pee7F+D9hWtiKTdLTQUODa+zlsZ+aOqniclw2biZpDNmFdNKwdc8CNvjrIK+qUeSsj0EsSACoJQaCjwF3AgMBY4HrlRKHR5nvXFaEKff/zaXjZvJt26e3PXdRU9N53v/fD0nHdMvzW3BzykFvwoiKSvsupD7KCRNHHo3UQsixDmrm9pCvQPVjv3RV4N7u2IUBDAaGKe1fkBrndZavwO8COwXZ6VxPqQ3538OwBIrH7+lvZO7J8/njXmruD+ge+i2V+ax80XjGff+4q7vsh3zR8vXcc34WSxZE23ef76CcO/Zsrdw5uK13P/GgqLbX0YVO7lx4pxIyukJ+NHRlTJgXbymmd0veZ5D//JSVbp4S6GmyiyIiplJrbV+Fzgp+7dlURwA3JuYUBFjbxBrA84duPy/MwE488F3ur6bt2I9mUyGo//2KuvbOnlx1nL+e+4B0QhLbkDND0feMAkwo8OzDt3W9ZjbJ83j2udm84+TvshB221Yqoi9mlyXUXKdTdB+7qEpnwCwaE0Lze2dNNRXTDcUO/aVX6pBOVaSBdGFUmoD4Gngbev/HkEq4hX831+4hrsnz2e95XqasXhtpOXX+vWFONr5s9MLb2d52biZNLd38uM7p5YgWXXjtx1cP2E2Fz013V9HEjiLKUGFYvsc9TtR6eTGIBIUxCcVpyCUUlsBk4FVwLe01vHmclY5cW7zWJu30F0VtOgewqermrh+whzunjyfSXNWRFJmfDoh4F4UtsNLjce0tHcyee4KWjuqI55hf6WSTtX2Q0UpCKXUnsAUYDxwrNZa1iRIEL8rocokrehpbOmeZ7FiXWvR4/08gZ74nM5+aBrfv30Kf3j8w6RF8YXdgqiGp1ExCkIpNRJ4FrhWa/3rnmg5VNsLGjbNNUm3QZJ3uKMzzWXPzODRtxd6HufntgbeIS7BIHXgSXq2a/to+TpOumMKz3m4Jb3IzvR/pMg9rxRyLIgqiEFUUnToVGBDYKxSaqzt+xu01hcmJFOkVEF7yMHvVgqVel3lnsbx0NRPuP3VjwH49hc3i6zcpBRKXNjrPuvBaXy8Yj2T5qxg/lVfT06oMlFtE+UqRkFora8ArkhajrBc+5xm2dpWrvzWLgVdM5kCn4XoSDL4qpc2RlZWqUtntHZ08trcFey15TAG9+tjjolItij5eMX6pEUoK6kqsyAqxsVUzSxe08zfXpjLw299ynM5i9vlUu7O65OVTSxvLO6/LoRfcVeua+PxaZVh4if5zvmNOQa1bAq57HJ2iHP8dsW4mZxy91ucdPuU7mNiujmBtzu1fa6rwB3/4qSmyrKYKsaCqGbsy/auWt9W8LhytoeFnzdx4DUvAqAvO4K+dbWBy/Ar7/89+n7gsnsi9g44k8nkBCTLWTfAPa+biZjv2ZZjKVf7W9bYQltHms2GNrj+bhe1f31tTkC+p1MTcJtYL1asa2X4gPpY25lYEBHQ1tH9oL0GROUc3b5o27Vs6Zp8K8JP43SawKXI//pHK7n62VnhC/BJohaELa2i3HKUEiiOkua2TsZcPpH9//yirx39+vUJPnCpZnInyoUv55n3F7HXZRP4wxPxZm+JgiiRZY0tHHXjpK6/PZV5OTsNx5T+NU3dM7db2js58oZJ/PjOqZ6KIspO7sTb3uDmlz5ylJ9h8twVLItwa9AoOr6W9k5enLUs8IYudoXqJUXQLK+oBoj25xmXAluwqjumUGj+hv0ZNdR3K4h0OhOJGyyTyVTuHIOIYhBnPTgNgAemfMKRN0zizfmrSpXMFVEQJXLDhNz1gLxe/qTSXC99Zga7XfIcE6z4yFPvLWLWkkZenr28yKYlufJGbck++e4ivn/7FPa5cmJkZUbR8f3u0ff5yd1vcvr9bwc6z94nRRmA/NNT0/nKdS/nuS9TLqPRto40d1qZVPnEFIPIca35Ob77c3+bBXHzS3PZ+/IJvL0gfGeXTmc47ubJHP7Xlyty8lwcMYiZi9dy/D9ej6YwB71eQRyx06iSzl/r8J8W6kS/fuMkWtrLN7XDbspOtNxNp1mb/9gX0+vwaKVxu0ketNbkKeVFmfrxqpxVbqMQ+Yl3FwGFR8CFCNpRepfV/Xl1Uztzlq3j8nEzCx5jjmvjtknzuOQZ99n1ORaETzmWN7aGHtX7GU/U13V3QX95bjYr1rVx8p1vsqyxhVPufpP/2OY3ZDeL6vRoMHppI+9+upp5y9czfnrhhBG/OJXMjEVrOeehacxYFG5Zm2pbrK/XK4ijd98EMKZuc1snf31+NlPmrSxyVjfNbU4F4f5aTF+0lllLuhtV3G2jkCVjgqf2vwuXEZWIhTqY/KU8gvPdW1/nJ3d3r+uUZJprZ46LqbAcYS0xr7W2npuxhN0veZ5rxutwhbvw7IeL2fvyCZz/H+8khJz0bcdlT/vk87xMOvszcntfGls7GPvEh7wwaxm/eeS9ru/vePVjDr/uZf74ZGG/u73T7UyXNiCbNGc5u1z0HNdP6F5W/qgbJ/HUe4ty3MpByH33REFUPNnnlcmYxdFumDiHE/75hu/z17fmjjC83v1yjhgKdUJX/HemQ3kUlikqP25rgeW/62qj8Vm9Ma/bJRH1Hf5sdTNjn/iQDz8rvjmT/XbF8aiXN7ZYZWdYsqYl5xlf5rAu3LCL5GcQdKbl5/5PgFnKdsX40uzlHHfzZPa+fELuMTZB3JpAKgUfLc+fH5G9xgcsy9PPEiSlcNIdU2nrSHP9hHiWla/UMIkdURBWA82Q4dW5wRdFa3JseuK1RHaJA5pAFOp6b5v0cewWhFM5FdroqE9t9M0v6o751Lvf5L43FvCNv70K4OnesC8XUbKLyeXuZ8u86n+z2OfKiYG2nG3vTHPx09O7/n5RLy+aZVSq+n76vUVFj3FbzsXvKsL/fGVeYJmSpliiwLrWDh59e2Hsys8vvV5BZF+DTMb75S9EkyPTxStIXQkWRBCc4oYVv7nAzmF2F1N05na093jWktzZ0Vf8t/BI3W4plZqQ4HU7brU6xiCTIB+c8gmvzc21GrwTFILsKOhbDHO8vQ4XN2Mq5a89uG1MVWlemwkzlnLSHVP4aHn+/u5u/cFv//M+5z3yHicG8GLESa9XEN0WRP4DW7a2hZ/e+5bnSMjpPvF6p8ppUnopKvuObl4iRZV1VUhB2GfRhlHOdm592aTQltpB2NOB3bijYIZQLl5y+JnY5HZ6KUp/zrL8ZUCKFRfX/Cv7vXELQ/md+BX1gKu5rZP7Xp/PPJfOPCyn3fsWk+as4Cd3vQnkKj43+cd9YHaMnLMsOhlKQRRE9oOLBXHBYx/w/IylnP3QtEjqyslyiTnl1fMd85tJEZGI7QW2GK2zuZjaStyG9Mr/mUl4pYg8Zd5KdrvkuYK/n/fv9wr+5qTCBrKhFKfv/aJKyN5ys1JqUv7uX7FBRdD5Jlf9byZjn5zOode+7Ov4TMb/vA03a63SrB03REFYDTRDJm+EH2bxNS+z3F7+31/8qOBxUeA5CvP5Qkdl8RSKvfSxDR+P/ftrkQTFS3npznvEWwE8+o7/YK1T8ba0d7JgZakL00U7pC+mAMIu9+7GH5/8kGmfmD3a7YMjt0w2vzGIqC3yx975LNDxR1w/iR8VmWzqxH6kpLlWAfYspo4QUWRnW/Z2MflvEO99uppfPPA2M0NuI+r5ivlccjislZMt/e0Fqzj+H5N5rUDw3945zF66jg8XFc8UKkYSkxEzmUxOOqb5LveY426ezEHXvMTEmUuLdvMdnenIUyDdS/OWxHcMwscx976+gONunmyOz3ExuVkQPhVEEQ0RtC3Yq13nYxa9XtrIpDkreOydz3zHg4LMZn/4zU98lRknoiDsMYgIsozCBKnXtXbkdQjH/P01/vvBEo79+2vh5PB4x3KTXOObKPftW17nzfmfc3mBwG6dI4up1DgElCZz2HNnLF6bnwrqKCur6E+9562i5RVKCy5lQO92bcXKC1Pd3ZPnBzrezYJIpfCldTpdLqqUe2Sf9Bpk3bDzHnmP/a4KvhpAsQHjbx/9IHCZUSMKIqsgMplITL6gi/VN/XgVu138XMFFtwp1FoXIrh/k33/s8Zvj76iDln0imgdhJwmr3S2Nt5SJci0FgvqlEfzG+H7etqL9zJmwD4ZcLQjHS3TaPW/lLTPy5LufRbreUmNLbnLC1I+DLffR3ulPFnu7CNNWZy1ZyzPvF08fjgpRENY4KZ2BxWty88L9vCDOQ4K6mM588B0605muyT+lsvflE2jt6Ixk28+4faTO0WOpyxZnMpmK2da1lL6rtSPtnsUUvkhXimcxxZPGlJPm6lJFTSqVc8yEmUu567Xc7LFz//VupO3T+e6XgzDyH3H9pK6F+spBr1cQpbx1974+n/krvXPJ7ZQjzbWprRO9pDESC8LZS0WtL5ybxZTqd+9MZwLJmE5nYhqtl3YtQa3GYphsm/zviymAcuzlUyiLyYlbTMDnoD0UUQboc7DJLDOpq4BSGsIfn5ye951nVlCZWkSxa8qZSU2GxWua+dOT+ctJxD0ar3WZdr6mqZ2/vziXOSEyyILe3hNve4Mxl08oade9Qnjq3SJyFlJaYZtqoftSrLi4Okn79btV4VbvsIb6vO/Cvk8dPlKqvVZEKIUwWUyTP1rBQ1M/SWTtpl6/o1zUr4DXIyzvTGr/s2B/ft/bvL9wDfe8viBn4/jQ4vqsOz8DLMVvH32fZ6cv4ZrxOvAm9ukCI2XXY9MZplh+5tsmzeP3R+0QqC47bpdbyvSSlvbOSK21jnQ61nkQQbEPPNzeiZpUKq8zdNMFYZIa3pq/ih/dOZXT9t+KX39VdcvkKMrpon1o6iecOGYLX3X88ckPaWzpYKPBfQPL58b3bzPbxg7qV/7uWhRExC+BlxJwa89x6IyaVMq3eyADvL/QPb00b6mN0sTKw03EZ23rGQUl7YhBeN2CnJGc9WCiHKFlMhme/XAJkz9awTmHbZv3mxemDblk94QczqTT7tZg0Swm2wGfrW5m0yH9u/7uTGe6Ykil3DW3d8Kt7bplLBUbcLn9fOo9b9HU1smNL8zNURBOPvhsjbXkh/n7gsc+YOsRAzzrA5iztJF7Xy+8RlbQJma3pCfOXOZxZDxUpItJKTVGKVXefN1KAAAgAElEQVSWUH0UwVw7XoOaYg36rtc+5o0AS40XoqbG+7pyNrv3VGi5v1X68sRBYhDRKoP872YvXcfp97/Nva8v4PT7cjceKtYOjAfE/Zgwcrt1rlC87ds76j/asuzue2MBu140nokzw+23UGwuQE1NKu/q3a47jEUeZBMhZ/H/tZbB8KKpwMKUruX7OOZ/H3YPmJJYwK+iFIRSKqWUOgV4Dsh3OsZA1BaE9xae3k3i4qdn8L0IFulata6NJ9/1Nys0QIw68qBaFPMe7JiRsj8yBT5HxVTbFpBvLfg857di820KdXxL1rbkbVDlh0KKs/g8iO4DPm/qTjMd+8SHrG/r7JrTEVRp5dx7n3MZ3Kpwaz/FlF4pA0I/btticZugcT37OmbL1vZyBQH8HjgXuLxcFfptLs9+aEYP7Z1pFq9pLnhcOZauKMb3b5/CczP8je48feWO316evbwEqfK51bFcc6m6Oh1gbRw72VOiDLlsskG/gscXG/l6/X7BY96b97iWF7Lh2S2IKDZ3cqNQDMKJ22KDYS4rrriK3/KDNk/7HJum9uCDg1KpNAVxJ7A78GbZavTZYE6//x3WNLfzw9unsO+VL/BKgc7S6+WOesQcDUFsiHgptbbOjP/xWY6bo4SaC+2pYN/rYsTAXGO4WG3pdOGO5L8fBI/RFLovc4usGGofMXuNnoPevdx7n48JUud+57Z9aBgXU2zpqxhFHHX59oy2zjjzegtQUQpCa71Ya13WuxDE5Gxsae/KfHHuD5yllBhEufDaMGjO0sauUUuFiOubdCbD3a/N93Wss8vs6EyHmix11oPvuH5v9/tvN3JQzm/FRvRRt5O7XvvYdYbzLx9+1/M8ezuJ0oAotvqr37rmLM1VcCffNbWolVvKZRTr+1s70jmuODdyLtfHc7a7mArFkuJEsphCtphC53nHIMLVFScn3pYb8/jKX19hmw0HMPG8g8s+kafUPmjx6hae8rGLGeQHSn//eLh1b96c/zluktuVQIfjRhYNUmcyTF8UbpFGN7xWDr7vjQW0tndy2gFb5/1mHw17uZhKadcF01wpPldhicN6e0kv5yVdxA0ao4vp2L+/FmgF6OyVt7R3sryxlc2HNeQdY7cgOhKwIERBRFye18tSaNG6JHFbQya7H3ClLFtx3+vzeX3eSq7+zm6exz0+LdhyzXb+/Zb/pbzzccuw6f7sdC0WU7yZTIY/PZU/CTMOxlrZSXtsMZQvjh6a81tNjgUR3ZtSTKHE6QYqhWLeBj/Kwc16OvqmV5m9dB3/OX3fvONb2rsVpVgQCeDlW/VqEIWeVXD3QP7xcaeT+l6DLawYEcs/1pqxPnr4gJzcdCcbDopmYlIUXDO+ezXQtx1ZTEXnQZRx7/Is81esz1MQdqIMUhedKBdRXa6ZWz5kKkQk2/jaPqczGdo60sy2XGU3TJyTd3yzxCCSJexDL7R3xNsLPmf+itI2hqmUWHZYMaJaS6i1o5PHp3WP7Bd+3uw5ugy0cUuOiyn6G/65x9alxZ5vEiPFFpf5AXY5Pe97gJby8Juf5FhrhTrxuG5BVMonCi5+egb7//mFrr/ddl60ryrrdFWWg4pUEFrrl7TWI8pRV9jmknXDOHnk7YUc/JeXCm6z6QenS+LXRYKJcRG24wyrIJx90A0T5vCrh7s34qmr8Tby//LcbN915Sy77HFcHP1J8RnA5e8IWi1XxoxFa/n2LZN59sMlOXJGdR+cexy4WxDR1PX0+4vy1l0qKUhdmjgGx+Uus60DVizGkEQWZEUqiHISxIJ4PMCWhM7169247/X5rFiXf5zzpXmsBN96KYTtp/wsSeCHB6fmLoFeW5OKLI+92GzeLH4mRz0SMH5R7D1PwoLMWhA/uP0N3l7wOaff/3ZZkircrrXWJc01DC/p5XnLXpSyhPkdjiXHo6bYoDLMjpel0usVRKFxwaQ5+dkQ1z7vf4T6vI+JamNdVoOF+LOd/M82LvM8CEd1zifzn7cXRrY0it8r8zNy/tebnwaru+hSGwm4mCwLwu4aK4scbi6mCIPUzsSFUkqO4r30eqfeK7AmWpYkBg69XkEUaosn31XaXL1+fWpDnxv3fIkLHvOX0hlWjLDy+zorBpeP9+5v0VdYykzquFjvsu1tbiyknGmu4csrhv1xTo9gD/SgVGKquxeiIAJ+75eBfStXQfglvIIIt22ms4Ny65yj6jvsdb04q3DufBx9VXEXU/mf/x2vfpw3KCpHLMSthigzpvKbUPcXX7/x1cjq8cNP732L9z5dXdY6S0UUREw51/alFoJS7VlM6UyG7/xjcuDzrnO48NyeTNhMjs50hmmffE6bFUC3l/LZ6sJra8WRk1/UgkggzRXy19ry62Iqpbm6LtYXoVpesLKJY256lfveMLEI5+O8ZvysgjvuRc3zM5Zy4wtz468oQkRBxFRuKSlplbKsdikj2Q8/Cz4TeNKcFTl/u/XNYf3iVz87i+NunsyZ1tIYfi8tDndHsborxYK03+tUyqQdn1lgaZGwuD1OM9clmnuwprmd9xau6ZoQ6Hycf3/xIybOXFZ1rp9yIQoiJg1RSoCvUiyIsEPD6Dq46B5OduVYP8kDORL0khiEG04x/vHSPMa977InQgniuimCOGdSuxW9rLG1YlYNqDREQXhtrFNCOy3FgvCTIlsOwr40UblIYl1xIcH+oNjgoVIGCE5F9fq8FQWODE/cMQgnbu/74P51YkEUQBRETG3Rz8bohTj8upcjlCQ85c5iihovKZIcMVaNi6kMcrhdayoVn/520z2D+/URBVGAXq8g4iKJafFRk/RSTKXq7kJyjHt/caIdQvG1mCqj7Xxtp1Fdn7fZcGBBa7sUZetmbcblYvrr87NdXYZ9amuq3sUU10C31ysIrxu7YGVT6HKTWJo3asKOZKMaAZfa6AtJceaD7yTaHVTiTGo3BvbtXsuztsb9eXxSwjsChV1McSjwGybOcd3XOUN5spji5MrjdomlXFEQBUZFpVoAnUnlKkbEg1M+SdzFtLTEPXi95mIU2gmuHBRz3VTKzoN2KdIZdwXxypzlJXWuMxfnZ7ulyN/rISpc1wnLJBqSioS4rC5REHHFICrkJQ/L7x//gGUhX9JKufRrxmv+5rKEMsBJd0wtszTdFHMxrW4uvBJsOXGuV+XWCa1rjX6f5HK/OxkqJ7U8NOJiiodKTHOtFFYV2T6xEDNcRoVJUWj9LDdXQ7ko1hfdWECpJUmhDnRdS0fko+/mELPwSyGTqZxBTVjiyvsSBRHTra12CwLCx1HaItoPohKIY4ScxH4P4eiWM53JuAZ4Y7EgSsgADIMJUFfLM3FHXEwxUYlprpVCT1BybvQvYSHFKKgW69LpYnJ7VVpdNhoqFbdtcOMkk4F/TQ22Im+lIVlMMRGXafbs9CUxlVw+qqUjK8aYyyfk/J30pmKNLdGPuuPgRb2s63M6437f4hi5ltsCzWA2+qpmREHERFw3NsxaRJVGT7Eg7Lt2Cf6xZ5FlcHcxxaEg4rBKvKj6ADXiYoqRytmjttKo9lTdQlRPDKByKORiMnMWor2fUe1p7hdpDYXp9Qoi1vV+qpxy+4LLRXb3NME/Jkid/30cI9fFa8o7R0UsCI9yYym1ihD9UJieEoMQSieTcV/ZtibGdZPKRQ/IJ4ltoFtX/JDyoZTaA7gV2AmYA5yutX4jzjrj2jCoJ9BTYhBC6dz3xgK2HzUo7/s4V14tF5WyOGIpxJWu78uCUEqNUEptr5TaOBYpTB39gKeBu4AhwI3AY0qp+rjqBKgVBVGQnhqDEMIxa0lj3nc1Ma2bVE56hospnnI9LQil1HbAPcAY23frgNeAh4EHtdZRrQtwCJDWWt9i/X2nUupXwDeBRyOqI48BJewd3dMRF1P18dUdR/JcwE2RSqEHGBBVP4sakktzvQfoAL4BfAnTifcBWoFrgNlKqf0ikmV7YIbjO41xN8XGwH4V5WWrKHrCirS9jXIbxD1g8N0jXExxRVOL9Y67A7tprbsWtFFKdQLnAZ8AZwFPK6WO0Fq/WaIsAwDn2sFNQEOJ5XrSt04siEK094ShVRXRpzZVcuZYnNt1uvHBZ2vKWl8c9ARLOS5LrpgFMQ/Yyu0HrXWH1vp64HfAXyKQpQno7/iuAVgXQdlCCCQGUV7qakpPKiy3gpg0J/ptSMvNuf96N2kRSiauZJtiLfJqTCzgMI9jxgN7RyDLTEA5vlPku50SZcvhsRo0FYW4mMpLXQTDwJqeEBQQApOIBaG1vge4HnhWKfWWUuoC65waAKVUHXAGsCgCWV4A+iqlzlZK9VFKnQKMxCigWDl8h418H3v6QdvEKEllIUtUlJfGCFZGFfXQO0lsLSat9TWYWMQHwG8xbqCZSqkVQCNwJvCrUgXRWrcCRwInAquAs4GjtdbrSy27GLeetJfvY4/cObZM34pj1fpw+0EIySEGRM9iv22G+zournkQvlJ4tNbTgZ8opU4FdsVkHA0BlgMvaK0/j0IYrfX7QFRZUb7xO9lnUL86NmjoE7M0ghCecscghHjx2zdVxExqrXUaeNf61+vYcGDfpEUQBG9EP/RKkgpSCzYkAChUOmJB9Cz8Pk/ZcrQCEP0gVDqydEzP4ecHbe3bxSSruVYAMjoTKh1poj2HbUYM9G0ZyI5yFYCs/CpUOtJGew5mDw5xMVUNtXK3hAqnJ7pBD9h2RGLXdcROo9hmwwGJ1N2ezvi+bglSVwDiYhIqnZ7YRI/bY9NIliEJwxbDG9hqxMDY6/n2npvlfdfRmfYfpBYXU/KI+S5UOj11EJPUxkSn7r9VWZRu3z75XXFHZwa/elGC1BVATzTfhZ5FT1QQqVQ061SFqXfk4H5lee/dduvrSAeIQYgFkTzlTCH82k4jy1aX0HPokQqCFLW15b+uDfqbVRPKcU9rUimGDcjdPDOQiykOoRAFEQivh7XJBv147BfRrRLy4323jKwsL76yoyiiJNhsqHNl+2jogfoBgN8dsT3gPtL+2YFbx1Jn9n0vxz11q0OC1FWG1zM4ePuN2HRIdC99uWZtX3/C7mWpZ4th+cuk/+OHe/KFjeIPAFYiu20+JJZy42g2A+pr2W5kss/phL035+mz9ufRM/IHYXG9KdlyyxF7rEml8vbG7uhMyzyIasLLgshkvBtq/z7Bdq4rl8+1XMG/vUYPzfvuiJ03ZsKvD2J0L9pjI0tsLoEYeopUKsXhOyRnaaZSRoZdNtuAAX3Lt0Vw9laW4w1xq6MjnUncxSQbMhdhu5EDmb3UbGq3/cb55m033hpiUL86mts7fdebVNZGbHhczn/POYD5K9ezzYYDeeTthYx94sPyydXTcNznQf3qaGwpfZ+Jim6PMYmWtRzKFYNwYrKYZKmNiuaLo4dxw/d25/gvbsZ5X3VueJeL15rs2YCXX8r1Qibls/79Udt3fR7Qt46dNtmAfn1qGWpbTj0uP70fThyzOQdtt2Fi9YfF3lEMbejDW384vOQyU8Ci1S0ll5NTZoTtLq69ELKvYFleRZc6OtLiYqp4alJwzO6bcs3xuzHQw7zNZLwfUsUqiDKtD+2sp7ZAgrf9uH4B3XJRsu82I7jnlDGxlR9bUNFRbN+6CO5hCt5buLr0cmwEGfEWu1epFPzy8G1LFSm/XLJB6mQsiPbODJu7xO78nh8FoiCK4PfGpzMZ7xhEfbAXtVotiPraGj6+8ijX3w7bvntr13Tafb9ruzx965JrntnbP9yRehgV5YhBRNmxrYvATWUnSPMudmgK9ySIUsnKWJYsJsD5RnR0pvnZgVt7DkzjRhREEYI0Dq8XsrU9HajecgWpo67lxDGbF7wPlx+3S9fndKaAgrB9TlZBGEmeOPPLZaszio7Ifu+jerYpoLGlPaLSDO2d7s8/DHF14Nl7WQ4r282gPn6vzenXp5Zbfrhn0fPFxZQQfi2IAv1dF03twUZg5ZrwFOUo8+pv78oFR+1QoB4Y2K97JLTZ0OIjvvoKsCDclkCImi9sNJCT99uS4QNK37EwjnFFKpVizFbDoi/Yd/1FfidV9P0rhXKM1ZzX8MSZX+66534UlLiYEsLvfc/gPWJravOfwQSUbXGyKJvVd/fe3DNuMLBvHZceuzM/O3Brjtx5lLs8OS6m5GIQYRXnX0/YzWf53Z+vPX43Ljp6p2gsiJhGu7870l3xl4NyxcmcZF/BcgzWnFXstMnggr/5OT8qREEUIYgF4XVoc0AFUa7FK8uVxZSt5qR9RvP7o3bwSN/r/r5PyPXVnzl7/1Dn2akJ4V7YesMBHLdH/qqcbthL7a6rdOy31e3ZHrDtCNfzvGbUp1IwpCFYkkU5SaXy/fff23vzksst70zq3EpyY0k+zo8rkyuWUktEKXWDUuovScsB/l/aDBnPh1SxFkTA1n/pMTuFrCf4cWHjMFtEMPEuTIAybGZOykddW/vck6BY3nwYV0yKZJfwKDZOcIp27yljuOrbu5Zcb7bc8mQxuddtPvtxMUUrT1e58RQbDqXUcKXU3cA5ScuSJdCSFxFaEJU6MWnvrYbx75/vy0u/OTjQeX5HOPajwtyDsw75QiQugV023SBPnmKEXczRj7xbj/CnIHKL8i+Pl+JIpVKJLAK4/ahBbDm8gUNs2W+uOJapiCo9OlVOC8LxrOx1JuliqrSZ1K8CrwGPJi1IFt83voiLqdPxBp62/1bc/urHBY+vVAVRU8aAZZj1qHbedHDg0dTE8w5i8twVfHWnUaxpbiedybDR4H5AsNFjoIw3l/MO2m5D/v3WQtfj/Y787R15lE0oCQXx33MOIJ3JUFfEhHCmiEYlaraccryKNSlylFzQbLS4rJyyKgilVB3gtupXWmu9FjhMa73IsiIqAr8j32LzIJxpnecfofIURCrV3RGUc2nxIMT9stgbut8Vnq8/YXd++fC72RIC+2NHDOjLSdbquSMtxdAlT4ByAnWitkOz5/3hGzsWVhC+ZbBVEVVHSXk6SSc1NSlqfDyBlENDRCVr9rkEea6D+tbR2Bp8zohXFX46/7geT7ldTAcDn7v8ex9Aa72ozPIUJZCHqchifnb61tVyxsHbdP197mHbsqvl1gBIVZTzz064pug7BmH77PfFbLBNQqxJhegYPV9O/8W4hY38LNeRbWOD+xUOBDtX+ixcVsr1c3EKl59dLK9SuO67uZli+QOCaGRNOf4PdFLQujzur58+qEdYEFrrCcSn7GLBdxYTwS/M/uAH9atz/FaZtyl+C8L+2V9lztnDQe+d5+gtwFN1q/fsQ7/Ay7OXe5br5zoLTDwPLE/Gty3iLKv78x5bDGHpmhYWrYl2fSa/fGvPzairreGch6Z1fWe/Lj+Pf/iAelaub8v7fvwvD+Rr178C2LOYArgZfR/p/zw/1feKIHUl4vfGF0tzdcPeSThf5spUD+EVV5jT/G+WYvscoi7Pw0u0IAvJEjQIWWjmuZOwz6dY8fZyMxn/Lq+4cIvhZPFzDwrJb5+97ye7LE+uEhIVPGy44vX2pjTXiiKABVFK0c7OsFItiLjFspfvN1DvzBmP8t4FS3N1LcG93JzzopM3jhiE2+AnzpnLfnAOCuzy+LnsQi47t2eRxEQ5v78FOSYMoiCK4N+C8J4H4YZ9tFFTk4olEyNqwr8s/s4L6noxx+V+DhyC8Kin1CB1VO4B31lMtsKiVTw2C4LwrqqocA4KMgV+C0ruoC0bpA53fhC8ZPZzPb0lzRUArfXJScuQJVAMogTXhrOTqlQFEbtcHlZVIXI6C1KRupiCuAzc5C10tltH5IXfDjln0BHhs8oRMZNJ3oKwf47w3clRsDXu5XueH7byFAXdEH5KlD2pE8L3bQ8zQ7UqXUzxyhXG9eI8LKiM3kHqIOUECWYGkzHtczHgXBdTfh2FOvZizTffgkgWL5eXn8dQSP6c+9e1H0QAufwf6vs8Xy6mkPUWQxREEYJM1grad3qlJFamegg/Ko0zSF2qMvXqrEuPQRQv108b821B5Ljo/MvjXbdLkDppCyJV+Dr9KN9C8rutgZTUYoFdcshqrpWL3/tebC0m17Jtn6vFggidxeT3OEfKatCyo/aNl5rm6uca/CgWv2mu9rKCPCuveRaZTCan3KTjD+BwMZHKkamUZczc0qzLEYMotUwJUieE75nU6RAxCA93QIXqh/AmtM8Tw7iYSh3glZpBkiWs8vR1XoilNqJsQimHBZE03okFfiyIQhtW5cdwgj3X6F9cURAVTJRbI+Yd73Ax5fpRK1NDxB6DSLl/9kuSnVfWVRT0hfZzmf6D1O6fS8FZs7nHSWcxdX9OpXKfeynB+RqX5xJsHkT4ugvdUT8KT+ZBJIT/LKZMScHRCl2bL4/QMYgQaa6+J8rZzgm1nHVkFkSgWm11hPeZ58tQOK7lVU6Q25YJIE9c5M2DKPBbIQoHqfPvX1lmUntZRGJBVC6+YxCh1tn3fpkrkXJaNmHmQUTdb4WJQfg5I+jgwPdMatsbHeezStrLlBeML7ASakECBakDyJVQDEKC1Anh9yUrdS2mKtEPsWcxBXW9REHUWUxBO+ZSlobwKiuIFF76x/lbJpPxvXhgXOS2E0f8roRycxQs2SB1EAsihhiELxdTPIiCKEKsazHljCKrQ0OUcx5EuayVqOdB+LIgfNafJdRifTG+3YlbEI6RftCZ1H5cTNmHVK4spoKBc3ExVS5ht5H0dbzDxVQJ6YPFiL3PDmFV2Q8LM7L1qibMTOrgAwU/LpHgi/VFNZp13tOk4w9Q7JkVP99PZxzmeSY1zJOZ1Anh/76H6JiqMkgdTlC/Z4UKUpc6Uc4zZdI/3TGIYCNYX2sxBZShULmFBiHBgtTJL7XhNXGvpLWYXOKCwXYVTOZFFgsiIXzHIEJlzwTLZKkEyqnIwoyAIw9SB4pB+I9qBu3Q/AapvebWhMUtzbWSYhAQXJ6CKaUuFmy5Xs1CMvm5NAlSJ4TvGESIsu1FV4sFEdZtEXZl1qDnhMsm8yo7yOixeHle53nhP801WLnd5fu/cRWxFpPj7zhWQg6z5WhS4zwJUieE73kQIXqmsMsiJEncYqYKfI61zog7FDuF24VtaYiY5kHEdf8yFbBjUM7MbocwpVhObqcmPXiTIHUF4/e+h9swyJGJkfSwzAdx7yiXhNstqnqCLM0QdPVR/y6mfB96yVRgu8x3MQU7v9DxrjGIEtfjKgfiYkoI/xZE8LKrM821fHX5Xu4756/kerPuoGbxY6MKqubL4F1uFIOQSnQx+f0tS6FgfclZTCU8yoKz3JNcPia5qquDODtEryWLK5W4FVkSMYio6FqLKeB5fo4PtSd1bD6mmMoNgH2J9Cg9XjUuFmxZltoIeV7X+eJiSoYgM6kDl237XC0WRPxrMbl/9n9WcgSZSZ2JKQYRh37Iy2KiArKYbJ+d8pQS9Hdrf9UQgxAXU0IE2ZM6KFHto1tOKn411+hECUyw/Re6P/vq0ELIEFae/N+cE+WSn9LpPfs9miB19nOck2XtlHJXJYspIcrVcddUSZA6LP5vY/AgdaXo1iBLbeSmZXqfkUqFczEFWWojSOdUGau55rqYglJ4HoSbiymIXMFliQKZSZ0Qfu97qBcm6AqUvYxqsaqyhF1qoxgpCLBhkP08lyB1BPIYn3/luJjyfguqoQsQZvHFOFqsv4lyMVSMKIii+I9BlOpiCnx6VeH38kp2MSUZpA4Q1Awqp18LIsxaVsVw1pzOJL/UhvMexyFP92qu0ZftRinXENcAsy6WUkOilPoD8DNgMPAucJbW+sMkZQqymmtJ9fR0DeGTMEHq3IBlkmmu/o8NErNKpVKxxyCCUAHz5Fz2yLa5hnyc76edZF10QWIacXTUSRrSFWNBKKVOBn4EHAyMACYA45RSicrot3GE8oPm5MIHP7+aCDdRLvg5SY5sg/isg4iZIuRM6kAKK4BAFYDzvYxjYBDGguhpr3HFKAiMUrhcaz1Pa90B3ABsAWyWrFj+COVi6kUxCN8BZ9vncBPlkiPIjnJBSKX8ty/7LauNqE25Ko/EXUzdn52ruQZdLLFYHeUIUlfq619WBaGUqlNKDXH5N1hr/Ret9T22w48GVgILyyHbOYdt6y7zqIF53/38wK0BOOPgbdh+1CAAfnn4dgBsMawh59jj9tgUgPO+sh1nHrKNOf8gc/6eo4cCUF9Xw+hhDfz2yO0B2G3zIQCctM9oAI7aZRSD+xlv4FYjBuTJU19bw6Hbb+TrOu0cufMoAPay5CjED760BQCbbNAv5/tLj9mp6/Op+2+V89tlx+6c8/f39t7cl0xbDGugb51pll8cPbRLtq1t112TgutP2B2AL39hOMp6BgAHbDsCgC9tNcxXfX7Y3XoeWf56wm5do8pTvtx93dm2cKl17QdsO4LtRw12LfO0A8x5znu6z9b5cl/1rV0Z+/UdARjS0Cfnt6037L4vg/rWsd3IQYwc3BeAY622l70nAL86fDs2dtS5xxZDutqvG1d8y1xP9p7+8Zs78ufv7FLw+GN33yTvu4uP7m4rfetqOP2gbQqeD7n31Y3NhzbQr49pJ/tsPYyjrTr71tUwrKEegNP2L1zGVd/elZP32xLoblsbDTL3bfNh/QH4+YFGxu1GDsovwIUrv7ULFxy1Q973J+0zmvpa7652982H8Odv7wqQ9y6PHt7gdgp3/WRvAEYMrPclXxhS5ZzwopQ6HHje5acFWustbccdCIwDfq61ftBHuVsCH0+cOJHNNgtncGQyGRasbGJoQz3zVqxj5OB+ZIBNh/QveOzo4Q20daZZsa6t67iW9k5WrGtl1fo2NhvawNCGPl3HAl2fsyPqT1c10bdPDRsN6tf198jB/aivq6EznWHm4rWoUYNo70zT2NLB0IZ6Zi9tZIvhDbS0d5IiRb8+NTTU1zF/5Xrqa2uorUlRV5Ni0ZoWvrDRQNa1dNDemWbYgHpWN7fT0ZkmRStuCAoAAA2CSURBVIpNh/antiZFe2eaJWtaaGzp6Frrf+iAej77vJlhA/rwhY0GsXhNM0Mb6unXpzbnXixYuZ7OdIatRgzIsxIWrFzP8IF9Wd/awcjBuZ2SF8sbW2lu62SL4Q1dsm06pD8zFq9lxMC+NPStZXC/Pny6qomNN+hHXW0Na5ra6UinGT7QvOTtnWkWrGxidVMb6Qxss+EAOtIZGlvaGdSvD5+saiIF7LTJBvSvr/WUp6W9k3nL19NQX0s6k2HrDQeytqWdto40wwfUM33RWvrU1uQoKqdsyxpb2KB/Hz79vJmdNx1M37pa13uavd7OdIbNhvZn0eoWtrDazqermhi1QT9WrmujJgULVzez66YbsKyxlcaWDjYc1JdhA+pZ29LOsrUtfGEjI09HZ5r5K5voX1/LpkP6096Z5oPP1pBOZxjSUM/o4Q30qa3hs9XN1NWk6FtnPq9v7WT4wHq22XBgjmybW4OgT1Y2UVdr2tqa5nY+b2qnob6WnTYZzKr1baxc30b/PuaejR4+gKa2DmYubmSXTTegT22KBSubGNLQh5Xr29h4g36sbe6gqa2DxpYOdtl0g6JxuRXrWmlq7ey6P8saW2ior2NgXzOYymQyfLxiPbVWOZsO6c8nq5qoq6lhi+ENpNMZZljv14p1rQzpX0//+lpa2jtZtb6NTWzv/icrm1jX2kFLRydq5CBWrW9j2IB6Zi1pRI0axMp1rYwePqDrObV1phk9rIGFnzczengDy9e10tqeBqCtM01LeydNbZ1sMayB+toahg6o76pnkyGm3djJtrf1rR0MG1BPc1snGw3ux5I1pl0Va8OFWLhwIYcddhjAVlrr+c7fy6og/KCUOgm4GThba323z3O2pEQFIQiC0NsopiAqLYtpLPBL4Bit9QtJyyMIgtCbqRgFoZT6CfArYD+t9ayk5REEQejtVIyCAC4ABgFvKaXs3++ttZ6ZjEiCIAi9l4pREFrrwmkUgiAIQtmppHkQgiAIQgUhCkIQBEFwRRSEIAiC4IooCEEQBMEVURCCIAiCK6IgBEEQBFdEQQiCIAiuiIIQBEEQXBEFIQiCILgiCkIQBEFwRRSEIAiC4IooCEEQBMEVURCCIAiCK6IgBEEQBFdEQQiCIAiuiIIQBEEQXBEFIQiCILgiCkIQBEFwRRSEIAiC4IooCEEQBMEVURCCIAiCK3VJC5BFKdUX+CtwPFAPvAT8Qmv9WZJyCYIg9FYqyYIYC+wIKGBDYCXwt0QlEgRB6MVUkoL4E3Ck1noVMBIYDKxIViRBEITeS1ldTEqpOmCgy09prfVaoFkpdRHwR2ARcGAZxRMEQRBslNuCOBj43OXf+7ZjrgIGAI8C45VSfcosoyAIgkCZLQit9QQgVeSYFgCl1PnAGcAuwDvxSycIgiDYqZgYhFLqTqXUGbav6jDyrU5IJEEQhF5NxaS5AlOB85VS/wOWATcAk7TW85IVSxAEoXdSSQriVmAj4DXMPIjnMHMiBEEQhASoGAWhtc4Al1j/BEEQhISpmBiEIAiCUFmIghAEQRBcEQUhCIIguCIKQhAEQXBFFIQgCILgiigIQRAEwRVREIIgCIIroiAEQRAEVypmolyJ1AIsWbIkaTkEQRCqBlufWev2e09REBsD/OAHP0haDkEQhGpkY+Aj55c9RUG8CRwALAY6E5ZFEAShWqjFKIc33X5MZTKZ8oojCIIgVAUSpBYEQRBcEQUhCIIguCIKQhAEQXBFFIQgCILgiigIQRAEwRVREIIgCIIroiAEQRAEV3rKRLlQKKX2AG4FdgLmAKdrrd9IVqrSUUrtD1wLbA+sAK7WWt+qlBoK3AkcCqwBLtZa32Gd0xe4GTgOaAdu1FpfnoT8paCUGgl8AJyitX5GKbUlcAcwBjOR8tda62esYwvej2pAKbUZ8A/gQGAt5jnf2JOfs1JqP+BGYDvM87xYa/1gT7xmpdQY4Amt9SbW36GuUSmVAq4ATsP0+fdi3oOik4p7rQWhlOoHPA3cBQzBNLrHlFL1iQpWIlYjegpzPUOB44ErlVKHA7cB64CRwHeAq5VSu1qnXg6MBrYC9gdOU0p9s8ziR8EdwHDb348AU4FhwLnAg0qpDa3fvO5HRWO99E8AMzHX+zXgIqsD7ZHPWSlVi7nmq7TWgzEd3j3WIKDHXLNSKqWUOgV4DrD3R2Gv8Uzg68CuwA7Al4Ff+JGl1yoI4BAgrbW+RWvdrrW+E1gJVGzD8cloYJzW+gGtdVpr/Q7wIrAfcCzwJ611i9Z6KvAg8FPrvB8CV2it12it5wA3AT9LQP7QKKVOB9YDn1p/7wDsAlxiPeP/AS8DP1JKDcT7flQ6XwI2AX5nXdt0YF/gM3rucx4CbAjUWQoyDbRhltfpSdf8e8xgpsvK8dFeva7xJOB6rfVirfUS4Ep8Xn9vVhDbAzMc32mMu6lq0Vq/q7U+Kfu3ZVEcAKSAdq31PPvhwE7WMSPJvR9VdS+UUtsC5wFn2L7eHpivtW62fZe9rm0pcD/iljUi9gSmY0aRS5RSs4F9MJZSj3zOWuuVGDfKQxg3yiTgLGAEPeua7wR2J3d9pILt1cc1Ovs6DexoKVlPerOCGAA0Ob5rAhoSkCUWlFIbYNxob2OsiGbHIdnrHWD72/lbxaOUqgPuB87VWq+y/eT1jAdQ+H5UA8MwVvAKYAvgZOBvwEB67nOuwch7PEbmbwLXA4PpQddsjfSdi+R5tddi1+h8D5owfX/fYrL05iB1E9Df8V0DxsdX9SiltgKewSzhewLG91joerONpz8m2Gn/rRoYC7yrtf6v43uvZ1ztz78VWKW1vtL6e7JS6lHgYnruc/4W8CWt9fnW3+OUUuOAi+i515ylWFuGwtfoPLcB6NBatxSrtDdbEDMB5fhOke92qjqUUnsCU4DxwLGWi2UO0EcptYX9UGCGNepeRu79qKZ7cQLwPaXUaqXUasyI+l+Ya9jSyvDIkr2ugvejTDKXigYGWNZTllpgGj33OW9B/qi3HXiHnnvNWUp5f519nbK+K0pvtiBeAPoqpc7GpAqehPHjjU9UqhKx0jyfBa7VWv85+73WulEp9SQmo+mnGP/k94GjrEPux2TBfAeTFXMW8H9lFT4kWuvt7X8rpeYDZ1lprt8BLlVKjcWkBx4M/MLH/ah0ngc+B65SSv0Ok8Z7HPAVYEt64HPGXPOVSqmfAHdj0nuPwzzXLemZ1wyU/P7eD5yvlHoBo1AvAO7zU2+vtSC01q3AkcCJwCrgbOBorfX6RAUrnVMxmR5jlVLrbP8ux2Q89AEWAo8C52utp1jn/QGYDcwCXgVu01o/Un7xI+dbwG6YEdb1wIla60+t37zuR0VjWYUHAztjru1B4BxrHk+PfM5a6w8w6Z3nYuYB/B34sdb6LXroNTsIe403A09i0r1nAK8B1/mpUDYMEgRBEFzptRaEIAiC4I0oCEEQBMEVURCCIAiCK6IgBEEQBFdEQQiCIAiuiIIQBEEQXOnNE+UEIQdrVvL/YdY1Gg2sxkycHKu1XpCAPC8Bb2mtf1PuugUBxIIQBDtXYCYa/gqzHMExwMbAy0qpil3cTRDiQiwIQejmVOCXWutx1t/zlVLHY2YqHwX8JzHJBCEBREEIQjdp4ECl1IPZ7Ri11quVUjsDiy0X1OWY5Vk2BpYDd2utfw+glLobs2zLEMwCgssw+1NsAlyCWXb5Lq31r23HtwEbAEcDC4ALtNaPuwlnrUH0e2BTzF4QF2itJ0R8DwShC3ExCUI312G2sVyglPqnUupEpdQwrfVsrXUj8Dvgu5hF0rYFLgUuUEodaCvjTMxKmbtiVhn9F0ahfA04H/iVUurLtuNPxiiSPTB7BT+ilNrdKZhS6ijMPuMXWGXfBzyjlNolqosXBCeiIATBwtpb4buYRc9+glkAb7FS6kpr960PgZO11q9qredrrW8BlpC7O5nWWl+jtf4Isz/2BsB5WuvpWuvbMcpgR9vxH2EW2Zultb4CmIz7tqe/A67RWv9Haz1Xa30j8Bhm4TpBiAVxMQmCDWsFzEeUUoMwO7adgumcP9Na36SUOkQpdQ0miL07MAqzD0MW+5aQ2Y1c5tu+ayZ3T4PXHLuHvQns5SLajsAYpdSFtu/qMft+CEIsiIIQBEAptSvwU6312WDW3weeAp5SSj0DHKGUGoFZFv5O4GHgl8DLjqI6XIpPe1TtPL4W6HQ5rg7jXnra8X2rR9mCUBKiIATBUAucpZR6WGv9quO3tZiO+EzgN1rruwCUUkMwm0wV3fzdgz0cf4/B3SqYCYzWWs/NfqGUugLj4rqxhPoFoSCiIAQB0FpPU0o9DvzHcuO8BAzGbCp1NLA/sC9wlFLqFcyuXVdgNnApuvm7B2Os3e4exgSz98QErp1cDTyklJoBTAS+CfwW+HoJdQuCJxKkFoRuTsTsvvVr4APgFeAg4DCt9buYjvsLmGD1v4G3MTt7Oa2AIDxnnf8eZmLeEVrr2c6DrNTXczCZUDOAnwMnaa2fLaFuQfBEdpQThISw5kEM1Fp/J2lZBMENsSAEQRAEV0RBCIIgCK6Ii0kQBEFwRSwIQRAEwRVREIIgCIIroiAEQRAEV0RBCIIgCK6IghAEQRBc+X+bq355jOJ3JQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "bad_chain = norm().rvs(1000)\n", "bad_chain[bad_chain<-1] = -2\n", "\n", "plt.plot(np.arange(1000),bad_chain)\n", "plt.ylim(-3.5,3.5)\n", "plt.xlabel('Sample')\n", "plt.ylabel('$\\\\theta$')\n", "plt.title('A Chain that never explores below -2');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Without getting too much in the weeds at this point, it may be possible to transform the parameter in the above case to sidestep the problem we can see above.\n", "\n", "### Without ergodicity you are screwed\n", "\n", "If we can't assume ergodicity, then we can't use our chain for inference. You might be able to resort to practices like *thinning* to whinnow your chain down to something more like iid (and then rely on SLLN), but it will be difficult to convince your audience and it will require running very long chains. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## MH MCMC Inference\n", "\n", "Having established these three properties and having used random-walk MH, we are assured that our chain represents **somewhat** correlated random draws from the posterior that can be used for all manner of inference:\n", "\n", "* Means\n", "* Credible Intervals\n", "* Medians\n", "* Standard deviations\n", "\n", "> Provided that our chain has converged to the limiting distribution, which is a numerical issue." ] } ], "metadata": { "hide_input": false, "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.8.12" }, "latex_envs": { "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 0 } }, "nbformat": 4, "nbformat_minor": 4 }