{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hxdc7U-1kKL0"
   },
   "source": [
    "$$\\newcommand{\\arr}[1]{\\underline{\\underline{#1}}}$$    \n",
    "$$\\newcommand{\\vec}[1]{\\underline{#1}}$$   \n",
    "$$\\require{mhchem}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Non-linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "n63PGfye9pig"
   },
   "source": [
    "## Recap on linear regression\n",
    "\n",
    "Last class, we talked about how we could turn linear regression into a linear algebra problem\n",
    "- You can calculate this yourself\n",
    "- You can also do linear regression using statistical packages like statsmodel\n",
    "\n",
    "Today we will discuss two ways of solving non-linear regression problems\n",
    "- Turn a non-linear problem into a linear one and solve\n",
    "- Non-linear curve fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NDQBwO8M2Tmt"
   },
   "source": [
    "\n",
    "## Turn a non-linear regression problem into a linear regression problem\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PXofbsJg2Tmu"
   },
   "source": [
    "-   Rate constants and reaction orders are determined by using models that are fit to experimental data\n",
    "-   A common case is to monitor concentration vs. time in a constant volume, batch reactor\n",
    "-   We consider the disappearance of $A$\n",
    "-   From the mole balance we know:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dN_A}{dt} = r_A V\n",
    "\\end{equation}\n",
    "\n",
    "-   Let us assume the rate law is of the form: $r_A = k C_A^\\alpha$ and a constant volume so that:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dC_A}{dt} = -k C_A^\\alpha\n",
    "\\end{equation}\n",
    "\n",
    "-   Let us be loose with mathematics, rearrange the equation, and take the log of both sides.\n",
    "    -   By loose I mean we take logs of quantities that are not dimensionless\n",
    "\n",
    "\\begin{equation}\n",
    "\\ln(-\\frac{dC_A}{dt}) = \\ln{k} + \\alpha \\ln C_A\n",
    "\\end{equation}\n",
    "\n",
    "-   This suggests that if we could numerically compute $\\frac{dC_A}{dt}$ from our data of $C_A(t)$ then a plot of the log of the negative derivative vs the log of concentration would have\n",
    "    -   an intercept equal to the log of the rate constant, $k$\n",
    "    -   and a slope equal to the reaction order $\\alpha$\n",
    "\n",
    "-   Given the following data, determine the reaction order in A and the rate constant with 95% confidence intervals.\n",
    "\n",
    "<table id=\"orgfd5ae74\" border=\"2\" cellspacing=\"0\" cellpadding=\"6\" rules=\"groups\" frame=\"hsides\">\n",
    "\n",
    "\n",
    "<colgroup>\n",
    "<col  class=\"org-right\" />\n",
    "\n",
    "<col  class=\"org-right\" />\n",
    "</colgroup>\n",
    "<thead>\n",
    "<tr>\n",
    "<th scope=\"col\" class=\"org-right\">time (min)</th>\n",
    "<th scope=\"col\" class=\"org-right\">C\\_A (mol/L)</th>\n",
    "</tr>\n",
    "</thead>\n",
    "\n",
    "<tbody>\n",
    "<tr>\n",
    "<td class=\"org-right\">0</td>\n",
    "<td class=\"org-right\">0.0500</td>\n",
    "</tr>\n",
    "\n",
    "\n",
    "<tr>\n",
    "<td class=\"org-right\">50</td>\n",
    "<td class=\"org-right\">0.0380</td>\n",
    "</tr>\n",
    "\n",
    "\n",
    "<tr>\n",
    "<td class=\"org-right\">100</td>\n",
    "<td class=\"org-right\">0.0306</td>\n",
    "</tr>\n",
    "\n",
    "\n",
    "<tr>\n",
    "<td class=\"org-right\">150</td>\n",
    "<td class=\"org-right\">0.0256</td>\n",
    "</tr>\n",
    "\n",
    "\n",
    "<tr>\n",
    "<td class=\"org-right\">200</td>\n",
    "<td class=\"org-right\">0.0222</td>\n",
    "</tr>\n",
    "\n",
    "\n",
    "<tr>\n",
    "<td class=\"org-right\">250</td>\n",
    "<td class=\"org-right\">0.0195</td>\n",
    "</tr>\n",
    "\n",
    "\n",
    "<tr>\n",
    "<td class=\"org-right\">300</td>\n",
    "<td class=\"org-right\">0.0174</td>\n",
    "</tr>\n",
    "</tbody>\n",
    "</table>\n",
    "\n",
    "\n",
    "- We can get the derivatives by first fitting a spline through the data. The spline is essentially just a smoothing function  \n",
    "-   We will use the `splev` function to numerically compute derivatives from spline fits of the function. \n",
    "-   This works best when the $x$ points are evenly spaced, and they should be monotically increasing or decreasing\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 297
    },
    "executionInfo": {
     "elapsed": 811,
     "status": "ok",
     "timestamp": 1651147976426,
     "user": {
      "displayName": "Zachary Ulissi",
      "userId": "07633171379186475882"
     },
     "user_tz": 240
    },
    "id": "bidKH_Yo2Tmv",
    "outputId": "f5f97a6d-38ce-40d8-df9c-cba0dfe3dafe"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Conc. [mol/L]')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAbgElEQVR4nO3df5BV5Z3n8ffH5ocdo3YGO440kMbAkMWwK3qDyYxJJbEMmJ0RYqgNbnblDzeM41A7birswGYrZdyZdZSZmMrETcqMbIiziWYYIJ0f2iaDMbOpldAdMA06nWnRlDRORBQMkx4F/O4f92m9XG5339P06Xtv9+dVdavvec5zzv0+nIYP58c9RxGBmZlZtc6qdQFmZtZYHBxmZpaJg8PMzDJxcJiZWSYODjMzy2RKrQsYDxdccEG0t7fXugwzs4bS3d39QkS0lrdPiuBob2+nq6ur1mWYmTUUSb+o1O5DVWZmlomDw8zMMnFwmJlZJg4OMzPLxMFhZmaZ5BockpZJ6pXUJ2l9hfnTJT2Q5u+U1J7a2yUNSNqTXl8uWeZyST1pmS9IUh61b9/dz+/82Q7mrv8uv/NnO9i+uz+PjzEzazi5BYekJuBu4BpgIXC9pIVl3W4EXoqIecBdwB0l856KiEvT66aS9i8BnwDmp9eysa59++5+Nmztof/IAAH0Hxlgw9Yeh4eZGfnucSwB+iJif0S8CtwPLC/rsxzYnN5vAa4abg9C0kXAeRHxWBTvB/81YMVYF76xs5eB4ydPaRs4fpKNnb1j/VFmZg0nz+BoA54tmT6Q2ir2iYgTwFFgRpo3V9JuSY9Kem9J/wMjrBMASWskdUnqOnToUKbCDx4ZyNRuZjaZ1OvJ8eeAORGxGPgk8HVJ52VZQUTcExGFiCi0tp72jflhzWxpztRuZjaZ5Bkc/cDskulZqa1iH0lTgPOBwxHxSkQcBoiIbuAp4LdS/1kjrPOMrVu6gOapTae0NU9tYt3SBWP9UWZmDSfP4NgFzJc0V9I0YBXQUdanA1id3q8EdkRESGpNJ9eRdDHFk+D7I+I54GVJ707nQm4AvjXWha9Y3Mbt1y2iraUZAW0tzdx+3SJWLK54VMzMbFLJ7SaHEXFC0lqgE2gCNkXEPkm3AV0R0QHcC9wnqQ94kWK4ALwPuE3SceA14KaIeDHNuxn4KtAMPJheY27F4jYHhZlZBSpenDSxFQqF8N1xzcyykdQdEYXy9no9OW5mZnXKwWFmZpk4OMzMLBMHh5mZZeLgMDOzTBwcZmaWiYPDzMwycXCYmVkmDg4zM8vEwWFmZpk4OMzMLBMHh5mZZeLgMDOzTBwcZmaWiYPDzMwycXCYmVkmDg4zM8vEwWFmZpk4OMzMLBMHh5mZZeLgMDOzTHINDknLJPVK6pO0vsL86ZIeSPN3Smovmz9H0jFJnyppe0ZSj6Q9krryrN/MzE6XW3BIagLuBq4BFgLXS1pY1u1G4KWImAfcBdxRNv9zwIMVVv+BiLg0IgpjXLaZmY0gzz2OJUBfROyPiFeB+4HlZX2WA5vT+y3AVZIEIGkF8DSwL8cazcwsozyDow14tmT6QGqr2CciTgBHgRmS3gz8MfDZCusN4GFJ3ZLWjHnVZmY2rCm1LmAItwJ3RcSxtANS6sqI6Jf0VuD7kv4hIn5U3imFyhqAOXPm5F2vmdmkkeceRz8wu2R6Vmqr2EfSFOB84DBwBXCnpGeAW4D/JmktQET0p5/PA9soHhI7TUTcExGFiCi0traO1ZjMzCa9PINjFzBf0lxJ04BVQEdZnw5gdXq/EtgRRe+NiPaIaAc+D/zPiPiipHMknQsg6RzgQ8DeHMdgZmZlcjtUFREn0l5CJ9AEbIqIfZJuA7oiogO4F7hPUh/wIsVwGc6FwLZ0+GoK8PWIeCivMZiZ2ekUEbWuIXeFQiG6uvyVDzOzLCR1V/rag785bmZmmTg4zMwsEweHmZll4uAwM7NMHBxmZpaJg8PMzDJxcJiZWSYODjMzy8TBYWZmmTg4zMwsEweHmZll4uAwM7NMHBxmZpaJg8PMzDJxcJiZWSYODjMzy8TBYWZmmTg4zMwsEweHmZll4uAwM7NMHBxmZpaJg8PMzDLJNTgkLZPUK6lP0voK86dLeiDN3ympvWz+HEnHJH2q2nWamVm+cgsOSU3A3cA1wELgekkLy7rdCLwUEfOAu4A7yuZ/Dngw4zrNzCxHee5xLAH6ImJ/RLwK3A8sL+uzHNic3m8BrpIkAEkrgKeBfRnXaWZmOcozONqAZ0umD6S2in0i4gRwFJgh6c3AHwOfHcU6AZC0RlKXpK5Dhw6NehBmZnaqej05fitwV0QcG+0KIuKeiChERKG1tXXsKjMzm+Sm5LjufmB2yfSs1FapzwFJU4DzgcPAFcBKSXcCLcBrkv4F6K5inWZmlqM8g2MXMF/SXIr/uK8C/n1Znw5gNfD/gJXAjogI4L2DHSTdChyLiC+mcBlpnWZmlqPcgiMiTkhaC3QCTcCmiNgn6TagKyI6gHuB+yT1AS9SDILM68xrDGZmdjoV/4M/sRUKhejq6qp1GWZmDUVSd0QUytvr9eS4mZnVqWEPVUn6QhXreDki/vsY1WNmZnVupHMcy4HPjNBnPeDgMDObJEYKjrsiYvNwHSS9ZQzrsRxs393Pxs5eDh4ZYGZLM+uWLmDF4orfmzQzG9Gw5zgi4vNDzZN0y0h9rPa27+5nw9Ye+o8MEED/kQE2bO1h+25//cXMRudMTo5/csyqsNxs7Oxl4PjJU9oGjp9kY2dvjSoys0Z3JsGhMavCcnPwyECmdjOzkZxJcEz8L4BMADNbmjO1m5mNZNjgkPQrSS+X/Bx8/QqYOU412hlYt3QBzVObTmlrntrEuqULalSRmTW6Ya+qiohzx6sQy8fg1VO+qsrMxspIXwDsBv4vxafw/TAi/mVcqrIxtWJxm4PCzMbMSOc4rgC2Ae8HHpX0PUl/JOm3cq/MzMzq0kiHqk4AP0wvJM0ElgF/Imke8FhE3JxzjWZmVkcy3VY9Ig4Cm4BNks4C3pNLVWZmVrdGOsfxbYa57DYirh3ziszMrK6NtMfx5+NShZmZNYyRznE8Ovhe0jRg8KR4b0Qcz7MwMzOrT1Wd45D0fmAz8AzFW43MlrQ6In6UX2lmZlaPqj05/hfAhyKiFyBdjvsN4PK8CjMzs/pU7b2qpg6GBkBE/ByYmk9JZmZWz6rd4+iS9FfAX6fpjwNd+ZRkZmb1rNo9jj8AngD+c3o9kdqGJWmZpF5JfZLWV5g/XdIDaf5OSe2pfYmkPen1uKSPlCzzjKSeNM/hZWY2zqra44iIV4DPpVdVJDUBdwNXAweAXZI6IuKJkm43Ai9FxDxJq4A7gI8Be4FCRJyQdBHwuKRvp2+yA3wgIl6othYzMxs7Ve1xSPpdSbslvVh6m/URFlsC9EXE/oh4FbgfWF7WZznFq7UAtgBXSVJE/LokJM7Gz/4wM6sb1R6q+jywGpgREedFxLkRcd4Iy7QBz5ZMH0htFfukoDgKzACQdIWkfUAPcFNJkATwsKRuSWuG+nBJayR1Seo6dOhQdaM0M7MRVRsczwJ7I2Lc/ucfETsj4hLgXcAGSWenWVdGxGXANcAfSnrfEMvfExGFiCi0traOU9VmZhNftVdV/Vfge5IeBV4ZbIyI4c559AOzS6ZnpbZKfQ5ImgKcDxwu7RART0o6BrwT6IqI/tT+vKRtFA+J+YuIZmbjpNo9jj8Ffk3xfMO5Ja/h7ALmS5qbbleyCugo69NB8RAYwEpgR0REWmYKgKS3Ae8AnpF0jqRzU/s5wIconkg3M7NxUu0ex8yIeGeWFacrotYCnUATsCki9km6jeKeQwdwL3CfpD7gRYrhAnAlsF7SceA14OaIeEHSxcA2SYO1fz0iHspSl5mZnRlVc9pC0p3ADyLi4fxLGnuFQiG6uvyVDzOzLCR1R0ShvD3LFwAfkjSQ4XJcMzObgKr9AuBI5zPMzGySGHaPQ9JvjrSCavqYmdnEMdKhqu9VsY5q+piZ2QQx0qGqfzPCuQwBPtdhZjaJjPTo2KbxKsTMzBpDtVdVmZmZAQ4OMzPLyMFhZmaZODjMzCyTUQeHpO+MZSFmZtYYzmSP4xNjVoWZmTWMah8de46ks0qmz6L4tD4zM5tkqt3j+DvgTSXTbwJ+MPblmJlZvas2OM6OiGODE+n9m4bpb2ZmE1S1wfHPki4bnJB0OTCQT0lmZlbPqn0C4C3A30g6SPH+VL8JfCy3qszMrG5V+zyOXZLeASxITb0RcTy/sszMrF5Vu8cB8C6gPS1zmSQi4mu5VGVWwfbd/Wzs7OXgkQFmtjSzbukCVixuq3VZZpNOVcEh6T7g7cAe4GRqDsDBYeNi++5+NmztYeB48dev/8gAG7b2ADg8zMZZtXscBWBhRESexZgNZWNn7+uhMWjg+Ek2dvY6OMzGWbVXVe2leEI8E0nLJPVK6pO0vsL86ZIeSPN3SmpP7Usk7UmvxyV9pNp12sR08Ejli/iGajez/FS7x3EB8ISknwCvDDZGxLVDLSCpCbgbuBo4AOyS1BERT5R0uxF4KSLmSVoF3EHxaq29QCEiTki6CHhc0rcpHh4baZ02Ac1saaa/QkjMbGmuQTVmk1u1wXHrKNa9BOiLiP0Aku4HlgOl/8gvL1n3FuCLkhQRvy7pczbFwKh2nTYBrVu64JRzHADNU5tYt3TBMEuZWR6qOlQVEY8C/wCcm15PprbhtAHPlkwfSG0V+0TECYr3v5oBIOkKSfuAHuCmNL+addoEtGJxG7dft4i2lmYEtLU0c/t1i3x+w6wGqr2q6t8BG4EfUvwC4F9KWhcRW/IqLCJ2ApdI+lfAZkkPZlle0hpgDcCcOXNyqNDG24rFbQ4KszpQ7aGqTwPviojnASS1UrzJ4XDB0Q/MLpmeldoq9TkgaQpwPnC4tENEPCnpGPDOKtc5uNw9wD0AhULBV4OZmY2Raq+qOmswNJLDVSy7C5gvaa6kacAqoKOsTwewOr1fCeyIiEjLTAGQ9DbgHcAzVa7TzMxyVO0ex0OSOoFvpOmPAcMeOkpXRK0FOoEmYFNE7JN0G9AVER3AvcB9kvqAFykGAcCVwHpJx4HXgJsj4gWASuuscgxmZjYGVO13+iRdR/EfdIC/j4htuVU1xgqFQnR1ddW6DDOzhiKpOyIK5e3D7nFImgdcGBE/joitwNbUfqWkt0fEU/mUa2Zm9Wqk8xSfB16u0H40zTMzs0lmpOC4MCJ6yhtTW3suFZmZWV0bKThahpnnez2YmU1CIwVHl6RPlDdK+k9Adz4lmZlZPRvpctxbgG2SPs4bQVEApgEfGXIpMzObsIYNjoj4JfDbkj5A8ZvbAN+NiB25V2ZmZnWp2meOPwI8knMtZmbWAKq95YiZmRng4DAzs4wcHGZmlomDw8zMMnFwmJlZJg4OMzPLxMFhZmaZODjMzCwTB4eZmWXi4DAzs0yqfea4mY2h7bv72djZy8EjA8xsaWbd0gWsWNxW67LMquLgMBtn23f3s2FrDwPHTwLQf2SADVuLz0tzeFgj8KEqs3G2sbP39dAYNHD8JBs7e2tUkVk2Dg6zcXbwyECmdrN6k2twSFomqVdSn6T1FeZPl/RAmr9TUntqv1pSt6Se9PODJcv8MK1zT3q9Nc8xmI21mS2Vn7o8VLtZvcktOCQ1AXcD1wALgeslLSzrdiPwUkTMA+4C7kjtLwC/FxGLgNXAfWXLfTwiLk2v5/Mag1ke1i1dQPPUplPamqc2sW7pghpVZJZNnnscS4C+iNgfEa8C9wPLy/osBzan91uAqyQpInZHxMHUvg9oljQ9x1rNxs2KxW3cft0i2lqaEdDW0szt1y3yiXFrGHleVdUGPFsyfQC4Yqg+EXFC0lFgBsU9jkEfBX4aEa+UtP1vSSeBvwX+JCKi/MMlrQHWAMyZM+cMh2I2tlYsbnNQWMOq65Pjki6hePjq90uaP54OYb03vf5jpWUj4p6IKEREobW1Nf9izcwmiTyDox+YXTI9K7VV7CNpCnA+cDhNzwK2ATdExFODC0REf/r5K+DrFA+JmZnZOMkzOHYB8yXNlTQNWAV0lPXpoHjyG2AlsCMiQlIL8F1gfUT8eLCzpCmSLkjvpwK/C+zNcQxmZlYmt+CIiBPAWqATeBL4ZkTsk3SbpGtTt3uBGZL6gE8Cg5fsrgXmAZ8pu+x2OtAp6WfAHop7LF/JawxmZnY6VTivPOEUCoXo6uqqdRlmZg1FUndEFMrb6/rkuJmZ1R8Hh5mZZeLgMDOzTBwcZmaWiYPDzMwycXCYmVkmDg4zM8vEwWFmZpn4meNmdka27+5nY2cvB48MMLOlmXVLF/jOvxOcg8PMRm377n42bO15/Rnq/UcG2LC1B8DhMYH5UJWZjdrGzt7XQ2PQwPGTbOzsrVFFNh4cHGY2agePDGRqt4nBwWFmozazpTlTu00MDg4zG7V1SxfQPLXplLbmqU2sW7qgRhXZePDJcTMbtcET4L6qanJxcJjZGVmxuM1BMcn4UJWZmWXi4DAzs0wcHGZmlomDw8zMMnFwmJlZJrkGh6Rlknol9UlaX2H+dEkPpPk7JbWn9qsldUvqST8/WLLM5am9T9IXJCnPMZiZ2alyCw5JTcDdwDXAQuB6SQvLut0IvBQR84C7gDtS+wvA70XEImA1cF/JMl8CPgHMT69leY3BzMxOl+cexxKgLyL2R8SrwP3A8rI+y4HN6f0W4CpJiojdEXEwte8DmtPeyUXAeRHxWEQE8DVgRY5jMDOzMnl+AbANeLZk+gBwxVB9IuKEpKPADIp7HIM+Cvw0Il6R1JbWU7rOit88krQGWAMwZ86cMxiGmU0WfrZIder65LikSygevvr9rMtGxD0RUYiIQmtr69gXZ2YTyuCzRfqPDBC88WyR7bv7a11a3ckzOPqB2SXTs1JbxT6SpgDnA4fT9CxgG3BDRDxV0n/WCOs0M8vMzxapXp7BsQuYL2mupGnAKqCjrE8HxZPfACuBHRERklqA7wLrI+LHg50j4jngZUnvTldT3QB8K8cxmNkk4WeLVC+34IiIE8BaoBN4EvhmROyTdJuka1O3e4EZkvqATwKDl+yuBeYBn5G0J73emubdDPwV0Ac8BTyY1xjMbPLws0Wqp+LFSRNboVCIrq6uWpdhZnWs/PnpUHy2yO3XLZq0J8gldUdEobzdt1U3M8PPFsnCwWFmlvjZItWp68txzcys/jg4zMwsEweHmZll4uAwM7NMHBxmZpaJr6oyM5tg8r5Zo4PDzGwCKf8i4+DNGoExCw8fqjIzm0DG42aNDg4zswlkPG7W6OAwM5tAxuNmjQ4OM7MJZN3SBTRPbTqlrXlqE+uWLhizz/DJcTOzCWQ8btbo4DAzm2DyvlmjD1WZmVkmDg4zM8vEwWFmZpk4OMzMLBMHh5mZZaKIqHUNuZN0CPjFKBe/AHhhDMuppYkylokyDvBY6tVEGcuZjuNtEdFa3jgpguNMSOqKiEKt6xgLE2UsE2Uc4LHUq4kylrzG4UNVZmaWiYPDzMwycXCM7J5aFzCGJspYJso4wGOpVxNlLLmMw+c4zMwsE+9xmJlZJg4OMzPLxMExBEnLJPVK6pO0vtb1ZCXpGUk9kvZI6kptvyHp+5L+Mf18S63rrETSJknPS9pb0laxdhV9IW2nn0m6rHaVn26IsdwqqT9tmz2SPlwyb0MaS6+kpbWp+nSSZkt6RNITkvZJ+qPU3nDbZZixNOJ2OVvSTyQ9nsby2dQ+V9LOVPMDkqal9ulpui/Nbx/VB0eEX2UvoAl4CrgYmAY8DiysdV0Zx/AMcEFZ253A+vR+PXBHrescovb3AZcBe0eqHfgw8CAg4N3AzlrXX8VYbgU+VaHvwvS7Nh2Ym34Hm2o9hlTbRcBl6f25wM9TvQ23XYYZSyNuFwFvTu+nAjvTn/c3gVWp/cvAH6T3NwNfTu9XAQ+M5nO9x1HZEqAvIvZHxKvA/cDyGtc0FpYDm9P7zcCKGtYypIj4EfBiWfNQtS8HvhZFjwEtki4an0pHNsRYhrIcuD8iXomIp4E+ir+LNRcRz0XET9P7XwFPAm004HYZZixDqeftEhFxLE1OTa8APghsSe3l22Vwe20BrpKkrJ/r4KisDXi2ZPoAw/9i1aMAHpbULWlNarswIp5L7/8JuLA2pY3KULU36rZamw7hbCo5ZNgQY0mHNxZT/N9tQ2+XsrFAA24XSU2S9gDPA9+nuEd0JCJOpC6l9b4+ljT/KDAj62c6OCauKyPiMuAa4A8lva90ZhT3VRvyWuxGrj35EvB24FLgOeAvaltO9SS9Gfhb4JaIeLl0XqNtlwpjacjtEhEnI+JSYBbFPaF35P2ZDo7K+oHZJdOzUlvDiIj+9PN5YBvFX6hfDh4uSD+fr12FmQ1Ve8Ntq4j4ZfrL/hrwFd447FHXY5E0leI/tP8nIram5obcLpXG0qjbZVBEHAEeAd5D8dDg4KPBS+t9fSxp/vnA4ayf5eCobBcwP12ZMI3iSaSOGtdUNUnnSDp38D3wIWAvxTGsTt1WA9+qTYWjMlTtHcAN6SqedwNHSw6d1KWyY/0fobhtoDiWVenKl7nAfOAn411fJek4+L3AkxHxuZJZDbddhhpLg26XVkkt6X0zcDXFczaPACtTt/LtMri9VgI70p5iNrW+KqBeXxSvCvk5xeOFn651PRlrv5jiVSCPA/sG66d4LPPvgH8EfgD8Rq1rHaL+b1A8VHCc4vHZG4eqneJVJXen7dQDFGpdfxVjuS/V+rP0F/mikv6fTmPpBa6pdf0ldV1J8TDUz4A96fXhRtwuw4ylEbfLvwZ2p5r3Ap9J7RdTDLc+4G+A6an97DTdl+ZfPJrP9S1HzMwsEx+qMjOzTBwcZmaWiYPDzMwycXCYmVkmDg4zM8vEwWFmZpk4OMzKSJpRcmvtfyq51fYxSf8rh8/7qqSnJd2UcbnvDX75a5g+G9MYPnVmVZq9YcrIXcwml4g4TPF+RUi6FTgWEX+e88eui4gtI3d7Q0R8uIo+6yT98+jLMjud9zjMqiTp/ZK+k97fKmmzpL+X9AtJ10m6U8WHZz2U7oWEpMslPZruUtxZza3F0x7IlyQ9Jml/+txNkp6U9NWSfs9IukBSe5r3lfQwn4fT7SfMcuHgMBu9t1N87sG1wF8Dj0TEImAA+LcpPP4SWBkRlwObgD+tct1voXizuv9C8fYXdwGXAIskXVqh/3zg7oi4BDgCfHTUozIbgQ9VmY3egxFxXFIPxadGPpTae4B2YAHwTuD76Vk5TRTvW1WNb0dEpHX/MiJ6ACTtS+veU9b/6YgYbOtOfcxy4eAwG71XACLiNUnH440bv71G8e+WgH0R8Z7Rrjut65WS9sF1D9Uf4CTgQ1WWGx+qMstPL9Aq6T1QfAaEpEtqXJPZGXNwmOUkis+rXwncIelxioeXfru2VZmdOd9W3azG0pVS38l6OW6G9d/K+FxSbJOE9zjMau8o8D+yfgGwGpI2Av8B8Hc5bMx4j8PMzDLxHoeZmWXi4DAzs0wcHGZmlomDw8zMMvn/zcyV08oftBMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data=np.array([[0,0.05],\n",
    "               [50,.038],\n",
    "               [100,.0306],\n",
    "               [150,.0256],\n",
    "               [200,.0222],\n",
    "               [250,.0195],\n",
    "               [300,.0174]])\n",
    "\n",
    "plt.plot(data[:,0],data[:,1],'o')\n",
    "plt.xlabel('Time [min]')\n",
    "plt.ylabel('Conc. [mol/L]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oCafkMlO2Tmy"
   },
   "source": [
    "So, we need to convert the list of numbers to a numpy array so we can do the analysis.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 297
    },
    "executionInfo": {
     "elapsed": 582,
     "status": "ok",
     "timestamp": 1651148064962,
     "user": {
      "displayName": "Zachary Ulissi",
      "userId": "07633171379186475882"
     },
     "user_tz": 240
    },
    "id": "kVYnxN6K2AVN",
    "outputId": "60f804e2-eff2-4ff7-9df9-70d7aae60e49"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f57372fcc50>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU9bnH8c+TDSIgVEArm0FZREABY7RVXLAC2ipUaQFthdYr16q9tb3SirYWqftSvVXbXlwqLlRaBEUFsRZE7EUgCBYxjSJSATdEgoKBbM/945yEIUySGTKTyfJ9v17zysw5vznzHCfm4bebuyMiIhKrtFQHICIiTYsSh4iIxEWJQ0RE4qLEISIicVHiEBGRuGSkOoCG0KlTJ8/JyUl1GCIiTcqqVas+dffO1Y+3iMSRk5NDfn5+qsMQEWlSzOzf0Y6rqUpEROKixCEiInFR4hARkbi0iD4OEWm+SktL2bx5M7t37051KE1W69at6datG5mZmTGVV+IQkSZt8+bNtGvXjpycHMws1eE0Oe7Otm3b2Lx5Mz179ozpPUltqjKzkWZWaGbrzeyaKOdbmdms8PxyM8sJj+eYWbGZrQkff4x4z/FmtjZ8z+8sSb8pT6/ewsm3LqLnNc9z8q2LeHr1lmR8jIjU0+7du+nYsaOSxgEyMzp27BhXjS1pNQ4zSwfuB84CNgMrzWyeu78VUewSYLu79zKzccBtwNjw3LvuPijKpf8AXAosB+YDI4EFiYz96dVbmDJnLf3KChiVXsBrO/oxZU4JAKMHd03kR4lIAihp1E+8//2SWePIA9a7+wZ3LwGeBEZVKzMKmBE+nw2cWVsNwswOBw5299c8WA/+UWB0ogO/Y2Eh/coKeCLrZn6W8VeeyLqZfmUF3LGwMNEfJSLS5CQzcXQFNkW83hwei1rG3cuAHUDH8FxPM1ttZkvMbGhE+c11XBMAM5tkZvlmlr9169a4Av+gqJiT0grIopQMqyCTMk5KK+CDouK4riMiLUN6ejqDBg2qetx66631vubUqVPp2rUrgwYNonfv3px//vm89dZbdb7vkUce4YMPPqj359emsXaOfwj0cPdtZnY88LSZ9Y/nAu4+HZgOkJubG9duVV06ZPPajn6UkU46ZTjGaxX96NIhO57LiEgj9PTqLdyxsJAPiorp0iGbySP61rsJOjs7mzVr1iQowr1++tOfcvXVVwMwa9Yshg0bxtq1a+nceb9VQKo88sgjDBgwgC5duiQ8nkrJrHFsAbpHvO4WHotaxswygPbANnff4+7bANx9FfAu0Ccs362Oa9bb5BF9Kcjox/iSX1JQ0Y1y0tiWcRiTR/RN9EeJSAOq7L/cUlSMA1uKipkyZ21SBr/s2LGDvn37UlgYNHGPHz+eBx54AIC2bdvy05/+lP79+3PmmWcSS6vI2LFjGT58ODNnzgRg2rRpnHDCCQwYMIBJkybh7syePZv8/HwuuugiBg0aRHFxcdRy9ZXMxLES6G1mPc0sCxgHzKtWZh4wIXw+Bljk7m5mncPOdczsSKA3sMHdPwQ+N7OTwr6Qi4FnEh346MFdueX8gXzc/jgmlV5NmjmP5bygjnGRJu6OhYUUl5bvc6y4tLze/ZfFxcX7NFXNmjWL9u3bc9999zFx4kSefPJJtm/fzqWXXgrArl27yM3NZd26dZx22mnccMMNMX3OkCFD+Ne//gXAlVdeycqVK3nzzTcpLi7mueeeY8yYMeTm5vLEE0+wZs0asrOzo5arr6Q1Vbl7mZldCSwE0oGH3X2dmU0D8t19HvAQ8JiZrQc+I0guAKcC08ysFKgALnP3z8JzlwOPANkEo6kSOqKq0ujBXfcmir9tpMey+6BoE3ToXvsbRaTRqqmfsr79lzU1VZ111ln89a9/5YorruCNN96oOp6WlsbYscEA0u9973ucf/75MX1OZG1h8eLF3H777Xz55Zd89tln9O/fn3PPPXe/98RaLh5J7eNw9/kEQ2Yjj10f8Xw38J0o73sKeKqGa+YDAxIbaR2G/jcc+10lDZEmrkuHbLZESRLJ6r+sqKigoKCAgw46iO3bt9OtW7eo5WIdDrt69Wpyc3PZvXs3l19+Ofn5+XTv3p2pU6dGnYcRa7l4aa2qWLQ+GA4L++Z3f57aWETkgE0e0ZfszPR9jmVnpiet//Luu++mX79+zJw5kx/84AeUlpYCQUKZPXs2ADNnzuSUU06p81pPPfUUL774IuPHj6/649+pUyd27txZdS2Adu3a8cUXXwDUWq4+GuuoqsZp0U3wz1lwxQrIbJ3qaEQkTpXNz4keVVXZx1Fp5MiR/OAHP+DBBx9kxYoVtGvXjlNPPZUbb7yRG264gTZt2rBixQpuvPFGDj30UGbNmhX1unfffTePP/44u3btYsCAASxatKhqRNWll17KgAED+OpXv8oJJ5xQ9Z6JEydy2WWXkZ2dzbJly2osVx+WiB72xi43N9cTspHThiXw6Hlw5q9h6M/qfz0RqbeCggL69euX6jDi0rZtW3bu3JnqMPYR7b+jma1y99zqZVXjiMeRp0Hfc2DJ7bDnC+h7NnTPS3VUIiINSn0c8RrwHSgrhld/CzPOg00rUh2RiDQxja22ES8ljngVvQeEIyDKS2Dj0pSGIyLS0JQ44pUzFDJagaVDelbwWkSkBVEfR7y658GEZ4OaRo+T4aCOdb9HRKQZUY3jQHTPCyYFFjwDD54Ju7alOiIRkQajxFEfQy4OJgT+PbZ1ZkSkebrpppvo378/xx57LIMGDWL58uW1lp86dSp33nknANdffz0vvfRSvWPYuHEj2dnZ+6yZ9f777zNmzBgA1qxZw/z58+u4SmzUVFUfh/aDk34Ey+4Pkki3/YY7i0gzt2zZMp577jlef/11WrVqxaeffkpJSUnM7582bVrCYjnqqKP2WzOrcrb4mjVryM/P55xzzqn356jGUV+nXwPtvgrP/wwqyusuLyKpt2kFLL0rIcPpP/zwQzp16kSrVq2AYHmPyr0wcnJy+PnPf87AgQPJy8tj/fr1+71/4sSJVX/cc3Jy+PWvf82QIUMYOHBg1Uq4u3bt4oc//CF5eXkMHjyYZ56JbVHwjRs3MmDAAEpKSrj++uuZNWtW1eq99aEaR321agcjboKXpsKOTfCVnFRHJNKy/emb+x/rPxryLoWSL+HhEfDxm+AVYGlw2AA48TIYfFHQX/mXi/d97w+er/Xjhg8fzrRp0+jTpw/f+MY3GDt2LKeddlrV+fbt27N27VoeffRRrrrqqjqXNe/UqROvv/46v//977nzzjt58MEHuemmmxg2bBgPP/wwRUVF5OXl8Y1vfIM2bdrs89533323aumTk08+mcmTJwOQlZXFtGnTyM/P57777qv182OhGkci9D8frlgJOz9J2L9iRCRJdu8IkgYEP3fvqNfl2rZty6pVq5g+fTqdO3dm7NixPPLII1Xnx48fX/Vz2bJldV6vcon1448/no0bNwLw4osvcuuttzJo0CBOP/10du/ezfvvv7/feyubqtasWcP9999fr/uqjWociWAGH/0TZpwLZSXBPI8J87QciUgq1FZDyDoILngwWPWhvCSYi3XBg3v/X23Tsc4aRjTp6emcfvrpnH766QwcOJAZM2YwceJEYN8l02NZPr2yySs9PZ2ysjIg2Ifjqaeeom/fxrELqWocibJxKZTtASqgfI9mlIs0Vt3zgn/YDbsuIf/AKyws5J133ql6vWbNGo444oiq15X9CbNmzeJrX/vaAX3GiBEjuPfee6s2clq9enXc14hcbr2+VONIlMoZ5WW7wR26qbYh0mh1z0tYi8DOnTv58Y9/TFFRERkZGfTq1Yvp06dXnd++fTvHHnssrVq14s9//vMBfcavfvUrrrrqKo499lgqKiro2bNn3FvAnnHGGVXNXVOmTKnagfBAaFn1RNq0AlY9AmuegFMnw7BfJv8zRVq4xrysek5ODvn5+XTq1CnVodRJy6qnSuW/YrwCXr0b+n97786BIiLNhBJHMoy4GYq3Q5r+84q0ZJWjopob/WVLhoMOgQvrN8FGRGLn7jGNWJLo4u2ySOqoKjMbaWaFZrbezK6Jcr6Vmc0Kzy83s5xq53uY2U4zuzri2EYzW2tma8ysATou6mHXNnjmCijaf7y1iCRG69at2bZtW9x//CTg7mzbto3WrVvH/J6k1TjMLB24HzgL2AysNLN57v5WRLFLgO3u3svMxgG3AZFd/b8FFkS5/Bnu/mmSQk+c0l3w5lz44mO46K/BfA8RSahu3bqxefNmtm7dmupQmqzWrVvTrVu3mMsns6kqD1jv7hsAzOxJYBQQmThGAVPD57OB+8zM3N3NbDTwHrAriTEmV4cecOav4IVr4OXbICMzGLariYEiCZOZmUnPnj1THUaLksymqq7ApojXm8NjUcu4exmwA+hoZm2BXwDR1it34EUzW2VmkxIedaLlTYJOfWHJLbDoRu1TLiJNXmOdOT4VuNvdo+3ofoq7DwHOBq4ws1OjXcDMJplZvpnlp7QKm5YOR50RPPcK7VMuIk1eMhPHFqB7xOtu4bGoZcwsA2gPbANOBG43s43AVcC1ZnYlgLtvCX9+AswlaBLbj7tPd/dcd8/t3Llzou7pwAy4IFgTR/uUi0gzkMw+jpVAbzPrSZAgxgEXViszD5gALAPGAIs8GBpR9ZfVzKYCO939PjNrA6S5+xfh8+FA4nZBSZbueTDx+b37lB/aOGe5iojEImmJw93LwlrCQiAdeNjd15nZNCDf3ecBDwGPmdl64DOC5FKbw4C54XjtDGCmu7+QrHtIqO550O0EePyCYE2rcTM1ykpEmiStVdXQlt0PC6+F0X+EQeNTHY2ISI1qWquqsXaON18n/giOOBkW/AJ2bE51NCIicVPiaGhpaTDqfqgog2euDJZgFxFpQpQ4UuGQnjD8N/DZBihcoO1mRaRJ0SKHqZL7Q+jYC2aO3buFpbabFZEmQDWOVDGDLflB0vByTQwUkSZDiSOVcoYGM8shSCSaGCgiTYASRypVTgw8fFDQWV5ekuqIRETqpMSRapXJ45CjYM4k+PKzVEckIlIrJY7GoFVbGPMQ7PwEFvw81dGIiNRKo6oaiy6Dg/kdXQanOhIRkVopcTQmx4WbH7rDns+hdfvUxiMiEoUSR2O04Bfwzt/guHHBXh6a2yEijYj6OBqj9t1h+wZ4+RbtGCgijY4SR2NUUQIY4FC2RxMDRaRRUeJojHKGBnt2AFABHXJSGY2IyD6UOBqj7nkw4Vk4+So4uCu06ZTqiEREqqhzvLHqnhc8zrx+77IkIiKNgGocjV1aejA899V7YNWMVEcjIqIaR5PgFfDeEtj4Khx+rCYJikhKqcbRFKSlw/kPQpvOMHMcLLpRQ3RFJGWUOJqKNh3h1F/Azo/glTs0v0NEUkaJoykp/pRgfgdQrvkdIpIaSU0cZjbSzArNbL2ZXRPlfCszmxWeX25mOdXO9zCznWZ2dazXbNZyhkJGa7A0SG+ljZ9EJCWSljjMLB24HzgbOAYYb2bHVCt2CbDd3XsBdwO3VTv/W2BBnNdsvrrnBfuSD/tl8LN4O3y0NtVRiUgLk8xRVXnAenffAGBmTwKjgLciyowCpobPZwP3mZm5u5vZaOA9YFec12zeKud3lO6Ge48Pah+TFmuSoIg0mGQ2VXUFNkW83hwei1rG3cuAHUBHM2sL/AK44QCuCYCZTTKzfDPL37p16wHfRKOV2RrGPgo7P4a/TIDy0lRHJCItRGPtHJ8K3O3uOw/0Au4+3d1z3T23c+fOiYusMel6PIy6D/79KrzQsrp7RCR1ktlUtQXoHvG6W3gsWpnNZpYBtAe2AScCY8zsdqADUGFmu4FVMVyzZTn2u0E/x//9DspLYPD3tX+HiCRVMmscK4HeZtbTzLKAccC8amXmARPC52OARR4Y6u457p4D3APc7O73xXjNlqfvOZCeBauf0PwOEUm6pCWOsM/iSmAhUAD8xd3Xmdk0MzsvLPYQQZ/GeuBnQK3tLTVdM1n30GS8/39QUQ5eHtQ61s1NdUQi0oyZu6c6hqTLzc31/Pz8VIeRPJtWBDWN8pJgXas2h8JlS6HdYamOTESaMDNb5e651Y831s5xiUfV/I7r4Lx7oeQL+PNYKNlV93tFROJUa+e4mf0uhmt87u6/TFA8cqAq53cAHNQRnrwQnroUxj6m/TxEJKHqqnGMIhjJVNvjgmQGKAfg6HPg7Nug8HlY+WCqoxGRZqau4bh3u3utuweZ2VcSGI8kyon/CdmHwDHnseTvz1OwbD4v7urFx+2PY/KIvoweHHXepIhInWpNHO5+T03nzOwqd7+ntjKSYsd+hyV/f54TX5nAKZQxISuTi3Zcy5Q5JQBKHiJyQOrTOf6zhEUhSVOwbD5ZlJJuThalnJRWQHFpOXcsLEx1aCLSRNUncVjCopCkeXFXL/aQiTuk4Xzi7QH4oKg4xZGJSFNVn8TR/CeANAMftz+Oi0qu4w9l57KV9kzJ/DNH2gd06ZCd6tBEpImqazjuFwQJwtg3URigvzxNwOQRfZkyp4TXS/vwl4rTmZ11Azdn/YmPRvw11aGJSBNVV+d4u4YKRJKjsgP8joWF/LvocH7W6gbGnzFEHeMicsDqqnGsAl4l2IXvZXff3SBRSUKNHtx1/0RRXgbzr4a2h0GvM7WirojErK4+jhOBucDpwBIzm29mPzGzPkmPTJJr9eOw6k+w5FaYca5W1BWRmNXVVFUGvBw+MLMuwEjgRjPrBbzm7pcnOUZJhuJtVHVdle2Gdxer1iEiMYlrVJW7f+DuD7v7d4Fc4InkhCVJlzMUMlpT9Suw/iUo25PSkESkaairj+NZahl26+7n1XROGrnKFXU3LoXdO2DlQ7C1EA4/NtWRiUgjV9daVXc2SBSSGpEr6n7tx9A23JvdHUzzO0Ukurr6OJZUPg+3aq3sFC9099JkBiYNrDJpLJ8OhfPhiJPhyNPU7yEi+4mpj8PMTgfeAe4Hfg+8bWanJjEuSZVPC2HDYlh8o/YvF5Go6mqqqnQXMNzdCwHC4bh/Bo5PVmCSIgd3Ye9oq2LY8LJqHSKyj1hHVWVWJg0Ad38byExOSJJSVaOtwj6O9X+HioqUhiQijUusNY58M3sQeDx8fRGQn5yQJKUiR1sVb4d2h0OatqYXkb1iTRw/Aq4A/it8vZSgr6NWZjYS+B8gHXjQ3W+tdr4V8ChBk9c2YKy7bzSzPGB6ZTFgqrvPDd+zEfgCKAfK3D03xnuQWEWOtqr0wWro2AtaafkykZYupsTh7nuA34aPmJhZOkFn+lnAZmClmc1z97ciil0CbHf3XmY2DrgNGAu8CeS6e5mZHQ68YWbPhjPZAc5w909jjUXqafcOeHR0UPvody70Pkv9HiItWKyjqr5lZqvN7DMz+9zMvjCzz+t4Wx6w3t03uHsJ8CQwqlqZUUDlnuazgTPNzNz9y4gk0Rrt/ZFardvDyT+BrQXwyu1a20qkhYu18foeYALQ0d0Pdvd27n5wHe/pCmyKeL05PBa1TJgodgAdAczsRDNbB6wFLotIJA68aGarzGxSTR9uZpPMLN/M8rdu3RrbXUotnKpfl7LdUPBsSqMRkdSJNXFsAt509wb7l7+7L3f3/sAJwBQzax2eOsXdhwBnA1fUNJ/E3ae7e66753bu3LmBom7GcoZCRiuw8FemaFPt5UWk2Yq1c/znwHwzWwJUrYTn7rX1eWwBuke87hYei1Zms5llAO0JOsmruHuBme0EBgD57r4lPP6Jmc0laBJ7Jcb7kAMVOdrqkCPh6G8FxyvKIS09tbGJSIOKtcZxE/AlQX9Du4hHbVYCvc2sZ7hcyThgXrUy8wiawADGAIvc3cP3ZACY2RHA0cBGM2tjZu3C422A4QQd6dIQuufB0P+G/t+G9EzYuRX+OFTNViItTKw1ji7uPiCeC4cjoq4EFhIMx33Y3deZ2TSCmsM84CHgMTNbD3xGkFwATgGuMbNSoAK43N0/NbMjgbkWLMCXAcx09xfiiUsSKC0dsg6CWd+Ho78ZdKBrtJVIs2exdFuY2e3AS+7+YvJDSrzc3FzPz9d8xaTYsAQeGw1eAWkZMHE+9Dgx1VGJSAKY2apoc+Vibar6EfCCmRXHMRxXWoIt+VQtT1JRBi/fktJwRCT5Yp0AqOnCEl3OUEjPgvKSYMTVSdpJWKS5q7XGYWZfresCsZSRZqxytNWw6+AH86HPcCgvhXk/hk/Xpzo6EUmCumoc84EhCSgjzVn1ta22b4R/PR+Mthr2K9hdFNRM1HEu0izUlTiOq6MvwwD1dci+OvWG/3gJHjkXnv8ZYMFS7RPmKXmINAO1NlW5e3q4xEhNj3buXn0ZEZFgkuCg8eELD5Yp2bg0pSGJSGJoowVJnt7D924KlZ4VNFeJSJMX6wRAkfh1z4MJzwY1jco+jlfuDJYrOfToVEcnIgdIiUOSK7LjfOdWWP6/8Oo9cNrkYN6HOs1FmhwlDmk4bTvDpJfh0VHwt+tRp7lI03TAfRxm9lwiA5EWon1XGPid8IVDWbE6zUWamPp0jl+asCikZTnqDMjIRp3mIk1TTE1V4RLmxe5eEb5OI9itTyR+kXt7VPZxrHokWLJkyMWpjk5E6hBrH8ffgW8AO8PXBwEvAl9PRlDSAkR2mrsHM83feRH+9Rx0OT6olajfQ6RRirWpqrW7VyYNwucHJSckaXHMYPyTcNyF8PZCePlmmPEt2LQi1ZGJSBSxJo5dZla1HpWZHQ8UJyckaZHS0qFTL6p+Jcv2wNvao0ukMYq1qeoq4K9m9gHB+lRfBcYmLSppmXKGQkarvUu09xkZHNe+5iKNSqz7caw0s6OBvuGhQncvTV5Y0iJF6zTf+A+Yf3WwLe3nWzRhUKQRiGcC4AlATvieIWaGuz+alKik5aq+RDsOX3wIc/8Tx9hDJhfuuZaP2x/H5BF9GT1Ya2yKNLSY+jjM7DHgTuAUggRyArDfPrQiCZdzCuT+EAcMJ8tLOCNtDVuKipkyZy1Pr96S6ghFWpxYaxy5wDHu7skMRiSqPiPZs/ReMr2UNJzWtgeA4tJy7lhYqFqHSAOLdVTVmwQd4nExs5FmVmhm683smijnW5nZrPD8cjPLCY/nmdma8PGGmX071mtKM9Q9jwv3XMtdZd/lx6VXcnPZRQAcYxv5tEjzUEUaWqw1jk7AW2a2AthTedDdz6vpDWaWDtwPnAVsBlaa2Tx3fyui2CXAdnfvZWbjgNsIRmu9CeS6e5mZHQ68YWbPAh7DNaUZ+rj9cfy+qE/V67Z8yeNZN/NZWkd4owQ+36yOc5EGEmvimHoA184D1rv7BgAzexIYBUT+kR8Vce3ZwH1mZu7+ZUSZ1gQJI9ZrSjM0eURfpsxZS3FpOQA7OYgpFVfyu6zfw9xJaKVdkYYTU1OVuy8B/gW0Cx8F4bHadAU2RbzeHB6LWsbdywjWv+oIYGYnmtk6YC1wWXg+lmtKMzR6cFduOX8gXTtkY0DXDtmcff7FtDqpcq3NcKXddxelMkyRFiHWRQ6/C9wBvEwwAfBeM5vs7rOTFZi7Lwf6m1k/YIaZLYjn/WY2CZgE0KNHjyREKA1t9OCu+3eEbxoJy34fzDQ3g6OGpSY4kRYk1qaq64AT3P0TADPrDLxE0LxUky1A94jX3cJj0cpsNrMMoD2wLbKAuxeY2U5gQIzXrHzfdGA6QG5urkaDNVfRJg1uLQw2ihr8ffi0UH0fIgkWa+JIq0waoW3U3cy1EuhtZj0J/riPAy6sVmYeMAFYBowBFrm7h+/ZFHaOHwEcDWwEimK4prQ01ScNfvo2vPtyuNaV+j5EEi3W4bgvmNlCM5toZhOB54Fam47CPokrgYVAAfAXd19nZtPMrHI01kNARzNbD/wMqBxeewrBSKo1wFzgcnf/tKZrxnqz0kL0OxfyqvV9FDyb0pBEmhOLdU6fmZ1P8AcdYKm7z01aVAmWm5vr+fn5qQ5DGtKmFTDjPCjbDTgcOxbOn57qqESaFDNb5e77rRJSa1OVmfUCDnP3f7j7HGBOePwUMzvK3d9NTrgi9RTZ99Gp795O8835sP09KHpffR8iB6iuPo57gClRju8Iz52b8IhEEqV634c7PH150GGOBUu4T3hWyUMkTnX1cRzm7murHwyP5SQlIpFkMYN+3wpfeNCMtfrxlIYk0hTVlTg61HIuO5GBiDSIPiMhIzvYKAqD12fAhrrmsopIpLqaqvLN7FJ3fyDyoJn9B7AqeWGJJElk30fXE2BrQdDXAfDmHPjsPeipvg+R2tSVOK4C5prZRexNFLlAFvDtGt8l0phF9n0ceWrw8+2/wewfBM8zWsGE55Q8RGpQa1OVu3/s7l8HbiCYgLcRuMHdv+buHyU/PJEG8tEbBKvpECxfsvCXUFyU0pBEGqtY9xxfDCxOciwiqdPz1GCGeXlJ8HrzcrgvF360DNp2Tm1sIo1MPHuOizRf1de8Ss8Klixp2zmYTFi4APqereYrEZQ4RPaqPu+jy6BwBvq3guarf/wPjHsiSCAiLVisa1WJtEwbl0J5WfDcy2HW9+GVO6C0OLVxiaSQEodIbSqbrSw9GG3VPQ8W3Qh/+HpQC9m0ApbeFfwUaSHUVCVSm2j7fby3FD4pgA/fCBZSLN8D6a20dLu0GEocInWp3vfRc2jwWHpXUOugIli6/c05ShzSIihxiByonKGQkRUmD4flf4SSnXD6FPh8y761FJFmRIlD5EB1zwtW1924FA47FjYsghUPQOELQQIpLwn6R9SEJc2MEodIfUQ2Y/U5C078T1hyO7zxZDAKq2w3rH9JiUOaFY2qEkmkr+TA8RPDkVhpgMOy38Or90DJrhQHJ5IYqnGIJFrkSKw2neGtZ+ClX8Oy+2DABXBQJ5aUHcO1K7P5oKiYLh2ymTyiL6MHd0115CIxUeIQSYbIJqwhF8P7r8H8n8Py/6UCI88zOKzkWrbQhy1FxUyZE+yXpuQhTYGaqkQaQo+ToP8oMCONCjIp44HMu7gsfR5tKKa4tJw7FhamOkqRmKjGIdJQcoZCeivKSvdQRjrv+6Fck/kkl2U8y4LyPD75ogNsaquOdGn0klrjMDsjHAoAABJISURBVLORZlZoZuvN7Joo51uZ2azw/HIzywmPn2Vmq8xsbfhzWMR7Xg6vuSZ8HJrMexBJmLDv48HMC7mw5Dq+XfobztvzG96u6Mb4jMX8V+ZcmHGuli+RRi9picPM0oH7gbOBY4DxZnZMtWKXANvdvRdwN3BbePxT4Fx3HwhMAB6r9r6L3H1Q+PgkWfcgknDd8/jqN6+lIKMfAP/0o3i5YhDlbsE2UuWlQaf66ifgk3+lNFSRmiSzqSoPWO/uGwDM7ElgFPBWRJlRwNTw+WzgPjMzd18dUWYdkG1mrdx9TxLjFWkQlR3gdyws5IOiYta3GYSXPQNeGgzj7XoC/OV7sHsHdD8JDj0ajrsQepyY4shFAslMHF2BTRGvNwPVf/Oryrh7mZntADoS1DgqXQC8Xi1p/MnMyoGngBvd3at/uJlNAiYB9OjRo563IpJYowd33XcE1abcfZco+fHqYAjv6sdg02uwagac97tghJZIijXqznEz60/QfDU84vBF7r7FzNoRJI7vA49Wf6+7TwemA+Tm5u6XWEQaleoLKbbpCIf0DJZz93LAYWs46uqz96B1e9i2XuthSUokM3FsAbpHvO4WHotWZrOZZQDtgW0AZtYNmAtc7O7vVr7B3beEP78ws5kETWL7JQ6RJq9yL5DKNa+OGRUcf2EKvPt3qAgTipZ0lwaWzFFVK4HeZtbTzLKAccC8amXmEXR+A4wBFrm7m1kH4HngGnf/R2VhM8sws07h80zgW8CbSbwHkdSpnIE+7Lp9E8OZ18Nh/YOaiIdLuuf/KbWxSouStBpH2GdxJbAQSAcedvd1ZjYNyHf3ecBDwGNmth74jCC5AFwJ9AKuN7Prw2PDgV3AwjBppAMvAQ8k6x5EUq56ExbAYcfA2bcHQ3fLwq6/9MzgZ3lpsKjiJ2+pCUuSxqL0Kzc7ubm5np+fn+owRBJr04qgj6PH16HLIMjMhlfuhEW/Cc6nZ8HFz8IRJ6U2TmmyzGyVu+dWP96oO8dFpBbRaiO7dwAGeNA38uQ4OOWncOKP4MM16kyXhFDiEGlO+p0bbCZVXgJpaXBwV3j90WA+yKOjtD+6JIQSh0hzErmke2XNYvfnsDJMJpWd6YtvhlH3Q3utxivxU+IQaW6qN2G1Pjgc2psJZR60ZG1YDPcMgKPODIb57vpETVgSMyUOkZYgcn/0nKHBBlNrngiG8b73ClSUQXoGfPO3MOgiMEt1xNKIaVSVSEv2yp1Bs5WXU9Wpfmh/GDQeOvaGT9apJtKCaVSViOyv56lB8igvCZqy8i6Ffy+DF3+5t0xGtjrTZR9KHCItWbTOdAiWNXntD1QN631vKRTMg97D4YhTYEu+hva2YEocIi1dtPkg/b8d9H9UrpPV8Sh45gr4v3vhoM6we3swQiu9FUu+9hDXrszmg6JiunTIZvKIvto7vZlT4hCR/UWrifQeDoXzYfEt8OVWACrK9rBqyTy2lJwLGFuKipkyZy2AkkczpsQhItFVr4lkHQQDx0CHHlXrZJWQwaulR/Nf6XM5Of1NXq/oRUlFJvMXbGD04MtTF7sklRKHiMQnYmjvhfPTeN370IfNdOVTTsz4F+5QsedpeGUXnDo51dFKEiRzWXURaa6658HQ/+bj9scB8GT5MGaWDwv2TjdIM4e1s/eWXz4dXr4tWJhRmjwlDhE5YJNH9CU7Mx2A1yqOoYRMyjyN8rTWMPymoNBbz8CCyfDyzfDwSPi/+4Ll36XJUlOViBywyg7wOxYWsrqoDz9pdQM/Oepj+n/9m3v7R7a+TfBv1IpgouGL18HSO2Hs48GILQ3rbXKUOESkXkYP7lr7CKojT4Old+2dZHjGtfDRm1DyJfxlDJTtBkuDkbdA7iXB0ifSqOkbEpHkqmmSYWUywYOayIKfw5Lbg6XhB5wPGa1VG2mklDhEJPmiTTLMGRo0VVXWRE6dDB+vg3/Ogvdfg+0bw31FMuDip+GIr6ckdNmfEoeIpEZNNZGSL2HxTcGSJ14O5eXw2Leh//nBEvCt2sKm5aqJpJASh4ikTrSaSNZBQYJY+VBY40gPksS/noc3ZoaF0iBDOxmmihKHiDQ+0WojZSXw7FXwxhNARZBUVj4U9JUc/S1o+1X4+J+qiTSApCYOMxsJ/A+QDjzo7rdWO98KeBQ4HtgGjHX3jWZ2FnArkAWUAJPdfVH4nuOBR4BsYD7wE28Jm4qItDTVayMZWZA7EdbN2bv44ldygg2p3n5hb7m0zCDpqE8kaZI2AdDM0oH7gbOBY4DxZnZMtWKXANvdvRdwN3BbePxT4Fx3HwhMAB6LeM8fgEuB3uFjZLLuQUQamcqayLDrgp9nTIGr1sIJlxJsRAVUlAZ7igAUPAuvPwav3KVZ6wmUzBpHHrDe3TcAmNmTwCjgrYgyo4Cp4fPZwH1mZu6+OqLMOiA7rJ0cAhzs7q+F13wUGA0sSOJ9iEhjUr0mYgbHfhdWP753hFbPoeAO834MxdvDcunwjRsg7z8gMzs1sTcTyUwcXYFNEa83AyfWVMbdy8xsB9CRoMZR6QLgdXffY2Zdw+tEXjPqzCMzmwRMAujRo0c9bkNEGr2aRmjlXhL0gVTOFfnbL+HzLXD2rUGt5J0Xoe/ZVeWfXr2FOxYWam+ROjTqznEz60/QfDU83ve6+3RgOgR7jic4NBFpbKKN0OozApbdv7cmMux66DUsaLZ6dBSU74FX74ZB3+Plg7/FtYsr+LLUGWJvc9IXBcyaMwC4QMmjmmQmji1A94jX3cJj0cpsNrMMoD1BJzlm1g2YC1zs7u9GlO9WxzVFRAK1zVqvqFxo0WHNY5zOYyxJa8+0tO9xe+YDZFJGKXP5yYIM7S1STTJXx10J9DaznmaWBYwD5lUrM4+g8xtgDLDI3d3MOgDPA9e4+z8qC7v7h8DnZnaSmRlwMfBMEu9BRJq6cAn4fWojOUMhvVXQ75GRDRfO5qcll/NKxUBy7GMyKSPDKmhNCRcVz4QP/xn0mWxaESSdFt7RnrQaR9hncSWwkGA47sPuvs7MpgH57j4PeAh4zMzWA58RJBeAK4FewPVmdn14bLi7fwJczt7huAtQx7iIxCtKTWTFwenMLTqFIfY2lzMP81IMOC19LfzvUMj+Cuz5HJxgKHALnnxoLWEKRG5urufn56c6DBFpxJ5evYUpc9ZSXFoe9HGkFbA6bQAXn3MqZ2e/Fewj8sm6oLClQ+e+Qcf6UWcGxza91uwmH5rZKnfPrX68UXeOi4g0lOp7i3zc7jgmj+jL2YO7AoOgYy+Ycd7ejnZLg1fvCUdthdKzYOLzzSp5RKPEISISqnVvkWgd7bt3wPNXw9q/Ag7lZcH5Np3g5VuhQ49gt8Ojv9mskokSh4hIrKoP+W3dHvIuDWaoVy6DkjMUPnsPChcEfSIA//gf6HsOnHMHtO8adK434b1GlDhEROqjpiG/J/8kWB7eKwCH9S8Fm1NtWgGPnBPUTtIzYewT0CfuqWoplczhuCIiLUO0Ib89T913yO/Fz0CbjkGCKS8jaNoqgZnfgftPhAXXNJkhv6pxiIgkQ001kZyhQc2jcnfDwRdB0fuw7R3YvDLogC8rDhLOyT+Bwd+DQ44MzjWS5i0NxxURaWg19XEsvQsW3RSsqxWpdQco2RnUSNKz4PtzGmTZeA3HFRFpLKKtqwXV9mHPglH3QsmuYMOqj9ZS1bw141zomgvdT4DsQ2DPFw26WKMSh4hIY1FT89ahx+w7h+Tob0LRpmBf9oqyoMyye2HifOZvyWb28y/RvayUUenreW1HP6bMKQFIWPJQ4hARaUyi1UZqSihLboeXbwlGblVUwMalvP7Kdh5PvxdPgwqMEjK5qORa7liYlbDEoVFVIiJNQbSRW0eevnfkVjiHZPYXA3mq7JRgSS1zMinjpLQCPigqTlgoqnGIiDRVUWoibTrs5Ikd3+Cc9BVkehmlZPBaRT+6dEjcrodKHCIiTVm1pq3JI/oyZU4JF5Vcy0lpBbxW0Y+CjH7cMqJvwj5SiUNEpBnZu1hjFn8o6kOXDtncolFVIiJSm1oXa0wAdY6LiEhclDhERCQuShwiIhIXJQ4REYmLEoeIiMSlRayOa2ZbgX8f4Ns7AZ8mMJxUai730lzuA3QvjVVzuZf63scR7t65+sEWkTjqw8zyoy0r3BQ1l3tpLvcBupfGqrncS7LuQ01VIiISFyUOERGJixJH3aanOoAEai730lzuA3QvjVVzuZek3If6OEREJC6qcYiISFyUOEREJC5KHDUws5FmVmhm683smlTHEy8z22hma81sjZnlh8cOMbO/mdk74c+vpDrOaMzsYTP7xMzejDgWNXYL/C78nv5pZkNSF/n+ariXqWa2Jfxu1pjZORHnpoT3UmhmI1IT9f7MrLuZLTazt8xsnZn9JDze5L6XWu6lKX4vrc1shZm9Ed7LDeHxnma2PIx5lpllhcdbha/Xh+dzDuiD3V2Pag8gHXgXOBLIAt4Ajkl1XHHew0agU7VjtwPXhM+vAW5LdZw1xH4qMAR4s67YgXOABYABJwHLUx1/DPcyFbg6Stljwt+1VkDP8HcwPdX3EMZ2ODAkfN4OeDuMt8l9L7XcS1P8XgxoGz7PBJaH/73/AowLj/8R+FH4/HLgj+HzccCsA/lc1TiiywPWu/sGdy8BngRGpTimRBgFzAifzwBGpzCWGrn7K8Bn1Q7XFPso4FEPvAZ0MLPDGybSutVwLzUZBTzp7nvc/T1gPcHvYsq5+4fu/nr4/AugAOhKE/xearmXmjTm78XdfWf4MjN8ODAMmB0er/69VH5fs4Ezzczi/Vwljui6ApsiXm+m9l+sxsiBF81slZlNCo8d5u4fhs8/Ag5LTWgHpKbYm+p3dWXYhPNwRJNhk7iXsHljMMG/bpv091LtXqAJfi9mlm5ma4BPgL8R1IiK3L0sLBIZb9W9hOd3AB3j/UwljubrFHcfApwNXGFmp0ae9KCu2iTHYjfl2EN/AI4CBgEfAnelNpzYmVlb4CngKnf/PPJcU/teotxLk/xe3L3c3QcB3QhqQkcn+zOVOKLbAnSPeN0tPNZkuPuW8OcnwFyCX6iPK5sLwp+fpC7CuNUUe5P7rtz94/B/9grgAfY2ezTqezGzTII/tE+4+5zwcJP8XqLdS1P9Xiq5exGwGPgaQdNg5dbgkfFW3Ut4vj2wLd7PUuKIbiXQOxyZkEXQiTQvxTHFzMzamFm7yufAcOBNgnuYEBabADyTmggPSE2xzwMuDkfxnATsiGg6aZSqtfV/m+C7geBexoUjX3oCvYEVDR1fNGE7+ENAgbv/NuJUk/tearqXJvq9dDazDuHzbOAsgj6bxcCYsFj176Xy+xoDLAprivFJ9aiAxvogGBXyNkF74XWpjifO2I8kGAXyBrCuMn6Ctsy/A+8ALwGHpDrWGuL/M0FTQSlB++wlNcVOMKrk/vB7Wgvkpjr+GO7lsTDWf4b/Ix8eUf668F4KgbNTHX9EXKcQNEP9E1gTPs5pit9LLffSFL+XY4HVYcxvAteHx48kSG7rgb8CrcLjrcPX68PzRx7I52rJERERiYuaqkREJC5KHCIiEhclDhERiYsSh4iIxEWJQ0RE4qLEISIicVHiEKnGzDpGLK39UcRS2zvN7PdJ+LxHzOw9M7sszvfNr5z8VUuZO8J7uLp+UYrslVF3EZGWxd23EaxXhJlNBXa6+51J/tjJ7j677mJ7ufs5MZSZbGa7Djwskf2pxiESIzM73cyeC59PNbMZZrbUzP5tZueb2e0WbJ71QrgWEmZ2vJktCVcpXhjL0uJhDeQPZvaamW0IP/dhMysws0ciym00s05mlhOeeyDczOfFcPkJkaRQ4hA5cEcR7HtwHvA4sNjdBwLFwDfD5HEvMMbdjwceBm6K8dpfIVis7qcEy1/cDfQHBprZoCjlewP3u3t/oAi44IDvSqQOaqoSOXAL3L3UzNYS7Br5Qnh8LZAD9AUGAH8L98pJJ1i3KhbPuruH1/7Y3dcCmNm68NprqpV/z90rj60Ky4gkhRKHyIHbA+DuFWZW6nsXfqsg+H/LgHXu/rUDvXZ4rT0RxyuvXVN5gHJATVWSNGqqEkmeQqCzmX0Ngj0gzKx/imMSqTclDpEk8WC/+jHAbWb2BkHz0tdTG5VI/WlZdZEUC0dKPRfvcNw4rj+VhhlSLC2EahwiqbcD+E28EwBjYWZ3AN8DNJdDEkY1DhERiYtqHCIiEhclDhERiYsSh4iIxEWJQ0RE4vL//zVCZtIXa74AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import interpolate\n",
    "\n",
    "# data will be a 2d list, which we convert to an array here\n",
    "data = np.array(data)\n",
    "t = data[:, 0]   # column 0\n",
    "Ca = data[:, 1]  # column 1\n",
    "\n",
    "# calculate a spline through the data\n",
    "tck = interpolate.splrep(t, Ca)\n",
    "\n",
    "t_eval = np.linspace(0,300)\n",
    "Ca_spline = interpolate.splev(t_eval, tck)\n",
    "plt.plot(data[:,0],data[:,1],'o', label='Exp Data')\n",
    "plt.plot(t_eval, Ca_spline,'--.',label='Spline Fit')\n",
    "plt.xlabel('Time [min]')\n",
    "plt.ylabel('Conc. [mol/L]')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 375
    },
    "executionInfo": {
     "elapsed": 1620,
     "status": "ok",
     "timestamp": 1651148121933,
     "user": {
      "displayName": "Zachary Ulissi",
      "userId": "07633171379186475882"
     },
     "user_tz": 240
    },
    "id": "WVNctlap2Tmz",
    "outputId": "f9494984-cb10-4e46-8976-d34dbe15bf17"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.7/dist-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.\n",
      "  import pandas.util.testing as tm\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha = 2.0354816446001127, conf interval [1.92418422 2.14677907]\n",
      "k = 0.1402128334966632, conf interval [0.09372748 0.20975319]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZfrG8e+TECAoEKQpSLOAXdGIXVEQsIKrKPbO6ro/RBZU7GJZlLXg2mB11VV3WVcR6QEUsSMgKlJFBCUoRQygBEh5fn+cAQMOMElm5kyS+3NduZhyJnOfi5Cbc+Y972vujoiISHmlhR1AREQqBxWKiIjEhQpFRETiQoUiIiJxoUIREZG4UKGIiEhcVAs7wI6Y2WHAs0BNoBD4k7t/GmW7y4E7Infvd/eXdva9GzRo4C1btoxjWhGRym/GjBmr3L1htOcsla9DMbMJwGPuPs7MTgdudvf222yzGzAdyAYcmAEc4e4/7+h7Z2dn+/Tp0xMTXESkkjKzGe6eHe25VD/l5UCdyO26wLIo23QGJrr76kiJTAS6JCmfiIhEpPQpL6A3kGNmfyMov2OjbNMU+L7E/aWRx0REJIlCLxQzmwTsHuWp24EOwE3u/oaZnQ88D3Qsx3v1BHoCNG/evKzfRkREogi9UNx9uwVhZv8Cbozc/R/wXJTNcoH2Je7vCby7nfcaCgyF4DOU0qcVEZHtSfXPUJYBJ0VunwJ8HWWbHKCTmdUzs3pAp8hjIiKSRKEfoezEtcBgM6sGbCByusrMsoHr3P0ad19tZvcB0yKvGeDuq8OJKyJSdaX0sOFE0rBhEamSfl0FXgy7NirTyyvysGEREYkHd/j83/BkNoztl5C3UKGIiFR2qxfBv7rCiOv5ckNjTp15LMcNfIcRM3Pj+jap/hmKiIiUVVEBfPwkvDuQAqrxYPHVvLjpZJw0yMun//BZAHRrG59L91QoIiKVUe4MGHkjLJ8F+53Jed9244tNtbbaJL+giEE58+NWKDrlJSJSmWz8Bcb3h+c6wvpVcMEr0ONVvlxTK+rmy/Ly4/bWOkIREaksFkyAMX1gzfeQfTV0vBtq1gWgSVYmuVHKo0lWZtzeXkcoIiIV3S8r4PWr4N/dofoucFUOnPnoljIB6Ne5DZkZ6Vu9LDMjnX6d28Qtho5QREQqKneY+QpMuAMK1kP72+D43lCtxu823fw5yaCc+SzLy6dJVib9OreJ2+cnoEIREamYVi2E0b1h8fvQ/Fg4azA0bL3Dl3Rr2zSuBbItFYqISEVSuAk+GgxTBkG1mkGRtL0M0sL/BEOFIiJSUXw/DUb1ghVz4IBucNpDUDva6h/hUKGIiKS6jevg7QHw6T+gThO4cBi0OS3sVL+jQhERSWXzxsLYvrB2GbTrCR3uhBq1w04VlQpFRCQVrfsRxt0Mc96CRgdA95eg2ZFhp9ohFYqISCopLobPXoKJd0PhBjjlTjjuRkjPCDvZTqlQRERSxcoFMOpG+O4jaHlCMIKr/t5hp4qZCkVEJGyFm+CDx+D9v0FGLTj7SWh7CZiFnaxUVCgiImH67pPgqGTlPDjoXOgysMyrKYZNhSIiEoYNa2DSvTD9eajbDC76H7TuFHaqclGhiIgk29xRwTK8vyyHo/8EJ98ONXYNO1W5qVBERJJl7bKgSOaNhsYHQ49XoekRYaeKGxWKiEiiFRcHp7Ym3QvFBdDxXjjmhgoxFLg0VCgiIom0Ym7wofv3U2Gv9nDmY7DbXmGnSggViohIIhRsgPcfCYYD16gN3Z6FQ3tUuKHApaFCERGJt8UfBkclP30Nh1wAnR+EXRqEnSrhVCgiIvGSnwcT7wqmTslqDpe8Aft0DDtV0qhQRETKyx3mjIBxt8CvK+HY/4P2/YP13asQFYqISHmsWQpj+sKCcbDHoXDRa9DksLBThUKFIiJSFsVFMO25YOErL4ZO98NR10N61f21WnX3XESkrJbPhpG9IHc67N0BznwU6rUMO1XoVCgiIrEq2ADvPQwfDoaadeEP/4CDu1fqocCloUIREYnFt+/BqN6w+hs49CLo/ADU2i3sVClFhSIisiPrV8PEO2HmK8FprUtHwN4nh50qJalQRESicYev3oDxtwalclxvOOkWqF4r7GQpS4UiIrKtvO9gdB9YOBGaHA6Xvgm7Hxx2qpSnQhER2ay4CKY+C+/cD1iwemK7npCWHnayCiGlC8XMDgOeBWoChcCf3P3TKNsVAbMid79z97OTl1JEKoUfvoRRvWDZTNi3E5zxSDB9isQspQsFeBi4193Hmdnpkfvto2yX7+5V89JUESmfTethykD46Mlg1NZ5/4QD/6ChwGWQ6oXiQJ3I7brAshCziEhl881kGN0bfl4MbS+FUwdoKHA5pHqh9AZyzOxvQBpw7Ha2q2lm0wlOiw109xHJCigiFdCvP8GE2+GL/8Bue8Plo6HVCWGnqvBCLxQzmwTsHuWp24EOwE3u/oaZnQ88D0SbC7qFu+ea2V7AO2Y2y92/ifJePYGeAM2b69yoSJXjDl++Bjn9YcMaOKEvnNgPMmqGnaxSMHcPO8N2mdkaIMvd3cwMWOPudXbymheB0e7++o62y87O9unTp8cvrIiktp8XB0OBv3kbmmbD2U9A4wPDTlXhmNkMd8+O9lxassOU0jLgpMjtU4Cvt93AzOqZWY3I7QbAccCcpCUUkdRWVAgfPgFPHR2s637aILh6gsokAUI/5bUT1wKDzawasIHI6Sozywauc/drgP2BIWZWTFCQA91dhSIisOxzGPl/8OOX0OZ0OH0Q1N0z7FSVVkoXirt/ABwR5fHpwDWR2x8BuoRVRH6z6VeY/CB88jTs0hC6vwQHdNVQ4ARL6UIREdnWiJm5DMqZz7K8fJpkZdKvcxu6tW362wYLJ8Hom4LpU464AjreC5lZoeWtSlQoIlJhjJiZS//hs8gvKAIgNy+f/sODSTK67Vs9GL01639Qf1+4Yiy0PC7MuFWOCkVEKoxBOfO3lMlm+QWFfDX2GbpVewU2/gIn3Qon9IFqNUJKWXWpUESkwliWl7/V/Rb2Iw9Ue57jC2bD7kfDWYOh0X4hpZNUHzYsIrJFk6xMAKpRyPXpI8mpfguHpC1iULU/wpXjVCYhU6GISIXRr3Mbjsz4lpHV7+SWjGFMLj6Ms4ofZd8zboQ0/ToLm055iUjFsPEXuv34d7qmD2EVWfxx0018VefE34/yktCoUEQk9S2YAGP6wJql2JFX07DDXQypWTfsVLINFYqIpK5fVsC4W2D2cGi4H1yVA82PCjuVbIcKRURSjzvMfBkm3AEF+XDy7XBcb6hWPexksgMqFBFJLasWBoteLX4fmh8bDAVu2DrsVBIDFYqIpIbCTfDRYJgyCKrVDIqk7WUavVWBqFBEJHzfT4NRvWDFHDigG5z2ENSOtu6epDIVioiEZ8NaeHsATHsO6jSBC4dBm9PCTiVlpEIRkXDMGwtj/gLrfoCj/gin3AE1aoedSspBhSIiybXuRxjbD+aOhEYHwgUvw55RV5SVCkaFIiLJUVwMn70EE++Gwg3Q4S44thekZ4SdTOJEhSIiibdyAYy6Eb77CFqeEIzgqr932KkkzlQoIpI4hRvhg8fg/UcgoxZ0fQoOu1hL8VZSKhQRSYwlHwdHJavmw0HnQZe/wq6Nwk4lCaRCEZH42rAGJt0D0/8JdZvBRf+D1p3CTiVJoEIRkfiZMzIYwfXrCjj6Bjj5Nqixa9ipJElUKCJSfmuXBUUybzQ0Phgu/Dc0PSLsVJJkKhQRKbviYpj+PEy6F4oLoOO9cMwNGgpcRalQRKRsVsyFkb1g6aewV3s48zHYba+wU0mIVCgiUjoFG4JhwB88FkyV0u1ZOLSHhgKLCkVESmHxBzCqN/z0NRzSAzo/ALs0CDuVpIhSF4qZ7QJscPeiBOQRkVSU/zNMvAs++xdktYBLhsM+HcJOJSlmp4ViZmlAD+Bi4EhgI1DDzFYBY4Ah7r4woSlFJBzuMPvNYF339T8Fc2+1vxWq7xJ2MklBsRyhTAYmAf2Br9y9GMDMdgNOBh4yszfd/ZXExRSRpMv7Hsb2hQXjYY9D4ZLXgz9FtiOWQuno7gVm1nJzmQC4+2rgDeANM9MYQZHKorgIPv0HvHMfeDF0egCOug7S9ZGr7NhOf0LcvSByczhweMnnzOxod/+kxDYiUpH9+FWwFG/uDNi7A5z5KNRrGXYqqSBi+QzlfIIiqW1m+wPzSxypDAUOSWA+EUmGgnyY8jB89ATUzII/PAcHn6ehwFIqsRzDfgjUBK4BHgXamFkesAzIT2A2EUmGRVNgdG9YvSiYWr7T/VBrt7BTSQUUS6EcSvCh/Nnu/hGAmdUHWgLzEhdNRBJq/WqYcCd8/grUawWXvRVc8S5SRrEUyjnAAKCxmc0DvgA+j/y5IYHZRCQR3GHW6zD+1uD6kuNvgpNugYzMsJNJBZe2sw3c/Vp3zwaeARYAiwiGC08FliQynJkdamYfm9ksMxtlZnW2s10XM5tvZgvN7NZEZhKp0H5eAq+eB8Ovgazm8Mcp0PEelYnERWnGAV7g7lsGoZvZ00C/+EfaynNAX3efYmZXRd7vzpIbmFk68BRwKrAUmGZmI919ToKziVQcRYUw9VmY/ABg0GUgtOsJaelhJ5NKZKdHKCWsNbMtCxy4+wygdfwjbaU18F7k9kTg3CjbtAMWuvsid98EDAO6JjiXSMXxwxfwXAeYcDu0PAFumApHX68ykbgrzRHK1cBwM5sGzAAOBhJ9/clsgnIYAXQHmkXZpinwfYn7S4GjEpxLJPVtWg/v/hU+fgpq1YfzXoADz9FQYEmYmAvF3ReY2eFAN4IymQvcVt4AZjYJ2D3KU7cDVwFPmNmdwEhgUznfqyfQE6B58+bl+VYioRoxM5dBOfNZlpdPk6xM+nVuQ7e2TX/b4Jt3glmB85bA4ZfBqQMgs154gaVKiOXCxmOATzywCXgt8hUX7t5xJ5t0iuRoDZwR5flctj5y2TPyWLT3GkpwMSbZ2dle6rAiKWDEzFz6D59FfkEw4XduXj79h88CoFvrmpBzG3w5DOrvA5ePhlYnhBlXqpBYPkO5DJhhZsPM7Aozi3Y0kRBm1ijyZxpwB/BslM2mAfuaWSszq04wM/LIZGUUSbZBOfO3lMlm+QWFfDl2CDyZDV+9Dif2g+s+VJlIUsUyl9f1AGa2H3Aa8KKZ1SWYhXg88GEC10a50MxuiNweDrwQydIEeM7dT3f3QjP7M5ADpAP/dPfZCcojErpleVtPUNHMlvNAtX9yYsEsaHwknPUEND4gpHRSlZl76c/8mFkmwbUopwHHRK5TqVCys7N9+vTpYccQKbXjBr5Dbl4+6RRxTfpYeld7g0LSGZpxCX+57SGN3pKEMrMZ2/udX5phw1u4e767jyW4LuRP5QknIqXTr3MbsjMWM7L6HfTP+A8fFB/MWcWPsPcZN6lMJFRlWQK4D8EoryzgAOAV4NM45xKRaDb+QrflT9E1/Rl+oi7Xb+rNl7VPpF+X/bYe5SUSgrKsmNMOGOXur5rZM+5+X7xDiUgUX0+E0X1gzXfYEVfSoOM9PJOZFXYqkS1KXSju3sPMzjazlwEtLC2SaL+sDCZy/Op1aNAarhwHLY4NO5XI75T1M5SRBOujfGZm/4hvJBEBglmBZ74SDAWe8xacdCtc94HKRFJWLBc23gmsdvenSj7u7huB+80s2nQoIlIeP30TLHr17XvQ7Gg4azA02i/sVCI7FMspr/OBI7d90MyuARq6+1/jnkqkqioqCJbhnfIwpFeHMx+Dw6+AtDKdTBBJqlgKpcDdoy2k9TLwGaBCEYmHpdNhZC9YMRv2PwtOGwR19gg7lUjMYimUTWa2h7v/UPJBd99oZomebVik8tu4Dt65H6YOgdp7wAWvwv5nhp1KpNRiKZRHgLfMrLu7b1mhMTLPliZYFCmP+eNhzF9gbS4ceQ10uAtqRl2YVCTlxTKX1//MrBbBBJGfEKwnn0awPsk9iY0nUkmtWw7jboY5I6Dh/nBVDjTXMj5SscUyyut0gtUSWxGshXIQ8CtwkbtPS2w8kUqmuBhm/gsm3gUF+XDyHXDcjVCtetjJRMotllNe5wADgMbAPOALYBGwwczSEzjTsEjlsuprGHUjLPkQWhwXDAVusG/YqUTiJpZTXtcCmNltBMvtLiKYaXgosJpgQSsR2Z7CTfDh4/DeIMjIDKaXb3uphgJLpVOaqVcucPdDN98xs6cJZhsWke35/tNgKPDKucF67l0egtqNw04lkhClKZS1ZnaEu88AcPcZkWV5RWRbG9bC2/fCtOehThO4cBi0OS3sVCIJVZpCuRoYbmbTgBkEU9jrOhSRbc0dDWP7wbof4Kg/wil3QI3aYacSSbiYC8XdF5jZ4QQjvQ4G5gK3JSqYSIWz9gcY1w/mjoJGB8IFL8OeFW4xU5EyK9X09e6+yczy3f3ORAUSqXCKi2HGCzDpHijcGFyceGwvSM8IO5lIUpVlga0HgFHxDiJSIa2cHwwF/u5jaHUinPk41N877FQioShLoVjcU4hUNIUb4f1H4f1HoPou0PVpOOwiMP3zkKqrLIWi+bukalvyMYzqBasWwEHnQZeBsGvDsFOJhK4shSJSNeXnBZ+TzHgB6jaHi1+HfU8NO5VIylChiOyMO8wdCWNvhl9XwDF/hvb9ocauYScTSSllKZTlcU8hkqrW5AbXlMwfA7sfDBf+B5oeHnYqkZRU6kJxdx3jS+VXXBRc5f72ACguhFMHwNF/0lBgkR2IZfp6c/cdfhAfyzYiFcbyOcGH7kunwV7tg3Xdd9sr7FQiKS+WI5TJZvYG8Ja7f7f5QTOrDhwPXA5MBl5MSEKRZCnYEMwI/OHjUKMOnDMEDrlAQ4FFYhRLoXQBrgL+Y2atgDygJpAOTAAed/eZiYsokgSLPwguUPxpIRzSAzo/ALs0CDuVSIUSy3ooG4CngafNLANoAOS7e16iw4nE04iZuQzKmc+yvHyaZGXSr3Mbuu1XCybcCTNfhqwWcMlw2KdD2FFFKqTSzuVVAPyQoCwiCTNiZi79h88ivyBYYDQ3bz3vvTmELpmvUHNTXjD3Vvtbg6veRaRMYi4UM6sBnAu0LPk6dx8Q/1gi8TUoZ/6WMmnCKu7LeIEOaTOZt2kv9uv5Juxx6E6+g4jsTGmOUN4C1hCshbIxMXFEEmNZXj5pFHN5eg59q72GAfcVXMxLRV1YqDIRiYvSFMqe7t4lYUlEEuiEOsvps+FJDkv7hilFh3B74VUs9UY0zcoMO5pIpVGaQvnIzA5291kJSyMSbwX5MOUhXiz4Oz9bLXptuoGRxccCRmZGOv06twk7oUilUZpCOR640swWEZzyMsDd/ZCEJBMpr0Xvwqje8PO3pB12MZ/scQMzJi/HSo7yats07JQilUZpCqULkRJJUJbfMbNDgWeBXYHFwMXuvjbKdouBdUARUOjuWne1Klu/GibcAZ+/CvVawWUjYa+TOAM446gDw04nUmnFMvXKOqKXyOZyqRPvUCU8B/R19ylmdhXQD9je8sMnu/uqBGaRVOcOs16H8bfChjw4vg+cdDNk6HMSkWSI5cLG2skIsh2tgfcitycCOWy/UKQq+3kJjOkDCydBk8Ph7Ldg94PCTiVSpaT6eiizga7ACKA70Gw72zkwwcwcGOLuQ5OUT8JWVAhTn4XJDwAGXR6CdtdCWnrYyUSqnNALxcwmAbtHeep2gjnEnjCzO4GRwKbtfJvj3T3XzBoBE81snru/t+1GZtYT6AnQvHnzuOSXEP3wBYzsBT98Dvt2hjMegazt/Z9DRBIt9EJx94472aQTgJm1Bs7YzvfIjfy5wszeBNrx26myktsNBYYCZGdna7r9imrTenj3Qfj4aahVH857AQ48R7MCi4Qs9ELZETNrFCmJNOAOghFf226zC5Dm7usitzsBmg6mslr4Noy+CfKWwOGXBQtfZdYLO5WIAGlhB9iJC81sATAPWAa8AGBmTcxsbGSbxsAHZvYF8Ckwxt3Hh5JWEufXVTC8J7zyh2DVxCvGwNl/V5mIpJCUPkJx98HA4CiPLwNOj9xeBGgypsrKHb4YBjm3wca1cGI/OKEvZNQMO5mIbCOlC0WquNXfBqe3Fk2GPdvBWYOh8QFhpxKR7VChSOopKoSPn4R3B0JaNTj9b5B9NaSl+hlakapNhSKpJfczGNULfpwFbc6A0wdBXc23JVIRqFAkNWz8BSY/CFOfgV0awfkvw/5naSiwSAWiQpHwfT0RRveBNd/BEVdCx3sgMyvsVCJSSioUCc8vK4OJHL96HRq0gSvHQ4tjwk4lImWkQpHkcw+mls+5HQrWQ/v+cPxNUK1G2MlEpBxUKJJcP30Do3vDt+9B82OCocANtWqiSGWgQpHkKCqAj56AKQ9DenU48zE4/AoNBRapRFQoknhLpwezAq+YDfufDac9DHX2CDuViMSZCkUSZ+M6eOd+mDoEau8BPf4N+0WdMFpEKgEViiTG/PEw5i+wNheOvAY63AU1E7latIiETYUi8bVuOYy7GeaMgIb7w9UToFm7sFOJSBKoUCQuRnz2PfPGPs31BS9S0wr4Zr9eHHDenVCtetjRRCRJNMRGym3Sex/Q9K3u3Fr4NHO9BV02DuTc2ccxYtbKsKOJSBLpCEXKrnATfPg4J05+iHyqc3PBtbxW1B4wKChiUM58urXVxI4iVYUKRcrm+0+DocAr55JTdDQDCi5jJVvPv7UsLz+kcCISBhWKlM6GtfD2vTDteajTFC56jYHDM1gZpTyaZGWGEFBEwqLPUCR2c0fDU0cFZXLUdXDDVGjdmX6d25CZkb7VppkZ6fTrrClVRKoSHaHIzq39Acb1g7mjoPFBcMErsOcRW57e/DnJoJz5LMvLp0lWJv06t9HnJyJVjApFtq+4GGa8AJPugaJN0OFuOPb/ID3jd5t2a9tUBSJSxalQJLoV82DUjfD9J9DqRDjzcai/d9ipRCSFqVBka4Ub4f1H4f1HoMau0PVpOOwiLcUrIjulQpHfLPkoOCpZtQAO7g6d/wq7Ngw7lYhUECoUgfw8mHQ3zHgR6jaHi1+HfU8NO5WIVDAqlKrMHeaOhLE3w68r4Jg/w8m3QfVdwk4mIhWQCqWqWpMLY/vC/LGw+yFw0TBo0jbsVCJSgalQqpriouDCxLcHQHEhnDoAjr4B0vWjICLlo98iVcnyOTCqFyydBnudHKzrvlursFOJSCWhQqkKCjbAe4Pgw8ehRh04ZwgccoGGAotIXKlQKrvFHwRDgX9aCIf0gM4Pwi71w04lIpWQCqWCGTEzN7Y5s9avhol3wcyXIasFXPom7H1K8gOLSJWhQqlARszMpf/wWeQXFAGQm5dP/+GzgN8maMQdZg+HcbcEpXJsL2jfH6rXCiu2iFQRKpQKZFDO/C1lsll+yZUR876HMX+Br3Ngj8Pgkjdgj0NDSisiVY0KpQLZ3gqIP+b9Cp88A2/fBzh0eiBYr0RDgUUkifQbpwJpkpVJ7jalsr8t4ZGaz8P4hbBPRzjjUajXIqSEIlKVpcSKjWbW3cxmm1mxmWVv81x/M1toZvPNrPN2Xt/KzKZGtvuvmVVPTvLkKrkyYg02cXO1YYyqfjt7Z6yGc58P5uBSmYhISFLlCOUr4A/AkJIPmtkBQA/gQKAJMMnMWrt70Tavfwh4zN2HmdmzwNXAM4mPnVybP3h/Z+z/6LPxaVqmLWdJs260uPAxqLVbyOlEpKpLiUJx97kA9vsL7boCw9x9I/CtmS0E2gEfb97AghedAlwUeegl4B4qYaGwfjXdFt9Pt4J/Q/1WcNZIWux1UtipRESAFCmUHWgKfFLi/tLIYyXVB/LcvXAH21Rs7jDrdRh/K2zIg+P7wEk3Q0Zm2MlERLZIWqGY2SRg9yhP3e7ubyUpQ0+gJ0Dz5s2T8Zbl9/MSGNMHFk6CpkfAWW/B7geFnUpE5HeSViju3rEML8sFmpW4v2fksZJ+ArLMrFrkKCXaNpszDAWGAmRnZ3sZ8iRPUSFMfQYmPwgYdHkI2l0LaelhJxMRiSolRnntwEigh5nVMLNWwL7ApyU3cHcHJgPnRR66HEjKEU/C/PAFPNcBJtwBrU6EG6bC0depTEQkpaVEoZjZOWa2FDgGGGNmOQDuPht4DZgDjAdu2DzCy8zGmlmTyLe4BegT+dC+PvB8svchLjatD0pk6Mmwdhmc9wJcOAyymu38tSIiIbPgP/hVT3Z2tk+fPj3sGL9Z+DaMvgnylsDhlwULX2XWCzuViMhWzGyGu2dHey7VR3lVfr+ugpzb4Mv/Qv194Iox0PL4sFOJiJSaCiUs7vDFsKBMNq6DE2+GE/4CGTXDTiYiUiYqlDCsXhSc3lr0LuzZDs5+AhrtH3YqEZFyUaEkU1EBfPwUvDsQ0qrB6X+D7KshLSXGRoiIlIsKJVlyP4ORvWD5LGhzBpw+COpWrgv6RaRqU6Ek2sZfgosTpz4DuzSC81+GA84OO5WISNypUBLp64kwug+s+Q6yr4KO90DNumGnEhFJCBVKIvyyIpjI8as3oEEbuHI8tDgm7FQiIgmlQoknd/j8Vci5HQrWQ/vb4PjeUK1G2MlERBJOhRIvP30Do26Exe9D82PgrMHQsE3YqUREkkaFUl5FBfDhYJjyMFSrCWc+DodfrqHAIlLlqFDKY+n0YCjwitlwQFc47WGoHW3JFxGRyk+FUhYb18Hb98GnQ6H2HtDjP7Df6WGnEhEJlQqltHJnwH8vg7W5wYJXp9wJNeuEnUpEJHQqlNKq2wzqtYDuL0CzdmGnERFJGfrkuJRGfF3Accv70uqplRw38B1GzIy62rCISJWjI5RSGDEzl/7DZ5FfUARAbl4+/YfPAqBbW83LJSJVm45QSmFQzvwtZbJZfkERg3Lmh5RIRCR1qFBKYVlefqkeFxGpSlQopdAkK7NUj4uIVCUqlFLo17kNmRnpWz2WmZFOv86aYkVERB/Kl8LmD94H5cxnWV4+TbIy6de5jT6QFxFBhVJq3do2VYGIiEShU14iIhIXKhQREYkLFYqIiMSFCkVEROJChSIiInFh7h52hlCY2UpgSQhv3QBYFcL7Jpv2s/KoCqXCDTgAAAXGSURBVPsI2s9YtXD3htGeqLKFEhYzm+7u2WHnSDTtZ+VRFfYRtJ/xoFNeIiISFyoUERGJCxVK8g0NO0CSaD8rj6qwj6D9LDd9hiIiInGhIxQREYkLFYqIiMSFCiVJzOwvZuZm1mA7z19uZl9Hvi5Pdr7yMrP7zOxLM/vczCaYWZPtbPewmc02s7lm9oSZWbKzlkcp9rN55Pm5ZjbHzFomN2nZxbqPkW3rmNlSM3symRnjIZb9NLPDzOzjyM/sl2Z2QRhZy6MUP7Pl/x3k7vpK8BfQDMghuJCyQZTndwMWRf6sF7ldL+zcpdzHOiVu9wKejbLNscCHQHrk62OgfdjZ472fkefeBU6N3N4VqBV29njvY+T5wcC/gSfDzp2I/QRaA/tGbjcBfgCyws6egP2My+8gHaEkx2PAzcD2RkB0Bia6+2p3/xmYCHRJVrh4cPe1Je7uQvR9daAmUB2oAWQAyxOfLn5i2U8zOwCo5u4TI6/5xd3XJyliucX4d4mZHQE0BiYkI1e8xbKf7r7A3b+O3F4GrACiXiWeqmL8+4zL7yAtsJVgZtYVyHX3L3Zwdqcp8H2J+0sjj1UoZvYAcBmwBjh52+fd/WMzm0zwvzwj+F/t3OSmLL+d7SfB/2rzzGw40AqYBNzq7kXJS1k+O9tHM0sDHgEuATomN138xPB3WXLbdgT/GfomCdHiKob9jMvvIB2hxIGZTTKzr6J8dQVuA+4KO2M87GQ/cffb3b0Z8Crw5yiv3wfYH9iT4If1FDM7IZn7EIvy7ifBf9ROAPoCRwJ7AVckKX5M4rCPfwLGuvvSZOYurTjs5+bvswfwMnCluxcnJ33s4rWf5Rb2+b3K/AUcTHCIvDjyVQh8B+y+zXYXAkNK3B8CXBh2/nLsd3PgqyiP9wPuLHH/LuDmsPMmYD+PBqaUuH8p8FTYeeO8j69GfpYXE0w0uBYYGHbeeO9n5Lk6wGfAeWHnTODfZ1x+B+kIJYHcfZa7N3L3lu7ekuAw8nB3/3GbTXOATmZWz8zqAZ0ij1UYZrZvibtdgXlRNvsOOMnMqplZBnASUKFOecW4n9OALDPbfK79FGBOorPFSyz76O4Xu3vzyM91X+Bf7n5rkiLGRSz7aWbVgTcJ9u/1ZGWLpxh/ZuPyO0iFEhIzyzaz5wDcfTVwH8EvomnAgMhjFcnAyCH2lwQ/jDfC1vsJvE5w/nkW8AXwhbuPCiVt2e10Pz34rKQv8LaZzSL4vOgfYQUug1j+LiuDWPbzfOBE4IrIsNvPzeywkPKWVSw/s3H5HaSpV0REJC50hCIiInGhQhERkbhQoYiISFyoUEREJC5UKCIiEhcqFBERiQsVikiCmdkvMWyTaWZTzCy9xGPdzWxq5NqH2WZ293ZeW93M3jMzzc0noVKhiKSGq4DhkYsiiaxHcQtwrrsfRjAnWNQLzdx9E/A2UOHW6pDKRRc2iiRY5AjlIGAc8AHBujC5QFd3z49s8xFwkbsvNrM6wLfAke6+KMb3OBT4q7ufnoh9EImFjlBEkmdfgkkiDwTygHNhy3xRe7n74sh23YCpsZZJxFcERzEioVGhiCTPt+7+eeT2DKBl5HYDgoLZ7CDgc3bAzP5rZn0334+cKttkZrXjF1ekdFQoIsmzscTtIn5b4C6fYCXLzX5lB/82I2tcjCZYHqGkGsCG8scUKRsVikjIPFhyNd3MNpfKOKC7mTUGMLMaZnZt5HZNoLu7vwzU3fw9zKw+sMrdC5KbXuQ3KhSR1DABOB7A3T8F7gFyIlOOfw40imzXD9jVzJ4FDjSzzMjjJwNjkppYZBsa5SWSAszscOAmd790B9s0B+5296sj9+8Gxrv71Mj69be6+4LkJBb5PRWKSIows6uAlzZfi1KK11UHerj7vxKTTCQ2KhQREYkLfYYiIiJxoUIREZG4UKGIiEhcqFBERCQuVCgiIhIXKhQREYkLFYqIiMSFCkVEROLi/wGCf+xogANSVwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import interpolate\n",
    "import statsmodels.api as sm\n",
    "\n",
    "# data will be a 2d list, which we convert to an array here\n",
    "data = np.array(data)\n",
    "t = data[:, 0]   # column 0\n",
    "Ca = data[:, 1]  # column 1\n",
    "\n",
    "# calculate numerical derivatives\n",
    "tck = interpolate.splrep(t, Ca)\n",
    "dCadt = interpolate.splev(t, tck, der=1)\n",
    "\n",
    "# do the transformation\n",
    "x = np.log(Ca)\n",
    "y = np.log(-dCadt)\n",
    "\n",
    "# setup and do the regression\n",
    "# column of ones and x:  y = b + mx\n",
    "X = np.column_stack([x, x**0])\n",
    "\n",
    "mod = sm.OLS(y, X)\n",
    "res = mod.fit()\n",
    "\n",
    "intercept = res.params[1]\n",
    "alpha = res.params[0]\n",
    "\n",
    "confidence_intervals = res.conf_int(0.05)\n",
    "intercept_error = confidence_intervals[1]\n",
    "alpha_error = confidence_intervals[0]\n",
    "\n",
    "\n",
    "print('alpha = {0}, conf interval {1}'.format(alpha, alpha_error))\n",
    "print('k = {0}, conf interval {1}'.format(np.exp(intercept), \n",
    "                                          np.exp(intercept_error)))\n",
    "\n",
    "\n",
    "# always visually inspect the fit\n",
    "plt.plot(x, y,'o')\n",
    "plt.plot(x, res.predict(X))\n",
    "plt.xlabel('$\\ln(C_A)$')\n",
    "plt.ylabel('$\\ln(-dC_A/dt)$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 495
    },
    "executionInfo": {
     "elapsed": 420,
     "status": "ok",
     "timestamp": 1588246512328,
     "user": {
      "displayName": "Zachary Ulissi",
      "photoUrl": "",
      "userId": "07633171379186475882"
     },
     "user_tz": 240
    },
    "id": "qZpIkBiThUO0",
    "outputId": "f286dfd3-b56a-4bc9-dcd1-3b6bd7ca43b2"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.6/dist-packages/statsmodels/stats/stattools.py:71: ValueWarning: omni_normtest is not valid with less than 8 observations; 7 samples were given.\n",
      "  \"samples were given.\" % int(n), ValueWarning)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.998</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.997</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   2210.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Thu, 30 Apr 2020</td> <th>  Prob (F-statistic):</th> <td>8.22e-08</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>11:35:11</td>     <th>  Log-Likelihood:    </th> <td>  13.785</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>     7</td>      <th>  AIC:               </th> <td>  -23.57</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>     5</td>      <th>  BIC:               </th> <td>  -23.68</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    2.0355</td> <td>    0.043</td> <td>   47.013</td> <td> 0.000</td> <td>    1.924</td> <td>    2.147</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>   -1.9646</td> <td>    0.157</td> <td>  -12.539</td> <td> 0.000</td> <td>   -2.367</td> <td>   -1.562</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>   nan</td> <th>  Durbin-Watson:     </th> <td>   2.377</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td>   nan</td> <th>  Jarque-Bera (JB):  </th> <td>   0.735</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.181</td> <th>  Prob(JB):          </th> <td>   0.692</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 1.454</td> <th>  Cond. No.          </th> <td>    40.4</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.998\n",
       "Model:                            OLS   Adj. R-squared:                  0.997\n",
       "Method:                 Least Squares   F-statistic:                     2210.\n",
       "Date:                Thu, 30 Apr 2020   Prob (F-statistic):           8.22e-08\n",
       "Time:                        11:35:11   Log-Likelihood:                 13.785\n",
       "No. Observations:                   7   AIC:                            -23.57\n",
       "Df Residuals:                       5   BIC:                            -23.68\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "x1             2.0355      0.043     47.013      0.000       1.924       2.147\n",
       "const         -1.9646      0.157    -12.539      0.000      -2.367      -1.562\n",
       "==============================================================================\n",
       "Omnibus:                          nan   Durbin-Watson:                   2.377\n",
       "Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.735\n",
       "Skew:                          -0.181   Prob(JB):                        0.692\n",
       "Kurtosis:                       1.454   Cond. No.                         40.4\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 20,
     "metadata": {
      "tags": []
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "J5P1bVrs2Tm9"
   },
   "source": [
    "\n",
    "-   You can see there is a reasonably large range of values for the rate constant and reaction order (although the confidence interval does not contain zero)\n",
    "\n",
    "-   The fit looks ok, but you can see the errors are not exactly random\n",
    "    -   There seems to be systematic trends in a sigmoidal shape of the data\n",
    "    -   That suggests small inadequacy in the model\n",
    "\n",
    "-   Let us examine some methods of evaluating the quality of fit\n",
    "\n",
    "-   First we examine the residuals, or the errors between the data and the model.\n",
    "\n",
    "-   In a good fit, these will be randomly distributed\n",
    "\n",
    "-   In a less good fit, there will be trends\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 283
    },
    "executionInfo": {
     "elapsed": 472,
     "status": "ok",
     "timestamp": 1588251490801,
     "user": {
      "displayName": "Zachary Ulissi",
      "photoUrl": "",
      "userId": "07633171379186475882"
     },
     "user_tz": 240
    },
    "id": "iqKO9K0I2Tm-",
    "outputId": "cf03be0c-dba2-48a3-e209-5e9a8530528a"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEKCAYAAAA8QgPpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU5bnA8d+TPSEbSQAhgIEk7C4oohCsu2Ctglav2k2trdfb2tZaabGuV7FSua21alttrcW2V229iMhiimtlUUBR2TNhEQiQlYTs63v/mDNhiBMymczMmZk8389nPpk5cybnORry5Jz3eZ9XjDEopZRSvRVldwBKKaXCkyYQpZRSPtEEopRSyieaQJRSSvlEE4hSSimfaAJRSinlkxi7AwimrKwsk5OTY3cYSikVNrKysigsLCw0xszq+l6/SiA5OTls3LjR7jCUUiqsiEiWp+16C0sppZRPNIEopZTyiSYQpZRSPtEEopRSyieaQJRSSvmkX1Vhhaslm0pYWLiTg9WNDEtPZO7MscyZnG13WEqpfk4TSIhbsqmEuxdvprG1HYCS6kbuXrwZQJOIUspWegsrxC0s3NmZPFwaW9tZWLjTpoiUUspJE0iIO1jd2KvtSikVLJpAQtyw9MRebVdKqWDRBBLi5s4cS0yUHLctMTaauTPH2hSRUko5aQIJcXMmZ3P6iDRcOWRAfDSPXn2KDqArpWynVVhhoLnNUJCXhTFQXtusyUMpFRL0CiTEdXQYisvqyB+cwvS8THaW1lJe22x3WEopZW8CEZFZIrJTRIpFZJ6H9+NF5GXr/Q9FJKfL+yNFpE5E7gpWzMFWUt1IY2s7Y4YkU5Dr7Ki8dleFzVEppZSNCUREooGngcuACcANIjKhy263AEeMMXnA48Avu7z/a2BloGO1k6OsFoD8IclMyk4jNSGGtcWVNkellFL2XoFMBYqNMbuNMS3AS8DsLvvMBhZZz18BLhIRARCROcAeYGuQ4rVFUWkdAHmDU4iOEs4ZnckavQJRSoUAOxNINrDf7fUBa5vHfYwxbUANkCkiycDPgP8OQpy2cpTWMSQ1nrTEWAAK8rI4cKSRfZUNNkemlOrvwnUQ/UHgcWNMXU87isitIrJRRDaWl5cHPjI/c5TVMmZISufrgrxMAFYX61WIUspediaQEmCE2+vh1jaP+4hIDJAGVAJnA4+JyF7gDuDnInK7p4MYY541xkwxxkwZNGiQf88gwDo6DI7SOvIGJ3duyx2UzJDUeL2NpZSynZ3zQDYA+SIyCmeiuB74Wpd9lgI3AuuAa4C3jTEGONe1g4g8CNQZY54KRtDBdKwC69gViIhQkJvFu0XldHQYorrMUldKqWCx7QrEGtO4HSgEtgP/MMZsFZGHRORKa7fncI55FAN3Al8o9Y1knRVYblcgANPzsqiqb2HH4Vo7wlJKKcDmmejGmBXAii7b7nd73gRc28P3eDAgwYUAh1WBlT845bjtrnGQtbsqmDAsNehxKaUUhO8ger9QVFrH4JR40pJij9s+NC2R0VkDWKMD6UopG2kCCWHFXSqw3E3Py+TDPVW0tncEOSqllHLSBBKiOjoMjrI68ocke3y/IDeLhpZ2PtlfHeTIlFLKSRNIiDpY00hDS/sXxj9cpuVmIoLexlJK2UYTSIhyDaCP6eYKJD0pjknD0rQvllLKNppAQlRRqauE1/MVCDjHQTbtP0JDS1uwwlJKqU6aQEKUo8xzBZa7gtwsWtsN6/dUBTEypZRy0gQSohyltd0OoLuclZNBXHQUa3fpbSylVPBpAglBxlgVWCe4fQWQGBfN5JHpOpCulLKFJpAQVFJtVWD1cAUCzvbuWw8epaq+JQiRKaXUMZpAQpCjzFWBdeIrEHAmEIB1ehtLKRVkmkBCkKPUcxNFT04bnkZyfIy2d1dKBZ0mkBDkKK1jUEo86UlxPe4bEx3F2aMyWKvjIEqpINMEEoKKyuq6nUDoyfS8LPZWNlBS3RjAqJRS6niaQEKMMYbi0toeK7Dcudq7azWWUiqYNIGEmIM1TdR7WYHlMnZIClnJcXobSykVVJpAQow3LUy6EhGm52axZlclzhV/lVIq8DSBhJjizlUIvb8CAedtrPLa5s4SYBV5lmwqoWDB24yat5yCBW+zZFOJ3SGpfs7WJW3VFxWV1pKVHM/AAT1XYLmbnuucD7KmuMKr+SMqvCzZVMLdizfT2NoOOCeb3r14MwBzJmfbGVrALdlUwsLCnRysbmRYeiJzZ46N+HMOF5pAQoyjlxVYLiMykhiZkcSa4kpuLhgVgMiUnRYW7uxMHi6Nre3c8+pmNu07QlSUEC1CdJQc99z1iBIhOgrr67FtMV32P/b82L5RUc79osV6v/P7ubZBtLVP52ek67Gt7xvt9hlru4h0e979OXGGA00gIcQYQ3FZHV89w7d/GAV5mSz79BBt7R3EROvdyUhysJsS7fqWdpZ8cpCODkO7MbR3WA9jCJfhMBGOJacuia+6oYWOLufR2NrOwsKdmkBCgCaQEHKopom65jbyfbwFNT03ixfX72dzSQ2TRw70c3TKLss+O9jte9npiayZd6HH94w5lkw6OuhMMB0dhrYOQ4dbwnE9d36l83lb1/fdEtWxfTucX431fucxj309/vvg8djtHj77tw/2eTy37hKqCi5NICHEVYHl6xjG9FznfJC1uyo1gUSA1vYOFqzcwXOr93ByRiKlR5tpauvofD8xNpq5M8d2+3kRISZawvof+Ts7yj1OkB2WnmhDNKorvc8RQorLfKvAcslMjmf80FRWO3Q+SLgrq23i63/6kOdW7+Gm6TmsuvN8Fnz1VLLTExGcVx6PXn1KxN/GmTtzLImx0cdti4mSEyZOFTzh/MdJxPG1AstdQW4mL3zwOU2t7SR0+YenwsPGvVV87+8fc7SplcevO42rJg8HnIPGkZ4wunKdr6sKKz42ivZ2wzTralvZS69AQkhRaZ3PVx8uBXlZtLR1sHHvET9FpYLFGMNf1uzh+mc/IDEumle/V9CZPPqzOZOzWTPvQvYsuJx/3XEeIsLCwp12h6XQBBIyXBVYvpTwups6KoOYKNH27mGmoaWNO17+hAdf38b5Ywex9PYZjB+aandYIWdkZhI3F+Twfx8fYEtJjd3h9HuaQEKEqwIrr4+TAAfEx3D6iHTtixVG9lbUc/Xv1rL004PcdekYnv3mFNISY+0OK2R9/8I8MpLieGjZNm3dYzNNICGicxXCPt7CAmd7980lNdQ0tvb5e6nAWrWtlCueWs3ho00sunkqt1+YT1RU9xPrFKQmxPLjS8awfk8VhVsP2x1Ov6YJJER0rkLohzYkBbmZdBj4YLcucxuq2jsM/1O4k+++sJGczAG8fvsMvjRmkN1hhY3rzxrBmCHJPLpyB81t7T1/QAWEJpAQ4SitIys5jow+VGC5TB45kMTYaL2NFaKq6lu46fn1PPVOMddNGcE/b5vGiIwku8MKKzHRUdxz+QQ+r2zghbWf2x1Ov6UJJEQUldWS54fbVwBxMVFMHZXBak0gIeezA9Vc8eRqPtxdxYKrT+GX15yq5dY+Om/MIM4fO4jfvu2gsq7Z7nD6JU0gIcC5CmGdX7voFuRlsqu8nsM1TX77nqpvXlq/j2t+vw6Af942jeunjrQ5ovB37+XjaWhp5zdvOuwOpV/SBBICDh9torYPPbA8cbV3X6vlvLZram3nZ698xrzFmzl7dAav/2AGp41ItzusiJA3OIWvnz2S/12/r3McUQWPrQlERGaJyE4RKRaReR7ejxeRl633PxSRHGv7JSLykYhstr567iYXJop8XETqRCYMTWVgUixrinUg3U4HjjRw7R/W8fLG/dx+QR5/uXmqX8a51DF3XDyGpLhoHlmx3e5Q+h3bEoiIRANPA5cBE4AbRGRCl91uAY4YY/KAx4FfWtsrgCuMMacANwJ/DU7UgeHoYxNFT6KihGm5mazdVaG18jb5d1E5X3lyNXsr6vnjt6Zw18yxRGuJrt9lDIjjhxfm8+7Oct4rKrc7nH7FziuQqUCxMWa3MaYFeAmY3WWf2cAi6/krwEUiIsaYTcYYV4/rrUCiiMQHJeoAcJTWkTnAPxVY7qbnZnGopok9FfV+/b7qxDo6DE+97eDG59czJCWBpT+YwSUThtgdVkT71vSTOTkzifnLttHW3tHzB5Rf2JlAsoH9bq8PWNs87mOMaQNqgK5d1L4KfGyM8ViGISK3ishGEdlYXh6af504ymrJ72MLE09m5FnL3O7S21jBUtPYyq1/3cj//KuIK08bxqvfn86orAF2hxXx4mOiufuy8TjK6nhxw/6eP6D8IqwH0UVkIs7bWv/Z3T7GmGeNMVOMMVMGDQq9iVrGGByldeQP9v865idnJjkXHNL27kGx4/BRZj+1mnd3lvPgFRP4zXWnkxSnDa+DZebEIZw9KoPHVxVpF4YgsTOBlAAj3F4Pt7Z53EdEYoA0oNJ6PRx4FfiWMWZXwKMNkNKjzdQ2t/W5iaInIsL03EzW7a6kveu6oMqvlmwqYc7Ta2hoaeelW8/hpoJRJ1zrW/mfiHDfVyZwpKGF371TbHc4/YKdCWQDkC8io0QkDrgeWNpln6U4B8kBrgHeNsYYEUkHlgPzjDFrghZxALhWIcwLwBUIONu71zS2su3g0YB8//6upa2DB5du5Y6XP+HU7HSW/XAGU3Iy7A6r35qUncZXzxjO82v2sq+ywe5wIp5tCcQa07gdKAS2A/8wxmwVkYdE5Eprt+eATBEpBu4EXKW+twN5wP0i8on1GBzkU/CLziaKAbgCgWPL3Gp7d/8rPdrEDX/8gL+s3cstM0bx9++ezeCUBLvD6vfmzhxLTLTw6Eot6w00W2/QGmNWACu6bLvf7XkTcK2Hz80H5gc8wCBwlNaSMSCOzOTAFJENTk0gf3Aya4oruO283IAcoz/6YHclt//vJhpa2njyhslccdowu0NSliGpCdx2Xi6/XlXEh7srOXu0rl4YKGE9iB4Jikpr/TqB0JOCvCw27K3SrqV+YIzhT+/v5ut/+pDUhBiWfL9Ak0cI+u65oxmalsD85dvp0PG/gNEEYiNjDI4y//bA8qQgL4um1g427asO6HEiXX1zG7e/uIn5y7dz8fjBvHZ7QcD/3ynfJMZF89NZY9lcUsOrm7rW5ih/0QRio9KjzdQ2tQVkDoi7s0dnECWwRrvz+mxXeR2zn17Dys2HmHfZOP7wjTNJSdBVA0PZ7NOyOW14Go8V7qChpc3ucCKSJhAbOcqsRaQCVIHlkpoQy6nD0zWB+OiNLYeY/dQaqupb+NstZ3PbeblaohsGoqKcZb2lR5t55r3ddocTkTSB2KiziWKAr0DA2d790wM11DbpBCtvtbV38OjK7dz2t4/JHZzMsh/MYLo1u1+Fhyk5GVx+6lCe+fcuDtU02h1OxNEEYqPiMmcFVlaAKrDcFeRm0d5hWL+nKuDHigQVdc1887n1PPPebr5+9kj+8Z/nMCw90e6wlA/mzRpHh4GFb+y0O5SIownERkWldX5bhbAnZ5w8kPiYKG3v7oVN+45wxZOr+XjfERZecyqPXHUK8TG6amC4GpGRxLcLRrF4UwmfHdBCEn/SBGITZw+s2oBNIOwqITaaKTkDdYGpEzDG8LcPPuc/nllHTLTwf/81nWunjOj5gyrkff+CXLKS43h42TZd3sCPNIHYpKy2maNNbQEfQHdXkJfFjsO1lNfq+tFdNbW2c9c/P+PeJVsoyMvi9dtnMCk7ze6wlJ+kJMRy5yVj2bD3CCu3HLY7nIihCcQmrh5YwRhAdynQZW492lfZwNW/W8viTQf40UX5/PnGs0hP0lUDI811Z41g3EkpPLpyO02tOqnWHzSB2MRR6uqBFbwrkEnZaaQmxLBWx0E6vbOjjCueWs2BIw38+caz+PElY4jSVQMjUnSUcO/lE9hf1chf1u61O5yIoAnEJo6yWgYmxZIZxPWxo6OEc0ZnamNFnKsGPr6qiG8v2sCw9ESW/eBcLhgXlv04VS/MyM/iwnGDefrtYirq9FZuX2kCsYmjtI78ISlBn5BWkJfFgSON/brVdXVDC7cs2sATbzm4anI2i/9rOiMzk+wOSwXJz788nsbWdh5fVWR3KGFPE4gNjDFBaaLoSUFe/27vvqWkhiueWs3q4grmz5nEr649jcQ4LdHtT/IGJ/ONc07mxfX72Hm41u5wwpqut2mDcqsCy45GfLmDkhmcEs+a4gpumDoy6McPpiWbSlhYuJOD1Y0MS0/kS2OyWPxxCQOT4vjHf05j8siBdoeobPKji/J5dVMJ85dv44VvT9XWND7SKxAbdLYwseEKRESYkZfFul2VEd3mesmmEu5evJmS6kYMUFLdyIvr9zNiYCLLfjhDk0c/N3BAHD+8KJ/3HRW8u7Pc7nDCliYQG3Q2UbSpFfj0vCwq61vYEcGX7wsLd9LooVSzobU9KK1jVOj75jknMyprAPOXb6O1vcPucMKSJhAbFJXWkZ4US1ayPXMNXOMgkTwf5GC158Z5h6qbghyJClVxMVHcfdk4dpXX8+L6fXaHE5Y0gdiguKyWMYODX4HlMjQtkdFZAyK6vXt3jQ+1IaJyd8mEIUwbncnjq4qoadBO1b2lCSTInBVYdeQFcQa6J9PzMlm/pypiL91/cOEX139PjI1m7syxNkSjQpWIcO9XxlPd2MqTbzvsDifseJVARKRARAZYz78hIr8WkZMDG1pkKq9tpqaxlTE2DKC7K8jNor6lnU/3R2Z30h2HnYUKg1LiESA7PZFHrz6FOZOz7Q1MhZyJw9L4jzNHsGjdXvZU1NsdTljx9grk90CDiJwG/ATYBbwQsKgimKMs+C1MPJmWm4kIEdnefUtJDS+s28u3pp3MhnsuZs+Cy1kz70JNHqpbP7l0DLHRUSxYud3uUMKKtwmkzTh7IM8GnjLGPA3Y+xswTLmaKNp9Cys9KY5Jw9IibkJhe4fhnlc3kzEgnp9cqrerlHcGpybwvfNzKdxayrpdkfdHVaB4m0BqReRu4BvAchGJAmIDF1bkcpQ5K7AGhUAp6fS8TDbtO0JDS5vdofjN/67fx6cHarjvK+NJS9QfUeW975w7muz0ROYv30Z7BM+R8idvE8h1QDNwizHmMDAcWBiwqCKYw2phEgozXwtys2htj5xlbstqm3jsjR0U5GVy5WnD7A5HhZmE2Gh+OmssWw8e5f8+PmB3OGHBqwRijDlsjPm1MeZ96/U+Y4yOgfSSqwLLrgmEXZ2Vk0FcdBRrI+SS/RfLt9Pc2sHDsyeFRIJW4efK04YxeWQ6Cwt3Ut8cOVfmgXLCBCIitSJy1MOjVkSOBivISFFe56zAsqOFiSeJcdFMHpkeEfNB1hZXsOSTg9x23mhGDwqN/74q/Ig41wwpr23mmfd22R1OyDthAjHGpBhjUj08UowxqcEKMlIU27CIVE8K8rLYdugoR+pb7A7FZ81t7dy7ZAsjM5L43gV5doejwtyZJw/kitOG8ez7u7vtaKCcejWRUEQGi8hI1yNQQUWqzmVsQ+QKBJxtTYyBdbvD9zbWs+/tZndFPQ/NnkhCrLZmV333s1ljMQYee2OH3aGENG8nEl4pIg5gD/AesBdYGcC4IlJRWR1pibEMSrG/Asvl1OHpJMfHhO1trH2VDTz1TjGXnzKU88fqioLKP4YPTOI7545iyScH+SRCJ9v6g7dXIA8D5wBFxphRwEXABwGLKkIVl9aFTAWWS2x0FGePygjLgXRjDPe9toWYKOG+r0ywOxwVYf7r/DyykuN5eNk2nNPgVFfeJpBWY0wlECUiUcaYd4ApAYwr4hhjKCqrDZkKLHfT87LYU1FPSZjd71255TDvFZXzk0vHclJagt3hqAiTHB/DXZeO4aPPj7B88yG7wwlJ3iaQahFJBv4N/F1EngC0aUwvVNS1UN3QyhibZ6B70rnMbRjdxqprbuOh17cxYWgq35qmbdlUYFw7ZQTjh6ayYOUOmjysL9PfeZtAZgONwI+BN3D2wrqirwcXkVkislNEikVknof340XkZev9D0Ukx+29u63tO0VkZl9jCTRH5wB66F2BjB2SQlZyHGvDKIE8vqqI0tomHrlqEjHR2lRaBUZ0lHDf5eM5cKSRP6/ZY3c4IcfbiYT1xph2Y0ybMWaRMea31i0tn4lINPA0cBkwAbhBRLreyL4FOGKMyQMeB35pfXYCcD0wEZgF/M76fiHrWBPF0LsCERGm5WaxZldlWNzr3XqwhufX7OFrU0fq0rQq4KbnZXHx+CH87p1dlNc22x1OSPG2Cst9QmGTiLT7YSLhVKDYGLPbGNMCvITzSsfdbGCR9fwV4CJxjkDPBl4yxjQbY/YAxdb3C1lFpbWkJsSEVAWWuxl5mZTXNlNsJbpQ1dFhuHfJFjIGxPHTmePsDkf1Ez//8jiaWtv59aqddocSUry9AumcUAgkAl8FftfHY2cD+91eH7C2edzHGNMG1ACZXn4WABG5VUQ2isjG8vLyPobsO0dZHWOG2LcKYU+m52YBoT8O8tKG/WzaV83PvzyetCRtlqiCY/SgZL457WRe3rCf7Ye0CYdLr28eG6clQMiPOwAYY541xkwxxkwZNGiQXTE4myiG4O0rlxEZSYzMSGJNCJfzVtQ1s2Dlds4ZncFVuraHCrIfXZRPSkIsjyzfHha3eoPB21tYV7s9rhGRBUBTH49dAoxwez3c2uZxHxGJAdKASi8/GzIq6lo40tAakgPo7gryMvlgVyVtIbrM7S9WbKextZ35c7RZogq+9KQ47rg4n9XFFby9o8zucEKCt1cgV7g9ZgK1fHG8orc2APkiMkpE4nAOii/tss9S4Ebr+TXA29bCVkuB660qrVFAPrC+j/EEjKPMqsAK4SsQcN7Gqm1uY3NJjd2hfMG6XZUs/riEW780mrwQT8Qqcn3jnJMZPWgAj6zYTmuI/qEVTDHe7GSMudnfBzbGtInI7UAhEA382RizVUQeAjYaY5YCzwF/FZFioApnksHa7x/ANqAN+L4xJmSLtB0h2ETRk+m5zvkga3dVhlR1U0tbB/e9toXhAxO5/YJ8u8NR/VhsdBT3fHk8tyzayN8++JybC0bZHZKtTphARORJoNubfcaYH/bl4MaYFcCKLtvud3veBFzbzWcfAR7py/GDxVHmrMAaHKIVWC6ZyfGMOymFNcUVfD+Eutr+8f3dFJfV8fxNZ5EYF9LV2qofuHDcYAryMvnNmw6umpxNelKc3SHZpqdbWBuBj4AE4AzAYT1OB/rvf7Veci0iFQ737WfkZbHx8yMhM+t2f1UDT77tYNbEk7hgnDZLVPZzrRlS29TKb98qtjscW/W0HsgiY8wi4FTgfGPMk8aYJ3E2Uzw9GAFGguKyupCcQOhJQV4WLW0dfPT5EbtDwRjDA0u3EiXC/Vdos0QVOsYPTeW6s0bwwrq97C4P7blTgeTtIPpAwH0BqWRrm+pBRV0zVfUtYTPwO3VUBjFREhLzQQq3lvL2jjLuvGQMw9IT7Q5HqePceclYEmKj+cWK/rtmiLcJZAGwSUT+IiKLgI+BXwQurMhxbAA9PK5ABsTHcPoI+5e5rW9u479f38q4k1K4aXqOrbEo5cmglHi+d0Eub24vDas+cv7k7Uz054GzgVeBxcA069aW6kFnCW+YXIGAs/fP5pIaahpbbYvhibccHKrRZokqtH27YBTZ6Yk8vHw77R39b3LhCf9lisg46+sZwDCc7UP2A8OsbaoHRaW1pCTEMCQ1tCuw3BXkZtJh4AOblrndfugoz63ew/VnjeDMkzNsiUEpbyTERjPvsnFsP3SUVz7a3/MHIkxP80DuBG4FfuXhPQNc6PeIIowjBFch7MnkkQNJjI1mbXEFMyeeFNRju5olpiXG8rNZ2ixRhb6vnDqUv6zdy8LCIi4/dRjJ8V5Nr4sIPVVh3Wp9vcDDQ5OHF1xNFMNJXEwUU0dl2NIX658f7eejz49w92XjGDhAK8VV6BNxLqlcUdfM79/tX2W93vbCulZEUqzn94rIYhGZHNjQwl9lZwVWeAyguyvIy6S4rI7So31teea9qvoWHl25g6k5GVxz5vCgHVepvjp9RDpzTh/GH9/fw4EjDXaHEzTejk7eZ4ypFZEZwMU4W4z8IXBhRYaiMGlh4omrvfvaXcGrLnl0xXbqmtqYf5U2S1ThZ+6scQjw2Bv9Z80QbxOIa1ry5cCzxpjl6Ez0HhVbFVjhmEAmDE1lYFIsqx3BuY21fk8V//zoAN85d3RY/vdSKjs9kVu/NJqlnx7k4332T8QNBm8TSImIPANcB6wQkfhefLbfKiqtIyU+vCqwXKKihGm5mazdVRHwtQ9a2zu4d8lmstMT+eFFodODS6neuu28XAanxPPwsm39Ys0Qb5PAf+DsmjvTGFMNZABzAxZVhHCUOReRCtfbMdNzszhU08SeivqAHue51XsoKq3jv6+cSFJc/6lgUZFnQHwMd80cy6Z91Sz99KDd4QSctxMJG4AyYIa1qQ1nU0V1As4S3vC9HVOQZy1zG8BqrANHGnjiTQeXTBjCxROGBOw4SgXLNWcMZ+KwVH65ckfINCUNFG+rsB4AfgbcbW2KBf4WqKAiQWVdM5X1LSG/iNSJ5GQmkZ2eGNA2Df/9+jYAHrxyYsCOoVQwRUU5u/UerGniudV77A4noLy9hXUVcCVQD2CMOQiE75/WQeAoc1Zg5YfxgLCIMD03k3W7K+kIQJuGVdtKWbWtlB9dnE+2NktUEWRabiaXThjC794ppqw2eKXwweZtAmmxlpI1ACIyIHAhRQZHqasCK3yvQMB5G6u6oZVth4769fs2tLTx4NKtjBmSzC0z+veqbioy/fzL42lp7+BXhUV2hxIwPSYQcY4AL7OqsNJF5LvAm8AfAx1cOHOUOSuwTkpNsDuUPnEtc7vaz7exfvtWMSXVjTxy1SnEarNEFYFysgZw47Qc/vHRfrYerLE7nIDo8V+udeVxLfAK8H/AWOB+a2Ep1Y2i0lrywrgCy2VwagL5g5P92t595+Fa/vT+bq49czhn5WizRBW5fnBRPumJscxftj0iy3q9/dPvY6DaGDPXGHOXMWZVIIOKBMVldYwJ49G6JqwAABQfSURBVAosdwV5WWzYW0VzW98rSowx3LdkC8kJMdz95fF+iE6p0JWWGMsdF49h3e5K3txeZnc4fudtAjkbWCciu0TkM9cjkIGFs6r6FirqwrsCy9303EyaWjvYtK+6z9/rlY8OsH5vFXdfNo4MbZao+oGvnT2S3EED+MWK7bS0ddgdjl95m0BmArk427df4fZQHrgG0MO5AsvdObmZRAl9Luc9Ut/CL1Zs58yTB3LtmSP8FJ1SoS02Oop7L5/Anop6/vrB53aH41feTiT83NMj0MGFqyJXCW8YduH1JDUhllOHp/d5QuEv39jB0aY2HrlqElFR4T02pFRvnD92EOfmZ/HEm0UcqW+xOxy/0fKXACgurSU5PoahaeFdgeWuIC+TT/ZXU9vk2zK3H31exUsb9nPLjFGMOynVz9EpFdpEnJML65rbeOKtyGnioQkkAIpK68gLs1UIe1KQm0V7h2H9nqpef7a1vYN7Xt3CsLQEfnRRfgCiUyr0jT0pheunjuSvH3xOsXWXItxpAgkAR1lt2E8g7OqMkwcSHxPFmuLe38b6y5q97Dhcy/1XTGRAP1ruU6mu7rxkDEmx0Ty6YrvdofiFJhA/66zAipASXpeE2Gim5Azs9QJTB6sbefzNIi4aN5iZE7VZourfspLj+f6Feby1o4zVjuAt1hYomkD87FgFVmRdgYCzvfuOw7VU1DV7/ZmHXt9GhzE8eOXEiLqlp5Svbi7IYURGIvOXb6M9AD3mgkkTiJ9FQhPF7szIcy1z691trLd3lPLG1sP84MJ8RmQkBTI0pcJGfEw0d182nh2Ha3l5w367w+kTTSB+5rAqsIZFUAWWy6TsNFITYryaD9LY0s79r20lb3Ay3z13dBCiUyp8XDbpJM7KGcivV+30ubIxFGgC8TNHWeRVYLlERwnnjM5kjRfjIE+94+DAkUbmz5lEXIz+mCnlzlXWW1HXwu/e3WV3OD7Tf9l+VlRaFzETCD0pyMtif1Uj+yobut2nuKyWZ/+9m6vPyOac0ZlBjE6p8HHaiHSunpzNc6v3sL+q+39PoUwTiB8dqW+hoq6ZMRE4/uFSkOdMCN1dhRhjuHfJFpLiYvi5NktU6oTmzhpLlMCCN3bYHYpPbEkgIpIhIqtExGF9HdjNfjda+zhE5EZrW5KILBeRHSKyVUQWBDf67rkG0PMisALLJXdQMoNT4rtt7/7qphI+2F3Fz2aNIys5PsjRKRVehqYlcuuXcln+2SE27u39JF272XUFMg94yxiTD7xlvT6OiGQAD+DsBDwVeMAt0fyPMWYcMBkoEJHLghP2iRV1rkIYuVcgIkJBXhbrdn1xmdvqhhYeWb6dySPTuf4sbZaolDduO280Q1LjeXj59oAsHR1IdiWQ2cAi6/kiYI6HfWYCq4wxVcaYI8AqYJYxpsEY8w6AMaYF51olw4MQc4+Ky+oYEBcdkRVY7grysqisb2GnlTBdHivcyZGGFubP0WaJSnkrKS6GuTPH8en+apZ+etDucHrFrgQyxBhzyHp+GPA0RTkbcC+SPmBt6yQi6Tjbyr8ViCB7y7kKYUpEVmC56xwHcbuN9fG+I7y4fh83F4xi4rA0u0JTKixdPTmbU7LT+OUbO2hs6fvCbcESsAQiIm+KyBYPj9nu+1lL5vb6uk1EYoAXgd8aY3afYL9bRWSjiGwsLy/v9Xn0hqMssiuwXIamJTI6a0DnhMI2q1nikJQEfnzJGJujUyr8REUJ914+nkM1Tfzx/W5/nYWcgCUQY8zFxphJHh6vAaUiMhTA+upprccSwP1G+nBrm8uzgMMY85se4njWGDPFGDNl0KBBfTupE6huaKG8tjnimih2Z3peJh/urqS1vYNF6z5n+6Gj3H/FBJK1WaJSPjl7dCazJp7E79/dRenRJrvD8Ypdt7CWAjdaz28EXvOwTyFwqYgMtAbPL7W2ISLzgTTgjiDE6pXOFiYR1kSxO7FRQn1LO/n3rGT+sm2MPymFyyadZHdYSoW1u788jvYOw/8U7rQ7FK/YlUAWAJeIiAO42HqNiEwRkT8BGGOqgIeBDdbjIWNMlYgMB+4BJgAfi8gnIvIdO07CXVEEN1HsasmmEl506+FjgN0V9bz2SXgNACoVak7OHMBNBTm88vEBtpTU2B1Oj2xJIMaYSmPMRcaYfOtWV5W1faMx5jtu+/3ZGJNnPZ63th0wxogxZrwx5nTr8Sc7zsOdo9RZgZWdnmh3KAG3sHAnTa0dx21rbutgYZj81aRUKLv9wjwGJsXx8LJtOIeIQ5fORPcTR1ltxPbA6upgdWOvtiulvJeaEMuPLxnDh3uqKNxaanc4J6QJxE8cpXUR2cLdk2HdXGV1t10p1Ts3nDWC/MHJPLpyO81toVvWqwnED2oaWimrbe4XJbwAc2eOJTE2+rhtibHRzJ051qaIlIosMdFR3HP5eD6vbOCv6z63O5xuaQLxg6KyyG9h4m7O5GwevfoUstMTESA7PZFHrz6FOZOze/ysUso7548dzHljBvHEWw6q6lvsDscjLdr3A0ep1USxn1yBgDOJaMJQKrDuvXw8s554n9+8WcRDsyfZHc4X6BWIHxSV1pLUTyqwlFLBkz8kha9NHcnfP9yHo0vvuVCgCcQPiq1VCLWBoFLK3+64OJ+kuGgeWbHd7lC+QBOIHxSV1vabGehKqeDKTI7nBxfm8e7Oct4rCmw/v97SBNJHrgqs/tIDSykVfDdOz+HkzCQeWb6NtvaOnj8QJJpA+shR1n9amCil7BEfE83dl42jqLSOlzbs7/kDQaIJpI/6WxNFpZQ9Zk48iamjMnh8VRFHm1rtDgfQBNJnRaW1JMZqBZZSKrBEhPsun0BVQwtPv11sdziAJpA+c7Yw0QospVTgnTI8jasnD+f5NXvZV9lgdziaQPrK1URRKaWC4aezxhIdJSx4w/6yXk0gfVDT2Erp0eZ+08JEKWW/IakJ3HZeLis2H2b9nipbY9EE0gfFrgosvQJRSgXRrV8azUmpCTy8bBsdHfatGaIJpA+KrB5YegWilAqmxLhofjprLJtLanh1U4ltcWgC6QNHaZ1WYCmlbDHn9GxOHZ7GY4U7aGhpsyUGTSB94BpA1wospVSwRUUJ931lAqVHm3n237vticGWo0YIVwmvUkrZ4aycDC4/ZSjPvLebwzVNQT++JhAf1TS2cvhok85AV0rZat5l42jvMDxWuCPox9YE4qPizlUI9QpEKWWfERlJ3Dwjh8Ufl/DZgeqgHlsTiI9cqxDqFYhSym63X5BH5oA4Hl62DWOCV9arCcRHRaV1JMRGMXygVmAppeyVkhDLnZeOYcPeI7yx5XDQjqsJxEdagaWUCiXXTRnB2CEpPLpyB81t7UE5piYQHzlK6xijt6+UUiEiJjqKe78ynn1VDfxlzd6gHFMTiA+ONjkrsPJ0AF0pFULOzR/EBWMH8dTbxVTUNQf8eJpAfOAaQNcrEKVUqLnn8vE0tLbz+KqigB9LE4gPinUZW6VUiMobnMI3zh7Ji+v3UVRaG9BjaQLxgasCa8TAJLtDUUqpL7jj4jEkx8cwf3lg1wzRBOKDolKtwFJKha6BA+L44UX5/LuonHd2lgXsOJpAfFBcVqcTCJVSIe1b03LIGhDLdxdtZNS85RQseJslfm79HuPX79YPHG1q5VBNk45/KKVC2orNhzja1EabteBUSXUjdy/eDMCcydl+OYZegfRScZm2MFFKhb6FhTtpaT++rUljazsLC3f67Ri2JBARyRCRVSLisL4O7Ga/G619HCJyo4f3l4rIlsBHfIyjVJsoKqVC38Hqxl5t94VdVyDzgLeMMfnAW9br44hIBvAAcDYwFXjAPdGIyNVAXXDCPcZRWkd8TBTDtQJLKRXChnWzUmp3231hVwKZDSyyni8C5njYZyawyhhTZYw5AqwCZgGISDJwJzA/CLEep6isjrzByURrBZZSKoTNnTmWxNjo47YlxkYzd+ZYvx3DrgQyxBhzyHp+GBjiYZ9sYL/b6wPWNoCHgV8BDT0dSERuFZGNIrKxvLy8DyE7FZfWkj9Yb18ppULbnMnZPHr1KWSnJyJAdnoij159it8G0CGAVVgi8iZwkoe37nF/YYwxIuJ1A3sROR3INcb8WERyetrfGPMs8CzAlClT+tQov7aplYM1TeQP0QF0pVTomzM5268Jo6uAJRBjzMXdvScipSIy1BhzSESGAp5mupQA57u9Hg68C0wDpojIXpzxDxaRd40x5xNgrgqsMZpAlFLKtltYSwFXVdWNwGse9ikELhWRgdbg+aVAoTHm98aYYcaYHGAGUBSM5AHuqxDqLSyllLIrgSwALhERB3Cx9RoRmSIifwIwxlThHOvYYD0esrbZpqi0lviYKEZkaAWWUkrZMhPdGFMJXORh+0bgO26v/wz8+QTfZy8wKQAheuQoqyN3kFZgKaUU6Ez0XnGU1uoEQqWUsmgC8ZJWYCml1PE0gXjpWA8svQJRSinQBOI1hyuB6BWIUkoBmkC85iitJS4mipFagaWUUoAmEK9pBZZSSh1PE4iXHKV1WoGllFJuNIF4oa65jZLqRm1hopRSbjSBeMFVgZWnFVhKKdVJE0gPlmwq4abn1wNw/2tb/L4ovVJKhStbWpmEiyWbSrh78WYaW9sBKD3a7PdF6ZVSKlzpFcgJLCzc2Zk8XPy9KL1SSoUrTSAnEIxF6ZVSKlxpAjmBYCxKr5RS4UoTyAkEY1F6pZQKVzqIfgKugfKFhTs5WN3IsPRE5s4cqwPoSimFJpAeBXpReqWUCld6C0sppZRPNIEopZTyiSYQpZRSPtEEopRSyieaQJRSSvlEjDF2xxA0IlIOfG7DobOAChuOG0z94RxBzzOS9IdzhL6fZwWAMWZW1zf6VQKxi4hsNMZMsTuOQOoP5wh6npGkP5wjBPY89RaWUkopn2gCUUop5RNNIMHxrN0BBEF/OEfQ84wk/eEcIYDnqWMgSimlfKJXIEoppXyiCUQppZRPNIEEkIj8RESMiGR18/6NIuKwHjcGO76+EJGHReQzEflERP4lIsO62e8xEdkqIttF5LciIsGOtS96cZ4jrfe3i8g2EckJbqR94+15WvumisgBEXkqmDH2lTfnKCKni8g662f2MxG5zo5Y+6IXP7N9//1jjNFHAB7ACKAQ58TFLA/vZwC7ra8DrecD7Y67F+eX6vb8h8AfPOwzHVgDRFuPdcD5dsfu7/O03nsXuMR6ngwk2R17IM7Tev8J4H+Bp+yO29/nCIwB8q3nw4BDQLrdsQfgPP3y+0evQALnceCnQHdVCjOBVcaYKmPMEWAV8IWZnqHKGHPU7eUAPJ+nARKAOCAeiAVKAx+d/3hzniIyAYgxxqyyPlNnjGkIUoh+4eX/T0TkTGAI8K9gxOVP3pyjMabIGOOwnh8EyoBBwYnQP7z8f+mX3z+6oFQAiMhsoMQY8+kJ7thkA/vdXh+wtoUNEXkE+BZQA1zQ9X1jzDoReQfnX3GC8y/W7cGNsu96Ok+cf7VWi8hiYBTwJjDPGNMevCj7rqfzFJEo4FfAN4CLgxudf3jx/9J936k4//jZFYTQ/MqL8/TL7x+9AvGRiLwpIls8PGYDPwfutzvGvurhHDHG3GOMGQH8Hbjdw+fzgPHAcJw/nBeKyLnBPAdv9PU8cf4hdi5wF3AWMBq4KUjhe80P5/k9YIUx5kAw4+4NP5yj6/sMBf4K3GyM6QhO9N7z13n2md336yLtAZyC87J3r/VoA/YBJ3XZ7wbgGbfXzwA32B2/j+c8EtjiYftc4D631/cDP7U73gCc5znAe26vvwk8bXe8ATjPv1s/y3txNtg7CiywO15/nqP1XirwMXCN3XEG8P+lX37/6BWInxljNhtjBhtjcowxOTgvDc8wxhzusmshcKmIDBSRgcCl1rawICL5bi9nAzs87LYPOE9EYkQkFjgPCKtbWF6e5wYgXURc98ovBLYFOjZ/8uY8jTFfN8aMtH6u7wJeMMbMC1KIfebNOYpIHPAqznN7JVix+ZOXP7N++f2jCSSIRGSKiPwJwBhTBTyM85fPBuAha1u4WGBdMn+G84fvR3D8OQKv4Lx/vBn4FPjUGPO6LdH6rsfzNM6xjruAt0RkM87xnj/aFbCPvPn/Ge68Ocf/AL4E3GSVwX4iIqfbFK+vvPmZ9cvvH21lopRSyid6BaKUUsonmkCUUkr5RBOIUkopn2gCUUop5RNNIEoppXyiCUQppZRPNIEo5WciUufFPoki8p6IRLttu1ZEPrTmHmwVkQe6+WyciPxbRLSXnbKVJhCl7PFtYLE1CRFrPYafAV81xpyOs6eWx4ldxpgW4C0g7NaqUJFFJxIq5WfWFcgkYCWwGue6KCXAbGNMo7XPWuBrxpi9IpIK7AHOMsbs9vIYpwGPGmO+HIhzUMobegWiVODk42yqOBGoBr4Knf2WRhtj9lr7zQE+9DZ5WLbgvEpRyjaaQJQKnD3GmE+s5x8BOdbzLJwJxWUS8AknICIvi8hdrtfWra8WEUnxX7hK9Y4mEKUCp9nteTvHFnBrxLlSo0s9J/i3aK3xsAznUgHu4oGmvoeplG80gSgVZMa5hGi0iLiSyErgWhEZAiAi8SLyXet5AnCtMeavQJrre4hIJlBhjGkNbvRKHaMJRCl7/AuYAWCMWQ88CBRaLbg/AQZb+80FkkXkD8BEEUm0tl8ALA9qxEp1oVVYStlARM4AfmyM+eYJ9hkJPGCMucV6/QDwhjHmQ2v99XnGmKLgRKzUF2kCUcomIvJtYJFrLkgvPhcHXG+MeSEwkSnlHU0gSimlfKJjIEoppXyiCUQppZRPNIEopZTyiSYQpZRSPtEEopRSyieaQJRSSvlEE4hSSimfaAJRSinlk/8H94gZeE7FkK4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light",
      "tags": []
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "residuals = y - res.predict(X)\n",
    "\n",
    "# always visually inspect the fit\n",
    "plt.plot(x, residuals, 'o-')\n",
    "plt.xlabel('$\\ln(C_A)$')\n",
    "plt.ylabel('residuals')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "renlH8Vs2TnB"
   },
   "source": [
    "\n",
    "-   You can see there are trends in this data\n",
    "    -   That means the model may not be complete\n",
    "\n",
    "-   There is uncertainty in the data\n",
    "    -   In each concentration measurement there is uncertainty in the time and value of concentration\n",
    "    -   You need more data to reduce the uncertainty\n",
    "    -   You may also need better data to reduce the uncertainty\n",
    "\n",
    "-   Derivatives tend to *magnify* errors in data\n",
    "    -   The method we used to fit the data contributed to the uncertainty\n",
    "\n",
    "-   We also *nonlinearly* transformed the errors by taking logs and exp of the data and results, which may have skewed the confidence limits\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "n-69aqdp2TnJ"
   },
   "source": [
    "\n",
    "## Nonlinear regression\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wvJkiE7x2TnK"
   },
   "source": [
    "-   Nonlinear models are abundant in reaction engineering\n",
    "    -   $r = k C_A^n $ is linear in the $k$ parameter, and nonlinear in $n$\n",
    "\n",
    "-   Nonlinear fitting is essentially a non-linear optimization problem\n",
    "\n",
    "-   Unlike linear regression, where we directly compute the parameters using matrix algebra, we have to provide an initial guess and iterate to the solution\n",
    "\n",
    "-   Similar to using fsolve, we must define a function of the model\n",
    "    -   The function takes an independent variable, and parameters, f(x,a,b,&#x2026;)\n",
    "    -   The function should return a value of $y$ for every value of $x$\n",
    "    -   i.e. it should be vectorized\n",
    "\n",
    "-   It is possible to formulate these problems as nonlinear minimization of summed squared errors. See [this example](http://jkitchin.github.io/blog/2013/02/18/Nonlinear-curve-fitting/).\n",
    "\n",
    "-   The function [scipy.optimize.curve_fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) provides nonlinear fitting of models (functions) to data. \n",
    "\n",
    "\n",
    "Let's say we want to fit some other data to the function $$y=ax/(b+x)$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 297
    },
    "executionInfo": {
     "elapsed": 371,
     "status": "ok",
     "timestamp": 1651148540478,
     "user": {
      "displayName": "Zachary Ulissi",
      "userId": "07633171379186475882"
     },
     "user_tz": 240
    },
    "id": "aesf39C9Nudl",
    "outputId": "3123d135-df25-40b4-828b-b5ed5a2e318a"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'y')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAPwElEQVR4nO3df6xfd13H8edrbeeuOCjSami30UFKdQKmeFmIMzpE0jLN1gAxmyLOAIvKkESsboGAGcYNmpBoHOJEgpjA+JGlqdmwUbeFKA57Z2FlWwplgOudujLpTOTiuvH2j/vtuLvc237ves/39Hs/z0dyk+/5nE+/5/3ptz2vnPM5389NVSFJatcZfRcgSeqXQSBJjTMIJKlxBoEkNc4gkKTGre67gKVat25dbdq0qe8yJGms3H333d+sqvUL7Ru7INi0aRNTU1N9lyFJYyXJNxbb560hSWqcQSBJjTMIJKlxBoEkNc4gkKTGjd1TQ5LUmt37p9m19yAPHZ1hw9oJdm7bwo6tG5ft/Q0CSaeFrk9242r3/mmuveUAM8eeAGD66AzX3nIAYNn+frw1JKl3x09200dnKL53stu9f7rv0nq3a+/BJ0PguJljT7Br78FlO4ZBIKl3ozjZjauHjs4sqf3pMAgk9W4UJ7txtWHtxJLanw6DQFLvRnGyG1c7t21hYs2qp7RNrFnFzm1blu0YBoGk3o3iZDeudmzdyPWveTEb104QYOPaCa5/zYt9akjSynL8pOZTQwvbsXVjp38XBoGk00LXJzstzltDktQ4g0CSGmcQSFLjDAJJalxnQZDkw0keTvKlRfb/apJ7khxI8rkkP9lVLZKkxXV5RfARYPsJ9n8N+LmqejHwHuCmDmuRJC2is8dHq+qzSTadYP/n5mzeBZzTVS2SpMWdLt8jeCPwmcV2JrkKuArgvPPOG1VNkksjqwm9B0GSVzAbBD+zWJ+quonBraPJyckaUWlq3CjWgZdOB70+NZTkJcCHgMuq6pE+a5Hmc2lktaK3IEhyHnAL8GtV9eW+6pAW49LIakVnt4aSfBy4GFiX5DDwbmANQFV9EHgX8BzgA0kAHq+qya7qkZZqw9oJphc46bs0slaaLp8auuIk+98EvKmr40unaue2LU+ZIwCXRtbK1PtksXS6cmlktcIgkE7ApZHVAtcakqTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcX6PQItyCWapDQaBFuQSzFI7vDWkBbkEs9QOg0ALcglmqR0GgRa02FLLLsEsrTwGgRa0c9sWJtasekqbSzBLK5OTxVqQSzBL7TAItCiXYJba4K0hSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1rrMgSPLhJA8n+dIi+5PkT5McSnJPkpd2VYskaXFdXhF8BNh+gv2vBjYPfq4C/rzDWiRJi+gsCKrqs8B/n6DLZcBHa9ZdwNokz+2qHknSwvqcI9gIPDhn+/Cg7fskuSrJVJKpI0eOjKQ4SWrFWEwWV9VNVTVZVZPr16/vuxxJWlH6DIJp4Nw52+cM2iRJI9RnEOwB3jB4eujlwKNV9R891iNJTVrd1Rsn+ThwMbAuyWHg3cAagKr6IHAbcAlwCPg28Btd1SJJWlxnQVBVV5xkfwFv6er4kqThjMVksSSpOwaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEnDYIkb03y7FEUI0kavWGuCH4U2Jfkk0m2J0nXRUmSRuekQVBV7wQ2A38FXAl8JckfJ3lBx7VJkkZgqDmCqirgPwc/jwPPBj6d5H0d1iZJGoHVJ+uQ5G3AG4BvAh8CdlbVsSRnAF8Bfr/bEiVJXRrmiuCHgddU1baq+lRVHQOoqu8Cv3SiPziYUziY5FCSaxbYf16SO5LsT3JPkkue1igkSU/bMHME766qbyyy7/7F/lySVcCNwKuBC4Arklwwr9s7gU9W1VbgcuADwxYuSVoeXX6P4ELgUFU9UFWPATcDl83rU8AzB6+fBTzUYT2SpAV0GQQbgQfnbB8etM31h8DrkxwGbgPeutAbJbkqyVSSqSNHjnRRqyQ1q+9vFl8BfKSqzgEuAf5mMAn9FFV1U1VNVtXk+vXrR17kcbv3T3PRDbdz/jW3ctENt7N7/3RvtUjScjnpU0OnYBo4d872OYO2ud4IbAeoqn9JchawDni4w7qelt37p7n2lgPMHHsCgOmjM1x7ywEAdmydf6EjSeOjyyuCfcDmJOcnOZPZyeA98/r8O/BKgCQ/DpwFnJb3fnbtPfhkCBw3c+wJdu092FNFkrQ8OguCqnocuBrYC9zP7NNB9ya5Lsmlg25vB96c5IvAx4ErB19eO+08dHRmSe2SNC66vDVEVd3G7CTw3LZ3zXl9H3BRlzUslw1rJ5he4KS/Ye1ED9VI0vLpe7J4bOzctoWJNaue0jaxZhU7t23pqSJJWh6dXhGsJMcnhHftPchDR2fYsHaCndu2OFEsaewZBEuwY+tGT/ySVhxvDUlS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqXKdBkGR7koNJDiW5ZpE+v5zkviT3JvlYl/VIkr7f6q7eOMkq4EbgVcBhYF+SPVV135w+m4FrgYuq6ltJfqSreiRJC+vyiuBC4FBVPVBVjwE3A5fN6/Nm4Maq+hZAVT3cYT2SpAV0GQQbgQfnbB8etM31QuCFSf45yV1Jti/0RkmuSjKVZOrIkSMdlStJbep7sng1sBm4GLgC+Mska+d3qqqbqmqyqibXr18/4hIlaWXrMgimgXPnbJ8zaJvrMLCnqo5V1deALzMbDJKkEekyCPYBm5Ocn+RM4HJgz7w+u5m9GiDJOmZvFT3QYU2SpHk6C4Kqehy4GtgL3A98sqruTXJdkksH3fYCjyS5D7gD2FlVj3RVkyTp+6Wq+q5hSSYnJ2tqaqrvMiRprCS5u6omF9rX92SxJKlnBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGtdpECTZnuRgkkNJrjlBv9cmqSSTXdSxe/80F91wO+dfcysX3XA7u/dPd3EYSRpLq7t64ySrgBuBVwGHgX1J9lTVffP6nQ28Dfh8F3Xs3j/NtbccYObYEwBMH53h2lsOALBj68YuDilJY6XLK4ILgUNV9UBVPQbcDFy2QL/3AO8FvtNFEbv2HnwyBI6bOfYEu/Ye7OJwkjR2ugyCjcCDc7YPD9qelOSlwLlVdeuJ3ijJVUmmkkwdOXJkSUU8dHRmSe2S1JreJouTnAG8H3j7yfpW1U1VNVlVk+vXr1/ScTasnVhSuyS1pssgmAbOnbN9zqDtuLOBFwF3Jvk68HJgz3JPGO/ctoWJNaue0jaxZhU7t21ZzsNI0tjqbLIY2AdsTnI+swFwOfArx3dW1aPAuuPbSe4Efq+qppaziOMTwrv2HuShozNsWDvBzm1bnCiWpIHOgqCqHk9yNbAXWAV8uKruTXIdMFVVe7o69nw7tm70xC9Ji+jyioCqug24bV7buxbpe3GXtUiSFuY3iyWpcQaBJDXOIJCkxhkEktS4VFXfNSxJkiPAN07SbR3wzRGUc7px3G1pddzQ7thPZdzPq6oFv5E7dkEwjCRTVdXJSqanM8fdllbHDe2Ovatxe2tIkhpnEEhS41ZqENzUdwE9cdxtaXXc0O7YOxn3ipwjkCQNb6VeEUiShmQQSFLjxjYIkmxPcjDJoSTXLLD/B5J8YrD/80k2jb7Kbgwx9p9N8m9JHk/yuj5q7MIQ4/7dJPcluSfJPyZ5Xh91Lrchxv2bSQ4k+UKSf0pyQR91LreTjXtOv9cmqeX+XSZ9GeLzvjLJkcHn/YUkbzrlg1bV2P0wu6z1V4HnA2cCXwQumNfnt4EPDl5fDnyi77pHOPZNwEuAjwKv67vmEY77FcAPDl7/1kr4zIcc9zPnvL4U+Lu+6x7FuAf9zgY+C9wFTPZd94g+7yuBP1vO447rFcGFwKGqeqCqHgNuBi6b1+cy4K8Hrz8NvDJJRlhjV0469qr6elXdA3y3jwI7Msy476iqbw8272L2t+KNu2HG/T9zNp8BrIQnQIb5Pw7wHuC9wHdGWVyHhh33shrXINgIPDhn+/CgbcE+VfU48CjwnJFU161hxr4SLXXcbwQ+02lFozHUuJO8JclXgfcBvzOi2rp00nEneSlwblXdOsrCOjbsv/PXDm6BfjrJuQvsX5JxDQJpUUleD0wCu/quZVSq6saqegHwB8A7+66na0nOAN4PvL3vWnrwt8CmqnoJ8Pd8787H0zauQTANzE3BcwZtC/ZJshp4FvDISKrr1jBjX4mGGneSXwDeAVxaVf83otq6tNTP+2ZgR6cVjcbJxn028CLgziRfB14O7FkBE8Yn/byr6pE5/7Y/BPzUqR50XINgH7A5yflJzmR2Mnj+70DeA/z64PXrgNtrMNMy5oYZ+0p00nEn2Qr8BbMh8HAPNXZhmHFvnrP5i8BXRlhfV0447qp6tKrWVdWmqtrE7JzQpVU11U+5y2aYz/u5czYvBe4/5aP2PUt+CrPrlwBfZnaG/R2DtuuY/ccAcBbwKeAQ8K/A8/uueYRjfxmz9xb/l9mroHv7rnlE4/4H4L+ALwx+9vRd84jG/SfAvYMx3wH8RN81j2Lc8/reyQp4amjIz/v6wef9xcHn/WOnekyXmJCkxo3rrSFJ0jIxCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDALpFCV52WABsLOSPCPJvUle1Hdd0rD8Qpm0DJL8EbPfZp8ADlfV9T2XJA3NIJCWwWBdmH3Mrov/01X1RM8lSUPz1pC0PJ4D/BCzq2Ke1XMt0pJ4RSAtgyR7mF0C+nzguVV1dc8lSUNb3XcB0rhL8gbgWFV9LMkq4HNJfr6qbu+7NmkYXhFIUuOcI5CkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXH/D609HksVeqeYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "x = np.array([0.5, 0.387, 0.24, 0.136, 0.04, 0.011])\n",
    "y = np.array([1.255, 1.25, 1.189, 1.124, 0.783, 0.402])\n",
    "\n",
    "\n",
    "plt.plot(x,y,'o')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 297
    },
    "executionInfo": {
     "elapsed": 500,
     "status": "ok",
     "timestamp": 1651148681053,
     "user": {
      "displayName": "Zachary Ulissi",
      "userId": "07633171379186475882"
     },
     "user_tz": 240
    },
    "id": "B4XHykHh6JRa",
    "outputId": "c32a4555-d6f6-4f63-809c-0c95a1d1e6cf"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'y')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAdqElEQVR4nO3de3RW9Z3v8fc3N5JwC5qAQri2iI2Cg0ar1ar1UkQt0OuRczrWVstMq52utoeprrpsq106rWvsnHNqaxn1tE5brToW6SkaWxWrVRSUm8AEEEFyQcIlF0hCSPI9f+wn8BADPIFnPzvJ/rzW2mtfk+e7E9if7N9vX8zdERGR+MqKugAREYmWgkBEJOYUBCIiMacgEBGJOQWBiEjM5URdQG8VFxf7hAkToi5DRKRfefPNN3e6e0lP6/pdEEyYMIHly5dHXYaISL9iZluPtE5NQyIiMacgEBGJOQWBiEjMKQhERGJOQSAiEnMKAhGRmOt3l4+KiMTNwhXV3FtRSU19C6OLCpg/Ywpzpo9J2/dXEIhInxL2QS9S7tBxANpb4EBrz+P2/XAgMW5vYfWWHbyzeguf79zPm1mn8XL9NG57ag1A2n4uCgIR6TMWrqjmtqfW0HKgA4Dq+pa0H/SOqLMTDjQHQ9u+4GB8cL750PTB5S1JQ/Ph0+2tSctak5a1gHf0qqxpwLQsIAt+3j6Llzun0XKgg3srKhUEIjLw3FtReTAEunzgoOceHFDb9kFbUzDevzcx3zXuYbrrAN81HDzA7wvG7S29Lzh7EOQVQk4B5BZAbiHk5gfTg05JLEsMOQWH1nVNJ49zBiXWdW0zCHIK+Lu7/0oLebSRgyd169bUH0e9R6AgEJHM6GiH/Y2JoQlak6b3N8L+vczdu4LBOS0MtRaG0MJgWhhirQxpboH7SBzwm8A7U/tMy4K8IcEBOm9wcNDOGwIFI2B4abAstzBYnptYn1t4+LLcgqTlBYeW5RZAVnaoPzKAwUUl1Pdw0B9dVJC2z1AQiEhqOtqhtR5a6oNxaz20NgRDS9J0a0NwYG9tOHSwb20M/vI+hq/lZLHX82mikH2ezz7yafIC6nNLmDxpUnAQHzQkMR6aOLgPCcYH55OW5eSDWQZ+OOGZP2PKYc1lAAW52cyfMSVtn6EgEImb9jZo2QPNu6BlNzTvDuZ7HOqDcWt90LxyNFm5UFAE+cNh0DDIHwbDRiemk5YNGhpMDxqamO9aNpQ/vr2H2/7w9gcOevd8aioMlA7jXupqEtNVQyLSM/egaWVfXXBg31cH+3ZC887gAN+8KzG/KzHsDppWjiQ7L2g26RqGl8IpUxMH+KIPjvOHHzr4p+Gv7zlnDwazgXvV0HGaM31MqD8DBYFIX+MeNKfs3QFN22Hv+8EBfu/7sLcO9u0I1u2rC4aOtp6/T24hFBZD4UlQeDIUT4aCk4L5ghHBssKTgmUFI4Lp3MLIm1LCPujJB4UWBGb2MHAtsMPdz+xh/f8AvgsY0AR8zd1XhVWPSJ+wvwkaa6GpJjjIN9Um5msPHfT37uj5CpasHBhcEgxDRsLIMhhcfGjZ4OJgKCwODvJ5hZnfP+mXwjwj+BXwM+CRI6x/F7jE3feY2UxgAfDREOsRCVdbMzRUQcO2YNxYA43ViaEmGPY3fvDrBg2DoafC0FEw9qPBeEj3YWTQFJOlp8JI+oUWBO7+VzObcJT1rybNLgVKw6pFJC3a9vH8a8uoeOV1BjdXMSV/DxeXNDPadgYH/uZd3b7AgoP4sNFw8odh4iXB9LAxMPSUxMH/lOAqGJEI9ZU+ghuBZ4600szmAfMAxo0bl6maJG7cgzb33e/Cnndh9+Zgevdm2LMFmndyOXA5QC60tOdRXVvC+6MmMuojn4LhYxNDaTAMPRVy8qLdJ5EURB4EZvYJgiC46EjbuPsCgqYjysvLPUOlyUC1fy/s2gi73oFdm4JhZ2I++Yoay4JhpXDSRDj9Gn6xuoN1zUVs85FUeQk7GQYYYxoL+NunLotsd0ROVKRBYGbTgAeBme7e/bxa5MS07IG6DVD3X1BXeWjcWJW0kUHR2KDpZux5cNKH4KRJwVA07rC/6H/y6p/o6a+QdN7qLxKFyILAzMYBTwF/7+4boqpDBoD2tuAv/PfXHhp2rAs6abvkFEDJaTDhQig+LTFMhhETg2e9pGB0UQHVId/qLxKFMC8ffRS4FCg2syrg+0AugLs/ANwBnAz83ILrltvdvTysemSAaGsODvS1K6F2VTDsWA+dB4L1WblQcjpMuCi4vHJkGZRMCdruT/CKm0zc6i8ShTCvGpp7jPU3ATeF9fkyAHS0w461ULUcqt+E6rdgZ+WhB44VnASj/w4uuDm4+3VkWfBXfnZuKOVk4lZ/kShE3lksctDeOnjvVdj2RnDgr1l56MaqwmIYcw6UzYJTzwqGYWMyfhes7nqVgUhBINFwDy7JfO812PpqMN61KViXnRcc6M+5AUrLg6FofOSPPhAZqBQEkjlN22HzS7B5Cbz70qHO3PzhMO4CmP73MP5jQQjkDIq0VJE4URBIr/TqfbJtzbDl5eDA/86LULc+WF4wAiZeDBO+FQTAyDI9OkEkQgoCSVlK75NtqIINFcHw7kvBe1pz8oMD/lnXwaRL4ZRpOvCL9CEKAklZz++TbWfhM88wZ8922PAsbA+CgaLxQRv/5E/C+AtTvlZfRDJPQSApO3QHrXO6beOa7KVck7WUSW3b4eWs4K/+K++E064KbthS565Iv6AgkJSdP2wX5ze/wLVZS/lQVi0dbrzaeQaP53yaW781HwafHHWJInIcFARydK2NsPYPsOI3PNr2Bh3ZxtLOMh46cDXPdpxLS+4I7rlmqkJApB9TEMgHucPWv8GK38C6p+FAMxRPgSvv4rmsi/nRS7sPXjV0h+6sFen3FARyyP69sOpReP2B4OauvKEw7QvB9f1jzgEzZgIzL4i6UBFJJwWBQP02eGMBvPVraG2A0WfDp38JH5ml996KxICCIM62vQGv3Q/r/xjMl82C878Opefqih+RGFEQxNG2ZbDkbnjnheDxDhfcDOfNC17QIiKxoyCIk+q3YMk9sPE5KDwZrrwLyr+il6eLxJyCIA5qV8GSf4HKxcFzfi7/fnAGoAAQERQEA9veOvjLD2Dlb4ImoE/cDh/9B8gfFnVlItKHKAgGoo52WPYgvHh3cA/Ax/4JPv4dKCiKujIR6YMUBAPNlr/B4vnBKx4nfQJm/iR4abuIyBEoCAaKvTug4nuw5vHgRe3/7Tdw+rW6DFREjklBMBBUPgtP3wz7m+Dif4aLvqUbwUQkZQqC/qytGZ67HZY/BKOmwg3/D0Z+JOqqRKSfCe01UWb2sJntMLO3j7DezOx/m9kmM1ttZmeHVcuAVLMSFlwShMDHvgFffV4hICLHJcz3Bf4KuOoo62cCkxPDPOAXIdYycHR2wMv3wYOXBw+Ju34RfPJHetm7iBy30JqG3P2vZjbhKJvMBh5xdweWmlmRmZ3q7rVh1dTvtTbAE1+Gd56Hsjlw7U+h8KSoqxKRfi7KPoIxwLak+arEsg8EgZnNIzhrYNy4cRkprs+pfw9++wXYtRGu/bfgfcC6IkhE0iDMpqG0cfcF7l7u7uUlJSVRl5N5VW/Cv18OjTXwxf+E8i8rBEQkbaI8I6gGkh93WZpYJsnWPQ1P/QMMGRlcFVQyJeqKRGSAifKMYBFwfeLqofOBBvUPJHGHV/4NHr8eTpkKNz2vEBCRUIR2RmBmjwKXAsVmVgV8H8gFcPcHgMXA1cAmoBn4cli19DudnbD4O7D8YTjjMzDnF5CbH3VVIjJAhXnV0NxjrHfg5rA+v99yh2dvDULgwm/C5T+ArH7RlSMi/ZTuLO5rXrgL3vglXHALXPFDdQqLSOj0p2Zf8vJ98PK/BpeGfvJHCgERyQgFQV/x+gJ4/ocw9fNwzX0KARHJGAVBX7Dit/DMfJhyTdAxnJUddUUiEiMKgqitXQiLboFJl8LnHobs3KgrEpGYURBEqepN+M+boPQ8uO53ukRURCKhIIhKyx544gYYeirMfRTyBkddkYjElC4fjYI7LLwZmmrgKxV6gqiIREpBEIWlv4DKP8GMu6G0POpqRCTm1DSUaVVvwp/vCK4QOv/rUVcjIqIgyKjm3Yf6Bebcr3sFRKRPUNNQprjD0zdDU23QL1AwIuqKREQABUHmLP05VC6GGfdA6TlRVyMicpCahjKhdlXQL3D6tXD+16KuRkTkMAqCsLnDM9+F/CKY/TP1C4hIn6OmobCt/QO891rwwnn1C4hIH6QzgjAdaAmahEZNhbOvj7oaEZEe6YwgTK/+H2jYBp9+QE8UFZE+S2cEYWmohld+CmWzYcJFUVcjInJECoKw/OUH0NkBV94VdSUiIkelIAjDtjdgzePwsW/AiPFRVyMiclQKgnTr7AwuFx1yClz0rairERE5plCDwMyuMrNKM9tkZrf2sH6cmb1oZivMbLWZXR1mPRmx+jGoeQuu/CEMGhJ1NSIixxRaEJhZNnA/MBMoA+aaWVm3zW4HHnf36cB1wM/Dqicj9jfBX34IY86BqV+IuhoRkZSEeUZwHrDJ3Te7exvwGDC72zYODEtMDwdqQqwnfEsfgL3b4aofQ5Za3USkfwjzaDUG2JY0X5VYluwHwBfNrApYDHyjp29kZvPMbLmZLa+rqwuj1hP29JvvsvPF+1nScRYX/raJhSuqoy5JRCQlUf/ZOhf4lbuXAlcD/2FmH6jJ3Re4e7m7l5eUlGS8yGNZuKKal57+vxSzh193fJLq+hZue2qNwkBE+oUwg6AaGJs0X5pYluxG4HEAd38NyAeKQ6wpFPdWVDKXZ9jaOZIlnWcB0HKgg3srKiOuTETk2MIMgmXAZDObaGZ5BJ3Bi7pt8x5wOYCZfYQgCPpm289RFDWs59ysDTzScSWe9COtqW+JsCoRkdSEFgTu3g7cAlQA6wmuDlprZnea2azEZt8Bvmpmq4BHgRvc3cOqKSz/WPg8zT6IJzouPWz56KKCaAoSEemFUB865+6LCTqBk5fdkTS9DrgwzBpC17ybq/0VnvCP08jgg4sLcrOZP2NKhIWJiKQm6s7i/u+tR8ju3M/Iy29hTFEBBowpKuCez0xlzvTuF0mJiPQ9egz1iejsgGUPwYSPc9kln+CyS6IuSESk93RGcCIqn4GG9+C8eVFXIiJy3BQEJ+KNBTCsFKb0/0ckiUh8KQiO147/gndfgnO/AtlqYROR/ktBcLzeWADZg+DsL0VdiYjICVEQHI/WBlj1GJz5WRjc726EFhE5jILgeKz8HRzYBx9VJ7GI9H8KguOx6jEYPT0YRET6OQVBbzXWQu1KOP3aqCsREUkLBUFvbXwuGJ92VbR1iIikiYKgtzY+F9w7MOqMqCsREUkLBUFvHGiFd16E02aAWdTViIikhYKgN7a+ElwtdNqMqCsREUkbBUFvbKiAnAKYeHHUlYiIpI2CIFXuQRBMugRy9cIZERk4FASpqquE+q1qFhKRAUdBkKoNzwbjyZ+Mtg4RkTRTEKRq43MwaioML426EhGRtFIQpKJ5N7y3VM1CIjIgHTMIzOwbZjYiE8X0We+8AN6hIBCRASmVM4JRwDIze9zMrjKL4Z1UGyqg8GQYc07UlYiIpN0xg8DdbwcmAw8BNwAbzexuM/tQyLX1DR3tsOnPQSdxVnbU1YiIpF1KfQTu7sD2xNAOjACeNLOfHO3rEmcQlWa2ycxuPcI2XzCzdWa21sx+18v6w1e1DFr2qFlIRAasY75s18y+CVwP7AQeBOa7+wEzywI2Av98hK/LBu4HrgSqCJqXFrn7uqRtJgO3ARe6+x4zG3miO5R2GysgKwc+dFnUlYiIhCKVt66fBHzG3bcmL3T3TjM72kP5zwM2uftmADN7DJgNrEva5qvA/e6+J/E9d/Sm+IzYUAHjLoD84VFXIiISilT6CL7fPQSS1q0/ypeOAbYlzVclliU7DTjNzP5mZkvNrMeH/JvZPDNbbmbL6+rqjlVy+tS/BzvW6d0DIjKgRX0fQQ5BR/SlwFzg382sqPtG7r7A3cvdvbykpCRz1W2oCMYKAhEZwMIMgmpgbNJ8aWJZsipgkbsfcPd3gQ0EwdA3bKiAkyZB8YejrkREJDRhBsEyYLKZTTSzPOA6YFG3bRYSnA1gZsUETUWbQ6wpdZ0dsPVV+PAVUVciIhKq0ILA3duBW4AKYD3wuLuvNbM7zWxWYrMKYJeZrQNeJLgiaVdYNfXKrneCl9CMnh51JSIioUrlqqHj5u6LgcXdlt2RNO3AtxND37J9dTA+ZWq0dYiIhCzqzuK+a/tqyM6D4ilRVyIiEioFwZFsXwMlp0NOXtSViIiESkHQE3eoXQ2nTou6EhGR0CkIetJUC8074ZSzoq5ERCR0CoKebF8TjNVRLCIxoCDoSe1qwOCUM6OuREQkdAqCnmxfFdxRPGho1JWIiIROQdCT7WvULCQisaEg6K61AfZs0RVDIhIbCoLuDnYUKwhEJB4UBN0pCEQkZhQE3dWuhiGjYOioqCsREckIBUF36igWkZhRECRr3w9169UsJCKxoiBItmM9dLbriiERiRUFQTJ1FItIDCkIkm1fDXlDYMTEqCsREckYBUGy2tUw6kzI0o9FROJDR7wunZ3w/tvqHxCR2FEQdNnzLrTtVf+AiMSOgqCLXlYvIjGlIOhSuxqycmDkR6KuREQko0INAjO7yswqzWyTmd16lO0+a2ZuZuVh1tOThSuqufBfXmDJS8+z0UtZuGZnpksQEYlUaEFgZtnA/cBMoAyYa2ZlPWw3FPgm8HpYtRzJwhXV3PbUGqrrWyjL2sqq9nHc9tQaFq6oznQpIiKRCfOM4Dxgk7tvdvc24DFgdg/b3QX8GGgNsZYe3VtRScuBDkqoZ6TVs7ZzPC0HOri3ojLTpYiIRCbMIBgDbEuar0osO8jMzgbGuvufjvaNzGyemS03s+V1dXVpK7CmvgWAM7K2ALCuc8Jhy0VE4iCyzmIzywLuA75zrG3dfYG7l7t7eUlJSdpqGF1UAECZbQVgnY8/bLmISByEGQTVwNik+dLEsi5DgTOBJWa2BTgfWJTJDuP5M6ZQkJtNWdYWtnaOpIlCCnKzmT9jSqZKEBGJXE6I33sZMNnMJhIEwHXAf+9a6e4NQHHXvJktAf6nuy8PsabDzJketFSdteg93m4fz5iiAubPmHJwuYhIHIQWBO7ebma3ABVANvCwu681szuB5e6+KKzP7o05ZcPg6VrGXnEjMy+5LOpyREQyLswzAtx9MbC427I7jrDtpWHWckR7tgTj4smRfLyISNR0Z3FjbTAepuYgEYknBUFjov962KnR1iEiEhEFQVMtWBYMGRV1JSIikVAQNNbA4JGQnRt1JSIikVAQNNaoWUhEYk1B0FSrjmIRiTUFQWM1DNUZgYjEV7yDoG0ftDbAsNFRVyIiEpl4B8HBewgUBCISX/EOgqaaYKwgEJEYi3cQNCaCYKiCQETiS0EAunxURGJNQZA/HPIGR12JiEhk4h0EuodARCTmQdBYo3sIRCT2FATqHxCRmItvEHQcgL3vq2lIRGIvvkGw933A1TQkIrEX3yDQm8lERIBYB4HeTCYiAnEOgiadEYiIQJyDoLEasgdBwYioKxERiVSoQWBmV5lZpZltMrNbe1j/bTNbZ2arzex5MxsfZj2HaawNHjZnlrGPFBHpi0ILAjPLBu4HZgJlwFwzK+u22Qqg3N2nAU8CPwmrng9oqtVTR0VECPeM4Dxgk7tvdvc24DFgdvIG7v6iuzcnZpcCpSHWc7jGagWBiAjhBsEYYFvSfFVi2ZHcCDwTYj2HuAdNQ7qHQESEnKgLADCzLwLlwCVHWD8PmAcwbty4E//A5t3QsV9nBCIihHtGUA2MTZovTSw7jJldAXwPmOXu+3v6Ru6+wN3L3b28pKTkxCvTm8lERA4KMwiWAZPNbKKZ5QHXAYuSNzCz6cAvCUJgR4i1HE5vJhMROSi0IHD3duAWoAJYDzzu7mvN7E4zm5XY7F5gCPCEma00s0VH+Hbp1agzAhGRLqH2Ebj7YmBxt2V3JE1fEebnH1FjDVgWDBkVyceLiPQl8byzuKkmCIHsPtFXLiISqXgGgS4dFRE5KKZBUKP+ARGRhHgGQZOCQESkS/yCoG0ftDYoCEREEuIXBF1vJtM9BCIiQCyDQG8mExFJFr8g0JvJREQOE78g6Doj0OWjIiJALIOgFvKLIK8w6kpERPqE+AWB3kwmInKY+AWB3kwmInKYGAaBHi8hIpIsXkHQcQD2vq8rhkREksQrCPa+D7juIRARSRKvIDj4QhqdEYiIdIlnEKiPQETkoHgGga4aEhE5KF5B0FQDOflQMCLqSkRE+ox4BUHXpaNmUVciItJnxCwIatRRLCLSTbyCoKlGl46KiHQTnyBwD5qG1FEsInKYUIPAzK4ys0oz22Rmt/awfpCZ/T6x/nUzmxBaMc27oWO/3kwmItJNaEFgZtnA/cBMoAyYa2Zl3Ta7Edjj7h8Gfgr8OKx6Dr2ZTEEgIpIszDOC84BN7r7Z3duAx4DZ3baZDfw6Mf0kcLlZSJf0HHwzmYJARCRZmEEwBtiWNF+VWNbjNu7eDjQAJ3f/RmY2z8yWm9nyurq646smfzicfi0UjTu+rxcRGaByoi4gFe6+AFgAUF5e7sf1TcadHwwiInKYMM8IqoGxSfOliWU9bmNmOcBwYFeINYmISDdhBsEyYLKZTTSzPOA6YFG3bRYBX0pMfw54wd2P7y9+ERE5LqE1Dbl7u5ndAlQA2cDD7r7WzO4Elrv7IuAh4D/MbBOwmyAsREQkg0LtI3D3xcDibsvuSJpuBT4fZg0iInJ08bmzWEREeqQgEBGJOQWBiEjMKQhERGLO+tvVmmZWB2w9zi8vBnamsZz+QPscD9rneDiRfR7v7iU9reh3QXAizGy5u5dHXUcmaZ/jQfscD2Hts5qGRERiTkEgIhJzcQuCBVEXEAHtczxon+MhlH2OVR+BiIh8UNzOCEREpBsFgYhIzA3IIDCzq8ys0sw2mdmtPawfZGa/T6x/3cwmZL7K9Ephny82s7fMrN3MPhdFjemWwj5/28zWmdlqM3vezMZHUWc6pbDP/2hma8xspZm90sN7wvudY+1z0nafNTM3s359SWkKv+MbzKwu8TteaWY3nfCHuvuAGggeef0OMAnIA1YBZd22+TrwQGL6OuD3UdedgX2eAEwDHgE+F3XNGdrnTwCFiemvxeT3PCxpehbwbNR1h73Pie2GAn8FlgLlUdcd8u/4BuBn6fzcgXhGcB6wyd03u3sb8Bgwu9s2s4FfJ6afBC43M8tgjel2zH129y3uvhrojKLAEKSyzy+6e3NidinBW/L6s1T2uTFpdjDQ368GSeX/M8BdwI+B1kwWF4JU9zetBmIQjAG2Jc1XJZb1uI27twMNwMkZqS4cqezzQNPbfb4ReCbUisKX0j6b2c1m9g7wE+CfMlRbWI65z2Z2NjDW3f+UycJCkuq/688mmjyfNLOxPazvlYEYBCKHMbMvAuXAvVHXkgnufr+7fwj4LnB71PWEycyygPuA70RdSwb9EZjg7tOAP3OodeO4DcQgqAaSE7I0sazHbcwsBxgO7MpIdeFIZZ8HmpT22cyuAL4HzHL3/RmqLSy9/T0/BswJtaLwHWufhwJnAkvMbAtwPrCoH3cYH/N37O67kv4tPwicc6IfOhCDYBkw2cwmmlkeQWfwom7bLAK+lJj+HPCCJ3ph+qlU9nmgOeY+m9l04JcEIbAjghrTLZV9npw0ew2wMYP1heGo++zuDe5e7O4T3H0CQV/QLHdfHk25JyyV3/GpSbOzgPUn/KlR95KH1PN+NbCBoPf9e4lldxL8AwHIB54ANgFvAJOirjkD+3wuQXvjPoKzn7VR15yBff4L8D6wMjEsirrmDOzz/wLWJvb3ReCMqGsOe5+7bbuEfnzVUIq/43sSv+NVid/x6Sf6mXrEhIhIzA3EpiEREekFBYGISMwpCEREYk5BICIScwoCEZGYUxCIiMScgkBEJOYUBCInyMzOTTwALN/MBpvZWjM7M+q6RFKlG8pE0sDMfkRwx3oBUOXu90RckkjKFAQiaZB4Lswygufhf8zdOyIuSSRlahoSSY+TgSEET8PMj7gWkV7RGYFIGpjZIoLHPk8ETnX3WyIuSSRlOVEXINLfmdn1wAF3/52ZZQOvmtll7v5C1LWJpEJnBCIiMac+AhGRmFMQiIjEnIJARCTmFAQiIjGnIBARiTkFgYhIzCkIRERi7v8DnwEJQ/BxzncAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "def func(x, a, b):\n",
    "    'nonlinear function in a and b to fit to data'\n",
    "    return a * x / (b + x)\n",
    "\n",
    "popt, pcov = curve_fit(func, x, y, p0=(3,3))\n",
    "\n",
    "\n",
    "xrange = np.linspace(0,0.5)\n",
    "fitted_y = func(xrange, *popt)\n",
    "\n",
    "plt.plot(x,y,'o')\n",
    "plt.plot(xrange,fitted_y)\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 168,
     "status": "ok",
     "timestamp": 1651148614348,
     "user": {
      "displayName": "Zachary Ulissi",
      "userId": "07633171379186475882"
     },
     "user_tz": 240
    },
    "id": "XO46nvzjNfKT",
    "outputId": "1ec2e2d9-b5e9-4339-bc49-de03212c1fa2"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.32753142 0.02646156]\n"
     ]
    }
   ],
   "source": [
    "print(popt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yZNeInig2TnP"
   },
   "source": [
    "-   We also need to estimate uncertainties in nonlinear parameters\n",
    "\n",
    "-   `lmfit` provides a nice way to do this\n",
    "\n",
    "[lmfit](https://lmfit.github.io/lmfit-py/)\n",
    "\n",
    "Read the [lmfit](https://lmfit.github.io/lmfit-py/) documentation to see how the confidence intervals are computed\n",
    "\n",
    "Here is an example usage of lmfit.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 10920,
     "status": "ok",
     "timestamp": 1651148799014,
     "user": {
      "displayName": "Zachary Ulissi",
      "userId": "07633171379186475882"
     },
     "user_tz": 240
    },
    "id": "9hUADZ3w6xEu",
    "outputId": "3ed3fd3a-1f16-435c-d279-98abfee0f57d"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collecting lmfit\n",
      "  Downloading lmfit-1.0.3.tar.gz (292 kB)\n",
      "\u001b[K     |████████████████████████████████| 292 kB 7.7 MB/s \n",
      "\u001b[?25hCollecting asteval>=0.9.22\n",
      "  Downloading asteval-0.9.26.tar.gz (40 kB)\n",
      "\u001b[K     |████████████████████████████████| 40 kB 4.1 MB/s \n",
      "\u001b[?25hRequirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.7/dist-packages (from lmfit) (1.21.6)\n",
      "Requirement already satisfied: scipy>=1.4 in /usr/local/lib/python3.7/dist-packages (from lmfit) (1.4.1)\n",
      "Collecting uncertainties>=3.0.1\n",
      "  Downloading uncertainties-3.1.6-py2.py3-none-any.whl (98 kB)\n",
      "\u001b[K     |████████████████████████████████| 98 kB 7.1 MB/s \n",
      "\u001b[?25hRequirement already satisfied: importlib_metadata in /usr/local/lib/python3.7/dist-packages (from asteval>=0.9.22->lmfit) (4.11.3)\n",
      "Requirement already satisfied: future in /usr/local/lib/python3.7/dist-packages (from uncertainties>=3.0.1->lmfit) (0.16.0)\n",
      "Requirement already satisfied: zipp>=0.5 in /usr/local/lib/python3.7/dist-packages (from importlib_metadata->asteval>=0.9.22->lmfit) (3.8.0)\n",
      "Requirement already satisfied: typing-extensions>=3.6.4 in /usr/local/lib/python3.7/dist-packages (from importlib_metadata->asteval>=0.9.22->lmfit) (4.2.0)\n",
      "Building wheels for collected packages: lmfit, asteval\n",
      "  Building wheel for lmfit (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
      "  Created wheel for lmfit: filename=lmfit-1.0.3-py3-none-any.whl size=84402 sha256=8328dc5ec51b17be5d13b3ab47d4b414f6a54deab3dd98f55708967677867dff\n",
      "  Stored in directory: /root/.cache/pip/wheels/b9/7a/d1/236aa0f8196b264fda481a112f7cfb1bfde7bfb20235f8e331\n",
      "  Building wheel for asteval (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
      "  Created wheel for asteval: filename=asteval-0.9.26-py3-none-any.whl size=17648 sha256=7d66fc626e5eeb588f3a81b94e0560796c7f9f0602982830b1c2e5983d01719d\n",
      "  Stored in directory: /root/.cache/pip/wheels/4c/e9/f0/bc343d5b77d2fded45177f424a6b0b9224b92ff6e7c150bad4\n",
      "Successfully built lmfit asteval\n",
      "Installing collected packages: uncertainties, asteval, lmfit\n",
      "Successfully installed asteval-0.9.26 lmfit-1.0.3 uncertainties-3.1.6\n"
     ]
    }
   ],
   "source": [
    "!pip install lmfit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 583
    },
    "executionInfo": {
     "elapsed": 344,
     "status": "ok",
     "timestamp": 1651148861738,
     "user": {
      "displayName": "Zachary Ulissi",
      "userId": "07633171379186475882"
     },
     "user_tz": 240
    },
    "id": "FOBWKAJzFNHe",
    "outputId": "0990defa-db2c-4c34-f0d5-367e4c4edf84"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Model]]\n",
      "    Model(func)\n",
      "[[Fit Statistics]]\n",
      "    # fitting method   = leastsq\n",
      "    # function evals   = 36\n",
      "    # data points      = 6\n",
      "    # variables        = 2\n",
      "    chi-square         = 6.9885e-04\n",
      "    reduced chi-square = 1.7471e-04\n",
      "    Akaike info crit   = -50.3470350\n",
      "    Bayesian info crit = -50.7635160\n",
      "[[Variables]]\n",
      "    a:  1.32753139 +/- 0.00972276 (0.73%) (init = 2)\n",
      "    b:  0.02646155 +/- 0.00102789 (3.88%) (init = 1)\n",
      "[[Correlations]] (unreported correlations are < 0.100)\n",
      "    C(a, b) = 0.711\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'y')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAdqElEQVR4nO3de3RW9Z3v8fc3N5JwC5qAQri2iI2Cg0ar1ar1UkQt0OuRczrWVstMq52utoeprrpsq106rWvsnHNqaxn1tE5brToW6SkaWxWrVRSUm8AEEEFyQcIlF0hCSPI9f+wn8BADPIFnPzvJ/rzW2mtfk+e7E9if7N9vX8zdERGR+MqKugAREYmWgkBEJOYUBCIiMacgEBGJOQWBiEjM5URdQG8VFxf7hAkToi5DRKRfefPNN3e6e0lP6/pdEEyYMIHly5dHXYaISL9iZluPtE5NQyIiMacgEBGJOQWBiEjMKQhERGJOQSAiEnMKAhGRmOt3l4+KiMTNwhXV3FtRSU19C6OLCpg/Ywpzpo9J2/dXEIhInxL2QS9S7tBxANpb4EBrz+P2/XAgMW5vYfWWHbyzeguf79zPm1mn8XL9NG57ag1A2n4uCgIR6TMWrqjmtqfW0HKgA4Dq+pa0H/SOqLMTDjQHQ9u+4GB8cL750PTB5S1JQ/Ph0+2tSctak5a1gHf0qqxpwLQsIAt+3j6Llzun0XKgg3srKhUEIjLw3FtReTAEunzgoOceHFDb9kFbUzDevzcx3zXuYbrrAN81HDzA7wvG7S29Lzh7EOQVQk4B5BZAbiHk5gfTg05JLEsMOQWH1nVNJ49zBiXWdW0zCHIK+Lu7/0oLebSRgyd169bUH0e9R6AgEJHM6GiH/Y2JoQlak6b3N8L+vczdu4LBOS0MtRaG0MJgWhhirQxpboH7SBzwm8A7U/tMy4K8IcEBOm9wcNDOGwIFI2B4abAstzBYnptYn1t4+LLcgqTlBYeW5RZAVnaoPzKAwUUl1Pdw0B9dVJC2z1AQiEhqOtqhtR5a6oNxaz20NgRDS9J0a0NwYG9tOHSwb20M/vI+hq/lZLHX82mikH2ezz7yafIC6nNLmDxpUnAQHzQkMR6aOLgPCcYH55OW5eSDWQZ+OOGZP2PKYc1lAAW52cyfMSVtn6EgEImb9jZo2QPNu6BlNzTvDuZ7HOqDcWt90LxyNFm5UFAE+cNh0DDIHwbDRiemk5YNGhpMDxqamO9aNpQ/vr2H2/7w9gcOevd8aioMlA7jXupqEtNVQyLSM/egaWVfXXBg31cH+3ZC887gAN+8KzG/KzHsDppWjiQ7L2g26RqGl8IpUxMH+KIPjvOHHzr4p+Gv7zlnDwazgXvV0HGaM31MqD8DBYFIX+MeNKfs3QFN22Hv+8EBfu/7sLcO9u0I1u2rC4aOtp6/T24hFBZD4UlQeDIUT4aCk4L5ghHBssKTgmUFI4Lp3MLIm1LCPujJB4UWBGb2MHAtsMPdz+xh/f8AvgsY0AR8zd1XhVWPSJ+wvwkaa6GpJjjIN9Um5msPHfT37uj5CpasHBhcEgxDRsLIMhhcfGjZ4OJgKCwODvJ5hZnfP+mXwjwj+BXwM+CRI6x/F7jE3feY2UxgAfDREOsRCVdbMzRUQcO2YNxYA43ViaEmGPY3fvDrBg2DoafC0FEw9qPBeEj3YWTQFJOlp8JI+oUWBO7+VzObcJT1rybNLgVKw6pFJC3a9vH8a8uoeOV1BjdXMSV/DxeXNDPadgYH/uZd3b7AgoP4sNFw8odh4iXB9LAxMPSUxMH/lOAqGJEI9ZU+ghuBZ4600szmAfMAxo0bl6maJG7cgzb33e/Cnndh9+Zgevdm2LMFmndyOXA5QC60tOdRXVvC+6MmMuojn4LhYxNDaTAMPRVy8qLdJ5EURB4EZvYJgiC46EjbuPsCgqYjysvLPUOlyUC1fy/s2gi73oFdm4JhZ2I++Yoay4JhpXDSRDj9Gn6xuoN1zUVs85FUeQk7GQYYYxoL+NunLotsd0ROVKRBYGbTgAeBme7e/bxa5MS07IG6DVD3X1BXeWjcWJW0kUHR2KDpZux5cNKH4KRJwVA07rC/6H/y6p/o6a+QdN7qLxKFyILAzMYBTwF/7+4boqpDBoD2tuAv/PfXHhp2rAs6abvkFEDJaTDhQig+LTFMhhETg2e9pGB0UQHVId/qLxKFMC8ffRS4FCg2syrg+0AugLs/ANwBnAz83ILrltvdvTysemSAaGsODvS1K6F2VTDsWA+dB4L1WblQcjpMuCi4vHJkGZRMCdruT/CKm0zc6i8ShTCvGpp7jPU3ATeF9fkyAHS0w461ULUcqt+E6rdgZ+WhB44VnASj/w4uuDm4+3VkWfBXfnZuKOVk4lZ/kShE3lksctDeOnjvVdj2RnDgr1l56MaqwmIYcw6UzYJTzwqGYWMyfhes7nqVgUhBINFwDy7JfO812PpqMN61KViXnRcc6M+5AUrLg6FofOSPPhAZqBQEkjlN22HzS7B5Cbz70qHO3PzhMO4CmP73MP5jQQjkDIq0VJE4URBIr/TqfbJtzbDl5eDA/86LULc+WF4wAiZeDBO+FQTAyDI9OkEkQgoCSVlK75NtqIINFcHw7kvBe1pz8oMD/lnXwaRL4ZRpOvCL9CEKAklZz++TbWfhM88wZ8922PAsbA+CgaLxQRv/5E/C+AtTvlZfRDJPQSApO3QHrXO6beOa7KVck7WUSW3b4eWs4K/+K++E064KbthS565Iv6AgkJSdP2wX5ze/wLVZS/lQVi0dbrzaeQaP53yaW781HwafHHWJInIcFARydK2NsPYPsOI3PNr2Bh3ZxtLOMh46cDXPdpxLS+4I7rlmqkJApB9TEMgHucPWv8GK38C6p+FAMxRPgSvv4rmsi/nRS7sPXjV0h+6sFen3FARyyP69sOpReP2B4OauvKEw7QvB9f1jzgEzZgIzL4i6UBFJJwWBQP02eGMBvPVraG2A0WfDp38JH5ml996KxICCIM62vQGv3Q/r/xjMl82C878Opefqih+RGFEQxNG2ZbDkbnjnheDxDhfcDOfNC17QIiKxoyCIk+q3YMk9sPE5KDwZrrwLyr+il6eLxJyCIA5qV8GSf4HKxcFzfi7/fnAGoAAQERQEA9veOvjLD2Dlb4ImoE/cDh/9B8gfFnVlItKHKAgGoo52WPYgvHh3cA/Ax/4JPv4dKCiKujIR6YMUBAPNlr/B4vnBKx4nfQJm/iR4abuIyBEoCAaKvTug4nuw5vHgRe3/7Tdw+rW6DFREjklBMBBUPgtP3wz7m+Dif4aLvqUbwUQkZQqC/qytGZ67HZY/BKOmwg3/D0Z+JOqqRKSfCe01UWb2sJntMLO3j7DezOx/m9kmM1ttZmeHVcuAVLMSFlwShMDHvgFffV4hICLHJcz3Bf4KuOoo62cCkxPDPOAXIdYycHR2wMv3wYOXBw+Ju34RfPJHetm7iBy30JqG3P2vZjbhKJvMBh5xdweWmlmRmZ3q7rVh1dTvtTbAE1+Gd56Hsjlw7U+h8KSoqxKRfi7KPoIxwLak+arEsg8EgZnNIzhrYNy4cRkprs+pfw9++wXYtRGu/bfgfcC6IkhE0iDMpqG0cfcF7l7u7uUlJSVRl5N5VW/Cv18OjTXwxf+E8i8rBEQkbaI8I6gGkh93WZpYJsnWPQ1P/QMMGRlcFVQyJeqKRGSAifKMYBFwfeLqofOBBvUPJHGHV/4NHr8eTpkKNz2vEBCRUIR2RmBmjwKXAsVmVgV8H8gFcPcHgMXA1cAmoBn4cli19DudnbD4O7D8YTjjMzDnF5CbH3VVIjJAhXnV0NxjrHfg5rA+v99yh2dvDULgwm/C5T+ArH7RlSMi/ZTuLO5rXrgL3vglXHALXPFDdQqLSOj0p2Zf8vJ98PK/BpeGfvJHCgERyQgFQV/x+gJ4/ocw9fNwzX0KARHJGAVBX7Dit/DMfJhyTdAxnJUddUUiEiMKgqitXQiLboFJl8LnHobs3KgrEpGYURBEqepN+M+boPQ8uO53ukRURCKhIIhKyx544gYYeirMfRTyBkddkYjElC4fjYI7LLwZmmrgKxV6gqiIREpBEIWlv4DKP8GMu6G0POpqRCTm1DSUaVVvwp/vCK4QOv/rUVcjIqIgyKjm3Yf6Bebcr3sFRKRPUNNQprjD0zdDU23QL1AwIuqKREQABUHmLP05VC6GGfdA6TlRVyMicpCahjKhdlXQL3D6tXD+16KuRkTkMAqCsLnDM9+F/CKY/TP1C4hIn6OmobCt/QO891rwwnn1C4hIH6QzgjAdaAmahEZNhbOvj7oaEZEe6YwgTK/+H2jYBp9+QE8UFZE+S2cEYWmohld+CmWzYcJFUVcjInJECoKw/OUH0NkBV94VdSUiIkelIAjDtjdgzePwsW/AiPFRVyMiclQKgnTr7AwuFx1yClz0rairERE5plCDwMyuMrNKM9tkZrf2sH6cmb1oZivMbLWZXR1mPRmx+jGoeQuu/CEMGhJ1NSIixxRaEJhZNnA/MBMoA+aaWVm3zW4HHnf36cB1wM/Dqicj9jfBX34IY86BqV+IuhoRkZSEeUZwHrDJ3Te7exvwGDC72zYODEtMDwdqQqwnfEsfgL3b4aofQ5Za3USkfwjzaDUG2JY0X5VYluwHwBfNrApYDHyjp29kZvPMbLmZLa+rqwuj1hP29JvvsvPF+1nScRYX/raJhSuqoy5JRCQlUf/ZOhf4lbuXAlcD/2FmH6jJ3Re4e7m7l5eUlGS8yGNZuKKal57+vxSzh193fJLq+hZue2qNwkBE+oUwg6AaGJs0X5pYluxG4HEAd38NyAeKQ6wpFPdWVDKXZ9jaOZIlnWcB0HKgg3srKiOuTETk2MIMgmXAZDObaGZ5BJ3Bi7pt8x5wOYCZfYQgCPpm289RFDWs59ysDTzScSWe9COtqW+JsCoRkdSEFgTu3g7cAlQA6wmuDlprZnea2azEZt8Bvmpmq4BHgRvc3cOqKSz/WPg8zT6IJzouPWz56KKCaAoSEemFUB865+6LCTqBk5fdkTS9DrgwzBpC17ybq/0VnvCP08jgg4sLcrOZP2NKhIWJiKQm6s7i/u+tR8ju3M/Iy29hTFEBBowpKuCez0xlzvTuF0mJiPQ9egz1iejsgGUPwYSPc9kln+CyS6IuSESk93RGcCIqn4GG9+C8eVFXIiJy3BQEJ+KNBTCsFKb0/0ckiUh8KQiO147/gndfgnO/AtlqYROR/ktBcLzeWADZg+DsL0VdiYjICVEQHI/WBlj1GJz5WRjc726EFhE5jILgeKz8HRzYBx9VJ7GI9H8KguOx6jEYPT0YRET6OQVBbzXWQu1KOP3aqCsREUkLBUFvbXwuGJ92VbR1iIikiYKgtzY+F9w7MOqMqCsREUkLBUFvHGiFd16E02aAWdTViIikhYKgN7a+ElwtdNqMqCsREUkbBUFvbKiAnAKYeHHUlYiIpI2CIFXuQRBMugRy9cIZERk4FASpqquE+q1qFhKRAUdBkKoNzwbjyZ+Mtg4RkTRTEKRq43MwaioML426EhGRtFIQpKJ5N7y3VM1CIjIgHTMIzOwbZjYiE8X0We+8AN6hIBCRASmVM4JRwDIze9zMrjKL4Z1UGyqg8GQYc07UlYiIpN0xg8DdbwcmAw8BNwAbzexuM/tQyLX1DR3tsOnPQSdxVnbU1YiIpF1KfQTu7sD2xNAOjACeNLOfHO3rEmcQlWa2ycxuPcI2XzCzdWa21sx+18v6w1e1DFr2qFlIRAasY75s18y+CVwP7AQeBOa7+wEzywI2Av98hK/LBu4HrgSqCJqXFrn7uqRtJgO3ARe6+x4zG3miO5R2GysgKwc+dFnUlYiIhCKVt66fBHzG3bcmL3T3TjM72kP5zwM2uftmADN7DJgNrEva5qvA/e6+J/E9d/Sm+IzYUAHjLoD84VFXIiISilT6CL7fPQSS1q0/ypeOAbYlzVclliU7DTjNzP5mZkvNrMeH/JvZPDNbbmbL6+rqjlVy+tS/BzvW6d0DIjKgRX0fQQ5BR/SlwFzg382sqPtG7r7A3cvdvbykpCRz1W2oCMYKAhEZwMIMgmpgbNJ8aWJZsipgkbsfcPd3gQ0EwdA3bKiAkyZB8YejrkREJDRhBsEyYLKZTTSzPOA6YFG3bRYSnA1gZsUETUWbQ6wpdZ0dsPVV+PAVUVciIhKq0ILA3duBW4AKYD3wuLuvNbM7zWxWYrMKYJeZrQNeJLgiaVdYNfXKrneCl9CMnh51JSIioUrlqqHj5u6LgcXdlt2RNO3AtxND37J9dTA+ZWq0dYiIhCzqzuK+a/tqyM6D4ilRVyIiEioFwZFsXwMlp0NOXtSViIiESkHQE3eoXQ2nTou6EhGR0CkIetJUC8074ZSzoq5ERCR0CoKebF8TjNVRLCIxoCDoSe1qwOCUM6OuREQkdAqCnmxfFdxRPGho1JWIiIROQdCT7WvULCQisaEg6K61AfZs0RVDIhIbCoLuDnYUKwhEJB4UBN0pCEQkZhQE3dWuhiGjYOioqCsREckIBUF36igWkZhRECRr3w9169UsJCKxoiBItmM9dLbriiERiRUFQTJ1FItIDCkIkm1fDXlDYMTEqCsREckYBUGy2tUw6kzI0o9FROJDR7wunZ3w/tvqHxCR2FEQdNnzLrTtVf+AiMSOgqCLXlYvIjGlIOhSuxqycmDkR6KuREQko0INAjO7yswqzWyTmd16lO0+a2ZuZuVh1tOThSuqufBfXmDJS8+z0UtZuGZnpksQEYlUaEFgZtnA/cBMoAyYa2ZlPWw3FPgm8HpYtRzJwhXV3PbUGqrrWyjL2sqq9nHc9tQaFq6oznQpIiKRCfOM4Dxgk7tvdvc24DFgdg/b3QX8GGgNsZYe3VtRScuBDkqoZ6TVs7ZzPC0HOri3ojLTpYiIRCbMIBgDbEuar0osO8jMzgbGuvufjvaNzGyemS03s+V1dXVpK7CmvgWAM7K2ALCuc8Jhy0VE4iCyzmIzywLuA75zrG3dfYG7l7t7eUlJSdpqGF1UAECZbQVgnY8/bLmISByEGQTVwNik+dLEsi5DgTOBJWa2BTgfWJTJDuP5M6ZQkJtNWdYWtnaOpIlCCnKzmT9jSqZKEBGJXE6I33sZMNnMJhIEwHXAf+9a6e4NQHHXvJktAf6nuy8PsabDzJketFSdteg93m4fz5iiAubPmHJwuYhIHIQWBO7ebma3ABVANvCwu681szuB5e6+KKzP7o05ZcPg6VrGXnEjMy+5LOpyREQyLswzAtx9MbC427I7jrDtpWHWckR7tgTj4smRfLyISNR0Z3FjbTAepuYgEYknBUFjov962KnR1iEiEhEFQVMtWBYMGRV1JSIikVAQNNbA4JGQnRt1JSIikVAQNNaoWUhEYk1B0FSrjmIRiTUFQWM1DNUZgYjEV7yDoG0ftDbAsNFRVyIiEpl4B8HBewgUBCISX/EOgqaaYKwgEJEYi3cQNCaCYKiCQETiS0EAunxURGJNQZA/HPIGR12JiEhk4h0EuodARCTmQdBYo3sIRCT2FATqHxCRmItvEHQcgL3vq2lIRGIvvkGw933A1TQkIrEX3yDQm8lERIBYB4HeTCYiAnEOgiadEYiIQJyDoLEasgdBwYioKxERiVSoQWBmV5lZpZltMrNbe1j/bTNbZ2arzex5MxsfZj2HaawNHjZnlrGPFBHpi0ILAjPLBu4HZgJlwFwzK+u22Qqg3N2nAU8CPwmrng9oqtVTR0VECPeM4Dxgk7tvdvc24DFgdvIG7v6iuzcnZpcCpSHWc7jGagWBiAjhBsEYYFvSfFVi2ZHcCDwTYj2HuAdNQ7qHQESEnKgLADCzLwLlwCVHWD8PmAcwbty4E//A5t3QsV9nBCIihHtGUA2MTZovTSw7jJldAXwPmOXu+3v6Ru6+wN3L3b28pKTkxCvTm8lERA4KMwiWAZPNbKKZ5QHXAYuSNzCz6cAvCUJgR4i1HE5vJhMROSi0IHD3duAWoAJYDzzu7mvN7E4zm5XY7F5gCPCEma00s0VH+Hbp1agzAhGRLqH2Ebj7YmBxt2V3JE1fEebnH1FjDVgWDBkVyceLiPQl8byzuKkmCIHsPtFXLiISqXgGgS4dFRE5KKZBUKP+ARGRhHgGQZOCQESkS/yCoG0ftDYoCEREEuIXBF1vJtM9BCIiQCyDQG8mExFJFr8g0JvJREQOE78g6Doj0OWjIiJALIOgFvKLIK8w6kpERPqE+AWB3kwmInKY+AWB3kwmInKYGAaBHi8hIpIsXkHQcQD2vq8rhkREksQrCPa+D7juIRARSRKvIDj4QhqdEYiIdIlnEKiPQETkoHgGga4aEhE5KF5B0FQDOflQMCLqSkRE+ox4BUHXpaNmUVciItJnxCwIatRRLCLSTbyCoKlGl46KiHQTnyBwD5qG1FEsInKYUIPAzK4ys0oz22Rmt/awfpCZ/T6x/nUzmxBaMc27oWO/3kwmItJNaEFgZtnA/cBMoAyYa2Zl3Ta7Edjj7h8Gfgr8OKx6Dr2ZTEEgIpIszDOC84BN7r7Z3duAx4DZ3baZDfw6Mf0kcLlZSJf0HHwzmYJARCRZmEEwBtiWNF+VWNbjNu7eDjQAJ3f/RmY2z8yWm9nyurq646smfzicfi0UjTu+rxcRGaByoi4gFe6+AFgAUF5e7sf1TcadHwwiInKYMM8IqoGxSfOliWU9bmNmOcBwYFeINYmISDdhBsEyYLKZTTSzPOA6YFG3bRYBX0pMfw54wd2P7y9+ERE5LqE1Dbl7u5ndAlQA2cDD7r7WzO4Elrv7IuAh4D/MbBOwmyAsREQkg0LtI3D3xcDibsvuSJpuBT4fZg0iInJ08bmzWEREeqQgEBGJOQWBiEjMKQhERGLO+tvVmmZWB2w9zi8vBnamsZz+QPscD9rneDiRfR7v7iU9reh3QXAizGy5u5dHXUcmaZ/jQfscD2Hts5qGRERiTkEgIhJzcQuCBVEXEAHtczxon+MhlH2OVR+BiIh8UNzOCEREpBsFgYhIzA3IIDCzq8ys0sw2mdmtPawfZGa/T6x/3cwmZL7K9Ephny82s7fMrN3MPhdFjemWwj5/28zWmdlqM3vezMZHUWc6pbDP/2hma8xspZm90sN7wvudY+1z0nafNTM3s359SWkKv+MbzKwu8TteaWY3nfCHuvuAGggeef0OMAnIA1YBZd22+TrwQGL6OuD3UdedgX2eAEwDHgE+F3XNGdrnTwCFiemvxeT3PCxpehbwbNR1h73Pie2GAn8FlgLlUdcd8u/4BuBn6fzcgXhGcB6wyd03u3sb8Bgwu9s2s4FfJ6afBC43M8tgjel2zH129y3uvhrojKLAEKSyzy+6e3NidinBW/L6s1T2uTFpdjDQ368GSeX/M8BdwI+B1kwWF4JU9zetBmIQjAG2Jc1XJZb1uI27twMNwMkZqS4cqezzQNPbfb4ReCbUisKX0j6b2c1m9g7wE+CfMlRbWI65z2Z2NjDW3f+UycJCkuq/688mmjyfNLOxPazvlYEYBCKHMbMvAuXAvVHXkgnufr+7fwj4LnB71PWEycyygPuA70RdSwb9EZjg7tOAP3OodeO4DcQgqAaSE7I0sazHbcwsBxgO7MpIdeFIZZ8HmpT22cyuAL4HzHL3/RmqLSy9/T0/BswJtaLwHWufhwJnAkvMbAtwPrCoH3cYH/N37O67kv4tPwicc6IfOhCDYBkw2cwmmlkeQWfwom7bLAK+lJj+HPCCJ3ph+qlU9nmgOeY+m9l04JcEIbAjghrTLZV9npw0ew2wMYP1heGo++zuDe5e7O4T3H0CQV/QLHdfHk25JyyV3/GpSbOzgPUn/KlR95KH1PN+NbCBoPf9e4lldxL8AwHIB54ANgFvAJOirjkD+3wuQXvjPoKzn7VR15yBff4L8D6wMjEsirrmDOzz/wLWJvb3ReCMqGsOe5+7bbuEfnzVUIq/43sSv+NVid/x6Sf6mXrEhIhIzA3EpiEREekFBYGISMwpCEREYk5BICIScwoCEZGYUxCIiMScgkBEJOYUBCInyMzOTTwALN/MBpvZWjM7M+q6RFKlG8pE0sDMfkRwx3oBUOXu90RckkjKFAQiaZB4Lswygufhf8zdOyIuSSRlahoSSY+TgSEET8PMj7gWkV7RGYFIGpjZIoLHPk8ETnX3WyIuSSRlOVEXINLfmdn1wAF3/52ZZQOvmtll7v5C1LWJpEJnBCIiMac+AhGRmFMQiIjEnIJARCTmFAQiIjGnIBARiTkFgYhIzCkIRERi7v8DnwEJQ/BxzncAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from lmfit import Model\n",
    "\n",
    "gmodel = Model(func, independent_vars=['x'],param_names=['a','b'])\n",
    "params = gmodel.make_params(a=2., b=1.0)\n",
    "result = gmodel.fit(y, params, x=x)\n",
    "\n",
    "print(result.fit_report())\n",
    "xrange = np.linspace(0,0.5)\n",
    "\n",
    "fitted_y = result.eval(x=xrange)\n",
    "\n",
    "\n",
    "plt.plot(x,y,'o')\n",
    "plt.plot(xrange,fitted_y)\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_SZfcjd02TnU"
   },
   "source": [
    "-   Here the two intervals are relatively small, and do not include zero, suggesting both parameters are significant.\n",
    "\n",
    "-   More importantly, the errors are not skewed by a nonlinear transformation.\n",
    "\n",
    "-   Note you have to provide an initial guess.\n",
    "    -   This will not always be easy to guess.\n",
    "    -   There may be more than one minimum in the fit also, so different guesses may give different parameters.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "itc32Ju085QS"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "26-bonus-nonlinear-regression.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}