ols-video/main.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ordinary Least Squares\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The Ordinary Least Squares (OLS) is an important method in machine learning and statistics for several reasons. The OLS is a straightforward and easy-to-understand method for fitting linear models. \n",
    "\n",
    "The model minimizes the sum of squared differences between the observed $y$ and predicted values $\\hat{y}$, making it intuitive to understand. \n",
    "\n",
    "Additionally, the OLS is the foundation for linear regression, one of the most widely used machine learning and statistics techniques. Linear regression is valuable for modeling relationships between variables when you suspect a linear relationship exists.\n",
    "\n",
    "While OLS is valuable in many scenarios, it's essential to acknowledge its limitations, especially when dealing with nonlinear relationships or complex data structures. \n",
    "\n",
    "More advanced machine learning techniques like decision trees, neural networks, or support vector machines may be more appropriate in such cases. However, OLS remains a foundational method with enduring relevance in machine learning and statistics.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Obtaining the OLS general equation\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(\\theta_0, \\theta_1)}{\\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": 3,
   "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": 3,
     "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": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAgn0lEQVR4nO3dfbBcZWE/8OfeG5JQCghUgUzCEt31DRBfEEbtjMlCzdBMhD+sylDNKJ3ETipQHK2ZqS8MthHHUbRkYFdXaesL2tpExzekuJgqUElCWmor7rW5mIiY6YwlBGrq7O5vnvP73fvLhrzcm+x9ds+5n8/M8eGcey73eTh37/n6PM85z0i32+0GAIBERlP9IACASPgAAJISPgCApIQPACAp4QMASEr4AACSEj4AgKSEDwAgqXlhyHQ6nfDYY4+Fk08+OYyMjAy6OgDANMR3lj755JNh0aJFYXR0NF/hIwaPJUuWDLoaAMAx2LVrV1i8eHG+wkfs8Zis/CmnnDLo6gAA07B3796s82DyPp6r8DE51BKDh/ABAPkynSkTJpwCAEkJHwBAUsIHAJCU8AEAJCV8AABJCR8AQFLCBwCQlPABACQlfAAASQkfAEBSwgcAzCG7d+8OzWYzKwdF+ACAOaLRaIRSqRSq1WpWxv1BGOl2u90wZKvinXrqqeGJJ56wsBwA9Ens6YiBo9PpTB0bGxsLExMTYfHixUnv33o+AGAOaLVaPcEjarfbYXx8PHldZhw+tmzZElatWhUWLVqULZu7efPmw577jne8IzvnlltuOd56AgDHoVKphNHR3tt+7Pkol8th6MPHU089FS688MKwcePGI563adOm8MADD2QhBQAYrMWLF4d6vZ4FjiiWtVqtL0MuMzVvpt9w+eWXZ9uR/PznPw/vfOc7w1133RVWrlx5PPUDAPrkmmuuCStWrMiGWmKPxyCCxzGFj6OJ40lvectbwrvf/e5w3nnnHfX8/fv3Z9uBE1YAgNkRA8egQsesTTi9+eabw7x588K11147rfM3bNiQzY6d3JYsWdLvKgEAQ6Sv4WPbtm3hE5/4RLjjjjuyiabTsX79+uyxnMlt165d/awSAFDk8PFP//RPYc+ePeGcc87Jej/i9uijj4Z3vetd4dxzzz3k9yxYsCB7HvjADQAorr7O+YhzPS677LKeY3FiSzz+tre9rZ8/CgCYK+Fj3759PS8k2blzZ9ixY0c4/fTTsx6PM844o+f8E044IZx11lnhBS94QX9qDADMrfCxdevWsHz58qn9G264IStXr16dzfUAAOhr+Fi2bFmYyXIw8Z3xAACTrO0CACQlfAAASQkfAEBSwgcADJHdu3eHZrOZlUUlfADAkGg0GqFUKoVqtZqVcb+IRrozeXQlgbiwXFzjJb5q3dtOAZgrdu/enQWOuEDrpLjsfXxqdNALwfX7/q3nAwCGQKvV6gkeUbvd7nmxZ1EIHwAwBCqVShgd7b0tx56Pcrkcikb4AIAhsHjx4lCv17PAEcWyVqvlYshlpsz5AIAhm/sxPj6e9XjkKXjM5P7d11VtAYDjs3jx4lyFjmNh2AUASEr4AACSEj4AgKSEDwAgKeEDAEhK+AAAkhI+AICkhA8AICnhAwBISvgAAJISPgBghmuvNJvNrOTYCB8AME2NRiOUSqVQrVazMu4zc1a1BYBpiD0dMXB0Op2pY3HZ+4mJicIvBNfv+7eeDwCYhlar1RM8ona7HcbHxwdWp7wSPgBgGiqVShgd7b1txp6Pcrk8sDrllfABANMQh1bq9XoWOKJY1mo1Qy7HwJwPAJjh3I841BJ7PASPY7t/zzviVwGAHjFwCB3Hx7ALAJCU8AEAJCV8AABJCR8AQFLCBwCQlPABACQlfAAASQkfAEBSwgcAkJTwAQAkJXwAUMj1V5rNZlYyfIQPAAql0WiEUqkUqtVqVsZ9hotVbQEojNjTEQNHp9OZOjY2NhYmJiYsBjdE9289HwAURqvV6gkeUbvdDuPj4wOrE88kfABQGJVKJYyO9t7aYs9HuVweWJ14JuEDgMKIQyv1ej0LHFEsa7WaIZchY84HAIWc+xGHWmKPh+AxfPfveYnqBADJxMAhdAwvwy4AQFLCBwCQlPABACQlfAAASQkfAEBSwgcAkJTwAQAkJXwAAMMdPrZs2RJWrVoVFi1aFEZGRsLmzZunvvab3/wm/Nmf/Vm44IILwkknnZSd89a3vjU89thj/a43ADBXwsdTTz0VLrzwwrBx48ZnfO3pp58O27dvD+973/uy8h/+4R/CI488El7/+tf3q74AFOT1581mMyuZe45rbZfY87Fp06Zw5ZVXHvacBx98MFx88cXh0UcfDeecc85R/53WdgEotkajEdasWRM6nU62Am1cCO6aa64ZdLU4TjO5f8/6nI9YiRhSnvWsZx3y6/v3788qfOAGQDHFno7J4BHFcu3atXpA5phZDR+//vWvszkgV1111WFT0IYNG7KkNLktWbJkNqsEwAC1Wq2p4DGp3W5nK9Ayd8xa+IiTT9/4xjeGOKpz2223Hfa89evXZ70jk9uuXbtmq0oADFilUsmGWg40NjYWyuXywOpEQcLHZPCI8zzuvvvuI479LFiwIPv6gRsAxRSXuY9zPGLgiGJZq9Wy48wd82YreMSutTiT+Ywzzuj3jwAgx+Lk0hUrVmRDLbHHQ/CYe2YcPvbt29czNrdz586wY8eOcPrpp4ezzz47vOENb8ges/3617+ejeM9/vjj2Xnx6/Pnz+9v7QHIpRg4hI65a8aP2t57771h+fLlzzi+evXq8MEPfjAsXbr0kN8Xe0GWLVt21H+/R20BIH9mcv+ecc9HDBBHyivH8doQAGAOsLYLAJCU8AEAJCV8AABJCR8AQFLCBwCQlPABACQlfAAASQkfAEBSwgcAz7B79+7szdSxhH4TPgDo0Wg0QqlUCtVqNSvjPgx0bZfZZm0XgMGJPR0xcHQ6naljcdn7iYkJC8HRt/u3ng8AprRarZ7gEcUVyg9czRyOl/ABwJRKpRJGR3tvDbHno1wuD6xOFI/wAcCUOLRSr9ezwBHFslarGXKhr8z5AOCQcz/iUEvs8RA86Pf9e960/o0AzCkxcAgdzBbDLgBAUsIHAJCU8AEAJCV8AABJCR8AQFLCBwCQlPABACQlfAAASQkfAEBSwgcAkJTwAZDTtVeazWZWQt4IHwA502g0QqlUCtVqNSvjPuSJVW0BciT2dMTA0el0po7FZe8nJiYsBEdu7t96PgBypNVq9QSPqN1uh/Hx8YHVCWZK+ADIkUqlEkZHe/90x56Pcrk8sDrBTAkfADkSh1bq9XoWOKJY1mo1Qy7kijkfADmd+xGHWmKPh+BB3u7f85LVCoC+iYFD6CCvDLsAAEkJHwBAUsIHAJCU8AEAJCV8AABJCR8AQFLCBwCQlPABACQlfAAASQkfALP4CvRms5mVwP8nfADMgkajEUqlUqhWq1kZ94H/y8JyAH0Wezpi4Oh0OlPH4uqzExMT1mOhsGZy/9bzAdBnrVarJ3hE7XY7W4UWED4A+q5SqYTR0d4/r7Hno1wuD6xOMEyED4A+i0Mr9Xo9CxxRLGu1miEX+H/M+QCYxbkfcagl9ngIHhTd3hncv+clqxXAHBMDh9ABz2TYBQBISvgAAJISPgCApIQPAGC4w8eWLVvCqlWrwqJFi8LIyEjYvHlzz9fjwzPvf//7w9lnnx1OPPHEcNlll2Uv3AEAOKbw8dRTT4ULL7wwbNy48ZBf/8hHPhI++clPhttvvz388z//czjppJPCihUrwq9//Wv/xQGAmT9qe/nll2fbocRej1tuuSX8+Z//ebjiiiuyY3/zN38TzjzzzKyH5M1vfvPx1xgAyLW+zvnYuXNnePzxx7OhlknxhSOXXHJJuP/++w/5Pfv3789eTHLgBgAUV1/DRwweUezpOFDcn/zawTZs2JAFlMltyZIl/awSADBkBv60y/r167NXsU5uu3btGnSVAIC8hI+zzjorK3/5y1/2HI/7k1872IIFC7J3wB+4AaRae6XZbGYlkNPwsXTp0ixk3HPPPVPH4hyO+NTLq171qn7+KIDj0mg0QqlUCtVqNSvjPjCk4WPfvn1hx44d2TY5yTT+889+9rPsvR/XX399+NCHPhS+9rWvhYcffji89a1vzd4JcuWVV85G/QFmLPZ0rFmzJnQ6nWw/lmvXrtUDAsP6qO3WrVvD8uXLp/ZvuOGGrFy9enW44447wnve857sXSDxg/3f//3f4Xd/93fDt7/97bBw4cL+1hzgGMUXH04Gj0ntdjuMj49bhRYSGOnGl3MMkThME596iZNPzf8AZkPs4YhDLQcGkLGxsTAxMSF8QIL798CfdgFILQaMer2eBY4olrVaTfCARPR8AHO6ByQOtZTLZcEDEt6/ZzznA6AoYuAQOiA9wy4AQFLCBwCQlPABACQlfAAASQkfAEBSwgcAkJTwAQAkJXwAAEkJHwBAUsIHAJCU8AEM9dorzWYzK4HiED6AodRoNLJl76vValbGfaAYrGoLDJ3Y0xEDR6fTmToWl72fmJiwEBwMqZncv/V8AEOn1Wr1BI+o3W6H8fHxgdUJ6B/hAxg6lUoljI72/nmKPR/lcnlgdQL6R/gAhk4cWqnX61ngiGJZq9UMuUBBmPMBDPXcjzjUEns8BA8YbjO5f89LViuAGYqBQ+iA4jHsAgAkJXwAAEkJHwBAUsIHAJCU8AEAJCV8AABJCR8AQFLCBwCQlPABACQlfADH9frzZrOZlQDTJXwAx6TRaIRSqRSq1WpWxn2A6bCwHDBjsacjBo5OpzN1LK48OzExYS0WmKP2zuD+recDmLFWq9UTPKJ2u52tQAtwNMIHMGOVSiWMjvb++Yg9H+VyeWB1AvJD+ABmLA6t1Ov1LHBEsazVaoZcgGkx5wM4rrkfcagl9ngIHjC37Z3B/XtesloBhRMDh9ABzJRhFwAgKeEDAEhK+AAAkhI+AICkhA8AICnhAwBISvgAAJISPgCApIQPACAp4QMASEr4gDmw/kqz2cxKgGEgfECBNRqNUCqVQrVazcq4DzBoVrWFgoo9HTFwdDqdqWNjY2NhYmLCYnDAQO/fej6goFqtVk/wiNrtdhgfHx9YnQAi4QMKqlKphNHR3o947Pkol8sDqxNAJHxAQcWhlXq9ngWOKJa1Ws2QCzBw5nzAHJj7EYdaYo+H4AEUcs5HHFN+3/veF5YuXRpOPPHE8LznPS/cdNNNYcgyDswZMXAsW7ZM8ACGxrx+/wtvvvnmcNttt4W//uu/Duedd17YunVreNvb3paloWuvvbbfPw4AmOvh47777gtXXHFFWLlyZbZ/7rnnhi9+8Yvhhz/8Yb9/FACQQ30fdnn1q18d7rnnnvCTn/wk2/+Xf/mX8P3vfz9cfvnlhzx///792TjRgRsAUFx97/l473vfmwWIF77whdns+jgH5C/+4i/C1VdffcjzN2zYEG688cZ+VwMAmCs9H1/+8pfD5z//+fCFL3whbN++PZv78dGPfjQrD2X9+vXZzNjJbdeuXf2uEgBQ5EdtlyxZkvV+rFu3burYhz70ofC5z30u/PjHPz7q93vUFgDyZ6CP2j799NOHfKviwa95BgDmpr7P+Vi1alU2x+Occ87JHrV96KGHwsc+9rHw9re/vd8/CgDIob4Puzz55JPZS8Y2bdoU9uzZExYtWhSuuuqq8P73vz/Mnz//qN9v2AUA8mcm92+vVwcA8j3nA5j52ivNZjMrAeYC4QMGqNFohFKpFKrValbGfYCiM+wCAxJ7OmLgOPBJsPhk2MTEhEXggNwx7AI50Gq1nvEIenwj8Pj4+MDqBJCC8AEDUqlUDvlOnHK5PLA6AaQgfMCAxKGVer2eBY4olrVazZALUHjmfMAQzP2IQy2xx0PwAPJqJvfvvr/hFJiZGDiEDmAuMewCACQlfAAASQkfAEBSwgcAkJTwAQAkJXwAAEkJHwBAUsIHAJCU8AEAJCV8wDRef95sNrMSgOMnfMARNBqNUCqVQrVazcq4D8DxsbAcHEbs6YiBo9PpTB2LK89OTExYiwXgOO7fej7gMFqtVk/wiNrtdrYCLQDHTviAw6hUKmF0tPcjEns+yuXywOoEUATCBxxGHFqp1+tZ4IhiWavVDLkAHCdzPmAacz/iUEvs8RA8AI7//j3viF8FssAhdAD0j2EXACAp4QMASEr4AACSEj4AgKSEDwAgKeEDAEhK+AAAkhI+AICkhA8AICnhAwBISvigEOL6K81mMysBGG7CB7nXaDRCqVQK1Wo1K+M+AMPLqrbkWuzpiIGj0+lMHRsbGwsTExMWgwMY0vu3ng9yrdVq9QSPqN1uh/Hx8YHVCYAjEz7ItUqlEkZHe3+NY89HuVweWJ0AODLhg1yLQyv1ej0LHFEsa7WaIReAIWbOB4WZ+xGHWmKPh+ABMNz373nJagWzKAYOoQMgHwy7AABJCR8AQFLCBwCQlPABACQlfAAASQkfAEBSwgcAkJTwAQAkJXwAAEkJHwBAUsIHydZeaTabWQnA3CZ8MOsajUYolUqhWq1mZdwHYO6alfDx85//PPzhH/5hOOOMM8KJJ54YLrjggrB169bZ+FEMudjTsWbNmtDpdLL9WK5du1YPCMAc1vdVbX/1q1+F17zmNWH58uXhW9/6Vnj2s58dWq1WOO200/r9o8iBeO0ng8ekdrsdxsfHrUILMEf1PXzcfPPNYcmSJeGzn/3s1LGlS5f2+8eQE5VKJYyOjvYEkLGxsVAulwdaLwAKNOzyta99LVx00UXhD/7gD8JznvOc8LKXvSx86lOfOuz5+/fvD3v37u3ZKI7Yu1Gv17PAEcWyVqvp9QCYw0a63W63n//ChQsXZuUNN9yQBZAHH3wwXHfddeH2228Pq1evfsb5H/zgB8ONN974jONPPPFEOOWUU/pZNQYozvGIQy2xx0PwACie2Hlw6qmnTuv+3ffwMX/+/Kzn47777ps6du2112Yh5P777z9kz0fcDqx8HLYRPgCgmOGj78MuZ599dnjxi1/cc+xFL3pR+NnPfnbI8xcsWJBV8sANACiuvoeP+KTLI4880nPsJz/5SfZ+BwCAvoePP/3TPw0PPPBA+Mu//MtsjP8LX/hCNuFw3bp1/f5RAEAO9T18vPKVrwybNm0KX/ziF8P5558fbrrppnDLLbeEq6++ut8/CgDIob5POE05YQUAGA4DnXAKAHAkwgcAkJTwAQAkJXwAAEkJH0y9/rzZbFrqHoBZJ3wQGo1G9hK4arWalXEfAGaLR23nuNjTEQPHwUveT0xMWAAOgGnzqC3T1mq1eoJH1G63s7fTAsBsED7muEqlEkZHe38NYs9HuVweWJ0AKDbhY46LQytx7Z0YOKJY1mo1Qy4AzBpzPpia+xGHWmKPh+ABwGzev+fN+N9OIcXAIXQAkIJhFwAgKeEDAEhK+AAAkhI+AICkhA8AICnhAwBISvgAAJISPgCApIQPACAp4QMASEr4yNHaK81mMysBIM+EjxxoNBqhVCqFarWalXEfAPLKqrZDLvZ0xMDR6XSmjsVl7ycmJiwEB0Au7996PoZcq9XqCR5Ru90O4+PjA6sTABwP4WPIVSqVMDrae5liz0e5XB5YnQDgeAgfQy4OrdTr9SxwRLGs1WqGXADILXM+cjT3Iw61xB4PwQOAYTOT+/e8ZLXiuMTAIXQAUASGXQCApIQPACAp4QMASEr4AACSEj4AgKSEDwAgKeEDAEhK+AAAkhI+AICkhA8AICnho8/rrzSbzawEAA5N+OiTRqMRSqVSqFarWRn3AYBnsqptH8Sejhg4Op3O1LGxsbEwMTFhMTgA5oS9M7h/6/nog1ar1RM8ona7HcbHxwdWJwAYVsJHH1QqlTA62vufMvZ8lMvlgdUJAIaV8NEHcWilXq9ngSOKZa1WM+QCAIdgzkef537EoZbY4yF4ADCX7J3B/XteslrNATFwCB0AcGSGXQCApIQPACAp4QMASEr4AACSEj4AgKSEDwAgKeEDAChW+Pjwhz8cRkZGwvXXXz/bPwoAmOvh48EHH8xeM/6Sl7xkNn8MAJAjsxY+9u3bF66++urwqU99Kpx22mlhWF5/3mw2sxIAKFj4WLduXVi5cmW47LLLjnje/v37s/fBH7jNhkajEUqlUqhWq1kZ9wGAgoSPO++8M2zfvj1s2LDhqOfGc+JCNJPbkiVL+l6f2NOxZs2a0Ol0sv1Yrl27Vg8IABQhfOzatStcd9114fOf/3xYuHDhUc9fv359tgLe5Ba/v99ardZU8JjUbrezFWgBgLT6vqrttm3bwp49e8LLX/7ynhv9li1bwq233poNs4yNjU19bcGCBdk2myqVShgdHe0JILEO5XJ5Vn8uAJCg5+PSSy8NDz/8cNixY8fUdtFFF2WTT+M/Hxg8UonL3Nfr9amfHcv4FE48DgDkvOfj5JNPDueff37PsZNOOimcccYZzzie0jXXXBNWrFiRDbXEHg/BAwAKEj6GWQwcQgcAzIHwce+996b4MQBADljbBQBISvgAAJISPgCApIQPACAp4QMASEr4AACSEj4AgKSEDwAgKeEDAEhK+AAA5vbaLt1uNyv37t076KoAANM0ed+evI/nKnw8+eSTWblkyZJBVwUAOIb7+KmnnnrEc0a604koCXU6nfDYY4+Fk08+OYyMjPQ9lcVQs2vXrnDKKaeEoil6++ZCG7Uv/4reRu3Lv72z1MYYJ2LwWLRoURgdHc1Xz0es8Gwvex//Yxf1l2outG8utFH78q/obdS+/DtlFtp4tB6PSSacAgBJCR8AQFJzKnwsWLAgfOADH8jKIip6++ZCG7Uv/4reRu3LvwVD0Mahm3AKABTbnOr5AAAGT/gAAJISPgCApIQPACCpwoWPjRs3hnPPPTcsXLgwXHLJJeGHP/zhEc//u7/7u/DCF74wO/+CCy4I3/zmN0NR2nfHHXdkb4k9cIvfN6y2bNkSVq1alb0dL9Z18+bNR/2ee++9N7z85S/PZm2Xy+WszcNspm2M7Tv4Gsbt8ccfD8Nmw4YN4ZWvfGX2duLnPOc54corrwyPPPLIUb8vT5/BY2ljnj6Ht912W3jJS14y9fKpV73qVeFb3/pWYa7fsbQxT9fvUD784Q9ndb7++uvDMF3HQoWPL33pS+GGG27IHiHavn17uPDCC8OKFSvCnj17Dnn+fffdF6666qpwzTXXhIceeij7QxK3f/u3fwtFaF8UP1y/+MUvprZHH300DKunnnoqa1MMWNOxc+fOsHLlyrB8+fKwY8eO7MP1R3/0R+Guu+4KRWnjpHiDO/A6xhvfsPne974X1q1bFx544IFw9913h9/85jfhda97Xdbmw8nbZ/BY2pinz2F8u3S8WW3bti1s3bo1VKvVcMUVV4Qf/ehHhbh+x9LGPF2/gz344IOhVqtlYetIBnIduwVy8cUXd9etWze13263u4sWLepu2LDhkOe/8Y1v7K5cubLn2CWXXNJdu3Zttwjt++xnP9s99dRTu3kUfzU3bdp0xHPe8573dM8777yeY29605u6K1as6Baljc1mMzvvV7/6VTdv9uzZk9X9e9/73mHPydtn8FjamOfPYXTaaad1P/3pTxfy+k2njXm9fk8++WS3Uql077777u5rX/va7nXXXXfYcwdxHQvT8/G///u/WZK97LLLetaJifv333//Ib8nHj/w/Cj2JBzu/Ly1L9q3b18olUrZIkJHS/d5k6frd7xe+tKXhrPPPjv83u/9XvjBD34Q8uCJJ57IytNPP72w13A6bczr57Ddboc777wz69WJQxNFvH7TaWNer9+6deuynuGDr8+wXMfChI//+q//yn6RzjzzzJ7jcf9w4+Px+EzOz1v7XvCCF4TPfOYz4atf/Wr43Oc+l60Y/OpXvzrs3r07FMHhrl9csfF//ud/QhHEwHH77beHr3zlK9kW//gtW7YsG3YbZvF3LQ6DveY1rwnnn3/+Yc/L02fwWNuYt8/hww8/HH77t387m0f1jne8I2zatCm8+MUvLtT1m0kb83b9ohio4t+IOEdpOgZxHYduVVv6Jyb5A9N8/MC86EUvysYAb7rppoHWjemJf/jiduA1/OlPfxo+/vGPh7/9278Nw/z/uuJ48fe///1QVNNtY94+h/H3Lc6hir06f//3fx9Wr16dzXU53M05j2bSxrxdv127doXrrrsum5M0zBNjCxM+fud3fieMjY2FX/7ylz3H4/5ZZ511yO+Jx2dyft7ad7ATTjghvOxlLwvj4+OhCA53/eLksBNPPDEU1cUXXzzUN/U/+ZM/CV//+tezJ3vi5L4jydNn8FjbmLfP4fz587Mnx6JXvOIV2aTFT3ziE9nNtijXbyZtzNv127ZtW/YQQnwKcFLsNY+/q7feemvYv39/di8Z9HUszLBL/GWKv0T33HPP1LHYPRb3DzeWF48feH4U0+KRxv7y1L6DxV/A2N0Yu/KLIE/Xr5/i/2MbxmsY59DGm3Lswv7ud78bli5dWrhreCxtzPvnMP6diTesIly/Y2lj3q7fpZdemtUv/p2Y3C666KJw9dVXZ/98cPAY2HXsFsidd97ZXbBgQfeOO+7o/vu//3t3zZo13Wc961ndxx9/PPv6W97ylu573/veqfN/8IMfdOfNm9f96Ec/2v2P//iP7gc+8IHuCSec0H344Ye7RWjfjTfe2L3rrru6P/3pT7vbtm3rvvnNb+4uXLiw+6Mf/ag7rLOzH3rooWyLv5of+9jHsn9+9NFHs6/HtsU2TvrP//zP7m/91m913/3ud2fXb+PGjd2xsbHut7/97e6wmmkbP/7xj3c3b97cbbVa2e9lnLE+Ojra/cd//MfusPnjP/7j7KmAe++9t/uLX/xianv66aenzsn7Z/BY2pinz2Gsd3xyZ+fOnd1//dd/zfZHRka63/nOdwpx/Y6ljXm6fodz8NMuw3AdCxU+or/6q7/qnnPOOd358+dnj6Y+8MADPRdg9erVPed/+ctf7j7/+c/Pzo+PbX7jG9/oFqV9119//dS5Z555Zvf3f//3u9u3b+8Oq8nHSg/eJtsUy9jGg7/npS99adbG5z73udljccNspm28+eabu8973vOyP3ann356d9myZd3vfve73WF0qHbF7cBrkvfP4LG0MU+fw7e//e3dUqmU1fXZz35299JLL526KRfh+h1LG/N0/aYbPobhOo7E/5m9fhUAgILO+QAA8kH4AACSEj4AgKSEDwAgKeEDAEhK+AAAkhI+AICkhA8AICnhAwBISvgAAJISPgCApIQPACCk9H8ApPDWCIYUyfUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "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": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAny0lEQVR4nO3dD3BU1fXA8ZMQg9QKChGBSUio9U/8h4aoU3WmgrZOSql2prU6qSaaTtShorUjykzVMtpGHEdBpERtakxQE20LdWxTihqkiqAEqdS2CkpiqmKSjiWANRiyvzm3ffvbDZuQhLf3/ft+ZtbHvt2Q93jZvOO955ybEYvFYgIAAGBJpq1vBAAAoAg+AACAVQQfAADAKoIPAABgFcEHAACwiuADAABYRfABAACsIvgAAABWZYnP9PX1yYcffihHHnmkZGRkeH04AABgCLRn6e7du2XKlCmSmZkZrOBDA4+8vDyvDwMAAIxAe3u75ObmBiv40BEP5+DHjh3r9eEAAIAh6O7uNoMHzn08UMGHM9WigQfBBwAAwTKUlAkSTgEAgFUEHwAAwCqCDwAAYBXBBwAAsIrgAwAAWEXwAQAArCL4AAAAVhF8AAAAqwg+AACAVQQfAABEyKZNm2TWrFlm6xWCDwAAIqSurk6am5ulvr7es2Pw3douAADAXW1tbdLV1WXWXWlsbDT7GhoapKysTGKxmOTk5Eh+fr7YQvABAEDIFRQUHLDwW2dnp8yYMSO+X4MQW5h2AQAg5FasWCFZWVlJQYaz1f36uk2MfAAAEHKlpaVSWFiYNNLh2LhxoxQVFVk9HkY+AACIkMzMzKStJ8fg2XcGAADWTJw4USZNmmRGP6qrq81Wn+t+2zJiNjNMhqC7u1vGjRsnu3btkrFjx3p9OAAAhEZPT49kZ2ebpFO9/e/bt09Gjx5t/f5NzgcAABExOiHQ0ADErcBjuJh2AQAghF1E/YzgAwCAEHYR9TOmXQAACGEXUT8j+AAAIIRdREM17bJu3TqZM2eOTJkyxfzjrlq1asD3XnfddeY9ixcvPtTjBADA1/zWRTRUwcfevXtl+vTpsmzZskHft3LlStmwYYMJUgAAiEIXUe0Wmoru19cxwmmXkpIS8xjMBx98IDfccIOsXr1aZs+ePdxvAQBAoGn30L6+vvgWac750H/kK6+8Um655RY55ZRThtTwRB+JTUoAAAhyF9G8vDypqKiQmpoaaW9v96SLaKSCj0WLFpm5rXnz5g3p/VVVVbJw4UK3DwMAAOtyc3OltbU13kW0srLS1S6iYeFqn4+WlhZZsmSJ1NbWxjN9D2bBggWmFavz0AgRAICg0kDDuQd62UU0MsHHn//8Z+no6JCpU6ea0Q99aN3zj3/846QSpER6UbQHfOIDAACEl6vBh+Z6vPnmm7Jly5b4Q6tdNP9Dk08BAIi6TbReH37Ox549e2T79u3x5zt27DBBxvjx482Ix4QJE5Lef9hhh5nkmxNPPNGdIwYAICSt14uLiyWKhh18aKQ2c+bM+PObb77ZbLV9rOZ6AACAZLReT5YR81mvVy21HTdunEk+Jf8DABAGiUUYGRkZJuBwtg6f3Y7Tev9mVVsAANKcm0Hr9WQEHwAAHCQ341DRej0Zq9oCAGAxNyOT1usEHwAAOBJ7Ujl5Gp2dnTJjxoxDzs2g9fr/I/gAAOB/NPeivLxcent7U+ZmHEpVJ63X/x/BBwAA/6O5F4WFhUkjHYm5GUVFRYf09ycGGhkRbr1OwikAACloTkbiNqg2+bCjarD/RQEAcJmTm6GjH9XV1Warz4Oam1HnYtWOW2gyBgBAPz09PfHcDL1NBi03oy2haqekpMQs+qrBU1NTU9o6qg7n/k3OBwAAIcvNKEhj1Y4bmHYBACBkVvi8oyojHwAAhExpmqt2DhUjHwAAhFimD6t2/HMkAAAgElU7VLsAABBSPRardqh2AQAA4teqHaZdAACAVQQfAADAKoIPAABgFcEHAACwiuADAABYRfABAACsIvgAALhu06ZNMmvWLLMF+iP4AAC4rq6uTpqbm6W+vt7rQ4EPEXwAAFzR1tYmLS0tsnnzZmlsbDT7GhoazHPdr68jWVRHiAg+AACuKCgokOLiYrOGSGdnp9mnW32u+/X1KN9wU6mL6AgRwQcAwBUrVqyQrKz/rtrhLBvmbHW/vh7lG66jjREiFpYDALhHb6A60tHfc889Z1ZU1fVFSkpKpKOjw6yu2tTUZAKUnJwcyc/PlyjIyMhI+rOev7N1+OzWPCQsLAcA8FRmZqb09fXFt9/85jcPuPk6UzJBvuGOxIoVK6S8vFx6e3tTjhDV1tZK2DHtAgBwjY5m6AiHBhXV1dVmq88ffPDBIU3JREFpaals3Lgx5Wu6X18PO0Y+AACuyc3NldbWVsnOzjYjHJWVlbJv3z6zlPt5552XckpGb7hFRUUSRZn9RoiigpEPAICrNNBwplZ0q88T6Y02ceslrypvJg4wQqT7o4CRDwCA1RtuXl6eVFRUSE1NjbS3t3t6w02svNFyYD+MEEUB1S4AAGt6enriN1y9/Xhxw9VS1q6uLipvXEa1CwDAlxIDjVRTMjY4zc6cY4hy5Y1XvJ9wAwDAh83QkD6MfAAAIkVLWQsLC6m88RAjHwCAyPJT5U2U8K8NAIicqJe6eo1qFwBAJPmh8iZMqHYBACAAlTdRxbQLAACwiuADAABYRfABAPAtr9ZeQXoRfAAAfCtx7RUHAUnwkXAKAPCVxLVXGhsbzb6GhgYpKyszVSnLli3zZDE4uIdSWwCArzjrrTh/Hug2xWJw/jKc+zfTLgAA36+9koqzGJyOfiQuFgf/G3bwsW7dOpkzZ45MmTLFRKSrVq2Kv/b555/LrbfeKqeddpocccQR5j1XXXWVfPjhh24fNwAgxGuv6BorB8NicBEKPvbu3SvTp083c279ffrpp7J582a5/fbbzfa3v/2tvP322/Ktb33LreMFAETIUNZe0UBFAxaEOOG0pKTEPFLRuZ41a9Yk7XvooYfk7LPPlvfff1+mTp068iMFAERu7ZW8vDypqKiQmpoaee+99+Rf//qXCUT6+vriWwRP2qtdNPFEp2eOOuqoAXvr6yMxYQUAEG25ubnS2toaX3ulsrLSBB/nn39+UkDS3t7OYnABlNbg47PPPjM5IFdcccWAma9VVVWycOHCdB4GACAEa68cd9xxBwQkLAYXTGmrdtHk08suu8wkBC1fvnzA9y1YsMCMjjgPjWIBAEhFAw2nFJfF4IIrM52BhzaK0RyQwep99QdHX098AACiha6l0ZKZrsBj27Zt8vzzz8uECRPc/hYAgJAFGanaqCO8hp3zsWfPHtm+fXv8+Y4dO2TLli0yfvx4mTx5snznO98xZbbPPfec7N+/X3bu3Gnep6/rPB0AILoSg4xjjjlm0DbqdC0Nr2G3V1+7dq3MnDnzgP36w/LTn/5Upk2blvLr9IftggsuOOjfT3t1AAjvWi3aqqGjo8NUqOi2fxv1/u3UfbYCCFy6fw975EMDiMF+GPhBAQAkSmx97iSLamv0VPeOxK6ltbW1Vo8T9rC2CwDA+lotznbUqFEpv4aupeFG8AEA8Gytlscff3zIbdQRHlxlAIA1/YMMLUbQNuq6Om11dbXZ6nO6loYbwQcARJit/hrOWi39gwxdBV27lurIyLXXXmu2+lzbqyO80r62CwAgGKWvxcXFVtdqSdUana6l0UDwAQARLn212V+j/1otBBnRRfABABEzUOmrToU4aJuAdCLnAwAiZrDSV92vrwPpRPABABEzWOnrwfprsAAc3EDwAQARNtz+GiwABzcQfABAyAxldGKg0tdU/TU0QbWlpcUsGpqYoKrPdb++DqR1Ybl0Y2E5ADg08+bNk6VLl5rtkiVLBnxfT09PvPRVbwWpSl8Tk1KdP7MAHA71/s3IBwB4zI08ipGMTmig4QQWg5W+kqAKt1FqCwAhaPSVzvJZTUAtLCxM+rsSE1SLiopG9Pciuhj5AAAPuJ1HYWt0ggXg4AZGPgDAA26PVKR7dMJJUM3Ly5OKigqpqamR9vZ2FoDDiBB8AIAHdCSivLxcent7U45U1NbWjvjv1lGJvr6++Nbm2izAUDBuBgABa/TlRvnsSAw1QRU4GIIPAPDYUPMoDlYV44xOsDw9/I7gAwA8MtyRiqF0F3VzdIJW6kgXcj4AwCNDyaPQqpeuri7zemJVTFlZmckRycnJkfz8fN+WAAOp0OEUAHzMdnfRxGCnpKREOjo6zEhMU1NT2oMdBNtw7t+MfABARKtibDcrAxzkfABAxKpiBkMrddhA8AEAAWGju6jtYAfRRPABAD6X7v4dA6GVOtKFnA8A8Dnb3UVppY50o9oFAHCAnp6eeLCjtwlaqeNgqHYBABySxECDVupwGxN5AADAKoIPAABgFcEHAHiAdVMQZQQfAOCBoSwSB4QVCacAYImXi8QBfkLwAQCWsG4K8F9MuwCAJaybAvwXIx8AYImui1JYWJg00pG4bkpRUZEnxwXYxsgHAHiAdVMQZfzUA4AlWlZ7ww03yIQJE6wvEgf4CdMuAGCxvHb9+vUyd+5cWbp0qZVF4gA/IvgAAMvltc8884xcc801lNcisgg+ACCNKK8FDkTOBwCksR065bXAgQg+ACCN7dC1vFbLaFPR/fo6EDVMuwCApXboWlbb19cX3wJRRfABAGnO19AyWi2nzcvLk4qKCqmpqZH29nbKaxFZTLsACAQvlqB3K18jNzdXWltbzTTLtddea7b6XPcDUTTs4GPdunUyZ84cmTJlivk/gVWrViW9rh/MO+64QyZPnixjxoyRiy66SLZt2+bmMQOIIC+WoHczX0P7eDijJ7qlrweibNjBx969e2X69OmybNmylK/fe++98uCDD5rOffrhPOKII+Tiiy+Wzz77zI3jBRCxnIuWlhbZvHlzUs6FPtf9+rqtEZaB2qF7MSIDRC7no6SkxDxS0VGPxYsXy09+8hO55JJL4v+3cuyxx5oRkssvv/zQjxhAZHjRIyNxhKW4uPig+Rr93w/Acs7Hjh07ZOfOnWaqxTFu3Dg555xz5NVXX3XzWwGIAFs9MgYbYfn444/lpZdeSsrX0Oe63+aIDBAmrla7aOChdKQjkT53Xuuvp6fHPBzd3d1uHhKAALO1BP1wRlj09RNPPHHI7wfgw2qXqqoqMzriPHRoEwBsLkE/3BEWupYCh8bVT7HOiyodjkykz53X+luwYIHs2rUr/tC5VABwODkX6VyCfrhVLXQtBXw07TJt2jTzS+GFF16QM844Iz6Noh/G66+/PuXXaLkZJWcADtYjIzs728oS9HQhBXw48rFnzx7ZsmWLeThJpvrn999/3/xiuOmmm+Tuu++WZ599VrZu3SpXXXWV6Qly6aWXpuP4AUTAYD0y3Cp1Hc4IiyaUah7bhAkT4tMvutXnup+EU2BwGbFhZkWtXbtWZs6cecB+XeugtrbWzHveeeed8sgjj8i///1vOf/88+UXv/iFnHDCCUP6+3WkRHM/dApm7Nixwzk0ABE0b948Wbp0qdkuWbLkkP4uTX53Rlj0d9lAIyxOIOT8Wd/rbB0knCJquodx/x528JFuBB8AhrPgm/Yd6ujoMCMUTU1NI1rwbbieeOIJKS8vl97e3gNe0xEQ/R8x8j4QNd0EHwDCzA8jD9rTI1UJsPb5cKsEGAiS4dy/PS+1BYDh8lOpazpLgIGw4tMCIHD8UOpqowQYCCtXS20BwDavSmNtlwADYcLIB4BA8sPIw2AlwAAGRsIpgMAaamksgPQj4RRAWrjV0MstjDwAwUTwAWDI6urqpLm5Werr60MVxACwi4RTAENu6NXY2Gj2NTQ0mK7GI23olRjEFBcXp+nIAfgVOR8ArDT08rorKQD/3L8Z+QAwKG3Y5bQST9XQS1uJD0VBQcEBAU1nZ2dSl1Cf/b8QgDQh5wOAlYZefupKCsBbBB8ArLQS90NXUgD+QPABwHpDLy/XQ6HSBvAeOR8ArLUSd4KYvLw8qaiokJqaGmlvb7falZRKG8B7VLsACH1XUiptgPSj2gWAbyUGGra6klJpA/gLOR8AQo9KG8BfGPkAEHpaSVNYWJg00pFYaVNUVOTJcQFRxcgHgEjxstIGwH/x6QMQCW6XCwMYOapdAESGF5U2QFR0U+0CAP6otAFwIKZdAACAVQQfAADAKoIPAABgFcEHAOtY3A2INoIPANYlLu4GIHqodgFgReLibo2NjWZfQ0ODlJWVsbgbEDEEHwCsYHE3AA6mXQBYweJuAByMfACwgsXdADgY+QBgHYu7AdHGJx+ANSzuBkCxsBwAq1jcDQgnFpYD4Fss7gaAaRcAAGAVwQfgs7bffjiGKB47AHsIPgCftf32wzFE8dgB2EPCKSIvse13SUmJdHR0mOqLpqYma22//XAMUTx2AN7cvwk+EHlOq2/nz/qRcLaOdH9M/HAMUTx2AN7cv5l2QeT5oe23H44hiscOwBuMfAAisnnz5pRtv1taWqy1/fbDMUTx2AG4g5EPIMBtv/1wDFE8dgD28BsC8Enbbz8cQxSPHYB9TLsAPmr77YdjiOKxAzh0tFcHAtr22w/HEMVjB2AX0y4AACDYwcf+/fvl9ttvl2nTpsmYMWPkuOOOk7vuuos6fwAAkJ5pl0WLFsny5cvl8ccfl1NOOcWs8XD11VebeaB58+a5/e0AAEDUg4/169fLJZdcIrNnzzbPCwoK5KmnnpLXXnvN7W8FAAACyPVpl3PPPVdeeOEFeeedd8zzv/zlL/Lyyy+bNR8GypDXDNnEBwAACC/XRz5uu+02E0CcdNJJMmrUKJMD8rOf/UxKS0tTvr+qqkoWLlzo9mEAAICojHw8/fTT8sQTT8iTTz5pWi5r7sd9991ntqksWLDA1AQ7j/b2drcPCQAAhLnJWF5enhn9mDt3bnzf3XffbRaX+sc//nHQr6fJGAAAwePp2i6ffvrpAes66PRLX1+f298KAAAEkOs5H3PmzDE5HlOnTjWltm+88Ybcf//9cs0117j9rQAAQAC5Pu2ye/du02Rs5cqV0tHRIVOmTJErrrhC7rjjDrPuw8Ew7QIAQPAM5/7NwnKILG2AN3/+fLn33nuluLjY68MBgEDzNOcDCIq6ujppbm6W+vp6rw8FACKFVW0RKW1tbdLV1WVWXW1sbDT7GhoapKyszKw/lJOTI/n5+V4fJgCEGsEHIkXb/Ts0AFGdnZ0yY8aM+H6fzUQCQOgw7YJI0X4zWVlZSUGGs9X9+joAIL0Y+UCkaJv/wsLCpJEOx8aNG6WoqMiT4wKAKGHkA5HlNMPr3xTP6wqcWbNmmS0AhJV/fusClkycOFEmTZpkRj+qq6vNVp/rfq9RgQMgCujzgUjq6ekxTe806VQ/Avv27ZPRo0d7XoFTUlJimvNpINTU1EQFDoDAGM79m5wPRFJioKE3fa8CD0UFDoCoYdoF8DhvgwocAFFD8AF4nLehFThaaZOK7tfXASBMmHYBfNQ5VStv+vr64lsACCOCD8AHeRtOBU5eXp5UVFRITU2NtLe3+6ICBwDcRvABuEDzMsrLy6W3tzdl3kZtbe2gX5+bmyutra3xCpzKykpPK3AAIJ0IPgCfdE71UwUOAKQTCadABDqnAoCf8NsRiEDnVADwEzqcAiHtnAoANtHhFPAIeRsAcHBMuwAAAKsIPgAAgFUEHwAAwCqCDwAAYBXBBwAAsIrgAwAAWEXwAd/btGmTzJo1y2wBAMFH8AHfq6urk+bmZqmvr/f6UAAALqDJGHypra1Nurq6TKOuxsZGs6+hoUHKyspM59CcnBzJz8/3+jABACNA8AFfKigoiP9ZAxDV2dmZtGqsz1YGAAAMEdMu8GUux4oVKyQrKyspyHC2ul9fBwAEE8EHfJnLUVpaKhs3bkz5dbpfXwcABBPTLvB9LkdmZqb09fXFtwCAYCP4gG9zOSZOnCiTJk2SvLw8qaiokJqaGmlvbzf7AQDBRfABqzRXo7y8XHp7e1PmctTW1sbfm5ubK62trZKdnW0ClcrKStm3bx/L1ANAwBF8wCrN1SgsLEwa6UjM5SgqKkralxhoaABC4AEAwUfCKTyjORyJWwBANPBbH9Y5uRw6+lFdXW22+pxcDgCIhoyYzzo1dXd3y7hx42TXrl0yduxYrw8HadLT0xPP5dAfQa9yObTHyPz58+Xee++V4uJi698fAMJiOPdvRj7gCQ00nGoXL3M5WDcGAOwj4RSRw7oxAOAtgg9EDuvGAIC3mHZB5LBuDAB4i5EPRM5we40AANzFyAcijV4jAGAfv3ERSfQaAQDv0OcDkeWXXiMAEAbDuX+T84HIYt0YAPAG0y4AACD4wccHH3wg3//+92XChAkyZswYOe2000wbawAAANenXT755BM577zzZObMmdLU1CTHHHOMbNu2TY4++mi3vxUAAAgg14OPRYsWSV5enjz22GPxfdOmTXP72wAAgIByfdrl2WefNauDfve73zVli2eeeaY8+uijg1YcaIZs4gMAAISX68HHe++9J8uXL5fjjz9eVq9eLddff73MmzdPHn/88ZTvr6qqMqU5zkNHTQAAQHi53udD+yboyMf69evj+zT4eP311+XVV19NOfKhD4eOfGgAQp8PAADC2efD9ZGPyZMny8knn5y0T9fReP/991O+X3sr6EEmPgAAQHi5Hnxopcvbb7+dtO+dd96R/Px8t78VAAAIINeDjx/96EeyYcMG+fnPfy7bt2+XJ598Uh555BGZO3eu298KAAAEkOvBx1lnnSUrV66Up556Sk499VS56667ZPHixWYZcwAAABaWAwAAwU44BQAAGAzBBwAAsIrgAwAAWEXwAQAArCL4AAAAVhF8AAAAqwg+AACAVQQfAADAKoIPAABgFcEHAACwiuAj4DZt2iSzZs0yWwAAgoDgI+Dq6uqkublZ6uvrvT4UAACGJGtob4OftLW1SVdXl2RkZEhjY6PZ19DQIGVlZaLrBObk5Eh+fr7XhwkAQEoEHwFUUFAQ/7MGIKqzs1NmzJgR3++zxYoBAIhj2iWAORgrVqyQrKyspCDD2ep+fR0AAL8i+AhgDkZpaals3Lgx5Wu6X18HAMCvmHYJeA5GZmam9PX1xbcAAPgdwUdAczAmTpwokyZNkry8PKmoqJCamhppb283+wEA8DOCDxdojkV5ebn09vamzMGora11/Xvm5uZKa2urZGdnm4CnsrJS9u3bJ6NHj3b9ewEA4CaCDxdojkVhYWHSSEdiDkZRUVFavm9ioKEBCIEHACAISDh1meZeJG4BAEAy7pAu52Do6Ed1dbXZ6nNyMAAASJYR81k3qu7ubhk3bpzs2rVLxo4dK0HS09MTz8HQf1ZyMAAAUdE9jPs3OR8uIgcDAICDY9oFAABYRfARQDZauAMAkC4EHwG8sdto4Q4AQLqQ8zHIjb24uFj8wssW7gAAuIlql3439pKSEuno6DAlsk1NTb65sTtt250/63E5W4fPLiUAIEK6qXbx/9osQWjhDgBAOpDz8b8bu97AVaobu77uhxbu2qo9Fd2vrwMAEAQEHwG8sdPCHQAQZNy9AnRjp4U7ACAMyPnod2PPy8uTiooKqampkfb2dl/d2HNzc6W1tTXewr2yspIW7gCAwKHaJQFrswAAkP77t//mFjxsLLZ169Z4tYtba7MEoWkZAAA2EXykuWMo3UgBAEgW2WmXdDYWC0LTMgAAvLp/Rzb4SGfHULqRAgCippucD28biwWhaRkAAF6J7MiH2rx5c1ILdUdLS4sUFRX59u8GAMBvGPnwUWMxPzctAwDAC5G+I6azYyjdSAEASC1S0y7aa2P+/Ply7733SnFxcdobi9G0DAAQFd1Muwy954YGA243FrPxdwMAEFShDz6054YmeWoCaGNjo9nX0NBgnut+ff1Q0cUUAIChC/3CcgUFBfE/O6MQnZ2dSZUohzrzlDii4kznAACAiI58pKvnho0RFQAAwijtCaf33HOPLFiwQG688UZZvHixJwmn6ei5QRdTAAB8mHD6+uuvy8MPPyynn366+IGbPTfoYgoAwMikLfjYs2ePlJaWyqOPPipHH320eCkdPTf03DZu3JjyNd2vrwMAAIsJp3PnzpXZs2fLRRddJHffffegvTD0kThs47bc3FxpbW2N99yorKx0teeGjqT09fXFtwAAwPLIh5N4WVVVddD36nt0jsh55OXlpeOQDtpzYyTlsnQxBQDAByMf7e3tJrl0zZo1cvjhhx/0/ZqMevPNNyeNfKQrAHG7XDbdIyoAAISR69Uuq1atkm9/+9syatSo+L79+/ebm7NOS+gUS+JrXq5qq+WwXV1d5thKSkqko6PDjFo0NTWZ5NGcnBzJz89P6zEAABAGw7l/ux587N69+4AeF1dffbWcdNJJcuutt8qpp5466NfbDD4olwUAwB3DuX+7Pu1y5JFHHhBgHHHEETJhwoSDBh62aTlseXm59Pb2piyXra2t9fgIAQAIn9C3Vx+MlsMWFhambECm5bIjbUAGAAA8Dj7Wrl0rfke5LAAAdoR+bZeDoVwWAICQre0yXDYTTh1ageOUy+o/B+WyAAAEKOE0iBIDjVQNyAAAgHsiP+0CAADsIvgAAABWEXwAAACrCD4AAIBVBB8AAMAqgg8AAGAVwQcAALCK4AMAAFhF8AEAAKwi+AAAAFYRfAAAAKsIPgAAgFUEHwAAwCqCDwAAYBXBhws2bdoks2bNMlsAADA4gg8X1NXVSXNzs9TX13t9KAAA+F6W1wcQVG1tbdLV1SUZGRnS2Nho9jU0NEhZWZnEYjHJycmR/Px8rw8TAADfIfgYoYKCgvifNQBRnZ2dMmPGjPh+DUIAAEAypl1GaMWKFZKVlZUUZDhb3a+vAwCAAzHyMUKlpaVSWFiYNNLh2LhxoxQVFXlyXAAA+B0jHy7IzMxM2gIAgIFxtzwEEydOlEmTJpnRj+rqarPV57ofAACklhHzWVZkd3e3jBs3Tnbt2iVjx44Vv+vp6ZHs7GyTdKr/lPv27ZPRo0d7fVgAAPj2/k3OxyFKDDQ0ACHwAABgcEy7AAAAqwg+AACAVQQfAADAKoIPAABgFcEHAACwiuADAABYRfABAACsIvgAAABWEXwAAACrCD4AAIBVvmuv7iw1oz3iAQBAMDj37aEsGee74GP37t1mm5eX5/WhAACAEdzHdYG5QK1q29fXJx9++KEceeSRZqE2t6IxDWba29sDsVLuSHCOwRf281OcYziE/RzDfn7pOkcNJzTwmDJlimRmZgZr5EMPODc3Ny1/t/4Dh/UHycE5Bl/Yz09xjuEQ9nMM+/ml4xwPNuLhIOEUAABYRfABAACsikTwMXr0aLnzzjvNNqw4x+AL+/kpzjEcwn6OYT8/P5yj7xJOAQBAuEVi5AMAAPgHwQcAALCK4AMAAFhF8AEAAKwKTfCxbNkyKSgokMMPP1zOOeccee211wZ9/zPPPCMnnXSSef9pp50mf/jDHyRM51hbW2s6xCY+9Ov8at26dTJnzhzTGU+PddWqVQf9mrVr10pRUZHJ1v7yl79sztnPhnuOen79r6E+du7cKX5UVVUlZ511lulOPHHiRLn00kvl7bffPujXBemzOJJzDNpncfny5XL66afHm0995StfkaamptBcw5GcY9CuYX/33HOPOeabbrpJ/HIdQxF8NDY2ys0332zKhjZv3izTp0+Xiy++WDo6OlK+f/369XLFFVdIRUWFvPHGG+YXiD7++te/SljOUemH6qOPPoo/2traxK/27t1rzkkDrKHYsWOHzJ49W2bOnClbtmwxH6of/OAHsnr1agnLOTr05pZ4HfWm50cvvfSSzJ07VzZs2CBr1qyRzz//XL7+9a+b8x5I0D6LIznHoH0WtcO03qxaWlpk06ZNMmvWLLnkkkvkrbfeCsU1HMk5Bu0aJnr99dfl4YcfNsHWYKxfx1gInH322bG5c+fGn+/fvz82ZcqUWFVVVcr3X3bZZbHZs2cn7TvnnHNi1157bSws5/jYY4/Fxo0bFwsi/bFcuXLloO+ZP39+7JRTTkna973vfS928cUXx8Jyjs3NzeZ9n3zySSyIOjo6zPG/9NJLA74niJ/F4Z5jkD+LjqOPPjr2y1/+MpTXcCjnGNRruHv37tjxxx8fW7NmTeyrX/1q7MYbbxzwvbavY+BHPvbt22ei14suuihpfRh9/uqrr6b8Gt2f+H6lowgDvT+I56j27Nkj+fn5ZvGgg0X1QRO0a3gozjjjDJk8ebJ87Wtfk1deeUWCYteuXWY7fvz40F7HoZxjkD+L+/fvl4aGBjOyo1MTYbyGQznHoF7DuXPnmhHi/tfHD9cx8MFHV1eX+eE59thjk/br84HmxnX/cN4fxHM88cQT5Ve/+pX87ne/kxUrVpjVgs8991z55z//KWEw0DXUlRr/85//SBhowFFdXS2/+c1vzEN/6V1wwQVm2s3v9OdNp8LOO+88OfXUUwd8X9A+iyM5xyB+Frdu3Spf/OIXTT7VddddJytXrpSTTz45VNdwOOcYxGvY0NBgfldontJQ2L6OvlvVFu7QCD4xitcPSmFhoZn7u+uuuzw9NgyN/sLTR+I1fPfdd+WBBx6Q+vp68fv/celc8csvvyxhNdRzDOJnUX/uNJdKR3Z+/etfS1lZmcl3GejmHETDOcegXcP29na58cYbTV6SXxNjAx985OTkyKhRo+Tjjz9O2q/PJ02alPJrdP9w3h/Ec+zvsMMOkzPPPFO2b98uYTDQNdSksDFjxkhYnX322b6/of/whz+U5557zlT3aGLfYIL2WRzJOQbxs5idnW0qyNSMGTNM0uKSJUvMzTYs13A45xi0a9jS0mKKEbQa0KGj5/rz+tBDD0lPT4+5p3h5HQM/7aI/QPqD88ILL8T36ZCYPh9o/k73J75faYQ42Hxf0M6xP/3B02FGHcoPg6BdQ7fo/6n59RpqHq3elHX4+sUXX5Rp06aF7jqO5BzD8FnU3zd6wwrDNRzJOQbtGl544YXm+PT3hfMoLi6W0tJS8+f+gYcn1zEWAg0NDbHRo0fHamtrY3/7299ilZWVsaOOOiq2c+dO8/qVV14Zu+222+Lvf+WVV2JZWVmx++67L/b3v/89duedd8YOO+yw2NatW2NhOceFCxfGVq9eHXv33XdjLS0tscsvvzx2+OGHx956662YX7Oy33jjDfPQH8v777/f/Lmtrc28ruem5+h47733Yl/4whdit9xyi7mGy5Yti40aNSr2xz/+MeZXwz3HBx54ILZq1arYtm3bzM+mZqpnZmbGnn/++ZgfXX/99aYiYO3atbGPPvoo/vj000/j7wn6Z3Ek5xi0z6Ieu1bv7NixI/bmm2+a5xkZGbE//elPobiGIznHoF3DVPpXu3h9HUMRfKilS5fGpk6dGsvOzjZlqRs2bEj6Ry8rK0t6/9NPPx074YQTzPu1ZPP3v/99LEzneNNNN8Xfe+yxx8a+8Y1vxDZv3hzzK6estP/DOSfd6jn2/5ozzjjDnOOXvvQlUw7nZ8M9x0WLFsWOO+4480tu/PjxsQsuuCD24osvxvwq1bnpI/G6BP2zOJJzDNpn8Zprronl5+eb4z3mmGNiF154YfymHIZrOJJzDNo1HErw4fV1zND/pGdMBQAAIIQ5HwAAIFgIPgAAgFUEHwAAwCqCDwAAYBXBBwAAsIrgAwAAWEXwAQAArCL4AAAAVhF8AAAAqwg+AACAVQQfAADAKoIPAAAgNv0fgsQEZaYP8DwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.56763627]\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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.99216691]\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": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.06512069])"
      ]
     },
     "execution_count": 8,
     "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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.5676362738874774\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": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(3.992166911745958)"
      ]
     },
     "execution_count": 10,
     "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": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAABD8ElEQVR4nO3dC5zVc/7H8fdM6YKKrG6mG+VapFJt609GtLQta3eRkGQnhKLcJckqLYpSU4muhKXcJZoQSlfkkrLNKJWKdNVFnf/j8/t2pplp7nNuv3Nez8fj9Ot3zpk5v9+cmfl95vP9fD/fpEAgEBAAAECEJEfqhQAAAAzBBwAAiCiCDwAAEFEEHwAAIKIIPgAAQEQRfAAAgIgi+AAAABFF8AEAACKqvGLMvn37tGbNGlWpUkVJSUnRPhwAAFAM1rN069atqlOnjpKTk/0VfFjgUbdu3WgfBgAAKIVVq1YpJSXFX8GHZTyCB1+1atVoHw4AACiGLVu2eMmD4HXcV8FHcKjFAg+CDwAA/KU4JRMUnAIAgIgi+AAAABFF8AEAACIq5mo+ijud5/fff9fevXujfSi+U65cOZUvX55pzACAqPFd8LF7926tXbtWO3bsiPah+Nahhx6q2rVrq0KFCtE+FABAAvJV8GENyFauXOn99W5NTOziyV/wJcsYWfC2YcMG7+vYuHHjIhvBAACQ0MGHXTgtALF5xPbXO0qucuXKOuSQQ5SVleV9PStVqhTtQwIAJBhf/tnLX+tlw9cPABBNXIUAAEBEEXwAAICIIviIkk2bNmnAgAHezB0AAMJt9erVysjI8LbRRvARJb169dJnn32mG264IdqHAgCIc+PGjVP9+vWVmprqbW0/mgg+ouDNN9/U1q1bve0RRxyhKVOmRPuQAABxavXq1UpLS/Nmixrb9ujRI6oZEF9NtY0XHTt29G5m/Pjx0T4cAEAcW758eXbgEWQdwlesWKGUlJSoHBOZDwAA4ljjfBpKWrPORo0aRe2YCD4AAIhjKSkpGjNmjBdwmHLJyRrdqJFS1q2L2jEldvCxYIGUmuq2EfoGGDlyZK77PvnkE69bq3UcBQAgHLp3767MzExvtktm167qvmyZNGmSoiWxg4+JE6WMjIi9Aa1bt9b8+fNzrbXSu3dv3XrrrV71MQAAYZGVpZSfflK7qlWV8uab7r6pU6VFi6SFC73HI8n/BaeBgFSSFW5/+EH6+WfJFqSzL7x5/nnp0kvd5zrqKKleveJ9LltfpgQL27Vp00YTJkzI3p80aZJWrVqlu+++u/jHDwBASTVocPB9GzZILVoc2LdrYIT4P/iwwOPww8v2OewNOPPMkn/ctm3SYYeVKPi46667tG3bNm813nvuuUcPPfSQDi/r8QMAUBjL8HftavNsDw42ype3qZeKJP8HHz7SokULr+J40aJFeu+993T00UerW7du0T4sAEAIWf8Mm95qs0yiNZU1l1WrpMmTcwceOc2bJzVvrkjyf/BhQx+WgSiJJUvyz3TMmSM1a1ay1y4BKyxt2rSpXn75ZY0dO1ZvvfUWK8wCQByxzqHBhl72+91mmVixZ1QEAtKzz0q33ipt2SJVqCDt3m1Lm7tAJLiNAv9f+azmwoY+SnKrXNl9bPDCH9za/SX5PCWo98g59DJ8+HB16NBB7dq1C+VXAgAQRTHVSfTHH62jpU1zcYFHmzbSu+9KtWq5Oo/0dLe1/Ro1Yj/4+PDDD9WpUyfVqVPHq1uYPn16gc+9/vrrvecMGzZMMcW+0FF6A0477TQdcsgh+s9//hP21wIAxEYn0YhmOyZMkE45RXr7baliRWnIEJfZP/tsKTPTDbP06OG2th+FoaESD7ts377du4Bee+21uuSSSwp83rRp0zR37lwvSIk59oW2L7iloCx7kZbmUlH2JoXZ1KlTddNNN0W1sxwAIHydRHMGIBHtJLpmjbueBafStmrlCklPOunAc3Je5+z6F4HrXkiCjwsuuMC7FebHH3/UzTffrBkzZmSvYRJzIvgG2Dfihg0bvLFAi4xfffXVsL0WACC6nURtqMUyHhZ4jB49OvxFp4GAKyi95Rbp11/dH9YPPij16eNmssSg8uG40F511VW6/fbbdYqlfYqwa9cu7xa0xcam4owNVdkyxieeeKJXbFq1atVoHxIAIAysuNRq+myoxTIeYQ881q51Qyivv+72W7Z02Y5iXH/jKvh45JFHVL58ed1iEVgxDBo0SAMGDFA8s8LSvOOAAID4ZAFHRLIdzz0n3XyztGmTdMgh0gMPSHfcEbPZjrDNdlm4cKGeeOIJb5l4KzQtDuvuuXnz5uybdfwEAAAF+OknyWour7zSBR7Wo8NapN9zjy8Cj5AHHx999JHWr1+vevXqedkPu9mCaX369FGD/Fq7eqUXFb1hiJw3AACQT7bDlgWxIRWbaWrZjoEDpblzpaZN5SchDZGs1qN9+/a57rOxL7ufTp4AAJTS+vXSjTdKL7/s9q0hpk2pPfVU+VGJgw9blyTnnOWVK1dqyZIlql69upfxOMoWZsvBelrUqlVLJ5xwQmiOGACAeG69vmCBq92w/hxWQPrSSy7w2LjRDav062c1Cy7z4VMlHnZZsGCBTj/9dO9mbrvtNu//999/fziODwCAuGDtFurXr+/NfrSt7edr4kQpI0MaM8atuG43CzwsyzF/vmTXWx8HHiYpEIjgGrrFYFNtq1Wr5hWf5q3/2Llzp5dpadiwoSpVqhS1Y/Q7vo4AEPmMhwUceRuQZWZmugxIVpYLMGyyhvXSsmEW+79dom0JkJtukqwztvXwiFGFXb/jb20XAAD83nq9QQM3xGLLfVjgYYK5Afu4J5+M6cCjpAg+AADIJ1ORkZERskXhgq3XcyqXs/W6dSgtVy7/D7Y6D3s8jhB8AABQmtqMUrRet4DD5Gq9/vPP0ltvWSpE+bIF4Lp0UTwh+Ihy59PevXtH+zAAAPtZpiMtLS17iMS2tlZLKDIg1nrdajwyMjK8re3rtdekJk1ct9Jgc85ghiRPpiSexO+ZxZnZs2d7XWN/tUWDAADRqc0oI8t0tGvXTimHHSZdfbV00UXSunVu5VlbdLRWLVf3kZ7utrZfo4bijT/6sEZzvjUAIGEEazPyzkrJrs0IhTfekNLS3KJwlt24/Xa3LovNPszMdIWllgWx5+zeHbVl78MpYTMf4RjTK8z27dt19dVX6/DDD1ft2rX12GOP5Xp80qRJatmypapUqeI1Zbviiiu8VvXG0nPnnHOO9/8jjzzSy4Bcc8013v4777yjM888U0cccYTX4O0vf/mLvv/++7CeCwDEq0JrM8rKMtf2u7tTJxd4WPPNjz+WBg92gYexQCM4/GLbOAw8Ejb4COeYXkFuv/12ffDBB3r11Vf17rvvesMoixYtyn58z549GjhwoD7//HNNnz7dCziCAUbdunX18v6WusuWLdPatWu9BfyCQY01erPmb++//74Xsf/tb39jFV0ACGVtRllZQamtyWIt0ZOSpL59pcWLpTZt5JcZOyEViDGbN2+2ic3eNq/ffvst8PXXX3vbspg1a5b3GnlvGRkZgXDYunVroEKFCoEXX3wx+76ff/45ULly5UCvXr3y/Zj58+d7x2Qfa+zYbH/Tpk2FvtaGDRu853355ZcFPidUX0cAQBF+/TUQ6NbNOna42/HHBwIffxz2l3366acDycnJ3vXAtrYfzet3XgmZ+ShyvnWI2TDI7t271bp16+z7bC2cnOvdLFy4UJ06dfLWx7Ghl7PPPtu7/4cffij0c1vNSufOnXXsscd6HeWCqwcX9XEAgDCbMcPNZHn2WZftuO02ackSqW3buMvul1RCBh9hHdMrBRs6sdV/LXiYMmWK5s+fr2nTpnmPWdBSGAtYfvnlF40dO1bz5s3zbsX5OABAmGzZIv3rX9Kf/2yRgGR/2H74oWS1fpUr+37GTigk7GwXG8OzC769GZbxCGfgcdxxx3mr+1pgYJkNs2nTJn333XdehuPbb7/Vzz//rMGDB3v1HcZqOHKqsL+trn0DBdnHWA2IBR7/93//5903Z86csJ0HAKAIM2faBUZatcplO265RXr4YenQQ+Nrxk4ZJWTm46D51mHOeNgMFwt2rOh01qxZWrp0qVdMGhz6sYDEgovhw4frf//7n1577TWv+DQnm5Fjs1zeeOMNbdiwQdu2bfNmvtgMF8viWBBln9uKTwEAEbZ1q9Sjh3T++S7wOPZYa9AkDRsW0cAjFrP7+Uno4COS/vOf/3jZCRsmad++vTc9toU1kJF09NFHa/z48XrppZd08sknexmQRx99NNfHH3PMMRowYIDuuusu1axZUzfddJMXvEydOtWrF2nSpIluvfVW73UAAGFkmenUVLc1778vNW0qjRnj9m++WfriC+mss+Jrxk4IJVnVqXyyJC9LwYcGX0cAKAMbShk+XLrhBrc/apTbNmwoPfOMrZ2hRLSlkOt3Xglb8wEAQLFlZUkbN7o6jhdecPeNHu2WuzdXXSWNHGnj7FE9TL8g+AAAoCj72xjkknNGyaRJ0sSJET0kP6PmAwCAokyebJWb+T9Wvrx7HMVG5gMAgMJs3y599pn1Osj/ceuv1Lx5pI/K13wZfMRYjazv8PUDgGL66COpWzdrVX3gPmuTYEMuwS3ie9jFGnWZHTt2RPtQfC349Qt+PQEAedjvyVtvlWypCws8rAGkDa3UqiVZm4T0dLe1/Ro1on20vuOrzIc1SrGl44NLzR966KFe4y0UP+NhgYd9/ezrGGxAAwDIwZa5t2zH8uVu33pkWGv0atWkf/zDWk67WS9pabaWRdwuex9Ovgo+TC2LMqXsAAQlZ4FH8OsIANjvt9+k++6Thg51a9Aec4z09NNujZagnIGGBSAEHokRfFimo3bt2qpRo4b27NkT7cPxHRtqIeMBIJxs9VRb3MzWGImllt6F+vRT6ZprpO++c/uW+Xj8cftrLdpHFpd8F3wE2QWUiygAxJZx48ZlL+duS0DYGiOx1tr7oGzH/fe7QMOKR+vUkcaOlS68MNpHFtd81V4dABDbGQ9bBDPvaqq2tkhMZkBsiqxlO7791u1ffbVbCO7IIyPy8qv9mCEK0fXbV7NdAACxyy6kOQMPs3fvXm/V7ZwXXFvszLZRs3OndNddUtu2LvCoXVt67TVpwoSIBR7jxo3zArXU1FRva/uJhOADABAS9he8DbXkZJmPRo0axc4Fd/58N0X2kUfcMMuVV0pLl0qdOkXsEFavXp09NGVs26NHj+gGZBFG8AEACAkbOrAaj2A9nm1Hjx7t3R/1C+6uXdI990ht2khffy3VrClNn+7WZKleXbGWIYp3vi04BQDEHisu7dChg3chtYxHsJahsAtuWOodFiyQ7rhDGjLE7Vttx1dfuf9fcYX05JPSUUcpmhmifXlqY4IZokRA5gMAEFIWTLRr1y5XUFHUkEzI2QqzGRnSjTe6bIcFHtaJ9JVXpClTohZ4FJUhShTMdgEARITVeNhQi2U8ghfckE7DzcqSNm50zb/at5c2bTrw2HnnuTqP009XrFi9evVBGSI/K8n1m+ADABAfF9ziLLeR45IXb1Ndo42ptgAA3wzJhMy//13wY+XLu4Xh9ouJmTcJjMwHAMDfbKmNhx+WHnpI+v33/J+zcKHUvLk/m6H5BJkPAEBi+OILqXVr6YEHXOBxzjnu/mBxa54iV8NU1+gj+AAA+DPbMXCg1LKltHix69Xx/POuS6mt2m2NxNLT3db2baZLtGbe4CD0+QAA+It1JO3aVVq0yO1ffLE0apQLMkxmplShgitATUuTdu+WKlY8aKpr3pk3DLlEDjUfAAB/sGEVaxpmQyyW+bB1WEaMkDp3Lt5Mlzif6uqn6zeZDwBA7LMmYdal1DqXGluLZfRotyhcKVnAQdARHdR8AABiO9sxeLCbqWKBxxFHuO6lr75apsAD0UXmAwAQm775xmU7PvvM7XfsKI0ZI9WpE+0jQxmR+QAAxJa9e11th7VCt8CjWjVp/Hjp9dcJPOIEmQ8AQOz49lupWzdp7ly3f8EF0tix0jHHRPvIEEIEHwCA2Mh2DBsm3XuvtGuXVLWqVvfrp+XNm6txICDKQuMLwy4AgMiz4tHUVLf97jvprLOkvn1d4HH++Rp3332qf+edSj333Fxrr9j02IyMDG8L/yLzAQCIPJuxkpEh3XabNH++tHOnVKWK9PjjWt2hg9IaNMhugW5bawi2adMm3Xnnnd6+dSi1RmHdu3eP9pmgFGgyBgCIjKwsaeNG1xDs/POln38+8FirVtLQoVLbtl5mw1abzcsCDhaDS9CF5T788EN16tRJderUUVJSkqZPn5792J49e7yotGnTpjrssMO851x99dVas2ZN6c4EABA/GjRwa7HYeis5Aw9js1r+9KcC117JG3gYFoPzrxIHH9u3b9dpp52mp5566qDHduzYoUWLFqlfv37e9pVXXtGyZcv017/+NVTHCwDwq8ceK7gNevny0uTJudZescyGse3gwYNZDC6Raz4uuOAC75YfS7fMnDkz130jRoxQq1at9MMPP6hevXqlP1IAgD9ZxmLkSKlfP6mgkf5581wX0/2slqNDhw651l6pXr06i8HFibAXnNrYjw3PHGEtcfOxa9cu75ZzzAgAECdWrpSuvVaaPdvt25DLwoU2juKCkuC2GGuv5BeQwJ/COtV2586dXg1I586dCyw+GTRokJcxCd7q1q0bzkMCAESCBRS2zH3Tpi7wOPRQtwLtK69ItWq5ICQ93W1tv0aNYn1aCzjatWtH4JHIs10sozFt2jRdfPHFBz1mxad///vfvbnYs2fPLjD4yC/zYQEIs10AwKcyMy1NIc2a5fath8czz0jHHef27Xd+hQqu/sMuQbt3SxUrRvWQEdnZLmEZdrHA49JLL1VWVpZmzZpV6EFUrFjRuwEAfM4CCVv4zZqFbdsmVa7sVqS96SY3vBKU83e+BSBcAxJO+XAFHsuXL/fmah911FGhfgkAQKz54Qfpuuuk4KSDM8+Unn1WKsZsFMuQ2zXDptgynJIYSlzzsW3bNi1ZssS7mZUrV3r/t9ksFnj84x//0IIFCzRlyhSvInndunXebbel1QAA8ZftePppqUkTF3hYtsOahVmdRz6BR9726NY23dqnW1OxnG3UEecCJZSRkWE1IgfdunbtGli5cmW+j9nNPq44Nm/e7D3ftgCAGPbDD4FAhw4Wfrhb27aBwLJlBT796aefDiQnJ3u/4207ZMiQ7P3grVy5coFVq1ZF9DQQGiW5ftNeHQBQMnbZsCGVW2+1X9pSpUrSv/8t9eplnb/y/RDLdFhmI2eX0vy6lhrLjNiMFvhL1AtOAQBx6scfpX/9S3r7bbffpo00frx0wgmFfpjVdOQNNIILxOVdr4WupfEvrH0+AABxlO2YMEE65RQXeNgMlSFDpDlzigw8ClqvxQKNRx55JFcbdbqWJgaCDwDAwRYskGxlWdva4qCdOknXXGNtq90KtIsXS7ffXuAwS175rddigUbfvn29lWltqMW21sUU8Y+aDwDAwW65RRo+XDr/fLfi7K+/usZgDz4o9enjFoIrBav9oD16fKLmAwBQ8v4aWVnSxo2u8dfzz7snvfuu2550kvTEE9J555XpNfOu14LERPABAAnI+mmkpaVlF33akEh3axJWkG++cVmQ2EqWw6eo+QCABMx4BAMPY1tbqn71Qw+5rEd+bJhl8uTIHijiFpkPAEgw+U17tY7UK4YMUUpBmY1586TmzSNzgIh7ZD4AIMHkO+1VUiNrGHb88e6O4ON5ngeEAt9VAJBgsqe97g8sLPAYnZyslAEDpHfekWrVklq0kNLT3db2a9TId20WoDQYdgGARLNhg7rPmKEO+/ZphWU8TjxRKTa7pVkz93hmpptWa/UfaWmSLQxasWL+Rar05UAp0OcDAOJMoUvUv/yydMMNXgDiNQi79153s2CjiM+Zd20WaxRmjcGYOouSXr8ZdgGAOFLgEvXWv+Pyy6V//MMFHk2buuZhNtRSROBRaJHqCsudACXDsAsAxHKmIgRTaG14JeW++6T161224+67pX79ihV05C1SZRE4hAKZDwCIxUxFKRSYnbC6DQs8bFG4uXOlgQNLFHgUtjYLQy4oDWo+ACBKQl1Hke/ns/rRpCSl3HWX1L+/W422jMfM2izIDzUfAOADoa6j8LITQ4eq3P4upd4U2tq1lWLZjocfLnPgEXyNdu3aEXigTAg+ACCWmn2VpY7i9dfVffBgZQYCykhKUuYNN6j7//4ntWoVmgMGQoTgAwCipEx1FAsWSKmpbrtpk9S1q/TXv0pr1yrlxBPV7tNPlTJypFSpUvhPBCghaj4AIMqKW0eRa1bMkCHS8OEu4LAAZM0a1xSsb183fbZy5YieA7ClBNdvptoCQJRZwFFUtuOg7qKHHiqvt+hrr7kn1KsnDR0qXXJJzE0BBvJi2AUAYly+/Tu2bVOu1VV++EH6+99jcgowkBfBBwD4cVaM5K3Lkq18eWny5JC8XkHNylhMDqFC8AEAfpwVYwvC5bxj3jypS5eQvB6t1BFuBB8AEMu2bFHKgw9qzL59XsCR3b/DGofZTp6gJCanAAN5EHwAQKx67z23ANzYsV5xaWa3bsqYOFGZRx+t7i1bSunpUosWUq1aUo0aIXtZWqkj3JhqCwCxZutW6Y47XHBhGjaUnnlGatfO7e/a5dZmsam19it89+6QdC/Ni1bqKAmm2gKAX82aJV17rZSV5fZ79pQGD5YOP/zAc3IGGhaAhCHwKO4UYKA0CD4AIBZs2ybdeadkXUlNgwYu23HOOdE+MiDkCD4AINo++EDq1k1audLtX3+9ZB1Mq1SJ9pEBYUHBKQBEy/bt0i23uFoOCzysS+nMmdKoUQQeiGtkPgAgGj76yGU7vv/e7aelSf/5j0ShPRIAmQ8AiKQdO6TevbX6rLOU8f33Wl27tjRjhjR6NIEHEgbBBwBEypw50mmnadwTT6i+pFRJ9X/6SeNWrYr2kQERRfABAOFiS92nprqgo08f6ayztHrFCqXZein7n8K6KUhEBB8AEC4TJ0oZGVKnTtLjj3sNwZb/+c/ZgUcQ66Yg0VBwCgChZM3BNm50XUiDy9D/+qv0hz9I992nxq1bK/ndd3Mt3Ma6KUg0BB8AEErWHCw/FpD07u0tBjfm6ae9oRbLeLBuChIRwQcAhMrOnVLHjtKbb+b/ePny0vjx6t6lizp06MC6KUhYBB8AEArz50tdu0rffFPgU1a/9pqWV6qkxqtXs24KEhoFpwBQFlbbcc89Ups2LvCoWVN67DH3WPL+X7HJybLqj/p/+YtSU1NVv359jQvWgwAJiOADAMoylbZFC2nQIJszK11xhfTVV9Kll0q1arnH0tO1umlTN712f5Ep02uR6Bh2AYDSZDsGDnRL3e/dK9Wo4QUZ+tvfDjwnM1OqUMFb8n5548bad+65+U6vZegFiYjgAwCKwbIUy5cvV+OdO5Vyxx3S0qXugcsvl4YPd1Npc6pYMfu/jY8/XsnJyUyvBfZj2AUAimD1GVan4dVrXHihxlngcfTR0n//Kz3//MGBRx6W3RgzZowXcBim1yLRJQUCgYBiyJYtW1StWjVt3rxZVVlkCUAMZDws8MiVtUhKUubixUo57bQSfy6m1yJeleT6zbALABRkzx4tv+++XIGH2RsIaMWmTV7DsJJgei1QymGXDz/8UJ06dVKdOnWUlJSk6dOn53rcEin333+/ateurcqVK6t9+/beOCkAlIVlDTIyMiI3Q+Tzz6VWrdR4woSDflFSrwFEOPjYvn27TjvtND311FP5Pj5kyBA9+eSTSk9P17x583TYYYd5nfx2Wuc/AChrzUW4e2Ts2eNmspxxhrRkiVKqV9eYtDTqNYBYqfmwzMe0adN08cUXe/v2qSwj0qdPH/Xt29e7z8Z+atasqfHjx+tyqwovAjUfAIqsuShXTpmZmWUPAKxPh81cGTJEatnSm8GyunNnLV+6VI1tmMR+t40a5fXsyK9eI3sGTOPGBCNIeFtKcP0O6WyXlStXat26dd5QS5AdSOvWrfXpp5+G8qUAJAi7uB9UcxGqJeiDS95PmCA9/LDGNWum+kuXKtW6kSYlaZyt02LNwvbXa7Rr1y47yIhoNgaIMyEtOLXAw1imIyfbDz6W165du7xbzsgJAIIsqxDSHhnBJe+TkqQXXnD3WRfS3393XUj3P21fIKAe11+vDn/+80FZDct4pKWlHdSx1IaYyYAAPujzMWjQIC87ErzVrVs32ocEIIaEvEeGLXlvQyzW+nz9enff77/LyuJz51cKzrCENRsDJICQBh+19qcnf/rpp1z3237wsbzuvvtub3woeFu1alUoDwlAHOjevbtX42GzXWxr+6U2ebJFMAfdbTUexZ3VEszGFOe5AMIcfDRs2NALMt5///1cwyg26+WPf/xjvh9TsWJFrzAl5w0A8spbc1Eqtg7LmjUHVpvN+fkljenXr1gZlmA2JhiA2JYZMEAYaz62bduWK7VoRaZLlixR9erVVa9ePfXu3VsPPfSQ95eBBSP9+vXzZsAEZ8QAQFQsWyZdc400d+6B+6zuwyb8WRCxb5+6X3yxOqSl0YUUiLWptrNnz9Y555xz0P1du3b1ptPap+vfv7/3V8Gvv/6qM888UyNHjtTxxx9frM/PVFsAJVHkdFfLdgwbJt13n2T9huz3Sr9+0mOPSVZjZkM4NlPFhnznz7e0RnSn/wI+VZLrN2u7APAtm94anHViQx/2R0+uepDvvpO6dZM++cTtd+ggjR3rgg6bZbd/yXsv+7F7d66VaAtjtSc2xTa/+21oCEhEW6LV5wMAIqWg6a5e+3XLdgwdKtnCbxZ4VKnigo6333aBh7FAwwIPY9tiBh6GglOgbAg+APhSgdNdP/hAsuzDbbe5YRZrerh0qXTddQeCjVib/gskGIZdAMTXUvcVKyrFgo7DD3d1Hf/6V8iCjvyOgeJUwGHYBUDcOyj7IGl0IOACD6vH+PJLKS0tbIFHyKb/AgmI4AOAb3Xv1k2Z/fsro2JFZdr+YYdJI0dKM2e6TqYA4n9tFwCImJUrpWuvVcrs2V6DMFkLAJsy27BhtI8MQBHIfADwF6vxsGXumza1xkPSoYdKI0ZI771H4AH4BMEHgBIVWFovC286ayQsWODqN2xrMjOl886TbrxR2r5dOuss6YsvpJ49822ZDiA28dMKoNgNvWx2iTXXsq3thz2ImTjROne57ejRLtsxa5ZUubL0xBPuseOOK/VxAIgOptoCiGg78SK7kmZlSRs3ulkqF1zglr0/5BBpzx73eMuW0vPPSzT0AmIK7dUBhFSo2okXK4gpztTY2Pq1BUD0+QAQYqFqJ15gV9IcK2Vr8mSpfAET8ex+exyArxF8AIhYO/EigxjLaNiCbwWtszJvntSlSynPAkCsIPgAUCxWl2HDIzbUYttcdRqhCGJ+/FHq2NEtcW8zWXIOwTCTBYgr1HwAiLhca6Icc4w0YYLUu7e0ebPLevTp4xqG1avnghH7/6pV0vz5FsGU+bVt+MeyMLRFB0KHglMA/rBmjVt/5c033X6rVtL48dJJJ7nhlwoVXPbDfk3t3l2iZe9LNdMGQKkRfACIbfZrxwpHb7lF+vVXF2Q8+KDLeBRUbBpD04UBlO36zdouACJr7VqpRw/p9dcP9O2wbMcpp4T1ZQubaUPwAUQWVVwAIpftmDLFBRkWeFjjsH//W/r007AHHqGcLgyg7Ag+AITfunXSJZdIV14pbdokNW8uLVwo3XNP2IZZwjVdGEDZUfMBIHzs18vUqdJNN0m//OKyHfffL915p/t/tGfaEHgAIUPNB4DoszVZbrhBeuUVt9+smZtSe+qpUT0sCzgIOoDoYtgFQOi9+KKr47DAw4ZVBgyQPvss6oEHgNhA8AGg9BYskGzBOduaDRukSy+VLrvMrUxrwYY1BrOhligNswCIPQy7ACi9iRNtaVtp0iQpM1O68UYXgFhR5733upv18ACAHCg4BVAyWVkuq2GdRy+4wNV2WOdR60hqTjhBeu45N6MFQMLYQsEpgLBp0ODg+4KBh1m2jMADQKGo+QBQMtYWvaDeHHa/PV6M6a62Oq5tASQegg8AJXP44VJBKdV586QuXYpc3M3WWElNTfW2tg8gsRB8ACgeaxJ21VXSxRe7/xur+zB52pYXxDIdwVVljW179OhBBgRIMAQfAIr2xhtSkyZuSMUCDWseVrOmWxQuPV1q0UKqVUuqUaPUi7sBSBwUnAIomK3D0ru3m1IbnMliK9C2aSMNHeqm0Vr2Iy1N2r3bzXopxuJueZe1Z3E3ILGQ+QCQv7fectkOCzwswOjbV1q82AUexgKN4LCLbYsIPAyLuwEw9PkAkNuvv0q33SY9+6zbP/549/+2bUP2EizuBsQf+nwAKJ133pGuu0768UeXzbj1Vumhh6TKlUP6MizuBiQ2gg8A0ubNUp8+Ng/W7VsNhmU7zjwz2kcGIA5R8wEkupkzpaZNXeBh2Y5evaTPPyfwABA2ZD6ARLV1qysiHTPG7R97rMt2nHVWtI8MQJwj8wHEWNvviBzD+++7bEcw8Lj5ZumLL8oceMTC1w9A7CP4AGKo7XdYjmHBAik11W0t22ENwtq3d6vTNmwoZWRITz4pHXZY7B07gLjEVFtg/1/sdsHM2/wqMzMzYrMywnYMt9wiDR8uXXKJtGiRlJnp7u/ZUxo82K3VEqvHDsA3SnL9JvMBxEjb75Aeg2U1Fi50wcbUqe6+V15xgYe1QX/uOWnEiJAEHiE/dgBxj4JTIEbafof0GBo0KPixdeukK66QOndWPH39APgHmQ8gRtp+h/QYrN6ioJVmy5d3C8TF2dcPgH9Q8wHEWNvvMh/DRx9J3bpJ33+f/+M2HNO8ueL16wcgOmivDvi47Xepj2HHDunee6UnnpDsbwpb8v6nn1wGxIZDgts4//oBiH0MuwDx4OOPpWbNpGHDXODRvbs0a5YrLm3RQkpPd1vbr1Ej2kcLIMGR+QD87LffpPvuk4YOdUHHMcdIY8dKF1zgHrfZLRUquLbpaWnS7t1SxYrRPmoACS7kmQ+bXtevXz81bNhQlStX1nHHHaeBAwcqxkpLAP/79FOX7Xj8cRd4WJ3H0qUHAg9jgYYFHsa2BB4A4jHz8cgjj2jUqFGaMGGCTjnlFC1YsEDdunXzilBusWZHAMqe7bj/fhd0WA1HnTou23HhhdE+MgCITvDxySef6KKLLlLHjh29/QYNGuj555/XZ599FuqXAhLPvHnSNddI337r9q++2tV5HHlktI8MAKI37NK2bVu9//77+u6777z9zz//XHPmzNEFOVPBOezatcubnpPzBiCPnTulu+6yHzAXeNSuLb32mjRhAoEHAN8Jeebjrrvu8gKIE0880Ws0ZDUg//73v9WlS5d8nz9o0CANGDAg1IcBxI/581224+uv3f6VV7rptNWrR/vIACA2Mh8vvviipkyZoueee06LFi3yaj8effRRb5ufu+++22tIErytWrUq1IcE+NOuXdI990ht2rjAw/p2TJ8uTZpE4AHA10Le4bRu3bpe9qOnrZi530MPPaTJkyfr2+A4dSHocApIWrDAZTu++srt21ostuz9UUdF+8gAIPZWtd2xY4e3wFRONvySd8VLADkCjdRUt7U+HP36uWyHBR7WEMxWo50yhcADQNwIec1Hp06dvBqPevXqeVNtFy9erMcff1zXXnttqF8KiA8TJ0oZGdJjj7mA48sv3f2XXy4NHy794Q/RPkIAiO1hl61bt3pNxqZNm6b169erTp066ty5s+6//35VsE6LRWDYBQkhK0vauNE1/vrzn6UNGw48dsQRVoktXX99NI8QAEqkJNdvVrUFoiHYdbQwsfWjCQCxW/MB+IktAZ+RkeFtI2r8eLfKbH7Kl5cmT47s8QBABBF8IGGNGzdO9evXV2pqqre1/Yj44gvXp6OgImzrYlpAXxwAiAcEH0hIlulIS0vLnoVl2x49eoQ3A7JnjzRwoNSypbR4sRRMSwYzIAVlQgAgzvDbDglp+fLlB03/tm68K1asCM8L2mqzNn3WFoSzIOTii90Ml1q1pBYtpPR0t7V9m14LAHEs5FNtAT9o3Lix148mZwBi/WgaNWoU2hf6/XdpyBDpgQdc0GHrsIwYIXXu7IpOMzMlmwVm/09Lc30+WPYeQJwj84GElJKSojFjxngBh7Ht6NGjvftDxlqi//GP0r33usCjUyfXx8O6lQZnu1igEfy/bQk8ACQAptoioVmNhw21WMYjZIGHZTsefVTq399lMqxvh7VGtwXhiphia8djQ0KWmQlpIAQAMXT9ZtgFCc0u8CG9yH/zjVuT5bPP3H7HjtKYMVKdOkV+qM22CRbB2pCQZWa6d+8eumMDgBhB5gMIhb17pccfd+uy2Gq01apJw4ZJXbsWq6GYZTxsum/eGpTMzEwyIAB8gcwHEEnLlrlsx9y5bv+CC6SxY6VjjgnJ7BuCDwDxhoJToCzZDlsMrlkzF3hYpG+Nyt58U6sDgRJ1Tg3OvskpLLNvACAGEHwApfHdd9JZZ0l9+0o7d0rnn+96eVx7rcY980yJO6dGZPYNAMQIaj6AoixYIN1xh+vXcfrpbubKPfe4oKNKFVfrYYWhSUllrt0Iy+wbAIgAaj6AUJo40XUjHT5c+t//pDlz3P3t27thlnr1Qla7EfLZNwAQgwg+gPxkZUkbN7qZKi+8cCAIMZUrS/fdJ91990EzWSLWORUAfIyaDyA/DRq4BeBsvZX163M/9ttvrmtpPlNoqd0AgKIRfAD5sSxHQavMli8vTZ5c4IdaYzCr8bDZLralURgA5MawC5DXypXSM89IeWo3ss2bJzVvXuinoHYDAApG5gMIsmBj1CipaVNp9mypUiV3fzADUlAmBABQIvw2BYwtbX/eedKNN0rbt0v/93/SzJlSrVqu7iM93W1tv0aNaB8tAPgawy5IbNbmxhZ+s2Zh27a5mSyDB0s33eQyHRaUVKjgikvT0twqtSx7DwBlQvCBxJ5Oe9110nvvuf0zz5SefVbKOS02Z6BhAQiBBwCUGcMuSMxshy38ZrUdFnhYbcfQoa7Og34cABB2ZD6QWFatkv71L2nGDLfftq3Ldhx/fLSPDAASBpkPJE62w6bPNmniAg8bPnn0UenDDwk8ACDCyHwg/tmy9lYs+vbbbr9NG5ftOPHEaB8ZACQkMh/wBVvt1TqG2rZE2Y7x4122wwIPy3bYyrS2MByBBwBEDcEHYt64ceO8ZepTU1O9re0Xac0aqVMnqVs3afNmqVUrafFi6fbbbcGVSBw2AKAABB+IaZbpSEtLy14l1rY9evQ4OAOyYIGUmirNny9NmiSdcor05puuR8egQdLHH0snnRSdkwAA5ELNB2La8uXLcy1Pb/bu3asVK1bkXjvFFoLLyJAuu8ytzWKsI+mECS4QAQDEDIIPRJVlMCzAaNy4cb4Lsdn9ycnJuQIQW6a+kfXjsCZhGze6Oy3IMBZ42KqzVmB6223SccdF7FwAAMXDsAtiupbDApIxY8Z4AYex7ejRo12g0qCB1LKlu23ZcuCDfv9dGjmShmEAEKOSAgGbEhA7tmzZomrVqmnz5s2qWrVqtA8HYcx4WMCRN6ORmZmZbwbEnm9DLZbx8B63b1tbf8WCjPxY9sNmunTpEs7TAACU4vrNsAtiu5ZjP7sv+/7166UbbpBeeaXgF5g3T2rePOTHDQAoO4ZdEBXBWo6csms5CvPii66A1AKPYG2HCX6uPJ8TABB7+E2NqCi0liM/GzZI//ynm81iRaannuqm1fbrJ9Wq5Wa2pKe7re3XqBHZEwIAFBs1H4iqg2o58vPf/0o33ugCEAtW7r3X3ayHh9m1y/3flry3b+fdu103UwBAxFDzAd/IVcuRl2U4rKj0hRfcftOmrog0by1HzkDDAhACDwCIaQy7IDZNm+ZqOyzwCGY7bJglxEWkpVozBgBQJgQfiC0//yxdcYV0ySVuVosFIHPnSg89FPKMRqnWjAEAlBk1H4gdr74q9egh/fSTm7Vy551S//5hGUYpaZ8RAEDort9kPhB9v/wiXXmldPHFLvCwBeA+/VR6+OGw1W8U1mcEABBeBB+Irtdfd0MrU6YcyHYsWiS1ahWbfUYAAGVG8IHo2LRJ6tpV+utfpXXrpBNPdMveDx4sVaoUe31GAAAhQ80HImPBAumOO6QhQ9zQinUmXbPGTY3t21caMECqXDk2+4wAAIpEnw/EnokTpYwM6eqrpW++cfcdf7z07LNS27ax2WcEABAWBB8In6ws1yjMshsWfJhg4GGrzVpr9BNOiOohAgAij2EXhI8FHUWJrW8/AIBfp9r++OOPuvLKK3XUUUepcuXKatq0qRbYmD8Si81cKYitSDt5ciSPBgAQr8MumzZt0p/+9Cedc845evvtt3X00Ud7PRWOPPLIUL8UYtWWLa6IdOzYgp8zb17IW6UDABI0+HjkkUdUt25dPWuFhPs1bNgw1C+DWPXee1L37tIPP7j9yy5z67NYTw1r6hXcAgASVsiHXV577TW1bNlS//znP1WjRg2dfvrpGlvIX8C7du3yxoly3uBDW7dK118vnXeeCzws4LTZLY8+KtWqJbVoIaWnu63t16gR7SMGAMRLwWml/Q2ibrvtNi8AmT9/vnr16qX09HR1taZSeTzwwAMaYD0e8qDg1EdmzZKuvdbNbjE9e7pmYYcf7vZ37ZIqVHAFqPbttns3y94DQAIXnIY8+KhQoYKX+fjkk0+y77vlllu8IORTW68jn8yH3XIevA3bEHz4wLZtrqh05Ei336CBLRUrpaZG+8gAAIk026V27do6+eSTc9130kkn6YdgDUAeFStW9A4y5w0+8MEH0qmnHgg8bMjliy8IPAAARQp58GEzXZYtW5brvu+++85bvhxxYPt2S2VJ7dpJK1dK9epJM2dKo0ZJVapE++gAAIkYfNx6662aO3euHn74YW/NjOeee85bwKun1QHA3z76SDrtNGn4cLdv67N8+aXUvn20jwwAkMjBxxlnnKFp06bp+eefV5MmTTRw4EANGzZMXaydNvxpxw6pd2/p7LOl77+X6taVZsyQRo+WGCYDAJQQ7dVRuDlzpG7dpBUr3L718HjsMalatWgfGQAghkS9vTriwG+/SX36SGed5QKPY46R3npLevppAg8AQJkQfMCxtXdspoptbUp0s2bS44+7vhyW+Vi6VLrggmgfJQAgDoS8vTp8ypa8t46kPXpIS5a4Fuh16rj1WS68MNpHBwCIIwQficw6km7c6DqPBleYXbTIbTt2lB5+2PXyAAAghCg4TWQWdBQltr49AAAxioJTFM+DDxb8WPnyB7IhAACEEMMuicjW0rHF/B55pODnzJsnNW8eyaMCACQIMh+Jxmaz2LL2gwa5otIOHdz9ycm5twAAhAlXmkRhy9j36ye1aSN99ZVUo4b0yiuub0etWi4gSU93W9u3xwEACAOGXRKBzWC55hq3Dou57DJpxAjpD39w+5mZUoUKrgDV1muxQKVixageMgAgfpH5iGcWRPTvL7Vq5QIPCzZeekmaOvVA4GEs0AjOfLEtgQcAIIzIfMSB1atXa/ny5WrcuLFSUlLcndYozLIdn3/u9v/5T+mpp6Sjj47qsQIAQPDhc+PGjVNaWpr27dun5ORkjRk1St3XrpUeekj6/XfpqKOkkSOlSy+N9qECAOChyZjPMx7169f3Ao+gclbCIcnLf1xyiQs8ataM5mECABLAFpqMJQYbaskZeJi9klZUqSI9/7z03/8SeAAAYg7DLpGswwgx+9w21JIr85GUpEazZ9MgDAAQs8h8hKkOw4ZDUlNTva3th9zvvytl4kSNSUryhlpMueRkjR4zRikEHgCAGEbNRyTqMMqVU2ZmZtkyINaZ9I47pCFDpMqV3UwWu89es317rbj+ejVq3TpsWRYAAEJ1/WbYJRJ1GHv3asWKFWULDCZOlDIypF69XNBhPTyOOEJ68kmlXHmlUoqzQi0AADGA4CMSdRjlyqlRo0Yl/2RZWdLGja7x15Qp7r5PPnHbM8+UHnvMNRADAMBHqPkIMctujBkzxgs4jG1Hjx5duqxHgwZSy5ZuvZVffsn92Jw5UuvWITpqAAAih+AjDLp37+7VeGRkZHhb2y8Vq+8oaDilfHlp8uQyHScAANHAsEuYWKaj1DUee/dKw4a5VWgLqgeeN4/ptAAAXyLzEWu++0466yypb19p1y6pTRt3f3Jy7i0AAD7FlSxWWLZj6FDptNNcUal1KR07VnrxRalWLVf3kZ7utrZfo0a0jxgAgFJh2CUWrFghdevmikhN+/bWqUyqV8/tZ2ZKFSq4+o+0NDfNlmXvAQA+ReYjmmw67hNPSKee6gKPww+XRo+W3n33QOBhLNAIFp7adn/gYQ3NrKjVtgAA+AXBRwHCfmH//nvpnHOk3r2l336TUlOlL790mY1iNAyLSAt3AADCgOAj0hd2y3aMGOGyHR9+KB12mFv2fuZM19ejGCwgSktLy25kZtsePXqQAQEA+ALBRyQv7CtXSueeK918s7Rjh8t8WLbjhhtKNIulsBbuAADEOoKPSFzY7fONGiU1bSrZcveHHuqyH++9JzVsWOoW7jmVuoU7AAARRvAR7gu7zVQ57zzpxhul7dtdD48vvpB69ix1z46QtnAHACDCCD5CeWG31WatcNS21pnUZq5YtmPWLKlyZTezxVamPe642GnhDgBAhCUFAgX1746OLVu2qFq1atq8ebOqVq0ateOwGg8barGMR7EzCrfcIg0fLl1zjX0CN6wSXIH22WclhkUAAHGqJNdvmoyVdW2WnMvev/CCu2/CBJf5sH4cd9wh9e9vKZSwHzMAAH5A5qOsitGTo8DF4QAAiBMluX5T81HWxmKTJhVcOLp/2Xs6kQIAcADBR1kai1kw8dxzbiptfubN07idO+lECgBADgy77M94WGCQs7+HzXKxWST51n3Yl8zqOqw1+ubN0iGHSHv2uAyIfY7929VvvaX6f/lL8T8vAAA+xbBLOBuLrVkjderkVqG1wKNVK2nGjHyXvV++dSudSAEAyIPZLjkai+XNUORqLGbZjsmT3XTaX391S9w/+KDUp4+r7chn2fvGGzYU/XkBAEgwZD6K01hs7Vrpooukq692gUfLltKiRdKdd7rAo4Bl7+lECgDAwaj5KKyxmH1prKDUFoLbtMnVdjzwgOvdEQw6SvN5AQCIMzQZC0VjsZ9+kq6/Xpo+3e1bHcf48VKTJmX7vAAAJLiEHHYptO+GZTuef146+WQXeFi2Y+BA6dNPSxV4AACABA8+Cu3nsX699I9/SFdcIf3yi9SsmVsk7r77XBACAADKLKFqPgrt5/HJJ26Ze1unxeo5+vWT7r67REGHfX6btmuzZxhmAQAkki30+ShhP48rr5Quu8wFHqeeKs2fL91/f4kCj1J1SAUAIAGFPfgYPHiwkpKS1Nu6gUZZ4x07DjphmwTb6IMPXLbDAg4LPGy4pQQs45GWlpYd2Ni2R48erOUCAECkg4/58+d7fS1OtWxCDEiZMUNjLODY34/DAo/Rdn/Tpt46LBowwDUKC2eHVAAAElzYgo9t27apS5cuGjt2rI488khFTVaWtHChawr2wgvqLikzEFCGbZOS1P2mm1xRafPmZe6QmhOdTAEAiHDw0bNnT3Xs2FHt27cv9Hm7du3yilRy3kKqQQPXkdT6dNhsFst0SGpnW6u1HTGiVNmOnOhkCgBA8YWlydjUqVO1aNEib9ilKIMGDdIAG+4IF1uP5ZprpN9/P/gxq/OwxmEh0L17d3Xo0IFOpgAARHqq7apVq9SyZUvNnDkzu9ajXbt2atasmYYNG5Zv5sNuQZb5qFu3bmin2tqQi2U+8rLhmDzDLUyXBQDAZ1NtFy5cqPXr16t58+YqX768d/vggw/05JNPev+3QsycKlas6B1kzlvYBOsy8tRnBDFdFgAAH2Y+tm7dqiwr8syhW7duOvHEE3XnnXeqSREtysPSZMymvJ5xhlS3ro2PWJRhKRo3rXZ/dqPQBmRkQAAAiN2F5apUqXJQgHHYYYfpqKOOKjLwCBsLHjIzXWGpTbNNS5N27/aWvS/OdFmCDwAAQidxOpxaoLG/v4e3zRF4GKbLAgDg49kuec2ePVuxLjhd1jqTWsaD6bIAAIRHQi0sVxxW+8F0WQAAfFTz4XcWcBB0AAAQPolT8wEAAGICwQcAAIgogg8AABBRBB8AACCiCD4AAEBEEXwAAICIIvgAAAARRfABAAAiiuADAABEFMEHAACIKIIPAAAQUQQfAAAgogg+AABARBF8AACAiCL4AAAAEUXwAQAAIorgI0RWr16tjIwMbwsAAApG8BEC48aNU/369ZWamuptbR8AAOQvKRAIBBRDtmzZomrVqmnz5s2qWrWqYp1lOizg2LdvX/Z95cqVU2ZmplJSUqJ6bAAAxOL1m8xHGS1fvjxX4GH27t2rFStWRO2YAACIZQQfZdS4cWMlJ+f+Mlrmo1GjRlE7JgAAYhnBRxnZ0MqYMWO8gMPYdvTo0Qy5AABQAGo+Qlj7YUMtlvEg8AAAJJotJbh+l4/YUcU5CzgIOgAAKBrDLgAAIKIIPgAAQEQRfAAAgIgi+AAAABFF8AEAACKK4AMAAEQUwQcAAIgogg8AABBRBB8AACCiCD4AAEBEEXwAAICIirm1XYLr3NkCNQAAwB+C1+3irFcbc8HH1q1bvW3dunWjfSgAAKAU13Fb3bYwSYHihCgRtG/fPq1Zs0ZVqlRRUlJSyKMyC2pWrVpV5HK/fhTv55cI58j5+V+8nyPn539bwnSOFk5Y4FGnTh0lJyf7K/NhBxzupentix2v31SJcH6JcI6cn//F+zlyfv5XNQznWFTGI4iCUwAAEFEEHwAAIKISKvioWLGi+vfv723jUbyfXyKcI+fnf/F+jpyf/1WMgXOMuYJTAAAQ3xIq8wEAAKKP4AMAAEQUwQcAAIgogg8AABBRcRd8PPXUU2rQoIEqVaqk1q1b67PPPiv0+S+99JJOPPFE7/lNmzbVW2+9pXg5v/Hjx3tdYnPe7ONi1YcffqhOnTp53fHsWKdPn17kx8yePVvNmzf3qrYbNWrknXMsK+k52vnlfQ/ttm7dOsWaQYMG6YwzzvC6E9eoUUMXX3yxli1bVuTH+elnsDTn6Kefw1GjRunUU0/Nbj71xz/+UW+//XbcvH+lOUc/vX/5GTx4sHfMvXv3Viy9j3EVfLzwwgu67bbbvClEixYt0mmnnaYOHTpo/fr1+T7/k08+UefOndW9e3ctXrzY+0Vit6VLlyoezs/YD9fatWuzb1lZWYpV27dv987JAqziWLlypTp27KhzzjlHS5Ys8X64rrvuOs2YMUPxco5BdoHL+T7ahS/WfPDBB+rZs6fmzp2rmTNnas+ePTr//PO9cy6I334GS3OOfvo5tO7SdrFauHChFixYoNTUVF100UX66quv4uL9K805+un9y2v+/PkaPXq0F2wVJirvYyCOtGrVKtCzZ8/s/b179wbq1KkTGDRoUL7Pv/TSSwMdO3bMdV/r1q0DPXr0CMTD+T377LOBatWqBfzIvjWnTZtW6HPuuOOOwCmnnJLrvssuuyzQoUOHQLycY0ZGhve8TZs2Bfxm/fr13rF/8MEHBT7Hbz+DpTlHP/8cmiOPPDLw9NNPx+X7V5xz9Ov7t3Xr1kDjxo0DM2fODJx99tmBXr16FfjcaLyPcZP52L17txfJtm/fPtc6Mbb/6aef5vsxdn/O5xvLJBT0fL+dn9m2bZvq16/vLSJUVHTvN356/8qqWbNmql27ts477zx9/PHH8oPNmzd72+rVq8fte1icc/Trz+HevXs1depUL6tjQxPx+P4V5xz9+v717NnTywznfX9i5X2Mm+Bj48aN3jdSzZo1c91v+wWNj9v9JXm+387vhBNO0DPPPKNXX31VkydP9lYMbtu2rVavXq14UND7Zys2/vbbb4oHFnCkp6fr5Zdf9m72y69du3besFsss+81Gwb705/+pCZNmhT4PD/9DJb2HP32c/jll1/q8MMP9+qorr/+ek2bNk0nn3xyXL1/JTlHv71/xgIq+x1hNUrFEY33MeZWtUXoWCSfM5q3H5iTTjrJGwMcOHBgVI8NxWO/+OyW8z38/vvvNXToUE2aNEmx/FeXjRfPmTNH8aq45+i3n0P7frMaKsvq/Pe//1XXrl29WpeCLs5+VJJz9Nv7t2rVKvXq1curSYrlwti4CT7+8Ic/qFy5cvrpp59y3W/7tWrVyvdj7P6SPN9v55fXIYccotNPP10rVqxQPCjo/bPisMqVKytetWrVKqYv6jfddJPeeOMNb2aPFfcVxk8/g6U9R7/9HFaoUMGbOWZatGjhFS0+8cQT3sU2Xt6/kpyj396/hQsXepMQbBZgkGXN7Xt1xIgR2rVrl3ctifb7GDfDLvbNZN9E77//fvZ9lh6z/YLG8uz+nM83Fi0WNvbnp/PLy74BLd1oqfx44Kf3L5TsL7ZYfA+thtYuypbCnjVrlho2bBh372FpztHvP4f2e8YuWPHw/pXmHP32/p177rne8dnvieCtZcuW6tKli/f/vIFH1N7HQByZOnVqoGLFioHx48cHvv7660BaWlrgiCOOCKxbt857/Kqrrgrcdddd2c//+OOPA+XLlw88+uijgW+++SbQv3//wCGHHBL48ssvA/FwfgMGDAjMmDEj8P333wcWLlwYuPzyywOVKlUKfPXVV4FYrc5evHixd7Nvzccff9z7f1ZWlve4nZudY9D//ve/wKGHHhq4/fbbvffvqaeeCpQrVy7wzjvvBGJVSc9x6NChgenTpweWL1/ufV9axXpycnLgvffeC8SaG264wZsVMHv27MDatWuzbzt27Mh+jt9/Bktzjn76ObTjtpk7K1euDHzxxRfeflJSUuDdd9+Ni/evNOfop/evIHlnu8TC+xhXwYcZPnx4oF69eoEKFSp4U1Pnzp2b6w3o2rVrrue/+OKLgeOPP957vk3bfPPNNwPxcn69e/fOfm7NmjUDF154YWDRokWBWBWcVpr3Fjwn29o55v2YZs2aeed47LHHetPiYllJz/GRRx4JHHfccd4vu+rVqwfatWsXmDVrViAW5Xdedsv5nvj9Z7A05+inn8Nrr702UL9+fe9Yjz766MC5556bfVGOh/evNOfop/evuMFHLLyPSfZP+PIqAAAAcVrzAQAA/IHgAwAARBTBBwAAiCiCDwAAEFEEHwAAIKIIPgAAQEQRfAAAgIgi+AAAABFF8AEAACKK4AMAAEQUwQcAAIgogg8AAKBI+n9KnOerYRVnPgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "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": 12,
   "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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "x,y = DataGen(9, 100, 1, 0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAja0lEQVR4nO3dCYxeVfk/8DNt6SIphCVsmbb0B43s+xKKwVqJI0ECxiAYSBAbIAiyNC6AAiEsBSJNBRFaUgEDFUgIghgwpDQgUmQTIkGgFUtaCFuCtIAWbeefc/1PnQ4z03ln7nvvufd+Psmbl/ftlPdOh3K+9znPOaeju7u7OwAAFGRUUR8EABAJHwBAoYQPAKBQwgcAUCjhAwAolPABABRK+AAACiV8AACFGhMSs379+vDWW2+FiRMnho6OjrIvBwAYgrhn6Zo1a8JOO+0URo0aVa3wEYPHpEmTyr4MAGAYVq5cGTo7O6sVPmLFo+fit9hii7IvBwAYgtWrV2fFg55xvFLho2eqJQYP4QMAqmUoLRMaTgGAQgkfAEChhA8AoFDCBwBQKOEDACiU8AEAFEr4AAAKJXwAAIUSPgCAQgkfAEChhA8AqKFVq1aFJUuWZM+pET4AoGYWLlwYpkyZEmbOnJk9x9cp6eju7u4OiZ2Kt+WWW4YPP/zQwXIA0KJY6YiBY/369RveGz16dFixYsUmj7ovavxW+QCAGlm2bNlGwSNat25dWL58eUiF8AEANTJt2rQwatTGw3usfOy6664hFcIHANRIZ2dnWLBgQRY4ovg8f/78tk65tErPBwDUtPdj+fLlWcWjiODRyvg9pu1XAwAUrrOzM6lqR2+mXQCggntlVJnwAQAV3CujUeHj8ccfD8ccc0zYaaedQkdHR/jNb36z0a/HFpJLLrkk7LjjjmHChAnhyCOPzJb9AEBVxErH6aefvmHJanw+44wzVEDKCh8ff/xx2HfffcONN97Y769fe+214frrrw8333xz+NOf/hQ233zz0NXVFf71r3/lcb0A0HZV2CujylpuOD3qqKOyR39i1WPevHnhJz/5STj22GOz9371q1+F7bffPquQnHjiiSO/YgAoaK+MvruEprRXRpXl2vPx97//Pbz99tvZVEuPuOzm0EMPDUuXLu3396xduzZbntP7AQBlqsJeGVWW61LbGDyiWOnoLb7u+bW+5syZEy677LI8LwMARmzWrFlZ20CRe2U0RemrXS688MJsQ5Kex8qVK8u+JADIxMAxY8YMwSPl8LHDDjtkz++8885G78fXPb/W17hx47Kd0Ho/AID6yjV8TJ06NQsZixcv3vBe7OGIq14OO+ywPD8KAGhKz8dHH3200VKj2GT6wgsvhK233jpMnjw5nHfeeeGKK67IOoVjGLn44ouzPUGOO+64vK8dAGhC+Hj22WfDl770pQ2vZ8+enT2fcsop4bbbbgs//OEPs71A4uYs//jHP8IXvvCF8PDDD4fx48fne+UAQCU51RYAKHT8Ln21CwDQLMIHAFAo4QMAKJTwAQAUSvgAAAolfAAAhRI+AIBCCR8AQKGEDwCgUMIHAFAo4QMAKJTwAQAUSvgAAAolfAAAhRI+ACBhq1atCkuWLMme60L4AIBELVy4MEyZMiXMnDkze46v66Cju7u7OyRk9erVYcsttwwffvhh2GKLLcq+HAAoxapVq7LAsX79+g3vjR49OqxYsSJ0dnaG1LQyfqt8AECCli1btlHwiNatWxeWL18eqk74AIAETZs2LYwatfEwHSsfu+66a6g64QMAEtTZ2RkWLFiQBY4oPs+fPz/JKZdW6fkAgMR7P5YvX55VPFIOHq2M32MKuyoAoGUxcKQcOobDtAsANGR/jVQIHwDQkP01UqHnAwAqvL9GKuzzAQDDnE6p8/4aqRA+AGiMoUynFLm/xqqG9pUIHwA0QhzgTz/99A1Vjfh8xhlnfGbgL2p/jYUN7ivR8wFAI8QKQxzo+3t/xowZhe6vsaqGfSX2+QCAAaZT+g74A02ntHN/jWWD9JVUNXy0wrQLAI2Q0nbl02p8bstQCB8ANMasWbOyqY041RKf4+t22FQjaWdCQagMej4AaIQYBOJ0R6w6tHOQj42jPY2tsboRQ8ZAIWdVRc5tyXv8Fj4AqL1WAsFI1LGRdKhsMgYALS6xzYMNyoZG+ACg1ooMBE1vJB0q4QOAWisyEKTUSLoq4d1ThQ8Aaq3oQFDUipoq756q4RSARqjTypIUm17tcAoABe5YmpJlFdg91bQLANTItAo0vQofAFCjJs/OhJpeByJ8AFD4YJvCIF3nJs9ZCTS9Dkb4AKDQwTalQbqKG5kNVax0zJgxI6mKRw/hA4DCBtsUB+k82Nm0NcIHAIUNtnUdpKvQ5JkS4QOAwgbbug7SVWjyTInwAUBhg22dB+nUmzxTYodTAArfNbQpu402yWo7nAKQ8q6hTdltlIKmXWLj0MUXXxymTp0aJkyYEHbZZZdw+eWXh8QKLABASXKvfFxzzTXhpptuCrfffnvYc889w7PPPhtOPfXUrBRzzjnn5P1xAEDTw8eTTz4Zjj322HD00Udnr3feeefw61//Ojz99NN5fxQAUEG5T7tMnz49LF68OLz22mvZ6xdffDE88cQT4aijjur369euXZs1qfR+AAD1lXvl44ILLsgCxG677ZYtoYo9IFdeeWU46aST+v36OXPmhMsuuyzvywAAmlL5uOeee8Kdd94ZFi1aFJ5//vms9+OnP/1p9tyfCy+8MFuW0/NYuXJl3pcEQM7qeDAcFd7nY9KkSVn146yzztrw3hVXXBHuuOOO8Morr2zy99vnAyBt8SC4nvNZ4m6lcdMwG2qxuoXxO/fKxyeffNLv1rl99/IHoHrqcDBc76pNXSo4qyr2feQePo455pisx+N3v/tdtr3sfffdF+bOnRu+/vWv5/1RABSs6gfDxarNlClTwsyZM8PkyZOzR/zn+F78tSpa2Ot7qsr3kfu0y5o1a7JNxmLoePfdd8NOO+0UvvWtb4VLLrkkjB07dpO/37QLQLrinXUc4HoHkFjdjjebZe5YGq8rBqN4cN1A19HftfeWwvdR5Z9HqdMuEydODPPmzQtvvPFG+Oc//xn+9re/ZT0fQwkeAKQtxYPhhnrn31/VpqoVnKpXohwsB0DLUjkYrpU7/9QrH6uGUL3p7/eofADQCHFgmzFjRulTFK3c+fet2sTFER0dHUlUcBYOs28jxUrUUKh8AFBZw7nz7121icqu4KzKoXqRQiWqlfE79x1OAaAoPXf+cblvrHgM5c4//lrvX0+5etM5xGvr+z2lTvgAoNLiBmddXV2l3/kP17Rp07IpoL6Vj57KTB3p+QCg8lLpQRmOzor2bYyEng8ASMCqBPo2RkLPBwClLwOlNZ0V69sYCdMuANRu+27SPgdG+ACg1gfJkV6QFD4AGnr3Weftu0k7SAofAA29+2xlGWhvdV8GWlfLEgqSwgdAxe4+i6yaNHEZaF1NSyhICh8AFbr7LKNqEjfxilt9x8ATn+NrqqczoSBpnw+AipzjkdIJplVlyXBo234iTrUFqOHdZ0pz9lVrnK1ar03dd4NV+QCoyN3ncCsf7brbj4N3T/9K7CWIoSrVKRlVo/ZT+QCo4d3ncKom7brbT2nZ5lCkWDVqMuEDoEJaaf5sZ0Co2mCe0koPhA+A2lZN2hkQqjaYp7TSA+EDoLbaGRCqOJhbMpwODacAFTGcxtHY4xGnWmLFoycg5DnoVv0YeMoZv4UPgAoYycoSAYEiCB8ANWKZKFVgqS1AjbSjcbRKm4NRP8IHQMMaR+30SdmED4DE5bmypGqbg1FPY8q+AIAyVO2Asdhc2tXVNeLG0cGmcKrw50A9qHwAjVPVaYc8DgSr2uZgdbWq4T03wgfQKE2fdqji5mB1s7Ci4TdPltoCjRLvNuP/9Pt7P1YVmsLeH+Wo87Lp1S2M33o+gEbpmXbo+z//pk07xIGu6oNdFem5+S/TLkCjmHagTHpu/kv4ABrHAWOURfj9Lz0fACWr2rJfRm5VDXtubK8OUJEllFY+NFNnDsumq0z4AMhJq0Gi6ct+aS7hAyAHwwkS7TgwDqpA+ADIwXCChJUPNJXwAZCD4QQJKx9oKuEDIAfDDRKW/dJEltoC5KiOSyhhKGyvDlAS25bDppl2AQAKJXwAAIUSPgCAQgkfAEChhA8AoFDCBwC11uphf1Q0fLz55pvh5JNPDttss02YMGFC2HvvvcOzzz7bjo8CgAE5Nbgh4eODDz4Ihx9+eNhss83CQw89FF5++eVw3XXXha222irvjwKAATk1OF25bzJ2zTXXhEmTJoVbb711w3tTp07N+2MAqKkYDuJBffG8nJFs2DbYYX82gqtZ5eOBBx4IBx10UDj++OPDdtttF/bff/9wyy23DPj1a9euzbZk7f0AoJk9GXlOkzg1uEHh4/XXXw833XRT9kP//e9/H84888xwzjnnhNtvv73fr58zZ062F3zPI1ZNAGiG3mFj8uTJ4bTTTsttmsSpwQ06WG7s2LFZ5ePJJ5/c8F4MH88880xYunRpv5WP+OgRKx8xgDhYDqDeYqiIwaPv1EhfsSoyY8aMEX2Ow/5qfrDcjjvuGPbYY4+N3tt9993Dvffe2+/Xjxs3LnsA0Cz99WT0lcc0icP+GjDtEle6vPrqqxu999prr2XpFgAG68no6OgwTdIAuYeP888/Pzz11FPhqquuyspcixYtyubczjrrrLw/CoAKb6LVX09GXKCwYsWK7PuJz7NmzSr7MqlCz0f04IMPhgsvvDArqcVltrNnz86aiPKeMwKo8lLQvBo2e/ayiFWEOJhXbcDWk1EPrYzfbQkfIyF8AClLabDvr2EzVg9ixcAgTsrjt7NdACq6Y+Zgm2hByoQPgIoO9jbRoqqED4CKDvY20aKqhA+Aigz2/a1qif0mVodQNRpOASqwOiOlRlfoj9UuQKOltBQ2D1a1UAVWuwCNleepqKlIrdEVRkr4AGojtaWwdW10hZESPoDaqGuFoOxGV8hb7qfaAuWpW6/DcCsEfXsj6lAhiM2lXV1dtiGnFlQ+oMJ6L72sY69Dq+peIYjfx4wZM2rz/dBcVrtARfVdehn/Kvf+69zk1RAOKoO0x2/TLlCTxsq+enodmjj4xu95qN9306eqoAymXaAmjZV91aXXoZ1MVUE5hA+oydLL+Lrnvbr1OrRDXZflQhUIH1CTxsr4+o033nDGR8OX5UIVaDiFCtNYOXy2LId82V4dGsLSy+Gr+7JcSJnKB9BoqkeQD0ttAdqwLHcoLN2FTTPtApATS3dhaEy7AMmrQjVBAytNt1rDKVAXVakmWLoLQyd8AMmq0kZg/W38ZpdZ6J/wAYz4RN12qVI1wdJdGDrhA0h2KqRq1YS4q2zs8bDLLAxO+ACSnQoZTjWhiIrMYGz8BpsmfEADbGpAbmXALnoqpJVqQlWaU6HphA+ouU0NyK0O2GVMhQylmlCl5lRoOuEDamxTA/JwBuxUGyur1JwKTSd8QI1takAe7oCdYmPlSCoyZfeJQNMIH1BjmxqQRzJgp9ZYOdyKjD4RKJ7wATW2qQE51SmU4Wq1IqNPBMrhbBdogE0dG1/EsfIpns8SQ0qsePT3fqzqpHjNkCpnuwAtTZG0ewol1amNwaadUr3mgehboUqED6CtUp7aGGjaKUr1mvtTtaAEwgfQ6CWw/fWJpH7NVQl3MJAxA/4KQI5TG70H89TOZ4kVkN5TTkVcc179JIMFJX0qpErlA2irKq6oafc15zlNUrXD9yCy2gUoRBEraqpwzfHfGQNH36pKnPIZ7mfE8BKnWmLFoycopbDxG82yuoXx27QLUMrURlOvuR3TJDFodHV1VS7c0VzCB0CB2tVPUsVwR3Pp+QDaxt4T9eiBgbzp+QDaIvYh9CwBjXf6ccDVh1DtHhjIa/wWPoDctaOpEkib7dWBUlVpky6geMIHkLuy957QawJpEz6AWjVVOucE0qfnA2qs7CPhi26q1GsC5dHzASRRAYgD/owZMwob+PWaQDW0PXxcffXVoaOjI5x33nnt/iig4Sedlt1rAiQQPp555plsnnefffZp58cAJVUAUmvstIEXNDx8fPTRR+Gkk04Kt9xyS9hqq63a9TFASRWAFKZ1+hM3Mos9HjEUxWcbm0GDwsdZZ50Vjj766HDkkUcO+nVr167NmlR6P4C0KwAjndZpd8Wk6F4TIIHwcdddd4Xnn38+zJkzZ5NfG78mdsf2PCZNmtSOS4LGaWcFYCTTOqlWTIAKL7VduXJlOOigg8Ijjzyyodcj3oHst99+Yd68ef1WPuKjR6x8xABiqS2ka7hLWi2Fhfoqdantc889F959991wwAEHhDFjxmSPxx57LFx//fXZP8e7o97GjRuXXWTvB1DPaR1LYYG2VD7WrFkT3njjjY3eO/XUU8Nuu+0WfvSjH4W99tpr0N9vkzGojlY3EVP5gPpqZfwek/eHT5w48TMBY/PNNw/bbLPNJoMHpKbsHUJTF/9MWvlz6amYxObUWPGwFBaayQ6nULHGyNT21miVpbCAs12gQtMDMQD1LHGN+3jEKkLP4K1KA5TJ2S4wQgM1Rv7sZz8r7ZoG21sj1SoNQH+EDxjiDqHR3LlzS5vuGCgQLV26tJHnuADVJXxAP+K0xezZsz/zfhzYy1oWOtCW6XHmtKrLV6vevwIMj/ABAzj33HOzE5lTOSF1oL01pk+fXsmTXE0VQXMJHzDIYB8PRkzphNT+VopU8STXkZ4NA1Sb1S6Q80ZaZanKdUYxPMWKR3/vx+MYgOopdZMxaPpGWkXpu7Q21escrH+l71Lm1KeKgHyYdoEKGm6/RCoNnlWcKgLyY9oFGrIB2mAblJWlSlNFwOBsMlZBqdyRkr7hnAybaoNnDByxx0PwgGYRPhJgySF57PcxWL+Eo+yBlAgfJUv1jrROUqoq5XEtw+mXGE5gAWgX4aNk7kibU1XK81paPRlWgyeQEg2nJUv19NQ6SOnPNpVr0eAJtIuG0wpxR9qMqlIq16LBE0iBTcYSEEvmXV1d7khrvJFVStcCUDaVj0S4I613VSl+5tVXX72h6VOFC2gyPR/U3lD6HPpuVZ63vht8XXPNNeH73/9+7p8DUIXxW/ig8dq982cqzaYA7aThFNq8z0or+3Wk0mwKkArhg0YbTjAYbL+O/kKJDb4ANiZ80GitBoPBKiUDhZKUGl8BUiB80GitBoOBKiVLly4ddPqm1R1JAerMPh80Xiv7rAy0X0fs2x5o+qbn3xefVTsAVD4aJ6VD1qq4z8pAlZLp06fr6wAYIuGjQVI6ZK3K+ptC0dcBMHT2+aiIkW6CZa+JYji4DWiq1fb5qJc8Khb2miiGbfIBNk34qOkmWH3Za6I+9O0AVSd8JC6vioWehHrQtwPUgZ6PxOXdq6Enobr07QAp0/NRI3lXLPQkVJe+HaAubDJWs02wqK+BNjjTtwNUjcpHRYykYqFBsR707QB1oeej5mJDYs9qmXjXHAcv54pUm74dIEWtjN/CR41pUASgKBpOyWhQBCBFwkeN2VgMgBQJHzWmQXF4NOgCtJfw0cATWOsqj9BgB1GA9tNwSi3ksapHgy7A8Gk4pVHyOnxPgy5AMYQPKi+v0KBBF6AYwkeFaYzMNzRo0AUohvBRURoj2xMaYp/I0qVLw9y5c7PnOjfoApRFw2kFaYxs37bjtqMHGB4NpzWnMTL/w/fybFwFYHDCRwVpjGwPoQ6gGMJHBRXRGNnEZlahDqCi4WPOnDnh4IMPDhMnTgzbbbddOO6448Krr76a98c0Xjt3Lm1qM6vVLgAVbTj96le/Gk488cQsgPznP/8JF110UXjppZfCyy+/HDbffPNN/n4Np+XSzJpP4ypA06xuYfwek/eHP/zwwxu9vu2227IKyHPPPReOOOKIvD+OAvsemjIQx+9zoO81BpP4ZxSnaJry5wFQuZ6PmICirbfeut9fX7t2bZaWej8oj76HgTV1OgqgUuEj3kGfd9554fDDDw977bXXgD0isUzT85g0aVI7L4lN0PfQP8twASqyydiZZ54ZHnroofDEE08MOHjFykd89IiVjxhA9HyUS9/DxmJjb6x49Pd+3FsEoOlWl9nz0ePss88ODz74YHj88ccHHbzGjRuXPZosxT6Cwfoemjwd1bcR13QUQALTLrGQEoPHfffdFx599NEwderUvD+iVpraR1C1fURMRwEkPO3y3e9+NyxatCjcf//94fOf//yG92MpZsKECZv8/U1aatvUZa1VPj/FdBTAyMfv3MNHR0dHv+/feuut4dvf/vYmf3+TwkcT+wjyClwpTlUBNNnqMg+Wi1mmv8dQgkfTNHFZax7npzR1qgqgLpztUuM+ghT7KkYauCx5Bag+4aOmZ7SkWh0YaeBy8ixA9bV1n4/haFLPR5MbWYfbuFmF7w2giVaX2fNB+apQHYhBITbVthoYLHkFqL62bTJGeeq+IVacmurq6rLkFaCiVD5qKNXqQJ4NsMOtnABQPj0fNZbShlhV3lgMgMQ3GRsp4aN8eW/gpUkUoP5WazglpSW6VWiABaA4wgdt38Ar5Z1cU9yIDaDuhI8aGu6A2q4KRaoNsKluxAZQd3o+amYkjZ3t7s1IqQFWHwpAvvR8NNRIp03aVaHoqcREqSyP1YcCUB7ho0byGFDzPmsm1amNlPtQAOpO+KhRE2NeA2peG3ilfAJtqn0oAE0gfOSszDv91AbU1Kc22nWiMACD03BawybGVBo7U/nzAKD9NJyWZKA7/aVLlxZ6HWWee9J7yim1SgwAaXCqbZtPk41OOOGELBHWvaw/0DJfJ9AC0JtplzYMwLGpMlY8eqv7dIMpFoBmW23apTzxTn/RokWfeT+lRst2rNZJvbkUgHQIH20wffr0yu4hMdzVOvbNAGCohI82qGqj5Uj25ajq9wxA8TSc5iwO1HEKIjZZxn6HKjVaDjZ1MpTr11wKwFAIH4kc6pbqap1Wp05i4BA6ABiMaZecpLyV+FCZOgGgCCofiUxZpMLUCQDtJnwkNGWRClMnALSTaZecmLIAgKGxw2lND3UDgFTHb9MuOTNl0f/S4zgt5c8FgMi0C8ntlgpAvZl2oS0cNAfQLKsdLEfZHDQHwECED9rCQXMADET4oC0sPQZgIHo+aCtLjwGaYbWltqTC0mMA+jLtAgAUSvgAAAolfAAAhRI+AIBCCR8AQKGEDwCgUMIHAFAo4QMAKJTwAQAUSvhIbCvyJUuWZM8AUFfCRyIWLlwYpkyZEmbOnJk9x9cAUEcOlktArHTEwLF+/foN78VTYFesWOFcFAAqoZXxu22VjxtvvDHsvPPOYfz48eHQQw8NTz/9dKijPKZKli1btlHwiNatW5edBgsAddOW8HH33XeH2bNnh0svvTQ8//zzYd999w1dXV3h3XffDXXqqchrqmTatGlh1KiNfxSx8hGPoQeAumlL+Jg7d2447bTTwqmnnhr22GOPcPPNN4fPfe5z4Ze//GWoS09FDDCnn376hopFfD7jjDOGFWzi1MqCBQuywBHF5/nz55tyAaCWcg8fn376aXjuuefCkUce+b8PGTUqe7106dLPfP3atWuzeaLej7zlGRTaNVUya9asrMcjVmbic3wNAHWUe/h4//33s0F4++233+j9+Prtt9/+zNfPmTMna1DpeUyaNCnvS2pLT0U7pkpipWPGjBkqHgDUWulLbS+88MKsM7bnsXLlytw/o11BwVQJALRuTMjZtttumw3E77zzzkbvx9c77LDDZ75+3Lhx2aOdeoJCnGqJFY+8gkKcGomNtLGCEoOM4AEAJVQ+xo4dGw488MCwePHiDe/FKY/4+rDDDgtlaVdPhakSACi58hHFZbannHJKOOigg8IhhxwS5s2bFz7++ONs9UuZYkAQEgCghuHjhBNOCO+991645JJLsibT/fbbLzz88MOfaUIFAJrH9uoAQD22VwcA6I/wAQAUSvgAAAolfAAAhRI+AIBCCR8AQKGEDwCgUMIHAFAo4QMAKJTwAQAUSvgAAAolfAAAhRI+GmjVqlVhyZIl2TMAFE34aJiFCxeGKVOmhJkzZ2bP8TUAFKmju7u7O1T0SF5aEysdMXCsX79+w3ujR48OK1asCJ2dnaVeGwDV1sr4rfLRIMuWLdsoeETr1q0Ly5cvL+2aAGge4aNBpk2bFkaN2vhHHisfu+66a2nXBEDzCB8NEqdWFixYkAWOKD7Pnz/flAsAhdLz0dDejzjVEiseggcARY/fY3L5RColBg6hA4CymHYBAAolfAAAhRI+AIBCCR8AQKGEDwCgUMIHAFAo4QMAKJTwAQAUSvgAAAolfAAAhRI+AIBCJXe2S885d/GAGgCgGnrG7aGcV5tc+FizZk32PGnSpLIvBQAYxjgeT7cdTEf3UCJKgdavXx/eeuutMHHixNDR0ZFLEotBZuXKlZs84pfi+Lmkyc8lXX42afJz+Z8YJ2Lw2GmnncKoUaOqVfmIF9yO497jfxRN/w8jRX4uafJzSZefTZr8XP5rUxWPHhpOAYBCCR8AQKFqHz7GjRsXLr300uyZdPi5pMnPJV1+Nmnycxme5BpOAYB6q33lAwBIi/ABABRK+AAACiV8AACFqn34uPHGG8POO+8cxo8fHw499NDw9NNPl31JjTZnzpxw8MEHZzvYbrfdduG4444Lr776atmXRR9XX311tsPweeedV/alNN6bb74ZTj755LDNNtuECRMmhL333js8++yzZV9Wo61bty5cfPHFYerUqdnPZJdddgmXX375kM40oQHh4+677w6zZ8/OlkE9//zzYd999w1dXV3h3XffLfvSGuuxxx4LZ511VnjqqafCI488Ev7973+Hr3zlK+Hjjz8u+9L4/5555pkwf/78sM8++5R9KY33wQcfhMMPPzxsttlm4aGHHgovv/xyuO6668JWW21V9qU12jXXXBNuuumm8POf/zz89a9/zV5fe+214YYbbij70iqj1kttY6Uj3mXH/0B6zo2Je/B/73vfCxdccEHZl0cI4b333ssqIDGUHHHEEWVfTuN99NFH4YADDgi/+MUvwhVXXBH222+/MG/evLIvq7Hi/6f++Mc/hj/84Q9lXwq9fO1rXwvbb799WLhw4Yb3vvGNb2RVkDvuuKPUa6uK2lY+Pv300/Dcc8+FI488cqNzY+LrpUuXlnpt/M+HH36YPW+99dZlXwohZFWpo48+eqO/N5TngQceCAcddFA4/vjjs5C+//77h1tuuaXsy2q86dOnh8WLF4fXXnste/3iiy+GJ554Ihx11FFlX1plJHewXF7ef//9bF4uptPe4utXXnmltOvif2IlKvYUxLLyXnvtVfblNN5dd92VTU/GaRfS8Prrr2fl/Th9fNFFF2U/m3POOSeMHTs2nHLKKWVfXqMrUvE029122y2MHj06G2uuvPLKcNJJJ5V9aZVR2/BBNe6yX3rppeyOgXLF48DPPffcrA8nNmeTTkCPlY+rrroqex0rH/HvzM033yx8lOiee+4Jd955Z1i0aFHYc889wwsvvJDdSMWj5P1cGh4+tt122yyRvvPOOxu9H1/vsMMOpV0X/3X22WeHBx98MDz++OOhs7Oz7MtpvDhFGRuxY79Hj3g3F38+sWdq7dq12d8nirXjjjuGPfbYY6P3dt9993DvvfeWdk2E8IMf/CCrfpx44onZ67gC6Y033shW8wkfDe/5iGXJAw88MJuX630XEV8fdthhpV5bk8X+5hg87rvvvvDoo49mS9Uo35e//OXwl7/8JbuD63nEO+5YRo7/LHiUI05J9l2KHvsMpkyZUto1EcInn3yS9RD2Fv+OxDGGhlc+ojhPGlNo/J/oIYccknXtxyWdp556atmX1uiplliqvP/++7O9Pt5+++3s/S233DLrFKcc8WfRt+9m8803z/aW0I9TnvPPPz9rbozTLt/85jezfYoWLFiQPSjPMccck/V4TJ48OZt2+fOf/xzmzp0bvvOd75R9adXRXXM33HBD9+TJk7vHjh3bfcghh3Q/9dRTZV9So8X/5Pp73HrrrWVfGn188Ytf7D733HPLvozG++1vf9u91157dY8bN657t912616wYEHZl9R4q1evzv5uxLFl/Pjx3f/3f//X/eMf/7h77dq1ZV9aZdR6nw8AID217fkAANIkfAAAhRI+AIBCCR8AQKGEDwCgUMIHAFAo4QMAKJTwAQAUSvgAAAolfAAAhRI+AIBCCR8AQCjS/wNNoBCspAV8jgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,y,'.k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.99109086])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "the0, the1 = MyOLS(x,y)\n",
    "the1"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv (3.9.6.final.0)",
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}