Amplitude Marginaliser: Amplitude Marginaliser library

Author:: Matthew Pitkin

Description

This library (amplitude-marginaliser) provides functions to analytically marginalise a Gaussian likelihood ratio (i.e. the ratio of a Gaussian likelihood given a model to a Gaussian likelihood for pure noise) over the amplitudes of an arbitrary number of independent components of a signal model. The algorithm used is described in the Appendix of Pitkin, Williams, Fletcher & Grant, arXiv:1406.1712 (if you use this library for work in a publication I would be grateful if you cite that paper). E.g., for the five amplitudes $A$ , $B$ , $C$ , of the five model components $Af(x_1)$ , $Bg(x_2)$ , $Ch(x_3)$ , $Di(x_4)$ and $Ej(x_5)$ ), the integral performed would be:

$\mathcal{L}(x_1, x_2, x_3, x_4, x_5) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty, 0}^{\infty} \exp{-\sum_{n=1}^N\frac{1}{2\sigma^2} ((Af(x_1)_n + Bg(x_2)_n + Ch(x_3)_n + Di(x_4)_n + Ej(x_5)_n)^2 - 2d_n(Af(x_1)_n + Bg(x_2)_n + Ch(x_3)_n + Di(x_4)_n + Ej(x_5)_n) ) } {\rm d}A {\rm d}B {\rm d}C {\rm d}D {\rm d}E.$

In the above equation we leave it open as to whether the final model amplitude is marginalised between $-\infty$ and $\infty$ , or 0 and $\infty$ (i.e. is only allowed to be positive). When expanded out and terms grouped such that summations on pairs of models, and models on data, are done separately the equation for the first integral e.g. over $A$ the equation can be placed into the format

$\mathcal{L}_A = \int_{-\infty}^{\infty} \exp{\left(-\frac{1}{2\sigma^2}(XA^2 + YA + Z)\right)} {\rm d}A,$

which has the solution (via e.g. Mathematica)

$\mathcal{L}_A = \sqrt{\frac{2\pi\sigma^2}{X}}\exp{\left(-\frac{1}{2\sigma^2}(Z-Y^2/(4X))\right)}.$

For the case that the integral is between 0 and $\infty$ this is instead

$\mathcal{L}_A = \sqrt{\frac{\pi\sigma^2}{2X}}\exp{\left(-\frac{1}{2\sigma^2}(Z-Y^2/(4X))\right)} {\rm erfc}{\left(\frac{Y}{2\sqrt{2\sigma^2X}}\right)}.$

Mathematica can in fact solve the full equation (for five amplitudes at least), presumably by repeated iterations of the above formula, but when the number of models is greater than 3 the output is so long that it starts to become unusable. However, it is possible to iterate the above equation numerically as demonstrated in these library functions.

As this is a likelihood ratio you can get back to the likelihood by multiplying by:

$(2\pi\sigma^2)^{-N/2}\exp{\left(-\sum_{n=1}^N \frac{d_n}{2\sigma^2} \right)}.$

Functions are also provided to analytically marginalise over all component amplitudes bar a final component (which itself can be a composite of multiple components). This requires repeated iterations of the above equation, but stopping for the final terms.

These functions do not calculate a Bayes factor as no model parameter priors are used. The functions integrate the model amplitudes over an infinite range and they don't explicitly include prior values on them. However, if one can assign a flat prior (or wide Gaussian prior) on the amplitudes, and the likelihood falls of to zero well within that prior range, then amplitude priors can be subsequently applied to the output. Priors on any other (non-independent amplitude) model parameters can be applied as required as this likelihood ratio must still be calculated over those parameter ranges. However, if actually comparing Bayesian evidences in which the same models have been used then using this likelihood ratio is fine as the priors would cancel out.

It should also be noted that as in Bretthorst the marginalisation integral could be simplified by converting the independent model components into a new set of orthogonal models, which makes the model component cross terms become zero. In the orthogonal model space the the likelihood model is much simpler, but does require pre-calculating the orthogonal models. Another issue with this is that it requires that all the model amplitudes are integrated between $-\infty$ and $\infty$ .

Example

An example of using these functions might be if you want to calculate the Bayesian evidence for a signal with an unknown amplitude in data which also contains a background unwanted polynomial variation. The signal amplitude and the polynomial coefficients (up to whatever order required) can be analytically marginalised out. This means that the data does not need detrending to remove the polynomial background, which could effect the signal by creating unwanted artifacts and causing a loss of signal power.

An example of the algorithm being used for such a case is given in Pitkin, Williams, Fletcher & Grant, arXiv:1406.1712, with accompanying code found here.