Fitting

After downloading/simulating data, creating a transient object, specifying a model, and creating a prior we now come to the exciting part; Fitting the model to data!

To fit our model to data we have to specify a sampler and sampler settings. The likelihood is set by default depending on the transient/data but one can use a different one or write their own as explained in the likelihood documentation.

Installing redback with minimal requirements will install the default sampler dynesty. Installing optional requirements will also install nestle. We generally find dynesty to be more reliable/robust but nestle is much faster. We note that dynesty has checkpointing, as do many other samplers.

Samplers

As we use bilby under the hood, we have access to several different samplers. Cross-checking results with different samplers is often a great way to ensure results are robust so we encourage users to install multiple samplers and fit with different samplers.

Nested samplers

Dynesty:
Nestle
CPNest
PyMultiNest
PyPolyChord
UltraNest
DNest4
Nessai

MCMC samplers

bilby-mcmc
emcee
ptemcee
pymc3
zeus

A full up to date list of samplers can be found in the bilby documentation. This page also provides guidance on how to install these samplers, while the bilby API provides information on the sampler settings for each sampler.

Fitting with redback

In redback, having created a transient object, specified a model, priors, fitting is done in a single line.

result = redback.fit_model(transient, model=model, sampler='dynesty',
                           nlive=200, transient=afterglow, prior=priors, **kwargs)

Here

transient: Is the transient object created we want to fit
model: is a string referring to a function implemented in redback. Or a function the user has implemented.
sampler: is a string referring to the sampler. It could be a string referring to any name of a sampler implemented in bilby.
nlive: is the number of live points to sample with. Higher = better. Typically we would use nlive=1000/2000 but this depends on the sampler.
transient: the transient object
prior: the prior object
data_mode: type of data to fit.
kwargs: Additional keyword arguments to pass to fit_model, such as the likelihood, or things required by the sampler, label of the result file, directory where results are saved to etc.

Please see the bilby documentation for more information on the sampler settings. As well as the redback API.

Fitting models with extinction, phase, or additional effects

In general most redback models work on the assumption that the time provided to the model is a time since the explosion/burst etc. I.e., time = 0 is when the transient starts. However, sometimes users will not know this and the time they will have is the times of observation in MJD or some other time system.

In this case, we must ensure both the model and the transient object are set up correctly. In particular, you must ensure that the transient object is set up with time_mjd as an attribute instead of time and that the model is set up to take t0 as an input.

# Create a transient object
sn = redback.transient.Supernova(name=name, time_mjd=time_mjd, magnitude=flux_density,
                                magnitude_err=mag_err, bands=bands, use_phase_model=True)

# Create a model
model = 't0_base_model' # This model is a general workhorse which just adds a t0 to the underlying model which we set by
base_model = 'arnett' #this could be any other redback model.

# Create a prior
priors = redback.priors.get_priors(model=base_model)

# we must add a prior for t0
priors['t0'] = bilby.core.prior.Uniform(minimum=sn.x[0]-100, maximum=sn.x[0]-1 name='t0', latex_label=r'$t_{0}~\mathrm{MJD}$')

# We must also make sure the model kwargs not include a pointer to the base model we want to use.
model_kwargs = {bands: sn.filtered_bands, output_format:'magnitude', base_model: base_model}

This is just one such example of a base model, there are models which include extinction, or additional physical effects. Please look at the API for more information.

Fitting bolometric luminosities of optical transients

Most redback models can also output bolometric luminosities of optical transients, which is often returned in erg/s. However, the transient objects assume luminosity in units of 10^50 erg/s. A simple workaround for this is to write a small wrapper function that takes the luminosity in erg/s from the function and divides by 10^50. And you then use this wrapper function to fit instead. This is shown in the examples,but we also provide a simple example here.

def luminosity_wrapper(tt, **kwargs):
    func = redback.model_library.all_models_dict['arnett_bolometric'] # or any other model
    luminosity = func(tt, **kwargs)
    return luminosity / 1e50

result = redback.fit_model(name='GRB', model=luminosity_wrapper, sampler='dynesty',
                           nlive=200, transient=afterglow, prior=priors,
                           data_mode='luminosity', **kwargs)

Fitting with your own/different likelihood

For users which wish to change the likelihood, we provide an easy way to do this. Many likelihoods are implemented but we can also write our own likelihood. Once written this likelihood can be passed to the fit_model function as follows:

result = redback.fit_model(name='GRB', model=model, sampler='dynesty',
                           nlive=200, transient=afterglow, prior=priors,
                           data_mode='luminosity', likelihood=likelihood, **kwargs)

Please look at the documentation for how to use the likelihoods correctly.

Using MPI

We note that some samplers have multiprocessing, which you can see how to use here.

To use MPI with redback, most samplers will work out of the box and all you need to do is pass a npool (or similar argument, please see the bilby docs) to the fit_model function. This is especially true for linux, and the general scripts we provide in the examples folder will work out of the box as soon as you pass npool or the relevant argument. For Mac users, there is some finicky behaviour with MPI which requires that you wrap the fit_model function in a __main__ block and to add a line that explicitly sets the multiprocessing start method to “fork” at the begining of your script.

This is a known issue with MPI on Mac and not specific to redback or bilby.

An example of how to do this is shown below:

import bilby
import redback
import multiprocessing
multiprocessing.set_start_method("fork", force=True)

# some code here

if __name__ == "__main__":
    result = redback.fit_model(name='GRB', model=model, sampler='nessai',
                           nlive=1000, transient=afterglow, prior=priors,
                           n_pool=4, likelihood=likelihood, **kwargs)

We will soon implement some GPU models and JAX functionality for more rapid inference workflows.