{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "#### EECS 16A Discussion 14B" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Problem 1: Polynomial Fitting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this discussion, we are trying to fit data (observations) of the form $\\{(x_i,y_i),i=1,2,...,n\\}$ to a polynomial that we know looks like this:\n", "\n", "$$y = f(x) = a_0 + a_1x + a_2x^2 + a_3x^3 + a_4x^4$$\n", "\n", "In other words, we want to find $a_0$, $a_1$, $a_2$, $a_3$, and $a_4$ that best fit the data.\n", "\n", "More generally, we might want to fit the data to a polynomial of different degree (for instance, if we do not know that the polynomial looks like as above), so we could try to solve for some $a_0,a_1,\\ldots,a_d$ that define a $d$-degree polynomial.\n", "\n", "Note that the observations are not perfect -- they are *noisy*, which means that $y_i \\neq f(x_i)$ in general! That is what makes this problem interesting.\n", "\n", "This first block of code contains functions that will help us set up some useful objects -- the polynomial curve for a set of parameters $a_0$, $a_1$, $a_2$, $a_3$, and $a_4$ and a \"data\" matrix that we will use to compute the least squares estimate later." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initialization" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import numpy.matlib\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import Axes3D\n", "from numpy import pi, cos, exp\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Construct Data Matrix\n", "Use this block to make the data matrix. The input data is given in the array `x_data`, and the output data is given in the array `y_data`.\n", "\n", "To construct the data matrix, which we call `DataMat`, you only need to specify what degree polynomial you will use to fit the data.\n", "\n", "If `x_data` has the form `x_{data}`$ =\\begin{bmatrix}x_1\\\\x_2\\\\ \\vdots \\\\x_n \\end{bmatrix}$, then `DataMat` has the form `DataMat` $=\\begin{bmatrix}1 &x_1&x_1^2& \\cdots &x_1^d\\\\1 &x_2&x_2^2& \\cdots &x_2^d\\\\ \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\1 &x_n&x_n^2& \\cdots &x_n^d \\end{bmatrix}$, where $d$ is the degree of the polynomial. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Function that defines a data matrix for some input data" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "\"\"\"You'll understand why it's constructed like this after\n", "doing the worksheet!\"\"\"\n", "def data_matrix(input_data,degree): \n", " # degree is the degree of the polynomial you plan to fit the data with \n", " Data=np.zeros((len(input_data),degree+1))\n", " \n", " for k in range(0,degree+1):\n", " Data[:,k]=(list(map(lambda x:x**k ,input_data)))\n", " \n", " return Data\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Function that helps with plotting " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "\"\"\"Function that defines the polynomial curve for a set of parameters and a range. The set of parameters defines the degree of the\n", "polynomial.\"\"\"\n", "def poly_curve(params,x_input):\n", " # params contains the coefficients that multiply the polynomial terms, in degree of lowest degree to highest degree\n", " degree=len(params)-1\n", " x_range=[np.min(x_input), np.max(x_input)]\n", " x=np.linspace(x_range[0],x_range[1],1000)\n", " y=x*0\n", " \n", " for k in range(0,degree+1):\n", " coeff=params[k]\n", " y=y+list(map(lambda z:coeff*z**k,x)) \n", " return x,y\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Printing the data matrix" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00]\n", " [1.000000e+00 5.000000e-01 2.500000e-01 1.250000e-01 6.250000e-02]\n", " [1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00]\n", " [1.000000e+00 1.500000e+00 2.250000e+00 3.375000e+00 5.062500e+00]\n", " [1.000000e+00 2.000000e+00 4.000000e+00 8.000000e+00 1.600000e+01]\n", " [1.000000e+00 2.500000e+00 6.250000e+00 1.562500e+01 3.906250e+01]\n", " [1.000000e+00 3.000000e+00 9.000000e+00 2.700000e+01 8.100000e+01]\n", " [1.000000e+00 3.500000e+00 1.225000e+01 4.287500e+01 1.500625e+02]\n", " [1.000000e+00 4.000000e+00 1.600000e+01 6.400000e+01 2.560000e+02]\n", " [1.000000e+00 4.500000e+00 2.025000e+01 9.112500e+01 4.100625e+02]]\n" ] } ], "source": [ "x_data = [0.0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5]\n", "y_data = [24.0, 6.61, 0, -0.95, 0.07, 0.73, -0.12, -0.83, -0.04, 6.42]\n", "\n", "D = data_matrix(x_data, 4)\n", "print(D)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Least Squares\n", "Recall the system of equations we are trying to solve here:\n", "\n", "$$\\texttt{D}\\vec{a} = \\vec{y},$$ where $\\vec{a} = \\begin{bmatrix}a_0\\\\a_1\\\\:\\\\a_d\\end{bmatrix}$ and $\\vec{y} = \\begin{bmatrix}y_1\\\\y_2\\\\ : \\\\y_n \\end{bmatrix}$.\n", "\n", "In the next block, you will implement code to compute $\\vec{a}$, the set of optimal polynomial coefficients (optimal in a least squares sense) to fit the data. Put the coefficients in an array called params." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Finding the least squares solution" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 24.00958042 -49.99515152 35.0039627 -9.99561772 0.99841492]\n" ] } ], "source": [ "# params = np.linalg.solve(np.dot(np.transpose(D), D), np.dot(np.transpose(D), y_data))\n", "params = np.linalg.inv(np.transpose(D) @ D) @ np.transpose(D) @ y_data\n", "print(params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot Curve Fit\n", "Use the next block to compare your fitted polynomial to the data. All you need to do is enter the polynomial coefficients (in degree of lowest degree to highest degree) into an array called params. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the fitted polynomial and the original datapoints" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 1.0, 'Polynomial of Degree 4')" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEKCAYAAAAGvn7fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsw0lEQVR4nO3deXxU5dn/8c+VEBLCoogRomChsq9BogRRlIIW0F8VKxVc6lrrjra4PxV9tFYLijx9rNZqRVGsuJbHlUVcQKKGHWQJSpQgyCL7Tub6/XGfwBATMklm5sycXO/Xa16TmTlnznVOku/cc5/7nCOqijHGmOBI8bsAY4wx0WXBbowxAWPBbowxAWPBbowxAWPBbowxAWPBbowxAWPBnuBE5HIR0bDbNhGZLyI3ikidKr7XRyLyUYxKjQsRaelth8urMW9U119E/p+ILBSR3V5NR1Yw3X3l/A4LRWSCiPwyWvUkAxH5uYjs9LZDa7/rCaoqBYPx1RCgGGjk/fw34BjgXj+L8sEaoBfwtZ9FeB+qLwGfATcAe4Ftlcx2KlACZAKtgAuA90XkReAyVQ3FruKE8XdgC1DP70KCzII9ecxT1RXez5O91s5walmwq+oeIN/vOoDjgIbARFX9JMJ5PlfV/WGPnxWRW4HHgHnAo9EtsWIikgbs1zgeoSgiFwHdgb8AY+K13NrIumKS15dAIxE5BkBEBojILBHZJSJbROQtEWlX0cwi0kxE9orI8HJeu8/7utzYe/yRiMwQkf4iMsd7bZGIDC5n3krrCHu/ASIyz5t2roj0FJE6IvKQiKwRkR9FZJyI1A+b9yddMSJykoi8JiLF3nst896jWq1CEckWkRdEZIOI7BGRBSJySfj2AYq8h8969XxUnWWp6hhgLu5DOryGLBF5SkRWezUsFZFryqm1v7ftdovIChG52ttmRWHTlG6z60XkryLyPbAHONJ7/XwRyfd+r5tF5FUROb6cZV3jdQPu9rbNsyJyVCTr6f0tPQaMADZHun1M9ViwJ69WuK/120VkAPAOsB24ELgO6AzMEJHjyptZVdcCbwGHhIWIpAJX4Vqim8JeOgEYi/vnPB/XJfJqeD9pFetoDYwCHsZ1LaUDk4AngWzgcuC/gYuBkZVsi+NxLd5rgQFenVcCz1Uy3094HyIfAwOBu4HzgIXA+LBgfcarGeBBXNfQ9VVdVpj3gBalYSoijYAZwCDgPuBs4P+AJ0XkprBaO3Jwew/16h0O/KKC5dwDtMX9zgcDu0XkWuB14Ctc19Dvcb+zj0WkYdiyHgaeAKYCvwJuw23r97y/mcr8FViqquMjmNbUlKraLYFvuIBToB2u66wx7p+vBHjLm6YAKATqhM3XCtgHPBb23EfAR2GPz/De+7Sw537lPZdXZr59QJuw547xarg77Lmq1LEP+Hk5y51aZv3fAFaGPW7pTXd5BdtLvO10CRACmlS0/hXMf6P3/meUeX4qsA5I9R63PlwdZea9z5u2TgWv/957vaf3+E/A7vDt7T3/T2BD6fsAE4D1QGbYNNnevEXlbLM5gIQ93wDX3/2vMstphdtncEvY/CXAvWWm6+2973mVrP9puG8IHcv8Tbf2+/8rqDdrsSePpbgw/BG3A+ol4EqvhXki8IqG9d+q6kpgJnB6RW+oqh/hWmq/D3v698ACVS3bj12oqoVh867DBV1pK7OqdSxX1W/KrB/AB+Wsd3MRkYrWQ0QaicgjIvI1LkD2AeNxId+movkq0AdY7W2bcC8CWUDHKr5fJErXrbS/ewDwObDS65qq4+2s/QBoElZDHvCuqu4sfSNVXYPboVuet9RLVk8v3M74l8osZxVuu/fxpjsT9+2+7HSf43YY96ECIlIX+AcwRlW/qnRLmKiwnafJYzBuVMw24FtV3Q0gIs1xwbCmnHnWAj+r5H2fBEZ7fe0NcKFyYznT/VjOc3uADO/nxlWsY1OZx3sP83wdIBXYT/meA/rjdiTPA3YAJ+O6DjIqmKciR1HxOpS+Hm0tvPvS5R6D+0awr4Lpm3j32bgP17J+AH5ezvNl1+sY735qBcvZVGa6FRVM16SC5wFuwf1t/I8cHA6a6d03FJGGqlrZaCJTRRbsyWORHhwVE24TrqXXrJzXmlF+IId7ATdK4XLcP+BO3LeBqqppHdUiIhnAucB9qjo27Pku1XzLH3HdXmU1C3s92gYB36nqKu/xRlxg/2THtmeZd7+Gg6EbrmkF85UdAbPRu78cWFzO9NvKTHcWP/3gDX+9PB1x2251Oa/NAeYDOYeZ31SDBXuSU9UdIjIbGCIi96lqCYCI/Aw4BTfe/XDzbxWRl3BdMA2Al1V1a7zrqIF0XGu+bOv28mq+38e4deitqjPDnr8IF7ZR7U4QN9wxB9eyLfU+cBMu7MtrkZfKBwaJSGZpd4yIZOP6vsv71lHWZ7jwbq2qzx9muim4/RXHq+qUCN433MPAuDLPDQDuwO0HWVZ2BlNzFuzB8Cfc6Ii3ReTvuIC+H7djLJKx0X/nYD/7Uz7WUWWqukVE8oE/isga3M7FK3HjzKtjHK6l/IaI3IPr/roY18/8+9IPrGrqKSIluO6hn+NGoQwEngf+J2y6MbhRRZ+KyBhc+NUH2uN2dJ/rTfeg9x4fiMho3Ifcn3BdMZUe7OR9qN8GPCEiWbjROVtw2+503I7mCar6tYg8AvyvN3T1Y9wO2ha47fKMqk6vYBlLObj/BHDDL70fP6/gW6ipIQv2AFDV90XkbNywwIm4fumPgNtV9fsI5l8gIsuBrao6x686amAYbl/BE8Aub9nDgber+kbeN4/TccPzHsYdhLQMuFRVX6xhnTO8+524FvUXwABVPWSHsfdhdQpun8EduKDd7NXxeth0X3nbexRunVcDj+BaxC0jKUhV/yEiq3DDFy/CZcJq4FPc/orS6e4WkSW4o2xvwHXrrAKm4UZCmQQih+4kN7WR1wpbAvxOVZ/1ux5TfSLSALeT8x1Vvcrveow/LNhrMW9ETWtcd0lrXF/rLn+rMlUhIn/D9ZV/DxyL+6bSHThJVRf4WZvxT8Tj2EWkhYhMF5GvRGSxNzyu9PDz1eIODZ8nIoNiV66JsquBD3GjKC6yUE9KGbjul8nA07ihnv0t1Gu3iFvs3t72bFWd4x1qPBt3uPVvgO2qOjpmVRpjjIlYxDtPvSPa1ng/b/N2pFR35IExxpgYqVYfuzdc6RPcyYL+gBszvBV3rpA/6qEnjyqd5xq8E07Vr1+/R/v27atdtDHG1EazZ8/eoKpZlU1X5WD39rp/DPxZVd8Qkaa4scMKPIDrrrnycO+Rm5urBQUFVVquMcbUdiIyW1VzK5uuSicBE3dy/teBl1T1DQBV/UFVS9Rd/eWfuHN0GGOM8UlVRsUI8CywRFUfC3s+O2yywcCi6JVnjDGmqqpy5Glv4FJgoYjM8567GxgmIjm4rpgiDj0FbK2Snw9jx8Ly5dC2LQwfDnl5fldljKltqjIqZgYHzxsd7t3olZO8Ro6E0aNh1y5QhXnzYNIkGDEC7r/f7+qMqb59+/ZRXFzM7t27/S6l1sjIyKB58+akpaVVa347V0wU5Oe7UN+9s4QuLGITjVkVOp6dO93zAwday90kr+LiYho2bEjLli05zPVOTJSoKhs3bqS4uJhWrVpV6z3sCkpRMHasa6mnsY85nMhVHDzdyu7d7nVjktXu3btp0qSJhXqciAhNmjSp0TckC/YoWL7cdb/sIYOVtKIDSw68FgpBoZ37ziQ5C/X4qun2tmCPgrZtIcXbkkvocEiwp6S4140xJl4s2KNg+HDI8K6suYQOtGU5qd7lOTMy4OabfSzOmDjLz4dhw6BHD3efX/ay6NWQmppKTk4OnTp1olu3bjz66KOEQoe/lkhRURETJkyo+cIrcfXVV/PVV4e/sNZbb71V6TTRZMEeBXl5bvRLZiYskw6ks5cTZCWZme5523FqaouRI6FfP3jlFZgzByZOdI9HjqzZ+9arV4958+axePFipkyZwnvvvcf9lQw3i1ewP/PMM3Ts2PGw08Q72FHVuN969OihQTRrlup/9Z+lCjrqtP/orFl+V2RMzX311VcRTTdrlmpmpqrb43ToLTNTa/T/UL9+/UMef/3113rUUUdpKBTSlStX6qmnnqrdu3fX7t2768yZM1VVtWfPntqoUSPt1q2bPvbYYxVOF27lypXarl07veiii7R9+/b661//Wnfs2KGqqlOnTtWcnBzt3LmzXnHFFbp7925VVT399NP1yy+/PFDn3XffrV27dtWePXvq2rVrdebMmdq4cWNt2bKlduvWTVesWKFjx47VDh06aJcuXfTCCy8sd53L2+5AgUaQsRbs0bZ5s9usDz/sdyXGREWkwT50qKpI+cGekuJer66ywa6qesQRR+jatWt1x44dumvXLlVVXb58uZbmy/Tp0/Xss88+MH1F04VbuXKlAjpjxgxVVb3iiit01KhRumvXLm3evLkuW7ZMVVUvvfRSHTNmjKoeGuyATpo0SVVVb7vtNn3ggQdUVfWyyy7TV1999cBysrOzD3wwbNq0qdx1rkmwW1dMtB1xBGRnw9KllU9rTICUjg4rTyxHh+3bt4/f/e53dOnShSFDhlTY5RHpdC1atKB3794AXHLJJcyYMYNly5bRqlUr2nojIS677DI++eSTn8xbt25dzjnnHAB69OhBUVFRucvo2rUrF198MS+++CJ16kT/cCIL9ljo0AGWLKl8OmMCJHx0WFnRHh32zTffkJqayjHHHMOYMWNo2rQp8+fPp6CggL1795Y7T6TTlR1qWJWhh2lpaQemT01NZf/+/eVO984773DDDTcwZ84cTjrppAqnqy4L9lgoDXa7nqypRcJHh5UVzdFh69ev59prr+XGG29ERNiyZQvZ2dmkpKQwfvx4SkpKAGjYsCHbtm07MF9F05X13XffMWvWLAAmTJjAqaeeSrt27SgqKmLFihUAjB8/ntNPPz3imsNrCYVCrFq1ir59+/LII4+wZcsWtm/fXq1tUREL9ljo0AG2boU1a/yuxJi4CR8dVtpyT0khKqPDdu3adWC4Y//+/TnrrLMY6Q21uf7663n++efp1q0bS5cupX79+oDr7khNTaVbt26MGTOmwunKateuHU888QQdOnRg06ZNXHfddWRkZPDcc88xZMgQunTpQkpKCtdee23E9Q8dOpRRo0bRvXt3CgsLueSSS+jSpQvdu3fn5ptv5sgjj6z+xilHta6gVFOBv9DGtGnQvz9MnerGehmTxJYsWUKHDh0inr70LKeFhdCmTXKd5bSoqIhzzjmHRYv8P/t4eds90gtt2EnAYqH0l7FkiQW7qXXy8pInyIPKumJiITsbGjWyHajGJJmWLVsmRGu9pizYY0HERsYYY3xjwR4rFuzGGJ9YsMdKhw6wdi1s3ux3JcaYWsaCPVbCd6AaY0wcWbDHigW7MVHToEEDAL7//nsuuOCCqLzn0qVLycnJoXv37nz99deccsopQPzOChlLFuyx0qoVpKdDPE/VaUzAHXvssbz22mtRea+33nqLCy64gLlz53LCCSfw2WefARbs5nBSU6FjRwjA0CljEkVRURGdO3cGYNy4cZx//vkMGDCANm3acPvttx+YbvLkyfTq1YsTTzyRIUOG/OSQ/XfffZfHH3+cJ598kr59+wIHvxXceeedfPrpp+Tk5DBmzJg4rVl02QFKsdS5M3z4od9VGBM9t9wC8+ZF9z1zcuDxx6s167x585g7dy7p6em0a9eOm266iXr16vHggw8ydepU6tevzyOPPMJjjz3Gvffee2C+QYMGce2119KgQQNGjBhxyHs+/PDDjB49mrfffrsGK+UvC/ZY6twZxo+HTZugcWO/qzEmcPr168cRRxwBQMeOHfn222/ZvHkzX3311YFT7+7du5devXr5WWbcWbDHkveVkUWL4LTT/K3FmGioZss6VtLT0w/8XHqaXFXlzDPP5OWXX/axMn9ZH3ssdeni7q2f3Zi4ycvLY+bMmQdOsbtjxw6WL18e8fxlT/ebjCzYY6l5c3fOGAt2Y+ImKyuLcePGMWzYMLp27UqvXr1YWoUrmpU93W8yivi0vSLSAngBaAoo8LSqjhWRo4BXgJZAEfAbVd10uPcK/Gl7w/XuDXXqwMcf+12JMdVS1dP2muioyWl7q9Ji3w/8UVU7AnnADSLSEbgTmKaqbYBp3mNTqnNn12K3qykZY+Ik4mBX1TWqOsf7eRuwBDgOOBd43pvseeC8KNeY3Lp0gR9/tKspGWPiplp97CLSEugOfA40VdXS1FqL66oxpcJHxhiTpPy40lptVtPtXeVgF5EGwOvALaq6tUwxiut/L2++a0SkQEQK1q9fX61ik1KnTu7egt0kqYyMDDZu3GjhHieqysaNG8mo6MrgEajSOHYRScOF+kuq+ob39A8ikq2qa0QkG1hXQbFPA0+D23la7YqTTVYWNG1qwW6SVvPmzSkuLqZWNch8lpGRQfPmzas9f8TBLiICPAssUdXHwl6aBFwGPOzd/6fa1QRVly6wcKHfVRhTLWlpabRq1crvMkwVVKUrpjdwKfALEZnn3QbhAv1MESkE+nuPTbjOnWHxYgiF/K7EGFMLRNxiV9UZgFTwcr/olBNQnTvDrl2wciWccILf1RhjAs6OPI0HGxljjIkjC/Z4KB0ZY/3sxpg4sGCPhwYNXBfM/Pl+V2KMqQUs2OMlJyf6FygwxphyWLDHS04OrFgBSX46UGNM4rNgj5ecHHe/YIGvZRhjgs+CPV5Kg926Y4wxMWbBHi/HHQdNmliwG2NizoI9XkRsB6oxJi4s2OMpJ8eNZd+/3+9KjDEBZsEeTzk5sGcPLFvmdyXGmACzYI8n24FqjIkDC/Z4atcO0tMt2I0xMWXBHk9pae6EYBbsxpgYsmCPt9KRMXaZMWNMjFiwx1tODmzYAN9/73clxpiAsmCPN9uBaoyJMQv2eOva1d3PnetvHcaYwLJgj7dGjaB1a5gzx+9KjDEBZcHuh5NOgoICv6swxgSUBbsfcnNh1Sr44Qe/KzHGBJAFux9yc9397Nn+1mGMCSQLdj907+7O9mjdMcaYGLBg90PDhtC+vQW7MSYmLNj9kptrwW6MiQkLdr+cdBKsWWNHoBpjos6C3S+lO1Ct1W6MiTILdr906wapqfDll35XYowJmIiDXUT+JSLrRGRR2HP3ichqEZnn3QbFpswAysyETp2sxW6MibqqtNjHAQPKeX6MquZ4t3ejU1YtUboD1U7ha4yJooiDXVU/AX6MYS21T26uO4Xvd9/5XYkxJkCi0cd+o4gs8LpqGlc0kYhcIyIFIlKwfv36KCw2AGwHqjEmBmoa7E8CJwA5wBrg0YomVNWnVTVXVXOzsrJquNiA6NrVXS7PdqAaY6KoRsGuqj+oaomqhoB/AidHp6xaIj3dXXgjP9/vSowxAVKjYBeR7LCHg4FFFU1rKtCrl2ux79/vdyXGmICoynDHl4FZQDsRKRaRq4C/ishCEVkA9AVujVGdwZWXBzt3wsKFfldijAmIOpFOqKrDynn62SjWUjv16uXuZ81yZ300xpgasiNP/fazn0GzZtbPboyJGgt2v4m47phZs/yuxBgTEBbsiaBXL1ixAmx8vzEmCizYE0Fenrv//HN/6zDGBIIFeyLIzXVnerTuGGNMFFiwJ4LMTHcaXwt2Y0wUWLAnil694IsvoKTE70qMMUnOgj1R9OoFO3bAIjt41xhTMxbsiaL0QKWZM/2twxiT9CzYE0WrVnDssfDpp35XYoxJchbsiUIETjvNBbtdUckYUwMW7ImkTx9YvRpWrvS7EmNMErNgTySnneburTvGGFMDFuyJpFMnaNzYgt0YUyMW7IkkJQVOPRU++cTvSowxScyCPdH06QOFhbB2rd+VGGOSlAV7ointZ58xw986jDFJy4I90Zx4ojt3jHXHGGOqyYI90aSluaNQbQeqMaaaLNgT0Wmnwfz5sHmz35UYY5KQBXsi6tvXHX368cd+V2KMSUIW7ImoZ0+oVw+mTfO7EmNMErJgT0Tp6a475sMP/a7EGJOELNgTVb9+sHixjWc3xlSZBXui6tfP3Vur3RhTRRbsiSonx503xvrZjTFVZMGeqFJT4YwzXLDb+dmNMVUQcbCLyL9EZJ2ILAp77igRmSIihd5949iUWUv16wfffmvnZzfGVElVWuzjgAFlnrsTmKaqbYBp3mMTLaX97NYdY4ypgoiDXVU/AX4s8/S5wPPez88D50WnLANAu3buOqgW7MaYKqhpH3tTVV3j/bwWaFrRhCJyjYgUiEjB+vXra7jYWkLEtdqnTYOSEr+rMcYkiajtPFVVBSrcy6eqT6tqrqrmZmVlRWuxwTdgAGzYALNn+12JMSZJ1DTYfxCRbADvfl3NSzKHOOss13J/7z2/KzHGJImaBvsk4DLv58uA/9Tw/UxZRx8NJ59swW6MiVhVhju+DMwC2olIsYhcBTwMnCkihUB/77GJtoED4YsvXJeMMcZUoiqjYoaparaqpqlqc1V9VlU3qmo/VW2jqv1VteyoGRMNAwe6g5QmT/a7EmNMErAjT5NBbq7rkrHuGGNMBCzYk0FKituJ+sEHEAr5XY0xJsFZsCeLgQNh/Xob9miMqZQFe7L45S9t2KMxJiIW7MkiK8sNe3z7bb8rMcYkOAv2ZHLeefDll1Bc7HclxpgEZsGeTM49191PmuRvHcaYhGbBnkzat4e2beGtt/yuxBiTwCzYk4mI646ZPh02b/a7GmNMgrJgTzbnnQf799voGGNMhSzYk03PntC0qXXHGGMqZMGebFJS4Fe/gnffhT17/K7GGJOALNiT0Xnnwfbt8OGHfldijElAFuzJqF8/aNQIXn3V70qMMQnIgj0ZpafD4MHwxhvWHWOM+QkL9mQ1dChs2cLoMz+gRw8YNgzy8/0uyhiTCCzYk9T9M/qxgSZkf/oKc+bAxImuh2bkSL8rMyb48vNdYypRG1UW7EkoPx/+OiaN1/k15/If6rGTUAh27oTRoxPvj8yYIBk50jWiXnmFhG1UWbAnobFjYdcu+DdDacAOzuadA6/t3u1eN8ZEX36+azzt3OmuVgkkZKPKgj0JLV/u/qg+oQ9raMaFvHLgtVAICgt9LM6YACttVAG8zvn8gUcPvJZIjSoL9iTUtq07TilEKq8yhHN4m4ZsBdzzbdv6XKAxAVXaqGrFN5zPm6Sx78BridSosmBPQsOHQ0aG+3kCF5HBHi7gNcA9f/PNPhZnTICVNqp+ywuEEF7i4gOvJVKjyoI9CeXlwYgRkJkJX0pPltKOyxlHZqZ7Pi/P7wqNCabhw6FeeojLeJ5p9KOYFgdeS6RGlQV7krr/fpg2DX5zoTDluMvpw6fMGLeC++/3uzJjgisvD/5nyKe0oogX5HLAtdQTrVFVx+8CTPXl5Xl/SKsvhePvofuC52HIA36XZUygXZkyjv2ZDUkZNJgeK6FNG9eST5RQBxAtHbMTR7m5uVpQUBD35QbawIGweDEUFbkmhDEm+rZvh2bN3JHfzzwT98WLyGxVza1sOkuAoLj8cli1ys74aEwsvfEG7Njh/t8SWFSCXUSKRGShiMwTEWuK++Hcc+HII+Ff//K7EmOC65ln4IQToHdvvys5rGi22Puqak4kXxNMDGRkwKWXwuuvw7p1fldjTPAsWgSffgrXXOOuP5zArCsmSK67DvbutVa7MbHw1FNQty5ccYXflVQqWsGuwGQRmS0i15Q3gYhcIyIFIlKwfv36KC3WHKJDBzjjDPcHWFLidzXGBMf27fDCCzBkCGRl+V1NpaIV7Keq6onAQOAGEelTdgJVfVpVc1U1NysJNkzSuv56+PZbeP99vysxJjgmTIBt29y34iQQlWBX1dXe/TrgTeDkaLyvqYbzznPDsf7+d78rMSYYVOHJJ6FLFzjlFL+riUiNg11E6otIw9KfgbOARTV9X1NNaWnwu9/Be+/B11/7XY0xyW/WLJg3z7XWE3ynaalotNibAjNEZD7wBfCOqlo/gJ+uvRbq1IHHH/e7EmOS36OPQuPG8Nvf+l1JxGoc7Kr6jap2826dVPXP0SjM1MCxx8LFF7vRMRs3+l2NMclrxQp4803XWq9f3+9qImbDHYPqj390l3V56im/KzEmeY0Z47o3b7zR70qqxII9qDp3dueP+dvf3KVdjDFVs3EjPPccXHIJZGf7XU2VWLAH2YgR8MMP8OKLfldiTPJ58kl3Hbw//MHvSqrMgj3I+vaFHj3gkUdg/36/qzEmeWzb5gYfnH02dOrkdzVVZsEeZCLwpz+5HUAvveR3NcYkjyeecF0x997rdyXVYudjDzpVyM2FLVtg6VI3DNIYU7Ht26FlSzj5ZHj3Xb+rOYSdj904InDffe5gpfHj/a7GmMRX2lofOdLvSqrNWuy1gSqcdBL8+CMsW+aGbxljfiqBW+tgLXYTTsRd/XrlSvjnP/2uxpjENWqUa63fd5/fldSIBXttMWgQnH66+3q5ZYvf1RiTeFavdsF+4YWuxZ7ELNhrCxF3zosNG+Avf/G7GmMSz733uusYBOD/w4K9NunRw10+7/HHoajI72qMSRzz57ujTG+6CVq18ruaGrNgr23+/GfXer/jDr8rMSYxqMKtt7qLwd9zj9/VRIUFe23TogXceSdMnAgffOB3Ncb4b8IEmD4dHnzQnZ43AGy4Y220Zw907Qr79rkrr2dm+l2RMf7YtAnat3dDHD/7DFJT/a7osGy4o6lYejr84x9u+OODD/pdjTH+ueceN6DgqacSPtSrwoK9tjrjDLjiCje8a/58v6sxJv4++sidwfGmm6B7d7+riSrriqnNNm50521v0gQKCiAjw++KjImPrVtdd2RamrueaZJcHcm6YkzlmjSBceNg8WK46y6/qzEmfm69FVatghdeSJpQrwoL9trul790X0UffxymTPG7GmNi7/XX3fWA77gDevXyu5qYsK4Y464Sk5vrdiLNng3Nm/tdkTGxUVjoDtTr2BE++QTq1vW7oiqxrhgTuXr14LXX3MWvf/1ru0ZqHOTnw7BhLmOGDXOPTYzt3AkXXOD61SdOTLpQrwoLduN06OD6G7/4Am64wR2NZ2Ji5Ejo1w9eeQXmzHEZ069fUp/+O/GFQnD11bBwobua2PHH+11RTFmwm4MGD4b/+i/X//joo35XE0j5+TB6tGs8ln52hkLu8ejR1nKPmZEj4eWX3Sk1Bgzwu5qYs2A3h7r/fvjNb+C22+yKSzEwdqzbpVGe3bvd6ybKnnvOHYh31VXudBq1gF0A0xwqJcV1yWzYAFdeCUcd5a7UbiIXCsH69fDDD+62di2sWwfbtjHwwx2cqjupzw5SKSFECiFSKCGVfaE0Mqc3goeOgCOOcOctyc52O7OPO85O/VAdEye6Lpj+/d3BSCJ+VxQXNirGlG/rVujb151L5tVX4Ve/8ruimMrPd63l5cuhbVsYPhzy8iqZaft212e7YAGsWOFuhYXu+rIV7IDek1qPrSX12Ukm+6lDCiFSKSGFEOns4ciUraSF9pa/vKOOcuc0ad/e7RNp397d2rWzyx2W5803YcgQOOUUeO+9QIxXj3RUTFSCXUQGAGOBVOAZVX34cNNbsCeJTZtcf+ScOTBhAvkthlQ9/JLAyJGuf3vXLtfvnZLiDsIdMcL1TAGuE/yLL9wnwNy57mjFwsKDHeXp6XDCCdC6tbu1bOla202bHrw1bEj+Fyn06+ferqzMTJg2DfJydrurXG3aBN9/D8XFB2/ffANLl8K33x6cMT0dunVzQ2xyc92tY0eoc/ALebU+uJJI2fV7oN2LtP7zFW5bTJ4MDRv6XWJURBrsqGqNbrgw/xr4OVAXmA90PNw8PXr0UJMktmxR7d1bSxC9K22UCiEF1ZQU1cxM1Xvv9bvAmpk1y62HS+iDt6as0WF1X9PvL7xF9aSTVOvUOfhiq1aqgwer3n+/6n/+o7pypWpJScTLvPdet8yUFK3+tty+XXXOHNWXXlIdMUL1jDNUGzU6WGPDhqoDBqg+9JD+84qZekS9PSpSg+UlsNLt6dYvpCNklCro1y37qm7e7Hd5UQUUaCS5HMlEh30D6AV8EPb4LuCuw81jwZ5c8qfv1NdSh6iCPs3Vms6uA/mRmenCMVkNHeoCIZPtOoB39VFu1QV0PhCQe1IzVPv0Ub37btV33lH98ceoLHfWLLfsHj3cfVS2YUmJ6rJlqi++qHrddaqdOh1Yjx3U02n01bv4s55IgQolSf+7Uz30gzmDnfoMV6qC/pvf6JH1dif9+pUVz2C/ANf9Uvr4UuB/DzePBXtyGTpUNYUSfYB7VEHn0k3bsvRA62/oUL8rrIaSEtWCAv3bsQ/pNPrqbuqqgu4iXafQT2/nYe3JLO3ZfY/fldbI785bp4N5Qx/jFp1DzoGgX8fR+jJD9e8nP6daXOx3mdVW+sH8c1bobLqrgj7APSqUJO/f5mFEGuxxGxUjItcA1wAcH/CDA4Jm+XIIkcKfeJBZ9GIclzOHE7mbh3gidAOFhUkyuKq42F01asoUmDoVNm7kRmA+XfkbNzGZs5jBqezCjT5JSYEL2/tbck3N/i6LOQzmTQYDkMU6zmQKZzGZs5hM9hf/hua4s3yedZa79enjjkZOAoXLQlynT/IId7CPNM7h/3iHcwDQkNsNUitFkv6Hu2FdMYE3dOjB/mBQzWa1vsNAVdDZdNc/nZmg33d37lR97z3VW25R7djx4Ao0a6b629+qjh+vBW+vKbePPQjdTKo//d2F31IkpLcPmK86apRq//6q6enuhfR01TPPdM/Pn68aCvm9GuX78ktdcvSpqqDv8UttzneHrl+K6rBhfhcZXcSxK6YO8A3QioM7Tzsdbh4L9uRS/g7GkP6aV3W1HOueOPts1YKCqC5z6FDVE0+sQh90KKS6cKHq6NEumMoG1ejRqgsW/CSoorIzM0FVtHO43A+uHTtU339f9dZbD+mf12bNVC+9VHX8eNW1ayNebpV/f5FavFj14otVQfc2ztJr6z6reDv1g/bBXFbcgt0ti0HActzomHsqm96CPflUFH4P3rFV9aGHVBs3di/06eNGauzaVeNlVTqKY/9+NzLk8cfdKJWjjz74X92xowuo9993gVWJmOzMTBDV/uAqLlZ97jm3QZo0Obhtc3JUb7tNddIk1Q0bKlxeVEfh7NvnfpfnnOPetF491TvvVN2yJdAfzGVFGux2gJKJWOlY4cJCaNOmzFjorVvddSP/8Q831rpRIxg40B3YdMYZcOyxES+j/HHeSqeMb3jjntm03T7HnV74yy/deG+AVq3g9NPdrV8/aNEiSmsdDIf93UUiFHLj9ydPdreZM93F0MEdINW7N/TqxcLUHPre0JGNu356lOyBcfqRLnf7dndq3fffd2dMW7cOjj7aXT/g+uvdz9FavyQR1wOUqsqCPcBCIfjwQ/ePOGmS+2cEF7Q9ehw8gKdZM3fY/BFHuANp9u+H/fv58+1bWPTRBo5mPcexmjYU0poVtGYF9fHSPi0NunRxB5/06eNuFuTxtXOnu5ziZ58dvG3cCEAIYQWtWURnVtCaVbRgFS1YK8fSe2AjHn26oTsKtKQE9u6FPXvcqRfWrIGiooNH886f7z480tPdaS0uuQQGDXKPaykLduO/UMj988+a5W4LFrjW/J49Ec2+lzS+4eesoDWFtOErOrKzQw9emtupVv9zJyRVWLGCEb9cSP2VC+nCQjqziJYUkUFkv+8DsrIOHknbv7/7NpAko3RiLdJgT5JxaiYppaTAySe72/Dh7rlQCFavdifJ2rLF3UpKXKs9NZX/HtOQidOz+EGz+JGjCJF6yNtdmANYpiceEWjThtU92zDx2/MJhUpfULJYz/F8x3Gyhv5527npiu2um6VOHXexi7p14ZhjDp7wLCur1pysK1Ys2E18paS4bpMKuk7OOhoeqeBcKhkZcPPNMa7P1Mjw4a4H7uDvT1jPMaznGJbUg7seAwLY951o7HzsJqHk5bmTb2Vmus8AcPeZme75IO4QCxL7/SUG62M3Cam2jHIIKvv9xYbtPDXGmICJNNitK8YYYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwKmRsEuIveJyGoRmefdBkWrMGOMMdVTJwrvMUZVR0fhfYwxxkSBdcUYY0zARKPFfqOI/BYoAP6oqpvKm0hErgGu8R7uEZFFUVh2ojoa2OB3ETEU5PUL8rqBrV+yaxfJRKKqh59AZCrQrJyX7gHycRtRgQeAbFW9stKFihSoam4kBSYjW7/kFeR1A1u/ZBfp+lXaYlfV/hEu8J/A25FMa4wxJnZqOiomO+zhYCDI3SvGGJMUatrH/lcRycF1xRQBv49wvqdruNxEZ+uXvIK8bmDrl+wiWr9K+9iNMcYkFxvuaIwxAWPBbowxAeNbsIvIEBFZLCIhEQnE8CQRGSAiy0RkhYjc6Xc90SQi/xKRdUE9/kBEWojIdBH5yvu7HO53TdEkIhki8oWIzPfW736/a4o2EUkVkbkiErjReSJSJCILvVO3FFQ2vZ8t9kXA+cAnPtYQNSKSCjwBDAQ6AsNEpKO/VUXVOGCA30XE0H7cAXYdgTzghoD9/vYAv1DVbkAOMEBE8vwtKeqGA0v8LiKG+qpqTiTj2H0LdlVdoqrL/Fp+DJwMrFDVb1R1L/Bv4Fyfa4oaVf0E+NHvOmJFVdeo6hzv5224gDjO36qiR53t3sM07xaYkRMi0hw4G3jG71oSgfWxR89xwKqwx8UEKBhqExFpCXQHPve5lKjyuirmAeuAKaoapPV7HLgdCPlcR6woMFlEZnunZzmsaJwrpkKHOx2Bqv4nlss2pjpEpAHwOnCLqm71u55oUtUSIEdEjgTeFJHOqpr0+0xE5BxgnarOFpEzfC4nVk5V1dUicgwwRUSWet+iyxXTYI/0dAQBsRpoEfa4ufecSRIikoYL9ZdU9Q2/64kVVd0sItNx+0ySPtiB3sCvvOtBZACNRORFVb3E57qiRlVXe/frRORNXNdvhcFuXTHR8yXQRkRaiUhdYCgwyeeaTIRERIBngSWq+pjf9USbiGR5LXVEpB5wJrDU16KiRFXvUtXmqtoS93/3YZBCXUTqi0jD0p+Bs6jkA9nP4Y6DRaQY6AW8IyIf+FVLNKjqfuBG4APcjreJqrrY36qiR0ReBmYB7USkWESu8rumKOsNXAr8IqBXBMsGpovIAlwjZIqqBm5YYEA1BWaIyHzgC+AdVX3/cDPYKQWMMSZgrCvGGGMCxoLdGGMCxoLdGGMCxoLdGGMCxoLdGGMCxoLdGGMCxoLdGGMC5v8DBGVUBKd477wAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig=plt.figure()\n", "ax=fig.add_subplot(111,xlim=[-1,5],ylim=[-5,25])\n", "x,y=poly_curve(params,x_data) ###### Plotting the fitted data\n", "\n", "ax.plot(x_data,y_data,'.b',markersize=15) ###### Plotting the original data\n", "ax.plot(x,y,'r') \n", "ax.legend(['Data points','line fit'])\n", "plt.title('Polynomial of Degree %d' %(len(params)-1),fontsize=16)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Additionally, it is helpful to see what the error of the our polynomial fit is. Run the error function in the cell below to analyze how close the polynomial is to the true values.\n" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.05312027972027798" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def error(params, x_data, y_data):\n", " return np.sum((y_data - numpy.polyval(params[::-1], x_data))**2)\n", " \n", "error(params, x_data, y_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exploring other polynomials\n", "What if we didn't know what the degree of the polynomial was? What if instead we decided that our polynomial had degree 1, 2 or even 20? Change the degree variable to see how different degree polynomials fit the data. Try degrees of 1, 4, 8, 15, 20 to get started. You can also look at how the error changes as we change the degree" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 1.0, 'Polynomial of Degree 4')" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAouklEQVR4nO3deZgU1dn+8e8zM8g4gOCCSAAdVEAQBAIJi0ajgEE0aowoGP2BS9BEDb6Ju4loogaVgCSaxddE3EXjGnfgdQNBRURFdgRkECK7CIws8/z+OD3QDjPMMNNDdXfdn+vqq5eq6nqqB+4+fepUlbk7IiKS/XKiLkBERPYMBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMSEAj9DmdlgM/Ok23oz+8jMLjOzvN18rzfM7I1aKnWPMLPCxOcwuBrLpnT7zezHZvaJmRUnampUwXw3lfM3nGdmj5rZj1JVTyYws0PNbGPiczg86nqy1W4Fg6Sl/kARsE/i8V+AA4EboywqAsuAHsCCKItIfNk+ArwDXApsBtZXstgxwDagAGgJnAm8YmYPA4PcvaT2Kk4bfwXWAXtHXUg2U+BnvunuPj/x+LVE62goMQt8d/8GmBJ1HUAzoAHwhLu/VcVl3nX3rUnP/2lm/wOMBKYDf0ptiRUzszrAVt+DR2Sa2TlAZ+CPwKg9td44UpdO9nkf2MfMDgQws75mNtnMNpnZOjN71szaVLSwmR1kZpvNbGg5025K/OzeN/H8DTObaGa9zWxaYtoMM/tJOctWWkfS+/U1s+mJeT80s25mlmdmt5nZMjNbbWZjzKxe0rI7demY2ffM7N9mVpR4rzmJ96hWK9LMmprZg2a20sy+MbOPzezc5M8HWJR4+s9EPW9UZ13uPgr4kPDlnVxDYzP7u5ktTdQw28yGlFNr78RnV2xm883sosRntihpntLP7JdmdoeZfQF8AzRKTD/DzKYk/q5rzexJMzu4nHUNSXQnFic+m3+a2X5V2c7Ev6WRwJXA2qp+PlI9Cvzs05LQPfC1mfUFXgS+Bs4GfgG0ByaaWbPyFnb35cCzwLdCxMxygQsJLdc1SZMOA0YT/tOeQehaeTK5H3Y36zgcuBMYTuiiqgs8D/wNaAoMBn4P/AwYVslncTChhXwJ0DdR5wXA/ZUst5PEl8ubwEnA9cDpwCfAQ0mBe1+iZoBbCF1Mv9zddSV5GWhRGrJmtg8wEegH3AScDPwH+JuZXZ5Uazt2fN4DEvUOBU6oYD03AK0Jf/OfAMVmdgnwFDCT0MV0MeFv9qaZNUha13DgHmA8cCpwFeGzfjnxb6YydwCz3f2hKswrNeXuumXgjRB8DrQhdM3tS/hPuQ14NjHPVGAekJe0XEtgCzAy6bU3gDeSnv8w8d4/SHrt1MRr3csstwVolfTagYkark96bXfq2AIcWs56x5fZ/qeBhUnPCxPzDa7g87LE53QuUALsX9H2V7D8ZYn3/2GZ18cDXwK5ieeH76qOMsvelJg3r4LpFyemd0s8/x1QnPx5J17/X2Bl6fsAjwIrgIKkeZomll1Uzmc2DbCk1+sT+tP/VWY9LQn7JK5IWn4bcGOZ+Y5OvO/plWz/Dwi/KNqV+Td9eNT/v7L1phZ+5ptNCMnVhB1fjwAXJFqk3wXGelL/sLsvBCYBx1X0hu7+BqFld3HSyxcDH7t72X7yee4+L2nZLwkBWNoq3d065rr7Z2W2D+DVcra7uZlZRdthZvuY2e1mtoAQLFuAhwjh36qi5SpwLLA08dkkexhoDLTbzferitJtK+1P7wu8CyxMdHHlJXYSvwrsn1RDd+Ald99Y+kbuvoywI7k8z3oicRN6EAYBPFJmPUsIn/uxifn6EHoJys73LmFH9bFUwMz2Av4BjHL3mZV+EpIS2mmb+X5CGKWzHljs7sUAZtacEBjLyllmOXBIJe/7N2BEoi+/PiFsLitnvtXlvPYNkJ94vO9u1rGmzPPNu3g9D8gFtlK++4HehB3Y04ENwPcJXRD5FSxTkf2oeBtKp6dai8R96XoPJPyC2FLB/Psn7psSvnTL+i9waDmvl92uAxP34ytYz5oy882vYL79K3gd4ArCv40/245hqwWJ+wZm1sDdKxvdJLtJgZ/5ZviOUTrJ1hBahgeVM+0gyg/qZA8SRk0MJvzH3Ej49bC7alpHtZhZPnAacJO7j056vUM133I1ofusrIOSpqdaP+Bzd1+SeL6KEOQ77VBPmJO4X8aOME7WpILlyo7IWZW4Hwx8Ws7868vMdyI7fyEnTy9PO8Jnt7ScadOAj4BOu1heqkGBn6XcfYOZfQD0N7Ob3H0bgJkdAvQkjNff1fJfmdkjhK6c+sBj7v7Vnq6jBuoSWv9lW8ODq/l+bxK24Wh3n5T0+jmEEE5pt4SFYZmdCC3hUq8AlxO+BMprwZeaAvQzs4LSbh0za0roWy/vV0pZ7xBC/XB3f2AX840j7A852N3HVeF9kw0HxpR5rS9wDWE/y5yyC0jNKfCz2+8IozVeMLO/EoL7ZsIOuaqM7f4rO/rx/x5hHbvN3deZ2RTgN2a2jLBT8wLCOPnqGENoWT9tZjcQutF+RujHvrj0i6yaupnZNkI306GEUTEnAQ8Af06abxRhlNPbZjaKEIr1gCMIO9hPS8x3S+I9XjWzEYQvv98RunQqPYgr8WV/FXCPmTUmjBZaR/jsjiPs4H7U3ReY2e3A3Ykhtm8Sdgy3IHwu97n76xWsYzY79s8AYZho4uG7FfxqlRpS4Gcxd3/FzE4mDF98gtDv/QZwtbt/UYXlPzazucBX7j4tqjpqYCBhX8Q9wKbEuocCL+zuGyV+qRxHGEY4nHBw1RzgPHd/uIZ1TkzcbyS0wN8D+rr7t3ZUJ77EehL2SVxDCOC1iTqeSppvZuLzvpOwzUuB2wkt6MKqFOTu/zCzJYRhlucQsmIp8DZhf0jpfNeb2SzCUcWXErqHlgATCCOzJI3Yt3fOi+yQaLXNAn7u7v+Muh6pPjOrT9i5+qK7Xxh1PRINBb7sJDHC53BCt8vhhL7cTdFWJbvDzP5C6Iv/AvgO4ZdNZ+B77v5xlLVJdNSlI+W5iNBtMBc4R2GfkfIJ3ThNCF1o7wG9Ffbxpha+iEhM6EhbEZGYSKsunQMOOMALCwujLkNEJKN88MEHK929cWXzpVXgFxYWMnXq1KjLSLkpU2D0aJg7F1q3hqFDoXv3qKsSkWxhZourMp+6dGrZsGHQqxeMHQvTpsETT4Tnwyo7sa+ISIop8GvRlCkwYgRs3Ah5vpnvsJSSkvB8xIgwXURkT1Hg16LRo2FTYkDjRI5hTNJpXIqLw3QRkT0lrfrws83cuVA66nUSR/ML/kY+myhmb0pKYJ4OPJcMtmXLFoqKiiguLo66lNjIz8+nefPm1KlTp1rLK/BrUevWMH06lJTAOPrwP9xFT97h/+hFTk6YLpKpioqKaNCgAYWFheziOjSSIu7OqlWrKCoqomXLltV6D3Xp1KKhQyE/cZmNtziWLeTRh3AW2fx8+NWvIixOpIaKi4vZf//9FfZ7iJmx//771+gXlQK/FnXvDldeCQUFsCmnPpPpQW/GU1AQXtfQTMl0Cvs9q6aftwK/lt18M0yYAGedBTOb9ua7TOPNp1dx881RVyYicaPA3wO6d4fHHoNLnupDDk7X9eVeE0Ikq02ZAgMHQpcu4T4Vw5Jzc3Pp1KkTRx55JB07duRPf/oTJSW7vsbLokWLePTRR2u+8kpcdNFFzJy56wuhPfvss5XOk0oK/D3pe9+DBg1gfEXXhhbJTrV1AOLee+/N9OnT+fTTTxk3bhwvv/wyN1fy83lPBf59991Hu3btdjnPng583D1tbl26dPGsd+qp7oceGnUVIjU2c+bMKs03ebJ7QYF7GKT87VtBQZheXfXq1fvW8wULFvh+++3nJSUlvnDhQj/mmGO8c+fO3rlzZ580aZK7u3fr1s332Wcf79ixo48cObLC+ZItXLjQ27Rp4+ecc44fccQR/tOf/tQ3bNjg7u7jx4/3Tp06efv27f3888/34uJid3c/7rjj/P33399e5/XXX+9HHXWUd+vWzZcvX+6TJk3yfffd1wsLC71jx44+f/58Hz16tLdt29Y7dOjgZ599drnbXN7nDkz1KmRs5CGffItF4P/5z+FjX7Ag6kpEaqSqgT9ggLtZ+YGfkxOmV1fZwHd3b9iwoS9fvtw3bNjgmzZtcnf3uXPnemm+vP76637yySdvn7+i+ZItXLjQAZ84caK7u59//vl+5513+qZNm7x58+Y+Z84cd3c/77zzfNSoUe7+7cAH/Pnnn3d396uuusr/8Ic/uLv7oEGD/Mknn9y+nqZNm27/wlizZk2521yTwFeXzp7Wp0+4nzAh2jpE9pDkAxDLqs0DELds2cLPf/5zOnToQP/+/SvsOqnqfC1atODoo48G4Nxzz2XixInMmTOHli1b0jpxUM2gQYN46623dlp2r7324pRTTgGgS5cuLFq0qNx1HHXUUfzsZz/j4YcfJi8v9YdJKfD3tDZtoFkzGDcu6kpE9ojWrSGngqRJ9QGIn332Gbm5uRx44IGMGjWKJk2a8NFHHzF16lQ2b95c7jJVna/skMjdGSJZp06d7fPn5uaydevWcud78cUXufTSS5k2bRrf+973KpyvuhT4e5oZnHhiCPwU/zFF0lHyAYhlpfIAxBUrVnDJJZdw2WWXYWasW7eOpk2bkpOTw0MPPcS2bdsAaNCgAevXr9++XEXzlfX5558zefJkAB599FGOOeYY2rRpw6JFi5g/fz4ADz30EMcdd1yVa06upaSkhCVLlnD88cdz++23s27dOr7++utqfRYVUeBH4aSTYO1aePfdqCsRqXXJByCWtvRzckjJAYibNm3aPiyzd+/enHjiiQxLDP355S9/yQMPPEDHjh2ZPXs29erVA0K3SW5uLh07dmTUqFEVzldWmzZtuOeee2jbti1r1qzhF7/4Bfn5+dx///3079+fDh06kJOTwyWXXFLl+gcMGMCdd95J586dmTdvHueeey4dOnSgc+fO/OpXv6JRo0bV/3DKkVbXtO3atatn4wVQdrJ2LRxwAFxzDdx6a9TViFTLrFmzaNu2bZXnL70Q0Lx50KpVZl0IaNGiRZxyyinMmDEj6lLK/dzN7AN371rZsjXeK2BmLYAHgSaAA/e6+2gz2w8YCxQCi4Cz3H1NTdeXFRo1gqOPhpdeUuBLbHTvnjkBn61S0aWzFfiNu7cDugOXmlk74Fpggru3AiYknkupfv3CqTSXLo26EhGpRGFhYVq07muqxoHv7svcfVri8XpgFtAMOA14IDHbA8DpNV1XVjnppHD/yivR1iEisZHSnbZmVgh0Bt4Fmrj7ssSk5YQun/KWGWJmU81s6ooVK1JZTnrr0CEMz3z55agrEZGYSFngm1l94CngCnf/Knla6ZFm5S3n7ve6e1d379q4ceNUlZP+zEK3zmuvwZYtUVcjIjGQksA3szqEsH/E3Z9OvPxfM2uamN4U+DIV68oq/frB+vUwaVLUlYhIDNQ48C0cPvZPYJa7j0ya9DwwKPF4EPBcTdeVdXr1gjp1wmgdEdlt9evXB+CLL77gzDPPTMl7zp49m06dOtG5c2cWLFhAz549gT13ls3alIoW/tHAecAJZjY9cesHDAf6mNk8oHfiuSRr0ACOPRZeeCHqSkQy2ne+8x3+/e9/p+S9nn32Wc4880w+/PBDDjvsMN555x1AgQ+Au090d3P3o9y9U+L2kruvcvde7t7K3Xu7++pUFJx1Tj0VZs0KZ5gSkWpZtGgR7du3B2DMmDGcccYZ9O3bl1atWnH11Vdvn++1116jR48efPe736V///47nbrgpZde4q677uJvf/sbxx9/PLDjV8S1117L22+/TadOnRg1atQe2rLUSv3p2GT3nHZaOOTwuefgqquirkakeq64IhxXkkqdOsFdd1Vr0enTp/Phhx9St25d2rRpw+WXX87ee+/NLbfcwvjx46lXrx633347I0eO5MYbb9y+XL9+/bjkkkuoX78+V1555bfec/jw4YwYMYIXMvgXuQI/aoccAp07w7PPKvBFUqRXr140bNgQgHbt2rF48WLWrl3LzJkzt5/iePPmzfTo0SPKMvc4BX46OP10uOkmWL4cDjoo6mpEdl81W+K1pW7dutsfl56O2N3p06cPjz32WISVRUtny0wHp58erhDxn/9EXYlI1urevTuTJk3afirjDRs2MHc39p2VPa1yJlLgp4MOHeDQQ0O3jojUisaNGzNmzBgGDhzIUUcdRY8ePZg9e3aVly97WuVMpNMjp4vf/AbuvhtWrgzDNUXS3O6eHllSoyanR1YLP12cfjps3qyTqYlIrVHgp4uePcNFUZ55JupKRCRLKfDTRW5uGJP/wgtQXBx1NSJVkk5dwnFQ089bgZ9OzjornExNp0yWDJCfn8+qVasU+nuIu7Nq1SryK7oifBVoHH46OeGE0K0zdiz85CdRVyOyS82bN6eoqIhYXcciYvn5+TRv3rzayyvw00leHvz0p/DQQ7BhA9SrF3VFIhWqU6cOLVu2jLoM2Q3q0kk3Z58NGzfqlMkiknIK/HRz7LHh9Apjx0ZdiYhkGQV+usnNhTPPhBdfDDtwRURSRIGfjs4+OwzN1Ll1RCSFFPjpqGdPaNYMHn886kpEJIso8NNRTg4MGBBOs7ByZdTViEiWUOCnq//3/2DLFojxubtFJLUU+OnqqKPCJd4efDDqSkQkSyjw09mgQTB1KsycGXUlIpIFFPjpbODAMExTrXwRSQEFfjpr0gROOimcamHbtqirEZEMp8BPd4MGwRdfwIQJUVciIhlOgZ/ufvxjaNQIHngg6kpEJMMp8NNd3bpwzjnw9NOwenXU1YhIBlPgZ4IhQ8KpFrTzVkRqQIGfCTp2hG7d4B//AF1dSESqSYGfKS6+GGbPhrffjroSEclQCvxMcfbZ0LBhaOWLiFSDAj9TFBTAeefBv/+tE6qJSLUo8DPJxRfD5s0wZkzUlYhIBlLgZ5L27eEHP4B77tGRtyKy2xT4meaKK2DRInj++agrEZEMk5LAN7N/mdmXZjYj6bX9zGycmc1L3O+binXF3mmnQWEhjBoVdSUikmFS1cIfA/Qt89q1wAR3bwVMSDyXmsrNhcsvD8MzP/gg6mpEJIOkJPDd/S2g7HH/pwGlJ4B5ADg9FesS4MILoX59GD066kpEJIPUZh9+E3dflni8HGhSi+uKl4YN4YILwkXOly2rfH4REfbQTlt3d6DccwKY2RAzm2pmU1esWLEnyskOv/pVGKmjVr6IVFFtBv5/zawpQOL+y/Jmcvd73b2ru3dt3LhxLZaTZQ47DPr3h7/+FdasiboaEckAtRn4zwODEo8HAc/V4rri6brrYP16uPvuqCsRkQyQqmGZjwGTgTZmVmRmFwLDgT5mNg/onXguqdSxI5xyCtx1F3z9ddTViEiay0vFm7j7wAom9UrF+8suXH899OwJ994Lv/511NWISBrTkbaZrkcPOP54GDEiXCRFRKQCCvxs8NvfhuGZOnWyiOyCAj8bnHBCuN16q/ryRaRCCvxsceutsGKFxuWLSIUU+Nmie3c49VS4805YXfYsFyIiCvzscsst8NVXLL3iTgYOhC5dYOBAmDIl6sJEJB0o8LNJhw581P4c9nvoLt55/HOmTYMnnoBevWDYsKiLE8l+U6aQ1o0tBX4WmTIFzpp3G44xnGsAKCmBjRvDqM10+8cnkk2GDQuNq7FjSdvGlgI/i4weDfO+OZg7uYqBPM7RTNw+rbhY+3NFasuUKaFRtXEjeOI0kenY2FLgZ5G5c8M/tju4miKaMZqhGCVA+Mc3b17EBYpkqdGjYdOm8Pg8HuR3/J5ctgLp1dhS4GeR1q0hJwc2Uo+ruYMuTON87gfC661bR1ygSJYqbWztxypG8mt6M55t5ALp1dhS4GeRoUMhPz88foyBTORobucaDmAF+fnhFPoiknqlja0/ch2NWMul3AMYkF6NLQV+FuneHa68EgoKICfHuJh/sA9f8efcX3PllWG6iKTe0KHwwzqTGML/chdXMIMO26elU2NLgZ9lbr4ZJkyAs86CvbscyQvtrmHgtoe5+ejXoi5NJGt177iJsfUvYLEdwu/tJiC07AsKSKvGVkpOjyzppXv3pH9gxTdApyfhkkvgk0+gXr1IaxPJSjffzAGr5vLl6HGcPLk+8+ZBq1ah5Z8uYQ9g7uVeajYSXbt29alTp0ZdRvZ56y047ji4/HL485+jrkYku0ydCt26wfnnw333RVKCmX3g7l0rm09dOnFw7LGhE/Evf4FXX426GpHsUVwcgv6gg8KA+zSnwI+L4cPhyCNh8GBYuTLqakSyw1VXwYwZoWXfqFHU1VRKgR8Xe+8NjzwSzqT585/vOBxQRKrn+efh7rvhf/4HTjop6mqqRIEfJx07hvPmP/ts+IcqItWzdGnoyuncGf74x6irqTIFftz8+tfw4x+H+0mToq5GJPNs3gwDBsA338Djj0PdulFXVGUK/LjJyYEHH4TCQujfH5Yvj7oikcwydChMnBj67dPlENoqUuDHUaNG8PTTsG4dnHlmaKmISOXuvRf+/ne4+urQys8wCvy46tAB/vWv0K0zeHA4w5OIVOyNN+Cyy+BHP4Lbbou6mmrRkbZxdvbZsGgRXHstHHJIGLopIjv7+GM47bRw+Oxjj0FubtQVVYsCP+6uvhoWL4bbb4cWLeDSS6OuSCS9LF4MfftCgwbwyiuw775RV1RtCvy4MwunW/jii/Bzda+9wjh9EQnDL/v0CZeumjgxNIoymPrwBfLywoU4+/WDIUNC375I3C1dCscfH0ayvfQStG8fdUU1psCXoG5deOqp8NP1oovCSASRuPr88x1h/8or0LNn1BWlhAJfdsjPh2eeCS39X/wCbrxRp2CQ+Pnoo3BO4y+/zKqwBwW+lJWfH069cMEF8Ic/wIUXhiMLReJg/Hj4wQ/CKJyJE7Mq7EGBL+XJywtHEd54I9x/f/hp+8UXUVcVG1OmwMCB0KVLuJ8yJeqKYsAdRo0KXZqFhTB5clb02ZelwJfymYXrJT72WPiJ27lzOPBEatWwYdCrV9iHPm0aPPFEeD5sWNSVZbH168MxKaXnmXr7bWjePOqqaoUCX3ZtwAB4770w9rhXL7jmmnDRB0m5KVPCNTQ2btyx66SkJDwfMUIt/Vrx9tuhMfPUU3DHHeGUIw0bRl1VrdE4fKlcu3bw/vvhasx33AH/+U/o6unWDQhBNHo0zJ0bziWVbtfxrKmUbp97+MJcuzacy2jduvB8yxbe/MNm+mzcQh5byGMr28hlC3XYQh22barD+Ovz6H5bAeyzTwilhg3DNYrNUrm5WaXCv93GjfDb38Jdd4UunNdfD1eGy3K1fk1bM+sLjAZygfvcvcLj93VN2wzw6qth2GZREQwezJ373sZN/2jKpk0hy3Jywn7fK68MPUKZbtiw0LqudPvc4b//DclSVATLlu18W7UqBP2WLakrMDc3fAE0bhwus3fQQdC06Y7HLVpAy5bhvk6dnRbP5i/rcv92dZ0xJz1O//evgSVLwpHlw4dD/fpRl1sjVb2mba0GvpnlAnOBPkAR8D4w0N1nlje/Aj9DfPUV3HYbJSNHsWHLXozk14xmKGvYb/ssBQUwYUJmh8eUKaEXa+PGHa/VYTOHM5+j9prDHRfM5uCNs2HOHJg9O7TWk+Xnh/AtvTVuHFrljRp9+z4/H/bai9/+vg6vjK9Dse/FVvLIY2uifb+FuraFE4/fwu+u3LTjl8G6deFvsXYtrFgRxowvXx6+XNav/3Ytubk7wr9lSzj0UMZOb82IF9rySXErvqFuVn1Z7/y3c47ndW7ht/RkMhtad6befaPDiJwskC6B3wO4yd1/lHh+HYC7l3uJGAV+Zhl6ygJ+8OI1nMlTrKc+93Ap93ApRbQgJwfOOivs881IK1dy04DZFE2YQxtmcwSzacMcDuUz8ti2Y75mzaBNGzjiiHBr0yYEa9OmIcx3o7ulvC+YUrv9BbphQwj/xYth4cKdb0nXQdhGDp9xaGIrj+Czvdpy6V+OoH3/thl73piBA8OOb/NtnMTLXM9t9GQyX9CUG+0WNvYfxKNjM/MEaOVJl8A/E+jr7hclnp8HdHP3y5LmGQIMATj44IO7LF68uNbqkdTq0iWMJDmSGdzArZzNWBzjJfpxL0NY2flEJk9L/dWAUtYNsXUrfPZZaJ2XttJLH69atX22Yuoyl9bsiP4jmEMbCjq14a0PG6Ruw9jRDVFcHHbY1lare3D/DXz073kcwSzaMosjmE1bZtGaudQl6biLJk3CPpy2bcOt9PFBB1Vr38Ge6kI6qUMRXWaM4SLuo5DFLOZghnMt93M+35BPly6QTW3LqgY+7l5rN+BMQr996fPzgLsrmr9Lly4umWPAAPecHPfQQ+p+CAv9Fq73ZTRxB99QZx/3c85xf/JJ95UrU7LOG290LyhwNwvrzMkJz2+8sYIFSkrcV6xwf+cd9zFj3K+91v30093btnWvU2dH8eDepIn7sce6DxniPnKk//HYl/wwW+A5bP3WbKXrHTgwJZu0k8mTw2fbpUu4nzw59ev47nd9p20C9xy2+mHM86GH/cf99tvdBw9279bNfZ99vj1jo0buPXq4X3ih+4gR7i++6L5wofu2bRWuc7f/drujpMR91iz3O+5w7959e52v0dv7M9bz2LxH/nZRAaZ6FTJZXTpSbRV1QeSxhR/XHcff+zzFgZOf29Fabt8+jITo2DFcgKV9+3DK2Rqury7FHJ6/lLF/KuLIfZbAggUwb15oRs6bF/q4txeXF85pXrYbpk2bnbovUtrFkmYGDgxj/Mu77k1OThiW/uijSS+6h4PvZs2CmTPDfenjFSt2zFdQED7T0l8EhYXQogUfrjqYH/6sGV9t2nnHcbU+y//+F2bMgE8+gXfegbfeCq8BdOnC513P4OQHzmJG8eGpWV+aS5cunTzCTttewFLCTttz3P3T8uZX4GeeSrsgtm4Nyfnmm+E2eTJ8/fWON9h//3CQS/PmofugXr0wYqJevbCjcdu2cNu6lRee2MiXc9bQkLU0Yi37s4pmLKUxK79dlFnoR2/dOoR76a1167DDspzRKtXevgyV0i+zlSt3fAEkfyEsWfKt2UowvuA7LKUZq9if1ezHavZjDftxWKf6nDsoL3wh5yVGi2/cGPZFbNgAq1eHs1cuXRred/XqHW/cogUcd1xoTPTpE75kyN6/XXnSIvAThfQD7iIMy/yXu99a0bwK/MxU2i87b17I1V32y5aUhB2Jn3wSWmhLloRhjEuXhpNVbdgQvhC2bt1p0U05Bawq2Ze1NGIN+7KGfVlKM4poThHNqdeqGX99vnn4D5+fH832ZZBaD8QNG8Lf9/PP+cPPP6fk8yUczOc0Yyn7sob9WM2+rKERa8lhFzm0115hRFOzZuHWvHn4FdG+ffileOCBFS6arX+7stIm8HeHAl+2++abkEK5ueGWk8PAc2z3uiGkUnsqEHfVhZRn2zj3p5u4/3+3hi/6rVvDjPXqhZ8bu/GLLK4U+JJ1srlPPdvpb1e7qhr4OpeOZIzu3UN3Q0FBaNFDuC8oCK8rMNKX/nbpQS18yThx6ZfNRvrb1Q516YiIxIS6dERE5FsU+CIiMaHAFxGJCQW+iEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmFPgiIjGhwBcRiQkFvohITCjwRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhMKPBFRGJCgS8iEhM1Cnwz629mn5pZiZl1LTPtOjObb2ZzzOxHNStTRERqKq+Gy88AzgD+kfyimbUDBgBHAt8BxptZa3ffVsP1iYhINdWohe/us9x9TjmTTgMed/dv3H0hMB/4fk3WJSIiNVNbffjNgCVJz4sSr4mISEQq7dIxs/HAQeVMusHdn6tpAWY2BBgCcPDBB9f07UREpAKVBr67967G+y4FWiQ9b554rbz3vxe4F6Br165ejXWJiEgV1FaXzvPAADOra2YtgVbAe7W0LhERqYKaDsv8iZkVAT2AF83sVQB3/xR4ApgJvAJcqhE6IiLRqtGwTHd/Bnimgmm3ArfW5P1FRCR1dKStiEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmFPgiIjGhwBcRiQkFvohITCjwRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmFPgiIjGhwBcRiQkFvohITCjwRURiQoEvIhITNQp8M7vTzGab2cdm9oyZNUqadp2ZzTezOWb2oxpXKiIiNVLTFv44oL27HwXMBa4DMLN2wADgSKAv8Fczy63hukREpAZqFPju/pq7b008nQI0Tzw+DXjc3b9x94XAfOD7NVmXiIjUTCr78C8AXk48bgYsSZpWlHhtJ2Y2xMymmtnUFStWpLAcERFJllfZDGY2HjionEk3uPtziXluALYCj+xuAe5+L3AvQNeuXX13lxcRkaqpNPDdvfeuppvZYOAUoJe7lwb2UqBF0mzNE6+JiEhEajpKpy9wNXCqu29MmvQ8MMDM6ppZS6AV8F5N1iUiIjVTaQu/EncDdYFxZgYwxd0vcfdPzewJYCahq+dSd99Ww3WJiEgN1Cjw3f3wXUy7Fbi1Ju8vIiKpoyNtRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhM2I7L0EbPzNYDc6KuoxYdAKyMuohapO3LXNm8bZD929fG3RtUNlNNL3GYanPcvWvURdQWM5uq7ctc2bx92bxtEI/tq8p86tIREYkJBb6ISEykW+DfG3UBtUzbl9myefuyedtA2wek2U5bERGpPenWwhcRkVqiwBcRiYm0C3wz629mn5pZiZllxTAqM+trZnPMbL6ZXRt1PalmZv8ysy/NbEbUtaSambUws9fNbGbi3+XQqGtKJTPLN7P3zOyjxPbdHHVNtcHMcs3sQzN7IepaUs3MFpnZJ2Y2vbLhmWkX+MAM4AzgragLSQUzywXuAU4C2gEDzaxdtFWl3Bigb9RF1JKtwG/cvR3QHbg0y/5+3wAnuHtHoBPQ18y6R1tSrRgKzIq6iFp0vLt3quxYg7QLfHef5e7ZdLTt94H57v6Zu28GHgdOi7imlHL3t4DVUddRG9x9mbtPSzxeTwiNZtFWlToefJ14Widxy6qRHGbWHDgZuC/qWqKWdoGfhZoBS5KeF5FFgREnZlYIdAbejbiUlEp0d0wHvgTGuXtWbR9wF3A1UBJxHbXFgdfM7AMzG7KrGSM5tYKZjQcOKmfSDe7+3J6uR6QyZlYfeAq4wt2/irqeVHL3bUAnM2sEPGNm7d09K/bHmNkpwJfu/oGZ/TDicmrLMe6+1MwOBMaZ2ezEr+6dRBL47t47ivVGZCnQIul588RrkiHMrA4h7B9x96ejrqe2uPtaM3udsD8mKwIfOBo41cz6AfnAPmb2sLufG3FdKePuSxP3X5rZM4Ru5HIDX106te99oJWZtTSzvYABwPMR1yRVZGYG/BOY5e4jo64n1cyscaJlj5ntDfQBZkdaVAq5+3Xu3tzdCwn/9/4vm8LezOqZWYPSx8CJ7OLLOu0C38x+YmZFQA/gRTN7NeqaasLdtwKXAa8Sdvg94e6fRltVapnZY8BkoI2ZFZnZhVHXlEJHA+cBJySGvU1PtBazRVPgdTP7mNA4GefuWTd0MYs1ASaa2UfAe8CL7v5KRTPr1AoiIjGRdi18ERGpHQp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhM/H/2Hw1dpgo4cgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "degree = 4 # use this to change the degree of your polynomial\n", "D = data_matrix(x_data, degree)\n", "params = np.linalg.lstsq(D, y_data, rcond=None)[0]\n", "#plotting the figure\n", "fig=plt.figure()\n", "ax=fig.add_subplot(111,xlim=[-1,5],ylim=[-25,25])\n", "x,y=poly_curve(params, x_data) ###### Plotting the fitted data\n", "\n", "ax.plot(x_data,y_data,'.b',markersize=15) ###### Plotting the original data\n", "ax.plot(x,y,'r') \n", "ax.legend(['Data points','line fit'])\n", "plt.title('Polynomial of Degree %d' %(len(params)-1),fontsize=16)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.05312027972027535" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "error(params, x_data, y_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Manipulating Data Points.\n", "What happens if we remove data points or the noise of our data is larger? Run the cells below to look at what happens if we remove data points or made our observations a lot noisier. You can change how many data points to drop or how much noise you want. Additionally, see how noise and dropping data points interacts with the degree of the polynomial and the error it produces." ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [], "source": [ "\"\"\"NOTE: We only have 10 data points, so choosing num_points greater than 10 or less than will break this function!!!\"\"\"\n", "def reduce_points(x_data, y_data, num_points=0):\n", " x = np.array(x_data)\n", " y =np.array(y_data)\n", " indexes = np.random.choice(len(x_data), size= num_points, replace=False)\n", " return x[indexes], y[indexes] " ] }, { "cell_type": "code", "execution_count": 94, "metadata": {}, "outputs": [], "source": [ "\"NOTE: try ranges of 0.01 to 10 to see the differences\"\"\"\n", "def noisy_points(y_data, noise_factor = 1):\n", " return np.array(y_data) + noise_factor * np.random.randn(len(y_data))" ] }, { "cell_type": "code", "execution_count": 109, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.05312027972027508\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAouklEQVR4nO3deZgU1dn+8e8zM8g4gOCCSAAdVEAQBAIJi0ajgEE0aowoGP2BS9BEDb6Ju4loogaVgCSaxddE3EXjGnfgdQNBRURFdgRkECK7CIws8/z+OD3QDjPMMNNDdXfdn+vqq5eq6nqqB+4+fepUlbk7IiKS/XKiLkBERPYMBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMSEAj9DmdlgM/Ok23oz+8jMLjOzvN18rzfM7I1aKnWPMLPCxOcwuBrLpnT7zezHZvaJmRUnampUwXw3lfM3nGdmj5rZj1JVTyYws0PNbGPiczg86nqy1W4Fg6Sl/kARsE/i8V+AA4EboywqAsuAHsCCKItIfNk+ArwDXApsBtZXstgxwDagAGgJnAm8YmYPA4PcvaT2Kk4bfwXWAXtHXUg2U+BnvunuPj/x+LVE62goMQt8d/8GmBJ1HUAzoAHwhLu/VcVl3nX3rUnP/2lm/wOMBKYDf0ptiRUzszrAVt+DR2Sa2TlAZ+CPwKg9td44UpdO9nkf2MfMDgQws75mNtnMNpnZOjN71szaVLSwmR1kZpvNbGg5025K/OzeN/H8DTObaGa9zWxaYtoMM/tJOctWWkfS+/U1s+mJeT80s25mlmdmt5nZMjNbbWZjzKxe0rI7demY2ffM7N9mVpR4rzmJ96hWK9LMmprZg2a20sy+MbOPzezc5M8HWJR4+s9EPW9UZ13uPgr4kPDlnVxDYzP7u5ktTdQw28yGlFNr78RnV2xm883sosRntihpntLP7JdmdoeZfQF8AzRKTD/DzKYk/q5rzexJMzu4nHUNSXQnFic+m3+a2X5V2c7Ev6WRwJXA2qp+PlI9Cvzs05LQPfC1mfUFXgS+Bs4GfgG0ByaaWbPyFnb35cCzwLdCxMxygQsJLdc1SZMOA0YT/tOeQehaeTK5H3Y36zgcuBMYTuiiqgs8D/wNaAoMBn4P/AwYVslncTChhXwJ0DdR5wXA/ZUst5PEl8ubwEnA9cDpwCfAQ0mBe1+iZoBbCF1Mv9zddSV5GWhRGrJmtg8wEegH3AScDPwH+JuZXZ5Uazt2fN4DEvUOBU6oYD03AK0Jf/OfAMVmdgnwFDCT0MV0MeFv9qaZNUha13DgHmA8cCpwFeGzfjnxb6YydwCz3f2hKswrNeXuumXgjRB8DrQhdM3tS/hPuQ14NjHPVGAekJe0XEtgCzAy6bU3gDeSnv8w8d4/SHrt1MRr3csstwVolfTagYkark96bXfq2AIcWs56x5fZ/qeBhUnPCxPzDa7g87LE53QuUALsX9H2V7D8ZYn3/2GZ18cDXwK5ieeH76qOMsvelJg3r4LpFyemd0s8/x1QnPx5J17/X2Bl6fsAjwIrgIKkeZomll1Uzmc2DbCk1+sT+tP/VWY9LQn7JK5IWn4bcGOZ+Y5OvO/plWz/Dwi/KNqV+Td9eNT/v7L1phZ+5ptNCMnVhB1fjwAXJFqk3wXGelL/sLsvBCYBx1X0hu7+BqFld3HSyxcDH7t72X7yee4+L2nZLwkBWNoq3d065rr7Z2W2D+DVcra7uZlZRdthZvuY2e1mtoAQLFuAhwjh36qi5SpwLLA08dkkexhoDLTbzferitJtK+1P7wu8CyxMdHHlJXYSvwrsn1RDd+Ald99Y+kbuvoywI7k8z3oicRN6EAYBPFJmPUsIn/uxifn6EHoJys73LmFH9bFUwMz2Av4BjHL3mZV+EpIS2mmb+X5CGKWzHljs7sUAZtacEBjLyllmOXBIJe/7N2BEoi+/PiFsLitnvtXlvPYNkJ94vO9u1rGmzPPNu3g9D8gFtlK++4HehB3Y04ENwPcJXRD5FSxTkf2oeBtKp6dai8R96XoPJPyC2FLB/Psn7psSvnTL+i9waDmvl92uAxP34ytYz5oy882vYL79K3gd4ArCv40/245hqwWJ+wZm1sDdKxvdJLtJgZ/5ZviOUTrJ1hBahgeVM+0gyg/qZA8SRk0MJvzH3Ej49bC7alpHtZhZPnAacJO7j056vUM133I1ofusrIOSpqdaP+Bzd1+SeL6KEOQ77VBPmJO4X8aOME7WpILlyo7IWZW4Hwx8Ws7868vMdyI7fyEnTy9PO8Jnt7ScadOAj4BOu1heqkGBn6XcfYOZfQD0N7Ob3H0bgJkdAvQkjNff1fJfmdkjhK6c+sBj7v7Vnq6jBuoSWv9lW8ODq/l+bxK24Wh3n5T0+jmEEE5pt4SFYZmdCC3hUq8AlxO+BMprwZeaAvQzs4LSbh0za0roWy/vV0pZ7xBC/XB3f2AX840j7A852N3HVeF9kw0HxpR5rS9wDWE/y5yyC0jNKfCz2+8IozVeMLO/EoL7ZsIOuaqM7f4rO/rx/x5hHbvN3deZ2RTgN2a2jLBT8wLCOPnqGENoWT9tZjcQutF+RujHvrj0i6yaupnZNkI306GEUTEnAQ8Af06abxRhlNPbZjaKEIr1gCMIO9hPS8x3S+I9XjWzEYQvv98RunQqPYgr8WV/FXCPmTUmjBZaR/jsjiPs4H7U3ReY2e3A3Ykhtm8Sdgy3IHwu97n76xWsYzY79s8AYZho4uG7FfxqlRpS4Gcxd3/FzE4mDF98gtDv/QZwtbt/UYXlPzazucBX7j4tqjpqYCBhX8Q9wKbEuocCL+zuGyV+qRxHGEY4nHBw1RzgPHd/uIZ1TkzcbyS0wN8D+rr7t3ZUJ77EehL2SVxDCOC1iTqeSppvZuLzvpOwzUuB2wkt6MKqFOTu/zCzJYRhlucQsmIp8DZhf0jpfNeb2SzCUcWXErqHlgATCCOzJI3Yt3fOi+yQaLXNAn7u7v+Muh6pPjOrT9i5+qK7Xxh1PRINBb7sJDHC53BCt8vhhL7cTdFWJbvDzP5C6Iv/AvgO4ZdNZ+B77v5xlLVJdNSlI+W5iNBtMBc4R2GfkfIJ3ThNCF1o7wG9Ffbxpha+iEhM6EhbEZGYSKsunQMOOMALCwujLkNEJKN88MEHK929cWXzpVXgFxYWMnXq1KjLSLkpU2D0aJg7F1q3hqFDoXv3qKsSkWxhZourMp+6dGrZsGHQqxeMHQvTpsETT4Tnwyo7sa+ISIop8GvRlCkwYgRs3Ah5vpnvsJSSkvB8xIgwXURkT1Hg16LRo2FTYkDjRI5hTNJpXIqLw3QRkT0lrfrws83cuVA66nUSR/ML/kY+myhmb0pKYJ4OPJcMtmXLFoqKiiguLo66lNjIz8+nefPm1KlTp1rLK/BrUevWMH06lJTAOPrwP9xFT97h/+hFTk6YLpKpioqKaNCgAYWFheziOjSSIu7OqlWrKCoqomXLltV6D3Xp1KKhQyE/cZmNtziWLeTRh3AW2fx8+NWvIixOpIaKi4vZf//9FfZ7iJmx//771+gXlQK/FnXvDldeCQUFsCmnPpPpQW/GU1AQXtfQTMl0Cvs9q6aftwK/lt18M0yYAGedBTOb9ua7TOPNp1dx881RVyYicaPA3wO6d4fHHoNLnupDDk7X9eVeE0Ikq02ZAgMHQpcu4T4Vw5Jzc3Pp1KkTRx55JB07duRPf/oTJSW7vsbLokWLePTRR2u+8kpcdNFFzJy56wuhPfvss5XOk0oK/D3pe9+DBg1gfEXXhhbJTrV1AOLee+/N9OnT+fTTTxk3bhwvv/wyN1fy83lPBf59991Hu3btdjnPng583D1tbl26dPGsd+qp7oceGnUVIjU2c+bMKs03ebJ7QYF7GKT87VtBQZheXfXq1fvW8wULFvh+++3nJSUlvnDhQj/mmGO8c+fO3rlzZ580aZK7u3fr1s332Wcf79ixo48cObLC+ZItXLjQ27Rp4+ecc44fccQR/tOf/tQ3bNjg7u7jx4/3Tp06efv27f3888/34uJid3c/7rjj/P33399e5/XXX+9HHXWUd+vWzZcvX+6TJk3yfffd1wsLC71jx44+f/58Hz16tLdt29Y7dOjgZ599drnbXN7nDkz1KmRs5CGffItF4P/5z+FjX7Ag6kpEaqSqgT9ggLtZ+YGfkxOmV1fZwHd3b9iwoS9fvtw3bNjgmzZtcnf3uXPnemm+vP76637yySdvn7+i+ZItXLjQAZ84caK7u59//vl+5513+qZNm7x58+Y+Z84cd3c/77zzfNSoUe7+7cAH/Pnnn3d396uuusr/8Ic/uLv7oEGD/Mknn9y+nqZNm27/wlizZk2521yTwFeXzp7Wp0+4nzAh2jpE9pDkAxDLqs0DELds2cLPf/5zOnToQP/+/SvsOqnqfC1atODoo48G4Nxzz2XixInMmTOHli1b0jpxUM2gQYN46623dlp2r7324pRTTgGgS5cuLFq0qNx1HHXUUfzsZz/j4YcfJi8v9YdJKfD3tDZtoFkzGDcu6kpE9ojWrSGngqRJ9QGIn332Gbm5uRx44IGMGjWKJk2a8NFHHzF16lQ2b95c7jJVna/skMjdGSJZp06d7fPn5uaydevWcud78cUXufTSS5k2bRrf+973KpyvuhT4e5oZnHhiCPwU/zFF0lHyAYhlpfIAxBUrVnDJJZdw2WWXYWasW7eOpk2bkpOTw0MPPcS2bdsAaNCgAevXr9++XEXzlfX5558zefJkAB599FGOOeYY2rRpw6JFi5g/fz4ADz30EMcdd1yVa06upaSkhCVLlnD88cdz++23s27dOr7++utqfRYVUeBH4aSTYO1aePfdqCsRqXXJByCWtvRzckjJAYibNm3aPiyzd+/enHjiiQxLDP355S9/yQMPPEDHjh2ZPXs29erVA0K3SW5uLh07dmTUqFEVzldWmzZtuOeee2jbti1r1qzhF7/4Bfn5+dx///3079+fDh06kJOTwyWXXFLl+gcMGMCdd95J586dmTdvHueeey4dOnSgc+fO/OpXv6JRo0bV/3DKkVbXtO3atatn4wVQdrJ2LRxwAFxzDdx6a9TViFTLrFmzaNu2bZXnL70Q0Lx50KpVZl0IaNGiRZxyyinMmDEj6lLK/dzN7AN371rZsjXeK2BmLYAHgSaAA/e6+2gz2w8YCxQCi4Cz3H1NTdeXFRo1gqOPhpdeUuBLbHTvnjkBn61S0aWzFfiNu7cDugOXmlk74Fpggru3AiYknkupfv3CqTSXLo26EhGpRGFhYVq07muqxoHv7svcfVri8XpgFtAMOA14IDHbA8DpNV1XVjnppHD/yivR1iEisZHSnbZmVgh0Bt4Fmrj7ssSk5YQun/KWGWJmU81s6ooVK1JZTnrr0CEMz3z55agrEZGYSFngm1l94CngCnf/Knla6ZFm5S3n7ve6e1d379q4ceNUlZP+zEK3zmuvwZYtUVcjIjGQksA3szqEsH/E3Z9OvPxfM2uamN4U+DIV68oq/frB+vUwaVLUlYhIDNQ48C0cPvZPYJa7j0ya9DwwKPF4EPBcTdeVdXr1gjp1wmgdEdlt9evXB+CLL77gzDPPTMl7zp49m06dOtG5c2cWLFhAz549gT13ls3alIoW/tHAecAJZjY9cesHDAf6mNk8oHfiuSRr0ACOPRZeeCHqSkQy2ne+8x3+/e9/p+S9nn32Wc4880w+/PBDDjvsMN555x1AgQ+Au090d3P3o9y9U+L2kruvcvde7t7K3Xu7++pUFJx1Tj0VZs0KZ5gSkWpZtGgR7du3B2DMmDGcccYZ9O3bl1atWnH11Vdvn++1116jR48efPe736V///47nbrgpZde4q677uJvf/sbxx9/PLDjV8S1117L22+/TadOnRg1atQe2rLUSv3p2GT3nHZaOOTwuefgqquirkakeq64IhxXkkqdOsFdd1Vr0enTp/Phhx9St25d2rRpw+WXX87ee+/NLbfcwvjx46lXrx633347I0eO5MYbb9y+XL9+/bjkkkuoX78+V1555bfec/jw4YwYMYIXMvgXuQI/aoccAp07w7PPKvBFUqRXr140bNgQgHbt2rF48WLWrl3LzJkzt5/iePPmzfTo0SPKMvc4BX46OP10uOkmWL4cDjoo6mpEdl81W+K1pW7dutsfl56O2N3p06cPjz32WISVRUtny0wHp58erhDxn/9EXYlI1urevTuTJk3afirjDRs2MHc39p2VPa1yJlLgp4MOHeDQQ0O3jojUisaNGzNmzBgGDhzIUUcdRY8ePZg9e3aVly97WuVMpNMjp4vf/AbuvhtWrgzDNUXS3O6eHllSoyanR1YLP12cfjps3qyTqYlIrVHgp4uePcNFUZ55JupKRCRLKfDTRW5uGJP/wgtQXBx1NSJVkk5dwnFQ089bgZ9OzjornExNp0yWDJCfn8+qVasU+nuIu7Nq1SryK7oifBVoHH46OeGE0K0zdiz85CdRVyOyS82bN6eoqIhYXcciYvn5+TRv3rzayyvw00leHvz0p/DQQ7BhA9SrF3VFIhWqU6cOLVu2jLoM2Q3q0kk3Z58NGzfqlMkiknIK/HRz7LHh9Apjx0ZdiYhkGQV+usnNhTPPhBdfDDtwRURSRIGfjs4+OwzN1Ll1RCSFFPjpqGdPaNYMHn886kpEJIso8NNRTg4MGBBOs7ByZdTViEiWUOCnq//3/2DLFojxubtFJLUU+OnqqKPCJd4efDDqSkQkSyjw09mgQTB1KsycGXUlIpIFFPjpbODAMExTrXwRSQEFfjpr0gROOimcamHbtqirEZEMp8BPd4MGwRdfwIQJUVciIhlOgZ/ufvxjaNQIHngg6kpEJMMp8NNd3bpwzjnw9NOwenXU1YhIBlPgZ4IhQ8KpFrTzVkRqQIGfCTp2hG7d4B//AF1dSESqSYGfKS6+GGbPhrffjroSEclQCvxMcfbZ0LBhaOWLiFSDAj9TFBTAeefBv/+tE6qJSLUo8DPJxRfD5s0wZkzUlYhIBlLgZ5L27eEHP4B77tGRtyKy2xT4meaKK2DRInj++agrEZEMk5LAN7N/mdmXZjYj6bX9zGycmc1L3O+binXF3mmnQWEhjBoVdSUikmFS1cIfA/Qt89q1wAR3bwVMSDyXmsrNhcsvD8MzP/gg6mpEJIOkJPDd/S2g7HH/pwGlJ4B5ADg9FesS4MILoX59GD066kpEJIPUZh9+E3dflni8HGhSi+uKl4YN4YILwkXOly2rfH4REfbQTlt3d6DccwKY2RAzm2pmU1esWLEnyskOv/pVGKmjVr6IVFFtBv5/zawpQOL+y/Jmcvd73b2ru3dt3LhxLZaTZQ47DPr3h7/+FdasiboaEckAtRn4zwODEo8HAc/V4rri6brrYP16uPvuqCsRkQyQqmGZjwGTgTZmVmRmFwLDgT5mNg/onXguqdSxI5xyCtx1F3z9ddTViEiay0vFm7j7wAom9UrF+8suXH899OwJ994Lv/511NWISBrTkbaZrkcPOP54GDEiXCRFRKQCCvxs8NvfhuGZOnWyiOyCAj8bnHBCuN16q/ryRaRCCvxsceutsGKFxuWLSIUU+Nmie3c49VS4805YXfYsFyIiCvzscsst8NVXLL3iTgYOhC5dYOBAmDIl6sJEJB0o8LNJhw581P4c9nvoLt55/HOmTYMnnoBevWDYsKiLE8l+U6aQ1o0tBX4WmTIFzpp3G44xnGsAKCmBjRvDqM10+8cnkk2GDQuNq7FjSdvGlgI/i4weDfO+OZg7uYqBPM7RTNw+rbhY+3NFasuUKaFRtXEjeOI0kenY2FLgZ5G5c8M/tju4miKaMZqhGCVA+Mc3b17EBYpkqdGjYdOm8Pg8HuR3/J5ctgLp1dhS4GeR1q0hJwc2Uo+ruYMuTON87gfC661bR1ygSJYqbWztxypG8mt6M55t5ALp1dhS4GeRoUMhPz88foyBTORobucaDmAF+fnhFPoiknqlja0/ch2NWMul3AMYkF6NLQV+FuneHa68EgoKICfHuJh/sA9f8efcX3PllWG6iKTe0KHwwzqTGML/chdXMIMO26elU2NLgZ9lbr4ZJkyAs86CvbscyQvtrmHgtoe5+ejXoi5NJGt177iJsfUvYLEdwu/tJiC07AsKSKvGVkpOjyzppXv3pH9gxTdApyfhkkvgk0+gXr1IaxPJSjffzAGr5vLl6HGcPLk+8+ZBq1ah5Z8uYQ9g7uVeajYSXbt29alTp0ZdRvZ56y047ji4/HL485+jrkYku0ydCt26wfnnw333RVKCmX3g7l0rm09dOnFw7LGhE/Evf4FXX426GpHsUVwcgv6gg8KA+zSnwI+L4cPhyCNh8GBYuTLqakSyw1VXwYwZoWXfqFHU1VRKgR8Xe+8NjzwSzqT585/vOBxQRKrn+efh7rvhf/4HTjop6mqqRIEfJx07hvPmP/ts+IcqItWzdGnoyuncGf74x6irqTIFftz8+tfw4x+H+0mToq5GJPNs3gwDBsA338Djj0PdulFXVGUK/LjJyYEHH4TCQujfH5Yvj7oikcwydChMnBj67dPlENoqUuDHUaNG8PTTsG4dnHlmaKmISOXuvRf+/ne4+urQys8wCvy46tAB/vWv0K0zeHA4w5OIVOyNN+Cyy+BHP4Lbbou6mmrRkbZxdvbZsGgRXHstHHJIGLopIjv7+GM47bRw+Oxjj0FubtQVVYsCP+6uvhoWL4bbb4cWLeDSS6OuSCS9LF4MfftCgwbwyiuw775RV1RtCvy4MwunW/jii/Bzda+9wjh9EQnDL/v0CZeumjgxNIoymPrwBfLywoU4+/WDIUNC375I3C1dCscfH0ayvfQStG8fdUU1psCXoG5deOqp8NP1oovCSASRuPr88x1h/8or0LNn1BWlhAJfdsjPh2eeCS39X/wCbrxRp2CQ+Pnoo3BO4y+/zKqwBwW+lJWfH069cMEF8Ic/wIUXhiMLReJg/Hj4wQ/CKJyJE7Mq7EGBL+XJywtHEd54I9x/f/hp+8UXUVcVG1OmwMCB0KVLuJ8yJeqKYsAdRo0KXZqFhTB5clb02ZelwJfymYXrJT72WPiJ27lzOPBEatWwYdCrV9iHPm0aPPFEeD5sWNSVZbH168MxKaXnmXr7bWjePOqqaoUCX3ZtwAB4770w9rhXL7jmmnDRB0m5KVPCNTQ2btyx66SkJDwfMUIt/Vrx9tuhMfPUU3DHHeGUIw0bRl1VrdE4fKlcu3bw/vvhasx33AH/+U/o6unWDQhBNHo0zJ0bziWVbtfxrKmUbp97+MJcuzacy2jduvB8yxbe/MNm+mzcQh5byGMr28hlC3XYQh22barD+Ovz6H5bAeyzTwilhg3DNYrNUrm5WaXCv93GjfDb38Jdd4UunNdfD1eGy3K1fk1bM+sLjAZygfvcvcLj93VN2wzw6qth2GZREQwezJ373sZN/2jKpk0hy3Jywn7fK68MPUKZbtiw0LqudPvc4b//DclSVATLlu18W7UqBP2WLakrMDc3fAE0bhwus3fQQdC06Y7HLVpAy5bhvk6dnRbP5i/rcv92dZ0xJz1O//evgSVLwpHlw4dD/fpRl1sjVb2mba0GvpnlAnOBPkAR8D4w0N1nlje/Aj9DfPUV3HYbJSNHsWHLXozk14xmKGvYb/ssBQUwYUJmh8eUKaEXa+PGHa/VYTOHM5+j9prDHRfM5uCNs2HOHJg9O7TWk+Xnh/AtvTVuHFrljRp9+z4/H/bai9/+vg6vjK9Dse/FVvLIY2uifb+FuraFE4/fwu+u3LTjl8G6deFvsXYtrFgRxowvXx6+XNav/3Ytubk7wr9lSzj0UMZOb82IF9rySXErvqFuVn1Z7/y3c47ndW7ht/RkMhtad6befaPDiJwskC6B3wO4yd1/lHh+HYC7l3uJGAV+Zhl6ygJ+8OI1nMlTrKc+93Ap93ApRbQgJwfOOivs881IK1dy04DZFE2YQxtmcwSzacMcDuUz8ti2Y75mzaBNGzjiiHBr0yYEa9OmIcx3o7ulvC+YUrv9BbphQwj/xYth4cKdb0nXQdhGDp9xaGIrj+Czvdpy6V+OoH3/thl73piBA8OOb/NtnMTLXM9t9GQyX9CUG+0WNvYfxKNjM/MEaOVJl8A/E+jr7hclnp8HdHP3y5LmGQIMATj44IO7LF68uNbqkdTq0iWMJDmSGdzArZzNWBzjJfpxL0NY2flEJk9L/dWAUtYNsXUrfPZZaJ2XttJLH69atX22Yuoyl9bsiP4jmEMbCjq14a0PG6Ruw9jRDVFcHHbY1lare3D/DXz073kcwSzaMosjmE1bZtGaudQl6biLJk3CPpy2bcOt9PFBB1Vr38Ge6kI6qUMRXWaM4SLuo5DFLOZghnMt93M+35BPly6QTW3LqgY+7l5rN+BMQr996fPzgLsrmr9Lly4umWPAAPecHPfQQ+p+CAv9Fq73ZTRxB99QZx/3c85xf/JJ95UrU7LOG290LyhwNwvrzMkJz2+8sYIFSkrcV6xwf+cd9zFj3K+91v30093btnWvU2dH8eDepIn7sce6DxniPnKk//HYl/wwW+A5bP3WbKXrHTgwJZu0k8mTw2fbpUu4nzw59ev47nd9p20C9xy2+mHM86GH/cf99tvdBw9279bNfZ99vj1jo0buPXq4X3ih+4gR7i++6L5wofu2bRWuc7f/drujpMR91iz3O+5w7959e52v0dv7M9bz2LxH/nZRAaZ6FTJZXTpSbRV1QeSxhR/XHcff+zzFgZOf29Fabt8+jITo2DFcgKV9+3DK2Rqury7FHJ6/lLF/KuLIfZbAggUwb15oRs6bF/q4txeXF85pXrYbpk2bnbovUtrFkmYGDgxj/Mu77k1OThiW/uijSS+6h4PvZs2CmTPDfenjFSt2zFdQED7T0l8EhYXQogUfrjqYH/6sGV9t2nnHcbU+y//+F2bMgE8+gXfegbfeCq8BdOnC513P4OQHzmJG8eGpWV+aS5cunTzCTttewFLCTttz3P3T8uZX4GeeSrsgtm4Nyfnmm+E2eTJ8/fWON9h//3CQS/PmofugXr0wYqJevbCjcdu2cNu6lRee2MiXc9bQkLU0Yi37s4pmLKUxK79dlFnoR2/dOoR76a1167DDspzRKtXevgyV0i+zlSt3fAEkfyEsWfKt2UowvuA7LKUZq9if1ezHavZjDftxWKf6nDsoL3wh5yVGi2/cGPZFbNgAq1eHs1cuXRred/XqHW/cogUcd1xoTPTpE75kyN6/XXnSIvAThfQD7iIMy/yXu99a0bwK/MxU2i87b17I1V32y5aUhB2Jn3wSWmhLloRhjEuXhpNVbdgQvhC2bt1p0U05Bawq2Ze1NGIN+7KGfVlKM4poThHNqdeqGX99vnn4D5+fH832ZZBaD8QNG8Lf9/PP+cPPP6fk8yUczOc0Yyn7sob9WM2+rKERa8lhFzm0115hRFOzZuHWvHn4FdG+ffileOCBFS6arX+7stIm8HeHAl+2++abkEK5ueGWk8PAc2z3uiGkUnsqEHfVhZRn2zj3p5u4/3+3hi/6rVvDjPXqhZ8bu/GLLK4U+JJ1srlPPdvpb1e7qhr4OpeOZIzu3UN3Q0FBaNFDuC8oCK8rMNKX/nbpQS18yThx6ZfNRvrb1Q516YiIxIS6dERE5FsU+CIiMaHAFxGJCQW+iEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmFPgiIjGhwBcRiQkFvohITCjwRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhMKPBFRGJCgS8iEhM1Cnwz629mn5pZiZl1LTPtOjObb2ZzzOxHNStTRERqKq+Gy88AzgD+kfyimbUDBgBHAt8BxptZa3ffVsP1iYhINdWohe/us9x9TjmTTgMed/dv3H0hMB/4fk3WJSIiNVNbffjNgCVJz4sSr4mISEQq7dIxs/HAQeVMusHdn6tpAWY2BBgCcPDBB9f07UREpAKVBr67967G+y4FWiQ9b554rbz3vxe4F6Br165ejXWJiEgV1FaXzvPAADOra2YtgVbAe7W0LhERqYKaDsv8iZkVAT2AF83sVQB3/xR4ApgJvAJcqhE6IiLRqtGwTHd/Bnimgmm3ArfW5P1FRCR1dKStiEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmFPgiIjGhwBcRiQkFvohITCjwRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmFPgiIjGhwBcRiQkFvohITCjwRURiQoEvIhITNQp8M7vTzGab2cdm9oyZNUqadp2ZzTezOWb2oxpXKiIiNVLTFv44oL27HwXMBa4DMLN2wADgSKAv8Fczy63hukREpAZqFPju/pq7b008nQI0Tzw+DXjc3b9x94XAfOD7NVmXiIjUTCr78C8AXk48bgYsSZpWlHhtJ2Y2xMymmtnUFStWpLAcERFJllfZDGY2HjionEk3uPtziXluALYCj+xuAe5+L3AvQNeuXX13lxcRkaqpNPDdvfeuppvZYOAUoJe7lwb2UqBF0mzNE6+JiEhEajpKpy9wNXCqu29MmvQ8MMDM6ppZS6AV8F5N1iUiIjVTaQu/EncDdYFxZgYwxd0vcfdPzewJYCahq+dSd99Ww3WJiEgN1Cjw3f3wXUy7Fbi1Ju8vIiKpoyNtRURiQoEvIhITCnwRkZhQ4IuIxIQCX0QkJhT4IiIxocAXEYkJBb6ISEwo8EVEYkKBLyISEwp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhM2I7L0EbPzNYDc6KuoxYdAKyMuohapO3LXNm8bZD929fG3RtUNlNNL3GYanPcvWvURdQWM5uq7ctc2bx92bxtEI/tq8p86tIREYkJBb6ISEykW+DfG3UBtUzbl9myefuyedtA2wek2U5bERGpPenWwhcRkVqiwBcRiYm0C3wz629mn5pZiZllxTAqM+trZnPMbL6ZXRt1PalmZv8ysy/NbEbUtaSambUws9fNbGbi3+XQqGtKJTPLN7P3zOyjxPbdHHVNtcHMcs3sQzN7IepaUs3MFpnZJ2Y2vbLhmWkX+MAM4AzgragLSQUzywXuAU4C2gEDzaxdtFWl3Bigb9RF1JKtwG/cvR3QHbg0y/5+3wAnuHtHoBPQ18y6R1tSrRgKzIq6iFp0vLt3quxYg7QLfHef5e7ZdLTt94H57v6Zu28GHgdOi7imlHL3t4DVUddRG9x9mbtPSzxeTwiNZtFWlToefJ14Widxy6qRHGbWHDgZuC/qWqKWdoGfhZoBS5KeF5FFgREnZlYIdAbejbiUlEp0d0wHvgTGuXtWbR9wF3A1UBJxHbXFgdfM7AMzG7KrGSM5tYKZjQcOKmfSDe7+3J6uR6QyZlYfeAq4wt2/irqeVHL3bUAnM2sEPGNm7d09K/bHmNkpwJfu/oGZ/TDicmrLMe6+1MwOBMaZ2ezEr+6dRBL47t47ivVGZCnQIul588RrkiHMrA4h7B9x96ejrqe2uPtaM3udsD8mKwIfOBo41cz6AfnAPmb2sLufG3FdKePuSxP3X5rZM4Ru5HIDX106te99oJWZtTSzvYABwPMR1yRVZGYG/BOY5e4jo64n1cyscaJlj5ntDfQBZkdaVAq5+3Xu3tzdCwn/9/4vm8LezOqZWYPSx8CJ7OLLOu0C38x+YmZFQA/gRTN7NeqaasLdtwKXAa8Sdvg94e6fRltVapnZY8BkoI2ZFZnZhVHXlEJHA+cBJySGvU1PtBazRVPgdTP7mNA4GefuWTd0MYs1ASaa2UfAe8CL7v5KRTPr1AoiIjGRdi18ERGpHQp8EZGYUOCLiMSEAl9EJCYU+CIiMaHAFxGJCQW+iEhM/H/2Hw1dpgo4cgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "\"\"\"Removing points\"\"\"\n", "num_points = 10 ### number of points to keep, set it to 10 if you want to keep all points\n", "new_x, new_y = reduce_points(x_data, y_data, num_points)\n", "new_y = noisy_points(new_y, 0) ### how much extra noise do you want, set it to 0 if you do not want any noise\n", "degree = 4 # use this to change the degree of your polynomial\n", "D = data_matrix(new_x, degree)\n", "params = np.linalg.lstsq(D, new_y, rcond=None)[0]\n", "#plotting the figure\n", "fig=plt.figure()\n", "ax=fig.add_subplot(111,xlim=[-1,5],ylim=[-25,25])\n", "x,y=poly_curve(params, new_x) ###### Plotting the fitted data\n", "\n", "ax.plot(new_x,new_y,'.b',markersize=15) ###### Plotting the new data\n", "ax.plot(x,y,'r') \n", "ax.legend(['Data points','line fit'])\n", "plt.title('Polynomial of Degree %d' %(len(params)-1),fontsize=16)\n", "print(error(params, new_x, new_y))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.3" }, "vscode": { "interpreter": { "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6" } } }, "nbformat": 4, "nbformat_minor": 1 }