Generalisations of Fisher Matrices

Abstract

Fisher matrices play an important role in experimental design and in data analysis. Their primary role is to make predictions for the inference of model parameters—both their errors and covariances. In this short review, I outline a number of extensions to the simple Fisher matrix formalism, covering a number of recent developments in the field. These are: (a) situations where the data (in the form of ($x,y$) pairs) have errors in both x and y; (b) modifications to parameter inference in the presence of systematic errors, or through fixing the values of some model parameters; (c) Derivative Approximation for LIkelihoods (DALI) - higher-order expansions of the likelihood surface, going beyond the Gaussian shape approximation; (d) extensions of the Fisher-like formalism, to treat model selection problems with Bayesian evidence.

1. Introduction

Fisher information matrices are widely used for making predictions for the errors and covariances of parameter estimates. They characterise the expected shape of the likelihood surface in parameter space, subject to an assumption that the likelihood surface is a multivariate Gaussian when viewed as a function of the model parameters. Diagonal terms are the inverse variances of the parameters, conditional on all others being known, and non-zero off-diagonal terms indicate correlations between inferred parameters. Diagonal terms of the inverse Fisher matrix yield the variances of parameters when all others are marginalised over. The Cramér–Rao inequality shows that the variances deduced from the Fisher matrix are lower bounds.

Fisher matrices have been extensively used in cosmology, where future experiments have been designed in order to deduce as precisely as possible the parameters of the standard cosmological model, so-called ΛCDM (Cold Dark Matter, with a cosmological constant Λ), and are routinely used to give “figures-of-merit” [1] for the power of each experiment. Normally, these studies are standard applications of Fisher matrix theory, often simplified by an approximation (which is very good for observations of the Early Universe) that the data are Gaussian-distributed.

In this article, I review a number of generalisations of the Fisher matrix approach. In Section 2 the derivation of the Fisher matrix for Gaussian data is sketched out; in Section 3 we consider Fisher matrices for data pairs that have errors in both x and y; in Section 4 we show how Fisher matrices may be used to estimate biases when some parameters are fixed at incorrect values; in Section 5 we explore better approximations for the likelihood surface (“DALI”), from expansions to higher order in derivatives, and in Section 6 we generalise the use of Gaussian likelihood surfaces to model selection and Bayesian evidence.

2. Gaussian Fields

In cosmology, one is very often dealing with Gaussian random fields, which are characterised statistically entirely by their mean and covariance. A pedagogical derivation for the Fisher matrix when the data $\overrightarrow{y}$ are Gaussian appears in [2]. The negative log-likelihood $\mathcal{L}\equiv -\ln L$ is

$$2\mathcal{L}\left(\overrightarrow{\theta }\right)=lndet2\pi \mathsf{C}+{(\overrightarrow{y}-\overrightarrow{\mu })}^T{\mathsf{C}}^{-1}(\overrightarrow{y}-\overrightarrow{\mu }),$$

(1)

where in general both the mean vector $\overrightarrow{\mu }$ and the covariance matrix $\mathsf{C}=\langle (\overrightarrow{y}-\overrightarrow{\mu }){(\overrightarrow{y}-\overrightarrow{\mu })}^T\rangle$ depend on the model parameters $\overrightarrow{\theta }$. If $\overrightarrow{y}$ represents 1-point statistics, such as Fourier coefficients, then typically $\overrightarrow{\mu }=\overrightarrow{0}$, and all the parameter dependence is in $\mathsf{C}$. If $\overrightarrow{y}$ represents 2-point statistics, then for Gaussian fields they have only approximately a Gaussian distribution, and the analysis is only approximately correct. In this case, the covariance matrix has some parameter dependence through the 4-point function, which for Gaussian fields can be written as products of the 2-point function.

The Gaussian assumption is widely applicable in cosmology, since the quantum fluctuations that are thought to give rise to the density and radiation fields should ensure this, and limits on departures from gaussianity are very tight [3]. Defining the data matrix $\mathsf{D}\equiv (\overrightarrow{y}-\overrightarrow{\mu }){(\overrightarrow{y}-\overrightarrow{\mu })}^T$ and using the matrix identity for positive definite square matrices $lndet\mathsf{C}=\mathrm{Tr}ln\mathsf{C}$, where Tr indicates trace, we can re-write Equation (1) as

$$2\mathcal{L}=\mathrm{Tr}\left[ln\mathsf{C}+{\mathsf{C}}^{-1}\mathsf{D}\right].$$

(2)

Using standard comma notation for partial derivatives, $Z_{,\alpha }=\partial Z/\partial {\theta }_{\alpha }$, and using the matrix identities $\left({\mathsf{C}}^{-1}\right){,}_{\alpha }=-{\mathsf{C}}^{-1}\mathsf{C}{,}_{\alpha }{\mathsf{C}}^{-1}$ and $(ln\mathsf{C}){,}_{\alpha }={\mathsf{C}}^{-1}\mathsf{C}{,}_{\alpha }$, we find after taking two derivatives and then the expectation value,

$${\mathsf{F}}_{\alpha \beta }=\langle \mathcal{L}{,}_{\alpha \beta }\rangle =\frac{1}{2}\mathrm{Tr}[{\mathsf{C}}^{-1}{\mathsf{C}}_{,\alpha }{\mathsf{C}}^{-1}{\mathsf{C}}_{,\beta }+{\mathsf{C}}^{-1}(\overrightarrow{\mu }{,}_{\alpha }\overrightarrow{\mu }{,}_{\beta }^T+\overrightarrow{\mu }{,}_{\beta }\overrightarrow{\mu }{,}_{\alpha }^T)],$$

(3)

The great advantage of the Fisher matrix approach is seen in this example: no data (real or simulated) are required to compute the expected log-likelihood surface, only the statistical properties of the data. This can be a big advantage if simulation is computationally expensive.

3. Fisher Matrix with Errors in x as Well as y

The previous section gives the standard analysis where only the covariance of the y values is considered. Let us now consider the fairly general case where the data consists of data pairs $(X,Y)$, where we have errors in both X and Y. We can compute the Fisher matrix via the application of a Bayesian hierarchical model, provided that the errors in X are small (this will be defined later). The full analysis is given in [4].

We assume $\overrightarrow{X}$ and $\overrightarrow{Y}$ are length m and n vectors (for data pairs, $m=n$, but in fact the analysis is more general), and have Gaussian errors, around true values $\overrightarrow{x}$, $\overrightarrow{y}$, with a covariance matrix $\mathsf{C}$, which also allows correlations between $\overrightarrow{X}$ and $\overrightarrow{Y}$. $\overrightarrow{x}$ and $\overrightarrow{y}$ are not observed, being latent variables, and are essentially nuisance parameters. In fact the $\overrightarrow{y}$ are not independent nuisance parameters as they are assumed to be related precisely to $\overrightarrow{x}$ through a deterministic theoretical model $\overrightarrow{y}=\overrightarrow{\mu }\left(\overrightarrow{x}\right)$ (however, a stochastic element could easily be included). Given the observed data, $\overrightarrow{X},\overrightarrow{Y}$, we seek the posterior $p\left(\overrightarrow{\theta }\right|\overrightarrow{X},\overrightarrow{Y})$. With a uniform prior for $\overrightarrow{\theta }$, this is proportional to the likelihood $L=p(\overrightarrow{X},\overrightarrow{Y}|\overrightarrow{\theta })$. We write this as the marginalised distribution over $\overrightarrow{x}$ and $\overrightarrow{y}$ as

$$\begin{array}{ccc}\hfill L& =& \int p(\overrightarrow{X},\overrightarrow{Y},\overrightarrow{x},\overrightarrow{y}|\overrightarrow{\theta })d\overrightarrow{x}d\overrightarrow{y}=\int p(\overrightarrow{X},\overrightarrow{Y}|\overrightarrow{x},\overrightarrow{y},\overrightarrow{\theta })p(\overrightarrow{x},\overrightarrow{y}|\overrightarrow{\theta })d\overrightarrow{x}d\overrightarrow{y}\hfill \\ \hfill & =& \int p(\overrightarrow{X},\overrightarrow{Y}|\overrightarrow{x},\overrightarrow{y},\overrightarrow{\theta })p\left(\overrightarrow{y}\right|\overrightarrow{x},\overrightarrow{\theta })p\left(\overrightarrow{x}\right|\overrightarrow{\theta })d\overrightarrow{x}d\overrightarrow{y}.\hfill \end{array}$$

(4)

A deterministic $\overrightarrow{y}\left(\overrightarrow{x}\right)$ relation gives a delta function, $p(\overrightarrow{y}|\overrightarrow{x},\overrightarrow{\theta })=\delta (\overrightarrow{y}-\overrightarrow{\mu }(\overrightarrow{x}))$, and assuming a uniform prior for $\overrightarrow{x}$ (a more general prior is considered in [4]), integration over $\overrightarrow{y}$ gives

$$L\propto \int p(\overrightarrow{X},\overrightarrow{Y}|\overrightarrow{x},\overrightarrow{\mu }\left(\overrightarrow{x}\right),\overrightarrow{\theta })d\overrightarrow{x}.$$

(5)

We now assume that the errors in $\overrightarrow{X}$ are small, for which we require that we can truncate at the linear term of the Taylor expansion of $\overrightarrow{\mu }$:

$$\overrightarrow{\mu }\left(\overrightarrow{x}\right)=\overrightarrow{\mu }\left(\overrightarrow{X}\right)+\mathsf{T}\left(\overrightarrow{X}\right)(\overrightarrow{x}-\overrightarrow{X}),$$

(6)

where

$${\mathsf{T}}_{ij}\equiv {\left.\frac{\partial {\mu }_i}{\partial x_j}\right|}_{\overrightarrow{x}=\overrightarrow{X}}.$$

(7)

$\mathsf{T}$ is diagonal for data consisting of $\overrightarrow{X},\overrightarrow{Y}$ pairs.

We assume a multivariate Gaussian for $p(\overrightarrow{X},\overrightarrow{Y}|\overrightarrow{x},\overrightarrow{y})$ (independent of $\overrightarrow{\theta }$), and write the covariance matrix of the data in block form as

$$\begin{array}{ccc}& & \overrightarrow{X}\overrightarrow{Y}\\ \mathsf{C}=& \begin{array}{c}\overrightarrow{X}\\ \overrightarrow{Y}\end{array}& \left(\begin{array}{cc}{\mathsf{C}}_{\mathrm{XX}}& {\mathsf{C}}_{\mathrm{XY}}\\ {\mathsf{C}}_{\mathrm{XY}}^T& {\mathsf{C}}_{\mathrm{YY}}\end{array}\right).\end{array}$$

(8)

Note that ${\mathsf{C}}_{\mathrm{XY}}\equiv \langle (\overrightarrow{X}-\overline{\overrightarrow{X}}){(\overrightarrow{Y}-\overline{\overrightarrow{Y}})}^T\rangle$ is not symmetrical, nor invertible or even square in general; although ${\mathsf{C}}_{\mathrm{XX}}$ and ${\mathsf{C}}_{\mathrm{YY}}$ are. The covariance matrix may include a number of elements, such as intrinsic scatter and measurement noise, with individual covariance matrices adding to give the final $\mathsf{C}$. We also assume that the function $\overrightarrow{\mu }\left(\overrightarrow{x}\right)$ is linear across the width of the Gaussian error distribution of $\overrightarrow{x}$, in which case the likelihood may be integrated analytically, as follows. We write

$$L\propto \int \frac{1}{\sqrt{det\mathsf{C}}}exp\left(-\frac{Q}{2}\right)d\overrightarrow{x}$$

(9)

where $Q\equiv {(\overrightarrow{Z}-\overrightarrow{z})}^T{\mathsf{C}}^{-1}(\overrightarrow{Z}-\overrightarrow{z})$, and $\overrightarrow{z}$ and $\overrightarrow{Z}$ are $m+n$-dimensional vectors: $z_i=x_i$ and $Z_i=X_i$ for $i\le m$, $Z_{m+j}=Y_j$ and $z_{m+j}={\mu }_j(\overrightarrow{X})+{[\mathsf{T}(\overrightarrow{X})(\overrightarrow{x}-\overrightarrow{X})]}_j$. The inverse of $\mathsf{C}$ in block form is

$${\mathsf{C}}^{-1}=\left(\begin{array}{cc}\mathsf{G}& -\mathsf{H}\\ -{\mathsf{H}}^T& \mathsf{E}\end{array}\right)$$

(10)

where

$$\begin{array}{ccc}\hfill \mathsf{G}& =& {\mathsf{C}}_{\mathrm{XX}}^{-1}+{\mathsf{C}}_{\mathrm{XX}}^{-1}{\mathsf{C}}_{\mathrm{XY}}\mathsf{E}{\mathsf{C}}_{\mathrm{XY}}^T{\mathsf{C}}_{\mathrm{XX}}^{-1}\hfill \\ \hfill \mathsf{H}& =& {\mathsf{C}}_{\mathrm{XX}}^{-1}{\mathsf{C}}_{\mathrm{XY}}\mathsf{E}\hfill \\ \hfill \mathsf{E}& =& {({\mathsf{C}}_{\mathrm{YY}}-{\mathsf{C}}_{\mathrm{XY}}^T{\mathsf{C}}_{\mathrm{XX}}^{-1}{\mathsf{C}}_{\mathrm{XY}})}^{-1}.\hfill \end{array}$$

(11)

Defining $\tilde{\overrightarrow{x}}\equiv \overrightarrow{X}-\overrightarrow{x}$, and $\tilde{\overrightarrow{Y}}\equiv \overrightarrow{Y}-\overrightarrow{\mu }(\overrightarrow{X})$, we find that Q has the quadratic form

$$Q={\tilde{\overrightarrow{x}}}^T\mathsf{A}\tilde{\overrightarrow{x}}-{\mathsf{B}}^T\tilde{\overrightarrow{x}}-{\tilde{\overrightarrow{x}}}^T\mathsf{B}+Q^{\prime },$$

(12)

where

$$\begin{array}{ccc}\hfill \mathsf{A}& =& \mathsf{G}+{\mathsf{T}}^T\mathsf{E}\mathsf{T}-\mathsf{H}\mathsf{T}-{\mathsf{T}}^T{\mathsf{H}}^T\hfill \\ \hfill \mathsf{B}& =& (\mathsf{H}-{\mathsf{T}}^T\mathsf{E})\tilde{\overrightarrow{Y}}\equiv \mathsf{P}\tilde{\overrightarrow{Y}}\hfill \\ \hfill Q^{\prime }& =& {\tilde{\overrightarrow{Y}}}^T\mathsf{E}\tilde{\overrightarrow{Y}}.\hfill \end{array}$$

(13)

With the definition of Q in Equation (12), the Gaussian integral of Equation (9) can be performed, using

$$\int exp\left(-\frac{1}{2}{\tilde{\overrightarrow{x}}}^T\mathsf{A}\tilde{\overrightarrow{x}}+{\mathsf{B}}^T\tilde{\overrightarrow{x}}\right)d\tilde{\overrightarrow{x}}=\frac{{\left(2\pi \right)}^{n/2}}{\sqrt{detA}}exp\left(\frac{1}{2}{\mathsf{B}}^T{\mathsf{A}}^{-1}\mathsf{B}\right),$$

(14)

and noting that $Q^{\prime }$ is independent of $\tilde{x}$. The likelihood then simplifies to

$$L\propto \frac{1}{\sqrt{det\mathsf{R}}}exp\left(-\frac{1}{2}{\tilde{\overrightarrow{Y}}}^T{\mathsf{R}}^{-1}\tilde{\overrightarrow{Y}}\right),$$

(15)

where the inverse of the marginal covariance matrix of $\tilde{\overrightarrow{Y}}$ is ${\mathsf{R}}^{-1}=\mathsf{E}-{\mathsf{P}}^T{\mathsf{A}}^{-1}\mathsf{P}$. This is obtained using the Woodbury formula [5] ${(\mathsf{K}+\mathsf{U}\mathsf{W}\mathsf{V})}^{-1}={\mathsf{K}}^{-1}-{\mathsf{K}}^{-1}\mathsf{U}{({\mathsf{W}}^{-1}+\mathsf{V}{\mathsf{K}}^{-1}\mathsf{U})}^{-1}V{\mathsf{K}}^{-1}$, giving

$$\mathsf{R}={\mathsf{C}}_{\mathrm{YY}}-{\mathsf{C}}_{\mathrm{XY}}^T{\mathsf{T}}^T-\mathsf{T}{\mathsf{C}}_{\mathrm{XY}}+\mathsf{T}{\mathsf{C}}_{\mathrm{XX}}{\mathsf{T}}^T,$$

(16)

This is a key result. We see that this looks just like a standard Gaussian (in terms of data) likelihood, but with the covariance matrix $\mathsf{C}$ (${\mathsf{C}}_{\mathrm{YY}}$ in our current notation) replaced by $\mathsf{R}$. Hence to compute the Fisher matrix, we can use the standard formula found in Equation (3) and Equation (15) of [2], and simply replace $\mathsf{C}$ by $\mathsf{R}$:

$${\mathsf{F}}_{\alpha \beta }=\frac{1}{2}\mathrm{Tr}\left[{\mathsf{R}}^{-1}{\mathsf{R}}_{,\alpha }{\mathsf{R}}^{-1}{\mathsf{R}}_{,\beta }+{\mathsf{R}}^{-1}({\overrightarrow{\mu }}_{,\alpha }{\overrightarrow{\mu }}_{,\beta }^T+{\overrightarrow{\mu }}_{,\beta }{\overrightarrow{\mu }}_{,\alpha }^T)\right].$$

(17)

Note that $\mathsf{R}$ depends not only on the standard covariance, but also on the covariance in the independent variable, ${\mathsf{C}}_{\mathrm{XX}}$, the meta-covariance, ${\mathsf{C}}_{\mathrm{XY}}$, and the first partial derivatives of the model function $\overrightarrow{\mu }$. In the case of uncorrelated data pairs, the result reduces to that found in [6]. For the simple case of no correlations between $\overrightarrow{X}$ and $\overrightarrow{Y}$ values $\mathsf{R}={\mathsf{C}}_{\mathrm{YY}}+{\mathsf{T}}^T{\mathsf{C}}_{\mathrm{XX}}\mathsf{T}$, and with diagonal covariance matrices ${\mathsf{C}}_{\mathrm{YY}}$ and ${\mathsf{C}}_{\mathrm{XX}}$ we recover the propagation of error result that the variance of $f\equiv Y-\mu \left(X\right)$ for each data point is effectively

$${\sigma }_f^2={\sigma }_{\mathrm{Y}}^2+{\mu }^{\prime }{\left(X\right)}^2{\sigma }_{\mathrm{X}}^2,$$

(18)

where ${\mu }^{\prime }=\partial \mu /\partial x$ and $\mathsf{C}$ can be replaced in the standard Fisher expression, Equation (3), by a diagonal $n\times n$ matrix with these enhanced entries.

Generalising Still Further

The analysis above is applicable not just to the simple case of data with errors in x as well as y, but to any system where the “data” $\overrightarrow{y}$ depend (in a locally linear way) on any parameters $\overrightarrow{x}$ that have some error associated with them.

4. Systematic Errors, or Errors from Simplified Nested Models

The Fisher matrix can also be useful to determine the errors in parameter inference that arise if one parameter is fixed at an erroneous value. This could arise in a number of contexts, such as a nuisance parameter (e.g., a calibration setting) being fixed at an incorrect value, or when considering nested models. An example of the latter would be cosmological models where the Universe is assumed to be flat. This is an example of a nested model, being a subset of a more general model, but with the curvature parameter (usually given the symbol ${\Omega }_k$) set to zero. In these cases, the maximum likelihood values of all the other parameters are, in general, shifted from their maximum likelihood values in the more general model. See Figure 1 for an illustration of this in two dimensions. With the usual Fisher assumption that the likelihood surface is a Gaussian function of the parameters, these shifts can be computed using the Fisher matrix.

We consider two models, M, which has more ($n^{\prime }+p$) parameters than a simpler nested model $M^{\prime }$, which has $n^{\prime }$. The extra parameters are designated ${\psi }_{\zeta }$, and these are fixed in $M^{\prime }$ at values that are $\delta {\psi }_{\zeta }$ from their maximum likelihood values in M. In this case, the maximum likelihood values of all other parameters of $M^{\prime }$, $\overrightarrow{\theta }$, are systematically shifted by [7,8]

$$\delta {\theta }_{\alpha }=-{\left({\mathsf{F}}^{\prime -1}\right)}_{\alpha \beta }{\mathsf{G}}_{\beta \zeta }\delta {\psi }_{\zeta }\alpha ,\beta =1\dots n^{\prime },\zeta =1\dots p$$

(19)

where

$${\mathsf{G}}_{\beta \zeta }=\frac{1}{2}\mathrm{Tr}\left[{\mathsf{C}}^{-1}{\mathsf{C}}_{,\beta }{\mathsf{C}}^{-1}{\mathsf{C}}_{,\zeta }+{\mathsf{C}}^{-1}({\overrightarrow{\mu }}_{,\zeta }{\overrightarrow{\mu }}_{,\beta }^T+{\overrightarrow{\mu }}_{,\beta }{\overrightarrow{\mu }}_{,\zeta }^T)\right],$$

(20)

which we recognise as a subset of the Fisher matrix.

5. Beyond the Gaussian Approximation—DALI

The Fisher matrix approach assumes that the likelihood surface is a multivariate Gaussian, which will be asymptotically true near the peak, but may not be a good approximation over the range of parameter values of interest. A generalisation of the Fisher matrix is DALI, Derivative Approximation for LIkelihoods [9], which expands the likelihood surface to include higher-order derivatives than the second. This is a rather elegant expansion, in derivatives rather than parameters, that ensures that the approximate distribution is a genuine probability distribution—i.e., it is non-negative and normalisable, non-divergent and asymptotically approaches the true likelihood.

The starting point is a Taylor expansion of the likelihood:

$$\begin{array}{cc}\hfill -\mathcal{L}\equiv lnL\approx & ln{L}_0+\frac{1}{2!}{\mathsf{F}}_{\alpha \beta }\mathsf{\Delta }{\theta }_{\alpha }\mathsf{\Delta }{\theta }_{\beta }+\frac{1}{3!}{\mathsf{S}}_{\alpha \beta \gamma }\mathsf{\Delta }{\theta }_{\alpha }\mathsf{\Delta }{\theta }_{\beta }\mathsf{\Delta }{\theta }_{\gamma }\hfill \\ & +\frac{1}{4!}{\mathsf{Q}}_{\alpha \beta \gamma \delta }\mathsf{\Delta }{\theta }_{\alpha }\mathsf{\Delta }{\theta }_{\beta }\mathsf{\Delta }{\theta }_{\gamma }\mathsf{\Delta }{\theta }_{\delta }+\dots ,\hfill \end{array}$$

(21)

where $L_0$ is a normalization constant and ${\mathsf{F}}_{\alpha \beta }={\mathcal{L}}_{,\alpha \beta }$, ${\mathsf{S}}_{\alpha \beta \gamma }={\mathcal{L}}_{,\alpha \beta \gamma }$ and ${\mathsf{Q}}_{\alpha \beta \gamma \delta }={\mathcal{L}}_{,\alpha \beta \gamma \delta }$.

If the expansion is arranged in order of derivatives, the expressions are normalisable and positive-definite. For example, to second order in the $\overrightarrow{\mu }$ derivatives, and assuming $\mathsf{C}$ is independent of $\overrightarrow{\theta }$, we have

$$\begin{array}{cc}\hfill L=L_{0}exp[& -\frac{1}{2}{\overrightarrow{\mu }}_{,\alpha }^T{\mathsf{C}}^{-1}{\overrightarrow{\mu }}_{,\beta }\mathsf{\Delta }{\theta }_{\alpha }\mathsf{\Delta }{\theta }_{\beta }-(\frac{1}{2}{\overrightarrow{\mu }}_{,\alpha \beta }^T{\mathsf{C}}^{-1}{\overrightarrow{\mu }}_{,\gamma }\mathsf{\Delta }{\theta }_{\alpha }\mathsf{\Delta }{\theta }_{\beta }\mathsf{\Delta }{\theta }_{\gamma }\hfill \\ & +\frac{1}{8}{\overrightarrow{\mu }}_{,\delta \gamma }^T{\mathsf{C}}^{-1}{\overrightarrow{\mu }}_{,\beta \alpha }\mathsf{\Delta }{\theta }_{\alpha }\mathsf{\Delta }{\theta }_{\beta }\mathsf{\Delta }{\theta }_{\gamma }\mathsf{\Delta }{\theta }_{\delta })+\mathcal{O}\left(3\right)].\hfill \end{array}$$

(22)

This is apparently true at every order (see [9] for the third-order expansion, and [10] for the case where the parameter dependence is in $\mathsf{C}$). Figure 2 shows the improvement in the expected likelihood surfaces for a supernova cosmology experiment.

6. The Expected Bayesian Evidence—Generalising Fisher Matrices to Model Selection

At the root of the Fisher matrix formalism is the Laplace approximation, i.e., the assumption that the likelihood surface is a multivariate Gaussian when viewed as a function of the model parameters. We can generalise this to the higher-level question of model selection, where we compute the posterior probabilities of different models, given the data collected, but regardless of the model parameters $\overrightarrow{\theta }$. The ratio of these probabilities is the ratio of the prior model probabilities, multiplied by the “Bayes factor”, which is the ratio of the marginal likelihoods (or Bayesian evidence) of the models, where the evidence for a model M is

$$p\left(\overrightarrow{y}\right|M)=\int p\left(\overrightarrow{y}\right|\overrightarrow{\theta },M)p\left(\overrightarrow{\theta }\right|M)d\overrightarrow{\theta }.$$

(23)

With the Laplace approximation for the first, likelihood term, and a uniform prior (which can be generalised to a Gaussian prior), we can compute the expected evidence (conditional on some fiducial set of parameters) by performing Gaussian integrals. For nested models (with $n^{\prime }$ and $n=n^{\prime }+p$ parameters respectively), the considerations of Section 4 on the locations of the peak likelihood is relevant, and the result depends on the shifts of the fiducial parameters away from the values that are fixed in the lower-dimensional model, $\delta {\psi }_{\zeta }$. If we further approximate that the expected Bayes factor is the ratio of the expected evidences, then the expected Bayes factor is (see [7] for details)

$$\langle B\rangle ={\left(2\pi \right)}^{-p/2}\frac{\sqrt{det\mathsf{F}}}{\sqrt{det{\mathsf{F}}^{\prime }}}exp\left(-\frac{1}{2}\delta {\theta }_{\alpha }{\mathsf{F}}_{\alpha \beta }\delta {\theta }_{\beta }\right)\prod _{q=1}^p\mathsf{\Delta }{\theta }_{n^{\prime }+q}$$

(24)

where $\mathsf{\Delta }{\theta }_{\alpha }$ are the prior ranges of the additional p parameters in the extended model, and the offsets $\delta {\theta }_{\alpha }$ are given by Equation (19). Note that $\mathsf{F}$ is an $n\times n$ matrix, ${\mathsf{F}}^{\prime }$ is $n^{\prime }\times n^{\prime }$, and $\mathsf{G}$ is an $n^{\prime }\times p$ block of the full $n\times n$ Fisher matrix $\mathsf{F}$, given by Equation (20). The expression we find is a specific example of the Savage-Dickey density ratio [11]; here we explicitly use the Laplace approximation to compute the offsets in the parameter estimates which accompany the wrong choice of model.

Figure 3 shows the ratio of expected evidences, assuming the Laplace approximation (as the Fisher matrix does), for nested cosmological models. Details are in the caption, but essentially one parameter is fixed in the simpler model, but allowed to vary in the more complex model. If the more complex model applies, then the data will favour the simpler model if the parameter is close to the fixed value. This is shown in the figure by the cusp in the graph. $lnB$ is positive to the left of the cusp, and negative to the right.

7. Discussion

This article reviews some recent developments in Fisher matrix theory, which have been motivated by cosmology. The Fisher matrix for data consisting of pairs that have errors in both x and y is derived, as a specific example of a general result where the data can depend on arbitrary variables x that may be measured with some error. The Fisher matrix is shown to be able to determine biases in some parameters when others are set to fixed values (such as in nested models where the simpler model does not allow some parameters to vary). DALI, which goes beyond the Laplace approximation by using higher-order derivatives, is found to allow much more accurate predictions for the expected shape of the likelihood surface. Finally, the concept of expected probabilities in the Laplace approximation is generalised to model selection, by computing the expected Bayesian evidence.