Cosmological Parameter Inference with Bayesian Statistics

Abstract

Bayesian statistics and Markov Chain Monte Carlo (MCMC) algorithms have found their place in the field of Cosmology. They have become important mathematical and numerical tools, especially in parameter estimation and model comparison. In this paper, we review some fundamental concepts to understand Bayesian statistics and then introduce MCMC algorithms and samplers that allow us to perform the parameter inference procedure. We also introduce a general description of the standard cosmological model, known as the $$\mathsf{\Lambda }$$CDM model, along with several alternatives, and current datasets coming from astrophysical and cosmological observations. Finally, with the tools acquired, we use an MCMC algorithm implemented in python to test several cosmological models and find out the combination of parameters that best describes the Universe.

1. Introduction

The beginning of the standard cosmology as it is known today emerged after 1920, when the Shapley-Curtis debate was carried out [1]. This debate was held between the astronomers Harlow Shapley and Heber Curtis, resulting in a revolution for astronomy at that time by reaching an important conclusion: “The Universe had a larger scale than the Milky Way”. Several observations at that epoch established that the size and dynamics of the cosmos could be explained by Einstein’s General Theory of Relativity. In its childhood, cosmology was a speculative science based only on a few datasets, and it was characterized by a dispute between two cosmological models: the steady state model and the Big Bang (BB) theory. It was not until 1990, when the amount of data increased enough to discriminate and rule out compelling theories, that the BB model awarded was the most accepted. During the same decade, David Schramm heralded the “Golden Age of Cosmology” at a National Academy of Sciences colloquium [2].

Once the new age of cosmological observations arrived with a large variety of data, it was necessary to confront the cosmological models with such data, usually done through statistics. It is important to stress out that, since we have a unique Universe, we cannot rely on a frequentist interpretation of statistics (we are not able to create multiple Universes and make a frequentist inference of our models). An alternative approach that will help in our task is the Bayesian statistics. In Bayesian statistics, the probability is interpreted as a “degree of belief”, and it may be useful when repetitive processes are complicated to reproduce.

The main aim of this work is to provide an introduction of Bayesian parameter inference and its applications to cosmology. We assume the reader is familiarized with the basic concepts of statistics, but not necessarily with Bayesian statistics. Then, we provide a general introduction to this subject, enough to work out some examples. This review is written in a generic way so that the parameter inference theory may be applicable to any subject, in particular we put into practice the Bayesian concepts into the field of cosmology.

The paper is organized as follows. In Section 2, we point out the main differences between the Bayesian and Frequentist approaches. Then, in Section 3, we introduce the basic mathematical concepts in Bayesian statistics to perform the parameter estimation procedure for a given model. Once we have the mathematical background, we continue, in Section 4, with some numerical resources that are able to simplify our task, especially for models with several parameters that need to be tested with many datasets. With these methods and tools in place, we provide the example of fitting a straight line in Section 5. In Section 6, we present an introduction to cosmology, and then, in Section 7, focus on some codes to compute the cosmological observables. In Section 8, we constrain the parameter space that describes the standard cosmological model, namely the $\mathsf{\Lambda }$CDM model, along with several extensions. Finally, in Section 9, we present our conclusions.

2. Bayesian vs. Frequentist Statistics

Fundamentally, the main difference between Bayesian and Frequentist statistics is on the definition of probability. From a Frequentist point of view, probability has meaning in limiting cases of repeated measurements

(1)$P=\frac{n}{N},$

where n denotes the number of successes, and N the total number of trials. Frequentist statistics defines the probability P as the limit of the number of independent trials going to infinity. Then, for Frequentist statistics, probabilities are fundamentally related to frequencies of events. On the other hand, in Bayesian statistics, the concept of probability is extended to cover degrees of certainty about a statement. For Bayesian statistics, probabilities are fundamentally related to our knowledge concerning an event.

Here, we introduce some key concepts to understand the consequences this difference entails; for an extended review see Reference [3,4,5,6,7,8,9,10,11,12,13] and references therein. Let x be a random variable related to a particular event and $P\left(x\right)$ its corresponding probability distribution, for both cases, the same rules of probabilities apply ¹:

(2)$P\left(x\right)\ge 0,$

(3)${\int }_{-\infty }^{\infty }dxP\left(x\right)=1.$

For mutually exclusive events, we have

(4)$P(x_1\cup x_2)=P\left(x_1\right)+P\left(x_2\right),$

but, in general

$P(x_1\cup x_2)=P\left(x_1\right)+P\left(x_2\right)-P(x_1\cap x_2).$

These rules are summed up as follows: the first condition (2) is necessary due to the probability of having an event is always positive; the second rule (3) is a normalized relation, which tells us that we are certain to obtain one of the possible outcomes; now, in the third point (4) we have that the probability of obtaining an observation, from a set of mutually exclusive events, is given by the individual probabilities of each event; finally, and in general, if one event occurs given the occurrence of another then the probability that both $x_1$ and $x_2$ happen is equal to the probability of $x_1$ times the probability of $x_2$ given that $x_1$ has already happened

(5)$P(x_1\cap x_2)=P\left(x_1\right)P\left(x_2\right|x_1).$

If two events $x_1$ and $x_2$ are mutually exclusive, then

(6)$P(x_1\cap x_2)=0=P(x_2\cap x_1).$

The rules of probability distributions must be fulfilled by both Frequentist and Bayesian statistics. However, there are some consequences derived by the fact that these two scenarios have a different definition of probability, as we shall see below.

2.1. Frequentist Statistics

Any frequentist inferential procedure relies on three basic ingredients: the data, the model and an estimation procedure. The main assumption in Frequentist statistics is that the data has a definite, albeit unknown, underlying distribution to which all inference pertains.

The data is a measurement or observation, denoted by X, that can take any value from a corresponding sample space. A sample space of an observation X can be defined as a measurable space $(x,\widehat{B})$ that contains all values that X can take upon measurement. In Frequentist statistics, it is considered that there is a probability function $P_0:\widehat{B}\to [0,1]$ in the sample space $(x,\widehat{B})$ representing the “true distribution of the data”

$X\sim P_0.$

For Frequentist statistics, the model Q is a collection of probability measurements $P_{\theta }:\widehat{B}\to [0,1]$ in the sample space $(x,\widehat{B})$. The distributions $P_{\theta }$ are called model distributions, with $\theta$ being the model parameters; in this approach $\theta$ is unchanged. A model Q is said to be well-specified if it contains the true distribution of the data $P_0$, i.e.,

$P_0\in Q.$

Finally, we need a point-estimator (or estimator) for $P_0$. An estimator for $P_0$ is a map $\widehat{P}:x\to Q$, representing our best guess $\widehat{P}\in Q$ for $P_0$ based on the data X. Therefore, the Frequentist statistics is based on trying to answer the following questions: “what the data is trying to tell about $P_0$?” or “considering the data, what can we say about the mean value of $P_0$?”

2.2. Bayesian Statistics

In Bayesian statistics, data and model are two elements of the same space [3], i.e., no formal distinction is made between measured quantities X and parameters $\theta$. One may envisage the process of generating a measurement’s outcome $Y=y$ as two draws, one draw for $\mathsf{\Theta }$ (where $\mathsf{\Theta }$ is a model with associated probabilities to the parameter $\theta$) to select a value of $\theta$ and a subsequent draw for $P_{\theta }$ to arrive at $X=x$. This perspective may seem rather absurd when thinking in a Frequentist way, but, in Bayesian statistics, where probabilities are related to our own knowledge, it results natural to associate probability distributions to our parameters. In this way, an element $P_{\theta }$ of the model is interpreted simply as the distribution of X given the parameter value $\theta$, i.e., as the conditional distribution $X|\theta$.

2.3. Comparing Both Descriptions

Table 1 provides a short summary of the most important differences between the two statistics. To understand these differences, let us review a basic example. Here, we present an experiment and, since we are interested in comparing both descriptions, we show only the basic results from both points of view: Frequentist and Bayesian.

Example 1.

Let us assume we have a coin that has a probability p to land as heads and a probability $1-p$ to land as tails. Our goal is to know whether this coin is fair ($p=0.5$) or not. In the process of trying to estimate p, we flip the coin 14 times, obtaining heads in 10 of the trials. Now, we are interested in the next two possible events. To be precise: “What is the probability that in the next two tosses we will get two heads in a row?”

• Frequentist approach. As mentioned previously, in Frequentist statistics probability is related to the frequency of events, then our best estimate for p is $P\left(head\right)=p=\frac{\#ofheads}{\#ofevents}=10/14$. So, the probability of having 2 heads in a row is $P\left(2heads\right)=P\left(head\right)P\left(head\right)\simeq 0.51$.

• Bayesian approach. In Bayesian statistics, p is not a value; it is a random variable with its own distribution, and it must be defined by the existing evidence (the 14 trials and 10 successes). Then, by considering that we do not know anything about p a priori and averaging over all possible values of p, we have that the probability of having two heads is

  (7) $P\left(2heads\right|D)=0.485.$

  This Bayesian example will be expanded in detail during the following section, but, for now, we just want to stress out that both approximations arrive at different results.

In the Frequentist approach, since we adopt the probability as a frequency of events (the probability of having a head was fixed by $p=10/14$); hence, the final result was obtained by only multiplying each of these probabilities (since we assume the events are independent of each other). On the other hand, in the Bayesian framework, it was necessary to average over all possible values of p in order to obtain a numerical value. However, in both cases, the probability differs from the real one ($P\left(2heads\right)=0.25$) because we do not have enough data for our estimations.

In this example, we have seen the forward application of statistics: using a mathematical model to relate measured quantities to an unknown quantity of interest, in this case, the probability of getting two heads in a row. This can also be applicable to the inverse problem. That is, by the having the data, we would like to obtain information about the parameters of a given model [4,14]. In Section 8, we will illustrate this point by using different cosmological observations and finding out the best-fit values of the parameters that describe a given model.

3. A First Look at Bayesian Statistics

Before we start with the applications of Bayesian statistics in cosmology, it is advisable to understand the important mathematical tools within the Bayesian procedure. In this section, we present a basic revision but encourage the reader to look for the formal treatment in the literature, cited in each section.

3.1. Bayes’ Theorem, Priors, Posteriors, and All That Stuff

When anyone is interested in the Bayesian framework, there are several concepts to understand before presenting the results. In this section, we review these concepts, and then we get back to the example of the coin toss given in the last section.

The Bayes’ theorem. The Bayes’ theorem is a direct consequence of the axioms of probability shown in Equations (2)–(5). From Equation (5), without loss of generality, it must be fulfilled that $P(x_1\cap x_2)=P(x_2\cap x_1)$. In such a case, the following relation applies:

(8)$P\left(x_2\right|x_1)=\frac{P\left(x_1\right|x_2)P\left(x_2\right)}{P\left(x_1\right)}.$

As already mentioned, in the Bayesian framework, data and model are part of the same space. Given a model (or hypothesis) H, considering $x_1\to D$ as a set of data, and $x_2\to \theta$ as the parameter vector of said hypothesis, we can rewrite the above equation as

(9)$P\left(\theta \right|D,H)=\frac{P\left(D\right|\theta ,H\left)P\right(\theta \left|H\right)}{P\left(D\right|H)}.$

This last relation is the so-called Bayes’ theorem and the most important tool in a Bayesian inference procedure. In this result, $P\left(\theta \right|D,H)$ is called the posterior probability of the model. $L\left(D\right|\theta ,H)\equiv P(D|\theta ,H)$ is called the likelihood, and it will be our main focus in future sections, $\pi \left(\theta \right)\equiv P\left(\theta \right|H)$ is called the prior and expresses the knowledge about the model before acquiring the data (this prior can be fixed depending on either previous experiment results or the theory behind), $\mathcal{Z}\equiv P\left(D\right|H)$ is the evidence of the model, usually referred to as the Bayesian Evidence.

The prior refers to the information one has a priori of the model. It can be defined in various ways; however, a common one is the uniform prior (also referred to as a flat prior):

(10)$\pi \left(\theta \right)\propto c,$

with c being a constant. This type of prior is telling us that every parameter value is equally probable a priori, as seen in Figure 1. Using this prior also means that the posterior probability will be proportional to the likelihood (since the Bayesian Evidence is a constant). Another convenient prior distribution is the beta distribution $B(\theta ;a,b)$ since it contains several statistical distributions by varying its parameters a and b (in particular the flat prior is obtained when $a=b=1$). It is defined as

(11)$B(\theta ;a,b)=\frac{\mathsf{\Gamma }(a+b)}{\mathsf{\Gamma }\left(a\right)\mathsf{\Gamma }\left(b\right)}{\theta }^{a-1}{(1-\theta )}^{b-1},$

where $\mathsf{\Gamma }$ is the gamma function. These are just two examples of useful priors, and it is evident that the choice of a prior will influence the posterior distribution, although its effect is reduced as more data are collected, as we shall see later in this section.

Now, regarding the Bayesian Evidence, we notice that it acts as a normalizing factor, and is nothing more than the average of the likelihood over the

(12)$P\left(D\right|H)=\int d^N\theta P\left(D\right|\theta ,H)P\left(\theta \right|H),$

where N is the dimensionality of the parameter space. This quantity is usually ignored, for practical reasons, i.e., when testing the parameter space of a unique model. Nevertheless, the Bayesian evidence plays an important role for selecting the model that best describes the data, this process being known as model selection. For convenience, the ratio of two evidences

(13)$K\equiv \frac{P\left(D\right|H_0)}{P\left(D\right|H_1)}=\frac{\int d^{N_0}{\theta }_{0}P\left(D\right|{\theta }_0,H_0)P\left({\theta }_0\right|H_0)}{\int d^{N_1}{\theta }_{1}P\left(D\right|{\theta }_1,H_1)P\left({\theta }_1\right|H_1)}=\frac{{\mathcal{Z}}_0}{{\mathcal{Z}}_1},$

or equivalently the difference in log evidence $ln{\mathcal{Z}}_0-ln{\mathcal{Z}}_1$ is often termed as the Bayes factor ${\mathcal{B}}_{0,1}$:

(14)${\mathcal{B}}_{0,1}=ln\frac{{\mathcal{Z}}_0}{{\mathcal{Z}}_1},$

where ${\theta }_i$ is a parameter vector (with dimensionality $N_i$) for the hypothesis $H_i$ and $i=0,1$. In Equation (14), the quantity ${\mathcal{B}}_{0,1}=lnK$ provides an idea on how well model 0 may fit the data when compared to model 1. Jeffreys provided a suitable guideline scale on which we are able to make qualitative conclusions (see Table 2 [15]).

We can see that Bayes’ theorem has an enormous implication with respect to a statistical inferential point of view. In a typical scenario, we collect some data and then interpret it with a given model; however, we usually do the opposite. That is, first we have a set of data, and then we confront a model considering the probability that our model fits the data. Bayes’ theorem provides a tool to relate both scenarios. Then, based on the Bayes’ theorem, we are able to select the model that best fits the data.

Example 2.

We go back to the example shown in the last section: the coin toss. We are interested in the probability of obtaining two heads in a row given the data $P\left(2heads\right|D)$ (D being the previous 14 coin tosses acting as data). First, let us assume that we have a model with a parameter p that defines the probability of obtaining the two heads, that is $P\left(2heads\right|p)$. This parameter p will have a probability distribution $P\left(p\right|D)$ depending on the data in place. Therefore, the probability can be obtained by averaging over all the possible parameters with its corresponding density distribution

(15) $P\left(2heads\right|D)={\int }_0^{1}P\left(2heads\right|p)P\left(p\right|D)dp.$

For simplicity, we do not update p between the two tosses, but we assume that both of them are independent of each other. With this last assumption, we have

(16) $P\left(2heads\right|p)={\left[P\left(head\right|p)\right]}^2,$

where $P\left(head\right|p)$ is the probability of obtaining a head given our model. We assume a simple description of $P\left(head\right|p)$ as

(17) $P\left(head\right|p)=p\Rightarrow P\left(2heads\right|p)=p^2.$

On the other hand, notice that we do not know a priori the quantity $P\left(p\right|D)$ but $P\left(D\right|p)$ (i.e., we know the probability of obtaining a dataset by considering a model as correct). A good choice for experiments that have two possible outcomes is the binomial distribution

(18) $P\left(x\right|p,n)=\left(\genfrac{}{}{0pt}{}{n}{x}\right)p^x{(1-p)}^{n-x},$

with n the number of trials (in this case, = 14) and x the number of successes (here =10). Hence, we have an expression for $P\left(D\right|p)$. Using the Bayes’ formula, we have

(19) $P\left(p\right|D)=\frac{P\left(D\right|p\left)P\right(p)}{P\left(D\right)}.$

For the prior, we will use the beta distribution (11), so

(20) $P\left(p\right)=B(p;a,b).$

In order to get the explicit form of $P\left(p\right|D)$, we still need to compute $P\left(D\right)$. That is, by plugging Equations (18) and (20) into the integral of Equation (12) yields to

(21) $P\left(D\right)=B(10+a,4+b)\equiv \frac{\mathsf{\Gamma }(10+a)\mathsf{\Gamma }(4+b)}{\mathsf{\Gamma }\left(\right(10+a)+(4+b\left)\right)};$

therefore,

(22) $P\left(p\right|D)=\frac{p^{10+a-1}{(1-p)}^{4+b-1}}{B(10+a,4+b)}.$

If we know nothing about p, then we can assume the prior is a uniform distribution; this means $a=b=1$. Notice from Figure 1 that our posterior result (Red figure) described by Equation (22) does not exactly agree with the real value of p (black dashed vertical line). We would expect the posterior distribution to be centered at $p=0.5$ with a very narrow distribution. Nevertheless, this value is recovered by increasing the experimental data.

Finally, solving the integral in Equation (15) using (17) and (22), we arrive at the result obtained in the previous section

(23) $P\left(2heads\right|D)=\frac{B(13,5)}{B(11,5)}=0.485.$

It is important to clarify that the inferred value for the parameter p is merely the probability of said value given our data D. This value can change and will generally be a better estimation if more data are collected (as will be seen in detail in the next section).

3.2. Updating the Probability Distribution

As seen in the coin example, we were not able to get the real value of p because of the lack of enough data. If we want to get closer, we would have to keep flipping the coin until the amount of data becomes sufficient. Let us continue with the example: suppose that after throwing the coin 100 times we obtain, let us say 56 heads, while after throwing it 500 times, we obtain 246 heads. Then, we expect to obtain a thinner distribution with its center close to $p=0.5$ (see Figure 2). Given this, it is clear that in order to confront a parameter and be more accurate about the most probable (or “real”) value, it is necessary to increase the amount of data (and the precision) in any experiment. That is, if we take into account the 500 tosses—with 246 heads—the previous result is updated to $P\left(2heads\right|D)=0.249$, much closer to the real value.

Then, we have some model parameters that have to be confronted with different sets of data. This can be done in two alternative ways: (a) by considering the set of all datasets; or (b) by taking each dataset as the new data, but our prior information updated by the previous information. The important point in Bayesian statistics is that it is indeed equivalent to choose any of these two possibilities. In the coin toss example, it means that it is identical to start with the prior given in Figure 2a, and then, by considering the 500 data points, we can arrive at the posterior in Figure 2d, or similarly start with the posterior shown in Figure 2c as our prior and consider only the last 400 data points to obtain the same posterior, displayed in Figure 2d.

In fact, if we rewrite Bayes’ theorem so that all probabilities are explicitly dependent on some prior information I [4]

(24)$P\left(\theta \right|DI,H)=\frac{P\left(\theta \right|I,H\left)P\right(DI|\theta ,H)}{P\left(D\right|I,H)},$

and then consider a new set of data $D^{\prime }$, letting the old data become part of the prior information $I^{\prime }=DI$, we arrive at

(25)$\begin{array}{rrr}P\left(\theta \right|D^{\prime }{I}^{\prime },H)& =& \frac{P\left(\theta \right|I,H)P\left(D{D}^{\prime }I\right|\theta ,H)}{P\left(D{D}^{\prime }\right|I,H)}\hfill \\ & =& P\left(\theta \right|\left[D{D}^{\prime }\right]I,H),\hfill \end{array}$

where we can explicitly see the equivalence of the two different options.

3.3. About the Likelihood

We mentioned that the Bayesian evidence is usually set apart when doing any inference procedure in the parameter space of a single model. Then, without loss of generality, we can fix it to $P\left(D\right|H)=1$. If we ignore the prior (a motivation for doing this will be expanded upon in the next section), we can identify the posterior with the likelihood $P\left(\theta \right|D,H)\propto L(D|\theta ,H)$; thus, by maximizing it, we can find the most probable set of parameters for a model given the data. However, having ignored $P\left(D\right|H)$ and the prior, we are not able to provide an absolute probability for a given model, but only relative probabilities. On the other hand, it is possible to report results independently of the prior by using the Likelihood ratio. The likelihood at a particular point in the parameter space can be compared with the best-fit value, or the maximum likelihood $L_{\max }$. Then, we can say that some parameters are acceptable if the likelihood ratio

(26)$\mathsf{\Lambda }=-2ln\left[\frac{L\left(D\right|\theta ,H)}{L_{\max }}\right],$

is bigger than a given value.

Let us assume we have a single-peaked distribution. We consider that $\widehat{\theta }$ is the mean of the distribution

(27)$\widehat{\theta }=\int d\theta \theta P\left(\theta \right|D,H).$

If our model is well-specified and the expectation value of $\widehat{\theta }$ corresponds to the real or most probable value ${\theta }_0$, we have

(28)$⟨\widehat{\theta }⟩={\theta }_0,$

then we say that $\widehat{\theta }$ is unbiased. Considering a Taylor expansion of the log likelihood around its maximum, we have

(29)$\begin{array}{rrr}lnL\left(D\right|\theta )& =& lnL\left(D\right|{\theta }_0)\hfill \\ & +& \frac{1}{2}({\theta }_i-{\theta }_{0i})\frac{{\partial }^{2}lnL}{\partial {\theta }_i\partial {\theta }_j}({\theta }_j-{\theta }_{0j})\hfill \\ & +& ...,\hfill \end{array}$

where ${\theta }_0$ corresponds to the parameter vector of the real model. In this manner, we have that the likelihood can be expressed as a multi-variable likelihood given by

(30)$L\left(D\right|\theta )=L\left(D\right|{\theta }_0)exp\left[-\frac{1}{2}({\theta }_i-{\theta }_{0i})H_{ij}({\theta }_j-{\theta }_{0j})\right],$

where

(31)$H_{ij}=-\frac{{\partial }^{2}lnL}{\partial {\theta }_i\partial {\theta }_j},$

is called the Hessian matrix. It controls whether the estimates of ${\theta }_i$ and ${\theta }_j$ are correlated, and if it is diagonal, these estimates are uncorrelated.

The above expression for the likelihood is a good approximation as long as our posterior distribution possesses a single-peak. It is worth mentioning that, if the data errors are normally distributed, then the likelihood for the data will be a Gaussian function as well. In fact, this is always true if the model is linearly dependent on the parameters. On the other hand, if the data is not normally distributed, we can resort to the central limit theorem. In this way, the central limit theorem tells us that the resulting distribution will be best approximated by a multi-variate Gaussian distribution [6].

3.4. Letting Aside the Priors

In this section, we present an argument for letting aside the prior in the parameter estimation. For this, we follow the example given in Reference [5]. In this example, there are two people, A and B, that are interested in the measurement of a given physical quantity $\theta$. A and B have different prior beliefs regarding the possible value of $\theta$. This discrepancy could be given by the experience, such as the possibility that A and B have made a similar measurement at different times. Let us denote their priors by $P\left(\theta \right|I_i)$, $(i=A,B)$, and assume they are described by two Gaussian distributions with mean ${\mu }_i$ and variance ${\Sigma }_i^2$. Now, A and B measure $\theta$ together using an apparatus subject to a Gaussian noise with known variance $\sigma$. They obtain the value ${\theta }_0=m_1$. Therefore, they can write their likelihoods for $\theta$ as

(32)$L\left(D\right|\theta ,HI)=L_{0}exp\left[-\frac{1}{2}\frac{{(\theta -m_1)}^2}{{\sigma }^2}\right].$

By using the Bayes formula, the posterior of the model A (and B) becomes

(33)$P\left(\theta \right|m_1)=\frac{L\left(m_1\right|\theta I_i)P\left(\theta \right|I_i)}{P\left(m_1\right|I_i)},$

where we have skipped writing explicitly the hypothesis H and used the notation given in Equation (24). Then, the posterior of A and B are (again) Gaussian with mean

(34)${\widehat{\mu }}_i=\frac{m_1+{(\sigma /{\Sigma }_i)}^2{\mu }_i}{1+{(\sigma /{\Sigma }_i)}^2},$

and variance

(35)${\tau }_i^2=\frac{{\sigma }^2}{1+{(\sigma /{\Sigma }_i)}^2},(i=A,B).$

Thus, if the likelihood is more informative than the prior, i.e., $(\sigma /{\Sigma }_i)\ll 1$, the posterior mean of A (and B) will converge towards the measured value $m_1$. As more data are obtained, one can simply replace the value of $m_1$ in the above equation by the mean $⟨m⟩$ and ${\sigma }^2$ by ${\sigma }^2/N$. Then, we can see that the initial prior ${\mu }_i$ of A and B will progressively be overridden by the data. This process is illustrated in Figure 3 where the green (red) curve corresponds to the probability distribution of $\theta$ for person A (B) and the blue curve corresponds to their likelihood.

3.5. Chi-Square and Goodness of Fit

We mentioned that the main aim of parameter estimation is to maximize the likelihood in order to obtain the most probable set of model parameters given the data. If we consider the Gaussian approximation given in Equation (30), we can see the likelihood will be maximum if the quantity

(36)${\chi }^2\equiv ({\theta }_i-{\theta }_{0i})H_{ij}({\theta }_j-{\theta }_{0j}),$

is minimum. The quantity ${\chi }^2$ is usually called chi-square and is related to the Gaussian likelihood via $L=L_0{e}^{-{\chi }^2/2}$. Then, we can say that maximizing the Gaussian likelihood is equivalent to minimizing the chi-square. However, as we mentioned before, there are some circumstances where the likelihood cannot be described by a Gaussian distribution, in these cases the chi-square and the likelihood are no longer equivalent.

The probability distribution for different values of ${\chi }^2$ around its minimum is given by the ${\chi }^2$ distribution for $v=n-M$ degrees of freedom, where n is the number of independent data points and M the number of parameters. Hence, we can calculate the probability that an observed ${\chi }^2$ exceeds, by chance, a value $\widehat{\chi }$ for the correct model. This probability is given by $Q(v,\widehat{\chi })=1-\mathsf{\Gamma }(v/2,\widehat{\chi }/2)$ [16], where $\mathsf{\Gamma }$ is the incomplete Gamma function. Then, the probability that the observed ${\chi }^2$ (even the correct model) is less than a given value ${\widehat{\chi }}^2$ is $1-Q$. This statement is strictly true if the errors are Gaussian and the model is a linear function of the likelihood, i.e., for Gaussian likelihoods.

If we evaluate the quantity Q for the best-fit values (minimum chi-square), we can have a measure of the goodness of fit. If Q is small (small probability), we can interpret it as:

• The model is wrong and can be rejected.

• The errors are underestimated.

• The error measurements are not normally distributed.

On the other hand, if Q is too large there are some reasons to believe that:

• Errors have been overestimated.

• Data are correlated or non-independent.

• The distribution is non-Gaussian.

3.6. Contour Plots and Confidence Regions

Once the best fit parameters are obtained, we would like to know the confidence regions where values could be considered good candidates for our model. The most logical election is to take values inside a compact region around the best fit value. Then, a natural choice are regions with constant ${\chi }^2$ boundaries. When the ${\chi }^2$ possesses more than one minimum, it is said that we have non-connected confidence regions, and for multi-variate Gaussian distributions (as the likelihood approximation in Equation (30)) these are ellipsoidal regions. In this section, we exemplify how to calculate the confidence regions, following Reference [6].

We consider a small perturbation from the best fit of chi-square $\mathsf{\Delta }{\chi }^2={\chi }^2-{\chi }_{\mathrm{best}}^2$. Then, we use the properties of ${\chi }^2$ distribution to define confidence regions for variations on ${\chi }^2$ to its minimum. In Table 3, we see the typical $68.3\%$, $95.4\%$, and $99.73\%$ confidence levels as a function of number of parameters M for the joint confidence level. For Gaussian distributions, these correspond to the conventional 1, 2, and 3 $\sigma$ confidence levels. As an example, we plot in Figure 4 the corresponding confidence regions associated with the coin example.

The general recipe to compute constant ${\chi }^2$ confidence regions is as follows: after finding the best fit by minimizing ${\chi }^2$ (or maximizing the likelihood) and checking that Q is acceptable for the best parameters, then:

• 1.. Let M be the number of parameters, n the number of data and p the confidence limit desired.

• 2.. Solve the equation:

  (37)$Q(n-M,min\left({\chi }^2\right)+\mathsf{\Delta }{\chi }^2)=p.$

• 3.. Find the parameter region where ${\chi }^2\le min\left({\chi }^2\right)+\mathsf{\Delta }{\chi }^2$. This defines the confidence region.

3.7. Marginalization

It is clear that a model may (in general) depend on more than one parameter. However, some of these parameters ${\theta }_i$ may be of less interest. For example, they may correspond to nuisance parameters, like calibration factors, or it may be the case that we are interested in only one of the parameter constraints rather than the joint of two or more of them simultaneously. Then, we marginalize over the uninteresting parameters by

(38)$P({\theta }_1,...,{\theta }_j,H|D)=\int d{\theta }_{j+1}...d{\theta }_{m}P(\theta ,H|D),$

where m is the total number of parameters in our model, and ${\theta }_1$,...,${\theta }_j$ denote the parameters we are interested in.

3.8. Fisher Matrix

Once we have a dataset, it is important to know the accuracy for which we can estimate parameters. Fisher suggested a way 70 years ago [17]. Let us start by considering again a Gaussian likelihood. As we notice, the Hessian matrix $H_{ij}$ has information on the parameter errors and their covariance. More specifically, when all parameters are fixed except one (e.g., the i-th parameter), its error is $1/\sqrt{H_{ii}}$. These errors are called conditional errors, although they are rarely used.

A quantity to forecast the precision of a model, that arises naturally with Gaussian likelihoods, is the so-called Fisher information matrix

(39)$F_{ij}=-⟨\frac{{\partial }^2\mathcal{L}}{\partial {\theta }_i\partial {\theta }_j}⟩,$

where

(40)$\mathcal{L}=lnL.$

It is clear that $F=⟨H⟩$, where the average is made with observational data.

As we can see from Equation (4), for independent datasets, the complete likelihood is the product of the likelihoods, and the Fisher matrix is the sum of individual Fisher matrices. A pedagogical and easy case is having one-parameter ${\theta }_i$ with a Gaussian likelihood. In this scenario,

(41)$\mathsf{\Delta }\mathcal{L}=\frac{1}{2}{F}_{ii}{({\theta }_i-{\theta }_{0i})}^2,$

when $2\mathsf{\Delta }\mathcal{L}=1$, and identifying the $\mathsf{\Delta }{\chi }^2$ corresponding to $68\%$ confidence level, we notice that $1/\sqrt{F_{ii}}$ yields the $1-\sigma$ displacement for ${\theta }_i$. In the general case,

(42)${\sigma }_{ij}^2\ge {\left(F^{-1}\right)}_{ij}.$

Thus, when all parameters are estimated simultaneously from the data, the marginalized error is

(43)${\sigma }_{{\theta }_i}\ge {\left(F^{-1}\right)}_{ii}^{1/2}.$

The beauty of the Fisher matrix approach is that there is a simple prescription for setting it up by only knowing the model and measurement uncertainties, and under the assumption of a Gaussian likelihood the Fisher matrix is the inverse of the covariance matrix. So, all we have to do is set up the Fisher matrix and then invert it to obtain the covariance matrix (that is, the uncertainties on the model parameters). In addition, its fast calculation also enables one to explore different experimental setups and optimize the experiment.

The main point of the Fisher matrix formalism is to predict how well the experiment will be able to constrain the parameters, of a given model, before doing the experiment and perhaps even without simulating it in any detail. We can then forecast the results of different experiments and look at trade-offs, such as precision versus cost. In other words, we can engage in experimental design. The inequality in Equation (42) is called the Kramer-Rao inequality. One can see that the Fisher information matrix represents a lower bound of the errors. Only when the likelihood is normally distributed the inequality is transformed into an equality. However, as we saw in Section 3.3, a Gaussian likelihood is only applicable to some circumstances, being generally impossible to be applied, so the key is to have a good understanding of our theoretical model in such a way that we can construct a Gaussian likelihood.

Constructing Fisher Matrices: A Simple Description

Let us construct Fisher matrices in a simple way. Suppose we have a model that depends on N parameters ${\theta }_1,{\theta }_2,...,{\theta }_N$. We consider M observables $f_1,f_2,...,f_M$ each one related to the model parameters by some equation $f_i=f_i({\theta }_1,{\theta }_2,...,{\theta }_N)$. Then, the elements of the Fisher matrix can be computed as

(44)$F_{ij}=\sum _k\frac{1}{{\sigma }_k^2}\frac{\partial f_k}{\partial {\theta }_i}\frac{\partial f_k}{\partial {\theta }_j},$

where ${\sigma }_k$ are the errors associated with each observable, and we have considered they are Gaussianly distributed. Here, instead of taking the real data values (which could be unknown), it is possible to recreate the data with a fiducial model. The errors associated with the mock data can be taken as the expected experimental errors, and then it is possible to calculate the above expression.

To complement the subject, there is also the Figure of Merit used by the Dark Energy Task Force (DETF) [18], which is defined as the reciprocal of the area in the plane enclosing the $95\%$ confidence limit of two parameters. The larger the figure of merit the greater accuracy one has measuring said parameters. The figure of merit can also be used to see how different experiments break degeneracies, and to predict accuracy in future experiments (experimental design).

3.9. Importance Sampling

We call Importance Sampling (IS) to different techniques of determining properties of a distribution by drawing samples from another one. The main idea is that the distribution one samples from should be representative of the distribution of interest (for a larger number of samples). In such case, we should infer different quantities out of it. In this section, we review the basic concepts necessary to understand the IS, following the Reference [19].

Suppose we are interested in computing the expectation value ${\mu }_f=E_p\left[f\left(X\right)\right]$, where $f\left(X\right)$ is a probability density of a random variable X and the sub-index p means average over the distribution p. Then, if we consider a new probability density $q\left(x\right)$ that satisfies $q\left(x\right)>0$ whenever $f\left(x\right)p\left(x\right)\ne 0$, we can rewrite the mean value ${\mu }_f$ as

(45)$\begin{array}{rrr}{\mu }_f=\int f\left(x\right)p\left(x\right)dx& =& \int f\left(x\right)\frac{p\left(x\right)}{q\left(x\right)}q\left(x\right)dx\hfill \\ & =& E_q\left[f\left(X\right)w\left(x\right)\right],\hfill \end{array}$

where $w\left(x\right)=p\left(x\right)/q\left(x\right)$, and now we have an average over q. So, if we have a collection of different draws $x^{\left(1\right)},...,x^{\left(m\right)}$ from $q\left(x\right)$, we can estimate ${\mu }_f$ using these draws as

(46)${\widehat{\mu }}_f=\frac{1}{m}\sum _{j=1}^{m}w\left(x^{\left(j\right)}\right)f\left(x^{\left(j\right)}\right).$

If $p\left(x\right)$ is known only up to a normalizing constant, the above expression can be calculated as a ratio estimate:

(47)${\widehat{\mu }}_f=\frac{{\sum }_{j=1}^{m}w\left(x^{\left(j\right)}\right)f\left(x^{\left(j\right)}\right)}{{\sum }_{j=1}^{m}w\left(x^{\left(j\right)}\right)}.$

For the strong law of large numbers, in the limit when $m\to \infty$, we will have that ${\widehat{\mu }}_f\to {\mu }_f$.

Another useful quantity to compute in Bayesian analysis is the ratio between evidences for two different models

(48)$\frac{P^{\prime }\left(D\right)}{P\left(D\right)}=E{\left[\frac{P^{\prime }(\theta ,D)}{P(\theta ,D)}\right]}_{P\left(\theta \right|D)}\simeq \frac{1}{N}\sum _{n=1}^N\frac{P^{\prime }\left(D\right|{\theta }_n)P^{\prime }\left({\theta }_n\right)}{P\left(D\right|{\theta }_n)P\left({\theta }_n\right)},$

where the samples $\left\{{\theta }_n\right\}$ are drawn from $P\left(\theta \right|D)$.

An important result for importance sampling is that, if we have a new set of data which is broadly consistent with the current data (in the sense that the posterior only shrinks), we can make use of importance sampling in order to quickly calculate a new posterior including the new data.

3.10. Combining Datasets: Hyperparameter Method

Suppose we are dealing with multiple datasets $\{D_1,\dots ,D_N\}$, coming from a collection of different surveys $\left\{S_1,\dots ,S_N\right\}$. Sometimes it is difficult to know, a priori, if all our data are consistent with each other, or whether there could be one or more that are likely to be erroneous. If we were sure that all datasets are consistent, then it should be enough to update the probability, as seen in Section 3.2, in order to calculate the new posterior distribution for the parameters we are interested in. However, since there is usually an uncertainty about this, a way to know how useful the data may be is by introducing the hyperparameter method. This method was initially introduced by Reference [20,21] in order to perform a joint estimation of cosmological parameters from combined datasets. This method may be used as long as every survey is independent of each other. In this section, we review the main steps necessary to understand the hyperparameter method.

The main feature of this process is the introduction of a new set of hyperparameters $\alpha$ in the Bayesian procedure to allow extra freedom in the parameter estimation. These hyperparameters are equivalent to nuisance parameters in the sense that we need to marginalize over them in order to recover the posterior distribution, i.e.,

(49)$P\left(\theta \right|D,H)=\frac{1}{P\left(D\right|H)}\int P\left(\theta \right|\alpha ,H)P\left(\alpha \right|D,H)d\alpha ,$

where we have used the Bayes’ theorem. Now, for the method, it is necessary to assume that the hyperparameters $\alpha$ and the parameters of interest $\theta$ are independent, i.e., $P(\theta ,\alpha ,H)=P\left(\alpha \right)P(\theta ,H)$, and it is also necessary to assume that each hyperparameter ${\alpha }_k$ is independent of each other, i.e., $P\left(\alpha \right)=P\left({\alpha }_1\right)P\left({\alpha }_2\right)\dots P\left({\alpha }_N\right)$. In this way, we can rewrite the above expression as

$P\left(\theta \right|D,H)=\frac{P(\theta ,H)}{P\left(D\right|H)}\left[\prod _{k=1}^N\int P\left(D_k\right|\theta ,{\alpha }_k,H)P\left({\alpha }_k\right)d{\alpha }_k\right].$

Here, the quantity inside the square brackets is the marginalized likelihood over the hyperparameters. We can identify the quantity inside the integration as the individual likelihood $L\left(D_k\right|\theta ,{\alpha }_k,H)$, for every ${\alpha }_k$ and the dataset $D_k$; $P\left(D\right|H)$ is the evidence and, similarly to a parameter inference procedure, it works as a normalizing function, i.e., $P\left(D\right|H)=\int d\theta P(\theta ,H\left)L\right(D|\theta ,H)$. Notice that, by considering $P\left({\alpha }_k\right)=\delta ({\alpha }_k-1)$, we rely on the standard approach, where no hyperparameters are used.

We add these ${\alpha }_k$ in order to weight every dataset and take away the data that does not seem to be consistent with other ones. Then, we would like to know whether the data supports the introduction of hyperparameters or not. A way to address this point is given by the Bayesian evidence K defined in Equation (13). If we consider a Gaussian likelihood with maximum entropy prior, and assuming that on average the hyperparameters’ weights are unity, we can rewrite the marginalized likelihood function $L\left(D\right|\theta ,H_1)$ for model $H_1$ as

(50)$P\left(D\right|\theta ,H_1)=\prod _{k=1}^N\frac{2\mathsf{\Gamma }(\frac{n_k}{2}+1)}{{\pi }^{n_k/2}{\left|V_k\right|}^{1/2}}{({\chi }_k^2+2)}^{-\left(\frac{n_k}{2}+1\right)},$

obtaining an explicit functional form for K, given by

(51)$K=\prod _{k=1}^N\frac{2^{n_k/2+1}\mathsf{\Gamma }(n_k/2+1)}{{\chi }_k^2+2}{e}^{-{\chi }_k^2/2}.$

Here, ${\chi }_k^2$ is given by (36) for every dataset, and $n_k$ is the number of points contained in $D_k$. In Equation (50), $V_k$ is the covariance matrix for the k-data. Suppose we have two models, one with hyperparameters, called $H_1$, and a second one without them, called $H_0$. The Bayesian evidence $P\left(D\right|H_i)$ is the key quantity for making a comparison between two different models. In fact, by using the Bayes factor K from Equation (51), we can estimate the necessity to introduce the hyperparameters to our model using the criteria given in Table 2. Notice that, if we have a set of independent samples for $H_0$, we can compute an estimate for K with the help of Equation (48).

4. Numerical Tools

In typical scenarios, it is very difficult to compute the posterior distribution analytically. For these cases, the numerical tools play an important role during the parameter estimation task. There exist several options to carry out this work; nevertheless, in this section, we focus on the Markov Chain Monte Carlo (MCMC) with the Metropolis Hastings algorithm (MHA). Additionally, we present some useful details we take into account to make our computation more efficient.

4.1. MCMC Techniques for Parameter Inference

The purpose of the MCMC algorithm is to build up a sequence of points (called chain) in a parameter space in order to evaluate the posterior of Equation (9). In this section, we review the basic steps for this procedure in a simplistic way; however, it is recommendable to check [22,23,24,25,26] for a more formal version of the MCMC theory.

A Monte Carlo simulation is assigned to algorithms that use random number generators to approximate a specific quantity. On the other hand, a sequence $X_1,X_2,\dots$ of elements of some set is a Markov Chain if the conditional distribution of $X_{n+1}$ given $X_1,\dots ,X_n$ depends only on $X_n$. In other words, a Markov Chain is a process where we can compute subsequent steps based only in the information given at the present. An important property of a Markov Chain is that it converges to a stationary state where successive elements of the chain are samples from the target distribution; in our case, it converges to the posterior $P\left(\theta \right|D,H)$. Hence, we can estimate all the usual quantities of interest out of the posterior (mean, variance, etc.).

The combination of both procedures is called an MCMC. The number of points required to get good estimates in MCMCs is said to scale linearly with the number of parameters, so this method becomes much faster than grids as the dimensionality increases.

The target density is approximated by a set of delta functions:

(52)$p\left(\theta \right|D,H)\simeq \frac{1}{N}\sum _{i=1}^N\delta (\theta -{\theta }_i),$

with N being the number of points in the chain. Then, the posterior mean is computed as

(53)$⟨\theta ⟩=\int d\theta \theta P(\theta ,H|D)\simeq \frac{1}{N}\sum _{i=1}^N{\theta }_i,$

where ≃ follows because the samples ${\theta }_i$ are generated out of the posterior by construction. Then, we can estimate any integrals (such as the mean, variance, etc.) as

(54)$⟨f\left(\theta \right)⟩\simeq \frac{1}{N}\sum _{i=1}^{N}f\left({\theta }_i\right).$

As mentioned before, in a Markov Chain, it is necessary to generate a new point ${\theta }_{i+1}$ from the present point ${\theta }_i$. However, as it is expected, we need a criterion for accepting (or rejecting) this new point depending on whether it turns out to be better for our model or not. If the new step is worse than the previous one, we may still accept it since, if we only accept steps with better probability, we could be converging into a local maximum in our parameter space and, therefore, just exploring a small region of the entire space. The simplest algorithm that contains all this information in its methodology is known as the Metropolis-Hastings algorithm.

4.1.1. Metropolis-Hastings Algorithm

In the Metropolis-Hastings algorithm [27,28], it is necessary to start from a random initial point ${\theta }_i$, with an associated posterior probability $p_i=p\left({\theta }_i\right|D,H)$. We need to propose a candidate ${\theta }_c$ by drawing from a proposal distribution $q({\theta }_i,{\theta }_c)$ used as a generator of new random steps. Then, the probability of acceptance of the new point is given by

(55)$p\left(acceptance\right)=min\left[1,\frac{p_{c}q({\theta }_c,{\theta }_i)}{p_{i}q({\theta }_i,{\theta }_c)}\right].$

If the proposal distribution is symmetric, the algorithm is reduced to the Metropolis algorithm

(56)$p\left(acceptance\right)=min\left[1,\frac{p_c}{p_i}\right].$

In this way, the complete algorithm can be expressed by the following steps:

• 1.. Choose a random initial condition ${\theta }_i$ in the parameter space and compute the posterior distribution.

• 2.. Generate a new candidate from a proposal distribution in the parameter space and compute the corresponding posterior distribution.

• 3.. Accept (or reject) the new point with the help of the Metropolis-Hastings algorithm.

• 4.. If the point is rejected, repeat the previous point in the chain.

• 5.. Repeat steps 2–4 until you have a large enough chain.

4.1.2. A First Example of Parameter Inference

In order to put into practice the numerical tools, we take again the coin toss example seen from Section 3.1. Let us try to estimate the value of p (or region of values for p) that best matches our data (we assume only the 14 times that the coin was thrown). To calculate the posterior distribution (20), we use the MHA.

As mentioned before, we consider a likelihood given by the binomial distribution (18) and a normal distributed prior (11) ($a=b=1$). As our first guess for p, we consider $p_i=0.1$. We generate a new candidate $p_c$ as $p_c=p_{cu}+G(p_{cu},\widehat{\sigma })$, where $G(p_{cu},\widehat{\sigma })$ is our proposed Gaussian distribution centered at $p_{cu}$ with variance $\widehat{\sigma }=0.1$; $p_{cu}$ is the current value of p, for our first step is $p_{cu}=p_i$. Then, we introduce the MHA in a Python code. Our final result (shown in Figure 4) is a posterior distribution that matches very well with the results calculated analytically (shown in Figure 1). Numerically, we obtained $p=0.{695}_{-0.107}^{+0.123}$, where the upper and lower values for p correspond to the $1\sigma$ standard deviation. Notice that we have plotted the width of our $1\sigma$, $2\sigma$ and $3\sigma$ confidence regions in the same figure. To complement the example (and Figure 4), we also show in Figure 5 the Markov Chains generated by our code. It is easy to see that the chains oscillate with small amplitude around the mean value.

Remark: In Figure 6, we include the MCMC algorithm using an explicit code for the MCMC process. However, in Python there are some modules that can simplify this task. For example, PyMC3 [29] is a Python module that implements statistical models and fitting algorithms, including the MCMC algorithm. We use this module at the end of this section.

4.1.3. Convergence Test

It is clear that we need a test to know when our chains have converged. We need to verify that the points in the chain are not converging to a false convergent point or to a local maximum. In this sense, we need that our algorithm takes into account this possible difficulty. The simplest way (the informal way) to know whether our chain is converging to a global maximum is by running several chains starting with different initial proposals for the parameters. Then, if we see by the naked eye that all chains seem to converge into a single region of the possible value for our parameter, we may say that our chains are converging to that region.

Taking yet again the example of the coins, we can run several chains for the above example and try to estimate whether the value (region) of p that we found is a stationary value. In Figure 5, we plot 5 different Markov chains with initial conditions $p=0.1,0.3,0.5,0.7,0.9$. As we expected from the analytical result, after several steps all the chains seem to concentrate nearby the same value.

The convergence method used above is very informal, and we would like to have a better way to ensure that our result is correct. The usual test is the Gelman-Rubin convergence criterion [30,31]. That is, by starting with M chains with very different initial points and N points per chain, if ${\theta }_i^j$ is a point in the parameter space of position i and belonging to the chain j, we need to compute the mean of each chain:

(57)$⟨{\theta }^j⟩=\frac{1}{N}\sum _{i=1}^N{\theta }_i^j,$

and the mean of all the chains

(58)$⟨\theta ⟩=\frac{1}{NM}\sum _{i=1}^N\sum _{j=1}^M{\theta }_i^j.$

Then, the chain-to-chain variance B is

(59)$B=\frac{1}{M-1}\sum _{j=1}^M{(⟨{\theta }^j⟩-⟨\theta ⟩)}^2,$

and the average variance of each chain is

(60)$W=\frac{1}{M(N-1)}\sum _{i=1}^N\sum _{j=1}^M{({\theta }_i^j-⟨{\theta }^j⟩)}^2.$

If our chains converge, W and $B/N$ must agree. In fact, we say that the chains converge when the quantity

(61)$\widehat{R}=\frac{\frac{N-1}{N}W+B(1+\frac{1}{M})}{W},$

which is the ratio of the two estimates, approaches unity. A typical convergence criterion is when $0.97<\widehat{R}<1.03$.

4.1.4. Some Useful Details

The proposal distribution: The choice of a proposal distribution q is crucial for the efficient exploration of the posterior. In our example, we used a Gaussian-like distribution with a variance (step) $\widehat{\sigma }=0.1$. This value was taken because we initially explored different values for $\widehat{\sigma }$, and we select the quickest that approaches the analytic posterior distribution of p. However, if the scale of q is too small compared to the scale of the target (in the sense that the typical jump is small), then the chain may take very long to explore the target distribution, which implies that the algorithm will be very inefficient. As we can see in Figure 7 (top panel), considering an initial step $p_i=0.6$ and a variance for the proposal distribution $\widehat{\sigma }=0.002$, the number of points are not enough for the system to move to its “real” posterior distribution. On the other hand, if the scale of q is too large, the chain gets stuck, and it does not jump very frequently (bottom panel of the figure corresponding to $\widehat{\sigma }=0.8$), so we will have different maxima in our posterior.

In order to fix this issue in a more efficient way, it is recommendable to run an exploratory MCMC, compute the covariance matrix from the samples, and then re-run with this covariance matrix as the covariance of a multivariate Gaussian proposal distribution. This process can be computed a couple of times before running the final MCMC.

The burn-in: It is important to notice that, at the beginning of the chain, we have a set of points far outside the stationary region. This early stage of the chain (called burn-in) must be ignored; this means that the dependence on the starting point must be lost. Thus, it is important to have a reliable convergence test.

Thinning: There are several Bayesian statisticians that usually thin their MCMC, which means that they do not prefer to save every step given by the MCMC; instead, they prefer to save a new step each time n steps have taken place. An obvious consequence of thinning the chains is that the amount of autocorrelation is reduced. However, as long as the chains are thinned, the precision for the estimated parameters is reduced [32]. Thinning the chains can be useful in other kind of circumstances, for example, if we have limitations in memory. Notice that thinning a chain does not yield incorrect results; it yields correct results, but less efficient ones, than using the full chains.

Autocorrelation probes: A complementary way to look for convergence in an MCMC estimation is by looking for the autocorrelation between the samples. The autocorrelation $lagk$ is defined as the correlation between every sample and the sample k steps before. It can be quantified as [33,34]

(62)$\begin{array}{rrr}{\rho }_k& =& \frac{Cov(X_t,X_{t+k})}{\sqrt{Var\left(X_t\right)Var\left(X_{t+k}\right)}}\hfill \\ & =& \frac{E\left[(X_t-X)(X_{t+k}-X)\right]}{\sqrt{E\left[{(X_t-X)}^2\right]E\left[{(X_{t+k}-X)}^2\right]}},\hfill \end{array}$

where $X_t$ is the t-th sample, and X is the mean of the samples. This autocorrelation should become smaller as long as k increases (this means that samples start to become independent).

More Samplers

Gibbs sampling: The basic idea of the Gibbs sampling algorithm [35] is to split the multidimensional $\theta$ into blocks and sample each block separately, conditional on the most recent values of the other blocks. It basically breaks a high-dimensional problem into low-dimensional problems.

The algorithm reads as follows:

• 1.. $\theta$ consists of k blocks ${\theta }_1,\dots ,{\theta }_k$. Then, at step i

• 2.. Draw ${\theta }_1^{i+1}$ from $p\left({\theta }_1\right|{\theta }_2^i,\dots ,{\theta }_k^i)$

• 3.. Draw ${\theta }_2^{i+1}$ from $p\left({\theta }_2\right|{\theta }_1^{i+1},{\theta }_3^i,\dots ,{\theta }_k^i)$

• 4.. …

• 5.. Draw ${\theta }_k^{i+1}$ from $p\left({\theta }_k\right|{\theta }_1^{i+1},{\theta }_2^{i+1},\dots ,{\theta }_{k-1}^{i+1})$

• 6.. Repeat the above steps for the wished iterations with $i\to i+1$.

The distribution $p\left({\theta }_1\right|{\theta }_2,\dots ,{\theta }_k)=\frac{p({\theta }_1,\dots ,{\theta }_k)}{p({\theta }_2,\dots ,{\theta }_k)}$ is known as the full conditional distribution of ${\theta }_1$. This algorithm is a special case of MHA where the proposal is always accepted.

Metropolis Coupled Markov Chain Monte Carlo ($M{C}^3$): It is easy to see that it could be a little problematic if our likelihood has local maxima. The $M{C}^3$ is a modification of the standard MCMC algorithm that consists of running several Markov Chains in parallel to explore the target distribution for different temperatures. The temperature controls the height of the peaks in the likelihood; this simplifies the way we sample our parameter space and helps us to avoid local maxima. If a distribution has a temperature $T<1$, it is said to be tempered.

We consider a tempering version of the posterior distribution $P(\theta ,T|D,H)$:

(63)$P(\theta ,T|D,H)\propto L{(\theta ,D)}^{1/T}P(\theta ,H),$

where L is the likelihood, and $P(\theta ,H)$ the prior. Notice that, for higher T, individual peaks of L become flatter, making the distribution easier to sample with an MCMC algorithm. Now, we have to run N chains with different temperatures assigned in a ladder $T_1<T_2<\dots <T_N$, usually taken with a geometrically distributed division, with $T_1=1$. The coldest chain $T_1$ samples the posterior distribution more accurately and behaves as a typical MCMC. Then, we define this chain as the main chain. The rest of the chains are running such that they can cross local maximum likelihoods easier and transport this information to our main chain.

The chains explore independently the landscape for a certain number of generations. Then, in a pre-determined interval, the chains are allowed to swap its actual position with probability

(64)$A_{i,j}=min\left\{{\left(\frac{L\left({\theta }_i\right)}{L\left({\theta }_j\right)}\right)}^{1/T_j-1/T_i},1\right\}.$

In this way, if a swap is accepted, chains i and j must exchange their current position in the parameter space, and then chain i has to be in position ${\theta }_j$ and chain j has to move to position ${\theta }_i$.

We can see that, since the hottest chain $T_{\max }$ can access easier to all the modes of $P(\theta ,H,T_{\max }|D)$, then it can propagate its position to colder chains, to be precise, it can propagate its position to the coldest chain $T=1$. At the same time, the position of colder chains can be propagated to hotter chains, allowing them to explore the entire prior volume. For an extensive explanation, or modification to make the temperature of the chains dynamical, see Reference [36,37,38].

Affine Invariant MCMC Ensemble Sampler: The main property of this algorithm relies on its invariance under affine transformations. Let us consider a highly anisotropic density:

(65)$p(x_1,x_2)\propto exp\left(\frac{-{(x_1-x_2)}^2}{2ϵ}-\frac{{(x_1+x_2)}^2}{2}\right),$

which is difficult to calculate for small $ϵ$. But, by making the affine transformation

(66)$y_1=\frac{x_1-x_2}{\sqrt{ϵ}},y_2=x_1+x_2,$

we can rewrite the anisotropic density into the easier problem

(67)$p(y_1,y_2)\propto exp\left(\frac{-(y_1^2+y_2^2)}{2}\right).$

An MCMC sampler has the form $X(t+1)=R\left(X\right(t),\psi (t),p)$, where $X\left(t\right)$ is the sample after t iterations, R is the sampler algorithm, $\psi$ is the sequence of independent identically distributed random variables, and p is the density. A sampler is said to be affine invariant if, for any affine transformation $Ax+b$,

(68)$R(AX\left(t\right)+b,\psi \left(t\right),p_{A,b})=AR(X\left(t\right),\psi \left(t\right),p)+b.$

There are already several algorithms that are affine invariant, one of the easiest is known as the stretch move [39]. An algorithm fully implemented in Python under the name EMCEE [40] is also affine invariant, and some others that can be found in Reference [41].

Even more samplers: The generation of the elements in a Markov chain is probabilistic by construction, and it depends on the algorithm we are working with. The MHA is the easiest algorithm used in Bayesian inference. However, there are several algorithms that can help to accomplish our mission, for instance, some of the most popular and effective ones, are the Hamiltoninan Monte Carlo (see, e.g., Reference [42,43]) or the Adaptative Metropolis-Hastings (AMH) (see, e.g., Reference [19]).

5. Fitting a Straight-Line

In this section, we carry out an standard example: fitting a straight-line. That is, we assume that we have a certain theory where our measurements should follow a straight line. Then, we simulate several datasets and focus on two different cases (Figure 8):

• 1.. Consider two datasets coming from the same straight-line but having different errors.

• 2.. Consider two datasets coming from different straight-lines and also having different errors.

In our analysis, we used the PyMC3 module, and the complete code can be downloaded from the Git repository [44].

5.1. Case 1

In this example, we start by considering that our measurements for a given theory (a straight-line $y=a+bx$) are given by the data shown in the upper panel of Figure 8. The two datasets, D1 and D2, were generated from the same line $y=3+2x$, adding a gaussian error to each point. For D1, we added an error with a standard deviation ${\sigma }_1=0.3$, while, for D2 we use ${\sigma }_2=0.2$. Then, we would like to estimate the parameters of the model, i.e., a and b. We will analyze this data with and without the hyperparameter method and discuss in detail our results.

Model $H_0$: without hyperparameters.

Before we make a Bayesian estimation, it is necessary to specify the priors. As we have seen, a good prior is a non-informative one. Suppose we only know some limits for a and b, as can be seen when plotting the data. Then, we consider the flat priors

(69)$a\propto U[0,5]\mathrm{and}b\propto U[0,3],$

where $U[\alpha ,\beta ]$ are uniform distributions with lower limit $\alpha$ and upper limit $\beta$.

From Equation (30), we can write our likelihood as

(70)$L(D;line)\propto exp\left[-\sum _d\frac{{(y_d-y)}^2}{2{\sigma }_d^2}\right],$

where $y_d$ is our data taken from the dataset $D=D_1+D_2$, and ${\sigma }_d$ its errors.

We use the MHA to generate the Markov chains. In our analysis, we ran 5 chains with 10,000 steps for each one. We ran each chain with temperature $T=2$ and we thinned them every 50 steps. The result corresponds to $a=2.982\pm 0.047$ and $b=1.994\pm 0.013$, and their posterior distributions are plotted in Figure 9. Notice that there are some regions where the frequency of events in our sample is increased. So, we can say that such parameter regions seem more likely to match the data. Additionally, we compute the Gelman-Rubin criterion for each variable in order to verify that our results converged, i.e., for a is $1.000017$, and for b is $1.000291$. We see that this number is very close to 1, so our convergence criterion is achieved. The bottom panel of Figure 9 displays the $1-4\sigma$ confidence regions. We also added a point in red to show the real value for our parameters. The real value for a and b are within the curve corresponding to one standard deviation of our estimations in the inferential method.

As we mentioned, we need the autocorrelation to be small as k increases in order to consider that our analysis is converging. We see in Figure 10 such plots and notice that our convergence criterion is fulfilled. Then, in Case 1, we can see that the model $H_0$ seems to be a very good estimation procedure.

Model $H_1$: with hyperparameters.

Now, let us consider the Hyperparameter method. In this case, our likelihood can be written as Equation (50). Similar to the last procedure, we compute the posterior with flat priors, using 5 chains with 10,000 steps for each one, and check for autocorrelations. Our results are as follows: $a=2.97\pm 0.038$ with Gelman-Rubin of 1.000113, and $b=1.995\pm 0.010$ with Gelman-Rubin 1.000155. Comparing both procedures, we observe they provide similar results. In fact, the confidence regions for both approximations, Figure 9 and the top panel of Figure 11 are similar as well. So, which method is better? We could say that the method with hyperparameters is as good as the one without them, but, in order to be sure, we compute the evidence ratio K between both models. We obtained from Equation (51)

(71)$K=3.$

Then, comparing with Table 2, we can say that the evidence for $H_1$ to be better than $H_0$ is weak. In such a case, it should be equally better to work with $H_0$ as to $H_1$, as explained before.

Finally, in order to exemplify our results, let us plot in the bottom panel of Figure 11 our data with the straight-line inferred by the mean parameters of both models. As we expected, our estimation fits well the data for both cases.

5.2. Case 2

Here, we consider that we have the same theory for the straight-line but different measurements. The data points are given in the lower panel of Figure 8. These correspond to our dataset $D_1$ and $D_2$, but now changing $D_2$ by 16 new points generated around the line $y=3.5+1.5x$ with a Gaussian noise and standard deviation $\sigma =0.5$. So, our datasets are not auto-consistent with each other. Let us make again the parameter estimation for the parameters a and b and look for the differences in both procedures.

Model $H_0$: without hyperparameters.

We follow the same procedure as in Case 1. We computed our posterior and verified that our results converged with the help of the Gelman-Rubin criterion and the autocorrelation plots. Our results are the following: $a=3.528\pm 0.056$ and $b=1.795\pm 0.014$. Then, we plotted $1-4\sigma$ confidence regions in the upper panel of Figure 12. It is easy to see that our estimation differs so much from the real parameters of the datasets (red points). This is because we are trying to fit a model with non-auto-consistent datasets; therefore, we arrive at incorrect results. Now, let us see what happens in the hyperparameters procedure.

Model $H_1$: With hyperparameters.

In the bottom-left panel of Figure 12, we plotted the posterior distribution. We see immediately that both approximations are very different. While, for model $H_0$, we obtained a single region far away of the real values of our data, for model $H_1$, we obtained two local maximum regions near the real values for our datasets (red dots).

Given the fact that we know a priori the real values of the parameters, we could immediately say that the method with hyperparameters is a better approximation than the case without them. However, we confirm this assumption by computing the ratio K between both models. We then obtain

(72)$K=37,$

which means that we have very strong evidence that $H_1$ is better that $H_0$. Finally, we can plot the straight-line inferred by model $H_0$ and the two lines inferred by model $H_1$. Considering the parameters inside the two regions in the bottom-left of Figure 12, we obtain the bottom-right panel of the same Figure.

6. Bayesian Statistics in Cosmology

6.1. Theoretical Background

Bayesian statistics is a very useful tool in Cosmology to determine, for instance, the combination of model parameters that best describes the Universe. In this section, we present the basics of Cosmology and apply the Bayesian statistics to perform the parameter inference. In our examples, we will focus on the background Universe—and avoid perturbations—since the main purpose of this article is the application of these techniques rather than the cosmology by itself. It should be clear, however, that the extension to consider perturbations is immediate, i.e., there is only an increment in the number of parameters, and the expressions turn out to be just a little more complicated.

6.1.1. Einstein Field Equations

In order to specify the geometry of the Universe, an essential assumption is the Cosmological Principle: for a particular time and on sufficiently large scales, the observable Universe can be considered homogeneous and isotropic, with great precision. For example, at scales greater than 100 Mega-parsecs, the distribution of galaxies observed on the celestial sphere justifies the assumption of isotropy. The uniformity observed in the temperature distribution (one part in ${10}^5$) measured through the Cosmic Microwave Background (CMB) is the best observational evidence we have in favor of a universal isotropy. Therefore, if isotropy is taken for granted and by taking into account that our position in the Universe has no preference—known as the Copernican Principle—, the homogeneity follows when considering isotropy in each point.

Homogeneity establishes that the Universe is observed equally in each point of space.

Isotropy establishes that the Universe is observed equally in all directions.

The formalism of General Relativity establishes the relationship between the geometry of space-time and the matter on it. That is, the curvature of the spacetime produces physical effects on the matter, and these effects are associated with the gravitational field. Additionally, the curvature is related to the matter, described by an energy-momentum tensor $T_{\mu \nu }$. The above expressions can be summarized by paraphrasing Wheeler: “matter tells space-time how to curve and, in turn, the geometry of this curvature tells matter how to move”. We can write this sentence down by the Einstein equations

(73)$G_{\mu \nu }=8\pi G{T}_{\mu \nu },$

where $G_{\mu \nu }$ is the Einstein tensor (geometry of the spacetime), and G is the gravitational Newton constant [45,46,47,48]. Throughout this review, we use natural units $c=\hslash =1$.

The distance between two points in a curved space-time can be measured as

(74)$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu },$

where $g_{\mu \nu }$ is the metric tensor that contains all the information about the geometry of the space-time. From now on, and unless stated otherwise, greek letters $\mu$, $\nu$, …denote spacetime indices ranging from 0 to 3, while latin letters i, j, …denote spatial coordinates ranging from 1 to 3.

The geometry that best describes a homogeneous, isotropic, and expanding Universe is given by the Friedmann-Lemaître-Robertson-Walker metric (FLRW), with a line element

(75)$d{s}^2=d{t}^2-a^2\left(t\right){\gamma }_{ij}d{x}^{i}d{x}^j,$

where

(76)${\gamma }_{ij}\equiv {\delta }_{ij}+\kappa \frac{x_i{x}_j}{1-\kappa \left(x_k{x}^k\right)}.$

In Equation (75), a represents the scale factor of the Universe which only depends on time, and by convention is normalized to today $a\left(t_0\right)\equiv 1$. Similarly, in expression (76), $x^i$ labeled the spatial coordinates (also called comoving coordinates), ${\delta }_{ij}$ is the Kronecker delta and $\kappa$ describes the curvature of the space-time.

6.1.2. Friedmann and Continuity Equations

The content of the Universe needs to satisfy homogeneity and isotropy, as well; hence, here it is described by the energy-momentum tensor of a perfect fluid

(77)$T_{\mu \nu }=(\rho +P)U_{\mu }{U}_{\nu }-P{g}_{\mu \nu },$

where $\rho$ is the energy density, P is the fluid pressure, and $U_{\mu }$ is the 4-velocity relative to the observer. If we take the velocity as $U^{\mu }=\left(1,0,0,0\right)$ (comoving observer), the energy-momentum tensor reduces to

(78)$T_{\nu }^{\mu }=g^{\mu \lambda }{T}_{\lambda \nu }=\left(\begin{array}{cccc}\rho & 0& 0& 0\\ 0& -P& 0& 0\\ 0& 0& -P& 0\\ 0& 0& 0& -P\end{array}\right).$

Using Equations (73) and (78), with the FLRW metric, we can obtain the Friedmann equations

(79)$H^2\equiv {\left(\frac{\dot{a}}{a}\right)}^2=\frac{8\pi G}{3}\sum _i{\rho }_i-\frac{\kappa }{a^2},$

(80)$\frac{\ddot{a}}{a}=-\frac{4\pi G}{3}\sum _i\left({\rho }_i+3P_i\right).$

In these expressions, H accounts for the rate of expansion/contraction of the Universe, named as the Hubble parameter. Subindex i labels all the components that believe the Universe is made of, as we will see in the next section. These equations describe the evolution of the Universe. By combining Equations (79) and (80), we can obtain the continuity equation given by

(81)$\dot{\rho }+3\frac{\dot{a}}{a}\left(\rho +P\right)=0.$

The meaning of (81) is the conservation of the energy-momentum tensor (${\nabla }_{\mu }{T}^{\mu \nu }=0$). In order to close the system, we need to include an equation-of-state that relates pressure and energy density for a given fluid. In particular, we are interested on barotropic fluids, which generally have the form of $P=\omega \rho$.

6.1.3. Content of the Universe

Once the equations that define the dynamics of the Universe are known, it is necessary to specify its content. The standard cosmological model, also known as $\mathsf{\Lambda }$ Cold Dark Matter ($\mathsf{\Lambda }$CDM), is one of the most accepted models to describe the Universe, with its content being:

• Dust: It has no pressure, and its energy density takes the form of $\rho \propto a^{-3}.$ Dust is conformed by baryons (ordinary matter).

• Dark matter: It is proposed to explain several astrophysical observations, like the dynamics of galaxies in the Coma cluster or the rotation curves of galaxies [49,50]. The $\mathsf{\Lambda }$CDM model assumes the dark matter only interacts gravitationally (and possibly by weak interaction) with the rest of the Universe, hence its name, Cold Dark Matter (CDM). Since it is proposed as interacting only via gravitational force, there can be several candidates to fulfill this requirement: it could be conformed by weakly interacting massive particles (WIMPs), by gravitationally-interacting massive particles (GIMPs); by axions (hypothetical elementary particles); or by sterile neutrinos, just to mention a few. For a short review about Dark Matter and possible candidates, see Reference [51].

• Radiation: This corresponds to relativistic particles that follow the relation $P=\frac{1}{3}\rho$, which implies a density with a behavior $\rho \propto a^{-4}$. We consider photons ${\rho }_{\gamma }$ and massless neutrinos ${\rho }_{\nu }$ as radiation, so the total radiation energy density in the Universe is given by

  (82)${\rho }_r={\rho }_{\gamma }+{\rho }_{\nu }.$

  The relation between these quantities is

  (83)${\rho }_{\nu }=N_{\mathrm{eff}}\times \frac{7}{8}\times {\left(\frac{4}{11}\right)}^{4/3}{\rho }_{\gamma },$

  where $N_{\mathrm{eff}}$ is the effective number of relativistic degrees of freedom, with standard value $N_{\mathrm{eff}}=3.046$ [13].

• Dark Energy: It is introduced to explain the current accelerated expansion of the Universe. In the $\mathsf{\Lambda }$CDM model, dark energy is given by the cosmological constant $\mathsf{\Lambda }$ or equivalently by an equation-of-state $\omega =-1$.

Each of these components can be described by its equation of state shown in Table 4, and defining the density parameter

(84)${\mathsf{\Omega }}_i\equiv \frac{{\rho }_i}{{\rho }_{\mathrm{crit}}},\mathrm{with}{\rho }_{\mathrm{crit}}=\frac{3H^2}{8\pi G},$

with ${\rho }_{\mathrm{crit}}$ being the condition to have a flat Universe or equivalently zero curvature, we can rewrite (79) as

(85)$\frac{H^2}{H_0^2}={\mathsf{\Omega }}_{r,0}{a}^{-4}+{\mathsf{\Omega }}_{m,0}{a}^{-3}+{\mathsf{\Omega }}_{k,0}{a}^{-2}+{\mathsf{\Omega }}_{\mathsf{\Lambda },0},$

where ${\mathsf{\Omega }}_{r,0}$ is the radiation density parameter, ${\mathsf{\Omega }}_{m,0}\equiv {\mathsf{\Omega }}_{b,0}+{\mathsf{\Omega }}_{DM,0}$ corresponds to the total matter, ${\mathsf{\Omega }}_{b,0}$ to baryons, ${\mathsf{\Omega }}_{DM,0}$ to dark matter, ${\mathsf{\Omega }}_k\equiv -\kappa /{\left(aH\right)}^2$ the curvature density parameter, and ${\mathsf{\Omega }}_{\mathsf{\Lambda }}\equiv \mathsf{\Lambda }/3H^2$ associated with the Cosmological Constant, and the subscript zero indicates they are evaluated by today’s values.

6.1.4. Alternatives to the $\mathsf{\Lambda }$CDM Model

The $\mathsf{\Lambda }$CDM model has had great success in modeling a wide range of astronomical observations. However, it is in apparent conflict with some observations on small-scales within galaxies (e.g., cuspy halo density profiles, overproduction of satellite dwarfs within the Local Group, amongst many others; see, for example, Reference [50,52]). In addition, all attempts to detect WIMPs either directly in the laboratory, or indirectly by astronomical signals of distant objects have failed so far. For some of these reasons, it seems necessary to explore alternatives to the standard $\mathsf{\Lambda }$CDM model. With this in mind, several alternatives have been suggested. For instance, the Scalar Field Dark Matter (SFDM) model proposes the dark matter as a spin 0 boson particle [53,54,55,56,57,58]; or the Self Interacting Dark Matter, as its name states, relies on the cold dark matter to be made of self interacting particles [59,60,61]. On the other hand, in order to explain the accelerated expansion of the Universe, there exist different modifications to the theory of General Relativity, i.e., $f\left(R\right)$ theories [62], braneworld models [63,64]. There are also some alternatives to the cosmological constant as Dark Energy, i.e., scalar fields (quintessence, K-essence, phantom, quintom, non-minimally coupled scalar fields [65,66,67,68,69]); or many more alternatives, i.e., anisotropic Universes [70,71,72]. Finally, if the dark energy is assumed to be a perfect fluid, one of the most popular time-evolving parameterizations for its equation of state consists of expanding $\omega$ in a Taylor series, for example, the Chevallier-Polarski-Linder (CPL) $\omega ={\omega }_0+{\omega }_a\left(1-a\right)$, with two free parameters ${\omega }_0,{\omega }_a$ [73,74]. It may also be expanded into Fourier series [75], or many more Bayesian approaches, as have been suggested to account for a dynamical dark energy [76,77,78].

6.2. Cosmological Parameters

6.2.1. Base Parameters

These parameters, also known as standard parameters, are the main quantities used in the description of the Universe. They are not predicted by a fundamental theory, but their values must be fitted to provide the best description of the current astrophysical and cosmological observables. To explain the homogeneous and isotropic Universe, we can use the density parameter of each component ${\mathsf{\Omega }}_{i,0}$ and the Hubble parameter $H_0$ related by (85). In particular, the radiation contribution is measured with great precision, so that ${\mathsf{\Omega }}_{\gamma }$ is pinned down very accurately, and, hence, there is no need to fit this parameter. Similarly, neutrinos, as long as they maintain a relativistic behavior, can be related to the density of the photons through (83).

On the other hand, the existence of strong degeneracies from different combinations of parameters is also notorious. In particular, the geometric degeneracy involving ${\mathsf{\Omega }}_m$, ${\mathsf{\Omega }}_{\mathsf{\Lambda }}$ and the curvature parameter ${\mathsf{\Omega }}_k=1-{\mathsf{\Omega }}_m-{\mathsf{\Omega }}_{\mathsf{\Lambda }}$. To reduce these degeneracies, it is common to introduce a combination of cosmological parameters such that they have orthogonal effects in the measurements.

6.2.2. Derived Parameters

The above standard set of parameters provides an adequate description of the cosmological models. However, this parameterization is not unique, and some others can be as good as this one. Various parameterizations make use of the knowledge of the physics or the sensitivity of the detectors and can, therefore, be interpreted more naturally. In general, other parameters could have been used to describe the Universe, for example: the age of the Universe, the current temperature of the neutrino background, the epoch of equality of matter-radiation, or the epoch of reionization. In the standard cosmological model ($\mathsf{\Lambda }$CDM), in order to decrease degeneracies, the physical energy densities ${\mathsf{\Omega }}_{DM,0}{h}^2$ and ${\mathsf{\Omega }}_{b,0}{h}^2$ are used as base parameters [13].

6.3. Cosmological Observations

In this section, we review some of the most common experiments and observables used to constrain the cosmological models on the background level.

Baryon Acoustic Oscillations (BAO): The BAO is a statistical property, a feature in the correlation function of galaxies or in the power spectrum. The best description of the early Universe considers that it was made of plasma of coupled photons and matter (baryons and dark matter). The interaction between the gravitational force due to matter and the radiation pressure formed spherical waves in the plasma. When the Universe cooled down enough, the protons and electrons were able to join together, forming hydrogen atoms; therefore, this process allowed photons to decouple from the rest of the baryons. The photons began to travel uninterrupted, while the gravitational field attracted matter towards the center of the spherical wave. The final configuration is an overdensity of matter in the center and a shell of baryons of fixed radius called sound horizon. This radius, used as a standard ruler, is the maximum distance that sound waves could have traveled through the primordial plasma before recombination. The sound horizon $r_d$ is given by

(86)$r_d={\int }_{z_d}^{\infty }\frac{c_s\left(z\right)}{H\left(z\right)}dz,$

where the sound speed (in terms of redshift z) in the photon-baryon fluid is $c_s\left(z\right)=3^{-1/2}c{\left[1+\frac{3}{4}{\rho }_b\left(z\right)/{\rho }_{\gamma }\left(z\right)\right]}^{-1/2}$, and $z_d$ is the redshift when photons and baryons decouple.

The BAO scale is determined by adopting a fiducial model to be able to translate the angular and redshift separations at comoving distances. The information of the measurement is found in the ratio ($\alpha$) of the measured BAO scale and that predicted by the fiducial model ($fid$). In an anisotropic fit, two ratios are used, one perpendicular ${\alpha }_{\perp }$ and one parallel ${\alpha }_{\Vert }$ to the line of sight. A measurement of ${\alpha }_{\perp }$ constrains the ratio of the comoving angular diameter distance to the sound horizon [79]:

(87)$\frac{D_M\left(z\right)}{r_d}={\alpha }_{\perp }\frac{D_{M,fid}\left(z\right)}{r_{d,fid}},$

where the comoving angular diameter distance is given by

(88)$D_M\left(z\right)=\frac{c}{H_0}{S}_k\left(\frac{D_c\left(z\right)}{c/H_0}\right).$

The line-of-sight comoving distance is defined as

(89)$D_c\left(z\right)=\frac{c}{H_0}{\int }_0^{z}d{z}^{\prime }\frac{H_0}{H\left(z^{\prime }\right)},$

and $S_k\left(z\right)$

(90)$S_k\left(x\right)=\left\{\begin{array}{cc}sinh\left(\sqrt{{\mathsf{\Omega }}_k}x\right)/\sqrt{{\mathsf{\Omega }}_k}\hfill & {\mathsf{\Omega }}_k>0,\hfill \\ x\hfill & {\mathsf{\Omega }}_k=0,\hfill \\ sin\left(\sqrt{-{\mathsf{\Omega }}_k}x\right)/\sqrt{-{\mathsf{\Omega }}_k}\hfill & {\mathsf{\Omega }}_k<0.\hfill \end{array}\right.$

The Hubble parameter can be constrained by measuring ${\alpha }_{\Vert }$ using an analogous quantity

(91)$\frac{D_H\left(z\right)}{r_d}={\alpha }_{\Vert }\frac{D_{H,fid}\left(z\right)}{r_{d,fid}},$

with $D_H\left(z\right)=c/H\left(z\right)$.

If redshift-space distortions are weak (which is valid for luminous galaxy surveys but not for the Ly-$\alpha$), an isotropic analysis measures an effective combination of (87) and (91), and the volume averaged distance $D_V\left(z\right)$ [79]

(92)$\frac{D_V\left(z\right)}{r_d}=\alpha \frac{D_{V,fid}\left(z\right)}{r_{d,fid}},$

with $D_V\left(z\right)={\left[z{D}_H\left(z\right)D_M^2\left(z\right)\right]}^{1/3}$.

The BAO measurements constrain the cosmological parameters through the radius of the sound horizon $r_d$, Hubble distance $D_H\left(z\right)$ and the comoving angular diameter distance $D_M\left(z\right)$; see Figure 13. The data used in this work to constrain $D_V/r_d$ is obtained from the 6dF Galaxy Survey (6dFGS [80]) from UK Schmidt Telescope, the Main Galaxy Sample (MGS [81]), and the BOSS LOWZ Sample [82] from SDSS. On the other hand, the data from BOSS CMASS Sample [82], BOSS Lyman $\alpha$ auto-correlation (Ly$\alpha$ auto, [83]), and BOSS Lyman-$\alpha$ cross correlation (Ly$\alpha$ cross [84]) are used to constrain $D_M/r_d$, $D_H/r_d$.

Supernovae type Ia (SNIa): This type of supernova occurs in binary star systems, one of which is a white dwarf that accretes matter from the star that accompanies it. When the white dwarf accumulates sufficient mass (≈1.4 solar masses), its core will start the ignition temperature for carbon fusion, and, within a few seconds, it releases enough energy to produce the supernova [86]. Since type Ia supernovae (SNIa) are hypothesized to occur near the same mass limit of 1.4 solar masses, commonly referred to as Chandrasekhar mass, their luminosity peaks are fairly consistent and can be standardized and, thus, be used as standard candles [13]. From several analyses of SNIa, the Supernova Cosmology Project and High-z Supernova Search Team both found evidence that the Universe is currently expanding at an accelerated rate [87,88,89].

These stars allow us to measure relative distances using the luminosity distance given by

(93)$D_L\equiv \sqrt{\frac{L}{4\pi S}},$

where L is the luminosity defined as the energy emitted per unit solid angle per second, and S is the radiation flux density defined as the energy received per unit area per second [49]. The observable quantity is the radiation flux density received, and it cannot be translated into the luminosity density unless the absolute luminosity of the object is known. Even if the luminosity is unknown, it will appear as a scaling factor [49]. The relation between $D_L$ and the cosmological parameters is given by

(94)$D_L=D_M(1+z),$

where $D_M$ is provided by Equation (88). Another important quantity in the observation of supernovae is the standardized distance modulus

(95)$\mu =m_B^{*}-M_B+\alpha X_1-\beta C,$

where $m_B^{*}$ is the observed peak magnitude in the rest frame of blue band (B), $\alpha$, $\beta$, and $M_B$ are parameters that depend on host galaxy properties [90]. $X_1$ is the time stretching of the light curve, and C is the supernova color at maximum brightness. The relation between the standardized distance modulus and the luminosity distance is

(96)$\mu =5lo{g}_{10}\left(\frac{D_L}{10pc}\right).$

The data used in this paper (shown in Figure 14) is obtained from the Joint Light-curve Analysis (JLA). It is a collaboration to analyze the data of 740 stars from the SDSS-II (previous version of BOSS), the SuperNova Legacy Survey (SNLS [91]) experiment that used the Canada-France-Hawaii telescope (CFHT), the Calán/Totolo Survey, the Carnegie Supernova Project, the Harvard-Smithsonian Center for Astrophysics, La Silla Observatory, Fred Lawrence Observatory, and the Hubble Space Telescope (HST) [90]. For simplicity, we compress all the information into a linear function fit over 30 bins (31 nodes), spaced evenly in $log\left(z\right)$ with a $31\times 31$ covariance matrix [79].

Although SNIa have been widely used to restrict cosmological models, there are still discussions about the way it is done. That is, in order to extract information from them, an a priori cosmological model has to be assumed, which may biased the outcomes [92].

Cosmic Microwave Background (CMB): Corresponds to the radiation that permeates all the Universe, discovered in 1965. Before recombination, baryons and photons were tightly coupled, and once photons decouple from the rest of the matter, they traveled uninterrupted until reach us. The temperature radiation measured at different parts of the sky contains information of the last scattering epoch, gravitational lensing, among others. Here, the CMB displays the primordial anisotropies studied in the angular power spectrum. One of the most important recent collaborations that studies the CMB corresponds to the Planck satellite, and previous probes include COBE [94] and WMAP [95]. It is a European Space Agency mission, in which the main objective is to measure the temperature, polarization, and anisotropies of the CMB over the entire sky. These results would allow to determine the properties of the Universe at large scales and the nature of dark matter and dark energy, as well as to test inflationary theories, determine whether the Universe is homogeneous or not, and obtain maps of galaxies in the microwave [96,97,98,99,100,101,102,103].

In this work, we use the CMB information as a BAO located at redshift $z=1090$, measuring the angular scale of the sound horizon $D_M\left(1090\right)/r_d$, in addition, to calibrating the absolute length of the BAO ruler through the determination of ${\mathsf{\Omega }}_b{h}^2$ and ${\mathsf{\Omega }}_{cb}{h}^2$ (the density parameters of baryonic and dark matter, respectively; more details about these parameters can be found in Reference [79]).

Large Red Galaxies: (Cosmic Chronometers) These are the most massive galaxies for each redshift z and contain the oldest star population. These kinds of galaxies are used to estimate the Hubble factor because they contain little stardust, which makes it easier to get their spectra. The way these chronometers work is by selecting two galaxies at different redshifts between $z\sim 0-2$ and compare their upper cut in its age distributions. By doing this, it is possible to obtain the difference of ages $\mathsf{\Delta }t$ and redshifts $\mathsf{\Delta }z$ such that the expression $\frac{dz}{dt}$ can be approximated as $\frac{dz}{dt}\simeq \frac{\mathsf{\Delta }z}{\mathsf{\Delta }t}$. This quantity is related with the Hubble factor via [104,105]

(97)$H\left(z\right)=\frac{-1}{1+z}\frac{dz}{dt}\simeq \frac{-1}{1+z}\frac{\mathsf{\Delta }z}{\mathsf{\Delta }t}.$

In this work, the data used to constrain $H\left(z\right)$ was obtained from References [106,107,108,109]. In Figure 15, a compilation of these data is shown.

Another important probe in the foreseen future corresponds to the Gravitational waves from merging black holes (standard sirens), as they directly measure the luminosity distance to the merger. Therefore, it is possible to constrain the distance-redshift relation and, hence, the cosmological parameters [110,111].

Even though this paper is focused on cosmological constraints on the background level, there exist several observations to pin down the parameters at the perturbation level, and some of them include gravitational lensing [112,113,114,115,116], clustering [117,118,119,120], the matter power spectrum [121,122], and primordial gravitational waves and the primordial power spectrum, both generated during inflation [123,124,125], among many others.

7. Parameters Inference in Cosmology

Here, we estimate the Hubble parameter $H_0$ at the present time and the matter density of the Universe ${\mathsf{\Omega }}_m$, assuming a $\mathsf{\Lambda }$CDM standard model. We use our own Python code, which can be found in Reference [44].

The data—For this particular example, we focus only on the Cosmic Chronometers, as shown in Figure 15.

The model—Our interest is to fit the density parameters for each component of the Universe, as well as the value of the Hubble parameter at present time. For this purpose, we use Equation (85), and, for simplicity, we assume a flat Universe with a well measured radiation content. Then, the model is given by

(98)$H^2=H_0^2[{\mathsf{\Omega }}_{m,0}{(1+z)}^3+{\mathsf{\Omega }}_{\mathsf{\Lambda },0}],$

constrained by the relation

(99)${\mathsf{\Omega }}_{m,0}+{\mathsf{\Omega }}_{\mathsf{\Lambda },0}=1.$

Notice that the above relation implies that we can get rid of a parameter—${\mathsf{\Omega }}_{m,0}$ or ${\mathsf{\Omega }}_{\mathsf{\Lambda },0}$—in the analysis. Then, the parameters we decided to estimate are $H_0$ and ${\mathsf{\Omega }}_{m,0}$. By assuming the Gaussian approximation, we construct the likelihood given by

(100)$L\propto exp\left[-\sum _i\frac{{(H_i-H\left(z_i\right))}^2}{2{\sigma }_i^2}\right],$

where $H\left(z_i\right)$ is described by Equation (98) evaluated at each redshift $z_i$, and $H_i$ is the value of the Hubble parameter measured at $z_i$ and ${\sigma }_i$ the error of the i-th measurement.

The only a priori information we have about the free parameters is that each component of the Universe ${\mathsf{\Omega }}_{i,0}$ must satisfy the relation $0\le {\mathsf{\Omega }}_{i,0}\le 1$, while, for the present Hubble parameter, its conservative prior can be obtained by observing the data at our disposal. In such cases, a good prior choice is a uniform distribution (flat prior) with limits ${\mathsf{\Omega }}_m\in [0,1]$ and $H_0\in [10,100]$. Hence, the priors

(101)${\mathsf{\Omega }}_{m,0}\sim U[0,1],H_0\sim U[10,100].$

Numerical estimation—We follow the same procedure as in the straight-line example. In the left panel of Figure 16, we have plotted the chains obtained for our estimations and its corresponding 1D posterior distribution. Similar to the above examples, we have also plotted the 2D posterior distributions with 1–$4\sigma$ confidence regions in the right panel of the same figure. Additionally, we have obtained the mean, standard deviation and the Gelman Rubin criteria for each parameter: $H_0=67.77\pm 3.13$ with convergence 1.00045, and ${\mathsf{\Omega }}_m=0.331\pm 0.0628$ with convergence 1.00044. We can see that our estimations are very similar to the values reported in the literature [85].

7.1. Cosmological and Statistical Codes

As we have shown, until now, during the process of parameter inference for a given model, there are three steps we have followed: first, obtain the data we would like to confront with the model parameters. Then, construct the likelihood associated with the theory we are working with. Depending on the nature of the data, the likelihood can depend on different ways of the parameters. For example, in the last exercise, the likelihood in terms of the parameters is via (98) and (100). Finally, it is necessary to program the numerical tools in order to obtain the parameter inference. Such programming can be done, for example, in PyMC3, as we saw before.

It is notorious that the above process can be a tremendous task for programming, for instance, when the theories involve numerous parameters, or some of them must be marginalized in order to ignore the non-interesting ones (like nuisance parameters), or/and when the models depend on the parameters of interest in a difficult way (like by solving differential equations, integrals, etc.). Then, we can proceed in two different ways: the first is by accepting the challenge and creating our own code. Developing a new code may be a good option rather than using others, as it would require a long time to learn the implementation and modification (specially if the theory we are dealing with is somehow simple). Otherwise, we need to rely on existing codes, principally when the theory of interest is quite complicated, as it is usually the case when perturbations are taken into account on the cosmological models. In this section, we present some cosmological and statistical codes to test cosmological models.

7.1.1. Cosmological Codes

Today, there are several cosmological Boltzmann codes available. Some of them include: CMBFAST (written in FORTRAN 77 [126,127,128]), CMBEASY (C++ [129]), CAMB (FORTRAN 90 [130,131]), CLASS (C [132,133]), and COSMOSIS [134] (written in Python, and it works as an interface between CLASS, CAMB, MontePython, CosmoMC, and more). All of them are used for calculating the linear CMB anisotropy spectra, based on integrations over the sources along the photon line of sight.

7.1.2. Statistical Codes

Once the cosmological model is established, we need a statistical code to estimate the free parameters of our model. There are several MCMC codes that can make this task easy to handle. Some of them are:

Monte Python [135,136,137]—It is a Monte Carlo code for Cosmological Parameter extraction that contains likelihoods of most recent experiments, and interfaces with the Boltzmann code Class for computing the cosmological observables. The code has several sampling methods available: Metropolis-Hastings, Nested Sampling (through MultiNest), EMCEE (through CosmoHammer), and Importance Sampling.

CosmoMC [138]—This is a Fortran MCMC engine for exploring cosmological parameter space. It contains Monte Carlo samples and importance sampling. It also has by default several likelihoods of the most recent experiments, and interfaces with CAMB.

SimpleMC—It is an MCMC code for cosmological parameter estimation where only the expansion history of the Universe matters. It was written by Anẑe Slosar and José A. Vázquez, initially released in Reference [79], and can be downloaded in Reference [139]. This code solves the cosmological equations for the background parameters in the same way as CLASS or CAMB, and it contains the statistical parameter inference of CosmoMC/MontePython. An advantage of this code is that it is completely written in Python with an interface to machine learning tools, such as artificial neural networks, genetic algorithms, as well as algorithms to compute the Bayesian evidence, i.e., Dynesty [140] or MCEvidence [141].

The main idea of MCEvidence is that, asymptotically, the number density n of the MCMC is proportional to the density of the likelihood multiplied by the prior, that is, the non-normalized posterior

$\tilde{P}\left(D\right|\theta ,H)=P\left(D\right|\theta ,H)P\left(\theta \right|H)=an\left(D\right|\theta ,H),$

where $n\left(D\right|\theta ,H)=NP(D|\theta ,H)$ with N the length of the chain. The code uses Bayesian inference to find a, the constant of proportionality (see Reference [142] for details). Once this constant has been found, the Bayesian evidence is given by $P\left(D\right|H)=aN.$

8. Examples with SimpleMC

The main interest of this section is to test several cosmological models through SimpleMC. Even though we selected this code, the results we present here can also be obtained in any of the aforementioned codes.

Throughout these examples, we consider Gaussian likelihoods for each dataset with the following form:

(102)$L\propto exp\left[-\sum _i\frac{{(T\left(z_i\right)-T_i)}^2}{2{\sigma }_i}\right],$

where $T\left(z_i\right)$ is the theoretical value related to the observation $T_i$; and ${\sigma }_i$ are the corresponding errors associated to each measurement. In our estimation, $T\left(z_i\right)$ is given by (94) for Supernovae; (87) for CMB; (87), (91), and (92) for BAO; and (97) for Cosmic Chronometers.

We use the BAO data mentioned in Section 6.3 (labeled as BBAO); for Supernovae, we use the Joint Light-Curve Analysis compressed data denoted as SN, the Planck data (denoted as Planck) for CMB, and the Cosmic Chronometers data (HD).

Models of the Universe

The base example corresponds to the flat $\mathsf{\Lambda }$CDM model, already explored in Section 7, but now including the full observations presented in Section 6.3 and using the SimpleMC code. Here, in order to test the Friedmann Equation (85) with the data, we consider as free parameters (along with flat priors): the total matter dimensionless density parameter ${\mathsf{\Omega }}_m\in [0.05,1.5]$, the baryon physical density ${\mathsf{\Omega }}_b{h}^2\in [0.02,0.025]$, and the dimensionless Hubble constant $h\in [0.4,1]$. Then, assuming the base $\mathsf{\Lambda }$CDM model, we let the curvature of the Universe be a free parameter (model o$\mathsf{\Lambda }$CDM), with its corresponding flat prior ${\mathsf{\Omega }}_k\in [-1.5,1.5]$. Moreover, because the cosmological constant is only a particular case for the dark energy equation-of-state $\omega =-1$, we let $\omega$ be a free parameter with flat priors $\omega \in [-2.0,0.0]$ and labeled it as model $\omega$CDM. We may combine the addition of curvature and constant $\omega$ to define the o$\omega$CDM model. In order to go even further and describe a dynamical dark energy, we use the CPL parameterization for the equation of state with flat priors on the parameters ${\omega }_0\in [-2.0,0.0]$ and ${\omega }_a\in [-2.0,2.0]$ and labeled it as ${\omega }_0{\omega }_a$CDM model. Again, we can incorporate the curvature of the Universe to the CPL parameterization, named as o${\omega }_0{\omega }_a$CDM.

By using the combined dataset BBAO+Planck+SN+HD, Table 5 shows the best fit values, along with 1$\sigma$ confidence levels. The first important result to highlight is how the constraints have shrunk once new information is taken into account. That is, in Section 7, for the $\mathsf{\Lambda }$CDM model and by using only Hubble data, we had $h=0.677\pm 0.313$ and ${\mathsf{\Omega }}_m=0.331\pm 0.0628$. Now, with the inclusion of BAO, Planck, and SN data, the constraints have improved considerably to $h=0.684\pm 0.006$ and ${\mathsf{\Omega }}_m=0.299\pm 0.007$. Figure 16 is updated with new data in order to get Figure 17. Here, the upper panel displays the chain for the parameter $H_0=100·h$, with 9000 steps. In the lower panel of the same figure, we plot the 2D posterior distribution, along with 1 and 2$\sigma$ confidence regions obtained from our estimations.

From Table 5, we also observe the best fit of most of the new parameters—additional to the base $\mathsf{\Lambda }$CDM model—remained well inside the 1$\sigma$ confidence level. However, the exception is for the o${\omega }_0{\omega }_a$CDM model, where ${\mathsf{\Omega }}_k$, ${\omega }_0$, and ${\omega }_a$ lay down right outside the 1$\sigma$ region. This main feature is better observed in Figure 18, where the 2D posterior distributions, along with 1 and 2$\sigma$ confidence regions, are shown. Here, the standard $\mathsf{\Lambda }$CDM values are marked with a dash line ( ${\mathsf{\Omega }}_k=0$, $\omega ={\omega }_0=-1$, ${\omega }_a=0$). The inclusion of extra parameters improves the fit to the data, observed through the minimum ${\chi }_{\min }^2$, before the last row of Table 5. However, it also carries out a penalization factor that affects directly the model selection, as seen in Reference [75,77,143].

In the last row of the same table, we show the Bayes factor for the extended models using the LCDM model as reference, where the Bayesian evidences were obtained using the MCEvidence code. We found significant evidence in favor of $\mathsf{\Lambda }$CDM compared to $\omega$CDM and ${\omega }_0{\omega }_a$CDM; strong evidence with respect to o$\mathsf{\Lambda }$CDM and o${\omega }_0{\omega }_a$CDM; and decisive evidence against o$\omega$CDM. This is in agreement with the results presented in Reference [144], where Planck 2013 found no evidence against $\mathsf{\Lambda }$CDM.

9. Conclusions and Discussion

The fact that the number of cosmological observations have increased impressively over the last decade allowed to obtain a better description of the Universe. However, since we still have the limitation of a unique Universe, a Frequentist approach may not be the best to rely on; hence, the Bayesian statistics came into consideration. In this work, we provide a review of the Bayesian statistics and present some of its applications to cosmology, mainly throughout several examples.

The Bayesian statistics rests on the rules of probability which yield to the Bayes’ theorem. Given a model or hypothesis H for some data D, Bayes’ theorem tells us how to determine the probability distribution of the set of parameters $\theta$. Bayes’ theorem states that

$P\left(\theta \right|D,H)=\frac{P\left(D\right|\theta ,H\left)P\right(\theta \left|H\right)}{P\left(D\right|H)},$

where the prior probability $P\left(\theta \right|H)$—the state of knowledge before acquiring the data—is upgraded through the likelihood $P\left(D\right|\theta ,H)$ when experimental data D are considered. The aim for parameter estimation is to then obtain the posterior probability $P\left(\theta \right|D,H)$ which represents the state of knowledge once we have taken into account the new information.

We noticed that, if the prior probability is constant, we can identify the posterior probability with the likelihood $P\left(\theta \right|D,H)\propto L(D|\theta ,H)$; thus, by maximizing it, we can find the most probable set of parameters for a model given the data. Moreover, if we assume a Gaussian approximation for the likelihood, then the chi-squared quantity is related to the Gaussian likelihood via $L=L_0{e}^{-{\chi }^2/2}$. Therefore, maximizing the Gaussian likelihood is equivalent to minimizing the chi-squared. Once the posterior distribution for a set of parameters, of a given model, is calculated, we show the results in the form of confidence regions of said parameters. In addition, for this particular case, in which likelihoods are Gaussian, Fisher’s matrix can be computed according to the Hessian matrix, where the latter contains information about the errors of the parameters and their covariances. The Fisher matrix gives information about the accuracy of the model and allows predicting how well an experiment will be able to constrain the set of parameters for a given model.

On the other hand, sometimes it is difficult to know, a priori, if multiple datasets are consistent with each other, or whether there could be one or more that are likely to be erroneous. Since there is usually this uncertainty, a way to know how useful a dataset is may be by introducing the hyperparameter method. These hyperparameters act as weights for every dataset in order to take away data that does not seem to be consistent with all of them. Here, the key quantity to estimate the necessity of introducing hyperparameters to our model is given by the Bayesian evidence $P\left(D\right|H)$.

The estimation of the posterior distribution is a very computationally demanding process, since it requires a multidimensional exploration of the likelihood and prior. To carry out the exploration of the cosmological parameter space, we focus on Markov Chain Monte Carlo methods with the Metropolis Hastings algorithm. A Markov process is a stochastic process (that aims to describe the temporal evolution of some random phenomenon) where the probability distribution of the immediate future state depends only on the present state. Any computational algorithm that uses random numbers is called Monte Carlo. Thus, the Metropolis Hastings algorithm uses a transition kernel to construct a sequence of points (called chain) in the parameter space in order to evaluate the posterior distribution of said parameter. The generation of the elements in a Markov chain is probabilistic by construction, and it depends on the algorithm we are working with. The MHA is the easiest algorithm used in Bayesian inference; however, to explore complex posterior distributions more efficiently, we provide a brief description of several samplers.

Finally, we show how Bayesian statistics is a very useful tool in Cosmology to determine, for instance, the combination of model parameters that best describes the Universe. In particular, we confront the standard cosmological model ($\mathsf{\Lambda }$CDM) to current observations and compare it to different models. We found the model that best fit the data, through ${\chi }_{\min }^2$, corresponds to a curved Universe with a dynamical dark energy, namely o${\omega }_o{\omega }_a$CDM. It is important to clarify that this does not mean that the o${\omega }_o{\omega }_a$CDM model is the final model; it merely shows that, for these particular datasets, there is an improvement in the fit when compared to $\mathsf{\Lambda }$CDM, but further analysis and more data is necessary to give a verdict. However, adding extra parameters brings up a penalized factor seen through the Bayesian evidence; hence, the favored model becomes $\mathsf{\Lambda }$CDM.

Author Contributions

All authors contributed equally. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

JAV acknowledges the support provided by FOSEC SEP-CONACYT Investigación Básica A1-S-21925, FORDECYT-PRONACES-CONACYT/304001/2020, and UNAM-DGAPA-PAPIIT IA104221. LEP, LOT, and LAE were supported by CONACyT México. LEP also acknowledges sponsorship from CONACyT through grant CB-2016-282569.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Note

1. These rules are defined for a continuous variable; however, the corresponding discrete definition can be given immediately by replacing $\int dx\to \sum$.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. The coin example: blue figure displays the prior distribution P(p) which is updated, when the data is taken into account, to get the posterior distribution P(p|D) (red). The vertical black line corresponds to the real value, p=0.5.

Figure 2. Posterior distributions P(p|D), when: (a) there is no available data, (b) after 14 coin tosses of which 10 were heads, (c) after 100 coin tosses of which 56 were heads, (d) after 500 tosses of which 246 were heads. Notice that, while we continue increasing the experimental results, the posterior distribution starts to be more localized near by the real value p=0.5.

Figure 3. Converging views in Bayesian inference (taken from Reference [5]). A and B have different priors P(θ|Ii) for a value θ (panel (a)). Then, they observe one datum with an apparatus subject to a Gaussian noise and they obtained a likelihood L(θ;HI) (panel (b)), after which their posteriors P(θ|m1) are obtained (panel (c)). After observing 100 data, it can be seen that both posteriors are practically indistinguishable (panel (d)).

Figure 4. One-dimensional posterior distribution for our example. We plot the prior distribution (red), true posterior (dashed-black), and the posterior calculated by the MHA (blue). We plot 1, 2 and 3σ confidence regions for the estimation of p.

Figure 5. The five Markov Chains used to estimate the posterior distribution. We use p=0.1,0.3,0.5,0.7,0.9 for starting points.

Figure 6. Here, we show our MCMC code, written in Python, from scratch.

Figure 7. Two Markov Chains considering different variance for our Gaussian proposal distribution. The upper panel corresponds to σ^=0.002, while the lower panel corresponds to σ^=0.8.

Figure 8. Datasets D1 and D2 measured by the straight-line theory. Top: case 1. Bottom: case 2.

Figure 9. Top panel: One-dimensional marginalized posterior distributions for our samples and the Markov chains for model H0. Bottom panel: Two-dimensional marginalized posterior distributions, along with 1–4 confidence regions for our parameters for model H0. The red point corresponds to the true value.

Figure 10. Autocorrelation plots for model H0.

Figure 11. Upper panel: Confidence regions for the parameters in model H1. Lower panel: The best straight-lines inferred, along with the data.

Figure 12. Upper panel: Confidence regions for the parameters in model H0. Bottom-left panel: Confidence regions for the parameters in model H1. Bottom-right panel: Best-fit values for the straight-lines for Case 2 inferred with the datasets shown.

Figure 13. BAO Hubble diagram. BAO measurements of DV/rd, DM/rd, and zDH/rd from the sources indicated in the legend. The scaling factor z is included for a better display of the error bars. Solid lines are plotted by using the best-fit values obtained by the Planck satellite [85]. The Lyα cross-correlation points have been shifted in redshift; auto-correlation points are plotted at the correct effective redshift.

Figure 14. Joint Light-curve Analysis data (JLA). Vertical axis is the standardized distance modulus μ (luminosity distance function), and the horizontal axis is the redshift z. Source: Reference [93].

Figure 15. Hubble data as a function of redshift z for Cosmic Chronometers. The solid line corresponds to the best-fit by using ΛCDM model.

Figure 16. Results for the ΛCDM model. Left panel: One-dimensional posterior distributions for the parameters H0 and Ωm, along with its chains. Right panel: Joint 2D posterior distributions with 1–4 confidence levels.

Figure 17. Top panel: Markov chain for the parameter H0, with 9000 steps. Bottom panel: Two-dimensional posterior distribution with confidence regions 1 and 2 σ for the joint parameters H0 and Ωm.

Figure 17. Top panel: Markov chain for the parameter H0, with 9000 steps. Bottom panel: Two-dimensional posterior distribution with confidence regions 1 and 2 σ for the joint parameters H0 and Ωm.

Figure 18. Two-dimensional posterior distributions with 1,2-σ confidence regions for different models. oΛCDM refers to ΛCDM with curvature, ωCDM is a flat Universe with the dark energy equation of state as a free parameter, oωCDM generalizes to non-zero curvature, wowaCDM uses the CPL parameterization, and owowaCDM generalizes to non-zero curvature. The dashed lines show the standard ΛCDM values.

Table 1

Main differences between the Bayesian and Frequentist interpretations.

Frequentist	Bayesian
Data are a repeatable random	Data are observed from
sample. There is a frequency.	the realized sample.
Underlying parameters remain	Parameters are unknown and
constant during this repeatable	described probabilistically.
process.
Parameters are fixed.	Data are fixed.

Table 2

Jeffreys guideline scale for evaluating the strength of evidence when two models are compared.

$|{\mathcal{B}}_{0,1}|$	Odds	Probability	Strength
<1.0	<3:1	<0.750	Inconclusive
1.0–2.5	∼12:1	0.923	Significant
2.5–5.0	∼150:1	0.993	Strong
>5.0	>150:1	>0.993	Decisive

Table 3

$\mathsf{\Delta }{\chi }^2$ for the conventional $68.3\%$, $95.4\%$, and $99.73\%$ as a function of the number of parameters (M) for the joint confidence level.

			$\mathsf{\Delta }{\mathit{\chi }}^2$
$\mathit{\sigma }$	$\mathit{p}$	$\mathit{M}=\mathbf{1}$	$\mathit{M}=\mathbf{2}$	$\mathit{M}=\mathbf{3}$
1	$68.3\%$	$1.00$	$2.30$	$3.53$
2	$95.4\%$	$4.00$	$6.17$	$8.02$
3	$99.73\%$	$9.00$	$11.8$	$14.20$

Table 4

Equation of state associated with each component of the Universe.

Component	$\mathit{\omega }$
Dust	0
Radiation	1/3
Cosmological Constant	−1

Table 5

Cosmological parameter constraints from BAO data combined with our compressed description of CMB from Planck, the JLA SN, and Hubble data (BBAO + Planck + SN + HD). Two-tailed distributions are shown, along with 1$\sigma$ C.L. Entries for which the parameter is fixed are marked with dash (${\mathsf{\Omega }}_k=0$, $\omega ={\omega }_0=-1$, ${\omega }_a=0$).

Parameter	$\mathit{\Lambda }$CDM	o$\mathit{\Lambda }$CDM	$\mathit{\omega }$CDM	o$\mathit{\omega }$CDM	${\mathit{\omega }}_0{\mathit{\omega }}_{\mathit{a}}$CDM	o${\mathit{\omega }}_0{\mathit{\omega }}_{\mathit{a}}$CDM
${\mathsf{\Omega }}_m$	$0.299\pm 0.007$	$0.298\pm 0.007$	$0.303\pm 0.009$	$0.299\pm 0.009$	$0.307\pm 0.010$	$0.306\pm 0.010$
${\mathsf{\Omega }}_b{h}^2$	$0.0224\pm 0.0002$	$0.0227\pm 0.0003$	$0.0224\pm 0.0003$	$0.0227\pm 0.0003$	$0.0224\pm 0.0003$	$0.0226\pm 0.0003$
h	$0.684\pm 0.006$	$0.679\pm 0.007$	$0.677\pm 0.108$	$0.676\pm 0.010$	$0.674\pm 0.011$	$0.670\pm 0.011$
${\mathsf{\Omega }}_k$	⋯	$-0.004\pm 0.002$	⋯	$-0.003\pm 0.003$	⋯	$-0.006\pm 0.003$
${\omega }_0$	⋯	⋯	$-0.96\pm 0.05$	$-0.97\pm 0.05$	$-0.91\pm 0.10$	$-0.83\pm 0.11$
${\omega }_a$	⋯	⋯	⋯	⋯	$-0.17\pm 0.41$	$-0.52\pm 0.51$
${\chi }_{\min }^2$	$73.57$	$71.59$	$73.13$	$71.5$	$72.8$	$69.7$
$|{\mathcal{B}}_{\mathsf{\Lambda }\mathrm{CDM},i}|$	0	3.61	1.48	5.59	1.23	4.3