---
abstract: '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.'
address: 'Department of Mathematics, College of Engineering, Mathematics and Physical Sciences, University of Exeter'
author:
- Louise Kimpton
- Peter Challenor
- Daniel Williamson
bibliography:
- 'REF1.bib'
title: Modelling Numerical Systems with Two Distinct Labelled Output Classes
---

Classification ,Uncertainty Quantification ,labelled outputs ,computer model ,two solutions

Introduction
============

In many areas of science, complex numerical models are used to represent real life physical systems [@Sacks1989]. In general, a mathematical simulator is used to approximate physical reality and the simulator parameters are estimated to specify those models that best represent the real world. We can reproduce data, make predictions and generally get a better understanding of these complex systems by using such models. For practical applications when making predictions, it is also important to include estimates of uncertainty.

In the majority of cases, the inner workings of the system produce similar output results no matter what small changes are made to the model parameters and we can represent the relationship between model inputs and outputs by a smooth, continuous function. In certain applications of scientific modelling, we find that this is not the case; different areas of input space create significantly different values or properties of the output. There will be regions in the model output space where the overall trend has no resemblence to other regions, whether this be in the shape, range or other properties of the output values. These occurences can produce discontinuities between regions in the output space; examples include tipping points and bifurcations. These discontinuities can create step functions at the transitions between regions, so it is important not to assume any continuity between the separate solutions. For example, in climate science, the Stommel model has a different solution for when the overturning circulation is turned on or off [@Sciences2005]. In other cases, the output may be in a binary or catergorical form such as computer code for a complex model that fails to run for certain input values. This corresponds to separate binary outcomes of ’runs’ and ’fails to run’.

A motivating example has been supplied by [@Voliotis2018] where the subject is the reproductive system in mammals. In particular how this is controlled by connections between the brain, the pituity gland, and the gonads. There are particular neurones in the brain that secrete a specific hormone known as the gonadotrophin-releasing hormone (GnRH). These are vital in regulating gametogenesis and ovulation. Signals are made by the pituitary gland which then simulate the gonads for this cycle to start. One of the regulators of the GnRH neurone is neuropeptide kisspeptin, of which two are located within areas of the hypothalamus (the arcuate nucleus (ARC) and the proptical area). Other research suggests that one of these areas (ARC) is the location of the GnRH pulse regulator of which the core are neurones (ARC kisspeptin or KNDy) that secretes two neuropeptides: neurokinin B (NKB) and dynorphin (Dyn). The object of the model presented is to understand the role of NKB and the firing rate of these neuropeptides on the regulation of GnRH, and subsequentially in controlling reproduction. To do this, the model identifies the population of the KNDy neurones where the GnRH pulse regulator is said to be found. The model consists of a set of coupled ordinary differential equations (ODEs) to describe the dynamics of $m$ synaptically connected KNDy neurones. There are several fixed parameters including the concentration of Dyn, rates at which Dyn and NKB are lost and those that describe the characteristic timescale for Dyn and NKB. The variables are the concentration of NKB secreted at the synaptic ends and the firing rate, measured in spikes/min. Using the population of KNDy neurones is shown to be critical for GnRH pulsatile dynamics and that this can stimulate GnRH secretion. Analysing the output of this model shows that the population can behave as a bistable switch so that the firing rate is either high or low. Hence, this causes us to have a system with two distinct solutions, and is an example of the type of system that we wish to model. This bistable system is coupled with negative feedback leading to sustained oscillations that drive the secretion of GnRH hormones that are involved in reproduction. In being able to model the system and locate the areas of low and high firing rates means that not only can we aide predictions on the repreduction rate but we can also have a better understanding of the specific input parameters that are associated with high rates of reproduction.

A common solution to the problem of uncertainty quantification with ’black box’ models is the use of Gaussian process emulation [@Kennedya]. A Gaussian process emulator is a statistical approximation to the simulator and is fast to run so that uncertainty estimates can be calculated. They are particularly useful when evaluation of the simulator is computationally expensive due to the complex nature of the underlying physical system. If we were dealing with simple and cheap to run models, then emulation would be redundant since the simulator could be run many times resulting in easy analysis of the separate output solutions and the system as a whole.

The main aim of this report is to emulate complex systems with multiple solutions, as in the motivating example above. Initially to simplify the problem, systems with exactly two output solutions are considered, generalising later to $n$ dimensions. We define these to be the output regions for the remainder of the paper. The model output can be either discontinuous or in a binary form, so it seems sensible to avoid current stationary methods of emulating data as a whole. When considering Gaussian processes, these types of regression models do not typically cope well with modelling discontinuities or step functions. This is stated by [@Neal1998], where it is said that Gaussian processes are not appropriate priors for models with discontinuities, or where smoothness varys as a function of the inputs. If applying a Gaussian process to a step function, as the height of the step increases in size, we find that the corresponding emulator becomes increasingly inaccurate over the whole function. It not only overshoots near the discontinuity, but also tends to induce fluctuation in the rest of input space as it tries to model this abrupt jump while still preserving the smoothness assumptions.

Some literature focus on non-stationary Gaussian processes that could be applied to this problem. A non-stationary Gaussian process has a covariance structure that varies throughout the input space, where there may be areas of higher variability. This is applicable to models with two ouput solutions since the two solutions are assumed to have different output trends and hence distinct underlying covariance structures. Examples of these include changes to the covariance function [@Schmidt2003], Composite Gaussian Processes [@Ba2012] and Treed Gaussian processes [@Gramacy2008]. Treed Gaussian processes work by partitioning up the input space to fit different models to data independently in each separate region. They specifically divide the input space up by making binary splits on the value based on a single variable so that the boundaries between regions are parallel to coordinate axes. This is an iterative process, such that new partitions are subpartitions of existing partitions. The main problem with treed Gaussian processes is that partitions are made on straight lines parallel to coordinate axes. This is similar to region partioning using Voronoi tessellation [@Gallier2008] introduced by [@Kim2005]. Input space is also partioned similarly with disadvantages due to straight line partions. Both of these methods result in loss of model flexibility and potential errors when boundaries between output regions are not linear.

We therefore conclude that it would be unrealistic to use information from one region to model the other and that current non-stationary Gaussian process models are not suitable for our exact specification. [@Neal1998] suggests using non-Gaussian models or only including a Gaussian process at the lowest level of the model. Taking this into consideration, we will consider estimation of the boundary between regions and modelling the region outputs separately.

First steps in this direction were made by [@Diggle1998] in using a logit transformation to map the domain of a Gaussian process on to the unit interval. With the two regions, identifying them with binary labels (0 for region 1 and 1 for region 2), they consider the probability of which region an outcome may lie on, hence it is appropriate to be able to model to the unit interval and for the data to be Bernoulli. The main aim of the paper is to address the assumption of data being Gaussian and instead concentrate on situations where the stochastic variation in the data is known to be non-Gaussian. Hence, it would seem appropriate to model the data as Bernoulli trials where a success is treated as being in the specified output region with the probability of success following a Gaussian process.

A similar method is mentioned by [@Chang2015] involving ice sheet models and binary data. They propose a novel calibration method for computer models whose output is in the form of binary spatial data. The approach is based on a generalised linear model framework with a latent Gaussian process. It follows the standard logistic regression framework corresponding to the probability for each observation. By assuming the elements in the model output are conditionally independent given the natural parameters, the likelihood function can be found. Construction of the Gaussian process element differs in that the likelihood maximised is now binomial.

The layout of the problem is also very similar to classification in machine learning mentioned by [@Seeger2004] and [@Nickisch2008]. In their formulation, the input data points, $\textbf{x}_{i}$, are associated with separate class regions with corresponding class labels, $y_{i} \in \{-1,1\}$. The process, $f(\textbf{x})$, becomes latent in the model, and is transformed using a sigmoid function, $\sigma$, so that the probability of being in one of the classes, $P(y = +1|\textbf{x})$, can be modelled. The class labels are assumed to be independently distributed Bernoulli random variables. A posterior distribution over the latent values, $f(\textbf{x})$, is found in terms of both the training and test latent values, $f(\textbf{x})$ and $f(\textbf{x}_{*})$. Note here that $\textbf{x}_{*}$ is a test point where the class memership probability is to be predicted. The predictive class membership probability, $P(y_{*} = 1|\textbf{x}_{*}, \textbf{y}, \textbf{X}, \theta)$, is obtained by averaging out the test set latent variables, $f(\textbf{x}_{*})$. The main disadvantage with the method outlined by [@Nickisch2008] is the fact that part of the posterior distribution is not analytically tractable. This is due to the observation likelihood no longer being Gaussian. The rest of this paper outlines ways to tackle this problem by describing different techniques to numerically approximate the posterior distribution for the predictive class membership. These include Laplace Approximation, Expectation Propagation and Kullback-Leibler divergence.

A similar approach to numerically approximate the posterior distribution is shown by [@Chan2013] that follows a Bayesian approach to generalised linear modelling (GLM). Different Gaussian process models are obtained by changing the form of the likelihood, which [@Chan2013] limit to the exponential family. The model is composed of a latent Gaussian process, a random component, $P(y|\theta)$, that models the output as an exponential family distribution and a link function that relates the mean of the output distribution with the latent function. A prior is placed on the latent function, adding in a Bayesian element to the model. A binomial distribution (with $n=1$ for GP classification) is used in the exponential family form and so the mean is related to the latent space through the logistic function.

All of these methods produce a posterior distribution for the predictive class membership of being in one of the two regions. When we sample from this, or use it to make predictions, we draw from an independent Bernoulli distribution where the 0/1 outputs correspond to either of the two regions. When classifying data into two specific regions, apart from directly at the boundary, we state that all points in the neighbourhood belong in the same region. In our previous example of a computer model that does not run to completion at certain input values, if we knew one point where the model is certain to crash, it is sensible to assume that other similar input values are also likely to cause a crash. Hence, correlation between neighbouring points is valuable and should be incorporated into our model. When we take draws from a Bernoulli distribution, each draw is independent from every other draw, so the classification for input values with the same probability of success is equivalent, regardless of any information in the neighbouring area.

Take a simple example with one input variable and four known classified data points, where it is known that there is exactly one change in region somewhere between the two centre points. Thus, we have two points known to be in region one, followed by two points known to be in region two. The change in region can happen anywhere between the two central points. The input space between each pair of points in the same regions, however, must be classified into the same respective region as the surrounding points. If we drew randomly between the two points in region one, we forfeit this known information, and it is thus important to include some correlation over distance in our model. This loss of information would result in random occurrances of points being classified into the wrong region; something that is not intended in the set up of the example.

Given that there are only two output regions, we assume a hard boundary. As we get closer to this boundary, the probability of being classified into the first region is going to become close to 50% since we are uncertain of where the exact boundary lies. Hence, the draws from the Bernoulli distribution become equally likely to fall on a 0 or a 1, so there will be a section (close to the boundary) where the classification may appear fairly random. Therefore, if we wish to pursue in a classification direction, we would require a classifier that included some correlation to help us obtain a clean cut boundary between regions and a ’smooth’ classification.

An interesting aspect of this method however, is the use of latent modelling with corresponding region class labelling. Each data point has two quantities attached to it; the function output and a class label that corresponds to which output region the point lies in. This may be something to consider when we work with models that have no associated output function or where the output is binary. If we have knowledge of which output region data points are in, we may be able to model only the class labels and ignore any corresponding system outputs.

[@Ranjan2008] propose an alternative method to classification and logistic regression by attempting to model the boundary between the two separate outout regions, specifically as a contour. They try to estimate the contour of a complex computer code based on an improvement function. A relatively small experimental design is performed and points are chosen sequentially based on the improvement function weighted towards choosing points on or near the estimate of the contour, or where the predicted variance is high. This process is aided with the use of Gaussian process emulation. Although this method appears to be an improvement in producing the uncertainty, it requires an underlying smoothness assumption. The whole output space is modelled by one single Gaussian process, where there is a simplifying assumption of the response surface being smooth in the form of the covariance function. Therefore, this method would be unsuitable for models with discontinuities. Also, it is likely to become difficult in higher dimensions.

A process known as history matching is used in a method developed by [@Caiado2015]. History matching is an iterative process designed to reduce the input space of the simulator such that input values that are not likely to result in the observed data are discarded [@Andrianakis2015; @Vernon2010]. Here, it is used to sort data into the separate output regions by discarding regions which are unlikely based on an implausibility criteria. Although this has no smoothness assumption, it may still be difficult in higher dimensions.

Overall, it is clear to see that there is a need for a model that can be used for systems with more than one class of output solution.

Gaussian Process Emulation
==========================

A mathematical simulator is often based on the solution of a set of physically justified PDEs and aims to mimic the behaviour of the complex system so that insight can be gained into the functioning of the system. The main disadvantages in many complex simulators are that run times tend to be lengthy and computationally expensive so it is usually not feasible to run large sets of inputs [@Craig2001]. Emulators are statistical approximations to simulators so encapsulate features of the simulator through a complete probability distribution. Consider observations of a simulator, $y$, assumed to be continuous. For a vector of inputs, $\textbf{x} = [x_{1},x_{2}, ..., x_{n}]$, this can be displayed as follows:

$$y(\textbf{x}) = f(\textbf{x}) + \epsilon \hspace{0.1cm},$$

where $f(\textbf{x})$ is the mean value of the output and $\epsilon$ is an error term. Emulators are a non-parameteric approach to regression such that they find a distribution over the possible functions, $f(\textbf{x})$, that are consistent with the observed data.

A Gaussian process is a generalisation of a Gaussian distribution over an infinite dimensional space and is fully defined by a mean function, $m(\textbf{x})$, and a covariance function, $v(\textbf{x}_{1},\textbf{x}_{2})$. [@Kennedya]. If a function, $f(\textbf{x})$, is distributed as a Gaussian process with mean function $m(\cdot)$ and covariance function $v(\cdot,\cdot)$ for any set of inputs $x_{i},...,x_{n}$, the associated finite set of random variables, $f(x_{1}), f(x_{2}), ..., f(x_{n})$, have distribution,

$$\begin{bmatrix} 
f(x_{1}) \\
\vdots \\
f(x_{n})
\end{bmatrix}
\sim
\mathcal{N} \left(
\begin{bmatrix} 
m(x_{1}) \\
\vdots \\
m(x_{n})
\end{bmatrix},
\begin{bmatrix} 
v(x_{1},x_{1}) & \dots & v(x_{1},x_{n}) \\
\vdots & \ddots & \vdots \\
v(x_{n},x_{1}) & \dots & v(x_{n},x_{n}) 
\end{bmatrix}
\right).$$

All marginal, joint and conditional distributions are Normal [@Seeger2004]. This can then be formally written as:

$$f(\textbf{x}) \sim GP(m(\textbf{x}), v(\textbf{x}_{i},\textbf{x}_{j})) \hspace{0.1cm},$$

We will restrict ourselves to Gaussian processes with linear prior mean functions, so the prior mean can be specified as $E[f(\textbf{x})|\boldsymbol{\beta}] = \textbf{h}(\textbf{x})^{T}\boldsymbol{\beta}$, where $\textbf{h}(\textbf{x})$ is a vector of basis functions of $\textbf{x}$ and $\boldsymbol{\beta}$ is a vector comprising of unknown coefficients. $v(\textbf{x}_{1},\textbf{x}_{2})$ is a covariance that can be defined as $\sigma^{2}c(\textbf{x}_{i},\textbf{x}_{j})$, where $c$ is a known correlation function of distance. The stationarity of the Gaussian process means that the covariance function does not change over the input space. A common choice of correlation function is the squared exponential; $c(\textbf{x}_{i},\textbf{x}_{j}) = \text{exp} \left\{ - \frac{| \textbf{x}_{i} - \textbf{x}_{j} | ^{2}}{\delta} \right\}$, where $\delta$ is the correlation length parameter and controls the wiggliness of the process.

Model Outline
=============

Initially, we include no knowledge of the function output of the system, simply a label on which region the output belongs to, i.e. whether the input produces an output in region 1 or region 2. This is important since the range of possible applications includes those where the output is simply binary and hence has no associated output value (for example, where the code does not run to completion). Instead of focusing on the function output, we will instead be using a class labelling (similar to that of classification [@Nickisch2008]) that we assign according to whether the input point lies in region 1 or region 2.

In order for the model to be applicable to cases that have either continous or discontinous transitions between regions, we aim to model the output regions separately. Therefore, we need to find where the boundary between these regions lies and be able to make predictions of region classification for other input values. This should include quantifying any uncertainty that is present with these estimates. Once the regions have been identified, if the actual system has real output values that we are interested in, then the regions can be modelled with separate Gaussian process emulators.

Due to limited initial information, the choice was made to model the boundary and region classification in latent space. We also need to ensure that correlation between data points is included in the model. Hence, a latent variable modelled as a Gaussian process is used to structure the two output solutions using our assigned class labelling. A latent variable is a quantity that is not observed directly, but rather inferred from the region classifications that are the observations. The variable is hidden behind the physical model and used to aid our prediction of the boundary. The values of the Gaussian process will not be measured themselves; only the values of the system and which region the points lie in.

As with logistic regression and Gaussian process classification, we are interested in adapting their use of class labels, $y$, for each input data point, $\textbf{x}$, associated with the separate class regions. In contrast to [@Nickisch2008] and [@Seeger2004], who use the values -1 and 1 to denote the different regions, we adapt this to account for minimal knowledge known about the form of the latent variable. We do not assign numerical values as class labels for the state, but simply let the labels be negative for one region and positive for the other. It is important to distinguish that for each input point we have a class label and a separate function output if the system is such that an output exists. We do not use the function output explicitly in the modelling of the two regions. A Gaussian process, $\eta(\textbf{x})$, is estimated with the constraints of having the correct sign for all initial known data points. This is different to a normal Gaussian process since the training data is usually the input variables and the corresponding function output for a few predesigned points. The latent aspect no longer trains the Gaussian process to lie at the function evaluations at the specific input points, but are allowed to take any value providing the sign class labels agree at these inputs. Once this is satisfied and the latent Gaussian process has been estimated, a threshold, $\psi$, is taken to split the input space into the separate regions. Due to the set up of the problem, the theshold is taken to be zero.

Metropolis Hastings
-------------------

For the estimation of the paarameters of the latent Gaussian process, a method is required to provide estimates that agree with the negative and positive region labels that have been defined. We take a Bayesian approach.

Bayesian inference derives a posterior probability distribution from a prior distribution and a likelihood function,

$$\pi(\theta|x) = \frac{\pi(x|\theta)\pi(\theta)}{\pi(x)} \hspace{0.1cm}.$$

$\pi(\theta|x)$ is the probability density function of the model parameters, $\theta$, given the data, $x$. This is made up of the prior knowledge of the model parameters before knowledge of the data, $\pi(\theta)$, along with the likelihood of the data, $\pi(x|\theta)$. What is of most interest is the denominator since this is where non-tractable problems lie. This quantity can be found by integrating $\pi(x,\theta)$ over all possible parameter values. However, this calculation is only analytically possible for very simple models, hence we turn to approximation methods by drawing samples from the posterior.

Markov chain Monte Carlo (MCMC) methods [@Brooks2012] refer to a collection of algorithms that are designed to approximate random samples from probability distributions that would cause problems if trying to sample directly. The result of such an algorithm tends to a Markov chain whose equilibrium distribution matches that of the distribution of interest. See [@Gelman2013] and [@Sivia2006] for more background on Bayesian methods.

One form of Monte Carlo sampling is known as rejection sampling or the acceptance-rejection method. A simple form of rejection sampling when starting with minimal information is Approximate Bayesian Computation (ABC). ABC works by simulating predicted model data and comparing with the known observations to estimate the posterior distribution for the model parameters [@Turner2012]. Hence, it will produce an estimate for the expected posterior of the latent Gaussian process assuming only the input values and some form of output comparison; we do not need to supply a likelihood. From [@Turner2012], model parameters, $\theta$, are accepted if $\rho(X,Y) \leq \delta$, where $\rho(.,.)$ is a measure of distance, usually taken to be the Euclidean distance, $\Vert X-Y \Vert$ (X being the simulated data from $\theta$, Y being the true observations and $\delta$ is a small quantity that specifies how close the approximation is to the true posterior distribution) [@Wilkinson2008].

We propose a modified version of the ABC MCMC algorithm for our problem. This produces an estimate to the posterior distribution for the latent GP parameters, $\theta = (\boldsymbol{\beta},\sigma,\delta)$. As mentioned previously, each data point will be given a class label of positive or negative. This concept is highly important in this algorithm since it will be used as the basis of the rejection criteria. So, instead of requiring samples to be close in value to the observations, we are forcing samples to have the same sign as the observation for comparison. Hence, it will be generating an estimate of the Gaussian process, $\eta$, that is negative in all input space of region 1 and positive in the input space of region 2.

Although using this method of ABC sampling proves to be effective, it is not very efficient and could be improved to include more information already known from the initial data. In a rejection sample such as ABC, the rejection rate tends to be very high, especially when working in higher dimensions. In alternative MCMC methods this is not always the case, and so we proceed by considering a version of the Metropolis Hastings algorithm.

The general Metropolis Hastings algorithm works by constructing Markov chains such that its stationary distribution is the distribution of interest. See [@Gilks1996] and [@Chib2017] for more information. Unlike ABC, Metropolis Hastings requires a likelihood. In this case it needs to reflect the knowledge of which points are definitely positive or negative in association with the different regions. Hence, the likelihood becomes:

$$\begin{split}
\mathcal{L}(\theta;x) & = P \left( \eta(x_{1})<0 , \eta(x_{2})<0 , \ldots , \eta(x_{j})<0 ,  \eta(x_{j+1})>0 , \ldots , \eta(x_{n})>0 \right) \\
& = \int_{-\infty}^{0} \dots \int_{-\infty}^{0} \int_{0}^{\infty} \dots \int_{0}^{\infty} \phi ( \eta(x_{1}), \eta(x_{2}), \dots, \eta(x_{j}), \eta(x_{j+1}), \dots, \\
& \hspace{1cm} \eta(x_{n}) )  \hspace{0.2cm} dx_{1} dx_{2} \dots dx_{n} .
\end{split}$$

This gives a joint distribution of the first $j$ points falling in negative space and the second $n-j$ points remaining in positive space. The main difference between sampling from an ordinary Gaussian process, to the latent one here is the use of the cumulative distribution function instead of the density function of a Normal distribution. By specifying this to be our likelihood, we are not putting any constraints on the specifc values of the generated Gaussian process, just their sign values. The process is not observed and so there is no concern for exactly what value it takes as long as the change in sign is estimated correctly to give us an estimation for the region boundary. It is also important to note the correlation assumption in using this likelihood. Since nearby points are more likely to fall in a similar region, the correlation must be accounted for through the use of a single multivariate draw of the distribution.

As well as the use of a likelihood, Metropolis Hastings also differs from ABC through the calculations to the final predictive distribution. Where ABC simulates many samples and rejects those that do not fit the criteria, Metropolis Hastings uses the given likelihood to downweight parameters that are less likely to give a sample that fits the comparison with the observations. So, as the iterations progress, the current parameter estimates in the chain progress closer to those that are more likely to result in accepted samples. In otherwords, a certain amount of rejection is done automatically with the use of the likelihood.

At the end of the algorithm, we are left with a chain of values that form an estimate to the posterior distribution of the Gaussian process parameters. We can then find a MAP (Maximum A Posteriori) estimate to give the parameters that maximise this posterior distribution [@Rice2007]. Using this, we can also produce a corresponding posterior distribution for the latent values themselves.

In proceeding with finding the MAP estimate of the parameters, this results in drawing many samples from a joint Gaussian distribution. This is inefficient and highly time consuming since not only do we have to make the draws, but it is also necessary to perform a process similar to a rejection sample to eliminate any of the remaining samples that do not follow the negative and positive region trend. One way in which the efficiency can be increased is to sample each latent value in $\eta(\textbf{x})$ in turn using the normal conditioning equations presented below:

$$\begin{split}
y_{j}|y_{i} & \sim \mathcal{N}(E[y_{j}|y_{i}],\text{var}[y_{j}|y_{i}]) \hspace{0.1cm} , \\
E[y_{j}|y_{i}]& = E[y_{j}] + \text{cov}[y_{j},y_{i}]\text{var}[y_{i}]^{-1}(y_{i} - E[y_{i}]) \hspace{0.1cm} , \\
\text{var}[y_{j}|y_{i}]& = \text{var}[y_{j}] - \text{cov}[y_{j},y_{i}]\text{var}[y_{i}]^{-1}\text{cov}[y_{i},y_{j}] \hspace{0.1cm} , 
\end{split}$$

where $i,j = {1,2,...,n}$. Now, to find a single draw for the latent process, we are only drawing $n$ times from a univariate Gaussian distribution instead of a vector of length $n$ from a multivariate distribution. Using the equations above ensures that the correlation between points remains included in the sample and computational time is greatly saved. On sampling each latent process point in turn, draws from a univariate Gaussian distribution are sampled until the value agrees with the sign of the corresponding region. Hence $\eta(x_{1})$ is first sampled and all values that are positive for region 1 are rejected until a negative value is sampled, stored as the first latent variable point and we move on to generating a value for the next point. The second, $\eta(x_{2})$, is then drawn conditional on the results from the first sample, $\eta(x_{1})$, and again values are rejected where appropriate. This then continues for all data values and ensures only one value is being sampled at a time, hence reducing computation time. Not only is time saved in this way, but also since we are able to reject points that do not concide with the correct region separately at each point. We never reject the whole of a sample, or points that have been shown to agree with the sign for that region. Instead we are resampling each point in turn until a valid value is found.

The main computational expense comes from sampling the first point over the boundary in region 2, i.e. the first time the latent process switches from negative to positive. This jump in changing sign can cause problems in the computation, especially if the resulting variance becomes very small and results in the Gaussian process finding it hard to make the initial jump. To ease the computation, it is useful to sample the boundary and make it obvious near the beginning. It is possible to do this since the numbering of the data points are arbitrary and so the ordering makes no difference to the resulting Gaussian process.

Examples 1
==========

1d Example
----------

To illustrate the concept, a simple toy example of one input variable and one output with two solution regions is presented in Figure \[fig1\]. The inputs are a vector of 12 values ranging between 0 and 20, with the undefined boundary situated in the region $[6,8]$. The points in the range $[0,6]$ are specified as region 1 and have associated class label negative, and the points in the range $[8,20]$ are in output region 2 with class label positive.

The latent Gaussian process is found by applying the Metropolis Hastings approach explained above. We find a posterior distribution for the parameters, $\theta = (\beta,\sigma^{2},\delta)$, by running our algorithm a large number of times to ensure convergence . Samples from the posterior distribution can be taken for the latent process at input points, $\eta(x_{1}), ..., \eta(x_{n})$. After an estimate for the latent process has been found, it is then thresholded at $\eta=0$ to give a value of $x$ for where the boundary between regions lies. This is $x = 7.15$ and is shown by the green dashed line in Figure \[fig1\]. This is a suitable value for the boundary since we set the boundary interval to be $[6,8]$ in the example. Due to the little information we have regarding the location of the boundary and lack of knowledge of the actual system output, we would expect there to be a high level of uncertainty in any results. As shown in Figure \[fig1\], the credible intervals for the estimate are very large and are roughly equal to the extreme bounds that we set the example up with, $[6,8]$. It is interesting to note that our estimate for the boundary is not that of the naive estimate of $x=7$. The estimate has favoured that of the right boundary and so gives us some level of confidence that our models are giving a more informative estimate than that of the centre value between data points. More investigation is needed to test this assumption.

![1 dimensional example with 2 output regions. The posterior mean of the latent Gaussian process (blue) is shown along with the prior mean (red) and boundary estimate (green). Both have 95% credible intervals included. Initial data points are shown in orange with size corresponding to misclassification.[]{data-label="fig1"}](pplot1.png)

Misclassification
=================

The method of model validation used in this example is based on a leave-one-out cross-validation. This usually involves leaving each training point out in turn, fitting a Gaussian process to the remaining points, and then using this to predict the point that was left out [@Seeger2004]. Given that the problem is set up in latent space, it is not possible to strictly follow this layout. Instead, we have looked into methods commonly used in Gaussian process classification to adapt the validation to our problem. As mentioned in [@Seeger2004], it is possible to look into the misclassification of the points during the validation. They provide a binary classification example of sorting images of digits 3 and 5 in the postal service. From a test set of data, they count the number of times a digit is wrongly classified in terms of the process’ standard deviation and length scale parameters.

Adapting this to suit our particular problem, we use a version of leave-one-out to find the misclassification rate. A leave-one-out cross-validation is performed on samples from the posterior distribution for the latent process points to predict the sign of each point left out in turn. We use the latent Gaussian process to predict the sign of the removed point only. From these samples, we calculate the proportion of times each point is classified into the wrong region. The output from this is shown in Figure \[fig1\], where the size of the data points corresponds to the rate of misclassification. As expected, the rate is largest for the two points either side of the boundary. In a 1d example such as this, these points are the most critical since they are the points that restrict the boundary to the precise region of input space. It is also interesting to note here that the remaining points have a misclassification rate of almost (but not quite) zero. From a more in-depth look, we can see that very occasionally the latent process crosses the axis. This is caused by the Gaussian process having a short correlation length parameter; leading to the latent process having the chance to bend quickly over the $x=0$ threshold between known points in the same region. We also identify that the prior mean function placed on the Gaussian process has some influence. The prior mean (as shown by the red line in Figure \[fig1\]) forces the posterior points in each region to follow the same pattern. Therefore, this linear effect appears to force the points to stay in the specified sign.

Based on this example, it has become clear that it is important to place suitable priors on the model parameters and the prior mean function. For the example shown in Figure \[fig1\], we found that a linear prior mean function was a particularly accurate choice when generating our Gaussian process.

One interesting aspect of Gaussian processes is their behaviour in the far field of the input space. They converge to the prior mean asymptotically, which would be a problem for, say, a constant prior mean. If a constant mean function is placed on the Gaussian process, then we start to observe the overall latent process tending towards the prior mean in the edges of our input space. Due to our model layout of negative and positive space, the prior mean would be estimated to be close to the $x=0$ axis and we find that it is very easy for the process to switch signs when the latent variable lies close to zero in the far field, forcing a misclassification in that area.

In a situation such as that in Figure \[fig1\], we have extra knowledge that there are only two output regions and so it should not be the case for the latent variable to slip back over the x-axis. This prior knowledge is incorporated into the method by the selection of the prior mean function. The linear placed on the example in Figure \[fig1\] forces the latent Gaussian process away from the $x=0$ axis in the edges of the input space since we are certain that there are only two distinct regions. If, for example, the two regions in the input space were separated by a circle, then it would be sensible to place a quadratic prior on the mean function, ensuring the Gaussian process would not return to the $x=0$ axis in the far field around the circle. Although this would appear to be a sensible choice, polynomials of a higher order come with a larger number of estimated parameters. Therefore, we should consider whether the classification in the edges of our input space is useful or not.

The prior mean in this example has a linear form of $ax + b$, where $b$ is the value of the latent process when it crosses the vertical axis. Since, it is known that the latent process must cross the $x$-axis at approximately the boundary between regions, we can encorporate this into the prior knowledge of our model. This will then help approximate where the latent crosses the vertical axis, $b$, more efficiently. A transformation is applied to the input points so that the boundary between regions in the $x$-axis now approximately lies at zero in the vertical axis. With this transformation, a tight prior can be placed over the axis intercept, $b$, ensuring the latent process crosses the axis at zero. If we contrast the plot in Figure \[fig2\] compared with that in Figure \[fig1\], we notice a significant difference in the resulting latent process. The prior means for each plot are shown in red. Figure \[fig1\] uses the transformed data and is shown to have an expected mean Gaussian process follow its prior by being close to a straight line through the boundary; something that is expected when there is a problem such as this with minimal information. Figure \[fig2\] does not include the transformation and is shown to differ by the posterior estimate in region 1 leveling out as it approaches zero. This is clearly not appropriate since with no information of the system input, we would expect both sides of the latent process to match. This shows that the transformation in the data greatly improves the estimate in the latent process and any predictions that would follow.

![Same example as of figure \[fig1\], but where the data are not transformed. The prior mean (red) crosses close to the origin (0,0).[]{data-label="fig2"}](pplot2.png)

When considering prior knowledge, it is also important to chose a correct form for the correlation length parameter, $\delta$. The correlation length parameter determines how much the Gaussian process is allowed to bend between each of the initial data points [@Seeger2004]. Particularly when considering the 1d example in Figure \[fig1\], we know that there is only one boundary where the latent Gaussian process is not expected to change signs between data points (apart from the boundary between regions). If the correlation lengths are allowed to become too small, then there is a chance that the Gaussian process would be able to curve round quickly and fall briefly in the wrong sign, causing a misclassification of regions in some input areas. To ensure this does not happen, particular inverse gamma priors are placed on the $\delta$’s so that they are forced away from zero and being too small. An inverse gamma prior is also placed on the variance, $\sigma^{2}$.

Examples 2
==========

2d Example
----------

The method is now expanded to include a simple 2 dimensional version of the 1d example. The general method used is very similar to that of 1 dimension, with the exception that the exact boundary is no longer so easily found. In the example in Figure \[fig3\], 20 input points are generated using a Latin hypercube [@Welch1992] over the region $[-1,7]$. The boundary between regions is defined as the line $x=3$. From Figure \[fig3\], the yellow points are those initial points in region 1 (input space $x < 3$) and the purple points are those in region 2 (input space $x > 3$). The latent GP has been applied to a grid of points over the input space to show where the estimated region lies.

![2 dimensional example where the two region are split by an $x_{1}=3$ plane (red). The dark blue region corresponds high probability of be classified into region 1, whilst light blue corresponds to high probability of being classified into region 2. A misclassification rate is also shown based on point size.[]{data-label="fig3"}](pplot3.png)

To show uncertainty within the 2d example, Figure \[fig3\] shows the probability of input points being in region 1 compared to region 2. The dark blue points represent high probability of being classified into region 1 and light blue represents high probability of being classified into region 2. As expected there is high uncertainty of the predicted region around the boundary. This is due to the minimal information known in this area of input space and so it is easy for an incorrect classification. A misclassification rate is calculated for each point in the same way as the validation performed on the 1d example and also shows high uncertainty in this boundary region. We can see from this plot that the points near the boundary have a larger rate of misclassification than anywhere else in the input space. Where the uncertainty increases here, there is a much higher chance of the latent Gaussian process to flip signs. What is also interesting to note here is the general slope of the estimated boundary. If our method was no more efficient than taking the naive approach to the problem, then we would expect the approximate boundary to misclassify equally across the boundary. This is not the case since the upper section is shown to curve far more into region 2 than that of the lower section. We can see that the latent process is in fact listening to the data.

Another 2d Example
------------------

Another example has been provided by Santner et. al, where their test function is shown in Figure \[fig8\] and has the following form:

$$y(x) = 
     \begin{cases}
       \infty &\quad \text{if} \hspace{0.2cm} x_{1}^{2} + x_{2}^{2} \le c_{1}^{2} \\
       \frac{e^{-(a'x + x'Qx)}}{(x_{1}^{2} + x_{2}^{2} - c_{1}^{2})} &\quad \text{if} \hspace{0.2cm}  c_{1}^{2} \le x_{1}^{2} + x_{2}^{2} \le c_{2}^{2} \\
       - \infty &\quad \text{if} \hspace{0.2cm} x_{1}^{2} + x_{2}^{2} \ge c_{2}^{2} ,\\
     \end{cases}$$

where,

$$a = [3,5] \hspace{0.5cm} 
Q = 
\left(\begin{array}{cc} 
2 & 1.5 \\
1.5 & 4  \\
\end{array}\right)
\hspace{0.5cm} c_{1}^{2} = 0.25^{2}, \hspace{0.3cm} c_{2}^{2} = 0.75^{2} \hspace{0.1cm}.$$

The space between the two circles is region one and the remainder is region two, both over the input space $[-1.25,1.25]$. There are function values associated with this example that is only feasible within region 1. We choose to neglect these and just focus on the classification side of our method, where the regions can be modelled separately after they have been fully classified.

![2 dimensional example with two regions. Region 1 lies within the two circles and region 2 is the remaining input space.[]{data-label="fig8"}](testfunction.png)

To set out this problem, 50 data points have been selected where they are given the class label of positive (region 1) if they lie between the two rings, and negative (region 2) otherwise. These are shown by the purple and orange points in Figure \[fig7\] along with the hard boundary (red). The classification after applying our method is also shown in the plot including uncertainty. As in the previous example, the light blue points represent high probability of being sorted into region 1 and the dark blue points shows high probability of being sorted into region 2. The largest areas of uncertainty correspond to the areas where our classification method performed the poorest.

Overall, our method is estimating the regions well with only a few larger deviations in the upper left and right sections of the doughnut region. This is likely to be caused by the lack of information in these areas. Due to the more complicated shape, we chose to fit a constant prior mean function. Overall, this has proven to be successful since no areas have been misclassified in the far corners of the input space, which may have been likely. Alternatively, a quartic polynomial could be used for the prior mean function, but this causes problems with a large number of parameters to be estimated, especially in areas of sparse data.

The two input points that are interesting to point out are those at the bottom of the larger circle; they are classified in different regions but are very close together. In this area, the latent Gaussian process has to change sign quickly but has proven to do so successfully; it is useful to confirm that our method can cope with cases such as this.

A misclassification rate is also included, where the points are more likely to misclassify in region 1 (between rings). This is likely to be due to a higher proportion of points being in region 2 and so the majority of the latent process is negative, making it more likely for points to be classified into region 2. This is supported by the constant mean function estimated to be -2.25; pulling the process to be overall more negative.

![Estimated regions for the 2d example shown in figure \[fig8\]. Initial data points are displayed (orange - region 1 and purple - region 2), with the actual region boundaries shown in red. Uncertainty on the estimate is included where light blue areas correspond to high probability of being classified into region 1 and dark blue areas correspond to high probability of being classified into region 2. A misclassification rate is also shown.[]{data-label="fig7"}](fplot2.png)

Motivating Example
------------------

Returning to the initial motivating example looking into reproduction rates in mammals, the inputs are NKB concentration and firing rate, where a Latin hypercube has been created over the input space of $[0.1,0.2]\times[10,200]$. The choice was made here to transform the data to a $0,1$ scale for computational simplicity. The system is bimodal, so for the 20 initial points, we labelled 5 of them as negative in region 1, and 15 of them as positive in region 2. As with the rest of our examples, we aim to predict the region of any points in the input space and model the system as a whole.

One of the most important choices to be made in this example was the form of the prior mean on the latent Gaussian process. To make the decision, we consulted the expert in the system as well as examining the initial points. From the initial points (yellow and purple) shown in Figure \[fig5\], we can notice that there is likely to be only one change in region, meaning that we have just two disjoint regions seprated by a fairly linear boundary. This agrees with prior knowledge collected from the expert. Therefore, it would be reasonable to either choose a constant or linear mean form. A constant mean would be preferable if we were unsure about the number or shape of regions, but would have a higher chance of misclassifying in sparse areas or in the far field of the input space. We therefore chose to proceed with a linear mean prior. The output of the predicted region boundary from a linear prior mean is shown in Figure \[fig5\] along with the actual boundary (red) and uncertainy. In general, our solution classifies correctly in most areas where, as expected, the area between the regions is the most uncertain. The main issue is that the region of highest uncertainty (where we expect our boundary to lie) does not appear to capture the correct curve of the actual boundary. This is down to lack of information in the area, but since the correct boundary still lies in our area of uncertainty, we can conclude that our model is doing a reasonable job. Misclassification is also shown through the size of the points where it is easier to misclassify points near the boundary. The method is the same as described in previous examples.

![2 dimensional example looking at the effects of hormone release on mammal reproduction. The system has two regions of high and low rates of hormone release where the actual boundary is shown in red. Initial points are displayed (orange - region 1 and purple - region 2), with predicted region classification and uncertainty. Dark blue areas correspond to high probability of being classified into region 1 and light blue areas correspond to high probability of being classified into region 2. A misclassification rate is also shown.[]{data-label="fig5"}](pplot7.png)

Discussion
==========

We have developed a new method for classifying the output of numerical models into one of two classes. Our method is suitable for modelling systems with two output solution classes that are either stationary across the input space, or where there are discontinuities. This includes systems where the output may be in a binary or catagorical form. A major disadvantage of most common methods of classification like those shown by [@Chang2015], [@Seeger2004] and [@Nickisch2008], involves assuming that class labels associated with input points are independently distributed Bernoulli random variables. This is something that causes concern due to any correlation between nearby points being ignored. Neighbouring input points are more likely to result in the same output region, so it is vital that we include this information in our model.

Keeping this in mind, we use aspects of classification from [@Nickisch2008] in the form of class labelling. We use a form of this for limited knowledge and define all input points in one region to have the class label negative and all inputs in the other region to have class label positive. To ensure that correlation between data points is included, a latent variable modelled as a Gaussian process is used to structure the two output solutions using our assigned class labelling. The latent Gaussian process is estimated using a version of the Metropolis Hastings algorithm. As a form of model validation, we have calculated a version of the misclassification rate as shown by [@Seeger2004]. This is based on a leave-one-out cross-validation.

We feel that this method will be applicable to a wide range of applications in computer science, climate science and biology. Our main motivating example is based on assessing reproduction rates in mammals [@Voliotis2018]. We have successfully modelled this bimodel system, where it can be used for class prediction for other inpur points with estimates of uncertainty included.

There are some obvious extensions to the work presented in this paper. One would be to now expand the method to cope with situations when there are more than two output solution classes. This would then increase the numer of applications where it is suitable. There is also room for research in areas of experimental design where we can improve the accuracy of our class classification and boundary estimation with limited intial data.

Acknowledgements {#acknowledgements .unnumbered}
================

We would like to thank T. Santner et al. and M. Voliotis et al. (2018) for providing us with the given examples. I would also like to thank EPSRC for my studentship, allowing me to carry out this research.

References {#references .unnumbered}
==========
