Created
March 12, 2014 23:55
-
-
Save PBarmby/9519225 to your computer and use it in GitHub Desktop.
ipython notebook for demonstrating least-squares fitting and chi^2
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "metadata": { | |
| "name": "" | |
| }, | |
| "nbformat": 3, | |
| "nbformat_minor": 0, | |
| "worksheets": [ | |
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "setup and load packages" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "%matplotlib inline" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 1 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "import scipy.stats as sps" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 2 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "import numpy as np" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 3 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "import matplotlib.pylab as plt" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 4 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Generate a set of x values" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "x = np.arange(0,5,0.1)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 37 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Now some Gaussian noise" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "gauss_noise_1 = sps.norm(loc=0,scale=0.3)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 38 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "random_gauss_1 = gauss_noise_1.rvs(len(x))" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 39 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "random_gauss_2 = gauss_noise_1.rvs(len(x))" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 40 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "random_gauss_1" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "metadata": {}, | |
| "output_type": "pyout", | |
| "prompt_number": 41, | |
| "text": [ | |
| "array([-0.01527424, 0.90483845, -0.14772177, 0.48751598, -0.58252677,\n", | |
| " 0.1727146 , -0.07703082, 0.17155032, 0.16298625, 0.30335679,\n", | |
| " -0.19924831, -0.13559424, -0.28148125, -0.17440548, 0.05189903,\n", | |
| " 0.34949444, -0.41032898, -0.44136499, -0.00962964, -0.68953672,\n", | |
| " -0.07359837, 0.12830748, 0.20711083, 0.48446765, 0.11520149,\n", | |
| " 0.1257288 , -0.51378533, 0.35586096, 0.10773178, 0.49506007,\n", | |
| " 0.24545989, 0.24043002, 0.64948792, 0.1179628 , 0.22867986,\n", | |
| " 0.30404989, -0.15652229, 0.25107982, 0.50090187, -0.53797168,\n", | |
| " 0.00661192, -0.12039878, 0.05337592, 0.06687503, 0.04433742,\n", | |
| " -0.46955137, -0.17347405, 0.35325241, -0.32480156, 0.35988961])" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 41 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "random_gauss_2" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "metadata": {}, | |
| "output_type": "pyout", | |
| "prompt_number": 10, | |
| "text": [ | |
| "array([ 0.10161237, -0.0728139 , -0.21675509, 0.07686902, 0.23972285,\n", | |
| " 0.0657801 , -0.10659676, 0.14038602, 0.22224802, 0.05252895,\n", | |
| " -0.2215578 , -0.00246837, 0.03466819, -0.29465311, 0.04524024,\n", | |
| " -0.22324192, 0.1550886 , 0.10407618, 0.18105872, 0.0231177 ])" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 10 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Define y=x+noise" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "y = x+ random_gauss_1" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 42 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "plt.plot(x,y)\n", | |
| "plt.errorbar(x,y,random_gauss_2,ls='None',marker='None',color='k')\n", | |
| "plt.plot(x,y1(x),marker=None,lw=1,ls='solid',label='unweighted',color='k')\n", | |
| "plt.plot(x,y2(x),marker=None,lw=3,ls='dashed',label='weighted')\n", | |
| "plt.plot(x,x,marker=None,lw=1,ls='solid',label='y=x')\n", | |
| "plt.legend(loc='upper left')" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "metadata": {}, | |
| "output_type": "pyout", | |
| "prompt_number": 53, | |
| "text": [ | |
| "<matplotlib.legend.Legend at 0x10799cf50>" | |
| ] | |
| }, | |
| { | |
| "metadata": {}, | |
| "output_type": "display_data", | |
| "png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEACAYAAAC57G0KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdcVfX/wPHXuWxBRBAURVkqmntrLtKGlV/NmeWe3yxH\n0z2w4fpmWmrDEWpmWZJZZlq/ErdirlIRFypJKooDBZmf3x/EFRAQ8FwuXN7Px+M85F7OOZ/3PV7e\n93M/5zM0pZRCCCGExTGYOwAhhBCmIQleCCEslCR4IYSwUJLghRDCQkmCF0IICyUJXgghLJTuCf7a\ntWu8+uqr+Pr6Ymdnh7u7O4899hiHDx/WuyghhBB5sNbzZNeuXaNFixacPXsWa2tratSogbW1NQcP\nHuTs2bM0bNhQz+KEEELkQdcEP2XKFM6ePYuXlxdbt27F398fgLS0NBITE/UsSgghxAPo1kSjlOKb\nb74BwNfXl169euHo6EidOnVYsmQJDg4OehUlhBAiHzS9piq4cuUKlSpVMj728PDA0dGRyMhIABYu\nXMgrr7yiR1FCCCHyQbcEHx0djZeXFwAVKlQgMjISBwcH2rZty549ewgICCA8PPxewZqmR7FCCFHq\n5Ddt69ZE4+7ujo2NDQA1a9bE0dERg8FA48aNATh//nyOQcqmmD59utljKC6bXAu5FnIt8t4KQrcE\nb2Njw2OPPQZAREQEd+7cIS0tjUOHDgHpSV8IIUTR0bUf/LvvvoudnR3Xrl3Dz88Pf39/du/ejaZp\nTJ8+Xc+ihBBCPICuCb5p06aEhobSsWNHEhISiIuLIzAwkN9//51u3brpWZRFCQwMNHcIxYZci3vk\nWtwj16JwdLvJWuCCNa3A7UlCCFHaFSR3ylw0QghhoXQdyaoXV1dXrl+/bu4wRAlTvnx5YmNjzR2G\nEMVGsWyikeYbURjyvhGlgTTRCCGEkAQvhBCWShK8EEJYKEnwQghhoSTBlyJBQUEYDAYMhoL/twcG\nBmIwGIzTURSlFStWGOO+cOFCkZcvREklCb4UqVq1Ki1btqRly5aFPseDZgF9mA+Rhy1bCJFVsewH\nL0xj6NChDB06tEjKkmQshPlJDV4nPj4+GAwGBg8ebHxu0KBBGAwGfH19s+wzcOBApk+fjqenJ+XL\nl6d///7cvn0buNccUaZMGVJSUgDo3bs3BoPBmJzj4uKwtrbGYDCwfv16AOLj45k4cSL+/v7Y2tri\n4eHBsGHDuHbtmjGenGrXSUlJjB49GhcXF9zc3Bg7diyTJk3KtRaulGLZsmX4+vpStmxZ/vOf/3D5\n8mUgvRnn7bffNu6XcY6VK1cCcPXqVV5++WWqVq2Kra0tXl5evPHGGyQkJGQpY8aMGXh4eODk5MSA\nAQO4efPmQ/zPCFGKKTPJq2gzhlVo3t7eStM0NXjwYONzAwcOVJqmKV9f3yz72NraKmdnZ+Xv7680\nTVOapqnJkycrpZQ6f/680jRNGQwGtW/fPqWUUp6enkrTNBUQEKCUUmrz5s1K0zRlZWWlrl27plJT\nU1WbNm2M527YsKEqV66c0jRN1alTRyUkJCillJo+fbrx3BneeustYwy+vr7K3d1dOTk53bdf+/bt\nlaZpqkyZMsrBwUEFBAQYj+vbt69SSqmXX35ZeXl5GZ9v1aqVatWqldq0aZOKi4tTNWvWVJqmKUdH\nR9WwYUNVpkwZpWma6tChg7GcxYsXG4+vUqWK8vLyyhLP+fPnc/0/KInvGyEKqiDv8xKV4IEi3/Kr\nIAm+XLlyKjo6WqWlpammTZsqTdNUy5Ytjcf5+voqTdPU+++/r06fPm08xmAwqCtXrqjJkycrTdNU\ngwYNlFJKff/998YEGBYWppRSKjo6Wjk4OChN09Rnn32mlLo/wd+5c0fZ29srTdPUc889p5RS6vbt\n28ZEnFOCt7a2Vn/++adSSqnu3bsrTdOUp6encb+goKD7jlVKqQULFihN05STk5MxSR8+fNiYzLds\n2aKUUqpatWrG65GcnKxSUlLUY489JgleiH8V5H1eoppoVDFfPSW/OnTogKenJ5qmERAQAKSvaZsh\no6fKrl272LlzJwAjR45EKcXOnTvZsWMHcG8K1bCwMOP1adGiBQaDgSpVqnD37l0A9u3bl2Mcp0+f\nJjExEYDnn38eAEdHRzp37pxr7PXq1aNevXoA1K5dG8DYRJMRQ04yYrxz546xqapRo0bG3+/bt4+4\nuDiioqIAeO6557C2tsbKyoru3bvnGo8QlkDTtDy3wpKbrDrJ+E9ITU01Ppdb27GLi4vxZ2vr9P+C\nzIkxMDCQ4OBgdu3ahZubG/b29owdO5a5c+fy+++/G5Nl9i6LBoOB5s2b31eep6dngV5LXh9sOcVe\nEA4ODjRo0OC+511dXXONwVQftEJYOknwOqlYsSLnz5/n9OnTAMTGxrJt27ZCnSsjccfExLBu3Tqa\nNm1KpUqVqFOnDitXriQxMRGDwUD79u0BjEk9LS2NGTNm8MQTTxgfh4aGGhdDz6569erY29tz9+5d\nQkJC6NOnD7dv32bjxo2FihugTJkyxp/j4+ONj5s3b85XX31FamoqS5YsoW7dugAkJyezefNmmjdv\nTtmyZalatSpRUVH88MMPvPHGG2iaZryRLERJkVut+8GVley/f7jeaCWqiaY469ixIwB79uyhRYsW\n1K9fn1u3bhXqXF5eXvj5+QHp3wJat24NQJs2bYy9berXr2+sTXfp0sW4z1NPPUWtWrWoU6cOzs7O\nPP7445w4cSLHcsqUKcPo0aMBCAkJwdfXFz8/P6KjowsVN9xrtlFK8cgjj9CqVSsiIyMZOnQoNWvW\nJCkpiYYNG1KvXj0CAgJwdnama9euxmaecePGAbB37158fHzw9fVlz549xnMKIfJPErxOJk6cSL9+\n/ShfvjwXLlygf//+9OnTJ8s+ObWn5fZJn1GL1zSNNm3aANC2bVvjc5mXMNM0jS1btjBx4kRq1KjB\nuXPnuHz5MrVr12bixInGGn5OZb377ruMGjUKFxcXbt68Se/evRkyZAgA9vb2BY69c+fODB8+HDc3\nN6KioggLCyMhIQEnJyd27dpl7CZ56tQpbt68SZMmTXj77beNXUlfeeUVpk2bRoUKFbh16xaPPvoo\n7733Xp7XSojiJuMe3qhRo/joo49Mek8vLzIffCl35coV7O3tcXZ2BiAhIYFmzZpx/PhxWrVqxa5d\nu8wcYf7J+0YUN6NHj6ZmzZrGb8q5uVd5ybmJJvP7uiDvc2mDL+V2795Nv379aNasGWXLlmX//v1c\nvnwZGxsb3nnnHXOHJ0Qpo++3VKnBl3J//vknY8aM4dixY9y8eRNXV1datWrFhAkTaNGihbnDKxB5\n34ji4EFNiTm9RwtyTEHe57om+KCgIONQ9exSUlKyDH2XBC/0Ju8bURwUJsHndVz2/c3eROPu7o6/\nv78pTi2EECVEwbo8mqJyYpIE/+yzz/L555+b4tRCCCHyySTdJNetW4eDgwOenp507tyZw4cPm6IY\nIYQQedC1Bq9pGlZWVnh6emJtbU14eDibNm3it99+Y8+ePTRs2DDL/kFBQcafAwMDs/TtFkKI0iQ2\nIZYFexdwPOY463qvMz4fGhpKaGhooc6p603WU6dO4e7ubhxh+csvv9CpUycgfbGJpUuX3itYbrIK\nncn7RhQHBb3Jei3+GvP3zuejfR8RlxQHwKtlX2X+6/NzPb9ZbrLWqFEjy+Mnn3wSV1dXYmNjjbME\nCiGESDd311ze2f4Ot5NuZ3l+X1LOM8AWlK5t8PPnz+eff/4xPv7111+JjY0F0lczEgUnC2ULUbJk\nnm48+1QF2WveSalJWZJ77aSKHLrUjV07y+sSi64J/sMPP8TLywtvb28eeeQRnnrqKQCcnJx49dVX\n9Syq1JCFsoWwPKmpqXzzzTesfX0tVklWNL7rybFznTn2aTINk8qjzftAl3J0baKZPHky3377LceP\nH+fcuXP4+vrSunVrpk6del/zjcgfWShbiJIny9+SE4xZOwZ2Q+LdRFavXs3s2bNxc3NjwUtjabzr\nN1w2bEHr7weHFkO1avoFku+1n3SWV9FmDKvQgoODlaZpysHBQSUnJyullOrVq5fSNE0NGTJEKaXU\nrVu3lJWVldI0TX333Xfqzp07asKECcrPz0/Z2Ngod3d3NXToUHX16lXjeTOW2dM0zfhcYmKiGjVq\nlCpXrpxydXVVY8aMURMnTrxvv4xl9gIDA9XSpUuVj4+PcnJyUp07d1aXLl3Ksk/2bcWKFUoppWJi\nYtTIkSOVl5eXsrGxUVWqVFGvv/66io+Pz/L6g4KClLu7u3J0dFT9+/c3LtH3oGX29FQS3zeWhIdc\n9tKSAAonFE+hmIwiCEUtVNWqVRWgaoNaBSoG1Lug3Auyzqqlrsk6fev09AuVbZu+dbpJ9s9tv5wU\ndLHsK1euyELZOiutyaS4kASf7u+bf6vRm0Yru3fssuSTcuPLqaPBweo7UP+AmgjKuRDXqSD7ynzw\nOqlWrRo+Pj4opdixYwdnzpzh0qVLODs7c+rUKWJiYoxrqdatW5c9e/awa9cuNE1j586dHDp0iPDw\ncOzt7Tl+/DirVq3KsZz4+HgWLlwIQNeuXTl79iyRkZFUrlw519iSkpLYt28fJ06coFu3bgD8/vvv\nACxevJhhw4YB6V8rd+/eze7du3n66adZvnw5p06dwtHRkePHj3Po0CF2794NwNatW/nll18AmDNn\nDgAtWrTg3LlznDt3jmbNmj3sJRUljPr3JuL48eOZNWuW2eZALww910LddGoTC8MWkpiavt4xCgbF\neHNyhw+PTJtGt48+otKdO6zx9uZIZKRJr5MkeB0VZLFsWShbCMsTFRXFkVVHMNwyoCkYFe1N7IYA\nPv/VDo9hY9FOn4bRoyHT0pamVKLmgw8KDCIoMKjY7J9dQRbL3r9/PyALZQtRHGS8X59//nm6d+9u\nrDzl5fyN83iW9cTWypaTJ08yZ84cvv/+e4YPHszvzsNptGYLZR1c0CZOhB49wMrK1C/jPiUqwRd3\nBV0sG2ShbCFKmsjrkczcOZMVh1cwod4ETn51kt9//52xL73EhalTcVy0CDw9Yf7H0KkTmLF3mjTR\n6Kggi2XLQtlClCxnYs8w9Ieh1FxUk2UHl5GSlsK7297lp3Xf0P/qVQa/+y7bXnsNgoNhxw54+mmz\nJneQBK+7/C6WLQtlC1FyHPrnEAGLAvj80OekpKUAUD4Bpm6DswpaAZ2BZwH+/XvPSeYbuOfPn8fX\n1/ehb+rmRZbsK4EsaaFsPcn7pniYMGECLi4uTJgwwdyh5MuDEmtqairr16+n/67+JJRLoFIcfHC8\nKj333MSmWw92tWlD0Fdf8euvvz50Wfl5/5p9RSdhWrJQthBFJ6PpdHqvAdT6eSPPHLjFxzej8Aei\ngoPTm2TImrwfnIALttpTYUkTTQnk5+dH06ZNOXr0KJs3b0YpRdeuXdm+fTsdOnQwd3hClEzux6DH\nC9BmBrDY+HTwm2+yt2ZNxs36lq4tB2Jz8jSvAiVhflypwZdA9evXL/QCAEKURnk2jXgA7YA6dUFT\nUF2D/Z1olggTgZaTJ8PYsbBoEZQrl+3goqmJF5YkeCFE6WQFdAfqZDyhQEHgJcUkmygCEuF/QLez\nZ4tsYJLepIlGCGFSuU0DYKqeI/mWCvzb6UxT0DkCdi9y4dPV1Vhz+3WqA4vARMldy7aZhtTghRAW\nIX8fFgo4A8wFvsUq9Dq942HC/0HqHZiZcoPvuEEaQ0wbbBGRGrwQwrJ5AsYplPoCLbDFleFMISIK\nRm6A8TehcQqsA9IKdPKC1cQzJhZTSuHt7U1kpsnGTNHFV2rwQogiUsQ3JKsA7YGaQApwGhzj1jAc\neIPZ/AkMAnam5D0XTUkerCcJXghhWarshcAZUONn41Plk2FURRgVB6FAF+CQiYovToPtJMELIYpI\nEdWEW4yDGtsBqBQHr+2BoX/ABgO0BU4W8HQZCdvBwYHY2FgcHByAklGzlwQvhCjxUlJS7j3YHoOP\nF7y1G/ocsWW1dRKNkuCCUrreOi1ONfXcSIIXQpQ4mqZBBeBq1udrAxOuhvPMx7DEBmolJxGTnM/z\n/eubb76hT58+QMlI4nmRXjQ6eeGFFzAYDLRs2TLL8+3atcNgMPDCCy8U6HwLFy7EYDBgb29PeHg4\nAKtWrcJgMGBjY5Prik9C6EHPPuuZe4nktBU0hq2RW9Pvjr5CepIHmgLfAb8D4YB/CkxOgJhCRZxz\nDHfv3qVMmTIPfT1ym03SFEyW4Hv37o3BYMBgMNCrVy9TFVNsjBo1CoD9+/dz/PhxAC5dumRcd3Xw\n4MGsWLHCeE1y21auXAnA6NGj6dSpE0lJSQwZMoSoqCjGjh0LwLRp02jRooV5XqgQRSDHhOcLDIYO\nqzqAT/pTXZqX5VqTJuyvVs241ukspbhZgA8SyP1DqKQzSRNNcHAw69atMz7W7dPJHDc18vmf3Lp1\na+rXr8+ff/7J8uXLmTdvHuvXr0cpRZUqVXjiiSf4+eef76vhZ+fh4WH8OTg4mHr16rFv3z6aNWtm\nXERkypQpD/WShHiQjOTm5ubGyZMncXNzM29ALYCn03/UFDx7EiZtB9cbcbx+5wBrgOQxY2DMGBMl\nZv26eBbpB4fS2enTp5WTk5Nq3bq1qlq1qtI0TfXq1eu+/fIq2gRhFYklS5YoTdNUxYoVVXJysurQ\noYPSNE1NmDCh0OfcsGGD0jRNaZqmnJ2d1blz53SM2LKU1PdNcebq6qquXr36wP1Iz4C5bgV179iz\nCkYqnMopq4kG1acH6ogH6kA5VI8yKIMOZeUvDpVt07+sgsSUX7o20aSkpNC3b1+sra1ZvXo1BkPp\nauLv27cv5cqVIyYmhuXLl7Nt2zY0TWPQoEEA/PTTT7Rs2TLPbdOmTVnOGRkZafw5Pj6ev//+uyhf\nkigBit0cLw9B3Ve7bYYtZRl2exonFjgzciuMuw1NbkJIvCItS24X2enaRDNjxgzCwsJYs2YNPj4+\nD9w/KCjI+HNgYKBxObuSqkyZMgwcOJCPPvqIN998k7S0NJo3b05AQAAAV69eJSwsLM8/vqtX73UL\nOHr0qHFVnMaNG3Pw4EEGDBjAkSNHcHJyMu2LEaJQCteUoZRi06lNjNs0jgoH0u+cOgIjeJ03WMyf\n1Gdwwg/sTGinb7glQGhoaOGnB9fra8P+/fuVlZWVGjBggPE5b29vpWma6tmz533751W0jmEVuYiI\nCGUwGIzNKp988kmhznP37l1Vv3594/WLi4tTfn5+StM0NXDgQH2DthAl+X1TWOjcNJJdwZtoCtaU\nkZaWpjac2KAC5gUoglAEoZpPqq6mgroMai29VEMO3ne+wpRVGKa+voWNKd/76lVocHCw0jRNOTg4\nKEdHR+Xo6GhMdNbW1srJyUndunUrX0GW9D/UJ554wngtrl+/XqhzvPbaa0rTNOXh4aFiYmKUUkpt\n375dWVlZKYPBoNatW6dnyBahpL9vCqMkJ/jwK+HKb46fMbFXegM1pzXqqgNquQ2q5gNemyT4B9Ot\nkTyj2SExMZGEhAQSEhKM7WmpqanEx8dbRLej/KhVqxYAXbp0wcXFpVDn+OCDD0hLS+Py5ctUqJD+\nlbVt27akpKSQmppKjx49dItXWILsuafwMrfjx8bGUqFCBV3b9lNTU/n666/p2akn52+ex+c6LP4J\njn0MLd0akPJHGEOT8zOlgOnnVFeZukza29sb81jGVuyZ4APGKKOJprT0ovnss89Uly5dlLW1tbKy\nslL79+83d0ilSkl93zwMHqI2ywNqp7lthTmXUunNjkuXLlXVq1dXrVu3Vts+/VQdebKhiimD+rVP\nc3X57F8mi08P9vb2Kj4+3mTnz6+CvEaTdnMpyXfzC2PPnj38+OOPeHp6smTJEpo2bWrukIR4CDp8\nK9CAOvDy+y9TvXp1vv32W9a+9RY73d1pN306AW26knbyJI9/tQ8P37pZS3/ACFilFO3bt2fr1q0l\nq1ZdhDRlpiuiaVqu/xl5/U6I3JTG9829ClT2153+fH6ux+DBg2nXrh2DBw/W5ZwTJkzAuZwzHh08\nmLR5EjHEUC7emb2+s6kV8h1ERMBbb8HQoYVeDi+3iqMp//+zzyZpLgV5n8tkY0JYhOLxTTk1LZWD\nKQfZdX0X8Zvj09c6PQmTdtyiivW7MPUd6NcPbG3NHWq+Zf4wKZPpA6kkVCakBi8sRml83zyoCTSv\nv7EHK1gN/vz588x5fw6f2nyKoayi9zGYsBPSrAyED3uOJ9/8BLeyHjkeW5yZ49tCXgryPi9dQ02F\nsDCZ255TUlIwGAw6tkfnr5fKiRMnGDx4MI0bN8bF1oH/Ha1BxEJ45YAVh197gWqnrvDCjJASmdwh\n93sBJYHU4IXFKK3vm8LUMB/Uzp4XpRTJqcn8uOtHvlr4Fdu2beP1ESMY4+BAmY8/JtzejuWPuTN5\n4WbKO5TP56sQ+VXi2+DLly9fqnrfCH2ULy/JRC8ZCST7bJJJqUmM+3oci44sIik5ifkNprI6IAC7\nJUsgMBA2bmTl2rVUcHHJV3Ivbs0flqZYJvjY2FhzhyCEyeid1IoiGSamJPLml2+y6PAicIGKKfD6\nXnhxbxBfpsKQEyfg3zmXWLvW5PGI/JE2eCFErpRSbNiwAY9RHiw6twgflT7q9PhicEiARrVhKNxL\n7oU4v1KKU6dO4e/vX6Lat0sCSfBCFLGMJHbx4kU8PT3NnNRyvpGalJREcHAwN2/e5JNPPuFtnz6s\nXA/7l8INA9SqC2OOQtTRf88iTarFUrFsohFC3K+gTTsPk3SrV69OQEAAbe3t2eLhge2CDUzyhNF1\n4dYhIKnQpxZFSBK8EEUop6Sb+bmirMlnlDV48GDqN6nPbzd+Y3PMZjpFPcWd9ZsYFxVFAPDmpk0s\nBxIuG4/MdiapvRdXkuCFKAHyqo0/TPfQ89Hn2XJrCysvrIQyimdj4JOIv7gDzIL0tU4LdeYHK04f\ndpZKErwQZpFzLbgo27K1RzR4FqzqQp9/R52mGOCDZzXmHUlmhbU1KzLvL+3sJY4keCFKFP2aR2zj\nYGAEjN8J0WVh3KOwJQoSPolAsy7K1CBNPqYivWiEKGUcgdeAxLTKDD3vxri+Hhxd9zEbvrmL2q+w\nt7Z/wBlMv9CG0IfU4IWwZHZAc+AwlI+DUaRvoQAbN1K1ekXWOLhhZ21nxiCFqUgNXgizyKsW/HAL\nbSil+P7n76Ed8CpUag5z3eA04AO0BZ4HaNSIymUr5zu553XTU26IFk+S4IUoUXL/YEhLS2P1utV4\n9fOix84e+DSCj39LH3VappIVycf2M0QpIordaFFp8jEVSfBCPITMi1HnZ2HqzNPN6jmS9YsvvqBu\n3brMXjSb8i7RrNiYxh9LINYBnprijeMny3ANaPBQZWSPv6RNnVsaSYIXoohlfAhUqVKFf/7556HX\nLvbx8eHzzz9nxSuv8JdrTXausee4OzwR5EfVhSvZPfE0gxoOwsbKRsdX8fAyf0Bkn4tGPjT0ITdZ\nhSgkUw0+ylczRRnS/3pvwY+vv07dH36AOXPgzTe5smAadW4c5Z26fbA2yJ94aSb/+0IUscLOHXP1\n6lVmfTSLjw99TErDJCbdaEjQfhtYtAgmTIC+fcHWlppAzWoNTRC5KGl0X9Fp2bJlfPrpp5w9e5a4\nuDhcXFx45JFHePnll3n++efvFVxKV98RxVfhJ/Mq2NqlhYnhpTdeIvhkMDRIoVt4KhN3QqoGHu/O\np8qg0WBlVeCy9PZwK0vl/5jSzqwrOu3evZuLFy/i4+NDWloax48fZ8eOHezYsQN/f3+aNm2qd5FC\nWDZrWGMXzAASGfdx+qjT8Y/D34/WYelTLalSDJK7KJ50T/CffPIJdnb3+tUuX76c4cOHo2kaN2/e\n1Ls4IXSTUSuqVq0aO3fupFq1amaJ4ejRo8yaNYstW7YwduhQXnN0JOWjuex2TWRwV4hr3oCp7abS\nrXY3DFrx6SdRmFq31NRNS/cEb2dnx/bt2xk3bhx37twhIiKCChUqMG7cODp27Kh3cUIUWw/qLpld\nWFgY7818j7B9YUwYMYLlPj7YL10KgYHEfP8d70YEMa71OLoEdClWiV0UXya5yXr9+nXCwsKMbUVX\nr17lwIEDJCcnY2Nzr6tWUFCQ8efAwEACAwNNEY4QJlb4Lo5KKbZu3cq096dxpOwRmtSuzAW/F7BZ\ntAi6doUdOyAgAHdgd5sn9QtZlBihoaGEhoYW6ljdb7JmdvnyZd577z0WLVoEwJIlSxg2bFh6wXKT\nVRRT+W2iyV/f9ZxvwKamprJx40amfzCdSK9I3Cvc5rXdqfQ5Cgl9elIlaB6YoYlIFH8FyZ0m/Z5X\nsWJF3nvvPePjw4cPm7I4IYpU5kE5M2bMYOrUqfkepNOgQQNGrB1Bat2/WPjXTfYuSeW6PdQaBXN7\nV5bkLnShaxNNQkICq1evpn///tjbp085+sMPPxh/7+7urmdxQtznQbVqvb81Zi/vnXfeyWXPRGCl\n8dFnw4fjtu5TXDZcZkFLGPUM1Kneii/aT+dJf2mKEfrQtYnmxo0buLq6Ymdnh5+fH4mJiZw9exYA\nV1dXDh06RNWqVdMLliYaYQJ6JPiC9KJ5cDNNHLAUmAfUI5DNTAKeqFqVhNdGU/32THyr1GF6++k8\n7ve4rJokHshs/eAdHBzo168f+/bt48KFCyQnJ+Pt7c1jjz3GlClTjMldCNMrmlWCHjgIyrUsWgN4\nditM4iKuwGzgidOncbC1Ze/NPng5e0liFyaha4K3s7Nj1apVep5SiBLl0qVLzJ8/H9zAqjX0toIJ\nuyHFDmYmwnogDQi2tQWgajmp9AjTkc60QuggMjKSl19+mYBHA/gp7TtGPKoRsQNeOpg+6nToggZ8\nm5ZGqsyUKIqQJHghyDqve1RUFN7e3vma3z08PJwBAwbQtGlTPMqUYUOnQLZ8dpquJxSDnoP2gyH1\nqSdY+MwiaYYRRc6k/eDzLFhusgoTKOxN1oIe98cffzBz5kx27tzJ+OHDeVkpHJYtI7V9O/7jvZuf\ny/7DU/5q2eC9AAAdzUlEQVRPMb39dFpVbVWwFyFEHgqSOyXBC4vy8Ak+95khlVJs376dmTNncuSf\nI0zvMYKhN+KwXbUqfdTp+PEQEMDm05txsXehpVfLh3w1QtxPErwoFfScavZBCX7jxo3MnDmTqKQo\nGnR05unfjzEkogz2g4bBG2/IwCRRZCTBC4uRVxIv/Pztecl8bCoZHc1qtK1Boxa2dP7lGE+fhs+a\nQMgTVdg76Sy2Vrb5OK8Q+ig2UxUIYUoZzSa2trbcvXtXx7U8k4DlQG0AmnrAnJhTfPjpMcLdwX8M\nTOkI3jWbcePuDR3KE8I0JMGLYi0jaffs2ZNvvvmmCBZk1gA7YBiBnOIXIOQKnGxQBb+xMKstPN6o\nO4f+e4j1z6/Hw9HDhLEI8XBkTVZRImVvbsmY+yjDw3wIdAYmgXHU6ZfAt0GL6fzXl0xtN5V6FesV\n+txCFCVJ8EIAVl4w+ckavBx6i5ioWCp++CHuI0YQbGVF8L/7dK3V1awxClFQ0kQjSjiVbSsY26ow\nrCWciIeOv53ilfZ36OTlTsKzzxaLhayFeBhSgxelkqM1jGgIb0TAEXsY3BV2eoNBi6d8ZbsHn0CI\nEkBq8MLiHT58mOeffx6A8sBU4GwKtLoCnV+EZ/vCzqowoMEAwl8Jp8yFMmaNVwi9SIIXFmvXrl08\n++yzPPvss3SoXZvEsWOJdXXl7SFD8Dhxgso7dvJXZSsGNRzEqbGnWPncSmq61TR32ELoRgY6iWLr\nYSbnateuHVF/RzGof3ueO3KS+tvCoX//+0adXrx1kSrOVbIcW5AFP4QoamZb8EOI4qJrhxo8svk4\nTf+3gmVNrfD44yCV/Orft1/25C6EJZEEL0qA3CcAS05OxsHBgerVq1OuXDmGPVOfGhtCqPX+cj5s\nAc+PgVv2qVw6+TkL/BYUfehCmJG0wYsSa/HixdSoUYO01FRWDR7MXmdnui1aQ0il6/iNhdltIaGM\nDSObjuT1Vq+bO1whipy0wYtiK+cZHm8B5dJ/T/qo04lkGnXqBimjNWysbBjeeDjjW48v8LJ4mdvg\n9VjEWwg9yWySwiJkTfBXgQ+BT7DiGr3Lw8TrkALMAkJIX+sUYMWhFTzu93iB29cLc1NX3sOiqEmC\nFxYhe8K1BQZWgvHxEO0C1eZ+hnfv4aDTUnh5J/jc7wMIUZTMNl3wvHnz6NChA1WqVMHOzg4vLy96\n9+7N0aNH9SxG6CjzuqP5WYO0qJw+fdr4syPwakU44wjPlYVBPaDdEJhoHapbcod7M1dm3u7Rsm1C\nFH+61uB9fHyIioqievXqGAwGIiIiAHB0dOTo0aN4e3vfK1hq8MWCnqsi6eGvv/5i1qxZ/PLLL7w+\naBBjra3hs0X8VOUOs9vAIc/0/crYlGF089HM6jjLpB9I0gYvihuz1eCHDRvG6dOniYiIIDw8nHnz\n5gFw584d1q9fr2dRQicZNdWAgADCw8OLYL71rLJ/c6hfvz5bv/qKcdeu8d9583C8coX43zbzQg8D\nhzzB0caR8a3Hc27sOWY/PrtYfNsQorjStR/8lClTsjx+/PHHjT9nn69bWKaCfCPI/pwP8BbQB1gN\nNAIufP457kDA8gACagewdMhSKpSpoG/QQlgokw50+uCDDwCoUKECvXr1uu/3QUFBxp8DAwMJDAw0\nZTiimEhLS+OHH35g5syZANTWYIILPJMAS2yh1i2IyXZMnSt16B3YW5K7KHVCQ0MJDQ0t3MHKBBIT\nE1X//v2VpmnKxcVF7dq16759TFS0KADun0w9y/Yw2rZtq7Zt25blueTkZLV69WpVp04d1bhxY/XL\n7Jnqu3Kof8qgJnREOU9AMaCjAnVfDD179lTffPPNQ8UkhCUoyN+m7jX4q1ev0q1bN3bt2kXlypX5\n6aefaNCggd7FiBLk7t27rFy5kjlz5lCtalWCBwyg1sbviJs9lVmPQt/GkGDz786eB6BsNMSZNWQh\nLIKuCT48PJzOnTsTGRlJo0aN+PHHH6lcubKeRQiTyLmPd0Flb39v3759ljMenTOHR9avh+XLSXjz\nVRrHnORy8vX0HRLKw57XYN8YSCxXqPKFEFnp2k2yVq1anDx5EoC6detSpsy9hROGDx/O0KFD7xUs\n3STNLuepAKCwg3hyusFqBfQifTqB+o0awaRJ0K0bWFnx3vb3mPLTFNgDhAGJ+S9L3juitDLbSFZf\nX18uXLiQ4++mT5/OtGnT7hUsCd7sTJfgXbGlCwM0b8bbfky0FsPMu7A5LS3LwKQ7SXdwKusESQWP\nXd47orQy23zwkZGRep5OlBCRkZHMnTsXSB91OtzwX96w+oQ/qyQw+LFEdtoBn3HfqFNHW0dUoiRq\nIUxFpgsWFHYY/rFjx+jfvz9NmzalioMDUzU4aw2tqs/nP0Nu8OygRHZ6A5WAANNELoTInSR4UWD7\n9++ne/fudOzYkWZeXkT368eUlSvxUdD2RXj+xbsc/ndKAW5XhC3AGXNGLETpJAm+FFOZJtXKPlVB\n9jY+pRRbt27liSeeoEePHnStX5+oLl0Y89ln2AEcOsRQ4GTGvHJxwGbgw8vpN1GTi/KVCSFAluwT\nD6CUYuPGjcycOZPY2Fhm9+/PMxVcsfpoAdYjX4ETJ8DD494BR0hv5TlC+mTtQgizkRp8KVCYKYFT\nUlL46quvaNCgAdOmTePtLl04GlCDJ+bP4X83NtLgTSeS3p6eJbkrpVApCvWHom2rtmzbti3XbwRC\nCNOTBT9KgbwmAMvtd/7+/lSqWJF5//kPjX/dQvyxw8xqkcRH9eKNo04/6/wZI5qMyHdZQoiHZ7bp\ngkXxlFGDfuutt5gzZ84Da9QasGHYMHampdEiOJgFfldwH3GDOY3vJXfvct642LsUzQsQQhSK1OAt\n3IOaYq5fv86iRYv46KOPCGzThtlNmuC3di1YWxtHnW44tZHn1j4HgI+LD5PbTmZAgwHYWtkWxUsQ\nQmQia7IKowcleFdXV7o98wzvVK9OpVWr0CpXhsmT4amnjAOTlFL0+KYHnWt2pn/9/thY2eR5TiGE\n6UiCL+H0bMe+fzqC88D/gMU4AhcmT8bl8+Wc8XJkUovbLPtfBOXsZbIvIYoraYMXOTgBDAIaUx4r\npgJngYs7Q2jf/RY1nz3DugqXWRi2ENB/Me7iuri3EJZMErwJFTapZdwEHTBgACtWrNCpm2E7KlKR\nOfTlFKvxKQ9tR0D9x06w0z3euNemU5seWJYkayFKBknwFmrHjh08/fTTQPpap4t5juMsxQFFIw4x\n1BFOZpqqv3aF2qzpvoYdg3cYvwIqpYiIiKBGjRoP/SGTcfyNGzdwdnaWvvFCFAFJ8CaUkcSSkpKw\ntrY2eVJTSvHzzz/Ttm1bBg8ezJCWLVkJ7Aeus5RaXGcMi4jCG/4GzkAd9zp83eNr/hr5Fy/UewEr\ng1WeZeRVU5davBDFi0xVYAFSU1P57rvvmDlzJikpKfyvd2+ePHAAtXgRU6rD6Etw63YOB4bAnyv/\nxKCZ9nM+p8Sf+TmpyQthGpLgTaQoklpSUhJffvkls2fPpryLC4t79qRVaChpSz5l/X9qMHxkHNcN\nwD7g5/uPV3dyjuFBsf97dPY9CvMShBAmJE00JVB8fDwLFy6kevXqfPXll3w7YAB7DAaarwxmTT2F\ny9Cr9Ky4jeuGf9fAawI4mTXkf6lsmxDClCTBm5x+Se3mzZvMmjULPz8/tv76K78PH84vly+TMmUK\nvfbuxe7aGfqV/Y3bWqY18KJhQ78NpN1KyzLxV/6+QUhCFqIkkwRfAsTExDB58mT8/f05+ddfHHzp\nJb47dozqv/4K//sfTYAQIC0WyJiP/SLwJbAEugR0kRugQpRCMpI1nwo6ulSPBa179uxJXFwc+/fv\nZ0D37kytWBG3FStIqFOTG6+NxPPp3veX5XoK3E7BqafJ+Pwu/OLZBW9nz6ksvRf3FqI0k5GsJUBe\ng4VOnTrFsGHD+PHHH3FOSeH8kCEs2LABu2OHefut5ji33clrd77L+cSxNeDUMxTPm56FW/tVCFE4\nuib47du307lzZypWrIjBYMBgMDBjxowHHpdbsitOIyQz2q0XL17MyJEjC9COXbCk1qdPHx599FFq\nu7iwzteX1fv2oS6dZ9Lb7Snf5Bem3/yelLQUvjn2Dcdjjj/kq9I3dqmJC1G86NpN8uDBg2zZsoWa\nNWsSExMDWP7gl8K+voxkGBISwsKFC3FycuLgwYPsW7uWt4E+8+bxBVDDBqKqrYMr2Y4/p6hTrw5c\nfrj4C6OgiTzz/jdv3qRatWrcvHlT77CEENnoWoMfMGAAcXFxhIWF6XnaEi2nkawZ2y+//ML06dPZ\nv38//Rs3JiowkD+AWCAAGAtEJQMHM50wEljx73ZfctenCSRzjNmnKpBauhAlh641eFdXVwBu385p\n2GR+WN7gmcw1/IyfQ0JCmDlzJgkJCQypV49HY2NptWQJqWNG4wfcArJci13/QIXKsJP02X5LoOzf\ndDIeyweGEKZj1pGsQUFB2Z4pngk9p2aYTz75pNDnmzVzJgu6dqX19u0k/PYbwdVc+eyNBkSnbvs3\nuWdz2zO9y+MDSLIUwvKEhoYSGhpaqGNN0k3y9u3bODs7A+lJfNq0afcXnKmrz4PasfUOsfBdHvOS\n87ePhIQEgoODmTt3Ln4+Pszr0IEGmzahxcby98v96Zv0E9vv7L332fY5cCH38+VVlp7XSRbPFqJ4\nKkg3SZPU4C33xmpeSTfn1+zn50fzxo3ZMngwNdetg/Xr6XnoEOufgbQbU+8/1J9/E7wQQjwck/SD\nz/zpYqoa38N0rcy4WdioUSMOHDhgfKx3N01bIGzYML4PD6fmr7/C3Llw4ED6qNPsbTEngaXAVuMr\nRI8bpoWV/Yaw3GAVouTRNcF/9913VK9enQYNGhif++ijj6hevTr9+vXL9bjckklJTSqOwGukL4nn\n9ccfsGIF7NgBnToZF7ImDIh3g4j/wJL9sEalTy8ghBA60TXBx8XFcfbsWc6dO2es/d64cYPIyEii\no6P1LMokMj5QoqKiqFKlSr4/YEaMGEH58uWZ9NJL3HjjDW67u/NBr17c3LKasWNqkNam9f0HJQEL\nI+CrHyC6aY5x5L5kn4wIFUI8mK5t8AMHDmTgwIF6njIfCta1MnuTS5MmTbKeLc+Enks7u4MDf7/4\nImXWrIHnnuOvdR8z8e8V/LQn/VtLu2rt6PFIj/sPTHDLM1YhhHgYMhfNQ/ABFgPjV62ijJUVf275\ngmc6XqL+1l78dOon437v7nhXl6amvM5REpuyhBCmZQEJvrDNFfmf6zwtLY3ffvuNjh074urqSt/G\njUl+8UUi3dx4edIkOHECPvyQMKt/+Pn0vaWTNDSer/M8X3T7IpebtQWP3VLuUwghTE+W7MtDWloa\nGzduZObMmdy4cYP3n3+e2jExVAgPx7pnT/j4YyhXzrj/gAYDeHf7u1y4eYE+dfswpd0UHnF/xIyv\nQAhRmkmCz0FKSgrx8fE0aNAAWxsb5nftStsdO9CCg9nRsiWTWldj1VuvY2dtl+U4WytbPu/6OZXL\nVqZWhVq5nj9zjTskJIQ1a9YQEhJistcjhCidSmyCN2X/ekhvMGl9/TqTANdDh9A+/5xtrb0Yvu4V\nTiWf4rHDwbzU9KX7ju/g28EkcQkhREGV2AT/8HJu87YCegMTgBRgJrDeB9K2D8kywnTklyMZ2XIk\nKkXav4UQxZNZE3xxmlHQFhgIjMePi1RhPBPZTCfwNsCgbDunkj6KyaaooxRCiPyzgF40BaOUIjo6\nmrfeegsrKyt6P/MMl8eP5yzwHDCIFbRnO5t5GtDSp+f9d4yWtcGa4Y2Hc/a1s6gNCpVg/g8mIYTI\njVlr8EVdcz937hxz587l66+/ZkTPnixwdeW/e/Zg4+REI+CwAUhre/+BofDfef9lYpuJeLt46xKL\nzI8uhDA1k0wXnK+CCzDl5cMKDw9n9uzZbNy4kTf79mVMaiqOX3/NBsB/yRLO17On8/ud0yf82mb6\nqXgh7+l4ZapeIURuCpI7LbqJ5o8//qBHjx4EBgbS3N2d6OeeY+Lq1ThaW6MOHuT15i48HzWNzl91\nhipAS8CuaOZ5kQFLQghTs7heNEoptm/fzsyZMzl+/Dgz+/Zlra0t1sHB8NJLcOIEWtWKkPhRekLP\nvPazDeAFnDFP7Bkk0Qsh9GAxNXilFD/99BNt2rRh2LBhvNysGecbN6Z/cDDWdevCmTPw3nvg4ZE+\nk+OdTAcnA3vg4riLqNP3zyYpCVcIURKV+Bp8amoq69atY9asWai0NBY89xyBu3ejrVoFb74JX30F\nZcpkOUYpxe6o3Ty+6nFeavoS41qPo5JTJTO9AiGEMI0Se5M1KSmJL774gjlz5lDBzY0FTzxBs19/\nRYuNhQkTSH2hD+tO/8CFmxd4q/VbOZ7jesJ1yjuUv+/5v//+m5YtW/L3338XOj4hhDAFi77JeufO\nHT788EP8/f1Z9/XXfN+nD7vu3KH5jz+ivf46qUf/Yk1TO+otb0KfkD5M2TqFi7dyXiopp+QuhBCW\nosQk+Bs3bvDee+/h5+fH7q1b2TFgAD9HRvLI1q1oc+fCwYN8VSuFOp/Vp+93fQm/Gg5AUmoS8/bM\nM3P0QghR9Ip9G/zly5dZsGABS5YsofuTT3Jk0CAqrV4NSUkQHAxt7w1M2nR6ExHXIoyPy9qWZWzL\nsbzW8jVzhC6EEGZVbGvwFy5cYMyYMdSuXZuUK1c4PXAgS3/7jUqRkbBxI2zalCW5A0xpOwWDZsDZ\nzplp7adx/tXzvPPYO7g6uJrpVQghhPkUu5usERERzJkzhw0bNvBanz7YfvwxQ4ENwBzgpAHwBiLT\n989+jpDjIXTw7VCo9nUZQSqEKO4KcpO12CT4Q4cOMWvWLEJDQ5nSty8jbt3Cfv16Prx+nXlAlBXQ\nAGgLuAAfAzH6Jl9J8EKI4s6svWi+/vprGjdujIODA66urvTq1YszZ3IeGqppmnHr3Lkzz/r4cLFj\nR8Z88QX2lSrBiRO8agVRTYDRQBegPOkzCLTTO3KZPkAIYVl0rcEvX76c4cOHA+Dn58e1a9e4efMm\nHh4eHDlyhIoVK94rOFNtuSmwt0sXrPbtg7Fj4eWXjWudao9pEJitoHg32HUNdkntWghRupilBp+U\nlMSECRMA6NmzJ6dPn+b48eOULVuWK1euMHPmzPuOCWQyvwDrAKvHH4ezZ2HixCwLWfMH6VMJANyp\nAL/OgQXnYFfe8WT+dpB5E0KI0kK3BL9//36uXbsGQI8ePQDw9PSkZcuWAGzevPm+Yz7lW74EagCJ\nI0fcN6UAALeBHcAvwIKrsGs8JJXVK2whhLBYuvWDj4qKAtJrzh4eHsbnM37O+H1mtemDMrwNlaHy\nmMqsH7+edt45NK5vL3g8GV9hgoKCsvwrhBAlSWhoKKGhoYU61uQDnfJqK1KtnKE14ASxxDJj2wx+\nG/CbqUMSQogSIzAwkMDAQOPjGTNm5PtY3RJ8tWrVgPSEfvnyZePzV65cyfL7LJ56M8vDE1dPcOXO\nFTwc730DKMxN1Jza2jNfFLkxK4QoDXRrg2/WrBlubm4AhISEABAdHc3evXsB6NSpU+4H34RFTy/i\nzJgzWZK7EEKIwtO1m+TSpUv573//C4CPjw/Xrl0jLi4Od3d3jhw5QqVK9+Zc1zSNqh9UZVLbSQxu\nOBg7azu9wshUgy+a9VWFEKKomHUk65o1a3j//fc5ceIE9vb2dOjQgdmzZ1O9evX7gkxMScTWylbP\n4o3nTicJXghhWUrkVAV6nzudJHghhGWx6AU/hBBC5E+xnw/+4cjIVSFE6SU1eCGEsFAWWYPP3D4l\nI1mFEKWVRd5kzTh/TuQGqxCiJJObrEIIISyziQakpi6EEFKDF0IICyUJXgghLJQkeCGEsFCS4IUQ\nwkKZNcHLOqlCCGE6UoMXQggLZdZuktKVUQghTEdq8EIIYaEkwQshhIWSBC+EEBZKErwQQlgoSfBC\nCGGhJMELIYSF0i3BL1myhPbt21O2bFkMBgMGg4Ft27bpdXohhBAFpFuC//nnnzlw4ACVKlUCZJRq\nQYSGhpo7hGJDrsU9ci3ukWtROLol+I8//phbt26xYMECvU5Zasib9x65FvfItbhHrkXh6DaS1dPT\nE5DRqUIIUVzITVYhhLBUKg+TJ09WmqbluYWGhmY55scff1SapimDwaC2bduW67kB2WSTTTbZCrHl\nV55NNE2aNGHQoEF57WJsmikoacoRQgjTyjPBd+vWjW7duhX65JLEhRDCfDSlUxYeP348ISEhxMfH\nc+nSJSC9du/g4MDYsWMZPXq0HsUIIYTIJ91usl65coXIyEguX75s7AN/6dIlIiMjuX79unG/r7/+\nmsaNG+Pg4ICrqyu9evXizJkzeoVRYmzfvp3OnTtTsWJF48CwGTNmmDusIjdv3jw6dOhAlSpVsLOz\nw8vLi969e3P06FFzh2YWy5Yto2nTpri6umJjY4O7uzvt27dn7dq15g7NrHr37m38O+nVq5e5wylS\nQUFBxteefUtLS8vzWN26SQYHBxMcHJznPsuXL2f48OEA+Pn5ce3aNUJCQtixYwdHjhyhYsWKeoVT\n7B08eJAtW7ZQs2ZNYmJiAErlwLCFCxcSFRVF9erVcXZ2JiIignXr1vHzzz9z9OhRvL29zR1ikdq9\nezcXL17Ex8eHtLQ0jh8/zo4dO9ixYwf+/v40bdrU3CEWueDgYNatW2d8XBr/TgDc3d3x9/cv2EH5\nvh37kBITE1WFChWUpmmqV69eSimloqOjlbOzs9I0TY0ZM6aoQikWrl27phISEtTt27eNPZJmzJhh\n7rCK3DvvvKPOnj1rfPzBBx8Yr8f8+fPNGJl53L17N8vjZcuWGXul/d///Z+ZojKf06dPKycnJ9W6\ndWtVtWrVLPmjtJg+fbrSNE0NHjy4wMcWWT/4/fv3c+3aNQB69OgBpLfRt2zZEoDNmzcXVSjFgqur\nK/b29qX+RvSUKVPw9fU1Pn788ceNP9vb25sjJLOys7Nj+/bttGzZknr16jFy5EgqVKjAnDlz6Nix\no7nDK1IpKSn07dsXa2trVq9ejcFQuoftrFu3DgcHBzw9PencuTOHDx9+4DFFdsWioqKA9K9XHh4e\nxuczfs74vSjdPvjgAwAqVKhQ6tpaM1y/fp2wsDCOHz9OSkoKV69e5cCBAyQnJ5s7tCI1Y8YMwsLC\n+OSTT/Dx8TF3OGajaRpWVlZ4enri5+fH5cuX2bRpE61atXpgkjf7R2Jpr8GKdElJSQwYMICVK1dS\nrlw5vv/+e9zc3Mwdlll07dqVtLQ0oqOjGTVqFABr165l5cqVZo6s6Pzxxx/MmjWL/v3706dPnyy/\nK20548UXXyQmJoaIiAiOHTtmbO1ITExk8eLFeR5bZAm+WrVqQPp/zuXLl43PX7lyJcvvS5vSesMo\ns6tXr9KxY0dWr15N5cqVCQ0N5dFHHzV3WGZXsWJF3nvvPePj/HwltxRHjx4lLS2Nb7/9FicnJ5yc\nnIzf8r///nvKli1LXFycmaMsGjVq1MDFxcX4+Mknn8TV1RV4cMtHkSX4Zs2aGWtkISEhAERHR7N3\n714AOnXqVFShFCuZayOlrWYCEB4eTosWLdi1axeNGjUiLCyMBg0amDsss0hISGDp0qXcvXvX+NwP\nP/xg/Nnd3d0cYZlFRsUnMTGRhIQEEhISjH8fqampxMfHl5q/l/nz5/PPP/8YH//666/ExsYCPLjp\nSsebvQ+0ZMkSYw8JX19fYw8aDw8P9c8//xRlKGYXEhKi/P39lZ+fn/GauLq6Kn9/f9W3b19zh1dk\nAgICjK+/Xr16qkWLFsZt2bJl5g6vSF2/fl1pmqbs7e3VI488ovz9/Y3Xxs3NTV24cMHcIZqVt7d3\nqexF4+3trQwGg6pWrZqqXbu28T1RtmxZFR4enuexRdoGP3z4cFavXk3Dhg25dOkSVlZWdO/enV27\ndhkXCikt4uLiOHv2LOfOnTMODLtx4waRkZFER0ebO7wik5iYaHz9x44dY//+/cbt4sWL5g6vSDk4\nONCvXz+qVavGhQsX+Pvvv/H29mbQoEGEhYVRtWpVc4doVqV1EaHJkyfTsWNHUlNTOXfuHL6+vvTr\n148DBw5Qq1atPI/VbaoCIYQQxYvZe9EIIYQwDUnwQghhoSTBCyGEhZIEL4QQFkoSvBBCWChJ8EII\nYaH+H5/CkCfv42kIAAAAAElFTkSuQmCC\n", | |
| "text": [ | |
| "<matplotlib.figure.Figure at 0x1077ec0d0>" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 53 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Do a least-squares linear fit" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "sps.linregress(x,y)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "metadata": {}, | |
| "output_type": "pyout", | |
| "prompt_number": 45, | |
| "text": [ | |
| "(1.0012842578118442,\n", | |
| " 0.047093021680977731,\n", | |
| " 0.97446382195943138,\n", | |
| " 8.5576353434068425e-33,\n", | |
| " 0.033302250906870397)" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 45 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Look up what the fit parameters mean" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "help(sps.linregress)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "Help on function linregress in module scipy.stats.stats:\n", | |
| "\n", | |
| "linregress(x, y=None)\n", | |
| " Calculate a regression line\n", | |
| " \n", | |
| " This computes a least-squares regression for two sets of measurements.\n", | |
| " \n", | |
| " Parameters\n", | |
| " ----------\n", | |
| " x, y : array_like\n", | |
| " two sets of measurements. Both arrays should have the same length.\n", | |
| " If only x is given (and y=None), then it must be a two-dimensional\n", | |
| " array where one dimension has length 2. The two sets of measurements\n", | |
| " are then found by splitting the array along the length-2 dimension.\n", | |
| " \n", | |
| " Returns\n", | |
| " -------\n", | |
| " slope : float\n", | |
| " slope of the regression line\n", | |
| " intercept : float\n", | |
| " intercept of the regression line\n", | |
| " r-value : float\n", | |
| " correlation coefficient\n", | |
| " p-value : float\n", | |
| " two-sided p-value for a hypothesis test whose null hypothesis is\n", | |
| " that the slope is zero.\n", | |
| " stderr : float\n", | |
| " Standard error of the estimate\n", | |
| " \n", | |
| " \n", | |
| " Examples\n", | |
| " --------\n", | |
| " >>> from scipy import stats\n", | |
| " >>> import numpy as np\n", | |
| " >>> x = np.random.random(10)\n", | |
| " >>> y = np.random.random(10)\n", | |
| " >>> slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)\n", | |
| " \n", | |
| " # To get coefficient of determination (r_squared)\n", | |
| " \n", | |
| " >>> print \"r-squared:\", r_value**2\n", | |
| " r-squared: 0.15286643777\n", | |
| "\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 16 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Define a function which returns the fitted value for any x" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "def y1(x):\n", | |
| " return(0.0471+x*1.0013)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 46 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Polyfit is another way to do these fits, can also weight by uncertainties" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "np.polyfit(x,y,1,w=1.0/random_gauss_2**2,cov=True)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "metadata": {}, | |
| "output_type": "pyout", | |
| "prompt_number": 49, | |
| "text": [ | |
| "(array([ 1.04069318, -0.1239218 ]),\n", | |
| " array([[ 0.00442845, -0.01308322],\n", | |
| " [-0.01308322, 0.04038001]]))" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 49 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "help(np.polyfit)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "Help on function polyfit in module numpy.lib.polynomial:\n", | |
| "\n", | |
| "polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)\n", | |
| " Least squares polynomial fit.\n", | |
| " \n", | |
| " Fit a polynomial ``p(x) = p[0] * x**deg + ... + p[deg]`` of degree `deg`\n", | |
| " to points `(x, y)`. Returns a vector of coefficients `p` that minimises\n", | |
| " the squared error.\n", | |
| " \n", | |
| " Parameters\n", | |
| " ----------\n", | |
| " x : array_like, shape (M,)\n", | |
| " x-coordinates of the M sample points ``(x[i], y[i])``.\n", | |
| " y : array_like, shape (M,) or (M, K)\n", | |
| " y-coordinates of the sample points. Several data sets of sample\n", | |
| " points sharing the same x-coordinates can be fitted at once by\n", | |
| " passing in a 2D-array that contains one dataset per column.\n", | |
| " deg : int\n", | |
| " Degree of the fitting polynomial\n", | |
| " rcond : float, optional\n", | |
| " Relative condition number of the fit. Singular values smaller than this\n", | |
| " relative to the largest singular value will be ignored. The default\n", | |
| " value is len(x)*eps, where eps is the relative precision of the float\n", | |
| " type, about 2e-16 in most cases.\n", | |
| " full : bool, optional\n", | |
| " Switch determining nature of return value. When it is\n", | |
| " False (the default) just the coefficients are returned, when True\n", | |
| " diagnostic information from the singular value decomposition is also\n", | |
| " returned.\n", | |
| " w : array_like, shape (M,), optional\n", | |
| " weights to apply to the y-coordinates of the sample points.\n", | |
| " cov : bool, optional\n", | |
| " Return the estimate and the covariance matrix of the estimate\n", | |
| " If full is True, then cov is not returned.\n", | |
| " \n", | |
| " Returns\n", | |
| " -------\n", | |
| " p : ndarray, shape (M,) or (M, K)\n", | |
| " Polynomial coefficients, highest power first.\n", | |
| " If `y` was 2-D, the coefficients for `k`-th data set are in ``p[:,k]``.\n", | |
| " \n", | |
| " residuals, rank, singular_values, rcond : present only if `full` = True\n", | |
| " Residuals of the least-squares fit, the effective rank of the scaled\n", | |
| " Vandermonde coefficient matrix, its singular values, and the specified\n", | |
| " value of `rcond`. For more details, see `linalg.lstsq`.\n", | |
| " \n", | |
| " V : ndaray, shape (M,M) or (M,M,K) : present only if `full` = False and `cov`=True\n", | |
| " The covariance matrix of the polynomial coefficient estimates. The diagonal\n", | |
| " of this matrix are the variance estimates for each coefficient. If y is a 2-d\n", | |
| " array, then the covariance matrix for the `k`-th data set are in ``V[:,:,k]``\n", | |
| " \n", | |
| " \n", | |
| " Warns\n", | |
| " -----\n", | |
| " RankWarning\n", | |
| " The rank of the coefficient matrix in the least-squares fit is\n", | |
| " deficient. The warning is only raised if `full` = False.\n", | |
| " \n", | |
| " The warnings can be turned off by\n", | |
| " \n", | |
| " >>> import warnings\n", | |
| " >>> warnings.simplefilter('ignore', np.RankWarning)\n", | |
| " \n", | |
| " See Also\n", | |
| " --------\n", | |
| " polyval : Computes polynomial values.\n", | |
| " linalg.lstsq : Computes a least-squares fit.\n", | |
| " scipy.interpolate.UnivariateSpline : Computes spline fits.\n", | |
| " \n", | |
| " Notes\n", | |
| " -----\n", | |
| " The solution minimizes the squared error\n", | |
| " \n", | |
| " .. math ::\n", | |
| " E = \\sum_{j=0}^k |p(x_j) - y_j|^2\n", | |
| " \n", | |
| " in the equations::\n", | |
| " \n", | |
| " x[0]**n * p[n] + ... + x[0] * p[1] + p[0] = y[0]\n", | |
| " x[1]**n * p[n] + ... + x[1] * p[1] + p[0] = y[1]\n", | |
| " ...\n", | |
| " x[k]**n * p[n] + ... + x[k] * p[1] + p[0] = y[k]\n", | |
| " \n", | |
| " The coefficient matrix of the coefficients `p` is a Vandermonde matrix.\n", | |
| " \n", | |
| " `polyfit` issues a `RankWarning` when the least-squares fit is badly\n", | |
| " conditioned. This implies that the best fit is not well-defined due\n", | |
| " to numerical error. The results may be improved by lowering the polynomial\n", | |
| " degree or by replacing `x` by `x` - `x`.mean(). The `rcond` parameter\n", | |
| " can also be set to a value smaller than its default, but the resulting\n", | |
| " fit may be spurious: including contributions from the small singular\n", | |
| " values can add numerical noise to the result.\n", | |
| " \n", | |
| " Note that fitting polynomial coefficients is inherently badly conditioned\n", | |
| " when the degree of the polynomial is large or the interval of sample points\n", | |
| " is badly centered. The quality of the fit should always be checked in these\n", | |
| " cases. When polynomial fits are not satisfactory, splines may be a good\n", | |
| " alternative.\n", | |
| " \n", | |
| " References\n", | |
| " ----------\n", | |
| " .. [1] Wikipedia, \"Curve fitting\",\n", | |
| " http://en.wikipedia.org/wiki/Curve_fitting\n", | |
| " .. [2] Wikipedia, \"Polynomial interpolation\",\n", | |
| " http://en.wikipedia.org/wiki/Polynomial_interpolation\n", | |
| " \n", | |
| " Examples\n", | |
| " --------\n", | |
| " >>> x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])\n", | |
| " >>> y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])\n", | |
| " >>> z = np.polyfit(x, y, 3)\n", | |
| " >>> z\n", | |
| " array([ 0.08703704, -0.81349206, 1.69312169, -0.03968254])\n", | |
| " \n", | |
| " It is convenient to use `poly1d` objects for dealing with polynomials:\n", | |
| " \n", | |
| " >>> p = np.poly1d(z)\n", | |
| " >>> p(0.5)\n", | |
| " 0.6143849206349179\n", | |
| " >>> p(3.5)\n", | |
| " -0.34732142857143039\n", | |
| " >>> p(10)\n", | |
| " 22.579365079365115\n", | |
| " \n", | |
| " High-order polynomials may oscillate wildly:\n", | |
| " \n", | |
| " >>> p30 = np.poly1d(np.polyfit(x, y, 30))\n", | |
| " /... RankWarning: Polyfit may be poorly conditioned...\n", | |
| " >>> p30(4)\n", | |
| " -0.80000000000000204\n", | |
| " >>> p30(5)\n", | |
| " -0.99999999999999445\n", | |
| " >>> p30(4.5)\n", | |
| " -0.10547061179440398\n", | |
| " \n", | |
| " Illustration:\n", | |
| " \n", | |
| " >>> import matplotlib.pyplot as plt\n", | |
| " >>> xp = np.linspace(-2, 6, 100)\n", | |
| " >>> plt.plot(x, y, '.', xp, p(xp), '-', xp, p30(xp), '--')\n", | |
| " [<matplotlib.lines.Line2D object at 0x...>, <matplotlib.lines.Line2D object at 0x...>, <matplotlib.lines.Line2D object at 0x...>]\n", | |
| " >>> plt.ylim(-2,2)\n", | |
| " (-2, 2)\n", | |
| " >>> plt.show()\n", | |
| "\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 19 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "def y2(x):\n", | |
| " return(1.0407*x-0.124)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 50 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Does scipy.stats already have a chi^2 function? Maybe, but chi here is the chi^2 probability distribution" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "help(sps.chi)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "Help on chi_gen in module scipy.stats.distributions object:\n", | |
| "\n", | |
| "class chi_gen(rv_continuous)\n", | |
| " | A chi continuous random variable.\n", | |
| " | \n", | |
| " | %(before_notes)s\n", | |
| " | \n", | |
| " | Notes\n", | |
| " | -----\n", | |
| " | The probability density function for `chi` is::\n", | |
| " | \n", | |
| " | chi.pdf(x,df) = x**(df-1) * exp(-x**2/2) / (2**(df/2-1) * gamma(df/2))\n", | |
| " | \n", | |
| " | for ``x > 0``.\n", | |
| " | \n", | |
| " | %(example)s\n", | |
| " | \n", | |
| " | Method resolution order:\n", | |
| " | chi_gen\n", | |
| " | rv_continuous\n", | |
| " | rv_generic\n", | |
| " | __builtin__.object\n", | |
| " | \n", | |
| " | Methods inherited from rv_continuous:\n", | |
| " | \n", | |
| " | __call__(self, *args, **kwds)\n", | |
| " | \n", | |
| " | __init__(self, momtype=1, a=None, b=None, xa=None, xb=None, xtol=1e-14, badvalue=None, name=None, longname=None, shapes=None, extradoc=None)\n", | |
| " | \n", | |
| " | cdf(self, x, *args, **kwds)\n", | |
| " | Cumulative distribution function of the given RV.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | x : array_like\n", | |
| " | quantiles\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | cdf : ndarray\n", | |
| " | Cumulative distribution function evaluated at `x`\n", | |
| " | \n", | |
| " | entropy(self, *args, **kwds)\n", | |
| " | Differential entropy of the RV.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information).\n", | |
| " | loc : array_like, optional\n", | |
| " | Location parameter (default=0).\n", | |
| " | scale : array_like, optional\n", | |
| " | Scale parameter (default=1).\n", | |
| " | \n", | |
| " | est_loc_scale(*args, **kwds)\n", | |
| " | `est_loc_scale` is deprecated!\n", | |
| " | \n", | |
| " | This function is deprecated, use self.fit_loc_scale(data) instead.\n", | |
| " | \n", | |
| " | expect(self, func=None, args=(), loc=0, scale=1, lb=None, ub=None, conditional=False, **kwds)\n", | |
| " | Calculate expected value of a function with respect to the distribution\n", | |
| " | \n", | |
| " | Location and scale only tested on a few examples.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | func : callable, optional\n", | |
| " | Function for which integral is calculated. Takes only one argument.\n", | |
| " | The default is the identity mapping f(x) = x.\n", | |
| " | args : tuple, optional\n", | |
| " | Argument (parameters) of the distribution.\n", | |
| " | lb, ub : scalar, optional\n", | |
| " | Lower and upper bound for integration. default is set to the support\n", | |
| " | of the distribution.\n", | |
| " | conditional : bool, optional\n", | |
| " | If True, the integral is corrected by the conditional probability\n", | |
| " | of the integration interval. The return value is the expectation\n", | |
| " | of the function, conditional on being in the given interval.\n", | |
| " | Default is False.\n", | |
| " | \n", | |
| " | Additional keyword arguments are passed to the integration routine.\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | expected value : float\n", | |
| " | \n", | |
| " | Notes\n", | |
| " | -----\n", | |
| " | This function has not been checked for it's behavior when the integral is\n", | |
| " | not finite. The integration behavior is inherited from integrate.quad.\n", | |
| " | \n", | |
| " | fit(self, data, *args, **kwds)\n", | |
| " | Return MLEs for shape, location, and scale parameters from data.\n", | |
| " | \n", | |
| " | MLE stands for Maximum Likelihood Estimate. Starting estimates for\n", | |
| " | the fit are given by input arguments; for any arguments not provided\n", | |
| " | with starting estimates, ``self._fitstart(data)`` is called to generate\n", | |
| " | such.\n", | |
| " | \n", | |
| " | One can hold some parameters fixed to specific values by passing in\n", | |
| " | keyword arguments ``f0``, ``f1``, ..., ``fn`` (for shape parameters)\n", | |
| " | and ``floc`` and ``fscale`` (for location and scale parameters,\n", | |
| " | respectively).\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | data : array_like\n", | |
| " | Data to use in calculating the MLEs.\n", | |
| " | args : floats, optional\n", | |
| " | Starting value(s) for any shape-characterizing arguments (those not\n", | |
| " | provided will be determined by a call to ``_fitstart(data)``).\n", | |
| " | No default value.\n", | |
| " | kwds : floats, optional\n", | |
| " | Starting values for the location and scale parameters; no default.\n", | |
| " | Special keyword arguments are recognized as holding certain\n", | |
| " | parameters fixed:\n", | |
| " | \n", | |
| " | f0...fn : hold respective shape parameters fixed.\n", | |
| " | \n", | |
| " | floc : hold location parameter fixed to specified value.\n", | |
| " | \n", | |
| " | fscale : hold scale parameter fixed to specified value.\n", | |
| " | \n", | |
| " | optimizer : The optimizer to use. The optimizer must take func,\n", | |
| " | and starting position as the first two arguments,\n", | |
| " | plus args (for extra arguments to pass to the\n", | |
| " | function to be optimized) and disp=0 to suppress\n", | |
| " | output as keyword arguments.\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | shape, loc, scale : tuple of floats\n", | |
| " | MLEs for any shape statistics, followed by those for location and\n", | |
| " | scale.\n", | |
| " | \n", | |
| " | fit_loc_scale(self, data, *args)\n", | |
| " | Estimate loc and scale parameters from data using 1st and 2nd moments.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | data : array_like\n", | |
| " | Data to fit.\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information).\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | Lhat : float\n", | |
| " | Estimated location parameter for the data.\n", | |
| " | Shat : float\n", | |
| " | Estimated scale parameter for the data.\n", | |
| " | \n", | |
| " | freeze(self, *args, **kwds)\n", | |
| " | Freeze the distribution for the given arguments.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution. Should include all\n", | |
| " | the non-optional arguments, may include ``loc`` and ``scale``.\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | rv_frozen : rv_frozen instance\n", | |
| " | The frozen distribution.\n", | |
| " | \n", | |
| " | isf(self, q, *args, **kwds)\n", | |
| " | Inverse survival function at q of the given RV.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | q : array_like\n", | |
| " | upper tail probability\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | x : ndarray or scalar\n", | |
| " | Quantile corresponding to the upper tail probability q.\n", | |
| " | \n", | |
| " | logcdf(self, x, *args, **kwds)\n", | |
| " | Log of the cumulative distribution function at x of the given RV.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | x : array_like\n", | |
| " | quantiles\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | logcdf : array_like\n", | |
| " | Log of the cumulative distribution function evaluated at x\n", | |
| " | \n", | |
| " | logpdf(self, x, *args, **kwds)\n", | |
| " | Log of the probability density function at x of the given RV.\n", | |
| " | \n", | |
| " | This uses a more numerically accurate calculation if available.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | x : array_like\n", | |
| " | quantiles\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | logpdf : array_like\n", | |
| " | Log of the probability density function evaluated at x\n", | |
| " | \n", | |
| " | logsf(self, x, *args, **kwds)\n", | |
| " | Log of the survival function of the given RV.\n", | |
| " | \n", | |
| " | Returns the log of the \"survival function,\" defined as (1 - `cdf`),\n", | |
| " | evaluated at `x`.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | x : array_like\n", | |
| " | quantiles\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | logsf : ndarray\n", | |
| " | Log of the survival function evaluated at `x`.\n", | |
| " | \n", | |
| " | moment(self, n, *args, **kwds)\n", | |
| " | n'th order non-central moment of distribution.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | n : int, n>=1\n", | |
| " | Order of moment.\n", | |
| " | arg1, arg2, arg3,... : float\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information).\n", | |
| " | kwds : keyword arguments, optional\n", | |
| " | These can include \"loc\" and \"scale\", as well as other keyword\n", | |
| " | arguments relevant for a given distribution.\n", | |
| " | \n", | |
| " | nnlf(self, theta, x)\n", | |
| " | \n", | |
| " | pdf(self, x, *args, **kwds)\n", | |
| " | Probability density function at x of the given RV.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | x : array_like\n", | |
| " | quantiles\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | pdf : ndarray\n", | |
| " | Probability density function evaluated at x\n", | |
| " | \n", | |
| " | ppf(self, q, *args, **kwds)\n", | |
| " | Percent point function (inverse of cdf) at q of the given RV.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | q : array_like\n", | |
| " | lower tail probability\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | x : array_like\n", | |
| " | quantile corresponding to the lower tail probability q.\n", | |
| " | \n", | |
| " | sf(self, x, *args, **kwds)\n", | |
| " | Survival function (1-cdf) at x of the given RV.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | x : array_like\n", | |
| " | quantiles\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | sf : array_like\n", | |
| " | Survival function evaluated at x\n", | |
| " | \n", | |
| " | stats(self, *args, **kwds)\n", | |
| " | Some statistics of the given RV\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | moments : str, optional\n", | |
| " | composed of letters ['mvsk'] defining which moments to compute:\n", | |
| " | 'm' = mean,\n", | |
| " | 'v' = variance,\n", | |
| " | 's' = (Fisher's) skew,\n", | |
| " | 'k' = (Fisher's) kurtosis.\n", | |
| " | (default='mv')\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | stats : sequence\n", | |
| " | of requested moments.\n", | |
| " | \n", | |
| " | ----------------------------------------------------------------------\n", | |
| " | Methods inherited from rv_generic:\n", | |
| " | \n", | |
| " | interval(self, alpha, *args, **kwds)\n", | |
| " | Confidence interval with equal areas around the median.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | alpha : array_like of float\n", | |
| " | Probability that an rv will be drawn from the returned range.\n", | |
| " | Each value should be in the range [0, 1].\n", | |
| " | arg1, arg2, ... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information).\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter, Default is 0.\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter, Default is 1.\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | a, b : ndarray of float\n", | |
| " | end-points of range that contain ``100 * alpha %`` of the rv's possible\n", | |
| " | values.\n", | |
| " | \n", | |
| " | mean(self, *args, **kwds)\n", | |
| " | Mean of the distribution\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | mean : float\n", | |
| " | the mean of the distribution\n", | |
| " | \n", | |
| " | median(self, *args, **kwds)\n", | |
| " | Median of the distribution.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | Location parameter, Default is 0.\n", | |
| " | scale : array_like, optional\n", | |
| " | Scale parameter, Default is 1.\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | median : float\n", | |
| " | The median of the distribution.\n", | |
| " | \n", | |
| " | See Also\n", | |
| " | --------\n", | |
| " | stats.distributions.rv_discrete.ppf\n", | |
| " | Inverse of the CDF\n", | |
| " | \n", | |
| " | rvs(self, *args, **kwds)\n", | |
| " | Random variates of given type.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information).\n", | |
| " | loc : array_like, optional\n", | |
| " | Location parameter (default=0).\n", | |
| " | scale : array_like, optional\n", | |
| " | Scale parameter (default=1).\n", | |
| " | size : int or tuple of ints, optional\n", | |
| " | Defining number of random variates (default=1).\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | rvs : ndarray or scalar\n", | |
| " | Random variates of given `size`.\n", | |
| " | \n", | |
| " | std(self, *args, **kwds)\n", | |
| " | Standard deviation of the distribution.\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | std : float\n", | |
| " | standard deviation of the distribution\n", | |
| " | \n", | |
| " | var(self, *args, **kwds)\n", | |
| " | Variance of the distribution\n", | |
| " | \n", | |
| " | Parameters\n", | |
| " | ----------\n", | |
| " | arg1, arg2, arg3,... : array_like\n", | |
| " | The shape parameter(s) for the distribution (see docstring of the\n", | |
| " | instance object for more information)\n", | |
| " | loc : array_like, optional\n", | |
| " | location parameter (default=0)\n", | |
| " | scale : array_like, optional\n", | |
| " | scale parameter (default=1)\n", | |
| " | \n", | |
| " | Returns\n", | |
| " | -------\n", | |
| " | var : float\n", | |
| " | the variance of the distribution\n", | |
| " | \n", | |
| " | ----------------------------------------------------------------------\n", | |
| " | Data descriptors inherited from rv_generic:\n", | |
| " | \n", | |
| " | __dict__\n", | |
| " | dictionary for instance variables (if defined)\n", | |
| " | \n", | |
| " | __weakref__\n", | |
| " | list of weak references to the object (if defined)\n", | |
| "\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 27 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "so we define our own chi^2 function" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "def chi2_lin(x,y,unc_y,f):\n", | |
| " chi = (y-f(x))/unc_y\n", | |
| " chi2 = chi**2\n", | |
| " return(chi2.sum())" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 28 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "and then evaluate it for the various fits" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "chi2_lin(x,y,random_gauss_2,y1)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "metadata": {}, | |
| "output_type": "pyout", | |
| "prompt_number": 58, | |
| "text": [ | |
| "1146.9897355522685" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 58 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "chi2_lin(x,y,random_gauss_2,y2)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "metadata": {}, | |
| "output_type": "pyout", | |
| "prompt_number": 59, | |
| "text": [ | |
| "1262.3276543106767" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 59 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "def y3(x):\n", | |
| " return(x)\n" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 56 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "chi2_lin(x,y,random_gauss_2,y3)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "metadata": {}, | |
| "output_type": "pyout", | |
| "prompt_number": 60, | |
| "text": [ | |
| "1159.1015927721865" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 60 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [] | |
| } | |
| ], | |
| "metadata": {} | |
| } | |
| ] | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment