{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [ "remove_input" ] }, "outputs": [], "source": [ "path_data = '../../../../data/'\n", "\n", "import numpy as np\n", "import pandas as pd\n", "from scipy import stats\n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "plt.style.use('fivethirtyeight')\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# HIDDEN \n", "galton = pd.read_csv(path_data + 'galton.csv')\n", "\n", "heights = galton[['midparentHeight', 'childHeight']]\n", "\n", "heights = heights.rename(columns={'midparentHeight':'MidParent', 'childHeight':'Child'})\n", "\n", "dugong = pd.read_csv('http://www.statsci.org/data/oz/dugongs.txt')\n", "\n", "dugong[['Age','Length']] = dugong['Age\\tLength'].str.split('\\t',expand=True)\n", "\n", "dugong = dugong.drop(columns=['Age\\tLength'])\n", "\n", "length = dugong.pop('Length')\n", "\n", "dugong.insert(0, 'Length', length)\n", "\n", "dugong['Length'] = dugong['Length'].astype(float)\n", "\n", "dugong['Age'] = dugong['Age'].astype(float)\n", "\n", "hybrid = pd.read_csv(path_data + 'hybrid.csv')" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [ "remove_input" ] }, "outputs": [], "source": [ "def standard_units(x):\n", " return (x - np.mean(x))/np.std(x)\n", "\n", "def correlation(table, x, y):\n", " x_in_standard_units = standard_units(table[x])\n", " y_in_standard_units = standard_units(table[y])\n", " return np.mean(x_in_standard_units * y_in_standard_units)\n", "\n", "def slope(table, x, y):\n", " r = correlation(table, x, y)\n", " return r * np.std(table[y])/np.std(table[x])\n", "\n", "def intercept(table, x, y):\n", " a = slope(table, x, y)\n", " return np.mean(table[y]) - a * np.mean(table[x])\n", "\n", "def fit(table, x, y):\n", " a = slope(table, x, y)\n", " b = intercept(table, x, y)\n", " return a * table[x] + b\n", "\n", "def residual(table, x, y):\n", " return table[y] - fit(table, x, y)\n", "\n", "def scatter_fit(table, x, y):\n", " #fig, ax = plt.subplots(figsize=(7,6))\n", " plt.scatter(table[x], \n", " table[y], \n", " color='darkblue',\n", " s=20)\n", " \n", " plt.plot(table[x], fit(table, x, y), lw=2, color='gold')\n", " plt.xlabel(x)\n", " plt.ylabel(y)\n", " \n", "def residual_plot(table, x, y):\n", " fig, ax = plt.subplots(figsize=(7,6))\n", " x_array = table[x]\n", " t = pd.DataFrame({x:x_array, 'residuals':residual(table, x, y)})\n", " ax.scatter(t[x], t['residuals'], color='r')\n", " xlims = np.array([min(x_array), max(x_array)])\n", " ax.plot(xlims, np.array([0, 0]), color='darkblue', lw=4)\n", " x_label = x\n", " y_label = y\n", " y_vals = ax.get_yticks()\n", " plt.ylabel(y_label)\n", " plt.xlabel(x_label)\n", " plt.title('Residual Plot')\n", " plt.show()\n", "\n", "def regression_diagnostic_plots(table, x, y):\n", " scatter_fit(table, x, y)\n", " residual_plot(table, x, y) " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "tags": [ "remove_input" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
MidParentChildFitted ValueResidual
075.4373.270.7123732.487627
175.4369.270.712373-1.512373
275.4369.070.712373-1.712373
375.4369.070.712373-1.712373
473.6673.569.5842443.915756
573.6672.569.5842442.915756
673.6665.569.584244-4.084244
773.6665.569.584244-4.084244
872.0671.068.5644672.435533
972.0668.068.564467-0.564467
\n", "
" ], "text/plain": [ " MidParent Child Fitted Value Residual\n", "0 75.43 73.2 70.712373 2.487627\n", "1 75.43 69.2 70.712373 -1.512373\n", "2 75.43 69.0 70.712373 -1.712373\n", "3 75.43 69.0 70.712373 -1.712373\n", "4 73.66 73.5 69.584244 3.915756\n", "5 73.66 72.5 69.584244 2.915756\n", "6 73.66 65.5 69.584244 -4.084244\n", "7 73.66 65.5 69.584244 -4.084244\n", "8 72.06 71.0 68.564467 2.435533\n", "9 72.06 68.0 68.564467 -0.564467" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "heights['Fitted Value'] = fit(heights, 'MidParent', 'Child')\n", "\n", "heights['Residual'] = residual(heights, 'MidParent', 'Child')\n", "\n", "heights.head(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Numerical Diagnostics ###\n", "\n", "In addition to visualization, we can use numerical properties of residuals to assess the quality of regression. We will not prove these properties mathematically. Rather, we will observe them by computation and see what they tell us about the regression.\n", "\n", "All of the facts listed below hold for all shapes of scatter plots, whether or not they are linear." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Residual Plots Show No Trend ###\n", "**For every linear regression, whether good or bad, the residual plot shows no trend. Overall, it is flat. In other words, the residuals and the predictor variable are uncorrelated.** \n", "\n", "You can see this in all the residual plots above. We can also calculate the correlation between the predictor variable and the residuals in each case.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.090556599598937e-16" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "correlation(heights, 'MidParent', 'Residual')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That doesn't look like zero, but it is a tiny number that is 0 apart from rounding error due to computation. Here it is again, correct to 10 decimal places. The minus sign is because of the rounding that above." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "round(correlation(heights, 'MidParent', 'Residual'), 10)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dugong['Fitted Value'] = fit(dugong, 'Length', 'Age')\n", "\n", "dugong['Residual'] = residual(dugong, 'Length', 'Age')\n", "\n", "round(correlation(dugong, 'Length', 'Residual'), 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Average of Residuals ###\n", "**No matter what the shape of the scatter diagram, the average of the residuals is 0.**\n", "\n", "This is analogous to the fact that if you take any list of numbers and calculate the list of deviations from average, the average of the deviations is 0.\n", "\n", "In all the residual plots above, you have seen the horizontal line at 0 going through the center of the plot. That is a visualization of this fact.\n", "\n", "As a numerical example, here is the average of the residuals in the regression of children's heights based on parents' heights in Galton's dataset." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "round(np.mean(heights['Residual']), 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The same is true of the average of the residuals in the regression of the age of dugongs on their length. The mean of the residuals is 0, apart from rounding error. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.0" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "round(np.mean(dugong['Residual']), 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### SD of the Residuals ###\n", "**No matter what the shape of the scatter plot, the SD of the residuals is a fraction of the SD of the response variable. The fraction is $\\sqrt{1-r^2}$.**\n", "\n", "$$\n", "\\mbox{SD of residuals} ~=~ \\sqrt{1 - r^2} \\cdot \\mbox{SD of }y\n", "$$\n", "\n", "We will soon see how this measures the accuracy of the regression estimate. But first, let's confirm it by example. \n", "\n", "In the case of children's heights and midparent heights, the SD of the residuals is about 3.39 inches." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.388079916395343" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.std(heights['Residual'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's the same as $\\sqrt{1-r^2}$ times the SD of response variable: " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.388079916395342" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r = correlation(heights, 'MidParent', 'Child')\n", "np.sqrt(1 - r**2) * np.std(heights['Child'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The same is true for the regression of mileage on acceleration of hybrid cars. The correlation $r$ is negative (about -0.5), but $r^2$ is positive and therefore $\\sqrt{1-r^2}$ is a fraction." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.5060703843771185" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r = correlation(hybrid, 'acceleration', 'mpg')\n", "r" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(9.43273683343029, 9.432736833430296)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hybrid['fitted mpg'] = fit(hybrid, 'acceleration', 'mpg')\n", "\n", "hybrid['residual'] = residual(hybrid, 'acceleration', 'mpg')\n", "\n", "np.std(hybrid['residual']), np.sqrt(1 - r**2)*np.std(hybrid['mpg'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let us see how the SD of the residuals is a measure of how good the regression is. Remember that the average of the residuals is 0. Therefore the smaller the SD of the residuals is, the closer the residuals are to 0. In other words, if the SD of the residuals is small, the overall size of the errors in regression is small.\n", "\n", "The extreme cases are when $r=1$ or $r=-1$. In both cases, $\\sqrt{1-r^2} = 0$. Therefore the residuals have an average of 0 and an SD of 0 as well, and therefore the residuals are all equal to 0. The regression line does a perfect job of estimation. As we saw earlier in this chapter, if $r = \\pm 1$, the scatter plot is a perfect straight line and is the same as the regression line, so indeed there is no error in the regression estimate. \n", "\n", "But usually $r$ is not at the extremes. If $r$ is neither $\\pm 1$ nor 0, then $\\sqrt{1-r^2}$ is a proper fraction, and the rough overall size of the error of the regression estimate is somewhere between 0 and the SD of $y$.\n", "\n", "The worst case is when $r = 0$. Then $\\sqrt{1-r^2} =1$, and the SD of the residuals is equal to the SD of $y$. This is consistent with the observation that if $r=0$ then the regression line is a flat line at the average of $y$. In this situation, the root mean square error of regression is the root mean squared deviation from the average of $y$, which is the SD of $y$. In practical terms, if $r = 0$ then there is no linear association between the two variables, so there is no benefit in using linear regression." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Another Way to Interpret $r$ ###\n", "We can rewrite the result above to say that no matter what the shape of the scatter plot,\n", "\n", "$$\n", "\\frac{\\mbox{SD of residuals}}{\\mbox{SD of }y} ~=~ \\sqrt{1-r^2}\n", "$$\n", "\n", "A complementary result is that no matter what the shape of the scatter plot, the SD of the fitted values is a fraction of the SD of the observed values of $y$. The fraction is $\\vert r \\vert$.\n", "\n", "$$\n", "\\frac{\\mbox{SD of fitted values}}{\\mbox{SD of }y} ~=~ \\vert r \\vert\n", "$$\n", "\n", "To see where the fraction comes in, notice that the fitted values are all on the regression line whereas the observed values of $y$ are the heights of all the points in the scatter plot and are more variable." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAEjCAYAAAC8fj/5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABfrUlEQVR4nO2deXhU1fnHvzOTyWSZhJCFSUIgbIGwJIGAgEYsqAWpiAIqWAtWS6JYFVqxrVShLAoqtmjLr0BwA1Eji6JVRFtBAZFFlIBsgZCQkISEJQlD1pm5vz+Gmcy9c9fZJ3k/z8Oj986557zn3Hvz3nPe5ahqa2sZEARBEEQQova3AARBEAThKqTECIIgiKCFlBhBEAQRtJASIwiCIIIWUmIEQRBE0EJKjCAIgghaSIkRBEEQQYvflJjZbMaSJUuQmZkJg8GAzMxMLFmyBCaTyV6GYRgsXboU6enpSExMxJ133onjx4/7S2SCIAgiwPCbEluxYgXWrl2Ll156Cfv378eyZcuQn5+Pv//97/Yyr732GlauXImXXnoJX3/9NRISEjBp0iRcvXrVX2ITBEEQAYTKXxk7pk6dis6dO2PVqlX2c4899hiuXLmCgoICMAyD9PR05ObmYu7cuQCAxsZGpKWlYfHixXj44Yf9ITZBEAQRQPhtJjZy5Ejs3r0bp06dAgCcOHECu3btwi9/+UsAQGlpKS5cuIBbb73Vfk14eDhuuukm7Nu3zy8yEwRBEIFFiL8anjNnDoxGI0aMGAGNRgOTyYS5c+di5syZAIALFy4AABISEljXJSQkoLKy0ufyEgRBEIGH32ZiW7ZswQcffIC1a9fim2++wapVq7B27VqsW7eOVU6lUrGOGYZxOidEUVGRx+QNJqjfHQvqd8eio/ZbCL/NxObPn48nnngCU6ZMAQAMHDgQZWVl+Mc//oEZM2bAYDAAAKqrq5GSkmK/7uLFi06zM4IgCKJj4reZWENDAzQaDeucRqOBxWIBAKSmpsJgMGDHjh3235uamrB3716MGDHCp7ISBEEQgYnfZmJ33HEHVqxYgdTUVKSnp6OwsBArV67EtGnTAFiXEWfNmoVXX30VaWlp6NOnD5YvX47IyEjce++9/hKbIAiCCCD8psRefvllvPDCC3j66adx8eJFGAwGPPTQQ/jTn/5kLzN79mw0NjbimWeeQW1tLYYOHYotW7YgKirKX2ITBEEQAYTf4sR8QVFREdLS0vwths+hfncsqN8di47abyEodyJBEAQRtJASI4h2QElJHcaO3Yrs7AKMHbsVpaX1/haJIHwCKTGCaAfk5e3E/v3VKC6ux/791cjN3SF9EUG0A0iJEUQ74OLFJtFjgmivkBIjiHZAfHyY6DFBtFdIiRFEOyA/fwyGD++CXr2iMXx4F+Tnj/G3SAThE/wWJ0YQhOdITY3Gl1/e7W8xCMLn0EyMIAiCCFpIiREEQRBBCykxgiAIImghJUYQBHEdChoPPkiJEQRBXIeCxoMPUmIEQRDXoaDx4IOUGEEQxHUoaDz4ICVGEARxHQoaDz4o2JkgCOI6FDQefNBMjCAIgghaSIkRBEEQQQspMYIgCCJoISVGEARBBC2kxAiCIIighZQYQRAEEbSQEiMIgiCCFlJiBEEQRNBCSowgCMLPUPZ81yElRhAE4Wcoe77rkBIjCILwM5Q933VIiREEQfgZyp7vOqTECIJwG7LpuAdlz3cdymJPEITb2Gw6AFBcXI/c3B2UDV4BlD3fdWgmRhCE25BNh/AXpMQIgnAbsukQ/oKUGEEQbkM2HcJfkE2MIAi3IZsO4S9oJkYQBEEELX5TYhkZGYiJiXH6d//99wMAZs2a5fTb7bff7i9xCYIgiADEb8uJO3bsgNlsth9XVVVh9OjRuOeee+znRo8ejdWrV9uPQ0NDfSkiQRAEEeD4TYnFx8ezjtevX4+oqCiWEtPpdDAYDD6WjCAIgggWAsImxjAM1q9fj6lTpyIiIsJ+fu/evejTpw+GDh2Kp556CjU1NX6UkiACF1cyZvgry4ZYu5T5g1BKQCixHTt2oLS0FNOnT7efu/3227Fq1Sps3boVS5YswQ8//ICJEyeiubnZj5ISRGDiShZ0f2VOF2uXsrkTSgkIF/t33nkH2dnZyMzMtJ+bMmWK/f8HDhyIwYMHIyMjA9u3b8fEiRMF6yoqKhI97ihQvzsWlZVXnY6lxsKVazyBWLtKZeqo97sj9TstLU30d78rsZqaGnz++edYvny5aLmkpCQkJyejuLhYtJxjh4uKiiQHoD1C/e5YFBUVISkpCmVljfZzSUlRkmORlHRM8TWeQKxdJTJ15PvdEfsthN+XE9977z3odDpMnjxZtNylS5dQWVlJjh4EwYMrGTP8lWVDrF3K/EEoxa8zMYZhsG7dOkyePBlRUVH280ajEcuWLcPEiRNhMBhw7tw5LFq0CAkJCZgwYYIfJSaIwMSVjBn+yrIh1i5l/iCU4lcltmvXLpw5cwZr1qxhnddoNDh27Bg++OAD1NXVwWAwYNSoUXjrrbdYyo4gCILo2PhVid1yyy2ora11Oh8eHo4tW7b4XiCCIBRTUlKHvLyduHixCfHxYcjPH4PU1Gh/i0V0EPxuEyMIIrght3jCn5ASIwg/ITew1xsBwCUldbjlls0wGN6EwfAGfvGLLS7XK7YhZrAFLwdiADghDikxgvATcmcw3pjp5OXtRGHhZTQ3m9HcbMHhw5dcrldsQ8xgm6UFYgA4IQ4pMYLwE2IzGFfKudO2O/WKucV7Q3Zv4i95g22cAglSYgThJ8RmMK6Uc6dtd+q1ucUfOjQVX355N8upwxuyexN/yRts4xRIkBIjCPjHJiE3sNcbAcD5+WOQmRkLnU4DnU6NrKw4rwQWB1vwsq/ltT13FRXXoNeHoHt3fVCMUyChqq2tZfwthLfoqOlZqN/KGTt2K/bvr7YfDx/eJWiCbul+By+uPHftod+ehGZiBAGySRD+gZ479yElRhAgm0Qg0NHczEtK6lBd3cA6156eO1/dT1JiBIHgs920Rzqam3le3k4YjSb7sV4f0q6eO1/dT79vxUIQ7uCY8igyEnj3XQPLO05OSiRPpk3ydQqmkpI6PPLIIVRXf48LF9qWosLCNNi8+Q7k5CSzZNJqgbKya2hpsSAsTIOCgnHIyUmW3ZZQ34R+UzIeFRXXRI/dkU/u9TNm/BcnT9YBYJCe3hnr1t3utfvH7Z/RaMLQoR8qvi+Biq+WSsmxox3SkfotZRiXYzj3pFOHrx1EuO05oteHoLz8YVllXGnLsW9CvykZj5SUt5xmJmKycZ9zd8eeb5y8ef+4/XVErO/B8n776l2g5UQiqJH62pPzNejJL0ZfG+rF6m9qMssu40pbjsdCvykZj9jYMNFjd+Rz5XpX6lCCWP+U3JdAxVdL9KTEiKBGyiFDjsOGJ506uNeGhqqRkvIW4uPXIiXlLezZU+Fy3XLacyQsTCO7jCttOR4L/aZkbBMTI0SP3ZHPletdqUMJMTGhgr8puS+BilgQvCchJUYENY5fexkZ0U5fe3K+Bj35xcitq7S0HkajCSYTA6PRhKlTt7tct1B7mZnRMBjYf2xtdhWuTP36RSMiQoOQEBX0+hB7GVf6JmdHZiVj6+598MT1vggAt6FSqTjHcOm+dHTIJtYOoX4HDnFxa2E2t71iGo0Kly7N9GgbRUVF0Gq7dLg9vQLxfishO7sAxcVtbuc6nRr7998ned+Cvd+ehmZiBOFFGIYRPfYUHc09vT3AXapsbrbQfXMBUmIE4UUSEsJEjz0FZX4IPvLzx0CnY9u+6L4ph5QYQXgR7tKQt5b4KONI8JGaGo2srDjWObpvyqFg5w6IrwNyvYlUsPPu3eWYNu0rNDWZ7c4OXbtGsvo/f/4wLFp0kPcYYFBcfBUAoFZbv56nTOnj1HZ8fBjuuacnnntuHywWa9u9e0ejtdUCvT4EMTE6xMbq0NxsRnZ2AeLjw5CX1x9z5uxhycYNTua7P9zf583rifnzh2HKlM/R3GwtU1h4CXv2VNgDZvnqZBgGeXk7UVXVgMuXm67LGAbA6oTi2Pbu3eWYMmWbvX7HYGqh+uU+U4H0PPpalvz8McjN3cFqj1AGOXa0Q6T6HcwZ27lI9YUvgHbAgFjWNXp9iFMZoSBUtRq4fDmXt20xhg/vAgCs8mo17ArP1i5fcLJUAHdmZjTCwsKcZHEMmOWrkyuPkNxffnk3b2CuVP1ynylXr/XG+x0M70ZH/bsmBM3EOiDtyX4i1Rdu0GhTk1lWGSEclY6SceMr61iXY7tKA7avXGmFVuvcpmM/XL3ntnJ8Y+KJ+t291tMEkiyEPMgm1gFpT/YTqb5wg0bDwjSyygihdnhjlIxbfHwY9Hq2puGECQkGJ0sdx8RoeWUJDW0Tlq8OOfLbyvCNieM5d56pQHoeA0kWQh6kxDog7Slju1Swc0HBOOj1IawgUm7/CwrGCR737t1mD7HZxPjaHj68C5YtG8lScjbaspOzV+579Yp2ko2vXqkA7sWL+1/3dGNrxe7d9YLX5OePsZ/r3l0PvT4EKSmRyMqKQ2ZmrFPbBQXjoNO11e0YTC1HZjEC6XkMJFkIeZBNrB3iTr+tmbz/h5MnrwBQoV+/Tli//peS2coDgR07CrF06Vm/y+YcxKpB166RKC83oqWlbQ0xIUGHxkazk2OHUmz3e9Cg91Be3pYZPSUlEkeP/tq9zriJN58Xer8JgGZiBIe8vJ0oLLyE5mYLmpvNKCy8zArADOSg2vnzjweEbM5BrGYUF9ezFBgA1NQ0ezQlVW1ts+ixPwjk54VoH5ASI1hIZfIOZMN3bW0r69hfsjkuSel08l8xdzOXu5sF3hsE8vNCtA9IiREspDJ5B5Lhm7v9eUQE2/lAiWy7d5d7LNu8Y/burKx4wXJc+5m7mcu5WdHFsqT7ikB6Xoj2CbnYEyzy88dgxoz/4sSJNpuYzbi9e3c5jh69BMDqWde3r7MjhTfh2leamkwoLLwMAHYblFarhloNpKd3tsvGvY4b3JyfPwb33/8lGhqsMyGj0YT77vsCFRWPiLbPV48tgNh2bsGCGzBv3vc4ceIKGAbQaICYGB2uXm1FZKQWFy402uu37uxbgC1bxtt3RZa703BJSd11O2YbLS3smd3u3eWYPPkLtLRYzeAqFbB27RgMHZpgl9nqQekc7Cxl2+ILKs/JSeYN5vWUnay8vBG///1Wp/tRUXENtbXNiI0NQ2JiBPLzx6CsrJ5XPqF7G0i2XkIccuxoh3ir30p33vU03EBUnU6D5mb+JTixXYe5wczDh3fhDfqtrc0VbZ+vHgCSQcViwdSOssvdabioqAi///0xp7I6nQYXLrQpYr6AZbUaGDaMv/9issgJKhd6NjwVUDxqVAGOHGlzoBEa1+HDu+DYscui8gVDkLONjvp3TQiaiRGyURIU7A2c7SnC319idjyu3BcvNkGlAhwTzHNjuOTWIy2z9LgJ7YosdE74PHt8+Nq1WMTtVHJ3aFbybHjKTsa1gQq1efFik6R8ZLsLXsgmRshGSVCwN+DaU8R2NRGz43Hljo8PQ9++7KUj7rGceqqrGxAaytZ+en2I03WO+4vxERWl5W1P6JzQ+fT0zqLyAtaZmJidSu4OzUqeDU/ZyWJi2MHjQm3Gx4dJyke2u+CFlBghG77AYV9i8/qzbV9hc1lXq63/VCogOTlCctdhbnBzfv4YfPjheNa5Dz8cL9i+Yz16fdtihtFowrlzRtY1KpXKSW6b8tXp1Kxgahu2PceU7DTMV3bduttZZawBy21KVqUCK+i5V69owWBnqSBgJc+GpwKKFy/uz3tfU1IiodeHoHt3vb1+KfkoyNkzqMw1CGn+H0KNKxB+5RHoq4dBZa70apu0nEiwEDNw5+Qks+wINu9AXxnDbV5/3GDiHj2isWXLHZg+/QtcuyZ8nSN89g4xGwjXceGf/xyFigqjkw2GGwt2+PAl5ORsRkHBODz66E5WMHJCQjh++GGqU39sdaamRuPbb6cIygQAGzeeQl7eN3bFaDCEIzU1ivdeMIwFKpUagHUpLTQUSEyMsI9PW6B7LQAVS0HzjaGNkpI6LFx4AF26RNgdLBYuPCD4XHDrcvU54s7EU1L0gjIyDIMBA2LtbaSk6Fm/S/XP204fQedYwligNpdA3XoEmtZCaEzW/6otVU5FNa2FMGmSvCYKOXa0Q9zptxIDt7+M4XIysntaFj7HhYYGk1MSXyHnAptC4HMucGccY2PznWQQqkNpJnq5sshxeBGrw9X+cx07vPWs+uI5V9KGz/+uMS1Qm45Dw1JYR6FirjoXVelhDsmAWWv7lwlLSDqg0vFU7Bn8NhPLyMhAWVmZ0/mxY8fiww8/BMMwWLZsGd555x3U1tZi6NChWL58Ofr37+8HaTsOSgzc/jKG87ltT5q0zauy8DkG8CmPgoJxWLjwgJMyaGoyIzk5EkZj23KjLRjZnT2l+GQA5DuUiGWiFzonVUaOw4vY73LvnZLg9kDPsh8wjiWWOmhaj9pnVprWQqhNJ6FCq3NRdaJVUYVkwqzNgEWbCYumB6DyrZXKb0psx44dMJvbHvaqqiqMHj0a99xzDwDgtddew8qVK7Fy5UqkpaXh5ZdfxqRJk3DgwAFERUX5Ser2T3x8GGtpS8roL7esJ+AuuXz00Xj7kou3ZQkL07BmF2FhGqeZmFptXXLl238rLEyDxMQIls0sMTECgPhSlhTcPclscPtfUlIHi8V50YWbid5xDPnq4YN7HXespOrgbiPDt60MHzExWpSVtcXZeetZ9cVz7ut3CQwDlaXKrqg0piNQtxZCYy5xLgoVzJo+1lnV9dmVOSQDjKaLd2WUScAsJy5fvhyvv/46Tpw4gfDwcKSnpyM3Nxdz584FADQ2NiItLQ2LFy/Gww/Li02i5UTllJbWO80KhNbmlZT1BGJLLqWl9Zg+fRuMRnhFlj17KjB16nZWsGxVVQNyc3fAYnHe9ZmvfEqK3uPjtXnzabsMgLBNjG+pUKcDtmy50x70W1pa7xTo7pj8WQjuc7BgwQ2iNjEuyclv2gPNASAiQuMUaM7Hzp2FePFFeQmf3XlWffGcK2lD8fvNmKE2n7EvB6pbj0BjKoTactG5KEJh1g6A5frsyqqwBgJqPU/FgUFAKDGGYTB48GCMHTsWr7zyCkpKSjB48GB8/fXXyM7Otpe7//77ERsbi1WrVsmql5RYG3IMx+5uMS+U/d7VernXVVU1sGYyvXpF49ChqQCsDg65ud/Yf+vdOxpbtowHwzCYOvULnDpVD4ax/oF8/fVbsGbNz/Z68/L648knd6Gx0QKVCujXLwYFBeN4Zc/NHYA//GE3mprMCAlRg2HMaL2+0pKYaHVsUKlUuHq11V73nDl77MpsxYqbsWbNMXt2jMuXr6G8nL10FB2thcXCIDxcg0uXrEl8IyJC7FkmHGXS67VoampCXZ0ZFy6w60lPb+sH13mke3c9EhMjrtcRwpLZlT/SNseXhgaTfSxSUvSynrODB6tZs8mQEBUuXpwp2abQ+y2UPSTQcPW9EP27xjRB03oMaoflQI3pZ6iYBueiqk4Oispmv+oLqGROhQOEgFBiX3/9NSZPnoxvv/0WmZmZ2LdvH8aNG4cjR46gW7du9nK///3vUVlZiS1btsiql5RYG3IMx540fjte72q9ShwG+Bwc+Bw+AOclOL4lOSHZhZbvhJDTllyEHEHEEOqHWNYQVxwX+JxGhOqSGlO5mWCE3m9/Z5aRi6vvha3fKssV6xKg3dniCNSmU1CBJ6hdnQKzdpBVYdmXA7vzR/UHGQHhYv/OO+8gOzsbmZmZrPMqzgAzDON0jktRUZHocUeB2+/KyqtOx66UEYJ7reP1rtbLvU6v16BXrwhcudKKmBgt5s3raa+HTzHwycRXlu/aH36oxqhRBaipaZYsK4actuTS2GjiHU8xbGM9b15PzJ/fZB872zYwYtcolU1uXVz54+NDYTSa0dxsRmioBsuXD5TdPl85riy2cQs05L8XDEJVVYhQn0KE5iR6h51C2PmT0Kmd3dkZRo1GS080WPqhwdIXDZZ+aDT3hQkxnJItAE57qiteRWoi4nclVlNTg88//xzLly+3nzMYDACA6upqpKSk2M9fvHgRCQkJovU5dphmYm0kJR1jGcGTkqJcKiME91rH612tl3td9+6dBL9U1eqdTgoiKcnqAMSVS87syGwGjhypZ8VKCZUVw5MzsfDwEN7xFMM21mlpwLfftn0kjh27FVVV/LM5Jfe9TbbdvEpRznPWo0eMS44tQu83VxbbuAUavO9Fn55Qm05Z7VemQuvsqvUI1MwVp+sZhMOsHXh9GfD6sqB2AKAKRwiA6Ov/2jt+z9jx3nvvQafTYfLkyfZzqampMBgM2LGjbQO9pqYm7N27FyNGjPCHmEGPnIwE7m4xn5UVB51ODZ1Og8zMWMyfPwxjx25FaWk91GpAo7EG0C5YcIOsOufPH8bKsiB2HVfW3r2j7dko+vWLtq+aRERoePupEwhjiYjQ2rdMUauBp58ebJeJL82RVgtkZcWx6nbsAzc7RkoKvxeaSmXd+dmx7ejoUIwduxULFtzAqqNfv0gYDHxpp2JEM3zY6sjMjHWSWQklJXXo2jXCPsZqNX/mFL62vZEdw9+ZZeSyNn8EHnmgHs898SM2/vNTvDF/AcLOJSLq4k2IqHsUumsrEdLyLdTMFVhUsWgNHY3myKdQ3LQEV+P3oT6xHNfi/4umTn9HS+TDMIcOBVTh/u6Wz/GrTYxhGAwbNgw5OTl4/fXXWb+tWLECr776KlauXIk+ffpg+fLl+O677xS52NNMzL8I2W5ctYlJXefJIG8bUvYaTwTCSgUZS7Xh7/vtr6B3f/dbCSrzRatHYGthm5eg+TRUPEmsi8tiUVrVCzeOGXt9dpUBRp1st18FU799gV+XE3ft2oUzZ85gzZo1Tr/Nnj0bjY2NeOaZZ+zBzlu2bKEYsSBCWcZ16XLeDADl1h0SokJ2dgIOHnQOWnbEnUBlobZ1Og2rnoAJhBUg0OXzKQxzPR1Tm7OFNR2Tc/5ABiEwh/SDWZuJl14z4b+743D4RDLqroZbPW8nTvVDB4IPUSXWuXNnSUcKPi5fviyr3C233ILa2lre31QqFZ599lk8++yzitsnAgO+4FnbeTFsrsfnz7OT6er1IV7L1ciVNTs7QTBo2RFGLJU+ByGXam7bWVlxrH75PBBWALnyK5Uv6PIG2mBaoDadcEjHVAhN689QMc7PvDUd0yBOOqb+9nRM/9mzFfsPtH0wtYcs+r66r6LLiUuXLnVSYp999hmOHz+O2267DX369AHDMDh9+jS+/vprDBgwAL/61a/wl7/8xeOCukJHnXYHSr9tAZx8O+2KPcx8m19mZcWhudmMw4cv2c97cllNKNiUL2jZMebIE7kmpQJdpX731f12VX5X65XCp8+5pR4a09HrsVfXZ1emEwLpmAwO6ZisWS4smp6i6Zi8GuzsJ3y1zCw6E+POgtavX4/q6mrs3bvXaRBPnjyJu+66i+VNSPgXT38JcQOae/SIRFiYVjBI1jE7uk0OObKeO8f+km1psS7hXbnCdne/eLGJdV1kJPDuuwaWDFaZ/4uTJ+sAMEhP74wXXhiBRYsOsra1f/rp3faA6KqqaygvNyI1NZqVuZ/b/65dI1BczHaT/vHHahgMb/K2xZ1ZVlU1YOzYrSwlz4fYbG/37nLcd98uNDbudPotLEyDzZvvcAqQduVZKCmpw+HD7AwPtvuZmhqNNWtGIy9vJyoqriEnZ7PsDxbHemwcOlSDsWO3uvy8uhXsbE/HdISTjuksb3GzpndbOqbrWS4YjUFWU0E7A5VJaWm96LGnUOTYccMNN2Dq1Kn2VFBcXn75ZWzcuBEHDhzwmIDuECxfLJ7G1m9PfwlJBdkK1e9KoDUffMHOgHj2er56ufXwBf7yBcgqCTIWq1vqN6k+OP4uFGTM7Ye7z4JS5xMl7bh6rdvBzozFno5JzdpOpMa5KEJhCelvd7RoS8fkuo0+KGagbsBNQKBWA5cv53q8HUWOHWVlZdAJ+SIDCA8P581MT/gHTxvclWYjVyKHHNliYnSsPaHkZK+Xk81dKru7Ehml6tHpNOjaNRLx8WGoqLjGq4Ck+uB4zCcnX/vuPgtKnU+UtGNzjjl0qAYmU9s3tavPK+/9ZZqgaT3Ok47JeQM6RhXtoKiuZ2gP6efxdEzt3SHGaopiOMeeR5ES6927N959913MmDEDnTp1Yv1WW1uL9evXo0+fPh4VkHAdTzsECDlqSNUvRw6pugEgOTnS6UtVqm6+ernZ3LlZ123nlMrIF8jMrVurbXuRY2N1rE0y5fZBr9faHVz4MtM7YmvP3WdBqfOJknZsy6UajRomU5sCcvV5TYhtRL8e5RjS/zwG969A9oAKRFf9WSAdU1eWs4U1HVOqT9IxBYrDjrcID9dwgs6d3ylPoGg5cfv27XjwwQfRuXNnPPDAA3aFVVRUhA8++ABXrlzBhg0bMG5cYAQXBsu0WylSa+m2fns6+zY3y3mPHpFQqVQ4e9YImw1o3brbndrYvfs8pk37UtRG4SirVguUlV1Dc7MFDMOgc+dQNDdbEBVltb852lsA2K+z2sTGOzlFTJ/+ld0mBjBobm575MPD1di0aTyeeGInzp61KhRuRnqh/nftGoHKymtobWUQFqbBa6+NwmuvHWbZ3558MgNz5lgTBjvOMgAgPb0ToqN1KC2tR01NExjG+reTmzyXex+bmkwoLGzzABbLAtKtWySOHPm1/R44JuiNjw8HwMBoNIk+HyUldXjoof+JZrfnq18oATAXZ0ceNbKy4iWv3bHjMN5c/QMSO5/GsIwKTL27EVHaY1CbnVeDGKhhCUmzz6xsdixGHcfqpy9tVHLeCz6C5e+alFOUp1Ac7PzNN99g/vz5KCwsZJ3PzMzEwoULMXr0aE/K5xbBcrOVEkjBr95OLMx3vVA9cvptMLyJ5ua2L3KdToMLFx7xmieVmB1NqG0pGbgZ6bt1C8eRI7/hrUelAq5cyZVlz/OkTVPJ+HH747g7gR3GBLWpiLVZo/HSD4iJcp7JNjVrcfhEIn483hU/nUhGdX1fvPXeU4AqQlQOXwdtt3ebmK9QHOz8i1/8At988w2qq6tx7tw5MAyD1NRUdOkSGBukdQQCaS3dFXuXu/YYV+uxwv1mY3jr8tSYitfD37bUtdxlqJgYrWB5m2OjKzs0C5339D3m9qdrsgqalgPX3dhtSusYVGDXGRMFXLwSYVVWx5Nx/lIfLH45F+PuOoa937NDMaQUmLt9cIVAeo+DGZczdnTp0oUUl58IpLV0V+xd7tpjXK0HANLTO7NizdLTOwOAU6Jf7rGriNnRbG0r7R83S8i8eT0F64mI0EjKIdWeN++xynwR76/V4LMtB5BqOIuMvuXokXQBqkvO66MWTXfWcuA900rw350aAFb71fDhXWDRpmPV6mSXsqj4+r0KpPc4mBFdTtyzZ49Llebk5LgskCdpr9PuQAl+lSOL3DJy2pAKmpbTbyFZfvGLLSzllpUVh2++mSxSkzLZxTafdDUonNvv0tJ63H//NtYGoBs3WuPE+Nro3FkHhpG2iXnkHjMMVOZSh8wW1/e/slQ4tcdAA8v1dExme3b2DEAdwyqnZGdnOfh6p3JX22uvf9dcRfRzc8KECYrcIm37fclNO9UekGMMlgq+VGpQ5gYRT5q0zWMvnZJAUcdAYovFjKqqa7jrrs+c/gDb5BVrz9EhQKdTobT0GiyWNicLoaBp27kzZ2px6VILgJ2CjhlistTUsLc2OXbsCrKzC6DXayHl/GCToazsKiorrfU47hBta89W7upVa5aHsrKryM3d4SC7FYNBzWpn48ZTePTRb1jj4di3lStP4+23dzr1qanJjKoq646+Ng/A0FANBgyw7jDgGPCdmzsAOTmbee+72P3jHVem1ZqOqcExHdNRgXRMkdfTMTk4XISkAyo5Xo38511x0OAGxUdFaVm/ecPhQ0nKMkIY0ZnY7t27Xar05ptvdlkgT+KLLxY5xlmp4EtP7Xxsu86dfivZFVfMWUBuH6SCdYG2IEm+/gLOOzc7XiMXvp2h+VAaqC0WECx3d2WpoNGYmHxBeYXGjtu2q7srW9Mx/cyTjqnFuai6i312ZbmutCyaXqLpmMQYNaoAR460KUZ3dhL3xs7kUpBjh2cQnYkFijIKZDwRXOuqgdcbhmE5gcBy2pMri1SwLtD2x1VJf23XePorWmmgtuNv3HJyx9ad3aGFxo7bNrdOPtlU5ir7Zo1qW1omczFvu2ZNL046pkzZ6ZjkUlvLzlto66Mr74XYffWV0w85driG33d2DnbkGGe5Aa/cQFpXDbzeMAxLySrWPvc3V9oTg9tedXUDevXiV0i2zSTz8nbav3aLi602CL6v3YgI4VmRI3yZ9OWOA9dZRKtVOcWO8V3Htzu0XGxluTJyx92xDZXKgkF9a6Ft/Ajq1kK0XD2E1vqfEBdT51Q/Ay1POqZBbqVjkktMjJa1M7JtzFx5L/juoTv1yYEcOzyDqBJ76aWXoFKpMHfuXKjVarz00kuSFapUKvzpT3/ymICBjpz9pAoKxjkF/Smtw9W2lSIlK7d9WyCxxWKGVqtBfHw4KxBZbnuONrFLl5rQ3Nz2V7t372h7ezk5m+1/fI1GE1QqFYYP78KyK9nsRoD8r13Hfmu1KvTo0QnNzWZERWlZzg/NzWYnpWi7D3w2Mcdx4NqXe/TohKgorZNNzLYrteM45+buYNnEHHnkkRS8+Wa5U58cy3KflQULbsDChQdQX3sVwwdfRO50E34+tBMZaeXI7FcBfUQLUGutJ0wNIAaorQ/DTyeSUXWlDybef49DOqZQ3jH1NosX93dy7ODrq5xn0fFZtgWqu1OfHLxVb0dD1CZm20+sqqoKoaGh6Ny5s3SFAeTY0VHXjuX2OxCyaPPJMGnSNsHgV7HAWL5+y7E7OMog5EFoQ6x9ofG0Oa9wZ3rdu+uRmBjByqbv6Gwh937s2FGIpUtleOlZalFT9j0+3bQVqV3OIqNfOXqnVEGtcp6BVl/ujLCYbIRGD8Efnr2K7Ts742x5LAAVfzAyz1h6+5mi95sAJGZiV65cET0mghu5S22+lkFsmUXpEoycr11HGRzhGxOx9oXGk0+BAcDly004d85oLz9t2pf2ckrux/z5x+0ODtbrvsZX24a1bSdis2OZz6FTKPCHX7dda7GocLy4C348noyfTiTbA4cvXtHbFf6hU1txtlzeho2B8EwRHQuyiXVgfGlYFvpC55Pho4/GCyoepUswctzDlTioiLUvNJ58ThLDh3dxymIv5fDDC2NGQvQ5/HpCMQanV2DIgPPIHlCJ6Gqe7OwIQ+HJRHz/k8GekunYmSRcNfIvB9rat/XZFmNm2weNb5ZFzgqEryEl1oHxpWFZ6AudTwYxxSNHKSlFiYOKWPtC48l1otDrQ+yu245Z7LnlnO4H0wBN6zGndEz/W8uOcQMAi6qz3Y3d7tYekoa8hz5jzTrFnERs7dv6bJPXaDTi3Dkj7yyLnBUIX6NYia1btw7vvPMOSkpKeJcXVSoVLl26xHNl+8Qfwc6eYv78YazM4+XlRtxyy2ZRm5AUQn2tqGDPDMrKrtp3NdbrQ5wy0/PBHSdHGxLDmK9nod8JwOpYkZTEn0mdK+Nrr41CXV2TPdMFAGi1aqjVQHOzGaWl9WAYRvAe2eRy7Et4uBpHj15CXNxa1nYpYWEadO0agezsAkRFaZGZGYsLFxpQU9OEa9es90GtBrolNeH1Za24VvYSCvfvRO+UEvRKqYJG7WzCLq3ojEPHkvHTcetyYOGpZMxfNAlTpjjbTXJzB+DgwWq7JyLXtV6rVUGtVgFQoanJhNLSens/xWZZjgHZAGAwhCM1NUqRs4LS96C8vBG///1Wt22Kjm27mjWF8B+KstgvWrQIK1aswMCBA5GTk4OYmBjecn/5y188JZ9btPdgZyHk9tudzOZCCPWVe57rNu5KQKpYsLBYH/hkBCAaeAwI7yDNd/+OHbvMW5/z7tQJuFZ7Gmndz2FI/wrrHljpFeiW5OzObjKpUVKZjO59c+w5BPsPPoLySuflQKFgbzk7QXN3zxbrp9yAbDkofQ+4wc5isitt25U6fAU5drBRNBN799138atf/Qrvvvuut+QJOvwZ7Owu7mQ2F0Kor7GxYTAajfbz3F1fXZFFTqC0nOuk6pEKcOa7f3x1hoSYkda9Ehlp1s0ah/S/br/SOy8HGq+F4mRpNxwrTsE3+xLw0/FkHC1KRNeUOJZnYFXNKThn5hcOiBbr6/DhXVBV1cC6T459E7MHuhOQzdcW3zEXbrCzSzZFibJk0wt8FCmxa9eu4fbbb/eWLEGJP4Od3cWdzOZCCPU1MTHC7okHOO/66kpAqtxAaW7dQjIK1WW7XonHZIPxMnp1LcOQ/hV2hTUorQq6UGclcuGSHj8eS8ZPJ7pavQSPJ+P0uXgMG5YIgD0DlOqLDSFbl1B5Wwhb58461n1ybE/MHuhOQLZjW0reA26ws5hNUWqpUulOAkTgoGg5ccqUKejVqxdeeeUVb8rkMXwx7ZaTiVpqh1NPZ8+W229uZvOWFgtaWtr+Eun1IdizZ4oiWYT6yu2jLdhWSZ/F6miziVkRs4nxyQiAN9jZccYhdI/Ol57G6tffR0r8GQzuX4FbhtcgFCVQqZxfrWut3fD9j4n48Xgyzl/sjdwnH8CUqT/ixIm25UOtVo0hQ+Il27X15b77vkBDQ5tyFEuAbOt7Y6MJZp5JWWZmLMLCQhQ/i5s3n3YKyOZrXwyl7wE3i73YMyW1VOnuTgK+hJYT2ShSYsXFxbjnnnvwpz/9CQ8++KCiDPf+wBc3W6kx2hdOHK4GO1dUXGN5y4kFtfoSuWMmp99y6uJmNE9P74x1625Hanc91OazVu/A67FX1u1ELji109KqwamSZJyt7IHsm25HJ8MNOFvRHTPzfnBqW9bOxgJ90Ou1uHTpKioqWsAwVgWSmMivvJUGlnMRc1ASGle54+3KO6Hk/VY6xoEMKTE2okpsxIgRTufq6upQXV2NsLAwJCcnQ6NhL42pVCp8//33npfUBQLFscOd8q7gqmOHO4ZxbyJ3zOT0W05dY8duxU8/VmBgWpV1OTC9ArcMr0Zmv0qoGCO3ShgbwvHDz0nXvQOTr8dfGdDaGsJqQ6htd58hIbj1yNkFQKxtMQcld/rmi2zuvnjvfAUpMTaiNrH4+Hin2VZCQgL69FG2TNCeUWqMDqRgUG7bsbFhGDAggvVFHAh4csx467LUQtN61D6zenvBLvTpVgGtluOdwAAWdVJbstuQDJi1WcgesRdnzjgrN26bQv1QGsDtqV0OpALLuYg5xAi1Jefe+eKdoDyF7RdRJfbZZ5/5So6gRakxOpCCQbmyJCZGBOTXqcfGjGEwIK0B/bsfw5D+5zGkfwVGZF1ApwuPsor172VNx3SiOAE/Hrc6WzQxA7DklUfBaBKcqo2LOyyqxKSyoSsN4JbjkANYs/w7xnkpDSznIuagJNQ3OffOF++EN4LkicCAMnZcx9V1eVvAsM1OsGDBDaLlA+GL0HE3ZUcbSm7uAKSkvCVrV2dA/pi5awcUGzPn5Lo7AQApKWHo3/sSRmRV4ck8Bp3DT0JjOoJPVzgH4pstOhw5ZcD+wiT7cmDhyWQ0NFrjr1Qq4MUXR+CX479jyWALgLYFOWu1Kly50upU/9Gjl7BnTwXuuacna0lr8uReouMjZJtzHI+oKC2qqupw4YKzx6HRaGJl1cjPH4Np07bjxIlaMIx1GbF373fQ0GB18lCp0Gb/47FbxcSwl5ufe26YXcbCQnbS7zNn6jB27FZeZwu+YHMAHnsnPG13VrLbuT/l7KhIOnZcuHABEyZMwMSJE/H8888Lllu8eDE+/fRTfP7554iPj/e4oK7gizXzQFxrl+q33IBkqd195fbdm2OUkvIWzKYGDEqrsgYK96/AkP4VyOxbiYhwZ4ViUcXgh6PJ+HZ/PH68nuHifE0S6urE/Zv4grMBOKVwEoqP0utD0NBg4g0IlmtPcvzNEW7QryNcBwY5O2nL3b1arw/BgAGxovY5V4L/5SL0nHv6efOUvFxoZ2fPIDkTW7VqFS5fvow5c+aIlps9ezbefvttrF69Gn/96189JZ/PcHVdPpBsXHIRsm24GwQsN2DUnTFSWS637SrcWojv3/8W6T2rodHwpWOKwY/HuqL0Qk/MfOIB6+7C6hRM+9WHrOWrEBnrEVzlxNcHsQDfpiazYECwXHuS0Dlu0K8j3KU5OQHictttajIrtgHzySA3aF0unn4nvSVvMP7tCEQkX98vv/wSkydPRlSU+E6t0dHRmDJlCrZt2xaUSszVdXnubr3c40BEyLahZFdnQP6YuTS2DAOV+Rw0pjaFpWk9CrWFvfnjwD7WdExHThmsW4lcDxr+6UQyrtRFALB+4f527p2C8sgJmubOsvR6LYqL60TLOBIWpuGdifHJI2RPcvzNEW7Qr1oNpKToeXNRyumrHLuVrS4p+5wrwf/u4mkbm7fkDST7eDAjGVd/9uxZDBo0SFZlAwcORHFxsdtC+YP8/DEYPrwLevWKxvDhXWSvy3O9NwM9dg6w7mKs14cgJEQFvT7EHuwrdF4IuWMmWY5phbr1Z2gb3kdY/TxEXpqA6As9EF2Thcgrv0GY8WVom7+A2lIORhUBk3Y4miN+h4ZOr2HP+QIkjV6KzLufxow/P4B/vPML7NjXB+GR8YLtceUpKBiHlBT2HxDHjBNhYRqnawDGaYkpP3+Mffx0OmtWEsexzM8fY6+Xu+sy3/jk549BZmYsdDoNdDo1srLieMd48eL+rI8ni6XNSYdrY7HdY41GBbUaiIsLRXi4GqGhatE28vPHoF+/aHtmj4gIjb1PmZmx4D72BkO44DOh9DlTiqvvshDektfTcnZUJKcNKpUKFpmJ0CwWS1D8EefDVe+lS5eaRI8DkZycZN41/ZycZOzePdlubF648ICosVnumLHKWYzQmH6G5tr1QGFTITStx6FCs9N1FnX89US3mbBoM3CuuicezitF+fnG61kVQpGY2IT/ff0AGIbB9Olf4No1SBrJGYa99JiSoseqVWPszi4WC6DRqKHVtjk62K65dq0FJSX1TjOuLl0iMGVKHwwdmsAy1ufl9cecOXtw992fIzRUjR49onD+/DUAKrz22mEMG9bFSR7Hcfv22yksB4Dc3B1O2doffjgJra1sgfbvr8bmzadZWTNKSuqwcOEBdOkSodiRIDU1Gvv28QcHr1t3u6iDAt/uAwMGxNqPU1L0smSQi+PzVlJS53Y2HKH3xRNyzp8/FNOmfYVz564iJ2ezx5xGOhKSjh033ngjhg4din/961+SlT355JM4ePAg9u7d6zEB3cEXBlBvGX3dwZ1+e9IorjJXW5cBTUegvr7DsNp8BiqehLVmTU9YtBn27OxW+1UiHD/xxTKNA/KDdpVknReqn68MnzOG2BKjXNmlgtIjIjSstFOObTtmkveWg41Uvd4KqvdUcLs/ceXvBzl2sJGciY0bNw6rV6/GU089hb59+wqWO3XqFDZt2oTHHnvMowIGOtzs7LGxwb2u7ZKxmbFAbS6xK6q2dExVzkWhhTkknbVZo1k7CFB3Uiyb2HkxueVmnZdTFwDWUhC3rNQihivBv1xZm5v5ZZdyRvGUI4GcoGpH3Mk2r5RAd57wtpNLR0DSJvbEE08gMjISd911FzZt2gSTif21ajKZsGnTJkycOBFRUVF44oknZDdeVVWFxx57DL1794bBYMCIESOwe/du+++zZs1CTEwM61+gZdFPTIwQPQ42uMZl592Fm6FuPQxtw3qE1f0JkRfHI/pCKqJqshFZ+1uEXXsV2uavoLZUgVFFwRR6I5oj8tDQ6V+4Gv8N6hPLYUzYhcaY/0NL5GMw63JkKTBeWRzOC8ldUmKNV8rOLsDYsVtRWlrPW1bMWM9Xv42ICOt1kyZtw9ixWxEVpWX9LpXNna9uW5Ayty82uLKGhvLLzm1b8t66iLvye9OhwVt99hTcsfC0k0tHQFYC4B9//BEPPvggqqqqEBYWhj59+kCv18NoNOL06dNoampCUlISNmzYgMGDB8tquLa2Fr/4xS8wcuRI5OXlIS4uDqWlpUhMTES/fv0AWJVYZWUlVq9ebb8uNDQUnTt3ltVGoGSx9zXu9NuxP6ndLMj/ZzyS405bZ1ithVCbTkIFnvgrdeL1dEzW2ZVFmwmLJhVQubAnh4RsfJnGAWD69G0wGtk2Mb7lJL7g6fJyI6ZO3W63idl2drbZxABrNvnS0nrU1DRBpVIhPNy6S/PJk21/sLlZ4B99dCBmz9513cWeYc2OIiI02Lv3XgBATs5mwSU2qR0AHnkkCW++WYkzZ2px6VILAP5M8t56VktL692S31U55Dzngfh+OiK1wwUftJzIRpY/+JAhQ7B371689dZb+OKLL3Dy5ElcvXoVUVFRyMzMxPjx4/Hb3/4WnTrJ+6IGgNdffx2JiYksBdWjRw+ncjqdDgaDQXa9vqasrB7Hjl1GU5MZ1dUNKC83Kn5JXM0I4LglvJLtL5wzBYxGj5Rr0LQWIkF9AMue+Aa9UkqQmnTReoGDBzUDFcyaNPtyoOV6DkFG00WRrI880g9vvnmSV3ah8RBzJCkpqYPNP6Kurgk33rgJLS0WmM3sbzTrti3O322HD9ew/gjbHCUOH76EESM2oVMnLa5dMyE2NgzDhnVhZaB3xGg0Yd262zF16hc4cKAa+/dXIyJCg61bf4Unn9zFSfMVibKyekyatA2tnO+Cw4cvYtOmU3jyyV1obLTKotWqEBWlZcnf1GTCSy+dwtmz7I01uUuJ7maHELu+rKwe166xV2gcl+2492337nK33xm5BHq6KW85jXQkFG3F4klGjBiB2267DZWVldi1axcSExMxY8YM5Obm2j0cZ82ahc8++wyhoaHo1KkTcnJy8PzzzyMhwTl/HR/B4tjhah1CW8KL9psxI/e3a6BpPWLfrHFYRhVio686FW1qDkHx+W7oM3CUwyxrAKBW7k3GlZWLoxOCK+OhJLM74OxIIedabj1CWdq59QP82S2kHEqEnEK4jhFCOI6puw4OYtfzZQFxNRu+EjrqjKSj9lsIv0XmlpSU4I033sDjjz+OOXPm4MiRI/jzn/8MAMjLywMA3H777bjrrruQmpqKc+fOYcmSJZg4cSJ27twJnU7HW29RUZHosadpbDQ5HStt09U6+DJA2K4rKiqCCk0IV59BhPokIjSnEKE+hXB1ET5c5mzcNjFRaDD3Q8HWKHy734Afj3fFybMJSErSY8sWxy15KhX1TUhWvt9tsrsyHpWVzkrYEY3Gajt65JEkvPDCKUXXCrVXVFSEefN6Yv78Jly50oqYGC3mzeuJJ58sdCrf2GjiLTtlirDyFBoz7viIXW8bN24fbfLLRex6PnnmzespWL8n3hkb3n6/A5WO1G8phe03JWaxWDBkyBAsWLAAAJCVlYXi4mKsXbvWrsSmTJliLz9w4EAMHjwYGRkZ2L59OyZOnMhbr2OHffHFEh6+m/VVGR4eorhN1+vYaf+/zp0aMDi9AgOTmtBweQ9iwkuhNp2CCs7eTpUX47D/sFVR/XQiGUxoJta99wigUmHtJ+wv7qSkKI+MoVq9U3ImZmvHlfFISjrGylrBxWwGGhvNePPNSiQlRbHKco/lYBuXtDTg228zOb+ddaovPDwEo0dnOpUND/9B8UxMIKyM93rbuHHHR+l9Fbuee7/0emtfhfDEOwN03BlJR+23EJ6zuivEYDDYHThs9O3bF+Xl5QJXAElJSUhOTg6orCCeiOZXVAfDQGU6h5Cm/+DlP/0PH//rLZT87wVc/n4Bvn57NcKvPo847RfQmI7DYrGgydIXLWH3ozFqMYyxW1FvKEZN9I94cd2fsX7bvai8eisWL7vPHo+Vnz8GWVlx0OnU0Ok0aGoysTzNXMUxWwUAdO7M9uJ78cWR9v9fsWIUK7OFLcu5GPPnD0NEhMYpWwY39r60tB51dc2szBMLFtyAZctGOld6nbAwDQyGMOj1IejeXW93Dtm48RRiY/MRE5OP2Nh8bN582t5Xx+wWajUQGRli9450pKBgHLRabotAeLgaL744EuHhbYOm1aqQlRWHuDj2KoSQB6RjBghu9o/mZrNdlt27y5GS8hZiY639GDBgg5OsQtklSkrq0LUr2yM3JUUv+sxwn/cVK2528iCVQ3l5o0vXEe0Lv9nEZs6cifPnz2Pbtm32c0uWLMGnn36Kffv28V5z6dIlpKen4/XXX8cDDzwg2UbQf7EwJqhNp+wBw5rWQqhbj0DN1DoVvdagRXFFd5w61x3bd8bix+NdcbQoEZlZXRUbtn0VFOuIYxuutC90Ddf+wje7cbVtIZukHLmk5Bcry+1TRIQGarVK0s4kd4yk2peSW+61UnJJwc3eH2iBzN4i6P+ueRi/LSc+/vjjGDt2LJYvX47JkyejsLAQa9assW/3YjQasWzZMkycOBEGgwHnzp3DokWLkJCQgAkTJvhLbO9hMUJjOmZPdqs2HYGm9ZhAOqY4mEMycaWxH/65Bvj+pyQYm7tj9Zrb8KfZ21gecK4Ed/oqKFboN1faF7qGG4xudRpy9lh0pW2hrPRy5JIqJ1aW26dOnUJQU9PCKsMXNCski1CArSvjruRaKbmk4GbvD7RAZsI3+E2JZWdnY8OGDVi0aBFeeeUVpKSkYN68eZg5cyYAQKPR4NixY/jggw9QV1cHg8GAUaNG4a233pLMqB/oqFuPQWM6DpW57LrSOgK1+bRAOqYePOmYkgCVCjoAcxexy3siM7a3smuLZTx3bMOV9oWuSUyMwLlzbX/ww8Ods7i72jZ3Vse3rOdqpn+xstw+xcfrcPWqWTLTupAsQpntue3zudkL3VMlz4yrzxs3e3+gBTITvsFvy4m+wO/TbsaMkJZvoG3YgNCmzeJFEQILbzqmGEVNlpbW8wb9Kq3DW0GxjjsSM4w1Ezy3DVfaF+o3X6DtvHnf48SJKwBU6NevE9av/6VLbW/efBq5uTtE4/Tk1ldaWo8ZM/4rKJdYnfPm9YRWGy8ZNCskiy3gtqHBqshsO31zZRULHOcLQJf7zLj6vO3cWYgXXzwbsIHM3sLvf9cCDFJibmL7Om0w1uK+O47gj787Ar2a36bHpTkiry1oOCQdUPGHDQghFBQst998X9b791fZg5IBoHfvaDCMsELk1pGbOwB/+MNu2YHbntqifffuctx335dobbUo3kK+pKQOU6d+gVOn6sEwVhvTxo13ICcnmSWfVguUlV1DS4u1jRUrbsbrrxfi5Mk6AIw9u0dqarT93nD3EAsNVcNiYezXv/LKIVbGD5UKSEriVyJC2O63JwKaZ8z4H06e5Fek2dkFrBlTr17R2LLlDqcM9Y4Z9vlkcDW4X6jffP3wxDPlS5TITEqMDSkxF1CZzyO04X1oGzdAYz4r+7qWsEloDX8QJt0YQOV+jjShoFG5/eb7sj54sFrQFV7OtvXOm0eKB7J6yonEnQBaPucE2/VizihiTiJCzhJS1/PVJYXtfns6oJlbh5zAbjkZ6r0d7Bzomev5UCIzKTE2gb8NsT9hGGhaf4K2cQNCGzZABXmxRIwqGi0RD6Il/EFYtPI2FHUFdzNg8xnUxf6oynE+4F4vJZOnnEjcGQu+Nm3XK3HscCwvp31XMtwrKe/u9dxzfDknJ03axiovJ0O9tzO3B3rmej6CUeZAwW9xYoGIpnkXIi4/iE6VMdZ/VZ2hvzQGuoa1ggqs+Hwy5r58JxJHPQ9V/1cw4qF3UJ94Dk3RS72qwAD3M2DzZfgWy7outG29I9zrpWTyVJZxd8aCr03b9VKOHVxsGdw9kY1c6Vi4O5ZS95cv56QrGeq9nbk90DPX8xGMMgcKpMSuE9L0JfSX74K2+TPBMq2ho9EQk4+6xArUJdWiLqkWxoTvsevn+xEZneLzLcbdDbTmC2DlBiX37h0tun06Xx1KZPLUFu0FBePswc5Kx4IbnBwRobFfb5NPp2P/obU5cmRmxrICqo1GE3Jzd9jvDVfRcYOvhf6A6/UhisfC3bHkBrpnZsay6sjL24n9+6tRXFyP/furkZu7w6nNgoJxkjJ4IkGAVD888Uz5kmCUOVAgm9h11K0/I+pijv24JXwaWsJ/A3PoTR7dTsQTSBmBvbVmztcuwzDIy9vpsneaEoO2WNmSkjpMn/4Frl1jO6FYHT62o7HRApUK6NcvBgUF41htONar14dApVLh6tVWVj18Tg2HDk1FSUkdRozYzNqYMiUlEsnJkdedQVQoLTUKLpmFhqqhUqmcNrbU6dTYv/8++/iKjY8r99sV5wehMXAXVx0xOqptqKP2WwhSYkGIlBHYW/2Wm7GdTyYldQpdJ1ZWSTYKbhtyMokI1S/kFCIny7wUQhn3uePjyv32ZFYUd3G13vb6fkvRUfstBDl2BCH+MgIraVeuTO7UKSfThpLMFWJt8jk1AEBFxTVWea3WOaOGqyjJ4mFD7qzGlWdIaAzchZwaCHcgJRaEeCujhqvtupOxQUlfxMoqyUbBZ0SXyiQitLlibS07LZhOF+KUUcNV+MZXalxtdivbdbm5O3jlduUZ8tYGk/56non2AS0nuok/AittGQ6qqhpw+XITYmJ0SE6OtLct1W8hmW1BqI2NZjAMg4SEMKSmRgtmv8jPH4OysquYNu1LyWwP3DZtQbFnztTi0qW2vH8GQzhSU6NY9jaxbe1tZUpL61Fd3QS1WgW1moHJZN2yRKWyOmGYzRC0iZWW1uPOO7eivLxtBtCpkxZmM4Po6FDU17dAp1Oz5ASs9cXH61BT06bIunfX49NP77SPU2ioCiUlwjaxkBAVTKa2V5AbXA5ANJvFxo2nkJf3jeD2LFy7Vds9NsF8XSS1Wvi+ycWd98DVjB07dhRi6VLK2NHRISXmJv4MrBRqW6rfnsxm7mqWdim7kVx7kJwdnfmyy3OJickX/V2sbsd4LyVjxXe90sBfqR2zufJIBWG7+vz64z2gLPYEQMuJbuPP9XxX2/ZkNnNXs7QrDYJWco6LVFCxO6hUKgwfniDLTsSVVadTo7WVLZzSwF+xvul0Gid5XBl3OfjjPaAs9gRAcWJu488gRVfbFrpOKGZJyk4lp6xUUCxfeTl1y+mzWAC3u4SHa/Dll3fj0KGp+PLLu0WXs7iyZmXFIyKC/R2pNPBXrG/p6TFO8kjVf/78NZc2mPTHexATw95NlGxpHRNSYm7izyBFV9sWus4WhKrRqKBWAwZDmGS9cmUQCoo1GKxZQmy2q+TkCFbgtVTdtjIpKZGIiNCge3c9evaMZO0MLWdcuDs7P/PMYAwf3gXdu+uh14cgLi7U6RqVCoqDqrn98UTAupAi48uwwb3HcXGh0OtDEBpqraS52WwPZFYqh6/fg8WL+1OAMEE2sfaIr/rtaMxvaWlhOUbExenQu3cnnxrb+frNzVAPAFqtCmq1ipV1Xo5jgpBziu04L68/5szZw8rO3rVrpOIs70qdJIqKinDhQjgmTGDnMTQYwnDy5HRZYxdogcxyoPebAEiJtUt81W85ThW+NLbz9VtKRqlgZrG6uM4pfE4aAwbEKs7yrtRJoqioCGPG7HZy2JDj0OJqm3LxpsMHvd8EQMuJhBu4s329r5Bq3/a7HMcEKecUvgz+UtfIaUfOGPI5bKi4iRpF8NZyIDcYnHtMEO5CSoxwGTmGdH8b26Xat/3uihMJ10mCL4M/9xqtlq1Y9HpnB2FXnCT4HDbCw+U7idgCmeU4qCiBGwzOPSYIdyEXezdxxZbiKbuAN+p13HVXq1UhNTUaLS0W3vod0xCJ2cTE5C0rq3fa5ZdhLKJB13y0Bf3utDtzTJnSB/PnD8OkSZ+j1cEbW6NRISREhZ49o9HUZMKgQe/hwgX2DGH//mp07pxvD45mGAY1NW1l1GrgueeGYcuWYpSW1qOmxnm2ZDSa7EtpcXGh6NpVj1On6lhlVCqV07jwBXVL8de/DsWzz7btKB4aqrY7iUg9J1I7OrsDNwVXbKz3P2pceS+CcTdowgrZxNzEFVuKt+0N7vRbLBjWXbn55D127LLTLr8AFAddc4N+bfYgsd2K5dj0bOUB50THtsBkubs4DxvWxamOXr2iER8f5tbzwWcTcwyalnr+pHZ0dgd/2MQCKbmxNyCbGBuaibmJK7YUT9mJvFGvWDCsu/Xzyatkl18luyzbjsXuh6vB4TZsssrdxZmvnvj4MI/cR7FxlKrf1SByOXgrabAYrownJSEOXsgm5iau2FI8ZSfyRr1iwbDc+ktK6jB27FZkZxfICpDlk5dvl1+xoGuhNrn2KNux2P1QEhwutkOxnABltdq5PdvGl564j2K7JUvV72oQuRy8ZWsTg69/Us8q7awcvJAScxMlAbme9vzyRr2Owbfh4Wr0799ZsH6+nX6VyssX7CsWdC3UpmPQr2OAs233ZZ1OA51OjaysONZvtkBpnc5512WVypr1whZ4LbT7M9+YOQZH2+Th9n/Pnil2W5+791EsaFqqfqkdnYMNvv5KPau0s3LwQjaxdoiv+u2tAFlX2hTa2dmbBnupnablOLHk5CS7LYec+90eHReUPOf+eFa9RUf9uyYEzcQIl/HHEoxQm3l5O3HkSL3Tl7bS2aISxOrm+23atK9gNJpgMjEwGk2YOnW7x2RxR9aOAC0Xtl9IiREu448lGKE2hQzz3jTYK91pWokTi6fp6I4LtFzYfiHvRMJlvLXTryttCu0O7M1dg5XuNF1d3cBygxdyCPHG0l9H3z3ZH88q4RtoJka0C6wOHNFOX9re/AIXq9udjPXeWPqjmQjRXqGZmJv402DubtuO2Tn4sq7r9SFQqVS4erXVp31zpV8Mw9iz1Dc1mTBjxn8Vye3YJl+/GYbhlenLL++2Xztp0jan37gy9urVCSdPXkFrK4Nnn93Lyoxhq+fQoRrWdZ5Y+isrq8exY5fR0GBCSUk9xo//FCkp+nbh4EF0bMg70U38GenvbsYObqYJvqzrjviqb57IuOCIJ64HICiTXHmlMmMIySAni73U/RbKKhLImSmk6Kheeh2130LQTMxN/Gkwd7dtPkcDsTp81TdPZFzw9vVSThyeqCckRIXs7ASPLP0JOZF0NAcPov1BNjE38afrrrtt82V5EKvDV31zpV/uyi11vZhMcuWVyozB/T07O8FjWS7EsqAQRDBDMzE38UduOE+1XVAwDlOnbmfZxFJS9PY6+WxDjvDZroRsR0qYP38Ypk370i7XggU3yBqL6dO3wWiEqNx8dsCcnGTWWApdLzTWcu9Dfv4YzJjxX5w40ZYt3pV6XMF2rxsarEuKiYkRdpsYQQQzZBNrh/hrZ2cp25E79cqpwxXbkGO292CFnvOORUfttxA0EyNcxlvZ0L1pZ/RnwDFBEJ7HrzaxqqoqPPbYY+jduzcMBgNGjBiB3bt3239nGAZLly5Feno6EhMTceedd+L48eN+lJhwhM8W5AkboTftjGLZ3gmCCD78psRqa2sxbpx1x9wPP/wQ+/btw8svv4yEhAR7mddeew0rV67ESy+9hK+//hoJCQmYNGkSrl696i+xCQf4Amg9EVTrzcBcuQHH/kTpFjcE0ZHxm01s0aJF2LNnD7Zv50+CyjAM0tPTkZubi7lz5wIAGhsbkZaWhsWLF+Phh6XtGErWjl0NHPZEsLOnA6bl9vv//u8w5s3bbz9etmwk7rije9BmO7f125Pjya1r/vxhWLTo4HUHEC0AazJfbjtSgeRi5S0WhrXJp5RNcMeOQixdejYg7pkvg/87qm2oo/ZbCL8psREjRuC2225DZWUldu3ahcTERMyYMQO5ublQqVQoKSnB4MGD8fXXXyM7O9t+3f3334/Y2FisWrVKsg0lN9tVZwJPBDt7OmBabr9jYvKdzg0f3iVotmnnYuu3J8eTW5deH8IbNMxtR04guVh5R6S2DRk1qgBHjrTN1vx5z3wZ/N9R/5h31H4L4TfHjpKSErzxxht4/PHHMWfOHBw5cgR//vOfAQB5eXm4cOECALCWF23HlZWVgvUWFRWJHgtRWXnV6VjOta5e5+k6uLh6vTdk8SVFRUUe7QO3rsZGfkXDbYdbrrHRJCqXWL2RkeL3s7a2VbBeX+Pr5yeYnk1P0pH6LaWw/abELBYLhgwZggULFgAAsrKyUFxcjLVr1yIvL89eTsXZbpdhGKdzjjh2WMkXS1LSMZSVNTocR8m61tXrPF2HI/L7vZNHliiPyuJLbP325Hhy6woPF56JObYTHr6bVS48PER0bLnl1WqgR49oWUtyMTGHAuaeefpZFqOjzkg6ar+F8Jtjh8FgQL9+/Vjn+vbti/LycvvvAFBdzc4ld/HiRafZmSdw1Zkg0B0ZxFi2bKTTcXvIdu7JPnDrKigYZz/OyopDZmYsbzt8DiRicnHLf/rpnTh0aKqsjB2LF/cPmHvWHp4fIrjwm01s5syZOH/+PLZt22Y/t2TJEnz66afYt2+f3bEjLy8PTz/9NACgqakJaWlpWLRokccdO9oT1O+OBfW7Y9FR+y2E32Zijz/+OA4cOIDly5ejuLgYH3/8MdasWYOZM2cCsC4jzpo1CytWrMAnn3yCY8eO4fHHH0dkZCTuvfdef4lNEARBBBB+s4llZ2djw4YNWLRoEV555RWkpKRg3rx5diUGALNnz0ZjYyOeeeYZ1NbWYujQodiyZQuioqL8JTZBEAQRQPg17dS4ceMwbpxwsKlKpcKzzz6LZ5991odSEQRBEMECbcVCBByUsYIgCLmQEiMCjry8ndi/vxrFxfXYv78aubk7/C0SQRABCikxIuDw527ZBEEEF6TEiIDDn7tlEwQRXJASIwIOCpglCEIutCkmEXCkpkYHTdJhgiD8C83ECIIgiKCFlBhBEAQRtJASIwiCIIIWUmJuQoG5BEEQ/oOUmJtQYC5BEIT/ICXmJhSYSxAE4T9IibkJBeYSBEH4D1JibkKBuQRBEP6Dgp3dhAJzg4uSkjrk5e3ExYtNiI8PQ37+GKSmRvtbLIIgXIRmYkSHghxxCKJ9QUqM6FCQIw5BtC9IiREdCnLEIYj2BSkxN1Ea7EzB0f6FHHEIon1Bjh1uYrOxAEBxcT1yc3eIOnooLU94FnLEIYj2Bc3E3ESpjYVsMgRBEJ6DlJibKLWxkE2GIAjCc5AScxOlNhayyRAEQXgOsom5iVIbC9lkCIIgPAfNxAiCIIighZQYQRAEEbSQEiMIgiCCFlJiBEEQRNBCSqwDQ9lDCIIIdkiJdWAooztBEMEOKbEODGUPIQgi2CEl1oGh7CEEQQQ7pMQ6MJQ9hCCIYIcydnRgKHsIQRDBDs3ECIIgiKDFb0ps6dKliImJYf3r27ev/fdZs2Y5/X777bf7S1yCIAgiAPHrcmJaWhr+85//2I81Gg3r99GjR2P16tX249DQUJ/JRhAEQQQ+flViISEhMBgMgr/rdDrR3wmCIIiOjV9tYiUlJejfvz8yMzPxyCOPoKSkhPX73r170adPHwwdOhRPPfUUampq/CMoQRAEEZCoamtrGX80/NVXX8FoNCItLQ0XL17EK6+8gqKiInz//feIjY3F5s2bER4ejtTUVJw7dw5LliyBxWLBzp07odPp/CEyQRAEEWD4TYlxMRqNGDx4MObMmYMnnnjC6ffKykpkZGTgzTffxMSJE/0gIUEQBBFoBIyLvV6vR3p6OoqLi3l/T0pKQnJysuDvBEEQRMcjYJRYU1MTioqKBB05Ll26hMrKSnL0IAiCIOz4TYk999xz2L17N0pKSnDw4EE89NBDaGhowAMPPACj0YjnnnsO+/fvR2lpKXbt2oVp06YhISEBEyZM8JfIBEEQRIDhNxf7iooKzJw5E5cuXUJ8fDyGDRuGr776Ct27d0djYyOOHTuGDz74AHV1dTAYDBg1ahTeeustREVF+UtkgiAIIsDw20zszTffxIkTJ1BTU4Pjx49j/fr1SE9PBwCEh4djy5YtOH36NGpqanD06FH8+9//RkpKilM9VVVVeOyxx9C7d28YDAaMGDECu3fv5m1z9uzZiImJwT//+U+v9s0XyOn36dOn8Zvf/Abdu3dHUlISbrnlFpw8edJPEnsGqX4bjUY888wzGDBgABITEzFs2DCsXLnSjxK7T0ZGhlP2mpiYGNx///0AAIZhsHTpUqSnpyMxMRF33nknjh8/7mep3Ues362trViwYAFuuukmJCcno1+/fpg5cybKysr8LbbbSN1vR9rT3zRXCeoEwLW1tRg3bhxGjhyJDz/8EHFxcSgtLUVCQoJT2a1bt+LQoUNISkryg6SeRU6/S0pKMG7cOEybNg2ffPIJYmJicOrUKURGRvpRcveQ0++//vWv2LlzJ1atWoXU1FR89913mD17NuLi4jBt2jQ/Su86O3bsgNlsth9XVVVh9OjRuOeeewAAr732GlauXImVK1ciLS0NL7/8MiZNmoQDBw4E9cqFWL8bGhpw+PBhzJ07FxkZGaivr8dzzz2He++9F3v27EFISPD+aZO63zba0980dwjeOw3g9ddfR2JiIis1VY8ePZzKnTt3Dn/5y1/w8ccf49577/WhhN5BTr+XLFmCW2+9FS+88IJgmWBDTr/379+PqVOn4pZbbgEApKamYv369fjhhx+CVonFx8ezjtevX4+oqCjcc889YBgG//73vzFnzhzcfbd1R4J///vfSEtLw6ZNm/Dwww/7Q2SPINbviIgIfPzxx6zf//GPf2DkyJE4efIkBg4c6ENJPYtYv220t79p7hAw3omu8Nlnn2Ho0KF4+OGH0adPH9x8881Ys2YNGKYt9M1kMmHmzJmYO3cu+vXr50dpPYdUvy0WC7744gv069cPU6ZMQe/evTFmzBhs2bLFz5K7h5z7PXLkSHzxxRcoLy8HAOzbtw9Hjx7Fbbfd5i+xPQrDMFi/fj2mTp2KiIgIlJaW4sKFC7j11lvtZcLDw3HTTTdh3759fpTUs3D7zcfVq1cBADExMT6UzLvw9bs9/k1zh6BWYiUlJXjjjTfQo0cPbN68GY899hgWLlyI/Px8e5mlS5eic+fO+N3vfudHST2LVL9rampgNBrx97//HWPGjMFHH32EKVOmIDc3F1988YWfpXcdOff7pZdeQkZGBgYNGoT4+Hjceeed+Nvf/oY77rjDj5J7jh07dqC0tBTTp08HAFy4cAEAnJbQExISUF1d7XP5vAW331xaWlrw3HPP4Y477kDXrl19LJ334Ot3e/yb5g5BvZxosVgwZMgQLFiwAACQlZWF4uJirF27Fnl5edi9ezfee+897Nq1y8+SehapflssFgDAr371K3v2k8zMTPz0009Yu3Zt0P5Bl+o3AKxevRr79u3D+++/j27duuG7777D888/j+7du7eLrXzeeecdZGdnIzMzk3VepVKxjhmGcToXzAj1G7DOTPLy8lBXV4f333/fD9J5D26/2+vfNHcI6pmYwWBwmk737dvXvpS0a9cuVFVVoV+/foiLi0NcXBzKysqwYMECDBgwwB8iewSpfsfFxSEkJES0TDAi1e/GxkYsWrQICxcuxPjx4zFo0CDk5eVh8uTJ7cJ7q6amBp9//jkeeugh+zlb8D931nXx4kVeB6dghK/fNkwmE373u9/h559/xtatWxEbG+sHCb0DX7/b6980dwjqmdjIkSNx+vRp1rnTp0+jW7duAICZM2fajd02pkyZgilTpvC+EMGCVL9DQ0ORnZ2NoqIiwTLBiFS/W1tb0dra6rQvnUajsc9Og5n33nsPOp0OkydPtp9LTU2FwWDAjh07kJ2dDcCa/Wbv3r1YtGiRv0T1KHz9Bqz3+5FHHsHx48fxn//8p91l8+Hrd3v9m+YOQa3EHn/8cYwdOxbLly/H5MmTUVhYiDVr1uD5558HYLULcL9GbXuYpaWl+UNkjyDVbwB46qmn8PDDD+Omm27CLbfcgl27dmHLli3YsGGDHyV3D6l+R0dHIycnBwsXLkRkZCS6deuGPXv24IMPPsDChQv9LL17MAyDdevWYfLkySy3eZVKhVmzZuHVV19FWloa+vTpg+XLlyMyMrJdeK0J9dtkMuGhhx7Cjz/+iPfffx8qlcpuH4yOjkZ4eLi/RPYIQv1ur3/T3CFgsti7yvbt27Fo0SKcPn0aKSkpyM3NxaOPPipoD8jIyEBeXh6efPJJH0vqWeT0e8OGDfj73/+O8+fPo1evXvjjH/8Y9H/YpPp94cIFLFy4EDt27MCVK1fQrVs3zJgxA0888URQ24i+/fZbTJw4Ef/73/8wdOhQ1m8Mw2DZsmV4++23UVtbi6FDh2L58uXtYnlJqN+lpaXIysrivWblypV48MEHfSWiVxC731zay980Vwl6JUYQBEF0XILasYMgCILo2JASIwiCIIIWUmIEQRBE0EJKjCAIgghaSIkRBEEQQQspMYIgCCJoISVGEB4kIyMDs2bNkixXWlqKmJiYoA4+J4hAgJQYQfCwYcMG+4663377LW+ZW2+9FTExMbjhhhs80ubSpUtZO/kmJCQgMzMTzz77LOrq6jzShrc4duwYli5ditLSUn+LQnQwgjrtFEF4m7CwMGzcuNG+yaaNM2fO4NChQwgLC2OdP3jwINRq974NX3nlFURHR+PatWv4+uuv8e9//xs//vgjtm3bFrBZR44fP46XXnoJN998M1JTU/0tDtGBoJkYQYgwduxYbN26Fc3NzazzBQUF6NKlC4YMGcI6r9PpoNVq3Wpz4sSJmDp1Kh555BG8++67uOuuu/D999/j4MGDbtXLMAyamprcqoMgAg1SYgQhwpQpU2A0Gp02E920aRMmT57sNOvis4lVVVXhoYceQkpKCnr06IEnnnjCvguxHGyzwJKSErS0tOCFF17A6NGjkZqaisTERNx22234/PPPna6LiYnBH/7wB3z88ce46aab0KVLF2zevBkAUF9fj+eeew4ZGRno0qULBg0ahL/97W9OytpWx1dffYVRo0bBYDAgOzsbmzZtspfZsGGDfYPGu+66y74cSvY+whfQciJBiJCcnIycnBxs3LjRvgXGwYMHUVxcjPvvvx9HjhwRvb6pqQl33303iouLkZubi9TUVHz66ad47LHHZMtw9uxZAEBsbCyuXr2Kt956C5MmTcJvfvMbNDY2YuPGjXjwwQexadMm3Hbbbaxr9+7di61btyI3NxcGgwF9+/ZFY2MjJkyYgNLSUvz2t79Fz549ceTIEfzrX//CqVOn8N5777HqOHDgAD777DM8/PDDmD59OtatW4e8vDxkZGSgX79+yMnJQW5uLvLz8/H000+jb9++AIARI0bI7iNBuAopMYKQ4L777sPcuXNRW1uLmJgYFBQUoHfv3vb9u8R45513cPLkSaxatQrTpk0DYN0TasKECYLXXLlyBSEhIXab2JtvvgmDwYAbb7wROp0OP//8M3Q6nb38o48+ilGjRuGf//ynkxI7efIkvvnmG9aOyK+++iqKioqwc+dO1iaj/fv3x9y5c/Hdd9/hpptusp8/ceIE9uzZYy97zz33YNCgQXj33XexePFi9OjRAyNHjkR+fj5Gjx6NUaNGSY4LQXgKWk4kCAnuvvtuqFQqbN26FSaTCR9//DHuu+8+Wddu374dCQkJrPIajQaPPvqo4DUjR45E7969kZmZiTlz5iArKwubNm1CREQENBqNXYG1tLTgypUruHr1KnJycvDTTz851TVixAiWAgOAjz76CCNGjEB8fDwuXbpk/zd69GgAcPLGHDVqFEvZdenSBWlpaSgpKZE1BgThTWgmRhASdOrUCWPHjsWHH36I5ORk1NTUyFZiZWVl6Nmzp9Nu03369BG85u2330ZMTAwiIiLQrVs3JCUlsX5ft24d/u///g8nT54Ew7TtpMTnudijRw+nc2fOnMHRo0fRu3dv3vYvXrzIOubbDTwmJgZXrlwR7ANB+ApSYgQhg/vuu8++/fvQoUMFFQAXhmF4lYuj8uFy4403wmAw8P62adMmPPXUUxg/fjxmz56NhIQEhISEYMOGDdi4caNTeb4dji0WC2655Rb88Y9/5G0jOTmZdcxVwHL6QBC+gpQYQchg3LhxiI6Oxp49e7Bs2TLZ13Xv3h1HjhyB2WxmKYMzZ864JMeWLVvQo0cPvPfee067eMulZ8+eMBqN9uVDTxCo8WtE+4dsYgQhA51Oh1dffRV//vOfce+998q+buzYsaipqWHNksxmM1avXu2SHDZF6DgLKikpwX/+8x/ZdUyePBmHDh3idctvbGyE0WhULFdERAQAoLa2VvG1BOEONBMjCJkoUV42HnroIbzxxht48sknUVhYiB49euCTTz5RFCfmyPjx4/Hpp5/igQcewPjx41FRUYE33ngDvXv3xtGjR2XV8eSTT+LLL7/E9OnTcf/992Po0KFobm7G6dOn8dFHH2Hjxo2KU2llZWVBrVbjH//4B+rq6hAeHo6hQ4fy2uQIwpOQEiMILxIeHo6tW7fiL3/5C9555x1otVpMmDABjz32GG6++WbF9f3617/GxYsX8cYbb2Dnzp3o1asXXnzxRRQXF8tWYuHh4fjkk0/w2muvYcuWLdi8eTMiIyPRo0cPzJo1C2lpaYrlSkpKwooVK7BixQrMnj0bZrMZK1euJCVGeB1VbW0tWWcJgiCIoIRsYgRBEETQQkqMIAiCCFpIiREEQRBBCykxgiAIImghJUYQBEEELaTECIIgiKCFlBhBEAQRtJASIwiCIIIWUmIEQRBE0EJKjCAIggha/h+iVctsLymhuQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "scatter_fit(heights, 'MidParent', 'Child')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fitted values range from about 64 to about 71, whereas the heights of all the children are quite a bit more variable, ranging from about 55 to 80.\n", "\n", "To verify the result numerically, we just have to calculate both sides of the identity." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.320949896063959" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "correlation(heights, 'MidParent', 'Child')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is ratio of the SD of the fitted values and the SD of the observed values of birth weight:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.32094989606395896" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.std(heights['Fitted Value'])/np.std(heights['Child'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The ratio is equal to $r$, confirming our result.\n", "\n", "Where does the absolute value come in? First note that as SDs can't be negative, nor can a ratio of SDs. So what happens when $r$ is negative? The example of fuel efficiency and acceleration will show us." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.5060703843771185" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "correlation(hybrid, 'acceleration', 'mpg')" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.5060703843771186" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.std(hybrid['fitted mpg'])/np.std(hybrid['mpg'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The ratio of the two SDs is $\\vert r \\vert$. \n", "\n", "A more standard way to express this result is to recall that \n", "\n", "$$\n", "\\mbox{variance} ~=~ \\mbox{mean squared deviation from average} ~=~ \\mbox{SD}^2\n", "$$\n", "\n", "and therefore, by squaring both sides of our result,\n", "\n", "$$\n", "\\frac{\\mbox{variance of fitted values}}{\\mbox{variance of }y} ~=~ r^2\n", "$$" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.12" } }, "nbformat": 4, "nbformat_minor": 2 }