{ "metadata": { "name": "", "signature": "sha256:d15af7fbf18219f2a9a5d886cb016c9b88c4277c4c37d4bf49874ad614185435" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Introduction to Scipy: Fitting data" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "We have talked about the Numpy and Matplotlib libraries, but there is a third library that is invaluable for Scientific Analysis: [Scipy](http://www.scipy.org). Scipy is basically a very large library of functions that you can use for scientific analysis. A good place to start to find out about the top-level scientific functionality in Scipy is the [Documentation](http://docs.scipy.org/doc/scipy/reference/)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Examples of the functionality include:\n", "\n", "* Integration (scipy.integrate)\n", "* Optimization/Fitting (scipy.optimize)\n", "* Interpolation (scipy.interpolate)\n", "* Fourier Transforms (scipy.fftpack)\n", "* Signal Processing (scipy.signal)\n", "* Linear Algebra (scipy.linalg)\n", "* Spatial data structures and algorithms (scipy.spatial)\n", "* Statistics (scipy.stats)\n", "* Multi-dimensional image processing (scipy.ndimage)\n", "\n", "and so on." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "This week, we will take a look at how to fit models to data. When analyzing scientific data, fitting models to data allows us to determine the parameters of a physical system (assuming the model is correct)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "There are a number of routines in Scipy to help with fitting, but we will use the simplest one, ``curve_fit``, which is imported as follows:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "from scipy.optimize import curve_fit" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The full documentation for the ``curve_fit`` is available [here](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html#scipy.optimize.curve_fit), and we will look at a simple example here, which involves fitting a straight line to a dataset." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "We first create a fake dataset with some random noise:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.random.uniform(0., 100., 100)\n", "y = 3. * x + 2. + np.random.normal(0., 10., 100)\n", "plt.plot(x, y, '.')" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "pyout", "prompt_number": 3, "text": [ "[]" ] }, { "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGG1JREFUeJzt3X/MnWV9x/H3hx+dFFYfCaQtUGnjSqALW8GCbtNxNrWD\nZqH4RxHNtgrMx8RNmYktLcnWJ5ooYnAuMZIwAStZuxXcSHHjR+t6RJ0UUQqV0rUQynzUPgylFsai\nZf3uj/s+PKeH85zfv6/PKznhPvevc5879Huu53t9r+tWRGBmZqPtuH5fgJmZdZ+DvZlZAhzszcwS\n4GBvZpYAB3szswQ42JuZJaBmsJf0Bkk7Je2StEfSZ/L1E5ImJT2Wvy4rO2a9pP2S9kpa3u0vYGZm\n9alenb2k2RHxiqQTgG8DnwDeBbwUEZ+v2HcJsAm4CDgT2A6cExFHu3HxZmbWmLppnIh4JV+cBRwP\nvJi/V5XdVwKbI+JIRBwAngYu7sB1mplZG+oGe0nHSdoFTAE7IuLJfNNHJT0u6TZJY/m6M4DJssMn\nyVr4ZmbWR4207I9GxFLgLOD3JRWAW4BFwFLgp8DNtU7Rges0M7M2nNDojhHxC0n/CiyLiGJpvaQv\nA/fmb38MLCg77Kx83TEk+QfAzKwFEVEthV5XvWqc00opGkknAe8BHpM0r2y39wK78+WtwFWSZkla\nBCwGHpnhgv2KYMOGDX2/hkF5+V74Xvhe1H61o17Lfj6wUdJxZD8Md0bENyR9VdJSshTNs8CH8wC+\nR9IWYA/wKvCRaPcKzcysbTWDfUTsBi6ssv7PahzzaeDT7V+amZl1ikfQ9lmhUOj3JQwM34tpvhfT\nfC86o+6gqq58qOTsjplZkyQR3eigNTOz0eBgb2aWAAd7M7MEONibmSXAwd7MLAEO9mZmCXCwNzNL\ngIO9mVkCHOzNzBLgYG9mlgAHezOzBDjYm5klwMHezCwBDvZmZglwsDczS4CDvZlZAhzszcwS4GBv\nZpYAB3szswQ42JuZJcDB3swsATWDvaQ3SNopaZekPZI+k68/VdI2SfskPShprOyY9ZL2S9oraXm3\nv4CZmdWniKi9gzQ7Il6RdALwbeATwOXACxFxk6TrgTdFxDpJS4BNwEXAmcB24JyIOFpxzqj3uWZm\ndixJRIRaObZuGiciXskXZwHHAy+SBfuN+fqNwBX58kpgc0QciYgDwNPAxa1cmJmZdc4J9XaQdBzw\nA+AtwC0R8aSkuRExle8yBczNl88AHi47fJKshW9mZi0YH4d9+2D27PbOUzfY5ymYpZLeCDwg6Q8q\ntoekWjmZqtsmJiZeWy4UChQKhUau18wsGcVikQcfLPLcc+2fq27O/pidpb8G/hf4c6AQEQclzQd2\nRMS5ktYBRMSN+f73AxsiYmfFeZyzNzNrwIoVcN99sGwZPPpol3L2kk4rVdpIOgl4D/AYsBVYne+2\nGrgnX94KXCVplqRFwGLgkVYuzMzMYNMmWLUKtm1r7zw1W/aSzifrgD0uf90ZEZ+TdCqwBXgzcAC4\nMiIO5cfcAFwDvApcFxEPVDmvW/ZmlrTyXPymTTA2Vv+YdqpxmkrjdIqDvZmlrlCAb34zW161CrZs\nqX9MV0svzcys80rVNcuWwa23dv/z3LI3M+uDQ4eyVM6ttzaWwgGncczMktBOsK9bZ29mZv3Ts0FV\nZmbWulaqbsrt2zfdkdsOd9CamXVRKVjfd18W+JtV3pHbDgd7M7MuarfqpieDqrrFHbRmlopWqm5m\n4mocM7MEeFCVmZnV5GBvZpYAB3szswQ42JuZ5cbHswnKVqzIOlZHiYO9mVmu0Zr4YfxRcLA3s+TM\nFKwbrYlvd6BUPzjYm1lyZgrW5QOYatXE93p64k5wnb2ZJaf8ua71Ans1nRwo1QwPqjIza0K/gnW7\nHOzNzBLgEbRmZlaTg72ZWQIc7M3MEuBgb2aWgJrBXtICSTskPSnph5I+lq+fkDQp6bH8dVnZMesl\n7Ze0V9Lybn8BM0tTJ0axDuNI2FbVrMaRNA+YFxG7JJ0CfB+4ArgSeCkiPl+x/xJgE3ARcCawHTgn\nIo5W7OdqHDNrS6Ew/WzWVatgy5b+nKOXulaNExEHI2JXvvwy8BRZEAeo9oErgc0RcSQiDgBPAxe3\ncmFmZrV0YhTrMI6EbVXDOXtJC4ELgIfzVR+V9Lik2ySVhiWcAUyWHTbJ9I+DmVnHNDq1QbfPMSxO\naGSnPIVzN3BdRLws6Rbgk/nmTwE3A9fOcHjVfM3ExMRry4VCgUKh0NgVm5kBa9fC88/DBz6QBe1W\ngvXY2GCnborFIsVisSPnqjuCVtKJwNeB+yLiC1W2LwTujYjzJa0DiIgb8233AxsiYmfFMc7Zm1lb\nhi3f3gldy9lLEnAbsKc80EuaX7bbe4Hd+fJW4CpJsyQtAhYDj7RyYWY2Wjpd+ZJSvr0T6lXjvAN4\nCHiC6XTMDcD7gaX5umeBD0fEVH7MDcA1wKtkaZ8HqpzXLXuzxHSqJT4+nk1RfOKJcPLJ8JWvjH6+\nvcQToZnZwCoF5yefhBdeaH1a4ZIU0zcl7QT7hjpozcxaVXpQCMBZZ7Vf+eL0TWs8XYKZdVV5cN69\nu/2US0rlkp3kNI6ZdVUnHhRSSgXNnt16meUocM7ezPqq28E45Tx9OT+8xMz6aqYHeHdKK3n6lCY5\na4Q7aM2sbZ3oNK3118GmTc2ngso7hsfHs+NSTgW5ZW9mbetEp2mtvw5K0xo0c+7KH6Bu//Ux6Bzs\nzaxtrQTjSp0uqaz8AUq9ZNMdtGbWca102Haiaqef5+8FD6oys4FSmS+vVz1T/uPQLYM+w2W3OY1j\nZh01Pg7f/W62PGcOfO5z9Y9JPZ/eCw72ZlZTsyWM+/bBr36VLR8+DGvW1D8m9Xx6LzjYm1lNM7W6\nZ/oRKE/FLF3aWPD2FAjd5w5aM6tpxYos0FfOVjnTqNZDh+DqqyEiremHe8HTJZhZ18xUxTLTj4B1\nj4O9mfVMyg8P6TfPjWNmPVPK4W/fDrNmOdAPCwd7s0S1OlFYNypnPGlZ9znYmyWq1dr2blTOuM6+\n+zyC1ixRrbbQuzES1XX23ecOWrMEVJurZpDmihmkaxlknhvHzICZJyCrNbf7IEh93ppecM7ebITM\nlPtuZG53d5KOtprBXtICSTskPSnph5I+lq8/VdI2SfskPShprOyY9ZL2S9oraXm3v4CZTZsp993I\n3O7uJB1tNXP2kuYB8yJil6RTgO8DVwBXAy9ExE2SrgfeFBHrJC0BNgEXAWcC24FzIuJoxXmdszfr\ngkZz39X2qxwRu3Zt2o/xG0Q9G0Er6R7gi/nrkoiYyn8QihFxrqT1wNGI+Gy+//3AREQ8XHEeB3uz\nBrTyEJBWVf4AzDT3jfVPT0bQSloIXADsBOZGxFS+aQqYmy+fAUyWHTZJ1sI3sxZ0OrVSKy9f+WhB\nl0OOloaqcfIUzteA6yLiJWn6hyUiQlKtZnrVbRMTE68tFwoFCoVCI5dilpROB9xmniC1aZPLIfut\nWCxSLBY7cq66aRxJJwJfB+6LiC/k6/YChYg4KGk+sCNP46wDiIgb8/3uBzZExM6KczqNY9aATtef\ne6bK4da1nL2yJvxG4GcR8fGy9Tfl6z6bB/ixig7ai5nuoP2NysjuYG/WHx68NNy6GezfATwEPMF0\nOmY98AiwBXgzcAC4MiIO5cfcAFwDvEqW9nmgynkd7M36oJcdvtZ5ns/eLDGtBm1X2Aw3z2dvlphW\nq3RcYZMuB3uzIVQK2qecAi++2Pj0Bn6wd7qcxjEbQocOweLF8MIL2XunZNLgNI5ZYsbG4KKLsmWn\nZKwRbtmbDSmXUabH1Thm1hSXYA4np3HMrCmezjg9DvZmPTJIDwdxCWZ6HOzNemSQWtMuwUyPn0Fr\n1iH18uCD1Jr2M1/T45a9WYfUa7l3szU9SCkiG0wO9mYdUm9Ua+XDQTppkFJENpgc7C1p9VrEzbSY\nN22C00+Hl1+G7dt7G3QHKUVkg8nB3pJWr0XcTIt5bCwLttD7oOsOV6vHwd6SVq9F3GyLuV9Bt5sp\nIhsNHkFrSas35YCnJLBB4ukSzMwS4OkSzMysJgd7M7MEONibmSXAwd6S1a9Rpx7tav3gYG/J6teo\nU492tX5wsLdk9WvUqUe7Wj/UDfaSbpc0JWl32boJSZOSHstfl5VtWy9pv6S9kpZ368LN2tWvAVAe\n7Wr9ULfOXtI7gZeBr0bE+fm6DcBLEfH5in2XAJuAi4Azge3AORFxtGI/19mbmTWpq3X2EfEt4MVq\nn1tl3Upgc0QciYgDwNPAxa1cmNkgcueqDat2cvYflfS4pNsklf4YPQOYLNtnkqyFb9Y1nQrAjZzH\nnas2rFp9UtUtwCfz5U8BNwPXzrBv1XzNxMTEa8uFQoFCodDipVjqSgEYsgDc6hOYGjmPO1etl4rF\nIsVisSPnamhuHEkLgXtLOfuZtklaBxARN+bb7gc2RMTOimOcs7eOWbEia2kvW9Zep2cj5/HEaNZP\nPZ8bR9L8srfvBUqVOluBqyTNkrQIWAw80spnmDVifBwOH4Z58+Duu9sLwI1UyXgqYRtWjVTjbAYu\nAU4DpoANQAFYSpaieRb4cERM5fvfAFwDvApcFxEPVDmnW/bWEYXCdOpl1So/RNtGm6c4tmTVSr2M\nj2d5+Nmzs1a7W+M27DzFsSWrVurFlTNm09yyt5G1YAFMTsKcOfDEE3D22f2+IrP2uGVvVkUpuB8+\nDGvW9PdazPrNwd5G1pw52X9dE2/mNI6NMNfE26hxNY6NPFfWmDlnbwlwZY1ZexzsbSh4Thqz9jiN\nY0PB+Xcz5+xtyDj/btYa5+xtqDj/btZ7rc5nb4lppzVeeWyn8u/+C8GscW7ZW0PaaY1XHls5n02r\nT5ryXwhmjXPL3hrSTmu88tjSnPAlrT5pyhU6Zo1zy94a0siDPVo9ttWg3c41maXG1TjWdy6rNGuM\nSy/NzBLg0ksbSq12zJpZ8xzsrW9cTWPWOw721jeupjHrHefsrW/cMWvWHHfQWt+URrE+80z2GMA5\nczya1axbHOytb+bPh4MHj123atX0wChPaWDWOV2txpF0u6QpSbvL1p0qaZukfZIelDRWtm29pP2S\n9kpa3spF2fD45S+PfV+Zf7/33ulO2A9+sKeXZmZlGumgvQO4tGLdOmBbRJwDfCN/j6QlwPuAJfkx\nX5LkTuAR9ta3Zv89/3xYufL1o1nLfwzUUnvEzDqhbiCOiG8BL1asvhzYmC9vBK7Il1cCmyPiSEQc\nAJ4GLu7MpVovNFv7ftddWdrmoYfgnnten6Yp/RhccAHccUfHL9fMGtTqRGhzI2IqX54C5ubLZwAP\nl+03CZzZ4mdYHzQ6KVmjufi77nLFjdkgaHvWy4gISbV6W6tum5iYeG25UChQKBTavRTrgEZr3xv9\nUaic4dLMGlcsFikWix05V0PVOJIWAvdGxPn5+71AISIOSpoP7IiIcyWtA4iIG/P97gc2RMTOivO5\nGmdAlde+r107c+t9xYqs03XZMs86adYr/ZgbZyuwOl9eDdxTtv4qSbMkLQIWA4+0+BnWB6WW+NhY\n7ekMPL2w2XCpm8aRtBm4BDhN0o+AvwFuBLZIuhY4AFwJEBF7JG0B9gCvAh9xE3541UrpOD1jNlw8\nqGpI9WKwkqczMBssHkGboEJhuoO0fMSqmY0uz2efIM8YaWbNcMt+SDWSYvG8NGajxWkcq8qpHrPR\n4jSOVeVUj5mVuGU/IqqlbFxNYzZanMaxllI2zumbDRencayllI0f+G2WDgf7EdHK9AXO6Zulw2mc\nhDmnbzZcnLMfAc6fm1k9ztmPgGr582afGlWunWPNbPS0/fAS64xq+fPyB4Scdx489VTjLf5GHy5i\nZmlwy35AVOtgLf0AABw8eGzFTL2Wuztfzaycc/YD7NChrEV/8ODrnwhVr67ena9mo8c5+xE1Npal\nbqqVVNZruZc/ccrMzC37IeWWu1l6XHo5YFxGaWbd4DTOgPE0BGY2aBzsu8CVMGY2aJzG6QLn082s\nG5yzHwLO45tZu5yzHwLO45tZP7U1XYKkA8Bh4P+AIxFxsaRTgX8CzgYOAFdGxEjNzlJqpT/zDJx9\nNsyZk7XW166dufXuPL6Z9VNbaRxJzwJvjYifl627CXghIm6SdD3wpohYV3HcUKdxykevlqxaBc8/\nP/OoVufxzaxd/U7jVH7w5cDGfHkjcEUHPqPrmpklstRKnzMn+2+ptV6r9e4RrWbWT+0G+wC2S3pU\n0ofydXMjYipfngLmtvkZPdFMTr00adkTTxw7lUErT4syM+uFdqc4/r2I+Kmk04FtkvaWb4yIkFQ1\nXzMxMfHacqFQoFAotHkp7Wkmpz42lr1Wrz52ZspS693MrBOKxSLFYrEj5+pY6aWkDcDLwIeAQkQc\nlDQf2BER51bsO3A5+2Zz6vVmnTQz67S+5OwlzZb06/nyycByYDewFVid77YauKfVz+imyhx9szl1\nV9eY2TBpJ40zF/gXSaXz/ENEPCjpUWCLpGvJSy/bvsouaOVJTuUDo265BdascXWNmQ2HloN9RDwL\nLK2y/ufAu9u5qF5opWVe/gOxZo1TN2Y2PJIdQVtZOdNI6aVTN2Y2rDw3Tq5eh+v4OOzZk42affjh\nbOSsmVkv9XtQ1UCp1kKvXFdtn/JW+0knvX77vn3wne9kz4Nds6a338nMrF1DEeybGd1abXBU5bpq\n+5SndZ577vXbncIxs2E2FMG+mdGt1YJy5bpq+5SXXlbb7tGxZjbMhiLYN9qqHh+Hw4dh3jy4++7p\noFwZqE8/HU47beagXS2we24bMxtmQ9FB2+jo1kZHtXr0q5kNo3Y6aNudG6cnGp1zptG/AJx/N7PU\nDEXLvlGN/gXgueXNbBgl8wxaP8fVzFKWTJ29n+NqZtaaoQr2zrWbmbVmqNI4zrWbWcqSydmbmaVs\n6HP2zUyHYGZmzRuIYO+OVzOz7hqIYO+OVzOz7hqInL07Xs3M6hvpDloPpDIzywx9B20t1fL57tA1\nM2vOQAX7ek+QKuXz3aFrZtacgZr1shTEIQviY2PV56d3h66ZWXO60rKXdKmkvZL2S7q+0eMqg/hM\nz331U6PMzJrT8WAv6Xjgi8ClwBLg/ZLOa+TYyiA+Uwt+lJ4aVSwW+30JA8P3YprvxTTfi87oRsv+\nYuDpiDgQEUeAfwRWNnJgZRBPoQXv/5Gn+V5M872Y5nvRGd3I2Z8J/Kjs/STwtlZO1OgTqszMrLZu\ntOw9w5mZ2YDp+KAqSW8HJiLi0vz9euBoRHy2bB//IJiZtWBgRtBKOgH4T+BdwE+AR4D3R8RTHf0g\nMzNrWMdz9hHxqqS/BB4Ajgduc6A3M+uvvsyNY2ZmvdXT6RJaHWw1CiQtkLRD0pOSfijpY/n6UyVt\nk7RP0oOSRrTI9PUkHS/pMUn35u+TvBeSxiTdLekpSXskvS3he7E+/zeyW9ImSb+Wyr2QdLukKUm7\ny9bN+N3ze7U/j6nL652/Z8G+ncFWI+II8PGI+E3g7cBf5N9/HbAtIs4BvpG/T8V1wB6mK7hSvRd/\nB/xbRJwH/BawlwTvhaSFwIeACyPifLI08FWkcy/uIIuP5ap+d0lLgPeRxdJLgS9JqhnPe9myb3mw\n1SiIiIMRsStffhl4imxMwuXAxny3jcAV/bnC3pJ0FrAC+DJQqi5I7l5IeiPwzoi4HbI+r4j4BQne\nC+AwWaNodl7oMZusyCOJexER3wJerFg903dfCWyOiCMRcQB4mizGzqiXwb7aYKsze/j5AyNvwVwA\n7ATmRsRUvmkKmNuny+q1vwXWAEfL1qV4LxYB/y3pDkk/kPT3kk4mwXsRET8Hbgb+iyzIH4qIbSR4\nL8rM9N3PIIuhJXXjaS+DvXuCAUmnAF8DrouIl8q35U90Gfn7JOmPgecj4jGmW/XHSOVekFXEXQh8\nKSIuBP6HijRFKvdC0luAvwIWkgWzUyT9Sfk+qdyLahr47jXvSy+D/Y+BBWXvF3DsL9PIk3QiWaC/\nMyLuyVdPSZqXb58PPN+v6+uh3wUul/QssBn4Q0l3kua9mAQmI+J7+fu7yYL/wQTvxTLgPyLiZxHx\nKvDPwO+Q5r0omenfRGU8PStfN6NeBvtHgcWSFkqaRda5sLWHn99XkgTcBuyJiC+UbdoKrM6XVwP3\nVB47aiLihohYEBGLyDrg/j0i/pQ078VB4EeSzslXvRt4EriXxO4FWcf02yWdlP97eTdZB36K96Jk\npn8TW4GrJM2StAhYTDaAdWYR0bMXcBnZ6NqngfW9/Ox+v4B3kOWndwGP5a9LgVOB7cA+4EFgrN/X\n2uP7cgmwNV9O8l4Avw18D3icrDX7xoTvxVqyH7vdZB2SJ6ZyL8j+yv0J8Cuy/s2ra3134IY8lu4F\n/qje+T2oyswsAQP1DFozM+sOB3szswQ42JuZJcDB3swsAQ72ZmYJcLA3M0uAg72ZWQIc7M3MEvD/\ntylhZXB+OsoAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Let's now imagine that this is real data, and we want to determine the slope and intercept of the best-fit line to the data. We start off by definining a function representing the model:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def line(x, a, b):\n", " return a * x + b" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The arguments to the function should be ``x``, followed by the parameters. We can now call ``curve_fit`` to find the best-fit parameters using a least-squares fit:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "popt, pcov = curve_fit(line, x, y)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The ``curve_fit`` function returns two items, which we can ``popt`` and ``pcov``. The ``popt`` argument are the best-fit paramters for ``a`` and ``b``:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "popt" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "pyout", "prompt_number": 6, "text": [ "array([ 2.90060552, 8.48416222])" ] } ], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "which is close to the initial values of ``3`` and ``2`` used in the definition of ``y``.\n", "\n", "The reason the values are not exact is because there are only a limited number of random samples, so the best-fit slope is not going to be exactly those used in the definition of ``y``. The ``pcov`` variable contains the *covariance* matrix, which indicates the uncertainties and correlations between parameters. This is mostly useful when the data has uncertainties." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Let's now try and fit the data assuming each point has a vertical error (standard deviation) of +/-10:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "e = np.repeat(10., 100)\n", "plt.errorbar(x, y, yerr=e, fmt=\"none\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "pyout", "prompt_number": 7, "text": [ "" ] }, { "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGh9JREFUeJzt3X2wHXV9x/H3h6cqMmNgqEkIgcsoGY2ljVSibbU5aa1C\n7QBOZxA71uDTMGMHKKMtSf7ovVdn5KEDpdqRqQISGEEZaS1ogYDNiXU6BB9AKSEKUy4SJEGBIJQ/\nmpBv/9g93s3mnHP3nLPncT+vmTt3z57dPXs3N5/93t/+9reKCMzMrBoOGfYOmJnZ4Dj0zcwqxKFv\nZlYhDn0zswpx6JuZVYhD38ysQgqFvqRXSdom6UFJ2yVdms6fkbRT0gPp1xmZdTZIelTSDknv7tcP\nYGZmxaloP31JR0bEy5IOA74LfAr4Y+DFiLgqt+xK4GbgNGAZcC+wIiL2l7nzZmbWmcLNOxHxcjp5\nBHAo8Hz6Wk0WPwu4JSL2RsQc8Biwuof9NDOzEhQOfUmHSHoQ2A1siYiH07cukPQjSddJWpTOOw7Y\nmVl9J0nFb2ZmQ9RJpb8/IlYBxwN/KKkGXAOcBKwCngaubLeJHvbTzMxKcFinK0TEC5K+Bbw1IuqN\n+ZKuBe5IXz4FLM+sdnw6j8zyPgmYmXUhIpo1qxdStPfOsY2mG0mvBv4EeEDSksxi7wMeSqdvB86V\ndISkk4CTgfub7Li/Ipienh76PozKl4+Fj4WPRfuvXhWt9JcCmyQdQnKiuCkivi3pRkmrSJpuHgfO\nT8N8u6Rbge3APuATUcbemplZTwqFfkQ8BJzaZP6H2qzzWeCz3e+amZmVzXfkjoBarTbsXRgZPhbz\nfCzm+ViUp/DNWaV/sOQWHzOzDkki+n0h18zMJoND38ysQhz6ZmYV4tA3M6sQh76ZWYU49M3MKsSh\nb2ZWIQ59M7MKceibmVWIQ9/MrEIc+mZmFeLQNzOrEIe+mVmFOPTNzCrEoW9mViEdPxjdzMwGp15P\nvhrTvfJDVMzMxoQEMICHqEh6laRtkh6UtF3Spen8YyTdI+mnkjZLWpRZZ4OkRyXtkPTubnfQzMzK\nU7jSl3RkRLws6TDgu8CngDOBX0bEFZIuAY6OiPWSVgI3A6cBy4B7gRURsT+zPVf6ZmYdKKPSL9ym\nHxEvp5NHAIcCz5OE/pp0/iagDqwHzgJuiYi9wJykx4DVwH3d7qiZ2STJt9U3nv1eq81P90Ph0Jd0\nCPBD4PXANRHxsKTFEbE7XWQ3sDidPo4DA34nScVvZmYcGO5SORdpi+ik0t8PrJL0WuBuSWtz74ek\ndu01bssxMxuyjrtsRsQLkr4F/C6wW9KSiNglaSnwTLrYU8DyzGrHp/MOMDMz8+vpWq1GrZ9/05iZ\njaF6vU69xD8DCl3IlXQssC8i9kh6NXA3MAu8B3g2Ii6XtB5YlLuQu5r5C7lvyF659YVcM7OEBK3i\nMNv2PzsLvV7ILRr6p5BcqD0k/bopIv5e0jHArcAJwBxwTkTsSdfZCHwE2AdcFBF357bp0DeziVfk\ngm270M8qo/eOb84yMxuQVuE+yND32DtmZhXisXfMzEZYtnmoDG7eMTMbEDfvmJnZQLnSNzMbkGxF\n380wDO69Y2ZWkkGMhVO0Gafd+g59M7OSFb1ZqtOTg0PfoW9mA1b2zVKdxJhD36FvZkNURm+afod+\n/iS1datD38ysK4MI/bKvFUgOfTOzrgy60i9Dr6HvfvpmZhXi0DczqxCHvplZhTj0zcwqxKNsmtlY\n6rVXTGPdmZmD159k7r1jZmOv2140jfXy609y7x2HvpmNvTJCf8uWYn85DGKMnnYc+mZWeZ2EfrPQ\nnp1NQn8cmnYGEvqSlgM3Aq8DAvhiRHxO0gzwMeAX6aIbI+LOdJ0NJA9GfwW4MCI257bp0DezrmXD\ne3YWpqeT6W4q7mE003RrUKG/BFgSEQ9KOgr4AXA2cA7wYkRclVt+JXAzcBqwDLgXWBER+zPLOPTN\nKqSfzSJlDGQ2LnHUa+gX6r0TEbuAXen0S5IeIQlzgGYffhZwS0TsBeYkPQasBu7rdkfNbLzlR7Hs\n9rmvrU4e2WlrreM2fUlTwFbgzcAngQ8DLwDfBz4ZEXskfR64LyK+kq5zLXBnRNyW2Y4rfbMJ1yqg\nZ2fLqaxb9b7pdjvjYCCVfubDjgK+DlyUVvzXAJ9O3/4McCXw0RarH3RIZ2Zmfj1dq9Wo+TRtNlFa\nVfezs0PaoTFUr9epd/tnUROFK31JhwPfJKnYr27y/hRwR0ScImk9QERclr53FzAdEdsyy7vSN6uQ\nbDXdaWW90F8MrvQ7WL/ghVwBm4BnI+LizPylEfF0On0xcFpE/EXmQu5q5i/kviGb8g59s8nR6dOo\negnZZttx6HewfsHQfwfwHeDHzDfTbAQ+AKxK5z0OnB8Ru9N1NpJ02dxH0hx0d26bDn2zCVRkjPoy\nQr9eh7Vrk66a+S6bsPBJaNg3WXXLN2eZ2UjpNvSLhnCnlX72JDGOIZ/n0DezkVJGpV8kxLPTRZcv\nsv1RN9DeO2ZmRbgv/ehypW9mpWo1YmW7Sj9/kti6NWmjb9b00qxNv2hzULt548LNO2Y2MhpBvG5d\nMj01lQT4unWwaVPyfWpq4ZCGzpp32nHo59Z36JtZO51eAM1X9p12qywS+o2eOp3sz0LzxoVD38wG\nppPKujFGfZEmmPz60NmF2U732aE/BA59s/GTDcuid8l2E9Lg0G/FoW9mA9NJd0yHfn+4y6aZDVSr\nCv+885KLtGvWJA8bb7xf9mdmtz9uN1aNAlf6ZlZYke6YRd7Lyz8FC1p32SxioYvPVa70HfpmFdTt\nkAT9Cv3sPs3OJtV82cMkeBiGdH2Hvlm1dVL19jP0u9mfKnKbvpktqF2Va9Xi0DebIO3CvXHxU0qm\ns8uWfWE0v92Ftp/dl37sj81z847ZhCrSvTL7utMHoWRfu3lncNymb2ZNdRr6vaybP3HMzc2PvTM3\nl5w0pqaKVe4O/fYc+mbWVK+hv2XLwgG+dm05jyssst+WcOibTaheuxiWWem3etJVp+PqFOHQb8+h\nb1YB3QRh0aCenU2q+nZj0Bc9gZTBod/eoB6Mvhy4EXgdyUPQvxgRn5N0DPA14ERgDjgnIvak62wg\neTD6K8CFEbE5t02HvllBZYV+dn6Ri7DNXvf7JieHfnuDCv0lwJKIeFDSUcAPgLOBDwO/jIgrJF0C\nHB0R6yWtBG4GTgOWAfcCKyJif2abDn2zgkYp9PvNod/eQG7OiohdwK50+iVJj5CE+ZnAmnSxTUAd\nWA+cBdwSEXuBOUmPAauB+7rdUTPrXaNCP++85Ht+sLRh9Yv3gGqD03GbvqQpYCvwW8DPIuLodL6A\n5yLiaEmfB+6LiK+k710L3BkRt2W240rfrKB89dtNn/r8tka10rf2BjoMQ9q0cxtwUUS8KM1/bkSE\npHa/Fge9N9M4nQO1Wo2aT+lmheTDvawhjG301Ot16iX+Axeu9CUdDnyTpGK/Op23A6hFxC5JS4Et\nEfFGSesBIuKydLm7gOmI2JbZnit9s4LaVfqzs/PPjF20CPbsmV+mcWLIz5+amn9QeaPvPbQeetiV\n/ugY1IVckbTZPxsRF2fmX5HOuzwN+kW5C7mrmb+Q+4Zsyjv0zYprF7rdBnIn23Toj45Bhf47gO8A\nP2a+mWYDcD9wK3ACB3fZ3EjSZXMfSXPQ3bltOvTNChpE6Le7TtC489aGzzdnmVVA2aFf5G7aSXno\nyKRx6JtNqKKh22ul76ab8eKHqJiNiU4r535U1K3Grc/uj002V/pmQzAKo1K60h9PrvTNRoTbwG0c\nOPTNSpK9ALp16/y8xomgn8Ff9ITTargDN+9Uh5t3rNIWemB4t5V742b1srtZNrS6Oaub7pVu3hkv\n7r1jVpIyu0X2O/TbbcvPpJ1sbtM3q6B2o1KateNK3yw1rpV+r9t2pT9eeq30DylzZ8zMbLS50jdL\nVanSd/fS8eULuWYlqVLo2/jyhVyzLuWrXUguiPar2m183twcnHhi8hmNseyH+ahCqxZX+maUX5kv\ntL1ut1v0s/1fa3K50jcbEflulFu39vcvB7NuOPTNSpK9izc/9EH2/X5o12/fJxzLcvOOVVojLGdn\nk7AsaxiGotwUY51y7x2zEkiwZcvCIV92V0eHvnXKoW+VV0YQdzN+TRmB7dC3Tg3qwejXA+8FnomI\nU9J5M8DHgF+ki22MiDvT9zaQPBT9FeDCiNjcZJsOfStdWQ8Jd+jbqBpU750vA58HbszMC+CqiLgq\nt0MrgfcDK4FlwL2SVkTE/m530qyVQfe1Nxt3hZt3JE0Bd2Qq/WngpYi4MrfcBmB/RFyevr4LmImI\n+3LLudK3UhXtG99qXVf6Ng6G3U//AkkfAr4PfDIi9gDHAdmA30lS8ZsNlMeXMTtYL6F/DfDpdPoz\nwJXAR1ss27SWmWl0JgZqtRo1/0+0Ptm69eCgH9aY9O5Tb52o1+vUG78wJei6eafVe5LWA0TEZel7\ndwHTEbEtt46bd6xU7Zp3+jHGvJtmbBiGNp6+pKWZl+8DHkqnbwfOlXSEpJOAk4H7u/0cMzMrT6Hm\nHUm3AGuAYyU9CUwDNUmrSJpuHgfOB4iI7ZJuBbYD+4BPuKQ3MxsNvjnLxkKRi7Ju3rEq8B25Vjmt\nwrbX0C9yYnGPIBs2h75VTj7A2w2alv0rwL9uNgkc+lY57Sr9YTye0GyQHPo2dnptIskGeNFtOfRt\nUjj0baQtFMrdhHHRYZDz6/jXzSaBQ98Goh/DFzfm9RrgrQI9u8+zszA93fk+m40ah74NXFnDFzeb\n101XySKh7542Nikc+jZw4xL6ZpNoaMMwmJnZ+HHom5lVSK/j6Zt1JdveXqSNfVjDIJtNGrfpW8fK\nbtOH+fndduH0r5JVhdv0zcysMDfvWE8azS5zc8n3qalkulZLpt090my0OPStFFNT8MQTcN55yY1Q\n552XhH29Pt/+Dsn03Nz8OmvWJI8ynJnxycFsENymbx1baMCzxvd6HW64IQn5ubnkpLBmTRL2jZNC\nfntF2+d945VVlW/OslK1CtNFi2DPnoPnNxu+uMgF2/x7jSESHOBm7Tn0rW+6fdpUN6HvXwWzYnoN\nfbfpW8da/TVgZqOvUKUv6XrgvcAzEXFKOu8Y4GvAicAccE5E7Enf2wB8BHgFuDAiNjfZpiv9EVav\nw9q1SbNL0XHqXemb9d9AmnckvRN4CbgxE/pXAL+MiCskXQIcHRHrJa0EbgZOA5YB9wIrImJ/bpsO\n/RHXLsTzyyy0vEPfrBwDa9OXNAXckQn9HcCaiNgtaQlQj4g3plX+/oi4PF3uLmAmIu7Lbc+hP+L6\nGfrufWPWnWG26S+OiN3p9G5gcTp9HJAN+J0kFb+VaNxDc1z202zSlHIhNyJCUruyvel7M5m7dmq1\nGjWnQGH5rpKNE4CZTZZ6vU69xP/gvTbv1CJil6SlwJa0eWc9QERcli53FzAdEdty23PzTkn61SZe\ntHlnejq5+WrTpuTmq2bDMLRr0zez4obZvHM7sA64PP3+jcz8myVdRdKsczJwfw+fY10aZBNQY0iF\n/Gdkh2HID4nsP+zMBq9o751bgDXAsSTt938H/BtwK3ACB3fZ3EjSZXMfcFFE3N1km670S9LNIwY7\n2W7RC7l5437dwWwU+Y5cKzSMQS9j2szOwpYtzYPa3S3NBsuhb4UGLGs2v2gl3m2lb2blc+hb16Ff\n9P38e262MRseh771JfQd7GajyaFvfa/0zWx0eJTNCVF2ZZ3fHsw/ncqVull1udIfQVLSWyb71ClI\n+sI3njoFzU8Ss7OdDXaWXcb/HGajz5X+hMp3t4T58AcPu2Bm3XGlP4LyVfdClfpCbfqu9M0mhyt9\naynbru8hEMwMXOmPpH5U+v0YpsHMBs9dNkdQrz1xugn9dsMwOPTNJodDf8T1MtBZ9jV0PxRCq/d9\nA5bZ+HHoj7hRDn0zGz++kDsmXFWb2Shwpd9n3VxYdaVvZq240h+CRtU+N5d8n5qCHTtgyRJYtAj2\n7IGzz+5um3Bw98pe9rHZ9vyXhVl1udJPddv80uzpUmV0oazX2w/DUHS8ezObLL6Q2wedhGi/Qr9f\n+2tm483NO0OS/cugUXU3RrEc9Oe7+cbMiuq50pc0B/wKeAXYGxGrJR0DfA04kdxD0zPrTUyl3zCs\nSt/MqqPXSv+QEvYhgFpEvCUiVqfz1gP3RMQK4NvpazMzG7KymnfyZ50zgTXp9CagzogFf7sLt52u\n00yjucVNL2Y2Sspo3vkf4AWS5p1/jogvSXo+Io5O3xfwXON1Zr2Rad7Jj12zdWvyul1IZy/gNiz0\nyML8umZmnRqFC7l/EBFPS/pN4B5JO7JvRkRIahpxM40SGKjVatSGWAY3dqUR4pldMzMbmnq9Tr3E\npyaV2mVT0jTwEvBxknb+XZKWAlsi4o25ZYdW6eebabKV/dq1yfzsrjVr1pmdTR5p2Fg+u06zSt7D\nMJhZGYbaT1/SkcChEfGipNcAm4FZ4F3AsxFxuaT1wKKIWJ9bt1DodxKWzZbN3tjUathhODCwITkJ\ntBumuJfmHTOzbg079E8C/jV9eRjwlYi4NO2yeStwAiV22ezmpqki8+Dg0M/fcJVdvlHhr1uXnBie\neAIWL24+DIMreTMr01Db9CPicWBVk/nPkVT7A5Wt9KF1r5n8cpAMbdCQ7XHTrJdO4/UNNyTfJdi1\nq6tdNjMbqLEahmFQlX6RG6sWWs7MrB8qNfbOqIV+u0cUmpn1Q+VCv2jQDrLSNzMblMqFfi83QDn0\nzWzcjcLNWWOp1SiVC60DybIeXsHMxpEr/SaPJvQImWY2qlzpl6Cbqt/MbByNRKVf9K7bsir97OfN\nzs5fHG71eR5CwcxGxcRdyM2HdavAnZ1tP75NI8yLDMNQ5PMc8GY2CiY+9FuvX87wxW6jN7NxMvFt\n+p08uMTMzNobq0q/MdhZvZ4036xZc/ATr4o0zbgJx8zG1cQ07zSCuFmYN6bzF3LBTTNmVi0TE/rz\n85PvzXrNbN2anBCmppILtEUea2hmNkkmOvQXmu+LsGZWNb2G/iFl7oyZmY02h76ZWYU49M3MKqRv\noS/pdEk7JD0q6ZJ2y9bryXg3jdErIZnOP9LQzMx605cLuZIOBX5C8pzcp4DvAR+IiEcyy7QcZbPI\nuPetljMzm2SjeiF3NfBYRMxFxF7gq8BZffosMzMrqF/DMCwDnsy83gm8rZsNtRr22P3yzcw616/m\nnT8HTo+Ij6evPwi8LSIuyCzTcfOOH2hiZlU3qgOuPQUsz7xeTlLtH2Am87SSRYtq7NlTA9o/irBd\n5e/q38wmTb1ep15ir5Z+VfqHkVzI/WPg58D9dHAht/k2XdWbmY3sMAySzgCuBg4FrouIS3PvLxj6\nHg3TzOxAIxv6C35wFw9GNzOrulHtsmlmZiPIoW9mViEOfTOzCnHom5lViEPfzKxCHPpmZhXi0Dcz\nqxCHvplZhTj0zcwqxKFvZlYhDn0zswpx6JuZVYhD38ysQhz6ZmYV4tA3M6sQh76ZWYU49M3MKsSh\nb2ZWIQ59M7MK6Tr0Jc1I2inpgfTrjMx7GyQ9KmmHpHeXs6tmZtarXir9AK6KiLekX3cCSFoJvB9Y\nCZwOfEGS/6Joo16vD3sXRoaPxTwfi3k+FuXpNYybPZH9LOCWiNgbEXPAY8DqHj9novkXep6PxTwf\ni3k+FuXpNfQvkPQjSddJWpTOOw7YmVlmJ7Csx88xM7MStA19SfdIeqjJ15nANcBJwCrgaeDKNpuK\n8nbZzMy6pYje81jSFHBHRJwiaT1ARFyWvncXMB0R23Lr+ERgZtaFiGjWtF7IYd2uKGlpRDydvnwf\n8FA6fTtws6SrSJp1Tgbuz6/fy06bmVl3ug594HJJq0iabh4HzgeIiO2SbgW2A/uAT0QZf06YmVnP\nSmneMTOz8TCU/vOSTk9v3HpU0iXD2IdhkbRc0hZJD0v6b0kXpvOPSS+c/1TS5kxvqIkn6dD0Br87\n0teVPBaSFkn6uqRHJG2X9LYKH4sN6f+RhyTdLOk3qnIsJF0vabekhzLzWv7snd4MO/DQl3Qo8E8k\nN26tBD4g6U2D3o8h2gtcHBFvBt4O/FX6868H7omIFcC309dVcRFJc2Djz86qHot/BP49It4E/Daw\ngwoei7RjyMeBUyPiFOBQ4Fyqcyy+TJKPWU1/9m5uhh1Gpb8aeCwi5iJiL/BVkhu6KiEidkXEg+n0\nS8AjJBe8zwQ2pYttAs4ezh4OlqTjgT8FrmX+Zr/KHQtJrwXeGRHXA0TEvoh4gQoeC+BXJMXRkZIO\nA44Efk5FjkVE/CfwfG52q5+945thhxH6y4AnM68re/NWWtG8BdgGLI6I3elbu4HFQ9qtQfsH4G+A\n/Zl5VTwWJwG/kPRlST+U9CVJr6GCxyIiniO57+dnJGG/JyLuoYLHIqPVz97xzbDDCH1fOQYkHQXc\nBlwUES9m30t7O038cZL0Z8AzEfEAzYf0qMyxIOlJdyrwhYg4Ffhfcs0XVTkWkl4P/DUwRRJqR0n6\nYHaZqhyLZgr87G2PyzBC/ylgeeb1cg48U008SYeTBP5NEfGNdPZuSUvS95cCzwxr/wbo94EzJT0O\n3AL8kaSbqOax2AnsjIjvpa+/TnIS2FXBY/FW4L8i4tmI2Af8C/B7VPNYNLT6P5HP0+PTeS0NI/S/\nD5wsaUrSESQXIW4fwn4MhSQB1wHbI+LqzFu3A+vS6XXAN/LrTpqI2BgRyyPiJJILdf8REX9JNY/F\nLuBJSSvSWe8CHgbuoGLHguQC9tslvTr9//Iukgv9VTwWDa3+T9wOnCvpCEkn0eJm2ANExMC/gDOA\nn5BcdNgwjH0Y1hfwDpL26weBB9Kv04FjgHuBnwKbgUXD3tcBH5c1wO3pdCWPBfA7wPeAH5FUt6+t\n8LH4W5KT3kMkFy4Pr8qxIPmr9+fA/5Fc//xwu58d2Jhm6Q7gPQtt3zdnmZlViB9uYmZWIQ59M7MK\nceibmVWIQ9/MrEIc+mZmFeLQNzOrEIe+mVmFOPTNzCrk/wFarZ3EYIos6QAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 7 }, { "cell_type": "code", "collapsed": false, "input": [ "popt, pcov = curve_fit(line, x, y, sigma=e)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "prompt_number": 8 }, { "cell_type": "code", "collapsed": false, "input": [ "popt" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "pyout", "prompt_number": 9, "text": [ "array([ 2.90060551, 8.48416251])" ] } ], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Now ``pcov`` will contain the true variance and covariance of the parameters, so that the best-fit parameters are:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print(\"a =\", popt[0], \"+/-\", pcov[0,0]**0.5)\n", "print(\"b =\", popt[1], \"+/-\", pcov[1,1]**0.5)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "a = 2.90060551367 +/- 0.0330874141721\n", "b = 8.48416250909 +/- 1.804681477\n" ] } ], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "We can now plot the best-fit line:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.errorbar(x, y, yerr=e, fmt=\"none\")\n", "xfine = np.linspace(0., 100., 100) # define values to plot the function for\n", "plt.plot(xfine, line(xfine, popt[0], popt[1]), 'r-')" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "pyout", "prompt_number": 11, "text": [ "[]" ] }, { "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOXZ//HPxabiBhZFRDRWoEpdcEOtVOLyVK0WrK2o\nrYpLbX30cUGU3QfShlVBFLdaN1xAcQd34MegqIj6iKgIghIUFXABBRVluX5/nBkySWYmM5PJTJLz\nfb9eeWVy5ix3jvidK/e5z33M3RERkXBoVOgGiIhI/ij0RURCRKEvIhIiCn0RkRBR6IuIhIhCX0Qk\nRNIKfTPb2sxeN7N5ZrbAzEZElw81s+Vm9nb066S4bQaY2WIzW2hmv6utX0BERNJn6Y7TN7Pm7v6D\nmTUBZgNXA8cBa919bKV1OwETgcOAtsB0oKO7b85l40VEJDNpd++4+w/Rl82AxsDq6M+WYPUewCR3\n3+DuZcASoEsN2ikiIjmQduibWSMzmwesBGa6+/vRty4zs3fM7C4zaxFdthuwPG7z5QQVv4iIFFAm\nlf5md+8M7A4cbWbFwG3AXkBn4AtgTKpd1KCdIiKSA00y3cDdvzWzZ4BD3T0SW25mdwJToz9+BrSL\n22z36DLi1teHgIhIFtw9Ubd6WtIdvdMq1nVjZtsA/wW8bWa7xq32R+Dd6OspwJlm1szM9gI6AHMT\nNFxf7gwZMqTgbagrXzoXOhc6F6m/airdSr8NMMHMGhF8UNzv7jPM7D4z60zQdbMU+Ec0zBeY2WRg\nAbARuMRz0VoREamRtELf3d8FDk6w/NwU2wwHhmffNBERyTXdkVsHFBcXF7oJdYbORTmdi3I6F7mT\n9s1ZOT+wmXp8REQyZGZ4bV/IFRGRhkGhLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJEYW+\niEiIKPRFREJEoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCZGMH4wuIiL5E4kEX7HXNaWH\nqIiI1BNmAHl4iIqZbW1mr5vZPDNbYGYjost3MrNpZvahmb1oZi3ithlgZovNbKGZ/S7bBoqISO6k\nXembWXN3/8HMmgCzgauB7sBX7j7azPoBLd29v5l1AiYChwFtgelAR3ffHLc/VfoiIhnIRaWfdp++\nu/8QfdkMaAysJgj9btHlE4AI0B/oAUxy9w1AmZktAboAc7JtqIhIQ1K5rz727Pfi4vLXtSHt0Dez\nRsD/AXsDt7n7+2bW2t1XRldZCbSOvt6NigG/nKDiFxERKoa7WW4u0qYj7SGb7r7Z3TsDuwNHm9kx\nld53IFV/jfpyRERqoDEba7yPjIdsuvu3ZvYMcAiw0sx2dfcVZtYGWBVd7TOgXdxmu0eXVTB06NAt\nr4uLiymuzb9pRETqoUgkQmTGDHjnHS7hZcbXcH9pXcg1s1bARndfY2bbAC8AJcAJwNfuPsrM+gMt\nKl3I7UL5hdz28VdudSFXRCRgBgnj8McfWdzvTnaZcB1f7tyJ8z8axGyOzsuF3DbAhGi/fiPgfnef\nYWZvA5PN7EKgDOgJ4O4LzGwysADYCFyihBeRMMrqgu3atXDbbXDDDXTo0gWmP8aOhx3G7Kyjvpxu\nzhIRyZNkFf2W5atXw/jxwdfxx8PAgbD//hXWy8vNWSIiUnta8WUQ8O3bw9KlMHs2TJpUIfBzRXPv\niIgUymefwfXXs4gJsOZMeOstKCqqsEp891AuqHtHRCRPtnTjlJXBqFHw8MNw/vnsNrYPn/tuaW2v\n7h0RkXqiI4vgvPPgkEOgZUtYtAjGjOELqg/8XFH3johIbZs/H4YPZzYzYO/L4aOPiMxrQeSW4O1u\n3SB221JtT8Og7h0REWppLpw33oDSUpg7F/r0YbtrLmadb5d1G3PRvaPQFxGpJOnNUqT54fDSS0HY\nL1wIffvChRfCNtuk3G+67VLoi4hkIJ3QTjecK6znDtOmBWH/+ecwYACccw40a5Z4/Swo9EVEaqDa\nm6XS2X7TZpg6FYYNg3XrgvH2Z54JTapeMs0m9Ct/SM2apdAXEclKjUJ/0ybOaPIoD+8/LAj4wYPh\n1FOhUcVBkbm+VmCm0BcRyUpWob9hA0ycCMOH8+qHv+A3zwyGk06K9b3UOoW+iEiWMgr9n36Ce+4J\nbqraay+49lrs2GJqkL9ZqWnoa5y+iEgq338Pd9wBY8bAgQfCgw/Cb35T6FZlTaEvIpLIt9/CrbfC\nuHHw29/ClClw8MGFblWNKfRFROJ9/TUl3Ah73xr01c+cCZ06FbpVOaPQF5F6qaajYmLbDh0avD7l\n0BUc+dpYurx3F204DebMCaY6bmB0IVdE6r1sb3oyA//kU8bvMZrLWj4If/0r9O2L7dEu85uz8kSj\nd0Qk9LIK348+4j/tR3LRTo8z+psL6frYVbw4f1cg9V8OtTJHTwYU+iISepmE/tx7F7DVDSPosPg5\nHtn5ElaccQX9r/sFM2fmJ7RrKi+hb2btgPuAXQAH7nD3m8xsKPA34MvoqgPd/bnoNgOAC4BNwOXu\n/mKlfSr0RSRr8RV3SQkMGRK8Tlpxv/12MFXCyy/DlVfCJZfAjjsChemmyVa+Qn9XYFd3n2dm2wFv\nAacCPYG17j620vqdgInAYUBbYDrQ0d03x62j0BcJkdrsFkkZ2q+9FoT922/D1VfD3/8O226b/vZ1\nTF5uznL3FcCK6Ot1ZvYBQZgDJDp4D2CSu28AysxsCdAFmJNtQ0Wkfqs8i2W2z31N9uER/xr3YEFp\nKXz0EfTrB48+Cltvnd1BG5CM+/TNrAiYBfwa6AOcD3wLvAn0cfc1ZjYemOPuD0a3uRN4zt0fi9uP\nKn2RBi5ZQJeU5KayjlXoWyp1d3juuaCy//LLYMbLv/4VmjZNaz/1QV6nYYh27TwKXBGt+G8D/hl9\n+1/AGODCJJtXOaVDY88HA4qLiymuD1dRRCRtyar7kpLcHsfYDI8/GYT9zz/DoEFw+unQuHFuD1QA\nkUiESLZ/FiWQdqVvZk2Bpwkq9nEJ3i8Cprr7/mbWH8DdR0bfex4Y4u6vx62vSl8kROKr6Uwr62R/\nMZSWbGTjAw/z/tnD+fWhzYPpjf/whyrTG2fStrouXxdyDZgAfO3uveOWt3H3L6KvewOHuftf4i7k\ndqH8Qm77+JRX6Is0HJk+jaomIWsG/tPPcP/9LPnbCNp3bcPvZl/Li5v/K+vpjRX6VQ/SFXgJmE95\nN81A4Cygc3TZUuAf7r4yus1AgiGbGwm6g16otE+FvkgDlM50xVmH7I8/cmnzu7ml3Si+ab0Pf3xz\nEMcM6VZlyCZU/yFU6JussqWbs0SkTsk29FOG8KHr4PbbYcwYpqw4jO6vD4YuXapeyE3Rnvoa8pUp\n9EWkTslFpb/l/TVrYPz44OvYY2HgQOzAA6rsJ53QT6eN9YEeoiIidU5aY+lTaMWXMGhcUN2fcgq8\n9BLss09tNDV0VOmLSE5VrqITVeOV14l9SGy/9nN2nzyG45ffy9JDTmdjn34ccdZeCfcficAxxwR9\n+am6a1TpV6RKX0RyJlbdn3de8LqoqPzn2PeiIujWLZjHHqIhXVRG8arR8NBD3LC6Fwcwn8/ebEsq\nsWCPu91H0qBKX0RSyvQCaOXKPmW/+4cfwogRwaMIL7oIrroKa70LkLqPPjZSJ5P2VLesvtCFXBHJ\nm3TCMj7kZ85M3AVz8h7vcti04TB9Olx2WfDVsuWW7SGzC7OZtlmhXwAKfZH6Jz4sq5tXJ2GF/8Yb\nwVQJc+ZA797B9Mbbb1/lGKDQT0ahLyJ5k8lwzArfX54dzHj5/vvQty9ceCE0b570GKDQT0YXckUk\nr5JV+FUv0jrvjJnBTEqh16fQvz889RRstVWNjlnlInBxVr9GaKnSF5G0pTMcE3d4+mle717K4ft8\nx9kLB/HAhjOhSfIas/JTsCC4DpBtqFd38TnMlb5CXySEsp2SIGXob9wEjz0W9Nk3asSf5g3msU1/\nxBo3SitgY20qKQmq+VxPk6BpGKLbK/RFwi2Tqjdh6P+8gXObTeK+Xw2HFi2C6Y1PPhlrZNVOkVDT\n9oSR+vRFpFqpqtys/fQTf+de+NUozmdPuOWWYH6cLKc3lvxQ6Is0IKnCPXbx0yx4Hb9uRhdGf/gB\n/vMfuO46enAA3H8/x3Y9Cj+uYjvi91vd/uPbknF7JCPq3hFpoNIZXhn/c7V93t99R78db2NU6xvg\nyCNh8GDs0EOSduGkMwNmJu2WgPr0RSShTEM/6TrffAM33QS33MLEr/6Lv7w7EPbbr8J6lT84ysrK\n594pKws+NIqK0qvcFfqpKfRFJKGahv7sx1fhY8bS+c3/8EyTU5nSqT+vrOpQIcCPOSa7aj6bdktA\noS/SQNV0iGHWob98OTe2u44rWt4PZ50FfftiRXsmfNJVOlMbZ0qhn5pCXyQEsgnCZI8krBzUJSXB\nxGjFe3wMo0bBI49w/eoLuPrzPtCmTcrj10ZAK/RTy9eD0dsB9wG7EDwE/Q53v8nMdgIeBvYEyoCe\n7r4mus0AggejbwIud/cXK+1ToS+SplyFfvzy2Pd97QM+OGcEPPssXHwxXHkltnOrpH8N1PZNTgr9\n1PIV+rsCu7r7PDPbDngLOBU4H/jK3UebWT+gpbv3N7NOwETgMKAtMB3o6O6b4/ap0BdJU22Efmeb\nx7zTh7PykVm0Lr0cLr00uLkqwbb5DGKFfmp5uTnL3VcAK6Kv15nZBwRh3h3oFl1tAhAB+gM9gEnu\nvgEoM7MlQBdgTrYNFZGai0SgC68zr10pz/IWt7x1Nf24mz8v3o6icYUbF68J1fIn4z59MysCZgH7\nAZ+4e8vocgO+cfeWZjYemOPuD0bfuxN4zt0fi9uPKn2RNCV7pmzsdaIuliqToL30EpSWsmz6h+x5\nSz+2vvQC1vvWaU09rOq77sjrNAzRrp3HgCvcfa3F3W7t7m5mqf5ZVHlvaNzDLYuLiynWR7pIWiqH\ne/zdrBW4wwsvBJOgrVgBAwfSYfpf+fmSZvx0aX7aKjUTiUSIJP0PnLm0K30zawo8TVCxj4suWwgU\nu/sKM2sDzHT3fcysP4C7j4yu9zwwxN1fj9ufKn2RNKWq9EtKyp8Z26IFrFkD5pvZ8NhTXPLNMJps\nXM9bJw5iblFPvFHjLTdNTZgAvXqV3zwFyaceVqVfd+TrQq4R9Nl/7e6945aPji4bFQ36FpUu5Hah\n/EJu+/iUV+iLpC9V6FZ4b9MmmDw5qOy32gquvRa6d4dGjbLfZzXrSn7lK/S7Ai8B8ynvphkAzAUm\nA3tQdcjmQIIhmxsJuoNeqLRPhb5ImqoN6J83wAMPwIgRsMsuwfTGJ5yQcsbLTK4TxO68lcLTzVki\nIZA09Nev55Jt7ubWPUZBhw5BZX/00dVOb5zO3bQN5aEjDY1CX6SBShm6h66Df/8bxoxh6heH8IfX\nBsERR2S0/2xnwZTC0kNUROqJTCvnhMu//RZuvhl63hi8+eyzdD+oM55m3iebtz6+PdKwqdIXKYCM\nq+uvvoJx4+D22+Hkk6F/f9h33+z2VakNqvTrF1X6InVErfSBf/EFjBkDd98Np58Oc+fCL39Zw5ZK\nmCn0RXIk/gLorFnly2IfBBkF/7JlMHo0TJoEZ58N8+fD7rsnXT3dD5xk0x2oeyc81L0joVbdA8Oz\nrdxjg2fSHQe/xeLFMHIkPPkkXHQR9O4NrVunbHf8zVnZDK9U9079otE7IjmSyc1K6ewLMtjfe+/B\n8OEwbVow2+Xll8NOO2XVbj2TtmFTn75IffbWW8Hds6++GlT1t98OO+xQ7WapZqUUSUWVvkhUPiv9\no+wVXjlpWNBXf801QVdO8+aZNTjFsVXpN1yq9EXqC3eYMQOGDeN+lkGPfvDEE8EcOSJ5okpfJKrW\nKn13eOYZKC0NpsAcOJCmvc5igzetcZuzbaumWKi/dCFXJEdyHfrGZjY/8ngQ9u4waBD86U/QuHGt\ndqmou6ZhU/eOSJYqV7sQXBCtcbW7cSN/5SEGMhyu2z4I/ZNPJjLLiPwrmL9+zz2DY8Tmsi8qUpUt\n+aFKX4QajKuP99NPcN99MHIkkY/bUcpgpm8+LqPpjXNBlX7DVtNKv+qTFUQkMz/+COPHs75de5aM\neoy7j76Xod0izOB4hpZY8kcZihSAundEsrV2Ldx2G9xwAxxxBFs/+wTL1x3KJxEo3rN8taymYchQ\nqnH76jKSeOrekVCLhWVJSRCW6UzDcPwhq+n69ngYPx6OPx4GDoT998/q+OqKkUxp9I5IDpjBzJnV\nDGNctYplvW9g5yfuYNGvujOm6QDa/75jxXWyOK7+N5BMKPQl9HIx5jzl/DWffQbXXw8TJsCZZ0Lf\nvlBUlJPAVuhLpvL1YPS7gZOBVe6+f3TZUOBvwJfR1Qa6+3PR9wYQPBR9E3C5u7+YYJ8Kfcm5mj5Q\npMLPS8tg1Ch4+GE47zzo0wfatq3xsXLRXgmvfI3TvwcYD9wXt8yBse4+tlKDOgFnAJ2AtsB0M+vo\n7puzbaRIMrUy1n7RIu5hBBwyFS6+GBYtgp13rmFLReqGtELf3V82s6IEbyX6tOkBTHL3DUCZmS0B\nugBzsm2kSDLx4R4bDh8buZKx+fODGS9nzmQJl8NHH0GLFjVvpEgdUtMhm5eZ2bnAm0Afd18D7EbF\ngF9OUPGL5FXaff1z5/Ikw+CEuXDVVXDXXQzbfjtKlffSANUk9G8D/hl9/S9gDHBhknUT9loOjSvJ\niouLKdaAYqkls2ZV7e55+8aX2OGmUlp9tZCl7ftR2vMhNn6/DcVv1m5bNKZeMhGJRIjk8A6/tEfv\nRLt3psYu5CZ7z8z6A7j7yOh7zwND3P31StvoQq7kVKqpFLZcMHUPnk5VWgqffw4DBsA550CzZonX\nr+Z4+ics+VawCdfMrI27fxH98Y/Au9HXU4CJZjaWoFunAzA32+OI5IKxGZ6aGoT9998HM16ecQY0\n0U3pEi5p/Ys3s0lAN6CVmX0KDAGKzawzQdfNUuAfAO6+wMwmAwuAjcAlKumlYDZtgkcfZR7D4J9N\ng7A/9VRopGmnJJx0c5bUC+lclK3QvbNhAzz4IIwYAa1a8ftXB/Hs5pNSzngZT907UlfpjlwJnWRh\nawZbsZ71t94T3FTVvj0MHgzdumGNLCdPk9ITp6TQFPoSOpVDPxKBV178nq9G3MHAZtezes+DeOm3\ng2h/zpEV/grQPzdpCPTkLAm3b7+l+LVbKb5rHI/yW3ae8zQ7H3QQHQvdLpE6SpW+5F1Nu0jMwL/6\nGm68kQ033sqCopOY3XUAj7zfKem+VOlLQ6HuHanTqgv4jMN4xQpGtxlL7x3uYv7ep/FK1348Pr99\ntR8cCn1pKBT6khe1MX1xbFm189gDfPIJXHcdPPggN60+m8s/uQbatUu638ptLimBIUMyb7NIXaPQ\nl7zL1fTFiZZVWWfJEhg5Ep54Av72N7jqKmzX1qm3idJIG2mIFPqSd3kJ/QULYPhweP55uOQSuOIK\n+MUvUm8jEgI1DX3dlih1Smfehj//GY45BvbbDz7+GP75zy2BLyI1o9CXumHOHDjlFJ7mFDjqqCDs\n+/eHHXYodMtEGhSN05eCiEQgMtMpJsJbLUvZ2z9idtf+/JlHWd9768TrR4LXlacjFpH0qU9fMlbj\nPn13eO45KC1l0WtfM5yBTPj5L9C0aVb7Vp++hIku5EreZRuyjWwzmx97MpjeeONGGDiQxmedzmYa\nb9mfQl8kNYW+5F18yMa6XcrKgu9FRcHr4uLgdXExFHfdCA8/zPtnD+fXh20bTIJ2yinQqFGVB58o\n9EVS09w7UicUFcGyZXDeecGNUOedB8W/+ZmFg+7nm9NGsHb73biKsRx50u8oe9zg8WCbbt2CRxkO\nHar+eZF8UKUvGUs1tbE7bGM/8uP4u1j/r9EsabwPd+wymClrjmbZsiDki4qiHwrFVfeXbtWuG68k\nrNS9IzmVLExbtIA1a6oujw/Z7Wwd6667nc+vGctu3Q8LnlLVpQtQ/fNrY1MkKMBFUlPoS61Ju698\nzRoYP55V/zueXc44lgMfHsA7fmCVfUE1Dy0XkWop9KXWVDenTfPvv6TN5HH86cvb+bDjHzjjnQEs\n8l8lnW4BFPoiNZWXC7lmdjdwMrDK3fePLtsJeBjYEygDerr7muh7A4ALgE3A5e7+YrYNlMKIdfEM\nHZqgy6Xj5xRPHQP33MPtq3vS/OM36bzXXnyY9T9DEcmXtCp9M/stsA64Ly70RwNfuftoM+sHtHT3\n/mbWCZgIHAa0BaYDHd19c6V9qtKv42IV+JZKvKwMRo+Ghx6Cc8+Fa67Bdm9b5SKsKn2R2pOXSt/d\nXzazokqLuwPdoq8nABGgP9ADmOTuG4AyM1sCdAHmZNtIKawOfAjnj4ApU+Dvf4eFC2GXXWq0z1TT\nKujirUjtqck4/dbuvjL6eiXQOvp6NyoG/HKCil9yKB9DFvfjXThrOK8wHYr+J5jbvmXLnOxb4S5S\nGDm5Ocvd3cxS/YGe8L2hsfIOKC4uplgpkLb40DQr/wDIiTffhGHDeJE5cFBvfvnQHawdsn0ODyAi\n6YpEIkRy+D942qN3ot07U+P69BcCxe6+wszaADPdfR8z6w/g7iOj6z0PDHH31yvtT336OZKzPvHZ\ns4N5cd5/H/r1Y5vLLuRH3ybl/mNj7MvKYMKEoKumyjQMxan79EUkfYWchmEK0AsYFf3+ZNzyiWY2\nlqBbpwMwtwbHkSyl1QXkDtOnB2G/fHkwh/2UKdCsGesvS/9YsSkVKh8jEinvr1ffvUjhpTt6ZxLB\nRdtWBP33/ws8BUwG9qDqkM2BBEM2NwJXuPsLCfapSj9H0qn0q6zjDk8/HYT92rUwcCCceSY0aVJl\nm+oq/WTvaaoEkdzTzVmS1jQGW8J50yZ47DEYNgwaNQqmSjjttOA1iYO6pARmzkwc1BpuKZJfCn1J\na8KypraBDfdODB423rIlXHstkea/JzIr+LeTqhLPttIXkdxT6Evq0P/pJ7j3Xj6+eBS/PKYoqOyP\nPbb8ymqCfaTaP6jbRqSQFPqSOPR/+AHuuAOuvx4OOICjnhvEK35UWvsABbtIXaXQlwqBvYN9x3cj\nboVx4+Coo4LK/uCDq+2GUTeNSP2gJ2c1EDWurL/5Bm66iY+4Bd79HW+MmM4zy/aDKRC5Klgl9nQq\nVeoi4aVKvw4yC0bL3HtvcKNTWVmwvKio/KlTEHw4bLtuJW0nj+XUr+5k4a9Opee8ASzx9lX2B6r0\nRRoCde80QJUDOFFovzb5UxqNuY4D5j/A8zv9hbKeffl2xz0oKclshstkxxSRukmh3wClDP2PP4aR\nI+HRR+GCC6BPH2y3NimHbCr0RRoO9emHxD58AOeOgGefhf/+b/jwQ2jVKuU2mr5YRCpTpV8HVai6\n33mHyZ2H0Y1ZtB52BVx6Key4Y9L1k1X6GU/TICJ1krp36qCajsQxA39tTjBVwltv0eeLPtzOxXzv\n2yZdP9U0DAp9kYZDoV/HZRSm7jBrFtOPKeX4PRZDv35wwQXYNltveTubY1T3gPPYa92AJVL3KfTr\nuLRC3x1eeCGY8XLlSi5YMoC7fzobmjXbso/YatkcQ1W8SMOh0K/jYoGbqKo238yfmjzFfk8Ng/Xr\ng7tne/bEmjSudshmomNU1wYRqf8U+nVcosBtbJvYNHFy0Ge/9dYweDB0775leuN0xulXd4xM3heR\n+kNDNgsgVrWXlQXfi4pg4ULYdVdo0QLWrIFTT02w4c8/wwMP8AEj4ZZdgsnQTjgBzFIOr6xJGxPt\nT/31IuGlSj8q24uaiZ4uVeX1j+vh7rth9Gjo0IFu0wcza/PRVaY3rtyeVNMwpDvfvYg0LOreqQWZ\nhGjK0F+3jj7b/5sxbcbAoYcGffaHH57zkFboi4SHuncKJP4vg1jVHZvFckfWwLBb4MYbOYJuwV20\nnTvX2vHVfSMi6apxpW9mZcB3wCZgg7t3MbOdgIeBPan00PS47RpMpR/jX34F48bx1bDbaXXuydC/\nP9Zp36zukBURSaSmlX6jHLTBgWJ3P8jdu0SX9QemuXtHYEb05wZrV77gOq6Gjh3hyy/pwlyYMAH2\n3bfQTRMRqSBX3TuVP3W6A92irycAEepY8Ke6cJvuNj06L+NmRnMWk7iPc2H+fNh9d5beUd7doq4X\nEalLctG98zHwLUH3zr/d/T9mttrdW0bfN+Cb2M9x29WZ7p3Kc9fMmhX8nDSkFy/mro4juXCnJxn5\nzUXcQG9W0TrlpGeVj1dHfnURqWfqwoXco9z9CzPbGZhmZgvj33R3N7OEETc0VgIDxcXFFBewDI41\nJdZHH9e0cu+9B8OHw7RpfMqlsHgxA36xU55aKCJhFIlEiMS6GHIgp0M2zWwIsA64iKCff4WZtQFm\nuvs+ldYtWKVfuZsmvrI/5phgeXzT3vz3W2w7bhjtPn2VB3fpzTdn/DcDR+7AzJnl68dvk6iS1+Rm\nIpILBR2nb2bNgcbuvtbMtgVeBEqA44Gv3X2UmfUHWrh7/0rbphX6mYRlonXjb2xKNu0wVAxsCD4E\nvnrqFa5YV0rrle+y6tyraT/679i2zSuMy68weifN7h0RkWwVOvT3Ap6I/tgEeNDdR0SHbE4G9iCH\nQzazuWkqnWUQH9jOccxgevEwPo4s45f/7g+9esFWW21ZP1bh9+oVfJAsWwatWyeehkGVvIjkUkH7\n9N19KVDlriN3/4ag2s+r+Eofko+aqbwewHm9nAOXP8NrlNKm+bc83mIgPTmL6R2bULxVxXVj+7r3\n3uC7GaxYkZNfQUSkVtWraRhqo9JvZJs5jcd59MBScOf0+YN5ZONp0Lhxyn1U94hCEZHaEKq5d3Ia\n+hs3wqRJLDh3BN+xA0dMHQwnn4w1srSeN1vdIwpFRGpD6EI/3aBNGtjrfwrulh01CvbYg+Mjg5jB\nccTOYboPGVd1LyKFELrQz/oGqB9+4PJt7+Sm3a+DX/86eHBJ164JR+8o9EWkrqoLN2fVbWvXwm23\nwdixHMOR8MQTcOihwcXcocEq8VMlpBK7+Dt0qKZXEJH6qcFW+i1tNauH3gQ33wzHHw+DBmH771ft\nxd10K32T3jVoAAAHfElEQVQRkUJQpV/ZqlVwww0s4Q5Y1gNeeSWY/TKFZHPTi4g0NHWi0k/3rtuU\nVfjy5cEzZ++7D846iz1v7csy37PSMYNt4o9XUlJ+cTjZ8TSFgojUFQ3uQm7lME8WuCUl0fWWLg1G\n4kyezKfHncfkPa5m7fa7bQnzdKZhSOd4CngRqQsafOgn8ytbxKJeI+Dpp+Ef/4Arr4Sdd854P+qj\nF5H6pMH36VeuvM/Y5x1+O3s4LzMT9r4MliwJJrwREZFq1Z9Kf+5cphxeygk7vcGsQ/pw2rSLObTb\ndlWeeJVO14y6cESkvmow3TuxIC4pCUbQxMK3R8uXOOiZUli4kP/5tC83/3AhbLNNlZuqRETCoMGE\nfvlyAMeff5E1fYex+dPPebnrAMavOYcZLzejW7fgomxZWRqPNRQRaWAaVuhv3kyPxlMZTCmHdfoB\nBg2Cnj2hSZOElb0uwopI2DSM0N+0CR55BIYN4//ea0opg3l806nQqFHc+sF3hb6IhFn9Dv2ff4YH\nHoARI6BVK7j2Wuz3JwKWcPoDUOiLSLjV7yGbHTpA+/Zwxx3B1VvL+vcQEZE01Fqlb2YnAuOAxsCd\n7j6q0vvur74KRx5ZZQhlogu0qvRFROpo946ZNQYWETwn9zPgDeAsd/8gbp2ks2ym8wDzZOuJiDRk\nNQ39RtWvkpUuwBJ3L3P3DcBDQI9aOpaIiKSptvr02wKfxv28HDg8mx0lm/ZY4/JFRDJXW907fwJO\ndPeLoj+fDRzu7pfFrZNx944eaCIiYVdXR+98BrSL+7kdQbVfwdC4p5W0aFHMmjXFQOpHEaaq/FX9\ni0hDE4lEiMRCLwdqq9JvQnAh9zjgc2AuGVzITbxPVfUiInVy9A6AmZ1E+ZDNu9x9RKX3qw19zYYp\nIlJRnQ39ag+cxYPRRUTCrq4O2RQRkTpIoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGF\nvohIiCj0RURCRKEvIhIiCn0RkRBR6IuIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIgo9EVEQkShLyIS\nIgp9EZEQyTr0zWyomS03s7ejXyfFvTfAzBab2UIz+11umioiIjVVk0rfgbHuflD06zkAM+sEnAF0\nAk4EbjUz/UWRQiQSKXQT6gydi3I6F+V0LnKnpmGc6InsPYBJ7r7B3cuAJUCXGh6nQdM/6HI6F+V0\nLsrpXOROTUP/MjN7x8zuMrMW0WW7Acvj1lkOtK3hcUREJAdShr6ZTTOzdxN8dQduA/YCOgNfAGNS\n7Mpz12QREcmWudc8j82sCJjq7vubWX8Adx8Zfe95YIi7v15pG30QiIhkwd0Tda2npUm2G5pZG3f/\nIvrjH4F3o6+nABPNbCxBt04HYG7l7WvSaBERyU7WoQ+MMrPOBF03S4F/ALj7AjObDCwANgKXeC7+\nnBARkRrLSfeOiIjUDwUZP29mJ0Zv3FpsZv0K0YZCMbN2ZjbTzN43s/fM7PLo8p2iF84/NLMX40ZD\nNXhm1jh6g9/U6M+hPBdm1sLMHjWzD8xsgZkdHuJzMSD6/8i7ZjbRzLYKy7kws7vNbKWZvRu3LOnv\nnunNsHkPfTNrDNxMcONWJ+AsM9s33+0ooA1Ab3f/NXAEcGn09+8PTHP3jsCM6M9hcQVBd2Dsz86w\nnosbgWfdfV/gAGAhITwX0YEhFwEHu/v+QGPgTMJzLu4hyMd4CX/3bG6GLUSl3wVY4u5l7r4BeIjg\nhq5QcPcV7j4v+nod8AHBBe/uwIToahOAUwvTwvwys92B3wN3Un6zX+jOhZntCPzW3e8GcPeN7v4t\nITwXwHcExVFzM2sCNAc+JyTnwt1fBlZXWpzsd8/4ZthChH5b4NO4n0N781a0ojkIeB1o7e4ro2+t\nBFoXqFn5dgNwDbA5blkYz8VewJdmdo+Z/Z+Z/cfMtiWE58LdvyG47+cTgrBf4+7TCOG5iJPsd8/4\nZthChL6uHANmth3wGHCFu6+Nfy862qnBnyczOwVY5e5vk3hKj9CcC4KRdAcDt7r7wcD3VOq+CMu5\nMLO9gSuBIoJQ287Mzo5fJyznIpE0fveU56UQof8Z0C7u53ZU/KRq8MysKUHg3+/uT0YXrzSzXaPv\ntwFWFap9efQboLuZLQUmAcea2f2E81wsB5a7+xvRnx8l+BBYEcJzcSjwqrt/7e4bgceBIwnnuYhJ\n9v9E5TzdPbosqUKE/ptABzMrMrNmBBchphSgHQVhZgbcBSxw93Fxb00BekVf9wKerLxtQ+PuA929\nnbvvRXCh7v+5+zmE81ysAD41s47RRccD7wNTCdm5ILiAfYSZbRP9/+V4ggv9YTwXMcn+n5gCnGlm\nzcxsL5LcDFuBu+f9CzgJWERw0WFAIdpQqC+gK0H/9Tzg7ejXicBOwHTgQ+BFoEWh25rn89INmBJ9\nHcpzARwIvAG8Q1Dd7hjic9GX4EPvXYILl03Dci4I/ur9HPiZ4Prn+al+d2BgNEsXAidUt3/dnCUi\nEiJ6uImISIgo9EVEQkShLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJkf8PDkxgYpLgjk4A\nAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "You should now be able to fit simple models to datasets! Note that for more complex models, more sophisticated techniques may be required for fitting, but ``curve_fit`` will be good enough for most simple cases." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Note that there is a way to simplify the call to the function with the best-fit parameters, which is:\n", "\n", " line(x, *popt)\n", "\n", "The * notation will expand a list of values into the arguments of the function. This is useful if your function has more than one or two parameters. Hence, you can do:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.errorbar(x, y, yerr=e, fmt=\"none\")\n", "plt.plot(xfine, line(xfine, *popt), 'r-')" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "pyout", "prompt_number": 12, "text": [ "[]" ] }, { "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOXZ//HPxabiBhZFRDRWoEpdcEOtVOLyVK0WrK2o\nrYpLbX30cUGU3QfShlVBFLdaN1xAcQd34MegqIj6iKgIghIUFXABBRVluX5/nBkySWYmM5PJTJLz\nfb9eeWVy5ix3jvidK/e5z33M3RERkXBoVOgGiIhI/ij0RURCRKEvIhIiCn0RkRBR6IuIhIhCX0Qk\nRNIKfTPb2sxeN7N5ZrbAzEZElw81s+Vm9nb066S4bQaY2WIzW2hmv6utX0BERNJn6Y7TN7Pm7v6D\nmTUBZgNXA8cBa919bKV1OwETgcOAtsB0oKO7b85l40VEJDNpd++4+w/Rl82AxsDq6M+WYPUewCR3\n3+DuZcASoEsN2ikiIjmQduibWSMzmwesBGa6+/vRty4zs3fM7C4zaxFdthuwPG7z5QQVv4iIFFAm\nlf5md+8M7A4cbWbFwG3AXkBn4AtgTKpd1KCdIiKSA00y3cDdvzWzZ4BD3T0SW25mdwJToz9+BrSL\n22z36DLi1teHgIhIFtw9Ubd6WtIdvdMq1nVjZtsA/wW8bWa7xq32R+Dd6OspwJlm1szM9gI6AHMT\nNFxf7gwZMqTgbagrXzoXOhc6F6m/airdSr8NMMHMGhF8UNzv7jPM7D4z60zQdbMU+Ec0zBeY2WRg\nAbARuMRz0VoREamRtELf3d8FDk6w/NwU2wwHhmffNBERyTXdkVsHFBcXF7oJdYbORTmdi3I6F7mT\n9s1ZOT+wmXp8REQyZGZ4bV/IFRGRhkGhLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJEYW+\niEiIKPRFREJEoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCZGMH4wuIiL5E4kEX7HXNaWH\nqIiI1BNmAHl4iIqZbW1mr5vZPDNbYGYjost3MrNpZvahmb1oZi3ithlgZovNbKGZ/S7bBoqISO6k\nXembWXN3/8HMmgCzgauB7sBX7j7azPoBLd29v5l1AiYChwFtgelAR3ffHLc/VfoiIhnIRaWfdp++\nu/8QfdkMaAysJgj9btHlE4AI0B/oAUxy9w1AmZktAboAc7JtqIhIQ1K5rz727Pfi4vLXtSHt0Dez\nRsD/AXsDt7n7+2bW2t1XRldZCbSOvt6NigG/nKDiFxERKoa7WW4u0qYj7SGb7r7Z3TsDuwNHm9kx\nld53IFV/jfpyRERqoDEba7yPjIdsuvu3ZvYMcAiw0sx2dfcVZtYGWBVd7TOgXdxmu0eXVTB06NAt\nr4uLiymuzb9pRETqoUgkQmTGDHjnHS7hZcbXcH9pXcg1s1bARndfY2bbAC8AJcAJwNfuPsrM+gMt\nKl3I7UL5hdz28VdudSFXRCRgBgnj8McfWdzvTnaZcB1f7tyJ8z8axGyOzsuF3DbAhGi/fiPgfnef\nYWZvA5PN7EKgDOgJ4O4LzGwysADYCFyihBeRMMrqgu3atXDbbXDDDXTo0gWmP8aOhx3G7Kyjvpxu\nzhIRyZNkFf2W5atXw/jxwdfxx8PAgbD//hXWy8vNWSIiUnta8WUQ8O3bw9KlMHs2TJpUIfBzRXPv\niIgUymefwfXXs4gJsOZMeOstKCqqsEp891AuqHtHRCRPtnTjlJXBqFHw8MNw/vnsNrYPn/tuaW2v\n7h0RkXqiI4vgvPPgkEOgZUtYtAjGjOELqg/8XFH3johIbZs/H4YPZzYzYO/L4aOPiMxrQeSW4O1u\n3SB221JtT8Og7h0REWppLpw33oDSUpg7F/r0YbtrLmadb5d1G3PRvaPQFxGpJOnNUqT54fDSS0HY\nL1wIffvChRfCNtuk3G+67VLoi4hkIJ3QTjecK6znDtOmBWH/+ecwYACccw40a5Z4/Swo9EVEaqDa\nm6XS2X7TZpg6FYYNg3XrgvH2Z54JTapeMs0m9Ct/SM2apdAXEclKjUJ/0ybOaPIoD+8/LAj4wYPh\n1FOhUcVBkbm+VmCm0BcRyUpWob9hA0ycCMOH8+qHv+A3zwyGk06K9b3UOoW+iEiWMgr9n36Ce+4J\nbqraay+49lrs2GJqkL9ZqWnoa5y+iEgq338Pd9wBY8bAgQfCgw/Cb35T6FZlTaEvIpLIt9/CrbfC\nuHHw29/ClClw8MGFblWNKfRFROJ9/TUl3Ah73xr01c+cCZ06FbpVOaPQF5F6qaajYmLbDh0avD7l\n0BUc+dpYurx3F204DebMCaY6bmB0IVdE6r1sb3oyA//kU8bvMZrLWj4If/0r9O2L7dEu85uz8kSj\nd0Qk9LIK348+4j/tR3LRTo8z+psL6frYVbw4f1cg9V8OtTJHTwYU+iISepmE/tx7F7DVDSPosPg5\nHtn5ElaccQX9r/sFM2fmJ7RrKi+hb2btgPuAXQAH7nD3m8xsKPA34MvoqgPd/bnoNgOAC4BNwOXu\n/mKlfSr0RSRr8RV3SQkMGRK8Tlpxv/12MFXCyy/DlVfCJZfAjjsChemmyVa+Qn9XYFd3n2dm2wFv\nAacCPYG17j620vqdgInAYUBbYDrQ0d03x62j0BcJkdrsFkkZ2q+9FoT922/D1VfD3/8O226b/vZ1\nTF5uznL3FcCK6Ot1ZvYBQZgDJDp4D2CSu28AysxsCdAFmJNtQ0Wkfqs8i2W2z31N9uER/xr3YEFp\nKXz0EfTrB48+Cltvnd1BG5CM+/TNrAiYBfwa6AOcD3wLvAn0cfc1ZjYemOPuD0a3uRN4zt0fi9uP\nKn2RBi5ZQJeU5KayjlXoWyp1d3juuaCy//LLYMbLv/4VmjZNaz/1QV6nYYh27TwKXBGt+G8D/hl9\n+1/AGODCJJtXOaVDY88HA4qLiymuD1dRRCRtyar7kpLcHsfYDI8/GYT9zz/DoEFw+unQuHFuD1QA\nkUiESLZ/FiWQdqVvZk2Bpwkq9nEJ3i8Cprr7/mbWH8DdR0bfex4Y4u6vx62vSl8kROKr6Uwr62R/\nMZSWbGTjAw/z/tnD+fWhzYPpjf/whyrTG2fStrouXxdyDZgAfO3uveOWt3H3L6KvewOHuftf4i7k\ndqH8Qm77+JRX6Is0HJk+jaomIWsG/tPPcP/9LPnbCNp3bcPvZl/Li5v/K+vpjRX6VQ/SFXgJmE95\nN81A4Cygc3TZUuAf7r4yus1AgiGbGwm6g16otE+FvkgDlM50xVmH7I8/cmnzu7ml3Si+ab0Pf3xz\nEMcM6VZlyCZU/yFU6JussqWbs0SkTsk29FOG8KHr4PbbYcwYpqw4jO6vD4YuXapeyE3Rnvoa8pUp\n9EWkTslFpb/l/TVrYPz44OvYY2HgQOzAA6rsJ53QT6eN9YEeoiIidU5aY+lTaMWXMGhcUN2fcgq8\n9BLss09tNDV0VOmLSE5VrqITVeOV14l9SGy/9nN2nzyG45ffy9JDTmdjn34ccdZeCfcficAxxwR9\n+am6a1TpV6RKX0RyJlbdn3de8LqoqPzn2PeiIujWLZjHHqIhXVRG8arR8NBD3LC6Fwcwn8/ebEsq\nsWCPu91H0qBKX0RSyvQCaOXKPmW/+4cfwogRwaMIL7oIrroKa70LkLqPPjZSJ5P2VLesvtCFXBHJ\nm3TCMj7kZ85M3AVz8h7vcti04TB9Olx2WfDVsuWW7SGzC7OZtlmhXwAKfZH6Jz4sq5tXJ2GF/8Yb\nwVQJc+ZA797B9Mbbb1/lGKDQT0ahLyJ5k8lwzArfX54dzHj5/vvQty9ceCE0b570GKDQT0YXckUk\nr5JV+FUv0jrvjJnBTEqh16fQvz889RRstVWNjlnlInBxVr9GaKnSF5G0pTMcE3d4+mle717K4ft8\nx9kLB/HAhjOhSfIas/JTsCC4DpBtqFd38TnMlb5CXySEsp2SIGXob9wEjz0W9Nk3asSf5g3msU1/\nxBo3SitgY20qKQmq+VxPk6BpGKLbK/RFwi2Tqjdh6P+8gXObTeK+Xw2HFi2C6Y1PPhlrZNVOkVDT\n9oSR+vRFpFqpqtys/fQTf+de+NUozmdPuOWWYH6cLKc3lvxQ6Is0IKnCPXbx0yx4Hb9uRhdGf/gB\n/vMfuO46enAA3H8/x3Y9Cj+uYjvi91vd/uPbknF7JCPq3hFpoNIZXhn/c7V93t99R78db2NU6xvg\nyCNh8GDs0EOSduGkMwNmJu2WgPr0RSShTEM/6TrffAM33QS33MLEr/6Lv7w7EPbbr8J6lT84ysrK\n594pKws+NIqK0qvcFfqpKfRFJKGahv7sx1fhY8bS+c3/8EyTU5nSqT+vrOpQIcCPOSa7aj6bdktA\noS/SQNV0iGHWob98OTe2u44rWt4PZ50FfftiRXsmfNJVOlMbZ0qhn5pCXyQEsgnCZI8krBzUJSXB\nxGjFe3wMo0bBI49w/eoLuPrzPtCmTcrj10ZAK/RTy9eD0dsB9wG7EDwE/Q53v8nMdgIeBvYEyoCe\n7r4mus0AggejbwIud/cXK+1ToS+SplyFfvzy2Pd97QM+OGcEPPssXHwxXHkltnOrpH8N1PZNTgr9\n1PIV+rsCu7r7PDPbDngLOBU4H/jK3UebWT+gpbv3N7NOwETgMKAtMB3o6O6b4/ap0BdJU22Efmeb\nx7zTh7PykVm0Lr0cLr00uLkqwbb5DGKFfmp5uTnL3VcAK6Kv15nZBwRh3h3oFl1tAhAB+gM9gEnu\nvgEoM7MlQBdgTrYNFZGai0SgC68zr10pz/IWt7x1Nf24mz8v3o6icYUbF68J1fIn4z59MysCZgH7\nAZ+4e8vocgO+cfeWZjYemOPuD0bfuxN4zt0fi9uPKn2RNCV7pmzsdaIuliqToL30EpSWsmz6h+x5\nSz+2vvQC1vvWaU09rOq77sjrNAzRrp3HgCvcfa3F3W7t7m5mqf5ZVHlvaNzDLYuLiynWR7pIWiqH\ne/zdrBW4wwsvBJOgrVgBAwfSYfpf+fmSZvx0aX7aKjUTiUSIJP0PnLm0K30zawo8TVCxj4suWwgU\nu/sKM2sDzHT3fcysP4C7j4yu9zwwxN1fj9ufKn2RNKWq9EtKyp8Z26IFrFkD5pvZ8NhTXPLNMJps\nXM9bJw5iblFPvFHjLTdNTZgAvXqV3zwFyaceVqVfd+TrQq4R9Nl/7e6945aPji4bFQ36FpUu5Hah\n/EJu+/iUV+iLpC9V6FZ4b9MmmDw5qOy32gquvRa6d4dGjbLfZzXrSn7lK/S7Ai8B8ynvphkAzAUm\nA3tQdcjmQIIhmxsJuoNeqLRPhb5ImqoN6J83wAMPwIgRsMsuwfTGJ5yQcsbLTK4TxO68lcLTzVki\nIZA09Nev55Jt7ubWPUZBhw5BZX/00dVOb5zO3bQN5aEjDY1CX6SBShm6h66Df/8bxoxh6heH8IfX\nBsERR2S0/2xnwZTC0kNUROqJTCvnhMu//RZuvhl63hi8+eyzdD+oM55m3iebtz6+PdKwqdIXKYCM\nq+uvvoJx4+D22+Hkk6F/f9h33+z2VakNqvTrF1X6InVErfSBf/EFjBkDd98Np58Oc+fCL39Zw5ZK\nmCn0RXIk/gLorFnly2IfBBkF/7JlMHo0TJoEZ58N8+fD7rsnXT3dD5xk0x2oeyc81L0joVbdA8Oz\nrdxjg2fSHQe/xeLFMHIkPPkkXHQR9O4NrVunbHf8zVnZDK9U9079otE7IjmSyc1K6ewLMtjfe+/B\n8OEwbVow2+Xll8NOO2XVbj2TtmFTn75IffbWW8Hds6++GlT1t98OO+xQ7WapZqUUSUWVvkhUPiv9\no+wVXjlpWNBXf801QVdO8+aZNTjFsVXpN1yq9EXqC3eYMQOGDeN+lkGPfvDEE8EcOSJ5okpfJKrW\nKn13eOYZKC0NpsAcOJCmvc5igzetcZuzbaumWKi/dCFXJEdyHfrGZjY/8ngQ9u4waBD86U/QuHGt\ndqmou6ZhU/eOSJYqV7sQXBCtcbW7cSN/5SEGMhyu2z4I/ZNPJjLLiPwrmL9+zz2DY8Tmsi8qUpUt\n+aFKX4QajKuP99NPcN99MHIkkY/bUcpgpm8+LqPpjXNBlX7DVtNKv+qTFUQkMz/+COPHs75de5aM\neoy7j76Xod0izOB4hpZY8kcZihSAundEsrV2Ldx2G9xwAxxxBFs/+wTL1x3KJxEo3rN8taymYchQ\nqnH76jKSeOrekVCLhWVJSRCW6UzDcPwhq+n69ngYPx6OPx4GDoT998/q+OqKkUxp9I5IDpjBzJnV\nDGNctYplvW9g5yfuYNGvujOm6QDa/75jxXWyOK7+N5BMKPQl9HIx5jzl/DWffQbXXw8TJsCZZ0Lf\nvlBUlJPAVuhLpvL1YPS7gZOBVe6+f3TZUOBvwJfR1Qa6+3PR9wYQPBR9E3C5u7+YYJ8Kfcm5mj5Q\npMLPS8tg1Ch4+GE47zzo0wfatq3xsXLRXgmvfI3TvwcYD9wXt8yBse4+tlKDOgFnAJ2AtsB0M+vo\n7puzbaRIMrUy1n7RIu5hBBwyFS6+GBYtgp13rmFLReqGtELf3V82s6IEbyX6tOkBTHL3DUCZmS0B\nugBzsm2kSDLx4R4bDh8buZKx+fODGS9nzmQJl8NHH0GLFjVvpEgdUtMhm5eZ2bnAm0Afd18D7EbF\ngF9OUPGL5FXaff1z5/Ikw+CEuXDVVXDXXQzbfjtKlffSANUk9G8D/hl9/S9gDHBhknUT9loOjSvJ\niouLKdaAYqkls2ZV7e55+8aX2OGmUlp9tZCl7ftR2vMhNn6/DcVv1m5bNKZeMhGJRIjk8A6/tEfv\nRLt3psYu5CZ7z8z6A7j7yOh7zwND3P31StvoQq7kVKqpFLZcMHUPnk5VWgqffw4DBsA550CzZonX\nr+Z4+ics+VawCdfMrI27fxH98Y/Au9HXU4CJZjaWoFunAzA32+OI5IKxGZ6aGoT9998HM16ecQY0\n0U3pEi5p/Ys3s0lAN6CVmX0KDAGKzawzQdfNUuAfAO6+wMwmAwuAjcAlKumlYDZtgkcfZR7D4J9N\ng7A/9VRopGmnJJx0c5bUC+lclK3QvbNhAzz4IIwYAa1a8ftXB/Hs5pNSzngZT907UlfpjlwJnWRh\nawZbsZ71t94T3FTVvj0MHgzdumGNLCdPk9ITp6TQFPoSOpVDPxKBV178nq9G3MHAZtezes+DeOm3\ng2h/zpEV/grQPzdpCPTkLAm3b7+l+LVbKb5rHI/yW3ae8zQ7H3QQHQvdLpE6SpW+5F1Nu0jMwL/6\nGm68kQ033sqCopOY3XUAj7zfKem+VOlLQ6HuHanTqgv4jMN4xQpGtxlL7x3uYv7ep/FK1348Pr99\ntR8cCn1pKBT6khe1MX1xbFm189gDfPIJXHcdPPggN60+m8s/uQbatUu638ptLimBIUMyb7NIXaPQ\nl7zL1fTFiZZVWWfJEhg5Ep54Av72N7jqKmzX1qm3idJIG2mIFPqSd3kJ/QULYPhweP55uOQSuOIK\n+MUvUm8jEgI1DX3dlih1Smfehj//GY45BvbbDz7+GP75zy2BLyI1o9CXumHOHDjlFJ7mFDjqqCDs\n+/eHHXYodMtEGhSN05eCiEQgMtMpJsJbLUvZ2z9idtf+/JlHWd9768TrR4LXlacjFpH0qU9fMlbj\nPn13eO45KC1l0WtfM5yBTPj5L9C0aVb7Vp++hIku5EreZRuyjWwzmx97MpjeeONGGDiQxmedzmYa\nb9mfQl8kNYW+5F18yMa6XcrKgu9FRcHr4uLgdXExFHfdCA8/zPtnD+fXh20bTIJ2yinQqFGVB58o\n9EVS09w7UicUFcGyZXDeecGNUOedB8W/+ZmFg+7nm9NGsHb73biKsRx50u8oe9zg8WCbbt2CRxkO\nHar+eZF8UKUvGUs1tbE7bGM/8uP4u1j/r9EsabwPd+wymClrjmbZsiDki4qiHwrFVfeXbtWuG68k\nrNS9IzmVLExbtIA1a6oujw/Z7Wwd6667nc+vGctu3Q8LnlLVpQtQ/fNrY1MkKMBFUlPoS61Ju698\nzRoYP55V/zueXc44lgMfHsA7fmCVfUE1Dy0XkWop9KXWVDenTfPvv6TN5HH86cvb+bDjHzjjnQEs\n8l8lnW4BFPoiNZWXC7lmdjdwMrDK3fePLtsJeBjYEygDerr7muh7A4ALgE3A5e7+YrYNlMKIdfEM\nHZqgy6Xj5xRPHQP33MPtq3vS/OM36bzXXnyY9T9DEcmXtCp9M/stsA64Ly70RwNfuftoM+sHtHT3\n/mbWCZgIHAa0BaYDHd19c6V9qtKv42IV+JZKvKwMRo+Ghx6Cc8+Fa67Bdm9b5SKsKn2R2pOXSt/d\nXzazokqLuwPdoq8nABGgP9ADmOTuG4AyM1sCdAHmZNtIKawOfAjnj4ApU+Dvf4eFC2GXXWq0z1TT\nKujirUjtqck4/dbuvjL6eiXQOvp6NyoG/HKCil9yKB9DFvfjXThrOK8wHYr+J5jbvmXLnOxb4S5S\nGDm5Ocvd3cxS/YGe8L2hsfIOKC4uplgpkLb40DQr/wDIiTffhGHDeJE5cFBvfvnQHawdsn0ODyAi\n6YpEIkRy+D942qN3ot07U+P69BcCxe6+wszaADPdfR8z6w/g7iOj6z0PDHH31yvtT336OZKzPvHZ\ns4N5cd5/H/r1Y5vLLuRH3ybl/mNj7MvKYMKEoKumyjQMxan79EUkfYWchmEK0AsYFf3+ZNzyiWY2\nlqBbpwMwtwbHkSyl1QXkDtOnB2G/fHkwh/2UKdCsGesvS/9YsSkVKh8jEinvr1ffvUjhpTt6ZxLB\nRdtWBP33/ws8BUwG9qDqkM2BBEM2NwJXuPsLCfapSj9H0qn0q6zjDk8/HYT92rUwcCCceSY0aVJl\nm+oq/WTvaaoEkdzTzVmS1jQGW8J50yZ47DEYNgwaNQqmSjjttOA1iYO6pARmzkwc1BpuKZJfCn1J\na8KypraBDfdODB423rIlXHstkea/JzIr+LeTqhLPttIXkdxT6Evq0P/pJ7j3Xj6+eBS/PKYoqOyP\nPbb8ymqCfaTaP6jbRqSQFPqSOPR/+AHuuAOuvx4OOICjnhvEK35UWvsABbtIXaXQlwqBvYN9x3cj\nboVx4+Coo4LK/uCDq+2GUTeNSP2gJ2c1EDWurL/5Bm66iY+4Bd79HW+MmM4zy/aDKRC5Klgl9nQq\nVeoi4aVKvw4yC0bL3HtvcKNTWVmwvKio/KlTEHw4bLtuJW0nj+XUr+5k4a9Opee8ASzx9lX2B6r0\nRRoCde80QJUDOFFovzb5UxqNuY4D5j/A8zv9hbKeffl2xz0oKclshstkxxSRukmh3wClDP2PP4aR\nI+HRR+GCC6BPH2y3NimHbCr0RRoO9emHxD58AOeOgGefhf/+b/jwQ2jVKuU2mr5YRCpTpV8HVai6\n33mHyZ2H0Y1ZtB52BVx6Key4Y9L1k1X6GU/TICJ1krp36qCajsQxA39tTjBVwltv0eeLPtzOxXzv\n2yZdP9U0DAp9kYZDoV/HZRSm7jBrFtOPKeX4PRZDv35wwQXYNltveTubY1T3gPPYa92AJVL3KfTr\nuLRC3x1eeCGY8XLlSi5YMoC7fzobmjXbso/YatkcQ1W8SMOh0K/jYoGbqKo238yfmjzFfk8Ng/Xr\ng7tne/bEmjSudshmomNU1wYRqf8U+nVcosBtbJvYNHFy0Ge/9dYweDB0775leuN0xulXd4xM3heR\n+kNDNgsgVrWXlQXfi4pg4ULYdVdo0QLWrIFTT02w4c8/wwMP8AEj4ZZdgsnQTjgBzFIOr6xJGxPt\nT/31IuGlSj8q24uaiZ4uVeX1j+vh7rth9Gjo0IFu0wcza/PRVaY3rtyeVNMwpDvfvYg0LOreqQWZ\nhGjK0F+3jj7b/5sxbcbAoYcGffaHH57zkFboi4SHuncKJP4vg1jVHZvFckfWwLBb4MYbOYJuwV20\nnTvX2vHVfSMi6apxpW9mZcB3wCZgg7t3MbOdgIeBPan00PS47RpMpR/jX34F48bx1bDbaXXuydC/\nP9Zp36zukBURSaSmlX6jHLTBgWJ3P8jdu0SX9QemuXtHYEb05wZrV77gOq6Gjh3hyy/pwlyYMAH2\n3bfQTRMRqSBX3TuVP3W6A92irycAEepY8Ke6cJvuNj06L+NmRnMWk7iPc2H+fNh9d5beUd7doq4X\nEalLctG98zHwLUH3zr/d/T9mttrdW0bfN+Cb2M9x29WZ7p3Kc9fMmhX8nDSkFy/mro4juXCnJxn5\nzUXcQG9W0TrlpGeVj1dHfnURqWfqwoXco9z9CzPbGZhmZgvj33R3N7OEETc0VgIDxcXFFBewDI41\nJdZHH9e0cu+9B8OHw7RpfMqlsHgxA36xU55aKCJhFIlEiMS6GHIgp0M2zWwIsA64iKCff4WZtQFm\nuvs+ldYtWKVfuZsmvrI/5phgeXzT3vz3W2w7bhjtPn2VB3fpzTdn/DcDR+7AzJnl68dvk6iS1+Rm\nIpILBR2nb2bNgcbuvtbMtgVeBEqA44Gv3X2UmfUHWrh7/0rbphX6mYRlonXjb2xKNu0wVAxsCD4E\nvnrqFa5YV0rrle+y6tyraT/679i2zSuMy68weifN7h0RkWwVOvT3Ap6I/tgEeNDdR0SHbE4G9iCH\nQzazuWkqnWUQH9jOccxgevEwPo4s45f/7g+9esFWW21ZP1bh9+oVfJAsWwatWyeehkGVvIjkUkH7\n9N19KVDlriN3/4ag2s+r+Eofko+aqbwewHm9nAOXP8NrlNKm+bc83mIgPTmL6R2bULxVxXVj+7r3\n3uC7GaxYkZNfQUSkVtWraRhqo9JvZJs5jcd59MBScOf0+YN5ZONp0Lhxyn1U94hCEZHaEKq5d3Ia\n+hs3wqRJLDh3BN+xA0dMHQwnn4w1srSeN1vdIwpFRGpD6EI/3aBNGtjrfwrulh01CvbYg+Mjg5jB\nccTOYboPGVd1LyKFELrQz/oGqB9+4PJt7+Sm3a+DX/86eHBJ164JR+8o9EWkrqoLN2fVbWvXwm23\nwdixHMOR8MQTcOihwcXcocEq8VMlpBK7+Dt0qKZXEJH6qcFW+i1tNauH3gQ33wzHHw+DBmH771ft\nxd10K32T3jVoAAAHfElEQVQRkUJQpV/ZqlVwww0s4Q5Y1gNeeSWY/TKFZHPTi4g0NHWi0k/3rtuU\nVfjy5cEzZ++7D846iz1v7csy37PSMYNt4o9XUlJ+cTjZ8TSFgojUFQ3uQm7lME8WuCUl0fWWLg1G\n4kyezKfHncfkPa5m7fa7bQnzdKZhSOd4CngRqQsafOgn8ytbxKJeI+Dpp+Ef/4Arr4Sdd854P+qj\nF5H6pMH36VeuvM/Y5x1+O3s4LzMT9r4MliwJJrwREZFq1Z9Kf+5cphxeygk7vcGsQ/pw2rSLObTb\ndlWeeJVO14y6cESkvmow3TuxIC4pCUbQxMK3R8uXOOiZUli4kP/5tC83/3AhbLNNlZuqRETCoMGE\nfvlyAMeff5E1fYex+dPPebnrAMavOYcZLzejW7fgomxZWRqPNRQRaWAaVuhv3kyPxlMZTCmHdfoB\nBg2Cnj2hSZOElb0uwopI2DSM0N+0CR55BIYN4//ea0opg3l806nQqFHc+sF3hb6IhFn9Dv2ff4YH\nHoARI6BVK7j2Wuz3JwKWcPoDUOiLSLjV7yGbHTpA+/Zwxx3B1VvL+vcQEZE01Fqlb2YnAuOAxsCd\n7j6q0vvur74KRx5ZZQhlogu0qvRFROpo946ZNQYWETwn9zPgDeAsd/8gbp2ks2ym8wDzZOuJiDRk\nNQ39RtWvkpUuwBJ3L3P3DcBDQI9aOpaIiKSptvr02wKfxv28HDg8mx0lm/ZY4/JFRDJXW907fwJO\ndPeLoj+fDRzu7pfFrZNx944eaCIiYVdXR+98BrSL+7kdQbVfwdC4p5W0aFHMmjXFQOpHEaaq/FX9\ni0hDE4lEiMRCLwdqq9JvQnAh9zjgc2AuGVzITbxPVfUiInVy9A6AmZ1E+ZDNu9x9RKX3qw19zYYp\nIlJRnQ39ag+cxYPRRUTCrq4O2RQRkTpIoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGF\nvohIiCj0RURCRKEvIhIiCn0RkRBR6IuIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIgo9EVEQkShLyIS\nIgp9EZEQyTr0zWyomS03s7ejXyfFvTfAzBab2UIz+11umioiIjVVk0rfgbHuflD06zkAM+sEnAF0\nAk4EbjUz/UWRQiQSKXQT6gydi3I6F+V0LnKnpmGc6InsPYBJ7r7B3cuAJUCXGh6nQdM/6HI6F+V0\nLsrpXOROTUP/MjN7x8zuMrMW0WW7Acvj1lkOtK3hcUREJAdShr6ZTTOzdxN8dQduA/YCOgNfAGNS\n7Mpz12QREcmWudc8j82sCJjq7vubWX8Adx8Zfe95YIi7v15pG30QiIhkwd0Tda2npUm2G5pZG3f/\nIvrjH4F3o6+nABPNbCxBt04HYG7l7WvSaBERyU7WoQ+MMrPOBF03S4F/ALj7AjObDCwANgKXeC7+\nnBARkRrLSfeOiIjUDwUZP29mJ0Zv3FpsZv0K0YZCMbN2ZjbTzN43s/fM7PLo8p2iF84/NLMX40ZD\nNXhm1jh6g9/U6M+hPBdm1sLMHjWzD8xsgZkdHuJzMSD6/8i7ZjbRzLYKy7kws7vNbKWZvRu3LOnv\nnunNsHkPfTNrDNxMcONWJ+AsM9s33+0ooA1Ab3f/NXAEcGn09+8PTHP3jsCM6M9hcQVBd2Dsz86w\nnosbgWfdfV/gAGAhITwX0YEhFwEHu/v+QGPgTMJzLu4hyMd4CX/3bG6GLUSl3wVY4u5l7r4BeIjg\nhq5QcPcV7j4v+nod8AHBBe/uwIToahOAUwvTwvwys92B3wN3Un6zX+jOhZntCPzW3e8GcPeN7v4t\nITwXwHcExVFzM2sCNAc+JyTnwt1fBlZXWpzsd8/4ZthChH5b4NO4n0N781a0ojkIeB1o7e4ro2+t\nBFoXqFn5dgNwDbA5blkYz8VewJdmdo+Z/Z+Z/cfMtiWE58LdvyG47+cTgrBf4+7TCOG5iJPsd8/4\nZthChL6uHANmth3wGHCFu6+Nfy862qnBnyczOwVY5e5vk3hKj9CcC4KRdAcDt7r7wcD3VOq+CMu5\nMLO9gSuBIoJQ287Mzo5fJyznIpE0fveU56UQof8Z0C7u53ZU/KRq8MysKUHg3+/uT0YXrzSzXaPv\ntwFWFap9efQboLuZLQUmAcea2f2E81wsB5a7+xvRnx8l+BBYEcJzcSjwqrt/7e4bgceBIwnnuYhJ\n9v9E5TzdPbosqUKE/ptABzMrMrNmBBchphSgHQVhZgbcBSxw93Fxb00BekVf9wKerLxtQ+PuA929\nnbvvRXCh7v+5+zmE81ysAD41s47RRccD7wNTCdm5ILiAfYSZbRP9/+V4ggv9YTwXMcn+n5gCnGlm\nzcxsL5LcDFuBu+f9CzgJWERw0WFAIdpQqC+gK0H/9Tzg7ejXicBOwHTgQ+BFoEWh25rn89INmBJ9\nHcpzARwIvAG8Q1Dd7hjic9GX4EPvXYILl03Dci4I/ur9HPiZ4Prn+al+d2BgNEsXAidUt3/dnCUi\nEiJ6uImISIgo9EVEQkShLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJkf8PDkxgYpLgjk4A\nAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Important Note:** the way ``curve_fit`` determines the uncertainty is to actually renormalize the errors so that the reduced $\\chi^2$ value is one, so the magnitude of the errors doesn't matter, only the relative errors. In some fields of science (such as astronomy) we do *not* renormalize the errors, so for those cases you can specify ``absolute_sigma=True`` in order to preserve the original errors." ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Exercise 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the following code, we generate some random data points:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.random.uniform(0., 10., 100)\n", "y = np.polyval([1, 2, -3], x) + np.random.normal(0., 10., 100)\n", "e = np.random.uniform(5, 10, 100)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit a line and a parabola to it and overplot the two models on top of the data:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "\n", "# your solution here\n" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 14 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Exercise 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As before, we use the [data/munich_temperatures_average_with_bad_data.txt](data/munich_temperatures_average_with_bad_data.txt) file, which gives the temperature in Munich every day for several years:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# The following code reads in the file and removes bad values\n", "import numpy as np\n", "date, temperature = np.loadtxt('data/munich_temperatures_average_with_bad_data.txt', unpack=True)\n", "keep = np.abs(temperature) < 90\n", "date = date[keep]\n", "temperature = temperature[keep]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 15 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit the following function to the data:\n", "\n", "$$f(t) = a~\\cos{(2\\pi t + b)} + c$$\n", "\n", "where $t$ is the time in years. Make a plot of the data and the best-fit model in the range 2008 to 2012. What are the best-fit values of the parameters? What is the overall average temperature in Munich, and what are the typical daily average values predicted by the model for the coldest and hottest time of year? What is the meaning of the ``b`` parameter, and does its value make sense?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "\n", "# your solution here\n" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 16 } ], "metadata": {} } ] }