{ "cells": [ { "cell_type": "markdown", "id": "controversial-sally", "metadata": {}, "source": [ "# Variance of fit parameters\n", "\n", "We use the bootstrap and the jackknife to compute the uncertainties of a non-linear least-squares fit. The bootstrap is generally superior to the jackknife, which we will also see here. We use `scipy.optimize.curve_fit` to perform the fit, which also estimates the parameter uncertainties with asymptotic theory. For reference, we also doing a Monte-Carlo simulation of the experiment with a large number of tries, to have a reference for the parameter uncertainties.\n", "\n", "In this case, the asymptotic theory estimate is very accurate, while the bootstrap and the jackknife estimates are similar and off. The accuracy of the non-parametric methods improves with the sample size." ] }, { "cell_type": "code", "execution_count": 1, "id": "major-companion", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import curve_fit\n", "from resample import bootstrap, jackknife\n", "\n", "rng = np.random.default_rng(1)\n", "\n", "# generate some random data, each y value scatters randomly\n", "x = np.linspace(0, 1, 100)\n", "y = 1 + 10 * x ** 2\n", "ye = 0.5 + x\n", "y += rng.normal(0, ye)" ] }, { "cell_type": "code", "execution_count": 2, "id": "intermediate-currency", "metadata": {}, "outputs": [], "source": [ "def model(x, a, b, c):\n", " return a + b * x + c * x ** 2\n", "\n", "def fit(x, y, ye):\n", " return curve_fit(model, x, y, sigma=ye, absolute_sigma=True)\n", "\n", "# fit original data and compute covariance estimate from asymptotic theory\n", "par, cov = fit(x, y, ye)" ] }, { "cell_type": "code", "execution_count": 3, "id": "successful-inquiry", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABWAElEQVR4nO3deXwTdf4/8NckzdFCWyjQNmC5KgK1KIKACChIkct6rXIoCqwLirAu8P25gIqlqIDHuvVAUFYBRcQTlmurBQUWL5BSVyiCQkGUttwtUNpc8/ujJDRpjkkyk6uv5+PRXWYymflkqMw7n+P9FkRRFEFEREQUJKpQN4CIiIgaFgYfREREFFQMPoiIiCioGHwQERFRUDH4ICIioqBi8EFERERBxeCDiIiIgorBBxEREQVVTKgb4MxqteLYsWOIj4+HIAihbg4RERFJIIoizp07h5YtW0Kl8ty3EXbBx7Fjx5CWlhbqZhAREZEfjh49iiuuuMLjMWEXfMTHxwOobXxCQkKIW0NERERSVFZWIi0tzf4c9yTsgg/bUEtCQgKDDyIioggjZcoEJ5wSERFRUDH4ICIioqBi8EFERERBFXZzPqQQRRFmsxkWiyXUTYk4arUaMTExXMZMREQhE3HBh9FoRGlpKaqqqkLdlIgVFxcHg8EArVYb6qYQEVEDFFHBh9VqRUlJCdRqNVq2bAmtVstv8D4QRRFGoxEnTpxASUkJOnTo4DURDBERkdwiKvgwGo2wWq1IS0tDXFxcqJsTkWJjY6HRaHDkyBEYjUbo9fpQN4mIiBqYiPzay2/rgeH9IyKiUOJTiIiIiIKKwQcREREFFYOPIBFFERMnTkRSUhIEQUCTJk0wderUUDeLiIgo6CJqwmkky8/Px7Jly7Blyxa0b98eKpUKsbGx9tfbtm2LqVOnMiAhIqKox+AjSA4ePAiDwYAbb7wx1E0hIiIKqcgPPkQRMIUo4ZgmDpCQZ2TcuHFYvnw5gNpqf23atEHbtm3RtWtX5OXloX///jhy5AimTZuGadOmAagdpiEiIpJTldGMjKc/BwAUzx2MOG1owoDIDz5MVcC8lqG59hPHAG0jr4e98sorSE9Px1tvvYWdO3dCrVbj3nvvtb/+2Wef4dprr8XEiRMxYcIEJVtMRERhLlwCBCVF3ycKQ4mJiYiPj4darUZqamq915OSkqBWqxEfH+/ydSIiomgS+cGHJq62ByJU1yYiIiKfRH7wIQiShj6IiIgoPDDPR5jQarWwWCyhbgYREZHiGHyEibZt22Lbtm34448/cPLkyVA3h4iISDEMPsLE3LlzcfjwYaSnp6NFixahbg4REZFiIn/OR4Rwzl66ZcsWh9dvuOEG/Pjjj8FtFBERUQiw54OIiIiCisEHERERBRWDDyIiIgoqBh9EREQUVAw+iIiIKKgYfBAREVFQMfggIiKioGqwwUeV0Yy2Mzeg7cwNqDKaQ90cIiKiBqPBBh/hoH///g6Jx4iIiBqCBht8WKyi/c87Sk47bIejLVu2QBAEnD17NtRNISIiCkiDDD7y95Qi6+Wt9u1xS3ei7/NfIn9PaQhbRURE1DD4HHxs27YN2dnZaNmyJQRBwJo1a+yvmUwmzJgxA126dEGjRo3QsmVLPPjggzh27JicbQ5I/p5STFpRiPLKGof9ZRXVmLSiULEA5MKFC3jwwQfRuHFjGAwG/OMf/3B4/b333sP111+P+Ph4pKam4r777sPx48cBAIcPH8aAAQMAAE2bNoUgCBg3blzt58nPR9++fdGkSRM0a9YMt912Gw4ePKjIZyAiIpKDz8HHhQsXcO2112LhwoX1XquqqkJhYSFmz56NwsJCfPbZZ9i/fz9uv/12WRobKItVRO66YrgaYLHty11XrMgQzOOPP46tW7fi3//+N7744gts2bIFhYWF9tdNJhOeeeYZ/Pjjj1izZg0OHz5sDzDS0tLw6aefAgD279+P0tJSvPLKKwBq/z6mT5+OH374AZs3b4ZKpcJdd90Fq9Uq+2cgIqLwFUkLKXyuajt06FAMHTrU5WuJiYkoKChw2Pf666+jZ8+e+O2339C6dWv/WimTHSWnUVpR7fZ1EUBpRTV2lJxG7/Rmsl33/PnzePvtt7FixQoMHDgQALB8+XJcccUV9mP+/Oc/2//cvn17vPrqq+jRowfOnz+Pxo0bIykpCQCQnJyMJk2a2I/905/+5HCtd955By1atEBxcTEyMzNl+wxERERyUXzOR0VFBQRBcHhg1lVTU4PKykqHH6UcP+c+8PDnOKkOHjwIo9GIXr162fclJSWhY8eO9u1du3YhOzsbrVu3Rnx8PG6++WYAwG+//ebx3L/88gtGjx6N9u3bIyEhAW3btpX0PiIiolBRNPiorq7GjBkzMHr0aCQkJLg8Zv78+UhMTLT/pKWlKdae5Hi9rMfJ5cKFCxg8eDASEhLw/vvvY+fOnVi9ejUAwGg0enxvdnY2Tp8+jSVLluD777/H999/L+l9REREoaJY8GEymTBixAiIoohFixa5PW7WrFmoqKiw/xw9elSpJqFnuyQYEvUQ3LwuADAk6tGzXZKs101PT4dGo7EHBgBw5swZHDhwAADw888/49SpU1iwYAH69euHTp062Seb2mi1WgCAxWKx7zt16hT279+Pp556CgMHDkTnzp1x5swZWdtOREQkN0WCD1vgceTIERQUFLjt9QAAnU6HhIQEhx+lqFUCcrIzAKBeAGLbzsnOgFrlLjzxT+PGjfHQQw/h8ccfx5dffok9e/Zg3LhxUKlqb3/r1q2h1Wrx2muv4dChQ1i7di2eeeYZh3O0adMGgiBg/fr1OHHiBM6fP4+mTZuiWbNmeOutt/Drr7/iyy+/xPTp02VtOxERSRNJEz5DTfbgwxZ4/PLLL9i0aROaNZNv4qYchmQasGhMNyQn6Bz2pybqsWhMNwzJNChy3RdffBH9+vVDdnY2srKy0LdvX3Tv3h0A0KJFCyxbtgwff/wxMjIysGDBArz00ksO72/VqhVyc3Mxc+ZMpKSkYMqUKVCpVFi1ahV27dqFzMxMTJs2DS+++KIi7SciIpKLz6tdzp8/j19//dW+XVJSgqKiIiQlJcFgMOCee+5BYWEh1q9fD4vFgrKyMgC1EyxtQwehNiTTgD5XNkeXOV8AAJaN74F+HVrI3uNRV+PGjfHee+/hvffes+97/PHH7X8ePXo0Ro8e7fAeUXRc8jt79mzMnj3bYV9WVhaKi4s9vo+IiCic+Bx8/PDDD/aEVwDs3fxjx47FnDlzsHbtWgBA165dHd731VdfoX///v63VGZ1A42e7ZIUDTyIiIjoMp+Dj/79+3v8Zh0p37rjtDE4vGB4qJtBRETkwLn2mNI986HQIGu7EBERhSOla4+FS1FVBh9ERERhwJfaY/6srAmnoqoMPoiIiEJM6dpjoSqq6k5EBh+RMq8kXPH+ERGFF19qj/kqlEVV3Ymo4EOj0QCorZ5L/rPdP9v9JCKi0FKy9piSgY2/fF7tEkpqtRpNmjSxpx6Pi4uDIETXDGAliaKIqqoqHD9+HE2aNIFarQ51k4iICMrWHgtVUVVPIir4AIDU1FQAqFf7hKRr0qSJ/T4SEVHo2WqPlVVUuxweEVCbiduf2mPhWFQ14oIPQRBgMBiQnJwMk8kU6uZEHI1Gwx4PIqIwY6s9NmlFIQTAIQAJtPaYkoGNvyIu+LBRq9V8iBIRUdSw1R7LWbvXYVVKaqIeOdkZftceUzKw8VdETTglIiKKZkMyDdg0/Wb79rLxPbB9xi0BFz0NVVFVdyK254OIiCgaKVV7LBRFVd1hzwcREZEC/MlCqrRwKarK4IOIiIiCisEHERERBRWDDyIiIh+F45BKJGHwQUREREHF4IOIiIiCisEHERFRA6KCNdRNYPBBRETUYBgv4APtsxit3hzSZjDJGBERkRdVRjMynv4cAFA8d7Cs56m7HadV8LFsroHu07HopfoZnYTfgKoZgDZFuet5wJ4PIiKiKGCxXq7asqPktMM2LGbg04egLvkKVaIO441/B+KahaCVtRh8EBERRbj8PaXIenmrfXvc0p3o+/yXyN9TClitwNq/AvvWQVRrMcE0HYXiVSFsLYddiIiIIlpBcTmmripyqFYLAGUV1Zi0ohCLMvZiyKGVgKCG8c5/4ev3Q//oZ88HERFFpYaSCGzexn31Ag8Al/aJyC1OgUUUgDvfgKXj8OA2zg0GH0RERDLwOOdCQeWVNW5fEyGgFM2xo+crwLWjgtIeKULf90JERBQkzqtN5Fpdkr+nFDlr99q3xy3diZQEnSznlsPxVlmhboIDBh9EREQBcDfn4riHHolgS47Xh7oJDjjsQkREFADPcy5qKTkEk5Kgg+DmNQGAIVGPnu2SFLu+Pxh8EBERBcDTnAubXUfOKHb9J4Z1BgAITiGQLSDJyc6AWuUuPAkNBh9EREQKO3FOuSGYQRkpWDQwBinCaYf9qYl6LBrTDUMyDYpd21+c80FERKSwFvHKTT5VHf4vhnz/IAZpa/CK+W68ZrkLS8f3Qr8OLcKux8OGPR9EREQB8DTnwqZ7m6ayXtM2h6QtSrH7g1xYTDUQrxqMNyx3QIQKPdslhW3gATD4ICIiCsjlOReO6m7LGQjUTaV+GAaMvvg4+lrexIZOC2COkAENBh9EREQBGJSRgkVjuiHZKa9HSoL8y1vz95Ri0orCepNcy8yNMfXjYtmvpxQGH0RERAEakmnApuk327eXje+Bguk3yXoNi1VE7rpir8t6I0Fk9M8QERGFubpDK0rk1dhRchqlFdVuX4+kAMTnno9t27YhOzsbLVu2hCAIWLNmjcProiji6aefhsFgQGxsLLKysvDLL7/I1V4iIqIG6XjZ76Fugmx8Dj4uXLiAa6+9FgsXLnT5+gsvvIBXX30Vixcvxvfff49GjRph8ODBqK52H60RERGRB2ePInn70z69JZyr+vo87DJ06FAMHTrU5WuiKCIvLw9PPfUU7rjjDgDAu+++i5SUFKxZswajRoVPRT0iIiJXlCo+J1WcNgaHFwy/vKPiD2B5NnpWHYZBNRJl1kSXQywCImfoRdYJpyUlJSgrK0NW1uXqeYmJiejVqxe+/fZbOS9FREQU/SpLgeXZwJkSqJPaIOfOawF4XtYbCWQNPsrKygAAKSkpDvtTUlLsrzmrqalBZWWlww8REVGkq1tMbkfJad+Ly50rrw08Th8EmrQGxq7HkJ5Xu1zWm5qoR96orjK0OjhCvtpl/vz5yM3NDXUziIiIZFNQXI55G/fZt8ct3QlDoh4zh3aSdoLzJ2oDj1O/AIlpwNj1QJM0ALXLevtc2Rxd5nwBoHZZb78OLVBjtng9bb0hnRCRtecjNTUVAFBeXu6wv7y83P6as1mzZqGiosL+c/ToUTmbREREJDtvvRpTVxXVTwRWUY2pq4q8n/zCSeDd24GT+4GEVsDYtUDTNg6HOC/rDedU6q7IGny0a9cOqamp2Lx5s31fZWUlvv/+e/Tu3dvle3Q6HRISEhx+iIiIwlVBcbk9vTlQ26vR9/kvUVB8+Yu334nAqk4D794BHC8G4g3A2HVAUvuA2xxufB52OX/+PH799Vf7dklJCYqKipCUlITWrVtj6tSpePbZZ9GhQwe0a9cOs2fPRsuWLXHnnXfK2W4iIopSoV5t4s3UVUX1AgmpvRoeA5CLZ4AP7gLK9wCNU2oDj2bpAbQ0fPn8N/rDDz9gwIAB9u3p06cDAMaOHYtly5bh73//Oy5cuICJEyfi7Nmz6Nu3L/Lz86HXy5/jnoiIKNjc9WoEMvCRgAvQffAnoOwnoFGL2sCjeYcAzhjefA4++vfvD1F0H7sJgoC5c+di7ty5ATWMiIgokvibYyMeVXhXuwDqsoNAXLPawKNFR1nbFm5YWI6IiMgLX5fJSukF2VFyGpaqs3hPOw9dVQchxiYBD64Fkjv718gIEl4DaURERGEmf08pctbu9fl93jKOjlu6E6nqSsxRNcNp8Thi71uN2NRMv9sZSdjzQURE5EZBcTkmrSist2zWlbq9HXmjutZLBOZKmaUxHjFNxcCal/B9VUvfE5FFKAYfREQUlQLOMApg3sZ9kuZyOA+zDMpIwabpN9u3m8Zp3LxTBUDAGSTYl+zm7yn1uZ2RhsEHERFFnfw9pV5zcUghpccDcJ3evG7irzNVJknnKauoxqQVhVEfgDD4ICKiqJK/p9TlUElZRTX+VicXh7+9Ic4euak9ts+4BYMyUrwf7IWtNbnriqN6CIbBBxERRQ2LVUTuumJJGUblGua4Ib2ZrOnNRQClFdXYUXJatnOGGwYfREQUNXaUnEZpRbXk470Nc6Qk6Lwum+3epqnH15vjDFJVFRBgldwuADh+TvrniDQMPoiIKGr4+sD2NszxxLDanBvOAUjdbU+9Hsk4gw+1z2KO+h0Agk9ZUJPjozczOIMPIiKKGv48sD0NcwzKSMGiMd3qLZtNSfB+HeHcMazSPoN0VSlubVqKRXemSVp+KwAwJOrRs12S1I8QcRh8EBFR1OjZLgmGRL1fdVbc9ZoMyTQ4LJtdNr4HCqbf5PFcQuUf0K24He1VZThqbYGaMesw5IZrHc4DuO9RycnOkHUeSbhh8EFERGGvymhG25kb0HbmBlQZzW6PU6sE5GRnAPC90JunXpO6gUDPdkkeA4M0oRy6926D6kwJfrO2wCjjUxCbtK53nldcJCJLTdRj0ZhuGJJp8LH1kYXBBxERRZUhmQaXQyXuyDnM0V44ho+0z0BV8RusTdtjpPFp/IEWLo91TkS2bHwPbJ9xS9QHHgCDDyIiikLOQyVTBqRDgLLDHMLxYnyonQuDcBrW5h1R/cA6lKKZx/f40qPiTZw2BocXDMfhBcMRpw3v0m0MPoiIyGdSh0FCqe6D/NEBV7rsDZFtmOOPQujfvx0thErstbZB9f1rgcapgZ0zioV3aERERCSTIZkG9LmyObrM+QJA7TBHvw4tAu7xUP2+A/hwJISaSqDV9bh6zCdAbNOwDcrCAXs+iIiowXA1zBFIL05v1V7oPrgHqKkE2vQBHlwDxHpOOkbs+SAiIvJLf1URFmv+CcFkAtoPAEatBLRxoW5WRGDPBxERkY/U+9fjLc0/oBdMMHcYAoxexcDDBww+iIiI6qibZt1V5Vv13k+h/ezP0AoWrLf0gvHuZYAmelOhK4HBBxER0SUFxeXIenmrfdtW+baguBwAcK96C7T/fhiCaMGnln74m2kKoNaEqLWeeQuiQonBBxER0SVTVxWhvLLGYV9ZRTWmrirCzaoivKh5CwJEmK4bh/9nehgWqEPUUs/y95S6DKLcVe8NNgYfREREl7jqGxAv/e8B6xWwiAJMPSfBNOQliGH6CM3fU4pJKwpdBlGTVhSGRQASnneOiIgarHAcLhAhoBTNMcM0AaaBzwBCeBZ9s1hF5K4r9hBEAbnrikN+T7nUloiIwkb+nlLkrN1r3x63dCcMiXrMHNophK267BNrf8wNg8DDlkrd2Y6S0yitcF2dF6gNQEorqrGj5DR6p3tO/a4kBh9ERFGsymhGxtOfAwCK5w4O65oftuEC5+/ktjkX5N3xc+4DD3+OUwqHXYiIKOSkDBcEg7s+jdD3dUiTHC9tya/U45TC4IOIiByEomiclOGCYBGcrhYpgQdQmzLekKj3GEQZEvXo2S4pmM2qh8EHERGFXKiHAQCgBc5iYbOPkYrTDvtTE/XIG9U1NI3ykVolICc7A0D9oMm2nZOdEXAxvUAx+CAiopAL9TDAFcJxfKzNxbALq/HfpOeQhtqkYsvG98D2GbdgUEZKSNvniyGZBiwa0w3JCTqH/amJeiwa0w1DMg0hatllDD6IiCjkpAwXKEU48TM+0eairaoc1iZtYBq7HkeRYm9XqHsJ/DEk04BN02+2b9uCqHAIPAAGH0REFAakDBco4o9d0K/IRqpwBgesrVDzwAaITdspecWgqRs0hVsQxeCDiIjCgqfhAiXmXKhKtgDLb4dw8TSKrOkYYXwaYrx8PQO2XByHFwwPeIlzKCYBKyl8F3wTEVHU8ZZ3ZEimAX2ubI4uc74AUDtc0K9DC9SYLbK24zbVt9B9uBiwmmBp0w/37x+HC4iV/H53Sb5IGvZ8EBFRWFF6uGCs+nO8qnkdgtUEZNyJmpEf+hR4UOAYfBARUcMgitBseRa5muVQCSJM3R8C7nkHiNF5fy/JisEHERFFP6sZWDsFmm/+CQB40TQCplufB1TqEDesYZI9+LBYLJg9ezbatWuH2NhYpKen45lnnoEohr4qIRER+SfUEx4DqXSrRw20n44Fdq+AKKgwwzQBCy13hm1l2oZA9gmnzz//PBYtWoTly5fj6quvxg8//IDx48cjMTERjz32mNyXIyKiKOeu0m1OdobXvBWJOI9/aV9CzC8HgBg9jHcuwYcrfO/t4ARTeckefHzzzTe44447MHx47V9S27Zt8cEHH2DHjh1yX4qIiELEuSeiX4cWfk8MdV4BU1dBcTmmripyWel20opCjxk7hco/8JF2LjqqfoeoT4Qw+kNYDD0AfO5XO+XgHMREw7JZf8g+7HLjjTdi8+bNOHDgAADgxx9/xPbt2zF06FCXx9fU1KCystLhh4iIwldBcTmyXt5q3x63dCf6Pv8l8veUyn6teRv3eax0m7uu2PUQzIn9iH13GDqqfgfiDRDG5wNtesvePvKP7D0fM2fORGVlJTp16gS1Wg2LxYLnnnsO999/v8vj58+fj9zcXLmbQURECvG3J8If5ZU1bl8TAZRWVGNHyWn0Tm92+YWjO4GV9wIXzwDNOgAPfAY0aS1bm+TUUIdzZO/5+Oijj/D+++9j5cqVKCwsxPLly/HSSy9h+fLlLo+fNWsWKioq7D9Hjx6Vu0lERCQjv3oiFORQEffnDcDy7NrAo9X1wJ8/D9vAoyGTvefj8ccfx8yZMzFq1CgAQJcuXXDkyBHMnz8fY8eOrXe8TqeDTsc11kREkc5tT4TC7BVxdywB/vN3QLQCHW4F7l0GaBs5HNtQexrCjew9H1VVVVCpHE+rVqthtVrlvhQREYUhh56IAKUk6DxWujUk6tGzTROg4Glg4/+rDTy6jQVGfVAv8HAnkGW85B/Zez6ys7Px3HPPoXXr1rj66quxe/duvPzyy/jzn/8s96WIiCgM2XsiZPDEsM6YuqoIAhyHe2wBSc6wq6BeMwHY82ntjltmA/3+T3IOD1fLeFMS2BuvNNmDj9deew2zZ8/Go48+iuPHj6Nly5Z4+OGH8fTTT8t9KSIiCgHnQKDu/tREPXq2S5LtWoMyUrBoTDfkrN3rMPk0NVGPnMFtMKTwYeDI14AqBrhjIXDtKMnnzt9TikkrCut9luMeJrmSPGQPPuLj45GXl4e8vDy5T01ERGHCbU9EdobsheBcVrptcRHqlfcCJ/cDugRg5HuouqIvMmZuAOC6Ym5dFquI3HXFHifP2o4j+bG2CxFRFFNiPkPeqK5IdhqaSE3UI29UVzyyolCRFOx1A5pesUehfmdQbeAR3xL4cz7Qvr9P59tRchqlFd7npuw6csbXppIEsvd8EBGR/JyzgHr6Vm8TSFpyTwZlpOCWTsmOPREdWqDGbPH7nFLdrPoR+vcmAKYLQEomcN9HQGIrn88jdVLsiXMcglECez6IiKKQbT6Dc5IuWzKwQLOR1u2J6NkuSfahFldGqL/C25oXIZguAO1uBsZv9CvwAKRPim0Rz8mnSmDwQUQUZaTMZ/CUDCzslp6KIjRb5+MFzRLECFaYu4wE7v8E0Cf6fcqe7ZJgSNS7XcZr071NU7+vQe4x+CAiijLe5jPUTQbmLH9Pqcu6LQXF5Uo01TtzNfDpX6D5+iUAwCvmu2C8bSEQow3otGqVgJzsDACoF4AITseR/Bh8EBFFGanzGZyP8zRUM3VVkVzNkywJldC9fyew5xOIqhg8bpqIf5rvlZzDw5shmQYsGtOt3uTZlAT58pSQaww+iIiijNT5DHWPk7r0NFjShT+wRjsb6j92AvpE1Iz6BB9b+st+nSGZBmyafrN9e9n4HiiYfpPs1yFHDD6IiKKMt/kM9rTkdZKBSRmqCRZVyVas1uagteoErE3bAX/ZDGvbfopdLxSTZxs6Bh9ERFGiymhG25kbkP7ERswc2gmA+/kMzsnA5KzHEgj17neh+3AEEoQq7LReheqxnwPNO4S6WSQzBh9ERFHIlpbcVTKwRWO61cvzIWc9Fn8IsGJmzEro/jMNgtWM1ZY+uN/4JBAXvOq4NrbKt4cXDJeUT4V8x7tKRBSlXKYl79DC5bCCbaimrKLabd0WOYZenJfx9uvQAjBVYZHmFQxR7wQAGPvNwLSCa1C/34aiBYMPIqIoJnU+g23p6aQVhW7rtgTKZcbVeA1mx36EYeqdqBFjgDteh67baBweKNNFKSxx2IWIiAC4X3pqq9sSCLfLeM8ZMfn4nfjYfBPuMz4JS+a9AV2HLgvn4aPwag0REYWUu6GaQOq2eF7GK0CAiJnmCbBA7fc1KLKw54OIiBy4GqoJJOW692W8AgMPL8Iu5X2A2PNBREQeuZqrkZIgveBauCzjjVRKVScOJfZ8EBGRWwXF5S7nahyvlF5qXollvLacJm1nbkCV0Sz7+cOF0tWJQ4XBBxFRA+LrQ3vexn1eU657GwLoqTkEg+oMBFi9Xi8ahhTkEmh14nDG4IOIiNxy/sbtyq4jZ9y/+L+PoF4+HDnqZQAEr8t2bVV0I/UbvZwCqU4c7hh8EBGFmUgbUjhxzkWAIlqBTXOAzyYAlhoM6dwci0ZeXW8ZryuRPqQgF3+rE0cCBh9ERBSQFvGOAUUjXIT2kweA7f+s3dF3OjBqJYZc186hgmzTOI3L80X6kIJc/KlOHCkYfBARkVspCTqvQyXd2zS1/zlNKMen2jmI+SUfUOuAu5cAWTmAqvZxU3cZ75kqk9tzRvKQglz8qU4cKRh8EBGRW08M6wzAfXVc4HJAoTr0FdZpn0In1VFYG6cA4/8DXDMioOtH4pCCXGwp7wHp1YkjBYMPIiJyy1113JSEOl39oghsz4PuwxFoIlxAkTUdNeM3AVd0D/j6kTikICdPKe9dVSeOFAw+iIgiQCgzXA7JNDjM1Vg2vgcKpt8EAIhFNbRr/gJsyoEgWvGhuT9GGmdDjG/p9byehnQieUhBbq7u//YZt0Rs4AEw+CAiCnv5e0qR9fJW+3YolqO6SrmeJpTjM20OYvatAVQxuHjri5hhnoAaaCUFSN6GdHwdUoi2FOR1Sa1OHCkYfBARhTF3GUZDvRw17ret+G9iLjqrjgKNkpHf7zP039IettBBSoDkbkjH1ZCCtwqt4RCgkXQMPoiIgkjuDKPBX44qIubbV4H37wGqzwKtrkf+zf/GpM/P+xUgyTGkEK4BGrnH4IOIKIx5yjAq13JUqQFRLKrxuuY1aL/KrU0idt0DsIzdgNzNpQEFSIEOKYRfgEbeMPggIopwwViOKpw5jM+0ObhN/R1ElQYY/jJw+2vYcfRCyFOAByNAI3kx+CAiinCKL0f9pQD6pQPRWXUUJ8RE1Nz/b6DHQ4AgSA58Ri/5LqSp4htyvpBwxOCDiCiMhXI5qgpWaLYtAN6/F0L1Wey2Xonbap6DNa2X/ZhIycMRKe1sKBh8EBGFMV+Wo8q51LQpKrFM8zw0218EIMLU7c8YaZyNcjgGOlJSgCuN+UIiD4MPIqIwJnU5qrulpgXF5T5fU/XHD1ivexI3qX+CGBML3PUWTENehBH1C8FJSQGutEDzhXhbxkvyY/BBRBTmvC1Hzd9T6nap6dRVRVg8ppvEB6uImF1vQ/febWglnMIhayqqx30BXDvSa/vcBUh5o7pK/pz+8iVfCIUHhnhEFNWqjGZkPP05AKB47mBZv9kqeW5n7pajWqwictcVu11qKqB2qemgjFSP3/5jUY15mreh/fxrAMB/LD3wuOlh7EjOkNS+IZkG9LmyObrM+QJAbYDUr0ML1Jgtkt4fKHfXj/RMoNFKkZ6PP/74A2PGjEGzZs0QGxuLLl264IcfflDiUkREDdqOktMBL3UVTv2CNdqncZf6a4iCGsaBczHJNBXnEedTW0KdAjzU1yfpZA/Tz5w5gz59+mDAgAH4z3/+gxYtWuCXX35B06ZN5b4UEVGDJ3UJqdvjiv8N/ZpH0VF1HsfFJkgYswL6K/vhcD8ZG0nkRPbg4/nnn0daWhqWLl1q39euXTu5L0NERJC+hNT5uBiYodk0G9jxBgQA31s7YYrxMWxt3VuBVhI5kn3YZe3atbj++utx7733Ijk5Gddddx2WLFni9viamhpUVlY6/BARkTRSlro6LzVthRP4SDsXmh1vAABMN/wV9xmfxAk0Uby9RIACwcehQ4ewaNEidOjQAZ9//jkmTZqExx57DMuXL3d5/Pz585GYmGj/SUtLk7tJRERRS8pS17pLTdUHNmKD7gl0U/0KUZcAjFwB0y1zYIE6aG0mkj34sFqt6NatG+bNm4frrrsOEydOxIQJE7B48WKXx8+aNQsVFRX2n6NHj8rdJCKiqOZpqat9qanZCOTPQszHD2KftTVeNd2JbbdugKXjbSFqNTVkss/5MBgMyMhwXJrVuXNnfPrppy6P1+l00Ol0Ll8jIvIkmEtdleT8OZwzlfbr0MLrOTwuNT1zGPh4PPKPqpFrehWlaFb7po+PwvDFCcwc2kn2z0Tkiew9H3369MH+/fsd9h04cABt2rSR+1JERFGnoLjc70ylLpeaFq8FFt+E/KNqTDJNRalTenRbIrJII2cqeQo+2b8mTJs2DTfeeCPmzZuHESNGYMeOHXjrrbfw1ltvyX0pIqKgCVYvy9RVRfUShvkVIJhrgE1PADvehEUUkGv9C0QX01JticgiSf6eUuSs3WvfHrd0JwyJevbgRBDZ/+vp0aMHVq9ejVmzZmHu3Llo164d8vLycP/998t9KSKiqOMpU6lUrYVy6N4dCpT9CADYkfEkSnfH+3TNcGVLJS9LgEYho0joftttt+G22ziJiYhILlIDhOGq7zBfswTqsotAbBJw12Icr8oEdhcp2TzJbEXc/CEllTxFhsicnUVERI5qzkO74XEs1K4EAFiu6AX1vUuBxFZIPngqxI2Th5RU8hQZWNWWiCjSHdsNvHkTYv63EhZRwCvmu1AzZi2Q2AqAtERkkUBqKnkKf+z5ICIKIwJcf4N3ud9qBb59Hdg8F7CaYI1vidEnH8IOsTMmqC7/825LRDZpRWG980gNPAIZLpGL1FTyFP7Y80FEFERSloi6y1Tq4FwZsOJuoGA2YDUBnbNR/Zdt2CF2dnldT4nI8kZ19e1DhEi09OAQgw8ioqDJ31PqNYdH3qiuXgME1a9fAIv6AIe+AmJigexXgBHvAbGeq4cPyTRg0/Sb7dvLxvfA9hm3YFBGSoCfLDikpJKnyMDgg4jCUpXRjLYzN6DtzA2oMppD3ZyA2ZaIllfWOOx3XiI6KCPFbYCggxE5Mcuh/2g0UHUSSOkCPLwV6D4OEKQ9fl0lIoukhF3R0INDnPNBRKQ4X5eIugoQhBM/Y432aXRW/Vb7wg2PAgNzAE1g8yDcJezKyc6orQkTBL7OJ3GXSr7GbFGqiSEXDnNu5MSeDyIihQW0RFS0At8thn7pQHRW/YaTYgKqR6wChswPOPAoKC532xszaUUh8veUBnR+JblMJU8Rg8EHEUU1V0MKwR7S8XeJaCpOQffBPUD+DAjmamy1XIOhNQtgvXKQLO2at3Gf294YAMhdVxzWQzAUuRh8EJGiQjl3Q8oEz2DwZ4nobapv8bluBtSHtwIxsTAOfgFjTTNwAk1ka5dzj0ddIoDSimrsKDkd0DUiaT4JBQ/nfBBRVAq0BohzITlXnB+s9hL2TmxLRMsqqr3n8Lh4FtovZuB17Se11zBcB/WflkDbvAMO9/babNkFktgrHOaTUHhizwcRRR1vEzzl4K5XxdU8CalLRG9U7YH+X/0Qs/cTmEUV8sx3o+bB/wDNO8jUat/5m9grkueTkPIYfBBR1FG6Bog/D1aPS0TvycDsmPewUjsPqnPHYG3aHvcY5yDPfA+g1gTYWvdSEnQeE3YZEvXo2S7Jr3NzPgl5wuCDiKKO0jVA/H2wukzy9WAzZH83Eg/F/AcAYLpuHKof2oIi8Ur5G+7kiWG12VDd9cbkZGf4tIrEthz0gwk3BGU+CUUuBh9EFHWUrgESyIPV9jCPgRk3/vEO1G9nQXVyP06IiRhvfBymof8AtI2UaHY9gzJS3PbGLBrTze95GVKDPxaKa7g44ZSIIpa7CZ8+TfBUiKcHa0fhN7ykWQzttsMAAPNVwzH4f9k4jQSFW1Wfu4RdgeTNkBr8yR0kRlsirmjGng8iikieJnyGQw0Qlw9WixkxX7+Mddon0UV1GKK+CXD3v2D80/KQBB42cifsklIALpD5JBT5GHwQUcTxVCfFNuFTyRogfk3UPP4z8HYWtFufg1awoMDSDRcnfg1cc6/kuiyRQkrw5+t8EoouDD6IKKJIWUZrm/CpVBVXnyZqWszA9n8Cb/YDju2GqE/ENOMkTDD9H9A41f6ZbHxJxOWcwM027HB4wXDEaUM7qu4p+AtkPglFBwYfRCQrpTOaSllGW3fCp1xDCnUDgsRYDRbeV//BmpKghwjgkRWFtZ/9xH5Y3h4EbJoDWIywXHkrqid8jdXWfrCFKq6Gj+puRzJ3wR8DD2LwQURBF0iAEoqVFAXF5fUChGc2FGPmkE72fcvG90DB9JsAACpYEfPda8DiflAfK0SlGIf/Z3oYNfeuhBhvcDivq+Gj4x5W00QaFoAjV7jahYgiSihWUkxdVeQyTfv0j360b9vmeFwlHMXzmiXQfvkrAMDSfiBuLb4TZWiGuU5zO7zlCwHARFwUldjzQUQRJRQrKSSlaTfXQLNtAdZrn8B1ql8h6uKB219HzcgPUYZmLs/rKV+Iza4jZ3xuL1G4Y/BBRBElnFZS2AKQDjgK/TsDoNn+4qWVLN1RPfEboNsDAa9kOXEueoZgiGwYfBBRWPK0AiTcVlJMjvk3VCf3Q4xrgcnGxzDBNB1ifEtZzt0iXuf9IKIIw+CDiMKOlIqx4bSSIkU4A3OXkbg48RtssN4AqanMPOULsenepmnA7SMKNww+iAKg9LLShkhKAjGbYK2kcD+/xIpknMEb5tthzH4DiPNtnom3fCEAuDqEohKDDyIKG74kEAu2+gGCFYCAs2iE/4rX+nVOd4XdUhKULYxHFGpcaktEinJV/M0dXxKI9U53vYJECXm3GTA//xeUmRvb96U2jsHM267B31YVBXRuJQq7ecLiaxQOGHwQkWIKissxb+M++/a4pTthSNRj5tBOLo8PVgIxd9VwnWlgxkT1ety+ZS1uU9fga+FqrDQPxKgHJqJfp1aoMVsCaocNE3FRQ8Pgg4gU4y4511Q3vQXBSCCWv6cUOWv32rdtAVFOdobDZFXVka+xUTsLHVR/AGYA7W5Czv47UCIa8PKVqQwQiALAOR9El3DyqPwkJeeqQ+kEYu7SmTtMZr1wClg9CZoVd+CkmIAV5luw7YYlqBr5KUpE1iQhkgN7Pogo6NwFILYEYpNWFEJwOs62XVpRjRqzxa+qrZ7SmQsAcj/diUG6v6HgQjpyTa+i1JaZdAuQUrjN5+uRcjh3JbKx54OIwoq7BGJyrADxlM5cBFB6UY3Xz92MSaaplwOPS6Kp2BtRqDH4IKKw4yqBmK1irNKWqu6G6GLgh8XeiOTDYRciUozz0Im3/XU5rwAJlrMmtddjdh05g/4dk/06P4cLiNjzQUQK85S9M9g8pTMXADSJ1Ug6D4u9EQVG8eBjwYIFEAQBU6dOVfpSRKQgqauB6g5JTB6Q7rL4W96orko1060EXMDslj8AEC9lJ73M1hNz9qJJ0rlY7I0oMIoOu+zcuRNvvvkmrrnmGiUvQ0QKqDKakfH05wCA4rmDJb3HOYfG618dRHK81r5ty94pV3IuKeJiVDg84jSwaQ5w+CRiND0wBw+jzBRnPyYlQY+yyupLf9bheGWNx2EhFnsjCoxiPR/nz5/H/fffjyVLlqBpU/6HShTt3BWEO3HOaP+zLXunc4ZRxSZw/lEIvJ0FrJ0CVJ0Eml+FIeOeRMGTd9oPcZ7M2hCLvdnmoRxeMNyvJcxEvlIs+Jg8eTKGDx+OrKwsj8fV1NSgsrLS4YcoUgTtIRrmpBSEsx2Xv6cUWS9vte8bt3Qn+j7/pUO12oCdPw6s/Suw5Bbgj12AtjFw67PAI18D6QM8pjMPpNgbfx+IpFEkxF21ahUKCwuxc+dOr8fOnz8fubm5SjSDSFFS03Q3BN4Kwtm8ufUgFn510GXK9UkrCrFoTLeA7l1V1Xm89tz/4dGYfyNeuFi785qRQFYukCD9vK6KvXVv09S+bVN35Yq734e/D+5o3+epjgxRQyJ7z8fRo0fxt7/9De+//z70eu/fFGbNmoWKigr7z9GjR+VuEpHs3A0xOKTpbkCkFnp777vfPPaO5K4r9rO3QIR631ro3+yNGZpViBcuwmLoCozPB+5+y6fAw8aXYm/ufh9KK6ox7aMf7dty9/Kwp4UilezBx65du3D8+HF069YNMTExiImJwdatW/Hqq68iJiYGFovjRDOdToeEhASHH6JwJmWIwf+HaGSSWuitwsNqElvq9B0lp326dqZwCB9qn4Fu9XioKn5DmdgU042PoGZcAdCmt0/n8oen3wdX5ApQgzJ8RaQQ2YddBg4ciJ9++slh3/jx49GpUyfMmDEDarX3BD5E4czbEEPdh2jv9GZuj4smtoJwZRXVkh/C7kjtRUFlKbSbcrFWuwoqQYQYEwvzDX/FgE2dcBF6PCsEJ42R1CEnG3sdmXXFGJThX3VcW0+LUsNXcmJSNXJF9v864+PjkZmZ6fDTqFEjNGvWDJmZmXJfjkg2UruwpT4cJT9Ew0CgFX1tBeGAwJOKJXpJ9KWDETHbXwJe646Y/30AlSACXUZA+OsPMN00AxcReA0YX/jz9+xvLw/AnjeKDsxwSgTfurClDjFIPS5aSCkI5ynDqI3bHBpWC+5WbcOXuv+Ddtt8wHQBuKIn8JfNwJ+WAIlXBPgJ/BPI37M/gYsvPW9E4SoowceWLVuQl5cXjEsR+czXyaO2IQZPaboNifqg1iMJFeceE28F4fzLoSFCdXAT9O8MwMvaxWglnII14QrgT28DD30BXHG9vB/KR95+HzzxJ3CJxp43anjY80ENmj9d2FKGGHKyMxrscko5c2h0EQ5hpeY56D8cCdXxvagU4zDfNBrVD38HdLkHEEJ/jz39PrgTSIDKnjeKBgw+qEHztwvb3RBDaqK+3mS/QOdTBCrU13fmrXcEAHC6BNo1E7BO9xRuVBdDVGth6jUZ/Wry8KYlG9DEBrnVnrn7fXAl0ACVPW8UDZhHlxq0QLqwXSWiYgIpaZx7R2ySUAnNF7OAwqWIsZpgFQWstvbBsMmvQ0xMQ8XWz0PRXADeV224+n2oqrEgd/1ehyG9VBeJ6HxZEWLraZm0otBeEM+GPW8UKRh8UIMWaBe2qyEGfwqyNXjG85isXoNHYtZB80NtZlJLuwHI/nkQisW2GJqYVm81UjgGes6/D3HaGPS7Sv4A1dbTkrPWe2BDFI4YfFCD5i0/hYDaf9DDpQvbObCJ9CJgOhgRs2MxNN/k4XHNCQCANeUaqG6di5q0fii+9FkLissxb+M++/siKZW9L5lSfcGeN4pknPNBDVo4Tx4Nt7kasrKYMFq9GV/ppkO76UkIVSdw2JqCx4yTUf3nzUD6APtQxOIx3TB1VRFT2bugVGBDpDQGH9Tg+TJ5lAJktQA/fgj9W70xX/M2WgqnYY1viZqh/0SW8UWstfYB6mQmZUItoujE4IMIrldgbJ9xCwOPOgIrYiZisGoH9P/qB6yeCNWZEpwQE5BregDVk3bCct2DMLsYBWZCLaLoFNkDxkQyYhe2o7rBxRtf/YqPd/1u37bNuZg5tJPnk4giVIc2Y632KVyjKgFOAtAnwnjDX3FTfntchB6Px7if9MuEWkTRicEHEdWTv6cUOWv32rdf/+pgvWPKKqoxdVWRmzOIGKAqgm75S1AfK8Q1KgDaxsANk4DeU2BWN8bFfO/LZplQiyg6MfggIgcFxeWYuqrIa3VaW3VWx50i1Af+g7Xa2bU9HccAxMQCPR4C+k4DGjWvPU7iBNpgrUZytYyXiJTD4IOIHMzbuM9r4GFz+Tgr1Ps3AF+/BF3Z/3CNCqgSddDcMAGafn8DGif71RalEmrVTerl3MsjeUiJiPzGCadE5MB5SasUM9QfQvfpg0DZ/yBqGmGRORt9a16BaWCu34GHjZKrkTwVFXQ/pEREgWLwQSER1TksGqCuql8hahsD/f4PFyfvxvPm0TiNBNnOL3U1ki8rcqQs4yUiZXDYhYj8zpMhwIoUnMF31k64dvI6xCW2kDyfw1feViO5Gz5xlwVVyjJeAPhgwg0Rn0mWKNyw54OogcvfU4qsl7f6/D7h0pTTCsThFcu9QGxT+RsnkafhE3dZULmMlyh0GHxQ2IjWoZjAknMpy91D2zXHdqcm6JE36jpcRODl7QO5R/5mQeUyXqLQYV8iUQC8lUIPRUE0qVV1PT20azkupk0VKpCEM9gvtsbb429Avw4tHIJEf5eoBnqPfMmC2ju9mX1/pBUVJIom7Pkgklndb9h/C4OCaK56FaqMZqQ/sdHjQ9sWeNyr+grvDjShYPafUCy2gwVq9GyXhILiMofhmnFLd6Lv819i24ETOLxgOA4vGC5prkSgReP8HT4J56KCRNGOwQeRjKTMnwhmQbSC4nKXAUJBcbnkc3xs7Y/rb74N6pjLgURBcbnPcyzcCbRoXCDDJywqSBQaDD6IZOLL/Al/C6L5OjfCXa/C1FVFUEPqvJr63/zdJSKTGjTEaWPwwYQbPF5V6j2yDZ+4658QABg8DJ+wqCBR8DH4oAbJ1eRW2/wNqcMFdXmfP+GapyED50Bj4/9KXfZieOplcB8giGiGSqTiFARYXb7X02CDpwBLatAg12oTOYZPWFSQKLgYfBDJwNukR3fcDRk4D9+MW7oTj66UZ5gDAEQIOI4kDFN9D0Bw+9AOhLegQc7VJhw+IYosDD4o6gVjCa+vuSBsD/fRS76r1yZfh28A/+ePLLUOQd6o61w+tPNGdfX5fHV5CxoCHS5xxuETosjB4INIBr7kgvDUq+DP8I2rYY5EnJf4XhUGZaS4fGgPykhx+76UBF3AQYMSq01CPXwSyNAdUUPC4INCIpwTb/nD27f4ujz1Kvg7fAMAxyurgd++h3bNRHynfRQGH+Zz+PrQfmJYZ5fn8TVoaAjDJQxIiOpj8EFB52o+g6/LP8ONp2/xdXnrVQgklXfylr8D79yKmOJPEasy4/6YTVBqPsegjBTZggYOlxA1PAw+KKgivYS5p/kjbr/FJ1wekvHWq+BPKm8BIgw4iZ5nNwAxepivvR/Da57DS+aRis3nAOQNGkI9XBKp2KtCkYrBBwVNQyhh7uqBXDD9Jsnv92X4BoB9WCWnyX+gvjUXmL4PxuGvYq/YDgD8ms/hCwYNROQPBh8UNFJLmNuOjdR5IIE8kD0P34ioV9xNcxGLBmox5O8rgT5/A+LqT/JkgEBE4YbBBwWNL/MZpCTQilb24Zt4rcN+A07hjZhX8M+Y1zFQ2IVlo67EF0/ehUc2m9D2if9EVSVgIopuHCSkoPF1PoMtgVa0rHqQxGoFDm/DkAMrkSWux05NaxxHEySrzqPn1R1gunYmOi+tgggVXstID3VriYj8wuCDgsZbCXNntoLuueuKMSgjNbqHC04dBIpWAj+uAip/B1D7H2eycBZfWq7DtOk5UDdNRY3RDBGfh7atREQBYvBBQWObzzBpRSEESJtkWjeBVu/0ZvVerzKakfF07cO4eO7giJrxH48qDFd/B927ecDvOy6/oEsEuvwJ1VePwMA3TwEQMK1Rc8nndc6h0r1NU/kaHUK2lR1EFPki519qklWoHtq2+Qw5a/dKSh9uE0j+i7BiMQEHvoS2aBV26tZCL5iA3wEIKiB9IND1PqDjMECjh9VoBnzs5cjfU4qctXvt2+OW7kSK01JbIqJQY/BBQTck04A+VzZHlzlfSH6PP/kvPHHuHejXoYVywzqiFT2FfbhD/Q1iX50CXDyNGAAxArDfegXaZf0F2utGA/GpAV2moLgcU1cV1etROu4iyPOnd4Q9D0QkF9mDj/nz5+Ozzz7Dzz//jNjYWNx44414/vnn0bFjR7kvRRGs7oM+JUGH45U1LodhBNQmxZJaXAyo36vj/KCtqrEgd71j74AhUY+c7AwZJ7aKyBRKoNn8NNTFq/GR7ljt7osAGiXD1PkO3P11G/wktkPxDUOglaHnad7GfV5zqFisouTeEQYbRKQU2YOPrVu3YvLkyejRowfMZjOeeOIJ3HrrrSguLkajRo3kvhxFgSeGdcbUVUX15oH4W1ysroLicszbuM++PW7pTpfHybay5uQv0BR9hM3a95CuKgW+r91dKcYh39IDtz/wGPRX9ofJAvy0Xd6Jo1KGsd7cehALvzooqXeEiEgpsgcf+fn5DtvLli1DcnIydu3ahZtukp7pkSKXr/NJbHVCnOeBpMrQG+FqGMIVv1fWiCJw4meg+N9A8Vrg+F5oAKSrgGpRg5jOw2C5+k/o8b4VNdDitnb9AXUMLCaT/RSKD/vU8d53v0nqHSEiUpLicz4qKioAAElJrrvNa2pqUFNz+YFTWVmpdJMCEurVFaG+vlKc54EsG99Dlgeyv6XpXa2sqXvk1cJhaLY8C+xfD5z65fJLghqW9gPwf/uuQoG1O3befScAoKbOxFFXwx7yD/u4VnHR5PWYXUfOoH/HZEXbQUQNm6JPLqvViqlTp6JPnz7IzMx0ecz8+fORm5urZDMoQoRLGnCXK2usVuCPH6DZswb/1X6ENNUJ4JtLr6m1QPsBQMbtQMdhqIlJwJqnXQ+puJsUKsewj6e5M744cY5DMESkLEWDj8mTJ2PPnj3Yvn2722NmzZqF6dOn27crKyuRlpamZLNCIlp7LKRy/vzhzL6yxnQRKNkGHMgH9ucD545BAyBNBVwUtdB0GoyYzDuBDrcC+oTLJ3BKc26buGmxiuj7/Jduhz0CTajmae6MLwFJi/iGtzSXk2uJgkuxJ+CUKVOwfv16bNu2DVdccYXb43Q6HXS64P9j19CDAVca+j0RAKQmaNDzzHrgg3zg0BbAVHX5AG08aq4cjDG7O6NIvBJLut3o09CQlMJ60oZ9XHM3dyYlQY+yyupLf/a+sqhfhxY+X5uIyBeyP11EUcRf//pXrF69Glu2bEG7du3kvgSFkXAPWKR+6xcuHZVz8QWo119eEXNMTMJmSzfce99fsMWUiZwN+1Eu1j7Y3c3VcPctWmqitEASqrmaO9O9TVP7tpIri8IBezCIIoPsVW0nT56MFStWYOXKlYiPj0dZWRnKyspw8eJFuS8V0aqMZrSduQFtZ26wVyN1tU8pzrkvonmFg+tHqVNpepzCIk0ehqh3Aq26AwOewsWHtuDGmtcw2/xnfG68BpM++F+95ay2uRpSqu9KTZSWGKuRdJw7nubO2HpHkp3yeqQm6htWAT8iCinZv6YuWrQIANC/f3+H/UuXLsW4cePkvlxYCPdv/87crbaYObRTCFuljLxRXTFv4z6HoMGAU5gd8x6aCudqK8bGVKNnh1ZQd3oQ6PA+EJ8CABCNZgC1ycE8JfCSOldDamE9uWuxOPcGyLmySK6eBvZYEDUsigy7UPjK31OKSSsKXa62mLqqyGFfJARVrtKkQxSRJpSjn2oPhhe/j+H4L3ZqWtUGGjiLnqqfIbToiHfK2mGb9RosnvUY1HGNPV7HUwIvqXM1PBXWq7vtKQhwfkj720MWLiuLiKhhCr+nSZhz+bBTiL8Pf3fvs1hF5K4r9vgNPljk+KbrsgdHexGz9R/hv7pLS10P1P5fJ5UVJ62J6Hb7JKivykJVbAqetS2HjZGnboyUuRruCuvVnRRqI+UeBfP3kYhILg0m+JCy1NPbwz4chisCKYgmZbVFpMjfsQeTPjt8qc2XP3+ZUYfJxgfxWsxZJKvOomv/u2FtPwDdFx+HFSoUXzsY0MbUWw4rB6lzOrxNCpUqHH4fiYj80WCCD2e+fmP0ZbhCKYFmxgxWWXrZv41brcDJ/cBv3wFHv4flyPfILZsGEUlw7q8RL82hnmJ+DICAZa16oHtqU1jh24PdmdzF75yHPXwVDr+PRET+apDBh6tiY56+McoxXOFPCXPnNgeaGVPusvSuyPFtXI8aqI5sB0p3Ar99D/y+A6iusL++w9IZpfCWB0OwX99VxVZfhdMS1XAaPiMi8keDDD7cPcTdfWMMdLhCaglzTwJdbVFlNGP0ku8AuM994WsmTGeeAiS338atVuDUr8CxQsT9UYjD7XcCZf8D3ncaFtHE1S6BTeuF46YewBaL5HbJUbFVyeJ3voqm4TMiapgaZPDh6zfGQIYr3HWP132ASRmakGO1RV3eUnD70zvjLUACRLTCSaj3/RsoLwKO7QaOFQHGc/XfFN8SaN0LSLv0k9oFUNfmv0g+eArY8p3kdslVsVWp4ne+CtbwGRGRUhpk8OGOu8eSL8MVdSeBeuoer8s2NPH3wR0dzuPrw1/qQ8lV7ovEOA0gAmcvVT111zvjaeKutwAJEHCX+r8o/GQveqp+hlq4dGdiYgHDtUCrbkDLbrVBR2IaILh+qEvNl+FKoBVbw2GJajCGz4iIlMTgQwJfHnZ1J4Emxmo9do/XVVpRjWkf/ehwHl+HZqQ+lAZlpOCWTsn2b/BTBqRj4VcH6302qcMVCTiPFqjAQbTyeuzrlrvxuuVuGLTVyOluwpCemUCLToBa+q+ip3wZ3kRDxVZvv4+BDp8RESlN9vTq4crX7va6KcdtDztA2mQ+2yTQTcVlvjbTwXGHPBA6t9cWABgCWG3x8a7f3Q6X2FgsVuBcOVSH/4sx6gLMiVkG3cq7oHklE29q/omh6u8lXxsAyox6TPomHvknm/kUeNjY8mU4pwn3Jhoqtnr6feRkUyKKBA0i+MjfU4qsl7f69J5xS3ei7/Nf2mt2+PKwsz20Vxf94WtTXZ4HAO7pVlsZ2N3DJpDVFp6GS2z2vDQU+MdV0K+8E89qlmJczBcoOHgRN516AqNNs/G65W4XrXbPdlTuumK/52EMyTRg0/Sb7dtN4zReH75ypy4PFXe/j6mJeuSN6hqaRhERSRT1wYdtwqeUB6wz56Jhzg87T0QApy+YJD0QpVi45SAS4zT1io7ZHjaPrCi0F6TzVDROBSuEc6VQ/b4Dd6n+iyEqaT0WX1e3hUVUwdq0HTZbrsN04yOYZJrqYsmr9E9bd6Ksv+oGXHNuv9plCwQ3x3sT7sX3nH8fl43vge0zbsGgjJQQtoqIyLuonvMhdcKnO66Wsfrau5B9bUu89+0RWcbhK6pMDuewrbaoMdcuO1XBioLCA5hfcMR+zLilO2HQVGF2s6+wXfc5UnAGmtdqj/+nFvjW0hn51l5er/265W582mg0Zt7cGX+TOYmVXKs33C2H9Sd1eaRkDw2HCbBERL6K6p4Pb/kQpAj02/ktnZL9mpvgri0CRDRCFWaqV6LvT09AveIu6P91E3bqJuH1mFcwdc0hlF1wzIFRZtJjctkw7LG2g0awQBTUsCZcgW8tGSixpiJVZ4IgITQqq6xRJHumnKs3XPUGFEy/yadzuOstY/ZQIiJ5RHXPh5z5EI6fPQ9YEgCLGWpYYIWA5IRYtym3ARExsKCH5hAaaY24KbsaD63ci1NiAg4g7VJPiOBwvJQhCxECLiAO16oOImZPbZZWFYAkCHjG/GC9Wie171FBgIi/mqYgyVSJr+aMAFQxGH1pyewr93Z1mb2z/rV9n9Aod1pyKZx7A3wpYheq7KGB1OxxxvL0RBTuojr4kPMbdfKakcC6fYgDcPDSafMv9sAkTIUA0V5TBAAEWAEIeF3zGhq9uxMAEAfgg0udH/mWHsg1PehivoT0x9uHlgHonjUC2sRUVOuaYcjy3z2mHBchwAQNytEMUDn+tbsbrnB9Ht+EU1pyKUKRPdRduv9gZ04lIgqWqB52seVDcP9oE5GM00gVTl8KGOoTYIUBJ9FT9XO914aod2KRJg+pOOOwP1U4ixdiFuM61S+wJqYBLTrBYrgO31oysNlyHbIyW2FLz+8wVpWP21Tf4L3rD+GNnqeQGiv90bbG2hfm3o8BXe+DNX0gDiNV8ntdTZ70ZTKtjaelvza2wMbVqgwp9WiCLRTZQ6euKnI5xFN3srPcbL0jhxcMr1e9mYhIaVH9r46nZFS2YY/jSMIrI10PO9Q+RFXIGdEX6s6HAdRm+Ow9fzMEiPjmiUEYoo9FHxPQ5ZkvAQh1yqMnAWagePJgxGljUGM024c5iu+qzQy6fEft9gu31x7Tb5jJa1l1OSauuktg5k8PhLsejbrCJS25FKHIHhpIzR4iokgU1T0fgPt8CCkJlx8yXr+dd0sHYpvYfyrQGGcRD+gTAW0c1BotbI/dQFYc1H2fAM9LRt2RemVv2UulnGfygHTJeSYiZVWGt94yuVodp43BBxNu8HiMHEuRiYjCUdQHH4C0FRDuciY4DwtI6a6WIz9E3qiufieQkvKAtH2zNiTqoYtR+3We1786CFG8/NmCkWeiymhG25kb7DlN5BbM7KFSh3hYSI6Iok1UD7vUJWUFhBzfzt1NHvQ1P4Rz/RXnnB511Q1uJg9Ix8e7fpeUVM1TNVxXxedcOXHOaP9zOPdo+MLWW+Y8ATdV5kmgUod4WEiOiKJNg+j5CCZ3kwf9yQ8hJRhyTh3v3Bshhatv1oMyUuqlLndFrnL14UZqT1ggpAzx+Fqzh4goEjD48MLXbn5vBdoCVXfYZ9uBEy6TYdXtjZDC3TfrusHOmSqT1/PsOnLG6zGRROl5KlKGeMJtKTIRkRwYfASJLQD5YMINsixt9JYMy0buariehKJcfaQvGfVUIC4clyITEcmBwUeQyTV5UGrq+Hu7K1cN11k0lKsPhWAM8RARhRMGHzKQsmzSxjbEEeg3dqlBTJtmjWT5Zu2pB8XGVq4+0nsjQiFSliITEcmhwTwZpNS7CKQmhm3yYFlFtd91TJyv72mOidQVEC3idejfMTngJF+e0qTbtsPlgcnaJkRE4Y09HzIJ9uRB76nja9l6I6R8s/bUY+EuEVvdZG0UfOxlIqJIxOBDRsGcPCg1GZacvRFylKsnIiLiVyUvfC11LmcdE2/DB+6SYaUk6FFWqUxWTOceFKk4FEJERDbs+fDAOYHXuKU70ff5L71WGg3m5MFQ90aw25+IiHzFp4Ub+XtKMWlFYb3Jo7ZS5+GUg0FK6ngiIqJwwZ4PF6Qk8MpdVxxV6cRdCcdeDTmK9hERUWgx+HDBWwIvljoPDX+HwYiIKLww+HCBpc7DQ926Ov8u+sNlHRvbMBgDECKiyMHgwwWWOg8/8zbua/DDYERE0SI8BvLDjBzZSsNNpE9Cde7xqKvuMFjv9GayXTPS7xkRUbhiz4cLLHUemTgMRkQUGRQLPhYuXIi2bdtCr9ejV69e2LFjh1KXUkQg2UrDcZVIQ8BhMCKiyKDIk/HDDz/E9OnTsXjxYvTq1Qt5eXkYPHgw9u/fj+TkZCUuqQg5s5VGi1ANRaQk6HC8siZqhsGIiBoyRXo+Xn75ZUyYMAHjx49HRkYGFi9ejLi4OLzzzjtKXE5RkVDqvCH0tDwxrDOA6B0Gawh/h0RENrIHH0ajEbt27UJWVtbli6hUyMrKwrffflvv+JqaGlRWVjr8EDlzV1VXiaJ9RESkLNmDj5MnT8JisSAlJcVhf0pKCsrKyuodP3/+fCQmJtp/0tLS5G4SRQlXdWy2z7iFgQcRUYQJ+WqXWbNmoaKiwv5z9OjRUDeJwlgkDIMREZFnsg8uN2/eHGq1GuXl5Q77y8vLkZqaWu94nU4HnU5Xbz8RERFFJ9l7PrRaLbp3747Nmzfb91mtVmzevBm9e/eW+3JEREQUYRSZVj99+nSMHTsW119/PXr27Im8vDxcuHAB48ePV+JyREREFEEUCT5GjhyJEydO4Omnn0ZZWRm6du2K/Pz8epNQIwFTbIdO3VotO0pOo1+HFiFsDRERyUWxhAJTpkzBlClTlDo9Rbn8PaXIWbvXvj1u6U4YEvWYObRTCFtFRERyYDYjCjv5e0oxaUVhvWymZRXVmLqqKBRNIiIiGYV8qS1RXRariNx1xS7TqLvaR0REkYfBB4WVHSWnUVrhvjotAxAiosjH4IPCyvFz7gMPIiKKDgw+KKwkx+tD3QQiIlIYgw8KKz3bJcGQqK9XvdaGydSJiCIfgw8KK2qVgJzsDAD1Aw3b9uIx3Vh6nogogjH4oLAzJNOARWO6ITnBseZPaqIei8Z0YxVbIqIIx6+OFJaGZBrQ58rm6DLnCwDAsvE90K9DC1axJSKKAuz5oLBVN9Do2S6JgQcRUZRg8EFERERBxeCDiIiIgorBBxEREQUVgw8iIiIKKgYfREREFFQMPoiIiCioGHwQERFRUDH4ICIioqBi8EFERERBxeCDiIiIgoq1XShsxWljcHjB8FA3g4iIZMaeDyIiIgoqBh9EREQUVAw+iIiIKKgYfBAREVFQMfggIiKioGLwQUREREHF4IOIiIiCisEHERERBRWDDyIiIgoqBh9EREQUVAw+iIiIKKgYfBAREVFQMfggIiKioGLwQUREREHF4IOIiIiCKibUDXAmiiIAoLKyMsQtISIiIqlsz23bc9yTsAs+zp07BwBIS0sLcUuIiIjIV+fOnUNiYqLHYwRRSogSRFarFceOHUN8fDwEQZD13JWVlUhLS8PRo0eRkJAg67npMt7n4OB9Dg7e5+DhvQ4Ope6zKIo4d+4cWrZsCZXK86yOsOv5UKlUuOKKKxS9RkJCAn+xg4D3OTh4n4OD9zl4eK+DQ4n77K3Hw4YTTomIiCioGHwQERFRUDWo4EOn0yEnJwc6nS7UTYlqvM/BwfscHLzPwcN7HRzhcJ/DbsIpERERRbcG1fNBREREocfgg4iIiIKKwQcREREFFYMPIiIiCqqoCz4WLlyItm3bQq/Xo1evXtixY4fH4z/++GN06tQJer0eXbp0wcaNG4PU0sjmy31esmQJ+vXrh6ZNm6Jp06bIysry+vdCtXz9fbZZtWoVBEHAnXfeqWwDo4Sv9/ns2bOYPHkyDAYDdDodrrrqKv7bIYGv9zkvLw8dO3ZEbGws0tLSMG3aNFRXVweptZFp27ZtyM7ORsuWLSEIAtasWeP1PVu2bEG3bt2g0+lw5ZVXYtmyZYq3E2IUWbVqlajVasV33nlH3Lt3rzhhwgSxSZMmYnl5ucvjv/76a1GtVosvvPCCWFxcLD711FOiRqMRf/rppyC3PLL4ep/vu+8+ceHCheLu3bvFffv2iePGjRMTExPF33//Pcgtjyy+3mebkpISsVWrVmK/fv3EO+64IziNjWC+3ueamhrx+uuvF4cNGyZu375dLCkpEbds2SIWFRUFueWRxdf7/P7774s6nU58//33xZKSEvHzzz8XDQaDOG3atCC3PLJs3LhRfPLJJ8XPPvtMBCCuXr3a4/GHDh0S4+LixOnTp4vFxcXia6+9JqrVajE/P1/RdkZV8NGzZ09x8uTJ9m2LxSK2bNlSnD9/vsvjR4wYIQ4fPtxhX69evcSHH35Y0XZGOl/vszOz2SzGx8eLy5cvV6qJUcGf+2w2m8Ubb7xR/Ne//iWOHTuWwYcEvt7nRYsWie3btxeNRmOwmhgVfL3PkydPFm+55RaHfdOnTxf79OmjaDujiZTg4+9//7t49dVXO+wbOXKkOHjwYAVbJopRM+xiNBqxa9cuZGVl2fepVCpkZWXh22+/dfmeb7/91uF4ABg8eLDb48m/++ysqqoKJpMJSUlJSjUz4vl7n+fOnYvk5GQ89NBDwWhmxPPnPq9duxa9e/fG5MmTkZKSgszMTMybNw8WiyVYzY44/tznG2+8Ebt27bIPzRw6dAgbN27EsGHDgtLmhiJUz8GwKyznr5MnT8JisSAlJcVhf0pKCn7++WeX7ykrK3N5fFlZmWLtjHT+3GdnM2bMQMuWLev9wtNl/tzn7du34+2330ZRUVEQWhgd/LnPhw4dwpdffon7778fGzduxK+//opHH30UJpMJOTk5wWh2xPHnPt933304efIk+vbtC1EUYTab8cgjj+CJJ54IRpMbDHfPwcrKSly8eBGxsbGKXDdqej4oMixYsACrVq3C6tWrodfrQ92cqHHu3Dk88MADWLJkCZo3bx7q5kQ1q9WK5ORkvPXWW+jevTtGjhyJJ598EosXLw5106LKli1bMG/ePLzxxhsoLCzEZ599hg0bNuCZZ54JddNIBlHT89G8eXOo1WqUl5c77C8vL0dqaqrL96Smpvp0PPl3n21eeuklLFiwAJs2bcI111yjZDMjnq/3+eDBgzh8+DCys7Pt+6xWKwAgJiYG+/fvR3p6urKNjkD+/D4bDAZoNBqo1Wr7vs6dO6OsrAxGoxFarVbRNkcif+7z7Nmz8cADD+Avf/kLAKBLly64cOECJk6ciCeffBIqFb87y8HdczAhIUGxXg8gino+tFotunfvjs2bN9v3Wa1WbN68Gb1793b5nt69ezscDwAFBQVujyf/7jMAvPDCC3jmmWeQn5+P66+/PhhNjWi+3udOnTrhp59+QlFRkf3n9ttvx4ABA1BUVIS0tLRgNj9i+PP73KdPH/z666/24A4ADhw4AIPBwMDDDX/uc1VVVb0AwxbwiSxJJpuQPQcVnc4aZKtWrRJ1Op24bNkysbi4WJw4caLYpEkTsaysTBRFUXzggQfEmTNn2o//+uuvxZiYGPGll14S9+3bJ+bk5HCprQS+3ucFCxaIWq1W/OSTT8TS0lL7z7lz50L1ESKCr/fZGVe7SOPrff7tt9/E+Ph4ccqUKeL+/fvF9evXi8nJyeKzzz4bqo8QEXy9zzk5OWJ8fLz4wQcfiIcOHRK/+OILMT09XRwxYkSoPkJEOHfunLh7925x9+7dIgDx5ZdfFnfv3i0eOXJEFEVRnDlzpvjAAw/Yj7cttX388cfFffv2iQsXLuRSW3+89tprYuvWrUWtViv27NlT/O677+yv3XzzzeLYsWMdjv/oo4/Eq666StRqteLVV18tbtiwIcgtjky+3Oc2bdqIAOr95OTkBL/hEcbX3+e6GHxI5+t9/uabb8RevXqJOp1ObN++vfjcc8+JZrM5yK2OPL7cZ5PJJM6ZM0dMT08X9Xq9mJaWJj766KPimTNngt/wCPLVV1+5/PfWdm/Hjh0r3nzzzfXe07VrV1Gr1Yrt27cXly5dqng7BVFk/xUREREFT9TM+SAiIqLIwOCDiIiIgorBBxEREQUVgw8iIiIKKgYfREREFFQMPoiIiCioGHwQERFRUDH4ICIioqBi8EFERERBxeCDiIiIgorBBxEREQUVgw8iIiIKqv8PC+A1dk48l3QAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.errorbar(x, y, ye, fmt=\"o\", label=\"data\")\n", "xm = np.linspace(np.min(x), np.max(x), 1000)\n", "plt.plot(xm, model(xm, *par), label=\"fit\")\n", "plt.legend();" ] }, { "cell_type": "code", "execution_count": 4, "id": "passive-cowboy", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "a = 1.10 +/- 0.18 jackknife=0.13 bootstrap=0.13 MC=0.18\n", "b = -0.72 +/- 1.09 jackknife=0.93 bootstrap=0.88 MC=1.04\n", "c = 10.52 +/- 1.22 jackknife=1.11 bootstrap=1.05 MC=1.19\n" ] } ], "source": [ "# now only return fit parameters\n", "def fit2(x, y, ye):\n", " return fit(x, y, ye)[0]\n", "\n", "# jackknife and bootstrap\n", "jvar = jackknife.variance(fit2, x, y, ye)\n", "bvar = bootstrap.variance(fit2, x, y, ye, size=1000, random_state=1)\n", "\n", "# Monte-Carlo simulation for reference\n", "mvar = []\n", "for itry in range(1000):\n", " y2 = 1 + 10 * x ** 2 + rng.normal(0, ye)\n", " mvar.append(fit2(x, y2, ye))\n", "mvar = np.var(mvar, axis=0)\n", "\n", "for n, p, e, ej, eb, em in zip(\"abc\", par,\n", " np.diag(cov) ** 0.5,\n", " jvar ** 0.5,\n", " bvar ** 0.5,\n", " mvar ** 0.5):\n", " print(f\"{n} = {p:5.2f} +/- {e:1.2f} \"\n", " f\"jackknife={ej:1.2f} \"\n", " f\"bootstrap={eb:1.2f} \"\n", " f\"MC={em:1.2f}\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.9" } }, "nbformat": 4, "nbformat_minor": 5 }