Parameter Estimation and Model Comparison

The posterior probability distributions of section 5.3.4 constitute a complete description of the author's beliefs about the models and their parameters, once the data (chapter 5) are taken into account. However, it is a very unwieldy description, consisting as it does of two scalar fields defined on a $47$-dimensional space. A summary is required for presentation. The author believes that the most useful summary will be to give, for each parameter $Q_m$, estimates of its posterior marginal expectations

$$<Q_m\vert D, M_N> = \int_{\textrm{All }\mathbf{Q}\textrm{ space}}Q_mP(\mathbf{Q}\vert D, M_N)\mathrm{d}^{47}\mathbf{Q}$$ (5.50)

and

$$<Q_m\vert D, M_M> = \int_{\textrm{All }\mathbf{Q}\textrm{ s... ..._mP(\mathbf{Q}\vert D, M_M)\mathrm{d}^{47}\mathbf{Q}\textrm{,}$$ (5.51)

and of its posterior marginal standard deviations,

$$\sigma{}(Q_m\vert D, M_N) = \sqrt{<Q_m^2\vert D, M_N>-<Q_m\vert D, M_N>^2}$$ (5.52)

and

$$\sigma{}(Q_m\vert D, M_M) = \sqrt{<Q_m^2\vert D, M_M>-<Q_m\vert D, M_M>^2}\textrm{,}$$ (5.53)

where

$$<Q_m^2\vert D, M_n> = \int_{\textrm{All }\mathbf{Q}\textrm{ ... ...^2P(\mathbf{Q}\vert D, M_n)\mathrm{d}^{47}\mathbf{Q}\textrm{.}$$ (5.54)

The summary should also include the estimates of the model posterior probabilities $P(M_n\vert D)$ (equations 5.47, 5.48.)

If a means can be found of drawing samples $\mathbf{Q}^{(\textrm{pos}, n, i, j)}$ from the distribution $P(\mathbf{Q}\vert D, M_n)$, where $i$ runs from $1$ to $I$, $j$ runs from $1$ to $J$, and the total number of samples is, therefore, $IJ$, $<Q_m\vert D, M_n>$ can [48] be estimated by

$$\bar{Q_m}^{(\textrm{pos}, n)} = \frac{\sum_{i=1}^I\bar{Q_m}^{(\textrm{pos}, n, i)}}{I}\textrm{,}$$ (5.55)

where

$$\bar{Q_m}^{(\textrm{pos}, n, i)} = \frac{\sum_{j=1}^JQ_m^{(\textrm{pos}, n, i, j)}}{J}\textrm{.}$$ (5.56)

Similarly, $\sigma{}(Q_m\vert D, M_n)$ can [48] be estimated by

$$\sigma_m^{(\textrm{pos}, n)} =\sqrt{\frac{IJ(\bar{Q_m^2}^{(\... ..., i)}-(\bar{Q_m}^{(\textrm{pos}, n, i)})^2)}{IJ-1}}\textrm{,}$$ (5.57)

where

$$\bar{Q_m^2}^{(\textrm{pos}, n)} = \frac{\sum_{i=1}^I\bar{Q_m^2}^{(\textrm{pos}, n, i)}}{I}\textrm{,}$$ (5.58)

and

$$\bar{Q_m^2}^{(\textrm{pos}, n, i)} = \frac{\sum_{j=1}^J((Q_m^{(\textrm{pos}, n, i, j)})^2)}{J}\textrm{.}$$ (5.59)

Similarly, if a means can be found of drawing samples $\mathbf{Q}^{(\textrm{pri}, n, i, j)}$ from the distribution $P(\mathbf{Q}\vert M_n)$, where $i$ runs from $1$ to $I$, $j$ runs from $1$ to $J$, and the total number of samples is, therefore, $IJ$, $P(D\vert M_n)$ can [48] be estimated by

$$\bar{L}^{(\textrm{pri}, n)} = \frac{\sum_{i=1}^I\bar{L}^{(\textrm{pri}, n, i)}}{I}\textrm{,}$$ (5.60)

where

$$\bar{L}^{(\textrm{pri}, n, i)} = \frac{\sum_{j = 1}^jP(D\vert\mathbf{Q}^{(\textrm{pri}, n, i, j)})}{J}\textrm{.}$$ (5.61)

All of these estimators are [90] non-Bayesian, i.e. they coincide exactly with Bayesian posterior expectations of the quantities being estimated, given the samples, only for specific, but not specified, prior probability distributions over moments of the distribution $P(\mathbf{Q}\vert D, M_n)$. The estimators will, therefore, differ by an amount $\delta$ from the posterior expectations that would be obtained with explicit statements of plausible prior probability distributions over these moments, if this were technically feasible. Fortunately, however, the size of $\delta$ tends [10] to decrease with decreasing marginal likelihood, i.e. with increasing number of samples. Equally fortunately, when the inference process aims to estimate moments of a probability distribution, from samples from that distribution, a critical $\delta$ for each moment is made available, by the degree of variation implied by the higher moments, such that significantly smaller values of $\delta$ can be safely ignored. All this means that, for large numbers of samples, the non-Bayesian estimators are acceptable.

The means of obtaining samples is the leapfrog method. This is explained in detail in Information Theory, Inference and Learning Algorithms [90], where it is attributed to Skilling. Each iteration of the method consists of the application of a leapfrog proposal density, followed by the application of a Metropolis decision algorithm. For sampling from the posterior probability distribution, after the $i$th iteration, the method's state is a set of $J$ parameter vectors $\mathbf{Q}^{(\textrm{pos}, n, i, j)}$, which are the samples to be used in the estimators. The leapfrog proposal density sets

$$\mathbf{Q}'^{(\textrm{pos}, n, i, j)} = 2\mathbf{Q}^{(\textrm{pos}, n, i, k)}-\mathbf{Q}^{(\textrm{pos}, n, i, j)}\textrm{,}$$ (5.62)

where, for each $j$, a $k$ is chosen at random from the other $J-1$ integers from $1$ to $J$, with each of these being equi-probable. The Metropolis decision algorithm then works as follows:

• if $\frac{P(\mathbf{Q}'^{(\textrm{pos}, n, i, j)}\vert D, M_n)}{P(\mathbf{Q}^{(\textrm{pos}, n, i, j)}\vert D, M_n)} \geq{} 1$, then $\mathbf{Q}^{(\textrm{pos}, n, i+1, j)}$ is set to $\mathbf{Q}'^{(\textrm{pos}, n, i, j)}$. Otherwise,
• $\mathbf{Q}^{(\textrm{pos}, n, i+1, j)}$ is set at random, either to $\mathbf{Q}'^{(\textrm{pos}, n, i, j)}$, with probability $\frac{P(\mathbf{Q}'^{(\textrm{pos}, n, i, j)}\vert D, M_n)}{P(\mathbf{Q}^{(\textrm{pos}, n, i, j)}\vert D, M_n)}$, or to $\mathbf{Q}^{(\textrm{pos}, n, i, j)}$, with probability $1-\frac{P(\mathbf{Q}'^{(\textrm{pos}, n, i, j)}\vert D, M_n)}{P(\mathbf{Q}^{(\textrm{pos}, n, i, j)}\vert D, M_n)}$.

An analogous process is used to draw samples $\mathbf{Q}^{(\textrm{pri}, n, i, j)}$ from the prior probability distribution.

As part of the preparation of the present thesis, the author has written Perl code that implements the leapfrog method; in the transparent copy of this thesis, the code for the null model $M_N$ is in the file null, and the code for the main model $M_M$ is in the file main. As far as the author has been able to discover, the present thesis is the first use of the leapfrog method to infer parameters and model likelihoods from real experimental data. Some further details are needed to define the particular version of the leapfrog method used in this thesis:

• Initial vector sets $\mathbf{Q}^{(\textrm{pri}, n, 0, j)}$ and $\mathbf{Q}^{(\textrm{pos}, n, 0, j)}$ are chosen, centred around the analytically-derived prior expectation values $<Q_m\vert M_n>$ of the parameters. The displacements of the $J$ initial vectors from the prior expectation are picked randomly from a perfect top-hat distribution out to a maximum size of each displacement component, which is the product of a constant factor $\alpha_{pri}$ or $\alpha_{pos}$, and the prior standard deviation $\sigma{}(Q_m\vert M_n)$ of the relevant parameter. That is to say, $Q_m^{(\textrm{pri}, n, 0, j)}$ is drawn from a distribution
  
  $$T_0(Q_m^{(\textrm{pri}, n, 0, j)}; <Q_m\vert M_n>-\alpha_{pr... ..., <Q_m\vert M_n>+\alpha_{pri}\sigma{}(Q_m\vert M_n))\textrm{,}$$
  
  
  and $Q_m^{(\textrm{pos}, n, 0, j)}$ is drawn from a distribution
  
  $$T_0(Q_m^{(\textrm{pos}, n, 0, j)}; <Q_m\vert M_n>-\alpha_{po... ..., <Q_m\vert M_n>+\alpha_{pos}\sigma{}(Q_m\vert M_n))\textrm{.}$$
  
  
  It has [90] been shown that the leapfrog method rapidly grows an initial vector set that is narrower than the distribution, from which one wishes to sample. However, the author does not know of any analogous proof that the algorithm can shrink an initial vector set that is broader than the distribution, from which one wishes to sample. Therefore, the initial vector set, for sampling from the prior probability distribution is made a factor of $\exp{}(1)$, in each dimension, narrower than the prior probability distribution. An attempt is made to achieve the same for the posterior probability distribution, using a guess at the width of the posterior probability distribution. This gives $\alpha_{pri} = \frac{\sqrt{3}}{\exp{}(1)}$, and $\alpha_{pos} = \frac{\sqrt{3}}{(\vert D\vert-A)\exp{}(1)}$, where $\vert D\vert = 1018$ is the number of data points, and the number of adjustable parameters is $A = 39$ for the null model $M_N$, and $A = 46$ for the main model $M_M$, because a parameter whose prior probability distribution is a Dirac delta function is not truly adjustable. One further complication is that it is necessary to ensure that, in each dimension, the width of the initial vector set is comfortably greater than the quantum of floating point numbers of the order of magnitude of $<Q_m\vert M_n>$. Mantissae have 53 bits `on typical hardware' [98]. One bit is assumed to be taken up by the sign, and four bits are left, to provide 16 possible starting ordinates, i.e. if the half-width of the top-hat distribution is less than $2^{-48}\vert<Q_m\vert M_n>\vert$, it is adjusted to $2^{-48}\vert<Q_m\vert M_n>\vert$.
• The author believes that the number $J$ of state vectors maintained simultaneously, in the Monte Carlo algorithm, must be at least the number $A$ of adjustable parameters. This is because every state vector ever obtained will be some linear combination of the $J$ initial state vectors, and to obtain a representative sample, the $J$ initial state vectors must, therefore, span the parameter space. In this thesis, therefore, $J = 2A$ is chosen, rather than the $J = 6$ or $J = 12$ suggested in Information Theory, Inference and Learning Algorithms [90].
• The author also believes that having a large value of $J$ will ameliorate the problem, noted by MacKay [90], of the state vectors in the leapfrog algorithm becoming stuck permanently at one vector set, after it has converged, because, as long as the parameter of interest varies smoothly in $\mathbf{Q}$ space, as the $Q_m$ must, if $J \gg{} 1$, there will be a reasonable sampling of the distribution even if all the individual state vectors are stationary. $P(D\vert\mathbf{Q}, M_n)$, however, might not vary smoothly is $\mathbf{Q}$ space. Therefore, the quantitative estimates of $P(D\vert M_n)$ obtained herein should be treated with some scepticism.