coursework/homeworks/vu_astr8070_s21_hw7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ASTR 8070: Astrostatistics\n",
    "***S. R. Taylor***\n",
    "___\n",
    "\n",
    "# Homework 7\n",
    "### Due: Saturday, Apr 10th at 11.59pm CST\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Only one problem this week\n",
    "\n",
    "This is a regression problem. Note of warning: the data I have provided have associated uncertainties, `yerr`. `sklearn` and `astroML` have different formats for using uncertainties; in `sklearn` fitting you can provide `sample_weight=1/yerr**2`, while in `astroML` fitting you can directly provide `y_error=yerr`.\n",
    "\n",
    "1. Read in `hw7_data_1.npy`. The dataset should consist of $500$ samples. The first two columns are the two data features, `X`. The third column is the target labels, `y`, and the final column are the heteroscedastic uncertainties on the labels, `yerr`. Make a 2-panel scatter plot of the labels versus each feature.\n",
    "\n",
    "\n",
    "2. Use multivariate linear regression on this data with its uncertainties (`astroML` rather than plain `sklearn` is the better option here). Print the best-fit intercept and slopes for the features. Make a prediction `Xgrid` over $50$ values from the min to the max value of each feature, and plot the best-fit model on a duplicate of the $2$-panel plot from (1); does this look like a good description of the data?\n",
    "\n",
    "\n",
    "3. Make a new two-dimensional prediction meshgrid `Xgrid` over the $2$-D feature space, use your best-fit model to predict the labels, and make a $3$-D interactive plot showing the data and your best-fit over the entire space. You may find the following code useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Makes the plot interactive. \n",
    "# You can pan the image.\n",
    "%matplotlib notebook  \n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "ax.plot(___, ___, y, color='k', \n",
    "        zorder=15, linestyle='none', \n",
    "        marker='o', alpha=0.5)\n",
    "ax.scatter(___.flatten(), ___.flatten(), ___, \n",
    "           facecolor=(0,0,0,0), s=20, \n",
    "           edgecolor='#70b3f0')\n",
    "ax.set_xlabel('X[0]', fontsize=12)\n",
    "ax.set_ylabel('X[1]', fontsize=12)\n",
    "ax.set_zlabel('Target', fontsize=12)\n",
    "\n",
    "ax.view_init(elev=28, azim=120)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To turn off interactive plotting\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Sticking with a linear model, now attempt a `GridSearchCV` procedure on the data to find the best `Ridge` regression hyper-parameters. Use the following `param_grid`:\n",
    "```\n",
    "ridge_params = {'alpha': [0.05,0.1,0.2,0.5,1.0], 'solver': ['svd', 'lsqr']}\n",
    "```\n",
    "What are the `best_params_` and the `best_score_` from this? Fit the data with the best options, and visualize interactively. Try `RidgeCV` also; what cross-validated `alpha_` does it return, and does this agree with `GridSeachCV`? *(Hint: remember to revert to inline matplotlib when you're done each time you make an interactive plot.)*\n",
    "\n",
    "\n",
    "5. Perform Nadaraya-Watson regression on the data, using a Gaussian kernel with a bandwidth of $0.2$. Visualize interactively as before. How does this look compared to previous models?\n",
    "\n",
    "\n",
    "6. Now model the data using Gaussian process regression and the default kernel. Check what the fit looks like with and with and without uncertainties accounted for. *(Hint: GPR in `sklearn` takes an `alpha` parameter equal to `(yerr/y)**2`.)*\n",
    "\n",
    "\n",
    "7. Let's look finally at polynomial fitting. Fit a $4$th degree polynomial to the data. Print the best-fit regression coefficients, predict on the 2D meshgrid, and visualize interactively as before. How does this look compared to the linear model and GPR models?\n",
    "\n",
    "\n",
    "8. Sticking with polynomial fitting, let's find the optimal hyper-parameters through cross-validation. \n",
    "- Using `sklearn.model_selection.train_test_split`, first do a train-test split with a `test_size` of $20\\%$. Then split out a $20\\%$ cross-validation set from the training set (leaving a somewhat smaller set to use for training). Note that you'll need to simultaneously split out `yerr` as well by just providing this as another array to `train_test_split`.\n",
    "- Compute the training and cross-validation rms errors (and from those, the BIC values) as a function of polynomial degree from $1$ to $10$ inclusive. Plot the rms errors and BIC values as a function of the polynomial degree. \n",
    "- Choose the optimal polynomial degree, and train on all training data from the first `train_test_split` operation. Predict the labels of the held-out test set, and finally compute the testing rms error."
   ]
  }
 ],
 "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": 2
}