Likelihood

By default the likelihood is determined by the type of transient/data being used. However, users can choose a different likelihood. We note that there is typically only one correct choice of likelihood but there may be edge cases such as errors in time, or non-detections, or uncertain y errors which requires users to use a different likelihood. Many different simple to more complicated likelihoods are included in redback, these should cover most of the cases seen in transient data but if not, users can write their own likelihoods. We encourage users to add such likelihoods to redback. Please check the API for an up-to-date list of the likelihoods available in redback and their usage.

Regular likelihoods

  • Gaussian likelihood - general Gaussian likelihood

  • Gaussian likelihood on magnitude flux (GaussianLikelihoodOnMagnitudeFlux) - accepts AB magnitude data and a model that returns magnitudes, but evaluates the Gaussian likelihood after converting both to flux density

  • GRB Gaussian likelihood - a GRB specific Gaussian likelihood

  • Poisson likelihood - For a poisson process

More advanced likelihoods

  • Gaussian likelihood with additional noise - When you want to estimate some additional uncertainty on your model

  • Gaussian likelihood with uniform x errors - When you have x errors that are bin widths

  • Gaussian likelihood with upper limits (GaussianLikelihoodWithUpperLimits) - A Gaussian likelihood where non-detection data points are treated as upper limits via a CDF. See below for details.

  • Gaussian likelihood with non detections and quadrature noise - Same as above but with an additional noise source added in quadrature

  • StudentT likelihood (StudentTLikelihood) - A StudentT likelihood, useful for data with some outliers, heavier tails than a Gaussian means its less sensitive to outliers. The degrees of freedom nu default to 3; add a prior on nu to sample it alongside the model parameters.

  • MixtureLikelihood - A mixture likelihood with two Gaussian components, assumes each data point either comes from one Gaussian that is consistent with the model or an outlier Gaussian. Provides a probabilistic estimate of each data points probability of being an outlier. Please look at the examples for more details.

Write your own likelihood

If you don’t like the likelihoods implemented in redback, you can write your own, subclassing the redback likelihood for example,

class GaussianLikelihoodKnownNoise(redback.Likelihood):
       def __init__(self, x, y, sigma, function, kwargs):
           """
           A general Gaussian likelihood - the parameters are inferred from the
           arguments of function

           Parameters
           ----------
           x, y: array_like
               The data to analyse
           sigma: float
               The standard deviation of the noise
           function:
               The python function to fit to the data. Note, this must take the
               dependent variable as its first argument. The other arguments are
               will require a prior and will be sampled over (unless a fixed
               value is given).
            kwargs: dictionary of additional keywords for the model
           """
           self.x = x
           self.y = y
           self.sigma = sigma
           self.N = len(x)
           self.function = function

           # These lines of code infer the parameters from the provided function
           parameters = inspect.getargspec(function).args
           parameters.pop(0)
           super().__init__(parameters=dict.fromkeys(parameters))


       def log_likelihood(self):
           res = self.y - self.function(self.x, **self.parameters, **self.kwargs)
           return -0.5 * (np.sum((res / self.sigma)**2)
                          + self.N*np.log(2*np.pi*self.sigma**2))

Gaussian likelihood for magnitude data in flux space

Photometric measurement noise is often approximately Gaussian in flux density rather than in magnitude. GaussianLikelihoodOnMagnitudeFlux allows a transient and model to remain in Redback’s magnitude data mode while evaluating the residuals in AB flux density. The observed magnitudes and model magnitudes are converted to mJy internally, and magnitude errors are propagated as

\[\sigma_f = \frac{\ln 10}{2.5} f_\nu \sigma_m.\]

Construct this likelihood explicitly and pass it to fit_model:

from redback.likelihoods import GaussianLikelihoodOnMagnitudeFlux

likelihood = GaussianLikelihoodOnMagnitudeFlux(
    x=transient.x,
    y=transient.y,
    sigma=transient.y_err,
    function=model,
    kwargs=model_kwargs,
)

result = redback.fit_model(
    transient=transient,
    model=model,
    likelihood=likelihood,
    prior=prior,
    model_kwargs=model_kwargs,
    ...
)

The data and model must use AB magnitudes, and the magnitude uncertainties must be finite and positive. The uncertainty conversion is a local linear propagation, so this likelihood is most appropriate when the supplied magnitude errors were themselves derived from approximately Gaussian flux-density errors. It currently handles detections only; use GaussianLikelihoodWithUpperLimits for datasets containing upper limits.

Non-detections / upper limits

GaussianLikelihoodWithUpperLimits handles datasets that mix genuine detections with non-detection upper limits. Each data point is flagged as a detection or upper limit via a boolean detections array (True = detection, False = upper limit).

For detections the standard Gaussian log-likelihood is used. For upper limits the likelihood contribution is the CDF probability that the true value lies on the correct side of the limit:

  • Flux / flux_density / luminosity data: upper limit means the true value is below the limit, so the contribution is \(\log \Phi\!\left(\frac{y_{\rm lim} - \mu}{\sigma}\right)\).

  • Magnitude data: upper limit means the true value is fainter (larger magnitude) than the limit, so the sign is reversed.

The simplest way to use this is to load simulated data with include_upper_limits=True, which automatically populates the detections array and substitutes limiting magnitudes for NaN values in the non-detection rows:

transient = redback.transient.Transient.from_simulated_optical_data(
    name='my_transient',
    data_mode='magnitude',
    include_upper_limits=True,
    upper_limit_sigma=3.0,   # sigma level of the limits, default 3
)

When fitting, redback automatically selects GaussianLikelihoodWithUpperLimits if the transient has upper limits:

result = redback.fit_model(
    transient=transient,
    model='arnett_bolometric',
    ...
)

You can also construct the likelihood manually:

import numpy as np
from redback.likelihoods import GaussianLikelihoodWithUpperLimits

detections = np.array([True, True, False, True, False])  # False = upper limit
likelihood = GaussianLikelihoodWithUpperLimits(
    x=time, y=magnitude, sigma=mag_err,
    function=my_model,
    detections=detections,
    upper_limit_sigma=3.0,   # sigma level the limits were reported at
    data_mode='magnitude',   # or 'flux', 'flux_density', 'luminosity'
)

Important: upper limit y-values must be finite. For simulated data, the limiting magnitude is substituted automatically. For manually constructed datasets, replace any NaN upper limit values with the actual limiting magnitude or flux before passing to the likelihood.

The measurement-error array must have one entry per observation. Errors for detected points must be finite and positive because they enter the Gaussian likelihood directly. Errors attached to upper-limit points are not used by GaussianLikelihoodWithUpperLimits: their effective uncertainty is calculated from the finite limit value and upper_limit_sigma. Consequently, an upper-limit error may be NaN or infinite without changing this likelihood. A finite placeholder is nevertheless convenient when the same transient will also be passed to plotting code or a custom likelihood that expects every error to be finite.

For example, when constructing a transient from a table with missing upper-limit uncertainties:

detections = np.asarray(detections, dtype=bool)
errors = np.asarray(errors, dtype=float)

invalid_detection_errors = detections & (~np.isfinite(errors) | (errors <= 0))
if np.any(invalid_detection_errors):
    raise ValueError("Detected data require finite, positive uncertainties")

# This value is ignored for rows where detections is False.
errors[~detections] = 1.0

transient = redback.transient.Supernova(
    name="my_supernova",
    data_mode="flux_density",
    time=time,
    flux_density=flux_density,
    flux_density_err=errors,
    bands=bands,
    detections=detections,
    upper_limit_sigma=5.0,
)

For flux-like data, Redback uses \(\sigma_{\rm lim}=y_{\rm lim}/N\) for an \(N\)-sigma limit. For magnitudes it uses \(\sigma_{m,\rm lim}=(2.5/\ln 10)/N\). These are derived inside the upper-limit likelihood; the placeholder in errors does not replace or modify them. Every value supplied through upper_limit_sigma must be finite and strictly positive.

Joint likelihoods

Any likelihood can be combined with another likelihood to form a joint likelihood. This is useful when you want to jointly fit two different types of data. For example, GW and EM data. Or GRB prompt and afterglow data. Or a spectrum and photometry of your favourite transient.