Computational issues: normalization, fit summaries, priors and approximations

At this point it is important to understand that in Bayesian approach the full result of the inference is given by final distribution, that in our case is - we rewrite it here:

$$f(m,c,\sigma_v\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I) = k\, \prod_i \frac{1}{\sqrt{\sigma^2_v + \sigma_{y_i}^2+m^2\,\sigm... ...+ \sigma_{y_i}^2+m^2\,\sigma_{x_i}^2) } \right]}\, f(m,c,\sigma_v\,\vert\,I)\,,$$

where $k$ is `simply' a normalization factor. (This factor is usually the most difficult thing to calculate and it is often obtained approximately by numerical methods. But this is, in principle, just a technical issue.) Once we have got $k$ we have a full knowledge about $f(m,c,\sigma_v\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I)$ and therefore about our uncertainty concerning the model parameters, the distribution of each of which can be obtained by marginalization:

$$f(m\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I) = \int\! f(m,c,\sigma_v\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I)\,\,dc\,d\sigma_v$$ (52)

$$f(c\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I) = \int\! f(m,c,\sigma_v\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I)\,\,dm\,d\sigma_v$$ (53)

$$f(\sigma_v\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I) = \int\! f(m,c,\sigma_v\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I)\,\,dm\,dc\,.$$ (54)

Similarly the joint distribution of $m$ and $c$ can be obtained as

$$f(m,c\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I) = \int\! f(m,c,\sigma_v\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I)\,d\sigma_v\,,$$ (55)

from which we can easily see that we recover Eq. (30) in the case we think the extra variability discussed in the previous section is absent. This limit case corresponds to a prior of $\sigma _v$ sharply peaked around zero, i.e. $f(\sigma_v\,\vert\,I) = \delta(\sigma_v)$.

Other interesting limit cases are the following.

• Errors only on the $y$ axis and no extra variability.
  Making the limit of Eq. (30) for $\sigma_{x_i}\rightarrow 0$ and neglecting irrelevant factors we get
  

  $$f(m,c\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I) \propto \prod_i \exp{ \left[ -\frac{(y_i-m\,x_i-c)^2} {2\, \sigma_{y_i}^2} \right]}\, f(m,c\,\vert\,I)$$ (56)

  $$\propto \exp{ \left[-\frac{1}{2}\sum_i \frac{(y_i-m\,x_i-c)^2} {\sigma_{y_i}^2} \right]}\, f(m,c\,\vert\,I) \,.$$ (57)

  
  
  This is the best known and best understood case.
• Errors only on the $y$ axis and extra variability.
  Making the limit of Eq. (53) for $\sigma_{x_i}\rightarrow 0$
  

  $$f(m,c,\sigma_v\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I) \propto \prod_i \frac{1}{\sqrt{\sigma^2_v + \sigma_{y_i}^2}}\, \exp{ \lef... ...2} {2\, (\sigma^2_v + \sigma_{y_i}^2) } \right]}\, f(m,c,\sigma_v\,\vert\,I)\,.$$ (58)

  
  

• Scattering of data point around the hypothesized straight line only due to `extra variability'.
  

  $$f(m,c,\sigma_v\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I) \propto \sigma_v^{-N} \prod_i \exp{ \left[ -\frac{(y_i-m\,x_i-c)^2} {2\,\sigma^2_v } \right]}\, f(m,c,\sigma_v\,\vert\,I)$$ (59)

  $$\propto \sigma_v^{-N} \exp{ \left[ -\frac{1}{2\,\sigma_v^2} \,\sum_i (y_i-m\,x_i-c)^2 \right]}\, f(m,c,\sigma_v\,\vert\,I) \,.$$ (60)

  
  
  This case corresponds to the joint determination of $m$, $c$ and $\sigma _v$ made by the method of the `residuals', that can be considered a kind of approximated solution of Eq. (61), achieved by iteration. [Indeed, if there are `enough' data points the `best estimates' achieved by the residual method are very close to the expected values of $m$, $c$ and $\sigma _v$ evaluated from $f(m,c,\sigma_v\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I)$ if we assumed a flat prior distribution for the parameters.]

Although, as it has been pointed out above, the full result of the inference is provided by the final pdf, often we do not need such a detailed description of our uncertainty, and we are only interested to provide some `summaries'. The most interesting ones are the expected values, standard deviations and correlation coefficients, i.e. $E(m)$, $E(c)$, $E(\sigma_v)$, $\sigma(m)$, $\sigma(c)$, $\sigma(\sigma_v)$, $\rho(m,c)$, $\rho(m,\sigma_v)$ and $\rho(c,\sigma_v)$. They are evaluated from $f(m,c,\sigma_v)$ using their definitions, that are assumed to be known [hereon we often omit the conditions on which the pdf depends, and we write $f(m,c,\sigma_v)$ instead of $f(m,c,\sigma_v\,\vert\,{\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}},I)$, and so on]. Obviously, these are not the only possible summaries. One might report in addition the mode or the median of each variable, one-dimensional or multi-dimensional probability regions [i.e. regions in the space of the parameters that are believed to contain the true value of the parameter(s) with a well defined probability level], and so on. It all depends on how standard or unusual the shape of $f(m,c,\sigma_v)$ is. I just would like to stress that the most important summaries are expected value, standard deviation and correlation coefficients, because these are the quantities that mostly matter in subsequent evaluations of uncertainty. Giving only `most probable' values and probability intervals might bias the results of further analyzes [15].

The prior $f(m,c,\sigma_v\,\vert\,I)$ has been left on purpose open in the above formulas, although we have already anticipated that usually a flat prior about all parameters gives the correct result in most 'healthy' cases, characterized by a sufficient number of data points. I cannot go here through an extensive discussion about the issue of the priors, often criticized as the weak point of the Bayesian approach and that are in reality one of its points of force. I refer to more extensive discussions available elsewhere (see e.g. [2] and references therein), giving here only a couple of advices. A flat prior is in most times a good starting point (unless one uses some packages, like BUGS [11], that does not like flat prior in the range $-\infty$ to $+\infty$; in this case one can mimic it with a very broad distribution, like a Gaussian with very large $\sigma$). If the result of the inference `does not offend your physics sensitivity', it means that, essentially, flat priors have done a good job and it is not worth fooling around with more sophisticated ones. In the specific case we are looking closer, that of Eq. (53), the most critical quantity to watch is obviously $\sigma _v$, because it is positively defined. If, starting from a flat prior (also allowing negative values), the data constrain the value of $\sigma _v$ in a (positive) region far from zero, and - in practice consequently - its marginal distribution is approximatively Gaussian, it means the flat prior was a reasonable choice. Otherwise, the next-to-simple modeling of $\sigma _v$ is via the step function $\theta(\sigma_v)$. A more technical choice would be a gamma distribution, with suitable parameters to `easily' accommodate all envisaged values of $\sigma _v$.

The easiest case, that happens very often if one has `many' data points (where `many' might be already as few as some dozens), is that $f(m,c,\sigma_v)$ obtained starting from flat priors is approximately a multi-variate Gaussian distribution, i.e. each marginal is approximately Gaussian. In this case the expected value of each variable is close to its mode, that, since the prior was a constant, corresponds to the value for which the likelihood ${\cal L}(m,c,\sigma_v\,;\, {\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}})$ gets its maximum. Therefore the parameter estimates derived by the maximum likelihood principle are very good approximations of the expected values of the parameters calculated directly from $f(m,c,\sigma_v)$. In a certain sense the maximum likelihood principle best estimates are recovered as a special case that holds under particular conditions (many data points and vague priors). If either condition fails, the result the formulas derived from such a principle might be incorrect. This is the reason I dislike unneeded principles of this kind, once we have a more general framework, of which the methods obtained by `principles' are just special cases under well defined conditions.

The simple case in which $f(m,c,\sigma_v)$ is approximately multi-variate Gaussian allows also to approximately evaluate the covariance matrix of the fit parameters from the Hessian of its logarithm.⁶This is due to a well known property of the multi-variate Gaussian and it is not strictly related to flat priors. In fact it can easily proved that if the generic $f({\mbox{\boldmath$\theta$}})$ is a multivariate Gaussian, then

$$(V^{-1})_{ij}({\mbox{\boldmath$\theta$}}) = \left.\frac{\partial^2 \varphi} {\partial\theta_i\,\partial\theta_j} \right\vert _{{\mbox{\boldmath$\theta$}}={\mbox{\boldmath$\theta$}}_m}$$ (61)

where

$$\varphi({\mbox{\boldmath$\theta$}}) = -\log f({\mbox{\boldmath$\theta$}})\,,$$ (62)

$V_{ij}({\mbox{\boldmath$\theta$}})$ is the covariance matrix of the parameters and ${\mbox{\boldmath$\theta$}}_m$ is the value for which $f({\mbox{\boldmath$\theta$}})$ gets its maximum and then $\varphi({\mbox{\boldmath$\theta$}})$ its minimum.

An interesting feature of this approximated procedure is that, since it is based on the logarithm of the pdf, normalization factors are irrelevant. In particular, if the priors are flat, the relevant summaries of the inference can be obtained from the logarithm of the likelihood, stripped of all irrelevant factors (that become additive constants in the logarithm and vanish in the derivatives). Let us write down, for some cases of interest, the minus-log-likelihoods, stripped of constant terms and indicated by $L$, i.e. $\varphi({\mbox{\boldmath$\theta$}}\, ;\, {\mbox{\boldmath$x$}},{\mbox{\boldmath... ...theta$}}\, ;\, {\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}} ) + \mbox{\it const}$.

• Simplest case: linear fit with only known errors on the $y$ axis [from Eq. (58)]:
  

  $$L(m,c\,;\, {\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}}) = \frac{1}{2}\sum_i\frac{(y_i-m\,x_i-c)^2} {\sigma_{y_i}^2} = \frac{1}{2}\, \chi^2(m,c\,;\, {\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}})\,,$$ (63)

  
  
  where we recognize the famous chi-squared. Applying Eq. (63) we get then the covariance matrix of the fit parameters as
  

  $$(V^{-1})_{m,c} = \left.\frac{1}{2} \frac{\partial^2 \chi^2(m,c\,;\, {\mbox{\boldma... ...rtial m\,\partial c} \right\vert _{\begin{array}{l} m=m_m \\ c=c_m\end{array}}$$ (64)

  
  
  (See Ref. [2] for the fully developed example yielding analytic formulas for the expected values and covariance matrix of the $m$ and $c$.) Note that the often used (but also often misused![15]) ` $\Delta\,\chi^2=1$ rule' to calculate the covariance matrix of the parameters comes from the same Gaussian approximation of the final pdf and prior insensitivity. [And, because of the factor $1/2$ between Eqs. (63) and (66), there is an equivalent `$\Delta\,$minus-log-likelihood = $1/2$' rule, applicable under the same conditions].
• Errors also on the $y$ axis:
  

  $$L(m,c\,;\, {\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}}) = \frac{1}{2}\sum_i \log {(\sigma_{y_i}^2+m^2\,\sigma_{x_i}^2)} + \frac{1}{2}\sum_i\frac{(y_i-m\,x_i-c)^2} {\sigma_{y_i}^2+m^2\,\sigma_{x_i}^2}\,.$$ (65)

  
  
  In this case expected values and covariance matrix cannot be obtained directly in closed form. Nevertheless, one can use iteratively the formulas for $\sigma_{x_i}=0$ in which the estimate of $m$ is used to evaluate the terms $\sigma_{y_i}^2+m^2\,\sigma_{x_i}^2$ (having the meaning of effective $y$-error) in the likelihood of the next iteration. Instead it is wrong to simply replace the denominator of the $\chi^2$ of Eq. (65) with $\sigma_{y_i}^2+m^2\,\sigma_{x_i}^2$, because this approximation does not take into account the first term of the r.h.s. of Eq. (67) and the slope $m$ will be underestimated (as a consequence, the intercept $c$ will be over- or under-estimated, depending on the sign of the correlation coefficient between $m$ and $c$, a sign that depends on the sign of the barycenter of the $x$ points.)
• Dispersion on the $y$ axis only due to $\sigma _v$ [from Eq. (61)]:
  

  $$L(m,c,\sigma_v\,;\, {\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}}) = N\,\log\sigma_v + \frac{1}{2\,\sigma_{v}^2}\sum_i\,(y_i-m\,x_i-c)^2\,.$$ (66)

  
  

• The most complete case seen here [from Eq. (53)]:
  

  $$L(m,c,\sigma_v\,;\, {\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}}) = \frac{1}{2}\sum_i \log {(\sigma^2_v + \sigma_{y_i}^2+m^2\,\sigma_... ...i \frac{(y_i-m\,x_i-c)^2} {\sigma^2_v + \sigma_{y_i}^2+m^2\,\sigma_{x_i}^2 }\,.$$ (67)

  
  

• As the previous item, but for the general $\mu_y(\,)$ [from Eq. (34)]:
  

  $$L({\mbox{\boldmath$\theta$}},\sigma_v\,;\, {\mbox{\boldmath$x$}},{\mbox{\boldmath$y$}}) \!\!\! \approx \frac{1}{2}\sum_i \log {[\sigma^2_v + \sigma_{y_i}^2 +{\mu_y^{\,\... ...+{\mu_y^{\,\prime}}^2 (x_i; {\mbox{\boldmath$\theta$}})\cdot\sigma_{x_i}^2 }\,.$$ (68)