Skip to content

Commit

Permalink
update obs priors tutorial
Browse files Browse the repository at this point in the history
  • Loading branch information
Clarissa Rizzo authored and Clarissa Rizzo committed Mar 1, 2024
1 parent 5e23a70 commit 5f19714
Showing 1 changed file with 100 additions and 29 deletions.
129 changes: 100 additions & 29 deletions docs/tutorials/ONeil-ObsPriors.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <a href = \"https://ui.adsabs.harvard.edu/abs/2019AJ....158....4O/abstract\"> O'Neil et al 2019</a>, 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": [
Expand All @@ -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",
Expand Down Expand Up @@ -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": [
"<div><i>Table length=9</i>\n",
"<table id=\"table140558710839472\" class=\"table-striped table-bordered table-condensed\">\n",
"<thead><tr><th>epoch</th><th>object</th><th>quant1</th><th>quant1_err</th><th>quant2</th><th>quant2_err</th><th>quant12_corr</th><th>quant_type</th><th>instrument</th></tr></thead>\n",
"<thead><tr><th>float64</th><th>int64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>bytes5</th><th>bytes5</th></tr></thead>\n",
"<tr><td>57298.0</td><td>1</td><td>-253.72</td><td>2.98</td><td>92.35</td><td>2.85</td><td>nan</td><td>radec</td><td>defrd</td></tr>\n",
"<tr><td>57606.0</td><td>1</td><td>-236.63</td><td>9.77</td><td>127.94</td><td>9.18</td><td>nan</td><td>radec</td><td>defrd</td></tr>\n",
"<tr><td>57645.0</td><td>1</td><td>-234.52</td><td>1.79</td><td>123.39</td><td>1.03</td><td>nan</td><td>radec</td><td>defrd</td></tr>\n",
"<tr><td>57946.0</td><td>1</td><td>-210.76</td><td>1.94</td><td>152.09</td><td>1.88</td><td>nan</td><td>radec</td><td>defrd</td></tr>\n",
"<tr><td>58276.0</td><td>1</td><td>-167.49</td><td>1.61</td><td>180.87</td><td>16.97</td><td>nan</td><td>radec</td><td>defrd</td></tr>\n",
"<tr><td>58287.0</td><td>1</td><td>-177.67</td><td>1.67</td><td>174.6</td><td>1.67</td><td>nan</td><td>radec</td><td>defrd</td></tr>\n",
"<tr><td>58365.0</td><td>1</td><td>-165.7</td><td>3.28</td><td>185.33</td><td>3.66</td><td>nan</td><td>radec</td><td>defrd</td></tr>\n",
"<tr><td>58368.0</td><td>1</td><td>-170.38</td><td>2.52</td><td>185.94</td><td>2.74</td><td>nan</td><td>radec</td><td>defrd</td></tr>\n",
"<tr><td>58414.0</td><td>1</td><td>-161.64</td><td>13.6</td><td>176.21</td><td>14.31</td><td>nan</td><td>radec</td><td>defrd</td></tr>\n",
"</table></div>"
],
"text/plain": [
"[ObsPrior, ObsPrior, Sine, Uniform, Uniform, ObsPrior, ObsPrior, ObsPrior]"
"<Table length=9>\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<ipython-input-8-1c048fde3cbd>\u001b[0m in \u001b[0;36m<module>\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",
Expand All @@ -130,7 +199,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"### Run MCMC!"
"### Run MCMC!\n",
"\n",
"Let's run this orbit fit!"
]
},
{
Expand Down Expand Up @@ -173,7 +244,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "python3.11",
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
Expand All @@ -187,7 +258,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.4"
"version": "3.8.8"
}
},
"nbformat": 4,
Expand Down

0 comments on commit 5f19714

Please sign in to comment.