coursework/homeworks/vu_astr8070_s21_hw5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ASTR 8070: Astrostatistics\n",
    "***S. R. Taylor***\n",
    "___\n",
    "\n",
    "# Homework 5\n",
    "### Due: Friday, Mar 27th 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 `hw5_data_1.npy`. This is a (10 x 3) numpy array, containing some time-series measurements. It could be flux measurements from a star or some other transient. The first column is time (arbitrary units), the second column is the flux (arbitrary units), and the third column are homoescedastic measurement uncertainties on the flux. Plot this data with the uncertainties.\n",
    "\n",
    "2) One possible model for this data is a sharp burst with an exponential tail rising above a nominal background level. This model can be described as:\n",
    "\n",
    "$y = b$, for $t<t_0$\n",
    "\n",
    "$y = b + A\\exp[-\\alpha(t-t_0)]$, for $t\\geq t_0$\n",
    "\n",
    "where $b$ is the background amplitude, $A$ is the burst amplitude, $t_0$ is the burst epoch, and $\\alpha$ encodes how rapidly the burst dies off.\n",
    "\n",
    "Code this burst model, and code a likelihood function for this data with the burst model. *Hint: This is similar in spirit to the likelihood when fitting polynomials to data, except now $x$ is time, and the polynomial model is now the burst model.* \n",
    "\n",
    "3) Using one of the MCMC codes that we looked at in class (`emcee`, `pymc3`, `PTMCMCSampler`) perform a Bayesian analysis that recovers an MCMC chain of posterior probability samples for the joint parameter space of $\\{b,A,t_0,\\alpha\\}$. You must decide on how many samples are appropriate for running your MCMC. For your priors, make them $b\\in \\mathrm{Uniform}[0,50]$, $A\\in \\mathrm{Uniform}[0,50]$, $t_0\\in \\mathrm{Uniform}[0,100]$, and $\\ln\\alpha\\in \\mathrm{Uniform}[-5,5]$.\n",
    "\n",
    "4) Perform an appropriate check to burn-in your chain, and thin the chain by selecting only points every auto-correlation length. If you are left with a small number of samples after thinning then you may need to run your MCMC for longer in (3).\n",
    "\n",
    "5) Make a corner plot of the $4$-D parameter space, with labels, $[68,95]\\%$ credible levels and quantiles, and titles.\n",
    "\n",
    "6) Now choose $100$ samples randomly from your burned, thinned chain and plot the $100$ corresponding burst time-series on a figure along with the original data. You may want to choose a transparency for these $100$ burst curves that allows you to see everything nicely. These are showing the posterior spread of model solutions that are consistent with the data.\n",
    "\n",
    "7) Another possible model that may explain the data is a Gaussian profile described by:\n",
    "\n",
    "$y = b + A\\exp[-(t-t_0)^2 / 2\\sigma_W^2]$\n",
    "\n",
    "The priors for this model are the same as above for the equivalent parameters, and with $\\ln\\sigma_W\\in\\mathrm{Uniform}[-2,2]$. Following the examples from `Lecture 12` and from relevant web tutorials/documentation, use one of the nested sampling packages (`dynesty` or `UltraNest`) to compute the Bayesian evidence for this Gaussian model, and for the original burst model. Compute the Bayes factor, and state which model is the favored description of the data. *Hint: Be careful with how you transform your parameter space into a unit hypercube for the nested sampling priors. An appropriate check on whether you have done things correctly is that the posterior samples you get from the burst analysis here should give an equivalent corner plot to the MCMC you did above.*"
   ]
  },
  {
   "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
}