{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Minimization</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import chi2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b> Minimization, not using least squares</b><br><br>\n",
    "\n",
    "Let's try a general minimization function, for Poisson uncertainties on a Gaussian profile. First write out obligatory routine to return a Gaussian. Here, we'll add another constant background/continuum value. We will also add an option to return the vector of derivatives of the Gaussian with respect to each parameter.\n",
    "<br>\n",
    "if \n",
    "$$f = A \\exp{-0.5(x-x_0)/\\sigma^2} + c$$\n",
    "then the derivatives are:\n",
    "$${df\\over dA} = \\exp{-0.5(x-x_0)^2/\\sigma^2}$$\n",
    "$${df\\over dx_0} = f {(x-x_0)\\over \\sigma^2}$$\n",
    "$${df\\over d\\sigma} = f {(x-x_0)^2\\over \\sigma^3}$$\n",
    "$${df\\over dc } = 1$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12f8b5390>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1fn48c/JSkJISCBhT4JCtVY0SkDcU7GKfN3w64IdEQSNa39Ctf0iWGtrU7V1Qa2KQRC0I4o7WBURG2y1ARGBIGJBSULY1yQQyDI5vz9mgkmYLbPcuXPzvF+vvMjcuTP34c69T+6c+5xzlNYaIYQQ0S8m0gEIIYQIDUnoQghhEZLQhRDCIiShCyGERUhCF0IIi4gzcmM9e/bUubm5Rm5SCCGi3ldffbVHa53paz1DE3pubi4rV640cpNCCBH1lFIV/qwnTS5CCGERktCFEMIiJKELIYRFSEIXQgiLkIQuhBAWIQldCCEsQhK6EEJYhCR0IYSwCEnoIjoVFDh/hBBHSUIXQgiLkIQuhBAWIQldCCEsQhK6EEJYhCR0IYSwCEnoQghhEZLQhRDCIiShCyGERUhCF0IIi5CELkR70gtVRCmfCV0pNUAp9U+l1LdKqW+UUne7lmcopZYopTa6/k0Pf7hCCCE88ecKvQm4R2v9U2AEcKdS6iRgKrBUaz0YWOp6LIQQIkJ8JnSt9Xat9SrX77XAt0A/4Apgnmu1ecCV4QpSCCGEbx1qQ1dK5QKnAcuBXlrr7eBM+kCWh9cUKqVWKqVW7t69O7hohRBCeOR3QldKpQBvAZO11jX+vk5rXay1ztda52dmZgYSoxBCCD/4ldCVUvE4k7lda/22a/FOpVQf1/N9gF3hCVEIIYQ//KlyUcBs4Fut9ROtnloIjHf9Ph54L/ThCSGE8FecH+ucDYwDypRSq13LpgGPAAuUUpOASuCa8IQohBDCHz4Tutb634Dy8PTI0IYjhBAiUNJTVFhSwdwCCuYWRDoMIQwlCV0IISxCEroQQliEJHQRfex2KC2FZcsgN9f5WAghCV1EGbsdCguhvt75uKLC+ViSuhCS0EWUmT4d6uraLqurcy4XopOThC6iS2Vlx5YL0YlIQhfRJTu7Y8uF6EQkoYvoUlQEycltlyUnO5cL0clJQhfRxWaD4mJITHQ+zslxPrbZIhuXECbgz1guQpiLzQazZjl/LymJaChCmIlcoQshhEVIQhdCCIuQhC6EEBYhCV0IISxCEroQrck4MSKKSUIXooWMEyOinCR0IVrIODEiyklCF6KFjBMjopwkdCFayDgxIspJQheihYwTI6KcdP0X1rR6dcdf0zIezKRJzhujOTnOZC7jxIgoIQldiNZknBgRxaTJRQghLEISujCnggLnTwDsZXZKU2tY1r2a3Bm52Mukjlx0DtLkIizFXmancFEh9bEagIrqCgoXFQJgGxKatvCCuQUAlEwoCcn7CREqcoUuLGX60unUNbbtHFTXWMf0pdI5SFifJHRhKZXV7jsBeVoeFkE0FwkRDEnowlKy09x3AvK0XAgrkYQuLKVoZBHJ8W07ByXHJ1M0UjoHCeuThC4sxTbERvFlxSQ6FGjIScuh+LLikN0QFcLMpMpFWI5tiI1Zs+8EoOTB8sgGI4SBJKGL6CS9OIU4hs8mF6XUHKXULqXUulbLHlRKbVVKrXb9jA5vmEKETsHcgqO15EJYiT9t6HOBUW6WP6m1znP9fBDasIQQQnSUz4Sutf4M2GdALEIIIYIQTJXLXUqpta4mmXRPKymlCpVSK5VSK3fv3h3E5oQQQngTaEJ/HjgeyAO2A497WlFrXay1ztda52dmZga4OSEESPu/8C6ghK613qm1dmitm4FZwPDQhiWEEKKjAkroSqk+rR6OAdZ5WlcIIYQxfNahK6XmAwVAT6VUFfB7oEAplQdooBy4NYwxCiGE8IPPhK61vt7N4tlhiEWI6BDIfKVCGEDGchFCCIuQhC6EwaRSRYSLJHQhhLAISejCeux2XptTw6dPVUNuLtg7OEl0SYkM/iWikiR0YT52O5SWwrJlHU/IdjsUFtK7VjsP7ooKKCzseFIXIgpJQhfm4krI1Nc7H3c0IU+fDnVtJ4mmrs65XAiLk4QuzCXYhFzpYTJoT8uFsBBJ6MJcgk3I2R4mg/a0vIPsZXZKU2tY1r2a3Bm52MuiqylHKmysTRK6MJdgE3JRESS3nSSa5GTn8iDZy+wULiqkPlaDgorqCgoXFbZN6sG0/wsRJEnowlyCTcg2GxQXs6ObohkgJweKi53LgzR96XTqGts2B9U11jF9qas5KNj2fyGCJAldmIsrIZOY6HwcSEK22Rg7MZUL7k6D8vKQJHOAymr3zT5Hl4fghqw0iYhgyCTRwnxsNpg1y/m7ierBs9OyqaiucLsckBuyIuLkCl1ERLivREtW51GyOi+k71k0sojk+LbNQcnxyRSNdDUHhfmGrBC+SEIXwk+2ITaKLysm0aFAQ05aDsWXFWMb4mrSCeMNWSH8IU0uQnSAbYiNWbPvBKDkwfJ2T7oS+6RJzhujOTnOZB6iNnwhfJGELkQombT9X3QO0uQiOhV7mZ3SqlKWVSyLuo5B0Ry7MIYkdNFpHO0Y5HDWibvtGGRS0Ry7MI4kdBEeBQXOHxPx2THIxKI5dmEcSeii0/DZMcjEojl2YRxJ6KLTONoByM/lZhLNsQvjSEIXnYbPjkEmFs2xC+NIQhedxtGOQbHOcWKO6RhkYtEcuzCO1KGLTsU2xMasr5x14iUTSjr+Bq75SrNqNbyba2jHoaBjF5YnV+hC+CsE85VGspZc6titTxK6EP4KcnjcSNaSSx175yBNLqLTa3I0s/XAYX7Yc4jyPYfYvOcQNYcbyc5IZmBmVwb2TGFgz66kBTk8rrda8nC3hUdy28I4ktBFp1N/JItDNSdy87yVbN5zkMp9dTQ69NHnuyXGkZoUz8I122j+cTFfpGbSt3rXsW/o5/C4kawllzr2zkESuugUmhzNLP5mJ3O/2Mz28vGgmtgSX8fgrG5c9LPeDOzZ9ehPj64JKKWob3KwZd9hNu85xOY9B/nsyGSufP4PdGmoP/q+DYldqJ32e3r4EYPPCTLCKJLbFsaRhC4sbd+hBuavqOTvpRVsrz5CdkYy6Y2vkOL4jMVTtnl9bWJcLIOyUhiUlQL0gvPug2HZ7Lh9HFm1mv09e/PQWTfw3uYsLnx5JTedlcuZx/dAeXi/opFFFC4qbNP0YVQteSS3LYwjCV1Y0rfJmbzU+3TefXgpDU3NnDOoJw9dcTI/PzGLkb8eE/gb22yM/dI1HvqM7fxf9WH6l1by6opKlqzfyU96pTAh6xSu2r2eLu1f6mqrnvTGOOpjNDndcygaWWRIG/bRbb83iXpHPTlpxm1bGEcSurCUww0O/rJ4Ay+dMoEkRwPXDO3PhLNyGdyrW1i21yctiXsvPoG7LhjEojXbmPtFOdOOu5gX+g7n8fJ95OdmtFnf6wQZYSZ17Nbns2xRKTVHKbVLKbWu1bIMpdQSpdRG17/p4Q1TWIlf9dAlJR2eIGJV5X7+5+l/8dLn5YzfsYrSVTMpGjMkbMm8tS7xsVyTP4D3f3UOf1+/gGYU17zwH4r+sZ4jjQ6/3kPqxEWw/KlDnwuMardsKrBUaz0YWOp6LIRP4aiHrm9y8OhHG7j6+S+ob2rm1ZvP4A9zf0fa0sUBvV8wE1grpbj/uPdQ/JpfDs9m1r82c+kz/2bNlgNeXyd14iIUfCZ0rfVnwL52i68A5rl+nwdcGeK4hEWFelzvdVurufyZz3m+5HuuzR/AR5PP5axBPUMRalBiqKdozBBenjicQ/VNXPX8Fzz+8Xc0NDW7Xd+w8c5Xr3b+CEsKtA29l9Z6O4DWertSKsvTikqpQqAQINvPel1hXaGqh250NPPsPzfxt083kdE1gZcmDOPnJ3o8DCPmvJ9k8tHk8/jjovU88+kmPvl2Fw1qAAl6S5v1pE5chELYu/5rrYu11vla6/zMzMxwb06YXCjG9d5Ve4RrZv6HGZ9s5NJT+vDxlPNMmcxbpCXF8/i1pzLrxnx219azLeFP1MZe0GYdGe9chEKgCX2nUqoPgOtfN93nRKdlt0NpKSxbBrm5bQavCnZc7027ahnz7Bd8t6OW52ynM2PsaXRPTghl9GHzi5N6sWTKeSQ1l7E3fhKPfrSBZldXVBnvXIRCoE0uC4HxwCOuf98LWUQiurlGJKTe1ZuyZURCAJstqHro0h/2UvjyShLjY1lw65kM6Z8Wrv+FVyWr87yvkOf5+fSuCWQ1PsE+fSPPl0DV/sM8ds0pIasTb7mZK2WJnZPPhK6Umg8UAD2VUlXA73Em8gVKqUlAJXBNOIMUUcTbiISuccMDqYd+b/VWfvPGWrJ7JDP3pmH0T0/2/SKTUjST0TSXOy/7FY98uIGdNUcoHjdU6sRF0HwmdK319R6eGhniWIQVBDkiYXtaa54r+Z6/Lv6OEcdl8MIN+aQlxwcRoDko4Lbzj6df9yTuWbCGq57/gnk3DY90WCLKyXjoIrQ8VTIFUOHU5Ghm2jtl/HXxd1yR15d5E4dbIpm3dtmpffn7zWew92ADY577nPrDvSMdkohiktBFaBUVQXK75pDkZOfyDjhY38SkeSuZv2ILd/18EDOuyyMxLjaEgZrH8IEZvHX7WSQlxLKj8jrqao/3uG7JhBLvzTFSZ96pSUIXoWWzQXExJDonMyYnx/m4A/NuVh9uxDarlH9v2sPDVw3h3otPQClPYxhaw6CsFN6+/WziE/eya+uVvL2qyv2KBQXOHyHckMG5ROjZbDDLeXOvo+Ox1Bxp5MbZy1m/vYaZNwzlFyf1Cn18JpXZLZHe2a+zq2oM97wRg1Iw5rT+kQ5LRBG5QhemUXOkkXGzV7B+ew3P2zpXMm8RE9NIVv+3GTGwB/csWMN7q7dGOiQRReQKXZiC88p8Beu3VfPsL0/nwjAmc5915BEWE9PE7An5TJz7JVNed7aHX5HXL8JRiWggV+giLAryVlOQ59/NudojjYyfs4J1W53J/KKfdeJKD9dNzeSEOOZMGMbwgRlMeX21XKkLv0hCFwErmNydgsndg3qP2iON3DhnBWVV1Txr6wTJ3G7ntTk1fPpU9THDIrTXktSH5TqT+sI13qfM80fJ6jzTf0MRgZOELiKm2ZHAeFcy/9svT+fiTpDMKSykd612nngtwyL4SOov3TSM/NwMJr/2NYt6nGBYuCL6SEIXEdHsiGdn1dWsrarmb788jVEnWzyZg/dhEbxITojjpQnDyM/JYPKgS3k/Q5K6cE8SujDckUYHO6uuov5wH565/jRGndwn0iEZI4hhEbomOq/Uh9Zu5e7Bl/LJ+p0hDk5YgSR0YahGRzN32FdRf3gAPfv+g0uGdJJkDkEPi9A1MY45G97i5EM7uOPVVXzx/Z6Obd/LsMbCGiShC8M4mjX3LFjDpxt2kdFrCSmpG4wPIpJJLQTDIqQ0NzJ3w1vk9kjmlnkrfc5VepSnYY0lqVuKJHRhCK01D7y3joVrtvHbUSeQmr7G+CBCldRKSjrcAxYIybAIAOlNR3hl0hlkpCQw/qUV/Hdnre8XBdh+L6KLJPQoFszs9Eb76+LvsC+v5Nbzj+OOgkGRCcIMSc1mgxEj4Pzzoby8w8m8Ra/ULtgnjSAhNoZxs5ezZV+d9xeEeFhjYU6S0EXYzVz2Pc+VfM/1w7OZOurEyAVisaSW3SOZVyadwZHGZmwvLmdXzREvK4duWGNhXpLQRVi9urySRz7cwKWn9OFPV54c2VET/Uhq9jI7pVWlLKtYRu6MXOxlHWuOCfb1HXVC727MmzicvQfrGfl8EaWpB1nWvfrYbYdoWGNhbpLQhUfBNukcihnB9HfLKDghkyeuzSM2JsJD4PpIavYyO4WLCql3ONvYK6orKFxU6HdSDvb1gcob0J1rzi3nm8N/pT7WAcrNtkPUfi/MTRK6CIu6mFPZHX87w3IyeN42lIS4toeaz4kawsFHUpu+dDp1jW3bousa65i+1L829mBfH4y53zyCVvXetx2i9nthXjLaogi5FZv3sTt+Mgl6Cy9OGE1SgolmGvIyVntltfu2dE/L/V3P39dD4CNBhmLbIvrJFboIqXVbq5k090ti9R56NTxKapfomQM0O819G7un5aF+fTA8bWNA2oCwb1uYhyR0ETKbdh3kxjkrSE2Kp/fghcTmeZ4b04yKRhaRHN+2jT05Ppmikf7dOAz29cEoShxNcmPbZUkNcOX2M9Fah337whwkoYuQqNpfx7jZy4lR8PebzyAu3o/OLiZjG2Kj+LJiEmOdbew5aTkUX1aMbYh/bc3Bvt4nL71cbY9+QPFCyDkASjv/nbUIfjPzU578ZGNoti9MT9rQTa6lysTwG4gdsLu2nnGzV3CwvonXC89kYM+ukQ4pYLYhNmZ95WxjD2SfB/t6jzz1cgXnfYHKSmwabGVtX6bVHp5eupG0pHgmnTPQ93ZaJqAOpCesiDi5QheeuWbPcavVRA3NOTkM/fwD5t40jJP6phobY2fhq5erpxr7AQO45OTePPT+ehas3BLeGEXESUKPUkZ3YGm7cTv2J29ixCRN3O/hjPE7ObViBkP//WHkY7MqX71cPdTYqz//mRlj8zh3cE+mvrWWDzN+Et44RURJQo9CfndgKSj48St0KLf/4t0UXtxIRXfQCiq6w22XNGJ/8e6Ida6xPF+9XF019ju6KZqhTY19YlwsL4wbymnZ6fy/QZfyWVquQUELo0lCj0KR7MACMC1vL3UJbZfVJcD0vL2GxRZNA5OFhD9d9202xk5M5YK7047pOJScEMec8cMYdHgvt5xwJV9s6uBY6iIqSEKPQpHsRNLoaGZLmvvnKtOkg0tQvI3VHoKu+2nJ8fz92wXkHDnApHkrKf1hb8fiC9M3PhE6ktCjUKQ6sDQ5mrn7ta9JaejmfvvxPSLauSaq+TNWewi67vdoOoz92wX0S09i4twv+bJ8X2jiF6YgCT0KRaIDS5OjmSkL1vBB2Q6uPe0PJKu2bS7JKoGiy5+KaOeaqGbgWO2ZjXW8essZ9E7rwoQ5K/iqYn/ItyEiQxJ6FDragcWhQIehA0s7jmbNvW+sYdGabdx3yYm8eM0UisfMabv9MXOwDbGFv3ONVRk8VntWty7Mv2UEmd0SmTBnBav9ncpOmJok9ChlG2JjRE0q5x9Io3xyedgSpkbx2zfX8u7qbfzm4hO49fzjfW7fNsTGiP4jOD/nfM+xSXtsWxGYgKJXahfmF44gvWsC42Yvp6yqOmzbEsYIKqErpcqVUmVKqdVKqZWhCko4RbKe215mpzS1hs+6H+Cpb37Beadu4M6fR2jquM4gQhNQ9ElLYn7hCNKS4rlh9nK+Sc4K6/ZEeIXiCv3nWus8rXV+CN5LuESynvvotmM1KHDE7OatzQ9ILXk4GTUBhZsJrvt1T2L+LSPomhDLDT+9xpJJvbOUuUqTi0lFstZ82ifTIlrn3mlFcAKKARnJzC8cQVJzE2NPuo4Vm6X6JRoFm9A18LFS6iulVKG7FZRShUqplUqplbt37w5yc51HpOq56xqaqKxxP+aHkbXklr6icnOVbAY5Pbryxjevktl4iHGzl/Pphp2RDkl0ULAJ/Wyt9enAJcCdSqnz2q+gtS7WWudrrfMzMzOD3FznEYl67gN1Ddzw4nJim3savm1hDv0aannjm/n8pFc3bnn5K975usrv11r6j3CUCCqha623uf7dBbwDDA9FUML4WvOdNUe47oVS1m2t4ddn/F5qyTuxHk2HmV84guG5GUx5fQ0vfb7ZmA1L5VPQAk7oSqmuSqluLb8DFwHrQhVYZxd0rbm3buTtlO85xNUzv6Bqfx0v3TSMv4y+k+L08fSvdk2WcDCW4vTxUktuda2OmZQTBjGvy0YuOqkXf1i0nieW/BeZ9ygIBv2xCmaCi17AO0qplvd5VWv9UUiiEoBrsoTZdwJQ8mC5/y/0NRlCK+u31XDjnBU4mpt59ZYRnDqgO9jt2O6dh+3ofVEHJM+D1LNlpnircnPMJNx+G8/PfIH78ofw9NKN7M+9kD+UfxK2SoqCPOfY+yUeVyhw/mvC+w9mEXBC11r/AJwawlhEqHjrRt4qIa+4zMakE64ipUd3Xis8i0FZKR16vbAQD5957O/u59HNm0lPTuCFz+BAXBcea3KQGBcbmTgD0NKfo95RT+6MXIpGFln226ZMQRfFSlbnuX/Cj27k81dU8sBPr2NA/QFeuX00/bondej1wpw8HhO+ePnMlVLcN/qnpL/yEo/knM/W4lJm3jCUrNQugQdqEE/9OYC2Sd0iV/9Shx6tvLWRe+lG3tDUzLR3yrjv7TLOrKnknXX2tsncx+stwaRlgxHl6zO327ntnafZ/OilPDPtKp6Z9AdWVbYb1MvblIUREum5A4wmCd3MWs3b2SZp+xpq1UM38urfPcgvZ5Xy6vJKbjv/eF7a8BZpriuXNiLUDV1EkLfPvNXxpoB+1buY9t6T/P1XD/P6lx341haBKpbONj6/JHSzcp1EvWu180NqnbR9DbXqpht5+cMzuHjnANZtq+aZ609j6iUnEuupbsHLdGbC5AL99uHtM3dzvCU11jPt36/wf2+V8bt319HQ1ByC4EPPFOPzd6DiLFjShm5W3pK2P23cNhvMmgXAgsde4f5315HVTfH27WdzUt9U39u32Rj7pavCZkZ5AP+BzqtkQkmkQwiMp8/cw/HWY/9OCs87juLPfuC7HbU4SCWWGgMC9V/RyCIKFxW2aXYxtE9FByrOQiE6rtB9fVWL5g4JnmL3lrT9bONuVDE8mHMBv31zLfk56Sy86xz/knknVzKhJHqTcjh4ON5UdjbTRv+Up8bmsXbrAbYl/ol6NdDj2xTkrT5amngMT82LQYr4+PwGTlwC0ZLQOyNvSduPNu5Nu2q5+me/ZG6fodx8zkBenjicjK7tZnYOUsmMA5TMkIkRLM/H8XZFXj/evO0sQLM94fc8X/I9TY4ONMF4a15svY6vZgsPF0c+x+cPZ5OIwRVjktDNyttJ5GWoVUez5oVl3zP66X9TmZjGs/99j/svPYm4WPmoLSESFTp+DO17cr80+tZPJ7n5ax79aANXz/wP3+8+6N/7+7qK9We+1UCF873B8IoxOct9iVRzjq+TyM1Qqz/sPsg1M7/g4Q83UPCTTD5e+xL/s++/xsdOeJstIjnxR6flY2hfe5mdL1O38W3XP1KXcSur9ixi9FP/4sV//YCj2cegAZWV2IdA7mSI+b3zX/sQfryKDWezRaje21OeMLhiTG6KBqlldLmwJK9WNza9XZU1N2vmflHOXxZvICE2hhnX5XFFXl/UnDqPr4lWfncUEYZpMyEKsPvwVpLinmF4z6786R/NLP5mB3+92nOncvv5GRSetZc6V4tgRXcovAzokYENwttsEe4mkZY/fJMmOb8F5OT8+C07DKL/Cj0E7V/ehv30eiPHBCoT0xg7q5Q/vr+eM4/rwZJfn8+Vp/XDNcaOOQXxmXW2jiJm4ulccPeZHG6qY3Pjizx+zals2FHLqKc+oyb2IjTHHpfTL+RoMm9Rl+BcDoS32cKIJhEDJy4xf0K329mx6jOa3Z384W7/MrG6hib+1ncEo06ZwLfbavjL1acwZ8Iwepm9O3aQn1ln6ygSDTzt+y3VW/jfof35eMp5nDGwB/vix7Mj4X6+btfDtLLJ/exIR5eHs9nCYp3ozJ3Qfd39NrgkyAyaHM28tqKSgr+W8Fj2uZxdXcHiKedxbf6AY67KTfntIsjPzBQdRUQbvj6TPmlJzL1pGD0ai2lUfRjz3BfcYf+KzXsO+fX6UHR083hPx6i5XA1i7oTu6+TvRINIaa1Zsn4no576F1PfLqN/ehJvrnuVWf99l77tx2IJkZLVeYEP9uRJkJ+Z0RN/CN+KRhaRrNq2mSSrhDafiVKKbo5l9K//NXePHEzJd7v5xRPLeOC9dUw9+w++P1ObjbETU7ng7rTAmi28FTdEcC7XUDN3Qvd18oeq/cvToEJh6uzQ+v39aUteldKHa1/4D7e8vJLmZs3MG4by1u1nkX9wa9i3HXJBfmYR7ygSblE4cJhtLRQv1OQccE2IcsD52La21Uquc6nkqZ1MueUi/pOzg+uGDcC+vJJnFvXiqtw/kuiICWwyF9f7R+R4jvS22zF3lUt2trOZxd1ycLZzFRa2vYoPVftXS3NPnavkKtRddn10CdZa8/WWAxQPvoKPevyEnnvq+NOVJ3PdsAHE+1tTnufh6trg7shthOAzsw2xMesrZ/WP9Og0genTsVU0Yvuq9cLGH8fPd3Mupf2/OykqLmbilCv4y0cbWLz2RPrpeaQ2vc+yqa+Q2iXe/+1H8nj2d9sG/ZE29xV6URH2ofFt61OHxv948oeg/cteZqc0tYZl3avb1jRPn479+Lq22z6+bVtvUPXQHpqTmqdN442VW7j8b59z1XNf8HlaDpO3fM6y3xRww4gc/5N5ANs25N6DH5+Z1JlHhsdzwRdf36S9nEvHZ6bwwrh83rr9TOL0DvbHj2PEn5dy/7tlbNxZ619skTyeTXYfz9RX6PZToPByxdE/7N2djzkFjp7+ftZqu33/dvWzrWuaSa2g8DKOrY1dVIGt9Wt91EN7rFP3eBJs4TdvrmVwVgoPXXkyV027ma7NjZAYwo8q0vcevHxmUmceGd7OBZ/73cc3abuPcwlgaE4GfRr+RL3K5eyh81iwsoq/l1Zy1vE9GNB/FX/7+jeeY4vk8Rzpc6kdU1+hT186nTrd0GZZnW4IWc2xt5rm6RfHuq+NvTjW52v94qHNeH+PXrx68xl8POU8xo3IoeunSzz/oQq0vdXEE1hInXlkBLXffZT++TqXWkvU5Tx2zan8Z+oF/HbUCZTvOcTjpQ96j82A49ljXxWTnUumTujhrjn29v6VKQ73z7mW+x1buxuuRxodfLJ+J69efiuH4xPbrNqclESPpx7jrEE9w9sxyN/a2wjcoJM688gIar/7KCv0dS650yMlkTsKBvHZb39Oc8we77GFopY80GPdz21767wYSqZO6OGuOfb2/tlpOR6ey+lwbP1MX14AAA6WSURBVM10YdGabdz56iqGPrSEm19eycMZp/PmrQ9wOCnFOc1ETg4xs2YZUzJl4gkspM48MoLe717KCn2dS97ExcZ4jCGmuSejZnzGk1nDqHr0Kc/HczirUEx2Lpk6oYe75rgocTTJjW2XJTc6l/vatrfnD9U38fmmPcz45L/siJ9KZeJMfjX/a5b/sJfL8/oxb+Jwvrr/F4x7ZhpJw4eiIlH/Gmxdb5hInXlk+LXf8/I8V04F+95e3t/d65Pikrjp5Gmkdonn6U83ck5VH865bTan3HMvC978F5svvhKttTG9yU10Lpn6pmjLzZhJ702i3lFPTloORSOLjr1JE2CzgO3RDyAVpo+EyjTIroaipWCr+QBuf8657TfGUR+jyenedttHY3M9n5ncj3OzfsVrJTk8MP9jHM0apSBOpZLqWMLsux7i9Ox0YmNMPMaKCfj9mYuQCud+b3+utD+Xgo1td209n3y7kz++MZe62Hx++6azAL5nSgKLZ9xLD09VKCa5iAklUyd0CHPNcWUlNg22snbLVeWP257tmpLrwXL2Hqyn9Ie9bNx1kE07a9m463j61b+MQ6XDEfimNoa8AbHcUXA8Q3PSOT0nncunZgIwLPfZ0MZuYVJnHhnh3O/tz6WAXu8htsxuiVw/PJsXXn0CjeLF+yr5snw/Kyv2kb53p9v305WVfFi2ncFZKeT06EpCnKkbK/xm+oQeVh7KrQ736cdbpRVsrz7M7vhbaVJZnP7QEvYd+rHipmtCLBMr/8Njb8+gd20NTf36E/Pww8SNu+ToOi21s/UxmtwZuW6vSlrGWikJ8X+tpZa73lHvcdtC+MuI4ymYPyKtz7VfzB9C0cginrjWBr91f45v7daTO+yrAIiLUeT27MrAnl3pm9aFPt2T6JPWhT5pzn8/rXyb0h8+83ge+zrPjTwXLZfQGx3N1NU7ONTQRF1DEwfrHdTVN3HgcCP76xo4UNfI/kMN7K9r5ISLJ3HjS0V0aaw/+vq6uESmnn4dC99d52weiTmJOL2Xi07qxaCsFAb36sbgrBT6/ONt1OOPHe1UkLC1Cm67FWIU2GzB1fUGSWq5RSiZ/Xjyeq556Jmc+fTjvH/hOWzcVcvGnQfZuOsgFXsPUfrDXmqPNB1d9WDsP9kX/zd0q/ee8M7NvPHlFs7tP4YN1R/yynfT2mz7loW3UHu4iRvzbLyz4XUK3zdu30VFQq89MITDB4/jppdW0OBopr6xmfqmZhqamqlvctDQ1MzhRgeH6h00+DGXYVJ8LOnJ8aw/qYCGG5sZ88aT9Kmp5lCvvmy6exoTbL/kvrQuZKYkcuE9GQA88r+/bvsm3nqI2Wxe63rDfRJEctvCevw9nkI+kJufvMY3udy5oN0EE4k2GyfjnDoPXB0A06FsSgm1RxrZUX2EbdVHuPqdW9FH6tu8d5M+wgdbnmTtxiFUJvwZR8zhNs8fbjrMXf/4DY+83ZOqxCk4Yow7F6MioTuaUmhsSGfPwQYS42JIiIuhW5c4EuJiSIyL5dPyxaiEJiYOu5quCbEkJ8aRkhhLckIcXRNjmfbpPcTE1PP29a/QPTmeLvGtOzScS0HKk0AaJTO2clr7jXu6q++jh1gk66mllluEktmPJ5/xdbA3ebcu8XTrEs/gXt3Ye2Sb23Wa2M2motHEPeS+Rt4Rs4epl5zI7SU+auhDLCoS+up7H/b6fMHcKQBMveR+t893Ka0CoHea+8kfArqy8NHdOTsug4qmvcc+HZfR8W11UHZaNhXVx8YmtdwiEGY/nsIZn7f3jolRHp/PScvmtvOP55Gvjd131ri1G0YeB8b30UOs6BNIbmj3dINzeRsB1vZ6I7XcIpTMfjyFM75g+qOEOzZ3JKEHyseogbZl+yheRNsxohc5l4c9NKuPGS4MZfbjKZzx2YbYKE4fT/9q13l8MJbi9PFt+qN427bR+y4qmlyC5m7yilDw1jaXnY2trOLYGvccY76mSi23CKWQHE9BfhP1tt2wHe92O7Z752E7el/TAcnzIPXsHy/efGzbyHNRrtAhPINQWWzyWSF8isLZlnwy2XjnvgSV0JVSo5RS3ymlNimlpoYqqI7wNRlCwIP2B8uPQXsiOZFDxPaL6LQ83o8yZOMlHv/YeD0PTTbeuS8BN7kopWKBZ4FfAFXAl0qphVrr9aEKzhdfHR4i2bnHuREbY790dXeeUd6h2MMp4vtFCJPweR76mgbTZIK5Qh8ObNJa/6C1bgBeA64ITVj+8TUov5knS4hkbGbeL0IYyee5EGVNp8HcFO0HbGn1uAo4o/1KSqlCoBAgO8R/1Xx1KDBzhwjpeCRE5PnVKQmO6WnafqRGX01JRjU1BZPQ3Y0Dq49ZoHUxUAyQn59/zPPB8NWhwJAOEQHeBIpkZw1TdBQJ8uaZVO6IUPDrXAhi3mKjBdPkUgUMaPW4P+C+n2yY+Czq9zKBhVFKVue57Ykayc4aZu8oIiInojcugxRI7FY7F4JJ6F8Cg5VSA5VSCcBYYGFowvKPz6L+Rz+geGG7zj0LXRNbRFgkO2uYvaOIEEax2rkQcJOL1rpJKXUXsBiIBeZorb8JWWR+8lq072MCi0iLZOcf6XgkhJOVzoWgeopqrT8AIn+564kZSo5M3uYmhLAOa/cUjbKSIyGESUVJL1hrj+XiKi3acfs4smo1MR5Kjqwq2r8+CiE6xhIJ3Wvi8tJbU5iX/DESouMskdB9idTUWEJYiZX/yFrl/2btNnQhhOhEOsUVuplZ5cpACBF5coUuhBAW0Tmu0KOg3EgIIYIlV+hCCGERktCFEMIiJKELIYRFdI42dOGWVNgIYS1yhS6EEBYhCV0IISxCEroQQliEJHQhhLAISehCCGERktCFEMIiJKELIYRFSEIXQgiLkIQuhBAWobTWxm1Mqd1ARYAv7wnsCWE4oSSxBUZiC4zEFphoji1Ha53p600MTejBUEqt1FrnRzoOdyS2wEhsgZHYAtMZYpMmFyGEsAhJ6EIIYRHRlNCLIx2AFxJbYCS2wEhsgbF8bFHThi6EEMK7aLpCF0II4YUkdCGEsAhTJXSl1DVKqW+UUs1Kqfx2z92nlNqklPpOKXWxh9cPVEotV0ptVEq9rpRKCFOcryulVrt+ypVSqz2sV66UKnOttzIcsbjZ5oNKqa2t4hvtYb1Rrn25SSk11aDY/qqU2qCUWquUekcp1d3DeobtN1/7QSmV6Pq8N7mOrdxwxtNquwOUUv9USn3rOifudrNOgVKqutVn/YARsbm27fUzUk5Pu/bbWqXU6QbFdUKr/bFaKVWjlJrcbh3D9ptSao5SapdSal2rZRlKqSWuPLVEKZXu4bXjXetsVEqN92uDWmvT/AA/BU4ASoD8VstPAtYAicBA4Hsg1s3rFwBjXb/PBG43IObHgQc8PFcO9DR4Hz4I3OtjnVjXPjwOSHDt25MMiO0iIM71+6PAo5Hcb/7sB+AOYKbr97HA6wZ9jn2A012/dwP+6ya2AuB9I48vfz8jYDTwIaCAEcDyCMQYC+zA2SknIvsNOA84HVjXatlfgKmu36e6Ow+ADOAH17/prt/TfW3PVFfoWutvtdbfuXnqCuA1rXW91nozsAkY3noFpZQCLgDedC2aB1wZznhd27wWmB/O7YTBcGCT1voHrXUD8BrOfRxWWuuPtdZNroelQP9wb9MHf/bDFTiPJXAeWyNdn3tYaa23a61XuX6vBb4F+oV7uyF0BfCydioFuiul+hgcw0jge611oL3Tg6a1/gzY125x62PKU566GFiitd6ntd4PLAFG+dqeqRK6F/2ALa0eV3Hswd0DONAqYbhbJ9TOBXZqrTd6eF4DHyulvlJKFYY5ltbucn3NnePh65w/+zPcJuK8gnPHqP3mz344uo7r2KrGeawZxtXMcxqw3M3TZyql1iilPlRK/czAsHx9RmY4xsbi+WIrUvsNoJfWejs4/3ADWW7WCWj/xYUkvA5QSn0C9Hbz1HSt9XueXuZmWft6S3/W8ZufcV6P96vzs7XW25RSWcASpdQG11/soHiLDXgeeAjn//0hnE1CE9u/hZvXhqR+1Z/9ppSaDjQBdg9vE5b95i5cN8vCelx1lFIqBXgLmKy1rmn39CqczQkHXfdK3gUGGxSar88o0vstAbgcuM/N05Hcb/4KaP8ZntC11hcG8LIqYECrx/2Bbe3W2YPza12c60rK3Tp+8xWnUioOuAoY6uU9trn+3aWUegfnV/ygE5O/+1ApNQt4381T/uzPgPix38YDlwIjtaux0M17hGW/ueHPfmhZp8r1madx7FfosFBKxeNM5nat9dvtn2+d4LXWHyilnlNK9dRah30AKj8+o7AdY366BFiltd7Z/olI7jeXnUqpPlrr7a5mqF1u1qnC2dbfoj/Oe4teRUuTy0JgrKviYCDOv6YrWq/gSg7/BK52LRoPeLriD4ULgQ1a6yp3TyqluiqlurX8jvOG4Dp364ZSu3bKMR62+SUwWDmrghJwfjVdaEBso4D/Ay7XWtd5WMfI/ebPfliI81gC57H1qac/RKHkaqefDXyrtX7Cwzq9W9rzlVLDcZ7Pew2IzZ/PaCFwo6vaZQRQ3dLMYBCP354jtd9aaX1MecpTi4GLlFLprmbTi1zLvDPiTm8H7giPwfmXqR7YCSxu9dx0nBUJ3wGXtFr+AdDX9ftxOBP9JuANIDGMsc4Fbmu3rC/wQatY1rh+vsHZ5GDEPnwFKAPWug6cPu1jcz0ejbNy4nsDY9uEs11wtetnZvvYjN5v7vYD8Eecf3QAuriOpU2uY+s4g/bVOTi/Yq9ttb9GA7e1HHfAXa59tAbnTeazDIrN7WfULjYFPOvar2W0qlozIL5knAk6rdWyiOw3nH9UtgONrtw2Cec9mKXARte/Ga5184EXW712ouu42wTc5M/2pOu/EEJYRLQ0uQghhPBBEroQQliEJHQhhLAISehCCGERktCFEMIiJKELIYRFSEIXQgiL+P9Bhj0WXKLMQQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# here we define the function. For use with curve_fit, we need each parameter to \n",
    "# be a separate argument (i.e., not have all arguments as a single list/array)\n",
    "def model(x,amp,cent,sigma,const,deriv=False) :\n",
    "    \"\"\" Gaussian profile function with amplitude, center, sigma, plus constant\n",
    "        with deriv=False, return array of model values at location of input x\n",
    "        with deriv=True, also return array of derivatives\n",
    "    \"\"\"\n",
    "    delta=(x-cent)**2/sigma**2\n",
    "    f= amp*np.exp(-0.5*delta)\n",
    "    if not deriv :\n",
    "        return f+const\n",
    "    else :\n",
    "        # extra code here for computing derivatives if requested\n",
    "        tmp=f*(x-cent)/sigma**2\n",
    "        d = np.array([f/amp, tmp, tmp*(x-cent)/sigma,np.ones(len(x))])\n",
    "        return f+const,d\n",
    "\n",
    "# simulate some data. Note that I've defined a par array for convenience, but pass it as\n",
    "#   individual arguments using *par\n",
    "x=np.arange(-10,10,0.5)\n",
    "par=[10.,0.,3.,1.]\n",
    "\n",
    "# here's what we would get if we assumed Gaussian uncertainties\n",
    "sig=np.sqrt(model(x,*par))\n",
    "y=model(x,*par)+np.random.normal(0.,sig,size=len(x))\n",
    "plt.errorbar(x,y,sig,fmt='ro')\n",
    "\n",
    "# here's using Poisson uncertainties\n",
    "y=np.random.poisson(model(x,*par))\n",
    "plt.errorbar(x,y,sig,fmt='go')\n",
    "plt.plot(x,model(x,*par))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that, for a Poisson distribution, the probability of getting an observed value, $y_i$, given an underlying model value, $f_i$, at each $x_i$ is: \n",
    "$$P(x_i|f_i) = {\\exp(-f_i) f_i^{x_i} \\over x_i!}$$\n",
    "\n",
    "what is the expression/function for the log(likelihood) of this data set? Remember the probability of the data set is just the product of the probabilties of the individual points.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b> ANSWER: </b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To maximize the likelihood, we will minimize the log(likelihood). For a Poisson distribution, the probability of getting an observed value, y, given an underlying model value, f, is: \n",
    "$$P(x|f) = {\\exp(-f) f^x \\over x!}$$\n",
    "Take the log to get:\n",
    "$$ ln(P(x)) = -f(x) + x ln(f) - ln(x!)$$\n",
    "Multiplying the probabilities of each individual point is taking the sum of the logs. Since the last term is independent of the model parameters, we can neglect it when looking for a maximum with respect to parameter values.\n",
    "$$ln(L) = \\sum -f_i + y ln(f_i)$$\n",
    "where $f_i$ is the model value at each $x_i$\n",
    "<br>\n",
    "Write a routine to routine the -ln(likelihood), which will use the data values, as well as the model values given the independent variable values. Include an option to return the vector of derivatives with respect to each parameters as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-132.1069679423997\n",
      "      fun: -140.25780751783714\n",
      " hess_inv: array([[ 1.25663502, -0.00775656, -0.17409553, -0.00186432],\n",
      "       [-0.00775656,  0.08412766, -0.01932584,  0.024792  ],\n",
      "       [-0.17409553, -0.01932584,  0.16306609, -0.13126844],\n",
      "       [-0.00186432,  0.024792  , -0.13126844,  0.1817842 ]])\n",
      "      jac: array([ 0.00000000e+00, -1.90734863e-06,  1.90734863e-06,  1.90734863e-06])\n",
      "  message: 'Optimization terminated successfully.'\n",
      "     nfev: 70\n",
      "      nit: 11\n",
      "     njev: 14\n",
      "   status: 0\n",
      "  success: True\n",
      "        x: array([9.51996749, 0.312664  , 2.77652864, 1.18851471])\n",
      "-140.25780751783714\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12fc4e850>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD6CAYAAAC4RRw1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU9bnH8c9D2BdZLiA7KKVcbXHBXKzikhYX0CIuuMYqym1ciutt1cp1N23VqnUrNhZFba5o3RCXKlKxLhUNKgYEyyJBBAFFQAhCIL/7x29SQ5hJJrOdOZPv+/WaVyZnfXLOmSdnfuf3nGPOOUREJHyaBR2AiIgkRglcRCSklMBFREJKCVxEJKSUwEVEQkoJXEQkpBpM4GbW18xeM7MFZjbfzC6JDO9iZjPMbFHkZ+f0hysiIjWsoX7gZtYT6Omce9/MOgBzgOOBccA659zvzOwqoLNz7sr6ltW1a1c3YMCAlAQuItJUzJkz50vnXLe6w5s3NKNzbhWwKvL+GzNbAPQGxgAFkckeBmYB9SbwAQMGUFZW1qjARUSaOjOriDa8UW3gZjYA2B+YDeweSe41Sb57ciGKiEhjxJ3Azaw98BRwqXNuYyPmKzKzMjMrW7t2bSIxiohIFHElcDNrgU/epc65pyODV0fax2vayddEm9c5V+Kcy3fO5XfrtksTjoiIJCieXigGTAYWOOfuqDXqOeDsyPuzgWmpD09ERGJp8CImMBz4GVBuZh9Ghl0N/A54wszGA8uBk9MTooiIRBNPL5Q3AYsxekRqwxERkXipElNEJKSUwEVEQkoJXHJCwZQCCqYUBB2GSEYpgYuIhJQSuIhISCmBi4iElBK4iEhIKYGLiISUEriISEgpgYuIhJQSuIhISCmBi4iElBK4iEhIKYGLiISUEriISEgpgYuIhJQSuIhISCmBi4iElBK4iEhIKYGLiISUEriISEgpgYuIhJQSuIhISCmBi4iElBK4iEhIKYGLiISUEriISEgpgYuIhJQSuIhISCmBS+iVlpfyzop3eL3idQb8YQCl5aVBhySSEUrgEmql5aUUTS9i646tAFRsqKBoelFKk3jBlAIKphSkbHkiqaIELqE2ceZEKqsqdxpWWVXJxJkTA4pIJHOUwCXUlm9Y3qjhIrlECVxCrV/Hfo0aLpJLlMAl1IpHFNO2RdudhrVt0ZbiEcUBRSSSOUrgEmqFQwopGV1Cq7xWAPTv2J+S0SUUDikMODKR9GsedAAiySocUsgDcx4AYNa4WcEGI5JBOgMXEQmpBhO4mT1oZmvMbF6tYdeb2edm9mHkdUx6wxQRkbriOQOfAoyMMvxO59x+kdeLqQ1LJLVUjCO5qMEE7pz7B7AuA7GIiEgjJNMGPsHMPoo0sXROWUQiIhKXRBP4JGAgsB+wCrg91oRmVmRmZWZWtnbt2gRXJyIidSWUwJ1zq51zO5xz1cADwLB6pi1xzuU75/K7deuWaJwiTY7a7aUhCSVwM+tZ69cTgHmxphURkfRosJDHzB4DCoCuZrYCuA4oMLP9AAcsA85LY4wiIhJFgwncOXd6lMGT0xCLiIg0gioxRURCSglcRCSklMBFREJKCVwkjdQVUNJJCVxEJKSUwEVEQkoJXEQkpJTARURCSglcRCSklMBFcpR6wOQ+JXCRepSWl/LOind4veJ1BvxhAKXlpUGHJPJvSuAiMZSWl1I0vYitO7YCULGhgqLpRUrikjWUwEVimDhzIpVVlTsNq6yqZOLMiQFFJLIzJXCRGJZvWN6o4YlQO7UkQwlcJIZ+Hfs1arhIpimBS8ak82xz1rhZzBo3K6XLLB5RTNsWbXca1rZFW4pHFKd0PSKJUgIXiaFwSCElo0toldcKgP4d+1MyuoTCIYUBRybiNfhEHpGmrHBIIQ/MeQAg5Wf4IsnSGbiISEgpgUvOC2MxThhjlsxTApecFsZinDDGLMFQApecFsZinDDGLMFQApeclolinFQLY8wSDCVwyWlhLMYJY8wSDCVwyWlhLMYJY8wSDCVwyWlhLMYJY8wSDBXySM4LYzFOGGOWzNMZuIhISCmBi4iElBK4SJoEWU2pSs6mQQlcJA2CrKZUJWfToQQukgZBVlOqkrPpUAIXSYMgqylVydl0qBuhSMSmTbBsGXz2GXzzDVRXQ5s2sHHJXrTuurpRy+rXsR8VGyqiDk+3INctmaUELk3Wtm3w2mswbRq8+SbMmwfORZtyEgC9boFhw+DYY+GnP4WePWMvu3hEMUXTi3ZqyshUNWWQ65bMUgKXJmfJErj3XpgyBdavh/btYfhwOOEE2Htv6NsXOnWCZs1g82YY95er2LK6Nwc3v4g33vAJ3wxGjYILLvA/8/J2XkdN1eT4aePZumMr/Tv2p3hEcUaqKYNct2SWErg0GZVf9OaUU+DJJ33CHTsWzjgDjjwSWreOPd9/lL8DwCPjLsI5mD8fHn8cJk+G0aN90r/+ejjpJJ/0awRZTalKzqahwYuYZvagma0xs3m1hnUxsxlmtijys3N6wxRJ3Ndfw6JHL+G9iY/w4ovw619DRQU89phPwPUl77rM4Ic/hJtu8suYOtU3u5xyim9eKStL398hUlc8vVCmACPrDLsKmOmcGwTMjPwuElNQhSVPP+3PkFfOOo5eBc+xZAkUF0OvXt9NUzClgIIpBY1edosWcOqpUF4OjzwCK1f6JH7RRf6CaENUbCPJajCBO+f+AayrM3gM8HDk/cPA8SmOS3JIEIUlGzf65pGTToIePWDotecx6Gd3sfvuqV/XiEcLmLyjgAULfPK+7z4YOhTmzIk9j4ptJBUS7Qe+u3NuFUDkZ/fUhSS5JtOFJXPnQn4+PPGEb+p4913o0H9xWtZVW8eOcNddvmfLli1w0EHw2csnR+3ZomIbSYW0X8Q0syKgCKBfP/VDbYoyWVhSWgr//d/QuTP8/e9w2GEpX0WDDj/c/xP5+c/h6am/YFPFIL49fee2dhXbSCokega+2sx6AkR+rok1oXOuxDmX75zL79atW4KrkzDLxCPCnIMbboAzz4QDD4QPPwwmedfo0sX3dhlw4p9Z88+jOPxw30ZeQ49Nk1RINIE/B5wdeX82MC014UguSvcjwrZtg7PP9l35zj4bXnkFumdBo54Z9B/9F35w0f8yf77va75okR+nx6ZJKsTTjfAx4J/AYDNbYWbjgd8BR5rZIuDIyO8iUaXzEWFbtvgCnEcf9e3dDz0ELVsmvdidzBo3K6m+1F2HvsmsWb5nyiGHwAcfpGabJNp7RnJHg23gzrnTY4wakeJYJIelo7Bk0yYYM8ZfNPzTn6CoKCWLTYv8fF+uf+SRUFAAzz8PhYeq2EaSo7sRSiht2AAjR8KsWb4PdjYn7xqDB8Pbb/s+6KNG+YQukgwlcAmdzZvhmGNg9mxf0n7mmUFHFL8+ffw3hj59fBLfuGTvoEOSEFMCl1D59ls4/nh45x1fxj52bNARNV6PHr6LY48e8NHtt7Jx6eCgQ5KQ0s2sJDSqquC00+DVV/2dBE86KeiIEterl0/i3z9gI+V33Mb8U+EHP9h5GrWLS0N0Bi6hUF0N55zjb+V6zz2+u2DY9e0L+15xOc1abGPkSP8gCZHGUAKXULjiCl9lWVwMEyYEHU3qtO76BUMuv4KNG+Hoo2Fd3bsOidRDCVyy3n33we23wy9+4W8Fm2va913KtGn+QROjR0NlZcPziIASuKRYqotLpk+Hiy/2ie2uu3x1YyKSLcZJt4IC/w3jn/+EwkLfZCTSECVwyVplZf6i5dCh/uELdR9blmvGjoU774Rnn4Wrrw46GgkD9UKRRqk5u0732eyyZf7Bwd26+bPwdu3SurqscfHFsHAh3HKLL/w555zEl5XN3zgkNZTAJets2gTHHQdbt/qilx49go4oc8zg7rt9e/h558Gee/rb04pEoyYUySrV1XDWWf7BwU88AXvtFXREmdeihf/bBw6EE0/87g6GInUpgUtWuekmeOYZ+P3v/Y2fmqpOnfwNr8x8U9KGDUFHJNlICVyyxtNPf3dP70svDTqa4A0c6LfJ0qXqmSLRKYFLVigv900nBx4I99+feHfBXHPYYb775AsvwHXXBR2NZBtdxJTAffmlv2i5227+jLP2syPDLhU9QS64wD/h/uabYf/9fbu4COgMPLRy5WksVVVwyimwapXv/9yrV9ARZR8zX406bJhvXpo/P+iIJFsogUug/ud/fFfBkhKfoCS61q39t5N27fztdL/+OuiIJBsogUtgHnnE31nwsst8+7fUr3dveOopqKjwFzVdtT6+TZ2OAAnE++/7QpWCArj11qCjCY/hw/0/vZdegmXPJFGmKTlBCVx2ke729S+/9E+S79rVPxKteYCX0kvLS3lnxTu8XvE6A/4wgNLy0ozOn4jzzoMfX1zK8u9dw+vLMrdeyT7qhSIZM2vcLLZv9w8jXr0a3ngDuncPLp7S8lKKphexdcdWACo2VFA03T8duXBIYdrnT1RpeSmzuxfB9sqMrleyj87AJaOuvhpmzoQ//hH+67+CjWXizIlUVu188+3KqkomzpyYkfkTNXHmRCq3Z369kn10Bi4Z88QTcNttcP75cO65QUcDyzcsb9TwVM8PifUTT8V6JTfoDFwyYt48n7QPOshXFmaDfh37NWp4qudPVKzld3DpXa9kHyVwSbv16/1Fyw4d4MknoWXLoCPyikcU07ZF252GtW3RluIRxRmZP1HR1tu8ui0bny7mySfTumrJMkrgklbV1fCzn/kHNPz1r9lVaVk4pJCS0SW0ymsFQP+O/SkZXRL3hcBk509UtPX+eUwJB7UvZNw4/21Hmga1gWepTD35Jt1uvNHfFvXee+GQQ4KOZleFQwp5YM4DQGLbOtn5ExVtvUc+CQcc4Cs133sPOnfOWDgSEJ2BS9pMnw433ODv33HhhUFHk/t69fKVmsuXwxlnwI4dQUck6aYEHkJBFI/Es97a43vfNoBTbipl6FCYNEm3h82Ugw/233b+9je45pqgo5F0UxNKyARZPFLfeuuOX1lZAUcVcdZPoE0bFZdkUlERlJXBb38LQ4f6p91LbtIZeMgEWjxSz3qjjadFJXd+lLq4cuUWuplwzz2+y6YuauY2JfCQCaqIo6H1qrgku7Rq5btsduig28/mMiXwkMm24pGa4UHFJbHpombuUwIPmWwqHqm93nP7F8O2zMcl9dNFzdymBB4y2VQ8UrPelSvh/l8U0vWfJbRq1iajcUnDior867e/hb/8JehoJJXUCyWEsql4ZMsW38b6zTfw9l2FXBRAXNKwe+6BRYtg/Hjo3x8OPTToiCQVkjoDN7NlZlZuZh+aWVmqgpJwqK72RTrvvefP7IYMCToiiaVlS98evsce/r40ixcHHZGkQiqaUH7snNvPOZefgmUJwRXqNHbdV13l729y220wZkzGQpQEde4ML7zg3x97LKxbF2w8kjy1gWeZWAUzmUjijVn3pEk+cV94oX+yvITDwIHw7LP+5mInnQTbtgUdkSQj2QTugFfMbI6ZFaUioKYuqEKdxqz7q7k/YsIEfxZ3110qkw+bQw6BBx+EWbP88zWdCzqi1GlqxV7JXsQc7pxbaWbdgRlmttA594/aE0QSexFAv37qE9yQIAti4ln3N8u+z8eTrmO//WDq1NQ+kDhX7sAYBoWFvh38+ut9f/Fi9fYMpaTOwJ1zKyM/1wDPAMOiTFPinMt3zuV369YtmdU1CUEWxDS07iVLYN4ffkuL9ht4/nlo3z7tIUkaXXut7174m99kz1OSpHESTuBm1s7MOtS8B44CdNeFJAVVqNPQuleuhCOPhOodzRly2ZX07Jn2cCTNzPzDpU88ES69FEobcZmlqTVVZKtkzsB3B940s7nAu8ALzrm/pSaspiuoQp361j2qdyFHHQVr18KQy66gXe+KtMcimZGX5xN3QYG/8dXf9AkOlYRbMJ1zS4F9UxiLRARVqBNt3Zs2wRFH+PbSl16CGyo+yWg8kn6tW/ueKYcf7numzJjhS/Al+6kbocT07be+6KOsDB5/HH7846AjknTp2NGffffqBSNHwuzZ6VuXml9SR6X0ElV1VUuOPx5mzoSHHlKhTjZK9bezHj3gtdf8mfjRR8Orr0K+yvOyms7AQ2rWuFlpa17Zsa0l8+66mVdegT//2ZfLS26JdRbcp49P4p07+4vW77+f+dgSFWQFc1CUwGUnW7bA/Ltv5uuP85k8Gc49N+iIJNP69fNJfLfdfBL/8MOgI2pYkBXMQVICl3/btAmOOw6+/jifwefcyjnnBB2RBGXAAJ/E27XzPVTeeivoiOoXZAVzkJTABYCvvoIRI/yHdvC5t9DjUPUna+r23BPeeAO6d/dn4i+/HHREsTXVR/opgQsrVvj7Q8+d62852uOQLP6kSkb17++T+ODBMHq0v/tkNmqqj/RTAm/iFi6E4cN9En/5ZfU2SUQ6Lyhng91399/Mhg2D006D++4LOqJdBVnBHCQl8Cbs1VfhRz/y/b1ff913H2uKcj0Bp0KnTvDKK/DTn8KECbC4dAKuOnr6CKKfd5AVzEFSP/AmatIkuOgi2GsveP55/1U5Hkp0TVfbtvD003DFFXDHHWPZsrY335wMHToEHZkXZAVzUHQGnqXSdVZYVQUXX+wfxDByJLz9dvzJWyQvD26/HQaddQfryocxfHhqH8+mKs3GUQJvQj7/HH7yE/+A28svh2nTUn/2lK5/PE2xSCNblZaXsnzQ1XBNc+YfOYAhZ5Ty7LPxz6v9mDpK4EkI09nCjBmw337wwQfwf//nz6Ly8oKOKj5NtUgjG+20Lwyqd6tg29FFnHBtKVdeCdu3xzkv2o+p0OQTeENJOExJOppt2+Caa/y9Lbp390+QP/30oKNqnKZapBGkWMd9tH1R3byS9mMmcuut/kL4ljXRbxav/Zh6oUjgYU+iQZk/3/cyuflmfz+Td9/1Fy3DpqkWaWSjWNt8c/PllJb6Y67s2gdZ9fqxuzxrU/sx9UKRwKVxDp88goGnTuKAA3z/7mee8XcUbNcu6MgS01SLNLJRffvijDOgvBx22/Nj/jXlV4wZ44+/eOaVxCiB55jZs+H9Gyex9IkLOOYYmDcPjj++ccvItn7RTbVIIxs1tC/69oV9fvlLBp5+L6++6r/x3XWXbxvPxH7MtmM33ZTAc8S6dXD++XDQQbDtm87sfeF1PPWUb/cOu6ZapJGN4tkX1szR56gnmTfP36Lh0kt9FefAzdqPqaZCnpDbsgXuvdc/Wfybb+Cyy2D2986ieZstmAUdXeo0xSKNbBXvvthzT3jhBX9/nUsu8ScXY8cWss8BL9K2z+fajymgM/CQqqqCBx+E73/fV8YdfLC/b/Ptt0PzNluCDk8EADMYOxY++QSuv94/V7Xsfx9m0V8uZuXKoKMLv6xP4Ons+N/QsusbH1RBwpYt/mZCgwbB+PH+GYavvebPdH74w4yEkBYq8MicID5T7dvDddf5qs0eh73AytfGsMcevtlv6dL0x5WrsroJJVbHfyDpdrOGll3feCCuuGq6Pqbiq+IXX/jHm91zD6xZ47+O3nOPv7lQ2JtK0rmfZWdBfqbAP3fz+2fdSd+RU9l/6WM89BA88AAMG1/KB32L2Fqduc9ULsjqM/B0dvxvaNn1jc9UQUJ1tX+o8Mkn+6v711wDQ4f6Owe+9Za/P3PYkzeowCOTgvxM1dam+yruvx8+/dRf5JzdYSJbq4M9BsJYb5LVZ+Dp7Pjf0LITWXcq4nLOP1hh6lT/qqiALl38RaDzzvNNJ8nKtrMYFXhkTpCfqWh69fLXbe68Ifa8334LrVsnHV5SGvrMBPWZyuoz8HR2/G9o2fWNT3Vc27f7p55cfbXvN7v//v6g3ntvePRRfxOq3/8+Nck7G6nAI3OC/EwlMq9b34/dd4czz4THHoOqTVly79oskdUJPJ0d/xtadn3jk42ruho+/ti3/Z16KnTrBocdBrfd5s9I/vQnWLUKXnzRH7hBn32kmwp1MifIz1Sj523elivziznxRP+0qDPOgLcvfpYPfnMPN97oL95XVsZYYBOR1U0oNRcvxk8bz9YdW+nfsT/FI4pTcmGroWXHs+5441q92jeLfPCBb7t+6y1feAPQsyecdBKMGgVHHAEdOyb9p4VOOvez7CzIz1TC8xbCjh1QVgan/uZRvvroIK6/3jc3tmgBBxwAhxzif+63n/+mGpY7bSYrqxM4pLeAo6Fl1ze+7rgdO2DZMliyxHeVWrQIPvrbrWz67Hv0OOe7+QYPhhNO8Afc8OHwve/lxoXIZKlQJ3OC/EwlOm9eHhx4IAw4YQoDTpjCs8fP4u23fdPjm2/C3Xf7O28CtGkD++wDQ4b4ZF7zGjjQj8slWZ/Ag1a9PY/tlR341798972VK33zxqpVsPCNq9i6viv/+Tvfl7Wq6rv5WraElj060WWf2Vx10ij23Rf23ddfkBRpylLxT6NTJzjmGP8Cn7wXLvTFbDWvadNg7drv5jGD3r2hTx/fVNm793evXr1g0/KBNG+3kc2b/ePjwnBiFYoEvn1LG6q3tWHFCv9Vavv2XV/RhtcMq6ryBTCVlTv/3LIFFs25hOqqlpz+MmzeDF9/DevXf/dz8+aZAAyuE1OLFtBst/1o2ekrhuzvbxg1cKA/ox440B8UIx71/VgvHTdql78p23qCiIRZy5b+rHuffeCss74bvn79d9+IFy3yJ1qffw4LFviHem/cWHspkwFo/0u/vC5d/KtTJ1iw4VbyWn3Lz2b6u3q2b+9/1rxWlR2LNa9iauvIyVutV6tW/uegQbDbbqn9u0ORwD/963msfO14+l6a2uW2bg3b835Cs5bbKFvld0Tnzj4Jd+7sX09/OpnmbTdx3chL6N7dt1n36uV37I8fPg2Av0ZJwjVVZVt3bGXAHwZkrE03qPVK7sjWYyiRuDp1gvx8/4pm0yaf0B+dW8qt8y6kqtlGOjXrzyFbi+mxppB162Bhi1LWH/wrXLsveGxzP9q+XUz13EI2b44sZEgpjLgPOi7n9NklMLMYyneN66WX/HNoUykUCbzbgX+nXZ+l/OrQy2nenJ1eeXnENaxNG/+1qOZnq1bQrBkUTBkDxD4TnjPlUQDOPPOSuOMNqrJQFY2SrHiPoUx/c0zFsR3tW2/79lC2rZQ7FxdRlee7tKx3Ffy9XRElF/tpiqYX4SIFSjvaV7Dj2CJK7ofTf1DIlDmlTHiliC3bI91hOlXQ+pQirr4JRnQrZOtW37yzbZsvwku1UCTwToM/otPgj/j5uMujjq+vOSKIpor6KtLSmUiDWq/kjmw9htIZV0MVpPWt98a3Jn6XvCO+3VHJ5E8ncs2Y9G+vUCTwdCbfdCw7qMpCVTRKsrL1GMq2CtJkKrZTKasLebJBfU/4iDUuqMpCVTRKsrL1GAqqgjSZiu1MUAJPg6AqC1XRKMnK1mMoqArSZCq2MyEUTShhE1RloSoaJVlBH0Oxvu0GXUGaTMV2OimBp0lQlYWqaJRkZesxFFQFaTIV2+mWVBOKmY00s0/MbLGZXZWqoEREpGEJJ3AzywPuA0YBewOnm9neqQosXtn42LOGBBlXtm4TyR31XfgPSq4e98k0oQwDFjvnlgKY2VRgDPBxKgKLRyoee5ZpQRbbqNBHmqJcPu6TaULpDXxW6/cVkWEZkw2PPWusIOPK1m0ikk65fNwncwYe7V5dbpeJzIqAIoB+/VLbNzKox54lI8iO/0EXHYgEIRXHfUNNQkE1GSWTwFcAfWv93gdYWXci51wJUAKQn5+/S4JPRr+O/ajYUBF1OFDvuKA0FHOurjsVsq1dVcIh3uM+jMdXMk0o7wGDzGwPM2sJnAY8l5qw4pPOx56lQrSLOUHGlQ3bRLJDNl5ojEcicefycZ/wGbhzbruZTQBeBvKAB51z81MWWRxS+dizTAmy43/QRQciQcjl4z6pQh7n3IvAiymKJSHJdMAPSpBxZes2EUmnXD3udS8UEZGQUgIXEQkpJXARkZDSzaxyUC618YlIbDmRwJWwcoP2o0jj5EQCr4+Sgkjj5eLnJhf/JrWBi4iElBK4iEhI5XwTSrbKxa9zIpJZOgMXEQkpJXARkZBSAhcRCSklcBGRkFICFxEJKfVCaYLUA0YkN+gMXEQkpJTARURCSglcRCSklMBFREJKCVxEJKSUwEVEQkoJXEQkpJTARURCSglcRCSkzDmXuZWZrQUqEpy9K/BlCsNJFcXVOIqrcbI1Lsje2HIxrv7OuW51B2Y0gSfDzMqcc/lBx1GX4mocxdU42RoXZG9sTSkuNaGIiISUEriISEiFKYGXBB1ADIqrcRRX42RrXJC9sTWZuELTBi4iIjsL0xm4iIjUklUJ3MxONrP5ZlZtZvl1xv3azBab2SdmdnSM+fcws9lmtsjMHjezlmmI8XEz+zDyWmZmH8aYbpmZlUemK0t1HFHWd72ZfV4rtmNiTDcysg0Xm9lVGYjrNjNbaGYfmdkzZtYpxnQZ2V4N/f1m1iqyjxdHjqUB6Yql1jr7mtlrZrYgcvxfEmWaAjPbUGv/XpvuuCLrrXe/mHd3ZHt9ZGZDMxDT4Frb4UMz22hml9aZJmPby8weNLM1Zjav1rAuZjYjkotmmFnnGPOeHZlmkZmd3eiVO+ey5gXsBQwGZgH5tYbvDcwFWgF7AEuAvCjzPwGcFnl/P3BBmuO9Hbg2xrhlQNcMbrvrgV82ME1eZNvtCbSMbNO90xzXUUDzyPtbgFuC2l7x/P3AhcD9kfenAY9nYN/1BIZG3ncA/hUlrgLg+UwdT/HuF+AY4CXAgB8BszMcXx7wBb6fdCDbCzgMGArMqzXsVuCqyPuroh33QBdgaeRn58j7zo1Zd1adgTvnFjjnPokyagww1Tm31Tn3KbAYGFZ7AjMz4CfAk5FBDwPHpyvWyPpOAR5L1zrSYBiw2Dm31Dm3DZiK37Zp45x7xTm3PfLrO0CfdK6vAfH8/WPwxw74Y2lEZF+njXNulXPu/cj7b4AFQO90rjOFxgCPOO8doJOZ9czg+kcAS5xziRYIJs059w9gXZ3BtY+jWLnoaGCGc26dc+5rYAYwsjHrzqoEXo/ewGe1fl/Brgf4f21OpVQAAANnSURBVADrayWLaNOk0qHAaufcohjjHfCKmc0xs6I0xlHbhMjX2AdjfGWLZzum07n4s7VoMrG94vn7/z1N5FjagD+2MiLSZLM/MDvK6IPMbK6ZvWRmP8hQSA3tl6CPqdOIfRIVxPaqsbtzbhX4f9BA9yjTJL3tMv5QYzN7FegRZdRE59y0WLNFGVa3+0w808QlzhhPp/6z7+HOuZVm1h2YYWYLI/+pE1ZfXMAk4Cb833wTvnnn3LqLiDJv0t2Q4tleZjYR2A6UxlhMyrdXtFCjDEvbcdRYZtYeeAq41Dm3sc7o9/HNBJsi1zeeBQZlIKyG9kuQ26slcBzw6yijg9pejZH0tst4AnfOHZHAbCuAvrV+7wOsrDPNl/ivb80jZ07RpklJjGbWHDgROKCeZayM/FxjZs/gv74nlZDi3XZm9gDwfJRR8WzHlMcVuTjzU2CEizT+RVlGyrdXFPH8/TXTrIjs547s+vU45cysBT55lzrnnq47vnZCd869aGZ/NLOuzrm03vMjjv2SlmMqTqOA951zq+uOCGp71bLazHo651ZFmpTWRJlmBb6tvkYf/PW/uIWlCeU54LRID4E98P9J3609QSQxvAaMjQw6G4h1Rp+sI4CFzrkV0UaaWTsz61DzHn8hb160aVOlTrvjCTHW9x4wyHxvnZb4r5/PpTmukcCVwHHOucoY02Rqe8Xz9z+HP3bAH0t/j/VPJ1UibeyTgQXOuTtiTNOjpi3ezIbhP7tfpTmuePbLc8BZkd4oPwI21DQdZEDMb8FBbK86ah9HsXLRy8BRZtY50uR5VGRY/DJxlbYRV3NPwP9X2gqsBl6uNW4ivgfBJ8CoWsNfBHpF3u+JT+yLgb8CrdIU5xTg/DrDegEv1opjbuQ1H9+UkO5t9yhQDnwUOXh61o0r8vsx+F4OSzIU12J8O9+Hkdf9dePK5PaK9vcDN+L/wQC0jhw7iyPH0p4Z2EaH4L86f1RrOx0DnF9znAETIttmLv5i8MEZiCvqfqkTlwH3RbZnObV6j6U5trb4hNyx1rBAthf+n8gqoCqSv8bjr5vMBBZFfnaJTJsP/LnWvOdGjrXFwDmNXbcqMUVEQiosTSgiIlKHEriISEgpgYuIhJQSuIhISCmBi4iElBK4iEhIKYGLiISUEriISEj9P8rMTcLE5sK6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def logl(pars,x,y,deriv=False) :\n",
    "    \"\"\" log(likelihood) function\n",
    "        Poisson p(x|f) = exp(-f)*f(y) / y!\n",
    "        ln(p) = -f + y*f - ln(y!)\n",
    "        but last term is independent of pars\n",
    "    \"\"\"\n",
    "    mod=model(x,*pars,deriv=deriv)\n",
    "    if not deriv :\n",
    "        tmp = np.sum(-mod+y*np.log(mod))\n",
    "        bd=np.where(mod <=0)[0]\n",
    "        if len(bd) > 0 : print('bad: ',pars)\n",
    "        return -tmp\n",
    "    else :\n",
    "        tmp = np.sum(-mod[0]+y*np.log(mod[0]))\n",
    "        return -tmp,-np.sum(-mod[1] + y*mod[1]/mod[0],axis=1)\n",
    "\n",
    "\n",
    "from scipy.optimize import minimize\n",
    "p0=[8,1,3,0]\n",
    "p0=[10,0,3,2]\n",
    "print(logl(p0,x,y))\n",
    "# default BFGS algorithm, numerical derivatives\n",
    "out=minimize(logl,p0,args=(x,y))\n",
    "print(out)\n",
    "print(logl(out['x'],x,y))\n",
    "\n",
    "xfit=np.arange(-10,10,0.01)\n",
    "plt.errorbar(x,y,sig,fmt='go')\n",
    "plt.plot(xfit,model(xfit,*out.x),color='b')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will use scipy.optimize.minimize() to find the minimum of our -ln(likelihood) function. You'll need to supply a starting guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      fun: -135.02762072896473\n",
      " hess_inv: array([[ 1.31411877, -0.01418979, -0.14358744, -0.03101242],\n",
      "       [-0.01418979,  0.09643986,  0.0179437 , -0.01432622],\n",
      "       [-0.14358744,  0.0179437 ,  0.08565652, -0.05610299],\n",
      "       [-0.03101242, -0.01432622, -0.05610299,  0.10216768]])\n",
      "      jac: array([0.00000000e+00, 0.00000000e+00, 3.81469727e-06, 1.90734863e-06])\n",
      "  message: 'Optimization terminated successfully.'\n",
      "     nfev: 100\n",
      "      nit: 17\n",
      "     njev: 20\n",
      "   status: 0\n",
      "  success: True\n",
      "        x: array([ 9.82479012, -0.22915828,  2.76639097,  0.91959623])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12352f850>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD5CAYAAAA+0W6bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXwU9f3H8deHcIkoIAJqFaIt3lSl0VbrEcUDUUSpRTT1RGOrVG219chD5YdNPVr9eR94gRoRC943oqD2gQfwQ0HBGpFLELAKiGgg5Pv747uRGHaTTXZmZ2fzfj4e+8hmdnbmM5PvfjL7nfnM15xziIhI/LSKOgAREWkeJXARkZhSAhcRiSklcBGRmFICFxGJKSVwEZGYat3YDGa2A/AQsA1QA4xyzt1iZlsB44BCYD4wxDn3dUPL2nrrrV1hYWGGIYuItCzTp0//0jnXrf50a+w6cDPbFtjWOTfDzLYApgPHA2cAXznnrjOzy4AuzrlLG1pWUVGRmzZtWnO3QUSkRTKz6c65ovrTG+1Ccc4tdc7NSDz/BpgD/AQYBIxJzDYGn9RFRCRLmtQHbmaFwD7AO0AP59xS8Eke6B50cCIiklraCdzMOgITgIucc6ub8L5SM5tmZtNWrFjRnBhFRCSJtBK4mbXBJ+8K59wTicnLEv3jtf3ky5O91zk3yjlX5Jwr6tZtkz54ERFppkYTuJkZcD8wxzl3U52XngFOTzw/HXg6+PBERCSVRi8jBH4NnArMMrOZiWlXANcBj5vZMGAh8NtwQhQRkWQaTeDOubcAS/Fyv2DDERGRdKkSU0QkppTApUUoHl1M8ejiqMMQCZQSuIhITCmBi4jElBK4iEhMKYGLiMSUEriISEwpgYuIxJQSuIhITCmBi4jElBK4iEhMKYGLZEhVnhIVJXARkZhSAhcRiSklcBGRmFICFxGJKSVwEZGYUgIXEYkpJXARkZhSAhcRiSklcBGRmFICFxGJKSVwEZGYUgIXEYkpJXARkZhSAhcRiSklcBGRmFICFxGJKSVwEZGYUgKXvFcxq4K3F7/NlAVTKLy5kIpZFVGHJBIIJXDJaxWzKih9tpSqDVUALFi1gNJnS5XEJS8ogUteK5tUxtr1a380be36tZRNKosoIpHgKIFLXlu4amGTpovEiRK45LWenXo2abpInCiBS14r71dOhzYdfjStQ5sOlPcrjygikeAogUteK+lTwqiBo2hX0A6AXp16MWrgKEr6lEQcmUjmWkcdgEjYSvqUcO/0ewGYfMbkaIMRCZCOwEVEYqrRBG5mD5jZcjObXWfaCDP73MxmJh4Dwg1TRETqS+cIfDTQP8n0/3XO7Z14vBBsWCLxoCpPiVKjCdw59wbwVRZiEYkVVXlK1DLpAx9uZh8kuli6BBaRSEyoylOi1twEfhfwU2BvYClwY6oZzazUzKaZ2bQVK1Y0c3UiuUdVnhK1ZiVw59wy59wG51wNcC+wXwPzjnLOFTnnirp169bcOEVyjqo8JWrNSuBmtm2dX08AZqeaVyRfqcpTotZoIY+ZjQWKga3NbDFwNVBsZnsDDpgPnBtijCI5qbaac9jTw6jaUEWvTr0o71euKk/JmkYTuHPu5CST7w8hFpHYUZWnREmVmCIiMaUELtKI4tHFFI8ujjoMkU0ogYuIxJQSuIhITCmBi4jElBK4iEhMKYGLiMSUEriISEwpgYuIxJQSuIhITCmBi4jElBK45AVVS0pLpAQuIhJTSuAiIjGlBC4iElNK4CIiMdXogA4i+UCDLUg+0hG4BEpXg4hkjxK4iEhMKYGLiMSUEriISEwpgYuIxJQSuIhITCmBi4jElBK4iEhMqZBHJEMqEpKo6AhcRCSmlMBFRGJKCVxEJKaUwKXFa+j+LRWzKnh78dtMWTCFwpsLqZhVkd3gRBqgBC6SQsWsCkqfLaVqQxUAC1YtoPTZUiVxyRlK4CIplE0qY+36tT+atnb9WsomlUUUkciPKYFLzsi1W9EuXLWwSdNFsk0JXCSFnp16Nmm6SLapkEeypvboOpuFLzU1MG0avPkmvPceVFbCF1/A6tVQUAAdOsCazW6nffcl3PEtHHQQ9OkDZlDer5zSZ0t/1I3SoU0HyvuVZy1+kYYogUte+vBDuPtumDABli7103r2hN13h5//HDp1gg0bYM0aePLtKlbO2Yfhw/18O+4IQ4fC2WeXMGogDHt6GFUbqujVqRfl/cop6VMS3YaJ1KEELnnljTfgf/4HXnsN2rWDY4+FE06Aww+HHj2Sv2fe6ItxDsYcOplJk+Dxx+GGG+D662Hw4BL22GMyWxR+opJ5yTlK4JIXvl3Sk3nj/sAhH8B228G118LZZ8PWW6f3fjMoLIRhw/xjyRK47Ta46y5YNf5eehzwEosPh+23D3UzRJqk0ZOYZvaAmS03s9l1pm1lZhPN7JPEzy7hhimS3Lp1MHIkTL/6PlZX7sn11/t+7ssuSz95J1P7T2DBAthhwKMsf/cwdt3VJ/SamuDiF8lEOlehjAb615t2GTDJOdcbmJT4XVq4bFctzpkDRUVw9fgK7OIdqP7rltzZtpAnKoNbb6dOsNNvR7Hv30/jgAPgvPPgiCNg0aLAViHSbI0mcOfcG8BX9SYPAsYkno8Bjg84LomZbFctPvII7LsvzN+ignYnllKz2fJQ17tZty94+WW45x54913o2xdefTXQVYg0WXOvA+/hnFsKkPjZPbiQJI6yVbVYXQ3Dh8Opp/okuuXgMqpcdqolzaC01F+W2KMHHHWU72ZxLvBViaQl9EIeMys1s2lmNm3FihVhr04iko2qxW++gUGD4I474OKL/ZUmS9Zkv1pyl13g7bfhpJPgiiv8Sc/165u3rFyrPpV4aW4CX2Zm2wIkfi5PNaNzbpRzrsg5V9StW7dmrk5yXdhVi0uXwsEHw8sv++u7//lPaN06umrJjh2hogKuvhoefBAGDvT/YESyqbkJ/Bng9MTz04GngwlH4qq8Xzkd2nT40bSgqhYXLfLJu7ISnn8ezj03O+ttjBmMGAH33ef7ww87DL7+OvTVivwgncsIxwJTgV3MbLGZDQOuA44ws0+AIxK/SwtW0qeEUQNH0a6gHQC9OvVi1MBRGVctfvaZT94rVsDEib7fORvrbYphw+Cpp+CDD3zB0Ff1T/mLhKTRQh7n3MkpXuoXcCwScyV9Srh3+r1AMPc7qayEQw+Fb7+FSZPgF7/Iznqb49hjfRI/4QTo188fkXftGkko0oLoboTSJNk66bZ4sT+a/e47eP311Mk7lxx9NDz9NMyd62NftSraeHSCNP8pgUvOWbHCF8t89ZU/abnXXlFHlL6jjoInn/Q30zruOP8PSCQsSuCSU6q/68DRR8P8+fDss/E48q6vf3946CF/C9uTTvLXrouEQTezkpxRU92aD2/9G2sqfX/yIYdEHVHzDR3qv0Gcfz5sX/wSS9/sj1nUUUm+0RG45ATn4D9jLmbl3L7cfz8cc0zUEWXuvPOg1/EPsOzf/RkxIupoJB8pgUtOuPZaWPbW0fQ6bjSnnRZ1NMHpddxDbHPQ84wc6Qt/RIKkBC6RGzcOysqg+68m0uv40VGHEygz6H3a/1JcDGedBW+9FXVEkk+UwCVSU6fC6afDgQfCLmfdkJf9xK1aVzNhAvTq5a8TnzcvmOXqMkFRApfIfPEFDB7sR7l58klo1aaZd4SKga228rcBqKnx/furV0cdkeQDJXCJRHU1nHyyL3Z58snMRs+Ji969/SDLn3ziv3VoZB/JlBK4ROKqq2DyZH9nwT59oo4me4qL/Z0Un3rKn7gVyYQSuGRF3eHWuv+9kGufq+Ccc/jhipNsD8eWLK5srffCC6GkBK68Ev77wX6hr0/ylxK4hK7+cGsr1i/ABpVywLkVSV8Pezi2VHFla71mMGqUv0XA3Huu5LtlPwl1fZK/lMBboGxfvZBsuDXXei0j3ixL+XpYw6I1Flc21gvQoQM88QRgNXx42zWsWRP6KptEV7jEgxK4hK6x4dayMRxbU5Yf9npr7bgj7PaHkXy7pBelpRpbU5pOCVxC19iwZ1ENixbVeuvaao/pFA5+gLFjfbeKSFMogUvoft+7HNanHvYsqmHR0l3v5DMmhzpQRM8Bj9K/P1xwAcyYEdpqJA8pgcdUXPooV6+GB/9cwpZTRtHWNgM2HfYsqmHRcmE4NgBr5Xj4YejeHX77W1i5MqurlxjT7WQlNM7B2WfDp5/Ca/eWcNW81MOeRTUsWi4Mxwa+kOnxx/34n2ed5Qt+8vG2AhIsHYFLaG69Ff71L1+wcvDBUUeT+/bfH264wVem3nJL1NFIHCiBS9rSKXqp7S+eOhUuuQQGDfI/ZVPJ9udFF8Hxx8Nf/uJv9BXksiX/KIFLWppS9LJiBQwZAj17wujR6gpIJtX+fHR2BQ8+6PfdkCHw5ZfBLVtJPP8ogUta0i162bDBl4mvWAHjx0PnzsHFEJcTt+loaH927uy7npYv97caaOpNr6IsUJLsUgKXtKRb9DJyJEycCHfcAfvsk43I4qmx/dm3r+8Hf/FFuO66YJct+UMJXNKSTtHLSy/BNdfAmWfCsGHZiiye0tmf554Lp5zib3r1+uvBLlvygxK4pKWxopeFC33XSZ8+cPvt2Y8v7GKboKVTRHTomGIWHnQ0O+/s752+dGlwy5b8oAQuaWmo6GXdOl+AUl3t+707dGhkYZJ2EVFB++8YPx6++cYn8erq4JYt8adCHklbqqKXiy+Gd9/1xSe9e0cUXAylW0S0xx5+4IvTTvMDYfz978EtW+JNR+CSkcce810mF1/sx7eUcJx6Kpxzji+KeuGFqKORXKEEHkOZFGkEWeAxZ44vlT/wQA0PlkzQxTS33gp77+2T+a2vh1eooyKg+FACj5lMijSCLPD49ls48UTYfHMYNw7atGnyIvJaGMU07dv768PX/rSCP70WTqGOioDiRQk8ZjIp0giqwMM5/3V+7lx49FHYbrsmvb1FCKuY5mc/g47Hl1HTOpxCHRUBxYtOYsZMJkUaQRV4fP7qYN4Y67tN+vVr0ltbjDCLaf67PrxlqwgoXnQEHjOZFGkEUeCx6pM9mTfuPAYNgksvTfttLU6YxTRxXbYETwk8ZjIp0si0wOOLL+CjO6+mfdcvGDNGN6lqSJjFNMmWvVnr8JatIqDcpQQeM5kUaWTy3upqGDoUqtduwe7Dr6JTp8y2I9+FWUzzo2U7YGUv9l06ilP2DHjZqAgo16kPPIYyKdJo7nsvvxymTIFdz7mRjjvMa9I6W6owi2nqLvugTyfzt7/BPXvC738f7LJVBJTbdAQujRo/Hv75TzjvPOhxwMSow5F6RoyAAQPgj3+EN9+MOhrJpowSuJnNN7NZZjbTzKYFFZTkjlmz4Iwz4Je/hJtuijoaSaagACoqYKed/LX5ixZFHZFkSxBH4Ic65/Z2zhUFsCzJIV9+CccdB1tuCU88Ae3aRR2RpNK5Mzz1FHz3HZxwgv/ZkHwaHKMlUx+4JLVunT+aW7oU3ngjmGId9aeGa7fd/JH4ccdBaSm4w3SlUL7L9AjcAa+Y2XQzKw0iIMkNF17oT1refz/st1/U0Ui6Bg70oyI98gh8/sqJUYcjIcv0CPzXzrklZtYdmGhmc51zb9SdIZHYSwF69lQxQBzcdZe/femll/pBGrJFR+jBKCuDmTPhiXF/YLNtVUGZzzI6AnfOLUn8XA48CWxyrOacG+WcK3LOFXXr1i2T1UmaMunffO01uOACOOYYKFftRiy1agVjxkDHHT5lzl0jeP/9qCOSsDQ7gZvZ5ma2Re1z4EhgdlCBSfbNnu1PgO2yi79JVUFB1BFJc3XsCHtedAUFm33LMcfA559HHZGEIZMj8B7AW2b2PvAu8Lxz7qVgwpJs+/xzOPpo/8F/8UV/5UlLEbfxNNPVrsuX9LnoMlatgmOP9cOySX5pdh+4c24esFeAsUhEVq/2XSYrV/pCkB12SD1vPia6fFX7t3ppP5/Ahw6Fp5+G1hFfe1bbvae2lDlVYrZw69f7AYlnz/YVl3vvHXVEErT+/eGOO/xQbOef7+/nLvlB14HHVBBHLzU1MGwYvPIK3HcfHHVU5nGFpXaYr6oNVRTeXEh5v/K8vsFS0Nt77rmwcKEfELlrV2Dn4GKV6OgIvIVyzt874+GH4ZprfCKvlWtVei1tmK+wtvdvf/OJ/NprYdGLJwURatblWtuMWotP4C21QXw24WzuvBMuucRfN5zLWtowX2Ftr5nvSjnpJJj3+B9YOuWYjJYn0YtFF4pOegRr4fMns+j531FaCjfckPvl1vk8zFeyNh3m9hYUwEMPwSsfvcN/xvyZx4+GIUMyXqxEpMUfgbc0t9wCn40/l5NPhjvvzP3kDS1vmK+wt7dtW9jj/Kvo1Hs2p5wC48ZtOk++XlqZb5TAW5B//AMuugh+8xtfqReXQp2WNsxXNra3oF0Vff50GQccAKecAmPHBrZoySIl8BaivBz++lff/zl2LLRpE3VE6Wtpw3xla3sL2n/HCy/AQQfB737n72Qo8RKLPnBpPuf8iC0jR/oP6YMPRl/I0RwtbZivbG1vx47w/PP+LoannQbff//jK5Ikt+kIPI9t2OALN0aOhDPPhNGj45m8JVybbw7PPQdHHglnn+2vFVexTzwogeepDevaMmSIvzXsX//qC3Xi0uct2dehAzzzjP+WVlbm7wdfUxN1VNKYnE/gtRVpUxZMofDmwkCLN8JcdlQqZlXw9qJ3eWvJqzzRs5DfXV/B9df7W4xKfmnulSKp2n2bNv7k9p//DLfdBief3PjQbEGtO+z35quc/liHWYGXj9V9FbMqOPvpUqpqvgNz0HkBT1THe5skWI21+1at4MYb/RVLjz8OhxwS3K1oM/nMBfF5zceivZxO4GFW4OVjdd+fnivj+w35tU0SrHTb/SWX+EGSP/oI9t0X3n03e+sO+r35LKcTeJgVaflU3VdT4+9vsaIqf7ZJwtGUdj9oEEydCu3awcEH+yuYMjm5mclnLp8+r0HK6QQeZkVavlT3LV8OAwbAFVdAh+r82KZ8kmsVjU1t9336+KPvAw6As86C00+HNWua1x2RyWcuXz6vQcvpBB5mRVo+VPdNmgR77eVHj7/7bhg1JP7bJOFqTrvv1g0mTvT1BI88Ar/4BaxZ+NOsrDuI9+aznE7gYVakxbm6b+1a30d5xBHQpYs/Qjr3XCj5eXy3SbKjue2+oACuvtofNHzzDcy45i4WPncK1dXhrzvT9+aznC/rCLMiLZer+1LdgfGNN3ylXGWlT9o33ugLMWrl8jZJbsikjRx6KMycCXscM5XPJpSy//6+b3zPPcNft9r2pnL6CFw2+vprX1V5yCG+wnLSJN9tUjd5i2RD9+6wx/lXs/t5VzN/PvTtC1dd5b8ZZiofL/ULkxJ4jtuwwVdT9u7tE/aFF8KsWXDYYeGsT8USkq5u+07ho4/8mKrXXAO77gr/+lduluHna7tWAs9RzsFXs/alb1847zz/FXXGDLj55vCOuvOxuEnC1a2bv4vhlCn+fMyQIf7g4u23o4kn2RF8PrdrJfAc9NprMPPa25h10z9Ytcof1bz+ur/iJEwqlpDmOvhgmD4dbr8dPvwQ9t8fjj0W/u//oo4s/XYdx+6bWCTwXLuWNgw1Nf5mQoccAv36wfdfbkPv027iP/+BE0/Mzsg5KpaQTLRu7c/TzJvn72j473/7/vGBA2Hy5Oi6VoJo17ma3GORwMMU9T+HNWv8QLO77uor3z77zA979svrS9ju0Gdo2zZ7sahYQoLQsSNcfrlvyyNG+O6UQw+FoiLf3VKzPouNmvxu1y0+gUfBOX854JlnwjbbwPDhsNVW8Nhj8OmncMEF0KrNuqzHpWIJCVLnzv7a8YUL4Z574Ntv/e1qp/5pPJUVf2TmzOzEkc/tOuevA88Xzvk+wief9Hd5q6z0RypDh/qb6P/qV1FHyA9FEcOeHkbVhip6depFeb/yFl8sIU1Tv4Zhs82gtNS380mToOTK91gyeSD77OPP65x4IgweDLvv3viym/NtOZ/btRJ4iNas8UfaL73k7+y2aJGvaCsu9tfNDh6ce9dxq1hCwtKqla8e3v3za1i/ZgtO5lkefRSuvNI/dt0Vvt3pHLrsMY3vv4f27YNbd762ayXwAK1aBdOm+aQ9aRK88w5UV/uGeOSR/lrZY4+Frl2jjlQkOnUT6PDhsGSJ/2Y6YQK8/tJQFr1Qwla3+8GW+/XzV7T07Zt7Bzu5QAm8GZyDZctgzhxfVPPee/7x8cf+9Vat/A1/LrnEN8Bf/9p/jawrVam8SEuz3Xb+6pXzz4cD7xrIqo/3oh/X8uqrcOmlfp6CAl8Lsd9+PpnvsYfvcmnpB0NK4CmsX+9HIlm4cOPj00990p4zB1au3Djvttv6hnXqqf7nvvv6EzjFo4uZuhgmbza5SeuurRqr2lBF4c2FedNfJ/khrPZZMauCaf99jaouL/JNp7GUjy3niB4lvPee/zb77rswfjzce+/G9/To4RP5zjvDV9tX8O+qmVQXrGL7Gwu57vByfrdXenE1tE25/HnM+wReVeXvnrZmzY9/fvMNfPklrFix8Wft82XLYOnSTa9b7dEDdtvNjxW4225w38KL6bDdfKb+aUJg8aaqGgPSbjQ6qs8Pufh3DKJ9NmW5owZCyTElHHOMn885fy7po498wVDto+KDCtZ0LYW2vmDn8zULOPXxUv54Aey0poTu3WHuN5fRZsuv+eeX/gCrUyf/mLqmgus+Sr5NQCjbG5RYJPBHHoFXX4V16zY+qqo2Pp+1ZC411W3YfvOf/jBt7VqfrNevb3z5W2wBW2/ty4K3286fGe/Z0z9umXsx7bsuY8qFj2zSDTJh9PTAt7WhqrFcaDDSsoXVPtNdrtnGz2b//hvnLby5jDWr6t1Nq+1a1h9cxjYzSli+HFbO24d1q7vwlxfrrfyiMui86bpPHV2GGdRsselrZ48tY1xZCa1b+wKmyQsnYa02MGCXI2nd2nf51L5W+/tZZ8EuuzR7FyUViwQ+Z44vJW/b1g/v1LbtxkeHDtBmy6+x1tXs9bOf/vB6+/Y+MT/2yb0UtPuOssMvoGNHP632Z9euPnE3dLb74USSrp+8w6JqSMllYbXPTJebar61rRfy/PP+efHok3AOnv3NZFauhNWr/YUHB726kGRFom6L5NMBvm+7kMWL/UUK1dWw5qudcTWteO3zjdM2bNj4vLoaDj+8hSbw8nL/SKV49OUAjEvylXPqaH/DmjPPuCCM0ALXs1NPFqxakHR6S5eLXQotTVjtM9PlpvP+uu1nyy3rzPNe8vf26uzfm+q1GTM2/l48+tRN1pENqsTMMflcNSa5o7m3kAirfWa63LCGa8v1z2MsjsBbknyuGpP4C6t9ZrrcTN6fzntz9fOoBJ6D8rVqTPJDWO0z0+WGNVxbLn8e1YUiIhJTGSVwM+tvZh+bWaWZXRZUUE3R0FBJmQ6jlMmyw4xLJFe1xLYd5TY3uwvFzAqAO4AjgMXAe2b2jHPuo6CCa0xDRQWQ2QX4mSw7zLhEclVYRT65LOptzuQIfD+g0jk3zzm3DngMGBRMWOlp6OL/TIcHy2TZYcYlkqtaYtuOepszOYn5E2BRnd8XA7+sP5OZlQKlAD17Bnstc3Mu/s+0MCCdZYcZl0iuyucitFQnL6Pe5kyOwJON0rhJ4ZJzbpRzrsg5V9StW7cMVrephoZKynQYpUyWHWZcIrkq7m27OdfGR73NmSTwxcAOdX7fHliSWThNE+YF+JksO86FARD9OKEST7nQtrPddqPe5ky6UN4DepvZjsDnwFDglECiSlOYF+Bnsuw4FwaINFdLLEKLepvN1b9nalPebDYAuBkoAB5wzjX4b6eoqMhNmzat2etLpaHBETIdOCGTZYf1XpGohfmZi6Owt9nMpjvniupPz6gS0zn3AvBCJsuQ5FpS4xeR5lElpohITCmBi4jElG5mJSJNpi6+3KAjcBGRmMqLI/CGjgZ0pCAiYYsqz+RFAg+T/jmISK5SAg+REryIhEkJPCJK7iKSKZ3EFBGJKSVwEZGYUheKiARK3YPZoyNwEZGYUgIXEYkpJXARkZhSAhcRiSklcBGRmFICFxGJKSVwEZGYUgIXEYkpJXARkZjKaFT6Jq/MbAWwoJlv3xr4MsBwgqK4mkZxNU2uxgW5G1s+xtXLOdet/sSsJvBMmNk051xR1HHUp7iaRnE1Ta7GBbkbW0uKS10oIiIxpQQuIhJTcUrgo6IOIAXF1TSKq2lyNS7I3dhaTFyx6QMXEZEfi9MRuIiI1JFTCdzMfmtmH5pZjZkV1XvtcjOrNLOPzeyoFO/f0czeMbNPzGycmbUNIcZxZjYz8ZhvZjNTzDffzGYl5psWdBxJ1jfCzD6vE9uAFPP1T+zDSjO7LAtx/cPM5prZB2b2pJl1TjFfVvZXY9tvZu0Sf+PKRFsqDCuWOuvcwcxeN7M5ifZ/YZJ5is1sVZ2/71Vhx5VYb4N/F/NuTeyvD8ysbxZi2qXOfphpZqvN7KJ682Rtf5nZA2a23Mxm15m2lZlNTOSiiWbWJcV7T0/M84mZnd7klTvncuYB7AbsAkwGiupM3x14H2gH7Ah8ChQkef/jwNDE87uBP4Qc743AVSlemw9sncV9NwK4pJF5ChL7biegbWKf7h5yXEcCrRPPrweuj2p/pbP9wHnA3YnnQ4FxWfjbbQv0TTzfAvhPkriKgeey1Z7S/bsAA4AXAQN+BbyT5fgKgC/w10lHsr+Ag4G+wOw6024ALks8vyxZuwe2AuYlfnZJPO/SlHXn1BG4c26Oc+7jJC8NAh5zzlU55z4DKoH96s5gZgYcBoxPTBoDHB9WrIn1DQHGhrWOEOwHVDrn5jnn1gGP4fdtaJxzrzjnqhO/vg1sH+b6GpHO9g/Ctx3wbalf4m8dGufcUufcjMTzb4A5wE/CXGeABgEPOe9toLOZbZvF9fcDPnXONbdAMGPOuTeAr+pNrtuOUuWio4CJzrmvnHNfAxOB/k1Zd04l8Ab8BFhU5/fFbNrAuwIr6ySLZPME6SBgmXPukxSvO+AVM5tuZqUhxlHX8MTX2AdSfJNSwNEAAANESURBVGVLZz+G6Sz80Voy2dhf6Wz/D/Mk2tIqfNvKikSXzT7AO0le3t/M3jezF81sjyyF1NjfJeo2NZTUB1FR7K9aPZxzS8H/gwa6J5kn432X9UGNzexVYJskL5U5555O9bYk0+pfPpPOPGlJM8aTafjo+9fOuSVm1h2YaGZzE/+pm62huIC7gGvw23wNvnvnrPqLSPLejC9DSmd/mVkZUA1UpFhM4PsrWahJpoXWjprKzDoCE4CLnHOr6708A99NsCZxfuMpoHcWwmrs7xLl/moLHAdcnuTlqPZXU2S877KewJ1zhzfjbYuBHer8vj2wpN48X+K/vrVOHDklmyeQGM2sNTAY+EUDy1iS+LnczJ7Ef33PKCGlu+/M7F7guSQvpbMfA48rcXLmWKCfS3T+JVlG4PsriXS2v3aexYm/cyc2/XocODNrg0/eFc65J+q/XjehO+deMLM7zWxr51yo9/xI4+8SSptK09HADOfcsvovRLW/6lhmZts655YmupSWJ5lnMb6vvtb2+PN/aYtLF8ozwNDEFQI74v+Tvlt3hkRieB04MTHpdCDVEX2mDgfmOucWJ3vRzDY3sy1qn+NP5M1ONm9Q6vU7npBife8Bvc1frdMW//XzmZDj6g9cChznnFubYp5s7a90tv8ZfNsB35ZeS/VPJyiJPvb7gTnOuZtSzLNNbV+8me2H/+z+N+S40vm7PAOclrga5VfAqtqugyxI+S04iv1VT912lCoXvQwcaWZdEl2eRyampS8bZ2mbcDb3BPx/pSpgGfByndfK8FcQfAwcXWf6C8B2iec74RN7JfAvoF1IcY4Gfl9v2nbAC3XieD/x+BDflRD2vnsYmAV8kGg829aPK/H7APxVDp9mKa5KfD/fzMTj7vpxZXN/Jdt+YCT+HwxA+0TbqUy0pZ2ysI8OxH91/qDOfhoA/L62nQHDE/vmffzJ4AOyEFfSv0u9uAy4I7E/Z1Hn6rGQY+uAT8id6kyLZH/h/4ksBdYn8tcw/HmTScAniZ9bJeYtAu6r896zEm2tEjizqetWJaaISEzFpQtFRETqUQIXEYkpJXARkZhSAhcRiSklcBGRmFICFxGJKSVwEZGYUgIXEYmp/wdeJGpexjxZmgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.optimize import minimize\n",
    "p0=[8,1,3,0]\n",
    "# default BFGS algorithm, numerical derivatives\n",
    "out=minimize(logl,p0,args=(x,y))\n",
    "print(out)\n",
    "\n",
    "xfit=np.arange(-10,10,0.01)\n",
    "plt.errorbar(x,y,sig,fmt='go')\n",
    "plt.plot(xfit,model(xfit,*out.x),color='b')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an example of specifying bounds on a parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      fun: -135.0276207289479\n",
      " hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>\n",
      "      jac: array([ 0.00000000e+00, -1.13686838e-05, -1.70530258e-05,  0.00000000e+00])\n",
      "  message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
      "     nfev: 85\n",
      "      nit: 15\n",
      "     njev: 17\n",
      "   status: 0\n",
      "  success: True\n",
      "        x: array([ 9.8247912 , -0.22915942,  2.7663895 ,  0.91959721])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1235ffd10>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD5CAYAAAA+0W6bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3gU1f3H8feXqyAKqIBahWhLvVIVo63WSxQviCJKrYKpoqKxVaq22nrJo/LTpl5a/Xm/4A3UiFjwfkcU1D6gAj8UFKwRuSNglZtoIOT8/jgbiWE32WRndnY2n9fz7JPN7OzMdyZnv5k9M9855pxDRETip0XUAYiISNMogYuIxJQSuIhITCmBi4jElBK4iEhMKYGLiMRUq4ZmMLOdgUeB7YFqYIRz7nYz2wYYAxQA84BTnXPf1Les7bbbzhUUFGQYsohI8zJt2rSvnHNd6k63hq4DN7MdgB2cc9PNbCtgGnAScBbwtXPuRjO7AujsnLu8vmUVFha6qVOnNnUbRESaJTOb5pwrrDu9wS4U59xS59z0xPM1wGzgJ8AAYFRitlH4pC4iIlnSqD5wMysA9gPeA7o555aCT/JA16CDExGR1NJO4GbWARgHXOKcW92I95WY2VQzm7pixYqmxCgiIkmklcDNrDU+eZc7555OTF6W6B+v6Sdfnuy9zrkRzrlC51xhly6b9cGLiEgTNZjAzcyAh4DZzrlba730PDAk8XwI8Fzw4YmISCoNXkYI/Bo4A5hpZjMS064CbgSeMrOhwALgt+GEKCIiyTSYwJ1z7wKW4uU+wYYjIiLpUiWmiEhMKYFLs1A0soiikUVRhyESKCVwEZGYUgIXEYkpJXARkZhSAhcRiSklcBGRmFICFxGJKSVwEZGYUgIXEYkpJXARkZhSAhfJkKo8JSpK4CIiMaUELiISU0rgIiIxpQQuIhJTSuAiIjGlBC4iElNK4CIiMaUELiISU0rgIiIxpQQuIhJTSuAiIjGlBC4iElNK4CIiMaUELiISU0rgIiIxpQQuIhJTSuAiIjGlBC55r3xmOVMWTWHS/EkU3FZA+czyqEMSCYQSuOS18pnllLxQQuXGSgDmr5pPyQslSuKSF5TAJa+VTihl3YZ1P5q2bsM6SieURhSRSHCUwCWvLVi1oFHTReJECVzyWveO3Rs1XSROlMAlr5X1KaN96/Y/mta+dXvK+pRFFJFIcJTAJa8V9ypmRP8RtG3ZFoAeHXswov8IinsVRxyZSOZaRR2ASNiKexXzwLQHAJh41sRogxEJkI7ARURiqsEEbmYPm9lyM5tVa9pwM1tsZjMSj37hhikiInWlcwQ+EuibZPr/Ouf2TTxeDjYskXhQladEqcEE7px7G/g6C7GIxIqqPCVqmfSBDzOzjxJdLJ0Di0gkJlTlKVFragK/F/gpsC+wFLgl1YxmVmJmU81s6ooVK5q4OpHcoypPiVqTErhzbplzbqNzrhp4ADiwnnlHOOcKnXOFXbp0aWqcIjlHVZ4StSYlcDPbodavJwOzUs0rkq9U5SlRa7CQx8xGA0XAdma2CLgWKDKzfQEHzAPODzFGkZxUU8059LmhVG6spEfHHpT1KVOVp2RNgwncOTc4yeSHQohFJHZU5SlRUiWmiEhMKYGLNKBoZBFFI4uiDkNkM0rgIiIxpQQuIhJTSuAiIjGlBC4iElNK4CIiMaUELiISU0rgIiIxpQQuIhJTSuAiIjGlBC55QdWS0hwpgYuIxJQSuIhITCmBi4jElBK4iEhMNTigg0g+0GALko90BC6B0tUgItmjBC4iElNK4CIiMaUELiISU0rgIiIxpQQuIhJTSuAiIjGlBC4iElMq5BHJkIqEJCo6AhcRiSklcBGRmFICFxGJKSVwafbqu39L+cxypiyawqT5kyi4rYDymeXZDU6kHkrgIimUzyyn5IUSKjdWAjB/1XxKXihREpecoQQukkLphFLWbVj3o2nrNqyjdEJpRBGJ/JgSuOSMXLsV7YJVCxo1XSTblMBFUujesXujpotkmwp5JGtqjq6zWfhSXQ1Tp8I778AHH0BFBXz5JaxeDS1bQvv2sLbdXWzRdQl3fwuHHgq9eoEZlPUpo+SFkh91o7Rv3Z6yPmVZi1+kPkrgkpc+/hjuuw/GjYOlS/207t1hzz3hF7+Ajh1h40ZYuxaemVLJytn7MWyYn2+XXWDQIDj33GJG9Iehzw2lcmMlPTr2oKxPGcW9iqPbMJFalMAlr7z9NvzP/8Cbb0LbtnDCCXDyyXDUUdCtW/L3zB15Kc7BqCMmMmECPPUU3Hwz3HQTDBxYzF57TWSrgs9UMi85Rwlc8sK3S7ozd8wfOPwj2HFHuOEGOPdc2G679N5vBgUFMHSofyxZAnfeCffeC6vGPkC3g19l0VGw006hboZIozR4EtPMHjaz5WY2q9a0bcxsvJl9lvjZOdwwRZJbvx6uuw6mXfsgqyv25qabfD/3FVekn7yTqfknMH8+7NzvCZa/fyS77+4TenV1cPGLZCKdq1BGAn3rTLsCmOCc6wlMSPwuzVy2qxZnz4bCQrh2bDl26c5U/XVr7mlTwNMVwa23Y0fY9bcjOODvZ3LwwXDBBXD00bBwYWCrEGmyBhO4c+5t4Os6kwcAoxLPRwEnBRyXxEy2qxYffxwOOADmbVVO21NKqG63PNT1tuvyJa+9BvffD++/D717wxtvBLoKkUZr6nXg3ZxzSwESP7sGF5LEUbaqFquqYNgwOOMMn0S3HlhKpctOtaQZlJT4yxK7dYNjj/XdLM4FviqRtIReyGNmJWY21cymrlixIuzVSUSyUbW4Zg0MGAB33w2XXuqvNFmyNvvVkrvtBlOmwGmnwVVX+ZOeGzY0bVm5Vn0q8dLUBL7MzHYASPxcnmpG59wI51yhc66wS5cuTVyd5LqwqxaXLoXDDoPXXvPXd//zn9CqVXTVkh06QHk5XHstPPII9O/v/8GIZFNTE/jzwJDE8yHAc8GEI3FV1qeM9q3b/2haUFWLCxf65F1RAS+9BOefn531NsQMhg+HBx/0/eFHHgnffBP6akV+kM5lhKOBycBuZrbIzIYCNwJHm9lnwNGJ36UZK+5VzIj+I2jbsi0APTr2YET/ERlXLX7xhU/eK1bA+PG+3zkb622MoUPh2Wfho498wdDXdU/5i4SkwUIe59zgFC/1CTgWibniXsU8MO0BIJj7nVRUwBFHwLffwoQJsP/+2VlvU5xwgk/iJ58Mffr4I/Jtt40kFGlGdDdCaZRsnXRbtMgfzX73Hbz1VurknUuOOw6eew7mzPGxr1oVbTw6QZr/lMAl56xY4Ytlvv7an7TcZ5+oI0rfscfCM8/4m2mdeKL/ByQSFiVwySlV37XnuONg3jx44YV4HHnX1bcvPPqov4Xtaaf5a9dFwqCbWUnOqK5qxcd3/I21Fb4/+fDDo46o6QYN8t8gLrwQdip6laXv9MUs6qgk3+gIXHKCc/CfUZeyck5vHnoIjj8+6ogyd8EF0OOkh1n2774MHx51NJKPlMAlJ9xwAyx79zh6nDiSM8+MOprg9DjxUbY/9CWuu84X/ogESQlcIjdmDJSWQtdfjafHSSOjDidQZtDzzP+lqAjOOQfefTfqiCSfKIFLpCZPhiFD4JBDYLdzbs7LfuIWraoYNw569PDXic+dG8xydZmgKIFLZL78EgYO9KPcPPMMtGjdxDtCxcA22/jbAFRX+/791aujjkjygRK4RKKqCgYP9sUuzzyT2eg5cdGzpx9k+bPP/LcOjewjmVICl0hccw1MnOjvLNirV9TRZE9Rkb+T4rPP+hO3IplQApesqD3cWte/F3DDi+Wcdx4/XHGS7eHYksWVrfVefDEUF8PVV8N/Pzow9PVJ/lICl9DVHW5txYb52IASDj6/POnrYQ/HliqubK3XDEaM8LcImHP/1Xy37Cehrk/ylxJ4M5TtqxeSDbfmWq1j+DulKV8Pa1i0huLKxnoB2reHp58GrJqP77yetWtDX2Wj6AqXeFACl9A1NNxaNoZja8zyw15vjV12gT3+cB3fLulBSYnG1pTGUwKX0DU07FlUw6JFtd7attlrGgUDH2b0aN+tItIYSuASut/3LIMNqYc9i2pYtHTXO/GsiaEOFNG93xP07QsXXQTTp4e2GslDSuAxFZc+ytWr4ZE/F7P1pBG0sXbA5sOeRTUsWi4MxwZgLRyPPQZdu8JvfwsrV2Z19RJjup2shMY5OPdc+PxzePOBYq6Zm3rYs6iGRcuF4djAFzI99ZQf//Occ3zBTz7eVkCCpSNwCc0dd8C//uULVg47LOpoct9BB8HNN/vK1NtvjzoaiQMlcElbOkUvNf3FkyfDZZfBgAH+p2wu2f685BI46ST4y1/8jb6CXLbkHyVwSUtjil5WrIBTT4Xu3WHkSHUFJJNqfz4xq5xHHoGdd/b78Kuvglu2knj+UQKXtKRb9LJxoy8TX7ECxo6FTp2CiyEuJ27TUd/+7NTJdz0tX+5vNdDYm15FWaAk2aUELmlJt+jluutg/Hi4+27Yb79sRBZPDe3P/feH226DV16BG28MdtmSP5TAJS3pFL28+ipcfz2cfTYMHZqtyOIpnf35+9/7W+5efTW89Vawy5b8oAQuaWmo6GXBAt910qsX3HVX9uMLu9gmaOkUER0xqogFhx5Hz54+kS9dGtyyJT8ogUta6it6Wb/eF6BUVfl+7/btG1iYpF1E1Krdd4wd6wuiBg/2+zioZUv8qZBH0paq6OXSS+H9933xSc+eEQUXQ+kWEe29tx/4YsgQPxDG3/8e3LIl3nQELhl58knfZXLppX58SwnHmWf6qtYbboCXX446GskVSuAxlEmRRpAFHrNn+6RyyCEaHiyZoItp7rgD9t0XzjgD7ngrvEIdFQHFhxJ4zGRSpBFkgce338Ipp8CWW8KYMdC6daMXkdfCKKZp185fH77up+X86c1wCnVUBBQvSuAxk0mRRlAFHs7BeefBnDnwxBOw446NenuzEFYxzc9+Bh1OKqW6VTiFOioCihedxIyZTIo0girwWPzGQN4e7btN+vRp1FubjTCLaf67IbxlqwgoXnQEHjOZFGkEUeCx6rO9mTvmAgYMgMsvT/ttzU6YxTRxXbYETwk8ZjIp0si0wOPLL+GTe65li22/ZNQo3aSqPmEW0yRbdrtW4S1bRUC5Swk8ZjIp0sjkvVVVMGgQVK3bij2HXUPHjpltR74Ls5jmR8t2wMoeHLB0BKfvHfCyURFQrlMfeAxlUqTR1PdeeSVMmgS7n3cLHXae26h1NldhFtPUXvahn0/kb3+D+/f2908JctkqAsptOgKXBo0dC//8J1xwAXQ7eHzU4Ugdw4dDv37wxz/CO+9EHY1kU0YJ3MzmmdlMM5thZlODCkpyx8yZcNZZ8Mtfwq23Rh2NJNOyJZSXw667+mvzFy6MOiLJliCOwI9wzu3rnCsMYFmSQ776Ck48EbbeGp5+Gtq2jToiSaVTJ3j2WfjuOzj5ZP+zPvk0OEZzpj5wSWr9en80t3QpvP12MMU66k8N1x57+CPxE0+EkhJwR+pKoXyX6RG4A143s2lmVhJEQJIbLr7Yn7R86CE48MCoo5F09e/vR0V6/HFY/PopUYcjIcv0CPzXzrklZtYVGG9mc5xzb9eeIZHYSwC6d1cxQBzce6+/fenll/tBGrJFR+jBKC2FGTPg6TF/oN0OqqDMZxkdgTvnliR+LgeeATY7VnPOjXDOFTrnCrt06ZLJ6iRNmfRvvvkmXHQRHH88lKl2I5ZatIBRo6DDzp8z+97hfPhh1BFJWJqcwM1sSzPbquY5cAwwK6jAJPtmzfInwHbbzd+kqmXLqCOSpurQAfa+5CpatvuW44+HxYujjkjCkMkReDfgXTP7EHgfeMk592owYUm2LV4Mxx3nP/ivvOKvPGku4jaeZrradv6KXpdcwapVcMIJsGZN1BFJ0JrcB+6cmwvsE2AsEpHVq32XycqVvhBk551Tz5uPiS5f1fytXj3QJ/BBg+C556BVxNee1XTvqS1lTpWYzdyGDX5A4lmzfMXlvvtGHZEErW9fuPtuPxTbhRf6+7lLftB14DEVxNFLdTUMHQqvvw4PPgjHHpt5XGGpGearcmMlBbcVUNanLK9vsBT09p5/PixY4AdE3nZb4OfBxSrR0RF4M+Wcv3fGY4/B9df7RF4j16r0mtswX2Ft79/+5hP5DTfAwldOCyLUrMu1thm1Zp/Am2uD+GLcudxzD1x2mb9uOJc1t2G+wtpeM9+VctppMPepP7B00vEZLU+iF4suFJ30CNaClwaz8KXfUVICN9+c++XW+TzMV7I2Heb2tmwJjz4Kr3/yHv8Z9WeeOg5OPTXjxUpEmv0ReHNz++3wxdjzGTwY7rkn95M3NL9hvsLe3jZtYK8Lr6Fjz1mcfjqMGbP5PPl6aWW+UQJvRv7xD7jkEvjNb3ylXlwKdZrbMF/Z2N6WbSvp9acrOPhgOP10GD06sEVLFimBNxNlZfDXv/r+z9GjoXXrqCNKX3Mb5itb29tyi+94+WU49FD43e/8nQwlXmLRBy5N55wfseW66/yH9JFHoi/kaIrmNsxXtra3Qwd46SV/F8Mzz4Tvv//xFUmS23QEnsc2bvSFG9ddB2efDSNHxjN5S7i23BJefBGOOQbOPddfK65in3hQAs9TG9e34dRT/a1h//pXX6gTlz5vyb727eH55/23tNJSfz/46uqoo5KG5HwCr6lImzR/EgW3FQRavBHmsqNSPrOcKQvf590lb/B09wJ+d1M5N93kbzEq+aWpV4qkavetW/uT23/+M9x5Jwwe3PDQbEGtO+z35quc/liHWYGXj9V95TPLOfe5EiqrvwNz0Gk+T1fFe5skWA21+xYt4JZb/BVLTz0Fhx8e3K1oM/nMBfF5zceivZxO4GFW4OVjdd+fXizl+435tU0SrHTb/WWX+UGSP/kEDjgA3n8/e+sO+r35LKcTeJgVaflU3Vdd7e9vsaIyf7ZJwtGYdj9gAEyeDG3bwmGH+SuYMjm5mclnLp8+r0HK6QQeZkVavlT3LV8O/frBVVdB+6r82KZ8kmsVjY1t9716+aPvgw+Gc86BIUNg7dqmdUdk8pnLl89r0HI6gYdZkZYP1X0TJsA++/jR4++7D0acGv9tknA1pd136QLjx/t6gscfh/33h7ULfpqVdQfx3nyW0wk8zIq0OFf3rVvn+yiPPho6d/ZHSOefD8W/iO82SXY0td23bAnXXusPGtasgenX38uCF0+nqir8dWf63nyW82UdYVak5XJ1X6o7ML79tq+Uq6jwSfuWW3whRo1c3ibJDZm0kSOOgBkzYK/jJ/PFuBIOOsj3je+9d/jrVtveXE4fgcsm33zjqyoPP9xXWE6Y4LtNaidvkWzo2hX2uvBa9rzgWubNg9694Zpr/DfDTOXjpX5hUgLPcRs3+mrKnj19wr74Ypg5E448Mpz1qVhC0tXlgEl88okfU/X662H33eFf/8rNMvx8bddK4DnKOfh65gH07g0XXOC/ok6fDrfdFt5Rdz4WN0m4unTxdzGcNMmfjzn1VH9wMWVKNPEkO4LP53atBJ6D3nwTZtxwJzNv/QerVvmjmrfe8lechEnFEtJUhx0G06bBXXfBxx/DQQfBCSfA//1f1JGl367j2H0TiwSea9fShqG62t9M6PDDoU8f+P6r7el55q385z9wyinZGTlHxRKSiVat/HmauXP9HQ3//W/fP96/P0ycGF3XShDtOleTeywSeJii/uewdq0faHb33X3l2xdf+GHPfnlTMTse8Txt2mQvFhVLSBA6dIArr/Rtefhw351yxBFQWOi7W6o3ZLFRk9/tutkn8Cg45y8HPPts2H57GDYMttkGnnwSPv8cLroIWrRen/W4VCwhQerUyV87vmAB3H8/fPutv13t5D+NpaL8j8yYkZ048rld5/x14PnCOd9H+Mwz/i5vFRX+SGXQIH8T/V/9KuoI+aEoYuhzQ6ncWEmPjj0o61PW7IslpHHq1jC0awclJb6dT5gAxVd/wJKJ/dlvP39e55RTYOBA2HPPhpfdlG/L+dyulcBDtHatP9J+9VV/Z7eFC31FW1GRv2524MDcu45bxRISlhYtfPXwnouvZ8ParRjMCzzxBFx9tX/svjt8u+t5dN5rKt9/D1tsEdy687VdK4EHaNUqmDrVJ+0JE+C996CqyjfEY47x18qecAJsu23UkYpEp3YCHTYMlizx30zHjYO3Xh3EwpeL2eYuP9hynz7+ipbevXPvYCcXKIE3gXOwbBnMnu2Laj74wD8+/dS/3qKFv+HPZZf5BvjrX/uvkbWlKpUXaW523NFfvXLhhXDIvf1Z9ek+9OEG3ngDLr/cz9Oypa+FOPBAn8z32st3uTT3gyEl8BQ2bPAjkSxYsOnx+ec+ac+eDStXbpp3hx18wzrjDP/zgAP8CZyikUVMXgQT201s1LprqsYqN1ZScFtB3vTXSX4Iq32Wzyxn6n/fpLLzK6zpOJqy0WUc3a2YDz7w32bffx/GjoUHHtj0nm7dfCL/+c/h653K+XflDKparmKnWwq48agyfrdPenHVt025/HnM+wReWenvnrZ27Y9/rlkDX30FK1Zs+lnzfNkyWLp08+tWu3WDPfbwYwXusQc8uOBS2u84j8l/GhdYvKmqxoC0G42O6vNDLv4dg2ifjVnuiP5QfHwxxx/v53POn0v65BNfMFTzKP+onLXblkAbX7CzeO18zniqhD9eBLuuLaZrV5iz5gpab/0N//zKH2B17Ogfk9eWc+MnybcJCGV7gxKLBP744/DGG7B+/aZHZeWm5zOXzKG6qjU7bfnTH6atW+eT9YYNDS9/q61gu+18WfCOO/oz4927+8ftcy5li22XMenixzfrBhk3clrg21pf1VguNBhp3sJqn+ku12zTZ7Nv303zFtxWytpVde6m1WYdGw4rZfvpxSxfDivn7sf61Z35yyt1Vn5JKXTafN1njCzFDKq32vy1c0eXMqa0mFatfAHTxAUTsBYb6bfbMbRq5bt8al6r+f2cc2C33Zq8i5KKRQKfPduXkrdp44d3atNm06N9e2i99TdYqyr2+dlPf3h9iy18Yn7yswdo2fY7So+6iA4d/LSan9tu6xN3fWe7H0sk6brJOyyqhpRcFlb7zHS5qeZb12oBL73knxeNPA3n4IXfTGTlSli92l94cOgbC0hWJOq2Sj4d4Ps2C1i0yF+kUFUFa7/+Oa66BW8u3jRt48ZNz6uq4KijmmkCLyvzj1SKRl4JwJgkXzknj/Q3rDn7rIvCCC1w3Tt2Z/6q+UmnN3e52KXQ3ITVPjNdbjrvr91+tt661jwfJH9vj07+valemz590+9FI8/YbB3ZoErMHJPPVWOSO5p6C4mw2memyw1ruLZc/zzG4gi8OcnnqjGJv7DaZ6bLzeT96bw3Vz+PSuA5KF+rxiQ/hNU+M11uWMO15fLnUV0oIiIxlVECN7O+ZvapmVWY2RVBBdUY9Q2VlOkwSpksO8y4RHJVc2zbUW5zk7tQzKwlcDdwNLAI+MDMnnfOfRJUcA2pr6gAMrsAP5NlhxmXSK4Kq8gnl0W9zZkcgR8IVDjn5jrn1gNPAgOCCSs99V38n+nwYJksO8y4RHJVc2zbUW9zJicxfwIsrPX7IuCXdWcysxKgBKB792CvZW7Kxf+ZFgaks+ww4xLJVflchJbq5GXU25zJEXiyURo3K1xyzo1wzhU65wq7dOmSweo2V99QSZkOo5TJssOMSyRXxb1tN+Xa+Ki3OZMEvgjYudbvOwFLMgunccK8AD+TZce5MACiHydU4ikX2na2227U25xJF8oHQE8z2wVYDAwCTg8kqjSFeQF+JsuOc2GASFM1xyK0qLfZXN17pjbmzWb9gNuAlsDDzrl6/+0UFha6qVOnNnl9qdQ3OEKmAydksuyw3isStTA/c3EU9jab2TTnXGHd6RlVYjrnXgZezmQZklxzavwi0jSqxBQRiSklcBGRmNLNrESk0dTFlxt0BC4iElN5cQRe39GAjhREJGxR5Zm8SOBh0j8HEclVSuAhUoIXkTApgUdEyV1EMqWTmCIiMaUELiISU+pCEZFAqXswe3QELiISU0rgIiIxpQQuIhJTSuAiIjGlBC4iElNK4CIiMaUELiISU0rgIiIxpQQuIhJTGY1K3+iVma0A5jfx7dsBXwUYTlAUV+MorsbJ1bggd2PLx7h6OOe61J2Y1QSeCTOb6pwrjDqOuhRX4yiuxsnVuCB3Y2tOcakLRUQkppTARURiKk4JfETUAaSguBpHcTVOrsYFuRtbs4krNn3gIiLyY3E6AhcRkVpyKoGb2W/N7GMzqzazwjqvXWlmFWb2qZkdm+L9u5jZe2b2mZmNMbM2IcQ4xsxmJB7zzGxGivnmmdnMxHxTg44jyfqGm9niWrH1SzFf38Q+rDCzK7IQ1z/MbI6ZfWRmz5hZpxTzZWV/NbT9ZtY28TeuSLSlgrBiqbXOnc3sLTObnWj/FyeZp8jMVtX6+14TdlyJ9db7dzHvjsT++sjMemchpt1q7YcZZrbazC6pM0/W9peZPWxmy81sVq1p25jZ+EQuGm9mnVO8d0hins/MbEijV+6cy5kHsAewGzARKKw1fU/gQ6AtsAvwOdAyyfufAgYlnt8H/CHkeG8Brknx2jxguyzuu+HAZQ3M0zKx73YF2iT26Z4hx3UM0Crx/Cbgpqj2VzrbD1wA3Jd4PggYk4W/3Q5A78TzrYD/JImrCHgxW+0p3b8L0A94BTDgV8B7WY6vJfAl/jrpSPYXcBjQG5hVa9rNwBWJ51cka/fANsDcxM/OieedG7PunDoCd87Nds59muSlAcCTzrlK59wXQAVwYO0ZzMyAI4GxiUmjgJPCijWxvlOB0WGtIwQHAhXOubnOufXAk/h9Gxrn3OvOuarEr1OAncJcXwPS2f4B+LYDvi31SfytQ+OcW+qcm554vgaYDfwkzHUGaADwqPOmAJ3MbIcsrr8P8LlzrqkFghlzzr0NfF1ncu12lCoXHQuMd8597Zz7BhgP9G3MunMqgdfjJ8DCWr8vYvMGvi2wslaySDZPkA4FljnnPkvxugNeN7NpZlYSYhy1DUt8jQiS/y0AAANHSURBVH04xVe2dPZjmM7BH60lk439lc72/zBPoi2twretrEh02ewHvJfk5YPM7EMze8XM9spSSA39XaJuU4NIfRAVxf6q0c05txT8P2iga5J5Mt53WR/U2MzeALZP8lKpc+65VG9LMq3u5TPpzJOWNGMcTP1H3792zi0xs67AeDObk/hP3WT1xQXcC1yP3+br8d0759RdRJL3ZnwZUjr7y8xKgSqgPMViAt9fyUJNMi20dtRYZtYBGAdc4pxbXefl6fhugrWJ8xvPAj2zEFZDf5co91cb4ETgyiQvR7W/GiPjfZf1BO6cO6oJb1sE7Fzr952AJXXm+Qr/9a1V4sgp2TyBxGhmrYCBwP71LGNJ4udyM3sG//U9o4SU7r4zsweAF5O8lM5+DDyuxMmZE4A+LtH5l2QZge+vJNLZ/pp5FiX+zh3Z/Otx4MysNT55lzvnnq77eu2E7px72czuMbPtnHOh3vMjjb9LKG0qTccB051zy+q+ENX+qmWZme3gnFua6FJanmSeRfi++ho74c//pS0uXSjPA4MSVwjsgv9P+n7tGRKJ4S3glMSkIUCqI/pMHQXMcc4tSvaimW1pZlvVPMefyJuVbN6g1Ol3PDnF+j4Aepq/WqcN/uvn8yHH1Re4HDjRObcuxTzZ2l/pbP/z+LYDvi29meqfTlASfewPAbOdc7emmGf7mr54MzsQ/9n9b8hxpfN3eR44M3E1yq+AVTVdB1mQ8ltwFPurjtrtKFUueg04xsw6J7o8j0lMS182ztI24mzuyfj/SpXAMuC1Wq+V4q8g+BQ4rtb0l4EdE893xSf2CuBfQNuQ4hwJ/L7OtB2Bl2vF8WHi8TG+KyHsffcYMBP4KNF4dqgbV+L3fvirHD7PUlwV+H6+GYnHfXXjyub+Srb9wHX4fzAAWyTaTkWiLe2ahX10CP6r80e19lM/4Pc17QwYltg3H+JPBh+chbiS/l3qxGXA3Yn9OZNaV4+FHFt7fELuWGtaJPsL/09kKbAhkb+G4s+bTAA+S/zcJjFvIfBgrfeek2hrFcDZjV23KjFFRGIqLl0oIiJShxK4iEhMKYGLiMSUEriISEwpgYuIxJQSuIhITCmBi4jElBK4iEhM/T/hT2pgY8i1jAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nobounds=(None,None)\n",
    "out=minimize(logl,p0,args=(x,y),bounds=[(0,None),nobounds,(0,None),(0,None)])\n",
    "\n",
    "print(out)\n",
    "\n",
    "xfit=np.arange(-10,10,0.01)\n",
    "plt.errorbar(x,y,sig,fmt='go')\n",
    "plt.plot(xfit,model(xfit,*out.x),color='b')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some examples for timing different algorithms, with and without analytical derivatives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# default BFGS algorithm, numerical derivatives\n",
    "print('BFGS, numerical derivatives')\n",
    "%timeit -n100 out=minimize(logl,p0,args=(x,y))\n",
    "out=minimize(logl,p0,args=(x,y))\n",
    "print(out)\n",
    "\n",
    "# BFGS algorithm, supplying derivatives: jac=True means that function returns derivatives as well as value\n",
    "print('\\nBFGS, analytical derivatives')\n",
    "%timeit -n100 out=minimize(logl,p0,args=(x,y,True),jac=True)\n",
    "out=minimize(logl,p0,args=(x,y,True),jac=True)\n",
    "print(out)\n",
    "\n",
    "# Downhill simplex (Nelder-Mead) algorithm\n",
    "print('\\nDownhill simplex (Nelder-Mead)')\n",
    "%timeit -n100 out=minimize(logl,p0,args=(x,y),method='Nelder-Mead')\n",
    "out=minimize(logl,p0,args=(x,y),method='Nelder-Mead')\n",
    "print(out)\n",
    "\n",
    "# here's treating the problem with least squares\n",
    "print('\\nleast squares')\n",
    "from scipy.optimize import curve_fit\n",
    "%timeit -n100 fit=curve_fit(model,x,y,p0=p0,sigma=sig)\n",
    "fit=curve_fit(model,x,y,p0=p0,sigma=sig)\n",
    "print('least squares: ',fit[0])\n",
    "\n",
    "xfit=np.arange(-10,10,0.01)\n",
    "plt.errorbar(x,y,sig,fmt='go')\n",
    "plt.plot(xfit,model(xfit,*out.x),color='b')\n",
    "plt.plot(xfit,model(xfit,*fit[0]),color='r')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do the speeds of the different algorithms compare? Are they what you would expect, given what you know about the algorithms?\n",
    "<br><b>ANSWER HERE: </b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}