I am in a context of forecasts in astrophysics. Don't be too rude if questions seem to you stupid or naive but rather indulgent, I am just looking for better undertsand all these numerical methods of Monte-Carlo alone/ Monte-carto coupling with Markov-Chain and the difference between a sampler and an estimator. This is little the mess in my head to grasp all the subtilities.

1. Using Covariance matrix at each step

In the following figure below below, especially in the central box I don't understand why I have to use the Covariance matrix at each call of a point that will be or not accepted in the distribution of the posterior : Is it done to compute the $\chi^2$ at each time and accept/reject it relying on some threshold, but on which criterion ?

In my code, I generate Power matter spectrum (in Cosmology at the upper left of the figure). Up to this, there is not random process. For me, this is in the central box that there is random with the computation of a posterior distribution with the formula :

$P(\Theta | data)=\dfrac{P(data | \Theta) \times P(\Theta)}{P(data)}$

As you can see, I need the Likelihood which directly depends of the theorical model, doesn't it ?

Then, I generate a sample of the Likelihood by taking random data in this likelihood ? I am a bit lost as you can see, mixing the 2 concepts and where the random processes occur.

Principe of "Inverse Physics problem

2. Monte-Carlo and Metropolis Hastings

Have I got to consider the term "Monte-Carlo" as a general way to generate distributions (or samples, I don't know which one of the two terms I must use (even if, with Monte-Carlo, I can compute and so estimate the expectation of a random variable knowing the PDF with an integral ?

And coupled with Metropolis-Hasting, the result is that we have a distribution of the posterior, from we can extract for example the mean (peak of the distribution) ?

3). Link between Likelihood and chi-squared : which is the deep link between Likelihood and chi-squared into Monte-Carlo Markov-Chain ?

4. Fisher formalism :

A last question : I heard that Fisher formalism could be only applied under the assumption that posterior/likelihood must be Gaussian.

Could anyone explain why ? and mostly, how to demonstrate it from a mathematical point of view ?

And if by lack of chance, the likelihood produced by a theorical model is not Gaussian, which other alternatives are possible to estimate a set of parameters ? Are there only Monte-Carlo-Markov-Chain methods which could circumvent the non-existence of Gaussian property of Likelihood ?

PS : I have asked different questions but all of them is linked in the sense they have connections between themselves from estimations and sampling method point of view.

So don't be too rude, I am just looking for trying better understand and grasp all the subtilities of all these concepts.

Even if I could have only one answer about one of my questions, I would be grateful.

  • 5
    $\begingroup$ Maybe you would get answers if you like your posts to one question. $\endgroup$
    – nicoguaro
    Jul 19, 2020 at 15:24
  • $\begingroup$ @nicoguaro Thanks but all the questions are linked between themselves, If I can't get an answer for each question, this doesn't matter, even partial or only one answer about a question would be enough for my understanding. $\endgroup$
    – user14831
    Jul 20, 2020 at 6:17
  • 1
    $\begingroup$ I understand. Was trying to help you. $\endgroup$
    – nicoguaro
    Jul 20, 2020 at 6:21
  • $\begingroup$ I opened a bounty, I tried to summarize my lack of understanding about all my questions. Sorry if I have been a little long in the detailed description of my bounty. Thanks $\endgroup$
    – user14831
    Jul 20, 2020 at 6:25

2 Answers 2


As I understand, your ultimate goal is to solve an inverse problem (i.e., infer some parameters from given data / observations). To this end, you want to apply Bayesian Inference, which relates the posterior (i.e., the probability distribution of the unknown parameters) to the likelihood (i.e., the probability model of observing some values given the parameters) and the prior (i.e., the probability distribution of your belief that the parameters attain some values). The evidence is only used to normalize in order to obtain a valid probability distribution (there are more use cases, e.g., model selection).

Since you're mentioning $\chi^2$, I'm supposing the likelihood looks like $$ p(\vec{y} | \vec{p}) \simeq \exp\left( -\frac{1}{2} (\vec{y} - \vec{p})^T \Sigma^{-1} (\vec{y} - \vec{p} ) \right), $$ which means that the data / observations $\vec{y}$ follow a normal distribution $\vec{y} \sim \mathcal{N}(\vec{p}, \Sigma)$ where the parameters $\vec{p}$ are the mean and the covariance $\Sigma$ is fixed. Note that the likelihood is just some function that can be (numerically) evaluated given the inputs $\vec{y}$ and $\vec{p}$.

Now, to infer the parameters, we are often interested in some functionals of the posterior. For example, mean, mode, standard deviation, quantiles, highest-posterior-density regions etc. Note that, for appreciating the Bayesian framework, the parameter inference should not be reduced to a single value (e.g., the mean of the posterior).

In this context, the Monte Carlo method essentially means to draw samples from the posterior and use a statistical estimator to infer some quantity (functionals such as mean, quantiles, etc.) from the distribution. That is, using the Monte Carlo method, we would simply need to draw random samples from the posterior and use this to estimate the parameters (i.e., take the sample mean to approximate the mean of the distribution). However, directly sampling from the posterior is usually not possible. In the example above (Likelihood is normal distribution) it depends on the choice of the prior distribution whether we obtain some known distribution for the posterior that can be sampled from directly (see conjugate priors).

As the name implies, Markov Chain Monte Carlo methods are a subset of Monte Carlo methods. It is a special method to generate samples from the posterior distribution, which can subsequently be used in a Monte Carlo estimator. The "standard" MCMC method is Metropolis-Hastings which works like this:

Given some initial state $\vec{p}_i$, perform the following steps:

  1. Draw a proposal $\vec{x} \sim Q(\vec{p}_i)$, where $Q$ is a probability distribution that may depend on $\vec{p}_i$.
  2. Calculate acceptance probability $$ \alpha_i = \min\left\{1, \frac{p(\vec{x} | \vec{y}) q(\vec{p}_i | \vec{x})}{p(\vec{p}_i | \vec{y}) q(\vec{x} | \vec{p}_i)} \right\}, $$ where $q(\cdot | \vec{a})$ is the density of $Q(\vec{a})$.
  3. Draw a random sample $u_i$ from the uniform distribution $U([0,1])$ on $[0,1]$ and set $$ \vec{p}_{i+1} = \begin{cases} \vec{x} & \text{if } u_i \leq \alpha_i \\ \vec{p}_i & \text{otherwise}. \end{cases} $$

In this algorithm, the posterior density $$p(\vec{p} | \vec{y}) \simeq p(\vec{y} | \vec{p}) p(\vec{p}) $$ without normalization is used. This involves the computation of the likelihood and prior at the proposed point $\vec{x}$, which, in turn, requires multiplication by the covariance matrix in the evaluation of the likelihood.

This should answer your first two questions.

  1. Link between Likelihood and chi-squared

This really depends on the modeling assumptions and the form of the likelihood. In the model used above, it is assumed that $$ \vec{y} = \vec{p} + \varepsilon, \qquad \varepsilon \sim \mathcal{N}_{\vec{0}, \Sigma}. $$ If the errors are not assumed to be Gaussian, the $\chi^2$ term would not appear in the likelihood.

  1. Fisher formalism

As far as I know, the maximum likelihood theory and the Fisher information do not depend on Gaussian distributions. They are fully generic.

And if by lack of chance, the likelihood produced by a theorical model is not Gaussian, which other alternatives are possible to estimate a set of parameters ?

Besides Monte Carlo methods (including MCMC), you can still apply maximum likelihood estimators for the model parameters.

  • $\begingroup$ Thanks for your anwser, I begin to see clearer but just one question : Sorry, I am stuck on a common english expression : Does "draw a sample" simply mean "generate a distribution" ? For example, I have seen many times "draw a sample", "draw a posterior", etc : it is just saying that we generate histogram, i.e a numerical distribution, isn't it ? $\endgroup$
    – user14831
    Jul 22, 2020 at 14:07
  • $\begingroup$ "Draw a sample" means to generate a random sample from a distribution. For example, in many languages, drawing a sample from a Gaussian $N(0,1)$ is done by randn() . It means that a random number (or vector) is generated such that, if repeated for a large number of times, the collection of samples approximates the probability distribution (i.e., same density or cumulative distribution function). Some distributions (e.g., normal, uniform, beta etc.) can be sampled easily using (pseudo) random number generators. Others require methods such as MCMC or importance/rejection sampling. $\endgroup$
    – cos_theta
    Jul 23, 2020 at 14:00
  • $\begingroup$ Fine to explain quickly. I would like to get further information about how you get the expression : $p(\vec{y} | \vec{p}) = \exp\left( -\frac{1}{2} (\vec{y} - \vec{p})^T \Sigma^{-1} (\vec{y} - \vec{p} ) \right),$ . I know the classical formula of $\chi^2=\sum_i\dfrac{(y_i-p)^2}{\sigma_{i}^2}$ but I have difficulties to pass to matricial form of $\chi^2$. Have you got a simple demonstration of this matricial expression ? I mean, involving the inverse of covariance matrix $\Sigma^{-1}$. Regards $\endgroup$
    – user14831
    Jul 24, 2020 at 13:05
  • $\begingroup$ Let's assume that the errors are independent such that $y_i\sim\mathcal{N}_{p_i, \sigma_i^2}$ for all observations $i=1,\dots,N$, where $\sigma_i^2$ is the variance. Then, the covariance matrix is diagonal:$$\Sigma = \mathrm{diag}\left(\sigma_1^2,\dots,\sigma_N^2\right). $$ Hence, we have $$\begin{aligned}\exp\left(-\frac{1}{2}(\vec{y}-\vec{p})^T\Sigma^{-1}(\vec{y}-\vec{p})\right)&=\exp\left(-\frac{1}{2}\sum_{i=1}^N \left(y_i-p_i\right)\frac{1}{\sigma_i^2}\left(y_i-p_i\right)\right) \\ &= \exp\left(-\frac{1}{2}\sum_{i=1}^N \frac{\left(y_i-p_i\right)^2}{\sigma_i^2} \right).\end{aligned}$$ $\endgroup$
    – cos_theta
    Jul 24, 2020 at 14:20
  • $\begingroup$ so the link between $\text{log}\mathcal L$ and $\chi^2$ would be some kind : $\text{log}\mathcal L= \sum_{i=1}^{N}\,p(y_i|p) \propto -\chi^2$, wouldn't it ? Or maybe simply we can write : $$\text{log}\mathcal L = -\dfrac{1}{2}\,\chi^2$$. So minimizing (or maximizing the Likelihood, I don't know since MLE says that derivate must be equal to zero) the $\text{log}\mathcal L$ is the same thing that minimizing the $\chi^2$ : could you confirm that ? $\endgroup$
    – user14831
    Jul 24, 2020 at 14:45

The previous answer pretty sums up my understanding on this problem. I just want to add 2 solid references on this regard (Both are from an astrophysics context).

The paper by Hogg et al provides a pretty hands-on approach while the the survey of Sharma is more of a survey of MCMC analysis usage in astrophysics.

I am not from the astrophysics community, but I learned a lot with Bayesian inference with MCMC from these two. Hope this can be helpful.

  • $\begingroup$ thanks for your links, it seems to be very instructive. I have to do deeper in the reading. I will contact you as soon as possible if I have issues of understanding about these documentations. $\endgroup$
    – user14831
    Jul 22, 2020 at 14:09
  • $\begingroup$ I read the abstract of Hogg's et al paper and I don't fully understand the sentence in bold : "Above all we emphasize the importance of having a “generative model” for the data, even an approximate one. Once there is a generative model, the subsequent fittingis non-arbitrary because the model permits direct computation of thelikelihood of the parameters or the posterior probability distribution.Construction of a posterior probability distribution is indispensible if there are “nuisance parameters” to marginalize away." . $\endgroup$
    – user14831
    Jul 23, 2020 at 0:05
  • $\begingroup$ I know what is a nuisance parameter, that is, by definition, which brings "noise" to the estimation of parameters. But "Marginalize them away" : what does it mean technically in my case of MCMC method ? How can I get rid of them when I apply the sampler (like Mettropolis-Hastings) ? $\endgroup$
    – user14831
    Jul 23, 2020 at 0:06
  • $\begingroup$ Or maybe "Marginalize away" means that we include them in the estimations of the wanted parameters, that is, actually we don't estimate these nuisance parameters but if we don't estimate them, then we fix them, and so how to find the fixed (fiducial) values of these nuisance parameters ? Hoping you will understand my issue about this notion of marginalization. Thanks $\endgroup$
    – user14831
    Jul 23, 2020 at 0:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.