#  Singular 2-webs and Polar Curves

Fernando Etayo111Departamento de Matemáticas, Estadística y Computación.
Facultad de Ciencias. Universidad de Cantabria. Avda. de los Castros, s/n,
39071 Santander, SPAIN. e-mail<EMAIL_ADDRESS>
###### Abstract

A 2-web in the plane is given by two everywhere transverse 1-foliations. In
this paper we introduce the study of singular 2-webs, given by any two
foliations, which may be tangent in some points. We show that such two
foliations are tangent along a curve, which will be called the polar curve of
the 2-web, and we study the relationship between the contact order of leaves
of both foliations and the singularities of the polar curve.

AMS Classification: 53A60, 53C15.

Keywords: Singular 2-web, Foliations, Pfaffian forms, Paracomplex structure.

##  1 Introduction

A 2-web in the real plane is given by two everywhere transverse 1-dimensional
foliations, i.e., two families of curves such that in any point the curves
passing through it have different tangent lines. These 2-webs are always
locally diffeomorphic to that of vertical and horizontal lines, and then they
have no local invariants. The case of 3-webs is completely different because
of the curvature of the Blaschke-Chern connection, which measures how far the
web is to be hexagonal [2, 6].

A 2-web can be defined by a (1,1) tensor field $F$ of maximum rank such that
$F^{2}=I$, where $I$ denotes the identity tensor field. The eigenspaces
associated to the eigenvalues $\pm 1$ define two distributions, which are
involutive in the case of the real plane. $F$ is said to be a paracomplex
structure.

One could think that 2-webs have no interest. Nevertheless, we want to go into
an unknown landscape, which has not been explored yet, as far as the author
knows. Let us consider two different 1-dimensional foliations in the plane.
What can you say about the set of points in which both foliations are tangent?
If this set is non-empty, we shall say that both foliations define a _singular
2-web_ and the set will be called the _polar curve_ of the singular 2-web. We
use this terminology following a similar idea [3] developed in the context of
a holomorphic foliation in $\mathbb{C}^{2}$, where the polar curve is defined
from the foliation and a direction in the plane, i.e., it is defined from
2-web defined by the holomorphic foliation and the foliation of lines parallel
to the direction.

Polar curves should not be confused with polar foliations. A polar foliation
[1] is a singular foliation in a complete Riemannian manifold such that for
each regular point $p$, there is an immersed submanifold $\Sigma_{p}$, called
section, that passes through $p$ and that meets all the leaves and always
perpendicularly.

In the previous paper [4] we have shown some explicit examples. In the present
one we state some results about the polar curve. The following example is
introduced as a motivating case.

###### Example 1

Let us consider the punctured plane $\mathbb{R}^{2}-\\{(0,0)\\}$ and the
foliations given by

  * •

${\cal F}_{1}$ is the set of circles with center $(0,0)$.

  * •

${\cal F}_{2}$ is the set of vertical lines.

Obviously, this is a singular 2-web having the horizontal axis as polar curve.
We can prove it by using different techniques: analyzing the associated
distributions, the paracomplex structure and the Pfaffian forms. As is well
known 1-dimensional distributions are always integrable, and then we can work
interchangeably with distributions and foliations. We shall show carefully
these three techniques, because we should choose the best one in order to
ahead more complex situations.

  * •

Associated distributions: The tangent vector of the foliation ${\cal F}_{1}$
at the point $(x,y)$ is $X=-y\frac{\partial}{\partial
x}+x\frac{\partial}{\partial y}$. The tangent vector to ${\cal F}_{2}$ is
$Y=\frac{\partial}{\partial y}$. These vectors are colinear in the polar curve
$\\{y=0\\}-\\{(0,0)\\}$.

  * •

The paracomplex structure [5] is the $(1,1)$-tensor field $F$ such that their
eigenspaces are the distributions tangent to the foliations. A straightforward
calculation shows that

$F=\left(\begin{array}[]{cc}1&0\\\ -\frac{2x}{y}&-1\end{array}\right)$

One can esaily check that $F^{2}=id;F(X)=X;F(Y)=-Y$. The paracomplex structure
$F$ is well defined in all the punctured plane unless the polar curve
$\\{y=0\\}-\\{(0,0)\\}$. Then the polar curve appears as the set of points
where the structure cannot be defined.

  * •

The Pfaffian forms defining ${\cal F}_{1}$ and ${\cal F}_{2}$ are
$\omega=x\,dx+y\,dy$ and $\eta=dx$. The points where these 1-forms are
dependent are those where $0=\omega\wedge\eta=(x\,dx+y\,dy)\wedge
dx=-ydx\wedge dy$, which are those of the polar curve.

The best option is to work with Pfaffian forms, because one has only to check
a product of 1-forms.

##  2 Algebraic foliations and paracomplex structure

As is well known, a foliation is said to be algebraic if it is given by a
1-form $\omega=\omega_{1}(x,y)dx+\omega_{2}(x,y)dy$, where $\omega_{i}(x,y)$
are polynomials. This doesn’t mean that the algebraic curves of the foliation
have to be algebraic curves. For example, the 1-form $\omega=ydx-dy$ is
algebraic and the curves of this foliation are the exponential
$C_{k}=\\{y=ke^{x}\\},k\in\mathbb{R}$.

In the same way, we can say that a (1,1)-tensor field
$F=\frac{\partial}{\partial x^{i}}\otimes F^{i}_{j}dx^{j}$ is algebraic (resp.
rational) if the functions $F^{i}_{j}$ are polynomial (resp. rational).

Then we have,

###### Theorem 2

Let $\omega=\omega_{1}(x,y)dx+\omega_{2}(x,y)dy$ and
$\eta=\eta_{1}(x,y)dx+\eta_{2}(x,y)dy$ be two 1-forms in an open subset of the
real plane. Then,

(1) The polar curve of the singular 2-web defined by $\omega$ and $\eta$ is
the curve $\\{\omega_{1}\eta_{2}-\omega_{2}\eta_{1}=0\\}$. Besides, if the
coefficient functions $\omega_{1},\omega_{2},\eta_{1},\eta_{2}$ are polynomial
of degrees $p,q,r,s$ then the polar curve is an algebraic curve of degree less
or equal to ${\rm max}(p+s,q+r)$.

(2) The paracomplex structure $F$ defined by $\omega$ and $\eta$ is

$F=\frac{1}{-\omega_{2}\eta_{1}+\omega_{1}\eta_{2}}\left(\begin{array}[]{cc}-\eta_{1}\omega_{2}-\omega_{1}\eta_{2}&-2\omega_{2}\eta_{2}\\\
2\omega_{1}\eta_{1}&\eta_{2}\omega_{1}+\omega_{2}\eta_{1}\end{array}\right)$

(3) If $\omega$ and $\eta$ are algebraic, then $F$ is rational.

_Proof_.

(1) Observe that
$\omega\wedge\eta=(\omega_{1}\eta_{2}-\omega_{2}\eta_{1})dx\wedge dy$, thus
showing that both 1-forms are dependent on the curve
$C=\\{\omega_{1}\eta_{2}-\omega_{2}\eta_{1}=0\\}$.

(2) In the points $p\in\mathbb{R}^{2}-\\{C\\}$, where $c$ denotes the polar
curve, one has $T_{p}\mathbb{R}^{2}=(ker\omega)_{p}\oplus(ker\eta)_{p}$. A
basis of $ker\omega$ is $-\omega_{2}\frac{\partial}{\partial
x^{1}}+\omega_{1}\frac{\partial}{\partial x^{2}}=(-\omega_{2},\omega_{1})$. A
basis of $ker\eta$ is $(-\eta_{2},\eta_{1})$. Then, a vector field
$v=v_{1}\frac{\partial}{\partial x^{1}}+v_{2}\frac{\partial}{\partial x^{2}}$
can be decomposed as

$v=\left(\begin{array}[]{c}v_{1}\\\
v_{2}\end{array}\right)=\alpha\left(\begin{array}[]{c}-\omega_{2}\\\
\omega_{1}\end{array}\right)+\beta\left(\begin{array}[]{c}-\eta_{2}\\\
\eta_{1}\end{array}\right)$

which produces, by using the Cramer rule,

$\alpha=\frac{\mid\begin{array}[]{cc}v_{1}&-\eta_{2}\\\
v_{2}&\eta_{1}\end{array}\mid}{\mid\begin{array}[]{cc}-\omega_{2}&-\eta_{2}\\\
\omega_{1}&\eta_{1}\end{array}\mid}\hskip 14.22636pt;\hskip
14.22636pt\beta=\frac{\mid\begin{array}[]{cc}-\omega_{2}&v_{1}\\\
\omega_{1}&v_{2}\end{array}\mid}{\mid\begin{array}[]{cc}-\omega_{2}&-\eta_{2}\\\
\omega_{1}&\eta_{1}\end{array}\mid}$

Then, the paracomplex structure $F$ defined as $F|_{ker\omega}=I$,
$F|_{ker\eta}=-I$ is that given by

$F(v)=\alpha\left(\begin{array}[]{c}-\omega_{2}\\\
\omega_{1}\end{array}\right)-\beta\left(\begin{array}[]{c}-\eta_{2}\\\
\eta_{1}\end{array}\right)=\frac{1}{-\omega_{2}\eta_{1}+\omega_{1}\eta_{2}}\left(\begin{array}[]{cc}-\eta_{1}\omega_{2}-\omega_{1}\eta_{2}&-2\omega_{2}\eta_{2}\\\
2\omega_{1}\eta_{1}&\eta_{2}\omega_{1}+\omega_{2}\eta_{1}\end{array}\right)\left(\begin{array}[]{c}v_{1}\\\
v_{2}\end{array}\right)$

(3) It is a direct consequence of (2).

###### Remark 3

Observe that the polar curve corresponds to the locus where the paracomplex
structure $F$ cannot be defined. Observe that formula (2) does not depend on
the algebricity of foliations.

###### Remark 4

If one changes the 1-forms $\omega$ and $\eta$ by proportional 1-forms
$f\omega$ and $g\eta$, f,g beings functions, then the 2-web is the same, and
the paracomplex structure $F$ remains invariable in the above theorem. This is
important: the 1-forms are not uniquely defined, but the paracomplex structure
is uniquely defined, up to sign.

##  3 Singular points of the polar curve

Let us assume that $\omega$ and $\eta$ are two algebraic foliations with polar
curve $C$. As $C$ is an algebraic curve, $C$ may have singular points. The
following examples suggest that singular points of $C$ correspond to the
points where the curves of each foliations have contact of order greater than
two. Remember the classical definitions:

  * •

Two curves $\alpha$ and $\beta$ are said to have _contact of order_ $k$ at a
point $p$ if their derivatives of order $0,1,\ldots,k$ coincide at the point
and the derivatives of order $k+1$ are different. We shall denote it as
$ord_{\alpha\beta}(p)=k$.

  * •

The multiplicity $mult_{C}(p)$ of a curve C at p is the order of the first
non-vanishing term in the Taylor expansion of f at p, where
$C=\\{f(x,y)=0\\}$. The point is said to be a _regular_ point if
$mult_{C}(p)=1$, and _singular_ if $mult_{C}(p)\geq 2$.

  * •

A 1-form $\omega=\omega_{1}(x,y)dx+\omega_{2}(x,y)dy$ is said to define an
_exact differential equation_ if $\frac{\partial\omega_{1}(x,y)}{\partial
y}=\frac{\partial\omega_{2}(x,y)}{\partial x}$. In this case, a curve
$f(x,y)=0$ is an integral curve iff $\frac{\partial f}{\partial
x}=\omega_{1}(x,y)$ and $\frac{\partial f}{\partial y}=\omega_{2}(x,y)$. As it
is well known, multiplying by an integrating factor $\mu(x,y)$ any 1-form can
be transformed into an exact differential equation, although in many cases
obtaining the integrating factor is not easy.

###### Example 5

Let us consider the foliations ${\cal F}_{\omega}$ and ${\cal F}_{\eta}$ given
by the 1-forms and their dual vector fields:

$\omega=y^{2}dx-dy,\;X_{\omega}=\frac{\partial}{\partial
x}+y^{2}\frac{\partial}{\partial y};$

$\eta=-x^{3}dx+dy,\;X_{\eta}=-\frac{\partial}{\partial
x}-x^{3}\frac{\partial}{\partial y}$.

The first one is the set of hyperbolas $\\{y=\frac{1}{k-x},k\in\mathbb{R}\\}$
and the horizontal axis, and the second one that of quartics
$\\{y=k+\frac{x^{4}}{4},k\in\mathbb{R}\\}$. The polar curve is the cubic
$C=\\{y^{2}=x^{3}\\}$, which has a singular point of multiplicity two (a cusp
point) in the origin. The curves of ${\cal F}_{\omega}$ and ${\cal F}_{\eta}$
through the origin are $\\{y=0\\}$ and $\\{y=\frac{x^{4}}{4}\\}$. The last one
has an inflection point in the origin, thus proving that both curves have
contact of order three.

###### Example 6

Let us consider the foliation ${\cal F}_{\alpha}$ given by the 1-form
$\alpha=dy$, and the foliation ${\cal F}_{\eta}$ given by the same 1-form
$\eta$ of the above example. Then de polar curve is $C=\\{x^{3}=0\\}$, which
is the vertical axis, being all of its points singular of multiplicity three,
equal to the contact order of tangent curves of each foliation.

We can state:

###### Theorem 7

Let ${\cal F}_{\omega}$ and ${\cal F}_{\eta}$ be two algebraic foliations
given by two exact differential equations $\omega$ and $\eta$ and let $C$ be
their polar curve. If the curves $\gamma_{\omega}$ and $\gamma_{\eta}$ of
${\cal F}_{\omega}$ and ${\cal F}_{\eta}$ through $p$ have contact of order
$\geq 2$ then $p\in c$ is a singular point of the polar curve. In this case,
$ord_{\gamma_{\omega}\gamma_{\eta}}(p)\geq mult_{C}(p)$.

_Proof_

First we introduce some notation. Let
$\omega=\omega_{1}(x,y)dx+\omega_{2}(x,y)dy$ and
$\eta=\eta_{1}(x,y)dx+\eta_{2}(x,y)dy$ the 1-forms corresponding to ${\cal
F}_{\omega}$ and ${\cal F}_{\eta}$ and let

$C=\\{\omega_{1}\eta_{2}-\omega_{2}\eta_{1}=0\\}=\\{f(x,y)=0\\}$

be their polar curve. Then a point $p\in C$ is singular iff

$\frac{\partial f}{\partial x}(p)=0=\frac{\partial f}{\partial y}(p)$

Thus, $p$ is a singular point iff the following equations hold at $p$:

(1) $\left\\{\begin{array}[]{cc}\frac{\partial\omega_{1}}{\partial
x}\eta_{2}+\omega_{1}\frac{\partial\eta_{2}}{\partial
x}-\frac{\partial\omega_{2}}{\partial
x}\eta_{1}-\omega_{2}\frac{\partial\eta_{1}}{\partial x}&=0\\\ &\\\
\frac{\partial\omega_{1}}{\partial
y}\eta_{2}+\omega_{1}\frac{\partial\eta_{2}}{\partial
y}-\frac{\partial\omega_{2}}{\partial
y}\eta_{1}-\omega_{2}\frac{\partial\eta_{1}}{\partial
y}&=0\end{array}\right\\}$

Let $\gamma_{\omega}$ and $\gamma_{\eta}$ the integral curves of ${\cal
F}_{\omega}$ and ${\cal F}_{\eta}$ through the point $p$. As the 1-forms
$\omega$ and $\eta$ are exact, we have:

(2) $\left\\{\begin{array}[]{cc}\frac{\partial\gamma_{\omega}}{\partial
x}&=\omega_{1}\\\ &\\\ \frac{\partial\gamma_{\omega}}{\partial
y}&=\omega_{2}\end{array}\right\\}\hskip 28.45274pt{\rm and}\hskip
28.45274pt\left\\{\begin{array}[]{cc}\frac{\partial\gamma_{\eta}}{\partial
x}&=\eta_{1}\\\ &\\\ \frac{\partial\gamma_{\eta}}{\partial
y}&=\eta_{2}\end{array}\right\\}$

Then, $ord_{\gamma_{\omega}\gamma_{\eta}}\geq 2$ iff the following equations
hold at $p$:

(3) $\left\\{\begin{array}[]{ccc}\frac{\partial\gamma_{\omega}}{\partial
x}=\frac{\partial\gamma_{\eta}}{\partial
x};&\frac{\partial\gamma_{\omega}}{\partial
y}=\frac{\partial\gamma_{\eta}}{\partial y}&\\\ &&\\\
\frac{\partial^{2}\gamma_{\omega}}{\partial
x^{2}}=\frac{\partial^{2}\gamma_{\eta}}{\partial
x^{2}};&\frac{\partial^{2}\gamma_{\omega}}{\partial x\partial
y}=\frac{\partial^{2}\gamma_{\eta}}{\partial x\partial
y};&\frac{\partial^{2}\gamma_{\omega}}{\partial
x^{y}}=\frac{\partial^{2}\gamma_{\eta}}{\partial x^{y}}\end{array}\right\\}$

Then, taking into account equation (2) one easily check that
$(3)\Rightarrow(1)$, thus proving the first statement of the theorem. The
second one follows from a similar reasoning.

The condition in the above theorem of being $\omega$ and $\eta$ exact
differential equations is necessary, as the following example in
$\mathbb{R}^{2}-\\{(\pm 1,0)\\}$ shows:

###### Example 8

Let ${\cal F}$ be the vertical axis and the set of circles $C_{a}$, with
center in the horizontal axis passing through the points $(a,0)$ and
$(1/a,0)$, where $a\neq\pm 1$.

Points $(a,0),(1/a,0),(1,0),(-1,0)$ define a harmonic quadruple of points,
i.e., their cross ratio is -1. The center of $C_{a}$ is the point
$\frac{a^{2}+1}{2a}$ and the radius is $\frac{a^{2}-1}{2a}$ thus obtaining the
equation

$C_{a}=\left\\{\left(x-\frac{a^{2}+1}{2a}\right)^{2}\,+y^{2}\,=\,\left(\frac{a^{2}-1}{2a}\right)^{2}\right\\}=\left\\{x^{2}-\frac{a^{2}+1}{a}\,x+y^{2}+1=0\right\\}$

Let us denote by $f(x,y)=0$ is the equation of the curve $C_{a}$. The tangent
line in $p=(x,y)\in C_{a}$ is $f_{x}(x-p_{1})+f_{y}(y-p_{2})$, whose direction
is generated by the vector $(-f_{y}(p),f_{x}(p))$, and the corresponding dual
form will be $f_{x}(p)dx+f_{y}(p)dy$. In our case,

$\omega=\left(2x-\frac{a^{2}+1}{a}\right)\,dx+\,2y\,dy$

Then, for the equation of $C_{a}$, we can deduce that

$\frac{a^{2}+1}{a}=\frac{x^{2}+y^{2}+1}{x}$

thus showing

$\omega=\frac{x^{2}-y^{2}-1}{x}\,dx\,+\,2y\,dy\>;\>{\rm if}\,x\neq 0$

We can replace the Pffaf form by another one obtained multiplying by the
function $\mu(x,y)=x$, and then we would have coefficients of degree two.
Then, we can take

$\omega=(x^{2}-y^{2}-1)\,dx\,+\,2xy\,dy$

which is not an exact differential equation.

Let us consider the singular 2-web given by the foliations:

  * •

${\cal F}_{\omega}$ is the vertical axis and the set of circles $C_{a}$, with
center in the horizontal axis passing through the points $(a,0)$ and
$(1/a,0)$.

  * •

${\cal F}_{\eta}$ is the set vertical lines, whose Pfaffian form is $\eta=dx$.

Then the polar curve is the reducible curve $\\{xy=0\\}$ given by both axis
and has a unique singular point. Nevertheless, both foliations have in common
the leaf $\\{x=0\\}$, and then in all of its points have contact of order
infinite, thus showing that there are regular points of the polar curve
corresponding to higher order contact of both foliations.

On the other hand, by using Theorem 2 one can obtain the paracomplex structure
associated to this 2-web:

$F=\left(\begin{array}[]{cc}1&0\\\
\frac{-(x^{2}-y^{2}-1)}{xy}&-1\end{array}\right)$

thus showing again that the polar curve is $\\{xy=0\\}$. But one cannot obtain
information about the contact order of leaves of both foliations.

In the present case, as one can easily check, $\mu(x,y)=\frac{1}{x^{2}}$ is an
integrating factor of $\omega$, and then we can re-write $\omega$ as

$\omega=\frac{x^{2}-y^{2}-1}{x^{2}}\,dx\,+\,\frac{2y}{x}\,dy$

Then, the equation of the polar curve is $C=\\{-\frac{2y}{x}=0\\}$, showing
the special property of the vertical axis $\\{x=0\\}$. The points where the
polar curve cannot be defined when it is obtained from exact differential
equations corresponds with those of a common leaf of both foliations. But this
example also shows that multiplying the 1-form by an integrating factor can
add points where the 1-form is not defined: in the example those of the
vertical axis $\\{x=0\\}$.

##  4 Conclusions

We write down some global conclusions.

  * •

The following mathematical objects are equivalent in the real plane: a 2-web,
two differential equations, two vector fields, two 1-forms, two distributions.
Assume that each of them has no singularities nor zeros (restricting to an
open subset of the plane, if necessary). Then, we have two families of curves,
and we ask where they are tangent. This is a natural question and, as far as
the author knows, there was not yet any answer about it.

  * •

The points where the foliations are tangent define a curve, called the polar
curve of the 2-web. If one defines the 2-web by means of two 1-forms it is
very easy to find an equation for the polar curve. If the 1-forms are
algebraic, the polar curve is an algebraic curve.

  * •

The 1-forms associated to a 2-web are not unique. If one takes exact
differential equations for them, then we have proved that higher order contact
points of the foliations are singular points of the polar curve.

  * •

There exist integrating factors which allow to obtain an exact differential
equation for any 1-form. In general, integrating factors are difficult to be
calculated. Besides they can exclude points of the plane.

  * •

The paracomplex structure is unique up to a sign and the points where it is
not defined define the polar curve. One can derive the expression of the
paracomplex structure from those of the 1-forms, but there is no a general way
for the reverse.

## References

  * [1] Alexandrino, M. M.: On polar foliations and fundamental group. _Results in Mathematics_ , 60, Issue 1 (2011), 213–223
  * [2] Blaschke, W.: _Geometrie der Gewebe_. Springer, 1938.
  * [3] Corral, N.: Infinitesimal adjunction and polar curves. _Bull. Braz. Math. Soc._ (N.S.) 40 (2009), no. 2, 181 -224.
  * [4] Etayo, F: Singular 2-webs. An Introduction. _Contribuciones matemáticas en homenaje a Juan Tarrés_. Ed. Complutense, 2012, 141-147. ISBN: 978-84-695-4421-1.
  * [5] Etayo, F. ; Santamaría, R; Trías, U. R.: The geometry of a bi-Lagrangian manifold. _Differential Geom. Appl_. 24 (2006), 33–59.
  * [6] Grifone, J.; Salem, E.: _Web theory and related topics_. World Scientific Publishing Co. Pte. Ltd. (2001)


# Extreme data compression for Bayesian model comparison

Alan F. Heavens,11footnotetext: Corresponding author. Arrykrishna Mootoovaloo
Roberto Trotta  Elena Sellentin

###### Abstract

We develop extreme data compression for use in Bayesian model comparison via
the MOPED algorithm, as well as more general score compression. We find that
Bayes factors from data compressed with the MOPED algorithm are identical to
those from their uncompressed datasets when the models are linear and the
errors Gaussian. In other nonlinear cases, whether nested or not, we find
negligible differences in the Bayes factors, and show this explicitly for the
Pantheon-SH0ES supernova dataset. We also investigate the sampling properties
of the Bayesian Evidence as a frequentist statistic, and find that extreme
data compression reduces the sampling variance of the Evidence, but has no
impact on the sampling distribution of Bayes factors. Since model comparison
can be a very computationally-intensive task, MOPED extreme data compression
may present significant advantages in computational time.

##  1 Introduction

Astronomical datasets can be very large, with galaxy catalogues and weak
lensing surveys having up to hundreds of millions of objects [e.g., 1, 2, 3],
and Planck microwave background data [4] containing billions of time-ordered
data. Extracting scientific knowledge from such large datasets almost
inevitably requires some form of data compression, typically into summary
statistics such as point estimates of two-point summaries, either correlation
functions or power spectra. Compression without loss of information is
generally possible since the datasets are noisy and correlated, and
furthermore the information may be captured by a small number of summary
statistics, depending on the statistical properties of the underlying data.
However, even with compression to two-point summaries, the number of
statistics can still be large: for example, cosmic shear surveys such as
proposed for the Euclid mission [5] may generate $\sim 10^{4}$ two-point
statistics, and this can present problems for likelihood-based scientific
analysis if the covariance of the summaries needs to be simulated, or indeed
if their distributions are not known [6]. Both of these considerations
motivate the use of more radical data compression techniques. This can also be
highly beneficial when the resultant summaries are coupled with techniques
such as simulation-based inference (SBI, also known as likelihood-free
inference), which is undergoing rapid development thanks in part to novel
machine learning techniques [e.g., 7, 8, 9, 10, 11]. Typically, the most
extreme compression possible without losing information is to reduce the
dataset down from its original size to the number of model parameters to be
inferred. Such extreme compression can have a number of benefits, in terms of
speed of analysis, in reducing the number of simulations required to determine
the sampling distribution of the statistics [6], and in making SBI feasible
[12, 13], whilst still allowing for new models to be considered [14]. Heavens
et al. [15] derived the optimal data compression algorithm MOPED, for
gaussian-distributed data whose mean signal depends on the parameters, showing
that extreme data compression can preserve the Fisher matrix, so the
compression is, in this sense, lossless. Techniques, based on Karhünen-Loève
compression [16] also exist for data whose covariance depends on the
parameters, but not the mean, in contrast to MOPED, which applies when the
mean depends on the parameters, but not the covariance. However, the former
situation does not lead to lossless extreme compression, nor does it
generalise neatly to multiple parameters. MOPED, whether in its original Gram-
Schmidt orthogonal form, or in the simpler non-orthogonalised version (which
we use in this paper) is a special case of using the score function (the
gradient of the likelihood with respect to the model parameters) to compress
the data, as shown by [17]. Alternative approaches to extreme data compression
include information-maximising neural networks (IMNNs) to find maximally-
informative summaries, as introduced by [7], and which have in certain cases
been shown to be lossless at Fisher level [9]. Typically, the neural network-
based methods operate on the fields themselves [e.g. 10, 9, 11], whereas the
MOPED and score compressions operate on some intermediate, physics-based
summaries [e.g. 6, 10, 18, 12].

The main motivation for this paper is to determine whether the qualities of
optimal extreme data compression which apply to parameter inference also
extend to Bayesian Evidence and Bayes factors, i.e., to model comparison. The
answer to this question is yes, provided that the extreme data compression is
performed in each of the models separately, and the compressed data then
combined, such that the combined compressed dataset is by construction optimal
in all models. This applies to both nested and non-nested models.

With this key result, committed Bayesians may stop reading, but we also take
the opportunity to investigate the sampling distribution of the Bayesian
Evidence and the related Bayes Factor, to look at their use as frequentist
statistics under extreme compression. The motivation for this comes from a
suggestion [19] that data compression may help in reducing the sampling
variability in the Bayesian Evidence. The answer is that it does, but the
sampling variance of the Bayes Factor, which is a ratio of evidences and the
relevant quantity for model comparison, remains virtually unchanged. To
rephrase this same point, one can use extreme MOPED data compression as
effectively as the full dataset for model comparison as well as parameter
inference. Since model comparison is typically a computationally expensive
exercise, extreme data compression, which often results in faster likelihood
evaluation, can be very valuable.

Note that the distinction between the sampling distribution of the Bayesian
Evidence and of the Bayes Factor is important, especially since it is tempting
to conclude –erroneously– that the sampling distribution of the Bayesian
Evidence is sufficiently wide to compromise its usefulness in model selection,
but this ignores the fact that much of the variation is common to both models
(the data are the same), and cancels out when the ratio of Bayesian Evidence
is computed, leaving the Bayes factor as an effective discriminant, as we show
in this paper.

The layout of the paper is as follows. In section 2 we review the MOPED data
compression algorithm, and show how it works effectively for model comparison
as well as parameter inference, and consider more general score compression.
In section 3.1 we explore the frequentist properties of Bayes factors, showing
that the sampling distribution of the Bayes factor is unaltered whether the
full or compressed data are used. In section 3.2 we apply MOPED compression to
the Pantheon+SH0ES dataset [20]) and show that both the cosmological
parameters and the Bayes Factor can be determined from just 3 numbers (the
MOPED coefficients), as precisely as from the original 1590 supernovae. We
conclude in section 4.

##  2 Extreme data compression: Bayesian Evidence with MOPED

###  2.1 Bayesian Evidence and Bayes Factors for the uncompressed dataset

Let us first consider gaussian-distributed data, with fixed covariance matrix
and arbitrary dependence of the mean on the parameters. Data compression with
MOPED [15] allows extreme compression to $p$ summary statistics, where $p$ is
the number of model parameters, whatever the length of the original data
vector. Although the resulting MOPED coefficients are not strictly sufficient
statistics except in the case of linear models, they can preserve the entire
Fisher matrix, so the likelihood surface near the peak is no broader than when
using the full, uncompressed dataset. As a result, the extreme compression can
be used for parameter inference without increasing errors, and has been used
in many contexts, such as star formation histories of SDSS galaxies [e.g. 21,
22, 23, 24], the cosmic microwave background [e.g. 25, 26, 27], galaxy power
spectrum and bispectrum [e.g. 28, 29], gravitational waves [30], weak lensing
[e.g. 12, 31], combined probes [e.g. 32, 33], planetary transits [34], and
extended to parameter-dependent covariance matrices [6] and mis-specified
models [14]. In this section, we explore how MOPED and related score
compression can be used in a new context: model comparison. To begin we review
the derivation of the unorthogonalised version of MOPED, following not the
original analysis [15], but rather the derivation in [17].

We consider a model with parameters
${\boldsymbol{\theta}}=\\{\theta_{1},\dots,\theta_{p}\\}$, data
${\boldsymbol{x}}\in\mathbb{R}^{n}$ (where $n\gg p$) generated from a gaussian
sampling distribution with covariance matrix ${\sf{C}}\in\mathbb{R}^{n\times
n}$ independent of ${\boldsymbol{\theta}}$ and expectation value of the data
${\boldsymbol{\mu}}({\boldsymbol{\theta}})$. The sampling distribution for the
uncompressed dataset is therefore

$p({\boldsymbol{x}}|{\boldsymbol{\theta}})=\frac{1}{\sqrt{|2\pi{\sf{C}}|}}\exp\left[-\frac{1}{2}({\boldsymbol{x}}-{\boldsymbol{\mu}})^{T}{\sf{C}}^{-1}({\boldsymbol{x}}-{\boldsymbol{\mu}})\right].$
(2.1)

Letting ${\boldsymbol{X}}\equiv{\boldsymbol{x}}-{\boldsymbol{\mu}}_{*}$, where
${\boldsymbol{\mu}}_{*}$ is the expectation value of the data at some fiducial
set of parameters ${\boldsymbol{\theta}}_{*}$, the posterior is

$p({\boldsymbol{\theta}}|{\boldsymbol{x}})\propto\frac{1}{\sqrt{|2\pi{\sf{C}}|}}\exp\left[-\frac{1}{2}({\boldsymbol{X}}-\tilde{\boldsymbol{\mu}})^{T}{\sf{C}}^{-1}({\boldsymbol{X}}-\tilde{\boldsymbol{\mu}})\right]\,\pi({\boldsymbol{\theta}})$
(2.2)

for a prior density $\pi({\boldsymbol{\theta}})$. Here

$\tilde{\boldsymbol{\mu}}({\boldsymbol{\theta}})\equiv{\boldsymbol{\mu}}({\boldsymbol{\theta}})-{\boldsymbol{\mu}}_{*},$
(2.3)

and a linear Taylor expansion around the fiducial point gives

$\tilde{\boldsymbol{\mu}}\simeq{\boldsymbol{\mu}}_{,\alpha}\tilde{\boldsymbol{\theta}}_{\alpha}\quad{\rm
where}\quad\tilde{\boldsymbol{\theta}}\equiv{\boldsymbol{\theta}}-{\boldsymbol{\theta}}_{*}.$
(2.4)

The comma indicates $\partial/\partial{\boldsymbol{\theta}}_{\alpha}$ and the
summation convention is used. Hence to this order,

$\displaystyle p({\boldsymbol{\theta}}|{\boldsymbol{x}})\propto$
$\displaystyle\frac{1}{\sqrt{|2\pi{\sf{C}}|}}\exp\left(-\frac{1}{2}{\boldsymbol{X}}^{T}{\sf{C}}^{-1}{\boldsymbol{X}}\right)\exp\left(\tilde{\boldsymbol{\mu}}^{T}{\sf{C}}^{-1}{\boldsymbol{X}}-\frac{1}{2}\tilde{\boldsymbol{\mu}}^{T}{\sf{C}}^{-1}\tilde{\boldsymbol{\mu}}\right)\,\pi({\boldsymbol{\theta}}).$
(2.5)

We see that the only coupling of the data ${\boldsymbol{x}}$ to the parameters
is through the $\tilde{\boldsymbol{\mu}}^{T}{\sf{C}}^{-1}{\boldsymbol{X}}$
term, and only via the $p$ MOPED coefficients

$Y_{\alpha}\equiv{\boldsymbol{\mu}}_{,\alpha}^{T}{\sf{C}}^{-1}{\boldsymbol{X}}\equiv{\boldsymbol{b}}_{\alpha}^{T}({\boldsymbol{x}}-{\boldsymbol{\mu}}_{*}),\qquad\alpha=1\ldots
p.$ (2.6)

and we identify the $p$ (unorthogonalised) MOPED weight vectors [15])

${\boldsymbol{b}}_{\alpha}^{T}\equiv{\boldsymbol{\mu}}_{,\alpha}^{T}{\sf{C}}^{-1}.$
(2.7)

The posterior may then be written

$\displaystyle p({\boldsymbol{\theta}}|{\boldsymbol{x}})$
$\displaystyle\propto\frac{1}{\sqrt{|2\pi{\sf{C}}|}}\exp\left(-\frac{1}{2}{\boldsymbol{X}}^{T}{\sf{C}}^{-1}{\boldsymbol{X}}\right)\exp\left[({\boldsymbol{b}}_{\alpha}^{T}{\boldsymbol{X}})\tilde{\theta}_{\alpha}-\frac{1}{2}({\boldsymbol{b}}_{\alpha}^{T}{\sf{C}}{\boldsymbol{b}}_{\beta})\tilde{\theta}_{\alpha}\tilde{\theta}_{\beta}\right]\,\pi(\tilde{\boldsymbol{\theta}}),$
(2.8)

or in matrix form

$\displaystyle p({\boldsymbol{\theta}}|{\boldsymbol{x}})$
$\displaystyle\propto\frac{1}{\sqrt{|2\pi{\sf{C}}|}}\exp\left(-\frac{1}{2}{\boldsymbol{X}}^{T}{\sf{C}}^{-1}{\boldsymbol{X}}\right)\exp\left({\boldsymbol{Y}}^{T}\tilde{\boldsymbol{\theta}}-\frac{1}{2}\tilde{\boldsymbol{\theta}}^{T}\Lambda\tilde{\boldsymbol{\theta}}\right)\pi(\tilde{\boldsymbol{\theta}}),$
(2.9)

where ${\boldsymbol{Y}}$ is a vector of only $p$ elements, consisting of the
MOPED coefficients:

${\boldsymbol{Y}}={\sf{B}}^{T}{\boldsymbol{X}}={\sf{B}}^{T}({\boldsymbol{x}}-{\boldsymbol{\mu}}_{*}),$
(2.10)

and ${\sf{B}}^{T}\in\mathbb{R}^{p\times n}$ is made up of $p$ rows of MOPED
vectors ${\boldsymbol{b}}_{\alpha}^{T}$, defined by equation (2.7). In matrix
form,

${\sf{B}}^{T}=\Phi^{T}{\sf{C}}^{-1},$ (2.11)

where $\Phi^{T}\in\mathbb{R}^{p\times n}$ is a matrix with rows given by
${\boldsymbol{\mu}}_{,\alpha}^{T}$. The symmetric matrix
$\Lambda\in\mathbb{R}^{p\times p}$ is

$\Lambda={\sf{B}}^{T}{\sf{C}}{\sf{B}}=\Phi^{T}{\sf{C}}^{-1}\Phi.$ (2.12)

Note that Heavens et al. [15] applied a Gram-Schmidt orthogonalisation to the
MOPED vectors, to decorrelate the elements of ${\boldsymbol{Y}}$, but the
information content is the same if this is not applied, and we omit it as it
simplifies the analysis. We have also subtracted ${\boldsymbol{\mu}}_{*}$ from
${\boldsymbol{x}}$ for convenience, and to ensure that a Taylor expansion to
linear order is accurate, following [6].

To normalise the posterior,
$p({\boldsymbol{\theta}}|{\boldsymbol{x}},M)=p({\boldsymbol{x}}|{\boldsymbol{\theta}},M)\pi({\boldsymbol{\theta}})/p({\boldsymbol{x}}|M)$,
we need the Bayesian Evidence $Z({\boldsymbol{x}})=p({\boldsymbol{x}}|M)$,
where we have added explicitly the dependence on the choice of model, $M$. The
Evidence is

$\displaystyle p({\boldsymbol{x}}|M)$
$\displaystyle=\frac{1}{\sqrt{|2\pi{\sf{C}}|}}\exp\left(-\frac{1}{2}{\boldsymbol{X}}^{T}{\sf{C}}^{-1}{\boldsymbol{X}}\right)\int\,\exp\left({\boldsymbol{Y}}^{T}\tilde{\boldsymbol{\theta}}-\frac{1}{2}\tilde{\boldsymbol{\theta}}^{T}\Lambda\tilde{\boldsymbol{\theta}}\right)\pi(\tilde{\boldsymbol{\theta}})\,d\tilde{\boldsymbol{\theta}}.$
(2.13)

We see that the data-only prefactor cancels in the normalized posterior:

$\displaystyle p({\boldsymbol{\theta}}|{\boldsymbol{x}},M)$
$\displaystyle=\frac{\exp\left({\boldsymbol{Y}}^{T}\tilde{\boldsymbol{\theta}}-\frac{1}{2}\tilde{\boldsymbol{\theta}}^{T}\Lambda\tilde{\boldsymbol{\theta}}\right)\pi(\tilde{\boldsymbol{\theta}}|M)}{\int\,\exp\left({\boldsymbol{Y}}^{T}\tilde{\boldsymbol{\theta}}-\frac{1}{2}\tilde{\boldsymbol{\theta}}^{T}\Lambda\tilde{\boldsymbol{\theta}}\right)\pi(\tilde{\boldsymbol{\theta}}|M)\,d\tilde{\boldsymbol{\theta}}}.$
(2.14)

From this it is apparent that, to the extent that the linear expansion is
applicable, we may compute the posterior from the MOPED coefficients alone,
without considering ${\boldsymbol{x}}$ in full — we only need ${\bf Y}$. This
is a very powerful feature of MOPED compression.

Note that for a Gaussian prior density,
$\pi(\tilde{\boldsymbol{\theta}}|M)={\cal{N}}({\bf 0},\Sigma)$ the integrals
may be evaluated analytically (see also [35, 36, 13]) giving

$\displaystyle Z({\boldsymbol{X}})$
$\displaystyle=\frac{1}{\sqrt{|2\pi{\sf{C}}||\Sigma\Lambda+\mathbb{I}_{p}|}}\exp\left(-\frac{1}{2}{\boldsymbol{X}}^{T}{\sf{C}}^{-1}{\boldsymbol{X}}\right)\exp\left[\frac{1}{2}{\boldsymbol{Y}}^{T}\left(\Lambda+\Sigma^{-1}\right)^{-1}{\boldsymbol{Y}}\right].$
(2.15)

####  2.1.1 Demonstration that the Bayes factor depends only on MOPED
coefficients

Before we derive the general expression for the Bayes factor in the next
section, we motivate our approach by demonstrating that, for nested models,
the Bayes factor from the full dataset depends only on the set of MOPED
coefficients defined in the extended model, to linear order in expansion. The
posterior odds between models $M_{0}$ and $M_{1}$ are

$\frac{p(M_{0}|{\boldsymbol{x}})}{p(M_{1}|{\boldsymbol{x}})}=\frac{\pi(M_{0})p({\boldsymbol{x}}|M_{0})}{\pi(M_{1})p({\boldsymbol{x}}|M_{1})}=B_{01}({\boldsymbol{x}})\equiv\frac{Z_{0}({\boldsymbol{x}})}{Z_{1}({\boldsymbol{x}})}$
(2.16)

where $B_{01}$ is the Bayes factor, and the last equality applies if we assume
equal prior model probabilities, $\pi(M_{0})=\pi(M_{1})$. We have defined
$Z_{i}\equiv p({\boldsymbol{x}}|M_{i})$.

Model $M_{0}$, with $q$ free parameters, is said to be ‘nested’ in model
$M_{1}$, with $p>q$ parameters, if the first $q$ parameters of Model 1 are
common to Model 0, and there is a choice of value for the additional
parameters in Model 1 such that the latter reverts to Model 0. In this
setting, we choose a fiducial parameter set within the prior support of Model
0 (and, therefore, a fortiori in Model 1); in this case,
${\boldsymbol{\mu}}_{*}$ is common to both models, and ${\boldsymbol{X}}$ is
the same in both. in this case, the prefactor in the Evidence (equation 2.13),
cancels out in the Bayes factor, for any choice of prior densities
$\pi(\tilde{\boldsymbol{\theta}}_{i}|M_{i})$

$B_{01}({\boldsymbol{x}})\equiv\frac{Z_{0}({\boldsymbol{x}})}{Z_{1}({\boldsymbol{x}})}=\frac{\int\,\exp\left({\boldsymbol{Y}}^{T}\tilde{\boldsymbol{\theta}}_{0}-\frac{1}{2}\tilde{\boldsymbol{\theta}}_{0}^{T}\Lambda_{0}\tilde{\boldsymbol{\theta}}_{0}\right)\pi(\tilde{\boldsymbol{\theta}}_{0}|M_{0})\,d\tilde{\boldsymbol{\theta}}_{0}}{\int\,\exp\left({\boldsymbol{Y}}^{T}\tilde{\boldsymbol{\theta}}_{1}-\frac{1}{2}\tilde{\boldsymbol{\theta}}_{1}^{T}\Lambda_{1}\tilde{\boldsymbol{\theta}}_{1}\right)\pi(\tilde{\boldsymbol{\theta}}_{1}|M_{1})\,d\tilde{\boldsymbol{\theta}}_{1}}.$
(2.17)

Note that for the simpler model the extra parameters are set at
$\tilde{\boldsymbol{\theta}}={\bf 0}$, which we accommodate by defining
$\Lambda_{0}$ as the top $q\times q$ block of $\Lambda_{1}$.

Just like the posterior for the parameters, the Bayes factor for the full
dataset depends on the data only via ${\boldsymbol{Y}}$. This motivates
adopting MOPED compression for model comparison in the general case, with the
potential for a large acceleration of what can be a computationally expensive
task for some applications. We develop this in the next section.

It is worth pointing out that the requirements on the accuracy of the Taylor
expansion are more stringent for model comparison than for parameter
inference, since we are evaluating an integral over the whole parameter space,
well beyond the region where the posterior density is appreciable. For some
applications this could limit the accuracy of the approach, but for typical
cosmological posteriors, their unimodal and relatively simple structure means
that the MOPED-derived posteriors are often very close to those obtained from
the full dataset, even far from the posterior peak. This is apparent even in
the first (astrophysical) example, presented in [15], where the posterior
contours are similar between MOPED and full data set at values smaller than
the peak by a factor $e^{100}$, even though this is not guaranteed by the
method - the compression does far better than one has any right to expect. We
therefore compress the data only once for each model, each at a single
fiducial parameter set, and we illustrate later in a cosmological supernova
study that the Bayes factor computation from the compressed data still agrees
with the full dataset within sampling errors.

###  2.2 Bayesian Evidence and Bayes Factor with MOPED compression for non-
nested models

We now develop the formalism for using MOPED extreme data compression for
computing Bayes factors for non-nested models, with the nested case requiring
only a minimal adjustment. In order for the Bayes factor to give the correct
Bayesian posterior odds ratio (for equal model priors), the data entering in
each of the Bayesian Evidence calculations have to be the same, or else the
$p({\rm data})$ denominator in Bayes theorem will not cancel in the ratio.
Since the MOPED vectors defined within each model are optimal for parameter
inference in their respective models, we construct a compressed dataset from
the union of the MOPED coefficients from both (or, more generally, all)
models. In this way the augmented compressed dataset will be optimal for
parameter inference in both models, since we will not lose information by
adding additional data beyond what is required. Of course we must take care to
include correctly the correlations between all the compressed data. We thus
compare the analysis with the full dataset with the data compressed to

${\boldsymbol{Y}}^{T}=\begin{pmatrix}{\boldsymbol{Y}}_{0}^{T}\\\
{\boldsymbol{Y}}_{1}^{T}\\\
\end{pmatrix}=\begin{pmatrix}{\sf{B}}_{0}^{T}({\boldsymbol{x}}-{\boldsymbol{\mu}}_{*0})\\\
{\sf{B}}_{1}^{T}({\boldsymbol{x}}-{\boldsymbol{\mu}}_{*1})\\\
\end{pmatrix}\in\mathbb{R}^{p+q},$ (2.18)

where ${\boldsymbol{\mu}}_{*0}$ and ${\boldsymbol{\mu}}_{*1}$ are fiducial
means in the two models. In the nested case, these would be chosen to be
equal, and since ${\boldsymbol{Y}}_{0}$ is already contained in
${\boldsymbol{Y}}_{1}$, it is discarded, with ${\boldsymbol{Y}}_{1}$ alone
being retained. Here we assume non-nested models when both
${\boldsymbol{Y}}_{0}$ and ${\boldsymbol{Y}}_{1}$ are present, and specialize
to the nested case later.

The Bayesian Evidence for each model, $Z_{i}({{\boldsymbol{Y}}})$ is obtained
by very similar arguments to the preceding analysis. For Gaussian data
${\boldsymbol{x}}\sim{\cal N}({\boldsymbol{\mu}},{\sf{C}})$,
${\boldsymbol{Y}}\sim{\cal N}(\bar{\boldsymbol{Y}},\Lambda)$ is also Gaussian-
distributed, with mean

$\bar{\boldsymbol{Y}}(\tilde{\boldsymbol{\theta}})=\begin{pmatrix}{\sf{B}}_{0}^{T}({\boldsymbol{\mu}}-{\boldsymbol{\mu}}_{*0})\\\
{\sf{B}}_{1}^{T}({\boldsymbol{\mu}}-{\boldsymbol{\mu}}_{*1})\\\
\end{pmatrix}.$ (2.19)

and covariance
$\Lambda=\Phi^{T}{\sf{C}}^{-1}\Phi\in\mathbb{R}^{(p+q)\times(p+q)}$, where in
block form $\Phi=(\Phi_{0},\Phi_{1})$. The Bayesian Evidence for model $i$ is
then

$Z_{i}({\boldsymbol{Y}})=\frac{1}{\sqrt{|2\pi\Lambda|}}\int\exp\left[-\frac{1}{2}\left({\boldsymbol{Y}}-\bar{\boldsymbol{Y}}\right)^{T}\Lambda^{-1}\left({\boldsymbol{Y}}-\bar{\boldsymbol{Y}}\right)\right]\,\pi(\tilde{\boldsymbol{\theta}}_{i}|M_{i})\,d\tilde{\boldsymbol{\theta}}_{i}.$
(2.20)

With a Taylor expansion for ${\boldsymbol{\mu}}$ as before (equation (2.4)),
$\bar{\boldsymbol{Y}}={\sf{B}}^{T}\Phi\tilde{\boldsymbol{\theta}}\equiv\Lambda\tilde{\boldsymbol{\theta}}$,
we find

$\displaystyle Z_{i}({\boldsymbol{Y}})$
$\displaystyle=\frac{1}{\sqrt{|2\pi\Lambda|}}\exp\left[-\frac{1}{2}({\boldsymbol{Y}}-\bar{\boldsymbol{Y}}_{*i})^{T}\Lambda^{-1}({\boldsymbol{Y}}-\bar{\boldsymbol{Y}}_{*i})\right]\times$
$\displaystyle\int\exp\left[({\boldsymbol{Y}}-\bar{\boldsymbol{Y}}_{*i})^{T}\Lambda^{-1}{\sf{A}}_{i}\tilde{\boldsymbol{\theta}}_{i}-\frac{1}{2}\tilde{\boldsymbol{\theta}}_{i}^{T}{\sf{A}}_{i}^{T}\Lambda^{-1}{\sf{A}}_{i}\tilde{\boldsymbol{\theta}}_{i}\right]\pi(\tilde{\boldsymbol{\theta}}_{i}|M_{i})d\tilde{\boldsymbol{\theta}}_{i},$
(2.21)

where ${\sf{A}}_{i}\equiv\Phi^{T}{\sf{C}}^{-1}\Phi_{i}$ and

$\bar{\boldsymbol{Y}}_{*i}=\begin{pmatrix}{\sf{B}}_{0}^{T}({\boldsymbol{\mu}}_{*i}-{\boldsymbol{\mu}}_{*0})\\\
{\sf{B}}_{1}^{T}({\boldsymbol{\mu}}_{*i}-{\boldsymbol{\mu}}_{*1})\\\
\end{pmatrix},$ (2.22)

which will have zeros in one or other block depending on whether $i=0$ or $1$.

####  2.2.1 Nested models with MOPED compression

There is a significant simplification for nested models, where we dispense
with the top blocks of ${\boldsymbol{Y}}$ and $\Phi^{T}$, and take
${\boldsymbol{\mu}}_{*0}={\boldsymbol{\mu}}_{*1}$. As a consequence
$\bar{\boldsymbol{Y}}_{*i}={\bf 0}$. ${\sf{A}}_{i}$ can be replaced by
$\Lambda$ (which reduces to $\Lambda_{1})$, since it always comes in
combination with $\tilde{\boldsymbol{\theta}}_{i}$, and we extend the
$\tilde{\boldsymbol{\theta}}_{0}$ vector to include zeros from the extended
parameter space. With these simplifications, the Evidence becomes

$\displaystyle Z_{i}({\boldsymbol{Y}})$
$\displaystyle=\frac{1}{\sqrt{|2\pi\Lambda|}}\exp\left(-\frac{1}{2}{\boldsymbol{Y}}^{T}\Lambda^{-1}{\boldsymbol{Y}}\right)\int\exp\left({\boldsymbol{Y}}^{T}\tilde{\boldsymbol{\theta}}_{i}-\frac{1}{2}\tilde{\boldsymbol{\theta}}_{i}^{T}\Lambda\tilde{\boldsymbol{\theta}}_{i}\right)\pi(\tilde{\boldsymbol{\theta}}_{i}|M_{i})d\tilde{\boldsymbol{\theta}}_{i}.$
(2.23)

This is very similar to the Evidence computed from the full data set (equation
2.13), but with a different data-dependent prefactor to the integral.

When we set the prior on the extended model parameters to Dirac delta
functions in Model 0, $\Lambda$ effectively becomes $\Lambda_{0}$, and when we
cancel the common prefactors, we recover equation (2.17), and find that the
Bayes factor for the MOPED compressed data for nested models is exactly the
same as the Bayes factor for the original dataset. In other words, to the
extent that a linear Taylor expansion is accurate, no information is lost in
performing Bayesian model comparison using MOPED coefficients rather than the
full dataset.

####  2.2.2 Gaussian prior

With a gaussian prior on the parameters, the Bayesian Evidence (equation 2.21)
for the different models (nested or otherwise) reduces, after application of
the Woodbury formula and the matrix determinant lemma, to a very simple
expression

$Z_{i}({{\boldsymbol{Y}}})=\frac{1}{\sqrt{|2\pi
Q_{i}|}}\exp\left[-\frac{1}{2}({\boldsymbol{Y}}-\bar{\boldsymbol{Y}}_{*i})^{T}Q_{i}^{-1}({\boldsymbol{Y}}-\bar{\boldsymbol{Y}}_{*i})\right]$
(2.24)

where

$Q_{i}=\Lambda+{\sf{A}}_{i}^{T}\Sigma_{i}{\sf{A}}_{i}.$ (2.25)

###  2.3 Score compression

We can generalise the result further to arbitrary likelihood functions, using
the score compression of Alsing & Wandelt [17], of which MOPED is a special
case. Here they expand the log likelihood as

${\mathcal{L}}={\mathcal{L}}_{*}+\tilde{\boldsymbol{\theta}}^{T}{\boldsymbol{\nabla}}{\mathcal{L}}_{*}-\frac{1}{2}\tilde{\boldsymbol{\theta}}^{T}{\boldsymbol{J}}_{*}\tilde{\boldsymbol{\theta}}$
(2.26)

where
${\boldsymbol{J}}\equiv-{\boldsymbol{\nabla}}{\boldsymbol{\nabla}}^{T}{\mathcal{L}}$,
and $*$ indicates that quantities are evaluated at the fiducial parameters.

${\mathcal{L}}_{*}$ depends on the data, but we see again that the
${\mathcal{L}}_{*}$ (which is a source of variability in the Bayesian
Evidence) cancels in the Bayes factor if the models are nested and
${\boldsymbol{\theta}}_{*}$ lies in the domain of both $M_{0}$ and $M_{1}$.
The only coupling between data and parameters at linear order is through the
score function

${\boldsymbol{t}}\equiv{\boldsymbol{\nabla}}{\mathcal{L}}_{*}.$ (2.27)

${\bf J}$ is strictly constant in the Gaussian case if ${\sf{C}}$ is
independent of ${\boldsymbol{\theta}}$, and is given by the Fisher matrix

$\langle{\boldsymbol{J}}\rangle={\boldsymbol{F}}={\boldsymbol{\nabla}}{\boldsymbol{\mu}}^{T}{\sf{C}}{\boldsymbol{\nabla}}^{T}{\boldsymbol{\mu}}.$
(2.28)

If ${\sf{C}}$ is parameter-dependent, the expression is complicated, and
couples data and parameters (see [17], equation (15)), and replacing
${\boldsymbol{J}}$ by ${\boldsymbol{F}}$ would be an approximation, but
possibly a good one.

We see that the preceding arguments hold in this case as well, that for nested
models, the ${\mathcal{L}}_{*}$ term cancels out in both the posterior and
also the Bayes factor, if the fiducial parameters are common to both (all)
models. As a result, general score compression, using the trick of
concatenating the scores from all models considered, can effectively be used
for model comparison, with the additional approximation (2.28).

##  3 Example Applications

###  3.1 Distribution of Bayes factors

Bayesians may wish to turn away now, as we consider the frequentist properties
of Bayesian Evidence and Bayes factors, as investigated by Jenkins & Peacock
[37] and Joachimi et al. [19]. These papers note that the Bayesian Evidence
and the Bayes factor are noisy statistics, and consider their sampling
variability and its impact on the outcome of model comparison. Note that we
remain of the Bayesian view that the Bayes factor can indicate a preference
for one model over another, but not falsify models in absolute terms, and we
are not seeking frequentist methods to do this, which typically rely on tail
probabilities of $p({\rm data}|{M})$ as opposed to posterior model
probabilities, $p({M}|{\rm data})$, cf. [38, 39, 40].

To the committed Bayesian, the sampling distribution is not relevant, since
the Bayesian view is that there is one set of data on which the inference is
based, and it does not matter what other realisations of the data might have
yielded. Nevertheless, $Z$ is a statistic, and its sampling distribution could
be of some interest. There is a contribution to the variability of $Z$ from
the data-only prefactor in equation (2.13) in the gaussian case (and we come
to a similar conclusion in the more general score compression case: see
section 2.3). The log of the prefactor is distributed as a $\chi^{2}$
distribution with $n$ degrees of freedom with a variance of $2n$. However,
this main source of variability (the prefactor) cancels in nested models when
the Bayes factor is computed, so the quantity on which model comparison
depends is not subject to this source of large variability.

This behaviour is apparent in Fig. B.1 of [19], where the Bayes factor is very
stable even though $Z_{0}$ and $Z_{1}$ are highly variable, which would invite
the naïve interpretation that the sampling variance of the Evidence may hinder
model comparison. However, even if the sampling distribution of $\ln Z_{0}$
and $\ln Z_{1}$ is large in comparison with $\ln B_{01}({\boldsymbol{X}})$,
the Bayes factor favours one model consistently, because $Z_{0}$ and $Z_{1}$
are not independent, as they use the same data realisation.

We illustrate the behaviour of a non-nested case with the following toy non-
linear models:

  * •

$M_{0}:\ x=E\sin(Ft);\quad{\boldsymbol{\theta}}_{0}=(E,F)$

  * •

$M_{1}:\ x=at^{2}+bt;\quad{\boldsymbol{\theta}}_{1}=(a,b)$

where $t$ has $20$ values uniformly spaced between $0$ and $\pi$ (an example
is shown in Fig. 1).

The data are generated with $E=2.5$ and $F=1.0$, and the fiducial values are
$E_{*}=3.0$, $F_{*}=1.2,a_{*}=-1.5,b_{*}=4.7$). The distribution of Bayesian
Evidence values for each model from $4\times 10^{4}$ realisations is shown in
Fig. 2 for the full dataset, and from MOPED-compressed data in Fig. 3. Notice
that the considerable sample variance exceeds the difference of the mean
Evidence values. As suggested by [19], the width of the sampling distribution
of the compressed data Evidence is reduced compared to the full dataset.

Fig. 4 shows that, although the Evidence values differ markedly for the full
and compressed datasets, the Bayes factors are almost identical, with
deviations being due to the linear Taylor expansion not being perfect in the
nonlinear case.

Figure 1: Example dataset, with noisy data drawn from the mean curve of Model
0. Parameters for Model 1 are arbitrary. Figure 2: The distribution of
Bayesian Evidence for $50000$ noise realisations of non-nested toy models
$M_{0}$ and $M_{1}$ with details given in the text, obtained from the full
dataset. The normalisation of the $y$ axis is arbitrary. Figure 3: As Fig. 2,
but for a compressed dataset consisting of the four MOPED coefficients, two
from each model. Figure 4: Comparison of Bayes factor sampling distribution
for the nonlinear, non-nested toy models $M_{0}$ and $M_{1}$ with details
given in the text, from the full dataset and the MOPED-compressed dataset. The
distributions overlay almost perfectly.

Our conclusion from this study is that MOPED compression reduces the variance
of the sampling distribution of the Bayesian Evidence, so it would appear
advantageous to use such a data compression. However, the Evidence values from
an individual noise realisation are highly correlated, as they are based on
the same data, and the log of the Bayes factor has a much smaller variance
than the variance of the individual $\ln Z$ values. We find that the
uncompressed data give the same Bayes factor as the MOPED compressed data
(inasmuch as the neglect of higher-order derivatives of ${\boldsymbol{\mu}}$
is valid). As a result MOPED compression can be used for accurate model
comparison, with in some cases a large increase in speed of computation. We
also note that with a single full likelihood evaluation and one compressed
likelihood evaluation at the fiducial parameters, the full dataset Evidence
can be computed, if it were ever needed, from the compressed Evidence:

$Z({\boldsymbol{X}})=\frac{\sqrt{|2\pi\Lambda|}}{\sqrt{|2\pi{\sf{C}}|}}\frac{\exp(-\frac{1}{2}{\boldsymbol{X}}^{T}{\sf{C}}^{-1}{\boldsymbol{X}})}{\exp(-\frac{1}{2}{\boldsymbol{Y}}^{T}\Lambda^{-1}{\boldsymbol{Y}})}\,Z({\boldsymbol{Y}}).$

###  3.2 Cosmological example — the Pantheon+SH0ES dataset

In this section, we re-analyse the Pantheon+SH0ES data presented in [20]. This
consists of 1701 supernovae, with Cepheid calibration. We use the distance
moduli from the catalogue222https://github.com/PantheonPlusSH0ES (with a
fiducial SNIa magnitude determined from SH0ES 2021 Cepheid host distances),
and the supplied covariance matrix, including systematics, to define a simple
gaussian likelihood function for 1590 supernovae with $z>0.01$. While not
nearly as sophisticated as a Bayesian hierarchical model [e.g. 41], this
simplified setup serves to demonstrate parameter inference and Bayesian
Evidence comparisons for the full and MOPED-compressed datasets. The two
models being compared are:

  * •

a flat universe, parameterised by ${\boldsymbol{\theta}}_{0}=(\Omega_{\rm
m},h)$;

  * •

a curved universe, which allows $\Omega_{\Lambda}$ to vary away from
$1-\Omega_{\rm m}$, with parameters ${\boldsymbol{\theta}}_{1}=(\Omega_{\rm
m},h,\Omega_{\Lambda})$,

where $h$ is the Hubble-Lemaître constant in units of $100\,{\rm km}\,{\rm
s}^{-1}\,{\rm Mpc}^{-1}$ and $\Omega_{\rm m}$ is the matter density parameter,
and $\Omega_{\Lambda}$ the vacuum energy density parameter. Priors on all
varying parameters are taken to be uniform over the range $[0,1]$. We choose a
fiducial model with $\Omega_{\rm m}=0.3$, $h=0.7$ and $\Omega_{\Lambda}=0.7$.
The compression concatenates the MOPED coefficients from both models, but
since they are nested, this is identical to using the MOPED coefficients from
the extended, 3-parameter model, yielding

${\bf Y}=(1304.68,\ 23857.17,\ -1341.94).$ (3.1)

All the cosmological information about $\Omega_{\rm m},h$ and
$\Omega_{\Lambda}$ is contained in just these three numbers. The covariance
matrix of these coefficients is

${\mathrm{Cov}}({{\boldsymbol{Y}}})=\begin{pmatrix}3120.0&35028.9&-2632.9\\\
35028.9&668580.5&-36364.1\\\ -2632.9&-36364.1&2559.2\end{pmatrix}.$ (3.2)

We use the nested sampling code dynesty333DOI 10.5281/zenodo.3348367 [42] with
5000 live points, and no linear approximation is made for the parameter
dependence of the mean distance modulus in the likelihood function. We find
virtually identical posteriors for full and compressed data, shown in Fig. 5.
The natural log of the Bayesian Evidence for the full dataset is $838.88\pm
0.035$ for the flat model, and $837.96\pm 0.037$ for the extended model. For
the MOPED-compressed data, the log Bayesian Evidence values are $-22.35$ and
$-23.29$, with the same errors. Hence the Bayes factor are

$\displaystyle-0.92$ $\displaystyle\pm$ $\displaystyle
0.05;\qquad\mathrm{Full\ dataset\ (from\ 1590\ numbers)}$ $\displaystyle-0.94$
$\displaystyle\pm$ $\displaystyle 0.05;\qquad\mathrm{MOPED\ compression\
(from\ 3\ numbers),}$ (3.3)

with in both cases the flat model being slightly favoured by the data. The
error is an estimate from the values of about $0.036$ for each, since we do
not know the extent to which the errors are correlated.

Fig. 5 and the Bayes factors of eq.(3.3) show the strength of MOPED for both
parameter inference and now for model comparison as well. The latter, along
with the sampling distributions of the Bayes factor in the compressed case
(Fig. 4) are the main new results of this paper.

Figure 5: Marginal posteriors for cosmological parameters $(\Omega_{\rm
m},h,\Omega_{\Lambda})$, using the full dataset of size 1590 (red) and only
the three MOPED-compressed coefficients (blue). The posteriors are almost
identical, resulting in a purple shade in the 1D marginals.

##  4 Discussion and Conclusions

In this paper we have shown how the extreme data compression algorithm MOPED
[15] can be used to compute Bayes factors with essentially no loss of
accuracy, thus extending its power from parameter inference to now include
model comparison. The key to this is to perform MOPED compression in all the
model spaces, and to concatenate the MOPED coefficients into a slightly larger
compressed data space, of dimension equal to the sum of the number of
parameters if the models are not nested. This is still an extreme compression,
if the number of data far exceed the number of parameters, which is typically
the case with cosmological survey data. In the nested case the MOPED
coefficients come from the more complex model. Many of the advantages extend
to other forms of score compression when MOPED is not applicable.

The computation of Bayes factors is generally expensive, requiring integrals
over the parameter space or sophisticated nested sampling techniques, so where
MOPED leads to much faster likelihood evaluations from the massively reduced
dataset, the process may be significantly accelerated.

Finally, we have investigated the sampling distribution of the Bayesian
Evidence and of Bayes factors, for uncompressed and MOPED compressed data. We
found that while extreme data compression reduces the variance of the Evidence
for each individual model, as previously claimed, the sampling distribution
for the Bayes factor remains unaffected by extreme compression of linear
models, and virtually identical for nonlinear models. Thus treating the Bayes
factor as a frequentist statistic, it is just as effective to compute it from
the compressed dataset as the full set. In a cosmological application, we have
shown that, assuming a simple gaussian data model for the Pantheon+SH0ES
dataset [20], the Bayes factor for the flat vs non-flat model can be computed
from only 3 compressed numbers, rather than the full dataset of 1590
supernovae, with no loss of accuracy.

## Acknowledgements

We are grateful to Andrew Jaffe, David van Dyk, Louis Lyons, and Alex
Geringer-Sameth for useful discussions. RT acknowledge co-funding from Next
Generation EU, in the context of the National Recovery and Resilience Plan,
Investment PE1 – Project FAIR “Future Artificial Intelligence Research”. This
resource was co-financed by the Next Generation EU [DM 1555 del 11.10.22]. RT
is partially supported by the Fondazione ICSC, Spoke 3 “Astrophysics and
Cosmos Observations”, Piano Nazionale di Ripresa e Resilienza Project ID
CN00000013 “Italian Research Center on High-Performance Computing, Big Data
and Quantum Computing” funded by MUR Missione 4 Componente 2 Investimento 1.4:
Potenziamento strutture di ricerca e creazione di “campioni nazionali di R&S
(M4C2-19 )” - Next Generation EU (NGEU). This publication was made possible
through the support of an LSSTC Catalyst Fellowship to AM, funded through
Grant 62192 from the John Templeton Foundation to LSST Corporation. The
opinions expressed in this publication are those of the authors and do not
necessarily reflect the views of LSSTC or the John Templeton Foundation.

## Data Availability

No new data are used in this publication, except for those generated randomly
by the processes described fully within.

## References

  * [1] K. Kuijken, C. Heymans, A. Dvornik, H. Hildebrandt, J.T.A. de Jong, A.H. Wright et al., _The fourth data release of the Kilo-Degree Survey: ugri imaging and nine-band optical-IR photometry over 1000 square degrees_ , _A & A_ 625 (2019) A2.
  * [2] C. Doux, B. Jain, D. Zeurcher, J. Lee, X. Fang, R. Rosenfeld et al., _Dark energy survey year 3 results: cosmological constraints from the analysis of cosmic shear in harmonic space_ , _MNRAS_ 515 (2022) 1942.
  * [3] R. Dalal, X. Li, A. Nicola, J. Zuntz, M.A. Strauss, S. Sugiyama et al., _Hyper Suprime-Cam Year 3 Results: Cosmology from Cosmic Shear Power Spectra_ , _arXiv e-prints_ (2023) arXiv:2304.00701 [2304.00701].
  * [4] Planck Collaboration, N. Aghanim, Y. Akrami, F. Arroja, M. Ashdown, J. Aumont et al., _Planck 2018 results. I. Overview and the cosmological legacy of Planck_ , _A & A_ 641 (2020) A1.
  * [5] R. Laureijs, J. Amiaux, S. Arduini, J.L. Auguères, J. Brinchmann, R. Cole et al., _Euclid Definition Study Report_ , _arXiv:1110.3193_ (2011)  [1110.3193].
  * [6] A.F. Heavens, E. Sellentin, D. de Mijolla and A. Vianello, _Massive data compression for parameter-dependent covariance matrices_ , _MNRAS_ 472 (2017) 4244 [1707.06529].
  * [7] T. Charnock, G. Lavaux and B.D. Wandelt, _Automatic physical inference with information maximizing neural networks_ , _PRD_ 97 (2018) 083004 [1802.03537].
  * [8] F. Gerardi, S.M. Feeney and J. Alsing, _Unbiased likelihood-free inference of the Hubble constant from light standard sirens_ , _Phys. Rev. D_ 104 (2021) 083531.
  * [9] T.L. Makinen, T. Charnock, J. Alsing and B.D. Wandelt, _Lossless, scalable implicit likelihood inference for cosmological fields_ , _JCAP_ 2021 (2021) 049 [2107.07405].
  * [10] N. Jeffrey, J. Alsing and F. Lanusse, _Likelihood-free inference with neural compression of DES SV weak lensing map statistics_ , _MNRAS_ 501 (2021) 954 [2009.08459].
  * [11] T.L. Makinen, T. Charnock, P. Lemos, N. Porqueres, A.F. Heavens and B.D. Wandelt, _The Cosmic Graph: Optimal Information Extraction from Large-Scale Structure using Catalogues_ , _The Open Journal of Astrophysics_ 5 (2022) 18 [2207.05202].
  * [12] K. Lin, M. von Wietersheim-Kramsta, B. Joachimi and S. Feeney, _A simulation-based inference pipeline for cosmic shear with the Kilo-Degree Survey_ , _arXiv:2212.04521_ (2022)  [2212.04521].
  * [13] F. Leclercq, _Simulation-based inference of Bayesian hierarchical models while checking for model misspecification_ , _arXiv:2209.11057_ (2022)  [2209.11057].
  * [14] A.F. Heavens, E. Sellentin and A.H. Jaffe, _Extreme data compression while searching for new physics_ , _MNRAS_ 498 (2020) 3440 [2006.06706].
  * [15] A.F. Heavens, R. Jimenez and O. Lahav, _Massive lossless data compression and multiple parameter estimation from galaxy spectra_ , _MNRAS_ 317 (2000) 965 [astro-ph/9911102].
  * [16] M. Tegmark, A.N. Taylor and A.F. Heavens, _Karhunen-Loève Eigenvalue Problems in Cosmology: How Should We Tackle Large Data Sets?_ , _ApJ_ 480 (1997) 22 [astro-ph/9603021].
  * [17] J. Alsing and B. Wandelt, _Generalized massive optimal data compression_ , _MNRAS_ 476 (2018) L60 [1712.00012].
  * [18] G. Jung, D. Karagiannis, M. Liguori, M. Baldi, W.R. Coulton, D. Jamieson et al., _Quijote-PNG: Quasi-maximum Likelihood Estimation of Primordial Non-Gaussianity in the Nonlinear Dark Matter Density Field_ , _ApJ_ 940 (2022) 71 [2206.01624].
  * [19] B. Joachimi, F. Köhlinger, W. Handley and P. Lemos, _When tension is just a fluctuation. How noisy data affect model comparison_ , _A & A_ 647 (2021) L5 [2102.09547].
  * [20] D. Scolnic, D. Brout, A. Carr, A.G. Riess, T.M. Davis, A. Dwomoh et al., _The Pantheon+ Analysis: The Full Data Set and Light-curve Release_ , _ApJ_ 938 (2022) 113 [2112.03863].
  * [21] A. Heavens, B. Panter, R. Jimenez and J. Dunlop, _The star-formation history of the Universe from the stellar populations of nearby galaxies_ , _Nature_ 428 (2004) 625 [astro-ph/0403293].
  * [22] B. Panter, A.F. Heavens and R. Jimenez, _Star formation and metallicity history of the SDSS galaxy survey: unlocking the fossil record_ , _MNRAS_ 343 (2003) 1145 [astro-ph/0211546].
  * [23] B. Panter, R. Jimenez, A.F. Heavens and S. Charlot, _The star formation histories of galaxies in the Sloan Digital Sky Survey_ , _MNRAS_ 378 (2007) 1550 [astro-ph/0608531].
  * [24] B. Panter, R. Jimenez, A.F. Heavens and S. Charlot, _The cosmic evolution of metallicity from the SDSS fossil record_ , _MNRAS_ 391 (2008) 1117 [0804.3091].
  * [25] S. Gupta and A.F. Heavens, _Fast parameter estimation from the cosmic microwave background power spectrum_ , _MNRAS_ 334 (2002) 167 [astro-ph/0108315].
  * [26] A. Zablocki and S. Dodelson, _Extreme data compression for the CMB_ , _Phys. Rev. D_ 93 (2016) 083525 [1512.00072].
  * [27] H. Prince and J. Dunkley, _Data compression in cosmology: A compressed likelihood for Planck data_ , _Phys. Rev. D_ 100 (2019) 083502 [1909.05869].
  * [28] Y. Lai, C. Howlett and T.M. Davis, _Faster cosmological analysis with power spectrum without simulations_ , _arXiv e-prints_ (2023) arXiv:2306.00388 [2306.00388].
  * [29] D. Gualdi, H. Gil-Marín, R.L. Schuhmann, M. Manera, B. Joachimi and O. Lahav, _Enhancing BOSS bispectrum cosmological constraints with maximal compression_ , _MNRAS_ 484 (2019) 3713 [1806.02853].
  * [30] P. Graff, M.P. Hobson and A. Lasenby, _An investigation into the Multiple Optimised Parameter Estimation and Data compression algorithm_ , _MNRAS_ 413 (2011) L66 [1010.5907].
  * [31] L. Thiele, G.A. Marques, J. Liu and M. Shirasaki, _Cosmological constraints from HSC Y1 lensing convergence PDF_ , _arXiv e-prints_ (2023) arXiv:2304.05928 [2304.05928].
  * [32] R. Ruggeri and C. Blake, _Compressing combined probes: redshift weights for joint lensing and clustering analyses_ , _MNRAS_ 498 (2020) 2948 [2006.11297].
  * [33] L. Wenzl, C. Doux, C. Heinrich, R. Bean, B. Jain, O. Doré et al., _Cosmology with the Roman Space Telescope - Synergies with CMB lensing_ , _MNRAS_ 512 (2022) 5311 [2112.07681].
  * [34] P. Protopapas, R. Jimenez and C. Alcock, _Fast identification of transits from light-curves_ , _MNRAS_ 362 (2005) 460 [astro-ph/0502301].
  * [35] R. Trotta, _Bayes in the sky: Bayesian inference and model selection in cosmology_ , _Contemporary Physics_ 49 (2008) 71 [0803.4089].
  * [36] F. Leclercq, W. Enzi, J. Jasche and A. Heavens, _Primordial power spectrum and cosmology from black-box galaxy surveys_ , _MNRAS_ 490 (2019) 4237 [1902.10149].
  * [37] C.R. Jenkins and J.A. Peacock, _The power of bayesian evidence in astronomy_ , _MNRAS_ 413 (2011) 2895.
  * [38] R.L. Wasserstein and N.A. Lazar, _The asa statement on p-values: Context, process, and purpose_ , _The American Statistician_ 70 (2016) 129.
  * [39] H. Koo, R.E. Keeley, A. Shafieloo and B. L’Huillier, _Bayesian vs frequentist: comparing Bayesian model selection with a frequentist approach using the iterative smoothing method_ , _JCAP_ 2022 (2022) 047 [2110.10977].
  * [40] R.E. Keeley and A. Shafieloo, _On the distribution of Bayesian evidence_ , _MNRAS_ 515 (2022) 293 [2111.04231].
  * [41] H. Shariff, X. Jiao, R. Trotta and D.A. van Dyk, _BAHAMAS: New Analysis of Type Ia Supernovae Reveals Inconsistencies with Standard Cosmology_ , _ApJ_ 827 (2016) 1 [1510.05954].
  * [42] J.S. Speagle, _DYNESTY: a dynamic nested sampling package for estimating Bayesian posteriors and evidences_ , _MNRAS_ 493 (2020) 3132 [1904.02180].