Bayesian Inference with CANDID & Pymask¶

In this notebook, we will show how to go from a calibrated .oifits file to posterior estimates of binary parameters, using two different packages for Bayesian inference:

• CANDID, developed by Antoine Merand & A. Gallenne, and
• Pymask, developed by Benjamin Pope & Anthony Cheetham.

With AMICAL, we provide some easy interface between these codes and the outputs of our extraction pipeline. We give below some example to analyze and extract the quantitative values of a simulated binary.

In [ ]:

from matplotlib import pyplot as plt
import amical

Your inputdata has to be a calibrated oifits file or list of oifits files. We'll use the output of our tutorial for calibrating NIRISS data.

In [ ]:

inputdata = "Saveoifits/example_fakebinary_NIRISS.oifits"

CANDID¶

First we have to initialize the parameters of CANDID.

In [ ]:

param_candid = {
    "rmin": 20,  # inner radius of the grid
    "rmax": 250,  # outer radius of the grid
    "step": 50,  # grid sampling
    "ncore": 1,  # core for multiprocessing
}

Now we apply the AMICAL interface to the CANDID binary grid search.

It will search on a grid defined in circles from rmin to rmax (in milliarcseconds), in stepsize step, and has support for multi-core threading.

If you want to save the figure locally as .pdf, use save=True (new feature June 2021).

In [ ]:

fit1 = amical.candid_grid(
    inputdata, **param_candid, diam=20, doNotFit=[], save=False
)

# Plot and save the fitted model
amical.plot_model(inputdata, fit1["best"], save=False)

We can also estimate contrast limits for non-detections.

To do this in CANDID, we Monte Carlo over the errorbars in the oifits data to estimate, in the absence of bias, the fraction of false positives that would be contributed purely by independent, identically distributed noise. Note that this can be quite different from the results of an end-to-end injection-recovery test, which can probe bias and correlated noise!

In [ ]:

cr_candid = amical.candid_cr_limit(
    inputdata, **param_candid, fitComp=fit1["comp"], save=False
)

Pymask¶

First, we define parameters for Pymask: we need priors in

• separation,
• position angle, and
• contrast

as well as flags for multiprocessing, and additional error.

Pymask proposes to add some extra_error_cp on the CP. This allows to take into account a possibly understimated uncertainties on the data. Indeed, some bias due to mismatch between the calibrator and the science spectral type, or some systematic temporal effect could produce additional errors not properly retrieved by the covariance matrix.

In addition, we can also add some scaling parameter (err_scale) on the CP uncertainties to deal with the number of independent closure phases, which is the number of linearly independent triangles (N(N-1)(N-2)/6), and is smaller than simply the set of triangles ((N-1)(N-2)/2). If you consider the full CP set (35 for a 7 holes mask), you double count your data, so you have to scale your uncertainties by the factor of additional CP, which is sqrt(N/3).

Note that if you used only a subset of CP (by selecting one common hole to save the oifits, see amical.save for details), this additional err_scale is unusable.

In [ ]:

param_pymask = {
    "sep_prior": [100, 180],  # Prior on the separation
    "pa_prior": [20, 80],  # Prior on the position angle
    "cr_prior": [230, 270],  # Prior on the contrast ratio
    "ncore": 1,  # core for multiprocessing
    "extra_error_cp": 0,
    "err_scale": 1,
}

Now we can do a grid search followed by local optimization, which proceeds similarly to CANDID:

In [ ]:

fit2 = amical.pymask_grid(inputdata, **param_pymask)

The main reason to use Pymask is its interface to Bayesian inference tools like nested sampling and MCMC. Here, let's use MCMC from emcee to infer proper Bayesian posteriors for the binary parameters.

In [ ]:

param_mcmc = {
    "niters": 800,
    "walkers": 100,
    "initial_guess": [146, 47, 244],
    "burn_in": 100,
}

fit3 = amical.pymask_mcmc(inputdata, **param_pymask, **param_mcmc)

As with CANDID, we can also Monte Carlo over the errorbars to estimate limits for non-detections - and the same caveat applies that this is no substitute for proper end-to-end injection-recovery simulations:

In [ ]:

cr_pymask = amical.pymask_cr_limit(
    inputdata,
    nsim=500,
    ncore=1,
    smax=250,
    nsep=100,
    cmax=5000,
    nth=30,
    ncrat=60,
)

In [ ]:

plt.figure()
plt.plot(cr_candid["r"], cr_candid["cr_limit"], label="CANDID", alpha=0.5, lw=3)
plt.plot(cr_pymask["r"], cr_pymask["cr_limit"], label="Pymask", alpha=0.5, lw=3)
plt.ylim(plt.ylim()[1], plt.ylim()[0])  # -- reverse plot
plt.xlabel("Separation [mas]")
plt.ylabel(r"$\Delta \mathrm{Mag}_{3\sigma}$")
plt.legend(loc="best")
plt.grid()
plt.tight_layout()