{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\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": [ "