{"nbformat":4,"nbformat_minor":0,"metadata":{"kernelspec":{"display_name":"Python 3 (Anaconda)","language":"python","name":"anaconda3"},"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.5.3"},"colab":{"name":"CurveFitting-Working.ipynb","provenance":[],"collapsed_sections":[]}},"cells":[{"cell_type":"markdown","metadata":{"id":"OmsFis6Mtcdt"},"source":["# Curve Fitting"]},{"cell_type":"markdown","metadata":{"id":"nHkH3jjKtcd0"},"source":["Suppose that you want to fit a set of data points $(x_i,y_i)$, where\n","$i = 1,2,\\ldots,N$, to a function that can't be linearized. For example, the function could be a second-order polynomial, $f(x,a,b,c)=ax^2+bx+c$. There isn’t an analytical expression for finding the best-fit parameters ($a$, $b$, and $c$ in this example) as there is for linear regression with uncertainties in one dimension. The usual approach is to optimize the parameters to minimize the sum of the squares of the differences between the data and the function. How the difference are defined varies. If there are only uncertainties in the y direction, then the differences in the vertical direction (the gray lines in the figure below) are used. If there are uncertainties in both the $x$ and $y$ directions, the orthogonal (perpendicular) distances from the line (the dotted red lines in the figure below) are used."]},{"cell_type":"markdown","source":["
"],"metadata":{"id":"MSNQLaybLHkA"}},{"cell_type":"markdown","metadata":{"id":"5J7N0WY_tcd4"},"source":["For the case where there are only uncertainties in the y direction, if the uncertainty in $y_i$ is $\\sigma_i$, then the difference squared for each point is weighted by $w_i=1/\\sigma_i^2$. If there are no uncertainties, each point is given an equal weight of one and the results should be used with caution. The function to be minimized with \n","respect to variations in the parameters is\n","$$\n","\\chi^2 = \\sum_{i=1}^N w_i \\left[y_i - f\\left(x_i,a,b,c\\right)\\right]^2.\n","$$ \n","For the case where there are uncertainties in both $x$ and $y$, the function to be minimized is more complicated. \n","\n","The **`general_fit`** function that performs this minimization is defined in the file \"fitting.py\", which must be in the same drectory as the Python program. If there are only uncertainties in the $y$ direction, it uses the **`curve_fit`** function from the \"scipy\" library to find the best-fit parameters. If there are uncertainties in both $x$ and $y$, the **`odr`** package from the \"scipy\" library is used."]},{"cell_type":"markdown","metadata":{"id":"5IL4jbhstcd5"},"source":["An example of performing a fit with uncertainties in only the $y$ direction is shown below. The first command imports the **`general_fit`** function. Second, the function (**`fitfunc`**) to be fit is defined. The function could have more than 3 parameters. Inital guesses at the parameters are placed in the list $p0$. The name of the fitting function, arrays containing the data points ($x$ and $y$), the initial guesses at the parameters, and an array of uncertainties ($yerr$) are sent to the **`general_fit`** function. If it succeeds, the **`general_fit`** function returns arrays with the optimal parameters and estimates of their uncertainties (in lists called $popt$ and $punc$), the reduced chi squared ($rchi2$), and the degrees of freedom ($dof$). "]},{"cell_type":"code","source":["from google.colab import drive\n","drive.mount('/GoogleDrive')"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"PnDCvjvnMb_q","executionInfo":{"status":"ok","timestamp":1653063821022,"user_tz":420,"elapsed":18183,"user":{"displayName":"Alan DeWeerd","userId":"09086135147919405400"}},"outputId":"e66baa9f-60f7-4b35-f875-b5f795ecb8e5"},"execution_count":1,"outputs":[{"output_type":"stream","name":"stdout","text":["Mounted at /GoogleDrive\n"]}]},{"cell_type":"code","source":["import sys\n","sys.path.append('/GoogleDrive/My Drive/Colab Notebooks/ComputationalTutorials')"],"metadata":{"id":"tpK7p_qPMi7m","executionInfo":{"status":"ok","timestamp":1653063822745,"user_tz":420,"elapsed":8,"user":{"displayName":"Alan DeWeerd","userId":"09086135147919405400"}}},"execution_count":2,"outputs":[]},{"cell_type":"code","metadata":{"id":"4Fefywnhtcd5","outputId":"3fac4076-da9b-40a5-a9d5-3bf6961c16c4","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1653063827524,"user_tz":420,"elapsed":1260,"user":{"displayName":"Alan DeWeerd","userId":"09086135147919405400"}}},"source":["from fitting import *\n","import pylab as pl\n","\n","# Define the fitting function\n","def fitfunc(x, a, b, c):\n"," return (a*x**2 + b*x + c)\n","\n","x = pl.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])\n","y = pl.array([-1.78, 4.09, 8.85, 17.9, 26.1, 35.2])\n","yerr = pl.array([1.24, 1.46, 1.05, 1.68, 1.18, 1.56])\n","\n","p0 = [1.0, 3.0, -2.0]\n","popt, punc, rchi2, dof = general_fit(fitfunc, x, y, p0, yerr)\n","print('optimal parameters: ', popt)\n","print('uncertainties of parameters: ', punc)"],"execution_count":3,"outputs":[{"output_type":"stream","name":"stdout","text":["results of general_fit:\n"," degrees of freedom = 3\n"," reduced chi squared = 0.3960236768747382\n","optimal parameters: [ 0.59883841 4.47445632 -1.7310894 ]\n","uncertainties of parameters: [0.13760793 0.69552672 0.72598714]\n"]}]},{"cell_type":"markdown","metadata":{"id":"c4LAqmlXtcd6"},"source":["For this example, the optimal parameter are\n","\n","$$ a = 0.59883841\\\\\n"," b = 4.47445632\\\\\n"," c = -1.7310894 $$\n","\n","and their uncertainties are\n","\n","$$ \\sigma_a = 0.13760793 \\\\\n"," \\sigma_b = 0.69552672 \\\\\n"," \\sigma_c = 0.72598714 $$\n","\n","If the initial guesses at the parameters are not reasonable, the optimization may fail. It is often helpful to plot the data first to help make a good guesses at the parameters."]},{"cell_type":"markdown","metadata":{"id":"8ateVfb5tcd6"},"source":["If your performing a fit with uncertainties in both the $x$ and $y$ direction, arrays containing the data points ($x$ and $y$) and their uncertainties ($yerr$ and $xerr$) are sent to the **`general_fit`** function as follows: \n","
\n","Note the order of the uncertainties! The uncertainty in $x$ is optional, so it is second. This is also consistent with the **`errorbar`** function."]},{"cell_type":"markdown","metadata":{"id":"OcKm4mbztcd7"},"source":["## Intrepeting the Results"]},{"cell_type":"markdown","metadata":{"id":"IuM2lBBwtcd8"},"source":["Plotting data with error bars and a best-fit function together gives some idea of whether or not the fit is good. If the curve passes within most of the error bars, the fit is probably reasonably good. The first line below makes a list of 100 points between the minimum and maximum values of $x$ in the data. In the second line below, all of the parameters are sent to the fitting function at once using a pointer (using an asterisk in front of the name)."]},{"cell_type":"code","metadata":{"id":"VAyVWUoCtcd9","outputId":"2baf0af9-3954-40b2-d89e-cc861851fb0e","colab":{"base_uri":"https://localhost:8080/","height":279},"executionInfo":{"status":"ok","timestamp":1653063840930,"user_tz":420,"elapsed":810,"user":{"displayName":"Alan DeWeerd","userId":"09086135147919405400"}}},"source":["xf = pl.linspace(min(x),max(x),100)\n","yf = fitfunc(xf,*popt)\n","\n","pl.figure()\n","pl.scatter(x,y,s=15,label='Data')\n","pl.errorbar(x, y, yerr, ls='None', capsize=2)\n","pl.plot(xf,yf,\"r-\",label='Best-Fit Curve')\n","pl.xlabel('x')\n","pl.ylabel('y')\n","pl.legend(loc='upper left')\n","pl.show()"],"execution_count":4,"outputs":[{"output_type":"display_data","data":{"text/plain":[""],"image/png":"iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5zOdf7/8cfLEEIlDTulopYV1QwmHaRdOaTVlk2S6PyLQgdbvumwHVbbdpChNklRKopUmxxikpIcZzRqQhRqHWKizUhhZt6/P95XcpipwXVdn+vwvN9u122uz+f6XHO9rjTv1+fz/rzfr7c55xARkeRTIegAREQkGEoAIiJJSglARCRJKQGIiCQpJQARkSRVMegA9sdRRx3l6tWrF3QYIiJxJTc391vnXOre++MqAdSrV4+cnJygwxARiStm9lVp+9UFJCKSpJQARESSlBKAiEiSiqt7AKXZuXMna9as4aeffgo6FClDlSpVqFu3LpUqVQo6FBHZTdwngDVr1lCjRg3q1auHmQUdjuzFOcemTZtYs2YN9evXDzocEdlN3HcB/fTTT9SqVUuNf4wyM2rVqqUrNJEYFPcJAFDjH+P07yMSmxIiAYiIyP5TAgiDlJQUMjIySE9Pp1mzZsyZM+eAfs+QIUPYtm1bqa+9//77HH744WRkZJCRkUHbtm3Jycnh5ptv3vX6r33u1KlTyczMpHHjxjRt2pTbbrvtgGIUkejKyl5OvQGTdz2yspeH7XfH/U3gWFC1alXy8vIAmDZtGnfeeScffPDBfv+eIUOG0KNHDw499NBSX2/VqhWTJk3aY19mZibgE0D16tU566yz9nlffn4+ffv2ZfLkyTRq1Iji4mJGjBhR7riKioqoWFH/q4gEoV+7hsxbuQmAcb3ODOvv1hVAmG3ZsoWaNWvu2n7sscc47bTTOPXUU7nvvvsA+OGHH+jYsSPp6emcfPLJjBs3jieeeIJ169bRunVrWrduXa7Pev/997ngggtYvXo1w4cPJysri4yMDD788MM9jnv00Ue5++67adSoEeCvWG688UYArr76aiZMmLDr2OrVq+/63a1ateLCCy+kcePGDBgwgKeeemrXcffffz+DBg0q8zuKSOxLrNO6W2+F0Jl42GRkwJAhv3rIjz/+SEZGBj/99BPr16/nvffeA2D69OmsWLGCBQsW4JzjwgsvZNasWRQUFHD00UczefJkAL7//nsOP/xwBg8ezMyZMznqqKNK/ZwPP/yQjIwMALp06ULLli0BXyPphhtuoHr16tx+++37vC8/P/+AunwWLVpEfn4+9evX5+OPP+bWW2+lT58+AIwfP55p06aV+R3POeec/f48EYmuxEoAAdm9C2ju3LlceeWV5OfnM336dKZPn07Tpk0B2Lp1KytWrKBVq1bcdttt3HHHHVxwwQW0atWqXJ+zdxfQ+++/H/bvsrsWLVrsGrvftGlTNm7cyLp16ygoKKBmzZoce+yxDB06tNTvqAQgEj4ViosoSQl/c51YCeA3ztSj4cwzz+Tbb7+loKAA5xx33nknvXr12ue4RYsWMWXKFO655x7atGnDvffeu8frb775Jg888AAAzz333EHF1KRJE3Jzc0lPT9/ntYoVK1JSUgJASUkJO3bs2PVatWrV9ji2S5cuTJgwgW+++YauXbsC/Op3FJEwyMnh8Qe6M+T6gUCc3QMwsypmtsDMFpvZZ2b2QGj/C2a2yszyQo+MSMcSDcuWLaO4uJhatWpx3nnnMWrUKLZu3QrA2rVrd51FH3roofTo0YP+/fuzaNEiAGrUqEFhYSEAf/3rX8nLyyMvL2/Xjd5fs/t799a/f38eeughli/3owdKSkoYPnw44LuPcnNzAZg4cSI7d+4s8zO6du3Kq6++yoQJE+jSpQtAmd9RRMLgjTdw55xDyY4dLP7mB7oMn0NB4faw/fpoXAFsB851zm01s0rAbDObGnqtv3Nuwq+8Ny78fA8A/Bnx6NGjSUlJoX379ixdupQzz/RZu3r16rz88st88cUX9O/fnwoVKlCpUiWefvppAHr27EmHDh04+uijmTlz5n7F8Je//IVLLrmEt956iyeffHKPbqVTTz2VIUOG0K1bN7Zt24aZccEFFwBw/fXXc9FFF5Genk6HDh32OevfXZMmTSgsLOSYY44hLS0NoMzvWLt27f2KX0R24xw89hjccQcr6jehW8c72VTtCFK+/o7eY3J57YZ9R/sdCHPOheUXlevDzA4FZgM3hh6T9icBZGZmur0XhFm6dCknnXRSWOOU8NO/k0g57dgBN94Io0ZB1640O6Ebm0t+OVevVjmFzx7osF+/0sxynXP7dCVEZRiomaWYWR6wEch2zs0PvfRPM/vEzLLMrHIZ7+1pZjlmllNQUBCNcEVEgrF5M5x3nm/877kHxo7lxON+WckxpQI0TjssbB8XlQTgnCt2zmUAdYEWZnYycCfQCDgNOBK4o4z3jnDOZTrnMlNT91nSUkQkMaxYAWecAXPmwIsvwsCBUKECw7o3p0blilQwaHZcTYZ1bx62j4zqKCDn3P/MbCbQwTk3KLR7u5k9D+w7gF1EJBnMnAmdO0OFCjBjBpx99q6XUmtUpvHR/qw/7mYCm1mqmR0Rel4VaAcsM7O00D4DOgH5kY5FRCTmPPcctG8Pv/sdLFiwR+MfadG4AkgDRptZCj7hjHfOTTKz98wsFTAgD7ghCrGIiMSG4mIYMAAGDfIJYPx4OPzwqIYQ8QTgnPsEaFrK/nMj/dkiIjFpyxbo3h0mTYK+fSErCwIouKhicGHwcznoJk2akJ6ezuOPP75rdm1ZVq9ezdixY6MUoYjEjNWroWVLmDoVhg2DJ58MpPEHJYCw+LkW0GeffUZ2djZTp07dVcahLEoAIklo9mxo0QLWrIF33vHj/X9DVvZy5q/azPxVm8O+HgDOubh5NG/e3O1tyZIl++yLtmrVqu2x/eWXX7ojjzzSlZSUuFWrVrmzzz7bNW3a1DVt2tR99NFHzjnnTj/9dHfYYYe59PR0N3jw4DKPSxSx8O8kEqjnn3euUiXnGjRwbtmyqH40kONKaVMTqxhcORQUbqf3mFyWrN9C47TDGNa9Oak1Sp2DdsBOOOEEiouL2bhxI7Vr1yY7O5sqVaqwYsUKunXrRk5ODg8//DCDBg3aVd1z27ZtpR4nInGuuBjuuAMefxzatvU3e3dbMyRISZcAeo/JZdHX31FcAovCXFejNDt37qRv377k5eWRkpKyqyDbgR4nInHk++/h8sthyhTo3dtXLK5UKeiodkm6BLBk/RaKQ/dni0v8dritXLmSlJQUateuzQMPPECdOnVYvHgxJSUlVKlSpdT3ZGVlles4EYkTX34Jf/mLn+H79NNwQ+yNdE+6m8CN0w4jJfStw11XA6CgoIAbbriBvn37YmZ8//33pKWlUaFCBV566SWKi4uBfcs3l3WciMShGTPgtNNgwwaYPj0mG39IwgQwrHtzmh1Xk2qVU8JWV+PnctBNmjShbdu2tG/fftfauL1792b06NGkp6ezbNmyXeWWTz31VFJSUkhPTycrK6vM40QkjjgH//63L+iWluZn9pZzje8gRLUc9MFSOej4pX8nSXg7dvhJXc8+67t+Xn4ZDgtvD8OBCrQctIhIQtu4Edq08Y3/nXfCf/4TM43/r0m6m8AiImGVlwcXXeSTwNix0K1b0BGVW0JcAcRTN1Yy0r+PJKxx4+Css6CkxM/yjaPGHxIgAVSpUoVNmzapkYlRzjk2bdqkYa2SWIqL4a674LLLoFkzWLgQmodvoZZoifsuoLp167JmzRq0XGTsqlKlCnXr1g06DJHw+P57X8lz8mS4/no/6ueQQ4KO6oDEfQKoVKkS9evXDzoMEUkGy5ZBp05+ktdTT/libmZBR3XA4j4BiIhExaRJ/sy/cmU/0eucc4KO6KBFY0nIKma2wMwWm9lnZvZAaH99M5tvZl+Y2Tgzi89rKBFJbCUl8OCDcOGF8PvfQ05OQjT+EJ2bwNuBc51z6UAG0MHMzgAeAbKcc78HvgOui0IsIiLlV1gIl1wCf/+7P/ufPRuOOy7oqMIm4gkgVI56a2izUujhgHOBCaH9o/ELw4uIxIYVK+CMM2DiRBg8GF58EapWDTqqsIrKMFAzSzGzPGAjkA18CfzPOVcUOmQNcEwZ7+1pZjlmlqORPiISFZMn/1LMbdo06Ncvrm/2liUqCcA5V+ycywDqAi2ARvvx3hHOuUznXGZqamrEYhQRoaQEBg70tXzq1/f9/W3aBB1VxER1FJBz7n9mNhM4EzjCzCqGrgLqAmujGYuIyB6+/x6uugreesv3948YAYceGnRUERWNUUCpZnZE6HlVoB2wFJgJXBI67CrgrUjHIiJSqqVL/WLtkybB0KHw0ksJ3/hDdK4A0oDRZpaCTzjjnXOTzGwJ8KqZPQh8DIyMQiwiInt6/XW4+mrf4M+YAX/8Y9ARRU3EE4Bz7hOgaSn7V+LvB4iIRF9REdxzDzzyiD/7f/11SLKSJZoJLCLJ59tvfeXOd9+FXr18t0/lykFHFXVKACKSXBYuhM6dff3+kSPh2muDjigwcV8OWkSkXJzzK3adfTZUqAAffZTUjT8oAYhIMvjxR7juOujZE/70J8jNjcv6/eGmBCAiiW3lSmjZEp5/3t/0nTIFatUKOqqYoHsAIpK4pkyBHj1898+kSdCxY9ARxRRdAYhI4iku9hU8O3aE44/3XT5q/PehKwARiRtZ2csZOmPFru1b2jSgX7uGex5UUACXX+6HeF57rV+yMcGqeIaLEoCIxI1+7Royb+UmAMb1OnPfA+bOhS5d/Dj/557zN36lTOoCEpH455yfzHXOOX6B9jlz1PiXgxKAiMS3LVuga1e49Vb48599f3+zZkFHFRfUBSQicaOgcDtL1m3hhx1FdBk+h2fSD+HIqy73Qz0feQT690/IhVsiRQlAROJG7zG5FG73CwnWnzSe6jc/DUfVgvfeS5iF2qNJCUBE4saS9VuosvMnBk4fTpf8d5lfL53T502DOnWCDi0uKQGISNxoa99x44t30vDbr3my5WXMvrw349T4HzAlABGJD2PHMmRQT75zFbnm0vvZ1rotw7qrns/BiMaSkMea2UwzW2Jmn5nZLaH995vZWjPLCz3+HOlYRCQO/fijr9nfvTvWtCl3/300P53bjtduOIvUGslXwz+cojEMtAi4zTnXGDgD6GNmjUOvZTnnMkKPKVGIRUTiyfLlcOaZfoH2AQNg5kw216wddFQJIxpLQq4H1oeeF5rZUuCYSH+uiMS5sWP9mX/lyjB5sh/jL2EV1YlgZlYPvz7w/NCuvmb2iZmNMrOaZbynp5nlmFlOQUFBlCIVkcBs2+br9nfvDunpkJe3q/HPyl7O/FWbmb9qM/UGTCYre3nAwcY3c85F54PMqgMfAP90zr1hZnWAbwEHDATSnHO/ujxPZmamy8nJiXywIhKMJUv8rN78fN/lM3AgVNRYlYNlZrnOucy990flv6yZVQJeB8Y4594AcM5t2O31Z4FJ0YhFRGKQczB6NPTpA9WqwTvvwHnnBR1VwovGKCADRgJLnXODd9ufttthfwXyIx2LiMSgwkK48kq45hpo0cJ3+ajxj4poXAG0BK4APjWzvNC+u4BuZpaB7wJaDfSKQiwiEkvy8uDSS+HLL+GBB+DuuyElJeiokkY0RgHNBkqrzqRhnyLJyjm/UMvtt8NRR/laPn/8Y9BRJR3dXRGR6Nq0ya/UNXGiX6bxhRd8EpCo03oAIhI9s2ZBRgZMnQpZWfD222r8A6QEICKRV1QE990HrVtDlSp+6cZbb1Xt/oCpC0hEIuvrr/2krtmz/Wiff/8batQIOipBCUBEImnCBLj+en8F8NJL0KNH0BHJbtQFJCLh98MPvpxDly7QsKEf7qnGP+YoAYhIeH38MWRmwnPP+XIOs2fDiScGHZWUQglARMKjpMSP7DnjDNiyBd59F/71L6hUKejIpAy6ByAiB2/9erj6apg+HTp18mf/tWoFHZX8Bl0BiMjBmTgRTj0VPvwQnn4a3nhDjX+cUAIQkQOzbRvceCNcdBHUrQu5uXDDDRrbH0eUAERk/+XmQrNmMHy4r+czbx6cdFLQUcl+UgIQkfIrLoaHH/Y3erduhRkz4LHH/LKNEnd0E1hEyuerr/xM3lmz4JJL4Jln4Mgjg45KDoKuAETk1zkHY8b49Xk//tiv3DV+vBr/BBCNFcGONbOZZrbEzD4zs1tC+480s2wzWxH6Weqi8CISoM2boVs3P4v35JNh8WJ/FaAbvQkhGlcARcBtzrnGwBlAHzNrDAwAZjjnGgAzQtsiEiuys/3wztdfh3/+Ez74AOrXDzoqCaOIJwDn3Hrn3KLQ80JgKXAMcBEwOnTYaKBTpGMRkXLYtg1uvhnat/dVO+fNg7vu0lKNCSiq9wDMrB7QFJgP1HHOrQ+99A1QJ5qxiEgpcnOheXN48kmfBBYt8tuSkKKWAMysOvA6cKtzbsvurznnHH5x+NLe19PMcswsp6CgIAqRiiShnTvhH//4ZXhndjYMHQpVqwYdmURQVBKAmVXCN/5jnHNvhHZvMLO00OtpwMbS3uucG+Gcy3TOZaampkYjXJHksmwZtGzpV+y67DL45BNo2zboqCQKIj4PwMwMGAksdc4N3u2licBVwMOhn29FOhaRRJKVvZyhM1bs2r6lTQP6tWtY/l9QUuK7egYMgEMPhdde8+P7JWlEYyJYS+AK4FMzywvtuwvf8I83s+uAr4BLoxCLSMLo164h81ZuAmBcrzP3781ffQXXXAMzZ0LHjvDss5CWFoEoJZZFPAE452YDZQ0abhPpzxeR3TgHo0ZBv37++ciRPhFoXH9SUikIkWSxbp1fn3fKFGjd2ieCevWCjkoCpFIQIonu51IOJ5/su3yeeMKv1qXGP+kpAYgksg0boHNnX8qhUSO/OPtNN0EF/emLEoBI3Coo3M6SdVtYuHozXYbPoaBw+54HjB/vz/qnTIFHH/UrdjXcj1FCkvCUAETiVO8xuRRuL6LEwaKvv6P3mFz/wsaN0KULdO3qa/csWgT9+6uUg+xDCUAkTi1Z/8uE+uKS0PZrr0GTJn6d3n/9C+bMgcaNA4xSYpkSgEicapx22K7ntX/8jpFvPwKXXupv7ubm+gleFTXQT8qmBCASp4Z1b06NQ1K4aOkHvDuyL6d/Nsef9c+d6/v+RX6DTg9E4lTq1s2MmvgQpy3+EE4/3Y/rV3eP7IffvAIws5u0WpdIDHHOL8vYpAnpS+bz0sV94aOP1PjLfivPFUAdYKGZLQJGAdNC5ZtFJNq+/hp69YJ33oGWLfm/9jexvs5xXKERPnIAfvMKwDl3D9AAX9HzamCFmT1kZidGODYR+VlJCQwb5kf4fPihr+I5axbr6xwXdGQSx8p1Ezh0xv9N6FEE1AQmmNmjEYxNRACWL4c//Qn69PELtuTnQ9++ms0rB6089wBuMbNc4FHgI+AU59yNQHOgc4TjE0leO3fCI49Aejp8+ik8/zxMn76rhk9W9nLmr9rM/FWbqTdgMlnZy4ONV+JOee4BHAlc7Jz7avedzrkSM7sgMmGJJLlFi+C663ztnosvhn//e596/f3aNdy/BWBE9lKeewD37d347/ba0vCHJJLEtm2DO+6AFi3gm2/g9df9Q4u1SAREvBPRzEaZ2UYzy99t3/1mttbM8kKPP0c6DpGYN2MGnHKKL9x2zTWwdKk/+xeJkGjcRXoB6FDK/iznXEboMSUKcYjEpk2bfIPftq0v2DZzpl+i8Ygjgo5MElzEE4BzbhawOdKfIxJ3nIOxY+Gkk+Dll+HOO2HxYj/iRyQKghxH1tfMPgl1EZU509jMeppZjpnlFBQURDM+kchZtQrOPx+6d/clm3Nz4aGHoGrVoCOTJBJUAngaOBHIANYDj5d1oHNuhHMu0zmXmZqaGq34RCLj56GdTZr48g1Dh/qSzaeeGnRkkoQCKQbnnNvw83MzexaYFEQcIlE1b54v4/DJJ9Cpk5/NW7du0FFJEgvkCsDMdh/T9lcgv6xjReLe99/7WbxnneVv+L75pn+o8ZeARfwKwMxeAf4EHGVma4D7gD+ZWQbggNVAr0jHIRJ1zvl1eW+91S/TePPNMHAg1KgRdGQiQBQSgHOuWym7R0b6c0UC9eWX0Lu3L93QvDm8/TZkZgYdlcgeVE1KJJy2b/dn+U2a+JW5nngC5s9X4y8xSSuCiYTLe+/BjTf66p2XXgqDB8MxxwQdlUiZdAUgcrDWr4fLL4c2baC42C/WMm6cGn+JeUoAIgeqqMh38TRqBG+8Affd58s2n3de0JGJlIu6gEQOxNy5vrtn8WJo396Xa27QIOioRPaLrgBE9kdBga/T//OY/gkTfJePGn+JQ0oAIuVRXAxPPw0NG8KLL0L//r5cc+fOYBZ0dCIHRF1AIr9l7ly/Bu+iRXDuub6756STgo5K5KDpCkCkLBs2+Dr9Z53ln7/6Krz7rhp/SRhKACJ727kThgzx3T1jxsCAAbBsGXTtqu4eSSjqAhLZ3XvvwU03wZIlfnTPE0/AH/4QdFQiEaErABGA1avhkkv8ZK4ff4T//MeP7lHjLwlMVwCS3LZt84uwP/IIVKjg6/jcfjtUqRJ0ZCIRpwQgyck5eO0139j/979w2WU+ERx7bNCRiUSNuoAk+eTl+YXXu3aFWrXggw/glVfU+EvSUQKQ5LFhA1x/PTRr5m/yPvMM5OTAOecEHZlIICKeAMxslJltNLP83fYdaWbZZrYi9LNmpOOQJLZ9Ozz2mC/X8MIL0K8frFgBPXtCSkrQ0YkEJhpXAC8AHfbaNwCY4ZxrAMwIbYuEl3N+7d3GjeH//s+f6efnw+OPwxFHBB2dSOAingCcc7OAzXvtvggYHXo+GugU6TgkyXz8sS/bcPHFULUqTJsGkyZpWKfIboK6B1DHObc+9PwboE5ZB5pZTzPLMbOcgoKC6EQn8Wv9erj2Wr8Ob34+PPWUv+nbvn3QkYnEnMBvAjvnHOB+5fURzrlM51xmampqFCOTePLkxDwGt+rBtuPrs2P0S+R0vsb38/fuDRU12lmkNEH9ZWwwszTn3HozSwM2BhSHxLviYnjpJW66+25Yt455zVpzxvhnyTzxxKAjE4l5QV0BTASuCj2/CngroDgknr37ru/queYaqFuXe28fTlbPf4Iaf5FyicYw0FeAucAfzGyNmV0HPAy0M7MVQNvQtkj5fPYZdOwI7drB//4Hr7xCwfT3efPQeixcvZkuw+dQULg96ChFYl7Eu4Ccc93KeKlNpD9bEsy6dXDvvfD881Cjhi/dcNNNUKUKvYfPoXB7EQCLvv6O3mNyee2GswIOWCS26e6YxL4tW/xErscfh6IiuPlmuOceX8YhZMn6LbueF5fsuS0ipQt8FJBImXbs8Msv/v738OCDcOGFfh3erKw9Gn+AxmmH7XqeUmHPbREpnRKAxJ6SEhg3zs/gvekmaNIEFizwSzKWcYN3WPfm1KhckQoGzY6rybDuzaMctEj8UReQxJYZM/wSjDk5cMopMGUKdOjwm0sxptaoTOOj/Vn/uF5nRiNSkbinKwCJDYsW+dm6bdvCxo0werQv53D++VqHVyRClAAkWMuX+7r8zZv7JDB4MHz+OVx55X5V6szKXs78VZuZv2oz9QZMJit7eQSDFkkM5isxxIfMzEyXk5MTdBgSDmvXwj/+ASNH+uUX//Y3uO02OPzwoCMTSThmluucy9x7v+4BSHR9+y08/LAv0lZc7Gv13H031CmzHqCIRIgSgETHli2+e2fwYPjhB7jiCrj/fqhXL+jIRJKWEoBE1rZt/mz/4Ydh82bo3Nl3/TRuHHRkIklPN4ElMrZv95O4TjzRr8bVooUf2jlhghp/kRihKwAJr507/bq7Dz4IX3/tl2EcPx5atQo6MhHZi64AJDyKivzY/UaN/GLrv/sdZGfD+++r8ReJUUoAcnCKi2HMGN+tc/XVfrH1SZNg3jw/qUuTuERilhKAHJjiYnjlFTj5ZOjRwy+8/uabvp+/Y0c1/CJxQAlA9k9xsS/KdsopcPnlfrbuhAm+bEOnTmr4ReJIoDeBzWw1UAgUA0WlzVST/ZeVvZyhM1bs2r6lTQP6tWt4cL+0uNhX6HzwQV+SuUkTf3O3c2eooPMIkXgUC6OAWjvnvg06iETSr11D5q3cBIShMmZRkT/jf/BBX6Pn5JPV8IskCP0FS+l+Hs550kl+1m7lyr6rZ/Fi6NJFjb9IAgj6r9gB080s18x6lnaAmfU0sxwzyykoKIhyeElo+3YYMQL+8Ae45ho47DB/c/fjj3XWL5Jggv5rPts51ww4H+hjZufsfYBzboRzLtM5l5mamhr9CJPFtm3wxBN++cVevSA1Fd5+24/q6dRJDb9IAgr0r9o5tzb0cyPwJtAiyHiS0pYt8MgjUL8+3HKL/zl9uh/Hf8EFGtUjksACSwBmVs3Mavz8HGgP5AcVTyIpKNzOknVbWLh6M12Gz6GgcPu+B337Ldx7Lxx/vF+CsWlTmDXLP9q1U8MvkgSCvAKoA8w2s8XAAmCyc+6dAONJGL3H5FK4vYgSB4u+/o7eY3J/eXHNGr/4yvHHw8CBcO65sHAhvPOOSjaIJJnAhoE651YC6UF9fiJbsn7LrufFJaHtzz+HRx+Fl16CkhI/iWvAAFXmFElisTAPQMKscdphLFz9HQBNv/mcOxZPhIGz/FDOnj3h9tu1EIuIKAEkomGXN+OeGx/nqtnjOevrTyg54gi46y64+WaoXTvo8EQkRigBJJKdO2H8eFIffZRnPvmEzYcfBYMGUaFnT6hRI+joRCTGKAEkgsJCePZZGDIE/vtfaNyYp666h49Oa8fYPvtMrRARAYKfCCYHY+1auOMOOPZYuO02P4b/7bfh00+ZdeafKa5YKegIRSSG6QogHi1eDIMHw9ixfkTPxRdD//5+3V0RkXLSFUC8KCmBKVP8KlsZGb4wW+/e8MUX8NprezT+WdnLmb9qM/NXbabegMlkZS8PMHARiVXmnAs6hnLLzMx0OTk5QYcRXdu2+bH7Q4bAsmVw9NF+NE/PnlCzZreuB5EAAAfYSURBVNDRiUgcMLPc0tZbURdQrFq7FoYNg2eegU2boHlznwguvRQOOSTo6EQkASgBxJoFC/zZ/muv+VW4LrrIl244+2zV5xGRsFICiAU7d/o+/aFDYf58P2b/ppugb1844YSgoxORBKUEEKQNG/ziK08/DevXQ4MGvib/1Vdr4paIRJwSQBAWLIAnn/Rr6+7YAR06wHPP+Z9aeEVEokQJIFp++gnGjYOnnvLll2vU8Ctv9enjl18UEYkyJYBIW7XKd/GMGuVH85x0kk8CV1yhbh4RCZQSQCQUF8PUqb7hnzrVd+t06uQnbrVurdE8IhITAk0AZtYBGAqkAM855x4OMp6D9s03/kx/xAj46itIS4O//x2uvx7q1g06OhGRPQSWAMwsBXgKaAesARaa2UTn3JKgYjogJSUwc6afsPXmm1BUBG3awKBBfgx/JRVkE5HYFOQVQAvgi9DSkJjZq8BFQHwkgI0bYfRof7b/xRdw5JG+REOvXtCwYdDRiYj8piATwDHAf3fbXgOcvvdBZtYT6Alw3HHHRSeysvx8tj9ihD/b37nTz9C97z645BKoUiXY+ERE9kPM3wR2zo0ARoAvBhdIEOvWwQsvwMiRsHKlP9vv08f37WtRdRGJU0EmgLXAsbtt1w3tiw07d8Lkyb7RnzrVj+xp3RoGDvT193W2LyJxLsgEsBBoYGb18Q3/ZcDl4f6QrOzlDJ2xYtf2LW0a0K/dr/TRL1sGzz/v+/c3bPAjefr3h2uv9aUaREQSRGAJwDlXZGZ9gWn4YaCjnHOfhftz+rVryLyVmwAY1+vM0g/assWXZRg1CubOhYoVoWNHuO46OP98vy0ikmACbdmcc1OAKYF8eHGxv6H7wgvwxhvw449+lu5jj/lZunXqBBKWiEi0JN+p7eefw4sv+seaNXD44XDVVXDNNXDaaZqlKyJJI+ETQEHhdtZ8uZY2i99j+cP/j4arl/jSDOed5ydrXXghVK0adJgiIlGX8Amg95hc/jb1aTrnv8ey1Hq8eHEfrnzyLr+2rohIEkv4BLBk/RaGnd6F507rxNLaJ1CtcgpXqvEXEUn8BNA47TAWbvfTDVIq+G0REYGEX35qWPfm1KhckQoGzY6rybDuzYMOSUQkJiT8FUBqjco0Ptqf9Zc5D0BEJAkl/BWAiIiUTglARCRJJXwCyMpezvxVm5m/ajP1BkwmK3t50CGJiMQEcy6YCssHIjMz0+Xk5AQdhohIXDGzXOdc5t77E/4KQERESqcEICKSpJQARESSlBKAiEiSUgIQEUlSgSQAM7vfzNaaWV7o8ecg4hARSWZBloLIcs4NCvDzRUSSmrqARESSVJBXAH3N7EogB7jNOfddaQeZWU+gZ2hzq5l9foCfdxTw7QG+N17pOycHfefkcDDf+fjSdkZsJrCZvQv8rpSX7gbm4b+IAwYCac65ayMSyC/x5JQ2Ey6R6TsnB33n5BCJ7xyxKwDnXNvyHGdmzwKTIhWHiIiULqhRQGm7bf4VyA8iDhGRZBbUPYBHzSwD3wW0GugVhc8cEYXPiDX6zslB3zk5hP07x1U1UBERCR8NAxURSVJKACIiSSopEoCZdTCzz83sCzMbEHQ8kWZmo8xso5klxc11MzvWzGaa2RIz+8zMbgk6pkgzsypmtsDMFoe+8wNBxxQtZpZiZh+bWVKMHjSz1Wb2aahsTlhXxEr4ewBmlgIsB9oBa4CFQDfn3JJAA4sgMzsH2Aq86Jw7Oeh4Ii00qizNObfIzGoAuUCnBP83NqCac26rmVUCZgO3OOfmBRxaxJnZ34BM4DDn3AVBxxNpZrYayHTOhX3iWzJcAbQAvnDOrXTO7QBeBS4KOKaIcs7NAjYHHUe0OOfWO+cWhZ4XAkuBY4KNKrKctzW0WSn0SOyzOcDM6gIdgeeCjiURJEMCOAb4727ba0jwxiGZmVk9oCkwP9hIIi/UFZIHbASynXMJ/52BIcD/ASVBBxJFDphuZrmh0jhhkwwJQJKEmVUHXgdudc5tCTqeSHPOFTvnMoC6QAszS+juPjO7ANjonMsNOpYoO9s51ww4H+gT6uINi2RIAGuBY3fbrhvaJwkk1A/+OjDGOfdG0PFEk3Puf8BMoEPQsURYS+DCUJ/4q8C5ZvZysCFFnnNubejnRuBNfLd2WCRDAlgINDCz+mZ2CHAZMDHgmCSMQjdERwJLnXODg44nGsws1cyOCD2vih/ksCzYqCLLOXenc66uc64e/u/4Pedcj4DDiigzqxYa2ICZVQPaE8bSOQmfAJxzRUBfYBr+5uB459xnwUYVWWb2CjAX+IOZrTGz64KOKcJaAlfgzwiTZZW5NGCmmX2CP8nJds4lxbDIJFMHmG1mi4EFwGTn3Dvh+uUJPwxURERKl/BXACIiUjolABGRJKUEICKSpJQARESSlBKAiEiSUgIQEUlSSgAiIklKCUDkIJjZaWb2Sag+f7VQbf6ErskjiUMTwUQOkpk9CFQBqgJrnHP/CjgkkXJRAhA5SKEaUwuBn4CznHPFAYckUi7qAhI5eLWA6kAN/JWASFzQFYDIQTKzifjyxPXxS1P2DTgkkXKpGHQAIvHMzK4EdjrnxobWn55jZuc6594LOjaR36IrABGRJKV7ACIiSUoJQEQkSSkBiIgkKSUAEZEkpQQgIpKklABERJKUEoCISJL6/3OjXA5CL7LxAAAAAElFTkSuQmCC\n"},"metadata":{"needs_background":"light"}}]},{"cell_type":"markdown","metadata":{"id":"xVdSGkRRtcd-"},"source":["The reduced chi squared and the degrees of freedom can also be used to judge the goodness of the fit. If $N$ is the number of data points and $C$ is\n","the number of parameters (or constraints) in the fit, the number degrees of freedom is \n","$$\n","d = N - C.\n","$$\n","In the example, $C = 3$ because there three parameters in the function. The reduced chi squared is defined as\n","$$\n","\\tilde{\\chi}^{\\, 2} = \\frac{\\chi^2}{d}.\n","$$\n","According to Taylor (p. 271), “If we obtain a value of $\\tilde{\\chi}^{\\, 2}$ of order one or less, then we have no reason to doubt our expected distribution; if we obtain a value of $\\tilde{\\chi}^{\\, 2}$ much larger than one, our expected distribution is unlikely to be correct.” \n","For an observed value (from fitting data) of the reduced chi square ($\\tilde{\\chi}^{\\, 2}_o$), you can look up the probability of randomly getting a larger $\\tilde{\\chi}^{\\, 2}$ with $d$ degrees of freedom on the table below (from Appendix D of Taylor’s book). A typical standard is to reject a fit if \n","$$\n","Prob_d\\left(\\tilde{\\chi}^{\\, 2} \\ge \\tilde{\\chi}^{\\, 2}_o \\right) < 5\\%.\n","$$\n","In other words, if the reduced chi squared for a fit is unlikely to occur randomly, then the fit is not a good one. \n","In the example above, six data points are fit and $\\tilde{\\chi}^{\\, 2}_o = 0.39$. Since $d = 6 - 3 = 3$, the table gives\n","$$\n","Prob_d\\left(\\tilde{\\chi}^{\\, 2} \\ge \\tilde{\\chi}^{\\, 2}_o \\right) \\approx 75\\%,\n","$$\n","and there is no reason to reject the fit."]},{"cell_type":"markdown","source":["
\n","From Error Analysis by John Taylor
"],"metadata":{"id":"7a9lLmhlLuJG"}},{"cell_type":"markdown","metadata":{"id":"U-nUR2Smtcd_"},"source":["## Additional Documentation"]},{"cell_type":"markdown","metadata":{"id":"OXpKndqftcd_"},"source":["More information is available at: \n","http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html \n","https://docs.scipy.org/doc/scipy/reference/odr.html"]}]}