coursework/homeworks/vu_astr8070_s21_hw6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ASTR 8070: Astrostatistics\n",
    "***S. R. Taylor***\n",
    "___\n",
    "\n",
    "# Homework 6\n",
    "### Due: Saturday, Apr 3rd at 11.59pm CST\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1\n",
    "\n",
    "This problem has to do with density estimation and clustering. \n",
    "\n",
    "1. Read in `hw6_data_1.npy`. The dataset should consist of $1000$ samples with $2$ features. Note that the last column are the labels, not another feature. Make a $2$-D scatter plot of the data colored by their labels. \n",
    "\n",
    "\n",
    "2. Start with kernel density estimation on the data set. First, perform a grid search $5$-fold cross-validation to find the best bandwidth, testing $30$ bandwidths between $0.1$ and $1.0$. What is this best bandwidth? Using this best bandwidth, plot the two-dimensional kernel density estimate of the distribution. \n",
    "\n",
    "\n",
    "3. Now try a nearest neighbors approach to estimating the density. Use the Bayesian nearest neighbors option. Plot the estimated distribution. What value of $k$ neighbors do you need to make the plot look similar to your KDE distribution?\n",
    "\n",
    "\n",
    "4. Now fit a Gaussian mixture model. Compute the AIC and BIC for the number of GMM components between $1$ and $10$. Plot these AIC and BIC values as a function of number of components. What is the optimal number according to the BIC? Fit the data set with this optimal number of GMM components, and create a new plot that shows the data set, the GMM component mean locations, and the $2$-sigma GMM component ellipses. *(Use the code below to help with plotting the ellipses.)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Kludge to fix the bug with draw_ellipse in astroML v1.0\n",
    "from matplotlib.patches import Ellipse\n",
    "\n",
    "def draw_ellipse(mu, C, sigmas=[1, 2, 3], ax=None, **kwargs):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "\n",
    "    # find principal components and rotation angle of ellipse\n",
    "    sigma_x2 = C[0, 0]\n",
    "    sigma_y2 = C[1, 1]\n",
    "    sigma_xy = C[0, 1]\n",
    "\n",
    "    alpha = 0.5 * np.arctan2(2 * sigma_xy,\n",
    "                             (sigma_x2 - sigma_y2))\n",
    "    tmp1 = 0.5 * (sigma_x2 + sigma_y2)\n",
    "    tmp2 = np.sqrt(0.25 * (sigma_x2 - sigma_y2) ** 2 \n",
    "                   + sigma_xy ** 2)\n",
    "\n",
    "    sigma1 = np.sqrt(tmp1 + tmp2)\n",
    "    sigma2 = np.sqrt(tmp1 - tmp2)\n",
    "\n",
    "    for scale in scales:\n",
    "        ax.add_patch(Ellipse((mu[0], mu[1]),\n",
    "                             2 * scale * sigma1, 2 * scale * sigma2,\n",
    "                             alpha * 180. / np.pi,\n",
    "                             **kwargs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Now try Kmeans clustering with $5$ clusters. \n",
    "    - Use the `StandardScalar` to scale the data.\n",
    "    - Fit the Kmeans model.\n",
    "    - Find the cluster centers, and then inverse transform thse cluster centers back to the original unscaled coordinates.\n",
    "    - Predict the labels for the scaled data.\n",
    "    - Finally, make a two-panel side by side plot showing (a) the original data colored by its labels, then (b) the data colored by the Kmeans predicted labels, with cluster centers superposed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Problem 2\n",
    "\n",
    "This problem has to do with dimensional reduction. We're going to load in a sample of SDSS Imaging data. \n",
    "\n",
    "1. Execute the cell below to read in the data, print out the feature names, and create a data matrix out of a subset of features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('ra', 'dec', 'run', 'rExtSFD', 'uRaw', 'gRaw', 'rRaw', 'iRaw', 'zRaw', 'uErr', 'gErr', 'rErr', 'iErr', 'zErr', 'uRawPSF', 'gRawPSF', 'rRawPSF', 'iRawPSF', 'zRawPSF', 'upsfErr', 'gpsfErr', 'rpsfErr', 'ipsfErr', 'zpsfErr', 'type', 'ISOLATED')\n"
     ]
    }
   ],
   "source": [
    "from astroML.datasets import fetch_imaging_sample\n",
    "data = fetch_imaging_sample()  \n",
    "data.shape  # number of objects in dataset\n",
    "\n",
    "print(data.dtype.names)\n",
    "\n",
    "keylist = ['ra', 'dec', 'rExtSFD', 'uRaw', \n",
    "           'gRaw', 'rRaw', 'iRaw', 'zRaw', \n",
    "           'uErr', 'gErr', 'rErr', 'iErr', \n",
    "           'zErr', 'uRawPSF', 'gRawPSF', \n",
    "           'rRawPSF', 'iRawPSF', 'zRawPSF', \n",
    "           'upsfErr', 'gpsfErr', 'rpsfErr', \n",
    "           'ipsfErr', 'zpsfErr']\n",
    "\n",
    "X = np.column_stack([data[key] for key in keylist]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Use PCA (with randomized svd_solver for speed if necessary) to reduce the data matrix down to $2$ features. What is the explained variance of the data encapsulated in these eigen-features?\n",
    "\n",
    "\n",
    "3. Access the `type` key of the `data` structure and make an array of labels out of these. Do some research in astroML documentation to find out what these integer types correspond to, and state that here.\n",
    "\n",
    "\n",
    "4. Choose $5000$ random integers between 0 and the number of samples in the data matrix. Record these integers because you'll use them later. Make a scatter plot of the PCA-reduced data for these $5000$ random samples, colored by their corresponding `type`. *(You may want to set the transparency to be lower than 1 to see the mixing of samples.)*\n",
    "\n",
    "\n",
    "5. Now try some non-linear dimensional reduction. These algorithms are slower than PCA, so you will operate only on the $5000$ random samples identified in the previous part. \n",
    "    - Try `LocallyLinearEmbedding`, `Isomap`, and `TSNE` algorithms, setting the number of components to be $2$ in all cases. \n",
    "    - As in the PCA case, make scatter plots of the dimensionally-reduced data, color coded by their `type`. For LLE and Isomap, experiment with the number of nearest neighbors between $5$ and $100$ to see what visually gives the best separation in `type` populations. For TSNE, do the same for the perplexity attribute. \n",
    "    - Which algorithm gives the cleanest way to visually see the two populations of sources? *(This will be subjective according to the samples you trained on, and even the randomness of the algorithms.)*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}