Chapter 5
Methodological Effort and Epistemological Reflection

We have stressed earlier that ISD studies were essentially empirical, because of the complexity of their subject. Most of the time, they consist in interpreting data with formulae and models. Yet, comparing observations to models is a very wide methodological topic. It also has epistemological ¹ consequences. All the knowledge we derive about ISD depends on the way it was inferred. The question of the methods we use, and the way we articulate different results to build a comprehensive picture of the ISM, is thus of utmost importance.

5.1 Understanding the Opposition between Bayesians and Frequentists

Historically, two competitive visions of the way empirical data should be quantitatively compared to models have emerged, the (i) Bayesian, and (ii) frequentist methods. We personally follow the Bayesian method and will try to give arguments in favor of its superiority. An efficient way to present the Bayesian approach is to compare it to its alternative, and to show how both methods differ. There is a large literature on the subject. The book of Gelman et al. (2004) is a reference to learn Bayesian concepts and techniques. The posthumous book of Jaynes (2003) is more theoretical, but very enlightening. Otherwise, several reviews have been sources of inspiration for what follows (Jaynes, 1976; Loredo, 1990; Lindley, 2001; Bayarri & Berger, 2004; Wagenmakers et al., 2008; Hogg et al., 2010; Lyons, 2013; VanderPlas, 2014). A good introduction to frequentist methods can be found in the books of Barlow (1989) and Bevington & Robinson (2003).

5.1.1 Two Conceptions of Probability and Uncertainty

Bayesian and frequentist methods differ by: (i) the meaning they attribute to probabilities; and (ii)  the quantities they consider random. Their radically different approaches and the various bifurcations the two methods take to address a given problem originate from these sole conceptions.

5.1.1.1 The Concept of Conditional Probability

The concept of conditional probability is central to what follows. As we will see, Bayesian and frequentist approaches differ on this aspect.

The meaning of conditional probabilities. If $A$ and $B$ are two logical propositions, the conditional probability, noted $p\left(A|B\right)$, is the probability of $A$ being true, knowing $B$ is true. To give an astronomical example, let’s assume that:

• $A$ is the probability per unit time to observe a SN Ia, when pointing a telescope at a random star; and
• $B$ is the probability to observe a binary system, when pointing a telescope at a random star.

Since most stars are LIMS, and that they have a lifetime $⟨\tau ⟩\simeq 10$ Gyr, we can estimate the probability to observe a SN Ia, knowing we are observing a binary system:

On the contrary the probability to observe a binary system, knowing we are observing a SN Ia is:

because SN Ia happen only in binary systems. We see that $p\left(A|B\right)\ne p\left(B|A\right)$. In our example, the two quantities do not even have the same units.

All probabilities are conditional. In practice, all probabilities are conditional. In the previous example, we have implicitly assumed that our possibility space, that is the ensemble of cases we can expect out of the experiment we are conducting, contained all the events where we are actually observing a star, when we are pointing our telescope at its coordinates. However, what if there is suddenly a cloud in front of the telescope? We would then need to account for these extra possibilities, which is equivalent to adding conditions. For instance, if we are conducting the same experiment from a ground-based telescope in London, during winter, we will get:

where the symbol $\wedge$ denotes the logical “and”.

Bayes’ rule. Thomas BAYES derived, in the middle of the XVIII${}^{th}$ century, a formula to reverse the event and the condition, in conditional probabilities. Fig. 5.1 shows a Venn diagram, that is a graphic representation of a possibility space. We have shown an arbitrary event A, in red, and another event B, in blue. The intersection of A and B, numbered (2), can be written $A\wedge B$. The probability of A, knowing B, is the probability of A, when B is considered as the new possibility space. In other words:

We therefore have: $p\left(A\wedge B\right)=p\left(A|B\right)p\left(B\right)$. By symmetry of A and B, we have $p\left(A\wedge B\right)=p\left(B\wedge A\right)$, thus $p\left(A|B\right)p\left(B\right)=p\left(B|A\right)p\left(A\right)$, which gives Bayes’ rule:

5.1.1.2 The Bayesian and Frequentist Assumptions

We review here the assumptions of both methods, in a general, abstract way. We will illustrate this presentation with concrete examples in Sect. 5.1.2. Let’s consider we are trying to estimate a list of $n$ physical parameters, $\overrightarrow{x}={\left(x_i\right)}_{i=1,\dots <!--FUNCTION\; APPLICATION-->,n}$, using a set of $m$ observations, $\overrightarrow{d}={\left(d_j\right)}_{j=1,\dots <!--FUNCTION\; APPLICATION-->,m}$. Let’s also assume that we have a model, $f$, predicting the values of the observables, ${\overrightarrow{d}}_{mod}$, for any value of $\overrightarrow{x}$: ${\overrightarrow{d}}_{mod}=f\left(\overrightarrow{x}\right)$.

The Bayesian viewpoint. The Bayesian approach considers that there is an objective truth, but that our knowledge is only partial and subjective. Bayesians thus assume the following.

Probabilities

quantify the plausibility of a proposition, when we only have incomplete knowledge. With this assumption, probabilities can be assigned to parameters and hypotheses. More precisely, the true value of a parameter is considered fixed and an hypothesis is either true or false, but our knowledge of the value of this parameter, or of the plausibility of this hypothesis, can be described by random variables. Probabilities are therefore used to quantify our partial, subjective knowledge.

The method

then consists in sampling the probability distribution of the parameters, conditional on the data, using Bayes’ rule:

$$\underset{posterior}{\underbrace{p\left(\overrightarrow{x}|\overrightarrow{d}\right)}}\propto \underset{likelihood}{\underbrace{p\left(\overrightarrow{d}|\overrightarrow{x}\right)}}\times \underset{prior}{\underbrace{p\left(\overrightarrow{x}\right)}}.$$

(5.6)

Compared to Eq. ($5.5$), Eq. ($5.6$) misses the denominator, $p\left(\overrightarrow{d}\right)$. This is because this distribution is independent of our variables, that are the physical parameters. If we were to explicit it, it would be:

$$p\left(\overrightarrow{d}\right)=\int \; <!--nolimits-->p\left(\overrightarrow{d}|\overrightarrow{x}\right)p\left(\overrightarrow{x}\right)d^n\overrightarrow{x},$$

(5.7)

which is simply the normalization factor of $p\left(\overrightarrow{x}|\overrightarrow{d}\right)$. This factor can thus be estimated by numerically normalizing the posterior, hence the proportionality we have used in Eq. ($5.6$). The three remaining terms are the following.

The posterior distribution,

$p\left(\overrightarrow{x}|\overrightarrow{d}\right)$, is what we are interested in. It is literally the PDF of the physical parameters (what we want), knowing the data (what we have).

The likelihood,

$p\left(\overrightarrow{d}|\overrightarrow{x}\right)$, is the probability of the data, for a fixed value of the parameters. It can be computed using our model, $f\left(\overrightarrow{x}\right)$. For instance, assuming that our observations are affected by uncorrelated, Gaussian noise, with standard deviations $\overrightarrow{\sigma }={\left({\sigma }_j\right)}_{j=1,\dots <!--FUNCTION\; APPLICATION-->,m}$, we can write:

$$p\left(\overrightarrow{d}|\overrightarrow{x}\right)=\frac{1}{{\left(2\pi \right)}^{m∕2}\prod _{j=1}^m{\sigma }_j}\exp <!--FUNCTION\; APPLICATION-->\left(-\sum _{j=1}^m\frac{{\left(f_j\left(\overrightarrow{x}\right)-d_j\right)}^2}{2{\sigma }_j^2}\right).$$

(5.8)

The prior distribution,

$p\left(\overrightarrow{x}\right)$, is a unique feature of the Bayesian approach. It literally quantifies all our prior knowledge about the values of $\overrightarrow{x}$. We will give concrete examples of what that could mean in the following sections.

The result

of the Bayesian approach is the posterior, as it contains all the information we want on the parameters, informed by the observations and our prior knowledge. We can then decide to synthesize this information by, for instance:

• quoting moments of the posterior: means, $⟨x_i|\overrightarrow{d}⟩$, standard deviation, $\sigma \left(x_i|\overrightarrow{d}\right)$, etc.;
• quoting their correlation coefficient: $\rho \left(x_i,x_j|\overrightarrow{d}\right)$, etc.;
• testing hypotheses: $P\left(x_i>\text{const}|\overrightarrow{d}\right)$, etc.;
• and so on.

The frequentist viewpoint. The frequentist approach also considers that there is an objective truth, but it differs with the Bayesian viewpoint by rejecting its subjectivity. Frequentists thus assume the following.

Probabilities

are the limit to infinity of the occurrence frequency of an event, in a sequence of repeated experiments, under identical conditions ². This repeated event can be, for instance, the measure of a quantity. The rejection of the use of probabilities as a quantification of knowledge forbids frequentists to consider parameters or hypotheses as random variables. In other words, physical quantities have a single, true value, and hypotheses are either true or false. Only the data, which are tainted with uncertainties, can be considered as random variables.

The method

then consists in using the likelihood, $p\left(\overrightarrow{d}|\overrightarrow{x}\right)$, to perform several tests, the most well-known being the maximum likelihood. We can see here that frequentists consider the probability distribution of the data, given the physical parameters. The data are thus considered as variables, and the physical parameters, fixed.

The results

consist in describing what values of the physical parameters we would find if we were to repeat the experiment in the same conditions. Frequentists then compute the uncertainties on their estimate of the parameters by simulating new data that could have been obtained in the same conditions. We thus end-up with a distribution of parameter values, that is not a probability distribution. For instance, hypothesis testing can not be performed the Bayesian way, because we dot not have a conditional probability of the parameter. We will see in Sect. 5.1.4.1 that we need to resort to the infamous significance tests, instead.

5.1.2 Comparison of the Two Approaches on Simple Cases

We now compare the two approaches on a series of simple examples, in order to demonstrate in which situations the two approaches may differ.

5.1.2.1 Simple Case: When the Two Approaches Agree

Let’s assume we are trying to estimate the flux of a star, $F_{\star }$, in a given photometric filter, with the following assumptions.

• The true flux of the star is $F_{true}=42$ (in arbitrary units).
• We have performed $m=3$ repeated measures, $d_j=F_j$ $\left(j=1,\dots <!--FUNCTION\; APPLICATION-->,m\right)$.
• Each flux has Gaussian independent, identically distributed (iid) uncertainties: ${\sigma }_j={\sigma }_F=14$ $\left(\forall <!--FUNCTION\; APPLICATION-->j=1,\dots <!--FUNCTION\; APPLICATION-->,m\right)$.
• The physical quantity we want to estimate is simply the value of the flux, $F_{\star }$. This corresponds thus to the simplest case, where we have only $n=1$ parameter and the model is the identity: $F_{mod}=f\left(x\right)=x$.

This is represented in Fig. 5.2. There is an analytic solution to this simple case:

The Bayesian solution. The Bayesian solution is obtained by sampling the posterior distribution in Eq. ($5.6$): $p\left(F_{\star }|F_1,F_2,F_3\right)$.

The likelihood term

is exactly Eq. ($5.8$), in our case, as we have Gaussian iid uncertainties.

The choice of the prior

is the arbitrary part of the Bayesian model. For this first example, let’s assume that we have no idea what the flux should be. We thus take a flat, uninformative prior, $p\left(F_{\star }\right)\propto 1$. This is a subjective choice that has the following consequences.

1.

This particular choice is called an improper prior, meaning it can not be normalized, as: ${\mathrm{\int \; <!--nolimits-->}<!--FUNCTION\; APPLICATION-->}_{-\infty }^{\infty }p\left(F_{\star }\right)d{F}_{\star }=\infty$. In practice, we thus need to choose lower and upper bounds, $\left[F_{-},F_{+}\right]$, beyond which the prior is 0. Since we are measuring a positive quantity ³, we can take $F_{-}=0$ as the lower bound. The upper bound could be taken as several times the flux of a 120 $M_{\odot }$ star at the distance of the source. We thus have:

$$p\left(F_{\star }\right)=\left\{\begin{array}{cc}\frac{1}{F_{+}-F_{-}}\hfill & \text{for }{F}_{-}\le F_{\star }\le F_{+}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}\right.$$

(5.10)

2.

This prior is also subjective, as it depends on the choice of the variable. If we decide to study $\ln <!--FUNCTION\; APPLICATION-->F_{\star }$, instead of $F_{\star }$, it will lead to a different result. We will address this issue in Sect. 5.1.4.1. We can summarize it the following way.

• If the choice of the prior matters in the final solution, it means that the weight of evidence brought by the data is weak. It is therefore natural that the way we decide to quantify our prior knowledge is important. This is an aspect of the Bayesian method we need to embrace.
• On the opposite, if the weight of evidence brought by the data is large, the choice of the prior will not matter significantly. In other words, if the width of the likelihood is much narrower than a dex, the difference between multiplying by $p\left(F_{\star }\right)$ or $p\left(\ln <!--FUNCTION\; APPLICATION-->F_{\star }\right)$ will be negligible.

The solution is represented in Fig. 5.3.a. We have sampled the posterior using a Markov Chain Monte-Carlo method (MCMC; using the code emcee by  Foreman-Mackey et al., 2013). We will come back to MCMC methods in Sect. 5.1.3. The estimated value is $F_{\star }\simeq 53.9\pm 8.1$ ($1\sigma$ uncertainty). This is exactly the analytical solution in Eq. ($5.9$). The $95\%$ credible range, that is the range centered on $⟨F_{\star }⟩$ containing a $95\%$ probability, is ${CR}_{95\%}\left(F_{\star }\right)=\left[37.8,69.8\right]$.

The frequentist solution. There are different ways to approach this problem in the frequentist tradition. The most common solution would be to use a Maximum-Likelihood Estimation (MLE) method.

Maximum-likelihood value.

There are numerical tools to compute the MLE of complex models. We have used a Levenberg-Marquardt algorithm (e.g. Markwardt, 2009). The MLE value of $F_{\star }$ we derive is $F_{ML}=53.9$, which is in agreement with the analytical solution in Eq. ($5.9$). It however does not provide uncertainties.

Uncertainty estimates.

A common solution to estimate the uncertainties on $F_{ML}$ would be to perform bootstrapping. Following the frequentist conception of probabilities, we randomly draw new observations around the MLE value, a large number of times. We thus obtain a set of synthetic repeated measures, assuming the population distribution has $F_{ML}$ for mean. For each new set, we derive a new MLE value, ${\bar{F}}_{ML}^{\left(k\right)}$ $\left(k=1,\dots <!--FUNCTION\; APPLICATION-->,3000\right)$. The first four draws are represented in Fig. 5.3.b-e. The standard deviation of this sample is $8.1$, in agreement with Eq. ($5.9$). We can also compute the $95\%$ confidence interval, from this sample: ${CI}_{95\%}\left(F_{\star }\right)=\left[37.8,69.8\right]$.

A few remarks. We can see that both methods give the same exact result, which is also consistent with the analytical solution (Eq. $5.9$). This is because the assumptions were simple enough to make the two approaches equivalent:

1.

by assuming a flat prior, we removed the effect of this Bayesian peculiarity;

2.

the symmetry and the iid nature of the noise made the sampling of the likelihood as a function of the parameters, or as a function of the data, identical.

Note also our subtle choice of terminology: (i) we talk about credible range in the Bayesian case, as this term designates the quantification of our beliefs; while (ii) we talk about confidence interval in the frequentist case, as it concerns our degree of confidence in the results, if the experiment was repeated a large number of times, assuming the population mean is the maximum likelihood.

5.1.2.2 Benefits of Using an Informative Prior

A first way to find differences between the Bayesian and frequentist approaches is to explore the effect of the prior. To that purpose, let’s keep the same experiment as in Sect. 5.1.2.1, but let’s assume now that the star we are observing belongs to a cluster, and we know its distance.

The Bayesian improvement. Contrary to Sect. 5.1.2.1, where we had to guess a very broad, flat prior, we can now refine this knowledge, based on the expected luminosity function, scaled at the known distance of the cluster. This is represented on Fig. 5.4. The posterior distribution (yellow) is now the product of the likelihood (green) and prior (blue). The frequentist solution has not changed, as it can not account for this kind of information. We can see that the maximum a posteriori is now closer to the true value than the maximum likelihood. This can be understood the following way.

The true flux

has been drawn from the luminosity function, because we have randomly targetted a star in this cluster. This is the definition of the luminosity function, which we happen to have chosen as the prior. This is often noted $F_{true}\sim p\left(F_{\star }\right)$, the $\sim$ symbol meaning “distributed as”.

The observed flux

has then been drawn from a normal law centered on $F_{true}$ with variance ${\sigma }_F^2$, that can be written: $F_j|F_{true}\sim \mathcal{N}\left(F_{true},{\sigma }_F^2\right)$. This is equivalent to saying that the observed flux has been drawn from the posterior $p\left(F_{\star }\right)\times \mathcal{N}\left(F_{true},{\sigma }_F^2\right)$. Sampling this distribution is thus the best choice we can make, considering the information we have. This is why we get an advantage over the frequentist result.

If we perform several such measures, there will be some Bayesian solutions that will get corrected farther away from the true flux. This is a consequence of stochasticity. On Fig. 5.4, keeping our value of $F_{true}$, this will be the case if the noise fluctuation, $-\delta F$ is negative, that is if the observed flux is lower than $F_{true}$ ($F_j=F_{true}-\delta F<F_{true}$). However, this deviation will be less important than the correction we would benefit from if the same fluctuation was positive, because the prior would be higher: $p\left(F_{true}-\delta F\right)>p\left(F_{true}+\delta F\right)$ (cf. Fig. 5.4). The prior would therefore correct less the likelihood on the left side, in this particular case. Thus, on average, taking into account this informative prior will improve the results.

Accumulation of data. In this case and the previous one (Sect. 5.1.2.1), we had three observations of the same flux. The posterior distribution we sampled was:

This is because we considered the three measures as part of the same experiment. However, we could have chosen to analyze the data as they were coming. After the first flux, we would have inferred:

This posterior would have been wider (i.e. more uncertain), as we would have had only one data point. Note that such an inference would have not been possible with a frequentist method, as we would have had one parameter for one constraint (i.e. zero degree of freedom). What is interesting to note is that the analysis of the second measure, can be seen as taking into account the first measure in the prior:

and so on. For the third measure, the new prior would be $p\left(F_{\star }|F_1,F_2\right)$:

which is formally equivalent to Eq. ($5.11$), but is a different way of looking at the prior. Notice that the original prior, $p\left(F_{\star }\right)$, appears only once in the product. The more we accumulate data, the less important it becomes.

5.1.2.3 Case Where the Two Approaches Differ: Non Gaussianity and Few Data

The other reason why the two approaches might differ is because Bayesians sample $p\left(\overrightarrow{x}|\overrightarrow{d}\right)$, whereas frequentists produce a series of tests based on $p\left(\overrightarrow{d}|\overrightarrow{x}\right)$. The difference becomes evident when we consider non-Gaussian errors with small data sets.

Flux with a non-linear detector. Let’s assume that we are measuring again the flux from the same star ($F_{\star }=42$), with $m=3$ repetitions, but that our detector is now highly non-linear. This non-linearity translates into a heavily-skewed split-normal noise distribution (cf. Sect. C.2.2):

with $\lambda =0.3$ and $\tau =50$. This noise distribution is the red curve in Fig. 5.5.a. In practice, this could for example be a very accurate detector suffering from transient effects. The measured value would be systematically higher than the true flux ⁴. We have simulated three such measures in Fig. 5.5.a, in blue. This problem was adapted from example 5 of Jaynes (1976).

The solutions. Knowing that the measured flux is always greater than or equal to the true flux, it is obvious that the solution should be lower than the lowest measured flux: $F_{\star }\le \underset{j}{\min <!--FUNCTION\; APPLICATION-->}{F}_j=47.6$, in our particular case. This flux range, corresponding to inconsistent values, has been hatched in grey, in Fig. 5.5.

The Bayesian solution

is obtained the same way as before, by sampling the posterior. We again assume a flat prior, and take the likelihood in Eq. ($5.15$). The posterior is shown in yellow, in Fig. 5.5.b. It has zero probability in the range that we qualified as “inconsistent”, and the true value falls in a high probability domain. The mean and standard deviation of the posterior give us $F_{\star }\simeq 43.8\pm 3.6$, with ${CR}_{95\%}\left(F_{\star }\right)=\left[34.6,47.9\right]$.

The frequentist solution

is obtained the same way as before. The maximum likelihood and the whole bootstrapping sample however falls in the “inconsistent” domain. We get $F_{\star }\simeq 51.9\pm 4.0$, with ${CI}_{95\%}\left(F_{\star }\right)=\left[47.4,62.5\right]$. We see here that the frequentist solutions fails at inferring the true flux. On top of that, it gives only “inconsistent” solutions. In this particular case, this is because of the asymmetry of the noise, which breaks the symmetry between $p\left(F_j|F_{\star }\right)$ and $p\left(F_{\star }|F_j\right)$ that we had in Sect. 5.1.2.1. This can be seen in Fig. 5.5.b. The frequentist solution is the mirror symmetric of the Bayesian posterior for that reason.

When $m$ increases, the frequentist solution gets closer and closer to the true flux. However, a bootstrapping analysis will reject the true flux in $100\%$ of the cases.

The reason of the frequentist failure. The failure of the frequentist approach is a direct consequence of its conception of probability (cf. e.g. VanderPlas, 2014, for a more detailed discussion and more examples). The frequentist method actually succeeds in returning the result it pretends to give: predicting a confidence interval where the solution would fall 95 % of the time, if we repeated the same procedure a large number of times. This is however not equivalent to giving the credible range where the true value of the parameter has a 95 % probability to be (the Bayesian solution). With the frequentist method, we have no guarantee that the true flux will be in the confidence interval, only the solution. We can see that the main issue with the frequentist approach is that it is difficult to interpret, even in a simple problem such as that of Fig. 5.5. “Bayesians address the question everyone is interested in by using assumptions no-one believes, while Frequentists use impeccable logic to deal with an issue of no interest to anyone” (Lyons, 2013). In the previous citation, the “assumption no-one believes” is the subjective choice of the prior, and the “issue of no interest to anyone” is the convoluted way frequentists formulate a problem, to avoid assigning probabilities to parameters.

5.1.3 Numerical Methods to Solve Bayesian Problems

Bayesian problems are convenient to formulate as they consist in laying down all the data, the model, the noise sources and the nuisance variables to build a posterior, using Bayes’ rule. The Bayesian results are also convenient to interpret as they all consist in using the posterior, which gives the true probability of the parameters. However, in between, estimating the average, standard deviation, correlation coefficients of parameters, or testing hypotheses can be challenging, especially if there are a lot of parameters or if the model is complex. Fortunately, several numerical methods have been introduced to make these tasks simpler. Most of these methods are based on Markov Chain Monte-Carlo (MCMC ⁵), which are a class of algorithms for sampling PDFs.

5.1.3.1 Sampling the Posterior Distribution

Markov Chains Monte-Carlo. A MCMC draws samples from the posterior. In other words, it generates a chain of $N$ values of the parameters, ${\overrightarrow{x}}_k$ ($k=1,\dots <!--FUNCTION\; APPLICATION-->,N$). These parameter values are not uniformly distributed in the parameter space, but their number density is proportional to the posterior PDF. Consequently, there are more points where the probability is high, and almost none where the probability is low. It has several advantages.

• Contrary to grid-based sampling or standard Monte-Carlo techniques, we spend most of our computing time estimating our model (which can be costly) where it matters, and not much time where it does not.
• From the MCMC, it is very simple and easy to estimate moments of a given parameter, marginalizing over the other ones. For instance, if we have two parameters, $x$ and $y$, the average of $x$, marginalizing over $y$ would be:

  $$⟨x⟩\equiv \iint x.p\left(x,y\right)dxdy\simeq \frac{1}{N}\sum _{k=1}^N{x}_k,$$

  (5.16)

  where the second equality is simply the average of the sample. We could do the same for the standard deviation, or the correlation coefficient.

• It also makes hypothesis testing very easy. For instance, if we want to know the probability that $x>y$, we literally compute the fraction of points in the MCMC, where $x_k>y_k$.

All these operations would have been much more expensive, in terms of computing time, if we had to numerically solve the integral. In particular, computing the normalization of the whole posterior would have been costly. From a technical point of view, a MCMC is a random series of values of a parameter where the value at step $k$ depends only the value at step $k-1$. We briefly discuss below the two most used algorithms. A good presentation of these methods can be found in the book of Gelman et al. (2004) or in the Numerical recipes (Press et al., 2007).

The Metropolis-Hastings algorithm. To illustrate this method and the next one, let’s consider again the measure of the flux of our star, with the difference, this time, that we would be observing it through two different photometric filters.

• Let’s call $F_{\star }$ and $G_{\star }$ the fluxes in these two bands, with uncertainties ${\sigma }_F=12$ and ${\sigma }_G=10$.
• Let’s assume that the uncertainties in the filters are correlated, with a correlation coefficient, $\rho =0.8$, because, for instance, of the way the calibration was performed.
• Let’s assume that we make only $m=1$ measure in each band, with observed fluxes $F_1=42$ and $G_1=36$. The posterior, assuming a flat prior, is simply a bivariate normal distribution centered on the observed flux, with covariance matrix:

  $$V=\left(\begin{array}{cc}\hfill {\sigma }_F^2\hfill & \hfill \rho {\sigma }_F{\sigma }_G\hfill \\ \hfill \rho {\sigma }_F{\sigma }_G\hfill & \hfill {\sigma }_G^2\hfill \end{array}\right).$$

  (5.17)

  The posterior is thus, noting $\overrightarrow{x}=\left(F_{\star },G_{\star }\right)$ and $\overrightarrow{d}=\left(F_1,G_1\right)$:

  $$p\left(\overrightarrow{x}|\overrightarrow{d}\right)\propto \exp <!--FUNCTION\; APPLICATION-->\left(-\frac{1}{2}{\left(\overrightarrow{x}-\overrightarrow{d}\right)}^T{V}^{-1}\left(\overrightarrow{x}-\overrightarrow{d}\right)\right).$$

  (5.18)

  Contours of this distribution are represented in Fig. 5.6.a-b.

The algorithm proposed by Metropolis et al. (1953) and generalized by Hastings (1970) is the most popular method to sample any PDF. This is a rejection method, similar to what we have discussed for MCRTs, in Sect. 3.1.1.4 (cf. also Appendix C.2.3.1).

A proposal distribution,

$p\left({\overrightarrow{x}}_k|{\overrightarrow{x}}_{k-1}\right)$, first needs to be chosen. The choice of this distribution is instrumental in the sampling efficiency: (i) if it is too wide, a lot of draws will be rejected; (ii) if it is too narrow, the sampling steps are going to be small, and more iterations are going to be necessary to sample the posterior. For our present example, we choose a bivariate normal distribution, centered on ${\overrightarrow{x}}_{k-1}$, whose width is the noise of our data, $s_F={\sigma }_F$ and $s_G={\sigma }_G$:

$$p\left({\overrightarrow{x}}_k|{\overrightarrow{x}}_{k-1}\right)\propto \exp <!--FUNCTION\; APPLICATION-->\left[-\frac{1}{2}{\left(\frac{F_{\star ,k}-F_{\star ,k-1}}{s_F}\right)}^2-\frac{1}{2}{\left(\frac{G_{\star ,k}-G_{\star ,k-1}}{s_G}\right)}^2\right]$$

(5.19)

The method

then consists in the following steps, iterated $N$ times.

1.

At each iteration, $k$, we draw a new set of parameters, ${\overrightarrow{x}}_k$, from the proposal distribution: ${\overrightarrow{x}}_k\sim p\left({\overrightarrow{x}}_k|{\overrightarrow{x}}_{k-1}\right)$.

2.

With this new value, we compute the acceptance probability, defined as:

$${\alpha }_k=\min <!--FUNCTION\; APPLICATION-->\left(1,\frac{p\left({\overrightarrow{x}}_k|\overrightarrow{d}\right)p\left({\overrightarrow{x}}_{k-1}|{\overrightarrow{x}}_k\right)}{p\left({\overrightarrow{x}}_{k-1}|\overrightarrow{d}\right)p\left({\overrightarrow{x}}_k|{\overrightarrow{x}}_{k-1}\right)}\right).$$

(5.20)

The $\min <!--FUNCTION\; APPLICATION-->$ function is here to make sure we get a result between 0 and 1. In our case, we have chosen a symmetric proposal distribution. Eq. ($5.19$) thus implies that:

$$\frac{p\left({\overrightarrow{x}}_{k-1}|{\overrightarrow{x}}_k\right)}{p\left({\overrightarrow{x}}_k|{\overrightarrow{x}}_{k-1}\right)}=1\Rightarrow {\alpha }_k=\min <!--FUNCTION\; APPLICATION-->\left(1,\frac{p\left({\overrightarrow{x}}_k|\overrightarrow{d}\right)}{p\left({\overrightarrow{x}}_{k-1}|\overrightarrow{d}\right)}\right).$$

(5.21)

We just need to estimate our posterior at one position (assuming we saved the value of $p\left({\overrightarrow{x}}_{k-1}|\overrightarrow{d}\right)$, after the previous iteration). In addition, since we need only the ratio of two points in the posterior, we do not need to normalize it. This is the reason why this algorithm is so efficient.

3.

To update ${\overrightarrow{x}}_k$, we draw a random variable, ${𝜃}_k$, uniformly distributed between 0 and 1.

• If ${𝜃}_k\le {\alpha }_k$, we accept the new value, ${\overrightarrow{x}}_k$.
• If ${𝜃}_k>{\alpha }_k$, we reject the new value and keep the old one, ${\overrightarrow{x}}_k={\overrightarrow{x}}_{k-1}$.

The initial value of the chain has to be a best guess.

The sampling of Eq. ($5.18$) with the Metropolis-Hastings method is shown in Fig. 5.6.a.

Gibbs sampling. The number of parameters, $n$, determines the dimension of the posterior. The higher this number is, the smaller the support of the function is (i.e. the region where the probability is non negligible). Metropolis-Hastings methods therefore will have a high rejection rate, if $n»1$, requiring longer chains to ensure convergence. Gibbs sampling (Geman & Geman, 1984) provides an alternative MCMC method, where all draws are accepted. Its drawback is that it requires normalizing, at each iteration, the full conditional distribution:

that is the posterior fixing all parameters except one. This is only a one dimensional PDF, though, much less computing-intensive than the whole posterior. The method consists, at each iteration $k$, to cycle through the different parameters, and draw a new value from Eq. ($5.22$):

Since this distribution has an arbitrary form, random drawing can be achieved numerically using the CDF inversion method (Appendix C.2.3.2). Fig. 5.6.b represents the Gibbs sampling of Eq. ($5.18$). The squared pattern comes from the fact that we alternatively sample each parameter, keeping the other one fixed.

5.1.3.2 Post-Processing MCMCs

Assessing convergence. One of the most crucial questions, when using a MCMC method, is how long a chain do we need to run. To answer that question, we need to estimate if the MCMC has converged toward the stationary posterior. Concretely, it means that we want to make sure the sampling of the posterior is homogeneous, and that the moments and hypothesis testing we will perform will not be biased, because some areas of the parameter space have only been partially explored. The reason why the sampling of the parameter space might be incomplete is linked to the two following factors.

The burn-in

refers to the first drawn values, before the MCMC could find the support of the posterior. This can be seen in Fig. 5.6.a-b (highlighted in cyan). The arbitrary starting value (magenta star) is outside the $5\sigma$ contour of the PDF. The MCMC thus walks a few steps before finding the probable region. This burn-in phase is also highlighted in cyan, in Fig. 5.6.c-d. This burn-in phase could actually be significantly longer, for several non-exclusive reasons: (i) a larger number of parameters; (ii) a more degenerate model, with several local maxima; (iii) a less lucky choice of initial conditions; or (iv) a very well-constrained model, resulting in a very small support over the whole parameter space. There is no universal method to identify burn-in, it needs to be investigated carefully, most of the time. Running several MCMCs, starting from initial conditions distributed over the whole parameter space, is usually efficient.

The autocorrelation

of the MCMC results from the fact that the parameter value at step $k+1$ depends on step $k$. If several successive iterations stay in the same region of the posterior, this will create a portion of correlated values. The AutoCorrelation Function (ACF; e.g.  Sokal, 1996) is an essential tool to determine the correlation length of a MCMC. The ACF, $\widehat{\rho }$, of a given parameter, depends on the lag, $\tau$, that is the number of steps between two arbitrary iterations:

$$\widehat{\rho }\left(\tau \right)\equiv \frac{N}{N-\tau }\frac{\sum _{k=1}^{N-\tau }\left(x_k-⟨x⟩\right)\left(x_{k+\tau }-⟨x⟩\right)}{\sum _{k=1}^N{\left(x_k-⟨x⟩\right)}^2}.$$

(5.24)

This is the correlation coefficient of the parameter with itself, shifted by $\tau$ steps. The ACFs of our example are displayed in Fig. 5.6.e-f. We can see that the ACF starts at 1, for $\tau =0$. It then drops over a few steps and oscillates around 0. The typical lag after which the ACF has dropped to 0, corresponds to the average number of steps necessary to draw independent values. This typical lag can be quantified, by the integrated autocorrelation time, $t_{int}$ ⁶:

$$t_{int}\equiv 1+2\sum _{i=\tau }^N\widehat{\rho }\left(\tau \right).$$

(5.25)

It is represented in Fig. 5.6.e-f. It corresponds roughly to the average number of steps needed to go from one end of the posterior to the other. Different parameters of a given MCMC can in principle have very different $t_{int}$ (e.g. Galliano, 2018). To make sure that our posterior is properly sampled, we thus need to let our MCMC run a large number of steps, times $t_{int}$, after burn-in. The effective sample size, $N_{e\mathrm{ff}}\equiv N∕t_{int}$, quantifies the effective number of steps that can be considered independent. We need $N_{e\mathrm{ff}}»1$.

With the Metropolis-Hastings algorithm, the integrated autocorrelation time will depend heavily on the choice of the proposal distribution. We have explored the effect of the width of this distribution on $t_{int}$. In Eq. ($5.19$), instead of taking $s_F={\sigma }_F$ and $s_G={\sigma }_G$, we have varied this parameter. Fig. 5.7.a represents the mean rejection rate as a function of $\sqrt{s_F{s}_G∕\left({\sigma }_F{\sigma }_G\right)}$.

• When $\sqrt{s_F{s}_G∕\left({\sigma }_F{\sigma }_G\right)}«1$, the proposal is much narrower than the posterior. The proposed steps are thus very small, and a lot of them are necessary to cross the posterior. This is why $t_{int}$ is very large in this case (cf. Fig. 5.7.b).
• When $\sqrt{s_F{s}_G∕\left({\sigma }_F{\sigma }_G\right)}»1$, the proposal distribution is much larger than the posterior. Most proposals therefore falls outside the support of the posterior. They are therefore rejected close to 100 % of the times. This is why $t_{int}$ is also very large in this case.

Fig. 5.7.b shows that the only range where $t_{int}$ is reasonable is when the width of the proposal distribution is comparable to the width of the posterior.

Parameter inference. Numerous quantities can be inferred from a MCMC. We have previously seen that the average, uncertainties, and various tests can be computed using the posterior of the parameters of a source. This becomes even more powerful when we are analyzing a sample of sources. To illustrate this, let’s assume we are now observing $N_{\star }=5$ stars, through the same photometric bands as before. Fig. 5.8.b shows the posterior PDF of the two parameters of the five stars. It is important to distinguish the following two types of distributions.

The posterior of individual stars

are represented in Fig. 5.8. The error bars in panel (b) correspond to the mean and standard deviation of the marginal distributions in panels (a) and (c). They represent the uncertainty on the measured fluxes of each individual star. They are the moments of a given parameter, over the whole MCMC. This is what we have focussed on, until now.

The statistic distribution across the sample

is represented in Fig. 5.9. In panel (a), we have shown the distribution of the standard deviation of the sample, at each step $k$ in the MCMC:

$$\sigma {\left(F_{\star }\right)}_k\equiv \frac{1}{N_{\star }-1}\sum _{i=1}^{N_{\star }}{\left(F_{i,k}-\frac{1}{N_{\star }}\sum _{j=1}^{N_{\star }}{F}_{j,k}\right)}^2,$$

(5.26)

where $F_{i,k}$ is the observed flux of the star $i$, at the MCMC iteration $k$. It is how we can quantify the dispersion of the sample. We can thus quote the sample dispersion as (cf. Fig. 5.9.a):

$$\sigma \left(F_{\star }\right)\simeq ⟨\sigma \left(F_{\star }\right)⟩\pm \sigma \left[\sigma \left(F_{\star }\right)\right].$$

(5.27)

In our example, we have: $\sigma \left(F_{\star }\right)\simeq 13.8\pm 1.6$ and $\sigma \left(G_{\star }\right)\simeq 6.8\pm 1.0$. We can see that these values correspond roughly to the intrinsic scatter between individual stars, in Fig. 5.8.b, but they are larger than the uncertainty on the flux of individual stars. We can do the same for the correlation coefficient, as shown in Fig. 5.9.b: $\rho \left(F_{\star },G_{\star }\right)\simeq -0.63\pm 0.15$. Notice that it is negative, because the correlation between stellar fluxes, in Fig. 5.8.b, points toward the lower right corner. However, the correlation between the uncertainties on $F_{\star }$ and $G_{\star }$ is positive: the individual ellipses point in the other direction, toward the upper right corner. We have deliberately simulated data with these two opposite correlations to stress the difference between the individual likelihood properties and those of the ensemble. Finally, we could have done the same type of estimate for the mean of the sample: $⟨F_{\star }⟩\simeq 52.2\pm 1.5$ and $⟨G_{\star }⟩\simeq 26.0\pm 1.0$.

Quantifying the goodness of a fit. It is important, in any kind of model fitting, to be able to assess the quality of the fit. In the frequentist approach, this is done with the chi-squared test, which is limited in its assumptions to normal iid noise, without nuisance parameters. In the Bayesian approach, the same type of test can be done, accounting for the full complexity of the model (non-Gaussian errors, correlations, nuisance parameters, priors). This test is usually achieved by computing posterior predictive $p$-values (ppp; e.g. Chap. 6 of  Gelman et al., 2004). To illustrate how ppps work, let’s consider now that we are observing the same star as before, through four bands (R, I, J, H) and are performing a blackbody Bayesian fit to this SED, varying the temperature, $T_{\star }$, and the dilution factor, ${\Omega }_{\star }$. This is represented in Fig. 5.10.b. The principle is the following.

1.

We generate a set of replicated data, ${\overrightarrow{d}}_{rep}$, from our posterior:

$$p\left({\overrightarrow{d}}_{rep}|\overrightarrow{d}\right)\equiv \int \; <!--nolimits-->p\left({\overrightarrow{d}}_{rep}|\overrightarrow{x}\right)p\left(\overrightarrow{x}|\overrightarrow{d}\right)d\overrightarrow{x}.$$

(5.28)

If we sampled our posterior with a MCMC, this integral can simply be computed by evaluating our model (the blackbody, in the present case), for values of our drawn parameters: ${\overrightarrow{d}}_{rep}={\left\{f\left({\overrightarrow{x}}_k\right)\right\}}_{k=1}^N$.

2.

We evaluate the mean and standard deviation of this replicated data set, $⟨{\overrightarrow{d}}_{rep}|\overrightarrow{x}⟩$ and $\sigma \left({\overrightarrow{d}}_{rep}|\overrightarrow{x}\right)$. These quantities are the average and the dispersion of the predicted flux in the different bands. They will serve as position and scale references, when comparing model and observations.

3.

We compute a discrepancy metric, $T\left(\overrightarrow{d}|\overrightarrow{x}\right)$. Several choices are possible, but the most common is to adopt a chi-squared equivalent:

$$T\left(\overrightarrow{d}|\overrightarrow{x}\right)\equiv \sum _{j=1}^m\frac{{\left[d_j-⟨d_j|\overrightarrow{x}⟩\right]}^2}{\sigma {\left(d_j|\overrightarrow{x}\right)}^2}.$$

(5.29)

We compute this quantity both for the replicated set, $T\left({\overrightarrow{d}}_{rep}|\overrightarrow{x}\right)$, which is the blue distribution in Fig. 5.10.e, and for the observations, $T\left(\overrightarrow{d}|\overrightarrow{x}\right)$, which is the red line in Fig. 5.10.e (it is a single value). To be clear, only the data term in Eq. ($5.29$) changes between $T\left({\overrightarrow{d}}_{rep}|\overrightarrow{x}\right)$ and $T\left(\overrightarrow{d}|\overrightarrow{x}\right)$. The quantities $⟨d_j|\overrightarrow{x}⟩$ and $\sigma \left(d_j|\overrightarrow{x}\right)$ are identical in both cases.

4.

The quality of the test is assessed by comparing both quantities. To that purpose, we compute the following probability:

$$p_B\equiv p\left(T\left({\overrightarrow{d}}_{rep}|\overrightarrow{x}\right)\ge T\left(\overrightarrow{d}|\overrightarrow{x}\right)|\overrightarrow{d}\right).$$

(5.30)

If the difference between the replicated set and the data is solely due to statistical fluctuations, we should have on average $p_B\simeq 50\%$. The fit is considered bad, at the $95\%$ level, if $p_B<2.5\%$ or $p_B>97.5\%$.

We have illustrated this test in Fig. 5.10, varying the number of parameters and observational constraints, in order to explore the different possible cases.

A good fit

is shown in Fig. 5.10.b. We have varied both $T_{\star }$ and ${\Omega }_{\star }$ to fit the 4 fluxes. Fig. 5.10.e shows that $T\left(\overrightarrow{d}|\overrightarrow{x}\right)$ falls in the high probability range of $T\left({\overrightarrow{d}}_{rep}|\overrightarrow{x}\right)$. In other words, the average deviation of the replicated data, relative to the reference we have chosen, $⟨{\overrightarrow{d}}_{rep}|\overrightarrow{x}⟩$, is comparable to the deviation of the actual data relative to the same reference. The observations could thus have likely been drawn from our posterior.

A poor fit

is shown in Fig. 5.10.a. We have intentionally fixed the temperature of the fit at $T_{\star }=7000$ K, while the true value is $T_{\star }=6000$ K. We see that $p_B\simeq 1$, in Fig. 5.10.d. The difference between the observations and the model, can thus not be explained by the scatter of the model. It is the sign of a bad fit. In our case, it is because the model is bad (wrong choice of fixed temperature).

An overfit

is shown in Fig. 5.10.c. This time we fit only two fluxes. With a chi-squared fit, we would have 0 degrees of freedom. We see that $p_B\simeq 0$, in Fig. 5.10.f. The average model therefore gets too close to the observations. The ppp tells us this is very unlikely. It is however not an issue in terms of derived parameters. We can see that the true model (green) is well among the sampled model (blue), and the inferred parameters are consistent with their true values.

There is no issue with fitting even fewer constraints than parameters, with a Bayesian approach. The consequence is that the posterior is going to be very wide along the dimensions corresponding to the poorly constrained parameters. However, the results will be consistent, and the derived probabilities will be meaningful.

5.1.4 Decision Making and Limitations of the Frequentist Approach

We now resume our comparison of the Bayesian and frequentist approaches, started in Sect. 5.1.2. We focus more on the interpretation of the results and synthesize the advantages and inconveniences of both sides.

5.1.4.1 Hypothesis Testing

Until now, we have seen how to estimate parameters and their uncertainties. It is however sometimes necessary to be able to make decisions, that is to choose an outcome or its alternative, based on the observational evidence. Hypothesis testing consists in assessing the likeliness of a null hypothesis, noted $H_0$. The alternative hypothesis is usually noted $H_1$, and satisfies the logical equation $H_1=\neg <!--FUNCTION\; APPLICATION-->H_0$, where the $\neg <!--FUNCTION\; APPLICATION-->$ symbol is the logical negation. The priors necessarily obey $p\left(H_0\right)+p\left(H_1\right)=1$. To illustrate this process let’s go back to our first example, in Sect. 5.1.2. We are observing a star with true flux $F_{true}$, $m$ times, with an uncertainty ${\sigma }_F$ on each individual flux measurements, $F_i$. This time, we want to know if $F_{true}\le F_{test}$, for a given $F_{test}$.

Bayesian hypothesis testing. Bayesian hypothesis testing consists in computing the posterior odds of the two complementary hypotheses:

The posterior odds is the ratio of the posterior probabilities of the two hypotheses. It is literally the odds we would use for gambling (e.g. a posterior odd of 3 corresponds to a 3:1 odd in favor of $H_1$). The important term in Eq. ($5.31$) is the Bayes factor, usually noted $B{F}_{10}\equiv p\left(\overrightarrow{d}|H_1\right)∕p\left(\overrightarrow{d}|H_0\right)$. It quantifies the weight of evidence, brought by the data, against the null hypothesis. It tells us how much our observations changed the odds we had against $H_0$, prior to collecting the data. Table 5.1 gives a qualitative scale to decide upon Bayes factors. We see that it is a continuous credibility range going from rejection to confidence. The posterior of our present example, assuming a wide flat prior, is:

where $⟨F⟩={\sum <!--FUNCTION\; APPLICATION-->}_{i=1}^m{F}_i∕m$. The posterior probability of $H_0=\left(F_{true}\le F_{test}\right)$ is then simply:

It is represented in Fig. 5.11.a. This PDF is centered in $⟨F⟩$, since, as usual in the Bayesian approach, it is conditional on the data. Fig. 5.11.a shows the complementary posteriors of $H_0$ (red) and $H_1$ (blue), which are the incomplete integrals of the PDF. When we vary $F_{test}$, the ratio of the two posteriors, $B{F}_{10}$, changes. Assuming we have chosen a very wide, flat prior, such that $p\left(H_1\right)∕p\left(H_0\right)\simeq 1$, the Bayes factor becomes:

Fig. 5.11.b represents the evolution of the Bayes factor as a function of the sample size, $m$. In this particular simulation, $H_0$ is false. We see that, when $m$ increases, we accumulate evidence against $H_0$, going through the different levels of Table 5.1. The evidence is decisive around $m\simeq 60$, here.

Frequentist hypothesis testing. Instead of computing the credibility of $H_1$ against $H_0$, the frequentist approach relies on the potential rejection of $H_0$. It is called Null Hypothesis Significance Test (NHST; e.g.  Ortega & Navarrete, 2017). It consists in testing if the observation average, $⟨F⟩$, could have confidently been drawn out of a population centered on $F_{test}$. In that sense, it is the opposite of the Bayesian case. This is demonstrated in Fig. 5.11.c. We see that the distribution is centered on $F_{test}$, which is assumed fixed, and the observations, $⟨F⟩$, are assumed variable. Since frequentists can not assign a probability to $F_{test}$, they estimate a $p$-value, that is the probability of observing the data at hands, assuming the null hypothesis is true. It is based on the following statistics:

The $p$-value of this statistics is then simply:

Notice it is identical to Eq. ($5.33$), but different from the Bayes factor (Eq. $5.34$). It however does not mean the same thing. A significance test at $p_0=0.05$ does not tell us that the probability that the null hypothesis is $5\%$. It means that the null hypothesis will be rejected $5\%$ of the time ⁷. NHST decision making is represented in Fig. 5.11.c. It shows one of the most common misconceptions about NHST: the absence of rejection of $H_0$ does not mean that we can accept it. It just mean that the results are not significant. Accepting $H_0$ requires rejecting $H_1$, and vice versa. Fig. 5.11.d shows the effect of sample size on the $p$-value, using the same example as we have discussed for the Bayesian case. The difference is that the data are not significant until $m\simeq 70$.

The Jeffreys-Lindley’s paradox. Although the method and the interpretations are different, Bayesian and frequentist tests give consistent results, in numerous applications. There are however particular cases, where both approaches are radically inconsistent. This ascertainment was first noted by Jeffreys (1939) and popularized by Lindley (1957). Lindley (1957) demonstrated the discrepancy on an experiment similar to the example we have been discussing in this section, with the difference that a point null hypothesis is tested: $H_0=\left(F_{true}=F_{test}\right)$. Lindley (1957) shows that there are particular cases, where the posterior probability of $H_0$ is $1-\alpha$, and $H_0$ is rejected at the $\alpha$ level, at the same time. This “statistical paradox”, known as the Jeffreys-Lindley’s paradox has triggered a vigorous debate, that is still open nowadays (e.g. Robert, 2014). The consensus about the paradox is that there is no paradox. The discrepancy simply arises from the fact that both approaches answer different questions, as we have been illustrating at several occasions in this chapter, and that these different interpretations can sometimes be inconsistent.