{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy\n", "from matplotlib import pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear regression\n", "\n", "Data from Example 6.1 of Seborg, Edgar, Melichamp and Doyle (3rd edition)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import pandas" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "df = pandas.read_excel('../../assets/example_6_1.xlsx')" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "x = df['Fuel Flow Rate']\n", "y = df['Power Generated']" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAQ2klEQVR4nO3dfWydZ3nH8e+FY1Sno3JHDTROs8A0vJcWLZXFeNG6jdC5jKrNqv1BpaKOVcqEECswzMiQhvYHKpvRBtIkpqgtLaJLx4rJJoTqVjDWIZUyp27nlNZD46XEKcRVZxjg0TS99oePg22S2Oec55zHt/39SJHt26fn+Slyfj2+n+c8V2QmkqTyvKDuAJKk1ljgklQoC1ySCmWBS1KhLHBJKtS2bh7soosuyt27d3fzkJJUvCNHjjydmQOr17ta4Lt372ZycrKbh5Sk4kXEt8+07haKJBXKApekQlngklQoC1ySCmWBS1KhunoViiRtJYenZhmbmOH4/AI7+vsYHRli357Byp7fApekDjg8NcuB8WkWTp4CYHZ+gQPj0wCVlbhbKJLUAWMTM6fLe8nCyVOMTcxUdgwLXJI64Pj8QlPrrbDAJakDdvT3NbXeCgtckjpgdGSIvt6eFWt9vT2MjgxVdgxPYkpSByydqPQqFEkq0L49g5UW9mpuoUhSoSxwSSqUBS5JhVqzwCPi9og4ERFHV62/MyJmIuKxiPjrzkWUJJ3Jel6B3wFctXwhIn4HuBZ4VWb+GvCR6qNJks5lzQLPzAeAZ1Ytvx34cGb+pPGYEx3IJkk6h1YvI3wl8JsR8SHg/4D3ZuZ/nOmBEbEf2A+wa9euFg8nSZ3R6TsGdlKrJzG3ARcCrwFGgU9HRJzpgZl5MDOHM3N4YOBnhipLUm2W7hg4O79A8tM7Bh6emq072rq0WuDHgPFc9FXgeeCi6mJJUud1446BndRqgR8G3gAQEa8EXgg8XVUoSeqGbtwxsJPWcxnhIeBBYCgijkXETcDtwCsalxbeDdyYmdnZqJJUrW7cMbCT1jyJmZnXn+VbN1ScRZK6anRkaMXUHKj+joGd5M2sJG1Z3bhjYCdZ4JK2tE7fMbCTvBeKJBXKApekQlngklQoC1ySCmWBS1KhLHBJKpQFLkmFssAlqVAWuCQVygKXpEJZ4JJUKAtckgrlzawkbUglz6rsFgtc0oazNKty6T7dS7MqAUt8GbdQJG04pc+q7BYLXNKGU/qsym5Zz0zM2yPiRGP+5ervvTciMiKcSC+pMqXPquyW9bwCvwO4avViRFwCXAk8WXEmSVvc6MgQfb09K9ZKmlXZLWsWeGY+ADxzhm/9LfA+wGn0kiq1b88gt1x3GYP9fQQw2N/HLddd5gnMVVq6CiUirgFmM/PRiFjrsfuB/QC7du1q5XCStqCSZ1V2S9MnMSNiO/AB4C/W8/jMPJiZw5k5PDAw0OzhJEln0cpVKL8IvBx4NCK+BewEHo6Il1UZTJJ0bk1voWTmNPCSpa8bJT6cmU9XmEuStIb1XEZ4CHgQGIqIYxFxU+djSZLWsuYr8My8fo3v764sjSRp3XwnpiQVygKXpEJZ4JJUKAtckgplgUtSoSxwSSqUBS5JhbLAJalQFrgkFcoCl6RCWeCSVCgLXJIKZYFLUqFaGqkmqRyHp2YZm5jh+PwCO/r7GB0ZclTZJmGBS5vY4alZDoxPs3DyFACz8wscGJ8GsMQ3AbdQpE1sbGLmdHkvWTh5irGJmZoSqUoWuLSJHZ9faGpdZVnPSLXbI+JERBxdtjYWEU9ExH9GxGcjor+zMSW1Ykd/X1PrKst6XoHfAVy1au1+4NLMfBXwX8CBinNJqsDoyBB9vT0r1vp6exgdGaopkaq0ZoFn5gPAM6vW7svM5xpffgXY2YFsktq0b88gt1x3GYP9fQQw2N/HLddd5gnMTaKKq1D+CPjHs30zIvYD+wF27dpVweEkNWPfnkELe5Nq6yRmRHwAeA6462yPycyDmTmcmcMDAwPtHE6StEzLr8Aj4kbgamBvZmZ1kSRJ69FSgUfEVcCfAb+VmT+uNpIkaT3WcxnhIeBBYCgijkXETcDfAS8C7o+IRyLi7zucU5K0ypqvwDPz+jMs39aBLJKkJvhOTEkqlAUuSYWywCWpUBa4JBXKApekQlngklQoC1ySCuVINakmzqpUuyxwqQbOqlQV3EKRauCsSlXBApdq4KxKVcECl2rgrEpVwQKXauCsSlXBk5hSDZZOVHoVitphgUs1cVal2uUWiiQVygKXpEJZ4JJUqPXMxLw9Ik5ExNFlaz8fEfdHxNcbHy/sbExJ0mrreQV+B3DVqrX3A1/IzF8CvtD4WpLURWsWeGY+ADyzavla4M7G53cC+yrOJUlaQ6t74C/NzKcAGh9fcrYHRsT+iJiMiMm5ubkWDydJWq3jJzEz82BmDmfm8MDAQKcPJ0lbRqsF/r2IuBig8fFEdZEkSevRaoH/C3Bj4/MbgX+uJo4kab3WcxnhIeBBYCgijkXETcCHgSsj4uvAlY2vJUldtOa9UDLz+rN8a2/FWSRJTfCdmJJUKAtckgplgUtSobwfuLacw1OzDlLQpmCBa0s5PDXLgfHp0xPhZ+cXODA+DWCJqzhuoWhLGZuYOV3eSxZOnmJsYqamRFLrLHBtKcfnF5palzYyC1xbyo7+vqbWpY3MAteWMjoyRF9vz4q1vt4eRkeGakoktc6TmNpSlk5UehWKNgMLXFvOvj2DFrY2BbdQJKlQFrgkFcoCl6RCWeCSVCgLXJIKZYFLUqEscEkqVFsFHhHvjojHIuJoRByKiPOqCiZJOreWCzwiBoE/AYYz81KgB3hLVcEkSefW7hbKNqAvIrYB24Hj7UeSJK1HywWembPAR4AngaeA72fmfasfFxH7I2IyIibn5uZaTypJWqGdLZQLgWuBlwM7gPMj4obVj8vMg5k5nJnDAwMDrSeVJK3QzhbKG4FvZuZcZp4ExoHXVRNLkrSWdu5G+CTwmojYDiwAe4HJSlJp03GQsFS9lgs8Mx+KiHuAh4HngCngYFXBtHk4SFjqjLauQsnMD2bmL2fmpZn51sz8SVXBtHk4SFjqDN+JqY5zkLDUGRa4Os5BwlJnWODqOAcJS53hTEx1nIOEpc6wwNUVDhKWqucWiiQVygKXpEJZ4JJUKAtckgplgUtSoSxwSSqUBS5JhbLAJalQFrgkFcoCl6RCWeCSVCgLXJIK1VaBR0R/RNwTEU9ExOMR8dqqgkmSzq3duxF+DLg3M/8gIl4IbK8gkyRpHVou8Ii4ALgC+EOAzHwWeLaaWJKktbSzhfIKYA74RERMRcStEXH+6gdFxP6ImIyIybm5uTYOJ0larp0C3wZcDnw8M/cAPwLev/pBmXkwM4czc3hgYKCNw0mSlmunwI8BxzLzocbX97BY6JKkLmi5wDPzu8B3ImJpMu1e4GuVpJIkrandq1DeCdzVuALlG8Db2o8kSVqPtgo8Mx8BhivKohYdnpp14ru0BTmVvnCHp2Y5MD7NwslTAMzOL3BgfBrAEpc2Od9KX7ixiZnT5b1k4eQpxiZmakokqVss8MIdn19oal3S5mGBF25Hf19T65I2Dwu8cKMjQ/T19qxY6+vtYXRk6Cz/haTNwpOYhVs6UelVKNLWY4FvAvv2DFrY0hbkFookFcoCl6RCWeCSVCgLXJIKZYFLUqEscEkqlJcRdol3DJRUNQu8C7xjoKROcAulC7xjoKROsMC7wDsGSuoEC7wLvGOgpE5ou8AjoicipiLic1UE2oy8Y6CkTqjiJObNwOPABRU816bkHQMldUJbBR4RO4E3Ax8C3lNJok3KOwZKqlq7WygfBd4HPH+2B0TE/oiYjIjJubm5Ng8nSVrScoFHxNXAicw8cq7HZebBzBzOzOGBgYFWDydJWqWdV+CvB66JiG8BdwNviIhPVZJKkrSmlgs8Mw9k5s7M3A28BfhiZt5QWTJJ0jl5HbgkFaqSe6Fk5peAL1XxXJKk9fEVuCQVygKXpEJZ4JJUKAtckgplgUtSoSxwSSqUBS5JhbLAJalQFrgkFcoCl6RCWeCSVCgLXJIKZYFLUqEscEkqlAUuSYWywCWpUBa4JBXKApekQrU8Ui0iLgE+CbwMeB44mJkfqyrYcoenZhmbmOH4/AI7+vsYHRli357BThxKkorRzkzM54A/zcyHI+JFwJGIuD8zv1ZRNmCxvA+MT7Nw8hQAs/MLHBifBrDEJW1pLW+hZOZTmflw4/P/BR4HKm/UsYmZ0+W9ZOHkKcYmZqo+lCQVpZI98IjYDewBHjrD9/ZHxGRETM7NzTX93MfnF5pal6Stou0Cj4ifAz4DvCszf7D6+5l5MDOHM3N4YGCg6eff0d/X1LokbRVtFXhE9LJY3ndl5ng1kVYaHRmir7dnxVpfbw+jI0OdOJwkFaOdq1ACuA14PDP/prpIKy2dqPQqFElaqZ2rUF4PvBWYjohHGmt/npmfbz/WSvv2DFrYkrRKywWemV8GosIskqQm+E5MSSqUBS5JhbLAJalQFrgkFSoys3sHi5gDvt3GU1wEPF1RnKpsxExgrmaZqznmak67uX4hM3/mnZBdLfB2RcRkZg7XnWO5jZgJzNUsczXHXM3pVC63UCSpUBa4JBWqtAI/WHeAM9iImcBczTJXc8zVnI7kKmoPXJL0U6W9ApckNVjgklSoDV/gEXF7RJyIiKN1Z1kuIi6JiH+NiMcj4rGIuLnuTAARcV5EfDUiHm3k+su6My0XET0RMRURn6s7y5KI+FZETEfEIxExWXeeJRHRHxH3RMQTjZ+z126ATEONv6elPz+IiHfVnQsgIt7d+Jk/GhGHIuK8ujMBRMTNjUyPVf13teH3wCPiCuCHwCcz89K68yyJiIuBi5cPdQb2VT3UuYVcAZyfmT9sDNz4MnBzZn6lzlxLIuI9wDBwQWZeXXceWCxwYDgzN9QbQCLiTuDfM/PWiHghsD0z5+vOtSQieoBZ4Dcys5036FWRZZDFn/VfzcyFiPg08PnMvKPmXJcCdwOvBp4F7gXenplfr+L5N/wr8Mx8AHim7hyrdWuoc7Ny0Q8bX/Y2/myI/0tHxE7gzcCtdWfZ6CLiAuAKFoemkJnPbqTybtgL/Hfd5b3MNqAvIrYB24HjNecB+BXgK5n548x8Dvg34PerevINX+AlONdQ5zo0tikeAU4A92fmhsgFfBR4H/B83UFWSeC+iDgSEfvrDtPwCmAO+ERjy+nWiDi/7lCrvAU4VHcIgMycBT4CPAk8BXw/M++rNxUAR4ErIuLFEbEd+D3gkqqe3AJv01pDneuQmacy89eBncCrG7/G1SoirgZOZOaRurOcwesz83LgTcA7Gtt2ddsGXA58PDP3AD8C3l9vpJ9qbOlcA/xT3VkAIuJC4Frg5cAO4PyIuKHeVJCZjwN/BdzP4vbJo8BzVT2/Bd6Gbgx1bkfjV+4vAVfVHAUWR/Bd09hvvht4Q0R8qt5IizLzeOPjCeCzLO5X1u0YcGzZb0/3sFjoG8WbgIcz83t1B2l4I/DNzJzLzJPAOPC6mjMBkJm3ZeblmXkFi9vBlex/gwXesm4NdW5WRAxERH/j8z4Wf7CfqDcVZOaBzNyZmbtZ/NX7i5lZ+yukiDi/cRKaxhbF77L4a2+tMvO7wHciYqixtBeo9QT5KtezQbZPGp4EXhMR2xv/NveyeF6qdhHxksbHXcB1VPj31s5Q466IiEPAbwMXRcQx4IOZeVu9qYAuDnVu0sXAnY0rBF4AfDozN8wlexvQS4HPLv6bZxvwD5l5b72RTnsncFdju+IbwNtqzgNAYy/3SuCP686yJDMfioh7gIdZ3KKYYuO8rf4zEfFi4CTwjsz8n6qeeMNfRihJOjO3UCSpUBa4JBXKApekQlngklQoC1ySCmWBS1KhLHBJKtT/A3FiZj1ZgpW0AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.scatter(x, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That resembles a straight line! Let's say the line is described by\n", "\n", "$$y = ax + b$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In regression, we are trying to find $a$ and $b$ given lots of values for $y$ and $x$, so we have something like\n", "\n", "\\begin{align}\n", "y_0 &= ax_0 + b \\\\\n", "y_1 &= ax_1 + b \\\\\n", "y_2 &= ax_2 + b \\\\\n", "\\vdots &= \\vdots \\\\\n", "\\end{align}\n", "\n", "Which we factor as \n", "\n", "$$\\underbrace{\\begin{bmatrix}y_0\\\\y_1\\\\\\vdots \\end{bmatrix}}_Y = \\underbrace{\\begin{bmatrix}x_0 & 1 \\\\x_1 & 1\\\\\\vdots & \\vdots \\end{bmatrix}}_X \\underbrace{\\begin{bmatrix}a\\\\b\\end{bmatrix}}_\\beta$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, when there are many points, $X$ is taller than it is wide and as such we know that there is no solution to the exact equations. Instead, we may try to get \"close\" to the solution. We can write the residuals as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\\epsilon = Y - X\\beta$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, $\\epsilon$ is a vector of errors, one for each data point. The (euclidian) length of this vector is given by \n", "\n", "$$||\\epsilon|| = \\sqrt{\\sum_i \\epsilon_i^2}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is common to focus on minimising the part without the square root to make calculations simpler. This leads to the popular way of determining the error of a fit called the \"sum of square errors\", sometimes expressed as $||\\epsilon||^2$.\n", "\n", "This minimisation is written in mathematical notation as \n", "\n", "$$\\min_\\beta \\sum_i \\epsilon_i^2$$\n", "\n", "which is read in words as as \"minimise (_with respect to_ or _by changing_ $\\beta$) the sum of the square error\"." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create the design matrices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The matrices $Y$ and $X$ are sometimes called the design matrices. We can build them using basic numpy functions as follows:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "Y = numpy.asmatrix(y).T\n", "X = numpy.bmat([numpy.c_[x], numpy.ones_like(Y)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "_Note_ `numpy.c_` produces a two dimensional array from a one dimensional one in a *c*olumn." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the case of polynomials, $X$ is a special matrix called a [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix), and numpy has a function which generates them more easily." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "X = numpy.asmatrix(numpy.vander(x, 2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is also a library called [patsy](https://patsy.readthedocs.io) which supplies a simplified syntax to construct these matrices:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "import patsy" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "Y, X = patsy.dmatrices(\"Q('Power Generated') ~ Q('Fuel Flow Rate')\", df)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "Y, X = map(numpy.asmatrix, (Y, X))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Pseudoinverse solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, let's apply the pseudoinverse method directly (note **you should never do this for production code**, as calculating inverses is computationally expensive)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The solution minimising the sum of the squares of the residual $||\\epsilon||^2$ can be shown to be \n", "\n", "$$ \\hat{\\beta} = (X^TX)^{-1}X^T Y $$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "matrix([[0.07854918],\n", " [1.85932397]])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Excel: =MMULT((MINVERSE(MMULT(TRANSPOSE(_X); _X))); MMULT(TRANSPOSE(_X); _Y))\n", "betahat = (X.T * X).I * X.T * Y\n", "betahat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Note:* The code above is as close as possible to the equation above, as I have made `X` and `Y` matrices. **The numpy developers advise against using the `numpy.matrix` class**.\n", "\n", "Matrix properties: \n", "* `.T`: Transpose\n", "* `.I`: Inverse\n", "* `.A`: Array form of matrix\n", "\n", "Normal `numpy.array`s don't have an `.I` property and don't multiply matrix-fasion but rather element-wise. Here is how we would have to write the code if we used arrays:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "Y = Y.A\n", "X = X.A" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.07854918],\n", " [1.85932397]])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "numpy.linalg.inv(X.T @ X) @ X.T @ Y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `@` operator always does matrix multiplication and expects two-dimensional arrays on both sides." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Dedicated solvers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculating the inverse is not a numerically well behaved operation, so you should rather use a dedicated routine if you are solving this kind of problem:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.07854918],\n", " [1.85932397]])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beta, residuals, rank, s = numpy.linalg.lstsq(X, Y, rcond=None)\n", "beta" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, polynomial fits are such a common operation that there are nicer routines to do this fit. There are a whole range of functions in numpy which start with poly. Press tab to see them." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "numpy.poly" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1.85932397, 0.07854918])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "poly = numpy.polyfit(x, y, 1)\n", "poly" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that we could just pass the data directly, and the routine handled building the various matrices and the fitting itself." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is useful to plot the regression with the data points, but we should sample on a finer grid." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "smoothx = numpy.linspace(min(x), max(x), 1000)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deVxU9f7H8deXRUFRMDUVl1BTU1EQZljUtLLUFk3NUnPJ2+160zJbJLP7K8vsZmlaaWpmZe65EO1hm5klCAiIG26hArkLyibLfH9/gFw1F2QGzgx8no+HD2HmcM5blreH75zz/SqtNUIIIRyPk9EBhBBClI8UuBBCOCgpcCGEcFBS4EII4aCkwIUQwkG5VObBGjRooH18fCrzkEII4fDi4uJOaK0bXvp4pRa4j48PsbGxlXlIIYRweEqpg5d7XIZQhBDCQUmBCyGEg5ICF0IIB1WpY+CXU1BQQGpqKnl5eUZHEcIqbm5uNGvWDFdXV6OjiGrC8AJPTU2lTp06+Pj4oJQyOo4Q5aK15uTJk6SmptKyZUuj44hqwvACz8vLk/IWDk8pRf369Tl+/LjRUYQdiYhPY0ZkMukZuXh7uRPWpx0DujS12f4NL3BAyltUCfJ9LC4UEZ/G5PAkcguKAEjLyGVyeBKAzUpcXsQUQogKMCMymdyCIjqpA/zb+SsAcguKmBGZbLNjSIEDHh4eAKSnpzN48GCD09iHrl27Gh1BCId2IiOTF1xWEFHjJUa7RFKHHADSM3Jtdgy7GEKxF97e3qxdu7ZCj1FYWIiLy+U/7Vd7rqyKiopwdna2ah8Af/zxh9X7EKK62rn5O76v+QIt1RFWFt7OG4UPc5ZaAHh7udvsOHIGfoGUlBR8fX0BWLx4MYMGDaJv3760adOG559/vnS79evXExoaSkBAAA8++CBZWVkATJ06FbPZjK+vL2PGjOH8ake33XYbL774Ij179uTdd9+96JivvPIKY8aMoXfv3owaNYqioiLCwsIwm8107tyZDz74AACLxcK4cePo2LEj9913H/fcc0/pfzY+Pj5MnTqV7t27s2bNGvbv30/fvn0JDAzk1ltvZffu3QCsWbMGX19f/Pz86NGjBwA7duwgKCgIf39/OnfuzN69e4H//VaitSYsLAxfX186derEZ599BsCGDRu47bbbGDx4MLfccgvDhw9HVncS1d3ZzFNEz3mEDpFDcVUWHin8D5ML/8UZagPg7upMWJ92NjueXZ2Bv/rVDnamn7HpPjt412VKv47l+tiEhATi4+OpWbMm7dq1Y/z48bi7uzNt2jR+/PFHateuzZtvvsmsWbN4+eWXefLJJ3n55ZcBGDlyJF9//TX9+vUDICMjg19//fWyx4mLi2PTpk24u7uzcOFCPD09iYmJ4dy5c3Tr1o3evXsTFxdHSkoKSUlJHDt2jPbt2/Poo4+W7sPNzY1NmzYB0KtXLxYsWECbNm2Ijo5m3Lhx/Pzzz0ydOpXIyEiaNm1KRkYGAAsWLGDChAkMHz6c/Px8ioqKLsoWHh5OQkICiYmJnDhxArPZXFr+8fHx7NixA29vb7p168bvv/9O9+7dy/W5FsLRJf6yhka/voBZnySq0RA6j5rBwL1Z7KvqV6HYq169euHp6QlAhw4dOHjwIBkZGezcuZNu3boBkJ+fT2hoKAC//PILb731Fjk5OZw6dYqOHTuWFviQIUOueJz+/fvj7l78a9X69evZtm1b6dl1ZmYme/fuZdOmTTz44IM4OTnRuHFjbr/99ov2cX7/WVlZ/PHHHzz44IOlz507dw6Abt26MXr0aB566CEGDRoEQGhoKK+//jqpqakMGjSINm3aXLTfTZs2MWzYMJydnWnUqBE9e/YkJiaGunXrEhQURLNmzQDw9/cnJSVFClxUOxknjrB36XjMmes56NScPfesI8TUC4ABXTxtWtiXsqsCL++ZckWpWbNm6dvOzs4UFhaiteauu+5i5cqVF22bl5fHuHHjiI2NpXnz5rzyyisX3V1au3btKx7nwue01syZM4c+ffpctM0333xz1azn92GxWPDy8iIhIeFv2yxYsIDo6Gi++eYb/P39SUhI4OGHHyY4OJhvvvmGPn36sGjRIu64446L8lzJ5T4/QlQX2mJha+Sn+ERPwV9nsbn5owSMeJ2abrUqLYOMgV+nkJAQfv/9d/bt2wdATk4Oe/bsKS3rBg0akJWVVe4XQ/v06cP8+fMpKCgAYM+ePWRnZ9O9e3fWrVuHxWLh6NGjbNiw4bIfX7duXVq2bMmaNWuA4gJOTEwEYP/+/QQHBzN16lQaNGjA4cOHOXDgAK1ateKpp56if//+bNu27aL99ejRg88++4yioiKOHz/Oxo0bCQoKKte/TYiq4kT6QRLe7kdg9NOccm7IocHfEvrY7Eotb7CzM3BH0LBhQxYvXsywYcNKhyamTZtG27Zt+de//kWnTp3w8fHBbDaXa/+PPfYYKSkpBAQEoLWmYcOGRERE8MADD/DTTz/h6+tL27ZtCQ4OLh3eudTy5csZO3Ys06ZNo6CggKFDh+Ln50dYWBh79+5Fa02vXr3w8/Nj+vTpLFu2DFdXVxo3blw6hn/ewIED2bx5M35+fiileOutt2jcuHHpC6NCVCfaYiHmi7nckjidDjqfqNZPYRr2Ei6uNQzJo6515YBS6mPgPuCY1tr3gsfHA08ChcA3Wuvnr7CLUiaTSV+6oMOuXbto3759OaJXP1lZWXh4eHDy5EmCgoL4/fffady4sdGxxAXk+7nqSk9J5sSqsXTOi2OXa0c8HppP8zZ+lXJspVSc1tp06eNlOQNfDMwFllyws9uB+4HOWutzSqkbbRVUXNl9991HRkYG+fn5vPTSS1LeQlQCS1ERW1a/Sefd7+CJIrrDi5gHT8TJBvdbWOuaBa613qiU8rnk4bHAdK31uZJtjtk+mrjUlca9hRAV42ByAjlrxxJSsJNt7iYaDptH8E22u47bWuUdA28L3KqUeh3IAyZqrWMut6FSagwwBqBFixblPJwQQlSMy80YeG/HBsSufJWAAwvJUzWI8f8vpv5jUU72dd1HeQvcBagHhABmYLVSqpW+zIC61nohsBCKx8DLG1QIIWztcjMGLl73Be2//JBQfYCtdXrQYsT7mBvb58lneQs8FQgvKewtSikL0ACQyZCFEA7j/IyBADXJZ7zL5zzu/BWnLXWI7/oeAX0eMTjh1ZW3wCOAO4ANSqm2QA3ghM1SCSFEJTg/M2CA2sNbrgu52SmdNYU9mFY4gsQ+V7572l5cc0BHKbUS2Ay0U0qlKqX+CXwMtFJKbQdWAY9cbvjEUch0sn9XEdPJvvLKK8ycOdPm+y2Pl19+mR9//PGq22zYsEFmZaziWtbVTHH5lLU1XsVN5TMqfxJhhY/j4dXQ6GhlUparUIZd4akRNs5iOJlO9n+qenFNnTr1mtts2LABDw8PmRu9ikra+DlL88No6nKcxYW9mVE4hGzcbT5jYEWyr5dUDSbTyVbedLK2yHghDw8PnnvuOQICAujVq1fp2pQJCQmEhITQuXNnBg4cyOnTpwEYPXr0RZ+/KVOmEBAQQKdOndi9ezcpKSksWLCA2bNn4+/vz2+//XbZbMLxZJ46wZZ3htHp59FYlCurfBfyocdYcnCnqZc7bwzqVKETUNmSfd1K/90LcCTJtvts3Anunl6uD5XpZCtuOtkxY8bYLCNAdnY2AQEBvP3220ydOpVXX32VuXPnMmrUKObMmUPPnj15+eWXefXVV3nnnXf+9vENGjRg69atzJs3j5kzZ7Jo0SIef/xxPDw8mDhxIgCdOnX6WzbhWOLXL6PZH/9HgM5kc9NRdBk5naHutRlqdLBysq8CtzMynWzFTCdr64wATk5OpZ+DESNGMGjQIDIzM8nIyKBnz54APPLIIxcd80LnjxUYGEh4ePhlt7lcNuEYThw5zKFlTxCQ9Sv7nVuS2X8ZoX6OP/WxfRV4Oc+UK4pMJ1sx08naOuPlXO8K8ef/LVf7d1wuW/369a/rOKJyaYuF2K8+oE38NHx1HptbjsX08Ku41qh57Q92ADIGfp1kOtlrTyc7efJkPv/88ys+b+uMUPyfwvnP+YoVK+jevTuenp7Uq1eP3377DYClS5eWno2XRZ06dTh79mzp+5fLJuzXkcP72DajD+b4Fzji0py/hv1A6OjpVaa8wd7OwB2ATCd77elkk5KS6N+//1X/nbbMCMW/gezYsYPAwEA8PT1LX2z99NNPefzxx8nJyaFVq1Z88skn1/oSlOrXrx+DBw/miy++YM6cOcyePftv2YT9sRQVsWXdLHx3vE1dLES1C8P80As4W3mFlz265nSytiTTyVrHUaaT7dOnD5GRkZV6TA8Pj9KrgYwk38/GOrwviTOrx9IxP4ntNf25YegHeLe8xehYVrNmOllhJxxlOtnKLm8hCgvyiV01Df998/BUrsR0noppwHi7m3zK1qTAHYhMJ3tl9nD2LYxxYHs0RRFPEFK4l/jaXWk2Yj5mbx+jY1UKKXAhhEM6l5fD1uUvYTr0CWdVbeLMbxNw96NV/qz7QlLgQgiHkxz7MzW/nUCo5RCxnnfReuQcAhs2MTpWpZMCF0I4jNzssyQuCcN8ZBUn1A0k9vgA0x2Oeh+l9aTAhRAOYfvvX1Hvx+cI0UeJbnA/HUa9g5/nDUbHMlT1GSwqo2tNeRoREcHOnTsrMZEQ1duZjJNEvzcS3x9GoFHsuGsFweOXUKealzdIgV83KXAhKkdEfBoTpr5B9mwTppNf8V3dB6k/MZaO3e41OprdcLgCj4hPo9v0n2n5wjd0m/4zEfFpVu/z9ddfp127dtx5550kJycD8OGHH2I2m/Hz8+OBBx4gJyeHP/74gy+//JKwsDD8/f3Zv3//ZbcTQlhnza/xOIU/xruW6WTq2gzKf5VnTw8mcs8Zo6PZFYcq8PMLkKZl5KIpXoB0cniSVSUeFxfHqlWriI+PJzw8nJiYGKB4drqYmBgSExNp3749H330EV27dqV///7MmDGDhIQEWrdufdnthBDloy0WYr/5kDt+7kdfp2hmFzxAv/zXSdQ3k1tQxIzIZKMj2pWyLKn2sVLqWMnyaZc+N1EppZVSDSom3sUuXID0PGu/qL/99hsDBw6kVq1a1K1bt3QOj+3bt3PrrbfSqVMnli9fzo4dOy778WXdTghxdcfS/iRh5r2YYiZyWDfkvvz/8m7RAxRccK3F+TUsRbGynIEvBvpe+qBSqjlwF3DIxpmu6EpfPGu/qJebenT06NHMnTuXpKQkpkyZctHUsOXZTghxedpiYcvaWbgvDKV9dgxRNz/Dk25vskc3/9u23l7uBiS0X9cscK31RuDUZZ6aDTwPVNpsWFf64lnzRe3Roweff/45ubm5nD17lq+++gqAs2fP0qRJEwoKCli+fHnp9pdOMXql7YQQ15Z2YBc73rydoO2vcqjmzZwYuYGQEa8w8e6OuLtevLarI61VWVnKNQaulOoPpGmtE22c56rC+rSz+Rc1ICCAIUOG4O/vzwMPPMCtt94KwGuvvUZwcDB33XUXt9zyv9nMhg4dyowZM+jSpQv79++/4nZCiCsrKiwkasVr1Pu0Jz55yUR3fIn2k36l2c3Fa9IO6NKUNwZ1oqmXOwocbq3KylKm6WSVUj7A11prX6VULeAXoLfWOlMplQKYtNYnrvCxY4AxAC1atAg8ePDgRc9f7/SbEfFpzIhMJj0jF28vd8L6tJMvqrAbMp3stR3cFUfeunG0K9xNonsQjR6eT+PmNxsdy67ZcjrZ1kBLILFk7LgZsFUpFaS1PnLpxlrrhcBCKJ4PvBzHu8iALk2lsIVwQAX554hd/jKBKYvIVu7EBrxJ4H1jqtXkU7Z23QWutU4Cbjz//rXOwIUQYm/Cbzh/+SShlhTi6t6Oz4i5mBo1MzqWw7tmgSulVgK3AQ2UUqnAFK21TS921lpf9yK0QtibylzdylHk5WQRv3QSQenLOaW8iO/6PoG9Rxgdq8q4ZoFrrYdd43kfawK4ublx8uRJ6tevLyUuHJbWmpMnT+Lm5mZ0FLuxM+p76kQ+Q6hOZ8sN99Ju1Ht0qVcpt4xUG4bPRtisWTNSU1M5fvy40VGEsIqbmxvNmsmwQNaZ0+xY8izBJ8JJV43Y3msJQbfeb3SsKsnwAnd1daVly5ZGxxBC2MC2X9Zy46+TMOuTRDUaQudRM/D28DQ6VpVleIELIRxf5smj7FkyHnNmJAedmrPn7rWEmO80OlaVJwUuhLDK1u8+4aboKfjrLDY3f5SAEa9T062W0bGqBSlwIUS5nEg/yKHlTxCQ/Rv7nFuTcf8qQjt3NTpWtSIFLoS4LtpiIeaL97kl8Q066nw2t34K87CXcHGtYXS0akcKXAhRZukpyZxYNY6gvFh2uXak9oPzCG3rb3SsaksKXAhxTZaiImLWvEWnXbPxRBHdYTLmwWE4OTtf+4NFhZECF0Jc1aE9CWSvGUtwwU62uZtoOGwewTfJtK72QApcCHFZBQX5xK54lYADH5CnahDj/zqm/uNk8ik7IgUuhPib/dv+QH/xJKFF+9nq0YMWI9/H3LiF0bHEJaTAhRCl8nKziV/2H0ypS8hUddga8i4BfUcbHUtcgRS4EAKA3TE/4v7dBEItqcR49aXtqPcIqN/I6FjiKqTAhajmss9mkLQkjKBjazimGpB428eYb3vA6FiiDKTAhajirrYMYdLGL6j/y0RC9DGiGw6i46hZNK5bz+DEoqykwIWowiLi05gcnkRuQREAaRm5TA5PoiD7FDfF/pegjG85rLzZ2fczgkP6GpxWXC8pcCGqsBmRyaXlfV73omh6/PAx9TnDZu9RdBn5Bs1reRiUUFhDClyIKiw9I7f07fpk8orrp/RzjmKn5SbODFpGqP+tBqYT1rrmFflKqY+VUseUUtsveGyGUmq3UmqbUupzpZRXxcYUQpSHt5c7oLnfaRM/1Ayjt1MsMwse5N9uM2gj5e3wynJL1WLg0sGxHwBfrXVnYA8w2ca5hBA2MMHszseuM3m3xjxSdGPuzf8vHzkN5rm7fY2OJmygLIsab1RK+Vzy2PoL3o0CBts2lhDCGpaiImLWzeaeHTNxcrLwJqP5IP9OmnjV5o0LrkIRjs0WY+CPAp9d6Uml1BhgDECLFnIrrhAV7fC+7ZxZ/TjB+Ulsd/On3pAFTGrVnklGBxM2Z1WBK6X+AxQCy6+0jdZ6IbAQwGQyaWuOJ4S4ssKCAmI/m4b/3vfxwoUtnV/BPHCCTD5VhZW7wJVSjwD3Ab201lLMQhjowM4tFIY/QUjhHuJrd6Xp8HkENW1pdCxRwcpV4EqpvsAkoKfWOse2kYQQZZV/Lo+4Zf9H4KGPyVK1iTO/TcDdj8pZdzVxzQJXSq0EbgMaKKVSgSkUX3VSE/hBKQUQpbV+vAJzCiEukbz1V2p8PZ5Qy0FiPe+k9ci5BDZsYnQsUYnKchXKsMs8/FEFZBFClEFu9lkSl4RhPrKKk6oeCbd+gKnXUKNjCQPInZhCOJAdv3+D54/PEqKPEN3gftqPnI2/V32jYwmDSIEL4QDOZJ5i16dPE3zqC1JVY7bftYzgbv2MjiUMJgUuhJ1L/PkzGm+cjEmfIqrxMPxGzaBZ7TpGxxJ2QApcCDt1+vhf7F/6JKYzP5Li1IL9935ESODtRscSdkQKXAg7oy0W4r77mFYxr+Kns9nc4l8EDn+NGm7uRkcTdkYKXAg7cjw9hdRlYzHl/MFelzZkDnyf0I7BRscSdkoKXAg7oC0WYj5/j1uS3qS9LiCqzdOYhvwHF9caRkcTdkwKXAiDnF+rUmUe5C3XRXR12s7OGp2o89B8Qm7uZHQ84QDkflshDBARn8aL6xLofTac9TUm0Unt5+XCR9ndZwXNpbxFGckZuBAGWPHNDyx1mkug015+KfLjxYLH+Iv6/PTDPgYFyrTLomykwIWoRAX554hdPoWlBR+So9x4On8cEZZugAIuXsNSiGuRAheikuxL3IT68klCi/7kOx3KS/mjOIHnRdsUr2EpRNlIgQtRwfJysohfOhlz+jIyVF3iu77PuYa3kx2eBAVFpdu5uzoT1qedgUmFo5ECF6IC7YyKpE7k04TqdGLq3UPbUe/R5YaGdCl5fkZkMukZuXh7uRMma1WK6yQFLkQFyDpzmu1LniPoeDhHnBqSdMdizD0GXrTNgC5NpbCFVaTAhbCxbRvWceOGSQTpE2y5cTCdRs3Eu46X0bFEFSQFLoSNZJ48yp4lT2HO/J5DTk3Z03c1IUG9jY4lqjApcCFsYOv3i2kRNYUu+gybm42my4j/4uZe2+hYooory5qYH1O8+vwxrbVvyWM3AJ8BPkAK8JDW+nTFxRTCPp04cohDS58gIHsj+51bkXn/CkI7dzM6lqgmynIr/WKg7yWPvQD8pLVuA/xU8r4Q1Ya2WIiJmIvrghA6Zm0mquUTtJgURWspb1GJyrKo8UallM8lD99P8Ur1AJ8CG4BJNswlhN3662Ayx1eOw5wXy27XDrgPnk9IO3+jY4lqqLxj4I201n8BaK3/UkrdeKUNlVJjgDEALVrIHA/CcVmKiohZOwPfnbPxRBPdfhLmByfh5OxsdDRRTVX4i5ha64XAQgCTyaQr+nhCVIRDexPJXj2W4IIdJLkFUH/YAoJ95K5JYazyFvhRpVSTkrPvJsAxW4YSwl4UFuQTs/I1AvbP55yqwRa/aZjvfwLlJDMxC+OVt8C/BB4Bppf8/YXNEglhJ/YnRaEjxhFatJ94j+40Hz6PIO+bjI4lRKmyXEa4kuIXLBsopVKBKRQX92ql1D+BQ8CDFRlSiMp0Li+Hrcv+g+nwp5xRHmwNfoeAvqNBKaOjCXGRslyFMuwKT/WycRYhDLc75kfcv3uaUMthYrz60GbkewQ0aGx0LCEuS+7EFALIycpk25KJBB1dwzFVn8SeizDfLr9YCvsmBS6qve2/fcENP4cRoo8S3XAQHUe9TeO6NxgdS4hrkgIX1Vbm6RMkL5lA0OmvOay82dlnFcGhdxsdS4gykwIX1U5EfBobvvqUF4oWEkgG33oO4Y7HZ9G8lofR0YS4LnIxq6hWVv+6FZfwf/KO5U1Oaw8G5L/Gc6cH8X1yptHRhLhucgYuqgVtsRD3zYfcGTsVD6dcZhUMZn5RfwpwgYIiZkQmy+o4wuFIgYsq72jqfo4sH4spN5p4fTPPF4xhr2520TbpGbkGpROi/KTARZVlKSoiJvwdOmyfQVuKiGr7HBMPhpCaX/C3bb293A1IKIR1pMBFlZS6bzuZq8cSnL+N7W7+1BuygJBW7ZkYn8bk8CRyC4pKt3V3dSasj0xMJRyPFLioUooKC4lZNQ2/ve/jiTNbOk3BPOjp0smnzo9zz4hMJj0jF28vd8L6tJPxb+GQpMBFlfHnzhgKwscRUriHhNqheA+fT1DTln/bbkCXplLYokqQAhcOL/9cHnHLXyLw4Edkq9rEmmcSePc/ZcpXUeVJgQuHtmfrr7h8/RShlhTiPHvRcsQcTDfK2bWoHqTAhUPKzT5L4tJJmP9awUlVj4TuCwi880oTZwpRNUmBC4ez449v8fzhWUL0X2yp349bRr2Lv1d9o2MJUemkwIXDOJt5ip1LniH4ZARpqhHb71xKUPf+RscSwjBS4MIhJP68msYbX8CkTxHVeBh+I9+iqUddo2MJYSirClwp9QzwGKCBJOAfWus8WwQTAuD0iSPsW/Ik5jM/kOLUnH33LCLEdIfRsYSwC+W+zkop1RR4CjBprX0BZ2CorYKJ6k1bLMR9+xF6rhn/zJ+JavZPmjy/hXZS3kKUsnYIxQVwV0oVALWAdOsjieruePpBUpeNJTDnd/Y630zGwHmE+AYbHUsIu1PuAtdapymlZlK8Kn0usF5rvf7S7ZRSY4AxAC1atCjv4UQ1oC0WYiLmcMu26bTXBUTdPAHT0P/DxbWG0dGEsEvWDKHUA+4HWgLeQG2l1IhLt9NaL9Ram7TWpoYNG5Y/qajS0lOS2f5mL4K2vUxqjVYcH/ETISOnSnkLcRXWDKHcCfyptT4OoJQKB7oCy2wRTFQPRYWFxKx5k06738UTRXTH/2B+4DmcnJ2NjiaE3bOmwA8BIUqpWhQPofQCYm2SSlQLB3dvJXfdOEIKdrHN3cyND88nuEUbo2MJ4TCsGQOPVkqtBbYChUA8sNBWwUTVEhGfVjqFa7O6LkysE8ndJz4lR7kR0+UNTP0el8mnhLhOVl2ForWeAkyxURZRRUVcsIhCR5XCW3kf0DH/IL+7dafdPxZgbtzc6IhCOCS5E1NUuBmRyVgKcglzCeffzl9zirr8O/8Ztrv14HcpbyHKTQpcVLjGmQksqbGQ1k5/saawB68VjuAMHihZSFgIq0iBiwqTdTaD7UueY02NdaRTn5H5L/CbpXPp87KQsBDWkQIXFSLp13Aa/PI8QfoEP9bpzwsZAzll+d813bKQsBDWkwIXNpV56jh7lozHnPEdh5yasqfPanoH9ybngqtQZCFhIWxDClzYzNbIJbTY/BJd9Bmimj2C/4g3cHOvDchCwkJUBClwYbXjRw5xaOkTBGZvZL9zKzL6ryDEr5vRsYSo8qTARblpi4XYLxfQJuF1Ouk8olo+QeDDU3CtUdPoaEJUC1LgolyOHNrL0RVjMefFsNu1A+6D5xPSzt/oWEJUK1Lg4rpYioqIWTsT352zqIsm6pZJmB98HmcX+VYSorLJT50os0N7t5G1ehzBBUlsd+vCDcM+IMRHLgUUwihS4OKaCgvyiV35Gv775+OlarDF7zXM9z8pk08JYTApcHFV+5OisUSMI6RoH/G1u9F8xHyCvG8yOpYQAilwcQXn8nLYuuz/MB1ezBnlQVzQOwT0fUTOuoWwI1Lg4m92x/6E27cTCLUcJsarN21GziGwQWOjYwkhLiEFLkrlZGWybcnzBB39jOPqBhJ7fIj5joeMjiWEuAIpcAHA9k1fUu+niYToo0Q3GECHUbPx87zB6FhCiKuQAq/mMjNOkvzpUwSd/ppU1YSdvVcS3PUeo2MJIcrAqgJXSnkBiwBfQAOPaq032yKYqHjxP66k6aYXCdSniWoyHP9Rb9GslsUCaGUAAA3FSURBVIfRsYQQZWTtGfi7wPda68FKqRpALRtkEhXs1LE0Diwdj+nsT/zp5ENmv8WEdOlpdCwhxHUqd4ErpeoCPYDRAFrrfCDfNrFERdAWC3HfLqJ17FQ66xw23/RvAodPpUZNN6OjCSHKwZoz8FbAceATpZQfEAdM0FpnX7iRUmoMMAagRYsWVhxOWONo2gH+WjYWU24Ue1zakvnAfELbm4yOJYSwgjV3ZbgAAcB8rXUXIBt44dKNtNYLtdYmrbWpYcOGVhxOlIe2WIheO4taC7vSLmcrUW2epfULm/GR8hbC4VlzBp4KpGqto0veX8tlClwYJ+3ADjJWjSU4P5EdNf3wGjqfkFYdjY4lhLCRche41vqIUuqwUqqd1joZ6AXstF00UV5FhYXEfPY6fnvmUhdntnSagnnQ03IbvBBVjLVXoYwHlpdcgXIA+If1kYQ1UnbFkr9uLCGFe0ioFUKT4fMIatba6FhCiApgVYFrrRMAGUw1WER8GrO+S2JA9hqedPmcLGoRa55B4D2PyVm3EFWY/HQ7uIj4NJat+5wP8ibyrOtavrcEcW/R26Q2vVfKW4gqTm6ld2B5OVmciniBz5y/4jhePJb/HD9aAgGYEZnMgC5NDU4ohKhIUuAOasfm7/Bc/wyPqr9YUXg70wsf5gy1S59Pz8g1MJ0QojJIgTuYs5mn2LnkGYJPRpCuGvG40xS+L/z7upTeXu4GpBNCVCYZJHUgib+sIXu2GfOJL4i6cQhez8XQt98Q3F2dL9rO3dWZsD6y2LAQVZ2cgTuAjBNH2Lt0PObM9aQ4NWfPPesIMfUCYEAXT6B4zDs9IxdvL3fC+rST8W8hqgEpcDumLRa2fr8Yny2v4K+ziGr+T7qMmEZNt4snfRzQpakUthDVkBS4nTqRfpBDy8YRmLOJvc43kzFgNSGdQoyOJYSwI1LgdkZbLMR8MZdbEqfTQecT1fopTMNewsW1htHRhBB2RgrcjqSnJHNy5eMEndvKTldf6jw0j5A2fkbHEkLYKSlwO2ApKmLL6jfpvPsdPFFEd3gR8+CJODk7X/uDhRDVlhS4wQ4mx5OzdhwhBTvZ5m6m4bD3Cb5JLgEUQlybFHgliYhPu+hSv+d6taTJzg8JOLCQPFWDGP//Yuo/VuYvEUKUmRR4JYiIT2NyeBK5BUUAeGXuou1Xz+DrlMLWOj1oMWIe5sbNDU4phHA0UuCVYEZkMrkFRdQkn/Eun/O481ecpg4T1URmTnzJ6HhCCAclBV4J0jNyCVB7eMt1ITc7pbOmsAevFY7gLB7MNDqcEMJhSYFXsOyzGUytuZThfE869RmVP4mNluJLA5vKhFNCCCtYXeBKKWcgFkjTWt9nfaSqI2nj59T/5XmGc5xllrt4s2AI2RSXtkw4JYSwli3OwCcAu4C6NthXlZB56jjJS54iKONbDjk1Jfmez6hbwxevyGRyZMIpIYSNWFXgSqlmwL3A68CzNknk4LZGLqPF5v8QoM+wuekouoycjpt7bdqDFLYQwqasPQN/B3geqGODLA7txJHDHFr2BAFZv7LfuRUZ/ZcT6tfd6FhCiCqs3HeNKKXuA45preOusd0YpVSsUir2+PHj5T2c3SqefGoeLgtC8D37O1E+T9BiUhQ3S3kLISqYNWfg3YD+Sql7ADegrlJqmdZ6xIUbaa0XAgsBTCaTtuJ4dufI4X0cXf445rwYdru0x33wPEJuCTA6lhCimih3gWutJwOTAZRStwETLy3vqspSVETMurfpuGMWdbEQdcvzmB+chLOLXJUphKg80jjX6fC+JM6uHktwfhJJbl2oP3QBIS1vMTqWEKIaskmBa603ABtssS97VViQT+yqafjvm4enciWm81RMA8bL5FNCCMPIGXgZHNgeTVHEE4QU7iW+dleajZiP2dvH6FhCiGpOCvwqzuXlsHX5S5gOfcJZVZu4oFkE9P2HnHULIeyCFPgVJMf+TI1vJxBqOUSs5120HjmHwIZNjI4lhBClpMAvkZt9loQlYQQdWcUJdQOJPT7AdMdQo2MJIcTfSIFfYPvvX1Hvx+cI1UeJbjCADqNm4+d5g9GxhBDisqTAgTMZJ9m15GmCT31JqmrMjrtWENztXqNjCSHEVVX7Ak/4aRXev03GpE8T1WQ4fiPfpFntaj+1ixDCAVTbAj91LI0DS8djOvsTfzrdRMZ9iwkJ6Gl0LCGEKLNqV+DaYiHu249oHTuVzjqbzTeNIXD4a9So6WZ0NCGEuC7VqsCPpf1J+vKxmHI2s8elLRmD5hHawWx0LCGEKJdqUeDaYiEm/B3aJ71FO4qIavMM5qH/J5NPCSEcWpVvsLQDOzm9aixB+QnsqNkZz4fmE3Kzr9GxhBDCalW2wIsKC4n57L/47ZmDJ85E+76MedDTODk7Gx1NCCFsokoWeMquOPLWjSOkcDeJtYJpPHw+wc1aGx1LCCFsqkoVeP65POJWTCEwZRHZyp3YwLcIvPdfMvmUEKJKqjIFvjd+I85fjSfUkkJc3TvwGTEHU6NmRscSQogK4/AFnpeTRfzSSQSlL+eU8iKh23wC73rY6FhCCFHhHLrAd27+jjrrnyVUp7PlhntpN+o9/Os1MDqWEEJUinIXuFKqObAEaAxYgIVa63dtFexCEfFpzIhMJj0jF28vd57p0YTmW98i+EQ46aoR23stIejW+yvi0EIIYbesOQMvBJ7TWm9VStUB4pRSP2itd9ooG1Bc3pPDk8gtKAKgzZnNhEZ+RBNOEdVoCJ1HzcDbw9OWhxRCCIdQ7gLXWv8F/FXy9lml1C6gKWDTAp8RmUxuQRFenOUl16U84LyJvZam/NN1Gp+Me9KWhxJCCIdikzFwpZQP0AWIvsxzY4AxAC1atLjufadn5BKg9vBBjVl4kc17hQOYWziQgnxX60ILIYSDs7rAlVIewDrgaa31mUuf11ovBBYCmEwmfb379/Zy51BGI3ZZbuKNwofZpW8CoKmXu3XBhRDCwVl1h4tSypXi8l6utQ63TaSLhfVpR7brDYwqmFxa3u6uzoT1aVcRhxNCCIdhzVUoCvgI2KW1nmW7SBcb0KUpwEVXoYT1aVf6uBBCVFfWDKF0A0YCSUqphJLHXtRaf2t9rIsN6NJUClsIIS5hzVUomwBlwyxCCCGug8zyJIQQDkoKXAghHJQUuBBCOCgpcCGEcFBS4EII4aCU1td9c2T5D6bUceCgFbtoAJywURxbscdMILmul+S6PpLr+lib6yatdcNLH6zUAreWUipWa20yOseF7DETSK7rJbmuj+S6PhWVS4ZQhBDCQUmBCyGEg3K0Al9odIDLsMdMILmul+S6PpLr+lRILocaAxdCCPE/jnYGLoQQooQUuBBCOCi7L3Cl1MdKqWNKqe1GZ7mQUqq5UuoXpdQupdQOpdQEozMBKKXclFJblFKJJbleNTrThZRSzkqpeKXU10ZnOU8plaKUSlJKJSilYo3Oc55SyksptVYptbvk+yzUDjK1K/k8nf9zRin1tNG5AJRSz5R8z29XSq1USrkZnQlAKTWhJNMOW3+u7H4MXCnVA8gClmitfY3Oc55SqgnQRGu9VSlVB4gDBmitbbqoczlyKaC21jqrZMWkTcAErXWUkbnOU0o9C5iAulrr+4zOA8UFDpi01nZ1A4hS6lPgN631IqVUDaCW1jrD6FznKaWcgTQgWGttzQ16tsjSlOLv9Q5a61yl1GrgW631YoNz+QKrgCAgH/geGKu13muL/dv9GbjWeiNwyugcl9Ja/6W13lry9llgF2D4qhO6WFbJu64lf+zif2mlVDPgXmCR0VnsnVKqLtCD4lWv0Frn21N5l+gF7De6vC/gArgrpVyAWkC6wXkA2gNRWuscrXUh8Csw0FY7t/sCdwRKKR+gCxBtbJJiJcMUCcAx4AettV3kAt4BngcsRge5hAbWK6XilFJjjA5TohVwHPikZMhpkVKqttGhLjEUWGl0CACtdRowEzgE/AVkaq3XG5sKgO1AD6VUfaVULeAeoLmtdi4FbiWllAfFCzs/rbU+Y3QeAK11kdbaH2gGBJX8GmcopdR9wDGtdZzRWS6jm9Y6ALgbeKJk2M5oLkAAMF9r3QXIBl4wNtL/lAzp9AfWGJ0FQClVD7gfaAl4A7WVUiOMTQVa613Am8APFA+fJAKFttq/FLgVSsaY1wHLtdbhRue5VMmv3BuAvgZHgeI1VPuXjDevAu5QSi0zNlIxrXV6yd/HgM8pHq80WiqQesFvT2spLnR7cTewVWt91OggJe4E/tRaH9daFwDhQFeDMwGgtf5Iax2gte5B8XCwTca/QQq83EpeLPwI2KW1nmV0nvOUUg2VUl4lb7tT/I2929hUoLWerLVuprX2ofhX75+11oafISmlape8CE3JEEVvin/tNZTW+ghwWCnVruShXoChL5BfYhh2MnxS4hAQopSqVfKz2Yvi16UMp5S6seTvFsAgbPh5s2ZV+kqhlFoJ3AY0UEqlAlO01h8ZmwooPqMcCSSVjDcDvKi1/tbATABNgE9LrhBwAlZrre3mkj071Aj4vPhnHhdghdb6e2MjlRoPLC8ZrjgA/MPgPACUjOXeBfzb6Cznaa2jlVJrga0UD1HEYz+31a9TStUHCoAntNanbbVju7+MUAghxOXJEIoQQjgoKXAhhHBQUuBCCOGgpMCFEMJBSYELIYSDkgIXQggHJQUuhBAO6v8Bc5NYkUd8GjEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def regplot(poly):\n", " smoothy = numpy.polyval(poly, smoothx)\n", " plt.scatter(x, y, label='data')\n", " plt.plot(smoothx, smoothy, label='linear regression')\n", " plt.plot(x, numpy.polyval(poly, x), label='linear regression, less points')\n", " plt.legend(loc='best')\n", "regplot(poly)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There's obviously no difference between the two for a linear fit, but what about higher orders?" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxM1//H8dfJngiC2BIiqKISJGJfizaoqrXoYuniW11/bUWri63aKlpU+60qSlvUUrSqaqnaq4hE7RRBEjuJRBLJzJzfH4l80agkM8mdST7Px8MDMzf3vjOJt5sz956jtNYIIYRwPE5GBxBCCJE/UuBCCOGgpMCFEMJBSYELIYSDkgIXQggH5VKYB/P19dWBgYGFeUghhHB4kZGRF7XW5W9/vFALPDAwkF27dhXmIYUQwuEppU7m9LgMoQghhIOSAhdCCAclBS6EEA6qUMfAc5KRkUFsbCxpaWlGRxHCKh4eHlSpUgVXV1ejo4hiwvACj42NpWTJkgQGBqKUMjqOEPmitebSpUvExsZSvXp1o+OIYsLwAk9LS5PyFg5PKUW5cuW4cOGC0VGEHVkeFcfE1YeJT0jFz8eTiPDadA/xt9n+DS9wQMpbFAnyfSxutjwqjhFL95KaYQYgLiGVEUv3AtisxOVNTCGEKADjVx0iNcNMBa7Q2elPAFIzzExcfdhmx5ACB7y9vQGIj4+nd+/eBqexDy1atDA6ghAO5+CZq4z+aT/tJv7O5atJPOf8E7+7v8ZHrjPwJgWA+IRUmx3PLoZQ7IWfnx9Lliwp0GOYTCZcXHJ+2f/tudwym804OztbtQ+Abdu2Wb0PIYoDrTVb/r7I1HVH2XXyCm7OTgz1P07Pa9OoxhnWmkN5z/QkyXgB4OfjabNjyxn4TWJiYggKCgJgzpw59OzZk06dOlGrVi2GDx+evd2aNWto3rw5oaGh9OnTh+TkZADGjh1L48aNCQoKYsiQIdxY7ahdu3a89dZbtG3blqlTp95yzNGjRzNkyBAefPBBBgwYgNlsJiIigsaNG1O/fn2+/PJLACwWC88//zz16tWja9eudOnSJfs/m8DAQMaOHUurVq1YvHgxx44do1OnTjRq1IjWrVtz6NAhABYvXkxQUBANGjSgTZs2AOzfv58mTZrQsGFD6tevz9GjR4H//VSitSYiIoKgoCCCg4NZuHAhABs2bKBdu3b07t2bOnXq8PjjjyOrO4ni5vDZJPrN2M6Ts3ZwJjGNCe292VdnNq+ef4uyJdwZYn6TZzOGcUpXBMDT1ZmI8No2O75dnYGPWbGfA/FXbbrP+/xKMerhevn62OjoaKKionB3d6d27dq89NJLeHp6Mm7cONatW0eJEiX46KOP+OSTTxg5ciQvvvgiI0eOBODJJ5/k559/5uGHHwYgISGBjRs35nicyMhItmzZgqenJzNmzKB06dLs3LmT69ev07JlSx588EEiIyOJiYlh7969nD9/nrp16/LUU09l78PDw4MtW7YA0KFDB6ZPn06tWrX4888/ef7551m/fj1jx45l9erV+Pv7k5CQAMD06dN55ZVXePzxx0lPT8dsNt+SbenSpURHR7Nnzx4uXrxI48aNs8s/KiqK/fv34+fnR8uWLdm6dSutWrXK12sthCNJTTczbf1RZmw6TkkPF97vUp2+1xfjsv0zcHaDjmMo2ex5uuy9wP6ifhWKverQoQOlS5cG4L777uPkyZMkJCRw4MABWrZsCUB6ejrNmzcH4Pfff2fChAmkpKRw+fJl6tWrl13gffv2veNxunXrhqdn5o9Va9as4a+//so+u05MTOTo0aNs2bKFPn364OTkRKVKlbj//vtv2ceN/ScnJ7Nt2zb69OmT/dz169cBaNmyJYMGDeLRRx+lZ8+eADRv3pz333+f2NhYevbsSa1atW7Z75YtW+jfvz/Ozs5UrFiRtm3bsnPnTkqVKkWTJk2oUqUKAA0bNiQmJkYKXBR5e2MTefn7KE5cvEbvUH9G1ziM98ZXICke6veFjmOgVGUg82oTWxb27eyqwPN7plxQ3N3ds//s7OyMyWRCa80DDzzAggULbtk2LS2N559/nl27dlG1alVGjx59y92lJUqUuONxbn5Oa820adMIDw+/ZZuVK1f+a9Yb+7BYLPj4+BAdHf2PbaZPn86ff/7JypUradiwIdHR0Tz22GM0bdqUlStXEh4ezsyZM2nfvv0tee4kp9dHiKLKYtHM2nKCCasPUa6EO8t6lSZk31vw81aoVB/6fA0BzQo1k4yB51GzZs3YunUrf//9NwApKSkcOXIku6x9fX1JTk7O95uh4eHhfPHFF2RkZABw5MgRrl27RqtWrfjhhx+wWCycO3eODRs25PjxpUqVonr16ixevBjILOA9e/YAcOzYMZo2bcrYsWPx9fXl9OnTHD9+nBo1avDyyy/TrVs3/vrrr1v216ZNGxYuXIjZbObChQts2rSJJk2a5OtzE8JRXbmWzqA5O3n/l4M8dI8HG+r9TMgvD8P5g9B1MgzZUOjlDXZ2Bu4Iypcvz5w5c+jfv3/20MS4ceO49957efbZZwkODiYwMJDGjRvna//PPPMMMTExhIaGorWmfPnyLF++nF69evHbb78RFBTEvffeS9OmTbOHd243b948hg4dyrhx48jIyKBfv340aNCAiIgIjh49itaaDh060KBBA8aPH893332Hq6srlSpVyh7Dv6FHjx788ccfNGjQAKUUEyZMoFKlStlvjApR1B05l8TTc3dy4WoqC0MP0uTE56hTCRD2NNz/FniVNSybutuVA0qp2UBX4LzWOuimx18CXgRMwEqt9fA77CJbWFiYvn1Bh4MHD1K3bt18RC9+kpOT8fb25tKlSzRp0oStW7dSqVIlo2OJm8j3c9Gy7sA5Xvk+imaufzPNZwFel/ZBtZbQ+SOoFFxoOZRSkVrrsNsfz80Z+BzgM+Cbm3Z2P/AIUF9rfV0pVcFWQcWdde3alYSEBNLT03n33XelvIUoIFprvth4jDmrt/NZySXcn/47pPtBr1kQ1AvsZNqEuxa41nqTUirwtoeHAuO11teztjlv+2jidnca9xZC2E5ahpm3lkTiu282mzyW4242Q+vXM3+53fliBCPkdwz8XqC1Uup9IA0YprXemdOGSqkhwBCAgICAfB5OCCEKxs0zBlYs5UFrFcULqTOp6XoGXasTqtOHULaG0TFzlN8CdwHKAM2AxsAipVQNncOAutZ6BjADMsfA8xtUCCFs7eYZAwPUOd5N/ZYHnHdzybMq9F6CqvWA0RH/VX4LPBZYmlXYO5RSFsAXkMmQhRAOY+Lqw6RmmBjs/CtvunxPBs58mNGfX917sNHOyxvyX+DLgfbABqXUvYAbcNFmqYQQohAkJ1xghuuXPOgcyVpzI97OeIrzlEElOsZNaXe9kUcptQD4A6itlIpVSj0NzAZqKKX2Ad8DA3MaPnEUMp3sPxXEdLKjR49m0qRJNt9vfowcOZJ169b96zYbNmyQWRmLsGvH/uAX97do5xTN2IwneTbjNc5TBrDtjIEFKTdXofS/w1NP2DiL4WQ62f8p6sU1duzYu26zYcMGvL29ZW70osZi4fK6jym17UPMlKWfaQy7zf97k9LWMwYWJLmV/iYynWzhTSdri4w38/b25vXXXyc0NJQOHTpkr00ZHR1Ns2bNqF+/Pj169ODKlSsADBo06JbXb9SoUYSGhhIcHMyhQ4eIiYlh+vTpTJ48mYYNG7J58+YcswkHc+0Sl2b2pOy2cWwgjDP91jCgd0/8fTxRgL+PJx/2DC7QCahsyb5upV/1Jpzda9t9VgqGzuPz9aEynWzBTSc7ZMgQm2UEuHbtGqGhoXz88ceMHTuWMWPG8NlnnzFgwACmTZtG27ZtGTlyJGPGjGHKlCn/+HhfX192797Nf//7XyZNmsTMmTN57rnn8Pb2ZtiwYQAEBwf/I5twHPrkNlLmD8Q77TKfez1Ht2dGUrVc5nXdjlLYt7OvArczMp1swUwna+uMAE5OTtmvwRNPPEHPnj1JTEwkISGBtm3bAjBw4MBbjnmzG8dq1KgRS5cuzXGbnLIJB2CxkLHpE5w2fMAFiy/fBUzj1ScfpYS749effX0G+TxTLigynWzBTCdr64w5yesK8Tc+l3/7PHLKVq5cuTwdRxSy5AtcX/wM7ic3sMLcjNMtP+StB0NwcrKPW+GtJWPgeSTTyd59OtkRI0awbNmyOz5v64yQ+Z/Cjdd8/vz5tGrVitKlS1OmTBk2b94MwLfffpt9Np4bJUuWJCkpKfvvOWUTduzEZtI/bwEntzLK8iwufb7m+U6hRaa8wd7OwB2ATCd79+lk9+7dS7du3f7187RlRsj8CWT//v00atSI0qVLZ7/ZOnfuXJ577jlSUlKoUaMGX3/99d2+BNkefvhhevfuzY8//si0adOYPHnyP7IJO2QxY9k4ETZ+xGlLRcZ7T2L4wF7UqljS6GQ2d9fpZG1JppO1jqNMJxseHs7q1asL9Zje3t7ZVwMZSb6fDZZ0DtPip3E5tZml5lZsrfMWY3s3dfjxbmumkxV2wlGmky3s8hYCgGO/k7HkWcypibxj+g91Ow9lUovAPL8f4kikwB2ITCd7Z/Zw9i0MYjZh2fAhavPHxFj8GO0+gWFPdSckoIzRyQqcFLgQwnFdjef6wsG4x21nkakt22oP5/NeTfDxcjM6WaGQAhdCOKaj67i++BnM6Sm8YXmBRt2HMrlRlSI9ZHI7KXAhhGMxm0hbOxaP7VM5bqnKtLIfMvyJbgT62tdqOYVBClwI4TgSY0n49kl8Lu5mgbkD55qPYmp4MK7OxfOWluL5Wf+Lu015unz5cg4cOFCIiYQQAMl//cy1qc1xuXCADzwjqPef2fxflwbFtrxBzsDzbPny5XTt2pX77rvP6ChCFGk31qq8kJDEOx6LGcAKDliqsT3sE4Z1aY+bS/Et7hsc7hVYHhVHy/Hrqf7mSlqOX8/yqDir9/n+++9Tu3ZtOnbsyOHDhwH46quvaNy4MQ0aNKBXr16kpKSwbds2fvrpJyIiImjYsCHHjh3LcTshhHVurFWpEk+x0G0sA1jBN6YH2NhmPk916yjlncWhXoUbX9S4hFQ0EJeQyoile60q8cjISL7//nuioqJYunQpO3fuBDJnp9u5cyd79uyhbt26zJo1ixYtWtCtWzcmTpxIdHQ0NWvWzHE7IYR1Pl61j8ctP7HK7U1qqjiGpr/CSNNgvtt13uhodsWhhlAyFyC9dS7o1AwzE1cfzvd8vps3b6ZHjx54eXkBZM/hsW/fPt555x0SEhJITk7+x+yAN+R2OyFE7hzcuoLZaW9QyzWO38whjDINJFZXACA+IdXgdPYlN2tizlZKnc9a//L254YppbRSyrdg4t3qTl88a7+oOV03OmjQID777DP27t3LqFGjbpkaNj/bCSH+3aX4E0R/0p26a5/ATWXwVPowns6IyC5vcJy1KgtLboZQ5gCdbn9QKVUVeAA4ZeNMd3SnL541X9Q2bdqwbNkyUlNTSUpKYsWKFQAkJSVRuXJlMjIymDdvXvb2t08xeqfthBC5cz0the3fvIPnl02pk7iFzVWGsKPzL/zhfOuMno60VmVhyc2ixpuUUoE5PDUZGA78aONMdxQRXpsRS/feMoxi7Rc1NDSUvn370rBhQ6pVq0br1q0BeO+992jatCnVqlUjODg4u7T79evHs88+y6effsqSJUvuuJ0Q4t9prYn8bQkVto6kmY5nd4mWlOv1Ma1rZs7m6OruxcTVh4lPSMXPx5OI8NoOu/RZQcnVdLJZBf6z1joo6+/dgA5a61eUUjFAmNb64h0+dggwBCAgIKDRyZMnb3k+r9Nv3ri0SL6owh7JdLK5c/Twfq4si6BJ2lZinfy40uY9gtv1NjqW3bLZdLJKKS/gbeDB3GyvtZ4BzIDM+cDzerzbdQ/xl8IWwkFdTEgkasEYWp/9hirKid21XqZ+77eo4i5j2/mRn6tQagLVgT1Zb/5VAXYrpZporc/aMpwQomi4bjKz/sdvCPrrQx5Q59hftgNV+35CaKVAo6M5tDwXuNZ6L5D9tvDdhlByuc9iNYOYKJoKc3UrR6G1ZvOOnbisGUFn8y7i3QKI7/w99UI7Gx2tSLhrgSulFgDtAF+lVCwwSmtts7tVPDw8uHTpEuXKlZMSFw5La82lS5fw8PAwOord2H/yLIcWj6Fr0mLMyoVjIW9S86HXwaV4zNVdGHJzFUr/uzwfaE2AKlWqEBsby4ULF6zZjRCG8/DwoEqVKkbHMNz5xFRWLZlJh1NT6KUucsKvC1X7TqKmj7x3ZWuG34np6upK9erVjY4hhLBSWoaZxb/+TuCusQxUezjnVZPkHrOpfm9bo6MVWYYXuBDCsWmtWRn5N5dXfUA/04+YnT241HIsFdu9AM5SMQVJXl0hRL5Fxlxm/Q9f8vjVGfipy1y4pxfle4zH07vC3T9YWE0KXAiRZ7FXUpiz/FfuPz6JCOf9XCldB0uvBZSv1szoaMWKFLgQItfSMszMWf8Xblsn8ob6FbObF9c7TKBMs2fAydnoeMWOFLgQ4q601vx24Bx//DidIde/prxTIqlBj1Gi81goUSiTkYocSIELIf7ViYvXmP3DCrrGTeZdp0MklwvGqedSSlRpZHS0Yk8KXAiRo5R0EzPXRlF6+yRGOa3B5F4K04NT8G40QIZL7IQUuBDiH1bvi2fn8s/5T8Y3lHNOIq3BQLzCR4FXWaOjiZtIgQshsp1JTGXWomV0Pv0J7zgdJblCCE49puDl19DoaCIHUuBCCMwWzcJNe3D6/X3eYi3X3X0wdf4c74aPgZNDrX1erEiBC1HMHYhLYP2Cj3ksaTY+KoVrDZ6iZKeR4OljdDRxF1LgQhRTqelmFi5fRsi+D3jR6RgXfRuh+kylZKVgo6OJXJICF6KIy2kZwgCPFOJ/GMGAjHUku5XhWvgX+Ib1B5nS2aFIgQtRhC2PirtlIfAzCdfY/cNE2jkvor5K42y9p/HrNgo8ShmcVOSHFLgQRdjE1YezyztUHeE916+p53SSP3U96j87Az//IIMTCmtIgQtRhMUnpOJOOm+5zGOgy1rO6LK8kP4yv1iackLK2+HlZkm12UBX4LzWOijrsYnAw0A6cAwYrLVOKMigQoi8a1LiDGMzJlPbKZbZpk5MMj1KCh74+8gq8EVBbi7wnAN0uu2xtUCQ1ro+cAQYYeNcQggrmExm1s8dyzemNymrkhiQ/gZjTQNIwQNPV2ciwmsbHVHYQG7WxNyklAq87bE1N/11O9DbtrGEEPkVF3uKs988Rfv0nRws2YzDTcdzbOsV1E1XoXQPkfUpiwJbjIE/BSy805NKqSHAEICAgAAbHE4IcSfb1yzinm3DCNIp7G3wNsE9IqirFN1bG51MFASrClwp9TZgAubdaRut9QxgBkBYWJi25nhCiJylpFxj16z/o82lRZxyDsDUdynB94YZHUsUsHwXuFJqIJlvbnbQWksxC2GQo/t2wdJnaGM5QVTFPgQNnoqrRwmjY4lCkK8CV0p1At4A2mqtU2wbSQiRG9piYeuij2l0cALXlTsH231JSLt+RscShSg3lxEuANoBvkqpWGAUmVeduANrVeatt9u11s8VYE4hxE0uXTjDidlP0yp1Kwe8QvEbNJe6FeU9puImN1eh9M/h4VkFkEUIkQt7Nv9Epd9eob5OJLL2a4T2ewclK+QUS3InphAOIv36dXbNGUaz+G+Jc/Yjrsd3NApuaXQsYSApcCEcwOm/95H6/WBamI6wq9zDBD39XzxKyARUxZ0UuBB2TFss7Prpv9wX9R5m5Ux086mEhQ8yOpawE1LgQtippISLHJn1LI2T1nPAPZiyT86hYdV7jI4l7IgUuBB26MjOtZT8ZSgNLJf4I3AoTZ4ch7OL/HMVt5LvCCHsiMWUQeS3bxEa8xVnnSrwd9clNG/cwehYwk5JgQthJy6cPsLlbwfSOP0Af5Z6gDpPz8Dfp6zRsYQdkwIXwiA3r1XZ13MHb1u+xB/YHjKepo88h5L1KcVdSIELYYAba1U6ZSQz0XUuvfUmdut7ONRiMo91amN0POEgpMCFMMDE1YepazrIJ25fUFWdZ6qpB5+aelIp2sRjty+fIsQdSIELUci06Tr9k+cw1O0nzlCOfunvslPXATLXsBQit6TAhShESaf2cWXeQF50+ZvFpjaMMQ0gGa/s5/1krUqRB1LgQhQGi4WTqz6h0s7xeGsPZvq/x6RTtUjDkr2JrFUp8koKXIgCZrp8kvi5T1EtcRfbnMMo3Xc6z9xbC9+brkKRtSpFfkiBC1FQtObSH9/isfZNyllMLPSL4KGBb+Dt4QpA9xB/KWxhFSlwIQpCymXiv3sOv/jV7Na1ufjAVPq2am50KlHESIELYWOpB34lfenz+GYk8K33INoNHkeob0mjY4kiSApcCFtJv8aFpRGUPzSPU5Yq7AiZSv9uD+Hi7GR0MlFE3fU7Syk1Wyl1Xim176bHyiql1iqljmb9XqZgYwph38yndpAwuRnlDs5nvvMjJA1Yx5M9HpbyFgUqN99dc4Db7w17E/hNa10L+C3r70IUP+YMrq4ajZodTnJKCtMCPuGhYbMIu6ey0clEMZCbRY03KaUCb3v4ETJXqgeYC2wA3rBhLiHs34XDJMwbjE/CfpbrtqjOH/Fy0zoyCZUoNPkdA6+otT4DoLU+o5SqcKcNlVJDgCEAAQEB+TycEHbEYiF16xc4rx+NxeLGh6Xf5rGBL1CtXAmjk4lipsDfxNRazwBmAISFhemCPp4QBSoxlivzn6XMuW2st4Twd7MPGBbeDFcZ6xYGyG+Bn1NKVc46+64MnLdlKCHsjtZcj16IZcXruJkzmOz5PB0eH077qvL+vTBOfgv8J2AgMD7r9x9tlkgIe5NymcuLXqRszEoiLbXYXv8Dhj7SAQ9XZ6OTiWLurgWulFpA5huWvkqpWGAUmcW9SCn1NHAK6FOQIYUwSvrhNVxf8jwl0y/zpevjNOg3ihfuqWh0LCGA3F2F0v8OT8lKq6LoSk/h0vI3KXdgLjEWf1bVmsFTfR6hZNY8JkLYA7kTU4jbXI/ZwbXvn6Zc2inmOXWlcq8PeSVYrqAS9kcKXIgbzBnEr3iPitHTSNVl+KL6FB579AlKe8lZt7BPUuBCAKnxB7j83WD8Uw6xyqkdJXt+wtCgmkbHEuJfSYGLYmf5TQsp+Jd2Z1TlP2gT8yle2o3vq4+ja/+heLvLPw1h/+S7VBQry6PiGLF0L6kZZipymQ9SvqRNzF62qhA8e39Bv6C6RkcUItekwEWxMnH1YVIzTDzs9AfjXGfjipm3Mp5mQ4kubJPyFg5GClwUK86JMcx2nUt752h2W+7htYyhxOjKqKvXjY4mRJ5JgYtiISnpKlHzR7HW7VvSceG9jCeYYw7HTObdlH4+ngYnFCLvpMBFkaa1Zsfq+VTZPpo2nGerVztGJD3KKbNP9jaers5EhNc2MKUQ+SMFLoqsmL/3c3HJazRN284p56r8/eACWjbtwms3XYXi5+NJRHhtWR1eOCQpcFHkpKQks2v+GJqc/poKyonI2q/SsPcInF3dAege4i+FLYoEKXBRZGitifxtMRW3jqSNPsOe0vdTtf9kGlWubnQ0IQqEFLgoEk6dOMy5Ra/ROHULsU7+HOk4lwYtuxsdS4gCJQUuHFpaWio7579Ho5MzKY9m1z0v0fDRt3Fxl6tKRNEnBS4cVuTvy/Dd9DatdRx7SrbGv99kwqrUMjqWEIVGClw4nPhTx4hb+CqNr20kTlXiYPvZNGjTy+hYQhQ6KXDhMNLS0ti58ANCj39JWczsqvEcDfqNwt/dy+hoQhhCClw4hOjNKyi9fgSt9Wn2ejen4qNTCKtWx+hYQhjKqgJXSr0KPANoYC8wWGudZotgQgCcjYvh5ILXaZq8jjOqAvvbfknw/f2MjiWEXch3gSul/IGXgfu01qlKqUVAP2COjbKJYiw9PZ0diz6iwdHPCSGDndWeoX7/0VT2LGl0NCHshrVDKC6Ap1IqA/AC4q2PJIq7vdt+pcS64bSynGSfV2PK9ZlC4xpBRscSwu7ku8C11nFKqUnAKSAVWKO1XnP7dkqpIcAQgIAAWRhW3NmFM6c4vmAYTa+u5pzy5a9Wn1G/wxOglNHRhLBLTvn9QKVUGeARoDrgB5RQSj1x+3Za6xla6zCtdVj58uXzn1QUWaaMDLbO/xD36U0JSVzHjiqDKD0sivodn5TyFuJfWDOE0hE4obW+AKCUWgq0AL6zRTBRPBza+Rsuvw6jpfk4+z1D8Ok1hSa1GhodSwiHYE2BnwKaKaW8yBxC6QDsskkqUeQsv20K11dblCUgagJNrqzkAmWJbjqZBuGDUE75/qFQiGLHmjHwP5VSS4DdgAmIAmbYKpgoOm5eSNgJC22Tfqbjb99TgjS2V36c4Mfep2GpMkbHFMLhWHUVitZ6FDDKRllEEZW5kLCZuuok412/ooHTcbZb6jLFbQjfP/eU0fGEcFhyJ6YocGcTknnBeQWvuPxAIt68nP4CP1laoNLlDUohrCEFLgrUlj+2stRtFA2cjrPC3Ix3MwaTQObNOLKQsBDWkQIXBeJ8YgobvxnLwxdnct3Jnf8zvcxyU7Ps52UhYSGsJwUubEprzc8bt1F5w+v04SAxvq3xH/AV7Y5b2CkLCQthU1LgwmZiLiSz9rvxPJYwA+XkzIX2kwlsNRiUonsIUthC2JgUuLCayWxh/trtVN/2Bs86/cVZ32ZUeOIrvMrI1AlCFCQpcGGVfbEJ/LpgKkOSp+PubOHq/eOp1Oo/IDfkCFHgpMBFvpjMFuas2UHAH28zzGkXV3xDcX9sJu7lahodTYhiQwpc5NnJS9dYNPcznkqcRinnNFLbjqZMm5fBydnoaEIUK1LgIte01izdug/3tW8QobaSUKYero/NxrWCLG0mhBGkwEWuXEi6zvxvZ9Dv3CTKqSSuNh+OT8fh4OxqdDQhii0pcHFX6/cc4+ryCF7Rv3HZ+x6cHvuRUv4y5asQRpMCF3eUfN3E/AXf0uXEOCqry1wOeYGyD40CF3ejowkhkAIXd7D77ziOf8ENJv8AABAeSURBVB/BENNKLntWxdJvFWUDm939A4UQhUYKXNwi3WRh8fIltPjrXXo7neVc3UFU7PEhuHkZHU0IcRspcJHt7/iLRH3zBv1SfyDRvSIpvZdRsXZ7o2MJIe5AClxgsWh+Xr2K2tuH00edJrZGH6r0/QQ8ShkdTQjxL6wqcKWUDzATCAI08JTW+g9bBBOF48zlq2yd8zaPJM7jmosPid3mU6XBQ0bHEkLkgrVn4FOBX7XWvZVSboAMlDqQ3zdvpMJv/0dvjhPj14VqT36O8iprdCwhRC7lu8CVUqWANsAgAK11OpBum1iiICUmp/H73FF0Pj+LNCdPzod/RWCzR42OJYTII2vOwGsAF4CvlVINgEjgFa31tZs3UkoNAYYABATI9KJGi4zahctPL9BdH+KYbzuqDZxB6VIVjY4lhMgHa+b8dAFCgS+01iHANeDN2zfSWs/QWodprcPKly9vxeGENdLSM/hl1ljqLu9CTX2aU20nU/PF5bhIeQvhsKw5A48FYrXWf2b9fQk5FLgw3uHDB7i26Dm6mPdwtFRTqgycSYCv/DQkhKPLd4Frrc8qpU4rpWprrQ8DHYADtosmrGU2W9i4aAqND03AWWmONBnHvZ1fBKWMjiaEsAFrr0J5CZiXdQXKcWCw9ZGELcSeOsHZ74bQPn0HR73qU+GJWdzrf6/RsYQQNmRVgWuto4EwG2UR+bQ8Ko6JWSu+Vy7twYsV/qLLqUkEqevsDXqToJ7DUbLYghBFjtyJ6eCWR8UxYuleUjPMlOEqb6V8StfT2znofC9ln5hFcI36RkcUQhQQKXAHN3H1YVIzzHR0iuRD15mUJpkJGY+yokQfNkt5C1GkSYE7ONfE43zlOo8HnHdzwFKNARlvclBXQ13NMDqaEKKASYE7qrSrnFo+hjVuX5OOKx9m9Ge2uTMZWV9SPx9PgwMKIQqaFLijsVhI2fkN5rVjCDBdZoVze8anP0qc+X8zB3q6OhMRXtvAkEKIwiAF7khO/UnislcpfWU/kZZa7A2aRL8e3THvO5t9FYqfjycR4bXpHuJvdFohRAGTAncEiXGkrnoHz0NLSdFl+dJ7GJ37vcSgqj4AdA/xl8IWohiSArdnGalYtn6KZfNknMwm/mvpgWvbYbza7j5cna2ZxkYIURRIgdsjreHActJ/eRu3a3GsMjdhVaXneb3vg1T3LWF0OiGEnZACtzdn/iJj5XBcY//gb0s1prqOpdMjvZjW0B8lc5gIIW4iBW4vrl3EtG4szlHfkIw3H5uexqvZYCZ1rENJD1ej0wkh7JAUuNFM6Vh2zMC0fjxOphRmmzqxo9ozRDzSlHsqlDQ6nRDCjkmBG8hyZC0pKyLwTjrBH+b6LCw3lMe7PsjT9/gaHU0I4QCkwAvJzTMGNil5mTEe86mTtI3zlkq87/UurTo/xuf1K8s4txAi16TAC8GNGQNdMpIY4bKMwem/kpruxqcuA6n28Ku817AaLnJZoBAij6TAC8EHvxykjfkP3nOfgy+JLDS342PTo7h5VGRbo+pGxxNCOCgp8AJyMfk6v+w9w4Zdexmd9ild3HawzxLIUxnD2KdrAKAS0wxOKYRwZFLgNnT+ahq/Hz7PL3vPsuXvC/RQG5nq9h1uzumMz+jHV+aHMPO/lXFkxkAhhDWsLnCllDOwC4jTWne1PpLjSMswsy8ukY1HLvD74fPsi7sKQKNSSawrP4vqiTuganPW3fMOc9cmYzabsz9WZgwUQljLFmfgrwAHgVJ329DRnb+aRuTJK+w+dYXIk1fYF3eVdLMFJwWhAWWIeLAWvcyrqLjjI1Sqgi6TIOxpOjo58WHJOJkxUAhhU1YVuFKqCvAQ8D7wmk0S2QmT2cKhs0nZZR158gqxV1IBcHNxor5/aQa3DCQkoAxNq5elTMoJ+OlFOP0n3NMRuk4Bn6rZ+5MZA4UQtmbtGfgUYDhwx1sGlVJDgCEAAQEBVh6u4KSkm4g6lcDOmMvsisk8y05JzxzyqFDSnUbVyjCoRSCh1cpQz68U7i5ZY9nmDNg6BTZOALcS0ONLqN8X5HpuIUQBy3eBK6W6Aue11pFKqXZ32k5rPQOYARAWFqbze7yCcPpyCusOnmP9ofP8efwy6WYLSkGdSqXo3agKjaqVoVG1Mvj7eOZ8g018NPz4IpzbC/V6QOcJ4F2h8D8RIUSxZM0ZeEugm1KqC+ABlFJKfae1fsI20QpGarqZFXviWRIZy46YywDULF+CgS2q0eIeXxpVK0Opu00elZEKG8bDtmlQojz0nQd1i9X7t0IIO5DvAtdajwBGAGSdgQ+z5/K+mpbB3K0xfL0thsvX0qnhW4LhnWrzUHBlqpXLwxzbMVvhp5fg8jEIeRIeHAeePgUXXAgh7qDIXwdusWgWR55mwq+HuXQtnfZ1KvCfNjVoUr1s3uYdSbsKv42BnTPBpxoM+BFqtCuo2EIIcVc2KXCt9QZggy32ZUuxV1J4dWE0O2OuEFatDHMGNyG4Sum87+jIGvj5VbgaB81egPZvZ75hKYQQBiqyZ+BrD5zjtUXRaA2T+jSgV2g+VrS5dglWj4C/FkL5OvD0WqjauGACCyFEHhXJAp+95QTvrTxAsH9pPn8slKplvfK2A61h/zL4JQLSEqDtG9D6dXBxL5jAQgiRD0WuwKesO8KUdUcJr1eRKX1D8HRzvvsH3ezqGVj5OhxeCX4h0O1HqBRUMGGFEMIKRarAp288xpR1R+kVWoUJvevj7JSHIROtYfc3sOZdMF/PvLqk6VBwLlIvkRCiCCky7fRjdBzjVx3i4QZ+eS/vyydgxctwYhNUawXdPoVyNQsurBBC2ECRKPB9cYm88cNfNAksy8d9GuS+vC1m+HM6/PYeOLtmzl8SOhCcZHUcIYT9c/gCv5qWwX++jaSMlxv/fSIUN5dclu+5A5k35MTtgns7wUOfQGmZbEoI4TgcvsDH/XyAM4mpLBnaAl/vXFwlYkqHLZ/ApkngUQp6zYKgXjL5lBDC4Th0ga8/dI5Fu2J5vl1NQgPK3P0D4iIzJ586fwCC+0Cn8VDCt+CDCiFEAXDYAk9NN/POsn3UrliSVzrW+veN01Pg9/dh+3/BuxL0Xwi1OxVOUCGEKCAOW+BfbjpGfGIai/qF/G9u7pyc2JQ51n0lBhoNhgfGgEc+bqcXQgg745AFHp+QyvSNx3gouDJNqpfNeaO0RFg7EiLnQJnqMPBnqN66UHMKIURBcsgC/2TtESwa3uxcJ+cNDq/KnHwq+Ry0eBnajQC3PN5OL4QQds7hCvzUpRSWRcUxoHm1f85xcu0irBoO+36ACvWg33zwDzUmqBBCFDCHK/AvNh7DWSn+0+amOyW1hr1LMsv7ehLc/za0/D9wcTMuqBBCFDCHKvD4hFSWRJ6mb+OqVCrtkflgYiz8/BocXQ1VGkO3aVChrrFBhRCiEDhUgc/ZFoNFw3Nta4LFApFfw9pRoM2Z13Q3GQJOeZx9UAghHJTDFHhahplFu04TXq8iVSxnYO7LcHILVG8LD0+FstWNjiiEEIUq3wWulKoKfANUAizADK31VFsFu9nyqDhG/7SfpNQ0av/9NebPF+Hs6gHdPoOQJ+Q2eCFEsWTNGbgJeF1rvVspVRKIVEqt1VofsFE2ILO8RyzdS6DpON+4zaC+PsE6S2NMHSfSKTTElocSQgiHku8C11qfAc5k/TlJKXUQ8AdsWuATVx8mzBzFbLeJJFCC59Nf5hdLU/w3XaFTC1seSQghHItNxsCVUoFACPBnDs8NAYYABAQE5Hnf8QmpXKQOs82d+MLUjQRKZj8uhBDFmdUrFyilvIEfgP/TWl+9/Xmt9QytdZjWOqx8+fJ53r+fjyfXceND0+PZ5X3jcSGEKM6sKnCllCuZ5T1Pa73UNpFuFRFeG0/XWy8N9HR1JiK8dkEcTgghHIY1V6EoYBZwUGv9ie0i3ap7SOYqORNXHyY+IRU/H08iwmtnPy6EEMWVNWPgLYEngb1Kqeisx97SWv9ifaxbdQ/xl8IWQojbWHMVyhZALsAWQgiDyPLrQgjhoKTAhRDCQUmBCyGEg5ICF0IIB6W01oV3MKUuACet2IUvcNFGcWzFHjOB5MoryZU3kitvrM1VTWv9jzshC7XAraWU2qW1DjM6x83sMRNIrrySXHkjufKmoHLJEIoQQjgoKXAhhHBQjlbgM4wOkAN7zASSK68kV95IrrwpkFwONQYuhBDifxztDFwIIUQWKXAhhHBQdl/gSqnZSqnzSql9Rme5mVKqqlLqd6XUQaXUfqXUK0ZnAlBKeSildiil9mTlGmN0ppsppZyVUlFKqZ+NznKDUipGKbVXKRWtlNpldJ4blFI+SqklSqlDWd9nze0gU+2s1+nGr6tKqf8zOheAUurVrO/5fUqpBUopD6MzASilXsnKtN/Wr5Xdj4ErpdoAycA3Wusgo/PcoJSqDFS+eVFnoLutF3XORy4FlNBaJ2ctuLEFeEVrvd3IXDcopV4DwoBSWuuuRueBzAIHwrTWdnUDiFJqLrBZaz1TKeUGeGmtE4zOdYNSyhmIA5pqra25Qc8WWfzJ/F6/T2udqpRaBPyitZ5jcK4g4HugCZAO/AoM1VoftcX+7f4MXGu9CbhsdI7baa3PaK13Z/05CbixqLOhdKbkrL+6Zv2yi/+llVJVgIeAmUZnsXdKqVJAGzIXTUFrnW5P5Z2lA3DM6PK+iQvgqZRyAbyAeIPzANQFtmutU7TWJmAj0MNWO7f7AncE/7aosxGyhimigfPAWq21XeQCpgDDAYvRQW6jgTVKqcisRbjtQQ3gAvB11pDTTKVUCaND3aYfsMDoEABa6zhgEnAKOAMkaq3XGJsKgH1AG6VUOaWUF9AFqGqrnUuBW+luizobQWtt1lo3BKoATbJ+jDOUUqorcF5rHWl0lhy01FqHAp2BF7KG7YzmAoQCX2itQ4BrwJvGRvqfrCGdbsBio7MAKKXKAI8A1QE/oIRS6gljU4HW+iDwEbCWzOGTPYDJVvuXArdCYSzqbI2sH7k3AJ0MjgKZS/B1yxpv/h5or5T6zthImbTW8Vm/nweWkTleabRYIPamn56WkFno9qIzsFtrfc7oIFk6Aie01he01hnAUqCFwZkA0FrP0lqHaq3bkDkcbJPxb5ACz7fCWtQ5r5RS5ZVSPll/9iTzG/uQsalAaz1Ca11Fax1I5o/e67XWhp8hKaVKZL0JTdYQxYNk/thrKK31WeC0Uqp21kMdAEPfIL9Nf+xk+CTLKaCZUsor699mBzLflzKcUqpC1u8BQE9s+LpZs6hxoVBKLQDaAb5KqVhglNZ6lrGpgEJc1DmPKgNzs64QcAIWaa3t5pI9O1QRWJb5bx4XYL7W+ldjI2V7CZiXNVxxHBhscB4AssZyHwD+Y3SWG7TWfyqllgC7yRyiiMJ+bqv/QSlVDsgAXtBaX7HVju3+MkIhhBA5kyEUIYRwUFLgQgjhoKTAhRDCQUmBCyGEg5ICF0IIByUFLoQQDkoKXAghHNT/A8VXebpwbf3DAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "poly9 = numpy.polyfit(x, y, 8)\n", "regplot(poly9)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we had just plotted the connecting lines, we would have missed the bit sticking up on the left!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Nonlinear regression\n", "\n", "We can apply the same principles to fit nonlinear functions as well. The `scipy.optimize.curve_fit` function can be used to fit an aribitrary function to data" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "import scipy.optimize" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's start by reproducing the results from the linear fit" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "def f(x, a, b):\n", " \"\"\"fitting function linear in coefficient\"\"\"\n", " return a*x + b" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1.85932396, 0.07854919])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beta, _ = scipy.optimize.curve_fit(f, x, y, [1, 0])\n", "beta" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's build some data which is obviously nonlinear" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "x = numpy.arange(1, 10)\n", "y = 2*numpy.sin(3*x) + x + 1" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "def f2(x, a, b, c, d):\n", " \"\"\" A nonlinear fitting function\"\"\"\n", " return a*numpy.sin(b*x)+ c*x + d" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "def fit_and_plot(beta0):\n", " beta, _ = scipy.optimize.curve_fit(f2, x, y, beta0)\n", "\n", " plt.scatter(x, y)\n", " plt.plot(smoothx, f2(smoothx, *beta))\n", " return beta" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.32658273, 1.45786791, 1.0854963 , 0.6763682 ])" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXyU1aH/8c8hCSRhC0uALISAQNgCBCOgqFXRIuCCom29SHvrQmvdWpVWaluv92r7s3T72fbS0qq1LlhEQKxIwA0ULVtCTAKELRCyQUJISMJknXP/SLBIsSSZmTwzk+/79eJF8jDO833F5JtnzpznHGOtRUREAk8XpwOIiEj7qMBFRAKUClxEJECpwEVEApQKXEQkQIV25Mn69+9vExMTO/KUIiIBb8eOHWXW2uizj3dogScmJrJ9+/aOPKWISMAzxhw+13ENoYiIBCgVuIhIgFKBi4gEKBW4iEiAOm+BG2OeM8YcM8Zkn3FssTFmjzHmU2PMKmNMlG9jiojI2VpzBf4X4Nqzjm0AxllrxwN7gUVeziUiIudx3mmE1tpNxpjEs46tP+PTfwC3eDeWiEjgW51RyOK0XIoqXMRGRbBwRhJzUuK89vzemAd+B/C3L/pHY8wCYAFAQkKCF04nIuL/VmcUsmhlFq6GJgAKK1wsWpkF4LUS9+hNTGPMY0Aj8PIXPcZau9Ram2qtTY2O/pcbiUREgtLitNzPyvs0V0MTi9NyvXaOdl+BG2O+AVwHTLfaFUJE5HOKKlxtOt4e7boCN8ZcC/wAuMFae8praUREgkRsVESbjrdHa6YRLgM+AZKMMQXGmDuB3wE9gQ3GmJ3GmD94LZGISBBYOCOJiLCQzx2LCAth4Ywkr52jNbNQbjvH4We9lkBEJAidfqPS32ehiIjIOcxJifNqYZ9Nt9KLiAQoFbiISIBSgYuIBCgVuIhIgFKBi4gEKBW4iEiAUoGLiAQoFbiISIBSgYuIBCgVuIhIgFKBi4gEKBW4iEiAUoGLiAQoFbiISIBSgYuIBCgVuIhIgFKBi4gEKBW4iEiAUoGLiAQoFbiIiA9V1zWy7VC5T55bBS4i4iN5ZTXc9PvN3PH8NipdDV5/fu1KLyLiA5v2lnLfK+l06WJYcvuF9I4I8/o5VOAiIl5kreXZj/L46drdjBjQkz99PZWEfpE+OZcKXETES2obmnhsVTavpxfw5TED+dVXJ9Kjm+9qVgUuIuIFx07WsuDFHew8UsED00fw3ekj6NLF+PScKnAREQ9lHqlgwYvbOelqZMm8ScxMjumQ86rARUQ8sCqjgB+8nkV0j268fs8ljInt1WHnVoGLiLRDk9vy9Lo9LN10kClD+/K/8ybRr0e3Ds2gAhcRaaNKVwMPLMtg495S5k8dwk+uH0NYSMffVqMCFxFpgwOl1dz9wnbyy0/x1E3jmDdliGNZzvsrwxjznDHmmDEm+4xjfY0xG4wx+1r+7uPbmCIizns/9xhzfreZClcDL981xdHyhtbdSv8X4Nqzjj0KvGutHQG82/K5iEhQstbyx40HuOMv24jvG8ma+6YxZVg/p2OdfwjFWrvJGJN41uEbgStaPn4B+AD4gRdziYj4hdqGJh59/VNW7yxidnIMi28dT2RX/xh9bm+KgdbaYgBrbbExZsAXPdAYswBYAJCQkNDO04mIdLySyloWvLidTwsqefiakdx31XCM8e3NOW3h818j1tqlwFKA1NRU6+vziYh4w47DJ/j2Szs4VdfI0vkX8uWxg5yO9C/aW+BHjTExLVffMcAxb4YSEXHS8u1H+NGqbAb1DuelO6eQNKin05HOqb0TF9cA32j5+BvAG96JIyLinMYmN0+8mcP3V3zK5KF9WXPfNL8tb2jFFbgxZhnNb1j2N8YUAI8D/w9Yboy5E8gHbvVlSBHxD6szClmclktRhYvYqAgWzkhiTkqc07G84kRNPfctS2fz/uPcMW0oP5w1ilAHbs5pi9bMQrntC/5pupeziIgfW51RyKKVWbgamgAorHCxaGUWQMCXeG5JFXf/dTsllbX8/JbxfCV1sNORWsW/f72IiN9YnJb7WXmf5mpoYnFarkOJvCMtp4Sb/3czroYmXv3W1IApb9Ct9CLSSkUVrjYd93dut+W37+3n1+/sZUJ8b/44P5VBvcOdjtUmKnARaZXYqAgKz1HWsVERDqTxTE1dI4+8lsnb2SXcnBLHT29OJjwsxOlYbaYhFBFplYUzkog4q+QiwkJYOCPJoUTtc6C0mjm/30xaTgmPzRrNL78yISDLG3QFLiKtdPqNykCehbIuu5hHXvuUrqFdePHOKUwb3t/pSB5RgYtIq81JiQuowj6tscnN4vW5/HHjQSYMjmLJvEkBOfRzNhW4iAS1suo67n8lg08OHuf2qQn8+LoxdAsNzCGTs6nARSRopeef4DsvpXPiVD2/uHUCt1wY73Qkr1KBi0jQsdby/OZD/Ozt3QzqHc7K71zC2NjeTsfyOhW4iASVEzX1LFzxKe/sPsrVowfwy1sn0jsyzOlYPqECF5Ggse1QOQ8sy6Csuo6fXDeGb05L9Kv1u71NBS4iAa/JbVnywX5+/c4+4vtEsPKeaSTHB9+QydlU4CIS0I6drOWh5Zl8tL+MGybE8tRN4+gZHpxDJmdTgYtIwFqbVcxjq5pXSHx6bjJfSR0c1EMmZ1OBi0jAqXQ18F9rcliVUciE+N786qsTuSC6h9OxOpwKXEQCysf7y3jktUyOVtXx4PQR3HfVcML8fOMFX1GBi0hAOFXfyOK0XJ7ffIih/bvz+j2XMHFwlNOxHKUCFxG/99G+Mhat+pQj5S7mTx3ColmjiOyq+tJXQET8VuWpBp58axev7ShgaP/uvLpgKlOH9XM6lt9QgYuI37HW8nZ2CY+vyaG8pp57rriAB6ePCNh1u31FBS4ifuVgaTWPr8nhw31ljIvrxfP/eRHj4oL/ppz2UIGLiF+oqWvkd+/v588fHiQ8NITHrx/D/KlDCO2kM0xaQwUuIo6y1vJWVjFPvbWb4spa5k6K59GZo4ju2c3paH5PBS4ijtl+qJyfrt1Nen4FY2J68bv/SOHCIX2djhUwVOAi0uEOlFbz83V7SMs5yoCe3fjZzc23wYd06Ty3wXuDClxEOszRk7X89r19LNt6hPDQLjx8zUjuvGyo5nS3k75qIuJzJZW1/GHjAV7Zmo/bbZk3JYEHpo+gfw+Nc3tCBS4iPnNmcTe5LXMnxXHflSNI6BfpdLSgoAIXEa87UFrNsx/lsWJHAW63Ze6keO69criK28tU4CLiFdZatuaV86cP83hn91G6hnZh7qR4vnPFBQzuq+L2BRW4iHikvtHN29nFPPtRHp8WVNK3e1cemD6Cr188RGPcPuZRgRtjvgfcBVggC/imtbbWG8FExL8dKT/FK1vzeW37Ecqq6xnavztP3TSOuZPitWZJB2l3gRtj4oAHgDHWWpcxZjnwNeAvXsomIn6mscnN+7mlvPSPw2zaV4oBpo8eyLwpCVw+IpoumsfdoTwdQgkFIowxDUAkUOR5JBHxN7klVaxML2BVRiHHquoY0LMb9181gq9dNJjYqAin43Va7S5wa22hMeYXQD7gAtZba9ef/ThjzAJgAUBCQkJ7TyciHay0qo41mUWsTC8gp+gkoV0MVyRFc8uF8UwfPbDTbmPmTzwZQukD3AgMBSqA14wxt1trXzrzcdbapcBSgNTUVOtBVhHxsaraBt7dfYw1mUVs3FtKk9uSHNebx68fw/UTYvWmpJ/xZAjlaiDPWlsKYIxZCVwCvPRv/ysR8Ssnaxt4Z9dR1maVsGlfKfWNbgb1Cufuy4Zx86Q4Rg7s6XRE+QKeFHg+MNUYE0nzEMp0YLtXUomIT1W6Tpd2MR/uK6O+yU1M73BunzKEWcmDmJTQR29IBgBPxsC3GGNWAOlAI5BBy1CJiPifylMNrN9VwtqsYj7aX0ZDkyW2dzjzLx7CrOQYUgZHqbQDjEezUKy1jwOPeymLiHhZxal61uccZW12MZtbSjsuKoL/vCSRWckxTBwchTEq7UClOzFFgkx5TT3rc0pYm13Cx/vLaHRb4vtEcMe0ocxKjmF8fG+VdpBQgYsEgePVdaTlNI9pf3LwOE1uS0LfSO66bBizk2MYF9dLpR2EVOAiAaq0qo51OSW8nVXMPw4ex21haP/ufPtLw5g5LoaxsSrtYKcCFwkg1XWNrMsu4Y2dhWzeX4bbwrDo7tx75XBmJccwalBPlXYnogIX8XMNTW427S1lVUYh7+w+Sm2Dm4S+kdx75XCuGx/LyIE9VNqdlApcOrXVGYUsTsulqMJFbFQEC2ckMSclzulYAOQUVfK3bUd4M7OIE6ca6BMZxldSB3PjxDgmJWj2iKjApRNbnVHIopVZuBqaACiscLFoZRaAYyVeVdvAmswiXt16hKzCSrqFduGaMQO5KSWOy0dGa/0R+RwVuHRai9NyPyvv01wNTSxOy+3wAs88UsHLWw7zZmYxroYmRg3qyX9dP4abUuLpHRnWoVkCkT+/kvIlFbh0WkUVrjYd97aGJjfrskt4fnMe6fkVRHYN4caJsXxtcgITNFe71fzxlVRHUYFLpxUbFUHhOcra1+tbn6ipZ9m2fF785DDFlbUM6RfJ49eP4ZYL4+kZrqvttvKnV1IdTQUundbCGUmfu3IDiAgLYeGMJJ+cr7jSxdJNB1m2NZ/aBjfThvfjyTnjuDJpgNYg8YDTr6ScpAKXTuv01Zmvx04PH6/hDxsPsGJHAdbCjRPjuPvyoYwa1Mur5+msnHol5Q9U4NKpzUmJ89nL7IOl1fz2vf28sbOQ0JAufO2iBBZcPozBfSN9cr7OqqNfSfkTFbiIlxVXunjm3X0s315At9Au3HnpUO6+bBgDeoU7HS0oddQrKX+kAhfxkhM19SzZeIC/fHwILMyfOoT7rhqubcg6gC9fSfkzFbiIh2obmnhucx5L3j9ATX0jN6XE892rR2ioRHxOBS7STtZa1u86ylNv7Sa//BRXjx7A968dpT0kpcOowEXaIbekiv/+ew6b9x9nxIAevHjnZC4bEe10LOlkVOAibVBV28Av1+/lr58comd4GE/cMJZ5UxII1Rol4gAVuEgrpeWU8PgbORytqmXelAQeviaJPt27Oh1LOjEVuMh5FFe6ePyNHNbvOsqoQT1ZcvskUhL6OB1LRAUu8kXcbstLWw7z83W5NLrdPDpzFHdeOlRLuorfUIGLnMOR8lM88lomW/LKuWxEf56ak0xCP00LFP+iAhc5g7WWZVuP8ORbuwgxhp/fMp5bL4zX0q7il1TgIi2KK1384PUsNu0t5dLh/Xn6lvHEdYIFkSRwqcBFgDd2FvKj1dk0Nln+58axzJsyREu8it9TgUunVl3XyE/eyGZleiEXDunDL2+dQGL/7k7HEmkVFbh0WlkFlTzwagaHj9fw4PQR3H/VcN2QIwFFBS6djttteW5zHk+v20P/Ht145e6pTB3Wz+lYIm2mApdOpay6joeXZ7Jxbykzxg7k6bnjiYrU3ZQSmFTg0mlsO1TOvS+nU+lq4Mk545g3JUHTAyWgeVTgxpgo4M/AOMACd1hrP/FGMBFvsdby7Ed5/OztPST0jeSFOyYzOkb7UUrg8/QK/P8D66y1txhjugK6VU38SlVtA99f8SlvZ5dw7dhB/PzW8fQKD3M6lohXtLvAjTG9gMuB/wSw1tYD9d6JJeK5PSUnueeldPLLT/HYrNHcddlQDZlIUPHkCnwYUAo8b4yZAOwAHrTW1pz5IGPMAmABQEJCggenE2m9lekF/HBVFj3Dw3jlrilM0SwTCUKeTHoNBSYBS6y1KUAN8OjZD7LWLrXWplprU6OjtWOJ+FZDk5sfr87moeWZjI+P4q0HLlV5S9Dy5Aq8ACiw1m5p+XwF5yhwkY5yvLqOe15OZ2teOQsuH8b3ZyTpxhwJau0ucGttiTHmiDEmyVqbC0wHdnkvmkjr7So6yd1/3U5ZdR2/+epE5qTEOR1JxOc8nYVyP/ByywyUg8A3PY8k0jZvZxXz0PJMekWEsvxbFzNhcJTTkUQ6hEcFbq3dCaR6KYtIm7jdlt+8u49n3t1HSkIUf7z9Qgb0Cnc6lkiH0Z2YEpCq6xp5ePlO0nKOcsuF8Tx10zi6hYY4HUukQ6nAJeDkHz/F3X/dzr5jVfz4ujHcMS0x6OZ3r84oZHFaLkUVLmKjIlg4I0nj+vIvVOASUD7eX8Z3XknHWnjhjslcNiL4pqauzihk0cosXA1NABRWuFi0MgtAJS6fozlWEhCstbzw8SHmP7eV/j268ca904KyvAEWp+V+Vt6nuRqaWJyW61Ai8Ve6Ahe/V9/o5idvZPPqtiNcPXoAv/7qRHoG8XomRRWuNh2XzksFLn6ttKqOe17awfbDJ7jvyuE8dM3IoN+rMjYqgsJzlHWsNliWs2gIRfxWdmElN/7uI7KLKvntbSk8MiMp6MsbYOGMJCLCPj+jJiIshIUzkhxKJP5KV+Dil97MLGLhikz6RnZlxbcvYVxcb6cjdZjTb1RqFoqcjwpc/Irbbfnlhlx+//4BLkrsw//Ou5Dont2cjtXh5qTEqbDlvFTg4jeqahv43t928s7uY9w2eTBP3DCOrqEa5RP5Iipw8Qt5ZTXc/dft5JXV8N83jmX+1CFBd3OOiLepwMVxm/aWct8r6YR0Mbx05xQuvkDrd4u0hgpcHHN6s+Gfrt3NyIE9+dPXUxncV9uqirSWClwcUdvQxA9XZbEyvZCZ4wbxi1sn0L2bvh1F2kI/MdLhjp6sZcGLO8g8UsFD14zkviuHd4r53SLepgIPMv6+il1G/gm+9eIOauoa+eP8C5kxdpDTkUQClgo8iPj7KnavbT/CY6uyGdi7Gy/eOY2kQT2djiQS0DTJNoj46yp2dY1N/Gh1FgtXfEpqYh/W3HupylvEC3QFHkT8cRW7ogoX97ycTuaRCr71pWEs/LJ2ihfxFhV4EPG3Vew27y/j/mUZ1De6WTJvEjOTYxzJIRKsdCkURPxlFTtrLUs+OMD8Z7fQt3tXVt87TeUt4gO6Ag8i/rCKXVVtA4+8lklazlFmj4/h53PHa363iI/oJyvIOLmKXXZhJfcvyyC//BQ/mj2aOy8dqvVMRHxIBS4eO71f5U/X7qFP9zBeuWsKU4ZpPRMRX1OBi0cqTzWwcEUm63cd5apRA/jFrRPo272r07FEOgUVuLTbjsMneGBZBseqavnR7NHcMW2obokX6UAqcGmzJrflj5sO8Mv1e4mNCue1b1/CxMFRTscS6XRU4NImR8pP8fDyTLYeKmd2cgw/m5tMr/Awp2OJdEoqcGkVay2vbS/giTdz6GIMv/rKBG5KidMsExEHqcDlvMqq61i0MosNu44ydVhffnHrBOL7aOMFEaepwOXfWpddzI9WZ3PS1ag3KkX8jMcFbowJAbYDhdba6zyPJP7g2MlafvJGDutyShgb24uX75qoFQRF/Iw3rsAfBHYDvbzwXOKw02PdT761i9pGNz+4dhR3XTaUMK0gKOJ3PCpwY0w8MBt4CnjIK4nEMYeP1/DDVVls3n+cyYl9+dncZC6I7uGV5/b3nYJEApGnV+C/Ab4PfOFra2PMAmABQEJCgoenE19w1Tex5IP9/GHTQbqGdOHJOeP4j8kJXhvr9vedgkQCVbtfFxtjrgOOWWt3/LvHWWuXWmtTrbWp0dHR7T2d+IC1lrScEq7+1UaeeW8/144dxDsPfYnbpw7x6huV/rpTkEig8+QKfBpwgzFmFhAO9DLGvGStvd070cSXDpZW88Sbu9i4t5SRA3uw7O6pXHyBbxag8sedgkSCQbsL3Fq7CFgEYIy5AnhE5e3/yqrreObdfbyyJZ+IsBB+fN0Yvn7xEJ++SelvOwWJBAvNA+8kTtU38uyHefxh4wFqG93cNnkwD04fSXTPbj4/98IZSZ8bAwdndgoSCTZeKXBr7QfAB954LvGuhiY3K3YU8OsNezlWVceXxwzkBzNHeW12SWv4w05BIsFIV+BBqr7RzevpBfz+/f0UnHCRkhDF7+dN4qLEvo7kcXKnIJFgpQIPMqeL+3fv7aewwsWEwVH8z43juCIpWgtPiQQZFXiQqK5r5G/bjvDcR3mfFfeTN43jipEqbpFgpQIPcCWVtTz/cR6vbMmnqraRyYl9VdwinYQKPEBlFVTy/Md5rNlZhNtaZibHcPdlw7QzjkgnogIPIKfqG/l7ZjEvbzlMZkElkV1DmH/xEO6YNpTBfbU+t0hnowIPAPuOVvHylnxeTy+gqraREQN68MQNY5mTEkfvCG1nJtJZqcD91Imaev7+aRGvpxey80gFXUO6MDN5EPOmDOGixD4a3xYRFbg/qW9080HuMVamF/LunqM0NFmSBvbkh7NGMXdSPP16+P6uSREJHCpwhzW5LdsPlfNWVjFvZhZx4lQD/Xt05esXJ3LzpDjGxPTS1baInJMK3AFNbsvWvHLWZhWzLqeE0qo6uoV24ZoxA5k7KZ7LRvQnVDvgiMh5qMA7SGOTm62HWko7+yhl1XWEh3XhqlEDmJUcw5VJA+jeTf87RKT11Bg+1NjkZkte8/BIWnYJx2vqiQgL4arRA5g1LoYrR0UT2VX/C0SkfdQeXtbY5OaTg8dZm1VMWs5RymvqiewawlWjBjA7OYYrkgYQ0TXE6ZgiEgRU4F7Q0OTm4wPHWftpMet3lXDiVAPdu4YwffRAZiXHcEVSNOFhKm0R8S4VeDvVN7rZvL+MtVnFrN91lEpXAz26hXL16AHMTI7hSyNV2iLiWyrwNqhrbOLDvWWszS5mw66jVNU20jM8lGvGDGTWuBguHdFfpS0iHUYFfh5ut2VLXjmrMwpZm11MVW0jvcJDmTF2ELOTY7hkeD+6haq0RaTjqcC/wJ6Sk6zKKGTNziKKK2vp3jWEGeMGccOEWC65oD9dQzVPW0ScpQI/Q3lNPSvTC1ixo4A9JVWEdjF8aWQ0i2aN5prRAzV7RET8SqcvcLfb8vGB4yzbls/6nBIamiwTB0fxPzeOZfb4WPp27+p0RBGRc+q0BV5WXcffth3h1W35HCl3ERUZxu1Th/C1ixJIGtTT6XgiIufV6Qo8p6iS5zcfYs3OIuqb3Fw8rB+PfDmJGWMHaQaJiASUTlHgTW7Lhl1HeW5zHlvzyokIC+GrFw3mG5ckMnxAj3Y95+qMQhan5VJU4SI2KoKFM5KYkxLn5eQiIl8sqAu8vtHNyvQClmw8wOHjp4iLiuCxWaP5ykWDPdrJZnVGIYtWZuFqaAKgsMLFopVZACpxEekwQVngrvomXt2Wz9JNBymurGV8fG+WzJvENWMGemWZ1sVpuZ+V92fnbGhicVquClxEOkxAFHhrhytc9U288Mkh/rTpIMdr6pk8tC9Pzx3PZSP6e3VThKIKV5uOi4j4gt8XeGuGK+ob3fxtWz7PvLef0qo6Lh8Zzf1XDeeixL4+yRQbFUHhOco6NirCJ+cTETkXvy/wfzdccf2EWNZkFvLrDfvILz/F5MS+LJk3iVQfFfdpC2ckfe6XCkBEWAgLZyT59LwiImfy+wL/omGJwgoXs5/5kD0lVYyJ6cXz37yIK0ZGd8j+kaev/DULRUSc5PcF/kXDFQA19Y08c1sK1yXH0KVLx278OyclToUtIo5q95QMY8xgY8z7xpjdxpgcY8yD3gx22sIZSUSc4wab2ckxbPjel7hhQmyHl7eIiD/w5Aq8EXjYWptujOkJ7DDGbLDW7vJSNqD5SjevrIbfvrcPt20ea3505ii+cUmiN08jIhJw2l3g1tpioLjl4ypjzG4gDvBqgQN8/eIhpOef4KFrRpKS0MfbTy8iEpC8MgZujEkEUoAt3ni+s/Xr0Y0X75zii6cWEQlYHt+WaIzpAbwOfNdae/Ic/77AGLPdGLO9tLTU09OJiEgLjwrcGBNGc3m/bK1dea7HWGuXWmtTrbWp0dHRnpxORETO4MksFAM8C+y21v7Ke5FERKQ1PLkCnwbMB64yxuxs+TPLS7lEROQ8PJmF8hGgCdgiIg7R1uoiIgFKBS4iEqBU4CIiAcpYazvuZMaUAoc9eIr+QJmX4niLP2YC5Wor5Wob5WobT3MNsdb+yzzsDi1wTxljtltrU53OcSZ/zATK1VbK1TbK1Ta+yqUhFBGRAKUCFxEJUIFW4EudDnAO/pgJlKutlKttlKttfJIroMbARUTknwLtClxERFqowEVEApTfF7gx5jljzDFjTLbTWc7UUXuCtpUxJtwYs9UYk9mS6wmnM53JGBNijMkwxvzd6SynGWMOGWOyWhZk2+50ntOMMVHGmBXGmD0t32cX+0GmpDMWr9tpjDlpjPmu07kAjDHfa/mezzbGLDPGhDudCcAY82BLphxvf638fgzcGHM5UA381Vo7zuk8pxljYoCYM/cEBeZ4e0/QduQyQHdrbXXLeu0fAQ9aa//hZK7TjDEPAalAL2vtdU7ngeYCB1KttX51A4gx5gXgQ2vtn40xXYFIa22F07lOM8aEAIXAFGutJzfoeSNLHM3f62OstS5jzHJgrbX2Lw7nGge8CkwG6oF1wD3W2n3eeH6/vwK31m4Cyp3OcTZrbbG1Nr3l4yrg9J6gjrLNqls+DWv54xe/pY0x8cBs4M9OZ/F3xphewOU0r7mPtbben8q7xXTggNPlfYZQIMIYEwpEAkUO5wEYDfzDWnvKWtsIbARu8taT+32BBwJf7wnaVi3DFDuBY8AGa61f5AJ+A3wfcDsd5CwWWG+M2WGMWeB0mBbDgFLg+ZYhpz8bY7o7HeosXwOWOR0CwFpbCPwCyKd5s/VKa+16Z1MBkA1cbozpZ4yJBGYBg7315CpwD51vT1AnWGubrLUTgXhgcsvLOEcZY64Djllrdzid5RymWWsnATOBe1uG7ZwWCkwCllhrU4Aa4FFnI/1Ty5DODcBrTmcBMMb0AW4EhgKxQHdjzO3OpgJr7W7gaWADzcMnmUCjt55fBe6B1uwJ6qSWl9wfANc6HAWad3C6oWW8+VWad3J6ydlIzay1RS1/HwNW0Txe6bQCoOCMV08raC50fzETSLfWHnU6SIurgTxrbam1tgFYCVzicCYArLXPWmsnWWsvp2KsoXAAAAEjSURBVHk42Cvj36ACbzd/3RPUGBNtjIlq+TiC5m/sPc6mAmvtImttvLU2keaX3u9Zax2/QjLGdG95E5qWIYov0/yy11HW2hLgiDEmqeXQdMDRN8jPcht+MnzSIh+YaoyJbPnZnE7z+1KOM8YMaPk7AbgZL37d2r2lWkcxxiwDrgD6G2MKgMettc86mwr4556gWS3jzQA/tNaudTATQAzwQssMgS7Acmut30zZ80MDgVXNP/OEAq9Ya9c5G+kz9wMvtwxXHAS+6XAeAFrGcq8BvuV0ltOstVuMMSuAdJqHKDLwn9vqXzfG9AMagHuttSe89cR+P41QRETOTUMoIiIBSgUuIhKgVOAiIgFKBS4iEqBU4CIiAUoFLiISoFTgIiIB6v8AGa66DFrw0DwAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fit_and_plot([1, 1, 1, 1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What went wrong? \n", "\n", "The initial values we chose were not sufficiently close to the \"correct\" values. This shows the first main problem with nonlinear regression: there may be multiple solutions which are returned based on the initial guess.\n", "\n", "| Linear regression| Nonlinear regression |\n", "|--|--|\n", "| single correct solution | multiple solutions possible depending on initial guess |\n", "| requires no initial guess | requires initial guess | \n", "| Never returns the \"wrong\" local minimum solution | Sometimes claims success with \"wrong\" answer |\n", "| less flexible in functional form | more flexible functional form |\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's try a different starting point. Remember the initial data were generated with $\\beta=[2, 3, 1, 1]$" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2., 3., 1., 1.])" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXzU1b3/8dfJvpGEkAWSEEhYwg4BBARBBRUQFdTaatWqVelib3tra6v1tr29m/2VLrfLrRXR4ooroFVJQEUEkT1AEiAQQshCJgmE7PvM+f2RxAIGSGa+M9/vTD7Px4MHZDLzPR9j8s6Z8z2L0lojhBDC+/iZXYAQQgjnSIALIYSXkgAXQggvJQEuhBBeSgJcCCG8VIAnG4uNjdXDhw/3ZJNCCOH19u7de1prHXfh4x4N8OHDh7Nnzx5PNimEEF5PKXWyp8dlCEUIIbyUBLgQQngpCXAhhPBSEuBCCOGlLhvgSqnnlVKVSqnccx5boZQ6opQ6qJRap5SKdm+ZQgghLtSbHvhqYNEFj20CJmitJwFHgScMrksIIcRlXHYaodb6U6XU8Ase23jOhzuArxhblhBCeL/12WWsyMrnVE0zidGhPLYwnWUZSYZd34h54N8EXr/YJ5VSy4HlACkpKQY0J4QQ1lZR18IfNh3lrb2ldDg6t+wuq2nmibU5AIaFuEsBrpR6EugAXrnYc7TWK4GVANOnT5fNx4UQPquupZ3/ef/wecF9ruZ2Oyuy8g0LcKdnoSil7gNuAu7WciqEEKKfKzrdyJI/beXNvaXcM2vYRZ93qqbZsDad6oErpRYBPwWu1lo3GVaNEEJ4oeIzTXz1mc/pcGje+NaVTBs2kE2HKijrIawTo0MNa7c30wjXAJ8D6UqpUqXUg8BfgAHAJqXUfqXU3wyrSAghvEh9SzsPvrCb1g4Hry2fxbRhAwF4bGE6oYH+5z03NNCfxxamG9Z2b2ah3NXDw88ZVoEQQnixX76TR+HpRl765gxGJwz44vHucW6rz0IRQoh+6cNDFazNLuP780cye2Tslz6/LCPJ0MC+kCylF0IIJzS0dvCzdTmMGTyA780fZUoN0gMXQggnPLPlOJX1rTxz7zSCAszpC0sPXAgh+qi8tplntxZy8+REMlIGmlaHBLgQQvTRnz4qwOGAnxg4o8QZEuBCCNEHttoW3t5byh3TkxkaE2ZqLRLgQgjRB6u2FmLXmm9fPcLsUiTAhRCit842tvHKzmKWTk40vfcNEuBCCNFrr+4qprndzrevMb/3DRLgQgjRK3aH5tWdxcweMei8FZdmkgAXQohe2HykkrKaZu69xE6DniYBLoQQvfDSjpMkRAZz3bgEs0v5ggS4EEJcxskzjWw5WsVdM1II9LdObFqnEiGEsKjXdpfg76e48wprHQspAS6EEJfgcGjWZ5cxb1Qsg6NCzC7nPBLgQghxCTsKz1Be28JtU5PNLuVLJMCFEOIS3t5XxoDgAK630M3LbhLgQghxEU1tHWzILWfJpCGEXHA8mhVIgAshxEVk5dloarNbcvgEJMCFEOKi1u4rY2hMKNOHmbfn96VIgAshRA+qG9vYfvwMt0xOxM9PmV1OjyTAhRCiBx8eqsDu0CyeMMTsUi5KAlwIIXqQmWcjeWAo4xMjzS7loiTAhRDiAvUt7Ww7dppF4wejlDWHT0ACXAghvmRzfhVtdgeLJgw2u5RLkgAXQogLZOXaiBsQzFQTT5zvjcsGuFLqeaVUpVIq95zHYpRSm5RSx7r+tvZ/pRBC9FJLu53N+ZXcMC7BsrNPuvWmB74aWHTBY48DH2mtRwEfdX0shBBe79OjVTS12S09+6TbZQNca/0pUH3Bw0uBF7r+/QKwzOC6hBDCFJl5NqJCA5mZFmN2KZfl7Bh4gta6HKDr7/iLPVEptVwptUcptaeqqsrJ5oQQwv3a7Q4+PFTBdWMTLHVww8W4vUKt9Uqt9XSt9fS4uDh3NyeEMIjWmsq6Fsprm3E4tNnleMSOwjPUtXRYfvZJtwAnX1ehlBqitS5XSg0BKo0sSghhnrqWdp7beoI1u4qprG8FIDoskGVTkhgRF87fthRyqqaZxOhQHluYzrKMJJMrNk5mro2wIH/mjoo1u5RecTbA3wXuA37d9fc7hlUkhDDNgZIaHnl1H6Vnm1kwJp55o+Pw91PsKDzDSztOYj+nJ15W08wTa3MAfCLE7Q5NVl4F16bHW3Lr2J5cNsCVUmuAa4BYpVQp8Es6g/sNpdSDQDFwhzuLFEK43+6iau5/fhfRYUG8/Z3ZTDtnB757Zg1jxn9/+EWPvFtzu50VWfk+EeDZxWc53dDKQi8ZPoFeBLjW+q6LfGqBwbUIIUxSUNnAN/++m4SoENY8PIuEyC+f/Vh1QXh3O1XT7O7yPGJDro0gfz/mj7nonAzLsf5tViGEWzW1dbD8pT0EBfjx0oMzewxvgMTo0D497k201mTm2pg7KpaIYGdHlj1PAlyIfu43mfkUVjXy57sySLpEGD+2MJ3QC8aGFfDItSPcXKH75Z2qo6ym2auGT0ACXIh+bU9RNS98XsR9Vw5j9shLz7xYlpHEU7dNJCk6FAXERgShFOwvqfFIre6UmWvD309x3VjrHVx8Kd7zXkEIYSi7Q/Pzd/JIjArlJ4vG9Oo1yzKSzrth+f8yj/D0J8dZOiWJOZf5BWBlmXk2ZqbGEBMeZHYpfSI9cCH6qXXZZRwur+Oni8cQ7uS47w8WjGL4oDCeXJdDW4fD4Ao9o6CynoLKBq9ZvHMuCXAh+qGWdju/25jP5OQobpro/KZNIYH+/GrpBIrONPHi50WG1edJmbk2ABaOlwAXQniBN/eUUF7bwk8XjXF5y9SrR8cxb3Qcf/roGGcb2wyq0HMy82xMTYm+6OwbK5MAF6Kf6bA7eObTQjJSorlyxCBDrvnkjWNpaO3gzx8XGHI9TympbiK3rM4rh09AAlyIfue9g+WUnm3mu9eMNOy8x/TBA7hj2lBe3nkSW22LIdf0hKw87x0+AQlwIfoVrTV/23Kc0QkRLDB4xeH35o/E4ei8vrfIyrMxdkgkwwaFm12KUyTAhehHdhRWc8RWz0Nz0ww/LmxoTBi3T03m1V3FXtELr6xvYc/Jsyzy0t43SIAL0a+8vPMkUaGB3DI50S3Xf+Ra7+mFb8yrQGu8dvwbJMCF6Dcq61vIyrXxlWnJbtsuNWVQGLdNTWLNrmLONPS8+ZVVZOXZSI0NZ3RChNmlOE0CXIh+4o3dJXQ4NHfPTHFrO8vnpdHa4eDFz0+6tR1XnG1sY/vxMyyeMNiwG7lmkAAXoh+wOzRrdpUwZ+Qg0uLc2+McGT+A68bG89KOkzS32d3alrM2Ha7A7tDc6MIiJiuQABeiH9h+/DRlNc18fcYwj7T38Nw0qhvbeGtfqUfa66sNOeUkDwxlfGKk2aW4RAJciH5g7b4yIkMCWDDWM4cVzEiNYfLQaFZtLTzvGDYrqG1uZ1vBaW6cOMSrh09AAlwIn9fY2kFmro0lkxI9dtajUorlc9M4eaaJTYdsHmmztz46XEG7XbPYi2efdJMAF8LHZeXZaG63c9tUz55buWjCYIbGhPLMp4UebfdyNuTaGBIVwuTkaLNLcZkEuBA+bl12GUNjQpl+ziHFnuDvp3joqjSyi2vYe7Lao21fTENrB1uOVrFowmDDFzKZQQJcCB9mq21hW8Fpbs1INmW8947pyUSFBvLspyc83nZPPj5SSVuHw+tnn3STABfCh717oAyt4dYMzw6fdAsLCuDrM1PIOmTj5JlGU2o414accuIHBDMtxbPvRtxFAlwIH/ZBjo2JSVGkxpq3WdP9s4cT4Kf4+2dFptUAUNfSzkdHKrlx4hCfGD4BCXAhfFZ5bTP7S2pM3+sjITKEmycn8saeEmqb2k2rIzPXRluHg6VT3LMPjBkkwIXwURvzKgBrbNb00FVpNLXZeWWXecvr39lfxrBBYUwZ6v2zT7q5FOBKqR8qpfKUUrlKqTVKKe87k0gIH5WZa2NUfAQj3Lx0vjfGJUYyZ+QgXtheZMrhxxV1LWw/foalU5K8fvHOuZwOcKVUEvB9YLrWegLgD9xpVGFCCOdVN7ax88QZS/S+uz00N42KulbeO3jK423/48AptIZlPjR8Aq4PoQQAoUqpACAM8Pz/GSHEl3x4qAKHttZRYdeMjmNUfATPbj2B1p5dXr9+fxmTkqPcvpGXpzkd4FrrMuC3QDFQDtRqrTde+Dyl1HKl1B6l1J6qqirnKxVC9Fpmns1ymzUppXhobiqHy+vYfvyMx9o9VlFPblkdS6eYM5XSnVwZQhkILAVSgUQgXCl1z4XP01qv1FpP11pPj4uLc75SIUSv1Le0s+3YaRaNt95e10unJBEbEcSqrZ5bXv/a7hIC/ZXPDZ+Aa0Mo1wEntNZVWut2YC0w25iyhBDO2pxfRZvdYanx724hgf7cO2s4m/OrKKisd3t7Le123t5Xyg3jBjMoItjt7XmaKwFeDMxSSoWpzl/zC4DDxpQlhHBWVq6NuAHBTLXoasN7ZqUQHODHqq3uX16flWejpqmdu2a49xQiswQ4+0Kt9U6l1FvAPqADyAZWGlWYEK5o63Dw7oFTbD5SSUVdCzHhQVw1KpbbpyYTHuz0t73ltbTb2Zxfya0ZSZZdbTgoIpjbpyXz1t5SfnRDOnED3NczXrOrmKExocweMchtbZjJpVkoWutfaq3HaK0naK3v1Vpb+xRT0S/sPXmW636/hR+/eYDs4rME+vtxtKKeX7yTx9UrNvNBTrnZJbrN1mOnaWqzW3L45FzfnJNKW4eDFz8vclsbJ043sqOwmjuvSLHsLzNX+W5XRPRL7x44xaOv72dIdAh/v/8KrkmP++JG3r7is/zq3Ty++8o+vnftSH50w2jL3eRzVWaujciQAGalWbvHOTI+gkXjB7P6syIevCqV6LAgw9tY/dkJAv0Vd0xLNvzaViFL6YXP2HSogh++vp9pwwby3r/M5dox8ecF9NSUgbz57dncecVQ/rK5gBVZ+SZWa7x2u4MPD1dw3bgEAv2t/6P9w+tH09DWwUo3HPhQ09TGG3tKWTolifhI310gbv3/y0L0QmFVAz98fT8TEiN5/v4riAoN7PF5QQF+PHXbRO6akcJfPznOW3uteeiuM3YUnqG2uZ3FE7xjr+v0wQO4aVIiq7cXcabB2NHXV3YW09xu56G5qYZe12okwIXXa+2w891X9hHor/jrPdMue5NSKcV/Lh3P7BGD+Nm6HHJKaz1UqXtl5toIC/Jn7qhYs0vptX+9bhQt7Xb+tuW4Ydds7bCzensRc0fFMmawdRYyuYMEuPB6z2wp5Iitnt99dTJJ0aG9ek2Avx9/+fpUwgL9ufWvnzH88feZ8+uPWZ9d5uZq3cPu0GTlVXBterzHDi42woi4CJZlJPHi5ycpq2k25JprdhZTVd/Kt68eYcj1rEwCXHi1E6cb+cvmAm6aNIT5YxL69NpPj1bR1Ganw9G5L0dZTTNPrM3xyhDPLj7L6YZWFlp89klPfnRDOkrB/3zg+jKS5jY7f9l8nJmpMT47dfBcEuDCq/3PB4cJ8vfjFzeN6/NrV2Tl02Y/f2vT5na7V97czMy1EeTvx7Xp3rddRVJ0KN+5eiTvHyzncyf3SFmfXcacX3/M2F9kcrqhlZmpMT43w6gnEuDCa+0vqWHToQq+NS/NqZkGpy7ylv1ij1uV1prMPBtXjYplQEjPN2+t7ltXp5EUHcov382ltcPep9euzy7jibU55w3BPLv1hFe+k+orCXDhtX63MZ+Y8CAeuMq5mQaJFxkvv9jjVpV3qo7Ss80sstDWsX0VEujPf906gaMVDfx+09E+vXZFVj7N7eeHvre+k+orCXDhlXadqGbrsdN85+oRRDi5NP6xhemEXnDDTwGPXj/agAo9JyvPhp+C68b17R6A1VybHs9dM4ay8tPCPg2l+Mo7KWdIgAuv9MyW4wwKD+LeK4c5fY1lGUk8ddtEkqJDUUBMWBAa8OxRA67LzLUxM3UQMeHGr2b0tCeXjCM1NpxHXt1HSXVTr14TH9nzXire9k7KGRLgwuscr2rgoyOV3HvlMJenzC3LSOKzx+dz4tdL2Pvz65iQFMmfPjpGu93z5zY6o6CygWOVDZbf+6S3IoIDWPWN6bTbHTz4wm5OX2aBT01TW4+Phwb689jCdHeUaCkS4MLrPL/tBEEBftwzy/ned0+UUjx6/WiKq5u8ZoVmVp4NgBvGe/fwybnS4iJ45p5pFFc3cefKHRRWNfT4vPLaZu56didnG9t55NoRX7yTSooO5anbJrIsw/dO4LmQbGYlvMrZxjbe3lfKbRlJxLphg/5r0+OZMjSav3xcwFemJVt+T5GsPBtThkYzJMq3hgtmj4xl9QMz+M7Le7n5z9v41tUj+NoVQ0mIDKG2uZ11+0r5w4fHsDs0q+6bzrzRcTy2cIzZZXucBLjwKm/vK6Wl3cH9c4a75fpKKb6/YCTfXL2H9w6e4tYM6+5kV3q2iYOltTy+2DeDa1baID74wVx+vj6P3286yu83HSUiOICG1o6uz8fwX8smMjLetw4q7gsJcOE1tNa8truEjJRot+5xcc3oeEYnRPDMlkKWTUmy7IKQzNzO4ZPFPjL+3ZMhUaGsum86Ryvq+fRoFWU1zcSEBTFnVCwZQ6Mt+//GUyTAhdfYe/IsBZUN/Ob2SW5tx89PsXzeCH785gE+OVrFtenxbm3PWRtybYwbEsmwQeFml+J2oxMGMDphgNllWI61B/iEOMeaXSWEB/mzZJL7t0u9ZXIiQ6JCeMbAXfKMZKttYe/Jsz7d+xaXJwEuvEJdSzvv55zililJHjnTMijAjwevSmVHYTX7S2rc3l5fdc8+WTzRO/b+Fu4hAS68wnsHymlpd3DnFUM91uadM1KIDAmwZC/8g5xyRsVH9OsbeEICXHiJd/aXMSIunEnJUR5rMyI4gLtnDSMrz0bxmd6tCvSEqvpWdhdVS+9bSIAL6yuvbWZXUTW3TPb8jJD7Zw/H30/x/GcnPNrupWw8ZMOhfXv2iegdCXBhee8dKEdruGVKosfbTogM4ZbJSby+u+Siy7Y9LTPXRmpsOGMGy6yM/k4CXFjeOwfKmJwcRWqsOdPlHp6XSnO7nVd2FpvS/rmqG9vYfvwMiyYM7vdzoIUEuLC4gsoGcsvquGWKeftajBkcybzRcazeXtTnwwaM9kFOOXaH5iYPTKUU1icBLizt3QOnUApuNjmwHp6bSlV9K+/sP2VqHe/uP8Wo+AjGDfHt09ZF77gU4EqpaKXUW0qpI0qpw0qpK40qTAitNe/uL2P2iEFOHZlmpKtGxjJm8ABWbS1Ea3N2DC8928SuomqWTkmU4RMBuN4D/yOQqbUeA0wGXD9WWogueafqKDrTxM2TPH/z8kJKKR6em8bRiga2HK0ypYZ/HCgHYKmJw0nCWpwOcKVUJDAPeA5Aa92mtbbekjXhtTJzbfj7KW6wyFmPN09OJCEymFVbzZlS+M7+MqamRDM0JsyU9oX1uNIDTwOqgL8rpbKVUquUUl+aJqCUWq6U2qOU2lNVZU7PRXinzDwbM1NjLHNUWFCAH/fPTmVbwWnyTtV6tO18Wz1HbPX94pAC0XuuBHgAMBV4WmudATQCj1/4JK31Sq31dK319Li4OBeaE/1JQWU9BRY8KuzrM1MID/LnOQ/3wtdll+Hvp7hRVl+Kc7gS4KVAqdZ6Z9fHb9EZ6EK4rHuv6xvGWSvAo0ID+doVKbx74BTltZ459bzd7uCtvaXMHxPvllOIhPdyOsC11jagRCnVfXLoAuCQIVWJfi8zz0ZGSjSDo8ydfdKTB+YMx6E1qz8r8kh7Hx+p5HRDq0c38hLewdVZKP8CvKKUOghMAf7H9ZJEf1dS3URuWZ1l9/oYGhPGjROH8OrOYupb2t3e3mu7ikmIDObq0TIEKc7nUoBrrfd3jW9P0lov01qfNaow0emIrY439pTw+u5i9pfUmDYH2ZO697peaJHZJz15eG4a9a0dvL67xK3tnKppZsvRKr46fSgBFj9gWXieHKlmUdnFZ/n3d/M4UHr+bIdR8RH8dNEYrhuXYFJl7peZa2OsxY8Kmzw0mhmpMfz9syLunz3c8HBdn13Giqx8ymo6x9mjQwMNvb7wDfIr3YKe33aC25/eTlV9K/+xdDyf/Pgatv7kWlZ8ZRJKwUMv7uHn63OxO3yvN15Z18Le4rMssnDvu9vyuWmU1TTzQdcNV6Oszy7jibU5X4Q3wG83HmV9dpmh7QjvJz1wi/nrJwX8JjOfheMT+O0dkxkQ8s+e19CYMJZOSWJF1hGe3XqC6sY2/nRXBv5+vrOsOutQBVrD4onWD/D5Y+JJiwvnmS3HuXnSEMOWt6/Iyqe5/fxNs5rb7azIypd54OI80gO3kPcOnuI3mfncMjmRv9497bzw7hYU4MeTS8bx5I1jeT+nnF/9I8+nxsWzcm2kxYYzyguOCvPzUzxyzUjyTtWxwcBe+KmanqcnXuxx0X9JgFtESXUTP33rINOHDWTFHZMu26t+eF4a35qXxoufn+TNPaUeqtK9apra+LzwDAu9aK/rZRlJjE6I4Lcb8+mwOwy5ZmJ0aJ8eF/2XBLgF2B2aH71xAD+l+ONdGQQH+PfqdT9ZNIYr0wbxy3fzKKhscHOV7vfh4UrsDu0V49/d/P0UP74hncKqRt7eZ8wv0h/fMJoLf3+HBvrz2ML0nl8g+i0JcAtYs6uYXUXV/OLmcST1oZfl76f4w9emEBLox4/e2O/1NzWz8mwMiQrx6MHFRrh+XAIZKdH8YdMxmto6XL5eQmQIDt256lMBSdGhPHXbRBn/Fl8iAW6y2uZ2fr/pKDNTY/jKtOQ+v35wVAi/vHk8B0preXWX+Ud+OauprYNPj1Zxw7gErxk+6aaU4t+WjMVW18KfPipw6Vpaa3636SiDI0PY+bMFnPj1Ej57fL6Et+iRBLjJ/m9zAWeb2vj5TeOcDq6lUxKZPWIQv8k8QlV9q8EVesanR6to7XBYevHOpUwbFsMd05JZtbWQYxX1Tl/nHwfL2XvyLN9fMIqQwN4NpYn+SwLcRJX1LazeXsRtGclMSHJ+2EApxX8snUBLu53fZuUbWKHnbMyrICo0kBmpMWaX4rTHF48hPDiAn63LcWo4q6G1g/9+/xATkiL5mux7InpBAtxEq7aeoMPu4F/mj3T5WiPjI7h31nDe3FtCQaXzPUAztNsdfHi4ggVj4716ufigiGB+cdM4dhed5elP+j6U8tusfCrqWvmPpRN8am6/cB/v/WnxctWNbby84yS3TE5keKwxS8a/N38kYUEB/CbTu3rhOwurqWvp8Nrhk3PdNjWJmycn8ocPj7Gz8EyvX7c5v5LV24t4YM5wpqYMdGOFwpdIgJtk9WcnaG6388i1rve+u8WEBzFvdCwbD1Uw/PH3mfPrj71i+XVWno2QQD/mjfL+3faUUvzXsgkMGxTG8pf29urdUEFlPT9Yk82YwQP46aIxHqhS+AoJcBO0dth5ZWcxC8YkMCphgGHXXZ9dxseHK7/4uKymmSfW5lg6xB0OzaZDFVw9Oo7QIN+4aRcVGsgLD8wg0N+Pe1bt4oit7qLPPV7VwH3P7yYowI9nvzFdblyKPpEAN8EHOeWcaWzj/tnDDb3uiqx8WjrOXw3YvYeGVR0sq8VW12K5k3dcNTQmjJcfmoFG85WnP+eVnSfPW6mptead/WXc/vR2WtrtrH5ghhxWLPpMNrMywYufnyQtLpw5IwcZel1v3EMjK6/z5PkFY+PNLsVwYwZHsu67c/jRGwd4cl0uf/m4gNkjYgnwU+w+WU1hVSOTk6P4450Zht0HEf2LBLiH5ZTWkl1cwy9vdn7e98UkRoeetwXpuY9bVVaejVlpMUSHWePkeaMlRofy6sMz2Xiogjf3lLCtoAq7A8YMHsAj14xkWUaSzDgRTpMA97CXdhQRFuTP7U6surycxxam88TanPO2IvVXyrJ7aBRU1lNY1Wj4UJLVKKVYOH6wT8yyEdYiAe5Bja0dvHewnFsmJxLZw1axrupebr0iK59TNc2EBfvT3Ga37OKYrLwKoHMvESFE30mAe9CGXBtNbXan9jzprWUZSV8EeUl1E1ev2MwLnxfxxOKxbmvTWRvzbExOjmJIlHWHeISwMpmF4kFr95UybFAY04Z5ZqHG0JgwFk/oPD29sdX1XfKMVFLdxIHSWhZa9OR5IbyBBLiHlNU083nhGW7LSPbobnsPzk2lvqWDN/e49/T0vtqQWw7ATRMTTa5ECO8lAe4h6/aVonXnUmtPmpoykKkp0Tz/WZGl9gt//2A5k5KjSBkkc5+FcJYEuAdorXl7XxkzUmNMWazx0Nw0iqub2HSowuNt96R7+OTGiUPMLkUIryYB7gEHSms5cbqR2z3c++52w7gEkgeG8ty2QlPav9D7OZ3DJ0skwIVwiQS4B2zILSfAT7FovDmBFeDvxzfnpLK76Cz7S2pMqeFc7x8sZ3JylCwdF8JFLge4UspfKZWtlHrPiIJ8jdaarFwbV44YRFSY8XO/e+urVwxlQHAAq7aa2ws/eaaRnLJalkyS3rcQrjJiHvgPgMNApAHX8jn5FfUUnWli+bwRptYRERzAXTNTeG7bCUrPNpE80LO93/XZZazIyv9iqX+An7z5E8JVLv0UKaWSgSXAKmPK8T0bcmwoZY3Vht1L1ld/VuTRdtdnl/HE2pzz9mlZkZVv6W1uhfAGrnaD/hf4CeC42BOUUsuVUnuUUnuqqqpcbM77ZOXZuGJYDHEDgs0uhcToUJZMHMJru0uoa2n3WLsrsvLP258FrL/NrRDewOkAV0rdBFRqrfde6nla65Va6+la6+lxcd5/4kpfnDjdyBFbvaVWGz48N42G1g7e2O25hT3euM2tEN7AlR74HOAWpVQR8BowXyn1siFV+YisPBsAiywU4BOTo5iZGsPfPys674ABd7rYdrZW3uZWCG/gdIBrrZ/QWidrrYcDdwIfa63vMawyH7Ah18ak5CiSLBZUD81No6ymmQ9ybR5p77GF6QRdcNp8aKC/Zbe5FcJbyFQANzlV08yBkhpL7gG9YEw8abHhrNpaiMdf1o0AAA8ZSURBVNbuX16/LCOJ8YmRdG8BkxQdylO3Tfxi10QhhHMM2U5Wa/0J8IkR1/IVGy04fNLNz0/xzatS+bf1uewuOuv2/cLPNLSSd6qO+64czr/fMt6tbQnRn0gP3E0y82yMTohgRFyE2aX06PapyQwMC+RZDyzseWtvKW12B3fPTHF7W0L0JxLgbnCmoZVdJ6pZZMHhk26hQf7cO2sYHx6u4GhFvdvacTg0r+4qZkZqDKMSBritHSH6IwlwN9h0qAKHxlLTB3vywJxUwgL9+dNHx9zWxmfHT3PyTJP0voVwAwlwN8jMs5ESE8a4IdbeXWBgeBDfmD2c93PKOeamXvgL24uICQ+y5L0AIbydBLjB6lra+azgNIsmDPboyTvOenhuGqGB/vz54wLDr11QWc+Hhyv5xpXDCA7wN/z6QvR3EuAG+/hwJe12bcnpgz2JCQ/i3iuH8Y+DpyiobDD02is/LSQk0I9vXDnc0OsKITpJgBssM9dGQmQwGUOjzS6l15bPTSMs0J8VWUcMu2ZFXQvrssv46vShxIQHGXZdIcQ/SYAbqLnNzidHK1k4fjB+ftYfPuk2KCKY71wzgqy8CnYUnjHkmk9/chyHhoeuSjPkekKIL5MAN9CWo1W0tDssPX3wYh68Ko0hUSH89/uHcbh4+HFJdROv7DzJV6cny6HFQriRBLiBsvJsRIcFun1lozuEBvnzk0Xp5JTVsmZ3sUvX+sOHR/FTiu8vGGVQdUKInkiAG6Stw8GHhyu4fmwCAf7e+WVdNiWJ2SMG8dQHR5ze6nVPUTVr95Vx/5zhDImy1iZeQvga70waC9p+/DT1LR0snuh9wyfdlFL8+rZJ2B2aJ9fl9Hmjq7YOB0+szSEpOpTvz5fetxDuJgFukKw8GxHBAcweEWt2KS5JGRTGTxalszm/ilVbT/TptX/86CjHKhv4z2XjCQ82ZJ80IcQlSIAbwO7QbMyr4Nox8YQEev+ClftnD2fxhME8teEw246d7tVrNh+p5P82H+dr04cyf4z5538K0R9IgBtgd1E1ZxrbvHL2SU+UUqy4YzIj4yNY/tIe9p6svuTz9548yyOv7mPckEh+tVS2ixXCUyTADZCZayM4wI9r0n3nzM+I4ABefnAmCZEhfP3Znby9t7THMfENOeV847mdxA8IZvUDV/jEOxAhvIUMVLpIa01Wno15o+N8btw3PjKEN799Jd99ZR8/evMAr+8u4Y7pyaTFhVNe28Kbe0rZcrSKKUOj+ds904iPDDG7ZCH6Fd9KHBMcKK2lvLaFH9/gm+c7xkYE8+pDM3ltdwl/3VzAY28d/OJzA8MCeWLxGO6fM1w2qxLCBBLgTlqfXcaKrHzKuuZLt3bYTa7IfQL8/bhn1jDunpnCEVs9troWBoUHMW5IpNfOeRfCF0iAO2F9dhlPrM2huf2fof2f7x0mLCjApw/qVUoxdkgkYy2+z7kQ/YV0n5ywIiv/vPAGaG63syIr36SKhBD9kVf0wLuHK07VNJMYHcpjC9NN7elebJm5s8vPhRDCGZbvgXcPV5TVNKOBsppmnlibw/rsMtNqSozueY+Piz0uhBDuYPkAt+JwxWML0wkOOP9LFxroz2MLfXMmihDCmiwf4FYcrliWkXTeop2k6FCeum2iT9/AFEJYj9Nj4EqpocCLwGDAAazUWv/RqMK6JUaHfjFV78LHzVRY1cistBheW36lqXUIIfovV3rgHcCPtNZjgVnAI0qpccaU9U+PLUwn9ILl2WYPVxyrqOdYZQM3ThxiWg1CCOF0D1xrXQ6Ud/27Xil1GEgCDhlUG8AXwxLdi2YU8G9Lxpo6XPGPg+UohdecPC+E8E2GTCNUSg0HMoCdRlzvQssykliWkcQRWx2L/ncrNc3t7mimV7TWrM8uY/aIQSTI3h9CCBO5fBNTKRUBvA38q9a6rofPL1dK7VFK7amqqnKprTGDI5k7KpbV24tMW7qeXVJDcXUTS6fIDUshhLlcCnClVCCd4f2K1nptT8/RWq/UWk/XWk+Pi3N9u9WH56ZRVd/Ku/tPuXwtZ7yTXUZwgB+LJsjwiRDCXE4HuFJKAc8Bh7XWvzeupEubOyqW9IQBPLftRJ/PbHRVu93BPw6Wc93YBCJDAj3athBCXMiVHvgc4F5gvlJqf9efGw2q66KUUjw4N5Ujtnq29vK4L6NsO3aa6sY2me8thLAEpwNca71Na6201pO01lO6/nxgZHEXs3RKInEDgnl2a6EnmvvC+v1lRIcFcvVo3zl5RwjhvSy/ErMnwQH+3HflMLYeO02+rd4jbdY0tbEh18ZNk4YQFOCVXzYhhI/x2iS6e+YwQgL9WOWhXvi67DLaOhzcNSPFI+0JIcTleG2ADwwP4o5pQ3ln/ykq61vc2pbWmjW7ipmcHMX4xCi3tiWEEL3ltQEO8OBVqbQ7HLy4/aRb29lXfJajFQ18fab0voUQ1uHVAT48Npzrxybw8s6TNLV1uK2dV3eWEBEcwE2TEt3WhhBC9JVXBzjA8nlp1DS188qOYrdc/0xDK+8dPMXSKYmEB3vFAUZCiH7C6wN8+vAYrhoZy9+2HKex1fhe+Es7TtLa4eCBOamGX1sIIVzh9QEO8OgNoznT2MYLnxcZet3mNjsvfn6S68bGMzI+wtBrCyGEq3wiwKemDOTa9DhWflpIfYtxOxW+va+U6sY2Hp6bZtg1hRDCKD4R4ACPXp9OTVM7T39y3JDrtXU4eObT40xOjmJGaowh1xRCCCP5TIBPTI7i1owkVm07QfGZJpev9/qeEkqqm/nh9aPp3LdLCCGsxWcCHODxxWMI8FP81/uuHQrU3Gbnzx8d44rhA2XfEyGEZflUgCdEhvDItSPZeKiCjw5XOH2dpz8poLK+lccWjpHetxDCsnwqwAEemptKesIAnlibQ21T329oFlQ28PSW4yybkihj30IIS/O5AA8O8Od3X51MdWMbT67P6dOhD3aH5mdrcwgJ9OfJJePcWKUQQrjO5wIcYEJSFI/eMJr3Dpb3ac/wP350jF1F1fz7zeOJGxDsxgqFEMJ1Prs2/DtXjyC3rJZfbzhCUnQYSyYNueTz39lfxp8+OsbtU5O5fVqyh6oUQgjn+WyAK6X47R2Tqapv5fuvZVPf0s6dF9nLe+2+Un769kFmpsbw37dO8HClQgjhHJ8cQukWFhTA3x+YwewRg3h8bQ7feXkv+bb6L8bFS6qbePSN/Tz6xgGmD4vh2fumExLob3LVQgjROz7bA+8WERzA6gdm8PQnBfz1k+NsyLURGxFMoL+ivLaFIH8/vnvNCH54/WgC/X3695kQwsf4fIAD+Pspvjd/FHfOSCErz8b+4hrsWjNm8ABumpRIYnSo2SUKIUSf9YsA7xYbEczdM4dx98xhZpcihBAukzEDIYTwUhLgQgjhpSTAhRDCS0mACyGEl3IpwJVSi5RS+UqpAqXU40YVJYQQ4vKcDnCllD/wf8BiYBxwl1JKdoASQggPcaUHPgMo0FoXaq3bgNeApcaUJYQQ4nJcCfAkoOScj0u7HjuPUmq5UmqPUmpPVVWVC80JIYQ4lysLeXo6quZLm29rrVcCKwGUUlVKqZMutBkLnHbh9e5gxZpA6uorqatvpK6+cbWuHlcfuhLgpcDQcz5OBk5d6gVaa5cOmFRK7dFaT3flGkazYk0gdfWV1NU3UlffuKsuV4ZQdgOjlFKpSqkg4E7gXWPKEkIIcTlO98C11h1Kqe8BWYA/8LzWOs+wyoQQQlySS5tZaa0/AD4wqJbeWOnBtnrLijWB1NVXUlffSF1945a6VF8O/RVCCGEdspReCCG8lAS4EEJ4KcsHuFLqeaVUpVIq1+xazqWUGqqU2qyUOqyUylNK/cDsmgCUUiFKqV1KqQNddf3K7JrOpZTyV0plK6XeM7uWbkqpIqVUjlJqv1Jqj9n1dFNKRSul3lJKHen6PrvSAjWld32duv/UKaX+1ey6AJRSP+z6ns9VSq1RSoWYXROAUuoHXTXlGf21svwYuFJqHtAAvKi1tsyR8UqpIcAQrfU+pdQAYC+wTGt9yOS6FBCutW5QSgUC24AfaK13mFlXN6XUo8B0IFJrfZPZ9UBngAPTtdaWWgCilHoB2Kq1XtU1VTdMa11jdl3duvZDKgNmaq1dWaBnRC1JdH6vj9NaNyul3gA+0FqvNrmuCXRuMzIDaAMyge9orY8ZcX3L98C11p8C1WbXcSGtdbnWel/Xv+uBw/SwlYCn6U4NXR8Gdv2xxG9ppVQysARYZXYtVqeUigTmAc8BaK3brBTeXRYAx80O73MEAKFKqQAgjMssLPSQscAOrXWT1roD2ALcatTFLR/g3kApNRzIAHaaW0mnrmGK/UAlsElrbYm6gP8FfgI4zC7kAhrYqJTaq5RabnYxXdKAKuDvXUNOq5RS4WYXdYE7gTVmFwGgtS4DfgsUA+VArdZ6o7lVAZALzFNKDVJKhQE3cv4KdpdIgLtIKRUBvA38q9a6zux6ALTWdq31FDq3N5jR9TbOVEqpm4BKrfVes2vpwRyt9VQ6t0Z+pGvYzmwBwFTgaa11BtAIWGbP/a4hnVuAN82uBUApNZDO3VBTgUQgXCl1j7lVgdb6MPD/gE10Dp8cADqMur4EuAu6xpjfBl7RWq81u54Ldb3l/gRYZHIpAHOAW7rGm18D5iulXja3pE5a61Ndf1cC6+gcrzRbKVB6zrunt+gMdKtYDOzTWleYXUiX64ATWusqrXU7sBaYbXJNAGitn9NaT9Vaz6NzONiQ8W+QAHda183C54DDWuvfm11PN6VUnFIquuvfoXR+Yx8xtyrQWj+htU7WWg+n8633x1pr03tISqnwrpvQdA1R3EDn215Taa1tQIlSKr3roQWAqTfIL3AXFhk+6VIMzFJKhXX9bC6g876U6ZRS8V1/pwC3YeDXzaWl9J6glFoDXAPEKqVKgV9qrZ8ztyqgs0d5L5DTNd4M8LOu7QXMNAR4oWuGgB/whtbaMlP2LCgBWNf5M08A8KrWOtPckr7wL8ArXcMVhcADJtcDQNdY7vXAt8yupZvWeqdS6i1gH51DFNlYZ1n920qpQUA78IjW+qxRF7b8NEIhhBA9kyEUIYTwUhLgQgjhpSTAhRDCS0mACyGEl5IAF0IILyUBLoQQXkoCXAghvNT/B3IkptXOzf3HAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "beta2 = fit_and_plot([2, 2.8, 1, 1])\n", "beta2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, now we know we're right, right? The error of this fit is essentially zero." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "def fiterror(beta):\n", " return sum((y - f2(x, *beta))**2)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.8617117363215879e-28" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fiterror(beta2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But wait, what about this:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2. , 9.28318531, 1. , 1. ])" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO29eXikZ3nme7+1LyotVSrtqtLS3erV7m6320sHbGzAYAhunGQCJ4GcTAKThGFgkjFjz5w5meQ6g32GXGQy18kw8cEkJHAMGWMaiAFjbIOD8dbuxb2oV+1V2qqkqpJqX97zx1dfqVSq5dvUpWo9v+vqq1sl6f3e3u56v+d7nvtmnHMQBEEQjYeu3hsgCIIglEECThAE0aCQgBMEQTQoJOAEQRANCgk4QRBEg2K4kRdrb2/nAwMDN/KSBEEQDc/bb78d4Jy7S1+/oQI+MDCAkydP3shLEgRBNDyMsclyr1MJhSAIokEhAScIgmhQSMAJgiAaFBJwgiCIBoUEnCAIokGpKeCMsa8xxhYYY+eLXvsSY+wSY+wdxth3GWOtm7tNgiAIohQpJ/C/A/CBktdeALCfc34LgCsAHtN4XwRBEA3PidM+HHviJQw++hyOPfESTpz2abp+TQHnnL8CYKnktZ9wzjP5D18H0KfprgiCIBqcE6d9+MIz78AXioMD8IXieOzZc5qKuBY18H8J4EeVPskY+zRj7CRj7OTi4qIGlyMIgtj6/JfnRpHK5ta9Fk9n8aXnL2t2DVUCzhj7jwAyAL5Z6Ws4509yzo9wzo+43RsmQQmCIG5KFleTZV/3h+KaXUPxKD1j7HcAfBjA/ZxifQiCIAosrpQXbwDoabVqdh1FAs4Y+wCAfw/gHs55TLPdEARB3AS8fGkBAGA26JDMrJVRrEY9HnlgRLPrSGkjfBrAawBGGGMzjLHfA/D/AHAAeIExdoYx9j812xFBEESD8/pYEE67CU88fAC9rVYwAL2tVjz+8AEcP9Sr2XVqnsA55x8v8/JTmu2AIAjiJuON8SXcMejERw/34aOHN69JjyYxCYIgNGQpmoIvFMfB/s2fbyQBJwiC0JCL/ggAYF9Py6ZfiwScIAhCQy7OhgEAe3uaN/1aJOAEQRAactEfQXeLBU67adOvRQJOEMS2g3OO2XAcmzHCMjq7gj3dm3/6BkjACYLYZnDO8e+/8w7uevwlfP7bZzRdO5fjGA9GMey2a7puJUjACYLYVrw2FsQ/npzBYLsd3zvjxytXtPNo8ofjSGVyGGxv0mzNapCAEwSxrXj6zWm0WI34wWd/Be1NJnzrrSnN1p4ICIPpA+02zdasBgk4QRDbhkw2h5cvLeDBA11oMhvwoQPd+OnoAuKprCbrjwdWAQBDdAInCILQlvP+CFaTGdw93A4AuHd3B1KZHN6eXNZk/fFADFajHp3NZk3WqwUJOEEQ24bXx4IAgDuHXACA2wecMOgYfnk9oMn644FVDLTbwRjTZL1akIATBLFteGMsiB0dTXA7hBNyk9mAvT3NODMd0mT98UAUQ+03pgMFIAEnCGIbccEfwa196z1K9vW04LwvrLonPJ3NYXo5fsMeYAIk4ARBbBOCq0ksrCSxp9ux7vUDvS2IJDKYWVaXlDMXTiCb4+hvIwEnCILQlEtzKwCwYUpyf6/w8XlfWNX6vnxUWm+bdok7tSABJwhiWzA6K7gE7u5afwLf1emAQcdw3q9SwPMn+F4NI9NqQQJOEMSWYTIYxSe/9ib+/AcXkc1p61MyOrsCt8MMV9P6Fj+LUQ+vy4brC1FV64thxVpmXtZCcagxQRCElnDO8e/+11m8NbGMV64sYmdnEz5+1KPZ+uOBVexwlx+wGXI34friqqr1faE42ptMsBj1qtaRA53ACYLYEpyZDuGtiWX8+UP7cGt/K/7nz69r6hY4EYxV7BAZdjdhIhhFJpsr+3kp+ELxG1o+AUjACYLYInz/rB8mgw7HD/Xifzvaj8lgDBfy6TZqCcfTWIqmMOAq36M97LYjneWqOlF8ofgNfYAJkIATBLFFeOXKIu4edqHZYsR793SCMeCno/OarD0ZFOrb3goCPpQvrSgto3DO4Q/F0dNCAk4QxDZjIZLA9cUo7sqPuLuazNjT1Yw3x5c0WX88IAj4YIUpSdG/W6mAB6MpJNI5OoETBLH9EM2kjg46C68dHXTi1NQy0irq0iKTQcHm1esqXwNvtZngspswtqisE6UeLYQACThBEFuAC/4I9Dq2bsjmsLcNiXQOV+fVdYcAwEQgiu4WS9UOkSG3XbGA16OFEJAg4IyxrzHGFhhj54teczLGXmCMXc3/3La52yQI4mbmgj+MHe6mdQK7L5/qfnFW/YPMiWC04gNMEa/LjqmlmKL1xSnMvi1YQvk7AB8oee1RAC9yzncCeDH/MUEQhCIu+CMFwRYZcNlhNepxUYNOlKmlWMXyiYjHacNcJIFEWn64w8xyHHaTHi1Wo9ItKqKmgHPOXwFQ+iThIQBfz//66wCOa7wvgiC2CYsrgsnU3hIB1+sYdnc7cHFW3Yh7Ip1FYDVVsz7tcQoCr6SV0B+Ko6fVesN8wEWU1sA7OeezAJD/uaPSFzLGPs0YO8kYO7m4qF14KEEQN5bXrgfxvTM+5DQecb+6IJhM7e5q3vC5vd3NuOiPqBro8Us0merPC/i0gjJKPXrAgRswSs85fxLAkwBw5MgRbf/mCYK4IfzyegC/9dU3wLnQ0fFv7t+p2dqFFj/3xhr1rk4HIokMFleS6Gi2KFq/4BIo8QSupA7uC8VxsL+19hdqjNIT+DxjrBsA8j8vaLclgiC2Gl/+yRX0tlrxrp3t+B8/u4ZwLK3Z2uOLUZgNOnSXEeihvKiPBZQbTYktfrU6RNqbTLAa9bIFPJrMIBRL3/AOFEC5gH8fwO/kf/07AL6nzXYIgthqjC2u4uTkMj5xpxe3DziRSOdw65//BMeeeAknTvtUrz8RjGKw3Q6dbmP9WJyQVNreBwinYx0Dulqqn+AZY/A4bbIF3F+nDhRAWhvh0wBeAzDCGJthjP0egCcAvI8xdhXA+/IfEwRxE/KLa0Lgr0HH8D9evlZ43ReK47Fnz6kW8bFA5Ra/7mYLLEYdxlQ4BfpCcXQ1W2DU1z6v9jttsmvgMxJLNJtBzRo45/zjFT51v8Z7IQhiC/La9SB6W6146hfjSGTWT0XG01l86fnLOH6oV9HamWwOU8EYHtjXVfbzOh3DgMteqJMrwbcs/QGjx2nDa9cD4JxL7igpTGFuxRM4QRDbl1yO4/WxIO4ccmE2nCj7NWIJQQm+UByZHK/oUQLkJyTVCHi+xU8KHqcV0VQWS9GU5PX9oTgMOoYOh7KHrGogAScIoiJjgSiWY2ncMeSsKIJqHt6N1TCZAoCh9iZMLcWQysj3RMnmOObCCcnlDY9LfieKLxRHV4sF+jI1/M2GBJwgiIqIOZL7eprxyAMjsJZ4iViNejzywIji9ScDos1r5SnJIbcd2RxX1N63sJJAJscllzfERHlZAr4s/YSvNSTgBEFU5NJcBAYdw46OJhw/1IvHHz5QOM0a9QyPP3xAcf0bEE6vZoMO7pKcymJED++pJfllFLkugX1t8od5fKF4XTpQABJwgiCqMDq7gmF3E8wG4eR9/FAvXn30PvzBPcMAgAcPdKtaX4whq/bA0FOYkJRfa5c6xCNiNenR4TAX7GdrkcrkMBdJoI9O4ARBKIFzjnBcu8GaYi7NRrCn27Hh9X09zUhnOa4tqAwCltAhonTABljzNZHTIeJ12TAp8Vpz4QQ4Xzu532hIwAmigeGc418/fRq3/tlP8PiPRjVdeyWRhj+cwK6ujQK+o0NdBJmIL5SoGUPGGEO/06pIwP2hONpsRthM0l1DPE47piSewGdCwtfVo4UQIAEniIbmhYvzeO6dWfS1WfHkK2O4Mr+i2doTAUGchsp0iAy228GYuglJwSUwKUn8PAoGbABlJlNel3RbWfGETzVwgiBk8/SbU+hpseDEZ47BbNDhH16b1Gzt8XwQ8EAZAbcY9ehttao6gftl1KfFCUm5roS+ZflBw2JHjJQ3DN9yHIwB3Tc4zFiEBJwgGpRwPI1Xrgbwqwd70N5kxr27OvD8hTnN7F4nxBY/Z+Ukd3UCLgwGSTkh97fZZA/YiEnxck/g4kNTKQ8yfaE4OhxmmAz1kVIScIJoUN4cX0I2x/GeEcGO/717O7GwksSlOW3KKGKOpNVUPkdyOJ8hqfQNwyfWjyWcwAudKDLCFsLxNKKprGyPErFtUcqDzJnlWN0eYAIk4ATRsLw+FoTZoMMhj+BDfXRASHQ/NbWsyfrjNXIkh9xNiKezmIuUH7GvhW9ZmksgoGxCUml9us1mRJPZgKlg7fq+2AZZL0jACaJBOTMdwi19LYUe7X6nFe1NJpya1EbAJ4MxDLRXPl0O5726lZZRfKEEOiW6BIoiLHfABpA/6i/aytY6gWdzHLOhRN0eYAIk4ATRkORyHKOzEezraSm8xhjDYU+bJifwcDyNpWiq+gm8XWglVOoU6AvFJIurzWRAe5NZnoDLnMIsxuuy1WwlnI/IG9PfDEjACaIBmQhGEUtlsbd7fY7kQU8rJoIxhGLSH/aVXT9QuQNFpMNhhtmgU9TeB8gvP3hk9oL7Q3FYjDo47SbZe/O4bJhejiFbpb4vd8pzMyABJ4gG5IJfMJkqTXLfkw8GvjKvbsBGrB/3V3lAp9Mx9CtIsAGKXAJlnF77nYKoSkXKmH4lvE470lmO2XDlh6Zil4r4gLUekIATRANycVYwmdrVuX5KciQ/NXlZ5UCPT+KEoRBBJt+jZHEliXSWyzyB2+APJZDOSrOVleMDXorYC16tjDIRiEKffxOrFyTgBLFJTASi+NTfn8Tf/Py65mtfX1jFQLt9Q/9xd4sFDosBl+ciqtb3LcfhMBvQYjVW/TqP0gEbGS2EIv1OW+HBoaRrLCt3CSz0gle5uxgPRtHXZpX0EHazIAEniE0gl+P4w2+ewgsX5/H4jy7h5csLmq4vBgGXwhjDSKcDV+bUe5RIGrBx2rCazGBZZkq9T8YQT+FaotWrhDJKIp1FMJpSXJ/uabXCqGdVh3kmqmR53ihIwAliE3j1egCjsxH811+7BR6nDV95WbtTeDbHMRGMlfUoAYBdXQ5cnl+RfSouRmr5QTypyq2Dix0ickoccnrBCw8YFZ7A9TqGvjZbRQ9yzjkmAuXfRG8kJOAEsQmcOO1Hi9WIhw714Ndv68ObE0uqsiOL8YfiSGVyFcVjqN2OcDwt+1RcjG85JmtCUraAh2JosQoDM1IRkuWZNAEX3yBUeJQMu+0V7XIXV5OIprIYqJIkdCMgAScIjeGc47XrAdw97ILZoMeDB4TE9Z9fWdRk/fEaLX7ibf2EhEnCcqwk0ogkMhJLKPIHbIC8D7jM8oZex9DbKq2VUO0JHAB2dTowthgtm8U5nndhHHQ3KV5fC0jACUJjppZi8IcTuHvYBQAYdjfB7TDjtetBTdYXhblSCUUU9gmFAzaiyZSU8oY4YCPVP7v4GkrE1euyY1LKiPtyHHodQ1ez8qT4kS4HMjledlBJ/DsYpBo4QdxciEJ9V17AGWO4a8iFtyaWNFl/bDEKu0kPt6N8jqTHaYOOKRdwuR0icgdsOOeKPUQGXDZMBmp3vfhDcXQ1W2BQ0SEitmiWa8kcW4zCpNehp1X5G4QWkIAThMacnQmj2WLAcNHt9S19LZgNJxBYTapefzwQxaDbXnFAxWTQobfNinGZp2IRn0wTKI/MYZ5IPIPVZEaZgLfbsZLM1LSVnQnFVYvrkNsOvY7hahkBvzgbwc7OJlVvEFqg6uqMsX/LGLvAGDvPGHuaMVbftyOC2AKMzkawp7t5ncCKniXnfWHV608vxWpO/w247IpP4DOhOIx6VjUpvhiP04bZcFzygI2aGDKp9f2ZJfU2r2aDHoPtdlwuY887OhvZYGNQDxQLOGOsF8C/AXCEc74fgB7Ax7TaGEE0Irkcx+W5Fewp+c8tjryLI/BKkVp+GGwXBFxJK6E/lEB3ixU6nbQR9H6nDTm+dnKvhdwTfjHihKQY91aORDoLfzihSY/2SKdjQwllYSWBwGpqw99xPVB7/jcAsDLGDABsAPzqt0QQjcvkUgzx9EaTqRarEV6XTfUJPBhNIZnJ1RTwAZdQagjKSLARkdpCKCK3lVCNCVRfm1Dfr/YgU+yIqWaFK5X9vS2YDMYQLCp9Xcy/CTe0gHPOfQD+AsAUgFkAYc75T0q/jjH2acbYScbYycVFbdqoCEItb44v4ZUri6qGXcpxabbyf+49Xc1lb8flULBIrVEeEMVLSsfGhmvIjCGTG7Yws6zcJVBKfX8i/zmvBifw2wfaAABvF3msn5sJgzE0fAmlDcBDAAYB9ACwM8Z+u/TrOOdPcs6PcM6PuN1u5TslCI14+dIC/sXfvIZPfu1NfPutaU3XFgc/dnRs7A8e7rBjaikmuVZcDqmnV08+x1JKrmMxqUwOCytJWafjTocFJoNO+gl8OY6+Npsil0BAuLuo9sYkfk6LIZsDfS0wGXQ4WSTgb04sYaTTgRZbdZ+YG4GaEsp7AYxzzhc552kAzwK4W5ttEcTmwDnHF384iiG3Hbf0teAvf3oFyUxWs/XHg1F0NZfPkRx2NyGT47JFtRipSe79TisYky/gc+EEOJdX3tDpGPrbrJJ7wdXGkA247BivUt+fCEbRajOi1Sb/hF+K2aDHrX0teGNMaA1NZrJ4e3IZdww6Va+tBWoEfArAnYwxGxPeSu8HMKrNtghiczjvi+Dqwio+9a4hfO7+nZiPJPHqtYBm608EohVrr0P5tsIxFUnuM8txNJkNaLZWH0E3G/TobrbIHnFX2iEip5VwZjmmakLS67JhJZFBqIJVwGQwpkn5ROTekQ6cnQnDH4rjF1cDiKWyuDcfJF1v1NTA3wDwDIBTAM7l13pSo30RxKbwwsU56HUMD+7vxq/sbIfDbMDz5+c1W38iGKvsUVLIkFTW3gfICynwuGyya+BKTKYAod48JcFWNpYSnAvV5EiKf47XKrwRjgeimnqUfOhANwDg229N49tvTaPZYsCxHe2ara8GVV0onPM/5Zzv5pzv55x/gnOufkqBIDaRX14P4kBvC1psRpgNety9w4XXxrQZca+VI9lsMcLtMKs6gfuWpT9g9DrtCkymBAHvlpAUX4xUW1k1OZUi4oTklTIDNslMFv5QHF4NQxYG2u34wL4u/NWLV/GTi/P4vV8Z2uDDXi+2xi4I4gaQSGdxZjqEO4dchdeOeJ2YWophcUX92aPw8KyKxeiw2644xR0A/GHpE4Yelw2B1RRWkxnp64ficDvMsBg31vCrXsspretlJqS8B1ykt9UKu0mPK2U6eq4vRJHjwM6SpCK1/NlD+/CeETcePtSLf3XPkKZrq4EEnNg2XJpbQSbHcbB/Lcn9sLcVADRJchdNj6p5RA+5m3B9UdmATTQp1H17W6WdLqXEgpWi9AGjV2Ir4UxhiEf5CZkxVvA8L+XyvNDGubtLWwHvbLbgb3/3KL78mwdlv7ltJiTgxLbhgl8YohHH2sVfm/Q6nJpUL+ATgRgYqx5yK3p1h+PyvbrlWqR6862ElUIJyl5DRommmEJaTg0B9y3HYdLrJI/pV2Kk04GrZYKbL82uwKTXVb0LupkgASe2DRf8ETRbDOtu3398fg4cHH/zyhiOPfESTpz2KV5/ckloIax2QutXGIAAFNePpZdQAOmthLkchz+cUHQCt5r06HCYa15rZjmG7laL5DH9SuzqdCAYTWEhsj4f89LcCnZ0NNU1p/JGsj1+lwQBQcD39qyZTJ047cNjz55DOiuUM3yhOB579pxiEfctxwsn0UqIp/NpBUnua0M80soPLVYjWm3GqsG8xQSiSaQkjOlXQkor4UQwqkmL3639YukrBED4uzz2xEv4+ZVFjAeiqt6IGwkScGLLoWZSsRKCyVRk3Yj7l56/jHh6/RBPPJ3Fl56/rOgaPgkWpqpO4HmXwI4KPuDl8Dpt0gdsVHaIeFzVBVzIkYxhUIMWv/29zTDpdTg9tVx4Ixbf4OLprKo34kaCBJzYMuRyHJ/55ins/k8/xj+8NqHp2nORBBLp3DqP7koZlUqyK7M5jrlw7ZSZJrMBLrtJcQlFjksgAHhcdkxKrIHLSeIpey2nLf/nXH6yVeyI0aI+bTbosbenGaemljV/I24kSMCJLcN3T/vw3LlZtFqN+C8/HMVcOFH7myQidogUx5BVEiolAjYfSSCT45LKG/1Om+wMSUDaCb8Ur9MGfygh6a7Gp8KnGxCsAjgX0mrKMSGhzVIOt3nb8M5MuHDyLkWrEOmtDAk4sWX4xhuT2NHRhGf/6G4k0jl859SMZmsXWvzca+LxyAMjsJY8cLQa9XjkgRHZ6/tldIj0y0ywKb6G1Pq3iMdlQzbHJXl1+5bjcJgNaLEqM2na2Snc3VSbkAS0y5G8Z5cbyUyuoquh0juJRoIEnNgSzEcSOD0VwsOHe+F12XHE24YfnNXOXn48EIXVqEenY+0Ee/xQLx5/+ECh5ms26PD4wwdw/FCv7PXXHjDWPiF7nFb4QnFkZNT609kc5iPyg4DFiUQpDzLl2siWMthuh44B18r0ZwOCT4xBx1QN8RRz55ALDrMBXpcNppKuE6VvxI0GCTixJXjliuAVf+8uwSTovXs7cWluRZMJSSDvj9Fu31A/Pn6oF68+eh+OH+xBe5NZkXgDawMqUk59HqdwKp6VUSKaCyeQ40CfAo8SAJiSkuQeSqg6tZoNenhd9oon8IlgFB6nTbMcSZNBh48d7cfZ6RBcTSaI9jC9rVbFb8SNBgk4sSV4Y3wJTrupMEEnjru/rpFPyXgguq7+XcpAux3+cLziA7ha+EJxOO0m2EzVXQIBZZ0o4glfrsB2OMwwG3SSesHlJvGUY9jdVPBEL2U8ENN8wOYP7hlGT6sVc5EEvvjRA5h44kPCG/I2EG+ABJzYIpybCeOWvpbCCXl/TzOsRv26JBSlpLM5TC1VdgkEhNt/zpW19wFCfVqyR4kSAV+WN4UpotMxeJy2miWUlUQakURGVQkFEOrg44HohvJQJpvD9cVVDLu1FXBXkxk//eN78Ppj9+PjRz2art0IkIATdSeeyuLqwgpu6V0bcTfoddjb01zIH1TDzHIc2RyvevoTHQTHFSa5+5ale4h0t1hh0DFFJ3C5LoGA4FNSqxdc6Qm/lB3uJqSzvBBrJjIWiCKVyRXCnbXEYtSjs1n+n8vNAAk4UXcuzoaR40KAbDH7e5pxwR9GLqcut3JmWRCT/iqnS1HcJxQI+FpSvLQOEX3+QZ6cVkLfsjKXQECIV6vl1S0KvFob1t3dQglsdHb9G+9olaxQQjkk4ETdOTcjmEwd6Fsv4Pt6WxBNZQv9w0rxSzhdtliNcNpNiq4ViqURS2VllR/kthKqiSHzumyIp7NVHwiLNfJKXuZS2dXpgMWow9np0LrXL/ojMOl16wapCPWQgBN15+rCKpotBnSV3Abvy99un1dZRvEtx6FjQFeN8oPXZVNUQpHTQigiJ4IMUBdDVjC1qnK98XyOpNqgXqNeh/09LTg7UyLgs5FtZTJ1o6A/TUISV+ZX8NmnT2+Kv8TYYhRD7qYNMWE7OxxV+4qlMhOKo6vZUlM8Blx2Wd7ZInJNpgBBwEMxabayuRyHP5RQ3D9d6AWv8nub1MhkChCMps75woUHmdkcx5npEG7tb6nxnYRcSMCJmiQzWfzhN97GD8768flvn8F5X1jT9ccCq4Wcw2JMBh08ThuuK3ywKCLV49rjtGE2kpCdUq+kQ0QMQJBSBw+sJpHK5mT3gIv0tdmgY9V7wScCMc1yJA/2tyKRzhXunC7NRbCSyOCOQVeN7yTkQgJO1OSFi/O4vhjFl//FrWi2GPCVn13XbO3VZAbzkWTF2uiQuwnXK/QVS0WIIastfl6XDZzLt3r1heKwGvVok1F+6JdwKhaZVthCKGIy6NDdYq1YQklmsvCH45qdwI/taAdjwM8uLwAAXr0WAAAcHXRqsj6xBgk4UZMTp/3obDbjoYO9+OihXrx4aR7xlLKBl1LG88ZHlfqDh912jAeiijtRsjmO2ZC0kIK1WDD5Se49rRZJSfEicnrBlZRoSvG6bBXfLKaX4uAcmp3AnXYTDvW3CmEZnOO5c3PY39u8LbxJbjQk4ERV0tkcfnk9gPfv7YJex/DAvi4k0jn889VFTdYfCwin66EqJ/BkJlfRca4WCyt5l0BJJRThTURqgo2IPxxHr8yMR4dF6HqR8mahdIinGG8Vr24pYcxy+Y0j/bg0t4K/fOEKzk6H8JFbezRbm1iDBJyoyjszIcRSWRzbIdQvbxtog8mgw1sTS5qsf30xCh1bO/2WIo6/j6kYsAGkhRS0N5lgM+llC7icIZ5ipHaizCzH0Gozoslce0y/8rXsWIqmyj40Fe1f1bYQFnP8YC+G2u347y9dQ2ezGb99p1eztYk1SMCJqrw1IYyyH80/gDIb9DjQ21KIslLLeCCKvjYbzIbyAyrDHcLJXGkdfK38UFtgGWOy2/viqSyC0ZSiDhGPs3JZoxg1PeAiO/N/jlfLdPRcmluB22GuaMuqBKtJj3/4/TvwyAMjePpTd0ryiCHko0rAGWOtjLFnGGOXGGOjjLG7tNoYsTW44I+gt9W67j/3YU8rzs2EZXdrlGNqKVbx9A0ALrsJzRZDodQiF9lJ7i5boaQgZ325QQvitfyheM2wBaUn/GL25HvqSyckAeDyfKRgIqYlva1WfOY9OyqWxwj1qD2B/xWAH3POdwO4FcCo+i0RW4kL/vAG/4rbvG1IZXO4oIFPiW85VvX0yhjDYLtddlljbf042mxGySdAr8uO6eW45IemfhUPGPudNuQ4qoYtcM4xsxxHn8waeyk9LRa0WI24OLv+BJ7J5nB1fhUjndoLOLH5KBZwxlgzgHcDeAoAOOcpzrk299XEliCazGA8EMX+npIR9/zHl2bVDdjEU1kEVlM1xUlpBBkgP6TA47QhlclhLiLNq1vuCb8Yr4ROlGA0hXg6qzoEgTGGvd3NuFhyAp8IxpDM5DCyCSdwYvNRcxazBQwAACAASURBVAIfArAI4G8ZY6cZY19ljG14CsIY+zRj7CRj7OTiojadC8SNYXQ2As7XRtpFelutsJv0uKJyQlLMYKwlTh6nreAoKPsay3H0tMgfsJF64vctx6HXMXTKSIpfu1a+66WKgIvmWtWscKWyp7sZl+ci66xeyWSqsVEj4AYAhwF8hXN+CEAUwKOlX8Q5f5JzfoRzfsTtdqu4HFGJtyeX8c6M9jc/V/MPDktPZzodw64uBy7PqRNwcUClloD3O23I5Dhmw/JaCTnnQo6kjNOr2IkhtRfclx/TV5Iy0+Eww2TQVb27EL1ZtGjxO+xdPyEJAKemlmEx6ugE3qCoEfAZADOc8zfyHz8DQdCJG8hz78zi177ySzz016/i51e0vcMZD0RhMujKDmCMdDpUn8BnCgJevYQiDr3InZAMx9OIprKyHgB2t1hg0DHpJ3AVHSI6HYPXWd1AayIYLdjPqkUcZS9OOTo1uYxb+lrJZKpBUfy3xjmfAzDNGBOTQ+8HcFGTXRGSyOY4vvjDUezucqC31YovPjda1fNZLmOLUQy4bNDrNk4Y7up0IBhNqcqsnFmOwaTXwd1UvfywJuDy6uAzEk/4xRj0OvS1VR47L0WcwlTKcA2rgIlADP1tVk0E1u0wY0dHU0HAw7E0zvsjODpAI+6Nitp/FZ8F8E3G2DsADgL4ovotEVJ59VoAvlAcn71vJz7znh24PL+iSWeIyHhgtWLtdWdnvj+7QoCtFGbyJlOlQcOldLdYoJeZYAMoT5nxSHQlzGSFh51qJiR3djZhcilWsSVTDGPWimPDLrwxtoTVZAY/u7KAbI7jvj0dmq1P3FhUCTjn/Ey+vn0L5/w451x9gCEhmZcuLcBi1OH+PR344P4uGHQMPzo/q8namUKOZPkeXrFWrCTBRkRoj6stfga9Dr2tVtkC7pcxxFOM1ymtF3xhJYlsjqvyKNnR0YRsjmMisPH3xjnHRDCq6YTkRw72IJ7O4gdn/fjWm9PoarbgYF+rZusTNxYqfDUwr1xdxB2DLliMerTaTLi1vxWvXdcmxd0XiiOd5RWT3HtarTDpdRhXkZZTqwe8GLkTksL6cViMOtkThl6XDZFEBqFYqvr6KloIRXbkJyRLk9xPnPbhrsdfQiyVxXdP+zTzYT/sacOtfS34P06cx2tjQfz+uwZr3gERWxcS8AYlFEthbDGKO4bW6pd3DjnxzkwY0WRG9fqi98hgBZdAvY6h32lVfAKX2gMuoqQXXHzAKMclEFirudd6kLnms6KuBs4YcHVh7YHwidM+PPbsuUIvejiexmPPntNExBlj+PJvHsSB3hY8fKgXn7xrQPWaRP0gAW9QzuVDFW4tuv094nUik+Oa1MGl9B8PttvL3vpLQWoPuIjHaUMwmsKqjDcnX0iaD3gpUvqzxfUBdUnuFqMe/W22dSfwLz1/GfH0+pp4PJ3Fl56/rPg6xQy7m3DiM8fw5d88CJOBJKCRob+9BkUU8OIpyUKGpAaJOZPBGJrMBriqlB8GXHZMBJV5dc/IcAkElHWi+CTW2Ctdq1qCDSB00TjtJtVGTTs6mnB1fk3A/RWscyu9TmxfSMAblHMzYXhdtnUhtB3NFrgdZk1O4DPLtcsPA+12JGWMnRfjDwnfI7V+LCcAAVhzCVTSo2016dHhMNcsoUwGqxtxSWVvdzOuLa4WQjIqnegpEIEohQS8Qbngj2zwKAGEU/gFv/oTuJQJRrG8oqQO7g8JI+gdDmn1Y7kncLG8odQEyuuy1SyhTAZjmnSI3NrfimyO43z+7+2RB0ZgKSltWI16PPLASLlvJ7YxJOANSCKdxfRyrOCVXcz+nhZcXVhVbfUq1I+ri6vYn6ykE0UcQS83JFSOFpsRzRaD5BO42g4Rj9NetZUwkRZzJNWfwMW09rPTgh3C8UO9+Nx7dxY+39tqxeMPH8DxQ72qr0XcXJDL+iaRyebwB994G6+PLeEvfuMWfGB/t2ZrTwZj4Lx8juTOTqGveDIYwy6FFqGryQzC8XTN/uauZgtMep2koZdSlIyge6rkOpYysyx8ndIxd6/Lhu+cSiKRzsJi3Bg2Mb0Uy+dIqj+Bdzgs6G214vT0mp9Ns1Uojf3s392r6SAPcXNBJ/BN4ptvTOGnowsw6Bn+w3fPI5ZS39onMpaffiyX5C6+pibJXWoGo17H0Oe0KvLq9ks44Zfiddqll1CW4zDoGDqblbX4rQUcl7/eRP73rMUJHBA81t8YWyo8EH5zfAluh1mz9YmbExLwTYBzjm+8PolDnlb8v588gqVoCt8749ds/bEqLX6DKjMkgbUWP0lJ7s7ateJSsjmOubD8EfR+pw3TyzFJtrK+UBzdrdJLNKXU6gUvBAFrNCV574gbgdUkLvgjSGdz+PmVRRwbdsnuYSe2FyTgm8CV+VVcXVjFw4f7cMTbBo/ThucvzGm2/vWFVXQ1W2AvE3JrNxvQ1WxR5VHiEztEpAi4y46pYFSWidbiShKZHJfdVeF12ZDOckldL2pjyAq94BXq4BPBKJotBrQWdQGp4Z5dbhh0DN897cOr1wIIxdKalt2ImxMS8E3gF9cCAID7dneAMYb37e3EL68FkUirz5AEgOuBKIYqTEieOO3DciyFZ0/5cOyJlxRN7/mW4zDqGTokhBR4XTZE8y17ktfPn/Blm0wVTsW17y6ENkjl5Yc2mxEOc+WHppPBGAba7ZqdkF1NZnzwQDe+/dYU/vT7F9DhMOPeEfLPJ6pDAr4JvHY9iAGXrXACPDroRCqb02TAhnOOscXVsgIujmAnM0Liii8UVzSC7QvF0d1S2yUQkJ9gI6wv/YRfjNRWwlQmh/kVdS6BjDF42yt7dV+dX63oE6OULzwwAqtJj8lgDP/xQ3vKPjwliGJIwDeBszMhHPa2FT4+7BF+fWpKvVnjUjSFlUSmrEugViPYvuWYjAlJeQk2wNpEYXeLvAeMUsMW5sIJcA70qRx82dVZPnUoHEtjLpLASJe2MWT9Tht++sf34J+/8B48dJBaBonakIBrzEIkgcWV5LohG7fDjH6nFacm1ceeiSPo4mm0GK1GsOV4iPQ7rWBM3gncH4qj2WKAwyKvfiw1bEFsIVSbYrO3uxkLK0kEVteHVlyaEyZdd3drH0PWajOhv8zfLUGUgwRcY8Qx9tIg4MOeNpyZ1k7Ay4mTFiPYqUwOCytJyeUHs0GP7maLrF5wYcpTmUhJcSWc0cDmFRAEHFgL/hW5nI+S2005kkSdIQHXGHGMfW+JgO/tbsZcJFHTY7oW0+KAShlxeuSBEVhL6qZyR7CVlB88EsbOi/GFEootWL0ShnmmgjHodUy1d8ieCgJ+aW4FzRah24cg6gkJuMZcnI1gwGXbUB7YlT+tXZlX3t4HCOWBFqsRzWXKD8cP9eLxhw8U6tdWo172CPZMqPIbRCW8Tru8h5jLMcXi6nHaEI6nEY6lK37NeDCKPg1yJNvsJnS3WPDOzPqHzxd8YezpbqYebaLubEsBvzq/gi88cxb/fFXbFHdACALe0bHx1nokP9Z+WYMk92q13eOHevHqo/fhriEX9vY0y/bPEKcw5Qisx2VDYDUpKUhiJZFGJJFRIeDiQ9PKbxgTAe1iyG4fcOLN8aVCn3s8lcUFfwS3FT2kJoh6se0EPJbK4H//27fwjydn8Km/P4nZsHYey7kcx3iFHu3uFgscZgOulOlqkIPUHEklEWTAms2rnA6RWmPnxcyGhfWVCvhAu3CtsUD5OxnOOSYC0apBFHK4c8iFhZVkYbL17EwImRwnASe2BNtOwL93xg9fKI7/+uu3IJcD/ubnY5qt7Q/HkczkyooHYwy7usq3pUmFc46Z5Rj6JTwA9LhsWFxJFjympeILxeB2mGX1IHud4tRibQEvuAQqrIEPttuhY5W9XhZXk4imshjQyEPkrmEXABSyRn92eREGHcORAWe1byOIG8K2E/DvnvJhV2cTfuO2Pty3uwPPnZuV5K0hhfEaMWTDbrsqj5JgNIVEOifpBC62ookPPaWi1CUQkNYLvpYUr0xgzQY9vC47rlYQcDHiTSsHvwGXDf1OK358fg6cc7xwcQ63DzjRYtVmhJ4g1LCtBDySSOPtqWW8b28nGGN48JZuLK4kcXZGfXsfsCbglSb0BtrtCKwmsZKo/ACuGmsthLXFrz8v8nKtXn3LtYMcSmmxGtFqM0o7geddAt0SxvQrsaOjaUOKu8hEsHaWpxwYY/i1w3149XoAT/1iHNcXozh+qEeTtQlCLdtKwN8cW0I2x/GunYLHxN352+M3xpY0WX9sMQq7SV9RnAZd0ksN5SgMqDil1cABeSfwXI7DH0ooMoHySqy5+0JxdLUodwkEBAEfD0SRzuY2fG4iEIVBx1QZWZXyiTu9cNlN+L+eG0Vvq5WCFYgtw7YS8LMzIeh1rJDk3t5kxs6OJrwxHtRk/fFAFIPuygZHhQQbhWUUOUHATrsJdpNe1oPMQDSJVDanSPw8LmmthNNLsbJTpHLY2dGETD60opRrC6vwuGwwqGwhLMbVZMbX/+VR/NG9w/jm798Bs4E8Soitgep/5YwxPWPsNGPsn7TY0GbyzkwYOzuaYDWt/Qc82N+KczNhWXaolRiv0b62ZvykVMCFHnApI+iMMUlTi8X4ZCbFF+N12uALxcueiouZ0kDAd+Sj5K4tbHwgPDoXwR6NPUoAYF9PC77wgd2UjkNsKbQ4pnwOwKgG62wqnHOc84VxS9/6IOB9Pc0IRlOYjyQrfKc0sjkOfyheVZxsJgM6m80YDygrofhDCVn+HnJbCdXkSHpctsKfQTlOnPbhrsdfRGA1hR+em1Vkcysipg5dnltfB19JpDG9FMeeTfAoIYitiCoBZ4z1AfgQgK9qs53NwxeKYymawoF8+URkX68g6GqT3OcjCWRyvKb4DbjshQdtcvHLMJkCRN+QuOS7C6lRauXwVkmwEW1uxR7wSCKjyOZWxG42YNhtxzslD58v5Vs0xRF4grjZUXsC/28AvgCg4n0zY+zTjLGTjLGTi4vaTz5KRey/3ltyOhNGooHzvki5b5PMWntcdfEbbLdjQmENXG7KjMdpQzydRWBVmv+KLxSHw2IoO6Zfi0KCTZkTv1Y2t8Xc2t+KszOhdW9OomcJCTixXVAs4IyxDwNY4Jy/Xe3rOOdPcs6PcM6PuN31SxgZWxREszQIuMlswIDLjouz6k7gPokCPtBuRzCaQkRmK2EkkcZKMiMrCFgs50gtowhTnsrq0x0OM8wGHabK3F1oZXNbzKH+VgRWU4UHuwBwdjoMZ96/hCC2A2pO4McAfIQxNgHgWwDuY4x9Q5NdbQJjgVU47Sa02kwbPrejo6kg8EqZkVh+EEsNcvuzRbGTW0IBpIctqMmR1OkYPM7yToFa2NyWcrBfGGU/XWTR+9bEEo5428hkitg2KBZwzvljnPM+zvkAgI8BeIlz/tua7Uxjri9EMVwhR3LILdSlMzU6KKrhC8XRZjPCZtoYNFyMR4ZvSDFKBLyvTQhbmArWPulyzuELSfNZqUQlq1ctbG5L2d3tgMNiwKtXhfxRfyiOqaUYjg7SiDuxfdg2feBjgVUMlYkhA4Dh9iaks3zd7bhchJACaSnugPxhHjFHUo5Pt8UohC1IaVuMxDNYTWZUDcB4nHZMLcU2PDQVbW7F4Z3eVqtsm9tSjHod3r3TjZcvLyCb4/jJhTkAwL0jHYrXJIhGQxMB55z/jHP+YS3W2gzC8TQCq6mKSe7DHcLrlRzupCC1/NBkNsBlN8nKkASENwijnqG9Sd4IutSwBdEHXO0JPJ7OYnFlY0vmBw90gXOOz963A68+ep8m04wPHujGwkoSL19awLOnfdjZ0VToESeI7cC2OIGPLQrCPOQu/59bPJkrrYOL5Qep5Q2PhFSZUnzL0pPii5EatqCmhVBELA+Ve8OYCMSQ49BUYN+/rxO9rVb8wTfexjszYXzqXUOarU0QjcA2EfC8yVSFE3ib3YQ2mxHXF5WdwEOxNGKprOTyg7fCw75qCD3g8rsrvO1C2MJqjbAFOWP6lRC9XsbLvBFezU9NaingRr0Of/1bh+F12fDwoV48fJg8SojtRfUnbjcJk0sxMIaqPtrD7iZcV3gCF1sIpZYfPC47vn/Wj2QmK9lXwx+K4868+ZYcRK/uqWBsQ05nMb5QHFajHk77xi4dqfQ7bbCZ9Lg4u7Gn/trCKhjb2MaploP9rXjxT+7VdE2CaBS2xQncH4qj02GByVD5tzvYbldsMuWT6XHtddqQ45D80DSTzWEuotAlUKL/imgjq6YFT69jGOlylBXwqwur6G+zyQqKIAiiOttCwKV4XHsVJtiI6wOQXOIoRJBJLKPMrySR48rKG9Xq0sXMhGKaWLDu7W7G6GxkQyfK6GwEuzrJo4QgtGR7CLiElBlx6GVGZoKNuL7FqJNcfvDIdCVUEjQs0mwxwmk3Vb2WkCMZ0yQEYW9PM1YSmXV3F+F4GmOLURzsb6nynQRByOWmF/BcjmM2XPsE3i9z7LwYsYVQavnB3WSGzaSX1N4HKBviKabSgI1IYDWF1WRGkxzJ/T2CSJ8pmpA8NyPYFNza31r2ewiCUMZNL+ALK0mks7zmCVyub0gx/nAcvTI8RBgTxs6lllB8IXklmlJqdb2I7oheDU7g+3qa0WQ24PWxtZCMM9PLAIBb+kjACUJLbnoB9+UHVGoJuMtugk1mgk3hGgo8RLwSB2wA4QQuZUy/Eh6XHf5wHMlM+fp+IYy5ShiFVAx6HW4faMNrRQJ+cnIZQ247BQEThMbc9AIu1WRKPBVPL8kbp4+nsghGU+iVeTr2uoSx81yutle3XB/wUgZcNnCOir83MUdSzRRmMXcNuzC2GMVsOI54KovXrgfx7p31c6IkiJuVm17Apdq8AkLau5wIsnXryxQ/j9OGVCaH+ZVEza9VGjQsIk6gVkpynwzG0O/ULkfy/j2dAIDvnfHjp6PzSGZyeG/+NYIgtOOmF3B/KI5WmxF2c+3ygxhBJicf0y+zB1xkrT+7+huG3DH9cuzqFAT8yvzGDElAzPJU/wBTZNjdhKMDTjz1i3H85QtXMNhux10KhpAIgqjOTS/gcurTHqdVVoINoPwEXjwhWY1wPK3aJdBmMsDjtOFyGQHnnGMyGC24JGrF//mrexGOpzEejOILD4wUnAgJgtCOm36U3heKV02KL0ZsJZxejsHtkOb651uOQ69j6JT49SI9rRYYdAyTNVwJxYeqHpUn5F2djkKsXDEzy3FEU1ns7NR2xH1/bwte+pN7kMzkNB+fJwhC4KY+gXPOJU1hioithHLq4L5QHF3NFtn1Y4Neh942a80SSkHAq6TdS2F3lwPjgeiGTpTNzJHsa7OReBPEJnJTC3g4nkZUhkugmAcpJ+5MTQyZWHOvhlYCvqvLgWyOb7DMHZ1dAWPACI25E0TDcVMLuFyLVKtJjw6HWVYvuJoYsloTkoBwN9DeZJL0ELYau7sEgS4to4zORuB12lSvTxDEjeemFnAlDxj7nTZMS/BDOXHah7sffxG+UBwvXJzHidM+2fvzOu0Ix9MIxyon1IstfmoZarfDZtLj9NTyutdH5yKbUj4hCGLzubkFXEFIgZRhnhOnfXjs2XPwh4Ue7pVkBo89e062iK85BVZ+kDm1FFNdPgGEmvthTxvenFgT8MBqEpPBGI24E0SDclMLuF+mSyAgnMCrjZ0DwJeev4x4ev3n4+ksvvT8ZVn7q9ULns7m4A/FNRFwALh9wIlLcxGE48KJ/+TEEgDg6GCbJusTBHFjaQgBP3Hah2NPvITBR5/DsSdeknzSFW1k5YQUeJ3C2Hm1sAVxeEfq65WoZaDlD8WR49CkhAIAtw+2gXPgzXFBuF+9FoTFqMP+XrJ5JYhGZMsLuFiu8IXi4BBEWWq5wheS5xIISAtbqDQVKXda0mYywO0wY6JCEpBoMiW1j70Wt3nb0Gwx4EfnZpHLcfz4whzu3dUhOdaNIIitxZYXcDXlCmUugYJYVgtAeOSBEVhLosGsRj0eeWBE1rWAvNVrhRO46F2iVRCw2aDHgwe68aPzc/i7X05gcSWJD9/arcnaBEHceLa8gCstVyh1CWxvEmxlJ6qcwI8f6sXjDx+AOZ+x2dtqxeMPH8DxQ/JT0T2uyr7g1xZW4bKbVAUNl/JH9+5AlnP8+T9dxK7OJnxwPwk4QTQqigWcMdbPGHuZMTbKGLvAGPuclhsTUVquUOpRwhgrWL1W4/ihXnS3WPChW7rx6qP3KRJvQDB+mosksJLY2Ep4dWEVwxqdvkU8Lhu+9ju343ePDeCrn7ydPEoIooFRcwLPAPgTzvkeAHcC+AxjbK8221pDablCblJ8MV6nrZBSU4lcjsMfSqBPZRDwnm5hwOZSyYAN5xxX51ewU2MBB4Bf2dmOP/3Vfar9VQiCqC+KBZxzPss5P5X/9QqAUQDKjqFVEMsVJpnlCp/EIIdyeNttmFmKI1slbCGwmkQqm1O0fjF7u4UOkIv+yLrXF1eSiCQymtW/CYK4+dBkfpoxNgDgEIA3ynzu0wA+DQAej0fR+scP9eK160G8eGkBrz56n6Tv8YeUuQQCwoRkKpvDXKRykIJYYumX2eVSSmezGU67aYOAn/cLQcB7aUqSIIgKqH6IyRhrAvAdAJ/nnEdKP885f5JzfoRzfsTtVh6r5XHZEFhNIpbKSPp6pS6BAArhBpMV2vuAohY/lUHAjDHs7W7GhdnwutfPTIWgY8CBPurRJgiiPKoEnDFmhCDe3+ScP6vNlsojNzVejo3shmsVRtyrJ7nrNcqRvKWvBZdmVxBNrr05nZkJY1enQ3GQMUEQNz9qulAYgKcAjHLOv6zdlspTEHCJVq8zyzHF4trdYoVJr6v6IHMiEEN/mxVGDXIk7xp2IZPjeCs/2s45x9npEA72k0cJQRCVUaM+xwB8AsB9jLEz+R8ParSvDRQmJCWcwFOZHGYjiYK/t1z0OoY+p7Xqm8V4QLsYsiNeJ4x6htfGggAEj+5wPI3DXvIoIQiiMmq6UH7BOWec81s45wfzP36o5eaKabEa4bAYJAn4bDgOzoF+FeWNAZe94jAP5xwTwSgGVda/RawmPQ71t+HnlxcBAC9fXgAA3Dui/JkBQRA3P1t+ElOEMSYpwQZAwQ5W6QkcEMbXry+uIpPNbfjc4moSsVRW0yT3D93SjUtzKzg7HcIzb8/gsKcVHQ55U6QEQWwvGkbAAaGMIkXAZ/KBDP1O5SfwXZ0OpDK5sqfwiYDwmtoOlGKOH+xFs8WA3/rqGxgPRPG7xwY1W5sgiJuThhLwfmftARtASJXX6xi6mpWfYMUIsivzG5PcxxYFkymtSigA0GIz4i9+41ZYjDr8+m19+NAB8ighCKI6DdWj5nHakMrmMB9JVPVCmVmOo7tFWQ+4yI6OJuiYMOL+YImYXppbgdWoVz3EU8r793Xh/fu6NF2TIIibl4Y6gXudwom3VhlleimmWlwtRj0GXHZcmdt4Ah+djWCkywEdGUERBFFHGkrApfaCzywrT4ovZlenA5fnN5pMXZpboSBggiDqTkMJeE+rBXodq3oCT6SzWFhJahJDNtLlwEQwum5Cci6SQDiext68iyBBEES9aCgBN+h16G21VhVw0UZWixP4QU8rOAfOTIcKr533CXYvu+kEThBEnWkoAQdQsxd8Ov85NT3gIrd528DYWggwICS5G/UMBygImCCIOtN4Al6jF3zNJVC9gDdbjNjb3VzwKAGANyeWcEtfKyxGCgImCKK+NJyAD7rsWIqmEIqlyn5+PBCFw2yAu0m+D3g57hh04e3JZUSTGQRXkzg7HcLdwy5N1iYIglBDwwn4jk4hoebK/GrZz48tRjHotkMwS1TPBw90IZnJ4YWL83j+wjxyHHiAerUJgtgCNNQgDyC09gHChOTRQeeGz48trpZ9XSm3edow4LLhr168ilQmhz3dzdjXQw8wCYKoPw13Au9psaDJbMDVMiPu8VQW/nACQ27tciR1Oob/9OG9mAhG4Q/H8dgHd2t2uicIglBDw53AGWPY0dG0YcAGAMYCQlllyK2dRwkA3L+nE8999l0w6FnhDoAgCKLeNJyAA8Cuzia8OLqw4fXRWUHURSMqLdlLZROCILYYDVdCAYQ6eDCaQnA1ue71i/4ILEYdBtu1K6EQBEFsVRpSwHd3CafhC/7Iutcv+MPY3dUMPZlMEQSxDWhIAb+1vwWMAaemlguvcc5xcTZCHSIEQWwbGlLAHRYjRjodeHtyTcAngzGsJDLY10Mj7gRBbA8aUsABwafkzFSokM7zej7R/eggJbkTBLE9aFgBPzroxEoyg3dmBKfAX14Pwu0wY1jDHnCCIIitTMMK+L27OmDQMfz4whwS6SxeurSAe3a5aciGIIhtQ0P2gQNCCPC7d7nxzMkZdDdbsJrM4PjB3npviyAI4obRsCdwAPjMe3YgGE3hP//gIg57WsklkCCIbYUqAWeMfYAxdpkxdo0x9qhWm5LKbd42/NXHDuLjR/vx3z9+iEKGCYLYViguoTDG9AD+GsD7AMwAeIsx9n3O+UWtNieFhw724iEqnRAEsQ1RcwI/CuAa53yMc54C8C0AD2mzLYIgCKIWagS8F8B00ccz+dfWwRj7NGPsJGPs5OLioorLEQRBEMWoEfByBWe+4QXOn+ScH+GcH3G73SouRxAEQRSjRsBnAPQXfdwHwK9uOwRBEIRU1Aj4WwB2MsYGGWMmAB8D8H1ttkUQBEHUQnEXCuc8wxj71wCeB6AH8DXO+QXNdkYQBEFURdUkJuf8hwB+qNFeCIIgCBk09CQmQRDEdoZxvqFxZPMuxtgigEkVS7QDCGi0Ha3YinsCaF9yoX3Jg/YlD7X78nLON7Tx3VABVwtj7CTn/Ei991HMVtwTQPuSC+1LHrQveWzWvqiEQhAE0aCQgBMEQTQojSbgT9Z7A2XYPVb5BQAABE9JREFUinsCaF9yoX3Jg/Ylj03ZV0PVwAmCIIg1Gu0EThAEQeQhAScIgmhQtryAM8a+xhhbYIydr/deimGM9TPGXmaMjTLGLjDGPlfvPQEAY8zCGHuTMXY2v68/q/eeimGM6Rljpxlj/1TvvYgwxiYYY+cYY2cYYyfrvR8RxlgrY+wZxtil/L+zu7bAnkbyf07ijwhj7PP13hcAMMb+bf7f/HnG2NOMMUu99wQAjLHP5fd0Qes/qy1fA2eMvRvAKoC/55zvr/d+RBhj3QC6OeenGGMOAG8DOH6jE4nK7IsBsHPOVxljRgC/APA5zvnr9dyXCGPsjwEcAdDMOf9wvfcDCAIO4AjnfEsNgDDGvg7gnznnX80bxtk456F670skn8rlA3AH51zNgJ4We+mF8G99L+c8zhj7RwA/5Jz/XZ33tR9C2M1RACkAPwbwh5zzq1qsv+VP4JzzVwAs1XsfpXDOZznnp/K/XgEwijKBFjcaLrCa/9CY/7El3qUZY30APgTgq/Xey1aHMdYM4N0AngIAznlqK4l3nvsBXK+3eBdhAGBljBkA2LA17K33AHidcx7jnGcA/BzAR7VafMsLeCPAGBsAcAjAG/XdiUC+THEGwAKAFzjnW2JfAP4bgC8AyNV7IyVwAD9hjL3NGPt0vTeTZwjAIoC/zZecvsoYs9d7UyV8DMDT9d4EAHDOfQD+AsAUgFkAYc75T+q7KwDAeQDvZoy5GGM2AA9ifY6CKkjAVcIYawLwHQCf55xH6r0fAOCcZznnByGEbBzN38bVFcbYhwEscM7frvdeynCMc34YwAcBfCZftqs3BgCHAXyFc34IQBTAo/Xd0hr5ks5HAPyveu8FABhjbRAyeQcB9ACwM8Z+u767AjjnowD+bwAvQCifnAWQ0Wp9EnAV5GvM3wHwTc75s/XeTyn5W+6fAfhAnbcCAMcAfCRfb/4WgPsYY9+o75YEOOf+/M8LAL4LoV5Zb2YAzBTdPT0DQdC3Ch8EcIpzPl/vjeR5L4Bxzvki5zwN4FkAd9d5TwAAzvlTnPPDnPN3QygHa1L/BkjAFZN/WPgUgFHO+ZfrvR8RxpibMdaa/7UVwj/sS/XdFcA5f4xz3sc5H4Bw6/0S57zuJyTGmD3/EBr5EsX7Idz21hXO+RyAacbYSP6l+wHU9QF5CR/HFimf5JkCcCdjzJb/v3k/hOdSdYcx1pH/2QPgYWj456Yq0OFGwBh7GsC9ANoZYzMA/pRz/lR9dwVAOFF+AsC5fL0ZAP5DPuSinnQD+Hq+Q0AH4B8551umZW8L0gngu8L/eRgA/H+c8x/Xd0sFPgvgm/lyxRiA363zfgAA+Vru+wD8q3rvRYRz/gZj7BkApyCUKE5j64zVf4cx5gKQBvAZzvmyVgtv+TZCgiAIojxUQiEIgmhQSMAJgiAaFBJwgiCIBoUEnCAIokEhAScIgmhQSMAJgiAaFBJwgiCIBuX/B8EWYPqXByfIAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "beta3 = fit_and_plot([2, 9.2, 1, 1])\n", "beta3" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.6505726415425997e-28" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fiterror(beta3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is clear that there are multiple solutions to this problem which are equally good! This is a property of nonlinear regression. It is often impossible to recover the \"right\" values of the parameters. You should therefore be careful of interpreting a good fit as evidence of correctness of your model." ] } ], "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.3" } }, "nbformat": 4, "nbformat_minor": 4 }