# Neural networks and standard cosmography with newly calibrated high redshift
GRB observations

Celia Escamilla-Rivera  0000-0002-8929-250X Maryi Carvajal  Cristian Zamora
and Martin Hendry  0000-0001-8322-5405

###### Abstract

Gamma-ray bursts (GRBs) detected at high redshift can be used to trace the
cosmic expansion history. However, the calibration of their luminosity
distances is not an easy task in comparison to Type Ia Supernovae (SNeIa). To
calibrate these data, correlations between their luminosity and other observed
properties of GRBs need to be identified, and we must consider the validity of
our assumptions about these correlations over their entire observed redshift
range. In this work, we propose a new method to calibrate GRBs as cosmological
distance indicators using SNeIa observations with a completely model-
independent deep learning architecture. An overview of this machine learning
technique was developed in [1] to study the evolution of dark energy models at
high redshift. The aim of the method developed in this work is to combine two
networks: a Recurrent Neural Network (RNN) and a Bayesian Neural Network
(BNN). Using this computational approach, denoted RNN+BNN, we extend the
network’s efficacy by adding the computation of covariance matrices to the
Bayesian process. Once this is done, the SNeIa distance-redshift relation can
be tested on the full GRB sample and therefore used to implement a
cosmographic reconstruction of the distance-redshift relation in different
regimes. Thus, our newly-trained neural network is used to constrain the
parameters describing the kinematical state of the Universe via a cosmographic
approach at high redshifts (up to $z\approx 10$), wherein we require a very
minimal set of assumptions that do not rely on dynamical equations for any
specific theory of gravity.

##  1 Introduction

The accelerated expansion of the Universe can be identified using the
luminosity of Type Ia Supernovae (SNeIa), comprehensive measurements of which
have been compiled over the last few years [2, 3]. However, due to their lack
of sufficiently luminous observable at higher redshift, combined with
constraints on the redshifts at which they form, most of these SNeIa are
observed at $z\approx 1$ – which further restricts our ability to distinguish
between different cosmic expansion histories in a model-independent way. In
this context, the SNeIa Pantheon compilation [4] is the only available sample
that reports an observed SNIa at high redshift ($z=2.3$), with only six
further SNeIa observed at $z\approx 1.4$. Hence, SNeIa will not be sensitive
to differences between the luminosity distance redshift relation at higher
redshifts in different cosmological models, and it is therefore important to
identify cosmological probes that can be extended to higher redshifts.
Furthermore, the robustness of the SNeIa calibration has also been questioned
due the tension on the local value of the Hubble constant, $H_{0}$, with calls
for a revision of the Cepheid period-luminosity relationships [5] alongside a
deep study of the latest Gaia Data Releases [6] and other parallax
measurements from the Hubble Space Telescope (HST).

On the other hand, Gamma-ray bursts (GRBs) are the most energetic cosmic
explosions observed to date and are detectable up to very high redshift [7].
Since their characteristics include a narrow range of intrinsic luminosities,
GRBs can be viewed as standard candles – and hence calibrated to trace the
Hubble expansion at redshifts that extend well beyond those probed by SNeIa.
However, calibrating the luminosity distances of GRBs is not a easy task due
there being several different observables and physical parameters associated
with their afterglow emission. Among the calibrating methods that have been
employed in the literature are correlations between: the GRB spectrum lag and
isotropic peak luminosity ($\tau_{\text{lag}}-L$ relation) [8]; the peak
energy of the $\nu F_{\nu}$ spectrum and the isotropic equivalent energy
($E_{p}-E_{\text{iso}}$ relation) [9]; the time variability and isotropic peak
luminosity ($V-L$ relation) [10], the peak energy of the $\nu F_{\nu}$
spectrum and the isotropic peak luminosity ($E_{\rm p}-L$ relation) [11]; the
minimum rise time of light curve and isotropic peak luminosity ($\tau_{\rm
RT}-L$ relation) [12] and the peak energy and collimation-corrected energy
$(E_{\rm p}-E_{\gamma})$ [13]. All of these correlations depend on a specific
cosmological model that can lead to circularity issues when using the
calibrated GRBs to constrain cosmological parameters. However, several model-
independent methods have been proposed to avoid this issue, e.g. by
calibrating the GRB relations first using a low-redshift subsample – where
there exists a corresponding subsample of SNeIa the luminosity distances of
which can be fitted in a model-independent way – and then extrapolating the
correlations to high redshift.

Following these ideas, in this work we set out to use GRBs to understand the
physics of the cosmic acceleration based on a model-independent description.
Whichever path we take, we seek to avoid the explicit assumption of a prior
cosmological model since this may lead to a biased reconstruction of the
dynamics – affecting directly our inference of the cosmological parameters as
a result of degeneracy problems or overfitting convergence using surveys.
These problems are related to the cosmological tensions, like those of
$\sigma_{8}$ and $H_{0}$ [15, 14].

An optimal strategy towards a model-independent panorama requires only the
postulation a priori of homogeneity and isotropy, where the dynamics are
independent of the energy densities. One path, then, is the so-called
cosmography approach, which has been shown to be a good method for
reconstructing the kinematic evolution of the cosmic acceleration without
assuming a cosmological model (see [17, 16, 19], and references therein). To
combine this approach and work with homogeneous data samples at low/high
redshift, in [1] a new computational tool based on machine learning (ML) for
supernovae, called Recurrent-Bayesian Neural Networks (RNN+BNN), was proposed.
A deep learning architecture inside this ML can produce a trained homogeneous
sample, e.g of luminosity distances. Here, therefore, we extend that previous
work and present a methodology to train supernovae and GRB data up to high-$z$
in order to test the viability of this cosmographic approach to reconstructing
the distance-redshift relation. While this cosmographic approach has been
previously studied from several points of view, setting limits on the standard
cosmography using Taylor series [20], Padé and Chebyshev polynomials [21], all
of them have inherent convergence problems at higher redshifts – and hence are
vulnerable to bias and systematic errors. Also, this approach for late time
observations has been considered in extended theories of gravity [22] and
corrected versions of such theories as studied in [23], where there are
significant consequences for the reconstructions beyond the third order
derivative in $H(z)$.

One interesting proposal to solve these issues was presented in [24], where a
generic EoS depending solely on the form of a dark energy-like term and its
derivative was found. This approach can be adapted to any polynomial, to test
the viability of cosmography since this result can be tested with a supernovae
sample trained via a RNN+BNN network that solves problems with over-fitting at
lower redshifts and increases the density of data in the $z$-region of
interest. Furthermore, an as extension, in this work we design a new neural
network to add a new GRB calibrated dataset, which can therefore help to
compare the cosmographic $\mu(z)$ relation inferred from supernovae
observations (assuming flatness) with the cosmographic $\mu(z)$ relation
inferred from supernovae and GRBs, extending to higher forecast redshifts. In
the near future, high-redshift spectroscopic surveys [25, 26] and/or
gravitational-wave missions and projects [27, 28] will provide accurate data
in the redshift range of $2<z<5$. Meanwhile, our trained RNN+BNN method will
allow us to acquire a higher density of precise distance modulus data, for
both SNeIa and GRBs, over the wide redshift range of $0.01<z<10$.

This paper is organised as follows: In Sec.2 we present the theory behind the
cosmographic approach in order to write a general expression for the
luminosity distance in terms of the cosmographic parameters. In Sec.3 we
describe the current cosmological observations: SNeIa and GRB samples, which
will be used in the RNN+BNN methodology. Also we display the GRB values for
the $\mu(z)$ function calibrated with the current Pantheon supernovae sample.
In Sec.4 we explain the basics behind our RNN+BNN method to train the
cosmological data at hand, in order to the constrain the cosmographic
parameters that appear in the definition of the luminosity distance. We
consider an extension of the Bayesian architecture to solve problems of e.g
overfitting. Also, the method to include cosmography is detailed. In Sec.5 we
present the results for the cosmographic analysis at high-$z$, and their
behaviour according to different $H_{0}$ values. A Bayesian Information
Criteria approach is then introduced to test the viability of the cosmography
obtained. Finally in Sec.6 we summarise and discuss our main results and
conclusions.

##  2 Cosmography background

According to the standard Cosmography paradigm, the only ingredient that we
need to take into account a priori in such an approach is a FLRW space-time
derived from kinematical requirements

$ds^{2}=-c^{2}dt^{2}+a^{2}(t)\left[\frac{dr^{2}}{1-kr^{2}}+r^{2}d\Omega^{2}\right],$
(2.1)

where $c$ is the speed of light, $a(t)$ is the scale factor and $k$ is a
curvature constant. With this metric at hand it is possible to write the
luminosity distance $d_{L}$ as a series expansion in the redshift parameter
$z$, where the coefficients of the expansion are functions of $a(t)$
multiplying its higher order derivatives. With this expression, we can study
the luminosity distance at high redshift in order to infer and constrain the
cosmographic parameters that best describe the model for our data.

At this point, the use of supernovae data (up to $z\approx 2$) is convenient;
moreover the density of the data in luminosity distance can be enhanced if we
consider also including GRBs. In particular, data samples coming from
observations of both supernovae and GRBs can be used to fit directly a
cosmographic expression for the luminosity distance $d_{L}$.

To construct this function, $d_{L}$ can be defined from the relation between
the apparent luminosity, or flux, $l=L/4\pi d^{2}_{L}$ of an object, where $L$
is its absolute luminosity. In a FLRW space-time this can be described for
small distances by the expression

$\displaystyle d_{L}(z)$ $\displaystyle=$
$\displaystyle\frac{c}{H_{0}}[z+\frac{1}{2}(1-q_{0})z^{2}-\frac{1}{6}(1-q_{0}-3q_{0}^{2}+j_{0}+\frac{kc^{2}}{H_{0}^{2}a^{2}(t_{0})})z^{3}+$
(2.2)
$\displaystyle+\frac{1}{24}[2-2q_{0}-15q_{0}^{2}-15q_{0}^{3}+5j_{0}+10q_{0}j_{0}+s_{0}+\frac{2kc^{2}(1+3q_{0})}{H_{0}^{2}a^{2}(t_{0})}]z^{4}+\ldots]\,,$

where the cosmographic parameters are defined as

$\displaystyle H_{0}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a(t)}\frac{da(t)}{dt}|_{t=t_{0}}\equiv\frac{\dot{a}(t)}{a(t)}|_{t=t_{0}}~{},$
(2.3) $\displaystyle q_{0}$ $\displaystyle\equiv$
$\displaystyle-\frac{1}{H^{2}}\frac{1}{a(t)}\frac{d^{2}a(t)}{dt^{2}}|_{t=t_{0}}\equiv-\frac{1}{H^{2}}\frac{\ddot{a}(t)}{a(t)}|_{t=t_{0}}~{},$
(2.4) $\displaystyle j_{0}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{H^{3}}\frac{1}{a(t)}\frac{d^{3}a(t)}{dt^{3}}|_{t=t_{0}}\equiv\frac{1}{H^{3}}\frac{a^{(3)}(t)}{a(t)}|_{t=t_{0}}~{},$
(2.5) $\displaystyle s_{0}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{H^{4}}\frac{1}{a(t)}\frac{d^{4}a(t)}{dt^{4}}|_{t=t_{0}}\equiv\frac{1}{H^{4}}\frac{a^{(4)}(t)}{a(t)}|_{t=t_{0}}~{}.$
(2.6)

Depending on the degree of the degree of Taylor series adopted, the
cosmographic quantities derived may produce ambiguous or degenerate results.
Moreover, this may lead to a misunderstanding in the definition of the proper
distance – although (as we discuss below) we may use the properties of the
data themselves to limit the impact of this ambiguity through appropriate
truncation of the Taylor series. In the rest of this work, we will refer to
luminosity distance as the most direct choice in the measure of distance for
SNeIa and GRBs.

As a first attempt to fit the luminosity distance redshift relation using data
samples, and to address the ill-conditioned behaviour at high-$z$ via
appropriate convergence and truncation of the series, it is useful to recast
$d_{L}$ as a function of the transformed variable $y=z/(1+z)$. In this manner,
we can map $z\in(0,\infty)$ into the range $y\in(0,1)$, and retrieve an
optimal behaviour for the series at any distance (i.e. at any redshift). The
consideration of this new variable $y$ does not change the definition of the
cosmographic parameters but does change the equation expressing the luminosity
distance as a series expansion in terms of $y$. The luminosity distance, at
fourth order in $y$, can be written in terms of the cosmographic parameters as

$\displaystyle d_{L}(y)=$
$\displaystyle\frac{c}{H_{0}}\left\\{y-\frac{1}{2}(q_{0}-3)y^{2}+\frac{1}{6}\left[12-5q_{0}+3q^{2}_{0}-(j_{0}+\Omega_{0})\right]y^{3}+\frac{1}{24}\left[60-7j_{0}-\right.\right.$
(2.7)
$\displaystyle\left.\left.-10\Omega_{0}-32q_{0}+10q_{0}j_{0}+6q_{0}\Omega_{0}+21q^{2}_{0}-15q^{3}_{0}+s_{0}\right]y^{4}+\mathcal{O}(y^{5})\right\\}~{},$

where $\Omega_{0}=1+kc^{2}/H_{0}^{2}a^{2}(t_{0})$, is the total energy
density. Using the latter expression, we can write the luminosity distance
logarithmic Hubble relation as [29]

$\displaystyle\ln{\left[\frac{d_{L}(y)}{y~{}}\right]}\rm Mpc^{-1}=$
$\displaystyle\frac{\ln{10}}{5}\left[\mu(y)-25\right]-\ln{y}=\ln{\left[\frac{c}{H_{0}}\right]}-\frac{1}{2}(q_{0}-3)y+\frac{1}{24}\left[21-4(j_{0}+\Omega_{0})+\right.$
(2.8)
$\displaystyle\left.+q_{0}(9q_{0}-2)\right]y^{2}+\frac{1}{24}\left[15+4\Omega_{0}(q_{0}-1)+j_{0}(8q_{0}-1)-5q_{0}+2q^{2}_{0}\right.$
$\displaystyle\left.-10q^{3}_{0}+s_{0}\right]y^{3}+\mathcal{O}(y^{4})~{},$

and the corresponding expression for the distance modulus is given by

$\displaystyle\mu(y)=$ $\displaystyle
25+\frac{5}{\ln{10}}\left\\{\ln{\left[\frac{c}{H_{0}}\right]}+\ln{y}-\frac{1}{2}(q_{0}-3)y+\frac{1}{24}\left[21-4(j_{0}+\Omega_{0})+q_{0}(9q_{0}-2)\right]y^{2}+\right.$
(2.9)
$\displaystyle+\left.\frac{1}{24}\left[15+4\Omega_{0}(q_{0}-1)+j_{0}(8q_{0}-1)-5q_{0}+2q^{2}_{0}-10q^{3}_{0}+s_{0}\right]y^{3}+\mathcal{O}(y^{4})\right\\}~{}.$

The higher the order in the expansion of $d_{L}$, the better the fit to the
data we can expect since there will be more free parameters – but of course
this comes at the cost of greater degeneracy in these parameters. However, for
a given data sample there will be an upper bound on the order of the series
which gives a statistically significant fit to those data. With this bound,
and our RNN+BNN approach, we will be able to objectively define an optimal
truncation to the expansion without having to simply cut the series by hand.

##  3 The data

###  3.1 Pantheon Type Ia Supernovae

In this paper we select the recent Type Ia supernovae (SNeIa) compilation
known as the Pantheon sample [4]. This sample consists of 1048 SNeIa
compressed in 40 redshift bins. It is the largest spectroscopic SNeIa sample
compiled to date.

As already noted, SNeIa can provide estimates of the distance modulus, $\mu$,
the value of which is related to the luminosity distance $d_{L}$ as follows

$\mu(z)=5\log{\left[\frac{d_{L}(z)}{1\text{Mpc}}\right]}+25,$ (3.1)

where the luminosity distance is given in Mpc. In this distance modulus
expression we should include the nuisance parameter, $\mathcal{M}$, as an
unknown offset of the supernovae absolute magnitude (and including other
possible systematics), which can also be degenerate with the value of $H_{0}$.
As is standard, we assume spatial flatness and suppose that $d_{L}$ can be
related to the comoving distance $D$ using
$d_{L}(z)=\frac{c}{H_{0}}(1+z)D(z),$ where $c$ is the speed of light. We
obtain

$D(z)=\frac{H_{0}}{c}(1+z)^{-1}\,10^{0.2\mu(z)-5}.$ (3.2)

The normalised Hubble function $H(z)/H_{0}$ is derived from the inverse of the
derivative of $D(z)$ with respect to $z$, i.e.

$D(z)=\int^{z}_{0}\frac{H_{0}\,d\tilde{z}}{H(\tilde{z})},$ (3.3)

where $H_{0}$ is the value of the Hubble parameter we consider as a prior to
normalise $D(z)$. For our sample, we calibrated the data by using a value
obtained from late universe measurements. corresponding to $H_{0}=73.24\pm
1.4\,\,\mbox{km s}^{-1}\mbox{Mpc}^{-1},$ from SH0ES + H0LiCOW, with the
corresponding value for the nuisance parameter ${\cal M}=-19.24\pm 0.07$ [30].

###  3.2 Updated GRB dataset

After the reconstruction/prior calibration of $d_{L}$ from SNeIa, we can use
them to calibrate the luminosity correlations of GRBs. These correlations can
be written by considering a general exponential form $R=AQ^{b}$, which can be
re-expressed in linear form as $y=a+bx$, with $y\equiv\log R$, $x\equiv\log Q$
and $a=\log A$. Also, we can write the six luminosity correlations measured in
the comoving frame as follows

$\displaystyle\log\left(\frac{E_{\text{iso}}}{\text{erg}}\right)$
$\displaystyle=$ $\displaystyle a_{1}+b_{1}\log\left(\frac{E_{{\rm
p},i}}{300\,\text{keV}}\right),$ (3.4)
$\displaystyle\log\left(\frac{L}{\text{erg}\,\text{s}^{-1}}\right)$
$\displaystyle=$ $\displaystyle a_{2}+b_{2}\log\left(\frac{\tau_{{\rm
RT},i}}{0.1\,\text{s}}\right),$ (3.5)
$\displaystyle\log\left(\frac{E_{\gamma}}{\text{erg}}\right)$ $\displaystyle=$
$\displaystyle a_{3}+b_{3}\log\left(\frac{E_{{\rm
p},i}}{300\,\text{keV}}\right),$ (3.6)
$\displaystyle\log\left(\frac{L}{\text{erg}\,\text{s}^{-1}}\right)$
$\displaystyle=$ $\displaystyle a_{4}+b_{4}\log\left(\frac{E_{{\rm
p},i}}{300\,\text{keV}}\right),$ (3.7)
$\displaystyle\log\left(\frac{L}{\text{erg}\,\text{s}^{-1}}\right)$
$\displaystyle=$ $\displaystyle
a_{5}+b_{5}\log\left(\frac{V_{i}}{0.02}\right),$ (3.8)
$\displaystyle\log\left(\frac{L}{\text{erg}\,\text{s}^{-1}}\right)$
$\displaystyle=$ $\displaystyle a_{6}+b_{6}\log\left(\frac{\tau_{\rm
lag},i}{0.1\,\text{s}}\right),$ (3.9)

with $V_{i}=V(1+z)$, $E_{{\rm p},i}=E_{\rm p}(1+z)$, $\tau_{{\rm
lag},i}=\tau_{\rm lag}(1+z)^{-1}$ and $\tau_{{\rm RT},i}=\tau_{\rm
RT}(1+z)^{-1}$, each of them obtained from the observations of GRB spectra.

To calibrate these six relations with our Pantheon sample, we consider that
GRBs radiate isotropically by computing their bolometric peak flux
$P_{\text{bolo}}=L/4\pi d^{2}_{L}$, where the uncertainty on $L$ propagates
from the uncertainties on $P_{\text{bolo}}$ and $d_{L}$. In the same manner,
we can compute the isotropic equivalent energy $E_{\text{iso}}$ and the
collimation-correlated energy $E_{\gamma}$, respectively as:

$\displaystyle E_{\text{iso}}=4\pi d^{2}_{L}\,S_{\text{bolo}}(1+z)^{-1},$
(3.10) $\displaystyle E_{\gamma}\equiv
E_{\text{iso}}(1-\cos\theta_{\text{jet}}),$ (3.11)

where $\theta_{\text{jet}}$ is the jet opening angle and $S_{\text{bolo}}$ is
the bolometric fluence of the GRBs. The GRB sample [32] used in this work
consists of 116 long GRBs in the redshift range $0.17\leq z\leq 8.2$. With
this sample we calibrated the luminosity correlations by maximizing the
likelihood [33]

$\mathcal{L}(\sigma_{\text{int}},a,b)\propto\prod_{i}\frac{1}{\sqrt{\sigma^{2}_{\text{int}}+\sigma^{2}_{yi}+b^{2}\sigma_{xi}}}\times\text{exp}\left[-\frac{(y_{i}-a-bx_{i})^{2}}{2(\sigma^{2}_{\text{int}}+\sigma^{2}_{yi}+b^{2}\sigma^{2}_{xi})}\right].$
(3.12)

The best-fitting parameters and their
uncertainties,$(a,b,\sigma_{\text{int}})$, obtained in this way are reported
in Table 1. With these fits we tabulated our full GRB sample, as reported in
Tables 2 and 3.

Correlation | $a$ | $b$ | $\sigma_{\text{int}}$
---|---|---|---
$E_{p}-E_{\text{iso}}$ | $52.778\pm 0.058$ | $1.546\pm 0.129$ | $0.461\pm 0.051$
$\tau_{RT}-L$ | $52.766\pm 0.078$ | $-1.250\pm 0.132$ | $0.420\pm 0.055$
$E_{p}-E_{\gamma}$ | $51.655\pm 0.065$ | $1.455\pm 0.156$ | $0.136\pm 0.010$
$E_{p}-L$ | $52.167\pm 0.059$ | $1.450\pm 0.138$ | $0.531\pm 0.049$
$V-L$ | $51.785\pm 0.140$ | $0.556\pm 0.131$ | $0.747\pm 0.069$
$\tau_{lag}-L$ | $52.366\pm 0.071$ | $-0.781\pm 0.116$ | $0.452\pm 0.058$
Table 1: Best-fit values, and their uncertainties, for the parameters
describing the linear correlations between the GRB observational data and
intrinsic properties. In each case we performed an MCMC analysis to calculate
the posterior probability density function of the parameter space with a flat
prior on all free parameters and a fixed constraint of
$\sigma_{\text{int}}>0$.

$z$ | $\mu$ | $\sigma_{\mu}$ | $z$ | $\mu$ | $\sigma_{\mu}$ | $z$ | $\mu$ | $\sigma_{\mu}$
---|---|---|---|---|---|---|---|---
0.03351 | 35.1195 | 0.933638 | 1.44 | 43.2858 | 0.95932 | 2.199 | 46.4216 | 1.06071
0.125 | 38.6856 | 1.84671 | 1.4436 | 43.8132 | 1.01313 | 2.2 | 44.9272 | 0.931098
0.17 | 38.4667 | 1.03564 | 1.46 | 44.1083 | 1.11878 | 2.20 | 45 46.8579 | 0.931104
0.25 | 39.9357 | 1.38443 | 1.48 | 43.9951 | 1.0662 | 2.22 | 45.2389 | 1.20041
0.3399 | 40.6749 | 0.95616 | 1.48 | 43.5582 | 0.929143 | 2.26 | 43.8782 | 0.929903
0.36 | 42.9717 | 1.07033 | 1.489 | 45.0134 | 1.043 | 2.27 | 45.029 | 0.935913
0.41 | 40.8059 | 1.1032 | 1.52 | 42.871 | 0.978916 | 2.296 | 45.502 | 1.15448
0.414 | 42.8313 | 1.21619 | 1.547 | 45.9243 | 0.946461 | 2.3 | 46.2112 | 1.2461
0.434 | 41.5047 | 0.985045 | 1.547 | 44.0644 | 0.973486 | 2.346 | 46.8641 | 1.16259
0.45 | 41.8452 | 1.0036 | 1.563 | 42.781 | 1.69995 | 2.433 | 46.867 | 1.09682
0.49 | 39.7308 | 0.946699 | 1.567 | 43.8473 | 0.927328 | 2.452 | 47.364 | 1.13965
0.5295 | 42.6995 | 1.03825 | 1.6 | 44.1755 | 1.0596 | 2.486 | 44.8287 | 0.944518
0.54 | 39.9036 | 0.925708 | 1.604 | 46.5316 | 0.962389 | 2.488 | 45.181 | 0.942494
0.542 | 41.4047 | 0.975989 | 1.608 | 46.8731 | 1.03718 | 2.512 | 45.5879 | 0.934876
0.543 | 42.2564 | 1.32376 | 1.613 | 44.6927 | 0.948442 | 2.58 | 44.4939 | 0.943808
0.544 | 42.0376 | 0.927353 | 1.619 | 44.3532 | 0.946002 | 2.591 | 46.1254 | 0.958853
0.55 | 42.5063 | 0.951189 | 1.64 | 45.0375 | 0.958753 | 2.615 | 45.9102 | 1.00661
0.606 | 41.6888 | 0.940851 | 1.6919 | 43.2837 | 1.03666 | 2.65 | 45.6525 | 1.00258
0.618 | 43.7538 | 1.07668 | 1.71 | 45.675 | 0.960875 | 2.671 | 45.3823 | 0.943181
0.6528 | 41.8718 | 1.00894 | 1.727 | 43.2599 | 0.973643 | 2.69 | 46.1013 | 1.16024
0.677 | 43.7533 | 0.983363 | 1.728 | 44.3889 | 0.946023 | 2.752 | 45.5087 | 1.07352
0.689 | 43.868 | 0.93225 | 1.728 | 44.3889 | 0.946023 | 2.77 | 46.0964 | 0.9302
0.69 | 42.2854 | 0.928094 | 1.77 | 42.9739 | 0.944986 | 2.821 | 46.6406 | 0.937652
0.695 | 43.1984 | 1.11582 | 1.798 | 44.5549 | 0.94189 | 2.83 | 46.6737 | 0.930181
0.706 | 41.2777 | 0.945272 | 1.8 | 45.3824 | 0.948821 | 2.893 | 47.598 | 1.08191
0.716 | 40.4242 | 0.95484 | 1.822 | 44.7661 | 0.925467 | 2.8983 | 45.3038 | 1.02857
0.716 | 41.2868 | 1.01881 | 1.858 | 45.1695 | 1.23345 | 3. | 46.1498 | 1.10789
0.736 | 44.393 | 0.9305 | 1.858 | 45.1695 | 1.23345 | 3.036 | 46.0503 | 0.959184
0.78 | 40.8528 | 0.974974 | 1.8836 | 46.2739 | 1.28842 | 3.038 | 44.8963 | 1.19331
0.8 | 42.3122 | 1.13123 | 1.9 | 45.7854 | 1.13272 | 3.075 | 46.7497 | 1.03846
0.8049 | 44.9536 | 0.94066 | 1.9229 | 44.6651 | 1.11463 | 3.2 | 45.7872 | 1.11734
0.82 | 42.816 | 0.942884 | 1.95 | 46.5089 | 1.08858 | 3.22 | 45.5887 | 1.00027
0.835 | 44.1715 | 1.08963 | 1.9685 | 44.6084 | 0.999689 | 3.35 | 47.5973 | 0.954116
0.842 | 42.7009 | 1.09192 | 1.98 | 44.5553 | 1.02204 | 3.36 | 47.7168 | 1.01118
0.846 | 43.8624 | 0.92956 | 2.05 | 46.8257 | 1.31834 | 3.37 | 47.4274 | 1.22111
0.859 | 43.1051 | 1.13207 | 2.07 | 43.9293 | 0.936188 | 3.42 | 47.0429 | 0.996013
0.8969 | 44.5543 | 0.999868 | 2.0858 | 44.6893 | 1.56663 | 3.424 | 47.0792 | 1.05247
0.937 | 42.9071 | 0.931116 | 2.088 | 45.8633 | 0.983381 | 3.5 | 46.7831 | 1.23155
0.947 | 42.1526 | 1.03421 | 2.106 | 46.3176 | 0.965919 | 3.5 | 45.2884 | 0.974999
0.958 | 42.9866 | 1.00157 | 2.1062 | 43.7433 | 0.925125 | 3.5328 | 46.978 | 1.1415
0.966 | 44.2503 | 0.957411 | 2.14 | 44.7932 | 0.957566 | 3.57 | 45.8457 | 0.986639
0.971 | 42.7535 | 1.19374 | 2.145 | 47.4073 | 0.984835 | 3.6 | 45.6286 | 1.09246

Table 2: Compilation of estimated GRB luminosity distance moduli, and their
errors, for redshifts $z<3.6$

.

$z$ | $\mu$ | $\sigma_{\mu}$
---|---|---
3.758 | 47.5053 | 0.963238
3.796 | 46.1613 | 1.12782
3.91 | 46.6447 | 1.11229
3.93 | 47.0881 | 1.17098
4.0559 | 48.3284 | 1.27013
4.109 | 46.899 | 1.19876
4.1745 | 47.4784 | 1.00188
4.35 | 47.145 | 0.960048
4.5 | 46.1465 | 1.22279
4.6497 | 49.1939 | 1.1254
5.46 | 47.1793 | 1.06039
6.295 | 49.5673 | 1.13213
6.695 | 49.788 | 1.31553
8.1 | 49.2684 | 1.0345
9.3 | 50.0158 | 0.959439
Table 3: Compilation of estimated GRB luminosity distance moduli, and their
errors, calibrated for redshifts $z>3.6$

##  4 RNN+BNN method for Cosmography

To combine our samples at low and high redshift, in [1] a new computational
tool was proposed for supernovae, based on machine learning (ML), called
Recurrent-Bayesian Neural Networks (RNN+BNN). In this work we extend this
approach to introduce our GRB sample, and apply our deep learning architecture
to obtain a trained homogeneous sample comprising both SNeIa and GRBs. The
method consists of a neural network (NN) implemented with a non-linear
regression process using the samples described in Sec.3. We adopt a real
target to train the NN for each data point. If the first point trained is far
away from the real one, then the algorithm penalizes this point and continues
the process until it reaches a true value. When the training is carried out
for the full sample, then the algorithm proceeds to minimize the loss
function111Here we use a Mean Squared Error (MSE) function combined with an
Adam optimizer.. The architecture for this NN consists of a cell where the
output data, $\mu(z)$, of the previous step is used to compute the next one.
Each cell is provided with the prior information of the output value $\mu$
using

$\displaystyle h^{<t>}$ $\displaystyle=$ $\displaystyle g(W_{h}\cdot
h^{<t-1>}+W_{x}\cdot x^{<t>}+b_{a}),$ (4.1) $\displaystyle y^{<t>}$
$\displaystyle=$ $\displaystyle g(W_{y}\cdot h^{<t>}+b_{y}),$ (4.2)

where $b$ is the bias, $g$ is the activation function, $y^{<t>}$ is the output
and $h^{<t>}$ and $h^{<t-1>}$, are the hidden state and its preceding value
respectively. Moreover, RNN fails in long sequences due to the loss of prior
information from the initial inputs. To improve our training performance we
use a modified version of RNN as Long Short Term Memory (LSTM) cells. These
cells are capable to forget or add information step by step by incrementing
the number of matrix expressions for each layer as [34]:

$\displaystyle i^{<t>}$ $\displaystyle=$ $\displaystyle\sigma(W_{xi}^{T}\cdot
x^{<t>}+W_{hi}^{T}\cdot h^{<t-1>}+b_{i}),$ (4.3) $\displaystyle f^{<t>}$
$\displaystyle=$ $\displaystyle\sigma(W_{xf}^{T}\cdot x^{<t>}+W_{hf}^{T}\cdot
h^{<t-1>}+b_{f}),$ (4.4) $\displaystyle o^{<t>}$ $\displaystyle=$
$\displaystyle\sigma(W_{xo}^{T}\cdot x^{<t>}+W_{ho}^{T}\cdot
h^{<t-1>}+b_{o}),$ (4.5) $\displaystyle g^{<t>}$ $\displaystyle=$
$\displaystyle A_{f}(W_{xg}^{T}\cdot x^{<t>}+W_{hg}^{T}\cdot
h^{<t-1>}+b_{g}),$ (4.6) $\displaystyle c^{<t>}$ $\displaystyle=$
$\displaystyle f^{<t>}\otimes c^{<t-1>}+i^{<t>}\otimes g^{<t>},$ (4.7)
$\displaystyle y^{<t>}$ $\displaystyle=$ $\displaystyle h^{<t>}=o^{<t>}\otimes
A_{f}(c^{<t>}),$ (4.8)

where $W$ are the weights of each layer, $\sigma$ is the sigmoid function that
takes values between 0 and 1. Here $\otimes$ is the direct product and the
superscript $T$ denotes the transpose of the quantity where it is indicated.
The feed process is explicitly given in the second r.h.s term of the Eqs.
(4.3)-(4.6).

To compare different trained reconstructions of $\mu(z)$ given by the Pantheon
SNeIa and GRB data, we use three different activation functions $A_{f}$
defined as

$\displaystyle A_{f_{\text{Tanh}}}$ $\displaystyle=$
$\displaystyle\tanh(x),\quad\text{in}\quad(-1,1),$ (4.9) $\displaystyle
A_{f_{\text{ELU}}}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\alpha(e^{x}-1)&\mbox{for $x\leq
0$},\\\ x&\mbox{for $x>0$},\end{array}\right.\text{in}\quad(-\alpha,\infty),$
(4.12) $\displaystyle A_{f_{\text{SELU}}}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\alpha\lambda(e^{x}-1)&\mbox{for
$x\leq 0$},\\\ x&\mbox{for
$x>0$},\end{array}\right.\text{in}\quad(-\alpha\lambda,\infty).\quad\quad$
(4.15)

Furthermore, NNs have several parameters that can lead to a high probability
of overfitting. To overcome this issue, we improve our NN by using a
regularization method combined with Variational Dropout (VD). These allow us
to perform a regularization that can ‘turn off’ specific neurons to avoid
overfitting given certain values for the size, epochs, layers, neurons and
batchsize inside the NN architecture.

Figure 1:  Representation of our NN architecture. The input $z_{i}$ denotes
the redshifts and the output are the distance moduli $\mu_{i}$ for our
observable supernovae and GRBs. This NN has one layer, and the white squares
inside each node denote the number of neurons $h_{i}$. In this work we take
$h_{i}=100$. As indicated by the small arrows, each neuron receives input
information from the previous neuron. After many steps, there is some
processing on the input and the computed result is passed to a second neuron
connected through a different path, and so. Every node is associated with an
activation function (e.g Tanh, ELU and SELU) which transforms the input to an
output result. The output value of these functions acts as an input value to
the next connected node. Final outputs can decide the class of input data.

At this point, we can have some additional overfitting problems due to the
transport of prior information. To solve this issue, a Bayesian Neural Network
(BNN) is adapted to the NN, which can compute the errors on the outputs in the
form of posterior distributions, making the network probabilistic. The prior
distribution on the weight function for the input point $x$ can be calculated
by integrating

$p(\mathbf{y}^{*}|\mathbf{x}^{*},\mathbf{X},\mathbf{Y})=\int
p(\mathbf{y}^{*}|\mathbf{x}^{*},\mathbf{\omega})p(\mathbf{\omega}|\mathbf{X},\mathbf{Y})\mathbf{d\omega},$
(4.16)

where $p(\omega|\mathbf{X},\mathbf{Y})$ is the posterior distribution over the
space of parameters and in this case X and Y are the redshift $z$ and the
modulus distance $\mu$, respectively.

A supporting feature of our RNN+BNN network is the following: our selected
$A_{f_{\text{Tanh}}}$ is bounded in comparison to the other two activation
functions. This is an ansatz that can be important in physically interpreting
and representing the behavior of an observable, in our case our supernovae and
GRB data. Due to the specifications of the SNeIa sample, and their use to
calibrate the GRB data, we can expect that the observed differences in peak
luminosities (and the six parameters for the GRBs) are closely correlated with
observed differences in the shapes of their light curves.

###  4.1 The method: processing the trained data

Figure 2: Left Top: Recurrent Neural Network architecture. The input $z_{i}$
are the redshifts and the output $\mu_{i}$ are the luminosity distance moduli
$\mu_{i}$. This example figure consists of one layer and the grey circles
inside each yellow node indicate the number of neurons $h_{i}$. The values for
layers and neurons are given in (1). Each neuron receives input results from
the previous neurons at each step and an activation function is associated
step by step. Left Bottom: Bayesian Neural Network architecture. This consists
of one hidden layer with 12 nodes and an output layer with one node. The
weight functions between input and hidden layers are given in (4.1), where the
resulting prior distribution is given by (4.16). Right: Schematic summarizing
the steps in training the network and carrying out the fit to the data, to
compute the best-fit cosmographic parameters and their credible regions.

Figure 3: RNN+BNN training, up to $z=10$, using three different activation
functions: top left: Tanh activation function; top right: SELU activation
function; bottom: ELU activation function. Black points denote the GRB data
and blue points denote the Pantheon SNeIa data, while the light blue regions
represent the 1-$\sigma$ confidence contours on the reconstructed $\mu-z$
relation. Notice that there is good agreement between the reconstructed
relation and the high-z GRB data for the tanh activation function, but not for
the other activation functions.

Now we describe the methodology used to train the observational data given by
the Pantheon supernovae and GRBs sample using our RNN+BNN network. We analyze
three different architectures that use different activation functions. The
process was performed using TensorFlow222https://www.tensorflow.org. Figure 1
presents a schematic that outlines the logistics of the full method.

To perform the deep learning training we follow the methodology given in [1]
adapted with the Cosmography described in [24]. We divide our description of
this process into two sections: firstly, we describe the RNN+BNN ML
architecture and secondly, we describe how we calculated the best-fit values
for the cosmographic parameters that characterise the $d_{L}-z$ relation.

####  4.1.1  Training the Pantheon and GRBs data with RNN+BNN architecture

  1. (1)

Design of the neural network. We choose the three stated activation functions
with the following 8 hyperparameters: Size = 4; Epochs = 1000; Layers = 1;
Neurons = 100; Bathsize = 10; a variational dropout with an input $z$; a
hidden state $h$; an output $\mu(z)$.

  2. (2)

Organising the full sample of SNeIa + GRB data. We ordered the data from high
to low redshift, and adopted a number of steps $n=4$. The reason for this
ordering is that a RN feeds itself during the training – i.e. a network with a
neuron will have a connection from the input and also from the output of its
prior step, from which it feeds itself. This approach has the advantage of
recording more information to train in regions where the density of data is
higher (i.e. at low redshifts, for this full sample) up to regions where there
is a lower density of data points (at high redshifts, for this full sample).

  3. (3)

Computing the confidence regions via BNN. Due to the tendency of neural
networks to overfit, it is important to apply a form of regularisation – i.e.
to allow the algorithm to compute errors via a regularisation method. After
testing several times, we found that our models could not be trained with
input dropout due to the loss of information. To relax this prior, we use an
Adam optimizer.

  4. (4)

Extending the training to higher redshifts. After the final training, we
recover the model and apply $500$ times the same dropout to our previous
model. The predicted data sample with the above characteristics consists of
1209 data points in two redshift ranges: for observable ($0.01<z<9.3$)
redshifts, and for the higher (forecast) range, $0.01<z<10$. The results of
our training are shown in Figure 3, in which we present the reconstructed
$\mu(z)$ relation, using the Pantheon SNeIa and GRB data, with three different
activation functions and for the extended redshift range to $z=10$. Concerning
our choice of activation function, due to the physical behaviour that
characterises type Ia supernovae and GRBs – where one sees correlations
between luminosity and the shapes of their light curves – the tanh activation
function of Eq. (4.9) gave the better evolution to perform our cosmographic
analysis, which we carry out in Sec.5. As expected, this better performance is
seen in Figure 3.

##  5 Results: SNeIa and GRB high-$z$ cosmography

We now use the trained $\mu(z)$ relation, obtained following the methods
described in Sec.4, to carry out a cosmographical analysis by using Eq.(2.7)
to compute the cosmographical parameters that best fit this relation. At this
point we increased the number of epochs up to 1000. By doing so we are
changing the weights of the network; when we increase the number of epochs to
be the same as the number of times weights are changed in the NN, our analysis
switches from underfitting to optimally overfitting.

We compute the best-fit cosmographic parameters by modifying the publicly
available codes emcee 333https://emcee.readthedocs.io/en/stable/. Corner plots
showing 2-$\sigma$ credible regions for the fitted parameters are presented in
the left-hand panels of Figure 4. The upper left panel shows the results
obtained using the Pantheon SNeIa + GRB observational data, while the lower
left panel is for the trained data up to $z=10$ using the Tanh activation
function. The numerical results of these parameter fits are also summarised in
Table 4. The right hand panels of Figure 4 also show the $\mu(z)$ relation,
together with its 1-$\sigma$ error band, that corresponds to the best-fit
cosmographical parameters derived from each data set, compared with the
$\mu(z)$ relation that corresponds to the $\Lambda$CDM cosmographic model.
This latter model can be derived from e.g. Planck 2018 data, for which
$\Omega_{m}=0.313$, and assuming flatness; this implies that
$\Omega_{\Lambda}=0.685$ and hence that $q_{0}=-0.523$. Moreover, as noted in
e.g. [18], the $\Lambda$CDM model enforces $j_{0}=1$. Once with the values of
$q_{0}$ and $j_{0}$ thus determined, we can then fix them and fit again the
value of $s_{0}$ from our trained $\mu(z)$ relation to found $s_{0}$.

From Table 4 and Figure 4, we see that the deceleration parameter, $q_{0}$ is
quite well constrained from the observational SNIe + GRB data (upper left
panel of Figure 4) but is rather less well constrained by the trained data up
to $z=10$. The jerk and snap parameters, $j_{0}$ and $s_{0}$ are much less
well constrained by the observational data (although both are strongly
correlated with $q_{0}$ and with each other) and the $2-\sigma$ credible
regions on these parameters are again substantially larger for the trained
data extended to $z=10$ (see bottom part of Table 4). This pattern is
consistent with the $\mu(z)$ graphs shown in the right-hand panels of Figure
4, where the blue and orange curves show good agreement over the full redshift
range of observational data (top right panel), but the blue curve is more
divergent from the orange and green curves for the trained data extended to
$z=10$ (bottom right panel).

Since we are working with a cosmographic approach via $\mu(y)$, that may
result in a degeneracy in the derived cosmographic parameters, a
straightforward procedure to compare two models and their parameters is the
likelihood ratio test. To compute this comparison, we can use the quantity
$2\ln{\cal L}_{{\rm simple}}/{\cal L}_{{\rm complex}}$, where ${\cal L}$ is
the maximum likelihood of the cosmographic model, and the subscripts ‘simple’
and ‘complex’ could denote, for example, predictions of $\mu(y)$ to lower and
higher order respectively. The properties of the chi-squared distribution and
Jeffreys’ scale can then be employed to assess the significance of any
increase in likelihood against the number of extra parameters (or number of
data points) introduced.

As a less computationally intensive adaptation of this approach, we can also
measure which models are better by taking into account how many parameters the
models require to fit the data, and how good is that fit to the data, by
making use of the Akaike Information Criterion (AIC) or the Bayesian
Information Criterion (BIC): these criteria can be thought of as
approximations to the Bayesian evidence for each model that are less
computationally intensive to determine. Here we apply the BIC to compare our
cosmographic $\mu(z)$ relation with the corresponding relation for
$\Lambda$CDM.

We can compute the BIC using the relation

$\mbox{BIC}=\chi^{2}_{\rm{min}}+d\ln N,$ (5.1)

where $N$ is the number of data points and $d$ is the number of model
parameters. The quantity $\Delta\mbox{BIC}_{AB}=\mbox{BIC}_{A}-\mbox{BIC}_{B}$
can be interpreted as a measure of the evidence against model $A$ compared to
model $B$. For $0\leq\Delta\mbox{BIC}_{AB}<2$ there is not enough evidence
against model $A$; for $2\leq\Delta\mbox{BIC}_{AB}<6$ there is some evidence
against model $A$ and for $6\leq\Delta\mbox{BIC}_{AB}<10$ there is strong
evidence against model $A$.

Table 5 presents the results of our BIC calculations. We see that for the SNIe
+ GRB sample over the observed range a value of $\Delta\mbox{BIC}=4.123$ the
$\mu(z)$ relation differs markedly from that for $\Lambda$CDM. When we
consider the trained sample, extended to $z=10$, however, we obtain a much
larger value of $\Delta\mbox{BIC}>35$, indicating much stronger evidence of a
significant difference between the extended cosmographic $\mu(z)$ derived from
our data and the corresponding relation for $\Lambda$CDM. This divergence,
which was already apparent in Figure 4, arises because the high-redshift GRBs
are still highly correlated with their low-redshift SNIe calibrators.
Moreover, the divergence points towards the following robust conclusions from
our analysis: 1) the BNN method is insufficient to correctly calibrate the GRB
data, and 2) the RNN+BNN architecture is reliable enough to capture correctly
the physical trend that dimmer SNeIa decline more rapidly after maximum
brightness, and to incorporate that trend when extending our training set up
to $z=10$. In order to improve its description, however, we need better
quality calibration data up to $z=2.3$.

Figure 4: Left: 2-$\sigma$ credible regions for the cosmographic parameters
evaluated at the present cosmic time. Right: Evolution of $\mu(z)$ for
different cosmographic parameters. The green curves show the $\mu(z)$
relations for the cosmographic parameters that correspond to the $\Lambda$CDM
model, while the light blue bands show the best fit $\mu(z)$ curves, and their
1-$\sigma$ errors, using the cosmographic parameter values reported in Table
4. The top panels show the analysis using the Pantheon SNeIa + GRB
observational data, while the bottom two panels are for the analysis using
Pantheon SNeIa + GRBs trained data up to $z=10$ with activation function tanh.

Observational $z$ | Best-fit$\pm\sigma$ |  $68\%$ upper |  $68\%$ lower
---|---|---|---
$q_{0}$ | $-0.819_{-0.094}^{0.094}$ | $-0.725$ | $-0.913$
$j_{0}$ | $3.754_{-1.937}^{1.977}$ | $5.731$ | $1.817$
$s_{0}$ | $23.576_{-22.103}^{24.706}$ | $48.282$ | $1.473$

Trained $z=10$ | Best-fit$\pm\sigma$ |  $68\%$ upper |  $68\%$ lower
---|---|---|---
$q_{0}$ | $-0.410_{-0.255}^{0.249}$ | $-0.160$ | $-0.665$
$j_{0}$ | $-3.519_{-3.578}^{3.962}$ | $0.443$ | $-7.096$
$s_{0}$ | $-50.287_{-12.583}^{25.340}$ | $-24.947$ | $-62.870$

Table 4: Cosmographic best fits from (2.9) obtained with the 1-$\sigma$ using
the Pantheon SNeIa + GRBs sample and RNN+BNN trained sample. Top Table:
Results reported up to the observational redshift. Bottom Table: Results
reported up to the trained redshift.

Model (Observational) | $\chi^{2}$ | BIC |  $\Delta$BIC
---|---|---|---
$\Lambda$CDM cosmography | $1117.398$ | $1131.594$ | -
Cosmography (2.9) | $1107.324$ | $1135.717$ | $4.123$

Model (Trained $z=10$) | $\chi^{2}$ | BIC |  $\Delta$BIC
---|---|---|---
$\Lambda$CDM cosmography | $184.930$ | $200.126$ | -
Cosmography (2.9) | $12.279$ | $35.080$ | $165.046$

Table 5: Bayesian analysis obtained using the Pantheon SNeIa + GRBs full
sample. First column: The cosmographic models in comparison to $\Lambda$CDM
for our two redshift range: observational ($z:[0.01,9.3]$) and RNN+BNN trained
($z:[0.01,10]$). Second column: $\chi^{2}$-statistics for each model. Third
column: BIC analysis. Fourth column: $\Delta$BIC analysis. Top Table: Results
reported up to the observational redshift. Bottom Table: Results reported up
to the trained redshift.

##  6 Summary and Conclusions

In this work we have implemented a new method that uses neural networks (NNs)
to calibrate Gamma-ray bursts (GRBs) detected at high-$z$ as cosmological
distance indicators, combining them with a current SNeIa sample in order to
trace the Hubble expansion using an approximate cosmographic relation. Using a
combination of two networks – a Recurrent Neural Network (RNN) and a Bayesian
Neural Network (BNN) – we were able to reconstruct the distance modulus-
redshift relation, $\mu(z)$, of SNeIa and extend it to include the full GRB
sample. Our analysis was performed over a range that extends $\mu(z)$ to high
redshift, i.e. RNN+BNN: $z\in[0.01,10]$. Specifically, we simulated two types
of data with the BNN: first using the same redshift range as the observational
data and second when extending the training up to $z=10$.

This cosmographic approach required a very minimal set of assumptions that
does not rely on the dynamical equations for a specific theory of gravity. Our
NN analysis considered three different activation functions to train our
models: ELU, SELU, and Tanh. We found that the trend obtained by the SELU
function (4.15) does not follow the physical behaviour with the full sample of
SNeIa + GRB at higher redshift; therefore we discarded this activation
function. By comparing the ELU function (4.12) with the Tanh function (4.9),
we obtained a significant $\chi^{2}$-difference in the performance of these
two activation functions – with the Tanh function clearly favoured. Due to the
random initialization of weights in the RNN, we employed a cross-validation
method to find the model with the lowest $\chi^{2}$ for the Tanh function.

From the results of our analysis, as summarised in Figure 4, we can draw the
following conclusions:

  * •

For the trained SNeIa + GRB using RNN+BNN only up to the observed redshifts
reported in the literature (i.e. for Pantheon up to $z=2.3$ and for the GRB
sample up to $z=9.3$), the best-fit cosmographic relation obtained with these
data indicates that the training is reliable up to $z=2.3$. Specifically, the
best-fit $\mu(z)$ is essentially indistinguishable from the $\Lambda$CDM
cosmographic relation up to $z=2.3$ and deviates at only about 1-$\sigma$ from
the latter relation at higher redshift. We note that for $z<2.3$ the density
of data points is larger compared with that for the observable points at
higher redshifts. Hence, for $z>2.3$ the slight deviation of our best-fit
$\mu(z)$ relation from the cosmographic relation derived from the $\Lambda$CDM
parameters is likely the result of two factors: 1) the lower density of data
points and the larger errors on the estimated luminosity distance moduli of
GRBs at higher redshift, and 2) the possibility that the cosmographic relation
derived from the $\Lambda$CDM model parameters is inadequate to explain the
observable GRB data for $z>2.3$.

  * •

For the trained SNeIa + GRB, using RNN+BNN to extend the training up to
$z=10$, again we see that the training is fully reliable for $z<2.3$. At
higher redshifts, however, the reconstructed $\mu(z)$ starts to deviate more
strongly from the cosmographic relation derived from the $\Lambda$CDM model
than was the case for our analysis using only the observable redshift range.
This larger deviation may be due in part to our choice of activation function,
which sets a bias on the training itself. Thus, while the trained data over
the observed redshift range does provide support for $\Lambda$CDM as the
preferred model, this result may be biased due to the low density of data
points above $z>3.6$ (or even the lack of data between $z=7$ and $z=8$, see
Table 3), and we can expect that bias to be more important when the RNN+BNN
training is extended up to $z=10$. Notwithstanding this, of course the best-
fit $\mu(z)$ relation for $z>3.6$ may also be telling us that the standard
$\Lambda$CDM cosmography is inadequate at high redshift.

The work presented here represents a powerful, new NN methodology for
combining and calibrating SNeIa and GRB datasets, but also emphasises the
requirement of more real data at higher redshifts that can help the NN to
identify the physics trends underlying the cosmographic relations traced out
by these datasets. With such additional calibration data, we believe that this
new NN approach can successfully calibrate GRB distance information over an
extended redshift range, up to $z=10$, and thus can help to compare the
cosmography inferred from supernovae observations with the cosmography
inferred at higher redshift.

## Acknowledgments

CE-R acknowledges the Royal Astronomical Society as FRAS 10147. CE-R, MACQ and
CDZM are supported by DGAPA-PAPIIT-UNAM Project IA100220. M. H. is supported
by the Science and Technology Facilities Council (Ref. ST/L000946/1). This
article is also based upon work from COST action CA18108, supported by COST
(European Cooperation in Science and Technology).

## References

  * [1] C. Escamilla-Rivera, M. A. C. Quintero and S. Capozziello, JCAP 03 (2020) no.03, 008 doi:10.1088/1475-7516/2020/03/008 [arXiv:1910.02788 [astro-ph.CO]].
  * [2] K. Barbary, G. Aldering, R. Amanullah, M. Brodwin, N. Connolly, K. S. Dawson, M. Doi, P. Eisenhardt, L. Faccioli and V. Fadeyev, et al. Astrophys. J. 745 (2012), 32 doi:10.1088/0004-637X/745/1/32 [arXiv:1010.5786 [astro-ph.CO]].
  * [3] M. Betoule et al. [SDSS], Astron. Astrophys. 568 (2014), A22 doi:10.1051/0004-6361/201423413 [arXiv:1401.4064 [astro-ph.CO]].
  * [4] D. M. Scolnic, D. O. Jones, A. Rest, Y. C. Pan, R. Chornock, R. J. Foley, M. E. Huber, R. Kessler, G. Narayan and A. G. Riess, et al. Astrophys. J. 859 (2018) no.2, 101 doi:10.3847/1538-4357/aab9bb [arXiv:1710.00845 [astro-ph.CO]].
  * [5] Breuval, Louise and Kervella, Pierre and Wielgórski, Piotr and Gieren, Wolfgang and Graczyk, Dariusz and Trahin, Boris and Pietrzyński, Grzegorz and Arenou, Frédéric and Javanmardi, Behnam and Zgirski, B. The Astrophysical Journal 913\. DOI 10.3847/1538-4357/abf0ae.
  * [6] Lindegren, L. and Klioner, S. A. and Hernández, J. and Bombrun, A. and Ramos-Lerate, M. and Steidelmüller, H. and Bastian, U. and Biermann, M. and de Torres, A. and Gerlach, E. and et al. Astronomy & Astrophysics 2021.
  * [7] A. Cucchiara, A. J. Levan, D. B. Fox, N. R. Tanvir, T. N. Ukwatta, E. Berger, T. Kruhler, A. K. Yoldas, X. F. Wu and K. Toma, et al. Astrophys. J. 736 (2011), 7 doi:10.1088/0004-637X/736/1/7 [arXiv:1105.4915 [astro-ph.CO]].
  * [8] J. P. Norris, G. F. Marani and J. T. Bonnell, Astrophys. J. 534 (2000), 248-257 doi:10.1086/308725 [arXiv:astro-ph/9903233 [astro-ph]].
  * [9] L. Amati, F. Frontera, M. Vietri, J. J. M. in’t Zand, P. Soffitta, E. Costa, S. Del Sordo, E. Pian, L. Piro and L. A. Antonelli, et al. Science 290 (2000), 953-955 doi:10.1126/science.290.5493.953 [arXiv:astro-ph/0012318 [astro-ph]].
  * [10] E. Fenimore and E. Ramirez-Ruiz, [arXiv:astro-ph/0004176 [astro-ph]].
  * [11] L. Amati, F. Frontera and C. Guidorzi, Astron. Astrophys. 508 (2009), 173-180 doi:10.1051/0004-6361/200912788 [arXiv:0907.0384 [astro-ph.HE]].
  * [12] B. E. Schaefer, Astrophys. J. 660 (2007), 16-46 doi:10.1086/511742 [arXiv:astro-ph/0612285 [astro-ph]].
  * [13] G. Ghirlanda, G. Ghisellini, D. Lazzati and C. Firmani, Nuovo Cim. C 28 (2005), 303-306 doi:10.1393/ncc/i2005-10046-0 [arXiv:astro-ph/0504184 [astro-ph]].
  * [14] E. Di Valentino, L. A. Anchordoqui, O. Akarsu, Y. Ali-Haimoud, L. Amendola, N. Arendse, M. Asgari, M. Ballardini, S. Basilakos and E. Battistelli, et al. [arXiv:2008.11283 [astro-ph.CO]].
  * [15] E. Di Valentino, L. A. Anchordoqui, O. Akarsu, Y. Ali-Haimoud, L. Amendola, N. Arendse, M. Asgari, M. Ballardini, S. Basilakos and E. Battistelli, et al. [arXiv:2008.11284 [astro-ph.CO]].
  * [16] Y. L. Bolotin, V. A. Cherkaskiy, O. Y. Ivashtenko, M. I. Konchatnyi and L. G. Zazunov, arXiv:1812.02394 [gr-qc].
  * [17] S. Capozziello, M. De Laurentis, O. Luongo and A. Ruggeri, Galaxies 1, 216 (2013) doi:10.3390/galaxies1030216 [arXiv:1312.1825 [gr-qc]].
  * [18] S. Capozziello, R. Lakosz and V. Salzano, Phys. Rev. D. 84, 12 (2011) doi:10.1103/PhysRevD.84.124061 [arXiv:1104.3096 [astro-ph.CO]].
  * [19] C. Escamilla-Rivera and S. Capozziello, Int. J. Mod. Phys. D 28, no. 12, 1950154 (2019) doi:10.1142/S0218271819501542 [arXiv:1905.04602 [gr-qc]].
  * [20] E. K. Li, M. Du and L. Xu, Mon. Not. Roy. Astron. Soc. 491 (2020) no.4, 4960-4972 doi:10.1093/mnras/stz3308 [arXiv:1903.11433 [astro-ph.CO]].
  * [21] A. Aviles, C. Gruber, O. Luongo and H. Quevedo, Phys. Rev. D 86 (2012), 123516 doi:10.1103/PhysRevD.86.123516 [arXiv:1204.2007 [astro-ph.CO]].
  * [22] S. Capozziello, R. D’Agostino and O. Luongo, Int. J. Mod. Phys. D 28 (2019) no.10, 1930016 doi:10.1142/S0218271819300167 [arXiv:1904.01427 [gr-qc]].
  * [23] S. Bahamonde, K. F. Dialektopoulos, C. Escamilla-Rivera, G. Farrugia, V. Gakis, M. Hendry, M. Hohmann, J. L. Said, J. Mifsud and E. Di Valentino, [arXiv:2106.13793 [gr-qc]].
  * [24] C. Z. Munõz and C. Escamilla-Rivera, doi:10.1088/1475-7516/2020/12/007 [arXiv:2005.02807 [gr-qc]].
  * [25] D. J. Schlegel et al., arXiv:1907.11171 [astro-ph.IM].
  * [26] Bowman, J., Rogers, A., Monsalve, R. et al. Nature 555, 67–70 (2018).
  * [27] M. Corman, C. Escamilla-Rivera and M. Hendry, [arXiv:2004.04009 [astro-ph.CO]].
  * [28] [LIGO Scientific and Virgo], [arXiv:2004.08342 [astro-ph.HE]].
  * [29] E. Lusso, E. Piedipalumbo, G. Risaliti, M. Paolillo, S. Bisogni, E. Nardini and L. Amati, Astron. Astrophys. 628 (2019), L4 doi:10.1051/0004-6361/201936223 [arXiv:1907.07692 [astro-ph.CO]].
  * [30] C. Escamilla-Rivera and J. Levi Said, Class. Quant. Grav. 37 (2020) no.16, 165002 doi:10.1088/1361-6382/ab939c [arXiv:1909.10328 [gr-qc]].
  * [31] H. W. Chiang, A. E. Romano, F. Nugier and P. Chen, JCAP 11 (2019), 016 doi:10.1088/1475-7516/2019/11/016 [arXiv:1706.09734 [astro-ph.CO]].
  * [32] J. S. Wang, F. Y. Wang, K. S. Cheng and Z. G. Dai, Astron. Astrophys. 585 (2016), A68 doi:10.1051/0004-6361/201526485 [arXiv:1509.08558 [astro-ph.HE]].
  * [33] G. D’Agostini, [arXiv:physics/0511182 [physics]].
  * [34] Aurelien, G. Hands-On Machine Learning with Scikit-Learn and Tensorflow: Concepts, Tools, and Techniques to Build Intelligent Systems. ISBN-10:1491962291. O’Reilly Media. 2017.