Leeney, Handley & Acedo (2023) proposed a fully Bayesian methodology for simultaneous anomaly flagging and excision, which is performed jointly with the primary model fit to the data, and can be readily folded into Bayesian analysis pipelines via a simple modification to their likelihood. The term likelihood, in the context of a single data point, defines the probability of observing that data point given some model and its constituent parameters . This assumes a predefined probability distribution for the noise. For instance, assuming uncorrelated Gaussian noise results in a likelihood of the form

where is the parameter(s) and is the noise amplitude. Alternatively, for a data set containing data points, and a corresponding model with a single parameter set for all points, assuming no correlations, the overall likelihood is the product of the probabilities of each datum (or alternatively, the sum of log probabilities), giving

where i indexes over the data points. In many cases, this will be an index over a variable x. For the examples in this work, this is taken to be observation frequency.

However, this likelihood cannot account for isolated data points that do not match the probability distribution of the rest of the data set. Throughout this work, we refer to such data points as ‘anomalous’. Contamination by RFI is a common source of such anomalous data points. Typically, anomalies are flagged and excised prior to the Bayesian fitting procedure. This can be problematic, as it leads to potentially useful information being thrown away. Leeney et al. (2023) showed that anomalies can be modelled simultaneously in a Bayesian fashion by using a piece-wise likelihood capable of modelling both the probability of abnormality and the probability that each datum fits to the model of interest.

This was achieved in Leeney et al. (2023) by first defining a model in which, if a data point is contaminated, the probability of observing that data becomes independent of the parameters. It is instead assumed to be uniform over a defined scale , where describes the scale of the contamination, in the same units as those of the data, and can typically be set to (data). The probability of observing a given data point given the model then depends on whether that point is contaminated or not, as

This can be expressed with a Boolean mask to identify each data point as contaminated or not:

Here, indicates uncontaminated data and indicates contaminated data.

As it is not necessarily known which data points are contaminated a priori, a probability can be assigned to each given epsilon mask. This gives the probability of observing the data point, given the model, as

Leeney et al. (2023) then makes the assumption that the probability of data points being contaminated are uncorrelated, with a probability that any individual point is contaminated, such that

and

Substituting in these values into equation (5) gives

Ideally, should be marginalized over. However, this is computationally impractical in most cases. Therefore, Leeney et al. (2023) makes the approximation that the that returns the highest likelihood dominates in the marginalization. This highest likelihood mask is given by

Therefore, the final likelihood can be expressed as

In practice, this likelihood operates by removing very low likelihood data points from the overall likelihood and replacing them with an Occam’s penalty. This is because the Bayesian evidence favours the simplest solution that best describes the data, so without the penalty the ‘most likely’ solution would be to flag all of the data as anomalous.

Leeney et al. (2023) showed that this likelihood enabled anomalous data points to be efficiently flagged and corrected for automatically in a Bayesian model fit of a one-dimensional (1D) data set. The focus of this work is extending this methodology for application to 2D, time-varying data sets.

Anstey, de Lera Acedo & Handley (2023) proposed a methodology for efficient inclusion of time-dependent data sets in Bayesian modelling. Typically, in the case of time-varying data, modelling each time bin separately is unfeasible as it requires a complete model with its own set of parameters for each bin, which can quickly result in the dimensionality of the fit becoming very large for more than a few time bins. As a result, time-varying data are typically modelled by fitting a model to the time averaged data set. For example, in the case of the uncorrelated Gaussian likelihood described in equations (1) and (2), a time-averaged likelihood would take the form

where j indexes time bins and is the total number of time bins.

However, in Anstey et al. (2023), the case was considered where a time varying model can be defined as a product of a parameter-dependent component and a parameter-independent component, where only the parameter-independent component has time dependence, of the form

If it is the case that the time-dependent component , despite being independent of the specific value of the parameters, is required to be different for the different parameters in , which will be the case if different parameters are describing different effects, this can be expressed more fully as

where k indexes over the parameters and expresses the individual elements of the vector .

An example of a case where a model of this form can be defined is for a set of parameters defining an astrophysical observable that is constant on the time-scales of the experiment, with time variance only being introduced by the rotation of the Earth. Another possibility would be a data set of just noise with transient signals, which would simply have a model of , which still satisfies this condition.

It should be noted that time dependence is not the only way a model can be factorized in the manner shown in equation (13). For example, if the frequency component of the model were parameter independent, it could be factorized instead. Provided their effect on the model is not dependent on the parameters being fit for, it may also be possible to factorize polarization, or even different instruments in this fashion. The fast anomaly mitigation techniques described in Sections 2.3 and 2.4 could then be equivalently applied in such cases. However, for the purposes of this work, we will focus on time-dependence as a test case.

If such a model can be defined, it becomes possible to implement simultaneous fitting of separate time bins of data to corresponding models without the dimensionality of the parameters increasing, as each time bin fits for the same parameter set. Such a fitting process can be implemented by modifying the likelihood shown in equation (11) to

It was demonstrated in Anstey et al. (2023), in the context of global 21-cm experiments, that using this full time-dependent likelihood in a Bayesian model fit enabled the time variance of the model to be exploited to constrain the model parameter more tightly than could be achieved for a time-averaged model.

Applying the Bayesian anomaly mitigation technique described in Section 2.1 to the time-dependent likelihood described in Section 2.2 is straightforward, requiring only that the flagged likelihood shown in equation (10) be extended into two dimensions as

where is the likelihood of a single data point in a single time bin

However, implementing this method in practice faces a challenge. The likelihood shown in equation (14) requires a summation over time bins. This means the computation time of the likelihood grows linearly with the number of time bins used in the data set. As a result, larger numbers of time bins can greatly slow the fitting procedure.

In Anstey et al. (2023), this issue was resolved by reformatting the likelihood such that all summations over time could be calculated once, outside of the likelihood, removing the dependence of the calculation time on the number of time bins. However, this solution is not possible to implement when the anomaly correcting procedure is also implemented. This is because the value of must be calculated for every time bin within the likelihood, in order to evaluate , as shown in equation (15). Therefore, an alternative method of speeding the likelihood evaluation is needed to make this process viable in practice. This can be achieved using likelihood reweighting.

The process of likelihood reweighting is an extension to importance sampling, pioneered in the context of gravitational waves (Payne, Talbot & Thrane 2019; Romero-Shaw, Lasky & Thrane 2019). It is a method for speeding the evaluation of a posterior and evidence in a Bayesian fit for the case of a complex model that is otherwise slow to evaluate. This process relies on several key criteria.

Firstly, two models, and , are required. One of these must be quick to evaluate, which will henceforth be assumed to be model F, and one is slower to evaluate, which we define as model S. These two models must be parametrized by the same parameter vector , with the same prior distribution . They must also have their posterior peak in approximately the same region of parameter space.

By definition, the posteriors of the two models can be expressed as

and

where and are the likelihoods calculated from the two models and and are the respective evidences.

However, given the aforementioned criterion that the two models have the same priors, the prior can be expressed in terms of model F, as

and substituted into equation (18) to give

Thus, the posterior of model S can be evaluated without having to perform a full Bayesian model fit with the slow-to-calculate likelihood. This is achieved by instead performing a model fit of the much faster evaluated model F. Given the criterion that the bulk of the two models’ posteriors occupy similar regions of the parameter space, the samples of this fast evaluated posterior will cover the same parameter volume as a hypothetical posterior of the slower model. Therefore, reweighting the posterior samples of model F by a factor of the ratio of the likelihoods

will convert the samples to samples of the posterior of model S, to within a constant factor given by the ratio of the evidence. By this method, a set posterior for model S can be evaluated with the slow likelihood only needing to be computed for the relatively small number of posterior samples. This can be significantly faster than sampling the complex likelihood across the entire prior volume. For a particularly complex and thus slow to calculate likelihood, this can make the fitting procedure significantly faster. This effect will be demonstrated in Section 3.3.3.

Therefore, the process of likelihood reweighting is well suited to achieving fast and efficient time-dependent anomaly mitigation. By taking an anomaly correcting likelihood that acts on precomputed time-summed likelihoods as the quickly evaluated model F and the full time separated likelihood as the slowly evaluated model S, the full time-dependent likelihood posterior can be evaluated quickly, allowing the Bayesian anomaly mitigation procedure to be implemented practically on time-dependent data sets as follows. For the slow likelihood, the likelihood of each data point is calculated according to equations (15) and (16).

For the fast likelihood, however, in order to preserve as much of the time-dependent data as possible, the following process is implemented.

First, the methodology described in Anstey et al. (2023) to evaluate a time-dependent likelihood in a time-independent fashion by expanding the likelihood given in equation (14) and precomputing the summations over time is implemented. As previously noted, precomputing time summations in this fashion, while greatly accelerating the likelihood evaluation time, prevents the evaluation of directly and thus prevents the implementation of the full time-dependent flagging method defined in equation (15). However, from equation (7) of Anstey et al. (2023), it is possible to evaluate the product of likelihoods (or equivalently the sum of log likelihoods) over all time bins in each data channel

from precomputed time summations.

Therefore, it is possible to perform the anomaly mitigation methodology in a time-independent fashion by flagging out entire channels where contamination occurs, based on the summed log likelihoods, rather than only the contaminated time bins of those channels.

In order to determine the appropriate flagging thresholds for this process, we repeat the derivation of Leeney et al. (2023), described in Section 2.1, applied to the summed log likelihood.

We first make the assumption that if any time bin in a particular channel is contaminated, that contamination will dominate, and the resulting probability of the product of all time bins becomes independent on the model, and is instead uniform over a defined scale. In this case, the resulting probability of a contaminated channel depends on how many contaminated points it contains:

where again describes the scale of the contamination, is the number of contaminated points in channel i, and .

This can again be expressed using a Boolean mask, , as

Here, refers to the case where no time bins in channel i are contaminated, and to the case where at least 1 time bin is contaminated.

As before, the probability of observing the data channel, given the model, can then be expressed as

It must be noted here that, as the value of in each case is not known, we sum over all possible values, weighted by their probabilities.

If, again, the probability that any given point is contaminated is , assuming again that there are no correlations, the probability that a channel contains no contaminated points is

where is the number of time bins.

In the case of , the probability of points being contaminated is a simple binomial distribution:

As the case of is defined to be where at least one point in the channel is contaminated, so , it can be seen that

for all .

Substituting equations (26)–(28) into equation (25) and noting that

gives

Making, as before, the assumption that the most likely model dominates the marginalization of the masks, and given that the most likely mask is given by

the final likelihood can be expressed as

It can be seen that in the case of , this reduces to equation (10), as would be expected.

Overall, this likelihood uses the methodology described in Anstey et al. (2023) to evaluate a time-dependent likelihood in a time-independent fashion, while simultaneously implementing the anomaly mitigation methodology in a time-independent fashion by flagging out entire channels where contamination occurs, rather than only the contaminated time bins of those channels. This results in a likelihood that takes the same parameters as the full flagging likelihood shown in equation (15), with a posterior peak expected in a similar position, but is much faster to evaluate, albeit with a lower constraining power due to utilizing a smaller fraction of the available data. Therefore, equation (32) serves as an ideal fast likelihood for a likelihood reweighting process, as required.

In the following sections, this complete process will be tested on simulated data to evaluate its efficacy.

In order to evaluate the efficacy of this method for correcting transient RFI in time separated data, a toy-simulated time-dependent data set was generated of the form

where x describes the data channel value, indexed by i. For this toy model, this is assumed to be frequency in MHz and runs from 50 to 200. , , , and are time-dependent variables that are chosen to vary gradually over the time bins. For each of the four values, a start point is randomly chosen uniformly from the range [0–5] for , [0–1] for and , and [110–130] for . The four variables are then iteratively defined according to

where indicates a value randomly drawn from a normal distribution of mean and standard deviation . and were both set to 0, to ensure that the absolute scale of the toy model did not significantly vary from time bin to time bin, with standard deviations of and . and were both randomly drawn from the range [−0.05 – 0.05], with , such that the sinusoidal distortion would vary over time.

is a realization of random Gaussian white noise added to the data. Unless otherwise specified, this was set to have a standard deviation of 0.25. A different noise realization is added to each time bin.

A toy model of this form was chosen as it approximates the form of data from a global 21-cm experiment, which is dominated by foregrounds of primarily diffuse synchrotron emission (Shaver et al. 1999). The power law with a 2.55 spectral index approximates the spectral variation of the diffuse emission from the sky and the time varying sinusoids approximate the convolution of the diffuse emission with a chromatic antenna beam, as the Earth rotates. This is described, e.g. in equation (18) of Anstey et al. (2021) and seen in Bowman et al. (2018). Global 21-cm experiments are an anticipated use case of this process, so a toy model of this form enables these experiments to be used as a test case of the methodology.

In addition, a simulated data set of this form has the required structure to apply the time-separated model fitting as specified in equation (12), with

and

with a ‘true value’ of .

Fig. 1 shows an example simulated data set generated for 20 time bins using this toy model. Once such a data set has been generated, any arrangement of anomalous points can then be added in order to test the proposed Bayesian anomaly mitigator.

In Leeney et al. (2023), the likelihood threshold value for determining if a point should be flagged as an anomaly or not, p, was set to a fixed value. However, doing so produces a challenge when applied to time-dependent modelling.

As shown in equation (15), the general method of operation for the Bayesian anomaly correction procedure is to calculate the likelihood of each separate data point, compare it with a predefined threshold value and if it exceeds the threshold, include it with an appropriate weighting and if it does not, flag it as an anomaly and instead add a fixed penalty to the total likelihood. This is acceptable under the assumption that only anomalous points will have likelihoods lower than the threshold.

However, in cases where the model used has relatively large priors and significant variability, it becomes possible for a given parameter sample to produce a model sufficiently different from the data set that significant numbers of data points have likelihoods below the threshold, as demonstrated by Fig. 2.

This becomes of concern if the priors are wide enough compared with the posterior peak to result in sections of the prior space in which every data point has a likelihood below the defined threshold set by p. In this case, the overall likelihood is a sum of only fixed penalty terms, and thus is constant. As a result, all sections of the prior space that satisfy this condition have the same likelihood. This manifests as the outer boundaries of the likelihood surface becoming flat.

If the majority of the prior space is not flattened in this manner, this has minimal impact on the process of performing a model fit and recovering the parameter posteriors. However, if the likelihood surface is particularly steep, it is possible for the majority of the space to become flat, except for a very narrow region around the posterior peak. In this case, performing a model fit becomes impossible, as there exists no variation to guide an algorithm towards the peak.

This is a challenge when performing time-dependent model fits, because the additional time information amplifies the effect of model disagreement. A parameter sample that gives a model that matches well to a data set with many time bins will give a high likelihood for each data point of each time bin, and thus give a larger overall likelihood than for an equally well matching data set with fewer time bins, and vice versa. As demonstrated in Fig. 3 for the proposed one-parameter toy model, defined by equations (36) and (35), with 2, 20, 200, and 2000 time bins, this has the effect of steepening the likelihood surface, the more time bins the data set has.

Fig. 4 shows the regions where the value of , averaged over time, falls above or below a threshold of 1e−3, which is a typical value as determined from Leeney et al. (2023). It can be seen that as the number of time bins increases, increasingly large portions of the space fall below the threshold, even in absence of anomalies, and so are flattened. By 2000 time bins, only a very small region of the parameter space is not flattened, which makes performing a model fit almost impossible.

This demonstrates that, especially in the case of time separated data, having an externally fixed threshold value can impede model fitting. Therefore, ideally the value of p should be dynamic, allowing it to be low in suboptimal regions of the parameter space to avoid overflagging and flattening the likelihood surface, and higher around optimal regions to avoid missing genuine anomalous points.

Therefore, as was suggested in Leeney et al. (2023), this issue can be resolved by fitting the value of p as a free parameter, which is estimated as a parameter, simultaneously alongside other parameters of interest by the chosen Bayesian numerical solver, rather than assigning it a fixed value. It is assigned a wide prior from its theoretical maximum of one down to effectively zero. This results in there being sections of the parameter volume where the threshold is low and so the variations in the likelihood surface with the other parameters are visible. As a result, the model fit can progress, optimizing the parameter values towards their posteriors while simultaneously optimizing the threshold towards the optimal posterior value of the probability that a point is anomalous for that data set. This resolves the issue described above and allows the fit to proceed while still accurately flagging anomalous points. The effect this additional nuisance parameter may have on the likelihood reweighting process can then be mitigated by marginalizing over the parameter before reweighting.

In the next section, tests of the entire process will be performed to demonstrate its efficacy. Throughout the rest of this work, p is fit as a free parameter with a logarithmically uniform prior in the range [1e−30 – 1]. The next section will demonstrate the functionality of this method.

In order to demonstrate the performance of this anomaly mitigation procedure, four test data sets were generated according to the model described in Section 3.1, with = 2, 20, 200, and 2000 time bins, respectively. These data sets will henceforth be referred to as the uncontaminated data. To each of these data sets, a random selection of anomalous peaks is added. In each case, the number of spikes added was equal to , such that every data set is contaminated in equal proportion. The amplitudes were uniformly sampled from the range 10–50. This scale was chosen to provide a trial case for the process, in which the anomalies are sufficiently above the noise to be large enough to impact the model fit, but still relatively small compared with the overall scale of the data, so they cannot be trivially removed. The time bins in which each anomaly was placed were uniformly randomly chosen and the x bin was uniformly set to one of 40 randomly chosen channels. This constrained the contamination to a subset of the data channels, in a manner more resembling RFI. This will be discussed in more detail in Section 3.3.4. Fig. 5 shows an example of the anomalous points added to the data for . These data sets will henceforth be referred to as contaminated. As the primary objective of this work is to present a method of fast correction for transient RFI, the contamination we test here is transient. However, this method will also compensate for time-constant and wideband RFI, which was demonstrated in Leeney et al. (2023).

In all model fits performed in this paper, is taken to be 7000, which is the maximum scale of the toy model data in all cases, as can be seen in Fig. 1. Furthermore, all fits were performed using the PolyChord nested sampling algorithm (Handley, Hobson & Lasenby 2015a, b), implemented in Python and run on a laptop on one CPU. As the core of this method is a modification to the likelihood and model of a Bayesian model fit; however, it can be implemented with any Bayesian fitting algorithm and will scale and parallelize as that algorithm does.

In order to test the ability of the method to correctly identify small transients in large data sets, a series of tests were run in which a single anomalous data point was added at random to the previously described 2000 time bin test data set, and the model fit with the flagging process implemented. This was repeated five times each for anomalous points with amplitudes of 0.5, 1.25, 2.5, 5, and 12.5, which correspond to SNRs of 2, 5, 10, 20, and 50 respectively.

The results of these fits are shown in Fig. 12. It can be seen that when acting as a flagger, this method correctly and reliably identifies the single anomalous point for the cases with SNR . Furthermore, it never erroneously identifies any data points as anomalous that were not anomalous in the data. However, the flagger does begin to fail to detect the anomalous point at lower SNRs, which can be expected given there is only one anomalous point out of 302 000 in this test case, meaning lower SNR anomalies can become indistinguishable from simple statistical fluctuations and so will not be flagged.

Although real transients such as FRBs and pulsars will often have additional structure and cover more data points than this simple, single contaminated data point test case, this minimal case demonstrates the potential for the proposed methodology to function as an efficient Bayesian transient detector. This will be explored in greater detail in a future work.