lectures/Lecture_14.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Density Estimation & Clustering\n",
    "\n",
    "*S. R. Taylor (2021)*\n",
    "\n",
    "This lecture and notebook are based on the \"DensityEstimation\" and \"DE2andClustering\" lectures of of G. Richards' \"Astrostatistics\" class at Drexel University (PHYS 440/540, https://github.com/gtrichards/PHYS_440_540).\n",
    "\n",
    "##### Reading:\n",
    "\n",
    "- [Textbook](http://press.princeton.edu/titles/10159.html) Chapter 6.\n",
    "\n",
    "---\n",
    "\n",
    "## Contents\n",
    "* [Density Estimation](#one)\n",
    "    * [Non-parametric density estimation](#onea)\n",
    "    * [Parametric density estimation](#oneb)\n",
    "* [Clustering Algorithms](#two)\n",
    "    * [$K$-means](#twoa)\n",
    "    * [Mean-shift](#twob)\n",
    "* [Correlation functions](#three)\n",
    "    \n",
    "---\n",
    "\n",
    "***Exercises required for class participation are in <font color='red'>red</font>.***\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Density Estimation <a class=\"anchor\" id=\"one\"></a>\n",
    "\n",
    "Inferring the pdf of a sample of data is known as ***density estimation***. Essentially we are smoothing the data to correct for the finiteness of our sample and to better recover the underlying distribution.\n",
    "\n",
    "Density estimation is useful because:\n",
    "- identifying low probability regions can help uncover rare sources. \n",
    "- if the data can be divided into sub-samples, one can estimate the pdf for each subsample and, in turn determine classifications for new objects.\n",
    "\n",
    "### Non-parametric Density Estimation <a class=\"anchor\" id=\"onea\"></a>\n",
    "\n",
    "*Nonparametric* density estimation is useful when we know nothing about the underlying distribution of the data, since we don't have to specify a functional form. This flexibility allows us to capture the shape of the distribution well, at the expense of more difficulty interpreting the results.\n",
    "\n",
    "#### Kernel Density Estimation (KDE)\n",
    "\n",
    "[*Kernel Density Estimation (KDE)*](https://en.wikipedia.org/wiki/Kernel_density_estimation) is the standard approach for non-parametric density estimation.\n",
    "\n",
    "Let's start by recalling some problems with making simple histograms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------\n",
    "# Draw some random data\n",
    "np.random.seed(1)\n",
    "x = np.concatenate([np.random.normal(-0.5, 0.3, size=14), \n",
    "                    np.random.normal(1, 0.3, size=7)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Modified from Ivezic, Figure 6.1, modified by GTR\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "#------------------------------------------------------------\n",
    "\n",
    "# First figure: silly histogram binning\n",
    "fig1 = plt.figure(figsize=(8, 4))\n",
    "fig1.subplots_adjust(left=0.12, right=0.95, wspace=0.05, \n",
    "                     bottom=0.15, top=0.9, hspace=0.05)\n",
    "\n",
    "FC = '#6666FF'\n",
    "XLIM = (-2, 2.9)\n",
    "YLIM = (-0.09, 1.1)\n",
    "\n",
    "ax = fig1.add_subplot(121)\n",
    "bins = np.linspace(-1.8, 2.7, 13)\n",
    "ax.hist(x, bins=bins, density=True, \n",
    "        histtype='stepfilled', fc='k', alpha=0.3)\n",
    "ax.plot(XLIM, [0, 0], '-k', lw=1)\n",
    "ax.plot(x, 0 * x - 0.05, '+k')\n",
    "ax.set_xlim(XLIM)\n",
    "ax.set_ylim(YLIM)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$p(x)$')\n",
    "\n",
    "#Shift bin centers by 0.25\n",
    "ax = fig1.add_subplot(122)\n",
    "ax.yaxis.set_major_formatter(plt.NullFormatter())\n",
    "ax.hist(x, bins=bins + 0.25, density=True, \n",
    "        histtype='stepfilled', fc='k', alpha=0.3)\n",
    "ax.plot(XLIM, [0, 0], '-k', lw=1)\n",
    "ax.plot(x, 0 * x - 0.05, '+k')\n",
    "ax.set_xlim(XLIM)\n",
    "ax.set_ylim(YLIM)\n",
    "ax.set_xlabel('$x$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The data points that make up the histogram are the same in both panels.  All we have done is shifted the locations of the bins by 0.25. ***The choice of number of bins and the location of bin centers can really change the histogram that we make.***\n",
    "\n",
    "The next panels are what happens if we center the bins on each point. This is an example of **kernel density estimation** using a \"***top-hat***\" kernel. It is a good description of the data, but pretty ugly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig1b = plt.figure(figsize=(8, 4))\n",
    "fig1b.subplots_adjust(left=0.12, right=0.95, wspace=0.05, \n",
    "                      bottom=0.1, top=0.95, hspace=0.05)\n",
    "\n",
    "ax = fig1b.add_subplot(111)\n",
    "ax.xaxis.set_major_formatter(plt.NullFormatter())\n",
    "binwidth = bins[1] - bins[0]\n",
    "x_plot = np.linspace(-4, 4, 1000)\n",
    "y_plot = (abs(x_plot - x[:, None]) <= 0.5 * binwidth).astype(float)\n",
    "y_plot /= (binwidth * len(x))\n",
    "ax.fill(x_plot, y_plot.sum(0), ec='k', lw=1, fc='k', alpha=0.3)\n",
    "ax.plot(x_plot, y_plot.T, '-k', lw=1)\n",
    "ax.plot(x, 0 * x - 0.05, '+k')\n",
    "ax.set_xlim(XLIM)\n",
    "ax.set_ylim(YLIM)\n",
    "ax.set_ylabel('$p(x)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can make it look nicer by choosing a different kernel, i.e. a different bin shape. The next plot shows a **KDE using a Gaussian kernel**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------\n",
    "# First figure: transition to KDE\n",
    "fig2 = plt.figure(figsize=(8, 8))\n",
    "fig2.subplots_adjust(left=0.12, right=0.95, wspace=0.05, \n",
    "                     bottom=0.0, top=1.0, hspace=0.05)\n",
    "\n",
    "## Just right\n",
    "ax = fig2.add_subplot(313)\n",
    "ax.yaxis.set_major_formatter(plt.NullFormatter())\n",
    "binwidth = bins[1] - bins[0]\n",
    "x_plot = np.linspace(-4, 4, 1000)\n",
    "y_plot = binwidth * stats.norm.pdf(x_plot, x[:, None], 0.2)\n",
    "y_plot /= (binwidth * len(x))\n",
    "ax.fill(x_plot, y_plot.sum(0), ec='k', lw=1, fc='k', alpha=0.3)\n",
    "ax.plot(x_plot, y_plot.T, '-k', lw=1)\n",
    "ax.plot(x, 0 * x - 0.05, '+k')\n",
    "ax.set_xlim(XLIM)\n",
    "ax.set_ylim(YLIM)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$p(x)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can think of KDE as replacing the points with \"clouds\". Each cloud is described by the kernel $K(u)$, where $K(u)$ can be any function that is:\n",
    "- smooth,\n",
    "- postive definite,\n",
    "- normalizes to unity, \n",
    "- has zero mean,\n",
    "- has positive variance. \n",
    "\n",
    "A common kernel is the **Gaussian kernel** that we just used above:\n",
    "\n",
    "$$ K(u) = \\frac{1}{(2\\pi)^{D/2}}\\exp{(-u^2/2)}$$\n",
    "\n",
    "where $D$ denotes the dimensionality of the data. Once a kernel is chosen the KDE at a point, $x$, is given by \n",
    "\n",
    "$$ \\hat{f}(x) = \\frac{1}{Nh^D}\\sum_{i=1}^N K\\left(\\frac{d(x,x_i)}{h}\\right),$$\n",
    "\n",
    "where $\\hat{f}$ is an ***estimator*** of our distribution.\n",
    "\n",
    "The argument of $K$ is just some measure of the distance between $x$ and each $x_i$. Normally $d(x,x_i) = (x-x_i)$. For the gaussian kernel that makes $h=\\sigma$. So, $h$ represents the \"width\" or what is usually called the **\"bandwidth\"** in this context.\n",
    "\n",
    "Here is a comparison of some different possible kernels.\n",
    "\n",
    "![Ivezic, Figure 6.2](http://www.astroml.org/_images/fig_kernels_1.png)\n",
    "\n",
    "**The Epanechnikov kernel is \"optimal\" because it minimizes the variance of the kernel density estimate**: \n",
    "\n",
    "$$K(x) = \\frac{3}{4}(1-x^2),$$\n",
    "\n",
    "for $|x|\\le 1$ and 0 otherwise. Below is the code that produces the plot above.  \n",
    "\n",
    "<font color='red'>Add the Epanechnikov kernel to it in the color and linestyle of your choice by modifying the cell below.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------\n",
    "# Compute Kernels.\n",
    "xgrid = np.linspace(-5, 5, 10000)\n",
    "dx = xgrid[1] - xgrid[0]\n",
    "\n",
    "# Gaussian kernel\n",
    "gauss = (1. / np.sqrt(2 * np.pi)) * np.exp(-0.5 * xgrid ** 2)\n",
    "\n",
    "# Exponential kernel\n",
    "exp = 0.5 * np.exp(-abs(xgrid))\n",
    "\n",
    "# Top-hat kernel\n",
    "tophat = 0.5 * np.ones_like(xgrid)\n",
    "tophat[abs(xgrid) > 1] = 0 # Range of the tophat kernel\n",
    "\n",
    "# Epanechnikov kernel \n",
    "ep = ____  # Add the Epanechnikov kernel function\n",
    "ep[____] = 0 # Set the range of the kernel\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot the kernels\n",
    "fig = plt.figure(figsize=(5, 3.75))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "ax.plot(xgrid, gauss, '-', c='black', \n",
    "        lw=3, label='Gaussian')\n",
    "ax.plot(xgrid, exp, '-', c='#666666', \n",
    "        lw=2, label='Exponential')\n",
    "ax.plot(xgrid, tophat, '-', c='#999999', \n",
    "        lw=1, label='Top-hat')\n",
    "# Add the Epanechnikov kernel to the plot\n",
    "ax.plot(____,____,____,____,label='Epanechnikov')  \n",
    "\n",
    "ax.legend(loc=1)\n",
    "\n",
    "ax.set_xlabel('$u$')\n",
    "ax.set_ylabel('$K(u)$')\n",
    "\n",
    "ax.set_xlim(-5, 5)\n",
    "ax.set_ylim(0, 0.8001)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### How do we determine the optimal kernel bandwidth?\n",
    "\n",
    "When we first discussed histograms and KDE we talked about theoretical computation of optimal bandwidths.  Let's now look at how we can empirically determine the optimal bandwidth through [cross validation](https://en.wikipedia.org/wiki/Cross-validation_(statistics)).\n",
    "\n",
    "To avoid fooling ourselves into thinking our model is better than it is, we set aside some of the data as the \"***test***\" set--that provides an independent check on the the fit from the ***training data***. \n",
    "\n",
    "That's fine if our model has no parameters, or has parameters that are all fixed ahead of time.  But what if we want to try to use the data to not only fit a model, but decide what the best model parameters are?  In that case, we need to set aside some of the training data as a \"***validation***\" set.  \n",
    "\n",
    "However, not only does that leave us even less data for training, but it also means that it matters which objects are in each of the training, validation, and test sets.  We can solve both of these problems with **cross validation**.\n",
    "\n",
    "With our KDE example, we are just trying to map the density, not predict the value, so we really only need two sets: training and validation. There are a number of different ways to do this.\n",
    "\n",
    "1. You could **randomly sample** to decide which data goes into the training or test sets. We don't just do this once, but rather many times so that each data point is treated both as a training point and as a test point at some stage.\n",
    "\n",
    "\n",
    "2. We could do this in a more ordered way (e.g., to make sure that each object gets sampled as training/test the same number of times) and do a **$K$-fold cross validation**. Here $K$ is the number of \"experiments\" that need to be done so that each data point appears in a test sample once.\n",
    "\n",
    "\n",
    "3. We can take $K$-fold cross validation to the extreme by having $K\\equiv N$, so that in each experiment we leave out just one object. This is called **\"Leave-One-Out\" cross validation**.\n",
    "\n",
    "Let's go back to our problem of trying to determine the best bandwidth for our Kernel Density Estimate. <font color='red'>We can do this in Scikit-Learn with [`GridSearchCV`](http://scikit-learn.org/0.17/modules/generated/sklearn.grid_search.GridSearchCV.html) and replot our histogram above as follows:</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Complete and Execute this cell to determine the bandwidth\n",
    "from sklearn.neighbors import KernelDensity\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "\n",
    "bwrange = np.linspace(____, ____, ____) # Test 30 bandwidths from 0.1 to 1.0\n",
    "K = ____ # Do 5-fold cross validation\n",
    "grid = GridSearchCV(KernelDensity(), {'bandwidth': ____}, cv=____) # Try each bandwidth with K-folds\n",
    "grid.fit(x[:, None]) #Fit the histogram data that we started the lecture with.\n",
    "h_opt = grid.best_params_['bandwidth']\n",
    "print(h_opt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Execute this cell to see the new \"histogram\"\n",
    "fig2 = plt.figure(figsize=(5, 5))\n",
    "ax = fig2.add_subplot(111)\n",
    "ax.xaxis.set_major_formatter(plt.NullFormatter())\n",
    "ax.yaxis.set_major_formatter(plt.NullFormatter())\n",
    "\n",
    "binwidth = bins[1] - bins[0]\n",
    "x_grid = np.linspace(-4, 4, 1000)\n",
    "y_plot = binwidth * stats.norm.pdf(x_grid, x[:, None], h_opt)\n",
    "y_plot /= (binwidth * len(x))\n",
    "\n",
    "ax.fill(x_grid, y_plot.sum(0), ec='k', lw=1, fc='k', alpha=0.3)\n",
    "ax.plot(x_grid, y_plot.T, '-k', lw=1)\n",
    "ax.plot(x, 0 * x - 0.05, '+k')\n",
    "\n",
    "ax.set_xlim(XLIM)\n",
    "ax.set_ylim(YLIM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 2-D Histograms\n",
    "\n",
    "Here is some sample code using [`sklearn.neighbors.KernelDensity`](http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html).  Play around with this and see how it works.  Make different variations of the plot.  (For example, try `bandwidth=0.01` and `bandwidth=1.0`.)  What we are doing here is using KDE to set the plot color to indicate the relative density of the points.  This is essentially a 2-D histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "# Two 2-D normal distributions with offset centroids\n",
    "# See what happens when you make changes to the next 2 lines.\n",
    "X = np.concatenate([np.random.normal([-1,-1],[0.75,0.75], size=(1000,2)), \n",
    "                    np.random.normal([1,1],[1,1], size=(500,2))]) \n",
    "kde = KernelDensity(kernel='gaussian', bandwidth=0.2)\n",
    "kde.fit(X) #fit the model to the data\n",
    "\n",
    "u = v = np.linspace(-4,5,80)\n",
    "Xgrid = np.vstack(map(np.ravel, np.meshgrid(u, v))).T\n",
    "dens = np.exp(kde.score_samples(Xgrid)) #evaluate the model on the grid\n",
    "\n",
    "plt.scatter(Xgrid[:,0], Xgrid[:,1], c=dens, \n",
    "            cmap=\"Purples\", edgecolor=\"None\")\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<font color='red'>Now copy the example from above to a new cell and splice in the cross validation code to produce a new plot with the \"optimal\" bandwidth.  Basically, splice in the lines of code for `GridSearchCV` between the lines setting `X` and instantiating `kde`.  Then replace the hardcoded bandwidth with the optimal value you computed.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "---\n",
    "\n",
    "#### Nearest-Neighbor Density Estimation\n",
    "\n",
    "Another very simple way to estimate the density of an $N$-dimensional distribution is to look to the nearest object (or the $K$ nearest objects) and compute their distances, $d_K$.  This is the [$K$-Nearest Neighbor](https://en.wikipedia.org/wiki/K-nearest_neighbors_algorithm) algorithm. The density at a given point, $x$ is estimated as\n",
    "\n",
    "$$\\hat{f}_K(x) = \\frac{K}{V_D(d_K)},$$\n",
    "\n",
    "where $V_D(d)$ is given generically by $\\frac{2d^D\\pi^{D/2}}{D\\Gamma(D/2)}$ where $\\Gamma$ is the complete gamma function, and this formula reduces to the usual equations for area and volume in 2 and 3 dimensions, respectively.\n",
    "\n",
    "We can simplify this to \n",
    "\n",
    "$$\\hat{f}_K(x) = \\frac{C}{d_K^D}$$\n",
    "\n",
    "since the constant, $C$ can be evaluated at the end.\n",
    "\n",
    "This estimator has some intrinsic bias, which can be reduced by considering *all* $K$ nearest neighbors:\n",
    "$$\\hat{f}_K(x) = \\frac{C}{\\sum_{i=1}^K d_i^D}$$\n",
    "\n",
    "See the [Scikit-Learn `neighbors` documentation](http://scikit-learn.org/stable/modules/neighbors.html) for more information.\n",
    "\n",
    "Ivezic, Figure 6.5 compares a Nearest Neighbor ($k=10$) with a KDE algorithm. <font color='red'>See what happens as you increase the number of neighbors used.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Based on Ivezic, Figure 6.5\n",
    "# Modified by GTR\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "\n",
    "from astroML.density_estimation import KNeighborsDensity\n",
    "from astropy.visualization import hist\n",
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Generate our data: a mix of several Cauchy distributions\n",
    "#  this is the same data used in the Bayesian Blocks figure\n",
    "np.random.seed(0)\n",
    "N = 10000\n",
    "mu_gamma_f = [(5, 1.0, 0.1),\n",
    "              (7, 0.5, 0.5),\n",
    "              (9, 0.1, 0.1),\n",
    "              (12, 0.5, 0.2),\n",
    "              (14, 1.0, 0.1)]\n",
    "true_pdf = lambda x: sum([f * stats.cauchy(mu, gamma).pdf(x) \n",
    "                          for (mu, gamma, f) in mu_gamma_f])\n",
    "x = np.concatenate([stats.cauchy(mu, gamma).rvs(int(f * N)) \n",
    "                    for (mu, gamma, f) in mu_gamma_f])\n",
    "np.random.shuffle(x)\n",
    "x = x[x > -10]\n",
    "x = x[x < 30]\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# plot the results\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "N = 5000\n",
    "k = 10 #Number of neighbors\n",
    "\n",
    "xN = x[:N]\n",
    "t = np.linspace(-10, 30, 1000)\n",
    "\n",
    "# Compute density with KDE\n",
    "kde = KernelDensity(bandwidth=0.1, kernel='gaussian')\n",
    "kde.fit(xN[:, None])\n",
    "dens_kde = np.exp(kde.score_samples(t[:, None]))\n",
    "\n",
    "# Compute density with Bayesian nearest neighbors\n",
    "nbrs = KNeighborsDensity('bayesian', n_neighbors=k)\n",
    "nbrs.fit(xN[:, None])\n",
    "dens_nbrs = nbrs.eval(t[:, None]) / N\n",
    "\n",
    "# plot the results\n",
    "plt.plot(t, true_pdf(t), ':', color='black', \n",
    "         zorder=3, label=\"Generating Distribution\")\n",
    "plt.plot(xN, -0.005 * np.ones(len(xN)), '|k')\n",
    "plt.plot(t, dens_nbrs, '-', lw=1.5, color='gray', \n",
    "         zorder=2, label=\"Nearest Neighbors (k=%i)\" % k)\n",
    "plt.plot(t, dens_kde, '-', color='black', \n",
    "         zorder=3, label=\"Kernel Density (h=0.1)\")\n",
    "\n",
    "# label the plot\n",
    "plt.ylabel('$p(x)$')\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlim(0, 20)\n",
    "plt.ylim(-0.01, 0.4001)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parametric Density Estimation <a class=\"anchor\" id=\"oneb\"></a>\n",
    "\n",
    "#### Gaussian Mixture Models (GMM)\n",
    "\n",
    "KDE centers each bin (or rather, kernel) at each point.  In a [**mixture model**](https://en.wikipedia.org/wiki/Mixture_model) we don't use a kernel for each data point, but rather we fit for the ***locations of the kernels*** in addition to the width. So a mixture model is sort of a hybrid between a tradtional (fixed bin location/size) histogram and KDE. \n",
    "\n",
    "- Using lots of kernels (maybe even more than the BIC score suggests) may make sense if you just want to provide an accurate description of the data (as in density estimation).  \n",
    "- Using fewer kernels makes mixture models more like clustering (later today), where the suggestion is still to use many kernels in order to divide the sample into real clusters and \"background\".\n",
    "\n",
    "Gaussians are the most commonly used components for mixture models.  So, the pdf is modeled by a sum of Gaussians:\n",
    "\n",
    "$$p(x) = \\sum_{k=1}^N \\alpha_k \\mathscr{N}(x|\\mu_k,\\Sigma_k),$$\n",
    "\n",
    "where $\\alpha_k$ are the \"mixing coefficients\" with $0\\le \\alpha_k \\le 1$ and $\\sum_{k=1}^N \\alpha_k = 1$.\n",
    "\n",
    "We can solve for the parameters using maximum likelihood analyis as we have discussed previously.\n",
    "However, this can be complicated in multiple dimensions, requiring the use of [**Expectation Maximization (EM)**](https://en.wikipedia.org/wiki/Expectation%E2%80%93maximization_algorithm) methods (see textbook for details).\n",
    "\n",
    "Ivezic Figure 4.2 (next cell) provides an example in 1-D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ivezic v2, Figure 4.2\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.mixture import GaussianMixture\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# This function adjusts matplotlib settings for a uniform feel in the textbook.\n",
    "# Note that with usetex=True, fonts are rendered with LaTeX.  This may\n",
    "# result in an error if LaTeX is not installed on your system.  In that case,\n",
    "# you can set usetex to False.\n",
    "if \"setup_text_plots\" not in globals():\n",
    "    from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=8, usetex=False)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Set up the dataset.\n",
    "#  We'll create our dataset by drawing samples from Gaussians.\n",
    "\n",
    "random_state = np.random.RandomState(seed=1)\n",
    "\n",
    "X = np.concatenate([random_state.normal(-1, 1.5, 350),\n",
    "                    random_state.normal(0, 1, 500),\n",
    "                    random_state.normal(3, 0.5, 150)]).reshape(-1, 1)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Learn the best-fit GaussianMixture models\n",
    "#  Here we'll use scikit-learn's GaussianMixture model. The fit() method\n",
    "#  uses an Expectation-Maximization approach to find the best\n",
    "#  mixture of Gaussians for the data\n",
    "\n",
    "# fit models with 1-10 components\n",
    "N = np.arange(1, 11)\n",
    "models = [None for i in range(len(N))]\n",
    "\n",
    "for i in range(len(N)):\n",
    "    models[i] = GaussianMixture(N[i]).fit(X)\n",
    "\n",
    "# compute the AIC and the BIC\n",
    "AIC = [m.aic(X) for m in models]\n",
    "BIC = [m.bic(X) for m in models]\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot the results\n",
    "#  We'll use three panels:\n",
    "#   1) data + best-fit mixture\n",
    "#   2) AIC and BIC vs number of components\n",
    "#   3) probability that a point came from each component\n",
    "\n",
    "fig = plt.figure(figsize=(10, 3.4))\n",
    "fig.subplots_adjust(left=0.12, right=0.97,\n",
    "                    bottom=0.21, top=0.9, wspace=0.5)\n",
    "\n",
    "\n",
    "# plot 1: data + best-fit mixture\n",
    "ax = fig.add_subplot(131)\n",
    "M_best = models[np.argmin(AIC)]\n",
    "\n",
    "x = np.linspace(-6, 6, 1000)\n",
    "logprob = M_best.score_samples(x.reshape(-1, 1))\n",
    "responsibilities = M_best.predict_proba(x.reshape(-1, 1))\n",
    "pdf = np.exp(logprob)\n",
    "pdf_individual = responsibilities * pdf[:, np.newaxis]\n",
    "\n",
    "ax.hist(X, 30, density=True, histtype='stepfilled', alpha=0.4)\n",
    "ax.plot(x, pdf, '-k')\n",
    "ax.plot(x, pdf_individual, '--k')\n",
    "ax.text(0.04, 0.96, \"Best-fit Mixture\",\n",
    "        ha='left', va='top', transform=ax.transAxes)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$p(x)$')\n",
    "\n",
    "\n",
    "# plot 2: AIC and BIC\n",
    "ax = fig.add_subplot(132)\n",
    "ax.plot(N, AIC, '-k', label='AIC')\n",
    "ax.plot(N, BIC, '--k', label='BIC')\n",
    "ax.set_xlabel('n. components')\n",
    "ax.set_ylabel('information criterion')\n",
    "ax.legend(loc=2)\n",
    "\n",
    "\n",
    "# plot 3: posterior probabilities for each component\n",
    "ax = fig.add_subplot(133)\n",
    "\n",
    "p = responsibilities\n",
    "p = p[:, (1, 0, 2)]  # rearrange order so the plot looks better\n",
    "p = p.cumsum(1).T\n",
    "\n",
    "ax.fill_between(x, 0, p[0], color='gray', alpha=0.3)\n",
    "ax.fill_between(x, p[0], p[1], color='gray', alpha=0.5)\n",
    "ax.fill_between(x, p[1], 1, color='gray', alpha=0.7)\n",
    "ax.set_xlim(-6, 6)\n",
    "ax.set_ylim(0, 1)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel(r'$p({\\rm class}|x)$')\n",
    "\n",
    "ax.text(-5, 0.3, 'class 1', rotation='vertical')\n",
    "ax.text(0, 0.5, 'class 2', rotation='vertical')\n",
    "ax.text(3, 0.3, 'class 3', rotation='vertical')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try a simpler example so you can see exactly how to use GMM. We'll fit a two-dimensional distribution composed of two Gaussian distributions using GMMs with varying numbers of components. The BIC will tell us how many components are favored:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell\n",
    "import numpy as np\n",
    "from sklearn.mixture import GaussianMixture\n",
    "\n",
    "X = np.concatenate([np.random.normal(size=(500,2)),\n",
    "                   np.random.normal(loc=3.0, size=(500,2))]) # 1000  points in 2D\n",
    "BIC = []\n",
    "for ii in range(1,6):\n",
    "    gmm = GaussianMixture(ii) #three components\n",
    "    gmm.fit(X)\n",
    "    log_dens = gmm.score(X)\n",
    "    BIC.append(gmm.bic(X))\n",
    "    \n",
    "plt.plot(np.arange(1,6), BIC);\n",
    "plt.xlabel('Number of components');\n",
    "plt.ylabel('BIC Score');\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do the results make sense?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustering <a class=\"anchor\" id=\"two\"></a>\n",
    "\n",
    "> *[Clustering](https://en.wikipedia.org/wiki/Cluster_analysis) algorithms attempt to group together similar objects in a data set.* \n",
    "\n",
    "This process allows us to put new objects into the resulting classes and to identify rare objects that do not fit any particular mold. **Clustering is inherently an \"unsupervised\" process** as we do not know the classification of the objects. Since we have no metric for determining when we are right, it is a bit of a dark art, but it also can be very powerful. Scikit-Learn's clustering suite is summarized at [http://scikit-learn.org/stable/modules/clustering.html](http://scikit-learn.org/stable/modules/clustering.html).\n",
    "\n",
    "### $K$-Means Clustering <a class=\"anchor\" id=\"twoa\"></a>\n",
    "\n",
    "We start with [$K$-means clustering](http://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html), which is one of the simplest methods.  $K$-means seeks to minimize the following\n",
    "\n",
    "$$\\sum_{k=1}^{K}\\sum_{i\\in C_k}||x_i - \\mu_k||^2$$\n",
    "\n",
    "where $\\mu_k = \\frac{1}{N_k}\\sum_{i\\in C_k} x_i$\n",
    "\n",
    "In words, this says to\n",
    "  * Take every object in class $C_k$ (as determined by which centroid it is closest to, specifically $C_k(x_i) = \\arg \\min_k ||x_i-\\mu_k||)$\n",
    "  * Compute the mean of the objects in that class\n",
    "  * Subtract that mean from each member of that class and square the norm\n",
    "  * Do that for each class and sum\n",
    "  * Shift the centroids of the *pre-determined* number of classes until this sum is minimized\n",
    "  * Do this multiple times with different starting centroids and take the result with the minimum sum\n",
    "  \n",
    "A typical call will look something like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn import preprocessing\n",
    "\n",
    "#X = np.random.normal(size=(1000,2)) #1000 points in 2D\n",
    "X = np.concatenate([np.random.normal(size=(500,2)),\n",
    "                   np.random.normal(loc=3.0, size=(500,2))]) # 1000  points in 2D\n",
    "clf = KMeans(n_clusters=2) #Try 2 clusters to start with\n",
    "clf.fit(X)\n",
    "centers = clf.cluster_centers_ #location of the clusters\n",
    "labels = clf.predict(X) #labels for each of the points\n",
    "\n",
    "# plot the data color-coded by cluster id\n",
    "colors = ['C0', 'C1', 'C2']\n",
    "for ii in range(3):\n",
    "    plt.scatter(X[labels==ii,0], X[labels==ii,1], \n",
    "                color=colors[ii])\n",
    "\n",
    "# To get some information on these try:\n",
    "# KMeans?\n",
    "# help(clf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X[labels[0],0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example of **K-means clustering on some multi-dimensional data**. The example and data derive from G. Richards' lecture \"DE2andClustering.ipynb\", where he attributes this to Professor Cruz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute the next few cells\n",
    "from astropy.table import Table\n",
    "t = Table.read('./code/cruz_all_dist.dat', format=\"ascii\")\n",
    "\n",
    "# Turn these data into a properly formatted Scikit-Learn array\n",
    "X = np.vstack([ t['col2'], t['col3'], t['col4'], t['col5'] ]).T\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is multi-dimensional data, which is tough to visualize in terms of the raw data and the clusters. Let's do some quick dimensional reduction-- we'll look at this in detail next week."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Project onto 2 axes with PCA\n",
    "from sklearn.decomposition import PCA\n",
    "pca = PCA(n_components=2) # 2 components\n",
    "pca.fit(X) # Do the fitting\n",
    "\n",
    "X_reduced = pca.transform(X)\n",
    "\n",
    "plt.scatter(X_reduced[:,0], X_reduced[:,1], \n",
    "            marker=\".\", color='C0', \n",
    "            alpha=0.01, edgecolors='None')\n",
    "plt.xlabel('Eigenvalue 1')\n",
    "plt.ylabel('Eigenvalue 2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above image shows the data projected into **two linear combinations that encompass the directions of greatest variance**. We'll do $K$-means clustering in these dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the KMeans clustering\n",
    "n_clusters = 6\n",
    "scaler = preprocessing.StandardScaler()\n",
    "clf = KMeans(n_clusters)\n",
    "clf.fit(scaler.fit_transform(X_reduced))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make some plots\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "ax = fig.add_subplot()\n",
    "\n",
    "# Compute a 2D histogram  of the input\n",
    "H, xedges, yedges = np.histogram2d(X_reduced[:,0], X_reduced[:,1], 50)\n",
    "\n",
    "# plot density\n",
    "ax.imshow(H.T, origin='lower', interpolation='nearest', aspect='auto',\n",
    "          extent=[xedges[0], xedges[-1],\n",
    "                  yedges[0], yedges[-1]],\n",
    "          cmap='Blues')\n",
    "\n",
    "# plot cluster centers\n",
    "cluster_centers = scaler.inverse_transform(clf.cluster_centers_)\n",
    "ax.scatter(cluster_centers[:, 0], cluster_centers[:, 1],\n",
    "           s=40, c='w', edgecolors='k')\n",
    "\n",
    "# plot cluster boundaries\n",
    "x_centers = 0.5 * (xedges[1:] + xedges[:-1])\n",
    "y_centers = 0.5 * (yedges[1:] + yedges[:-1])\n",
    "\n",
    "Xgrid = np.meshgrid(x_centers, y_centers)\n",
    "Xgrid = np.array(Xgrid).reshape((2, 50 * 50)).T\n",
    "\n",
    "H = clf.predict(scaler.transform(Xgrid)).reshape((50, 50))\n",
    "\n",
    "for i in range(n_clusters):\n",
    "    Hcp = H.copy()\n",
    "    flag = (Hcp == i)\n",
    "    Hcp[flag] = 1\n",
    "    Hcp[~flag] = 0\n",
    "\n",
    "    ax.contour(x_centers, y_centers, Hcp, [-0.5, 0.5],\n",
    "               linewidths=1, colors='k')\n",
    "\n",
    "    H = clf.predict(scaler.transform(Xgrid)).reshape((50, 50))\n",
    "    \n",
    "ax.set_xlim(xedges[0], xedges[-1])\n",
    "ax.set_ylim(yedges[0], yedges[-1])\n",
    "\n",
    "ax.set_xlabel('Eigenvalue 1')\n",
    "ax.set_ylabel('Eigenvalue 2')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above figure displays the underlying density of samples, the cluster centers, and the cluster boundaries! Very powerful!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean-shift Clustering <a class=\"anchor\" id=\"twob\"></a>\n",
    "\n",
    "**Mean-shift clustering** works by finding the modes in a kernel density estimator of the distribution. Clustering is achieved by the ***mean-shift algorithm***:\n",
    "\n",
    "1. The KDE of the dataset is computed.\n",
    "2. This allows the gradient of the distribution to be calculated. Easy to do since it's a bunch of overlapping Gaussians.\n",
    "3. Each data point is shifted in the direction of increasing gradient, which drives the points toward the modes. \n",
    "4. This process is iterated until all points have converged with clusters of other points at each of several distinct modes.\n",
    "5. Each data point is then associated with a cluster of other points.\n",
    "\n",
    "\n",
    "Let's try this using the same dataset as in $K$-means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the Mean-shift clustering\n",
    "from sklearn.cluster import MeanShift\n",
    "\n",
    "scaler = preprocessing.StandardScaler()\n",
    "bandwidth = 0.4\n",
    "#bandwidth = estimate_bandwidth(X_reduced) # this takes a long time...beware\n",
    "ms = MeanShift(bandwidth=bandwidth, \n",
    "               bin_seeding=True, \n",
    "               cluster_all=False)\n",
    "ms.fit(scaler.fit_transform(X_reduced))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels_unique = np.unique(ms.labels_)\n",
    "n_clusters = len(labels_unique[labels_unique >= 0])\n",
    "print(labels_unique)\n",
    "print(bandwidth)\n",
    "print(\"number of estimated clusters :\", n_clusters)\n",
    "\n",
    "# Make some plots\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "ax = fig.add_subplot()\n",
    "\n",
    "# Compute a 2D histogram  of the input\n",
    "H, xedges, yedges = np.histogram2d(X_reduced[:,0], X_reduced[:,1], 50)\n",
    "\n",
    "# plot density\n",
    "ax.imshow(H.T, origin='lower', interpolation='nearest', aspect='auto',\n",
    "          extent=[xedges[0], xedges[-1],\n",
    "                  yedges[0], yedges[-1]],\n",
    "          cmap='Blues')\n",
    "\n",
    "# plot cluster centers\n",
    "cluster_centers = scaler.inverse_transform(ms.cluster_centers_)\n",
    "ax.scatter(cluster_centers[:, 0], cluster_centers[:, 1],\n",
    "           s=40, c='w', edgecolors='k')\n",
    "\n",
    "# plot cluster boundaries\n",
    "x_centers = 0.5 * (xedges[1:] + xedges[:-1])\n",
    "y_centers = 0.5 * (yedges[1:] + yedges[:-1])\n",
    "\n",
    "Xgrid = np.meshgrid(x_centers, y_centers)\n",
    "Xgrid = np.array(Xgrid).reshape((2, 50 * 50)).T\n",
    "\n",
    "H = ms.predict(scaler.transform(Xgrid)).reshape((50, 50))\n",
    "\n",
    "for i in range(n_clusters):\n",
    "    Hcp = H.copy()\n",
    "    flag = (Hcp == i)\n",
    "    Hcp[flag] = 1\n",
    "    Hcp[~flag] = 0\n",
    "\n",
    "    ax.contour(x_centers, y_centers, Hcp, [-0.5, 0.5],\n",
    "               linewidths=1, colors='k')\n",
    " \n",
    "    H = ms.predict(scaler.transform(Xgrid)).reshape((50, 50))\n",
    "    \n",
    "ax.set_xlim(xedges[0], xedges[-1])\n",
    "ax.set_ylim(yedges[0], yedges[-1])\n",
    "\n",
    "ax.set_xlabel('Eigenvalue 1')\n",
    "ax.set_ylabel('Eigenvalue 2')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mean-shift will not only estimate cluster centers and boundaries, but also the number of clusters!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Summary of Chapter 6 methods from Ivezic Table 6.1. Entries correspond to (L)ow, (M)edium, (H)igh.** \n",
    "\n",
    "|Method          |Accuracy|Interpretability|Simplicity|Speed|\n",
    "|----------------|--------|----------------|----------|-----|\n",
    "|K-nearest Neighbor| H | H | H | M |\n",
    "|Kernel Density Estimation| H | H | H | H |\n",
    "|Gaussian Mixture Models| H | M | M | M |\n",
    "|Extreme Deconvolution| H | H | M | M |\n",
    "||||||\n",
    "|K-Means| L | M | H | M |\n",
    "|Max-radius minimization| L | M | M | M |\n",
    "|Mean shift| M | H | H | M |\n",
    "|Hierarchical Clustering| H | L | L | L |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlation functions <a class=\"anchor\" id=\"three\"></a>\n",
    "\n",
    "> ***Correlation functions*** *tell us how far (and on what scales) a distribution of data samples differs from a random distribution.*\n",
    "\n",
    "They have been used extensively in astrophysics, e.g., \n",
    "- examining fluctuations and structure on varying scales of the galaxy density distribution in terms of luminosity, galaxy type, age of the Universe.\n",
    "- examining the two-point correlation function of temperature fluctuations of the cosmic microwave background to unveil the composition of the Universe.\n",
    "- searching for long-timescale correlations in time-series data to examine noise in AGN lightcurves or find GW signals in pulsar-timing data.\n",
    "\n",
    "\n",
    "![](https://www.astroml.org/_images/fig_corr_diagram_1.png)\n",
    "\n",
    "\n",
    "One of the most prominent is the **two-point correlation function** which characterizes the excess probability of finding pairs of points at varying separations when compared to a random distribution. It can be described in terms of the power spectrum of fluctuations, $P(k)$ where $k=2\\pi/\\lambda$ and $\\lambda$ is the scale/wavelength of the fluctuation:\n",
    "\n",
    "$$ \\xi(r) = \\frac{1}{2\\pi^2}\\int dk\\, k^2 P(k)\\frac{\\sin(kr)}{kr}.$$\n",
    "\n",
    "This correlation function can be used to describe the density fluctuations of sources by\n",
    "\n",
    "$$ \\xi(r) = \\left\\langle \\frac{\\delta\\rho(x)}{\\rho}\\frac{\\delta\\rho(x+r)}{\\rho}\\right\\rangle,$$\n",
    "\n",
    "where $\\delta\\rho(x)/\\rho = (\\rho-\\bar\\rho)/\\rho$ is the density contrast relative to the mean $\\bar\\rho$ at position $x$.\n",
    "\n",
    "In many situations in astronomy or cosmology, the spatial correlation function or angular correlation function is modeled as a power-law, e.g. $w(\\theta) = (\\theta/\\theta_0)^\\delta$ (this is astronomy after all, where power-laws rule supreme). Angular correlation functions are often used because we care about *projected structure* on different scales, rather than depth clustering, e.g., in the CMB this kind of angular clustering is indicative of fluctuations in the primordial density field that can unveil the composition of the Universe at the time of last scattering.\n",
    "\n",
    "Higher $n$-point correlation functions can be computed (see the image above), but the two-point function is the most common. Why? When dealing with large populations of sources spread across the Universe, we can often invoke the central limit theorem to describe the statistical distribution of sources-- hence **correlations can be approximated as obeying Gaussian statistics, which are entirely defined by the mean and two-point correlators (i.e. the variance)**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell to implement AstroML's two-point correlation function on some random data\n",
    "from astroML.correlation import two_point_angular\n",
    "\n",
    "RA = 40 * np.random.random(1000)\n",
    "DEC = 10 * np.random.random(1000) # ra and dec in degrees\n",
    "\n",
    "bins = np.linspace(0.1, 10.0, 11) # edges for the 10 bins to evaluate\n",
    "corr = two_point_angular(RA, DEC, bins, method='landy-szalay')\n",
    "print(corr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now try the example from Figure 6.17 of the textbook. This computes the angular correlation function for a subset of the SDDS spectroscopic galaxy sample in the range $0.08<z<0.12$. This may take a minute or so to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "\n",
    "from __future__ import print_function, division\n",
    "\n",
    "from astroML.utils.decorators import pickle_results\n",
    "from astroML.datasets import fetch_sdss_specgals\n",
    "from astroML.correlation import bootstrap_two_point_angular\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Get data and do some quality cuts\n",
    "data = fetch_sdss_specgals()\n",
    "m_max = 17.7\n",
    "\n",
    "# redshift and magnitude cuts\n",
    "data = data[data['z'] > 0.08]\n",
    "data = data[data['z'] < 0.12]\n",
    "data = data[data['petroMag_r'] < m_max]\n",
    "\n",
    "# RA/DEC cuts\n",
    "RAmin, RAmax = 140, 220\n",
    "DECmin, DECmax = 5, 45\n",
    "data = data[data['ra'] < RAmax]\n",
    "data = data[data['ra'] > RAmin]\n",
    "data = data[data['dec'] < DECmax]\n",
    "data = data[data['dec'] > DECmin]\n",
    "\n",
    "ur = data['modelMag_u'] - data['modelMag_r']\n",
    "flag_red = (ur > 2.22)\n",
    "flag_blue = ~flag_red\n",
    "\n",
    "data_red = data[flag_red]\n",
    "data_blue = data[flag_blue]\n",
    "\n",
    "print(\"data size:\")\n",
    "print(\"  red gals: \", len(data_red))\n",
    "print(\"  blue gals:\", len(data_blue))\n",
    "\n",
    "\n",
    "#-----------------------------------------------------------\n",
    "# Set up correlation function computation\n",
    "#  This calculation takes a long time with the bootstrap resampling,\n",
    "#  so we'll save the results.\n",
    "@pickle_results(\"correlation_functions.pkl\")\n",
    "def compute_results(Nbins=16, Nbootstraps=10,  method='landy-szalay', rseed=0):\n",
    "    np.random.seed(rseed)\n",
    "    bins = 10 ** np.linspace(np.log10(1 / 60.), np.log10(6), 16)\n",
    "\n",
    "    results = [bins]\n",
    "    for D in [data_red, data_blue]:\n",
    "        results += bootstrap_two_point_angular(D['ra'],\n",
    "                                               D['dec'],\n",
    "                                               bins=bins,\n",
    "                                               method=method,\n",
    "                                               Nbootstraps=Nbootstraps)\n",
    "\n",
    "    return results\n",
    "\n",
    "(bins, r_corr, r_corr_err, r_bootstraps,\n",
    " b_corr, b_corr_err, b_bootstraps) = compute_results()\n",
    "\n",
    "bin_centers = 0.5 * (bins[1:] + bins[:-1])\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot the results\n",
    "corr = [r_corr, b_corr]\n",
    "corr_err = [r_corr_err, b_corr_err]\n",
    "bootstraps = [r_bootstraps, b_bootstraps]\n",
    "labels = ['$u-r > 2.22$\\n$N=%i$' % len(data_red),\n",
    "          '$u-r < 2.22$\\n$N=%i$' % len(data_blue)]\n",
    "\n",
    "fig = plt.figure(figsize=(7, 3))\n",
    "fig.subplots_adjust(bottom=0.2, top=0.9,\n",
    "                    left=0.13, right=0.95)\n",
    "\n",
    "for i in range(2):\n",
    "    ax = fig.add_subplot(121 + i, xscale='log', yscale='log')\n",
    "\n",
    "    ax.errorbar(bin_centers, corr[i], corr_err[i],\n",
    "                fmt='.k', ecolor='gray', lw=1)\n",
    "\n",
    "    t = np.array([0.01, 10])\n",
    "    ax.plot(t, 10 * (t / 0.01) ** -0.8, ':k', linewidth=1)\n",
    "\n",
    "    ax.text(0.95, 0.95, labels[i],\n",
    "            ha='right', va='top', transform=ax.transAxes)\n",
    "    ax.set_xlabel(r'$\\theta\\ (deg)$')\n",
    "    if i == 0:\n",
    "        ax.set_ylabel(r'$\\hat{w}(\\theta)$')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The left panel shows red galaxies with $u-r>2.22$ and the right panel shows blue galaxies with $u-r<2.22$. Uncertainties are derived from $10$ bootstrap samplings. There is significantly greater correlation structure on small angular scales in red galaxies than in blue galaxies. \n",
    "\n",
    "With tens of thousands of galaxies in each sample, this code has been optimized using **Scikit-Learn's ball-tree method for computing the correlation function-- further details given in Chapter 2, Section 2.52**."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:astr8070] *",
   "language": "python",
   "name": "conda-env-astr8070-py"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}