From 5f19714ba6ab0b430623c44fee74b837be1712e4 Mon Sep 17 00:00:00 2001 From: Clarissa Rizzo Date: Fri, 1 Mar 2024 12:33:09 -0800 Subject: [PATCH] update obs priors tutorial --- docs/tutorials/ONeil-ObsPriors.ipynb | 129 +++++++++++++++++++++------ 1 file changed, 100 insertions(+), 29 deletions(-) diff --git a/docs/tutorials/ONeil-ObsPriors.ipynb b/docs/tutorials/ONeil-ObsPriors.ipynb index bd7f208b..01c29644 100644 --- a/docs/tutorials/ONeil-ObsPriors.ipynb +++ b/docs/tutorials/ONeil-ObsPriors.ipynb @@ -6,12 +6,26 @@ "source": [ "# Using the O'Neil (2019) Observation-based Priors\n", "\n", - "by Sarah Blunt (2024)" + "by Sarah Blunt and Clarissa R. Do Ó (2024)\n", + "\n", + "\n", + "### Background on Observable-base Priors \n", + "\n", + "Observable-based priors are a different set of priors from the standard \"parameter\" priors that most orbit fitting codes (including orbitize!) use. Here, rather than placing uniform priors on orbital parameters such as eccentricity and periastron passage, uniform priors are placed on the \"observables\" of the orbit (which are linearly related to the measured positions of the companion). The idea behind these priors is that there is an equal probability (i.e., uniform probability) of obtaining observations in the regions of parameter space that are possible to observe. For that reason, observable priors present uniformity in the orbital observables, which can be transformed back to orbital parameters.\n", + "\n", + "The paper where these priors are introduced is O'Neil et al 2019, so if you use this functionality please cite that paper.\n", + "\n", + "Note: The observable-based priors are set up in a different orbital basis from orbitize!'s basis. The main difference is that they use period P rather than semi-major axis a and periastron passage epoch (in years) rather than the unitless parameter $\\tau$.\n", + "\n", + "\n", + "### Initial Set up\n", + "\n", + "First, we import orbitize's driver and the \"priors\" functionality:\n" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -28,21 +42,27 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Set up MCMC run" + "### Set up MCMC run\n", + "\n", + "Now, we will set up our MCMC run. Here, we will use the object HD 206893 B as an example orbit for this tutorial. \n", + "\n", + "We begin but providing the astrometry data in CSV format and the system's mass and parallax. For this tutorial, and in general, one should use RA and Dec astrometry measurements rather than separation/P.A..\n", + "\n", + "We also set up the number of temps, walkers and threads for our MCMC run and set up the driver object." ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ - "filename = \"{}/GJ504.csv\".format(orbitize.DATADIR)\n", + "filename = \"{}/hd206893b.csv\".format(orbitize.DATADIR)\n", "\n", "# system parameters\n", "num_secondary_bodies = 1\n", - "system_mass = 1.75 # [Msol]\n", - "plx = 51.44 # [mas]\n", + "system_mass = 1.35 # [Msol]\n", + "plx = 24.5272375 # [mas]\n", "\n", "# MCMC parameters\n", "num_temps = 10\n", @@ -70,44 +90,93 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Modify Priors\n", - "\n", - "Define the priors on `sma`, `ecc`, `tau`, `plx`, and `mtot` to be the O'Neil observation-based prior." + "Let's check if the data is being read correctly!" ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 7, "metadata": {}, "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "{'sma1': 0, 'ecc1': 1, 'inc1': 2, 'aop1': 3, 'pan1': 4, 'tau1': 5, 'plx': 6, 'mtot': 7}\n" - ] - }, { "data": { + "text/html": [ + "
Table length=9\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
epochobjectquant1quant1_errquant2quant2_errquant12_corrquant_typeinstrument
float64int64float64float64float64float64float64bytes5bytes5
57298.01-253.722.9892.352.85nanradecdefrd
57606.01-236.639.77127.949.18nanradecdefrd
57645.01-234.521.79123.391.03nanradecdefrd
57946.01-210.761.94152.091.88nanradecdefrd
58276.01-167.491.61180.8716.97nanradecdefrd
58287.01-177.671.67174.61.67nanradecdefrd
58365.01-165.73.28185.333.66nanradecdefrd
58368.01-170.382.52185.942.74nanradecdefrd
58414.01-161.6413.6176.2114.31nanradecdefrd
" + ], "text/plain": [ - "[ObsPrior, ObsPrior, Sine, Uniform, Uniform, ObsPrior, ObsPrior, ObsPrior]" + "\n", + " epoch object quant1 quant1_err ... quant12_corr quant_type instrument\n", + "float64 int64 float64 float64 ... float64 bytes5 bytes5 \n", + "------- ------ ------- ---------- ... ------------ ---------- ----------\n", + "57298.0 1 -253.72 2.98 ... nan radec defrd\n", + "57606.0 1 -236.63 9.77 ... nan radec defrd\n", + "57645.0 1 -234.52 1.79 ... nan radec defrd\n", + "57946.0 1 -210.76 1.94 ... nan radec defrd\n", + "58276.0 1 -167.49 1.61 ... nan radec defrd\n", + "58287.0 1 -177.67 1.67 ... nan radec defrd\n", + "58365.0 1 -165.7 3.28 ... nan radec defrd\n", + "58368.0 1 -170.38 2.52 ... nan radec defrd\n", + "58414.0 1 -161.64 13.6 ... nan radec defrd" ] }, - "execution_count": 13, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], + "source": [ + "my_driver.system.data_table" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Modify Priors\n", + "\n", + "Define the priors on `sma`, `ecc`, `tau`, `plx`, and `mtot` to be the O'Neil observation-based prior." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "ename": "AttributeError", + "evalue": "module 'orbitize.priors' has no attribute 'ObsPrior'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 13\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 14\u001b[0m \u001b[0;31m# define the `ObsPrior` object\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 15\u001b[0;31m \u001b[0mmy_obsprior\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpriors\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mObsPrior\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mra_err\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdec_err\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepochs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 16\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 17\u001b[0m \u001b[0;31m# TODO: change the basis\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mAttributeError\u001b[0m: module 'orbitize.priors' has no attribute 'ObsPrior'" + ] + } + ], "source": [ "# convert input sep/PA measurements to RA/decl\n", - "sep = np.array(my_driver.system.data_table[\"quant1\"])\n", - "sep_err = np.array(my_driver.system.data_table[\"quant1_err\"])\n", - "pa = np.radians(np.array(my_driver.system.data_table[\"quant2\"]))\n", + "ra = np.array(my_driver.system.data_table[\"quant1\"])\n", + "ra_err = np.array(my_driver.system.data_table[\"quant1_err\"])\n", + "dec = np.radians(np.array(my_driver.system.data_table[\"quant2\"]))\n", "pa_err = np.radians(np.array(my_driver.system.data_table[\"quant2_err\"]))\n", + "[]\n", "\n", - "\n", - "ra_err = np.sqrt((np.cos(pa) * sep_err) ** 2 + (sep * np.sin(pa) * pa_err) ** 2)\n", - "dec_err = np.sqrt((np.sin(pa) * sep_err) ** 2 + (sep * np.cos(pa) * pa_err) ** 2)\n", + "#ra_err = np.sqrt((np.cos(pa) * sep_err) ** 2 + (sep * np.sin(pa) * pa_err) ** 2)\n", + "#dec_err = np.sqrt((np.sin(pa) * sep_err) ** 2 + (sep * np.cos(pa) * pa_err) ** 2)\n", "\n", "\n", "epochs = np.array(my_driver.system.data_table[\"epoch\"])\n", @@ -130,7 +199,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Run MCMC!" + "### Run MCMC!\n", + "\n", + "Let's run this orbit fit!" ] }, { @@ -173,7 +244,7 @@ ], "metadata": { "kernelspec": { - "display_name": "python3.11", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -187,7 +258,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.8.8" } }, "nbformat": 4,