{ "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", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "plt.style.use('fivethirtyeight')\n", "\n", "import functools\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from scipy import optimize\n", "\n", "def minimize(f, start=None, smooth=False, log=None, array=False, **vargs):\n", " \"\"\"Minimize a function f of one or more arguments.\n", " Args:\n", " f: A function that takes numbers and returns a number\n", " start: A starting value or list of starting values\n", " smooth: Whether to assume that f is smooth and use first-order info\n", " log: Logging function called on the result of optimization (e.g. print)\n", " vargs: Other named arguments passed to scipy.optimize.minimize\n", " Returns either:\n", " (a) the minimizing argument of a one-argument function\n", " (b) an array of minimizing arguments of a multi-argument function\n", " \"\"\"\n", " if start is None:\n", " assert not array, \"Please pass starting values explicitly when array=True\"\n", " arg_count = f.__code__.co_argcount\n", " assert arg_count > 0, \"Please pass starting values explicitly for variadic functions\"\n", " start = [0] * arg_count\n", " if not hasattr(start, '__len__'):\n", " start = [start]\n", "\n", " if array:\n", " objective = f\n", " else:\n", " @functools.wraps(f)\n", " def objective(args):\n", " return f(*args)\n", "\n", " if not smooth and 'method' not in vargs:\n", " vargs['method'] = 'Powell'\n", " result = optimize.minimize(objective, start, **vargs)\n", " if log is not None:\n", " log(result)\n", " if len(start) == 1:\n", " return result.x.item(0)\n", " else:\n", " return result.x" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [ "remove_input" ] }, "outputs": [], "source": [ "\n", "def standard_units(any_numbers):\n", " \"Convert any array of numbers to standard units.\"\n", " return (any_numbers - np.mean(any_numbers))/np.std(any_numbers) \n", "\n", "def correlation(t, x, y):\n", " return np.mean(standard_units(t[x])*standard_units(t[y]))\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", " \"\"\"Return the height of the regression line at each x value.\"\"\"\n", " a = slope(table, x, y)\n", " b = intercept(table, x, y)\n", " return a * table[x] + b" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Least Squares Regression\n", "In an earlier section, we developed formulas for the slope and intercept of the regression line through a *football shaped* scatter diagram. It turns out that the slope and intercept of the least squares line have the same formulas as those we developed, *regardless of the shape of the scatter plot*.\n", "\n", "We saw this in the example about Little Women, but let's confirm it in an example where the scatter plot clearly isn't football shaped. For the data, we are once again indebted to the rich [data archive of Prof. Larry Winner](http://www.stat.ufl.edu/~winner/datasets.html) of the University of Florida. A [2013 study](http://digitalcommons.wku.edu/ijes/vol6/iss2/10/) in the International Journal of Exercise Science studied collegiate shot put athletes and examined the relation between strength and shot put distance. The population consists of 28 female collegiate athletes. Strength was measured by the the biggest amount (in kilograms) that the athlete lifted in the \"1RM power clean\" in the pre-season. The distance (in meters) was the athlete's personal best." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "shotput = pd.read_csv(path_data + 'shotput.csv')" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "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", "
Weight LiftedShot Put Distance
037.56.4
151.510.2
261.312.4
361.313.0
463.613.2
566.113.0
670.012.7
792.713.9
890.515.5
990.515.8
\n", "
" ], "text/plain": [ " Weight Lifted Shot Put Distance\n", "0 37.5 6.4\n", "1 51.5 10.2\n", "2 61.3 12.4\n", "3 61.3 13.0\n", "4 63.6 13.2\n", "5 66.1 13.0\n", "6 70.0 12.7\n", "7 92.7 13.9\n", "8 90.5 15.5\n", "9 90.5 15.8" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "shotput.head(10)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAEfCAYAAADPxvgvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxE0lEQVR4nO3deVhTV94H8G8AUURoEDGiAm5UwKVVURAlWm1Rq9QFq1hrLYq4taP2RetuXUHHWq0P4l60BdShjHutM68b4vpaBetWFEFhRMESFFRkyftHhtQIhESy5/t5Hp5p7rnJ+eU3Pvw49557jkAikUhBRERkhCz0HQAREdGbYhEjIiKjxSJGRERGi0WMiIiMFosYEREZLRYxIiIyWixiRERktFjEiIjIaLGIqSktLU3fIRgE5kGGeZBhHmSYBxld5oFFjIiIjBaLGBERGS0WMSIiMlosYkREZLRYxIiIyGixiBERkdGy0ncARETGJuNJCcJOSZBXXIZGdS2xpZcQbnZ19B2WWeJIjIhITWGnJLiQ+xLpT8pwIfclJpyU6Dsks8UiRkSkprziMqWvSXdYxIiI1NSorqXS16Q7LGJERGra0kuIbk7WaGVviW5O1tjSS6jvkMwWJ3YQEanJza4Ojg5y0ncYBI7EiIjIiLGIERGR0WIRIyIio8UiRkRERosTO4iISON0taoJR2JERKRxulrVhEWMiIg0TlermrCIERGRxulqVRMWMSIi0jhdrWrCiR1ERKRxulrVhEWMiOgNVDX7jnRPr5cTk5OTERwcDE9PTwiFQsTGxiq0FxYWYubMmfDy8kKTJk3g7e2NqKgoPUVLRPQX7ilmGPQ6EisqKoKXlxdGjRqFSZMmVWqfN28eTpw4gY0bN8LNzQ1nzpzBtGnT4OjoiODgYD1ETEQkwz3FDINeR2IBAQFYuHAhBg8eDAuLyqFcuHABI0eOhFgshpubG0aNGgVvb29cunRJD9ESEf2Fe4oZBoOenejr64sjR44gKysLAHD+/Hn8/vvv6Nu3r54jIyJzxz3FDINBT+xYuXIlZsyYgfbt28PKShbqqlWr0L9/fz1HRkS6pqtljFRV1ey7tBw9BWPGDLqIbdq0CefPn0d8fDxcXFxw5swZLFiwAK6urnj//ferfE9aWprW49JFH8aAeZBhHmS0nYdxKXVx9anskl06yjDm1wfY9k6xVvt8E/z3IKOpPLi7uyttN9gi9vz5cyxZsgQxMTEYMGAAAKB9+/a4evUq1q9fX20Rq+kL11ZaWprW+zAGzIOMqeZB3VGPLvJQlJoD4K/JE4UCa7i7u2q1T3WZ6r8HdekyDwZ7T6ykpAQlJSWwtFS8WWppaYny8nI9RUVkHgxx+jgnUlBV9DoSKywsRHp6OgCgvLwcWVlZSE1NhYODA1xcXNCjRw8sXrwYtra2cHFxQXJyMnbt2oXFixfrM2wik2eI08e39BJiwkk+XEyK9FrELl++jMDAQPnriIgIREREYNSoUYiOjsb27duxePFihIWFIT8/Hy4uLpg3bx7CwsL0GDWR6WtU1xLpr1y6M4RRjzaWMTK0ySKkPr0WMX9/f0gkkmrbRSIRNmzYoLuAiAiA+Yx6Ki6bArLJIhNOSnSy3h9pjsFO7CAi/dHV4q36ZoiXTUk9Bjuxg4hI2zhZxPixiBGR2eKqG8aPlxOJyGyZy2VTU8aRGBERGS0WMSIiMlq8nEhElfD5KTIWao/E7t+/j127dmH9+vXyLVJKS0uRm5uL0tJSjQdIRLpniMtOEVVFrZHY3LlzsXnzZpSVlUEgEKBjx45o3rw5nj17hs6dO2P27NmYOnWqtmIlIh3h81P6w1GwelQeiX3//feIjo7G1KlTsXfvXkilUnmbvb09Bg4ciIMHD2olSCLSLT4/pT8cBatH5SK2Y8cOjBgxAosXL0aHDh0qtbdr1w537tzRaHBEpB98fkp/OApWj8pFLCsrC35+ftW229nZoaCgQCNBEZF+vXKhhTQs40kJAg7movPPOQg4mIvMpyUK7RwFq0flItawYUPk5FS/9/a1a9fg7OyskaCISL94SUt7asotR8HqUXliR0BAAHbs2IHQ0FAIBAKFtpSUFPz0008YN26cxgMkIt0zp0taup5IUVNuuYqIelQeic2dOxcWFhbw8/PDN998A4FAgNjYWIwbNw4ffPABmjZtipkzZ2ozViLSEXO6pKXrUac55VYXVC5iIpEIJ06cQP/+/XHgwAFIpVL84x//wL///W+MHDkSR48ehVAo1GKoRKQr5nRJS9ejTnPKrS6o9ZxYo0aNsG7dOqxbtw55eXkoLy9Ho0aNYGHB1auITIk5XdLS9S7W5pRbXXjjZacaNWqkyTiIiPTCXHaxNlUqD6HmzJmDzp07V9vepUsXLFiwQCNBERHpSsXI6LegJjg6yImrYxgZlYvY0aNHMWzYsGrbhw4diiNHjmgkKCIiIlWoXMSys7Ph6upabburqyuys7M1EhQREZEqVL4nZmdnh4yMjGrb7969i3r16mkiJiIyIhlPSjAupS6KUnO4YC3pnMojMbFYjO3bt1dZyDIyMvDDDz9ALBar1XlycjKCg4Ph6ekJoVCI2NjYSufcvn0bn376KVxdXeHs7AyxWIxbt26p1Q8RaU/YKQmuPrXk6h6kFyqPxObOnYt//etf6NGjBz755BN4eXlBIBDg2rVriI+Ph6WlJebNm6dW50VFRfDy8sKoUaMwadKkSu0ZGRno168fgoODsX//fgiFQvzxxx+wtbVVqx8i0h51nrPiNiOkaSoXsdatW+PXX39FeHg4tm7dqtDWo0cPrFq1Cu7u7mp1HhAQgICAAADAlClTKrUvW7YMffr0wfLly+XHWrRooVYfRKRd6jxnVbE6BgCkowwTTkr4zBTVilrPiXl6euLQoUN4/PgxMjIyIJVK0apVKzRs2FDjgZWXl+PIkSOYPn06goKCcOXKFbi6uuLLL79UOkuSiHRrSy8hxvz6AIUC6xqfszKnNRlJN97oYWdHR0c4OjpqOhYFubm5KCwsxJo1azB37lwsWrQIp06dwoQJE1C/fn30799fq/0TkWrc7Opg2zvFcHevfvZyBV2vjlEdXtY0HQKJRKLyzkFlZWU4duwYMjIykJ+fr7C7MwAIBALMmjXrjQJp1qwZVq1ahdGjRwMAHjx4AE9PTwwfPlzh8mVoaCgkEgkSEhKq/Jy0tLQ36p+ItC/7uQAL/7BGfokAwjpSLH37JZrZ6H7zsnEpdXH16V8FtKNdGba9U6zzOKhmNd2mUnkklpqaik8//RRZWVmVileF2hSx1zk6OsLKygpt27ZVOP72228jMTGx2vepe19OXWlpaVrvwxgwDzL6zkNtRhSaHI2omgd3AKc6vlEXGlWUmgO8MiL8s9wKU2/Z1DoX+v73YCh0mQeVp9iHh4ejsLAQP/74I+7evYv8/PxKP3/++afGArO2tkbnzp0rjaxu374NFxcXjfVDZMxqs42IOW98+fplTEkxzDYXxk6tkdicOXMwcOBAjXVeWFiI9PR0ALKJHFlZWUhNTYWDgwNcXFzwt7/9DSEhIfDz84NYLEZSUhISExOrfJ6MyBzVZqKEOU+yeH3R35znpSgs/OsKkznlwtipPBJr3LgxrKzeeNH7Kl2+fBlisRhisRjPnz9HREQExGIxVqxYAQAYNGgQ1q5di/Xr18PPzw+bNm3Cxo0b0a9fP43GQWSsarPBojlvzvj6or9NbBR/t5lTLoydylUpLCwMu3btQlhYGOrU0cwsHn9/f0gkEqXnjB49Wj7Zg4gU1WYbEW5B8hfmwnipXMSaNm0KKysrdO/eHZ9++imaN28OS8vKf60MHTpUowESUfWqmWOlEm1uzmhsU9i5UaXxUrmIjR8/Xv7fixcvrvIcgUDAIkakQ4a6AoahxkWmR+UiduDAAW3GQURvwFAnZxhqXGR6VC5iPXv21GYcRPQGDGUFjNcZalxkelSenUhEhmdLLyG6OVmjlb0lujlZG8yEBEONi0yPWnPmc3Nz8eOPP+LKlSsoKChAeXm5QrtAIMD+/fs1GiARVc9QJyQYalxkelQuYjdv3sTAgQNRVFSE1q1b48aNG/Dw8IBEIsGDBw/QsmVLNGvWTJuxEhERKVD5cuI333yDOnXq4Ny5c9i/fz+kUikiIiJw/fp1bNmyBRKJBEuXLtVmrERERApULmJnz55FSEgIWrRoAQsL2dsqFgIePnw4hg0bhgULFmgnSiIioiqoXMRKSkrg7OwMAKhXrx4AoKCgQN7eoUMHXL58WcPhERERVU/lIta8eXPcu3cPAGBjY4MmTZrgwoUL8vbr16/D1tZW8xESERFVQ+WJHf7+/jh8+DDmz58PAPj444+xYcMGPHnyBOXl5di9ezfGjBmjtUCJiIhep3IRmz59OsRiMV68eIF69eph3rx5ePLkCf75z3/C0tISI0eOxJIlS7QZKxEZCWNbO5GMl8pFzMXFRWEzyrp162Lt2rVYu3atNuIiIiPGtRNJV1S+JzZ16lT83//9X7Xtly5dwtSpUzUSFBEZN66dSLqichGLi4vD3bt3q23PzMxEfHy8RoIiIuNmzhtukm5pbO3EP//8E3Xr1tXUxxGREePaiaQrSu+JJScn4/Tp0/LXBw4cQHp6eqXzJBIJEhMT0b59e81HSERGh2snkq4oLWJJSUlYuXIlANnivgcOHKh2XzF3d3dERERoPkIiIqJqKC1iX375JcaNGwepVAoPDw+sXr0agYGBCucIBALUr1+fDzoTEZHOKS1itra28uKUkpKCRo0aoX79+joJjIiIqCYqPyfm5OSEoqIihSKWl5eHnTt3QiKRYPDgwejSpYtWgiQiIqqKyrMTZ8yYgaCgIPnroqIi9O3bF0uXLsX69evRv39/nDt3Tq3Ok5OTERwcDE9PTwiFQsTGxlZ77rRp0yAUCrF+/Xq1+iAiItOlchE7d+4cBgwYIH+dkJCAe/fuISEhAbdu3ULbtm2xevVqtTovKiqCl5cXIiMjYWNjU+15+/btw2+//SZfRZ+IiAhQo4g9fPhQYefmX375Bd26dUPfvn3RuHFjjB49GqmpqWp1HhAQgIULF2Lw4MHyPcped+/ePcyePRtbt26FlZXKVz+JiMgMqFzEbG1tIZFIAAClpaU4c+YMevfuLW+3sbHB06dPNRpcaWkpQkNDER4ejrZt22r0s4mIyPipPLTp1KkTfvzxR4jFYvzyyy8oLCxE//795e13795F48aNNRpcREQEHBwcMH78eI1+LpG54GryZOpULmLz58/H0KFD8d5770EqleKjjz5Cp06d5O0HDx6Ej4+PxgI7ffo04uLikJSUpNb70tLSNBaDPvswBsyDjCHnYVxKXVx9Klu3MB1lGPPrA2x7p1grfRlyHnSJeZDRVB7c3d2VtqtcxN555x1cvHgR58+fh52dHfz9/eVtEokEoaGh6NGjx5tH+pqkpCTk5OQoXEYsKyvDokWLEB0djevXr1f5vpq+cG2lpaVpvQ9jwDzIGHoeilJzAPy1gnyhwBru7q4a78fQ86ArzIOMLvOg1kwJR0dHfPjhh5WOC4VCTJ48WWNBAUBoaCgGDx6scCwoKAhBQUEYO3asRvsiMlWN6loi/ZUixtXkydTodbpfYWGhfEHh8vJyZGVlITU1FQ4ODnBxcYGTk+IColZWVhCJRPxLh0hFW3oJMeGk4j0xIlNSbRFzcHCAhYUFHjx4AGtrazg4OEAgECj9MIFAgMePH6vc+eXLlxXWYoyIiEBERARGjRqF6OholT+HiKrG1eTJ1FVbxGbNmgWBQCB/NqvitSb5+/vLp+2r4urVqxrtn4iIjFu1RWzOnDlKXxMREembxnZ2JiIi0jWVJnYUFxdj9+7dOH78OO7evYvCwkI0aNAArVq1Qp8+fTBixAhYW1trO1YiIiIFNRaxa9eu4ZNPPsH9+/chlUphb2+PBg0aIDc3FykpKdi7dy/WrFmD+Ph4Lg1FREQ6pfRyYmFhIUaNGoXc3FwsWLAA165dQ2ZmpsL/zp8/Hzk5OQgODkZRUZGu4iYiIlJexGJjY5GVlYXdu3djxowZaNq0qUJ706ZN8dVXXyE+Ph6ZmZmIi4vTarBkfDKelCDgYC46/5yDgIO5yHxaou+QiMiEKC1iR48eRZ8+fRSWmKpKr1698N577+HIkSMaDY6MX9gpCS7kvkT6kzJcyH2JCScl+g6JiEyI0iJ2/fp19OzZU6UPEovF1a5nSOYrr7hM6WsiotpQOrEjPz9f5e1VnJyckJ+fr5GgyHjUtNVHTWv3casQIqoNpSOx4uJi1Kmj2i8UKysrvHz5UiNBkfGo6XLhll5CdHOyRit7S3Rzsq60dh8vNxJRbdQ4xT4jIwOXLl2q8YPu3r2rkYDIuNR0ubCmtft4uZGIaqPGIlaxKG9NpFKpxtdWJMNX260+uFUIEdWG0iIWFRWlqzjISNV2qw9uFUJEtaG0iH3yySe6ioOMVG23+niT93MyCBFV4ALAZHQ4GYSIKrCIkdHhZBAiqsAiRkbn9ckfnAxCZL5YxMjo1PTsGRGZD5X2EyOqjj4mWdR2MgkRmQ6VR2LJycnIy8urtv3x48dITk7WSFBkPDjJgoj0SeUiFhgYiOPHj1fbfvLkSQQGBmokKDIenGRBRPqkchGTSqVK21++fAkLC95iMzecZEFE+qT0ntiTJ09QUFAgf/3nn3/i/v37lc6TSCT4+eef4ezsrPkIyaBxxQ0i0ielRWzDhg1YtWoVAEAgEGDOnDmYM2dOledKpVIsWLBArc6Tk5Oxfv16pKSk4MGDB4iKisLo0aMBACUlJVi2bBn+9a9/ISMjA3Z2dvD398eiRYvg4uKiVj+kPZxkQUT6pLSI9e7dG/Xq1YNUKsWSJUswbNgwdOjQQeEcgUCA+vXro1OnTvD29lar86KiInh5eWHUqFGYNGmSQtuzZ8+QkpKC8PBwdOjQAU+ePMH8+fMxfPhwJCcnw8qKEyuJiMyd0krg6+sLX19fALK9xQIDA9GuXTuNdR4QEICAgAAAwJQpUxTa3nrrLezdu1fh2HfffQdfX1/cunVLo3EQEZFxUnk4M3v2bG3GoZKnT58CAIRCoX4DISIigyCQSCTKpx3+18qVK2v+MIEAs2bNeqNAmjVrhlWrVsnvib3u5cuXCAwMhIODA3bt2lXt56Slpb1R/0REZHjc3d2Vtqs8EouMjKy2TSAQyDfFfNMipkxpaSnCwsJQUFCA+Ph4pefW9IVrKy0tTet9GAPmQYZ5kGEeZJgHGV3mQeUilp+fX+lYeXk57t27h02bNuH8+fNISEjQaHCArICNHz8e169fx8GDB9GwYUON90FERMapVk8nW1hYoEWLFoiIiICbm5vG75uVlJQgJCQE165dw4EDByASiTT6+UREZNw0Nk/d398fixcvVus9hYWFSE9PByAb1WVlZSE1NRUODg5wdnbG2LFjcfnyZcTHx0MgEODhw4cAAHt7e9jY2GgqdCIiMlIaWycqLS2txqWpXnf58mWIxWKIxWI8f/4cEREREIvFWLFiBbKzs3H48GE8ePAAvXv3Rtu2beU/iYmJmgqbiIiMmMojsepWqC8oKEBSUhK2bNmCIUOGqNW5v78/JBJJte3K2oiIiFQuYoMGDYJAIKh0XCqVwtLSEkFBQSpNwyciItIUlYvY/v37KxUxgUAAoVAIV1dX2NnZaTw4IiIiZVQuYv7+/tqMg4iISG01TuzYuXMnfHx8IBKJ4Onpiblz5+Lly5e6iI2IiEgppUXsH//4B6ZNm4bs7Gy0a9cO5eXl2LhxIxYuXKir+IiIiKqltIht3rwZrVq1wm+//YZjx47h2rVrGDp0KGJiYvDs2TNdxUhERFQlpUXsxo0b+Pzzz9G4cWMAgJWVFWbMmIHi4mJkZGToIj4iIqJqKS1iRUVFaNKkicKxpk2bAgDy8vK0FxUREZEKapzYUdWzYURERIagxin269atw+7du+WvS0pKAACLFy+utKK8QCDAnj17NBwiERFR1ZQWsebNm6OgoAAFBQUKx11cXJCbm4vc3FyF4xy1ERGRLiktYlevXtVVHERERGrT2Cr2REREusYiRkRERotFjIiIjBaLGBERGS0WMSIiMlosYkREZLRULmKBgYE4efJkte2nTp1CYGCgRoIiIiJShcpF7PTp03j06FG17Xl5eUhOTtZIUERERKrQ2OXE7Oxs2NraaurjiIiIaqR0xY5Dhw7h8OHD8tcxMTE4ceJEpfMkEglOnjyJLl26aDxAUi7jSQnCTkmQV1yGRnUtsaWXEG52dfQdFhGRTigtYjdu3MDPP/8MQLYu4sWLF3Hp0iWFcwQCAerXrw9fX19ERkaq1XlycjLWr1+PlJQUPHjwAFFRURg9erS8XSqVIjIyEjt27IBEIkGXLl2wevVqeHp6qtWPKQs7JcGF3JcAgHSUYcJJCY4OctJzVEREuqH0cmJ4eDhycnKQk5MDqVSKqKgo+euKnwcPHuDOnTtISEhAmzZt1Oq8qKgIXl5eiIyMhI2NTaX2devWISoqCitXrsSxY8fg5OSEoUOH4unTp+p9SxOWV1ym9DURkSmrcSuWCvn5+RrvPCAgAAEBAQCAKVOmKLRJpVJER0dj+vTpGDx4MAAgOjoa7u7uSEhIQEhIiMbjMUaN6loiHWUKr4mIzIXKRazC3bt3cfToUdy7dw8A4OrqioCAALRs2VKjgWVmZuLhw4fo06eP/JiNjQ38/Pxw/vx5FrH/2tJLiAknFe+JERGZC7WK2Lx587Bx40aUl5crHJ87dy4mTZqE5cuXayywhw8fAgCcnBTv7zg5OeHBgwfVvi8tLU1jMeizD3VEtf3rv1/mSJCWo5t+DS0P+sI8yDAPMsyDjKby4O7urrRd5SIWFRWFDRs2YNCgQfjb3/6Gtm1lvzlv3bqF9evXIzo6Gs2aNat0WbC2Xt9oUyqVKt18s6YvXFtpaWla78MYMA8yzIMM8yDDPMjoMg8qPye2c+dOBAQE4Mcff0TXrl1hb28Pe3t7dO3aFTt37sT777+PmJgYjQUmEokAoNID1nl5eZVGZ0REZJ5ULmIZGRnySRhVCQgIQGZmpkaCAgA3NzeIRCIcP35cfuzFixc4e/YsfHx8NNYPEREZL5UvJzo4OCi9xnn79m04ODio1XlhYSHS09MBAOXl5cjKykJqaiocHBzg4uKCyZMn49tvv4W7uzvatGmD1atXw9bWFsOHD1erHyIiMk0qj8Q+/PBDbNu2DbGxsZBKpfLjUqkUcXFx2L59OwYOHKhW55cvX4ZYLIZYLMbz588REREBsViMFStWAACmTZuGKVOmYObMmXjvvfeQk5ODxMRE2NnZqdUPERGZJoFEIpHWfJpsaalBgwbh+vXrcHR0ROvWrQEA6enpyM3NRfv27XHgwAEIhUJtxqt3vHErwzzIMA8yzIMM8yCjyzyofDlRKBTi2LFjiImJUXhOrGPHjujXrx8+++wz1K1bV2uBEhERvU6t58Ssra0RFhaGsLAwbcVDRESkMu7sTERERkutkdiJEyewY8cOZGRkID8/X2GCByB7MPnKlSuajI+IiKhaKhex6OhozJs3D40aNYK3tze3QyEiIr1Ta9mpHj164Oeff4a1tbU2YyIiIlKJyvfEHj9+jGHDhrGAERGRwVC5iL377rvyafVERESGQOUitnz5csTFxeHUqVPajIeIiEhl1d4T+/jjjysds7e3x5AhQ9C6dWu4uLjA0lJxF2GBQIA9e/ZoPkoiIqIqVFvEbt68WeW+Xc2bN0dxcTFu375dqU3ZPl9ERESaVm0Ru3r1qi7jICIiUhtX7CAiIqOl1oodr0pKSsKePXuQk5ODt99+G5MmTYKLi4smYyMiIlJK6UgsMjISTk5OePjwocLx2NhYDB48GD/99BP+/e9/Y8OGDejTpw+n4BMRkU4pLWJJSUno06cPRCKR/FhxcTHmzJkDe3t77Nu3D1lZWdi+fTsKCwuxZs0arQdMRERUQWkRS09Ph7e3t8KxkydP4unTp/jiiy8gFotha2uLoUOHYsSIEThx4oQ2YyUiIlKgtIjl5+ejSZMmCseSkpIgEAjQr18/hePvvvsucnJyNB8hERFRNZQWscaNG+M///mPwrGzZ8+iQYMGaN++veIHWVhwXUUiItIppUWsc+fOiIuLg0QiAQD8/vvvuHz5MsRicaUHm2/duoVmzZppLVAiIqLXKZ1iP3PmTPTp0wedO3eGh4cHfv/9dwgEAkybNk3hPKlUioMHD6JPnz5aDZaIiOhVSkdi7dq1w759++Dt7Y28vDx069YNiYmJ6Nq1q8J5SUlJaNCgAT766COtBktERPSqGh929vX1rXFRX7FYjDNnzmgsqAplZWWIiIjAnj178PDhQ4hEIowYMQKzZ8+GldUbP6dNREQmwqArwdq1a7F161ZER0fDy8sL165dw+TJk2FtbY1Zs2bpOzwiItIzgy5iFy5cQP/+/TFgwAAAgJubGwYMGIBLly7pOTIiIjIEBr0AsK+vL06fPo0//vgDgGx7mKSkJHzwwQd6joyIiAyBQCKRSPUdRHWkUimWLVuGNWvWwNLSEqWlpQgPD8f8+fOrfU9aWpoOIyQiIm1yd3dX2m7QlxMTExOxa9cubN26FR4eHrh69Spmz54NV1dXfPbZZ1W+p6YvXFtpaWla78MYMA8yzIMM8yDDPMjoMg8GXcQWLlyIL774AkFBQQBkU/7v37+P7777rtoiRkRE5sOg74k9e/YMlpaWCscsLS1RXl6up4iIiMiQGPRIrH///li7di3c3Nzg4eGB1NRUREVFITg4WN+hERGRATDoIrZq1SosX74c//M//4O8vDyIRCKMHTuWz4gREREAAy9idnZ2iIyMRGRkpL5DISIiA2TQ98SIiIiUYREjIiKjxSJGRERGi0WMiIiMFosYEREZLRYxIiIyWixiRERktFjEiIjIaLGIERGR0WIRIyIio8UiRkRERotFjIiIjBaLGBERGS0WMSIiMlosYkREZLRYxIiIyGixiBERkdFiESMiIqNlpe8ADEXGkxKEnZIgr7gMjepaYksvIdzs6ug7LCIiUoIjsf8KOyXBhdyXSH9Shgu5LzHhpETfIRERUQ1YxP4rr7hM6WsiIjI8LGL/1aiupdLXRERkeAy+iOXk5GDSpElo3bo1RCIRfHx8cPr0aY33s6WXEN2crNHK3hLdnKyxpZdQ430QEZFmGfTEDolEgn79+sHX1xd79uyBo6MjMjMz4eTkpPG+3Ozq4OggzX8uERFpj0EXse+//x5NmjTBpk2b5MdatGihv4CIiMigGPTlxEOHDqFLly4ICQlBmzZt0LNnT2zevBlSqVTfoRERkQEQSCQSg60IIpEIADBlyhQMGTIEV69exddff41FixYhLCysyvekpaXpMkQiItIid3d3pe0GXcScnJzQqVMnHD16VH5syZIlOHjwIC5cuKCXmNLS0mpMqjlgHmSYBxnmQYZ5kNFlHgz6cqJIJELbtm0Vjr399tvIysrSU0RERGRIDLqI+fr64vbt2wrHbt++DRcXFz1FVPPQ1lwwDzLMgwzzIMM8yOgyDwZdxKZMmYKLFy9i9erVSE9Px969e7F582aEhobqOzQiIjIABn1PDAB+/fVXLFmyBLdv30bz5s0xYcIETJw4EQKBQN+hERGRnhl8ESMiIqqOQV9OJCIiUoZFjIiIjBaLWA2+/fZbCIVCzJw5U35MKpUiIiICHh4eaNKkCQYOHIgbN27oMUrtqGnxZXPIQ1lZGZYtW4aOHTtCJBKhY8eOWLZsGUpLS+XnmGIekpOTERwcDE9PTwiFQsTGxiq0q/Kdi4uLMXPmTLRq1QpNmzZFcHAwsrOzdfk1ak1ZHkpKSrBo0SL4+fmhadOmaNu2LUJDQ3H//n2FzzD1PLxu2rRpEAqFWL9+vcJxbeWBRUyJixcvYseOHWjXrp3C8XXr1iEqKgorV67EsWPH4OTkhKFDh+Lp06d6ilTzKhZflkql2LNnD86fP49Vq1YpLL5sDnlYu3Yttm7dipUrV+LChQuIjIzEli1bsGbNGvk5ppiHoqIieHl5ITIyEjY2NpXaVfnOc+bMwYEDB7Bt2zYcPnwYT58+xciRI1FWZjx79SnLw7Nnz5CSkoLw8HCcPHkScXFxyM7OxvDhwxX+yDH1PLxq3759+O233+Ds7FypTVt54MSOahQUFKBXr15Yt24dVq1aBS8vL/z973+HVCqFh4cHJkyYgPDwcADA8+fP4e7ujqVLlyIkJETPkWvGkiVLkJycjF9//bXKdnPJw8iRI+Hg4ICNGzfKj02aNAn5+fnYvXu3WeShWbNmWLVqFUaPHg1Atf/vCwoK0KZNG0RFRWHEiBEAgKysLHTo0AEJCQno27ev3r7Pm3o9D1W5efMmfH19kZycjHbt2plVHu7du4d+/fph7969GD58OMLCwvDll18CgFbzwJFYNaZPn47BgwejV69eCsczMzPx8OFD9OnTR37MxsYGfn5+OH/+vK7D1JqaFl82lzz4+vri9OnT+OOPPwDIfkklJSXhgw8+AGA+eXiVKt/5ypUrKCkpUTinefPmaNu2rcnmBYB8JCoUCgGYTx5KS0sRGhqK8PDwSqssAdrNg0FvxaIvO3bsQHp6usIWMBUePnwIAJX2NHNycsKDBw90Ep8uZGRkYNu2bZgyZQqmT58uX3wZAMLCwswmD9OnT0dhYSF8fHxgaWmJ0tJShIeHyx+4N5c8vEqV7/zo0SNYWlrC0dGx0jmPHj3STaA69vLlS8yfPx/9+/dHs2bNAJhPHiIiIuDg4IDx48dX2a7NPLCIvSYtLQ1LlizBL7/8Amtr62rPe/1ha6lUalIPYJeXl6NTp05YtGgRAOCdd95Beno6tm7dqrCDgKnnITExEbt27cLWrVvh4eGBq1evYvbs2XB1dcVnn30mP8/U81CVN/nOppqX0tJShIWFoaCgAPHx8TWeb0p5OH36NOLi4pCUlKT2ezWRB15OfM2FCxfw+PFjdO/eHY6OjnB0dERycjK2bt0KR0dHNGzYEAAq/fWQl5enlR2n9aWmxZcrtskx9TwsXLgQX3zxBYKCgtCuXTsEBwdj6tSp+O677wCYTx5epcp3bty4McrKyvD48eNqzzEVpaWlGD9+PK5du4Z9+/bJf0cA5pGHpKQk5OTkoG3btvLfmffv38eiRYvg5eUFQLt5YBF7zcCBA3HmzBkkJSXJfzp16oSgoCAkJSWhTZs2EIlEOH78uPw9L168wNmzZ+Hj46PHyDWrpsWX3dzczCIPz549g6WlpcIxS0tLlJeXAzCfPLxKle/87rvvok6dOgrnZGdn49atWyaVl5KSEoSEhODatWs4cOCAvMBXMIc8hIaGIjk5WeF3prOzM6ZMmYJ9+/YB0G4eeDnxNUKhUH5TtkL9+vXh4OAg/6ti8uTJ+Pbbb+Hu7o42bdpg9erVsLW1xfDhw/UQsXZMmTIFAQEBWL16NYYNG4bU1FRs3rwZCxYsACC7lGQOeejfvz/Wrl0LNzc3eHh4IDU1FVFRUQgODgZgunkoLCxEeno6ANml5aysLKSmpsLBwQEuLi41fue33noLY8aMwcKFC+Hk5AQHBwfMmzcP7dq1Q+/evfX4zdSjLA/Ozs4YO3YsLl++jPj4eAgEAvn9Qnt7e9jY2JhFHlxcXCqNpqysrCASieSr2WszD5xir4KBAwfKp9gDsuu4kZGRiImJgUQiQZcuXbB69Wp5kTMVNS2+bA55ePr0KZYvX46DBw8iLy8PIpEIQUFBmDVrFurVqwfANPOQlJSEwMDASsdHjRqF6Oholb7zixcvsGDBAiQkJODFixcQi8X49ttv0bx5c11+lVpRlofZs2fjnXfeqfJ9UVFR8inopp6H6OjoSsc7dOigMMUe0F4eWMSIiMho8Z4YEREZLRYxIiIyWixiRERktFjEiIjIaLGIERGR0WIRIyIio8UiRqQBsbGxEAqFyMzMVPu9mZmZEAqF8qWsdK2i/9c3Orxy5QoGDBiA5s2bQygUvtHaeOqYPHkyOnTooNU+yPSwiJHJ2r9/P4RCIRISEiq1BQYGKm1zc3OTbztjSM6ePYuIiAhIJBKVzq8orhcvXlSrn7KyMoSEhOA///kPlixZgk2bNqFt27bYvXs3NmzY8AaRE2kHixiZrO7duwOQ/eJ/VWlpKS5dugQrK6tq23x9fdVaXTs4OBg5OTlwdXWtfeBKnDt3DitXrkRBQYHGPtPV1RU5OTnypbQA2YaFd+/excSJEzFu3DiMHDkSjRs3xp49e6pcoYFIX7h2IpksJycntG7dulKhSklJwbNnzzBixIhq23x9fdXqy9LSstJCwcZCIBDIl9CqkJeXB0C25h2RIeNIjExa9+7dcfPmTYXLb+fOnYOzszNGjhxZZVvF+yocP34cgwYNQvPmzdG0aVMMGjSo0m601d0Ti4mJQadOnSASidCzZ08cOXJE6b2f+Ph4dO3aFY0bN4afnx9OnDghb4uIiMDixYsByPZ3q1isurb3ql6/JzZ58mT5dvFTp06FUChEhw4dMHDgQPzv//4v7t+/L+/71cWypVIpNm/eDD8/P4hEIrRs2RITJkxAdnZ2pT5/+ukndOnSBSKRCD169MAvv/xSq+9A5osjMTJpvr6++Omnn3DhwgUEBAQAkBUqHx8fdO3aFQAqtdWrVw+dOnUCACQkJCAsLAz+/v6YN28eysvLERsbi48++giHDh2Ct7d3tX3HxMRg+vTp6Nq1K8LCwpCXl4eJEyfKd/193b59+/D48WOEhISgXr16iI6OxqeffoqrV6/CwcEBgYGBSEtLQ2JiIlasWCHfJbeq7eBrIyQkBG5uboiMjMTnn3+O7t27w9bWFra2tpBIJMjJycGKFSsqve+rr77Czp07MXLkSISGhuLhw4fYvHkzzp8/j1OnTskLXlxcHL744gt07twZoaGhyM3NxcSJE41qQVwyHCxiZNIqRlTnzp2TF6rz589jxowZsLe3h4eHR6W2Tp06oW7duigqKkJ4eDhGjhypcB8oJCQEvr6+WLJkCfbv319lvyUlJVi6dCnat2+PQ4cOyXcJF4vFGDx4sHxftlfdvXsXly5dQqNGjQAAPXv2hFgsRkJCAiZMmID27dujQ4cOSExMxMCBA+Hm5qa5RL2iW7duEAgEiIyMRNeuXTFy5Eh5W5MmTfDkyROFY4Asbz/88IPC6u2AbJJM7969sXnzZsyaNQulpaX45ptv4OHhgcOHD8svY/bs2RPDhg2rMi9EyvByIpm01q1bQyQSye993blzB48ePZLf8/L19a3U5ufnB0B2GVEikWDEiBF4/Pix/Of58+fo3bs3zp49i5KSkir7/e233/D48WN8/vnn8gIGAL169YKnp2eV7xkyZIi8gAFAx44dYW9vj4yMjFrnQdv++c9/okGDBggICFDIlbOzM1q3bo1Tp04BkOXl0aNH8tFmhT59+sDDw0Nf4ZMR40iMTJ6Pjw+OHj2Kly9f4ty5c6hfv778npSPjw/i4uLkbQDkBe7OnTsAgKFDh1b72QUFBQqFp8L9+/cByIro61q3bo2UlJRKx6sahbz11lvIz8+v6Svq3Z07d1BYWCjfBPF1FTM9K/JS1Xlt2rSpMi9EyrCIkcnz9fXF/v37cfnyZZw7dw5dunSBlZXsn76Pjw9evHghb7OwsEC3bt0AyHawBYANGzagadOmVX62vb292vFU9/xZdbMbDfF5tdeVl5ejYcOG2L59e5Xt9evXB/DXd6nq8QVj+J5keFjEyORVXB48d+4czp07hyFDhsjbWrRogSZNmsjb2rVrJ59W3rJlSwBAo0aN1N5CvWJUdefOHbz33nsKbRXbvL8JdZ5d04bq+m/ZsiWOHz+OLl26wM7Ortr3VzxH98cff1TKS8XIl0gdvCdGJq9Dhw5o0KABDh06hLS0tErPgPn4+FTZ1rdvX7z11ltYvXo1iouLK31uxbNUVenUqRMcHR0RExODly9fyo+fPHkSN27ceOPvUjGiUXXFDk2rX79+lQ9aDxs2DOXl5YiMjKzUJpVK8fjxYwCyvDg5OSEmJgYvXryQn3Ps2DHcvHlTe4GTyeJIjEyepaUlvL29ceLECVhYWFSaFu/j44O5c+cC+GvUBgB2dnZYt24dxo8fj549e+Ljjz+GSCRCdnY2kpKSYGtrW+WyVQBgbW2NefPm4auvvsLAgQMRFBSEvLw8bNmyBV5eXigsLHyj71Ix9X/p0qUICgqCtbU1xGIxnJyclL4vLi5O4ZmzCmPHjlW7//379+Prr7+Gt7c3LCwsEBQUBD8/P0ycOBFRUVH4/fff8f7776N+/frIzMzEwYMHMWbMGMyYMQN16tTBwoUL8eWXX+LDDz/Exx9/LM+Lp6fnG+eFzBeLGJmF7t2748SJE/D09Ky0CsWro6/XR2lDhgyBs7Mz1qxZgw0bNuD58+cQiUTw9vbGZ599prTPcePGAQC+//57LFq0CO7u7ti0aRPi4uLeeNTRtWtXzJ8/HzExMZg6dSrKy8tx4MCBGovYDz/8UOXxfv36qbUqR1hYGG7evIk9e/Zg8+bNkEqlCAoKAgCsXLkS7777LrZt24aIiAhYWFigadOm6Nu3LwYNGiT/jDFjxkAqlWLt2rVYtGgR2rRpg02bNmH//v04ffq0yrEQAYBAIpHwbiqRDvXo0QNOTk7Yu3evvkMhMnq8J0akJS9evKg04+7kyZO4du0axGKxnqIiMi0ciRFpSVJSEr7++mt89NFHaNKkCW7cuIGYmBg4OjrizJkzCusOEtGb4T0xIi1xdXWFm5sbfvjhB/z555+wt7fHoEGDsHDhQhYwIg3hSIyIiIwW74kREZHRYhEjIiKjxSJGRERGi0WMiIiMFosYEREZLRYxIiIyWv8PxdIN2BumRe4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "shotput.plot.scatter('Weight Lifted', 'Shot Put Distance');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's not a football shaped scatter plot. In fact, it seems to have a slight non-linear component. But if we insist on using a straight line to make our predictions, there is still one best straight line among all straight lines.\n", "\n", "Our formulas for the slope and intercept of the regression line, derived for football shaped scatter plots, give the following values." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.09834382159781994" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slope(shotput, 'Weight Lifted', 'Shot Put Distance')" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5.959629098373956" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "intercept(shotput, 'Weight Lifted', 'Shot Put Distance')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Does it still make sense to use these formulas even though the scatter plot isn't football shaped? We can answer this by finding the slope and intercept of the line that minimizes the mse.\n", "\n", "We will define the function `shotput_linear_mse` to take an arbirtary slope and intercept as arguments and return the corresponding mse. Then `minimize` applied to `shotput_linear_mse` will return the best slope and intercept." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def shotput_linear_mse(any_slope, any_intercept):\n", " x = shotput['Weight Lifted']\n", " y = shotput['Shot Put Distance']\n", " fitted = any_slope*x + any_intercept\n", " return np.mean((y - fitted) ** 2)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.09834382, 5.9596291 ])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "minimize(shotput_linear_mse)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These values are the same as those we got by using our formulas. To summarize:\n", "\n", "**No matter what the shape of the scatter plot, there is a unique line that minimizes the mean squared error of estimation. It is called the regression line, and its slope and intercept are given by**\n", "\n", "$$\n", "\\mathbf{\\mbox{slope of the regression line}} ~=~ r \\cdot\n", "\\frac{\\mbox{SD of }y}{\\mbox{SD of }x}\n", "$$\n", "\n", "$$\n", "\\mathbf{\\mbox{intercept of the regression line}} ~=~\n", "\\mbox{average of }y ~-~ \\mbox{slope} \\cdot \\mbox{average of }x\n", "$$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAApMAAAGTCAYAAAB0/plgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABExklEQVR4nO3de1iUZf7H8c8wZCiCoBGMp5GU8FyslocSzfbnYbNVU9PWDmsqHqpfJ8tcJ9mMTXSz9GdmipXtLraVa6gdtq7SXCS1pSytTcOkWTXwlCjHBGZ+f5BTI4Mww8AM8H5dl9el9/PM83y5kfr43M9934a8vDy7AAAAAA8E+LoAAAAANFyESQAAAHiMMAkAAACPESYBAADgMcIkAAAAPEaYBAAAgMcIkwAAAPBYkwmTWVlZvi6hwaHP3EefuY8+cx995j76DKg7TSZMAgAAwPsIkwAAAPAYYRIAAAAeI0wCAADAY4RJAAAAeIwwCQAAAI9VGyafeeYZ3XDDDerQoYM6d+6siRMn6j//+Y/TOXa7XYsWLVLXrl0VFRWlm266SV9//XWdFQ0AAAD/UG2Y3LFjh6ZOnar33ntPmzdvVmBgoMaMGaPTp087zlm+fLlWrlypxYsXa+vWrYqIiNDYsWOVn59fp8UDAADAtwKrO2Hjxo1Of169erU6duyoXbt2aeTIkbLb7Vq1apUeeOABjR49WpK0atUqxcTEaMOGDZoyZUrdVA4AAACfqzZMXqigoEA2m01hYWGSJKvVqmPHjmno0KGOc5o3b66BAwdq9+7dVYZJX+xGwA4I7qPP3EefuY8+cx995r666rOYmJg6uS7QULgdJh977DH16tVL1157rSTp2LFjkqSIiAin8yIiIpSTk1Plder7hy8rK4sfeDfRZ+6jz9xHn7mPPnPf+T4zlFkVlJ+kgPIc2YwmlYRYZA80+7o8oEFzazb3H/7wB+3atUt//etfZTQanY4ZDAanP9vt9kptAAD4iqHMquAfxqhZyRsKLN2hZiVvKPiHMTKUWX1dGqqRmpqqdu3a+boMSdKKFSvUq1evWl8nPT1dYWFhOnXqlMvjVqtVYWFh2rNnT5XXqMk59aHGYXLevHn6xz/+oc2bN6tTp06O9sjISEnS8ePHnc4/efJkpaeVAAD4SlB+kozl2U5txvJsBeUn+aiipuP48eOaO3eurr76al1++eXq1q2bxo8fr/fff9/XpcELahQm586dqw0bNmjz5s268sornY6ZzWZFRkZq27ZtjraSkhLt3LlT/fr18261AAB4KKDc9atXAeW59VxJ02K1WjV48GBt3bpViYmJysjIUFpamoYNG6aHHnqo3uo4d+5cvd2rOv5UizdUGybnzJmj9evXa+3atQoLC9OxY8d07NgxFRQUSKoY3p41a5aWLVumzZs36z//+Y9mz56t4OBgjR8/vs6/AAAAasJmNFXRHlXPlfiO1XpW06dv1ahRb2n69K2yWs/W+T3nzJkju92ubdu2aezYsYqJiVFsbKwSEhK0Y8cOSdLhw4c1efJktW/fXu3bt9ftt9+uo0ePXvS6L7/8suLi4hQREaG4uDi98sorTsfDwsKUkpKi22+/XW3bttXChQs9qn/58uW68sor1a5dO82YMcORf8777LPPNHbsWF1xxRXq0KGDRowYoU8++cTtWn788UdNnjxZ8fHxOnHihKP94MGDGjFihCIjI3XNNddo69atVdbqaujc1VD4/v37deutt6p9+/bq0qWLpk6d6pgD44lqw+TatWuVn5+v0aNHKzY21vFrxYoVjnPuv/9+zZ49W4888ohuuOEG5ebmauPGjQoJCfG4MAAAvKkkxKJyY7RTW7kxWiUhFh9VVL+s1rMaM+ZdvfHGt9qxI0dvvPGtxox5t04D5enTp/XBBx9o+vTpatmyZaXjYWFhstvtmjx5sk6cOKHNmzdry5Ytys3N1eTJk2W3211ed8uWLXrkkUc0a9Ys7dy5UzNnztTDDz+sd9991+m8xYsXa9iwYfr44481bdo0SVK7du0u+uuXD8LefPNNJSUlad68edq+fbtiYmL0/PPPO90jPz9fEydO1LvvvqsPP/xQvXr10oQJEyq9C+mqlvPOnj2rcePG6fTp03rrrbecXhNMTEzUjBkzlJ6eriFDhuh3v/udvv/++xr0vmu5ubn6zW9+o27duunDDz9UWlqaCgoKdNttt8lms3l0zWpnc+fl5VV7EYPBoHnz5mnevHkeFQEAQF2zB5pV2Drtp9ncubIZo5rUbO6kpExlZzsHx+zss0pKylRKytAqPlU7hw4dkt1ur/SK3C999NFH+vLLL7Vnzx6ZzRXfi7Vr1youLk7bt2/XkCFDKn3mueee08SJE5WQkCBJ6tKliz7//HMtX75cI0eOdJw3duxY3XnnnU6fTU9Pv2jNQUFBjt+vWrVKt912m2OZwzlz5ig9PV2HDh1ynDN48GCnzy9ZskSbN2/WBx98oIkTJ1ZZy+HDhyVVzDFJSEiQyWTSunXrnO4vSXfffbfGjh0rSY7NYV566SVZLJ79I+jFF19Uz5499cQTTzjaVq9erU6dOmnPnj3q06eP29d0e2kgAAAaKnugWcXhKb4uwydycopctufmum73hqqeLP7SgQMHZDKZHEFSkjp16iSTyaT9+/e7DJMHDhzQ5MmTndoGDBhQ6clkXFxcpc9eccUVNay+4j533HGHU9s111zjFCZPnDihP/3pT0pPT9eJEydUXl6u4uJiHTlypNpaJOmWW25R79699de//lWBgZVj2TXXXOP4fUBAgPr06aP9+/fX+Gu40BdffKGPP/7Y5ez47OxswiQAAHDNZGrhsj0qynW7N3Tu3FkGg0HffPNNledcbCnBiy0x6OrYhW3BwcGVzqluiaEBAwZow4YNFz3nl2bNmqXjx4/rqaeeUseOHXXppZfqt7/9baVJNq5qkaThw4frzTff1FdffaWrrrqqxvd1JSCg4u3FX4b4srIyp3NsNpuGDRumpKTKqxh4ugoPYRIAgCbAYumrzMwTTkPd0dGhslj61tk9w8PDdeONNyolJUUzZsyo9N5kXl6eunbtqu+//15Wq9XxdPK7775TTk6Ounbt6vK6sbGx2rVrl9NTw507d1Z5/i+5M8wdGxurzMxMp/tkZmY6nb9r1y4lJydr+PDhkiqWQXJnMsv8+fMVHh6uMWPGaNOmTerdu7fT8czMTMdQut1u12effebYvvpCl112maSK9yLP/37fvn1O51x11VV688031aFDB11yySU1rvNiCJMAADQBZnOo0tJGKikpU7m5RYqKaiGLpa/M5tA6ve/TTz+t4cOH64YbbtD8+fPVo0cP2e12paen69lnn9W+ffvUs2dPJSQkaPHixbLb7Xr00Ud11VVXKT4+3uU177vvPv3+97/X1VdfraFDh+qDDz7QG2+8ob/+9a/V1uPOMPfMmTM1c+ZM/epXv9L111+vTZs26dNPP3VsKS1VPH19/fXX1bdvXxUVFWnBggVq1qxZje8hSY8//rjsdrsjUP5yUfSXXnpJXbp0Uffu3bV27VodPnxYd999d5VfW/v27ZWcnKw//vGP+u9//6s///nPTudMmzZNr7zyiqZMmaIHHnhAl112mb777jvHZCNPJk8TJgEAaCLM5tA6m2xTlU6dOmn79u1aunSpEhMTlZOTo9atW6tnz5569tlnZTAYlJqaqrlz52rUqFGSKia1LFmypMph7lGjRmnJkiVasWKF5s2bpw4dOmjp0qVOk2+84ZZbbtF3332nJ598UsXFxRo5cqRmz56t9evXO8557rnn9MADD2jIkCGKiorSY489VuWuNhezYMEC2e12jR49Wps2bVJoaEXIT0xM1MqVK/XFF1+oQ4cO+tvf/lblUP0ll1yiF198UQ8//LCuv/569erVSwsWLHCaCGQymfTee+/piSee0Lhx4/Tjjz+qffv2uuGGG3TppZe6XbckGfLy8qp/O7YRYC9b99Fn7qPP3EefuY8+cx99BtQdt/bmBgAAAH6JMAkAAACP8c4kAPgpq7ViQemcnCKZTPUzWQIA3EWYBAA/dH7ru18u45KZeUJpad6dYAAAtcUwNwD4oYttfQcA/oQwCQB+yBdb3wGAJwiTAOCHfLH1HQB4gjAJAH7IYumr6GjnyTZ1vfUdAHiCMAkAfuj81ncTJnTWoEEmTZjQWWlpI5nNDcDvECYBwE+d3/puy5ZRSkkZSpBEk5WamlrlFoINQVhYmDZt2lTl8VOnTiksLEzp6en1WJX3ECYBAECdOn78uObOnaurr75al19+ubp166bx48fr/fff93Vp8ALWmQQAAHXGarVqxIgRatmypRITE9WzZ0/ZbDZt375dDz30kL788st6qePcuXNq1qxZvdyrqeHJJAAATYShzKrmp6cr+OQoNT89XYYya53fc86cObLb7dq2bZvGjh2rmJgYxcbGKiEhQTt27JAkHT58WJMnT1b79u3Vvn173X777Tp69OhFr/vyyy8rLi5OERERiouL0yuvvOJ0PCwsTCkpKbr99tvVtm1bLVy40KP6n3nmGcXExKhdu3aaMWOGkpOT1atXL8dxm82mJUuWqEePHrr88ss1cOBAvf322xe95meffabBgwcrMjJSgwYNUmZmw14/lieTAAA0AYYyq4J/GCNjeXZFQ6lkLM1UYes02QPNdXLP06dP64MPPpDFYlHLli0rHQ8LC5PdbtfkyZMVFBSkzZs3y2Aw6JFHHtHkyZO1bds2GQyGSp/bsmWLHnnkET311FMaOnSoPvzwQz388MO6/PLLNXLkz7tELV68WAsWLFBSUpKjrbp3LwcMGKANGzZIkv7xj39o8eLF+vOf/6yBAwdq8+bNWrZsmVq1auU4f9WqVVqxYoWeeeYZxcXF6bXXXtMdd9yhjz76SL179650/cLCQt1666267rrrtGrVKuXk5GjevHnVd6YfI0wCANAEBOUn/Rwkf2Isz1ZQfpKKw1Pq5J6HDh2S3W7XlVdeWeU5H330kb788kvt2bNHZnNFqF27dq3i4uK0fft2DRkypNJnnnvuOU2cOFEJCQmSpC5duujzzz/X8uXLncLk2LFjdeeddzp9trpJLkFBQY7fv/DCC/rd737nuMZDDz2k9PR0HTx40KmWe++9VxMmTJAkzZ8/Xx9//LGee+45rVmzptL133jjDZ07d04rV65Uy5Yt1b17dz388MOaMWPGRevyZ4RJAACagIDynCrac+vsnna7vdpzDhw4IJPJ5AiSktSpUyeZTCbt37/fZZg8cOCAJk+e7NQ2YMAAvfvuu05tcXFxlT57xRVX1LB66ZtvvqkURvv06eMIk2fPnlVOTo769+9fqZaqJhcdOHBAPXr0cHpSe+2119a4Jn9EmAQAoAmwGU1Sqav2qDq7Z+fOnWUwGPTNN99UeY7dbnc5lC2pyvaqjl3YFhwcXOkcd4a5q6vhYqr6XE0CdkNDmAQANBlW61klJWUqJ6dIJlMLWSx9m8z6nSUhFhlLM52GusuN0SoJsdTZPcPDw3XjjTcqJSVFM2bMqPTeZF5enrp27arvv/9eVqvV8XTyu+++U05Ojrp27eryurGxsdq1a5fuuOMOR9vOnTurPP+X3BnmvvLKK/XZZ5/p9ttvd7R99tlnjt+HhobKZDJp165dGjx4sFMtsbGxLq/ftWtXvfrqqyosLHSE3X//+9/V1u3PCJMAgCbBaj2rMWPeVXb2WUdbZuaJJrOzkD3QrMLWaQrKT1JAea5sxiiVhFjqbPLNeU8//bSGDx+uG264QfPnz1ePHj1kt9uVnp6uZ599Vvv27VPPnj2VkJCgxYsXy26369FHH9VVV12l+Ph4l9e877779Pvf/15XX321hg4dqg8++EBvvPGG/vrXv1ZbjzvD3DNnztQ999yjuLg4DRw4UG+99ZYyMzMVFhbmVMuiRYvUuXNnXX311Xrttde0c+dOffTRRy6vOX78eD355JO699579eijjyo3N1dLly6tcU3+iDAJAGgSkpIynYKkJGVnVzypTEkZ6qOq6pc90Fxnk22q0qlTJ23fvl1Lly5VYmKicnJy1Lp1a/Xs2VPPPvusDAaDUlNTNXfuXI0aNUqSNHjwYC1ZsqTKoeJRo0ZpyZIlWrFihebNm6cOHTpo6dKlTpNvvGHcuHH67rvv9MQTT6i4uFijRo3S3XffrXfeecdxzsyZM1VQUKDExEQdP35cMTEx+stf/uJyJrcktWzZUq+99poeeughDR48WDExMfrjH/+o2267zau11ydDXl5e4xu8dyErK0sxMTG+LqNBoc/cR5+5jz5zH33mvqysLD344AHt2FF5EsqgQSZt2TLKB1WhIZo8ebLKysr02muv+boUv8GTSQBAk2AytXDZHhXluh0oKirSiy++qF//+tcKDAzU5s2b9c477+gvf/mLr0vzK4RJAECTYLH0VWbmCaeh7ujoUFksfX1YFfyZwWDQBx98oGeeeUYlJSW64oortHr1at18882+Ls2vECYBoIkzlFl/mpSRI5vRVC+TMnzBbA5VWtpIJSVlKje3SFFRTWs2N9zXvHlzbdq0yddl+D3CJAA0Yb7YYs+XzObQJjPZBqgvAb4uAADgOxfbYg8AaoInkwDgp+pjgW1fbLEHoHEhTAKAH7rYAtve5Ist9gA0LgxzA4AfutgC2+4ylFnV/PR0BZ8cpeanp8tQZnUcKwmxqNwY7XR+XW+xB6Bx4ckkAPihnJwil+25ua7bq1LdBBtfbbEHoPGo0ZPJjIwMTZo0Sd26dVNYWJhSU1OdjhcUFOiRRx5R9+7dFRUVpb59+2rlypV1UjAANAXeWmC7JhNszm+xV3jZFhWHpxAkAbilRmGysLBQ3bt3V3Jyspo3b17p+Pz58/X+++/rhRde0O7du/Xwww/riSee0N///nevFwwATYHF0lfR0c6TbTxZYJsJNgDqWo3C5LBhw7RgwQKNHj1aAQGVP/LJJ59o4sSJio+Pl9ls1m233aa+ffvq008/9XrBANAUnF9ge8KEzho0yKQJEzorLW2k27O5bUZTFe1MsAHgHV55Z7J///765z//qTvvvFPt27fX7t279eWXX+p///d/vXF5APC6+lh2p7a8scB2SYhFxtJMp6FuJtgA8CZDXl6e3Z0PtGvXTkuWLNHkyZMdbefOndODDz6o1NRUBQZW5NMlS5bo7rvvrvI6WVlZHpYMALVz9Gix7r33Cx05UuJoa98+SM89d5Xatav8Kk9D18xwVO2avaBLDCdUao/Q0XMzdc7eztdlNRoxMTG+LgHwKa88mVy9erV2796tV199VR06dNDHH3+sxx9/XB07dtSvf/1rl5+p7x++rKwsfuDdRJ+5ryn3madP+nzRZ0uWbHUKkpJ05EiJUlNPNoit9tzvsxhJQ2RXxX/0m+L0mqb8swnUtVqHyeLiYi1cuFDr1q3TyJEVi+n27NlT+/bt04oVK6oMkwAaj4stsO1vQ8eS95bdAQB4YdHy0tJSlZaWymg0OrUbjUbZbLbaXh5AA+DNBbbrg7eW3QEA1PDJZEFBgQ4dOiRJstlsOnLkiPbu3avw8HB16NBB1113nZ544gkFBwerQ4cOysjI0N///nc98cQTdVo8AP/Q0J70WSx9lZl5wikAe7LsDgCghk8m9+zZo/j4eMXHx6u4uFiLFi1SfHy8nnrqKUnSSy+9pLi4OCUkJKh///5atmyZ5s+fr4SEhDotHoB/aGhP+ry17E5dCyjJUMtjvRWS01Etj/VWQEmGr0sCgEpq9GRy0KBBysvLq/J4ZGSknn/+eW/VBKCBaYhP+ryx7M7F1HbpoYCSDLU8PVoGlVU02M6q5enRKgjfJOnyuikaADzA3twAau38k76kpEzl5hYpKso/122sL96YkNTizKyfg+RPDCpTizOzJP3Dm+UCQK0QJgF4RV0/6WtILjYhqaZ9ZLDlVdF+prblAYBX1Xo2NwDAmTcmJNkDwqpob+VJSQBQZwiTAOBl3piQVNRqlewXDB7ZFaiiVqtqVRsAeBthEgC8zGLpq+ho53cj3Z2QZAu6TgXhm1Qe0FE2tVJ5QEcVhG+SLeg6b5cLALXCO5MA4GXempBkC7pOBUF766hKAPAOwiQA1IFfTkgylFkVlP+wAk7myGY0qSTEIntgU9whG0BjRJgEgDpiKLMq6Ow8XfLjhzLox4rGUslYmqnC1mkESgCNAu9MAkAdCCjJUMuTA9Xsx3d+DpI/MZZnKyg/yUeVAYB38WQSALzIUGZV0JnHdMm592SQrcrzAspz67EqAKg7hEkA8BJDmVXBP4yRsTy72nNtxqh6qAgA6h5hEoBX1HYv6sYgKD+pRkGy3BitkhBLPVQEAHWPMAmg1ryxF3VjEFCec9HjdgWotNlwlbRKZvINgEaDCTgAau1ie1E3JTajqcpjdgWrIHyLitu86ndB0mo9q+nTt2rUqLc0ffpWWa1nq/8QAPyEJ5MAas0be1E3BiUhFhlLM52Guu0KUumlQ1USusjvQqTEU2UAtceTSQC15o29qBsDe6BZha3TdC5ogsouGaRzQROUH7Fbxa3X+2WQlHiqDKD2eDIJoNYslr7KzDzhFErc3Yu6sbAHmlUcnuLrMmqMp8oAaoswCaDWvLUXNeofT5UB1BbD3AC8ym73dQW40MUm2FgsfRUd7Rz6m+pTZQCe4ckkgFpjEof/qu57w1NlALVFmARQaxebxJGSMtRHVUGq2ffGbA7l+wTAY4RJALXGJA7X/GFXIL43AOoaYRJArTXESRyGMquC8pMUUJ4jm9GkkhCLV5fv8Zeh/4b4vQHQsDABB0CtNbRJHIYyq4J/GKNmJW8osHSHmpW8oeAfxshQZvXaPfxl/caG9r0B0PDwZBJArTW0SRxB+UlOu9RIkrE8W0H5SV5bI9Jfhpcb2vcGQMNDmATgFQ1pEkdAeU4V7bleu4c/DS83pO8NgIaHYW4ATY7NaKqiPcpr92B4GUBTwZNJAE1OSYhFxtJMp6HucmO0SkIsXrsHw8sAmgrCJIAmxx5oVmHrtJ9mc+fKZozy+mxuieFlAE0DYRJAk2QPNHttsg0ANGW8MwkAAACPESYBAADgMYa5AfiM1XpWjz/+HxUUHPDZdoMAgNohTALwCX/ZbhAAUDsMcwPwCU+3GzSUWdX89HQFnxyl5qene3ULRACA+3gyCcAnPNlu0FBm1aXHf6tmhp8CZKlkL/5EP16+2evL+gAAaqZGTyYzMjI0adIkdevWTWFhYUpNTa10zsGDB3X77berY8eOMplMio+P14EDB7xeMIDGwZPtBstzE3WpwflJ5KUGq8pzE71aGwCg5moUJgsLC9W9e3clJyerefPmlY5/9913Gj58uMxmszZv3qydO3fKYrEoODjY6wUDaByq224woCRDLY/1VkhOR7U81lsBJRk6aj3o8lpHqmgHANS9Gg1zDxs2TMOGDZMkzZ49u9LxpKQkDR06VH/6058cbZ06dfJOhQAapfPbDc6du02FhUan7QYDSjLU8vRoGVRWcbLtrFqeHq2i4liX18o5HqqO9Vg7AOBntZ6AY7PZ9M9//lOxsbEaN26cOnfurBtuuEEbN270Rn0AGjGzOVRPPtldW7aMUkrKUMcs7hZnZv0cJH9iUJmi232vg9bWTu0Hra216eM76q3m6litZzV9+laNGvWWpk/fKqv1bPUfAoAGzJCXl2d35wPt2rXTkiVLNHnyZEnSsWPHFBsbqxYtWugPf/iD4uPj9a9//UuJiYlKTU3ViBEjXF4nKyur9tUDaJSubnGDAgMKKrWfK2+pwXc9rHsmblbbiLP6/kSoVr72Wz32+I1q167yKzj17ejRYt177xc6cqTE0da+fZCee+4qv6gPdSMmJsbXJQA+VevZ3DabTZL0m9/8Rvfee68kqXfv3vr888+1du3aKsNkff/wZWVl8QPvJvrMffSZewxlVpV+P1etWhTIZjSpJMQie6BZhmOtJVvlMGm8pLVWvzRVSUlXKTe3SFFRLbT6Jf9Z6HzJkq1OQVKSjhwpUWrqSaWkDPXaffh75j76DKg7tQ6Tbdq0UWBgoGJjnd9luvLKKxnqBvyE1VqxfmNOTpHHO8144xq/ZCizKviHMTJeki2VSiqVjKWZKmydpqJWq5zfmZRkV6CKWq2SOTLUq8HMm6pa7ig7u2Lo21t9BwD+pNZhslmzZvrVr35Vadj64MGD6tChQ20vD6CWvLHTTF3sVhOUnyRjebZTm7E8W0H5SSoOT1FB+KaKdydtZ2QPaKWiVqtkC7rOo3vVl6qWO/r669PKzDzh+DM7/QBoTGo0AaegoEB79+7V3r17ZbPZdOTIEe3du1eHDx+WJP3v//6v3nzzTa1bt06HDh3SK6+8oo0bN2ratGl1WjyA6nm604y3r3GhgPKcKtpzJUm2oOtUELlX+SarCiL3+n2QlFwvdxQcHKjCQufJRLXtOwDwJzUKk3v27FF8fLzi4+NVXFysRYsWKT4+Xk899ZQkadSoUVq2bJlWrFihgQMHavXq1XrhhRc0fPjwOi0eQPU82WmmLq5xIZvRVEV7lMfX9LXzyx1NmNBZgwaZNGFCZ3XrFu7y3Nr0HQD4kxoNcw8aNEh5eXkXPWfy5MmOGd4A/IcnO83UxTUuVBJikbE002mou9wYrZIQi8fX9Adms/M7ndOnb3Ua4j6vNn0HAP6k1utMAvBv1e00U1/XuJA90KzC1mk6VTpCZZcM0rmgCSpsndbo9tiui74DAH9S6wk4APzb+aHXpKRMx3I67s4m9sY1XLEHmpX945OK6dh4l2ypq74DAH9BmASaELtbWxQ4u3D49kLGwn+oxdn7ZNCPsutSFYWuUHnwOI/v5+2liHypur4DgIaMMAk0cnWxrM+FjIX/UPDZqTL89GeDihR8dqoKJY8CZX3UDADwDt6ZBBq5uljW50IVTySdGX5q90R91AwA8A7CJNDI1cWyPhcy6Ee32qtTHzUDALyDMAk0cnWxrM+F7LrUrfbq1EfNAADvIEwCjZy3lqYxlFnV/PR0BZ8cpeanp8tQZnUcKwpdoQvn9th/avdlzQCAuscEHKCR88bSNIYyq4J/GPPzAuOlkrE007EuZHnwOBVKXpvNzXI6ANBwECaBJqC2S9ME5Sc57VQjScbybAXlJ6k4PEVSxazt/FosBXQhltMBgIaBYW4A1Qooz6miPbeeKwEA+BvCJIBq2YymKtqj6rkSAIC/IUwCkHTxCTYlIRaVG6Odzi83RqskxFLfZQIA/AzvTAKodoKNPdCswtZpCspPUkB5rmzGKJWEWGQPNPu2cACAzxEmAdRogo090Oz4PQAA5zHMDYAJNgAAjxEmATDBBgDgMcIk0AQElGSo5bHeCsnpqJbHeiugJMPpOBNsAACeIkwCjVxASYZanh4to+2/CtBZGW3/VcvTo50C5fkJNueCJqjskkE6FzTBMfkGAICLYQIO0Mi1ODNLBpU5tRlUphZnZqkgaK+jzd8m2FitZ5WUlKmcnCKZTGynCAD+ijAJNHIGW14V7WfqtxA3WK1nNWbMu8rOPutoy8w8obS0kQRKAPAzDHMDjZw9IKyK9lb1W4gbkpIynYKkJGVnVzypBAD4F8Ik0MgVtVol+wWDEHYFqqjVKh9VVL2cnCKX7bm5rtsBAL5DmAQaOVvQdSoI36TygI6yqZXKAzqqIHyTbEHX+bq0KplMLVy2R0W5bgcA+A7vTAJNgC3oOqfJNv7OYumrzMwTTkPd0dGhslj6+rAqAIArhEkAfsdsDlVa2kglJWUqN7dIUVHM5gYAf0WYBOCXzOZQpaQM9XUZAIBq8M4kAAAAPEaYBAAAgMcIkwAAAPAYYRIAAAAeI0wCAADAY4RJAAAAeIwwCQAAAI8RJgEAAOAxwiQAAAA8RpgEAACAx2oUJjMyMjRp0iR169ZNYWFhSk1NrfLc+++/X2FhYVqxYoXXigQAAIB/qlGYLCwsVPfu3ZWcnKzmzZtXed6mTZv02WefyWQyea1AAAAA+K8ahclhw4ZpwYIFGj16tAICXH/kv//9rx577DGtXbtWgYGBXi0SAAAA/skrqa+srEzTpk3TnDlzFBsbW6PPZGVleePWbvHFPRs6+sx9DbXPjh4t1gsvZOvEiR8VEXGpZs6MVrt2VY9EeFND7TNfos/cV1d9FhMTUyfXBRoKr4TJRYsWKTw8XFOnTq3xZ+r7hy8rK4sfeDfRZ+5rqH1mtZ7Vgw++q+zss462AwdKlJY2UmZzaJ3eu6H2mS/RZ+6jz4C6U+vZ3Dt27ND69eu1cuVKb9QDwAeSkjKdgqQkZWefVVJSpo8qAgA0FLUOk+np6crNzVVsbKzatGmjNm3a6PDhw0pMTFT37t29USOAOpaTU+SyPTfXdTsAAOfVeph72rRpGj16tFPbuHHjNG7cON111121vTyAemAytXDZHhXluh0AgPNqFCYLCgp06NAhSZLNZtORI0e0d+9ehYeHq0OHDoqIiHC+aGCgIiMjeT8FaCAslr7KzDzhNNQdHR0qi6WvD6sCADQENRrm3rNnj+Lj4xUfH6/i4mItWrRI8fHxeuqpp+q6PgD1wGwOVVraSE2Y0FmDBpk0YULnepl8AwBo+Gr0ZHLQoEHKy8ur8UX37dvnaT0AfMRsDlVKylBflwEAaGDYmxsAAAAeI0wCAADAY4RJAAAAeIwwCQAAAI8RJgEAAOAxwiQAAAA8RpgEAACAxwiTAAAA8Fit9+YGzrNazyopKVM5OUUymVrIYunLDioAADRyhEl4hdV6VmPGvOu0t3Nm5gm25AMAoJFjmBtekZSU6RQkJSk7u+JJJQAAaLwIk/CKnJwil+25ua7bAQBA40CYhFeYTC1ctkdFuW4HAACNA2ESNWK1ntX06Vs1atRbmj59q6xW5yFti6WvoqOd342Mjg6VxdLXresAAICGhQk4qFZNJteYzaFKSxuppKRM5eYWKSqq8mxuJukAAND4ECZRrYtNrklJGepoM5tDnf7s6XUAAEDDwTA3quWtyTVM0gEAoPHhySSq5a3JNbW9jqHMqqD8JAWU58hmNKkkxCJ7oNmtGgAAgHfxZBLVqunkmrq8jqHMquAfxqhZyRsKLN2hZiVvKPiHMTKUWd2qAQAAeBdPJlGtmkyuqevrBOUnyVie7dRmLM9WUH6SisNT3KoDAAB4D2ESNVLd5Jq6vk5AeU4V7bm1LQkAANQCw9xoEPJLLquivU09VwIAAH6JMIkGwbJ8uA5aWzu1HbS2lmX5cB9VBAAAJIa50UDs/TpEv56aoKT731PbiLP6/kSoLMuHq+MVIb4uDQCAJo0wiQbBZGqhHTva6I5Hf+fUfu1A9v4GAMCXGOZGg+Ct5YkAAIB38WQS9crThce9tTwRAADwLsIkasRqrdhDOyenSCaTZ0Hu/MLjjvUiSyVjaaYKW6fVOFCyhzcAAP6FMIlqWa1nNWbMu8rOPutoy8w8obS0kW4FShYeBwCg8eGdSVQrKSnTKUhKUnZ2xZNKd7DwOAAAjQ9hEtXKySly2Z6b67q9KjajqYr2KLdrAgAA/oEwiWqZTK6X34mKcm9ZnpIQi8qN0U5t5cZolYRYPK4NAAD4FmES1arpsjyGMquan56u4JOj1Pz0dBnKrE7H7YFmFbZO07mgCSq7ZJDOBU2o8eQbAADgn5iAg2rVZFmems7UtgeamWwDAEAjQphEjVS3LA8ztQEAaJoY5oZXMFMbAICmqUZhMiMjQ5MmTVK3bt0UFham1NRUx7HS0lIlJiZq4MCBatu2rWJjYzVt2jQdPny4zoqG/2GmNgAATVONwmRhYaG6d++u5ORkNW/e3OlYUVGRvvjiC82ZM0fbt2/X+vXrdfToUY0fP15lZWV1UjTqX0BJhloe662QnI5qeay3AkoynI4zUxsAgKapRu9MDhs2TMOGDZMkzZ492+lYq1atlJaW5tT27LPPqn///jpw4IB69OjhnUrhMwElGWp5erQM+ukfB7azanl6tArCN8kWdJ2kn2dqV+y7nSubMarG+24DAICGq04m4OTn50uSwsLC6uLyqGctzsz6OUj+xKAytTgzSwVBex1tzNQGAKDp8XqYPHfunCwWi0aMGKF27dpVeV5WVpa3b10tX9yzocvKytLVLX5w+UKEvewH+tQF+sR99Jn76DP31VWfxcTE1Ml1gYbCq2GyrKxMCQkJOnPmjF599dWLnlvfP3xZWVn8wLvpfJ8ZjrWWbAWVjhsCW9OnF+DvmfvoM/fRZ+6jz4C647WlgcrKyjR16lR99dVX2rRpk1q3bu2tS8PHilqtkv2Cf3fYFaiiVqt8VBEAAPAXXgmTpaWlmjJlir766itt2bJFkZGR3rgs/IQt6DoVhG9SeUBH2dRK5QEdnSbfAACApqtGw9wFBQU6dOiQJMlms+nIkSPau3evwsPDZTKZdNddd2nPnj169dVXZTAYdOzYMUlSaGhopaWE0DDZgq5zmmwDAAAg1fDJ5J49exQfH6/4+HgVFxdr0aJFio+P11NPPaWjR4/qnXfeUU5OjoYMGaLY2FjHr40bN9Z1/QAAAPChGj2ZHDRokPLy8qo8frFjAAAAaLzYmxsAAAAeI0wCAADAY4RJAAAAeIwwCQAAAI8RJgEAAOAxwiQAAAA8RpgEAACAxwiTAAAA8BhhEgAAAB4jTAIAAMBjhEkAAAB4jDAJAAAAjxEmAQAA4LFAXxcAZ4Yyq4LykxRQniOb0aSSEIvsgWZflwUAAOASYdKPGMqsCv5hjIzl2RUNpZKxNFOFrdMIlAAAwC8xzO1HgvKTfg6SPzGWZysoP8lHFQEAAFwcYdKPBJTnVNGeW8+VAAAA1Axh0o/YjKYq2qPquRIAAICaIUz6kZIQi8qN0U5t5cZolYRYfFQRAADAxTEBx4/YA80qbJ3202zuXNmMUczmBgAAfo0w6WfsgWYVh6f4ugwAAIAaYZgbAAAAHiNMAgAAwGOESQAAAHiMMAkAAACPESYBAADgMcIkAAAAPEaYBAAAgMcIkwAAAPAYYRIAAAAeI0wCAADAY4RJAAAAeIy9uf2M1XpWSUmZyskpksnUQhZLX5nNob4uCwAAwCXCpB+xWs9qzJh3lZ191tGWmXlCaWkjCZQAAMAvMcztR5KSMp2CpCRlZ1c8qQQAAPBHhEk/kpNT5LI9N9d1OwAAgK/VKExmZGRo0qRJ6tatm8LCwpSamup03G63a9GiReratauioqJ000036euvv66Tghszk6mFy/aoKNftAAAAvlajMFlYWKju3bsrOTlZzZs3r3R8+fLlWrlypRYvXqytW7cqIiJCY8eOVX5+vtcLbswslr6KjnZ+NzI6OlQWS18fVQQAAHBxNQqTw4YN04IFCzR69GgFBDh/xG63a9WqVXrggQc0evRode/eXatWrVJBQYE2bNhQJ0U3VmZzqNLSRmrChM4aNMikCRM6M/kGAAD4tVrP5rZarTp27JiGDh3qaGvevLkGDhyo3bt3a8qUKbW9RZNiNocqJWVo9ScCAAD4gVqHyWPHjkmSIiIinNojIiKUk5NT5eeysrJqe2u3+eKeDR195j76zH30mfvoM/fVVZ/FxMTUyXWBhsJr60waDAanP9vt9kptv1TfP3xZWVn8wLuJPnMffeY++sx99Jn76DOg7tR6aaDIyEhJ0vHjx53aT548WelpJQAAABqXWodJs9msyMhIbdu2zdFWUlKinTt3ql+/frW9PAAAAPxYjYa5CwoKdOjQIUmSzWbTkSNHtHfvXoWHh6tDhw6aNWuWli5dqpiYGHXp0kVPP/20goODNX78+DotHgAAAL5VozC5Z88e3XzzzY4/L1q0SIsWLdJtt92mVatW6f7771dxcbEeeeQR5eXlqU+fPtq4caNCQkLqrHAAAAD4Xo3C5KBBg5SXl1flcYPBoHnz5mnevHneqqtahjKrgvKTFFCeI5vRpJIQi+yB5nq7PwAAALw4m7s+GcqsCv5hjIzl2RUNpZKxNFOFrdMIlAAAAPWo1hNwfCEoP+nnIPkTY3m2gvKTfFQRAABA09Qgw2RAuevF0APKc+u5EgAAgKatQYZJm9FURXtUPVcCAADQtDXIMFkSYlG5MdqprdwYrZIQi48qAgAAaJoa5AQce6BZha3TfprNnSubMYrZ3AAAAD7QIMOkVBEoi8NTfF0GAABAk9Ygh7kBAADgHwiTAAAA8BhhEgAAAB4jTAIAAMBjhEkAAAB4jDAJAAAAjxEmAQAA4DHCJAAAADxGmAQAAIDHCJMAAADwGGESAAAAHiNMAgAAwGOESQAAAHiMMAkAAACPESYBAADgMcIkAAAAPEaYBAAAgMcIkwAAAPAYYRIAAAAeI0wCAADAY4RJAAAAeIwwCQAAAI8RJgEAAOAxwiQAAAA8RpgEAACAxwiTAAAA8BhhEgAAAB4jTAIAAMBjhEkAAAB4jDAJAAAAj3klTJaXlyspKUm9e/dWZGSkevfuraSkJJWVlXnj8gAAAPBTgd64yLJly7R27VqtWrVK3bt311dffaVZs2apWbNmevTRR71xCwAAAPghr4TJTz75RCNGjNDIkSMlSWazWSNHjtSnn37qjcsDAADAT3llmLt///7asWOHvvnmG0nS/v37lZ6erv/5n//xxuUBAADgpwx5eXn22l7EbrcrKSlJzzzzjIxGo8rKyjRnzhxZLJYqP5OVlVXb2wIA4HMxMTG+LgHwKa8Mc2/cuFF///vftXbtWnXt2lX79u3TY489po4dO+rOO+90+Zn6/uHLysriB95N9Jn76DP30Wfuo8/cR58BdccrYXLBggW69957NW7cOElSjx49dPjwYT377LNVhkkAAAA0fF55Z7KoqEhGo9GpzWg0ymazeePyAAAA8FNeeTI5YsQILVu2TGazWV27dtXevXu1cuVKTZo0yRuXBwAAgJ/ySphcsmSJ/vSnP+nhhx/WyZMnFRkZqbvuuos1JgEAABo5r4TJkJAQJScnKzk52RuXAwAAQAPB3twAAADwGGESAAAAHiNMAgAAwGOESQAAAHiMMAkAAACPESYBAADgMcIkAAAAPEaYBAAAgMcIkwAAAPAYYRIAAAAeI0wCAADAY4RJAAAAeIwwCQAAAI8RJgEAAOAxwiQAAAA8RpgEAACAxwiTAAAA8BhhEgAAAB4jTAIAAMBjhEkAAAB4jDAJAAAAjxEmAQAA4DHCJAAAADxGmAQAAIDHCJMAAADwGGESAAAAHiNMAgAAwGOBvi7AU1brWSUlZSonp0gmUwtZLH1lNof6uiwAAIAmpUGGSav1rMaMeVfZ2WcdbZmZJ5SWNpJACQAAUI8a5DB3UlKmU5CUpOzsiieVAAAAqD8NMkzm5BS5bM/Ndd0OAACAutEgw6TJ1MJle1SU63YAAADUjQYZJi2WvoqOdn43Mjo6VBZLXx9VBAAA0DQ1yAk4ZnOo0tJGKikpU7m5RYqKYjY3AACALzTIMClVBMqUlKG+LgMAAKBJa5DD3AAAAPAPXguTubm5mjlzpjp37qzIyEj169dPO3bs8NblAQAA4Ie8Msydl5en4cOHq3///nr99dfVpk0bWa1WRUREeOPyAAAA8FNeCZP/93//p6ioKK1evdrR1qlTJ29cGgAAAH7MK8Pcb7/9tvr06aMpU6aoS5cuuv7667VmzRrZ7XZvXB4AAAB+ypCXl1frxBcZGSlJmj17tsaMGaN9+/Zp7ty5SkxMVEJCgsvPZGVl1fa2AAD4XExMjK9LAHzKK2EyIiJCcXFxev/99x1tCxcu1FtvvaVPPvmktpf3iqysLH7g3USfuY8+cx995j76zH30GVB3vDLMHRkZqdjYWKe2K6+8UkeOHPHG5QEAAOCnvBIm+/fvr4MHDzq1HTx4UB06dPDG5QEAAOCnvBImZ8+erX//+996+umndejQIaWlpWnNmjWaNm2aNy4PAAAAP+WVdyYl6b333tPChQt18OBBtW/fXtOnT9eMGTNkMBi8cXkAAAD4Ia+FSQAAADQ97M0NAAAAjxEmAQAA4DHCJAAAADxGmAQAAIDHCJMAAADwWKMNk0uXLlVYWJgeeeQRR5vdbteiRYvUtWtXRUVF6aabbtLXX3/twyp9Lzc3VzNnzlTnzp0VGRmpfv36aceOHY7j9Jmz8vJyJSUlqXfv3oqMjFTv3r2VlJSksrIyxzlNvc8yMjI0adIkdevWTWFhYUpNTXU6XpP++fHHH/XII4/oiiuuUNu2bTVp0iQdPXq0Pr+MenWxPistLVViYqIGDhyotm3bKjY2VtOmTdPhw4edrkGfpVZ57v3336+wsDCtWLHCqb2p9RlQVxplmPz3v/+tV155RT169HBqX758uVauXKnFixdr69atioiI0NixY5Wfn++jSn0rLy9Pw4cPl91u1+uvv67du3dryZIlioiIcJxDnzlbtmyZ1q5dq8WLF+uTTz5RcnKyUlJS9MwzzzjOaep9VlhYqO7duys5OVnNmzevdLwm/TNv3jxt2bJFL774ot555x3l5+dr4sSJKi8vr88vpd5crM+Kior0xRdfaM6cOdq+fbvWr1+vo0ePavz48U7/iKHPXNu0aZM+++wzmUymSseaWp8BdaXRrTN55swZDR48WMuXL9eSJUvUvXt3/fnPf5bdblfXrl01ffp0zZkzR5JUXFysmJgYPfnkk5oyZYqPK69/CxcuVEZGht577z2Xx+mzyiZOnKjw8HC98MILjraZM2fq9OnTeu211+izC7Rr105LlizR5MmTJdXs79SZM2fUpUsXrVy5Urfeeqsk6ciRI+rVq5c2bNigG2+80WdfT324sM9c2b9/v/r376+MjAz16NGDPquiz/773/9q+PDhSktL0/jx45WQkKD77rtPkpp8nwHe1OieTD7wwAMaPXq0Bg8e7NRutVp17NgxDR061NHWvHlzDRw4ULt3767vMv3C22+/rT59+mjKlCnq0qWLrr/+eq1Zs0Z2e8W/L+izyvr3768dO3bom2++kVTxP/X09HT9z//8jyT6rDo16Z/PP/9cpaWlTue0b99esbGx9OFPzj/FDQsLk0SfuVJWVqZp06Zpzpw5io2NrXScPgO8J9DXBXjTK6+8okOHDmn16tWVjh07dkySnIZwz/85JyenXurzN999951efPFFzZ49Ww888ID27dunuXPnSpISEhLoMxceeOABFRQUqF+/fjIajSorK9OcOXMc+9DTZxdXk/45fvy4jEaj2rRpU+mc48eP10+hfuzcuXOyWCwaMWKE2rVrJ4k+c2XRokUKDw/X1KlTXR6nzwDvaTRhMisrSwsXLtS7776rZs2aVXnehXuF2+32Jrt/uM1mU1xcnBITEyVJV111lQ4dOqS1a9cqISHBcR599rONGzfq73//u9auXauuXbtq3759euyxx9SxY0fdeeedjvPos4vzpH/ow4qnbQkJCTpz5oxeffXVas9vqn22Y8cOrV+/Xunp6W5/tqn2GVAbjWaY+5NPPtGpU6c0YMAAtWnTRm3atFFGRobWrl2rNm3aqHXr1pJU6V+cJ0+erPSUpKmIjIysNPxz5ZVX6siRI47jEn32SwsWLNC9996rcePGqUePHpo0aZLuuecePfvss5Los+rUpH8uv/xylZeX69SpU1We0xSVlZVp6tSp+uqrr7Rp0ybHf9Mk+uxC6enpys3NVWxsrOP/B4cPH1ZiYqK6d+8uiT4DvKnRhMmbbrpJH3/8sdLT0x2/4uLiNG7cOKWnp6tLly6KjIzUtm3bHJ8pKSnRzp071a9fPx9W7jv9+/fXwYMHndoOHjyoDh06SJLMZjN9doGioiIZjUanNqPRKJvNJok+q05N+ufqq6/WJZdc4nTO0aNHdeDAgSbbh6WlpZoyZYq++uorbdmyxRHKz6PPnE2bNk0ZGRlO/z8wmUyaPXu2Nm3aJIk+A7yp0Qxzh4WFOV5GP69FixYKDw93/Et01qxZWrp0qWJiYtSlSxc9/fTTCg4O1vjx431Qse/Nnj1bw4YN09NPP61bbrlFe/fu1Zo1a/T4449LqhiKpM+cjRgxQsuWLZPZbFbXrl21d+9erVy5UpMmTZJEn0lSQUGBDh06JKniVYojR45o7969Cg8PV4cOHartn1atWumOO+7QggULFBERofDwcM2fP189evTQkCFDfPiV1Z2L9ZnJZNJdd92lPXv26NVXX5XBYHC8exoaGqrmzZvTZy7+nl34dDEwMFCRkZGKiYmR1DT/ngF1pdEtDfRLN910k2NpIKniXZjk5GStW7dOeXl56tOnj55++mlH2GyK3nvvPS1cuFAHDx5U+/btNX36dM2YMcPxzhB95iw/P19/+tOf9NZbb+nkyZOKjIzUuHHj9OijjyooKEgSfZaenq6bb765Uvttt92mVatW1ah/SkpK9Pjjj2vDhg0qKSlRfHy8li5dqvbt29fnl1JvLtZnjz32mK666iqXn1u5cqVjORz6rML5v2cX6tWrl9PSQFLT6zOgrjTqMAkAAIC61WjemQQAAED9I0wCAADAY4RJAAAAeIwwCQAAAI8RJgEAAOAxwiQAAAA8RpgE6klqaqrCwsJktVrd/qzValVYWJhj28b6dv7+qampTu2ff/65Ro4cqfbt2yssLMyjvZDdMWvWLPXq1atO7wEAcA9hEk3a5s2bFRYWpg0bNlQ6dvPNN1/0mNlslt3uf8u07ty5U4sWLVJeXl6Nzj8fcv/973+7dZ/y8nJNmTJF33//vRYuXKjVq1crNjZWr732mp5//nkPKgcANESESTRpAwYMkFQRwH6prKxMn376qQIDA6s81r9/f8dOQTUxadIk5ebmqmPHjrUv/CJ27dqlxYsX68yZM167ZseOHZWbm+vYNlKSjhw5ouzsbM2YMUN33323Jk6cqMsvv1yvv/66yx1IAACNU6PZmxvwREREhDp37lwpMH7xxRcqKirSrbfeWuWx/v37u3Uvo9Eoo9FY65p9wWAwOLaLPO/kyZOSKvY4BgA0XTyZRJM3YMAA7d+/32lYeNeuXTKZTJo4caLLY+c/d962bds0atQotW/fXm3bttWoUaO0e/dup/tU9c7kunXrFBcXp8jISF1//fX65z//edF3A1999VVdc801uvzyyzVw4EB99NFHjmOLFi3SE088IUm66qqrFBYW5pV3GS98Z3LWrFm68cYbJUn33HOPwsLC1KtXL91000368MMPdfjwYce9w8LCHNex2+1as2aNBg4cqMjISEVHR2v69Ok6evRopXv+7W9/U58+fRQZGanrrrtO7777bq2+BgBA3eDJJJq8/v37629/+5s++eQTDRs2TFJFYOzXr5+uueYaSap0LCgoSHFxcZKkDRs2KCEhQYMGDdL8+fNls9mUmpqq3/72t3r77bfVt2/fKu+9bt06PfDAA7rmmmuUkJCgkydPasaMGWrXrp3L8zdt2qRTp05pypQpCgoK0qpVq3T77bdr3759Cg8P180336ysrCxt3LhRTz31lNq0aSNJio2N9Vp/SdKUKVNkNpuVnJys3//+9xowYICCg4MVHBysvLw85ebm6qmnnqr0uYceekh/+ctfNHHiRE2bNk3Hjh3TmjVrtHv3bv3rX/9yBM/169fr3nvv1a9+9StNmzZNJ06c0IwZM9S+fXuvfh0AgNojTKLJO/+EcdeuXY7AuHv3bj344IMKDQ1V165dKx2Li4vTpZdeqsLCQs2ZM0cTJ050ek9wypQp6t+/vxYuXKjNmze7vG9paamefPJJ9ezZU2+//baaNWsmSYqPj9fo0aPVoUOHSp/Jzs7Wp59+qssuu0ySdP311ys+Pl4bNmzQ9OnT1bNnT/Xq1UsbN27UTTfdJLPZ7L2O+oVrr71WBoNBycnJuuaaazRx4kTHsaioKJ09e9apTarot5dfflkrV67U5MmTHe0333yzhgwZojVr1ujRRx9VWVmZ/vjHP6pr16565513HMPr119/vW655RaX/QIA8B2GudHkde7cWZGRkY53I7/99lsdP37c8U5k//79Kx0bOHCgpIrh7by8PN166606deqU41dxcbGGDBminTt3qrS01OV9P/vsM506dUq///3vHUFSkgYPHqxu3bq5/MyYMWMcQVKSevfurdDQUH333Xe17oe69uabb6ply5YaNmyYU1+ZTCZ17txZ//rXvyRV9Mvx48cdT1/PGzp0qLp27eqr8gEAVeDJJCCpX79+ev/993Xu3Dnt2rVLLVq0cLyz2K9fP61fv95xTJIjaH777beSpLFjx1Z57TNnzjgFwPMOHz4sqSLMXqhz58764osvKrW7eirXqlUrnT59urov0ee+/fZbFRQUKCYmxuXx8zPjz/eLq/O6dOnisl8AAL5DmARUEQ43b96sPXv2aNeuXerTp48CAyt+PPr166eSkhLHsYCAAF177bWSJJvNJkl6/vnn1bZtW5fXDg0NdbueqtavrGo2uD+ud3khm82m1q1b66WXXnJ5vEWLFpJ+/lpcLbvUEL5OAGhqCJOA5Bi23rVrl3bt2qUxY8Y4jnXq1ElRUVGOYz169HAshxMdHS1JuuyyyzRkyBC37nn+KeO3336rG264wenYoUOHPPxKXIew+lTV/aOjo7Vt2zb16dNHISEhVX7+/Dqc33zzTaV+Of8kGADgP3hnEpDUq1cvtWzZUm+//baysrIqrSHZr18/l8duvPFGtWrVSk8//bR+/PHHStc9vxajK3FxcWrTpo3WrVunc+fOOdq3b9+ur7/+2uOv5fwTvprugONtLVq0cLlg+i233CKbzabk5ORKx+x2u06dOiWpol8iIiK0bt06lZSUOM7ZunWr9u/fX3eFAwA8wpNJQBXDx3379tVHH32kgICASsv59OvXT3/4wx8k/fwUU5JCQkK0fPlyTZ06Vddff70mTJigyMhIHT16VOnp6QoODna5HaMkNWvWTPPnz9dDDz2km266SePGjdPJkyeVkpKi7t27q6CgwKOv5fySRU8++aTGjRunZs2aKT4+XhERERf93Pr1653WrDzvrrvucvv+mzdv1ty5c9W3b18FBARo3LhxGjhwoGbMmKGVK1fqyy+/1K9//Wu1aNFCVqtVb731lu644w49+OCDuuSSS7RgwQLdd999+s1vfqMJEyY4+qVbt24e9wsAoG4QJoGfDBgwQB999JG6detWaVeXXz6NvPCp5ZgxY2QymfTMM8/o+eefV3FxsSIjI9W3b1/deeedF73n3XffLUn6v//7PyUmJiomJkarV6/W+vXrPX4Kd80118hisWjdunW65557ZLPZtGXLlmrD5Msvv+yyffjw4W7tcpOQkKD9+/fr9ddf15o1a2S32zVu3DhJ0uLFi3X11VfrxRdf1KJFixQQEKC2bdvqxhtv1KhRoxzXuOOOO2S327Vs2TIlJiaqS5cuWr16tTZv3qwdO3bUuBYAQN0z5OXl8UY74Geuu+46RUREKC0tzdelAABwUbwzCfhQSUlJpRnK27dv11dffaX4+HgfVQUAQM3xZBLwofT0dM2dO1e//e1vFRUVpa+//lrr1q1TmzZt9PHHHzvtaw0AgD/inUnAhzp27Ciz2ayXX35ZP/zwg0JDQzVq1CgtWLCAIAkAaBB4MgkAAACP8c4kAAAAPEaYBAAAgMcIkwAAAPAYYRIAAAAeI0wCAADAY/8P0GFNAlX1+yoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fitted = fit(shotput, 'Weight Lifted', 'Shot Put Distance')\n", "\n", "shotput['Best Straight Line'] = fitted\n", "\n", "fig, ax = plt.subplots(figsize=(7,6))\n", "\n", "ax.scatter(shotput['Weight Lifted'], \n", " shotput['Shot Put Distance'], \n", " label='Color=darkblue', \n", " color='darkblue')\n", "\n", "ax.scatter(shotput['Weight Lifted'], \n", " shotput['Best Straight Line'], \n", " label='Color=gold', \n", " color='gold')\n", "\n", "x_label = 'Weight Lifted'\n", "\n", "y_label = ''\n", "\n", "y_vals = ax.get_yticks()\n", "\n", "plt.ylabel(y_label)\n", "\n", "ax.legend(bbox_to_anchor=(1.04,1), loc=\"upper left\", frameon=False)\n", "\n", "plt.xlabel(x_label)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Nonlinear Regression\n", "The graph above reinforces our earlier observation that the scatter plot is a bit curved. So it is better to fit a curve than a straight line. The [study](http://digitalcommons.wku.edu/ijes/vol6/iss2/10/) postulated a quadratic relation between the weight lifted and the shot put distance. So let's use quadratic functions as our predictors and see if we can find the best one. \n", "\n", "We have to find the best quadratic function among all quadratic functions, instead of the best straight line among all straight lines. The method of least squares allows us to do this.\n", "\n", "The mathematics of this minimization is complicated and not easy to see just by examining the scatter plot. But numerical minimization is just as easy as it was with linear predictors! We can get the best quadratic predictor by once again using `minimize`. Let's see how this works." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recall that a quadratic function has the form\n", "\n", "$$\n", "f(x) ~=~ ax^2 + bx + c\n", "$$\n", "for constants $a$, $b$, and $c$.\n", "\n", "To find the best quadratic function to predict distance based on weight lifted, using the criterion of least squares, we will first write a function that takes the three constants as its arguments, calculates the fitted values by using the quadratic function above, and then returns the mean squared error. \n", "\n", "The function is called `shotput_quadratic_mse`. Notice that the definition is analogous to that of `lw_mse`, except that the fitted values are based on a quadratic function instead of linear." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def shotput_quadratic_mse(a, b, c):\n", " x = shotput['Weight Lifted']\n", " y = shotput['Shot Put Distance']\n", " fitted = a*(x**2) + b*x + c\n", " return np.mean((y - fitted) ** 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now use `minimize` just as before to find the constants that minimize the mean squared error. " ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-1.04004837e-03, 2.82708043e-01, -1.53182103e+00])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "best = minimize(shotput_quadratic_mse)\n", "best" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our prediction of the shot put distance for an athlete who lifts $x$ kilograms is about\n", "$$\n", "-0.00104x^2 ~+~ 0.2827x - 1.5318\n", "$$\n", "meters. For example, if the athlete can lift 100 kilograms, the predicted distance is 16.33 meters. On the scatter plot, that's near the center of a vertical strip around 100 kilograms." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "16.3382" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(-0.00104)*(100**2) + 0.2827*100 - 1.5318" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here are the predictions for all the values of `Weight Lifted`. You can see that they go through the center of the scatter plot, to a rough approximation." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "x = shotput.iloc[:,0]\n", "shotput_fit = best.item(0)*(x**2) + best.item(1)*x + best.item(2)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAApMAAAGTCAYAAAB0/plgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABBnklEQVR4nO3de1zUVeL/8fcwZHhD0BBG0ZGUwHusul5KNGu9bLZqatrazRS8ZJtblrVOshl52yxdMzOtbAvbyq+hdtn2UaaLprakpfVLw8RZMfCWqIAkMPP7g5iaGG7DwAzwej4ePh55PrczRybfnvM55xiys7PtAgAAANzg5+0KAAAAoO4iTAIAAMBthEkAAAC4jTAJAAAAtxEmAQAA4DbCJAAAANxGmAQAAIDbGkyYTEtL83YV6hzarOpos6qjzaqONqs62gyoOQ0mTAIAAMDzCJMAAABwG2ESAAAAbiNMAgAAwG2ESQAAALiNMAkAAAC3ESYBAADgNsIkAAAA3EaYBAAAgNsIkwAAAHCbv7crAABwzWq9oMTEVGVm5slkaiKLpbfM5kBvVwsAnNAzCQA+yGq9oNGjP9Dbb3+nnTsz9fbb32n06A9ktV7wdtWAWpeUlKS2bdt6uxqSpJUrV6p79+7Vvk9KSoqCgoJ09uxZl8etVquCgoK0f//+Mu9RmXNqA2ESAHxQYmKq0tOdg2N6enFPJVDXnDp1SnPnztW1116r1q1bq3Pnzho3bpz+/e9/e7tq8ACGuQHAB2Vm5rksz8pyXQ74KqvVquHDh6tZs2ZKSEhQt27dZLPZtGPHDj344IP66quvaqUely9fVqNGjWrlWRW5fPmyt6vgUfRMAoAPMpmauCwPC3NdDlSG1XpBcXHbNHLku4qL21Yrr03MmTNHdrtdn3zyicaMGaPIyEhFRUUpPj5eO3fulCQdP35ckyZNUnh4uMLDw3XHHXfoxIkT5d73lVdeUUxMjEJCQhQTE6NXX33V6XhQUJDWrl2rO+64Q23atNGCBQvcqv+KFSt0zTXXqG3btpo2bZpycnKcju/bt09jxozR1VdfrXbt2mn48OH67LPPqlyXH3/8UZMmTVJsbKxOnz7tKD9y5IiGDx+u0NBQ9enTR9u2bSuzrq6Gzl0NhR86dEi33XabwsPD1alTJ02ZMkUnT56sctuUIEwCgA+yWHorIsJ5sk1ERKAslt5eqhHqOm+8h3vu3Dl99NFHiouLU7NmzUodDwoKkt1u16RJk3T69Glt2bJFW7duVVZWliZNmiS73e7yvlu3btXDDz+sGTNmaPfu3Zo+fboeeughffDBB07nLVmyREOHDtWnn36qqVOnSpLatm1b7q9x48Y5rn/nnXeUmJioxx57TDt27FBkZKSef/55p2dcvHhREyZM0AcffKCPP/5Y3bt31/jx40u9C+mqLiUuXLigsWPH6ty5c3r33XcVEhLiOJaQkKBp06YpJSVFgwcP1h//+Ed9//33lWh917KysvT73/9enTt31scff6zk5GTl5OTo9ttvl81mc+ueDHMDgA8ymwOVnDxCiYmpysrKU1jYz7O509Lc70FAw1Xee7hr1w6pkWcePXpUdrtd11xzTZnnbN++XV999ZX2798vs9ksSVq3bp1iYmK0Y8cODR48uNQ1zz33nCZMmKD4+HhJUqdOnfTFF19oxYoVGjFihOO8MWPG6K677nK6NiUlpdw6BwQEOP579erVuv322zV58mRJxb2sKSkpOnr0qOOcQYMGOV2/dOlSbdmyRR999JEmTJhQZl2OHz8uSTpz5ozi4+NlMpm0fv16p+dL0r333qsxY8ZIKg6k27Zt08svvyyLxVLu5yjLSy+9pG7duumJJ55wlK1Zs0YdOnTQ/v371atXryrfkzAJAD7KbA6ssb/k0fB44z3csnoWf+nw4cMymUyOIClJHTp0kMlk0qFDh1yGycOHD2vSpElOZf379y/VMxkTE1Pq2quvvrqStS9+zp133ulU1qdPH6cwefr0aT311FNKSUnR6dOnVVRUpEuXLikjI6PCukjSrbfeqh49eui1116Tv3/pWNanTx/Hf/v5+alXr146dOhQpT/Dr3355Zf69NNPXc6OT09PJ0wCAADXvPEebseOHWUwGPTtt9+WeY7dbpfBYHB5rKzyso79uqxp06alzqloiaH+/ftr48aN5Z7zSzNmzNCpU6e0cOFCtW/fXldeeaX+8Ic/lJpk46oukjRs2DC98847+vrrr9WzZ89KP9cVP7/itxd/GeILCwudzrHZbBo6dKgSExNLXf/L4fWqIEwCANAAWCy9lZp62mmou6bfww0ODtaNN96otWvXatq0aaXem8zOzlZ0dLS+//57Wa1WR+/ksWPHlJmZqejoaJf3jYqK0p49e5x6DXfv3l3m+b9UlWHuqKgopaamOj0nNdV5ea49e/Zo8eLFGjZsmKTiZZCqMpll3rx5Cg4O1ujRo7V582b16NHD6XhqaqpjKN1ut2vfvn0aNWqUy3tdddVVkorfiyz574MHDzqd07NnT73zzjtq166drrjiikrXszyESQAAGoDy3sOtSU8//bSGDRumG264QfPmzVPXrl1lt9uVkpKiZ599VgcPHlS3bt0UHx+vJUuWyG6365FHHlHPnj0VGxvr8p7333+/7rnnHl177bUaMmSIPvroI7399tt67bXXKqxPVYa5p0+frunTp+s3v/mNrr/+em3evFmff/65goKCHOd07NhRb731lnr37q28vDzNnz+/yksQPf7447Lb7Y5A+ctF0V9++WV16tRJXbp00bp163T8+HHde++9ZX628PBwLV68WH/961/1v//9T3/729+czpk6dapeffVVTZ48WbNnz9ZVV12lY8eOOSYbNW/evEp1lwiTAAA0GN54D7dDhw7asWOHli1bpoSEBGVmZqply5bq1q2bnn32WRkMBiUlJWnu3LkaOXKkpOJJLUuXLi1zmHvkyJFaunSpVq5cqccee0zt2rXTsmXLnCbfeMKtt96qY8eO6cknn9SlS5c0YsQIzZw5Uxs2bHCc89xzz2n27NkaPHiwwsLC9Oijj5a5q0155s+fL7vdrlGjRmnz5s0KDCwO+QkJCVq1apW+/PJLtWvXTq+//nqZQ/VXXHGFXnrpJT300EO6/vrr1b17d82fP99pIpDJZNKHH36oJ554QmPHjtWPP/6o8PBw3XDDDbryyiurXG9JMmRnZ1f8dmw9kJaWpsjISG9Xo06hzaqONqs62qzqaLOqo82AmsM6kwAAAHAbYRIAAABu451JAECDYbUWL9KdmZknk6l2JqAA9R1hEgDQIJRsJ/jLpXFSU08rOXkEgRKoBoa5AQANQnnbCQJwH2ESANAgeGM7QaAhIEwCABoEb2wnCDQEhEkAQINgsfRWRITzu5E1vZ0g0BAwAQcA0CB4aztBoL6jZxIA0GCUbCe4detIrV07hCBZRyQlJZW5hWBdEBQUpM2bN5d5/OzZswoKClJKSkot1spzCJMAAKBGnTp1SnPnztW1116r1q1bq3Pnzho3bpz+/e9/e7tq8ACGuQEAQI2xWq0aPny4mjVrpoSEBHXr1k02m007duzQgw8+qK+++qpW6nH58mU1atSoVp7V0NAzCQA+ymq9oLi4bRo58l3FxW2T1Xqh4ouAchgKrWp8Lk5Nz4xU43NxMhRaa/yZc+bMkd1u1yeffKIxY8YoMjJSUVFRio+P186dOyVJx48f16RJkxQeHq7w8HDdcccdOnHiRLn3feWVVxQTE6OQkBDFxMTo1VdfdToeFBSktWvX6o477lCbNm20YMECt+r/zDPPKDIyUm3bttW0adO0ePFide/e3XHcZrNp6dKl6tq1q1q3bq0BAwbovffeK/ee+/bt06BBgxQaGqqBAwcqNbVur3VKzyQA+KDydmsB3GEotKrpD6NlLEovLiiQjAWpym2ZLLu/uUaeee7cOX300UeyWCxq1qxZqeNBQUGy2+2aNGmSAgICtGXLFhkMBj388MOaNGmSPvnkExkMhlLXbd26VQ8//LAWLlyoIUOG6OOPP9ZDDz2k1q1ba8SIn78jS5Ys0fz585WYmOgoq+jdy/79+2vjxo2SpP/7v//TkiVL9Le//U0DBgzQli1btHz5crVo0cJx/urVq7Vy5Uo988wziomJ0Ztvvqk777xT27dvV48ePUrdPzc3V7fddpuuu+46rV69WpmZmXrssccqbkwfRpgEAB9U3m4tjzzSzku1Ql0WcDHx5yD5E2NRugIuJupS8NoaeebRo0dlt9t1zTXXlHnO9u3b9dVXX2n//v0ym4tD7bp16xQTE6MdO3Zo8ODBpa557rnnNGHCBMXHx0uSOnXqpC+++EIrVqxwCpNjxozRXXfd5XRtRZNcAgICHP/9wgsv6I9//KPjHg8++KBSUlJ05MgRp7rMmjVL48ePlyTNmzdPn376qZ577jm9+OKLpe7/9ttv6/Lly1q1apWaNWumLl266KGHHtK0adPKrZcvI0wCgA9itxZ4ml9RZhnlWTX2TLvdXuE5hw8flslkcgRJSerQoYNMJpMOHTrkMkwePnxYkyZNcirr37+/PvjgA6eymJiYUtdeffXVlay99O2335YKo7169XKEyQsXLigzM1P9+vUrVZeyJhcdPnxYXbt2deqp/e1vf1vpOvmiSr0zuWvXLk2cOFGdO3dWUFCQkpKSnI7n5OTo4YcfVpcuXRQWFqbevXtr1apVNVJhAGgI2K0FnmYzmsooD6uxZ3bs2FEGg0HffvttmefY7XaXQ9mSyiwv69ivy5o2bVrqnLZt25b7a9y4cZWuQ3nKuq4yAbuuqVTPZG5urrp06aLbb79d06dPL3V83rx52r59u1544QWZzWZ9+umneuCBB9SqVStNnDjR45UGgPrOYumt1NTTTkPdJbu1XL580os1Q12V39wiY0Gq01B3kTFC+c0tNfbM4OBg3XjjjVq7dq2mTZtW6r3J7OxsRUdH6/vvv5fVanX0Th47dkyZmZmKjo52ed+oqCjt2bNHd955p6Ns9+7dZZ7/S1UZ5r7mmmu0b98+3XHHHY6yffv2Of47MDBQJpNJe/bs0aBBg5zqEhUV5fL+0dHReuONN5Sbm+sIu//9738rrLcvq1SYHDp0qIYOHSpJmjlzZqnjn332mSZMmKDY2FhJktls1muvvabPP/+cMAkAbihvt5a0NMIkqs7ub1Zuy2QFXEyUX1GWbMYw5Te31NjkmxJPP/20hg0bphtuuEHz5s1T165dZbfblZKSomeffVYHDx5Ut27dFB8fryVLlshut+uRRx5Rz549Hbni1+6//37dc889uvbaazVkyBB99NFHevvtt/Xaa69VWJ+qDHNPnz5d9913n2JiYjRgwAC9++67Sk1NVVBQkFNdFi1apI4dO+raa6/Vm2++qd27d2v79u0u7zlu3Dg9+eSTmjVrlh555BFlZWVp2bJlla6TL/LIO5P9+vXTv/71L911110KDw/X3r179dVXX+lPf/qTJ24PAA1SyW4tgKfY/c01NtmmLB06dNCOHTu0bNkyJSQkKDMzUy1btlS3bt307LPPymAwKCkpSXPnztXIkSMlSYMGDdLSpUvLHCoeOXKkli5dqpUrV+qxxx5Tu3bttGzZMqfJN54wduxYHTt2TE888YQuXbqkkSNH6t5779X777/vOGf69OnKyclRQkKCTp06pcjISP3jH/9wOZNbkpo1a6Y333xTDz74oAYNGqTIyEj99a9/1e233+7RutcmQ3Z2dpUG79u2baulS5c6vfh6+fJl/fnPf1ZSUpL8/Yvz6dKlS3XvvfeWeZ+0tDQ3qwwA1XfixCW98EK6Tp/+USEhV2r69Ai1bdvY29VCHRQZGentKqAWTZo0SYWFhXrzzTe9XRWf4ZGeyTVr1mjv3r1644031K5dO3366ad6/PHH1b59e910000ur6ntL19aWhpf+CqizaqONqs6b7SZ1XpBf/6z8xqOhw/nKzl5RJ3Yq5mfs6qjzeCOvLw8vfTSS7rpppvk7++vLVu26P3339c//vEPb1fNp1Q7TF66dEkLFizQ+vXrHd3L3bp108GDB7Vy5coywySA+sVqLV4DMTMzTybTz+/3+aLy1nBkWBlACYPBoI8++kjPPPOM8vPzdfXVV2vNmjW65ZZbvF01n1LtMFlQUKCCggIZjUancqPRKJvNVt3bA6gDytutxRcDJWs4AqiMxo0ba/Pmzd6uhs+r1DqTOTk5OnDggA4cOCCbzaaMjAwdOHBAx48fV2BgoK677jo98cQTSklJ0bFjx5SUlKR//vOfjhdpAdRv5fX0+SLWcAQAz6lUmNy/f79iY2MVGxurS5cuadGiRYqNjdXChQslSS+//LJiYmIUHx+vfv36afny5Zo3b55jmyMA9Vtd6+mzWHorIsK5x7RkDUcAQNVUaph74MCBys7OLvN4aGionn/+eU/VCUAdU9d6+spbwxH1m6HQ+tM6i5myGU21ss4iUN+xNzeAaitvtxZfVdNrONalCUkNhaHQqqY/jP55B5gCyViQqtyWyQRKoBoIkwCqjZ4+Z3VtQlJD6a0LuJjotJWgJBmL0hVwMbHWF/IG6hPCJACPYLeWn/na0kPlhcWG1FvnV5RZRnlWLdcEqF8IkwDgYb40IamisNiQeutsRpNU4Ko8rPYrA9QjlZrNDQCoPF+akFReWJQaVm9dfnOLiowRTmVFxgjlN7d4qUZA/UCYBAAP89TSQ4ZCqxqfi1PTMyPV+FycDIXWKtelorBoM5pcHq+PvXV2f7NyWybrcsB4FV4xUJcDxtfL4XygtjHMDQAe5okJSYZCq5qeuVlGe0ZxQYFk/HG3cq96r0p1qWhoN7+5RcaCVKfey/rcW2f3N9e74XvA2wiTAFADqjshKeD8oz8HyZ8Y7RkKOP+opAWVvk9FYbGkt654gk6WbMawejubG0DNIEwCgA/yL3C9FWVZ5WWpTFiktw5AdRAmAaCGOJbkKUyXwXZKdkOIbFdcXes9f4RFADWJMAkAHuQIkD/+Pxnth2RQ0S+O/k8q+rxS6zgWNuqjRj++77IcAHwJs7kBwENK1nRslP+2/O1f/ypI/uyXS/OUJT9wkYr8wp3KivzClR+4yGP1BQBPoGcSgEewF7XrNR3LUtE6jnZ/s3JbvVfGu45pHqgtAHgGYRJAtdW1vahrSllrOrpSmXUcedcRQF3AMDeAaitvL+qGpKwFwH+tPq/jCKDhoWcSQLX50l7U3uRqTccSdvmryBAt25WdfW4dR15RAFAdhEkA1eZLe1F7k9OajgXpMthPye4XKpt/B58LkCV4RQFAdREmAVSbxdJbqamnnQKJO3tR1wd17T3H8l5RqM4OPgAaDsIkgGrzxF7Utc2xHmRRpmxGk8/2HNY0XlEAUF2ESQAeUd29qGtTyXqQjncbC1SphcTrI15RAFBdzOYG4BFW6wXFxW3TyJHvKi5um6zWCxVf5CWu1oOszELidVV5fzYWS29FRDj3IDfUVxQAuIeeSQDVVtcmcZS1HmRFC4nXRRX92dTFVxQA+BbCJIBqq2uTOGxGk1TgqrzihcSrwheW3KnMn01dekUBgO8hTAKotro2icPVepCeXkjcV3pr69qfDYC6h3cmAVRbXZvEUbIe5OWA8Sq8YqAuB4z3+OQbX9kVqK792QCoe+iZBFBttbnOpKeW9Knp9SB9pUeQNUAB1DTCJIBqq61JHHVpSR9f6RFkgg2AmkaYBOARtTGJo7wlfXxt1xlf6hFkgg2AmkSYBOCT/PJ3qcn5GTLYsmX3C1Jei9V1akkfegQBNBSESQA+xy9/l5qdGyWDCosLbBfU7NwoFTS63uX5nl7Sx1PoEQTQEDCbG4DPaXJ+xs9B8icGFcp4+bCKjBFO5Z5e0gcAUDX0TALwGkOhVRFXPq6mZ3KcZmYbbNmuz1eeclr+66fZ3FmyGcPcns0NAPAMwiQAr3DMzL4ivXg3ml/OzPYLkmyl9/a2+7Wo8SV9AABVwzA3AK8ob2Z2XovVsv/q37p2+SuvxerarCIAoBIIkwC8oryZ2baA65QTvFlFfu1lUwsV+bVXTvBm2QKuq+VaAgAqwjA3AK+wGU3Fw9ulyotnZtsCrlNOwIFarhUAoKromQTgFd9deEDHvg9xKjv2fYi+u/CAl2oEAHAHPZMAPMKY+39qcuF+GfSj7LpSeYErVdR0bJnnJySe0mefTlHiAx+qTcgFfX86UJYVw/TbAae0tpz5NVbrBSUmpiozM08mEwuBA4C3ESYBVJsx9//U9MIUGX76vUF5anphinKlMgNlZmaerCda6c5H/uhU3j4rr8znWK0XNHr0B05bFKamnlZy8ggCJQB4SaWGuXft2qWJEyeqc+fOCgoKUlJSUqlzjhw5ojvuuEPt27eXyWRSbGysDh8+7PEKA/A9xT2Szgw/lZfFZGrisjwszHW5JCUmpjoFSUlKTy/uqQQAeEelwmRubq66dOmixYsXq3HjxqWOHzt2TMOGDZPZbNaWLVu0e/duWSwWNW3a1OMVBuB7DPqxSuWSZLH0VkSEc29iRESgLJbeZV6Tmem61zKrnN5MAEDNqtQw99ChQzV06FBJ0syZM0sdT0xM1JAhQ/TUU085yjp06OCZGgLweXZdKYNKBzq7rizzGrM5UMnJIzR37ifKzTUqLKzi9x/d6c0EANSsas/mttls+te//qWoqCiNHTtWHTt21A033KBNmzZ5on4AfISh0KrG5+LU9MxINT4XJ0Oh1XEsL3Cl7L863/5TeXnM5kA9+WQXbd06UmvXDqnwvUd3ejNrm9V6QXFx2zRy5LuKi9smq7X0Tj4AUJ8YsrOzf/13QLnatm2rpUuXatKkSZKkkydPKioqSk2aNNFf/vIXxcbG6j//+Y8SEhKUlJSk4cOHu7xPWlpa9WsPoFY0MpzQNY1nKcAvw1GWbwvXt5ee02V7W0lSC+O/FHHlU/IzXJbN3kjpP87T+SLX3//qOHHikl54IV2nT19WSEgjTZ8eobZtS79+4w0nTlzSrFlfKiMj31EWHh6g557r6TN1hOdFRkZ6uwqAV1V7NrfNZpMk/f73v9esWbMkST169NAXX3yhdevWlRkma/vLl5aWxhe+imizqquvbdb43FI1ys9wKgvwy1B0q6Rf7JMdqVz9POGm9U+/KlLVNouMlAYP7lHp82vT0qXbnIKkJGVk5Csp6YzWrh3isefU15+zmkSbATWn2sPcrVq1kr+/v6KiopzKr7nmGmVkZJRxFYDaVN2h1/K2PsTPypoglJ7O0DeA+qvaPZONGjXSb37zm1LD1keOHFG7du2qe3sA1VTZtRkNhVYFXEyUX1GmbEaT8ptbZPc3S5Iu5l+lYGPpe1/Mb8U2Wr9Q1gShb745p9TU047fszYmgPqkUn8P5OTk6MCBAzpw4IBsNpsyMjJ04MABHT9+XJL0pz/9Se+8847Wr1+vo0eP6tVXX9WmTZs0derUGq08gIpVZm1GQ6FVTX8YrUb5b8u/YKca5b+tpj+MdkyysawYpiPWlk73OGJtKcuKYTX/AeoQVxOEmjb1V25uoVMZa2MCqE8qFSb379+v2NhYxcbG6tKlS1q0aJFiY2O1cOFCSdLIkSO1fPlyrVy5UgMGDNCaNWv0wgsvaNgw/qIBvK0yazMGXEyUsSjd6bixKF0BFxMlSQe+aa6bpsTr9a0x2rano17fGqObpsTr4KHmNVfxOqhkuaPx4ztq4ECTxo/vqM6dg12ey9qYAOqLSg1zDxw4UNnZ2eWeM2nSJMcMbwC+ozJrM1b0TqTJ1EQ7d5be+vC3A1jf8dfM5kCnyTZxcduchrhLsDYmgPqC152Aeq4yazPajCaX19qMYZW+B1yj7QDUd9WegAPAt5nNgfowOViBORY1uTJXeT821YVmq9T6F5M/8ptbZCxIdRrqLjJGKL+5xXGP5OQRSkxMVVZWXqV2q0Ex2g5AfUeYBOo5v/xd6njFH+XXsngSSPOmeQqx/1G5+ZtlC7hOkmT3Nyu3ZfJPs7mzZDOGOc3mlkoP39Y0q7V4kkpmZp5MprodwGq77QCgNhEmgXqu0dlp8vNznk3sZyhUo7PTlN/2K0eZ3d/8iwXIvauyyxkBALyPdyaBeq7gxx+qVO4LKrOcEQDANxAmgXruQo7rWcNllfuCyixnBADwDYRJoJ5bnfwnXS5w/qpfLvDT6uQ/ealGFavMckYAAN9AmATqufF336NJj83W0Yxg/XA+QEczgjXpsdkaf/c93q5amVhOBwDqDibgAPWc2RyoxxfO1qOJ1zuWpnl8oW/PjGY5HQCoOwiTQANQF5emqYt1BoCGiGFuAAAAuI0wCQAAALcRJgEAAOA2wiQAAADcRpgEAACA2wiTAAAAcBthEgAAAG4jTAIAAMBthEkAAAC4jTAJAAAAt7GdIgCfZLVeUGJiqjIz82QysTc3APgqwiQAn2O1XtDo0R8oPf2Coyw19bSSk0cQKAHAxzDMDcDnJCamOgVJSUpPL+6pBAD4FsIkAJ+TmZnnsjwry3U5AMB7CJMAfI7J1MRleViY63IAgPcQJgH4HIultyIinN+NjIgIlMXS20s1AgCUhQk4AHyO2Ryo5OQRSkxMVVZWnsLCmM0NAL6KMAnAJ5nNgVq7doi3qwEAqADD3AAAAHAbYRIAAABuI0wCAADAbYRJAAAAuI0wCQAAALcRJgEAAOA2wiQAAADcRpgEAACA2wiTAAAAcBthEgAAAG4jTAIAAMBtlQqTu3bt0sSJE9W5c2cFBQUpKSmpzHMfeOABBQUFaeXKlR6rJAAAAHxTpcJkbm6uunTposWLF6tx48Zlnrd582bt27dPJpPJYxUEAACA76pUmBw6dKjmz5+vUaNGyc/P9SX/+9//9Oijj2rdunXy9/f3aCUBAADgmzyS+goLCzV16lTNmTNHUVFRnrglgFpmtV5QYmKqMjPzZDI1kcXSW2ZzoLerBQDwcR4Jk4sWLVJwcLCmTJlS6WvS0tI88egq8cYz6zrarOrqYpudOHFJs2Z9qYyMfEfZ7t0n9NxzPdW2bdmvtnhKXWwzb6PNqq6m2iwyMrJG7gvUFdUOkzt37tSGDRuUkpJSpetq+8uXlpbGF76KaLOqq6tttnTpNqcgKUkZGflKSjqjtWuH1Oiz62qbeRNtVnW0GVBzqr00UEpKirKyshQVFaVWrVqpVatWOn78uBISEtSlSxdP1BFADcvMzHNZnpXluhwAgBLV7pmcOnWqRo0a5VQ2duxYjR07VnfffXd1bw+gFphMTVyWh4W5LgcAoESlwmROTo6OHj0qSbLZbMrIyNCBAwcUHBysdu3aKSQkxPmm/v4KDQ1lSAGoIyyW3kpNPa309AuOsoiIQFksvb1YKwBAXVCpYe79+/crNjZWsbGxunTpkhYtWqTY2FgtXLiwpusHoBaYzYFKTh6h8eM7auBAk8aP76jk5BHM5gYAVKhSPZMDBw5UdnZ2pW968OBBd+sDwEvM5sAan2wDAKh/2JsbAAAAbiNMAgAAwG2ESQAAALiNMAkAAAC3ESYBAADgNsIkAAAA3EaYBAAAgNsIkwAAAHAbYRIAAABuq9QOOEBlWK0XlJiYqszMPJlMTWSx9GY7PgAA6jnCJDzCar2g0aM/UHr6BUdZauppj+7vbCi0KuBiovyKMmUzmpTf3CK7v9kj9wYAAO5hmBsekZiY6hQkJSk9vbin0hMMhVY1/WG0GuW/Lf+CnWqU/7aa/jBahkKrR+4PAADcQ5iER2Rm5rksz8pyXV5VARcTZSxKdyozFqUr4GKiR+4PAADcQ5iER5hMTVyWh4W5Lq8qv6LMMsqzPHJ/AADgHsIkKsVqvaC4uG0aOfJdxcVtk9XqPKRtsfRWRITzu5EREYGyWHpX6T5lsRlNZZSHVeFTAAAAT2MCDipUmck1ZnOgkpNHKDExVVlZeQoLKz2buzqTdPKbW2QsSHUa6i4yRii/ucVTHxMAALiBMIkKlTe5Zu3aIY4ysznQ6ffu3scVu79ZuS2Tf5rNnSWbMYzZ3AAA+ADCJCrkqck11b2P3d+sS8Frq/RMAABQs3hnEhXy1OSamp6kAwAAah9hEhWq7OSa2roPAADwHQxzo0KVmVxTm/cBAAC+gzCJSqlock1t3wcAAPgGwiTqDKu1eOZ3ZmaeTCZ6NQEA8AWESdQJ1VmjEgAA1Bwm4KBWGQqtanwuTk3PjFTjc3EyFFordV15a1QCAADvoWcStcZQaFXTH0b/vItNgWQsSFVuy+QKFx/31FqXAADAs+iZRK0JuJjotB2iJBmL0hVwMbHCa1mjEgAA30SYRK3xK8osozyrwmtZoxIAAN/EMDdqjc1okgpclYdVeC1rVAIA4JsIk6gUTyzLk9/cImNBqtNQd5ExQvnNLZW6njUqAQDwPYRJVMhTy/LY/c3KbZmsgIuJ8ivKks0Ypvzmlgon3wAAAN9FmESFyluWp6o9hXZ/sy4Fr/Vk9QAAgBcxAQcVYlkeAABQFnomUaHKLstjKLT+NISdKZvRxBA2AAANAGESFbJYeis19bTTUPevl+WpzoLkAACg7mKYGxUqWZZn/PiOGjjQpPHjO5aafFOdBckBAEDdRc8kKqWiZXmqsyA5AACou+iZhEfYjKYyyitekBwAANRdhEl4RH5zi4qMEU5lVVmQHAAA1E2VCpO7du3SxIkT1blzZwUFBSkpKclxrKCgQAkJCRowYIDatGmjqKgoTZ06VcePH6+xSsP3lCxIfjlgvAqvGKjLAeOZfAMAQANQqTCZm5urLl26aPHixWrcuLHTsby8PH355ZeaM2eOduzYoQ0bNujEiRMaN26cCgsLa6TS8E0lC5LnXrVVl4LXEiQBAGgAKjUBZ+jQoRo6dKgkaebMmU7HWrRooeTkZKeyZ599Vv369dPhw4fVtWtXz9QUAAAAPqdGZnNfvHhRkhQUFFTmOWlpaTXx6HJ545l1XUmbNfX7XBFXPiF/w0UV2psr/ccE5dp6ebl2vomfs6qjzaqONqu6mmqzyMjIGrkvUFd4PExevnxZFotFw4cPV9u2bcs8r7a/fGlpaXzhq6ikzfzyd6nZuVkyqPi1BX/lKLrJLOUEb5Yt4Dov19K38HNWdbRZ1dFmVUebATXHo7O5CwsLFR8fr/Pnz+v555/35K3hRU3Oz3AEyRIGFarJ+RleqhEAAPAVHuuZLCws1JQpU/T//t//07vvvquWLVt66tbwMoMtu4zy87VbEQAA4HM80jNZUFCgyZMn6+uvv9bWrVsVGhrqidvCR9j9gsoob1G7FQEAAD6nUj2TOTk5Onr0qCTJZrMpIyNDBw4cUHBwsEwmk+6++27t379fb7zxhgwGg06ePClJCgwMLLWUEOqevBar1ezcKKehbrv8ldditRdrBQAAfEGleib379+v2NhYxcbG6tKlS1q0aJFiY2O1cOFCnThxQu+//74yMzM1ePBgRUVFOX5t2rSppuuPWmALuE45wZtV5NdeNrVQkV97Jt8AAABJleyZHDhwoLKzs8s8Xt4x1A+2gOuUE3DA29UAAAA+hr25AQAA4DbCJAAAANxGmAQAAIDbCJMAAABwG2ESAAAAbiNMAgAAwG2ESQAAALiNMAkAAAC3ESYBAADgNsIkAAAA3EaYBAAAgNsIkwAAAHAbYRIAAABuI0wCAADAbYRJAAAAuI0wCQAAALcRJgEAAOA2wiQAAADcRpgEAACA2wiTAAAAcBthEgAAAG4jTAIAAMBthEkAAAC4jTAJAAAAtxEmAQAA4DbCJAAAANxGmAQAAIDbCJMAAABwG2ESAAAAbiNMAgAAwG2ESQAAALiNMAkAAAC3+Xu7AnBmKLQq4GKi/IoyZTOalN/cIru/2dvVAgAAcIkw6UMMhVY1/WG0jEXpxQUFkrEgVbktkwmUAADAJzHM7UMCLib+HCR/YixKV8DFRC/VCAAAoHyESR9yOS+jSuUAAADeRpj0IfsPun7roKxyAAAAbyNM+pC/vzFGR6wtncqOWFtq5T/HeKlGAAAA5atUmNy1a5cmTpyozp07KygoSElJSU7H7Xa7Fi1apOjoaIWFhenmm2/WN998UyMVrs8MjTropinxen1rjLbt6ajXt8bopinx0hUdvF01AAAAlyoVJnNzc9WlSxctXrxYjRs3LnV8xYoVWrVqlZYsWaJt27YpJCREY8aM0cWLFz1e4frMYuktv0YRuvORP+rGydN15yN/lF+jCFksvb1dNQAAAJcqFSaHDh2q+fPna9SoUfLzc77Ebrdr9erVmj17tkaNGqUuXbpo9erVysnJ0caNG2uk0vWV2Ryo5OQRGj++owYONGn8+I5KTh4hsznQ21UDAABwqdozO6xWq06ePKkhQ4Y4yho3bqwBAwZo7969mjx5cnUf0aCYzYFau3ZIxScCAAD4gGqHyZMnT0qSQkJCnMpDQkKUmZlZ5nVpaWnVfXSVeeOZdR1tVnW0WdXRZlVHm1VdTbVZZGRkjdwXqCs8tuaMwWBw+r3dbi9V9ku1/eVLS0vjC19FtFnV0WZVR5tVHW1WdbQZUHOqvTRQaGioJOnUqVNO5WfOnCnVWwkAAID6pdph0mw2KzQ0VJ988omjLD8/X7t371bfvn2re3sAAAD4sEoNc+fk5Ojo0aOSJJvNpoyMDB04cEDBwcFq166dZsyYoWXLlikyMlKdOnXS008/raZNm2rcuHE1WnkAAAB4V6XC5P79+3XLLbc4fr9o0SItWrRIt99+u1avXq0HHnhAly5d0sMPP6zs7Gz16tVLmzZtUvPmzWus4gAAAPC+SoXJgQMHKjs7u8zjBoNBjz32mB577DFP1QsAAAB1AHtzAwAAwG2ESQAAALiNMAkAAAC3ESYBAADgNsIkAAAA3EaYBAAAgNsIkwAAAHAbYRIAAABuI0wCAADAbYRJAAAAuI0wCQAAALcRJgEAAOA2wiQAAADcRpgEAACA2wiTAAAAcBthEgAAAG4jTAIAAMBthEkAAAC4jTAJAAAAtxEmAQAA4DbCJAAAANxGmAQAAIDbCJMAAABwG2ESAAAAbiNMAgAAwG2ESQAAALiNMAkAAAC3ESYBAADgNsIkAAAA3EaYBAAAgNsIkwAAAHAbYRIAAABuI0wCAADAbYRJAAAAuI0wCQAAALcRJgEAAOA2wiQAAADcRpgEAACA2wiTAAAAcJtHwmRRUZESExPVo0cPhYaGqkePHkpMTFRhYaEnbg8AAAAf5e+Jmyxfvlzr1q3T6tWr1aVLF3399deaMWOGGjVqpEceecQTjwAAAIAP8kiY/OyzzzR8+HCNGDFCkmQ2mzVixAh9/vnnnrg9AAAAfJRHhrn79eunnTt36ttvv5UkHTp0SCkpKfrd737nidsDAADARxmys7Pt1b2J3W5XYmKinnnmGRmNRhUWFmrOnDmyWCxlXpOWllbdxwIA4HWRkZHergLgVR4Z5t60aZP++c9/at26dYqOjtbBgwf16KOPqn379rrrrrtcXlPbX760tDS+8FVEm1UdbVZ1tFnV0WZVR5sBNccjYXL+/PmaNWuWxo4dK0nq2rWrjh8/rmeffbbMMAkAAIC6zyPvTObl5cloNDqVGY1G2Ww2T9weAAAAPsojPZPDhw/X8uXLZTabFR0drQMHDmjVqlWaOHGiJ24PAAAAH+WRMLl06VI99dRTeuihh3TmzBmFhobq7rvvZo1JAACAes4jYbJ58+ZavHixFi9e7InbVYqh0KqAi4nyK8qUzWhSfnOL7P7mWns+AAAAPBQma5uh0KqmP4yWsSi9uKBAMhakKrdlMoESAACgFnlkAk5tC7iY+HOQ/ImxKF0BFxO9VCMAAICGqU6GSb+izDLKs2q5JgAAAA1bnQyTNqOpjPKwWq4JAABAw1Ynw2R+c4uKjBFOZUXGCOU3L3v7RgAAAHhenZyAY/c3K7dl8k+zubNkM4YxmxsAAMAL6mSYlIoD5aXgtd6uBgAAQINWJ4e5AQAA4BsIkwAAAHAbYRIAAABuI0wCAADAbYRJAAAAuI0wCQAAALcRJgEAAOA2wiQAAADcRpgEAACA2wiTAAAAcBthEgAAAG4jTAIAAMBthEkAAAC4jTAJAAAAtxEmAQAA4DbCJAAAANxGmAQAAIDbCJMAAABwG2ESAAAAbiNMAgAAwG2ESQAAALiNMAkAAAC3+Xu7Au6yWi8oMTFVmZl5MpmayGLpLbM50NvVAgAAaFDqZJi0Wi9o9OgPlJ5+wVGWmnpayckjCJQAAAC1qE4OcycmpjoFSUlKTy/uqQQAAEDtqZNhMjMzz2V5VpbrcgAAANSMOhkmTaYmLsvDwlyXAwAAoGbUyTBpsfRWRITzu5EREYGyWHp7qUYAAAANU52cgGM2Byo5eYQSE1OVlZWnsDBmcwMAAHhDnQyTUnGgXLt2iLerAQAA0KDVyWFuAAAA+AbCJAAAANzmsTCZlZWl6dOnq2PHjgoNDVXfvn21c+dOT90eAAAAPsgj70xmZ2dr2LBh6tevn9566y21atVKVqtVISEhnrg9AAAAfJRHwuTf//53hYWFac2aNY6yDh06eOLWAAAA8GGG7Oxse3Vv0rdvX914443KzMxUSkqKwsLCdNdddykuLk4Gg8HlNWlpadV9LAAAXhcZGentKgBe5ZEwGRoaKkmaOXOmRo8erYMHD2ru3LlKSEhQfHx8tSvpCWlpaXzhq4g2qzrarOpos6qjzaqONgNqjkeGuW02m2JiYpSQkCBJ6tmzp44ePap169b5TJgEAACA53lkNndoaKiioqKcyq655hplZGR44vYAAADwUR4Jk/369dORI0ecyo4cOaJ27dp54vYAAADwUR55Z3Lfvn0aOnSoHn30Ud166606cOCA/vSnP+nxxx9XXFycJ+oJAAAAH+SRMClJH374oRYsWKAjR44oPDxccXFxmjZtWpmzuQEAAFD3eSxMAgAAoOFhb24AAAC4jTAJAAAAtxEmAQAA4DbCJAAAANxWb8PksmXLFBQUpIcffthRZrfbtWjRIkVHRyssLEw333yzvvnmGy/W0vuysrI0ffp0dezYUaGhoerbt6927tzpOE6bOSsqKlJiYqJ69Oih0NBQ9ejRQ4mJiSosLHSc09DbbNeuXZo4caI6d+6soKAgJSUlOR2vTPv8+OOPevjhh3X11VerTZs2mjhxok6cOFGbH6NWlddmBQUFSkhI0IABA9SmTRtFRUVp6tSpOn78uNM9aLOkMs994IEHFBQUpJUrVzqVN7Q2A2pKvQyT//3vf/Xqq6+qa9euTuUrVqzQqlWrtGTJEm3btk0hISEaM2aMLl686KWaeld2draGDRsmu92ut956S3v37tXSpUsVEhLiOIc2c7Z8+XKtW7dOS5Ys0WeffabFixdr7dq1euaZZxznNPQ2y83NVZcuXbR48WI1bty41PHKtM9jjz2mrVu36qWXXtL777+vixcvasKECSoqKqrNj1JrymuzvLw8ffnll5ozZ4527NihDRs26MSJExo3bpzTP2JoM9c2b96sffv2yWQylTrW0NoMqCn1bmmg8+fPa9CgQVqxYoWWLl2qLl266G9/+5vsdruio6MVFxenOXPmSJIuXbqkyMhIPfnkk5o8ebKXa177FixYoF27dunDDz90eZw2K23ChAkKDg7WCy+84CibPn26zp07pzfffJM2+5W2bdtq6dKlmjRpkqTK/UydP39enTp10qpVq3TbbbdJkjIyMtS9e3dt3LhRN954o9c+T234dZu5cujQIfXr10+7du1S165dabMy2ux///ufhg0bpuTkZI0bN07x8fG6//77JanBtxngSfWuZ3L27NkaNWqUBg0a5FRutVp18uRJDRkyxFHWuHFjDRgwQHv37q3tavqE9957T7169dLkyZPVqVMnXX/99XrxxRdltxf/+4I2K61fv37auXOnvv32W0nFf6mnpKTod7/7nSTarCKVaZ8vvvhCBQUFTueEh4crKiqKNvxJSS9uUFCQJNrMlcLCQk2dOlVz5sxRVFRUqeO0GeA5/t6ugCe9+uqrOnr0qNasWVPq2MmTJyXJaQi35PeZmZm1Uj9fc+zYMb300kuaOXOmZs+erYMHD2ru3LmSpPj4eNrMhdmzZysnJ0d9+/aV0WhUYWGh5syZo6lTp0ri56wilWmfU6dOyWg0qlWrVqXOOXXqVO1U1IddvnxZFotFw4cPV9u2bSXRZq4sWrRIwcHBmjJlisvjtBngOfUmTKalpWnBggX64IMP1KhRozLP+/X2jna7vcFu+Wiz2RQTE6OEhARJUs+ePXX06FGtW7dO8fHxjvNos59t2rRJ//znP7Vu3TpFR0fr4MGDevTRR9W+fXvdddddjvNos/K50z60YXFvW3x8vM6fP6833nijwvMbapvt3LlTGzZsUEpKSpWvbahtBlRHvRnm/uyzz3T27Fn1799frVq1UqtWrbRr1y6tW7dOrVq1UsuWLSWp1L84z5w5U6qXpKEIDQ0tNfxzzTXXKCMjw3Fcos1+af78+Zo1a5bGjh2rrl27auLEibrvvvv07LPPSqLNKlKZ9mndurWKiop09uzZMs9piAoLCzVlyhR9/fXX2rx5s+P/aRJt9mspKSnKyspSVFSU4++D48ePKyEhQV26dJFEmwGeVG/C5M0336xPP/1UKSkpjl8xMTEaO3asUlJS1KlTJ4WGhuqTTz5xXJOfn6/du3erb9++Xqy59/Tr109HjhxxKjty5IjatWsnSTKbzbTZr+Tl5cloNDqVGY1G2Ww2SbRZRSrTPtdee62uuOIKp3NOnDihw4cPN9g2LCgo0OTJk/X1119r69atjlBegjZzNnXqVO3atcvp7wOTyaSZM2dq8+bNkmgzwJPqzTB3UFCQ42X0Ek2aNFFwcLDjX6IzZszQsmXLFBkZqU6dOunpp59W06ZNNW7cOC/U2PtmzpypoUOH6umnn9att96qAwcO6MUXX9Tjjz8uqXgokjZzNnz4cC1fvlxms1nR0dE6cOCAVq1apYkTJ0qizSQpJydHR48elVT8KkVGRoYOHDig4OBgtWvXrsL2adGihe68807Nnz9fISEhCg4O1rx589S1a1cNHjzYi5+s5pTXZiaTSXfffbf279+vN954QwaDwfHuaWBgoBo3bkybufg5+3Xvor+/v0JDQxUZGSmpYf6cATWl3i0N9Es333yzY2kgqfhdmMWLF2v9+vXKzs5Wr1699PTTTzvCZkP04YcfasGCBTpy5IjCw8MVFxenadOmOd4Zos2cXbx4UU899ZTeffddnTlzRqGhoRo7dqweeeQRBQQESKLNUlJSdMstt5Qqv/3227V69epKtU9+fr4ef/xxbdy4Ufn5+YqNjdWyZcsUHh5emx+l1pTXZo8++qh69uzp8rpVq1Y5lsOhzYqV/Jz9Wvfu3Z2WBpIaXpsBNaVeh0kAAADUrHrzziQAAABqH2ESAAAAbiNMAgAAwG2ESQAAALiNMAkAAAC3ESYBAADgNsIkUEuSkpIUFBQkq9Va5WutVquCgoIc2zbWtpLnJyUlOZV/8cUXGjFihMLDwxUUFOTWXshVMWPGDHXv3r1GnwEAqBrCJBq0LVu2KCgoSBs3bix17JZbbin3mNlslt3ue8u07t69W4sWLVJ2dnalzi8Juf/973+r9JyioiJNnjxZ33//vRYsWKA1a9YoKipKb775pp5//nk3ag4AqIsIk2jQ+vfvL6k4gP1SYWGhPv/8c/n7+5d5rF+/fo6dgipj4sSJysrKUvv27atf8XLs2bNHS5Ys0fnz5z12z/bt2ysrK8uxbaQkZWRkKD09XdOmTdO9996rCRMmqHXr1nrrrbdc7kACAKif6s3e3IA7QkJC1LFjx1KB8csvv1ReXp5uu+22Mo/169evSs8yGo0yGo3VrrM3GAwGx3aRJc6cOSOpeI9jAEDDRc8kGrz+/fvr0KFDTsPCe/bskclk0oQJE1weK7muxCeffKKRI0cqPDxcbdq00ciRI7V3716n55T1zuT69esVExOj0NBQXX/99frXv/5V7ruBb7zxhvr06aPWrVtrwIAB2r59u+PYokWL9MQTT0iSevbsqaCgII+8y/jrdyZnzJihG2+8UZJ03333KSgoSN27d9fNN9+sjz/+WMePH3c8OygoyHEfu92uF198UQMGDFBoaKgiIiIUFxenEydOlHrm66+/rl69eik0NFTXXXedPvjgg2p9BgBAzaBnEg1ev3799Prrr+uzzz7T0KFDJRUHxr59+6pPnz6SVOpYQECAYmJiJEkbN25UfHy8Bg4cqHnz5slmsykpKUl/+MMf9N5776l3795lPnv9+vWaPXu2+vTpo/j4eJ05c0bTpk1T27ZtXZ6/efNmnT17VpMnT1ZAQIBWr16tO+64QwcPHlRwcLBuueUWpaWladOmTVq4cKFatWolSYqKivJYe0nS5MmTZTabtXjxYt1zzz3q37+/mjZtqqZNmyo7O1tZWVlauHBhqesefPBB/eMf/9CECRM0depUnTx5Ui+++KL27t2r//znP47guWHDBs2aNUu/+c1vNHXqVJ0+fVrTpk1TeHi4Rz8HAKD6CJNo8Ep6GPfs2eMIjHv37tWf//xnBQYGKjo6utSxmJgYXXnllcrNzdWcOXM0YcIEp/cEJ0+erH79+mnBggXasmWLy+cWFBToySefVLdu3fTee++pUaNGkqTY2FiNGjVK7dq1K3VNenq6Pv/8c1111VWSpOuvv16xsbHauHGj4uLi1K1bN3Xv3l2bNm3SzTffLLPZ7LmG+oXf/va3MhgMWrx4sfr06aMJEyY4joWFhenChQtOZVJxu73yyitatWqVJk2a5Ci/5ZZbNHjwYL344ot65JFHVFhYqL/+9a+Kjo7W+++/7xhev/7663Xrrbe6bBcAgPcwzI0Gr2PHjgoNDXW8G/ndd9/p1KlTjnci+/XrV+rYgAEDJBUPb2dnZ+u2227T2bNnHb8uXbqkwYMHa/fu3SooKHD53H379uns2bO65557HEFSkgYNGqTOnTu7vGb06NGOIClJPXr0UGBgoI4dO1btdqhp77zzjpo1a6ahQ4c6tZXJZFLHjh31n//8R1Jxu5w6dcrR+1piyJAhio6O9lb1AQBloGcSkNS3b1/9+9//1uXLl7Vnzx41adLE8c5i3759tWHDBscxSY6g+d1330mSxowZU+a9z58/7xQASxw/flxScZj9tY4dO+rLL78sVe6qV65FixY6d+5cRR/R67777jvl5OQoMjLS5fGSmfEl7eLqvE6dOrlsFwCA9xAmARWHwy1btmj//v3as2ePevXqJX//4q9H3759lZ+f7zjm5+en3/72t5Ikm80mSXr++efVpk0bl/cODAyscn3KWr+yrNngvrje5a/ZbDa1bNlSL7/8ssvjTZo0kfTzZ3G17FJd+JwA0NAQJgHJMWy9Z88e7dmzR6NHj3Yc69Chg8LCwhzHunbt6lgOJyIiQpJ01VVXafDgwVV6Zkkv43fffacbbrjB6djRo0fd/CSuQ1htKuv5ERER+uSTT9SrVy81b968zOtL1uH89ttvS7VLSU8wAMB38M4kIKl79+5q1qyZ3nvvPaWlpZVaQ7Jv374uj914441q0aKFnn76af3444+l7luyFqMrMTExatWqldavX6/Lly87ynfs2KFvvvnG7c9S0sNX2R1wPK1JkyYuF0y/9dZbZbPZtHjx4lLH7Ha7zp49K6m4XUJCQrR+/Xrl5+c7ztm2bZsOHTpUcxUHALiFnklAxcPHvXv31vbt2+Xn51dqOZ++ffvqL3/5i6SfezElqXnz5lqxYoWmTJmi66+/XuPHj1doaKhOnDihlJQUNW3a1OV2jJLUqFEjzZs3Tw8++KBuvvlmjR07VmfOnNHatWvVpUsX5eTkuPVZSpYsevLJJzV27Fg1atRIsbGxCgkJKfe6DRs2OK1ZWeLuu++u8vO3bNmiuXPnqnfv3vLz89PYsWM1YMAATZs2TatWrdJXX32lm266SU2aNJHVatW7776rO++8U3/+8591xRVXaP78+br//vv1+9//XuPHj3e0S+fOnd1uFwBAzSBMAj/p37+/tm/frs6dO5fa1eWXvZG/7rUcPXq0TCaTnnnmGT3//PO6dOmSQkND1bt3b911113lPvPee++VJP39739XQkKCIiMjtWbNGm3YsMHtXrg+ffrIYrFo/fr1uu+++2Sz2bR169YKw+Qrr7zisnzYsGFV2uUmPj5ehw4d0ltvvaUXX3xRdrtdY8eOlSQtWbJE1157rV566SUtWrRIfn5+atOmjW688UaNHDnScY8777xTdrtdy5cvV0JCgjp16qQ1a9Zoy5Yt2rlzZ6XrAgCoeYbs7GzeaAd8zHXXXaeQkBAlJyd7uyoAAJSLdyYBL8rPzy81Q3nHjh36+uuvFRsb66VaAQBQefRMAl6UkpKiuXPn6g9/+IPCwsL0zTffaP369WrVqpU+/fRTp32tAQDwRbwzCXhR+/btZTab9corr+iHH35QYGCgRo4cqfnz5xMkAQB1Aj2TAAAAcBvvTAIAAMBthEkAAAC4jTAJAAAAtxEmAQAA4DbCJAAAANz2/wH00GtSe2hzdwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "shotput['Best Quadratic Curve'] = shotput_fit\n", "\n", "fig, ax = plt.subplots(figsize=(7,6))\n", "\n", "ax.scatter(shotput['Weight Lifted'], \n", " shotput['Shot Put Distance'], \n", " label='Color=darkblue', \n", " color='darkblue')\n", "\n", "ax.scatter(shotput['Weight Lifted'], \n", " shotput['Best Quadratic Curve'], \n", " label='Color=gold', \n", " color='gold')\n", "\n", "x_label = 'Weight Lifted'\n", "\n", "y_label = ''\n", "\n", "y_vals = ax.get_yticks()\n", "\n", "plt.ylabel(y_label)\n", "\n", "ax.legend(bbox_to_anchor=(1.04,1), loc=\"upper left\", frameon=False)\n", "\n", "plt.xlabel(x_label)\n", "\n", "plt.show()" ] } ], "metadata": { "anaconda-cloud": {}, "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": 1 }