# Fully Adaptive Bayesian Algorithm for Data Analysis
FABADA

Pablo M. Sánchez Alarcón, Yago Ascasibar Sequeiros

###### Abstract

The aim of this paper is to describe a novel non-parametric noise reduction
technique from the point of view of Bayesian inference that may automatically
improve the signal-to-noise ratio of one- and two-dimensional data, such as
e.g. astronomical images and spectra. The algorithm iteratively evaluates
possible smoothed versions of the data, the smooth models, obtaining an
estimation of the underlying signal that is statistically compatible with the
noisy measurements. Iterations stop based on the evidence and the $\chi^{2}$
statistic of the last smooth model, and we compute the expected value of the
signal as a weighted average of the whole set of smooth models. In this paper,
we explain the mathematical formalism and numerical implementation of the
algorithm, and we evaluate its performance in terms of the peak signal to
noise ratio, the structural similarity index, and the time payload, using a
battery of real astronomical observations. Our Fully Adaptive Bayesian
Algorithm for Data Analysis (FABADA) yields results that, without any
parameter tuning, are comparable to standard image processing algorithms whose
parameters have been optimized _based on the true signal_ to be recovered,
something that is impossible in a real application. State-of-the-art non-
parametric methods, such as BM3D, offer a slightly better performance at high
signal-to-noise ratio, while our algorithm is significantly more accurate for
extremely noisy data (higher than $20-40\%$ relative errors, a situation of
particular interest in the field of astronomy). In this range, the standard
deviation of the residuals obtained by our reconstruction may become more than
an order of magnitude lower than that of the original measurements. The source
code needed to reproduce all the results presented in this report, including
the implementation of the method, is publicly available at
https://github.com/PabloMSanAla/fabada.

###### Index Terms:

Image denoising, Bayes’ methods, Data analysis, Astronomy.

##  I Introduction

The acquisition of any kind of experimental data is affected by several
sources of statistical error, which ultimately translate into a random noise
component in the measurements to be recorded. There are different types of
noise depending on their physical origin, both related to electronic (thermal
noise, fluctuations) and mechanical (non-perfect lenses, antennas, etc.)
devices. In astronomy, for example, errors can be produced in the acquisition
of the images due to defects in the optics of the telescopes and in the read-
out process of the detector (typically a CCD) transforming the light captured
by the telescope into an electrical signal. The noise introduced can sometimes
be comparable to or even larger than the signal, and different image
processing algorithms may be used to recover the information that is buried
deep in the data.

Smoothing, where measurements are weighted at nearby spatial or temporal
points (using different schemes to assign weights), is one of the most popular
techniques to mitigate the effects of random noise [1, 2]. Nowadays there is a
large number of smoothing algorithms, based on many different techniques, see
[3] for a review, such as central moving average, data grouping/segmentation
[4], fitting smooth functions, different types of statistical analysis,
partial differential equations, wavelength transformation filters, linear and
nonlinear filtering, sparse models and nonlocal self-similarity models [5] and
more recently artificial neural networks, see [6, 7]. Rest of these methods
rely on some explicit or implicit assumptions about the true (noise-free)
signal in order to separate it properly from the random noise. A common
assumption is that the signal being retrieved varies gradually and that the
data can be fit by a smooth function, see [8].

Many techniques analyze the probability that the data correspond to a random
Gaussian realization of the model that attempts to describe the underlying
signal plus random fluctuations of known amplitude [9]. In this work we use
Bayesian inference to evaluate and combine different candidate models that
iteratively attempt to improve the quality of the fit to the data. This new
Bayesian technique incorporates an automatic selection criterion based on the
statistical properties of the residuals, and therefore it yields a fully non-
parametric method. Although the motivation of our scheme is the application in
the field of astronomy, our new algorithm, Fully Adaptive Bayesian Algorithm
for Data Analysis FABADA is focused in an general way, and it is possible to
generate a smooth model for any type of data. Several algorithms have been
developed over the years for denoising data, and their ability to recover the
underlying signal from an experimental data set with its corresponding errors
has been evaluated and compared according to different standard metrics. Our
method is fully described in Section II, and we use a set of synthetic tests
to compare it with other prescriptions in the literature. Details of the
comparison procedure are provided in Section III, and results are discussed in
Section IV. Our main conclusions are succinctly summarized in Section V.

##  II FABADA

The goal of our algorithm, FABADA is to estimate an unknown signal $y=y(x)$ at
$N_{X}$ different locations $x\in X$ specified by a $1\mathrm{D}$ or
$2\mathrm{D}$ coordinate that belongs to the data domain. Its input are
$N_{X}$ independent measurements, contaminated with random Gaussian white
noise $\eta(x)\sim\mathcal{N}(0,e(x))$, and the associated errors $e=e(x)$.
The noisy observational data, $z:X\rightarrow\mathbb{R}$, have the form
$z(x)=y(x)+\eta(x)$, and FABADA returns an estimation $\hat{y}=\hat{y}(x)$ of
the original signal $y(x)$ that is statistically compatible with the
measurements $z(x)$.

For this purpose we apply the Bayes’ theorem in an iteratively way to generate
smooth models of the noisy measurements. Therefore, we must define a suitable
likelihood function to evaluate these models to be tested, and specify a prior
probability distribution for the signal $y$ to initiate the process. Our
likelihood is based on the statistics of a Gaussian process, and we start from
an improper, constant prior. Then, we evaluate different smooth versions of
the posterior probabilities until a certain condition is reached, and we
combine all the smooth models to produce the final estimation $\hat{y}(x)$,
taking into account both their Bayesian evidence and the $\chi^{2}$ of the
residuals.

###  II-A Iterative models

|  Symbol | Definition
---|---|---
DIMENSIONS |  $x\in X$ | Coordinates of the data set (either 1D or 2D coordinates)
|  $i$ | Iteration number $\\{1,2,,..,N_{i}\\}$
INPUTS |  $I=\\{z,e\\}$ | Set of inputs of the algorithm
|  $z=z(x)$ | Observational measurements
|  $e=e(x)$ | Errors associated to the measurements
MODELS |  $m_{i}=m_{i}(x)$ | Set of all possible models with shape $(N_{X},N_{i})$.
|  $\mu_{i}=\mu_{i}(x)$ | Expectation value of the element $x$ of the model $m_{i}$
|  $v_{i}=v_{i}(x)$ | Associate variance of the expectation value
|  $\mathcal{L}(z|m_{i},e)$ | Likelihood of the model $m_{i}$
|  $\mathcal{E}(m_{i})$ | Evidence of the model $m_{i}$
PRIORS |  $p(m_{i})$ | Prior probability distribution of the model $m_{i}$
|  $\tilde{\mu}_{i}$ | Central moving mean of $\mu_{i}$ along $N_{avg}$ neighbors
POSTERIORS |  $\mathcal{P}(m_{i}|z,e)$ | Posterior probability distribution of the model $m_{i}$
OUTPUT |  $\hat{y}=\hat{y}(x)$ | Estimation of the underlying signal
TABLE I: List of all symbols used to describe the FABADA algorithm in Section
II.

FABADA, as the name suggests, is a fully automatic algorithm that only takes
the data set $z=z(x)$ and its associated errors $e=e(x)$ as input,
$I=\\{z,e\\}$, where the domain $X$ may have one or two dimensions. Within the
field of astrophysics, $y$ may be, for example, the spectral energy
distribution as a function of wavelength $x$, or a broadband photometric image
of an arbitrary region of the sky. This is indeed the kind of data that we
will use to illustrate the performance of the algorithm, but the results can
be easily generalized to any other type of empirical measurement in one or two
dimensions affected by Gaussian random errors.

Under this assumption, the likelihood function $\mathcal{L}$ for a model
$m(x)$ of the underlying signal $y$ is given by

$p(z|m,e)=\mathcal{L}(z|m,e)=\prod_{x\in
X}\frac{e^{-\frac{[z(x)-m(x)]^{2}}{2\cdot e(x)^{2}}}}{\sqrt{2\pi}\cdot e(x)}.$
(1)

using the notation described in Table I. Since we do not assume any previous
knowledge about the signal, the prior probability distribution for our initial
model $m_{0}$ will be homogeneous in the range of all possible values, i.e.
$p(m_{0}(x))=1$ for all $x$, during the first iteration of the algorithm.
According to Bayes’ theorem, the posterior probability distribution

$\mathcal{P}(m|z,e)=\frac{p(m)\cdot\mathcal{L}(z|m,e)}{\mathcal{E}}$ (2)

is in this case the straightforward multivariate Gaussian

$\mathcal{P}(m_{0}|z,e)=\prod_{x\in
X}\frac{e^{-\frac{[z(x)-m_{0}(x)]^{2}}{2\cdot e(x)^{2}}}}{\sqrt{2\pi}\cdot
e(x)}$ (3)

centered on the empirical measurements $z$. The expected value of $m_{0}(x)$
is thus $\mu_{0}(x)\equiv\langle m_{0}(x)\rangle=z(x)$, its variance
$v_{0}(x)=e^{2}(x)$ is determined by the corresponding errors, and the
Bayesian evidence of the model $\mathcal{E}=\int p(m)\,\mathcal{L}(z|m,e)\
{\rm d}m$ reduces to unity.

To create smoothed versions of this first model, we iteratively apply a
central moving average filter

$\tilde{\mu}_{i}(x)=\frac{\sum_{x\in\rm avg}\mu_{i}(x)}{N_{\rm avg}}$ (4)

to the expected values $\mu_{i}(x)$ of the last iteration. We simply adopt
$N_{\rm avg}=3$ in 1D (including the two adjacent measurements) and $N_{\rm
avg}=9$ (a $3\times 3$ square) in 2D. The basis of our method is to use
information from neighboring points to update our priors regarding the correct
value of $m(x)$. For every iteration $i>0$, the prior probability distribution

$p(m_{i}(x))=\frac{e^{-\frac{[\tilde{\mu}_{i-1}(x)-m_{i}(x)]^{2}}{2\cdot
v_{i-1}(x)}}}{\sqrt{2\pi\cdot v_{i-1}(x)}}$ (5)

becomes a Gaussian centered on the smoothed expectation (4) of the posterior
distribution from the previous iteration, using its local variance $v_{i}(x)$
as a measure of our uncertainty.

We stress that we are forsaking the strict Bayesian philosophy (the prior
distribution should be, as its name indicates, totally independent from the
data) on purely practical grounds. Once we accept this premise, one may
compute the posterior probability distribution

$\mathcal{P}(m_{i}|z,e)=\frac{e^{-\frac{[\mu_{i}(x)-m_{i}(x)]^{2}}{2v_{i}(x)}}}{\sqrt{2\pi
v_{i}(x)}},$ (6)

where

$\mu_{i}(x)=v_{i}(x)\cdot\left[\frac{\tilde{\mu}_{i-1}(x)}{v_{i-1}(x)}+\frac{z(x)}{e^{2}(x)}\right]$
(7)

and

$\frac{1}{v_{i}(x)}=\frac{1}{e^{2}(x)}+\frac{1}{v_{i-1}(x)},$ (8)

as well as the evidence

$\mathcal{E}_{i}(x)=\frac{e^{-\frac{[\tilde{\mu}_{i-1}(x)-z(x)]^{2}}{2[e^{2}(x)+v_{i-1}(x)]}}}{\sqrt{2\pi[e^{2}(x)+v_{i-1}(x)]}}.$
(9)

As long as the model remains relatively close to the data (within an
environment of the order of $\sqrt{e^{2}(x)+v_{i}(x)}$), the evidence in its
favour will be high, but if it departs significantly, the exponential term
will indicate that we have reached the maximum smoothing that is statistically
compatible with the data.

###  II-B Stopping criteria

On the one hand, we evaluate the average evidence of model $m_{i}$ as

$\langle\mathcal{E}_{i}\rangle=\frac{1}{N_{X}}\sum_{x\in X}\mathcal{E}_{i}(x)$
(10)

and we ensure that we iterate until it reaches a maximum;
$\langle\mathcal{E}_{i}\rangle<\langle\mathcal{E}_{i-1}\rangle$. Since our
priors are based on the data themselves, it is not surprising that the
evidence increases rapidly during the first few iterations of the algorithm.
However, our model uncertainties $v_{i}(x)$ decrease monotonically, and the
smooth models quickly become less compatible with the data than the first
(self-tuned) estimates. Thus, the average evidence reaches its maximum at a
very early stage, and then it slowly declines as the number of iterations
increase.

In order to overcome the bias arising from the lack of independence of the
priors (i.e. overfitting), we made use of the chi-square statistics

$\chi^{2}_{i}=\sum_{x\in X}\frac{[\mu_{i}(x)-z(x)]^{2}}{e^{2}(x)}$ (11)

and iterate until $\chi_{i}^{2}>N_{X}-2$ (the absolute maximum of the
associated probability density function). Note that this condition is achieved
when, on average, the model has departed by about one sigma from the
observational measurements.

We always impose _both_ criteria, so that the algorithm stops on the largest
number of iterations. Usually this is set by the chi-square statistics,
although the evidence-based condition may dominate at very high signal-to-
noise ratio.

###  II-C Model selection

Once the algorithm stops, after $N_{i}$ iterations, we combine all the models
$m_{i}(x)$ to generate our final estimation $\hat{y}(x)$ of the real signal
$y(x)$ as a weighted sum over the expected values

$\hat{y}(x)=\frac{\sum_{i=0}^{N_{i}}w_{i}(x)\,\mu_{i}(x)}{\sum_{i=0}^{N_{i}}w_{i}(x)}$
(12)

where the adopted weights $w_{i}(x)$ are an important ingredient of our
algorithm.

We explored several phenomenological prescriptions in an attempt to optimize
the performance, stability, and reliability of the results. On the one hand,
our tests show the convenience of taking into account the local evidence
$\mathcal{E}_{i}(x)$ calculated at each iteration for every location $x$. This
quantity adapts to the local structure of the data, giving more weight to the
smoother models (larger number $i$ of iterations) in areas where the data are
indeed smooth, while preserving the information when sharp edges are present.

On the other hand, we use the overall $\chi^{2}_{i}$ statistic of each model
$m_{i}$ to gauge its compatibility with the measurements $z$, given the errors
$e$. The probability density function of the $\chi^{2}$ of a large number of
random variables can be approximated as a narrow Gaussian centered around
$N_{X}-2$. Since our priors are biased because they are based on the data
themselves, the maximum evidence occurs at values of $\chi^{2}$ that are
typically much lower than this value (i.e. they overfit the measurements).
Using the probability density function directly would give minimal weight to
any model that is not extremely close, $\chi^{2}\simeq N_{X}-2$, which we find
tends to yield overly smoothed models. In order to avoid this problem, we use
the actual value of $\chi^{2}$ to give more weight to the smoother models, but
not so much that the models favored by the Bayesian evidence are almost
completely ignored. We thus adopt

$w_{i}(x)=\mathcal{E}_{i}(x)\ \chi^{2}_{i},$ (13)

as our final prescription for $i>0$.

Of course, this expression is not valid for $i=0$, since $\chi^{2}_{0}=0$ and
the initial evidence is $\mathcal{E}_{0}(x)=1$ for every point. Note that, due
to the use of an improper prior, this value is in general very different from
$\mathcal{E}_{1}(x)$. To avoid these problems, we a posteriori set
$\chi^{2}_{0}=\chi^{2}_{1}$ based on the first iteration (i.e. the lowest of
all smooth models) and use

$\mathcal{E}_{0}(x)=\frac{e^{-\frac{1}{2}}}{\sqrt{2\pi}e},$ (14)

which is merely a reformulation of (9), assuming that our initial model based
on the measurements should be, on average, about one sigma away from the true
signal. Thus,

$w_{0}(x)=\frac{\chi^{2}_{1}\ e^{-\frac{1}{2}}}{\sqrt{2\pi}e}.$ (15)

##  III Synthetic tests

We develop a battery of synthetic tests based on real astronomical data to
assess the performance of the algorithm. More precisely, we apply FABADA to a
set of astronomical spectra and images, with different levels of Gaussian
random noise, and compare the quality of the reconstructed signal (in terms of
the peak signal-to-noise ratio PSNR and the structure similarity index measure
SSIM), as well as the execution time, with other methods available in the
literature.

###  III-A Other Algorithms

Over the last decades, lots of effort have been placed into the development of
several applications that help in the analysis of digital images in different
fields. Noise reduction is one of the basic problems in this context, and we
attempt to provide a fair comparison of our algorithm with other methods that
are representative of the current state of the art. Leaving aside the
techniques based on some kind of training, such as e.g. neural networks [6, 9,
7], we focus on a more traditional, statistical approach, closer in philosophy
to the formalism propose here. It is important to stress that many of these
standard methods involve a number of free parameters, and we optimize their
values according to the metrics used to compare. Note that, of course, such
optimization is only possible in a synthetic tests, since the true signal in a
real problem is unknown. Thus, our results represent an upper limit to the
performance of parametric methods.

We now briefly describe the main principles and free parameters, if any, of
all the techniques that we have considered. A succinct summary is provided in
Table II.

Method | Parameters | 1D | 2D
---|---|---|---
Median | Window size ($w$) | ✓ | ✓
SGF | Window ($w$) and order ($o$) | ✓ | ✓
GF | Radius ($R_{0}$) | ✓ | ✓
Wiener | Low Frequency ($k$) | ✓ | ✓
LOWESS | Fraction window ($f$) | ✓ | -
BM3D | —– | - | ✓
FABADA | —– | ✓ | ✓
TABLE II: List of all noise reduction methods used to compared with FABADA
along their parameters and space implementations, one or two dimensions.

####  III-A1 Median filter

One of the classical non-linear digital filtering techniques, it is still
often used to remove noise from an image or signal. The main idea of the
median filter is to run through the data, replacing each point with the median
of neighboring entries. The number of neighbors used in the median is called
the ”window”, which slides, entry by entry, over the entire signal. For each
data point $z(x)$, the region used to compute the median contains, for one
dimensional data, $(w-1)/2$ neighbors on each side, whereas for two dimensions
it corresponds to a square of size $w$ centered in $z(x)$. The optimization
procedure consisted on computing the estimation $\hat{y}^{w}_{M}(x)$ for
different values of $w$, starting with the smallest value of the window
length, $w_{0}=3$, and increasing it according to
$w_{i+1}=w_{i}+2(1+w_{i}//5)$, where $//$ denotes the integer division.

####  III-A2 Savitzki-Golay filter (SGF)

As first noted by Savitzky and Golay in [1], a smoothed version of the data
may be obtained by fitting successive sub-sets of adjacent points with a low-
degree polynomial using the method of least squares. When the data are equally
spaced, the solution of the least squares (i.e., the coefficients of the
polynomials) is analytical and independent of the data to be smoothed. Thus,

$\hat{y}_{SG}^{w,o}(x)=\sum_{i=\frac{1-w}{2}}^{\frac{w-1}{2}}C_{i}^{o}(w,o)\
z(x+i)$ (16)

for $\frac{w-1}{2}\leq x\leq N_{X}-\frac{w-1}{2}$, where the two free
parameters of the method are the window length $w$ of the data region, i.e.,
the number of data points to be fitted, and the order $o$ of the polynomial.
$C_{i}^{o}$ are the $w\geq o$ Savitzky-Golay coefficients, and
$\hat{y}_{SG}^{w,o}(x)$ is the smoothed result of the filter at position $x$.
We vary the window length according to $w_{i+1}=w_{i}+2(1+w_{i}//5)$ and
$o\leq 10$ for the order of the polynomial.

Figure 1: Example spectra used to compare the performance of the different
algorithms: a Kurucz stellar atmosphere model (left), supernova remnant N132D
in the Large Magellanic Cloud (center) and interacting galaxy Arp 256 (right).

####  III-A3 LOWESS

A popular variant of the Savitzky-Golay method is the locally weighted
scatterplot smoothing (LOWESS) [2], where the regions to be fitted are not
evenly spaced and the least squares procedure takes into account weighted
values of the data, according to their distance from the point to be
evaluated. This scheme involves computing the coefficients of the fitted
polynomial each time, producing a less efficient algorithm. To compare with
FABADA we use the implementation explained in [2], which uses a linear fit and
can only be used for one dimensional data. This implementation has only one
parameter, which is the fraction $f$ of data points used to accomplish the
linear regression at each point. For the optimization we begin with the
smallest fraction possible ($f=2/N_{X}$) and then it is incremented it by a
factor of $1.5$ until the best recovery is found.

####  III-A4 Gaussian Filter (GF)

Another classical technique of noise reduction consists in filtering the high
frequency components of the data using a Gaussian filter. The fast Fourier
transform (FFT) is the most computationally efficient way to convert the data
$z(x)$ to the frequency domain $\tilde{z}(w)$. Once we have the spectrum of
the image, defined as the amplitude of the FFT of the data, we can apply a
low-pass Gaussian filter to discard the highest modes:

$GF(\tilde{z}(w),W_{0})=e^{-|\tilde{z}(w)|^{2}/2W_{0}^{2}}$ (17)

where $|\tilde{z}(w)|$ is the distance from the center (zero frequency
component) in Fourier space, and $W_{0}$ is the radius of the filter,
equivalent to $W_{0}=2\pi/R_{0}$ in configuration space. This is again a one-
parameter method, in which we select the radius $R_{0}$ of the low-pass
Gaussian filter by evaluating in twenty logarithmic steps between $R_{0}=1$
and $R_{0}=630$.

####  III-A5 Wiener Filter

The Wiener filter minimizes the mean square error taking into account that the
measurements are a random process where the statistical properties (in
particular, the spectrum) of the noise are known. In this work we have used
the implementation in the Scipy Signal Tools library [10, 11], where the
output of the filter given the signal $z(x)$ is given by

$\hat{y}_{W}(x)=\left\\{\begin{array}[]{cl}\frac{\sigma^{2}}{\sigma_{z}^{2}}m_{z}+\left(1-\frac{\sigma^{2}}{\sigma_{z}^{2}}\right)z&\sigma_{z}^{2}\geq\sigma^{2}\\\
m_{z}&\sigma_{z}^{2}<\sigma^{2}\end{array}\right.$ (18)

where $m_{z}$ and $\sigma_{z}^{2}$ are local estimates of the mean and
variance within a window of size $w$, and
$\sigma^{2}=\langle\sigma_{z}^{2}(x)\rangle$ is the the average variance
across the data. We increment $w$ starting from 3 until the best recovery is
found.

####  III-A6 BM3D

Finally, we also include the Block-Matching and 3D filtering [4] algorithm,
which arguably represents the state of the art in the research field of image
analysis [9]. A detailed account of this method, where image denoising is
implemented as two-step process, can be found in [4]. For the first step, the
noisy image is divided into equal-size square blocks. For each block, a 3D
group is formed with similar regions (block-matching), and noise is attenuated
by hard-thresholding the coefficients of a 3D transform. The filtered image is
used to estimate the energy spectrum of the signal, and the process can be
repeated a second time using a Wiener filtering instead of hard thresholding.
The final smoothed result of the image is generated as a weighted average of
the denoised blocks in their original positions. This is a non-parametric
method, and no optimization is necessary.

Figure 2: Battery of images used for the comparison procedure. From left to
right, top to bottom, the objects shown in this figure are the Bubble nebula
(NGC 7635), a galaxy cluster (Abell S1063), the Crab nebula (M1), the Eagle
nebula (M16), a spiral galaxy pair (NGC 4302 & 4298), the Ghost nebula (IC
63), Saturn and a globular star cluster (NGC 1466). Labels on the bottom left
corner are used to identify each object throughout this work.

###  III-B Data Sample

All the methods explained in the previous section are applied to a set of test
data in one and two dimensions (astrophysical spectra and monochromatic
images, respectively) with different levels of Gaussian random noise.

An important aspect in the recovery of spectra is the conservation of their
features, such as the Balmer break or emission and absorption lines, after
noise reduction. For this purpose, we consider three different spectra
(represented in Figure 1) that show these characteristics in different
degrees. The first spectrum (left) is a Kurucz stellar atmosphere model [12]
with an effective temperature $T_{eff}=11500K$, metal abundance $\log Z=0.1$
and surface gravity $\log g=5.0$, typical of an O/B type star, with a
prominent Balmer break at $\sim 400$ nm and several strong absorption lines.
The spectrum of a supernova remnant, plotted on the middle panel, is a
composite of 5 different observations [13] from the Faint Object Spectrograph
(FOS) instrument of the Hubble Space Telescope (HST). This high-resolution
spectrum ($0.9$ Å/pixel) is characterized by very prominent emission lines,
useful for inferring different physical properties of these objects. The last
spectrum (right) is taken from the TRDS Brown Atlas [14] which consist on a
pair of interacting galaxies, Arp 256, in the constellation of Cetus, and it
contains a combination of emission and absorption lines with a stellar
continuum. The Kurucz and Arp 256 spectra can be found in the Synthetic
Photometry SYNPHOT [15] Python package that simulates photometric data and
spectra, observed or otherwise. The aim of using these different spectra is to
obtain a good representation of the possible features that can appear in one-
dimensional astrophysical data and see how the different algorithms perform in
digging up spectral features out of a noisy signal.

For astronomical images, we consider eight different targets, displayed in
Figure 2, that are intended to sample the wide range of features that may be
encountered in the field, including planets, stars, diffuse nebulae, and
galaxies, either alone or in potentially blended groups. Saturn is arguably
the target that is most similar to the ordinary test images (e.g natural
landscapes, human subjects) that are often used in the context of digital
image processing. In addition, our sample includes two examples of nebulae
(Crab and Bubble) dominated by the gaseous component, two with a more
significant contribution of the stellar population (Eagle and Ghost nebulae),
and a globular cluster full of stars with different brightness. There is also
an image with a galaxy pair (NGC 4302 & NGC 4298) in which we can see two
different orientations of the galaxies, as well as a galaxy cluster with a
wide variety of morphologies and apparent sizes. All of these images have been
taken from the Hubble Space Telescope gallery produced by NASA and the Space
Telescope Science Institute (STScI). All of them have been compressed to 8-bit
values, with a maximal dynamical range of $0-255$ counts and $512\times 512$
pixels size to lighten up the computational load. For simplicity, we have also
normalized the astronomical spectra to $255$ in order to have the same
dynamical range and represent the noise in terms of this value for both
dimensions.

###  III-C Test Statistics

We apply different levels of Gaussian random noise $\eta(x)$ with constant
variance $\sigma^{2}$ to the real signal $y(x)$:

$z(x)=y(x)+\eta(x)$ (19)

where $\eta(x)=\mathcal{N}(0,\sigma)$, the element $x\in X$ denotes
independent measurements (spectrum wavelengths or image pixels), and we assume
that statistical errors are correctly characterized in the input data. Once
$z(x)$ is computed, a softened estimation $\hat{y}(x)$ of the real signal
$y(x)$ is carried out using the different algorithms explained above.

In one dimension, noise levels vary from $\sigma=5$ counts to $\sigma=95$
counts, out of the $255$ maximal value that sets the dynamical range of our
data (i.e. of the order of $\approx 2-40\%$ relative errors). We extend the
values of $\sigma$ to even higher values in two-dimensional data, specifically
to $1024$ counts ($400\%$) to illustrate the challenging situation (not so
seldom encountered in astronomy) that the signal is actually fainter than the
background noise.

We evaluate the quality of the reconstruction in terms of the Peak Signal to
Noise Ratio (PSNR) and the Structure Similarity Index Measure (SSIM) of the
estimators $\hat{y}(x)$, following common practice in the signal processing
literature.

By definition, the PSNR (usually expressed in decibels, dB) is related to the
Mean Square Error (MSE)

$\mathrm{MSE}(\hat{y})=\frac{1}{N_{X}}\sum_{x\in
X}\left(\hat{y}(x)-y(x)\right)^{2}$ (20)

as

$\mathrm{PSNR}(\hat{y})=10\cdot\log_{10}\left(\frac{255^{2}}{\mathrm{MSE}}\right),$
(21)

where $255$ is the dynamical range in our data. In principle, a more faithful
recovery of the underlying signal should yield smaller values of the MSE and
higher values of the PSNR.

Figure 3: Results obtained for the Arp 256 spectrum (see figure 1) by all the
models explained in section III-A (rows) for three different noise levels
(from left to right columns). The real signal $y(x)$, the noisy input data
$z(x)$, and the estimation $\hat{y}(x)$ are represented as red, gray and black
lines, respectively. Numbers on each panel quote the PSNR and SSIM obtained by
each method, to be compared with the noisy data (column headers).

The Structural Similarity Index Measure (SSIM) is another typical metric used
in image restoration that takes into account properties such as luminance,
contrast and structure. It is defined by the expression [16]

$\operatorname{SSIM}(\mathbf{x},\mathbf{y})=\frac{\left(2\mu_{x}\mu_{y}+C_{1}\right)\left(2\sigma_{xy}+C_{2}\right)}{\left(\mu_{x}^{2}+\mu_{y}^{2}+C_{1}\right)\left(\sigma_{x}^{2}+\sigma_{y}^{2}+C_{2}\right)}$
(22)

where $x$ and $y$ denote the two images being compared, $\mu$ and $\sigma^{2}$
are their mean and variance, and $\sigma_{xy}$ their covariance. The constants
$C_{1}$ and $C_{2}$ are two variables to stabilize the division when the
denominator approaches zero, and they are usually set to $C_{1}=(0.01\,L)^{2}$
and $C_{2}=(0.03\,L)^{2}$, where $L$ is the dynamic range of the images. The
SSIM metric can adopt values from 0 (absolute lack of correlation) to 1 (high
structural similarity).

Another metric that we consider is the CPU time used to generate the
estimation of the real data on a 2.40 GHz Intel i9-9980HK CPU along with 16Gb
DDR4 2400 MHz RAM memory. Please note that this time corresponds to the final
execution time for the given noise level in the Python implementation of the
algorithms, once the optimal parameters have been found, but it does _not_
include the time invested in the optimization, which is considerably larger.

Figure 4: Example of the estimation results for the bubble image (see Figure
2), showing the recoveries for three different values of noise
($\sigma=15,95,255$), i.e. values of PSNR ($24.63$, $8.59$, $0.02$) dB and
SSIM ($0.33$, $0.16$, $0.002$). Rows and columns represent different
algorithms and noise levels, respectively, as indicated by the labels. The
last three columns show a close-up of the bubble structure to better gauge the
recovery of fine details.

##  IV Results

We now assess the ability of our algorithm to recover the underlying signal
for the synthetic test cases described in Section III-B. In order to
facilitate the comparison with previous results reported in the literature, we
use the Peak Signal-to-Noise Ratio (PSNR) defined in (21), which is just a
measure of the Mean Squared Error (MSE), expressed in decibel (dB), as well as
the Structural Similarity Index Measure (SSIM) defined in (22). All the
results shown from the parametric methods are optimize using the PSNR metric
(similar results are obtained if SSIM was used instead).

###  IV-A Recovery Examples

Figures 3 and 4 show two examples of the results obtained by the different
algorithms explained in section III-A: the median, Savitzky-Golay filter
(SGF), Gaussian Filter (GF), Wiener filter, and locally weighted scatterplot
smoothing (LOWESS) for one-dimensional spectra, and the median, SGF, GF,
Wiener filter and block-matching (BM3D) filters for the images. Each row of
both figures provides the recovered estimations of the signal $\hat{y}(x)$ for
different noise levels, represented in each column.

####  IV-A1 Spectra Exmaple

We represent in figure 3 the recoveries obtained for three random realizations
with high, low, and extremely low signal-to-noise ratios (SNR) of the Arp 256
spectrum (see figure 1). The PSNR achieved by each method is quoted within the
corresponding panel along with the SSIM index.

At high SNR ($\sigma=10$, original PSNR$=28.13$ dB, SSIM$=0.52$), all
algorithms display not only a similar performance, but actually converge to
very similar solutions. The highest value of PSNR/SSIM is obtained by the
optimized Wiener filter (33.53 dB / 0.82), a little better than FABADA (33.15
dB / 0.82). The difference is lower than the difference between these two
recoveries and any of the others. The improvement with respect to the
originally high quality data is necessarily modest in all cases, of the order
of $\approx 3$ dB, i.e. an increase of $\sim 50\%$ in MSE. All algorithms are
able to correctly trace the presence of the most prominent emission lines, as
well as the strong Lyman-$\alpha$ absorption line near the peak at the left
end of the spectrum. Nevertheless, it is important to note that, while the
Wiener filter provides the best recovery in terms of to the overall noise
reduction, FABADA tends to preserve the true intensity of these features
slightly better than any of the other algorithms.

This trend becomes more significant as the noise increases, and it is more
difficult to discriminate significant spectral features from Gaussian random
fluctuations. In the middle panels, where $\sigma=35$, all models are able to
reproduce the overall shape of the continuum. However, they fail to recover
even the strongest absorption and emission lines, although hints of the
brightest emission lines are still present in the Wiener, LOWESS, and GF
filters. Only our prescription is able to provide a good description of these
prominent features with this level of noise in the input data, albeit weaker
absorption and emission lines are completely lost. This reflects on the values
of the metrics, where FABADA obtains the highest values ($27.93$ dB/$0.77$); a
difference of $10.68$ dB, which is more than an order of magnitude of noise
reduction with respect to the original measurements.

If we now turn to the recovery of the most noisy spectrum ($\sigma=95$, on the
right column), one can see that FABADA obtains the highest value of PSNR,
while GF and LOWESS obtain higher values of SSIM and close values of PSNR.
While the results of FABADA and GF are actually similar, all the spectral
information in LOWESS is gone. At these high levels of noise, the MSE and the
SSIM metrics are rather inadequate to assess the quality of the reconstructed
emission line spectra, because they incur in minimal penalty for failing to
reproduce a handful of peaks that are barely statistically significant. The
criteria implemented in FABADA are more conservative, and a lot of random
fluctuations are kept (hence the slightly lower SSIM) together with the most
significant remains of the actual signal. It is somewhat remarkable that
FABADA manages to recover the brightest line even at this noise level, at
variance with the other methods, while still obtaining high values of the
noise reduction metrics. We stress once again that this does not necessarily
imply a failure of the methods, but of the MSE as a goodness-of-fit indicator.
On the other hand, it does highlight the robustness of FABADA in this respect,
although the results reported here suggest that perhaps the MSE and SSIM is
not the optimal metric to gauge the quality of the recovered solution or, more
likely, that it should be complemented with another test statistic that
quantifies information loss and/or gives more weight to informative features.

####  IV-A2 Image Example

Similar trends are observed in the results obtained for the 2D data. Figure 4
shows the recovery of the Bubble nebula image for the whole set of different
algorithms, represented in each row, with three different noise levels
($\sigma=15,\ 95,\ 255$) along the columns. The PSNR and SSIM values obtained
for each are also shown in each panel, and a close-up of some structures is
provided in the last three columns to illustrate whether the smoothing methods
can reproduce their shape and edges. All models yield fairly similar
reconstructions for the high signal-to-noise case ($\sigma=15$). The best
reconstruction in terms of the MSE and SSIM is provided by the state-of-the-
art BM3D method, that improves the PSNR from $24.63$ dB to $35.30$ dB, more
than one order of magnitude of noise reduction in terms of the MSE, and from
$0.33$ to $0.89$ for the SSIM index. This is around $20\%$ more than the
recovery with FABADA ($34.22$ dB), which is the second highest value of MSE
and slightly higher than, yet similar to, the optimized Wiener filter ($34.02$
dB). In terms of the SSIM, the Wiener filter actually yields a value that is
$2\%$ higher than FABADA. The other classical filters (median, SGF, GF) obtain
comparable results when their parameters are tuned to minimize the MSE, about
$\approx 1$ dB ($\sim 25\%$) below FABADA’s solution. Regardless of the MSE
and SSIM statistics, the Wiener filter virtually misses some of the filaments
in the shape of the zoomed structure, while the other classical filters
recover some of the shapes, although in a blurrier way than FABADA and BM3D.
Visually, the latter reproduces the shape of the zoomed structure better, both
in terms of sharpness and smoothness. In particular, several edges in FABADA
seem to be a little bit more blurry compared to BM3D, together with a clearly
visible salt and pepper noise component.

Once again, the advantages of our algorithm become more evident as the noise
increases. In the middle range, with a noise level of $\sigma=95$, one may see
that FABADA’s difference in MSE with respect to BM3D decrease to $0.6$ dB and
increases to $1.11$ dB with respect to the Wiener filter, achieving higher
values of SSIM in both cases. On the other hand, the GF reaches the highest
values of MSE and SSIM, whereas the median and SGF methods display
increasingly lower values than the other methods. Most of them are able to
bring up the large-scale gradients and structures, but they fail to recover
the smallest filaments (as shown by the zoomed panels) due to the high level
of noise. The most significant problem of FABADA is still the salt and pepper
noise. The Wiener filter, despite achieving optimal metrics, returns a
composition of smoothed and noisy regions that hardly reproduces the true
underlying signal. This again give us a hint of the MSE and SSIM metric may
not be the optimal way to quantify the recovery of this type of images. On the
other hand, it is remarkable how the state-of-the-art BM3D method starts to
add some artificial edges to the recovered image.

If we push the noise even further, as it is often the case in practical
astrophysical applications, the signal itself is comparable to or even lower
than the statistical uncertainties. In our last test, the noise level is
$\sigma=255$, equal to the dynamic range of the original data. The GF obtains
again the highest values of MSE ($26.23$ dB, $0.39$ dB higher than FABADA’s
solution), almost three orders of magnitude of noise reduction, and a SSIM
value of 0.74, equal to FABADA’s. It is noteworthy to see that these two
solutions are far above the rest, by more than $\approx 2.5$ dB in PSNR and
$\approx 0.3$ in SSIM.

Upon visual inspection of the different reconstructions, including the zoomed
region, one can see that almost every algorithm, except FABADA and GF,
produces some artifacts in the images either at large or small scales. The GF,
for instance, is very efficient at filtering the high-frequency noise, but, by
construction, highlights graininess at the characteristic scale of the filter.
A similar effect is also evident in the SGF and median filters. The Wiener
filter still mixes regions with very different levels of smoothing, while BM3D
is extremely successful in eliminating the high-frequency ’grain’, albeit the
smooth areas of the real image are transformed into more staggered gradients.
This feature might be inherent to the block-matching algorithm, whose aim is
to classify similar square sections of the image into groups, thus resulting
in patches with similar gradients and/or edges. The fully automatic method
developed in this work, FABADA, seems to offer a reasonable compromise
solution. At high levels of noise, it is able to reduce the MSE and increase
the SSIM at the cost of significantly blurring the image, but the introduction
of spurious patterns, apart from salt and pepper noise, is not as severe as in
the other methods.

1D |  2D
---|---
Method | PSNR (dB) | SSIM | Method | PSNR (dB) | SSIM
all $\sigma$ | all $\sigma$ | all $\sigma$ | $\sigma>95$ | all $\sigma$ | $\sigma>95$
$\\#_{max}$ | $\delta$ | $\\#_{max}$ | $\delta$ | $\\#_{max}$ | $\delta$ | $\\#_{max}$ | $\delta$ | $\\#_{max}$ | $\delta$ | $\\#_{max}$ | $\delta$
FABADA | $18$ | 0.511 | $2$ | $0.137$ | FABADA | $20$ | $1.457$ | $20$ | 0.475 | $6$ | 0.098 | $4$ | $0.074$
LOWESS | $0$ | $1.072$ | $9$ | $0.069$ | BM3D | 63 | 1.422 | $10$ | $2.838$ | 58 | $0.102$ | $9$ | $0.197$
WIENER | 19 | $0.535$ | $13$ | $0.139$ | WIENER | $0$ | $5.382$ | $0$ | $7.804$ | $2$ | $0.181$ | $2$ | $0.237$
GF | $14$ | $0.751$ | 23 | 0.034 | GF | $25$ | $1.501$ | 22 | $0.537$ | $39$ | $0.107$ | 34 | 0.061
SGF | $6$ | $1.241$ | $10$ | $0.081$ | SGF | $0$ | $4.400$ | $0$ | $5.301$ | $3$ | $0.179$ | $3$ | $0.199$
MEDIAN | $0$ | $1.737$ | $0$ | $0.100$ | MEDIAN | $4$ | $3.131$ | $4$ | $2.277$ | $4$ | $0.181$ | $4$ | $0.179$

TABLE III: Table showing the number of tries that a specific algorithm have
recovered the best solution ($\\#_{max}$) for each metric (PSNR, SSIM). Along
these values it is shown the average euclidean distance to the maximum value,
which we defined as $\delta=<PSNR_{max}(\sigma)-PSNR(\sigma)>$. Figure 5: Peak
Signal to Noise Ratios (PSNR), Similarity Structural Indexes (SSIM), and CPU
times obtained for all the one dimensional spectra samples and noise ranges
considered in the comparison procedure. In the top set of figures is shown the
PSNR in decibels (dB), in the middle is represented the SSIM, and in the
bottom one the CPU time in seconds (s). Each figure of both groups is labeled
with the reference name given at the top of the column. The dashed yellow line
represents the PSNR and SSIM of the noisy data. Solid lines with filled
symbols refers to the non-parametric methods. In 1D data the only automatic
method is the one presented in this work, FABADA, an it is represented with a
red solid square ($\blacksquare$). Dotted lines with unfilled symbols refers
to the optimize methods. The LOWESS algorithm is represented with the blue
diamond ($\Diamond$), the Wiener filter with the orange circle
($\mathbin{\scalerel*{\circ}{j}}$), GF with the purple cross (), SGF with the
brown hexagon ($\hexagon$), and the median filter with the green triangle
($\bigtriangleup$). Figure 6: Peak Signal to Noise Ratios (PSNR), Similarity
Structural Indexes (SSIM), and CPU times obtained for all the two dimensional
images samples and noise ranges considered in the comparison procedure. In the
top set of figures is shown the PSNR in decibels (dB), in the middle is
represented the SSIM, and in the bottom one the CPU time in seconds (s). Each
panel is labeled with the reference image name given. The dashed yellow line
represents the PSNR and SSIM of the noisy data. Solid lines with filled
symbols refers to the non-parametric methods. In 2D data there are two
automatic methods, FABADA, which it is represented with a red solid square
($\blacksquare$) and the block matching 3D filtering BM3D, with blue stars
(★). Dotted lines with unfilled symbols refers to the optimize methods. The
Wiener filter is represented with the orange circle
($\mathbin{\scalerel*{\circ}{j}}$), GF with the purple cross (), SGF with the
brown hexagon ($\hexagon$), and the median filter with the green triangle
($\bigtriangleup$).

###  IV-B Entire database results

We have computed the PSNR, SSIM and CPU time for all data samples, noise
levels, and methods as a function of $\sigma$, averaged over ten different
random realizations. The number of cases that a specific algorithm has
recovered the best solution, as well as the average difference with respect to
the optimal choice, are quoted in Table III.

####  IV-B1 One dimension - Spectra

In one dimension (Figure 5), FABADA behaves very similarly to the optimized
Wiener filter, particularly in terms of PSNR (top panel). In general, the best
results regarding this metric are achieved with the latter method; 19 out of
the 57 ($33\%$) test cases, followed by FABADA, who obtained the best PSNR for
18 ($32\%$) estimations, being especially successful in the low signal-to-
noise regime. Then, the optimized GF obtained the best reconstruction in 14
($25\%$) cases, and the remaining 6 ($10\%$) correspond to the optimized SGF.
Neither LOWESS nor the median filter obtained in any case the best values for
PSNR. Just taking into account these results, one can see how FABADA preforms
as well as the best possible solutions of the standard methods typically used
in astronomy. The average difference with respect to the highest PSNR achieved
by any algorithm is only $0.511$ dB. This supports the idea that FABADA
automatically achieves the limit of the standard methods, when their
parameters are tuned to the (unknown) optimal values.

In terms of the SSIM metric (middle panels), our algorithm obtains slightly
worse results, being the best option only in $2$ ($4\%$) test cases, with an
average distance with respect to the highest SSIM of $0.137$. The optimized GF
filter obtained $23$ ($40\%$) of the highest values and an average difference
of $0.034$. This result is due in part to the $SN132D$ spectra, where most of
the information is concentrated in the emission lines. At high values of the
noise level, almost all standard methods converge towards a horizontal line,
where all the information is lost (see e.g. LOWESS in Figure 3). However, the
SSIM (and, to some extent, PSNR) do not strongly penalize the misfitting of
the emission lines when their statistical significance becomes low, which
hints the limitations of these metrics to quantify the recovery of relevant
information. Discarding the SN132D spectra, FABADA is seldom the optimal
choice, but the average distance decreases to $0.048$ in terms of SSIM,
compared to $0.022$ for the GF and $0.056-0.085$ for the other algorithms.

As regards to the CPU time, it is easy to see that FABADA has a strong
dependence on the level of noise, at variance with most classical methods, due
to the increasing number of iterations required to fulfill our $\chi^{2}$
stopping condition. This behavior is also seen in the median, due to a similar
increase of the optimal window size, whether other methods are less sensitive
to the noise level. FABADA is, in general, significantly slower than the other
algorthims (except LOWESS, at high signal-to-noise ratios), although it must
be borne in mind that we have not taken into account the time consumed by the
optimization process of the standard algorithms, only the execution time once
the optimal parameters have been found.

####  IV-B2 Two dimensions - Images

For the two-dimensional images (figure 6) BM3D reaches the highest PSNR for
$63$ out of our $112$ ($56\%$) test cases with different target and noise
levels, followed by the optimized GF with $25$ ($22\%$), FABADA with $20$
($18\%$), and $4$ ($4\%$) that correspond to the optimized median filter.
Neither SGF nor Wiener filter ever obtained the highest values. Nevertheless,
BM3D is on average $1.42$ dB below the optimal choice, comparable to the
$1.46$ dB of FABADA and $1.50$ dB of GF, due to the significantly different
trends observed at different noise levels.

In general, BM3D stands over the other methods, including FABADA, at high SNR
($\sigma\lesssim 95$ dB), in particular for the Saturn image. Its
collaborative filter is particularly well suited for periodic data, or images
with repetitive patterns, which are virtually absent in other test cases. The
stars image would be a paradigmatic example, and the difference in this test
is insignificant.

On the other hand, FABADA and GF predominate in the low-SNR regime
($\sigma\gtrsim 95$ dB), where BM3D only achieves the best reconstruction in
$10$ out of the $56$ ($18\%$) test cases. The remaining $46$ are achieved by
GF with $22$ ($39\%$), FABADA with $20$ ($36\%$), and the median filter with
$4$ ($7\%$). Furthermore, BM3D is on average $2.84$ dB below the highest
value, while FABADA and GF are only $0.48$ and $0.54$ dB lower, respectively.
Essentially, FABADA and the optimized GF achieve remarkably similar results,
and they preform better than BM3D at low signal-to-noise ratios.

This behavior is also seen in the middle panels, where the SSIM is
represented. From the $58$ ($52\%$) of the highest values achieved by BM3D,
only $9$ ($20\%$) are above $\sigma>95$ counts, and the difference with
respect to the highest value of SSIM increases with noise level from $0.1$ to
$0.2$. Although in this metric FABADA only achieves the best recovery in $6$
($5\%$) cases, $4$ ($7\%$) at high noise levels, it is only $0.1$ apart, on
average, from the best solution and $0.07$ at the high noise regime. The
optimized GF produces better estimations than FABADA in $39$ ($35\%$) cases,
$34$ ($60\%$) in the low SNR regime. Despite this difference, GF is still
$0.1$ below the best solution ($0.06$ at low SNR), again very close to the
results obtained by our algorithm.

This change on the trend of the best solution is at the cost of CPU time.
While FABADA seems to converge to its solution in an efficient way in the high
SNR regime (always below the second), for high noise levels (low SNR) the time
rises up to $30$ s, as shown in the bottom panels. Once again, for the
optimized models we only take into account the execution time with the fine-
tuned parameters already known. The results are similar to the one-dimensional
case, but the time scale increases with the dimensionality of the problem. The
optimized Wiener filter is still the fastest algorithm, obtaining an average
time of $0.02$ seconds, while median filter seems to be the slowest, with
$9.34$ seconds. FABADA also features a high execution time ($8.08$ seconds on
average), whereas BM3D, GF, and SGF yield $6.7$, $0.87$, and $0.05$ seconds,
respectively. However, the performance of FABADA and the median filter vary
significantly as a function of SNR. Combining these results with the above
metrics for the quality of the reconstruction, we argue that our method
provides a competitive alternative both at high SNR, where it achieves
slightly less reliable results than the state-of-the-art algorithm BM3D at a
fraction of the computational cost, as well as in the low SNR regime, where it
provides a more faithful reconstruction after a significantly larger execution
time.

##  V Conclusions

In this work we present the theory and implementation of a novel automatic
algorithm for noise reduction: the Fully Adaptive Bayesian Algorithm for Data
Analysis (FABADA). Our method iteratively evaluates progressively smoother
models of the underlying signal, and then combines them according to their
Bayesian evidence and $\chi^{2}$ statistic. The source code is publicly
available at https://github.com/PabloMSanAla/fabada.

We compare FABADA with other methods that are representative of the current
state of the art in image analysis and digital signal processing. For this
comparison we used the most typical metrics, the Peak Signal to Noise Ratio
(PSNR), which is a measure of the Mean Square Error (MSE), the Structural
Similarity Index (SSIM) and the CPU time. One important advantage of our
method, shared by BM3D, over classical algorithms is the absence of free
parameters to be tuned by the user. Our results suggest that FABADA and BM3D
achieve values of PSNR or SSIM comparable to or better than the best possible
solution attainable by the classical methods.

Beyond the precise values of the global quantitative metrics, both FABADA and
BM3D are quite successful in adapting to the structures present in the input
data. Perhaps the most significant difference between them is that FABADA’s
priors assume that the signal is smooth, whereas BM3D uses block-matching to
look for repetitive patterns. This might be relevant when one must recover the
height and shape of the spectral features in 1D or the gradients and
boundaries in 2D. We argue that FABADA appears to offer a trade-off between
noise reduction, increasing the metric values significantly, in a way that is
statistically compatible with the data, keeping significant features without
introducing considerable artifacts.

Regarding execution time, non-parametric methods are more expensive than the
classical alternatives _once the optimal values of their parameters are known_
, something that is of course impossible in practice. FABADA is faster than
BM3D at high SNR, although it usually yields a poorer reconstruction, and the
trend reverses for low SNR.

## References

  * [1] A. Savitzky and M. J. E. Golay, “Smoothing and differentiation of data by simplified least squares procedures.” _Analytical Chemistry_ , vol. 36, no. 8, pp. 1627–1639, 1964. [Online]. Available: https://doi.org/10.1021/ac60214a047
  * [2] W. S. Cleveland, “Robust locally weighted regression and smoothing scatterplots,” _Journal of the American Statistical Association_ , vol. 74, no. 368, pp. 829–836, 1979. [Online]. Available: https://www.tandfonline.com/doi/abs/10.1080/01621459.1979.10481038
  * [3] B. Goyal, A. Dogra, S. Agrawal, B. Sohi, and A. Sharma, “Image denoising review: From classical to state-of-the-art approaches,” _Information Fusion_ , vol. 55, pp. 220–244, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1566253519301861
  * [4] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Image denoising by sparse 3-d transform-domain collaborative filtering,” _IEEE Transactions on Image Processing_ , vol. 16, no. 8, pp. 2080–2095, 2007.
  * [5] S. Gu, L. Zhang, W. Zuo, and X. Feng, “Weighted nuclear norm minimization with application to image denoising,” in _2014 IEEE Conference on Computer Vision and Pattern Recognition_ , 2014, pp. 2862–2869.
  * [6] A. E. Ilesanmi and T. O. Ilesanmi, “Methods for image denoising using convolutional neural network: a review,” _Complex & Intelligent Systems_, 2021. [Online]. Available: https://doi.org/10.1007/s40747-021-00428-4
  * [7] K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang, “Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising,” _IEEE Transactions on Image Processing_ , vol. 26, no. 7, pp. 3142–3155, 2017.
  * [8] V. Katkovnik, V. Katkovnik, K. Egiazarian, and J. Astola, _Local Approximation Techniques in Signal and Image Processing_ , ser. SPIE Press monograph. SPIE Press, 2006. [Online]. Available: https://books.google.to/books?id=YChTAAAAMAAJ
  * [9] M. El Helou and S. Susstrunk, “Blind universal bayesian image denoising with gaussian noise level learning,” _IEEE Transactions on Image Processing_ , vol. 29, p. 4885–4897, 2020. [Online]. Available: http://dx.doi.org/10.1109/TIP.2020.2976814
  * [10] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” _Nature Methods_ , vol. 17, pp. 261–272, 2020.
  * [11] J. S. Lim, _Two-dimensional signal and image processing_. Englewood Cliffs, N.J. ; Prentice Hall, 1990.
  * [12] F. Castelli and R. L. Kurucz, “New Grids of ATLAS9 Model Atmospheres,” in _Modelling of Stellar Atmospheres_ , N. Piskunov, W. W. Weiss, and D. F. Gray, Eds., vol. 210, Jan. 2003, p. A20.
  * [13] W. P. Blair, J. A. Morse, J. C. Raymond, R. P. Kirshner, J. P. Hughes, M. A. Dopita, R. S. Sutherland, K. S. Long, and P. F. Winkler, “Hubble Space Telescope Observations of Oxygen-rich Supernova Remnants in the Magellanic Clouds. II. Elemental Abundances in N132D and 1E 0102.2-7219,” _Astrophysical Journal_ , vol. 537, no. 2, pp. 667–689, Jul. 2000.
  * [14] M. J. I. Brown, J. Moustakas, J. D. T. Smith, E. da Cunha, T. H. Jarrett, M. Imanishi, L. Armus, B. R. Brandl, and J. E. G. Peek, “An Atlas of Galaxy Spectral Energy Distributions from the Ultraviolet to the Mid-infrared,” _Astrophysical Journal, Supplement_ , vol. 212, no. 2, p. 18, Jun. 2014.
  * [15] P. L. Lim, R. I. Diaz, and V. Laidler, “Pysynphot user’s guide (baltimore, md: Stsci),” 2015. [Online]. Available: http://synphot.readthedocs.io/en/latest/
  * [16] Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” _IEEE Transactions on Image Processing_ , vol. 13, no. 4, pp. 600–612, 2004.

##  VI Acknowledgments

Pablo M. Sánchez-Alarcón acknowledges support from the State Research Agency
(AEI-MCINN) of the Spanish Ministry of Science and Innovation under the grant
”The structure and evolution of galaxies and their central regions” with
reference PID2019-105602GB-I00/10.13039/501100011033, and from IAC project
P/300724, financed by the Ministry of Science and Innovation, through the
State Budget and by the Canary Islands Department of Economy, Knowledge and
Employment, through the Regional Budget of the Autonomous Community. Yago
Ascasibar acknowledges financial support from grant ”Starbursts throughout the
evolution of the Universe” (PID2019-107408GB-C42) from the AEI-MICINN, Spain.

|  Pablo Manuel Sánchez Alarcón was born in Madrid, Spain, in 1996. He
received the B.S. degree in physics from the Universidad Autónoma de Madrid,
Spain in 2019, the M.S. degree in theoretical physics from Universidad
Autónoma de Madrid and the Instituto de Física Teórica, Madrid, Spain in 2020.
He is currently pursuing the Ph.D. degree in astronomy, specializing in galaxy
evolution, with the Universidad de La Laguna, Tenerife, Spain. He is also a
researcher at the Instituto de Astrofísica de Canarias.
---|---
|  Yago Ascasibar Sequeiros was born in Madrid, Spain, in 1975. He received
the Ph.D. degree in Astrophysics and Cosmology from the Universidad Autónoma
de Madrid, Madrid, Spain, in 2003. He held postdoctoral research positions at
the University of Oxford, UK, in 2003, the Harvard-Smithsonian Center for
Astrophysics, USA, in 2004, the Leibnitz Institut für Astrophysik Potsdam,
Germany, in 2005, and the Universidad Autónoma de Madrid, Spain, in 2008. He
is currently Associate Professor at that institution, since 2019. His main
research interests are galaxy formation and evolution, dark matter, statistics
and data analysis. Dr. Ascasibar is member of the Spanish Astronomical
Society, European Astronomical Society, and International Astronomical Union.
---|---


# Defining a credible interval is not always possible with “point-null”
priors: A lesser-known correlate of the Jeffreys-Lindley paradox

Harlan Campbell and Paul Gustafson
Department of Statistics, University of British Columbia

###### Abstract

In many common situations, a Bayesian credible interval will be, given the
same data, very similar to a frequentist confidence interval, and researchers
will interpret these intervals in a similar fashion. However, no predictable
similarity exists when credible intervals are based on model-averaged
posteriors whenever one of the two nested models under consideration is a so
called “point-null”. Not only can this model-averaged credible interval be
quite different than the frequentist confidence interval, in some cases it may
be undefined. This is a lesser-known correlate of the Jeffreys-Lindley paradox
and is of particular interest given the popularity of the Bayes factor for
testing point-null hypotheses.

##  1 Introduction

Recently, several Bayesian tests using Bayes factors have been proposed as
alternatives to frequentist hypothesis testing; see Heck et al. (2022) for a
recent review. When using the Bayes factor (or the posterior model odds) for
testing, it is often recommended that researchers also report parameter
estimates and their credible intervals (e.g., Keysers et al. (2020)). Indeed,
following a controversial debate about the strict binary nature of statistical
tests, many now call for an additional focus on parameter estimation with
appropriate uncertainty estimation; see Wasserstein & Lazar (2016).

Campbell & Gustafson (2022) consider how Bayesian testing and estimation can
be done in a complimentary manner and conclude that if one reports a Bayes
factor comparing two models, then one should also report a model-averaged
credible interval (i.e., one based on the posterior averaged over the two
models under consideration). Researchers who follow this recommendation can
obtain credible intervals congruent with their Bayes factor, thereby obtaining
suitable uncertainty estimation.

In many familiar situations, a posterior credible interval will be, given the
same data, very similar to a frequentist confidence interval and researchers
will interpret these intervals in a similar fashion; see Albers et al. (2018).
However, when comparing two models, one of which involves a so-called “point-
null”, it is less clear whether or not such similarity can be assumed.

Previous work has examined the properties of Bayesian credible intervals and
how they relate to frequentist confidence intervals under various prior
specifications (e.g., Casella & Berger (1987), Datta & Ghosh (1995), Greenland
& Poole (2013), Held (2020) ). In this paper, on the basis of a few simple
examples, we will examine properties specific to model-averaged credible
intervals. We will show that, when one of the two models under consideration
is a point-null model, not only can the model-averaged credible interval be
quite different than the confidence interval, oftentimes, for a desired
probability level, it may be undefined. This is perhaps an unexpected
correlate of the Jeffreys-Lindley paradox, the most well known example of the
rift between frequentist and Bayesian statistical philosophies; see
Wagenmakers & Ly (2021). The limitations/particularities of working with
point-null models are of particular interest given the recent popularity of
the Bayes factor for testing point-null hypotheses.

We begin in Section 2 by re-visiting an example of two Normal models
considered previously by Wagenmakers & Ly (2021) in their discussion of the
Jeffreys-Lindley paradox. In Section 3, we extend this example to consider the
consequences of specifying a point-null model. We conclude in Section 4 with
thoughts on the consequences, with respect to parameter estimation, of
specifying point-null models.

##  2 A mixture of two Normals

Let $\theta$ be the parameter of interest for which there are two a priori
probable models: $M_{0}$ and $M_{1}$, defined by two different priors
$\pi_{0}(\theta)$ and $\pi_{1}(\theta)$. The posterior density which
appropriately acknowledges the uncertainty with regards to which of the two
models is correct is the mixture density:

$\displaystyle\pi(\theta|data)=$
$\displaystyle\textrm{Pr}(M_{0}|data)\pi_{0}(\theta|data)+\textrm{Pr}(M_{1}|data)\pi_{1}(\theta|data),$
(1)

where the model-specific posteriors, $\pi_{0}(\theta|data)$ and
$\pi_{1}(\theta|data)$, are weighted by their posterior model probabilities,
$\textrm{Pr}(M_{0}|data)$ and $\textrm{Pr}(M_{1}|data)$; see, for instance,
Campbell & Gustafson (2022). Note that this “mixture” posterior is obtained as
a result of specifying the “mixture” prior:

$\pi(\theta)=\textrm{Pr}(M_{0})\pi_{0}(\theta)+\textrm{Pr}(M_{1})\pi_{1}(\theta),$
(2)

where $\textrm{Pr}(M_{0})$ and $\textrm{Pr}(M_{1})$ are the a priori model
probabilities.

As an example, consider two a priori equally probable Normal models,
$M_{0}:\theta\sim N(0,g_{0})$ and $M_{1}:\theta\sim N(0,g_{1})$, such that
$\textrm{Pr}(M_{0})=\textrm{Pr}(M_{1})=0.5$. The prior density functions for
the two models are defined as:

$\displaystyle\pi_{0}(\theta)$ $\displaystyle=f_{Normal}(\theta,0,g_{0}),$ (3)

and

$\displaystyle\pi_{1}(\theta)$ $\displaystyle=f_{Normal}(\theta,0,g_{1}),$

where $f_{Normal}(x,\mu,\sigma^{2})$ is the Normal probability density
function evaluated at $x$, with mean parameter $\mu$ and variance parameter
$\sigma^{2}$. Let $y_{i}$ be the $i$-th data-point, for $i=1,...,n$; let
$\bar{y}=\sum_{i=1}^{n}{y_{i}}/n$ be the sample mean; and suppose these data
are normally distributed with known unit variance such that:

$\displaystyle\textrm{Pr}(data|\theta)$
$\displaystyle=\prod_{i=1}^{n}f_{Normal}(y_{i},\theta,1).$

Then the Bayes factor is:

$\displaystyle\textrm{BF}_{01}=\sqrt{\frac{1+ng_{1}}{1+ng_{0}}}\times\textrm{exp}\Big{(}\frac{(g_{0}-g_{1})nz^{2}}{2(1+ng_{0})(1+ng_{1})}\Big{)},$

where $z=\sqrt{n}\bar{y}$. The posterior model probabilities can be calculated
from the Bayes factor as:

$\displaystyle\textrm{Pr}(M_{0}|data)=\frac{\textrm{Pr}(M_{0})}{\textrm{Pr}(M_{1})/\textrm{BF}_{01}+\textrm{Pr}(M_{0})}\quad\textrm{and}\quad\textrm{Pr}(M_{1}|data)=1-\textrm{Pr}(M_{0}|data).$
(4)

Finally, the model specific posteriors are defined as:

$\displaystyle\pi_{j}(\theta|data)=$ $\displaystyle
f_{Normal}\Big{(}\theta,\frac{zg_{j}}{\sqrt{n}(\frac{1}{n}+g_{j})},{\frac{g_{j}}{1+g_{j}n}}\Big{)},$

for $j=0,1$.

Having established all the components of equation (1), let us now consider how
to define a credible interval based on the model-averaged posterior. An upper
one-sided $(1-\alpha)\%$ credible interval is defined as:

$\displaystyle\textrm{one-sided }(1-\alpha)\%\textrm{CrI}$
$\displaystyle=[\theta^{*},\infty),$

where $\theta^{*}$ satisfies the following equality:

$\displaystyle\textrm{Pr}(\theta<\theta^{*}|data)=\alpha.$ (5)

Let us define an equal-tailed two-sided $(1-\alpha)\%$ credible interval from
a combination of two upper one-sided intervals as:

$\displaystyle\textrm{two-sided }(1-\alpha)\%\textrm{CrI}$
$\displaystyle=[\theta^{l*},\theta^{u*}),$

where $\theta^{l*}$ and $\theta^{u*}$ satisfy:
$\textrm{Pr}(\theta<\theta^{l*}|data)=\alpha/2$ and
$\textrm{Pr}(\theta<\theta^{u*}|data)=1-\alpha/2$. Note that, in our example
of two Normal models, these posterior values are calculated as:

$\displaystyle\textrm{Pr}(\theta<\theta^{*}|data)$
$\displaystyle=\int_{-\infty}^{\theta^{*}}\pi(\theta|data)d\theta=\frac{\int_{-\infty}^{\theta^{*}}\Big{(}f_{Norm}((z-\theta\sqrt{n}),0,1)\times\pi(\theta)\Big{)}d\theta}{\int_{-\infty}^{\infty}\Big{(}f_{Norm}((z-\theta\sqrt{n}),0,1)\times\pi(\theta)\Big{)}d\theta},$

where $\pi(\theta)$ is defined as in equation (2),  and the integral in the
denominator ensures that the posterior density integrates to one.

Now suppose $g_{0}=0.02$, $g_{1}=1$ and that we observe data for which
$\bar{y}=1.645/\sqrt{n}$ which corresponds to a $p$-value of $p=0.05$ when
using these data to test against the null hypothesis
$\textrm{H}_{0}:\theta<0$. See Figure 1 which plots priors and posteriors for
this scenario with $n=10$. The lower bound of an upper one-sided $(1-A)$%
confidence interval (CI) will be equal to
$\textrm{CI}_{A}=\bar{y}-\frac{Q_{Norm}(1-A)}{\sqrt{n}}$, where $Q_{Norm}()$
is the Normal quantile function. For instance, for the observed data with
$\bar{y}=1.645/\sqrt{n}$, we have $\textrm{CI}_{0.10}=(1.645-1.282)/\sqrt{n}$,
such that an upper one-sided 90% CI will be $=[0.363/\sqrt{n},\infty)$. An
upper one-sided 95% CI for these data will be $[0,\infty)$, since
$\textrm{CI}_{0.05}=0$. How do these frequentist intervals compare to model-
averaged Bayesian credible intervals? While most literature describing the
asymptotic agreement of Bayes and frequentist inferences considers the regime
of a fixed true parameter value as $n$ increases, for our purposes it is
useful to consider the regime of a fixed $p$-value for a particular point null
hypothesis. Consider two observations.

Figure 1: For the “mixture of two normals” example ($g_{0}=0.02$ and
$g_{1}=1$), panels A, B, and C, plot the $M_{0}$ prior, the $M_{1}$ prior, and
the mixture-prior, respectively. For data with $\bar{y}=0.520$ and $n=10$,
panels D, E, and F, plot the $M_{0}$ posterior, the $M_{1}$ posterior, and the
model-averaged posterior, respectively.

First, setting $\theta^{*}=0$ in equation (5), we see that as $n$ increases
(and $p=0.05$ remains fixed), the corresponding value of $\alpha$ approaches
$p=0.05$: For $n=10$, we obtain $\alpha=0.160$, whereas for $n=10000$, we
obtain $\alpha=0.050$; see how the
$\textrm{Pr}(\theta<\textrm{CI}_{0.05}|data)$ curve approaches 0.05 as $n$
increases in Figure 2. Second, setting $\theta^{*}=\textrm{CI}_{A}$ in
equation (5), we see that as $n$ increases (and $p=0.05$ remains fixed), the
corresponding value of $\alpha$ approaches $A$. In Figure 2, we plot values of
$\alpha$ corresponding to $A=0.05,0.10,0.20$, and $0.30$. One can clearly see
that each $\textrm{Pr}(\theta<\textrm{CI}_{A}|data)$ curve tends
asymptotically towards $A$. One can verify this asymptotic behaviour by re-
expressing posterior expectations arising from the specified prior as
posterior expectations arising under an improper uniform prior. For
completeness, we give the necessary details in the Appendix.

Figure 2: Let $\textrm{CI}_{A}$ be the lower bound of a frequentist upper one-
sided (1-A)% confidence interval. We consider
$\textrm{Pr}(\theta<\textrm{CI}_{A}|data)=\alpha$ and data corresponding to
($n$, $p$), where $n$ is the sample size and $p$ is the frequentist $p$-value
obtained when testing the data against the null hypothesis
$\textrm{H}_{0}:\theta<0$. For the normal mixture example with $g_{0}=0.02$
and $g_{1}=1$, and $p=0.05$, we have $\textrm{CI}_{0.05}=0$ and see that, as
$n$ increases, $\alpha$ approaches $A$ for $A=0.05,0.10,0.20$ and $0.3$.

Based on the asymptotic behaviour of the posterior in this example, one might
reasonably conclude that, with a sufficiently large sample size, the model-
averaged credible interval will approximate the frequentist’s confidence
interval for any $A$ probability level. However, Wagenmakers & Ly (2021) argue
that, in this scenario, “the Jeffreys-Lindley paradox still applies”
indicating that there is indeed a conflict between Bayesian and frequentist
interpretations of the data.

Wagenmakers & Ly (2021) explain their reasoning as follows. From equation (4),
we calculate
$\textrm{lim}_{n\rightarrow\infty}\textrm{Pr}(M_{1}|data)=(1+\sqrt{g_{1}/g_{0}})^{-1}=(1+1/\sqrt{0.02})^{-1}=0.124$
and $\textrm{lim}_{n\rightarrow\infty}\textrm{Pr}(M_{0}|data)=0.876$.
Therefore, with sufficiently large $n$, we have that
$\textrm{Pr}(M_{1}|data)<\textrm{Pr}(M_{0}|data)$ regardless of the data
(i.e., regardless of the fixed value of $z=\sqrt{n}\bar{y}$); see Figure 3.

In this scenario, model selection (i.e., evaluating the relative values of
$\textrm{Pr}(M_{0}|data)$ and $\textrm{Pr}(M_{1}|data)$) is not addressing the
same question as estimation (i.e., evaluating $\textrm{Pr}(\theta|data)$ to
determine which values of $\theta$ are a posteriori most likely). The
posterior density of $\theta$ describes one’s belief in the probability of
different possible values of $\theta$, whereas the posterior model
probabilities describe the probability of different data generating processes
(DGP) (including the generation of $\theta$). As such, while it is true that
the Jeffreys-Lindley paradox still applies with regards to model selection
(i.e., with a sufficiently large sample size and fixed $z$, the Bayesian will
inevitably select $M_{0}$), the paradox does not apply when it comes to
parameter estimation (i.e., with a sufficiently large sample size and fixed
$z$, the Bayesian will inevitably agree with the frequentist when it comes to
estimating $\theta$, with their credible interval approximately equal to the
frequentist’s confidence interval). One way to think about this is to consider
the diminishing influence of the prior as the sample size increases and to
recall that the confidence interval and the credible interval will agree
exactly if one specifies the flat (albeit improper) reference prior,
$\pi(\theta)\propto 1$; see the worked examples in Held (2020).

In order for the Jeffreys-Lindley paradox to apply to parameter estimation, a
point-mass in the prior is required. We consider this situation in the next
Section.

Figure 3: For the normal mixture model example with $g_{0}=0.02$ and
$g_{1}=1$, the $\textrm{Pr}(M_{1}|data)$ (blue curve) increases towards 0.876
with increasing $n$, while the value of $\textrm{Pr}(\theta<0|data)$ (grey
line) approaches 0.05 (dashed black line).

##  3 Parameter estimation with a point null

Consider the same scenario as above but with the null model, $M_{0}$, defined
as a so-called “point-null” such that the prior density function under $M_{0}$
is:

$\displaystyle\pi_{0}(\theta)$ $\displaystyle=\delta_{0}(\theta),$ (6)

where $\delta_{0}()$ is the Dirac delta function at 0 which can be informally
thought of as setting $g_{0}=0$ in equation (3), or alternatively thought of
as a probability density function which is zero everywhere except at 0, where
it is infinite. Note that these are merely informal, intuitive
interpretations.

Figure 4: For the “point-null” example, panels A, B, and C, plot the $M_{0}$
prior, the $M_{1}$ prior, and the mixture-prior, respectively. For data with
$\bar{y}=0.520$ and $n=10$, panels D, E, and F, plot the $M_{0}$ posterior,
the $M_{1}$ posterior, and the model-averaged posterior, respectively.

We now have that $\textrm{Pr}(\theta=0|data)=\textrm{Pr}(M_{0}|data)$, or
equivalently, $\textrm{Pr}(\theta\neq 0|data)=\textrm{Pr}(M_{1}|data)$. As
such, model selection (selecting between $M_{0}$ and $M_{1}$) and null
hypothesis testing (selecting between $\textrm{H}_{0}:\theta=0$ and
$\textrm{H}_{1}:\theta\neq 0$) are equivalent in this scenario.

With the “point-null” prior for $M_{0}$ as defined in (6), and with $g_{1}=1$,
as defined previously, the “mixture” prior, $\pi(\theta)$, is recognizable as
a “spike-and-slab” prior (see van den Bergh et al. (2021)) and the Bayes
factor is equal to:

$\displaystyle\textrm{BF}_{01}=\sqrt{{1+n}}\times\textrm{exp}\Big{(}\frac{-nz^{2}}{2(1+n)}\Big{)},$

The posterior density is nonatomic with a spike (i.e., a discontinuity with
infinite density) at 0:

$\pi(\theta|data)=\textrm{Pr}(M_{0}|data)\delta_{0}(\theta)+\textrm{Pr}(M_{1}|data)f_{Normal}\Big{(}\theta,\frac{z}{\sqrt{n}(\frac{1}{n}+1)},{\frac{1}{1+n}}\Big{)},$

where the posterior model probabilities, $\textrm{Pr}(M_{0}|data)$ and
$\textrm{Pr}(M_{1}|data)$, can be calculated from the Bayes factor as in
equation (4).

Returning to our hypothetical data with $z=1.645$, we see that for
$\theta^{*}=0$, as $n$ increases, $\alpha$ (such that
$\textrm{Pr}(\theta<\theta^{*}|data)=\alpha$) does not approach $p=0.05$ and
instead approaches 0: For $n=10$, we obtain $\alpha=0.03$, and for $n=1000$,
we obtain $\alpha=0.005$; see trajectory of the grey curve in Figure 5.
Whatsmore, as $n$ increases and $\bar{y}=1.645/\sqrt{n}$ remains fixed, the
posterior probability on the “spike” at $0$ increases towards infinity such
that: $\textrm{lim}_{n\rightarrow\infty}\textrm{Pr}(M_{0}|data)=1$; as
famously emphasized by Lindley (1957) and originally demonstrated by Jeffreys
(1935).

Perhaps even more puzzling is that, for fixed $\alpha=0.05$, there is simply
no corresponding value of $\theta^{*}$ (such that
$\alpha=\textrm{Pr}(\theta<\theta^{*}|data)$) for any $n>2$. For $n=2$ we can
define $\theta^{*}=-0.0163$, such that
$\textrm{Pr}(\theta<-0.0163|data)=0.05$. However, for $n=3$, a precise value
of $\theta^{*}$ cannot be defined since, due to the discontinuity in the
posterior, we have: $\textrm{Pr}(\theta<0|data)=0.045<\alpha$, and
$\textrm{Pr}(\theta\leq 0|data)=0.465>\alpha$. For $n=10$ the gap is even
wider: $\textrm{Pr}(\theta<0|data)=0.030<\alpha$ and $\textrm{Pr}(\theta\leq
0|data)=0.522>\alpha$. Figure 5 plots these numbers for increasing values of
$n$. As a consequence, it is no longer the case that, with a sufficiently
large sample size, a Bayesian’s credible interval will approximate a
frequentist’s confidence interval. In fact, for certain values of $\alpha$ and
$n$, calculating a credible interval is not even possible.

Figure 5: For the hypothetical data with $z=1.645$, as $n$ increases along the
horizontal axis, values of $\alpha$ such that
$\textrm{Pr}(\theta<0|data)=\alpha$ (grey line) and $\textrm{Pr}(\theta\leq
0|data)=\alpha$ (red line) are plotted on the vertical axis. Figure 6: For the
hypothetical data with $z=1.645$ and $n=10$, the plotted line corresponds to
the cumulative distribution function of the posterior (i.e.,
$\textrm{Pr}(\theta<\theta^{*}|data)$) for increasing values of $\theta^{*}$.

In general, determining a specific value of $\theta^{*}$ for a given value of
$\alpha$ (such that $\alpha=\textrm{Pr}(\theta<\theta^{*}|data)$) is only
possible for values of $\alpha$ outside of the “incredibility interval”:

$\Big{[}\Big{(}\textrm{Pr}({\theta}<0|data,M_{1})\textrm{Pr}(M_{1}|data)\Big{)},\Big{(}\textrm{Pr}({\theta}<0|data,M_{1})\textrm{Pr}(M_{1}|data)+\textrm{Pr}(M_{0}|data)\Big{)}\Big{]}.$

The bounds of the “incredibility interval” are the limits of the “jump” in the
cumulative distribution function of the posterior, i.e., the values ranging
between $\textrm{Pr}(\theta<0|\textrm{data})$ and $\textrm{Pr}(\theta\leq
0|\textrm{data})$. In Figure 6, we plot the cumulative distribution function
of the posterior for hypothetical data with $z=1.645$ and $n=10$. In this
situation, the “incredibility interval” equals
$[\textrm{Pr}(\theta<0|data),\textrm{Pr}(\theta\leq 0|data)]=[0.03,0.522]$. In
Figure 5, the lower grey curve corresponds to the lower bound of the
incredibility interval and the upper red curve corresponds to the upper bound.
Notably, since $\textrm{lim}_{n\rightarrow\infty}\textrm{Pr}(M_{0}|data)=1$
and $\textrm{lim}_{n\rightarrow\infty}\textrm{Pr}(M_{1}|data)=0$, the width of
the incredibility interval increases as $n$ increases. As a result,
determining a precisely $\alpha$-level value of $\theta^{*}$ such that
$\alpha=\textrm{Pr}(\theta<\theta^{*}|data)$, becomes increasingly impossible
as $n$ grows large. This is true regardless of the data; see Figure 7 for
values of the lower bound obtained with data where $\bar{y}=2.575/\sqrt{n}$
(data for which one obtains a $p$-value of $p=0.005$ when testing against
$\textrm{H}_{0}:\theta<\theta_{0}$).

Figure 7: With data where $\bar{y}=2.575/\sqrt{n}$, as $n$ increases, the
lower bound of the incredibility interval (the solid line) decreases towards
zero. As a consequence, determining a value of $\theta^{*}$ such that
$\textrm{Pr}(\theta<\theta^{*}|data)=\alpha$, when $\alpha=0.005$ (the dotted
line) is only possible for $n<20$.

When $\alpha$ is inside the incredibility interval, there remains an
unconventional way for defining a $(1-\alpha)\%$ credible interval. In order
to establish a correct value for $\theta^{*}$ such that
$\textrm{Pr}(\theta<\theta^{*}|data)=\alpha$ (over repeated samples) one
defines $\theta^{*}$ stochastically such that

$\displaystyle\theta^{*}=\begin{cases}0,&\text{with probability
}\gamma;\textrm{ and}\\\ 0+\epsilon,&\text{with probability
}1-\gamma,\end{cases}$ (7)

where:

$\displaystyle\gamma=\frac{\alpha-\textrm{Pr}(\theta\leq
0|data)}{\textrm{Pr}(\theta<0|data)-\textrm{Pr}(\theta\leq 0|data)},$

and $\epsilon$ is an arbitrarily small number.

Returning to our example data with $\bar{y}=1.645/\sqrt{n}$, we note that, for
$n=10$, $\textrm{Pr}(\theta<0|data)=0.030$ and $\textrm{Pr}(\theta\leq
0|data)=0.522$. As such, for $\alpha=0.05$ (which is inside the incredibility
interval of [0.030, 0.522]), we define $\theta^{*}$ as:

$\displaystyle\theta^{*}=\begin{cases}0,&\text{with probability
}\gamma=0.959;\textrm{ and}\\\ 0+\epsilon,&\text{with probability
}(1-\gamma)=0.041.\end{cases}$

Defining $\theta^{*}$ in this way will guarantee that
$\textrm{Pr}(\theta<\theta^{*}|data)=0.05$. One way to think about this is to
consider the various values of $\theta$ that, over a researcher’s lifetime
give rise to the various datasets they analyse. Across all of these studies,
the average posterior probability content of the $[\theta^{*},\infty)$
interval will be 0.95. Moreover, if the model-averaged prior does in fact
correspond to the true data generating mechanism, we can be assured that,
amongst all of the researcher’s studies for which $z=1.645$, 95% of these were
the result of a $\theta$ value from inside of their interval. Furthermore,
since this is true for any arbitrary value of $z$ and any arbitrary value of
$\alpha$, then we have that
$\textrm{Pr}\Big{(}\theta_{j}\in[\theta^{*},\infty)\Big{|}z_{j}\Big{)}=1-\alpha$,
where $\theta_{j}$ and $z_{j}$ are values obtained from a joint draw from the
amalgamation of the prior and statistical model (i.e., the data generating
mechanism).

As another example, suppose $n=100$ and $\bar{y}=2.054/\sqrt{n}=0.2054$ which
corresponds to a $p$-value of $p=0.04$ when using the data to test against the
null hypothesis $\textrm{H}_{0}:\theta=0$, and a $p$-value of $p=0.02$ when
using the data to test against the null hypothesis $\textrm{H}_{0}:\theta<0$.
One can easily calculate an upper one-sided frequentist 95% confidence
interval for these data equal to:
$[\bar{y}-1.645/\sqrt{n},\infty)=[0.040,\infty)$, which clearly excludes 0.
However, one cannot calculate an upper one-sided 95% credible interval since
$\alpha=0.05$ is within the incredibility interval for this data: [0.009,
0.564]. The closest one can do is to calculate an upper one-sided 99.1%
credible equal to: $[0,\infty)$ which includes 0, or calculate an upper one-
sided 43.6% credible interval equal to $(0,\infty)$ which excludes 0. The only
way to define an upper one-sided interval with exactly 95% probability of
including the true value of $\theta$ (over repeated samples) is to do so
stochastically as equal to: $[\theta^{*},\infty)$, where $\theta^{*}=0$ with
probability $\gamma=(0.050-0.564)/(0.009-0.564)=0.926$, and
$\theta^{*}=0+\epsilon$ with probability $1-\gamma=0.074$.

We are not seriously suggesting that researchers define credible intervals in
this bizarre stochastic way. We simply wish to demonstrate that this is the
only way one can correctly define the credible interval from a posterior with
point masses. When model-averaged posteriors involve point-null models,
credible intervals must therefore be approached and interpreted with the
utmost caution. The issue only gets thornier as the sample size increases.

For a very very large $n$ it is possible that both $\alpha/2$ and
$(1-\alpha/2)$ are within the incredibility interval. In this case, the equal-
tailed two-sided $(1-\alpha)\%$ credible interval must be defined in an even
more bizarre way. When both $\alpha/2$ and $(1-\alpha/2)$ are both in the
incredibility interval, the credible interval must be defined stochastically
as either a single point or as an entirely empty interval:

$\displaystyle(1-\alpha)\%\textrm{CrI}=\begin{cases}[0],&\text{with
probability }\psi;\textrm{ and}\\\ \emptyset,&\text{with probability
}(1-\psi),\end{cases}$ (8)

where

$\displaystyle\psi=\frac{\textrm{Pr}(\theta=0|data)-\alpha}{2\times\textrm{Pr}(\theta=0|data)-1}.$

To be clear, the “stochastic credible interval” is not defined in equations
(7) and (8) to ensure that it has a certain (asymptotic) coverage. Rather it
is defined in the only possible way such that (over repeated samples) the
boundaries of the interval contain the correct amount of posterior mass (as
required by the definition in equation (5)). As such, it may not be
immediately obvious that, when we look at the asymptotic behaviour of these
stochastic credible intervals, we see that the Jeffereys-Lindley paradox
reduces the data to be entirely inconsequential (at least when assuming a
fixed $p$-value). Indeed, as $n$ increases, both $\gamma$ and $\psi$ approach
$1-\alpha$ since:

$\displaystyle\textrm{lim}_{n\rightarrow\infty}\gamma$
$\displaystyle=\textrm{lim}_{n\rightarrow\infty}\Big{(}\frac{\alpha-\textrm{Pr}(\theta\leq\theta_{0}|data)}{\textrm{Pr}(\theta<\theta_{0}|data)-\textrm{Pr}(\theta\leq\theta_{0}|data)}\Big{)}$
$\displaystyle=\Big{(}\frac{\alpha-1}{-1}\Big{)}$ $\displaystyle=1-\alpha,$

and:

$\displaystyle\textrm{lim}_{n\rightarrow\infty}\psi$
$\displaystyle=\textrm{lim}_{n\rightarrow\infty}\Big{(}\frac{\textrm{Pr}(\theta=0|data)-\alpha}{2\times\textrm{Pr}(\theta=0|data)-1}\Big{)}$
$\displaystyle=\Big{(}\frac{1-\alpha}{2-1}\Big{)}$ $\displaystyle=1-\alpha.$

Therefore, for sufficiently large $n$ and $z$ remaining constant, the
probability that one will exclude 0 from a $(1-\alpha)$%credible interval will
equal $\alpha$ regardless of the data; see Figure 8. While this may strike one
as paradoxical, it is entirely congruent with the wildly-known consequence of
the Jeffereys-Lindley paradox: As $n$ increases and $z$ is fixed, the
probability of selecting $M_{0}$ will go to 1.

Figure 8:  Each line corresponds to observing data corresponding to a
$p$-value of $p$ when testing against $\textrm{H}_{0}:\theta<0$.

##  4 Conclusion

We demonstrated that when one of the two models under consideration is a
point-null model, not only can  a model-averaged credible interval be rather
different than the frequentist confidence interval, oftentimes it will be
simply undefined (at least in a conventional sense). As a consequence, it may
be tempting to compare (e.g., using the Bayes factor) two a priori probable
models, $M_{0}$ and $M_{1}$, for the purpose of model selection, but then
simply report the uncertainty about $\theta$, conditional on $M_{1}$ being
unquestionably true. We caution that this strategy, while seemingly
straightforward, will lead to unavoidable inconsistencies between one’s priors
and posteriors. Campbell & Gustafson (2022) explain in detail why disregarding
$M_{0}$ “for the purpose of parameter estimation” (Wagenmakers & Gronau 2020)
is inadvisable; see also Tendeiro & Kiers (2019).

Some researchers may be happy to avoid model selection entirely and may see no
reason to entertain point-null priors (e.g., Gelman & Rubin (1995): “realistic
prior distributions in social science do not have a mass of probability at
zero” […] “we believe model selection to be relatively unimportant compared to
the task of constructing realistic models that agree with both theory and
data.”). However, if researchers truly believe that there is a non-zero prior
probability that the parameter of interest is precisely zero (and this prior
probability is equal to the value assigned to $\textrm{Pr}(M_{0})$), Bayesian
testing with a point-null will be optimal in the sense of minimizing the
expected loss (with respect to a joint distribution of the data and
parameters); see Berger (1985). These researchers should be aware that, while
perhaps optimal, Bayesian testing with a point-null can lead to rather
unexpected asymptotic behaviour. There will still be credible intervals; it is
just that, as a consequence of the discontinuity in the model averaged
posterior, certain specific credible intervals do not exist. Some researchers
might therefore wish to explore alternative means of conveying the uncertainty
surrounding the parameter of interest (e.g., Wagenmakers et al. (2022), Rice &
Ye (2022)).

One limitation of this work is that we only considered univariate models where
one wishes to define a credible interval for a single parameter of interest.
However, the ideas we discussed also apply to multivariate settings where one
wishes to define credible sets and where there may be several different nested
models under consideration. For instance, researchers using Bayes factors in
multiple regression models (Rouder & Morey 2012) should be aware that it may
be impossible to define certain model-averaged credible intervals/sets for the
regression coefficients. A second limitation is that we did not consider how
the undefinability of specific credible intervals will also occur in discrete
parameter models. In such cases, specific confidence intervals will also be
undefined (Tingley & Li 1993, Berger 1985), so while they may both be puzzled,
Bayesians and frequentists should at least agree in their inability to define
an uncertainty interval!

Finally, we note that the consequences of the Jeffereys-Lindley paradox on
model selection (and null hypothesis testing) are often understood as
“intuitive” and not necessarily unfavourable: When sample sizes are very
large, researchers might indeed prefer to sacrifice some power in order to
lower the probability of a type I error, a trade-off that occurs necessarily
when testing a point-null hypothesis with the Bayes factor; see Pericchi &
Pereira (2016) and Wagenmakers & Ly (2021). Indeed, the benefits of such a
trade-off are routinely discussed by frequentists and Bayesians alike (e.g.,
Leamer (1978): “from every reasonable viewpoint the significance level should
be a decreasing function of sample size”; and recently, Wulff & Taylor (2023):
“From a Neyman-Pearson perspective, it is logical that $\alpha$ should be a
decreasing function of the sample size.”). However, the consequences of the
Jeffereys-Lindley paradox on parameter estimation –specifically with regards
to model-averaged credible intervals and the inability to define these for
certain probability levels– were previously less well understood, and
certainly strike us as less intuitive.

## References

  * (1)
  * Albers et al. (2018) Albers, C. J., Kiers, H. A. & van Ravenzwaaij, D. (2018), ‘Credible confidence: A pragmatic view on the frequentist vs Bayesian debate’, Collabra: Psychology 4(1).
  * Berger (1985) Berger, J. O. (1985), Statistical decision theory and Bayesian analysis, Springer Science & Business Media.
  * Campbell & Gustafson (2022) Campbell, H. & Gustafson, P. (2022), ‘Bayes factors and posterior estimation: Two sides of the very same coin’, arXiv preprint arXiv:2204.06054 .
  * Casella & Berger (1987) Casella, G. & Berger, R. L. (1987), ‘Reconciling Bayesian and frequentist evidence in the one-sided testing problem’, Journal of the American Statistical Association 82(397), 106–111.
  * Datta & Ghosh (1995) Datta, G. S. & Ghosh, J. K. (1995), ‘On priors providing frequentist validity for Bayesian inference’, Biometrika 82(1), 37–45.
  * Gelman & Rubin (1995) Gelman, A. & Rubin, D. B. (1995), ‘Avoiding model selection in Bayesian social research’, Sociological methodology 25, 165–173.
  * Greenland & Poole (2013) Greenland, S. & Poole, C. (2013), ‘Living with $p$-values: Resurrecting a Bayesian perspective on frequentist statistics’, Epidemiology pp. 62–68.
  * Heck et al. (2022) Heck, D. W., Boehm, U., Böing-Messing, F., Bürkner, P.-C., Derks, K., Dienes, Z., Fu, Q., Gu, X., Karimova, D., Kiers, H. A. et al. (2022), ‘A review of applications of the Bayes factor in psychological research’, Psychological Methods .
  * Held (2020) Held, L. (2020), Bayesian tail probabilities for decision making, in ‘Bayesian Methods in Pharmaceutical Research’, CRC Press Taylor & Francis Group, pp. 53–73.
  * Jeffreys (1935) Jeffreys, H. (1935), Some tests of significance, treated by the theory of probability, in ‘Mathematical proceedings of the Cambridge philosophical society’, Vol. 31, Cambridge University Press, pp. 203–222.
  * Keysers et al. (2020) Keysers, C., Gazzola, V. & Wagenmakers, E.-J. (2020), ‘Using Bayes factor hypothesis testing in neuroscience to establish evidence of absence’, Nature Neuroscience 23(7), 788–799.
  * Leamer (1978) Leamer, E. E. (1978), Specification searches: Ad hoc inference with nonexperimental data, Vol. 53, John Wiley & Sons Incorporated.
  * Lindley (1957) Lindley, D. V. (1957), ‘A statistical paradox’, Biometrika 44(1/2), 187–192.
  * Pericchi & Pereira (2016) Pericchi, L. & Pereira, C. (2016), ‘Adaptative significance levels using optimal decision rules: balancing by weighting the error probabilities’, Brazilian Journal of Probability and Statistics 30(1), 70–90.
  * Rice & Ye (2022) Rice, K. & Ye, L. (2022), ‘Expressing regret: a unified view of credible intervals’, The American Statistician 76(3), 248–256.
  * Rouder & Morey (2012) Rouder, J. N. & Morey, R. D. (2012), ‘Default Bayes factors for model selection in regression’, Multivariate Behavioral Research 47(6), 877–903.
  * Tendeiro & Kiers (2019) Tendeiro, J. N. & Kiers, H. A. (2019), ‘A review of issues about null hypothesis Bayesian testing.’, Psychological Methods 24(6), 774.
  * Tingley & Li (1993) Tingley, M. & Li, C. (1993), ‘A note on obtaining confidence intervals for discrete parameters’, The American Statistician 47(1), 20–23.
  * van den Bergh et al. (2021) van den Bergh, D., Haaf, J. M., Ly, A., Rouder, J. N. & Wagenmakers, E.-J. (2021), ‘A cautionary note on estimating effect size’, Advances in Methods and Practices in Psychological Science 4(1), 2515245921992035.
  * Wagenmakers & Gronau (2020) Wagenmakers, E.-J. & Gronau, Q. F. (2020), ‘Overwhelming evidence for vaccine efficacy in the Pfizer trial: An interim Bayesian analysis’, PsyArXiv.
  * Wagenmakers et al. (2022) Wagenmakers, E.-J., Gronau, Q. F., Dablander, F. & Etz, A. (2022), ‘The support interval’, Erkenn 87, 589–601.
  * Wagenmakers & Ly (2021) Wagenmakers, E.-J. & Ly, A. (2021), ‘History and nature of the Jeffreys-Lindley paradox’, arXiv preprint arXiv:2111.10191 .
  * Wasserstein & Lazar (2016) Wasserstein, R. L. & Lazar, N. A. (2016), ‘The ASA statement on $p$-values: context, process, and purpose’, The American Statistician 70(2), 129–133.
  * Wulff & Taylor (2023) Wulff, J. N. & Taylor, L. (2023), ‘How and why alpha should depend on sample size: A Bayesian-frequentist compromise for significance testing’.

##  5 Appendix

To verify the behavior seen in Figure 2, consider the asymptotic regime with
$\bar{y}=a+bn^{-1/2}$, for some $a$, and some $b>0$. This corresponds to a
fixed $p$-value against $\textrm{H}_{0}:\theta<a$. We are interested in the
posterior probability content of the CI $(\bar{y}-kn^{-1/2},\infty)$ $=$
$(a+(b-k)n^{-1/2},\infty)$, with $k>0$ chosen to give the desired coverage,
i.e., a coverage probability of $\Phi(k)$, where $\Phi()$ is the standard
normal CDF. We presume the prior density $\pi()$ on $\theta$ to be continuous,
with $\pi(a)>0$. In what follows we let Let $\textrm{E}^{*}_{n}$ indicate
expectation with respect to $\theta\sim N(a+bn^{-1/2},n^{-1})$ (so the
posterior on $\theta$ that would arise under a locally uniform prior). And let
$\textrm{E}^{*}$ indicate expectation with respect to $S\sim N(0,1)$ (as one
could get by standardizing the distribution above).  With this set-up, the
posterior probability of interest is the complement of:

$\displaystyle\textrm{Pr}\left(\left.\theta<a+\frac{b-k}{n^{1/2}}\right|\bar{Y}=a+\frac{b}{n^{1/2}}\right)$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{a+(b-k)n^{-1/2}}\pi(\theta|\bar{Y}=a+{b}{n^{-1/2}})d\theta$
$\displaystyle=$
$\displaystyle\frac{\int_{-\infty}^{a+(b-k)n^{-1/2}}\Big{(}f_{Norm}(a+bn^{-1/2},\theta,n^{-1})\times\pi(\theta)\Big{)}d\theta}{\int_{-\infty}^{\infty}\Big{(}f_{Norm}(a+bn^{-1/2},\theta,n^{-1})\times\pi(\theta)\Big{)}d\theta}$
$\displaystyle=$
$\displaystyle\frac{\int_{-\infty}^{a+(b-k)n^{-1/2}}\Big{(}f_{Norm}(\theta,a+bn^{-1/2},n^{-1})\times\pi(\theta)\Big{)}d\theta}{\int_{-\infty}^{\infty}\Big{(}f_{Norm}(\theta,a+bn^{-1/2},n^{-1})\times\pi(\theta)\Big{)}d\theta}$
$\displaystyle=$
$\displaystyle\frac{\textrm{E}^{*}_{n}\left\\{I_{(-\infty,a+(b-k)n^{-1/2})}(\theta)\pi(\theta)\right\\}}{\textrm{E}^{*}_{n}\left\\{\pi(\theta)\right\\}}$
$\displaystyle=$
$\displaystyle\frac{\textrm{E}^{*}\left\\{I_{(-\infty,-k)}(S)\pi\left(a+(b+S)n^{-1/2}\right)\right\\}}{\textrm{E}^{*}\left\\{\pi\left(a+(b+S)n^{-1/2}\right)\right\\}}$
$\displaystyle\xrightarrow[n\to\infty]{}$
$\displaystyle\frac{\pi(a)\textrm{E}^{*}\left\\{I_{(-\infty,-k)}(S)\right\\}}{\pi(a)}$
$\displaystyle=$ $\displaystyle\Phi(-k),$

as claimed. Specifically, n Figure 2, we have $\pi(\theta)$ defined as the
“mixture” prior (defined in equations 2 and 3) and set $a=0$ and $b=1.645$,
with 4 different values of $k$: $k=1.645$, $k=1.282$, $k=0.842$, $k=0.524$,
such that $\Phi(-1.645)=0.05$, $\Phi(-1.282)=0.10$, $\Phi(-0.842)=0.20$, and
$\Phi(-0.524)=0.30$.