{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Least squares optimization\n",
    "\n",
    "Many optimization problems involve minimization of a sum of squared residuals. We will take a look at finding the derivatives for least squares minimization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In least squares problems, we usually have $m$ labeled observations $(x_i, y_i)$. We have a model that will predict $y_i$ given $x_i$ for some parameters $\\beta$, $f(x) = X\\beta$. We want to minimize the sum (or average) of squared residuals $r(x_i) = y_i - f(x_i)$. For example, the objective function is usually taken to be \n",
    "\n",
    "$$\n",
    "\\frac{1}{2} \\sum{r(x_i)^2}\n",
    "$$\n",
    "\n",
    "As a concrete example, suppose we want to fit a quadratic function to some observed data. We have\n",
    "\n",
    "$$\n",
    "f(x) = \\beta_0 + \\beta_1 x + \\beta_2 x^2\n",
    "$$\n",
    "\n",
    "We want to minimize the objective function\n",
    "\n",
    "$$\n",
    "L = \\frac{1}{2} \\sum_{i=1}^m (y_i - f(x_i))^2\n",
    "$$\n",
    "\n",
    "Taking derivatives with respect to $\\beta$, we get\n",
    "\n",
    "$$\n",
    "\\frac{dL}{d\\beta} = \n",
    "\\begin{bmatrix}\n",
    "\\sum_{i=1}^m f(x_i) - y_i \\\\\n",
    "\\sum_{i=1}^m x_i f(x_i) - y_i \\\\\n",
    "\\sum_{i=1}^m x_i^2 f(x_i) - y_i\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with matrices\n",
    "\n",
    "Writing the above system as a matrix, we have $f(x) = X\\beta$, with\n",
    "\n",
    "$$ \n",
    "X = \\begin{bmatrix}\n",
    "1 & x_1 & x_1^2 \\\\\n",
    "1 & x_2 & x_2^2 \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "1 & x_m & x_m^2 \n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "and \n",
    "\n",
    "$$\n",
    "\\beta = \\begin{bmatrix}\n",
    "\\beta_0 \\\\\n",
    "\\beta_1 \\\\\n",
    "\\beta_2\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "We want to find the derivative of $\\Vert y - X\\beta \\Vert^2$, so\n",
    "\n",
    "$$\n",
    "\\Vert y - X\\beta \\Vert^2 \\\\\n",
    "= (y - X\\beta)^T(y - X\\beta) \\\\\n",
    "= (y^T - \\beta^TX^T)(y - X\\beta) \\\\\n",
    "= y^Ty - \\beta^TX^Ty -y^TX\\beta + \\beta^TX^TX\\beta\n",
    "$$\n",
    "\n",
    "Taking derivatives with respect to $\\beta^T$ (we do this because the gradient is traditionally a row vector, and we want it as a column vector here), we get (after multiplying by $1/2$ for the residue function)\n",
    "\n",
    "$$\n",
    "\\frac{dL}{d\\beta^T} =  X^TX\\beta - X^Ty\n",
    "$$\n",
    "\n",
    "For example, if we are doing gradient descent, the update equation is\n",
    "\n",
    "$$\n",
    "\\beta_{k+1} = \\beta_k + \\alpha (X^TX\\beta - X^Ty)\n",
    "$$\n",
    "\n",
    "Note that if we set the derivative to zero and solve, we get\n",
    "\n",
    "$$\n",
    "X^TX\\beta - X^Ty = 0\n",
    "$$\n",
    "\n",
    "and the normal equations\n",
    "\n",
    "$$\n",
    "\\beta = (X^TX)^{-1}X^Ty\n",
    "$$\n",
    "\n",
    "For large $X$, solving the normal equations can be more expensive than simpler gradient descent. Note that the Levenberg-Marquadt algorithm is often used to optimize least squares problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "You are given the following set of data to fit a quadratic polynomial to:\n",
    "\n",
    "```python\n",
    "x = np.arange(10)\n",
    "y = np.array([  1.58873597,   7.55101533,  10.71372171,   7.90123225,\n",
    "                -2.05877605, -12.40257359, -28.64568712, -46.39822281,\n",
    "                -68.15488905, -97.16032044])\n",
    "```\n",
    "\n",
    "Find the least squares solution using gradient descent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.arange(10)\n",
    "y = np.array([  1.58873597,   7.55101533,  10.71372171,   7.90123225,\n",
    "                -2.05877605, -12.40257359, -28.64568712, -46.39822281,\n",
    "                -68.15488905, -97.16032044])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def f(x, y, b):\n",
    "    return (b[0] + b[1]*x + b[2]*x**2 - y)\n",
    "\n",
    "def res(x, y, b):\n",
    "    return sum(f(x,y, b)*f(x, y, b))\n",
    "\n",
    "# Elementary form of gradient\n",
    "def grad(x, y, b):\n",
    "    n = len(x)\n",
    "    return np.array([\n",
    "            sum(f(x, y, b)),\n",
    "            sum(x*f(x, y, b)),\n",
    "            sum(x**2*f(x, y, b))\n",
    "    ])\n",
    "\n",
    "# Matrix form of gradient\n",
    "def grad_m(X, y, b):\n",
    "    return X.T@X@b- X.T@y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  227.0657638 ,  1933.9094954 , 15758.14427298])"
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad(x, y, np.zeros(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  227.0657638 ,  1933.9094954 , 15758.14427298])"
      ]
     },
     "execution_count": 140,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.c_[np.ones(len(x)), x, x**2]\n",
    "grad_m(X, y, np.zeros(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.55079998,  7.31478229, -2.04118936])"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.linalg import solve\n",
    "\n",
    "beta1 = solve(X.T@X, X.T@y)\n",
    "beta1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "max_iter = 10000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.73391723,  7.23152392, -2.03359658])"
      ]
     },
     "execution_count": 144,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = 0.0001 # learning rate\n",
    "beta2 = np.zeros(3)\n",
    "for i in range(max_iter):\n",
    "    beta2 -= a * grad(x, y, beta2)\n",
    "beta2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.73391723,  7.23152392, -2.03359658])"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = 0.0001 # learning rate\n",
    "beta3 = np.zeros(3)\n",
    "for i in range(max_iter):\n",
    "    beta3 -= a * grad_m(X, y, beta3)\n",
    "beta3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtEAAAEICAYAAACZEKh9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4lFX6xvHvkwKEGhGUFhZEWUU0qKEqoCsiqD+xoRSx\nh8WKrgqia+9lratiw65YULGAAqsC6oKCioplpUkJKCqht5Dz++NMJERISD3vJPfnuuYymZnMPEFu\n3mfe9xRzziEiIiIiIjsvIXQBIiIiIiLxRk20iIiIiEgxqYkWERERESkmNdEiIiIiIsWkJlpERERE\npJjURIuIiIiIFJOaaCkzZtbCzJyZJYWuRSSqzOw6M3sudB0iEp/MbLaZHRq6DgE1OyIigpk5YC/n\n3JzQtYhURWb2FLDYOffPwp7nnNu3YiqSouhMtIiIlIp5Op6IlCNd5Y0e/aMnfzCz4Wa2xMxWm9kP\nZjbQzNabWf18zznAzH41s2QzSzSzu2LfzwOODli+SKSYWRMzG2Nmy81svpldtIPndTKzT8ws28xm\n5b9Ma2YfmtlNscfXmNlbZrarmT1vZqvM7DMza5Hv+Xub2UQz+z2W4ZPzPfaUmT1oZu/EMj7dzFrF\nHpsSe9qs2PucYma7mNnbsfpXxL5uVqC2m83sY2AdcKmZzSzwu/3DzMaW/k9TJLrMbIGZXW5mX5nZ\nWjN7wsx2N7PxsaxNMrNdYs99xcyWmdlKM5tiZvvG7h8MDASG5WU932sPN7OvgLVmlhS7r0fs8XFm\n9q98tYw2s1EV/odQRamJFgDM7K/ABUB751wd4EhgGvBf4MR8Tx0AvOqc2wxkAscABwAZwEkVWrRI\nRMXOyr4FzAKaAocDF5vZkQWe1xR4B7gJqA9cBowxs4b5ntYPGBR7nVb4TD4Ze/53wLWx16oFTARe\nAHaL/dxDZtamwGtdD+wCzAFuBnDOdYs9nu6cq+2cewl/fHgS+AvQHFgP/LvArzoIGAzUAe4HWprZ\nPgUef6bIPzCR+HcicATQGvg/YDxwJdAQn6W8D9Hjgb3wGf0ceB7AOfdo7Os7Yhn8v3yv3R9/kirV\nOZdT4H3PAgaZ2d/MbCDQARha9r+ebI+aaMmzBagOtDGzZOfcAufcXPwBuT/4S7b4g/ALsZ85GbjX\nObfIOfc7cGuAukWiqD3Q0Dl3g3Nuk3NuHvAYPj/5nQqMc86Nc87lOucmAjOAo/I950nn3Fzn3Er8\nAXiuc25S7GD6Cv5DLPgPtAucc08653Kcc18AY4C++V7rdefcp7GffR5ot6NfwDn3m3NujHNunXNu\nNb7h7l7gaU8552bH3m8j8FLsdyJ2hq0F8PZO/HmJxLsHnHM/O+eWAFOB6c65L5xzG4DXieXUOTfK\nObc6lpfrgHQzq1fEa98fO86uL/iAc24ZcC7wNHAfcFosr1IB1EQLALHJRBfjQ/1L7JJQE/xBuLOZ\nNQa6Abn4fyAAmgCL8r3MTxVXsUik/QVoEhuikW1m2fizUrtv53l9CzzvEKBxvuf8nO/r9dv5vna+\n1+pY4LUGAo3yPX9Zvq/X5fvZPzGzmmb2iJn9ZGargClAqpkl5nvaogI/9jQwIPaBexDwcqxZEKns\nisxpbAjkbWY2N5apBbHHGxTx2gVzVtBbQCLwg3Puo2LULKWkJlr+4Jx7wTl3CP5g7IDbnXMrgAnA\nKfihHKOdcy72I0uBtHwv0bwi6xWJsEXAfOdcar5bHefcUdt53rMFnlfLOXdbCd9zcoHXqu2cO7eE\nv8OlwF+Bjs65uvgP0QCW7zku/w8456YBm4Cu+H8vni3he4tURgOAPkAPoB7+Sg1szZTbzs8Udn+e\nm/FDuxqbWf9S1ijFoCZaAD8mOjamqjqwAf/JOTf28AvAafgxzy/k+7GXgYvMrFls0sQVFVmzSIR9\nCqyOTQhKiZ2Bamtm7Qs87zng/8zsyNhzapjZofkn8BXD20BrMxtkfuJvspm1LzBGuTA/A3vk+74O\n/t+BbPOTi6/dydd5Bj92erPOiolsow6wEfgNqAncUuDxghkskpl1A87EH6NPBx6IzbWQCqAmWvJU\nB24DfsVf8t0NGBF77E38RIhlzrlZ+X7mMeA9/OSpz4HXKqxakQhzzm3Bj1FuB8zH5+px/Nmn/M9b\nhD8zdSWwHH82+XJK8G9zbBxkT/y46yx8jm/HZ3tnXAc8HRsKcjJwL5ASq30a8O5Ovs6zQFv8BwQR\n2eoZ/LDHJcC3+Fzl9wR+XlK2mb1R1IuZWd3Ya17gnFvinJsae40nY0OqpJzZ1ivzIiIipWNmKcAv\nwIHOuR9D1yMiUl50JlpERMrSucBnaqBFpLLT7jciIlImzGwBfpLUcYFLEREpdxrOISIiIiJSTBrO\nISIiIiJSTHExnKNBgwauRYsWocsQiYyZM2f+6pxrWPQzw1BmRbYV5cwqryLb2tm8xkUT3aJFC2bM\nmBG6DJHIMLNI7w6pzIpsK8qZVV5FtrWzedVwDhERERGRYlITLSIiIiJSTGqiRURERESKSU20iIiI\niEgxqYkWERERESmmuFidI55lZa9n5OS5zFqUTXpaKkO6t6JJakroskRkB5RZkfihvEpIaqLLUVb2\nenrfN5W1G3PIyXXMzlrF2C+zGD+0q0IuEkHKrEj8UF4lNA3nKEcjJ8/9I9wAObmOdRtzGDl5buDK\nRGR7lFmR+KG8SmhqosvRrEXZf4Q7z+Zcx6xF2YEqEpHCKLMi8UN5ldDURJej9LRUkhJsm/uSE4z0\ntNRAFYlIYZRZkfihvEpoGhNdXpzj3PT6zH7vY+r99jPVNm3gq7Q2rN11N4Z0b1Xil9UkCpFysGED\nZGVxUfIyNv7vI2qt/J0fdk3jy+b7kpRSU5kViZItW+CXX7iwbjYr539KvRXLya5Wi+kt27GhfgPl\nVSqMmujtKDJEzsHvv8PixbBokf9v/lvsvsbr1jGmwGtv3mdfkrOPgCOOgG7doHbtYtWlSRQi2yo0\nr7m5sHw5LFkCWVn+v3m3/N///jsADYDb87325mrV2XJIV2rU/tpndr/9wOxPNRRWmzIrslWRx9dV\nq7af0fzfL1sGW7bQELivwOtv3m9/ktce6fN6yCGQsvM5U16luNREF5CVvZ7e906h6eK5NFmxDFvz\nK+Pu+p0BTROp+XPW1kZ5w4ZtfzAxEZo0gWbNID0djjnGf513S0yEDz8kedIkGDkS7r0XkpKgc2cf\n9h49oH17f98OFDaJ4oY+bcvzj0UkkrKy13PUPZM58PtP6fzzfBqt/Z3ZI36nQe1NVFu2FJYuhZyc\nbX8oIQF23x2aNoU99oCuXX12mzb1tyZNYNddYeZMkidOJHnCBLjsMv+zjRr5rPbs6f/buHGh9Smz\nIlvlNam7Z82ny9zPabz2Nz6/Jptd6+VQ/eelvkles+bPP7jLLlszuu++W3OaP7NLlkBeXu+9F+68\nE2rU8Pnu2dMfZ/ffv9APwcqrFJea6Pyc44M7n+DpJx+g3dL//XH35oRE1uy6OzVbt4SDDoLjjtu2\nQU5L8wflxMTCX79DBxg2zDfgH38MkybBxIlw7bVwzTVQty4cdpg/OB9xBLRuvU3gNYlCJJ8tW/j4\nxgd46dmH+OvynwBYVa0mv9TZlUXVmtLqsMO2Pcjmfb377oV+WP3D0Uf7G/gPzhMn+tt778Fzz/n7\n99vPZ7VnT3+wrllzm5dQZkW2euPJt7n9xfvo+cMnJODYmJjEL7V3ZVmTJvwlPR2OOurPmW3S5E+5\n2q7GjSEjA0aM8I34lCk+rxMmwOWX++fsvvvWD8FHHPGnD8HKqxSXmmjwl3xffx1uuomBX37Jonq7\nc/URQ5jVuDVL6zTk11r12D9tF8ZecEjZvF+NGnD44f52663w22/w/vu+qZ40CcaO9c9r1mxrQ334\n4aSnpTI7a9U2IdckCqlyNm2CZ5+F226j75w5/G/X5gw95lIm7tWJddX8Jdf0ZvXKLq/gs3jmmf6W\nmwuzZm09QD/4INx9N1Sr5hvpvKY6PV2ZFQF/0ujmmzlv/HhWVa/FvzufzAvterOszq5gVvZ5rV3b\nN+RHHeW/j52l/uP2/PP+/rZtt+a1WzflVYqtajfRW7bASy/BzTfDt9/CXnvx2oU3cmXK/mywrWeV\nyz1Eu+4Kffv6G8C8eVsb6jffhKeeAuDqNm3Zs15rPkxL579N9yWnRgo1qyeVahKFSNxYtw4ef9xf\npl28GA46iBevuJdraMVmt/WKTbnnNSEBDjjA34YN83VNnbq1qb7iCn9r2JAruh0GOc34T1o6S2rt\nSnKCKbNSNTjnj2E33wyTJ0PDhkw8dSiXNzyY7GpbzyxXSJPatCmccYa/5ebCV19tzetDD8E990C1\nalzVqQu7Ju/B+83T+aphS5ISE5VXKZQ554p+VmAZGRluxowZZfeCmzf7T6K33AI//ujHWF11FZx8\nMlmrN20zsSDvoBdsYsGWLfDllz7wkybhPvoI27iRVTXrMjFzOJ2vvZgmu+zEpS6pVMxspnMuI3Qd\nO1KmmV21auuB7pdf/Nneq66Cnj3JWrkhWnkFPw47b6jWxImwbBm5Zow/5Hi+OG8YZ/XaX5OUqqAo\nZ7ZM85qbC2+95Zvnzz7zDezll0NmJlmbLHp5Xb8ePvrIN9QTJ/qrTMAPLdow4R83c+KgI5XXKmin\n8+qci/ztoIMOcmViwwbnRo50rkUL58C5du2cGzPGuS1btnnakhXr3NVvfO2OfWCqu/qNr92SFevK\n5v3Lwrp1zr33nnNduvjf4cgjnZs/P3RVUsGAGS4C2dzRrUwy++uvzl19tXOpqf7veq9ezk2Z8qen\nRTqvubnOffWVcxde6JyZc2lpzo0bF7oqCSDKmS2TvG7e7NzzzzvXtq3P6x57OPfoo/64m0+k8+qc\nc8uW+T6hQQPnkpOdu+465zZuDF2VVLCdzWvw8O7MrdQBX7fOufvuc65pU/8rd+jg3Ftv+QNcvNqy\nxbkHHnCudm3natVy7t57ncvJCV2VVJAoH5BdaTObleXcpZf6v9fg3AknODdjRslfLyo++cS5ffbx\nv9Oppzq3fHnoiqQCRTmzpcrrxo3OPfaYc61a+b/bbdo499xzvqmOZ7/84tyAAf53atvWuWnTQlck\nFWhn81q5dyxcswbuugtatoShQ/1yVhMmwLRpfgm6Yqz3GjkJCXDBBTB7tl9v+uKL/ZqYs2eHrkyk\nZBYsgPPO83m99144/nj/93nMGL8qTrzr3Bm++MKvxDN6NLRp4//roj+kTuRP1q2D+++HVq0gMxNS\nU+G11+Drr2HgwJ1bASfKGjb0wz7ffhuys31+//EPWLs2dGUSIZWziV650o/HatHCj8Xabz8/sWHK\nFD8TN56b54KaN4d33vFhnzPHT3a67jrYuDF0ZSI75/vv/YSfPfeEJ56A00+H//3Pr8DRpk3o6spW\n9epw/fXw+ef+36f+/aFPHz9RUiQerFoFt93m//4OHeo/9L77rh//fPzx/gRPZXL00f7D/Lnn+nkZ\n++3n5zyIUNma6N9/92sut2gB//wndOoEn3ziJwt06xa6uvJjBgMG+BVGTj7ZH6QPPNCfcReJqi+/\n9CvStGkDL78MF17oV6Z55BF/1agy228/+O9/4V//8gfkfff1v3dubujKRLbvt9/8VZS//MWvxXzQ\nQf7E1JQpcOSRlevkVEF16/qlLCdPhuRkfzLu7LNhxYrQlUlglaKJXrogi6knns36Js3ghhtYf0g3\nmDnTX4bp3Dl0eRWnYUO/CcQ778Dq1dCliz9TsL0doEQCycpez1fdjoIDDmDDuHdZfcll8NNP/ixP\n06ahy6s4iYn+8vDXX/tNIoYMgb/9za8YJBIRWdnrGXPxLaxvkgY33sj6rt1hxgwYP96vlFOVdOvm\nV+8YMQKeftqfAHjttdBVSUBx30RnZa/nlIc/od07o5m4R3t6n/MgnTLOJ2uPfUKXFs5RR/nLT+ef\nDw884BeUf++90FWJ/LHt77uJu3Nn10F0/vsTHFK3B1nJtUOXFk6rVv5s9OOP+7Pz++8Pd9zx5+3K\nRSpYXl6fXVmL8Xt1olfmw3Q66DyyWlWyYVbFUaOGXx73s8/8jocnnuhvS5eGrkwCiPsmeuTkuWQl\n1uTgIaO46NhhfLfrX/7Y675Kq1PHN9AffQQpKdCrF5x2Gvz6a+jKpAobOXkuazfm8FDHk3iwyyms\nqFZLeQV/Kfzss/2QrF69YPhw6NjRN9UigeTl9cvd9+Qfx1zK9/XTlNc8BxwA06f78eHvvOPPSj/5\npCYKVzHBmmgz62VmP5jZHDO7oqSvk7fX/aoaW89kaa/7fLp08Qfiq6+GF1/0QX/xRQVdiqWs85qf\n8ppPkyb+8vArr/itijMy/MYyGzaErkziiPJaQZKT/Qfer77yV5DOOstvIT5vXujKpIIEaaLNLBF4\nEOgNtAH6m1mJrg+lp6WSlLDthAbtdV9A9epwww1+RYCWLf0kxGOPhUWLQlcmcUB5rWBmcNJJ/qz0\noEH+0nG7dv6qkkgRlNcAWreGDz6Ahx/2Z6f328/P8diyJXRlUs5CnYnuAMxxzs1zzm0CRgN9SvJC\nQ7q3olb1pD+CnreNqPa634799vOrldx9N7z/vl8R4KGHtCKAFEV5DaF+fX95+L33/JKVXbv6eQ6r\nVoWuTKJNeQ0hIcFPDp49Gw47zE8aPvhg+Oab0JVJOQrVRDcF8p8GXRy77w9mNtjMZpjZjOXLl+/w\nhZqkpjB+aFcGdGxOerN69O/YnPFDu2qv+x1JTIRLLvHB7tjRH5S7d/fr8opsX5F5hZ3LrPJaAj17\n+hU8hg71Z7ratvUrI4hsn/IaUloavPUWvPACzJ3rl5u97jrYtCl0ZVIOzAUYG2tmJwG9nHPnxL4f\nBHR0zl2wvednZGS4GTNmVGSJVYNzrHj4MaoPv5ycXHj2llEcf/pR+gcyDpjZTOdcRgW9V7HyCsps\neVk+4UPcOeew26K5vHH+dXS4aZjyGicqKrPKa4T8+ivrzruQmq+MZvp+h/DuNfeT2WNvZTYO7Gxe\nQ52JXgKk5fu+Wew+qUBZKzdw6PIWHH3q3axNSKb/iLO44KpnycpeH7o0iRblNQKystfT47+bOLTf\nXXzY8iCOffB6HjrzGuVVClJeIyIrqRad9z2Ta3ueS8evP6LjPy/g/+7+QJmtREI10Z8Be5lZSzOr\nBvQD3gxUS5WVt3zR/HqN6N//FjYnJvHI08N59fmJoUuTaFFeIyAvr+sSkvn78Vfy8V/SueGNf/Hx\nDfeHLk2iRXmNiLzMPn3A0Vx/eCa9fviE61+9nUff/yF0aVJGgjTRzrkc4ALgPeA74GXn3OwQtVRl\n+ZcvWlC/KQP63QzAqVeepV3T5A/KazTkz+vG5OpknvhPpjdvywn3XumXxBNBeY2S/Jl9MqMPNx96\nFsd8N4Xutw3Xyh2VRLB1op1z45xzrZ1zrZxzN4eqoyoruHzR3F3TOK3/LVQn128/rLUuJUZ5Da9g\nXjck12BI32tZvHc69O8Pr78esDqJEuU1Ggpm9rGOJ/Cv7qdx2GcTIDNTK2NVAnG/Y6GU3PaWL8pq\n1oq1b4+Hdev8Mj0//RS4ShGB7eeVOnWo9u54aN8eTjkF3n47cJUikmd7mX3msIGsHn6lX77y3HO1\n8VmcSwpdgISTt3zRyMlzmbUom/S0VIZ0b8VuqSkwcSIcfrhvpCdP9sv2iEgwO8pr49QUv+TdEUfA\niSfC2LF+63ARCWpHma1T7whIcHDrrX7Xwwce8JssSdxRE13FNUlN4YY+bf/8wIEHwoQJ0KOHH9ox\nebLfklhEgtlhXlNT/aYshx8Oxx3nz0j36FHxBYrINnaY2Ztvhs2b4a67fCN9991qpOOQhnPIjrVv\n7w/My5b5RnrZstAViciO1K/vryC1bg3HHus/+IpINJnBHXfARRfBvffCFVdoaEccUhMthevUyV8q\nXrzYn+X65ZfQFYnIjjRoAJMmQYsWcPTR8NFHoSsSkR0x8w30uef6hvqaa0JXJMWkJlqKdsgh8M47\nMH++v0T822+hKxKRHdltN/jPf6BpUzjqKJg2LXRFIrIjZvDvf8M558BNN8GNN4auSIpBTbTsnO7d\n4a23/PrRRxwBK1aErkhEdqRxY3j/fd9Q9+oF2tJZJLoSEuCRR+D00/3Z6NtuC12R7CQ10bLzDj8c\n3ngDZs+Gnj0hOzt0RSKyI02b+kZ6l118Xr/4InRFIrIjCQnwxBMwYACMGOEnGkrkqYmW4jnySHjt\nNZg1C3r3hlWrQlckIjvSvDl88AHUru2vIH39deiKRGRHEhPh6aehb1+49FK/9J1EmppoKb6jj/bb\nDM+Y4cdcrlkTuiIR2ZEWLXwjXb26v5r07behKxKRHUlKguef90tVXnQRjBwZuiIphJpoKZk+feDF\nF/2kpWOOgbVrQ1ckIjvSqpUf2pGY6BvpH34IXZGI7EhyMrz0kj+2nnuuH+YhkaQmWkrupJPguedg\n6lTfVK9fH7oiEdmRv/7Vr9qxZYtf933OnNAViciOVKsGr77qJwZnZsIzz4SuSLZDTbSUTr9+8NRT\n/izX8cfDhg2hKxKRHWnTxjfSGzf6RnrBgtAViciOVK/u5yD97W9w5pn+6q9EippoKb1Bg+Dxx/3u\nhied5A/QIhJN++3nN2RZswYOOwwWLgxdkYjsSEoKvPkmdO3qj7WvvBK6IslHTbSUjbPO8utcvvMO\n33XrzQn3fcg1Y78hK1tDPEQip107mDCB3N9X8FvHQzjzpteVV5GoqlkT3n6bjRkd2NJ/ALdk3qK8\nRoSaaCkzWScP4pajzmefTz+g76hbeWH6QnrfN1VBF4mgrD335bSTb6D6779y5b1DeX3q/5RXkYjK\nyknksMMu5+vdW3HZqGv5fsy7ymsEqImWMjNy8lxGpR/Fg5360v+rCfT47iPWbcxh5OS5oUsTkQJG\nTp7LtIZ7Mvj4q2j122JGTHpUeRWJqJGT5/KLVee0vteztG4D7h57JwkrVyqvgamJljIza1E2ObmO\new4ZyJeNW3Pbuw/QIPsXZi3SzoYiUZOX109atOORjicyYNZ7/O37j5VXkQjKy+uqGrW5+JjLaLT6\nV65970HlNTA10VJm0tNSSUowchKTGPp/l5GUu4V73rmbdk3qhC5NRArIyyvA3V0H8lWjPblt/AMc\nUlMTg0WiJn9ev2i6N/cd3J/jvp3MoLlTA1dWtamJljIzpHsralVPIinB+GmXJtx4xN/ptPBrLp31\nZujSRKSA/HndnJjMpccOo8aWTVz09E2Qmxu6PBHJJ39eAR49+BQ+b74vJzxxK8ybF7i6qktNtJSZ\nJqkpjB/alQEdm5PerB7VzzmL9cefSN1bboDPPgtdnojkUzCvnY/qzMZ/3UP1KR/CXXeFLk9E8imY\n11M6t6TpW2NISEiAU0+FnJzQJVZJ5pwLXUORMjIy3IwZM0KXISWxYoVfTqtaNfjiC6hdO3RFlYKZ\nzXTOZYSuY0eU2TjlHPTtC2PHwrRpcNBBoSuqNKKcWeU1jo0eDf37wzXXwPXXh66m0tjZvOpMtJSv\nXXbxW4PPmwcXXRS6GhEpjBk8+ig0agQDBsDataErEpHC9OsHp50GN90EH30UupoqR020lL+uXeHK\nK+HJJ7XbkkjU1a8Pzz4LP/4IF18cuhoRKcq//w0tWvhhHdlaraMiqYmWinHNNdCxIwwerG2GRaLu\n0ENh+HB4/HEYMyZ0NSJSmDp14IUXYPFiOO88PyxLKoSaaKkYycnw/PN+8sOgQbBlS+iKRKQw118P\nGRmQmekPziISXR07+sy++KIfQikVQk20VJxWreDBB2HKFLj99tDViEhhqlXzZ7c2bfJjLvXBVyTa\nrrgCunWD88+HudrJsCKUWxNtZnea2fdm9pWZvW5mqfkeG2Fmc8zsBzM7srxqkAgaNMhPhLj2Wvj0\n09DVSIzyKtu1115w//3wwQdw552hq5EY5VW2KzHRz2dITISBA2Hz5tAVVXrleSZ6ItDWObc/8D9g\nBICZtQH6AfsCvYCHzCyxHOuQKDGDhx+Gpk397P/Vq0NXJJ7yKtt35plw0klw9dVa7z06lFfZvubN\n4ZFHYPp0uOGG0NVUeuXWRDvnJjjn8lb/ngY0i33dBxjtnNvonJsPzAE6lFcdEkGpqX7M1vz5WvYu\nIpRX2aH8y94NHAhr1oSuqMpTXqVQJ58MZ5wBt9wCU7UteHmqqDHRZwHjY183BRble2xx7D6pSg45\nBK66Cp56Cl5+OXQ1si3lVbaVt977nDkwdGjoamRbyqv82f33Q8uWWvaunJWqiTazSWb2zXZuffI9\n5yogB3i+mK892MxmmNmM5cuXl6ZMiaprroFOnfyydz/9FLqaSq888xr7WWW2Muve3U9cGjUKXn01\ndDWVnvIqpZK37F1WFgwZomXvyklSaX7YOdejsMfN7AzgGOBwt3V/8SVAWr6nNYvdV/C1HwUeBb8l\naWnqlIhKSvLL3rVr5yccfvCBnxAh5aI88xp7fWW2srv+epg0yS9717EjpKUV/TNSIsqrlFqHDj6z\nV10FvXvD6aeHrqjSKc/VOXoBw4BjnXPr8j30JtDPzKqbWUtgL0DLNFRVe+zhl72bOhVuvTV0NVWW\n8io7JTnZn93avFnrvQekvMpOGz7cX0W64AI/HEvKVHmOif43UAeYaGZfmtlIAOfcbOBl4FvgXeB8\n55z+Ja7KTj0V+veH666DadNCV1NVKa+yc/bc028zPHky3HFH6GqqKuVVdk7esndJSVr2rhyYi4Nx\nMhkZGW7GjBmhy5DytHIlpKf7wH/xBdStG7qiSDOzmc65jNB17IgyW8k559d7f+01+OQTaN8+dEWR\nF+XMKq9VwCuv+FU7rroKbropdDWRt7N51Y6FEg316vnx0QsWwIUXhq5GRApjBiNHQuPGfr13LXsn\nEm19+/o132+5xe8aLGVCTbREx8EHwz//Cc88A6NHh65GRAqTt+zdvHla710kHtx/P7Rq5YdQrlgR\nuppKQU20RMvVV0Pnzn5JHi17JxJt3brBiBHw5JP+crGIRFft2n5i8NKlWvaujKiJlmjJW/bOOf9p\nOSen6J9HjWQyAAAgAElEQVQRkXCuvdYvdzd4MCxcGLoaESlM+/Zw441+k7Onnw5dTdxTEy3R07Il\nPPQQfPSRlr0TibrkZP/BNydHy96JxIPLL4dDD9Wyd2VATbRE08CB/nb99fDf/4auRkQK06qVX/Zu\nyhS4/fbQ1YhIYfKWvatWzU8M1rJ3JaYmWqLrwQf9jmgDB8KqVaGrEZHCnHYanHKKH97xqfb3EIm0\nZs3gscfgs8/8Hg1SImqiJbrq1WP5yFHk/vQTb/Y5h2vGfkNW9vrQVYnI9pix9PZ7yN6lIYv6nMx1\nr32pvIpEWNbhRzGzxwnk3norD93zqvJaAmqiJbKystfTY9pmRqcfSe8pr/HJuP/S+76pCrpIBGVl\nr6fXU19zRbezSVv2E+6RR5VXkYjKyl5P7/umcs5+p5Bdow4H3Hcjve+dorwWk5poiayRk+eydmMO\ndx88gA1J1Rj2wSjWbcxh5OS5oUsTkQLy8vpuq45MS2vLRVOfJ3FltvIqEkF5eV1RrRb3HDKAzj99\nRafZHyuvxaQmWiJr1qJscnIdv9bahYc6n0zPH6eRMX8WsxZlhy5NRArIyytm3Pi3c9hl/WoGf/yS\n8ioSQX/kFXgxvRdz6jdj+H+eYPb85YEriy9qoiWy0tNSSUowAEZl9GFx3d24+oPHadekTuDKRKSg\n/Hmd3WhPxrQ9nLNmvsmhyasDVyYiBeXPa05iEjf97Wz2WJHF2d+8F7iy+KImWiJrSPdW1KqeRFKC\nsTGpGncddgZtfp7HxUunhS5NRArIn1eA+w4dxJaERIaMfyxwZSJSUMG8frxne/67x4H0GvMo/P57\n4Orih5poiawmqSmMH9qVAR2bk96sHnXPOJVN7Tuwy83Xw5o1ocsTkXwK5vVvPQ5k86WXkTL2db9x\nkohERsG89u/0F1o9M5KEVSv9joayU8zFwd7pGRkZbsaMGaHLkCiYNg06d4ZrrvEbsVRRZjbTOZcR\nuo4dUWYFgLVr4a9/hSZNfHYTqu55myhnVnmVP/z97zBqFMyeDa1bh64mmJ3Na9X9F03iU6dO0K8f\n3HknLF4cuhoRKUytWnDLLX5DhxdeCF2NiBTlhhugRg0YNix0JXFBTbTEn1tvhdxcuOqq0JWISFFO\nPRUOPBBGjIB160JXIyKF2X13uPJKGDsWPvggdDWRpyZa4k+LFnDJJfDMM6BLkCLRlpAA99zjrxzd\nfXfoakSkKBdfDM2bwz/+AVu2hK4m0tRES3waMQIaNvQhj4Nx/SJVWrducMIJcNttsHRp6GpEpDAp\nKXD77fDll/5kleyQmmiJT3Xr+hnEU6fC66+HrkZEinL77bBpE/zzn6ErEZGinHKKn4N01VVaDasQ\naqIlfp19NrRt6ydAbNwYuhoRKcyee8KFF8KTT/ozXCISXWZ++NXSpX4iv2yXmmiJX0lJ8K9/wdy5\n8OCDoasRkaJcfTXUrw+XXqphWCJR17mzVsMqgppoiW89e0Lv3n5Znl9/DV2NiBQmNRWuuw7efx/e\nfjt0NSJSlLzVsK68MnQlkaQmWuLfXXf5MVtVePMVkbjx97/D3nvDZZf5MdIiEl15q2E9+6xWw9oO\nNdES/9q0gcGD4eGH4fvvQ1cjIoVJTvYffP/3Pxg5MnQ1IlKUESNgt918M61hWNtQEy2Vw/XX+93R\nLr88dCUiUpSjjoIePfzQjt9/D12NiBQmbzWsjz6C114LXU2kqImWyqFhQ7901ttvw6RJoasRkcKY\n+UnB2dn+4Cwi0XbWWVoNazvUREvlceGF0LKln/mvXZZEom3//f0ylf/+tx/aISLRlbca1rx5PrMC\nVEATbWaXmpkzswax783M7jezOWb2lZkdWN41SBVRo4bf0OGrr/xatFJsyqtUqBtv9LkdPjx0JXFJ\neZUK1bOnH4p1442wfHnoaiKhXJtoM0sDegIL893dG9grdhsMPFyeNUgVc9JJ0KWLH9qxenXoauKK\n8ioVrlEjP2npjTfgww9DVxNXlFcJ4s47tRpWPuV9JvoeYBiQfzpnH+AZ500DUs2scTnXIVWFGdxz\nD/z8sz8rLcWhvErFu+QSaN4c/vEPDcMqHuVVKl6bNn6ZypEj4bvvQlcTXLk10WbWB1jinJtV4KGm\nwKJ83y+O3Vfw5web2Qwzm7Fclw2kODp0gIED/fithQuLfr6UOq+x11BmpfhSUuC22+CLL/xatFIk\n5VWCuu46qF3br/VexZWqiTazSWb2zXZufYArgWtK+trOuUedcxnOuYyGDRuWpkypim65xf9Xuyz9\noTzzCsqslEK/ftCxo8/r2rWhq4kE5VUiK281rHHjYMKE0NUEVaom2jnXwznXtuANmAe0BGaZ2QKg\nGfC5mTUClgBp+V6mWew+kbLTvLlfpeP55+HTT0NXEwnKq0SWGdx9Nyxd6sdcivIq0XbhhbDHHlV+\nNaxyGc7hnPvaObebc66Fc64F/pLSgc65ZcCbwGmxWcSdgJXOuaXlUYdUccOHw+67+7GW2mVph5RX\niYQuXeCUU+COO2Dx4tDVRJbyKpFQvbrP6jffwBNPhK4mmBDrRI/Df5KeAzwGnBegBqkK6tSBm26C\njz+GV18NXU28Ul6l4tx2G+TmwlVXha4kXimvUnFOOAEOOQSuvhpWrQpdTRAV0kTHPjH/GvvaOefO\nd861cs7t55ybURE1SBV15pl+U4fhw2HDhtDVxAXlVYJp0QIuvhieeQZm6K/azlBeJZi8YVi//OI/\nAFdB2rFQKrfERL9Kx/z58MADoasRkaJceaWfuKRhWCLR1749nHqqb6Z/+il0NRVOTbRUfj16wDHH\n+KEdWspJJNrq1vU7ok2dCq+/HroaESnKLbdAQgJccUXoSiqcmmipGu680y+dde21oSsRkaKcfTbs\nuy8MGwYbN4auRkQKk5bm14wePRr++9/Q1VQoNdFSNey9N5x7Lu6RR3jggTfo8++PuGbsN2Rlrw9d\nmYgUlJTkh2HNncu75/5TeRWJumHD2NKoEYvOHEKfB6ZWmbyqiZYqY9nFw1idnMJ+997MrMUreWH6\nQnrfN7VKBF0k3mR17MbUPdvT5YWHWfi/hcqrSIRl5SRyfccBpP3wFc0nvlVl8qomWqqMh75ZyYMH\n9+PQeTPo/NNX5OQ61m3MYeTkuaFLE5ECRk6ey82HnUXNTes5/78vK68iETZy8lxG730os3fbg2FT\nnoHNm6tEXtVES5Uxa1E2Tx1wNMtrpXLutFcA2JzrmLUoO3BlIlLQrEXZfF8/jbFtujNg1rukrl+l\nvIpE1KxF2Wwigbu7DiRt5c8c8/3UKpFXNdFSZaSnpbKlWnWePOhYui34gn1/nktygpGelhq6NBEp\nID0tlaQE45GOJ1Jz80ZO+/wd5VUkovLy+n6r9vzQoDl/nz6GZKPS51VNtFQZQ7q3olb1JF486GhW\nV0thyKevUbN6EkO6twpdmogUkJfXebu3ZFKr9pwx8y3qW47yKhJBeXlNTEzkkY4nss/yBfRc+EWl\nz6uaaKkymqSmMH5oV/6vexsmde3D0d9PZcJxaTRJTQldmogUkJfXAR2bM+n/zqD++lW8V2+O8ioS\nQfnzuuCIY8lu0Ii7F7xX6fNqLg52hMrIyHAztAWslKUlS6BlSxg8GP7979DVFJuZzXTOZYSuY0eU\nWSlzBx8MWVnw449+Cbw4E+XMKq9S5u69Fy65xK8b3alT6GqKbWfzqjPRUjU1bQqDBsGoUdrFUCQe\nDB8OCxbAyy+HrkREinLOObDLLnD77aErKVdqoqXquvxy2LABHnggdCUiUpRjjoE2beCOOyAOrqCK\nVGm1a8MFF8DYsfD996GrKTdqoqXq2ntv6NPHD+dYsyZ0NSJSmIQE/8F31ix4773Q1YhIUS68EGrU\ngDvvDF1JuVETLVXb8OGwYgU8/njoSkSkKAMGQLNmcNttoSsRkaI0bAhnnQXPPuvnIVVCaqKlauvU\nCbp1g7vvhs2bQ1cjIoWpVs1PVpo8GaZPD12NiBTl0kshN9dPNKyE1ESLDB8OixbBiy+GrkREipKZ\nCamplX7Ckkil0LIlnHwyjBzpr/pWMmqiRXr3hrZt/YSl3NzQ1YhIYerUgfPPhzfegB9+CF2NiBRl\n2DA/7+jhh0NXUubURIuY+ZDPng3jxoWuRkSKctFFUL16pZ6wJFJptGsHRx4J990H69eHrqZMqYkW\nAejXD5o392ejRSTadttt64SlrKzQ1YhIUYYPh19+gaefDl1JmVITLQKQnAz/+AdMnep3WBKRaLv0\nUsjJqbQTlkQqlUMPhfbt4a67YMuW0NWUGTXRInnOOQfq19eEJZF4sMceWycsZWeHrkZECmPmz0bP\nnQtjxoSupsyoiRbJU6vW1h2WvvsudDUiUpRhw2D1at9Ii0i0HXcctG7tT1RVkl1H1USL5HfhhZCS\noglLIvHggAOgZ08/pGPDhtDViEhhEhP9rqOffw7/+U/oasqEmmiR/Bo0gLPPhueeg8WLQ1cjIkUZ\nPhx+/rnSTVgSqZQGDYLGjSvNsEk10SIFVfIdlkQqlcMOg4yMSjdhSaRSql4dLr4YJk2CmTNDV1Nq\naqJFCmrRAk45BR55pFLusCRSqeRNWJozB157LXQ1IlKUv/8d6tatFEvKqokW2Z5KvMOSSKVz/PGw\n116VasKSSKVVrx6cey68+qpfrSOOlWsTbWYXmtn3ZjbbzO7Id/8IM5tjZj+Y2ZHlWYNIiaSnQ69e\nlXKHpR1RXiVuJSbCZZf5y8Pvvx+6mgqhvEpcGzoUkpL8MKw4Vm5NtJkdBvQB0p1z+wJ3xe5vA/QD\n9gV6AQ+ZWWJ51SFSYsOGVcodlrZHeZW4d9pp0KhRpZmwVBjlVeJe48Zw+unw5JN+YnCcKs8z0ecC\ntznnNgI4536J3d8HGO2c2+icmw/MATqUYx0iJVNJd1jaAeVV4luNGn7C0sSJfgmtyk15lfh32WWw\naZO/4hunyrOJbg10NbPpZjbZzNrH7m8KLMr3vMWx+7ZhZoPNbIaZzVi+fHk5limyA/l3WKr8E5ZK\nlVdQZiUChgypNBOWiqC8Svxr3RpOOAEeeghWrQpdTYmUqok2s0lm9s12bn2AJKA+0Am4HHjZzGxn\nX9s596hzLsM5l9GwYcPSlClScscdV2kmLJVnXkGZlQioV8830q+8EvcTlpRXqRKGD4eVK+HRR0NX\nUiKlaqKdcz2cc223cxuL/wT8mvM+BXKBBsASIC3fyzSL3ScSPXk7LFWCCUvKq1QJF1/sJyz961+h\nKykV5VWqhPbt/Vrv99wDGzeGrqbYynM4xxvAYQBm1hqoBvwKvAn0M7PqZtYS2Av4tBzrECmdQYOq\nwoQl5VUqh8aN/STDOJ+wVATlVSqP4cMhKwuefz50JcVWnk30KGAPM/sGGA2cHvvUPBt4GfgWeBc4\n3zlX6WdtSRyrGhOWlFepPC6/3J/Vuv/+0JWUF+VVKo+ePaFdOz+XITc3dDXFYi4OxnlmZGS4GTNm\nhC5DqrKVK6F5c+jdG0aPDl0NZjbTOZcRuo4dUWYluBNP9EOwFi6EOnVCVxPpzCqvEtyLL8KAAfD6\n634uUmA7m1ftWCiyMyrRhCWRKmH4cMjOjtsJSyJVSt++0LJl3E3iVxMtsrMqyYQlkSqhQwe/1vs9\n9/i1aEUkupKS/LrR06bB1Kmhq9lpaqJFdlb+CUu//FL080UkrOHDYcmSuJywJFLlnHkmNGwYV5P4\n1USLFMdll1X2CUsilceRR0J6elxOWBKpclJS4KKLYNw4+Prr0NXsFDXRIsXx17/6SQ8PPghr1oSu\nRkQKk7fr6Pffw1tvha5GRIpy3nlQq1bc7DqqJlqkuPImLD32WOhKRKQocTphSaRKql8fBg/2q3X8\n9FPoaoqUFLoAkbjTsSN0786Wu/7FTWmH8vmytaSnpTKkeyuapKaErk5E8ktKgksvhQsu4Nfx/+H+\nzY2YtShbmRWJqksugQceYM2td3BH73MjnVediRYpgd/Ov5jErCWseepZZi1eyQvTF9L7vqlkZa8P\nXZqIFHTmmWzZtQHfXnwVL0xfqMyKRFlaGuv69iPxiScY/8HXkc6rmmiRErgvuRXf7daSzGljMJdL\nTq5j3cYcRk7WGtIikVOzJh/2PIVuP35Kq5/nAyizIhH2RJeTSMnZyMAZfi5DVPOqJlqkBGYtXsnI\nDifQ+reF/G3uZwBsznXMWpQduDIR2Z5R6UexNrkGf58+5o/7lFmRaJrErkzcswOnz3yblE0bgGjm\nVU20SAmkp6Xy7r7dWVx3N4ZM8wfl5AQjPS01cGUisj2t9m7Oy+2O5NhvJ9N0pV/nXZkViab0tFQe\n69yXXTas5pSvJgDRzKuaaJESGNK9FTVSqjOqw/G0X/It7X6eQ83qSQzp3ip0aSKyHUO6t+LFQ04C\nYNAX75CcYMqsSEQN6d6KH1rtz8ym+3DGzLeoZi6SeVUTLVICTVJTGD+0K0lnnMbG5OpctfRjxg/t\nGrmZwyLiNUlN4alrTuTH9t3p9+37nHpgY2VWJKLyjrHz+w6iRfZSRtT8OZJ5VRMtUkJNUlO4cmAX\nqg/oR/tP3qVJ0pbQJYlIIZqkprDPNZeSunoF19q8yB2QRWSrJqkpnHTLxZCaypnfTopkXtVEi5RW\nZiasXg0vvxy6EhEpSq9e0KwZPPpo6EpEpCgpKXDaafDaa/Drr6Gr+RM10SKl1aUL7LOPDsoi8SAx\nEc4+GyZOhAULQlcjIkXJzIRNm+CZZ0JX8idqokVKy8yHfPp0+Prr0NWISFHOOsv/94knwtYhIkVr\n2xY6d4bHHgPnQlezDTXRImVh0CCoVs2HXESirXlz6N0bRo2CnJzQ1YhIUTIz4fvv4aOPQleyDTXR\nImWhQQM48UR49llYH61tSUVkOzIzISsLxo0LXYmIFOXkk6Fu3cidqFITLVJWMjMhOxvGjCn6uSIS\n1tFHQ6NGkTsoi8h21KoFAwfCK6/AihWhq/mDmmiRsnLoobDnnppgKBIPkpP92Ohx42Dx4tDViEhR\nMjNhwwZ47rnQlfxBTbRIWTGDc86BqVP92C0Ribazz4bcXD82WkSi7YAD4KCDIjXBUE20SFk64wxI\nSoLHHw9diYgUZY894Igj/CodW7RZkkjkDR7sV8GaPj10JYCaaJGytfvu0KcPPP00bNwYuhoRKUpm\nJixcCBMmhK5ERIrSv78fHx2RuQxqokXKWmam31lp7NjQlYhIUfr0gYYNI3NQFpFC1KnjG+nRo2HV\nqtDVqIkWKXNHHAF/+YsmGIrEg2rV/DCst96CpUtDVyMiRcnMhHXr4IUXQleiJlqkzCUk+AlL//kP\nzJ0buhoRKco55/hNV556KnQlIlKU9u1h//0jcfWo3JpoM2tnZtPM7Eszm2FmHWL3m5ndb2ZzzOwr\nMzuwvGoQCeass3wzHSfbCiuvUqW1bg3du/sJwbm5oaspkvIqVZqZn2D4+ecwc2bQUsrzTPQdwPXO\nuXbANbHvAXoDe8Vug4GHy7EGkTCaNvWbOTz5JGzeHLqanaG8StU2eDDMmwcffBC6kp2hvErVNnAg\npKQEPxtdnk20A+rGvq4HZMW+7gM847xpQKqZNS7HOkTCyMyEZcvgnXdCV7IzlFep2k44AerXj5e5\nDMqrVG2pqdC3rx8XvWZNsDLKs4m+GLjTzBYBdwEjYvc3BRble97i2H3bMLPBsctUM5YvX16OZYqU\nk969oUmTeDkolyqvoMxKnKtRAwYNgtdfh+j//VVeRQYPhtWr4aWXgpVQqibazCaZ2TfbufUBzgUu\ncc6lAZcAxRoc6px71DmX4ZzLaNiwYWnKFAkjKcmPjX73Xb8ObWDlmVdQZqUSyMz0w6+efjp0Jcqr\nSFG6dIF99gk6pKNUTbRzrodzru12bmOB04HXYk99BegQ+3oJkJbvZZrF7hOpfM4+2/83AtsKK68i\nRdh3X39gfvzx4NsKK68iRTDzH3ynT4evvgpSQnkO58gCuse+/hvwY+zrN4HTYrOIOwErnXNanFMq\npxYtoGdP30RHe1th5VUE/EH5hx9g6tTQlRRGeRUBOO00v9Z7oLPR5dlEZwL/MrNZwC34mcIA44B5\nwBzgMeC8cqxBJLzMTFi0CN57L3QlhVFeRQBOPhnq1Yv6XAblVQRg113hxBPhuedg/foKf/uk8nph\n59xHwEHbud8B55fX+4pEzv/9H+y2mz8oH3VU6Gq2S3kVialZ0y+f9cQTcP/9fsWOiFFeRfIZPBhe\nfBFefdVPDq5A2rFQpLzlbSv89tvaVlgkHmRmwsaN/uyWiERb9+6w115Brh6piRapCOec48dEP/lk\n6EpEpCjt2vmthR99NPgEQxEpgpk/xn70EXz3XYW+tZpokYqw115w2GFxs62wSJWXmQmzZ8O0aaEr\nEZGinHEGJCf7Y2wFUhMtUlEyM2H+fPjPf0JXIiJF6dcPatUKvq2wiOyE3XaDPn38Gu8bN1bY26qJ\nFqkoxx/vJynpoCwSfXXqwIABMHo0rFwZuhoRKUpmJvz2m991tIKoiRapKDVq+DUt33gjHrYVFpHM\nTL9s1gsvhK5ERIrSo4ffm6ECJxiqiRapSBHaVlhEipCR4ScZaoKhSPQlJPgJhh98AHPmVMxbVsi7\niIjXpg0cfLAf0qGDski05W0r/OWXMHNm6GpEpChnngmJiRU2wVBNtEhFy8yE//0PpkwJXYmIFGXg\nQEhJ0VwGkXjQpAkcc4xfTnbTpnJ/OzXRIhWtb1+/rbAOyiLRV68enHKKHxe9Zk3oakSkKJmZ8Msv\n8NZb5f5WaqJFKlretsKvvgq//x66GhEpSmamb6BHjw5diYgUpVcvaNasQk5UqYkWCWHwYG0rLBIv\nOnf28xl09Ugk+hIT4eyzYcIEWLCgXN9KTbRICOnpflthTTAUiT4z/8H3009h1qzQ1YhIUc46y//3\niSfK9W3URIuEkpkJ33yjbYVF4sGgQVC9us5Gi8SD5s2hd28YNQpycsrtbdREi4SibYVF4kf9+nDi\niX4I1rp1oasRkaJkZkJWFowbV25voSZaJJQ6daB/f3jpJVi1KnQ1IlKUwYP9FuCvvhq6EhEpytFH\nQ6NG5XqiSk20SEiDB/uzWtpWWCT6unWD1q0rdFthESmh5GQ/NnrcOFi8uFzeQk20SEgZGX6SoYZ0\niESfmd9W+OOP4dtvQ1cjIkU5+2zIzfVjo8uBmmiRkPK2Ff78c20rLBIPTj/dn+HSB1+R6NtjDzji\nCL9Kx5YtZf7yaqJFQtO2wiLxY7fd4Ljj4JlnYMOG0NWISFEyM2HhQr9udBlTEy0SWmqq3wpc2wqL\nxIfMTL/b6Ouvh65ERIrSpw80bFguJ6rURItEweDBsHo1vPxy6EpEpCiHHw4tW2qCoUg8qFYNzjgD\n3noLli4t05dWEy0SBV26wD77aEiHSDxISPATDD/8EH78MXQ1IlKUc87xm6489VSZvqyaaJEoMGPl\nqWfAtGlcdOXTXDP2G7Ky14euSkR2YNkJ/diSkMiY865TXkWirnVrNh7cld/ve4jj7p9SZplVEy0S\nAVnZ6zk2uyWbEpM4YMIYXpi+kN73TdWBWSSCsrLXc+RLc/jPnh3o/sk7vPLJXOVVJMKystdz7e5d\nqP/zYmp9MrXMjrFqokUiYOTkuSxJqsW7rbtwwjfvk7BpE+s25jBy8tzQpYlIASMnz2Xtxhxe2L8n\nDdatpNO8L5RXkQgbOXkub7bqxIoadTh+9gfk5LoyyWxSGdUnIqUwa1E2ObmO+7r05+FOJ7EpKRly\nHbMWZYcuTUQKyMvrlJYHcsKpd/J5k72VV5EIm7Uom3UJyQw65UZ+3DUNgM1lkFmdiRaJgPS0VJIS\njLkN0vhutz0ASE4w0tNSA1cmIgXl5TU3IZHPm+4DZsqrSITlZfabRnuyMbk6UDbH2FI10WbW18xm\nm1mumWUUeGyEmc0xsx/M7Mh89/eK3TfHzK4ozfuLVBZDureiVvUkkhIM8OGuWT2JId1blen7KLMi\npae8isSX8spsaYdzfAOcADyS/04zawP0A/YFmgCTzKx17OEHgSOAxcBnZvamc+7bUtYhEteapKYw\nfmhXRk6ey6xF2aSnpTKkeyuapKaU9VspsyKlpLyKxJfyymypmmjn3HcAZlbwoT7AaOfcRmC+mc0B\nOsQem+Ocmxf7udGx5yrgUuU1SU3hhj5ty/U9lFmRsqG8isSX8shseY2Jbgosyvf94th9O7r/T8xs\nsJnNMLMZy5cvL6cyRSRGmRWJH8qrSAQUeSbazCYBjbbz0FXOubFlX5LnnHsUeBQgIyPDldf7iFQ2\nyqxI/FBeReJXkU20c65HCV53CZCW7/tmsfso5H4RKQPKrEj8UF5F4ld5Ded4E+hnZtXNrCWwF/Ap\n8Bmwl5m1NLNq+IkRb5ZTDSKy85RZkfihvIpEQKkmFprZ8cADQEPgHTP70jl3pHNutpm9jJ/MkAOc\n75zbEvuZC4D3gERglHNudql+AxHZacqsSPxQXkWizZyL/lAoM1sO/LQTT20A/FrO5ZRUlGsD1Vca\nIWr7i3OuYQW/507bycxG+f8pRLu+KNcGqm97IptZ5bVCRLm+KNcGEc5rXDTRO8vMZjjnMop+ZsWL\ncm2g+kojyrVFWdT/3KJcX5RrA9VXGUX9z0z1lVyUa4No16dtv0VEREREiklNtIiIiIhIMVW2JvrR\n0AUUIsq1georjSjXFmVR/3OLcn1Rrg1UX2UU9T8z1VdyUa4NIlxfpRoTLSIiIiJSESrbmWgRERER\nkXKnJlpEREREpJgqRRNtZr3M7Aczm2NmV4SuJz8zSzOzD8zsWzObbWZDQ9dUkJklmtkXZvZ26FoK\nMrNUM3vVzL43s+/MrHPomvIzs0ti/1+/MbMXzaxG6JriQVQzGw95BWW2pJTXkolqXiE+Mqu8lkw8\n5DXum2gzSwQeBHoDbYD+ZtYmbFXbyAEudc61AToB50esPoChwHehi9iB+4B3nXN7A+lEqE4zawpc\nBEDmVIEAAAKwSURBVGQ459ridwjrF7aq6It4ZuMhr6DMFpvyWjIRzyvER2aV12KKl7zGfRMNdADm\nOOfmOec2AaOBPoFr+oNzbqlz7vPY16vxf0Gbhq1qKzNrBhwNPB66loLMrB7QDXgCwDm3yTmXHbaq\nP0kCUswsCagJZAWuJx5ENrNRzysos6WkvBZfZPMK0c+s8loqkc9rZWiimwKL8n2/mAgFKD8zawEc\nAEwPW8k27gWGAbmhC9mOlsBy4MnYpbDHzaxW6KLyOOeWAHcBC4GlwErn3ISwVcWFuMhsRPMKymyJ\nKK8lFhd5hchmVnktgXjJa2VoouOCmdUGxgAXO+dWha4HwMyOAX5xzs0MXcsOJAEHAg875w4A1gKR\nGY9nZrvgz8i0BJoAtczs1LBVSVmIYl5BmS0N5bVyi2JmldeSi5e8VoYmegmQlu/7ZrH7IsPMkvHh\nft4591roevI5GDjWzBbgL9H9zcyeC1vSNhYDi51zeWcVXsUHPip6APOdc8udc5uB14AugWuKB5HO\nbITzCspsaSivJRPpvEKkM6u8llxc5LUyNNGfAXuZWUszq4YfeP5m4Jr+YGaGH2/0nXPu7tD15Oec\nG+Gca+aca4H/c3vfOReZT3rOuWXAIjP7a+yuw4FvA5ZU0EKgk5nVjP1/PpyITMqIuMhmNsp5BWW2\nlJTXkolsXiHamVVeSyUu8poUuoDScs7lmNkFwHv42ZujnHOzA5eV38HAIOBrM/sydt+VzrlxAWuK\nJxcCz8f+8Z4HnBm4nj8456ab2avA5/gZ4l8Q4e1JoyLimVVeSy+SmVVeSybieQVltrSU11LQtt8i\nIiIiIsVUGYZziIiIiIhUKDXRIiIiIiLFpCZaRERERKSY1ESLiIiIiBSTmmgRERERkWJSEy0iIiIi\nUkxqokVEREREiun/Af2HfoPw9BxPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107217f60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "titles = ['svd', 'elementary', 'matrix']\n",
    "\n",
    "plt.figure(figsize=(12,4))\n",
    "for i, beta in enumerate([beta1, beta2, beta3], 1):\n",
    "    plt.subplot(1, 3, i)\n",
    "    plt.scatter(x, y, s=30)\n",
    "    plt.plot(x, beta[0] + beta[1]*x + beta[2]*x**2, color='red')\n",
    "    plt.title(titles[i-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Curve fitting and least squares optimization\n",
    "\n",
    "As shown above, least squares optimization is the technique most associated with curve fitting. For convenience, `scipy.optimize` provides a `curve_fit` function that uses Levenberg-Marquadt for minimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logistic4(x, a, b, c, d):\n",
    "    \"\"\"The four paramter logistic function is often used to fit dose-response relationships.\"\"\"\n",
    "    return ((a-d)/(1.0+((x/c)**b))) + d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "nobs = 24\n",
    "xdata = np.linspace(0.5, 3.5, nobs)\n",
    "ptrue = [10, 3, 1.5, 12]\n",
    "ydata = logistic4(xdata, *ptrue) + 0.5*np.random.random(nobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "popt, pcov = curve_fit(logistic4, xdata, ydata) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Param\tTrue\tEstim (+/- 1 SD)\n",
      "a\t10.00\t10.40 (+/- 0.14)\n",
      "b\t 3.00\t 4.20 (+/- 1.16)\n",
      "c\t 1.50\t 1.38 (+/- 0.09)\n",
      "d\t12.00\t12.03 (+/- 0.10)\n"
     ]
    }
   ],
   "source": [
    "perr = yerr=np.sqrt(np.diag(pcov))\n",
    "print('Param\\tTrue\\tEstim (+/- 1 SD)')\n",
    "for p, pt, po, pe  in zip('abcd', ptrue, popt, perr):\n",
    "    print('%s\\t%5.2f\\t%5.2f (+/-%5.2f)' % (p, pt, po, pe))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFNW5x/HvO8MAI9uwiTCIgBoURQEnqCFGxSiIXgSi\nuWpiMJKLWbyJuUrEbHpzzQVDot4kRoMbagwuUXABg8SNLG6DIKsI4saADIsswjDM8t4/qkaaZpae\npbu6Z36f56mnq09Vdb9d0P1OnXPqHHN3REREsqIOQERE0oMSgoiIAEoIIiISUkIQERFACUFEREJK\nCCIiAighiIhISAlBREQAJQQREQm1ijqA+ujWrZv37ds36jBERDLKokWLtrh797r2y6iE0LdvXwoL\nC6MOQ0Qko5jZB4nspyojEREBlBBERCSkhCAiIoASgoiIhJQQREQEyLBeRiIiVeYsLmL6/NVs2F5C\nr7xcJo8cwNgh+VGHldGUEEQk48xZXMT1TyyjpKwCgKLtJVz/xDIAJYVGUJWRiGSc6fNXf5YMqpSU\nVTB9/uqIImoelBBEJONs2F5Sr3JJjBKCiGScXnm59SqXxNSZEMzsXjMrNrPlMWXTzextM1tqZrPN\nLK+a4w43sxfNbKWZrTCzH8Rsu9HMisxsSbiMbrqPJCLN3eSRA8jNyT6gLDcnm8kjB0QUUfOQyBXC\nTGBUXNkC4Hh3PwF4B7i+muPKgWvcfSBwCvA9MxsYs/1Wdx8cLvPqH7qItFRjh+Qzdfwg8vNyMSA/\nL5ep4wepQbmR6uxl5O4LzaxvXNlzMU9fBS6s5riNwMZwfZeZrQLygZWNiFdEBAiSghJA02qKNoQr\ngGdr2yFMKEOA12KKrwqrnO41s85NEIeIiDRCoxKCmf2EoGrooVr2aQ88Dlzt7jvD4juAI4HBBFcR\nv6nl+ElmVmhmhZs3b25MuCIiUosGJwQzuxw4H/iau3sN++QQJIOH3P2JqnJ33+TuFe5eCdwFDKvp\nfdx9hrsXuHtB9+51zu8gIiIN1KCEYGajgB8BY9x9Tw37GHAPsMrdb4nb1jPm6ThgOSIiEqlEup3O\nAl4BBpjZejObCPwe6AAsCLuN3hnu28vMqnoMDQcuA0ZU0730V2a2zMyWAmcCP2zizyUiIvVkNdT2\npKWCggLXFJoiIvVjZovcvaCu/XSnsoiIAEoIIiISUkIQERFACUFEREJKCCIiAmjGNBFJM5oaMzpK\nCCKSNjQ1ZrRUZSQiaUNTY0ZLCUFE0oamxoyWEoKIpA1NjRktJQQRSRuaGjNaalQWkbRR1XCsXkbR\nUEIQkbSiqTGjoyojEREBdIUgkvZ0o1by6RwHErpCMLN7zazYzJbHlE03s7fNbKmZzTazvBqOHWVm\nq81srZlNiSnvZ2avheWPmFnrxn8ckeal6katou0lOPtv1JqzuCjq0JoNneP9Eq0ymgmMiitbABzv\n7icA7wDXxx9kZtnA7cC5wEDgEjMbGG6+GbjV3Y8CPgEm1jt6kWZON2oln87xfglVGbn7QjPrG1f2\nXMzTV4ELqzl0GLDW3dcBmNnDwAVmtgoYAVwa7nc/cCNwRz1iF2n2dKNW8jX6HFeUQdkeKCsJlvK9\n+x/L90LZ3v3r5XuhvDR83Bc8VuyLKS+FitJgW0XMUl4Ko6fD4cOa8JMfrKnaEK4AHqmmPB/4KOb5\neuBkoCuw3d3LY8pbXoWdSB165eVSVM0Pk27UStz+9oE99O2UzY/O6Mm5R7eH0p1QuouvdlhK6e6d\ntLcSDmEv7ayUQ9hL9zbl8Phs2Lc7WMr2wL49ULZ7/49/2R6oLK87iJpk5UCrNpDdGlq1hVatIbvq\nebjeqi207QRZyW/ybfQ7mNlPgHLgocaHU+3rTwImAfTp0ycZbyGStiaPHHDAYG/Qwm/Ucod9n8Ke\nbVCyDUo+iVu2w97tsHcH7N3B9m1bGLp9K0+zmw5tSsgprYD5BEvoZoC4Fszd3pbsVu1hfQdo3R5a\nHxI8tu8BOYdATi60bhc8tsoNHquWVm33P7ZqCzltg31atQnKs1sH663aQtaBN+FFrVEJwcwuB84H\nznJ3r2aXIuDwmOe9w7KtQJ6ZtQqvEqrKD+LuM4AZAAUFBdW9h0iz1SJu1Kooh93F8Okm+LQ4WHYX\nw+4tsHtzuGyFPVthz5agCqUmrXIhNy/4i7ptHit25bK58kh2ejt2kcunfgi7yCXnkE7c8JWToU0H\naNOeBe/u4bf/+Jh1O6Bzp05cO+rY5nWOE9TghGBmo4AfAae7+54adnsDONrM+hH84F8MXOrubmYv\nErQ7PAxMAJ5saCwizVlG36i1bw/sLIIdH8GOIti5IXi+a2O4fBz88FPN33qtO0C7rtCuO3TKh54n\nBs9zu8AhXeGQLpDbOXie2zlIBK3aHPASX58yt7pXxnbBDceM/uz52b3g7NOa9qNnooQSgpnNAs4A\nupnZeuAGgl5FbYAFZgbwqrt/28x6AXe7+2h3Lzezqwgu0LKBe919Rfiy1wEPm9lNwGLgnib8XCJS\ngybtc19RBts/hE/e379s/xC2fxA87tl68DHtukOHntChF/QaGq73CKpj2vcItrc/NKheaSS1wdSP\nVV/Tk54KCgq8sLAw6jBEMlb8BDQQtElMHT+o5qTgHvwVv2U1bHkHtqyFrWtg67vBD39Mo+o+WlHa\nrjcdDusPeX2g0+Hhkg8d86Fjr4P+ik+mBn3eZsjMFrl7QV376U5lkRaktj73Y4fkQ+ku2LQSNi2H\nTSugeBVsfjtowK3SKhe6HgWHHc/qriN4YHU2a8q686EfyiY607Yih6lfTo8f3BbRBtOEdIUg0oL0\ni6lT78huBmWt4wR7j+Oz3uO87pth27r9O7fpCIceC92PCR67fS5YOuZDVnBP6/BpL1RbJZOfl8s/\np4xIwSeSROgKQUT2q6yA4pV8t/3L9C9dyWBby5FZGz/bXGQ9oMcwOPFSOGwQ9DgOOvWGoH2wRrpx\nrnlRQhBpjirKYeMSeG8hfPAv+Oh1KN3BZGBLVicWVx7F42WnsdSPZG32UUwZfyr5DahGUaNt86KE\nINIcuMOWNfDuC8Hywb9g365gW7cBcPw46HMqHH4y/3g/h+nPvfNZnfqURtSp68a55kUJQSRTlZXA\ne3+Hd/4KaxbAjg+D8i794YSLoO9pwdK++wGHje0CY4f2bpIQ1GjbvCghiGSSku3wznxY9VRwJVC2\nB3LaQf8z4LQfwpEjoHPflIaU0TfOyQGUEETSXekueHseLH88SAKVZcHNXIMvhQHnwhFfDMbLEWkk\nJQSRBKV0Vq2Kclj3Iiz5M6yeFwyP3LE3nHwlDBwL+Sd91vVTpKkoIYgkIP6O16pZtYCmTQrb3oM3\nHwgSwacfB2P0DPk6DLoIeg9TEpCkUkIQSUCdd/g2RmVF0DD8xt1BlZBlwdEjgyqhz41M6VAP0rIp\nIYgkICk3YO3dAW8+CK/PCMYE6pgPZ/w4uCLopEZaST0lBJEENOkNWLs+hlf/AIX3BbN29fkCnPM/\nMOA8yNZXUqKj/30iCWiSG7B2rIe//wYW/ykYIXTgWBj+feg1JAkRN15TN6KntFFeGkQJQSQB9bkB\nK/6H72df6sSobQ8FjcUAQ74Gw38Q3ECWppq6ET1ljfLSKBrtVKQJxf7wdWQ332n1FN/M/iuts5ys\noZfBaddA3uF1v1DEmnoUU42KGq1ERzutsw+bmd1rZsVmtjym7CIzW2FmlWZW7ZuY2QAzWxKz7DSz\nq8NtN5pZUcy20dW9hkimmT5/NWVlpVyR/Swvt/khV2Y/w7zKk7ko53b4t9syIhlA0zeia1TUzJBI\nldFM4PfAAzFly4HxwB9rOsjdVwODAcwsm2BO5dkxu9zq7r+uZ7wiaa3/zteZ2foBjs4qYmHFIKaV\nX8JK74uVRR1Z/TT1KKYaFTUz1HmF4O4LgW1xZavCH/xEnQW86+4f1DM+kcywowge/hoPtp5KDuVc\nse9avlE2hZXeF8i8H77JIweQm5N9QFljRjFt6teT5EhVo/LFwKy4sqvM7BtAIXCNu39S3YFmNgmY\nBNCnT5+kBilSb5UV8MY98PwvoLKcFcdezSXLC9hZuf9vrUz84WvqUUw1KmpmSKhR2cz6As+4+/Fx\n5S8B17p7jS29ZtYa2AAc5+6bwrIewBbAgf8Berr7FXXFoUZlSStb1sCc78D6N4JRRs+7Bbr0U/dK\nSTvpNIXmucCbVckAIHbdzO4CnklBHCJNo7ISXv8j/O1GyMmFcTPghK9+Nt2khoOWTJWKhHAJcdVF\nZtbT3asmdB1H0Egtkv52FMHsK+H9vwfjDY35LXQ4LOqoRJpEnQnBzGYBZwDdzGw9cANBI/PvgO7A\nXDNb4u4jzawXcLe7jw6PbQecDVwZ97K/MrPBBFVG71ezXST9rH42qCIq3wdjfgdDLqtzEnqRTFJn\nQnD3S2rYNDu+wN03AKNjnu8Gulaz32X1iFEkqeqs8y/fBwt+Dq/dAYcNggtnQrejIou3Jmq7kMbS\n0BXSotU5pMLOjfDYBPjoNRh2ZTAIXRoOR62hIaQpaLYNadFqm+eAD/4FM06Hj5fDhffB6F+lZTKA\nOj6HSIJ0hSAtWk1DJ5y2ax7cfx/k9YHL5kCPgSmOrH40NIQ0BV0hSIsWfwdxFpVc3+ohpuXcBf2+\nBP/xYtonA6j5TuhMu0NaoqWEIC1a7JAKbSnljpzbuLLVXNb1vQQufQxy8yKOMDEaGkKagqqMpEWr\nanC986+F/G/JLzkxax1Lj7+eE75yXUZ1KdXQENIUNB+CyI718OB4+OQ9+Mo9MHBM1BGJNKl0GrpC\nJH1tWQMPXAClu+DrT0C/06KOSCQySgjSchWvgvvHgFfC5XOh5wlRRyQSKSUEaZk+XhZcGWTlwOXP\nQHc1voooIUjLs3EpPDAGcg6BCU9D1yOjjkgkLSghSMtSvAoeHAut2wdXBp37Rh2RSNpQQpCM1KCB\n3La+u7+a6BtPKhmIxFFCkIzToIHctn8YNCBXVsA356VFNZFGJ5V0ozuVJePUeyC33VvhwXGwbxdc\nNjstGpCrklrR9hKc/UltzuKiqEOTFqzOhGBm95pZsZktjym7yMxWmFmlmdV4s4OZvW9my8xsiZkV\nxpR3MbMFZrYmfOzc+I8iLUW9BnLbtxv+/NXg5rNLHkmbrqUanVTSUSJXCDOBUXFly4HxwMIEjj/T\n3QfH3SU3BXje3Y8Gng+fiyQk4YHcKsrgscthw5vBHchHnJr84BKk0UklHdWZENx9IcGUmbFlq9y9\nMX/KXADcH67fD4xtxGtJC5PQQG7uMPe/YM1zcN4tcOz5KY6ydhqdVNJRstsQHHjOzBaZ2aSY8h7u\nvjFc/xjokeQ4pBkZOySfqeMHkZ+XiwH5eblMHT/owAbZf/0O3nwATrsWCr4ZWaw10eikko6S3cvo\ni+5eZGaHAgvM7O3wiuMz7u5mVuMIe2EimQTQp0+f5EYrGWPskPyae+SseiaYA3ngWDjzJ6kNLEEa\nnVTSUVITgrsXhY/FZjYbGEbQ7rDJzHq6+0Yz6wkU1/IaM4AZEIx2msx4pRnY+BY88R+QPxTG3QlZ\n6duRrtakJhKBpH1bzKydmXWoWgfOIWiMBngKmBCuTwCeTFYc0oLs3gKzLoXcLnDxLMjJZc7iIoZP\ne4F+U+YyfNoL6tYpUotEup3OAl4BBpjZejObaGbjzGw9cCow18zmh/v2MrN54aE9gH+Y2VvA68Bc\nd/9ruG0acLaZrQG+HD4XabiqHkV7tsDFf4IOPdTXX6Se6qwycvdLatg0u5p9NwCjw/V1wIk1vOZW\n4KzEwxSpw4Kfw/t/h3F/hF5DgNr7+quqRuRg6VvBKpKotx6BV/8AJ38bTrz4s2L19RepHyUEyWzF\nq+CZq+GI4XDOTQdsUl9/kfpRQpC0Uq9G4NJP4dEJwVDWF94L2TkHbFZff5H60WinkjbqNYqpOzzz\nQ9jyTjCUdYfDDno99fUXqR8lBEkb9WoEfvN+WPZocONZ/9NrfE319RdJnKqMJG0k3Ahc/DY8OwX6\nnwmnXZOCyERaBiUESRsJNQKX7YXHJ0LrdkEX06zsao8RkfpTQpC0kVAj8N9uhE3LYewd0EFjIoo0\nJbUhSNqosxH4nefgtTuC+w0+d06EkYo0T0oIklZqbATevQWe/C4cehx8+b9TH5hIC6CEIOnPHZ7+\nAezdAZfNgZy2UUck0iypDUHS39JH4O1ngi6mhx0fdTQizZYSgqS3Heth3mQ4/BT4wn9GHY1Is6aE\nIOnLHeZ8FyorYNwd6mIqkmRqQ5D0teg+eO9lOO8W6NI/6mhEmr1EJsi518yKzWx5TNlFZrbCzCrN\nrKCG4w43sxfNbGW47w9itt1oZkVmtiRcRjfNx5FmY/uH8NzPoN+XoOCKqKMRaRESqTKaCYyKK1sO\njCeYH7km5cA17j4QOAX4npkNjNl+q7sPDpd51b+EtEju8NT3g8cxvwOzqCMSaRESmTFtoZn1jStb\nBWC1fFHdfSOwMVzfZWargHxgZcPDlRZh8YOw7kUY/Wvo3DfqaERajJQ0KocJZQjwWkzxVWa2NKyS\n6pyKOCQD7NwI838KR3wRCiZGHY1Ii5L0hGBm7YHHgavdfWdYfAdwJDCY4CriN7UcP8nMCs2scPPm\nzckOV6L27GSoKIUxv4UsdYITSaWkfuPMLIcgGTzk7k9Ulbv7JnevcPdK4C5gWE2v4e4z3L3A3Qu6\nd++ezHAlaqueDpbTr4OuR0YdjUiLk7SEYEEDwz3AKne/JW5bz5in4wgaqaUlK9kOc6+FHoN0A5pI\nRBLpdjoLeAUYYGbrzWyimY0zs/XAqcBcM5sf7tvLzKp6DA0HLgNGVNO99FdmtszMlgJnAj9s6g8m\nGeZvN8LuYhjzfwfNjSwiqZFIL6NLatg0u5p9NwCjw/V/ANV2Q3L3y+oRozR3H74a3IR2yvcg/6So\noxFpsdRqJ9GqKINnfggde8OZP446GpEWTUNXSLRe+T0Ur4SLZ0Gb9lFHI9Ki6QpBovPJ+/DSzXDM\n+XCMRi8RiZoSgkTDPRjWOisbzr056mhEBCUEicqqp2HNc0G7QafeUUcjIighSBRKP4W/ToEex8Ow\nK6OORkRCalSW1Fv4K9hZBBfeB9n6LyiSLnSFIKlVvApeuR2GfB36nBx1NCISQwlBUscd5l4DbTrA\nl38RdTQiEkfX65I6Sx+FD/4J598G7bpGHY2IxNEVgqTG3h3w3E+DoSmGTog6GhGphhKCpMTaR39K\n5e7NjFk3juG/eok5i4uiDklE4ighSNK98NIL9H33T/y5fARLvT9F20u4/ollSgoiaUYJQZLLna4v\n/5idHML08n//rLikrILp81dHGJiIxFNCkORa+ign+ipuLr+EHRw4eN2G7SURBSUi1VFCkOTZuwMW\n/IwV9jkerTj9oM298nIjCEpEapJQQjCze82s2MyWx5RdZGYrzKzSzApqOXaUma02s7VmNiWmvJ+Z\nvRaWP2JmrRv3USTtvDQNPi1m82k30TbnwFnQcnOymTxyQESBiUh1Er1CmAmMiitbDowHFtZ0kJll\nA7cD5wIDgUvMbGC4+WbgVnc/CvgEmJh42JL2Nq2A1/4IJ13OGSNGMnX8IPLzcjEgPy+XqeMHMXZI\nftRRikiMhG5Mc/eFZtY3rmwVgFm1s2RWGQasdfd14b4PAxeY2SpgBHBpuN/9wI3AHYmHLmmramjr\nth3hrJ8DMHZIvhKASJpLdhtCPvBRzPP1YVlXYLu7l8eVS3Ow7C/BHcln3QCHdIk6GhFJUNo3KpvZ\nJDMrNLPCzZs3Rx2O1GXvzuCO5F5DYOg3oo5GROoh2QmhCDg85nnvsGwrkGdmreLKD+LuM9y9wN0L\nunfvntRgpQm8fDN8ugnO+00wG5qIZIxkJ4Q3gKPDHkWtgYuBp9zdgReBC8P9JgBPJjkWSbZNK+HV\nO+CkCcGYRSKSURLtdjoLeAUYYGbrzWyimY0zs/XAqcBcM5sf7tvLzOYBhG0EVwHzgVXAo+6+InzZ\n64D/MrO1BG0K9zTlB5MUc4d514YNyTdEHY2INECivYwuqWHT7Gr23QCMjnk+D5hXzX7rCHohSXOw\n7LH9Q1urIVkkI6V9o7JkgL07YP5PoNdQNSSLZDBNkCON98IvYfdm+NqjakgWyWC6QpDG2bAE3rgL\nPv+toKupiGQsJQRpuMpKtj32n2zzDpzw92EMn/aC5jgQyWBKCNJgi5/8P7p8spRf7PsaO2mniW9E\nMpwSgjTMp8Uc+dZ0XqkYyJzK4Z8Va+IbkcylhCANM//HtPFSflJ+BXDgAIea+EYkMykhSP2tfR6W\nPcafWn2Fdd7roM2a+EYkMykhSP2UlcDca6DrUXQ/dwq5OQd2M9XENyKZS/chSP0snA6fvAcTnmZM\nv/5UZrdh+vzVbNheQq+8XCaPHKB5D0QylBKCJO7jZfDP/4MTL4V+XwI08Y1Ic6IqI0lMRTk8eRXk\ndoaRv4w6GhFJAl0hSGJevR02LoGLZmrwOpFmSlcIUret78KL/wsDzoOBY6OORkSSRAlBaldZCU99\nH7Jbw3m/BrO6jxGRjFRnQjCze82s2MyWx5R1MbMFZrYmfOxczXFnmtmSmGWvmY0Nt800s/ditg1u\n2o8lTeb1P8IH/4BRU6HjwfcciEjzkcgVwkxgVFzZFOB5dz8aeD58fgB3f9HdB7v7YGAEsAd4LmaX\nyVXb3X1Jg6KX5NqyFv7233D0OTD4a1FHIyJJVmdCcPeFwLa44guA+8P1+4G6KpYvBJ519z31jlCi\nUVkBc74DrdrAv/1WVUUiLUBD2xB6uPvGcP1joEcd+18MzIor+6WZLTWzW82sTQPjkGT51+9g/esw\nejp07Bl1NCKSAo1uVHZ3B7ym7WbWExgEzI8pvh44Bvg80AW4rpbjJ5lZoZkVbt68ubHhSiI2vgUv\n3ATH/hsMuijqaEQkRRqaEDaFP/RVP/jFtez7VWC2u5dVFbj7Rg+UAvcBw2o62N1nuHuBuxd07969\ngeFKwvbtgce/Be26qapIpIVpaEJ4CpgQrk8Anqxl30uIqy6KSSZG0P6wvJrjJArzfwxb1sC4O3UD\nmkgLk0i301nAK8AAM1tvZhOBacDZZrYG+HL4HDMrMLO7Y47tCxwOvBz3sg+Z2TJgGdANuKnxH0Ua\n7e25sOg++MJ/Qv8zoo5GRFLMgiaAzFBQUOCFhYVRh5GW5iwuatyoo9s/hDtPg85HwMS/QavWyQtW\nRFLKzBa5e0Fd+2kso2ZgzuIirn9iGSVlFQCfzW0MJJYUykvh0QngHoxVpGQg0iJp6IpmYPr81Z8l\ngyr1mtt4/k9gw5sw9g/QpX8SIhSRTKArhGagpjmMqyuPr1q67bi1fH7RXUG7wbHnJztUEUljukJo\nBmqawzi+vKpqqWh7CQ7k7VjFcYU/ZUuXoXDWDSmIVETSmRJCMzB55ICE5jaOrVrqznbuav1rPqE9\n39x9FWTnpCxeEUlPqjJKc4n0Hqp6Xtd+VVVIbdjHjNa3kMduLtp3AytL26bmw4hIWlNCSGP16T2U\nyNzGvfJyKdq+h2k5dzEkay3f3nc1K7wv+TVUOYlIy6IqozTW6N5DcSaPHMDPWj/MuOx/8uuyi/hr\n5bBqq5ZEpGXSFUIaq0/voUSMLXkCsp7miexR3L53LPkNuYFNRJotJYQ0FlTxHPzjX1OvolotmQXP\n/RQGjmX8hfcyPiu77mNEpEVRlVEaS7T3UJ2WPw5Pfg/6nQ7jZ4CSgYhUQ1cIaSzR3kO1WvoYzJ4E\nfU6Fi/8czIAmIlINJYQ0l0jvoRq99XAwDeYRw+HSR6B1u6YNTkSaFVUZNVev/AFmfxv6fhEufVTJ\nQETqpCuE5qayIhis7rU7gikwx98FObrPQETqpoTQnJR+CrOvhLefgVO+C+fcpAZkEUlYQlVGZnav\nmRWb2fKYsi5mtsDM1oSPnWs4tsLMloTLUzHl/czsNTNba2aPmJkG4W+MzavhrhGweh6MnAqjpioZ\niEi9JNqGMBMYFVc2BXje3Y8Gng+fV6fE3QeHy5iY8puBW939KOATYGLiYcsBlv0FZpwJJdvgsjlw\n6ncb/ZJzFhcxfNoL9Jsyl+HTXmDO4qImCFRE0llCCcHdFwLb4oovAO4P1+8Hxib6pmZmwAjgLw05\nXkJ7tsHj34LHJ0LPE+DKhdD/9Ea/bPww2VVjKCkpiDRvjell1MPdN4brHwM9ativrZkVmtmrZlb1\no98V2O7u5eHz9UC1fSvNbFJ4fOHmzZsbEW4zs/IpuH0YrJgNp0+BCU9Dx15N8tJNPYaSiGSGJmlU\ndnc3M69h8xHuXmRm/YEXzGwZsKMerz0DmAFQUFBQ03u0HMWrgiEo1v4Nep4YVBEddnyTvkVTj6Ek\nIpmhMQlhk5n1dPeNZtYTKK5uJ3cvCh/XmdlLwBDgcSDPzFqFVwm9AdVH1GbnBlg4HRbNhDYd4Jxf\nwslXJmVimyYdQ0lEMkZjqoyeAiaE6xOAJ+N3MLPOZtYmXO8GDAdWursDLwIX1na8AFvWwJNXwW0n\nwJsPwLBJ8P0l8IXkzXLWZGMoiUhGSegKwcxmAWcA3cxsPXADMA141MwmAh8AXw33LQC+7e7fAo4F\n/mhmlQTJZ5q7rwxf9jrgYTO7CVgM3NNknyrTle8Luo8ufhDWPh+MP3TS5UES6Nw36W/fJGMoiUjG\nseCP9cxQUFDghYWFUYeRHOWl8N7fYfVcWPkk7NkKHfNhyNfh8/8B7btHHaGIZCgzW+TuBXXtpzuV\no1JZAZuWw/v/hA/+Cetehn27IKcdHH12kAiOHKGby0QkZZQQkq2yAnZ8BNvWwdZ3g15CHy+FTSug\nbE+wT+e+cPx4OOa8YM6CHE16LyKp1zISwu6twV/f1TmgysyD5+7hemWwVFZAZRlUlAeP5aXhUgL7\n9sC+T4Ol5BPY80lwx/Cuj4Pl003gMX3623SCwwbB0AmQPzQYmrqT6uZFJHotIyG89L/wxt3Jf59W\nbSG3CxzSBdofCocOhI49odPh0PVI6HIkdDgMzJIfi4hIPbWMhHDCxZB/Ui07xPxAW1b4g22QlRU+\nz4KsnKCGSBDEAAAGeUlEQVSbZ1Z28MPfqm3Q+6d1O2jdIXhUVY+IZLCWkRAO/3ywiIhIjTRjmoiI\nAEoIIiISUkIQERFACUFEREJKCCIiAighiIhISAlBRESAlnIfQhqas7hIw0uLSFpRQohA1ST2VfMW\nV01iDygpiEhk6qwyMrN7zazYzJbHlHUxswVmtiZ87FzNcYPN7BUzW2FmS83s32O2zTSz98xsSbgM\nbrqPlP40ib2IpKNE2hBmAqPiyqYAz7v70cDz4fN4e4BvuPtx4fG3mVlezPbJ7j44XJbUP/TMpUns\nRSQd1ZkQ3H0hsC2u+ALg/nD9fmBsNce94+5rwvUNQDGgab+oebJ6TWIvIlFqaC+jHu6+MVz/GOhR\n285mNgxoDbwbU/zLsCrpVjNr08A4MpImsReRdNTobqfuVbPJVM/MegIPAt9098qw+HrgGODzQBfg\nulqOn2RmhWZWuHnz5saGmxbGDsln6vhB5OflYkB+Xi5Txw9Sg7KIRKqhvYw2mVlPd98Y/uAXV7eT\nmXUE5gI/cfdXq8pjri5Kzew+4Nqa3sjdZwAzAAoKCmpMPJlm7JB8JQARSSsNvUJ4CpgQrk8Anozf\nwcxaA7OBB9z9L3HbeoaPRtD+sDz+eBERSa1Eup3OAl4BBpjZejObCEwDzjazNcCXw+eYWYGZVc1V\n+VXgS8Dl1XQvfcjMlgHLgG7ATU36qUREpN7MPXNqYQoKCrywsDDqMEREMoqZLXL3grr201hGIiIC\nKCGIiEgoo6qMzGwz8EEDD+8GbGnCcJqK4qofxVU/iqt+0jUuaFxsR7h7nTcGZ1RCaAwzK0ykDi3V\nFFf9KK76UVz1k65xQWpiU5WRiIgASggiIhJqSQlhRtQB1EBx1Y/iqh/FVT/pGhekILYW04YgIiK1\na0lXCCIiUotmlxDMbJSZrTaztWZ20MQ9ZtbGzB4Jt79mZn3TJK7LzWxzzDAf30pBTAfNhhe33czs\nt2HMS81saLJjSjCuM8xsR8y5+nmK4jrczF40s5XhTIA/qGaflJ+zBONK+Tkzs7Zm9rqZvRXG9d/V\n7JPy72OCcaX8+xjz3tlmttjMnqlmW3LPl7s3mwXIJphzoT/B/AtvAQPj9vkucGe4fjHwSJrEdTnw\n+xSfry8BQ4HlNWwfDTwLGHAK8FqaxHUG8EwE/796AkPD9Q7AO9X8O6b8nCUYV8rPWXgO2ofrOcBr\nwClx+0TxfUwkrpR/H2Pe+7+AP1f375Xs89XcrhCGAWvdfZ277wMeJpjdLVbsbG9/Ac4KR12NOq6U\n8+pnw4t1AcFote7B8OV5VSPVRhxXJNx9o7u/Ga7vAlYB8WOYp/ycJRhXyoXn4NPwaU64xDdapvz7\nmGBckTCz3sB5wN017JLU89XcEkI+8FHM8/Uc/MX4bB93Lwd2AF3TIC6Ar4TVDH8xs8OTHFMiEo07\nCqeGl/zPmtlxqX7z8FJ9CMFfl7EiPWe1xAURnLOw+mMJwZwpC9y9xvOVwu9jInFBNN/H24AfAZU1\nbE/q+WpuCSGTPQ30dfcTgAXs/ytADvYmwa34JwK/A+ak8s3NrD3wOHC1u+9M5XvXpo64Ijln7l7h\n7oOB3sAwMzs+Fe9blwTiSvn30czOB4rdfVGy36smzS0hFAGxmbx3WFbtPmbWCugEbI06Lnff6u6l\n4dO7gZOSHFMiEjmfKefuO6su+d19HpBjZt1S8d5mlkPwo/uQuz9RzS6RnLO64orynIXvuR14ERgV\ntymK72OdcUX0fRwOjDGz9wmqlUeY2Z/i9knq+WpuCeEN4Ggz62fBjG0XE8zuFit2trcLgRc8bKGJ\nMq64euYxBPXAUXsK+EbYc+YUYIfvn/40MmZ2WFW9qZkNI/h/nPQfkfA97wFWufstNeyW8nOWSFxR\nnDMz625meeF6LnA28Hbcbin/PiYSVxTfR3e/3t17u3tfgt+IF9z963G7JfV8NXRO5bTk7uVmdhUw\nn6Bnz73uvsLMfgEUuvtTBF+cB81sLUHD5cVpEtf3zWwMUB7GdXmy47JgNrwzgG5mth64gaCBDXe/\nE5hH0GtmLbAH+GayY0owrguB75hZOVACXJyCpA7BX3CXAcvC+meAHwN9YmKL4pwlElcU56wncL+Z\nZRMkoEfd/Zmov48JxpXy72NNUnm+dKeyiIgAza/KSEREGkgJQUREACUEEREJKSGIiAighCAiIiEl\nBBERAZQQREQkpIQgIiIA/D9q3jlJLz+b/QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11091a5f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, 4, 100)\n",
    "y = logistic4(x, *popt)\n",
    "plt.plot(xdata, ydata, 'o')\n",
    "plt.plot(x, y)\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}