least squares_NeklyudovSA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear (ordinary) least squares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we have a set of points $\\left\\{ (x_i, y_i), i = 1, \\ldots, n \\right\\}$ and a functional form $f(x, \\vec \\beta)$, where $\\vec \\beta$ is a vector of parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, *least squares approximation* is a problem of finding $\\vec \\beta$ by minimising function\n",
    "$$ \\xi (\\vec \\beta) = \\sum\\limits_{i = 1}^n \\left| y_i - f(x_i, \\vec \\beta) \\right|^2 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Least squares approximation becomes *linear (or ordinary)* when we assume that $f(x, \\vec \\beta)$ is a linear function of $\\vec \\beta$: \n",
    "$$f(x, \\vec \\beta) = \\sum \\limits_{j = 1}^m \\beta_j \\varphi_j (x),$$\n",
    "where $\\varphi_j (x)$ are called basis functions which do not need to be linear ($x^k$, $e^{i \\alpha_k x}$ etc.)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We minimize with respect to $\\vec \\beta$\n",
    "$$ \\xi (\\vec \\beta) = \\sum\\limits_{i = 1}^n \\left| y_i - \\left( \\beta_1 \\varphi_1(x) + \\ldots + \\beta_m \\varphi_m(x) \\right) \\right|^2 .$$\n",
    "It can be represented as \n",
    "$$\\xi (\\vec \\beta) = || \\vec y - A \\vec \\beta ||^2,$$\n",
    "if we suppose that a vector norm $||\\vec a|| = \\sqrt{\\sum\\limits_{i = 1}^n |a_i|^2}$ and *the design matrix*\n",
    "$$ A = \\begin{pmatrix} \n",
    "            \\varphi_1(x_1) & \\varphi_2(x_1) & \\ldots & \\varphi_m(x_1) \\\\\n",
    "            \\varphi_1(x_2) & \\varphi_2(x_2) & \\ldots & \\varphi_m(x_2) \\\\\n",
    "            \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "            \\varphi_1(x_n) & \\varphi_2(x_n) & \\ldots & \\varphi_m(x_n) \\\\\n",
    "        \\end{pmatrix}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To minimize $\\xi (\\vec \\beta)$ we set \n",
    "$$\\cfrac{\\partial}{\\partial \\beta_k} \\xi (\\vec \\beta) = 0, \\; k = 1, \\ldots, m, $$\n",
    "and obtain *the normal equation*\n",
    "$$A^T A \\vec \\beta = A^T \\vec y.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $A$ is poorly conditioned, then solving this matrix equation becomes a big problem. In this case QR decomposition of $A$ is often used. The design matrix can be represented as \n",
    "$$A = Q \\cdot R,$$\n",
    "where $Q$ is orthoginal ($Q^T Q = 1$), and $R$ is upper triangular."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since $m < n$ last $n-m$ rows of $R$ are zero:\n",
    "$$A = Q \\cdot \\begin{pmatrix} \n",
    "                    R_1 \\\\\n",
    "                    0 \\\\\n",
    "                \\end{pmatrix}.$$\n",
    "Recall tht matrix $Q$ is a rotation operator, and introduced norm of a vector is invariant to rotation. Then, we rotate the residual $\\vec y - A \\vec \\beta$:\n",
    "$$\\xi (\\vec \\beta) = || \\vec y - A \\vec \\beta ||^2 = || Q^T \\left( \\vec y - A \\vec \\beta \\right) ||^2 = || Q^T \\vec y - \\begin{pmatrix} \n",
    "    R_1 \\\\\n",
    "    0 \\\\\n",
    "\\end{pmatrix} \\vec \\beta ||^2. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us write \n",
    "$$Q^T \\vec y = \\begin{pmatrix} \n",
    "                    \\vec f \\\\\n",
    "                    \\vec r \\\\\n",
    "                \\end{pmatrix},$$\n",
    "where $\\mbox{dim} \\vec f = m$.\n",
    "\n",
    "This way,\n",
    "$$\\xi ( \\vec \\beta ) = ||\\vec f - R_1 \\vec \\beta||^2 + ||\\vec r||^2 \\rightarrow \\min \\Leftrightarrow R_1 \\vec \\beta = \\vec f.$$\n",
    "\n",
    "The last equation is solved by back substitution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let the basis function be $x^k, \\; k = 0, \\ldots, m$. Write function that finds $\\vec \\beta$ for $m = 1, 2, 3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 338,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_least_squares(m, xdata, ydata):\n",
    "    A = np.zeros((len(xdata), m+1))\n",
    "    for i in range(len(xdata)):\n",
    "        A[i] = np.array([np.power(xdata[i], j) for j in range(0,m+1)])\n",
    "    q, r = np.linalg.qr(A)\n",
    "    beta = np.linalg.inv(r)@q.T@ydata\n",
    "    return beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now implement your function to approximate real data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 339,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X90VPWd//HnmwQSm8RfUKJA+aEGpCAQjFC1rLAUi7ZgqK4CX7vaI/K11F3hu9927e5pQen3fD2t54C7FStSa3v8Qbu6KLYU5atgqVjll4oU+SE/lsDWCFZhAgSGvL9/ZIKTyfy4dzJ35t4778c5OcnMfHLz/sydvDP5zL2vEVXFGGNM8ehS6AKMMcbklzV+Y4wpMtb4jTGmyFjjN8aYImON3xhjiow1fmOMKTLW+I0xpshY4zfGmCJjjd8YY4pMaaELSKZHjx7av3//QpfRTlNTExUVFYUuw1Nhn6PNL/jCPsfOzG/jxo2HVPXzTsb6svH379+fDRs2FLqMdtasWcPYsWMLXYanwj5Hm1/whX2OnZmfiOxzOtaWeowxpshY4zfGmCJjjd8YY4qML9f4kzl16hQNDQ2cOHGiID//nHPOYdu2bQX52fniZo7l5eX06dOHrl27elxV8dl3uInH1u7m+c0HaWqOUlFWSn1tL+4ccxH9uof3hU2TP4Fp/A0NDVRVVdG/f39EJO8//+jRo1RVVeX95+aT0zmqKocPH6ahoYEBAwbkobLisXp7I7Oe3MSp0y1EW1rfKyPSHGXpW/t5buMBFt06knGDeha4ShN0gVnqOXHiBN27dy9I0zftiQjdu3cv2H9fYbXvcBOzntzE8VOnzzT9NtEW5fip08x6chP7DjcVqEITFoFp/IA1fR+xfZF7j63dzanTLWnHnDrdwpK1e/JUkQmrwCz1uGFrpCaInt98sMMz/UTRFmXZ5gPMrx+ap6pMGAXqGb8Tq7c3MnHhWpa+tZ9IcxTlszXSiQvXsnp7Y9bbbmho4IYbbqCmpoaLL76Ye+65h5MnTyYde/DgQW666aaM27z++uv55JNPsqpn3rx5PPjggxnHVVZWpr39k08+YdGiRVnVYHKnqTnqbNxJZ+OMSSVUjd/LNVJV5Rvf+Ab19fXs3LmTHTt2EIlE+Nd//dcOY6PRKL169eLZZ5/NuN0VK1Zw7rnnuq4nl6zx+0NFmbN/wCu6hfIfdZNHGRu/iDwuIo0i8l6K278rIm/HPt4TkdMicn7str0isiV2m+cZDF6ukb722muUl5fzrW99C4CSkhIWLFjA448/zrFjx3jiiSf4u7/7OyZNmsS1117L3r17GTq09d/xY8eOcfPNNzNs2DBuueUWRo8efSaSon///hw6dIi9e/cyePBg7rzzToYMGcK1117L8ePHW+f12GNcccUVDB8+nBtvvJFjx46lrXXPnj1ceeWVXHHFFfzgBz84c30kEmH8+PGMHDmSyy67jBdeeAGAe++9lw8++ICrr76a7373uynHGW/V1/aitEv6105KuwhTanvnqSITVk6e8T8BTEx1o6r+RFVHqOoI4PvAa6r6cdyQcbHb6zpXamZu1kjd2rZtG5dffnm7684++2z69u3Lrl27AHjjjTf45S9/yauvvtpu3KJFizjvvPN49913+cEPfsDGjRuT/oydO3fyne98h61bt3Luuefy3HPPAfCNb3yD9evX88477zB48GB+/vOfp631nnvu4dvf/jbr16/nggsuOHN9eXk5y5YtY9OmTaxevZp/+qd/QlV54IEHuPjii3n99df5yU9+knKc8dadYy6ia0n6X8muJV2YMcYOoTWdk7Hxq+ofgI8zjYuZBjzTqYo6wcs1UlVNeiRL/PUTJkzg/PPP7zDmj3/8I1OnTgVg6NChDBs2LOnPGDBgACNGjADg8ssvZ+/evQC89957jBkzhssuu4ynnnqKrVu3pq319ddfZ9q0aQB885vfbFfrv/zLvzBs2DC+8pWvcODAAT788MOkc3IyzuRWv+4VLLp1JGd1LenwzL+0i3BW1xIW3TrSDlAwnZazNX4R+Ryt/xk8F3e1Ai+LyEYRmZmrn5WKl2ukgwcP7pAYeuTIEfbv38/FF1/cut0UcapOny2XlZWd+bqkpIRotPUP1O23385Pf/pTtmzZwty5cx0dP5/sj9RTTz3FRx99xMaNG3n77beprq5Oui2n40zujRvUk5WzxzBtVF8qy0oRgcqyUqaN6svK2WPs5C2TE7l8lWgS8HrCMs/VqnpQRHoCq0Tk/dh/EB3E/jDMBKiurmbNmjXtbj/nnHM4evRo2gK+NuTzPPf2fxNNs8xf2gW+PvTzGbeVaMyYMcybN49HH32U6dOnc/r0aWbPnn3m6xMnTnDy5Mkz241EIrS0tHD06FGuuOIKnnrqKerq6nj//ffZsmULTU1NHD16FFUlEom0Gw/Q3NxMc3MzR48e5ciRI1RVVfHxxx/zq1/9igsvvJCjR4/S3NxM165dO8xl9OjR/OIXv2Dq1KlnloWOHj3Khx9+yLnnnsuJEyd4+eWX2bdvH5FIhMrKSo4cOcLp06fTjkv8OSdOnOiwn/wsEokEpt7x58L4cWVA25OBQ+zZcoh0r04FaX7ZCvsc8zW/XDb+qSQs86jqwdjnRhFZBowCkjZ+VV0MLAaoq6vTxEzqbdu2ZYwTmDV+EMu3NBJtOZ1yTNeSEr79t4OoqnL37/LRo0d54YUXmDVrFg8++CAtLS1cf/31PPjgg5SVlVFeXk63bt3O1FhZWUmXLl2oqqpizpw53HbbbVx99dXU1tYybNgwevXqRVVVFSJy5nDLtvHQ+uz/1KlTVFVV8aMf/Yjx48fTr18/LrvssjPRCmVlZZSVlXW4Xx5++GGmT5/O4sWLufHGGwGoqqrijjvuYNKkSYwbN44RI0Zw6aWXUllZSf/+/fnyl7/MVVddxde+9jX++Z//Oem4xJ9TXl5ObW2tq/uxkCzLPfhyOUc/nu+Tr30oTpYhRKQ/8FtVTXrWiIicA+wBvqCqTbHrKoAuqno09vUq4H5VXZnp59XV1Wnissq2bdsYPHhwxlqTZZ1A6xpp15IuWWeddCar5/Tp05w6dYry8nI++OADxo8fz44dO+jWrVtW2/OK2zk63Sd+EfbGGPb5Qe7m6FWf6KxOvhHLRqcH0WR8xi8izwBjgR4i0gDMBboCqOrPYsOmAC+3Nf2YamBZbK25FHjaSdPvrLY10iVr97Bs8wGaTkap6FbKlNrezBgzoCB/yY8dO8a4ceM4deoUqsojjzziu6ZvTLGIP98nUbRFiba0nu+zcvaY0L6QnrHxq+o0B2OeoPWwz/jrdgPDsy2sM/p1r2B+/VDfnNZeVVXlu7eSNKZYuTnfxy89JNcCdeauHUvuH7YvTFB5eb5PUASm8ZeXl3P48GFrOD7QlsdfXl5e6FKMcc0ykQKUztmnTx8aGhr46KOPCvLzT5w4EfpG52aObe/AZUzQVJSVEnHQ/MOciRSYmXXt2rWg7/a0Zs2aQB26mI1imKMx9bW9WPrW/rTLPWHPRArMUo8xxuSCZSJZ4zfGFBnLRLLGb4wpQsWeiRSYNX5jTEfrdh1i7vKt7GyMtF6x8nfU9KzkvslDuOqSHnmrw4/xB5n47XyffLJn/MYE1EOv7GD6kjc/a/oxOxsjTF/yJg+9siMvdXj5dqfGG9b4jQmgdbsOsWDVzrRjFqzaybpdhzytw8u3OzXescZvTADNXZ7+zXjazHvR2bhsefl2p8Y71viNCaDE5Z1UdnzobFy2LP4gmKzxG2OyZvEHwWSN3xiTNS/f7tR4xxq/MQFU07PS0biB1c7GZau+tleHk6AShT3+IIis8RsTQPdNHuJo3LxJzsZly+IPgskavzEBdNUlPZgzoSbtmDkTajw/icviD4LJGr8xAXXP+IE8PWN0h+WcgdWVPD1jNPeMH5iXOoo9/iCI7BUXYwLsqkt68PKca4DCvtl6MccfBJE1fhM4XuXC+CVvxk0d8WMjzVEqV7+Us5qzrcPJfeeHOXq5v/3yWEpFMr2VoYg8DnwdaFTVDn/ORWQs8ALQdmref6rq/bHbJgIPASXAElV9wElRdXV16rc3Jy/ks6l8CcIcV29vZNaTmzh1uqXdiUOlXYSuJV1YdOvIlEsL6ebXme3mkps6vKzZyzr8MEe/3HeJOvM7KCIbVbXOyVgna/xPABMzjFmrqiNiH21NvwR4GLgO+CIwTUS+6KQoY5LxKhfGL3kzburwsmYv6/DDHP1y3xVSxsavqn8APs5i26OAXaq6W1VPAkuBG7LYjjGAd7kwfsmbcVOHlzV7WYcf5uiX+66QcnVUz5Ui8o6I/F5E2g4c7g3sjxvTELvOmKx4lQvjl7wZN3V4WbOXdfhhjn657wopFy/ubgL6qWpERK4HngdqgGSn86W8R0RkJjAToLq6mjVr1uSgtNyJRCK+qynX/D7HiNNcmOZo0nmkml9nt5srbupI31raj3Vbs5d1+GGOXu5vrx6judbpxq+qR+K+XiEii0SkB63P8L8QN7QPcDDNdhYDi6H1xV2/vcgYhBc+O8vvc6xc/ZKjX6yKstKk80g1v85uN1fc1AHOmkw2NXtZhx/m6OX+9uoxmmudXuoRkQtERGJfj4pt8zCwHqgRkQEi0g2YCizv7M8zxcurXBi/5M24qcPLmr2sww9z9Mt9V0gZG7+IPAO8AQwSkQYRuUNE7hKRu2JDbgLeE5F3gH8DpmqrKHA38BKwDfiNqnr7rhAm1LzKhfFL3oybOrys2cs6/DBHv9x3heTkqJ5pqnqhqnZV1T6q+nNV/Zmq/ix2+09VdYiqDlfVL6nqurjvXaGqA1X1YlX9P15OxISfV7kwfsmbcVOHlzV7WYcf5uiX+66QMp7AVQh2AldhBGWO+w43sWTtHpZtPkDTySgV3UqZUtubGWMGpP2FyjS/bLeba27qaDc2doZormrOug4H950f5ujl/vbqMZqOmxO4rPE7FJSm2Blhn6PNL/jCPsd8NX7L6jGe8HtWSbFat+sQc5dvbfeevTU9K7lv8hDPI5zzwS/5QtnWnOu8pVQsltnk3OrtjUxcuJalb+0nEjseO9IcZelb+5m4cC2rtzcWusSi9NArO5i+5M0Ob9S+szHC9CVv8tArOwpUWW64edy5fYx69ZhO3C452m4m1vhNTgUlq6TYrNt1iAWrdqYds2DVTtbtOpSninLLL/lCXtWca9b4TU4FJauk2Mxd7uxI6nkvBvOIa7/kC3lVc65Z4zc5FZSskmKTuLyTyo4PnY3zG7/kC3lVc67Zi7smp5qcZpWcdDbOGCdcPe4cHsjY9hj16jFdyN8Ve8ZvcqotYyXjuG72nMPkjpvHndvHqFeP6UL+rljjNzkVlKySYlPTszLzIOjwxu1B4Zd8Ia9qzjVr/CangpJVUmzumzwk8yBg3iRn4/zGL/lCXtWca9b4TU4FJauk2Fx1SQ/mTKhJO2bOhJrAnsTll3whr2rONWv8JufGDerJytljmDaqL5VlpYhAZVkp00b1ZeXsMXl503LT0T3jB/L0jNEdlnMGVlfy9IzR3DN+YIEqyw03jzu3j1GvHtMdtkt+flfsFTbjiX7dK5hfP5T59UMLWkeHiIKVv0sZUbBsUwNzl2/lyInPjqI4u7yU+yYPYcrIPh227Xa8H/Q+7yxGDTifg5+cOBM7MGrA+fQ+76wOY4MYfwCgbYftaMLlBG4fo149puO3m68sIgtpcyjs4VAQvjk+9MqOtGerzplQc+ZZ7uxfb+b5zSnfII762l4svKX2zGW34/Mh0/5bvb2RWU9u4tTplnbHj5d2EbqWdGHRrSPPPMP0amw2453OsTPb9Yt8hbTZUo8JJTcRBcs2NaRt4tB6ss2yTQ0Arsf7gVeRBhZ/EEzW+E0ouYkocBtnEMT4A68iDSz+IJis8ZtQchNREL9Gn86nx1vHuR3vB15FGlj8QTDZi7vGFAHPIg0s/iCQ7Bm/MUXAq0gDiz8IJmv8JpTcRBScXe6sGZxzVus4t+P9wKtIA4s/CKaMjV9EHheRRhF5L8Xt/0NE3o19rBOR4XG37RWRLSLytoj46/hME2puIgrcxhkEMf7Aq0gDiz8IJifP+J8AJqa5fQ9wjaoOA+YDixNuH6eqI5weX2pMLriJKJgysg/1tb3Sjq2v7XXmpCy34/3Aq0gDiz8IpoyNX1X/AHyc5vZ1qvrX2MU/Af55tJui5iaiYOEttSy4eXiH5Zlzziplwc3DO5yM5Xa8H3gVaeDb+AOLCknJ0Zm7ItIf+K2qpj1XWUT+N3Cpqs6IXd4D/JXW1/4fVdXE/wbiv3cmMBOgurr68qVLlzqcQn5EIhEqK4MZWetU2Odo8wu+sM+xM/MbN26c4zN3c/bqk4iMA+4Avhx39dWqelBEegKrROT92H8QHcT+KCyG1sgGv0UHhC3OIJkwzjE+EybSLFSWNecsE8ZN3kyHzCBImRmULSf7z8uMnHwI42M0Xr7ml5OjekRkGLAEuEFVD7ddr6oHY58bgWXAqFz8PGOcWL29kYkL17L0rf1EYsd5R5qjLH1rPxMXrmX19sacbVvTbPuhV3YwfcmbHU4q29kYYfqSN3nolR1Z1+FVzSbcOt34RaQv8J/AN1V1R9z1FSJS1fY1cC2Q9MggY3LNy+wWN9t2kxnkJcuyMfGcHM75DPAGMEhEGkTkDhG5S0Tuig35IdAdWJRw2GY18EcReQd4C/idqq70YA7GdOBldoubbfsl18eybEy8jGv8qjotw+0zgBlJrt8NDO/4HcZ4z012i9t8dTfbjjiMEtjxobNsoWx5eX+Y4LEzd00oeZndEsRcmCDWbLxjjd+EkpfZLUHMhQlizcY71vhNKHmZ3eJm224yg7xkWTYmnjV+E0peZre42bZfcn0sy8bEs8ZvQsnL7BY323aTGeQly7Ix8azxm9DqkN1C7rJb3OTCuMkM8pJl2Zg29kqOCbV+3SuYXz+U+fVDc346fPy2M7nqkh68POeanP3sbLmp2YSXNX7jWBBzXtpn9USpXP1SypqDOD9jsmGN3ziyensjs57cxKnTLWdOBGrLeXlu4wEW3TrSd0sFbmoO4vyMyZat8ZuMgpjz4qbmIM7PmM6wxm8yCmLOi5uagzg/YzrDGr/JyE3Oi1+4qTmI8zOmM2yN32QUxJwXVzVnfhO6z8YaEwL2jN9kFMScFzc1B3F+xnSGNX6TURBzXtzUHMT5GdMZ1vhNRkHMeXFTcxDnZ0xnWOM3GQUx58VNzUGcnzGdYY3fOBLEnBc3WT1BnJ8x2bJXq0LGbezAul2HmLt8KzsbY2/9t/J31PSs5L7JQzokRgY150Vjh+1owuVEbufn5r72Kg7CTSSFMW1E1eGxbHlUV1enGzZsyDwwj3Id8OWFZLED0Lpc0bWkS4fYgYde2cGCVTtTbm/OhJq8JUd6we394dW2varDy/n5VRB+DzujM/MTkY2qWudkrKOlHhF5XEQaReS9FLeLiPybiOwSkXdFZGTcbbeJyM7Yx23OpmDcchs7sG7XobRNH2DBqp2s23XIs5q95GUMgx/iICxmwnSG0zX+J4CJaW6/DqiJfcwEHgEQkfOBucBoYBQwV0TOy7ZYk5rb2IG5y7c62u68F52N8xsvYxj8EAdhMROmMxw1flX9A/BxmiE3AL/SVn8CzhWRC4GvAqtU9WNV/SuwivR/QEyW3MYOnFnTz2DHh87G+Y2XMQx+iIOwmAnTGbl6cbc3sD/uckPsulTXdyAiM2n9b4Hq6mrWrFmTo9JyIxKJ+K6meBGnEQXNUdfz8PO8U/Hy/nCzbaevoLmtw8v5+Znffw87K1/zy1XjT3bao6a5vuOVqouBxdD64q7fXsDx+4tKlatfctQMKspKW+ex8neOt+3neafi+v7waNvgrEm7rcPL+fmZ338POytf88vVcfwNwBfiLvcBDqa53uSY29iBmp6Vace2SXyf2KDwMobBD3EQFjNhOiNXjX858Pexo3u+BHyqqv8NvARcKyLnxV7UvTZ2nckxt7ED900e4mi78yY5G+c3XsYw+CEOwmImTGc4PZzzGeANYJCINIjIHSJyl4jcFRuyAtgN7AIeA2YBqOrHwHxgfezj/th1Jsfcxg5cdUkP5kyoSbvNORNqOpzEFRRexjD4IQ7CYiZMZ9gJXA4FZW1x3+Emlqzdw7LNB2g6GaWiWylTanszY8yAlGfuzntxa7ujdwZWVzJvUsczd4Oo3f0RO2M23f2R9bYz3Ndu94sf5udHQfk9zFa+TuCyxu9Q2B9wEP452vyCL+xz9NWZu8YYY8LDGr8xxhQZa/zGGFNkrPEbY0yRscZvjDFFxhq/McYUGWv8xhhTZKzxG2NMkbHGb4wxRcYavzHGFBlr/MYYU2Ss8RtjTJGxxm+MMUXGGr8xxhQZa/zGGFNkrPEbY0yRscZvjDFFxhq/McYUGWv8xhhTZBw1fhGZKCLbRWSXiNyb5PYFIvJ27GOHiHwSd9vpuNuW57J4Y4wx7pVmGiAiJcDDwASgAVgvIstV9c9tY1R1Ttz4fwBq4zZxXFVH5K5kY4wxneHkGf8oYJeq7lbVk8BS4IY046cBz+SiOGOMMbknqpp+gMhNwERVnRG7/E1gtKrenWRsP+BPQB9VPR27Lgq8DUSBB1T1+RQ/ZyYwE6C6uvrypUuXZj0pL0QiESorKwtdhqfCPkebX/CFfY6dmd+4ceM2qmqdk7EZl3oASXJdqr8WU4Fn25p+TF9VPSgiFwGvisgWVf2gwwZVFwOLAerq6nTs2LEOSsufNWvW4Leaci3sc7T5BV/Y55iv+TlZ6mkAvhB3uQ9wMMXYqSQs86jqwdjn3cAa2q//G2OMyTMnjX89UCMiA0SkG63NvcPROSIyCDgPeCPuuvNEpCz2dQ/gauDPid9rjDEmfzIu9ahqVETuBl4CSoDHVXWriNwPbFDVtj8C04Cl2v5Fg8HAoyLSQusfmQfijwYyxhiTf07W+FHVFcCKhOt+mHB5XpLvWwdc1on6jDHG5JiduWuMMUXGGr8xxhQZa/zGGFNkrPEbY0yRscZvjDFFxhq/McYUGWv8xhhTZBwdx2/Ca9/hJh5bu5vnNx8k0hylcvVL1Nf24s4xF9Gve0WhyzPGeMAafxFbvb2RWU9u4tTpFqItrSdcR5qjLH1rP89tPMCiW0cyblDPAldpjMk1W+opUvsONzHryU0cP3X6TNNvE21Rjp86zawnN7HvcFOBKjTGeMUaf5F6bO1uTp1uSTvm1OkWlqzdk6eKjDH5Yo2/SD2/+WCHZ/qJoi3Kss0H8lSRMSZfrPEXqabmqLNxJ52NM8YEhzX+IlVR5ux1/Ypu9vq/MWFjjb9I1df2orRLsnfV/ExpF2FKbe88VWSMyRdr/EXqzjEX0bUk/e7vWtKFGWMG5KkiY0y+WOMvUv26V7Do1pGc1bWkwzP/0i7CWV1LWHTrSDuJy5gQssZfxMYN6snK2WOYNqovlWWlCFBZVsq0UX1ZOXuMnbxlTEjZK3dFrl/3CubXD2V+/VDWrFnD2LFjU46Nj3doao5SUVZq8Q7GBJA1fuOIxTsYEx6OlnpEZKKIbBeRXSJyb5LbbxeRj0Tk7djHjLjbbhORnbGP23JZvMkPi3cwJlwyNn4RKQEeBq4DvghME5EvJhn6a1UdEftYEvve84G5wGhgFDBXRM7LWfUmLyzewZhwcfKMfxSwS1V3q+pJYClwg8PtfxVYpaofq+pfgVXAxOxKNYVi8Q7GhIuTNf7ewP64yw20PoNPdKOI/A2wA5ijqvtTfG/SM4JEZCYwE6C6upo1a9Y4KC1/IpGI72rKtVRzjDiNd2iO+vo+Cvs+DPv8IPxzzNf8nDT+ZKd3Jj79exF4RlWbReQu4JfA3zr83tYrVRcDiwHq6uo03dElhZDpiJcwSDXHytUvOWr+FWWlvr6Pwr4Pwz4/CP8c8zU/J0s9DcAX4i73AQ7GD1DVw6raHLv4GHC50+81/mfxDsaEi5PGvx6oEZEBItINmAosjx8gIhfGXZwMbIt9/RJwrYicF3tR99rYdSZALN7BmHDJ2PhVNQrcTWvD3gb8RlW3isj9IjI5NuwfRWSriLwD/CNwe+x7Pwbm0/rHYz1wf+w6EyAW72BMuDg6gUtVVwArEq77YdzX3we+n+J7Hwce70SNxgfa4h2WrN3Dss0HaDoZpaJbKVNqezNjzABr+sYEiJ25axyLj3cwxgSXNf4iF5+/E2mOUrn6pZT5O+t2HWLu8q3sbIycua6mZyX3TR7CVZf0yFvNyzY1MHf5Vo6c+OxIo7PLS7lv8hCmjOyTtzosu8gElaVzFrHV2xuZuHAtS9/af+Zwzbb8nYkL17J6e+OZsQ+9soPpS95s1/QBdjZGmL7kTR56ZUdeap79683M+c077Zo+wJETUeb85h1m/3pzXupIvO+U1PedMX5jjb9IucnfWbfrEAtW7Uy7vQWrdrJu1yEvS2bZpgae35z+aODnNx9k2aYGT+uw7CITdNb4i5Sb/J25y7c62ua8F52Ny5Zf6rDsIhN01viLlJv8ncTlnVR2fOhsXLYSl3dS+fS4s3HZsuwiE3TW+ItUk9P8nZPeNtEgsvvOBJ01/iJVUebsgK6KbnbgVyK770zQWeMvUm7yd2p6Vjra5sBqZ+OydXa5s0Z6zlneNlzLLjJBZ42/SLnJ37lv8hBH25w3ydm4bPmlDssuMkFnjb9IucnfueqSHsyZUJN2e3Mm1Hh+EteUkX2or+2Vdkx9bS/PT+Ky7CITdNb4i1hb/s60UX2pLCtFgMqyUqaN6svK2WPavXn6PeMH8vSM0R2WcwZWV/L0jNHcM35gXmpeeEstC24e3mE555yzSllw83AW3lKblzo63HeS+r4zxm/s1aciF5+/k+lNIK66pAcvz7kmf8WlMGVkn7xGM6Ri2UUmqKzxp+EmxyZxfC6zW7zMhAliVo9fWFaPCSpb6knBTY5NsvG5ym7xMhMmiFk9fmFZPSbIrPEn4TaLxavsFi8zYYKY1eMXltVjgs4afxJus1i8ym7xMhMmiFk9fmFZPSborPEn4TaLxavsFi8zYYKY1eMXltVjgs4afxJus1i8ym5wQK/SAAAKOUlEQVTxMhPG8mayZ/edCTpr/Em4zWLxKrvFy0wYy5vJnt13JugcNX4RmSgi20Vkl4jcm+T2/yUifxaRd0XkFRHpF3fbaRF5O/axPJfFe8VtFotX2S1eZsIEMavHLyyrxwRdxsYvIiXAw8B1wBeBaSLyxYRhm4E6VR0GPAv8OO6246o6IvYxOUd1e8ptFotX2S1eZsIEMavHLyyrxwSdk2f8o4BdqrpbVU8CS4Eb4geo6mpVPRa7+Ceg8KdVdoLbLBavslu8zIQJYlaPX1hWjwk6UU1/dIKI3ARMVNUZscvfBEar6t0pxv8U+Iuq/ih2OQq8DUSBB1T1+RTfNxOYCVBdXX350qVLs5tRDjUea+GlvadYdzDKiahSXipc1auUr/bvSs/Pdfyb2X48lJeSdnx2deRuux23nX6O2w5HeXLbSQ5EPnvM9K4Ubh3cjcHd/b+eHYlEqKzM3XKUl/slG7menx+FfY6dmd+4ceM2qmqdk7FOfluTLWYm/WshIrcCdUB8oEtfVT0oIhcBr4rIFlX9oMMGVRcDiwHq6uo0XWZMZyzb1MDc5VvbvY3f2eWl3Dd5SNL8l5tjnzPl2LRte/1bW2l757/jUVjfCF//0kDGJmzbTfyBm+26jRH48e+38ep/7Y7tUOF4FF79ryiDBvTl5rGD243966YGjm7ZSuvf8FZHoyVcMODSTs3Pbd1utt1+rABNOY2ZuDnzECA/8Q5OHqNBF/Y55mt+Tp7xXwnMU9Wvxi5/H0BV/2/CuK8A/w5co6pJz1cXkSeA36rqs+l+Zl1dnW7YsMHpHByb/evNPL/5YMrb62t7pUx3zLRD3Gz7oVd2pD0Tds6EmjNpl262u3p7I7Oe3MSp0y3tjjMv7SJ0LenColtHtkuNnLjwNd7/S+pj7y+9oJKVs6/xdH5u63azbbd1eMXtfslW2JsihH+OnZmfiDh+xu/k/9H1QI2IDBCRbsBUoN3ROSJSCzwKTI5v+iJynoiUxb7uAVwN/NnZNHJr2aaGtI0LWk/MWbapwdNtu4k/cLNdtzECP/79trRNH+D9v0T48e+3eTY/8C46wi8xExbvYPwo41KPqkZF5G7gJaAEeFxVt4rI/cAGVV0O/ASoBP5DRAD+K3YEz2DgURFpofWPzAOqWpDG7yZ2wG3kr5tt96wqdzz2L5+ecDx20vBejmME5tcP5ZHXdjva9s9e202lw7c8dDu/l+dc4yr+4E+7DzvedoZ/ZDvU4RU387N4Z5Mvjn6jVXUFsCLhuh/Gff2VFN+3DrisMwXmSvyafjqfHnd/tqWbbX96PPfxB58ej7qKEZhfPzT5izRJtODt/NzUHXF4xqyb+87rmAm3+8WYfLAzd0MiqDECQa3bqbDPzwSTNf6QCGqMQFDrdirs8zPBVDSN/2yH69SJ7+Wa6227iT9ws123MQLpR36mC97ND7yLjvBLzITFOxg/KprG72XsgJttezXWbYzAt6+5yNG277rmIs9qBu+iI/wSM2HxDsaPiqbxTxnZh/raXmnH1Nf2yupNvN1s2038gZvtuo0R+N51g7n0gvTPdi+9oJLvXTfYs/mBd9ERfomZsHgH40cZT+AqBK9O4ILWY+7nvbi13dE755xVyrxJyc/cbeP0zF2n21636xDzXtza7qiSgdWVzJuU/Mxdp9vdd7iJJWv3sGzzAZpORqnoVsqU2t7MGDMg5Zm7P3ttN/EHHHah9Zn+965rf+auV/NzW7ebbbutwytu90s2wn5yE4R/jvk6gavoGn+2wv6Ag/DP0eYXfGGfY74afygOJchHDkohuZmf2ywiP9RsjMmvwDf+ZDkokeYoS9/az3MbD+QsB6VQ3MwvVZ7OkRNR5vzmHV7b+VHKLKJC1WyMyb9Av7gb9hwUN/PzMovIjbDvE2PCINCN300OShC5mZ+bvCAvhX2fGBMGgW78bnJQgsjN/LzMInIj7PvEmDAIdOMPew5KEOcXxJqNKTaBbvxhz0EJ4vyCWLMxxSbQjT/sOShu5udlFpEbYd8nxoRBoBt/2HNQvMqx8VLY94kxYRDoxh/2HBQ38/Myi8iNsO8TY8Ig0I0fYNygnqycPYZpo/pSWVaKCFSWlTJtVF9Wzh4T+BOF3Mxv4S21LLh5eIflnHPOKmXBzcPzcvKW25qNMfkXilfY+nWvYH790NC+dZ2b+U0Z2Sev0QyphH2fGBNkgX/Gb4wxxh1r/MYYU2Ss8RtjTJHxZR6/iHwE7Ct0HQl6AIcKXYTHwj5Hm1/whX2OnZlfP1X9vJOBvmz8fiQiG5y+yUFQhX2ONr/gC/sc8zU/W+oxxpgiY43fGGOKjDV+5xYXuoA8CPscbX7BF/Y55mV+tsZvjDFFxp7xG2NMkbHGn4SIlIjIZhH5bZLbbheRj0Tk7djHjELU2BkisldEtsTq35DkdhGRfxORXSLyroiMLESd2XIwv7Ei8mncPvxhIerMloicKyLPisj7IrJNRK5MuD3Q+w8czTGw+1BEBsXV/baIHBGR2QljPN2Hocjq8cA9wDbg7BS3/1pV785jPV4Yp6qpjhe+DqiJfYwGHol9DpJ08wNYq6pfz1s1ufUQsFJVbxKRbsDnEm4Pw/7LNEcI6D5U1e3ACGh9kgkcAJYlDPN0H9oz/gQi0gf4GrCk0LUU0A3Ar7TVn4BzReTCQhdlQETOBv4G+DmAqp5U1U8ShgV6/zmcY1iMBz5Q1cQTVj3dh9b4O1oIfA9oSTPmxti/X8+KyBfyVFcuKfCyiGwUkZlJbu8N7I+73BC7LigyzQ/gShF5R0R+LyLevjtNbl0EfAT8IrYcuUREEt/cIOj7z8kcIbj7MN5U4Jkk13u6D63xxxGRrwONqroxzbAXgf6qOgz4f8Av81Jcbl2tqiNp/XfyOyLyNwm3J3vvxCAd/pVpfptoPb19OPDvwPP5LrATSoGRwCOqWgs0AfcmjAn6/nMyxyDvQwBiS1iTgf9IdnOS63K2D63xt3c1MFlE9gJLgb8VkSfjB6jqYVVtjl18DLg8vyV2nqoejH1upHVtcVTCkAYg/j+ZPsDB/FTXeZnmp6pHVDUS+3oF0FVEeuS90Ow0AA2q+mbs8rO0NsnEMYHdfziYY8D3YZvrgE2q+mGS2zzdh9b446jq91W1j6r2p/VfsFdV9db4MQnrbJNpfRE4MESkQkSq2r4GrgXeSxi2HPj72JEFXwI+VdX/znOpWXEyPxG5QEQk9vUoWn8PDue71myo6l+A/SIyKHbVeODPCcMCu//A2RyDvA/jTCP5Mg94vA/tqB4HROR+YIOqLgf+UUQmA1HgY+D2QtaWhWpgWex3phR4WlVXishdAKr6M2AFcD2wCzgGfKtAtWbDyfxuAr4tIlHgODBVg3Um4z8AT8WWCnYD3wrR/muTaY6B3oci8jlgAvA/467L2z60M3eNMabI2FKPMcYUGWv8xhhTZKzxG2NMkbHGb4wxRcYavzHGFBlr/MYYU2Ss8RtjTJGxxm+MMUXm/wNbr/sTmdr9twAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn import datasets\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "diabetes = datasets.load_iris()\n",
    "X = diabetes.data\n",
    "y = diabetes.target\n",
    "\n",
    "ydata = X[:100,-1]\n",
    "xdata = X[:100,0]\n",
    "\n",
    "plt.plot(xdata, ydata, 'o', label='Original data', markersize=10)\n",
    "plt.grid(True)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 340,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-5.9188815 ,  1.74339767, -0.09341773])"
      ]
     },
     "execution_count": 340,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta = np.asarray(linear_least_squares(2, xdata, ydata))\n",
    "beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the resulting functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 341,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X90VPW57/H3kwQSmkRFkChQfqgBKYgEI1TbHOFQLHoLBrUKHHvtuaJLOd4C99ye2tPbgtq7rrd1LfAcpYrU4+nyB3prQTzHoiwkLS1W+S2mCESQEmiNYBUmQMgkz/0jE5wkk8neyeyZvWc/r7WyyMx8882zvxmeTPbs/dmiqhhjjAmPnEwXYIwxJr2s8RtjTMhY4zfGmJCxxm+MMSFjjd8YY0LGGr8xxoSMNX5jjAkZa/zGGBMy1viNMSZk8jJdQCL9+/fXYcOGZbqMbqmvr6ewsDDTZWSUrYGtQdi3H9K/Blu3bj2qqhc4GevLxj9s2DC2bNmS6TK6paqqikmTJmW6jIyyNbA1CPv2Q/rXQEQOOh1ru3qMMSZkrPEbY0zIWOM3xpiQ8eU+/kQaGxupra3l9OnTmS4lqXPPPZfdu3dnugzPFRQUMHjwYHr16pXpUrLOwWP1PLVxP6u3H6G+IUphfh6VZQO5q+JihvYL9xumJjUC0/hra2spLi5m2LBhiEimy+nUiRMnKC4uznQZnlJVjh07Rm1tLcOHD890OVllw5465j27jcamZqLNLdfKiDREWfnOIV7eephlt49n8sgBGa7SBF1gdvWcPn2afv36+brph4WI0K9fP9//9RU0B4/VM+/ZbZxqbDrb9FtFm5VTjU3Me3YbB4/VZ6hCky0C0/gBa/o+Yj+L1Htq434am5qTjmlsambFxgNpqshkq8Ds6nHD9pGaIFq9/UiHV/rtRZuVVdsP81DlmDRVZbJRoF7xO7FhTx3Tlm5k5TuHiDREUT7fRzpt6UY27Knr9ty1tbXceOONlJaWcskllzB//nzOnDmTcOyRI0e45ZZbupzzhhtu4NNPP+1WPYsXL+aRRx7pclxRUVHSxz/99FOWLVvWrRpM6tQ3RJ2NO+NsnDGdyarG7+U+UlXlpptuorKykn379rF3714ikQg/+MEPOoyNRqMMHDiQX/7yl13O+9prr3Heeee5rieVrPH7Q2G+sz/AC3tn5R/qJo26bPwi8rSI1InIe508/l0R2RH7eE9EmkTk/NhjH4rIrthjnmcweLmP9M0336SgoIC///u/ByA3N5clS5bw9NNPc/LkSZ555hm++c1vcuutt3Ldddfx4YcfMmZMy5/jJ0+e5NZbb2Xs2LHcdtttTJw48WwkxbBhwzh69Cgffvgho0aN4q677mL06NFcd911nDp1qmW7nnqKq666iiuuuIKbb76ZkydPJq31wIEDXH311Vx11VX88Ic/PHt/JBJhypQpjB8/nssvv5xXXnkFgPvvv58PPviAcePG8d3vfrfTccZblWUDyctJ/t5JXo4ws2xQmioy2crJK/5ngGmdPaiqP1XVcao6Dvg+8BtV/SRuyOTY4+U9K7VrbvaRulVdXc2VV17Z5r5zzjmHIUOGUFNTA8Bbb73FE088wZtvvtlm3LJly+jbty/vvvsuP/zhD9m6dWvC77Fv3z7+4R/+gerqas477zxefvllAG666SY2b97Mzp07GTVqFD//+c+T1jp//nzuvfdeNm/ezIUXXnj2/oKCAlatWsW2bdvYsGED//iP/4iq8vDDD3PJJZewY8cOfvrTn3Y6znjrroqL6ZWb/L9kr9wc5lbYIbSmZ7ps/Kr6W+CTrsbFzAZe6FFFPeDlPlJVTXgkS/z9U6dO5fzzz+8w5ne/+x2zZs0CYMyYMYwdOzbh9xg+fDjjxo0D4Morr+TDDz8E4L333qOiooLLL7+c5557jurq6qS1/v73v2f27NkAfOtb32pT6z//8z8zduxYvva1r3H48GE++uijhNvkZJxJraH9Cll2+3j69Mrt8Mo/L0fo0yuXZbePtwMUTI+lbB+/iHyBlr8MXo67W4E3RGSriNydqu/VGS/3kY4ePbpDYujx48c5dOgQl1xyScu8nUSwOn21nJ+ff/bz3NxcotGWX1Df/va3eeyxx9i1axeLFi1ydPx8ol9Szz33HB9//DFbt25lx44dlJSUJJzL6TiTepNHDmDtggpmTxhCUX4eIlCUn8fsCUNYu6DCTt4yKZHKd4mmA79vt5vnK6p6REQGAOtE5P3YXxAdxH4x3A1QUlJCVVVVm8fPPfdcTpw4kbSA/zL6Al7e8WeiSXbz5+XAN8Zc0OVc7U2YMIFIJMKTTz7JnDlzaGpqYsGCBWc/P336NGfOnKGpqYkTJ04QiURobm7mxIkTXHXVVTz33HOUl5fz/vvvs2vXLurr6zlx4gSqSiQSaTMeoKGhgYaGBk6cOMHx48cpLi7mk08+4Re/+AUXXXQRJ06coKGhgV69enXYlokTJ/Jv//ZvzJo16+xuoRMnTvDRRx9x3nnncfr0ad544w0OHjxIJBKhqKiI48ePn52ns3Htv8/p06c7/Jyg5b2ERPeHSU/XYMp5MGVyPtD6YuAoB3YdJShH8NtzwN9rkMrGP4t2u3lU9Ujs3zoRWQVMABI2flVdDiwHKC8v1/Y51rt37+4yCmHelJGs2VVHtLmp0zG9cnO5929HUlzs/s/lV155hXnz5vHII4/Q3NzMDTfcwCOPPEJ+fj4FBQX07t2b3NxciouLKSoqIicnh+LiYhYuXMgdd9zBV77yFcrKyhg7diwDBw6kuLgYETl7uGXreGh59d/Y2EhxcTE//vGPmTJlCkOHDuXyyy8/GwuRn59Pfn5+h3V5/PHHmTNnDsuXL+fmm28GoLi4mDvvvJPp06czefJkxo0bx2WXXUZRURHDhg3jq1/9KldffTXXX3893/ve9xKOa/99CgoKKCsr67BOlsVuaxCE7ff6fB8/r4E42Q0hIsOA/1DVhGeNiMi5wAHgi6paH7uvEMhR1ROxz9cBD6rq2q6+X3l5ubbfrbJ7925GjRrVZa2Jsk6gZR9pr9wcz7NOEmX1NDU10djYSEFBAR988AFTpkxh79699O7d27M60qGzn4mfn/DpEvY18Pv2p6NPZOBCLFudHkTT5St+EXkBmAT0F5FaYBHQC0BVn4gNmwm80dr0Y0qAVbF9zXnA806afk+17iNdsfEAq7Yfpv5MlMLeecwsG8TciuEZeWPs5MmTTJ48mcbGRlSVn/3sZ4Fv+sYEVfz5Pu1Fm5Voc8v5PmsXVGTtG+ldNn5Vne1gzDO0HPYZf99+4IruFtYTQ/sV8lDlGN+c1l5cXBzYS0kak23cnO/jlx6SaoE6c9eOJfcP+1mYoPLyfJ+gCEzjLygo4NixY9ZwfKA1j7+goCDTpRjjmmUiBSidc/DgwdTW1vLxxx9nupSkTp8+HYqG2HoFLmOCpjA/j4iD5p/NmUiB2bJevXoF4mpPVVVVCQ9xNMb4Q2XZQFa+cyjp7p5sz0QKzK4eY4xJBctEssZvjAkZy0Syxm+MCaGwZyIFZh+/MaajTTVHWbSmmn11kbP3lQ4o4oEZo7nm0v5pq6N9/EF+Ltzy6S5fX+7Ub+f7pJO94jcmoB5dv5c5K95u0/QB9tVFmLPibR5dvzctdSS63OnpJlJyuVPjDWv8xgTQppqjLFm3L+mYJev2sanmqKd1eHm5U+Mda/zGBNCiNckvxtNq8avOxnWXl5c7Nd6xxm9MALXfvdOZvR85G9ddFn8QTNb4jTHdZvEHwWSN3xjTbV5e7tR4xxq/MQFUOqDI0bgRJc7GdVdl2cAOJ0G1l+3xB0Fkjd+YAHpgxmhH4xZPdzauuyz+IJis8RsTQNdc2p+FU0uTjlk4tdTzk7gs/iCYrPEbE1Dzp4zg+bkTO+zOGVFSxPNzJzJ/yoi01JEo/qBPHqGJPwgie8fFmAC75tL+vLHw2kyX0SH+oOVC4+GLQggKa/wmcNrnwhTm51FZNrDHuTBezetlHV7W3JM6usrq8cM2er12v6hu4L4Nr2f0udQZ6epShiLyNPANoE5VO/wKF5FJwCtA66l5v1LVB2OPTQMeBXKBFar6sJOiysvLNagXJ295pTMp02VklJdrsGFPHfOe3UZjU3ObE4fycoReuTksu318t3YtpHre7q6Bmzq8Wguv6/DDNqZj7c5Em2iKa6+pmDsZEdmqquVOxjrZx/8MMK2LMRtVdVzso7Xp5wKPA9cDXwJmi8iXnBRlTCJe5cL4JW/GTR1e1uxlHX7YxnStXVO719R+yi7qsvGr6m+BT7ox9wSgRlX3q+oZYCVwYzfmMQbwLhfGL3kzburwsmYv6/DDNvpl7TIpVUf1XC0iO0Xk1yLSeuDwIOBQ3Jja2H3GdItXuTB+yZtxU4eXNXtZhx+20S9rl0mpeHN3GzBUVSMicgOwGigFEp3O1+mKiMjdwN0AJSUlVFVVpaC09ItEIoGtPVW8WoOI01yYhqir7+/FvN1ZAzd1JG8tbcf6qQ4/bKNXzyOv506lHjd+VT0e9/lrIrJMRPrT8gr/i3FDBwNHksyzHFgOLW/uBvUNUntz17s1KNrwuqP/WIX5ea6+vxfzdmcN3NQBzpqM27Xwug4/bKNXzyOv506lHu/qEZELRURin0+IzXkM2AyUishwEekNzALW9PT7mfDyKhfGL3kzburwsmYv6/DDNvpl7TKpy8YvIi8AbwEjRaRWRO4UkXtE5J7YkFuA90RkJ/AvwCxtEQXuA14HdgMvqaq3V4UwWc2rXBi/5M24qcPLmr2sww/b6Je1yyQnR/XMVtWLVLWXqg5W1Z+r6hOq+kTs8cdUdbSqXqGqX1bVTXFf+5qqjlDVS1T1f3u5ISb7eZUL45e8GTd1eFmzl3X4YRvTtXa57V74+ym7qMsTuDLBTuAKNq/X4OCxelZsPMCq7YepPxOlsHceM8sGMbdieI/P5EzVvD1ZAzd1eLUWPa2jIBduKR/aaR1+2Eav1+7BF3/H2x+R8rk74+YELmv8KWaN39YAbA3Cvv2Q/jVw0/gtq8d4wi+5N6atTTVHWbSmus01e0sHFPHAjNGeRzing5cZQEHMDOqMxTKblHv34yjTlm5k5TuHiMSOx440RFn5ziGmLd3Ihj11mS4xlB5dv5c5K97ucKH2fXUR5qx4m0fX781QZamxYU+d4+edm7HdGe9Fzalkjd+k1MFj9Ty2oyHjuTemrU01R1mybl/SMUvW7WNTzdE0VZRafskX8qrmVLPGb1LqqY376SKqxBdZJWGzaI2zI6kXvxrMI679ki/kVc2pZo3fpNTq7Uc6pBK254eskrBpv3unM3s/cjbOb/ySL+RVzalmb+6alKp3mlVyxtk4Y5xw9bxzeCBj63PUq+d0Jv+v2Ct+k1KtGStdjuttrzlM6rh53rl9jnr1nM7k/xVr/CalKssGdjhjsT0/ZJWETemAoq4HQYcLtweFX/KFvKo51azxm5S6q+Jiuogq8UVWSdg8MGN014OAxdOdjfMbv+QLeVVzqlnjNyk1tF8h943Lz3jujWnrmkv7s3BqadIxC6eWBvYkLr/kC3lVc6pZ4zcpN/aCPNYuqGD2hCEU5echAkX5ecyeMIS1Cyo8udC06dr8KSN4fu7EDrtzRpQU8fzcicyfMiJDlaXG5JEDHD/v3Iztzngvak4le4fNeGJov0IeqhzDQ5VjMlqHm4iCVdtqWbSmmuOnPz+K4pyCPB6YMZqZ4wd3mDvZ+L4ebEsqDOrbhwnDz+fIp6fPxgNMGH4+g/r26TC2J/EH+blwy6e70h5/AKCth+1ou9vtuH2OevWczsT/FQtpSzELp/LPGjy6fm/Ss1UXTi09+yp3wYvbWb290wvEUVk2kKW3lZ293dX4qy/K4YX513ejau9s2FPHvGe30djU3Ob48bwcoVduDstuH3/2FaZXY7sz3ovtSwc/h7TZrh6TldxEFKzaVpu0iUPLyTarttUCOBr/1p+bz473A68iDSz+IJis8Zus5CaiwG2cQRDjD7yKNLD4g2Cyxm+ykpuIgvh99Ml8dqplnNvxfuBVpIHFHwSTvblrTAh4Fmlg8QeBZK/4jQkBryINLP4gmKzxm6zkJqLgnAJnzeDcPi3j3I73A68iDSz+IJi6bPwi8rSI1InIe508/nci8m7sY5OIXBH32IcisktEdohIMI/PNIHkJqLAbZxBEOMPvIo0sPiDYHLyiv8ZYFqSxw8A16rqWOAhYHm7xyer6jinx5cakwpuIgpmjh9MZdnApGMrywaePYnLyfirL8pJeNJXpngVaWDxB8Hk6AQuERkG/IeqJj21TET6Au+p6qDY7Q+BclV1dT03O4Er2Py0BptqjrL41eo2FxgZUVLE4umJz9xd/Gp1m6Nxzu2Tx+LpnZ+529n4vsdrfLMG8Q4eq2fFxgOs2n6Y+jNRCnvnMbNsEHMrhic8G7e7Ywty4ZbyoQnHup3bq+3zmp9P4Ep14/+fwGWqOjd2+wDwV1re+39SVdv/NRD/tXcDdwOUlJRcuXLlSif1+04kEqGoKJjRtqlia2BrEPbth/SvweTJkx03/pS9+yQik4E7ga/G3f0VVT0iIgOAdSLyvqr+NtHXx34pLIeWV/x+fLXkhJ9e7WaKn9bAy0yYZHMf2LW5zRq4yQzykpfrEc9Pz4FM8fMapOSoHhEZC6wAblTVY633q+qR2L91wCpgQiq+nzFObNhTx7SlG1n5ziEiDVEUiDREWfnOIaYt3ciGPXWezf3ux5/v/nl0/V7mrHi7w0ll++oizFnxNo+u39vtOlJZc0/WwwRLjxu/iAwBfgV8S1X3xt1fKCLFrZ8D1wEJjwwyJtW8zG5xMvdjOxo4eKzeVWaQlyzLxsRzcjjnC8BbwEgRqRWRO0XkHhG5JzbkR0A/YFm7wzZLgN+JyE7gHeA/VXWtB9tgTAdeZrc4mbupGVZsPOCbXB/LsjHxutzHr6qzu3h8LjA3wf37gSs6foUx3nOT3eI2B93J3E0Kq7YfJuIwSiD+qCMveLkeJnjszF2TlbzMbgliLkwQazbescZvspKX2S1BzIUJYs3GO9b4TVbyMrvFydy5AjPLBrnKDPKSZdmYeNb4TVbyMrvFydy5OTC3Yrhvcn0sy8bEs8ZvspKX2S1O5r5vXD5D+xW6ygzykmXZmHjW+E3WmjxyAGsXVDB7whCK8vMQgaL8PGZPGMLaBRU9uvB2V3OPveDzfeXzp4zg+bkTO+zOGVFSxPNzJ5694LvXvFwPEyz2To7JakP7FfJQ5RhPDlFMNnf7o+GvubQ/byy8NuU1uOXlepjgsMZvHEtXzksquak5iNtnTHdY4zeObNhTx7xnt9HY1Hz2RKDWnJeXtx5m2e3jfberwE3NQdw+Y7rL9vGbLgUx58VNzUHcPmN6whq/6VIQc17c1BzE7TOmJ6zxmy65yXnxCzc1B3H7jOkJ28dvuhTEnBdXNXd9EbrPxxqTBewVv+lSEHNe3NQcxO0zpies8ZsuBTHnxU3NQdw+Y3rCGr/pUhBzXtzUHMTtM6YnrPGbLgUx58VNzUHcPmN6whq/cSSIOS9uag7i9hnTXfZuVZZxGzuwqeYoi9ZUs6/u80v/lQ4o4oEZozskRgY150VbD9vRdrfbcbt9ydbazdie/CVhMROmO6zxZxG3sQOPrt/LknX7Osyzry7CnBVvs3BqadqSI73gZQxDV3PfOzaPSR7XYTETprsc7eoRkadFpE5E3uvkcRGRfxGRGhF5V0TGxz12h4jsi33ckarCTVtuYwc21RxN2PTjLVm3j001Rz2r2UtexjA4mfuxHQ2exkFYzITpCaf7+J8BpiV5/HqgNPZxN/AzABE5H1gETAQmAItEpG93izWdcxs7sGhNtaN5F7/qbJzfeBnD4GTupmY8jYOwmAnTE44av6r+FvgkyZAbgV9oiz8A54nIRcDXgXWq+omq/hVYR/JfIKab3MYOxO/TT2bvR87G+Y2XMQxO5m5SPI2DsJgJ0xOp2sc/CDgUd7s2dl9n93cgInfT8tcCJSUlVFVVpai09IpEIhmpPeI0oqAh6ro+t+MztQZtavBwPdzM7TANwnUdXm5fKvjhOZBpfl6DVDX+RKc9apL7O96puhxYDlBeXq6TJk1KUWnpVVVVRSZqL9rwuqNmUJif11Lf2v90PLfb7cnUGsRzvR4ezQ3OmrTbOrzcvlTww3Mg0/y8Bqk6jr8W+GLc7cHAkST3mxRzGztQOqAo6dhW7a8TGxRexjA4mTtX8DQOwmImTE+kqvGvAf5r7OieLwOfqeqfgdeB60Skb+xN3eti95kUcxs78MCM0Y7mXTzd2Ti/8TKGwcncuTl4GgdhMROmJ5wezvkC8BYwUkRqReROEblHRO6JDXkN2A/UAE8B8wBU9RPgIWBz7OPB2H0mxdzGDlxzaX8WTi1NOufCqaUdTuIKCi9jGJzMfd+4fE/jICxmwvSEqDp9+yl9ysvLdcuWLZkuo1syvV/v4LF6Vmw8wKrth6k/E6Wwdx4zywYxt2J4p2fuLn61us3ROyNKilg8veOZu05leg3iuV2PVM19YNfmNmvgVR1ebl9P+Ok5kCnpXgMR2aqq5Y7GWuNPLXvC2xqArUHYtx/83fgtpM0YY0LGGr8xxoSMNX5jjAkZa/zGGBMy1viNMSZkrPEbY0zIWOM3xpiQscZvjDEhY43fGGNCxhq/McaEjDV+Y4wJGWv8xhgTMtb4jTEmZKzxG2NMyFjjN8aYkLHGb4wxIWON3xhjQsYavzHGhIw1fmOMCRlHjV9EponIHhGpEZH7Ezy+RER2xD72isincY81xT22JpXFG2OMcS+vqwEikgs8DkwFaoHNIrJGVf/YOkZVF8aN/+9AWdwUp1R1XOpKNsYY0xNOXvFPAGpUdb+qngFWAjcmGT8beCEVxRljjEm9Ll/xA4OAQ3G3a4GJiQaKyFBgOPBm3N0FIrIFiAIPq+rqTr72buBugJKSEqqqqhyU5j+RSCSwtaeKrYGtQdi3H/y9Bk4avyS4TzsZOwv4pao2xd03RFWPiMjFwJsisktVP+gwoepyYDlAeXm5Tpo0yUFp/lNVVUVQa08VWwNbg7BvP/h7DZzs6qkFvhh3ezBwpJOxs2i3m0dVj8T+3Q9U0Xb/vzHGmDRz0vg3A6UiMlxEetPS3DscnSMiI4G+wFtx9/UVkfzY5/2BrwB/bP+1xhhj0qfLXT2qGhWR+4DXgVzgaVWtFpEHgS2q2vpLYDawUlXjdwONAp4UkWZafsk8HH80kDHGmPRzso8fVX0NeK3dfT9qd3txgq/bBFzeg/qMMcakmJ25a4wxIWON3xhjQsYavzHGhIw1fmOMCRlr/MYYEzLW+I0xJmSs8RtjTMg4Oo7fZK+Dx+p5auN+Vm8/Qn1DlML8PCrLBnJXxcUM7VeY6fKMMR6wxh9iG/bUMe/ZbTQ2NRNtbjnhOtIQZeU7h3h562GW3T6eySMHZLhKY0yq2a6ekDp4rJ55z27jVGPT2abfKtqsnGpsYt6z2zh4rD5DFRpjvGKNP6Se2rifxqbmpGMam5pZsfFAmioyxqSLNf6QWr39SIdX+u1Fm5VV2w+nqSJjTLpY4w+p+oaos3FnnI0zxgSHNf6QKsx39r5+YW97/9+YbGONP6QqywaSl5Poqpqfy8sRZpYNSlNFxph0scYfUndVXEyv3OQ//l65OcytGJ6miowx6WKNP6SG9itk2e3j6dMrt8Mr/7wcoU+vXJbdPt5O4jImC1njD7HJIwewdkEFsycMoSg/DxEoys9j9oQhrF1QYSdvGZOl7J27kBvar5CHKsfwUOWYLsdavIMx2cEav3HE4h2MyR6OdvWIyDQR2SMiNSJyf4LHvy0iH4vIjtjH3LjH7hCRfbGPO1JZvEkPi3cwJrt02fhFJBd4HLge+BIwW0S+lGDoi6o6LvaxIva15wOLgInABGCRiPRNWfUmLSzewZjs4uQV/wSgRlX3q+oZYCVwo8P5vw6sU9VPVPWvwDpgWvdKNZli8Q7GZBcn+/gHAYfibtfS8gq+vZtF5G+AvcBCVT3UydcmPCNIRO4G7gYoKSmhqqrKQWn+E4lEAlt7ZyJO4x0aolRVVWXlGrgV9jUI+/aDv9fASeNPdHpn+5d/rwIvqGqDiNwD/Dvwtw6/tuVO1eXAcoDy8nKdNGmSg9L8p6qqiqDW3pmiDa87av6F+XlMmjQpK9fArbCvQdi3H/y9Bk529dQCX4y7PRg4Ej9AVY+pakPs5lPAlU6/1vifxTsYk12cNP7NQKmIDBeR3sAsYE38ABG5KO7mDGB37PPXgetEpG/sTd3rYveZALF4B2OyS5eNX1WjwH20NOzdwEuqWi0iD4rIjNiw74hItYjsBL4DfDv2tZ8AD9Hyy2Mz8GDsPhMgFu9gTHZxdAKXqr4GvNbuvh/Fff594PudfO3TwNM9qNH4QGu8w4qNB1i1/TD1Z6IU9s5jZtkg5lYMt6ZvTIDYmbvGMTfxDsYY/7LGH3Ju8nc21Rxl0Zpq9tVFzt5XOqCIB2aM5ppL+6et5lXbalm0pprjpz8/0uicgjwemDGameMHp62OZGtnjJ9ZOmeIbdhTx7SlG1n5ziEiDVGUz/N3pi3dyIY9dWfHPrp+L3NWvN2m6QPsq4swZ8XbPLp+b1pqXvDidha+tLNN0wc4fjrKwpd2suDF7Wmpo6u1e/dju2Sl8S9r/CHlJn9nU81Rlqzbl3S+Jev2sanmqJcls2pbLau3Jz8aePX2I6zaVutpHU7W7rEdDZZdZHzLGn9IucnfWbSm2tGci191Nq67/FKHk7Vrasayi4xvWeMPKTf5O+1373Rm70fOxnVX+907nfnslLe7WZysXZNi2UXGt6zxh1S90/ydM7avuj1bOxN01vhDqjDf2QFdhb3twK/2bO1M0FnjDyk3+TulA4oczTmixNm47jqnwFkjPbePtw3XydrlCpZdZHzLGn9IucnfeWDGaEdzLp7ubFx3+aUOJ2uXm4NlFxnfssYfUm7yd665tD8Lp5YmnW/h1FLPT+KaOX4wlWUDk46pLBvo+UlcTtbuvnH5FmNhfMsaf4i15u/MnjCEovw8RKAoP494ShX5AAAK/ElEQVTZE4awdkFFm4unz58ygufnTuywO2dESRHPz53I/Ckj0lLz0tvKWHLrFR1255zbJ48lt17B0tvK0lJHV2s39gLbv2/8y56dIecmf+eaS/vzxsJr01BVcjPHD05rNENnkq2dHcFv/MwafxJucmxax/+iuoH7NrzuaLxXdXg1t1+yevzCsnpMUIlq8hNRMqG8vFy3bNmS0Ro27Klj3rPbaGxqbnOyTl6O0Cs3h2W3j2+zK6R1/JloE01xS9rZeK/q8GruR9fvTRrbsHBq6dndPX6+5FyqdLV2947N4zvf/FoGK8ysMDwHupLuNRCRrapa7mSs7eNPwE2OTfvxTe1+jyYa71UdXs3tl6wev7CsHhN01vgTcJNj053xXtXh1dx+ycjxC8vqMUFnjT8BNzk23RnvVR1eze2XrB6/sKweE3TW+BNwm8XiVXaLl5kwljfTfbZ2Juis8SfgNovFq+wWLzNhLG+m+2ztTNA5avwiMk1E9ohIjYjcn+Dx/yEifxSRd0VkvYgMjXusSUR2xD7WpLJ4r7jJsenOeK/q8Gpuv2T1+IVl9Zig67Lxi0gu8DhwPfAlYLaIfKndsO1AuaqOBX4J/CTusVOqOi72MSNFdXvKTY5Nd8Z7VYdXc/slI8cvLKvHBJ2TV/wTgBpV3a+qZ4CVwI3xA1R1g6qejN38A5D50yp7wE2OTfvxue1eCCYa71UdXs3tl6wev7CsHhN0XZ7AJSK3ANNUdW7s9reAiap6XyfjHwP+oqo/jt2OAjuAKPCwqq7u5OvuBu4GKCkpuXLlypXd26IUqjvZzOsfNrLpSJTTUSjIg2sG5vH1Yb0Y8IWOvzPrTjbz6t6TbDkqjsZ7VYdXc+8+FuXZ3Wc4HPn8OTOoSLh9VG9G9ft8f3YkEqGoKPt3+yRbuy80nwzFGnQmLM+BZNK9BpMnT3Z8ApeTd58S7cxM+NtCRG4HyoH4QJchqnpERC4G3hSRXar6QYcJVZcDy6HlzF2vznhbta2WRWuq21zG75yCPB6YMTph/sutLufetmknrVf+OxWFzXXwjS+PYFK7ud3EH6zaVsvmd6odzes23uEnv97Nm3/af/YHeioKb/4pysjhQ7h10qg2Y/+6rZYTu6pp+R3e4kQ0lwuHX9amjqqqKnoPHuMq3sGr6AivYyY6e360P2vTy9gNP7Izd/29Bk5e8V8NLFbVr8dufx9AVf9Pu3FfA/4VuFZV6zqZ6xngP1T1l8m+p1eRDQte3M7q7Uc6fbyybGC30x3dzO0m/sDNvG7jHaYt/Q3v/6XzY+8vu7CItQuudV3HwhWvs6qm80MZ47fPbd1u1s7N2FSL/0/vZeyGX/m56aVL0CMbNgOlIjJcRHoDs4A2R+eISBnwJDAjvumLSF8RyY993h/4CvBHZ5uRWqu21SZtXNByYs6qbbWezu0m/sDNvG7jHX7y691Jmz7A+3+J8JNf73a9fcmafvz2gXfREX6JmfAydsOY7upyV4+qRkXkPuB1IBd4WlWrReRBYIuqrgF+ChQB/09EAP4UO4JnFPCkiDTT8kvmYVXNSON3EzvgNvLXzdwDigscj/3LZ6cdj51+xUDHEQwPVY7hZ7/Z72juJ36znyKHlzx0u31vLLzWVXTEH/Yfczy30+zB1jq84mb7nERjG5MKjv5Hq+prwGvt7vtR3OcJYwhVdRNweU8KTJX4ffrJfHbK/dmWbub+7FTq4w8+OxV1FcHwUOWYxG/SJNCMt9vnpu6IwzNm3ayd1zETbn8uxqSDnbmbJYIaIxDUup3K9u0zwWSNP0sENUYgqHU7le3bZ4IpNI3/HIf7qdtfyzXVc7uJP3Azr9t4h+QjP5eDd9sH3kVH+CVmwsvYDWO6KzSN38vYATdzezXWbbzDvdc6uzzgPdde7FnN4F10hF9iJryM3TCmu0LT+GeOH0xl2cCkYyrLBnbrIt5u5nYTf+BmXrfxDv90/SguuzD5q93LLizin64f5Xr7Zl6a/C+E+HgHr6Ij/BIz4WXshjHdFbpr7q7aVsviV6vbHL1zbp88Fk9PfOau27n/1692Uh/3Pl1nc2+qOcriV6vbHFUyoqSIxdMTn7nrtOaDx+pZsfEAq7Yfpv5MlMLeecwsG8TciuGdnrn7xG/2E3/AYQ4tr/T/6fq2Z+46raP1zF2n2+e2bjdr52ZsKiU6c9fNzyXo7AQuf5/AFbrG7zV7wtsagK1B2Lcf/N34s+JQgmzPQXGzfW6ziPxQszEmvQK/j3/DnjqmLd3IyncOEWmIokCkIcrKdw4xbelGNuxJGBsUGG62b8GL21n40s4OJ1wdPx1l4Us7WfDi9rTU/O7H0az+mRgTdIFu/Nmeg+Jm+7zMInLj4LF6HtvRkLU/E2OyQaAbv5sclCBys31u8oK89NTG/XRRcqB/JsZkg0A3fjc5KEHkZvu8zCJyY/X2IzR1cbxAkH8mxmSDQDf+bM9BCeL2BbFmY8Im0I0/23NQgrh9QazZmLAJdOPP9hwUN9vnZRaRG5VlAztccL69IP9MjMkGgW782Z6D4lWOjZfuqriYLkoO9M/EmGwQ6Maf7TkobrbPyywiN4b2K+S+cflZ+zMxJhsEuvEDTB45gLULKpg9YQhF+XmIQFF+HrMnDGHtgorAX8TazfYtva2MJbde0WF3zrl98lhy6xXdvpC8W2MvyMvqn4kxQZcV77AN7VfIQ5VjsvbSdW62b+b4wWmNZuhMtv9MjAmywL/iN8YY4441fmOMCRlr/MYYEzK+zOMXkY+Bg5muo5v6A0czXUSG2RrYGoR9+yH9azBUVS9wMtCXjT/IRGSL04shZCtbA1uDsG8/+HsNbFePMcaEjDV+Y4wJGWv8qbc80wX4gK2BrUHYtx98vAa2j98YY0LGXvEbY0zIWONPIRFZKCLVIvKeiLwgIgWZrimdRGR+bNurRWRBputJBxF5WkTqROS9uPvOF5F1IrIv9m/fTNbotU7W4Jux50GziPjyyJZU6mQNfioi74vIuyKySkTOy2SN8azxp4iIDAK+A5Sr6hggF5iV2arSR0TGAHcBE4ArgG+ISGlmq0qLZ4Bp7e67H1ivqqXA+tjtbPYMHdfgPeAm4LdpryYznqHjGqwDxqjqWGAv8P10F9UZa/yplQf0EZE84AvAkQzXk06jgD+o6klVjQK/AWZmuCbPqepvgU/a3X0j8O+xz/8dqExrUWmWaA1Udbeq7slQSWnXyRq8Efu/APAHIPPpiTHW+FNEVQ8DjwB/Av4MfKaqb2S2qrR6D/gbEeknIl8AbgC+mOGaMqVEVf8MEPvXcqjNfwN+nekiWlnjT5HYftwbgeHAQKBQRG7PbFXpo6q7gf9Ly5+3a4GdgF1R3YSeiPyAlv8Lz2W6llbW+FPna8ABVf1YVRuBXwHXZLimtFLVn6vqeFX9G1r+7N2X6Zoy5CMRuQgg9m9dhusxGSIidwDfAP5OfXTsvDX+1PkT8GUR+YKICDAF2J3hmtJKRAbE/h1Cyxt7L2S2ooxZA9wR+/wO4JUM1mIyRESmAd8DZqjqyUzXE89O4EohEXkAuI2WP+u2A3NVtSGzVaWPiGwE+gGNwP9Q1fUZLslzIvICMImWJMaPgEXAauAlYAgtLwi+qart3wDOGp2swSfAvwIXAJ8CO1T165mq0WudrMH3gXzgWGzYH1T1nowU2I41fmOMCRnb1WOMMSFjjd8YY0LGGr8xxoSMNX5jjAkZa/zGGBMy1viNMSZkrPEbY0zIWOM3xpiQ+f958P9elKasDQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(beta[1]*xdata, ydata, 'o', label='Original data', markersize=10)\n",
    "plt.grid(True)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Output the comma-separated list of coefficients before $x$. You will need this for Google Form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 342,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6968201364375617,1.7433976718964401,-30.004765706028103\n"
     ]
    }
   ],
   "source": [
    "coefs = []\n",
    "\n",
    "# assuming that beta[1] is always a coefficient before x\n",
    "for m in range(1, 4):\n",
    "    coefs.append(str(linear_least_squares(m, xdata, ydata)[1]))\n",
    "    \n",
    "print(\",\".join(coefs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Output the comma-separated list of average quadratic deviations. You will need this for Google Form.\n",
    "\n",
    "*Reminder*: Average quadratic deviation = $\\cfrac{1}{n} \\sum\\limits_{i = 1}^{n} \\left| y_i - f(x_i, \\vec \\beta) \\right|^2 $, where $n$ is the number of observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 299,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qd(x):\n",
    "    return np.sum(np.abs(np.power(x, 2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 337,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.11953278794286401,0.11769389348016475,0.1051873346328916\n"
     ]
    }
   ],
   "source": [
    "aqds = []\n",
    "for m in range(1, 4):\n",
    "    beta = linear_least_squares(m, xdata, ydata)\n",
    "    cur_aqd = qd(ydata - beta@[np.power(xdata, i) for i in range(m+1)])/len(xdata)\n",
    "    aqds.append(str(cur_aqd))\n",
    "    \n",
    "print(\",\".join(aqds))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}