Cosmological parameter inference with Bayesian statistics
Luis E. Padilla Luis O. Tellez Luis A. Escamilla J. Alberto Vazquez





  

Departamento de Física, Centro de Investigación y de Estudios Avanzados del IPN, A.P. 14-740, 07000 México D.F.,
  México.
  Department of Astronomy and Texas Cosmology Center, University of Texas, Austin, TX, 78712-1083, U.S.A.
  
Departamento de Física, Centro de Investigación y de Estudios Avanzados del IPN, A.P. 14-740, 07000 México D.F.,
  México.
  Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México, 
Apdo. Postal 48-3, 62251 Cuernavaca, Morelos, México.
      
Departamento de Física, Centro de Investigación y de Estudios Avanzados del IPN, A.P. 14-740, 07000 México D.F.,
  México.
  Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México, 
Apdo. Postal 48-3, 62251 Cuernavaca, Morelos, México.
   
javazquez@icf.unam.mx
Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México, 
Apdo. Postal 48-3, 62251 Cuernavaca, Morelos, México.
Departamento de Física, Centro de Investigación y de Estudios Avanzados del IPN, A.P. 14-740, 07000 México D.F.,
  México.





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 Λ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.




Cosmological parameter inference with Bayesian statistics
    J. Alberto Vazquez
    December 30, 2023
=========================================================




§ INTRODUCTION




The beginning of the standard cosmology as it is known today emerged after 1920, when 
the Shapley-Curtis debate was carried out <cit.>. 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 <cit.>. 

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 <ref>, we point out the
main differences between the Bayesian and Frequentist approaches. Then, in Section 
<ref>, 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 <ref>, 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 <ref>.
In Section <ref>, we present an introduction to cosmology, 
and then, in Section <ref>, 
focus on some codes to compute the cosmological observables. 
In Section <ref>, we constrain the 
parameter space that describes 
the standard cosmological model, namely the ΛCDM model, along with several extensions. 
Finally, in Section <ref>, we present our conclusions.




§ 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

    P=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 <cit.> and references therein. 

Let x be a random variable related to a particular event and P(x) its corresponding
probability distribution, for both cases, the same rules of probabilities 
apply [These rules are defined for a continuous variable; however, the 
corresponding discrete definition can be given immediately by replacing 
∫ dx →∑.]:




    P(x)≥ 0,


    ∫_-∞^∞ dxP(x)=1.



For mutually exclusive events, we have





    P(x_1∪ x_2)=P(x_1)+P(x_2),


but, in general



    P(x_1∪ x_2)=P(x_1)+P(x_2) - P(x_1 ∩ x_2).



These rules are summed up as follows:
the first condition (<ref>) is necessary due to the probability of having an
event is always positive; the second rule (<ref>) is a normalized relation, 
which tells us that we are certain to obtain one of the possible outcomes;
now, in the third point (<ref>) 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 

    P(x_1∩ x_2)=P(x_1)P(x_2|x_1).




If two events x_1 and x_2 are mutually exclusive, then 

    P(x_1∩ x_2) = 0 = P(x_2∩ 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.




 §.§ 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,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:B̂→ [0,1] in the sample space (x,B̂) representing the “true 
distribution of the data” 
    X∼ P_0.


For Frequentist statistics, the model Q is a collection of
probability measurements P_θ:B̂→[0,1] in the sample space (x,B̂). 
The distributions P_θ are called model distributions, with θ being the model parameters; in this approach θ 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∈ Q.


Finally, we need a point-estimator (or estimator) for P_0. An estimator 
for P_0 is a map P̂:x→ Q, representing our best guess 
P̂∈ 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?”





 §.§ Bayesian statistics




In Bayesian statistics, data and model are two elements 
of the same space <cit.>, 
i.e., no formal distinction is made between measured quantities X and 
parameters θ. One may envisage the process of generating a measurement's outcome Y=y as 
two draws, one draw for Θ (where Θ is a model with associated probabilities 
to the parameter θ) to select a value of θ and a subsequent 
draw for P_θ 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_θ of the model is interpreted 
simply as the distribution of X given the parameter value θ, i.e., as the conditional
distribution X|θ.





 §.§ Comparing both descriptions






Table <ref> 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.- 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(head)=p=#  of  heads/#  of  events=10/14. 
So, the probability of having 2 heads in a row is P(2heads)=P(head)P(head)≃ 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


    P(2heads|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(2heads)=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 <cit.>. In Section <ref>, 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.

 



§ 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.   




 §.§ 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 (<ref>) - (<ref>). From Equation (<ref>), 
without loss of generality, it must be fulfilled that P(x_1∩ x_2)=P(x_2∩ x_1). 
In such a case, the following relation applies:


    P(x_2|x_1)=P(x_1|x_2)P(x_2)/P(x_1).



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→ D as a set of data, 
and x_2→θ as the parameter vector of said hypothesis, we can rewrite 
the above equation as


    P(θ|D, H)=P(D|θ,H)P(θ|H)/P(D|H).



This last relation is the so-called Bayes' theorem and the most important tool 
in a Bayesian inference procedure. In this result, P(θ|D, H) is called the
posterior probability of the model. L(D|θ,H) ≡ P(D|θ,H) 
is called the likelihood, and it will be our main focus in future sections,
π(θ) ≡ P(θ|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), 𝒵≡ P(D|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):


    π(θ)∝ 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 <ref>. 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(θ;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

    B(θ;a,b)=Γ(a+b)/Γ(a)Γ(b)θ^a-1(1-θ)^b-1,

where Γ 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 

    P(D|H)=∫ d^Nθ P(D|θ,H)P(θ|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



    K≡P(D| H_0)/P(D| H_1)= 
    ∫ d^N_0θ_0 P(D|θ_0, H_0) P(θ_0| H_0)/∫ d^N_1θ_1 P(D|θ_1, H_1) P(θ_1| H_1)=
         𝒵_0/𝒵_1,


or equivalently the difference in log evidence ln𝒵_0 - ln𝒵_1
is often termed as the Bayes factor ℬ_0,1:



    ℬ_0,1=ln𝒵_0/𝒵_1,

where θ_i is a parameter vector (with dimensionality N_i) 
for the hypothesis H_i and i=0,1. 

In Equation (<ref>), the quantity ℬ_0,1=ln K 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 <ref> <cit.>).






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.- 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(2heads|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(2heads|p). 
This parameter p will have a probability distribution P(p|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

    P(2heads|D)=∫^1_0 P(2heads|p)P(p|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

    P(2heads|p)=[P(head|p)]^2 ,


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

    P(head|p)=p   ⇒   P(2heads|p)=p^2.



On the other hand, notice that we do not know a priori the quantity 
P(p|D) but P(D|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

    P(x|p,n)=nxp^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(D|p). 
 Using the Bayes' formula, we have

    P(p|D)=P(D|p)P(p)/P(D).



For the prior, we will use the beta distribution (<ref>), so


    P(p)=B(p;a,b).


In order to get the explicit form of P(p|D), we still need to compute P(D). 
That is, by plugging Equations (<ref>) and (<ref>) into the integral of Equation (<ref>) yields to

    P(D)=B(10+a,4+b)≡Γ(10+a)Γ(4+b)/Γ((10+a)+(4+b));


therefore, 

    P(p|D)=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 <ref> that our posterior result (Red figure) described by 
Equation (<ref>) 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 (<ref>) using (<ref>) and (<ref>), 
we arrive at the result obtained in the previous section

    P(2heads|D)=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). 

 



 






 §.§ 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 <ref>). 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(2heads|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 <ref>a, and 
then, by considering the 500 data points, we can arrive at the posterior in 
Figure <ref>d, or similarly 
start with the posterior shown in Figure <ref>c as our prior and consider only the 
last 400 data points to obtain the same posterior, displayed in Figure <ref>d. 


 



In fact, if we rewrite Bayes' theorem so that all probabilities are explicitly dependent on 
some prior information I <cit.>

    P(θ|DI,H)=P(θ|I,H)P(DI|θ,H)/P(D|I,H),


and then consider a new set of data D', letting the old data become part of the prior 
information I'=DI, we arrive at 

    P(θ|D'I',H)    =   P(θ|I,H)P(DD'I|θ,H)/P(DD'|I,H)
       =   P(θ|[DD']I, H),


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




 §.§ 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(D|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(θ|D, H)∝ L(D|θ,H); thus, by maximizing it, we can find the 
most probable set of parameters for a model given the data. However, 
having ignored P(D|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

    Λ=-2ln[L(D|θ, H)/L_ max],

is bigger than a given value.

Let us assume we have a single-peaked distribution. 
We consider that θ̂ is the mean of the distribution 

    θ̂=∫ dθθ P(θ|D, H).



If our model is well-specified and the expectation value of θ̂ corresponds to the real or most probable value θ_0, we have 

    ⟨θ̂⟩=θ_0,


then we say that θ̂ is unbiased. Considering a Taylor expansion of the 
log likelihood around its maximum, we have 


    ln L(D|θ)   =   ln L(D|θ_0) 
       +   1/2(θ_i-θ_0i)∂^2ln L/∂θ_i ∂θ_j(θ_j-θ_0j) 
       +   ...,

where θ_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 

    L(D|θ)=L(D|θ_0)exp[-1/2(θ_i-θ_0i)H_ij(θ_j-θ_0j)],

where 

    H_ij=-∂^2ln L/∂θ_i ∂θ_j,

is called the Hessian matrix. It controls whether the estimates of 
θ_i and θ_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 <cit.>.




 §.§ 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 <cit.>. In this example, there are two people, 
A and B, that are interested in the measurement of a given physical quantity θ. 
A and B have different prior beliefs regarding the possible value of θ. 
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(θ|I_i), (i=A,B), and assume they are described by two Gaussian distributions 
with mean μ_i and variance Σ_i^2. Now, A and B measure 
θ together using an apparatus subject to a Gaussian noise with known variance σ. 
They obtain the value θ_0=m_1. Therefore, they can write their likelihoods 
for θ as 

    L(D|θ, HI)=L_0exp[-1/2(θ-m_1)^2/σ^2].



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

    P(θ|m_1)=L(m_1|θ I_i)P(θ|I_i)/P(m_1|I_i),

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

    μ̂_i = m_1+(σ/Σ_i)^2μ_i/1+(σ/Σ_i)^2,


and variance 

    τ_i^2=σ^2/1+(σ/Σ_i)^2,   (i=A,B).



Thus, if the likelihood is more informative than the prior, i.e., (σ/Σ_i)≪ 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 σ^2 by σ^2/N. 
Then, we can see that the initial prior μ_i of A and B will progressively be overridden 
by the data. This process is illustrated in Figure <ref> where the green (red) 
curve corresponds to the probability distribution of θ for person A (B) and the 
blue curve corresponds to their likelihood.



 





 §.§ 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 (<ref>), we can see the 
likelihood will be maximum if the quantity

    χ^2≡(θ_i-θ_0i)H_ij(θ_j-θ_0j),


is minimum. The quantity χ^2 is usually called chi-square and is related 
to the Gaussian likelihood via L=L_0e^-χ^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 χ^2 around its minimum is given 
by the χ^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 χ^2 exceeds, by chance, a value χ̂ for the 
correct model.
This probability is given by Q(v,χ̂)=1-Γ(v/2,χ̂/2) <cit.>, where 
Γ is the incomplete Gamma function. Then, the probability that the observed χ^2 
(even the correct model) is less than a given value χ̂^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.





 §.§ 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 χ^2 boundaries. 
When the χ^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 (<ref>)) these are ellipsoidal regions. 
In this section, we exemplify how to calculate the confidence regions, 
following Reference <cit.>. 

We consider a small perturbation from the best fit of chi-square 
Δχ^2=χ^2-χ^2_ best. Then, we use the properties of χ^2 
distribution to define confidence regions for variations on χ^2 to its minimum. 
In Table <ref>, 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 σ confidence levels. As an example, we plot in Figure <ref> the corresponding confidence regions associated with the coin example.






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



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

  * Solve the equation:

    Q(n-M,min(χ^2)+Δχ^2)=p.


  * Find the parameter region where χ^2≤ min(χ^2)+Δχ^2. This defines the confidence region.





 §.§ Marginalization



It is clear that a model may (in general) depend on more than one parameter. 
However, some of these parameters θ_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

    P(θ_1,...,θ_j,H|D)=∫ dθ_j+1...dθ_mP(θ,H|D),

where m is the total number of parameters in our model, and θ_1,...,θ_j 
denote the parameters we are interested in.




 §.§ 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 <cit.>. 
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/√(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

    F_ij=-⟨∂^2 ℒ/∂θ_i ∂θ_j⟩,

where 

    ℒ=ln L.



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

As we can see from Equation (<ref>), 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 θ_i with a 
Gaussian likelihood. In this scenario, 



    Δℒ=1/2F_ii(θ_i- θ_0i)^2,

when 2Δℒ=1, and identifying the Δχ^2 corresponding 
to 68% confidence level, we notice that 1/√(F_ii) yields the 1-σ 
displacement for θ_i. In the general case,



    σ_ij^2 ≥ (F^-1)_ij.



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

    σ_θ_i≥ (F^-1)^1/2_ii.



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 (<ref>) 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 <ref>, 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 θ_1,θ_2,...,θ_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(θ_1,θ_2,...,θ_N). Then, the elements of the 
Fisher matrix can be computed as

    F_ij=∑_k 1/σ_k^2∂ f_k/∂θ_i∂ f_k/∂θ_j,

where σ_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) <cit.>, 
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).






 §.§ 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 <cit.>.
 

Suppose we are interested in computing the expectation value μ_f=E_p[f(X)], 
where f(X) 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(x) that satisfies q(x)>0 whenever f(x)p(x)≠ 0, we can 
rewrite the mean value μ_f as

    μ_f = ∫ f(x)p(x)dx    =   ∫ f(x)p(x)/q(x)q(x)dx 
       =   E_q[f(X)w(x)],

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



    μ̂_f = 1/m∑_j=1^m w(x^(j))f(x^(j)).



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

    μ̂_f=∑_j=1^m w(x^(j))f(x^(j))/∑_j=1^mw(x^(j)).



For the strong law of large numbers, in the limit when m→∞, 
we will have that μ̂_f→μ_f.

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

    P'(D)/P(D)=E[P'(θ,D)/P(θ,D)]_P(θ|D)≃1/N∑_n=1^NP'(D|θ_n)P'(θ_n)/P(D|θ_n)P(θ_n),

where the samples {θ_n} are drawn from P(θ|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.





 §.§ Combining datasets: Hyperparameter method




Suppose we are dealing with multiple datasets { D_1,…,D_N}, coming from a 
collection of different surveys { S_1,…,S_N}. 
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 <ref>, 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 <cit.> 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 
α 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.,

    P(θ|D,H)=1/P(D|H)∫ P(θ|α,H)P(α|D, H)dα,

where we have used the Bayes' theorem. Now, for the method, it is necessary to assume that the 
hyperparameters α and the parameters of interest θ are independent, i.e.,
P(θ,α,H)=P(α)P(θ,H), and it is also necessary to assume that each 
hyperparameter α_k is independent of each other, i.e., 
P(α)=P(α_1)P(α_2)… P(α_N). 
In this way, we can rewrite the above expression as

    P(θ|D,H)=P(θ,H)/P(D|H)[∏_k=1^N∫ P(D_k|θ,α_k,H)P(α_k)dα_k].



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(D_k|θ,α_k,H), for every 
α_k and the dataset D_k; P(D|H) is the evidence and, similarly to a parameter 
inference procedure, it works as a normalizing function, i.e., 
P(D|H)=∫ d θ P(θ,H)L(D|θ,H). 
Notice that, by considering P(α_k)=δ(α_k-1), we rely on the standard 
approach, where no hyperparameters are used. 

We add these α_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 (<ref>).

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(D|θ,H_1) for model H_1 as

    P(D|θ,H_1)=∏_k=1^N2Γ(n_k/2+1)/π^n_k/2|V_k|^1/2(χ_k^2+2)^-(n_k/2+1),


obtaining an explicit functional form for K, given by 

    K=∏_k=1^N2^n_k/2+1Γ(n_k/2+1)/χ^2_k+2e^-χ_k^2/2.



Here, χ_k^2 is given by (<ref>) for every dataset, and n_k is the number of 
points contained in D_k. In Equation (<ref>), 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(D|H_i) is the
key quantity for making a comparison between two different 
models. In fact, by using the Bayes factor K from Equation (<ref>), we can estimate the 
necessity to introduce the hyperparameters to our model using the criteria given in 
Table <ref>.
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 (<ref>). 






§ 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 . 




 §.§ 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 (<ref>). In this 
section, we review the basic steps for this procedure in a simplistic way; however, 
it is recommendable to check <cit.> 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,… of elements of some set is a Markov Chain 
if the conditional distribution of X_n+1 given X_1,…,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(θ|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:

    p(θ|D,H)≃1/N∑_i=1^N δ(θ-θ_i),


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

    ⟨θ⟩=∫ dθθ P(θ,H|D)≃1/N∑_i=1^Nθ_i ,


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

    ⟨ f(θ)⟩≃1/N∑_i=1^N f(θ_i).



As mentioned before, in a Markov Chain, it is necessary to generate a new point θ_i+1 
 from the present point θ_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. 





  §.§.§ Metropolis-Hastings algorithm




In the Metropolis-Hastings algorithm <cit.>, it is necessary to 
start from a random initial point θ_i, with an associated posterior 
probability p_i=p(θ_i|D,H). We need to propose a candidate θ_c 
by drawing from a proposal distribution q(θ_i,θ_c) 
 used as a generator of new random steps. Then, the probability of 
acceptance of the new point is given by

    p(acceptance)=min[1,p_cq(θ_c,θ_i)/p_iq(θ_i,θ_c)].



If the proposal distribution is symmetric, the algorithm is reduced to the 
Metropolis algorithm
  
    p(acceptance)=min[1,p_c/p_i].



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


  * Choose a random initial condition θ_i in the parameter space and compute 
the posterior distribution.

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

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

  * If the point is rejected, repeat the previous point in the chain.

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






  §.§.§ 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 <ref>. 
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 (<ref>), we use the MHA. 

As mentioned before, we consider a likelihood given by the binomial distribution 
(<ref>) and a normal distributed prior (<ref>) (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,σ̂), where 
G(p_cu,σ̂) is our proposed Gaussian distribution centered 
at p_cu with variance σ̂=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 <ref>) is a posterior distribution 
that matches very well with the results calculated analytically (shown in Figure <ref>). 
Numerically, we obtained p= 0.695^+0.123_-0.107, where the upper and lower 
values for p correspond to the 1σ standard deviation.
Notice that we have plotted the width of our 1σ, 2σ 
and 3σ confidence regions in the same figure. 
To complement the example (and Figure <ref>), we also show in Figure <ref> 
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 <ref>, 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 <cit.> 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.





  §.§.§ 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 <ref>, 
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 
<cit.>. 
That is, by starting with M 
chains with very different initial points and N points per chain, 
if θ_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:

    ⟨θ^j⟩ =1/N∑_i=1^N θ_i^j,

and the mean of all the chains

    ⟨θ⟩ =1/NM∑_i=1^N∑_j=1^Mθ_i^j.


Then, the chain-to-chain variance B is

    B=1/M-1∑_j=1^M(⟨θ^j⟩-⟨θ⟩)^2 ,

and the average variance of each chain is

    W=1/M(N-1)∑_i=1^N∑_j=1^M(θ_i^j-⟨θ^j⟩)^2 .



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

    R̂=N-1/NW+B(1+1/M)/W,


which is the ratio of the two estimates, approaches unity. A typical convergence 
criterion is when 0.97<R̂<1.03. 

 





  §.§.§ 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) σ̂=0.1. 
This value was taken because we initially explored different values for 
σ̂, 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 <ref> (top panel), 
considering an initial step p_i=0.6 and a variance for the proposal 
distribution σ̂= 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 σ̂= 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 <cit.>. 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 lag k is defined as the 
correlation between every sample and the sample k steps before. 
It can be quantified as <cit.>

    ρ_k    =   Cov(X_t,X_t+k)/√(Var(X_t)Var(X_t+k))
       =   E[(X_t-X)(X_t+k-X)]/√(E[(X_t-X)^2]E[(X_t+k-X)^2]),

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
<cit.> is to split the multidimensional θ 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:


  * θ consists of k blocks θ_1, …, θ_k. Then, at step i

  * Draw θ^i+1_1 from p(θ_1|θ^i_2, …,θ^i_k)

  * Draw θ^i+1_2 from p(θ_2|θ^i+1_1, θ^i_3, …,θ^i_k)

  * …
  * Draw θ^i+1_k from p(θ_k|θ^i+1_1, θ^i+1_2, …,θ^i+1_k-1)

  * Repeat the above steps for the wished iterations with i → i+1.


The distribution p(θ_1|θ_2, …,θ_k)=p(θ_1, …, θ_k)/p(θ_2, …, θ_k) is known as the 
full conditional distribution of θ_1. This algorithm 
is a special case of MHA where the proposal is always accepted.




Metropolis Coupled Markov Chain Monte Carlo (MC^3): It is easy to see that it could be a little problematic if our likelihood has local maxima. The MC^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(θ,T|D,H):

    P(θ,T|D,H) ∝ L(θ,D)^1/TP(θ,H),


where L is the likelihood, and P(θ,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<…<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

    A_i,j=min{(L(θ_i)/L(θ_j))^1/T_j-1/T_i,1} .


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 θ_j and chain j has to move to 
position θ_i. 

We can see that, since the hottest chain T_ max can access easier to all the 
modes of P(θ,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 <cit.>.


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: 

    p(x_1, x_2) ∝exp( -(x_1-x_2)^2/2ϵ-(x_1+x_2)^2/2),

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

    y_1=x_1-x_2/√(ϵ),    y_2=x_1+x_2,

we can rewrite the anisotropic density into the easier problem 

    p(y_1, y_2) ∝exp( -(y_1^2+y_2^2)/2).



An MCMC sampler has the form X(t+1)=R(X(t),ψ (t),p), where X(t) is the 
sample after t iterations, R is the sampler algorithm, ψ 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,

    R(AX(t)+b,ψ (t),p_A,b)=AR(X(t),ψ (t),p)+b.


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


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 <cit.>) or the Adaptative Metropolis-Hastings 
(AMH) (see, e.g., Reference <cit.>). 



 
 

§ 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 <ref>):





    
  * Consider two datasets coming from the same 
straight-line but having different errors.
    
  * 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 <cit.>. 








 §.§ 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 <ref>. 
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 σ_1 = 0.3, while, for D2 
we use σ_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

    a ∝ U[0,5]    and   b ∝ U[0,3],

where U[α,β] are uniform distributions with lower limit 
α and upper limit β.

From Equation (<ref>), we can write our likelihood as

    L(D;line)∝exp[-∑_d (y_d-y)^2/2σ_d^2],


where y_d is our data taken from the dataset D=D_1+D_2, and σ_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 ± 0.047 and b=1.994 ± 0.013,
and their posterior distributions are plotted in
Figure <ref>. 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 <ref> displays the 1-4 σ 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 <ref> 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 (<ref>). 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 ± 0.038
with Gelman-Rubin of 1.000113, and b=1.995± 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 <ref> and the top panel of Figure <ref> 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 (<ref>) 

    K = 3.


Then, comparing with Table <ref>, 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 <ref> 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.  






 §.§ 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 <ref>. 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 σ = 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 ± 0.056 and b = 1.795± 0.014.
Then, we plotted 1-4 σ confidence regions in the upper panel of Figure <ref>. 
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 <ref>, 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

    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 <ref>, we obtain 
the bottom-right panel of the same Figure.









§ BAYESIAN STATISTICS IN COSMOLOGY






 §.§ 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. 




  §.§.§ 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_μν. 
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

    G_μν = 8π G T_μν,

where G_μν is the Einstein tensor (geometry of the spacetime),
and G is the gravitational Newton constant <cit.>.
Throughout this review, we use natural units 
 c=ħ=1.
 
The distance between two points in a curved 
space-time can be measured as

    ds^2 = g_μνdx^μdx^ν,

where g_μν 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 μ, ν, …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

    ds^2 = dt^2 -a^2(t)γ_ijdx^idx^j ,

where

    γ_ij≡δ_ij + κx_i x_j/1 - κ(x_k x^k).


In Equation (<ref>), a represents the scale factor of the Universe 
which only depends on time, and by convention is normalized to today 
a(t_0) ≡ 1.
Similarly, in expression (<ref>), x^i labeled the spatial coordinates 
(also called comoving coordinates), δ_ij is the Kronecker delta 
and κ describes the curvature of the space-time. 





  §.§.§ 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 

    T_μν = (ρ + P)U_μU_ν - Pg_μν,
 
where ρ is the energy density, P is the fluid pressure, and U_μ is the 
4-velocity relative to the observer. If we take the velocity as 
U^μ = (1,0,0,0) (comoving observer), the energy-momentum tensor reduces to

    T^μ_ν = g^μλ T_λν = [  ρ  0  0  0;  0 -P  0  0;  0  0 -P  0;  0  0  0 -P;    ] .


Using Equations (<ref>) and (<ref>), with the FLRW metric,
we can obtain the Friedmann equations


    H^2  ≡(ȧ/a)^2 = 8 π G/3∑_iρ_i - κ/a^2,


    ä/a = -4 π G/3∑_i(ρ_i + 3P_i).



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 
(<ref>) and (<ref>), we can obtain the continuity equation given by

    ρ̇ + 3 ȧ/a(ρ + P) = 0.


The meaning of (<ref>) is the conservation of the energy-momentum tensor
(∇_μ T^μν=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=ωρ.





  §.§.§ 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 Λ 
Cold Dark Matter (Λ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 ρ∝ 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 <cit.>. The Λ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 <cit.>.
    
	
	
  * Radiation: This corresponds to relativistic particles that follow the relation 
	P = 1/3ρ, which implies a density with a behavior ρ∝ a^-4. 
	We consider photons ρ_γ and massless neutrinos ρ_ν
	as radiation, so the total radiation 
	energy density in the Universe is given by
	
    ρ_r = ρ_γ + ρ_ν.

	
	The relation between these quantities is
	
    ρ_ν = N_ eff×7/8×(4/11)^4/3ρ_γ,

	where N_ eff is the effective number of relativistic degrees of freedom, 
	with standard value N_ eff = 3.046 <cit.>.
	
	
  * Dark Energy: It is introduced to explain the current accelerated
	expansion of the Universe. In the ΛCDM model, dark energy is given 
	by the cosmological constant Λ or equivalently 
	by an equation-of-state ω=-1.








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

    Ω_i≡ρ_i/ρ_ crit,    with  ρ_ crit = 3H^2/8 π G,
 

with ρ_ crit being the condition to have a 
flat Universe or equivalently zero curvature, we can rewrite (<ref>) as

    H^2/H_0^2 = Ω_r,0 a^-4 + Ω_m,0 a^-3 + Ω_k,0 a^-2 + Ω_Λ,0,

where Ω_r,0 is the radiation density parameter, 
Ω_m,0≡Ω_b,0+Ω_DM,0 corresponds to the total matter, 
Ω_b,0 to baryons, Ω_DM,0 to dark matter, Ω_k ≡ -κ/(a H)^2 the curvature 
density parameter, 
and Ω_Λ≡Λ /3H^2 associated with the Cosmological Constant, and the 
subscript zero indicates they are evaluated by 
today's values. 






  §.§.§ Alternatives to the ΛCDM model
 



The Λ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
<cit.>). 
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 Λ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 <cit.>;
or the 
Self Interacting Dark Matter, as its name states, relies on the cold dark 
matter to be made of self interacting particles <cit.>.
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(R) theories <cit.>, 
braneworld models <cit.>.
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
<cit.>); or many more alternatives, i.e.,
anisotropic Universes <cit.>.
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 ω in a Taylor series, 
for example, the Chevallier-Polarski-Linder (CPL) 
ω = ω_0 + ω_a(1-a), with
two free parameters ω_0, ω_a <cit.>. It may also be 
expanded into Fourier series <cit.>, or many more Bayesian approaches, as
have been suggested to account for a dynamical dark energy <cit.>.






 §.§ Cosmological parameters






  §.§.§ 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 Ω_i,0 and the Hubble parameter H_0 related 
by (<ref>). 

In particular, the radiation contribution is measured with great precision, so that
Ω_γ 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 (<ref>). 

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





  §.§.§ 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 (ΛCDM), in order to decrease degeneracies, the 
physical energy densities Ω_DM,0h^2 and Ω_b,0h^2 are used as 
base parameters <cit.>.




 §.§ 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

    r_d = ∫_z_d^∞c_s(z)/H(z)dz,

where the sound speed (in terms of redshift z) in the photon-baryon fluid is c_s(z) = 3^-1/2c[1+3/4ρ_b(z)/ρ_γ(z)]^-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 (α) of the measured BAO scale and that predicted by the fiducial model (fid). In an anisotropic fit, two ratios are used, one perpendicular α_ and one parallel α_∥ to the line of sight. A measurement of α_ constrains the ratio of the comoving angular diameter distance to the sound horizon <cit.>:

    D_M(z)/r_d = α_D_M,fid(z)/r_d,fid,


where the comoving angular diameter distance is given by

    D_M(z) = c/H_0S_k(D_c(z)/c/H_0).


The line-of-sight comoving distance is defined as


    D_c(z) = c/H_0∫_0^z dz'H_0/H(z'),


and S_k (z)


    S_k(x) = 
    sinh(√(Ω_k)x)/√(Ω_k)   Ω_k > 0,
    
    x    Ω_k=0,
    sin(√(-Ω_k)x)/√(-Ω_k)   Ω_k < 0.


The Hubble parameter can be constrained by measuring α_∥ using an analogous quantity

    D_H(z)/r_d = α_∥D_H,fid(z)/r_d,fid,

with D_H(z) = c/H(z).


If redshift-space distortions are weak (which is valid for luminous galaxy surveys but not for the Ly-α), an isotropic analysis measures an effective 
combination of (<ref>) and (<ref>), 
and the volume averaged distance D_V(z) <cit.>

    D_V(z)/r_d = αD_V,fid(z)/r_d,fid,


with D_V(z) = [zD_H(z)D_M^2(z)]^1/3.

The BAO measurements constrain the cosmological parameters through the radius of 
the sound horizon r_d, Hubble distance D_H(z) and the comoving angular diameter 
distance D_M(z); see Figure <ref>.
The data used in this work to constrain D_V / r_d is obtained from the 6dF 
Galaxy Survey (6dFGS <cit.>) from UK Schmidt Telescope, the Main Galaxy Sample 
(MGS <cit.>), and the 
BOSS LOWZ Sample <cit.> from SDSS. On the other hand, the data from 
BOSS CMASS Sample <cit.>, BOSS 
Lyman α auto-correlation (Lyα auto, <cit.>), and BOSS Lyman-α cross 
correlation (Lyα cross <cit.>) 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 <cit.>. 
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 <cit.>. 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 <cit.>.
 
These stars allow us to measure relative distances  
using the luminosity distance given by

    D_L ≡√(L/4π 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 <cit.>. 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 <cit.>. The relation between D_L and the
cosmological parameters is given by

    D_L = D_M (1+z),


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

    μ = m_B^* - M_B + α X_1 -β C,


where m_B^* is the observed peak magnitude in the rest frame of blue band (B), 
α, β, and M_B are parameters that depend on 
host galaxy properties <cit.>. 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 

    μ = 5 log_10(D_L/10pc).


The data used in this paper (shown in Figure <ref>) 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 <cit.>) 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) <cit.>. 

For simplicity, we compress all the information into a linear 
function fit over 30 bins (31 nodes), spaced evenly 
in log(z) with a 31× 31 covariance matrix <cit.>.

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 <cit.>.






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 <cit.> and WMAP <cit.>.
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 <cit.>.

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(1090)/r_d, in addition, to calibrating the absolute length of the BAO ruler through the 
determination of Ω_b h^2 and Ω_cb h^2 (the density parameters of baryonic and dark matter, respectively; more details about these parameters can be found in Reference <cit.>).


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 ∼ 0-2 and compare their upper cut in its age distributions. By doing this, it is possible to obtain the difference of ages Δ t and redshifts Δ z such that the expression dz/dt can be approximated as dz/dt≃Δ z/Δ t. This quantity is related with the Hubble factor via <cit.>

    H(z) = -1/1+zdz/dt≃-1/1+zΔ z/Δ t.


In this work, the data used to constrain H(z) was obtained from Reference <cit.>. 
In Figure <ref>, 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 <cit.>.
 
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 <cit.>, clustering <cit.>, the matter power spectrum <cit.>, and primordial gravitational waves and the primordial power spectrum, both generated during inflation <cit.>, among many others. 



 


§ PARAMETERS INFERENCE IN COSMOLOGY














Here, we estimate the Hubble parameter H_0 at the present time and the matter density 
 of the Universe Ω_m, assuming a ΛCDM standard model. 
We use our own Python code, which can be found in 
Reference <cit.>. 


The data.-
For this particular example, we focus only on the Cosmic Chronometers, 
as shown in Figure <ref>.




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 (<ref>), 
and, for simplicity, we assume a flat Universe with a well measured radiation content.
Then, the model is given by


    H^2 = H_0^2[Ω_m,0 (1+z)^3+ Ω_Λ,0],


constrained by the relation

    Ω_m,0+Ω_Λ,0=1.


Notice that the above relation implies that we can get rid of a parameter—Ω_m,0 or Ω_Λ,0—in the analysis. Then, the parameters 
we decided to estimate are H_0 and Ω_m,0. By assuming the Gaussian 
approximation, we construct the likelihood given by

    L ∝exp[-∑_i(H_i-H(z_i))^2/2σ_i^2],


where H(z_i) is described by Equation (<ref>) evaluated at each redshift z_i, and
H_i is the value of the Hubble parameter measured at z_i and σ_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 Ω_i,0 must satisfy the 
relation 0≤Ω_i,0≤ 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 Ω_m∈ [0,1] and H_0∈ [10,100].

Hence, the priors

    Ω_m,0∼ U[0,1],      H_0∼ U[10,100].


Numerical estimation.- We follow the same procedure as in the straight-line
example. In the left panel of Figure <ref>, 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σ 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 ± 3.13 with convergence 1.00045, and 
Ω_m = 0.331 ± 0.0628 with convergence 1.00044. 
We can see that our estimations are very similar to
the values reported in the literature <cit.>. 






 §.§ 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 (<ref>) and (<ref>). 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.





  §.§.§ Cosmological codes



Today, there are several cosmological Boltzmann codes available. Some of them include: 
CMBFAST (written in FORTRAN 77 <cit.>), CMBEASY (C++ <cit.>), 
CAMB (FORTRAN 90 <cit.>), CLASS (C <cit.>), and COSMOSIS <cit.> (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.






  §.§.§ 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 <cit.>.-
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 <cit.>.- 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 <cit.>, and can be downloaded in Reference <cit.>.
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 <cit.> or MCEvidence <cit.>. 

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


    P(D|θ,H) = P(D|θ,H)P(θ|H) = an(D|θ,H),


where n(D|θ,H) = N P(D|θ,H) with N the length of the chain. The code uses Bayesian inference to find a, the constant of proportionality (see Reference <cit.> for details). Once this constant has been found, the Bayesian evidence is given by
P(D|H) = aN.






§ 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:

    L∝exp[-∑_i (T(z_i)-T_i)^2/2σ_i],


where T(z_i) is the theoretical value related to 
the observation T_i; and σ_i are the corresponding errors associated 
to each measurement. In our estimation, T(z_i) is given by (<ref>) for Supernovae; (<ref>) for CMB; (<ref>), (<ref>), and (<ref>) for BAO; and (<ref>) for Cosmic Chronometers.

We use the BAO data mentioned in Section <ref> (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 ΛCDM model, already explored in 
Section <ref>, but now including the full observations presented in 
Section <ref> and using the SimpleMC code.



Here, in order to test the Friedmann Equation (<ref>) with the data,
we consider as free parameters (along with flat priors):
the total matter dimensionless density parameter Ω_m ∈ [0.05,1.5], the 
baryon physical density Ω_b h^2 ∈ [0.02,0.025], and the 
dimensionless Hubble constant h ∈ [0.4,1]. 

Then, assuming the base ΛCDM model, we let the curvature of the Universe be a 
free parameter
(model oΛCDM), with its corresponding flat prior Ω_k ∈ [-1.5,1.5].
Moreover, because the cosmological constant is only a particular case for the 
dark energy equation-of-state ω=-1, we let ω be a free parameter
with flat priors ω∈ [-2.0,0.0] and labeled it as model ωCDM.
We may combine the addition of curvature and constant 
ω to define the oω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 
ω_0 ∈ [-2.0,0.0] and ω_a ∈ [-2.0,2.0]
and labeled it as ω_0 ω_aCDM model. 
Again, we can incorporate the curvature 
of the Universe to the CPL parameterization, named as oω_0 ω_aCDM.

By using the combined dataset BBAO+Planck+SN+HD, Table <ref> 
shows the best fit values, along with 1σ 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 <ref>, for the ΛCDM model and 
by using only Hubble data, we had 
h=0.677± 0.313 and Ω_m=0.331 ± 0.0628.
Now, with the inclusion of BAO, Planck, and SN data, the constraints have improved considerably to h=0.684± 0.006 and Ω_m=0.299 ± 0.007.

Figure <ref> is updated with new data in order to get Figure <ref>.
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σ confidence regions obtained from our estimations. 

From Table <ref>, we also observe the best fit of most of the new 
parameters—additional to the base ΛCDM model—remained well 
inside the 1σ confidence level. 
However, the exception is for the oω_0ω_aCDM model, 
where Ω_k, ω_0, and ω_a lay down right outside the 1σ
region. This main feature is better observed in Figure <ref>, where 
the 2D posterior distributions, along with 1 and 2σ confidence regions,
are shown. Here, the standard ΛCDM values 
are marked with a dash line ( Ω_k=0, ω= ω_0 =-1, ω_a=0).
The inclusion of extra parameters improves the fit to the data, 
 observed through the minimum χ^2_ min, before the last row of Table <ref>.
However, it also carries out a penalization factor that affects directly the model selection, as seen in Reference <cit.>. 






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 ΛCDM compared to ωCDM and ω_0ω_aCDM; strong evidence with respect to oΛCDM and oω_0ω_aCDM; and decisive evidence against oωCDM. This is in agreement with the results presented in Reference <cit.>, where Planck 2013 found no evidence against ΛCDM.
 













§ 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 θ.

 Bayes' theorem states that 


    P(θ|D, H)=P(D|θ,H)P(θ|H)/P(D|H),


where the prior probability P(θ|H)—the state of knowledge
before acquiring the data—is upgraded through the likelihood
P(D|θ,H) when experimental data D are considered.

The aim for parameter estimation is to then obtain the posterior probability
P(θ|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(θ|D, H)∝ L(D|θ,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_0e^-χ^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 useful a dataset is may be by introducing the hyperparameter method.

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(D|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 (ΛCDM) to 
current observations and compare it to different models. We found the model that best fit the data, through χ^2_ min, corresponds to a curved Universe with a dynamical dark energy, namely oω_o ω_aCDM. It is important to clarify that this does not mean that the oω_o ω_aCDM model is the final model; it merely shows that, for these particular datasets, there is an improvement in the fit when compared to Λ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 ΛCDM.




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.
































JHEP




