Skip to content

Commit

Permalink
add some more explanation of syntax to tutorial
Browse files Browse the repository at this point in the history
  • Loading branch information
@sblunt
sblunt committed Mar 1, 2024
1 parent 8d6f750 commit beb9347
Showing 1 changed file with 42 additions and 46 deletions.
88 changes: 42 additions & 46 deletions docs/tutorials/ONeil-ObsPriors.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,18 +9,13 @@
"by Sarah Blunt and Clarissa R. Do Ó (2024)\n",
"\n",
"\n",
"### Background on Observable-base Priors \n",
"### Background on Observable-based 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",
"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"
"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"
]
},
{
Expand All @@ -40,13 +35,14 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"### Set up MCMC run\n",
"### Set up Orbit\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",
"Let's set up our System class, which holds all the info about the planetary system we're modeling. 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",
"We begin but reading in the astrometry data, and providing 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."
"Importantly, to use the observation-based priors in orbitize!, **you must also use the \"ObsPriors\"\n",
"fitting basis** (as shown below)."
]
},
{
Expand All @@ -72,8 +68,6 @@
"\n",
"data_table = read_input.read_file(filename)\n",
"\n",
"# TODO (sarah): note that the fitting basis choice here is required (can't fit in any other orbitize! basis)\n",
"# if we want to use obspriors\n",
"mySystem = system.System(\n",
" 1, data_table, system_mass, plx, mass_err=0, plx_err=0, fitting_basis=\"ObsPriors\"\n",
")\n",
Expand Down Expand Up @@ -143,35 +137,26 @@
"source": [
"## Modify Priors\n",
"\n",
"Define the priors on `sma`, `ecc`, `tau`, `plx`, and `mtot` to be the O'Neil observation-based prior."
"We'll next modify the priors on `per1`, `ecc1`, and `tp1` to be the O'Neil observation-based prior. For more info on modifying priors in `orbitize!`, check out the [modifying priors tutorial](https://orbitize.readthedocs.io/en/latest/tutorials/Modifying_Priors.html). \n",
"\n",
"The observation-based prior requires information about the error bars and observation times of your data points, so we need to pass those in."
]
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/bluez3303/miniconda3/envs/python3.11/lib/python3.11/site-packages/erfa/core.py:154: ErfaWarning: ERFA function \"dtf2d\" yielded 1 of \"dubious year (Note 6)\"\n",
" warnings.warn('ERFA function \"{}\" yielded {}'.format(func_name, wmsg),\n",
"/Users/bluez3303/miniconda3/envs/python3.11/lib/python3.11/site-packages/erfa/core.py:154: ErfaWarning: ERFA function \"utctai\" yielded 1 of \"dubious year (Note 3)\"\n",
" warnings.warn('ERFA function \"{}\" yielded {}'.format(func_name, wmsg),\n",
"/Users/bluez3303/miniconda3/envs/python3.11/lib/python3.11/site-packages/erfa/core.py:154: ErfaWarning: ERFA function \"taiutc\" yielded 1 of \"dubious year (Note 4)\"\n",
" warnings.warn('ERFA function \"{}\" yielded {}'.format(func_name, wmsg),\n"
"ename": "NameError",
"evalue": "name 'mySystem' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[2], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m ra_err \u001b[38;5;241m=\u001b[39m \u001b[43mmySystem\u001b[49m\u001b[38;5;241m.\u001b[39mdata_table[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mquant1_err\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n\u001b[1;32m 2\u001b[0m dec_err \u001b[38;5;241m=\u001b[39m mySystem\u001b[38;5;241m.\u001b[39mdata_table[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mquant1_err\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n\u001b[1;32m 3\u001b[0m epochs \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray(mySystem\u001b[38;5;241m.\u001b[39mdata_table[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mepoch\u001b[39m\u001b[38;5;124m\"\u001b[39m])\n",
"\u001b[0;31mNameError\u001b[0m: name 'mySystem' is not defined"
]
},
{
"data": {
"text/plain": [
"[ObsPrior, ObsPrior, Sine, Uniform, Uniform, ObsPrior, 24.5272375, 1.35]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
Expand All @@ -188,14 +173,17 @@
" dec_err,\n",
" epochs,\n",
" system_mass,\n",
" period_lims=(0.0, 100),\n",
" period_lims=(\n",
" 0.0,\n",
" 100,\n",
" ), # these two \"limits\" on period and periastron passage time are optional, but aid convergence.\n",
" tp_lims=(\n",
" lower_tp_lim,\n",
" upper_tp_lim,\n",
" ),\n",
")\n",
"\n",
"# set the priors on `per`, `ecc`, and `tp` to point to this object\n",
"# set the priors on `per1`, `ecc1`, and `tp1` to point to this object\n",
"for i in [\n",
" mySystem.param_idx[\"per1\"],\n",
" mySystem.param_idx[\"ecc1\"],\n",
Expand Down Expand Up @@ -351,21 +339,29 @@
}
],
"source": [
"# create an MCMC object using our newly modified System object\n",
"# these are the MCMC parameters I'd use if I wanted to publish this orbit fit.\n",
"# num_temps = 5\n",
"# num_walkers = 50\n",
"# num_threads = 2\n",
"# n_steps_per_walker = 100000\n",
"# n_burn_steps = 4000\n",
"\n",
"# since this is a tutorial, let's only run it for a few steps\n",
"num_temps = 5\n",
"num_walkers = 50\n",
"num_threads = mp.cpu_count()\n",
"num_threads = 2\n",
"n_steps_per_walker = 100\n",
"n_burn_steps = 10\n",
"\n",
"# create an MCMC object using our newly modified System object\n",
"myMCMC = sampler.MCMC(\n",
" mySystem, num_temps=num_temps, num_walkers=num_walkers, num_threads=num_threads\n",
")\n",
"\n",
"# thes\n",
"n_steps_per_walker = 100000\n",
"n_burn_steps = 4000\n",
"\n",
"myMCMC.check_prior_support()\n",
"\n",
"myMCMC.run_sampler(n_steps_per_walker * num_walkers, burn_steps=10, examine_chains=True)"
"myMCMC.run_sampler(\n",
" n_steps_per_walker * num_walkers, burn_steps=n_burn_steps, examine_chains=True\n",
")"
]
},
{
Expand Down Expand Up @@ -407,7 +403,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.8"
"version": "3.12.2"
}
},
"nbformat": 4,
Expand Down

0 comments on commit beb9347

Please sign in to comment.