{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Robust Linear Models" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import statsmodels.api as sm\n", "import matplotlib.pyplot as plt\n", "from statsmodels.sandbox.regression.predstd import wls_prediction_std" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimation\n", "\n", "Load data:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "data = sm.datasets.stackloss.load(as_pandas=False)\n", "data.exog = sm.add_constant(data.exog)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Huber's T norm with the (default) median absolute deviation scaling" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-41.02649835 0.82938433 0.92606597 -0.12784672]\n", "[9.79189854 0.11100521 0.30293016 0.12864961]\n", " Robust linear Model Regression Results \n", "==============================================================================\n", "Dep. Variable: y No. Observations: 21\n", "Model: RLM Df Residuals: 17\n", "Method: IRLS Df Model: 3\n", "Norm: HuberT \n", "Scale Est.: mad \n", "Cov Type: H1 \n", "Date: Fri, 21 Feb 2020 \n", "Time: 13:56:28 \n", "No. Iterations: 19 \n", "==============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "var_0 -41.0265 9.792 -4.190 0.000 -60.218 -21.835\n", "var_1 0.8294 0.111 7.472 0.000 0.612 1.047\n", "var_2 0.9261 0.303 3.057 0.002 0.332 1.520\n", "var_3 -0.1278 0.129 -0.994 0.320 -0.380 0.124\n", "==============================================================================\n", "\n", "If the model instance has been used for another fit with different fit parameters, then the fit options might not be the correct ones anymore .\n" ] } ], "source": [ "huber_t = sm.RLM(data.endog, data.exog, M=sm.robust.norms.HuberT())\n", "hub_results = huber_t.fit()\n", "print(hub_results.params)\n", "print(hub_results.bse)\n", "print(hub_results.summary(yname='y',\n", " xname=['var_%d' % i for i in range(len(hub_results.params))]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Huber's T norm with 'H2' covariance matrix" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-41.02649835 0.82938433 0.92606597 -0.12784672]\n", "[9.08950419 0.11945975 0.32235497 0.11796313]\n" ] } ], "source": [ "hub_results2 = huber_t.fit(cov=\"H2\")\n", "print(hub_results2.params)\n", "print(hub_results2.bse)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Andrew's Wave norm with Huber's Proposal 2 scaling and 'H3' covariance matrix" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parameters: [-40.8817957 0.79276138 1.04857556 -0.13360865]\n" ] } ], "source": [ "andrew_mod = sm.RLM(data.endog, data.exog, M=sm.robust.norms.AndrewWave())\n", "andrew_results = andrew_mod.fit(scale_est=sm.robust.scale.HuberScale(), cov=\"H3\")\n", "print('Parameters: ', andrew_results.params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "See ``help(sm.RLM.fit)`` for more options and ``module sm.robust.scale`` for scale options\n", "\n", "## Comparing OLS and RLM\n", "\n", "Artificial data with outliers:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "nsample = 50\n", "x1 = np.linspace(0, 20, nsample)\n", "X = np.column_stack((x1, (x1-5)**2))\n", "X = sm.add_constant(X)\n", "sig = 0.3 # smaller error variance makes OLS<->RLM contrast bigger\n", "beta = [5, 0.5, -0.0]\n", "y_true2 = np.dot(X, beta)\n", "y2 = y_true2 + sig*1. * np.random.normal(size=nsample)\n", "y2[[39,41,43,45,48]] -= 5 # add some outliers (10% of nsample)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 1: quadratic function with linear truth\n", "\n", "Note that the quadratic term in OLS regression will capture outlier effects. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 5.11003863 0.5285205 -0.01374639]\n", "[0.45516517 0.07027136 0.00621793]\n", "[ 4.76637875 5.03591903 5.30087908 5.56125892 5.81705852 6.06827791\n", " 6.31491707 6.55697601 6.79445472 7.02735321 7.25567148 7.47940952\n", " 7.69856735 7.91314494 8.12314232 8.32855947 8.5293964 8.7256531\n", " 8.91732959 9.10442585 9.28694188 9.46487769 9.63823328 9.80700865\n", " 9.97120379 10.13081871 10.28585341 10.43630788 10.58218213 10.72347616\n", " 10.86018996 10.99232354 11.1198769 11.24285003 11.36124294 11.47505563\n", " 11.58428809 11.68894033 11.78901235 11.88450414 11.97541572 12.06174706\n", " 12.14349819 12.22066909 12.29325977 12.36127022 12.42470045 12.48355046\n", " 12.53782025 12.58750981]\n" ] } ], "source": [ "res = sm.OLS(y2, X).fit()\n", "print(res.params)\n", "print(res.bse)\n", "print(res.predict())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Estimate RLM:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 5.04611616e+00 5.11413247e-01 -2.69305305e-03]\n", "[0.13964091 0.02155867 0.00190761]\n" ] } ], "source": [ "resrlm = sm.RLM(y2, X).fit()\n", "print(resrlm.params)\n", "print(resrlm.bse)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Draw a plot to compare OLS estimates to the robust estimates:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAHSCAYAAADlm6P3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd1zVZRvH8c/vHKaKew/EiXvi3rtMc1SWWTYsSystU9OnYWWmhaamaKamWaaZq9zmtsQEV05c4cBcKArKPOf3/HGJSKIiAodxvV+v80J+5xy4z5M+fLnPdV+XYZomSimllFJKZQcWRy9AKaWUUkqp9KLhVymllFJKZRsafpVSSimlVLah4VcppZRSSmUbGn6VUkoppVS2oeFXKaWUUkplG07p+c0KFixoenl5pee3VEoppZRS2dDOnTsvmaZZ6L/X0zX8enl5ERgYmJ7fUimllFJKZUOGYZxM6rqWPSillFJKqWzjvuHXMIzvDMO4YBjG/v9cf8swjCDDMA4YhvFl2i1RKaWUUkqp1JGcnd/ZwCO3XzAMoxXQBahhmmZVYGzqL00ppZRSSqnUdd+aX9M0txiG4fWfy/2AMaZpRt98zIWULiA2NpYzZ84QFRWV0i+RKbi5uVGyZEmcnZ0dvRSllFJKqWwrpQfeKgLNDMMYBUQBg03TDEjJFzpz5gweHh54eXlhGEYKl5OxmaZJaGgoZ86coUyZMo5ejlJKKaVUtpXSA29OQD6gITAEWGDcJbkahtHXMIxAwzACL168eMf9UVFRFChQIMsGXwDDMChQoECW391WSimllMroUhp+zwCLTbEDsAMFk3qgaZrfmqbpY5qmT6FCd7RaA8jSwTdedniNSimllFIZXUrD71KgNYBhGBUBF+BSai3KkT7++GPGjr37+b2lS5dy8ODBdFyRUkoppZRKLclpdTYP8Ae8DcM4YxhGH+A7oOzN9mfzgRdM0zTTdqli6e4QmozZQJlhK2gyZgNLd4ekx7dN+P4afpVSSimlMq37hl/TNHuaplnMNE1n0zRLmqY50zTNGNM0nzNNs5ppmnVM09yQHotdujuE4Yv3ERIWiQmEhEUyfPG+hw7Ao0aNwtvbm7Zt2xIUFATA9OnTqVevHjVr1uSJJ57gxo0bbNu2jd9++40hQ4ZQq1Ytjh8/nuTjlFJKKaVUxpSpJrz5rgkiMtaW6FpkrA3fNUEp/po7d+5k/vz57N69m8WLFxMQIE0runfvTkBAAHv37qVy5crMnDmTxo0b8/jjj+Pr68uePXsoV65cko9TSimllFIZU0pbnTnE2bDIB7qeHFu3bqVbt27kyJEDgMcffxyA/fv388EHHxAWFkZERAQdOnRI8vnJfZxSSimllHK8TLXzWzyv+wNdT66kOjG8+OKLTJ48mX379jFixIi7tilL7uOUUkoppZTjZarwO6SDN+7O1kTX3J2tDOngneKv2bx5c5YsWUJkZCTh4eEsW7YMgPDwcIoVK0ZsbCxz58699XgPDw/Cw8NvfX63xymllFJKqYwnU5U9dK1dApDa37NhkRTP686QDt63rqdEnTp1ePrpp6lVqxalS5emWbNmAIwcOZIGDRpQunRpqlevfivwPvPMM7z66qt8/fXXLFy48K6PU0oppZRSGY+RTh3KAPDx8TEDAwMTXTt06BCVK1dOtzU4UnZ6rUoppZRSjmQYxk7TNH3+ez1TlT0opZRSSqlMwm539AqSpOFXKaWUUko9nMhICAiAvXvl8ytXIHduuJTxBgBnqppfpZRSSimVQUyaBH/9BXv2wOHDYLPB00/D/PmQLx+89RbExjp6lXfQ8KuUUkoppe5kmhAcLOF29275mCOHhFuAWbPg4kWoVQu6d5ePdesmPH/0aIcs+340/CqllFJKZXexsbJ7e+QIPPGEXOvRAxYulD9bLFCpEjRvnvCcP/6QMJzJaPhVSimllMqONm6EuXNlR3f/foiOluvXroGHBzz/PLRtKzu61avfGXQzYfAFDb+EhobSpk0bAM6dO4fVaqVQoUIA7NixAxcXF0cuTymllFIq5S5elJKF228rVkD58nDwICxdCrVrw4ABEnJr14acOeW5jz/u2LWnkWwffgsUKMCePXsA+Pjjj8mVKxeDBw9O9BjTNDFNE4tFm2MopZRSKgMyTTh5UsJtrVpQpgwsXw6dOyc8pnRpCbfxh9Beew369wfDcMyaHUTT3F0cO3aMatWq8frrr1OnTh1Onz5N3rx5b90/f/58XnnlFQDOnz9P9+7d8fHxoX79+mzfvt1Ry1ZKKaVUdhEaCoMHQ5s2UKCABN7u3WHZMrm/bl0YOxbWr5fHBgfDkiUQP3DLySnbBV/IYDu/b78tZSepqVYtmDAhZc89ePAgs2bN4ptvviEuLu6ujxswYABDhw6lYcOGBAcH06lTJ/bv35/CFSullFJK3RQdDQcOyI7url3y8dFH4cMPwdUVpk6FqlXhqadkV7dOHanPBShWDN5917Hrz4AyVPjNaMqVK0e9evXu+7h169YRFBR06/MrV64QGRmJu7t7Wi5PKaWUUlnJjRsyJOLGDdnNNU3w9IQLF+R+Dw8JuIULy+e5csHVq7KDq5ItQ/2vldId2rSSM77gG7BYLJimeevzqKioW382TVMPxymllFLqwf34I6xdK7u6hw7JSOAaNSQEGwaMHAl588qObtmy0nLsdhp8H5jW/CaTxWIhX758HD16FLvdzpIlS27d17ZtW/z8/G59vie1azeUUkoplXlduQIbNkj9bc+e0Lix7OoCrFwJ69aBlxd88IF0X1i+POG5fftKv93y5e8MvipF9NeFB/DFF1/wyCOP4OnpSZUqVYi+2Q/Pz8+Pfv36MWvWLOLi4mjVqlWiMKyUUkqpbOLKFdnFbdYMXFxk5/ajjxLu9/SUXdzISOmT+9134ObmuPVmQ8btb+WnNR8fHzMwMDDRtUOHDlE5/tRhFpedXqtSSimVLfzzDyxYADt3yu3ECbkeGCjdFrZulUlodetK6C1Y0LHrzUYMw9hpmqbPf6/rzq9SSiml1P1cvSo7uoGBcuvXD1q2lLA7bJiULdStC6++Kh+9veV5zZrJTWUYGn6VUkoppW53/bp0XChUCM6ehVat4MiRhPu9vOCJJ+TPTZvCpUvSZ1clcumSnOHLaNlfw69SSimlsrcdO2Q3NyBAPh48KNPPpkyBIkWk+8Lzz0O9erKre3vpgqur3NQtp0/L2b7p0yFPHjh1CpydHb2qBBp+lVJKKZU92Gxw+LCE3NhYKVEA6aZw8qTs9Pr4yJS0tm3lPqsVfvnFcWvORA4fhi++kO5tAL16wXvvZazgCxp+lVJKKZXVTZkiATYwECIi5FrVqgnhd/58KF4cSpXKluN+H1ZgIIweLZOT3dygf38ZLOfp6eiVJU3Dr1JKKaUyv8uXpXwh/nbgABw9KkMgjh+X1mIvvgj160v5QsWKCc9t2NBhy04NS3eH4LsmiLNhkRTP686QDt50rV0iTb+nacLGjRJ6162TORzvvw8DBsgGekam3ZKBM2fO0KVLFypUqEC5cuUYOHAgMTExbNq0iU6dOt3x+OXLl1O7dm1q1qxJlSpVmDZtmgNWrZRSSmVTMTFSunD9unw+caIcOHv0Ufj4Y2k/1rIlhIfL/ePGwfbtMGmS1O5WqpRlBkYs3R3C8MX7CAmLxARCwiIZvngfS3eHpMn3s9tlh7dhQ5nAvH8/fPmlVI2MHJnxgy9o+MU0Tbp3707Xrl05evQoR44cISIigvfffz/Jx8fGxtK3b1+WLVvG3r172b17Ny1btkzfRSullFLZyZUrUprwzjsyHS13btnB3bZN7m/WDD7/HNavh7Aw2fWdNQvy5XPsutOB75ogImNtia5FxtrwXROUqt8nNha+/x6qVZOS6EuX4Jtv5PeMIUPkP0lmke3LHjZs2ICbmxsvvfQSAFarlfHjx1OmTBlatWp1x+PDw8OJi4ujwM2WJq6urnjH9/JTSiml1MOJiJBdXX9/aSPWvLn00u3ZE9zdpdvCm29CgwYyNALkY/yfs5mzYZEPdP1B3bgBM2dK94ZTp6TxxU8/wVNPSUXJ3TiiFCO5Mlb4fftt2LMndb9mrVowYcJd7z5w4AB169ZNdC137tx4enpy7NixOx6fP39+Hn/8cUqXLk2bNm3o1KkTPXv2xJJF3j5RSiml0l10tGQAf3/Yt0/eWwf49FMJvzVqyICJatUyXusAByue152QJIJu8bzuD/V1w8LAz08qSi5ehCZN5Nxgx473PxMYX4oRvyMdX4oBZIgAnO0Tm2maGEn8V7zbdYAZM2awfv166tevz9ixY3n55ZfTeplKKaVU5hceLqUJn30GnTrJDi6Aiwts2gSFC8MHH8DKlRAaCh9+KPc7O0Pt2hp8kzCkgzfuztZE19ydrQzpkLJ3pc+dk/Zknp7yn8LHB7ZskQnNjz2WvGYY8aUY0ZZDXHVaQLTlUJqUYqRUxtr5vccObVqpWrUqixYtSnTt2rVrnD59mnLlyt31edWrV6d69eo8//zzlClThtmzZ6fxSpVSSqlMxDRly7BwYfn86adh4cKEXd3KlaFmTfmzYchgCW0z9sDid1IftsTgxAnw9ZVS6dhYKWsYNkzeQH9QZ8MiuW75k0suXwI2DFwoEjOKs2GVH/yLpYGMFX4doE2bNgwbNow5c+bQu3dvbDYb7777Li+++CI5cuS44/EREREEBgbeOuS2Z88eSpcunc6rVkoppTKYmBjYvVsOof35p3y8fl0Oq1ksUr9buTI0aiT1unnzJn6+Bt8U61q7RIrLCfbtgzFj5Dyhk5N0gxsyBMqXf/CvZbPbWHl0JWE5RxFm+0suGmCacURZ9lE2d8aoy8724dcwDJYsWUL//v0ZOXIkdrudjh078vnnn+Pv78/69espWbLkrcfPmzePL7/8ktdeew13d3dy5sypu75KKaWyn9BQCbht2kCOHFKfO2qU3OflBa1aSaFobKyM/33rLYcuVyW2bZv06F2+HHLmlEYagwbJrI8H9W/4v8zcPZNvd37L6WunyedahNxx7Qg3NmOacRg4kceoleJSjNSW7cMvQKlSpVi2bNkd11u2bElk5J1F5M2aNUuPZSmllFIZR2gorFolxZ9bt0qZAkitbosW0o2hVi1pRZaSBKXSnGnC6tUSerduldbIn3wipdf58z/o1zLZGLyRqYFTWXp4KXH2ONqWbcuERybQuWJnVvx9gY9WLeb0jUBK5fDh00e7Z4jDbqDhVymllFL/ZbfL9IKtW2UaWv36Mi3t+eeloWuTJtCrl5Qy1K8vz6laVW4qw7HZpNx6zBhpqlWgSBxenU5gr3CCNS4u1DjpTdf89w+m/qf9WXlsJVejrrLm+BqOhB4hv3t+BjYYyGt1X6NCgQq3HiulGBlzt1/Dr1JKKaWk3dj48RJ4//wTrl6V6x98IAG3bl1JTtWqgdV676+lMoToaJgzRyawHTsmE53fHHGF32N2EGWPwyB5bchM02TGrhn0W9EPmynty6oXrs4P3X7gySpP4ubkll4vKVVo+FVKKaWymxs3ZNzv5s3g4QGDB0u7sfHjoWBB6czQtKlMTos/1O3snNCdQWVoEREwbRp89RWcPSvzP375Bbp1g+a+u4mKikv0+Pg2ZP8NvxExEfy07yemBk5lz7mEOQxWw0rPaj15rsZz6fJ6UpuGX6WUUiq7mDgRFiyQCWqxsdKFoUsXuc8wpN9VzpyOXaNKsdBQ+PprmDRJmmy0agWzZ0PbtgnNNJIzEe7AhQNMDZzKnL1zCI8Jp2aRmgxtPJRJOyYRY4vBxepCS6+Waf+C0oiGX6WUUiqruXJFyhe2bIG9e2HtWkk/QUFy6mnQIDmk1rgx5MmT8DwNvpnSmTMwbhx8+61s6nfpIj16Gza887F3mwhXNI+VefvmMTVwKltPbcXV6kqPqj3o59OPhiUbYhgGXSt1ZVPwJlp6taRRqUbp8MrShoZfpZRSKqv47Tc5vr97t4RcV1dJQJcvy9F+Pz/tp5uFBAVJPe8PP8gZxWeflels9zp3OKSDN8MX7yPMtp8oyz6c7CWwOx/nCut5dnEo5fKV48u2X/JS7ZcomKNgouc2KtUoU4feeNl+vDGA1WqlVq1aVKtWjc6dOxMWFgZAcHAw1apVu+Px8QMwwsPDb10bOHAghmFw6dKldFu3UkqpbCo8XNqODR0qh9G2b5frLi7SjWHECKnnDQuTVmQFCsj9GnyzhJ07ZQJb5crw00/Qt68caJsz5/4NN7rWLsHzLSI57/o/wpzmcMl1NFcsv9CqTFNW91rNkbeOMKTJkDuCb1ai4Rdwd3dnz5497N+/n/z58+Pn53ff55QvX55ff/0VALvdzsaNGylRImP0r1NKKZVFBQfLTm6+fNCxI0yYAG5uMl0N4JFHYONGCb/Nm8t9KkswTfk9pkMH8PGRSpZhw+SvxOTJMlfkfi5cv8DoraPxDXwZk1gwwMBgSOMhLH1mKR3Kd8BiZP1omClfof9pf0ZvHY3/af9U/9qNGjUiJCTkvo/r2bMnP//8MwCbNm2iSZMmODlpFYlSSqlUEBMjwyQ++URC7EcfyfWiRSXQDhsGv/8uO7tbtshjVJZkt0s1S+PGcoBt717p13vqFHz+ORQpcu/nm6bJH6f+4NlFz1Lyq5L8b8P/KJe/HC5WF6yGFTcnN7pW6po+LyaDyFBp7e3VbydqpZGUq9FX+fv839hNOxbDQo0iNcjjmueuj69VtBYTHpmQrO9vs9lYv349ffr0ue9jK1SowK+//sqVK1eYN28ezz33HKtWrUrW91FKKaXu6tlnJe1cvy5lCnXqSOgFCb6bNjl0eSp9xMXB/PkSdA8ckJ3dKVPgxRfB3f3+zw+PDufHv39kSuAU9l/YTx7XPPSv15/XfV6nUsFK+J/2zxKH11IiQ4Xf5LgadRW7aQfAbtq5GnX1nuE3OSIjI6lVqxbBwcHUrVuXdu3aJet53bt3Z/78+fz1119MmzbtodaglFIqGzFNOH4c1q+X2/nzUqMLUKgQvPACtGkDLVs++NxZlalFRsKsWeDrKyUN1arBjz9K6+V7vcEcH2ZL5S7FtjPb+OHvH4iIiaB20drM6DyDZ6o9Q06XhG4eWeXwWkpkqPCbnB1a/9P+tJnT5lafubnd5z70f7z4mt+rV6/SqVMn/Pz8GDBgwH2f98wzz1CnTh1eeOEFLJZMWUGilFIqvc2YAZ99BidPyuclS0rQjY2VQRITJzp2fSpVLd0dgu+aIM6GRVI8rztDOngnOUnt6lWYOlXmjFy4IKXdX38Njz0m7ZjvZXPwZtr/2J4Ym9R+O1ucebb6s/Tz6Uf9EvUx9KBjIhkq/CZHo1KNWN97fZps1efJk4evv/6aLl260K9fv/s+3tPTk1GjRtG2bdtUW4NSSqksIiYGtm2T2ty1a+U97HLlZKJa7drSqaFtW6hQQbswZFFLd4cwfPE+ImNlJHBSo4QvXJBzi35+cO2aHGgbPlzKuO/31+Jk2Em+3fktE/+aeCv4Ghi81+Q9RrYemXYvLJO7b/g1DOM7oBNwwTTNav+5bzDgCxQyTTPdenyl5VZ97dq1qVmzJvPnz6dZs2YEBQVRsmTJW/ePHz8+0eNfe+21NFmHUkqpTOrIEXjnHSljuH4drFZo1EgOp4G8f/30045do0oXvmuCbgXfePGjhGvlK8HYsTBzJkRHw5NPyjnGOnXu/TXtpp21x9cyJWAKK46uAKBxqcYEhAQQZ4/DxepCxwod0+olZQnJ2fmdDUwG5tx+0TCMUkA74FTqLyt9RUREJPp82bJlt/4cGxt7x+OfeuqpJL9OcHBwqq5LKaVUBhcaKju7a9bIVt1LL0mN7vHjcjKpfXup282d29ErVQ6Q1CjhmIu52Lu8HOU/kHKG3r3lTYCKFe/9tUJvhDJrzyy+CfyG41eOUzhnYYY3HU7fun3xzOOZrQ+wPaj7hl/TNLcYhuGVxF3jgaHAr6m8JqWUUirjMk0YNQqWL4cdO+Tz/PnB21vuL1gQDh927BpVhnD7KOHos3m5ur0ckUeLYnG2MXCATJm+7c3lW24PslaLlSkBU5i/fz7RtmiaejZlZKuRdK/cHVcn11vPyc4H2B5Uimp+DcN4HAgxTXPv/YqoDcPoC/QFqZFVSimlMpV//5Wd3dOn4cMPpRDz999l227ECBks4eMj5Q1K3WZwe28GjjvLhT/KEH2qIBa3GAo0O8ZXn+Skd6tiST4n/mB/VFwUACYmOZ1z8lKtl+hXrx81itRIz5eQJT1w+DUMIwfwPtA+OY83TfNb4FsAHx8f80G/n1JKKZXudu+Gn3+G1atlqgBA6dJyEsnJSdqT6WAjdRc2GyxZAmPGlOD0zhK45I4mX6uDVGxxkWGPl6dr7aSD77HLxxj6+1Ai4xLKJTpX7MyP3X8kt6uWzqSWlPzLLQeUAeJ3fUsCuwzDqG+a5rmULMI0zSzfhsM0NfcrpVSG9e+/sGoVPPEE5MkjO7vjxkHTpjJl4JFHoEaNhOP3GnxVEmJipCfvF1/IuccKFWD6dHj+eVdcXask+Ryb3cbKoyvxC/BjzfE1WA0rVsOKiYmr1ZXhTYdr8E1lD/yv1zTNfUDh+M8NwwgGfFLa7cHNzY3Q0FAKFCiQZQOwaZqEhobipjPWlVIqY7DZICAAVqyAlSth1y65ni8fdOsGfftCv37Slkyp+7h+XULuuHFw5ox0sluwALp3v3s1zMXrF5m5eybfBH7DyasnKe5RnI9bfMyrdV/lZNhJPbyWhpLT6mwe0BIoaBjGGWCEaZozU2sBJUuW5MyZM1y8eDG1vmSG5ObmlqhlmlJKqXR2+TJERICnJxw9Ku3HLBZo3BhGj4aOHaF6dXls3ryOXavKFC5fhkmTZBjF5cvQooXMMGnfPukevaZp8lfIX/gF+LHgwAJibDG08mrF2PZj6eLdBWerMwDFPYpr6E1Dyen20PM+93s9zAKcnZ0pU6bMw3wJpZRS6k6mCQcOwLJlsru7bRv06gVz5khnhkWLdHywSpGQEPjqK5g2TXZ9O3eWcvBGd8mrG//ZyOQdk9l/YT9HLh/Bw8WDV+u8Sv96/alSKOlyCJV2tGhJKaVU1mGaCVtu7drJwTSQyQH/+x907SqfG4a8J63UAzh6FL78Un5/stngmWfgvfcS3jD4r2OXj/HBhg/4+cDPgExfG9xoMB+1+AgPVy2pcRQNv0oppTK30FA5rLZsGQQGykkjqxV69oQePaBTJyhe3NGrVJnYnj1SGbNwITg7Q58+MGQIJPXGtc1uY9WxVfgF+LH62GoshuXWfRbDQn73/Bp8HUzDr1JKqcxp40b4+GP44w+w26FIEXn/OTxcanb79HH0ClUmZpqwdauE3tWr5ezjkCHw9ttQtOidjw+9EcrM3TOZGjiV4LBgiuUqxogWI6hTrA7PLHyGGFsMLlYXWnq1TPfXohLT8KuUUirjs9ngr79g6VJ4+mmoW1cCb1iYlDN07iyDJiyW+38tpe7BNKUJyOjRUiZeqBB8/rk0/0jqHGTg2UD8AvyYt28e0bZoWpRuwRdtv6BbpW63DrCt771euzdkIBp+lVJKZUxxcbB2rQTe336D8+flPecKFST8tmmTMIBCqYcUFyftycaMgX37ZKbJ5Mnw8svg7p74sZuCN+EX4MeBCwc4dOnQrQlsb9R/g2qFq93xtXX0cMai4VcppVTGcfUqnDwpAyXsdqnbtdulDVm3bvDoozKEQqlUEhUFs2eDry+cOAFVqsiBtmeekd+1bncy7CQfbviQH/b9AMgBtncavMOIliPI46Z/LzMLDb9KKaUcKyREdnaXLpU63vLl4eBBcHGBzZuhcmVwdXX0KlUWc+0afPMNjB8P585BgwbSvqxz58TVM6Zpsu7EOvwC/Fh2ZFmiia0Ww0KhnIU0+GYyGn6VUko5zrBhMgsWpJzhnXekHVl8y7JatRy7PpXlXLwIEyeCn5+UjLdtC3PnQqtWiQdTXI26yvd7v2dKwBSCQoMomKMgw5oMw6e4D70W99IDbJmYhl+llFJpzzRlN3fRIrktWQJly0Lr1lLG0LUrVKqU9FgspVLByZMyfnjGDCl16N4dGna7wJLT++mzNpLiO9wZ0sGbCiXC8AvwY87eOVyPvU6DEg2Y03UOT1V9CjcnN0APsGV2xu3b92nNx8fHDAwMTLfvp5RSysEuXoQJEyTwBgVJuG3cWK75+Dh6dSobOHhQ3lz46Sf5/PnnYehQOBwZwvDF+wiz7SfSshdMg1invdyw/I2r1ZWe1XvyRr038Cmuf08zK8Mwdpqmecd/QN35VUoplXrsdmlJZrdDkyYSdseNkz8PGCCH1ooVc/QqVTawY4e0K1u6VLo1vPEGDBoEnp5yf58xQYTa/LnoMgqwgQEWez48ra+wc9BoCuYo6ND1q7Sj4VcppdTDsdulIeqCBbLDe/astCFbtw4KFpTdXw+daKXSnmnKROvRo2HDBunL++GH8ntXwYLxjzHZEbKDvddHct1lCxj2m3cYeNg6YYnuqsE3i9Pwq5RS6sHFH0gDKZ789Vdwc4NHHoEnn5SRwvE0+Ko0ZrfLDu/o0TLhulgxaV322msJf/2i4qL4ef/PTA6YTODZQKxOOcgR14gb1h2ADQMn3O01KJ7X/Z7fS2V+Gn6VUkolT3xJw4IF0pps507ZWnvlFZm61qmTBl2VrmJipJb3iy/g8GEoVw6mTYPeveV3MYBTV0/xTeA3TN81nUs3LlG5YGX8OvqRjzZ8+ts/hMXsJ8qyDzd7dfJaqzGkg7djX5RKcxp+lVJK3dupU3JAbeFCOH1a+u8+8ghcvizh9/ZdXqXSwY0b0rVh7Fj5K1mjBsybJ286BPzrz1cBG/Fw9WDjPxv5NehXALp4d+HN+m/SyqsVxs13LdydcuG7xoWzYZUpnle6PXStXcKRL02lA+32oJRSKjHThF27ZLxVjRrSpaFGDejQAXr0kCkAOmVNOcCVK9Kfd+JEuHQJmjaF4cNl8J9hwPoT63l07qPE2mMByOOah/71+vO6z+t45vF08OpVetNuD0oppe7t0CHZPps/H44elfmu8+aBt7ckDS1pUDct3QRNwjMAACAASURBVB2C75ogzoZFptmO6e3fo6AlDyVDarJusQcRETLtevhwCb8AR0KPMCVgCt8EfnMr+Fqw8G6jd/mwxYepui6V+Wn4VUopBV26SB2vYcioq6FD5SBbPA2+6qalu6U/bmSsDYCQsEiGL94HkGoBOP57XLvgyrUd1QjeV5JAu4Wm7W8w+Ysc1KwJdtPOiiOrmBwwmdXHVuNscaaVVys2n9xMnD0OF6sLbcu2TZX1qKxFw69SSmU3587JobVVqyTwOjtL3W6bNvDUU9qHV92T75qgW8E3XmSsDd81QakWfj+efYZTv1fnxuHiYLGTq9oZcjc4gb2MSWnvOoz3n4VfgB/HrxynWK5ifNLyE/rW7UvRXEXxP+2v09fUPWn4VUqp7ODqVTmwNm8ebNwonRtq1oSQEPDygldfdfQKVSZxNizyga4/iD/+gDFjYO+KBhguceSudwIPn39w8ogmxgjm7+vLKfHVFm7E3qBJqSaMaj2K7pW742x1vvU1GpVqpKFX3ZOGX6WUyqqio+H6dcifH/bskZZk5cvD++9LPW+VKo5eocqEiud1JySJoJvS/rimKW9CjB4t4bdgQSjV7gRmlaPE5thLmHUVMcYpYq0nsODCi1Wf4836b1K7WO2HfSkqm7I4egFKKaVSkd0OW7ZA375QtCh8+qlcb9ZMuv8fOSLXNPiqFBrSwRt3Z2uia+7O1gfuj2uzydnK2rXhscfg5Enp4hAcDJ99ep2I3JM47/Ie1502EWs5QV5bJ2Y9spOZXWZq8FUPRXd+lVIqq/j8c+nwf+oU5MwpB9a6dZP7LBaoW9ex61NZQnxdb0q7PURHw/ffw5dfwvHjUKkSzJoFzz4LBy/vYcC6Sfy0/yeirFEJTzIsdKpWjd4NqqXFS1LZjIZfpZTKrEJCYMUKqdc1DNkyq1ZNiiYff1wCsFJpoGvtEg98uC08HL79FsaNg3//BR8fWLQIHuscy69HltBm7iT+OPUHOZxz8GLNF2nq2ZRXl71KjC0GF6sL/Rs9nkavRmU3Gn6VUioziYiAxYthzhzYsEEKJps1g8qVZdf35uQqpTKKS5fg669h8mQZUtG6tfz1rd7wAjN2TWfApKmEhIdQNl9ZxrUfx0u1XiKfez4AyuYrq50bVKrTCW9KKZVZbN8ObdvKIbayZeG55+RWoYKjV6bUHU6fll3e6dNlHHHXrtDpdX922X7knyv/sP6f9cTYYmhXth0DGgzg0fKPYrVY7/+FlUomnfCmlFKZTVAQ/PADlCkDffrIiOHnn5fiyKZNdZdX3VV6TGC7m6Ag+OIL+PFHeWPi2Wdh0JBYfg0dTd/Nn2A37QA8UfkJPmv9GZUKVkqXdSkVT8OvUkplJJcvw88/y4mgv/6Sg2r9+sl9OXLA1KmOXZ/K8NJjAltSdu6UdmWLF4ObG7z+OvTuf55V57+l49pvOBt+9tZjrYaVusXqavBVDqGtzpRSytHs9oQ/9+4N/ftLaYOvL5w5I8WSSiXTvSawpTbTlNLzdu3kANu6dfC//8GvATu52voFmiz05KNNH1G9cHXGthuLu5M7VsOKi9WFll4tU309SiWH7vwqpZSjHDokPZ7mzoUdO6BECfjkExg5EmrV0rIGlSJpOYEtnt0uk7HHjJE3KIoUgc/HxFKk9SJm7pvEqIXbyOWSi751+vJm/TfxLig9gBuXaqwH2JTDafhVSqn0FBEBP/0koXf7dnBykg7/ERFyv/biVQ8ptSew3S42ViZkf/EFHLzmT77am+g/oQaFqu9m8p6pnF15lnL5yjGhwwRerPUiedzyJHq+jh5WGYGGX6WUSmt2u9TyFiwI165JDW/lynIU/rnnoHBhR69QZSFDOngnqvmFlE1gu11kJMycCWPHyiS2si38cX66FVfMaKaEAVuhQ7kOfNvpWx6t8CgWQ6sqVcal4VcppdLKyZNycG3WLGlHtnYtFC8u5Q4VKmhZg0qWB+3c8LAT2G4XFgZTpsCECXDxIjRsEsvTny7h58tDib0aDYCBwdsN3+arDl+l7AUqlc40/CqlVGpbvVp2ddevlxNBbdrASy8l3F+xouPWpjKVlHZuSMkEttudPw/jx0tzkWvXoHWnS5TrMZ1Vl6bw5T9nKO5RHGeLM3bTjovVhaeqPJXi76VUetPwq5RSqeHAAenHmyOH/PnoURgxAl54Aby8HL06lUndq3NDWrQt++cfaTLy3XdS39u2117cW33N6pC5bDgRTduybZnScQodK3RkR8gOPbymwGaTedWnT8vt1KnEHz08YONGR68yEQ2/SimVUhERsGCBjLDavl1KHHr3hjffhHfekR69Sj2E9OjcADDxl/N89rnJpb2FwRpHjecW4lb/W9Ze2EyOf3PwUq2XeLP+m1QtXPXWc/TwWjZx/bqUcAUHy+3UqcQB9+xZiItL/BwPD/D0lFvlyo5Y9T1p+FVKqQcVFQUDB0rXhogI+T/3r76Cjh3lfldXx65PZRlp2bkBwN8fBg6LJGBLEai4Amu/LzEL7ONv6xUKhZfEt50vfWr3IZ97vlT5fioDioxMCLa33/75Rz5evJj48S4uULIklCoFLVrIR0/PxB/z5Pnvd8lQNPwqpVRyXLkCu3ZJ/a6rK+zbB08+Ca++Co0a6eE1lSbSonODacrZy9GjYfNmsJYKwqn/YOIKrcdmAKZB3pjelHd5jsGN26XCq1AOZZpSxH38+J23f/6R+27n4gKlS0u5Vrdu8vH2W5Eimf5dLQ2/Sil1N6YJf/4J06bBwoVgtcK5c5Arl1zXwKvSWGp2brDZYNEiGUyxe4+d/A1XUvGziRyJWwfm7WHGAAP+vRqTSq9Cpbm4ONmlTSrgnjgBN24kPNYwZHe2XDno1EnOKsQH2zJloGjRTB9u70fDr1JKJWXzZnjjDTm8ljs3vPwy9OkjwRc0+Kp087CdG6Kj4Ycf4Msv4eipaxRqP4vCT0ziQtxx3N1L4HnjZSJveHLJ5XNMMw4DJ9zs1VOttEKlErtd6muPHJHb0aMJfz5xInHdrZsblC0rAbdtW/kY/7mXV7YvzdLwq5RSILu8O3bIQY0qVSB/fnB3l87+Tz8NOXM6eoVKPZCICPj2WylHD4k8RqHHJuH+/Cwu2sNpVLQRXzcYRffK3Vnx9wWGL96HNWYUUZZ9uNmrk9da7aFKK9RDCAuDw4fldnvIPXpU6nPjublJv/Bq1aB7d/lz+fIScIsVy/K7tw9Dw69SKnsLD4e5c+Gbb2DvXmlNNns2VK8OAQGOXp1SDyw0FCZNgq9+2UZ4+RnkfuoIRt5thFmc6FG1BwMbDKReiXq3Hp9QWuHC2bDKD1VaoZLJboczZyTgHjqUEHYPH5bSqnhOTrJjW7GinDeoWFFCbsWKUKKEBtwU0vCrlMq+/vc/+PpraeVTs6Z09H/2WUevSqkUOXNGdnmnfXeDG/U+gad8wTC5BrxU6yVGtR5FMY9iST73YUsr1F3ExsKxY1I+9d+Qe3sdbr580jWmY0eoVEn+7O0tNbhOGtVSm/4vqpTKPm7cgKVL4ZlnZMfExQV69IDXXoP69bWOV2VKR45IPe/3S05jq+uH84DpYL18636rYaVC/gp3Db4qFcTFyeGyAwcS34KCJADHK11agm3z5gkht1IlKFRI//8nHWn4VUplfUePSlnDrFnSsqxoUWjdGj7+2NErU9nQ0t0hqdK9Yfdu+Hy0ycLt/lgaT8R8axGGxaRzpW60KdOGd9e+S4wtBherCy29Wqb+C8mO7HbpqvD334lD7uHDEHNbdwwvL6haVXZyq1aVW6VKenYgg9Dwq5TKus6fl4lra9fKW4dPPAH9+0OzZo5emcqmlu4OSdS3NyQskuGL9wEkKwCbJmzZAqPGxPB7yC9YmkyAPoF4uOTl1brv8Eb9N/DK6wVAraK1dPzwwwgLk37ef/+dcNu/X04SxvP0lGDbvn1CyK1cOaErjMqQDNM00+2b+fj4mIGBgen2/ZRS2dD58/JWY/Pm8lZk8+bw6KMyjKJoUUevTmVzTcZsSHJiW4m87vw5rPVdn/fnKX+mrdlEwG81OXxtF5YGU7Dn/JcK+bx5u9EAetfsTS4XDVwpEhcn/5/x99+Jw+7p0wmPyZcPatSQW/XqcqtaVbrDqAzLMIydpmn6/Pe67vwqpTK/+GEUU6bIMIoCBeQHl5MTbNvm6NUpdcvZJILvva7HxcFns/359FQrTEsM1DXBgNZlOvBu4+9oX649FkNP/CdbRIQE2z17pG5kzx4JvNHRcr+Tk5QnNGuWOOyWKKE1uVnIfcOvYRjfAZ2AC6ZpVrt5zRfoDMQAx4GXTNMMS8uFKqVUkn7/HQYPlh9ouXNDv35y0xPSKgMqntc9yZ3f/w6UiIqCmd/Z+PTnFVyo9Q7kl3BmYDCwwUDGPzI+XdabqZ0/nxBw4z8ePSq/LIP08q5dG958U7q91KwpwdfFxbHrVmkuOT8dZgOTgTm3XfsdGG6aZpxhGF8Aw4H3Un95SimVhBMnpMF78eIyctgwpJv/s8/qgRKVoQ3p4J2o5hfA3dl6a6DEtWswYUo4Y9fPIrzK19D6OHmshblhOmM37bhYXehRtYejlp8xmab0edu5M+G2e3fifrleXhJ0e/WSj7VqQcmSupubTd03/JqmucUwDK//XFt726fbgSdTd1lKKfUfpgnr1kn3/uXLYcAAmDABWrWSH3T6Q0xlAgkDJRJ3e2hcogT93z/BzP2TiKnyHTS9RtXcjfmo/Wi6V+5GQEiAHl6DxEE3MDAh7F68KPdbrTKhsUMHCbi1a8uObt68jl23ylCSdeDtZvhdHl/28J/7lgE/m6b54/2+jh54U0qlyIwZMG6ctBMqXBj69oXXX5c6PKUc6GHblgUHmwyauIVfL0zAXuFXDKx0KNGDTzoOpH6J+mm48kwiJEQmLd4r6Pr4QN26cqtZU8aSK0UaHXgzDON9IA6Ye4/H9AX6Anh6ej7Mt1NKZScnT0pDeJDDbB4eMGeODKVwdXXs2pTi4dqWfb9pMyNWTOZk9F4ocBRXjwK8WHU4H3boT4nc2fSXurAwCbk7dkjg3bEDzp6V+6xW6a7QqZMGXfXQUrzzaxjGC8DrQBvTNG/c5amJ6M6vUuqe7HZYvVpKG1avlh+APj5y+sfNzdGrUyqRlLQtW7nlPP1/eZ+TBWaCAZgG/Wu8x9jOH+HunI2CXHQ07N0rATf+FhSUcH/FijJ1sX59qFdPg65KkVTd+TUM4xHkgFuL5AZfpZS6q8hI+O47mDhRTmMXLSrT1+LfLdLgqzKg5LYtM02Y9utePlk9kXOF5kLBGLi572S1WChZKHfWDr6mKYdUt28Hf38Junv2JIz9LVoUGjSA55+XsOvjI311lUojyWl1Ng9oCRQ0DOMMMALp7uAK/G7IIZPtpmm+nobrVEplRTEx0lbIZoP335fdnrlz4ckntd2QcogHqeG9X9uyOJudD+eswG/XeMILbsQomIPG7q/wartm9F/7ctYdPRwRIeUL/v5y2749oU43Vy4Jt4MGJezsag9dlc6S0+2hZxKXZ6bBWpRS2YFpyuCJCRPkbc69e+UH4r592npIOdSD1vDerW3Z681L8tzXk1lwaiKxHsdwci9Jd48vmPzyKxTLmx8A76Kls0b3BtOEY8cSQq6/v/Tcttvlfm9veOwxaNgQGjWSul2r1bFrVtmedoFXSqWPmBj45RcJvYGB8rbmq69KPa+7O5Qq5egVqmzOd01QoiALEBlrw3dNUJLhN/7aR6sWc/pGIIWcy+EUc5EXlv+A6XqVHLEN6V/iM8b07o6bi3Oi5zYq1Shzht6oKOm48Oef8McfEnYvXZL7cueW8oX335eg26CBDJJQKoPR8KuUSh9Ll8Jzz8kEpalTpb5PB1KoDORBRw8DFC5wkiP2IUQ7RxNmAk4WCl18iuGt3ubtJxtm/jcyLl2Sd2riw25goPwiC1Km1LkzNG4sO7uVK+uursoUNPwqpdLGgQMwfrz04Rw0CLp1gzVroG1bsFgcvTql7pDc0cMAsbZYZmxbxNB1Q4i2RN+8atCr7GB+/OSLNF5pGjFNOXD6558JYTe+A4OLi7QXGzgQmjSRwFuokGPXq1QKafhVSqWe+ClsX30lrcrc3GDIELnP2Rnat3fs+pS6h/uNHga4EnmF0WunMyVwEtedzkB4SYzczlgMOy7OLrzRuqsjlp4ydrvU2m/ZknC7cEHuy59fAu5LL0nY9fHRrisqy9Dwq5RKPW++CVOmQJEiMHKkTGErWNDRq1IqWe42erhr7RIcDT3Khysmsuj4bOIs1zHOtObRPFP5+q2OXHT5K3McXouNhV27EoLuH3/IYAmQgTIdOkCzZtC0qRxU03doVBaVrCEXqUWHXCiVxYSGwrRp0KuX/PDctk3eJu3ZU3eJVKZnmiabgjfxwcrxbLu4HGzOOB1+lufLv83nA2tStKijV3gfkZHSUzc+7G7bBjdutub39obmzeXWrFnCNEWlspA0GW+slMqmjh6Vrg2zZ8sP0/z5ZZe3cWO5KZWJbQ7ezKQdk9l2fA//xhyD64VwP/AhAxv3473pRcmb19ErvIuoKPjrL9i4ETZtkk4MMTHSPrBGDejTJyHsFini6NUq5TAafpVSyWe3w9NPw6JFUsPbq5ccZqtW7f7PVSqDu3j9IsPX/Y+Ze2bIBdPAY+9wPm71Ea9/7EaOHI5d3x1iYmRnd+NGufn7SwA2DKhTB956C1q0kDIGnZim1C0afpVS92azSW1gixZSA1i8uPTxfOMNMv77vkrd38GLBxn35wR+2PsDsUTJ6GEDLIaFoW97MKhFBinhiY2FgADZ1d24UToyREZK2K1ZE/r1g5YtZXc3w25PK+V4Gn6VUkmLjJSyhq++kglOe/fKW6cTJzp6ZUoly71GFZumyboT6/jyj69YF7wa4txgT2/KWltypk4fbMjo4TZlWzruBdjtMi1t/XrporJ1K1y/LvdVry5DYlq1krCrwySUSjYNv0qpxK5dk3reSZOkwX39+jKZrWpVR69MqWS726jiGFsUEdbNjP1jPIcu78e4XgT+GknzHK/z0eCCtG4N2894Oa57w4kTCWF3w4aE6WmVKsELL0Dr1hJ2tceuUimm3R6UUiImRhrZX70Knp7yA3bIEDkck+nHVKnspsmYDYSERRJtOUSUZR/O9jLEWI5y3WkVccYVjAs1MP8cxOPlnuGDYa7Uq+eghV68KCF33ToJvf/8I9eLF5eBMG3ayK3EneOVlVL3pt0elFJJCwgAX184flxGl+bJI7tPBQo4emVKJXKvMob/Onsz+J5z+R8QKxcN4EhHLDve5bkmrRg2z6By5XRbvoiMlPKFtWsl8O7dK9fz5JEShkGDJPR6e+svnUqlEQ2/SmVHpikT2L78Ug7P5M4th2Wio6U/rwZflcHcrYwBuCMAm6ZJjtwHOBM9EYybwdc0YPtAil0bxPb1pfD0TKeFm6ZMUVu7Vm5btsi/MxcXmZw2apTs7NatC076I1mp9KD/0pTKjhYuhB49oGRJGDtWDs7kzu3oVals5kF2cn3XBCUaOwwQGWvDd03QredExUUx9++5jN8+noOxByAmn7TkM+xgOlOiWlkmP2tJ++B7/rzs6sYH3nPn5HqVKtC/v4z5bt6cjNc7TansQcOvUtlBZCTMmgW5ckHv3tClC/zwgwRgFxdHr05lQw+ykwtSxpCUs2GRXLh+gakBU/EL8OPijYu4X60JG74n779P4tZ8BZHlVuKZuw6fdut+13D9UKKjpR1gfNjds0euFygA7dpJ2G3XTn7ZVEo5nB54Uyoru3IFpkyR9mQXL8KTT0rnBqUcLP5A2n+VyOvOn8NaJ+vxMcYpbDmWc9VYT7QtmlxnHyPi90GUMVrx3lCDF15IwynbwcGwapXcNmyQFmTOzlLK0L693GrXlt7YSimH0ANvSmU306fL4ZmICHj0URg2TDo3KJUB3GsnNylDOngzfPE+rtj2E2FdQ6xxmhhrEFabG+5BLxK97m3KFq/EsFHw1FNpUD4bFSX1uvGBNyhIrpcpIy3IHnlEBkx4eKTyN1ZKpTYNv0plJUeOyGSnwoWlXdnjj8PQoTL9SakMpHhe9yR3fovndU/y8Y9WL8jcQ6tYeHQSYIIJzn+/Quzq0dSsXZDhs6Fjx1RukHDiRELY3bgRbtwAV1eZdvj66/JLZcWK2pVBqUxGw69SWcGuXTB6NCxaJGF3zBjo0EFuSmVA8Tu5tx9ic3e2MqSDd6LHhd4I5ZvAb5gcMJlzEedujR7GbqVMvrLMWF0w9d7QiImR3d3lyyXwHjki18uVg5dflt3dVq30oJpSmZyGX6Uys61bpVXSmjXSJ3TYMBg40NGrUuq+4g+e3a3bw9HQo4zfPp7Ze2YTGRdJqegOWDcPw9ZqOIZTDC7OLswe0ZJGpR5yIefPw8qVEnh//x3Cw2V3t1UreOMN2d2tUOEhv4lSKiPR8KtUZmOaCW+zTp6csOvbr58EYKUyia61SyTqvmCaJltObmGc/ziWBS3DyeJMqSvP8c9P73D+SjX6vAhtH6vPsbiHGD1st8Pu3bBihQTegAC5XqIE9OwJnTrJCOGcOVPlNSqlMh7t9qBUZmGzweLFUtIwZw5UrSr9Q3Pn1rdhVablf9qf9f+sx2a3sfzocgLPBpLHuQCFT/bn6Nw3yGUU4fXX4Z13ZOJvikRESN/dFSvk9u+/8gtkgwYSdh97TOritXZXqSxFuz0olVnFxsLcuRJ6g4LkgM2lS3Jf0aKOXZtSNz3IwIp4vx//ncd+eoxYu0xhK+hSivKHv+HYoudxyp2DT9+TyoP8+VOwoJAQWLYMfvsN1q+Xet7cuaUOvlMnKWcoVCgFX1gpldlp+FUqI7PZZEfq0CGoVQsWLIDu3cFqdfTKlLrlQQdWnAw7ycS/JuIX4Hcr+GJauLTqddz+eY0JvvDKKw9YeWCa8PffEnZ/+w3i32UsV04SdOfO0LSp9OJVSmVrGn6VymjCw6VrwwsvSMh96y3w8pKT5vq2rMqAkjN6GCAgJIBx/uNYeHAhABWdW3M4ZgsmcRimC8OfbcWIlx9g6GB8d4b4wHvypPwbadhQ6uAffxwqV9Z/N0qpRDT8KpVRhIXBpEkwYQJcvgzVqoGPjxxkUyoDu9fACrtpZ1nQMsb5j2Prqa3kdslNU6d3ODx7AIeOlsK7rT91n9hEv0db0rR0Mg6whYVJG7Jff5WP166Bu7uMD/7wQ6nf1XIgpdQ9aPhVytHCw+GLLyT4Xrsmu1Xvvy/BV6lMIKmBFXaicPbYSqXJ73D08lFKenjS3v4VOyb2YfP53LRqBT/4Qdu2jTCM+4Tef/+VsLtkiYwSjouTQS5PPSX/Xtq21UOfSqlk0/CrlKPExckMVicnmDVLDuL8739S26tUJhI/sCLMtp8blh3YjMtEWQOwxV2jpqUej16fzybfJ1h73YkuXWD4cGm0cE/HjknYXbIEtm+Xmt7y5aXtQ9eu8gW09l0plQIafpVKb6dPg6+vnEDfs0fesj18GDw8HL0ypVKka+0SbDmzhPG73gPTDgaUy9EA7398WTuqKfvtBr16yfDBqlXv8kVMU/49xAfe/fvleu3a8Mkn0K2bPFnrd5VSD0nDr1Lp5cQJKW+YNUt+0PfuDdevQ968GnxVpmSaJpuCNzHWfywrj66UiwZgWjm+rAshAc14/TUYPBhKl07iC9hs8OefCYH35EmwWKQrw/jxssPr5ZWOr0gplR1o+FUqPezcmfA27SuvyBaY/lBXmVSsLZZfDv7C2G1j2X1uN4VyFOKxoq+wKuRH7MSC3YXezVviu1BKcxOJi4PNm2HhQhnacuGCtHeIP7D2+OPaf1cplaY0/CqVVo4ckf68XbrIW7effirty0rcu/G/UhnVtehrTN85nYl/TeT0tdNUKlCJN0pNZ8eM51jh70a+Gi9Tv8cm3u3eknaVbzvEFhMjB9UWLoSlSyE0VA6oPfYYPPEEdOyo734opdKNjjdWKrUdOgSffQbz50OxYhAcLIfalMpk/E/7syl4E5ULVuaP038wfdd0rkVfo4VnS+pEv8uayR05eMBCmTIwZAi8+KKUsAMQFQW//y49q3/9VVqUeXjIdLUnn5S+1dqhQSmVhnS8sVJp7dgxedv2558lAbz7rhQ7avBVmZD/aX9afd+KaFs0ABYsPFn5abz+fZefP6nL5pPSivrHH+Hpp2/+NY+MhCWrZYd32TJp45cnj7z78eSTUtrg5ubYF6aUyvb0p7JSD8tul0M6oaGwfDm89x4MGqR1iypTMk2T1cdW8+aqN28FXwODJpZ32fTWl1y4AI0aSVvqxx4DS0wUrFgjv/T99psc4syfH3r0kJKGNm0eYGSbUkqlPQ2/SqXU3r0wciQUKADTpsmBtpAQyJ3b0StT6oFFx0Uzb/88xm4by4GLByiYoyBOFmdsNjtmnAtbv+9Gh9rSo7d5wxiM39fCiwukpOHaNQm8zz4rgydatgRnZ0e/JKWUSpKGX6Ue1J490nd06VIJuoMHJ9ynwVdlMlcirzBt5zS+/utr/o34l+qFqzO26fccXfwMs1bvJK7EJlqXaYnvQh/qXFkNs3+GrkulhjdvXiln6NEDWrfWwKuUyhQ0/Cr1IKZOhf79pY5xxAgYOBDy5XP0qlQ2s3R3CL5rgjgbFknxvO4M6eBN19oP1kUkOCyYCdsnMGPXDK7HXqdd2XZ8XGs2W2a34723DCwWeOn5enzULIISf34HHTvB5cvyC17XrlLo27atljQopTIdDb9K3c++fdKft0oVePRRCb1vvy27Xkqls6W7Qxi+eB+RsTYAQsIiGb54H8A9A3B854YiuYqw9vhafjn4CxbDQs9qPWmX810W+tXktd8gVw47E578k9bXZ1Ds51/J910YN1zcCW3zCKX6vQTt24Ora7q8VqWUSgsafpW6m/37pbxhONIDNwAAIABJREFU4UIZrbp4sQym+PhjR69MZWO+a4JuBd94kbE2fNcE3TX8/nnqT1rPaU2MLQaAHM45GNTwXWpFD2DGVyXpvcmkRe7d/NF4Hg1Pzsf682kinVxZX74+yyo3Y1OZulhy5GB0yep01eCrlMrkNPwq9V8HDshAil9+gVy54IMP4J13HL0qpQA4GxaZ7OvRcdHM3TeX4euH3wq+BgYd8wxmw/BP+HXXEV7z+IRFheeR/0IQ7HCCDh34pMlzLChem+uut/XhvU/AVkqpzELDr1L/tWgRrFwpx9oHDZJuDkplEMXzuhOSRNAtntf91p//e4itXL5yhEWGEWe3QZwzTmMuM/uiD9XZiRlhYNRpDp++I4fXChRg9rAVJDX+6G7BWymlMhMNv0odPiwty554Arp3l13eN97Q0KsypCEdvBPV/AK4O1sZ0sGbk2En5RDb7hlExETQrmw7vnlkNmeW+rBt3hfkKLCAF4KDaXJmMqaPD/Qch/H003eM3E5OwFZKqcxKw6/Kvo4dk/KGuXNlIlvjxnLdw8Ox61LqHuLLDm7v9vBEw1gWnBjCk8sWYBgGz1R7hlfKv8W5z85g6fMNfWJX0J8YIkp6k/PVT6HnMxgVKtz1e9wrYCulVGan4VdlTyNGwKhR0pf0nXdg6FAoXNjRq1IqWbrWLkHhAieZvusH9l/Yz7ubA8nlkouBDQbw/OXGRI5YS6V9HWhBGFdcixD69BsUH9KLXHXqgGEk6+sDD91OTSmlMiLDNJOq7EobPj4+ZmBgYLp9P6USOXNGShnc3WW3d8cOGDYMihVz9MqUSrZYWywjt4zksy2fYd6szO1f9hkGbytG7h8XUSDiFBHkZHeZ7pR87znK9GkNTrrPoZTKfgzD2Gmaps9/r+v/I6qs79w5GD1aRhB//rkcYuvVS25KZRLXoq8xfed0Jvw1gTPXzty6brVDyW/nU+oPK+ss7TnbZjStJ3ahWdWctx6TGkMxlFIqq9Dwq7Kuixfhyy/Bzw9iYuDFF+VAm1KZyNnws0zcPpFpO6dxNfoqLd2rMDC4Ih+VOEKMBZzsBvvPDmTMm8N45f0iPFI08fNTOhRDKaWyqvuGX8MwvgM6ARdM06x281p+4GfACwgGepimeSXtlqlUCvTqBevXw3PPwYcfQvnyjl6RUsl28OJBxm4by49//4jNbuPJG6UZvDiWescOcs7FkyOFX2B15fw81uIpRu9sdNeBgykZiqGUUllZcnZ+ZwOTgTm3XRsGrDdNc4xhGMNufv5e6i9PqQdw7RpMmACvvQZFioCvr4xhrVTJ0StT6r78T/uzMXgjeV3zsvLYSlYcXYG76UTfIHcGrQmn1I0LrMjRg6H0JrhoCwYPtTDxZSlhv5cHGYqhlFLZwX3Dr2maWwzD8PrP5S5Ay5t//h74P3t3HV9l/QVw/PPsLhgwukd3ScgEB4KTFgUmISAi8SNFRRSEEYKS0o2UpKSUhKKU1OhuaRg5YAzW997n98chREJgcRfn/XrdF9vdje/V7d7znOd8z9mIBr/KUUJCpLThhx/g1i3Inh1atYKSJR29MqVeyJYLW6gyq8rDKWypoyx8twU67LISkbMcI50/4ceIeuTOl4Luo6FxY2lU8iK0Z69SSj3O6RXvl9k0zSsA9/99Zo8owzDaGoax2zCM3Tdu3HjFp1PqGcaPh3z5oFs3KFcOdu+WwFepBCAsKoyJO8dTd1ath4Gvkx2+PJqGt90HUjPFBXIc/5OteZsxd1kKDh2CZs1ePPAF6dnr7mJ57Drt2auUSspifcObaZqTgckgrc5i+/lUEmC3g9P947bNm6Ws4Zdf4K23HLsupV7QzdCbTFjTj7EHpnDDCKXIdbiXHmxOBk64MvzYr3x3vDzVqsF6P/DxeaH2vE+lPXuVUupxrxr8XjMMI6tpmlcMw8gKXI/JRSn1VDYbzJ8vU9kWL4bixWH6dEiW7NUjA6Xi0LmAI4xY8CXTgtYTarHz3t/wTVRZCvt8zRcnMrP40DaiTvlQp6w33WeD1xPdKV+Nb2lPDXaVUuq+Vw1+fwWaA4Pv/7s8xlak1L+ZJixdCt9+C0eOQIkSUucL/73bRylHM032/T6doRsGsND9DE4mND3nQZdCrXBp3ZUBMzyZ20lu2qzZ23wzQfdoKqVUbHqRVmfzkM1tGQzDuAT0QYLehYZh/A+4ADSMzUWqJMw05Zzvpk1QqBAsWAANGjwqe1AqHvK/6M+Gg8tJse8wKy+tZ23WMDycofO94nR69zsuv/cBvQYbLOsFyZNDzYb3uJLzABsI4sQyLUtQSqnY9CLdHpo840dVYngtSj2ybx+UKiXlDHXrQsuW0q9Xx7Sq+Cwqis3zh1D1VG8iDRMMSJ/emR/SNKRt05Hs3u9J80Gwfj2kSXO//bTPFQZtOKBDKJRSKo5o+kzFL/v2wXvvweuvw5o1ct1XX8l0Ng18VXx18iQh3Toz9t101DnSi0gnCXydMPii2rfkz7GQ6rU8qVYNjh2TFtQXLkj5+qSdx545hEIppVTM02hCxQ8nTkhN78KFkDat9OytVMnRq1Lq2UJDYfFiAmdOZJzNn3Fl4WZFeM09N6GRl7HZbVhw5adeVTm/VTryTZoEzZvL7JUHdAiFUkrFLQ1+leNZrVCtmgyo6NULvv6aZ85qVcrR9u2DqVM5u2I2w1+7y0/eBmHOUDd3Dbq+05vSGSrQ80d/pq/fyJ0DPqRJ780P86VU3WJ58uF0CIVSSsUtDX6VY1y/DhMmQM+e0rF/7lwoWBAyPXNeilIOsWxfABOW76Gs/xo+PvwnoebfDK3oxMJWJhYnZ5qVbEaX8l3J4lyE8ePBdzQEBnpTsaI3fpOhZs3nd+LrWqMQfksOPVb6oEMolFIq9mjwq+JWUBAMGwajRkF4uHRy8PHRARUq/jFN/pq1AvuIcfjd2cDsklYa1nXlcGZIbklOl3Kf0unNThj3sjFyKEycCPfuScm6nx9UqPBiT6NDKJRSCcatW3DxorQbDQ2Vf0NC4IMPpPXoX3/B2rWP/3zMGClnjEc0+FVxIyoKRo6EwYPh9m1o1Ai++07alykVn9y5Az//DJMmUeHwQfq+bWHgBzZMAzAjSWmtRVHXdrTNW4fvusCMGVK506gRdO8ubahflg6hUEo9sGxfQOwdDIeGwrlzEBgIN27AzZuSlGraFDw9Yd06GDFCrvvnZdcuKFoUZs2Czp2ffNxz5yBXLti6FQYMgBQp5JI8uQTAGvyqJMlikdKGN9+UP4zSpR29IqUeMU15c580CebPJywylOl1ctCzViqCkgXDw8HsThh3PTmyuAwFe0kDkpYtoWtX2dCmlFLRsWxfwGNlUP/Z+jAkBC5fhvTpIV06OH8epk6VwPafl/HjoXJl+PNP8PV98nG8vCT4jYiAa9dk303WrJA6tXydKhXL9gUw91J60vn2wD1tKupVKEDF0rklwM2WTR6nWzc59RXPp65q8Ktih90unRtGjIDff5c/yk2bIFUqR69MqUeCgyXLO3ky7N/PrfTJGd+uCGMznuFG5EVSGkVIE/EBd1zmYZpWsLlyd3FLnK5n5uuvJQGSNevjDxmrWRulVKI2dM0JwqJsuFkjyXz3Jpnv3STzvVssCb6Ib+k20iOxRQu4ckWC3uBgueOUKdC6tQS6AwdChgxyyZgRiheXLCxA2bIwb55cnyGDBM1p00oAC1Crllz+5WFQ7pwOCpUH4PdzFga9nhnfwv94f3vart54SINfFbNMU/rz+vnB/v3yR3f5sgS/Gviq+MA0YfduCXjnzYOQEC56F2XEkMpMidpBSNQe3sv1Ht9U+IbAwLx0GnaVkJNNiEq+F+NKedLnS8/IOdf5+O1sTzz0S2dtlFJJj2lKucGFC48uhQpBjRrcvXqDXVPbkzEk6LG7jKrQBGgjdbWRkfLZWr26ZFyzZn20yeD116XM8FlTULNmhcaNX3rJD4Lyf3rQjzwhvrdp8KtiTlgYvPuuFLznyQOzZ0OTJgnmSFAlcqGhEuxOnAh79kDy5Bz+pAZDy0Qw98ofEH6SJsWb0LV8V4pmeI1Fi6RE/dKB7LimDiedV3YKfHKD7rVT4Fv6ycAXEt8HhFLqFd29C6dPP7pkyyZTSk1Tsq43bz5++9atoUYNPDJn4M/8b3LFIz2XU2Xiqkd6rqdICzlz8iXIfbdsefbzPivojabE1o9cg18VfVevQpYsckRaoAA0bAht2oCrq6NXphQcPw4//ggzZ0JQENsq5WHmoPIcTmdj25WlJL+enI5vdOQr76/InCwnM2fCB0Pk86pwYZg+HT76KBmursX+86kS2weEUuoZTFM++x4Et05O0KyZ/KxsWdlD8E9160rwaxhSL5UyJeTMKZvEcuaU8gOga83C+IV2eqL14aA6r8XVK3uqxNaPXINf9erOnZOpbAsWwJEjkD+/1B0p5WhRUbB8uWR5168HFxfsDeoz9N3U9Dg7BXvEWbgCrUu3ZnDVwbja0jNpkpSoX7kiez8WL5Z9IS+TSElsHxBKJXn37sHJkxLoPqiF/d//YP58OZv0QKFCj4LfJk2gXj3ZBfvgkjr1o9v27PnMp4uvrQ8TWz9yDX7Vy7t+XTo2TJwoJQ2dOklNr1Jx7N+by3qX8qCm/0o5CLtyBXLmJHLA98ytkIohhydx7Myxh/e1GBYyu+Vl9OD0jBsnHfgqV5ZOPlWqvNpm5cT2AaFUkmC3S+/aXLnk+xkzYM4cOHECLl2S69zdJRB2coKSJWUPS/78kDevBLe5cz96vKe1AnsJ8bH1YXwNyl+VYZrmf98qhnh5eZm7d++Os+dTsSAkRN4ggoKgVSvo00faoygVxx5sLguPjKL8+YM027eKqn/vwIKJUbMmd9u1YEqGC4zcOZpLwZcokbkE9QrX44etPxBpi8Swu+I0Zx2Rp73x9ZU9mmXLxsy6EssHhFKJ0rFjsGqVnLE8fBiOHpUs7o0b0gFh2DD45RfJ5v7zUrRorNXUqthhGMYe0zS9nrheg1/1nyIiYPVqmeACUgTp7S0FkUo5SPW+v1J+yyqa7VtFvlsB3HRPxcIS1VldyYdSNa8yftd4boffxie3D90qdKNGvhqcOGHw9Sh/fj+2EeO8D03f9qZbN/lMU0olEqYpZ36OHHkU4B45Iv1vixWDadNkg1mWLNI1oVgxeRNo3Fi7EiUyGvyql2e3y+74Xr2kvnfPHmmjopQjHTkC48cTMnU6KaLCmV4qF1O9cnAmXXnuuR4hxPInGFH4FvalW4VulMtejt27YdAgWLoU3Nzkc69Ll0dnOZVSCVREhGRuDxyAcuWgSBEZ5FC9+qPbZMggAe6wYVLQf/eu7AvQcr1E71nBr9b8qieZJvz2m5wHPnhQprFNmqRT2ZTjWK3w668wbhxs2ABubmwo8Q6DyhRiW9YJwHlgC2Ahk6U6mzqMpGD6QmzYANVayqj51KmhfqtgzmU5wEprMHvnaUmCUgmK1SpjFW/elLra/fulhMFqlZ//8IMEv6VLw5gxj7K6mTI9/jgeHnG/dhWvaPCrnhQSAp98IiMN582DDz/UOiflGDduyKnKiRNlQ0rOnDB4MGarVhw4sIldmzuBESW3NSGdvR4Ta4/k2BZPPhkEO3dC5szymZjd+zL9/jioAyiUSgguXpRhNHv2wN69ktlt2BBGjZLgdfNmKVV4/33ZgFaqlGxAA8n0fv65Y9ev4jUNfpV40At1+HDpP7hunRxBa69edV+cbuTatUuyvPPnyzSjqlVh7Fjs79VixanVDP61DtsvbSe5ayqsVgumaeJkuFDF+X/0/sSTo0dlzsrEiTIJNFkyqDD4uA6gUCq+Mc1HgW5UFDRqJNeXLSvtxSwWCXKrVHk0xczVFc6eddyaVYKnwW9SFxAA330nGwBSpJAODiVKyJG0UvfFydjeyEhYtEhOV+7cKQdhbdvCp58SVTA/cw/N5YdJJTkWeIzcaXIzvtZ4WpZqyY7z+xmxdCO7FvqwaLc3r70GP/8sJyyc//EOpwMolIpHJk6UXtx79kBgoFxXuPCj4HfKFJlmVqKEtBlTKgZp8JtUhYVBv35yCslqlVNEPXvKm41S/xKrY3uvX5ezDhMnSqanUCEYOxY++YSQZBam7p3K8DE1uBh8kRKZSzC33lwaFmvIvWBnRg6FUaO8uXHDm/LlYcoKeO+9p/fo1QEUSsWx4GA5i7Njh1xOnpQNq05Osp/k6lWoUwfKlJFLiRKP7vv++45bt0r0NPhNakxTIgMXF1i2TKbQ9Osn54iVeoZYyZru2wejR0tdeWQk1KyJf6vqbMwcRums+dmxdyRjdo7hVtgtKuWqxKT3J1Ezf02uXTPo1UNi5eBgqFlT9mZWqvT8p9MBFErFIqsVDh2ScrlkyWRcYpcu8pkDULCgdGO4d0/aiU2Y8GqTZJSKARr8JhU2m5wLHjdO6nk9PKTGKnlyR69MJQAxljW1WuVU5+jRsmElRQrpO/b55/inuE3lWZWJOBqBiXxg1ilUh24VulE+R3nOnoWOHeGnn6Q0sEED6N79xZuQJLYJRUo5VFCQfJZs2yZZ3b175YziX3/JkWj58lJSV7asXNKmffz+GvgqB9LgN7F70Lase3c5Ki9TRk41eXho4KteWLSzprduSdeG8ePhwgUZBTp8uNSYp0nD8cDjfLbkM8Kt4QAYGHxR7gtG1RzF4cPwsZ/sfbNYoHlz6NoVChR4+dcRH8eGKhXv2Wzy+bFtG7zxhlyOHZMjUDc36f/erp1kdosVk/u8+aZclIqHNPhNzO7ehdq15Ug8Xz5YsEDerLRtmXpJr5w1PXZMsryzZklWyMdHvq9dGywWdgbsZPCawSw7vgwXiwvOTs6YpomrxZWiZiPq1IEVKyRB/OWX0tpTp2krFQfCw2UyzLZtsH27lCsA9O4twe/rr8v1pUtrVyCV4OiEt8QoKEh69AI0bSrtYVq31jcoFTdMU6ZKjBgBv/8umaGmTeGLL6BkSUzTZO2ZtQzeOpj1Z9eTJlkaPi/7OfncGzBs3V+cCziGfU9dQnfWIH16udtnn+kwJqVizaVLUoa0aZOM/O3TR/6Os2SBbNmkhKFCBfk3Vy4tWVAJhk54SwquXZMaq9mzZUdtzpxS56tUXAgPh7lzYeRIOHxYpkt8/z20bw8ZM2Kz21h69BcGbxnMnit7yOaRjWHVhtG2TFv+OBTMF4NuELj1EyKvpcaSMoxM1Y4xum9qGpfP5uhXplTi1LOnbDh90DPXwwOaNJGvDUNKlNzcHLc+pWKJBr+Jwd27Uj85bJjMOW/bVvsiqrhz/bq0XpgwQb4uUQKmT4cmTfC/vpe1R38kNCqUxccW8/etvymQrgBTak+hWYlmYHNj9mz4vLuF8JueOKe7R/p3D5CiWACGxWTsJncNfpWKDptNDkY3bZLs7vHjMi3NMKTLSunS0KkTVKwo/d0tlkf31cBXJVIa/CZ0ISHSGPzyZannHTjw1XYCKfWyjhyRLO+cOXLQVasWfPUVVK4MhsHaM2up9XMtouwyfrhg+oIsariIDwp/QFiohfFjpDIiIABcM0eRoe4ekhe8ivGPknQdQKHUS7LZJLB1cpJBEd98I6VwADlySCeGe/ckyzt0qGPXmkDF6bRLFSs0+E2ITFM2Gnh7y06grl3l63LlHL0yldiZJvzxhwS9a9ZIP88WLWQ3WuHCAASGBjJ2x1iGbBvyMPB1MpxoXrI572RuQL/vZYbFrVuy/+2nn6Dvnj1cvqMDKJR6aXa7dGLYsEEumzZJh58335SNzg0bSla3UiWp11XREifTLlWs0+A3odm0SY7kd+yQPr1lykjgoVRsioiQet7hwyXjmyUL9O8v7Y0yZADg4p2LDPcfzpS9UwiNCqVizorsDNiJ1W7FxcmVg7++Q6535WRFnToymOJBJ6TQjDqAQqkXYppSX+/uLuUMb78tR5IgwW6DBjIaHOQsTOXKjltrIhSr0y5VnNHgN6E4ckR69a5cCdmzS7qsVClHr0oldrdvy+jhMWOkP3SJEjBzJjRq9LAe8ETgCX7Y+gNzDs7BbtppWqIp35T/hmKZirFgmz/DFm1k3zIffrnozUcfybFb8eKPP40OoFDqGUwTTp+WgRLr1sHGjdK9Z+BACXZ9feUUio+PlDWoWBUr0y5VnNPgNyEIDZXTVnY7DB4svZ90Q5uKTefOwahRMpgiJASqVZNevVWrPmxztOfyHgZvHczio4txc3ajXZl2dCnfhVxpcrF3LzTsCIsXe+Pm5k37/8mk09y5n/2UOoBCqfvCw6WkyDRlQ9qBA3J99uwyz7tCBfne3R2mTXPcOuOh2K7HjbFpl8qhNPiNr+7ckWDjs89kEtvChfImmD69o1emErM9e6RryKJFEuQ2aQJffy27wIFtF7Yx/cB0Dl47yM6AnaR2S43fW350erMTGZNnYtMmaDtQyoJTpZKTFZ06SdczpdQzhIRIJ4a1a+XyYKKaYUDjxtLBp2pV2cysPXafKS7qcaM97VLFCxr8xjeRkXKa+fvv4eZN8PKSzWxVqzp6ZSqxsttlGMXQoXJK1cNDRql16iSZJsBu2hm2bRh+6/ywm3YAOnh1YFCVQXi4pmbVKhkG5e8PmTLJ1x06QOrUDnxdSsVXNtujlmL9+sklKkoGEb31lrzf2+3SsaF7d8euNQGJi3pcLdFKHDT4jS9MU7K7PXrAmTOySWHIENnQplRsiIiQISjDh8PRozI3eOhQaNPmYdRqtVtZcHgBg7YM4siNIw/vajEseKbMwYpfUvPDD7LvJnduGD8eWrbUqhylnnDxonRIWbNGsru7d0vNbqlScrBZtaqUMyRP7uiVJlhxVY+rJVoJnwa/8UVkpBzhe3hIm5oaNfT0loodwcFydmHUKLhyRUoaZs+GDz98OAI73BrOjP0zGLJ1CGeDzlIsYzH6VOrDkG1DiLRF4mS6Mq6rD1d3QbFicvdGjcDFxcGvTan4Zu9e+PhjOHZMvvf0hHr1JOEBULu2XFS0aT2uelEa/DrSkSOSdRs/XlJl69fLSOJ/TthRKqZcuQKjR8s0tuBgqFIFZsyQzWz3D7TuRtzlx90/MmL7CK7eu0o5z3KMqjmK9wu+z727TlzfXoPZWzZy77APubN7M2k5vP++nJ1VKkkzTanT/eMPye7Wqye1PzlyyKV1a0lqFC2qiY1YovW46kVp8OsIly9Dnz7SrszDQ3qllisHefI4emUqMfr7bylnmDkTrFaoX5+NdVvQ82IyLq8PI9veDbTzycipsEWM3TmWoPAgquWtxtx6c/HJ7cONGwa9esr04jt3vKle3Ru/2dJeVD/DVZJnmvIevmqVvLeDBLgPToNkzCjBsIp1Wo+rXpRhPjj1Ege8vLzM3bt3x9nzxTtRUbKRbfhwCUI6doRevbSDg4odu3ZJ3fjixVLO0KIFdOnCsrvu+C05RJDtMKFO27Ea1wiz7MI0IqhXpB5+b/nhlc2Lc+ek8cO0aVIeXL++VOZoGbpKskxTztitXg3Xr8sfCMjpj+TJJbNbvbr221UqnjAMY49pml7/vl4zv3HBNCVF5uws4yfr1pUG5ZrpVTHtwfjhH36Q37XUqSVi/eILmcoGDB28nkDbBgJdhwHSuSGZzYvC7p+y+MOWHDkCn3SXgW5OTtCsmQymKKRnDlVStWWLFLb/9ptsXAN4441HXRtWrnTs+pRSL0WD39hkmrBkCQwYIJmCLFlkl2+yZI5emYrHXqlJu80mvXl/+AH274ds2aTUoW1babh734GrB9gf0pdQ1y2ACQZgOpHMLMq14wXx9YXlyyWJ9fnn0uL3frezuHkdSjmaacLJk1LG0KoVpEkDW7fCvHlSH9+njwya8NTfZaUSKg1+Y8vWrdC1qzQ+LVpURsNmyaKBr3qul27SHhEhw1CGDIFTpyQ9O20aNG36cPwwwNYLWxm0ZRCr/l6FxTkFya3vEGbZgmlawXQh5PfGBO2pwKa08O23EvhmyBCHr0MpR4qKkiETK1ZIFvfUKbm+cGGoVUtK1Dp3ftgNRSmVsGnwG9OsVpnIs3gxZM0KU6ZIraWz/qdW/+2Fm7SHhMDkyVJzePmyFOIuXgy+vg9bL5imyR+n/2DgloFsOr+JDMkz0P+d/uRO9gH9ll/g6pmPuXf7LLZD72IGlaZF5zuM+S41Hh5x+DqUcpSbN+XvKGdO2RRapYocMFauLIFurVqP5nGnTOnQpSqlYpZGZDHl3j15g3R2lpTZ99/DV19BihSOXplKQP6zSfvt2zBunLQsu3kTfHykXVnVqg9bL9hNO0uPLWXgloHsvbIXTw9PRtUYRevXW+NCCubMgXtzC3DnnAvOaULI985FBnW7ScNyMReUxlWzeaVemGnC8eOS3V2xArZtkzMks2ZBkSJS5lCpkga6SiUBGvxGV0gIjBghHRw2bYISJWSAgFKv4FlN2l9zCpVdZxMnyoFW7drg5yejr+/bfH4zY3aMYdflXZy/c5786fIztfZUmpVsRlS4K1MmyK/ppUtQsqQLw+dDgwYpsFgKx9nr0GbzKk492GwM0oVh7Vr5unRp6NkTPvhAvjcMyfQqpZIEDX5flc0G06dLgeSVK/ImqhkDFU3/btKePegqHXcv5cNDa8EaJWPUuneXg6z7wqLC6L2hNyP8R2BiYmDwvc/39KjYgztBFgYPgDFjJFFcqZJU4sT2AEFtNq8c5u5d6au7fDns2CGju52dZYJh/frSluxVd3EqpRIFDX5fhd0O5cvDzp3w5puwcCG89ZajV6USgQf1sItmraHBn3OofewvDIszlpYtZANl/vwPb3s34i4Td09khP8IroVce3i9k+FEaIgz3b6xMGmSJIrff18SxeXLx+3r0G4PKs5s3Sqdddatk3Hx6dLJGZLgYPm6TRtHr1ApFU9Ea8iFYRidgdaACRwCWpqmGf6s2yf4IRfHjsnuX8OQusssWSSToGPY87l5AAAgAElEQVSuVEzZv18+wBcvlpHX7dtL7fg/2irdDL3JmB1jGLtzLLfDb1MtbzVyp3ibaYf6YTejwO6KZfY6zAvladxYEsWvvebA16RUTHtQv7t8ubQdK1VKgt62baWPet26UKGCbjRWKomL8SEXhmF4Al8ARU3TDDMMYyHQGJjxyquMr86dk/qwuXOlDc5778Fnnzl6VSox2bkT+veXjTipUsnvW6dOj/Ubu3L3CsP9h/Pj7h8JiQrBt7Avfm/5cfm6J50nnMPl2GtEOB2G85XwyJCdQWOu0r52Fge+KKVikN0uUwuXLIFly6QXL0j7yFKlpEvDqVOajIjHtPe3ii+ie1jsDLgbhhEFJAcuR39J8cjt2zKJbcwYaR/Vo4eWN6iYtXkz9OsHf/4pp2b79ZMDqzRpHt7k7O2zDNk6hOn7pxNlj6Jx8cb4veVH8UzF2bwZWrYOJOhkBQzXKFKVLkqq2uewpDzE7CPuGvyqhM1qlR2auXNLL95q1SAsDN55Rw4O69R5VL+rQW+8pr2/VXzyysGvaZoBhmEMAy4AYcAfpmn+EWMrczTThIoVZbNE8+YSlOgmCfUKnsh2VC+I781jkundtAkyZZIhFR06PNw06X/Rn4VHFnI88Dh/nvkTJ8OJFqVa0K1CN/Kmzcfq1dB+kJQ5Orl7kKbScTxKn8cpmfXh82pbMZUghYdLV4YlS+DXX+Xv4+hR6cG7ahUULw5p0zp6leolae9vFZ9Ep+whLVAXyAMEAYsMw/jYNM05/7pdW6AtQM6cOaOx1Dhgt8vptNq1wcVF+kJlzfrYznqlXsZj2Q7TpNDuv8gz7lMIOCF1vKNHQ+vWMk/4vhn7Z9D619bYTPmg+LDohwyvMZwsybOzaBHUGwwHD0pv/jFj4Ofb/lwNDXniubWtmEpwRo2C3r1ll2aqVPJeXK/eo5ZlFSs6eoXqFWnvbxWfOEXjvlWBs6Zp3jBNMwpYAjyxl9w0zcmmaXqZpumVMWPGaDxdLPvrL+ncUL++zHAH6Qelga+KhqFrThAeGUXNE1tZNaMTPy3+nnT3ghji+yWcPg1ffPEw8N1yYQvv/vwuLZe3fBj4WgwLxTOUYtW87BQqBB99JGd/Z8yQ8sbPP4futQvg7mJ57Hm1rZiK94KCZMBEnTpw5oxclycPNGkCv/0GN27AnDkS/GpJQ4L3rINxPUhXjhCdmt8LwJuGYSRHyh6qAAmvlcOxY9Ctm2w0yp5dooqmTR29KpUY2O2U2v4n07bOo3Dgec6kzUaXWl+yrKgPNosz37i5YZoma8+sZcDmAfx1/i8yJM9AuzLtmHVgFpG2SJxMV0Z19uHWAShbVk5G1KnzcIIxoG3FVAISEgKLFsnlzz/lSC5HDjh/HvLmfdSpQSU62vtbxSfRbXX2HdAIsAL7gNamaUY86/bxrtWZaUK5cnDihDRB7dRJ2kspFR12O/zyi4y4PnKE0+myM6Z8I1YUqYTdSTK02VK70e2DEAZsHsCuy7vI5pGNruW70ub1NoTeScE3Y/yZv2Mj4cd8qFrYGz8/2eOjCTCV4Ny8KVncwoUhMFBaRGbPDg0aQMOGclSnv9hJgnZ7UHHtWa3OohX8vqx4Efzeuyd1lh06yO76o0chY0a5KBUdNptktPr1k9+rIkXY1awjLe7lJuR+ssPEhtV1G85plnHh7gnypMlD97e607xkc65ddmP4cJnAFh4uQwO7d4c33nDsy1LqpQUGwtKlchC4bp3U6m7YID87cQIKFtSAVykV62K8z2+CY7U+Gkd89apkHpo3h6JFHb0yldDZbLBggXRvOHZMfqfmz4cGDXjDYmHAvgB6rV7I2fBFWJ3/JpJAirgVYXbV2TQu3phTJ51p30bKGwE+/lgqcQoXduzLUuqVfPWV7MS02SBfPplM2KDBo58X0tPcSinHShrB76pV8M03ko0rX15a6Hh7O3pVKqGz2STI7ddPslnFikkQ3KDBw6LcsKgwNl8fwVHrSExnEwODgZUH0u2tbuzd40SjhpIgS5YMOnaUuCG+N0VR6qGgIOmQs3ixHL2lTg2vvy5Hbw0bQsmSmuFVSsU7SSP4nTxZNlYsWQK+vvpmrKLHan0U9J48KbODFy2SXen3g967EXf5cfePDPcfzrWQaw/v6mQ4ceo0VP/OiXXrZJZFz57S9EErb1SCEBIiAe+CBbBmDURGyhCK06cl8P34Y0evUCmlnitpBL/TpklGwsXF0StRCdmD8obvvpOgt0QJyXj5+j4MeoPCgxi7YyyjdoziVtgtquatyreFvqXLn12IsEZiWl35qbcPWW0wdCi0awceHg5+XUr9l5AQyfJ6esLFixLgZs8u0wgbNZLCdE0qKKUSiKQR/GbI4OgVqITsQfeGvn2lpve1154Iem+E3GDk9pGM2zmOu5F3qV2wNj0r9uT1zOWYOxcy/V6a82zEM8qHb3t788knUuqgVLwVESH9dufPl1aQdepID/TChWHXLsnyOkWnVbxSSjlG0gh+lXoVDyb+9ekDhw/LRraFC2UQyv0P/YDgAIZtG8akPZMIt4bTsFhDerzVgwKpSjJtGnw4DC5cgBIlvJnn502DBuCsf3UqvuvVC8aNgzt3JHnQrNnj/c+9ntg8rZRSCYZ+DCv1b6Ypma4+fWD/fmnLNHcufPghWCz4X/Rn6fGl/H3zb1afWo3NbuPjEh/T/a3uZHEuzPjx0k3vxg2oUAEmTIBatfSssIqnTBN27pSdl/37y9GZu7uc2WjSBKpU0SM2pVSiou9oSj1gmrB6tQS9e/ZIm6ZZsyQAuP/hv+DwApouafpw/LBvIV9G1BiBe0QeRg6FiRPh7l0Jdv384K23HPmClHqOo0floG7ePBkv7OYm9bulS8suTKWUSqQ0+FVJ2rJ9AQz9/Tj5923lG/95FLt0XHau//STnOq9H/QeunaIAZsHsODIgof3tRgW8ruXZUiPPEyfLg1FGjaUwRSlSjnoBSn1PHa7lOxs3y7tHp2cJLPbu7dMVUmd2tErVEqpWKfBr0qylu0LYMmInxm1YSZvBBzlUqqM9H6vE17ffkndsrkB2H15N/039Wf5ieV4uHrQrEQzFh1dRJQ1CmyuDP/CB5er0KKF9PLPn9+hL0mpJ92+La345syRrgzDh8u/EydKaUOWLI5eoVJKxSkNflXStGMHnh9/yqxTe7maMh29qn/KghLViLK4sH79GTJkCaD/5v78fup30iZLS9+3+/JFuS84ujctZxZ0YGvARtyv+dDxQ286d4Zs2Rz9gpT6l9WrYepUGfITGSldGh4cnVks0L69Y9enlFIOosGvSlr275cR1ytWkCd5avpVbs2cUu8S4eKGiUm40372hi7gremHyJg8I4OrDKa9Vwe2bUhF3RqweTOkT+/N95286dgR0qVz9AtS6j67XVqQlS0ruysXLoRt22R04McfSy2v7rpUSikNflUScfy4bGRbuFDqGvv3p7G1FKfCINzpGPcsvxHhdAqr0wVcSc+oGqNoVaoNq5cnx+dziZmzZ4dRo6B1a0iRwtEvSKn7jh2Tkoaff4bz52HfPik6HzFCMr/aqUEppR6j74oqcTt7ViayzZ4t7Zt69oSvv4a0afl070XaLfuWQKeZYJhgQlpbfca9O4SQQ3l5vRmcOgWFCsn+t6ZNwdXV0S9IqfuOH5dfyr17ZeNatWrSqqxAAfm5npZQSqmn0uBXJU4BARIITJ0q9Y1ffiltGDJmxGa3sfDQPAZsH0Cg5cij+xhO5IkqRtcGebl8GcqUkcFuvr7yEEo5VHg4rFwpmVxfX8iRQ05BjBwJjRvrxjWllHpBGvyqxOXGDRg8GMaPB5sN2rSRbK+nJ1G2KH7eP4OBmwfy962/KZqxKH3f7svgrT8QERWJaXVl75yavFMIZsyAqlW1RFI5mGlKW7KZM2HBAggKkgyvr68Evps2OXqFSimV4GjwqxKH4GCpcRw+HEJDpUdvnz6QJw8R1ghm7J7E4K2DORd0jlJZSvFLw1/wSvkBI0c4Ya6qjpllIxVz+DD0F2/KlXP0i1Hqvlat5EjM3R3q1YPmzaFyZUevSimlEjQNflXCFh4u/UoHDoTAQAkQ+veHIkUIiwpj6o6xDNk2hEvBlyjrWZax744lv/09hgwxaDJHNsh/9JE33bp5U6yYo1+MStLu3oXFi6U+ffZs6Z/30UdQsSI0aACpUjl6hUoplSho8KsSJqtVRg/37QsXL0qNwsCB8MYbrDuzjuE/f832S9u5HX6bijkr8lOdn0h7uyqDexgsWSKTXNu1gy5dIFcuR78YlWTZ7VK6MH26FJiHhkov3vPnJfitVs3RK1RKqURHg1+VsJimZMd695bd7m+8IYFDlSrcCb9Dt5XtmbRnEgBOhhPj3x1PkZBPGdQB/vxTupz5+UGnTpApk4Nfi0q6rFbZuHbliowXTpFCOjc0bw7ly2uxuVJKxSINflXCYJqwdi306AG7d0ORIrBkCfj6cjs8iNEb+zJ6x2iCwoP+cR+DQaPucGkeZM4s++A6dIje2eNl+wIYuuYEl4PCyJbGna41CuFb2jP6r08lfmFhsHSpHKy5usrkNU9PWLNGAt7kyR29QqWUShKcHL0Apf7Tjh2SHateHa5fl+Dh0CECa1Skx/qe5BqVi+/++g6f3D5Mee8nXA13sFuwR7liO+3DhAlw7hx06xb9wNdvySECgsIwgYCgMPyWHGLZvoCYeqUqMdq/X0YJZ80q2d1Tp8DbWw7oQEp2NPBVSqk4o5lfFWdeOmt6/LhkepcuhYwZYfRoaNeOq1G3GbauGxN3TyQsKoyGxRry9Rs92bWyBAO6Q6S1MJnf3EiHmj703OodYwOuhq45QViU7bHrwqJsDF1zQrO/6nHXrkmNTbJkktmdNQvq14eWLcHHR4ZSKKWUcggNflWceJA1fRA8PsiaAk8Gjpcvy0a2n36SFk/ffQedOxNgBjNk/TdM3juZSFskTYo34YvSPVk3vwi1O0pSuHx5GOvnzXvvecd42eTloLCXul4lMVarBLpTp8owipkzpVtDhw6S+U2d2tErVEophQa/Ko68UNY0KAiGDIFRoySQ6NgR/9Y1WHp1I3+v+oTVp1ZjN+00K9GMNoX9+HVGAaq1lBa/NWvKRraKFWNvr1C2NO4EPCXQzZbGPXaeUCUMkZHw/fdSjnP5suyk7NwZ3nxTfq4typRSKl7R4FfFiedmTcPDYcIEGDAAbt2SbFm/fiwK20OTxXWwmRI01y1Uly7FRzFvYm4qN4aICGjYUKYWly4tjxebG9K61ij0WPYawN3FQtcahWLk8VUCEh4OBw9C2bLg4gIrVkCpUjBuHLz/vlynlFIqXtLgVwGx38XgaVlTJ7uNlue2QaFP4cIF2dA2eDAnc6ZgwObvmH1gNiayKcgJC+e3lcPn49w4OUlHqG++gQIFHn8NL1xa8QoePIZ2e0jCDhyAadNgzhyIipJWZSlTws6d0jxaKaVUvKfBr4r1oBH+lTU1TXzO7KH7ppkUvn4WypSBadM4WjIb/Tf1Z8GKBbhZ3Piw2IcsPbacSFsUdqsrx3/3oVMn+Oor6RD1b3GxIc23tKcGu0nRli1SyrB7t7Qpq18f/ve/R10aNPBVSqkEQ4NfFWdBI8DKqcv538of8b5wiHvZc8H8+RysVIj+Wwbyy4RfSO6SnK+9u+AV+TUTh2Ui8pQ/7kU38lF5H37w9yZ9+mc/h25IUzHGNCXQ9fCAwoUluxsZCWPGSLuydOkcvUKllFKvSINfFTdB4+nT+A7uge/ChdK2bNw4Ttb2op//YJZNXoaHqwfdKvhRILAzE/wyMHSPZHdHdvGmTRtvUqT476fQDWkq2u7cgblzYfJk6c/bsqV0HSlVSr7XyWtKKZXgafCrYjdovHkT+vWTDW0uLvh/25LZpSwcvD2PrdM/I7Vbanq91YeMZ75gfMd0nDwJBQtKWeXHH8sZ5helG9JUtHTtKr+noaES7E6cKJsvH9DAVymlEgUNflXsBI3h4XKKeOBAuHsXWrVicuMCdNjqh/2gHYCWr7Uh/7mhTGiemoAA6diwcCHUqwcWy8s/pW5IUy/lzh1YvBhatJChE+7uEuy2bQteXhrsKqVUImWYD0ZsxgEvLy9z9+7dcfZ86sXFWLcHux1+/hl69ZIODu+9x+YuH/L9xdmsPbP24c0MLCTb1o+wP/zw8ZEevdWqPR5vxHYHCpUEmaZ0Zpg0CebPh7Aw2cxWoYKjV6aUUiqGGYaxxzRNr39fr5lfBcRQF4N16+TU8b59mK+X5q9RX/Jd8K9s/Ks5mVJkokWRjsw5+hNWeySm3ZUy6X0Y6v9oFsA/xUUHCpXEnD8PH3wA+/ZBihTQrJlkecuUcfTKlFJKxSEdMK+i7/BhqFULqlbFvHWTtVP8ePuzlLxz8CuOBx6nx+sjqXXyLD9/PA779HWUCurH7Grr2DzX+6mBLzy/A4VSL+zwYVi9Wr729JTpaxMnSn/eSZM08FVKqSRIM7/q1V2+DL17w4wZmKk8+GNQa77PcJhtAYPI5pGNb4qP4dSi1gzq6o6rK7RpA126eJMnj/d/P7S2LVOvKiICfvlFgtytWyFvXnj3XXB2ht9/d/TqlFJKOZgGv+rlhYTAsGH4zxnEBs8o3Lq+xaLcIey4NpUcd3PQueAEDs9uyZDfkpEqFXTrBl9+CZkzv/hTaNsy9UrmzJFhFIGBkD8/DBsm4wB185pSSqn7NPhVL85uh1mzoGdPtjld5p1WTkQ62YFNZAnJQseck9kzrTkjt7qSMaM0eujQAdKkefmn0rZl6oVYrbBqFZQsCblzQ5YsUKkStG8PVapIFwellFLqH/STQb2Y9euhTBnMli1ZViYFDT7NcD/wBQMn7Ns7Mr5VG65ccmXcONlb5Of3aoEvyKa2QfVewzONOwbgmcadQfVe081uSty4AYMGQb584OsrgygAqlaV9mXVqmngq5RS6qk086ue78QJ6NoV+8oVLK2YgX6DcnEg4m883TyxRLpis9swba4kv1aFmTOhSRNwcYmZp46RDhQqcTFNaNcOZs6UccOVK8PIkVCnjqNXppRSKoHQ1Ih6usBA+Pxz7MWLsejyn5Tqm4UGlQMJSe5GI5dZWIedwzZ1I9n/7sfgous4vdGbTz6JucBXqYfCwuDXX+Vrw5Cxf23bwtGj0l6vXj3ZzKaUUkq9AB1yoR4XEQFjx2Ib0I9FOe/Sr04ajjrfpkCawhQN7M2GcY0IDrJQrZqUNfj46F4iFUvOnpWODdOmwa1bcOgQFC/u6FUppZRKIHTIhXo+08R/1kDWLxmG9U4QC9ql5Ji7SYE0Wal+YQKbhjTkVJiFevWge3e4ZAmg15oTXF6j09dUDDt/Hj77TDayOTnJYIqOHaFYMUevTCmlVCKgwa+CPXvY3KclVUsfIrI0YED25JmodHYaW/s34CxONGsG33wDhQvr9DUVC4KD4eJFCXDTpZNa8169pLwhe3ZHr04ppVQiosFvUnblCtaefszdO5OvahpEPvhtMJ24tOJ/3NrzIZ91hK+/hhw5Ht3tedPXNPhVL+Xvv2HcOJg+HXLmlNIGDw8JfrWeRimlVCzQ4DcpCgvDOmIYc37tz4A3Izn1AXi65sEIv4Rp2sDuSkufdxiyBDJkePLuOn1NRZu/P/TvL6OHXVygUSP4/PNHAa8GvkoppWKJBr9JiWkStXAec6Z8wYBiNzldC/JaipBv+yBOr6pD+lLbKdtwI1/V96FqoWePINbpa+qV3L0rQW3KlHD6NOzdC337SuuyLFkcvTqllFJJhLY6SyKidu3gp8YFKbytKa0q3oTU+fHctJwzvY/A8bpMmmQQsN2b1T38nhv4gkxfc3exPHadTl9Tz3TqlMy39vSECRPkukaNZGNbnz4a+CqllIpTmvlNxPwv+rPu4FJC1v7OAvMQZ4tCgbBcZPxjFKe31aVkSYPh86FBA7BY/vvxHnhQ1zt0zQkuB2m3B/UM69fDiBFS2uDsDB9+KBPYQBtCK6WUcpho9fk1DCMNMBUoDphAK9M0/Z91e+3zG3c2nfyTqvNqEmXawYBc99IRuu5HbuxrQMWKBn5+ULOmllaqGGazPTqSqlED9u2DDh2gfXvImtWxa1NKKZWkxFaf39HA76ZpNjAMwxVIHs3HU9EUaY1g5owv6XpmElFuJhiA3Ynz27/mvWwN8RsLFSo4epUq0bl6VQZSTJkC27dL54Zp02THZLJkjl6dUkop9dArB7+GYaQCKgEtAEzTjAQiY2ZZ6mVF2iKZsWoAA/x/4EKyCPIFuxKc1sQ07GC6kKVkOlr/L4AKWpqgYtKBAzByJMybB5GR8P77Mo4YtD+vUkqpeCk6md+8wA1gumEYJYE9QCfTNEP+eSPDMNoCbQFy5swZjadTTxNpi2TG1vEMWPstFyz38Lpp4bUNn/LbmeEke3s1zqV/I4VbIdxS5tA+vCpmBQaClxe4ukKbNvDFF1CwoKNXpZRSSj3XK9f8GobhBWwHKpimucMwjNFAsGmavZ91H635jTmRtkhm7JnGgDU9uWC/TdlLUGbj+6y6Op2gYnfw8DqLc8qIx+5jAGcHv+eYBauE7949mDFDsr1Tpsh1v/4KFStC2rQOXZpSSin1b8+q+Y1Oq7NLwCXTNHfc//4X4PVoPJ56AZG2SCbvmUyBoTlp99unZL1wm8GzS+D6ywE8W6xg/4UMFPc9/0TgC9qHV72iK1egRw+p4/38czhyBEJD5Wd16mjgq5RSKkF55bIH0zSvGoZx0TCMQqZpngCqAEdjbmnqnzad38TI7SPxP7eFa+GBvHkRvt+Yia1B43HvXZ81bQyS399u2LVGIfyWHHpsBLH24VWvZPVq8PUFqxXq1ZNZ197P7wOtlFJKxWfR7fbwOfDz/U4PZ4CW0V+S+qdIWyS9N/Rm6NahmJg42WHwGhfcTvfGMrAL41u4P9EyVfvwqldmmrB2LTg5QZUqUL48fPqpZHzz5XP06pRSSqloi1af35elNb8vLtIWyYz9M+j/1wAu3r0gXZQNcLIbfJK6K9O+/AEnnc+nYkpkJMyfD8OHw8GDUL06rFnj6FUppZRSryw2an5VLIi0RTJlzxTyjy5Iu5XtSHksiBG/QTKrgQUn3FyT0bahrwa+KubMmgV58kDz5jKk4qefZCObUkoplQjpeON4IsoWxYz9M/h+4wAu3TtPrktZWLHRibfPG9z8cjxvtn6NjRe34JPbB+8cWnOpounCBUiXDlKmBLsdihaVoLd6dR37p5RSKlHTsodoWrYvIFq1tVG2KGYemEnf9QMICDlHhoC8jN5wi0angghu2Ia04wdAxoyx+ApUknLwIAwdKkMphgyBr76SOl8NeJVSSiUysTXeOElbti/gsa4KAUFh+C05BPCfAfDm85sZtX0Um85sIzDyKskDijHir8J8efI4kaXKYdk1jrReT/z/UurVbNwowe5vv0GKFDKQokED+ZkGvkoppZIQDX6jYeiaE4+1EwMIi7I9d5JalC2KJvO/YvHf48EwwW7QcG055m3bgZEhE8b06bh98gn/LOqNbnZZJVH/zOh+95305+3fHzp0kJIHpZRSKgnS4DcaLgeFvfD1VruVmftn02Xl9wSZ5x5ebzFNSlp2cq7J/8g3YRikSfPY/aKTXVZJVHg4zJwJY8fCH39AtmwymS1TJnDXQSdKKaWStkTdM2DZvgAqDF5Pnu6rqDB4Pcv2BcTo4z9rYto/r7farUzbMxPPQYVpvaIVQQHpyL2xLW5WA4sNnE2D5cU680mJj54IfOH52WWlHnP7NgwcCLlyQfv2EujeuCE/y5VLA1+llFKKRBz8PsiYBgSFYfIoYxqTAXDXGoVwd7E8dt2DSWpWu5Wpu2aTbUBRWq9swfWLqSi2eRazf8vO2Y2TWbTQgzeuViKtdQjXU1d+qSzy865XSVRwMOTNCz17wuuvw/r1sHMnlCzp6JUppZRS8UqiLXt4lXrcl/W0SWpfVc/P9YiNZO33PYGchKslKXJlMfM8r1Fibydsd4L5sWw9xlRoQqirO273m208L4sc8JRA91m3V0nI6dPw++/QsSOkSgUDBkCFChrwKqWUUs+RaIPfuMqY+pb2JHOGC6w/60/QXSttl3fkltNxuPYape8sZkolT8rs/gyW7wYfH/767FtG7494LDB/kC1+mq41Cj1W8/tft1dJwMGDMHgwLFgArq5Qvz5kySJjiJVSSin1XIk2+I2rjOnWC1t5Z2ZlomyRYABBeSkfuogJH71Nyfm9oMMUCUzmzoXGjaliGAzK++LdG56WXdZuD0nU2bPSomzlShlO8fXX0Lmz/H4ppZRS6oUk2uA3tjOmdtPOyDVL6Lm1I1HO9wNf04nOb7dixK0gaFgEgoIkOOnTR05L3+db2vOlgteXvb1KREwTbt6EDBnkd+jQIfj+eyl10HZlSiml1EtLtMFvbGVMTdNk8PJlDPLvy93kByEkF06pXDGcbLg6OdNw8lxYcxQqVYLx46F48Zh4OSqpsdlg6VIYNEh6Pu/cCenTS52vxfLf91dKKaXUUyXa4BdiNmNqt5t8N/9Xhu/tS4jHfpwiClI/+Rwm9mjMqetr2Ti9Dz4LduAdeRNmz4amTXVylnp5VquMHh44EI4fhwIFoFs3sNsl6NXAVymllIqWRB38xgSr1aTnzFWMPdyXsDR7sNjz8ZH7LCZ81YTUKS0wZw4Zu3TBOzAQPvtCJmk9pV+vUi/k55+hRQsoUUI2tNWvrwGvUkopFYM0+H0K/4v+rD21gYP7krHi7HwiMuzCxTkPLdNNZ1y3j0mezBmOHZPd9Rs3wptvwpo1UKqUo5euEpqwMJg2Tep3P/oIGjeW8oZatR4bca2UUkqpmKHB77+sPbGNmvMrYzMjwADnFFlon20qI5t/QjJXFwgNhR7fwrBhsuN+0iRo3VoDFfVyQkLgxx/l9+jqVYCMx6AAABSzSURBVAl6P/oI3Nzg/fcdvTqllFIq0dLg977AQJOvxm5gzr1WmKkk8DVw4tt3O9L77f/JjVauhM8/h3PnoHlzGDIEMmVy6LpVAjRnDnz5pXRxqFIF5s+XDZJKKaWUinVJPvi9dAk6j9rEktvfYs/5Fxa3DNhMZ8COYbjgYi0BFy5Ap06wbBkULQp//aXBino5t27J2YE0aeRSrhz06gXe3o5emVJKKZWkJNng9+RJ+GrkVlaH9cHMsw539yy8m6U/hy+9RnDkScKdDpEiqiiWQb9i3doYZwOZqtW5s0zVUupFBAbCiBEwdqwMqBgwQMoatLRBKaWUcogkF/zu3QtdRu1gg9kH8q/B3Z6Jr71G0KN6e6oO9yciKgw3ilDhgp3+ayZQOPA8m4uUp+JvcyFXLkcvXyUUgYEwfDiMGyf1vR9+CE2aOHpVSimlVJKXJIJf05RKBb8xe9ju3gcKrsLdTE+XN4fQ7Z1PSeGaAoDLQWGkCQvGb8N0Gh36k0upMtK6Xm/WFSjHWQ181cv47DNYuBAaNYLevaVcRimllFIOl+iD30mr/Onzy1yucQBKbsadtHQtP5AulT7Dw83j0Q1Nk1ZnNvPpyomkigjhx3L1GV2+CWGuyfBM4+64F6AShhs3JNPbqhUULAj9+sG332rQq5RSSsUziTr49b/oz2d73saaOwqAliVaM6rWcFK5pXr8hidOQIcO9N6wgb3Zi+JX/VNOZMwNgLuLha41CsXxylWCceOGtCsbP17a4OXKJcFvgQKOXplSSimlniJRN6fdeG4jJnYALIaFAhnyPh74RkTIRLYSJWDfPpg0iQvLfudegSIYgGcadwbVey3GRiSrRKZPH8idG4YOhbp14cgR6NDB0atSSiml1HMk6syvT24fXC2uRNoicbW44pPb59EPN26Edu2k7UOTJrIjP0sWfAHfMjkctGIV7925A6lTP/r6gw+kZVnhwo5dl1JKKaVeSKIOfr1zeLPuk3VsPLcRn9w+eOfwll34XbrAzJmQN6+MJa5e3dFLVfFdUJAcII0aBb/+Cj4+MHIkGIajV6aUUkqpl5Cog1+QANg7h7e0fJg+XQLf4GDo0UMydu66mU09x927MHq0bGYLCoIGDSBbNvmZBr5KKaVUgpPog18Ajh+H9u2l31mFCjBpEhQr5uhVqfjOboeyZeX3p25dqQ8vWdLRq1JKKaVUNCT+4HfePGjeHFKkgMmT4X//kzGzSj1NWBjMmSMtyywWGDgQcuQALy9Hr0wppZRSMSDxB7/e3tC0qYwmzpzZ0atR8VVEBEydKsHu5csS8NasKRvalFJKKZVoJP4UaO7cUuurga96GqsVpk2T3ryffQb58kknkJo1Hb0ypZRSSsWCxJ/5Veq/DBsGWbNK5rdqVd3IppRSSiViiT/zq9Q/mSasWgXVqkknB2dn2LAB/P3lOg18lVJKqURNg1+VdGzaBBUrwvvvw5kzcPasXJ8liwa9SimlVBKhwa9K/EJC4N134e23JeidOFHal5Uo4eiVKaWUUiqOac2vSrxu3YJ06aTNXapUMGSIbGrTwSZKKaVUkqWZX5X4XLgg/Zxz5JCvARYsgK5dNfBVSimlkjjN/KrE4+ZN6dM7bpx8/+mnkvVVSimllLpPg1+VOAQHS6/eoCCZ6Ne3L+TM6ehVKaWUUiqe0eBXJVxWK6xfD9WrS03vgAHSzaFYMUevTCmllFLxlNb8qoTHNGHpUiheHGrUgAMH5Pr27TXwVUoppdRzafCrEpZNm6B8eahXT3rzLlumLcuUUkop9cK07EElHMHBULs2pEwpo4ibN5cJbUoppZRSL0gzvyp+u3hRNq/Z7VLX+/vv8Pff0spMA1+llFJKvSQNflX8FBwMPXpIB4dBg+DQIbne2xuSJ3fs2pRSSimVYGnwq+IXqxUmTID8+SXorV8fTp6EkiUdvTKllFJKJQJ63ljFLzYbDB0KRYvCsGHg5eXoFSmllFIqEYl25tcwDIthGPsMw1gZEwtSSdDu3dC4MYSFgZsb+PvDhg0a+CqllFIqxsVE2UMn4FgMPI5Kas6fh6ZN4Y03ZFjFsfu/RlmySBszpZRSSqkYFq3g1zCM7MB7wNSYWY5KEiIjoVs3KFQIliyRjW2nTsHrrzt6ZUoppZRK5KJb8zsK+AbwiIG1qKTCxUVKGxo1gv79IUcOR69IKaWUUknEK2d+DcN4H7humuae/7hdW8MwdhuGsfvGjRuv+nQqofvtNyhbFq5ckZKGtWth5kwNfJVSSikVp6JT9lABqGMYxjlgPlDZMIw5/76RaZqTTdP0Mk3TK2PGjNF4OpUgHT0KNWtCrVoQFCTBL4Crq2PXpZRSSqkk6ZWDX9M0/UzTzG6aZm6gMbDeNM2PY2xlKmGz2+Hzz6FECdixA0aOhMOHta5XKaWUUg6lfX5VzLLbwclJLvfuQfv2Mp44QwZHr0wppZRSKmYmvJmmudE0zfdj4rFUAmWasHIlFC8O+/fLdT/9xP/bu/tYK8o7gePfn1TYrqW8iStaWVeCWzEqAqEu7lYjFl+Copt11RpWfIm1i0Zb1qyJjW23TRqlmCBaLaj1NUrctquhNtaqrW5SXQsFr5YK2ErEgu76hpY2LPDsH8+wvT2cc7mHezlz7pnvJ5ncOTNzcn73l+fM+eWZZ57hllssfCVJUtvw8cbqu64umDEDzjgj9/xu2ZK3O1evJElqMxa/6psvfhEmToTly+Hmm3MhPG1a2VFJkiTV5ZhfNW/bNhg0KPfsDh8OV1wBX/4yjBxZdmSSJEk9sudXzXnqqdzTu2xZfn399bBwoYWvJEkaECx+1Tvr18M558D06fC738GQIWVHJEmS1DSLX+3ewoXwyU/CD34AX/tafnDFjBllRyVJktQ0x/yqvpTyss8+MGwYzJoFN94IY8eWHZkkSdIes+dXu3rpJTj55DxHL8CcOfDQQxa+kiRpwLP41R+9+y5cdVW+oe0Xv4ChQ8uOSJIkqV857EHZo4/CpZfC22/D5z6Xx/aOGlV2VJIkSf3K4rfqUsrz9Q4bBocfDosWwbHHlh2VJEnSXmHxW1XvvZfn6B0yBObPhxNOgGef9ZHEkiSpoznmt2pSgvvuy1OX3XILbN2at4GFryRJ6nj2/FbJmjV5XO+zz8KnPpXn7Z08ueyoJEmSWsbit0oi4NVXYckSuPjiPIevJElShVj8drKU4MEH4ac/hW9/G8aPh9/8BgYPLjsySZKkUtj116lWr4aTToILLoAVK2Dz5rzdwleSJFWYxW+n+f3v4UtfgmOOgVWr4Pbb4bnn4OMfLzsySZKk0jnsodNs2ZKHOJx7LixYAAccUHZEkiRJbcOe306waRNcdx1s356fyvbLX+bpzCx8JUmS/oTF70C2fTt861t5zt4FC2D58rx99Ohy45IkSWpTFr8D1cqVMG0azJ0LU6ZAVxdMnVp2VJIkSW3NMb8DUUoweza89Rbcfz989rM+nU2SJKkXLH4HkmXL4IQTYOhQWLoUxoyBESPKjkqSJGnAcNjDQPDGG3D22XDGGXDrrXnbhAkWvpIkSU2y+G1nO3bkeXonTIDHH4cbboB588qOSpIkacBy2EM7u+YauOkmmD49z907blzZEUmSJA1oFr/tZutW+PBDGDkSLr8cjjoKLrzQG9okSZL6gcMe2snPfgaTJsEll+TX48fDnDkWvpIkSf3E4rcdfPABXHklHH88bN4Ml15adkSSJEkdyWEPZVu+HM46K8/ocOWV8PWv56nMJEmS1O8sfsuSUh7OcOihcPjh8PDDcNxxZUclSZLU0Rz20Gopwb33wimnwLZtMGoUPPmkha8kSVILWPy20oYNMHNmnr1hyxZ4992yI5IkSaoUi99WSAmWLIEjj4Sf/AQWLoRnnoHRo8uOTJIkqVIc89sKW7fmh1VMngx33AGHHVZ2RJIkSZVkz+/esmMH3HlnfmDFkCHw9NPw4x9b+EqSJJXI4ndvWLsWTjwxz9d7991524EHwj6mW5IkqUxWY/1p+3ZYsACOPhq6uuA734G5c8uOSpIkSQXH/PanefPyzWxnngm33QYHHVR2RJIkSerG4revtm/P43qHDctPaJs6Fc4/Pz/AQpIkSW3F4rcv1qzJc/aOHg2PPALjxuVFkiRJbckxv3tix448vGHiRHjlldzTK0mSpLZnz2+zNmyA2bPzwypmzoTFi2HMmLKjkiRJUi9Y/Dbrox+FjRvhrrtgzhzH9kqSJA0gDnvojddfhy98AbZtg1Gj4OWX4aKLLHwlSZIGGIvfnqQE99wDRx0FS5bAqlV5+6BB5cYlSZKkPWLx28imTTBrVh7acPTRufCdPLnsqCRJktQHjvlt5Jxz4IUX8hPbrrrK3l5JkqQOsMfFb0QcAtwLHAjsABanlBb2V2CleO892Hdf2G8/WLQIhgyBI44oOypJkiT1k74Me9gGzEspHQEcB8yNiAn9E1YJnnoqj+295pr8euJEC19JkqQOs8fFb0ppY0ppRbH+AbAaOLi/AmuZP/wB5s2D6dPzNGZz5pQdkSRJkvaSfhnzGxGHAscCz9fZdxlwGcDYsWP74+P6z+rVcO650NUFn/88zJ+fhzxIkiSpI/W5+I2IjwHfBa5OKW2u3Z9SWgwsBpgyZUrq6+f1q8GDc8/vY4/BaaeVHY0kSZL2sj5NdRYR+5IL3wdSSt/rn5D2svXr4atfzXP4jhuXe38tfCVJkiphj4vfiAjgTmB1Summ/gtpL0kJ7r8/z9n7zW/CunV5u1OYSZIkVUZfen6PB2YDJ0XEymI5vZ/i6l/vvJPH9s6enWd0ePFFGD++7KgkSZLUYns85jel9J9A9GMse0dK8JnP5JvavvGNPJWZvb2SJEmV1PlPeIvIszgMHw6TJpUdjSRJkkrU+cUvwEknlR2BJEmS2kCfZnuQJEmSBhKLX0mSJFWGxa8kSZIqw+JXkiRJlWHxK0mSpMqw+JUkSVJlWPxKkiSpMix+JUmSVBkWv5IkSaoMi19JkiRVhsWvJEmSKsPiV5IkSZVh8StJkqTKsPiVJElSZVj8SpIkqTIsfiVJklQZFr+SJEmqDItfSZIkVUaklFr3YRH/Daxv2Qf+0f7A/5TwuQOV+WqeOWuO+WqO+WqO+WqO+WqO+WpOmfn6y5TS6NqNLS1+yxIRP08pTSk7joHCfDXPnDXHfDXHfDXHfDXHfDXHfDWnHfPlsAdJkiRVhsWvJEmSKqMqxe/isgMYYMxX88xZc8xXc8xXc8xXc8xXc8xXc9ouX5UY8ytJkiRBdXp+JUmSpM4qfiPi1Ih4JSLWRcS1dfYPiYilxf7nI+LQ1kfZHiLikIh4OiJWR8TLEXFVnWNOjIj3I2JlsVxfRqztIiJei4iuIhc/r7M/IuLmon29GBGTyoizHUTEX3drNysjYnNEXF1zTOXbV0TcFRFvRcRL3baNjIgnImJt8XdEg/deWByzNiIubF3U5WmQr/kR8aviO/f9iBje4L09fn87UYN8fSUi3uj2vTu9wXt7/D3tRA3ytbRbrl6LiJUN3lvF9lW3jhgQ57CUUkcswCDgVeAwYDCwCphQc8w/A7cX6+cBS8uOu8R8jQEmFetDgTV18nUisKzsWNtlAV4D9u9h/+nAD4EAjgOeLzvmdliK7+Ym8nyL3bdXvn0BnwYmAS9123YjcG2xfi1wQ533jQR+XfwdUayPKPv/KSlfM4CPFOs31MtXsa/H728nLg3y9RXgX3bzvt3+nnbiUi9fNfsXANc32FfF9lW3jhgI57BO6vmdCqxLKf06pbQVeAiYVXPMLOCeYv3fgekRES2MsW2klDamlFYU6x8Aq4GDy41qwJsF3Juy54DhETGm7KDawHTg1ZRSGQ+4aWsppWeAd2o2dz9P3QOcVeetpwBPpJTeSSm9CzwBnLrXAm0T9fKVUvpRSmlb8fI54BMtD6xNNWhfvdGb39OO01O+ilrhH4EHWxpUG+uhjmj7c1gnFb8HA693e72BXYu5/z+mOFm+D4xqSXRtrBj+cSzwfJ3dfxMRqyLihxFxZEsDaz8J+FFELI+Iy+rs700brKLzaPyDYfva1V+klDZC/nEBDqhzjG2tvovJV1/q2d33t0quKIaJ3NXgkrTta1d/B7yZUlrbYH+l21dNHdH257BOKn7r9eDWTmXRm2MqJSI+BnwXuDqltLlm9wrypepjgEXAf7Q6vjZzfEppEnAaMDciPl2z3/ZVIyIGA2cCD9fZbfvac7a1GhFxHbANeKDBIbv7/lbFbcA4YCKwkXwpv5bta1fn03Ovb2Xb127qiIZvq7OtZW2sk4rfDcAh3V5/Avhto2Mi4iPAMPbsklBHiIh9yQ32gZTS92r3p5Q2p5Q+LNYfA/aNiP1bHGbbSCn9tvj7FvB98qXB7nrTBqvmNGBFSunN2h22r4be3Dlcpvj7Vp1jbGvdFDfLzAQuSMWAwlq9+P5WQkrpzZTS9pTSDmAJ9fNg++qmqBf+Hlja6Jiqtq8GdUTbn8M6qfh9ARgfEX9V9DadBzxac8yjwM47Cv8BeKrRibLTFeOX7gRWp5RuanDMgTvHREfEVHJ7ebt1UbaPiNgvIobuXCffZPNSzWGPAv8U2XHA+zsv/VRYw94S21dD3c9TFwKP1DnmcWBGRIwoLlvPKLZVTkScCvwrcGZKaUuDY3rz/a2EmvsQzqZ+Hnrze1olJwO/SiltqLezqu2rhzqi/c9hrb47cG8u5Lvt15DvUr2u2PZv5JMiwJ+RL7+uA/4LOKzsmEvM1d+SLzG8CKwsltOBy4HLi2OuAF4m3+n7HDCt7LhLzNdhRR5WFTnZ2b665yuAW4v21wVMKTvuknP25+Ridli3bbavP83Rg+RLz/9L7gm5hHwfwpPA2uLvyOLYKcAd3d57cXEuWwdcVPb/UmK+1pHHDu48j+2c0ecg4LFive73t9OXBvm6rzg/vUguUsbU5qt4vcvvaacv9fJVbL9753mr27G2r8Z1RNufw3zCmyRJkiqjk4Y9SJIkST2y+JUkSVJlWPxKkiSpMix+JUmSVBkWv5IkSaoMi19JkiRVhsWvJEmSKsPiV5IkSZXxf0YpYtWYxtpuAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure(figsize=(12,8))\n", "ax = fig.add_subplot(111)\n", "ax.plot(x1, y2, 'o',label=\"data\")\n", "ax.plot(x1, y_true2, 'b-', label=\"True\")\n", "prstd, iv_l, iv_u = wls_prediction_std(res)\n", "ax.plot(x1, res.fittedvalues, 'r-', label=\"OLS\")\n", "ax.plot(x1, iv_u, 'r--')\n", "ax.plot(x1, iv_l, 'r--')\n", "ax.plot(x1, resrlm.fittedvalues, 'g.-', label=\"RLM\")\n", "ax.legend(loc=\"best\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 2: linear function with linear truth\n", "\n", "Fit a new OLS model using only the linear term and the constant:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[5.66410251 0.39105655]\n", "[0.39504002 0.03403825]\n" ] } ], "source": [ "X2 = X[:,[0,1]]\n", "res2 = sm.OLS(y2, X2).fit()\n", "print(res2.params)\n", "print(res2.bse)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Estimate RLM:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[5.13515772 0.48757594]\n", "[0.10878613 0.00937345]\n" ] } ], "source": [ "resrlm2 = sm.RLM(y2, X2).fit()\n", "print(resrlm2.params)\n", "print(resrlm2.bse)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Draw a plot to compare OLS estimates to the robust estimates:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFlCAYAAAAzqTv+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hUVRPA4d/dTaUGQpPQQpVOEFRAMBQFFKQIShULRUBRP0VAEFSkKCqKhCpVEVDEgCiiNEEIvQXpSgmhlwAhffd+fwwBwYSEtN1N5n0eHsnuJnui1517zpkzY5imiVJKKaWylsXRA1BKKaVyIg3ASimllANoAFZKKaUcQAOwUkop5QAagJVSSikH0ACslFJKOYBbVr5ZoUKFzDJlymTlWyqllFIOs3379gumaRZO6rksDcBlypRh27ZtWfmWSimllMMYhnE8ued0CVoppZRyAA3ASimllANoAFZKKaUcIEv3gJMSHx/PyZMniYmJcfRQMpWXlxclSpTA3d3d0UNRSinlBBwegE+ePEnevHkpU6YMhmE4ejiZwjRNLl68yMmTJ/H393f0cJRSSjkBhy9Bx8TE4Ovrm22DL4BhGPj6+mb7Wb5SSqnUc3gABrJ18E2UE35HpZRSqecUAdiZvPfee3zyySfJPh8cHMy+ffuycERKKaWyI5cLwME7w2kwdjX+g3+mwdjVBO8Mz9r31wCslFIqA7hUAA7eGc6QxaGER0RjAuER0QxZHJruIDxq1CgqVapEs2bNOHjwIADTp0+nbt261KxZk6effpqoqCg2btzI0qVLGThwILVq1eLvv/9O8nVKKaVUSlwqAI9bcZDoeNttj0XH2xi34mCaf+b27dtZsGABO3fuZPHixWzduhWA9u3bs3XrVnbv3k3lypWZMWMG9evX56mnnmLcuHHs2rWLcuXKJfk6pZRSKiUOP4Z0L05FRN/T46mxfv162rVrR65cuQB46qmnANi7dy/Dhg0jIiKCyMhImjdvnuT3p/Z1SimlnNylS1CwYJa9nUvNgIv7eN/T46mVVIby888/z8SJEwkNDWXEiBHJHiFK7euUUko5qbAwaN8e/Pzg9Okse1uXCsADm1fC291622Pe7lYGNq+U5p/ZqFEjfvzxR6Kjo7l27Ro//fQTANeuXeO+++4jPj6eefPm3Xx93rx5uXbt2s2vk3udUkopJxYTAzdyfvDxgd274Y03IAurFbrUEnTbAD9A9oJPRURT3Mebgc0r3Xw8LWrXrs2zzz5LrVq1KF26NA0bNgRg5MiRPPTQQ5QuXZrq1avfDLqdOnWiV69eTJgwgUWLFiX7OqWUUk4oPBwmTYJp06BIEdi7F/LmhcOHwZK1c1LDNM0se7M6deqYd/YD3r9/P5UrV86yMThSTvpdlVLKqezeDWPHwqJFYLNBmzYwYAAEBkImFkoyDGO7aZp1knrOpWbASimlVKrFxUFCAuTKBX/9BcuXS9B95RVwgrr8LrUHrJRSSqXo7Fl4/30oVQomTJDHOnaEkyfh00+dIviCzoCVUkplF9u3wxdfwMKFMvtt2RLq1ZPn3N2zNMEqNVKcARuGMdMwjHOGYey94/FXDcM4aBjGX4ZhfJx5Q1RKKaWSYbff+vvQofDjj9C7t2Q4//ILPPqo48aWgtTMgGcDE4G5iQ8YhtEYaAPUME0z1jCMIpkzPKWUUioJFy7A9OkwdSr88QeULi3Zzb6+kD+/o0eXKikGYNM01xmGUeaOh/sCY03TjL3xmnMZPzSllFLqDrt3y77uvHkQGwuPPQaRkfJc2bKOHds9SusecEWgoWEYo4AY4C3TNLcm9ULDMHoDvQFKlSqVxrfLPBcvXqRp06YAnDlzBqvVSuHChQHYsmULHh4ejhyeUkqpRBcuQJ064OEBL7wAr74KVao4elRpltYA7AYUAB4G6gLfGYZR1kziULFpmtOAaSDngNM60Mzi6+vLrl27AOkFnCdPHt56663bXmOaJqZpYsniQ9pKKZWjXboEX30lR4jmzIFChWDxYnjkEShQwNGjS7e0RpSTwGJTbAHsQKGMG5bjHTlyhGrVqvHyyy9Tu3ZtwsLC8PHxufn8ggUL6NmzJwBnz56lffv21KlThwcffJBNmzY5athKKeX69u6FPn2gRAkYNAhOnIDoG013WrfOFsEX0j4DDgaaAGsNw6gIeAAX0juY11+HG5PRDFOrFnz+edq+d9++fcyaNYspU6aQkJCQ7OsGDBjA22+/zcMPP8yxY8do1aoVe/fuTfb1SimlkrFwIXTqBF5e0K2bLDPXqOHoUWWKFAOwYRjzgUCgkGEYJ4ERwExg5o2jSXFAj6SWn11duXLlqFu3boqvW7lyJQcP3upJfPnyZaKjo/H2Tl+XJqWUyvYiImDmTCmO0a4dNG8OY8ZAz56y5JyNpSYLunMyT3XL4LGkeaaaWXLnzn3z7xaLhX/fY/y77aBpmpqwpZRS92L/fvjyS9nbjYqSs7vt2klnosGDHT26LKFZRalksVgoUKAAhw8fxm638+OPP958rlmzZgQFBd38eldGr6MrpVR28tprkr08cyY8+yzs2CHneXMYDcD34KOPPqJFixY0bdqUEiVK3Hw8KCiIDRs2UKNGDapUqcL06dMdOEqllHIyV65IiciICPm6USMYNQrCwiQIBwQ4dnwOou0Is1BO+l2VUooDB2DiRJg9G65fh2++ga5dHT2qLKXtCJVSSmWdmBho2xZWrJCiGZ07SzbzAw84emRORZeglVJKpd+VK9L8AOQIUaFCMHKkLDPPnq3BNwk6A1ZKKZV2Bw/eWmaOiYFTp6BwYVluVnelM2CllFL3bt8+aNEC7r8fpk2TI0QhIRJ8VaroDFgppVTqXLki9Zn9/SFPHgnCH3wgZ3iLFnX06FyOBmCllFJ3d/DgraIZjzwCy5dDqVJw7Bhok5o0039zwMmTJ2nTpg0VKlSgXLlyvPbaa8TFxbF27VpatWr1n9cvW7aMgIAAatasSZUqVZiaAw+QK6VygDVrpDTk/ffD9OnQvr0kViXS4JsuOf7fnmmatG/fnrZt23L48GEOHTpEZGQkQ4cOTfL18fHx9O7dm59++ondu3ezc+dOAgMDs3bQSimVWa5cgbg4+fvmzRAaeiubec4c6cerMkSOD8CrV6/Gy8uLF154AQCr1cr48eOZOXMmUVFR/3n9tWvXSEhIwNfXFwBPT08qVaqUpWNWSqkMt38/9O8Pfn6waJE8NmAAHD8Ow4ZBkSKOHV825Fx7wA7oR/jXX3/xwB3n0/Lly0epUqU4cuTIf15fsGBBnnrqKUqXLk3Tpk1p1aoVnTt3xqJLMUopV2O3w88/y/7u77/fKpqR2P4vVy7Hji+by/FRwzRNDMNI9eMAX331FatWreLBBx/kk08+4cUXX8zsYSqlVMaJj5d/GoZ0Htq3Dz788FbRjGrVHDq8nMK5ZsAO6EdYtWpVfvjhh9seu3r1KmFhYZQrVy7Z76tevTrVq1ene/fu+Pv7M3v27EweqVJKpdO+fTLbXbIEDh2So0TLlkGJEuDu7ujROVx0NJw8CRUqZM375fgZcNOmTYmKimLu3LkA2Gw23nzzTZ5//nlyJbH8EhkZydq1a29+vWvXLkqXLp1Vw1VKqXtjs0nAbdYMqlaFWbOkgEZkpDzv75/jg29CAnz1lQTedu1kZT4r5PgAbBgGP/74I99//z0VKlSgYsWKeHl5MXr0aABWrVpFiRIlbv7ZuXMnH3/8MZUqVaJWrVqMGDFCZ79KKecVGiqNEQ4ehNGjZYo3cyYUK+bokTmcacIPP8iKe69e4FM9hHpvj2FzeEiWvL+2I8xCOel3VUo5yJ49sszs7g6TJsljK1dCYCC4OdeuY2oF7wxn3IqDnIqIpriPNwObV6JtgF+6fuaqVbL9vW0bVKkCzw5exYfHWpJgT8DLzYtVz62iXsl66R773doR5vgZsFJKubyEBJnKBQZCzZowbx5YrTLFA1l+duHgO2RxKOER0ZhAeEQ0QxaHErwzPE0/b9s2eOwx+Vdy7hxMnRHNi9PH83FYG+Lt8ZiYxNniWHtsbYb+HknRAKyUUq7ugw+gQwcpDfnxx7LM/OWXkuXs4satOEh0vO22x6LjbYxbcfCefs7Bg9CxI9StK6ddx30Wx1sLpvBBRAXe+v1/FMtVDgN3MC2YphtGXNWM/DWS5Jq3REoplZPt2CEBtkcPmfW++KL0223VSma+2cipiOh7evxOJ0/K/cnMmeBRNoTGw1dT5+EYJh2ax9HfjlKvRD36VP+Mb9flo6htLzGWULzs1Zmz1oP7C4ane6n7bjQAK6WUK4iPl2XmL7+EjRulSEa9ehKAy5SRP9lQcR9vwpMItsV9vO/6fZcuwdix8q/LZoP2r20g2KcJa+xxrNkCFQtW5OcuP9OyfEse+WgN0fHReFIZT7vk6UTbZZadmQFYl6CVUsrZmaYE286d4exZGD8ewsOlDWA2N7B5Jbzdb5/Ve7tbGdg86RLA169LsnfZsvDJJ9DxGZMpq3/mz+IdibdLjWsLFnrU6sETFZ7AMIx0z7LTSgOwUko5G9OETZvglVckwcow4M03pWzkoUNSttfHx9GjzBJtA/wY0746fj7eGICfjzdj2lf/z8w0Ph4mT4by5WHoUHj0UZi+cjVHGjXgpVWtMAwDD4sHVsOKp5snjcs0vvm9yc2mU5plp5cGYKQBQ61atahWrRqtW7cmIiICgGPHjlEtiZJsiUU6rl27dvOx1157DcMwuHDhQpaNWymVzcTGwty58OCDMuOdOxf++kue69wZnngiR7YAbBvgx4bBTTg69kk2DG5yW/C12+HNMZfIVyyKfv0gyiuCvlOWE/l0U3qub8qJKyeY2moqx147xtrn1zKy8cj/HDG611l2Rsl5/yWT4O3tza5du9i7dy8FCxYkKCgoxe8pX748S5YsAcBut7NmzRr8/DJvr0Aplc0dOQIlS0pi1fXrEBQky8w1azp6ZE7JNGH5cihXOY7P3ilIfPGNeL36DDHd6zP5zBNsC9/N+ObjOTLgCL0f6I271Z16JesxpOGQ/5zvTe0sO6O5ZBJWSFgIa4+tJbBMYIYclP63evXqsWfPnhRf17lzZxYuXEi3bt1Yu3YtDRo0YPny5Rk6FqVUNmaasH69pOl26SKblu3by3Gipk2zxRGizBISIkU01q0Dz4I28nQfQ2TZYdgMO5iQJ6EllT368/rDT6b6Z7YN8Mv0gHsnpwrAr//6OrvO3L0d4ZXYK+w5uwe7acdiWKhRtAb5PfMn+/paxWrxeYvUNXmw2WysWrWKl156KcXXVqhQgSVLlnD58mXmz59Pt27dNAArpVIWFSWFMiZOlKpV5ctDp06ytDxliqNH59T++kv2d5csgaJF4f0JR/n4bC8i3VYDiVUdLbiZhTl7xZEjTR2XW4K+EnMFuymVsu2mnSsx6f+3HB0dTa1atfD19eXSpUs89thjqfq+9u3bs2DBAjZv3kzDhg3TPQ6lVDY3f740vE/MXp4+HXbvzpH7uvfi+HF4/nmoXh3WrIG3R4bTakpfRkZUJNp9PblsDTHwANOCgRte9uqZnkCVEZxqBpyamWpIWAhN5zYlzhaHh9WDee3npXsZOnEP+MqVK7Rq1YqgoCAGDBiQ4vd16tSJ2rVr06NHDyz6P5BS6k6mKXWYy5SRVjv+/lIH8dVX4ZFHdJk5BefOyZGiyZPBLBHCg4N/pni1v5lwNBhbqI1etXvxQIEX+WT5BSLibhXR8LFWy/QEqozgVAE4NeqVrMeq51Zlyh5w/vz5mTBhAm3atKFv374pvr5UqVKMGjWKZs2aZdgYlFLZwNWrksE8caLUQHz9dTm7+/DD8N13jh6d07t6FT77DD79VFbsH+/7GyuLPMlmMwEOwxMVnmBiy4n4F/AHwNc7nHErPDgVUTnDmjVkBZcLwCBBOKOTrxIFBARQs2ZNFixYQMOGDTl48CAlSpS4+fz48eNve32fPn0yZRxKKRc1ZIgE3shIeOgh+PprKUKsUhQbK9vgH34IFy5Am46R+HeawJT9I0lISADAalh5pOQjN4MvOCaBKiO4ZADOaJGJjalv+Omnn27+PT4+/j+v75jM/0zHjh3L0HEppVyAzQarV0t7HcOQKNK+vRTRqFvX0aNzCim1E7TZ4JtvYPhwOHECApvFUKvXZOYdH8OS0PM0KNmAbae2kWBPwMPqQWCZQMf9MhlIA7BSSqXFxYtS4X/SJOlCtGED1K8va6fqpsR2gokdjRLbCQK0qeXH0qXwzjuw72oIxRus4pl3r7Hh6jzW7g+nqX9TPmzyIQ+XeDhTj586SooB2DCMmUAr4JxpmtXueO4tYBxQ2DRNLQGllMr+zp+XQ6jffgsxMdCoEYwbJ9Wr1H8k105w2OSzfBTqx6ZNUKL+n7g925RTZhzfhUO1wtX4ut3XNPa/VS4yM7ceHSU1M+DZwERg7r8fNAyjJPAYcCLjh6WUUk4kLk5muRUrQt68sGqVVKzq31/Oxqhk3dnQIO5sPi7/UYnjR4tQ3M/OyxMW80NMPxKibjRKMCx0rt75tuCbXaUYgE3TXGcYRpkknhoPvA0sSe8gTNPEyObp+KZppvwipZRzOX0apk2TzKBcuaQRgpeXlI100x281EhsJxh/ORcR6ysStd8PwyuWIp2+onCjSUw5t5PS+UvjYfXAZrfhYfW4rVFCdpamK8gwjKeAcNM0d6cUOA3D6A30Bjm2cycvLy8uXryIr69vtg3Cpmly8eJFvLy8HD0UpVRq7NkDY8bAokXSjeiJJySpKvEzSoNvqr1UuzJvvhNPxM4SGFY7uZ78BnvtMZyz7iNPfFnmtp1Ll+pd2BK+Jdvt8abknq8iwzByAUOBx1PzetM0pwHTAOrUqfOfaWCJEiU4efIk58+fv9ehuBQvL6/bjjMppZxMdLT0tMuXD/75Ryr9v/oq9Osn5SLVPYmIgI8/hs8/v4+4eJO87T4hqvJ4otxO4+tVjC+aTeWFWi/gbnUHsuceb0rSchtXDvAHEme/JYAdhmE8aJrmmXv9Ye7u7vj7+6f8QqWUygxHj0qppRkzZE/3gw+gdWvpRJQ7t6NH53Kio+UY9JgxcPkytHx+N2fr9GfHhQ0AuFvcWfTs/GxzlCg97jkAm6YZChRJ/NowjGNAHc2CVkq5lJUrYcIEWLZMajG3bw8tW8pzVqsG33uUkACzZsF778GpU9Cw3QG8W45g+anv8IrwwsDAxMRu2gkJC9EATCqaMRiGMR8IASoZhnHSMIyUWwUppZQzioq69ffJk2HzZmmvc+yYlIisl7OWQDOCacL330PVqtJjouj9R2kx5QU21KrKxgu/MKzhMJZ0WoKXmxdWw5qtCmmkl5GV2bl16tQxt23blmXvp5RSgPSxCwqSspDbtkGlSnDmDBQoAJ6ejh6dy1q5Uo5Ebz8bQqEGSyn/8EG2X1uGxbDQv25/Bj8ymMK5CwOZ28fdmRmGsd00zTpJPaepfEqp7CkhAZYulQ3JNWsk0HbuDO6S9EOxYo4dnwvbulVKXq9aBYXr/4KlZxsukMCFK9Du/nZ82fJL/PLdXps5JyZZpUR76Cmlshe79AvnyhXo2lUymj/6CE6elE3KsmUdOz4XduAAdOggRb927o/g8THvcqVlO+zcapRQt3jd/wRflTSdASulXJ9pyn7uxIkQFgZ//AG+vhASIpWqrFZHj9ClhYXB++/L/Yt3/kiavDeBHZ7j+C02gib+TdgYtpF4W7zu794jDcBKKdcVHQ0LFsj+7vbtUiby+eeldKSHB9Sq5egRurSLF+GlN66xZN1uKLOKXM8dx1rhJ1bHX6BV6VaMbDySWsVq5dj93fTSAKyUcl3z5kGvXlClinQl6tZNgrBKl8hI+PxzGPORnaiCO+H5ZmDEE2WAZ2wFxjaawaCmT918fbbY3zVNyRU4cwa6dMmSt9Q9YKWUa7Db4bff4Kmn5AgRSFLV6tWwdy/07avBN53i4mQxoXx5eHe4DaPhZIwercASDwZgGnjb6rF0ax5HDzXjXLkCX34pN3FNm8LIkRKMs4AGYKWUc4uIkOnY/fdD8+ay15sod25o3PhWjWaVJna7LCbcfz+88qqdAo8soszH1bn+0CtYjXxguoFpwcAdL3v1/3Q4ckl79sDLL4OfHwwYICVI58yBHTuy7HrSJWillHN75hn4/Xdpdv/ee/D003p2N4OYppS8HjAuhL8T1nBfgCfl+s7jQNROKntXpmLscGKu1SHOcpAYSyhe9up42itT3Mfb0UNPm7g4+OEH2a7480/pbNW5s9T7rpPkUd1MpQFYKeU84uPhxx9h6lSYPx+KFJElwbFjoXZtR48uW9mwQc7yrj8aAi80Bmssp4Hi1uI3OxT9tPsMQxaHYsRXxtNeGQBvdysDm1dy7ODvVViYXFNffQVnz0K5cvDpp5KwV7Cgw4alAVgp5XinTknf3WnTpAevv7+c3y1SBB56yNGjy1ZCQ6X65k8/QcEamyjUuzsX7LEAWLDQt05futfsDkDbADnPO27FQU5FRFPcx5uBzSvdfNypmaZUCgkKkoIspgmtWsls9/HHpf63g2kAVko51oULEnDj4qQZwvTp0KKFnt0FgneGpzv4Jf6M48cgbktlLuwqRi7/PVR6fxgHzWX4uPvgHu+O3bTjYfWgqX/T276/bYCfawTcRBERspc7aRIcOgSFCsHbb0OfPlCmjKNHdxsNwEqprBUZKRk/hw7JMmChQlJAo3Fj7bv7L8E7wxmyOJToeBsA4RHRDFkcCpDqgBi8M5yBcw9ydr0/13aWhkIH8Ojdk+v3/cJZTx9G1R/FgIcGEHo21PXP8e7aJUF33jxpulGvntT+7tjRaXMGtBmDUiprHDokH5CzZsHVq7Knu3Gj0344OlqDsasJTyLb2M/Hmw2Dm6T4/VevQtV2RwlfXxKzfDCWx9/H7vMXBp74uXUg9K0v8fHyyYyhZ53YWFi0SK6rjRvB21vKj/brBwEBjh4doM0YlFKO9vXX8Nxz0gihY0dpfF+vnh4fuovkjvqkdAQoJkaOSY8aBRfjPLF2fwpbyeVIhWwLheIG4xZTx7WD7/Hjt5Kqzp+HChVg/Hjo0UM6XLkIDcBKqYx3/jzMmCGlIFu0gMcek2zmXr2gaFFHj84lFPfxTnIGnNwRoIQEuc8ZMQLCLl6gdNexXCo+ERtx8gIDMCHe8g/F8zXMxJFnErtdjqNNmgTLlsljTz0ls92mTZ0iqepeud6IlVLOKbEhwnPPQYkScsZl5Up5rlgxGDZMg+89GNi8Et7utyeiJXUEyDQhOBhq1IAX+0VgazQc7yH+hJUYT2CppyhuewcDjxuFNNzIb9RyrWNEly/DZ59JD+cWLaTBxuDBcPSoHFl77DGXDL6gM2ClVEbp0kUaI+TNKzPdfv2kvJ9Kk9QcAVq7FvqPDWFfzG/kr32GPF0Wcsp2mY6VOvJ+4PtULlyZ4J3hDF9elLCobZTMVYcPWrZ3jazm7dtltjt/vjTdaNBAWjJlo0IsmoSllEqbf/6RPbhhwyBXLgm+ly5B9+5ak5mMOUKUnJ07ZYFhxYG18NzjYI0HoH6J+kx8YiIB9zlHAtI9i4mB776TwLt5s1xX3brJzVzNmo4eXZpoEpZSKmPY7fDrr1LcYPlyWfpr0gSaNYNOnRw9OqeREUeIknLkCLz7Liz4Lp5cDWbh+cJAYpHgazEstKrYyjWD79GjklQ1Y4acC69UCb74QrYzfFw4WSwFrrlwrpTKehcuyDndJ5+UgvXvvivZqM2aOXpkTmfcioM3g2+i6Hgb41YcTNPPO3VKmj3dX8XG4iPfUODdykQ17UP5wiXxtHpiNax4Wj0JLBOYAaPPIna73MS1aiWlIceNg4YNJW9g/35pkJCNgy/oDFgpdTfbtkmrv+efl4IZjz8uBTPatZOG9ypJaT1CdKeICGmUMD9kDTYb+Ayax2W3fVQpVouvGy/jiQpPsOnkJtcqonHxopwFnzxZtjGKFpXamL17Q8mSjh5dltIArJS6XeI+XFAQbNkChQtLgpWHB0yZ4ujRuYR7PUJ0p6goKQ42cvZGIp9uDI3iwIB8+Usz7fHvaV+5PRZDFjDrlaznGoF361bZ212wQK6xhg3lsHL79jn2Zk6XoJVStyxdKrOQHj2klNKECXD4cI79gLyb4J3hNBi7Gv/BP9Ng7GqCd4bffC61R4juFB8v/SgqVIBBk9eS0LYLuEnwtWChV+1edKjS4WbwdXrR0TB7Njz4oPz5/ntZTdm9G9atk7yBHHxt6QxYqZzMbocVK+C++6RoRoUK0KiRZJ02aaKVqpKRUpLVvXYRstulouKwYXA4ajP5nx4GvivJk6sQtphbjRKa+KdcgtIp/POPLDHPnCmZ8ZUry83cc89B/vyOHp3T0GNISuVEly7Jh2PiPtxLL8mRIpUq6a3TnOjHHeEMm3ieI7/4E8cxvFsPJbrEzxTOVZghjwyhb92+7Dy90zX2eG22Wxnyv/4qGfJt20rZ0cDAHHszp8eQlFK3vP02fPml7MM98sitfTiVahmRZPXx3HO8P9yLqDwHoEM3KPwXsWZuutz/NlPbvUsejzyAC+zxXrggN3NTpshxomLFJEO+d2/wc4GCHw6kAVip7C4mRkr2PfOM9Nj18ZE93n79pH6humfpSbLav1+Sfn9cEwWt+0LZX+UJ04pv3CCOHW1wM/g6LdOUBL1Jk2DhQulK1KgRjB0rGfLu7o4eoUvQAKxUdvXPPzIrmTlTjn74+soxonfecfTIXN7A5pVu2wOGlJOsTpyQSoqzFp3CrckoGDAdjBvfbwCmSbzlCKciamfy6NMhKkqymCdNklKRefLI9kXfvlCtmqNH53I0ACuV3Zw/L5mmiZWq/r0PpzLEvSRZXbgAY8bAxJkXSHh4LJY3gjCtCRSlJbbomlz0+BjTTMDADS979VQfVcpSR45IvsCsWdIcoWpV2evVsqPpogFYqezgwgXYt0+WAQsWlA9J3YfLVP/Odk5KZCS88VkIczf8Spx3GG4DvgdrFN1qdGf4o8PZc8yTIYtDcYsbRYwlFC97dXys1ZynU5HNBr/8IoF2xQpwc5Ncgf795bshCMQAACAASURBVAxvDk2qykgagJVyVYn7cEFBUjgjXz4ID5f9t40bHT06l5fWZgpxcVLWePjM1US0bgH14sGARmWaMPGJiVQuXBmAsjf6xo9b4cGpiMoZ3rAhzRJ7OU+ZIqVGixeXtfNeveS4msowGoCVckVr1sBbb0lN5sR9uH79NPklg6SlmYLNJp3zhr4Xw4lCU7G0GnqzS5HVsNKsbLObwTdRSrPoLGOasGmT7O1+953cRTRpAp9+Kk3v9brKFBqAlXIVhw+Dl5dUqrJaJfN00iRp16b7cBnqbs0U7gyYpgk//wxDhsaz130W7h1GgvdJAoo9wN7ze0mwJ+Bh9XDORglRUfDtt3Id7dwpqyh9+khSVeXKKX+/ShcNwEo5s4QE+XSfNAl++w1eeUXO8DZsCKGhug+XSVJ7znfKzyGMmb+aE+FxuD82D/L+zQN+DzOq6Rya+DchJCzEOYtoHDokSVWzZ0vHh+rVZcm5a1dZUVFZIsUAbBjGTKAVcM40zWo3HhsHtAbigL+BF0zTjMjMgSqV43z2GXz+OYSFQYkSMHIk9Owpz2ngzVQpnfPdswf6jt7IxgqNoXwcVIDSBSrweQvpUGTc+O/jVEU0Em/mgoLg999lWblDB9m6aNBArykHSE1F79lAizse+x2oZppmDeAQMCSDx6VUzmOa0jEmsTzsgQNw//1SROPoUSkUXKyYY8eYQyTXTOG5alXo1t2kZoflhJToeFujhBcCnufJik/eDL5O4+xZGD0aypaVI2n798OHH8rB5G+/lWpozjbmHCLFGbBpmusMwyhzx2O//evLTUCHjB2WUjnItWvwzTeyzLx3L2zeLJ1jJk+WvV6VKe6W5XznOd9C1nwUORxA//nbsAd2gK4bKJKrGJdjPLCZNjysHjQu09iRv87tTFMy4SdNkg5E8fHQrBl88QW0bi1HipTDZcR/hReBhck9aRhGb6A3QKlSpTLg7ZTKJi5cgBEjYO5cOTQaECANEapWlec1+Gaa1GQ5tw3wo3FZPz79FMbN30xM/XbQfSVFcxXnvcaTeTHgRbaf2u5ce7yRkTBvngTePXuk81C/fpJUVclJzherm1LVDenGDHhZ4h7wvx4fCtQB2pup+EHaDUnleHFxcrayQgXJQC1fXspD9u0rs15dCswSKXUziomRGDb8h7lcDxgDhQ9Q0LMwwx4dwst1Xsbb3cmqVR04IAOeM0f6ONesKQUzunSB3LkdPbocLVO6IRmG0QNJzmqamuCrVI524oR0Wp8+XZohHDgAuXLBsWM5uiF5ZkmpiEZyWc7hl2KYOROGjj/ImYD+8PgqANwt7nz3zHyalm2aJeNPlYQEWLJEAu/q1ZJU1bGjBN569fRmzgWkKQAbhtECGAQ8appmVMYOSalsZMsWSYD56Sf5+sknZbZrmvIBqcE3w6VmefnOLGfThOhDxbi6x5uXqr8AT8/FzWLFZhqYmNhNO1vCtzhHAD59Wm7kpk2TymclS0pLyZ49oUgRR49O3YMUs6ANw5gPhACVDMM4aRjGS8BEIC/wu2EYuwzDmJLJ41TKdVy8KGcrQZabQ0Jg8GDpTrR0KbRsKU0SVKa4WxGNRP/Oco457svpH/w5HzeT2K61cK89n9cffo3gzj/i5eaF1bA6vpCGacK6dfDss1CqlOQOVK0qM+CjR6XDlQZfl5OaLOjOSTw8IxPGopTrSjxCNGmStGsbOlSaIbRrJ6X8PD0dPcIcIzVFNNoG+PH3PneGT1tDVLmx0HkLFqudl2q/xPDAYZTIVwKAVc+tcmyS1Z0Z8j4+8Oqr8PLLULFi1o9HZSjNRVcqvWbPlupUiXWZX3hBusaAHPfQIx9ZKqUiGocOwaARVwiOHgyNZ4EhtZoXdFhAhyq3n6h0WCGNffsk6M6dK0E4MUO+c2fJHVDZgn4yKJUWJ09KdSqAxYsluzkoSOoy58vn2LHlcAObV7ptDxikiMaLAZV58eXrzN4/AbPeOPC+fNv3Hb54OKuHerv4eFlSDgqCtWslP+CZZySp6qGHNKkqG9IArFRqJX5ATp4sH5CHD0t1oa+/lqCrH5BO4c4iGkU881IsrDIvf/gtCQ+PhibnaFbqSboGdKDfz/2Is8U5do/31KlbSVWnTkHp0jBmjHS4KlzYMWNSWUIDsFIpuXQJJkyQD8jTpyUJZuRI2Y8DKXagnErbAD/y5znBmIUb+WPTNeIqdoVmJ6lXrDGfPhF8c1m5km8lx+zxmib88YfMdn/8UXoZtmghDRGeeEKLsOQQqSrEkVG0EIdyGXa7VKoqUgTOnZOg27ixHCF68kn9gHSwu53zjY+HIZP/5NMLTcEitZrL5a3KtHYTaOLfxLEDv3pVVkwmTZJ93gIFJGegb18pyqKynUwpxKFUtnTpkiRVTZkCRYvC+vUShE+ehEKFHD06RfLnfO12iD1UnP9NX8yZB/pBnjgALIaFF+t2cWzw3btXgu7XX0u5yAcegJkzoVMn8Hayqloqy2gAVgpg925p/bdgAcTEQP36ctQjsWCGBl+ncec5X9OEiwd96NJ/E7E1R8OjOyjsUZorNg9sdgc2SoiLg+BgWWZet06OonXqJElVdetm/XiU09EArHKuqCgpiOHlBX/+KV1jevSQ4vU1ajh6dDna3ZaYE8/zxlr2Exl5hOiD/thK/AQtN1DIrQwft5xN91pd2Rq+1TH7u+Hhki8wbRqcOQP+/jBunCw1+/pm3TiU09M9YJXzHDggS8yzZ8NHH0GfPhKMExL0CFEWSKlO851LzCDHiMa0r07bAD8ajF3NwRMHuFjojZt7vMQVwN+zBweGfoSH1QHlPU0T1qyRZebgYMkheOIJme02b66Vz3Iw3QNWyjRh0SI5QrRmjRSuf/pp2YsDLW6QRVJTp/lupSRr+/qRsDM3FwtNBKvs8WIaFLI+yWet38r64HvlihTLmDRJbux8feHNN2X7wt8/a8eiXI4GYJW9Xbkix4QMAz79VJYER4+GF1+UJCuVpe4WXO9cYv43W5Q7uzfZ8d/ZCXuVhVhtubGbVkxMLIY7bz7a8bZZdKbbs0eC7jffwPXrUihjzhwpnOHllXXjUC5NA7DKfux2+O03me2uWiUNEXx9pWJV0aJ6hMiBUlOn+d+lJO2xViL2unMt79fw6Ne44U2/Gu/wQYu3OHDhQNbu8cbFwQ8/SOD9808JtF26yBGiOkmuMCp1VxqAVfZx+bJUFJo6VToPFSkCAwbcer54cceNTQEp12kGKSU56Lu9hP99mOs+X0O97YAbgYV7sfCFDyiSW7r+ZFmd5rAwuaamT5cz4eXKyWrK889DwYKZ//4q29IArFybacq5yrx54fx5GDQIHn1U+qO2b6/9dp1McnWaBzavBEhBqFO7PTkeNpvYKgtuvMLK4AeDGPPES1k3ULtdmtwHBUkLSdOEVq0kqeqxxzSpSmUIDcDKNV29KvtvU6ZAhQqyNFixovRGLVPG0aNTybizTnNiFnSbWn4sCL7CgPmfcr7ceCh/XbKbAasB+fKey5oBRkRIdvzkydI2qVAhePttyZTX60plMA3AyrWEhsqsZN48mfkGBMjMJJF+SDpcSseM2gb43fb1b2uvU6bbWE6U/BiqXObhvB3o37QdvZf1zLpGCbt23bquoqOhXj2pWtWxo/ZyVplGA7ByfjExt/rqfv+9ZJt26iTJL3XrahciJ5KaY0YAIWEhfL1xJb+sucjxvPOh4jmqeT7BjK4f8mDJAADKFfTP3CSr2Fi5niZNgpAQKQnZtatcV7VrZ/z7KXUHLcShnNfhw5L8MmsWzJgBbdtKrWbQ5Bcn1WDs6iSTrPx8vNkwWGoxLwhZT9cVTbETDwYUNwL4uvNEmlSonzWDPH5cti6++koablSoINXPevSQ5ghKZSAtxKFch80mPXenTIHff5dZb7t20o0INPA6ubsdMwo/ZeP5T+ez0u01yBUPgAUL/Rt3yPzga7fL9RQUBD//LI+1bi2Bt1kzTapSDqEBWDmHqCipRmUYMHCg9JQbOVKakt93n6NHp1IpqWNGtlgr0eHHKDm6Jmbhv8iXUJ4YSyQ2MwsaJSR2t5o8GY4ckQb3gwdLUlXiTZ1SDqIBWDlO4qxk8mTYtAmOHZPiBr//DqVLa8EMF5R4zCjCtpcocy9xJ4sSk+d3qLqTvLGV+PCRhbzSpAObT27O3P3d7dtlb/fbbyWHoEEDeP99KT+qSVXKSWgAVlnv4kXZ000smFG4sJSGjI2VAFy2rKNHqNKobYAf+85vZ+iGoWDEQXlwiy3Gu7Vm807rrrhZ5CMnU4poxMTAd99J4N28WVZUuneXZeZatTL2vZTKABqAVdYwTTnekSuXFK0fNAgaNYIPP5SCGTorcXmmCR99s4X3dr0J+aRRgoGFES36M6xRj8x746NHJWdgxgy5uatUCb74Ap57Dnx8Mu99lUonDcAqc125cqtgRsOGMjupX1+CcKVKjh6dyiDTl+7h7eXvElFsKRbP/FhxB8OOh9WDpv5NM/4N7XZYsUKSqn75RXIH2rSRSlVNmujRNOUSNACrzLFz5609uKgoafvXoIE8ZxgafLOJRWsP8sqi9zhbaCFGgXy0yz+Smb1fY//FvZmzx3vxohxLmzxZti+KFoWhQ6F3byhZMuPeR6ksoAFYZZzoaNnDNQzZ3503T7rFvPyydovJZr5cuZjhK0cT4bkD8ueimecQ5r7+Fvf5yDnaerkyeI9361a5oVuwQPZ6GzaUtpLt2mm9b5Uxrl6VVZWOHbPsLbUQh0q/ffsk4M6ZA8uXSxm/M2ckGOsenFNJqUxkSrYfOkWHGa9yzHsxABbcWNA2mI61nsz4wUZHw8KFEni3boXcuW8lVVWvnvHvp3Kuv/6Sns7Xr8vfq1TJsB+thThUxouPlwYIkyfDunXg7g4dOkC+fPJ8sWKOHZ/6j9SWiUzK4fALdJ30EVuNieAlCVYYYMfO0kN/ZmwA/uefW0lVly5B5cowcaIE38TrS6n0uHpVtsdsNskbqFxZzoY/84z8PYtoAFb3JrFgRuKF6+MDH30kvVGLFHH06NRdjFtx8LY2gADR8TbGrTiYbAA+dekK3YM+Y3X0ePCIpOjFTpjulTifbzSmmYCBG+tDixBcIfyeZtL/YbPBr7/KbHf5cqlM1a6dXGOPPqpJVSr9TBO2bYNp02D+fJntPvaYXGMWi/R4zmIagFXKEhJg2TKZlfzzj2Qwe3lJ8Yxy5bSMn4u4W5nIO/12eDWvf/8x+yM3guc17ovqwKQn32fc+jOER0RTNG4UMZZQvOzVwV7xrkH8ri5cgJkz5do6elSqng0fLklVxYvf+89TKjlDh8KYMTKB6NxZrrG6dR06JA3AKnmnTsnd4ldfQXg4+PnJRRsXJwG4QgVHj1Ddg6TKRCY+nigqLoYOk4ew/PLn0o/Xw8KQyrMY/czzALzx01EAPO2V8bTfWqpLLrgnyTRlTzcoSPZ4Y2MhMFBWUtq2le0MpdIjcbY7dSq89prkDLRvL5nyXbpA/vyOHiGgAVjdyW6XD0Rvb6km9MEH0KKFfFg++aQ0R1AuKbFM5L+Xob3drQxsXom4hHgGfjubyX99QHyukzeft1oM8t53+ubXqQniyYqKkizmSZOkVGSePNCzp7T/q1o1fb+cUiB7u/PmycRh1y6Z7TZuLAG4Th2nO42ha4dKnD0LY8dC+fIyEwHpFvP331LooE0bDb4urm2AH2PaV8fPxxsDaRH4Ybsq7Dq6lvxDKzPhaG8s1/0ItIzBYniCacE03TDibgXHgc0r4e1+e43uxCCerCNH4M03oUQJaa4REyNB+NQpSa7S4KsyQnw8VKwoWfKmeesa69rV0SNLln6i5nR//CGZzIsXywUcGHjrLtHNDfz9HTo8lbHaBvhRtNAJ1hzbxLlzBi98N4Br3n/hFlWDPkWX0qRFACOW7aWI7cObe7xz1npwf0FJskrc503xKJPNJjduQUFyttLNTZYA+/eXM7yaVKXS6+pVSaZavx6+/lq2Lj7+GO6/X/Z2XeAaS/EcsGEYM4FWwDnTNKvdeKwgsBAoAxwDnjFN83JKb6bngJ1EZKQs/4F0h1mzRpqR9+kjF6/Ktjae2EjjOU2Is8WCAcbVUjzjM45pr3cgX14LDcauTnKJ2c/Hmw2Dm6T8BufPy/GhKVOk8X3x4nJd9eqlbSVVxkjc203MZK5ZE1atAl9fR48sSXc7B5yaJejZQIs7HhsMrDJNswKw6sbXypmZJmzcKIG2SBE4fFge//JLSbAaP16Dbza3ZNc6WnzVhTi7BF9MC0Me682Cd58hX175KLiXTOmbTBNCQqBbN1lmHjJEsuMXLZIWk8OHa/BVGSM4WGa3334Lzz4reSo7dzpt8E1JikvQpmmuMwyjzB0PtwECb/x9DrAWGJSB41IZ5fp1qVA1ZQqEhkLevHJmN7F8nx71yPbWHNpCz2+H8Y/xO8T4YsnljmGx4+HmQauqt89q7ynJ6vp1mYUEBUnCS758Mtvt2zdLixmobGz7dkmoqlVLrqvmzWVv14kymdMjrXvARU3TPA1gmuZpwzC0AoOzuXJFLtCYGPjf/6S02rRpcv4tcflZZUshYSGsPbaWorlKMP7XRexNWApRhah+9RPmvtKP6Py7km2UcLdM6ZsOHZIPwdmz5TqrXl1u8Lp21WtLpV9kpNzYTZ0qAdjb+1ZlPW9vCcTZRKYnYRmG0RvoDVCqVKnMfrucLTJSjnlMmSIJCSEhsjRz4ACUKePo0aksEBIWQpO5TYhJiAETiMtN2dMjmdH7NQLr5b3xquQbJSSbZFW9qCz/BQXBypVyfT39tCRVNWjgEgkvykV07QpLl0K1arJF1q1btq0pn6pmDDeWoJf9KwnrIBB4Y/Z7H7DWNM0U+8tpElYm2b9fjnN8/TVcuyYXbp8+ko6vVapyjGOXj9Ns+tP8HbX9xh6vQdeSQ/nmpZFp/6Fnz0ohlqlTISxM9nj79JHzu1rvW6XX9esyafjqK8kZ8PODLVuk+l69etnixi4zmjEsBXoAY2/8c0kaf45Kq+gb+3Te3tIMYcYMKST+8svZ5sJVqXP62mn6LxxFcNg0TDtgcQPsWAx3qpVJQ6k904QNG2SZedEiOZ7WrBl88YWcDdfz4Cq9du+Wm7pvvpFJQ+XKcPKkBOAHH3T06LJMiv8nGYYxH0m4KmQYxklgBBJ4vzMM4yXgBJB1DRRzugMHbrX++/BDmeV27y6diFw0E1ClzcWoi7z540d8fXAiduLxPPAC+RKewVJxN7Fue/5zhjdFkZFSRWjSJNizR3II+vWTPbdKKS5wKZU6587BAw/IjVzHjrKikkO3MVKTBd05maeaZvBYVHJMU5Zppk6Vwhnu7lLUILFgRq5c8kflCL///TvDfx/L1lObsFmi8TjUjTdqj2Bt9TDOXL8O3I9Xghwpi7bfvdsRIDd1kybJTd3Vq3Kucto0yTTNnTtrfimVfe3dK59d585J7e8iRWRlpWHDHD9p0LUkZ3b+PBQuLHeGQUFw+rSUi3zhBW39lw0F7wy/a4Wp63HX6Rc8kLn7Jsser2HlqYRvmDOlCz4+4D/4QJI/N8kzvAkJkugSFASrV8uxtI4dZcarWxgqvaKjJchOnSrbGR4eskpns4HVKk03lAZgpxMfD0uWyIW7fr0kvhQuDD/8IP/UpKpsKXhn+G3Hf8IjohmyOBSAltUL8fmfUxm5djTXjbM3v8dqhYebHb+ZIJqqM7ynT99KqgoPh1KlYPRoqdGsN3UqvUxTbt4mT5b63xUqwLhxUgCocOFkvy2lm8/sSj/NncXZs/DOO9Iuq2NHOWs5bNithJeiRTX4ZmPjVhy87ewtQFR8LP/76XOKja7A4D9e4/qxytQ9MwUvN2+shhUPqweBZQJvvj7ZRgmPV5REvU6dJOAOHy7nwoODpb/zkCEafFXaxcZKZapHH5UcApCAu2qVbG+89VaKwXfI4lDCI6IxuXXzGbwzPGvG70A6A3akhAS4dEk+/K5flzvFli0lk7l5c5niqBwhcZk41rKfaMseTLuNSPuf2M0TcPIh6l6bxdRBTQgIMAgJq5FkIY07z/CW9zYZF7uTWj0Gyj6cjw+8+qokVWkvZ5Vehw9LrsDs2XDhApQte+szy9cXmqSidjhJ33xGx6cid+EeOeMsWwOwI5w4IcuAM2ZIItWSJXLxnj4NhQo5enTKAYr7ePP31R2c9RgCJMge78VyeP85n18+fJbAQNmTlQ+RaE5F1GCZTzQDm9+e4dw2wI+2nlckqWrSXDniERAg11vnzpqsp9IncYkZZKVu717Zz+3TB5o2TdMqXZrqj9+ju23xODIIawDOSqtXw2efwfLlciG3aCEFDRJp8M2RTNOkacApNq3/GIwEedBuIW9Cc+ZMaEhg7VvBN9kPkWpF5EYuKAjWrpWkl2eekUpVDz2kSVUqfY4ehenT4fvvYccOqSk/Y4bUkk9no417qj+eRlk1y75XGoAzW3i4LDG7u8Off0pt0yFDJPBqecgcb93xdQxYMpTdl/+E6GJgcQeLHYvFjXeebk672rc+HJL6EMl76RwXBg6F/b9L8/HSpSVT/sUX77rvplSKEhJg2TIpbfvbb3IT16qVbJvlzStneTNAquqPp1NWzLLTQgNwZrDZpAn51KlyAS9cKCn4b74pwdfd3dEjVJkgNXtMiY0SiuQuwqyt37HhzG9w7T5ybw/i/TY9qfPUdjaGJ90o4eaHhWnycFgo3Xf8TPNDIbiZdllNmTpVcgg0d0Clh90uS8n79kG7dlKdavhwyZQvWTLD3y7Z+uMZODPNill2WmgAzkixsZJI9dVX0oy8SBF4+23pXwla1CAbS80eU0hYCE3mNCHGFiPfFJMPj02fMDCwH4OWepM3L0A9Hi2bdKOE8l52Ht74C913/ELFiye47JWXGXXbsqpRO74b1y2zf0WVndls8OuvMtv19ZXEqho1ZNusYcNMLz/aNsAvU5eCs2KWnRYagNPLbpdswEqVZN/tm2+kGfm4cdCmza2+uypbS2mP6dDFQ/Rc8rJ0KTIAu8FD5hv89N2bKa8U790Lkybxy5y5uEddZ3exCgxs+RpLKzfCkisXY9pXz7TfS2Vzp0/LXu706ZIcWrSo5A0katzYcWPLQFkxy04LDcBpdeYMzJolF+7ly7L/5u0te7w6081xkttLOnHlOD1+eIlv9s7BnuAGhhsYJobhRtvHaicffOPi4McfJZt53Trw9MS9UyfWNuvI0JPeTvUholyM3S7/tFgkKfSTT6TZxqefyqQhm26RZfYsOy00AN+r0FB4/33JOE1IgMBAqSSUuO+mwTdHunOPKYFLXLF+R6TbCubuNmDrK+Q61RPv+uuw+W5OvlHCyZNyUzdtmtzklS0rqykvvAC+vgQCG7L8t1PZwvnzsrQ8dapkyzdvDm+8Ab1767lwB9EAnBrnzsmMpEQJiIqSYx6vvQa9emmXGBeR2YfwE/eYLto2c8Xte+I4DNhhx4s8GDuM6HJnuVrzOFAaEkoD/2qUUKs4rFkjH4pLlsgM5YknZCmweXOtgKbSzjSlpO2UKVLONi5O9nS9vOT54sUdO74cTgNwcux2+VCcOlVK9vXoITOTBx+Uo0Weno4eoUqlrDiE37hyHor7z+LAie/kAdNKoQ1fM/fNzrRoAWWH7PnP9+SNvc7jK5fCnP5Sss/XVzLl+/SRma9SaZWQIIlTNpt0tbp+XSqg9e4tZUiVU9AAnJRJk2D8eDhyBAoWhFdekQsX5CycBl+XkpmH8K/HXWfilomMXvcxV+MvyYMGWAx4Y/gxWjaSh/69RF353D903/ELbfetIVd8rBTKmDNHCmckzkyUulemCZs332rksn+/7OcuWwYVK2oVNCekARjkwt2w4VZT6P37pbrLe+/B00/rh6KLy4xD+LEJsUzbPo3314ziYuxZONySfGEdiW7SH7sRh4fVg8b+gTdf/3YTf9aNmUqnrT9RN3wfMW4e/Fw1kIIDX6dx15ZpHodSXLsmTRCmTIHduyFPHujaFSIjoUABqFXL0SNUycjZATgxKWHaNJnt/vEHNGoks99MPvemsk5GHsJff2I940PGs/7YBi7EnINjj5J78w8M7d6AAVNhz6X7b2+UcOIETJtGm+nTaXPuHCcLFmdU4xdZ/0hrXm5Xl8ZOlpWpXEh8vMxwN22S5eVatSQId+nCjUPlysnlzChz7hwMGACLF8tF/MgjUuklsWCGBt9sJSMO4dtNO10WDGThwfFgmGC3YF35Of97ZACD/zQoWFBeVy93Per5PSSt2F5tCz/9JCssrVpB//6UeOwxhmpSlUqrqCj47jsJtA0ayNGhpk1h61YpDak1v11Kzok0Fy7A33/LfpuPD+zaBf36SSZz1aqOHp3KROk5hG+aJsEHghnw8zucvH7gX88Y5H9kL/U7naJgwRs/JyJCVlQmT5Z+zoUKSSW0Pn207rdKn/37ZW93zhy5zu6//1YylcUiXdWUyzFM08yyN6tTp465bdu2LHs/TFOODE2bJrPd4sUlCFsst7fVUuoOpmny29+/8c6qYew4sw0uVYCd3eHR0WCNw8CNonGjKJuvNhtaFJQjRPPmQXQ0PPywHCHq2FET9lTaJS4xg5zCmD9fasq//LIcJXKSzy9n7LPrTAzD2G6aZpJ3SNl3BhwcLLOPw4dlxpuYgp+4/OckF69yHomNEvJ75Wd+6AL+DFuP9VppWDUTzystKNDwb7B5E2OGkieuMk/vP0P3Hf1hyAGpgta1q1xntWs7+ldRruyff2TSMHOmNHUJCICRI6UgS5Eijh7dbZy1z66ryD4B2DTl3G6lStK9w2qVi/Xdd+Wu0duxXS+UcwsJC6HxnMbE2mIBsMb6wsogahk9+egTD4ZvXU14RDR+l33psiuaZ/eMpVDUFU74+knSXo8eknGqVFr8u/XfihXy+dW69a18Pj9OCAAAIABJREFUlFKlHDu+ZDhrn11X4foB+Nw52RdJzGQeMUKOD7VqJRewUikIPRtKz5963gy+2A18Dr7K1EH9aN8eDNOO94HTRM+aSODhLQCsKv8gC+u2pvX/ulPqgYxv0aZyiLg4adgSHQ3PPQf58snn10svSeW9LJCeJWRn7bPrKlw3AJumXLALF97KZB4xQs7tgi4xqxQdvniYEWtHsGDvAiwJ3jcbJbhbPfjx08dpmO8SjJ8NkydT/8gRYgv48s2jnZhWuRmUKq17XSpt7HZYuVKS9Y4dgx075NjQn39KYlUWnsJI7xKys/bZdRWuG4ANQxof9O8vmcxaXk2l0okrJ/jgjw+YvWs22DwxNwwi776BdHv9IIXrruUxWzHqjfwKvv0WYmLkuMf77+P59NP08PSkh6N/AeWaLlyQfd2pU2Wft1AhmenGxUmyXo0aWT6k9C4hO2ufXVfhugEYZL9EqVRadmgZo9ePZmv4Vux2C+aWV/DaOoQ3Xy7KW/NiyL/iMAxZIuX8cuWSFZbEAgdKpYVpyv6uu7vs7Q4aJMV+PvwQ2rd3eJZ8epeQnbXPrqtw7QCs1A1328e6GHWR1359jXmh88AETCuWRQvp17gtI947SuEfBkHVGXDxoiTxffGFBF8fH8f+Usp1XbkCX38tk4Tu3SXwPv203Mw5Ud2BjFhCdsY+u65CA7ByecntY0XFX+NQ1AI+C/mMa3HXJPgaYACDegUzeuVX8PAvcjStTRspzNKkieYPqLTbsUP2dr/9VqpW1alzq9eul5dTBV/QJWRH0wCsXN6d+1h2YjhjLqPH8sUkcBWvo+1he3us7XqCJRYPm53Wn86B+KIwbJicD8+ijFOVDf27YMbQoVJTvksXKZjh5BWqdAnZsTQAK4dLbyWdxP2qGEsoV6w/Emvdh2lEYhx7HFaMpkshO6NKTuLoPDtri5sE+tSi3ieDoV07OQKSAWNQOdD+/bLEPG+ezHxLlYKJE6WvswttX+gSsuNoAFYOlRGVdIrld+ev69O55rZI1pftFtyWTeDFw+6MLdqXAqFb4WgeinV/iXp9+0L16hk+BpVDxMVJlb3Jk6XMrbu77O3Gx8vz5co5dHjKtWhbFuVQdzsGkRK7aefb0G857t6Ha+6Lbj5umCaD8wxk6tW+FHCLlFlJeDhMmvSf4JveMagcIjHAXrwoy8vHjsGYMXDypNRo1sCr0kBnwMqh0nIMwjRNlhxcwrtr3mXvub3kj65BpfUvcLTxHGwWOx52k7oFH4A1o+DRR1NMqtJqPipJNhv8+qvMdmNj4fff4b77YMsWyWbWtpIqnTQAK4dK7TGIkLAQ1hxbg7ebN9/u/ZZtp7bhE1eO1ku7MP6vjZQzZ7Hhui/r2tYksN1r1PvwqQwfg8ohzp2DGTOkYMbx41C0qBT7sdsl6GqzDZVB0hWADcN4A+iJHPAIBV4wTTMmIwamcobUHIMYu2op7/zZAdOMBwNyxRWg1/L6fLZ7G3nsfxPXIBBeH0eDNm1okJiNmsFjUNmcaUqAtVolqeqdd6BxY+lA1KbNzWQ9pTJSmgOwYRh+wACgimma0YZhfAd0AmZn0NhUDpDSMYhxq39m+Lq+mBYJvoYdBv95mf+FhmLv1gsG9cMjnWVI9ShGDnblCsydK9nMAwfC88/Diy9Cy5bS9N7Jafa+a0vvErQb4G0YRjyQCziV/iGpnCapYxB7z+3l3TXvEnwgGK94d0w3WWZxtxscyv0MbYd14/cRrTJ1DCobu7NgRt26ULiwPJc/v/xxcpq97/rSHIBN0ww3DOMT4AQQDfxmmuZvGTYylW2k5i49JCyEtcfWUq5AOYIP/MiCvxaSJ97K+39C/012JpdoxPyG+blUrBHra1TG0Pwoda9sNlliBujZEw4ckIzmvn3hgQccO7Y00F68ri89S9AFgDaAPxABfG8YRjfTNL+543W9gd4ApZy0qbTKPKm5Sw8JC6HJ3CbEJsRgAp4JMCgEumwsTLDbiwQ+Vo+rVSSZ2dOUn6sJUirVDh2SJebvv4e9e2V2O3euVD9zoYIZd9LsfdeXnjz6ZsBR0zTPm6YZDywG6t/5ItM0p5mmWcc0zTqFE5d4VI6R0hnbM9dO878fehETL8HXMKH9xtJc+2sR+6cfp+qaviTUst52kkgTpFSK4uPhhx+gWTNpsPHll1C/Ply9Ks9Xq+bSwReSvwnVm1PXkZ494BPAw4Zh5EKWoJsC2zJkVCrbSO5u/Pz5Ewz6tAVfRvxGrMXECthNA8PmSYk28xn1a70b5XX9sFg0QUqlUuJRoUOHoEMHKFlSWv+99BIUK+bo0WUozd53fenZA95sGMYiYAeQAOwEpmXUwFT2cOcZ25KXj1Dw2mTWlj7IuGvQ5rAPXr//j5WF61O72xYGdQukSYV6t/0MTZBSd2W3w6pVklSVO7e0AaxaFdatk1lv4r5vNqPZ+67PME0zy96sTp065rZtOknOSYJ3hjNs0S7KnJjPea+l7C907f/t3Xd8VMX6x/HPJBCKiDTpIqIIiCBqQCOKEa40kSL6Q0UELFzFzr0gXEVBUVH0SjGRjgW7ImAHgaA06UUEBOlBQAjXIAlpO78/ZkOJSQghycluvu/Xi1eSPWeTGc7ueXbmzHkeDpeEFttrkPDti/x8qAePPW4YMAAqVPC6tRJQ4uLgrbfc9d3Nm6FSJVeB6PnnvW6ZyDHGmBXW2kzLYikTluSfffto9+UE5v4ygjGXJfjv4zVU+PI1Fq56gvvug+mDoYY+sEtOpQ8YjIFXX3X5mJs3h2efdVPOJUp42z6R06BkppK3rIUFC0i983YmdajBxXFDGdMk4YTNIdRodpRffnEDFwVfyZGEBJceMjwcvvnGPfbII7BmDSxYAN27K/hKwFEAlrzx118wbhy+JpfxwYPXcUnlT7ivQxrFyzfinMUjIbUU+EIpUTyMcQMjufhirxssAWHjRnj8cahe3d27m5x8fFu1atC4sXdtEzlDmoKWM7NxI0RHs2jWJN68JIFFrcPYWgbOD2tAjTkv8NucjjRrZujRtBmHK8QQWTuSiPMiTv17peiy1k0x+3zQti3s2eOmlx98EK699pTVrUQChQKwZCvTLFaNqsDMmRAVhZ07l1HXhNDvdh/WgCGFGpuGsuPDp6lfL4Rp06BzZzAmAlDglWzs3g0TJsAXX8CSJa4Awvvvu1q7Vap43TqRPKcALFnKmMUqeXcsO594i8QNsym1fy8/Nq3MU0Nr86Pdfuw5Ni2Ew/HFmTQxhLvvhmJ6hUl2fD6YOxeio92HOp/PFUI4cMBNO1/zt9w+IkFDp0fJ0ojvNpGYnEqz3evpsfIr2v66iOK+NN66qgEfPlmN7w6v4twSVbk4th+/ln0TQpMpHhrGjFGRRF7odeslICxaBDfe6G4h+ve/4Z//hAsu8LpVIgVCAVgyd/gwN8z7jLtWfsWhkjv49qLiTK/XhJn10zhQajXlkytydfwIlr7Ul4Ripen1xK3UujaGtvV1jVeysXy5G+1WrgzDh7tbiD79FG66CUqW9Lp1IgVKAVhO9ssv7gT5zjsMO3yYqY1r0KtzKGkmBVgBlODsbU+Q8PkQViSWpe8/4emnoUoVXeOVLCQmwkcfudfVsmVQujT07eu2GQNdu3rbPhGPKACLS1w/YwZERUFMjLufsls3Pr7+eu7b8RJpJg0MYA0s6sdf379A9+6GoUOhTh2vGy+F3r/+5dJENmjgiiL06BEQ9XZF8pvuAy7K9uyBoUOhdm247TbYvh1efpm9m1bw6O3n0GPPg6SGbAdCwRcKqSWpH9aa1asN776r4CuZSE2F6dOhTRs33Qzw6KMwbx6sXw8PP6zgK+KnEXBRYy3Mn++mAz//3BUpb9sWxo3j4PXNeGXJa4yZ2pTktGRalO3N1imD2REXy/nXxzDwjkgeGKZpZsnE3r0wcSKMG+duJ6pZ0z0GUL+++yciJ1EALiri412VmOhod523QgWXYeiBB/jebGP4guEsHrOYxJREIivdyb4PhzBvwUU0agRfjq1F+/YRyn8gmUtJcRmp/vjDrWgeMwY6dNA9aCKnoHdIsPv5Zxd0333XpYsMD4cpU6BbNxKKWfrP6k/08mgAQgil3i/vMO/ju6hd2z3ljjtOXc0t02QdKokWvOLjYepU+P57V/S+eHEYPx4uuQTlGBXJOQXgYJSc7KaXo6NdTdQSJeD22+Ghh6BpU5JSk5iwcgIv/PgCe//ae+xpvjTYFb+LMWOgTx+XiOhUMibriP1fIoOmrQNQEA42GT/MXXGFG/VWruzSnYnIaVEADmAZR56Dm5xN20VfuHR+e/e6VVIjRkDv3lCxIqm+VN5eOYnnfniOnX/u5KoqLai9/RmWnPMvCEmmmD+JRqvTGMSM+G7TseCbLjEljRHfbVIADiazZ0Pr1sc/zPXtC02bKi+zyBlQAA5Qx0aeyalcs2MNPVZ9zT82L8FiMe3bu9FumzYQEsLCnQsZM28MC3cuZPfh3VxeuSnhRybwxWM3AlDxxqqkXvIFtcpdzuEjtU6rHXv+l3haj0uA2L3bTStXq+aKIFx/Pbz+uruFqGJFr1snEhQUgAPUm9NX0G3xV/RY9TUXxu0mrlRZJjbrwvctuvDpK90BsNbyyoKXGTRnkAvMGG5kOIufHMCaI4Ybbj7C9lrLSC0dBnQlPp7Tnj6uXq4UsZkE2+rlSuVZX6WAnJiXecYMt2K+Tx+3LSzMLdoTkTyjAOyhXC1eWrMGoqP5fMrblE5JYlW1evyr/RN82eA6koqFuXwZ1jJ762yenvs0y/YsO/ZU6wth9lwft9xoGDYM7pvxE6kZgufpTh/3b1PvpGvAAKWKh9K/Tb0c/z9IIfHgg27UW6kS9O8PDzzg7hEXkXyhAOyR01q8lJTkVptGR8PChVCyJHMb38DYS9vyc9WLTtq1TNktRL4dyQ87fqBW2Vp0OGcQXx0ciTXJhBDGuKciua+t23fP22c+fZzeVq2CDkCrVrnX1FNPuUDbsydcd52rvau8zCL5TgHYIzlavLRzpxuRTJgA+/e7uqivvQa9epG6I5Hfpq0jKe1njoasI8RXnqTiC9iRsoKqB6vSp+YYFoy+ny/XluDiVjfTvHsM97WK5JpaxxNp5NX0cefLayjgBoqjR13xg6goV3O3dGlX/q92bVf6T+X/RAqMArBHshpl/n7oiFtxGhXlCpNb6yrFPPSQW4Ua4rKHdq4AG+NW8J8F/8HaFDBQKrQMD9Z+hZXjHmL8gtLUrety4N96awQhIX/PYKXp4yImKcl9iNuzB+rWdYuqevWCcuW8bplIkaQA7JGMo8+yR//i1nVz6LX2G3hlt7sON2CAq4+ayXW4LXFb+HDzM1hc8DUYauzox5uD+1O9ussI2Ls3fPVzLNe9kvn0sKaPg5zPB7NmuXvBX3zR3UL073/DpZdCq1bHPsyJiDeMtbbA/lh4eLhdnp6gvYhLvwZ8we7N9Fj5JZ1/mU+p1CQONr6Siv0fz/I63M4/d/L8/OeZsnoKoSGhpPl8pPkspIZRZtocBveK4OGH3cxixuvM4Ea4L93SSEE2mMXFuWxnb74Jv/0GVaocTz8qIgXKGLPCWhue2TaNgL2QlETn9fO4buZIKq5dQWKxEsy+vBVnP/4IN9zZ9qRdF+9aTMz2GBpXacys32YxdsVYAHo1fIi0mEG8+8U2itWJ4Y6rIxm1MILy5Y8/V0kyiqDvv4ebb3bXeq+9FoYNg1tuyVlaMxEpUArABWnHDhg7FiZNgj/+oKL/Olypnj3peGLk9Fu8azGt3mnF0dSjWCwhhHDXpfdQdtVgJt9di6QkuP/+qgweHEH16n//c0qSUQQcPQoffwzly7vA27Qp3HOPu3TRuLHXrRORbCgA5zef7/iiqq++co917OhS+WVzHS4+KZ7B8waTmOqCpcFwbbHH+bLPa8TFuWyAzz3n1tJkRUkygti2bcc/zB086C5Z3Hyzq7UbFeV160QkB7QKI78cOgT//S/Uq+fq7f70Ewwa5E6cn3/uyrZlEnwTUhIYsXAEdUbVYc62OYSYUAyhkFqSH8beSrNmsHIlfPBB9sEX3CrnUsVPLmWkVc5BYMAAt5r51VehRQs37fzxx163SkROk0bAeW3FCpfc4IMPIDERmjd3Q9WuXbO9DpexQlHbC9sSyfNEv5nCzpAYGpaJ5I13I4iMzHlTtMo5SMTFwVtvuWXt5ctDs2YueUafPnDeeV63TkRySaug80L6dbjoaDfSPessuOsul9rvssuyfeqPO39k5JKRLNy5kH1H9tHi/BZ0KTuM94dfx7Jl0KCBu4OkUycVnilyVqxw08kffOBeY1OnQvfuXrdKRE6DVkHnl23b3A23kybBgQNQvz6MHg133+2uxWXDZ308P/95hs4feqxQQr9LXmfthMd4YrbhvPPgkSGHWF58FU8sSWTERo1ei4zERGjZ8nimqrvvdmsGTvFhTkQCiwLw6fL54Lvv3Gj3q6/cddxOnVymqhtuOOUw1VrLzE0zGTxvMOv2rzvh8RD+OzqRihsMr78O1a+O5dmv1pGYqEL3RcKOHbBggRvhlioFjRq5lXY9eypTlUiQ0iKsnDp40C16qVsX2reHZcvg6adh+3ZXKKFly2yDr7WWWb/N4qqJV9H5o84cTT1KvyuGEOorBWmhkBpG75aRbN3qqr6Nisn6Hl4JEumZqjp1gjp14N574X//c9vGj4fHHlPwFQliGgGfyvLl7jrchx+663AtWriLsl265Ci5weJdi3lrzVssjV3K6r2rqXVOLUa3msSOmXcTNaAYVG1NxG0xPHVXJDc1Pp6vWffwBrlFi1we5s2b4dxz4cknXfk/BVyRIuOMArAxphwwEbgUsMA91trFedEwTyUmukVVUVFupFumjFuB+uCDbmowh6asmsJ9X9yHz/oAePjKfpy75kUGdyhBfDz06AFDh0ZQu/bfCyXoHt4gtGaNG/VefjnUrAmVK8Mzz8Btt7k8zSJSpJzpCHgU8K219lZjTBhQOg/a5J2tW13+3MmT3a0fDRrAmDFuEUzZsjn+Nev3r+eZmGeYtmHa8QdtKBPHVODorBJ07OgyBGYXy1WpKEgkJ8O0ae7D3IIF0KGDq3JVq5b7WUSKrFwHYGNMWaAF0AvAWpsMJOdNswpQWhp8+607QX77rVtU1aWLW1R1/fWnde/PlrgtDIkZwvvr3qdksbMom9aGeBMDpEJaGL6D4bw05Q8G9jr3lL9L9/AGgagoeP552LfPXeN99VU3kyIiwpmNgOsAfwBTjDGXASuAx6y1R/KkZfntwAF3+9DYsW4hVbVqbjrw/vuhRs6CXHqhhPqV6vP15q+ZsnoKYaFh/Pua/kz/sB3bv7wcwn4hpOG3nFO5Fme3SuWLvesYSMsc/X4Vug8w1rrSf1df7aaUExMhPNx9mGvTRuX/ROQkuU7EYYwJB5YAza21PxljRgHx1trBGfbrA/QBqFWr1pU7duw4wyafAWtdoozoaHeNNynJjXIfegg6d4bixXP8qxbvWkzLd1qSlJqExVIspBgPhj9I61L/YcSzVfnhByhW7gjlrvuV0g32HBtIG2Db8Jvyp3/ijcOHXZKMqChYvx7eew/uvNO93pQ9RaRIy69EHLuB3dban/w/fwoMzLiTtXY8MB5cJqwz+Hu5l5DgsglFR7tEymefDffd5xZVNWx42r8uLjGOgXMGcjT1KOAKJdx9UT92jnuZm2e48qsXdNxE2sW/YUJP7rIWUQWRxES3evmtt1wQvuIKN6vSubPbruArItnIdQC21u41xuwyxtSz1m4CWgG/5F3T8sDmzW5R1ZQp7v7Khg1dEL7rLheET1N8Ujwjl4zktcWvEZ8UT6gJdYMcXxiTn+xM2T/d4qrHHoPvN5dh0LQQLaIKNqmpsGGDW0FXsqS7nSg9EctVVynoikiOnekq6EeA9/wroLcC3q8wSUtzGaqiolySg2LFXCGEvn3huutydYJMSEkgamkULy98mYOJB+lSvwuPXTaU6Il/8dnyGEJ2RfJEtwgGDoSKFd1ztIgqyPzxB0yc6NYMHDoEsbHuQ9xPP0Fo6KmfLyKSwRkFYGvtaiDTue0Ct3//8UVVO3e6hVTPPeemmqtVy9WvnL99Pq8tfo2FuxYSlxhHmwvbMOiqYcS8H06H+9zMdq9eEQx5N/OiNFpEFQR+/dVNa3z0kbulqGVLePhhly4SFHxFJNcCOxOWtS5hfVQUfPKJO0G2agWvv+6K3hfLunvTV8VmOTpN9aXyzLxnGL5gOBZLiAlhdOtofEsf5NZr3ALqrl3debl+/YLqrBSYo0chPt4lykhIgOnT3er4vn3hkku8bp2IBInADcDTp8PQobB6tUuS8cADblFVDiLi9FWxJyW5SC904LM+ksMW8mzMs/x68Ndj+/t8hkHPx3HkGxffX3oJYovFcu90TS8HlR073JqBiRPhppvg7behSRP4/XdXYlJEJA8FbgDes8el9Rs3zt3yUaZMjp864ruTCx1YLAfTFtLjy74k2G3UOrselVJ7cyDkfdKTaKT83owhbx7g2QcqZRnAQZWKAtL8+W7W5Isv3M+dOp2cMEPBV0TyQUBmBpi+Kpbr/leXC9q+QPO4i5i++c/Ten56QYOjIRs4UHwke0r05Y8Sw0hOS+L9W96nzKapHJ78GkyeR8iSpzlnaxRVWycz639rgb8HcFClooATH+8uYQDMnOlWMw8c6Go8T5sGkZGeNk9Egl/AjYBPGn0ak6vRZ/VypdhweBqHikeDsWChbMpt1I5/nHefvIZfvoHQMolUuPYsyjRqhglxJ+r0wK1KRQFs/Xq3ZuCdd9xljH/8AwYPdhWuVBBBRApQwAXg7EafOQnAK/asIPGc5zmUFOPqNwEQQtLGC1n72TXsKg/nt9uMr8EWQor7TnpuehINVSoKMGlpMGMGvPEGzJvnAu0ddxxPOaoSgCLigYCbgs7t6HP9/vV0/bgr4RPC2XF4Le1q98RQAnyhkFoC36oO/Oc/riDSyBdKc1bpk+8XPjGJRv829ShVPDTL7VJIpKS4rz4fPPKIO7jDh8Pu3S45S4MG3rZPRIq0gBsBn87oc/GuxXy24TPW7VvH7K2zKRNWhiHXD6F3gycYN7oscz/+JynVY+jcJJI3FkQcu134VEk0lGSjELMWli51o91Fi2DjRpfjOybGVSTSfbsiUkgEXADOaZ3c6Rumc+snt5Jm3X7dG3VneOQoPphUkSa3uGRGd9wRwXPPRXDRRX//O6dKoqEkG4XM0aMuWcYbb8Dy5S5LVa9e7j7ec86BunW9bqGIyEkCLgCfavS57699vPjji0QtizoWfENNKAnbG3J144rExkK7dm7NTZMmnnVD8kp6xaEffnABt0EDt8iqR49c5fsWESkouS5HmBvh4eF2+fLl+fK74xLjGLFwBKOXjiYpNYl2ddvx/dbvSU5NwaaGYd+aQ8R5Ebz0kqtAKAHMWpg7141269Vz13V9PliwINf5vkVE8kN+lSMsFA4nHWbkkpG8uvhVDicd5vZLb2dI5BB2rrqYLZMXs/FoDBeYSEZGRXDzzTo3B7TDh93tQ1FRriJRpUquAhG4YvctWnjbPhGR0xCwAThmWwyvLHqFhTsXEp8cT+f6nXku8jmO7mxE324wZw6ElbucSs1LUfWaOHznxWKMrtkGtMcfh8mT4corXQ3ebt1cSUARkQAUkAF48a7F3Dj1RlJ9qYSYECbePJHmpe/l6Yfgs8+gbPk0Krf+lZKNtmOK+dgTj1JFBpr0spJvvAEjRsBll0H//tCnDzRrpqkMEQl4ARmAY7bH4LMuSYbBMPbd/fQZDaVLw5Ah8A0L2Jv410nPOZ1kHeKhuDhXVjI6GrZvd8kyYmNdAFbpKREJIgGXiAMgsnYkYaElMDaUtOQwVs+I5NFHXZ6FZ5+FfRmCbzqliizkUlLcKuYBA+D8812JyW3boH17r1smIpLnAnIEHHFeBB3j5vDJshjaXxJJ9LwIzj//+HaligwQKSmu8MF337lRb/HiMHIkNGwIjRt73ToRkXwVkAEYYPSACJ45EEHDhn/fltNkHeKRvXth/HgYO9bV2q1TB/btg6pVXY5mEZEiIGADcJUq7l9mlCqyEPvxR2jVyo1+27aFCRNcZpSQgLwaIiKSawEbgE9FqSILiaNH4cMPISwM7rzTrWDu1w/uvVfpIUWkSNOwQ/LHzp0waBDUrAm9e8N777nHS5RwmasUfEWkiFMAlrw3bBhccAG88opLDTlnDnz5pdetEhEpVBSA5cwdOQLjxrn7dcFNM/fv7+4L+/xzaNlSiTNERDJQAJbc++03dz23Rg144AGXhgygdWs3zXzivWEiInKSoF2EJfnI54OuXWHGDFfg/tZb4eGH4ZprvG6ZiEjA0AhYciY+/vgINyQEqleHp5+GHTvggw+geXNNM4uInAaNgCV7Gze6gghvvw1//eWu615wgSsJKCIiuaYRsGRuyxZ3LbdBA5cso0sXWLrUBV8RETljGgHLcYcOwZ49LhdzhQquGtGwYXD//VC5stetExEJKgrAAj//DGPGwNSpLvguXeoC8KZNuq4rIpJPNAVdlM2dCzfcAI0awTvvuEII48cf367gKyKSbzQCLmoOHIBSpeCss9x9vFu3wssvu9zMFSt63ToRkSJDI+CiYvVqF2TPOw8mT3aP9erlgvCAAQq+IiIFTCPgYGYtfPopjB4NCxZA6dJw993wj3+47cWLe9s+EZEiTAE4GCUkuGBrjFtcFRsLr74K99wD5ct73ToRESEPArAxJhRYDsRaazuceZMk11ascAF32jTYvBmqVIGPP4Zzz3UpI0VEpNDIixHwY8AGoGwe/C45XSkpLuCOHg2LFrnFVT17unzNAFWrets+ERHJ1BktwjLG1ARuAibmTXMkx6x1X3ftcrcP7dsHI0e66eaoKKifMwaOAAAJr0lEQVRWzdv2iYhIts50BDwSGACcndUOxpg+QB+AWrVqneGfE5Ytc9PMCQlugVWdOvDTT3Dlla5IgoiIBIRcn7GNMR2A/dbaFdntZ60db60Nt9aGn3vuubn9c0VbcrKrOBQR4Yrdf/65q8GbPs3ctKmCr4hIgDmTs3ZzoKMxZjvwIdDSGDM1T1olJxs1Cu68Ew4edN/HxrqvCroiIgEr12dwa+0ga21Na21t4HZgrrX2rjxrWVG2bBn06OFGuuASZnz9tSsN+OijUFbr3UREAp2GUIVFcjK8/z5cfbWbZp4xw410wd1G1K6dRrwiIkEkTxJxWGtjgJi8+F1FVuvWMH8+XHyxu6WoZ0+NdEVEgpiGVF5ZtszV2T1yxP3cvz988w1s2ACPPKLgKyIS5BSAC1LG1cwffeSKJADcdBO0batpZhGRIkK5oAvKvn1w+eXw++9Qt66mmUVEijgF4Py0fDmsWwe9e7u8zN26uWu9bdpopCsiUsQpAOe1lBT47DM3wl282K1gvvNOKFECXn/d69aJiEghoWFYXvruO6hd2+Vm/uMPlyxjyxYXfEVERE6gEfCZWrnSBdiGDV3wbdQIJkzQgioREcmWIkRupKS4OrvXXuuKIAwb5h6vVw++/Rbat1fwFRGRbClKnK7oaFeBqFs3t6L5v/+FsWO9bpWIiAQYTUHnxLp1cMklEBoKe/dC/fouELdv7x4TERE5TQrAWUlNhZkz3Wrm+fNh+nTo1AmGDNH0soiInDFFkoySkmDECLjwQujaFbZvdz+3aOG2K/iKiEge0Ag43aFDUL48FCvmrunWqeNuI7r5Zk0zi4hInivaATgtzdXZHTUK1q6FnTuhZElYsQLKlfO6dSIiEsSK5nzqn3/CyJGu9F/HjrBpE/Tr5wIyKPiKiEi+K1ojYJ/PXcNduxaeeMLdxzt8OHTp4qaeRURECkjwRx2fD2bNctPMF10EY8a4wLt2rctaJSIi4oHgnYI+fBiiotz9u+3aubq755/vthmj4CsiIp4K3hHwwIEuWUbTpjB1Ktx2G4SFed0qERERIFhGwNbC3LnQuTMsWeIe69cPFi2CpUuhe3cFXxERKVQCOwAnJMDEiXDZZdCqFSxcCLt2uW0XXggREd62T0REJAuBOwXt80GTJrB5MzRuDJMmuTq8pUp53TIREZFTCtwAHBICzz4LNWu6NJHGeN0iERGRHAvcAAzu2q6IiEgACuxrwCIiIgFKAVhERMQDCsAiIiIeUAAWERHxgAKwiIiIBxSARUREPKAALCIi4gEFYBEREQ8oAIuIiHhAAVhERMQDCsAiIiIeUAAWERHxgAKwiIiIB4y1tuD+mDF/ADvy8FdWAg7k4e/zkvpS+ARLP0B9KayCpS/B0g/I+76cb609N7MNBRqA85oxZrm1NtzrduQF9aXwCZZ+gPpSWAVLX4KlH1CwfdEUtIiIiAcUgEVERDwQ6AF4vNcNyEPqS+ETLP0A9aWwCpa+BEs/oAD7EtDXgEVERAJVoI+ARUREAlJABGBjTFtjzCZjzBZjzMBMtpcwxnzk3/6TMaZ2wbfy1Iwx5xlj5hljNhhj1htjHstkn0hjzJ/GmNX+f8940dacMMZsN8as87dzeSbbjTFmtP+4rDXGXOFFO7NjjKl3wv/1amNMvDHm8Qz7FNpjYoyZbIzZb4z5+YTHKhhjZhtjNvu/ls/iuT39+2w2xvQsuFZnLou+jDDGbPS/fj43xpTL4rnZvhYLWhZ9GWKMiT3hddQ+i+dme74rSFn046MT+rDdGLM6i+cWtmOS6fnX0/eLtbZQ/wNCgd+AOkAYsAa4JMM+fYGx/u9vBz7yut1Z9KUacIX/+7OBXzPpSyTwpddtzWF/tgOVstneHvgGMMDVwE9et/kU/QkF9uLu2wuIYwK0AK4Afj7hsVeAgf7vBwIvZ/K8CsBW/9fy/u/LF8K+tAaK+b9/ObO++Ldl+1osJH0ZAvz7FM875fnO635k2P4a8EyAHJNMz79evl8CYQTcDNhird1qrU0GPgQ6ZdinE/C2//tPgVbGGFOAbcwRa+3v1tqV/u8PAxuAGt62Kl91At6xzhKgnDGmmteNykYr4DdrbV4mi8lX1tofgLgMD5/4fngb6JzJU9sAs621cdbaQ8BsoG2+NTQHMuuLtXaWtTbV/+MSoGaBNywXsjguOZGT812Bya4f/nPs/wEfFGijcimb869n75dACMA1gF0n/LybvwetY/v436x/AhULpHW55J8mvxz4KZPNEcaYNcaYb4wxDQu0YafHArOMMSuMMX0y2Z6TY1eY3E7WJ5NAOSYAVay1v4M76QCVM9kn0I4NwD24GZXMnOq1WFg87J9On5zFVGcgHZfrgH3W2s1ZbC+0xyTD+dez90sgBODMRrIZl27nZJ9CwxhTBvgMeNxaG59h80rcFOhlwBhgekG37zQ0t9ZeAbQDHjLGtMiwPWCOizEmDOgIfJLJ5kA6JjkVMMcGwBjzFJAKvJfFLqd6LRYGbwIXAk2A33HTtxkF0nG5g+xHv4XymJzi/Jvl0zJ57IyPSyAE4N3AeSf8XBPYk9U+xphiwDnkbvon3xljiuMO/nvW2mkZt1tr4621f/m//xooboypVMDNzBFr7R7/1/3A57jpsxPl5NgVFu2AldbafRk3BNIx8duXPtXv/7o/k30C5tj4F7x0ALpb/wW5jHLwWvSctXaftTbNWusDJpB5GwPiuPjPs7cAH2W1T2E8Jlmcfz17vwRCAF4G1DXGXOAfpdwOzMywz0wgfVXarcDcrN6oXvJfM5kEbLDW/jeLfaqmX782xjTDHaODBdfKnDHGnGWMOTv9e9ximZ8z7DYTuNs4VwN/pk/1FEJZfpoPlGNyghPfDz2BGZns8x3Q2hhT3j8V2tr/WKFijGkLPAl0tNYmZLFPTl6Lnsuw/qELmbcxJ+e7wuAfwEZr7e7MNhbGY5LN+de794vXK9Ny8g+3mvZX3OrAp/yPPYd7UwKUxE0dbgGWAnW8bnMW/bgWN22xFljt/9ceeAB4wL/Pw8B63OrHJcA1Xrc7i77U8bdxjb+96cflxL4YIMp/3NYB4V63O4u+lMYF1HNOeCwgjgnuQ8PvQAruU/q9uPUPc4DN/q8V/PuGAxNPeO49/vfMFqB3Ie3LFty1t/T3S/rdDtWBr7N7LRbCvrzrfx+sxZ30q2Xsi//nv53vClM//I+/lf7+OGHfwn5Msjr/evZ+USYsERERDwTCFLSIiEjQUQAWERHxgAKwiIiIBxSARUREPKAALCIi4gEFYBEREQ8oAIuIiHhAAVhERMQD/w88ixdYpeQgeQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "prstd, iv_l, iv_u = wls_prediction_std(res2)\n", "\n", "fig, ax = plt.subplots(figsize=(8,6))\n", "ax.plot(x1, y2, 'o', label=\"data\")\n", "ax.plot(x1, y_true2, 'b-', label=\"True\")\n", "ax.plot(x1, res2.fittedvalues, 'r-', label=\"OLS\")\n", "ax.plot(x1, iv_u, 'r--')\n", "ax.plot(x1, iv_l, 'r--')\n", "ax.plot(x1, resrlm2.fittedvalues, 'g.-', label=\"RLM\")\n", "legend = ax.legend(loc=\"best\")" ] } ], "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.7.6" } }, "nbformat": 4, "nbformat_minor": 1 }