{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook, we investigate nonlinear least squares, and goodness of fit" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The text.latex.preview rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The mathtext.fallback_to_cm rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: Support for setting the 'mathtext.fallback_to_cm' rcParam is deprecated since 3.3 and will be removed two minor releases later; use 'mathtext.fallback : 'cm' instead.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The validate_bool_maybe_none function was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The savefig.jpeg_quality rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The keymap.all_axes rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The animation.avconv_path rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The animation.avconv_args rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n" ] } ], "source": [ "import numpy as np\n", "import matplotlib\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Nonlinear least squares \n", "\n", "OK, how about a nonlinear fit? Let's do a Gaussian function, with 3 parameters: amplitude, center, and sigma" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD4CAYAAAD4k815AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVLElEQVR4nO3df6wlZX3H8c+nILSrWLRcFVl2F1NCpcREekNWbYW4WFdKRZu2gax2I5vcmNRWmzYVsok1MSS1tqZttJqLUGi7gTb+qMRAZbUupolQL3T51QVZdMGVlb3W+KMlWdz47R8zFw7H82PumWfO/Hq/ks2eMzN35slz7/meZ77znWccEQIAdMfP1N0AAEBaBHYA6BgCOwB0DIEdADqGwA4AHXNiHQc97bTTYsuWLXUcGgBa6+677/5uRCxM266WwL5lyxatrKzUcWgAaC3bjxXZjlQMAHQMgR0AOobADgAdQ2AHgI4hsANAxxDYAaBjCgd229fbPmr7gYFlH7b9kO37bH/W9qmVtBIAUNh6Ruw3SNo+tGyvpPMi4lWSvi7p6kTtAgDMqHBgj4ivSPre0LLbI+J4/vZOSRsTtg2o1kUXZf+AjkmZY79S0m3jVtpesr1ie2V1dTXhYYEJCN7ooSSB3fZuSccl7Rm3TUQsR8RiRCwuLEyd6gAAMKPSc8XY3inpUknbgufsAUDtSgV229slvU/ShRHxVJomAQDKWE+5402SvirpHNuHbe+S9FFJp0jaa3u/7U9U1E4AQEGFR+wRccWIxdclbAsAIAHuPAWAjiGwA1Wh1BI1IbADQMcQ2AGgYwjs6Kc9e6Q775TuuEPasiV7D3QEgR3dNS5479kjLS1Jx45l7x97LHtPcEdHENjRTZOC9+7d0lND99M99VS2HOgAAju6aVLwfvzx0T8zbjnQMgR2dNOk4L1p0+h145YDLUNgRzdNCt7XXCNt2PDc5Rs2ZMuBDiCwo5smBe8dO6TlZenkk7Plmzdn73fsmH87gQqUnrYXaKS1IL1rV3YBdfPmZ4P62vprr81e79tXSxOBqhDY0V11Bu+1Ustjx7JSy8EvFaBipGKA1KiTR80I7EBq1MmjZgR2IDXq5FEzAjuQGnXyqBmBHUiNOnnUbD3PPL3e9lHbDwwse7HtvbYfyf9/UTXNBFqEOnnUbD0j9hskbR9adpWkL0XE2ZK+lL8HsGOHtHWrdOGF0qFDBHXMVeHAHhFfkfS9ocWXSboxf32jpLemaRbQAjz6Dg1V9gall0bEEUmKiCO2XzJuQ9tLkpYkaRMXkdAE025aWgva3JmKlpnbnacRsSxpWZIWFxdjXsdFzxGU0UNlq2KetH26JOX/Hy3fJACkeVBG2cB+i6Sd+eudkj5Xcn8AgJLWU+54k6SvSjrH9mHbuyT9uaQ32n5E0hvz9wCAGhXOsUfEFWNWbUvUFgBAAtx5CgAdQ2AHgI7hQRtAVSi1RE0YsaPZKPsD1o3ADsxi7dF3d9yRPfqOpyOhQQjswCiTAjePvkPDEdiBYdMCN4++Q8MR2IFh0wJ31Y++I82DkgjswLBpgbvKR9+R5kECBHZg2LTAXeWj70jzIAECOzBsWuCu8tF3Vad50AvcoAQMWwvQu3ZlKZHNm7OgPhi4d+yQrr02e53yRqRNm7L0y6jlg3gICCZgxA6MUtczS6tM86A3COxAk1SZ5kFvkIpBc62V/R07lpX9DadDuqqqNA96gxE7mqkPZX/Mg4OKENjRTJT9lcOXRq8R2NFMlP0BMyOwo5mqvLsT6Lgkgd32H9l+0PYDtm+y/bMp9oseK1r2R8oB+CmlA7vtMyT9oaTFiDhP0gmSLi+7X/QcZX/AzFKVO54o6eds/1jSBklPJNov+qzpZX9NbBOgBCP2iPi2pL+U9LikI5J+EBG3D29ne8n2iu2V1dXVsocFqrdvXzODN9P6YooUqZgXSbpM0lmSXi7p+bbfPrxdRCxHxGJELC4sLJQ9LNBPfajvR2kpLp5eLOmbEbEaET+W9BlJr02wX6C/xp0tUN+PAlIE9sclbbW9wbYlbZN0IMF+ge6aNZ1CfT8KSJFjv0vSpyTdI+n+fJ/LZfcLdFaZdAr1/SggSR17RPxZRPxSRJwXEe+IiGMp9gt0Upl0SpH6fi6u9h53ngLzViadMq2+n4urEIEdmL+y6ZRJDwHh4ipEYEebtTXlUOVTkri4CvGgDTTduBuExqUcpOZPO1DkmaqzKvrMVHQaI3a0U9tTDlU9U5VnpkIEdrQVKYfRmDwNIhWDtiLlMF7TJ09D5Rixo51IOQBjEdhRraoehEHKARiLVAzai5QDMBKBHWgjvsgwAakYAOgYAjsAdAyBHQA6hhw7UBfy5KgII3YA6BhG7EAXcTbQa4zYAaBjkgR226fa/pTth2wfsP2aFPtFyxWZL73snan79jE6BYakSsX8jaR/i4jftn2SpA3TfgAd1+b50oGWKz1it/1CSa+XdJ0kRcTTEfH9svvFHFUxn0vb50sHWixFKuYVklYl/b3t/7L9SdvPH97I9pLtFdsrq6urCQ6LRmO+dKA2KQL7iZLOl/TxiHi1pP+TdNXwRhGxHBGLEbG4sLCQ4LBojFEj/rIPbAYwsxSB/bCkwxFxV/7+U8oCPfqM+dKB2pQO7BHxHUnfsn1OvmibpP8uu1+0HPOlA7VJVRXzB5L25BUx35D0zkT7RZsxXzpQiySBPSL2S1pMsS8AQDnceQoAHUNgB4COIbCjPkWmHEA1qnrIOBqBwN53dQXXcVMOENyB0gjsfVZncGXKAaAyBPY+qzO4MuUAUBkCe5+lCK6zpnKYcgCoDIG9z8oG1zKpHKYcACpDYO+zssG1SCpn3IMwmHIAqAzPPO2ztSC6a1c26t68OQvqRYNr2VQOUw4AlSCw912Z4LppU5Z+GbUcQG1IxWB25MmBRiKwY3bkyYFGIhWDcsiTt89aieqxY1mJ6nquq6AVGLEDfcJUDr1AYMd0TBjVHUzl0AsEdqBPmMqhFwjsQJ8wlUMvENiBPqFEtReSVcXYPkHSiqRvR8SlqfaLjqOSZr7K3m2MVkhZ7vgeSQckvTDhPjEPBNd+oUS185KkYmxvlPQbkj6ZYn8AgNmlyrH/taQ/lfSTcRvYXrK9YntldXU10WEBAMNKB3bbl0o6GhF3T9ouIpYjYjEiFhcWFsoeFk0ybmpetBf3LrRaihH76yS9xfYhSTdLeoPtf0qwXwDADEoH9oi4OiI2RsQWSZdL+veIeHvplgEAZkIdOyab9Zmm6C7SNI2XdHbHiNgnaV/KfaJG4yaMkqh7BhqMETvGY8IooJWYjx3jMWFUd1HF1GmM2DEeE0YBrURgx3hMGAW0EoEd4/FMU6CVyLEXtVbe1cbcZJm2M2FU//BM1NZjxJ4Ktb3oAp6J2gkEdgDPosS1EwjsbcEZAeaBEtdOILADeNa0ElemmGgFAjuAZ00qcSX/3hoEdgDPmlTiSv69NSh37II2l2KiecaVuJJ/bw1G7ACKYYqJ1mhXYKcypB48+g4SU0y0SLsCe12oBACYYqJFyLFP0/aHTXB7OFJq+xQTPbkexYh9miKVAE0d0VOeBvRS6cBu+0zbX7Z9wPaDtt+TomGNMa0SoMnBk/I0tBHX0kpLMWI/LumPI+KVkrZK+n3b5ybYbzNMqwSYR/Cc9YyA8jSgl0oH9og4EhH35K9/JOmApDPK7rcxplUCVB08y5wRUJ4G9FLSHLvtLZJeLemulPut1bRKgKqDZ5kzAsrTUAdSKbVLVhVj+wWSPi3pvRHxwxHrlyQtSdKmto0YJ1UCXHNNNoIeDL4pg2eZM4K1L59du7IR/+bNVMWgmI5XjXRdkhG77ecpC+p7IuIzo7aJiOWIWIyIxYWFhRSHbYaqa3vLzra3Y4e0dat04YXSoUMEdaBKDTlbSVEVY0nXSToQER8p36QWqjJ4MtsegHVKMWJ/naR3SHqD7f35v0sS7Pe5mlorXjVm20PTMMVE45XOsUfEf0hygraM1/a7P8titj2guHF3l/boLux23HnKyHQ0yhmBYnqWtmxHYGdkOhrljEAxPRsctiOwMzIdjdn20DVVXUvr2eCwHYGdkel4lDOiK6pMl/RscNiOwM7IFOi+KtMlPRsctmc+9rrngaa8C5iuTOVJlemSnt2F3Z7APg9lJuFvcuBvctvQHWXLkjdtyn5m1PIUqh4cNqicsh2pGADNVzaV0uZ0ScPKKQnsANIom0opei2tIfOxPEfDyilJxbQF6RQ0XdFUyqSUZ93X0qYZ1/aGlVMyYgeQRhNSKXXNKdWwckoCO4A06i5LrjPP3YQvtQGkYrqgiaes6Kc6UymT8txVf7k0rJySwA6gG+rOcxf5UitTUr0OpGIAdEPD8tx1aldgb/ME/00s0QK6JEWeuyOfU1IxALqhSJ67rQPDdSKwr2nQ7cAAZtT0Ovg5aVcqpioNux0YAMpIEthtb7f9sO2Dtq9Ksc9KjMufNeF24I7k9oDSJl1L6+tD7depdGC3fYKkj0l6s6RzJV1h+9yy+52rusukgL4oE5jrPrNu0ZdKihH7BZIORsQ3IuJpSTdLuizBfueHMikgnXEj7rKBuc4z67q/VNYpRWA/Q9K3Bt4fzpc9h+0l2yu2V1ZXVxMcNqGG3Q4MdFLZwFznmXUT0rXrkCKwe8Sy+KkFEcsRsRgRiwsLCwkOm1Ddc1wAfVA2MNd5Zl207Q251yZFYD8s6cyB9xslPZFgv/PFQ6GBapUNzHWeWaf4Upljjj5FYP+apLNtn2X7JEmXS7olwX4BdEnZwFznmXXZts85R186sEfEcUnvlvQFSQck/UtEPFh2vzNpaslgi66mA5VJEZirPLOe9Dkt2/Y55+iT3HkaEbdKujXFvjqn7AN+gS6Zx52hs+y3yOe0TNvnfOGXO0+r1rKr6UAvVf05nfOFXwJ71bj5CWi+qj+nc77w25/A3uRnIZKDB+pV9Yh6zhd++xHYm/wsxJbd0QZ00jxG1HMsqe5HYK8zzz3tm5ocPFC/jt2k2I/52Ndz11gVJl1NJwcPNEOH5nLvx4i9yZN8NbltAFqpH4G9yZN8NbltQBM1ZD6WJutHKqbIsxDr0uS2AXiulnyhdCewT3tmaZPzZ01uG5Aaf+OV60YqhpJBAHhGNwI7JYMA8IxupGIoGQTQBnNKQ3VjxE7JIAA8oxsj9muuyXLqg+mYppUMcsEIaL6OfE67MWLv2O3AAFBGN0bsEiWDAJDrxogdAPCMUiN22x+W9JuSnpb0qKR3RsT3E7SrfzjLAJBI2RH7XknnRcSrJH1d0tXlm1Qh5pgA0AOlAntE3B4Rx/O3d0raWL5JAIAyUubYr5R0W8L9AQBmMDXHbvuLkl42YtXuiPhcvs1uSccljZ2cxfaSpCVJ2sSNQwBQmamBPSIunrTe9k5Jl0raFhExYT/LkpYlaXFxcex2AIByylbFbJf0PkkXRsRT07YHAFSvbI79o5JOkbTX9n7bn0jQJgBACaVG7BHxi6kakgSljADAnacA0DUEdgDoGAI7AHQMgR0AOobADgAdQ2AHgI4hsANAxxDYAaBjCOwA0DGeMG9XdQe1VyU9NuOPnybpuwmbkxJtmw1tmw1tm02b27Y5Iham7aSWwF6G7ZWIWKy7HaPQttnQttnQttn0oW2kYgCgYwjsANAxbQzsy3U3YALaNhvaNhvaNpvOt611OXYAwGRtHLEDACYgsANAxzQysNv+HdsP2v6J7cWhdVfbPmj7YdtvGvPzL7a91/Yj+f8vqqid/5w/EnC/7UO294/Z7pDt+/PtVqpoy4hjfsD2twfad8mY7bbnfXnQ9lVzatuHbT9k+z7bn7V96pjt5tZv0/rBmb/N199n+/wq2zNw3DNtf9n2gfwz8Z4R21xk+wcDv+v3z6Nt+bEn/o5q7LdzBvpjv+0f2n7v0DZz6zfb19s+avuBgWWF4tRMn9GIaNw/Sa+UdI6kfZIWB5afK+leSSdLOkvSo5JOGPHzfyHpqvz1VZI+NIc2/5Wk949Zd0jSaXPuww9I+pMp25yQ9+ErJJ2U9+25c2jbr0s6MX/9oXG/n3n1W5F+kHSJpNskWdJWSXfN6fd4uqTz89enSPr6iLZdJOnz8/z7Kvo7qqvfRvx+v6Ps5p5a+k3S6yWdL+mBgWVT49Ssn9FGjtgj4kBEPDxi1WWSbo6IYxHxTUkHJV0wZrsb89c3SnprJQ3N2bak35V0U5XHqcAFkg5GxDci4mlJNyvru0pFxO0RcTx/e6ekjVUfc4oi/XCZpH+IzJ2STrV9etUNi4gjEXFP/vpHkg5IOqPq4yZUS78N2Sbp0YiY9W730iLiK5K+N7S4SJya6TPayMA+wRmSvjXw/rBG/5G/NCKOSNkHQ9JLKm7Xr0l6MiIeGbM+JN1u+27bSxW3ZdC789Pf68ec5hXtzypdqWxEN8q8+q1IP9TeV7a3SHq1pLtGrH6N7Xtt32b7l+fYrGm/o9r7TdLlGj/oqqvfpGJxaqb+OzFJ82Zg+4uSXjZi1e6I+Ny4HxuxrNJ6zYLtvEKTR+uvi4gnbL9E0l7bD+Xf4JW1TdLHJX1QWf98UFmq6MrhXYz42ST9WaTfbO+WdFzSnjG7qaTfRjV3xLLhfpj7395zDm6/QNKnJb03In44tPoeZWmG/82vpfyrpLPn1LRpv6O6++0kSW+RdPWI1XX2W1Ez9V9tgT0iLp7hxw5LOnPg/UZJT4zY7knbp0fEkfy07+gsbZSmt9P2iZJ+S9KvTNjHE/n/R21/VtnpVekAVbQPbV8r6fMjVhXtz3Ur0G87JV0qaVvkycQR+6ik30Yo0g+V9dU0tp+nLKjviYjPDK8fDPQRcavtv7N9WkRUPtFVgd9Rbf2We7OkeyLiyeEVdfZbrkicmqn/2paKuUXS5bZPtn2Wsm/X/xyz3c789U5J484AUrhY0kMRcXjUStvPt33K2mtlFw4fGLVtSkN5zLeNOebXJJ1t+6x8ZHO5sr6rum3bJb1P0lsi4qkx28yz34r0wy2Sfi+v8tgq6Qdrp9FVyq/fXCfpQER8ZMw2L8u3k+0LlH2u/2cObSvyO6ql3waMPZuuq98GFIlTs31G53FFeIYryG9T9k11TNKTkr4wsG63sqvED0t688DyTyqvoJH0C5K+JOmR/P8XV9jWGyS9a2jZyyXdmr9+hbIr2fdKelBZKmIeffiPku6XdF/+h3D6cNvy95coq7R4dI5tO6gsb7g///eJuvttVD9Ietfa71bZKfHH8vX3a6Baq+J2/aqyU+/7BvrrkqG2vTvvo3uVXYx+7ZzaNvJ31IR+y4+9QVmg/vmBZbX0m7IvlyOSfpzHtl3j4lSKzyhTCgBAx7QtFQMAmILADgAdQ2AHgI4hsANAxxDYAaBjCOwA0DEEdgDomP8HPqCiAr/XMBMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "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 func(x,amp,cent,wid) :\n", " return amp*np.exp(-0.5*(x-cent)**2/wid**2)\n", "\n", "# simulate some data. Note that I've defined a par array, but pass it as\n", "# individual arguments using *par\n", "x=np.arange(-10,10,0.5)\n", "par=[10,0,3]\n", "sig=1\n", "y=func(x,*par)+np.random.normal(0.,sig,size=len(x))\n", "plt.errorbar(x,y,sig,fmt='ro')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use [scipy.optimize.curve_fit()](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.htm) to do the nonlinear fit. curve_fit returns both the best fit paramters and the covariance matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Moments for initial parameters" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "11.03839521600013 19 -0.5 -0.20767008783541083\n" ] } ], "source": [ "amp=y.max()\n", "cent=x[y.argmax()]\n", "mean=(y*x).sum()/y.sum()\n", "sigma=(y*x**2).sum()/y.sum()\n", "print(amp,y.argmax(),cent, mean)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "parameters: [10.12727042 -0.02159847 2.85080193]\n", "covariance matrix: \n", " [[ 1.87759941e-01 -1.08577429e-06 -3.52419268e-02]\n", " [-1.08577429e-06 1.98371456e-02 6.70284425e-07]\n", " [-3.52419268e-02 6.70284425e-07 1.98416511e-02]]\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAArP0lEQVR4nO3dd3wUdfrA8c+TnlASQkJLIXRBpCsIKpw029nP44h3nKKcP+88PVFPxXKNExvenQUPy9mwnb0hIAqcCkoREKT30BJqIL18f3/sIiHsJpvs7M7O7vN+vfLK7szszJNJ8ux3v/PM9yvGGJRSSjlXlN0BKKWU8o8mcqWUcjhN5Eop5XCayJVSyuE0kSullMPF2HHQtLQ0k5OTY8ehlVLKsZYuXbrPGJNee7ktiTwnJ4clS5bYcWillHIsEdnmabl2rSillMNpIldKKYfTRK6UUg6niVwppRxOE7lSSjmcJnKllHI4nxO5iDwvIvkisqrGsodFZK2IrBSRd0UkJSBRKqWU8qohLfIXgPNqLZsD9DTG9ALWA3dZFJdSSikf+ZzIjTELgAO1ls02xlS6ny4CMi2MTSnrDRvm+lIqjFjZR34tMNPbShGZICJLRGRJQUGBhYdVqgZN1CoCWZLIRWQSUAnM8LaNMWa6MWaAMWZAevpJQwUopZRqJL/HWhGRccBFwHCj88YppVTQ+ZXIReQ84I/AUGNMsTUhKaWUaoiGlB++BiwEuolInoiMB54AmgFzRGS5iDwdoDiVUkp54XOL3BjzCw+Ln7MwFqWUUo2gd3YqpZTDaSJXyipa+qhsoolcKaUcThO5Uko5nCZyFTlmzIBFi2D+fMjJcT1XKgxoIlfho65EPWMGTJgAZWWu59u2uZ5rMldhQBO5Cg/1JepJk6C41j1rxcWu5Uo5nCZyFR7qS9Tbt3t+nbflSjmIJnIVHupL1NnZntd7W66Ug2giV+GhvkQ9eTIkJZ24LinJtVwph9NErsJDfYk6NxemT4f4eNfz9u1dz3NzgxunUgHg9zC2SoWEYwl5/HjXBc/27V1JvGaizs2FZ55xPZ43L+ghKhUomshV+LAzUR8rfSwrc5U+1n4TUSqAtGtFKX9pjbqymSZypfylNerKZprIlfKX1qgrm2kiV8pfWqOubKaJXCl/aY26sllD5ux8XkTyRWRVjWWpIjJHRDa4v7cITJhKBdbBonK27CtiY/4R8o+UUl1tfH+x1qgrmzWk/PAFXJMtv1Rj2Z3AXGPMFBG50/38j9aFp1Rg5B8pZdbqvcxds5fv8w6zv6j8hPVN4qI5tV0yQ7ulc8FpbemQ1qTuHWqNurJRQyZfXiAiObUWXwIMcz9+EZiHJnIVwjYmpPJUxkA+nPI5FVWGDmlNGN69FV1bNyO1SRzRUcLhkgo2FxSxbPtBHp61jodnrWNo13RueO9xzizcoYlahRx/bwhqbYzZDWCM2S0irbxtKCITgAkA2XoRSAXZ4eIKHvtsPS/3u474mChyB2QxdmA2XVo1RURcGx2bb7NGot51qIR3luXxwtdb+UWPMYw8sIH7DhSTlZp00jGUskvQ7uw0xkwHpgMMGDCgAR2QSjWAh9bywk37ufXN5ewtLOUXZ2QzcVQ3UpvE+bS7dimJ/O7cLlx3dkeeH3cXT2Scyeh/LODPF5/Klf0zj78JKGUjf6tW9opIWwD393z/Q1LKGsYY/vnZBsY+u4iE2GjevXEIky87zeckXlNCbDQ37vqWz1Y8T6/MZG5/ayW3vLGc0ooqa4IdNuz4JwKlGsjfRP4BMM79eBzwvp/7U8oSZZVV3PrmCh77bD2X9sngo5vOondWit/7bVd+hBnXDWLiyK68v3wXY59ZxP6jZf4HrJQfGlJ++BqwEOgmInkiMh6YAowUkQ3ASPdzpWxVXF7JNf9ZzLvf7eS2UV2ZelVvmsRb14sYHSXcNLwL03L7sXpXIZdP+5qdh0os279SDdWQqpVfeFk13KJYlPJbcXkl176wmG+3HGDqVb25vF9mwI51/mltaZ2cwLjnv2XM9IW8FteczPLCgB1PKW/0zk4VNkorqn5M4o/9vE9Ak/gx/bJb8Mr4gRwqrmBMj5+zN7aeenOlAkATuQoL1dWGif9dwaLNB5h6VR8u6ZMRtGP3zkrhlfEDORCbxK9PuZIjpRVBO7ZSoIlchYkpn67l45W7ufuCU7i0b/CS+DG9s1J46rohrG/WmhtnLKOiqjroMajIpYlchZ4GluK9+s12pi/YzK/ObM/1Z3cMWFj1GdatFQ9cfhr/27CPP3+42rY4VOTRRK4cbfmOQ/zpg9UM7ZrO/T89NXA36Bybym3+fNdUbl5m/7lqQBa/Oacjryzazn+X7AhMLErVoolcOdb+o2Xc+MpSWjWP559j+hAd5UcSrytRN3Aqt9tHd2NI55ZMem8VK/MONT4mpXykiVw5UlW14ebXl7OvqJynr+5PSlLD79b8UX2JuoFTucVER/GvMX1JbxrP/72yjEK9+KkCTBO5cqRn/7eZLzfu4y8Xn0rPjGT/dlZfom7EVG4tm8bz+Ni+7Cks5d73VnndDvC520YpbzSRK8dZveswj8xex/k92/Dz07P832F9ibqRU7n1y27BzcO78P7yXbz33U7PGzWw20YpTzSRK0cprajilteX0yIpjr9fdpo1FzfrS9R+TOX225905vScFtzz3ip2HCg+eYMGdtso5YkmcuUoU2auZUP+UR75WW9aNGIUQ4/qS9R+TOUWHSVMvaoPAtz23xUnTyHXiG4bpWrTRK4cY/HWA7zw9VZ+PTiHc7qmW7djXxJ1bi4MGgRDh8LWrQ2ajzMrNYl7LurON1sO8Oq3tRK0L902OsStqocmcuUIpRVV3Pn2SjJbJHLHed2sP4AfidoXVw3IYkjnlkyZuZZdNUdK9KPbRqljNJErR3jqi41sKihi8mWnkRQXtImtLCMiTLm8F1XVhrvf/R5j3F0sfnTbKHWMJnIVWjyU4q3dU8hT8zZxWd8MhlrZpRJkWalJ3D66G/PWFfDe8hpVLAH+NKDCnyZyFTo8lOKZCRP4+M5HaZ4Yy70X9bA3Pn+4+7nHDc6hb3YKkz9ew+ESvVFIWUMTuQodHkrxpLiYn787jXsu7N6ouTZDTXSU8NdLenKgqJzH5qy3Zqd6MTTiaSJXocNLyV27wn1cZsPQtIHSMyOZXw5qz0sLt7Jq52G7w1FhQBO5Ch1eSvEqMzIDN6qhTW4d1Y0WSXHc9/6qk2vLlWogSxK5iPxBRFaLyCoReU1EEqzYr4owHkrxyuMTiHvwgRO3C4OuhOTEWO66oDvLth/irWV5doejHM7vRC4iGcDvgQHGmJ5ANDDG3/2qCOQuxTPx8VQDu5JbUTnt32FbxXF53wz6t2/BgzPXUhjt/P5/ZR+rulZigEQRiQGSgF0W7VdFmtxcPrjwGjr+8SPmz15M0jW/sjui4+bNc31ZJCpK+PPFp3KguJyn2g2ybL8q8vidyI0xO4FHgO3AbuCwMWZ27e1EZIKILBGRJQUFBf4eVoWp0ooqpmSfw2lH93DVAAtGNmwIixO1L3pmJHN530yeb9ufHfEehuPVIW6VD6zoWmkBXAJ0ANoBTUTk6trbGWOmG2MGGGMGpKc796YOFVjPfbmF3fHNuXfb5/7N+OMgt4/uRpQxPJh19okrdIhb5SMrulZGAFuMMQXGmArgHWCwBftVEWbf0TKmzdvEyAMbOOOIl/G7w1Cb5AQmjOrBR2ndWbrt4PEVOsSt8pEViXw7MEhEksRVIzYcWGPBflWE+dfcDZRUVHHn9vl2h2ItH7pHfnNOR1o1i+evH/1wfBwWHeJW+ciKPvJvgLeAZcD37n1O93e/KrJsLjjKq99s5xdnZNGp9GD9L3AKH7tHmsTHcNvobizfcYgPV+52LWzkzEQq8lhStWKMud8Yc4oxpqcx5pfGmDIr9qsix4OfriU+Joqbh3e1OxRrNaB75Ip+mfRo25wHZ66lrLLKtyFu9WKoQu/sVCFg8dYDzFq9lxuGdiK9Wbzd4VirAd0j0VHCXRecws5DJcxYtL3+IW71Yqhy00SubGWMYfLHa2jdPJ7rzu5odzjWa2D3yNld0hnSuSVPfLGRo2WVdQ9xqxdDlZsmcmWrmav2sHzHISaO6kZiXHT9L3BaV0IjZgC6Y/QpHCgq59n/ba5733oxVLlpIle2qao2PDp7HV1aNeWKfpnHV3i7MceJXQmNmAGod1YK5/dswzMLNrP/aB2Xm/RiqHLTRK5s8+53O9lUUMStI7v6dvOPU7sSGjED0MRR3SipqOLJLzZ530jn+1RumsiVLcorq/nHZ+vpmdGc83q28e1FEdSV0LlVU37WP4tXFm0jL6655410vk/lpolc2eKNJTvIO1jCxFHdfB9rPMK6Em4Z2QUEHssa4n0jne9ToYlc2aC0ooonPt/AgPYtGNaQyZQjrCuhbXIivx6cwztpp7IuMc3ucFQI00SuAqOOyR9eXriNvYVl3Da6Aa1xiMiuhP8b2ommVeU8knWW3aGoEKaJXAXV0bJKps3fxNld0hjUsWXDdxBhXQktmsQxYfe3zEntwnfbw2joAmUpTeQqqJ7/cgsHisqZOKqb3aE4xrXP/42WTeJ4eNY6u0NRIUoTuQqaQ8XlPLNgMyN7tKZPVord4ThGk/gYfvuTzny9aT9fbdxndzgqBGkiV0Hz7wWbOVpeycRRYTYwVhDkDsqmXXICD81ad3yYW6XcNJGroMg/UsoLX23lp73acUobL3XRyqv4mGhuGdGVFTsOMfuHvXaHo0KMJnIVFE99sYnyqmr+MDJCW+MWzAd6eb8MOqY14dHZ66iq1la5Ok4TuQq4nYdKePWb7VzZL5MOaU3sDsexYqKjuHVUV9bvPcoHKyJnKjxVP03kKuAen7sBgN+P6GJzJM53Qc+2nNquOY/N2UB5ZbVroQWtfeVsmshVQG3ZV8R/l+YxdmA2GSmJdofjeFFRwm2ju7H9QDFvLNlhdzgqRFiSyEUkRUTeEpG1IrJGRM60Yr/KoWqMGZ7cvQuXrpnHjT/pdHx9HXd9+iTCW6DDuqZzek4LHp+7gZLyKrvDUSHAqhb5P4FPjTGnAL2BNRbtVzlNrTHDU/ftZsrMx2n1wds2BxY+RITbR59C/pEyXlq41e5wVAjwO5GLSHPgHOA5AGNMuTHmkL/7VUHgb8vYEw9jhseWlYb+mOEOc0aHVIZ1S2fa/E0UllbYHY6ymRUt8o5AAfAfEflORJ4VkZNKE0RkgogsEZElBQUFFhxWhaQIGjPcbreN6sah4gqeXVDPlHAq7FmRyGOAfsA0Y0xfoAi4s/ZGxpjpxpgBxpgB6ekNGLpUhSZvrfkIGzPcTj0zkrnwtLY8++UW9tU1JZwKe1Yk8jwgzxjzjfv5W7gSu4pEkydTlVirOiWMxwy3262julJaUcVTdU0Jp8Ke34ncGLMH2CEix4azGw784O9+lTOZsWN5Yswd7G6WhoGIGDPcTp3Sm3Jl/0xeWbSNnYdK7A5H2cSqqpWbgBkishLoA/zdov0qh1mwYR+PtTqdOedfjUTImOF2u3mEa9iDYzdeqchjSSI3xix393/3MsZcaozREfAjkDGGR2evIyMlkTH5K+wOJ2JkpCSSOyib/y7NY3PBUbvDUTbQOzuVZWb/sJeVeYe5eUQX4ky13eFElBuHdSY+Joqpc9bbHYqygSZyZYmqasPU2evpmN6Ey/tm2B1OxElvFs+1Qzrw0crdrN512O5wVJBpIleW+GjlLtbtPcIfRnQlJrqOP6sat++Tk+N6rixx/TkdSU6M5dHZtVrlgbjxS4UUTeSRysKEWllVzT8+28ApbZpx4Wlt6z5mjdv32bbN9VyTuSWSE2O5YWgnPl+bz5KtB+wORwWRJvJIZHFCfXtZHlv2FTFxVDeiosT7hh5u36e4WG/ft9C4we1JbxbPQ5/qlHCRRBN5JLIwoZZVVvGvuRvpnZXCiO6t6t5Yb98PuKS4GG46tzPfbj3Agg06UXOk0EQeifxNqDW6ZSoy29P/q0+4fVQ3ROpojYPevh8kY07PJrNFIg/PWqut8gihiTwS+ZNQa3XLNN27k4dnP8GQb2fV/9rJk12369ekt+9bLi4mij+M6MqqnYXMXLXH7nBUEGgij0T+JFQP3TLx5WVI7W4ZT5M/5Oa6btePj3c919v3A+bSvhl0adWUR2evo5J6Pikpx9NEHon8Saj+dsvk5sKgQaC37wdUdJQwcVQ3NhUU8W56D7vDUQGmiTxSNTahaj+3Y4w+tTW9MpP5R+YQyiTa7nBUAGkiVw0zeTLVOkytI7imhOvGzvhkXmvVy+5wVABpIlcNk5vLS9dMYmfzdB2m1gHO6pzGoMPbeSLjTIrLK+0ORwVIjN0BKGdZmXeIPzXrw4GfXs+teV9F9Gz2TiCvvsqLL95GbEkJR2e0hakP6ZtuGNJErnxmjGHKzLWkNonj+t2L7Q5H1cddKhpf4qoyap6/C3P9BFcNiybzsKJdK8q7WoMtLdiwj6837eemczvTrKrctrCUjzyUikqJDokQjjSRK59UVlXz94/XkJWayNiBWqHiCF5KQo0OiRB2NJErn7y+eAfr9h5h0gXdiY/RUjZH8FISeiitTZADUYGmiVzVq7C0gqlz1jOwQyqjT9Uk4Bge7uAtj0/gz2fmsn7vEZuCUoFgWSIXkWgR+U5EPrJqnyo0PPH5Rg4Wl3PvRT3qHxjLF55u31fW83AHb/lTT/N53xH89aMfdECtMGJl1crNwBqguYX7VIHkQzLdGp/Cf77awpX9MumZkRz4mJS1cnPhmWdcj+fNoynwh6+28OcPf2DumnxG9Ghta3jKGpa0yEUkE7gQeNaK/anQ8UD7ocRGR3H76G52h6IscvWg9nRu1ZS/ffwDZZVVdoejLGBV18o/gDsAr1Oni8gEEVkiIksKCgosOqwKpK+bZzErtSs3DutEq+YJdoejLBIbHcW9F/Vg6/5iXvhqq93hKAv4nchF5CIg3xiztK7tjDHTjTEDjDED0tPT/T2sCrDKqmr+0v5cMsoOc93ZHU/eQPu5HW1o13TOPaUVj3+4goKRF9gdjvKTFS3yIcDFIrIVeB04V0ResWC/ykYvL9rG2iatuGfbPBJitdwwHN1zYXdKo2J4JOtsu0NRfvI7kRtj7jLGZBpjcoAxwOfGmKv9jkzZJr+wlKmz13POoS2cd2C93eGoAOmY3pRr9izjzfTTWLb9oN3hKD9oHbk6yQMz1zJy+Vyee/EOZP58yMlxjduhws7NeV/Ruvwok95dRWWVh0tctYZpUKHJ0kRujJlnjLnIyn2q4Fq0eT9Vr7zCg7MeJ7a0xLVw2zbXPJ2azMNO0+oK7t/2OWt2F/LC11vtDkc1krbI1Y8qqqq57/1V3P3ly8SWlZ64slgHWwpX5x1Yz0+6pfPYnPXsPlxidziqETSRqx/956strN97lNaHvZSH6mBLzlRPhZEAf764J5XVhr98+EOwolIW0kSuANi2v4ipc9YzskdrROfljDjZLZP4/fAuzFy1hy/W5tsdjmogTeQKYwx3vv09sVFR/PWSnh4HW9J5OcPf9Wd3pFN6E+59fxVFZTotnJNoIle8sXgHCzfv5+4Lu9MmOcHjYEs6L2f4i4uJYsoVvdh5qISHPl1rdziqATSR+8qJZVg+xLy3sJTJn6zhzI4tGXN61vEVubkwaBAMHQpbt2oSD0czZsCiRVCjxPT0nFTGnZnDiwu3sWjzfrsjVD7SRG4FJyZ5XF0q97y3ivLKah64/DRrhqhVzuCez5OyMtfzGiWmd5zXjezUJP749kpKonRaXyfQRB7BPlixizk/7OXWkV3JSWtidzgqmDzM53msxDQpLoYHr+jFtv3FPKy37zuCJvJQF6DW/q5DJdzz3ir6t2/B+LM6WL5/FeK8lZK6l5/ZqSW/HNSe/7Tpz+JmGUEMTDWGJvIIVF1tmPjmCqqrDY9d1YeYaP0ziDg+lJhOOvwdC6f9mv4fvUp1dnu9szeE6X9wBHr+qy0s3Lyf+37ag+yWSfW/QIWf+kpMZ8wg4cYbaHNkP1FA1I7tOkxDCNNEHmHW7inkoU/XMbJHa64akFX/C1R4qq/EtI4+dBV69JK00x3rP/dhkoeS8ipufm05zRNjtEpFnTSf5wm89KGb7dvRv5rQoy3yCHLf+6tYn3+ER6/qQ1rTeLvDUaHMSx96QUorKjwNd6ts5axE7tB67VDw3yU7+O/SPH73k84M7erjVHs6nVvk8tCHXpmQyN+GXM2UmXrXZ6hxViK3i4c74Jxk3Z4j3Pv+Ks7s2JJbRnS1OxzlBB760GOefYYW1/2a577cwkcrd9kbnzqB9pHXx9sdcBDat62733xMWRnJ3btwxYhruPmO+4mO0h5O5SMPfeiTKqv5fudh7nhrJae0aUbnVs3si68uDbh2FA60RV6f+q7eh2JrvcabjwBtDu3lLx/9k1YfvG13ZMrh4mKieDK3H4mx0fzm5aUc1VESQ4LfiVxEskTkCxFZIyKrReRmKwILGXXdAVfHeBW28vDmE11aoqVjyhJtkxN5fGxftuwr4pbXv6Oq2jRuR3rNyzJWtMgrgYnGmO7AIOC3ItLDgv2GhrrugAt0rW1jW/v13H6tlL8Gd0rj/p+eymdr8nngkzV2hxPx/E7kxpjdxphl7sdHgDVA+AzOUNcdcIFMmH609svaeTn9OsOPstC4wTn8enAOz365hRnfbLM7nIhmaR+5iOQAfYFvrNyvreq6Ay6QU6I1srW/Mf8IfxqYS2lsrTpxneFHBcA9U3/HsIObue/91SxY72WuVxVwliVyEWkKvA3cYowp9LB+gogsEZElBQUO+4V7m2QhkFOiNaK1v/twCb967lvm9B3B0cen6Qw/qn5+3isQg+HxjR/SpVVTbnhlKcu2H7QsNOU7SxK5iMTiSuIzjDHveNrGGDPdGDPAGDMgPd3HG1JCXSCnRPOltV+jD706uz0v/vbvFJZW8sI1p5P2m2t0hh8VFM2qynnp2jNIbxbPNf9ZzNo9J7XjwkeIXqC1ompFgOeANcaYqf6H5DCBmhLNh9HpavahR+3Yzu/feJh3krfQMyPZmhiU8lGr5gm8Mn4gCbFR/PK5b9m2v8jukCKKFS3yIcAvgXNFZLn76wIL9nuiUKzXDqRGjE6XVFlG139NCXKgKqw1oOslKzWJl8cPpKKqmrHPfMP2/cX1v0hZwoqqlS+NMWKM6WWM6eP++sSK4H4UqvXagVZXa19LDFUI6tq6Ga+MH0hReSVX/XshmwuOBuZAdXVxRFqjD6fc2aljI5+kKtPLWOJaYqhs1jMjmdeuH0RFVTVX/XsRG/YeCd7BI7TR54xErq3PE+QdLOaBs35JSYyWGKrQ1L1tc974zSCiBK7690KWbgtSNUuENvqckcgDWa/tMGv3FHLFtK95s+vZ7Hr4X1piqEJW51bNePM3Z5KcGMvYZxbx6ardx1cGqvsjQht9zkjkgazXdpA5P+zliqe+BuDNG86k0y0TtMRQhbSctCa8/X+D6dGuOf83YxnP/m8zJpDdHxHa6HNGIg9kvbYDGODJLzYy4eUldGrVlPd/exantGlud1hK+aRl03heu34Qo3u04W8fr+HgLbcHrvsjQht9zhmPvK75BYPBpnGND0fHc1fH0Xwyax0X927HQ1f2IiE22pZYlDrBse6RsjJX98jkyV4bVwmx0TyZ249/zt1AyoN7PO/Piu6PY8cfP94VV/v2dcYVLpyTyAPNn4HoA5Tkl20/yE29xrE3til3nX8KE87p2LAJkyNkUH1lg0ZMuBIdJdw6sisl7TJI3JV38gZWdX8EqtHXgDeuYHNG10qEqaiq5skvNvKzpxcixvDmD6/xm6GddNZ7FTr8qA5JfGgK1YmJJywziSHe/RHiZY2ayEPMqp2HufTJr3h41jrO69mGj79/iX5Hd9f/QqWCyZ/qkNxcop55BhMfjwF2Nk/nTxfdzPzTRx3fJtTGNAnxskbtWgkRRWWVPP75Rp7532ZaJMUxLbcf55/WFsbOsjs0pU6Wne1qlXpafkxd3ZW5uYi7+6Pg5ff435vLefH5b7nwtLbcc1F32loecAN4ijvEyxq1RW6z6mrDm0t28JNH5vH0/E1c3jeDubcOdSVxpUKVhdUhfbJSmHnz2Uwc2ZXP1uxl+KPzmdbuDErFSzvTjlvwQ7ysURO5TYwxfLEun4uf/JI73lpJRotE3rlxMA//rDfJSbF2h6dU3SwuCY6Pieam4V347NahDO7UkgezhzK073W8+s12Kqqqj29oV191iJc1atdKkBlj+GxNPo9/voGVeYfJSEnkn2P6cHHvdo27mKmVKcouAagOyUpN4tlxp/PNT6/moexzuPvd75m+YBMTzunE5f0ySKirrzqQFSQhXtaoiTxIisoqefe7nby0cCvr9x4lOzWJh67oxWX9MoiN1g9GStU08Egeb61+lblPvcE/5q7n7ne/Z+qcdSzevh2PzZ1g9FX78sblTxmzHzSRB5AxhtW7CnlraR5vL83jSFklPTOa8+jPenNJn3bEaAJXyisBRvRozfDurVi4aT/T/7eZnc3SyCz0MFVkiPRV28VZidwh3Qhb9xXx4YpdvLd8J5sKioitruLC/ev41f3X0TcrRevBlWoAEWFw5zQGd05jd8kDlN/6O+LKSn9cXx6fQN4f7ian2hAV5cP/lk2t5kByViIPURVV1SzeeoAv1uYzd20+mwtc01yd0SGV8Wd15IK7J5BSVQrZt9scqVLO1vbG8ZCcgHH3Ve9v2YbJQ67m3d0ZtJ4yl1E92nB2lzQGfTqH5gmRUzSgiRwafOttcXkly7cf4tutB1i89QDfbT9EcXkVcdFRDOyYyq8GtWfkqW3ISHHfvVZV6nVfSqkGqlGDnjZvHn8urWDomnw+XbWHt5bm8fKibUQJ9MpMYUjnlvTLbkGvzBTSm8XXs2Pn0kRex5gRZuxY8o+UsW7PEX7YXcgPuwpZs7uQTQVHqTYgAqe0ac7P+mcyuHMaZ3VOo0m8nlKlgql5QiyX9s3g0r4ZlFdW8932g3y1cR9fbdrP0/M3U1VtAMhISaRPVgrd2w2ic8l+uhQcpX1qUlhcq7Ik64jIecA/gWjgWWNM6M0A7KFfrLKqGu66ixgP5Ux7b5rIsLUtKamo+nFxRkoi3ds247yebeiX3YJ+7VuQnOjnx7cw7K9Tyife/ub9GJwqLiaKgR1bMrBjS27F9el59a5CVuw4xHL318fZZ7s2fnQ+sdFCh7QmZKcmkZGSSEaLRDJSu5JRVkjbwlJSm8Q5oqrM70QuItHAk8BIIA9YLCIfGGN+8Hff3lRWVVNSUUVpRTWlFVWUVboeu5a5HheVVVJYWsHh4grX946jKYxO4PD0RRwqqaDgSBn7i8rYtMPDKGxAq4P5jB2YTU7LJDqlN6VHu+akJMUF6kdSKrw0Nhk3YlTFuiTFxXB6Tiqn56T+uOzouSPZlNiSjX95hA35R9mYf5S8g8V8s/kAP1k2hwsWvEe7wn3syvqcief8inn9R5LWNJ7UJnG07HoJqRUlNPtkDU3jY2gSH0Oz+BiaJrgeN23ajiZV5cTvKyI+Jsr1FRtNfEwUMVESsEIHK1rkZwAbjTGbAUTkdeASwPJE/qcPVvPKom1Uuj8q+SoxNprmKR1IriyleVU17ZIT6JOVTHqzBIpfaEfTPTtPeo20z+bei3pYFbpS4clTq9qfZByEG36aVlfQu2gPvftnnhS3eXgaUuI6fmZhAY/OeZL3u6TzRdeRHDhaztaEFixtmsHRhVsprag+eec93TE+Mu+kVVHiuoP137/szzld0y35WY6xIpFnADtqPM8DBtbeSEQmABMAshtZ83lmp5Y0iY8mISaahNhoEmJd73YJsdEkxES5l0WTGBtNYlw0yYmxNE+MIT4m+ngXxr/mnbjTRx50/ZHV/OMJoVtvlXIcf5KxnYNTTZr0YxI/JraslCvfmcaVU//oWjDsTtf3efOorKqmqKyKI2UVHC2rpKiskiO/v5Xi6DjK772fssoqyiqrKauoPv64spqMFolYzYpE7umzwklNZmPMdGA6wIABAxrWpHYbfWobRp/apjEv9S7Eb71VynH8Sca+jKoYKL7EXeMTSEx0FMlJUSeOjXR4q+t73wzLw6uLFb34eUBWjeeZwC4L9hs8ubk6ibFSVvFnpEA7B6fyd4RDO0ZldLMikS8GuohIBxGJA8YAH1iwX6WUE/mTjO2caN2fuG2eQcjvRG6MqQR+B8wC1gBvGmNW+7vfRgm1WUXA1ndppWzhbzIO5Cfkuv4f/Ynb5hmELKkjN8Z8Anxixb7CisWlVEo5RqAmQD6mMfv05f+xsXHbPINQ6Fe6O1mIz/OnVEQJ5P+jzTMIaSIPpBCf50+piBLI/0ebZxCKjERuVz91fe/S2n+uVPAEstVs50VaIiGR23k1ua53aZuvcisVcQLdaraxjDn8E7md/dR1vUtr/7lSwWVzqzmQwn/MVV/7xQI1+qC3q+Daf65U8AW6msYm4d8it/lqslehGpdSynHCP5HbfDXZq1CNS6lQMG9eWLWYAy38u1ZCdVCsUI1LqUjnwDeQ8EnkdQ1kH6r9YqEal1JW0L/poAmPrhUt5VNKRbDwSORayqeUimDh0bWipXxKqVBgU3dSeLTItZRPKRXBwqNFPnlyaM+7qRd9lAodYfj/GB4t8jC+9VYppeoTHi1y0FI+pVTECo8WuVJKRTC/WuQi8jDwU6Ac2ARcY4w5ZEFckUM/PSil/ORvi3wO0NMY0wtYD9zlf0gBomM3KKXClF+J3Bgz2xhT6X66CMj0PySllFINYWUf+bXATAv3p5RSygf19pGLyGdAGw+rJhlj3ndvMwmoBLwObiIiE4AJANl6o45SSlmm3kRujBlR13oRGQdcBAw3xpg69jMdmA4wYMAAr9sppZRqGH+rVs4D/ggMNcYU17e9Ukop6/nbR/4E0AyYIyLLReRpC2JSSinVAH61yI0xna0KxBJaXqiUikB6Z6dSSjmcJnKllHI4TeRKKeVwmsiVUsrhNJErpZTDaSJXSimH00SulFIOp4lcKaUcThO5Uko5nNQxzlXgDipSAGxr5MvTgH0WhmMVjathNK6GCdW4IHRjC8e42htj0msvtCWR+0NElhhjBtgdR20aV8NoXA0TqnFB6MYWSXFp14pSSjmcJnKllHI4Jyby6XYH4IXG1TAaV8OEalwQurFFTFyO6yNXSil1Iie2yJVSStWgiVwppRwuJBO5iPxMRFaLSLWIDKi17i4R2Sgi60RktJfXp4rIHBHZ4P7eIgAxvuGe3m65iGwVkeVettsqIt+7t1tidRwejvcnEdlZI7YLvGx3nvscbhSRO4MQ18MislZEVorIuyKS4mW7oJyv+n5+cfmXe/1KEekXqFhqHDNLRL4QkTXuv/+bPWwzTEQO1/j93hfouNzHrfP3YtP56lbjPCwXkUIRuaXWNkE7XyLyvIjki8iqGst8ykV+/z8aY0LuC+gOdAPmAQNqLO8BrADigQ7AJiDaw+sfAu50P74TeDDA8T4K3Odl3VYgLYjn7k/AbfVsE+0+dx2BOPc57RHguEYBMe7HD3r7nQTjfPny8wMXADMBAQYB3wThd9cW6Od+3AxY7yGuYcBHwfp78vX3Ysf58vA73YPrhhlbzhdwDtAPWFVjWb25yIr/x5BskRtj1hhj1nlYdQnwujGmzBizBdgInOFluxfdj18ELg1IoLhaIsBVwGuBOkYAnAFsNMZsNsaUA6/jOmcBY4yZbYypdD9dBGQG8nj18OXnvwR4ybgsAlJEpG0ggzLG7DbGLHM/PgKsATICeUwLBf181TIc2GSMaewd434zxiwADtRa7Esu8vv/MSQTeR0ygB01nufh+Q+9tTFmN7j+OYBWAYzpbGCvMWaDl/UGmC0iS0VkQgDjqOl37o+3z3v5KOfreQyUa3G13jwJxvny5ee39RyJSA7QF/jGw+ozRWSFiMwUkVODFFJ9vxe7/6bG4L0xZcf5OsaXXOT3uYtpdHh+EpHPgDYeVk0yxrzv7WUelgWsftLHGH9B3a3xIcaYXSLSCpgjImvd79wBiQuYBvwV13n5K65un2tr78LDa/0+j76cLxGZBFQCM7zsxvLz5SlUD8tq//xB/Vs74cAiTYG3gVuMMYW1Vi/D1X1w1H394z2gSxDCqu/3Yuf5igMuBu7ysNqu89UQfp872xK5MWZEI16WB2TVeJ4J7PKw3V4RaWuM2e3+eJcfiBhFJAa4HOhfxz52ub/ni8i7uD5G+ZWYfD13IvIM8JGHVb6eR0vjEpFxwEXAcOPuHPSwD8vPlwe+/PwBOUf1EZFYXEl8hjHmndrrayZ2Y8wnIvKUiKQZYwI6OJQPvxdbzpfb+cAyY8ze2ivsOl81+JKL/D53Tuta+QAYIyLxItIB1zvrt162G+d+PA7w1sL31whgrTEmz9NKEWkiIs2OPcZ1wW+Vp22tUqtf8jIvx1sMdBGRDu7WzBhc5yyQcZ0H/BG42BhT7GWbYJ0vX37+D4BfuasxBgGHj31EDhT39ZbngDXGmKletmnj3g4ROQPX//D+AMfly+8l6OerBq+fiu04X7X4kov8/38MxtXcRlz9vQzXu1QZsBeYVWPdJFxXeNcB59dY/izuChegJTAX2OD+nhqgOF8Abqi1rB3wiftxR1xXoFcAq3F1MQT63L0MfA+sdP8xtK0dl/v5BbiqIjYFKa6NuPoBl7u/nrbzfHn6+YEbjv0+cX3cfdK9/ntqVE8FMKazcH2kXlnjPF1QK67fuc/NClwXjQcHIS6Pvxe7z5f7uEm4EnNyjWW2nC9cbya7gQp3/hrvLRdZ/f+ot+grpZTDOa1rRSmlVC2ayJVSyuE0kSullMNpIldKKYfTRK6UUg6niVwppRxOE7lSSjnc/wPDocbg8XtGgwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from scipy.optimize import curve_fit\n", "init=[12,1,2]\n", "fit=curve_fit(func,x,y,p0=init,sigma=sig*np.ones(len(x)))\n", "print('parameters:' ,fit[0])\n", "print('covariance matrix: \\n',fit[1])\n", "\n", "# plot the data and the best fit\n", "plt.errorbar(x,y,sig,fmt='ro')\n", "xfit=np.arange(-10,10,0.01)\n", "plt.plot(xfit,func(xfit,*fit[0]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now calculate the $X^2$ statistic for the data set and fit:\n", "$$X^2 = \\Sigma {(y_i - f(x_i, parameters))^2 \\over \\sigma_i^2}$$" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "46.80136673096389\n", "0.8703309955431144\n" ] } ], "source": [ "x2=((y-func(x,*fit[0]))**2/sig**2).sum()\n", "print(x2)\n", "from scipy.stats import chi2\n", "print(chi2.cdf(x2,len(y)-len(fit[0])))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use scipy.stats.chi2() to determine the probability of getting a $\\chi^2$ smaller than and larger than your value. Would you rule this fit out as a good fit?\n", "
ANSWER HERE: " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Change the size of the uncertainties going into the simulation but not the uncertainties that you feed the fit (or vice versa). Redo the fit , calculate $X^2$ and the probability again.\n", "
DISCUSS YOUR RESULTS: " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now add a second Gaussian to your simulated data, with different parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x=np.arange(-10,10,0.5)\n", "par=[10,0,3]\n", "sig=1\n", "y=func(x,*par)+func(x,10,6,3)+np.random.normal(0.,sig,size=len(x))\n", "plt.errorbar(x,y,sig,fmt='ro')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit it with a single Gaussian" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "parameters: [13.26859434 3.13414992 4.67786386]\n", "covariance matrix: \n", " [[ 0.13120967 -0.00599569 -0.03835598]\n", " [-0.00599569 0.02497582 0.00693781]\n", " [-0.03835598 0.00693781 0.03102229]]\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD6CAYAAACiefy7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xUVfr48c9DKgECBEInoSlVaRGwguuKimJZdRWCoCJYUH+6X11RXOtiYdVdFZWiuBbAtaOICqsC6ooYOkgLPdTQS0IgyfP7YyY6hJlkkqmZed6vV16ZuffcOSd3Js/c+9xzzxFVxRhjTOSrFuoGGGOMCQ4L+MYYEyUs4BtjTJSwgG+MMVHCAr4xxkQJC/jGGBMlvA74IjJJRHaJyHKXZY+JyFYRWez86edh24tFZLWIZIvISH803BhjTMWIt/3wReQ84DDwtqp2ci57DDisqs+VsV0MsAa4EMgBfgEGqOqv5dVZv359bdGihVftM8YYAwsWLNitqqnu1sV6+yKqOldEWlSi/h5AtqquBxCR94ArgHIDfosWLcjKyqpElcYYE51EZJOndf7I4d8pIkudKZ+6btY3Bba4PM9xLnNLRIaLSJaIZOXm5vqhecYYY8D3gP8a0BroAmwHnndTRtws85hHUtUJqpqhqhmpqW7PSowxxlSCTwFfVXeqapGqFgMTcaRvSssBmrs8bwZs86VeY4wxFedTwBeRxi5PrwKWuyn2C3CKiLQUkXjgeuAzX+o1xhhTcV5ftBWRqUAfoL6I5ACPAn1EpAuOFM1G4FZn2SbA66raT1ULReRO4GsgBpikqiv8+lcYY4wpl9fdMkMhIyNDrZeOMcZ4T0QWqGqGu3V2p60xxkQJC/jGGBMlLOAbEyn69HH8GOOBBXxjjIkSFvCNqUrsKN74wAK+McZECQv4xpjy2ZlFRLCAb4wxUcICvjHGRAkL+MYYEyUs4BtjTJSwgG9MJJg8GebNgzlzoEULx3NjSrGAb0xVN3kyDB8OBQWO55s2OZ5b0DelWMA3pqrwdBQ/ahTk5Z1YNi/PsdwYF16Ph2+McSrpjz57dvDq9HQUD7B5s9tNdPNm/pe9m+xdh9l+4Cj7jhxjX94xihWqCcTFVqN+jXhSayXQpE512jVKpnWDGiTExgTpjzLBZgHfmKqgrKP4tDTHF0ApW2vVJ/P1nwGIixHqJsVTJymOmGrVUFWOFRaz+3ABB48W/rZNbDWhY5NkzmpTn3Pa1KdHyxTi3pvqOLMoKHCcWYweDZmZgfxrTYBYwDemKijjKP77R16gx9MjSTxW8NvygvhE1t3zEFNu6EmbBjVJrZWAnH++Y2WpM5Ojx4vYsjePVTsOsXL7QeZv2MvEuet5bfY6Bqz9nse/eJF4d2cW3gb9UJwRGbcs4Bvjb4EIcB6O4rclpzL46Cnceu19jPjkZWrlHYS0NBKeeoreXgbkxLgYTmlYi1Ma1qJ/5yYAHC4o5H/Zu8noM4z4gqMnblByZhGMo3z7svArC/jGBFNlA9jo0Y4ja5e0Tn5cAnOH3MP0u86hU9NLoc8PlXttN2omxNK3YyPYs8Ptet28maPHiqgeb/n+qsR66RhTBWy++ComDX6QnORUihEONmhC8fjxDHjxQTo1rR24itPS3C7eWqs+5475jrd/2khhUXHg6jd+ZQHfmDCWf6yI52eu5o//nMNzqRm8338Yh/5wIck7t1LjpiGBb8Do0ZCUdOKypCSOPf532jSowSPTVnDZyz8wf8PeytdhI3EGjdcpHRGZBFwG7FLVTs5l/wD6A8eAdcBNqrrfzbYbgUNAEVDoaUZ1Y8zvfli7mwc+WsrW/flc0aUJD17SnkZXPBOYyjylmkry9EOHOnrppKfD6NG0ysxkqipfr9jBk9NX8ufxP/Gnbk159LKO1E6KC0wbjc8qksP/NzAWeNtl2SzgQVUtFJFngQeBBzxsf76q7q5UK42JIkcKCnlqxkom/7yZVqk1+M/wXvRsVS90DcrMhIkTHY9dvhBEhIs7Nab3qQ0Y+91axs9Zz4/ZuxlzTWd6n5oamraaMnmd0lHVucDeUstmqmpJJ955QDM/ts2Y8BPgMWsWbd7HxS/OZcr8zQw7tyUz7j43tMHeC9XjY7j/onZ8csfZJCfGMWTSfEZ9soz8Y0WhbpopxZ85/JuBLz2sU2CmiCwQkeFlvYiIDBeRLBHJys3N9WPzjPFRAMesUVUm/bCBP4//CVV4/9YzGXVpBxLjqk4vmNOa1ebzu85h+HmtmPzzZq569Uc27D4S6mYZF34J+CIyCigEPH3yz1bVbsAlwAgROc/Ta6nqBFXNUNWM1FQ7LTRhJEBj1hw8epw7Ji/kiem/0qdtA76461zOaJHi02ueJEijaSbGxfBQv/a8edMZ7Dh4lHHDnyB//gIbxTNM+NwPX0SG4LiYe4GqqrsyqrrN+XuXiHwC9ADm+lq3MUHl4W7XE5aXBFYvhyHYtOcIN//7FzbuyeOhfu0Ydm4rRKRy7fPU/76scXgCdPPU+W0b8E2jrdSY/iKJxytZbwX3pSmfT0f4InIxjou0l6tqnocyNUSkVsljoC+w3Jd6jQkJD33Sf1teXsqn1FF29r8mcOUrP7LnyDHeHdqT4ee1rnywL0uIRtOs99Tjvwf7itZrQz4HhNcBX0SmAj8BbUUkR0SG4ui1UwuYJSKLRWScs2wTEZnh3LQh8IOILAHmA1+o6ld+/SuMCQYPfdIZPdrxuKzA6iaANbn/bq5eNZdP7jibM1t7eWF29uyK30nrzZlJIJQx/k+5bMjnwFDVsP3p3r27GhNW3n1XNSFBFVTT0x3PS4g4lpf+EXGUdbOuqHla4NvsoW5NT/+9TO/ejh9PyltfgXq312mg63MPO8p42p9l7Utf2xXhgCz1EFPtTltjKiIzE3r1gt69YePGE3PKZaV8PBzVVsvZ4v82llbemYk3F3Qrc2bhpt7i6tV56YKbuHbcT+SMfd1z2qa89JmpFAv4xvhLGYFVmzd3v00wAlhmJkyYAAkJjufp6Y7nmZmBzZW7qbfaxInc/OrDxMUI1R4uI21T3peUqRxPh/7h8GMpHROWykojuElRHCss0jduf1KPxCacmJ5ISjoxJRSKdnuT7glAvVv2HtEiyknblJU+K+tvinJYSseYICmV8jl+/QDunrqIJ5K7Mu/Bp90fZYdSiC7oNqubVP5ZT1npM1MpNh6+MQFSWFTMPe8t5svlO/jbZR34wzmXwtxpjpXhMqGHh4lVgpFqinn6KXTYcCTfJa1jaZuAsiN8YwKgEOHe95fwxbLtjOrXnqHntPRuw2APFRzKXHlmJjJxAkWJiRQDW2s3YM3fX7Aj+QCygG+MnxUh3Ne6H58v2caDl7Rj2HmtQt0kz8q6oBuk+mN69mTvBZdw4yPvc8XedN/G1jdlsoBvolOAjqRVlcdaXMCnqR24/6K23Nq7td/r8Dt/5Mp93J/1C/OYPKwnjeskMvTfv7Bi24FKv5bxzAK+MRVVRp/0F79ZyzuNunLrtvmMOL9NUJtV1TWolcg7Q3tSMzGWIZN+YdOeckbaDNKAcJHEAr4xfvLOTxv513/Xcu2uZYzcPCfUzamSmtapzjtDe1BYXMwNb8xnV1wN9wVtrJ1KsYBvjB98vmQbj3y2gj+2b8jT678mAEOg+a4yd8uGoN42DWrx5o1nsPtwAYPbXcOBmISTC9lYO5ViAd8YH/2YvZu/vL+YM1qkMHZgV2Jnf+c5wIUq6FYxXdPqMm5Qd9YlN2TYtY9RUFhq9qxQDQhXxVnAN6a0ClyAXLPzELe9s4DWqTV5fUhGlZqhKtydd2oqz13bmfkb9/LXD5fiuInUyduxdoLdzTXMWcA3kSkI/+i5hwq46c1fSIyP4Y0bzyA5MS6g9UWjK7o05b6+pzJt8Tb++d+1v6+wsXYqxe60NaYS8o8VccvbWew9coz3bz2TpnWq+/6i0TrDUzkprhHnt2HTnjxe+mYt6SlJXN292e/7ZehQx/5KT4+e/eUDC/jGVFBxsfKX9xezNGc/4wd157RmtX1/0RBMQ3iCML6uICKMvuo0tu7PZ+THS2lSp7pjwpjMTJg40VEojNsfTiylY6KPj/23x3y9mi+X72BUv/b07djIP22qyr1OgtAfPj62Gq8N6k56vRrc+k4W2bsO+72OaGAB30QXH/tvf5C1hXFz1jGoV5r34+N4o6r2Oglif/ja1eN488YziIupxrC3sziQf9zvdUQ6C/gmuvhwJL1w8z5GfbKcs9vU47H+Hf074XhVneEpyGcmzVOSGHdDd3L25XH31EUUhecdD2HLAr6JLpU8kt558Ci3vbOAhrUTGDugG7Exfv7Xqaq9TkJwZnJGixQev7wTc9bkMibt3IDVE4ks4JvoUokj6YLCIm57dwGHCwqZODiDujXi/d+uUI9aWVkhOjMZ2DONQb3SGN+kJ9PqtQ9oXZGkQgFfRCaJyC4RWe6yLEVEZonIWufvuh62HeIss1ZEhvjacGMqpYITeuvkyfzt0+Us2ryf56/tTLtGyYFrW1Wc4SmEZyaPXNaRHge38NfWF7Esx83omja42sk8zX3o7gc4D+gGLHdZNgYY6Xw8EnjWzXYpwHrn77rOx3XLq8/mtDWVUt5cqJ7Wv/uuY55Zl/lVjydU17su+z99/utVwWl7VZyj1Zu5ZwMk94JL9Kx7Jmuvp/6ruw4ePbFNpd7Lk+YQror72guUMaetqOvtyl4QkRbAdFXt5Hy+GuijqttFpDEwW1XbltpmgLPMrc7n453lppZVV0ZGhmZlZVWofSbKlfQaySs1bV7p9EjJXbiu/bdbtHA73d/ueo1I2bWNatWCcIHQXbuqghC2e/nWA1wz7n+c1rQ2U4b1Ii6mmsf3kvR0x9kTVN19XQ4RWaCqGe7W+SOH31BVtwM4fzdwU6YpsMXleY5z2UlEZLiIZIlIVm5urh+aZ6KKL71GPFxorLd3Z3CCvamUTk1r8+zVp/PLxn08PWOVY2FV7eYaYMG6aOvuv8XtqYWqTlDVDFXNSE1NDXCzTMTx5R/dw4VGCfeukYYrujTlxrNaMOnHDUxfuq3qdnMNMH8E/J3OVA7O37vclMkBmrs8bwZs80PdxpzIl3/0qto10gDwUL/2dEurwwMfLmXHA3+z99INfwT8z4CSXjdDgGluynwN9BWRus5ePH2dy4zxL1+CdmYmKx5/jpzkVMfpZ1XpGmkAx/ALr2Z2JzEuhkH5bTj66riq1801wCraLXMq8BPQVkRyRGQo8AxwoYisBS50PkdEMkTkdQBV3Qs8Cfzi/HnCucwY//KhP/v2A/kMymvNLYOfpaDPBVWna6T5TaPaibw8oCvrcw9zX1xHtKp1cw2wCo2WqaoDPKy6wE3ZLOAWl+eTgEkVap0xlVGJURQLi4q5a8oijhUW8+raz0jUwsC1rywR1mMkFM5qU5/7L2rHs1+tonujbty0Y2HlXigCe/HYnbbGAM/NXEPWpn089afTaHV0X6ibY3x0W+9WXNihIaPT+pBV022HwKhkAd9EJ5e5Zb9btYtxc9YxsGcaV3Sx4FApYTZXr4jw/J8706zgICNO7c/uwwUnFojSu3At4Juotm1/Pn95fzHtGyfzyGUdQt0c40fJiXG8tnYa+2MTufc/iykudvYED+KQzuHGAr6JWseLirlrqjNvn9nNJiCPQO3zcnl8wzd8v3Y3r87OdiysypPN+MgCvolaz81czYJN+3j66tNpWb9GqJtjAuS63GVc0aUJL8xaw8/r90T1XbgW8E1U+nbVTsbPWU9mzzQu79wk1M0xASTA6KtOo0W9Gtz93iKKmjV3XzAK7sK1gG8iUxkXEbfuz+cv7y+hQ+Nk/uYubx9mFyCND5zvZc2EWF7J7Mb+vOOMv+hmNErvwrWAb6LK8aJi7pqykMIi5RXL20eV9o2TebR/R8bUy2DWPU9G5V24Fbrxypiq7rmvV7Nw835eHtDV8vZRaECP5sxbv4fblsJ7F15Nj0Nb3Z/NlXTbLChwdNscPToivhDsCN9EjW9W7mT83PUM6pVGf8vbRyUR4ak/nUZ6vRrcdUp/9sRWP7lQBHfbtIBvokJJ3r5jk2QevtT620ezmgmxvDKwG/tiq3Nvm0t/759fwptum336/D70QhViAd9EvONFxdw5ZSFFxcorAy1vb6BDk2Qe3fgtc+u05LU5605cGcHdNi3gm4g35qtVLNq8n2evPp0Wlrc3TgN3LaH/7pW8MGsN8ze4DN4bwZOnWMA3EW3WrzuZ+P0GBp+ZzqWnNw51c0wYEeCpDTNJS0nirqkL2VMy3k4ET4RjAd9ErJx9edz3wRI6NU3moX7tQ90cE4ZqFR1j7MCu7Ms7zv99sMSRz/dhToVwZwHfRKRjhcXcOWURxZa3N+Xo2KQ2f7usA7NX5zJ+7nrHwsxMiMDJU6wfvolIY75axeIt+3k1sxvp9Sxvb9xw6X8/qGca89bt4bmZqzmjRV0yWqSErl0BZEf4JnQC1LVt5oodvP7DBoacmU6/0yxvb8onIjx99Wk0rVOdu6YuYt+RY6FuUkBYwDfhqxJfCFv2uuTtL7W8vfFecmIcrwzsxp7Dx7jvgyVo+ZtUORbwTcQ4VljMnVMXoQqvDOxGQqzl7U3FnNasNqMubc83q3bxeuMM94Wq8GxZlsM3EePZr1axZMt+XrO8vfHB4DPT+WndHp4tPo/uh7bSzXWlp2EXoEpc2PX5CF9E2orIYpefgyJyT6kyfUTkgEuZR3yt1xhXX6/YwRs/bODGs1pwieXtjQ9EhGevOZ3Gxw5x1yn92Z/nks+v4rNl+RzwVXW1qnZR1S5AdyAP+MRN0e9LyqnqE77Wa0yJLXvzuP+DJZzerDYP9msX6uaYCFC7ehxj137Orria3PfBUlSdGf0qPuyCv3P4FwDrVHWTn1/XGLcc/e0XosDYAZa3N/7T+Yv3GHn5afx35U4m/bjRsbCKD7vg74B/PTDVw7ozRWSJiHwpIh09vYCIDBeRLBHJys3N9XPzTKR5+suVLMk5wD+uOZ20eknlb2BMBdx8dgsu7NCQZ75cyeIt+6v8sAt+C/giEg9cDnzgZvVCIF1VOwMvA596eh1VnaCqGaqakZqa6q/mmQj01fIdvPnjRm48qwUXd7K8vfE/EeEf15xOg1qJ3DllIQf+9OcqPeyCP4/wLwEWqurO0itU9aCqHnY+ngHEiUh9P9Ztqhofu7Zt2ZvH/R8uobPl7U2A1UmKZ+zAruw4cJS/frgEHTiwyg674M+APwAP6RwRaSQi4nzcw1nvHj/WbaoSb2YUKuMLoaCwiBFTFgIw1vrbmyDomlaXBy5ux9crdvLW/zaGujmV5peALyJJwIXAxy7LbhOR25xPrwGWi8gS4CXgev3tsreJOuV1bSvnC+GpL1ayNOcAz13bmeYplrc3wXHLuS25oF0DnpqxiqU1Goa6OZUi4Rx3MzIyNCsrK9TNMP5WrRq4+9yJQHGx44h+k5uOXunpTP/sf9w5ZRG3nNOShy+zqQpNcO07coxLX/qe2J07mL7sLZK/mRnqJp1ERBaoqtvbhG1oBRN85XVt89CnWTdvZuRHy+iWVocHLrG8vQm+ujXieXlgV7YmJDOy1cWE8wGzOxbwTfCV17XNwxfCrjoNiIsRxg7sRlyMfXRNaHRPT+H+zd8zo15b3p1XtW45sv8aE3zlzSjk5gvhWHwio88exD+v60KTOtUdCwM0vLIx5Rm+fT599q3nyekrWb71QKib4zUL+CY0yppRqNQXwpFGTbmv7wjS7hxGn7YNQtNeY1xUA15YN4OUGvGMmLKQQ0ePh7pJXrGAb8KT8wth9cVX0/2WiezqfzX3/PGUULfKGIfZs0n575e8NKArOfvyefDjZVUin28B34StI9XiuOPUy6mZEMdL13cl1vL2Jsz0aJnCXy48lelLt/Puz+E/gJr9B5mwpKo81KovGxLr8tKALjRITgx1k4xx6/berenTNpUnPl/Bos37Qt2cMlnAN2Fp8s+bmVa/A/fm/MhZrW0UDhO+qlUT/nVdFxomJ3LH5IXsPlwQ6iZ5ZAHfhJ2Fm/fx+Ocr6LNvPSO2zgt1c4wpV52keMYN6s7eI8e4a8oiCouKQ90ktyzgm7CSe6iA299dQOPa1Xkxe7p9QE2V0alpbf5+ZSd+Wr+H52auCXVz3LL/JxM2jhcVM2LKQg7kH2fcoO7ULgrfU2Nj3Lk2ozkDe6Yxbs46vlq+PdTNOYkFfBM6s2c7fpye+XIV8zfs5Zk/nU6HJsllb+vj8MrGBMqj/TvQuXkd7vtgKetyD4e6OSewgG/CwrTFW3+bhPzKrk0dC0t9IfzGm+GVjQmRhNgYXsvsRnxsNW57ZwFHCgpPLhSiu8Qt4JuQW7XjICM/WsYZLeoy6tL25W9Q3vDKxoRYkzrVeXlAV9blHuaBj5aGzU1ZFvBNSB3IP86t7yygVmIsr3g7KJqH0TQ9LjcmBM5uU5/7L2rH9KXbef37DaFuDmAB34RQUbFy738Ws21/Pq8N6ub9zVXlDa9sTJi4rXcrLunUiKe/XMmcNbmhbo4FfBM6Y75exberdvFI/450T0/xfsPyhlc2JtScOXoR4fk/d6Zto2TunLKQ9SG+iGsB34TExwtzGD9nPZk907ihV3rFNi5veGVjwkhSfCwTB3cnLqYat7ydxYH80I2saQHfBN2izfsY+fEyerVK4bHLO1buRcoaXtmYMNOsbhKvZXZj85487p66iCIkJO2wgG+CaseBo9z6zgIaJifwamZ3m7nKRI2ererx+BUdmbMmlzFp54WkDbEhqdVEpaPHixj+ThZHCgp595azSakRH+omGRNUmT3TWbX9EOPnQdu8XP4U5Pr9dnglIhtFZJmILBaRLDfrRUReEpFsEVkqIt38VbcJf6rK/R8uZdnWA7x4fVdObVgr1E0yJiQePbSYX14ZzJUfjeNYs+ZBvWHQ30f456vqbg/rLgFOcf70BF5z/jZR4JXvsvl8yTb+enFb/tihYaibY0xoTJ5M7G23kuq8cTB+aw7Fw4Y5jryDcB0qmAnUK4C31WEeUEdEGgexfhMi0xZv5bmZa7iqa1Nu79061M0xJnTc3CVeLT+f4gcfCkr1/gz4CswUkQUiMtzN+qbAFpfnOc5lJxCR4SKSJSJZubmhv1HBlKOcMUF+Xr+H+z9YSs+WKTxz9WmIhKZ3gjFhwdPd4Fu2cKww8GPo+zPgn62q3XCkbkaISOnL0O7+008aYEJVJ6hqhqpmpKam+rF5JtjW5R5m+DsLaJ5SnQk3ZJAQGxPqJhkTeGWN5OrhbvBtyfV56JPAT4Tut4Cvqtucv3cBnwA9ShXJAZq7PG8GbPNX/Sa87D5cwE1v/kJcjPDvm3pQOynO/5V4Gk3TmFApbyRXD3eJL7n9AT5ckMPYb7MD2jy/BHwRqSEitUoeA32B5aWKfQYMdvbW6QUcUNXwmyHA+Ozo8SJueSuLnQePMnFwBs1TksrfyJhIUN5Irh7uEu/39F+4qmtTnp+1hmmLtwasef7qpdMQ+MSZn40FpqjqVyJyG4CqjgNmAP2AbCAPuMlPdZswUlSs3PPeYpbk7Oe1zG50HXyVY4UdiZto4M1IrpmZMHGi47Hz/0KAZ64+ja3787n/g6U0qJXIma3r+b15fjnCV9X1qtrZ+dNRVUc7l49zBnucvXNGqGprVT1NVU/qq2+qNlXl4U+X89WKHTx8aQcu7mSdsEyU8WEk14TYGCbekEF6vSRGTFnofuIUH9l97cZvnp+5hqnzN3N7n9YMPadlqJtjTPD5OJJr7aQ4/n1zD166vis1Evw/EIINrWD8YtIPGxj7XTbXn9Gcv17UNtTNMSY0Sm6eGjrUceE2Pd0R7CtwU1XTOtVpWqd6QJpnAd/47NN67Xli+q9c1LEhf7+yk/W1N9HNTY4+XFjAN5U3eTL587O4PH8uZ6d8Su1/jiE2JiPUrTLGeGAB31TO5MkU3TKM6kfzAUjduwNuvw1iqtnY9MaEKbtoGw3KGf6gMtsXPDCSGGew/41rf2NjTNixI3xTYQs376PLVg83h3jqh2yM+V2Icvt2hG8qZFnOAYZMms+uOh7GOXLtb1zWmCLGmKCzgG+89uu2gwx642dqV48jfswzZfc3Lm9MEWNM0FlKx3hleVIDBr0+j6T4GKYO60VKyh8gKd5zf+OyxhSxi7om0oVZd8wSFvDN7xdkPXxIF9Rswo3triY5PpYpw3r+PhhaWf2NvRlTxBgTVJbSiRS+9sTx4Kd1e7ih/bXUO57PB7edSXq9Gt5t6MOYIsaYwLCAbzyasyaXG9+cT9OCg7z/61SaVOR2bx/HFDHG+J8FfOPWV8t3MOytLFqn1uS9X9+jwfEjFXsBD+N+W/7emNCxgB/pKtE18t15m7hj8gI6Nk3mw+QN1Pvhu8p1rczMhF69oHdv2LjRgr0xIWYXbSOZp66R4Db4qiovzFrDy99mc0G7Brwmq4i/43avtzfGhDc7wo9k5U235qIQ4YGPlvLyt9lcl9Gc8Td0J/7Rv3m9vTEm/NkRfiTzsmvk4Wpx3H1Kf77NyuHuP7Th3gtPdQxxbF0rjYkoFvAjWVqaIw3jbrlTzr48bumYyZqkejx5ZSdu6JVeoe3D9QYTY8zJLKUTycrpGrlg015eG/4Ek966n3VjLueG63ufeFHWulYaE1Es4EcCTz1xyuga+dGCHCbf/Qx/m/YvmhzajcDJ491Y10pjIoqoqm8vINIceBtoBBQDE1T1xVJl+gDTgA3ORR+r6hPlvXZGRoZmZWX51L6IV9ITx/XialLSiYHZZeiEwqJixny9mglz1/PLxFscE5eUlp7u6EZZopyhF8rky7bGmAoTkQWq6nbqOX/k8AuB/1PVhSJSC1ggIrNU9ddS5b5X1cv8UJ9xVYFBynYdPMqdUxcxf8NeBvVKo/6Yne5f058XZS3QGxM2fA74qrod2O58fEhEVgJNgdIB3wSClz1pfq7VjBEv/cDhguP887rOXNW1mXcXZY0xEcOvOXwRaQF0BX52s/pMEVkiIl+KSMcyXmO4iGSJSFZubq4/mxeZyhmkrLhYGd/4DAZ2uI7kxFimjTjHEezBLsoaExV/mFcAAA7hSURBVGX8FvBFpCbwEXCPqh4stXohkK6qnYGXgU89vY6qTlDVDFXNSE31MKuS+V0ZQXvnwaMMeXM+T6f3oe/etUy782zaNqr1ezm7KGtMVPFLP3wRicMR7Cer6sel17t+AajqDBF5VUTqq+puf9Qf1UqCc6mJSL7ucgEj/zWX/ONFjF4/k4G7liCJf3G/vacx7Y0xEcXnI3wREeANYKWqvuChTCNnOUSkh7PePb7WbZxcBinLW5PNg9VP59Z3FtCkTnWm33UumbuWOLpdGmOimj+O8M8GbgCWichi57KHgDQAVR0HXAPcLiKFQD5wvfraH9Sc5MfkNEb+ay45+/K5tXcr/u/CtsTH2q0WxhgHf/TS+QHKPoBU1bHAWF/rMu4dPHqcp1r25b2GnWkhwnvDetGzVb3fC/iaqrFUjzERwcbSqeJm/bqThz9dRm6D07h123zuffIREuNiQt0sY0wYsoBfRW3ac4THP/+Vb1ftom3DWkz48XU6H9kBFuyNMR5YwK9i8o8V8ersbMbPWU9cjDCqX3uGnNWC+GmPhLppxpgwZwG/iiguVj5fuo0xX61m6/58ruralAcvaUeD5MRQN80YU0VYwA9zqsr3a3fz7FerWLHtIB0aJ/PP67rQo2WK/yqxi7LGRAUL+GFsac5+nv1qFT9m76FZ3er867ouXN65CdWquekUZUHbGFOOyAz4VXxI3gWb9jL222y+W51LSo14Hu3fgYE900iItQuyxpjKi8yAXwWpKj+t28PL32bz0/o9pNSI5/6L2jL4zHRqJcaFunnGmAhgAT/EjhUWM2PZdt78cQNLcg7QoFYCD1/anoE900iKd3l7qvhZizEm9Czgh0juoQKm/LyZd3/eRO6hAlql1uDvV3bimu7N7MYpY0xAWMAPouJi5cd1u3k/K4evl+/gWFExfdqmctPZLTm3TX33F2ONMcZPoi/ghyA1smVvHh8syOGjBTls3Z9P7epxDOyZxg1nptM6tWbQ2mGMiW7RF/CDZPuBfGYs28GMZdtZsGkfInDuKak82K8df2zf0NI2xpigs4DvRzn78pj1606+WLqdrE37AOjQOJn7L2rLVV2b0uTKS+Bt7MKrMSYkIi/gT54M8+Y5Zn9q0cIxBWCApuw7XlRM1sZ9zF69i+9W72LNzsMAtGtUi/v6nkq/0xrTylI2xpgwEVkBf/JkGD7cEewBNm1yPAe/BP3iYmXVjkP8vGEP89bv4X/ZezhUUEhcjNCjZQp/zmjOH9o1sCBvjAlLkRXwR42CvLwTl+XlOZZXIuAfPV7Eyu0HWbh5P/PW7+GXjXvZn3ccgOYp1bn09Mac364BZ7epT82EyNqVxpjIE1lRavPmii13cbyomHW5h1m65QBLcvazJGc/q3cc4niRYybG9HpJ9O3QkF6t6tGzVT2a1ql+8osEqgdQENNUxpjIFVkBPy3NkcZxt9ypGMhJqM3qX3eyZuchVu84xJqdh1iXe/i34F4rIZbTm9fmlnNb0blZbbo0r0uj2iEahjjAaSpjTPSIrIA/erQjGLqkdQoTqjPt6jv48q1f2LQnj8097qGgWhy8nQVA0zrVObVhTXq3TaVdo1qc3qwOLevVCJ+boPycpjLGRK/ICviZmRQrbB9xL40P7mZbcn3GnDeYWUmnk74vn+vWzOWad5+n9pEDHGvajOK/j6b6jYND3eqy+ZCmMsYYV34J+CJyMfAiEAO8rqrPlFqfgKMHendgD3Cdqm70R92lVRuUyZfTF1GnMJ/0557k4XpJvFgzAZkyBUY989vRcsLWHBhxu2MO2HA+UvYiTWWMMd6o5usLiEgM8ApwCdABGCAiHUoVGwrsU9U2wD+BZ32ttyy37Mjimt0rOKNFCg1qJSIiZadGgqHkwuucOY4Lr5Mne7fd6NGQlHTisqQkx3JjjKkAnwM+0APIVtX1qnoMeA+4olSZK4C3nI8/BC4QkeAmyUOZGvF04dWboJ+ZCRMmQEKC43l6uuN5OJ+VGGPCkj8CflNgi8vzHOcyt2VUtRA4ANRz92IiMlxEskQkKzc31w/Nc/KUAim9vE+f37tX+ouvZxeZmdCrF/TuDRs3WrA3xlSKPwK+uyN1rUQZx0LVCaqaoaoZqampPjfuN6FMjdiFV2NMGPBHwM8Bmrs8bwZs81RGRGKB2sBeP9TtPX+lRjydAZSVo/fm7CIQZxbGGOPCHwH/F+AUEWkpIvHA9cBnpcp8BgxxPr4G+FZV3R7hB1SgUiPl5ejtwqsxJgz4HPCdOfk7ga+BlcD7qrpCRJ4Qkcudxd4A6olINvAXYKSv9YaV8nL0duHVGBMG/NIPX1VnADNKLXvE5fFR4Fp/1OWVYI83702OPjMTJk50PLbx8I0xIeCPlI7xtgeQMcaEkAV8f7AcvTGmCrCA76qyd8MGI0c/e7algowxPomswdN84eswxJajN8aEOTvCLxHqsXaMMSbAoi/ge0qNeNPTprIpn/IE6nWNMcZF9AV8T8rraePLAGhlCdTrGmNMKRbwS5TX08YfKR93ZxeWSjLGBIkF/BLl9bQJ1ABoNrCaMSZIrJeOq7J62gRq5imb0coYEyR2hO+tQN1cZTdtGWOCxI7wvVWS2hk61HGBNT3dEZRd++hXpv+9N69rjDF+YAG/IgJ1c5XdtGWMCQJL6RhjTJSwgG+MMVHCAr4xxkQJC/jGGBMl7KJtaXbR1BgToewI3xhjooQFfGOMiRKW0qmoQKV8LJVkjAkwnwK+iPwD6A8cA9YBN6nqfjflNgKHgCKgUFUzfKnXGGNMxfma0pkFdFLV04E1wINllD1fVbtYsDfGmNDwKeCr6kxVLXQ+nQc0871JxhhjAsGfF21vBr70sE6BmSKyQESGl/UiIjJcRLJEJCs3N9ePzTPGmOhWbg5fRP4LNHKzapSqTnOWGQUUAp7m5TtbVbeJSANgloisUtW57gqq6gRgAkBGRoZ68TcYY4zxQrkBX1X/WNZ6ERkCXAZcoKpuA7SqbnP+3iUinwA9ALcB3xhjTGD4lNIRkYuBB4DLVTXPQ5kaIlKr5DHQF1juS73GGGMqztcc/ligFo40zWIRGQcgIk1EZIazTEPgBxFZAswHvlDVr3ys1xhjTAX51A9fVdt4WL4N6Od8vB7o7Es9xhhjfCce0u5hQURyATczfHulPrDbj83xF2tXxVi7KsbaVTGR2K50VU11tyKsA74vRCQrHG/ysnZVjLWrYqxdFRNt7bLB04wxJkpYwDfGmCgRyQF/Qqgb4IG1q2KsXRVj7aqYqGpXxObwjTHGnCiSj/CNMca4sIBvjDFRokoHfBG5VkRWiEixiGSUWvegiGSLyGoRucjD9i1F5GcRWSsi/xGR+AC08T/Ou5AXi8hGEVnsodxGEVnmLJfl73a4qe8xEdnq0rZ+Hspd7NyH2SIyMgjt+oeIrBKRpSLyiYjU8VAuKPurvL9fRBKc73G287PUIlBtcamzuYh8JyIrnZ///+emTB8ROeDy/j4S6HY56y3zfRGHl5z7a6mIdAtCm9q67IfFInJQRO4pVSYo+0tEJonILhFZ7rIsRURmOePQLBGp62HbIc4ya51jmFWcqlbZH6A90BaYDWS4LO8ALAESgJY4ZuOKcbP9+8D1zsfjgNsD3N7ngUc8rNsI1A/ivnsMuK+cMjHOfdcKiHfu0w4BbldfINb5+Fng2VDtL2/+fuAOYJzz8fXAf4Lw3jUGujkf18Ix+VDpdvUBpgfr8+Tt+4LjDvwvAQF6AT8HuX0xwA4cNycFfX8B5wHdgOUuy8YAI52PR7r7zAMpwHrn77rOx3UrWn+VPsJX1ZWqutrNqiuA91S1QFU3ANk4Ruj8jYgI8AfgQ+eit4ArA9VWZ31/BqYGqo4A6AFkq+p6VT0GvIdj3waMhtekOt78/Vfg+OyA47N0gfO9DhhV3a6qC52PDwErgaaBrNOPrgDeVod5QB0RaRzE+i8A1qlqZe/g94k6hoXfW2qx62fIUxy6CJilqntVdR+O2QYvrmj9VTrgl6EpsMXleQ4n/0PUA/a7BBd3ZfzpXGCnqq71sN7rSWL86E7nafUkD6eR3uzHQPLLpDo+8Obv/62M87N0AMdnKyicKaSuwM9uVp8pIktE5EsR6RikJpX3voT6M3U9ng+6QrG/ABqq6nZwfJkDDdyU8ct+82nwtGAQLyZgcbeZm2Wl+596U8YrXrZxAGUf3Xs9SYw/2gW8BjyJ429+Eke66ebSL+FmW5/78Xqzv8SPk+r40lQ3ywL2OaooEakJfATco6oHS61eiCNtcdh5feZT4JQgNKu89yWU+yseuBz3c2+Han95yy/7LewDvpYzAYsHOUBzl+fNgG2lyuzGcToZ6zwyc1fGL20UkVjgT0D3Ml7D75PEeLvvRGQiMN3NKm/2o9/bJeEzqY43f39JmRzn+1ybk0/Z/U5E4nAE+8mq+nHp9a5fAKo6Q0ReFZH6qhrQgcK8eF8C8pny0iXAQlXdWXpFqPaX004Raayq253prV1uyuTguM5QohmOa5cVEqkpnc+A6509KFri+Kae71rAGUi+A65xLhoCeDpj8NUfgVWqmuNupYRgkphSedOrPNT3C3CKOHozxeM4Hf4swO0Kp0l1vPn7P8Px2QHHZ+lbT19S/uK8RvAGsFJVX/BQplHJtQQR6YHjf31PgNvlzfvyGTDY2VunF3CgJJ0RBB7PskOxv1y4foY8xaGvgb4iUteZfu3rXFYxgb4qHcgfHIEqBygAdgJfu6wbhaOHxWrgEpflM4AmzsetcHwRZAMfAAkBaue/gdtKLWsCzHBpxxLnzwocqY1A77t3gGXAUucHrnHpdjmf98PRC2RdkNqVjSNXudj5M650u4K5v9z9/cATOL6QABKdn51s52epVRD20Tk4TueXuuynfsBtJZ8z4E7nvlmC4+L3WUFol9v3pVS7BHjFuT+X4dK7LsBtS8IRwGu7LAv6/sLxhbMdOO6MXUNxXPP5Bljr/J3iLJsBvO6y7c3Oz1k2cFNl6rehFYwxJkpEakrHGGNMKRbwjTEmSljAN8aYKGEB3xhjooQFfGOMiRIW8I0xJkpYwDfGmCjx/wF7VQvrv6f1qwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "init=[12,1,2]\n", "fit=curve_fit(func,x,y,p0=init,sigma=sig*np.ones(len(x)))\n", "print('parameters:' ,fit[0])\n", "print('covariance matrix: \\n',fit[1])\n", "\n", "# plot the data and the best fit\n", "plt.errorbar(x,y,sig,fmt='ro')\n", "x=np.arange(-10,10,0.01)\n", "plt.plot(x,func(x,*fit[0]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate the $X^2$ statistic and it's probability" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Modify your function to fit the sum of two Gaussians, and\n", "see how well you can recover the parameters. Experiment with different starting guesses." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " DISCUSS YOUR RESULTS: " ] }, { "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 }