05-Radial-Velocity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Radial Velocity Planet Search"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One area in astronomy where Bayesian approaches have been applied with great success is in the hunt for extrasolar planets.\n",
    "There are several benefits of Bayesian modeling in this situation. In particular:\n",
    "\n",
    "1. Bayesian modeling can account for the many nuisance parameters typical in a planet search. For example, when seeking evidence of a planet around a star, we need our model to account for things like the phase and the longitude of periastron, but we don't necessarily care about these parameters in the end.\n",
    "\n",
    "2. Often we have very important prior information – for example, we might have a very good constraint on the period from eclipse data, and use this to model orbital parameters with radial velocity data.\n",
    "\n",
    "3. The forward-modeling aspect of Bayesian approaches can be advantageous when dealing with detectors that have strong systematic uncertainties. For example, this idea is key to some of the recent analysis of K2 data.\n",
    "\n",
    "Here we'll take a look at a Bayesian approach to determining orbital parameters of a planet from radial velocity (RV) measurements.\n",
    "We'll start with some generated data in which we know the correct answer, and then take a look at some real RV measurements from the literature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preliminaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, we start with some imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn; seaborn.set() #nice plot formatting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Radial Velocity Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first important step is to define a mathematical (and computational) model of how the parameters of interest are reflected in our observations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some references relating to what we're going to compute below:\n",
    "\n",
    "- Balan 2009: http://adsabs.harvard.edu/abs/2009MNRAS.394.1936B\n",
    "- Exofit Manual: http://www.star.ucl.ac.uk/~lahav/ExoFitv2.pdf\n",
    "- Hou 2014: http://arxiv.org/pdf/1401.6128.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation for radial velocity is this:\n",
    "\n",
    "$$\n",
    "v(t) = V - K[ \\sin(f + \\omega) + e \\sin(\\omega)]\n",
    "$$\n",
    "\n",
    "where $V$ is the overall velocity of the system, and\n",
    "\n",
    "$$\n",
    "K = \\frac{m_p}{m_s + m_p} \\frac{2\\pi}{T}\\frac{a \\sin i}{\\sqrt{1 - e^2}}\n",
    "$$\n",
    "\n",
    "The true anomaly $f$ satisfies\n",
    "\n",
    "$$\n",
    "\\cos(f) = \\frac{\\cos(E) - e}{1 - e\\cos E}\n",
    "$$\n",
    "\n",
    "Rearranging this we can write\n",
    "$$\n",
    "f = 2 \\cdot{\\rm atan2}\\left(\\sqrt{1 + e}\\sin(E/2), \\sqrt{1 - e} \\cos(E/2)\\right)\n",
    "$$\n",
    "\n",
    "The eccentric anomaly $E$ satisfies\n",
    "$$\n",
    "M = E - e\\sin E\n",
    "$$\n",
    "\n",
    "and the mean anomaly is\n",
    "$$\n",
    "M = \\frac{2\\pi}{T}(t + \\tau)\n",
    "$$\n",
    "\n",
    "and $\\tau$ is the time of pericenter passage, which we'll parametrize with the parameter $\\chi = \\tau /  T$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the parameters needed to compute the radial velocity:\n",
    "\n",
    "- $T$: orbital period\n",
    "- $K$: amplitude of RV oscillation\n",
    "- $V$: secular offset of RV oscillation\n",
    "- $e$: eccentricity\n",
    "- $\\omega$: longitude of periastron\n",
    "- $\\chi$: dimensionless phase offset\n",
    "\n",
    "Additionally, we will fit a scatter parameter $s$ which accounts for global data errors not reflected in the reported uncertainties (this is very similar to the third parameter from the linear fit we saw earlier)\n",
    "\n",
    "For convenience, we'll store these parameters in a ``namedtuple``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "params = namedtuple('params', ['T', 'e', 'K', 'V', 'omega', 'chi', 's'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a function to compute the observed radial velocity as a function of these parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import optimize\n",
    "\n",
    "@np.vectorize\n",
    "def compute_E(M, e):\n",
    "    \"\"\"Solve Kepler's eqns for eccentric anomaly given mean anomaly\"\"\"\n",
    "    f = lambda E, M=M, e=e: E - e * np.sin(E) - M\n",
    "    return optimize.brentq(f, 0, 2 * np.pi)\n",
    "\n",
    "\n",
    "def radial_velocity(t, theta):\n",
    "    \"\"\"Compute radial velocity given orbital parameters\"\"\"\n",
    "    T, e, K, V, omega, chi = theta[:6]\n",
    "    \n",
    "    # compute mean anomaly (0 <= M < 2pi)\n",
    "    M = 2 * np.pi * ((t / T + chi) % 1)\n",
    "    \n",
    "    # solve for eccentric anomaly\n",
    "    E = compute_E(M, e)\n",
    "    \n",
    "    # compute true anomaly\n",
    "    f = 2 * np.arctan2(np.sqrt(1 + e) * np.sin(E / 2),\n",
    "                       np.sqrt(1 - e) * np.cos(E / 2))\n",
    "    \n",
    "    # compute radial velocity\n",
    "    return V - K * (np.sin(f + omega) + e * np.sin(omega))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to get a sense of whether we've done this correctly, let's use IPython's interactive features to see how the parameters change the observed RV curve (you may have to first ``pip install ipywidgets``)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from ipywidgets import interact\n",
    "\n",
    "def plot_RV(T, e, K, V, omega, chi):\n",
    "    t = np.linspace(0, 5, 200)\n",
    "    theta = [T, e, K, V, omega, chi]\n",
    "    plt.plot(t, radial_velocity(t, theta))\n",
    "    \n",
    "interact(plot_RV,\n",
    "         T=(0, 5.), K=(0, 2000.), V=(-2000., 2000.),\n",
    "         e=(0, 1.), omega=(0, 2 * np.pi), chi=(0, 1.));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model seems to be working as expected!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting Simulated Data\n",
    "\n",
    "Now let's generate some simulated data so that we can explore a Bayesian approach to modeling the radial velocity effects of a planet orbiting a star.\n",
    "We'll choose some reasonable parameters and create some data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "theta_sim = params(T=700, e=0.38, K=60, V=12,\n",
    "                   omega=3.10, chi=0.67, s=1)\n",
    "Nobs = 50\n",
    "\n",
    "rng = np.random.RandomState(0)\n",
    "t_sim = 1400 + 600 * rng.rand(Nobs)\n",
    "err_sim = 5 + 5 * rng.rand(Nobs)\n",
    "rv_sim = radial_velocity(t_sim, theta_sim) + err_sim * rng.randn(Nobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.errorbar(t_sim, rv_sim, err_sim, fmt='.k');\n",
    "xlim = plt.xlim()\n",
    "t_fit = np.linspace(xlim[0], xlim[1], 500)\n",
    "plt.plot(t_fit, radial_velocity(t_fit, theta_sim), color='gray', lw=8, alpha=0.2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian Fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next let's use ``emcee`` to do a Bayesian model fit to this data.\n",
    "We'll follow some of the references above and use a flat prior on most parameters, and a Jeffreys Prior on the scale factors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "theta_lim = params(T=(0.2, 2000),\n",
    "                   e=(0, 1),\n",
    "                   K=(0.01, 2000),\n",
    "                   V=(-2000, 2000),\n",
    "                   omega=(0, 2 * np.pi),\n",
    "                   chi=(0, 1),\n",
    "                   s=(0.001, 100))\n",
    "theta_min, theta_max = map(np.array, zip(*theta_lim))\n",
    "\n",
    "def log_prior(theta):\n",
    "    if np.any(theta < theta_min) or np.any(theta > theta_max):\n",
    "        return -np.inf # log(0)\n",
    "    \n",
    "    # Jeffreys Prior on T, K, and s\n",
    "    return -np.sum(np.log(theta[[0, 2, 6]]))\n",
    "\n",
    "def log_likelihood(theta, t, rv, rv_err):\n",
    "    sq_err = rv_err ** 2 + theta[6] ** 2\n",
    "    rv_model = radial_velocity(t, theta)\n",
    "    return -0.5 * np.sum(np.log(sq_err) + (rv - rv_model) ** 2 / sq_err)\n",
    "\n",
    "def log_posterior(theta, t, rv, rv_err):\n",
    "    ln_prior = log_prior(theta)\n",
    "    if np.isinf(ln_prior):\n",
    "        return ln_prior\n",
    "    else:\n",
    "        return ln_prior + log_likelihood(theta, t, rv, rv_err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the likelihood model in place, we can now use ``emcee`` to sample the posterior and view the resulting chains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import emcee\n",
    "\n",
    "ndim = len(theta_sim)  # number of parameters in the model\n",
    "nwalkers = 50  # number of MCMC walkers\n",
    "\n",
    "# start with theta near the midpoint of the prior range\n",
    "rng = np.random.RandomState(42)\n",
    "theta_guess = 0.5 * (theta_min + theta_max)\n",
    "theta_range = (theta_max - theta_min)\n",
    "starting_guesses = theta_guess + 0.05 * theta_range * rng.randn(nwalkers, ndim)\n",
    "\n",
    "# create the sampler; fix the random state for replicability\n",
    "sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=(t_sim, rv_sim, err_sim))\n",
    "sampler.random_state = rng\n",
    "\n",
    "# time and run the MCMC\n",
    "%time pos, prob, state = sampler.run_mcmc(starting_guesses, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_chains(sampler):\n",
    "    fig, ax = plt.subplots(7, figsize=(8, 10), sharex=True)\n",
    "    for i in range(7):\n",
    "        ax[i].plot(sampler.chain[:, :, i].T, '-k', alpha=0.2);\n",
    "        ax[i].set_ylabel(params._fields[i])\n",
    "\n",
    "plot_chains(sampler)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yikes! it's all over the place!\n",
    "\n",
    "The issue here is that our initialization was haphazard and the posterior is extremely multimodal (especially in *T*); given a number of steps approaching infinity, the MCMC algorithm would converge, but we don't have an infinite amount of time to wait! Instead we can more carefully initialize the walkers.\n",
    "\n",
    "First, let's use a Lomb-Scargle periodogram to find a suitable guess at the period.\n",
    "The [gatspy](http://www.astroml.org/gatspy/) package has a nice implementation (first you need to ``pip install gatspy``)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from gatspy.periodic import LombScargleFast\n",
    "model = LombScargleFast()\n",
    "\n",
    "model.fit(t_sim, rv_sim, err_sim)\n",
    "periods, power = model.periodogram_auto()\n",
    "plt.semilogx(periods, power)\n",
    "plt.xlim(0, 10000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can choose a sensible starting point with this period, and with other parameters estimated from the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def make_starting_guess(t, rv, rv_err):\n",
    "    model = LombScargleFast()\n",
    "    model.optimizer.set(period_range=theta_lim.T,\n",
    "                        quiet=True)\n",
    "    model.fit(t, rv, rv_err)\n",
    "\n",
    "    rv_range = 0.5 * (np.max(rv) - np.min(rv))\n",
    "    rv_center = np.mean(rv)\n",
    "    return params(T=model.best_period,\n",
    "                  e=0.1,\n",
    "                  K=rv_range,\n",
    "                  V=rv_center,\n",
    "                  omega=np.pi,\n",
    "                  chi=0.5,\n",
    "                  s=rv_err.mean())\n",
    "\n",
    "theta_guess = make_starting_guess(t_sim, rv_sim, err_sim)\n",
    "theta_guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sampler.reset()\n",
    "start = theta_guess * (1 + 0.01 * rng.randn(nwalkers, ndim))\n",
    "pos, prob, state = sampler.run_mcmc(start, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_chains(sampler)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much more reasonable! The trace appears to have stabilized by the end of this, so let's reset and get a clean 1000 samples from the posterior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sampler.reset()\n",
    "pos, prob, state = sampler.run_mcmc(pos, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the [corner.py](https://pypi.python.org/pypi/corner) package, we can take a look at this multi-dimensional posterior, along with the input values for the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import corner\n",
    "corner.corner(sampler.flatchain, labels=params._fields, truths=theta_sim);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can view the fit to our data by sampling from the posterior, and computing the model associated with each sample: because the trace already is a sample from the posterior, we can simply draw randomly from these points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rng = np.random.RandomState(42)\n",
    "\n",
    "rv_fit = [radial_velocity(t_fit, sampler.flatchain[i])\n",
    "          for i in rng.choice(sampler.flatchain.shape[0], 200)]\n",
    "\n",
    "plt.errorbar(t_sim, rv_sim, err_sim, fmt='.k')\n",
    "plt.plot(t_fit, np.transpose(rv_fit), '-k', alpha=0.01)\n",
    "plt.xlabel('time (days)')\n",
    "plt.ylabel('radial velocity (km/s)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can treat everything but the period and eccentricity as a nuisance parameter (i.e. marginalize over them) and take a look at our parameter constraints:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import corner\n",
    "corner.corner(sampler.flatchain[:, :2], \n",
    "              labels=params._fields[:2],\n",
    "              truths=theta_sim[:2]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the marginalized posterior distribution for the period and eccentricity (i.e. integrating over all the other parameters as nuisance parameters).\n",
    "We see that, as we might hope, the true value lies well withing the uncertainty implied by the marginalized posterior!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Breakout: 47 Ursae Majoris\n",
    "\n",
    "Your task now is to repeat this analysis using some real data, which we'll take from table 1 of [Fischer et al 2002](http://iopscience.iop.org/article/10.1086/324336/meta).\n",
    "\n",
    "If you're curious about one possible solution to this problem, see [Solutions-05](Solutions-05.ipynb). As usual, try to fight the temptation to peek at this until after you've given the problem a reasonable effort!\n",
    "\n",
    "The following IPython magic command will create the data file ``47UrsaeMajoris.txt`` in the current directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%file 47UrsaeMajoris.txt\n",
    "date      rv     rv_err\n",
    "6959.737  -60.48 14.00\n",
    "7194.912  -53.60 7.49\n",
    "7223.798  -38.36 6.14\n",
    "7964.893  0.60 8.19\n",
    "8017.730  -28.29 10.57\n",
    "8374.771  -40.25 9.37\n",
    "8647.897  42.37 11.41\n",
    "8648.910  32.64 11.02\n",
    "8670.878  55.45 11.45\n",
    "8745.691  51.78 8.76\n",
    "8992.061  4.49 11.21\n",
    "9067.771  -14.63 7.00\n",
    "9096.734  -26.06 6.79\n",
    "9122.691  -47.38 7.91\n",
    "9172.686  -38.22 10.55\n",
    "9349.912  -52.21 9.52\n",
    "9374.964  -48.69 8.67\n",
    "9411.839  -36.01 12.81\n",
    "9481.720  -52.46 13.40\n",
    "9767.918  38.58 5.48\n",
    "9768.908  36.68 5.02\n",
    "9802.789  37.93 3.85\n",
    "10058.079  15.82 3.45\n",
    "10068.980  15.46 4.63\n",
    "10072.012  21.20 4.09\n",
    "10088.994  1.30 4.25\n",
    "10089.947  6.12 3.70\n",
    "10091.900  0.00 4.16\n",
    "10120.918  4.07 4.16\n",
    "10124.905  0.29 3.74\n",
    "10125.823  -1.87 3.79\n",
    "10127.898  -0.68 4.10\n",
    "10144.877  -4.13 5.26\n",
    "10150.797  -8.14 4.18\n",
    "10172.829  -10.79 4.43\n",
    "10173.762  -9.33 5.43\n",
    "10181.742  -23.87 3.28\n",
    "10187.740  -16.70 4.67\n",
    "10199.730  -16.29 3.98\n",
    "10203.733  -21.84 4.92\n",
    "10214.731  -24.51 3.67\n",
    "10422.018  -56.63 4.23\n",
    "10438.001  -39.61 3.91\n",
    "10442.027  -44.62 4.05\n",
    "10502.853  -32.05 4.69\n",
    "10504.859  -39.08 4.65\n",
    "10536.845  -22.46 5.18\n",
    "10537.842  -22.83 4.16\n",
    "10563.673  -17.47 4.03\n",
    "10579.697  -11.01 3.84\n",
    "10610.719  -8.67 3.52\n",
    "10793.957  37.00 3.78\n",
    "10795.039  41.85 4.80\n",
    "10978.684  36.42 5.01\n",
    "11131.066  13.56 6.61\n",
    "11175.027  -3.74 8.17\n",
    "11242.842  -21.85 5.43\n",
    "11303.712  -48.75 4.63\n",
    "11508.070  -51.65 8.37\n",
    "11536.064  -72.44 4.73\n",
    "11540.999  -57.58 5.97\n",
    "11607.916  -43.94 4.94\n",
    "11626.771  -39.14 7.03\n",
    "11627.754  -50.88 6.21\n",
    "11628.727  -51.52 5.87\n",
    "11629.832  -51.86 4.60\n",
    "11700.693  -24.58 5.20\n",
    "11861.049  14.64 5.33\n",
    "11874.068  14.15 5.75\n",
    "11881.045  18.02 4.15\n",
    "11895.068  16.96 4.60\n",
    "11906.014  11.73 4.07\n",
    "11907.011  22.83 4.38\n",
    "11909.042  23.42 3.78\n",
    "11910.955  18.34 4.33\n",
    "11914.067  15.45 5.37\n",
    "11915.048  24.05 3.82\n",
    "11916.033  23.16 3.67\n",
    "11939.969  27.53 5.08\n",
    "11946.960  21.44 4.18\n",
    "11969.902  30.99 4.58\n",
    "11971.894  38.36 5.01\n",
    "11998.779  33.82 3.93\n",
    "11999.820  27.52 3.98\n",
    "12000.858  23.40 4.07\n",
    "12028.740  37.08 4.95\n",
    "12033.746  26.28 5.24\n",
    "12040.759  31.12 3.54\n",
    "12041.719  34.04 3.45\n",
    "12042.695  31.38 3.98\n",
    "12073.723  21.81 4.73"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An easy way to load data in this form is the [pandas](http://pandas.pydata.org) package, which implements a DataFrame object (basically, a labeled data table).\n",
    "Reading the CSV file is a one-line operation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "data = pd.read_csv('47UrsaeMajoris.txt', delim_whitespace=True)\n",
    "t, rv, rv_err = data.values.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this data in hand, you can now start to explore it and search for a planet in the radial wobbles of the star."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize the Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Fill-in code to visualize the data\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute the Periodogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Compute the periodogram to look for significant periodicity\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialize and run the MCMC Chain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Fill-in your code here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the chains: have they stabilized?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Fill-in your code here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### If necessary, reset and re-run the sampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Fill-in your code here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make a Corner Plot for the Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Fill-in your code here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the Model Fits over the Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Fill-in your code here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results: Report your (joint) uncertainties on period and eccentricity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Fill-in your code here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extra Credit\n",
    "\n",
    "If you finish early, try tackling this...\n",
    "\n",
    "The source of the above data is a paper which actually reports *two* detected planets. Build a Bayesian model which models both of them at once: can you find signals from both planets in the data?"
   ]
  }
 ],
 "metadata": {
  "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}