commit e93e41c96fa518e78b3c2c23d373763525b6de84 Author: Gerardo Marx Date: Sat Feb 17 20:20:13 2024 -0600 Script done diff --git a/main.ipynb b/main.ipynb new file mode 100644 index 0000000..ef87ae9 --- /dev/null +++ b/main.ipynb @@ -0,0 +1,410 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# OLS Ordinary Least Squares\n", + "\n", + "The OLS general model $\\hat{y}$ is defined by: \n", + "\n", + "$$ \\hat{y} = \\theta_0+\\theta_1 x_1 $$\n", + "\n", + "Applying the partial derivatives with rescpect $\\theta_0$ and equaliting to zero:\n", + "\n", + "$$\\frac{\\partial SSR}{\\partial \\theta_0}=0 $$\n", + "\n", + "here SSR is defined as:\n", + "\n", + "$$ \\sum_{i=1}^n (y^i - \\hat{y}^i)^2 $$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Resulting in:\n", + "\n", + "$$ \\theta_0 = \\frac{\\sum_{i=1}^n y^i}{n} - \\frac{\\theta_1 \\sum_{i=1}^n x^i}{n}$$\n", + "\n", + "or \n", + "\n", + "$$ \\theta_0 = \\bar{y} -\\theta_1 \\bar{x} $$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In a similar way, the partial derivative of SSR with respect of $\\theta_1$ will result in: \n", + "\n", + "$$\\theta_1 = \\frac{\\sum_{i=1}^n x^i(y^i-\\bar{y}) }{\\sum_{i=1}^n x^i(x^i-\\bar{x})}$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Implementing OLS in Python" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 3.9654 , 4.50131579, 5.03723158, 5.57314737, 6.10906316,\n", + " 6.64497895, 7.18089474, 7.71681053, 8.25272632, 8.78864211,\n", + " 9.32455789, 9.86047368, 10.39638947, 10.93230526, 11.46822105,\n", + " 12.00413684, 12.54005263, 13.07596842, 13.61188421, 14.1478 ])" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import numpy as np\n", + "x = np.linspace(0,4,20)\n", + "theta0 = 3.9654\n", + "theta1 = 2.5456\n", + "y = theta0+theta1*x\n", + "y" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAj5klEQVR4nO3df4wcZ33H8c/uXO6cBt8Rp8T2ac4TwywGnMQQQiInlWJv3Fqu5cZ/8CtyjQWubCrT2AQBtkQIUaCXVAgCxUp26RCn5UdKARvEr9QN67iQBGyfrzUUwg54zT0Ex6pEb20Hjmh3+sfVq659d7497z6zs/d+SfPHzT17833yZP18NM/sPqkoiiIBAABYko67AAAAMLsQPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABY1RV3AeerVqt6/vnnNXfuXKVSqbjLAQAA0xBFkU6fPq3+/n6l01Pf22i78PH8889rYGAg7jIAAMAMjIyMyHXdKdu0XfiYO3eupPHie3t7Y64GAABMR7lc1sDAQG0en0rbhY9zSy29vb2EDwAAEmY6j0zwwCkAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAADCLGGNUKBRkjImtBsIHAACzRBAE8jxP2WxWnucpCIJY6khFURTFcuVJlMtl9fX1aXR0lI3lAABoEmOMPM9TtVqtnXMcR6VSSa7rXvLfb2T+5s4HAACzQLFYrAseklSpVBSGofVaGg4fBw8e1Lp169Tf369UKqV9+/ZN2vbd7363UqmUHnrooUsoEQAAXKpMJqN0un7adxxHvu9br6Xh8HH27FktW7ZMu3fvnrLd3r179eyzz6q/v3/GxQEAgOZwXVf5fF6O40gaDx65XK4pSy6N6mr0BWvWrNGaNWumbPPrX/9af/M3f6MnnnhCa9eunXFxAACgeTZv3qzVq1crDEP5vh9L8JBmED4uplqtauPGjXr/+9+vpUuXXrT92NiYxsbGaj+Xy+VmlwQAAP6P67qxhY5zmv7A6YMPPqiuri7ddddd02o/ODiovr6+2jEwMNDskgAAQBtpavg4cuSIPvWpT2nPnj1KpVLTes2uXbs0OjpaO0ZGRppZEgAAaDNNDR///u//rlOnTmnRokXq6upSV1eXTpw4ofe973265pprJnxNT0+Pent76w4AANC5mvrMx8aNG7Vq1aq6c6tXr9bGjRv1zne+s5mXAgAACdVw+Dhz5kzdF5IcP35cw8PDmjdvnhYtWqSrrrqqrv1ll12mBQsWaMmSJZdeLQAASLyGw8fhw4e1cuXK2s933323JGnTpk3as2dP0woDAACdqeHwsWLFCjWyHUypVGr0EgAAoIOxtwsAALCK8AEAAKwifAAAAKsIHwAAtBFjjAqFgowxcZfSMoQPAADaRBAE8jxP2WxWnucpCIK4S2qJVNTIR1csKJfL6uvr0+joKN92CgCYNYwx8jxP1Wq1ds5xHJVKpdg3gpuORuZv7nwAANAGisViXfCQpEqlUvfFnp2C8AEAQBvIZDJKp+unZcdx5Pt+TBW1DuEDAIA24Lqu8vm8HMeRNB48crlcIpZcGsUzHwAAtBFjjMIwlO/7iQoejczfTd3VFgAAXBrXdRMVOmaCZRcAAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAKABxhgVCgUZY+IuJbEIHwAATFMQBPI8T9lsVp7nKQiCuEtKJHa1BQBgGowx8jxP1Wq1ds5xHJVKpY7fCG46Gpm/ufMBAMA0FIvFuuAhSZVKRWEYxlRRchE+AACYhkwmo3S6ftp0HEe+78dUUXIRPgAAmAbXdZXP5+U4jqTx4JHL5VhymQGe+QAAoAHGGIVhKN/3CR7/TyPzd5elmgAA6Aiu6xI6LhHLLgAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEA6DjGGBUKBRlj4i4FEyB8AAA6ShAE8jxP2WxWnucpCIK4S8J52NUWANAxjDHyPE/VarV2znEclUolNoNrsUbmb+58AAA6RrFYrAseklSpVBSGYUwVYSKEDwBAx8hkMkqn66c2x3Hk+35MFWEihA8AQMdwXVf5fF6O40gaDx65XI4llzbDMx8AgI5jjFEYhvJ9n+BhSSPzd5elmgAAsMZ1XUJHG2PZBQAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGBVw+Hj4MGDWrdunfr7+5VKpbRv377a71566SV98IMf1HXXXacrrrhC/f39esc73qHnn3++mTUDAIAEazh8nD17VsuWLdPu3bsv+N2LL76ooaEh3XPPPRoaGtLXvvY1Pffcc/qLv/iLphQLAOgMxhgVCgUZY+IuBTG4pL1dUqmU9u7dq/Xr10/a5tChQ7rpppt04sQJLVq06KJ/k71dAKCzBUGgLVu2qFqtKp1OK5/Pa/PmzXGXhUvUyPzd8mc+RkdHlUql9PKXv3zC34+NjalcLtcdAIDOZIypBQ9Jqlar2rp1K3dAZpmWho/f//73+uAHP6g777xz0hQ0ODiovr6+2jEwMNDKkgAAMSoWi7XgcU6lUlEYhjFVhDi0LHy89NJLeutb36ooivTwww9P2m7Xrl0aHR2tHSMjI60qCQAQs0wmo3S6fupxHEe+78dUEeLQkvBxLnicOHFC+/fvn3Ltp6enR729vXUHAKAzua6rfD4vx3EkjQePXC4n13Vjrgw2dTX7D54LHsViUYVCQVdddVWzLwEASLDNmzdr9erVCsNQvu8TPGahhsPHmTNn6tbmjh8/ruHhYc2bN08LFy7Um9/8Zg0NDemb3/ymKpWKTp48KUmaN2+euru7m1c5ACCxXNcldMxiDX/U9sCBA1q5cuUF5zdt2qSPfOQjWrx48YSvKxQKWrFixUX/Ph+1BQAgeRqZvxu+87FixQpNlVcu4WtDAADALMDeLgAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AgAsYY1QoFGSMibsUdCDCBwCgThAE8jxP2WxWnucpCIK4S0KHaXhvl1ZjbxcAiI8xRp7nqVqt1s45jqNSqcRGcJhSI/M3dz4AADXFYrEueEhSpVKp280cuFSEDwBATSaTUTpdPzU4jiPf92OqCJ2I8AEAqHFdV/l8Xo7jSBoPHrlcjiUXNBXPfAAALmCMURiG8n2f4IFpaWT+7rJUEwAgQVzXJXSgZVh2AQAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwBIIGOMCoWCjDFxlwI0jPABAAkTBIE8z1M2m5XneQqCIO6SgIawqy0AJIgxRp7nqVqt1s45jqNSqcRGcIhVI/M3dz4AIEGKxWJd8JCkSqWiMAxjqghoHOEDABIkk8kona7/p9txHPm+H1NFQOMIHwCQIK7rKp/Py3EcSePBI5fLseSCROGZDwBIIGOMwjCU7/sED7SFRubvLks1AQCayHVdQgcSi2UXAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwBaxBijQqEgY0zcpQBthfABAC0QBIE8z1M2m5XneQqCIO6SgLbBxnIA0GTGGHmep2q1WjvnOI5KpRL7saBjNTJ/c+cDAJqsWCzWBQ9JqlQqCsMwpoqA9kL4AIAmy2QySqfr/3l1HEe+78dUEdBeCB8A0GSu6yqfz8txHEnjwSOXy7HkAvwfnvkAgBYxxigMQ/m+T/BAx2tk/u6yVBMAzDqu6xI6gAmw7AIAAKwifAAAAKsIHwAAwCrCBwAAsKrh8HHw4EGtW7dO/f39SqVS2rdvX93voyjShz/8YS1cuFCXX365Vq1apWKx2Kx6AQBAwjUcPs6ePatly5Zp9+7dE/7+7/7u7/TpT39ajzzyiH74wx/qiiuu0OrVq/X73//+kosFAADJ1/BHbdesWaM1a9ZM+LsoivTQQw/pQx/6kO644w5J0j/+4z9q/vz52rdvn97+9rdfWrUAACDxmvrMx/Hjx3Xy5EmtWrWqdq6vr08333yznnnmmQlfMzY2pnK5XHcAAIDO1dTwcfLkSUnS/Pnz687Pnz+/9rvzDQ4Oqq+vr3YMDAw0syQAANBmYv+0y65duzQ6Olo7RkZG4i4JAAC0UFPDx4IFCyRJL7zwQt35F154ofa78/X09Ki3t7fuAAAbjDEqFAoyxsRdCjCrNDV8LF68WAsWLNCTTz5ZO1cul/XDH/5Qy5cvb+alAOCSBEEgz/OUzWbleZ6CIIi7JGDWaDh8nDlzRsPDwxoeHpY0/pDp8PCwfvWrXymVSmnHjh366Ec/qm984xs6duyY3vGOd6i/v1/r169vcukAMDPGGG3ZskXValWSVK1WtXXrVu6AAJY0/FHbw4cPa+XKlbWf7777bknSpk2btGfPHn3gAx/Q2bNntWXLFv3P//yP/uRP/kTf/e53NWfOnOZVDQCXoFgs1oLHOZVKRWEYsgstYEEqiqIo7iL+v3K5rL6+Po2OjvL8B4CWMMbI87y6AOI4jkqlEuEDmKFG5u/YP+0CALa5rqt8Pi/HcSSNB49cLkfwACzhzgeAWcsYozAM5fs+wQO4RI3M3w0/8wEAncJ1XUIHEAOWXQAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QNA2zLGqFAoyBgTdykAmojwAaAtBUEgz/OUzWbleZ6CIIi7JABNwq62ANqOMUae56lardbOOY6jUqnERnBAm2pk/ubOB4C2UywW64KHJFUqFYVhGFNFAJqJ8AGg7WQyGaXT9f88OY4j3/djqghAMxE+ALQd13WVz+flOI6k8eCRy+VYcgE6BM98AGhbxhiFYSjf9wkeQJtrZP7uslQTADTMdV1CB9CBWHYBAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AMyYMUaFQkHGmLhLAZAghA8AMxIEgTzPUzabled5CoIg7pIAJAQbywFomDFGnuepWq3WzjmOo1KpxF4swCzVyPzNnQ8ADSsWi3XBQ5IqlYrCMIypIgBJQvgA0LBMJqN0uv6fD8dx5Pt+TBUBSBLCB4CGua6rfD4vx3EkjQePXC7HkguAaeGZDwAzZoxRGIbyfZ/gAcxyjczfXZZqAtCBXNcldABoGMsuAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAXQ4Y4wKhYKMMXGXAgCSCB9ARwuCQJ7nKZvNyvM8BUEQd0kAwK62QKcyxsjzPFWr1do5x3FUKpXYDA5A0zUyf3PnA+hQxWKxLnhIUqVSURiGMVUEAOMIH0CHymQySqfr3+KO48j3/ZgqAoBxhA+gQ7muq3w+L8dxJI0Hj1wux5ILgNjxzAfQ4YwxCsNQvu8TPAC0TKzPfFQqFd1zzz1avHixLr/8cr3qVa/S/fffrzbLOMCs4bquVqxYQfAA0Da6mv0HH3zwQT388MN67LHHtHTpUh0+fFjvfOc71dfXp7vuuqvZlwMAAAnT9PDx9NNP64477tDatWslSddcc42+9KUv6Uc/+lGzLwUAABKo6csut9xyi5588kn9/Oc/lyT9x3/8h77//e9rzZo1E7YfGxtTuVyuOwAAQOdq+p2PnTt3qlwu6zWveY0cx1GlUtHHPvYxbdiwYcL2g4ODuu+++5pdBgAAaFNNv/Px5S9/WV/4whf0xS9+UUNDQ3rsscf08Y9/XI899tiE7Xft2qXR0dHaMTIy0uySAABAG2n6R20HBga0c+dObdu2rXbuox/9qD7/+c/rZz/72UVfz0dtAQBInlg/avviiy9O+K2K53/NMwAAmJ2a/szHunXr9LGPfUyLFi3S0qVLdfToUX3iE5/Qu971rmZfCgAAJFDTl11Onz6te+65R3v37tWpU6fU39+vO++8Ux/+8IfV3d190dez7AIAQPI0Mn/z9eoAAOCSxfrMB4DGGGNUKBRkjIm7FACwgvABxCgIAnmep2w2K8/zFARB3CUBQMux7ALExBgjz/PqPgnmOI5KpRKbwAFIHJZdgAQoFosXfAS9UqkoDMOYKgIAOwgfQEwymcyE34nj+35MFQGAHYQPICau6yqfz8txHEnjwSOXy7HkAqDj8cwHEDNjjMIwlO/7BA8AidXI/N30bzgF0BjXdQkdAGYVll0AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQP4CKMMSoUCjLGxF0KAHQEwgcwhSAI5HmestmsPM9TEARxlwQAicfGcsAkjDHyPE/VarV2znEclUol9mIBgPM0Mn9z5wOYRLFYrAseklSpVBSGYUwVAUBnIHwAk8hkMkqn698ijuPI9/2YKgKAzkD4ACbhuq7y+bwcx5E0HjxyuRxLLgBwiXjmA7gIY4zCMJTv+wQPAJhEI/N3l6WagMRyXZfQAQBNxLILAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8oCMYY1QoFGSMibsUAMBFED6QeEEQyPM8ZbNZeZ6nIAjiLgkAMAV2tUWiGWPkeZ6q1WrtnOM4KpVKbAYHABY1Mn9z5wOJViwW64KHJFUqFYVhGFNFAICLIXwg0TKZjNLp+v+NHceR7/sxVQQAuBjCBxLNdV3l83k5jiNpPHjkcjmWXACgjfHMBzqCMUZhGMr3fYIHAMSgkfm7y1JNQEu5rkvoAICEYNkFAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+YIUxRoVCQcaYuEsBAMSM8IGWC4JAnucpm83K8zwFQRB3SQCAGLUkfPz617/WX/7lX+qqq67S5Zdfruuuu06HDx9uxaXQ5owx2rJli6rVqiSpWq1q69at3AEBgFms6bva/va3v9Wtt96qlStX6jvf+Y5e8YpXqFgs6sorr2z2pZAAxWKxFjzOqVQqCsOQXWgBYJZqevh48MEHNTAwoEcffbR2bvHixc2+DBIik8konU7XBRDHceT7foxVAQDi1PRll2984xu68cYb9Za3vEVXX3213vCGN+izn/3spO3HxsZULpfrDnQO13WVz+flOI6k8eCRy+W46wEAs1gqiqKomX9wzpw5kqS7775bb3nLW3To0CFt375djzzyiDZt2nRB+4985CO67777Ljg/Ojqq3t7eZpaGGBljFIahfN8neABAByqXy+rr65vW/N308NHd3a0bb7xRTz/9dO3cXXfdpUOHDumZZ565oP3Y2JjGxsZqP5fLZQ0MDBA+AABIkEbCR9OXXRYuXKjXve51dede+9rX6le/+tWE7Xt6etTb21t3AACAztX08HHrrbfqueeeqzv385//XJ7nNftSAAAggZoePt773vfq2Wef1d/+7d8qDEN98YtfVD6f17Zt25p9KQAAkEBNDx9vetObtHfvXn3pS1/Stddeq/vvv18PPfSQNmzY0OxLAQCABGr6A6eXqpEHVgAAQHuI9YFTAACAqRA+AACAVYQPAABgFeEDAABYRfiApPGvPy8UCmx1DwBoOcIHFASBPM9TNpuV53kKgiDukgAAHYyP2s5yxhh5nnfBlvelUokN4AAA08ZHbTFtxWKxLnhIUqVSURiGMVUEAOh0hI9ZLpPJKJ2u/9/AcRz5vh9TRQCATkf4mOVc11U+n5fjOJLGg0cul2PJBQDQMjzzAUnjz36EYSjf9wkeAICGNTJ/d1mqCW3OdV1CBwDACpZdAACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhIyGMMSoUCjLGxF0KAACXhPCRAEEQyPM8ZbNZeZ6nIAjiLgkAgBljV9s2Z4yR53mqVqu1c47jqFQqsREcAKBtNDJ/c+ejzRWLxbrgIUmVSkVhGMZUEQAAl4bw0eYymYzS6fphchxHvu/HVBEAAJeG8NHmXNdVPp+X4ziSxoNHLpdjyQUAkFg885EQxhiFYSjf9wkeAIC208j83WWpJlwi13UJHQCAjsCyCwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifDSRMUaFQkHGmLhLAQCgbRE+miQIAnmep2w2K8/zFARB3CUBANCW2NW2CYwx8jxP1Wq1ds5xHJVKJTaDAwDMCo3M39z5aIJisVgXPCSpUqkoDMOYKgIAoH0RPpogk8kona7/T+k4jnzfj6kiAADaF+GjCVzXVT6fl+M4ksaDRy6XY8kFAIAJ8MxHExljFIahfN8neAAAZpVG5u8uSzXNCq7rEjoAALgIll0AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWNXy8PHAAw8olUppx44drb4UAABIgJaGj0OHDimXy+n6669v5WUAAECCtCx8nDlzRhs2bNBnP/tZXXnlla26TEOMMSoUCjLGxF0KAACzVsvCx7Zt27R27VqtWrVqynZjY2Mql8t1RysEQSDP85TNZuV5noIgaMl1AADA1FoSPh5//HENDQ1pcHDwom0HBwfV19dXOwYGBppejzFGW7ZsUbValSRVq1Vt3bqVOyAAAMSg6eFjZGRE27dv1xe+8AXNmTPnou137dql0dHR2jEyMtLsklQsFmvB45xKpaIwDJt+LQAAMLWm72p75MgRnTp1SjfccEPtXKVS0cGDB/WZz3xGY2Njchyn9ruenh719PQ0u4w6mUxG6XS6LoA4jiPf91t6XQAAcKGm3/m4/fbbdezYMQ0PD9eOG2+8URs2bNDw8HBd8LDFdV3l8/natR3HUS6Xk+u61msBAGC2a/qdj7lz5+raa6+tO3fFFVfoqquuuuC8TZs3b9bq1asVhqF83yd4AAAQk6aHj3bmui6hAwCAmFkJHwcOHLBxGQAAkADs7QIAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACr2m5vlyiKJEnlcjnmSgAAwHSdm7fPzeNTabvwcfr0aUnSwMBAzJUAAIBGnT59Wn19fVO2SUXTiSgWVatVPf/885o7d65SqVRT/3a5XNbAwIBGRkbU29vb1L/dDjq9f1Ln95H+JV+n95H+JV+r+hhFkU6fPq3+/n6l01M/1dF2dz7S6XTLt73v7e3t2P+ppM7vn9T5faR/ydfpfaR/ydeKPl7sjsc5PHAKAACsInwAAACrZlX46Onp0b333quenp64S2mJTu+f1Pl9pH/J1+l9pH/J1w59bLsHTgEAQGebVXc+AABA/AgfAADAKsIHAACwivABAACs6rjwsXv3bl1zzTWaM2eObr75Zv3oRz+asv2//Mu/6DWveY3mzJmj6667Tt/+9rctVTozjfRvz549SqVSdcecOXMsVtuYgwcPat26derv71cqldK+ffsu+poDBw7ohhtuUE9Pj3zf1549e1pe56VotI8HDhy4YAxTqZROnjxpp+AGDA4O6k1vepPmzp2rq6++WuvXr9dzzz130dcl6T04kz4m6X348MMP6/rrr699+dTy5cv1ne98Z8rXJGn8pMb7mKTxm8gDDzygVCqlHTt2TNnO9jh2VPj453/+Z91999269957NTQ0pGXLlmn16tU6derUhO2ffvpp3Xnnndq8ebOOHj2q9evXa/369frxj39sufLpabR/0vg32P3mN7+pHSdOnLBYcWPOnj2rZcuWaffu3dNqf/z4ca1du1YrV67U8PCwduzYob/6q7/SE0880eJKZ67RPp7z3HPP1Y3j1Vdf3aIKZ+6pp57Stm3b9Oyzz2r//v166aWX9Gd/9mc6e/bspK9J2ntwJn2UkvM+dF1XDzzwgI4cOaLDhw8rm83qjjvu0E9+8pMJ2ydt/KTG+yglZ/zOd+jQIeVyOV1//fVTtotlHKMOctNNN0Xbtm2r/VypVKL+/v5ocHBwwvZvfetbo7Vr19adu/nmm6OtW7e2tM6ZarR/jz76aNTX12epuuaSFO3du3fKNh/4wAeipUuX1p1729veFq1evbqFlTXPdPpYKBQiSdFvf/tbKzU106lTpyJJ0VNPPTVpm6S9B883nT4m+X0YRVF05ZVXRv/wD/8w4e+SPn7nTNXHpI7f6dOno0wmE+3fvz+67bbbou3bt0/aNo5x7Jg7H3/4wx905MgRrVq1qnYunU5r1apVeuaZZyZ8zTPPPFPXXpJWr149afs4zaR/knTmzBl5nqeBgYGLpvukSdL4XarXv/71Wrhwof70T/9UP/jBD+IuZ1pGR0clSfPmzZu0TdLHcDp9lJL5PqxUKnr88cd19uxZLV++fMI2SR+/6fRRSub4bdu2TWvXrr1gfCYSxzh2TPj47//+b1UqFc2fP7/u/Pz58yddHz958mRD7eM0k/4tWbJEn/vc5/T1r39dn//851WtVnXLLbfIGGOj5JabbPzK5bJ+97vfxVRVcy1cuFCPPPKIvvrVr+qrX/2qBgYGtGLFCg0NDcVd2pSq1ap27NihW2+9Vddee+2k7ZL0HjzfdPuYtPfhsWPH9LKXvUw9PT1697vfrb179+p1r3vdhG2TOn6N9DFp4ydJjz/+uIaGhjQ4ODit9nGMY9vtaovmWb58eV2av+WWW/Ta175WuVxO999/f4yVYbqWLFmiJUuW1H6+5ZZb9Itf/EKf/OQn9U//9E8xVja1bdu26cc//rG+//3vx11Ky0y3j0l7Hy5ZskTDw8MaHR3VV77yFW3atElPPfXUpJNzEjXSx6SN38jIiLZv3679+/e39YOxHRM+/viP/1iO4+iFF16oO//CCy9owYIFE75mwYIFDbWP00z6d77LLrtMb3jDGxSGYStKtG6y8evt7dXll18eU1Wtd9NNN7X1pP6e97xH3/zmN3Xw4EG5rjtl2yS9B/+/Rvp4vnZ/H3Z3d8v3fUnSG9/4Rh06dEif+tSnlMvlLmib1PFrpI/na/fxO3LkiE6dOqUbbrihdq5SqejgwYP6zGc+o7GxMTmOU/eaOMaxY5Zduru79cY3vlFPPvlk7Vy1WtWTTz456Vre8uXL69pL0v79+6dc+4vLTPp3vkqlomPHjmnhwoWtKtOqJI1fMw0PD7flGEZRpPe85z3au3evvve972nx4sUXfU3SxnAmfTxf0t6H1WpVY2NjE/4uaeM3man6eL52H7/bb79dx44d0/DwcO248cYbtWHDBg0PD18QPKSYxrFlj7LG4PHHH496enqiPXv2RP/1X/8VbdmyJXr5y18enTx5MoqiKNq4cWO0c+fOWvsf/OAHUVdXV/Txj388+ulPfxrde++90WWXXRYdO3Ysri5MqdH+3XfffdETTzwR/eIXv4iOHDkSvf3tb4/mzJkT/eQnP4mrC1M6ffp0dPTo0ejo0aORpOgTn/hEdPTo0ejEiRNRFEXRzp07o40bN9ba//KXv4z+6I/+KHr/+98f/fSnP412794dOY4Tffe7342rCxfVaB8/+clPRvv27YuKxWJ07NixaPv27VE6nY7+7d/+La4uTOqv//qvo76+vujAgQPRb37zm9rx4osv1tok/T04kz4m6X24c+fO6KmnnoqOHz8e/ed//me0c+fOKJVKRf/6r/8aRVHyxy+KGu9jksZvMud/2qUdxrGjwkcURdHf//3fR4sWLYq6u7ujm266KXr22Wdrv7vtttuiTZs21bX/8pe/HL361a+Ouru7o6VLl0bf+ta3LFfcmEb6t2PHjlrb+fPnR3/+538eDQ0NxVD19Jz7WOn5x7k+bdq0KbrtttsueM3rX//6qLu7O3rlK18ZPfroo9brbkSjfXzwwQejV73qVdGcOXOiefPmRStWrIi+973vxVP8RUzUL0l1Y5L09+BM+pik9+G73vWuyPO8qLu7O3rFK14R3X777bVJOYqSP35R1HgfkzR+kzk/fLTDOKaiKIpad18FAACgXsc88wEAAJKB8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMCq/wWk8NYI2FLn9gAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt \n", + "plt.plot(x,y, '.k')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAq20lEQVR4nO3df3DT933H8dfXKBZZF7shjjA+C0PbjDg/SrAMOcjukhBWjnEs7LZ16anUdN6RdnSEpQslvSW5HutouF7TrONCttI4WKR2txbWS0YzljiwNomKpLCRppeGBqtOKDHeMZuQVQT82R9UqmVkY0lffb/6Ss/HnU+29LX0+fJF6MXnx/tjGWOMAAAAHFLjdgMAAEB1IXwAAABHET4AAICjCB8AAMBRhA8AAOAowgcAAHAU4QMAADiK8AEAABzlc7sB442Ojur48eO64oorZFmW280BAABTYIzR6dOn1dTUpJqayfs2yi58HD9+XMFg0O1mAACAAgwMDKi5uXnSY8oufFxxxRWSLjS+rq7O5dYAAICpGBkZUTAYzHyOT6bswkd6qKWuro7wAQCAx0xlygQTTgEAgKMIHwAAwFGEDwAA4CjCBwAAcBThAwAAOIrwAQAAHEX4AAAAjiJ8AAAARxE+AACAowgfAAB4QCwW09KlSxWLxdxuStEIHwAAeMCuXbvU19en7u5ut5tStLLb2wUAAFyQTCY1NDQky7LU29srSerp6VFHR4eMMWpoaFBLS4vLrcwf4QMAgDI1Z86czPfpDdtOnjypUCiUud8Y43SzisawCwAAZSoSicjnu9BPkA4Z6Vufz6dIJOJa24pBzwcAAGUqHA6rtbU1q6cjLRqNqq2tzYVWFY+eDwAAPKCmpibr1su8fwYAAFSwQCCgxsZGhUIh7dixQ6FQSI2NjQoEAm43rWCWKbOZKiMjI6qvr9fw8LDq6urcbg4AAK5LpVKqra2VZVkyxujs2bPy+/1uNytLPp/fzPkAAKDMjQ0almWVXfDIF8MuAAB4kJcrnhI+AADwIC9XPGXYBQAAj6iUiqeEDwAAPKJSKp4y7AIAgEdUSsVTej4AAPCISql4Ss8HAAAe5OWKp95rMQAAVawSKp5S4RQAAI8px4qnVDgFAKCCeb3iKcMuAADAUYQPAADgKMIHAABwFOEDAAA4Ku/wcfDgQa1atUpNTU2yLEt79+6d8NjPfOYzsixLX//614toIgAAqCR5h48zZ85o/vz52r59+6TH7dmzRy+//LKampoKbhwAAKg8eS+1XbFihVasWDHpMW+//bb+8i//Us8++6xWrlxZcOMAAEDlsb3Ox+joqNasWaP77rtP119//SWPT6VSSqVSmZ9HRkbsbhIAACgjtk84ffjhh+Xz+bRhw4YpHb9161bV19dnvoLBoN1NAgAAZcTW8BGPx/Xoo4+qq6tLlmVN6Xfuv/9+DQ8PZ74GBgbsbBIAACgztoaP//zP/9Tg4KBmz54tn88nn8+nZDKpz3/+85ozZ07O3/H7/aqrq8v6AgDAbrFYTEuXLlUsFnO7KVXP1jkfa9as0bJly7LuW758udasWaNPf/rTdr4UAAB52bVrl/r6+tTd3a329na3m1PV8g4f7777ro4ePZr5+dixYzp8+LBmzJih2bNn66qrrso6/rLLLlNjY6PmzZtXfGsBAMhDMpnU0NCQLMtSb2+vJKmnp0cdHR0yxqihoUEtLS0ut7L65B0+YrGYbr/99szP9957rySpo6NDXV1dtjUMAIBijR3yT89FPHnypEKhUOZ+Y4zTzap6eYeP2267La8L1d/fn+9LAABgi0gkorVr1+rcuXOZz670rc/n4z/NLrG9zgcAAOUiHA6rtbU1q6cjLRqNqq2tzYVWgY3lAABVoaamJusW7uEKAAAqWiAQUGNjo0KhkHbs2KFQKKTGxkYFAgG3m1a1LFNmM21GRkZUX1+v4eFhan4AAGyRSqVUW1sry7JkjNHZs2fl9/vdblZFyefzmzkfAICKNzZoWJZF8HAZwy4AAMBRhA8AAOAowgcAAHAU4QMAgCpQThvrET4AAKgCYzfWcxvhAwCAMmRHT0UymVQ8HlcikchsrNfd3a2FCxequ7tbyWTSrubmhaW2AACUobE9Fe3t7QU9R66N9U6dOqVYLKZPfepTktzZWI+eDwAAykSunoqenh4lEgnF4/G8eyoikYh8vgv9DLlCxpYtWwp63mJR4RQAgDKR7p1If2+Mydym5fuxnUgkcm6sN16xcSCfz296PgAAKJDdK0hy9VSkb30+nyKRiC2vM1apnncyhA8AAApk9wqScDisaDSa87FoNKpwOJz3c6Y31rvuuutsfd5iMOEUAIA8JJNJDQ0NybKsrHkZHR0dMsaooaFBLS0tRb9OTU2NRkdHM7eFam5uVn9/v1599VW1t7dnhnGKfd5iED4AAMhDrhUkJ0+ezJpXUcz8iXRPRTAYVGdnp3bu3KmBgQEFAoGCn9Pv92vmzJm2P2+hmHAKAEAedu/erbVr1+rcuXMXPebz+dTV1VX0MEYqlVJtbW2ml+Ls2bO27MRbqueV8vv8pucDAIA8hMNhtba25lxBEo1G1dbWVvRrjA0ElmXZFhBK9bz5YsIpAAAFqqmpybrF1PCnBQBAntLzMkKhkHbs2KFQKKTGxkZX5k94EXM+AAAoQCnnT3gRcz4AACixcpk/4UUMuwAAAEcRPgAAgKMIHwAAwFGEDwAA4CjCBwAAcBThAwCAMWKxmJYuXapYLOZ2UyoW4QMAUNacDgO7du1SX1+furu7HXm9akSdDwBAWRsbBtrb20vyGslkUkNDQ7IsS729vZKknp4edXR0yBijhoYGtbS0lOS1qxEVTgEAZWdsGFixYoUGBwcVCAS0b9++koQBy7KyvjfGZG7TyuzjsuxQ4RQA4Glz5szJfJ8OBidPnszaSdbOMBCJRLR27VqdO3cu87zpW5/Pp66uLtteC8z5AACUoUgkIp/vwv+Pc4WBSCRi6+uFw2FFo9Gcj0WjUYXDYVtfr9oRPgAAZcfNMFBTU5N1C/vxJwsAKGvpEJAefnnttddK8jqBQECNjY0KhULasWOHQqGQGhsbFQgESvJ61YzwAQAoS+PDwNVXXy1J6uvrK8nrNTc3q7+/X9FoVHfffbei0aj6+/vV3NxckterZkw4BQCUTCwW06ZNm7Rt27a8l8k2NzfrwIEDGhkZyRoCefrpp5VIJEqy6sXv92e+tywr62fYh/ABACiZYmt0zJs3L/O9E6te4AyGXQAAtkomk4rH40okElkFuxKJhOLxuJLJ5JSfy+lVL3AGRcYAALayu2BXIpHI6ulIi8fjamtry7t9xQwFYWL5fH7T8wEAsFWpeivsWgLL3i3uY84HAMBW4XBYra2tOXsrotFo3r0V6VUvwWBQnZ2d2rlzpwYGBvJaAsveLeWF8AEAKJmamhqNjo5mbguRXgJbW1sry7K0bt06nT17Nq+VKE6Xa8fkGHYBANjO7oJdfr8/ExoKWQLLxNXywoRTAEBJpFKpTG+FMSbv3gq72T1xFdmYcAoAcF2xvRWlwt4t7uNPHgBQFdi7pXww7AIAqBrlNhRUSfL5/Ga1CwCgarB3S3lg2AUAADiK8AEAABxF+AAAAI4ifAAAAEflHT4OHjyoVatWqampSZZlae/evZnH3n//fX3hC1/QjTfeqA984ANqamrSpz71KR0/ftzONgMAAA/LO3ycOXNG8+fP1/bt2y967L333lMikdADDzygRCKh733ve3r99df1B3/wB7Y0FgAAeF9RdT4sy9KePXu0evXqCY85dOiQFi1apGQyqdmzZ1/yOanzAQCA95RVefXh4WFZlqUPfvCDpX4pAADgASUtMvarX/1KX/jCF/SJT3xiwhSUSqWUSqUyP4+MjJSySQAAwGUl6/l4//339fGPf1zGGD322GMTHrd161bV19dnvoLBYKmaBAAAykBJwkc6eCSTSe3fv3/SsZ/7779fw8PDma+BgYFSNAkAAJQJ24dd0sHjjTfeUF9fn6666qpJj/f7/dTWBwCgiuQdPt59910dPXo08/OxY8d0+PBhzZgxQ7NmzdIf//EfK5FI6Omnn9b58+d14sQJSdKMGTNUW1trX8sBAIAn5T3sEovFtGDBAi1YsECSdO+992rBggV68MEH9fbbb+v73/++3nrrLd10002aNWtW5uvFF1+0vfEAAO+JxWJaunSpYrGY202BS/Lu+bjttts0WWmQIsqGAACqwK5du9TX16fu7m61t7e73Ry4oKRLbQEAkKRkMqmhoSFZlqXe3l5JUk9Pjzo6OmSMUUNDg1paWlxuJZxC+AAAlNycOXMy31uWJUk6efKkQqFQ5n56zqsHu9oCAEouEonI57vw/910yEjf+nw+RSIR19oG5xE+AAB5y3fSaDgcVjQazflYNBpVOBy2s3koc4QPAEDexk4azVdNTU3WLaoPcz4AAFNS7KTRQCCgxsZGBYNBdXZ2aufOnRoYGFAgEHDqFFAmLFNmM3zy2ZIXAOCc9ETR9PfGmMxt2qU+UlKplGprazO/d/bsWapcV4h8Pr/p8wIATMlkk0Yl6WMf+9gln8Pv92dCjGVZBI8qRfgAAEzJZJNGJenw4cNKJBKKx+NKJpMOtgxew5wPAIAtqNuBqaLnAwCqSLH7qqQnjX7oQx+6aLUKdTswVYQPAKgixSyRlaTm5mb19/fr6NGjOnToUM5jqNuBSyF8AECFSyaTisfjSiQSWUtkC52fMXbSqHRx3Y7XXnuNXWsxKeZ8AECFK9W+KhPV7Xj++efZtRaTInwAQIWLRCJau3atzp07l3Nfla6uroKeNz0EU1tbq1/84hcKhUI6d+6c7rzzTknsWouJUWQMAKpAIpHI6ulIi8fjamtrK/r57ShABm+jyBgAIKdS7avCrrXIB+EDAKpAen5GKBTSjh07FAqF1NjYaNu+Kuxai3ww5wMAqsDY+RmWZWndunUl21elpqZGo6OjmVtgPHo+AKBKlHpflVL3rqByMOEUAGAbdq2tXvl8fjPsAgCwzdigwa61mAjDLgAAwFGEDwDwsGI3igPcQPgAAA8rdqO48QgzcAJzPgDAY5LJpIaGhmRZVtZGcXaUMh8bZtiXBaVC+AAAj7F7o7iJwszChQv16KOP6sEHH9SqVavsaTwgltoCgOfs3r07s1HceOmN4vKpKMq+LLADe7sAQAWzu5T5ZPuySFJdXZ0SiYTi8biSyWSBrQZ+g2EXAPAwO0qZh8Nhtba25tz1VpJOnz5d8JAOkAs9HwDgQaUqZT52CCaN3WlhN8IHAORQ7ktO0xvFRaNR3X333YpGo+rv71dzc3NBz5cOM+3t7friF7+Y8xh2p4VdCB8AkIPd9TPyMdXgY+dGcWPDzB/90R9JujCkM/YWsAt/owDg15LJpOLxuBKJRNaSU6cnW7oVfNJhht1pUWostQWAX3NzyenYWhsrVqzQ4OCgAoGA9u3bV3ThsEKwOy3yxa62AFCASCSSqZ8xfslpun5GqdhdOKxY7E6LUmLYBQB+ze76GfmYrNYGq0xQaQgfAJCD05Mt3Qw+gNMIHwAwRjlMtmSVCSodf7MBYAy762fkw63gU+41TVB5CB8AMM6RI0d0xx13KBaLOTrZ0q3g42ZNE1QnVrsAwDhjP4zb29sdfW2nVpmMXdo7tqZJR0eHK0t7UV2o8wEAKr86G6XmZk0TVKZ8Pr8ZdgEAXaiz0d7erlAopJMnT0r6TZ2N9vb2rDocaV6eK8HSXriJ8AEAKuzD2MtzJVjaCzcx5wMAdOHDuLW1NauiaFo0GlVbW5ukypwrUVNTo9HR0cwtUGqEDwAYZ7IP43Irg16M9NLeYDCozs5O7dy5UwMDA2wgh5IjfADAr03lw9jN/V/sll7am95Abt26dWwgB0ew2gUAxpjKbq6JRCLn8Ew8Hs8MzwDVhtUuAFAgv9+fGU65VJ0NyqADheEdAwB5Kof9XwAvY9gFAAowleEZoJrk8/nNhFMAKIBTZdCBSsSwC4CK4OVqo0C1IXwAqAherjYKVBuGXQB4ViVWGwWqQd49HwcPHtSqVavU1NQky7K0d+/erMeNMXrwwQc1a9YsXX755Vq2bJneeOMNu9oLABmFbAYHwH15h48zZ85o/vz52r59e87Ht23bpr//+7/Xjh07FI1G9YEPfEDLly/Xr371q6IbCwBjsTMr4E1FLbW1LEt79uzR6tWrJV140zc1Nenzn/+8/vqv/1qSNDw8rJkzZ6qrq0t33XXXJZ+TpbYAxovFYtq0aZO2bdum9vb2rMeoNgqUB9cqnB47dkwnTpzQsmXLMvfV19fr5ptv1ksvvZTzd1KplEZGRrK+AGCsqUwmpdoo4B22vktPnDghSZo5c2bW/TNnzsw8Nt7WrVtVX1+f+QoGg3Y2CYBHJZNJxeNxJRKJrMmkiURC8XhcyWRSEtVGAS9yfbXL/fffr3vvvTfz88jICAEEwJS3rmdnVsB7bO35aGxslCS98847Wfe/8847mcfG8/v9qqury/oCgHwmk+azGVwxKGQG2MPW8DF37lw1Njbqueeey9w3MjKiaDSqxYsX2/lSACpcOBxWNBrN+Vg0GlU4HHa4RRQyA+yS97DLu+++q6NHj2Z+PnbsmA4fPqwZM2Zo9uzZ2rhxo/72b/9W11xzjebOnasHHnhATU1NmRUxAJCvmpoajY6OZm6dRCEzwH55h49YLKbbb78983N6vkZHR4e6urq0adMmnTlzRuvWrdP//u//6nd/93f1gx/8QNOnT7ev1QCqQnoyaTAYVGdnp3bu3KmBgQFHJ5NOde4JgKkrqs5HKVDnA8BYbm9dv3v3bq1du1bnzp276DGfz6euri5XhoCAcpPP57frq10AYDJub10fDofV2tqas5BZNBqlkBlQAKrxAMAUUcgMsAfvIAC4BAqZAfZizgcATIHbc0+AcsecDwCwmdtzT4BKwrALgItUaiXPSj0vwGsIHwAuUqmVPCv1vACvYc4HAEnZlTxXrFihwcFBBQIB7du3z9OVPCv1vIByk8/nN+EDgKTfVO9Mf2+Mydymldk/F1NSqecFlJt8Pr8ZdgEgKb9dZL2kUs8L8DJ6PgBkJBKJnJU84/G4pyt5TnRe3d3d+uQnP+lCi4DKQ88HgKJUaiXP8efzzDPPuNQSoLpV1r8sAIpSqZU8A4GArr76al177bX64he/mBmG2b9/vxKJhOLxuJLJpMutBKoHwy4AslRqJU8mngKlxbALgIL5/f7MB3UlVfJk4ilQPggfAGxTzhVEw+GwotFozsei0ajC4bDDLQKqF+EDgG28UkG0UifUAl7BxnIAijK2gmhvb68kqaenRx0dHWVXQTQ9oTYYDKqzs1M7d+7UwMCA5yfUAl7DhFMARfHaRM5KnVALuI0JpwAc47WJnJU6oRbwEoZdABQlHA6rtbU1ZwXRaDTq6cqoAEqDng8AtmEiJ4Cp4F8IAEWr1MqoAEqDCacAbMFETqC65fP5zZwPALYYGzSYyAlgMgy7AAAARxE+AACAowgfAADAUYQPAADgKMIHAABwFOEDAAA4ivABAAAcRfgAAACOInwAAABHET4AAICjCB8AAMBRhA8AAOAowgcAAHAU4QMAADiK8AEAABxF+AAAAI4ifAAAAEcRPgAAgKMIHwAAwFGEDwAA4CjCBwAAcBThAwAAOIrwAQAAHEX4AAAAjiJ8AAAARxE+AACAowgfAADAUYQPAADgKMIHAABwFOEDyEMsFtPSpUsVi8XcbgoAeJbt4eP8+fN64IEHNHfuXF1++eX68Ic/rC1btsgYY/dLAY7btWuX+vr61N3d7XZTAMCzfHY/4cMPP6zHHntMTz75pK6//nrFYjF9+tOfVn19vTZs2GD3ywEll0wmNTQ0JMuy1NvbK0nq6elRR0eHjDFqaGhQS0uLy60EAO+wPXy8+OKLuvPOO7Vy5UpJ0pw5c/Ttb39bP/7xj+1+KcARc+bMyXxvWZYk6eTJkwqFQpn76dkDgKmzfdhlyZIleu655/Szn/1MkvRf//Vf+uEPf6gVK1bkPD6VSmlkZCTrCygnkUhEPt+FnJ4OGelbn8+nSCTiWtsAwIts7/nYvHmzRkZGdO2112ratGk6f/68vvzlLyscDuc8fuvWrfrSl75kdzMA24TDYbW2tmb1dKRFo1G1tbW50CoA8C7bez6+853vaPfu3XrqqaeUSCT05JNP6qtf/aqefPLJnMfff//9Gh4eznwNDAzY3STANjU1NVm3AID82d7zcd9992nz5s266667JEk33nijksmktm7dqo6OjouO9/v98vv9djcDsFUgEFBjY6OCwaA6Ozu1c+dODQwMKBAIuN00APAc28PHe++9d9H/CqdNm6bR0VG7XwpwTHNzs/r7+1VbWyvLsrRu3TqdPXuW4AwABbC973jVqlX68pe/rGeeeUb9/f3as2ePvva1r+kP//AP7X4pwFF+vz+z2sWyrJzBgyJkAHBptvd8fOMb39ADDzygv/iLv9Dg4KCampp0991368EHH7T7pYCyM7YIWXt7u9vNAYCyZJkyK1AwMjKi+vp6DQ8Pq66uzu3mAJc0tgjZihUrNDg4qEAgoH379lGEDEDVyOfz2/aeD8DLYrGYNm3apG3btk2554IiZACQH9YLAmMUsncLRcgAID8Mu6Dq2TFskkgkchYhi8fjFCEDUBUYdgHyYOewSU1NjUZHRzO3AICLMewCz7JrWasdwybpImShUEg7duxQKBRSY2MjRcgAIAd6PuBZdi1rtWPvFoqQAcDUET7gKWPnZ/T29kqSenp61NHRYcuy1mKGTcYGjYmKkAEACB/wmFIta2XvFgBwDuEDnhKJRLR27VqdO3cu5/yMrq6ugp6XYRMAcA7hA55ix/yMiTBsAgDOYLULPCu9e/L4XZQnw8ZvAOA+wgc8p5hlrYVUMB2PAAMAxWHYBZ6T7/wMu1fIsHMtABSH8uqoeOlVMenvjTGZ27RLvQ3YuRYAJkd5dWAMO1bIsHMtANiHOR+oeOFwWNFoNOdj0WhU4XD4ks/BzrUAYB/CB6pKIStkJHsCDADgAsIHqoKdG78VGmAAABcw5wNVwY4KppRgBwB7sNoFyEMqlcoEGGMMJdgB4NdY7QKUCCXYAaB4DFoDAABHET4AAICjCB8AAMBRhA9UPDaCA4DyQvhAxbNjJ1sAgH0IH6gI43s3ksmk4vG4EolE1k62iURC8XhcyWTSzeYCQFVjqS0qwvht7tkIDgDKFz0f8KzJeje2bNnCRnAAUKbo+YBnTaV3I5doNKq2trZSNg0AMAl6PlA2cq1KmWylyqW2ud+yZYskNoIDgHLDv8YoG7lWpUy2UuVS29yvXbvWtp1sAQD2YdgFrkomkxoaGpJlWZl5G7t379aiRYskSU899ZSkC3M5Ojo6ZIzR8ePH9cgjj2jbtm1ZvRqjo6OZW8menWwBAPZjV1u4Kj1XI/39RH8dcz22YcMG3XfffVq4cOFF29wfOnRIzc3NJW07AOA38vn8ZtgFrso1byOX9GPTpk3L/KXu6enR4OCgvve97+lv/uZv1Nvbq+3bt6u/v5/gAQBljPCBkphqSfPJ5m3kcv78eZ0+fVrSb1a2LFmyRHfeeaf6+voUiUQYVgGAMkf4QEkUUtI816qU9Pdjh2fGr2wZ3xtCFVMAKG/M+YBtxk4eXbFihQYHBxUIBLRv3z4ZY9TQ0KCWlpaLfu+tt966aN5Gf3+/pAu1PNL3vfnmm/qf//mfnK+dnhMyfm5Imf31BoCKlc/nN+EDtsk1eXR8GDh06JDa29sv+t1UKpVZlWKM0dmzZyUp675oNKrFixdnVrRMNkHV5/Opq6tL4XDY5rMEAOTChFO4YrKiX+lgMtEwjN/vzxxjWZb8fv9F9zU3N2fV7Whvb9dVV12V8/mi0SjBAwDKFD0fsFUikchZ3vzKK6/UqVOnpjQMM5nxPSTje0PSt/F4nBLqAOCgfD6/KTKGkhhb7EuSTp06Jan4nWXHrmQZ2xsyvs4HVUwBoHwRPmCrQCCQCQM33HCDnnjiiazHx+690tXVVfTrUcUUALyHYRfYbuzQSDwezznBlGERAKgsTDiFq8ZPFJXYWRYA8Bt8EqCk0sMw7CwLAEhj2AUll6uGB3MyAKCyMOyCspKrhsdEpronDADAuwgfKCuF7AkDAPAWltrCdWP3hOnt7ZV0YYO4jo6OgouRAQDKF+EDrpszZ07m+/TwTLHFyAAA5YthF7husj1hfD6fIpGIa20DANiPng+4LhwOq7W1NeeeMNFolGJkAFBh6PlAWaEYGQBUPv6FR1mgGBkAVA+KjKFsUIwMALzL9SJjb7/9tj75yU/qqquu0uWXX64bb7yRolG4pHyKkQEAvMv2CaenTp3SLbfcottvv1379u3T1VdfrTfeeENXXnml3S8FAAA8yPbw8fDDDysYDOqJJ57I3Dd37ly7XwYAAHiU7cMu3//+99Xe3q4/+ZM/USAQ0IIFC/RP//RPEx6fSqU0MjKS9QUAACqX7eHjzTff1GOPPaZrrrlGzz77rD772c9qw4YNevLJJ3Mev3XrVtXX12e+gsGg3U0CAABlxPbVLrW1tWpvb9eLL76YuW/Dhg06dOiQXnrppYuOT6VSSqVSmZ9HRkYUDAZZ7QIAgIe4utpl1qxZuu6667Lua21t1S9+8Yucx/v9ftXV1WV9AQCAymV7+Ljlllv0+uuvZ933s5/9jF1JAQCApBKEj7/6q7/Syy+/rL/7u7/T0aNH9dRTT+kf//EftX79ertfCgAAeJDt4WPhwoXas2ePvv3tb+uGG27Qli1b9PWvf13hcNjulwIAAB5EeXUAAFA018urAwAATITwAQAAHEX4AAAAjiJ8AAAARxE+AACAowgfAADAUYQPAADgKMIHAABwFOEDAAA4ivABAAAcRfgAAACOInwAAABHET4AAICjCB8AAMBRhA8AAOAowgcAAHAU4QMAADiK8AEAABxF+AAAAI4ifAAAAEcRPgAAgKMIHxUgFotp6dKlisVibjcFAIBLInxUgF27dqmvr0/d3d1uNwUAgEvyud0AFCaZTGpoaEiWZam3t1eS1NPTo46ODhlj1NDQoJaWFpdbCQDAxQgfHjVnzpzM95ZlSZJOnjypUCiUud8Y43SzAAC4JIZdPCoSicjnu5Ad0yEjfevz+RSJRFxrGwAAk6Hnw6PC4bBaW1uzejrSotGo2traXGgVAACXRs9HBaipqcm6BQCgnPFp5WGBQECNjY0KhULasWOHQqGQGhsbFQgE3G4aAAATskyZzUocGRlRfX29hoeHVVdX53Zzyl4qlVJtba0sy5IxRmfPnpXf73e7WQCAKpPP5zdzPjxubNCwLIvgAQAoewy7AAAARxE+AACAowgfAADAUVUVPtiADQAA91VV+LBrAzZCDAAAhav41S6l2IBtbIhpb28vRbMBAKhYFV/nI73pWvp7Y0zmNm0qfwRjQ8yKFSs0ODioQCCgffv2sYssAKDqUedjjEgkorVr1+rcuXM5N2Dr6uqa0vOwiywAAPao+Dkf4XBY0Wg052PRaFThcHhKz8MusgAA2KPiw8dYxWzAZleIAQCg2lVF+LB7AzZ2kQUAoHAVP+dDkpqbm9Xf35/ZgG3dunUFbcCWDjHBYFCdnZ3auXOnBgYG2EUWAIA8VPxqF7uxiywAABdjtUsJsYssAADFYdICAABwFOEDAAA4ivABAAAcRfgAAACOInwAAABHET4AAICjCB8AAMBRhA8AAOAowgcAAHAU4QMAADiq7Mqrp7eaGRkZcbklAABgqtKf21PZMq7swsfp06clScFg0OWWAACAfJ0+fVr19fWTHlN2u9qOjo7q+PHjuuKKK2RZ1oTHjYyMKBgMamBgoCx3v7UL51lZOM/KwnlWFs6zOMYYnT59Wk1NTaqpmXxWR9n1fNTU1Ki5uXnKx9fV1VX0X5I0zrOycJ6VhfOsLJxn4S7V45HGhFMAAOAowgcAAHCUZ8OH3+/XQw89JL/f73ZTSorzrCycZ2XhPCsL5+mcsptwCgAAKptnez4AAIA3ET4AAICjCB8AAMBRhA8AAOCosg4f27dv15w5czR9+nTdfPPN+vGPfzzp8f/8z/+sa6+9VtOnT9eNN96of/u3f3OopcXJ5zy7urpkWVbW1/Tp0x1sbWEOHjyoVatWqampSZZlae/evZf8nRdeeEFtbW3y+/36yEc+oq6urpK3s1j5nucLL7xw0fW0LEsnTpxwpsEF2Lp1qxYuXKgrrrhCgUBAq1ev1uuvv37J3/Pa+7OQ8/Ti+/Oxxx7TRz/60UzBqcWLF2vfvn2T/o7XrqWU/3l68Vrm8pWvfEWWZWnjxo2THuf0NS3b8NHb26t7771XDz30kBKJhObPn6/ly5drcHAw5/EvvviiPvGJT6izs1OvvPKKVq9erdWrV+vVV191uOX5yfc8pQtV6X75y19mvpLJpIMtLsyZM2c0f/58bd++fUrHHzt2TCtXrtTtt9+uw4cPa+PGjfrzP/9zPfvssyVuaXHyPc+0119/PeuaBgKBErWweAcOHND69ev18ssva//+/Xr//ff1sY99TGfOnJnwd7z4/izkPCXvvT+bm5v1la98RfF4XLFYTEuXLtWdd96pn/zkJzmP9+K1lPI/T8l713K8Q4cO6fHHH9dHP/rRSY9z5ZqaMrVo0SKzfv36zM/nz583TU1NZuvWrTmP//jHP25WrlyZdd/NN99s7r777pK2s1j5nucTTzxh6uvrHWpdaUgye/bsmfSYTZs2meuvvz7rvj/90z81y5cvL2HL7DWV8+zr6zOSzKlTpxxpUykMDg4aSebAgQMTHuPV9+dYUznPSnh/GmPMlVdeab75zW/mfKwSrmXaZOfp9Wt5+vRpc80115j9+/ebW2+91dxzzz0THuvGNS3Lno+zZ88qHo9r2bJlmftqamq0bNkyvfTSSzl/56WXXso6XpKWL18+4fHloJDzlKR3331XLS0tCgaDl0zuXuXF61mMm266SbNmzdLv/d7v6Uc/+pHbzcnL8PCwJGnGjBkTHlMJ13Mq5yl5+/15/vx59fT06MyZM1q8eHHOYyrhWk7lPCVvX8v169dr5cqVF12rXNy4pmUZPoaGhnT+/HnNnDkz6/6ZM2dOOBZ+4sSJvI4vB4Wc57x58/Stb31L//qv/6pIJKLR0VEtWbJEb731lhNNdsxE13NkZET/93//51Kr7Ddr1izt2LFD3/3ud/Xd735XwWBQt912mxKJhNtNm5LR0VFt3LhRt9xyi2644YYJj/Pi+3OsqZ6nV9+fR44c0W//9m/L7/frM5/5jPbs2aPrrrsu57Fevpb5nKdXr6Uk9fT0KJFIaOvWrVM63o1rWna72mJyixcvzkrqS5YsUWtrqx5//HFt2bLFxZahEPPmzdO8efMyPy9ZskQ///nP9cgjj6i7u9vFlk3N+vXr9eqrr+qHP/yh200pqamep1ffn/PmzdPhw4c1PDysf/mXf1FHR4cOHDgw4QezV+Vznl69lgMDA7rnnnu0f//+sp4gW5bho6GhQdOmTdM777yTdf8777yjxsbGnL/T2NiY1/HloJDzHO+yyy7TggULdPTo0VI00TUTXc+6ujpdfvnlLrXKGYsWLfLEh/nnPvc5Pf300zp48KCam5snPdaL78+0fM5zPK+8P2tra/WRj3xEkhQKhXTo0CE9+uijevzxxy861svXMp/zHM8r1zIej2twcFBtbW2Z+86fP6+DBw/qH/7hH5RKpTRt2rSs33HjmpblsEttba1CoZCee+65zH2jo6N67rnnJhyfW7x4cdbxkrR///5Jx/PcVsh5jnf+/HkdOXJEs2bNKlUzXeHF62mXw4cPl/X1NMboc5/7nPbs2aPnn39ec+fOveTvePF6FnKe43n1/Tk6OqpUKpXzMS9ey4lMdp7jeeVa3nHHHTpy5IgOHz6c+Wpvb1c4HNbhw4cvCh6SS9e0ZFNZi9TT02P8fr/p6uoyr732mlm3bp354Ac/aE6cOGGMMWbNmjVm8+bNmeN/9KMfGZ/PZ7761a+an/70p+ahhx4yl112mTly5IhbpzAl+Z7nl770JfPss8+an//85yYej5u77rrLTJ8+3fzkJz9x6xSm5PTp0+aVV14xr7zyipFkvva1r5lXXnnFJJNJY4wxmzdvNmvWrMkc/+abb5rf+q3fMvfdd5/56U9/arZv326mTZtmfvCDH7h1ClOS73k+8sgjZu/eveaNN94wR44cMffcc4+pqakx//Ef/+HWKVzSZz/7WVNfX29eeOEF88tf/jLz9d5772WOqYT3ZyHn6cX35+bNm82BAwfMsWPHzH//93+bzZs3G8uyzL//+78bYyrjWhqT/3l68VpOZPxql3K4pmUbPowx5hvf+IaZPXu2qa2tNYsWLTIvv/xy5rFbb73VdHR0ZB3/ne98x/zO7/yOqa2tNddff7155plnHG5xYfI5z40bN2aOnTlzpvn93/99k0gkXGh1ftJLSsd/pc+to6PD3HrrrRf9zk033WRqa2vNhz70IfPEE0843u585XueDz/8sPnwhz9spk+fbmbMmGFuu+028/zzz7vT+CnKdX6Ssq5PJbw/CzlPL74//+zP/sy0tLSY2tpac/XVV5s77rgj84FsTGVcS2PyP08vXsuJjA8f5XBNLWOMKV2/CgAAQLaynPMBAAAqF+EDAAA4ivABAAAcRfgAAACOInwAAABHET4AAICjCB8AAMBRhA8AAOAowgcAAHAU4QMAADiK8AEAABxF+AAAAI76fwwW71kKZFohAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x = 4*np.random.rand(50, 1)\n", + "y = theta0 + theta1*x+0.5*np.random.randn(50, 1)\n", + "plt.plot(x,y, '*k')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Implementing with `for` \n", + "$$\\theta_1 = \\frac{\\sum_{i=1}^n x^i(y^i-\\bar{y}) }{\\sum_{i=1}^n x^i(x^i-\\bar{x})}$$" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[2.4717291]\n" + ] + } + ], + "source": [ + "# for implementation for computing theta1:\n", + "xAve = x.mean()\n", + "yAve = y.mean()\n", + "num = 0\n", + "den = 0\n", + "for i in range(len(x)):\n", + " num = num + x[i]*(y[i]-yAve)\n", + " den = den + x[i]*(x[i]-xAve)\n", + "theta1Hat = num/den\n", + "print(theta1Hat)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[4.18459936]\n" + ] + } + ], + "source": [ + "# for implementation for theta0:\n", + "# $$ \\theta_0 = \\bar{y} -\\theta_1 \\bar{x} $$\n", + "theta0Hat = yAve - theta1Hat*xAve\n", + "print(theta0Hat)\n", + "#real values are\n", + "#theta0 = 3.9654\n", + "#theta1 = 2.5456" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([2.27654582])" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "total = 0\n", + "for i in range(len(x)):\n", + " total = total + x[i]\n", + "total/len(x)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Implementing OLS by numpy methods" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.4717291029649546\n" + ] + } + ], + "source": [ + "# For theta1:\n", + "# $$\\theta_1 = \\frac{\\sum_{i=1}^n x^i(y^i-\\bar{y}) }{\\sum_{i=1}^n x^i(x^i-\\bar{x})}$$\n", + "num2 = np.sum(x*(y-y.mean()))\n", + "den2 = np.sum(x*(x-x.mean()))\n", + "theta1Hat2 = num2/den2\n", + "print(theta1Hat2)\n", + "\n", + "# Efficacy --> time\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "4.184599360470533" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "theta0Hat2 = yAve-theta1Hat2*xAve\n", + "theta0Hat2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Comparing Model and Data" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCwUlEQVR4nO3df3yT5b3/8Xea0iIIRTyWUlODW6tOZSCggNOBFe1XORydP7Z5EFB61upggG4IuE3kOAXcjhMUpXWRUmBTEWF7bFPmXMGBWmkrjk3U1pFCBoigowW0xeb+/nHRlpa0Tdr87uv5ePQRcuduct2ENh+u6/p8PjbLsiwBAACESUKkBwAAALoXgg8AABBWBB8AACCsCD4AAEBYEXwAAICwIvgAAABhRfABAADCiuADAACEVWKkB9Ca1+vV3r171adPH9lstkgPBwAA+MGyLNXW1io9PV0JCe3PbURd8LF3715lZGREehgAAKAT9uzZI4fD0e45URd89OnTR5IZfN++fSM8GgAA4I+amhplZGQ0fY63J+qCj8allr59+xJ8AAAQY/zZMsGGUwAAEFYEHwAAIKwIPgAAQFhF3Z4Pf1iWpS+//FINDQ2RHkrMsdvtSkxMJI0ZABAxMRd81NfXa9++fTp27FikhxKzevXqpYEDByopKSnSQwEAdEMxFXx4vV7t2rVLdrtd6enpSkpK4n/wAbAsS/X19frkk0+0a9cuZWVldVgIBgCAYIup4KO+vl5er1cZGRnq1atXpIcTk0477TT16NFD1dXVqq+vV8+ePSM9JABANxOT/+3lf+tdw98fACCS+BQCAABhRfABAADCiuAjQj777DMtWLBA+/bti/RQAAAxwuPxqKSkRB6PJ9JD6RKCjwiZOXOm3n77bd19992RHgoAIAa4XC45nU5lZ2fL6XTK5XJFekidRvARAX/4wx9UW1urP/zhD+rXr5/WrFkT6SEBAKKYx+NRXl6evF6vJFN6Ij8/P2ZnQGIq1TZejB8/XuPHj5ckFRUVRXYwAICoV1lZ2RR4NGpoaFBVVZUcDkeERtV5zHwAABDlfBWFtNvtyszMjNCIuobgAwCAKOdwOFRYWCi73S7JBB4FBQUxOeshdffgo6xMys42t2HgcDj01FNPtTj2xhtvqFevXqqurg7LGAAAsSk3N1dut1slJSVyu93Kzc2N9JA6rXvv+SgulkpKpFWrpBEjQv5yI0eO1LZt25ruW5alWbNm6Z577pHT6Qz56wMAYpvD4YjZ2Y6TxX7wYVlSIB1ud++WDh2SbDbpuefMsd/8Rvr2t81znXmmdM45/j1Xr17mefw0atQorVy5sun+qlWrtGfPHs2bN8//8QMAEONiP/g4dkw6/fSuPccnn0hXXBH49x05IvXu7ffpo0aN0ty5c3XkyBHZbDbdf//9+tnPfqbTuzp+AEC34/F4VFlZqaysrJibDeneez7CbPjw4UpISFBFRYUWL16ss846S3feeWekhwUAiDGxXnAs9mc+evUyMxCB2L7d90zHli3S0KGBvXYAevXqpcGDB2vdunV65pln9Mc//pEOswCAgLRVcCwnJydmZkBiP/iw2QJa+pAknXaauU1IkLze5tvTTgv8uQI0atQoPfHEE7rhhhs0duzYkL4WACD+xEPBse753+7UVCktTRo+XFq+3NympZnjITZkyBD16NFDP//5z0P+WgCA+BMPBce6Z/DhcEhut1RaKuXnm1u32xwPseeee07Tp0+PqX8kAIDoEQ8Fx2J/2aWzkpOb/2yztbwfZF6vV5988olcLpcqKyv129/+NmSvBQCIf7m5ucrJyVFVVZUyMzNjKvCQunPwEUavv/66srOzdcEFF2jdunXq27dvpIcEAIhxsVxwjOAjDMaOHXvK5iAAALqr7rnnAwAARAzBBwAACCuCDwAAupMwd3T3heADAIDu5OSO7hEScPDx+uuva8KECUpPT5fNZtOGDRvaPPeuu+6SzWbT448/3oUhAgCALqmulsrLpYoKafVqc+y558z98nLzeBgFnO1y9OhRDRkyRFOnTtVNN93U5nnr16/XW2+9pfT09C4NEAAAdNGgQace++QTU+G7kWWFbTgBBx/XXXedrrvuunbP+de//qUf/OAH2rhxo8aPH9/pwQEAgCB49FHpvvtaHmsMNhITpaKisA4n6HU+vF6vJk2apNmzZ+uiiy7q8Py6ujrV1dU13a+pqQn2kAAA6J4sS3rqKWn+/LbPKS2Vhg0L35gUgg2nixcvVmJiombMmOHX+QsXLlRKSkrTV0ZGRrCHBABA97N/vzR+vDR9uvT559LIkeZ4Y1O6hMjlnAT1lcvLy7VkyRIVFRXJZrP59T3z5s3T4cOHm7727NkTzCFFtbFjx2rWrFmRHgYAIN6sXy9dfLH08sumd9mSJdILL0Sso3trQV12+etf/6oDBw7onHPOaTrW0NCgH/7wh3r88cfldrtP+Z7k5GQlh7CpW7zYtGmTrrrqKn322Wfq169fpIcDAIhGtbXSzJnSihXm/tChJrulcRuE2y0lJZmGqnl5Un19SBurtiWowcekSZM0bty4FsdycnI0adIk3XnnncF8qaDweDyqrKxUVlZWzDbnAQAEJm5/97/xhjRpkvTPf5rg4r77pAULWgYXYezo3p6Al12OHDmi7du3a/v27ZKkXbt2afv27dq9e7fOPPNMXXzxxS2+evToobS0NJ1//vnBHnuXuFwuOZ1OZWdny+l0yuVyhfT1jh49qsmTJ+v000/XwIED9X//938tHl+1apVGjBihPn36KC0tTf/93/+tAwcOSJLcbreuuuoqSdIZZ5whm82mO+64Q5L0yiuv6IorrlC/fv105pln6j//8z/10UcfhfRaACBWhft3f1gcPy795CfSlVeawMPplDZtkhYtilhw0ZGAg4+ysjJdcskluuSSSyRJ9957ry655BI98MADQR9cqHg8HuXl5TV1mvV6vcrPz5fH4wnZa86ePVubN2/Wb3/7W/3pT3/Spk2bVFFR0fT48ePH9dBDD+ndd9/Vhg0b5Ha7mwKMjIwMrVu3TpL0wQcfaN++fVqyZIkkE9Tce++9Kisr02uvvaaEhAR961vfoosuALQSid/9Iff++9Lo0dLDD0ter5n5ePdd6ZvfjPTI2hXwssvYsWNlBVCIxNc+j0irrKw85cO5oaFBVVVVIZmCO3LkiFwul1avXq2rr75akrRy5coWrzV16tSmP3/lK1/R0qVLdemll+rIkSM6/fTT1b9/f0lSampqiz0fN998c4vXevbZZ3XWWWfpvffe08UXXxz0awGAWBXu3/0hZVnS009LP/qRyWQ54wypoEC69dZIj8wv3bK3S1ZWlhJapRjZ7XZlZmaG5PU++ugj1dfXa2RjmpOk/v37t1iKKi8v14QJE3TOOeeoT58+GjNmjCRp9+7d7T53ZWWlbrvtNn3lK19R3759NehEFbuOvg8Auptw/+4PmcYU2mnTTOBxzTXSjh0xE3hI3TT4cDgcKiwslN1ul2T+8RUUFEQs8j169KhycnLUt29frVmzRtu2bdP69eslSfX19e1+74QJE/Tpp5/qmWeeUWlpqUpLS/36PgDobqLtd3+n+EqhfeUV6eyzIz2ygAS9wmmsyM3NVU5OjqqqqpSZmRnSf3xf/epX1aNHD5WWljalIX/22Wf68MMPNWbMGL3//vs6dOiQFi1a1FRkraxVq+OkpCRJZoqw0aFDh/TBBx/omWee0ZVXXilJ2rJlS8iuAwBiXTh/9wdVRym0MabbBh+SiYLD8Q/v9NNPV25urmbPnq0zzzxTqamp+vGPf9w0/XfOOecoKSlJTzzxhO666y79/e9/10MPPdTiOZxOp2w2m37/+9/r+uuv12mnnaYzzjhDZ555pgoLCzVw4EDt3r1bc+fODfn1AEAsC9fv/k4rKzNpso8+Ko0Y4V8KbYzplssukfDzn/9cV155pSZMmKBx48bpiiuu0PAT3QTPOussFRUVae3atbrwwgu1aNEi/eIXv2jx/WeffbYWLFiguXPnasCAAZo+fboSEhL03HPPqby8XBdffLHuuece/fznP4/E5QEAgqW4WCopkVaujLkUWn/ZrEBSV8KgpqZGKSkpOnz4sPr27dvisS+++EK7du3Sueeeq549e0ZohLGPv0cAiDLV1dLBg2Zm47rrpAMHTLfZL780j3/rW2bJJSUlsuNsR3uf361162UXAACiwolMxRYaAw/JbDR96aWwDSfUWHYBACDSVq82Mx2+JCaax+MIMx8AAERar15S797S4cOnPlZaKg0b1uWXiKaeNsx8AAAQKbW10tSp0k03NQceNpu5TQjeR3S09bSJyeAjyvbIxhz+/gAgCmzdKg0ZYjaS2mzS3XdLAwaY9Nrly6Xhw6W0NCk1tUsvE409bWIq+OjRo4ck6dixYxEeSWxr/Ptr/PsEAIRRYxfab35T2rWrOYX2qadM1ktpqZSfL8+6dSopKlJXQ4S2etq8+eabXXzmzoupPR92u139+vVrajXfq1cv2Rqnp9Ahy7J07NgxHThwQP369WsqMQwACJP335duv10qLzf3J02SnniiOYX2RP0Ol8vVNFuRkJCgwsJC5ebmduolG3vatA5AvvOd76impqbTz9sVMVXnQzIfoPv379e///3v8A8uTvTr109paWkEbgAQLpZlZjZmz+6wC63H45HT6WwRLNjtdrnd7k5vFHW5XMrPz2/RoiMYz3uyuK7zYbPZNHDgQKWmpur48eORHk7M6dGjBzMeABBk7WaS7NtnNpW+8oq5f801Zp9HG83g2lomqaqq6nSQkJubqz59+ug73/lOUJ+3s2Iu+Ghkt9v5EAUARFy7SyTr10vf+5506JBZUnn0UWn69HYzWXwtk9jtdmVmZnZpnJdffnlInrczYmrDKQAA0aTNTJL3329OoT10yHShLS+XZszoMIXW4XCosLCw6T/YdrtdBQUFXZ6dCNXzdkbM7fkAACBalJSUKDs7+9TjaWkau39/cxfa//1fKSkpoOf2eDyqqqpSZmZmUAOEUD1vXO/5AAAgWvhcIpGUuX+/SaEtLjYptZ3gcDhCMisRqucNBMsuAAB0UtNSxomlFLukAkmOSZOkd9/tdOAR7wg+AAAIRFmZlJ1tbi1LuV98IXdSkkokuVNSlPvCC2bGo7F2B07BsgsAAIEoLpZKSkwJ9H/9S3rlFTkkOTpIoUUzgg8AADpSXS0dPGg2kD7/vDn27LOmeFiPHtL990sPPBDUZnDxjGwXAAA64k9F6Oj6OA27QD6/CdEAAOjI6tVSW4UtExPN4/Abyy4AALTn+HFp506pVcnzJqWl0rBh4R1TjCP4AACgLa270EpmCcayzP6OtgIStItlFwAAWrMsadkyM6NRXm660D79tJSWJo0YYTJdhg8391NTIz3amMOGUwAATtZeF9q6OlMmvXH2o77eNIwDG04BAPHH4/GopKREHo8ndC+yfr00eLD0yivyJCWpZPp0eX71q+baHcnJzZkvNhuBRycRfAAAop7L5ZLT6VR2dracTqdcLldwX6C2tkUXWldGhpxffqnsJ5+U89xzg/963RzLLgCAqObxeOR0Ols2b7Pb5Xa7g9MgbetWadIkadcuyWaT56675CwoCN3rxSmWXQAAcaOysrJFICBJDQ0Nqqqq6toTHz8u/eQnpvnbrl2mC+2mTaq89dbQvB6akGoLAIhqPtvW2+3KzMzs/JO2TqGdPFlaulRKSVGWxxP810MLzHwAAKJaU9v6ExVG7Xa7CgoKOrcE4iuF9oUXpJUrm7rQBvX14BN7PgAAMcHj8aiqqkqZmZmdCwTaS6ENxet1M4F8frPsAgCICQ6HoykI8Hg8qqysVFZWln+Bwfr10ve+Jx06ZNJjH31Umj693S60J78egotlFwBAWASrTodfabdlZVJ2trR5c4sUWg0dapZbZsxoN/BAaPE3DwAIuWDV6fB4PMrLy2vaDOr1epWfn39qQFNcLJWUSP/1X2ZpxWaT5swxTeAuuqirl4MuYtkFABBSbQUMOTk5AS9rtJd262hokA4elBoapMbgpqbG9F/53/+Vrr3WlEZHxBF8AABCqt2AIcDgo92024wM39+0f7+Ul2f+fCLHIuA9Iwgqll0AACHVGDCcrLN1M9pMgz37bGnKlLa/MTFRWr1aUhhKtaNDpNoCAELO5XIpPz9fDQ0NTQFDbm5up5+vRRqs3d4yhdaX8nJp2LDQl2rvxiivDgCIKrm5uXK73SopKZHb7e5S4CGZGZCxY8fKsW1bUxdaJSdLP/yhOaFxpqXVjEvISrUjIOz5AACERVDrZtTWSjNnmkwWyaTQrl5tqpSuWWP2f+Tmmo2ne/ZIqamSQlSqHQEj+AAAxJZWXWh1330mm6Uxk8XtNn+22cxG0/p6Myui5j0jrZeAWHIJL/Z8AABiw/Hj0oIF0sKFktdrutAWF5uutAGidHrwUV4dABBfWnehnTRJeuKJpmZwgaJ0emSx4RQAEL3a6kJbXNzpwAORx8wHACA6BdiFFrGDmQ8AQPRZv75lCu2SJebPBB5xgZkPAED0aCuFlmZwcYWZDwBAdNi6VRoyhC603UDAwcfrr7+uCRMmKD09XTabTRs2bGh67Pjx45ozZ44GDx6s3r17Kz09XZMnT9bevXuDOWYAQKwrK5Oys83t8ePST35iUmZ37TIptJs2SYsW0YU2TgUcfBw9elRDhgzRsmXLTnns2LFjqqio0E9/+lNVVFTopZde0gcffKD/+q//CspgAQBxorhYKimRli6VRo+WHn7Y1O6YNEl6991O1e5A7OhSkTGbzab169frxhtvbPOcbdu26bLLLlN1dbXOOeecDp+TImMAEKeqq6WDB82SynXXSQcOND/Wp4+Z6fj+9yM3PnRJVBUZO3z4sGw2m/r16+fz8bq6OtXV1TXdr6mpCfWQAACRMGhQ24/V1krTphF8dBMh3XD6xRdfaM6cObrtttvajIIWLlyolJSUpq+MjIxQDgkAECmrV0t2u+/HEhPN4+gWQhZ8HD9+XN/+9rdlWZaefvrpNs+bN2+eDh8+3PS1Z8+eUA0JABAptbXSa69JDQ2+Hy8tlSZODO+YEDEhWXZpDDyqq6v1l7/8pd21n+TkZCWf6DYIAIhDJ3ehbZSQYDaYNt6iWwn6zEdj4FFZWak///nPOvPMM4P9EgCAWOArhXbtWiktTRo+XFq+3NympUmpqZEeLcIo4JmPI0eOqKqqqun+rl27tH37dvXv318DBw7ULbfcooqKCv3+979XQ0OD9u/fL0nq37+/ksjXBoDuob0utBMmmPodNpuUlyfV15sS6ug2Ak613bRpk6666qpTjk+ZMkUPPvigzj33XJ/fV1JSorFjx3b4/KTaAkAMsyzpqaek2bOlzz83XWgLCqRbb/V5usfjUWVlpbKysmhxH+NCmmo7duxYtRevdKFsCAAglgXYhdblcikvL09er1cJCQkqLCxUbm5uGAeMSOlSkbFQYOYDAGLQ+vXS974nHTpkllAefVSaPt1sKPXB4/HI6XTKe9JmU7vdLrfbzQxIjArk85vGcgCAzqutNbMdN91kAo+hQ80+jxkz2gw8JKmysrJF4CFJDQ0NLfYUIn4RfAAAOqcLXWizsrKU0Co4sdvtyszMDNVoEUUIPgAAgamvl37846YUWs/ZZ6vkscfkmT7d7y60DodDhYWFsp+oeGq321VQUMCSSzfBng8AgP9apdC6Ro9WXmlppzeNejweVVVVKTMzk8AjxgXy+U3wAQDoWOsU2v795XnkETm//302jUISG04BAF1VViZlZ5vbffuk66832Suff25SaHfsUOV55/ncNLp27Vp5PJ4IDRyxgOADAHCq4mKppESaP18aPNjU7ujZU1qyxPw5Pd3nplFJuvfee+V0OuVyuSIwcMQCll0AAEZ1tXTwoMlc+X//T/rkk+bHzjtPWrZMGjeuxbe4XC7l5+erwUe3WpZguheWXQAAgRs0SBoxwjR7OznwkKQPPzTLLa3k5ubK7XbrscceO+Ux6nagLQQfANCNeTwelZSUmD0aRUVtFwZLTJRWr/b5kMPh0K233krdDviN4AMAuimXyyWn06ns7GyzR+OBB6RWG0iblJZKEye2+Vy+6nYsWrRIlZWVbD7FKdjzAQDdkM/eKpLcp58ux5EjZgbE622+LS+Xhg3z63mrqqq0bds2zZ07l6Zx3Qh7PgAA7fLZW0VS1S9+IaWlmX0fy5eb27Q0KTXVr+d1OBzKzMxsCjwkyev1Kj8/nxkQNEmM9AAAAOGXVVWlBEknhx92u12Z48dLd9xhyqTbbFJenimnnpzs93O31zSOzBdIzHwAQPdSUyPdeacceXkqlFlqkVr1VklONoGHZG4DCDwkmsahYwQfANBdbN1qWt4XFUk2m3LnzpX7o49UUlIit9sdtD0ZNI1DR9hwCgDxrr5eWrBAWrTIbB51OqVVq6Qrrwzpy9I0rnsJ5PObPR8AEM/ef9+kyFZUmPuTJ0tLl0opKSF/aYfDQdABn1h2AYB4ZFmmHPqwYSbw6N9fWrtWWrkyLIEH0B5mPgAg3uzbJ02dahrASaYselGRlJ4e0WEBjZj5AIB48tJLbXahBaIFMx8AEA9qa6WZM+VZsUKVkrIuvFCOtWulCy+M9MiAUzDzAQCxqKxMys42t1u3SkOGyLVihZySsiU5339frjffjPQoAZ8IPgAgFhUXSyUl0t13S9/8pjy7dilPzRVLg1HSvEXHWyCICD4AIFZUV5sGbxUV0po15lhZmeT1qnLkSLXuR9tY0rwzTul463IRjCBoKDIGALGiseS5Dx5JTp3aq8Xtdgdca8NXx1ubzSabzUaXWrSJrrYAEI+efLLNAMSRmKjC//mfoJQ099UYzrIsutQiaMh2AYBYsH69NH++KR7mS2mpcocNU878+V0uad7YGK51AHIyutSiK5j5AIBoVltrCobddJN06JB03nnmeGPX2FbdYx0Oh8aOHduloKB1Y7iEhATZWs240KUWXUHwAQB+CvuGyxMptFqxwiy3zJ0rvfyylJYmDR8uLV9ubtPSpNTUoL50bm6u3G63SkpKVF1drWeeeYYutQgaNpwCgB9cLpfy8vLCs+GyVRdaz9lnq/JHP1LWLbeYD/y6OikpyQQklmXOT04OzVhOQpdatCeQz2+CDwDogK/sj85mknTo/fel2283KbWSXKNHK6+0lCwTRD2yXQAgiHxlf3SlhoZPJ3ehLS+XzjhDnuXLmwIPiSwTxA+CDwDoQGP2x8mCuuFy3z7p+uul6dOlzz83XWh37FDleeeFPugBIoDgAwA60Dr7I6gbLtvqQnv22aEPeoAIIfgAAD+cnP3hdru7vu+iMYX25ptNCu3QoWa5ZcaMpvTZkAY9QASx4RQAwm3rVmnSJGnXLpOxMmeOyW5JSvJ5ejizTDwejyorK5WVlUWQg4Cw4RQAgiwoNT7q66Uf/1j65jdN4OF0Sps2SQsXthl4SMEpHOYPX83kgFAg+ACADnT6Q7msTMrONrfvvy9dfrn0yCOS1ytNniy9+64JRKKAx+NpqmMikVmD0CL4AIB2dOlDubhYKimRZs9uTqHt319au1ZauVJKSQnx6P0XlnRi4AQaywFAO9r7UPa1DOJ56y1Vvvuuss45R47f/MYc3LTJ3I4cKS1dKl12WYhHHThfzeTIrEGoMPMBAO0IJN3V5XLJOXq0su+6S87rr5fr4MGWJ5SWmgAkCpFZg3Ai+ACAdvj7ody0PHPivldSvqQWizOJidLq1WEYdecEPZ0YaAPLLgDQgdzcXOXk5LSb7upzeUZSlaSms0tLzd6PKOZwOJjtQMgRfACAH9r9UK6vV9ZLLylB0snhh11SpmSKhrUKTIDujGUXAOiK99+XRo+W48knVSjJbrNJMoFHwaBBcixfLg0fLqWlSampER0qEC2Y+QCAzrAs6amnTBrt559L/fsrt6BAOaNGmeWZjAw5vvIVU8E0L88UGEtOjvSogahA8AEAgdq3z/RleeUVc//aa6UVK6T0dDmkU5dnbDYCD+AkLLsAQCBad6FdulR6+WUpPT3SIwNiBjMfAOCPmhpp5kypqMjcHzpUWrNGuvDCSI4KiEnMfACIW0FpBieZLrRDh5rAw2aT5s41abMEHkCnEHwAiEtB6dDqqwvt5s0ddqEF0L6Ag4/XX39dEyZMUHp6umw2mzZs2NDiccuy9MADD2jgwIE67bTTNG7cOFVWVgZrvADQoaB0aN25Uxo9urkL7ZQp0t/+Jl15ZYhGDXQfAQcfR48e1ZAhQ7Rs2TKfjz/66KNaunSpli9frtLSUvXu3Vs5OTn64osvujxYAPBHlzq0Wpa0bJmpRFpR0dyFtqhI6ts3NAMGupmAN5xed911uu6663w+ZlmWHn/8cf3kJz/RDTfcIEkqLi7WgAEDtGHDBn33u9/t2mgBwA+d7tDaTgotgOAJ6p6PXbt2af/+/Ro3blzTsZSUFI0cOVJvvvmmz++pq6tTTU1Niy8A8EdbG0r97tBaViZlZ5tbUmiBsAlq8LF//35J0oABA1ocHzBgQNNjrS1cuFApKSlNXxkZGcEcEoA41dGGUr86tBYXSyUl0uTJ0s03S4cOmayW8nLpBz8wPVkABF3Ef7LmzZunw4cPN33t2bMn0kMCEOX83VDqcDg0duzYljMe1dUmuKioaG5vv3Onub3jDumFF0ihBUIsqEXG0tLSJEkff/yxBg4c2HT8448/1tChQ31+T3JyspIpOwwgAO1tKO2wHfygQW0/VlRkvizL77F4PB5VVlYqKyuLVvSAn4I683HuuecqLS1Nr732WtOxmpoalZaWavTo0cF8KQDdWOOG0pP5taFUkhYvbvuxxMTm2RA/BKWWCNANBRx8HDlyRNu3b9f27dslmU2m27dv1+7du2Wz2TRr1iz97Gc/0+9+9zvt2LFDkydPVnp6um688cYgDx1Ad+X3htKTWZb05JPS/Pltn1NaKk2c6NcYglJLBOimAl52KSsr01VXXdV0/95775UkTZkyRUVFRbrvvvt09OhR5eXl6d///reuuOIKvfLKK+rZs2fwRg2g28vNzVVOTo5pX5+Z2X7g0TqFdtQo6a23zIZSr7f5NgBdWvoBujmbZQWwuBkGNTU1SklJ0eHDh9WXgj4Auuqll6S8PJPJ0rOn9Oij0g03SCNHShkZUm6u5HJJe/ZI27ZJfgYOHo9HTqfzlFoibreb4APdUiCf33S1BRCfOupC63ab/iw2mwlO6uulADa/Ny795Ofnq6Ghwb+lHwCSmPkAEI+2bJEmTTIBhs0mzZkjLVgQkmZwHo/Hv6UfIM4x8wGge6qvN0HGokVmD4fTKa1aFdJmcA6Hg6ADCBDBB4D4sHOndPvtpniYZLrQLl1KMzggCkW8wikAdAldaIGYQ/ABICBtNXOLiH37pOuvl6ZPl774wnSh3bFDuuUWv58iqq4H6CYIPgD4LaoqegahC21UXQ/QjZDtAsAvUVPXoqMUWj9FzfUAcSKQz29mPgD4pb2KniFTViZlZ5tbSdq61QQbRUUmhXbuXFMSvRNdaCNyPQAkke0CwE+NzdxazxT41cyts4qLpZISE2ysXx/UFNqIXA8AScx8APBTp5q5dUZ1tVRebjJXnn/eHCsokB55xAQeN90kvftul2t3tL6ehIQE3XPPPV0dPQA/sOcDQEBCXtHTZuv4nCD+2vJ4PFqyZIkee+wxeb1eJSQkqLCwULm5uUF7DaA7COTzm+ADQHRZs0a64w7pyy9PfSwx0SzB+Nn23h9sPAWCgw2nAGLXaadJvXv7fqy0NKiBh8TGUyASCD4ARIfaWmnqVOnmm6XDh82xxiWYhND9qmrceHoyNp4CoUXwASAs2q0kunWrNGSItGKFCTi+/31pwABpxAhp+XJp+HApLU1KTQ36uMK2kRZAE/Z8AAg5l8ulvLy8Uzd0tteFtq5OSkoywYhlmXOTk0M2xpBvpAXiHBtOAUSNNjd0/vnPcvzwh81daCdPNiXSU1IiNFIAXRHI5zdFxgCEVJsbOnNy5KivN11oCwoCagYHILYRfAAIKZ+VRCVl1tdL11xjUmcDaAYHIPax4RRASDVt6DyRUWKXVJCYKMeSJaYjLYEH0O0w8wEgtGpqlLtli3K8XlVJyrzwQjnWru1UMzgA8YHgA0DobN0qTZok7dolh80mx5w5JrslKSnSIwMQQQQfAIKvvRRaAN0ewQeA4Hr/fVMCnRRaAG1gwymAzisrk7Kzza1lScuWScOGmcCjf39p7Vpp5UoCDwAtMPMBoPOKi6WSElMC/V//MtkrknTttaZUOpksAHwg+AAQmOpq6eBBU/b8+efNsWefNTMfSUnS/fdLP/1pSJvBAYhtlFcHEJjGTrPtia5fKwDCIJDPb/5rAiAwq1dLJzrAniIx0TwOAO1g2QWA/+rrpffeM+mzvpSWmg2nANAOgg8A/tm5U7r99uYUWsns6/B6m28BwA8suwBon68U2uXLpbQ0afhw8+fhw8391NRIjxZADGDDKYC27dsnTZ3qO4W2rs5kt9hsJkCpr5eSkyM7XgARw4ZTAF330kvS4MEm8OjZ01Qpffnl5todycnNmS82G4EHAL+x5wNASzU10syZUlGRuT90qLRmDV1oAQQNMx8Amm3ZIg0ZYgIPm02aO9dksBB4AAgiZj4A0IUWQFgRfADdXesU2ilTzP4ONnwDCBGWXYDuyrKkJ588tQttURGBB4CQYuYD6I7aS6EFgBBj5gPobjpKoQWAECP4AOJdWZmUnS1t2iTdead0883SoUPSJZdI5eXSD35gyqMDQJiw7ALEu+JiqaTEBCG1tc0ptA8+aCqUAkCYEXwA8ai6Wjp4UPryS8nlMsdqa03/lYcekq65hsADQMQQfADxaNAg38f375e+9z3z5+hq6wSgG2GhF4g3liVNntz244mJ0urV4RsPALTCzAcQT/buNSm0Gze2fU5pqantAQARwswHEC8aU2g3bjQptLNnm+ONmSxktACIEvw2AmJdTU1zCu2nnzan0M6YYTaYDh8uLV9ubtPSpNTUSI8YQDfHsgsQy7ZskSZNktxu3ym0brf5s80m5eWZBnLJyREcMAAQfACxqb7eBBmLF7ffhfbkQMNmI/AAEBWCvuzS0NCgn/70pzr33HN12mmn6atf/aoeeughWaT1AcGxc6c0erS0cKEJPKZMkf72t1MDDwCIUkGf+Vi8eLGefvpprVy5UhdddJHKysp05513KiUlRTNmzAj2ywFRxePxqLKyUllZWXI4HMF9csuSli0zG0m/+MJ0oS0sNHs9ACCGBD34eOONN3TDDTdo/PjxkqRBgwbpN7/5jd5+++1gvxQQVVwul/Ly8uT1epWQkKDCwkLl5uYG58lbp9Dm5EjPPkszOAAxKejLLpdffrlee+01ffjhh5Kkd999V1u2bNF1110X7JcCoobH42kKPCTJ6/UqPz9fHo+n60++bl3LFNonnqALLYCYFvSZj7lz56qmpkYXXHCB7Ha7Ghoa9PDDD2vixIk+z6+rq1NdXV3T/ZqammAPCQi5ysrKpsCjUUNDg6qqqjq//FJTI82cKRUVmfuXXGIqk154YdcGCwARFvSZjxdeeEFr1qzRr3/9a1VUVGjlypX6xS9+oZUrV/o8f+HChUpJSWn6ysjICPaQgJDLyspSQqsiXna7XZmZmZ17wi1bpCFDTOBhs0nz5klvvUXgASAuBD34mD17tubOnavvfve7Gjx4sCZNmqR77rlHCxcu9Hn+vHnzdPjw4aavPXv2BHtIQMg5HA4VFhbKbrdLMoFHQUGBf7MeZWVSdra5ra+X7r9fGjPG1OhwOqXNm6VHHqELLYC4EfRll2PHjvn8H2DrKelGycnJSqb2AOJAbm6ucnJyVFVVpczMTP+XW4qLpZISackS6b33pIoKc3zKFGnpUqlv39ANGgAiIOjBx4QJE/Twww/rnHPO0UUXXaR33nlHjz32mKZOnRrslwKijsPh8C/oqK6WDh40SyrPP2+ONXaa7dvXFA+7667QDRQAIshmBbn6V21trX76059q/fr1OnDggNLT03XbbbfpgQceUJIf08Y1NTVKSUnR4cOH1Zf/8SHONNUByc5WhyEKhfkAxJBAPr+DHnx0FcEH4lWLOiA2mwol5fr68UtMNBtN28gQA4BoRPABRBmPxyOn09li75Ndkls6dQakvFwaNix8gwOAIAjk8zvo2S5AvPN4PCopKQmogJjPOiCSqiSpcYN2Aj+OALoHftsBAXC5XHI6ncrOzpbT6ZTL5fLr+7KcTiXYbC2O2RMSlHnWWdLw4dLy5eY2LU1KTQ3F0AEgarDsAvjJ59KJ3S63291+hsvOndLtt8tVUaF8mRmPxjogubffbup32Gxmg2l9PW3vAcSkQD6/g55qC0STYHaZDbiEeqsutLn9+ytn4UJVnXee7zogNhuBB4BugeADcSvYXWYbS6i3nvnwWUK9dRfaa6+VVqyQIz294xRbAIhz7PlAXApFl1m/S6i/9FLLLrRLl9KFFgBOwswH4lJIusyqgxLqdKEFAL8QfCAuBbREEiCfJdS3bJEmTTLN4Gw2ae5c6cEHaQYHAD6w7IK41KUus4Gor5d+/GO60AJAAEi1RVzzeDwBd5n1O0PmRAotXWgBgAqnQBOHw6GxY8f6HXj4VUTMsqQnnzQl0CsqpP79pbVrzV4PAg8A6BDBB3BCmxkyf/iDlJ0tlZVJ+/ZJ118v/eAH0hdfmBTaHTukW2455bkCLcEOAN0FwQdwQpsZMi6XVFIizZ9vUmhfeaXdFNrOlmAHgO6CPR/ACT7LpyckyN2vnxyfftp84nnnmcql48b59xz+lGAHgBjHng+gE07JkJFU4PW2DDwk6cMPpWuu8fkc7dUXAQAYBB/ASXJzc+V2u1UyZ47cNpt8FmNPTDTFw3xorC9ysmDVFwGAeEHwAbTiqK3V2FdflaOtFcnSUmniRN/fG676IgAQw6hwCjRq1YVWffuakukJCZLX23zbgXZLsAMAmPkAJPlOoX3tNSktTRo+XFq+3NympUmpqR0+XaD1RQCgO2HmA3jpJSkvTzp0yKTQPvqoNG2amelwu02ZdJvNnFNfLyUnR3rEABDTCD7QffnThfbkQMNmI/AAgCBg2QXd05Yt0pAhJvCw2aR586S33moZeAAAQoKZD3Qv9fXSggXSokVm8+igQdKqVdIVV0R6ZADQbRB8oPtopwut351sAQBdxrIL4p+vLrQvvtjUhZZeLAAQXvR2QXzbu1eaOlXauFEeSZUjRijr6aflGDFCEr1YACBY6O0CSCaFdvBgaeNGuRIT5bTZlF1WJufIkU2zG/RiAYDwI/hAfCgrk7KzzW1NjXTnndLNN0uffirPRRcpz+uV98Qkn9frVX5+vjweD71YACACCD4Qczwej0pKSuTxeJoPFhdLJSXS4sWnpNBWPvZYm7Mb9GIBgPBjzwdiisvlUl5enrxerxISElR4//3K/da35LnmGlV++qmyJDkkaeBAk8lyyy1+7evweDz0YgGALgjk85vgAzHDZxAhaZGkOZK8MlN5hZJyJbPB9C9/UVZWljZu3Kj8/Hw1NDQ0zW7k5uZG4CoAID4F8vlNnQ/EDJ+bQ9UceOjEbb6kT202zZXkzc42MySFhXK73cxuAEAUYM8HYoavzaEJag48GjUFJK02mO7bt09RNtEHAN0SwQciyufm0TY0bQ49EYDYJS1KSDjlH3GCdEqQ0dDQoFGjRlFIDACiAMEHIibgyqI1NcrdskVur1clktwXXaTZf/6zCvv2lf3EKXZJi08//ZQZEklNSzYnp9oCAMKP4AMR4fF4mrJWJD8CgpO60DoSEjR23jw5Kiqkq65S7oEDcu/erZKSErl379aPDh5skT7rKxChkBgARA4bThER7VUWbbEZtL5eevBBU7+jrS60yclyZGTIkZHRdCg3N1c5OTmqqqpS7969NWrUqFNSbSkkBgCRwcwHIqKtyqK9e/du3gOyc6c0erS0cKEJPO64Q3r33ZaBRzscDofGjh2rSy+9lEJiABBFqPOBiHG5XC1qb9x+++1atWqVKSBms6nQblful1+aLrSFhaZcehdQSAwAQociY4gZjQGBz6URSe4xY+T49a+l9PTIDRIA0CG62iJmNC6NHDlyxGcBsar58wk8ACDOEHwg8mpqlPXUU6f8Y7Tb7crMyorIkAAAoUPwgcg6kULrePFFFUqy22yS2BQKAPGMVFtEho8U2txVq5QzaBCbQgEgzhF8IDzKyqT77pMefVTq3Vu6/XaposI8dscd0pIlUt++ckgEHQAQ5wg+EB7FxVJJifSjH0mlpdIXX3Qphdbj8aiyslJZWVkEKwAQY9jzgdCprpbKy80Mx69/bY5t3mwCj1GjpJdf7lTgEXBPGABAVKHOB0LnxObRdgX4z8/j8cjpdJ5SKt3tdjMDAgARRJ0PRIdnnmk7AElMlFavDvgp2+sJAwCIDez5QGhs2SI9/HDbMxulpdKwYQE/bWNPGJrEAUDsYuYDwVVfL91/vzRmjOR2SwMHmuONTeR8tLcPhMPhoEkcAMQ4Zj4QPDt3tkyhnTJFmjNHys6WMjKk3FzJ5ZL27JFSUzv9Mrm5ucrJyaEeCADEqJBsOP3Xv/6lOXPm6OWXX9axY8eUmZmpFStWaMSIER1+LxtOY5BlScuWSbNn+06hrauTkpLM/g/LMrMjycmRHTMAIKgC+fwO+szHZ599pm984xu66qqr9PLLL+uss85SZWWlzjjjjGC/FKLB3r3S1KnSxo3mfk6O9OyzLZvBnRxo2GwEHgDQzQU9+Fi8eLEyMjK0YsWKpmPnnntusF8G0WDdOikvT/r0U6lnT+nnP5emTfMvxRYA0G0FfcPp7373O40YMUK33nqrUlNTdckll+iZZ54J9ssgkmpqpDvvlG65xQQel1xiiolNn07gAQDoUNCDj3/+8596+umnlZWVpY0bN+ruu+/WjBkztHLlSp/n19XVqaampsUXotiJLrQqKjKBxrx50ltvSRdeGOmRAQBiRNA3nCYlJWnEiBF64403mo7NmDFD27Zt05tvvnnK+Q8++KAWLFhwynE2nEaZ1l1onU5p1SrpyisjPTIAQBSIaIXTgQMH6sJW/wv+2te+pt27d/s8f968eTp8+HDT1549e4I9JHTVzp3S6NHSwoUm8JgyRfrb3wg8AACdEvQNp9/4xjf0wQcftDj24Ycfyul0+jw/OTlZyWQ/RKeOUmgBAOiEoAcf99xzjy6//HI98sgj+va3v623335bhYWFKiwsDPZLIZRap9Bee620YkXLFFoAADoh6Msul156qdavX6/f/OY3uvjii/XQQw/p8ccf18SJE4P9UgimsjJTibSsTHrpJWnwYBN49OwpLV0qvfwygQcAIChCUuG0K6hwGiEzZkhPPCF97Wtmj4dkUmhXryaTBQDQoYhWOEUMqa6WDh40KbON7e0bA4877jAN4rKyIjY8AEB8YuajO/OnIFh0/fMAAESpiKbaIoYsXtz2Y4mJzbMhAAAEEcsu3VFjCu38+W2fU1oqDRsWvjEBALoNZj66m717peuuk37wA1O7Y9QoczwhoeUtAAAhwidNd7Ju3akptM8/L6WlScOHS8uXm9u0NCk1NdKjBQDEKZZduoOaGpNK29jcr3UKrdstJSWZDah5eaaPC1VnAQAhwsxHvGvsQrtyZdtdaJOTmzNfbDYCDwBASDHzEa9ad6EdNMh0ob3iikiPDADQzRF8xKOdO6Xbb5cqKsz9O+6QliyRqJsCAIgCLLvEE8syJdKHDTOBR//+0osvmoZwBB4AgCjBzEe82LtXuvNO6U9/MvdzcqRnn6UZHAAg6jDzEQ8aU2j/9CeTQvvEE3ShBQBELWY+YpmvFNo1a0xnWgAAohQzH7GirEzKzja3kvTXvzan0CYkNKfQEngAAKIcMx+xorhYKimRiorMMsvixWaDKSm0AIAYQ/ARzaqrpYMHTeGv5583xwoKpC+/NH++5RbJ5SKTBQAQU2yWZVmRHsTJampqlJKSosOHD6tvd/9Qbaw62p7oevsAAN1UIJ/f7PmIZqtXS4ltTE4lJprHAQCIMSy7RLPkZKlXL5PV0lppqSkmBgBAjGHmIxrV1JiS6Lfe2hx4NC7BJPCWAQBiG59k0aZ1Cu20adKAAdKIEdLy5dLw4VJampSa2qmn93g8KikpkcfjCfLAAQDwD8FHtKivN7U6xoyR3G6TQrt5s/TkkybrpbRUys83t2635HAE/BIul0tOp1PZ2dlyOp1yuVzBvgoAADpEtks0eO8904X2nXfM/RB0ofV4PHI6nfJ6vU3H7Ha73G63HJ0IZAAAOBnZLrGisQvt8OEm8AhhF9rKysoWgYckNTQ0qKqqKqivAwBAR8h2iZQwd6HNyspSQkLCKTMfmZmZIXk9AADawsxHJLz4Yti70DocDhUWFsput0sygUdBQQFLLgCAsGPPRzhFQRdaj8ejqqoqZWZmEngAAIImkM9vll3C5a9/lSZPNpkqCQnSnDnSgw9KSUlhHYbD4SDoAABEFMFHqNXXS/Pn04UWAIATCD5CaedOaeLEkKbQAgAQa9hwGkxlZVJ2trRtm9lEOmxYyFNoAQCINcx8BFNxsVRSIt1yi7R7tzkW4hRaAABiDcFHV1VXSwcPmsZvxcXm2O7dZiPpzJnS979P4AEAwElItW2Hx+NRZWWlsrKy2s4Qaew2257o+isGACDoKK8eBH43YfvJT9p+ksREafVqOskCAHASgg8fPB6P8vLymkqRe71e5efntwwe6uul+++XHn647ScqLZXriy/oJAsAwEkIPnzosAnbzp3S6NHSwoVmSWXCBHM8IaHFrefjjzsOYgAA6GYIPnxobMJ2MrvdrsyvfrU5hbaiojmF9qmnpLQ00512+XJzm5amytpaOskCANAK2S4+NDZhy8/PV0NDg2nCtnixHN/7nrRxozmpdQqt220yXGw2KS9Pqq9X1ief0EkWAIBWmPloQ25urtxut0pKSuR+6inlPvKICTza6kKbnNyc+WKzScnJdJIFAMAHUm3bH0xQutDSSRYAEO/oahsMW7ZIkyaZ5RSbTZo7t9NdaOkkCwBAM4KP1urrTZCxeLHk9dKFFgCAICP4ONnOndLtt5tMFokutAAAhAAbTiVTq8NXCi1daAEACDpmPvbulaZObTuFFgAABFX3mvkoK5Oys82tJK1bJw0e3H4KLQAACKruNfNRXCyVlEgul/Tkk11OoQUAAIGL/+Cjulo6eNCkyz7/vDlWWGgyWSTp+9+XfvnLTqXQAgCAwMV/kbHGqqPtia6/AgAAYk4gn9/xv+dj9WopsY0JnsRE8zgAAAib+F92mTjR7OUYPvzUx0pLTXotAAAIm5DPfCxatEg2m02zZs0K9Ut1LCGh5S0AAAi7kH4Kb9u2TQUFBfr6178eypfpWGqqlJZmZj+WLze3aWnmOAAACKuQBR9HjhzRxIkT9cwzz+iMM84I1cv4x+EwDeJKS6X8fHPrdpvjAAAgrEIWfEybNk3jx4/XuHHj2j2vrq5ONTU1Lb5CIjm5OfPFZjP3AQBA2IVkw+lzzz2niooKbdu2rcNzFy5cqAULFoRiGAAAIAoFfeZjz549mjlzptasWaOePXt2eP68efN0+PDhpq89e/YEe0gAACCKBL3I2IYNG/Stb31Ldru96VhDQ4NsNpsSEhJUV1fX4rHWgl5kDAAAhFwgn99BX3a5+uqrtWPHjhbH7rzzTl1wwQWaM2dOu4EHAACIf0EPPvr06aOLL764xbHevXvrzDPPPOU4AADofqi2BQAAwios5dU3bdoUjpcBAAAxgJkPAAAQVgQfAAAgrAg+AABAWIVlz0cgGsuOhKzMOgAACLrGz21/yodFXfBRW1srScrIyIjwSAAAQKBqa2uVkpLS7jlBr3DaVV6vV3v37lWfPn1ka2wEFyQ1NTXKyMjQnj174rJ6arxfnxT/18j1xb54v0auL/aF6hoty1Jtba3S09OVkND+ro6om/lISEiQI8St7vv27Ru3/6ik+L8+Kf6vkeuLffF+jVxf7AvFNXY049GIDacAACCsCD4AAEBYdavgIzk5WfPnz1dycnKkhxIS8X59UvxfI9cX++L9Grm+2BcN1xh1G04BAEB861YzHwAAIPIIPgAAQFgRfAAAgLAi+AAAAGEVd8HHsmXLNGjQIPXs2VMjR47U22+/3e75a9eu1QUXXKCePXtq8ODB+uMf/ximkXZOINdXVFQkm83W4qtnz55hHG1gXn/9dU2YMEHp6emy2WzasGFDh9+zadMmDRs2TMnJycrMzFRRUVHIx9kVgV7jpk2bTnkPbTab9u/fH54BB2DhwoW69NJL1adPH6WmpurGG2/UBx980OH3xdLPYGeuMZZ+Dp9++ml9/etfbyo+NXr0aL388svtfk8svX9S4NcYS++fL4sWLZLNZtOsWbPaPS/c72NcBR/PP/+87r33Xs2fP18VFRUaMmSIcnJydODAAZ/nv/HGG7rtttuUm5urd955RzfeeKNuvPFG/f3vfw/zyP0T6PVJpoLdvn37mr6qq6vDOOLAHD16VEOGDNGyZcv8On/Xrl0aP368rrrqKm3fvl2zZs3S//zP/2jjxo0hHmnnBXqNjT744IMW72NqamqIRth5mzdv1rRp0/TWW2/p1Vdf1fHjx3Xttdfq6NGjbX5PrP0MduYapdj5OXQ4HFq0aJHKy8tVVlam7Oxs3XDDDfrHP/7h8/xYe/+kwK9Rip33r7Vt27apoKBAX//619s9LyLvoxVHLrvsMmvatGlN9xsaGqz09HRr4cKFPs//9re/bY0fP77FsZEjR1r5+fkhHWdnBXp9K1assFJSUsI0uuCSZK1fv77dc+677z7roosuanHsO9/5jpWTkxPCkQWPP9dYUlJiSbI+++yzsIwpmA4cOGBJsjZv3tzmObH2M9iaP9cYyz+HlmVZZ5xxhvWrX/3K52Ox/v41au8aY/X9q62ttbKysqxXX33VGjNmjDVz5sw2z43E+xg3Mx/19fUqLy/XuHHjmo4lJCRo3LhxevPNN31+z5tvvtnifEnKyclp8/xI6sz1SdKRI0fkdDqVkZHRYXQfa2Lp/euqoUOHauDAgbrmmmu0devWSA/HL4cPH5Yk9e/fv81zYv099Ocapdj8OWxoaNBzzz2no0ePavTo0T7PifX3z59rlGLz/Zs2bZrGjx9/yvvjSyTex7gJPg4ePKiGhgYNGDCgxfEBAwa0uT6+f//+gM6PpM5c3/nnn69nn31Wv/3tb7V69Wp5vV5dfvnl8ng84RhyyLX1/tXU1Ojzzz+P0KiCa+DAgVq+fLnWrVundevWKSMjQ2PHjlVFRUWkh9Yur9erWbNm6Rvf+IYuvvjiNs+LpZ/B1vy9xlj7OdyxY4dOP/10JScn66677tL69et14YUX+jw3Vt+/QK4x1t4/SXruuedUUVGhhQsX+nV+JN7HqOtqi+AZPXp0i2j+8ssv19e+9jUVFBTooYceiuDI4K/zzz9f559/ftP9yy+/XB999JF++ctfatWqVREcWfumTZumv//979qyZUukhxIy/l5jrP0cnn/++dq+fbsOHz6sF198UVOmTNHmzZvb/HCORYFcY6y9f3v27NHMmTP16quvRvXG2LgJPv7jP/5DdrtdH3/8cYvjH3/8sdLS0nx+T1paWkDnR1Jnrq+1Hj166JJLLlFVVVUohhh2bb1/ffv21WmnnRahUYXeZZddFtUf6tOnT9fvf/97vf7663I4HO2eG0s/gycL5Bpbi/afw6SkJGVmZkqShg8frm3btmnJkiUqKCg45dxYff8CucbWov39Ky8v14EDBzRs2LCmYw0NDXr99df15JNPqq6uTna7vcX3ROJ9jJtll6SkJA0fPlyvvfZa0zGv16vXXnutzbW80aNHtzhfkl599dV21/4ipTPX11pDQ4N27NihgQMHhmqYYRVL718wbd++PSrfQ8uyNH36dK1fv15/+ctfdO6553b4PbH2HnbmGluLtZ9Dr9eruro6n4/F2vvXlvausbVof/+uvvpq7dixQ9u3b2/6GjFihCZOnKjt27efEnhIEXofQ7aVNQKee+45Kzk52SoqKrLee+89Ky8vz+rXr5+1f/9+y7Isa9KkSdbcuXObzt+6dauVmJho/eIXv7B27txpzZ8/3+rRo4e1Y8eOSF1CuwK9vgULFlgbN260PvroI6u8vNz67ne/a/Xs2dP6xz/+EalLaFdtba31zjvvWO+8844lyXrsscesd955x6qurrYsy7Lmzp1rTZo0qen8f/7zn1avXr2s2bNnWzt37rSWLVtm2e1265VXXonUJXQo0Gv85S9/aW3YsMGqrKy0duzYYc2cOdNKSEiw/vznP0fqEtp09913WykpKdamTZusffv2NX0dO3as6ZxY/xnszDXG0s/h3Llzrc2bN1u7du2y/va3v1lz5861bDab9ac//cmyrNh//ywr8GuMpfevLa2zXaLhfYyr4MOyLOuJJ56wzjnnHCspKcm67LLLrLfeeqvpsTFjxlhTpkxpcf4LL7xgnXfeeVZSUpJ10UUXWX/4wx/CPOLABHJ9s2bNajp3wIAB1vXXX29VVFREYNT+aUwrbf3VeE1TpkyxxowZc8r3DB061EpKSrK+8pWvWCtWrAj7uAMR6DUuXrzY+upXv2r17NnT6t+/vzV27FjrL3/5S2QG3wFf1yWpxXsS6z+DnbnGWPo5nDp1quV0Oq2kpCTrrLPOsq6++uqmD2XLiv33z7ICv8ZYev/a0jr4iIb30WZZlhW6eRUAAICW4mbPBwAAiA0EHwAAIKwIPgAAQFgRfAAAgLAi+AAAAGFF8AEAAMKK4AMAAIQVwQcAAAgrgg8AABBWBB8AACCsCD4AAEBYEXwAAICw+v8L8jQQA8bB0AAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "xNew = np.linspace(0,4,20)\n", + "yHat = theta0Hat + theta1Hat*xNew\n", + "plt.plot(xNew, yHat, '-*r', label=\"$\\hat{y}$\")\n", + "plt.plot(x,y,'.k', label=\"data\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Functions for data and OLS" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "def DataGen(xn: float,n: int, disp,theta0=3.9654,theta1=2.5456):\n", + " x = xn*np.random.rand(n, 1)\n", + " #theta0 = 3.9654\n", + " #theta1 = 2.5456\n", + " y = theta0+theta1*x+disp*np.random.randn(n,1)\n", + " return x,y" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [], + "source": [ + "x,y = DataGen(9, 100, 1, 0,1)" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAoUElEQVR4nO3db4xcVf3H8c/OlG4bsrsRtC3NDENhJ0H5D62l1OiwNqwEiBiDYiAWnHQ3pEBrE7VVgQcWlgJqf7RQumQCGKiKMShigJA6gJUipYiByJ9ZYWBvmxY0ulOqLmRmfg/aHbq7M7NzZ++9596571cyD3q7M/d0BvZ85pzvOaetXC6XBQAAYEDEdAMAAEB4EUQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAxBBAAAGDPDdAPqKZVK2rt3rzo6OtTW1ma6OQAAoAHlclkHDhzQ/PnzFYnUH/PwdRDZu3ev4vG46WYAAIAmDA8PKxaL1f0ZXweRjo4OSYf+IZ2dnYZbAwAAGlEoFBSPxyv9eD2+DiJj0zGdnZ0EEQAAAqaRsgqKVQEAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAgICyLEvZbFaWZZluStMIIgAABFAmk1EikVBPT48SiYQymYzpJjWlrVwul003opZCoaCuri6NjIxw6B0AAIdZlqVEIqFSqVS5Fo1Glc/nFYvFDLbsEDv9NyMiAAAETC6XGxdCJKlYLGpoaMhQi5pHEAEAIGCSyaQikfFdeDQaVXd3t6EWNY8gAgBAwMRiMQ0ODioajUo6FEK2bt3qi2kZu6gRAQAgoCzL0tDQkLq7u30VQuz03zM8ahMAAHBYLBbzVQBpBlMzAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAYJBlWcpms7Isy3RTjGg6iDz77LO65JJLNH/+fLW1tek3v/nNuL8vl8u68cYbddxxx2n27NlatmyZcrncdNsLAEDLyGQySiQS6unpUSKRUCaTMd0kzzUdRA4ePKgzzjhDd911V9W/v+2223TnnXfqnnvu0Z///GcdffTR6u3t1f/+97+mGwsAQNDUGvGwLEt9fX2VU3RLpZL6+/tDNzLSdBC58MILtX79en3lK1+Z9HflclkbN27UD3/4Q335y1/W6aefrp/97Gfau3fvpJETAABaVb0Rj1wuVwkhY4rFooaGhiSFZ8rGlRqRt99+W/v27dOyZcsq17q6urR48WLt3LnTjVsCAOArU414JJNJRSLju+FoNKru7u5QTdm4EkT27dsnSZo7d+6463Pnzq38XTWjo6MqFArjHgAABNFUIx6xWEyDg4OKRqOSDoWQrVu3SlKopmx8tWpmYGBAXV1dlUc8HjfdJAAAmlJvxGNMOp1WPp9XNptVPp9XOp2eMsC0GleCyLx58yRJ+/fvH3d9//79lb+rZt26dRoZGak8hoeH3WgeAACuqzXiEYvFJv1cKpWqXG8kwLQSV4LIggULNG/ePG3fvr1yrVAo6M9//rOWLFlS83nt7e3q7Owc9wAAIKiqjXhMpdEA0ypmNPvEDz74YNww0dtvv62XX35ZxxxzjI4//nitXr1a69evVzKZ1IIFC3TDDTdo/vz5uvTSS51oNwAAgRCLxWyHiHQ6rd7eXg0NDam7u7tlQ4g0jSDy4osv6vzzz6/8ec2aNZKk5cuX6/7779d3v/tdHTx4UH19ffr3v/+tz33uc3riiSc0a9as6bcaAIAW10yACaK2crlcNt2IWgqFgrq6ujQyMsI0DQAADrMsS7lcTslk0tHQY6f/9tWqGQAAMDUnNjvzy14lBBEAAALEiQDhp+3lCSIAAASEUwHCT3uVEEQAAAgIpwKEn/YqIYgAAOAQtw+qcypA+GmvEoIIAAAO8KL408kA0cxma25g+S4AANNkWZYSicS4aZNoNKp8Pu/KKINlWb7e7MxO/930hmYAAOCQerUbbgSFVtrsjKkZAACmyU/Fn0FDEAEAYJr8VPwZNNSIAADgEL/XbniFGhEAAAxopdoNrzA1AwAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAANAC3D5wzy0EEQAAAs6LA/fcwoZmAAAEmNcH7jXCTv/NiAgAAAFW78C9ICCIAAAQYEE/cI8gAgBAgAX9wD1qRAAAaAF+OnCPQ+8AADDMsizlcjklk0lPgkFQD9xjagYAAIcFeTmt15iaAQDAQX5cTus1lu8CAGBI0JfTeo0gAgBAFdW2TG9kG/WgL6f1GkEEAIAJqtV4NFr3EfTltF6jRgQAgCNUq/EYG+GwU/fhp+W0XmP5LgAATapW4zHxz9LHdR+1QkZQl9N6jakZAEBFUI+Sd1K1Go9IJELdh0sIIgAASex9MaZajcfg4CB1Hy6hRgQAwN4XVVSr8Qhz3Ycd1IgAAGypt/eF3Q7X663N3VKtxoO6D+cxNQMAcGzvC6Z3YBdBBADgyN4XlmWpr6+vMrJSKpXU398f6sJXTI2pGQCAJCmdTqu3t7fpGggnp3ca1SrTQGHm2ohIsVjUDTfcoAULFmj27Nk66aST9KMf/Ug+ro0FgNCLxWJKpVJNdepeb23ONFBrcC2IbNiwQVu2bNHmzZv12muvacOGDbrtttu0adMmt24JADDIy63NmQZqHa5NzTz33HP68pe/rIsuukiSdMIJJ+jnP/+5XnjhBbduCQAwbLrTO40yMQ0Ed7g2InLeeedp+/btevPNNyVJf/3rX7Vjxw5deOGFNZ8zOjqqQqEw7gEACJbpTO80qto0kCTt2rXLtXvCHa4FkbVr1+ryyy/XySefrKOOOkpnnXWWVq9erSuuuKLmcwYGBtTV1VV5xONxt5oHAAiwWCymDRs2TLq+bt06pmcCxrUg8vDDD+uhhx7Stm3b9NJLL+mBBx7QHXfcoQceeKDmc9atW6eRkZHKY3h42K3mAQA84tb5Neecc86ka2PTMwgO12pEvvOd71RGRSTptNNO0zvvvKOBgQEtX7686nPa29vV3t7uVpMAAB7LZDKVotJIJKLBwUGl02lHXntsembitvRhP4guaEuaXRsR+c9//lN1GVe1o5QBIMg4sbY6t1e2eLlKJyiCuKTZtSByySWX6Oabb9bvf/975fN5PfLII/rJT36ir3zlK27dEgA8F8Rf/F6pt7KlHjvBLp1OK5/PK5vNKp/POzbaEkRBXdLs2um7Bw4c0A033KBHHnlE7733nubPn69vfOMbuvHGGzVz5syGXoPTdwH4GSfW1tfM++PmVE6ry2az6unpqXo9lUp52hY7/bdrIyIdHR3auHGj3nnnHf33v//V3//+d61fv77hEAIAftfsN/6wsDt1EtRv9H7h9c62TuHQOwBoUlB/8XvJztQJwW56glozw6F3ANCksV/8/f39KhaLgfnF77VYLNbQe9LMKpigrRBxm1c72zqJEREAmAaKJZ0z8Rt9JBLRrbfeWreehELhybzY2dZJrhWrOoFiVQAIn9tvv13f+973VC6XaxasUijsb74oVgUAwC7LsrR27VqNfUeuVbBKPUnrIIgAAHyj0YBBoXDrIIgAAHyj0YAR1BUimIwgAgAwbmw3VUkNB4ywFAq3+hECBBEAaFCrdwimTFz9IqnhgBG0FSJ2hWFlEKtmAKABbD3uDla/1Bbk94ZVMwDgILYet6/R0SNWv9QWlveGIAIAUwhLh+AUO9MJrH6pLSzvDUEEAKYQlg7BCXZHj1j9UltY3huCCABMISwdghOaGT0Ky+qXZoThvaFYFQAaZFlWoA4TMyHIBZZwDsWqAOCCVl8q6gRGj2AXIyIAAMcxehRudvrvGR61CQAQIrFYzLEAYlmWcrmckskkoaYFMTUDAPCtMOwsGnZMzQBAi3Bq5MAvIxAUvgYXxaoAEDJOjRz4aQSCjeTCgRERAAg4OyMH9UY7/DYC4bf2oHGMiABAiDQ6cjDVaIepEYha59KwFDgcGBEBgICZOKrRyMiBUz/jtEZONWYpcPAwIgIALaraqEYjIweNjHZ4PQLR6Lk01TaSa/R0X/gfIyIAEBBTjVjUGzmwW0fixQhENptVT09P1eupVKrm8xoZRYFZjIgAQAuaalSj3hb0dkY7vNrKvplTje2e7gv/I4gAQEA003EfyW8nuTYzFcSS3tbDFu8AcAS/bOZVzVjH3d/fr2Kx2FQNh5NbrzshnU6rt7e34amgsTA2cYqp0TAG/2FEBAAO82Izr+kWWfptVMMJdqaC6o2iUMAaTBSrAoC8WbpKkaVzJhbU8t76i53+myACAGp+BUej2CXUPby3/sOqGQCwabqFoFOhyNI9vLfBRhABALm/mZfbQSfMeG+DjSACAIe5WQjabNChAHNqnEkTbNSIAICH7OxaSgGmPZxJ4x8UqwJAwFGAiSCjWBUAAo4CTIQFQQQAfIgCTISFq0Fkz549uvLKK3Xsscdq9uzZOu200/Tiiy+6eUsAaAkUYCIsXDtr5l//+peWLl2q888/X48//rg+9alPKZfL6ROf+IRbtwSAlmL3HBYgiFwLIhs2bFA8Htd9991XubZgwQK3bgcALclvh9QBTnNtaubRRx/VwoULddlll2nOnDk666yzdO+999Z9zujoqAqFwrgHAABoXa4FkbfeektbtmxRMpnUk08+qWuuuUbXX3+9HnjggZrPGRgYUFdXV+URj8fdah4AAPAB1/YRmTlzphYuXKjnnnuucu3666/Xrl27tHPnzqrPGR0d1ejoaOXPhUJB8XicfUQAAAgQX+wjctxxx+kzn/nMuGuf/vSn9e6779Z8Tnt7uzo7O8c9AABA63ItiCxdulRvvPHGuGtvvvmmEomEW7cEADSB82xgkmtB5Nvf/raef/553XLLLRoaGtK2bds0ODiolStXunVLAIBNmUxGiURCPT09SiQSymQyppuEkHH1rJnHHntM69atUy6X04IFC7RmzRqtWLGi4edz1gwAuIfzbOAWO/23a/uISNLFF1+siy++2M1bAIAnLMtSLpdTMpl0rJN24zXtqHeeDUEEXuGsGQCYghvTF36YEuE8G/gBQQQAarAsSw8//LBWrFhRGTkolUrq7++fVmGnZVnq6+tz9DWbwXk28ANXp2YAIKgymcy4sHCk6U5f+GlKhPNsYBpBBAAmmDhiMVEz0xdH1oOMTYlMLBI1NSXCeTYwiakZAJig2ojFmGamLybWgzz55JNMiQCHubp8d7pYvgtMn+mVGUFUa1nrz3/+cy1ZssTW+1hviawkpkTQknyxxTsA8/ywMiOIahVxXnbZZbYDw1T1IKlUihCCUGNEBGhRbFY1fZZlTXvEgs8BYcSICIC638TRGCdGLFgiC9THqhmgRfltZUaYsUQWqI0REaBF8U3cX6gHAaqjRgRocU7UOQCAHb459A6AeWxWFVwsvUYYMDUDYEqWZSmbzXp+FkqYsfQaYUEQAVAXHaL3/HIoHuAFggiAmugQzWDpNcKEIAKgJj90iGGcFhpben0kll6jVRFEAIe1UsdpukMM67SQH5Zet9J/x/A3ggjgoFbrOE12iGGfFurt7dWmTZt04403aufOnUqn057du9X+O4a/sY8I4BC3zxQxuZTTxF4k2WxWPT09Va+nUilP2mBKJpPRihUrNPbrua2tTffee68nYYSzceAEzpoBJvBimNnNegrT31BN7ApqelrIlLGRoCO/I5bLZfX19XkyGuSHuiCEC0EELc+rTtytjtMvUxRe1wz4oU6iEU6/L9WCgHToc/ciDIQ1AMIcgghampeduFsdpx++oZoakUmn08rn88pms8rn857WSTTCjfelWhCQpEgk4kkYCEoAROugRgQtzUSdgdP1FF7P2U+sRaFmoDo335dMJjMuQHtZIzKGM4owHZw1Axw29u1yYmfh5jdLp892GfuG2t/fr2Kx6Oo31CM7wEgkosHBQZ144ok1R2TC3EHVG6ma7vuSTqfV29urnTt3SpKWLFni+XvNGUXwCiMiaHmZTGZSJ+63If5GuP0NtdY3/J07d+rcc89lRGQCRoqA2lg1AxzB73UGjXJ75Uqtb/gHDx6kZqAKaikAZzAiAkDS1N/wqRmojvcFmIwaEQC2TVWLQs1AdbwvwPQwIgJgHL7hA5guRkQANI1v+AC8RLEqAAAwhiACAACMIYgAAABjCCIAMIHXB/wBYUYQARBIdsNCoz9v6oA/IKwIIgACx25YaPTnvTytGcAhBBEAgWI3LNj5+XoH2QFwB0EEQKDYDQt2fn7stOYjuX1aMxB2ngWRW2+9VW1tbVq9erVXtwTQguyGBTs/z0F2gPc8CSK7du3S1q1bdfrpp3txO6DlhXlVh92wYPfnW+W0ZiAoXD9r5oMPPtDZZ5+tu+++W+vXr9eZZ56pjRs3NvRczpoBJstkMpWah0gkosHBwVB2lnbPxOEMHcA7dvpv14PI8uXLdcwxx+inP/2pUqlU3SAyOjqq0dHRyp8LhYLi8ThBBL5jWZZyuZySyaSnnZplWUokEuNqHqLRqPL5PJ0rAN+wE0RcnZr5xS9+oZdeekkDAwMN/fzAwIC6uroqj3g87mbzgKaY3GfCzVUdYZ7uAWCOa0FkeHhYq1at0kMPPaRZs2Y19Jx169ZpZGSk8hgeHnareUBTTO8z4daqDjbxAmCKa0Fk9+7deu+993T22WdrxowZmjFjhp555hndeeedmjFjhorF4qTntLe3q7Ozc9wD8BPT+0y4sarDdLgCEG4z3HrhL37xi3rllVfGXbv66qt18skn63vf+17lFykQJGMjEhNrNLzcZyKdTqu3t9exwst64crpuhNTtTUA/Mu1INLR0aFTTz113LWjjz5axx577KTrQFCMjUj09/erWCwa22ciFos5dk+vwhWrfQBUw86qgE2tts+EF5t4Mf0DoBbXRkSqefrpp728HeAaJ0ck/MDp6Z6JnJj+YVoHaE2MiACQdChcpVIpVzr56a72aWRVD8uPgWAiiABw3XSmfxqZ1mH5MRBcru+sOh1s8Q60lma2Wc9ms+rp6al6PZVKsdss4EN2+m9Pa0TgLebU4TfN1NZMtarHy+XHAJzH1EyLYqgarWKqaR23dpsF4A2mZloQQ9VoRfWmdTKZzKS9XYK+rBoIMqZmQo6haviNE9OE9aZ13F5+DMA9TM20IIaq4SdeTRO6ufwYgHsIIi3Ii50yncT+D62LHVUBTIUg0qKCsg05RbWtzfRpxQD8j2JVGENRbevjMwbCyU7/zYgIjOHbsj1BnMIK2jQhAO8RRGAMRbWNC/IUVlCmCQGYwdQMjGL/h6kxvQEgaNhHBIHB/g9TY18YAK2MIALjmjl/JEymOmsFAIKMGhHA5yj4BNDKqBEBAqLeWSsA4CfUiAAecOL8FDuYwgLQipiaAZoQ5OW0AOAnTM0ADThy9EMSy2kBoA52VgUcNHH04//+7//YERYAHMKICFBHtc3ExnaDZUQEAKpjRARwSLXNxEqlktasWcNy2jqCeC4OADMIIkAdtc7DWbVq1aTzU+h8D6GQF4AdBBGgjnqbicViMaVSKcViMTrfwyzLUl9fX2UUqVQqqb+/P/ThDEBt1IgADai3mRiH0n0sm82qp6en6vVUKuV9gwAYwYZmgMPqbSbGoXQf41wcAHYxNQNMU606kjB2vpyLA8AugggwTXS+46XT6UmFvABQCzUigEM4lA4ADqFGBDCAQ+kAwD6mZgAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAY42oQGRgY0KJFi9TR0aE5c+bo0ksv1RtvvOHmLQEAQIC4GkSeeeYZrVy5Us8//7yeeuopffTRR7rgggt08OBBN28bOJZlKZvNclR6i+LzBYDaXA0iTzzxhK666iqdcsopOuOMM3T//ffr3Xff1e7du928baBkMhklEgn19PQokUgok8mYbhIcxOcLAPV5etbM0NCQksmkXnnlFZ166qmT/n50dFSjo6OVPxcKBcXj8ZY9a8ayLCUSiUlHpufzebYKb5BlWcrlckomk757z/h8AYSVnbNmPCtWLZVKWr16tZYuXVo1hEiHakq6uroqj3g87lXzjMjlcuM6KUkqFosaGhoy1KJg8ftoA58vAEzNsxGRa665Ro8//rh27NhR89sgIyJ8Y26U1+9dMyMvfL4Awsp3IyLXXnutHnvsMWWz2bq/gNvb29XZ2Tnu4RduFBzGYjENDg4qGo1KOtRJbd26lU6qAV6ONjQ78sLnCwBTc3VEpFwu67rrrtMjjzyip59+Wslk0tbz7SQqN2UyGfX19alUKikSiWhwcFDpdNqx17csS0NDQ+ru7qaTapBXow1O3IfPF0DY+GZEZOXKlXrwwQe1bds2dXR0aN++fdq3b5/++9//unlbR1mWVQkh0qFal/7+fsdHRlKpFJ2UDV6NNjgx8sLnCwC1zXDzxbds2SJJSqVS467fd999uuqqq9y8tWPqdUR0LGal02n19va6OtqQTCYViUQmjYh0d3c7fi8ACCNXg4iHK4NdQ0fkb7FYzNVAODby0t/fr2KxSJ0HADiMs2amQMEh0um08vm8stms8vm8o/VBABB2nm5oZpdfilUlCg79xM+bmAEAfFSs2kpaveAwKOeh+H0TMwCAPQQRD/m1sw9K5+7FCiYAgLcIIh7xa2cfpM6dLdMBoPUQRDzg587eTuduekRnbAXTkVjBBADBRhDxgJ+/yTfaufthRIcVTADQelg14wG/H36WyWQm7ZNx5BJVv7WfFUwA4G+smvEZv3+Tn2qfDL+N6LT6CiYACBNGRDwU1G/yJkdE2DMEAIKHERGfCtI3+SMLU02N6PihLgUA4C5GRHzCT9/8M5lMZZVPJBLR4OCg0um0pyM6fqtLAQA0jhGRgPHTN/96S429HNHxW10KAMAdBBHD/LbHiF8CAHuGAEA4EEQM80vHP8YvAcDvK40AAM4giBjml45/jJ8CwFTLigEAwUexqg9MtaGYCUFdagwAMM9O/00Q8Qk6fgBAq7DTf8/wqE2YQiwWI4AAAEKHGhFIqn2yrukTdwEArY0ggpr7mPhpfxMAQGuiRiTkau1gunPnTp177rnsbAoAsI2dVdGwWvuY7Nixw1f7mwAAWhNBxAA/1V3U2sfkc5/7nK/2NwEAtCaCiMf8VndRawOzRYsW+WZjMwBA66JGxEN+PlG21j4m7G8CALCLfUR8qt65MqY7+Vr7mLC/CQDATUzNeMhv58oAAGAaQcRDfjpQDgAAP6BGxADqLgAArYwaERdZlqVcLqdkMtl0iJhYd+HEa6I63lsA8DemZmxwY+mt35bzthLeWwDwP6ZmGuTG0ls/L+cNOt5bADCHLd5dUG/prcnX9NMurX7ixucFAHAeQaRBbiy9ne5rMvVQG0ulASAYCCINcmPp7XRe07Is9fX1Vb71l0ol9ff3MzJyGEulASAYqBGxyY2lt828ZjabVU9PT9XrqVTKkXa1ApZKA4D37PTfBJGAohgTAOBXvipWveuuu3TCCSdo1qxZWrx4sV544QW3bxkKTD0AAFqBqyMiv/zlL/XNb35T99xzjxYvXqyNGzfqV7/6ld544w3NmTNnyuczIjI1ph4AAH7jm6mZxYsXa9GiRdq8ebOkQwWV8Xhc1113ndauXTvl8wkiAAAEjy+mZj788EPt3r1by5Yt+/hmkYiWLVumnTt3Vn3O6OioCoXCuAemh31GAAB+5loQ+cc//qFisai5c+eOuz537lzt27ev6nMGBgbU1dVVecTjcbeaFwrsMwIA8Dtf7SOybt06jYyMVB7Dw8OmmxRY7DMCAAgC107f/eQnP6loNKr9+/ePu75//37Nmzev6nPa29vV3t7uVpNCpd4W5xS1AgD8wrURkZkzZ+qcc87R9u3bK9dKpZK2b9+uJUuWuHVbHMYW5wCAIHB1ambNmjW699579cADD+i1117TNddco4MHD+rqq69287aOC2LBJ/uMAACCwLWpGUn6+te/rvfff1833nij9u3bpzPPPFNPPPHEpAJWP8tkMpVai0gkosHBQaXTadPNakg6nVZvby/7jAAAfIst3utgG3UAAOzzxT4iraBewWctQZzGAQDAFIJIHXYLPtm3AwAAewgiddgp+GTfDgAA7HO1WLUVNFrw6Zd9OyzLUi6XUzKZpI4FAOB7BJEaJnboU3XqY9M4Ewtbvdy3I8grfAAA4cTUTBXN1HqY3rfDL1NDFOsCAOwgiEwwnQ49nU4rn88rm80qn897OhrRzAofp1GsCwCwiyAywXQ79FgsplQq5Xl9hukt3f0yIgMACJbQBpFaUwimO/RmmZ4a8sOIDAAgeEIZROpNIZju0KfD5NRQUAMcAMCs0G3x3ui27ZZlcUaLTZlMRv39/SoWi5UAx6odAAgfO/136JbvNrrfRyNLdjEeh+wBAOwKXRBpZr8PNglrHAEOAGBH6GpE7NaAsCQVAAD3hK5GZEwjNSCN1pMAAICPUSPSgEamEPxyfgwAAK0qdFMzdrAkFQAAdxFE6gjyniIAAARBaGtE7GBPEQAAGkeNiMNYkgoAgDuYmgEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgjCtBJJ/PK51Oa8GCBZo9e7ZOOukk3XTTTfrwww/duB0AAAioGW686Ouvv65SqaStW7equ7tbr776qlasWKGDBw/qjjvucOOWtlmWpVwup2QyqVgsZro5AACEUlu5XC57caPbb79dW7Zs0VtvvdXwcwqFgrq6ujQyMqLOzk7H2pLJZNTX16dSqaRIJKLBwUGl02nHXh8AgDCz0397ViMyMjKiY445pu7PjI6OqlAojHs4zbKsSgiRpFKppP7+flmW5fi9AABAfZ4EkaGhIW3atEn9/f11f25gYEBdXV2VRzwed7wtuVyuEkLGFItFDQ0NOX4vAABQn60gsnbtWrW1tdV9vP766+Oes2fPHn3pS1/SZZddphUrVtR9/XXr1mlkZKTyGB4etv8vmkIymVQkMv6fHY1G1d3d7fi9AABAfbZqRN5//33985//rPszJ554ombOnClJ2rt3r1KplM4991zdf//9kwLAVNysEenv71exWFQ0GtXWrVupEQEAwCF2+m/XilX37Nmj888/X+ecc44efPBBRaNR26/hVhCRDtWKDA0Nqbu7m1UzAAA4yE7/7cry3T179iiVSimRSOiOO+7Q+++/X/m7efPmuXFL22KxGAEEAADDXAkiTz31lIaGhjQ0NDSps/dotTAAAAgAV1bNXHXVVSqXy1UfAAAAYzhrBgAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRCawLEvZbFaWZZluCgAALY8gcoRMJqNEIqGenh4lEgllMhnTTQIAoKW1lX18El2hUFBXV5dGRkbU2dnp6r0sy1IikVCpVKpci0ajyufzk04QBgAAtdnpvxkROSyXy40LIZJULBY1NDRkqEUAALQ+gshhyWRSkcj4tyMajaq7u9tQiwAAaH0EkcNisZgGBwcVjUYlHQohW7duZVoGAAAXUSMygWVZGhoaUnd3NyEEAIAm2Om/Z3jUpsCIxWIEEAAAPMLUDAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGN8fdbM2Hl8hULBcEsAAECjxvrtRs7V9XUQOXDggCQpHo8bbgkAALDrwIED6urqqvszbeVG4oohpVJJe/fuVUdHh9ra2hx5zUKhoHg8ruHh4SmPJoZ3+Fz8ic/Fn/hc/InP5WPlclkHDhzQ/PnzFYnUrwLx9YhIJBJRLBZz5bU7OztD/x+KH/G5+BOfiz/xufgTn8shU42EjKFYFQAAGEMQAQAAxoQuiLS3t+umm25Se3u76abgCHwu/sTn4k98Lv7E59IcXxerAgCA1ha6EREAAOAfBBEAAGAMQQQAABhDEAEAAMaEKojcddddOuGEEzRr1iwtXrxYL7zwgukmhdrAwIAWLVqkjo4OzZkzR5deeqneeOMN083CBLfeeqva2tq0evVq000JvT179ujKK6/Uscceq9mzZ+u0007Tiy++aLpZoVYsFnXDDTdowYIFmj17tk466ST96Ec/auiMFRwSmiDyy1/+UmvWrNFNN92kl156SWeccYZ6e3v13nvvmW5aaD3zzDNauXKlnn/+eT311FP66KOPdMEFF+jgwYOmm4bDdu3apa1bt+r000833ZTQ+9e//qWlS5fqqKOO0uOPP66//e1v+vGPf6xPfOITppsWahs2bNCWLVu0efNmvfbaa9qwYYNuu+02bdq0yXTTAiM0y3cXL16sRYsWafPmzZIOnWMTj8d13XXXae3atYZbB0l6//33NWfOHD3zzDP6/Oc/b7o5offBBx/o7LPP1t13363169frzDPP1MaNG003K7TWrl2rP/3pT/rjH/9ouik4wsUXX6y5c+cqk8lUrn31q1/V7Nmz9eCDDxpsWXCEYkTkww8/1O7du7Vs2bLKtUgkomXLlmnnzp0GW4YjjYyMSJKOOeYYwy2BJK1cuVIXXXTRuP9vYM6jjz6qhQsX6rLLLtOcOXN01lln6d577zXdrNA777zztH37dr355puSpL/+9a/asWOHLrzwQsMtCw5fH3rnlH/84x8qFouaO3fuuOtz587V66+/bqhVOFKpVNLq1au1dOlSnXrqqaabE3q/+MUv9NJLL2nXrl2mm4LD3nrrLW3ZskVr1qzR97//fe3atUvXX3+9Zs6cqeXLl5tuXmitXbtWhUJBJ598sqLRqIrFom6++WZdccUVppsWGKEIIvC/lStX6tVXX9WOHTtMNyX0hoeHtWrVKj311FOaNWuW6ebgsFKppIULF+qWW26RJJ111ll69dVXdc899xBEDHr44Yf10EMPadu2bTrllFP08ssva/Xq1Zo/fz6fS4NCEUQ++clPKhqNav/+/eOu79+/X/PmzTPUKoy59tpr9dhjj+nZZ59VLBYz3ZzQ2717t9577z2dffbZlWvFYlHPPvusNm/erNHRUUWjUYMtDKfjjjtOn/nMZ8Zd+/SnP61f//rXhloESfrOd76jtWvX6vLLL5cknXbaaXrnnXc0MDBAEGlQKGpEZs6cqXPOOUfbt2+vXCuVStq+fbuWLFlisGXhVi6Xde211+qRRx7RH/7wBy1YsMB0kyDpi1/8ol555RW9/PLLlcfChQt1xRVX6OWXXyaEGLJ06dJJy9vffPNNJRIJQy2CJP3nP/9RJDK+K41GoyqVSoZaFDyhGBGRpDVr1mj58uVauHChPvvZz2rjxo06ePCgrr76atNNC62VK1dq27Zt+u1vf6uOjg7t27dPktTV1aXZs2cbbl14dXR0TKrTOfroo3XsscdSv2PQt7/9bZ133nm65ZZb9LWvfU0vvPCCBgcHNTg4aLppoXbJJZfo5ptv1vHHH69TTjlFf/nLX/STn/xE3/rWt0w3LTjKIbJp06by8ccfX545c2b5s5/9bPn555833aRQk1T1cd9995luGib4whe+UF61apXpZoTe7373u/Kpp55abm9vL5988snlwcFB000KvUKhUF61alX5+OOPL8+aNat84oknln/wgx+UR0dHTTctMEKzjwgAAPCfUNSIAAAAfyKIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMOb/AWncnKDfPxO4AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(x,y,'.k')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [], + "source": [ + "def MyOLS(x,y):\n", + " # for implementation for computing theta1:\n", + " xAve = x.mean()\n", + " yAve = y.mean()\n", + " num = 0\n", + " den = 0\n", + " for i in range(len(x)):\n", + " num = num + x[i]*(y[i]-yAve)\n", + " den = den + x[i]*(x[i]-xAve)\n", + " theta1Hat = num/den\n", + " theta0Hat = yAve - theta1Hat*xAve\n", + " return theta0Hat, theta1Hat" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([1.12539439])" + ] + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "the0, the1 = MyOLS(x,y)\n", + "the1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# TODO - Students\n", + "- [ ] Efficacy --> time: For method Vs. Numpy" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "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.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}