Search

Motivated by big data and the vast parameter spaces in modern machine learning models, optimisation approaches to Bayesian inference have seen a surge in popularity in recent years. In this paper, we address the connection between the popular new methods termed generalised Bayesian inference and Bayes linear methods. We propose a further generalisation to Bayesian inference that unifies these and other recent approaches by considering the Bayesian inference problem as one of finding the closest point in a particular solution space to a data generating process, where these notions differ depending on user-specified geometries and foundational belief systems. Motivated by this framework, we propose a generalisation to Bayes linear approaches that enables fast and principled inferences that obey the coherence requirements implied by domain restrictions on random quantities. We demonstrate the efficacy of generalised Bayes linear inference on a number of examples, including monotonic regression and inference for spatial counts. This paper is accompanied by an R package available at github.com/astfalckl/bayeslinear., Comment: Submitted to the Journal of the Royal Statistical Society: Series B

Computer model calibration involves using partial and imperfect observations of the real world to learn which values of a model's input parameters lead to outputs that are consistent with real-world observations. When calibrating models with high-dimensional output (e.g. a spatial field), it is common to represent the output as a linear combination of a small set of basis vectors. Often, when trying to calibrate to such output, what is important to the credibility of the model is that key emergent physical phenomena are represented, even if not faithfully or in the right place. In these cases, comparison of model output and data in a linear subspace is inappropriate and will usually lead to poor model calibration. To overcome this, we present kernel-based history matching (KHM), generalising the meaning of the technique sufficiently to be able to project model outputs and observations into a higher-dimensional feature space, where patterns can be compared without their location necessarily being fixed. We develop the technical methodology, present an expert-driven kernel selection algorithm, and then apply the techniques to the calibration of boundary layer clouds for the French climate model IPSL-CM., Comment: 50 pages

Whether and how data scientists, statisticians and modellers should be accountable for the AI systems they develop remains a controversial and highly debated topic, especially given the complexity of AI systems and the difficulties in comparing and synthesising competing claims arising from their deployment for data analysis. This paper proposes to address this issue by decreasing the opacity and heightening the accountability of decision making using AI systems, through the explicit acknowledgement of the statistical foundations that underpin their development and the ways in which these dictate how their results should be interpreted and acted upon by users. In turn, this enhances (1) the responsiveness of the models to feedback, (2) the quality and meaning of uncertainty on their outputs and (3) their transparency to evaluation. To exemplify this approach, we extend Posterior Belief Assessment to offer a route to belief ownership from complex and competing AI structures. We argue that this is a significant way to bring ethical considerations into mathematical reasoning, and to implement ethical AI in statistical practice. We demonstrate these ideas within the context of competing models used to advise the UK government on the spread of the Omicron variant of COVID-19 during December 2021., Comment: 26 pages, 2 figures

Modern scientific problems are often multi-disciplinary and require integration of computer models from different disciplines, each with distinct functional complexities, programming environments, and computation times. Linked Gaussian process (LGP) emulation tackles this challenge through a divide-and-conquer strategy that integrates Gaussian process emulators of the individual computer models in a network. However, the required stationarity of the component Gaussian process emulators within the LGP framework limits its applicability in many real-world applications. In this work, we conceptualize a network of computer models as a deep Gaussian process with partial exposure of its hidden layers. We develop a method for inference for these partially exposed deep networks that retains a key strength of the LGP framework, whereby each model can be emulated separately using a DGP and then linked together. We show in both synthetic and empirical examples that our linked deep Gaussian process emulators exhibit significantly better predictive performance than standard LGP emulators in terms of accuracy and uncertainty quantification. They also outperform single DGPs fitted to the network as a whole because they are able to integrate information from the partially exposed hidden layers. Our methods are implemented in an R package $\texttt{dgpsi}$ that is freely available on CRAN.

Any experiment with climate models relies on a potentially large set of spatio-temporal boundary conditions. These can represent both the initial state of the system and/or forcings driving the model output throughout the experiment. Whilst these boundary conditions are typically fixed using available reconstructions in climate modelling studies, they are highly uncertain, that uncertainty is unquantified, and the effect on the output of the experiment can be considerable. We develop efficient quantification of these uncertainties that combines relevant data from multiple models and observations. Starting from the coexchangeability model, we develop a coexchangable process model to capture multiple correlated spatio-temporal fields of variables. We demonstrate that further exchangeability judgements over the parameters within this representation lead to a Bayes linear analogy of a hierarchical model. We use the framework to provide a joint reconstruction of sea-surface temperature and sea-ice concentration boundary conditions at the last glacial maximum (19-23 ka) and use it to force an ensemble of ice-sheet simulations using the FAMOUS-Ice coupled atmosphere and ice-sheet model. We demonstrate that existing boundary conditions typically used in these experiments are implausible given our uncertainties and demonstrate the impact of using more plausible boundary conditions on ice-sheet simulation., Comment: Accepted with major revisions in the Journal of the American Statistical Association, resubmission is in review

Deep Gaussian processes (DGPs) provide a rich class of models that can better represent functions with varying regimes or sharp changes, compared to conventional GPs. In this work, we propose a novel inference method for DGPs for computer model emulation. By stochastically imputing the latent layers, our approach transforms a DGP into a linked GP: a novel emulator developed for systems of linked computer models. This transformation permits an efficient DGP training procedure that only involves optimizations of conventional GPs. In addition, predictions from DGP emulators can be made in a fast and analytically tractable manner by naturally utilizing the closed form predictive means and variances of linked GP emulators. We demonstrate the method in a series of synthetic examples and empirical applications, and show that it is a competitive candidate for DGP surrogate inference, combining efficiency that is comparable to doubly stochastic variational inference and uncertainty quantification that is comparable to the fully-Bayesian approach. A $\texttt{Python}$ package $\texttt{dgpsi}$ implementing the method is also produced and available at https://github.com/mingdeyu/DGP., Comment: This article has been accepted for publication in Technometrics, published by Taylor & Francis

In many real-world applications, we are interested in approximating black-box, costly functions as accurately as possible with the smallest number of function evaluations. A complex computer code is an example of such a function. In this work, a Gaussian process (GP) emulator is used to approximate the output of complex computer code. We consider the problem of extending an initial experiment (set of model runs) sequentially to improve the emulator. A sequential sampling approach based on leave-one-out (LOO) cross-validation is proposed that can be easily extended to a batch mode. This is a desirable property since it saves the user time when parallel computing is available. After fitting a GP to training data points, the expected squared LOO (ES-LOO) error is calculated at each design point. ES-LOO is used as a measure to identify important data points. More precisely, when this quantity is large at a point it means that the quality of prediction depends a great deal on that point and adding more samples nearby could improve the accuracy of the GP. As a result, it is reasonable to select the next sample where ES-LOO is maximised. However, ES-LOO is only known at the experimental design and needs to be estimated at unobserved points. To do this, a second GP is fitted to the ES-LOO errors and where the maximum of the modified expected improvement (EI) criterion occurs is chosen as the next sample. EI is a popular acquisition function in Bayesian optimisation and is used to trade-off between local/global search. However, it has a tendency towards exploitation, meaning that its maximum is close to the (current) "best" sample. To avoid clustering, a modified version of EI, called pseudo expected improvement, is employed which is more explorative than EI yet allows us to discover unexplored regions. Our results show that the proposed sampling method is promising.

Classification is a vital tool that is important for modelling many complex numerical models. A model or system may be such that, for certain areas of input space, the output either does not exist, or is not in a quantifiable form. Here, we present a new method for classification where the model outputs are given distinct classifying labels, which we model using a latent Gaussian process (GP). The latent variable is estimated using MCMC sampling, a unique likelihood and distinct prior specifications. Our classifier is then verified by calculating a misclassification rate across the input space. Comparisons are made with other existing classification methods including logistic regression, which models the probability of being classified into one of two regions. To make classification predictions we draw from an independent Bernoulli distribution, meaning that distance correlation is lost from the independent draws and so can result in many misclassifications. By modelling the labels using a latent GP, this problem does not occur in our method. We apply our novel method to a range of examples including a motivating example which models the hormones associated with the reproductive system in mammals, where the two labelled outputs are high and low rates of reproduction., Comment: arXiv admin note: substantial text overlap with arXiv:1901.07413

Calibration of expensive computer models with high-dimensional output fields can be approached via history matching. If the entire output field is matched, with patterns or correlations between locations or time points represented, calculating the distance metric between observational data and model output for a single input setting requires a time intensive inversion of a high-dimensional matrix. By using a low-dimensional basis representation rather than emulating each output individually, we define a metric in the reduced space that allows the implausibility for the field to be calculated efficiently, with only small matrix inversions required, using projection that is consistent with the variance specifications in the implausibility. We show that projection using the $L_2$ norm can result in different conclusions, with the ordering of points not maintained on the basis, with implications for both history matching and probabilistic methods. We demonstrate the scalability of our method through history matching of the Canadian atmosphere model, CanAM4, comparing basis methods to emulation of each output individually, showing that the basis approach can be more accurate, whilst also being more efficient.

The use of emergent constraints to quantify uncertainty for key policy relevant quantities such as Equilibrium Climate Sensitivity (ECS) has become increasingly widespread in recent years. Many researchers, however, claim that emergent constraints are inappropriate or even under-report uncertainty. In this paper we contribute to this discussion by examining the emergent constraints methodology in terms of its underpinning statistical assumptions. We argue that the existing frameworks are based on indefensible assumptions, then show how weakening them leads to a more transparent Bayesian framework wherein hitherto ignored sources of uncertainty, such as how reality might differ from models, can be quantified. We present a guided framework for the quantification of additional uncertainties that is linked to the confidence we can have in the underpinning physical arguments for using linear constraints. We provide a software tool for implementing our general framework for emergent constraints and use it to illustrate the framework on a number of recent emergent constraints for ECS. We find that the robustness of any constraint to additional uncertainties depends strongly on the confidence we can have in the underpinning physics, allowing a future framing of the debate over the validity of a particular constraint around the underlying physical arguments, rather than statistical assumptions.

In many real-world applications we are interested in approximating costly functions that are analytically unknown, e.g. complex computer codes. An emulator provides a fast approximation of such functions relying on a limited number of evaluations. Gaussian processes (GPs) are commonplace emulators due to their statistical properties such as the ability to estimate their own uncertainty. GPs are essentially developed to fit smooth, continuous functions. However, the assumptions of continuity and smoothness is unwarranted in many situations. For example, in computer models where bifurcations or tipping points occur, the outputs can be discontinuous. This work examines the capacity of GPs for emulating step-discontinuous functions. Several approaches are proposed for this purpose. Two special covariance functions/kernels are adapted with the ability to model discontinuities. They are the neural network and Gibbs kernels whose properties are demonstrated using several examples. Another approach, which is called warping, is to transform the input space into a new space where a GP with a standard kernel, such as the Matern family, is able to predict the function well. The transformation is perform by a parametric map whose parameters are estimated by maximum likelihood. The results show that the proposed approaches have superior performance to GPs with standard kernels in capturing sharp jumps in the true function.

We present a new method of modelling numerical systems where there are two distinct output solution classes, for example tipping points or bifurcations. Gaussian process emulation is a useful tool in understanding these complex systems and provides estimates of uncertainty, but we aim to include systems where there are discontinuities between the two output solutions. Due to continuity assumptions, we consider current methods of classification to split our input space into two output regions. Classification and logistic regression methods currently rely on drawing from an independent Bernoulli distribution, which neglects any information known in the neighbouring area. We build on this by including correlation between our input points. Gaussian processes are still a vital element, but used in latent space to model the two regions. Using the input values and an associated output class label, the latent variable is estimated using MCMC sampling and a unique likelihood. A threshold (usually at zero) defines the boundary. We apply our method to a motivating example provided by the hormones associated with the reproductive system in mammals, where the two solutions are associated with high and low rates of reproduction., Comment: Old version with errors

Ice sheet models are used to study the deglaciation of North America at the end of the last ice age (past 21,000 years), so that we might understand whether and how existing ice sheets may reduce or disappear under climate change. Though ice sheet models have a few parameters controlling physical behaviour of the ice mass, they also require boundary conditions for climate (spatio-temporal fields of temperature and precipitation, typically on regular grids and at monthly intervals). The behaviour of the ice sheet is highly sensitive to these fields, and there is relatively little data from geological records to constrain them as the land was covered with ice. We develop a methodology for generating a range of plausible boundary conditions, using a low-dimensional basis representation of the spatio-temporal input. We derive this basis by combining key patterns, extracted from a small ensemble of climate model simulations of the deglaciation, with sparse spatio-temporal observations. By jointly varying the ice sheet parameters and basis vector coefficients, we run ensembles of the Glimmer ice sheet model that simultaneously explore both climate and ice sheet model uncertainties. We use these to calibrate the ice sheet physics and boundary conditions for Glimmer, by ruling out regions of the joint coefficient and parameter space via history matching. We use binary ice/no ice observations from reconstructions of past ice sheet margin position to constrain this space by introducing a novel metric for history matching to binary data.

Existing methods for diagnosing predictability in climate indices often make a number of unjustified assumptions about the climate system that can lead to misleading conclusions. We present a flexible family of state-space models capable of separating the effects of external forcing on inter-annual time scales, from long-term trends and decadal variability, short term weather noise, observational errors and changes in autocorrelation. Standard potential predictability models only estimate the fraction of the total variance in the index attributable to external forcing. In addition, our methodology allows us to partition individual seasonal means into forced, slow, fast and error components. Changes in the predictable signal within the season can also be estimated. The model can also be used in forecast mode to assess both intra- and inter-seasonal predictability. We apply the proposed methodology to a North Atlantic Oscillation index for the years 1948-2017. Around 60% of the inter-annual variance in the December-January-February mean North Atlantic Oscillation is attributable to external forcing, and 8% to trends on longer time-scales. In some years, the external forcing remains relatively constant throughout the winter season, in others it changes during the season. Skillful statistical forecasts of the December-January-February mean North Atlantic Oscillation are possible from the end of November onward and predictability extends into March. Statistical forecasts of the December-January-February mean achieve a correlation with the observations of 0.48., Comment: 16 pages, 10 figures

Weakly stationary Gaussian processes (GPs) are the principal tool in the statistical approaches to the design and analysis of computer experiments (or Uncertainty Quantification). Such processes are fitted to computer model output using a set of training runs to learn the parameters of the process covariance kernel. The stationarity assumption is often adequate, yet can lead to poor predictive performance when the model response exhibits nonstationarity, for example, if its smoothness varies across the input space. In this paper, we introduce a diagnostic-led approach to fitting nonstationary GP emulators by specifying finite mixtures of region-specific covariance kernels. Our method first fits a stationary GP and, if traditional diagnostics exhibit nonstationarity, those diagnostics are used to fit appropriate mixing functions for a covariance kernel mixture designed to capture the nonstationarity, ensuring an emulator that is continuous in parameter space and readily interpretable. We compare our approach to the principal nonstationary GP models in the literature and illustrate its performance on a number of idealised test cases and in an application to modelling the cloud parameterization of the French climate model.

The calibration of complex computer codes using uncertainty quantification (UQ) methods is a rich area of statistical methodological development. When applying these techniques to simulators with spatial output, it is now standard to use principal component decomposition to reduce the dimensions of the outputs in order to allow Gaussian process emulators to predict the output for calibration. We introduce the `terminal case', in which the model cannot reproduce observations to within model discrepancy, and for which standard calibration methods in UQ fail to give sensible results. We show that even when there is no such issue with the model, the standard decomposition on the outputs can and usually does lead to a terminal case analysis. We present a simple test to allow a practitioner to establish whether their experiment will result in a terminal case analysis, and a methodology for defining calibration-optimal bases that avoid this whenever it is not inevitable. We present the optimal rotation algorithm for doing this, and demonstrate its efficacy for an idealised example for which the usual principal component methods fail. We apply these ideas to the CanAM4 model to demonstrate the terminal case issue arising for climate models. We discuss climate model tuning and the estimation of model discrepancy within this context, and show how the optimal rotation algorithm can be used in developing practical climate model tuning tools.

We develop Bayesian state space methods for modelling changes to the mean level or temporal correlation structure of an observed time series due to intermittent coupling with an unobserved process. Novel intervention methods are proposed to model the effect of repeated coupling as a single dynamic process. Latent time-varying autoregressive components are developed to model changes in the temporal correlation structure. Efficient filtering and smoothing methods are derived for the resulting class of models. We propose methods for quantifying the component of variance attributable to an unobserved process, the effect during individual coupling events, and the potential for skilful forecasts. The proposed methodology is applied to the study of winter-time variability in the dominant pattern of climate variation in the northern hemisphere, the North Atlantic Oscillation. Around 70% of the inter-annual variance in the winter (Dec-Jan-Feb) mean level is attributable to an unobserved process. Skilful forecasts for winter (Dec-Jan-Feb) mean are possible from the beginning of December., Comment: 27 pages, 8 figures, 4 tables

In this paper, we are concerned with attributing meaning to the results of a Bayesian analysis for a problem which is sufficiently complex that we are unable to assert a precise correspondence between the expert probabilistic judgements of the analyst and the particular forms chosen for the prior specification and the likelihood for the analysis. In order to do this, we propose performing a finite collection of additional Bayesian analyses under alternative collections of prior and likelihood modelling judgements that we may also view as representative of our prior knowledge and the problem structure, and use these to compute posterior belief assessments for key quantities of interest. We show that these assessments are closer to our true underlying beliefs than the original Bayesian analysis and use the temporal sure preference principle to establish a probabilistic relationship between our true posterior judgements, our posterior belief assessment and our original Bayesian analysis to make this precise. We exploit second order exchangeability in order to generalise our approach to situations where there are infinitely many alternative Bayesian analyses we might consider as informative for our true judgements so that the method remains tractable even in these cases. We argue that posterior belief assessment is a tractable and powerful alternative to robust Bayesian analysis. We describe a methodology for computing posterior belief assessments in even the most complex of statistical models and illustrate with an example of calibrating an expensive ocean model in order to quantify uncertainty about global mean temperature in the real ocean., Comment: Published at http://dx.doi.org/10.1214/15-BA966SI in the Bayesian Analysis (http://projecteuclid.org/euclid.ba) by the International Society of Bayesian Analysis (http://bayesian.org/)

In this paper we tackle the problem of generating uniform designs in very small subregions of computer model input space that have been identified in previous experiments as worthy of further study. The method is capable of producing uniform designs in subregions of computer model input space defined by a membership function that consists of a continuous function passing a threshold test, and does so far more efficiently than current methods when these subregions are small. Our application is designing for regions of input space that are not ruled out by history matching, a statistical methodology applied in numerous diverse scientific applications whereby model runs are used to cut out regions of input space that are incompatible with real world observations. History matching defines a membership function for a region of input space that is not ruled out yet by observations in the form of a distance metric called implausibility. We use this distance metric to drive a new type of Evolutionary Monte Carlo algorithm with a uniform distribution on the not ruled out yet region as its target distribution. The algorithm can locate and generate uniform points within extremely small subspaces of the computer model input space with complex and even disconnected topologies. We illustrate the performance of the technique in comparison to current methods with a number of idealised examples. We then apply our algorithm to generating an optimal design for the not ruled out yet region of a galaxy simulation model called GALFORM following 4 previous waves of history matching where the target region is 0.001% the volume of the input space.