LinearFit

{
 "metadata": {
  "name": "",
  "signature": "sha256:91bec98ee2b02f23a903259ca2069b693a2ea62a054a55e43a0cb48959b18e0a"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from numpy.random import rand,choice\n",
      "import matplotlib.pyplot as plt\n",
      "from astropy.io import ascii"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dat=ascii.read('data_err.dat')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 5,
     "metadata": {},
     "source": [
      "First assume the same error for all the data points."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Follow Wall and Jenkins (p. 129) - assume $\\sigma_i = \\sigma$ is the same for all of the points.\n",
      "\n",
      "Define $Sx=\\sum_i(\\rm X_i), Sy = \\sum_i(\\rm Y_i) , \\rm Sxy = \\sum_i(\\rm X_i*Y_i), \\rm Sxx = \\sum_i(\\rm X_i^2)$ then\n",
      "\n",
      "$$a = {{N\\, \\rm Sxy - \\rm Sx\\, \\rm Sy} \\over {N \\, \\rm Sxx -(\\rm Sx)^2}}$$\n",
      "\n",
      "$$b = {(Sy - a\\,Sx) \\over N } $$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.array(dat['xi'])\n",
      "y = np.array(dat['yi'])\n",
      "yerr = np.array(dat['err'])\n",
      "Sx = np.sum(x)\n",
      "Sy = np.sum(y)\n",
      "Sxy = np.sum(x*y)\n",
      "Sxx = np.sum(x**2)\n",
      "#print Sx, Sy, Sxy, Sxx "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = len(x)\n",
      "#print n\n",
      "a = (n*Sxy - Sx*Sy) /(n*Sxx - (Sx)**2)\n",
      "b = (Sy - a*Sx)/n\n",
      "print a, b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.942153821198 -0.128180870694\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plot the data, solution to the linear equation and the linear model with the assumed values for the simulated data points."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.errorbar(x, y, yerr, fmt='o')\n",
      "plt.plot(x, b + a*x, 'r')   # solution\n",
      "plt.plot(x, -0.1+0.8*x, 'g') # actual model\n",
      "plt.xlabel(\"X\")\n",
      "plt.ylabel('Y')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "<matplotlib.text.Text at 0x1065229d0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEPCAYAAABV6CMBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPWh9/HPIYGwqCzKKlHbqNe6VXG3lgwqCRpbucWr\nUrSt9VZ9LAG1ULZAEnGp+FgEaq1Vr9uV2oq2tURJQJnkeRBFUWlFqBgeK8sFFEHcAiQ5zx9ntjOZ\nmZxZz8zk+3698mLOmTNnfgfIfOe3HhAREREREREREREREREREREREREREZEcNQbYCGwCpsU47iyg\nFfhBJgolIiLZqQD4ADgG6A68A3wrynGvAEuBcZkqnIiIRNbNxfc+Gys4PgQOAs8Al0c4rhJYAnyc\nsZKJiEhUbgbHkcCWkO2tvn3hx1wOPOjbNjNQLhERicHN4HASAvcD033HGr4fERFxUaGL770NKA7Z\nLsaqdYQ6A6sJC+AI4BKsZq0XQg8qKSkxm5ub01RMEZG81Qwc63Yh4lGIVehjgB5E7xz3e4zoo6rM\nfFZdXe12EdImn6/NNHV9uS7fr48Em//drHG0AhOBeqyRU48CG4Abfc8/5FK5REQkBjeDA+Al30+o\naIFxXZrLIiIiDrjZOS4OeTwet4uQNvl8baDry3X5fn2JypdRSr7mOhERccowDEggB1TjEBGRuCg4\nREQkLgoOERGJi4JDRETiouAQEZG4KDhERCQuCg4REYmLgkNEROKi4BARkbi4vVaViEjW8nqtH/9j\n/wokHk/wcVekJUdERBwwDMi3jxktOSIiIhmh4BARkbgoOEREJC4KDhGRZOzaZXWA/OIXbpckYxQc\nIiKJqqyEwYOtxzfd5G5ZMkjDcUVEYqira2LhwgagkPLyViZNKqNi2KEwYoR1wNSpMG+eq2XMNA3H\nFRGJoq6uicmT62luvhOAAlr5e4/hnHhgp3XA7t0wYICLJUxOosNxVeMQkbyR6gl7Cxc2BELjSv7I\nH7kaDsCvTh3L9HV/Trq8uUrBISJ5IzQgDCMYIonav7+Qb9JMM8cC8DancRZvcEH/O5ie3KlzmjrH\nRUSi8DbWBkLjLNYwgrdpo5CePdtsxz27/lnWbFvjRhFdoeAQEQm3ZIlVZfExMHmTswAoKZlJZeVo\nAFZsXoFRa3DlkitZsXmFK0V1g5qqRET8Wluhe/fg9rp11G3ZS/mi2dTXF1Be3kZl5RiKz+yHUWsF\nS4FRwJZbtzD00KEuFTrzNKpKRPKKf/hsQ0MhZWW+4bMVIzt/4c9+Bo88Yj0+/3xYtcr2tGHAls+2\nUjy/OLBv/c3rOXHgiaksfkZpVJWIdHnhw2cbGqC5eRZA9PDYtg2GDw9uf/45HHKI7ZDPWj6DW0+h\neP4WALw/9lJ6TGnqLyBHqI9DRPJG6PBZv+bmO1m0aHnkF/TqFQyNe++11k0PCY0vD3yJUWvQ755+\n0HcLfxj3B8xqs0uHBrgfHGOAjcAmYFqE5ycA64C/A6uAUzNXNBHJNfv3R25EaWkpsO9Ytsxqe2pp\nsbbb22HKlMDT7WY7Rq3BIXdbITKozyCoMbn65KvTUu5c42ZTVQHwG+BiYBvwBvACsCHkmM3ASOAz\nrJD5PXBuZospIrmiqKg14v7A8Nn2digICZHVq+Fc+0dKrzt70dLaEthun9OOYRgYU1Ne3JzlZo3j\nbOAD4EPgIPAMcHnYMauxQgPgdWA4IiJRTJpURknJLNu+wPDZadOCoXHccVazVEholD1VhlFrBELj\nq5lfYVab/g5kCeFmjeNIYEvI9lbgnBjHXw+8mNYSiUhO83eALwoZPnvbj86l7LKQPomw9aVmvjyT\nu//v3YHt7bdtDwytDV3CpLQUamqsx/4Z6l31nuRuRuk4rOann/m2r8EKjsoIx44CHgC+A+yJ8LyG\n44pkuUx/yBoGmMVHwRbf99OqKpg7N/D84n8sZsLzEwLba29Yy4ihI5J7vxz7GMrF4bjbgOKQ7WKs\nWke4U4GHsUImUmgAUOP/KgB4PB48+Rz3Ijko1etIxfTss5hcGWzTaG8PzAR/betrnPfoeYFDl/zH\nEsadOC6NhckeXq8Xbwr+4t2scRQC/wQuArYDa4Dx2DvHjwJewaqNvBbjXKpxiOSQtH07N03oFtJ1\n+/LLcOGFAGz5bAtH3X9U4KnbPbczu3R20m+Z8ITDLJCLNY5WYCJQjzXC6lGs0LjR9/xDwBygP/Cg\nb99BrE51ERG7Sy6xhtn6GJiYF1pzMfzDagHGnjCWP1+VmiXRE5pwmAfyZbiAahwiOSAt384//hgG\nDQpuf/ABlJRgdGuH6uDQ24G9B7Jr6q7k3itMeXkVDQ13RNg/m2XL5kZ4RXbJxRqHiHQhafl2HjpU\n9pvfxPtoM96noNYwoDr41Csj2xk1KvXfkx1POMwzbs8cF5EuIu7lQGJ58UV7aBw8CM3NjGo0rNDw\n8c/FSEdogIMJh3lKwSEiGZGyb+eGARUV1uN588A0Oefx7wSWOQd4f+L7mNUmvbr3SrS4jsSccJjH\n1FQlIhmR9Lfz666Dxx8Pbpsm01dM556QwFg6fikVx1ckUcr4RJpwWFk5Jq87xkGd4yKSIeF9HGB9\nO1+woJMP2n37oG/f4Pa6dTxXuIkrnr0isGvuqLlUjaxKR7Ed60oTABUcIpIxdXVNLFq0POTb+ejY\noRHaj9GzJ+/+6w1OefCUwK6Lv3kxy69NoI8kDRQcuUfBIZJDOv2QXbUKLrggsLlnz3YGLBhmO8Ss\nzq7feQVH7lFwiOSQmB+yIbWMtmlTKex1r+3pbAqMXF/kUMGh4BDJap1+yP7yl9Zd+HyMGvvr/ffF\nkNRRcCg4RHJTS4t1C1ef8MDYN30fhxYdmtkydRGaOS4iuadHD2vyHnDOf8KakFu1vT/xfY47/DiX\nCiaxaAKgiGTeunVWX8bBg0y/2Kpl+ENj6filmNWmQiOLqcYhIpnl66d47ltwxVXB3dkwF0OcUXCI\nSEReL/z+9028/noDu3cXcvjhrZxzThk33DAysRFDP/0pPPYY6wbDaf8ruDub5mKIMwoOEYnoyy+b\nWLOmns2brZnen30GhjGLCRMA4lhS4+BB6NGDnX1gSI39qWwaWivOaVSViESUkntNGAat3aD7HPtu\nBUZ20KgqEUmppFazfecdOP30DkNr2+a00c3QmJxcp+AQkYgSXs3WMKzAqAnu+vSXn9K/V/9UFU1c\npugXkYjivtfErFkYtYatlrH2hrWY1aZCI8+oj0NEonK0mm17OxXXFvDi8cFdj3zvEa4fcX1mCytx\n05IjCg6RlKmra6KmpoGPPy6koKCVXr3KuOIKKzBCF/BbcK7BLZcEXzf+5PEsHrc44+WVxKhzXERS\n4u67m7j33nr27AnecKl//1kUFcGMGVZ4vPras3yn/krwhUaPbt3ZP/uAG8UVFyg4RMTG622whQbA\nnj130tg4mxtvOZnD5x1ue05Da7seBYeI2EQchmu0UX/eHRw+Lzivw6xqhQIHQ3Ml7yg4RMSmwzDc\nGnsT+P9++hIavjGCumWrYt/2VfKWgkNEbDyeMt54YxZ7Jt9l27/vLjjsgMkUgE3Q3GwN1VV4dD0a\nVSUiNkat/WNh7UNwe8EV/HXrsx2OPeus2axZ43D5Eck6GlUlkgfq6ppYuLCB/fsLKSpqZdKksox9\no69YXMGLm14MbD/yV7j+bcA02eupga0dX9O7t/o4uiIFh0iWqKtrYvLkepqbgyOaMtEctOC1BdxS\nf0tge/w/YPFzwP791h36SGL5EclLbi85MgbYCGwCpkU5ZqHv+XXA6Rkql0jG1dQ02EIDoLn5Tmpr\n03Ovile3vIpRawRCo0crmDWwuOcEMM1AaIDV79G/v335kf79Z1JaGmX5EclrbtY4CoDfABcD24A3\ngBeADSHHXAocCxwHnAM8CJyb2WKKZMbBg5F/HQ8cSG1z0M4vdjLkviG2fWaN/0HkvsIZM0Zy6qmw\naNFsWloK6NmzjcrKMeoY76LcDI6zgQ+AD33bzwCXYw+O7wNP+B6/DvQDBgM7M1NEkcwZPDhyc9CQ\nIfbmoET7QQ62HaTHHT1s+wKBsWIFXHRRzNdXVIxUUAjgbnAcCWwJ2d6KVavo7JjhKDgkD02aVEZz\n8yxbc5W1Gu2YwHai/SDhI6XaaqGbv3KhEYkSJzeDw+n/1vChYhFfV1NTE3js8XjwJHRTZBH3+D/4\nYzUHLVwYuR9k0aLZEYMjPDB2zYOBX/k2Pv8cDjkktRchWc3r9eL1epM+j5vzOM7FutWL/+vUDKAd\nuCfkmN8BXqxmLLA60kvpWOPQPA7pEo49tobm5poO+0tKavjgg+D+8MBofAxG/su34fHAypXpKqLk\nkETncbg5qupNrE7vY4AewFVYneOhXgB+5Ht8LrAXNVNJF1ZSErkf5NhjrX6Q4vnFttC4450BmDUh\noWGaCg1JmpvB0QpMBOqB94A/YnWM3+j7AXgR2IzVif4QcHPmiymSPaLdla/32E0YtQZb91mz9E4f\n9G3MGpj1l0+tg377W/VlSMpoyRGRHOO/K19LSwGfD/knb33rGdvzgZFSgR363ZDIdAdA/XJIF7J1\n31aK5xfb9pnnLoMxwRFYbN8OQ4dmuGSSSxQcCg7pAiLOxag2wQj7VdbvgzigRQ5F8lyHuRhz2uhW\nVm4PDQWGZICCQyTLhQfGjl/sYHDRAOgWshTJnDlQW5vhkklXpeAQyVLhgbH82uVc/M2L1SwlrnN7\ndVwRCVN0R5EtNKZ9ZxpmtcnF/9PLHhrr1ys0xBWqcYhkiSkNU7hv9X2B7eGHDWfLrb6l2lTLkCyi\n4BBx2UubXuLSxZfa9pnVvmAYNw6efz74RHt7xxARyTAFh4hLIs7F8AWG95V2PBcFO7/r+nyLp7//\nO25oNND6neI2BYdIhkWdi+FnGHhCnjMw4UsoWTOLCRMAdE8McVe+1Hk1AVByQsS5GIZvjMo//wkn\nnBB47nxWsZrzbceXl89m2bK5aS+ndA2aACiSxTrcF2PKLgb2GRhygP15T2k1qxvtoQHQ0pLYbWS9\nXuvH/9jf3OXxoKYviZuG44qkkVFr2EKj8SeNmNVmMDTmzLGHxsGDYJoUFUVePr1nz7aI+zvj8cBZ\nZzWxenUVjY01rF5dxVlnNSk0JCGqcYikQfH84sAS5wBzR82lamSV/aDQwCgrg/r6wKaT28jGI/yW\nsw0Nzm45KxKJgkMkhcLnYpw25DTevvFt+0EO5mT06TOSs88GmM0nnxRwxBFtnH32GPr0SexDPt5b\nzorEouAQSYGYczH8tm+HI48MbtfXWzWNCKy+h5GkagTV/v2Rf9UT7TORrk3BIZKEWHMxbFye+Z3q\nPhPp2hQckvfSMaKo07kYfr/9Lfz858HtL7+E3r0Te9MkpLrPRLo2zeOQLsUwkv+yH3MuRvib+f3b\nv8HGjcm9cZL8t5ytry+gvLyNysrR6t/o4nQHQAWHOJBMcIQHxs4pOxnUZ1DHA/v2hX37gttZ8H9T\n8zgkEgVHFvxySvaqq2ti4cIGGhoKKStrZdKkMsfftsMDo/EnjYw8OsJr9+6F/v2D208/DT/8YTLF\nFkkrzRwXiSLROQyO5mL4adlz6UI0c1zyXvQ5DMsjHj+lYQpGrREIjdOGnIZZbUYOjSVL7KGxe7dC\nQ/KeahyS95zOYXA0FyNUaGD06AH79ydcRpFcouCQnNdZx29ncxgcz8XwO/10eOedkINVw5CuRZ3j\nklcijZoK7+MAaw7DffMvZuxbF9mOjRkYLS3Qq1dw+9e/hltvTUWxRVyhUVUKDiH6cNvwOQz1591h\nez7qXIzQE4fS/zfJAwoO/SILnc/TcDwXw2/lSrjwwuD2Rx9BcXH040VySC4Oxx0A/BE4GvgQuBLY\nG3ZMMfAkMAgwgd8DCzNXRMkV/nkaUEh5uX2ehtcLoxrtvxs/MRv5sWckg/rEOKlqGSIRuVnjmAd8\n4vtzGtAfmB52zBDfzzvAIcBaYCywIew41Ti6sMh9GLNYsKCcy94stR1785k380DFA7FPOHYs/PWv\nwe329o4hIpIH0tFU9RJwM/D/EixTZzYCpcBOrHDwAifEegHwF2AR8HLYfgVHHkh0WYzy8ioaGux9\nFowbD6c8E9gsKiiipaoldgHa2qAwpBL+y1/CPfc4KrtILkpHU9V/AfXAE1i1goMJlSy6wVihge/P\nwZ0cfwxwOvB6isshWSI0IAwjGCKdsc3TOPkZuGK87fmYI6X81Cwl4lis4HgWq9YxB3gTeAqrnwHf\nn792cP7lWLWJcLPCts2Qc0dyCLAEmAx8EemAmpqawGOPx4NHK7d1GUVFrdB/M0wuse0vX13FsmVz\nY7/4nXeseRl+774LJ52UhlKKuM/r9eJ1+o0shs6qKEVY/Q8TgGeA9pDnapN8742AB9gBDAVWErmp\nqjuwFCvE7o9yLjVV5ahozVO1tc6+9Ee6LwY1JiUlM1mwYEzshQxVy5AuLh19HGOwahV/wwqJrxIq\nWXTzgN3APVid4v3o2DluYDWV7QZizbRScOSB0KG0TpY/Dx9aW7Z6Fg31hZ3fa2LiRHggpIO8rQ26\nadk26XrSERz/B7gJWJ9gmTozAPgTcBT24bjDgIeBCuACoAn4O8GmrBnAsrBzKTjygNPgCA+Mbbdt\nY9ihwzp9HaZpD4gf/tBa+lyki0pHcBjE7nfIJgqOPGAYsHRp9PtmhAfG0vFLqTi+osM5Iv5XULOU\nSAeaOa4PgpxnGE2UlHScj9F87V2242LNxegQHBs2wIknBrfr66GsLJXFFslZCg4FR84zjCogZD5G\nAnMxbMGhWoZITLm45IhIGN9/x1Ofgh/8yPZMrLkYoSOzSkvh9bMrOeeN3wQPaGmBoqLUFlWkC1Nw\niOsC60wd/hZUhn35qTEpL58d8/W2meWhtYyBA2HXrlQWVURQcIjL6uqamHTLi2y+5h44P+SJGquG\nUVIyk8rKMZ2fSM1SIhmjPg5xVfhIKWpfAfMV4CPKy4+KPR8DYNs2GD48uP3UU3DNNWkpq0i+Uee4\ngiOndAiM+7bB58NCdtRgmjWdnES1DJFkJBocmi4rGWXUGrbQOH3DlVazlC00ANqin2TePHto7N2r\n0BDJIAWHZER4YNx85s2Y1SZzr/05JSX2NS9LSmYCo6OcyIBp04Lbpgl9+6ahxCISjZqqJK3GPzee\nZ94NzsXoUdCD/VX7bceE3w+8snI0l1020l6JULOUSMqpj0MfJFnlpU0vceniS237OrsvRsS1qvbu\nhf79gwfdey9MmZLi0op0TZoAKFlh676tFM8vtu1zdCMlH/9tVUpLUS1DJEupxiEpEem+GPEEBoTU\nMh5/HK67LvjEtm0wLLzzXESSpRqHuCZ8aG3bnDa6GQmOu1AtQyTrKTgkYeGBsXPKTgb1GZTYyXr3\nxuTr4LYCQyRraTiuxC18aG3jTxoxq82EQqOpocWqZXxthca7J11FTbVJCm6LLCJpoj4Ocax4fjFb\n920NbM8dNZeqkVWJn1DNUiKu0sxxSZupDVMxao1AaJw25DTMajPx0Kivt4fGhg0KDZEcoj4OiSqR\nuRidUi1DJOcpOKSDZOdiRHTeefDaa8Ht9vaOISIiOUHBIQGpmIvRQVsbFIb8N7v8cvjLX5I7p4i4\nSsEhQIrnYgROqmYpkXyk4OjinM7FCL2vt9cbvFWr7batfuvXw8knB7dffdVqqhKRvJAvjcwajhun\n8MBo/EkjI4+Ocae90NcaMSoPqmWI5AwNxxVHLn36UltozB01F7PadBwaUU2dag+N1laFhkieUlNV\nF3H/a/dza/2tge3xJ49n8bjFyZ/YNKFbyPePf/93eP755M8rIllLwZHnVn20igseuyCwHelGSglT\ns5RIl6TgSLO4OpVTaOcXOxly3xDbvqSH1voU8xEYRwd3rFyZ3osRkaziVuf4AOCPwNHAh8CVwN4o\nxxYAbwJbge9FOSYnOsdjdiqnSFt7G4Vz7d8HUhUYgGoZInkk1zrHpwPLgeOBl33b0UwG3gP0CdUJ\no9awhUbbnLbUhcb8+fbQ+PprhYZIF+VWjWMjUArsBIYAXuCECMcNBx4H7gRuQzWOiEY/NZoVm1cE\ntvdN38ehRYem7g1CA2PECIy31iozRPJArt0BcDBWaOD7c3CU4+YDU4HDMlGoXDNjxQx+tepXge3t\nt21n6KFDU/cG3btbw2r9/GlhBO8NHqnfxq1+HRHJjHTWOJZj1SbCzQKeAPqH7PsUq98j1GXAJcDP\nAQ/wC2LUOKqrqwMbHo8HTxZ+QqWqxvH035/mmj9fE9h+aMRatq8dAXT+Qe3oQ333bjjiiOCLliyB\nceMiXkdn15SJfh0Rccbr9eINuUtabW0tJJADbjZVeYAdwFBgJR2bqu4CrgVagZ5YtY7ngB9FOF+X\naKp6betrnPdocOmOJf+xhHEnjrMdE897RDzWQee3gkMkPyTaVOVWcMwDdgP3YHWM9yN2B3kpMIUu\n2sfx0WcfcfT9weGvt3tuZ3bp7Kjv4a98dVb7sJXn6afhmmAthr17oW/fqO+h4BDJfbnWx/Er4E/A\n9QSH4wIMAx4GKiK8pst9/Hxx4AsOvTvYyT32hLH8+ao/d/o6f/+DYeDs3t2htYx+/WDPnoiH1dU1\nsXBhA1BIeXkrkyaVAUkuVSIiOUeLHDqUig5fp9++2812Cm4vCGwP7D2QXVN3xf0enb3f+8bxHM+m\n4I4YB9fVNTF5cj3NzXcG9pWUzKK5uRzTjB4eqnGIZK9ca6pKtYw2VSX6Yejkdb3u7EVLa0tgu31O\nu/8fN+73iPp+X38NvXsHtx98EG66KeZ5y8uraGi4I8IzszHNuY7KIyLZJdeaqiRM+FyMr2Z+Ra/u\nvVL/RgnO/N6/P9p/lYIo+0UkX2lZdZfNWDEDo9YIhMb227ZjVpupD43Vq+2h8fHHGHF0GxUVtUZ5\npi25colIzlGNwyXhczHW3rCWEUNHpOfNQgOjpAQ++CDuU0yaVEZz86ywPo6ZNDePSUUJRSSHKDgy\nIHQ00jnj/sWaUx8PPBdpLkaqLGIiGA8EdyTR2VBRYXWAL1o0m/r6AsrL26isHMNll2lUlUhXo87x\nBMTT4RsYjfTJjXBrcC7GhGE/5b9/9mh6ytbaBoUh3wl+9zu48cbIxybZya95HCK5K9dWx+0y7l/0\nEs3fXRUMjQ1jocbkk+eGpeX93uUke2iYZsTQEBFJlJqq0sQ0Ta5/4XpWnPeYtaNhHrw6NfB8S0uK\nRyN9/DEMGsRJ/u0dOzCGDO56syZFJO0UHGlwR9MdzF5pLQlS/D9nsOWhNwivDfbsmcLRSEOHwo4d\nALzMhVxkvpy6c4uIhFFTVQo9ue5JjFqD2StnU15SzoGqAzz4vV9TUlJlO66kZCaVlaOTf8PGRqsT\nwRcatLdzMZkJjbq6JsrLq4AaysurqKtrysj7ioj7VONIgRWbVzD6KSsIjhtwHGtvWBu4kVK00Uj+\n/QkxTegWkvleL3VfGCwcMxsn60hFWnMqvvLYlx9paIDm5llA8HqTfw8RkfQyM8n/dut2rDOpwaQG\ns6C2wNy+b7uj1yWlttY6EZjm4MGmaZrm0qWNZknJzMBuMH3bjR1eHu3YpUs7Hhv9OmbZXu//KS+v\nStl7iEj60QUXjw2V2b/sw7YEAoMazPW71jt7XTLF3LPH/im9a1fgqbKyyB/kUNXhNNGO9X/oO7uO\n6ojnKC2tTtl7iEj6kWBwqI8jThNfnAi3FQPg/bEXs9rkxIEnpvdNTz0V+vtumHjLLdbn8MCBgaej\nryO1pUP/Q7Rj4xvlFXn5EX+Hf2reQ0Sylfo44vTdo77LA9MuwPzH1el/s7Vr4cwzg9ttbfa+DZ/o\n60gV09Aw19b/EO3Y+EZ5lfmWVLcvP1JZOSZmeVI6kkxEJEkZrt5l4HWhbTxLl8Y8NFKfAswwodFh\n/8OMOPs4rPOUl1eZUG2Wl1fZXp+K9xCR9CPBpiotORIH/0ihhoZCysriHynkaPmNhQth8uTgC9rb\nHZdt0aLl1NdvAYqB0YSOqiotrcHrrQk71j/Ka3TC1xHtmpJ9DxFJP93IKc3BEe0OeAsWlDv+QIwZ\nHF9+CYccEtzesgWGD4+7nIZRBXS84VJ5+WyWLZsbdqzzdaSi3QGxtlZrVYnkqkSDI1+kvUqXmtFI\nUZ4YNSp4wh//OKlygvNmolT8tXV2jgy3IopIHEiwqUqd4w6lZaTQe+/BSScFtw8cgO7dEz8fACNZ\nsCDFEw5FREIoOBxK+Uih0Jsr/eEPcHXqRmlVVIykomIkhgHLlqXstCIigOZxODZpkjUENVRCa049\n/rg9NEwzpaEhIpJu+dIp4muuS6+kRgodOABFRcHt99+H445LeRkzfZMl3chJJHdpVFWW3gEQgD/9\nCa66ynr8ve/BCy+kpVyg4BAR5xQc2Rgce/bAgAEAfNFnMD8Y8SHnX9gTsIaz+oe0pqtsCg4RiUXB\nkW3BMXMm3H239fjNN+GMM9JeLlBwiIhziQaHRlWl2oYNcKJv0cObboIHH3S3PCIiKabgSJX2dhg1\nCpp8K9Hu2AGDB7tbJhGRNHBrOO4AYDnwPtAA9ItyXD9gCbABeA84NyOli9ff/gYFBVZoPPKI1Taj\n0BCRPOVWH8c84BPfn9OA/sD0CMc9ATQC/4VVO+oDfBbhOHf6OD7/3Or8bm2Fb3wDNm6EHj0yVo6Y\nZUN9HCISW6J9HG7VOL6PFQr4/hwb4Zi+wHexQgOsuwdFCg133HUXHHaYFRqrVsHmza6HhohIJrhV\n49iDVcvwl+HTkG2/04CHsJqovg2sBSYDX0U4X+ZqHM3NcOyx1uNrroEnn7TPBHdZJmoc0VbK9Q8x\n7ux5EckO2TgcdzkwJML+WVi1jNCg+BSr3yPUmcBq4HzgDeB+YB8wJ8I5MxccZ5wBb70FH30ExcWZ\nec9OxLPkuT7URcQvG4Mjlo2AB9gBDAVWAieEHTMEKzi+4du+AKsf5LII5zOrq6sDGx6PB08aPwVz\npd0+V8opIpnh9Xrx+r85ArW1tZBDwTEP2A3cgxUG/YjcOd4E/CfW6KsaoBdWZ3q47JsAmAVypZwi\n4o5cq3EGVCwHAAAGBklEQVQMAP4EHAV8CFwJ7AWGAQ8DFb7jvg08AvQAmoHryKZRVVlITVEi4lSu\nBUeqKThEROKUa8NxRUQkRyk4REQkLmqqckh9ByKSb9THoU4HEZG4qI9DREQyQsEhIiJxUXCIiEhc\nFBwiIhIXBYeIiMRFwSEiInFRcIiISFwUHCIiEhcFh4iIxEXBISIicVFwiIhIXBQcIiISFwWHiIjE\nRcEhIiJxUXCIiEhcFBwiIhIXBYeIiMRFwSEiInFRcIiISFwUHCIiEhcFh4iIxEXBISIicVFwiIhI\nXBQcIiISF7eCYwCwHHgfaAD6RTluBrAe+AewGCjKSOlERCQqt4JjOlZwHA+87NsOdwzwM2AEcApQ\nAFydofJlFa/X63YR0iafrw10fbku368vUW4Fx/eBJ3yPnwDGRjhmH3AQ6A0U+v7clpHSZZl8/s+b\nz9cGur5cl+/Xlyi3gmMwsNP3eKdvO9ynwH3AR8B2YC+wIiOlExGRqArTeO7lwJAI+2eFbZu+n3Al\nwC1YTVafAc8CE4CnU1dEERGJl+HS+24EPMAOYCiwEjgh7JirgNHAf/q2rwXOBX4e4XwfYAWNiIg4\n1wwc63YhnJoHTPM9ng78KsIx3wbeBXphBdwTRA4NERHpAgZg9VeED8cdBtSFHPdLgsNxnwC6Z7CM\nIiIiIiLSVTmdQNgPWAJsAN7D6iPJdk6vDay5LW8Df8tAuVLFyfUVY/V7rcdqrpyUsdIlbgxW390m\ngs2w4Rb6nl8HnJ6hcqVKZ9c3Aeu6/g6sAk7NXNGS5uTfDuAsoBX4QSYKlUJOrs+D9VnyLuDNSKlc\nMA+rGQusv4hIfSRgNW/91Pe4EOib5nKlgtNrA7gNa5TZC+kuVAo5ub4hwGm+x4cA/wS+lf6iJawA\na4DGMVjNqe/QsbyXAi/6Hp8DvJapwqWAk+s7j+Dv1xhy5/qcXJv/uFeApcC4TBUuBZxcXz+sL2nD\nfdtHZKpwmbaR4NyPIb7tcH2BzRkrUeo4uTaw/pFXAKPIrRqH0+sL9RfgorSVKHnnActCtqfTcTWE\n32GNFPQL/XvIdk6uL1R/YGtaS5Q6Tq/tFuBm4DFyKzicXN/NwO3xnDRXFzl0MoHwG8DHWP/QbwEP\nY80+z3ZOrg1gPjAVaM9EoVLI6fX5HYPVrPN6GsuUrCOBLSHbW337OjtmOLnByfWFup5g7SrbOf23\nuxx40Lcdad5ZtnJyfcdhNSGvBN7EmvoQUzonACYr2QmEhVjrXE0E3gDux0raOSksY6KSvbbLgF1Y\nbZKelJYsNZK9Pr9DsPqoJgNfpKZoaeH0gyR83lSufADFU85RWM3D30lTWVLNybX5PztMrH9Dt+a/\nJcLJ9XXH+qy8COvL9WqspsZN0V6QzcExOsZzO7E+mPwTCHdFOGar7+cN3/YSYlevMynZazsfa72v\nS4GewGHAk8CPUlvMhCV7fWD9Z34O+G+spqpstg2rQ9+vmI5NNeHHDCd31l5zcn1gdYg/jNXHsScD\n5UoFJ9d2BvCM7/ERwCVY6+jlQt+ik+vbAnwCfO37acKaRxc1OHKVkwmEYP0FHO97XAPck95ipYTT\na/MrJbf6OJxcn4EVhPMzVagkFWLNwD0G6EHnnePnkjudx+Ds+o7C6oTNhZGLoZxcW6jHyK1RVU6u\n7wSs/tICrBrHP4ATM1fEzHE6gfDbWDWOdcDz5MaoKqfX5ldKbnzz8XNyfRdg9d28g9Uc9zbWt9hs\ndgnW6K8PsO4jA3Cj78fvN77n12E1DeSSzq7vEWA3wX+vNZkuYBKc/Nv55VpwgLPrm0JwsnUuDH8X\nEREREREREREREREREREREREREREREfErxlpgs79vu79v+yjXSiSSQgVuF0AkD+3DmqV7JdYy3Auw\nJj3mysJ/IiLigkKsGeK3YM3G1Zc0ERHpVDnW0inZfC8Rkbjl6v04RHLBJcB24BS3CyIiItnvNKz7\nNxcD/yLy/UlEREQAa1n41QSbqCZi3VdEREQkohuAP4RsdwPWAt91pzgiIiIiIiIiIiIiIiIiIiIi\nIiIiIiIiIiIiIiIikvf+P4G8/qL8DntgAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1063b3910>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "err = np.array(dat['err'])\n",
      "w = 1./err**2\n",
      "s = np.sum(w)\n",
      "sx= np.sum(w*x)\n",
      "sy =np.sum(w*y)\n",
      "sxx = np.sum(w*x**2)\n",
      "sxy = np.sum(w*x*y)\n",
      "delta= s*sxx - (sx)**2\n",
      "b = (sxx*sy - sx*sxy)/delta\n",
      "a = (s*sxy - sx*sy)/delta\n",
      "print a, b\n",
      "plt.errorbar(x, y, err, fmt='o')\n",
      "plt.plot(x, b + a*x)\n",
      "#a = 0.8\n",
      "#b = -0.1\n",
      "#plt.plot(x, b + a*x)\n",
      "siga = np.sqrt(s/delta)\n",
      "sigb = np.sqrt(sxx/delta)\n",
      "print siga, sigb\n",
      "plt.plot(x, 1.1*x +0.1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.508490528215 -0.0439528852037\n",
        "0.0312698028105"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0.00268271474156\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "[<matplotlib.lines.Line2D at 0x10671b7d0>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcFNXd7/FPM2yCkUWIgBs6akATlygIapg26jSCEW9i\nXC5GEnzUqwgkbiiLDHFF4wK4L48Sn+ujEQ0xDMtAQjvqRYGIOyiM8ARQVJRFwRmYoe8fp5veqme6\nu6qrevm+X69+TZ2umqpTDPPrM7+zFIiIiIiIiIiIiIiIiIiIiIiIiIiISBEbDKwCVgPjUhzjB1YA\nHwBBV2olIiI5UwasAXoDbYB3gL4Jx3QGPgQOCpe7uVU5ERGx1srm9/fHBP91wG7geWBYwjH/G3gJ\n2BAub7Z5TRERsclu8D8QWB9T3hB+L9aRQFdgMbAc+I3Na4qIiE2tbX5/KI1j2gA/Bc4AOgBLgDcx\nfQQiIuIBu8F/I3BwTPlgoumdiPWYVM/34VctcBwJwb+8vDxUV1dnszoiIiWnDjgi02+ym/ZZjknr\n9AbaAhcCryQc8zfgNEzncAfgZOCjxBPV1dURCoWK9jV58mTP66B70/3p/pp51dQQevttz+ub6Qso\nzyZ42w3+jcA1wAJMQH8BWAlcGX6BGQY6H3gPeAt4AovgLyLiiW3bwOeDykqYPdvr2rjGbtoHYF74\nFeuxhPKfwi8RkfwxZQpUVZntJUtgwABPq+MmJ4K/pMHv93tdhZwp5nsD3V+hs7y/1avhqKPM9ogR\n8MwzblYpL/i8rkCMUDh/JSKSG6EQDBkC8+eb8oYNcGDi6PTC4vP5IItYbjfnLyJSGObNg1atTOB/\n6CHzQVDggd8OpX1EpLh99RX88Idmu3Nn+PxzaN/e2zrlAbX8RaR4HXhgNPDfcgts2aLAH6aWv4gU\nn9deg0GDouU9e8xwTtlLLX8RKR6hkAnykcD/z39G35M4Cv4iUhzuuMN06ILJ7YdCcPrp3tYpjynt\nIyKFbft26NQpWv78c+jRw7v6FAi1/EWkcPXrFw38o0aZ1r4Cf1rU8heRwvPOO3DCCdFyYyOUlXlX\nnwKklr+IFBafLxr4//pX09pX4M+Ygr+IFIZHHokftRMKwXnneVefAqe0j4jkt++/hw4douW1a6F3\nb8+qUyzU8heR/DVkSDTwX3CBae0r8DtCLX8RyT9vvgkDB0bLDQ3Qtq139SlCavmLSH7x+aKB/4Yb\nTGtfgd9xCv4ikh8mTUru0L37bu/qU+SU9hERbzU0xK+0GQxCRYVn1SkVCv4i4p0uXWDr1mhZT/Nz\njdI+IuK+jz4yKZ5I4N+6VYHfZQr+IuIunw+OOcZsX3yxCfqxC7OJKxT8RcQd06cnd+g+95x39Slx\nTgT/wcAqYDUwrpnj+gGNwC8duKaIFIqmJhP0x4415dmzleLJA3Y7fMuAB4EzgY3AMuAVYKXFcVOB\n+YAeqSNSKo47Dt57L1pW0M8bdlv+/YE1wDpgN/A8MMziuNHALOArm9cTkULw73+b1n4k8H/+uQJ/\nnrEb/A8E1seUN4TfSzxmGPBIuKz/ASLFzOeDQw812xUVesBKnrKb9kknkD8A3BQ+1kczaZ+qqqq9\n236/H7/fb692IuKe556D4cOj5T179OD0HAgGgwSDQdvnsfuTGQBUYTp9AW4G9mDy+xGfxlynG7AT\nuBzTNxArFNKfhSKFJxSKPjgd4KmnYORI7+pTYnzmAzbjWG43+LcGPgbOAD4DlgIXk9zhG/E08Hfg\nZYt9Cv4ihWboUJg7N1rW77Drsg3+dtM+jcA1wALMiJ6nMIH/yvD+x2yeX0Ty0ebN0L17tLxmDZSX\ne1cfyVg+JeTU8hcpBLF5/MMPh7o67+oiWbf8NcNXRNLzzDPxgX/XLgX+AqZVPUWkZbFBf+RI06kr\nBU3BX0RSO/54ePfdaFmp2aKhtI+IJNuyxbT2I4E/GFTgLzJq+YtIvMSJWQr6RUktfxExqqvjA/93\n3ynwFzG1/EUkPuifcQYsWuRdXcQVavmLlLLzzkt+wIoCf0lQ8BcpRd9/b4L+3/5myrNmKcVTYpT2\nESk1rVrFB/oiD/rBoHlFtiOLBfv90e1SpOUdRErFRx9FH5wO8NVX0K2bd/XxgM9XfJ91Wt5BRFLz\n+aKB/+KLTQQsscAv8RT8RYrZmDHJHbrPPeddfSRvKOcvUoyamqB1zK/3ww/DVVd5Vx/JOwr+IsVG\nM3STVFfXMn16DdCaQKCRMWMqGTp0kNfV8pSCv0ixWLUK+vaNlteuhd69PatOvqiurmXs2AXU1d0O\nQE0N1NVNACjpDwCN9hEpBkXU2nd6aGYgMJGamtss3p/E/Pm3ZlXHfOLVYxxFxEtTp8JNN0XLe/Yk\nfxAUmNgg7/NFPwiy1dBgHebq68vsnbjAabSPSCEKhUxkjAT+SZOi70mcdu0aLd9v377J5ZrkF7X8\nRQrNfvvBt99GywWc4nHDmDGV1NVN2JvzBygvH8/o0YM9rJX38qmZoJy/SHPWrjUPTI9491049ljv\n6uMCp2bkVlfXMmPGQhYsKCMQaGL06LOKprM325y/gr9IISiiDt10RIZm1tS0prLSuaGZWt4hSmkf\nkXyW2KHb0ABt23pXHxdoaKY7nOjwHQysAlYD4yz2DwfeBd4D3gCK++9UEafEduieeqppshZ54AeY\nPr0mLj8PUFd3OzNmLPSoRsXJbsu/DHgQOBPYCCwDXgFWxhzzKTAI2Ib5oHgcGGDzuiLFq8RSPIk0\nNNMddoN/f2ANsC5cfh4YRnzwXxKz/RZwkM1rihSnzZuhe/doefZsGDbMu/p4xOmhmbGTxioqoKrK\nbJf6ev52g/+BwPqY8gbg5GaOvwyYa/OaIsWnxFv7sZwemtlSkC/Vh73YDf6Z/A89HRgJnJrqgKrI\nRzLg9/vxF/O/vAjAs8/CpZdGy9u2mXH8ecqNQBnp1J0xY1LM0MzBOevsdXpGca4Fg0GCDlTS7lDP\nAUAVJpcPcDOwB5iacNyxwMvh49akOJeGekppiW3td+8OX37pXV2y4MawSbeGZuZqaKkbvBrquRw4\nEugNfAZcCFyccMwhmMB/CakDv0jp6NQJtm+Plgus0VNsyyNraGn2zgY+xgT2m8PvXRl+ATwJfA2s\nCL+WpjhPSKSoffttKGRCvXndd5/XNcrYnDmvhsrLx8fdRnn5+NCcOa/m5Hq5DguNTY2hvr8ZHOKW\nViGqCHHM83vvKxCYmNuLO4TM0u97aYaviBuKpEPXjeWRc92v8D9b/4cbFt7Aix+9GL+jdgLUToTG\n9gBUVFQRDFbZv2COaYavSD76+9/h3HOj5fXr4aDCHe3sxhh8p0fZhEIhzn/xfF5e+XLc+3269eHB\nsx/k7qsWW36gFfuqnwr+IrlSJK39WIWyPPLKr1Zy9MNHJ71/6XGXcm/lvXTr0G3ve/Vj2pTkqp8K\n/iJOGzIE5s2Llosg6Efk8/LI4xaO4+7/d3fS+xcecyH//av/jqRHkrg9tDRfKOcv4pTGRmjTJloe\nOxYeeMC7+uRIviyPvOm7TfS8t6flvtd/9zqnHpJySlFKhbjqp5Z0FvFSEaZ4WuJFoHx42cOMmjsq\n6f2Tep3EGyPfoG2ZvYXvFPy9oeAvheftt+HEE6PllSuhTx/v6uMiNwLld7u+47Bph7F55+akfbN+\nPYtfHf0r29co9OUdFPxF3FaCrX03AuU9b9zDjYtuTHq/6z5dWTt2Lfu1y9/lL7yg4C/ilquvhkce\niZabmqCVE4/GKE27mnbR7rZ2lvuuH3g991Te43KNCovG+YvkWigUH+QvucQszCYZe3nly/zqL9Yp\nm1WjVvGjbj9yuUalR8FfJB0lmOJxUigUosMdHahvrE/ad0z3Y/jg6g88qFVpU/AXac7atXD44dHy\nG2/AKad4V58MBYPw+OO1vPVWDV9/3Zr992/k5JMrueKKQTnvzFy6cSknP2n9eI/FIxbj753jCkiz\nFPxFUimC1v6OHbUsXbqATz81k7K2bQOfbwLDh4N5uqqzBj09iNf+/ZrlvqZbmmjlU99IvlCHr0ii\nqVOjD04HaGgo2Aen53ohtnVb13HYtMMs9z0z7BlGHD/C9jWkeerwFXFCbGv/1FPh9de9q4sDcrEQ\n28CnBvLmhjct9+0Yv4MObTpkfW5xj4K/CBRFiseKEwuxbd65me73dLfcN/FnE7n1584s5SzuUvCX\n0rZ5s3mEYsTs2TBsmHf1cVi2C7GNmTeGGUtnWO5bO3YtvTv3drKa4gEFfyldRdraj6iurqWqqobt\n2zezzz4X4vP1pGfPH9C//2A6dozv7G1uotUBHQ9g0/Wb3KiyuEjBX0rPs8/CpZdGy9u2wX7FtWTA\nnXfWcs89C9iyJdri79JlApdddhY332wC/xP/eoIr5lxh+f1vjHyDUw4unCGtkjmN9pHSEtva794d\nvvzSu7rkkPUonxBtJnRkd5vvLb8nNFm/f4VIo31EmnPkkbBmTbRc5A2NvaN8egfht6fvfX93zDE/\nXn0ud118XdE/tESsqeUvxW3nTujYMVp+9FG48krv6uOCfe/Ylx27d1jua3/35dTvfHxvubx8AtOm\nBfQBUMC0qqdIoiLv0I1Y880ajpxxpPXO12+ERVMBaN36QhobX0g6xKkJX+INpX1EImpqIBCIlj/7\nDHpaP+4vF6qra5k+vYaGhta0a9fImDGVjresK5+tZOGnCy33fXPjN3TZp4t53OKShdRXVNG+fRMf\nf9yTdeuSj9+0KfsJX1K4nAj+g4EHgDLgSWCqxTHTgbOBncBvgRUOXFckmcet/erqWsaOXRA3rr6u\nbgKArQ+A5iZanXHYGSy6dFHS+0OHDoq7Zr9+Ey2Df9u26U/4kuJhN+1TBnwMnAlsBJYBFwMrY44Z\nAlwT/noyMA0YYHEupX0ke0OGwLx50bJH/5f69ZvI8uXJa+n06zeJpUszS61cNOsiXvgwOU0D8PE1\nH3PU/kdldD6rD6by8vFMmzZYOf8C5lXapz+wBlgXLj8PDCM++J8LzAxvvwV0Bg4AvrB5bRFobIQ2\nbaLliRPhVu/y17t3W/9K7doVn1qxSg1VDh5I29tSLyBndyhmJMDPmDGJ+voy2rdvYvRoBf5SZTf4\nHwisjylvwLTuWzrmIBT8xa487NA94ADrtXR69IimVuJa4APvhdOvp2b57bA8/nucekB5rMRUkJQu\nu8E/3d+2xD9JLL+vqqpq77bf78ef66dNSGF6+2048cRoeeVK6NPHu/rESGctnXOWV8BvAO5I+v49\nt+yJ/BkvYikYDBIMBm2fx+7/sgFAFabTF+BmYA/xnb6PAkFMSghgFVBBcstfOX9pWR629hNVV9cy\nY8bCvamVfpfsw23hTt8kr4+DRXcBUF5exZo1Ve5VVIqCVzn/5cCRQG/gM+BCTIdvrFcwHb7PYz4s\ntqKUj2TqF7+AOXOi5aam+Iep55GhQweZ1n3Ygrr4/We8eSP/mJ88KO6IIzTqRtxjN/g3YgL7AszI\nn6cwnb2RKZSPAXMxI33WADuA39m8ppSSUCg+yJ92Grxm/ZhAL328+WP6PGSdejqp10ksu3zZ3nJ1\ndS3rVme+zLKIk/Ipuai0j8TL8xTPjx78EZ98/Ynlvo3XbqTXD3ql/N7E1NDo0WepI1ayouUdpHis\nWgV9+0bLNTVw1lne1Sdsa/1WukztYrmvzFdG4y3WI31EcknLO0hxyLPW/uWvXM6TK5603PfWf7xF\n/wP7u1wjEWco+Et++MMf4IEHouX6emhn/WSpXGra00TrW1P/WmjNeykWCv7ivdjWfrdu8NVXrl7+\nkWWPcPXcqy33Pfu/nuWSYy9xtT4iblDwF+94mOLxTUmdIm26pYlWvvwcRiriFAV/cd+mTfFLLD/z\nDIwYkdNLBtcFOX3m6Zb7rht4HX+q/FNOry+SbxT8xV0utvaba91/e/O37Nt235xdGyAYhMcfr+Wt\nt2r4+uvW7L9/IyefXMkVVwxCK5eI1xT8xR3TpsHvfx8tf/MNdLEeNpmttVvWcvj0wy33DThoAEsu\nW+Lo9VqyY0ctS5cu4NNPzWSubdvA55vA8OEAGtMv3tI4f8m9HLb2+z3Rj+WfLbfct27sOg7tfKhj\n18pUIDCRmprktf312ERxksb5S/7JQdDfVr+NzlM7p9yfT0MxGxqsf73q67N7bKLSSOIkBX9x3rZt\n0DkmQN96q3nISpbO/8v5vLTyJct9r/72VQYdmp8plHbtrGf8tm+f3QJuSiOJkxT8xVkOtPb3hPZQ\n9sfUreN8at03J521/TMxfXpN3LkA6upuZ8aMSVoXSDKm4C/OeOklOP/8aHnNGigvT/vb71tyH9fV\nXGe577bTb2PCoBTr4eexjh0H0b8/wCQ2by6jW7cm+vcfTMeO2QVqp9NIUtoU/MW+LFv7zQ3F3D1p\nN61bFfZ/T78f/P5BOJWScTqNJKWtsH+7xFtHHmla+BEtBP1Fny7irGetV+c8+4izmTt8bsrvDQbN\nK7Id6eA0ATbN+hY4p9NIUto01FMyt3s3tG0bLY8cCU89ZXloc637r2/8mq77dM348j6f54t9eiby\nHIAFC8oIBPQcANFQT3FLCyme5p5oBfY6a6ura5k+vQZoTSDQyJgxlSUV+IJBWLZsEAMGDKK+HgYM\ngGXLoGPH0vnrR5yjlr+k5/XX4Wc/i5ZXrIDjjwfgkPsPYf329Zbf9t7/eY+fHPAT25evrq5l7NgF\nCSmPCUybFiipDwCRRHqSl+ROQmt/R8N37Htn6nVxcjEUU7NlRawp7SPOGzUKHn54b7HVLRBqBSQE\n/hd//SLnH30+uaRhjiLOUvCXZKEQtIquZz+rL/z6woRDHG7dtzSaR8McRZyltI/s9eKHL/LrH18Q\n956vynwdd+o47jrzLlfqYTWaxzrnP55p0wYr5y8lTTl/yUpkKGbvLbB2WvT9gZdB7WO7aFPWxv06\npRjKqWGOIsm8Cv5dgReAQ4F1wAXA1oRjDgb+DPwQCAGPA9MtzqXg74IVn6/gp4//NO69UFXCQR7/\nHFoax1/K4/xFEnkV/O8GNoe/jgO6ADclHNMj/HoH2Bf4F3AesDLhOAX/HOl0Vye2N2xPev/G12Hq\nopg3GhriJ2+5LDKOv6amNZWV8eP4NcNXxJpXwX8VUAF8gQnwQSD1DB9jNjAD+EfC+wr+Dvniuy/o\ncW8Py329ftCLjddujB++ecop8MYbLtXOmsbxi2THq+C/BdPaj5zrm5iyld7Aq8AxwHcJ+xT8bTjv\n+fP428d/s9z3/lXv8+Mf/tgUcvwM3Wxb6BrHL5KdXI7zX4hp1SdKXGM3FH6lsi8wCxhLcuAHoKqq\nau+23+/Hr7/nU9rVtIt2t7VLuT9pKObmzdC9e7Q8ezYMG+Z4vWKDvM8X/SBoicbxi6QnGAwSTPcX\nqxnpBH/rZRiNSLpnE9AT+DLFcW2Al4D/wqR9LMUGf0n2nyv+k8teucxy3+wLZzOsT4pg7lFrPxMa\nxy+SnsSG8ZQpU7I6jxMdvl8DUzEdvZ1J7vD1ATPDx/2hmXMp7ZMgFArR6o+tUu7fc8ueyJ981p59\nFi69NFretg3228/BGiaLHYmTyagcjeMXyY6XQz3/AhxC/FDPXsATwFDgNKAWeI9oWuhmYH7CuRT8\ngfe/eJ9jHz3Wct+Tv3iSy35q3fJPEvuh0L07fJnqjzJn+XwwZ07qUTvN0Th+kcxpklcBm/r6VG76\nR+IfTEbDxAbalmUw/DLDB6w4zeerpbzc3qgdjeMXSZ+CfwHZ8v0Wetzbg11Nu5L2zRs+j8FHZPFk\npp07zcLuEY8+CldeaaOW2fH5JgL2Ru0o+IukT6t65rnn3n+O4S8PT3r/sM6H8cHVH9ChTYfsT57j\nDt10RB+0ssFyf0ujdmI7jSsqINL3r0lcIrmh4J8j9Y31nPDYCazavCpp31PnPsXIE0bav0hNDQQC\n0fJnn0HPnvbPm6H4ztqJlse0NGpHQV7EXQr+Dlq8djE///PPLfd9dcNXdOvQzbmL5UFrP2L69JqY\nHH8lZgqIHjIuks8U/G3YE9rDsOeHMeeTOUn7JldMpspf5fxFR4+GBx+MlvMgOR4/QSvSqTsJWE8g\ncDCjR2u4pki+UfDP0AdffsBPHrF+Ju2a0Wso71qemws3NUHrmB/Xww/DVVfl5loZSp6gNSj80tIM\nIvlKwT8N81bPY9TcUazdujbu/RHHjeDpYU83P9HKCXmU4rEyZkwldXUTkiZo1dUp1SOSrzTU08L2\nhu3c+uqt/GnJn+LeP6nXSTw05CH6H9jfnYqsWgV9+0bLa9dC797uXDtDVhO0zjlnUL59TokUHY3z\nt2nJ+iVcM+8a3v787b3v7dN6Hx4a8hAjjh9BK1/qZRZyIs9b+6lku7yDiGRHwT9D9Y313LfkPib8\nM35x0qFHDuW+wH0ctf9RrtUlznPPwfCY+QB79iR/EOQxnw8mTzbbeuiKSO4p+Gcgce37eyvv5Zr+\n12S2jILTQiFoFfPXxVNPwUgH5gK4TK19EXcp+Gfg3U3v0rinkRN7nejK9Vp0550wfrzZ7tIFvvnG\n2/rYoOAv4i4F/0K0fTt06hQtb9oEBxzgXX2ypOfrinhHwb/Q9OsHy5eb7VGj4iduiYikSQu7FYp3\n3oETToiWGxuhTI8qFBF3uTx+scT5fNHA/9e/muS4Ar+IeEDB3w2PPRY/XDMUgvPO864+IlLylPbJ\npfp62GefaDmPZ+iKSGlR8M+Vc86B6mqzfcEF8MIL3tYHjcoRkSiN9nHa6tVwVMzs4IYGaOvh5LEU\nNB5fpDhkO9pHOX8n+XzRwD9zpomueRj4RUSU9nHCCy/ARRdFy2pSi0ieU/BPQ6pc+ek/a6TizDbR\nAz/8EI4+2t3KiYhkwU7OvyvwAnAosA64ANia4tgyYDmwAfhFimMKIue/N1f+u9/BM8+YN/1+WLzY\nw1plTjl/keLgxfIOdwObw1/HAV2Am1Icey1wIvAD4NwUxxRE8D/It4ENHBx9Y8cO6NDBuwplScFf\npDh40eF7LjAzvD0TSDVr6SBgCPAk+TW6KHPdukUD/wMPmOhZgIFfRMROzv8A4Ivw9hfhspX7gRuA\n/Wxcy1uJwzeLpMlcVWW+Wo3515wAkeLWUvBfCPSweH9CQjkUfiU6B/gSWAH4W6pMVSQaAX6/H7/X\nUSYUgiFDYP58U964Ed+BvSxv1K5Mgq1TgTnyz+3zRc8XEXsuq/0i4o1gMEjQgV9IO2mYVZiAvgno\nCSwG+iQccwfwG6ARaI9p/b8EXGpxvvzK+dfUQCBgth96CK6+GnAnV57JNbKtTybP2lX/gEj+8mJJ\n51eAEcDU8NfZFseMD78AKoDrsQ78+WPnTujZ0zxopUcPsx5P+/auV6O5lIwd1dW1TJ9eA7QmEGhk\nzJhKYJC9k4pIwbE71PMvwCHED/XsBTwBDE04vgK4jnwe7XP//XDttWZ78WLLSOt2y9/JVnl1dS1j\nxy6gru72ve+Vl0+gri5AKJT6A0Atf5H8pSd5pSFVrjzQ998MvOhQU/jlL2HWrPglmGMUcvAPBCZS\nU3ObxZ5JhEK3OnINEXGXnuSVhqROzMUhuPhimBJecfPTT+Gww7yqXs41NKT6ceuBMiKlpmQXdjuN\n16BVK7Muz113maZtEQd+gHbtGlPsaXK1HiLivdIM/oMH8xqDzINWvv0Wxo1r8Vuqq2sJBCYCVQQC\nE6murs19PR02Zkwl5eXxo3TLy8cDZ3lTIRHxTEmlffY65hiGLBjD3J1D0jo8saO0pgbq6kwQHTq0\ncEbKROo6Y8YkFiwoIxBoYvTowZxzTuHcg4g4o6Q6fGM50VEaCExi/vzUHaXZcmMMvsb5ixQHPcwl\nh1J1lNbXq6NURAqTgn8aUnWUtm+vjlIRKUwK/mlI1VE6erSzHaWZdCo70QHd0jmKoZNbRPJfyE2Z\nXm7OnFdDgcDEEEwOBQITQ3PmvOpofebMeTVUXj4+ZLLr5mXKyddJdWwmdYLmz+HENUQk97BeVLOg\nuPwP5u73taSyckJcoI2+JqZ9bCCQfGwq0Pw5nLiGiOQeWQZ/pX3yROrZt+uTUi7OdEA3fw51cosU\nt9Ic55+HUs++PZiamlvj5hU40wHd/DnUyS1S3Equ5Z+vnZhWncpmNWzTqVxXdzszZixMeWzmHdDN\nn8OtTm4R8UZJtfzzeaZu/Ozb9cDBwGBi19qPpFxSzdTN7B4GMW1a6nM4cw0RyVclNcPXiZm67izp\nPBFIr56Z1CfVktZTpmiGr0ih0pLOaSicTszK8ENWYh+6Mp7RowfbOmuqJ4FNmWLrtCJSgEoq+BdO\nJ2bzKRkREbtKKu1j/RjD8Uybln5gLeQneWV7DqV9RPKXHuOYpurqWmbMWBjToj4roxa1gr+I5BMF\n/ww5sRRyrij4i0i6FPwz5MQomVQdqHYp+ItIuhT8M5TPAU3BX0TSpYe5iIhI2uwE/67AQuAToAbo\nnOK4zsAsYCXwETDAxjVFRMQBdoL/TZjgfxTwj3DZyjRgLtAXOBbzISAiIh6yk/NfBVQAXwA9gCDQ\nJ+GYTsAK4PA0zqecf5gbOf9MOrHz+d9KpNR50eG7BegSc55vYsoRxwOPYdI9xwH/AsYCOy3OV9LB\nP5N1d9wYfeT2CCcRyU6ugv9CTKs+0QRgJvHB/htMP0Csk4AlwCnAMuABYDtwi8U5Q5MnT95b8Pv9\n+HMYZfIt+KdSKPUUEXcEg0GCkZYZMMUszuVqy38V4Ac2AT2BxSSnfXpggv9h4fJpmL6BcyzOV9It\n/1QKpZ4i4g0vVvV8BRgBTA1/nW1xzCZgPaZT+BPgTOBDG9csCbEpl4oKqKoy20q5iIhT7LT8uwJ/\nAQ4B1gEXAFuBXsATwNDwcccBTwJtgTrgd8A2i/Op5S8ikiHN8M2Qgr+IFAMF/zRoBIuIFBsFfxGR\nEqS1fUREJG0K/iIiJUjBX0SkBCn4i4iUIAV/EZESpOAvIlKCFPxFREqQgr+ISAlS8BcRKUEK/iIi\nJUjBX0SkBCn4i4iUIAV/EZESpOAvIlKCFPxFREqQgr+ISAlS8BcRKUEK/iIiJUjBX0SkBCn4i4iU\nIDvBvyvlVZHEAAAEeElEQVSwEPgEqAE6pzjuZuBD4H3gOaCdjWuKiIgD7AT/mzDB/yjgH+Fyot7A\n5cBPgZ8AZcBFNq5ZsILBoNdVyJlivjfQ/RW6Yr+/bNkJ/ucCM8PbM4HzLI7ZDuwGOgCtw1832rhm\nwSrm/4DFfG+g+yt0xX5/2bIT/A8AvghvfxEuJ/oGuBf4N/AZsBVYZOOaIiLigNYt7F8I9LB4f0JC\nORR+JSoHfo9J/2wDXgSGA/83o1qKiIijfDa+dxXgBzYBPYHFQJ+EYy4EzgL+I1z+DTAAGGVxvjWY\nDwsREUlfHXCEmxe8GxgX3r4JuMvimOOAD4B9MB80M7EO/CIiUiC6YvL3iUM9ewHVMcfdSHSo50yg\njYt1FBERERERr6U7SawzMAtYCXyE6TPId+neG5i5DyuAv7tQL6ekc38HY/qBPsSk/sa4VrvsDcb0\nZa0mmtJMND28/13gBJfq5ZSW7m845r7eA94AjnWvaral87MD6Ac0Ar90o1IOSuf+/JhY8gEQdKVW\nWbobkxICczNWfQZgUkUjw9utgU45rpcT0r03gGsxo59eyXWlHJTO/fUAjg9v7wt8DPTNfdWyVoYZ\ndNAbk5p8h+T6DgHmhrdPBt50q3IOSOf+BhL9/RpM4dxfOvcWOe6fwBzgV25VzgHp3F9nTEProHC5\nm1uVy8YqonMDeoTLiToBn7pWI+ekc29gflCLgNMprJZ/uvcXazZwRs5qZN9AYH5M+SaSZ60/ihnB\nFhH775Dv0rm/WF2ADTmtkXPSvbffA1cDT1NYwT+d+7sa+GMmJ/VyYbd0JokdBnyF+WG9DTyBmSWc\n79K5N4D7gRuAPW5UykHp3l9Eb0yK5K0c1smuA4H1MeUN4fdaOuYgCkM69xfrMqJ/5eS7dH92w4BH\nwmWreUn5Kp37OxKTjl0MLMcMq29WS5O87LI7Saw1Zl2ga4BlwAOYT7xbHKxjtuze2znAl5gcnd/R\nmjnD7v1F7IvpsxkLfOdM1XIi3WCQODemUIJIJvU8HZNqPTVHdXFaOvcWiR0hzM/Qzhwnt6Vzf20w\nsfIMTAN5CSZttzrVN+Q6+J/VzL4vMMElMknsS4tjNoRfy8LlWTT/p6qb7N7bKZj1kYYA7YH9gD8D\nlzpbzazZvT8w/yFfAv4Lk/bJZxsxndQRB5Oc9kg85iAKZ62qdO4PTCfvE5ic/xYX6uWEdO7tROD5\n8HY34GzMumOF0NeWzv2tBzYD34dftZh5VimDv5fSmSQG5iaOCm9XAVNzWy1HpHtvERUUVs4/nfvz\nYT7M7nerUja1xsyU7A20peUO3wEUTocopHd/h2A6FgthRF2sdO4t1tMU1mifdO6vD6b/sAzT8n8f\nONq9KmYm3Ulix2Fa/u8CL1MYo33SvbeICgqjBRKRzv2dhunLeAeT2lqBaU3ms7Mxo5LWYJ5DAXBl\n+BXxYHj/u5g/swtJS/f3JPA10Z/XUrcraEM6P7uIQgv+kN79XU90Qm0hDK0WERERERERERERERER\nERERERERERERERERESld/x9OFJkNiIeNbwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106616910>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sig = 0.05\n",
      "lnlike0 = (1/sig**2)*np.sum((x - (1.1*x+0.1))**2)\n",
      "lnlike1 = np.sum((x - (a*x+b))**2/err**2)\n",
      "print lnlike0, lnlike1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "90.7296 8350.74328904\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    }
   ],
   "metadata": {}
  }
 ]
}