sigma.clipping.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "from astropy.io import fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# reset defalult plotting values\n",
    "plt.rcParams['figure.figsize'] = (10, 5)\n",
    "plt.rc('font', family='sans-serif')\n",
    "plt.rc('axes', labelsize=14)\n",
    "plt.rc('axes', labelweight='bold')\n",
    "plt.rc('axes', titlesize=16)\n",
    "plt.rc('axes', titleweight='bold')\n",
    "plt.rc('axes', linewidth=2)\n",
    "plt.rc('xtick',labelsize=14)\n",
    "plt.rc('ytick',labelsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# $\\sigma$-Clippinig\n",
    "## removing outliers from your data\n",
    "\n",
    "![]()\n",
    "\n",
    "### Prof. Robert Quimby\n",
    "&copy; 2019 Robert Quimby"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## In this tutorial you will..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Identify outliers from a sample population\n",
    "- Compare sample mean and STD with and without the outliers\n",
    "- Employ iterative techniques to identify outliers\n",
    "- Check if \"outliers\" are actually consistent with the expected distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Lets use a real bias image for our sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# load a bias frame\n",
    "bias = fits.getdata('media/bias.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# plot a histogram of the pixel values\n",
    "bmin = ????\n",
    "bmax = ????\n",
    "plt.hist(bias.flatten(), range=(bmin, bmax), bins=bmax-bmin+1, width=10);\n",
    "plt.xlabel('Counts'), plt.ylabel('Number'), plt.title('Count Distribution (Bias Frame)');\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Checking for Outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "med = ????\n",
    "std = ????\n",
    "\n",
    "# locate outliers\n",
    "w = ????\n",
    "print(\"found {} outliers\".format(w.sum()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# show the outliers\n",
    "j, i = np.indices(bias.shape)\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(bias, vmin=med-std, vmax=med+5*std, origin='lower', cmap='gray');\n",
    "plt.plot(i[w], j[w], 'ro');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## $\\sigma$-clipping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* one (often dangerous!) way to get rid of outliers in a distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# calculate the sample mean and STD (including outliers)\n",
    "mean = np.mean(bias)\n",
    "std = bias.std(ddof=1)\n",
    "print(mean, std, bias.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# identify the (non)outliers\n",
    "w = ????\n",
    "\n",
    "# now remove outliers and re-calculate mean, STD\n",
    "sample = ????\n",
    "print(sample.mean(), sample.std(), sample.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Iterative $\\sigma$-clipping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* can be very (dangerous) useful"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sigclip(sample, nsigma=5.0, niter=100):\n",
    "    # check for outliers\n",
    "    med = np.median(sample)\n",
    "    std = np.std(sample, ddof=1)\n",
    "    isok = ????\n",
    "\n",
    "    # are we done sigma clipping?\n",
    "    if (niter <= 0) or (isok.sum() == sample.size):\n",
    "        return sample\n",
    "\n",
    "    # recursive call\n",
    "    return ????"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# test it out\n",
    "sample = sigclip(bias)\n",
    "print(\"removing {} outliers:\".format(????))\n",
    "print(\"  bias mean is {:.4f}\".format(sample.mean()))\n",
    "print(\"  bias std is {:.4f}\".format(sample.std()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## $\\sigma$-clipping is a common task"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from astropy.stats import ????\n",
    "\n",
    "sample = ????\n",
    "print(\"bias mean is {:.4f}\".format(sample.mean()))\n",
    "print(\"bias std is {:.4f}\".format(sample.std()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Should I clip at $3\\sigma$, $4\\sigma$, or $5\\sigma$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Think about your data\n",
    "- How many outliers would you expect by random chance?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Consider the cumulative distribution of the parent population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-5, 5, 100)\n",
    "plt.plot(x, ????)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# how many bias pixels do we expect to be n-sigma larger than the mean?\n",
    "nsigma = ????\n",
    "ndata = ????\n",
    "expect = ????\n",
    "print(\"expect {:.2f} sample points {}-sigma from the mean\".format(expect, nsigma))"
   ]
  }
 ],
 "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.6.10"
  },
  "livereveal": {
   "scroll": "True",
   "start_slideshow_at": "selected"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}