{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# EECS16A Discussion 13B" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 3: 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": "code", "execution_count": 1, "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 notebook\n", "\n", "\"\"\"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=[x_input[1], x_input[-1]]\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", " \n", "\"\"\"Function that defines a data matrix for some input data. 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", " \n", "np.random.seed(10)\n", "\n", "# Parameters corresponding to the polynomial we want: $y = 1600 + 40x -102x^2 -x^3 + x^4$\n", "smarap=np.array([1600,40,-102,-1,1])/100.0\n", "x_data=np.linspace(-11,11,50)\n", "y_data=np.dot(data_matrix(x_data,4),smarap)+(np.random.rand(len(x_data))-0.5)*10\n", "\n", "y_data_training = y_data[0::2]\n", "x_data_training = x_data[0::2]\n", "y_data_test = y_data[1::2]\n", "x_data_test = x_data[1::2]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Input and output data\n", "\n", "`x_data` is the input (set of $x_i$'s), and `y_data` is the output (set of $y_i$'s). Here we plot our test points. Does it look like $y$ should be a quartic polynomial in $x$?" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "50\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAegAAADCCAYAAACCEFHoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAW+klEQVR4nO3df5AU5Z3H8c+XFd0USykiEOLCQaoMAW/5ufH01NOcHDHGiklKUuTUI5KEwh+Vyx/nHR5V+XFWquJRHnVWzKU4Y0JMjMZLVCrqCfG8MpRoXCiQXyKrbuJSiAQ5AhEU5Xt/dA+ZDDO7DdM983TP+1W1NTPdvdPPMz09336efn6YuwsAAIRlSLMTAAAAjkeABgAgQARoAAACRIAGACBABGgAAAJEgAYAIECn1PsGZjZO0g8ljZHkkpa7+7+b2ZmSHpA0QVKfpM+6+76B3uuss87yCRMm1JskAAByY926db9z91GVy63eftBmNlbSWHdfb2bDJa2T9ClJn5f0prt/y8wWSxrh7v800Ht1d3d7T09PXekBACBPzGydu3dXLq+7itvdd7n7+vj5AUnbJJ0t6SpJK+LNVigK2gAAIIFU70Gb2QRJMyQ9J2mMu++KV72uqAocAAAkkFqANrMOST+T9BV3/335Oo/q0avWpZvZQjPrMbOePXv2pJUcAAByLZUAbWZDFQXnH7v7z+PFu+P706X71G9U+193X+7u3e7ePWrUcffIAQBoSXUHaDMzSd+TtM3d/61s1UpJ8+Pn8yU9Uu++AABoFXV3s5J0oaTrJG0ysw3xsn+W9C1JPzWzL0j6jaTPprAvAABaQt0B2t3XSLIaqy+r9/0BAGhFjCQGAECACNAAAASIAA0AQIAI0AAABIgADRTIwYPS2rXRY5H2BeRNGudHGt2sAATg4EFpyhRp3z5pxAhp61apoyP/+wLyJq3zgxI0UBCbNkU/CAcPRo+bNhVjX0DepHV+UIIGCqKrK7pal6LHrq5i7AvIm7TODwI0UBAdHVFV2qZN0Q9CllXOjdwXkDdpnR8EaKBAOjqkCy4o3r6AvEnj/GjZe9C0QAXCwLkIVNeSAbrUwm7OnOiRHwa0kiQBMa2gOdj7cC4CtbVkFXd5C7vSa6rq0AqSdP9Iq4tIkvfhXARqa8kSdKmFXUcHLVDRWpJ0/0iri0iS9+FcBGpryRI0LVDRqpJ0/0iri0iS9+FcBGozd292Go7p7u72np6eZicDKLSDBwcPiEm2SWtfQKszs3Xu3l25vCVL0EArS9L9I60uVI3uisUFAYqkJe9BA8iXpC3PaRGORmhU10BK0ACClrRVOS3C0QiNnCiGEjSAoCVtVU6LcDRCIyeKoQQN5ESr3l9N2qqcFuFohEZOFEOABnKgledfPpHAy/jgyFojLwSp4gZyoNXnXy4F3la5KEHYGvV9pAQN5ADzLwOthwAN5AD3V4HWQxU3kKE0+0tSzQu0FgI0UEUaUzIycAaAelDFDVRIa0pGBs4AUA9K0ECFtKZkZOCM/GrUUI7AQChBAxXSmpKRhl351Mp9zhGWVErQZnaPmb1hZpvLlp1pZqvNbEf8OCKNfQFZKwXWVatq/zgn2aa0HQ278qXV+5wjHGlVcf9A0uUVyxZLetLdz5H0ZPwayIUkgZXgW0zcmkAoUgnQ7v60pDcrFl8laUX8fIWkT6WxLwDIUtLaESBrWd6DHuPuu+Lnr0saU20jM1soaaEkjR8/PsPkAEAyjOmNEDSkFbe7uySvsW65u3e7e/eoUaMakRwgOLQaBlApywC928zGSlL8+EaG+wJyiwFNGouLIeRFlgF6paT58fP5kh7JcF9AbtFquHG4GEKepNXN6ieS1kqaZGb9ZvYFSd+S9DdmtkPS7Pg1gAq0Gm4cLoaQJ6k0EnP3z9VYdVka7w8UGQOaNA7TdiJPGEkMCACthhuDiyHkCQEaQEvhYgh5wWQZAAAEiAANAECACNAAAASIAA0AQIAI0AAABIgADQBAgAjQAAAEiAANAECACNADYNYbAECzMJJYDaVZb/bti8bs3bqVYQEBAI1DCboGZr0BADQTJegamPUGANBMlKBrKM16s2oV1dt5QrsBAEVBCXoAzHqTL7QbAFAklKBRGLQbAFAklKBRGLQbAFAkBGgURqndwKZNUXCmehtAnhGgUSi0GwBQFNyDBgAgQARotBy6YgHIA6q40VLoioW0HDxIewdkixI0WgpdsZCG0oXenDnRI7UxyAIBGi2l1BWro4OuWDh5XOihEajiRkuhKxbSQJ97NAIBGi2HrlioFxd6aAQCNACcBC70kDXuQSM36B4FoJUQoBGEwYIvrWYBNEJIBYHMA7SZXW5m282s18wWZ70/5E+S4EurWQBZC60gkGmANrM2SXdJ+rikKZI+Z2ZTstwn8idJ8KV7FICshVYQyLoEfZ6kXnd/xd3fkXS/pKsy3idyJknwLbWaXbWK0b8AnJzBqq9DKwhk3Yr7bEmvlb3ul/QX5RuY2UJJCyVp/PjxGScHIUraZYVWswBOVpJhfkPrPtf0RmLuvtzdu929e9SoUc1ODpqkFHybfUIAKKak1dch/RZlXYLeKWlc2evOeBkAAA2Tx9Hfsg7Qz0s6x8wmKgrM8yT9bcb7BADgT4RWfZ1EpgHa3d81s5slPSGpTdI97r4ly30CAFBN3tqxZD7Up7s/JumxrPcDAECRNL2RGAAAOB4BGgCAABGgASAjjRzXOaQxpJEOppsEgAwkGRijtF29LYuT7gv5Qgm6Tly1AqgmycAYaU3OENoY0kgHJeg6cNUKoJYkA2OUB9bS65PpBpTHQTgwOErQdSjyVSs1A0B9kkzwktbkDEwmU0yUoOtQ1KtWagaAdAw2MEaao1vlbRAODI4SdB2KetVa5JoBIDQhTc6AsFCCrlMRr1qLWjMAAHlCCToQId3zTbtmIKS8AUBeUIIOQIj3fNOqGQgxbwCQB5SgA1Dke75FzhsAZIkSdACKfM+3yHkDgCwRoAOQx4nEkypy3gA0RhrDoeYRAToQRWwNXlLkvAHIViu3Y+EeNAAgWK3cjoUSNAAgWK3cjoUADQAIViu3YyFAAwCC1qrtWLgH3YIY2QvIF87Z1kQJugFC6iLQyi0igTzinG1dlKAzVjq55syJHpt9BdzKLSKBPMrzOUvJvz6UoDNWfnKVXjfzXkort4gE8iiv5ywl//pRgs5Y6eTq6Ajj5CrqHNZAUTX6nE2r1Nvokn8RS+uUoDMWYheBVm0RCeRVo87ZNEu9jSz5F7W0Tgm6AUonVxG+MACKK81SbyNL/nm+Tz8QStAAAEnpl3obVfLP6336wRCgAQCSwrwll0Re0z0YAjQA4Ji8tlHJa7oHUtc9aDOba2ZbzOyomXVXrLvVzHrNbLuZfay+ZAIA0FrqbSS2WdJnJD1dvtDMpkiaJ+lcSZdL+o6ZtdW5LwBAIIrYrSk0dVVxu/s2STKzylVXSbrf3d+W9KqZ9Uo6T9LaevYHAGi+onZrCk1W3azOlvRa2ev+eNlxzGyhmfWYWc+ePXsySg4AIC1F7dYUmkFL0Gb2S0nvr7Jqibs/Um8C3H25pOWS1N3d7fW+HwAgW0Xt1hSaQQO0u88+iffdKWlc2evOeBkAIOeK2q0pNFlVca+UNM/MTjOziZLOkfTrjPYFAGgwRkjMXr3drD5tZv2SLpD0qJk9IUnuvkXSTyVtlfTfkm5y9/fqTSzCQ0tOAMhGva24H5L0UI1135T0zXreH2GjJScAZIfJMnIktNIqLTkBIDsM9ZkTIZZWackJ5MvBgzTsyhMCdE6Ul1ZLr5s97iwtOYH8CPEiHwOjijsnSqXVjo6wSqu05ATygVtS+UMJOicorQKoB7ek8ocAnSNFnE4NQGNwkZ8/VHEXTGgtvQGEI8RbUvxm1UYJukBoBAIgT/jNGhgl6AKhEQiAPOE3a2CUoAuERiAA8oTfrIERoAuERiAA8oTfrIERoAuGlt4A8oTfrNq4Bw0AQIAI0AAABCj4Ku4jR46ov79fhw8fbnZSWkZ7e7s6Ozs1dOjQZicFAFpW8AG6v79fw4cP14QJE2RmzU5O4bm79u7dq/7+fk2cOLHZyQGAlhV8Fffhw4c1cuRIgnODmJlGjhxJjQUANFnwAVoSwbnB+LwBoPlyEaCbra2tTdOnT9e5556radOm6Y477tDRo0cH/J++vj7dd999mafti1/8orZu3TrgNg8//PCg2wAAwkKATuB973ufNmzYoC1btmj16tV6/PHH9Y1vfGPA/2lUgL777rs1ZcqUAbchQANA/hQyQGc5O8ro0aO1fPlyffvb35a7q6+vTxdffLFmzpypmTNn6plnnpEkLV68WL/61a80ffp0LVu2rOZ25fr6+vThD39Y11xzjSZPnqyrr75ab731liTpySef1IwZM9TV1aUFCxbo7bffliRdeuml6unpkSR1dHRoyZIlmjZtms4//3zt3r1bzzzzjFauXKlbbrlF06dP18svv6w777xTU6ZM0dSpUzVv3rz0PyQAQP3cPZi/WbNmeaWtW7cet2wgBw64jxvn3tERPR44cEL/XtWwYcOOW3b66af766+/7n/4wx/80KFD7u7+0ksveSkPTz31lH/iE584tn2t7cq9+uqrLsnXrFnj7u7XX3+9L1261A8dOuSdnZ2+fft2d3e/7rrrfNmyZe7ufskll/jzzz/v7u6SfOXKle7ufsstt/htt93m7u7z58/3Bx988Nh+xo4d64cPH3Z393379lXN84l+7gCAkyOpx6vExMKVoBs9O8qRI0f0pS99SV1dXZo7d27NquSk240bN04XXnihJOnaa6/VmjVrtH37dk2cOFEf+tCHJEnz58/X008/fdz/nnrqqbryyislSbNmzVJfX1/VfUydOlXXXHONfvSjH+mUU4LvaQcALalwAbo0O0pHR3azo7zyyitqa2vT6NGjtWzZMo0ZM0YbN25UT0+P3nnnnar/k3S7yhbUJ9KieujQoce2b2tr07vvvlt1u0cffVQ33XST1q9fr4985CM1twMANE/hAnRpdpRVq7KZ/HvPnj1atGiRbr75ZpmZ9u/fr7Fjx2rIkCG699579d5770mShg8frgMHDhz7v1rbVfrtb3+rtWvXSpLuu+8+XXTRRZo0aZL6+vrU29srSbr33nt1ySWXJE5zeVqOHj2q1157TR/96Ed1++23a//+/TqYxc16AEBdChegpT/OjpJWcD506NCxblazZ8/WnDlz9LWvfU2SdOONN2rFihWaNm2aXnzxRQ0bNkxSVI3c1tamadOmadmyZTW3qzRp0iTdddddmjx5svbt26cbbrhB7e3t+v73v6+5c+eqq6tLQ4YM0aJFixKnf968eVq6dKlmzJihHTt26Nprr1VXV5dmzJihL3/5yzrjjDPq/owAAOmy6P50GLq7u73UIrlk27Ztmjx5cpNS1Fh9fX268sortXnz5mYnpaU+dwBoJjNb5+7dlcsLWYIGACDvCNABmTBhQhClZwBA8xGgAQAIUF0B2syWmtmLZvaCmT1kZmeUrbvVzHrNbLuZfaye/YR0n7wV8HkDQPPVW4JeLenP3X2qpJck3SpJZjZF0jxJ50q6XNJ3zKztZHbQ3t6uvXv3EjQaxOP5oNvb25udFABoaXUNI+Xuq8pePivp6vj5VZLud/e3Jb1qZr2SzpO09kT30dnZqf7+fu3Zs6eepOIEtLe3q7Ozs9nJAICWluY4jwskPRA/P1tRwC7pj5cdx8wWSlooSePHjz9u/dChQzVx4sQUkwkAQPgGDdBm9ktJ76+yaom7PxJvs0TSu5J+fKIJcPflkpZLUT/oE/1/AACKaNAA7e6zB1pvZp+XdKWky/yPN4p3ShpXtllnvAwAACRQbyvuyyX9o6RPuvtbZatWSppnZqeZ2URJ50j6dT37AgCgldQ11Gfc+Os0SXvjRc+6+6J43RJF96XflfQVd388wfvtkfSbk07Q8c6S9LsU36+ZyEuYipKXouRDIi+hKkpessjHn7n7qMqFQY3FnTYz66k2vmkekZcwFSUvRcmHRF5CVZS8NDIfjCQGAECACNAAAASo6AF6ebMTkCLyEqai5KUo+ZDIS6iKkpeG5aPQ96ABAMiropegAQDIpdwHaDOba2ZbzOyomXVXrBt0Ri0zm2hmz8XbPWBmpzYm5QOL07Ih/uszsw01tuszs03xdj0NTmYiZvZ1M9tZlp8ramx3eXyses1scaPTmcRAM7hVbBfkcRnsM47HLnggXv+cmU1oQjIHZWbjzOwpM9san/9/X2WbS81sf9n37qvNSGsSg31fLHJnfFxeMLOZzUjnQMxsUtlnvcHMfm9mX6nYJthjYmb3mNkbZra5bNmZZrbazHbEjyNq/O/8eJsdZjY/tUS5e67/JE2WNEnS/0rqLls+RdJGRf20J0p6WVJblf//qaR58fPvSrqh2XmqksY7JH21xro+SWc1O42DpP/rkv5hkG3a4mP0QUmnxsduSrPTXiWdcySdEj+/XdLteTkuST5jSTdK+m78fJ6kB5qd7hp5GStpZvx8uKLZ9CrzcqmkXzQ7rQnzM+D3RdIVkh6XZJLOl/Rcs9M8SH7aJL2uqH9vLo6JpL+SNFPS5rJl/yppcfx8cbXzXdKZkl6JH0fEz0ekkabcl6DdfZu7b6+y6tiMWu7+qqTSjFrHmJlJ+mtJ/xUvWiHpUxkm94TFafyspJ80Oy0ZO09Sr7u/4u7vSLpf0TEMiruvcvd345fPKhrGNi+SfMZXKToPpOi8uCz+DgbF3Xe5+/r4+QFJ21RjQp6CuErSDz3yrKQzzGxssxM1gMskvezuaQ48lSl3f1rSmxWLy8+HWvHhY5JWu/ub7r5P0TTMl6eRptwH6AGcLem1stfVZtQaKen/yn5wa8661UQXS9rt7jtqrHdJq8xsXTwzWKhujqvm7qlRTZTkeIVmgaJSTTUhHpckn/GxbeLzYr+i8yRYcTX8DEnPVVl9gZltNLPHzezcxqbshAz2fcnb+TFPtQsVeTkmkjTG3XfFz1+XNKbKNpkdmzSnm8yMJZhRK48S5utzGrj0fJG77zSz0ZJWm9mL8ZVgQw2UF0n/Iek2RT9Ctymqsl/QuNSdmCTHxQafwS2I41J0ZtYh6WeKhhP+fcXq9YqqWA/G7R4eVjQvQIgK832J2/F8UtKtVVbn6Zj8CXd3M2tot6dcBGgfZEatGpLMqLVXUVXRKXFpoaGzbg2WLzM7RdJnJM0a4D12xo9vmNlDiqoxG35iJz1GZvafkn5RZVUwM6AlOC6f1/EzuFW+RxDHpUKSz7i0TX/8/TtdfxxrPyhmNlRRcP6xu/+8cn15wHb3x8zsO2Z2lrsHNx50gu9LMOdHAh+XtN7dd1euyNMxie02s7Huviu+pfBGlW12Krq3XtKpqE1U3YpcxT3ojFrxj+tTkq6OF82XFFKJfLakF929v9pKMxtmZsNLzxU1YNpcbdtmqrhX9mlVT+Pzks6xqFX9qYqqyFY2In0nwmrP4Fa+TajHJclnvFLReSBF58X/1LoIaab4vvj3JG1z93+rsc37S/fPzew8Rb93wV1sJPy+rJT0d3Fr7vMl7S+reg1NzVq/vByTMuXnQ6348ISkOWY2Ir59NydeVr9mt5yr90/RD36/pLcl7Zb0RNm6JYparW6X9PGy5Y9J+kD8/IOKAnevpAclndbsPJWl8weSFlUs+4Ckx8rSvjH+26KoCrbp6a6Sj3slbZL0gqIv/NjKvMSvr1DUGvflgPPSq+h+04b4r9TiORfHpdpnLOlfFF1wSFJ7fB70xufFB5ud5hr5uEjRLZMXyo7FFZIWlc4ZSTfHn/9GRQ36/rLZ6a6Rl6rfl4q8mKS74uO2SWU9VkL6kzRMUcA9vWxZLo6JoouKXZKOxDHlC4raXzwpaYekX0o6M962W9LdZf+7ID5neiVdn1aaGEkMAIAAFbmKGwCA3CJAAwAQIAI0AAABIkADABAgAjQAAAEiQAMAECACNAAAASJAAwAQoP8HASg+1qXhFAQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# x_data and y_data have already been defined\n", "\n", "fig=plt.figure(figsize=(8,3))\n", "ax=fig.add_subplot(111,xlim=[-11,11],ylim=[-21,21])\n", "x,y=poly_curve(smarap,x_data)\n", "ax.plot(x_data,y_data,'.b',markersize=5)\n", "ax.legend(['Data points'])\n", "print(len(x_data))" ] }, { "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": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "degree=4\n", "DataMat=data_matrix(x_data,degree)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Least Squares\n", "Recall the system of equations we are trying to solve here:\n", "\n", "$$\\texttt{DataMat}\\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", "Why does this make sense? Think about it.\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 a $(d+1)$ dimensional numpy array called params." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "D = DataMat\n", "\n", "# Least Squares\n", "params = np.linalg.solve(np.dot(np.transpose(D), D), np.dot(np.transpose(D),y_data))" ] }, { "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. \n", "\n", "Once you're done, re-run all the blocks of code for different degrees $d$! What do you observe?" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "19.275906555946552\n" ] }, { "data": { "text/plain": [ "Text(0.5, 1.0, 'Polynomial of Degree 4, cost=19.275907')" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAegAAADUCAYAAABXn7BuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABBQElEQVR4nO2dd5xT1fLAv0NzlUXpiICCilKkydopVkBEbKCo+EAUxV55iv6eos9neaI8u2JDUBRBUVQQEFgFwQKIiChFRQXpIL0s7Pn9MTcQQrLJ7ia5SXa+n08+Se499565dc6ZM2dGnHMYhmEYhpFalPJbAMMwDMMw9sUUtGEYhmGkIKagDcMwDCMFMQVtGIZhGCmIKWjDMAzDSEFMQRuGYRhGCmIK2kdEpKeIuKDPRhH5XkRuFJEyhdxXrojkJkjUpCAidb3z0LMI28b1+EXkXBH5QUS2eTJVjFCuf5hruFBEholI+3jJkw6IyOEissU7D0f6LU80RORWEbkwjvvrISLvicjv3jkYXEDZ60TkZxHZLiJ/iMi/RaRsDHV0Capjq4jMF5FHRKRCSLnBIfdl8OfnkLKRyjUPKVdVRF4TkVVe3V9HusdFpHfQ8c0XkT4h6+sWUK8TkW7RzkVJoFBKwEgYXYElwIHe72eA6sB9fgrlA8uAk4Bf/BTCaxy9BUwDbgB2ABujbNYK2AUcANQDugCfisibQA/nXH7iJE4ZngfWA/v7LUiM3ApMBd6P0/66A9WACehzHBYR6Qf8BxgIfAo0Bx4AagJXR6njTuAP4B70ndEC6A+cJiInB91n/wZeDNm2LvA2MDrMfgcDL4UsWxAk837AJKAq8E9gOXAV8LGInOWcyw0q29vb1yPAZ8AZwPMiIs65F7xigWc9lIfQZ2lcuIMvcTjn7OPTB+gJOODIkOWTgfWF3FcukOv3Mfl4LuN2/MBh3nXpFUPZ/l7ZMmHW3eatuyPJ56IsIEmu8zJgBar09rmnU/EDLAbejOP+SgX9XgIMDlMmC23sDQ5ZfieQDzSOUke1MMv+4Z3z06Ns+y+vXOOQ5Q54KMq23b1ypwYtE2AO8E3QsjLASuCNkO1fA1YDZQuo4wBgAzDC73sjVT5m4k5NvgUOFJHqACLSQUSme2al9SLygYgcHWljETlYRHaIyC1h1vX3zJCVvP+5IjJVRM4UkVneurkickGYbaPKEbS/DiIy2yv7nYicICJlRORhEVkmIms9M1z5oG33MXGLyHEiMlJElgSZ9B4WkSL10kSkpogMEZHVnvltjoh0Dz4/6Isb4FVPntyi1OWcGwh8B+x1HUSkmoi8KCJLPRl+FpFrwsh6pnfutonIIhG52jtni4PKBM7Z9SLyXxH5C9gOVPTWXygiX3nX9W8RGSEih4ap6xrR4ZVt3rl5VUQqx3Kc3r30JKpk/o71/BSwv2YiMkpE1gRd835B60VEbvOW7/Dup2dF5MCQ/dwiIj95+1gnIjMC97V3Dg8DLpc9ZtXBxZHbxWYlOQbIBsaGLP8UVXjnR6ljVZjF33rftaLU/Q9gpnPux+hi7sOJwFYX1FN2qlXHA8eJSKDuk1Arwpsh2w8FqqC940hcCFQA3iiCfBmJKejUpB5qLt0kIh2AT4BNwCXAdehDPjXoodgL59xy4ANgr5e+iJRGzVLvOufWBa06AngKfcleiJqfRkjQOGIh5TgSeBx4FDX17Yea1V5AzXg9gQeBy4H7o5yLQ4HZQB+ggydnL+D1KNvtg9cY+Bw4GzURng/8AAwNUpCvsMc8+RD6wrm+sHUFMRaoE1CKnhKZCnREe9/nAB8BL4jITUGyNmLP+e7myXsLcHqEeu4FjkKv+QXANtFxv/eAeajJ/Vr0mn0uQWOWIvIo8BxqjuwM9EXP9VjvnonGf4GfnXNDYyhbICJyPDAdvSdvQ8/Pk0DtoGL/8ZZNAM716u8JfCIipbz9XA48gZp0O6L32kgg0Oi4ADXTjkOv8UmoWTjQACgTwyeWcxPKLu97R8jy7d73MUXYZ1vv+6dIBUTkFPS5jKT8rvMai1tEZJKItA5ZvwvIC7NdqNyNve+5IeUCjYJGkWQEeqC9708LKFOy8LsLX5I/7DFxH42ahiqhL9FdwAdemRnAQoJMqKgCzwOeDFqWS5CJFzjV23froGWdvWUnhmyXB9QPWlbdk+GeoGWFkSMPODxMvZ+FHP/7wG9B/+t65XpGOF/inafuqDmwSqTjj7D9jYSY6bzln6EvhtLe/yMLkiNk2/5EMHF766/11p/g/f8XsC34fHvLX0ZNgGW8/8OAVcABQWVqetsuDnPOZhFk1kZ7aeuB10LqqYcqh1uDtt8F3BdS7hRvv+dHOf7W6Eu6Ucg9XSQTN/AF8GfwcYesr+zVNzhkecAE29n7/ywwK0pdiwlj4mbPsxPtE/F+I7KJO9s734+FLA+YqccV8nzV8u7dCVHKveRd96ph1g1FG92tvfP4PfoMnxpU5npPvoYh207yll/q/b/H+58VUq6Mt/xfBRzHLoLeJfYxE3eq8DP6QKxFHW3eAnp5Pb5jgeHOuZ2Bws6534Av2dNy3genpqh5qIIIcC0wxzn3VUjxhc65hUHbrkQf+kCvr7ByLHDO/RpyfLCv48fPQG0RkUjHISIHishjIvIL+mLOQ18oAtSPtF0E2gBLXZCZzuNN1CxXUOu+qASOLZCVpgPwNfBbcG8MPTdVgmQ4ERjjnNsS2JFzbhnquBaOD5z3pvM4CXU6fCuknj/R897GK3cWakkLLfc1OlbahgiISDn0xT/QOTcv6pmIgogcgDYM3go+7hBOBMqxrwn1HWAne+7Fb4HmIvKMN1RwQCFEmQkcF8Pn2kg7iIRzbhM6HnujiHQTkYoichrwMKqgYnYmFJFs4EP0uK8soFwWcDHwsXNudRiZrnDODXfOTXHOvYmaof9CLUgBhqENyDdEpImoR/c97Lk/iusEeQV6Hw4u5n4yCvPiTg0uQFvcG4HfnXPbAESkNvqCXxZmm+XoGFpBvAAMEB2LzkaVw41hyq0Ns2w76tAC2rMvjBzrQv7vKGB5GaA0+pIJx+vAmahH+2xgM3A8apLNirBNJCoT+RgC6+NNHe87UG91tIcezlwIqqRBe8srw6xfARweZnnocVX3vj+LUM+6kHKLosgTjlvRe+Np2TMNLaAIK4hIBedcNO/3YCqhL+klBZQJXKO9jtc5t1NE1gStH4LeH1ehvb88ERkD3O6cWxxFjk3ovRYNF71IWO5Az+sw9Lnaht7f/yT8/bkPoj4YH6H3QlvnXEHnrDPqkxDT2K5zbqOIfIKeu8Cyv0WnpL2BOoaBzrbojw4NBOQO3FeVQo4lcF3CvWtALQiznXNzIqwvkZiCTg3mOufCvSDXoS+Bg8OsO5jIN3uAIehUh57oA7MF7Z0XluLKUSS8lv95QH/n3FNBy5sUcZdr0eGEUA4OWh9vOgJ/OOf+9P6vQRXvPg58HvO972XsUZ7B1IiwXaiyWON992TP+F8wG0PKtWPfBlTw+nA0Qs/d0jDrZqGm0uYFbB/KOrQnVpCzU+AaHUzQcXm9/iqB9Z414SXgJc+JrR06Jj0cOCGKHG3RmRTR+Bw1hxcK59wG4EIRqYYex2K0YfNf1D+hQETnS48EcoCznHM/RNmkB9r7HVNYUUPkniIiR6ANzNLoNKy+wFbU6gB7rklj9lbQAcvQPpYWETkOaIj6HBhBmIJOYZxzm0VkJtBVRPo753YBiMhhwMnofOmCtt8gIm+hprhs4G3v5ZBUOYrBfuiLILS32bOI+/scPYZTnHNfBi2/DFWaxTbTBiMit6EK6tagxZ8CN6FKO1wPOcBXQEcROSBg7hWRmqgJOJZe1jRUCR/pnCuo5zQBVYqHOucmxLDfYB5lX5NkB+AudCxzfugGBeGc2yIiU4HuIvKgc25rmGJfoZaXbsDEoOWXoO+z3DD7XQcMF5ET2NssvZ3wc7YDJu5oFMY6sA9OPbJXAYjIvagSHVHQNp4T3Fuos2CnMMNVoeVrAO2B55xzkaw2odscCHQCvgkjs0N9UQIm9t7AUOfcZq/IdO84Lmdv6013tPEU/NwF6IFa0IbFIl9JwhR06vMv1Jv3YxF5HlW0D6AOQE/EsP3z7HkphQYuSKYchcY5t15EvgLuEJFl6IPfi+jTSSIxGO25vu+9EJegL5KzgGsDDY8icoKI7ELNqoejXtNnoybBp4PKDUSVyRQRGYgqsfJAA9Sh7zyv3EPePsaJyAC0sfIv1MQddbzPa5z1BZ7zempj0WtVC+0h5jrnhjnnfhGRx4BnRafMfY6aXOug5+UV51zY3qRz7mf2+BcAOu3L+/l1sFVIdOrc68BpYXwAgrnTk2G6iDyBXqPDgebOuZucc2u95f1EZDPaK2zona+p6D2KiAxCFeh0tPF1FDrOOT6ornlAaxHphA5zrHbOLfbM8jMKkDEsnud9oKe4P3CYiHTx/n/uKWRE5BLU5DsftWxdiI4RXxQ8JCAi/0DHq89wzn3uLX4OnWXwH2CziJwYJMKSMKbuy9FGbthGmojciVqVJqPjzoeh1+Bgb9vgso+gjZfVaC+6L9p43j0FzjmXJyL/QgOTLEWV9Onoc3uTc25HyD7LoY2tsVEarCUTv73USvKHGD1e0V7JdNSUtB51DDk6pEwuEbxK0RfBtxHW5QJTwyxfzL6esrHKMTVkWV3vOK8OWd6fIA9ownhxe8vGoi/blah37jnsGzQh4vGH1FkTdTJbjfag5gDdQ8oUxYs78NmMjucOA9pH2KYSqqh/Q3uDK4EpeJ7VQeXOQsdCtwO/og2tUcB30c5t0PqO6Mt3AzrEsRB96TcKKXcF2jvdjI7B/uSd69rxuKfRiGz7eAFH2EcLdHz1b+9e+xm4K2i9oObQ+d75W4YqrgODyvTw7omV3vn7zTvnwWUaeOd9iyfb4MIcawz3QvAn+F69GJ3et8W7LuOBUwo4l8HbLi6gjv5h9vE98EMBMp+L9mpXo8p2DTol8vgwZV9DG0w7vO9ngMoR9nstagLf7t1z10cod4En+0XFOfeZ+hHvJBkZitcr+gno7Zx71W95jKLjmRQXAZ84566KVj6VEJFhQEXnXEe/ZTGMdMFM3BmK5wF+JGqGXoaN76QdIvIMOpb8F3AIap6vhAZrSTfaoD1HwzBixBR05nI1OnVjAXCZC+9wY6Q2WcBjqOf2DtRp50yXhlNRnHO1o5cyDCMYM3EbhmEYRgpikcQMwzAMIwUxBW0YhmEYKUixx6BFpA4asaoG6i4/yDn3lGiquuHoNJDFwMVu7wxK+1C1alVXt27d4ooE69fDokVw1FFQoUL08oZhGIaRCGLQRzNnzlztnKsWurzYY9BedKOazrlZoinsZqJp/HoCa51zj4rI3UAl59xdBe0rJyfHzZhR6PgA+7JmDVStCg8/DP36RS9vGIZhGIngvvvgP/+BDRugfPmwRURkpnMuJ3R5sU3czrllzrlZ3u+N6JzbWmgM5UD0mjeIkog8rlSpAg0awNSoYW0NwzAMI3F89RU0bRpRORdEXMegvTB/LdBUdTWcpscDDaMXKch/YmjTBr78EnYVJ3qjYRiGYRSR/Hz4+ms48cToZcMQNwXtRTl6Dw1XuFdCBqd29LC2dBG5RkRmiMiMVatWxUscaN1abf9z58Zvn4ZhGIYRKz/9pKZtPxW0l/7sPTTR+vve4hXe+HRgnDpsIHTn3CDnXI5zLqdatX3GyItO69b6/cUX8dunYRiGYcTKV16yMb8UtIgI8Crwk3PuyaBVo9GA9XjfHxa3rkJx2GFw6KEwZUpSqzUMwzAMQIdZK1eG+vWLtHk8Qn2egmbC+UFEZnvL7kFzxb4rIlcBv+NHHN7WreGzz8A5EEl69YZhGEYJZsoUaNUKShWtL1xsBe2cm4qmfwvHGcXdf7Fo0wbeekvnoBWxBWMYhmEYhWb5ctU911xT5F1kdiQxG4c2DMMw/ODLL/U7oIeKQGYr6AYNNGCJjUMbhmEYyWTqVMjKgmOPLfIuMltBi2jrxXrQhmEYRjKZMgVOOAHKlSvyLjJbQYMq6N9+gyVL/JbEMAzDKAls3AjffVcs8zaUBAXdpo1+m5nbKAFs2gTTp+t3JtVlGGnF119Dfj7zKrcq1vOR+Qq6WTPNIGJmbiPD2bQJGjWCdu30O5GKM5l1GUa6sWPiFHZRijP/dVKxno/MV9Blyug8tNxcvyUxjITyww+wbp2+DNat0/+ZUJdhpBtbJ0zlh1LNWLb5wGI9H/EIVJL6nH469O0LS5dCrVp+S2MYCaFJE6hUSX9XqqT/CyQ/X+dq/vGHatn16+HvvyEvTwMriGgDt3JlzRBXtao+P5UrF74uwygp5OVx4E9f8d0BV5FN8Z6PkqGgz/DipUyeDN27+yuLYSSI7GyYN09b602a6P/dLF8OM2fCrFn6+flnWLwYtm0rfEWVK5Ndvz6/tDqKJdVbUqPzCRxQpjmQFZ8DMYx0ZvZsZMsWLn2jFQ3qh3kWC0HJUNDNmmkzZtIkU9BGRpOdDSedhPaGR03SULcTJ8L8+VpABI46Cho3hnPOgXr1NG59lSpQsSIcdBDst5/2rvPztTe9bh2sXq2fP/6AhQth4ULK5n5GvWVD4SmgbFk47jho3x46dICWLaF0aR/PhGH4hOeQnHVmK046pHi7KhkKulQpOO00fVF5cbk3bYrQ0zCMdGXdOvjwQxgxAiZMUOVavjy0bQu9e8Pxx0Pz5uo0WRhq1468bulS+OYb9VqdNAn694f771dz+IUXQrduOpOiAGVtz6KRUeTmwhFHwCHF1M6UFAUNOg79/vvw229sqn44jRrp+6xSJTUL2ovBSEvy87WX/PLLqpzz8rRHfPPN0LmzprkLCZQQi0KMWWnWqgUXXKAfgFWrtHHw8ccaB3/QIHZUPQQuv5xyN12rL66QeuxZNDKGXbt0xlDXrnHZXeZ7cQc4/XT9njTJPFCN9Ofvv+HRR1XhtW+v/hU33qi92d9+gwEDtOcaRjlHmx5VrClU1arBZZfBsGFs+mUFN1R5hwnrcij11JNw5JFq/h49WhsWmDe4kWF8950OLwX0TTEpOQq6QQOoWRMmTdrtgZqdbR6oRpqxZAnceSfUqQP9+ukY8jvvqKn5ySd1HLiA1KqxKMR4Kc0ffi3PkO2X0GnXhzQ64Hf+vKq/7uy883QM/I03aNIgz55FI3OYPFm/Tz01LrsrOQpaRFs1kyaRXd4xbx6MH28mNSNN+PNPHUeuVw/+9z9VcrNn67jvJZeoY1cMxNI4jVcDNng/26rUotL/7lfP8bff1p59z55kNz+SBXe8yIQxefYsGunP5Ml7OoNxoOQoaFAFvWIFzJu329vVXghGSrNmjfaY69eHIUOgTx/NMfvmmzo7oZAEpmIV1DiNpUyR6ypbVh3HZs/Wcepatci69TpOvPoYsse9p06chpGO5OWpB/dpp8VtlyVPQYP2Ogwjldm+XceYDz8cBg7Ucd0FC+CZZ6Bu3WLtOpbGabwasBH3I6LTvL78Useky5SBLl20cCCPbhGw+OCGb8ycqTeeKegiUreumgg/+8xvSQwjMuPGqX24Xz8dy5ozB157Tb2zMw0ROPdcPcZXX9Ux9latoFcv9Qj3iEXxWnxwI1mEvR/jPP4MJU1Bgz69kyerOcIwUok//tC5wx066P9PP9WpU40b+ytXMihdWpXy/Plw990wdCgcfTQMGsSmDfkxKV7zCDeSQcSG4OTJcMwxOpMhTpQ8Bd2+vebqnD7db0kMQ3EOBg3Sh3vcOHj4YdUu7dv7LVnyKV8eHnkEvv8emjaFa68lv3UbKq1ZFFXx2uwMIxmEbQju2KFDM3E0b0NJVNCnn67jXZ9+6rckhgG//65N8Wuv1SlSP/6opu0wXtklany1USPtkQweTIU/fmT61mbcXu5ZKlfMj6h44+XcZhgFEbYh+M03sGWLKehic9BB6ogybpzfkhglGefgpZe01/zVV/Dii+obEcEBrESOr4pAjx7I3LmUO6stT+y4iV8OP5PsNb9H3MRmZxiJJmxDcPJkvV/bto1rXSVPQYOO8c2apVOuDCPZrF2rY819+mgozrlztQddzAAjGUutWpT59BN45RXKfDdD44mPGuW3VEYJZp+G4Gef6X1ZuXJc6ymZCjowtjdhgr9yGCWPKVP0Qf7kE3jiCbXkxOCdXeLHV0Xgqqt0/nT9+trAuemmoqXLNIx4snEjTJuWEJ+RkqmgW7RQTzsbhzaSxa5d8OCDOgWjXDl9oG+/XTOtxYCNr3ocfjhMnarn7tlntRuzYIHfUhklmdxc2LlTx5/iTMlU0KVKaWtn/PjdQfsNIxFs2gTfjlvLzvYdNQ3jpZfq8EpOTqH3ZeOrHuXKqfXho490alpOjk5HMww/GD8eDjgATj457rsumQoaVEGvWqXZRwwjhFgDYxRUZtMmuPDIOVQ9O4ddE3PZ9vQgnd974IGJEbqk0amTmrwbNIDzz1cLhTW4jWQzfrxaxmKMh18YSq6CDpgjzJvbCCFeKRn/evIdRq04iXJuO2fv/znf5fQu0BHMKAJ16mj+3X/8Qy0UF10EGzb4LZVRUli8WIdYEhSzoOQq6OrV4dhjYexYvyUxUoxip2TctQv69uWo+y/lx3LH0uaAmSyqemLJc+xKFllZMHiwZvn66CP1jP/ll2LtskTNOTeKzvjx+p2A8WcoyQoa1EQ2bZpmDDIMj2KlZNy0CS64AAYMgOuvp9Gyibz52cEl27ErGYjALbfozIwVK1RJFzFaYImcc24UjfHj1Ypz9NEJ2X1cFLSIvCYiK0VkbtCyyiIyQUQWet+V4lFXXDn3XB2zGjPGb0mMFKLIKRn/+gvatNEpVM89B889R3blcubYlUxOO00Dvxx0kEYNfO+9Qu+iRM85N2Jn506YOFFbcgkauopXD3ow0CFk2d3AROdcfWCi9z+1OPZYOOQQTXdnGEEUOiXj99/DCSfAwoVqZr3++qTJaoRQv772nlu0gK5d1eO7EHmmS/yccyM2ZsyAv/9OmHkb4qSgnXNfAGtDFp8HvOH9fgM4Px51xZVSpdTM/emnmn/XMIrC2LGaIhF0jm7Hjv7KY2icg4kT1WnszjvhxhvVNyAGbM65ERPjx2vP+YwzElZFIsegazjnlnm/lwM1whUSkWtEZIaIzFgVlP81aXTurLaszz9Pft1G+jN4sA6V1K8PX38NzZr5LZERYP/9YfhwVdDPP69z0GNsiNuccyMqY8dqgpsqVRJWRVKcxJxzDghrY3LODXLO5TjncqrFMY9mzJx+uj7IZuY2CsuAAXDllTru+fnnOlxSRMxrOEGUKgWPP65m7hEjtDFlJ9koLitXaoO8U6eEVpNIBb1CRGoCeN8rE1hX0dl/fzjrLB03LMQ4lVGCcQ7uugv69tUxzo8/hgoVirw78xpOArffDq++ChMnsuv0s/h23Fo7z0bRGTNG3wNprKBHAz283z2A1I3F17mzhgycM8dvSYxUZ+dOuPpq+O9/NRvV228XO4KQeQ0niV692DpkBDu/ncUBHdty2tF/mZI2isbHH6vFrHnzhFYTr2lWbwPTgaNFZImIXAU8CpwlIguBM73/qck55+j3Rx/5K4eR2mzbpj3m116D++7Tcc3SpYu9W/MaTh6zD7+Qi7LGcFj+b7y7rBXzx0fOLW0YYdmxQyNQduqU8MiA4lLIrJuTk+NmzJjhT+Unnqhent9+60/9RmqzZYsGIBk/Hp5+WlMdxpFNm7Tn3KSJOSYlksBwwhGrv+aDbR2oUOcgSuVOhnr1/BbNSBc++2zPsGicTNwiMtM5t08GnZIdSSyYCy7QeW2/W4vaCGHzZn0QJ0zQ3nOclTOY13CyCEyhenjiCZT9/DNKbdwAbdvCokV+i2akCx9/rOFlTz894VWZgg7QpYt+FyHykJHBbNgAHTqol/aQIeq1baQ1gcbQAa1bwqRJah059VTLK21ExzntOZ9xhqaYTDCmoAMccYRGHhoxwm9JjFQhECVo+nR1Buve3W+JjHjTvDlMnqzjiqeeCj//7LdERiozfz78+mvCvbcDmIIOpksXjeP7559+S2L4zdq1cOaZMGsWjBwJF1/st0RGomjSBHJzNS6/KWmjIAKOxAHH4gRjCjqYgJn7/ff9lcPwl9WrNfjI3LkwahScf77fEhmJplEjVdKgY4s2Jm2E46OPNFpgnTpJqc4UdDBHHQVNm2qPySiZrF+vZu0FCzS6XJJaykYK0KCBeuju2KFjjOYwagSzfLnG2r/ggqRVaQo6lC5d4MsvNXWgUbLYvFkV8ty5akVJYJYaI0U55hj11t+wQXvSS5f6LZGRKnzwgTqJXXRR0qo0BR1Kly56EczMXbLYvl1bxtOnw1tvwdln+y2R4RctWmggilWrtCe9YoXfEhmpwHvvaVKcxo2TVqUp6FAaNtQLYN7cJYedO6FbN+05vfKKRgszSjbHH6/xlv/8U50FV6/2WyLDT9auVW//iy5KePSwYExBh+Pii2HKFPPmLgnk5+vc5g8+gKeesnnOxh5atVI/hIULdbjj77/9lsjwi9GjNdJkEs3bYAo6PJdfDs6xffDblgIwk3EObrwR3nwTHnoIbr7Zb4mMVOOMM9STf+5cnfu6ZYvfEhl+8N57cOih0LJlUqs1BR2OI45g1/En8du/h9LuLGcpADMR5+Duu+GFF+Cf/4R77vFbIiNVOfts9UuYNk2HP/Ly/JbISCYbNmgM/gsvTKp5G0xBR+T3NlfQIG8uh2+eYykAM5FHHtmTMvLRR5P+4BlpRteu8OKLOi7ds6cOjRglgzFjdOpdks3bYAo6ItVvvJg8ytCr7JuWAjDTeOYZuPdeDd353HOmnI3YuOYaePhhGDYMbrlFrTBG5vPee3DwwXDyyUmv2hR0BLIPq4J07Mh1Bw1j3g+7LMtQmrBpEwX7DQwerGPN558Pr78OpewRMArB3XfD7bfDs8/Cgw/6LY2RaDZu1OxVF17oy7vC3k4FUKZnd8qt/ovsGbl+i2LEQCDXb7t2hPcbGDkSrrpKc7m+8w6UKeOLnEYaIwIDBqiZu39/tcYYmcsHH8C2bXDZZb5Ubwq6IM49Fw48EIYO9VsSIwZ++AHWrVPFvI/fwNix+pCdeKJ65e63n29yGmmOCLz8Mpx3nlpj3nrLb4mMRPH223DYYZqf1AdMQRdEVpY6h7z3nrlxpwFNmkClSprvdy+/gc8/VxPVMcfAJ59A+fK+ymlkAGXKqBWmbVvo0UPvKyOzWLVKvbcvvdS3oTBT0NHo1UuV87vv+i2JEYXsbJg3T5+pefP0P99+q5aQunU1fGPFij5LaWQMWVkawKJZM23If/WV3xIZ8WTECA1O4pN5G0BcCnki5uTkuBkzZvgtxt44p6E/DzpIvY+M9GHuXO3hHHSQRoarVctviYxMZMUK9fBdv14T7Rx9tN8SGfGgdWsdK5s7N+FVichM51xO6HLrQUdDBHr31tZxEi6UEScWLVJnsP320xSCppyNRFGjhlpnSpWC9u1h2TK/JTKKyx9/aGpJH3vPYAo6Nq64AsqV00QKRuqzZIkmOMjLU+V8+OF7rY46FcswCsuRR+o49OrVGnlswwa/JTKKQ8Dxr1s3X8UwBR0LVatqKsKhQ9Xl3khdVq5U5bxunfZqGjXaa3XUqViGESP7NPSOO06n8v34ozolbt/uq3xGEXFOYyS0abNP4z7ZmIKOlauv1pRjo0b5LYkRiXXrVPP+8Yf2ZsIEti9wKpZhxEjEhl6HDmppmzjRQoKmK9OmaQazFMhsZwo6Vk4/XVtTL7zgtyRGODZtgnPOgZ9+0uACrVqFLRZxKpZhFIICG3o9emis93fegb59fZPRKCKvvaZTMbt08VsSU9AxU6oU3HCDegPPnu23NEYw27Zp6M5vvtGXYrt2EYuGnYplGIUkakPvrrs0lemTT8ITT/gio1EENm/WKbUXX5wSLwdT0IWhVy844AAL75dK5OXBJZeoSfH119VXIArZ2RoYKAWePyNNidrQE4H//U97YXfeqQk2jNRn5Eg1i/Tq5bckgCnowlGxonp0DxsGa9b4LY2xa5eO840erVmprrjCb4mMEkTUhl7p0upY2rat3qcTJyZTPKMovP461K8Pp5zitySAKejCc+ONalK1KVdJZy+vWefg+uu1sfToo/rbMFKNrCz1iTj6aLXufPed3xIZkViwQMMC9+yZMiloTUEXlmOOUYex55+HnTv9liZjiDY3eS+v2YaOHbf0hUGD4J57dLzPMFKVihU1WUvFijpH+rff/JbICMeLL0KZMmy+pFfKxElIuIIWkQ4iMl9EFonI3YmuLyncfLNO5Xn/fb8lyQhimZsc7DXbe+VDlHvmCbVmPPRQ8gU2jMJSuzZ8+ins2KHRxlav9lsiI5gtW2DwYPI6X0jD0w5OmTgJCVXQIlIaeA44G2gEXCoijQreKg0491w1WT32mJpajWIRy9zkgNfsP8v9j3/tuI+8y3rAU0+ljCnKMKLSqJH6S/z5p04J3LzZb4mMAMOHw7p1LDjj+pSKk5DoHvTxwCLn3K/OuR3AO8B5Ca4z8ZQqpfMbZ83SUJJGsYhlbnJ2Niy461Ue23EbO8+7iLJvvOJbCjjDKDKtWulUwBkzNANWXp7fEpUoIg6lvfACNGrEYVe0Sak4CYl+w9UC/gz6v8RbthsRuUZEZojIjFWrViVYnDjSvTsccoj2oo1iEdPc5OHDybqpN3ToQJnhb2k+XsNIR847T8c7x47VCIVmhUsKEYfSZszQtLTXXUd2BUmpOAm+d0Gcc4OccznOuZxq1ar5LU7s7Lcf3HabTp1ItRSZaUiBU1ZGj9YGUatW8N57eu4NI53p3RsefBCGDIG7M8M1J9WJOJT23HMa38KbpplKcRISraCXAnWC/tf2lmUG11yjuYYffdRvSTKX8ePVFHjssfDxx/ogGUYm8H//p9MD//tfGDjQb2kynrBDacuX61TNnj31XZ5iJFpBfwvUF5F6IlIO6AaMTnCdyePAA+Gmm7RXN2eO39JkHl98oSE8GzRQc+CBB/otkWHEDxF4+mm46CK4/fY9KQ6NhBB2KO3ZZ9UP4Lbb/BYvLAlV0M65ncCNwDjgJ+Bd59yPiawz6dx+u7a87r/fb0kyi2++UU/Xww6DCROgcmW/JTKM+FO6NLz55p5oY+PH+y1RRrOX+XrzZnUOO+88zeedgiR8DNo5N8Y5d5Rz7gjn3H8SXV/SqVRJlfQHH8DMmX5LkxnMnq1zRatXVy/56tX9lsgwEkdWFnz4ITRurHmkzaclObzxhqYQvuMOvyWJiLgU8iDMyclxM9Lx5tywAerVgxNP1DzERtH56SdNlJ6VpZnD6tb1WyLDSA7LlsHJJ2vP7ssvNSa0kRh27dJYFlWqwFdf+R5PQURmOudyQpf77sWdERx4oM6LHjNGJ9kZReOXX+CMM9TsN3GiKWejZFGzJowbp9Ou2rdXByYjMYwape+bO+7wXTkXhCnoeHHTTVCjhirqFLJKpA2//AKnngrbt6tZ+6ij/JbIMIpNtBjz+3DUUdrQX7lS43Zv2JC4ukoq+fnw73/rub7oIr+lKRBT0PGifHmNC/3ll5pT1IidRYtUOW/dCpMmaUISw0hzYokxHyi3l2I97jidGTJ3rs5i2LYtbnUZwEcf6aybe+9Va10KYwq6mOz1cF15JTRtqtmVYnioDFQ5n3aaKueJE6FZM78lMoy4EEuM+YiKtX17GDwYJk+GSy6JGhI0lroM1Lr54INwxBFw2WV+SxMVU9DFYJ+Ha2tpeOIJTSf39NN+i5f6BPecTTkbGUYsMeYLVKyXX65RrkaP1ilYu3YVqy4DdeKdNUt7z2kQLtgUdDEI+3CdeSZ06gT/+Y96ZaYpCR/PCijn7dvVrG3K2cgwYokxH1WxXn+9RiocNkx/R/BviSmefUkn0HuuW1dDB6cBqd+ESGECDxeEPFxPPql/brtNM9ekGQHLwLp1elxxf+Dnz1dv7YBytua+kaEEAmMUtH7ePG3cN2kS4Tm76y51Fnv4YahQAR5/PKzncbS6SjzvvadJMV59FcqW9VuamLAedDGI2GqtX19NKMOHa4jKNCOh41mzZ0Pr1jqmZsrZMGJLzvDQQzpT5Ikn1APZKBx5eXDPPRoMpkcPv6WJGetBF5OIrdZ//hPeflvNUnPnqpd3mhDRMlBcpk2Djh113rhNpTKM2BGB//0PNm7UsMIVKqRs/OiU5OWXYeFC9eBOcc/tYKwHnSj22w9eegkWL44pTncqzWGM93jWpk0w76kJuLPO0rCdU6eacjaMwlKqlCqaLl00vPCgQX5LlB5s3AgPPKARCs85x29pCoX1oBNJ69Zw7bU6Jn3uuRoQPwwJH/MtAvEaz9q0CW6rO4pn13RjXtmG1P10HOUPrVH8HRtGSaRMGc16tWWLvltKlYKrr/ZbqtTm8cc18Mvo0SkdNSwc1oNONE88oZlS/vEP+PvvsEUyeQ7jygdf5MU1XZhJSzqUm8ycFaacDaNYlCunDk9nnw29e8Mrr/gtUeqyaJHm2770UjjhBL+lKTSmoBNN+fKaTm7pUrjxxrBFMnIOY34+9OvH4Y9fx+SsjlxQfgJSuVJmHJth+E1WFrz/vinpgnAObr5ZGzQDBvgtTZEwBZ0Mjj9ex6HfeguGDNlndcbNYdy+Ha64Qudv9unDictG8cGE8plxbIaRKgSUdMeOGa+ki+Sj8+GHOovmgQfgkEMSJlsisXSTyWLnTg05Nn26xus+9li/JUoM69ZpTtvcXHjkEZ3DmWbjPoaRVmzbpkkfxoxRx7Hevf2WKK4UyUdn82bd6MADNXJYis97tnSTflOmjM6LrlZNFdjq1X5LFH/mzVNrwZdfqrXg7rtNORtGognuSV9zDQwc6LdEcaVIPjr9+sEff8Dzz6e8ci4IU9DJpFo1fZCWL9cA+Dt2+C1R/PjwQ3XC2LhRA/ynQSB6w8gY9ttPcxx37apTsO67L2PS3hbaR2fyZHjmGR1/bt06KTImClPQySYnR81QkyZBr17qTJXO5OdrfNvzz4cGDWDGDDjlFL+lMoySR7lyGhzpqqs02tjNN6f/+4VC+uhs3Kjv1SOP1CG2NMfmQfvBP/4BS5ZoONAaNXQqVjqyapWGzRs7Vp3CXnoJ9t/fb6kMo+RSurQGM6lUST2X//4bXnstrc28UIi4DHfeCb//DlOmwAEHJFyuRGMK2i/69VNT95NPQuXKqqyTxKZNUYLzx8Lnn6sZe80aTYl33XU23mwYCaJQz6yIzv2tVEnfK+vWadKeTJ9C8e67ap3s2zdjrHimoJNA2IcrEFt37Vr4v//T8ej+/ROu5IodtSwvTwP3P/SQmpE++QSaN0+UuIZR4inSMyuiySGqVNF8AG3aaBzqWrWSInPSWbhQI6qddJKm+s0QbAw6wQQernbt9HuveXylSsEbb+iYyYMPwh13JHzMqFhRy2bPVi/tBx/UZPIzZ5pyNowEU6xn9tpr4eOPVYGdcAJ8/33C5AxHUnIMbNkCF1+sZvx33kl7c34wpqATTNSHKzBmdPPNOj3ioot0Dl+CKFLUsu3btXd/3HGwbJl6iw4ZkvkmM8NIAYodafDsszVBjQi0aqXzpZNAgZ2TeJGfr34w338PQ4fCoYcmoBL/MAWdYGJ6uEqVUnP3009rQPdTToEFCxIiT6Gjlo0Zo0I/8IBODZs3Tz22DcNICnGJNNisGXz9teaqP/dcjfIXYRpWvHq9SckxcP/9MHIkDBjApjYdUyYjYNxwzqXMp2XLli4T2bjRuWnT9DsqY8Y4V7myc+XLO/f6687l5ydavPD8+KNz55zjHDh31FHOjR3rjxyGYcSPjRudu+QSfa4vuMC5v//eZ3WdOs5lZ+t3TO+sAqqK177CMmiQHkevXm7jhvzE1pVggBkujE60HnQSCEwRiKnle/bZaq7JyYErr9T/ixYlXMbdzJ+v48vHHANffKGp2n74ATp0SJ4MhmEkhuxsnSs9cKBa644/Hn78cffqePZ6E5pjYPhwHV8/+2x44QV+mCsZmRHQFHQqUrs2TJyoJu9p01RZ3nGH5jRNBM7ptKmuXXWw6IMP4J//hF9/1XmF5colpl7DMJKPCNx6qwZLWr9elfRLL4Fzcc+sV6jOSayMGgXdu+t4+siRUK5cZmYExJJlpD7Llul0iSFDNOZuz57acmzatPj7/vNPvcFff12bnJUra6D922+H6tWLv3/DMFKbv/5SJ6vPPoNOneCVV9hUvkbx4yQkiiFDdNZLTg6MGwcHHbR7VVziO/hEpGQZpqDThQULdH7f8OHqVd28OXTurKbnFi1UeUdj82Z1FMnNhQkT4KuvdHnLljpX8tJLLRKYYZQ08vM1dvVdd2n2pxde0IQ+qRR4yDkN6nTnnXDGGWrlSzctXAAJUdAi0hXoDzQEjnfOzQha1w+4CtgF3OycGxdtf6agY2DNGp1OMHKkulrm52umrMaN4bDD1Dx+wAFqlt68WQdkVqyAn37S7C6gXuMtW6o3dteu6tlpGEbJZu5cDdk7e7Z6ej/zjL5T/GbbNujTR2NGXHQRvPlmbB2SNCJRCrohkA+8BNwZUNAi0gh4GzgeOAT4DDjKOberoP2Zgi4kq1erI9eMGTBnjirgpUth61btZZcvrwMyVatqIouGDTUPdatW2lI2DMMIJi8PnnpKpy+BRjm85Zawca2TYlL++Wd1Wp01S2Mx/Otf2sHIMBJq4haRXPZW0P0AnHOPeP/HAf2dc9ML2o8paMMwjBTg999VMX/4oYYHffBBHasuXRqIQ8jgaOzapbmc77pLGwevvgrnnRfHClKLSAo6UU2RWsCfQf+XeMvCCXaNiMwQkRmrVq1KkDiGYRhGzBx2mI7zfvGFDptddZVa4F58EbZuTWwQkilT1Ans5puhbVvdeQYr54KIqqBF5DMRmRvmE5cz5pwb5JzLcc7lVKtWLR67NAzDMOJB69bq6/L++1CxomatO+wwWo76P1pkL4zftCbn1Hn1rLM0sceaNZqdaswYqFkzDgeSnkTNZuWcO7MI+10K1An6X9tbZhiGYaQTInDBBepU+sUXMGAA5Z54hC/y/8OGJiez38Xns9/is9VRtTCe385pEKb334fBg3W8uUYNeOwxuOEG9aEp4SRqDLoxMIw9TmITgfrmJGYYhpEBLF0Kb72lnzlzdNnBB+vskBYtdGZIzZqqcAPZpdat02BLCxboNl9+CYsX67pWrTTGw2WXlcipnony4r4AeAaoBvwNzHbOtffW3Qv0AnYCtzrnxkbbnyno9COdgwMYhhEHlizRoCG5ufDddzqlM1ra3Nq1dZy5XTuN5VCvXlJETVUsUIkRdxLuyZmG5OXlsWTJErZt2+a3KCWCrKwsateuTdkMygGc9mzdqkp72TLtMe/apebsihWhWjVVxpUr+y1lShFJQUcdgzZSh1TrrQZ7cgb+n3SSvzL5zZIlS6hQoQJ169ZFUikSUwbinGPNmjUsWbKEeiW8B5ZS7L+/mrgtAFKxybwZ3xlKUpKfF5JMDVBfHLZt20aVKlVMOScBEaFKlSpmrSgE8cr1bCQH60GnCanYWw2kk0ulXn0qYMo5edi5jh0bkko/rAedJqRqbzUh6eSMYlG6dGmaN29O48aNadasGU888QT5UZx2Fi9ezLBhwxIu29VXX828efMKLPPBBx9ELWMUnoQGFzESginoNCGhyc+NjGL//fdn9uzZ/Pjjj0yYMIGxY8fywAMPFLhNshT0K6+8QqNGjQosYwo6MaRqI9+IjCnoNMJ6q5lJIscFq1evzqBBg3j22WdxzrF48WJat27Nsccey7HHHsu0adMAuPvuu5kyZQrNmzdn4MCBEcsFs3jxYho0aMDll19Ow4YN6dKlC1u2bAFg4sSJtGjRgiZNmtCrVy+2b98OwKmnnkpgpkZ2djb33nsvzZo148QTT2TFihVMmzaN0aNH07dvX5o3b84vv/zC008/TaNGjWjatCndunWL/0kqIVgjPw1xzqXMp2XLls4oHhs3Ojdtmn4byWfevHmFKr9xo3N16jiXna3f8bhu5cuX32fZQQcd5JYvX+42b97stm7d6pxzbsGCBS7wzE2ePNmdc845u8tHKhfMb7/95gA3depU55xzV155pXv88cfd1q1bXe3atd38+fOdc85dccUVbuDAgc4559q2beu+/fZb55xzgBs9erRzzrm+ffu6f//7384553r06OFGjBixu56aNWu6bdu2OeecW7du3T5yFPacG6mFvbOcA2a4MDrRetAZRCp6ehsFk+xxwby8PHr37k2TJk3o2rVrRFNyrOXq1KnDKaecAkD37t2ZOnUq8+fPp169ehx11FEA9OjRgy+++GKfbcuVK0enTp0AaNmyJYsDUaVCaNq0KZdffjlvvvkmZcqYX2smYe+sgjEFnUGYE0j6kYxxwV9//ZXSpUtTvXp1Bg4cSI0aNfj++++ZMWMGO3bsCLtNrOVCvagL41VdtmzZ3eVLly7Nzp07w5b75JNPuOGGG5g1axbHHXdcxHJG+mHvrIIxBZ1BmBNI+pHoccFVq1bRp08fbrzxRkSE9evXU7NmTUqVKsXQoUPZtUvD41eoUIGNGzfu3i5SuVD++OMPpk/XNO/Dhg2jVatWHH300SxevJhFixYBMHToUNq2bRuzzMGy5Ofn8+eff3Laaafx2GOPsX79ejZZNytjsHdWwZiCziDMCSQ9ibfz39atW3dPszrzzDNp164d999/PwDXX389b7zxBs2aNePnn3+mvJcxqGnTppQuXZpmzZoxcODAiOVCOfroo3nuuedo2LAh69at47rrriMrK4vXX3+drl270qRJE0qVKkWfPn1ilr9bt248/vjjtGjRgoULF9K9e3eaNGlCixYtuPnmm6lYsWKxz5GRGtg7q2AsFrdhxJGffvqJhg0b+i1GUli8eDGdOnVi7ty5vspRks65kZlEisVtPWjDMAzDSEFMQRuGUSTq1q3re+/ZMDIZU9CGYRiGkYKYgjYMwzCMFMQUtGEYhmGkIKagDcMwDCMFMQVtGEbcGDFiBA0bNuS0005jxowZ3HzzzQDk5uaGTbhhGEZkLLCtYRgA7Ny5s9ixrl999VVefvllWrVqBUBOjk7tzM3NJTs7m5NPPrnYchpGScF60IaRYQwZMoSmTZvSrFkzrrjiCgB69uzJyJEjd5fJ9kI25ebm0rp1azp37kyjRo24++67ee6553aX69+/PwMGDADg8ccf57jjjqNp06a7I5MF8+CDDzJ16lSuuuoq+vbtS25uLp06dWLx4sW8+OKLDBw4kObNmzNlypREHr5hZAzWgzaMRHHrrTB7dnz32bw5/O9/EVf/+OOPPPTQQ0ybNo2qVauydu3aqLucNWsWc+fOpV69enz33Xfceuut3HDDDQC8++67jBs3jvHjx7Nw4UK++eYbnHN07tyZL774gjZt2uzez3333cekSZMYMGAAOTk55ObmAjpfuk+fPmRnZ3PnnXcW5+gNo0RhCtowMohJkybRtWtXqlatCkDlypWjbnP88cdTr149AFq0aMHKlSv566+/WLVqFZUqVaJOnTo89dRTjB8/nhYtWgCwadMmFi5cuJeCNgwjvpiCNoxEUUBPN9mUKVOG/Px8QDNEBaePDE2E0bVrV0aOHMny5cu55JJLAHDO0a9fP6699trkCW0YJRwbgzaMDOL0009nxIgRrFmzBmC3ibtu3brMnDkTgNGjR5OXlxdxH5dccgnvvPMOI0eOpGvXrgC0b9+e1157bXeqx6VLl7Jy5cqY5QpNZ2kYRnRMQRtGBtG4cWPuvfde2rZtS7Nmzbj99tsB6N27N59//jnNmjVj+vTpEdNHBvaxceNGatWqRc2aNQFo164dl112GSeddBJNmjShS5cuhVK45557LqNGjTInMcMoBJZu0jDiiKU+TD52zo10x9JNGoZhGEYaYQraMAzDMFKQYiloEXlcRH4WkTkiMkpEKgat6ycii0Rkvoi0L7akhmEYhlGCKG4PegJwjHOuKbAA6AcgIo2AbkBjoAPwvIiULmZdhpEWpJJfR6Zj59rIZIqloJ1z451zO72/XwG1vd/nAe8457Y7534DFgHHF6cuw0gHsrKyWLNmjSmOJOCcY82aNWRlZfktimEkhHgGKukFDPd+10IVdoAl3jLDyGhq167NkiVLWLVqld+ilAiysrKoXbt29IKGkYZEVdAi8hlwcJhV9zrnPvTK3AvsBN4qrAAicg1wDcChhx5a2M0NI6UoW7bs7rCZhmEYxSGqgnbOnVnQehHpCXQCznB77HpLgTpBxWp7y8LtfxAwCHQedHSRDcMwDCPzKa4Xdwfgn0Bn59yWoFWjgW4isp+I1APqA98Upy7DMAzDKEkUdwz6WWA/YIKIAHzlnOvjnPtRRN4F5qGm7xucc7uKWZdhGIZhlBhSKtSniKwCfo/jLqsCq+O4Pz+xY0lNMuVYMuU4wI4lVcmUY0nEcRzmnKsWujClFHS8EZEZ4eKbpiN2LKlJphxLphwH2LGkKplyLMk8Dgv1aRiGYRgpiClowzAMw0hBMl1BD/JbgDhix5KaZMqxZMpxgB1LqpIpx5K048joMWjDMAzDSFcyvQdtGIZhGGlJ2itoEekqIj+KSL6I5ISsi5ryUkTqicjXXrnhIlIuOZIXjCfLbO+zWERmRyi3WER+8MrNSLKYMSEi/UVkadDxdIxQroN3rRaJyN3JljMWCkqxGlIuJa9LtHPsBRca7q3/WkTq+iBmVESkjohMFpF53vN/S5gyp4rI+qD77j4/ZI2FaPeLKE9712WOiBzrh5wFISJHB53r2SKyQURuDSmTstdERF4TkZUiMjdoWWURmSAiC73vShG27eGVWSgiPeImlHMurT9AQ+BoIBfICVreCPgeDaRSD/gFKB1m+3eBbt7vF4Hr/D6mMDI+AdwXYd1ioKrfMkaRvz9wZ5Qypb1rdDhQzrt2jfyWPYyc7YAy3u/HgMfS5brEco6B64EXvd/dgOF+yx3hWGoCx3q/K6DpbkOP5VTgY79ljfF4CrxfgI7AWECAE4Gv/ZY5yvGUBpaj83vT4poAbYBjgblBy/4L3O39vjvc8w5UBn71vit5vyvFQ6a070E7535yzs0PsypqykvR8GenAyO9RW8A5ydQ3ELjyXgx8LbfsiSY44FFzrlfnXM7gHfQa5hSuMgpVtOBWM7xeehzAPpcnOHdgymFc26Zc26W93sj8BOZnTHvPGCIU74CKopITb+FKoAzgF+cc/EMPJVQnHNfAGtDFgc/D5H0Q3tggnNurXNuHTAB6BAPmdJeQRdALeDPoP/hUl5WAf4OeuGmYlrM1sAK59zCCOsdMF5EZnqZwVKVGz3T3GsRzESxXK9UoxfaqwlHKl6XWM7x7jLec7EefU5SFs8M3wL4Oszqk0TkexEZKyKNkytZoYh2v6Tb89GNyJ2KdLkmADWcc8u838uBGmHKJOzaxDMfdMKQGFJepiMxHtelFNx7buWcWyoi1dGY6D97LcGkUtCxAC8A/0ZfQv9GTfa9kidd4Yjlukj0FKspcV0yHRHJBt4DbnXObQhZPQs1sW7y/B4+QBP3pCIZc794fjydgX5hVqfTNdkL55wTkaROe0oLBe2ipLyMQCwpL9egpqIyXm8hYlrMRBDtuESkDHAh0LKAfSz1vleKyCjUjJn0BzvWayQiLwMfh1kVc4rSRBPDdenJvilWQ/eREtclhFjOcaDMEu/+Owh9TlIOESmLKue3nHPvh64PVtjOuTEi8ryIVHXOpVw86Bjul5R5PmLgbGCWc25F6Ip0uiYeK0SkpnNumTeksDJMmaXo2HqA2qhPVLHJZBN31JSX3st1MtDFW9QDSKUe+ZnAz865JeFWikh5EakQ+I06MM0NV9ZPQsbKLiC8jN8C9UW96suhJrLRyZCvMEjkFKvBZVL1usRyjkejzwHoczEpUiPET7xx8VeBn5xzT0Yoc3Bg/FxEjkffdynX2IjxfhkN/MPz5j4RWB9kek01Ilr90uWaBBH8PETSD+OAdiJSyRu+a+ctKz5+e84V94O+8JcA24EVwLigdfeiXqvzgbODlo8BDvF+H44q7kXACGA/v48pSM7BQJ+QZYcAY4Jk/977/IiaYH2XO8xxDAV+AOagN3zN0GPx/ndEvXF/SeFjWYSON832PgGP57S4LuHOMfAg2uAAyPKeg0Xec3G43zJHOI5W6JDJnKBr0RHoE3hmgBu98/896tB3st9yRziWsPdLyLEI8Jx33X4gaMZKKn2A8qjCPShoWVpcE7RRsQzI83TKVaj/xURgIfAZUNkrmwO8ErRtL++ZWQRcGS+ZLJKYYRiGYaQgmWziNgzDMIy0xRS0YRiGYaQgpqANwzAMIwUxBW0YhmEYKYgpaMMwDMNIQUxBG4ZhGEYKYgraMAzDMFIQU9CGYRiGkYL8P/2aSDlXm+ehAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "cost=np.linalg.norm(y_data-np.dot(DataMat,params))\n", "print(cost)\n", "\n", "fig=plt.figure(figsize=(8,3))\n", "ax=fig.add_subplot(111,xlim=[-11,11],ylim=[-21,21])\n", "x,y=poly_curve(params,x_data)\n", "\n", "ax.plot(x_data,y_data,'.b',markersize=5)\n", "ax.plot(x,y,'r') \n", "ax.legend(['Data points','curve fit'])\n", "plt.title('Polynomial of Degree %d, cost=%f' %(len(params)-1,cost),fontsize=16)" ] }, { "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.7.4" } }, "nbformat": 4, "nbformat_minor": 1 }