coursework/homeworks/vu_astr8070_s22_hw3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ASTR 8070: Astrostatistics\n",
    "***S. R. Taylor***\n",
    "___\n",
    "\n",
    "# Homework 3\n",
    "### Due: Saturday, Feb 12th at 11.59pm CST\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Only one problem this week\n",
    "\n",
    "This problem uses a dataset in `/coursework/homeworks/hw_data/`.\n",
    "\n",
    "1) Read in `hw3_data_1.npy`. This is a (50 x 2) numpy array, with measurements in the first column and uncertainties in the second column. Using the analytic results for heteroscedastic Gaussian data from lectures, compute the sample mean and the standard error on the sample mean from for this data.\n",
    "\n",
    "2) Reusing some approaches and tools from `Lecture_6`, write a ln-likelihood function for heteroscedastic Gaussian data, and use it in a fitting algorithm to find the best-fit mean. *Remember that scipy optimizers are set up to minimize functions.*\n",
    "\n",
    "3) Using the same numerical technique from `Lecture_5`, compute the Fisher uncertainty estimate on the mean.\n",
    "\n",
    "4) Using the bootstrap method, generate $1000$ bootstrap realizations of this dataset. *DO NOT use the `astroML` code. Write your own bootstrap function from scratch. Also recall that when resampling data, measurements and uncertainties should stay paired together.*\n",
    "\n",
    "5) Repeat (2) with all $1000$ boostrap datasets to find the distribution of the sample mean. Plot a normalized histogram of these bootstrap means, and overplot a Gaussian pdf with the mean and std found in (1). Do these agree?\n",
    "\n",
    "6) While we have fitted a heteroscedastic Gaussian to this data, let's try something else. Write some code to define a ln-likelihood for a Laplace distribution evaluated on this data. Fit simultaneously for the Laplace location parameter $\\mu$ and scale parameter $\\Delta$.\n",
    "\n",
    "7) Compute the AIC values for the heteroscedastic Gaussian model and the Laplacian model. Which model is favored by the data?\n",
    "\n",
    "8) Using the $1000$ bootstrap datasets from before, fit for the Laplacian $\\mu$ and $\\Delta$ for each. Make a nice `corner` plot of the distributions of $\\mu$ and $\\Delta$ that shows both the marginal $1$D distributions and the joint $2$D distribution. Make sure the plot has labels, shows the titles on each $1$D marginal panel, and has $68\\%$ and $95\\%$ levels.\n",
    "\n",
    "9) Let's finish with a Fisher uncertainty estimate of the Laplacian parameters. Use the following code to install `numdifftools` which provides a simple way to compute derivatives. We can then compute the Hessian matrix, which is the matrix of the second derivatives of the user's function. This should be computed at the best-fit Laplacian parameters $\\mu$ and $\\Delta$. To finish, invert the matrix, and then take the square root. The diagonal entries will then be the Fisher uncertainties on $\\mu$ and $\\Delta$. How does these compare to the bootstrap distribution widths found in (8)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install numdifftools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numdifftools as nd\n",
    "H = nd.Hessian(f_lnlaplace)([beta_laplace[0], beta_laplace[1]])\n",
    "sigma_laplace = np.linalg.inv(H)**0.5"
   ]
  },
  {
   "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
}