{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "# Specifying Sample Variables\n", "Author(s): Paul Miles | Date Created: July 19, 2019\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": [ { "name": "stdout", "output_type": "stream", "text": [ "1.9.0\n" ] } ], "source": [ "# Import required packages\n", "import numpy as np\n", "from pymcmcstat.MCMC import MCMC\n", "import matplotlib.pyplot as plt\n", "import pymcmcstat\n", "print(pymcmcstat.__version__)\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": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hUZfr/8fdNAgiigtgoZsECIsWyUYi6Cou66teuX8VVcW3ID4KAiyKsIIJflCJGwAa6ClhXRcWCy9LEXQcUsCBFZS00C6CgUtLm+f1xJjoZUibJydTP67q4Mpk5zDznQj95cp/7eY455xARkeRXJ94DEBERfyjQRURShAJdRCRFKNBFRFKEAl1EJEVkxuuDDzjgANeqVat4fbyISFJatmzZFufcgWW9FrdAb9WqFUuXLo3Xx4uIJCUz+7q811RyERFJEQp0EZEUoUAXEUkRCnQRkRShQBcRSREKdBGRFBFVoJtZfzP7xMxWmtmAMl43M5toZmvN7GMzO97/oYqIJL9AAO65x/vqt0r70M2sA3AjcCJQALxlZq8759aGHXY2cGToT2fg4dBXEREJCQSge3coKIB69WDePMjJ8e/9o5mhtwOWOOd2OueKgLeBiyOOuQCY7jyLgcZm1sy/YYqIJL+F84NcvXsq7Ys/oqAAFi709/2jCfRPgD+YWVMzawicAxwacUwLYH3Y9xtCz5ViZr3MbKmZLd28eXN1xywi4pvaLIGEf8aTfd/n+se68KjrxXX2JPXqQdeu/n5OpSUX59xqMxsDzAF2AB8CxdX5MOfcFGAKQHZ2tm6VJCJxVdslEID339rKqnOHcm3xVL7nYOZf9xQ7D/8z87r5/1lRXRR1zj3unPu9c+5U4Efgs4hDNlJ61t4y9JyISMJauNAL8+Ji/C+BFBfDo4/S4eI2XFP8OHkMoF2dT1lyxJUMGWq+hzlE3+VyUOhrFl79/JmIQ2YBPUPdLl2A7c65b3wdqYiIz7p29WbmGRn4WwJ57z3o3Bl696awbQe61P+Q2zImkF9/X9/LLOGi3W3xJTNrChQCfZ1z28ysN4Bz7hHgTbza+lpgJ3BtbQxWRMRPOTlemWXhQi/Mazxr3rIFhgyBxx+HQw6Bp59m3yuuYNJi8+8zKmDOxaeUnZ2d7bR9roikhOJimDoVhg6Fn3+G/v1h+HDYd1/fP8rMljnnsst6LW77oYuIpIQlS6BPH1i+HLp1g0mToH37uAxFS/9FRKpj82a44Qbo0gW+/Raee86r38QpzEGBLiJSqVK96sXF8OCD0KYNTJsGt94Ka9bA5ZeDWVzHqZKLiEgFwnvVT8kI8HqrvjT67AP44x9h8mRo1y7eQ/yVZugiIhVYuBD2y/+eKcXXsbDgJILffg/PPw9z5yZUmIMCXUTkV3tsA1BURI8tk1kVbMtVPMX4zMGsnrkGLrus0vJKLLYUiKSSi4gIe24DsCTvXTo+3JfWH37IthNO54mTJ3HyZUfROYo+8lhsKVAWBbqICL9tA9C0+DvG7h5Mx5umQcuW8MILNL7kEm6qwgXPsrYUiEWgq+QiIgJ0PaWI/nUm8iltucI9w8arb/e6Vy69tMrdK7W2pUAlNEMXEXnnHXJyc8kp/JgvjjyTL0dO5Lgebav9dr5vKRAlBbqIpK9vv4XbboMZMyArC156icMuusiXfvKcnNgFeQmVXEQk/RQVQV4etG3rtSAOHQqrV8PFF8d9cVBNaIYuIikrECij7LFoEfTtC598An/6k7f3ypFHVus9IfZllYoo0EUkLsoMW5/fP7x1cNHz35D9/K3w9NNeeWXmTLjwwirNyMPfMyPD+6tFRbFtTayIAl1EYi4WfdolrYNWXEif3ZPo+L8jwOXDHXd4e5Y3bFjt9ywuhmDQe8652LYmVkSBLiIxF4s+7a5doXvGQu4rzqWDW8m6I8/mzTMf4JhzjoSPqvfbQUk7Ylkz9Fi1JlZEgS4iMRcejLUShhs3kjNxEP8seI5tTVrxzx6vcNET51Ow2siYXP1SSWQ7IpT+wVDbZaTKKNBFJOZqrU+7sBAeeADuust7PHw4jW+/neV5DSgo9KdUEtmOWPI4Xsv9wynQRSQufO/TXrAAcnNh1Sr4n//xgv3ww4HYlEritdw/nAJdRJLbxo0waJB3x6DWrWHWLDjvvFKHVFYq8UOtl5GioJtEi0hyKij4rbxSVAS33w6DB0ODBr68fXXq4bGooesm0SKSWubP98orq1d7s/G8PDjsMN/evrr18Hgs9w+npf8ikjw2bPDu3dm9O7u35/NCz9cIDJnla5hD2fXwZKBAF5HEV1AAY8bAUUfBrFmsv+Eumv2wkiuePpfu3f2/K1C8tr+tKQW6iCS2f/0LOnXyauSnnw6rVvHUYcP5uXCvWptBl1xEHTUqMZb0R0s1dBFJTOvWwS23wEsv8eP+h/Pt+Ddo99dzgNh0lMS7Hl4dmqGLSGLJz/furtyuHcHX3mBk5khabvuE3w8759fSSrLOoGubZugiEjOVtvXNmcOuG/vRYN1n/HDqhTzX+X5GTmhFcRAyCmD69NJ/P3yVZiJtYxsvCnQRiYkKWwHXrYOBA2HmTDbaEfSvM5sF759F3pWlV3g+8cSee7AkwpL7RKGSi4jERJmtgPn5MHq0170yezYLz/w/OtknvBk8i4IC2Lr1t9LKddd5YR55ITRZWwxrg2boIhITkRcyL9jrn9CxH3z+OVxyCUyYQP2NWfCOV14pudhZUloJBGDatD0vhCbCkvtEoaX/IhIzgQB88PJX9FgykP0XvQJt2ni3gDvzzFLHlFcPL++1dKqhV7T0X4EuIpXyJTB374Zx47wSS506MGyYVzevXz+tArmmtJeLiFSbLxcd33gD+veH//4XLr0UJkyAQw/17/0F0EVREalEjS46fvklP/zhAjj3XHYVZnqrPl944dcwr/H7SykKdBGpULX2Ndm9G0aOJHjU0dT79zxutzEc8v3HBPY+3Z/3lzKp5CIiFary7eJef90rr3zxBWs6XsbZK+9jXbAlGYVl38Wn1m5Hl4aiCnQzGwjcADhgBXCtc2532Ot/AcYBG0NPTXbOPebvUEXED9W5ABnVviZffOEF+euve33lc+eyvWF3Nncv3YZY7feXSlUa6GbWArgZONo5t8vM/gH0AJ6MOPR551yu/0MUEb/UygXIXbu8rW3vvRcyM2HsWC/Y69UjB82+Yynakksm0MDMCoGGwKbaG5KI1BZfb2TsHLz2GgwYAF9+CT16wPjx0KJFqcNiPftO5xbISgPdObfRzMYD64BdwBzn3JwyDr3EzE4FPgMGOufWRx5gZr2AXgBZWVk1GriIVJ1vqyr/+1+4+WZ48004+mjvlnDduvk40upJ9xbISrtczKwJcAHQGmgO7G1mV0Uc9hrQyjnXCfgXMK2s93LOTXHOZTvnsg888MCajVxEqqzG287u3AnDh0P79rBokTcj//DDhAhzUAtkNCWX04EvnXObAcxsJnAS8FTJAc65rWHHPwaM9XOQIuKfapVAnINXX/XKK19/zcpOV7Br1Hiyz29eK2OsrnTf1yWaQF8HdDGzhngll+5AqTX7ZtbMOfdN6NvzgdW+jlJE4mftWq+8Mns2O1u356J6C5i3siv1ekBenrcjYqLUq9O9BTKaGvoSM3sRWA4UAR8AU8xsJLDUOTcLuNnMzg+9/gPwl9obsojExM6d3r4r48ZB/fpw331M2tmPeSPqUlzs7XybmwvBYGLVq9O5BTKqLhfn3J3AnRFPDw97fQgwxMdxiUi8OAevvOKVV9atgyuv9EK9WTNODUC90V5Jw8yrVQeDe3bMpHOnSTxppahImisVvk0/g379YM4c6NgR3n4bTj3112PDSxpNm3qZH1mvTvdOk3hSoIuksZLwzczfQWad/+NEN57CzAZ80z+P1uP7eguFIoSXNDp2/O2HAXj3dl63zsded6kSBbpIGlu4wHFO/stMCA4gK7iep+xqBheO5ccphzDv8sqDOPxuQiWz8oyM334OpGOnSTwp0EXS1Wef8f9e7Ufj4Bw+ohM9M57hHXcKwaC390pVZtbh/d8AN94IWVmqoceaAl0k3ezYAXffDffdR+MGDfhywAPMbtqHPx+UyXtl1MSjEdn/3bOngjweFOgiKaTC7hLn4MUX4ZZbYMMGuOYaGDOG1gcfzO2hQ8Jr4lUJ5HTv/04UCnSRJFcS4pFdJ6W6S9as8bpX5s6FY4+F556Dk0/e471q0sOdzv3fiUKBLpLEwi9Gmnk94aX6wjv+AqNGEZxwPwWZDdn018kcNqa3d+VSUo5uQSeSxMIvRgaDXk5nZEC9uo5Liv/h3Whi7FieclfRuuAzOjzUl8B70YV5IOC1IQYCtXsO4h/N0EWSUHiZJfxiZF4esHo1l7+Ty37D5sOxxzL9vBe4bmoOxVXoXtHioOSkQBdJMpFhW7JB1h9P+JnOc0bB5PuhUSN48EG46SaOfC+DetOq1r3i640wJGYU6CJJJjJst25xDDnsebjmr7BpE1x3nXc7uNA9B6rTgZLu29AmKwW6SJIJD9tjMlfS56VcWL4Qjj8eXnoJunTZ4+9UtQNFbYjJSYEukmRycmDBrJ/hrrs4IfAAdb7cBx56CHr18rV7RW2IyUeBLpJMnINnn6XzoEHw7bdw/fVeK8oBB8R7ZJIAFOgiyeKTT7w7Srz9NmRne3uWn3hivEclCUSBLhIn4cv0oYJ69fbtMGIETJoE++0Hjz7qzcy1OEgiKNBF4iByu1kzKCqK6Pl2js9HPE2z+29l71++w2680bslXNOm8R6+JCitFBWJg/DWw8LC0m2I06fDYzd/zDdtTuPIkVez+udD+UO99wj85VGFuVRIgS4SByWthxkZULfub4+b1NnO0VMH8JdJx1N37Spusil0ZjGLi7JZuDDeo5ZEp5KLSBxE9nnjHFvun8Gpb9zGPru+Zwq9uMNG81Pm/tQJanGPREeBLhInv/Z5f/SR173y73/zc/vOnLr2dRYXZVOvHkwOLevX4h6JhgJdJF62bYPhw709V/bfHx57jH2uvZZxS+pohaZUiwJdJNaCQZgxA267DbZsgd69YdQoL9TRCk2pPgW6SCx9+CH07Qvvvuul9uzZ3h4sUajw9nIiKNBFasUe4fvjjzBsGDz8sNd6+MQT3p2U60TXaKb9ySUaCnQRn4WHb/26QT66ZRpHTB3sXd3s0wdGjoQmTar0ntqfXKKhQBcJ8aukURK+nYqX81BxX44YvRhOOgnmzPFu0FwN2p9coqFAFyH6kkY0oX/68T+wP8O4kYfZwgGs/dsTHDEy+vJKWbQ/uURDgS5CdCWNSkM/GIS//50Thgwh2/3AspxcgiNGcuKZjX0Zo7pfpDJa+i9C6aX45ZU0ygr9Eh8/sYyNrU6CG2+Etm2x5cvJfneib2EuEg0Fugi/lTRGjSq/3BIZ+k2bQt6wraw6rTcdrjuBjPVfcUPdaQTGvgPHHBPzcxBRyUUkpLKSRngdu2mTIB/lPs5dhUNozDYmcTPDuYsdwf04/G3IOSlmwxb5lQJdpAzlXfzMyYGczPfZdHFfehW+zyL+QD97kNWZHQlqEy2JMwW6SIRyL35u3QpDh8LUqRyw/8FcW/cpZhT/mXr1TZtoSUJQoItEiLz4OWNakG1jH+P0+UOou2M79O9PvREj6LVqP9osrFqIx3L5vrYKSD9RBbqZDQRuABywArjWObc77PX6wHTg98BW4HLn3Fe+j1akmqoSbuGLeDrbe1w3pS/Zbinv1DmVfaY9yLFXdQCq3kYYy+X72iogPVXa5WJmLYCbgWznXAcgA+gRcdj1wI/OuSOA+4Exfg9UpLpKwm3YMO9rIFDx8Tk58PZLW3j/+F68U9SF5m4jf+ZpurGQ2es7VHscFbU9+i2WnyWJI9qSSybQwMwKgYbApojXLwBGhB6/CEw2M3POOV9GKVINJbPydeuqsA9KcTFMncoJQ4fCTz+x6YpbOO7l4Wwt3LfGFzxjuXxfWwWkp0oD3Tm30czGA+uAXcAc59yciMNaAOtDxxeZ2XagKbAl/CAz6wX0AsjKyqr56EXKEV5yyMiAzNB/6SX94/fcU0b5ZckSb/Os5cu9FydPpnn79rxSw1p0eLknVsv3tVVAeqo00M2sCd4MvDWwDXjBzK5yzj1V1Q9zzk0BpgBkZ2dr9i61JrzkAN4CzqwsL8wHDIioLR+xGW6/Hf7+d2jeHJ59Fi6/HMyAmi25L6uWPWSIP+dYGW0VkH6iWSl6OvClc26zc64QmAlELpvYCBwKYGaZwH54F0dF4iJyVWfPnl6Qbt36W9AX5Rfz070PQZs2MH063HorrFkDPXr8GuY1pVq2xFI0NfR1QBcza4hXcukOLI04ZhZwDRAALgXmq34u8VReyaEk6I/PDzDZ9eXYWR9At24weTIcfbTv41AtW2Ipmhr6EjN7EVgOFAEfAFPMbCSw1Dk3C3gcmGFma4Ef2LMLRiQmItsTI0sOOYd/z1fdBnPQm0+Sf2ALmPQcXHaZbzPySKplSyxZvCbS2dnZbunSyIm+SPVV2HtdVOTd/m3YMNixA265xXvcqFFcxyxSVWa2zDmXXdZr2m1RElYg4HWjlNU3XtZr5dar//MfyM6Gm2+GE06AFStgzBiFuaQcLf2XhFTRbDvytbzQPipNm5auV5/R6Tv4y2CYNg1atoQXXoBLLilVXtHyeEklCnRJSBXdQSj8tfx8yM3l150O8/Lgh++LuHzrQ7S+cjjs3Om1JN5xB+y9d6nP0PJ4STUKdElIFXWHhL9m5gV7MOh932DZv7l9cV/4+GM44wyYNAnati3zM6K57ZxIMlGgS0KqqDuk1I0mQguFmuR/y1hu48opM+DQQ+Gll+CiiyrsXlFLoaQadblIcisq4stBD9LskeHUK95FnVsHwd/+tkd5pTyqoUuyqajLRTN0SV6LFkFuLq1XrICzzoIHHvBWfVaBlsdLKlHbosRMRW2IVfLNN3DVVXDaafDTT/Dyy/Dmm1UOc5FUoxm6xERNO0oCAVg0r5DLN0+m1RN3eu0td9zhbdDSsGHtDVwkiSjQJSZq0lESCMCdXd9mQkFfWrGSH3POocn0B+CII2pzyCJJRyUXiYnI3Q+j7ijZtIl9/9+VzCnoyt7s4MI6r/LIua8rzEXKoBm6xESVN6kqLISJE2HECNoVFDI6czj3BAdTXL8hg7vFYMAiSUiBLjETdUfJggXe8s9Vq+Ccc6gzcSLdvj8cW6j2QpGKKNAlcWzcCIMGwXPPQatW8OqrcN55YEbO4Qpykcqohi7xV1AA48bBUUd5LYh33unNzs8/v9b2KRdJRZqhS3zNm+eVV9as8WbjeXlw2GHxHpVIUtIMXeJi2asbWNXpcjj9dG+G/tprMGuWwlykBhToKcS3lZi1qaCAr/uMoe2FR9F6xSxGZd7F4sdXwrnnAklyDiIJSiWXFFGbe3v7toHV3LmQm8vvPv2UWZxPf/JY71qTGYAuXbU/uUhNaYaeIsq9/VoNlYTssGHe12rNnNevh//9X29/8qIiVo9/gx4NXmV9RutSi4xq6xxE0oVm6Cmitvb2rtFNIPLz4f77YdQocM77OmgQ7fbai3kn7Tnr1/7kIjWjQE8RVV6JGaXqhGwgAF9PncMF8/rRYN1ncOGFXrC3alVqvJFjrK1zEEkXusGFVKoqNfRlL69j/aUDuTA4k7V2BIX3TaLdwLNiMUyRtKAbXEiNRLVkPz8f7ruPTiPupl0Q/sbd3G+DGLa7Pu1iMkoRUaBLtZXM3C/c6y3aPXwzfP45P3W9mJMCE/hv0e9UBxeJMQW6VEsgAH/p9jX35g+kHS+zK6sNDd56i6Z/+hNP6j6dInGhQJeq272bguHj+SB/NA5jqN1D4xsGctuf6gO6T6dIvKgPXapm9mzo0IHT5g7jn3XOoUOd1eTtdTt/OL1+td9Sq0NF/KEZukTnq69gwABvS9u2bWHOHA5pdAa9FtastKLVoSL+UaDLr8psT9y9G8aO9abQGRlw770wcCDUq0cONQ/fGi1cEpFSFOhpLDzAofRMOS8Pmvzndc6b25+9Nn3hLd2fMAFatvR1DFodKuIfBXqaiix1XHPNbzPl5ru/oHnvAZzrXmONHUXxxLm079e9Vsah1aEi/lGgp6nIUgfAvnV30T84lsHuHorI5FbGMtn6M/yXerSv5P1qsiOjumJE/KFAT1ORpY7+h73G/fsPYK9NX7D2hB6ctWI8XxW2iKoMogubIolBbYtpqqTUMWngF2w4/jyOuu189tqvPsyfzxHvPcuM+S0YNSq6cK5o21u1JIrEjmboSciXG07s2kXOW/eSM2kM1K0L48fDzTd7j6laGaS8C5uauYvElgI9ydQ0JAPvOjY+PItz5w1gr2++giuu8MK8efNqj6m8C5tqSRSJrUoD3czaAs+HPXUYMNw5lxd2TFfgVeDL0FMznXMjfRynhEQbkmXN4j94YS3be/Tn0uCbrLT2MHkB7ft29WVcZc3o1ZIoEluVBrpz7lPgWAAzywA2Ai+Xceg7zrlz/R2eRIomJCNn8Qve2Enn+ffQ8Z6xHB6sz0Am8LDlcudPdSvtXqkJtSSKxFZVSy7dgf86576ujcFI5aIJyd9m8Y6z8l+l7cUDYNvX/HjmlXReNI51hc1iNmNWS6JI7FQ10HsAz5bzWo6ZfQRsAgY551ZGHmBmvYBeAFlZWVX8aCkRHpJllVa6doV2mZ8zrvhmzgq+xc79O8ArCznwtNN4WlvbiqSsqG9BZ2b18MK6vXPuu4jX9gWCzrlfzOwc4AHn3JEVvZ9uQVdzZV4g7bQDRo8mOG48BXX24pved9F6XN9fu1dEJLlVdAu6qvShnw0sjwxzAOfcT865X0KP3wTqmtkB1RqtRK3UBdJ8x6bJM+Hoo2H0aOr0uJy9vvqU1nkDfA1z9ZWLJK6qlFyuoJxyi5kdAnznnHNmdiLeD4qtPoxPKlBygfR3+Z8xkX6c8cwc6NQJnn4aTjnF989TX7lIYotqhm5mewNnADPDnuttZr1D314KfBKqoU8EerhoazlSbTmddvD5pUP4xDrQreESmDgRli2rlTCHileEikj8RTVDd87tAJpGPPdI2OPJwGR/hyblcg5efBFuuYUWGzZ4WyWOGQMHH1yrH6u+cpHEppWiyWbNGujXD+bOhWOOgeefh5NOislHq69cJLEp0JPFL7/A3Xd7N5lo2BAmTYLevSEztv+E6isXSVwK9AQWCMDCBY5Lgi/Q5tG/woYNcO213m3gDjoo3sMTkQSjQE9QgQD06baa8fn9aMM8fmlzHI3e/Ue1pse+7M4oIglPgZ5gAgF4958/0/GVUbyXfz+/0Ihce5CWPW/i9pyMar2fWg1F0oMCPYEE3nU81PUf3Ft4Cy3YxBN1rmcI9/BT/QOZ98fqvae2sBVJHwr0ONmjDLJyJc2u7seMwgUs43guq/MSnXp1oX9WzUolajUUSR8K9DgIL4PsX/dnPrzoLpq/8AAtG+xDv7qP8GjxDWTWz2B8z5rPptVqKJI+FOhxsHCht/fKZcFnGV88iGbPfQs33EDm6NH8+fMDaL7Q3/BVq6FIelCg+yy8lAJlz4zPPvQTTiaXU3mbZZbN1qmv0PH6EwHIOUDhKyLVo0D3UXgpJSMDzKCoKKy7pP1PMGIEx06cSGGj/Zj9x0dp/Nfr+f0pVe9eERGJpED3UXhHSTDoPeecV17ZnPcMLBoE330HN95I3dGjObtp0wrfL5x6yUWkMgp0H4V3lJTM0I8qXMEkcvnDPxbBCSfArFne1ypQL7mIRKMqN7hIS1W5oUNJR8moUfDO69v54sKBLOc4uuyzEqZMgcWLqxzmoG1rRSQ6mqFXoDoz45wujpzPZ8DVt8H338NNN1Hn7ruhCuWVSOolF5FoKNArUOVVlh99BLm58O9/Q+fO8PrrkF3mrf9Kqaw+rl5yEYlGWgR6dS8oRj0z3rYNhg+HBx+E/feHxx7zdkWsU3lFK9rfAtRLLiKVSflAr8kFxUpnxsEgzJgBt90Gmzd7+5PffbcX6lHSXisi4peUDfSSWfm6dTULzHJnxh99BH37wn/+A126wOzZcPzxVR6n6uMi4peUDPTIBT4lN/XxJTC3bYNhw+Chh7wLnX//u3dPzyjKK2VRfVxE/JKSgR5exgC48UbIquGuhQSDMG0aDB4MW7dCnz4wciQ0aVLj8ao+LiJ+SMlAjyxj9KzproXLl3vllcWLvTeaM4fArmNZ+EjVfkhotaeI1KaUDHTfyhg//gh33AGPPOKVV558Eq6+msCSOlW+0KrVniJS21Iy0KGGZYxg0AvvwYPhhx+82fnIkdC4MVC9zhR1s4hIbdPS/0jLl8PJJ8P118NRR3nfT5z4a5jDbyWdjIzoL7RW5++IiFRFys7Qq+yHH+Bvf4NHH4WDDoLp0+Gqq7wdtiJUp6SjbhYRqW3mnIvLB2dnZ7ulS5fG5bNLCQbh8cdhyBCvJTE3F+66C/bbL94jExHZg5ktc86VuadI2s3QS3Wa1F3q1cffew9OOcVbut+pU7ndKOpSEZFEllaBXtJp0ih/KwfZULoEp2IHH+wt37/ySjArtxtFXSoikujS6qLo2/OL6bl7CquDbbim+HHmdRzAhF5rCBz+W628vL3HtSe5iCS69Jmhv/8+uU/1oZFbyiJOpX/mg6z+tANF/wf1xkFenrcAtGnTsvdW0Z4rIpLoUj/Qt2yBoUPhscdodMghfD7iaf5T9wq6rDdWTPVm3Pn53rXQYNAL65JwD6+Vq0tFRBJd6gZ6cbG3L/nQobB9OwwcCHfeyZH77ssQvJr4tGnejNvstxs7FxR4YT5kyJ5vqT1XRCSRpWagL1niTbmXLvWm05MnQ/v2pQ4Jn3E3bQoDBqicIiLJLbUCffNmb2r9+OPQvDk88wz06FHm4iAoPePu2LH65RS1M4pIIkiNQC8uhilTvJWeP/8MgwZ5t4TbZ5+o36K65RS1M4pIoqi0bdHM2prZh2F/fjKzARHHmJlNNLO1ZvaxmVX91j3VFQjAiSd6+5Mfd5x3J6Fx46oU5jWhdkYRSRSVBrpz7lPn3LHOuWOB3wM7gZcjDjsbODL0pxfwsE+rxM0AAAYtSURBVN8D3cP338N118FJJ8F338Hzz8PcuXD00bX+0eG06ZaIJIqqlly6A/91zn0d8fwFwHTnbQyz2Mwam1kz59w3vowyXHGxtz/5HXfAL794N2geNgwaNfL9o6KhdkYRSRRVDfQewLNlPN8CWB/2/YbQc/4H+hNPeB0sp58OkyZ5W9zGmdoZRSQRRB3oZlYPOB8oo0M76vfohVeSISsrq3pv0rMnHHwwnHtuud0rIiLpqCp7uZwNLHfOfVfGaxuBQ8O+bxl6rhTn3BTnXLZzLvvAAw+s2khDAsvqcc8n5xFYrDAXEQlXlZLLFZRdbgGYBeSa2XNAZ2B7bdTP1SIoIlK+qGboZrY3cAYwM+y53mbWO/Ttm8AXwFpgKtDH53ECahEUEalIVDN059wOoGnEc4+EPXZAX3+HtifteCgiUr6kWilaWYugluCLSDpLqkCH8lsEVV8XkXSXMncsUn1dRNJdygS6luCLSLpLupJLebQEX0TSXcoEOmgJvoikt5QpuYiIpLukDvRAAO65x/sqIpLukrbkojZFEZHSknaGrjZFEZHSkjbQ1aYoIlJa0pZc1KYoIlJa0gY6qE1RRCRc0pZcRESkNAW6iEiKUKCLiKQIBbqISIpQoIuIpAgFuohIijDvdqBx+GCzzcDX1fzrBwBbfBxOMtA5pwedc3qoyTn/zjl3YFkvxC3Qa8LMljrnsuM9jljSOacHnXN6qK1zVslFRCRFKNBFRFJEsgb6lHgPIA50zulB55weauWck7KGLiIie0rWGbqIiERQoIuIpIiEDXQzO9TMFpjZKjNbaWb9yzjGzGyima01s4/N7Ph4jNUvUZ7zlaFzXWFm75rZMfEYq1+iOeewY08wsyIzuzSWY/RbtOdsZl3N7MPQMW/Hepx+ivK/7f3M7DUz+yh0zLXxGKtfzGwvM3sv7HzuKuOY+mb2fCjDlphZqxp9qHMuIf8AzYDjQ4/3AT4Djo445hxgNmBAF2BJvMcdg3M+CWgSenx2Opxz6LUMYD7wJnBpvMcdg3/nxsAqICv0/UHxHncMznkoMCb0+EDgB6BevMdeg3M2oFHocV1gCdAl4pg+wCOhxz2A52vymQk7Q3fOfeOcWx56/DOwGmgRcdgFwHTnWQw0NrNmMR6qb6I5Z+fcu865H0PfLgZaxnaU/ory3xmgH/AS8H0Mh1crojznPwMznXPrQscl9XlHec4O2MfMDGiEF+hFMR2oj0K59Evo27qhP5FdKBcA00KPXwS6h86/WhI20MOFfg05Du8nXLgWwPqw7zdQdhgknQrOOdz1eL+hpITyztnMWgAXAQ/HflS1q4J/5zZAEzNbaGbLzKxnrMdWWyo458lAO2ATsALo75wLxnRwPjOzDDP7EG8i8i/nXLkZ5pwrArYDTav7eQl/Czoza4Q3MxvgnPsp3uOJhWjO2cy64QX6KbEcW22p5JzzgMHOuWANJi8Jp5JzzgR+D3QHGgABM1vsnPssxsP0VSXn/CfgQ+CPwOHAv8zsnWT+/945Vwwca2aNgZfNrINz7pPa+ryEnqGbWV28f/ynnXMzyzhkI3Bo2PctQ88lrSjOGTPrBDwGXOCc2xrL8dWGKM45G3jOzL4CLgUeMrMLYzhE30VxzhuAfzrndjjntgCLgGS/AF7ZOV+LV2Zyzrm1wJfAUbEcY21xzm0DFgBnRbz0a4aZWSawH1Dt/6cTNtBDdaTHgdXOuQnlHDYL6BnqdukCbHfOfROzQfosmnM2syxgJnB1ss/WILpzds61ds61cs61wqsz9nHOvRLDYfoqyv+2XwVOMbNMM2sIdMarOyelKM95Hd5vJJjZwUBb4IvYjNB/ZnZgaGaOmTUAzgDWRBw2C7gm9PhSYL4LXSGtjkQuuZwMXA2sCNWgwLsKngXgnHsEr+PhHGAtsBPvJ3wyi+ach+PV2B4KlR+KXHLvVBfNOaeaSs/ZObfazN4CPgaCwGO1+at6DETz7zwKeNLMVuB1iAwO/XaSrJoB08wsA2/y/A/n3OtmNhJY6pybhfdDboaZrcW7CNyjJh+opf8iIikiYUsuIiJSNQp0EZEUoUAXEUkRCnQRkRShQBcRSREKdBGRFKFAFxFJEf8fIiUGDfDYH34AAAAASUVORK5CYII=\n", "text/plain": [ "