ipython notebook: least-squares special cases

Least-squares special cases.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.optimize as spo"
     ],
     "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": "code",
     "collapsed": false,
     "input": [
      "xi = np.array([25,16,11.11,8.16,6.25,4.94,4,2.78,1.78,1])\n",
      "ci = np.array([44,18,17,6,8,9,9,11,3,3])\n",
      "obs_unc = np.sqrt(ci)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Standard least-squares fit. Note that this is the correct use of w (I had been squaring the number before, which isn't necessary). Polyfit returns coefficient of highest power first, so it's 1.22*x+1.52"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.polyfit(xi,ci,1,w=1/obs_unc)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "array([ 1.22327698,  1.51995345])"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.errorbar(xi,ci,obs_unc,marker='s',ls='None',color='k')\n",
      "xfit = np.arange(0,30,0.5)\n",
      "plt.plot(xfit,1.52+1.22*xfit,marker=None,ls='solid',label='Std least-sq')\n",
      "plt.plot(xfit,2.54+1.35*xfit,marker=None,ls='solid',label='Gaussian, analytic')\n",
      "plt.plot(xfit,2.10+1.32*xfit,marker=None,ls='solid',label='Poisson, analytic')\n",
      "plt.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "<matplotlib.legend.Legend at 0x105724d90>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VNX28PHvmYQUqYEQWigJBEKvSgIqQYqIiBTlSgel\niQiiiHr1J3Av+nq90kEjzQsYLIggSBEViZRgkA5BIBSlS5UeUtb7x3YOqSQhCQlkfZ5nHpKZM+fs\nM0zW7Nln77UsERGUUkrd8xy53QCllFJ3hgZ8pZTKJzTgK6VUPqEBXyml8gkN+EoplU9owFdKqXwi\nQwH/7NmzvPTSS/j5+eHu7k7JkiVp3rw527ZtAyA2NpYxY8bg7++Pm5sbvr6+DB8+nMuXL+do45VS\nSmWcld48/LNnz9K4cWMOHjyIq6srAQEBuLq6cvjwYT755BM6depEz549CQsLw8XFhYCAAA4cOEBs\nbCzNmjVj9erVWJZ1p85HKaVUGtLt4b/11lscPHgQX19f9uzZw+7du9m+fTvnz5/nscceY8uWLYSF\nhQEwadIkoqKiWLhwIQDh4eEsXrw4Z89AKaVUhtwy4IsIX375JQB+fn48/fTTFCxYkJo1azJ9+nQ8\nPT1ZsWIFAJZl0blzZwDatm2Lu7s7ACtXrszJ9iullMog11s9ePr0ac6fPw/A2rVr8fHxoXTp0uzZ\ns4fBgwcTHx/P0aNH7e19fHwAcDgceHt7c+zYMY4cOZKDzVdKKZVRtwz4cXFx9s/e3t4cOHAAT09P\nHnroISIiIpg6dSohISGpPvdWlwZ0TF8ppW5PVtKf3XJIp2TJkhQoUACAqlWrUrBgQRwOBw0aNADg\n8OHDlC9f3t7+1KlTACQkJHD27FkAKlSokGaj79XbqFGjcr0Nen56fvnx/O7lcxPJep7LWwb8AgUK\n0Lx5cwD27t3LlStXSEhIYOvWrQBUq1aNNm3a2AHcebF22bJlxMTEANiPK6WUyl3pztIZO3Ys7u7u\nnD17Fn9/fypXrsyGDRuwLItRo0bRoEEDunbtCsCwYcOoXr26ffH24YcfpkOHDjl7BkoppTIk3YDf\nqFEj1qxZQ4sWLbh27RqXLl0iJCSE1atX07FjRwDmzJnD22+/TYUKFTh06BA+Pj4MHTqUb7/9NsdP\nIC9K67rGvULP7+52L5/fvXxu2SHdhVc5clDLypbxKKWUyk+yGjs1l45SSuUTt5yWqdSdVLx4cXvd\nh1L5lZeXF+fOncuRfeuQjsoz9H2h1K3/DnRIRymlVIZowFdKqXxCA75SSuUTGvCVUiqf0ICvVB5V\nqVIlHA4Hffv2TXObNWvW4HA4cDgchIeH38HWqbuRBnylsig2Npb33nuPGjVqULBgQYoUKULlypXp\n0KEDmzZtsrcbPXq0HZwzIyPZZS3LyrUstBn5YFJ5gwZ8pbJo5MiR/POf/+S3337D19cXf39/zp49\ny9KlS9mzZ0+K7TU9uMotGvCVyqLPPvsMgLfffpu9e/eybds2Lly4QEREBA888ABgcrz861//Akxm\nWWdPf86cOQDs3r2bBx98EA8PD6pXr57l0qDh4eG0bt2aokWL4uHhQcOGDe1stk6vvvoqNWvWpFix\nYhQoUICyZcvSp08fTp48aW9z6tQpevbsSdmyZXF3d8fHx4eHHnqIsLAwDh8+jMPh4I8//gBMTq2M\nfIOJjIykVatWeHt74+7uTvny5WnXrh2bN2+2t1m0aBFVq1bFw8ODZs2asXz58hSvmboNkgty6bAq\nj7tb3xc+Pj5iWZYEBwfL0qVL5cSJEym2GTx4sPj6+oplWfa2wcHBsnz5crl27ZqUL19eLMsSNzc3\nqVWrltx3333i4eEhlmVJ37590zz2Tz/9JJZlicPhkPDwcBERWbFihTgcDrEsS8qVKyfVq1e3jztn\nzhz7ubVq1RIvLy+pU6eO1KhRw37OAw88YG/TqVMnsSxLChcuLI0aNRI/Pz9xdXWV/v37y4kTJyQo\nKEjc3d3FsiwpWbKkfV5piY+PF29vb7EsS0qXLi0NGzaUMmXKiMPhkLCwMBER2bFjh7i4uIhlWVK0\naFGpXr26FCpUKNVzuBfd6u8gq38jGvBVnpHa+wLu/C2zRo8ebQcj561atWry9ttvy9WrV1Ns53A4\nkjx/9uzZ9v1LliwREZEffvjB3ldmA369evXEsixp3769xMfHi4jIsGHD7A8Apx07dkhCQoL9+4wZ\nM+x9HTx4UETMh4JlWfLpp5/a2505c0Z27Nhh/16xYsV025n4uc5jHD161L4/Ojpa/vjjDxER6dWr\nl1iWJUWKFJFjx46JiMhbb72lAT+dxzJCh3RUnpYbIT+zRo0axddff82TTz5J0aJFsSyLffv28e9/\n/5tevXolOpfUd75z504A3N3deeKJJwBo0aIFXl5emW7LjRs32L59OwBLly7F1dUVh8PB5MmTAThx\n4oRdZ3rbtm00atSIQoUK4XA4GDBggN3O48ePA9C+fXsAevfuTeXKlWnbti0ff/wxZcqUSbctHTt2\nJCgoyL4tX76cEiVKEBwcjIgQEBBArVq16NKlCz/99JO9T+fr0aRJE8qWLQvAP/7xj0y/FiolTZ6m\nVDbo0KGDXexny5YtDBgwgC1btuRqTYgKFSrYATOx+Ph41q1bR+/evQFTr7pWrVpcunTJvsgcHx8P\nwDvvvEPTpk1ZtWoVu3btYv369axcuZIFCxbYle/Ssm3bNnt8H+DMmTMA/Pjjj8yfP58NGzYQFRXF\n4sWL+eqrr9i1axcTJ060t0/8AZnWh6XKHO3hK5VFb731lt2rBmjQoAHVqlUDoHDhwvb99913n/3z\n1atX7Z9r164NwPXr11m6dCkAq1evvq3MoW5ubtStWxeAKlWq8NNPP7FhwwY2bNjAwoULeeWVV6hU\nqRK//PILYGYM7dy5k40bN9KzZ0/7Pqd169bRrFkzJk6cyA8//MDHH38MwI4dO+z2Oc/r8uXLSdpy\n6NAh4uPj7Zvz28769evp06cPM2fOZMOGDTz77LMA/PTTTwDUqVMHgA0bNnDixAkAFixYkOnXQqUi\nSwNCtymXDqvyuLv1fVGqVCmxLEu8vb2lQYMGUqFCBXu8eejQofZ2S5Ysse+vWLGiBAUFycGDB+X6\n9ev2RdsCBQpIzZo1xdPTU9zc3DI8hm9ZVpKLts6LnsWKFZP69etLuXLlxLIsqV+/voiIfP/99/bz\nSpQoIYGBgVKiRIkU+2ratKm4ublJ5cqVpUGDBuLp6SmWZUmFChXsNjgv7Lq4uEiDBg1u2d7Y2Fh7\nfL5GjRpSq1Yt+2Jxjx49RERk586d4urqal8sDgwM1Iu2GXgsI7SHr1QWvfPOO/b4/b59+/jzzz+p\nWrUqb775Jh988IG9Xbt27ejfvz8lSpTgyJEjREZGcu3aNdzd3Vm+fDlNmjTBxcWFmJgYPvnkE8qW\nLZvunH3n44m3a9OmDatXr+bRRx/Fsiz27NmDm5sbnTp14q233gKgZcuW/Oc//6Fs2bLExMRQo0YN\nPvrooxT7euaZZ2jcuDFXrlxh9+7dFC1alA4dOrBixQp7m7FjxxIUFIS7uzvbtm1j165dabbX1dWV\nQYMGUblyZU6ePEl0dDQVKlTg+eefZ9q0aQDUqlWLBQsWEBAQQFxcHN7e3nz++ecZ/e9Qt6D58FWe\noe8LlZbDhw/j7+8PwP/+978kF8PvNZoPXymlVJZpwFdK3TU0LUXW6JCOyjP0faGUDukopZTKBhrw\nlVIqn9CAr5RS+YQGfKWUyifSDfiJq/QkvyUkJACm4s+YMWPw9/fHzc0NX19fhg8fnmKptVJKqdyT\n4eRpJUuWpHLlyknuc06RevbZZwkLC8PFxYWAgAAOHDjApEmT2LZtG6tXr9apVEoplQdkeEjn8ccf\nt5MwOW+WZbFlyxbCwsIAmDRpElFRUXZlnfDw8CxX7lFK5SznN/YxY8bkdlPuuD59+uBwOPDz88uW\n/TmrgOXVylwZDvhfffUVnp6elClThnbt2rFt2zYAO6eGZVl07twZgLZt2+Lu7g7AypUrs7vNSuU5\nN27cYNq0aTz88MMUL14cNzc3ypUrx4MPPsjYsWO5dOlSbjcxTY0bNyYoKIjy5cvndlNyTWZHIdIq\nSO/h4WG/nj4+PtnZxGyR7pCOZVm4uLhQpkwZXF1d2bNnD8uXL+fHH38kIiLCLqYA2CfocDjw9vbm\n2LFjSR5X6l504cIFWrZsyZYtWwDw9PQkMDCQa9eusXnzZiIiIujUqRM1atTI5ZamLiIiIrebcNdK\n/kFRunTpPP16ptvD79atG6dPn2bv3r3s3r3b7rHHxMQwbdq0ND8Z01sNNnr0aPu2Zs2azLdcqTxi\nyJAhdrAfMmQIZ8+eZceOHezfv5+//vqLhQsXUrJkScDkwe/QoQN+fn4ULFgQd3d3qlatyqhRo4iN\njbX3GRISgsPhoHnz5vZ9qfUq0ysIfuXKFV544QUqVKiAh4cHJUqUoHHjxkyYMMHeR/IhnStXrmS6\njdOmTaNSpUoULlyYJ554glOnTmX6dVy5ciUPPfQQPj4+uLm5UaRIER5++OEkowSJh0zGjx9Pjx49\nKFy4ML6+vrzzzjtJ9tezZ08CAgIoXLgwbm5uVKxYkWHDht3y21bXrl1xOBwEBQUluf/hhx/G4XDQ\ntWvXWxakT2tIJzo6mh49elCmTBnc3NwoW7YsgwYNSvc1WbNmTZJYmWW3k1PZmTe7devW8s4779g1\nKp3Fm+Pj4+0CzAMHDkzx/Ns8rLrH3Y3viwsXLti52+vVq5fu9qdPnxbLsqRMmTIpcue/+uqr9nbN\nmjUTy7KkefPm9n2jRo1KUhM3IwXBX375ZbEsSzw8PKRhw4ZSpUoVcXNzk1atWtn7dR5/zJgxt9VG\nNzc38fT0lGrVqtnbde/ePdOv5X//+19xc3OTKlWqSMOGDaVIkSJ2jYDt27eLiMihQ4fsY7i5uUm5\ncuXsIvKWZcn3339v769QoUJSsmRJqV+/vlSpUsXe5umnn7a36d27t1iWJX5+fiIism7dOvs13r17\nt4iInDhxQhwOhzgcDvnuu+9uWZD+8OHD9vOdefv3798vxYoVE8uyxNXVVWrUqCG+vr72MZO71d9B\nVv9G0n32+PHj5fjx4/bvq1atsk904MCBsnnzZvv3qVOnikjSQg+LFi3K9kare1Nq7wtGc8dvmREZ\nGZlqsZOBAwcmKWo+YsQIERG5ceOGREVFJdlHjx49xLIsKV++vH1fRgJ+RgqCt2vXTizLkrFjx9qP\nX7x4UX799Vf79+QBP7NtdHV1tYuaO4uhlClTJjMvo4iIHD58WC5cuGD/fu7cOSlcuLBYliX/93//\nJyJJA37Tpk0lNjZWzpw5IwUKFBDLsuT111+3n79t27Yk+3/zzTftD4qYmBgRSRnwRUTq1q0rlmXJ\nyy+/LCIiH374oViWJb6+vnbR97QK0idunzPg9+3b1z7u2rVr7W23bt2a6uuQqwG/YsWK4nA4pEKF\nClK9enX7ZAoXLix79uwREZFu3brZFW8CAwPtF79Zs2Y50mh1b7ob3xdpBfz3339fgoKCUvSM4+Li\n5I033pCAgAC7opXz5urqaj8/IwFfRKRJkyZiWZZ4enpKzZo15emnn5YZM2ZIbGysiIhMnz7d3n/5\n8uWlRYsW8tZbb8nhw4ftfSQP+Jlto7OKlsjNoJo8EGbE3r175cknn5SSJUvaVbCct379+olI0oA6\nadIk+7nOil7PPvusfd8HH3wgNWvWtEcbnLfEH5CpBXzna1aqVCmJjY2VRx55JMWHSWr/F8nb5wz4\nNWrUEMuy5JFHHsnQ65CTAT/dMfw333yTFi1aEB8fz+HDh/Hz86NHjx5s3ryZwMBAAObMmcPbb79N\nhQoVOHToED4+PgwdOjRXCzgrdScEBgbi6mrmPqxfv96+/9VXX0314t17773He++9R3R0NOXKlSMo\nKIhy5coB2AsZ4ebFQGcxcYC//vorxf5+/PFHZs6cSbdu3ShSpAiLFy9mwIABjBgxAoD+/fsTHh7O\niBEjqFmzJtu3b7cLk1+7di3Vc8poG52KFStm/+x8LW7H448/zpIlS/jrr7+oW7cuQUFBuLm5pXgd\nbnVc+fvaYVhYGK+++ipRUVH2dQtnAZW09ufUvXt3ihYtyunTp5k1axbh4eFYlkWfPn1u+9wSty03\npRvw+/fvz6pVqzh69ChXr17lwIEDzJ07l4CAAHsbV1dXRo8ezcGDB7l+/TpHjx5l4sSJFCpUKEcb\nr1RuK1y4MM888wwAW7Zs4fXXX09yYTO5jRs3AlCtWjUOHjzI2rVr7aLdiZUqVQowFykTEhK4fv16\nkrKCTukVBI+MjKRmzZq8//77rFixwi6Sfvz4cX777bcstTGj0prCmNjZs2c5cOAAAP/617/YsmUL\n8+fPv+1jOs+hSJEiHDp0iIiICFq1apWh595333307t0bEWHEiBEkJCRw//3324Xpnds4JS5Inxrn\nBeB169axYcMG+37n1PY7SXPpKJVFU6ZMoWHDhgC8//77FC9enHr16lGxYsUU29atWxeAvXv34ufn\nR8WKFfnll1+ApD3AFi1aAHDkyBHq169PrVq1iI6OTrJdXFwcrVq1wsvLi5o1a1K7dm1mzJgBYAfo\nyZMnU6pUKfz8/GjYsCGPPvooAIUKFUqxcj6zbcysW811L1GiBL6+vgC8/fbb1K5dm4YNG1KgQIHb\nOpbzHC5evEilSpXw9/dnwYIFQMbOYfDgwViWxZUrVwBS9O6rV69u76tGjRoEBwdz6NChVPf1z3/+\nk2LFihEXF8fDDz9MjRo1qFChAp06dbqtc8sKDfhKZVHRokVZv349EyZMICgoCBcXF/bu3YtlWYSE\nhPDBBx8wfPhwwPzx9+7dm2LFinH58mW6devG4MGDgaQBsW/fvgwdOhRvb2+OHj1Ky5YtGTZsWJLt\nMlIQvF27doSEhBAbG8vu3btxd3enVatWrFixgiJFitjHS3zsjLbRsqwUQTy1oH7u3DnAFCe/lYUL\nF3L//ffj6uqKiDB//ny8vb1T3Wdqx01833PPPcfLL7+Mt7c3V65c4ZFHHrGnUiY/h9RUrVrV/tD1\n8PCwv8U5pVWQPrX9Vq5cmU2bNtGtWzd8fHw4cOAACQkJ9ofvnaQVr1Seoe+Le1O9evXYtWsX69ev\np3HjxrndnAwbOnQoU6dOpUuXLnz++ed37Li3+jvI6t/I7V9hUUqpdJw7d46dO3cycODAuybYT58+\nnWXLlrF8+XIcDod9AfxeoAFfKZVjihcvfssZMXlRREQES5cuxdfXl9GjR9OoUaPcblK20SEdlWfo\n+0IpLWKulFIqG2jAV0qpfEIDvlJK5RMa8JVSKp/QgK+UUvmEBnyllMonNOArlcfk56Li2SG1amFZ\nsWbNGvv/JDw8PFv2mVs04CuVRc4A47y5uLhQrlw52rdvf1v1TbWoePa4VbK21PTp0weHw4Gfn1+S\n+4sWLWr/nxQtWjQ7m3jH6UpbpbKJm5sbDRo0ICYmhp07d/Ltt9+ycuVK1q9fz/3335/h/eTlItj5\nQfIPivr1698z/yfaw1cqm5QtW5YNGzawefNmFi9eDJgUxonzui9ZsoSHHnqIQoUK4eHhQb169QgN\nDU2yn+RDOvHx8bz11ltUqVIFT09PvLy8qFevHm+88Yb9nJUrV9K0aVO8vLzw9PTEz8+PTp06cfjw\nYXubdevW8eijj1K0aFHc3d0JDAzk3XffJS4uzt6mUqVKOBwOevfuzahRoyhTpgxeXl707NmTy5cv\nZ/o1mTdvHg888ADe3t4UKFAALy8v2rRpw6ZNm+xtEg+ZzJkzh3bt2nHffffh7+/P7Nmz7e0yWgA+\nueDgYBwOR4qMlxUrVsThcPDGG29QqVIl5s6dCyQtlB4eHp7mkM7mzZvp0KGDXUC+UqVKdkbOPCtL\n9bJuUy4dVuVxd+v7wlnqL3GZvG+//dYudTd8+HAREZk3b559X+nSpcXf39/+/Y033rCfm7zk4OTJ\nk+1i3vXq1ZPAwEDx9PSUgIAAETFFx52lCCtWrCj169cXb29vcTgcsn79ehER+emnn+xi68WLF5fA\nwED7OF27drWPXbFiRbv+apEiRaRy5cr2dm+++WamX5shQ4bIfffdJ4GBgVK/fn3x9PQUy7KkSJEi\ncvLkSbttzmO4ubmJv7+/XfTbxcVFfvvtN/s8nfVyM1MA/tNPP7XLQJ47d05ERH755Re7ROHevXul\nY8eOUrJkSbEsS9zd3e3C5Fu3brXb53A4JDw8XERE1q9fb7/m7u7uUrt2bSlVqlSSkpS361Z/B1n9\nG9GAr/KMVN8XcOdvmeQMMO7u7tK4cWOpV6+eHVzd3NwkMjJSRMQOUA888IBdRLtr1672dufPnxeR\nlAF/yJAhSeq6iohcu3ZNIiIiRETk119/FcuypGjRonL16lV7mx07dsjp06dFROThhx+2PxCchcJf\nf/11+1i7du0SkZsBv2jRonL8+HFJSEiQRo0aiWVZEhQUlOnXZt++fUnatH//fvuYs2bNEpGkAb9L\nly522533hYaGisjtF4CPiYkRHx8fsSxLJk+eLCIiI0eOFMuyJDg42H5enz59UnxwJ25f4oDfvHlz\nsSxLvLy87A8kkbQLk2dGTgZ8HdJReVtuhPzbdOPGDTZt2sSuXbsoVaoU7dq1Izw8nPvvv58///yT\nI0eOANCxY0e7VmvXrl0BM/Szffv2VPf7xBNPYFkWs2bNokyZMjRr1ow333yTwoULA6awiL+/Pxcv\nXsTHx4f69evTo0cPoqKi8Pb2BrCHUNq0aWNfeOzWrZt9jF9//TXJMR955BHKlCmDZVl2ab8///wz\n06/J+fPnad++PcWLF8fhcFC1alX7sRMnTqTYvnv37sDNilKJj+twOJg3bx5Vq1bF3d0dh8NBWFhY\nmvtycnNzo1+/fgB88skngCm2AtC7d297O0nn/z7x484KYB07dkxS+rBevXq33Edu04u2SmWTSpUq\ncfDgwXS3Sxw4nD/fKti0bt2aLVu28NVXX7F9+3a2bt3K2rVrmTFjBlFRUfj6+rJ582bmzZtHZGQk\nUVFRfPbZZ8yfP58TJ07w0ksv3fLYqblVgfCMunz5Mo8++ih//fUXnp6edslCZ73Z5GmTLcuyj5u4\nGLrzuM7i6mBe69KlS3PkyBGOHTuWanH1xAYOHMh7773H9u3bmTVrFgcPHky1klVmZfY1yW3aw1fq\nDvDx8aFChQoALFq0iJiYGESEzz77DDC9UGcd1uR27NiBt7c3//73v1myZIndW798+TKbNm3i0qVL\nREVFMWTIEObOncuvv/5qF+x2FjN/4IEHAHNx98KFCwD2xWTLsjKd8z2tKYyJ7d27l7/++guA2bNn\ns2nTJsaPH5+p4ySWleLqFSpU4IknnkBE7A/A9u3bJ5lm6SxMnlZR8sSzd5yFyRcvXsz+/fvt+3fs\n2JGJM7rzNOArdYe88847gBk+qVixIv7+/nzxxRdYlsUrr7ySpFed2Jdffkn58uWpUKECDRs2pHbt\n2oDpBdesWZNTp07RtGlTSpQoQZ06dQgMDGTVqlXAzWLmY8aMwdXVlSNHjuDn50e1atV4//33AXjm\nmWeoWbPmbZ3Trea6V65cmYIFCwLw7LPPUqdOHTp27Jjm9un1lrNaXP2FF14ASLcw+Z9//km1atUI\nDg7m+vXrqbZv7NixuLm5ceHCBWrVqkXt2rUpXbq0XXc4q74/8H227Cc5DfhKZVFqxbxT0717d775\n5huaNm3KlStXOHXqFHXr1uXDDz+0PwwS79MpJCSExx9/HMuyiIqKIiEhgaZNm7Jw4UKqVq2Kt7c3\nffv2pWzZsvzxxx8cOXKEKlWq8PrrrzNq1CgAmjVrxk8//UTr1q0B+P3336lWrRpjx461pyOmdS63\nKkzu/PBJTbFixViwYAE1atRARPDw8GDp0qVp7jO9+7JSXB2gVatWBAQEAFCmTJkURcSfffZZOnfu\nTLFixdi/fz+RkZFJhooS7zM4OJgNGzbw5JNPUrRoUfbv34+np2e2re4dHT6av67/lS37SkwrXqk8\nQ98Xd4eEhAS8vb25ceMGu3fvpmLFirndpAxr37493377LSNHjrSvB+Q1WsRcKZVn7NixgwsXLvDe\ne+/dNcH+nXfeYe3ataxatQpPT0+GDh2a203KFRrwlVKZUq9evXRnxeQ1P/zwA2vXrqVKlSqMHz+e\nsmXL5naTcoUO6ag8Q98XSv39d/DZZ/CPf0Aq11Oy8jeSqYu2Xbp0sXNKPP300/b9sbGxjBkzBn9/\nf9zc3PD19WX48OG3lXtDKaXyvfnz4dKlbN9thod0PvnkE7766iv798RXrJ999lnCwsJwcXEhICCA\nAwcOMGnSJLZt28bq1asznaZUKaXytSVLcmS3GerhHzhwgKFDh9KkSRN8fX2TPLZlyxZ7efOkSZOI\nioqyly2Hh4fbWQOVUkrlrnQDflxcHN27d8fV1ZVPP/0UhyPpU1asWAGYHn/nzp0BaNu2Le7u7oBZ\n2aeUUir3pTukM2bMGCIjI5k/fz6VKlVK8bgzIRSY5eNgkhx5e3tz7NixJI8nNnr0aPvnkJAQQkJC\nMtdydc/x8vLS4T+V73l5edk/r1mzhjVr1mTfzm+VSnPTpk3i4uIivXr1su9zpk99+umnRURk4MCB\ndurQhIQEe7ty5cqJZVny2GOPpdhvOodVSqlM27Fjh9SqVSvHj3P0r6Py2vevSYn/lJDOX3SWDX9s\nSH3DhASR8HCRjh1FSpQQefVVkd9/z9Kxsxo7bzmks2vXLhISEliwYAGFChWiUKFCdo990aJFFC5c\nOMl81lOnTgFmJd7Zs2cB7IRRSil1N9txage9F/em9ke1uRp7lcj+kXzV5SuCywcn3fDGDfj0U2jU\nCPr1g5Yt4fBheP99yOV4eMuA7/x6HRMTw7Vr17h27Zo9BzQ+Pp6rV6/Srl07wCQWcl6sXbZsGTEx\nMYDJv62UUncjEWFl9EpazWtFm0/bEFgikOih0Ux+bDL+Xv5JNz5zBt59F/z84JNPYMwY+O03GDwY\nChXKnROwNqJaAAAgAElEQVRILrNfCZIP6YiIdOvWzS5HFhgYKAUKFBDLsqRZs2ap7uM2DquUUreU\nnUM612Ovy6wts6TmtJpS+8Pa8snWT+R67PXUN46KEhkwQKRYMZG+fUW2b8+WNqQmq7Ez06kVUstE\nN2fOHAICApg7dy6HDh3Cx8eHp556irFjx2bTx5JSSuW8s1fPEvprKFM3TaVuqbpMeHQCLf1bppxM\nIAKrVsGECbBtGzz/vOnNlyqVOw3PIE2toJS6J+zcuZNu3bqxc+fOTD83+lw0EzZOYP7O+XQI7MDL\nQS9Tu1QqqZ+vXYN582DiRChQAF56Cbp2BQ+PbDiD9Gm2TKWUug0iwoYjG/gg4gPW/bGOAQ0HEDU4\nijKFy6Tc+PhxmDYNZsyAoCDzc0hIilw3eZ0GfKVUvhKXEMeiPYsYFzGOM1fPMDxoOJ92/JSCbgVT\nbrx5sxm2Wb4cuneH9evh7yIqOSUhATZtgsaNs3/fOqSjlLqrpbVYL3mMuXzjMrO3zmbixomUKVyG\nEcEjaF+tPS4Ol6RPjI+Hb74xwza//w4vvgjPPQeJFkTlhLNn4X//g48/NiNEa9dCopK7gA7pKKXU\nLR27eIwpkVOYuWUmzf2aM7/zfIJ8g1JuePEizJ4Nkyebi6/Dh0OnTuCac2FSBCIi4KOPYOlSaN/e\nBP3g4JwZLdKAr5S6q6XV491xagfjIsaxdO9SetbtSWT/yJRz5wEOHTJBfu5caNXKpCYOSuUDIRtd\nvGjWZoWGwvXrMGiQ+UJRokSOHlYDvlLq3iEifHfgO8ZFjCPqdBQvPvAiEx+diJenV/INYd06Mz7/\n889myGbbNihfPkfbt2WLCfILFpgFuBMnQvPmd+7arwZ8pdRdLyYuhvk75zN+43gsLEY0GcEztZ7B\nzcUt6YY3bsCXX5pIe/GimVY5d26OroS9ehW++MIM25w6BQMGQFQUlEllMlBO04u2Sqm7lnOh1LRN\n06hTqg6vBL+S+kKpM2fM1dAPP4Tq1c34/GOPgSNTRf8yZfduc8iwMGjSBAYONId0cUn/uWnRi7ZK\nqXwn+lw0EzdOJGxnGB0CO/Bdj+9SXyi1Z4/pzX/5JXTsCCtWQJ06OdaumBj4+mvTm9+/3+RO27IF\nKlbMsUNmigZ8pdRdwblQalzEOH7+/WcGNhqY+kKpXEh7cOAATJ9ucqbVrg1Dh8KTT5rFuHmJBnyl\nVJ6WeKHU6auneTnoZeZ1nJdyodTVq2bqizPtwfDhsHhxjqU9iIszUylDQ00vvndvcx24atUcOVy2\n0ICvlMqTki+UGtl0JE9WezLlQiln2oPp0810yqlTc3Tqy5EjMHOmuVWqZL5AfPPNHUunkyUa8JVS\neUqGF0o50x4sW2bSHmzYkGNpDxISzChRaKiZxdmtG6xcaYZv7iYa8JVSeUKGFkrFx8OSJSbQHz5s\n0h5MmZJjaQ9OnTLj8tOnm0M8/7wZNcor9UwySwO+UirXZHih1MWLMGuWWRFbunSOpj0QgfBw05v/\n7jvo3NnMo7///mw/1B2nAV8pdcdleKGUM+3BnDnQujV89lmOpT04f94cJjTUzJV//nnzc7FiOXK4\nXKEBXyl1x5y9epaPN3/M1Mip1ClVh/Gtx6dcKJU87cGzz5rplTlQAFwEIiNNYF+0CB5/3KS8f/DB\nuy7VfYZowFdK5bjkFaVSXSiVPO3BsGE5lvbg0iWTIy001Pw8cCC8/z6ULJnth8pTNOArpXJEhitK\nnTljropOm2bSHoweDW3b5kjag+3bTZD/4gtTsOr996FFixzNsJCnaMBXSmWr1BZKpVpRKioKJk3K\n8bQH166Z7JShoWYOff/+sHMnlCuX7YfK8zTgK6WyhXOh1ISNEyhbuGzqC6VEzNSXiRNzPO3B3r0m\nedncuWaGzeuvmy8OOVjPJM/Lx6eulMoOiRdKhVQKYX6n+QSXD066UfK0By+9lCNpD27cMLsNDTVf\nIJ591tSH9fPL1sPctTTgK6Vuy/aT2xkXMY5v931Ljzo9+KXfL1QuXjnpRs60BzNm5Gjag8OHzSFm\nzzaXAQYNgg4dwM0t3afmKxrwlVIZlnih1O4/dzO08VAmtZmUcqHU5s2mN+9Me7B+fbanPYiPh+XL\nTW/+l1+gZ0/46ScIDMzWw9xTtACKUipdzoVS4yLG4bAcvBL8Cl1rd026UCo+3mQRmzjxZtqDfv2y\nPe3BiRMmcdmMGebC66BB0KULeHpm62HyJC2AopTKMckrSk14dELKhVJ3IO1BQgKsXm168z/+aAL8\nkiVQr162HSJfyNDs05kzZ9KoUSOKFy9OgQIFKFmyJM2aNeOLL76wt4mNjWXMmDH4+/vj5uaGr68v\nw4cP5/LlyznWeKVUzog+F82Q5UOoMqUK0eej+a7Hd6zssZJWlVvdDPYHD5rgXqmSGVP57DOIiDDR\nOJuC/Zkz8MEHUK0ajBhhCn//8YeZfaPB/jZIBvTt21dKly4t9evXl7p160qBAgXEsiyxLEs2bdok\nIiI9evQQy7LE1dVVqlevLm5ubmJZloSEhEhCQkKS/WXwsEqpOyghIUHW/b5OOn7eUbzf95Z//vhP\nOX7xePKNRH7+WaRjR5ESJURGjhT5449sbofI2rUi3buLFCsm0quXSESEuT+/y2rszNCzr1+/nuT3\nmTNnimVZ4nA45IcffpDNmzfbHwDTpk0TEZGlS5fa93399dfZ2milVPaJi4+TBbsXSOMZjaXypMoy\n9ZepcjnmctKNYmJE5s0TadBAJCBAZNo0kUuXsrUdFy6ITJkiUrOmSLVqIhMmiJw9m62HuOvdkYAv\nIhIeHi6NGzeWWrVqSYECBaRkyZLy3//+V0RExo4da38AnDx5UkRE4uPjxcPDQyzLkgEDBmRro5VS\nWXcp5pJM2jhJ/Cb6SZNZTeTrqK8lLj4u6UanT4uMHStStqxIixYi334rEh+fre349VeRfv1Mb75L\nF5HVq7U3n5asxs4MD7SdP3+eyMhI+yrxmTNn2Lx5M7GxsRw5csTezsfHBwCHw4G3tzfHjh1L8rjT\n6NGj7Z9DQkIICQm5nREppVQmHb90nCmRU5ixeUbaFaWiosxsmwULciTtwZUr8Pnn5iLs6dMmeVkO\n1xm/K61Zs4Y1a9Zk3w4z+wlx8uRJefHFF+3hmhkzZsigQYPsHn7i8fpy5cqJZVny2GOPZeunlFIq\n87af3C69FvUSr/e85MXlL8qBcweSbpCQILJihUjr1iKlSomMGiXy9zf27LJrl8iQISLFi4u0by+y\nfLlIXFz6z1NGVmNnpi+llypVinfeeYepU6cCsHXrVsqXL28/furUKUqXLk1CQgJnz54FoEIO5LFW\nSqVPklWUGnL/ECYOTVZRKnnag+HDs7Uqd0wMLFxoevPR0WZq/tatOZLeXqUj3WmZ165dY8aMGVy/\nft2+b8mSJfbPPj4+tGnTBjBvroULFwKwbNkyYmJiAOzHlVJ3RkxcDLO3zqb2R7UZ+f1IetbpyaFh\nh3jjoTduBvvjx+HNN820ym+/NSkQtm2DPn2yJdhHR8PIkVC+PPzvfyZ9zu+/w7/+pcE+16T3FeD8\n+fNiWZZ4eHhIjRo1pHLlyvZwTokSJeSPv6dkdevWTSzLEhcXFwkMDLSnbjZr1izbv5YopVJ35soZ\nGRs+Vkp/UFpaz2stq6JXpZgWLb/+KtKjh4iXlxlf2bcv245/44bIwoUirVqJlCwpMmKEyP792bb7\nfC+rsTPdZ1+/fl169uwpVatWlUKFCom7u7tUqlRJ+vbtKwcO3BwDjI2NlVGjRomfn5+4u7tLuXLl\nZNiwYXIplalbGvCVyl77z+6XwcsGS7H3ikmfxX1kx8kdSTeIixP5+muRhx4SKV9e5P33Rc6dy7bj\n//67yP/9n0iZMiIPPigSFiZy7Vq27V79LauxU3PpKHWXkr8rSo2LGMfPv//MwEYDGXL/kKQVpS5e\nNCkkJ082U2Beegk6d86WlbDx8Sa1fWioKUHbrZvJa1OrVpZ3rdKguXSUymeSV5QaHjSceR3nJa0o\ndfAgTJliqn+0bm0KuAYFpb3TTDh50nyGTJ9uasAOGmSyKhQsmP5zVe7SgK/UXcJZUWrixomUKVwm\nZUUpEdPVnjABfv4ZnnvOXIRNNIvudonAmjWmN79qFTz1FHz1FTRqlOVdqztIA75SeVziilKpLpS6\nccPUhZ040QzhvPQSzJuXLV3uc+dgzhwT6N3cTG9++nQoWjTLu1a5QAO+UnnUjlM7GBcxjqV7l9Kz\nbk8i+0fi7+V/c4MzZ0z0nTbNlHkaPdoUbXVkKAlumkRg40YT5L/5Bp54wgzhNGmS7YWq1B2mAV+p\nPEREWHVglakodXo3Lz7wIhMfTbZQKioKJk0yvfpsTHtw8SKEhZlAf/WqSXcwbhx4e2d51yqP0ICv\nVB7grCg1fuN4LKyUFaVEzJSYiRPNuPzzz2db8plt20yQ/+ILaNHCBPlHHsnyFwWVB2nAVyoXJa8o\nNb71+KQVpa5dM+PxidMeLF6c5ZWwV6+aLwihoWbBbf/+sHs3lC2bDSel8iwN+Erlguhz0UzcOJH5\nO+fzZOCTfNfjO2qXqn1zg+PHzdj8jBlmOuW0aRASkuVB9D17TLWoTz81u33rLXjsMXBxydr5qLuD\nBnyl7pDEC6XW/rGWAQ0HsHvw7qQLpTZvNr35ZcvMSqb16yEgIEvHjYmBRYtMb/6338xszV9/NSl0\nVP6iK22VymHJF0q9HPQyfer1ublQKj7eTIeZOBEOH4YXXzQpJb28brnf9Bw6ZCbxzJ5tVr8OGgRP\nPmmmV6q7k660VSqPci6UmrBxAmULl025UOriRZg1y6Q9KFPGzJ/v1ClLaQ/i4syXg9BQ04vv1cus\nwapWLZtOSt3VNOArlc0SL5QKqRTC/E7zCS4ffHMDZ9qDOXNM2oPPP4fGjbN2zGMwc6YZ8q9Y0fTm\nv/4aPD2zeDLqnqIBX6lssv3kdsZvHM+SvUvoWSfZQqnU0h5s356ltAcJCfDDD6Y3v2YNPPOM6d3X\nrZs956PuPTqGr1QWSKKKUrv/3M3QxkMZ2HDgzYVSqaU96NULChW67WOePg2ffGJm2xQpYnrz3bpB\n4cLZdFIqz9IxfKVygXOh1LiIcTgsR8qFUmfOmIj84YfZkvbA+QXho4/MwtqOHU2Gyvvv13QHKuM0\n4CuVCYkXStUuVZvxj46nlX+rmwuloqJMb37BgmxJe3Dhgll3FRpqhnAGDTJT8rM4gUflUxrwlcqA\nA+cOMGHjBMJ2htEhsEPShVIisHKlGZ/fvt1E5SykPRAxM2xCQ82F1zZtzBeFhx/W3rzKGg34St3C\nhiMb+GDDB/ZCqajBUTcXSqWW9uCbb2477cHly2aYJjQUzp83ycv27gUfn2w8IZWv6UVbpZKJT4hn\n0W9/L5S6YipKJVkolTztwfDhWUp7sHOnCfKffQbNmpkvCK1aafIylZJetFUqmySvKPVqk1eTLpTa\nvNkM2yxbBt27ZyntwfXrpmJUaKhZEdu/P+zYAb6+2XhCSiWjPXyV7x2/dJwpkVOYsXkGIZVCGNFk\nxM2KUtmc9mD/fjN5Z84caNjQDNs88US21BRX+YD28JW6TYkrSvWo0yPpQqnEaQ9Kl76Z9qBAgZsz\ncpJJ6w8xNtZ8ZoSGmuGbvn1NRanKlXPqzJRKnQZ8la+kW1Hq0CET5J1pDz77zIzT34Y//jDD/LNm\nmZGfgQOhc2dwd8/GE1IqEzTgq3whJi6GsJ1hjI8Yn3KhlAisXZvhtAfOnvxzzz1HkyZNeO655+zH\n4uPNDM3QUNiwAXr0MOkPatS4I6ep1C1pwFf3NOdCqambplK3VF0mPDrhZkWpGzfgs09vpj0YNgzm\nzr2ttAcnTpg0xNOnm+n3zz9vSgbed18OnJRSt0kDvsoRmR3nzm7R56KZsHEC83fOp0NgB1b1WHVz\noVTitAeBgbed9kDEIiqqLE8/bXrxXbqYQiMNGmT/+SiVHdJ9h48bN45HHnmEcuXK4e7ujq+vL126\ndGHXrl32NrGxsYwZMwZ/f3/c3Nzw9fVl+PDhXL58OUcbr1RiIsL6P9bT6YtOBM8KpphHMaIGR/HJ\nk5+YYB8VZQbSAwLgwAFYvhx+/BHatctUsD971hT6XrhwLF980Zjmzc0Eno8/1mCv8jhJR8WKFcXh\ncEjVqlUlMDBQLMsSy7KkUKFCcvjwYRER6dGjh1iWJa6urlK9enVxc3MTy7IkJCREEhISUuwzA4dV\n94itW7dK3bp1c/QYsfGx8uWuL6XxjMbiP8lfpvwyRS7HXDYPJiSIrFwp8uijIqVKiYweLXLyZKaP\nkZAgsm6dSI8eIkWLivTsKfL44+/IjBkzs/lslEpbVmNnut2afv36ER0dzd69e9mzZw/jxo0D4MqV\nKyxatIgtW7YQFhYGwKRJk4iKimLhwoUAhIeHs3jx4hz7sFL52+Ubl5n8y2SqTqnKxF8mMrLpSPYN\n2ceQB4ZQMM4yA+o1a8LIkSZZ/OHDMGpUpnLc/PWXGfmpW9dMp6xf33w5mDsXSpU6oLlt1F0l3TH8\nt956K8nvLVu2tH/28PBgxYoVgBmz7dy5MwBt27bF3d2dmJgYVq5cSceOHbOzzSqfS1xRqrlfc+Z3\nnn9zodSxYyZCO9MeTJ0KzZtnOu3Bli1mps2CBSbNwYQJ8MgjmrxM3d0yfdF2/PjxAHh7e/PUU08l\n+UDw+TvLk8PhwNvbm2PHjnHkyJFU9zN69Gj755CQEEJCQjLbFJXP3HKhVDakPbhyxcysCQ2FP/+E\nAQNgzx6z7kqp3LBmzRrWrFmTbfvLcMC/ceMG/fr149NPP6Vo0aIsXrwYb2/vNLeXdGZjJA74SqVF\nElWUijodlXShVHy8yR88YQL8/rtJezBlSqbTHuzebS64hoVBkyZm1KdNG3BxyaGTUiqDkneGx4wZ\nk6X9ZSjgnzlzho4dO7J+/XrKli3LsmXLqPt34cwKFSrY2506dYrSpUuTkJDA2bNnUzyuVEY5K0qN\n3zgeC4sRTUbwTK1nzEKpixfh44kp0x5kIiFNTIz5rPjoI4iONmuttm4Ffbuqe1m6F2337NlD48aN\nWb9+PfXr1ycyMtIO9gBt2rQBTE/MebF22bJlxMTEJHlcqYw4d+0c7659F79Jfnyx+wvGtx7P9kHb\n6VW3F26/HzWpiP38TDKazz4zy1m7dMlwsD9wAF57zQT2WbPMWqvff4d//1uDvcoH0pvGU61aNXsq\nZu3ataVx48b2beZMMyWtW7duYlmWuLi4SGBgoBQoUEAsy5JmzZrlyNQidffI6LTM/Wf3ywvLXhCv\n97yk7+K+svPUTvNAQoJIeLhIx44iJUqIjBwp8scfmWpDbKzI11+LtG4t4u0tMmKEyL59t3M25r17\nq5tSOSmr77F0u0UxMTH2qsndu3cneaxt27YAzJkzh4CAAObOncuhQ4fw8fHhqaeeYuzYsdn2waTu\nPSLChiMbGBcxzq4otXvwblNR6sYN+DRraQ+OHoWZM82tUiVTWCQLBamUuutpPnyVo7Zt20afPn3Y\ntm2bfV9cQhyL9vxdUerqaV4OevlmRankaQ+GD89U2oOEBFi1ysy0+fln6NbNLK6tXTt7zietlBFO\n+r5WOUnz4au7RvKKUiObjrxZUSoqCiZNgi+/NBdgV6yAOnUyvO8//7yZvMzLyyQv+/TT28qDptQ9\nSwO+ynYpesGFwWplQQPgMGz4YAPB5YNNWuJVq8y0ym3bTJTORNVuEdOL/+gj+O47k2v+iy/g/vuz\n/ZSUuifokI7KdnbALwUEA9WAHcBG4DzI1aswb54Zny9QwAzbdO2a4cog58+b4fzQUDPSM2gQ9OwJ\nxYrlzPkkpkM6KjfpkI7KU0QEqmACvQ/wC/AdcA3KAIMBKlaExo0zlfZABCIjTZBftMgM60+fDg8+\nqOkOlMooDfgqWzgXSo2LGAetgA3ALiDejOQMB9oCYZCptAeXL5sVsKGhcOmSuQD7/vtQsmQOnYhS\n9zAd0lFZ4qwoNW3TNOqUqsMrwa/QukprHMCTwEtARWAKMAu4QMaGPXbsMEH+88+hWTMzvN+yZaZr\nlOSo1EocKpWTdEhH5Yroc9FM3DjxZkWpnquo5VMLLl5kGDAUOAlMBL4G4jOwz2vXTHbK0FA4cgT6\n94edO6FcuZw8E6XyDw34KsNuuVDq0CFz8XXuXIKArkBkBve7d68Zj58718ywef11M0afidQ4SqkM\nyENfkFVeFZ8Qz1dRXxE8K5jei3vTwq8Fh4cd5p3mYymzLdrMh7z/fnBzg61bMxTsb9wwvfkWLcyQ\njZsb/PKLqTrYvr0Ge6Vygv5ZqTQlXyj1WtPXaF+tPS5x8WaBlDPtwUsvwZw59iqnxGOMyVfaHj4M\nb75pFkkFBpoplR07moCvlMpZGvBVCscvHWdK5BRmbJ6RtKLUmTPw/94zaQ+qV4cxY+Cxx9K9kiri\nYOlSMzb/yy/Qowf89JMJ+EqpO0cDvrJtP7md8RvHp6woFRVl5kNmMu3BiRMwY0Yp9uxZxrvvmt78\nV1+Bp+cdOBmlVAoa8PM5SauilEexv9MeDM5U2oOEBPjxR9ObX70aWrQogJ/fC0REfHGHzkgplRYN\n+PlUTFwMYTvDGB8xHofluFlR6kZ8yrQHixenm1P49Gn43/9MosuCBc3nwyefwMGDR+nTZ++dOSml\n1C3li4CfVv6T/Lj4y7lQauqmqdQtVZcJj06gpX9LrBMn4O0xMGMGBAVlKO2BiFk0+9FHpnZ4hw4m\nQ2XjxpruQKm8KF8EfGUWSk3YOOHmQqkeq6hdqjZs3mwyjy1bBt27ZyjtwV9/mS8BoaEQF2fG5qdM\ngeLF79DJKKVuS76Yhy8iiAh79+4lICDA/v1eJyKs/2M9Hb/oSPCsYIp5FCNqcBSftJtJ7fX74eGH\nzZzIOnXg4EHTq79FsP/1V+jXz+Q+W7vWBPk9e8yszPwU7C3LwrIsZs+eTb9+/ezflcrrtId/D0pe\nUWp40HA+7fgpBa/Hw6zZMHkylCplxuc7dbrlKqcrV0yt8NBQMytz4ED47TcoXfoOnpBSKltowL+H\npFlR6vDv8NpbJndBq1Ywf74Zp7+FXbvMBdj5800K4n//G1q3BheXO3QyeVh++Hao7k0a8O8Bxy4e\nY0rkFGZumXlzoVS5xmbc5c2nTVmo556DrVuhQoU093P9OixcaHrzBw+ap2zbBuXLZ75NyYc4nL9r\nsFQq92jAv4vtOLWDcRHjki6UKuhrFkhNeMEkkB82zPTsb1HcNTra9ObnzIH69eHll6FdOzMrUyl1\n79CAf5dJc6HUlXj48GOT9iAw0KQ9aNs2zbQHsbGwZInpzW/fDn36wIYNUKVK9rVTKZW3aMC/hYzO\n378T8/ydFaXGbxyPhXVzodS+A/DS66ZX37FjumkP/vjDTLWfNcsE90GDTLLLDJaTVUrdxe75gJ9a\nME58X17viSavKDW+9Xha+rXA+v57eK39zbQHv/1mZt6kIj4evvvO9ObXr4du3eD776FmzTt8Mkqp\nXHXPB/yscH4YvPbaaxQvXpzXXnvtlttduXIFHx8frly5kuVjHzh3IOlCqZ6rqFXI3yxlnVgrQ2kP\nTp40aYinTzc1YJ9/3kyxLFgwy81TSt2FNODnIWlWlLqYAJM+zFDaAxFYs8b05letgqeeMhkqGzW6\n8+ejlMpb0l1p+/PPP9OuXTtKlSqFw+HA4XAwZsyYJNvExsYyZswY/P39cXNzw9fXl+HDh3P58uUc\na/i9JC4hjgW7FxA8K5hei3vxiN8jpqJU0U6Uef5VqFXL5DNYt85caX3kkRTB/tw5mDDBXK8dOtQs\noj182HxGaLBXSkEGevhbtmzhu+++o2rVqpw+fRpIOS7+7LPPEhYWhouLCwEBARw4cIBJkyaxbds2\nVq9ercvO03Ap5hJFmhWBIOASEAGOPfAjLzLkoS9NxH7xRZPDwMsrxfNFYONG05v/5ht44gkzhNOk\niSYvU0qllG4Pv1evXly6dInIyNSrlG7ZsoWwsDAAJk2aRFRUFAsXLgQgPDycxYsXZ2Nz7w3HLh7j\n9R9ex2+SH1QEFkLh2fDSHogGRgIMGWJWP736aopgf/GiyVBZrx706mUm5URHm4RmTZtqsFdKpS7d\ngF+8eHE8PDzSnM2yYsUKwPT6O3fuDEDbtm1x/3ue38qVK7OrrXe97Se302tRL2p/VJursVeJ7B+J\n/L8DyFMv8ed999EjIAC/iAiCRaBLlxQ5brZtM9MoK1Y0RUbGjTM1SV55Bby9c+mklFJ3jSxftD1y\n5Ij9s8/f1ZAcDgfe3t4cO3YsyeN3k9SGoV5//XX7Z+cHYHrTPhMSEuyFUrv/3M3QxkOZ9OhEvH7d\nBc+NgPBweO45Phs5ki1nztAwWY6bq1fNFPvQUDh+HPr3h927oWzZ7DpTpVR+kWOzdNKb3z569Gj7\n55CQEEJCQnK8Hfv27aNdu3bs27cvR46VhAtQB2p/VBuH5eCV4FfoWq0zbgsXw7BWKdIeXJ4yxaSj\n/NuePSbdgbOgyJtvmoWzmrxMqfxjzZo1rFmzJtv2l+GAn9aF1wqJknGdOnWK0qVLk5CQwNmzZ1M8\nnljigH9P8QQaAQ8ApzAVpYrUw5oxAx6vBtWrp5n2ID7ehS++ML35PXtM8rJff4VKlXLhPJRSuS55\nZzj5DMnMynABlMQ95cQ/t2nTxr7PebF22bJlxMTEJHn8nlccaAsM/fvneVD9U2j1nwVYVauaq6or\nVsAPP5jMZImC/aFDsHRpMHPm/Ivp02HwYJMC4Z13NNgrpbKRpGPhwoVSuXJl8ff3F8uyxLIsKV68\nuFSuXFm6d+8uIiLdunUTy7LExcVFAgMDpUCBAmJZljRr1izVfWbgsDli7969EhAQkKFtgVveREQS\nEhKE8gj/QHgV4RGEgsijICtBToDI6NEiJ0+m2H9srAi0F1gucFrgA4GqKY6hlFJOWY0L6Q7pXLp0\nicwKy1sAAA5rSURBVIMHDyYp43bhwgUuXLhgD9fMmTOHgIAA5s6dy6FDh/Dx8eGpp55i7Nix2f4B\nlSdY8FXUV4yLGAcdgQjwXAg942AYEAdMAJ4Ero8aleSpx46ZxGUzZgC8BoQCnYDr2dc8LdqulEqF\nJbkQBSzLypXgk5mLtqkGTTegPhAETeo04ZXgVxhSozMvAP2BjcBE4KdETxEREhLMSE5oqEl78Mwz\nplRgvXq3njB/u6+RBnyl7k1ZjZ2aSycjCmMuwjYEDgELYX3PyfDuREKAMKApZtFUUt68/76ZbVOk\niEleNmcOFC6cs811viGOHj1KUFAQR48ezdkDKqXuChm+aHs3cw5HVatWjf379ycZnkqLiNgLpTxe\n9qBJSBMOvLoPeWYh4veQKf5dpw7Fz53jRRH2i3D58mU8Pe/j55+Fbt2EYsVO89tvJkPlli0wYEDO\nB3ullEqL9vCTkVQqSr14oTePHjiG/+xHTc754cNNwE+0EvbCBZg1y5Xr1zcxYIBZETt1aqopcJRS\nKlfki4CfkTGvmLgYwnaGMT5i/M2FUoWCcJsWypWP5/BH1aowf75JT2zv18yTDw2Fr7+Gli1dcHMb\nTlTUd5rPRimV5+SLgH8rzopSUzdNpW6pukxoPZ6Wxz2w/jURfn4FnnuOyX364KhUiep/B/vLl80w\nTWgonD9vLsDu3QsFC8bg47MuV4P93V7hSymVc/LFGH5qos9F88LyF6gypQrR56P5vssyVtKDVs/8\nE6tfP2jZ0qQn/s9/uFCkCAA7d8ILL0CFCrB8Obz7rllP9dpr8HcaoUwREfs2efJkhgwZkuQ+pZTK\nTvmqhy/JKkr1b9CfPf9YR+mwxTD4CZP2YPToJGkPrl+H3bvrs2vXg0yebJKX7dgBvr65ey5KKZVZ\n+aKHn1pFqd8fXcm7X56ldP0H4cCBFGkP9u0DyxqHp+dpli0rxu+/v8jx466MGWOlCPbOWT+FChXi\n6tWrGZoFpJRSd9o93cO/fOMys7fOZuLGiZQpXIaRTV7lyd/vw+XtSbBtrJkY/9tvZuYNEBtrKkeF\nhprhG4jFlKM6mItnoZRS2eOeDPjHLh5jSuQUZm6ZSUilED57bBaN1+yHp/4PChQw0yoXLwYPDwB+\n/92kOpg9G6pWNRdhO3UCd/fXgddvfTD0QqhS6u5wTw3p7Di1g96Le9sVpTa3W8pX26vR+MF/wLff\nwrRppmxUnz7EF/Dg22/NCE6DBmbmzY8/mtQHXbvC3wW7lFLqnnHX59JJbaHUYKsxRUJnw7Jl0L07\nDB0KAQEAnDhhevLTp0Pp0mZUp0sXuO++bGlOpmV33puMrCBWSt2dsho782TAz0gQtFwtqA0Em98d\n66H9LljU9CEznfLFF6FfP/DyIiEBfvrJjM3/8IMJ8AMHmp59bsvJRGeaS0epe0u+S5527to5Qn8N\nhZeAU1B4BTx72NQdOQUwZIid9uDsWfjfOJO8zMPD9OZnzTKJzPIK7XErpe6UPDmG71x4FBsbi4uL\nCyJC9NlohiwfQpXJVdh/bj97uixHvF/izNnCdC5bFv+ICIJFkKe7sP4XV3r2hMqVYft2+N//zL/P\nP5+3gr1SSt1JeTLgJya+QucvOxM0K4ii7kXYV282n8z5i8C2PcHNjZ8mTOCDRo24WCOIDz+EunWh\nb1+oX99Mr587F5o0QXPbKKXyvTw3hm9ZlvkYCsSMzxeEAhugy1Z4KR4aBQQweP9+5gBXAVORZBDw\nNPAD8BEJCas1wKNj+Erda+6pi7aXb1ym8MOFzVqnS1B8LQzcDy8AezDVpL6Nj8dyKQQ8gwn0PsB0\nYDZ/j+Ln+3FxrXil1L3pnrhom3ihFBWg+ucw7BR0ARYBjwE7AajOsOEO4AgQAYwBVgIJudRypZS6\ne+TqGL6zolTtj2pz9cYVdvr/lxULYPUpOI4Z1XkON3byDLAG+OHvi64NgCeA5WiwTylxxk3NvqmU\ncsq1IZ2Wc1sSdTqK4XUG8vxvRSj44XRwdaXvzp18BsTgDwwA+gA7gFBgCSKxurhIKZUv3bVDOgNK\nt6PTnoa4dJtqqkhNnUrcQ835n1snzNh8A2AO8CCplQdXSimVObl30bZYMZP2YNgwjnoGMHMmzJwJ\nx46tw/TmvwJiUjxXRLSHr5TKl7Law8+1MfyE6IOsbDeVDq8GUKcOnDljUtLDQ0AYqQV7pZRSty/X\nevh+foKXl1n9+swzUKhQyu3i4uLw8PAgLi4uzX198803zJ49m2+++SYHW6yUUrnvrp2HHxkp3H//\n/2/vfEOa/No4/j3TOVdaNuemTU23MjVEjUItQsuKCCFIDdIM9iIsiP68qxBMoqAXud5IoAYJRglZ\nUdQvKEIEIRTJopT+WGLgfzPQctPa9Xvhs/O45/HPbHO39zwfOHCf+z54X1++4/LeznXOPfv1mXDa\nPE3UmgsEgmXGkvpJ5+7du9i8eTPUajU0Gg3y8/PR2dk549jZkr1AIBAIFgePJfybN2+ioKAAbW1t\nMBgMICLU19dj+/bt6O/vX9DfcqWOfCnXmjc0NEgdwqIi9MkbX9bny9o8gUcS/sTEBM6dm3oVYF5e\nHj5//oz29nYEBwdjYGAAV65c8cRtZIOvf+iEPnnjy/p8WZsn8EjCb2lpwfDwMAAgNzcXABAREYH0\n9HQAwLNnzzxxG4FAIBC4gUcS/rdv3wBMTSjodDp+3nHsuC4QCAQCCSEPcOfOHWKMkUKhoJcvX/Lz\nhYWFxBgjtVrtNB6AaKKJJppof9HcwSNbK0RHRwNTkThN0A4MDDhdd0BLZHJVIBAIlhMe+Uln69at\nCA0NBQDU19cDAHp6evDq1SsAwL59+zxxG4FAIBC4gccWXlVVVaG4uBgAEBMTg+HhYYyOjiIsLAxv\n3rxBeHi4J24jEAgEgr/EY3X4x44dQ21tLVJSUtDX1wc/Pz8cPHgQTU1NPNkvZGGWnLh48SIUCsWM\nzW6X1379jY2NyMnJgV6v5xrKysqcxkxOTqKsrAxGoxEBAQGIjIzE2bNnMTY2JlHUruOKvqysrBm9\n3LFjh0RRu861a9ewa9cuGAwGqFQqREZG4tChQ3j37h0fI2f/XNEnZ/+qq6uxZcsWaDQaKJVKhIWF\nITMzE3V1dXyMW/65NQOwAKqrq4kxRowxMplMFBISQowx0uv11NfX560wFoXS0lJijJFOp6OMjAyn\n9ufPH6nDWxAWi4X8/f0pMTGR+1VWVuY05siRI8QYI39/f0pISKCAgABijFFWVhbZ7XaJIncNV/Rl\nZmYSY4zWr1/v5OXx48clitp11q1bRwqFguLi4ig+Pp5rDAoKoq6uLiKSt3+u6JOzf2azmcLDwyk1\nNZWSk5NJqVRyjS0tLUTknn9eSfg2m420Wi0xxig/P5+IiHp6emjVqlXEGKNTp055I4xFw5HwzWaz\n1KG4zfDwMI2Pj9PY2NiMCbG1tZWfr6ioICKix48f83P379+XKnSXmE8f0X8TRk1NjURR/j2XLl2i\nL1++8H55eTnXabFYZO/fXPquX79ORPL2z2q1OvUdD8oKhYJevHjhtn9e2R55uSzMunfvHtRqNSIi\nIpCTk4O2tjapQ1owGo0GgYGBs1ZS/TO1hzUYY9zL/fv3Q6VSAVj6Xs6nbzpnzpyBSqWC0WhEcXEx\nrzpbypSUlCA2Npb3d+/ezY8DAwNl799c+hwaHMjRP5VKhcbGRqSnpyMpKQknTpyAVqvF1atXkZ2d\n7bZ/Xkn4vr4wizEGPz8/REREwGg0or+/H0+fPkVGRoYsk/5cTPfK4Z9CoYBWq/2/63KFMYYVK1Yg\nMjISer0eXV1dqKqqQkZGBn79+iV1eAuivLwcAKDVapGXl+dz/k3Xl5+fD0D+/o2MjKC5uRnt7e34\n/fs3hoaG0NraisnJSbf9k/Ql5q48ZcmBgoICDA4O4sOHD3j//j3/L2uz2VBRUSFxdN7BV7wEAIvF\ngpGREbx9+xbd3d04f/48AODr16948OCBxNG5xsTEBI4ePYqamhqsXr0aDx8+5ElhJuTm30z6HKXh\ncvfvwIEDsNvt6OnpwcmTJwEAdXV1qKmpcXtbeK8k/IUuzJIbGzZsQEhICO/v3bsXGo0GgPyemBzM\n9sGa7pXDS7vdzn+yk4uXc70mMyUlBUqlkvcPHz7Mj+Xg59DQELKzs1FbW4u1a9eioaEB27ZtA+Ab\n/s2lD5C/fw70ej0uX77M+69fv0ZUVBTv/41/Xkn4vr4wy2KxoLe3l/efP3+O79+/A5hakyBHpj8x\nTD92eEX/2f4aAJ48eQKbzeZ0fakzm77BwUHcuHHD6av/9JK4pe5nR0cH0tLS0NTUhNTUVDQ3NyM5\nOZlfl7t/8+mTs3/j4+OoqqqC1Wrl5x49esSPdTqd+/55cIJ5TiorK/lMcmxsLK/Q0el01Nvb660w\nFgVHqVh0dDQlJCRwncHBwdTR0SF1eAuivr6eTCYTGY1GrkOj0ZDJZKLCwkIiIiooKCDGGPn5+VF8\nfDwvHcvMzJQ2eBeYT19XVxcxxkipVFJ8fDxFRUXxcZs2bSKbzSa1hDnZuHEjjzcpKYnS0tJ4q66u\nJiJ5+zefPjn7NzIyQowxCgwMpMTERDKZTDz20NBQ6u7uJiL3/PNawiciun37NqWmppJaraY1a9ZQ\nbm4uffr0yZshLAqVlZW0Z88eMhgMpFaryWg0UlFREX38+FHq0BbMrVu3eBnY/7adO3cSEdHk5CSV\nlpZSbGwsqVQqMhgMdPr0aRodHZU4+vmZT9/Pnz+ppKSE0tLSSKvV0sqVKykxMZEuXLhAP378kDr8\neYmJiZlRm0Kh4OWncvZvPn1y9s9qtVJRURHFxcVRUFAQqVQqiomJIbPZTJ2dnXycO/5J8k5bgUAg\nEHgfSat0BAKBQOA9RMIXCASCZYJI+AKBQLBMEAlfIBAIlgki4QsEAsEyQSR8gUAgWCb8C33WyuW7\nLpiCAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1054fead0>"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now define a vector function representing the 2 equations to be solved for a and b, for the Gaussian analytic uncertainty case."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def alt_func1(a,xi,ci):\n",
      "    term2 = (ci/(a[0]+a[1]*xi))**2\n",
      "    eq1 = len(xi)-term2.sum()\n",
      "    term2 = xi*(ci/(a[0]+a[1]*xi))**2\n",
      "    eq2 = xi.sum()-term2.sum()\n",
      "    return(np.array([eq1,eq2]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And use scipy.optimize.root to find the roots of this equation. \"args\" are the extra arguments to the function, in this case our data."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "spo.root(alt_func1,np.array([1.5,1.2]),args=(xi,ci))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "  status: 1\n",
        " success: True\n",
        "     qtf: array([  7.05012347e-10,   9.96187576e-09])\n",
        "    nfev: 12\n",
        "       r: array([ -10.73741172, -101.08350504,  -11.6317085 ])\n",
        "     fun: array([ -4.27444746e-11,   1.12549969e-11])\n",
        "       x: array([ 2.53529043,  1.34635396])\n",
        " message: 'The solution converged.'\n",
        "    fjac: array([[-0.18532915, -0.9826765 ],\n",
        "       [ 0.9826765 , -0.18532915]])"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Same thing for the Poisson uncertainty case"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def alt_func2(a,xi,ci):\n",
      "    term2 = ci/(a[0]+a[1]*xi)\n",
      "    eq1 = len(xi)-term2.sum()\n",
      "    term2 = xi*(ci/(a[0]+a[1]*xi))\n",
      "    eq2 = xi.sum()-term2.sum()\n",
      "    return(np.array([eq1,eq2]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "spo.root(alt_func2,np.array([1.5,1.2]),args=(xi,ci))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "  status: 1\n",
        " success: True\n",
        "     qtf: array([  3.26631559e-08,  -1.65932458e-08])\n",
        "    nfev: 10\n",
        "       r: array([ -5.58603901, -51.17813446,  -6.35955669])\n",
        "     fun: array([  2.21742624e-11,   2.47695198e-11])\n",
        "       x: array([ 2.09937404,  1.32073883])\n",
        " message: 'The solution converged.'\n",
        "    fjac: array([[-0.21038093, -0.97761949],\n",
        "       [ 0.97761949, -0.21038093]])"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now what about the case with uncertainties in both parameters? First define some x uncertainties:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x_unc = np.log(xi)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.errorbar(xi,ci,yerr=obs_unc,xerr=x_unc,marker='s',ls='None')\n",
      "#plt.plot(xfit,1.16+1.41*xfit,ls='solid',marker='None')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 29,
       "text": [
        "<Container object of 3 artists>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHW9JREFUeJzt3X1QVNfBBvDnLh8LgomSZRVdEZasQRPH6NjxoxPYhKRD\nqdOmAWwVaaNJmqbNaPjjTdOWEXaiaZMZNfSNr63aJjhsrakok4jQpsOsdkw7WJAYhZgEYcSPEF2x\nIwbXhT3vH5u9fLiwC7vLftznN3Nn9+693HsORx/Onnv2riSEECAiooinCnYBiIhocjDwiYgUgoFP\nRKQQDHwiIoVg4BMRKQQDn4hIIbwKfKvVipdeegnp6elQq9VITk7Go48+ipaWFgCA3W6HyWSCXq9H\nbGwsdDodSkpK0NvbG9DCExGR9yRP8/CtViuWLVuG8+fPIzo6GgaDAdHR0ejs7MTbb7+Np556CsXF\nxTCbzYiKioLBYEB7ezvsdjuys7PR0NAASZImqz5ERDQKjz380tJSnD9/HjqdDm1tbTh79iw++ugj\n9PT04Nvf/jaam5thNpsBABUVFWhtbUV1dTUA4NixY6ipqQlsDYiIyCtjBr4QAu+++y4AID09HYWF\nhUhISMCDDz6I3bt3Iz4+HnV1dQAASZKQn58PAMjLy4NarQYA1NfXB7L8RETkpeixNl69ehU9PT0A\ngH/+85/QarWYOXMm2tra8LOf/QwDAwO4ePGivL9WqwUAqFQqaDQaXLp0CV1dXQEsPhEReWvMwO/v\n75efazQatLe3Iz4+Ho888gj+9a9/4a233oLRaHT7s2NdGuCYPhHRxPhy+7Mxh3SSk5MRExMDAJg3\nbx4SEhKgUqmwZMkSAEBnZyfmzJkj79/d3Q0AcDgcsFqtAIDU1NRRCx2pS1lZWdDLwPqxfkqsXyTX\nTQjf73M5ZuDHxMTg0UcfBQCcO3cOt27dgsPhwKlTpwAADzzwAHJzc+UAd12sra2thc1mAwB5OxER\nBZfHWTpbtmyBWq2G1WqFXq9HRkYGPvzwQ0iShLKyMixZsgRr1qwBAGzatAnz58+XL95mZWXhySef\nDGwNiIjIKx4Df+nSpbBYLMjJyUFfXx9u3rwJo9GIhoYGfP/73wcAVFZWYvPmzUhNTUVHRwe0Wi02\nbtyII0eOBLwCoWi06xqRgvULb5Fcv0iumz94/OBVQE4qSX4ZjyIiUhJfs5P30iEiUggGPhGRQjDw\niYgUgoFPRKQQDHwiIoVg4BMRKQQDn4hIIRj4REQKwcAnIlIIBj4RkUIw8ImIFIKBT0SkEAx8IiKF\nYOATESkEA5+ISCEY+ERECsHAJyJSCAY+EZFCMPCJiBSCgU9EpBAMfCIihWDgExEpBAOfiEghGPhE\nRArBwCciUggGPhGRQngM/PLycqhUKreLw+EAANjtdphMJuj1esTGxkKn06GkpAS9vb0BrwAREXkn\n2tsdk5OTkZGRMew1SZIAABs2bIDZbEZUVBQMBgPa29tRUVGBlpYWNDQ0yPsREfnbaPEixOSWIxx4\nPaTzne98Bx9++OGwRZIkNDc3w2w2AwAqKirQ2tqK6upqAMCxY8dQU1MTmJITEdG4eB34Bw8eRHx8\nPFJSUrBq1Sq0tLQAAOrq6gA4e/v5+fkAgLy8PKjVagBAfX29v8tMRCQTwrmcPg089NDgOt3NY+BL\nkoSoqCikpKRAr9eju7sbR48exYoVK9DS0oKuri55X61W6zyoSgWNRgMAw7YTEVHweBzDX7t2LTZt\n2oRp06YBAP7+978jNzcXNpsNO3fuRHS0+0MID39iy8vL5edGoxFGo9H7UhMRKYDFYoHFYvHb8STh\nKZnd0Gg0uH79Op544glkZ2ejtLQUkiTh0qVLmDlzJhwOBxISEmCz2fCTn/wEv//974efVJI8/kEg\nIhqPjz8G1q51PkYqX7PT45DOjh07cOXKFXn9gw8+wPXr1wEA6enpyM3NBeDs0bsu1tbW1sJmswGA\nvJ2IiILLYw8/LS0NXV1d0Ol0SEhIwCeffAIASExMRGNjIzIzM1FUVIT9+/dDpVLJ0zL7+/uRlZXl\n9u0Ie/hE5G/s4XvmsYf/61//Gjk5ORgYGEBnZyfS09Oxbt06NDU1ITMzEwBQWVmJzZs3IzU1FR0d\nHdBqtdi4cSOOHDky4YIREZF/TWgM3+eTsodPRH7GHr5nvJcOEZFCMPCJiBSCgU9EpBAMfCIihWDg\nExEpBAOfiEghGPhERArBwCciUggGPhGRQjDwiYgUgoFPRKQQDHwiIoXw+I1XREShSpIkN68NPudN\nGodjD5+IJkSShocrDReKvx/28IkoAozsyYdY0oYIBj4R+eTtt4NdgtGFctmCgYFPRD45fjzYJRhd\nKJctGPiNV0Q0Ia7x6WD+Vx68aOt+SCeYOROI34+v2ckePhFFAI7Ze4OBT0QTwjfpYwvF3w+nZRJR\n2BJCyMvp0wIPPSSGvUbDMfCJiBSCgU9EpBAMfCIihWDgExEpBAOfiEghGPhERAoxrsBfvXo1VCoV\nVCoVCgsL5dftdjtMJhP0ej1iY2Oh0+lQUlKC3t5evxeYiIgmxusPXr399ts4ePCgvD70PtQbNmyA\n2WxGVFQUDAYD2tvbUVFRgZaWFjQ0NLi9ZzUREU0ur3r47e3t2LhxI1auXAmdTjdsW3NzM8xmMwCg\noqICra2tqK6uBgAcO3YMNTU1fi4yERFNhMfA7+/vR1FREaKjo1FVVQWVaviP1NXVAXD2+PPz8wEA\neXl5UKvVAID6+np/l5mIiCbA45COyWRCY2Mj/vznPyMtLe2u7V1dXfJzrVYLAFCpVNBoNLh06dKw\n7UOVl5fLz41GI4xG4/hKTkQU4SwWCywWi9+ON2bg/+c//8FvfvMbFBcX44c//OG4DuzpPhZDA5+I\niO42sjNsMpl8Ot6YQzpnzpyBw+HAX//6VyQmJiIxMVHusR8+fBhTp07FrFmz5P27u7sBAA6HA1ar\nFQCQmprqUwGJiMg/xgx81+wam82Gvr4+9PX1yT33gYEBfPXVV1i1ahUAZ4/edbG2trYWNpsNAJCb\nmxuwwhMRkffGDPwf//jHcDgcGBgYkBdXj72goAADAwNYsmQJ1qxZAwDYtGkT5s+fL1+8zcrKwpNP\nPhngKhARkTfG/QUokiTdNa++srISBoMB+/btQ0dHB7RaLQoKCrBlyxa/FZSIiHzD77Qloojw8cfA\n2rXOx0jla3byXjpERArBwCciUgh+iTkRhbWRt+pyrXPU+G7s4RMRKQR7+EQU1tiT9x57+ERECsHA\nJyJSCAY+EZFCMPCJiBSCgU9EpBAMfCIihWDgExEpBAOfiEghGPhERArBwCciUggGPhGRQjDwiYgU\ngoFPRKQQDHwiIoVg4BMRKQQDn4hIIRj4REQKwcAnIlIIBj4RkUIw8ImIFIKBT0SkEF4F/t69e7F0\n6VIkJSUhJiYGycnJyM7OxoEDB+R97HY7TCYT9Ho9YmNjodPpUFJSgt7e3oAVnoiIvCcJIYSnnTZs\n2IC6ujqkpKTA4XCgtbUV/f39AIDGxkYsXboUxcXFMJvNiIqKgsFgQHt7O+x2O7Kzs9HQ0ABJkgZP\nKknw4rRERDSEr9npVQ9/165duHLlCpqbm9HS0oJdu3bJJ//vf/+L5uZmmM1mAEBFRQVaW1tRXV0N\nADh27BhqamomXEAiIvIPrwJfrVbj+PHjWL58ORYuXIgXXngBGo0Gr7/+OnJyclBXVwfA+QcgPz8f\nAJCXlwe1Wg0AqK+vD1DxiYjIW9He7tjT04PGxkb5LcW1a9fQ1NQEu92Orq4ueT+tVgsAUKlU0Gg0\nuHTp0rDtLuXl5fJzo9EIo9E48VoQEUUgi8UCi8Xit+N5NYY/VHd3N7Zu3Yq33noLALB79240NTXh\nD3/4AyRJQn9/vzxer9PpcPnyZeTm5uLo0aODJ+UYPhHRuE3KGP5QM2bMwNatW+X1U6dOYc6cOfJ6\nd3c3AMDhcMBqtQIAUlNTJ1xAIiLyD4+B39fXhz179uD27dvya++99578XKvVIjc3FwAghJAv1tbW\n1sJmswGAvJ2IiILH45DOjRs3kJSUBLVaDb1eD5vNhvPnzwMAkpKS5B5+UVER9u/fD5VKJU/L7O/v\nR1ZW1l1jUBzSISIav4AP6cTHx2PdunVITU3FhQsXcPHiRcydOxdPP/00Ghsb5eGcyspKbN68Gamp\nqejo6IBWq8XGjRtx5MiRCReOiIj8Z9wXbf1yUvbwiYjGbdIv2hIRUXhi4BMRKQQDn4hIIRj4REQK\nwcAnIlIIBj4RkUIw8ImIFIKBT0SkEAx8IiKFYOATESkEA5+ISCEY+ERECsHAJyJSCAY+EZFCMPCJ\niBSCgU9EpBAMfCIihWDgExEpBAOfiEghGPhERArBwCciUojoYBeAKNxIkvvXhZjcchCNF3v4REQK\nwcAnGichnMuGDcDevYPrRKGOgU9EpBAcw6eA4Dg3Uejx2MPftm0bHnvsMcyePRtqtRo6nQ6rV6/G\nmTNn5H3sdjtMJhP0ej1iY2Oh0+lQUlKC3t7egBaeiIi8Jwkxdp8rLS0NXV1duP/++6FSqXDu3DkA\nQEJCAs6cOYO5c+eiuLgYZrMZUVFRMBgMaG9vh91uR3Z2NhoaGiCN6O5JkgQPp6UI0dICPP208zHS\nPPMMsHKl85FoMvianR57+M8++yw+//xznDt3Dm1tbdi2bRsA4NatWzh8+DCam5thNpsBABUVFWht\nbUV1dTUA4NixY6ipqZlw4YiIyH88Bn5paSnS09Pl9ccff1x+HhcXh7q6OgDOvzz5+fkAgLy8PKjV\nagBAfX29XwtMREQTM+6Lttu3bwcAaDQaFBQUoLS0VN6m1WoBACqVChqNBpcuXUJXV5fb45SXl8vP\njUYjjEbjeItCRBTRLBYLLBaL347ndeDfuXMHzz77LKqqqnDvvfeipqYGGo1m1P09jTMNDXwiIrrb\nyM6wyWTy6XhezcO/du0acnJyUFVVhVmzZsFisWDlypUAgNTUVHm/7u5uAIDD4YDVar1rOxERBY/H\nwG9ra8OyZctw4sQJLF68GI2NjVi0aJG8PTc3F4CzR++6WFtbWwubzTZsOxERBZfHaZmZmZn49NNP\nAQAPPfQQpkyZIm977rnn8Mwzz6CoqAj79++HSqWSp2X29/cjKyvL7fgTp2UqR6RNyxw5xXgk/rum\nQPI1Oz2O4dtsNvkf+dmzZ4dty8vLAwBUVlbCYDBg37596OjogFarRUFBAbZs2TLhghF5w5W/SstZ\npdabfOOxhx+Qk7KHrxiB7uFPdvAN9vBHnlD6uhyTUxAGvjIFvIdPFA5+8INgl8ApVMpB5A4Dn/zO\n3Tj30JcC0Qt+6im/H9Ktd98Nj3IQucMhHQUK9HDAZF7YHFqXyRjm4JAOBROHdCiEuQ/FyBBJdSGl\nYOAr2L33hv9577nH/8cMB+zZ00RwSEeBXMMBN24E5vjTpo097HHjRmDaftq0r886Sf+0eHtkmmwc\n0qEJi4QePhF5j4FPAcRxbqJQwsBXoEgdTYvUehH5i1d3yyQaDyGEvJw6JbBokRj2GhEFBwOfiEgh\nGPhERArBwCciUggGPhGRQihils5ot3bh9UMiUhL28ImIFEIRgS+Eczl3DjAYBteJJkKSnMuf/gQ8\n++zgOlGoU0TgExGRQsbwifyJ7w4pXDHwKSBGDnHwCzuIgo9DOkRECsEePgUEe/JEoYeBPwZv5+9z\nnj8RhYOID3x3X6g99CXevZGIlCLsx/AnYw70yy8Dv/3t6PP3Xa/39gJTpnCePxGFpojv4Q9y//2q\nRERK4THwjx8/jjfeeAMnT57E1atXAQBlZWUoKyuT97Hb7XjttddQWVmJixcvQqvVorCwEK+++ioS\nExMDV/oh1q+f3J8jIgo3HgO/ubkZf/vb3zBv3jw58EeOi2/YsAFmsxlRUVEwGAxob29HRUUFWlpa\n0NDQ4HYc3d+ysty//s47E/s5b37WH3jBl4gmjfDAarWKvr4+0dvbKyRJEpIkCZPJJG9vamqSX9+5\nc6cQQoj3339ffu3QoUN3HdOL03rNNWI++nZ8vYgRCzyWw7Xvyy8L8dvfei5Lb68QU6ZMrPwjFyKi\nkXzNTo8XbZOSkhAXFzfqbJa6ujoAzl5/fn4+ACAvLw9qtRoAUF9f7/tfJb+QRiyhwRXxv/sd8OKL\nvOBLRIHj8yydrq4u+blWq3UeVKWCRqO5a3sgBCogncNQzuWNNyS88ooESRpchu7nWhITJXz1lfv9\niIiCLWCzdEZ7R+BSXl4uPzcajTAajQErhyt3z50DVq0CPv307v0Cda+XiRyX950hIgCwWCywWCx+\nO57XgT9abzU1NVV+3t3djZkzZ8LhcMBqtd61faihgR/avJ3OyWmfRORfIzvDJpPJp+N5HfhDe+xD\nn+fm5qK0tBRCCFRXV+PnP/85amtrYbPZ5O2h4pNPgDt3gNOnR9/n4EHvjuXtfi5jnRMALl8Grl3z\nvB8R0URJwsPYy6FDh/Dyyy9DCIGOjg4AwPTp0zF9+nQsX74cVVVVKCoqwv79+6FSqeRpmf39/cjK\nynL7dkSSpEm9pYHrzcm8eUBnJ/DAA3fv8/HHzse4OOfj7duuHrr7nntcnPBqP9frCxe6O+fY7wIm\n83dERKHP1+z02MO/efMmzp8/P+wi5I0bN3Djxg15uKayshIGgwH79u1DR0cHtFotCgoKsGXLlgkX\nLBDef985hu+uF+36o9DXN3x9NLdvj+/cY53T3zi3n4jc8djDD8hJg9TDH89FW88zbFzl924/d9Ud\nPIf7dwcT/R0x8IkiU8B7+JHA9ftxF/Qj9/G31193Pr7xxvh/diI/M/Scv/iF85FBT0SAQgJ/Il+3\nN/Sv6C9+ASQlDQbo0GP8z/8M7nfnDvB//we89NLgfteuTbTUvv0sEdFIihrSGcnbIowV+EOPcesW\noNU6H70rV2CGdMYqIxGFLw7peCH0A49z9oko8BQR+IEQqn9E+A1fRDSasP/Gq3AmhtzM83e/E3jx\nRQEhBpdQNxnfNkZE/sMefsTirR6IaDgG/hhG9l5fecX5OLLzPZFZQCMNDAAOB2C3j6+M4xXo4xNR\n6FLELJ2J8nZ2jy+zgFw/q1I594+K8r587vT3jz3zJzraf7/3/v6vzxT6TUkUEThLJ4C8/b36I/De\nfNP5wbD//V/fjuNpTN2fPXyO3xOFFwZ+xGIaE9FwHNIJMn/f98bTPYD4eycKX75mZ0hOy3RN9xu5\njHcfJRo6rbOrS2D27PCa6klEgcMhnSBjBhPRZAnJHr7ri8ntduesFXdfVO56raYG+O53A/dl5kRE\nkSIkA5+IiPwv5IZ0vLkXDO8XQ0Q0fiEX+P6k1NsD++OTv0QUeUI48L25FwzvF0NE5K0QDnz37r/f\nP/tEMvbkicidkPvgladvgfrsMwGDwfM+AGAwfL0XA5CIIoCvH7wKu8AfftF27K8G5Ng1EUWSCL55\nmjfj8RyzJyLyVggHvu/YsyciGhRyH7waet8Xu10gKurue8EMXa+pEfjud3m/GCIiT0Kyh+/NPHLO\nNSciGh+/9vD/8pe/YMmSJYiPj0dSUhIKCwvR3t7uz1MQEdEE+S3w//jHP2Lt2rVoaWnB7NmzIYRA\ndXU1vvnNb6K7u3tcx3LdCG3kMt59gsVisQS7CAHF+oW3SK5fJNfNH/wS+Hfu3MErX3/Dd0FBAT7/\n/HO0trZi6tSp+PLLL/Haa6/54zRhI9L/0bF+4S2S6xfJdfMHvwT+yZMnYbVaAQD5+fkAgJSUFCxf\nvhwAUF9f74/TEBGRD/wS+F1dXQCcHwrQarXy667nru1ERBREwg/2798vJEkSKpVKNDQ0yK8XFRUJ\nSZJEfHz8sP3h/IgsFy5cuHAZ5+ILv0zLTE1NBZwlGXaB9ssvvxy23UWEytVVIiIF8cuQzje+8Q3c\nd999AIDq6moAwOXLl/Hvf/8bAJCbm+uP0xARkQ/8dvO0PXv24PnnnwcApKWlwWq14ubNm0hOTsZH\nH32EmTNn+uM0REQ0QX6bh//cc8+hqqoKDz/8ML744gtERUXhqaeewokTJ+Swj9QPZpWXl0OlUrld\nHA5HsIs3LsePH8eqVaswY8YMuQ4mk2nYPna7HSaTCXq9HrGxsdDpdCgpKUFvb2+QSu09b+pnNBrd\ntuUjjzwSpFJ7b9u2bXjssccwe/ZsqNVq6HQ6rF69GmfOnJH3Cef286Z+4dx+e/fuxdKlS5GUlISY\nmBgkJycjOzsbBw4ckPfxqf18ugIwDnv37hWSJAlJkkRGRoaYNm2akCRJzJgxQ3zxxReTVYyAKCsr\nE5IkCa1WK1asWDFsGRgYCHbxxmXHjh0iOjpaLFiwQG4vk8k0bJ9169YJSZJEdHS0mD9/voiNjRWS\nJAmj0SgcDkeQSu4db+qXnZ0tJEkS999//7C2/OlPfxqkUntv7ty5QqVSiXnz5onMzEy5jomJiaKz\ns1MIEd7t5039wrn91q9fL2bOnCkWL14sFi1aJGJiYuQ6njx5UgjhW/tNSuDbbDah0WiEJEmisLBQ\nCCHE5cuXxT333CMkSRIbN26cjGIEjCvw169fH+yi+MxqtYq+vj7R29vrNhCbmprk13fu3CmEEOL9\n99+XXzt06FCwiu4VT/UTYjAwKisrg1TKiXv11VfF+fPn5fXt27fL9dyxY0fYt99Y9XvzzTeFEOHd\nfrdv3x627uooq1Qq8Y9//MPn9puUu2Uq5YNZBw8eRHx8PFJSUrBq1Sq0tLQEu0jjlpSUhLi4uFFn\nUtXV1QFwfubC1ZZ5eXlQq9UAQr8tPdVvqJdeeglqtRp6vR7PP/+8POsslJWWliI9PV1ef/zxx+Xn\ncXFxYd9+Y9XPVQeXcGw/tVqN48ePY/ny5Vi4cCFeeOEFaDQavP7668jJyfG5/SYl8CP9g1mSJCEq\nKgopKSnQ6/Xo7u7G0aNHsWLFirAM/bEMbStX+6lUKmg0mru2hytJkjBlyhTodDrMmDEDnZ2d2LNn\nD1asWIGvvvoq2MUbl+3btwMANBoNCgoKIq79htavsLAQQPi3X09PDxobG9Ha2or+/n5cu3YNTU1N\nsNvtPrdfUO+H700vKxysXbsWV69exblz53D27Fn5r6zNZsPOnTuDXLrJESltCQA7duxAT08PTp8+\njQsXLuCXv/wlAKCjowOHDx8Ocum8c+fOHfzoRz9CZWUl7r33XtTU1Mih4E64tZ+7+rmmhod7+33v\ne9+Dw+HA5cuX8eKLLwIADhw4gMrKyiFf7zqct+03KYE/3g9mhRuDwYBp06bJ69/61reQlJQEIPx6\nTC6j/cMa2lautnQ4HPKQXbi05Wj1A4CHH34YMTEx8vqaNWvk5+HQnteuXUNOTg6qqqowa9YsWCwW\nrFy5EkBktN9Y9QPCv/1cZsyYga1bt8rrp06dwpw5c+T1ibTfpAR+pH8wa8eOHbhy5Yq8/sEHH+D6\n9esAnJ9JCEdDewxDn7vaSnx9+2sAqK2thc1mG7Y91I1Wv6tXr2LXrl3D3voPnRIX6u3Z1taGZcuW\n4cSJE1i8eDEaGxuxaNEieXu4t5+n+oVz+/X19WHPnj24ffu2/Np7770nP9dqtb63nx8vMI9p9+7d\n8pXk9PR0eYaOVqsVV65cmaxiBIRrqlhqaqqYP3++XM+pU6eKtra2YBdvXKqrq0VGRobQ6/VyPZKS\nkkRGRoYoKioSQgixdu1aIUmSiIqKEpmZmfLUsezs7OAW3gue6tfZ2SkkSRIxMTEiMzNTzJkzR97v\nwQcfFDabLdhVGNMDDzwgl3fhwoVi2bJl8rJ3714hRHi3n6f6hXP79fT0CEmSRFxcnFiwYIHIyMiQ\ny37fffeJCxcuCCF8a79JC3whhDCbzWLx4sUiPj5eTJ8+XeTn54vPPvtsMosQELt37xZPPPGEmD17\ntoiPjxd6vV4UFxeLTz/9NNhFG7d33nlHngY2cnn00UeFEELY7XZRVlYm0tPThVqtFrNnzxabNm0S\nN2/eDHLpPfNUv1u3bonS0lKxbNkyodFoREJCgliwYIH41a9+JW7cuBHs4nuUlpbmtm4qlUqefhrO\n7eepfuHcfrdv3xbFxcVi3rx5IjExUajVapGWlibWr18v2tvb5f18aT+/3VqBiIhCW1Bn6RAR0eRh\n4BMRKQQDn4hIIRj4REQKwcAnIlIIBj4RkUL8P5atFVUvsqpOAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x105742a90>"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now define the chi^2 function to be minimzed. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def chi_errxy(a,x,unc_x,y,unc_y):\n",
      "    numer = (y-(a[0]+a[1]*x))**2\n",
      "    denom = unc_y**2 + (a[1]*unc_x)**2\n",
      "    chi2=numer/denom\n",
      "    return(chi2.sum())\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "and pass it to scipy.optimize.minimize"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "spo.minimize(chi_errxy,[1.5,1.2],args=(xi,x_unc,ci,obs_unc))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 27,
       "text": [
        "  status: 0\n",
        " success: True\n",
        "    njev: 9\n",
        "    nfev: 36\n",
        "     fun: 8.446592441153806\n",
        "       x: array([ 1.16276895,  1.41272227])\n",
        " message: 'Optimization terminated successfully.'\n",
        "    hess: array([[ 0.81804591, -0.10377374],\n",
        "       [-0.10377374,  0.03095001]])\n",
        "     jac: array([  1.19209290e-07,  -1.19209290e-07])"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}