{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# MCMC vs OFTI Comparison\n", "\n", "by Sarah Blunt, 2018" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Welcome to the OFTI/MCMC comparison tutorial! This tutorial is meant to help you understand the differences between OFTI and MCMC algorithms so you know which one to pick for your data. \n", "\n", "Before we start, I'll give you the short answer: **for orbit fractions less than 5%, OFTI is generally faster to converge than MCMC**. This is not a hard-and-fast statistical rule, but I've found it to be a useful guideline. \n", "\n", "This tutorial is essentially an abstract of [Blunt et al (2017)](https://ui.adsabs.harvard.edu/#abs/2017AJ....153..229B/abstract). To dig deeper, I encourage you to read the paper (Sections 2.2-2.3 in particular).\n", "\n", "**Goals of This Tutorial**:\n", "- Understand qualitatively why OFTI converges faster than MCMC for certain datasets.\n", "- Learn to make educated choices of backend algorithms for your own datasets.\n", "\n", "**Prerequisites**:\n", "- This tutorial assumes knowledge of the `orbitize` API. Please go through at least the [OFTI](https://orbitize.readthedocs.io/en/latest/tutorials/OFTI_tutorial.html) and [MCMC](https://orbitize.readthedocs.io/en/latest/tutorials/MCMC_tutorial.html) introduction tutorials before this one.\n", "- This tutorial also assumes a qualitative understanding of OFTI and MCMC algorithms. I suggest you check out at least Section 2.1 of [Blunt et al (2017)](https://ui.adsabs.harvard.edu/#abs/2017AJ....153..229B/abstract) and [this blog post](https://jeremykun.com/2015/04/06/markov-chain-monte-carlo-without-all-the-bullshit/) before attempting to decode this tutorial.\n", "\n", "**Jargon**:\n", "- I will often use **orbit fraction**, or the fraction of the orbit spanned by the astrometric observations, as a figure of merit. In general, OFTI will converge faster than MCMC for small orbit fractions. \n", "- **Convergence** is defined differently for OFTI and for MCMC (see the OFTI paper for details). An OFTI run needs to accept a statistically large number of orbits for convergence, since each accepted orbit is independent of all others. For MCMC, convergence is a bit more complicated. At a high level, an MCMC run has converged when all walkers have explored the entire parameter space. There are several metrics for estimating MCMC convergence (e.g. GR statistic, min Tz statistic), but we'll just estimate convergence qualitatively in this tutorial." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import astropy.table\n", "import time\n", "\n", "np.random.seed(5)\n", "\n", "from orbitize.kepler import calc_orbit\n", "from orbitize import system, sampler\n", "from orbitize.read_input import read_file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generate Synthetic Data\n", "\n", "Let's start by defining a function to generate synthetic data. This will allow us to easily test convergence speeds for different orbit fractions. I'll include the number of observations and the noise magnitude as keywords; I encourage you to test out different values throughout the tutorial!" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "mtot = 1.2 # total system mass [M_sol]\n", "plx = 60.0 # parallax [mas]\n", "\n", "\n", "def generate_synthetic_data(sma=30.0, num_obs=4, unc=10.0):\n", " \"\"\"Generate an orbitize-table of synethic data\n", "\n", " Args:\n", " sma (float): semimajor axis (au)\n", " num_obs (int): number of observations to generate\n", " unc (float): uncertainty on all simulated RA & Dec measurements (mas)\n", "\n", " Returns:\n", " 2-tuple:\n", " - `astropy.table.Table`: data table of generated synthetic data\n", " - float: the orbit fraction of the generated data\n", " \"\"\"\n", "\n", " # assumed ground truth for non-input orbital parameters\n", " ecc = 0.5 # eccentricity\n", " inc = np.pi / 4 # inclination [rad]\n", " argp = 0.0\n", " lan = 0.0\n", " tau = 0.8\n", "\n", " # calculate RA/Dec at three observation epochs\n", " observation_epochs = np.linspace(\n", " 51550.0, 52650.0, num_obs\n", " ) # `num_obs` epochs between ~2000 and ~2003 [MJD]\n", " num_obs = len(observation_epochs)\n", " ra, dec, _ = calc_orbit(\n", " observation_epochs, sma, ecc, inc, argp, lan, tau, plx, mtot\n", " )\n", "\n", " # add Gaussian noise to simulate measurement\n", " ra += np.random.normal(scale=unc, size=num_obs)\n", " dec += np.random.normal(scale=unc, size=num_obs)\n", "\n", " # define observational uncertainties\n", " ra_err = dec_err = np.ones(num_obs) * unc\n", "\n", " # make a plot of the data\n", " plt.figure()\n", " plt.errorbar(ra, dec, xerr=ra_err, yerr=dec_err, linestyle=\"\")\n", " plt.xlabel(\"$\\\\Delta$ RA\")\n", " plt.ylabel(\"$\\\\Delta$ Dec\")\n", "\n", " # calculate the orbital fraction\n", " period = np.sqrt((sma**3) / mtot)\n", " orbit_coverage = (\n", " max(observation_epochs) - min(observation_epochs)\n", " ) / 365.25 # [yr]\n", " orbit_fraction = 100 * orbit_coverage / period\n", "\n", " data_table = astropy.table.Table(\n", " [observation_epochs, [1] * num_obs, ra, ra_err, dec, dec_err],\n", " names=(\"epoch\", \"object\", \"raoff\", \"raoff_err\", \"decoff\", \"decoff_err\"),\n", " )\n", " # read into orbitize format\n", " data_table = read_file(data_table)\n", "\n", " return data_table, orbit_fraction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Short Orbit Fraction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's use the function above to generate some synthetic data with a short orbit fraction, and fit it with OFTI:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The orbit fraction is 2.0%\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEICAYAAAC0+DhzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVHklEQVR4nO3df7RlZX3f8ffH4ac18vOGkpnRQcEa4oqDjghFWwt1qbQJ6MIwtBHqIgtqsVF0JZF0rRqzSpckUVqSaMCCAjECQROpRSsqxrDWAGvAAWHwx1U0zDjCgPyQUmnAb/84z80c7nPnxx3vOWfG+36tddbZ+9nP3vNlc879nGefffZOVSFJ0rBnTboASdKux3CQJHUMB0lSx3CQJHUMB0lSZ49JF7AQDj744FqxYsWky5Ck3cptt932YFVNzbXsZyIcVqxYwdq1ayddhiTtVpJ8b2vLPKwkSeoYDpKkjuEgSeoYDpKkjuEgSeoYDpKkjuEgSeoYDpKkjuGwwE69eA2nXrxm0mVI0k/FcJAkdQwHSVJnbOGQZEmSryb5TJs/LMktSaaTXJ1kr9a+d5ufbstXjKtGSdLAOEcO7wDuGZq/ALiwqg4HHgbObO1nAg+39gtbP0nSGI0lHJIsA/4V8D/afIDjgWtbl8uBk9v0SW2etvyE1l+SNCbjumT3fwN+G/i5Nn8Q8EhVPdXmNwBL2/RS4D6AqnoqyaOt/4PDG0xyFnAWwPOe97ydLmyhzyxav+mxkWz36rOPXdDtSdK2jHzkkORfAw9U1W0Lud2quqSqVlXVqqmpOe9VIUnaSeMYORwH/GqSE4F9gOcC/x3YP8kebfSwDNjY+m8ElgMbkuwB7Ac8NKriFvoT+cyIwU/6knZnIx85VNV5VbWsqlYAq4EvVdW/BW4ETmndzgA+3aava/O05V+qqhp1nZKkLSb5O4ffAd6VZJrBdwqXtvZLgYNa+7uA90yoPklatMZ6D+mq+jLw5Tb9HeDoOfr8GHjzOOuSJD2Tv5CWJHUMB0lSZ6yHlRYDz1KS9LPAkYMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqTPycEiyT5Jbk9yR5O4k72vtH0tyb5J17bGytSfJRUmmk9yZ5GWjrlGS9Ex7jOHfeBI4vqoeT7IncFOSz7Zlv1VV187q/wbgiPZ4JfDh9ixJGpORjxxq4PE2u2d71DZWOQm4oq13M7B/kkNHXackaYuxfOeQZEmSdcADwA1VdUtbdH47dHRhkr1b21LgvqHVN7S22ds8K8naJGs3b948yvIladEZSzhU1dNVtRJYBhyd5CXAecCLgVcABwK/M89tXlJVq6pq1dTU1EKXLEmL2ljPVqqqR4AbgddX1aZ26OhJ4KPA0a3bRmD50GrLWpskaUzGcbbSVJL92/S+wGuBr898j5AkwMnAXW2V64DT21lLxwCPVtWmUdcpSdpiHGcrHQpcnmQJgzC6pqo+k+RLSaaAAOuAf9/6Xw+cCEwDTwBvHUONkqQhIw+HqroTOGqO9uO30r+Ac0ZdlyRp6/yFtCSpYzhot3HqxWs49eI1ky5DWhQMB0lSx3CQJHUMB0lSx3CQJHUMB0lSZxw/gtMitdBnFq3f9NiCb/fqs49dsG1JP0scOUiSOo4cNDIL/al8ZsTgp31p9Bw5SJI6hoMkqWM4SJI6hoMkqWM4SJI6nq2k3YZnKUnj48hBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJnZGHQ5J9ktya5I4kdyd5X2s/LMktSaaTXJ1kr9a+d5ufbstXjLpGSdIzjWPk8CRwfFW9FFgJvD7JMcAFwIVVdTjwMHBm638m8HBrv7D1kySN0cjDoQYeb7N7tkcBxwPXtvbLgZPb9Eltnrb8hCQZdZ2SpC3G8p1DkiVJ1gEPADcA3wYeqaqnWpcNwNI2vRS4D6AtfxQ4aI5tnpVkbZK1mzdvHvF/gSQtLmMJh6p6uqpWAsuAo4EXL8A2L6mqVVW1ampq6qfdnCRpyFjPVqqqR4AbgWOB/ZPM3E9iGbCxTW8ElgO05fsBD42zTkla7MZxttJUkv3b9L7Aa4F7GITEKa3bGcCn2/R1bZ62/EtVVaOuU5K0xTjuBHcocHmSJQzC6Jqq+kyS9cBVSf4L8FXg0tb/UuDKJNPAD4HVY6hRkjRk5OFQVXcCR83R/h0G3z/Mbv8x8OZR1yVJ2jp/IS1J6hgOkqSO4SBJ6hgOkqSO4SBJ6uxwOCS5fOb3Cm3+gCSXjaQqSdJEzWfk8MvtF84AVNXDzHGKqiRp9zefcHhWkgNmZpIcyHh+RCdJGrP5/HH/AHBzkmva/JuB8xe+JEnSpO1wOFTVFUnWMrgPA8Cbqmr9aMqSJE3SfL6QDvAy4MCq+hPg8STd5S8kSbu/+Xzn8CEGl9o+rc3/CPjTBa9IkjRx8/nO4ZVV9bIkX4XB2UpJ9hpRXZKkCZrPyOHv22W3Cwb3aQB+MpKqJEkTNZ9wuAj4K+CQJOcDNwH/dSRVSZImaj5nK308yX3AaxjctvPkqrpnVIVJkiZnu+HQzlJ6L/B2BiONAE8BBwG/P9LqJEkTsSOHlc4FjgNeUVUHVtUBwCuB45KcO9LqJEkTsSPh8BbgtKq6d6ah3eLz14HTR1WYJGlydiQc9qyqB2c3VtVmYM+FL0mSNGk7Eg7/byeXSZJ2UztyttJLkzw2R3uAfRa4HknSLmC74VBVS8ZRiCRp1+FtQiVJHcNBktQxHCRJnZ0OhySvSrLdS3YnWZ7kxiTrk9yd5B2t/feSbEyyrj1OHFrnvCTTSb6R5HU7W6MkaefM6x7QSY4C/g3wa8APgBcD52xntaeAd1fV7Ul+DrgtyQ1t2YVV9Uez/o0jgdXALwG/AHwhyYuq6un51CpJ2nnbHTkkeVGS9yb5BvAR4EHgNVX1SuCH21u/qjZV1e1t+kfAPcDSbaxyEnBVVT3ZfpU9DXjHOUkaox05rPR14ETglKpaVVUXDF1Ko+bzjyVZARwF3NKa3p7kziSXJTmgtS0F7htabQNzhEmSs5KsTbJ28+bN8ylDkrQdOxIObwLuBT6f5Mokv5Jk3pfNSPIc4JPAO6vqMeDDwAuBlcAm4APz2V5VXdLCatXU1NR8y5EkbcN2w6Gq/rqqVgOHA58FzgI2JPko8Nwd+UdamHwS+HhVfapt9/6qerqqfsLgcNXMoaONwPKh1Ze1NknSmOzw2UpV9X+q6i+q6lcYfBG9Brhze+u1+0FcCtxTVR8caj90qNsbgbva9HXA6iR7JzkMOAK4dUfrlCT99OZ1ttKMqnoYuKQ9tuc4Bpf9/lqSda3td4HTkqxk8L3Fd4Gz27bvTnINsJ7BmU7neKaSJI3XToXDfFTVTQwu0jfb9dtY53zg/JEVJUnaJn8hLUnqGA6SpI7hIEnqGA6SpI7hIEnqGA6SpI7hIEnqGA6SpI7hIEnqGA6SFq1TL17DqRevmXQZuyTDQZLUMRwkSR3DQZLUMRwkSR3DQZLUGfn9HCRpoSz0mUXrNz02ku1effaxC7q9SXDkIEnqOHKQtNtY6E/kMyOGn4VP+gvNkYMkqWM4SJI6hoMkqWM4SJI6hoMkqePZSpIWLc9S2jpHDpKkjuEgSeqMPBySLE9yY5L1Se5O8o7WfmCSG5J8qz0f0NqT5KIk00nuTPKyUdcoSXqmcYwcngLeXVVHAscA5yQ5EngP8MWqOgL4YpsHeANwRHucBXx4DDVKkoaMPByqalNV3d6mfwTcAywFTgIub90uB05u0ycBV9TAzcD+SQ4ddZ2SpC3G+p1DkhXAUcAtwCFVtakt+gFwSJteCtw3tNqG1jZ7W2clWZtk7ebNm0dXtCQtQmMLhyTPAT4JvLOqHhteVlUF1Hy2V1WXVNWqqlo1NTW1gJVKksYSDkn2ZBAMH6+qT7Xm+2cOF7XnB1r7RmD50OrLWpskaUzGcbZSgEuBe6rqg0OLrgPOaNNnAJ8eaj+9nbV0DPDo0OEnSdIYjOMX0scBbwG+lmRda/td4P3ANUnOBL4H/Fpbdj1wIjANPAG8dQw1SpKGjDwcquomIFtZfMIc/Qs4Z6RFSZK2yV9IS5I6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqTPycEhyWZIHktw11PZ7STYmWdceJw4tOy/JdJJvJHndqOuTJPXGMXL4GPD6OdovrKqV7XE9QJIjgdXAL7V1PpRkyRhqlCQNGXk4VNVXgB/uYPeTgKuq6smquheYBo4eWXGSpDlN8juHtye5sx12OqC1LQXuG+qzobV1kpyVZG2StZs3bx51rZK0qEwqHD4MvBBYCWwCPjDfDVTVJVW1qqpWTU1NLXB5krS4TSQcqur+qnq6qn4CfIQth442AsuHui5rbZKkMZpIOCQ5dGj2jcDMmUzXAauT7J3kMOAI4NZx1ydJC+3Ui9dw6sVrJl3GDttj1P9Akk8ArwEOTrIBeC/wmiQrgQK+C5wNUFV3J7kGWA88BZxTVU+PukZJ0jONPByq6rQ5mi/dRv/zgfNHV5EkaXv8hbQkqWM4SJI6hoMkqTPy7xwkaXe00GcWrd/02Ei2e/XZxy7o9mY4cpAkdRw5SNIcFvoT+cyIYVSf9BeaIwdJUsdwkCR1DAdJUsdwkCR1DAdJUsezlSRpDHaXs5RmOHKQJHUMB0lSx3CQJHUMB0lSx3CQJHUMB0lSx3CQJHUMB0lSx3CQJHVSVZOu4aeWZDPwvQXY1MHAgwuwnVGyxoWzO9RpjQtnd6hz3DU+v6qm5lrwMxEOCyXJ2qpaNek6tsUaF87uUKc1Lpzdoc5dqUYPK0mSOoaDJKljODzTJZMuYAdY48LZHeq0xoWzO9S5y9Todw6SpI4jB0lSx3CQJHUWTTgk+SdJ1g09HkvyzrbsPyb5epK7k/zB0DrnJZlO8o0kr5tUjUlWJrm5ta1NcnTrnyQXtRrvTPKyUdc4VOu5bX/dleQTSfZJcliSW1o9VyfZq/Xdu81Pt+UrJljjx9v/z7uSXJZkz9Z3IvtyrhqHll2U5PGh+V1pPybJ+Um+meSeJL/Z+u5qr8kTktze3js3JTm89Z3UvnxHq+/uob8/Bya5Icm32vMBrX1i+xKAqlp0D2AJ8APg+cC/AL4A7N2W/Xx7PhK4A9gbOAz4NrBkQjV+HnhDaz8R+PLQ9GeBAMcAt4yptqXAvcC+bf4a4N+159Wt7c+At7Xp/wD8WZteDVw9wRpPbPsrwCeGahz7vtxajW16FXAl8PhQ/11pP74VuAJ4Vmufed/saq/JbwK/OLT/PjbBffkS4C7g2Qxu0fwF4HDgD4D3tD7vAS6Y5L6ceSyakcMsJwDfrqrvAW8D3l9VTwJU1QOtz0nAVVX1ZFXdC0wDR0+oxgKe29r3A74/VOMVNXAzsH+SQ8dU3x7Avkn2YPBi3wQcD1zbll8OnDxU5+Vt+lrghCSZQI3fr6rr2/4q4FZg2VCNk9iXXY1JlgB/CPz2rL67zH5k8L75/ar6CXTvm13lNfl9tv3eGfe+/EUGf+CfqKqngL8B3jSrltnvm0nty0UbDqsZfGoEeBHw6ja0/Jskr2jtS4H7htbZ0NomUeM7gT9Mch/wR8B5rX0iNVbVxlbH3zEIhUeB24BH2ot+di3/UGdb/ihw0LhrrKrPzyxvh5PeAnxudo1z1D/uGt8OXFdVm2atsivtxxcCp2ZwmPOzSY6YXWMzsddkq/M3gOuTbGDw//v9s+sc175kMGp4dZKDkjybwchgOXDI0P/rHwCHzK6xGevfoEUXDu04+K8Cf9ma9gAOZDBs+y3gmjF9GtuqOWp8G3BuVS0HzgUunVRtAO2Y6EkMDrf9AvCPgNdPsqbZ5qoxya8PdfkQ8JWq+ttJ1AdbrfF04M3AH0+qrmHb2I97Az+uwaUePgJcNrkqt1nnucCJVbUM+CjwwUnVWFX3ABcwOEz8OWAd8PSsPsVgtDNxiy4cgDcAt1fV/W1+A/CpNnS7FfgJg4tfbWSQ6jOWtbZJ1HgG8Kk2/ZdsObw1qRr/JXBvVW2uqr9vtR3HYNi7xxy1/EOdbfl+wEMTqPGfthreC0wB7xrqP4l9OVeN72NwHHo6yXeBZyeZnl3jLrAfN7DlNflXwC/PrrGZ9GvypVV1S+tzNe01wGT2JVV1aVW9vKr+GfAwg+9E7p85XNSeZw7RTfJv0KIMh9PYcrgG4K8ZfClNkhcBezG4KuJ1wOp2VsNhwBEMjlFPosbvA/+8TR8PfKtNXwec3s5qOIbBUHr2oYhR+DvgmCTPbqOsE4D1wI3AKa3PGcCnh+o8o02fAnypfUIad433JPkN4HXAaTPHy4dqHPe+nKvGD1bVP66qFVW1Aniiqg4fqnGX2I8MvW8YvDa/OVTjrvSa3K+9rwFe22qfqXPc+5IkP9+en8fg+4a/mFXL7PfNJPblwDi//Z70g8Hhj4eA/Yba9gL+nMHxwNuB44eW/ScGZyl9g3a20IRqfBWDY/p3ALcAL2/tAf601fg1YNUY9+X7gK+3/XYlg8MML2AQoNMMRjgzZ4Dt0+an2/IXTLDGp9r+Wtce/3mS+3KuGmctHz5baVfaj/sD/6vtqzUMPqHviq/JN7Y67gC+PLPPJrgv/5ZBaN0BnNDaDgK+yOBD3xeAAye9L6vKy2dIknqL8bCSJGk7DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdpJyU5OUklefE2+jzd7iVwV5L/mWT/+awvTYrhIO2804C17Xlr/m9VrayqlwA/BM6Ztf5N21lfmgjDQdoJSZ4DvIbBJaF39I/7Gtoll9v6rwLOZHB5dmmXYjhIO+ck4AtVdQfweJKXb6tzu4HPCQwupjaz/ueq6pvAQ9tbXxo3w0HaOacxuBUl7Xlro4d9k6xjy01cbhha/6o2fdU21pcmwgvvSfOU5EAGV+pdVlVPJnkBg1s+Pq9mvaGSPF5Vz2l3/vrfDK4E+ucM7oewmcGNXZa05+fPXl+aFEcO0vydAlxfW+47/h0Gt6Z89dZWqKongN8E3t3Wv7Kqnl+D+zYsB+7d1vrSuO2x/S6SZjkNeGm7U9uMg1r7V7a2UlV9NcmdwMUM7vY37JPbW18aJw8rSZI6HlaSJHUMB0lSx3CQJHUMB0lSx3CQJHUMB0lSx3CQJHX+PwEwUgBU9otMAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# generate data with default kwargs\n", "short_data_table, short_orbit_fraction = generate_synthetic_data()\n", "print(\"The orbit fraction is {}%\".format(np.round(short_orbit_fraction), 1))\n", "\n", "# initialize orbitize `System` object\n", "short_system = system.System(1, short_data_table, mtot, plx)\n", "\n", "num2accept = 500 # run sampler until this many orbits are accepted" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Converting ra/dec data points in data_table to sep/pa. Original data are stored in input_table.\n", "497/500 orbits found\r" ] }, { "ename": "KeyboardInterrupt", "evalue": "", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", "Input \u001b[0;32mIn [4]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 4\u001b[0m short_OFTI_sampler \u001b[38;5;241m=\u001b[39m sampler\u001b[38;5;241m.\u001b[39mOFTI(short_system)\n\u001b[1;32m 6\u001b[0m \u001b[38;5;66;03m# perform OFTI fit\u001b[39;00m\n\u001b[0;32m----> 7\u001b[0m short_OFTI_orbits \u001b[38;5;241m=\u001b[39m \u001b[43mshort_OFTI_sampler\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun_sampler\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnum2accept\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 9\u001b[0m \u001b[38;5;28mprint\u001b[39m(\n\u001b[1;32m 10\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mOFTI took \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m seconds to accept \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m orbits.\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mformat(\n\u001b[1;32m 11\u001b[0m time\u001b[38;5;241m.\u001b[39mtime() \u001b[38;5;241m-\u001b[39m start_time, num2accept\n\u001b[1;32m 12\u001b[0m )\n\u001b[1;32m 13\u001b[0m )\n", "File \u001b[0;32m~/Documents/GitHub/orbitize/orbitize/sampler.py:593\u001b[0m, in \u001b[0;36mOFTI.run_sampler\u001b[0;34m(self, total_orbits, num_samples, num_cores, OFTI_warning)\u001b[0m\n\u001b[1;32m 588\u001b[0m OFTI_warning \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[1;32m 589\u001b[0m \u001b[38;5;28mprint\u001b[39m(\n\u001b[1;32m 590\u001b[0m \u001b[38;5;28mstr\u001b[39m(orbits_saved\u001b[38;5;241m.\u001b[39mvalue) \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m/\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mstr\u001b[39m(total_orbits) \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m orbits found\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m 591\u001b[0m end\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;130;01m\\r\u001b[39;00m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m 592\u001b[0m )\n\u001b[0;32m--> 593\u001b[0m \u001b[43mtime\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msleep\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m0.1\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 595\u001b[0m \u001b[38;5;28mprint\u001b[39m(\n\u001b[1;32m 596\u001b[0m \u001b[38;5;28mstr\u001b[39m(total_orbits) \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m/\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mstr\u001b[39m(total_orbits) \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m orbits found\u001b[39m\u001b[38;5;124m\"\u001b[39m, end\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;130;01m\\r\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 597\u001b[0m )\n\u001b[1;32m 599\u001b[0m \u001b[38;5;66;03m# join the processes\u001b[39;00m\n", "\u001b[0;31mKeyboardInterrupt\u001b[0m: " ] } ], "source": [ "start_time = time.time()\n", "\n", "# set up OFTI `Sampler` object\n", "short_OFTI_sampler = sampler.OFTI(short_system)\n", "\n", "# perform OFTI fit\n", "short_OFTI_orbits = short_OFTI_sampler.run_sampler(num2accept)\n", "\n", "print(\n", " \"OFTI took {} seconds to accept {} orbits.\".format(\n", " time.time() - start_time, num2accept\n", " )\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "start_time = time.time()\n", "\n", "# set up MCMC `Sampler` object\n", "num_walkers = 20\n", "short_MCMC_sampler = sampler.MCMC(short_system, num_temps=5, num_walkers=num_walkers)\n", "\n", "# perform MCMC fit\n", "num2accept_mcmc = 10 * num2accept\n", "_ = short_MCMC_sampler.run_sampler(num2accept_mcmc, burn_steps=100)\n", "short_MCMC_orbits = short_MCMC_sampler.results.post\n", "\n", "print(\n", " \"MCMC took {} steps in {} seconds.\".format(\n", " num2accept_mcmc, time.time() - start_time\n", " )\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.hist(\n", " short_OFTI_orbits[:, short_system.param_idx[\"ecc1\"]],\n", " bins=40,\n", " density=True,\n", " alpha=0.5,\n", " label=\"OFTI\",\n", ")\n", "plt.hist(\n", " short_MCMC_orbits[:, short_system.param_idx[\"ecc1\"]],\n", " bins=40,\n", " density=True,\n", " alpha=0.5,\n", " label=\"MCMC\",\n", ")\n", "\n", "plt.xlabel(\"Eccentricity\")\n", "plt.ylabel(\"Prob.\")\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These distributions are different because the MCMC chains have not converged, resulting in a \"lumpy\" MCMC distribution. I set up the calculation so that MCMC would return 10x as many orbits as OFTI, but even so, the OFTI distribution is a much better representation of the underlying PDF. \n", "\n", "If we run the MCMC algorithm for a greater number of steps (and/or increase the number of walkers and/or temperatures), the MCMC and OFTI distributions will become indistinguishable. **OFTI is NOT more correct than MCMC, but for this dataset, OFTI converges on the correct posterior faster than MCMC**. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Longer Orbit Fraction\n", "\n", "Let's now repeat this exercise with a longer orbit fraction. For this dataset, OFTI will have to run for several seconds just to accept one orbit, so we won't compare the resulting posteriors." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# generate data\n", "long_data_table, long_orbit_fraction = generate_synthetic_data(sma=10, num_obs=5)\n", "print(\"The orbit fraction is {}%\".format(np.round(long_orbit_fraction), 1))\n", "\n", "# initialize orbitize `System` object\n", "long_system = system.System(1, long_data_table, mtot, plx)\n", "num2accept = 500 # run sampler until this many orbits are accepted" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "start_time = time.time()\n", "\n", "# set up OFTI `Sampler` object\n", "long_OFTI_sampler = sampler.OFTI(long_system)\n", "\n", "# perform OFTI fit\n", "long_OFTI_orbits = long_OFTI_sampler.run_sampler(1)\n", "\n", "print(\"OFTI took {} seconds to accept 1 orbit.\".format(time.time() - start_time))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "start_time = time.time()\n", "\n", "# set up MCMC `Sampler` object\n", "num_walkers = 20\n", "long_MCMC_sampler = sampler.MCMC(long_system, num_temps=10, num_walkers=num_walkers)\n", "\n", "# perform MCMC fit\n", "_ = long_MCMC_sampler.run_sampler(num2accept, burn_steps=100)\n", "long_MCMC_orbits = long_MCMC_sampler.results.post\n", "\n", "print(\"MCMC took {} steps in {} seconds.\".format(num2accept, time.time() - start_time))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.hist(long_MCMC_orbits[:, short_system.param_idx[\"ecc1\"]], bins=15, density=True)\n", "plt.xlabel(\"Eccentricity\")\n", "plt.ylabel(\"Prob.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It will take more steps for this MCMC to fully converge (see the [MCMC tutorial](https://orbitize.readthedocs.io/en/latest/tutorials/MCMC_tutorial.html) for more detailed guidelines), but you can imagine that MCMC will converge much faster than OFTI for this dataset." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Closing Thoughts\n", "\n", "If you play around with the `num_obs`, `sma`, and `unc` keywords in the `generate_synthetic_data` function and repeat this exercise, you will notice that the OFTI acceptance rate and MCMC convergence rate depend on many variables, not just orbit fraction. **In truth, the Gaussianity of the posterior space determines how quickly an MCMC run will converge, and its similarity to the prior space determines how quickly an OFTI run will converge. In other words, the more constrained your posteriors are (relative to your priors), the quicker MCMC will converge, and the slower OFTI will run.**\n", "\n", "Orbit fraction is usually a great tracer of this \"amount of constraint,\" but it's good to understand why!\n", "\n", "**Summary**:\n", "- OFTI and MCMC produce the same posteriors, but often take differing amounts of time to converge on the correct solution.\n", "- OFTI is superior when your posteriors are similar to your priors, and MCMC is superior when your posteriors are highly constrained Gaussians." ] } ], "metadata": { "interpreter": { "hash": "e899b22145868d3cd465733d82c36c2ae3ac0d3591d6a0807ec2e5e577a9cf5c" }, "kernelspec": { "display_name": "Python 3.7.5 64-bit ('python3.7': conda)", "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.10.4" } }, "nbformat": 4, "nbformat_minor": 2 }