lsq-fits.with.a.shift.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Least-Squares Fitting with a Shift\n",
    "\n",
    "### Prof. Robert Quimby\n",
    "&copy; 2018 Robert Quimby"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## In this tutorial you will..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Quickly review how to fit a line to data with least-squares\n",
    "* Consider the uncertainty in the best-fit model\n",
    "* Find that model uncertainty can often be minimized with the right choice of coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Load some data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "data = np.genfromtxt('media/xy.dat', names='x, y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# setup plotting\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = (10, 7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# make a scatter plot\n",
    "plt.????\n",
    "plt.xlabel('x value', weight='bold', fontsize=14)\n",
    "plt.ylabel('y value', weight='bold', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recall how to least-squares fit a line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume the model:\n",
    "$$y_i = mx_i + b$$\n",
    "\n",
    "where $x_i$ and $y_i$ are the x and y values, respectively for the $i^{\\rm th}$ observation. Or in matrix form:\n",
    "\n",
    "$$\n",
    "\\left[ \\begin{array}{c}\n",
    "y_1  \\\\\n",
    "y_2  \\\\\n",
    "\\vdots \\\\\n",
    "y_N  \\end{array} \\right]\n",
    "= \n",
    "\\left[ \\begin{array}{cc}\n",
    "x_1 & 1 \\\\\n",
    "x_2 & 1 \\\\\n",
    "\\vdots & \\vdots\\\\\n",
    "x_N & 1 \\end{array} \\right] \n",
    "\\left[ \\begin{array}{c}\n",
    "m \\\\\n",
    "b \\end{array} \\right] $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "X = ????\n",
    "Y = ????\n",
    "p = ????"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plt.scatter(data['x'], data['y'], c='r')\n",
    "plt.xlabel('x value', weight='bold', fontsize=14)\n",
    "plt.ylabel('y value', weight='bold', fontsize=14);\n",
    "plt.plot(????, ls='dashed');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Determine the model parameter uncertainties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "modely = X * p\n",
    "dY = Y - modely\n",
    "M = X.shape[0] # number of data samples\n",
    "N = p.shape[0] # number of model parameters\n",
    "sample_var = (dY.T * dY) / (M - N)\n",
    "pvar = (sample_var * np.diagonal((X.T * X).I))\n",
    "psig = np.sqrt(pvar)\n",
    "\n",
    "# best fit parameters\n",
    "m, b = p.A1\n",
    "msig, bsig = psig.A1\n",
    "\n",
    "# now for the parameter errors\n",
    "var = (sample_var * np.diagonal((X.T * X).I))\n",
    "msig, bsig = np.sqrt(var).A1\n",
    "print(\"m is {:.2f} +/- {:.2f}\".format(m, msig))\n",
    "print(\"b is {:.2f} +/- {:.2f}\".format(b, bsig))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Randomly generate some possible models within the uncertainty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_models(data, m, msig, b, bsig, n=30):\n",
    "    # show n realizations of the model with errors\n",
    "    ms = ????\n",
    "    bs = ????\n",
    "    modelx = np.array([min(data['x'])-1, max(data['x'])+1])\n",
    "    for thism, thisb in zip(ms, bs):\n",
    "        plt.plot(modelx, thism * modelx + thisb, color='0.75')\n",
    "    plt.plot(data['x'], data['y'], 'ro')\n",
    "    plt.plot(modelx, modelx * m + b, '--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plot_models(data, m, msig, b, bsig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## One more time with a shifted origin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$(y - y_0) = m(x - x_0) + b$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "cx = np.mean(data['x'])\n",
    "cy = np.mean(data['y'])\n",
    "X = np.matrix( [[????, 1] for xi in data['x']] )\n",
    "Y = np.matrix( [[????] for yi in data['y']] )\n",
    "p = (X.T * X).I * X.T * Y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# predicted y values\n",
    "modely = X * p\n",
    "\n",
    "# residuals from observed values\n",
    "dY = Y - modely\n",
    "\n",
    "# variance of the data from the model\n",
    "M = X.shape[0] # number of data samples\n",
    "N = p.shape[0] # number of model parameters\n",
    "var = (dY.T * dY) / (M - N)\n",
    "\n",
    "# best fit parameters\n",
    "m2 = p[0, 0]\n",
    "b2 = p[1, 0]\n",
    "\n",
    "# now for the parameter errors\n",
    "err_p = np.sqrt(var * np.diagonal((X.T * X).I))\n",
    "\n",
    "msig2 = err_p[0, 0]\n",
    "bsig2 = err_p[0, 1]\n",
    "print(\"m is {:.2f} +/- {:.2f}\".format(m2, msig2))\n",
    "print(\"b is {:.2f} +/- {:.2f}\".format(b2, bsig2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def plot_models2(data, m, msig, b, bsig, n=30, cx=0, cy=0):\n",
    "    # show n realizations of the model with errors\n",
    "    ms = np.random.normal(m, msig, n)\n",
    "    bs = np.random.normal(b, bsig, n)\n",
    "    modelx = np.array([min(data['x'])-1, max(data['x'])+1])\n",
    "    besty = modelx * m + b\n",
    "    for thism, thisb in zip(ms, bs):\n",
    "        plt.plot(modelx, (modelx - cx) * thism + thisb + cy, color='0.75')\n",
    "    plt.plot(data['x'], data['y'], 'ro')\n",
    "    plt.plot(modelx, besty, '--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plot_models2(data, m2, msig2, b2, bsig2, cx=cx, cy=cy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fit comparision"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(\"m  is {:.2f} +/- {:.2f}\".format(m, msig))\n",
    "print(\"m2 is {:.2f} +/- {:.2f}\".format(m2, msig2))\n",
    "print(\"b  is {:.2f} +/- {:.2f}\".format(b, bsig))\n",
    "print(\"b2 is {:.2f} +/- {:.2f}\".format(b2, bsig2))\n",
    "print(\"value of y at x=0 is {:.2f} (shifted coordinates)\".format(cy - (cx * m2) + b2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  },
  "livereveal": {
   "scroll": "True",
   "start_slideshow_at": "selected"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}