Tutorial: Fit binary model to observed data using nested sampling via UltraNest#
In this tutorial we go through the procedure of fitting observations to a binary star model using the nested sampling code UltraNest.
Imports and Filter Setup#
from phitter import observables, filters
from phitter.params import star_params, binary_params, isoc_interp_params
from phitter.calc import model_obs_calc, phot_adj_calc, rv_adj_calc
from phitter.fit import likelihood, prior
import numpy as np
from phoebe import u
from phoebe import c as const
import matplotlib as mpl
import matplotlib.pyplot as plt
import pickle
from astropy.table import Table
import ultranest
import ultranest.stepsampler
%matplotlib inline
# The following warning regarding extinction originates from SPISEA and can be ignored.
# The functionality being warned about is not used by SPISEA.
/Users/abhimatlocal/software/miniforge3/envs/phoebe_py38/lib/python3.8/site-packages/pysynphot/locations.py:345: UserWarning: Extinction files not found in /Users/abhimatlocal/models/cdbs/extinction
warnings.warn('Extinction files not found in %s' % (extdir, ))
/Users/abhimatlocal/software/miniforge3/envs/phoebe_py38/lib/python3.8/site-packages/pysynphot/refs.py:124: UserWarning: No thermal tables found, no thermal calculations can be performed. No files found for /Users/abhimatlocal/models/cdbs/mtab/*_tmt.fits
warnings.warn('No thermal tables found, '
filter_153m = filters.hst_f153m_filt()
filter_127m = filters.hst_f127m_filt()
Prepare observations data#
Let’s read in the mock data generated in a separate example here. Let’s assume this data is for a red giant binary star belonging to an 8 Gyr old star population at the Galactic center.
with open('./mock_obs_table.pkl', 'rb') as in_file:
obs_table = pickle.load(in_file)
print(obs_table)
obs_times obs obs_uncs obs_types obs_filts
------------------ ------------------ -------- --------- ------------------------------------------------------
53800.714762925796 16.769110541124228 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53802.78887450278 16.569960382011555 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53804.08607971392 16.523859319835843 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53804.64136070461 16.54425287095397 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53805.203583377115 16.551835703948942 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53805.96328097839 16.559083907180494 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53806.315847597616 16.526208341236522 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53806.333495717256 16.549137149396568 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53807.488405082564 16.54357246902387 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53807.67541617193 16.512582058125357 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53809.991688679045 16.57630387875109 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53810.3952177884 16.571876780401222 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53810.56053852182 16.5629106720859 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53812.06192493336 16.79990290964248 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53814.764347052136 16.577347274346987 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53814.870873942076 16.584590478742268 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53815.35036219548 16.55250309097789 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53816.915279343724 16.55815171414512 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53817.22147451534 16.56401507633422 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53818.297738065274 16.527494721029015 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53822.21481869134 16.55859419544436 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53823.32369330329 16.60712912435184 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53823.33777948042 16.631108542728825 0.015 phot <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
... ... ... ... ...
53811.39251280399 19.466756200612434 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53812.36929572477 19.60648938621428 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53813.378127665645 19.529182837164768 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53814.38351959213 19.31805797850974 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53814.39102396979 19.342402712730358 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53815.47248701535 19.324506542981187 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53815.71319961314 19.301216716218455 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53818.81438132841 19.238703441352225 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53821.48902609574 19.29472566814667 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53822.61199465376 19.28945924040662 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53822.83753761046 19.346641516593092 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53824.797307830486 19.601928270465987 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53824.82761134125 19.606919387014845 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53824.827858456956 19.58491117884066 0.015 phot <phitter.filters.hst_f127m_filt object at 0x169a462e0>
53802.139836776325 130.228146222555 10.0 rv_pri <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53803.32979530006 114.02051863903269 10.0 rv_pri <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53804.46291200766 114.54568517998948 10.0 rv_pri <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53813.273397048906 163.65481126090992 10.0 rv_pri <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53819.03215627737 175.28898537909566 10.0 rv_pri <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53802.139836776325 164.33130628324724 10.0 rv_sec <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53803.32979530006 176.97897304376025 10.0 rv_sec <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53804.46291200766 190.80797816955624 10.0 rv_sec <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53813.273397048906 129.70047968985767 10.0 rv_sec <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
53819.03215627737 97.50308377106813 10.0 rv_sec <phitter.filters.hst_f153m_filt object at 0x1041bd6a0>
Length = 60 rows
We create two new phitter.observables objects, one for modeling containing only the times, filters, and types of observations, and another for computing the likelihood, which will also contain the observations and associated uncertainties as well.
# Model observables object, which only contains times and types of observations
model_observables = observables.observables(
obs_times=obs_table['obs_times'].data,
obs_filts=obs_table['obs_filts'].data, obs_types=obs_table['obs_types'].data,
)
# An observables object for the observations, used when computing likelihoods
observations = observables.observables(
obs_times=obs_table['obs_times'].data, obs=obs_table['obs'].data, obs_uncs=obs_table['obs_uncs'].data,
obs_filts=obs_table['obs_filts'].data, obs_types=obs_table['obs_types'].data,
)
Make stellar parameters and binary parameters objects for fitting#
Now we can make a stellar parameters object that we will use to derive the stellar parameters from an isochrone while fitting. Let’s use our previous assumptions about what type of star this is: red giant in an 8 Gyr old star population.
We also can create a binary parameters object that will store the binary parameters.
isoc_stellar_params_obj = isoc_interp_params.isoc_mist_stellar_params(
age=8e9,
met=0.0,
use_atm_func='merged',
phase='RGB',
ext_Ks=2.2,
dist=8e3*u.pc,
filts_list=[filter_153m, filter_127m],
ext_law='NL18',
)
# Make binary params objects
bin_params = binary_params.binary_params()
Set up model and likelihood objects#
In order to carry out our fitting, we will generate two more objects: a model object that will be used to compute observables and a likelihood object that will compute the likelihood from our observations.
# Set up a binary model object
binary_model_obj = model_obs_calc.binary_star_model_obs(
model_observables,
use_blackbody_atm=False,
print_diagnostics=False,
)
# Set up likelihood object for fitting parameters
log_like_obj = likelihood.log_likelihood_chisq(
observations
)
Set up function for UltraNest to return likelihood from fit parameters#
Now we need to write a function that will be used by Ultranest for its sampling. This function will provide Ultranest with a likelihood for a given set of fit parameters.
For our fitting, we’ll fit our data to eight parameters:
radius of star 1
radius of star 2
binary orbital period
binary inclination
system t0
binary CoM velocity
extinction in F153M passband
extinction law alpha
def un_evaluate(model_params):
(
star1_radius,
star2_radius,
bin_period,
bin_inc,
bin_t0,
bin_rv_com,
ext_153m,
ext_alpha,
) = model_params
# Obtain stellar params by interpolating along the isochrone
star1_params = isoc_stellar_params_obj.interp_star_params_rad(
star1_radius,
)
star2_params = isoc_stellar_params_obj.interp_star_params_rad(
star2_radius,
)
# Set binary params
bin_params.period = bin_period * u.d
bin_params.inc = bin_inc * u.deg
bin_params.t0 = bin_t0
# Run binary model
modeled_observables = binary_model_obj.compute_obs(
star1_params, star2_params, bin_params,
num_triangles=300,
)
# Check for situation where binary model fails
# (i.e., unphysical conditions not able to be modeled)
if np.isnan(modeled_observables.obs_times[0]):
return -1e300
# Apply distance modulus
# (We're assuming we know the distance, but this can be a fit parameter as well)
modeled_observables = phot_adj_calc.apply_distance_modulus(
modeled_observables,
8e3*u.pc,
)
# Apply extinction
modeled_observables = phot_adj_calc.apply_extinction(
modeled_observables,
2.2, filter_153m,
ext_153m,
isoc_red_law='NL18',
ext_alpha=ext_alpha,
)
# Add RV center of mass velocity
modeled_observables = rv_adj_calc.apply_com_velocity(
modeled_observables,
bin_rv_com * u.km / u.s,
)
# Compute and return log likelihood
log_like = log_like_obj.evaluate(modeled_observables)
return log_like
Let’s test out our un_evaluate() function on some mock binary parameters
test_params = (
18., 10.,
23., 85.,
53_801.,
140.,
4.4, 2.22,
)
log_like_test = un_evaluate(test_params)
print(f'log likelihood from test parameters = {log_like_test:.3f}')
# Since we simulated these mock data ourselves, we know what the "truth" is here. So we can check that as well
truth_params = (
15., 12.,
25., 75.,
53_800.,
150.,
4.5, 2.23,
)
log_like_truth = un_evaluate(truth_params)
print(f'log likelihood from "truth" parameters = {log_like_truth:.3f}')
log likelihood from test parameters = -43318.653
log likelihood from "truth" parameters = -64.610
Set up priors#
We need to set up a prior for each of our fit parameters. phitter.fit.prior provides some useful functions to create priors and then put them into a collection that can be passed to Ultranest or other sampling codes.
Using the function phitter.fit.prior.prior_collection.prior_transform_ultranest, we can transform the unit cube passed in by UltraNest into physical parameters to use during fitting.
star1_rad_prior = prior.uniform_prior(10.0, 25.0)
star2_rad_prior = prior.uniform_prior(8.0, 15.0)
bin_period_prior = prior.uniform_prior(24.0, 26.0)
bin_inc_prior = prior.uniform_prior(0.0, 180.0)
bin_t0_prior = prior.uniform_prior(53_795.0, 53_805.0)
bin_rv_com_prior = prior.uniform_prior(100., 200.)
ext_f153m_prior = prior.uniform_prior(4, 6)
# We have a constraint on the extinction law (from Nogueras-Lara et al., 2019 here),
# so we can set Gaussian prior on the extinction law alpha
ext_alpha_prior = prior.gaussian_prior(2.23, 0.03)
param_priors = prior.prior_collection([
star1_rad_prior,
star2_rad_prior,
bin_period_prior,
bin_inc_prior,
bin_t0_prior,
bin_rv_com_prior,
ext_f153m_prior,
ext_alpha_prior,
])
param_names = [
'star1_rad',
'star2_rad',
'bin_period',
'bin_inc',
'bin_t0',
'bin_rv_com',
'ext_f153m',
'ext_alpha',
]
# Ultranest supports wrapped / circular parameters, which inclination is in this example.
# So we need to make a list to pass to Ultranest for which parameters are wrapped so that
# it can appropriately handle it during sampling.
wrapped_params = [
False, # Star 1 Rad
False, # Star 2 Rad
False, # Binary period
True, # Binary inclination
False, # Binary t0
False, # Binary CoM RV
False, # F153M-band extinction
False, # Extinction law alpha
]
Set up UltraNest Sampler and Run Sampler#
sampler = ultranest.ReactiveNestedSampler(
param_names,
loglike = un_evaluate, # Function that we wrote for ultranest to calculate likelihood
transform = param_priors.prior_transform_ultranest, # Phitter's prior collection function to transform prior
log_dir='./un_out/',
resume='resume-similar',
warmstart_max_tau=0.25,
)
# Below is an example of how to set up a slice sampler in UltraNest.
# We won't use it here since it is less efficient for this example, but a
# "slice sampler" may be more efficient when fitting a large number
# of parameters (i.e., >~ 10).
# sampler.stepsampler = ultranest.stepsampler.SliceSampler(
# nsteps=(2**3)*len(param_names),
# generate_direction=ultranest.stepsampler.generate_mixture_random_direction,
# )
The syntax below shows how to run the sampler. However for speed, I wrote a separate script that can be parallelized with MPI. In order to run the example script, the following syntax can be used in the shell:
mpiexec -np [num_cores] ./fit_with_ultranest.py
Importantly, the script demonstrates how to parallelize Ultranest’s sampling and turn off parallelization in PHOEBE which can interfere with parallelization in Ultranest:
import os
# Turn off parallelisation in phoebe
# Needs to be done *before* phoebe is imported
os.environ["OMP_NUM_THREADS"] = "1"
os.environ["PHOEBE_ENABLE_MPI"] = "FALSE"
Since Ultranest can resume from previous runs, the code block below picks up on the samples from the parallelized script.
result = sampler.run(
show_status=True,
update_interval_volume_fraction=0.98,
frac_remain=0.1,
)
[ultranest] Resuming from 11417 stored points
[ultranest] Widening roots to 452 live points (have 400 already) ...
[ultranest] Widening roots to 512 live points (have 452 already) ...
[ultranest] Explored until L=-4e+01 [-39.2230..-39.2198]*| it/evals=10984/3029824 eff=inf% N=400 0
[ultranest] Likelihood function evaluations: 3029824
/Users/abhimatlocal/software/miniforge3/envs/phoebe_py38/lib/python3.8/site-packages/numpy/core/_methods.py:239: RuntimeWarning: overflow encountered in multiply
x = um.multiply(x, x, out=x)
[ultranest] Writing samples and results to disk ...
[ultranest] Writing samples and results to disk ... done
[ultranest] logZ = -64.05 +- 0.1572
[ultranest] Effective samples strategy satisfied (ESS = 2665.6, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.09 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.18, need <0.5)
[ultranest] logZ error budget: single: 0.24 bs:0.16 tail:0.10 total:0.18 required:<0.50
[ultranest] done iterating.
Examine the fit results#
UltraNest offers several outputs and plots to examine the results of our fit. Let’s run those here.
sampler.print_results()
logZ = -63.978 +- 0.274
single instance: logZ = -63.978 +- 0.237
bootstrapped : logZ = -64.051 +- 0.257
tail : logZ = +- 0.095
insert order U test : converged: True correlation: inf iterations
star1_rad : 11.66 │ ▁▁▁▁▁▁▁▁▁▁▂▃▃▃▃▅▅▆▆▆▇▇▆▆▅▄▄▃▂▁▁▁▁▁▁ ▁ │16.07 14.02 +- 0.59
star2_rad : 10.18 │ ▁▁▁▁▁▁▁▁▁▁▁▁▂▂▂▂▃▃▄▅▅▅▆▇▇▇▇▆▇▆▆▅▅▅▄▃▃▃│15.00 13.47 +- 0.80
bin_period : 24.61 │ ▁ ▁▁▁▁▁▂▂▄▃▅▅▅▇▇▇▇▇▇▇▄▄▄▃▂▂▁▁▁▁▁▁▁▁▁▁ │25.59 25.07 +- 0.13
bin_inc : 71 │ ▁▁▇▅▁ ▁▂▇▃▁ │109 89 +- 16
bin_t0 : 53799.664│ ▁ ▁ ▁▁▁▁▁▁▂▂▃▃▅▆▇▇▇▆▆▅▅▄▃▃▂▁▁▁▁▁▁▁▁▁ │53800.249 53799.976 +- 0.067
bin_rv_com : 132.6 │ ▁ ▁▁▁▁▁▁▁▂▂▃▃▄▅▆▆▇▇▇▆▆▆▅▃▃▂▂▁▁▁▁▁▁▁ ▁ │159.3 146.1 +- 3.1
ext_f153m : 4.389 │ ▁▁ ▁▁▁▁▁▁▂▂▂▂▃▅▅▆▇▇▇▇▇▇▆▅▅▄▃▂▂▁▁▁▁▁▁▁ │4.664 4.537 +- 0.036
ext_alpha : 2.166 │ ▁ ▁▁▁▁▁▁▂▃▃▄▄▅▇▇▇▆▆▆▆▅▄▃▂▂▂▂▁▁▁▁▁▁▁▁▁ │2.294 2.225 +- 0.017
sampler.log_to_disk = False
sampler.plot_run()
sampler.plot_trace()
The corner plot may be particularly useful to examine the constraints on the parameters from our fit. We’ll also put in the “truth” values in here (shown as blue lines), corresponding to the value of each parameter with which we generated the mock data.
The fit recovers the generated parameters. For the stellar radii and the center of mass RV, our “truth” value appears to be slightly out of 1 sigma, but still within the 95% confidence bounds. Also note the degeneracy in inclination, which is to be expected. The UltraNest fit recovers a peak at both the mock inclination, $i$, and a peak at $(180\degree - i)$.
ultranest.plot.cornerplot(sampler.results, truths=truth_params,)
