{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fitting problem 3\n",
"\n",
"This is a linear model. (It is not a linear function of $x$, but the model\n",
"is linear in the fit parameters $a_1$ and $a_2$.)\n",
"\n",
"The data for this problem is available in the file `fit_3.dat` which can be downloaded from http://www.eg.bucknell.edu/~phys310/hw/hw4.html. One option is to download the data file to a local directory, and then import it using `np.loadtxt()`. Another option is to look download it directly into the Jupyter notebook. I will use the second option here."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy import optimize\n",
"\n",
"import urllib # for importing data from URL\n",
"\n",
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Following is an Ipython magic command that puts figures in notebook.\n",
"#%matplotlib notebook\n",
"#for making pdf-file of this notebook, use instead following line\n",
"%matplotlib inline\n",
" \n",
"# M.L. modification of matplotlib defaults\n",
"# Changes can also be put in matplotlibrc file, \n",
"# or effected using mpl.rcParams[]\n",
"plt.style.use('classic')\n",
"plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures\n",
"plt.rc('axes', labelsize=16, titlesize=14)\n",
"plt.rc('figure', autolayout = True) # Adjusts supblot params for new size"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Options for getting data into notebook\n",
"+ Start by downloading the data files into the working directory. The `np.loadtxt` \n",
"function imports the content of the data file into a `numpy` array.\n",
"The `unpack = 'True'` option transposes the array so that each column is in \n",
"a separate array.\n",
"+ Download directly from a URL"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def f(x, a1, a2):\n",
" '''Linear model function with adjustable parameters a1 and a2'''\n",
" return a1*np.sin(2.*np.pi*x)+ a2*np.sin(4.*np.pi*x)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"g = urllib.request.urlopen('https://www.eg.bucknell.edu/~phys310/hw/assignments/fitting_3/fit_3.dat')\n",
"data = np.loadtxt(g)\n",
"\n",
"x, y, u = data.T"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"#Uncomment to look at the data.\n",
"#data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### (a) Is the suspected relationship relationship between $x$ and $y$ a linear model, or a nonlinear model?\n",
"\n",
"\\begin{equation}\n",
" y=f(x)=a_1 \\, \\sin(2 \\pi x) + a_2 \\, \\sin(4 \\pi x)\n",
"\\end{equation}\n",
"This is a linear model. (It is not a linear function of 𝑥 , but the model is linear in the fit parameters $a_1$ and $a_2$.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### (b) Perform a fit to these data using the generic optimize.curve_fit routine from scipy that we demonstrated in class, using the assumed functional form shown above."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Perform fit"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1.96501629 0.87923822]\n",
"[[6.39999997e-03 1.10687288e-11]\n",
" [1.10687288e-11 6.39999999e-03]]\n",
"0.07999999982280558\n",
"0.07999999993543226\n"
]
}
],
"source": [
"popt, pcov = optimize.curve_fit(f, x, y, sigma=u, absolute_sigma=True)\n",
"print(popt) # best fit parameters\n",
"print(pcov) # covariance matrix\n",
"for i in range(len(pcov)):\n",
" print(np.sqrt(pcov[i,i])) # uncertanties in fit parameters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Results\n",
"$a_1 = 1.96 \\pm 0.08\\quad\\quad\\mbox{and}\\quad\\quad a_2 = 0.88 \\pm 0.08$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAFbCAYAAACOHWQYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAxOAAAMTgF/d4wjAABG+0lEQVR4nO3dd3iN5xsH8O85GRKyJEKEBBVixogiqB2htNQqpWaL1Cgxa9RoqZWYkRq19yaCtDRoNT+rRltKY0VJrJCThOQkOef3R5oQTnL2ec+b8/1clwvvvPNeSe7zPO/z3I8kJiZGCSIiItKZVOgAiIiIxI7JlIiISE9MpkRERHpiMiUiItITkykREZGemEyJiIj0xGRKRESkJyZTIiIiPYk6mU6dOhWtWrXChQsXhA6FiIgsmGiT6ZEjR5CRkSF0GEREROJMpomJiVi/fj0mTJggdChERETiS6YKhQJz587FgAED4O7uLnQ4RERE4kumu3fvhr29PTp06CB0KERERAAAa6ED0Mbdu3exc+dOfP/99xodr1Ao8PTpU9jb20MikRg5OiIiEhOlUomXL1/Czc0NUql+bUuJmJZgO3r0KObPn58vMSoUCkilUrRq1QpTp07Nd/zjx4/Rs2dPU4dJREQisnPnTr1fG4qqZdqsWTP4+vrm2zZo0CCEhISgYcOGbx1vb28PALh69SrKlStnkhiLismTJ2POnDlChyFKfHa64XPTHZ+dbmQyGby8vPJyhT5ElUwdHBzg4ODw1nYPDw+VnypyW7COjo5wcnIyenxFia2tLZ+ZjvjsdMPnpjs+O/0Y4jWg6AYgERERmRtRtUxViYmJETqEIikoKEjoEESLz043fG6647MTHlumpBJ/OHXHZ6cbPjfd8dkJj8mUiIhIT0ymREREemIyJSIi0hOTKRERkZ6YTImIiPTEZEpERKQnJlMiIiI9MZkSERHpicmUiIhIT0ymREREehJ9bV5LEx2d8yc7Gzh7FmjYELCyAoKCcv4QEZHpMZmKTG7SlMkAZ+ecxMqVl4iIhMVuXiIiIj0xmRIREemJyZSIiEhPTKZERER6YjIlIiLSE5MpERGRnphMiYiI9MRkSkREpCcmUyIiIj0xmRIREemJyVSk5HI5gIT//iYiIiGxNq/IKBQKjBsXisjIeABeCAi4hw8+8MbChWMhlfKzERGREJhMRWbcuFBERLyL9PTxAIC4OCAi4gSAUISFjRc0NiIiS8WmjIjI5XJERcUjPb1lvu3p6S0RFXWXXb5ERAJhMhWRp0+fQibzUrkvJcULSUlJJo6IiIgAJlNRcXNzg5PTPZX7HB3/haurq4kjIiIiQITvTLdu3YqjR4/i0aNHKFasGGrVqoVhw4bBy0t1i60osbW1RceO3oiPP5Gvq9fO7gQ6dvSGra2tcMEREVkw0SVTT09PfPnll/D09ERaWho2bNiAr776Cps3bxY6NJNYuHAsgFBERu5GXJwXfHz+zRvNS0REwhBdMm3ZsmW+/w8cOBCDBw9GUlKSRXRzSqVShIWNx+TJcri7JyE21hWlSrFFSkQkJNEl09dlZGTg6NGj8PLygouLi9DhmFROl64H2LNLRCQ8USbT2NhYzJo1CxkZGShfvjzmzZvHggVERCQYUWagunXrYs2aNViyZAkqVKiAb775BllZWUKHRUREFkqULVN7e3uUK1cO5cqVQ7Vq1fDhhx/izJkzaNq0qcrjZ86cCUdHRwBAUFAQgoKCTBmuVqKjc/5kZwNnzwINGwJWVkBQUM4fIiLSXXR0NKKjowHAoIVuRJlM36RUKmFlZVXg/unTp6N8+fImjEh3uUlTJgOcnXMSq5OT0FERERUNrzeoZDIZwsPDDXJd0SXTlStXolmzZnBzc8OzZ8+wbds2ODs7o1atWkKHRkREFkp0yfTRo0eYMWMGkpOT4ezsDD8/P4SGhsLBwUHo0IiIyEKJLplOmzZN6BCIiIjyEeVoXiIiInMiupYpGY+6kcQcaUxEpBqTKeVRN5KYI42JiFRjMhWZ11uHjRsD06Zp3jpky5KIyDiYTEVGn8THliURkXFwABIREZGemEyJiIj0xGRqpnJqRiYYtHYkEREZB9+ZmhmFQoFx40IRGRkPwAsBAffwwQfeWLhwLJeZIyIyU0ymZmbcuFBERLyL9PTxAIC4OCAi4gSAUISFjRc0NiIiUo1NHTMil8sRFRWP9PSW+banp7dEVNRddvkSEZkpJlMz8vTpU8hkXir3paR4ISkpycQRERGRJphMzYibmxucnO6p3Ofo+C9cXV1NHBEREWmCydSM2NraomNHb9jZnci33c7uBDp29Iatra1J4lA3kpgjjYmI8uMAJDOzcOFYAKGIjNyNuDgv+Pj8mzea19jUjSTmSGMiItWYTM2MVCpFWNh4TJ4sh7t7EmJjXVGqlGFbpDktyqeQy90AvLq2upHEHGlMRKQak6mZyunS9YC2PbtKpRKXEi8h5k4M4pPjkZCagAcpD/Ds5TO42rkica8Ejy6WAWyqod67CejepQpCQ8chKyvrv5HE+ZNizkji3Zg1K7XQ/XPnyk3WDU1EZG6YTE3MGCu3KJQKHPnnCA5cP4DD/xxGckYyWldqjcolK6NB2QYoW7UsStqVxLJvduH2T32QJQ8EAPx7B1i8/AiO/NMaI8Z0h0xWXuX1U1K8cPPmTbUjjT08PHT7AoiIRI7J1MQMuXJLtiIbO//aiW9/+RbJ6cnoUaMH1ndZj+YVmsPWKn8rUS6XI+Ts0bxEmierAx5f3oBll5bhcVZlABPfuo+j47+oXLkynJzuITHx7Tg40piILB2TqQgplUrs/GsnZpycgReZLzC52WQMqDsAxayLFXhOYXNYi2XXQ8wnYRj+57c4uCUaisxXTeTckcQODg7o2NEb8fEn8hWVMPVIYyIic8RkKjKP0h5hSOQQnLl/BrNazkL/uv3faoWqkjuHtaCWZSm3Utjzw3KMKxmK3ft3416CK6QOV9AsyB3zF6wFoH6kMRcfJyJLxWQqIpHXI/FZ5Gd4z/s9/Bn8J9yKu2l8bu4cVnUty1cjiZ9i/cnTmPbbODRb1wzfd/oedT3qFjrSmIuPE5GlYjIVgczsTIw6Mgpb/9yK5R2Wo69fX0gkEq2vo+kc1pzEWhadq3VHt3rvY/ap2Wi6tinmtpmLEQ1H6DzSmC1XIiqqmEzNXKo8FT0PdkdiaiKuDLuCCi4VdL6WLnNYi9sUx+w2s9GhSgf02t0LJ+6ewKKWPwBw0fr+bLkSUVHFsjXmrMQjdNrdCpmKTJwccFKvRPq6Vy1LzZuWzbyb4dKwS0jPSsd7W+oBnucNEgsRUVHAZGqmbj2/CQxugndcfHD4k8NwtnMWOiSUKl4Kkb0jMchvGDCgJaJuHhA6JCIis8BkaobuPL+D93c2B65/iDUdthQ65cXUpBIpxrw7Edi3EZ8f6YPws+FCh0REJDi+MzUzOw49wrCz7eCd2RXeyaGY/rVE40E6Jh3gc60r9nX1QO+DHyI+OR7ftf0OUgk/mxGRZWIyNSOyDBnmP+iADnUaYHPXJZBqOWLX1AN8Gnk2wW+Df0OHLR2QmJaItR+uhZXUyng3JCIyU2xKCOTNNUHTs9LRZXsXuBd3x/ou60XTyqvqVhWnB53GufvnMODAAGQrsoUOiYjI5ETXMt28eTNOnTqFe/fuoXjx4mjYsCGGDh0KFxcXoUPTiKo1QTt18sKDgAt4kfkCkb0jNapoZE48HDwQ0z8GrTe2Rr/9/bCs9QaI8FuLiEhnovuN9+eff6JHjx7w9fVFWloali5dilmzZiEsLEzo0DSiak3Q5St+gt3va3Dr8GmUsC0hcIS6KeNQBj/3+xmtN7bG0KP9AOlGiPDbi4hIJ+LoS3zN3LlzERgYCG9vb1SvXh0jRozAxYsXkZqaKnRoasnl8v/WBG2Zb3uWPBCu/74HZxvhp7/oo4xDGcT0j8GfiZeAoF5Iz0gXOiQiIpMQfdMhOTkZtra2sLe3FzoUtQpbuSXzZVWzXxP09dHCjRsD06blHy2sUCgwd9oGpB9oASSWQNW6/TGwpz9CQ8dBKhXd5zYiIo2JOpnK5XJs3LgRQUFBsLIy/1Gk6lZuEXpNUHXJUt0Umze7sJNf5HRhSyShCAsbn3dczqCrp5DL3QCI6/0wEZEqok2m2dnZmDNnDgAgODi40GNnzpwJR0dHAEBQUBCCBKqqrunKLcaib7IszKsu7PH5tmfJA7F1z0bMnSuHtbX1W4Ovcgvts+VKRKYQHR2N6OhoAMibTWEIkpiYGKXBrmYiCoUCc+fORVxcHBYvXgynAiZTpqWloVOnTrh37x7Kly9v4ihVUygU+PjzL7HnwA0ok5rDp3KiwRNK7jzT5GTTFZJPSEhA/fqbkJg44a19EofR+OFwHfyx78l/LdeWefvs7E4gOPhcvpYrEZEpyGQyODs749ChQyhRQr/Bn6JrDiiVSixYsABXr17FwoULC0yk5io1MxXnakfiux3vA8rBiI3N6QIVe8sstwtblbKlnmPEsRHYc+DGW4Ov0tNbIirqrkE/IRIRmZrofoOHhYUhNjYWU6ZMAQAkJSUhKSkJ2dniKBYw5ugYVHGrgqENRkLblVvMWW4Xtp3diXzb7exO4OOPamKs/1j8+1j1B5+UFC8kJSWZIEoiIuMQ3TvTQ4cOAQC++OKLfNu3bdtm1iNhASDyeiT2XNuDP4L/MGqFI6EG+BS2+HhmZiaWOvZFcsrb55nD4CsiIn2ILpnGxMQIHYJOnrx4gs8jP8eS9kvg5ewFmczw91BVXcmUA3wKW3y8WLFi6N/DH8tXREOR+WqU0+uDr0xaqJ+IyIBEl0zNXUEJ4VLl0WhUvhH61elntHurqq4UEXECQKhJB/i8Wnw8//ZFYROQpZiD1duXIjO5Dny8U/KSPWD6Qv1ERIbCZGpgqhLC+ac/Y832g7j26TVItFwJRlMFTU3JGeCzG3PnygV/PyuVShG+dCoGjLqGhosD8GXPUIxoPljQmIiIDEF0A5DERp4txxdRX2Bmy5ko51TOaPcprLqSuQ3w8S1dHTiyGxN/GYUz/54ROhwiIr2xZfoGQ7+3W3phIYpZF8PIRiMNH+xrzL260ltutcW0prPx0Y6PcH7IeXg6egodERGRzphM32DQ93YutxF6ZjZ+6vcTrKXGfdRCV1fSRXC9L/H380vouqMrTg44iWLWxYQOiYhIJ0ymRqJUKoEOo9C92ido4tXEJPcsbGqKOZJIJPi+0/dosb4FRh4ZiVUfrBIsFo4kJiJ9MJkayeFbBwGvWMxsdt1k9yxsaoq5srO2w+4eu1F/VX00KtcIg+sLMyCJI4mJSB8cgGQE8mw5pp4aBxyfA1d7N5Pf/9XUFPNOpLm8nL2wo/sOjDo6CucfnBc6HCIirbFlagQR5yJQzKoYcHHQW/vUrdxiqVpXao3pLaaj285uONH7AoBSQodERKQxJlMDe/byGWadmoXV7begm+Ltx2vpSbMw45uMx5n7ZzDocG9AchSA+a9RS0QEMJlqTd1Aldm/zIZ/WX+0qcCMqS2JRIJ1ndehwcqGQMsZAL4ROiQiIo0wmWqpsIEqt57dQvi5cJz57IzRKh2ZM0N0YTsVc8LadtvxXlxTRN1ogN4NOhs3aCIiA2AyNaBJxyahT+0+8CvjhydPhFm5RUj6dmHnK9R/ZwL6bv0ex3qew+rls0S/3isRFW1MpgYSey8Wh/85jL+H/42QkAWCrdwiZm8W6lc8BdatOQoH23lYsvgrgaMjIioYk6kBKJVKfHX8K4QEhCBs5jazWLnFHBXWDdyqlepC/crM9tiw6wcsmC98oX4iooIwmRZAmwW2f779M648vIKdXXfivWHfmPXKLUIqrBs4IaHgQv2ylPJY/etqDG893IjRERHpjv2Ob1AoFAgJWYCAgLEANiEgYCxCQhZAoVCoPF6pVGJqzFSMbzIe2WnZolm5xdzkFupXxdPtOSaenoirj6+aOCoiIs0wmb4h971dXNwyABMQF7cMERHvYty4UJXH/3j7MG4m3cTIRiMLTQhmuXKLGckt1G9ndyLfdju7E+j5UQ182eRLdN/ZHanyVGECJCIqBLt5X6P1AtsSBb6NnYavmn0FB1sHABB05RaxV1cqrFC/AgrE/huLYYeGYdNHm4y6yLqljcImIv0xmb5GkwW2PTw8Xm2stg+P0x5iWINheZuEXLlFLEmzIIUV6pdCiq3dtqLeynpYdWEVhjYYatB755uWw1HYRKQlJtPXaLPAdrYiG2j1NcY1mgp7G/u87WJcucXcvCrUn3+7h4MHdnTfgQ5bOqCBZwP4e/rn7dN3CbU3p+VwFDYRaYPJ9DXaLLC998YOwDYN/WqpXjKsoIRA+mleoTmmt5iO7ru648KQC3C1z/mAo88Salp37xMRvYH9V29YuHAsgoPPwcdnBIB58PEZieDgc/m6aRVKBRac/hY4PhzIFi5WSzW+yXj4lfFDv339oFCqHmWtDU2694mICsNk+obcbtrY2DAA/REbm9PNl/veTKFQoPOAYYibUxH4Q6F26gwZnkQiwYYuG3D18VXM/XWu3tfjKGwi0heTaQEKWmB77NiFiNrWFdlPDgOYqHbqDBmHi50L9vTcg9m/zMbxW8f1ulZh03JMMQqbiMSPyVQLcrkcO/dfgzKzfb7tOe/W7v43rYJMpV7ZeljSfgl67+mNf2X/6nUtTbr3iYgKwgFIWnjy5AkeP3NXuU/l1BkyusH1BuO3e7+hx64eODngJHSdG8pR2ESkDyZTLVyRXUG2zR8q9/HdmvGpnv4iQZfAcFzOboaQ6BDMabZci3PfnjrDUdhEpAsmUy3MOzMPjVq54mKkMBWOLF3B01/s4fd8D/xX+cPPtRGAT7U4l4hIf6JLpqdOncL+/ftx48YNpKWl4dixY7CysjL6fU/Hn8aFBxdwZ/0dfDv1B0EqHFHBKrpUxJauW9B9Z3egjB+AOkKHREQWRHTJNCMjA/Xr14e/vz/WrFljsvvO/20+vnj3C7gWd+W7NTPV3qc9RjeYiNkfd0VS+jk4ObHbnYhMQ3TJNDAwEABw6dIlk93z2uNriI6LRkTHiLxtfLdmeIYo1D+u0RTMXncBg6J64cf+h2EtFd23OBGJEH/TaGDhbwvRp3YfeDp6qj1W7Cu3CMkQzygrMwvYOx/3G3bG5OOTMT9wvmGCIyIqBJOpGgmpD7Dljy24OPSiRsczaQoj36ovci+kL2uOpeV2wq+0H/rW6St0eERUxDGZqvH9xaUI8glCdffqQodChXhz1Zf424DNg+MYOHwgauyqgfpl6wscIREVZRaRTGfOnAlHR0cAQFBQEIIKaTq+3k3boKkMK85GoMuLw4guyRanuSpo1ZfMjDYoFd8UnTd3xrngc/BwYEENIksXHR2N6OhoADBo1TqLSKbTp09H+fLlNTr29W7ahb+tgu3ftbBjUFMjRkf6KmzVF5usOmjoKsdHOz5CTP8YAHamDY6IzMrrDSqZTIbw8HCDXFd0tXllMhni4uJw//59AEBcXBzi4uLw8uVLg95Hni3Hov8twoQmEwx6XTK8wld9uY91vddBAgkGHRgEpVJp4uiIyBKIrmX622+/Yd68eXn/HzZsGABg0aJFqFu3rsHus/3P7XCwdcAHvh8Y7JpkHOoWdXcq7oR9H+9DwzUNseDMtwCmvXUNfUdha1qukIiKJtEl0/bt26N9+/bqDyyAJr/0lEolFv1vEUY3Gg2pRHSN9yIv5z3HU8jlbsgtbJ9TgSq0wMpUZRzK4FDvQ2iyqglQtTTk8oF4vSi+vkmP5QqJLJvokqm+NPmld/LuScQnx6NfnX7CBEkq5Zv+Ai8EBNzLS5jqVn1RKBRYN+8oHPd2R+qTv1G3QTB6dq2Wdy4RkT4sLplqIiw2DEP9h6KEbQmhQ6HXvDn9JS4OiIg4ASAUYWE52wqqTPXmuffTgBUrfs53LhGRrviR/A3/PP0H0TejMfzd4Sr3R0cDISE579Ry362FhORsJ+N5Nf2lZb7tmizMXtC5GRmtcSDyJhd1JyK9sWX6hiVnlqBHjR4o51RO5X4OKBFGYdNf1C3MXti58Y8cEJ8QD58KPgaLlYgsD1umr3n28hnWXVqHMY3HCB0KvaHw6S+FL8xe2Lm2xeMw9OehyMjKMEicRGSZmExfs+rCKjTwbAB/T3+hQ6E35E5/sbM7kW+7JguzF3buwJ7vIjU7Fb329EKWIssIkRORJWA3738yszOx/NxyLOuwTOhQqADqpr/odu5XeJ4RjBbrW2DggYHY0GUDp0MRkdaYTP+z99pe2Eht8EFVFmkwV+qmv+h6rqu9K3769Cc0X9ccw6OGY0XHFZBIJMb8UoioiGEy/c/Ss0sxsuFIWEmthA6F1NBnYfaCzvVw8MCxfsfQbG0zjP9pPBYELhBNQmX1JSLhMZkCOP/gPC4nXkbUJ1FCh0IC8nb2xvF+x9F6Y2tkK7IRFhQmioTK6ktEwuPLIQDLzi5D/zr94WLnInQoJLAqblVwcsBJ7P17L0YdGaV1YfycOasJnLtKZGG0apnOmTMHxYsXh5+fH2rXrg13d3djxWUyj188wo4/d+Di0ItCh0Jm4p2S7+DkgJNovaE1MhWZWNFxhdpBSepKHRJR0aZVMvXz80NMTAyio6Mhl8vh4eGB2rVro3bt2vDz84OXl+qJ8eYot1j6D7+vRPMKzVHdvbrQIZEe9F315U0VXSrmJNSNrTHowCCs/mA1bKxsCjxek1KHRFR0SWJiYrRe4DE7Oxt///03/vjjD1y5cgUXL16EXC6Hu7s7unXrhm7dupnFp/G0tDR06tQJ9+7dy1sc/PUWRFxceViVOolOHStg79pws4iZ1Mt9N5icrP27QW3PTUhJQIctHeDp6ImdPXbCwdbhrWPkcjlq1x6LGzfenlZVteoI/PFHWKHzYA1Fn+dCZIlkMhmcnZ1x6NAhlCihXy12nbKHlZUVatasiV69emHOnDlYsWIFgoKC0LZtW+zfvx+TJ09Gdna2XoEZS24LIi5uGYCJyH5yGNE7emLcuFChQyMzVNaxLE4NPIWM7Ay02tAKj9IevXWMJqUOiaho0yqZymQy/Prrr3jy5Em+7ZUqVYKXlxc+++wzbNy4EX5+fti8ebNBAzWEgoult1JbLJ0sl1MxJxz+5DB8XH3QdG1TxCXF5duvSanD3AUSvvwSCAjI+fv1BRLU7Sci86bVO9NvvvkG9+/fx6NHj1CvXj20aNEC1apVQ2ZmJuLicn7BWFlZ4ZNPPkFERIRRAtaHPsXSybIVsy6GLV23YMJPE9BwdUPs6L4DgZUDAbwqVxgffyLfB7XXSx2qm75SlKe36DMPlnNoSSy0Sqa1a9fGggULcOvWLURHR2PTpk14/PgxrK2tMXZsTkm3//3vf0hOTkbJkiWNErA+clsQiYlv71NXLJ1IKpFiYbuFqF26Nrrs6IJvWn2DMY3HQCKR6FXqsKjT54NCUf6QQUWLVsm0WrVq2LRpEwIDAxEcHIzg4GCkpKTA2toa9vb2AICrV69i+/btGDp0qFEC1kduC+JufAwy0lvlbdekWDpRrv51+6O6e3V03dEVFxMvYlWnVbC3sde51CERiZ/Wo3llMhnOnz+P1q1bF3hMcnIynJ2d9Q5OXwWN5m3RuxfOn8xA+sMm+VoQHM1r3sytuzAxNRHdd3ZHckYytnXbhlqla6kdUavv/sKY+2jeovy1kTgZcjSv1uUEnZycCk2kAMwikRZECSXuBZxFxJh5GBjQgi0IEdEn8RnjHZuHgwdi+sfgm1PfoNGaRpjXdh4+9R0OQJgShLlzp+VyNwCvvqf53pHI+CyuNm/UP1HIVmbjw2rdAFjrVCydKJeNlQ1mtZqFtu+0RZ+9fRD1dzRQYg2AMiaLQV31Jb53JDI+i+vXXH52OYb5D4O11OI+R5ARNa/QHFeGXUFxmxLAiOrY+OcPWtf11VX+udMTEBe3DBER73LuNJEJWVQy/fvJ3zh19xQ+9/9c6FBIZDSZB1rSviTWddwO7NuAef+biVYbWuH6k+tGjavgudMtOXeayIQsqnm24twK9KjZA6VLlIZMJnQ0JCZadZXe+ABn+rXCgvPTUG9lPYxsOBJfvfeVUVYlMsTcab5TJdKfxSTTlIwUrL+0Hj9++qPQoZAFcLB1wKL2i9CvTj+M/2k83lnyDqa8NwWfVhsOwE6raxWW7Fq10n/uNN+pEunPYpLp5iubUdWtKhqVayR0KGRB6pWth58+/Qk/3vwRE49NxNIzy4AGk/Ayqz+cYK/RNQpPduqrLxGR8VnEO1OlUonwc+EY/u5wSCTCTFsgyyWRSBDkE4Tfh/6OyQ1nATWWoNbKCph9ajaevXyW71hdFhdfuHAsgoPPwcdnBIB58PEZieDgc6y+RGRCFtEyPXv/LB6kPECvWr2EDoUsVL7pK3EDYFPhOr4/eQBzms1Bb7/eGFR3EHYt+hWHDt2DtouLS6VSo1ZfMod3qgXNoSUyF6JNplu3bsXevXuRmpoKf39/jB07tsD3Q+svr8egeoNgb6NZtxqRob25eHjCXcDu4Ql08z4E23ov0OqTTyA/8z2QNQGAbouL53Tpehh87rSh3qnqkhDVzaElMhei/G48cuQINm3ahFGjRmH58uVIS0vDzJkzCzz+aNxRBDcINmGEZGlyp85MmwY0bpzzd+7UmcKmr5w7mY75LebD+1FHIKv9W/v3HriBjIwME34lhqdQKBASsgABAWMBbEJAwFiEhCyAQqFQey7n0JJYiLJlum/fPnTr1g3NmzcHAEyYMAF9+vRBXFwcfHx83jq+mVczVHatbOowyYIU1uWZkFD49JWbN28iNaWCyv13H5VAxe8qIrBOIAI8WgPOraBUekOokoWqqOsGfrNVrmmr+9WHkPzH5Myh3Y25c+UcYEVmQ3TJVC6X4+bNm/lWpfH09ISHhweuXr2qMpla/zUAISGvftnl/uDntiA4p46MSd3Sf5UrVy5wf5WyWVjWez1+uf8L1v2xEvjyM/iuckcjr4Zo6NkQDTwboHaZ2ijrUBbqEqyx3jsW1g2sT0Lk+sMkJqJLpjKZDAqF4q31Ul1cXPD8+XOV5ywa1RK1ar76P5MmmZK6xcMdHBwK3N+pUwUE+QYhyDcIExp8C+dSadhw+SL+en4WZ++fxbpL63Dr2S2UtC+J6q61gI41EfG7L/zKVUVVt6qo4FIBUkgFe++oT0Lk+sNkLLm9KS9eGO6aokumutQ7ffniBWQseUQGkPttpO2309dff46MjOU4enQHbt3yxjvv3EP79p74+usRkMlkavfn3TMTqOnoh4Byfvis5mcAgDR5Gq4/vY4L967i9MtrOH4jGmvPr8CtZ7eghBIlTvkg9cwSKDJfdbOGrziO1JczsXDeq+kz6r42Xfbb2NjAweGuyuNLlLgLa2vrQn8227Ytg/j4GKTnW384Bm3blkF6ejrS09MLPJeoIAEBQOPGSsw8PBdYaZhrar2eqdDkcjk6dOiA+fPnw9/fP29779690bt3b3z44Yd523LXMyWybIEAVFX+CgRwzAT39wCwFkCH17YdATAIgIpmp8rzawGoAeAagD80PI9IM4KsZyo0W1tbVK5cGZcuXcpLpgkJCUhMTESNGjVUnnP27FX4+pYzZZhUxCgUCkyZshxHjiTg9m0vVKp0Dx06lMXs2SO06iqVyQAvL+DevYIXBy9ovy7nJiYmonnzvXj48O3j3Uu3RPiuIUiySkJcUhyuPYxD1Jm/YVX6H9hZ26FaqWqo61EX/mX94e/pjzLWVVCxgpXWsec+u6NHD+Zrdc+efU3jZ/fkiRyVKz/DzZsluf4wGcSQyCGwSrfH1s/WG+R6okumANClSxcsX74cVatWRdmyZbFixQr4+fmpHHwEAA4OjnBisVHSQ0jIAqxd2yzvnebt28DatSdQrNhqjeeBvs7JqfC5mqr25w4gsrNzg5NTwQnl9XPt7Ozg7HxfZTIt6ZKIzo3H5w0AevJEDvfPnyI+0RFPcRt/PPoDvyf8jm3/bMO4U+NgJbEC+jbBuhtt8H611qjrURdWUiuNYg8Pn5pzffcknDmja1GJUihVinWDSf9CIo/THmP/7f048uEv2Ir1BolJlMn0/fffx7Nnz7B48eK8og3jxo0TOiwqooSeoqFP4QJ1g59sbW3fun6LZq+u/0ntTwAAWYosnLn9F5r1OYn/3T+OhWe/hZXECh2qdEDXal3R3qc9gMK7yYxVVIIsj76FRNZeXIum3k1RxdXXYDGJMpkCQJ8+fdCnTx+hwyALIPQUDV3naebKqdEbisjI3YiL84KPz795yVLT61tLrVHbvQ5wpg62/TgKxR2ycOHBBRy4fgBTfp6Cvvv6ok2FIKBGX8izPwRL/pGQCmu5tg3MRsT5CIQFhRn0nqJNpkSmIuQUDUO0igur3avN9V+fp+oktUWj8o3QqHwjzGkzB9ceX8Pmi7sQFTgB1Vd/gX51PsXg+oNRw131OAYiYyqs5Rp5/TCyFFn40PdDJD4w3NwYUZYTJDKl3K5SO7sT+babYpkzTVrFhZUyfN2rbtZX8WpyfU3KAVZ3r46Jjb8GlsZh7fvbkZiWCP9V/gjaHITjt47rNKWNyBjCz4VjqP9QWEsN25Zky5RIA+q6So1Fk1axPtW7NLm+Vt3MSilaeLfGB7Va43HaY4SfC0fP3T1RyaUSRtSbAEi6g5/hSSj/PP0HMXdisL7LeoNfm9/VRBrI7SqNjQ0D0B+xsTmJxNgVhIzdKlZ3fQAFFumPirqbb93VN9didS/hjhktZyB+dDz61+mPab+MA4b44+e7P+kVM5GuIs5HoGv1rvBwMPwYB4tomcrlmUKHQEWEECNSjd0qLuz6Dx8+VNsNXLp06UJHG5ewLYGRjUbi4yqfo0yncAyM+hgrLjXAvLbzUK9sPYN8DUTqpMnTsO7SOkT2jjTK9S0imXbuPAddu1blGogkSsZe/Luw6xuyG9jO2g6IHYtL+wch/PJ3aLq2KQbXG4zZbWbDqRgnj5Jxbf1jKyo4V0BTr6Z52zIz5YWcoR2LyCx374ZyDUQSPVUDiIx9fUN2A+cqaVcS8wPn4/Kwy/jr8V+oEV4D+//eb+CvhugVpVKJ8HPhGNFwBCQSSd6gurZtpxjsHhbRMgW4BiIJ6/V5b6qW/lO3X0j6dgMXNAe3ilsVHO93HBsvb8RnBz/DxssbsbLTSriXcDfiV0OW6H8PTuNu8t28IiSvelOGAlhlkHtYTDIFuAYiCUddUjSHpFkQfbuBCyORSNC/bn+8X+V9BEcFw+97P2zsshGBlQMBmPeHDBKPVZeWY1DdQShuU/yNudWGW03MopIp10Ak0p2qwVealCvUJCG6l3DHrh67sPbiWnTd2RXD/IdhdpvZCAqyZdKkAmm04L3jAxy6uQ9/Bf0FoPC51fqwmGRqign2RJZI3WhjTVuREokEg+sPRjPvZui9pzdifojB3o/3wtvZ27hfAImOVvWq/VehpXdb+LjmLIRSWG+KPiwimVasGIKPPvI1+gR7Iktk6NHGvqV8ETs4FqOPjkaDVQ2wu+duNK/Q3IARk9hpOoI89WUq4LsCA2usyduWvzelvsFisohkun//FNSpU0noMIgEYar3jrrMwS24IHkxRHSKQL2y9dBhSwcsCFyA4AbBkEgkhgu40Pvznay50qSetLW1NcaNC8XWPX8CT/ph6vEfcfaD63kt19zelAMHtuLWLcPEZRHJ1NbWRugQiARjzolB3VJaQ/yHoIZ7DXTf2R2XEy8jvGO4QWuq6ruUF5meJvWk58/flK/levONlmtub8rIkU/wzjurDRKXRcwzJSLxaubdDOeHnEfsv7Hosr0L0uRpQodEAsp956mKo+O/cHBw0Hjus42NAYufGOxKREWYpiuzkHGUdyqPUwNPIS0zDa035hTRJ8ukrpBISkqK2parMVhENy+Rvsy5q9RSuNi54Gifo+i3vx+arm2K6L7RqFSSYyEsUWEjyLOyslDc4bbK84w5PZLJlMjMGXsAkZgKIxSzLoZt3bZhzNExaLq2KX7u/zOqlaomdFhkYoWNILe1tYVHvce4czcaisxX38DGnh7JZEpk5oyd1MwxaRZGKpFicfvFcLB1QJNVLdEp6ThKZtbkaFwLpGoEeWZ2Jm43/A2dipfC1dORJlt/mMmUiERHIpHg29bfwtbKFuHnWmF/12NY6uPH0biEvdf2wt7WHnvXhuNZUrZRVlpShcmUiAplrt3AEokE01tOh7XUGp12twI8jgHg+qiWbunZpRjx7ghYSa1ga2sFU60/zGRKRIUSOmmqM6X5FGRlWmNGv0Bce3ISjZxqCh0SCeT8g/O4nHgZUZ9EmfzenBpDRKI35t2JwNkR6LI3EDeTbgodDglk2dll6F+nP1zsXEx+byZTIioaTkxHN99eaLupLf6V/St0NGRiCSkJ2PHnDoxsNFKQ+7Obl4j0Yoh3qhotpaWWBLObhyJDmYLATYE4OeAkSpcordf9WbtXPFacW4E277QRbKqURSTTBQsAFxf+ABAZgz4/V1otpaUBiUSC7zt9j777+qLDlg44OeAkHGwddL4/a/eKw8usl/j+wvfY1m2bYDFYRDKdNQsoX17oKIjoTZoupaUNK6kVNnTZgPe3vI8eu3rgYK+DsLFSvdiFMe5Pprfj2maUdSiLNpXaaHR8bo/DixeGi4HvTIlIEK+W0mqZb7uqguTasrWyxZ6ee/Ag5QGGHhoKpVJp0vuTKSkR8ftijG48WuMl+oKCgLAwYP58w0XBZEpEgtBkKS19ONs54/Anh3Hs1jHMPDnT5PcnE6n8I56mP8EntT8RNAwmUyIShLqltAxRkLycUzkc6XMES84swQ+//2Dy+5MJBIThM78vYGdtJ2gYokqmly9fxqRJk9ClSxe0atUK9+/fFzokItKRuqW0tClIntMlm6Cya7Zm6ZrY23MvRh0dhRN3Xt3LkPcn08pdEvGLGX9BUukUHh8JFnxJRFENQEpPT0fVqlXRrFkzhIaGCh0OEempsKW0NKFuNO6rqS2t4PFgCTpkdEPvtDP4ONAHQUH635+EkTvK+vODizFQ+QlWd84/BUqIEpiSmJiYt9/Mm7nExET07t0bmzdvRrly5Qo8Li0tDZ06dcK9e/dQnsN5iczWkyc5S2k9fqxdQfKQkAX/jcZtmbfNzu4EgoPP5RuNmzu1ZcTBsfjpbhRiB8eipH1Jje+fe35yMqfGmIvE1ERUWlIJ5z8/j5qldSshKZPJ4OzsjEOHDqFEiRJ6xSOqbl4iKppeLaWlXdeutqNxZ703H1XcqqDn7p7IzM7U6/4krKVnlqLtO211TqSGxmRKRKKky2hcK6kVtnbdioepDzEmeoyxQyQjSclIwYpzKzChyQShQ8ljFu9Mw8LCEBkZWeD+OnXqYPHixTpff+bMmXB0dAQABAUFIYhlkIhEL3c0bmLi2/sKG43rWMwRB3sfRINVDVC/bH0MqjfIyJGSoa3+fTWqlaqGZt7NtD43Ojoa0f+NVDLkXGKzSKZDhgzBp59+WuB+GxvV1Us0NX36dL4zJSpickfjxsefeOudqbrRuBVdKmJH9x34cPuHqOFeAzWcGhd6L8PUDiZDyMzOxKL/LcKS9ks0LtLwutcbVDKZDOHh4QaJyyySqYODAxwcCq6fSUSkij6jcdu80wazW89G1x1dcaL3BQBl3zrG0LWDSX/b/9wOe2t7dPbtLHQo+ZhFMtXUy5cvcf/+fTx58gQAcPfuXbx8+RKlS5eGE4fYEVkcqVSKsLDxmDw5ZzRubKx2o4G/bPQlfk/4HZ8e6gZYxQAolm8/a/cKo6DVetq1U2L+7fkY12QcrKRWQoeZj6iS6fXr1zFmzKtBA1OmTAEATJw4Ee3btxcqLCIS2KvRuNqdJ5FIsLLTSjT9oTnQYRSAlXn7Xo0Wzp80c0YL78bcuXKO/jWSglbrOfLPUTz68xH61ekndIhvEVUyrVu3LmJiYoQOg4iKEHsbe2zqtAe1btfH5r8a44uAgQA0Gy3s4eFhylAt3tzTczGq4SjBSweqIqpkSkRkDF5O3sDu7Rjv1BmNK9ZB/bL1dR4tTOrpsuj6qbuncDnxMg72OmjaYDXEZEpEBAC32mJco6notrMbLgy5AFd7V51HC1PhdFl0ffYvszGy4Ug42zmbJkgtMZkSUZGn6dSWMe9OxOUnZ9Bnbx8c6n2ItXvNxIXEszgdfxpbum4ROpQCMZkSUZGleSH8nILo07+WopzVBhwo8y6+OfUNZrScoddoYTKMhWdmI7hBMEoVLyV0KAViMiUiwRh7dQ91U1tU38cZwx7uQcAPAWhcvjHa+7TXebSwkHR5L2mWylxGTPxP+OGjleqPFRCTKREJxpi/2PWZ2lK7TG2s6LgCfff2xe9Df4eLxNs4QRqRLu8lzdJ7c/BprcHwcDDvkdMs4UFERZIuhfBf169OP3St3hU9d/WEPNtwNVxJczeS/gaqHcAof/MvkMFkSkRFUu7UFlU0ndqytMNSZCoyMeUUBxwJYcGZb4ErfXOmLpk5JlMiKpJyC+Hb2Z3It12bqS121nbY3WM3dlzbDNTabqRISZW/Hv2FA//sBk5NFToUjfCdKREVWYaY2lKpZCWsDNqIXrI+uJFUFw2cqhkvYMoz/cR0fOI7AOueF/tvapN5j/5iMiUi0VI3GljfQvi5OlT+ADj3Bfp5dse5IWdQwraEEb4aynX+/nnsX3IGXomdAGwSxWo9TKZEJFqajgY2yNSWn7+F68exCI4KxoYuG3RaS5M00/PzUcC5H3Ansx0AcazWY54pnojI3Cissfb97fjx5o9Y8/saoaMpEnK6bxP++zvHqZuncOesE7L/S6S5cqY03c13rDlhMiUi0pCHQ1ls67YNo6NH42LCRaHDES2FQoGQkAUICBiLnG7csQgJWQCFQoEpUVNQXNlA5XmaTGkSCrt5iYi00KpSK0x5bwq67+qOC0MuwMXOReiQRKegylT3kr/ElUpX4OFaGzdV5ExzXq2HyZSIqAAFDXAKbDcJ1Uv9hg6r+qPxnX1QZEtNXrJPrOUCC6tMFRkVislbx+J5kj0iIsS1Wg+TKRFRAQpOTFI0erkR/qv8UbbbQgyrPcHkJfs0LReo6Yo5plJYZarMF1XQt0pfVFxYEWJbrYfJlIhIB672rtjVYxdarG+BWiUbAWiRb7/QLUd1K+YIpbBF10u7PkH5MuUNNqXJlJhMiYh01MCzAcLahWFg1MeA4+8APPP2CV1oXt2KOULJrUz15qLrEpuj6PVRrXzduGJarYejeYmI9DDEfwjaVAgCephPQfxX7yVb5ttuLtNLFi4ci+Dgc/DxGQHgO0hLtUe3Tw8jNHScoHHpg8mUiEgPEokEYW0iANtUTD1lHslA3xVzjC23Gzc2NgxocwPt5ttg1w9Lzba6kSbEGzkRkZkoblMc2LEX269twuYrm4UOxyAr5pjCndRbQOPt+KbVAqFD0RvfmRJRkaWudq9BPXsHazpswYBDPeFXxg9+ZfwMfAPNFfRe0pymlyiVSoyPGQn8/jmquop/8QAmUyIqskw957Jdpfcxvsl4dN3RFWc/PwtXe+FagIZYMceYdv61E1ef/AH8vFvoUAyCyZSIyICmtZiGCwkX0Gt3LxzucxhC/Zo15+klsgwZxkSPwZwWYfgsw1nocAyC70yJiAxIKpFic9fN+Ff2LyYdmyR0OK9NLzGPRAoA036ehuru1dHdt7fQoRgMW6ZERAbmVMwJB3odQMM1DeHrXBdAX6FDMhsXEy5i9e+r8fvQ34vUMnZMpkREelJVsq+KWxVs77Yd3XZ2AzyrAVC9Eoox7i20gqo/BbZTYOa/wQgJCEG1UtUgkwkdqeEwmRIR6Uhdyb4gnyBMajwD0552QULqWTg5eaq9pqHuLaSCqj+Fn43Aw+sPMfm9yYLGZwyiSqaHDh3C0aNHcefOHVhbW8PPzw/Dhg2Dp6fhvkGJiDSlScm+obVHYtr8c+ixpyNih/6KErYlTHZvcxKXFIeJxybiYO+DOfNyC2HSKU0GIqpkevnyZbRr1w41a9YEAKxevRqTJk3C2rVrYW0tqi+FiESusKXEoqJ2Y86cdEyevCyn5RhXH/9cyoDfkSD8fTAGNtY2Rr333LlysxpwlK3IRv/9/TGo3iC0rtRa7fHmnDQLIqoMNGXKlHz/Hz9+PLp37467d++icuXKAkVFRJZIXcm+L7+ci40bW+YlvPSHwK2fjqBxz49wYe8ho947KSkJHh4eet3DkJZdCMWTF08wt+1coUMxGlEl0zclJycDAJxMuRQDEREKX0rMweEuYmKkbxWaR1YHXPplCZb9tgwjm4w0yr01LReozxJxWp1b+g/M/d8M/Nz/57zuXTF246oj2mSqVCrxww8/4N1334W7u7vQ4RCRhSmsZF/Llm6IjHRUeV5JaRNMjJwIL1cvdKnWxeD31rRcoD5LxGm8MHm2HPioH4LrjUbj8o3fOr8oMYtkGhYWhsjIyAL316lTB4sXL863bcWKFbh9+zaWLVtm5OiIiFQrqGTfnDlTcPLkeJUtRzeXx4jovQF99/bFoU8OoWXFlga9t7mUCwSA6b9MBCRKTGo8XehQjM4skumQIUPw6aefFrjfxib/y/rVq1fjxIkTWLp0Kdzc3NRef+bMmXB0zPmUGBQUhKCi9pGIiARRWMm+wlqOPfx6ICU7BZ23d8aJ/idQr2w9g97bHOz4cwe2XF0H7DyPYvOKCR1OnujoaERHRwOAQdd1lcTExCgNdjUT2LBhA/bv34/FixejQoUKhR6blpaGTp064d69eyhfvryJIiQiS5Pb3Zmc/Kq789U80LtvtRxz54HOPz0fobGhOD3oNHxcfQx2b0Dz95oFna/Pvf969Bca/9AYa9pvRa/6H+h0bVOQyWRwdnbGoUOHUKKEflOWzKJlqqmtW7di27ZtmDVrFhwdHfMWuHV0dHyr9UpEJCRNWo4Tmk7A47THaLuxLU4MOIGKLhUNdn+h3kvKMmTourMrRjcajQ6VPzB9AAIRVTI9ePAgMjIyMHHixHzbFy1ahLp16woTFBFRIV4Vmle9f37gfKRnpaPl+pYGT6implQqMfDAQFR0qYgZLWcgLVXoiExHVMl0+/btQodARGRQEokESzsshVQiRYv1LXCi/wlUKllJ6LB0mjoz5ecpuPDgAs4POQ8rqZVZ1g02FlElUyKioiJ/spLg3YaL4VZCgkbft8D/hp3AP2ff0XkeqCFoO3VmUewirP59NX4d+Ctc7VwRErLALOsGGwuTKRGRAN5OVhI4Oi7C2B+leG/de4j6JAphQXV1mgdqatuvbsL0E9Pxc/+f4VvKFyEhC0RVN9gQiuZHBCIiEZJIJAhtF4oR745Ai/Ut8NPNn4QOSb0qURhzfBj2frwXDTwbvFY3uGW+w3LqBt816HQUc8JkSkRkRiQSCb567yuEvx+OLju6YMtf64UOqUDH7hwFuvfCiqD1aPtOWwCa1Q0uitjNS0Rkhvr69YWnoye67ugKtLqJbMUMAFYGv4+ug4TWX1qPLw58AewNRceQznnbDVE3WIzYMiUiMlOtK7VGdM/TQM2d6LK3HRJTVWQoHSkUCoSELEBAwFgAmxAQMBYhIQugUCgKPU+pVOLbk99iyIjpKLmuN3A9Jd+5uXWD7exO5DtPm7rBYsSWKRGRGateqiaw6jzKbByGut/XxdZuWzVaE1QdXRYXz8zOxOijo7F+XjQk59fiQUYbleeKoW6woTGZEhHpyGRLickdsbr9Zuy6+QM+3PYhRjQcgQD51zh5rLhO99ZlcfEbSX9j2Pa+SM9IR5kHbXH7v0Ra0LnmXDfYGJhMiYh0ZMqSfRKJBJ/V/wyNyjXCkENDsC2lOpYMWYLOvp0hkUi0upY2i4srlUqgYThabp2E4AbB+KL6F2gyZ5dG56qr/lSUMJkSERmBsVqttcvUxulBp7H+0np8Hvk5Vl1YhdB2oajuXl3ja2g6SOhS4iWEHJkANLmOHZ0PoWPNlpDL5RY5wEgdDkAiIjKCoCAgLAxYsgSIjc35OyzMMC1ZqUSKQfUG4fqI66jkUgn1VtZDp62d8PPtn3NakmqoGyT0x5M/0Hl7ZzT5oQmqudUAIq7gPa+WGp1bVAcYqcOWKRGRSLnauyK8YzimNp+K8HPh6LmrJ8o5lcOAOgPQrnI71HCvUWAX8JuDhCq+cwc+AWn4q+5DNF8/A8P8h2Hllythne6KiIynkMvtkTt1xhIHGKkjuvVMtcH1TInI3KlbT/TJEznc3Z/i8WM3tYN4XmS+wJYrW7Dv7304efckXOxcEPhOIGqVroUyJcrAw8EDpUuURlpmGh6kPEBCSgKuP7iD8L0/QVrpLzSu2BidqnTCoHqD4F7c/b/1WOP/S5hv19fNiS0Jjx+rHmCkz1qppmCx65kSERU1BRVNeLW4uObF4ovbFMfn/p/jc//PkZGVgd/u/Yafbv2Es/fP4mHaQySmJuJR2iOUsCmBso5lUdahLEoV8wQuT8TNhe1RsYxb3rU0qa9rSQOM1GEyJSISgLpkqcs80NcVsy6GVpVaoVWlVoUeJ5MB67oBrvavtukydcbScQASEZEAcpNlXNwyABMQF7cMERHvYty4UMGLxVtqfV19MJkSEZmYumSZkJAgaDLLnTqjiiVPfykMkykRkYmpa/lJJBJBkxmnv2iP70yJiExMXdEEDw8PdOzojfj4E/lar6ZMZpz+oh0mUyIiE8tt+RWWLIVOZlKpVOf6uiarWWxGmEyJiASgLlnqk8wMSZfpL0U5aRaEyZSISACaJkvO5RQHJlMiIgExWRYNHM1LRESkJ7ZMiYgskD6DhCxxgJE6TKZERBZIn8RnyUmzIOzmJSIi0hOTKRERkZ6YTImIiPQkqnemR48exa5du5CQkACpVIqqVatiyJAhqFatmtChERGRBRNVMnV1dcVnn30Gb29vZGVlYc+ePZgwYQK2bdum9yrpREREuhJVN2/Dhg0REBCAcuXKoUKFChg2bBhSUlJw9+5doUMrcqKjo4UOQbT47HTD56Y7Yzy76GggJCRn2kvu9JeQkJzt9DZRtUxfl5WVhUOHDsHJyQne3t5Ch1PkREdHI4hj33XCZ6cbPrf8tJnLaYxnx+kv2hFdMr116xaGDx8OuVyOkiVLYv78+XBwcBA6LCIirahLlkxm4mIWyTQsLAyRkZEF7q9Tpw4WL14MAPDy8sKaNWsgk8kQFRWFWbNmISIiAk5OTm+dp1QqAQApKSmQyWRGib2oksvlfGY64rPTjaU9t4CAnD+qaPsYLO3ZGUruM8vNFfqQxMTE6H8VPaWmpuLly5cF7rexsYGLi4vKfZ9++im6dOmCbt26vbXv8ePH6Nmzp6HCJCKiImjnzp1wd3fX6xpm0TJ1cHDQuatWoVDAyspK5T43Nzfs3LkT9vb2kEgk+oRIRERFjFKpxMuXL+Hm5qb3tcwimWpq48aN8PPzg4eHB1JTU3Hw4EEkJyejcePGKo+XSqV6f9ogIqKiy1BjbkSVTFNSUjBv3jw8ffoUDg4O8PX1RWhoKDw8PIQOjYiILJhZvDMlIiISM1EVbSAiIjJHourmVWXr1q3Yu3cvUlNT4e/vj7Fjx8LV1VXlsS9fvsTSpUtx6tQpWFtbo127dhg2bFiBA5iKMk2fm0wmw9q1a3Hu3Dk8efIEpUqVQlBQEPr06WORzw3Q7nsuV1paGgYPHoyHDx/i2LFjfHYaPrvjx49j69atuHfvHpycnNC9e3f06tXLhBGbB22e2+3bt7FixQpcu3YNVlZWqFOnDoYPH44yZcqYOGphnTp1Cvv378eNGzeQlpam9udO3/wg6pbpkSNHsGnTJowaNQrLly9HWloaZs6cWeDxixcvxtWrV7FgwQJMnz4dMTEx2LBhgwkjNg/aPLenT5/i+fPnGDlyJNauXYvhw4dj37592Lx5s4mjNg/afs/lWrp0qcVX6tL22f34449YunQpevTogXXr1mHOnDnw9fU1YcTmQdvnNnXqVDg4OGDFihUIDQ1Famoqvv32WxNGbB4yMjJQv3599O7dW6Pj9c0Pok6m+/btQ7du3dC8eXP4+PhgwoQJuHLlCuLi4t46NiUlBceOHcPIkSNRo0YN1K9fH4MGDcKBAweQnZ0tQPTC0ea5VapUCTNmzEDjxo1Rrlw5NGnSBN27d8fp06cFiFx42jy7XL/88gvi4+Px8ccfmzBS86PNs8vKysL333+P4OBgtG/fHuXKlUPVqlVRr149ASIXljbP7fnz53jw4AH69OkDb29v+Pj4oHv37rhx44YAkQsrMDAQffv2Rc2aNdUea4j8INpkKpfLcfPmzXw/XJ6envDw8MDVq1ffOj73m6lu3bp52+rXrw+ZTIb79+8bPV5zoe1zUyU5ORmOjo7GCtFs6fLskpKSsHz5ckyaNMliu3YB3X5enz17huzsbAwcOBA9e/bEd999h+TkZFOGLThtn5uTkxPKly+PH3/8EXK5HC9fvsTx48fRoEEDU4YtOobID6JNpjKZDAqFAiVLlsy33cXFBc+fP3/r+GfPnsHBwQHW1tb5jgWg8viiStvn9qYHDx7g8OHD6Nixo5EiNF+6PLvQ0FB07doVFSpUMEGE5kvbZ5eYmAgg513h0KFD8fXXXyM+Pt7iuiu1fW5SqRQLFizA+fPn0aFDB3Ts2BEPHjzAV199ZaKIxckQ+UG0yVTbWoqqjrfEqkj61KB89uwZJk2ahNatW6N169YGjEoctH12R44cQXJyMnr06GGkiMRD22enUCgA5JQLbdy4MWrVqoWxY8fi/PnzePTokTFCNEu6PLfFixejQoUKCA8Px5IlS1C8eHGL+xCiLUPkB9GO5nV2doZUKsWzZ8/ybX/+/LnKOr6urq5ITU1FVlZW3qeP3HMLqvtbFGn73HIlJydj3Lhx8PX1xejRo40bpJnS9tldvnwZ165dQ2BgYL7t7dq1w+jRo/HBBx8YM1yzou2zy22JvT5oK/ffjx49QunSpY0XrBnR9rldvHgRFy9eRGRkJGxtbQEAX331FXr06IFbt27hnXfeMUXYomOI/CDalqmtrS0qV66MS5cu5W1LSEhAYmIiatSo8dbxVapUAZDzCy7XxYsX4eTkhHLlyhk9XnOh7XMDcl7Ojx8/HmXLlsWkSZMglYr220Yv2j67wYMHY82aNXl/xo0bBwBYuXIlWrZsaaKozYO2z87X1xfW1tb53lfl/tuSpnho+9zS09MhkUjy/Yzm/ju3tU9vM0R+EPVvxS5dumDPnj345ZdfEBcXhwULFsDPzw8+Pj54/Pgx+vXrh2vXrgHIeTHfpk0bLFu2DNeuXcPFixexdu1adO7c2eIGhmjz3NLS0jBhwgRYWVlh5MiRSE5ORlJSkkW9Z36dNs/O3d0dlSpVyvtTtmxZADkjpC1xAJc2z87BwQFBQUFYt24drly5gps3b2Lx4sVo1KiRxdXb1ua51axZEzY2NggNDUV8fDxu3ryJhQsXwtPT0+Le28tkMsTFxeV9CIuLi0NcXBxevnxplPwg2m5eAHj//ffx7NkzLF68OG8yc+6n/+zsbNy7dw8ZGRl5x48ZMwZLlizBuHHjYGVlhXbt2qF///5ChS8YbZ7bP//8g7///hsA8k2WL1OmDLZv32764AWm7fccvaLtsxs5ciTCw8MxZcoUWFlZoWHDhhgxYoRQ4QtGm+fm4uKCuXPnYtWqVfjiiy9gZWWFGjVq4LvvvoONjY2QX4bJ/fbbb5g3b17e/4cNGwYAWLRoETw8PAyeH1ibl4iISE+i7uYlIiIyB0ymREREemIyJSIi0hOTKRERkZ6YTImIiPTEZEpERKQnJlMiIiI9MZkSERHpicmUiIhIT0ymREREemIyJSIi0hOTKRERkZ6YTImIiPQk6iXYiEi9gwcPQiaT4c6dO2jXrh0eP36MpKQkxMXFYcSIERa3PiiRMbBlSlSEHTlyBJUrV0bfvn3RuXNnzJo1CzY2NvD19cWpU6dw69YtoUMkKhLYMiUqwmQyGWrWrAkAePr0KSQSCdq0aYPMzEwsWrQIdevWFTZAoiKCyZSoCPv444/z/n3lyhX4+fnBysoKVlZWTKREBsRuXiILcfHiRSZQIiNhMiUqorKzs3H+/HlkZ2fj2bNnuHPnDvz8/PL2b926VcDoiIoWJlOiIioqKgrjx4/Hw4cPERMTg2LFisHZ2RkA8Ouvv6JixYrCBkhUhEhiYmKUQgdBRIZ3+/ZtbN68Gd7e3qhcuTJevHiBCxcuwNPTEx4eHggKChI6RKIig8mUiIhIT+zmJSIi0hOTKRERkZ6YTImIiPTEZEpERKQnJlMiIiI9MZkSERHpicmUiIhIT0ymREREemIyJSIi0hOTKRERkZ6YTImIiPT0f+FYINjJADzZAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure()\n",
"plt.errorbar(x,r,u, fmt = 'ko')\n",
"plt.xlabel('x')\n",
"plt.ylabel('residuals ($y_i - f(x_i)$')\n",
"plt.axhline(0)\n",
"plt.xlabel('$x$')\n",
"plt.ylabel('residuals ($y_i - f(x_i)$)')\n",
"plt.xlim(-0.1, 1.1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### (d) Determine the goodness-of-fit parameter $\\chi^2$ for this data set and model."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Goodness of fit\n",
"\n",
"$$ \\chi^2 = \\sum_i\\left(\\frac{y_i - y(x_i)}{\\alpha_i}\\right)^2 $$ "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"36.61459756007561 0.7472366848995022 51\n"
]
}
],
"source": [
"chi2 = np.sum(r_norm**2) # Chi-square\n",
"chi2_nu = chi2/(len(r)-2) # reduced chi-square\n",
"print(chi2, chi2_nu, len(r))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###### Conclusions:\n",
"\n",
"#### (e) Given your analysis, does your fit to the data appear to be reasonable? Comment briefly. & (f) What are your resultant values for $a_1$ and $a_2$? (Include uncertainties.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The goodness of fit parameter $\\chi^2$ is less than the number of data points (but not much less), or, equivalently\n",
"the reduced chi-square $\\chi^2_\\nu$ is less than 1 (but not much less than one), so the data is consistent with the \n",
"assumed model. Also, in Fig.2 the residuals fluctuate around $0$ and show no systematic dependence on $x$.\n",
"\n",
"Therefore we can use the best-fit parameters found above:\n",
"\n",
"$$ a_1 = 1.96 \\pm 0.08\\quad\\quad \\mbox{and}\\quad\\quad a_2 = 0.88 \\pm 0.08 $$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Version information\n",
"`version_information` is from J.R. Johansson (jrjohansson at gmail.com); see Introduction to scientific computing with Python for more information and instructions for package installation.\n",
"\n",
"`version_information` is installed on the linux network at Bucknell"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"%load_ext version_information"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"application/json": {
"Software versions": [
{
"module": "Python",
"version": "3.7.8 64bit [GCC 7.5.0]"
},
{
"module": "IPython",
"version": "7.17.0"
},
{
"module": "OS",
"version": "Linux 3.10.0 1127.19.1.el7.x86_64 x86_64 with centos 7.9.2009 Core"
},
{
"module": "numpy",
"version": "1.19.1"
},
{
"module": "scipy",
"version": "1.5.0"
},
{
"module": "matplotlib",
"version": "3.3.0"
}
]
},
"text/html": [
"
Software
Version
Python
3.7.8 64bit [GCC 7.5.0]
IPython
7.17.0
OS
Linux 3.10.0 1127.19.1.el7.x86_64 x86_64 with centos 7.9.2009 Core