{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "
\n", "\n", "# Specifying Sample Variables\n", "Author(s): Paul Miles | Date Created: July 12, 2018\n", "\n", "Many models require a set of parameters as input; however, it is not always of interest to estimate those parameter's posterior distribution. The [pymcmcstat](https://github.com/prmiles/pymcmcstat/wiki) package will allow you to send static parameters to the model along with the parameters being sampled. We demonstrated this feature by considering a simple linear model." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Import required packages\n", "import numpy as np\n", "from pymcmcstat.MCMC import MCMC\n", "import matplotlib.pyplot as plt\n", "np.seterr(over = 'ignore');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define a linear model function and corresponding sum-of-squares function." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# define test model function\n", "def linear_model(xdata, theta):\n", " m = theta[0]\n", " b = theta[1]\n", " nrow, ncol = xdata.shape\n", " y = np.zeros([nrow, 1])\n", " y[:, 0] = m*xdata.reshape(nrow,) + b\n", " return y\n", "\n", "def ssfun(theta, data):\n", " xdata = data.xdata[0]\n", " ydata = data.ydata[0]\n", " # eval model\n", " ymodel = linear_model(xdata, theta)\n", " # calc sos\n", " ss = sum((ymodel[:, 0] - ydata[:, 0])**2)\n", " return ss" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Generate Synthetic Data\n", "We generate sythetic data using the linear model and adding observation errors\n", "$$y_i = 2x_i + 3 + \\epsilon_i, \\;\\; \\epsilon_i \\sim N(0, 0.1)$$\n", "We plot the data and linear model and see randomly distributed noise with respect to the model response." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xmc1XP///HHa6aF4gqVLRLqQhQyaGwXV7ZcyNL1la2LEBeupFJN2RKVUkI/dZUwJfuaJSKy5IgpJWSJrkYJFaK9mXn//vic4XQ6Z87nzJx9nvfbrducOfOZz3l/btXzvM/r817MOYeIiOSWvHQ3QEREEk/hLiKSgxTuIiI5SOEuIpKDFO4iIjlI4S4ikoMU7iIiOUjhLiKSgxTuIiI5qE66XrhJkyauRYsW6Xp5EZGsNGfOnJXOuaaxjktbuLdo0YKSkpJ0vbyISFYysyV+jlNZRkQkByncRURykMJdRCQHKdxFRHKQwl1EJAcp3EVEcpDCXUQkhQIBGDrU+5pMaRvnLiJS2wQC0KEDbNoE9erBjBlQWJic11LPXUQkRWbO9IK9vNz7OnNm8l5L4S4ikiLHH+/12PPzva/HH5+811JZRkQkRQoLvVLMzJlesCerJAMKdxGRlCosTG6oV1JZRkQkByncRURykMJdRCQHKdxFRJIsVROXQumGqohIDIFA9Ue4pHLiUihfPXczu87MPjWzz8ysZ4Sfm5nda2aLzOwTM2uX+KaKiKReZTjfdJP3tared6QeeionLoWK2XM3s4OAK4AjgE3Aq2b2knNuUchhHYFWwT9HAmODX0VEslqkcI7U847WQ6+cuFT5fDInLoXy03M/AJjtnFvnnCsD3gbOCTumEzDJeT4AdjCz3RLcVhGRlPM7qzRaD71y4tLgwSElGeeS3m4/4f4pcKyZNTazBsBpwJ5hxzQDvgv5fmnwORGRrBYxnCOo6k2gsBCKiqCw0edw2mnwyCNJb3fMsoxzbqGZ3QlMB9YC84Dy6ryYmXUHugM0b968OqcQEUk5P7NKq1xaYOVKuOUW+O9/Ybvt4Nxzk9haj6/RMs65icBEADMbgtczD7WMLXvzewSfCz/PeGA8QEFBQfI/l4iIJFH4KJqt3gQ2boQxY7xu/5o1cNVVXsg3bZr0tvkKdzPb2Tn3k5k1x6u3tw87ZCpwrZk9jncjdbVzbnlimyoikhg1GdoYeo6oQxydg+efhxtugG++gVNPhZEjoXXrBF1BbH7HuT9jZo2BzcA1zrlfzewqAOfcOOAVvFr8ImAdcGkyGisiUlOJGncedRTN3LnQqxe8/TYrdj6QFaNepfX1pyT4KmLzW5Y5NsJz40IeO+CaBLZLRCQp/A5tjCV8iONJB34Plw6E4mI2N2pM77pj+e/Ky8kfWIcZ7VMzcSmUZqiKSK3id9x5rNJN5Q3U96av47xlI2l+wZ2weTP06cN92w7k/jsaUV4B+TV4A6kJhbuI1Cp+NszwVbqpqKDwm0cpfKAIli71RsAMHw777ENhAOqNSP3EpVAKdxGpdWINbYxZunnvPa+u/tFH0K4dTJkCxx23xflTteNSNAp3EZEwUUs3ixdDv37w1FOw++5QXAwXXQR5W88HTdWOS9Eo3EVEwmzV8z7wN+g/BO6+G+rUgVtvhT59oGHDNLc0OoW7iNQK8Y5tLyyEwsPLYOJE6HQTrFgBXbsy59whTP+sGcd/kt6eeSwKdxHJedUa2/76615d/dNP4Zhj4OWXCZQdnpa12atDOzGJSFaKZ3cjv2uqBwIw/vqF/HLUP+Dkk2HtWq++/s47cPjhaVubvTrUcxeRrBNvT9zP2PaPpq1k3hmDuKJ8LGtpyJJrhrPXyB5Qv35c58kU6rmLSNaJtwdd5bK9mzbBqFEcdE4rupffzwSuYP+8r3m02Q1bBHvM82QY9dxFJOuE96AbN/ZKNFXdLN1qaKJz8MIL3uJeixax8chTOGbeSOaXHVhlrzzdQxz9UriLSNYJHarYuDH07BnnTc558+D6670THHAAvPIKO3TsyJgErBaZKVSWEZGsVLm70apVcZRoli+Hyy7zZpUuWAD/7//BJ59Ax45bnDPbgx3UcxeRDFCT9dV93eRcv95bT33YMO/AXr3gxhthhx1q3PZMpXAXkbQKH/kyerTXG49nslHUdVycg8ceg/794bvv4OyzvcW9WrZMzsVkEIW7iKRV6MiXjRvh2muhoiK+SUKhNzkrPwWc3jhAmwevh9mzvTLM5Mnwt78l81IyimruIpJWlWWV/Hxv/a3y8upPEgoE4JITlrDPwC60ufIo1n5RykvnPkTg3o9qVbCDeu4ikmZVjXyJa5LQb7/hioYxf+MoKsjjNm5m1Nq+rHm+IfVeif4pIBH7qWYihbuIpF1oWaVNmzjDtrwcHnoIbryRo378kUfzL2KAG8KyvD2pqPBKPNG200vUfqqZSOEuIhklrklCM2Z4I18++QSOPhpefJG9yw7nypn+PgUkaj/VTKRwF5GUSVgJ5MsvvZmlL74ILVrAk09C585gRiH+PwVk01ox8VK4i0hKJKQE8vPPMGgQ3H8/bLst3Hkn9OgB22wT8fBYnwIyYTu8ZFG4i0hK1KgEsmmTF+i33QarV8MVV3iPd965xu3KlrVi4qWhkCKSEqFDHn2XQJyDqVPhoIO8tWAKCrx1YcaNS0iw5zL13EUkJeIugcyf790sffNN2G8/eOklOO00MEtBa7Ofwl1EUsZXCeSHH7x1Xx58EHbcEe69F666CurWTUkbc4XKMiKSGdavhyFDoFUrKC72xjEuWkSg4D8Mvauur+305E/quYtIejkHTzwB/fpBaSl06gQjRkCrVjk9ySjZ1HMXkfT54AM46ig4/3zYaSevvv78817vnfi305M/KdxFJPVKS+GCC7xu+P/+BxMnQkkJnHDCFodVa4SNACrLiIgPCZtZ+vvv3sSjkSO972+80SvHbLddxMNzeZJRsincRaRKCal7l5fDww97Yf7DD3Dhhd7N0+bNY/5qrk4ySjaVZUSkSjWue7/5Jhx2GFx+Oey9t1dnf+QRX8Eu1adwF5Eqhde9GzeGoUOpcmhiIADjen/Nz8d28rr9v/4Kjz8Os2bBkUemrO21ma+yjJldD1wOOGABcKlzbkPIzy8BRgDLgk+Ncc49kNimikg6VLWZRqQSzYev/ULJP27jqvIxbGAbllw1hL1G9fQW+pKUidlzN7NmQA+gwDl3EJAPdIlw6BPOuUOCfxTsIjmksBCKiryNq6OWaDZvhvvuo83ZLbm6/F4e5hL2z/uaYRQxdPS2moSUYn5vqNYBtjWzzUAD4PvkNUlEMlXE9c+dg5dfhj594Msv2VjQgb8tGMXcsrbk53ubJJWVJWYSUq5uiZcMMcPdObfMzO4CSoH1wHTn3PQIh55rZscBXwHXO+e+S2xTRSTdwocmNlj0CYu79mbvRW94i3tNncoOp5/OPR8YM2d6w9knTEjMTkearRofP2WZHYFOwN7A7kBDM7so7LAXgRbOubbA60BxlHN1N7MSMytZsWJFzVouImlRWAhF3X5kn2HdOajroTRaNIfede/hgwkL4IwzvN2QgmWcrl0TNwlJs1Xj42e0zInAYufcCufcZuBZ4KjQA5xzq5xzG4PfPgAcFulEzrnxzrkC51xB06ZNa9JuEUmHDRtg2DBo1YomLz3EGOtBSxZxT0UPih+tu9Uomsqe/uDBNe9pa7ZqfPzU3EuB9mbWAK8s0wEoCT3AzHZzzi0PfnsmsDChrRSR9HLO26e0Xz9YsgTOPJNPLhhB0aV/ZdMmqqytJ2oSkmarxsdPzX22mT0NzAXKgI+B8WZ2G1DinJsK9DCzM4M//xm4JHlNFqmd0nYz8cMPvV2Q3n8f2raFN96ADh04FJjRnITX1qui2ar+mXMuLS9cUFDgSkpKYh8oIum5mVha6hXOH30UdtkFbr8dLr3U66ZnQvtqKTOb45wriHWc1pYRyQI12lw6XmvWeIt73XWXV44ZMAD694ftt4/6KyqZZB6Fu0gWiDi+PNHKy2HSJBg4EJYvZ+WJXXji0GG0O30vCqPn+h9UMsksCneRLJD0nvHMmV5dfd48aN+eBbc+w5E9C9n0FtQbozJLNtLCYSJZonLseEJD9uuv4eyzvU0yVq3y6uvvv89LqwqTPqY8EIi9AJlUn3ruIrXRL794g8/HjPHqPLffDr16/bG4V7LLQLoBm3zquYvUJps3e4HesiWMHg1du1Ly+CKG5g0kMO/PVRsTOfkoEs02TT713EWyUNxj3p2DadOgd2/44guvDDNqFIH1h0TtQSfzBmlKbhDXcgp3kQwWKcTjLml8+qlXcnn9dWjVCl544Y81YGYO3bIHPWlSaoYzauhk8incRTJUtBD3O+b9o5d/Iv+2mzm0ZALWqBHcfTdcfbV3sqDQHnSil+eNRUMnk0s1d5EMFa0uHXMBrQ0bWHLNcP56eivafDiR++1aim9axND1PQnMqbfFoaG19W7dvGBXHTw3qOcukqGi1aWjljScg2eegb592WvxYl7kDG5gOF+7/cnvBxUVkXvklT3oQACKi1UHzxUKd5EMVVVdOrSkEQjAl498xDmzevGX+e9BmzZ8fs/rnNf/RDZtgjzzeuMVFVWXcVQHzy1aOEwki815YSlfnDuAC8sn8yM7s6bf7ex7RzfIz//jZmz4ptajR3vzlRTg2UkLh4lkoIQt27t2LQwfTtuhIziwvIKh9Gd4XhF9G/2FouCijaG9+zZtIge9Jg/lLoW7SIokZFZmRQVMnuyt1Pj996zucB7HvDeMRWUtqqyTVwb90KEpXF1S0kqjZURSpMazMt95B444Ai65BPbYA2bNoskbj/PQWy3+mEkKVa/Xoq3qag/13EVSpNqzMr/5Bvr2hWefhT33hClToEsXyPP6ZqGjXWJ9MtBN09pD4S6SInEH66+/wh13wL33Qt263mD0Xr2gQYOIh/ud3KTJQ7WDwl0khXwFa1kZjB8Pt9ziDWu55BJv1cbdd6/y17Rei4RSuItkkmnTWHdNHxos/pzV7Y6n0WsjoV07X7+qkouEUriLZILPPvNWbHztNb63lvTNe45XP+/EjI1GPBmtkotU0mgZkTgkfPegFSu8xbwOPhhmz+aN00bRxj7juYqz2LTZtL6LVJt67iI+JXT3oI0b4b77vJuka9fCv/8Nt9xCw6+bYG9BvurmUkMKdxGf/I5GqZJz3pDGvn3h22/htNPgrrvggAMAKGyiurkkhsJdxKfqjkapXHLgH7vOoe3DvbzJSAcdBK+9BiefvNXxqptLIijcRXyqzmiUQAAu/vsybt44gLZuEpt3aErdcePgssugjv77SfLoX5dIHMKX2o0W9IEAzJq+lv1fuov5G4ZThzKGWz/q9Cii15WNtjpvwhYUEwlSuItUQ1U3VwOzKnjghCkM2lzEHizj6bx/0p87+b7+3sw4Nb5ziVSXwl2kGqLeXH3vPVp0uZ6Jm0v4iAIuzHuc1t2P4bLmW/fKK3vrpaVaqVEST+EuUg3hN1dbVHzLwjb9OODTp9mpaTMuqzuJSeUXUrd+HsO6Ri7bVPbW8/P/LL9r+KMkisJdpBoqb64GXl3N3z8YwgE3jqaMOgyuM4iTn+jD5ds0oOXM6DX00J4/wBVXQPMIvXuR6lK4i1RHWRmF8x+gcOzNuJUrKbZ/MdDdzo+uGXU+gKKiqkM6vOffNULvXqQmFO4i8Zo+3Vt697PP4LjjWHDp3Vx9dbu4xr9rkS9JNoW7iF8LF0KfPvDKK7DPPvDMM3D22bQ1Y8Z+8Qe1JitJMvkKdzO7HrgccMAC4FLn3IaQn9cHJgGHAauA85xz/0t4a0USzNf48pUr4dZbYdw4aNgQRoyA//wH6tf/4xAFtWSamOFuZs2AHkBr59x6M3sS6AI8HHLYZcAvzrmWZtYFuBM4LwntFUmYmOPLN22CMWPgtttgzRq48kov5Js2TVeTRXzzW5apA2xrZpuBBsD3YT/vBNwafPw0MMbMzDnnEtJKkSQIH6s+aVKwF/83R+GPz8MNN8A33/BLYUeePmIEB513IIXKdckSMcPdObfMzO4CSoH1wHTn3PSww5oB3wWPLzOz1UBjYGWC2yuSMKEjVvLz4aGHoM3muRxNL6h4G1q35vNRr1Iw8BQ2fQj1xmv2qGQPP2WZHfF65nsDvwJPmdlFzrlH4n0xM+sOdAdo3rx5vL8utUAq11gJHbGyeuH37P/IQLq6YlbRmEePHUvpyZez5Ks6mj0qWclPWeZEYLFzbgWAmT0LHAWEhvsyYE9gqZnVARrh3VjdgnNuPDAeoKCgQCUb2UJ111ipyRtC4cHrKHxjJOVPDaPMlTHK+nBnnYH8/mEjyt7X7FHJXn7CvRRob2YN8MoyHYCSsGOmAv8CAkBn4E3V2yVe1dkMo9qLblVUwKOPerONli4l/9xzmX/enWxetC+dS2HCBM0elezmp+Y+28yeBuYCZcDHwHgzuw0occ5NBSYCk81sEfAz3mgakbhUZzOMau2ONGuWNwnpww+hXTuYMgWOO452QDu8N4ziYs0elexm6epgFxQUuJKS8A8AUtvFW2IJ77mPHg2rVv35xrDFuRYvhn794KmnYPfdvZ2uL7oI8rbeJ17rq0umMrM5zrmCmMcp3CUdEhmeledq3Bh69vxz9IsZlJVB47q/MbfzEJo9NdoL8r59vWGODRsm4EpEUstvuGv5AUk5v3Xy0DcAiP5mUDk7dOjQP0s0FRWQ58rpxkQGl9/ELo/85NVXhgwhUNqMmfeqVy65TeEuCRerV+6nTh6+3nllL7yqN4PQmv1J9gbDy3vRxi1gVt4x/DjhZdp2K9CuR1JrbF1sFKmByvC86SbvayCw9TGVIZyfH/3GaegbwObNW78ZRFJYCLMmfsEXLU9nWtlJtNptDc+e/xR5775D224FW523qnOJZDv13CWh/PTK/Sx3Gz57NLTnHnEUzapVMGgQh44dCw0awPDhbNOjB+eELO4Vfl6NW5dcpnCXhPIbnrFWUQx/A4AobwabNsH998OgQfDbb9C9u/d45519nVclGclVGi0jCZeKYYSB9x3Lxk7lH2/fwLbffQ2nnAIjR8KBBybnBUUyhEbLSNoke23z+cXz2NitF50r3mKhHQAjX+GAXh2r/B2NW5faRuEu2WP5crjxRto+9BCr3E5cwxgmWndu2ViXA6r4NY2QkdpIo2Uk861fD3fcAa1aweTJLO/Si7bbfM1/868hr37dmDdFNUJGaiP13CWiZJQxop0z6ms5B489Bv37w3ffwdlnw/Dh7N6yJc/E0T6NkJHaSOEuW0lGGSPaOaO+1vvve4t7zZ7tLe41eTL87W9/nC+eur5GyEhtpLKMbCUZZYxo5wx//uPn/gddusDRR0Npqbc90kcfbRHs1VFY6K3uq2CX2kI9d9lKMsoY0c5Z+Xz9jb8xwIZy1T13Q34e3Hyzt7jXdtvV/MVFaiGFu2wlkWWM0Hp6pHMWHlHOgh4PsuuYG2m49idvCd4hQ2DPPWt6GSK1msJdIkrEWPVI9fSiopAD3ngDevVi3wULvDLM3S/B4YfX7EVFBFDNXZIoau3+yy/hzDPhpJPg99/hySfh3XcV7CIJpHCXpAlf/fHEdj/DddfBQQd5ST90KCxcCP/8p7cymIgkjMoykjSVtft3Zmzm/1bez97nD4LVq70dpwcNgl12SXcTRXKWwl2SxzkKV75E4eQ+8NVXcOKJMGoUtGmT7paJ5DyVZSQ55s/3aupnnumVXF56CaZPV7CLpIh67pIQlUMeT2rzAwUv3AQTJ8KOO8J998GVV0LduuluokitonCXaqsM9MaNod91G7h6493s54ZQkb+BvJ49vb32dtwx3c0UqZUU7lItf4xh3+g4jyf4uKI/LVjCC3Ri+XUjuGpkq3Q3UaRWU81dqmXmTDhk42zeqTiaKRXn8ys7cmLem5y/7fMc3FnBLpJuCneJX2kpV7x1Ae9XtGdvFnNV3Yl8NLaEDrefoI0wRDKEyjLi35o1MGwYjBxJE2DpvwbyRIt+/OuU7RXoIhlG4S6xlZdDcTEMHAg//AAXXABDh7JH8+b0TnfbRCQihbtU7a23vE0z5s3z6i3PPQft26e7VSISg2ruEtnXX8NZZ8Hf/w6//AKPPw6zZinYRbKEwl229MsvcP310Lo1vPmmt7jXF1/AeedpcS+RLKKyTI6Ke4PrzZtZ3G8cu467lW02/opddhkMHqzFvUSylMI9B8W1wbVz8PLLrL+mD3uXfskMOlBUfxT3XNqWQuW6SNZSWSYH+d7gesECOPlkOOMM1q2Ds/KmciKvM7esbdTfCQS8Sk0gkJy2i0hixAx3M9vPzOaF/PnNzHqGHXO8ma0OOebm5DVZYgnfJGOrDa5//NFbzOuQQ2DOHLjnHr5+dgHT659Bfr5F3RS78hPBTTd5XxXwIpkrZlnGOfclcAiAmeUDy4DnIhz6rnPu9MQ2T6oj6gbXGzaw5PrR7PrQEOqWrSevRw8vqXfaifbE3hQ70icCTV4SyUzx1tw7AN8455YkozGSOKEbXAfed/x435N0eL0fe61awlTO5Mb6I/jv//2Vwp0i/04klZ8IKmv5kXr3IpIZ4q25dwEei/KzQjObb2bTzOzAGrZLEmTBxA/h2GM46/EuLF7ViJPsDTrxAp+X/TV6LT6Kyk8EgwfHuEkrImnnu+duZvWAM4GiCD+eC+zlnFtjZqcBzwNbLQ1oZt2B7gDNmzevVoPFp+++g6Ii2kyZwg/swuVMoNguxerkk19R/Z53rN69iGSGeMoyHYG5zrkfw3/gnPst5PErZna/mTVxzq0MO248MB6goKDAVbPNUpU1a2D4cBgxApxjadci2j1ZxM+bt6dePRg9GlatimP8u4hkpXjC/XyilGTMbFfgR+ecM7Mj8Mo9qxLQvloh2oSjuCYiVVSw6KZidrlvINv/vhy6dPEW92rRgheuinNCk4hkPV/hbmYNgZOAK0OeuwrAOTcO6Az828zKgPVAF+eceuY+RJtwFNdEpJkzWXNlL1p+9TGzOZJ+9Z5haI9CClt4P45WSol7FquIZA1f4e6cWws0DntuXMjjMcCYxDatdog2vNDXsMNFi6BvX3juOcob7clFNoUp7nzyyy3mMMV43jz0JiCSfTRDNYUize6MNuEo/PnGjUN+99dfoXdvb3Gv6dPh9tv54vkveXabC6qchBTK7yxWTVwSyU5aWyZFovWUwyccgRfixx//5/ONG0PPnlCxcTMr8v/LofVupf66n/np9G7sMv522HVXjiT2JKRQfsesa+KSSHZSuKdIVSFZGfKR3gCKimDoEEeHjdMYXtGbAyq+4K2yE+hjo1j4xiHMWAyFu255Hj+izmINo4lLItlJ4Z4ifkIy4hvA9p/y7+d7UVTxOl/RirPzXmCqO4OKCiN/E0yaVP16uJ83A79vAiKSWSxdg1oKCgpcSUlJWl47XWLdmAztuTer+xOzT72ZXadOgL/8hf91vZknmlzDjrvUo2dP75j8fG//jLIyHyNqRCQnmNkc51xBrOPUc0+hWD3lwkJ4c9pG1g27h+PevYM6L66Fa66BW26hRePG9Ase16aN9yZRWgoTJqgeLiJbU7hnCufgmWdo37cvLF4Mp5/uzTLdf/+tDg2t0RcXqx4uIltTuGeCjz6CXr3gvfe8bvn06XDSSTF/TfVwEYlG4Z5EVdXYAwGY88JSzps/gKavToadd4bx46FbN6+Y7pMW8hKRSBTuSVLVDNDZb67lrVOG07NsBHlUsLRrEXvc1x/+8pf0NlpEcoZmqCZJxBmgFRXw8MO0PqsVA8puYypncmDeF0zef4iCXUQSSj33JAkf135Go3fg8Oth7lxofQQnLHqad8uP0o1QEUkKhXuChdbZZ8yAec98w//N6Uvja56FPfeERx5h+/PPZ8jsPN0IFZGkUbgnUGidvUnd1cw953YKn74X6tb19qbr1QsaNAB0I1REkks19wSaORPKN5ZxRflYFmxoyW6PjYQLL4SvvoIbb/wj2EVEkk3hHkGkpXn9OGubV/nYHcxYrubzvIP45ME58OCDsPvuyWmoiEgUKsuEiWsHpEqffw69e3PAq6+yodm+PH3CczT7dycOPspS0mYRkXAK9zB+1y8PBGD2SyvosvAWdp06HrbbDkaOZJtrr6VzvXqpbraIyBYU7mH8LM37wdsbmXriffQru53tWMPyzv9mt7G3QJMmqW6uiEhEqrmHqVyvZfDgCCUZ5+DZZ9nvnNYMLbuBWRzNIXkLeLjdfQp2Ecko6rlHEHGY4pw53lDGd96h/t4Hcsaa15hWfrImIYlIRlK4x7JsGQwc6G151KQJjB1Lg8svZ8BHdThqpiYhiUhmUrhHs26dt5768OHeVkc33AADBkCjRoAmIYlIZlO4h6uogEcf9XamXroU/vlPGDYM9tkn3S0TEfEtZ2+oxjsRKRCAyd3fZc2BR8LFF8Ouu8K778KTTyrYRSTr5GTPPd6JSHOf/pbl5/Xj4oqnWUYzlt80iVa3Xgh5eTE3tRYRyUQ5Fe6VQVxa6m8i0oevr6Zi8B0cNuse9quow80M4u68PgzYtgFFedWcrSoikgFyJtxDgzg/H+oEryziUMWyMr4d8AAtRtzMzqxgUt6/uLXeHZSWN9vi+PDZqpMmqRcvItkhZ8I9NIgBrrgCmjf/M6iHDg2G8u/ToVcv9vnsM97lWDoyjfl2GFd0+/P4yuAOna2anw8PPeQNnFEvXkQyXc6Ee/iyAV27euFb2aPfZ+NC2tEbKqbBPvvw5ZBnOOW2s9m02bY4PlTlbNXKUs+ECbFLPSIimSBnwj00iEN737NfXsmIDbdypRvHWhoyo+MIOjz3H/arX58Zx8cus1SOZw8EoLg48pozuukqIpnGnHNpeeGCggJXUlKSvBfYtAnGjKHslttgzRomWHeG1h/EE282rXYARwpx3XQVkVQysznOuYJYx+VMz/0PzsHzz3szSr/5hjqnnsq8i0fy65LWPHF8zYI30qxUv0sEi4ikUm6F+8cfe4t7zZwJrVvDtGlw6qkcAhySpJf0s0SwiEiqxZyhamb7mdm8kD+/mVnPsGPMzO61hEwTAAAGcklEQVQ1s0Vm9omZtUtekyNYvhy6dYPDDoNPP4X774f58+HUU5P+0lUuESwikiYxe+7OuS8JdnzNLB9YBjwXdlhHoFXwz5HA2ODX5Fq3DkaOhDvv9LrOvXt7KzjusEPSXzqUFhETkUwTb1mmA/CNc25J2POdgEnOuzv7gZntYGa7OeeWJ6SV4Soq4LHHoH9/b3Gvc87xVm/cd9+Yv6qRLSJSG8Qb7l2AxyI83wz4LuT7pcHnEh/uc+bA1VfDhx96ZZgpU+C443z9qka2iEht4XtVSDOrB5wJPFXdFzOz7mZWYmYlK1asqN5J1q3zNtAoLvYC3mewQ+SRLSIiuSienntHYK5z7scIP1sG7Bny/R7B57bgnBsPjAdvnHscr/2nY4+Fb7/1ut5x0sgWEakt4gn384lckgGYClxrZo/j3UhdnbR6O1Qr2CH6LFYRkVzjK9zNrCFwEnBlyHNXATjnxgGvAKcBi4B1wKUJb2lQTW+IamSLiNQGvsLdObcWaBz23LiQxw64JrFN25puiIqI+JNV2+zphqiIiD9ZFe6VN0Tz83VDVESkKlm1toxuiIqI+JNV4Q66ISoi4kdWlWVERMQfhbuISA7K6nAPBLyNrwOBdLdERCSzZF3NvZLGvIuIRJe1PXeNeRcRiS5rw11j3kVEosvasozGvIuIRJe14Q4a8y4iEk3WlmVERCQ6hbuISA5SuIuI5CCFu4hIDlK4i4jkIIW7iEgOMm+HvDS8sNkKYEk1f70JsDKBzckGuubaQddcO9TkmvdyzjWNdVDawr0mzKzEOVeQ7nakkq65dtA11w6puGaVZUREcpDCXUQkB2VruI9PdwPSQNdcO+iaa4ekX3NW1txFRKRq2dpzFxGRKmRsuJvZnmb2lpl9bmafmdl1EY4xM7vXzBaZ2Sdm1i4dbU0Un9d8YfBaF5jZ+2Z2cDramih+rjnk2MPNrMzMOqeyjYnm95rN7Hgzmxc85u1UtzORfP7bbmRmL5rZ/OAxl6ajrYliZtuY2Ych1zMowjH1zeyJYIbNNrMWCWuAcy4j/wC7Ae2Cj7cHvgJahx1zGjANMKA9MDvd7U7BNR8F7Bh83LE2XHPwZ/nAm8ArQOd0tzsFf887AJ8DzYPf75zudqfgmgcAdwYfNwV+Buqlu+01uGYDtgs+rgvMBtqHHXM1MC74uAvwRKJeP2N77s655c65ucHHvwMLgWZhh3UCJjnPB8AOZrZbipuaMH6u2Tn3vnPul+C3HwB7pLaVieXz7xngP8AzwE8pbF5S+LzmC4BnnXOlweOy+rp9XrMDtjczA7bDC/eylDY0gYK5tCb4bd3gn/CbnJ2A4uDjp4EOweuvsYwN91DBjyqH4r3zhWoGfBfy/VIiB0PWqeKaQ12G98klJ0S7ZjNrBpwNjE19q5Krir/nvwI7mtlMM5tjZl1T3bZkqeKaxwAHAN8DC4DrnHMVKW1cgplZvpnNw+uUvO6ci5phzrkyYDXQOBGvnfE7MZnZdng9tp7Oud/S3Z5U8HPNZnYCXrgfk8q2JUuMax4N9HPOVSSoU5MRYlxzHeAwoAOwLRAwsw+cc1+luJkJFeOaTwHmAX8H9gVeN7N3s/n/vXOuHDjEzHYAnjOzg5xzn6bitTO6525mdfH+IUxxzj0b4ZBlwJ4h3+8RfC5r+bhmzKwt8ADQyTm3KpXtSwYf11wAPG5m/wM6A/eb2VkpbGLC+bjmpcBrzrm1zrmVwDtAtt88j3XNl+KVopxzbhGwGNg/lW1MFufcr8BbwKlhP/ojw8ysDtAISMj/6YwN92DdaSKw0Dk3KsphU4GuwVEz7YHVzrnlKWtkgvm5ZjNrDjwLXJztvTjwd83Oub2dcy2ccy3w6pJXO+eeT2EzE8rnv+0XgGPMrI6ZNQCOxKtTZyWf11yK90kFM9sF2A/4NjUtTDwzaxrssWNm2wInAV+EHTYV+FfwcWfgTRe8u1pTmVyWORq4GFgQrFmBdze9OYBzbhzeyInTgEXAOrx3/mzm55pvxqvJ3R8sUZS57F50yc8155qY1+ycW2hmrwKfABXAA6n6OJ8kfv6eBwMPm9kCvJEm/YKfWrLVbkCxmeXjdaSfdM69ZGa3ASXOual4b3iTzWwR3g3kLol6cc1QFRHJQRlblhERkepTuIuI5CCFu4hIDlK4i4jkIIW7iEgOUriLiOQghbuISA5SuIuI5KD/D63oLHFYev1yAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nds = 100\n", "x = np.linspace(2, 3, num=nds).reshape(nds, 1)\n", "theta_data = np.array([2.0, 3.0])\n", "y = linear_model(x, theta_data) + 0.1*np.random.standard_normal(x.shape)\n", "plt.plot(x, y, '.b');\n", "plt.plot(x, linear_model(x, theta_data), '-r');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Initialize MCMC Object\n", "- initialize data structure\n", "- define simulation options" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Initialize MCMC object\n", "mcstat = MCMC()\n", "mcstat.data.add_data_set(x, y)\n", "mcstat.simulation_options.define_simulation_options(\n", " nsimu=int(5.0e3),\n", " updatesigma=True,\n", " method='dram')\n", "mcstat.model_settings.define_model_settings(sos_function=ssfun)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Specify Sample\n", "When defining the model parameters, we can turn the sampling on or off by specifying the value of the `sample` argument.\n", "- `sample = True`: Parameter will be treated as a random variable and be included in the sampling chain. This is the default behavior.\n", "- `sample = False`: Parameter will be held constant at the initial value defined by the user." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "mcstat.parameters.add_model_parameter(\n", " name='m',\n", " theta0=2.,\n", " sample=False)\n", "mcstat.parameters.add_model_parameter(\n", " name='b',\n", " theta0=5.,\n", " minimum=-10,\n", " maximum=100,\n", " sample=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Run Simulation\n", "Observe that the only parameters displayed are the ones where `sample = True`." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Sampling these parameters:\n", " name start [ min, max] N( mu, sigma^2)\n", " b: 5.00 [ -10.00, 100.00] N( 0.00e+00, inf)\n", " [-----------------100%-----------------] 5000 of 5000 complete in 1.6 sec\n", "---------------------\n", "name : mean std MC_err tau geweke\n", "b : 3.0184 0.0087 0.0004 5.9851 0.9995\n", "---------------------\n" ] } ], "source": [ "mcstat.run_simulation()\n", "# Extract results\n", "results = mcstat.simulation_results.results\n", "chain = results['chain']\n", "s2chain = results['s2chain']\n", "sschain = results['sschain']\n", "names = results['names']\n", "\n", "# define burnin\n", "burnin = 2000\n", "# display chain statistics\n", "mcstat.chainstats(chain[burnin:, :], results)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# MCMC Plots\n", "The plots only show the results of the sampled parameters `b`." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAEJCAYAAAAw1SpDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcFOW97/HPjxlmQGAGkV0dFkE22URkBAEVXCAaURljjvFcCcmJmuV4uPdEzU1OjjnmqDFRY8R44tW4ZiGBgJCgbKKIbAKygwoo6MgMi7OwzP67f1T32IzMMNNT3U919+/9etULpqu6+lfSfuepp6qeR1QVY4xxqYXrAowxxoLIGOOcBZExxjkLImOMcxZExhjnLIiMMc5ZEBljnLMgMsY4Z0FkjHHOgsgY45wFkTHGuXTXBbgkIgJ0B0pd12JMkinVJjzImtJBhBdCn7guwpgklA2UNHbjVA+iUoD9+/eTlZXluhZjEl5JSQnnnntuk9+X6kEEQFZWlgWRMQ5ZZ7UxxjkLImOMcxZExhjnLIiMMc5ZZ7VJClu2bGHDhg0cOnSIyy67jBEjRrguyTSBBZFJaKWlpdx777089dRTta+lp6fz17/+leuvv95hZaYp7NTMJKxdu3YxbNiw2hAaP348o0ePpqqqiry8PF577TXHFZrGsiAyCWnVqlWMHj2aPXv20KNHD5YsWcLy5ct58803ycvLo7Kykttvv53KykrXpZpGsCAyCeedd95h4sSJHDlyhIsvvph169YxYcIEwDste/nll+nSpQsFBQX84x//cFytaQwLIpNQNmzYwOTJkzl+/DhXXXUVy5Yto1OnTidtk5GRwW233QbA888/76BK01QWRCZh7Nu3j0mTJlFcXMzYsWP529/+Rps2bU657e233w7AggULKCwsjGOVJhoWRCYhHD9+nClTplBYWMjQoUNZsGABZ5xxRr3bDxo0iJEjR1JVVcUf/vCHOFZqomFBZAJPVfnWt77Fxo0b6dixI/PmzWvUQ8rTpk0D4OWXX451iaaZLIhM4M2cOZM//vGPpKenM3v2bHr06NGo94XvI9q4cSMlJY0eGsc4YEFkAm3t2rXMmDEDgF/84heMGzeu0e/t3r07PXv2pKamhrVr18aqROMDCyITWIcOHaq9J+jGG2/k7rvvbvI+Ro8eDcDKlSv9Ls/4yILIBFJ1dTW33nor+/bto0+fPjz33HN4Q4w3TTiI3nnnHb9LND6yIDKBdP/997No0SJat27NnDlzyM7Ojmo/Y8aMAWD16tVUV1f7WaLxkQWRCZzXX3+dBx54AIBnnnmGwYMHR72vCy64gLZt21JSUsK2bdv8KtH4zILIBMqnn37KN77xDVSVO+64g1tvvbVZ+0tPTyc3Nxew07MgsyAygaGqTJs2jUOHDjF06FAee+wxX/ZrHdbBZ0FkAmP27NksXryYzMxMZs2aRatWrXzZ7yWXXALAunXrfNmf8Z8FkQmEY8eO1d4vdM8993D++ef7tu8hQ4YA8MEHH3DixAnf9mv8kxBBJCL3ioiKyOMRr7USkZkiclhEjorIbBHp4rJOE72HH36Y/fv306NHD+655x5f992tWzc6dOhATU0NO3bs8HXfxh+BDyIRGQl8B9hcZ9VjwHVAHjAeb/roOfGtzvjh6NGjPPHEEwD88pe/bPBh1miISO2Vty1btvi6b+OPQAeRiLQFXgG+DXwe8Xo2MB2YoarLVHU9MA0YLSK5Too1UXvppZcoLi6mT58+3HjjjTH5DAuiYAt0EAEzgb+r6pI6r48AWgK1r6vqTmAfcEl9OxORTBHJCi9AuxjUbJqgpqamtjX0/e9/nxYtYvOVDPcTWRAFU2Bn8RCRW4ALgZGnWN0VqFDVojqvF4TW1ec+4Kf+VGj8sHjxYnbu3Em7du1qBzOLBWsRBVsgW0Qici7wa+BWVS3zcdcPAtkRyzk+7ttE4emnnwa8sYMaM8ZQtAYNGgTAZ599xuHDh2P2OSY6gQwivFOvzsAGEakSkSq8DukfhP5eAGSISPs67+sCHKhvp6parqol4QUojVH9phFKSkpqB7efPn16TD+rXbt29OrVC7BWURAFNYiWAoOBYRHLu3gd1+G/VwITwm8QkX5ADrAq3sWa6Lz66qtUVFTQv3//Zj1P1lh2ehZcgQwiVS1V1a2RC3AMOBz6uRh4FnhURC4XkRHA74FVqrraZe2m8WbNmgXAzTffHNUQH01lQRRcge2sboR/A2qA2UAm8Dpwl9OKTKMVFxfz+uuvA14QxYMFUXAlTBCp6mV1fi4DvhtaTIIJn5YNHDiwtiM51gYOHAjAjh07UNW4tMJM4wTy1MwkvzlzvJvg8/Ly4vaZ559/Pi1atKC4uJgDB+q9pmEcsCAycVdeXs6SJd69qF/96lfj9rmZmZmcd955APbMWcBYEJm4e/vttzl69Chdu3Zl2LBhcf3sAQMGALB9+/a4fq5pmAWRibvwvUOTJk2K2SMd9QkHkbWIgsWCyMRdOIgmT54c988Od1hbiyhYLIhMXO3Zs4edO3eSlpbGlVdeGffPtxZRMFkQmbgKt4YuvfTSqKcIao7+/fsDUFBQwJEjR+L++ebULIhMXC1atAjw+odcaNeuHeeeey5graIgsSAycVNZWcny5csBnJyWhdnpWfBYEJm4WbduHaWlpXTo0CHul+0jWYd18FgQmbhZunQpABMmTIj7ZftIdi9R8FgQmbgJ3009YcKE02wZWxdccAGATUEdIBZEJi6OHj3KqlXeUFETJ050Wkv4IdtPPvmEoqK6ow0bFyyITFysWLGCyspKevbsSe/evZ3Wkp2dXXvlzFpFwWBBZOJi2bJlgHdaFoThN8KnZ1u3bnVciQELIhMnb731FgCXXXaZ20JCLIiCxYLIxNzRo0dZv349AOPGjXNcjceCKFgsiEzMrV69murqanr06EFOTo7rcoAvgmjLli2oquNqjAWRibnwadnYsWMdV/KFAQMGICIcPnyYwsJC1+WkPAsiE3PhIArKaRlA69at6dOnD2CnZ0FgQWRiqry8nDVr1gDBahGB9RMFiQWRial3332XsrIyOnXqRL9+/VyXc5LIfiLjlgWRiakVK1YAXmsoCPcPRRo6dCgAGzdudFyJsSAyMbV6tTfx7pgxYxxX8mUXXXQR4LWIysrKHFeT2iyITMyoam0Q5ebmOq7my3JycujYsSOVlZVs3rzZdTkpzYLIxMzHH39MQUEBLVu2ZPjw4a7L+RIRYeTIkYDXl2XcsSAyMRNuDQ0bNozWrVs7rubUwqdnFkRuWRCZmAkH0ahRoxxXUj8LomCwIDIxE75/KIj9Q2HhINq2bRvHjx93XE3qCmwQicidIrJZREpCyyoRmRSxvpWIzBSRwyJyVERmi0gXlzWbL5SXl7NhwwYg2EHUvXt3unXrRk1NDe+9957rclKWL0EkIi1F5FwR6SciHfzYJ/AJcC8wArgIWAbME5FBofWPAdcBecB4oDswx6fPNs303nvvUVFRQceOHZ0PhHY6dnrmXtRBJCLtQq2WN4ES4CNgB3BQRD4WkWdEZGS0+1fV+ar6D1X9QFXfV9X/CxwFckUkG5gOzFDVZaq6HpgGjBaR4P76TSGR/UNBu5GxrnAQhYeyNfEXVRCJyAy84JkGLAGmAMOA84FLgPuBdGCRiLwmIn2bU6SIpInILUAbYBVeK6ll6LMBUNWdwL7Q59e3n0wRyQovQLvm1GXqF25dXHzxxY4rOb3ww7jLly+3IUEcSY/yfSOBcapa34C/a4HnROQOvLAaC3zQ1A8RkcF4wdMKrzV0g6puF5FhQIWq1h35vADo2sAu7wN+2tQ6TNOFgyh8n06Q5ebmkpmZyYEDB9i1a1fttNQmfqJqEanq1xsIocjtylX1aVV9LprPAXbhtbRGAb8FXhCRgVHuC+BBIDtiOacZ+zL1KCkpYdeuXQCMGDHCcTWn16pVK0aPHg3AG2+84bia1BTtqVlLEblSRMaKyFl+FxWmqhWq+qGqrlfV+4BNwL8CB4AMEWlf5y1dQuvq21+5qpaEF6A0VrWnsg0bNqCq5OTk0LlzZ9flNMrll18OWBC5Em1n9Ry8q1V/A1aLyKci8pp/ZdWrBZAJrAcqgdqZ+kSkH5CDdypnHAqfloU7gRNBeFB/6ydyI9o+ohxVvU5ELlbVYSLyXaCHn4WJyIPAQrwO6HbAPwGXAVerarGIPAs8KiJH8K7a/QZYpaqr/azDNN26deuAxOgfCrv44otp3bo1Bw8eZPv27bWTMJr4iLZFFB4zoUJEMlR1JnCpTzWFdQZexOsnWorXQX61qi4Orf83YAEwG3gL75TsRp9rMFFIxBZRZmZm7VAldnoWf9EG0ROhGxdnA0+LyHSgo39lgapOV9Weqpqpqp1VdWJECKGqZar6XVXtoKptVPVGVa23f8jEx5EjR9izZw+QGB3VkSZM8M70X3stHr0MJlLUfUSqekRVH8ZrjfQHpvpXlklU4fnL+vTpw5lnnum4mqa59tprAVi6dCnHjh1zXE1qiTaI3g7/RVWfV9V/Byr8KckksvBpWaK1hgAGDRpEz549KSsrY8mSJad/g/FNk4JIRK4TkXuAtiJybp3Vf/avLJOoErF/KExEuO666wCYP3++42pSS1NbRFvx7r3pCLwoIrtF5C0RmYV3Od2kuEQOIqA2iBYsWEBNTY3jalKHRHPPhIiMU9W3Qn8/G+/S/dbQTYIJI/S8WXFxcTFZWVmuy0l4Bw8erL2BsaioiOzsbMcVNV14xIDS0lLWrFmTEM/KBUlJSUn43z27KXkQbR/RqvCd1UCZqr6TaCFk/BfuqD7//PMTMoQAMjIyuPrqqwF49dVXHVeTOhLtzmoTYOEgSsSO6kjXX389APPmzXNcSeqINohyVPVfgE9UtS/w34DNx5LiEr1/KOwrX/kKaWlpbN26ld27d7suJyUE+c5qk2CSpUV05plnMn78eMBaRfES2DurTWIpLCxk//79iEgg5zBrqilTpgAwd+5cx5WkhqbeR3QegKq+YndWm0iRHdXJcAUy3E+0cuVKCgsLHVeT/JraInpaRPaGZtT4n9BT97uBB1TV+ohSWCLfUX0qOTk5DB8+nJqaGhYsWOC6nKTXpCBS1StVtRcwH+/p+LOBHwOHReTDGNRnEkQiDQ3bWOGbGxcuXOi4kuQXbR/Rzap6g6r+SFWvBiYDK32syySYZLliFmnSJG8avSVLllBVVeW4muQW9VWzyLGjVXURcIE/JZlEk5+fT35+Pi1atGDYsGGuy/HNyJEjOfPMMykqKmLt2rWuy0lq0QbRdODPIvIbEZkuIo8DNr5migp3VA8YMIC2bds6rsY/aWlpXHXVVYCNURRr0c7isQ1vbrEVQE/gY2BSQ+8xySsZT8vCrrnmGsCCKNaiHbMaVa0AZoUWk8KSOYjCz529++67HDx4kE6dOjmuKDk1Z8rpH4vIZBHp4mdBJrGoalIHUbdu3Rg6dCiqyuLFi0//BhOVqIMI+BneZfz80EOv80XkP0Xkq6GhQUwK+PTTTyksLCQtLY2hQ4e6LicmJk6cCMCbb77puJLk1ZwgWgd8CjwAPAQcwptF46/APhE5ICL/aH6JJsjCraFBgwbRunVrx9XExrhx4wB46623HFeSvJrTRzRKRG7He/J+HTBDVXeLSCbeNNEXAon/0JFpUPiKWTKeloVdeqn3PPfOnTspLCxMmNlrE0lzWkSo6vPA+Xhzj20ITYqYpqprVPW3oaFCTBLbsGEDABdeeKHjSmKnQ4cODB48GIAVK1Y4riY5NSuIAFT1qKr+ELgI76bGD0Xkn5tdmQk8VU2aoT9Ox07PYqvZQQQgIul4c9L/EfgE+H1omBCTxPLz8ykoKKBFixYMGTLEdTkxZUEUW1H3EYnIvcDg0NIfb7C0zcBa4H+AYj8KNMEVPi0bOHAgZ5xxhuNqYmvs2LEAbNq0iaKiItq3b++4ouTSnBbRfwOX4A2OdoGqZqnqpar6PVV9VlWr/SnRBFX4tCyZ+4fCunXrRt++fVFV3n777dO/wTRJc4JoBXAW8FNgvYisDD17Nk1EhopImj8lmqAKt4iSvX8oLHz1bNWqVY4rST5RB5GqjlfVbKAf8G28YUAGAL8CNgJHRcQeWU5iqdQiAsjNzQVgzZo1jitJPn5cNftAVf+kqj9U1Ymq2gE4D/hnIOoJxEXkPhFZJyKlIlIoInNFpF+dbVqJyEwROSwiR0Vktj1yEh8HDhwgPz8fEUmqoT8aEg6itWvXUl1tPQ9+8uWqWV2quldV/6KqP2rGbsYDM4Fc4EqgJbBIRNpEbPMYcB3eHGvjge54c66ZGNu4cSMA/fr1S6qhPxoyaNAg2rRpQ2lpKTt27HBdTlKJSRD5QVWvUdXnVXWbqm4Cbgdy8IYfQUSy8cZFmqGqy1R1PTANGC0iua7qThWbNm0CSJnWEHjjE4WnoF69erXjapJLVEEUGkB/TxTLD5pRa3gO4yOhP0fgtZJqT/9UdSewD+9q3qnqzhSRrPACtGtGPSktHETJ+qBrfcKnZxZE/or2PqLbo3zfR9G8SURaAI8DK1V1a+jlrkCFqhbV2bwgtO5U7sO7ymeayYLIgshPUQWRqsZ7PISZeI+PNHc22QeBRyN+bod3J7hpghMnTrBr1y4g9YJo1KhRAGzfvp2SkpKkmMMtCHztI5IQn/f5JHAtcLmqRobGASBDROre4toltO5LVLVcVUvCC1DqZ62pYtu2bdTU1NCxY0e6devmupy46tKlCz179kRVbUB9H/n1rNl0EdmK95hHmYhsFZFvNXOfEgqhG4ArVHVvnU3WA5XAhIj39MPr0LY7zmIo8rTM5987CcHuJ/Jfs4NIRH4G/BpvtMa80DIfeCy0LlozgW8A/wSUikjX0NIaQFWLgWeBR0XkchEZAfweWKWqdgIfQ6naPxQWPj2zIPJP1A+9RrgT+Laq/jHitVdFZDPwG+A/mrFfgOV1Xp8GPB/6+78BNXjPu2UCrwN3Rfl5ppFSPYgiO6xVNSVbhX7zI4haAu+e4vX1zdm/qp72X1dVy4DvhhYTB6qa8kE0bNgwWrZsycGDB/noo4/o1auX65ISnh99RC/xResl0r8Ar/iwfxMg+/bto7i4mJYtWzJgwADX5TjRqlWr2hs57TK+P6K9ofHR8II3w+u3Qh3U/y+0bMF7ELbGz2KNe1u2bAG8RzsyMjIcV+OOdVj7K9oW0fCIZTDeadhBvIddz8Ob0WMDMMiHGk2AbNu2DaB2DOdUZR3W/or2hsbL/S7EJIatW70b2y+44ALHlbgVbhFt2LCB8vJyMjMzHVeU2KI9Nctp4vY24WKSsCDy9O7dm44dO1JRUcF7773nupyEF+2p2ToR+R8RGVnfBiKSLSLfDt3oeFOUn2MCpKqqqnb4i1QPIhGpbRWtXLnScTWJL9ogGggcAxaHZnT9u4g8Exoq9mUR2QAUAt8EfqiqT/hVsHFn9+7dlJeXc8YZZ9CzZ0/X5TgXHlDf5jprvqiCSFUPq+oMoBvePTwfAB2BvqFNXgFGqOolqmrTTieJ8GnZwIEDadEisENZxU04iN5++21U1XE1ia1ZNzSq6onQ0Kz/rqqVPtVkAsr6h042YsQIWrduzaFDh9i5c2fK3lflBz9+rf0GuFZEOtVdISKnHKDMJKbwpXsLIk9GRkbtZXw7PWseP4JIgL8AB0TkMxFZJCK/FJHvAQt82L8JCGsRfVl4BlgLoubx41kzgN5AJ2AIMBRvGNebAbv/PUmUl5fz/vvvAxZEkcL9RDYVdfP4FURlocHr1/u0PxMwu3btorq6muzsbLp37+66nMDIzc0lLS2Nffv2sW/fPnJymnSLnQnx69LHABFp6dO+TABFnpbZsBdfaNu2LRdddBEAS5ZEPY1fyvMriJbhzey6RUT+ICL3ishkETnHp/0bx6x/qH7XXHMNAK+99prjShKXH0G0E29g+2uA3+Hd6HgDMAv42If9mwAIB1GqP+x6KuEgWrx4MVVVVY6rSUzN7iNS1YGhv+4A3gi/HhpE/7zm7t8Eg7WI6jdy5Eg6dOjAkSNHWLNmDWPGjHFdUsKJ2e2x6vkwVvs38XP06FH27vXmLhg0yEZ2qSstLY2rrroKgIULFzquJjHZffrmtLZv3w5A165d6dixo+NqgmnSpEmA9RNFy4LInJadlp3e1VdfDcD69espKChwXE3isSAyp2VBdHpdunSpvYw/d+5cx9UkHgsic1oWRI0zdepUAP7yl784riTxWBCZ07Igapy8vDwAli9fzsGDBx1Xk1gsiEyDDh8+zGeffQZ44xCZ+vXu3ZsLL7yQ6upqOz1rIgsi06DweMy9e/emXbt2jqsJPjs9i44FkWnQxo0bARg+fLjjShJDOIiWLVvGoUOHHFeTOCyITIMsiJqmb9++DB8+nOrqaubMmeO6nIRhQWQaZEHUdLfccgsAf/rTnxxXkjgsiEy9jh8/zq5duwALoqa4+eabAe/qWbij3zQssEEkIuNEZL6I5IuIisiUOutFRH4WGp72hIgsEZG+9e3PNN3mzZupqamhS5cudOvWzXU5CaNnz55ccsklqKp1WjdSYIMIaANswpuu6FR+CPwAuAMYhTf8yOsi0io+5SW/8BUzaw01nZ2eNU1gg0hVF6rqj1X1b3XXhYYYuRt4QFXnqepm4J+B7sCUutub6Fj/UPTy8vIQEVatWsXHH9uwXKcT2CA6jV5AV6B2bE5VLQbWAPVOYSQimSKSFV4AuzGmARZE0evWrVvtwPp29ez0EjWIuob+rPuYc0HEulO5DyiOWD7xv7TkUFVVxZYtWwALomjddNNNAMyePdtxJcGXqEEUrQeB7IjFxtSux3vvvUdZWRnt27end+/erstJSDfeeCMAK1euJD8/33E1wZaoQXQg9GeXOq93iVj3Japarqol4QUojVWBiS48YeCYMWNsnvsonXPOOeTm5gLwt799qavTREjUb9hevMCZEH4h1OczCljlqqhk8vbbbwNfTCBoohN+5MNOzxoW2CASkbYiMkxEhoVe6hX6OUdVFXgc+LGIfFVEBgMvAvmAPfbcTKpa2yKyIGqecD/Rm2++ac+eNSCwQQRcBGwMLQCPhv7+s9DPvwB+gzeF0TqgLXCNqpbFuc6k8/7773Pw4EEyMzMZMWKE63ISWs+ePRk+fDg1NTXMnz/fdTmBFdggUtXlqiqnWG4PrVdV/Q9V7aqqrVR1oqq+77jspBBuDY0aNYrMzEzH1SS+KVO8W9tsjKL6BTaIjDt2WuavcBAtWrSIY8eOOa4mmCyIzJdYEPlr8ODB9OrVi7KyMhYtWuS6nECyIDInef/999m7dy/p6elcckm9N6mbJhAROz07DQsic5LwZeYrrriCrKwsx9Ukj3AQzZ8/n8rKSsfVBI8FkTlJOIjC978Yf4wZM4bOnTvz+eefs2zZMtflBI4Fkam1d+9e1q9fT4sWLWp/gxt/pKWl1T7yMWvWLMfVBI8FkakVfkp8/PjxdOrUyXE1ySc8cuPcuXPt9KwOCyJT669//Svwxd3Axl/jxo2jc+fOHDlyxE7P6rAgMoA39tDq1atp0aIFN9xwg+tyklJaWlptyNvp2cksiAwADz/8MABf+9rX6N69u+Nqkld4Wuo5c+Zw4sQJx9UEhwWR4cMPP6wd5P3ee+91XE1yGzduHD169KCoqMhaRREsiAyPPPIINTU1TJ48mSFDhrguJ6mlpaXxne98B4Df/va3jqsJDvFG1EhNoTGMiouLi1P25r2NGzcycuRIqqureeutt+yxjjgoLCzknHPOobKykvXr13PhhRe6Lsk3JSUlZGdnA2SHBh9sFGsRpbCKigqmTZtGdXU1eXl5FkJx0rlz59pO66eeespxNcFgQZTCHnroITZt2sRZZ53Fk08+6bqclHLXXXcB8OKLL9bOppvKLIhS1MKFC7n//vsBeOKJJ+jcubPjilLLpZdeyuTJk6msrOR73/seqdxFAhZEKWnbtm187Wtfo6amhunTp/P1r3/ddUkpR0R44oknyMzMZMmSJSk/NbUFUYrZs2cPkyZNorS0lHHjxvHUU0/hTZxr4u28887jvvvuA2D69OkpPVaRBVEK2bt3L5dffjn79++nf//+zJ49m4yMDNdlpbR77rmHK664gqNHj/KVr3yFn/zkJ+zYsSPlTtXs8n2KXL5fsWIFU6dOpbCwkH79+vHGG2/QrVs312UZvKuX3/zmN3nllVdqX2vTpg05OTkMHz6csWPHMnXqVDp27OiwysaJ9vI9qpqyC5AFaHFxsSar48eP6wMPPKDp6ekK6JAhQzQ/P991WaaO6upqfemll3Ty5MnasmVLBU5asrKy9MEHH9QTJ064LrVBxcXF4ZqztAn/L1qLKElbRB988AGzZ89m5syZfPLJJ4D3HNmzzz5LmzZtHFdnGnLixAn279/P3r17Wb16NXPmzGHz5s0AjBw5knnz5gW2NRtti8iCKAmCSFV5//33Wbp0KcuXL2fVqlW14QOQk5PDz3/+c2699VbrmE5ANTU1vPLKK9x9990cOXKEs88+m7///e8MHTrUdWlfYkEUhUQPooKCAl544QVeeukltm7detK69PR0Lr/8cvLy8rjtttto1aqVoyqNX3bv3s21117Lzp07ad++PQsXLiQ3N9d1WSexIIpCIgZRdXU1S5cu5fe//z2zZ8+uHekvIyODMWPGcMUVVzBmzBhGjhxJ27ZtHVdr/FZUVMS1117LypUradOmDXPnzmXixImuy6plQRSFRAiikpIStm7dyqZNm1ixYgXLli2joKCgdv2oUaOYPn06eXl5tG/f3mGlJl6OHTvGlClTWLJkCenp6Tz33HPcdtttrssC7KpZ0lw1q6qq0sWLF+u3v/1tHTBggIrIl66gdOjQQe+8805dv36963KNI2VlZXrLLbfUfidmzJihZWVlrsuyq2bRCFKL6MSJEzzzzDM88sgjJ3U0A5x99tkMHjyY3Nxcxo4dy6WXXmo3Ihpqamr40Y9+VDu65rBhw3j66acZNWqUs5qsRZSgLaLKykr93e9+p927dz+pxXPHHXfoq6++qgUFBc5qM4lh3ry1E8S5AAAJF0lEQVR5etZZZ9V+f2666SZduXKl1tTUxL2WaFtEzsOguQvwXeAjoAxYA1zchPc6C6KKigp98cUXtX///rVfoJycHH366acD0cQ2iSU/P1+nTZt20qn8wIED9cc//rGuXLkybt+plDw1E5GvAS8Cd+CF0N1AHtBPVQsb8f64npp9/vnnbNy4kfnz5zNr1izy8/MBOOuss/jJT37CHXfcQWZmZszrMMlr69at/OpXv+LPf/7zSYPzZ2RkMGTIEAYPHsyAAQPo27cvffr0oWfPnr5eXU3Jq2YisgZYp6rfC/3cAtgP/EZVH2rE+08bRGvWrGHHjh3U1NRQU1NDVVUVVVVVVFRUUFFRQVlZGSdOnKhdysrKKCsro7y8vHZdUVERBQUFHDx48KR9d+rUiRkzZnDXXXc576MyyaWoqIhXX32VBQsW8MYbb3Do0KF6tz3rrLM455xz6N69O507d6ZDhw5kZWXRpk0bMjMzycjIID09naysrNpJIuuTckEkIhnAcWCqqs6NeP0FoL2qXn+K92QCkU2OdsAnDQXR97//fV9HL+zVqxejR4/m5ptv5uqrr7YWkIk5VWXv3r1s2LCBbdu2sXPnTj744AP27NnD559/3uj99OjRg48++qjBbaINovRGVxE8HYE0oKDO6wVA/3recx/w06Z8yKBBg5g8eTItWrRAREhPTyc9PZ2MjAwyMjJo3bo1rVq1OunPzMxMMjMzadWqFWeccQZZWVl06dKFc889lzPPPLPpR2pMM4gIvXv3pnfv3kydOvWkdUVFRezfv5/9+/dTUFBAQUEBRUVFFBcX17byw2cBsXz6P5FbRN2BT4HRqroq4vVfAONV9UvXMKNpERljGi8VW0SHgGqgS53XuwAHTvUGVS0HysM/2wOgxgRDwo7QqKoVwHpgQvi1UGf1BGBVfe8zxgRPIreIAB4FXhCRd4G1eJfv2wC/d1qVMaZJEjqIVPXPItIJ+BnQFXgPuEZV63ZgG2MCLKGDCEBVnwRsdkBjEljC9hEZY5JHwreI/FBS0viHhI0x9Yv2/6WEvY/IDyJyNvDJaTc0xjRVajzi4QfxbiTqDpSeYnU7vJA6p571icqOK7Ek6nGVahPCJaVPzUL/oT491bqImx1Lm5LsQWfHlViS9bjqss5qY4xzFkTGGOcsiOpXDtxPxLNpScKOK7Ek63GdJKU7q40xwWAtImOMcxZExhjnLIiMMc5ZEBljnEv6IBKR+0RknYiUikihiMwVkX6neU9LEfkPEdktImUisklErqmzTZqI/JeI7BWRE6FtfyJxHPZRRO4Ukc0iUhJaVonIpNO8J09EdoaOa4uITK6zXkTkZyLyWei4lohI39geyZdq9PW4Qv+eD4dePyYi+SLyYmi44biJxb9XnW2fFhEVkbv9rz7GmjIJWiIuwGvA7cAgYCjwd+BjoE0D73kY747ryUBv4E7gBDA8Ypsf4Q1X+xWgJzAV7xb8H8Tx2K4L1dgXOB/4OVABDKpn+9FAFfDvwADgv0LbXxCxzT1AEXA9MASYB+wBWiXqcQHZwGLgZqAfkIs3D967cf4u+v7vFbHtDXjjcX0K3B3P4/Llv43rAuJ+wNAJbybKcQ1skw98t85rs4GXI35eADzb0DaOju8IML2edX8GFtR5bTXwdOjvAnwG/J+I9dl4s+jekqjHVc97Roa+BzmJflxA+OHtQXizHidcECX9qdkpZIf+PNLANpl4//NFOgFcGvHzO8AEETkfQESGhtYv9KnOJgmdKt6CN1RufWN2XwIsqfPa66HXAXrhjXRZu42qFuO1Hi7BAZ+O61Sy8YKoqNlFRsGv4wqN0/4S8IiqbotFrfGQUg+9hv7RHgdWqurWBjZ9HZghIm8Bu/EG5L8Rbx61sIeALGCniFSH1v1fVX0lJsXXQ0QG432RWwFHgRtUdXs9m3fl1PPAdY1Yz2m2iQufj6vuvlvhnX7/UeP8IGkMjusevNO3J3wuNa5SKoiAmcAFnNyyOZV/BZ4BduL91tyNNyD/NyO2uRm4FfgnYBswDHhcRPJV9QWf627IrtBnZ+P1U70gIuMb+HInipgcl4i0BGbhnYbe2ewqm8634xKREXjf1Qs1dI6WqFLm1ExEngSuBS5X1QYHQ1PVg6o6Ba/Z3ANv5tijeJ22YY8AD6nqn1R1i6q+BDyGN5ts3Khqhap+qKrrVfU+YBPel/NUDtDwPHAHIl6rb5u48Pm4gJNCqAdwZbxbQ+D7cY0FOgP7RKRKRKrwju1XIvKR/9XHTtIHUehy9JN4VxWuUNW9jX2vqpap6qd4Lceb8K4ghZ0B1NR5SzXu/5u24OTZbCOtImIeuJAr+aKPYi/elzxyrrgsYBTu54prznFFhlBfYKKqHo5FkVFoznG9hHdlc1jEko/3S/Jq3yuNJde95bFegKfwOiTH451bh5fWEdu8CDwY8fMovD6h3ni/dZbitYbaR2zzPN6VivDl+xuAg8DDcTy2B4Fxoc8fHPq5Bu+3/amOazRQCfxvvFbef3Lqy/efA18N7XMu8b987+txAS3xfonsx7uFI/J7kJGox1XPZ3xEAl41c15AHP7xtZ7l9ohtlgPPR/w8HtiOd+XsUOgL0r3OftvhdXx/jHdFbTfwQJy/2M+GvnjlQCHeFZYr6zuu0Gt5eP0U5cBWYHKd9YI3T9yB0PEvAc6P87+Zr8cV+h+/vu/BZYl6XPV8RkIGkQ0DYoxxznV/hjHGWBAZY9yzIDLGOGdBZIxxzoLIGOOcBZExxjkLImOMcxZExhjnLIhMYInIchF53HUdJvYsiIwxzlkQGWOcsyAyQZcuIk+KSLGIHArNnBK3mVJMfFgQmaD7X3hDoV6MN4DYDOBbTisyvrOn701gichyvBEIB2l4jBKRh4CvqupAl7UZf1mLyATdaj35t+UqoK+IpNX3BpN4LIiMMc5ZEJmgG1Xn51zgA1WtdlGMiQ0LIhN0OSLyqIj0E5GvA98Hfu26KOOvVJvXzCSeF4HWwFq8WVJ+DfzOaUXGd3bVzBjjnJ2aGWOcsyAyxjhnQWSMcc6CyBjjnAWRMcY5CyJjjHMWRMYY5yyIjDHOWRAZY5yzIDLGOGdBZIxx7v8DxQHCq2Xi878AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAASMAAAEJCAYAAADfF0F9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnX+4H0V971+fcxRokcRrVahEoEpQBPxRjEhsAaF5KClYvCU00ssNt7H8LAVLlAKXH9pCrokEbisELClIb3/k9GmrDV6L5UmoPk3wKFZ+lAsFNXIOIqBIAigJkLl/zC7ZM2dmdmZ/fL/7/Z55P88833N2Z3dndnfe+/k1nxGlFAkJCQn9xki/G5CQkJAAiYwSEhI6gkRGCQkJnUAio4SEhE4gkVFCQkInkMgoISGhE0hklJCQ0AkkMkpISOgEEhklJCR0AomMEhISOoFERgkJCZ3Aq/rdgC5CRAR4E/Bsv9uSkDAA2AP4gao50TWRkR1vAib73YiEhAHCHOCxOidIZGTHswATExPMmjWr321JSOgstm7dypvf/GZoQItIZOTBrFmzEhklJPQIyYCdkJDQCSQySkhI6AQSGSUkJHQCiYwSEgYYk5OwYYP+HXQkMkpIGFCsWQP77gtHH61/16zpd4vqQVJC/ukQkVnAli1btiRvWkInMTmpCWjHjp3bRkdh82aYM6d37di6dSuzZ88GmK2U2lrnXEkySkgYQDz88FQiAnj5ZXjkkf60pwkkMkpIGEDMnQsjxugdHYX99+9Pe5pAIqOE1jBMxtWuYc4c+NznNAGB/r3xxt6qaE0jkVFCKwgxriayqoelS7WNaMMG/bt0ab9bVA/JgG1BMmDXQ4hxdc0aOP10XWdkRH/lB30wzUQkA3ZCp1FmXJ2c3ElEoH/POCNJSDMdiYwSGkeZcXUYPUEJ9ZHIaIDwjW/AqlX6t8soM64OoycooT76SkYicpaI3CsiW7OySUSOKzlmkYg8KCIviMh9IrLQU/cGEVEicn7zre8tTjsN3vc+uOAC/Xvaaf1ukR8+4+oweoIS6qOvBmwROQF4GXgYEGAJ8HHgPUqp/7DUnw98FbgIuA04BbgQ+GWl1P1G3Q8DlwNvAFYqpa6NaFenDNjf+IYmIBPj4zBvXu/b0xQmJ7Vqtvvu8L3v6W3z5ydSahOTk1pNnju3mfvcpAEbpVSnCvA0sNSxby1wm7HtLuAGY9ve6LSxBwGbgfNLrrkrMKtQ9gbUli1bVBdw9dVKwfRyzTX9bll93HSTUiI7+ySityU0j5tuUmpkRN/nkZFm7vOWLVsUoIBZqubY74zNSERGRWQxsDuwyVHtcOAOY9vt2fb8PCPAX6KloWnSlQMXAVsKpVN+nV/9Vfv2D3ygt+1oGpOT8Hu/p2koh1LJs9YGBsGD2XcyEpFDROQ5YBtwA/BhpdQDjup7AU8Y257Itue4EHgJ+NOIZiwHZhdKpxSFefNgyZKp25YsGWwVDbS6YLMSJM9a8xgED2YXcmA/BLwbTQInAZ8XkSM9hOSEiBwKnIe2IQUbw5RS29BkmJ8n9tKt45Zb4Jxz4N/+TUtEg05EoO0WItMJKXnWmkfuwTQDUbt0n/suGSmltiulHlFK3a2Uugi4B00oNvwQ2NPYtme2HeBXgTcCj4rISyLyErAvcLWIbG6+9b3FvHlw/vnDQUSgDah//ueakHKMjCTPWhsYBA9m56aDiMh64FGl1GmWfWuBn1dKnVDYthG4Vyl1poj8AvCLxmG3o21INyulHgpsQ6e8acOOyUnYlFkJDz+8WwNk2JB7MPffv3vetL6qaSKyHPgy8Ch6VcpTgKOAY7P9twKPZRITwP8G/lVELgC+BCwG3gucDqCU+jHwY+MaLwI/DCWihN5jzhxYtKjfrZgZmDOnu2Tfb5vRG4Fb0dLMFuBe4Fil1L9k+/cBXtFylVIbReQU4E+Aq9DxSScqI8YoISFhOpqOMWoanVPTuoAuqGldf3ESBgttZUlIs/aHHMOWaD2hvxiEGCNIZNQ5DMqLkzA4GIQYI0hk1DkMyovTdaQskjsxKFkSEhl1DIPy4nQZvVJze0F4TVxjEGKMgO5NlO1CQU+W7dtE2ZtuUmp0VE9oHB1NE0djMDGxczJoXkZH9fYm0cak07avMTGh1IYNzd6LJifKJm+aBV3xpjUZnDZTsGGDlohs2486qplr9GIBxbJrdMXbmrxpMwBz5ujB0+SLNhPsKL1Qc3th1/NdY1i9rYmMZgiG9QU20Qv7SC8Iz3WN3XcfXm9rIqMZAFu4wOmnw9jYcLzEJtpeT6wXhOe6xnPPDa+3NdmMLOiCzahJuOwokNYsq4Ne2PXMa/TCXhWDJm1GiYws6BcZtWWUtL3ARfTzZTbRFcNsl7FmjVbNXn55p8S0dKm+dxs36jp1c4mHPoehzoHdhUIfXPttu4qL4QK2smFDs9ergph7MDGh1Pr1zbvsBwETE0qtXavU2NjO/jeZSzzmOTTp2u/7wO9i6TUZ9So2ZmJCv8C9uFaVtoW2qxcxPk2gDcK09X1iYioR1Xmuse/iUCbknymwudebdhW7XPh53qAuRuOG3oNBmbvXhPfSfI6uvm/c2Fwu8b5OR6rLZsNYaEkycn3Rm5SMQqWGNqJx6yD0HqxfP10C6IqamaOJ52l7jq6+j40Nh2TU94HfxdIGGZU95CamgMS+SF2zu4Tcg16ptHVQlzBdfRwfd/fdtBnVUV9j3sWkpg0gysTfJmJjYkTsrgZBFlUQGwZh0mfdoEjXc3z+eXffly6FRx/VsWNjY/D971cP12g7TsuJumw2jIUeSUYjI9orYn7Vq0osoVJDF6WLKlJdl9RME3Uk3bJ70aW+JzVtAMlIqakvqMhOsbooUtf1FIUMgi7aXbrYJhdsHwvXtqqkMSiZGxIZDSgZKeV3r/tsArHX8A2CYZCMise5pMheudbbCjXokgTkQiKjASYjpdxSwKpV9u1tSAe9+PLGkkFsm3wk0AZBuFTtrhF7L5HIqINkFDPwqnhL2kAbX978PqxcWY0MQtvkk6TakvxcH5FBUS/LUEWSTGTUMTKq8hV2SQGDYiuwoXgfzFKHDGyDxGdjasv+NMySUfHZiSi1bFlYH4aGjICz0As3bs3KJuC4kmMWAQ8CLwD3AQuN/Vdk+58HfgLcARwW2a5gMqrzFTalgHzO0fXXT513NAiw3YcmyCA2UHR8XN/DtgjC9rGo6znrd6yX69mFfFiHiYxOABYCc4EDgCuB7cBBjvrzgZeAjwMHAn+c1T+4UOcU4NeAtwAHATehV6t9Q0S7gsmoqa+wGbQGgyUVlakwbUQDmySwZMnUr3v+d9MSpk2VrKLydmWOne/ZlT23oSEja4PgaWCpY99a4DZj213ADZ7zzcpu1jERbeiJZFQ8hy2cX6TaFzP2a2vW9x0/Pq7U1VfrX/McPhWtykALIfqcBGz2tpGR7kqYXfJolkm1vg/rUJIRMAosBrYB73DUeRQ439j2SeAeR/1dgGXAM8DrPdfeNSOgvOwdazNyiemuwVuE78s0NhbUhCltifnamvWL0oV5/JIlU9u2ZIn/PqxYEW6MtpFfzIAdlDilXBW/9NJutfemm9wfkhkjGQGHAM9l6tczpg3IqLsd+Iix7WzgCWPb8dk5dwCPAfNK2nBFdkOnlFhvmjnwygZv8dgmyKhKFHOZnSc/fnzcvt8mIcUQUJnXLdQe0yVJwwWbKt6V9rrehZUr/ccNGxntAuwPHAosB57ySEahZLR7ds73A2uA7wFv9LShlmRkQ+jgzbFixfS6IyNxk1xDpYP82LVr/URUPP7qq+37rrnG3yYbYr1uRVXMd/6mPZFNGpddqnix30XjfK+N2lUly6Eio2kN0t6vGx37otS0Qp2HgYsi2lA76DFk8JooSgm+wVQnFYnpwvUNkBjJKCZ1SRWvW5XUKLlKZJv/ZzuuSABNG5d9qvg117R33VDYnkuIzXLYyWg9cItj31pgnbFto8+AndX5DnBFRBtqk1GsZJSj7lQOn3TgeuFMj5TreJfa2YRtp0wyilXBYtKwmgSwYkXzKp9LMjInwLahasZIrGb7XKaFHENDRpladgSwX2Y7Wp7ZeRZk+28FlhfqzwdeBC4A3p7Zel5x7Wfq2VWZerZvpvr9BTomyRou4GhXI9NBQm1GMYjxMJkvny85lxnv5CLE8XH9JS+SaoyIX8XrFqtChAx8X3tc7atrXC7LOdSGET5G0qryAR0mMloDbEZ70J7MVLQFhf13mlISOujxoeyY+ykYvIHdgH9AG623AT8AvkiJAdvSrsbmptkGbx3YBlpoCECbX96Y85rqaJnXLeT8xa+/T/oyB7arbhPhGi4P4diYPeSg6ecTe74qpoWhIaOulibJKEdTRkkbGfkM3Sbamm4Set6qwZ2+89vUrFDJyDUfcOXKavdpYkJPpahq92ny+cRKWjNaMupqaZqMTKPxihXVz9WEKJ9/nV2G3YkJPSXl0kvtrntfyo4yCadOrmbb+V1f/5Urp9uMzPlWNq+e6dUK8eIVz2frX9nHwhZ02sQE5iqSVqxpIZHRAJGRL34jxtvjO18Tyd6L+8y25i9kXU9PjPpU95z5YB4bU+rMM6e322UrMsm3KS/hsmX249r2nlWRtGJMC4mMBoiMXIPFdKvHLLoX+4KNj2sp5/rr/WlKfMGX69a1N+3FRgIx5/SRiou8XTFWpiOgKS+hyHRJs1eBmk1JWjYkMmq5NC0ZlcXyxLyIuUg/Ph72gplit08y8Q2os88uH7wh8KkyTaxmURz8y5a5SSdkMcsYldgWtOoquaQ5KFNYfEhk1HJp2mYU86L63NWxhlGXQdJFgj7J6Nxz3cfFYmJCqdWrp5NSHanAJvHlpOS6TpmEGSq5hARymmV83H3cJZf0P61IKBIZtVza8KaZc7BiJCOXNFE2eF2u2mIbbJ4ps+5JJ1Wbt+RD01KBT6ozAzvNQFCfhFlngQMfQeXuct+HKlRabMpTWwWJjFoudcmoaKMxPT/5i18WAFc8pmp6B5dklNtUXB61XHK5/PKdnqQmicPVr6ZjaorFDOyMPXeVOKjxcX0fbe25+Wb9sbj+enebQ+5Jv3MiJTLqMBmZNhqfYboYAOdyH/u++CMj5Z44V3tiXmLfYIvJ+23WbTrmyTUBNyZ0oKqE4euL+Qze+lY/AYWSfqhHsE0kMuooGfkkEVe0sFJ+YnB98V3rrrnadfnl+iud24aqzPVyZVUsu35Z/2xSR1ViyG1rsSTXhIThk6Byd/nNN8cRke+5+Dy1vZKQEhl1lIx8Npr861ZlUmaRCEZGlDrjDL97vix2yTdHzYdiEGAomVUlvjaJoYk2VoXvHQGlfuM3wonUp5rGROXXQSKjjpLRwoX+F8MlVoeI5sXB5SKTZcvCYpdcL7E58F1ZKmNsSFUmubZJDDaJq5cu9jIPZ65mhRKpL2FbbJbQKkhk1HKpQkaul6xICE1MyswlHxupxXjcyuwrvmkBMYQRa28KzUpQRYVzSVx11qur0pay2K9YEnQZwfPJuG162hIZtVyqkJFL/L788p116k7KNOe4Fd3zy5bFv9yugECXB6goIcUYn2PsTWUZL6uqcC5SLIZcFLfXsTO5CKC4fXxcqRNOmN7XKlKgLbBWpPpCmjFIZNRyaVIysmVCtL30VdzHIyM7v36hxvOyc4JSxx9vP9fZZ7tDFYrbXAOxzN7kilbPpcs6KpyLrF3eqDKJIoTcigTgIi6zftVJ1CG2SJ/3taoElcio5VLVZhSTCTHWBVumvrj2n3yy/wU76ST74HdJWb4vbIjU4uuHLyd36HwyG1wk57LXhUS6x6jcZWpgUxJMiF2xSFY5+dRxGCQyarnUiTPKXbjr1pUn+4pN++GzvbimQ/heMNcgBaWOPdb9Msemhi1+dX31yhYICJlPZoNLKvqd33Hfryr9dd3/Vavcz7+JGK68Pfl7UPY+FItvukwIEhm1XOpGYLeVU7nM9lLMXx1CIL6ve676uSbIjo1NHSw+D1/RzrVihdve5CPH0PlkJnxSUejqKK4Ph9kW13P2ZTzwxQrl7Szro80EUPY+VOmvDYmMWi51yMhnS2gi2nh8XH9pXS94vj/kBXPZjIpzzmx1isbzfLC4JAXX+V02MtsM/Nj5ZEX4SDK0zS7JyJY9wVTV58+3XyO/vi9WyHf9smeYX2fdOqUuuyxMSgq5nolERi2XOmRUluyrTl4Z05vmu07ogCpGK7sMqGbQZegseJ/R2HcPisbu4rJDRTUkRn2xkUHRO2d+JGzbiqqmz4tWNuhtRB6S1cElqYTkUTKv63omNuIvQyKjlksbklHdOI/QF903yIowie2MM/xR2zlBuFSbVat2kkaRQGIHlw0uNaRqMrqLL54azOnyChYnNRfvlSsHeRkxuKSQMnWxqmRkIxxbOAHo51/lQ5nIqOXShM2oyQmgSoW96KbE4ZLEykR7UxKIPbZYt8rgMvNBu64XQ/L5vVi0aOo5Fi60Xzekv8XiUvtCytiY31ZW9v6EJtCrSng+JDJquTSRz6iuSmY7X8iL7lrKutiOENHeZ0C12XVsL7XPXmODTf0pa2veX1+cTL5v3Tr7OQ45xO3aDpV2RHQktGkbLDMiF1XcsgUE6rwXIaXK9JdERi2XNpKrNYFQ202xflXbRtlXc2Ki3FBeZq8xzxcao2MbyC4ycU17Ce1v7L3KPYbmopi2sATQ4QW5epzXs62p5kIV1dBVjjkmPv1IIqOWS1fJSKnptgxfBkOf7con3YR+NUPsY6Eqq8/w7/OwlcU3VZUaiv2NvVcxYRRFEquizq9bF55jPbTErHw8NGQEnAXcC2zNyibguJJjFgEPopesvo+pK8q+Gvh0tv159IqytwJvimxXK4s4xi5LVDw2dq2ymGWwi6qFq5hZK3OYq8O6givLVNYyYrN52Mr6WVVqcEmCoffKRt5VXPi+565UM7YiVwmVkIaJjE4AFgJzgQOAK4HtwEGO+vOBl4CPAwcCf5zVPzjbPxv4F+Bk4G3A+4GvA9+MbFfjiziWpfZwvXhNTwz1GY5dqkSxXHyxvW11F6fMzxdr+Pf10zVf7wMfmPr//PnxgZQbNrilEhH7YA6RsFx5r8w2+eYihpCNy6uWF9+S1kUMDRlZGwRPA0sd+9YCtxnb7gJu8JxvXnaz9vHU2TUjoLzs3RQZTUyUp/bw2XZ6vZR1yIBZtKjdEIZYw7+rnz6pyVyosKrDwSWd+HJJubxnPvWyeG/Hx5X60Ifs17388rCPyvi4lvKSZGQnhFFgMbANeIejzqPA+ca2TwL3eM77a8AO380Crshu6JTS5iKO+aDwGW4vvdR+XEzSrCqDzGeczosrh07ZAgHnnadL03maXbFCTRJmTNhB2bVMabn4EfKRaJlqVswS4fqo+CZw5x+bUAwVGQGHAM9l6tczRRuQpe524CPGtrOBJxz1dwPuBv6qpA2VJKMynT6v4/sKVlnmJiejkOs33e7iF9i2vThQc3f61VfrdKq+QVEnsb8PTcV8VQk7KCNnl/fMRaKu0ATzfhbPkxO0bcnqJiZwDxsZ7QLsDxwKLAee8khGwWSUGbP/CfhW7I0KsRnF2HJ8X8EqruMyd7YPse12tePmm+3bx8fj3OmXXBLXnqr2szoxX1XCDoofnCqwkagrgd+JJ1aTNJuQHIeKjKY1CO4AbnTsC1LTMiL6R+Ae4BcqtMFLRlUeoi+GxIwfcr3cZZNSy17Iqu0+8sipxyxZ4s9sGUuuIe2pOnCqSo9FuObY5WEHTS/XXWx7kUTLEvhVQV3JcdjJaD1wi2PfWmCdsW1j0YBdIKL7gTdUbIOXjFzi7WWXxblmiyi6rm1EUyQxn2pnxhqFpPnIxXJXAv58X1HMdw2M884LIyFfaSKxv1L1koYV72GssdnlTWsCvtzkVVFHchwaMsrUsiOA/TLb0fLM2Lwg238rsLxQfz7wInAB8PbM8Fx07b8a+CIwAbwL2KtQdoloV7RkVHwR85nfVQdC2dfKd31fVLJPuqjyktsWiIwhHduy2U1JRi7pMXbAuSax5lNa6thdqkptNvtPvzBMZLQG2Iz2oD2ZqWgLCvvvNKUkdNDjQ9kx9zM16HG/7MbYylER7Sq1GZWtkV4WwBdiQHZ9rSYmpiYtM4svI6KN6OqI//kCkSGq2cKFSn3sY7rkeaZjFlyMUSlcJHH88eGBpz4VzOdRC1Efm5DauoChIaOulhAyqhLdm9sY6ryEpnHYFodUliu6SHQTE0qddZa9fmjgW9m9yMnH1Q/XpFAbaYeqFGXeQFsc0Pi4Ur/7u0oddZRSixfbCdaVjsWX/8jWtjZitPqBREYtlxAyClmMz3zZ6qzPpZTbPmFGd4emuS3zeoWqAWUeQZPUbJHLoZN8Y+BbvqlI3DY11VVWr3bfA1v+I1vbqy4s0EV0iowAAaTuebpUysjIJb4Xi22KQd24jrJ8yS6Csg2IMgKJNYz67kmR1HyDvmy2//XXx83tiw2bCCllxFgm9bjuU5KMapARsDSz2Wwr2G8+WrdBXSg+MioT/4tl3brpqSRcElOIIdM1SEMGkHluF7Gdc850iSjU0DoxMX1VkSKp+aTJkDxIRbLNbV0u71+OKilEyoqPOFxtz72hoapfSN9i0USYg4m+kxHwKXTU9HLgQ1lZDjwLfKpuo/pdfGQUYyuy2VxM+4JvdVUbzONtKplP6sgRareooi65vD2u2CTTflNFovFJcrmhvMlUGy5p1tX2kRG32mhO72nDfd+WwbwLZPQURiR0tv0jwI/qNqrfpUnJqOrqqj6YRtyyhO6uc1YJIXCl1whJj+KSjNat87cttIQEfY6NabtPKNnFhB8U2+4ipLJztRHY2KbBvAtk9Aww17L9AOCZuo3qd2nKZlT2JWpibpDvPLZrm6K6zzsV0r6Q9ChFmF/9hQvd6khOHqFkFOr9y9ttC3H46EeV+uAH9a9t0mlolLIvPsl3Lpf0aOtbmceuziKioSpdF8joz4BVlu2fAa6r26h+lxBv2sSEUldd5ScB83+bLaaJL5ZLWjOvGSuql7WvLDrZhVyNMyfPFtWR4mAo84pVlR5iIo/NcIiQODHXvfNdN1QyCk21OzISv4hozHvSFzICVhXKn6IzM94P3JSV+4AtwJ/VbVS/SwgZKeVO8eEjKHOZnCZnlZsetBUrtAfq0kv9q5qWndfVvrL0KD74Bl3IYDJLE3YVEzb10xXZ7kqMFxN/lKPMZlRGdLZ9rkVEi+3J+xvznvSLjDYElvV1G9XvEkJGvhntZQvmmS9ZzBfah1ytGRvzJ86KIQ1f+6pKRkrFTbQdGdHrepnq4MUXtzctwqZ+2kjRtiijeY9C44+K8E358HnsYhYRNQNPfaYH13vSdzVt2EsZGfmM2CKaZEKN3K7c0nUwMRF27SaMmL70KD64JCNX0ra22u+K8napvVXb1KQR2eexsxFmqImgyn1OZNRyKSMjn3oS46kxX6KiuFwnHsTXvnyQ1VEJTRQlspjlp23qSKxbP//SV7lfLkkl5P6VtclEHWdFnvXzyit3qowuj52pkrk+EqEhKmXvSSKjlktVyWh0tHyp4rLiWkAxZsD5JCMzELNJNBWTFOrWzwdelfiZMrtLqGQbmo/JFn4RIhnZItZzj6XrXTvzTPf98gXgmuQVsn5bIqOWS6jNyKaeuKKkx8amJykL/TJVGXA2m1EVI29oJHDTsSwTE1rK9N2XWC9REWWSSkj4Bih17rnlDgjXwF+50t/Gsoj1skyTvv7lfczbXrR9xUjNiYxaLqHetKJ6UvziuGJJqkQW5wTk+8KZEpMpwp9wgk4TGzu9wPwqn3SSWzprKmaqCNeX/5prtITnyjYQapQvI7KcEC+7zJ00zraCbFP3xmXkLx4fmiXURdRmyEKs1JzIqOUSSkYmQtJixEYWu14wMx2JiF4q+bzzmnGD+77KoRNvi2pPrE3H5a0cHdWrV8QMON81fFJN6Ly2kNS9VSS4EMlo/fqdC1v6TARN2giLSGTUcqlCRjEvXC5RXX/9Tje8TVwWUeqUU9yTa+tMAC2TkMq+yraUHza1tYodyWezOffcZgecL2wh1pDuQ9V4MpvNaGRk6jQVn4kAtA2pDRuhUomMWi9VyChGFHcFzuWDYuVK92DM61ddujkvZVMnyvI1wc64Flfebt/2KvcS3KrZkiVTgxNtAYsxcLVh6dJwo7WJqvFkeTbNq67S9/ySS6a3y5bJM/fStolERi2XNiWjsnplnpw8HUjolzskx1B+XVOV8uUeMiU4Wx3XQpBlEpLPW+laO6xI7macV3G+XKjK6Lu/VY29TcD3ftgyebaNREYtlzo2ozJRvEyCCgkNsK3HbiurV+sXsmx6gU+Vyl3vl1wy9Ytb5mkqUyVDZr67gimrLARgeiVzw3Po87SRX4jru2n4VobpdVuUSmTUeqlKRkqVf5V8klGIOzmXDnLPWO69M+uZhOOaXhBr69qwQU/NcEkMJhH7CLPMzjIxoVWTU0+dnmYk709MxHYV9/rEhFvCC/Xa1VEXi+dxqcQh/WgLiYxaLnXIKASuCZQhrtn586duK85vW71a2xZi3PexbucyNdImcVTN2xSSZCxUXXXVCZEoqnrDYtOr5NfyhWrkamiTdqE6Ef+JjFoubZORUnqArlq1kzhcpJCTy4YNbntJLxNvhRjOXRkjfSrs+LgOSzjvPP13TJIxW/ZMWwaDMluLD7HeMJ/dy3VvXY4Nlze1CbtQSDiKD4mMWi69kIxCX7riixGTeCu2PcUvre+FDJVEVq+eHmTpUmFthvIDDojrq3luW0Bq1SkZOcwPiA+x6VVcz7/NlURcz9KcK+lDIqOWS5tkVGYzcn19Jybc9pEmUmjktqeQmKBLLpk6f67MzuULsgwJIfD1NTZxfdGIHeMJq5KYLkYy8qUFsa2N14SxOkTKLevr0JARcBZwb5aobSuwCTiu5JhFwIPAC1lCt4XG/v8KfAX4cXaT3l2hXa2RUZmNxiY9+IzAMfPNfLaB0Ahq2xQRc45TDGGWBVf6+lo1cX2s67spm5FvYLuuYYspAj29py5CPwQzIoUIcAKwEJiLzp99JbAdOMhRfz7hulaXAAAWGElEQVTwEvBx4EDgj7P6BxfqnApcBny0K2RU/HrHvti2+iLxhuqyL7uLJH3LaBeJJh/gruyXLvUqZEDYlk9qI3G9C3Xm3dnURRtcmTp9Uuf8+e30K6avQ0NG1gbB08BSx761wG3GtruAGyx19wslI2DXjIDysndTZGT7escYQ5uYgBpCgLY6IyNh9qEi0fhWAHFJZb7gShdR++xn+SRa26ojVVBVMqpz/pGRsJizOn0Mtf/NCMloSkNgFFiMXhDyHY46jwLnG9s+CdxjqRtDRldkdaeUumTk+3qHqgpNDIRQQjNJsmoifJNcQlZKGR9X6mMfU2rBgjCbjuvevu51069dhhDXdlO5ym3w2YvKyOKcc+pd2+a8cOXLtmGoyAg4BL0g5EvoJZAWeupux1ivDTgbeMJSt++SUVPer7oDIWaqytq1O1WKkC+ny06TByVWWQgglKh9ElWo9BBjmG5rmkWZU8P3DJqQ/mz9Cu3rsJHRLsD+wKHoVWmf8khGrZCR5dhGbEZ17RqmranOQCibQOkalDYi9CWLN9FGnqMcuYfxIx/xk5FLemhb/YpBmSfVJqXWtRk1gaEio2kNgjuAGx37WlHTLMc2ZsCu6vFpY4ljV8bIskFZhwjbGvChuYZ80kObRFkFIVOJrrpKvwtN2cPqYtjJaD1wi2PfWmCdsW1jXQO25djGvWkxS+r0eonjtgdl0/aWUMNrmfTQJckoFnWmcPjOEXveoSGjTC07IiOOQ7L/dwALsv23AssL9ecDLwIXAG/PDM+ma/91wLvRIQMK+O3s/70i2tX6dBAf2oi09hlJezEom7S3hLik99wzTHqoSpRNkEFVhNq5XGSzdq2e7Gyeo0oivGEiozXAZrQH7clMRVtQ2H+nKSWhgx4fyo65n+lBj6dlN8csV0S0qxYZ1X1ReyUZmS9iW96iplEmGV11Vfz5YoiyyqBtCqEfDlsbfVkhbGEcIR+koSGjrpY6ZNTUi9qGzci31lbuPetVUq66KOtLW+i3aheiUrsybIYuvxSjqicy6igZ+YLXqrysPltTVemriYmX/VRRzHa4VmJp4ty2PjZpX6tyH8vIMHSJpSQZDUipSkY+W4Yv2O/SS+OWuTalr9AZ1krV/7L3U0VxoWmJrjigzRxETUlGtmcYmoTNpVLHGPbN9ldV1RMZdZSMyl4G86W1pU8NmQ3uSxQWcnxVaaLfKkovMDExXbIQmdrHNoJQq7wHJgFXWaTh5JOrBTvmSGTUUTJSqnxdtFycdxmpzfQQpihf9sKFJu+qkkira3E5bcClxo6NTa0XM2hjn2FVko+VjJpIRdIkGY2Q0CiWLoXNm2FsDEaMuzs6Cvvvr//+2tfsx+/YAY88ov9eswb23ReOPlr/rlkDc+dOP28RL7+88/giJifh9NP1+UG/jtdcE9U167WLfRo0TE7q5zQ2pv+OwZw5cNRR+teHKs8Q3M+xrE2f+5x+JqCvIWKvOzqq65a1v6eoy2bDWGgozsgnzpdJRmXzlVzSV5XkXTFG1EEKAfDBlZ/apqZVlSCqPsO66m9RajOfV9lS3LFIatqAkJFSfnHeZzMKTcIWOsPaNjCK63/FGKMHKQTABhvhuIiiTsL7mGfoSsLWVKR1W88rkdEAkVEZ8tVC8zXOcsQYi0NfNnOQ2VZGzddvH1SiCUFIfmrXPL4iyogi9hmaSdi66Lk0kchoiMjIhzZUopy4XIbarr/8TaBMMgohkVCiqDPdZBA8l4mMZggZKdWeiB3ieeniy98UfPmpQ9Sr2NTBsc9wUDyXTZLRq3pmKU8oxeQkPPyw9rbkXo45c9rxeOSelzPO0J6bkZGdnrYcuUenUx6XhrB0KRx7LGzapP8//PCd/cy9XcX7UfQaPvxw+L2yPdMQlLVhKFGXzYax0AfJqF/2gfyrXXXV12FFWbKzqpNVm2pDXTQ1pSepaUNGRl2xD1RZNXWYjd0+9arsXjX1TNtQ05v88CUyGjIy6pJ9IMYzNxOM3T747lWXnmkRTX/4ks1oyNAl+0CIjcqM5t6xQ///znfCvHntt7Er8N2rLj3TImLsXb1Gmg7SAZhh/KOjcOON/X85XLC90Dt2wPvfr6c7JNif6fLl+t7FTj1pEl2e0iNKqyUJBYjILGDLli1bmDVrVs+uOzmpv1D7799dIgLdzn33nU5IoF/szZu73f5eIn+m3/wmXHihvmcjI5qoli7tT5vWrNnpRc0/fFXbsnXrVmbPng0wWym1tU67EhlZ0C8yGiSsWTNVVStiwwY9iTRBw0be/Sbtpj58TZJRUtMSKmHpUrjrru6K/F2Cz07jw+SkJvY21LrQrAO9RCKjhMqYN2+wbF39QhU7jS31SCzaJLNWUNcdN4yFDk0HGQQM+iz+MjQRTxUTw9WE+71XoRdNuvaTzciCZDNKyFG0jdU1PIfaaTZs0BKRbXuILa6XNqpkM0pI6AFs8VRnnFFd7Qm109R1v1e1UfUbfSUjETlLRO4Vka1Z2SQix5Ucs0hEHhSRF0TkPhFZaOwXEfmUiDwuIj8TkTtEZG67PUkYRvRrUNeNO+tyLJEP/ZaMJoE/Ag4F3gusB74oIgfZKovIfOBv0CvRvgf4AvAFETm4UO0TwB8AZwKHAc8Dt4vIbm11YlAwcAbNllF2P/o5qPNc6hs26N8Y1XDQgmhfQV2jU9MFeBpY6ti3FrjN2HYXcEP2twCPA8sK+2cDLwCLI9owdAbsNJdsKtpOjtYF9MKxMJQGbBEZBRYBnwfeo5R6wFLnUWCVUurawrZPAicqpd4lIm8BvpMd/+1CnX8Fvq2UOs9x7V2BXQub9gAmh8WA3cWgu34i9n4MSmR8P9CkAbvvE2VF5BBgE7Ab8BzwYRsRZdgLeMLY9kS2ncKvr44NFwGXh7Z50NDlyZH9QOz9aCvBXcJU9NtmBPAQ8G60fWc18HkReUeP27Acrc7lZahevUE1aLaFdD+6ib6TkVJqu1LqEaXU3Uqpi4B7AKs6BfwQ2NPYtme2ncKvr46tDduUUlvzAjwb1YmOY2ANmi0h3Y9uou9kZMEIU+03RWwCjjG2Lci2A3wPTTqv1MkCGA8r1BkoNOUBq+OdGQaY93Gm349Ooq4FvE5Bq0dHAPsBh2T/7wAWZPtvBZYX6s8HXgQuAN4OXAFsBw4u1LkQ+AnwoeycXwC+C+wW0a5OeNOSB6wZpPvYHoYm7Sw6XmgzsA14ErgjJ6Js/53ALcYxi9B2pm3A/cBCY78An0JLSC9k5zwgsl19J6Ou5MUedLR9H4c9D3gZhibtrFLKKxwrpY6ybPs74O88xyjgsqwMLJIHrBm0eR+bnLeW0E2bUQLJ49MU2rqPTc9bS0hk1Fkkj08zaOs+Dupk1C6jMxHYXUKXUoiERP9WXbW0bXSpXU1HUaeodo2UQmQGoSztRBMZAdtA19rVdJrVLkmuwzIBOklGFnRJMvKhq1/nrrarDfR73lq/jehJMkoAumu36Gq72kA/E9sPmxE9kdEAo6set662a9gwbKSfyGiA0SW7xSC0q5fohR1n2Eg/2YwsGBSbUY5+2y1c6Gq72kYv7ThNrg5bBWlF2ZYxaGSU0B30w3jfT9IfquRqCQnDhH5M4xmW5G/JZpSQ0CCGzY7TSyQySkhoEMl4Xx3JZmRBshkl1MVMMd4nm1FCQscxLHacXiKpaQkDh2GZi5UwFYmMEgYKXZuAm9Acks3IgmQz6iZm0gTcQUGaKJswIzFsc7FcmKlqaCKjhIHBTIjhmclqaCKjhIHBIMTw1JFqhi0lSCwSGSUMFLq8+GJRqtlnH1i5Mu74maKGupDijBIGDl2M4TGlGqXgE58AEVi2LOwcuRpqGuiHSQ31IUlGCQkNwCbVAFx4YbiaNQhqaJvoKxmJyEUi8g0ReVZEnhSRL4jI20qOebWIXCYi3xGRF0TkHhH5daPOHiJyrYh8X0R+JiIbRWReu71JmMmYO1dLQSZ27IhTs7qshraNfktGRwLXAe8HFgCvBr4iIrt7jvkT4AzgXOAdwA3AP4rIewp1bsrOdypwCPAV4A4R2bvxHiQkoKWXT396+vYqalY/82r3E50KehSRNwBPAkcqpb7qqPMD4Eql1HWFbX8P/Ewp9d9E5OeAZ4HfVEp9qVDnbuDLSqn/aTnnrsCuhU17AJMp6DEhFp/5jFbNduzoT+bFXmOYJ8rOzn6f9tTZFXjB2PYz4Feyv18FjJbUMXERcHl4MxMS7Fi2DBYvnhkz9ptGZyQjERkB/gl4rVLKRRqIyF8D7wJOBL4DHAN8ERhVSu2a1dkIbAdOAZ4APgJ8HnhEKTXNJpUko4SEahjW6SDXAQcDi0vqnQc8DDyIJpzPAjcDRV/GqYAAjwHbgD8A/sao8wqUUtuUUlvzglbzEhISeohOkJGIfBY4HvigUsrrCFVKPaWUOhHYHdgXeDvwHPDdQp3vKKWOBF4DvFkp9T60cfy7llMmJCR0AP127UtGRB8GjlZKfS/0WKXUC0qpx9A2ot9Cq2pmneeVUo+LyH8BjrXVSUhI6Ab6bcC+Dm3X+U3gWRHZK9u+RSn1MwARuRV4TCl1Ufb/YcDewLez3yvQpLoiP6mIHItW0x4C9gdWotW6m2Mat3VrLRU4IWHo0egYUUr1rQDKUU4r1LkTuKXw/5HAA2hv2Y+AW4E3Gec9GW3c3gY8jrYrzY5o196etqWSSirTy951+aAz3rQuQUQEeBN+Q/YewCQwp6TesGEm9nsm9hnC+70H8ANVk0z6raZ1EtlNfcxXR3bG/j9b16U5SJiJ/Z6JfYaofjdyTzrhTUtISEhIZJSQkNAJJDKqjm3AJ7PfmYSZ2O+Z2Gfocb+TATshIaETSJJRQkJCJ5DIKCEhoRNIZJSQkNAJJDJKSEjoBBIZVYSInCMim7M83F8Xkff1u02hEJEjRGSdiPxARJSInGjsFxH5lIg8nuUQv0NE5hp1XicifyUiW0XkGRFZIyKvMeq8U0S+lt2jCRH5RC/6Z0NIvnUR2U1ErhORH4vIcyLy9yKyp1FnHxH5koj8NDvPShF5lVHnKBH5lohsE5FHROS0HnRxGkTkLBG5N3tGW0Vkk4gcV9jfrf72c27aoBbgt9Huzv+BzsP9OeAnwBv73bbA9h+HziX+YfS8ohON/RcCz6AnML8Tne3gu8BuhTpfRk9WPgydQfNh4K8L+2cBPwT+D3AQOk/VT4HT+9TnfwZOy9ryLuBLwPeB3Qt1VgOPAkcDhwKbgH8r7B8F7gP+BXh3dh+fAq4q1Pkl4HngauBA4PeBl4Bj+9DnE4CFwFzgAOBKdA6wg7rY374PjEEswNeBzxb+H0FPH/mjfretQl+mkBE628HjwLLCttnoicmLs/8PzI57b6HOr6OT170p+/8sdPrgXQp1/hfwYL/7nLXlDVkfjij0cTtwUqHO27M678/+Pw54GdizUOdMYEveT+DTwP3Gtf4W+Od+9zlry9PA0i72N6lpkRCRXdBfkTvybUqpHdn/h/erXQ3il4C9mNq/LWgCzvt3OPCMUuqbhePuQJPRYYU6X1VKbS/UuR14W5Zfqt8w860fik7AV+z3g2jJodjv+5RSTxTOcztaCjyoUOcOpuJ2+vxuiMioiCxGJyXcRAf7m8goHq9Hi69PGNufQA/iQUfeB1//9kKv4vIKlFIvoQd2sY7tHMVr9AVZvvVr0SrJ/dnmvYDtSqlnjOpmv8v65KozK1u5pqcQkUNE5Dm0WeEG4MNKqQfoYH/TrP2EmYg837pz4YchwkNoe89s4CTg8yJyZH+bZEeSjOLxIzI92ti+J9pgO+jI++Dr3w+BNxZ3Zh6W1xl1bOcoXqPn8ORb/yGwi4i81jjE7HdZn1x1tqose2kvoZTarpR6RCl1t9LZUu9BL2rRuf4mMopEZgO5G71EEvCK2H8MWhcfdHwP/YIV+zcLbQvK+7cJeK2IHFo47mj0+/T1Qp0jROTVhToLgIeUUj9pqe1OBORbvxt4kan9fhuwD1P7fYiIFIl4ATqfzwOFOscwFQvozrsxgl6Wq3v97bd1fxAL2rX/ArAE7Vm6Ee3a37PfbQts/2vQovu70d6Tj2V/75PtvzDrz4fQy4N/Abtr/1vA+4APAP/JVNf+bDSp3Yo2dv422gXcL9f+9ehwhSPRdo68/Fyhzmq0u/+DaAPvRmBjYX/u6r4dHR5wLNp2ZnN1r0B7p86mf6795cARwH7Zc1yOdjIs6GJ/+z4wBrWg4ym+jzYMfh04rN9timj7UdjzGN+S7RfgUxmZvID2lhxgnON1wF+j05FuAf4CeI1R553A17JzTAIX9rHPrtzNpxXq7Ia2Jz2dDbB/APYyzrMv8H/RMVNPAZ8BXmW5v/+evRvfKV6jx31eA2zO2vFk9hwXdLW/KYVIQkJCJ5BsRgkJCZ1AIqOEhIROIJFRQkJCJ5DIKCEhoRNIZJSQkNAJJDJKSEjoBBIZJSQkdAKJjBISEjqBREYJQ4UsFfD5/W5HQjwSGSVUhojcIiJfyP6+U0Su7eG1TxMRMxcPwDx0GuCEAUPKZ5TQKYjILmpqdsgoKKWearI9Cb1DkowSakNEbkHPhj8vW21Eich+2b6DReTL2eoTT4jIX4rI6wvH3ikinxWRa0XkR+gZ4ojIH4rIfSLyfLayyPX56iMichRwMzC7cL0rsn1T1LRsdYsvZtffKiJjxRUwROQKEfm2iJyaHbtFRP5WRPZo964lmEhklNAEzkPnr/lz4BezMpEl7lqPntH9XnTS/j2BMeP4Jejk8B9AJ3wHneriD9DpR5ag8yWtyPZtBM5H59XJr/cZs1FZnqkvojMMHInOs/MWYK1R9a3Aieika8dndf8o6g4k1EZS0xJqQym1RUS2Az9VSr2SxVFEfh/4d6XUxYVtv4smqgOUUv+ZbX5YKfUJ45xF+9NmEfmf6BzOZyultovIFl1N+bJGHoPO4/NLSqmJ7Pr/HfgPEZmnlPpGVm8Enfbi2azOX2bHXhJ7LxKqI5FRQpt4F/DBLCG8ibeiE7KBzjo4BSLya8BF6IRds9Dv6m4i8vNKqZ8GXv9AYCInIgCl1AOZ4ftAICejzTkRZXgcI61uQvtIZJTQJl4DrENnjjTxeOHv54s7MnvTbehMhJegk3/9CjpZ2C7oRF9N4kXjf0UyYfQciYwSmsJ2dJrSIr4F/BZa8ngp4lyHosngAqXXpENETg64non/B7xZRN5cUNPeAbyWnTmcEzqCxP4JTWEzcJiI7Ccir8+Mx9ehjcd/IyLzROStInKsiNwsIj4ieQS9wOC5IvIWETmVnYbt4vVeIyLHZNf7ect57kDncP4rEfllEXkfOif3v6qpC1AmdACJjBKawmfQSzg9gM6VvI9S6gdoD9ko8BU0MVyLToy/w3UipdQ9wB+i1bv7gd9B24+KdTaiDdprs+t9wjgNSudU/k304gJfRZPTd9GLAyR0DCkHdkJCQieQJKOEhIROIJFRQkJCJ5DIKCEhoRNIZJSQkNAJJDJKSEjoBBIZJSQkdAKJjBISEjqBREYJCQmdQCKjhISETiCRUUJCQieQyCghIaET+P/U5B49CRuLwAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# generate mcmc plots\n", "mcpl = mcstat.mcmcplot # initialize plotting methods\n", "mcpl.plot_density_panel(chain[burnin:, :], names,\n", " figsizeinches=(3, 3))\n", "mcpl.plot_chain_panel(chain[burnin:, :], names,\n", " figsizeinches=(3, 3));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Plot Prediction/Credible Intervals" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Generating credible/prediction intervals:\n", "\n", "\n", "Interval generation complete\n", "\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAGoCAYAAAAJjpFOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XlYVXX+wPH3uZd9X0WUTVEEAQUFETW3SsW0cpqWmWaqqWydmWqW9sZS25tfU9PMNFlNNdXY4pIiqGUuWaFCoILgAiKCKPt6We5yfn+YJMpygXsR5fN6Hp8nLuec7zk8PXz4fs/n+/koqqoihBBCXGiaC30DQgghBEhAEkIIMUBIQBJCCDEgSEASQggxIEhAEkIIMSBIQBJCCDEgSEASQggxIEhAEkIIMSBIQBJCCDEg2FyogX18fNSQkJALNbwQQoh+kpGRUaGqqm93x12wgBQSEkJ6evqFGl4IIUQ/URTlmDnHyZKdEEKIAUECkhBCiAFBApIQQogB4YK9Q+qIXq+nuLiY5ubmC30rA4qDgwMBAQHY2tpe6FsRQgirGVABqbi4GFdXV0JCQlAU5ULfzoCgqiqVlZUUFxczYsSIC307QghhNQNqya65uRlvb28JRmdRFAVvb2+ZNQohLnkDKiABEow6ID8TIcRgMOACkhBCiMFJApIVhYSEUFFR0edjhBBiMJCAJIQQYkCQgHSOwsJCwsPDue222wgLC+Pmm2/mq6++YurUqYwePZrdu3dTVVXFtddey7hx45g8eTL79u0DoLKykjlz5hAZGcmdd96Jqqpt1/3www+ZNGkSMTEx3H333RiNxgv1iEIIMSAN7ICkKJ3/e+utn457662uj+2hI0eO8Mc//pG8vDzy8vL4+OOP2blzJ6+88grPPfccS5YsITY2ln379vHcc89xyy23APDMM88wbdo0cnJyWLRoEUVFRQDk5ubyySef8O2335KVlYVWq+Wjjz6yyI9ICCEuFQNqH9JAMWLECKKjowGIjIzk8ssvR1EUoqOjKSws5NixY6xatQqA2bNnU1lZSV1dHTt27GD16tUAXHXVVXh6egKwZcsWMjIyiI+PB6CpqYkhQ4ZcgCcTQoiBa2AHpLOWvLp0112n/1mIvb19239rNJq2rzUaDQaDoccVE1RV5dZbb+X555+32D0KIcSlZmAv2Q1Ql112WduS27Zt2/Dx8cHNzY3p06fz8ccfA5Camkp1dTUAl19+OZ9//jllZWUAVFVVceyYWdXYhRBi0BjYM6QB6umnn+b2229n3LhxODk58f777wOwZMkSfvGLXxAZGcmUKVMICgoCYOzYsSxfvpw5c+ZgMpmwtbXlH//4B8HBwRfyMYQQYkBRVHOXxSwsLi5OPbdBX25uLhERERfkfgY6+dkIIfqLqqp8tvoLgsNiSIgO6fP1FEXJUFU1rrvjZMlOCCFEG4PBQHltM/l1Lnz80X8xGAz9NrYs2QkhhABOz4z++vq/UfxicPIYyryF16PVavttfAlIQgghqK+vp7jKgHZYArb2DgA0Nzf16z3Ikp0QQgje/XAVn32Z0RaMAHZ/t61fq8pIQBJCiEGsurqafUdOgt9E3L2Htn1uMhrx8PLh+eefJzk5uV8CkwQkIYQYxFK+/p6VqbvblVkzGY188MpDfPTvV2hsbGTJkiUsWrTI6kFJAlI3nn76aV555ZVOv7927VoOHDjQj3ckhBB919TUxN7cQqptgvELGNnue7mZOzE2VZORkc4LL7xAWloaxcXFpKamWvWeLuqAZDQaSU5OZtmyZf02pTyXBCQhxMVof14B763ZhtF0/l7U4oJckpLmtpVJs7W1Ze7cuWRlZVn1ni7agGQ0Glm0aBFLlixBp9NZdEr57LPPEhYWxrRp0zh48CAAK1asID4+nvHjx3Pdddeh0+n47rvvWLduHX/+85+JiYkhPz+/w+OEEGKgMJlM5Bw8wv5TNoREdLxXNWBkBKmpm9Dr9QDo9Xo2bdpETEyMVe/tog1IqamplJSUkJaWxvPPP2+xKWVGRgYrV64kKyuLlJQU9uzZA8DPfvYz9uzZw969e4mIiOCdd95hypQpXH311bz88stkZWURGhra4XFCCDFQlFXW8u+PU6nXtXZ6TETsNLSOnsTFxfPoo4+SkJBAQEAASUlJVr23i3YfUmZmJnPmzOlwSrlgwYJeX/ebb75h0aJFODk5AXD11VcDkJ2dzZNPPklNTQ0NDQ3MnTu3w/PNPU4IIfpTY2Mjm77eSaU2iNCYWV0eq9FqueVPr3IyZwvOmkaWLl1KUlKS1TfJXrQzpNjYWDZv3txvU8rbbruNN954g/3797NkyRKam5v7dJwQQvQHVVUxGAwUnNSx63Adjc3mlQLSaLXY2dnx2GOPsWDBgn6p2HDRBqSkpCSGDx9OQkICjz32mMWmlNOnT2ft2rU0NTVRX1/P+vXrgdO7mP39/dHr9e26vbq6ulJfX9/2dWfHCSFEf/r08zXkHi1jw7YfePy5f7Bl70mGjYhA6UEX7bkLf46NTf8tpF20AUmr1bJmzRqWLl2Ks7MzS5cuZc2aNX2O4hMmTODGG29k/PjxJCUltXV5XbZsGQkJCUydOpXw8PC242+66SZefvllYmNjyc/P7/Q4IYToD6qqUlBaR9ZJG1L2FHOk2o7h0Zf3KBCd8e32L/s1e1naT1wk5GcjhDDH6rXryS6zwd03sM/XcjOe5NafzUaj6dvcxdz2ExdtUoMQQoj2VFWlySkEF0/LtIwYHR7V52DUExftkp0QQoif1NTU8Le3Pqa8AbQ2tha55kfvvNGWONYfJCAJIcQloMmgpVbx6dW7os7ceMvdktQghBDCfDk5B/hixwE8fYdb9Lonio/Rn3kGEpCEEOIitzf/FJV1LRa/bv6hXEwmk8Wv2xlJahBCiItY7uGjVOOHi3vfAofJaCQ3cyfFBbkEjIwgInYas+YulCW7M4rLG/nbik946+MNpOWW8fuHn2FLegFrvkrnqedeJy23jJf/+V/e+/xL0nLLuOfBx9mx9zifpHzHMy//i7TcMp5/7V0+/GIbabll/OHPj3c75muvvUZUVBSRkZH87W9/a/v86aefZvjw4cTExBATE0NKSgoA3377LePGjSMuLo7Dhw8Dp18uzpkzp9O/LGbOnMmYMWMYP348U6dObSvgOnPmTM5NhRdCiM7odDre+vALmpuayUnfzqZP3yQnfTumHu4dOtP/KC35baKDXUhLfpsPXnmIrZvXS8fYM/7+j3+h+ETS5BjMrrwyQhOvI/u4jqJ6B3wiZrMrrww7/1hqtf7syisjYsYvySyo5WSLGx6jLmNXXhlOQZOoNPmwK6+MG25/sMvxsrOzWbFiBbt372bv3r0kJydz5MiRtu8/9NBDZGVlkZWVxfz58wH461//SkpKCn/729948803AVi+fDmPP/54l+mSH330EXv37uXWW2/lz3/+swV+WkKIS11rayuVlZXA6U6vhWXNBEXP4sNX/3ReMOlJUGrrf5S+mxdfeIGM9N0YdFVUV1VYNEmiOwM6IE2bNcei1ztyqOu+Rbm5uSQkJODk5ISNjQ0zZsxg9erVXZ5ja2uLTqdDp9Nha2tLfn4+x48fZ+bMmWbd0/Tp09sFPSGE6My29EP839urWZ9WxL8/2cJH67Z2GkxyM3eed77JaOxwJtVR/6P58+dRXVkuAekMOzs7i17veGFBl9+Piorim2++obKyEp1OR0pKCsePH2/7/htvvMG4ceO4/fbbqa6uBuCxxx7jlltu4fnnn+e3v/0tTzzxBMuXLzf7ntavX090dHTvHkgIMSg0NjbyzoeryD2pMiQskYLSOhz9whkSOKbTYFJyNK/dNTpbljMZjR32P0pO3sCpsoKBtw9JUZQHFEXJVhQlR1GU89a9lNNeVxTliKIo+xRFmWCJm9uc3PXspKdmzem6LUVERASPPPIIc+bMYd68ecTExLTVxrv33nvJz88nKysLf39//vjHPwIQExNDWloaW7dupaCgAH9/f1RV5cYbb+RXv/oVp06d6nCsm2++mZiYGL799tsuW6QLIURdk4nCWlvoYLbSUTBJSdnI8BHta2l2NZM60/9owsR4HnnkUSZMjMPd1YkRw4MG1gxJUZQoYDEwCRgPLFAUZdQ5hyUBo3/8dxfwL0vc3PW/utMSl2mzKXlVt8fccccdZGRksGPHDjw9PQkLCwPAz88PrVaLRqNh8eLF7N69u915qqqyfPlynnrqKZ555hleeuklFi9ezOuvv97hOB999BFZWVmsXbuWwMC+15wSQlx6VFXlqy3bWLPjED7DRnZ4TLtg8uijTJgYj42TFxGx09od19VM6kz/o8SFd5JxOJ9rJsXy1YE8Hsnej3YgBSQgAtilqqpOVVUDsB342TnHXAN8oJ6WBngoiuLf15v7Yfe3fb1EO+GR47s9pqysDICioiJWr17NL3/5SwBKS0vbjlmzZg1RUVHtzvvggw+YP38+Xl5e6HQ6NBoNGo1GWpgLIXpNbzCReayRJkPnQeGnYLKY7KJGEhcu5pY/vYrmnM4H3c2kCvLSGTPOj3se+z2/37sfp7o6dCdKMfbj7zBzAlI2cJmiKN6KojgB84Fz/6QfDhw/6+viHz/rkyPZe7CtO0xC+BCK0tcRPtSGkV4mTuzdSEL4ENSKHBybjpEQPoT871cRFehEkGszFblfkxA+hNbSTNyNpSSEDyF3+8fEjeu+HcR1113H2LFjWbhwIf/4xz/w8PAA4OGHHyY6Oppx48axdetWXn311bZzdDod7733Hvfffz8Af/jDH5g/fz4PPvgg99xzT19/DEKIQeiL5BRWrP0eO89gbGy7fp+u0WqJjJvBnOvvJjJuxnnBCDqfSYWNS0CjMRIVOxT/AC9QFNKee5oDd97G+t/9Ho2zs7Ue8TxmtZ9QFOUO4D6gEcgBWlRVffCs7ycDL6iquvPHr7cAj6iqmn7Ode7i9JIeQUFBE48dO9ZuHGmx0Dn52QgxOBgMBuqajLyXvBtV64KNrWUKpcJPm19LjuYxfEQ4EbHTOLx/B6Gj3ZhxKI/Cq+aBVktjfT01Vc3UFp7i8Qfu7PPmWIu2n1BV9R3gnR8v/BynZ0BnK6H9rCngx8/Ovc5bwFtwuh+SOWMLIcSlzmQyUVtbi4urO//+z/+o0fjh6jXM4uOcmUlFxs2gsb4GW9smZs4YzdRnnmX4jp24FRSy78H7qSxrpeRYM/FR4/u1/YRZAUlRlCGqqpYpihLE6fdHk885ZB3wW0VRVgIJQK2qqqXnXkcIIUR7qqqyPf0QyalfMiJmNk2u4bjYOVh93KpTOfgaTFzxxN/xysml1dWFnMhwstIO4Ts8nqDRCr5+Lv2aZWfuPGyVoijegB64X1XVGkVR7gFQVfVNIIXT75aOADrgN9a4WSGEuJQUFhby2cY07PyiCBo3C6NJxc7e0aJjnFujztN3KP4BzlwxwpVpD/4Zl5ITNPoPZcerL3LcLhjbCj1wOgit/t9/mD7uybbMPGszd8nusg4+e/Os/1aB+y14X0IIcUlTVZVjNVpa7Idi2RIAPzmzGdbYVE1S0lxSk98GOzueuf5qZj33Eva1dVRFjOGzB37L/oO1jI72xsPnp/N/cdu9/RaMQKp9CyFEv6upqeGNd1biMmIqrmdFgI4qbneUMWeuszfD2trasnzZMibGx9Hw7vvY19ZRMjWRrx97FL1hOMEdVAs/mLOXKWOH9Nt7JAlIQgjRzyoaFfCKaPd+pqPZzJ4tqzrcU2SuDjfDzktiY9kJYmfNIDXxMrZ9tBFV1XQYAKurKqRB3xnNRh2rkz9j8/ZUalsrWf7iUoorC9l3KJN/vv0Gta2VrFz9EdvTtlDbWslflj1JWf0J0ven8fYH/6a2tZIPPnmP7374htrWSp78yxPdjtkf7Sf0ej2PPvooo0ePZsKECSQmJpKamgpASEgIFRUVff3RCSEGqMzMLNZsScfVw7vd5z0pkmqugJERbNiQ0n4z7MZUhowayaYZV/Lxv97jeM73nVYJn3zZ5W3l0/rDgA5Ib775L+JnjidiYih1rZXcev9NqPateA5zZuFNc6hrrWTqnHhGRgVS11rJ4oduoUVpwG+EF3N/Nou61kpmXTWFwDA/6lor+fOTf+hyvP5qP/HUU09RWlpKdnY2P/zwA2vXrqW+vt5CPzUhxEBW1qihVT3/rZG5RVLN1VBXTeSEcTi7uzJp3Dgee/hh4sePx8nDncDQGHZuSseo6zoApn7xKQaDoXcP2gsDOiDNu2quRa+Xk33h20/odDpWrFjB3//+d+zt7YHTdfJuuOGGXj2TEOLiUXjsOMdqtLi4eZ73PXOLpJqrtiIPp+pskisrWZaXh8Prr3FVfCwzrrkZrd0oWpp13QbACZOm9Os+pAEdkM78wraUI4fzu/x+f7SfOHLkCEFBQbi5uVnmoYQQF41P131JZXVNh98zt0gqdN7XCKC06DAmQxGXBzvxyyXL8M07yKxh/sS/vwL/6+7FyfV0bWxzAqCLq/vAqvZ9IX22svvq3D1xzaKFXX6/P9tPCCEGl6ZmPe4jE3Fyce/w++YWSe2qrxGoePkYGFW0l1l33IvziVKqxkbw0ZKnSM6twmj0wMPbDzAvAKZ+8Um/tjA3q5adNcTFxanp6e1K3Z1Xr622tZK61kqLjblp9VbuvOVus49//PHHCQgI4L777mv3eWFhIQsWLCA7O7vtM1VVmTt3LitXruR3v/sdzz33HIWFhWzevJlnn3227TidTkdgYCBHjx7tcJYUEhJCeno6Pj4+7T6XWnZCXLwMBgNPPPsaQ6Ou7PMSWE76dtKS325L5dbr9UyYGE/I+BlcNieBkRFBJP7pMYZv30nJZVPY8ueHqVODaNIZsbNvXwGio9p2ZwfAiaN9mBY1tE/3C+bXshvQM6Rvtvc+u6QjMbEXvv2Ek5MTd9xxBw888ACtra0AlJeX89lnn/Xp2YQQA0/ewSOs+HA1a78/jltQnEXex3SW/KDV1BAYejp47HnmKfb/9h4+vfe3pGWcwGi0PS8YQfdVwrOz0vt1hjSgA1LaN7vZ+30ubnbe/G/Fagz1Co0Vrax6fz1udt7s2baX3Ix83Oy8ef8fK1Fa7Kg+0cj6lZtxs/Pm2817KMg+jpudNyte/YDRoWO6HbM/2k8sX74cX19fxo4dS1RUFAsWLJB3SkJcQr5PS2PN1n2k7K2h3DSEkopGXM5J8+6tDt/9bNjA7Lpq7H5831NjNLJ5ymy0jlGEjTv/HZS5mpubLHLP5hrQS3biJ/KzEWLgMxgMNLaYeGfVDlo0Tjg5W/4PzTPvkPSNlVx1VRKpqan4l5eTUlpK0TVXkfHko5woqudYfi1Bo6K6v2AXZMlOCCEuUv/3xtu8uToNretQqwQj+Cn5IXzyTA7m5fBYVRUppaUY3N3ISkxgz44cVCWwz8EIYN3nH/brPiQpHSSEEH3Q0NDAx599gWdoItphk7rt7tpXpUWH8fWzYfHcSUx74DOcT5TSMHwY2199kRMOI3Eua7bYWFNnzunXSg0DLiCpqtqvee8Xgwu1rCqE6JzJZOJwYQn7jrdyvNkLXUVjj4NRT4qpnv7dCL5+KmOKs5m19Dns6huojBrLqgd/R+b+csJjY3B1dyYnfbtFCrT2ZzCCAbZk5+DgQGVlpfwCPouqqlRWVuLgYP2GXUII85WVl/PGf1ZzoqoJn6GBPf5Duuv9ROc7mpdOU91egkcNYfzmr7Crb6B4+jTWP/ss+iGTCY2c3uNrdmfbl8mDdx+SXq+nuLiY5mbLTTkvBQ4ODgQEBPRrXxIhROdaW1vJyK9md155r6/R2X6ixIWLiYyb0XacQa/Hxs6Eq0sljk422NrZoW1uZuTqL9iWMI28A5WEx0zv0TXNFRfmw9TI/ktqGFBLdra2towYMeJC34YQQnTp409Xk1/ris+w4F5fo7P9RNlH89oFj6LDuxk2TMP0LT/wWUAABUePMjQwgIMxiTjZj2XM+J5f0xxVZSeoca0BCwQkcw2oJTshhBjojCYT+I7D2z+o7bOuast1prtacg21VaimKqYnBnPz2++w5N0P2Pj+R4R4+PD1p2t5/6+vYTSY2i0VWrJAq5+vB2NDh/f4vL6QgCSEEN1oaWmhtbWVmto6nn7xn9Q0GtsCQW/f23RXS05XfxTnqn3Mvv8Bsr7bRUlAALv27uWFF14gI30PNmrLeb2SelKgtcvnbdYRFzG831esBtSSnRBCDDT1TXreWPE/nH0CwckP2yGR7WYlHbUJnzAxntzMnV0uk53ZT5SbuZPso3kkLlxMROw0ThXn4+WrYbqfhmkPLMf55Cl2ubly5aJF3S7FdXbNnmbZ1VeU0HBSBxH9O0OSgCSEEJ0oLCzkyx+KcQyYiKoooIK7l2+7Y/ry3uZMLbnIuBk/Zher+A3TEHbsALOWPo9dw+m07prZs9m6YQPLli1rS1ZISdlI4sLFXV6zt+bNTrRIMkNPyZKdEEJ04tipOo6XNXSZ0m2p9zaFeRk01vxA4EgfIrd/g11DA8dnXMYXy5cxetZt2Dj59HkpzhzHj2SjaSyx+HXNMaDSvoUQYqA4UXqSTftqqdPpuzzuzDskg66K+fPnkZKyEa2jJ3Gzf8aJwoPdbk416FuxsTXh5laNg6MWWzs7NK2tjFi7np3xUziQXU54zPRuW0VYSoCHhstjh7UVlrYEc9O+JSAJIcQ5VFVl+V/fxGboeLNq0p0dLPyDw0jfugZTUzVJSXNJTd2E1tGzw2Z7AIV5afj5mbgqI4MjN16H0cmJ+to6ThVX4eQ+HlW1buUaXWMdhtZW3Dx9KM7P5s6fz2BMiL9Fx7go9yEJIcSFcvz4cfz8/LCzs6NO14r3mBnoDSazzj37vU1O+nZMZiQ51NdU4uQC0xICmfLUEvy/24Vn3kHSXlxOk86WinJbAt2sG4w0GoXxAXbYKyqXTY3gyCENvm72Vh2zKxKQhBCDXkFpHf98fwMJUy7Dlla+/HIzYZOu6tW1zE1yaG0qZqi+mtlPvY7nwcO0eLizO2ku332VSfCYWQSGWnaWcq7SYweJC3Vn/swr2j4bO3asVcfsjiQ1CCEGtaKiIv63OoXAyGmcqFUprLFh1MR5vb5ed0kOpccO0dRwkGneJn75l2V4HjxMfWAAm//9BqfCZ+PlF22R5+qMqqpoFBN3LJrO3JkJVh2rpyQgCSEGNXsnN0wOQ4DT74IOZOzgy1UrzK64cK7ONqeGx0wFTPgH2jKuNJeZd96H06kyKqOjWLlsCeuzSmnVO+Pq4WPhJ2yv4kQBni1HGB3si6urq1XH6ilJahBCDFqtra1sy8gn96SxLVvOaGYyQlc6yog7djgT7yFGIsaPIfaFVwhd9QXHZ01n0+9+T7NNKC1NKjZWLqCsmoxcPWUEw73ssbfvv3dFkmUnhBDdKCkp4ZV3kwmJmtphpezYCXEMj0hEq7XpdW8hg74VrY0RD8867B002NraohgMhKxL4fu4yeTsPdlWrduampsaaSjYyTOPP4RG07+LY9LCXAghumPvQUjUVOD8ZASNRoPJqOd4zs4+9RYqzv+BzK3/IePhh8nZuh2j0Uhdo46No6MwaUb3SzAyGg1MjgriiT/d3+/BqCcG7p0JIYQVlZeX885/P237+txkhPXr16OqKnuzMnnxhRfISN+NQVd1XkHTztTXVKJvPsH3G/7D/q+24jLrclI+/IRXH3kKXaMtNdWOVt9jZNDrsdFqKM3+mugA+wHf6FPSvoUQg0JLq5FTZSfxcHPFzc0NJxdXbL1Gtn0/InYae7asYsLEeObPn8d7/3mfW2+7pde9hQz6Uo6lb6W1opJd+/dja2vL0mefZXxMLOs+XIVWa0NtZZnVKi7YtZaj1Bxn8R03o788GGdnZ4uPYWkyQxJCXPJ27NjJijXf8c9Pd7Ju214qa3Vs2p6Ok5t32zFnKmUnLlxMdlEjMTOuYePGzT2uUVdadJjG2gNMcWvF8/2PmHv11e2WAVFNfV4G7EpDbRVupjLuvn4Gv7vr19jZaC+KYAQSkIQQlziTSeVwpQad0R7/EZFUGj1ZsS6Tbbtzzzv2TMWFOdffzbwb7u1Rb6HTCWImhgfZEVN2iFl33sekxkY2rV9vsWVAcwT6OjMp3A+NRtNlUdiBSJbshBCXrPz8fDZ+l4vBJQR7x58+d3RxJzQ6sctze9pb6NjhTDw8W4icEMHoz/Zh29hI5MzpKE3NvVoGPJM6XlyQa1aGn8lopLb0EHfddS2O9tZNH7cWSfsWQlxyTCYTzc0tfJtzgu/2FuLh7We1sQz6VrRaA17ejdjanw4yGI0Eb9jIntgE9v5QjGpSKDmah8Ggp/Tgrnap5RMmxpO4cHG7gNSbPVE2GhOhjhXMvXLWgJsZSXFVIcSg9UXqVn44XIZPcLRVgxFA6bF9eLs3EZ+SxuFfXE+rhwf1DY2kBoXhoYQydsIoACLjZvwYaA61zZhSUjZ2uAzY0y60ddUVLJw2mpiwcVZ9VmuTgCSEuGSsT06lWhlCZZMP3kFDrDpWfU0FDg56pkwcypTHn2LornS8sg/wzRv/R0uLIw11rrh7t5+pmLsM2NMutPbGGpxotN7D9hNJahBCXDJq8KKs4fQv/s6WrUxGIznp29n06Zu9rld3+jqncK/LZfY9v2XornSavTz59qbr2Z66m5YW306z8c5OnIiMm9HhElxPutAaDa0s/kUSYWFhvXqOgUQCkhDiomc0Gknd8i0VLY7Y2nVeo+3Mu5m05Ld7nXZdWnSYuur9JDo3cfNfluNxOJ+64CA2/fM1ysbMYmhQfJ+fp7MCrWNjp2Ey/dSjyWQ0Urx3E442FyYXwNIkIAkhLnq6piZ2/JDf7cv8s9/N9DTtWlVVVNVI0EhH4suPMHPx/TiVlVMRM47PnlvKFxnFtLQ44uLu1efnOXdPVOLCxUxfcDOzxzrjWJ9HXVkhAB6uDjz31J8GfAUGc8k7JCHEgFdZWYmNjQ3u7u58tWULUxITcXJyAqCxsZG9BdUMHTWh2+v09N3M2YoOZ+HqriM6LpKQVQexbdRRdPlMUu6+B6NDOFGTxvf9Qc9ydhdaG61YBjNpAAAgAElEQVSGSSE2hI0MIHxUEAajia935XKiMBsH+zEWHfdCkoAkhBjwVn+5m2qdgrtfMIdyyjipluBkKGO4lz3NrQa+2HmY4DGx3V4nYGQEqclvs3zZsra065SUjSQuXNzpOfrWFrRaPbGJgdjYnl4ay77vLmpDR5IxLo79P5QQHqPFCtV/ADhReJApET7Ej/spYNoBs+LDaB0fYp1BLxBZshNCDGgnT5XT4hCAvWcAza1GgkZHU1nfwsEKhW0Hmyhs8iIoLMasa3X2bqaz6gsAZSUHaKrMYNI7K3CprQOgvr6BFJ8gDOoIwmMus8hzdubnc+KYnnB+OrezszOenp5WHbu/yQxJCDGgffpFKjXaQNy9fNt97uTsBkBNQ6vZG0F7Un2hrroCO/tmEmN8mPLY6/jtTsdrfw7b3nwdfaszzS2eVq3W3dKko6Ush7hrbxtwG12tRQKSEGJAGzFuJgWldRa73tnvZro8TqnAu/4Es55+DY8jBTR7ebLz1pvZuiGNMTHzGRbs3eX5feXj5cLkCVMGTTACCUhCiAFs+zc7+SFXh4fvsH4bs7ToME5ODSQ6G5n2yDIcyyuoCw5iy4vLOeUWQUB5s9Xv4WTRES5PimVsWLDVxxpI5B2SEGLAUu29cPhxac7qY/2Y1h0c6sykyqPMvOt+HMsrKI8dz6rnl7M+o5iWZgecXT2sfi+RIe74e7tYfZyBRmZIQogBqb6+nrIGLQ5O/fOL+Xj+Xpyc6xg/aRyByUXYNuo4dsUsNty5GNVpLNGTrF8nzmQ0om2t4pc3zUSjGTxLdWdIQBJCDEh79+eQtjuf4DGW3d9zLn1rCxqllQmTA9H+mNade8et1I0IZu/YWPanHyc8RgucLjlkbjuI3jC06hhmV8sgem3UjizZCSEGjLPb4bj5h1k9GAFUnMyjqSqD+Df/iWtFJQD1dfWkOA+h1RhMeMxlFik51J3mpkZmx4Vy48+vHVSJDGeTGZIQ4oL77vvv2bqvDC8XO4Z72jJu7CiSN/2As5/1CobWVVdga9tIQrQXUx59Fb89GXjvz+Hrt/+JweCC3jikLa27p+0gumLQ62mtLsQnMJx6XQta7elfwzXF2TjFugHWrVI+kMkMSQhxwRkd/dE6D6HJ3p+8Gje2H6jFoLVuMoONTQ0+DUeYfddv8duTQbOXFzvuvI0t679Fp/PCP2h027GdlRwqOZrX43EVDMSP8WFOtDtFe7cAoFHg4XtvJjy84wrhg4UEJCHEBbX9m+/JOFyBo7MrWq0N9o5O1LRocfcZapXxSosOU12WQZy2mpuXLMc9v4C6EcGk/ONVykbPIjhs+nnn9KQdRFeadQ0kRAzlipmXERw4jBeeuJ8ZUb5U5X2Nh4v9oF2qO0MCkhDiglFVlfSDpRhV6/8qUlUV1WRgRJgLiVXHfkrrnhDDquefZX16Mc3N9ji5uJ93bm9KDnVEV1WM0lDc9rWDgwNjgz255zc3DvpgBPIOSQhxgZhMJrIPFWHrHUZ/dPMpLtiPo0M14yePx//kKWx1TRTOuZzk236D4hTFuMToTs/tScmhzphMJm6/YR7Bfu3T2O3s7PDzs26b9YuFBCQhxAVRWVnJ2x+vIzR2tlXH0be2AM1MmByAxuZ0xYdDv7qJ+qBAssPHcWDPMcaM16LpZpJmbsmhjqiqSkF6Cp5X3teLJxg8JCAJIfpdU1MTh8uN3QYjk9FIbubOPu39qSo7jJNNOfFrtnPkputpDBhOfV09P9i4MVQfyJjxQX15lG6pqoqLoy1P/PEe3Nz6p+rExUreIQkh+t07H3zC5p37uzymr3t/aqvKaaw7SvxYN+565z+M/mQVkx/7C6rJhMHojqoZbtVq3WeUFOQwTHuKoT7WLzl0sZOAJIToN6qqUlPfhNY/Dk9f/y6P7Uu7cQAHxwb8dPnMXnw/QzIyafL2Ytv99/Dl2h3oGt0ZGhBqiUfqmqpy66JZTJ8y0fpjXQIkIAkh+k3ewYM8+9r7tJoxyent3p+Tx49QUbqLCWo5v/zLctyPFlI7IoQNf3+VstDpjBhr3XdWZzTUVtFYlEbkyCG4uAy+Qqm9IQFJCGEVLS0t5ObmUlXfwqrkr/hozZccqnJgyKg4s87v6d4fVVUxmfSMDHdlSlURjbffyysVlXwYOpJVzy8nOeP46bRuK1cPNxmNNDXUcHnCGBbfPHjLAPWGJDUIISxq//79KHaunKhTWb9hO6HjDLS2uKHRaLE91YCdvYNZ14mIncaeLauYMDGe+fPnkZKyEa2jJyaTiU2fvnlekkNJYQ52NuUMSxjHSx9/Sl1ICFdeey0vbdyIzWvvcPODf7XmYwM/JjBoGhlqX8bkiJ7tURKgnF3MsD/FxcWp6enpF2RsIYT1rPh8B6W1xvNajvfGmSy7kqN5+AeHkb51DaamapKS5pKaugmtoyc3/f55NJpWAoJA0RjZ+/0uUj74H+m792Bra4ter2fCxHgSFy7uVcq2uUoL9hMf7s91V8202hgXK0VRMlRV7XZqLDMkIYRFGI1GNm3fRaPGE3cvyyxTnb33Jyd9O6YOCpxmfrOa6KhhTPp8K0du/DlHDx4mae68du+e5s2bw9adqVZpHVFyNI8rp8dx84wk3FycLHLNwUreIQkhLKKhsZHt6Yetdv3Okhy0zaXc9e57hK76goTHlxA8KpSUjRvb3j01Nzfz/nvv01xV1OfWEaeO5xPmA8OdGqk6/A1RIZ7MiPQmfpQHPl4e2NnZWfSZBxuzluwURXkIuBNQgf3Ab1RVbT7r+7cBLwMlP370hqqqb3d1TVmyE+LS0dzczN6j1aTlVVptjJz07aQlv902Q9Lr9UyYMIGnyk9xw6lymny82bjsKdblV5H5zTcYdVXMnz+PlSs/xcnRgX379vZ5Cc9HW8mcxLH4+PhgMBjagqPomsWW7BRFGQ78HhirqmqToiifAjcB751z6Ceqqv62NzcrhLi4ZWTu5fOtuYREmJdB1xtnkhxiJ0zkqvnz2Zi6geCiIq6rq6N25Ag2L19ClUcU4a56xiX8rK3unG/wWC5PjDpvZpV9NM/sgGQyGlGaK/jFL2e3tRbXaDQkJyeTmZlJbGwsSUlJaC3cQXawMXfJzgZwVBTFBnACTljvloQQFxNVVWmyDyA43LqbPxWNhl899CKzbriJUxm7WZKbx+K6Op4Y5s8zN/yc9XuKaGqyw8HJpe3d05zr72bCtKQ+t44w6ZsIdNG1tRY3Go0sWrSIJUuWoNPpWLJkCYsWLcJowQ6yg1G3AUlV1RLgFaAIKAVqVVXd3MGh1ymKsk9RlM8VRQns6FqKotylKEq6oijp5eXlfbpxIcSFp6oqr7z2JukHjll9v82JogNUnvyOGUnTuSdxEu/Z2fF0RATqr2/hq5Qv2ffdVlST6bzz+to6orG+hrmTR3HdtQvbnjE1NZWSkhLS0tJ4/vnnSUtLo7i4mNTUVIs+82DTbUBSFMUTuAYYAQwDnBVF+dU5h60HQlRVHQd8Cbzf0bVUVX1LVdU4VVXjfH37nhIqhLiwGpr02A6biKMVN5u2tjTT3FRJbMJwxiWMB+ALL0+KQkPZvXdvW1khYydlhc60jkhcuJjsokYSFy7mlj+9anaWXUv5IVw1je0+y8zMZM6cOe2WAefOnUtWVlYfn3ZwM2fJ7grgqKqq5aqq6oHVwJSzD1BVtVJV1ZYfv3wbkMJNQlzidu3ezesfbECxtW6qc0PNMfRVWUx+8QW8CwoBOLgvm7lJSeeldv+wM5VNn755OkX8rOWzs5fwIuNmmB2M7LQqj/72FoKDg9t9Hhsby+bNm9stA27atImYmBgLPPHgZU5AKgImK4ripJyer14O5J59gKIoZ1dJvPrc7wshLj2VRi+0LtZpMw5QW1VGbcVBYkc5cve7/2HEFxtIeOJpVL2egNDxpKRaJ7X7jIa6aqoPbcPB7vzglZSUxPDhw0lISOCxxx4jISGBgIAAkpKSej2eMD/t+xngRsAAZHI6BfwJIF1V1XWKojzP6UBkAKqAe1VV7bICoqR9C3Fxqq6u5j8r16EMiUXTXVe7PjAZivEoy+W6V17DrfAYTb4+bF76F9YeqSA6YRH//esfMFghtfuM0cPduSLGr9O9RUajkdTUVLKysoiJiZEsuy6Ym/YtpYOEEGZpbm7mxKkKsk8YyMw+goePdWZHp4oLUI0lXOFly9SHHsahqpqaUSPZ/MwSqr2iqa3SY+/o1K6s0KmSo1yeGMWLL7zQdp1HHn2U7KJG5lx/d4/voSA7jbuuu4zIsWMs+WiDlrkBSSo1CCHMsn13Dv9cuZWjJxusEoxUVcVoaGX0WDcubzzFjLt/h0NVNacmxbHu5RdI+eE4TTob7B1Pv7OydGo3QGHuHvxc4I4b5jIqNMTSjyi6IbXshBDnObMclZmZybBhw1BcA6i3HU5gWKzVxjxVfBDVeJygEfG4KAo2zc0UJM1h9Q034aANI2ZKeKep5R1VBjc3tbu24iT2raXcdesNnAzXEhwcjJOT1KS7EGTJTgjRzplNnyUlJcyZM4fU1I0YtE7c8dg/LVaQ9GytLU0YjQ2EhNqhoEdrc/rvZO+9+zgaGEbmrqOERk7udp/T2Ut4w0eEd1tA1aDXY6MYWDBlNK62rQwZMsSizyV+ItW+hRBmMxqNvPXuh4RETWPX16s5fvw4u3efrhm3dOlSJkyMJzdzp1XaN+jqS7BpOcbkZzeR//NrqY4cS2N9Pd+V6QjyHMKoKPOWB8+uDG7ewKWMGW5P6HCPPty9sCQJSEIMciaTiRNVTVRpA2gp1bFjdw5z5p5fVbsntd/MUVtVhslYQcxIF6b9+V18svbhvXc/qf97D6PRE2ePsSiKdV5zG/Ut3HPTlTg7SHHUgUSSGoQYxFRV5eXX3mTll/tx8fRDURRGRcVbJEGgO+4eRoL1J5h95734ZO1DN8SXLU89yoZV26irc8ZnaJBFxztDVVWK9n1Fa1ODVa4vek9mSEIMUqqqcqysAY1/HHb2P73E70uCgDlOlRRgbD3OHG87pv5l6Y9p3aGkPvMUdV7jGOvcwoGMHVZppgegKAqP//F3eHpK4sJAIzMkIQapr7dt4+2Vm9oFI+h77bfOqKqKQd9M2FgP5jSVM+Oe359O606IZ/0rL7Lph2Ia6xVWvv44aclvW6ziwrn3UHZkNy72FrmcsDCZIQkxSGk8QnH3c+n4ez1NEDBDWckhDK2FBI9MwMHBAW1rKwXz5/L5dT/HSRlFzJQwDmTswNhBm3JLJlRMixuLo6OjRa4lLEsCkhCDjKqqrE3eSH6jD3b21v/F3NrShEFfx7j4oSiKDwClM6bx9btvUjR8FMVpBYz00KAoSqdtyi2VUKG01jF3RqLVW2WI3pElOyEGoYIKE1qbjmu0WVqzrhRTbTaJS5fit3c/AI0NDaw7XkVtnTehkQltASJgZITVEiqMRgO1RZkYDIZzPjeSnJzMsmXLSE5OliZ7F5AEJCEGEb1ez9a0fdh4BFtlk+vZaqvKqDy5n/FBWu5++x2CUzcT/8xzmJqbMRo9cfcdf15ad1+b6XWm8lQxikHHo3+4t12xVOn8OrBIQBJiECktq2TdV7utPo6qqnh6mwg1nWTWHffis3c/uiFD+GrpkySv2kZdrRPeQwLOO8+SCRVNjfWU5H7H1Eg/Lo/25meJw7GzaX8d6fw6sMg7JCEGAVVVycnJY+8pG0ZEJVp1rLITR2nVHSXJ154pTy3FobqGmlGhbFjyJA0+44l2NqKaTBzI3NlharelEioiRvoTMcmHsDBfCOu4Q3VXnV8XLFjQp/FFz8kMSYhBoKm5lQ/W7eBkZc83g5qMRnLSt3fYifVsqqqib20iLNKTpJZKpt/zAA7VNZycHM+G/3uJL7NO0KSzwcbGlg9eechqqd0mo5GCzK+YFe1DWNjoLo+Vzq8DiwQkIS5x+7Oz2ZBWyPCIqdjY9qxUjsloNDt4VJTmc+r4dpxdwNbDHY3BQP6CJP7zmzuoMo0kdurVKIpCbubOttTuF194gYz03Rh0VeRm7uzzs6qqip+3Mw/c/jNcXZy7PV46vw4sEpCEuIS1Gox8sXU/R0oqe3W+OcGjtaWJxroTRE30Je6yeABOTZ7ElvdX8M2DD1NaYYfR+NO7m3NTuzUaDWFho9i2/r9dzsA6Y6g+iqH8AONGeHEyewuTgm0ZGRJkVmq3VqtlzZo1LF26FGdnZ5YuXcqaNWuk8+sFIu0nhLgEqarKZ6vXUqMNokntvizBmdYN577T2fTpm0QHu3TZibWh9hhKw2F+lbyBwoXzKUuIR9fQyJ7tGYRGLzgvky4nfTtpyW+Tkb4bjUbDNddcQ35+PldffTWpGzejdfQ0K5HBZDTipNFx21UTUFUVBwcHmpubsbe3l31GA4x0jBVikDIajZysbiK/1pUGQ/d/6Xe1LNfVvqDaqjLKSjKJDlC5a8U7BG36itilz7Pnq69Z8/5qKssNqKbz/+A9O7X7hhtvJD8/n3379vHiiy/2aPmuRVeLq74Ee3t7HBwcAHBwcJBgdBGTLDshLjFvvPkfGl1G4+xpXh+hs5flzi3X01mh1fCYqdjbV+NSUsWsOx7HtaiY+iG+zAkMpO7Dz1iw4CpSUlaT9uXnjIm9jMDQsW2zrjOp3bmZO9m2/r8suvrqHldmaG1pImlaJONHTrfIz0wMDDJDEuISYTAYOFWtQx0Sg5Orl9nndVaup+RoXof7gpJufoDi/K+JqDzETY8vwbWomOqw0Tz9q1/QaDCRlfkDzz37LMFBAdhr9IzrIBniTGr3zIW/JnXj5h5XZjh1KA0vW10vf1JioJKAJMQl4osNG3njw00oNj1btuquXM+Z4HHFdYsZHR1HeLQ3C4x1zLj399jX1HIyMYENf32ZrbsOkDTvdGBLTU2ltLS026U4cysztDSfDj4l+ftx1uh46g+LCQwM7MuPSwxAEpCEuMhVVlZy/EQZlZpgvIaF9vh8c4NCdVkhp4q24+QCGl8fMKkcWZjE27fcSrUxhLgZ17YFtszMTK688soOZ11nM6cyQ11VOZTv4/a5Ydy9aBLXzxyDu4tDL35SYqCTLDshLkKqqlJeXo6tkzsfrd7MqQYFH/+QXl/vTJZdydE8ho8Ib1c5oaVZR2tTJaMj3VEUAxrN6b9j3fKPUjokmMxdRwkeHQuqygevPIRBV0VY2CgO5OSwb98+bG1t0ev1TJgYT+LCxT2qwGAyGpkWPYzxI9zb1aATFxdzs+wkqUGIi1BdfQMv/et/DIuYjtZ1BD6ufbteV+V6jIYqNLqDTHlyHUVJcymeOpldX29n5+btjJt8LZETp7cFqTPJCsX5BzBoCvrUdVZVVQoyUrljzn0SjAYJmSEJcZGpra3lh4J69hXWWHecqjKaGopIHOvN1D8+hvf+HBp8vLk8fAxN9Y0kzZtHauqmTvcNnT3r8g8OA+BE4UGz2pKrqoq3mwML4v3wdHez6nMK65N9SEJcolK/+obNOzr+Y87cunPdUVUTvkO1jLWtYfbt9+K9P4fGoX68ePNNNFTXkrFnT7dlf9qSIX52Jxlb17A75V2za9cdP5xFgH2FBKNBRgKSEBeRxqZmdC6j8Qs8P3mhJ3XnulJxsojCg5sIP3mAGx5/GpfiEqrDw1j14rPsrTSxYMGCbpMVztbT2nUKcPdN85gcJwVOBxsJSEJcJOrq6nh82d+oa2ztMK27r0VLTSYTLc0NhEW6c42xken3P4h9bS2lUyeT8teX2Z5TwZCAnnd07Wqf07lqq8rQn9hDWLAvjo7Wb68uBhYJSEIMcAaDgUOHDlNYoScoZk6ne4x68ou/IzUVxykr2o6jk4oa4I+KwuGrr2LFL35FlT6YcZOvInLi9HYp4rET4mhsMXE8/0CnS4TmtiVXVZWZ8RH85hfX9OTHIy4hEpCEGEB0Tc1s2fI1cHpGVFdXx4Gjlfz9o01sySzBxrbzbDNzf/Gfq6VZR23VMcaOd2fSrNPN+6qiIvnqo3f57qFHqGlwx2A4HQTP3je0v7CeFpMWFwct40JcO10iNGefk8lopDBzM/GjPfDyMr/KhLi0SJadEAOEXq9n3XdH+W5XBuFRsVSX5lPf0IBPUKRZ5595h2TQVbVLte6ucnZz4wlMNQf49eq1HJ9zOSVXzKJJp2Pn5u8Ij12ERtPxuWdX7e5ur1Fn+5wMej2NNSdZdGUCvo6tDB1qXv09cXGRfUhCXERUVeXZV97AIXAyQaOj0bUYsPcKxr4Hk4Wzi5ZmH80jceHiLtOra6vKaKw7yrTIIUx5cQXe2Qfw3refovgJ6PFl+IjpnQYj6HyJsKPCqJ3tcwr0sUdrqycqxFOqdAtZshNiINC1GPCNmI2jc992uJ75xT/n+ruJjJvRaTAymYwM8bcl2q6eWbffg3f2ARr9h/Llqy+xbu1OamsccPca0uVYvV0ihNPVumuPZXD9rAiuW3S1BCMByAxJiAuuoKCA/677Fo/g2H4Zr+LkcWor97HI35UpTyzBvraOqvAxfPHon2jxmcDEGRPbKi90pbPWFOZUY3BzcWTu7DgJRKIdCUhCXEAmk4kmxQ0bz5B+Gau1pZ6wKHcCdjYz/b6/oG1tpXRaIlsff4Jvtx8g0lljdvvuni4RnlFRWsSiGRGMixxhiccSlxAJSEJcIOUV1bz89xWETLwKFzdPq49XW1VMfeU+QsOmYQgKwmSjpWD+Aj6aPQevluFEJwT3+Jpd1cDrTLCPHX4eUptOnE8CkhD9QFVV6uvrcXNz4/td6ewvKMfoEoxvxGz0RutmurY069A1nCQi2guNdioAtWNGs/GDd9h6qIhdyTuImMDpIqlmzo56Q1VVTE3V3HrjbGxt5PW1OJ/8XyGElRkMBvYePsVjy1/n3U0H2XqwhRq8aTWYsLO3fl8fhXocWvKZ9uhjhGzYCEBTYyN/eOpZvlq5kmmxoexOebdXZYbOpW9tISLQjcoTBZy7pUSrtuBhOI5WfuuITsgMSQgr2rcvh5Wp3+EXlkhY4rXU6/Q4u7r3y9i1VeXUVR1ixjh/pi5fgdeBXLz351AwJZH0PYdwsLFp20O0fNkyJkyMJzdzZ4+W385VXnwYg3s54b5G6mxUWowKRoOe6tIj/PH2q/Fw6XYrihjEJCAJYQU6nY7jJyvJKNXiFTLBrHPObB4tLsg1q0VDl9cyGRk63JZRLc3Muv0eXE6U0jjMn69fWs669d9SVV7H/PnzzNpD1BPzLp/OzPH+KIpCRXUDL/zrE8aOn0TkWG/cnGx7fV0xOMjkWQgr2JN1gHdW7aROZ8DWzr7b4y1VqRugsqyYgpwURhfv47rH/4LLiVKqIsJZ+ezTlHjFMnHGjQSGju31HqLOnCo5ilZ3oi2V28vdiRvmxPKLK8KZe8Vss1LJxeAm/4cIYWEmk4kmWz+Gj4o2+5y+VupuG1dXQ1ikO9fZ6pl+/0PY19Zx4rKpbPq/v5J2sJ6mRg1arY1Z9eV6KmS4H6NChrV9rdFomBQ3AUd7mRkJ88iSnRAWtm5DKmn5zfgHjzH7nJ6U4elMfc0JasszGTXmMlpGjcRob0f+gnl8MH02Q5r8iZoU1HZsb/cQdaalWcf0CYEMG+bfq/OFAAlIQljc0NCJ+DaV9uicgJERpCa/zfJly9oKlaakbCRx4eJuz21p0tFQV0JUjA+K9vQMp35ECJs/+g/Vbv60pBdjUm04tyhCb/YQdaa2vJiqkkYYIwFJ9J4s2QlhQQcPHuS7rEPY2PZsmaovS2habQNOxmNM+/OjhK5aC0BLcwtrd2ZRXu5EwMhIq5fouXJGIpfPnmXVMcSlT2ZIQljQ8bIGGpoMuPew2WlvltBqq8qpqchlerQ/TX/6K+8UFxOZuRf9ZVNodQxmxNi5aLXWf39zsugIUV51gMyORN9IPyQhLKSlpYXU9BMcK2u0+lgmowFXNx2Oh9L4+C/LOOk7hCuvuYbUTRupbzZw39L/mpXd11uqqlKQ9TW/uH4RQb6OONpppbGe6JT0QxKinyWnbmZXQUufUqfP6GpPUnX5CcpP7Ob6YR5UP/wkJwMDScvOxtbWlmXPPsuEifEc2pdmkXdDHTHo9cydFIzHpOsYOnQoNjbya0RYhrxDEsJC/EbF4R8c1ufrdLYnyaDXo2uoInSMEwFHc9l+7wOsamriymuvPS87r+RoXofXzUnfzqZP3yQnfXuvywSVHNyDWn+CgIAACUbComTJTog+MplM/POt92jxiMbGvocvjzrQWWvw8bN+hq+flu83b6K5upZ5c+fy6cqV2Ds7s2/fvi7biJ8JcsamapKS5pKaugmto2e37c07cmXsMMYEupvdpkIIc5fsZIYkRB/V6vTonEehtbNModQO9yQlzcPQVIS9g5am2jr2pKfzwksvceDQISorK4mJndBldp4lNt4ClBxKx88VCUbCKiQgCdEH65JTeevzb7B38bJYanVHrcE3pGxg/rffUv/ZGpLm/lSDzsHBgVtvuw1H72CyixpJXLi43aznzDLdtnUfkDRvjllLe12ZER+Bu7ubRZ5TiHNJQBKil1r1Rk4ZfbFx8rbodc/sSYqdEMcjjzzCxLiJBJ48yc0H8rhifzapqantgtXGjZuZMC2JOdffTWTcjHbBqO1d1JhA1q1b16fadUZdBfNmxMt7I2E18n+WEL2QtXcvyTtycQ2IwsbCW300Wi03P/QSxw5t4+ierfzl+HGurq5htZcXq2dNp2JPFhMmxjN//jxSUjZ2uoH27GU6jUbDNddcQ3R0NNdccw0pqZt6VLvOaNBjrDyCyTQVrTQ0ElYiAUmIHjKZTByrd8bGs+ctv81RU1HKyaLvuTHAkynbdqCtb2CBry8lQ4eSNGoMXoeO0thiYv+xhi430J77LuqLL77ghhtvZG3qdmaYsfFWVVXsbODAt1/w2/vvIywgBo3GuhUfxOAmf+oI0Xz+HYIAACAASURBVAOlpaU8+cK/KCxvxtHZ1aLXNplMNNZXMCrCmZtcNFz2uz9iV9/Ah1FjOTFsGBkZGbz4wgv8kLEHZ3sNASMj2i3Rnevcd1Emk4lDh44wY+GvuzwPoPzEUVx1h7lnQRTPPv47woO8JBgJq5MZkhBdUFUVRVEoLi7m+MkKmrS+uAWZ13Cvp3QNZVSV7iFs7GU0REbQ6u5G0awZvF3bwLzY2B5XAo+IncaeLavMWt47m6oa+fXCKQx1t0GjUXB3758Ot0LIPiQhOqCqKq16A8+9/BpxsxZx+NhJqqpr8fUP6v7kHmpuaqSuspCYOH8UG4UzZbltq6qpd/Rh7X/XcWjXxvP2JZ2716gjZyo+lBzNY/iIcDPq453CreUYD953u0WfUQxu5u5DkoAkxDm2frOLPdkFuAfFUF1dY/GluXOppmpaT/3Arz9eSVVUJLl33oZer2fzqs1EJ1yPomhPV2rQVbWb6fRmU2tnaipOYksLNy24jGGedjg4WGZPlRAgAUmIHjt58iRHTzXy/aE6UBS0WuuuaNfXVFBWnMWc+BCmPvQInnmHaPFwZ817b6FzGUlZqb4tGPZ0pmMOk9FI3q4N3HvXHbja6lENzYwaNcoSjyZEO1JcVYge2ro7j+yieoYMH2H1sYyGVvwD7QlrVpj1m3twPnmK+sAAtr38HBu+3kfM1NE4Ov80S7FkM70znBzt/r+9Ow+PqrofP/4+M5nsITvZgYSdsGRjCYugyJKoIHXf9711r2hbtWJrsfZba7+1tUr7q361VaviggRQqiBqgEDY1wCShS2E7OvMnfP7I4EmIcskmQlJ+Lyex+dJZu6958x94nw4557z+fCLR+4gOqK/y+slCeEICUjivLdnzx6y9hVSZITQPyrE5e2VFZ+g4OBarosOIvWJX+BeUUHR6Hje/fF9mHzjSZ6e4PIAUXSigLFR7sREjnRpO0J0hCz7Fue9nfk1HCy0ubwdu91ORekJBo/w4jp/S/2y7ooK8i+czhevvMzOo5qqStUto5URsRFMGDvY5e0I0REyQhLnrYKCAj5ZvRGr3xD8g0Kddt3WahnVVJ7k1LH1DB89jfIx8VSHBJM3fRp/S0whsjiI4eNmOK0PbSk6cpDrr7+Q0EDXLtYQoqMkIInzVl6J5ni1N0FO/F4+q8zDsiVkfvE+adfdRcrEGIaOmAJATUgIX7y1hHLPYCxbTqCUu/M60QatNdF+VrzdZXJE9DwO/VUqpR5RSu1USu1QSv1LKeXZ7H0PpdR7SqkcpdR6pdQgV3RWCGcoLy/nlb+8yfq9xQT1j3TqtVsq86BrSji5fy3THnuC0X95AwCbzcbSjHWcOOpBWPTgbpmmO7j9W6YM8+Pum6/Ax8fH5e0J0VHtBiSlVBTwIJCitR4NmIFrmx12B1CstR4CvAy86OyOCuEspTVQYok6swHVmZrnjzOZTAwfNpScf7xNVuYGYj7+DGtuPlWVvoxImsferd93uYJre0qLjjM2LogHrp/FuGHODcBCOJOj43Y3wEsp5QZ4A0eavT8feLPh5w+AmUrWkYoe6MOln/L+6h0EBIe75PqN88cZhsH8+fPZuWMHqT/+Mc8MH86sYcP4/OttFBVq3v/TL84qU+7soGQ3DDi1m6mjQhk2JA539+6ZGhSiM9oNSFrrAuB3QC5wFCjVWq9qdlgUkNdwvA0oBZxbJEaILrLa7ORW+oBb18uMt2Zk4lS0xZeExESuufpqDuTksG37dhb/9res376d8uo6vHwj2b99g1MquLbH38+Tnz9+/5kRmxA9mSNTdoHUj4BigUjARyl1Y2caU0rdrZTKUkplFRYWduYSQrTLMAyWLVvG888/z7JlyzAMg3379rHkg/9g8glzWQYGu2FQWV7Iwpd/w40Tkin+6CPmz5vXJCnq3Lmzyf52hdMquLalvKSIwt1fYzbJAgbROzjyl3oxcEhrXai1tgIfAZObHVMAxAA0TOv5A0XNL6S1fl1rnaK1TgkNdd4yWyFOMwyDBQsW8Oyzz1JVVcWzzz7LpZfN47u9JRRWtv3nfrrcd2ef6dTVFlN8IgtPL4i//iruCQlm5eefnyn/UFNTw5v/eJOaU7lOqeDamMmksNvtTV4bGBPBvbc1f9wrRM/lSEDKBSYppbwbngvNBHY3O+ZT4JaGn68E/qPPVZI8cV7LyMigoKCAzMxMfvOb35CZmUleXj7fZW7EP6h/q+c1KffdwWc6NdWVFBzKJi5WkTojFYC6gADUO3/HEhxCYlIyC598kuEjRhEcHEz25k28/957DB48mDFjxrBw4UKSksd3qILraXW11eRu+w83XjiIvOwM6mprADh14gjh7mVSOkL0Ko48Q1pP/UKFzcD2hnNeV0otUkrNazjsb0CwUioHeBR40kX9FaJN2dnZzJ7ddCrskkvS2p0Ka2m5tqPPdHy8bYR7lzDtsZ8y7uX/Ba0xbDaWZXzLlff9lsnz7mFHbiWhA0cxb379FJ7ZbOaTTz4hfvRoPs5YQ+pld3Uoe7fWmqqKEi5OieXh2xcQ2M+H5xbezyWpg6koOUlEiB+Do4McupYQPYVDk8ta62e11iO01qO11jdprWu11s9orT9teL9Ga32V1nqI1nqC1vqga7stRMsSExNZtWpVh6fCmi/XduSZTkXpKfZkf85Ay3Fu//0rhGduJGbVamyHDlNZ6Uv8hCvw8PQhPmU6s6+6h6SpaZ2u4NpcbUUxwdYfSBwSQkxMDAB+fn6MiPEnwHqImSlxxMXFOXw9IXoCedop+pT+/cPA4kNS8ngWPvmkw1Nhzct9txfIrHU1RAzwYGaEJzNvv4+A/QcoHxDDF6+/Skbmfk6dVHh4ejc5Z2TiVMxegR3uW3M1VRX8aOY47r797LVFJpOJB+6+jaEDwzp0TSF6AqmHJHqt4uJivsncTMqESWRnZTJ29CjWbMnnh8Jqcvdv71DtoNPPkBwpgldZVsQPe77khuggJi98GktlJYVjx/DOPXdiGTAVa52l1cwLXa1rpLXm2I4vefLB2wgICHD4PCHOJSnQJ/q8A7kn+PunmfSPiqPwaC79AkLw8PJu/8RWtBUs7IbBzk1rObg7i6QpY0g3WZmx8GlMhkHexRfy3VM/Y923Bxg0fDxuFtdsPjVsVqJC/bhyWhxms0xuiN5DCvSJPi1rUzabj5jpH1X/nCQ0YkCXr9m8CN7pZeB5OTvJ2ZGJp5ud9LlzWPHO23zr58PnMdGcSJ3I6yPHEF3kz5DRU7rch7YU5u0hPjAMs1mquoq+SQKS6HW01nyRuYd9PxzleP7BJiUenKVx1u6hQwfjZq9mc9Y2LBYLv/rVr0gen8JLd91K/AVz8d1+CpPZs/2LdoHWmruuu4QBoZ0fAQrR08m4X/QqFRUVrNm4h/f/7w02rfq/Du8ZcnTza+Nl4AnjxjF//vymK/DmprFi2WqO5psJDh/g0mzddsPgh80riAx0x81N/g0p+i4JSKJXMAwDm81G5tYcXn719Q7tGTodhFa8+2f+/OxtDm1+bbwMPDExkS+++KLJCrzPM1aQMPVHuHt4dTnDQ3u8vdx57P6b8fR07ShMiHNNApLoUVrKQ6e15s9/+ycvv/0VO0+4Ya2zOrxnqHEGBm/jBG72aocCWVhMHBkrVmC1WklLSyMiIuJMVoVxCYmYPAMYO+GiLmV4cERx4VGCOUZMpCzjFn2fBCTRY7SUh27ylKn8v+XbsAaOxt0/AujYnqF2p95aCGRV5cWYOEmAgonx8fziqafI37+furo6duaXMmX+vdz601cwmc1dyvDgiLiYUFLiY51yLSF6OglIosdoKQ9dZWUl361bi9ntv+UTOrLBtL2pt8aBzG4YFBfmMXCIGzdHBrB623YW7d+PdemHzL7lZi65/UkunH83o8fPOLOAojMZHhxRXVHGiT3fsGBGPIMGDerStYToLSQgiR6j5Tx06Rw9vK/JcSazmZsff5nUy+5iR25lm3ngGo+mmk+9NQ9kNmsZlaXbsFhsFCcmUDF0MINvvA7bTbcSNmQWQ+Inn7XHqKMZHtpz7PBegswl3Ht5Ej++9XL8vKSOkTh/yMZY0WMsW7aMZ599lszMTCwWC1arlaTk8aRedteZvUEd1TwDw+efZ1BVpxmRdMGZ5eJ1dTUcy9vO5CmDMZs09tNVVSsqqDEHsn93MYEhUS2upOtIhofWKHsth7Z8xf333oXFVoZ/P18CAwM79XmF6IkkU4PodQzDYNasWZwoPMkll6R36su9JW1lYNBa4+5ehfVoFte+8f+oDQpiw6JfYGjN8neXM3riFbh7tL33pyvpgGrLC3nw2mlUlpcSHCxFlkXfJAFJ9Erbco7xuz+/RdGxvE7leuuIyvJSDu5azeVTRzHt4YX4HzhITVAQH738ItWRCRQXKSzuHi5pG+pTAdmPbuSh+26TEuOiT5PUQaLXOXDgAN9tKyJ5WrrL26qtqSQyxo0RFf7U3nAHfywrY3j/UNxf/QNfrj/ASPdxuHu4bt+P3W4nKtSPa66622VtCNHbyKIG0WPsyDnC8aIyl7dTVVHCvi2fMnDrWv6+8Be8EB5O5U9/yvNh4Tz3P68zZuI83D28XNqHwiM/YDu+1aVtCNHbyAhJ9AhVVVWUmcMJCKl1WRt2w6D4ZD4jxgaRWhJE2cNPcHTIEDJ37MBisfD8r39NUvJ4dmeva3cRxennRvkHd3cql96c6SlMGi7PjIRoTEZI4pyprKwkJycHq83O/77xDrkFx1zanmFUUFO+E4vFxqnEBL4PCWHW5Zd3eA9RZ7MzGIYNgNw9Gxjgb+Dh4brnU0L0RhKQRLerrKxkz6Ej/HvNfl774Fv+8tkuCB2Lbz/XLHWuqa7k4O5vGRhVx6TUZJRS2D09sD/yABkrO76HqDPZGXx0KTWHv+eGmUN46IaLiY4Md+pnFKIvkCk70e2+Wb+V1ZsOEz1kDDHDk9GAp5ePS9rSWtPPH3z6G1zwyOPYfH35fvEi7EpRcLQc3P1JSh7fZA9ReyXFW8vOsOPQnrOm+ux2O/ayAm69aRZm8+T6vUz9ul67SYi+SEZIotuUlZWxbuMODlQEED1kjMvbq6osY/v6D4mx53HTb35HWNZmgrbvoHbXHspKvUiccj23PvFHhzI+NNaR7AyGtZYYfxtms9mlJSqE6AtkhCS6zfb9BXz89XZiOhmM2lpI0Py92JEJRES5EV8Zwszb78XrZBFlsQNZ8z8v8tX6Awzz0Hg0pOVpXCXWESMTp7Jx9YftjqzqaqsZPzycGYnjO/V5hTjfyMZY4TSGYZCRkUF2djaJiYmkpaVhNpux2+18v3ELW497YDU69/fWuIJrWtocMjJWYvYK5ObHXwZo8t7yjBXUWOt4465bmPKL57BUVXMiKYF/3HQj3rEXYNi7vqTbkewMhbm7mJUQwfTp07rcnhC9mWRqEN3qdOmIgoICZs+ezapVq4iKimLp0qWUV1bxqz+9R9TwCZhMnZsl3pm1hsxlS9iUteGsPHfAWe+lJCXx/K5dzLPbyZ0zi++feooNG/KIHDj6rASpzlJbXYXZzYLZzY26mmquvHAEseF+MlUnznuSqUF0q8alIywWC4sWLWLixIm8/fbbBAyZSszISV26flsLCbTWZ72XdsklfHeqiCGzZ/LqgMHEHvViwJCkLn/O1vh4uqGKchg6dDhR0eF88sE7hPWLl2AkRAfIogbhFC2VjpgzZw7//ng5W/cc7vL121pIEBYVx+efL2/63ooM7D++l023/5jwgdNwc3fdKr6crFXMTQjmwTuvJ216MmMHh/L0wofx8XFNm0L0VRKQhFMkJiayatWqpkEhI4OoURfgH9S/y9dvrSjfiIQpjJ8+GXcPNyaMGcOTCxeSnJKCV4A/R49XcCRP0y8wzCUjlZqqCkYNDOSpH19PVHiI068vxPlGniEJpzj9DCk/P5/Zs2fz+fIMtMWPO5561WnZupsvJBg0fBz7t63iygvGMPmhx/n2h1w2entRefstDLvkVirK3ZtUmnUmrTX5W1byyyfuxdfX1yVtCNFXyKIG0e0Mw+CZZ59jy958ooYlu7R0RHVlGRExbgTt3cDFP3sWr6JTlMbFsvb3L/LZ+gMMHn2hyzbb2g2D8CBvrrwgFoubPIYVoj2yqEF0K7vdTmWtQXjCpcwcZ3HKyKS1fUe11eXk7FjGlJMBpP7sl7hVV3M8OZH/d8P1eOuBxI8f6YRP1Lpjh3czwi8Ui9tQl7YjxPlGApJwilVffsV3e4sJjh6BM8ZEZ+07WraEDas/YO71DzBufBTjyyOYev/DKLudw2mzyXzyKYo25mOpc8NFs3T1tObu69KICXFteQohzkeyqEF0idaa4rJKCqzhBIQPdtp1W0pgalSdImdbBmazlVPjxnB84nh23XYTv02ZxOECd6LjEly2xwigtqaKI9tXMSjMF3d317UjxPlKRkiiU+x2O0op/rNuIyvWbSN6xCSnLiBoad/RJelp5BUexWQyoU0mvvrNImwmfyIOVOPu4ee0tlvTPziAO35yO2YXPRcT4nwnIyTRYYbdzm9+/2f+unQDW4+7EzlsgtPbaGnf0cqPPyZtzTpMdXVorVn56dccybXj5x/i8g2ohUdzCXcvJjjINSUyhBASkEQH7dm7l/e+2o9HTCq1ygez2a3T6YDaMjJxKnj0Y1xCIk8++SQTRo8mJi+Py0tLqdq+k9JiL8amXo2nt/NHRnZbLda6ppVrhw0MY/SQKKe3JYT4LwlIwmFaaz5atYHco6ewuLu22mlNdQWPLl7MVTOm4vnHP/L8vn38IzyMr/7+GuvyaykqtGE2O3/GWduqSQyvxVr8w5nXSk4eIX3ycKKjo53enhDiv+QZknBIcXExOw8X4xOd6PK2amsqOLRrOVML/bjx72/iVlPD8fHJvHbt1fhYoxmZ1HZF184yDBsndn5F6ryHmDbZwidfbWXT3gKCTCVYTG2XJxdCdJ0EJOGQrG37WJV1mOi4US5rwzBsFB7JYez4KFIqIpl674Mou50f0uew4cmnKN10FA+bu8uWdft6e/LAzx/Fw1L/v8XEkf3x89BcODkdN7NMJgjhahKQRJuOHj3KnkNHOFTp79JgBKBULYo8zOb+nBo7miMXTKVkSByvRsQwMNdM5MDRLmu7qqKU2vy9eKT/d/QVFRVFVJQ8NxKiu0hAEi2qq6vDatN8s/0Y2XsLCI+Oc1lbtTVVHNr9HRfNHE1cwHCsDYsk1vzyZximAKIO1ODh5e+y9gGCAwOYM322S9sQQrRNApJo0f+99zH5lV4EhMW6NBjZrFZy939HwfavqHvnj0z39mbdX17B8HDni0/XMiwxHV//YJe1D/Wjo6EDvBk4QBYtCHEuSUASTeTl5VFQYqfCZyQBvq7d21NdWc4bv7oDXw8LaXPm8Ptde/hnfj7XbNlG7dCpjE292mXJWRtzw0qwl83l7Qgh2iYBSZyhtWbZ19kcq/IkqH+ky9qxGwbZ360gZ/tXeLuZyMrKwmKx8Pyvf01ycjLvfreLCwIm4eXj+nxxtTVVzJs+jhEDZMOrEOeaBKTznNYapRSfZ6zilBFAnW8cQS4s72M3DP7x0oMYlacIDAwgPT29aenx9HR25Fbi5eP6VEAAJw9lY0oOACQgCXGuSUA6j+3J+YF/fvAZIyakkXPEDU9v8PB0XXuGYePblf+E2lKyszezcuVKnn32WaxWKxaL5UxZ8tTL7nJdJxpRwKP33kRYoGTuFqInkM0V56mCo8f5dl8V/nFTKCytwT+oPx6e3i5t02S2cvLoTtLmzq0fDaWlERUVRWJiIgsXLjxTlnxk4lSXtF9VWcbpgpQ2q5W8rV8Q7GdxeR48IYRjJCCdhyoqKnjxT29RVm1zabmG0+pqq9m1aQVRgSVMn5BMxsoVWK1WzGYz77//PuXllazO3EnqZXdx8+MvO3UhQ11tNcVH9jMvdQD62CYCPO0YNitenu7cd9tVZ6YLhRDnnkzZnWcqKyvZnltJXPLcbhkZ2A0D/0AT44cFcOHDj3FhTQ2r+oeSnJJC2ty5fL58BUERsdzwk187fUWdYbUyNi4Yt7AKBoX58dRjP0YpxdJlX1JdayVuYJpT2xNCdI0EpPPMW+8u5Zg1hKAw1++5qautZnvmh1x9wTimP70Iv9w8qvqHct/9d7Ox2MT2PQeZPO/uM6XJnenU8Xwi3E8xZ8L1QNN9VJdfMpOamhqntieE6DoJSOeR8qpa6J9IkE27vq2SkwyI82bMqBhm3Hk/nsUllAwdzDf/81s2ZhcQNTiZoWNdlxlhxqSxTBwW0OJ7Sim8vGQhgxA9jTxDOk8UFBTwy5f+Sl03BCNrXRV5OauJ3fAlFz/4GJ7FJRybOJ7/vetO9teEMXTMhXj79HNZ+we2rGFEuBve3q5dpCGEcC4ZIZ0nSm3ehA9PdWkbhs3K0dzdJKYOIqEmhqm33YvSmkOXpbPhiYXUbCnEW7t4JZ9S3H71HML6h7i0HSGE80lA6mMMw6Cqqgo/v/9uLP3go0/IKQ/Aq59rc8K5uRl4eRVhNkdRMmokh+ZdQlV4GK8GRzDgEPSPGu7S9ivKign3qCAhXhYrCNEbSUDqY0rLKvn171/j/nvvxs/ditlkotQcgbu365Y319XWsG/rambPTiTOK4aahgUK6x56AGUOZECuDS9v103RFR45TET/YFJHRRDuY3VZO0II15KA1IecPHmS9ftKiElMY9mGfIqOHMTPy4R74CBcUO0bqM++4B8IqfEhXPjwY1gqKvl6yavU+fqyduX3xI6ciW+/IJe0bbNaGTM4lKTwWuIGhEvtIiF6OVnU0IesXb+FdVk7z+wvCo6Mwz1wkMvaq6utIfubfxFZvperfv5LQrbtwFJRScm2nZwqcid+/BUuC0YA+XsyifYqZ1pqsgQjIfoACUh9RFlFJaWWAYTHDO6e9ooLCYuwc2N8DLPuegC/3HxKhg5h9d9eY3OFN6dOGphMrv3zevju6xkxwrXPpYQQ3UcCUh9QU1PDMy+8QnFZVbe0Z62r5sjBr4nLXMXMBx/Ho6SUY5Mm8Mqdt7OvKoS4URe4dFm3zVrH0V1fExXi4/KgJ4ToPvIMqZez2WwUVRgMTErH1MkHRXbDYHf2OvIP7iY6bmSrmRMMm5WCH3aSPCWOcbZBTLrxdlYAq4cPw3L5AizewwAX1q5o4O3lyW1Xp+PmJn++QvQl8s/LXu5f7y9lyYdruhSM3vrdI2QuW8KYgb5kLlvCW797BLthnHWsxd2On18ZZrOVU4NjmRs7iKfjR+E+/3Iy3nmX5W//yeUjlvKSk0R5lTN8qOvKqgshzg0JSL2UzWajsKSScu9h+IfGdPo6u7PXYVQXsylrAy8uXsymrA3Yqk6xO3vdmWOsdbVsy/yMcP8ixkeHYTab2fL9ek7082dD9hZefPFFNm/aeNZ5rhAV7E1smGRgEKIvkoDUS639NpPfvbEUqzZ3KWt3/sHdpKXNaVK1NT19LgWH9gD103T+gTBtbBgXPfgYFzzwMO4lJRzcu5+0Oa2f15zdMNiZtYaV77/Gzqw1LY7A2lJbXUXV8b3cmJ7MmDFjOv15hRA9lwSkXsYwDI4cL+JwTQhhsaO7fL3ouJFkZKzEaq3fUHq6amtU7AisdbVkr3uXyNI9XPHUswRv34Gpzkrxzt2EhI8gY8WqFs9rriPTgq2JCvVlZtJALG7OzQouhOg5JCD1Mtt37eW3f/03ZVU2zE7Y7ToycSpmr0CSksez8Mknz1RtjRg4hP7hBjfFx3Dx3Q/gl5dP8fChrF7yGluq+hEZO6HF81qq9urItGBrrHW1HN2byYJpg5kwIaXLn1cI0XPJMqVewmazse/AYbKPmBk0eorTrmsym7n58ZfZnb2OHYf2kHrZXQwbO56DuzKYVQSpz72AubaOY6kTeW3+PHzKAxk0YgjAWee1tjqvtWnBHYf2EJ8yvcV+VVeWo4xq0i8Yh4r3xMvD9ZVthRDnloyQeolDecf467tfUFZldXqlV5PZTHzKdGYuuB3ffr5ERlu5fOQgpvz8Oait442J43lw1Gj25taglM9Z582+6h7iU6a3WmSvrWnB5uyGgdkEoyI9mDbUmzGxwYweHd8t1W2FEOeWBKQerqCggH+8t4wvdpQRN/bs6TBncne3ExBYjdnNRtmQwey68VpmDR3KK5XVxIaEsHPdZ7zz8hMdXpDQ2rRgS9N7J3I2khJpY95FSUyePMlZH00I0Qsordsu2KaUGg681+ilOOAZrfUfGh0zA/gEONTw0kda60VtXTclJUVnZWV1ps/nBZvNRllFFWu2HWfL7kME9Y90XVvWOnZsXE5a+nj8Skspiwhjy3eZfP35CsqPnyR782YsFgtWq5Wk5PGkXnZXq1NtrTm9+bbg0B6iYkc0md7TWpN/YAc/SptG8tBQ3N3dZUQkRB+ilNqktW73IXC7z5C01nuBhIaLmoECYGkLh36jtb60ox0VLduQlc3Ha3cRNTTZ5cGoX4DmwuQYZj74KJaTRcyNi6Wyqpp+vn6kp6V16NlPa05P77V0nkkpZoyLZMxAfzw8PJzyuYQQvU9Hp+xmAge01odd0RlRr6KyitzqQCLiElzajs3asKy7ZBc/+unPCd6xi1W1tVSUV5C1YSMLFy5k9erVDj376ayy4pPEBdm4dPb0JkUFhRDnn44GpGuBf7XyXqpSaqtSKkMpFd/SAUqpu5VSWUqprMLCwg42fX4oLy/nZ796haOnqlpdJNBVWmuKC48QGm7jltEDuPjun+CbX8CpkcNZekk6aWmXYLFYSEtLIyoqisTERBYuXNjms5/OCvVVxIa4rnigEKL3aPcZ0pkDlXIHjgDxWuvjzd7rB9i11hVKqXTgFa310LauJ8+Qzma1Wtl6qIS1W/Mxu7nuS9owajm4cznXuhmkLvoN5ro6jk5N5S+XXEreIl0THQAAIABJREFUSdj0xb/ZlLUBi8VCTU0Nw0eMInRgPElT57a6tLujtNZUlxzjoetn4O0pAUmIvsxpz5AaSQM2Nw9GAFrrskY/L1dK/VkpFaK1PtmB65+XCgsLKS0rxzsgjD+++hr9BozHx981Re0Mm5XD+7OZeMEw4nUsU665GaU1OQvmsemRxzHtKSc+tj9b1n1JUvJ40tPnsnz5CoLCB3HDT37llEBUevIYhq2aAQMG4WM6gYdFFnoKIep1JCBdRyvTdUqpcOC41lorpSZQPxVY5IT+9WnVtVY+WruXQ3nHiRg4nLBRF3f4S9/R0hEA7p6a0DArZotBeewgdt19O3aLhT/79ydyv5XAkAGA4xteO+JY7n4unz6SmInxVFVVMmTIEEBy0gkh/suhKTtVvxsyF4jTWpc2vHYvgNb6NaXUj4H7ABtQDTyqtf6urWue71N227Zv54MvswkelNTusa0FndM54ozqYtLS5pCRsRKzVyA3P/5ykwBis1rZvv4z0i+ZQL/iU1QMrA885SWlKHMwhcdNuHu4NoN2pFcFM5IHExoa6tJ2hBA9j6NTdg4/Q3K28zkgaa1ZtTGXLfuP4unddkG7toLO7ux1ZC5bcuZ5T0v7hKx1tQQEgeepnaQ/vxivE4X85++vUREawpt/+Du1NWYGj0p2eBTUkRHZ6eON8iM8cvMczGaZnhPifORoQJJviG5WXl7Oc7/9E7vyStsNRtB2YtL2SkfYrHVs+fY9wot2cPljTxG0czd2i4Wje/by4iNPk79rGynD+zucfbtTWbuNOoaEaEwm2egqhGibBKRucuzYMVav/Z7teVW4hyc4XFm1raDTWo64yEHDKTqeR/9wK7fFxzDr3gfxLTjCqZHDWb3kNZbvPoG1sqzD2bc7krXbbhhkff0ZOd9/iJuyY7fbO3jHhBDnGwlILpaTk8Oq7/fw0Xd5rN52ko17C/HpF+jw+W0lJm0tR9yIhAmUFG4mdu1yZjz4OB6lZRyZOpmXb76ZXSd9qKmykZ421+Hieqe1NyI77fRIastX7+LnCc8++ywLFizA6GAOPCHE+UUCkout35FL5q4jWJUXYVGxrR7XWkXVthKTni4dkXrZXezIrWR82q1EDh7J96v+SnB1CZOeeR5zXR05P5rH2hcWYwlNxt0zqEPZtxtz9LwdG7/GqC5m86YsFi9eTGZmJvn5+WRkZHTlVgoh+jgJSC5itVr5LmsHpwjDP6h/m8e29WymedBJveyuJqvo/ls64g62fPMJudu/YUhIKB9lfMHsYUPJ/vE9vByfyL59tfgHRWF2s3Qo+3Zjjp53eO8m0puNpObMmcOWLVu6cEeFEH2drLJzkcP5R/nDm8sZNGpiu8fuzFrT7mq5thg2K5+/8xInf9hJ1sas/14jJZmZ19xMVNx03D28mpzTVvbttrSXtbuuupwI01H+8LsXyMzMPNOXiRMnsmjRIi69VPLvCnG+cUWmBuGg/Px8Nh6qcygYQccqqjZfdj04PpngUAuBHrWMn930GnNnz2HVZyuJGvTDWUu028q+3ZaWzjMMG4bNSnSoL/mHs7nqvjv497/eZOLEicyZM4eVK1cSHR1NWlpah9oSQpxfZMrOyWpqbfzvm5+x//Axh89x9NlMS1N7ry+6ndATW7n2y6/44pNPzlyjpqaGN998C6PiuONLtDtBAebyHxgTWsV1M0fy00cewNPTk6VLl7Jo0SJ8fHxYtGgRS5cuxeyiZLFCiL5BpuycpK6ujrff/RAjeAxVdR0793SgsVWdOpM/zs076KyMCy1N7aUkJfHLAweYV13NpaEhFEREkDY3nXffew9vLy+2bdvaLEnqKJKmpjklHVBNVQUDA01ckzZeCuoJIVolU3bdSGvN3oIKDlX6E+yrO/zlfHrhQnv541qa2ktLS2PnSy8xcdoUkubOpf9J2JF3nNCB8cxMHY3FYsEwDK6++mr8fL2ZOSmejGVL2Lj6w7MCXke5U0eMv0WCkRDCKWSE1EVaa155dQk1AfF4+fRzaVvbN/yHzGVLyN686b+LBeLjuSc0GP8XXyJnfzXeviGY3SxNRlMrV67k2WefbbLIoLOlyE+rLCvmtkuTiAj2cfKnFEL0NZI6qJvUWe2Y+o/B09u11U611oybNBmTpzvJKSk8+eSTjE9MIMCk2HrRHPbsrsIvIOJMHaXGS7RffPFFZs6c2eGNsG31pfbYNrzMVqd9PiGEkBFSFxw6dIiPv9qGKXCwS9sxDBtbvl1KevoEAo4dY83Jk/ywL4fQiAjGTZpFSYk3bhaPs847vSJv87oMak7lNhlZdXaEVHrqBFH9g7glfQxuskhBCOEAGSG5wMZN2RQUFAD1o4RTtW6U2dtPkNoVtdVV9PM3SE8dyuyHHuOiBx5mRnQUV9xxC8oUQH5uXYvBCP67RPuGn/wai09IhzfCnmY3DIqO5xMV7M2IUDszRvWTYCSEcDpZ1OCgWqtBxqbj+PWrZnB0FVlrlhE68kICgsNc1qZhWNm58SOuTB7GRU/8HJ+jx6iIiiTvcC5VbgOIi5/t0HUcXTTRmkAfMwHVxfxo6iBMpriufCQhhGiVBCQHVFdXs/iVJQQMnY7JZOJQYS2+sVMwMOOK9WVaa04cOcSIMf25ZUQUU+9/CPfyCk7Fj2Td737Lwf1l+FWYcaB6xRmd2QhbWV5KzfFd3P/I7VjcRnfikwghhOMkIDnAzeJOv5jEJiUjPLxcWWHVSm3lHgZ/tZ5JL7yE2WqlYPpUXktLw+u4O2HRY13Ydj1bXS1zJg2lv9cALG4yPSeEcD0JSO2oq6vj0y/X4+bteMmIzrIbBvt3fEvqRaMY6j6Y1AXXYrLZ2HvlAjY/+AheB614+gS5vB+V5SWok9tIvibZ5W0JIcRpEpDaUV5eQdauw0QMdm5Aap6TbkTCFDy9NLFDfXBzs7Nufw6fTU0l3r8fe0aNI2R3JUGh0U7tQ4v9stuZmjiYScMnuLwtIYRoTFbZtcEwDI6W2YkYnODU67aUk+7Vp28gwJTLuLpqXvnZL1n+1r8ImDiZN37I59sVXxMQFOHUPrTm4KZVDAs14e7u3i3tCSHEaRKQ2rBv3z6WvPW+06/bUilwH3czpnvup+q+h6gpLiFrw0YWL17Mpo0bMaqL2y0v3lVaa7w9zPzi0TsIC2u7fpMQQriCBKQWbNu+kw8//5qdJz0ZMCrV6ddvsdzE3Lnk5B9hk4cHc+emOS2rQntqa6oAqDh+gGAjj7CQQMlNJ4Q4JyQgNdBa886//s3y7w+wclspe06aKThZeSYVjzNFxY7g8+UZTcpNrPr4Y+IGDkA/9jCfZ6zocHnxjrDb7dTV1hAe6EnRri+56aJYHrl5FpfM7lxeOyGEcIbzflFDaWkp+w78QLkplD0lPgTml+Hu3Q9XPkEZlTiJdZ+/QUpiImnp6az6+GNCqyrZ+cB9uIcm4ua9jqTk8U1KUTiaVcER9tLDDAoyc9WMFK6e/pSMiIQQPcJ5nctOa83aLYf4dNX3DBg2zuXt2erqWP3x63h51TEyOozk37zETptByJRU/H71Arm5Gg/vQBSqU+XF25O7bytXzxrH+IRRKKUkEAkhuoWjuezO24C0dfsOPlm9Gf+BSd3SnmGz8fbLj2CrPEV6WhoZK1cQYBj8/qLpvBUciX9oIsFhMS7tw5AQmDZuAP36ubZMhhBCNCYF+lpRWFhIZa2dzQUmPMOc91ymLXbD4NO3foNRdYpNWVlYLBaef/55klKSeStkAAOHX4zZzbXLrEtP5DJt1gz6+boyw4QQQnTeebGowW63c/DgQY4XV/PBF9m8/uF3lFXb8fB0/ZdzdWU5/fxthPnYSZvdfGVdGtuzdrs8GGmtCXIrw2KWKTohRM91XgSk8qo6/vetz3n7i93UeUUSNmBot7Rrtxvs3bKM/kc3c/XKL/jik0+arp7LWMmQ+HZHsV3m4w733XY1Xl5eLm9LCCE6q89P2RUVFbFm+3EGjbuw09donuanvUUGWmuO5e1jZEIkt8SFMvmBRzFXVLDE4k5ySjJpc9P5PCMDi3ewU1fPtaSqsozS/A2YLpFs3UKInq3PB6Rtu3PI3JzX6X08p9P8GNXFpKXNIWPZEjau/pCbH3+51aCklAH2XIZ8uZHUxf+DyWYj/8LpJM6aReARKzvyiph82d1OWz3XmtqaKobEhPGj634iK+qEED1enw5IhmFQbg4nKtan09donObHYrHwq+efJyl5PLuz151VW8hut7N3y9ekzhxNnPdgJi74OSabjd3XXEn2Aw/R77CdlKH9MJu74bZrTcmB75l20XWYpbqrEKIX6NMB6S9vvEmR28AuVXVtMc1P+lx2HNrTJCDZ7XY8Pe2MGBuKl7cbtb5BrH/+GXyOHuMfQeH47SglJHxglz9Ta2qqKlAmE+4eXuTv3chDd1xF1OU/aVLDSQgherI+/W0VNDgV/6CuJQqNjhtJRsbKNlP52O12Nq15nwDTD4w9cYzsb7/ngyX/4J/l5ayfeTkDhs1xWTDSWjMozJdoj1MkRsKNM4dyxYXxRAZ7SzASQvQqfXKEZLfb+ce/PqLEfXCXc9GNTJzKxtUftprKp7K8lLBIDy6fEc+FDz7OrUePkRcXx9z0dJZ9+BlrPlvDHU+96oyP1aLDezaREj6S+Vemn3ktZNJEl7UnhBCu0qcyNVRVVVFw9DjFVh8+WfkNoTHDnPIw//Qqu+apfOx2gx3rP+TykdHMfOpp/nP8BM8MH8767duxWCxYrVaSkseTetldZz1vcpZJw4MYGxsoS7qFED3WeZepwbDbWZt9kFVrsxg4IoX+A4Y77doms5n4lOlngorWmiM/7GFUYiQ3Dwxi8k8exVJZybdh/Zl9+eXtPm9yBq01R3M2MXLmFRKMhBB9Qp94yHDgwAFeeu199p6AgSNcv9FUKQOz6QhDv/iYaQ89jqWykvyZMzg4d+5ZZSWcXTqisRkpw/DxkVRAQoi+oU8EpCMVHhjekS5vx263szPrC/z9TzJpVCzjX3kVk2Gw69qrWPPMIhJm3o7ZO4ik5PEsfPJJkpLHO710BIBhs9JPlZF20WRZ0i2E6DN69TOkyspKXvt/72IKT3HpBlOof47k6aWx1uwlPCYCk8lEWOYGfHPz+L+gSLz9RhAaGdvq8yanqqsg2qOQKy6/zLnXFUIIF+jzz5C01uQXllHlNYh+rg5Gdjub1v6b9NlJxBTkc3xgFAC7YgehhqQQU9EPU8Nm1+bPm5zdj7ITP/DgjbMJ9PN0+vWFEOJc6rVTdlarlT+/9jf6BXV+06sjKsqK8etXy49mxJP+6BNMeXQh/ddvBODQ/iKO5NaeCUau5m7WjOiv6eft/LLqQghxrvXagOTu7s6VN9zp0ja0tvPDntX0z9vE3PseJmD/ASqio9hZVc3BfZVEDEwlMDTCZe2XHt1PbelR+nlbyNubxYIpsVxzxXx5biSE6JN67ZSdYRjs3LoJfJyfAUFrTcGhXcQnRXNjtD+TH6xfSVc4bizfvvQbCnMN3O2uq7paWV5K4rBwElNTcXd3JygoiC1h1QT383BZm0IIca712oCktaamugrV+byprTKZ7Hh4nGTYl5uY9OLvMRkGebMu4q8Xz8KUYyU4fLDzG23EXH2MaB9/whu1k5CQ4NI2hRDiXOu1AcnNzY3kSdPYnHOyzeM6UstIa82urC+ZNHM00cMGcPyJn/JrwyBo+lT8nvwF/Y65Y3H3dcXHAcBmrcPPYuP+OxZgceu1s6lCCNEpvfZbz2q1suyjf7Z5zOlaRpnLljBmoC+Zy5bw1u8ewW4YZx1rGDY8PA3GTojB29uN377wEs8NGULlE0/wl6JSXnzsZ1jcfV27vLymhGjPYglGQojzUq/95jObzaRecHGbxzSuZfTi4sVsytqAreoUu7PXNTlOa032ug/wNw6QeCCHrZkbqC4tY+OmTSx+8UWyN2/CZFSfOc9uGOzMWsPK919jZ9aaFgNcR9gNg+qSo9x91TTS0+Z06VpCCNFb9dqApJRqd7VZa7WMCg7tOXNMeUkRvn41zJsyCtM9D7D6pz9j0zvvkTZnbovndWTU5Sg3ZSXWrwJfz147gyqEEF3WawOS1WplzZeft3lMe7WMtLaTf2AtwYcyefOu+3nJz5+qxx9nX2ERn3zySYvnOTrqcoTWmsL8fVx3cTzXXLlAyowLIc5rvTYgubu7M/+qm896vfF0mt1ux+QVeFZuuREJU8jN2U5AQBk3hPugHvwpR8PCyNyxg8W//S179u6lqOgUCYlJZ+Wkaz7qMplMDBs2hK8/+78OT9+5u8G4aHf8vWVkJIQQvTYg2Ww2tmR93+S15tNpG5b/HbRm4iV3sCO3ktTL7uLmx1/GzaLw9S1l+KqlTHn4CbbX1jKrUdkIT09Pbrn1FryCBzY5z2Q2Nxl1GYbB/Pnz2bVzJwvSpjs8fae1Jm/fZi6dEMVVCy49064QQpzPem1AUkqdNcXV0nSaUV2MyWRi9lX3MCJhChnv/ZGvPv49RvExxv71b5gMg6Dp08hY2XRqb8WKVSRNTWP2VfcQnzL9zOq6kYlTMTeMuq6+5hoOHDjAtm3bePHFFx2evvP2cCM9dSjhQa5bQi6EEL1Nrw1IZrOZccmTmrzW1iIGa10Nb7/8CMf2ZTE8PILP/v0xc4YMZsPCRzmYNh+r8nKobITJbObmx18m9bK72L4vn3nz5rW5aKK5H3Z+x/SRvlx0QX0WBiGEEPV6bfmJuro6nlr0MgOSLznz2s6sNWQuW8KmrA1NSohPuvROcnZ8w6nDe8jKymr0XjLTLr+GYQmXoe26w2UjWmuvecny6rKT1BUd4M5brsWoOkVERARubvLcSAhxfujz5ScsFgsXzZ1HTuF/n9eMTJzKxtUfkpQ8nvT0uSxfvgLc/Rh/QQqVh9YycU6z0VNaOhu27md4ogllpsNlI1pqr/HIyma14mOu4/YrJ2KrG0dgoA8EuyDXkRBC9AG9doRkt9v55Ktt5JY1jamNC+RFDhqG2VRC+oBArD/9OS+Eh5O5Y0ebo5kO96ONgnxG5UmGBFQz79K0Tl9fCCF6uz4/QjIMgw3ffU346KbZGkxmM6OSL8DXP4DRydEM3ZXFpEeewFRVzRteniSnpJCWlsbnyzOweAd3ubx4SwX57HY7tWUnuOeKKXy7djXPP/88iYmJpKWlSekIIYRoRa8dIQF8s/1Yi8lVlTIoPrGB2Qd2M/F//ojJMMidM4u/zpzJvpwSqsorXFdeHMBWQ7+a/Xz0/j85cuQIs2fPZtWqVURFRbF06VIJSkKI80qfHyHZbDayMtdiChl15jWtNTs2rCT1otEMGD6QhF/8DJNhsP3Ga8m+8z4Cj7kzOdbLpQlSTx45yH3XXMTmzGMcOXKEzMxMLBYLixYtYuLEiWRkZHDppZe6rH0hhOitem1AUkrh7e1LTcPvNqsVb19FytQh+AV4YzOZWPeHlwjctZt/B0Vh2lpEeMxQl/bJhGZokI1AHzPZ2dnMnj27ySKKOXPmsGXLFglIQgjRgl69D2nU2CSgfmS09ful+NXuJWnzJkym+o+V4+dL5oRZRAya4fJgdCx3PxeODeHm636El5cXiYmJrFq1qslm25UrV0qhPSGEaEWvHSHV1dXx77ffIHjoZKIG+nHFBaOY8fhT+B88hFaKvLTZHM2vAZNBcH/XxV273Y6Xh4WkWD+ig73OvJ6Wlsbrr7/OxIkTmTNnDitXriQ6Opq0NFlxJ4QQLem1AclisZC+4FqWrfyQBJsPs59+Dq+TRZTGDiLLy4vibUUEh7t2NKK15mDWcn72yF1E9B/Z5D2z2czSpUvJyMhgy5YtLFq0SFbZCSFEG3rtKrva2lqeeeklLvPzZNJTz+BWXc2JpATWvfgCe4+YMZn74entmlxxWmvKCvO44bIphPUz4+srOemEEKI1fX6VHcCkzVlM/nQZJsPg8NxZ/G3WHKw7ywgfEO/ahu0GUV7lDAj1lkzdQgjhJO0GJKXUcOC9Ri/FAc9orf/Q6BgFvAKkA1XArVrrzU7uaxPuVitzN2TVL+u+6To233Y3/oXemEweDl/jdJaF/IO7iY4b6dC+pMryEi5IGMjUsdd09SMIIYRopN2n/VrrvVrrBK11ApBMfcBZ2uywNGBow393A39xdkebqzaZWJiUQNbTT/LPpCns3nYKN4u3w3uMOl2KvKKAQHOZEz6BEEKIxjo6ZTcTOKC1Ptzs9fnAW7r+gVSmUipAKRWhtT7qlF62wMPDg6FXXcV3vhGE6XCU6thKusa1kywWC796/nmSksezO3tdi7nt7IaBhToeueMKvDx69UynEEL0SB1dD30t8K8WXo8C8hr9nt/wmsuYzWb8vCM5eUJ3OBhB27WTGjt5LA+zthLjb8erYq8EIyGEcBGHv8mVUu7APODfnW1MKXW3UipLKZVVWFjY2csA9fuQNnz7DcH9Oxf3Gpcih/qNq8uXryAqdsSZY0IDPBkbCVdOjuDKWYncceuNXeqzEEKI1nVkaJEGbNZaH2/hvQIgptHv0Q2vNaG1fl1rnaK1TgkNDe1YT5uxWCxceeOdnT6/cSnylqrE5u7NZnSYnavmpxEeHt6lvgohhGhfR+afrqPl6TqAT4EfK6XeBSYCpa58fgT1e4GOFeQB/p06/3Qp8t3Z69hxaA+pl93VZJXdZRcmMXigS2cdhRBCNOJQQFJK+QCzgHsavXYvgNb6NWA59Uu+c6hfhXeb03vajGEY5OzbhX9caqev0VIto9rqKtyqC5ixYIEzuimEEMJBDgUkrXUlENzstdca/ayBB5zbtbZZLBYunH1Zi/WQusLf10LSsAFOvaYQQoj29cps34Zh8Mknn/DcwvvYmbWm/b1DDrBZreTt+p4F04YxPiXZCb0UQgjREb1uDbNhGCxYsICCggJmzZrF8mVL2Lj6Q25+/OVOF96zWeuIjw3lomETCPH3av8EIYQQTtfrRkgZGRkUFBSQmZnJ4sWL2ZS1AVvVKXZnr+vU9epqKinLWcOclGgSE8adqaUkhBCie/W6b9+WKrG2tKHVEYbNysXjB/P0Ez+hPh2fEEKIc6XXBaSWKrE239DqqB+2fU1/r1rc3HrdzKUQQvQ5vS4gpaWlERUVxcSJE3myhQ2tjnJ3M/H0I3cSGRnpop4KIYToiF43NDhdiXXZsmX86fU3z9rQ6ojSUyfwtxUQGjTKhT0VQgjREb0uIEF9UJo/fz5BcRM7tA9Ja01NVQWjhw1k2ogxLuyhEEKIjup1U3anWa1WVq/4xKFja6oqsNXVEupjx1S4hcunxNHVXHpCCCGcq1eOkKB+lDR81FhO1LX8fm11FcWF+aROSKZg7x5Sx44kYUw8evZYWVEnhBA9UK8NSEop+odHciK3/Oz3gIsSI6ksrGV66kBIHdjkPCGEED1Pr56y+/Cffz/r9braao7sWE3isAimXzDtHPRMCCFEZ/TagOTu7s51t9531uvhIf48eOc1mDuZRkgIIcS50WsDkmEY7N21rclreTk7iO+viYyQgnpCCNHb9NqApLWmuKi+DLrdMNix8WuO7/2GfbuyMZyQ/VsIIUT36rUByc3NjYlTL8JuGPz9hXtY//kbRIX48dJLL7FgwQIJSkII0cv02oBkGAaf/usv9Kvchbuysikri8WLF5OZmUl+fj4ZGRnnuotCCCE6oNcGJLPZzAvP/YIjuTlnZf+eM2cOW7ZsOcc9FEII0RG9NiBB/Uq7lrJ/r1y5koSEhHPcOyGEEB3RqwMSNM3+/dRTTzFx4kSio6NJS0s7110TQgjRAUprfU4aTklJ0VlZWU65lmEYZGRksGXLFhISEkhLS5N9SEII0UMopTZprVPaPa4vBCQhhBA9l6MBqddP2QkhhOgbJCAJIYToESQgCSGE6BEkIAkhhOgRJCAJIYToESQgCSGE6BEkIAkhhOgRJCAJIYToESQgCSGE6BEkIAkhhOgRJCAJIYToESQgCSGE6BEkIAkhhOgRzlm2b6VUIXC4i5cJAU46oTt9idyTs8k9aZncl7PJPTmbM+7JQK11aHsHnbOA5AxKqSxHUpqfT+SenE3uScvkvpxN7snZuvOeyJSdEEKIHkECkhBCiB6htwek1891B3oguSdnk3vSMrkvZ5N7crZuuye9+hmSEEKIvqO3j5CEEEL0ET0+ICmlYpRSXymldimldiqlHmrhGKWU+qNSKkcptU0plXQu+tpdHLwnNzTci+1Kqe+UUuPORV+7iyP3pNGx45VSNqXUld3Zx+7m6D1RSs1QSm1pOGZNd/ezOzn4/46/UuozpdTWhmNuOxd97S5KKU+l1IZGn/e5Fo7xUEq91/Adu14pNcglndFa9+j/gAggqeFnP2AfMKrZMelABqCAScD6c93vHnBPJgOBDT+nyT05c5wZ+A+wHLjyXPf7XN8TIADYBQxo+L3/ue53D7gnPwNebPg5FDgFuJ/rvrvwnijAt+FnC7AemNTsmPuB1xp+vhZ4zxV96fEjJK31Ua315oafy4HdQFSzw+YDb+l6mUCAUiqim7vabRy5J1rr77TWxQ2/ZgLR3dvL7uXg3wnAT4APgRPd2L1zwsF7cj3wkdY6t+G4Pn1fHLwnGvBTSinAl/qAZOvWjnajhu/NioZfLQ3/NV9cMB94s+HnD4CZDffHqXp8QGqsYZiYSH0EbywKyGv0ez4tfxn1OW3ck8buoH4EeV5o7Z4opaKABcBfur9X51YbfyfDgECl1NdKqU1KqZu7u2/nShv35E/ASOAIsB14SGtt79bOdTOllFkptYX6f6h9obVu9TtWa20DSoFgZ/fDzdkXdBWllC/1/7J9WGtddq770xM4ck+UUhdSH5CmdmffzpV27skfgIVaa7sL/nHXY7VzT9yAZGAm4AV8r5TK1Frv6+Zudqt27skcYAtwETAY+EIp9U1f/t7RWhtAglJ8p5aRAAAB2UlEQVQqAFiqlBqttd7R3f3oFSMkpZSF+j+ed7TWH7VwSAEQ0+j36IbX+iwH7glKqbHAEmC+1rqoO/t3LjhwT1KAd5VSPwBXAn9WSl3ejV3sdg7ck3xgpda6Umt9ElgL9PUFMO3dk9uon8bUWusc4BAwojv7eK5orUuAr4C5zd468x2rlHID/AGnf6f0+IDUME/5N2C31vr3rRz2KXBzw2q7SUCp1vpot3WymzlyT5RSA4CPgJv6+r92wbF7orWO1VoP0loPon4e/H6t9cfd2M1u5eD/O58AU5VSbkopb2Ai9c9V+iQH70ku9SNGlFJhwHDgYPf0sPsppUIbRkYopbyAWcCeZod9CtzS8POVwH90wwoHZ+oNU3ZTgJuA7Q1znFC/CmYAgNb6NepXTKUDOUAV9f/C6cscuSfPUD/H++eG6Smb7ttJIx25J+ebdu+J1nq3UmoFsA2wA0vOxVRNN3Lk7+R54B9Kqe3Ur0Bb2DB67KsigDeVUmbqBynva62XKaUWAVla60+pD+L/p5TKoX6Rx7Wu6IhkahBCCNEj9PgpOyGEEOcHCUhCCCF6BAlIQgghegQJSEIIIXoECUhCCCF6BAlIQgghegQJSEIIIXoECUhCCCF6hP8PigImVqt80vAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# generate prediction intervals\n", "def pred_modelfun(preddata, theta):\n", " return linear_model(preddata.xdata[0], theta)\n", " \n", "mcstat.PI.setup_prediction_interval_calculation(\n", " results=results,\n", " data=mcstat.data,\n", " modelfunction=pred_modelfun,\n", " burnin=burnin)\n", "mcstat.PI.generate_prediction_intervals(calc_pred_int=True)\n", "# plot prediction intervals\n", "data_display = dict(\n", " marker='o',\n", " markersize=5,\n", " markeredgecolor='k',\n", " markeredgewidth=1,\n", " markerfacecolor='w')\n", "model_display = dict(color='r', linestyle='--')\n", "mcstat.PI.plot_prediction_intervals(\n", " adddata=True,\n", " plot_pred_int=True,\n", " figsizeinches=(6, 6),\n", " model_display=model_display,\n", " data_display=data_display);" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": true, "user_envs_cfg": false }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }