---
author:
- 'S. D. P. Vitenti'
- 'M. Penna-Lima'
title: A general reconstruction of the recent expansion history of the universe
---

Introduction {#sec:introduction}
============

Many indications of the accelerated expansion of the universe come from distance measurements, such as the distance modulus of type Ia supernovae (SNe Ia) [@Riess1998; @Perlmutter1999]. In the last two decades, several models have been proposed in order to explain this phenomenon and, in general, they can be classified into dynamic and kinematic models. Assuming the general relativity, the first is described by adding a fluid, Dark Energy (DE), in which several propositions provide different DE equation of state (EoS) (for a review, see [@Joyce2015] and references therein). Other common dynamic approach is to modify the geometric setting of the gravitational theory instead of the energy-momentum tensor, such as the high-dimensional models [@Dvali2000] and $f(R)$ theories [@Sotiriou2010; @Nojiri2011]. These approaches are labeled as dynamic in the sense that there are differential equations of motion for the metric, whose modifications consist in altering the source term or the equation of motion itself.

In the context of kinematic models, the expansion history of the universe can be probed without assuming any theory of gravitation nor its matter content, and one only needs to define the space-time metric to study it. Considering the Friedmann-Lemaître-Robertson-Walker (FLRW) metric, the recent expansion of the universe is described in terms of the scale factor and its $n$-order derivatives with respect to time, such as the Hubble, deceleration and jerk functions [@Turner2002c; @Visser2004; @Rapetti2007; @Zhai2013], as well as in terms of the luminosity distance [@Daly2003; @Benitez-Herrera2012] (and references therein).[^1]

Since the only unknown in this metric is the scale factor (and possibly the spatial curvature), once one of the kinematic functions is determined, the others can be found by integrating and/or differentiating it. Therefore, all kinematic functions are related. Which function will be chosen to be reconstructed depends on the questions one wants to answer.[^2] For example, one can model the luminosity distance by a linear piecewise function and obtain a statistically sound and unbiased fit using observational data. However, this study will not contribute at all to the understanding of the recent accelerated expansion since, in this case, the deceleration function is assumed zero in the entire redshift interval.[^3]

Keeping the above idea in mind, the reconstruction of a kinematic function has been addressed using different methods. To understand some of these different methodologies, we divide the problem in two parts: (i) the theory underlying the observable quantities and (ii) the relation between the observables and the data, including the data probability distribution.

Regarding (i), the kinematic model provides a perfect description of the observables, not taking into account any noises nor errors in the measurements. For the sake of argument, suppose that (ii) is not part of the problem, i.e., the data is perfectly known. In this case, we could use a parametric function and adjust its parameters, such that the observable (hopefully) matches the data points, or we could use the data points to determine the observable function using, for example, interpolation. In this sense, we say that the analysis is model-dependent when a given parametric function is chosen a priori and model-independent when we use the data to determine it.

There are also two main procedures to treat (ii). We can assume which is the probability distribution of the data and, consequently, the only problem left is to determine the observable curve, which can be done in a model-dependent or independent way, as discussed above. In statistics texts, this is described as a parametric method. On the other hand, we can follow an even more conservative path and not impose a given probability distribution for the data. This way, known as non-parametric, also uses the data to reconstruct their own probability distribution.

In the model-dependent parametric approach, one assumes a priori a specific functional form of kinematic quantities, such as the deceleration function $q(z)$, and a probability distribution of the data [@Riess2004; @Shapiro2006; @Avgoustidis2009; @Nair2012]. A feature of this strategy is that its results have potentially smaller error bars when compared to the others. After all, one is introducing a reasonable set of assumptions which can lead to biased results. A natural improvement to this is to apply a model-independent approach, where one tries to reconstruct the curve when still using the assumed distribution for the data. Among these approaches is the Principal Component Analyses (PCA), in which the kinematic function is described in terms of a set of basis functions and the data is used to determine which subset of this basis is better constrained. Then, the function is reconstructed by using this subset [@Huterer2003; @Shapiro2006; @Clarkson2010; @Ishida2011; @Ruiz2012; @Benitez-Herrera2013].

Another possibility is to use smoothing methods [@Shafieloo2006; @Shafieloo2012]. In this model-independent non-parametric case, only mild assumptions are made about the data and, usually, no assumption is made about the model. This allows a direct translation of the data into a kinematic curve. Still in this context, we also have the Gaussian Process (GP), in which one chooses to model directly the probability distribution of the kinematic function itself [@Holsclaw2010; @Seikel2012]. For a more complete list of non-parametric methods see [@Montiel2014] and references therein.

Recovering both the probability distribution of the data and a reconstructed kinematic function require a large amount of data and, in practice, the current observational cosmology did not seem to have reached this level yet. This is evinced by the results obtained so far in the literature [@Ruiz2012; @Shafieloo2012; @Seikel2012; @Montiel2014]. Regarding the data, there is a good perspective to increasingly improve their probabilistic descriptions, since different error sources, such as the systematic ones, are being included in their modeling (e.g., [@Conley2011; @Betoule2014] ). This presents an additional challenge to the non-parametric methods, as they must incorporate all the error sources in their reconstruction.

Even in a model-independent and non-parametric approach, the estimated curves are not free from assumptions. Each method has some internal choices of parameters. Currently in the literature, these parameters are obtained using the observational data. However, as we usually have only one set of data, doing so will calibrate the method for this one particular realization of the data. In this case, there is no way to know if this calibration provides the best balance between bias and variance. This difficulty can be circumvented using different realizations of the **same** data set. For a given calibration, i.e., for a given choice of the internal parameters, the method is applied to a large number of simulations obtaining the bias and variance for this calibration. Then, repeating this process for different calibrations one can find the best suited one for the chosen data set. In other words, the internal parameters in these reconstructions must not be related to one particular realization of the data, but to their probability distribution.

This idea can be extended to the study of the statistical properties of the data. For example, in [@Montiel2014], among other results, the authors apply a bootstrap-like procedure to calibrate the smoothing parameter applied to the data. This kind of analysis can provide a insightful information about the statistical properties of the data when little is known about their relationships.

In this work, we use the current available observational data for small redshifts ($z \lesssim 2.3$) and their likelihoods, namely, the Sloan Digital Sky Survey-II and Supernova Legacy Survey 3 years (SDSS-II/SNLS3) combined Joint Light-curve Analysis (JLA) SNe Ia sample [@Betoule2014], baryon acoustic oscillation (BAO) data [@Beutler2011; @Padmanabhan2012; @Kazin2014; @Ross2014] and $H(z)$ measurements [@Stern2010; @Riess2011; @Moresco2012; @Busca2013]. Currently, there is not enough data to perform a full model-independent and non-parametric reconstruction of the recent evolution of the universe. Therefore, we use the usual likelihood for these data, but, to be conservative, we reconstruct $q(z)$ along with some astrophysical parameters of SN Ia, the drag scale (present in the BAO likelihood) and the Hubble parameter $H_0$.

Besides the above data, there is also a wealth of data concerning the large scale structure connected to the perturbations around a FLRW metric, such as the temperature fluctuations of the cosmic microwave background (CMB) [@Hinshaw2013; @Planck2015]. Since we assume no dynamic model, we would have to propose a kinematic one for the perturbations. Such model is not feasible as it would require a set of functions of both time and space. In principle, one could also directly use derived observables, as, e.g., the CMB distance priors [@Komatsu2011], to fit the background model. However, these parameters are obtained in the context of a specific model, e.g., $\Lambda$CDM. Thus this model would be indirectly reintroduced in the results. For this reason, we choose not to use these data and focus on the distance-like measurements. We reconstruct $q(z)$ using a model-independent parametric approach, which consists in describing $q(z)$ by a cubic spline. The choice of the deceleration function comes from the fact that one can use such function to directly test the energy conditions [@Lima2008; @Lima2008a]. Likewise, the deceleration function is related to the underlying dynamics of the metric, since it is a simple combination of first and second derivatives of the scale factor. Therefore, their knowledge is necessary to constrain models which chooses a different dynamic for the gravitation sector (Alves et al. to be submitted). In addition, the spline method allows us to vary the complexity of the functional form as a function of the knots number, for example. Here we introduce a novel method to continuously vary the complexity of the reconstructed function. As a result of this analysis, we also obtain the reconstruction of those SN Ia parameters and the drag scale.

The paper is organized as follows. In section \[sec:def\_q\] we review basic concepts, such as $q(z)$, and the few assumptions made here. In section \[sec:review\], we discuss some approaches used in the literature and some of their respective drawbacks. In section \[sec:reconst\] we present our reconstruction method and the tools to vary the function complexity and to select the best one. Following, section \[sec:data\], we specify the observational data sets used in our study and their respective likelihood functions. We then perform a Monte Carlo analysis in different scenarios calibrating the method (section \[sec:meth\_valid\]). Finally, we use this calibration to reconstruct the deceleration function using the observational data, in section \[sec:real\_reconst\], and we summarize our conclusions in section \[sec:conclusions\].

Deceleration function $q(z)$ {#sec:def_q}
============================

Measurements of the CMB [@Hinshaw2013; @Planck2015] show that the universe is nearly homogeneous and isotropic at large scales. Therefore, we assume that the universe follows the cosmological principle, restricting the metric to the FLRW metric, $$\label{eq:rw_metric} 
ds^2 = -c^2\,dt^2 + a^2 (t) \left [ dr^2 + S_k^2(r) (d\theta^2 + \sin^2 \theta  d\phi^2) \right ]\;,$$ where $S_k(r)=(r\,$, $\sin(r)$, $\sinh(r))$ for flat, spherical and hyperbolic spatial section ($k=0, 1, -1$), respectively, $c$ is the speed of light and $a(t)$ is the cosmological scale factor. In this case, the expansion history of the universe can be defined knowing $a(t)$ and $k$.

In practice, we do not measure $a(t)$ directly, but related quantities such as the distances to astronomical objects. Considering a null trajectory of photons emitted by a galaxy traveling along the radial direction to us, we have that $$\label{eq:com_dist1}
r = c \int_{t_e}^{t_0} \frac{{\mathrm{d}}t^\prime}{a(t^\prime)},$$ where $t_e$ and $t_0$ are the emitted and observed times, respectively. Expanding the scale factor to second order around $t_0$, gives $$a(t) = a_0 + H_0 (t - t_0) - \frac{q_0 H_0^2}{2} (t - t_0)^2,$$ where $a_0$ is the scale factor today and $H_0$ and $q_0$ are, respectively, the Hubble and deceleration functions at $t_0$, $$H_0 = \left(\frac{ \dot{a}}{a} \right)_{t_0} \quad \text{and} \quad q_0 \equiv - \left(\frac{\ddot{a}a}{\dot{a}^2}\right)_{t_0}.$$ Rewriting the comoving distance $D_c$ and the deceleration function in terms of the redshift, $1 + z = a_0/a$, we obtain $$\label{eq:com_dist}
D_c(z) = a_0 r = \frac{c}{H_0} \int_0^z \frac{{\mathrm{d}}z^\prime}{E(z^\prime)}$$ and $$\label{eq:qz}
q(z) = \frac{(1 + z)}{H(z)}\frac{\mathrm{d}H(z)}{\mathrm{d}z} - 1,$$ whose integral solution is $$\label{eq:Ez}
E(z) = \exp{\int_0^z \frac{1+q(z^\prime)}{1+z^\prime}{\mathrm{d}}z^\prime},$$ and where $E(z)$ is the normalized Hubble function, $E(z) = H(z)/H_0$.

Equations  – evince that we shall reconstruct $q(z)$ in order to access local information about the accelerated/decelerated phase [@Lima2008]. Besides, assuming that $q(z)$ is continuous guarantees that $E(z)$ and $D_c(z)$ are also continuous functions and at least once and twice differentiable, respectively, although the opposite is not true.

Review of other approaches {#sec:review}
==========================

In this section we briefly present some of the most used methods in reconstructing the expansion history of the universe. We discuss some intrinsic issues of these methodologies, which motivated us to develop a novel approach (presented in section \[sec:reconst\]).

Parametric models {#sec:PM}
-----------------

In general the Taylor series approaches have two related problems. The first is the convergence radius of the Taylor expansion itself, which can only be estimated since the real scale factor is naturally unknown. On the other hand, as we are fitting the coefficients of such expansion, the functional form for the scale factor (or distance) can be interpreted as a simple polynomial interpolation. In this sense, changing the time parameter can be useful [@Cattoeen2008] and provide a better polynomial interpolation. However, this leads us to the second problem, the Runge’s phenomenon. That is, after a given order, higher order polynomials provide worse and worse approximations. Therefore, when using the Taylor expansion one should stay on a small convergence radius, which would restrain the analysis to a very small but unknown redshift or use a polynomial interpolation keeping in mind its caveats.

More generally, the problem of finding a good kinematic description of the expansion history can be addressed using a parametric method. In this context, one assumes a priori a specific functional form of a kinematic function, like the polynomial form discussed above, then proceeds by fitting its parameters using observational data. For example, in Refs. [@Riess2004; @Shapiro2006; @Avgoustidis2009; @Nair2012] they fit different functional forms of the deceleration function $q(z)$. The drawback of this method is that the choice of a functional form introduces a form-bias in the estimates if the functional form is different from the true one. As we do not know it, the result of such fit can be misleading since, even if the parameters’ error bars are small, their form-biases can still be large. A more conservative approach is to use flexible functional forms. However, this translates in using many parameters and, consequently, obtaining larger error bars. In this way, there is a natural trade-off between variance and form-bias which should be evaluated to determine the optimal reconstruction.

An additional, less discussed, difficulty is the estimator-bias. The functional forms are usually fitted using a Least-Squares (LS) or Maximum-Likelihood (ML) approach. As it is well known, both approaches can provide biased estimators for the parameters (ML estimators are usually asymptotically unbiased). This means that, even if we knew the correct functional form, the fact that we have only a finite number of observations can lead to estimator-biases.[^4] Heuristically, when fitting a functional form with $n$ parameters using $m$ data points we will have $m/n$ observations per point. Hence, if the estimators are only asymptotically unbiased, then the higher the number of parameters higher the estimator-biases. Thus, the final trade-off must consider variance, form-bias and estimator-bias.

Principal Component Analysis {#sec:PCA}
----------------------------

Another popular methodology is the Principal Component Analysis (PCA) [@Huterer2003; @Shapiro2006; @Clarkson2010; @Ishida2011; @Ruiz2012]. In this approach the kinematic function is described by a flexible parametrization (e.g., in [@Huterer2003] they express the DE EoS as a constant piecewise function within $N=50$ bins). From the estimated covariance matrix of these parameters, the eigenvectors and eigenvalues are deduced. Say, for example, that we choose the first 5 eigenvectors, whose eigenvalues correspond to the smallest variance terms. This means that the original division in 50 bins is being described by a 5-dimensional parametrization.

The eigenvectors, whose eigenvalues correspond to the smallest variance terms, provide the parametrization which is better constrained by the data. The appeal of this method is that it provides a straightforward way to determine the curves better constrained by the data. The choice of how many eigenvectors should be used can be answered by looking the variance-bias trade-off [@Huterer2003]. On the other hand, the method as described above is subject to a potential drawback. The relationship between different kinematic functions are given in terms of their integral in time, for example, the relation between the deceleration and the Hubble functions is given by Eq.  while the distance is another integral of the Hubble function \[Eq. \]. Note that, to calculate the distance to a given SNIa at redshift $z_i$, the deceleration function is integrated twice between $[0, z_i]$. Therefore, all bins in this interval contribute to the final value of the distance. This leads to the following problem, the first bins contribute to the value of the distance for almost all supernovae, while the final bins affect only few objects. Besides, since the initial bins always contribute to the value of the distance at higher redshifts, they will be naturally correlated to them, any modification in their value will have to be compensated by the other. These facts have the following consequence, the best constrained modes will be strongly connected to the first bins while the worst will be related mostly to the last bins. Therefore, when one chooses to use only a few (better constrained) eigenvectors, the final parametrization will provide almost no power in the last bins. This is natural since it is equivalent to choose the coefficients of the last eigenvectors to be fixed at zero. This problem can be seen in [@Huterer2003], where the issue persists even when the authors consider a forecast with 3000 SNe Ia uniformly distributed in redshift. This problem was also noted in [@Ishida2011]. In the latter, they realized that the first modes have this deficiency and proposed the addition of a new parameter to circumvent this problem.

Similarly to what is commonly used in the PCA approach, the authors in [@Crittenden2012; @Zhao2012] reconstructed the EoS function describing it as a set of bins (discontinuous piecewise constant function of the scale factor). Instead of limiting the estimated curve variance by using a small set of eigenvectors (principal components), they introduced the correlated prior, in an approach similar to the GP, which treats the reconstructed curve as a set of random variables emerging from a multidimensional Gaussian distribution. This prior introduces a correlation between the bins controlled by the correlation length and prior strength. Once the prior is calibrated, the problem described in the above paragraph can be systematically resolved. The discussion above elucidates an important characteristic of the PCA, it is dependent on the initial description of the function. In many works the simple step function (or bins) were used. If another basis functions were used, with different behavior at large redshifts, we would end up with a similar problem, i.e., the large redshift behavior of the best constrained basis function would dominate at large redshifts, providing biased estimates at these points. We conclude that, in the specific context of constraining a function by its integral, the PCA approach has this potential problem.

Smoothing methods {#sec:SM}
-----------------

As we mentioned in section \[sec:introduction\], these approaches are model-independent and non-parametric. Therefore, they suit the study cases where the data probability distribution is unknown. However, as pointed out by Montiel et al. [@Montiel2014], these methods find difficulties due to the limited amount of observational data or even due to some features of the methods themselves, such as the size of the smoothing parameter and the assumptions on priors or fiducial models. Besides, since no assumption is made about the data distribution, one cannot use resampling to perform a self-validation, but only bootstrap like procedures, e.g., jackknife, which usually requires large samples.

Finally, the fact that this method makes such minimal assumptions is not necessarily useful. As one can only reconstruct the direct variable associated to the observable, i.e., cosmological distances and $H(z)$, any inference about the derived kinematic quantities is limited since it is highly dependent on the smoothing technique as, e.g., smoothing linear splines has always zero second derivative, and top-hat moving average filters are non-continuous. Therefore, any analysis about kinematic quantities, different of the reconstructed one, requires new assumptions [@Daly2004; @Shafieloo2012].

Gaussian Process {#sec:GP}
----------------

In this approach, instead of modeling a kinematic function, one chooses to model a probability distribution for the curve as a Gaussian probability distribution. This assumption dictates the data probability distribution by relating both the curve and observable probability distributions. In this sense, this approach unifies the two aspects of the model, the data distribution and the curve reconstruction. All the assumptions are comprised in the mean curve and the two-point covariance, which define the Gaussian distribution of the curve.

The drawback is similar to the parametric procedure. In general, one has to assume a function to describe the mean of the GP and a two-point function to describe the variance. If the considered mean function differs from the true one, it will impose a bias in the reconstruction. See, for example, references [@Holsclaw2010; @Seikel2012], where they assume a constant mean function. Since the GP determines the observable statistical distribution, it is also necessary to include the data distribution to calculate the joint probability distribution of both curve and data. In this case, one has a model-independent but parametric method in the sense that one is assuming a particular distribution for the data. The GP validation also has to be performed for a number of realizations of the data, since the indirect curve determination through a Gaussian distribution can lead to bias in an unpredictable way.

Reconstruction of $q(z)$ {#sec:reconst}
========================

Broadly speaking each strategy described in section \[sec:review\] has a better suited application. As a rule of thumb, to use less hypothesis it is necessary to have more data. Therefore, if the amount of data is limited, the generality of GP and soothing methods, for example, are restricted, leading to underdetermined problems. Besides that, including natural hypothesis can also be considerably difficult in those non-parametric studies, e.g., after the determination of the cosmological distances through a smoothing method it is necessary to add new assumptions to describe the Hubble function.

Therefore, in this work we adopted a model-independent parametric approach. Nonetheless, to be conservative we reconstruct the kinematic curve along with all the phenomenological parameters related to the modeling of each data set. Doing so, we minimize the assumptions on the data distribution bypassing any bias which could result from it. In a model-independent technique we need to use a set of functions to perform the reconstruction. To avoid the problems described above on the PCA approach, we use a cubic spline to reconstruct our observable, as described in section \[sec:spline\]. One advantage of the cubic spline is that it is continuous and twice differentiable on every knot, thus, all the parameters (the value of the function on each knot) are related through these conditions.[^5]

Another point on the model-independent approach is how to deal with the bias. In the PCA approach one uses the particular set of data to determine which basis functions form the minimal set. Depending on the data distribution, this tends to favor the intervals where the sample is denser, which usually creates biases on the other regions. To avoid this issue, we impose a global penalization (see section \[sec:penalty\]) on the curvature of the curve. Applying the same penalization in each interval, we evade the localization problem described before.

Piecewise deceleration function {#sec:spline}
-------------------------------

In this work, we avoid making arbitrary choices of the $q(z)$ form, and, consequently, *a priori* restricting it to specific functional forms, by approximating $q(z)$ by a piecewise third-order polynomial function, i.e., a cubic spline.

The first step to build an estimator of $q(z)$ \[denoted as $\hat{q}(z)$\] is to specify the redshift interval (domain $\mathrm{D}$) in which the function is defined. This interval is $\mathrm{D} = [z_{min}, z_{max}]$, where $z_{min}$ ans $z_{max}$ are the minimum and maximum redshifts of the used data. The next step is to choose the partition of the domain $\mathrm{D}$ into $n$ sub-intervals in which we define $\hat{q}(z)$ as a cubic polynomial function, namely, $$\hat{q}(z) = \left\{ \begin{array} {llllll}
             p_0(z) = a_0(z-z_0)^3 + b_0(z-z_0)^2 + c_0(z-z_0) + d_0 \quad & {z \in [z_0,z_1)}\\
             p_1(z) = a_1(z-z_1)^3 + b_1(z-z_1)^2 + c_1(z-z_1) + d_1 & {z \in [z_1,z_2)}\\
             \qquad\qquad\qquad\qquad \vdots & \quad \vdots\\       
             p_{n-1}(z) = a_{n-1}(z-z_{n-1})^3 + b_{n-1}(z-z_{n-1})^2 + c_{n-1}(z-z_{n-1}) + d_{n-1} & z \in [z_{n-1}, z_n],
     \end{array} \right.$$ where $z_0 = z_{min}$, $z_n = z_{max}$ and $p_i$ is the cubic polynomial defined in the $i$-th sub-interval. Note that each polynomial $p_i(z)$ in the segment $[z_i,z_{i+1})$ depends on 4 parameters ($a_i$, $b_i$, $c_i$ and $d_i$) and, consequently, we would need to estimate $4n$ parameters to define $\hat{q}(z)$ in the whole domain $\mathrm{D}$. However, imposing the following continuity conditions, $$\begin{aligned}
p_{i+1}(z_{i+1}) &=& p_i(z_{i+1}),  \\
p^{\prime}_{i+1}(z_{i+1}) &=& p^{\prime}_i(z_{i+1}), \\
\quad p^{\prime\prime}_{i+1}(z_{i+1}) &=& p^{\prime\prime}_i(z_{i+1}),\end{aligned}$$ on the $n-1$ internal knots $i \in (1,n-1)$, where ${}^\prime$ denotes the derivative with respect to $z$, and the two not-a-knot boundary conditions $$\hat{q}^{\prime\prime\prime}_{0}(z_{1}) = \hat{q}^{\prime\prime\prime}_1(z_{1}) \quad \text{and} \quad \hat{q}^{\prime\prime\prime}_{n-2}(z_{n-1}) = \hat{q}^{\prime\prime\prime}_{n-1}(z_{n-1}),$$ we end up with only $n + 1$ parameters to determine $\vec{Q} = \{d_0, ..., d_{n}\}$, i.e., the values of $\hat{q}(z)$ at each knot. As we have a one to one relation between the function in each knot $\hat{q}_i \equiv \hat{q}(z_i)$ and the parameter $d_i$ of each polynomial, from now on we rename the set as $\vec{Q} = \{\hat{q}_i\}$ for the sake of notational simplicity. Thus, using this cubic spline approximation, we fit the vector $\vec{Q}$ in order to obtain the estimates of the deceleration function (see description in section \[sec:meth\_valid\]).

Any interpolation method introduces an error source limiting the set of functions able to be reconstructed. The interpolation error for a cubic spline has an upper bound proportional to both the largest distance between nearby knots to the fourth power and the fourth derivative of the function with respect to $z$ (for details see [@Boor2001]). At first sight, larger the number of knots, smaller the interpolation error. But in practice, the number of knots is limited by the increasing number of parameters. Besides, there is also the overfitting that rises when fitting the parameters using data points. As a rule of thumb, one should choose $n$ and, consequently, the intervals between knots, such that the estimated function is expected to be well approximated by a cubic polynomial in these intervals. One can test the choice of $n = n_1$ applying the reconstruction for another one, e.g., $n = n_1 + 1$, and probing the results for any significant improvements on the fit. Notwithstanding the interpolation error, another important source of uncertainty comes from the statistical errors (bias/overfitting), as we will discuss in the next section. Finally, the use of cubic splines represents a large advantage in comparison to the step functions frequently used in the literature. For a constant piecewise function, the interpolation error is bounded by the first derivative times the largest distance between nearby knots. As a result, the number of knots (and, consequently, parameters) necessary to reconstruct functions with the same interpolation error bound is much larger in a binned approach.

Function complexity {#sec:penalty}
-------------------

Assuming a cubic spline to approximate $q(z)$, we are able to address both model-dependent and model-independent parametric methods by varying the number of knots. The simplest case, $n = 4$, is equivalent to consider that $q(z)$ is a third-order polynomial. On the other hand, we approach a model-independent case increasing $n$. The complexity of the function $\hat{q}(z)$ is, in principle, parameterized by the number of knots, as the number of knots goes to infinity any interpolation error drops to zero. Nonetheless, the choice of the domain partition is rather arbitrary and, at first, one would have to test different options in order to achieve, for example, a “model-independent limit” trying to minimize the over-fitting error. Another difficulty inherent of this approach is that the number of knots is a discrete variable and, as such, it is difficult to include it as another parameter in the analysis.

Instead of varying the number of knots by adding/removing actual knots to the function representation, we can fix the number of knots in some large value and penalize independent values of the parameters. For example, given that the parameters are just the values of the function at each knot, a penalty factor as a increasing function of $\vert\hat{q}_i - \hat{q}_{i+1}\vert$ will correlate all parameters $\hat{q}_i$. The correlation will be proportional to the weight of the penalization in the analysis, higher the weight more correlated are the parameters. In the strong correlation limit, all parameters would be equal, i.e., $\hat{q}_i = \hat{q}_{i+1}$. In this last case, even with a large number of knots, the effective number of degrees of freedom would be one. In short, varying the weight of the penalization, we can vary the effective number of degrees of freedom, circumventing the difficulties described above.[^6]

In practice, we include a set of penalty factors $P_i(\sigma_i)$ in our estimator. The initial likelihood is $L(\vec{D}, \vec{\theta})$, where the vector $\vec{D}$ represents the data set and the vector $\vec{\theta}$ all the parameters, including the spline parameters $\hat{q}_i$ and other parameters as described in sections \[sec:mc\_sneia\] and \[sec:mc\_alldata\]. To obtain the parameter estimators, we add to the likelihood $L(\vec{D}, \vec{\theta})$ the penalization $P_i(\sigma_i)$ defining the penalized likelihood $$\label{eq:likel_penalty}
-2 \ln\left[L_P (\vec{D}, \vec{\theta})\right] \equiv -2 \ln\left[L\left(\vec{D}, \vec{\theta}\right)\right] + \sum_{i = 2}^{n-1} P_i(\sigma_i),$$ where the penalty factor is given by $$\label{eq:penalty_fac}
P_i(\sigma_i) = \left(\frac{{\bar{\hat{q}}_i - \hat{q}_i}}{\sigma_i}\right)^2,\qquad \bar{\hat{q}}_i = \frac{(\hat{q}_{i-1} + \hat{q}_{i+1})}{2}, \qquad \sigma_i = \sigma_{abs} + \bar{\hat{q}}_i{\sigma_{\text{rel}}},$$ and we use $\sigma_{abs} = 10^{-5}$.[^7] The penalization factor is schematically illustrated in figure \[fig:basis\] showing the positions of $\hat{q}_1$ and $\bar{\hat{q}}_1$.

We control the complexity of $\hat{q}(z)$ by varying the value of the relative error ${\sigma_{\text{rel}}}$. For example, we are able to recover a high complexity function, in particular, a full $n+1$ knots spline for large ${\sigma_{\text{rel}}}$, and a straight line in the entire redshift interval when ${\sigma_{\text{rel}}}$ goes to zero. The former has many coefficients and can tend to fit the data noise, i.e., it is over-fitting dominated. The second naturally sharpen the constraints on $\{\hat{q}_i\}$, but they can be biased if the assumed functional form significantly differs from the true one. It is worth mentioning that this penalty factor allows us to explore a wide range of functional forms, since its simplest case is a linear function. Meanwhile, without using the penalty factor, the simplest model would be a third-order polynomial.

Finally, we emphasize that, in principle, one could use a large set of knots while constraining the allowed shape with the penalty factor. The restriction will be practical, the computational cost increases with the number of knots. Therefore, one should find the best balance between computational cost and flexibility of the method.

![\[fig:basis\] A descriptive example of the penalization on a cubic spline. For each three knots, the penalization is proportional to the distance between the straight line connecting the first and third knots, and the interpolation function.](penal_example.pdf)

Bias-variance trade-off {#sec:bv_trade}
-----------------------

Giving the penalty function which allows us to explore different complexity forms of $\hat{q}(z)$, we now have to introduce a criterion to determine the best ${\sigma_{\text{rel}}}$ value. We want to find the scenario where the combined error due to the biases (estimator- and form-bias) and the over-fitting error is minimized. For this, we decompose the error into bias and variance components as described below [@Bishop1996; @Keijzer2000; @Wasserman2001].

We create a controlled environment introducing a fidicual deceleration function $q^{{\text{fid}}}(z)$, which is determined by a given set of values for the parameters $\vec{\theta}^{{\text{fid}}}$. The idea is to use this function to generate a new data set $\vec{D}$, which is possible since we know the likelihood of the data $L(\vec{D},\vec{\theta})$. We define the ML estimators using the penalized likelihood, then, given a simulated sample $\vec{D}^{(l)}$, we obtain the estimates $\tilde{\vec{\theta}}(\vec{D}^{(l)})$ computing $$\frac{\partial L_P}{\partial \vec{\theta}}\left(\vec{D}^{(l)}, \tilde{\vec{\theta}}\right) = 0.$$ This provides an implicit definition of the function $\vec{\theta}^{(l)} \equiv \tilde{\vec{\theta}}(\vec{D}^{(l)})$. In principle, we could calculate the bias in the estimator integrating the function $\vec{\theta}^{(l)}$, i.e., $$\left\langle\tilde{\vec{\theta}}\right\rangle = \int{\mathrm{d}}\vec{D} \, \vec{\theta}^{(l)}L_P\left(\vec{D},\vec{\theta}^{{\text{fid}}}\right).$$ However, such integration is computationally unfeasible. It is a $N$ dimensional integration, where $N$ is the number of data points, and the function $\vec{\theta}^{(l)}$ is usually only determined numerically by maximizing the penalized likelihood.

Instead, we use the Monte Carlo (MC) approach to deal with such integrals.[^8] Since we know the probability distribution of the data, we can create a new realization of the data (mock catalog) $\vec{D}^{(l)}$ by resampling,[^9] i.e., using a (pseudo)random number generator to create a new data set. Given a large enough number of realizations $m$, we can approximate the expected value of a function of the data as $$\left\langle f\left(\vec{D}\right) \right\rangle \simeq \frac{1}{m} \sum_{l=1}^{m} f\left(\vec{D}^{(l)}\right).$$

Using these tools, we introduce the mean squared error (MSE) of approximating a fiducial function $q^{{\text{fid}}}(z)$ by $\hat{q}(z; {\sigma_{\text{rel}}})$. For a fixed ${\sigma_{\text{rel}}}$, it is $$\text{MSE} = \left\langle\left[\hat{q}^{(l)}(z; {\sigma_{\text{rel}}}) - q^{{\text{fid}}}(z)\right]^2\right\rangle \simeq \frac{1}{m} \sum_{l=1}^{m} \left[\hat{q}^{(l)}(z; {\sigma_{\text{rel}}}) - q^{{\text{fid}}}(z)\right]^2.$$ Given the estimate of the expected value, $$\label{eq:expec_val}
\langle \hat{q}(z, {\sigma_{\text{rel}}}) \rangle \simeq \frac{1}{m} \sum_{l=1}^{m} \hat{q}^{(l)}(z; {\sigma_{\text{rel}}}),$$ we have that (omitting ${\sigma_{\text{rel}}}$ for simplicity) $$\begin{aligned}
\left[\hat{q}^{(l)}(z; {\sigma_{\text{rel}}}) -  q^{{\text{fid}}}(z)\right]^2 &= \left[\hat{q}^{(l)}(z) - \langle \hat{q}(z) \rangle + \langle \hat{q}(z) \rangle -  q^{{\text{fid}}}(z)\right]^2 \nonumber \\
&= \left[\hat{q}^{(l)}(z) - \langle \hat{q}(z) \rangle \right]^2 + \left[\langle \hat{q}(z) \rangle -  q^{{\text{fid}}}(z)\right]^2 \nonumber \\
& + 2 \left[\hat{q}^{(l)}(z) - \langle \hat{q}(z) \rangle\right]\left[\langle \hat{q}(z) \rangle -  q^{{\text{fid}}}(z)\right],\end{aligned}$$ and, therefore, $$\label{eq:mse_1}
\text{MSE} \simeq \left(\langle \hat{q}(z) \rangle -  q^{{\text{fid}}}(z)\right)^2 + \frac{1}{m} \sum_{l=1}^{m} \left(\hat{q}^{(l)}(z) - \langle \hat{q}(z) \rangle \right)^2,$$ since $$\frac{1}{m} \sum_{l=1}^{m} \left(\hat{q}^{(l)}(z) - \langle \hat{q}(z) \rangle\right)\left(\langle \hat{q}(z) \rangle -  q^{{\text{fid}}}(z)\right) \simeq 0.$$ The first term of eq.  is the squared bias $b_{\hat{q}}(z)^2$ and the second is the variance. Since this variance estimator is biased, in this work, we evaluate the bias-variance trade-off computing $$\begin{aligned}
\label{eq:mse_2}
\text{MSE} &\simeq \left(\langle \hat{q}(z) \rangle -  q^{{\text{fid}}}(z)\right)^2 + \frac{1}{m -1} \sum_{l=1}^{m} \left(\hat{q}^{(l)}(z) - \langle \hat{q}(z) \rangle\right)^2 \nonumber \\
&\equiv b_{\hat{q}}(z; {\sigma_{\text{rel}}})^2 + \text{Var}\left(\hat{q}(z; {\sigma_{\text{rel}}})\right).\end{aligned}$$ Both the squared bias and the variance have the same weight in the above expression. Nonetheless, when applying the reconstruction for real data, we usually do not have access to an estimate of the bias. Therefore, minimizing MSE can potentially lead us to a methodology with a large bias (as we will see in sections \[sec:mc\_sneia\] and \[sec:mc\_alldata\]). To avoid this, in what follows we will minimize the MSE satisfying the constraint $$\label{eq:m_b}
\frac{b_{\hat{q}}(z; {\sigma_{\text{rel}}})}{\sigma(\hat{q}(z; {\sigma_{\text{rel}}}))} \leq m_b, \qquad \sigma(\hat{q}(z; {\sigma_{\text{rel}}})) \equiv \sqrt{\text{Var}(\hat{q}(z; {\sigma_{\text{rel}}}))},$$ where $m_b$ controls the maximum ratio between bias and variance. The variance of the reconstructed curve $\hat{q}(z; {\sigma_{\text{rel}}})$ can be written in terms of the covariance of the spline parameters $\text{Cov}\left(\hat{q}_i, \hat{q}_j\right)$. In turn this covariance can be estimated using the unbiased covariance estimator $$\label{eq:cov_q}
\text{Cov}\left(\hat{q}_i, \hat{q}_j\right) = \frac{1}{m - 1} \sum_{l = 1}^{m} \left(\hat{q}_{i}^{(l)} - \langle \hat{q}_i \rangle\right) \left(\hat{q}_{j}^{(l)} - \langle \hat{q}_j\rangle \right),$$ where $\hat{q}_i^{(l)}$ is the best-fitting value of the $i$-th spline parameter using the $l$-th mock catalog.

Observational data {#sec:data}
==================

As $q(z)$ is not a direct observable, we need to use other quantities to access $\{\hat{q}_i\}$. In this section, we present the samples of type Ia SNe, BAO and $H(z)$ measurements, and also their respective likelihood functions that we utilize to recover the deceleration function $q(z)$.

Type Ia supernova data
----------------------

We use the JLA sample [@Betoule2014] of 740 SNe Ia, whose likelihood is $$\label{eq:lik_snia}
-2\ln(L_{SNIa}) = \Delta{}\vec{m}^T\mathsf{C}_{SNIa}^{-1}\Delta{}\vec{m},$$ where the data covariance is a combination of the systematic and statistical errors $\mathsf{C}_{SNIa} = \mathsf{C}_{sys} + \mathsf{C}_{stat}(\alpha,\beta)$, and $$\begin{aligned}
\label{eq:snia_mb}
\Delta{}m_i &=&  m_{Bi} - m_{Bi}^{\text{th}}\\ 
&=& m_{Bi} -5\log_{10}(\mathcal{D}_L(z^{\text{hel}}_i, z^{{\text{cmb}}}_i)) + \alpha X_i - \beta \mathcal{C}_i - M_{h_i} + 5\log_{10}(c/H_0) - 25. \nonumber\end{aligned}$$ $m_{Bi}$ is the rest-frame peak B-band magnitude of the $i$-th SN Ia, and $z^{\text{hel}}_i$ and $z^{\text{cmb}}_i$ are its heliocentric and CMB frame redshits, respectively. The SN Ia astrophysical model contains four parameters $(\alpha, \beta, M_1, M_2)$, where the first two are related to the stretch-luminosity and colour-luminosity, respectively, and $M_1$ and $M_2$ are absolute magnitudes. The luminosity distance [@Davis2011] is $$\begin{aligned}
\label{eq:dis_SN}
D_L(z^{\text{hel}}, z^{\text{cmb}}) &= \frac{c}{H_0} \mathcal{D}_L(z^{\text{hel}}, z^{\text{cmb}}) \nonumber \\ 
&= (1 + z^{\text{hel}}) D_M(z^{\text{cmb}}),\end{aligned}$$ where the transverse comoving distance in a flat spatial sections universe is $D_M(z) = D_c(z)$.

Baryon acoustic oscillation
---------------------------

The peak position of the angular correlation function of the matter density can be measured by the distance ratio $D_V (z) / r_s(z_d)$ (see [@Thepsuriya2014; @Aubourg2014] and references therein). The volume-averaged-distance for perturbations along and orthogonal to the line of sight is defined as [@Eisenstein2005] $$\label{eq:bao_dv}
D_V(z) \equiv \left[ D_M(z)^2 \frac{cz}{H(z)}\right]^{1/3}.$$ The sound horizon $r_s(z_d)$ at the drag redshift $z_d$ (i.e., epoch at which baryons were released from photons) is $$\label{eq:bao_rd}
r_d \equiv r_s(z_d) = \frac{1}{H_0}\int_{z_d}^\infty dz \frac{c_s(z)}{E(z)},$$ where $c_s(z)$ is the sound wave speed in the photon-baryon fluid.

In this work we use 6 BAO data points as described in table \[tab:bao\] in appendix \[app:sampling\]. The first is measured by Beutler et al. [@Beutler2011] using galaxies from the 6dF Galaxy Survey. Padmanabhan et al. [@Padmanabhan2012] reported an improved data obtained with the reconstruction method [@Eisenstein2007] using Luminous Red Galaxy sample from SDSS Data Release 7 (DR7). Kazin et al. [@Kazin2014] give three points computing the power spectrum and correlation function of galaxies from the WiggleZ Survey in three correlated redshift bins. The last data is obtained by Ross et al. [@Ross2014] which used galaxies from SDSS DR7 with $z < 0.2$.

The BAO likelihood is $$\label{eq:lik_bao}
-2 \ln L_{\text{BAO}} = \left(\vec{b}^{\text{th}} - \vec{b} \right)^T\mathsf{C}_{\text{BAO}}^{-1}\left( \vec{b}^{\text{th}} -\vec{b} \right) - 2 \ln L_{\text{Ross}},$$ where $\vec{b}^{\text{th}}$ is the observable vector calculated using the theoretical model, i.e., the components are given by $b^{\text{th}}_i = D_V(z_i)/r_d$ calculated at each redshift (second column of table \[tab:bao\]). The vector $\vec{b}$ represent the observed version of these quantities and its components are provided by the third column of table \[tab:bao\]. The matrix $\mathsf{C}^{-1}_{\text{BAO}}$ is the inverse covariance matrix appearing in the BAO likelihood (see table \[tab:bao\]). Finally, Ross et al. [@Ross2014] pointed out that their data should be used considering their estimate of the likelihood distribution \[which we called $L_{\text{Ross}}$ in eq. \], since, in this case, the Gaussian distribution is not a good approximation.

References [@Beutler2011] and [@Padmanabhan2012] used the Eisenstein & Hu [@Eisenstein1998] (EH98) fitting function to compute $r_d^{{\text{fid}}}$ while [@Kazin2014; @Ross2014] used CAMB [@Lewis2000]. As mentioned in Ref. [@Kazin2014], the difference between $r^{{\text{fid}}}_{d,\text{EH98}}$ and $r^{{\text{fid}}}_{d, \text{CAMB}}$ is of order of $3\%$. So due to the current error magnitude of these data, this difference is relevant and, hence, we have to re-scale the data such that all measurements refer to the same method. In particular, we multiply Beutler and Padamanabhan’s data by $r_{d, \text{EH98}}^{{\text{fid}}} / r_{d, \text{CAMB}}^{{\text{fid}}} = 1.027$ and 1.025, respectively.

The BAO observable depends on the kinematic model through $D_V(z)$ and also requires the $r_d$ value. However, to calculate $r_d$ theoretically, we would need to extend the kinematic model to high redshifts and to compute the decoupling redshift $z_d$. We avoid this making $r_d$ a free parameter in the analysis. Therefore, throughout this work we fit $r_d$ along with the other parameters.

Hubble function
---------------

We work with 21 measurements of $H(z)$: 11 are provided by Stern et al. [@Stern2010] in the redshift range $0.1 \leq z \leq 1.75$ obtained from the spectra of red-enveloped galaxies; Riess et al. [@Riess2011] estimated the Hubble constant $H_0$ using optical and infrared observations of over 600 Cepheid variables; other 8 measurements are comprised within $0.1791 \leq z \leq 1.037$ as presented in Moresco et al. [@Moresco2012]. These last were obtained from the differential spectroscopic evolution of early-type galaxies with respect to $z$. In particular, we use their $H(z)$ values computed assuming the BC03 model for the differential evolution. The last datum is an estimate of $H(z)r_d / (1+z)$ at $z = 2.3$, and it is derived from the transmitted flux fraction in the Ly$\alpha$ forest of over 48,000 quasars combined with CMB observations as showed by Busca et al. [@Busca2013].

Assuming that $H(z)$ follows a Gaussian distribution and given that the error of each measurement is independent, we have that the likelihood is $$\label{eq:lik_hz}
-2 \ln L_{H} = \sum_{i = 1}^{20} \frac{\left(H(z_i) - H^{\text{obs}}_i\right)^2}{\sigma_i^2} + \frac{\left.\left[\frac{H(z)r_d}{1+z} - H_{r}^{\text{obs}}\right]^2\right\vert_{z=z_{21}}}{\sigma_{21}^2},$$ where $H_i^{\text{obs}}$ and $\sigma_i$ are the data points and their respective errors, $H(z) = H_0 E(z)$ and $E(z)$ is given by eq. . The observed values of the $20$ measures of the Hubble function are depicted in table \[tab:Hz\] in appendix \[app:sampling\]. The last observable $H^{\text{obs}}_r$ is described in the footnote of the same table.

Methodology and Validation {#sec:meth_valid}
==========================

The observational data set represent just one realization of their underlying data probability distributions. So by construction, we cannot access the bias of an estimated function, e.g., $\hat{q}(z)$, by fitting it using this data set. As we discussed in section \[sec:bv\_trade\], the $\hat{q}(z)$ bias is inferred knowing both the expected value and the true value of $q(z)$ \[see eq. \]. This last is the missing piece to compute the bias-variance trade-off. However, one can indirectly infer the bias without knowing the underlying $q(z)$.[^10] First, one fits the model obtaining the best fit for the real data. Then, using this best fit as the fiducial model, several mock catalogs are generated and, following the description in section \[sec:bv\_trade\], one calculates the bias and, consequently, the bias-variance trade-off.

The problem with this procedure is that one must choose the parameters of the reconstruction method to compute the best fit. In particular, we have to fix the number of knots and the complexity of the function $\hat{q}(z)$, through the penalization parameter ${\sigma_{\text{rel}}}$. Thus, if the initial best fit for a given ${\sigma_{\text{rel}}}$ is already significantly biased, so will be the bias-variance trade-off analysis. In other words, performing a MC analysis of the bias-variance trade-off around the best fit does not take into account the variance of the estimated curve. Therefore, a better approach consists in using not only the best fit curve but a set of curves inside some statistical significance, i.e., curves whose parameters are inside some confidence interval of the best fit. This means that the procedure should be capable of reconstructing not only the best fit curve but also every other curves inside some significance level.

There is a high computational cost to study the bias-variance trade-off for a given fiducial curve. It is necessary to resample from the model $m$ times and to find the best fit for each realization. The whole calculation must be performed for different values of ${\sigma_{\text{rel}}}$ until the best bias-variance trade-off is attained.

In this work, we address the problem by performing the MC analyses for three different fiducial curves and seven ${\sigma_{\text{rel}}}$ values (see sections \[sec:mc\_sneia\] and \[sec:mc\_alldata\]). The functional forms of these fiducial models were purposely defined to have quite different features as shown in figure \[fig:fiducial\_models\]. Given these distinct scenarios, we can explore the capability and efficiency of the proposed method in reconstructing $q(z)$ and also verify any dependence on the underlying model. Essentially, we want to find the ${\sigma_{\text{rel}}}$ value which best reconstructs all the three fiducial models.

We carry out this study considering: (i) only the SN Ia data and (ii) jointly the SN Ia, BAO and $H(z)$ measurements, in sections \[sec:mc\_sneia\] and \[sec:mc\_alldata\], respectively. Therefore, each realization is a (pseudo)randomly-generated catalog of SNe Ia $\{m_{B, a}, X_a, \mathcal{C}_a\}$ or random samples of all three observables, i.e., $\{\{m_{B,a}, X_a, \mathcal{C}_a\}, \{(D_V/r_d)_b\}, \{H_c\}\}$, where $a = 1,...,740$, $b = 1,...,6$ and $c = 1,...,21$. The methodology and algorithms to generate these samples are described in appendix \[app:sampling\].

We use the Monte Carlo object of the *Numerical Cosmology* Library (NumCosmo) [@DiasPintoVitenti2014], called . This object proceeds as describe in algorithm \[algo:NcmFitMC\]. The code generates a minimum of $\mathsf{prerun}$ mock samples and the catalog of the best fit for each sample.

In the list below, we summarize the steps to obtain $\hat{q}(z)$ using the MC method:

1.  Define the redshift domain $\mathrm{D} = [z_{min}, z_{max}]$, which is determined by the real data sample as described in appendix \[app:sampling\].

2.  Define the fiducial model $q^{\text{fid}}(z)$, where $z \in \mathrm{D}$.

3.  Choose a ${\sigma_{\text{rel}}}$ value and the number of knots $n$, which will be the same for all ${\sigma_{\text{rel}}}$ values to be tested.

4.  Run the MC algorithm for $m = \mathsf{prerun}$ minimum realizations.

5.  Finally, compute the MSE.

These steps are repeated for different values of ${\sigma_{\text{rel}}}$ and their respective results are compared in order to determine the best scenario which minimizes the MSE \[eq. \] for a given $m_b$ \[eq. \].

$l = 0$\

Monte Carlo analyses: SNe Ia {#sec:mc_sneia}
----------------------------

We carry out the first analysis using the SNe Ia data, such that the realizations are generated using the covariance matrix and redshifts of the JLA sample (see appendix \[app:sampling\]). Thus, $\mathrm{D} = [0.0, 1.3]$ which we divide in 7 equally spaced intervals, i.e., $n + 1 = 8$ knots.[^11] The first and second fiducial models, denoted as $q^{\text{fid1}}(z)$ and $q^{\text{fid2}}(z)$, are defined by the blue and green curves, respectively, in figure \[fig:fiducial\_models\]. The third one, $q^{\text{fid3}}(z)$, is the $\Lambda$CDM model (red curve), in which the cold dark matter and baryon density parameters are $\Omega_c = 0.3$ and $\Omega_b = 0.05$ and the DE EoS is $w = -1.0$. In all scenarios we assume a flat universe ($\Omega_k = 0$).[^12] Finally, these three fiducial models are completely defined by fixing the SN Ia astrophysical parameters, namely $\alpha^{\text{fid}}= 0.141$, $\beta^{\text{fid}}= 3.101$, $M_1^{\text{fid}}= -19.05$ and $M_2^{\text{fid}}= -19.12$. These values correspond to the best-fit obtained in [@Betoule2014] considering $\Lambda$CDM and SNe Ia data (including both systematic and statistical errors).[^13]

![The three $q^{\text{fid}}(z)$ fiducial models for which we study the reconstruction method via cubic spline. []{data-label="fig:fiducial_models"}](fiducial_models_new_z2.pdf)

![The top part of each panel shows the reconstructed curve of $q(z)$ using 8 knots, MC approach and sampling from SNe Ia data. The colored (blue, red and green) lines and shaded regions are the mean function $\langle\hat{q}(z)\rangle$ and their $1\sigma$ error bar, respectively, obtained for a given ${\sigma_{\text{rel}}}$ and $q^{\text{fid}}(z)$. The black lines correspond to $q^{{\text{fid}}1}$ (upper panel), $q^{{\text{fid}}2}$ (middle) and $\Lambda$CDM (lower). The bottom part of each panel shows the bias (dashed lines) and its $1\sigma$ error bar of the mean curve.[]{data-label="fig:qz_qs8"}](qz_plot_fiduc1_s8_jla.pdf "fig:") ![The top part of each panel shows the reconstructed curve of $q(z)$ using 8 knots, MC approach and sampling from SNe Ia data. The colored (blue, red and green) lines and shaded regions are the mean function $\langle\hat{q}(z)\rangle$ and their $1\sigma$ error bar, respectively, obtained for a given ${\sigma_{\text{rel}}}$ and $q^{\text{fid}}(z)$. The black lines correspond to $q^{{\text{fid}}1}$ (upper panel), $q^{{\text{fid}}2}$ (middle) and $\Lambda$CDM (lower). The bottom part of each panel shows the bias (dashed lines) and its $1\sigma$ error bar of the mean curve.[]{data-label="fig:qz_qs8"}](qz_plot_fiduc2_s8_jla.pdf "fig:") ![The top part of each panel shows the reconstructed curve of $q(z)$ using 8 knots, MC approach and sampling from SNe Ia data. The colored (blue, red and green) lines and shaded regions are the mean function $\langle\hat{q}(z)\rangle$ and their $1\sigma$ error bar, respectively, obtained for a given ${\sigma_{\text{rel}}}$ and $q^{\text{fid}}(z)$. The black lines correspond to $q^{{\text{fid}}1}$ (upper panel), $q^{{\text{fid}}2}$ (middle) and $\Lambda$CDM (lower). The bottom part of each panel shows the bias (dashed lines) and its $1\sigma$ error bar of the mean curve.[]{data-label="fig:qz_qs8"}](qz_plot_xcdm_s8_jla.pdf "fig:")

In order to obtain the features of the reconstruction procedure as a function of the penalty factor (through ${\sigma_{\text{rel}}}$) and $q^{\text{fid}}(z)$, we perform the MC analyses considering 7 different ${\sigma_{\text{rel}}}$ values, ${\sigma_{\text{rel}}}= \{5\%, 15\%, 30\%, 50\%, 75\%, 100\%, 150\% \},$ for each fiducial model. Independently of ${\sigma_{\text{rel}}}$ and $q^{\text{fid}}(z)$, we standardize our analyses fixing the number of realizations to $m= 42000$, with which we obtain small $\mathsf{lre} < 1\%$ (algorithm \[algo:NcmFitMC\]) in all cases. Then, for each mock catalog, we minimize the function $$\label{eq:f_snia}
-2 \ln (L_{\text{SNIa},P}) = -2 \ln (L_{\text{SNIa}}) + \sum_{i = 2}^{n-1} P_i({\sigma_{\text{rel}}}),$$ with respect to the free parameters $$\label{eq:set_p1}
\vec{\theta} \doteq \{\hat{q}_0, \, \hat{q}_1, \, \hat{q}_2, \, \hat{q}_3, \, \hat{q}_4, \, \hat{q}_5, \, \hat{q}_6, \, \hat{q}_7, \, \alpha, \, \beta, \, M_1, \, M_2\},$$ where the right-side terms of eq.  are given in eqs.  and , respectively. Therefore, the expected values $\langle \hat{q}(z; {\sigma_{\text{rel}}}) \rangle$ and the error bars $\sigma(\hat{q}(z; {\sigma_{\text{rel}}})) \equiv \sqrt{\text{Var}(\hat{q}(z; {\sigma_{\text{rel}}}))}$ of the reconstructed function $\hat{q}(z; {\sigma_{\text{rel}}})$ are estimated by computing eqs.  and .

Figure \[fig:qz\_qs8\] displays the reconstructed curves (colored solid lines) for ${\sigma_{\text{rel}}}= 5\%, 50\%$ and $150\%$, along with their respective fiducial models (black lines). As we discussed in section \[sec:penalty\], a small ${\sigma_{\text{rel}}}$ implies in a big constraint on the function complexity. Indeed, we see that ${\sigma_{\text{rel}}}= 5\%$ imposes $\hat{q}(z; 5\%)$ to be a linear function independently of the underlying fiducial model (upper-left part of the three panels).

Naturally, if the true model differs from a linear function, the result will be strongly biased and, consequently, the estimated curve will not be capable to recover the true form. In particular, $\hat{q}(z; 5\%)$ is outside the $1\sigma$ error bar in a large fraction of the redshift interval for $q^{{\text{fid}}1}(z)$ and $q^{{\text{fid}}2}(z)$. The bottom-left part of the three panels in figure \[fig:qz\_qs8\] show the bias (colored dashed lines) as a function of the redshift and its $1\sigma$ error bar.

Overall, until $z \lesssim 0.6$, the bias $b(z; {\sigma_{\text{rel}}})$ decreases as the $\hat{q}(z; {\sigma_{\text{rel}}})$ function complexity (i.e., ${\sigma_{\text{rel}}}$) increases, such that the reconstructed function approximates better and better the fiducial curve. On the other hand, due to the small amount of data at higher redshifts, mainly for $z \gtrsim 1.0$, the standard deviation $\sigma(\hat{q}(z; {\sigma_{\text{rel}}}))$ greatly increases, as evinced in figure \[fig:qz\_qs8\] by the colored shaded areas. As a direct result, we have that the constraints on the highest parameters $\hat{q}_i$ are degenerated causing an increment on $b(z; 150\%)$, for $z \gtrsim 0.6$, in comparison to $b(z; 50\%)$ for all three fiducial models, as shown in figure \[fig:qz\_qs8\].

![The top part of each panel shows the MSE of the reconstructed $\hat{q}(z)$ for 7 different ${\sigma_{\text{rel}}}$ values $\in [5\%, 150\%]$, and obtained using 8 knots, MC approach and sampling from SNe Ia data. The upper, middle and lower panels refer to $q^{{\text{fid}}1}(z)$, $q^{{\text{fid}}2}(z)$ and $\Lambda$CDM models, respectively. The MSE decomposition into variance (dotted lines) and squared bias (dashed lines) is displayed in the bottom part of each panel.[]{data-label="fig:mse_qs8"}](fiduc1_bias_var_qs8.pdf "fig:") ![The top part of each panel shows the MSE of the reconstructed $\hat{q}(z)$ for 7 different ${\sigma_{\text{rel}}}$ values $\in [5\%, 150\%]$, and obtained using 8 knots, MC approach and sampling from SNe Ia data. The upper, middle and lower panels refer to $q^{{\text{fid}}1}(z)$, $q^{{\text{fid}}2}(z)$ and $\Lambda$CDM models, respectively. The MSE decomposition into variance (dotted lines) and squared bias (dashed lines) is displayed in the bottom part of each panel.[]{data-label="fig:mse_qs8"}](fiduc2_bias_var_qs8.pdf "fig:") ![The top part of each panel shows the MSE of the reconstructed $\hat{q}(z)$ for 7 different ${\sigma_{\text{rel}}}$ values $\in [5\%, 150\%]$, and obtained using 8 knots, MC approach and sampling from SNe Ia data. The upper, middle and lower panels refer to $q^{{\text{fid}}1}(z)$, $q^{{\text{fid}}2}(z)$ and $\Lambda$CDM models, respectively. The MSE decomposition into variance (dotted lines) and squared bias (dashed lines) is displayed in the bottom part of each panel.[]{data-label="fig:mse_qs8"}](xcdm_bias_var_qs8.pdf "fig:")

In order to define which ${\sigma_{\text{rel}}}$ provides the best reconstructed curve, we compute the MSE \[eq. \] as a function of $z$ of all 21 reconstructed curves. Figure \[fig:mse\_qs8\] shows the MSE (solid lines) for each ${\sigma_{\text{rel}}}$ and $q^{{\text{fid}}1}(z)$ (upper panel), $q^{{\text{fid}}2}(z)$ (middle) and $\Lambda$CDM (lower) models. For these three fiducial models, $\hat{q}(z; 5\%)$ is the function with the smallest MSE for any $m_b$. However, these are also the most biased reconstructed curves and, even through visual inspection (figure \[fig:qz\_qs8\]), they do not give satisfactory reconstructions to their respective true $q^{\text{fid}}(z)$. Nevertheless, if we have only analyzed the fiducial curve closer to the best fit one (fiducial 3), we would wrongly conclude that ${\sigma_{\text{rel}}}= 5\%$ provides the smallest MSE with a insignificant bias (see the last panel in figure \[fig:mse\_qs8\]).

More importantly, when performing the analysis with real data, we will be able to estimate the error bars but not the bias. We can note in figure \[fig:mse\_qs8\] that the MSE for ${\sigma_{\text{rel}}}= 5\%$ has about the same contribution from bias and variance. Thus, if we choose the smallest MSE (${\sigma_{\text{rel}}}= 5\%$) disregarding $m_b$, the estimated error in the real data analysis would provide only half of the total uncertainty in the reconstruction. In a conservative approach one would estimate the variance and then double by hand the error bars.[^14] Notwithstanding, inspecting the ${\sigma_{\text{rel}}}= 5\%$ reconstructions in figure \[fig:qz\_qs8\], we note that this reconstruction looses all information about the shape of the curve. So even correcting the error bar for $\hat{q}(z)$ by doubling the computed error, for example, the form of the curve will always be close to a straight line.

The objective of this work is to reconstruct the form of the kinematic curve. We select the best reconstructed function (smallest MSE) requiring the bias to be at most $10\%$ of the total error, i.e., $m_b  = 0.1$. Making this imposition, the estimated variance for the fit using real data will provide a good approximation of the uncertainty of the reconstruction. Therefore, taking into account the three fiducial models, we find that the best bias-variance trade-off as a function of $z$ is achieved for ${\sigma_{\text{rel}}}= 30\%$. It is worth noting that the MC results for $\alpha$, $\beta$, $M_1$ and $M_2$ are unbiased (bias smaller than $0.1\%$) for all 21 cases that we studied.

Monte Carlo analyses: SNe Ia + BAO + $H(z)$ {#sec:mc_alldata}
-------------------------------------------

In this section, we perform the MC analyses combining the SNe Ia, $H(z)$ and BAO data. In this case, we equally divide the redshift interval, $\mathrm{D} = [0.0, 2.3]$, using $n + 1 = 12$ knots. The fiducial models and parameters are those defined in section \[sec:mc\_sneia\] and, additionally, $r_d^{\text{fid}}= 103.5 \, h^{-1}\text{Mpc}$ and $H_0^{\text{fid}}= 73.0 \, \text{km} \, \text{s}^{-1} \text{Mpc}^{-1}$. Similarly, the MC study is done for the same ${\sigma_{\text{rel}}}$ set and the three fiducial models. In view of the increased amount of data, the present $q(z)$ reconstruction is better constrained and $m= 20000$ mock catalogs (for each MC analysis) are sufficient to provide $\langle \hat{q}(z; {\sigma_{\text{rel}}}) \rangle$ with small $\mathsf{lre} < 1\%$ (see algorithm \[algo:NcmFitMC\]).

Since the SN Ia, BAO and $H(z)$ data sets are independent, the joint likelihood is $$-2\ln(L_{\text{SBH},P}) = -2 \left( \ln L_{\text{SNIa}} + \ln L_{\text{BAO}} + \ln L_{H} \right) + \sum_{i = 2}^{n-1} P_i({\sigma_{\text{rel}}}),$$ where $-2 \ln L_{SNIa}$, $-2 \ln L_{BAO}$ and $-2 \ln L_{H(z)}$ are given by equations , and , respectively. Then, for each realization, we compute the best-fitting values of the following 18 parameters, $$\label{eq:set_p2}
\vec{\theta} \doteq \{\hat{q}_0, \, \hat{q}_1, \, \hat{q}_2, \, \hat{q}_3, \, \hat{q}_4, \, \hat{q}_5, \, \hat{q}_6, \, \hat{q}_7, \, \hat{q}_8, \, \hat{q}_9, \, \hat{q}_{10}, \, \hat{q}_{11}, \, \alpha, \, \beta, \, M_1, \, M_2, H_0, r_d\},$$ and, with them, we calculate the expected value of each parameter estimator, $\langle \vec{\theta} \rangle$, and their covariance matrix.

![The top part of each panel shows the reconstructed curve of $q(z)$ using 12 knots, MC approach and sampling from SNe Ia, BAO and $H(z)$ data. The colored (blue, red and green) lines and shadow regions are the mean function $\langle\hat{q}(z)\rangle$ and their $1\sigma$ error bar, respectively, obtained for a given ${\sigma_{\text{rel}}}$ and $q^{\text{fid}}(z)$. The black lines correspond to $q^{{\text{fid}}1}$ (upper panel), $q^{{\text{fid}}2}$ (middle) and $\Lambda$CDM (lower). The bottom part of each panel shows the bias (dashed lines) and its $1\sigma$ error bar of the mean curve.[]{data-label="fig:qz_qs12"}](qz_plot_fiduc1_s12_jla.pdf "fig:") ![The top part of each panel shows the reconstructed curve of $q(z)$ using 12 knots, MC approach and sampling from SNe Ia, BAO and $H(z)$ data. The colored (blue, red and green) lines and shadow regions are the mean function $\langle\hat{q}(z)\rangle$ and their $1\sigma$ error bar, respectively, obtained for a given ${\sigma_{\text{rel}}}$ and $q^{\text{fid}}(z)$. The black lines correspond to $q^{{\text{fid}}1}$ (upper panel), $q^{{\text{fid}}2}$ (middle) and $\Lambda$CDM (lower). The bottom part of each panel shows the bias (dashed lines) and its $1\sigma$ error bar of the mean curve.[]{data-label="fig:qz_qs12"}](qz_plot_fiduc2_s12_jla.pdf "fig:") ![The top part of each panel shows the reconstructed curve of $q(z)$ using 12 knots, MC approach and sampling from SNe Ia, BAO and $H(z)$ data. The colored (blue, red and green) lines and shadow regions are the mean function $\langle\hat{q}(z)\rangle$ and their $1\sigma$ error bar, respectively, obtained for a given ${\sigma_{\text{rel}}}$ and $q^{\text{fid}}(z)$. The black lines correspond to $q^{{\text{fid}}1}$ (upper panel), $q^{{\text{fid}}2}$ (middle) and $\Lambda$CDM (lower). The bottom part of each panel shows the bias (dashed lines) and its $1\sigma$ error bar of the mean curve.[]{data-label="fig:qz_qs12"}](qz_plot_xcdm_s12_jla.pdf "fig:")

![The top part of each panel shows the MSE of the reconstructed $\hat{q}(z)$ for 7 different ${\sigma_{\text{rel}}}$ values $\in [5\%, 150\%]$, and obtained using 12 knots, MC approach and sampling from SNe Ia, BAO and $H(z)$ data. The upper, middle and lower panels refer to $q^{{\text{fid}}1}(z)$, $q^{{\text{fid}}2}(z)$ and $\Lambda$CDM models, respectively. The MSE decomposition into variance (dotted lines) and squared bias (dashed lines) is displayed in the bottom part of each panel.[]{data-label="fig:mse_qs12"}](fiduc1_bias_var_qs12.pdf "fig:") ![The top part of each panel shows the MSE of the reconstructed $\hat{q}(z)$ for 7 different ${\sigma_{\text{rel}}}$ values $\in [5\%, 150\%]$, and obtained using 12 knots, MC approach and sampling from SNe Ia, BAO and $H(z)$ data. The upper, middle and lower panels refer to $q^{{\text{fid}}1}(z)$, $q^{{\text{fid}}2}(z)$ and $\Lambda$CDM models, respectively. The MSE decomposition into variance (dotted lines) and squared bias (dashed lines) is displayed in the bottom part of each panel.[]{data-label="fig:mse_qs12"}](fiduc2_bias_var_qs12.pdf "fig:") ![The top part of each panel shows the MSE of the reconstructed $\hat{q}(z)$ for 7 different ${\sigma_{\text{rel}}}$ values $\in [5\%, 150\%]$, and obtained using 12 knots, MC approach and sampling from SNe Ia, BAO and $H(z)$ data. The upper, middle and lower panels refer to $q^{{\text{fid}}1}(z)$, $q^{{\text{fid}}2}(z)$ and $\Lambda$CDM models, respectively. The MSE decomposition into variance (dotted lines) and squared bias (dashed lines) is displayed in the bottom part of each panel.[]{data-label="fig:mse_qs12"}](xcdm_bias_var_qs12.pdf "fig:")

Analogously to the results presented in section \[sec:mc\_sneia\], in figure \[fig:qz\_qs12\] we show the reconstructed curves, $\langle \hat{q} (z; {\sigma_{\text{rel}}}) \rangle$ (colored solid lines), the biases (dashed lines) and their respective $1\sigma$ error bars (shaded areas) for ${\sigma_{\text{rel}}}= 5\%, 50\%$ and $150\%$ and the three fiducial models. As expected, combining SN Ia, BAO and $H(z)$ data decreases $\sigma(\hat{q}(z; {\sigma_{\text{rel}}}))$ and we note, in the three upper-left panels of figure \[fig:qz\_qs12\], that these improved constraints no longer restrict $\langle \hat{q}(z; 5\%) \rangle$ to be a linear function.[^15] Besides, as displayed in figure \[fig:mse\_qs12\], it also leads to a stronger dependency of $b(z; {\sigma_{\text{rel}}})$ on the fiducial model.

For example, the squared bias function of $q^{{\text{fid}}2}(z)$ (dashed lines on the middle panel of figure \[fig:mse\_qs12\]) presents a large variation with respect to ${\sigma_{\text{rel}}}$, similar to what we obtained considering only SN Ia data. On the other hand, the bias function of the $\Lambda$CDM fiducial model is nearly invariant with respect to ${\sigma_{\text{rel}}}$ (lower panel of figure \[fig:mse\_qs12\]). This highlights that, giving a sufficient amount of data, if the true model has low complexity form then, naturally, a satisfactory reconstruction will rapidly be reached even for small ${\sigma_{\text{rel}}}$. In particular, taking into account the bias-variance trade-off and $m_b = 0.1$, the best reconstructed curve for the $\Lambda$CDM fiducial model is obtained for ${\sigma_{\text{rel}}}= 15\%$.

As it was shown in figures \[fig:mse\_qs8\] and \[fig:mse\_qs12\], the MSE is dominated by the variance of the fitted parameters and, consequently, by the measurement errors of the observable data sets. In this way, we note that the MSE is independent of the fiducial model, depending only on the penalty factor, i.e., ${\sigma_{\text{rel}}}$. This aspect points out the fact that an statistical inference with relative small error bars can be completely misleading. For example, the well-constrained function $\hat{q}(z; 5\%)$, obtained for $q^{{\text{fid}}2}(z)$, is not at all a good reconstruction of the true model, diverging with more than $1\sigma$ in a large fraction of the redshift interval. Ultimately, following the discussion on section \[sec:PM\], the nature of the biases for ${\sigma_{\text{rel}}}= 5\%$ and $150\%$ can be classified as form-bias and estimator-bias, respectively.

As in section \[sec:mc\_sneia\], we expect that a good reconstruction is one which does not present a significant bias. Therefore, applying the same requirement $m_b = 0.1$, we obtain that ${\sigma_{\text{rel}}}= 30\%$ provides the best balance between bias and variance, given the current data sets, to reconstruct the $q(z)$ curve.

In summary, the MC outcomes (sections \[sec:mc\_sneia\] and \[sec:mc\_alldata\]) show that our method is efficient in reconstructing $q(z)$ and, given the current errors of the observational data, ${\sigma_{\text{rel}}}= 30\%$ is a safe and conservative choice to recover $q(z)$ imposing minimal assumptions and guaranteeing that the reconstructed curve will not be bias dominated.

Results {#sec:real_reconst}
=======

Defined the best estimator, $\hat{q}(z; 30\%)$, we now obtain the deceleration function given the real JLA, BAO and $H(z)$ samples (see section \[sec:data\]). As before, we reconstruct $q(z)$ (i) in the redshift interval $\mathrm{D} = [0, 1.3]$ using JLA SN Ia data and (ii) in $\mathrm{D} = [0, 2.3]$ combining those three observable data.

The likelihood function is now interpreted as the posterior distribution $$\mathcal{P}(\vec{\theta} | \vec{D}) = L_P(\vec{D}| \vec{\theta}),$$ since we are assuming flat priors on all parameters. We then use the NumCosmo algorithm ,[^16] which implements an ensemble sampler with affine invariance for Markov Chain Monte Carlo analysis, to compute the mean function $\overline{q}(z)$ given the JLA sample and its $68.27\%$, $95.45\%$ and $99.73\%$ confidence intervals (CI). We ran 100 chains, computing $8\times 10^5$ points in the 12-dimensional parametric space, shown in eq. , attaining a multivariate potential scale reduction factor (MPSRF) equal to $1.016$.[^17] The best-fitting values and the error bars of the SNe Ia parameters are $\alpha =0.141 \pm 0.007$, $\beta = 3.108 \pm 0.081$, $M_1 = -19.05 \pm 0.03$ and $M_2 = -19.12 \pm 0.03$.

![The model-independent reconstructed $q(z)$ (blue solid line) and its $68.23\%$, $95.45\%$ and $99.73\%$ confidence intervals (blue shaded areas), for ${\sigma_{\text{rel}}}= 30\%$, using SNe Ia (left panel) and SNe Ia + BAO + $H(z)$ data (right panel). The red line and contours are the $q(z)$ mean and its CI’s obtained, assuming XCDM model, fitting $\Omega_c$ (with $w = -1$) (left panel) and $(H_0, \Omega_c, w)$ (right panel) along with the SNe Ia nuisance parameters. The black lines correspond to the Planck+BAO+JLA+$H_0$ best-fit assuming $\Lambda$CDM model.[]{data-label="fig:real_qz"}](qz_plot_real_data_s8_planckbf.pdf "fig:") ![The model-independent reconstructed $q(z)$ (blue solid line) and its $68.23\%$, $95.45\%$ and $99.73\%$ confidence intervals (blue shaded areas), for ${\sigma_{\text{rel}}}= 30\%$, using SNe Ia (left panel) and SNe Ia + BAO + $H(z)$ data (right panel). The red line and contours are the $q(z)$ mean and its CI’s obtained, assuming XCDM model, fitting $\Omega_c$ (with $w = -1$) (left panel) and $(H_0, \Omega_c, w)$ (right panel) along with the SNe Ia nuisance parameters. The black lines correspond to the Planck+BAO+JLA+$H_0$ best-fit assuming $\Lambda$CDM model.[]{data-label="fig:real_qz"}](qz_plot_real_data_s12_planck.pdf "fig:")

The mean (blue line) and CI’s (blue shaded areas) are displayed in the left panel of figure \[fig:real\_qz\], where we note that the ESMCMC $68.27\%$ CI is in agreement to the respective MC error bar. Naturally, as in the MC study, the difficulty in constraining $q(z)$ with minimal assumptions is that the result is highly degenerated. So, besides the fact that $\overline{q}(z) < 0$ in the entire redshift, this is not statistically significant. The main results of this $q(z)$ reconstruction, using only SNe Ia data, are: the indication of a transition redshift at $z_{\tiny{T}} \simeq 0.53$ with $68.27\%$ significance level, and the evidence of an accelerated expansion phase $\geq 99.73\%$ within the redshift interval $\sim [0.04, 0.19]$.

To compare this result with the flat $\Lambda$CDM model, i.e., assuming GR and DE EoS given by $w = -1.0$, we carry out the ESMCMC in the 5-dimensional parametric space $$\vec{\theta} \doteq (\Omega_c, \alpha, \beta, M_1, M_2).$$ For this, we fixed the other cosmological parameters to the JLA best-fit. We obtain the following best-fitting values and standard deviations: $\Omega_c = 0.244 \pm 0.034$, $\alpha = 0.141 \pm 0.007$, $\beta = 3.103 \pm 0.081$, $M_1 = -19.05 \pm 0.023$ and $M_2 = -19.12 \pm 0.026$. The mean $\overline{q}^{\text{XCDM}} (z)$ (red line) and the CI’s (red shaded areas) are consistent with our model-independent reconstruction overall redshift range, as shown in the left panel of figure \[fig:real\_qz\]. In both analyses, the SN Ia nuisance parameters $\alpha$ and $\beta$ are weakly (anti-) correlated ($\lesssim 0.1$) to $\{\hat{q}_i\}$ $(i = 0,..., 7)$ and $\Omega_c$. $M_1$ and $M_2$ are also weakly correlated to most parameters excepted for $q_0$ and $\Omega_c$, in which there are moderate correlations $\sim 0.56 - 0.66$. At last, we also plot the Planck best-fitting curve (black line), $\overline{q}^{\text{Planck}}(z)$, which we computed using the Planck+BAO+JLA+$H_0$ best-fitting parameters obtained assuming $\Lambda$CDM [@Planck2015], namely, $H_0 = 67.74$, $\Omega_c = 0.259$ and $\Omega_b =0.049$. We note that $\overline{q}^{\text{Planck}}(z)$ is inside the $68.27\%$ CI in the entire redshift interval.

We follow the same procedure to compute $\overline{q}(z)$ and their CI in the redshift interval $[0.0, 2.3]$ given the JLA, BAO and $H(z)$ data. In this case, $8\times 10^5$ points are calculated in the 18-dimensional parametric space \[eq. \], getting $\text{MPSRF} = 1.033$. In this case, the best-fitting and standard deviation values of the non-$q_i$ parameters are $H_0 = 71.68 \pm 1.69 \, \text{km} \, \text{s}^{-1} \text{Mpc}^{-1}$, $r_d = 101.15 \pm 1.8 \, h^{-1}\text{Mpc}$, $\alpha = 0.141 \pm 0.007$, $\beta = 3.111 \pm 0.081$, $M_1 = -19.01 \pm 0.05$ and $M_2 = -19.07 \pm 0.05$. We note in the right panel of figure \[fig:real\_qz\] an expressive improvement on the constraints due to the combined data, regarding that we are assuming only FLRW metric and flat space, and given the high dimension of the parametric space. In this case, we obtain $z_{\tiny{T}} \simeq 0.58$ with $68.27\%$ CI and we can attest with significance $\geq 99.73\%$ that the universe is accelerating within $0.02 \lesssim z \lesssim 0.22$. At last, we also have an indication of a decelerated phase in the redshift interval $1.84 \lesssim z \lesssim 2.13$ with $68.27\%$.

We compare the model-independent $q(z)$ reconstruction with the XCDM model, i.e., constant EoS for the dark energy component $w = constant$, performing the ESMCMC on the 7-dimensional parametric space $$\vec{\theta} \doteq (\Omega_c, H_0, w, \alpha, \beta, M_1, M_2).$$ Their best-fitting values and error bars are $\Omega_c = 0.269 \pm 0.017$, $H_0 = 71.04 \pm 1.59 \, \text{km} \, \text{s}^{-1} \text{Mpc}^{-1}$, $w = -1.11 \pm 0.07$, $\alpha = 0.142 \pm 0.007$, $\beta = 3.108 \pm 0.081$, $M_1 = -19.03 \pm 0.050$ and $M_2 = -19.10 \pm 0.05$. In figure \[fig:real\_qz\] (right panel), we plotted the results, $\overline{q}^{\text{XCDM}} (z)$ (red line) and the CI’s (red shaded areas), showing that XCDM model is also compatible within $99.73\%$ confidence level with the reconstructed curve in the entire redshift interval. It is worth noting that our conservative analyses provide only lower bounds for $z_{\tiny{T}} \simeq 0.19$ (JLA) and $\simeq 0.22$ \[JLA+BAO+$H(z)$\] within $99.73\%$. On the other hand, model-dependent works can estimate better constrained intervals for $z_{\tiny{T}}$ [@Capozziello2014; @Santos2015], as evinced by the XCDM results. However, as we discussed in sections \[sec:introduction\] and \[sec:review\], these estimates are accurate as long as the form-bias is small.

Besides the results presented so far, we can also infer other kinematic quantities from the reconstructed function $q(z)$ (see sections \[sec:introduction\] and \[sec:def\_q\]), such as $H(z)$ and the jerk function $j(z) \equiv - \dddot{a}/(aH^3)$ (i.e., the third order term of the scale factor expansion). The functions $H(z)$ and $j(z)$ are estimated by integrating and differentiating, respectively, $q(z)$ with respect to z. Figure \[fig:real\_Hz\_jz\] shows the results obtained using only JLA[^18] (left panels) and JLA+BAO+$H(z)$ (right panels). The upper panels display the means $\overline{H}(z)$ (blue lines) and the $68.23\%$, $95.45\%$ and $99.73\%$ CI’s (blue shaded areas).

As the $H(z)$ estimate is related with $q(z)$ by a numerical integration, it is consequently better constrained and the enhancement due to the combined data is even more apparent, as shown figure \[fig:real\_Hz\_jz\] (upper right panel), where the 21 $H(z)$ measurements are also plotted. As before, we compare these results with the respective outcomes assuming, respectively, $\Lambda$CDM and XCDM models (red lines and shaded areas). The decrease in the CI’s of the reconstructed $H(z)$ highlights its concordance with those cosmological models within $99.73\%$ level. The Planck best-fit $\overline{H}^{\text{Planck}}(z)$ (black lines) also lies inside the CI’s over the entire redshift interval for JLA+BAO+H(z), and it is outside the $99.73\% $ CI only in $0.0 \leq z \lesssim 0.03$ when we use JLA.

Contrarily to the above scenario, $j(z)$ is obtained by numerical differentiation implying in high degenerated results, as displayed in the lower panels of figure \[fig:real\_Hz\_jz\]. The means $\overline{j}(z)$ and the CI’s are represented by the blue lines and blue shaded areas, respectively. The estimated jerk functions, obtained assuming $\Lambda$CDM (left panel) and XCDM models (right panel), correspond to the red lines and red shaded areas. Besides the redshift intervals $[0.14, 0.39]$ (left panel) and $[0.13, 0.45]$ (right panel), which provides $-6 \lesssim j(z) \lesssim 6$ with $\geq 99.73\%$ significance, the degenerated results do not provide relevant information about this kinematic function.

![$H(z)$ (upper panels) and $j(z)$ (lower panels) mean estimates (blue solid lines) and their $68.23\%$, $95.45\%$ and $99.73\%$ confidence intervals (blue shaded areas) obtained integrating and deriving, respectively, the reconstructed $q(z)$ using SNe Ia (left panels) and SNe Ia + BAO + $H(z)$ data (right panels). Similarly, the red line and contours are the $H(z)$ and $j(z)$ means and their CI’s obtained integrating/deriving the respective $q(z)$ where the $\Lambda$CDM (left panels) and XCDM (right panels) models were assumed. The black lines correspond to the Planck+BAO+JLA+$H_0$ best-fit assuming $\Lambda$CDM model. The black dots and error bars represent the 21 $H(z)$ measurements used in this work.[]{data-label="fig:real_Hz_jz"}](Hz_plot_real_data_s8_planckbf.pdf "fig:") ![$H(z)$ (upper panels) and $j(z)$ (lower panels) mean estimates (blue solid lines) and their $68.23\%$, $95.45\%$ and $99.73\%$ confidence intervals (blue shaded areas) obtained integrating and deriving, respectively, the reconstructed $q(z)$ using SNe Ia (left panels) and SNe Ia + BAO + $H(z)$ data (right panels). Similarly, the red line and contours are the $H(z)$ and $j(z)$ means and their CI’s obtained integrating/deriving the respective $q(z)$ where the $\Lambda$CDM (left panels) and XCDM (right panels) models were assumed. The black lines correspond to the Planck+BAO+JLA+$H_0$ best-fit assuming $\Lambda$CDM model. The black dots and error bars represent the 21 $H(z)$ measurements used in this work.[]{data-label="fig:real_Hz_jz"}](Hz_plot_real_data_s12_planck.pdf "fig:") ![$H(z)$ (upper panels) and $j(z)$ (lower panels) mean estimates (blue solid lines) and their $68.23\%$, $95.45\%$ and $99.73\%$ confidence intervals (blue shaded areas) obtained integrating and deriving, respectively, the reconstructed $q(z)$ using SNe Ia (left panels) and SNe Ia + BAO + $H(z)$ data (right panels). Similarly, the red line and contours are the $H(z)$ and $j(z)$ means and their CI’s obtained integrating/deriving the respective $q(z)$ where the $\Lambda$CDM (left panels) and XCDM (right panels) models were assumed. The black lines correspond to the Planck+BAO+JLA+$H_0$ best-fit assuming $\Lambda$CDM model. The black dots and error bars represent the 21 $H(z)$ measurements used in this work.[]{data-label="fig:real_Hz_jz"}](jz_plot_real_data_s8.pdf "fig:") ![$H(z)$ (upper panels) and $j(z)$ (lower panels) mean estimates (blue solid lines) and their $68.23\%$, $95.45\%$ and $99.73\%$ confidence intervals (blue shaded areas) obtained integrating and deriving, respectively, the reconstructed $q(z)$ using SNe Ia (left panels) and SNe Ia + BAO + $H(z)$ data (right panels). Similarly, the red line and contours are the $H(z)$ and $j(z)$ means and their CI’s obtained integrating/deriving the respective $q(z)$ where the $\Lambda$CDM (left panels) and XCDM (right panels) models were assumed. The black lines correspond to the Planck+BAO+JLA+$H_0$ best-fit assuming $\Lambda$CDM model. The black dots and error bars represent the 21 $H(z)$ measurements used in this work.[]{data-label="fig:real_Hz_jz"}](jz_plot_real_data_s12.pdf "fig:")

At last, in this work we also obtain a model-independent estimate of the sound horizon, since $r_d$ is fitted throughout the reconstruction procedure. Thus, we obtain $r_d = 101.15 \pm 1.8 \, h^{-1}\text{Mpc}$ given the combined analysis in the 18-dimensional parametric space using JLA+BAO+$H(z)$ data. Our result is in accordance with that presented in reference [@Heavens2014], where they obtained $r_d = 100.7 \pm 2.0 \, h^{-1}\text{Mpc}$ (SNe+BAO+$H(z)$ without Hubble prior) and $r_d = 101.9 \pm 1.9 \, h^{-1}\text{Mpc}$ (idem with Hubble prior) assuming a linear spline for $h^{-1}(z) = 100 \, \text{km}\,\text{s}^{-1}\text{Mpc}^{-1}/ H(z)$, 6 knots equally-spaced and curved universe. It is worth noting that their reconstruction was directly applied to the data and no test regarding the assumed curve was performed.

Conclusions {#sec:conclusions}
===========

In this work we presented a general model-independent approach to reconstruct any one-variable function. In particular, we applied it to estimate the deceleration function $q(z)$ using the JLA SN Ia, BAO and $H(z)$ data sets. We performed a conservative analysis in the sense that we reconstructed $q(z)$ considering minimal assumptions (FLRW metric and flat universe) and also fitting simultaneously $H_0$, $r_d$ and the SN Ia astrophysical parameters ($\alpha$, $\beta$, $M_1$, $M_2$).

The MC results presented in sections \[sec:mc\_sneia\] and \[sec:mc\_alldata\] show that it is crucial to estimate both the form- and the estimator-biases in order to validate any proposed reconstruction method, independently if it is parametric or non-parametric. They also reveal that the MSE’s are independent of the underlying fiducial model and are dominated by the variance of the current data sets, even for cases where the reconstructed functions are far to be a good representation of the true curves, i.e., highly biased. Nevertheless, a blindly minimization of the MSE could lead to a bias dominated reconstruction. Such reconstructions must be carefully analyzed and the bias added by hand or estimated in some way. Even so, the information about the shape of the curve is mostly lost. Thus, requiring the bias to be at most $10\%$ of the total MSE, we evaluated the bias-variance trade-off obtaining that, currently, the best penalty factor is given by ${\sigma_{\text{rel}}}= 30\%$.

Even though the MSE is approximately independent of the fiducial model, analyzing MSE of only one model can lead to misleading conclusions. If we had only used the $\Lambda$CDM fiducial model, the smaller ${\sigma_{\text{rel}}}= 5\%$ would prove to be the best choice, i.e., smaller MSE and insignificant bias. In other words, if a method is good in reconstructing a $\Lambda$CDM shaped curve, this does not mean that the same method will be able to reconstruct other curves close to it. In addition, capping the bias at $10\%$ is a straightforward way to obtain the total uncertainty of the estimates, since in this case we can use directly the variance analysis without requiring additional bias estimation.

Our main reconstructions are summarized in figure \[fig:real\_qz\], given that there are no assumptions on the gravitational dynamics nor on the matter content, and we also include the relevant astrophysical parameters into the study. The standard cosmological model agrees with the reconstructed results within $99.73\%$ CI in the entire redshift intervals. This is even further evident observing the estimate for the Hubble function in figure \[fig:real\_Hz\_jz\]. In this sense this work discards large deviations from the standard model with the caveat that we are assuming the FLRW metric and we are not testing this hypothesis. Notwithstanding, with more data available, we advocate that to measure departures from the standard model, one should repeat the process described in this work, evaluating the capability of the model to differentiate between a whole class of fiducial models.

SDPV thanks CAPES for the grant 2649-13-6. MPL thanks CNPq (PCI/MCTI/INPE program and grant 202131/2014-9) for financial support. This research was performed using the Mesu-UV supercomputer of the Pierre & Marie Curie University – France (UPMC) and the computer cluster of the State University of the Rio Grande do Norte (UERN) – Brazil. We also thank Pierre Astier and Marc Betoule for kindly providing the complete covariance matrix necessary for the SNe Ia analysis, and Nicolás Busca for useful comments.

Simulated data {#app:sampling}
==============

In this appendix we describe the methodology to generate the SN Ia, BAO and $H(z)$ mock catalogs which we used to perform the MC study in section \[sec:meth\_valid\].

SN Ia sampling
--------------

The SN Ia distance modulus $\mu$ is written in terms of its light-curve parameters $(m_B, X, \mathcal{C})$ as $$\label{eq:dist_mod}
\mu = m_B - (M_{h_i} - \alpha X + \beta \mathcal{C}),$$ where $m_B$ is the observed peak magnitude in rest-frame B, $X$ and $\mathcal{C}$ are the observed stretching and the color at maximum brightness, respectively, and $\alpha$, $\beta$ and $M_{h_i} = \{M_1, M_2\}$ are nuisance parameters [@Betoule2014]. Thus, in order to generate SNe Ia samples, we must provide fiducial values of $(m_{B}, X, \mathcal{C})$, for each SN Ia, and their respective probability distributions. The fiducial values are usually determined choosing a specific set of the observable model parameters. In this case, from eq.  and the theoretical distance modulus, $\mu^{\text{th}} = 5\log_{10}(\mathcal{D}_L(z^{\text{hel}}, z^{{\text{cmb}}}) H_0/ c) + 25$, we can write the fiducial magnitude as function of $X_i$ and $\mathcal{C}_i$, $$\label{eq:mb_sneia}
m_{B_i}^{\text{fid}}= 5\log_{10}\left[\mathcal{D}_L\left(z^{\text{hel}}_i, z^{{\text{cmb}}}_i, \hat{q}^{\text{fid}}\right)\right] - \alpha^{\text{fid}}X_i + \beta^{\text{fid}}\mathcal{C}_i + M_{h_i}^{\text{fid}}- 5\log_{10}(c/H_0^{\text{fid}}) + 25,$$ where $i$ denotes the SN Ia index and the SN Ia redshift values $\{z^{\text{hel}}_i, z^{{\text{cmb}}}_i\}$ are always kept fixed to their observable estimates. As described in section \[sec:mc\_sneia\], the fiducial parameters are $$\{\alpha^{\text{fid}}, \beta^{\text{fid}}, M_1^{\text{fid}}, M_2^{\text{fid}}, H_0^{\text{fid}}\} = \{0.141, 3.101, -19.05, -19.12, 73.0\},$$ and they correspond (besides $H_0^{\text{fid}}$) to the best-fitting values assuming $\Lambda$CDM model obtained in reference [@Betoule2014]. We consider the three different $q^{\text{fid}}(z)$ represented in figure \[fig:fiducial\_models\]. At this point, we still have to define the fiducial values for $X_i$ and $\mathcal{C}_i$.

However, there is no astrophysical models for $X$ and $\mathcal{C}$. Therefore, we define the SNe Ia fiducial model $(m_{B_i}^{\text{fid}}, X_i^{\text{fid}}, \mathcal{C}_i^{\text{fid}})$ as being the best-fitting values obtained by maximizing, with respect to all $X_i$ and $\mathcal{C}_i$, the multivariate Gaussian distribution [@Betoule2014] $$\label{eq:sn_all_dist}
G(\vec{\zeta}, \mathsf{C}_{cmp}) = \frac{1}{\sqrt{(2\pi)^{3N} \vert\mathsf{C}_{cmp}\vert}} e^{-\frac{1}{2}\vec{\zeta}^T \mathsf{C}_{cmp}^{-1} \vec{\zeta}},$$ where $\vert ... \vert$ denotes the determinant, $\mathsf{C}_{cmp}^{-1}$ is the complete inverse covariance matrix (with 2220 rows and columns), $N$ is the number of SNe Ia of the JLA sample (740) and $$\vec{\zeta} = \left(m_{B_1}, ..., m_{B_N}, X_{N+1}, ..., X_{2N}, \mathcal{C}_{2N+1}, ..., \mathcal{C}_{3N}\right),$$ with $\{m_{B_i}\}$ given by eq. .[^19]

Finally, writing the SNe Ia variables as $$m_{B_i} = m_{B_i}^{\text{fid}}+ \delta m_{B_i}, \quad X_i = X^{\text{fid}}_i + \delta X_i \quad \text{and} \quad \mathcal{C}_i = \mathcal{C}^{\text{fid}}_i + \delta\mathcal{C}_i,$$ we create a SN Ia sample $\{m_{B_i}, X_i, \mathcal{C}_i\}$ by randomly generating 2220 values from the normal distribution with variance $\mathsf{C}_{cmp}$ and zero mean corresponding to $$\overrightarrow{\delta\zeta} = \left(\delta m_{B_1}, ..., \delta m_{B_N}, \delta X_{N+1}, ..., \delta X_{2N}, \delta\mathcal{C}_{2N+1}, ..., \delta\mathcal{C}_{3N}\right).$$

  -------------------------- ------- -------------- -------------- ------- --------- --------- --------- --- --
          Reference             z     $D_V(z)/r_d$                                                           
   \[0.5ex\] Beutler et al.   0.106     $0.336$      $0.015^{-2}$        0         0         0         0   0 
      Padmanabhan et al.      0.35        8.88            0          34.60         0         0         0   0 
                              0.44       11.550           0              0    4.8116   -2.4651    1.0375   0 
                              0.60       14.945           0              0   -2.4651    3.7697   -1.5865   0 
                              0.73       16.932           0              0    1.0375   -1.5865    3.6498   0 
         Ross et al.          0.15       4.466            -              -         -         -         -   - 
           \[1ex\]                                                                                           
  -------------------------- ------- -------------- -------------- ------- --------- --------- --------- --- --

  : BAO data[]{data-label="tab:bao"}

In Ref. [@Beutler2011] the authors provide $r_d / D_V(z)$. Note that it is the reciprocal of the provided in the other references. As this data is not correlated with the others, we can include it directly in $L_{BAO}$ \[eq. \].

Refs. [@Kazin2014] and [@Ross2014] provide $D_V(z)(r_d^{{\text{fid}}}/r_d)$, where $r_d^{{\text{fid}}} = 148.6$ and $148.69 \, h^{-1}\text{Mpc}$, respectively. We consider these values to build $L_{\text{BAO}}$.

$H(z)$ and BAO sampling
-----------------------

Differently from the SN Ia sampling, the BAO and $H(z)$ mock catalogs can be directly generated, since we have theoretical models for the respective observables, i.e., $D_V(z)/r_d$ \[or $r_d/D_V(z)$\], $H(z)$ and $H(z)r_d /(1+z)$, as displayed in eqs. , and .

Each $H^\text{obs}$ catalog correspond to $\{H_i^\text{obs}\}$, where $i = 1,..., 20$, and the point $\{H_r^{\text{obs}}\}$. The 21 redshift values are not altered in the sampling procedure and they are equal to the observed data $z$ of table \[tab:Hz\]. We obtain $H_i^\text{obs}$ writing it as $$H_i^\text{obs} = H^{\text{fid}}(z_i) + \delta H_i,$$ and, then, we randomly generate the quantity $\delta H_i$ from a Gaussian distribution with mean 0 and standard deviation $\sigma_i$ (last column of table \[tab:Hz\]), since $H^\text{obs}_i$ follows a Gaussian distribution. Analogously, we generate the point $H^\text{obs}_r$ using a zero mean Gaussian with standard deviation $\sigma_{21}$ and add the result to $H(z_{21})r_d /(1+z_{21})$. The value of $H^{\text{fid}}(z_i)$ correspond to the fiducial value, which is defined by $q^{\text{fid}}(z)$, $H_0^{\text{fid}}= 73.0 \, \text{km} \, \text{s}^{-1} \text{Mpc}^{-1}$ and $r_d^{\text{fid}}= 103.5 \, h^{-1}\text{Mpc}$.

Applying the same methodology, we generate the first two points of the BAO set $\{D_{V_j}/r_d\}$, where $D_V(z_j)/r_d = D_V^{\text{fid}}(z_j)/r_d^{\text{fid}}+ \delta{}b_j$, the fiducial values are also defined by $q^{\text{fid}}(z)$, $H_0^{\text{fid}}$ and $r_d^{\text{fid}}$, and their respective inverse variances are presented in the first two rows of table \[tab:bao\]. The three correlated BAO points, related to Kazin et al. [@Kazin2014] data, are resampled randomly generating the 3 values $\delta{}b_j$ from a multivariate Gaussian distribution with means 0 and covariance matrix provided in reference [@Kazin2014], see table \[tab:bao\].

The last BAO point ($j = 6$) follows a different probability distribution, which is defined by the likelihood function $L_{Ross}$ provided by Ross et al. [@Ross2014] and centered in the observed value $(D_{V_j}/r_d)^{\text{obs}} = 4.466$. Thus, to produce a mock catalog from a fiducial model we perform the inverse transform sampling. First, we shift this function, $L_{Ross}^s$ such that the new center is $D_V^{\text{fid}}(z_j)/r_d^{\text{fid}}$. Then, we randomly generate a number $u$ from a uniform distribution $(0, 1)$ and, finally, we invert the equation $$u = \int_{\alpha_i}^\alpha {\mathrm{d}}\alpha^\prime L_{Ross}^s (\alpha^\prime)$$ obtaining $\alpha(u)$. The lower bound $\alpha_i$ corresponds to the lower bound described by the $L_{Ross}$ function.

  ------------------------------------- ------ -------- ---------- ----------- ------ -------- ----------
                Reference                 z     $H(z)$   $\sigma$   Reference    z     $H(z)$   $\sigma$
   \[0.5ex\] Riess et al. [@Riess2011]   0.0    $73.8$     2.4                  0.1      69        12
                                         0.18     75        4                   0.17     83        8
                                         0.20     75        5                   0.27     77        14
                                         0.35     83        14                  0.4      95        17
                                         0.59    104        13                  0.48     97        60
                                         0.68     92        8                   0.88     90        40
                                         0.78    105        12                  0.9     117        23
                                         0.88    125        17                  1.3     168        17
                                         1.04    154        20                  1.43    177        18
                                                                                1.53    140        14
        Busca et al. [@Busca2013]         -       -         -                   1.75    202        40
  ------------------------------------- ------ -------- ---------- ----------- ------ -------- ----------

  : $H(z)$ data[]{data-label="tab:Hz"}

The observable of this entry is related to $H(z_{21})r_d/(1+z_{21})$, instead of simply $H(z)$ as the other data points. In this entry the redshift is $z_{21} = 2.3$, the observable value and its standard deviation are $H^{\text{obs}}_r = 224 \, \text{km} \, \text{s}^{-1} \text{Mpc}^{-1}$ and $\sigma_{21} = 8 \, \text{km} \, \text{s}^{-1} \text{Mpc}^{-1}$, respectively.

[67]{}

Collaboration, A. G. Riess, A. V. Filippenko, P. Challis, A. Clocchiattia, A. Diercks, et al., [*Observational evidence from supernovae for an accelerating universe and a cosmological constant*]{}, [*Astron. J.*]{} [**116**]{} (1998) 1009–1038, \[[[astro-ph/9805201]{}](http://xxx.lanl.gov/abs/astro-ph/9805201)\].

Collaboration, S. Perlmutter, G. Aldering, G. Goldhaber, R. A. Knop, P. Nugent, et al., [*Measurements of omega and lambda from 42 high-redshift supernovae*]{}, [*Astrophys. J.*]{} [**517**]{} (1999) 565–586, \[[[ astro-ph/9812133]{}](http://xxx.lanl.gov/abs/astro-ph/9812133)\].

A. [Joyce]{}, B. [Jain]{}, J. [Khoury]{}, and M. [Trodden]{}, [*Beyond the cosmological standard model*]{}, [*Phys. Rep.*]{} [**568**]{} (2015) 1–98, \[[[arXiv:1407.0059]{}](http://xxx.lanl.gov/abs/1407.0059)\].

G. Dvali, G. Gabadadze, and M. Porrati, [*4d gravity on a brane in 5d minkowski space*]{}, [*Phys. Lett. B*]{} [**485**]{} (2000) 208–214, \[[[hep-th/0005016]{}](http://xxx.lanl.gov/abs/hep-th/0005016)\].

T. P. [Sotiriou]{} and V. [Faraoni]{}, [*f(r) theories of gravity*]{}, [*Rev. Mod. Phys.*]{} [**82**]{} (2010) 451–497, \[[[arXiv:0805.1726]{}](http://xxx.lanl.gov/abs/0805.1726)\].

S. [Nojiri]{} and S. D. [Odintsov]{}, [*Unified cosmic history in modified gravity: From f(r) theory to lorentz non-invariant models*]{}, [*Phys. Rep.*]{} [**505**]{} (2011) 59–144, \[[[ arXiv:1011.0544]{}](http://xxx.lanl.gov/abs/1011.0544)\].

M. S. Turner and A. Riess, [*Do sne ia provide direct evidence for past deceleration of the universe?*]{}, [*Astrophys. J.*]{} [**569**]{} (2002) 18, \[[[astro-ph/0106051]{}](http://xxx.lanl.gov/abs/astro-ph/0106051)\].

M. Visser, [*Jerk, snap, and the cosmological equation of state*]{}, [ *Classical Quant. Grav.*]{} [**21**]{} (2004) 2603–2616, \[[[gr-qc/0309109]{}](http://xxx.lanl.gov/abs/gr-qc/0309109)\].

D. Rapetti, S. W. Allen, M. A. Amin, and R. D. Blandford, [*A kinematical approach to dark energy studies*]{}, [*Mon. Not. R. Astron. Soc.*]{} [**375**]{} (2007) 1510–1520, \[[[ astro-ph/0605683]{}](http://xxx.lanl.gov/abs/astro-ph/0605683)\].

Z.-X. [Zhai]{}, M.-J. [Zhang]{}, Z.-S. [Zhang]{}, X.-M. [Liu]{}, and T.-J. [Zhang]{}, [*[Reconstruction and constraining of the jerk parameter from OHD and SNe Ia observations]{}*]{}, [*Phys. Lett. B*]{} [**727**]{} (2013) 8–20, \[[[arXiv:1303.1620]{}](http://xxx.lanl.gov/abs/1303.1620)\].

R. A. [Daly]{} and S. G. [Djorgovski]{}, [*A model-independent determination of the expansion and acceleration rates of the universe as a function of redshift and constraints on dark energy*]{}, [*The Astrophys. J.*]{} [**597**]{} (2003) 9–20, \[[[ astro-ph/0305197]{}](http://xxx.lanl.gov/abs/astro-ph/0305197)\].

S. [Benitez-Herrera]{}, F. [R[ö]{}pke]{}, W. [Hillebrandt]{}, C. [Mignone]{}, M. [Bartelmann]{}, et al., [*Model-independent reconstruction of the expansion history of the universe from type ia supernovae*]{}, [*Mon. Not. R. Astron. Soc.*]{} [**419**]{} (2012) 513–521, \[[[arXiv:1109.0873]{}](http://xxx.lanl.gov/abs/1109.0873)\].

D. Huterer and G. Starkman, [*Parameterization of dark-energy properties: a principal-component approach*]{}, [*Phys. Rev. Lett.*]{} [**90**]{} (2003) 031301, \[[[ astro-ph/0207517]{}](http://xxx.lanl.gov/abs/astro-ph/0207517)\].

R. A. [Daly]{}, S. G. [Djorgovski]{}, K. A. [Freeman]{}, M. P. [Mory]{}, C. P. [O’Dea]{}, et al., [*Improved constraints on the acceleration history of the universe and the properties of the dark energy*]{}, [*Astrophys. J.*]{} [**677**]{} (2008) 1–11, \[[[arXiv:0710.5345]{}](http://xxx.lanl.gov/abs/0710.5345)\].

V. [Sahni]{} and A. [Starobinsky]{}, [*Reconstructing dark energy*]{}, [*Int. J. Mod. Phys. D*]{} [**15**]{} (2006) 2105–2132, \[[[astro-ph/0610026]{}](http://xxx.lanl.gov/abs/astro-ph/0610026)\].

Collaboration, A. G. Riess, L.-G. Strolger, J. Tonry, S. Casertano, H. C. Ferguson, et al., [*Type ia supernova discoveries at z&gt;1 from the hubble space telescope: Evidence for past deceleration and constraints on dark energy evolution*]{}, [*Astrophys. J.*]{} [**607**]{} (2004) 665–687, \[[[astro-ph/0402512]{}](http://xxx.lanl.gov/abs/astro-ph/0402512)\].

C. Shapiro and M. S. Turner, [*What do we really know about cosmic acceleration?*]{}, [*Astrophys. J.*]{} [**649**]{} (2006) 563–569, \[[[astro-ph/0512586]{}](http://xxx.lanl.gov/abs/astro-ph/0512586)\].

A. [Avgoustidis]{}, L. [Verde]{}, and R. [Jimenez]{}, [*Consistency among distance measurements: transparency, bao scale and accelerated expansion*]{}, [*J. Cosmol. Astropart. Phys.*]{} [**6**]{} (2009) 12, \[[[arXiv:0902.2006]{}](http://xxx.lanl.gov/abs/0902.2006)\].

R. [Nair]{}, S. [Jhingan]{}, and D. [Jain]{}, [*Cosmokinetics: a joint analysis of standard candles, rulers and cosmic clocks*]{}, [*J. Cosmol. Astropart. Phys.*]{} [**1**]{} (2012) 18, \[[[ arXiv:1109.4574]{}](http://xxx.lanl.gov/abs/1109.4574)\].

C. [Clarkson]{} and C. [Zunckel]{}, [*Direct reconstruction of dark energy*]{}, [*Phys. Rev. Lett.*]{} [**104**]{} (2010) 211301, \[[[arXiv:1002.5004]{}](http://xxx.lanl.gov/abs/1002.5004)\].

E. E. O. [Ishida]{} and R. S. [de Souza]{}, [*Hubble parameter reconstruction from a principal component analysis: minimizing the bias*]{}, [*Astron. Astrophys.*]{} [**527**]{} (2011) A49, \[[[arXiv:1012.5335]{}](http://xxx.lanl.gov/abs/1012.5335)\].

E. J. [Ruiz]{}, D. L. [Shafer]{}, D. [Huterer]{}, and A. [Conley]{}, [*Principal components of dark energy with supernova legacy survey supernovae: The effects of systematic errors*]{}, [*Phys. Rev. D*]{} [**86**]{} (2012) 103004, \[[[arXiv:1207.4781]{}](http://xxx.lanl.gov/abs/1207.4781)\].

S. [Benitez-Herrera]{}, E. E. O. [Ishida]{}, M. [Maturi]{}, W. [Hillebrandt]{}, M. [Bartelmann]{}, et al., [*[Cosmological parameter estimation from SN Ia data: a model-independent approach]{}*]{}, [*Mon. Not. R. Astron. Soc.*]{} [ **436**]{} (2013) 854–858, \[[[ arXiv:1308.5653]{}](http://xxx.lanl.gov/abs/1308.5653)\].

A. [Shafieloo]{}, U. [Alam]{}, V. [Sahni]{}, and A. A. [Starobinsky]{}, [*Smoothing supernova data to reconstruct the expansion history of the universe and its age*]{}, [*Mon. Not. R. Astron. Soc.*]{} [**366**]{} (2006) 1081–1095, \[[[astro-ph/0505329]{}](http://xxx.lanl.gov/abs/astro-ph/0505329)\].

A. [Shafieloo]{}, [*Crossing statistic: reconstructing the expansion history of the universe*]{}, [*J. Cosmol. Astropart. Phys.*]{} [**8**]{} (2012) 2, \[[[arXiv:1204.1109]{}](http://xxx.lanl.gov/abs/1204.1109)\].

T. [Holsclaw]{}, U. [Alam]{}, B. [Sans[ó]{}]{}, H. [Lee]{}, K. [Heitmann]{}, et al., [ *Nonparametric dark energy reconstruction from supernova data*]{}, [*Phys. Rev. Lett.*]{} [**105**]{} (2010) 241302, \[[[arXiv:1011.3079]{}](http://xxx.lanl.gov/abs/1011.3079)\].

M. Seikel, C. Clarkson, and M. Smith, [*Reconstruction of dark energy and expansion dynamics using gaussian processes*]{}, [*J. Cosmol. Astropart. Phys.*]{} [**6**]{} (2012) 36, \[[[ arXiv:1204.2832]{}](http://xxx.lanl.gov/abs/1204.2832)\].

A. [Montiel]{}, R. [Lazkoz]{}, I. [Sendra]{}, C. [Escamilla-Rivera]{}, and V. [Salzano]{}, [*Nonparametric reconstruction of the cosmic expansion with local regression smoothing and simulation extrapolation*]{}, [*Phys. Rev. D*]{} [**89**]{} (2014) 043007, \[[[ arXiv:1401.4188]{}](http://xxx.lanl.gov/abs/1401.4188)\].

A. Conley, J. Guy, M. Sullivan, N. Regnault, P. Astier, et al., [*Supernova constraints and systematic uncertainties from the first three years of the supernova legacy survey*]{}, [*Astrophys. J. Suppl. Ser.*]{} [**192**]{} (2011) 1, \[[[arXiv:1104.1443]{}](http://xxx.lanl.gov/abs/1104.1443)\].

M. [Betoule]{}, R. [Kessler]{}, J. [Guy]{}, J. [Mosher]{}, D. [Hardin]{}, et al., [ *Improved cosmological constraints from a joint analysis of the sdss-ii and snls supernova samples*]{}, [*Astron. Astrophys.*]{} [**568**]{} (2014) A22, \[[[arXiv:1401.4064]{}](http://xxx.lanl.gov/abs/1401.4064)\].

F. [Beutler]{}, C. [Blake]{}, M. [Colless]{}, D. H. [Jones]{}, L. [Staveley-Smith]{}, et al., [*The 6df galaxy survey: baryon acoustic oscillations and the local hubble constant*]{}, [*Mon. Not. R. Astron. Soc.*]{} [**416**]{} (2011) 3017–3032, \[[[ arXiv:1106.3366]{}](http://xxx.lanl.gov/abs/1106.3366)\].

N. [Padmanabhan]{}, X. [Xu]{}, D. J. [Eisenstein]{}, R. [Scalzo]{}, A. J. [Cuesta]{}, et al., [*A 2 per cent distance to z = 0.35 by reconstructing baryon acoustic oscillations - i. methods and application to the sloan digital sky survey*]{}, [*Mon. Not. R. Astron. Soc.*]{} [**427**]{} (2012) 2132–2145, \[[[arXiv:1202.0090]{}](http://xxx.lanl.gov/abs/1202.0090)\].

E. A. [Kazin]{}, J. [Koda]{}, C. [Blake]{}, N. [Padmanabhan]{}, S. [Brough]{}, et al., [*The wigglez dark energy survey: improved distance measurements to z = 1 with reconstruction of the baryonic acoustic feature*]{}, [*Mon. Not. R. Astron. Soc.*]{} [**441**]{} (2014) 3524–3542, \[[[arXiv:1401.0358]{}](http://xxx.lanl.gov/abs/1401.0358)\].

A. J. [Ross]{}, L. [Samushia]{}, C. [Howlett]{}, W. J. [Percival]{}, A. [Burden]{}, et al., [*The clustering of the sdss dr7 main galaxy sample i: A 4 per cent distance measure at z=0.15*]{}, [[arXiv:1409.3242]{}](http://xxx.lanl.gov/abs/1409.3242).

D. [Stern]{}, R. [Jimenez]{}, L. [Verde]{}, M. [Kamionkowski]{}, and S. A. [Stanford]{}, [*Cosmic chronometers: constraining the equation of state of dark energy. i: H(z) measurements*]{}, [*J. Cosmol. Astropart. Phys.*]{} [**2**]{} (2010) 8, \[[[arXiv:0907.3149]{}](http://xxx.lanl.gov/abs/0907.3149)\].

A. G. [Riess]{}, L. [Macri]{}, S. [Casertano]{}, H. [Lampeitl]{}, H. C. [Ferguson]{}, et al., [*A 3 hubble space telescope and wide field camera 3*]{}, [*Astrophys. J.*]{} [ **730**]{} (2011) 119, \[[[ arXiv:1103.2976]{}](http://xxx.lanl.gov/abs/1103.2976)\].

M. [Moresco]{}, A. [Cimatti]{}, R. [Jimenez]{}, L. [Pozzetti]{}, G. [Zamorani]{}, et al., [*Improved constraints on the expansion rate of the universe up to z \~ 1.1 from the spectroscopic evolution of cosmic chronometers*]{}, [*J. Cosmol. Astropart. Phys.*]{} [**8**]{} (2012) 6, \[[[arXiv:1201.3609]{}](http://xxx.lanl.gov/abs/1201.3609)\].

N. G. Busca, T. Delubac, J. Rich, S. Bailey, A. Font-Ribera, et al., [ *Baryon acoustic oscillations in the ly[$\alpha$]{} forest of boss quasars*]{}, [*Astron. Astrophys.*]{} [**552**]{} (2013) A96, \[[[arXiv:1211.2616]{}](http://xxx.lanl.gov/abs/1211.2616)\].

G. [Hinshaw]{}, D. [Larson]{}, E. [Komatsu]{}, D. N. [Spergel]{}, C. L. [Bennett]{}, et al., [*Nine-year wilkinson microwave anisotropy probe (wmap) observations: Cosmological parameter results*]{}, [*Astrophys. J. Suppl. Ser.*]{} [**208**]{} (2013) 19, \[[[ arXiv:1212.5226]{}](http://xxx.lanl.gov/abs/1212.5226)\].

, P. A. R. [Ade]{}, N. [Aghanim]{}, M. [Arnaud]{}, M. [Ashdown]{}, et al., [*Planck 2015 results. xiii. cosmological parameters*]{}, [[ arXiv:1502.0158]{}](http://xxx.lanl.gov/abs/1502.0158).

E. Komatsu, K. M. Smith, J. Dunkley, C. L. Bennett, B. Gold, et al., [ *Seven-year wilkinson microwave anisotropy probe (wmap) observations: Cosmological interpretation*]{}, [*Astrophys. J. Suppl. Ser.*]{} [**192**]{} (2011) 18, \[[[arXiv:1001.4538]{}](http://xxx.lanl.gov/abs/1001.4538)\].

M. P. Lima, S. Vitenti, and M. J. Rebouças, [*Energy conditions bounds and their confrontation with supernovae data*]{}, [*Phys. Rev. D*]{} [**77**]{} (2008) 083518, \[[[ arXiv:0802.0706]{}](http://xxx.lanl.gov/abs/0802.0706)\].

M. P. Lima, S. D. P. Vitenti, and M. J. Rebouças, [*Energy conditions bounds and supernovae data*]{}, [*Phys. Lett. B*]{} [**668**]{} (2008) 83, \[[[arXiv:0808.2467]{}](http://xxx.lanl.gov/abs/0808.2467)\].

C. [Catto[ë]{}n]{} and M. [Visser]{}, [*Cosmographic hubble fits to the supernova data*]{}, [*Phys. Rev. D*]{} [**78**]{} (2008) 063501, \[[[arXiv:0809.0537]{}](http://xxx.lanl.gov/abs/0809.0537)\].

M. [Penna-Lima]{}, M. [Makler]{}, and C. A. [Wuensche]{}, [*Biases on cosmological parameter estimators from galaxy cluster number counts*]{}, [*J. Cosmol. Astropart. Phys.*]{} [**5**]{} (2014) 39, \[[[arXiv:1312.4430]{}](http://xxx.lanl.gov/abs/1312.4430)\].

R. G. [Crittenden]{}, G.-B. [Zhao]{}, L. [Pogosian]{}, L. [Samushia]{}, and X. [Zhang]{}, [*[Fables of reconstruction: controlling bias in the dark energy equation of state]{}*]{}, [*J. Cosmol. Astropart. Phys.*]{} [**2**]{} (2012) 48, \[[[arXiv:1112.1693]{}](http://xxx.lanl.gov/abs/1112.1693)\].

G.-B. [Zhao]{}, R. G. [Crittenden]{}, L. [Pogosian]{}, and X. [Zhang]{}, [ *[Examining the Evidence for Dynamical Dark Energy]{}*]{}, [*Phys. Rev. Lett.*]{} [**109**]{} (2012) 171301, \[[[ arXiv:1207.3804]{}](http://xxx.lanl.gov/abs/1207.3804)\].

R. A. Daly and S. G. Djorgovski, [*Direct determination of the kinematics of the universe and properties of the dark energy as functions of redshift*]{}, [*Astrophys. J.*]{} [**612**]{} (2004) 652–659, \[[[astro-ph/0403664]{}](http://xxx.lanl.gov/abs/astro-ph/0403664)\].

C. de Boor, [*A Practical Guide to Splines*]{}. Springer, 2001.

C. M. Bishop, [*Neural Networks for Pattern Recognition*]{}. Oxford University Press, 1996.

M. Keijzer and V. Babovic, [*Genetic programming, ensemble methods and the bias/variance tradeoff – introductory investigations*]{}, in [*Genetic Programming*]{} (R. Poli, W. Banzhaf, W. Langdon, J. Miller, P. Nordin, et al., eds.), vol. 1802 of [*Lecture Notes in Computer Science*]{}, pp. 76–90. Springer Berlin Heidelberg, 2000.

L. Wasserman, C. J. Miller, R. C. Nichol, C. Genovese, W. Jang, et al., [ *Non-parametric inference in astrophysics*]{}, [[astro-ph/0112050]{}](http://xxx.lanl.gov/abs/astro-ph/0112050).

S. [Dias Pinto Vitenti]{} and M. [Penna-Lima]{}, [*Numcosmo: Numerical cosmology*]{}, [*Astrophysics Source Code Library*]{} (2014) \[[[ascl:1408.013]{}](http://xxx.lanl.gov/abs/ascl:1408.013)\].

T. M. [Davis]{}, L. [Hui]{}, J. A. [Frieman]{}, T. [Haugb[ø]{}lle]{}, R. [Kessler]{}, et al., [*The effect of peculiar velocities on supernova cosmology*]{}, [ *Astrophys. J.*]{} [**741**]{} (2011) 67, \[[[arXiv:1012.2912]{}](http://xxx.lanl.gov/abs/1012.2912)\].

K. [Thepsuriya]{} and A. [Lewis]{}, [*Accuracy of cosmological parameters using the baryon acoustic scale*]{}, [[ arXiv:1409.5066]{}](http://xxx.lanl.gov/abs/1409.5066).

. [Aubourg]{}, S. [Bailey]{}, J. E. [Bautista]{}, F. [Beutler]{}, V. [Bhardwaj]{}, et al., [*Cosmological implications of baryon acoustic oscillation (bao) measurements*]{}, [[ arXiv:1411.1074]{}](http://xxx.lanl.gov/abs/1411.1074).

Collaboration, D. J. Eisenstein, I. Zehavi, D. W. Hogg, R. Scoccimarro, M. R. Blanton, et al., [*Detection of the baryon acoustic peak in the large-scale correlation function of sdss luminous red galaxies*]{}, [*Astrophys. J.*]{} [**633**]{} (2005) 560–574, \[[[astro-ph/0501171]{}](http://xxx.lanl.gov/abs/astro-ph/0501171)\].

D. J. [Eisenstein]{}, H.-J. [Seo]{}, E. [Sirko]{}, and D. N. [Spergel]{}, [ *Improving cosmological distance measurements by reconstruction of the baryon acoustic peak*]{}, [*Astrophys. J.*]{} [**664**]{} (2007) 675–679, \[[[astro-ph/0604362]{}](http://xxx.lanl.gov/abs/astro-ph/0604362)\].

D. J. Eisenstein and W. Hu, [*Baryonic features in the matter transfer function*]{}, [*Astrophys. J.*]{} [**496**]{} (1998) 605, \[[[astro-ph/9709112]{}](http://xxx.lanl.gov/abs/astro-ph/9709112)\].

A. Lewis, A. Challinor, and A. Lasenby, [*Efficient computation of [CMB]{} anisotropies in closed [FRW]{} models*]{}, [*Astrophys. J.*]{} [**538**]{} (2000) 473–476, \[[[ astro-ph/9911177]{}](http://xxx.lanl.gov/abs/astro-ph/9911177)\].

B. Efron and R. Tibshirani, [*An Introduction to the Bootstrap*]{}. Chapman & Hall/CRC Monographs on Statistics & Applied Probability. Taylor & Francis, 1994.

J. Goodman and J. Weare, [*Ensemble samplers with affine invariance*]{}, [ *CAMCoS*]{} [**5**]{} (2010) 65–80.

D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, [*emcee : The [MCMC]{} hammer*]{}, [*Publ. Astron. Soc. Pac.*]{} [**125**]{} (2013) 306–312.

S. P. Brooks and A. Gelman, [*General methods for monitoring convergence of iterative simulations*]{}, [*J. Comput. Graph. Stat.*]{} [**7**]{} (1998) 434–455.

S. [Capozziello]{}, O. [Farooq]{}, O. [Luongo]{}, and B. [Ratra]{}, [*[Cosmographic bounds on the cosmological deceleration-acceleration transition redshift in f(R) gravity]{}*]{}, [*Phys. Rev. D*]{} [**90**]{} (2014) 044016, \[[[arXiv:1403.1421]{}](http://xxx.lanl.gov/abs/1403.1421)\].

M. V. dos Santos, R. R. R. Reis, and I. Waga, [*Constraining cosmic deceleration-acceleration transition with type ia supernova, bao/cmb and h(z) data*]{}, [[arXiv:1505.0381]{}](http://xxx.lanl.gov/abs/1505.0381).

A. [Heavens]{}, R. [Jimenez]{}, and L. [Verde]{}, [*Standard rulers, candles, and clocks from the low-redshift universe*]{}, [*Phys. Rev. Lett.*]{} [**113**]{} (2014) 241302, \[[[ arXiv:1409.6217]{}](http://xxx.lanl.gov/abs/1409.6217)\].

[^1]: There are also works which explore the properties of the DE equation of state (i.e., assuming the general relativity) using kinematic quantities, e.g., [@Huterer2003; @Daly2008].

[^2]: For a discussion about the issues in deriving a kinematic function from a reconstructed one of another quantity, see [@Sahni2006].

[^3]: The second derivative of the distance, given by a linear spline, is zero everywhere.

[^4]: For a detailed discussion of this problem in the context of the cluster number counts see [@Penna-Lima2014].

[^5]: In practice, once the value of the function at each knot is defined, the coefficients of each cubic polynomial are obtained from the solution of a linear system including not only the knots on its neighborhood but all knots.

[^6]: One can also easily interpret, using the Bayesian point of view, the penalty factor as a prior on the fitted function.

[^7]: The $\sigma_{abs}$ factor guarantees that, even if some $\bar{\hat{q}}_i \simeq 0$ the denominator in the penalty factor does not go to zero.

[^8]: Note that here we are using the MC method to perform integrals in the data $\vec{D}$ given a set of parameters $\vec{\theta}^{\text{fid}}$. Therefore, we are sampling new data from the given likelihood. This procedure is similar but has different applications than the MC used to study the parameter space $\vec{\theta}$ given a data set $\vec{D}$.

[^9]: The likelihoods used in this work were all implemented on top of the data description objects of the Numerical Cosmology library [@DiasPintoVitenti2014]. These objects automatically provide the resample feature for any likelihood using them.

[^10]: There are methods to infer the bias without the knowledge of the true underlying model, e.g., bootstrap [@Efron1994]. In this work we do not explore these approaches, since their application is not straightforward for correlated data.

[^11]: Following the discussion presented in the end of section \[sec:spline\], we also performed the analyses using $n +1 = 6$ and $n+1 = 10$ knots. The results showed a small improvement from 6 to 8 knots, and a negligible variation from 8 to 10 knots.

[^12]: The densities $\Omega_c$ and $\Omega_b$ are the energy densities divided by the critical density today.

[^13]: In the reference, the authors use a different parametrization $(M_B, \Delta_M)$ such that $M_1 = M_B$ and $M_2 = M_B + \Delta_M$.

[^14]: In a more careful analysis, one could calculate the bias for several fiducial models and their respective upper limits. Then, this value could be added to the curve obtained from the real data.

[^15]: This is natural since, in this case, we have more knots and we would need a smaller ${\sigma_{\text{rel}}}$ value to constrain into a linear curve.

[^16]: The algorithm is describe in [@Goodman2010] and was implemented in C in the NumCosmo library [@DiasPintoVitenti2014]. Another unrelated implementation in Python is described in [@Foreman-Mackey2013].

[^17]: The MPSRF is a diagnose tool used to check the convergence of a Markov Chain, which was originally proposed in [@Brooks1998]. A value $<1.2$ indicates the convergence.

[^18]: In this case, we recover $H(z)$ \[see eq. \] considering $H_0 = 70.0 \, \text{km}\,\text{s}^{-1}\text{Mpc}^{-1}$, which is the reference value used in [@Betoule2014].

[^19]: The covariance matrix $\mathsf{C}_{cmp}$ also depends on the intrinsic standard deviation of each SN Ia subsample, i.e., SDSS, SNLS3, low-z and HST (Hubble Space Telescope), $\{\sigma^1_{int}, \sigma^2_{int}, \sigma^3_{int}, \sigma^4_{int}\}$, respectively. In particular, we used original values calibrated in [@Betoule2014].
