Non-Gaussian Process Regression
Authors: Yaman Kindap
Simon J. Godsill
Venue: Submitted to ICLR 2023 (2023)
Keywords: Non-parametric regression
Bayesian methods
Approximate Inference
Levy processes

Abstract: Standard GPs offer a flexible modelling tool for well-behaved processes. However, deviations from Gaussianity are expected to appear in real world datasets, with structural outliers and shocks routinely observed. In these cases GPs can fail to model uncertainty adequately and may over-smooth inferences. Here we extend the GP framework into a new class of time-changed GPs that allow for straightforward modelling of heavy-tailed non-Gaussian behaviours, while retaining a tractable conditional GP structure through an infinite mixture of non-homogeneous GPs representation. The conditional GP structure is obtained by conditioning the observations on a latent transformed input space and the random evolution of the latent transformation is modelled using a Lévy process which allows Bayesian inference in both the posterior predictive density and the latent transformation function. We present  Markov chain Monte Carlo inference procedures for this model and demonstrate the potential benefits compared to a standard GP.

Under review as a conference paper at ICLR 2023
NON-GAUSSIAN PROCESS REGRESSION
Anonymous authors
Paper under double-blind review
ABSTRACT
Standard GPs offer a flexible modelling tool for well-behaved processes. How- ever, deviations from Gaussianity are expected to appear in real world datasets, with structural outliers and shocks routinely observed. In these cases GPs can fail to model uncertainty adequately and may over-smooth inferences. Here we extend the GP framework into a new class of time-changed GPs that allow for straightforward modelling of heavy-tailed non-Gaussiag behaviours, while re- taining a tractable conditional GP structure through an infinite mixture non- homogeneous GPs representation. The conditional GP structure is obtained by conditioning the observations on latent transformed input space and the random evolution of the latent transformation is modelled using Lévy process which allows Bayesian inference in both the posterior predictive density and the latent transformation function. We present Markov chain Monte Carlo inference proce- dures for this model and demonstrate the potential benefits compared to a standard GP.
INTRODUCTION
Gaussian processes (GPs) are stochastic processes which are widely used in nonparametric regres- sion and classification problems to represent probability distributions over functions (Rasmussen & Williams (2006)). They allow Bayesian inference in a space of functions such that consistent uncer- tainty measures over predictions are obtained rather than only point estimates. In its simplest form GP defines a distribution over functions through its particular mean and covariance (kernel) func- tions which determine the smoothness, stationarity and periodicity of random realisation in the function space. As a prior distribution in Bayesian inference, using a zero mean GP reflects the lack of information the values and trend of the function. In this case the covariance function, which defines the similarity between any two points in the input space, fully characterises the properties of the random function space.
The design of kernel functions that are able to represent wide range of characteristics and make consistent generalisations a fundamental area of research. Some recent work in this area include modelling the kernel via spectral densities that scale-location mixtures of Gaussians (Wilson & Adams (2013)) and similarly using Lévy process priors over adaptive basis expansions for the spectral density (Jang al. (2017)). Furthermore, extensions to the standard GP model can be made by directly modelling the covariance matrix as stochastic process (Wilson & Ghahramani 201I)), assuming heteroscedastic noise on the observations and carrying out variational inference LázaroGrredilla & Titsias (2011)), or learning nonlinear transformations of the observations such
that the latent transtormed observations are modelled well by GP (Snelson et al. (2003); Lázaro-
Gredilla 2012)). Nonstationarity in the measurement process can be expressed as product of multiple GPs (Adams& Stegle (2008) and heavy-tailed observations may be modelled through the
Student-t process (Shah et al. 20).).. Particularly relevant extensions of GP models are presented
in (Rasmussen & Ghahramani (200T) where the input space is locally modelled by separate GPs,
and string GPs (Samo & Roberts (2016)) introduce link functions between local GPs such that the global process is still a GP and provides efficient inference methods on large data sets. In (Schmidt) & O'Haaa (2003); Snoek et al. (2014)) latent space is defined between the inputs and observations through a separate GP and a class of bounded functions in [o. 1]. respectively.
By designing expressive covariance functions stacking multiple GPs in structured arrangements, the GP framework produces accurate predictive models in numerous application domains. How- ever, these models are limited by their Gaussianity assumption such that the local patterns learned

---

Under review as a conference paper at ICLR 2023
through these models are highly dependent on particular observations instead of learning the overall dynamics of the data generating system. A more natural and interpretable way to define complex relationships may be assume that the underlying random function is Gaussan which yields more sparse representations (Unser & Tafti (2014) as discussed in Section 4
In this work, we present a novel approach to modelling non-Gaussian dynamics by constructing non-Gaussian process (NGP) such that the observations form a conditional GP that is conditioned on a latent input transformation function that is separately modelled as Lévy process. Building on the definition of a stationary kernel, the latent layer between the input and output spaces represent the random distances between any two points on an input space, In order to define the distribution of random distances without referring to specific origin, and in order to maintain monotonicity of the input space transformation, the latent space of transformation functions modelled by a special class of Lévy process called a subordinator that is non-negative and non-decreasing. Such process characterised by the distribution of its stationary and independent increments which as result defines probability distribution over the distance between any two input values. Making random monotonic transformations of input values allow the kernel to adapt to the local characteristics of an input space in other words to the varying rate of change in the observations and the learned subordinator provides uncertainty estimates over the variation of the observed process everywhere on its domain. In this paper we focus principally upon 1-dimensional GPs for the sake of brevity, but we emphasise that our approach can be readily extended to multiple dimensions, as described throughout the text and illustrated in the experimental results.
NGPs are related to continuous-time stochastic volatility models studied in the mathematical fi- nance literature to model the behaviour of a stochastic process which has randomly distributed variance (Ghysels et al. (1996). The time-change operation defined for continuous- processes a standard approach to building stochastic volatility models. A common example the time-changed Brownian motion where the time-change is chosen to be a subordinator and the time- changed motion produces a Lévy process (Veraart Winkel| (2010). In such model, the process is conditionally a Brownian motion i.e. the integral of nhite-noise GP. Similarly, our construction of a stationary NGP follows a GP conditioned on the latent values of a subordinator, thus it is a time-changed GP. Particular non-Gaussian behaviour can be expressed through the characterisation of a subordinator, examples include the stable law, normal-tempered stable, and generalised hyper- bolic (including Student-t) processes. Hence, NGPs provide flexible and expressive probabilistic framework for nonparametric learning of functions.
stochastic
In Section we briefly review the GP regression framework, introduce the time-change operation and define NGPs. An inference method for NGP regression is presented in Section 3 following an introduction shot-noise simulation methods for Lévy processes. In Section 4 present the results of applying NGP regression on representative synthetically generated non- Gaussaan data sets to visually highlight their dynamics and compare the results to alternative GPs. Furthermore, multidimensional example using a data set available in TensorFlow is presented.
2 MODELS
In this section, we briefly present the standard GP regression framework and introduce the time- change operation which results in a non-Gaussian process (NGP). The series representation of a Lévy process (Rosiáski (200T) reviewed in Section 2.2 is central to the inference methodology studied in Section 3
2.1 GAUSSIAN PROCESS REGRESSION
A stochastic process (f(x) R; E X²} is defined by the probability distribution of all possible finite subsets of its values, where E Rdd is d-dimensional input space. In the case of GPs, for any finite set of inputs {x; in the corresponding values of the function (f+))]] has multivariate
Gaussian distribution (MacKay (2003) characterised by its mean m(x) = [z)( and covariance kernel functions K(x,, x) Cov (f(x). fz)) where x', E X. Given set of inputs {2+} the mean function forms a vector m and the kernel function forms a positive-definite covariance matrix
Σ. The resulting multivariate Gaussian distribution can be extended to any input xx E X follow-
2

---

Under review as a conference paper at ICLR 2023
ing the Kolmogorov extension theorem (Oksendal (2014) which produces an interpretation of the stochastic process as a random function f such that p(f) ~ (m(m()), (xx,,))
In the standard GP regression setting, it is assumed that noisy observations of function faz are made such that Vi.n = f + Ei:m where 61:m ~ NN S2) are a sequence of independent and identically distributed Gaussian noise. Following Bayesian approach, prior distribution on the function space is defined such that p(f) ~ GP(O, K(x', x)) where any marginal f; = fx;; has a Gaussian distribution. In general referring to any particular marginal of f is not necessary since GP is defined for all points in X' and the finite distribution is understood from the context. Since the likelihood p( (f)1n.l.) is a product of Gaussians, the posterior distribution over the function space can be analytically found to be a GP with particular mean m(-) and kernel function (-..)) where both functions are defined for any finite set of inputs as shown in (Rasmussen & Williams 2006)). The posterior GP is denoted as p(f(vi:n) GP (m, K).
2.2 TIME-CHANGE
In this section, the classical time-change operation is introduced in one dimension of time and the operation is generalised multidimensional input spaces using subordinated Gaussian fields (Do- brushin [979); Merkle & Barth (2022a b)).
Let (2( be an isotropic stochastic process which has uniformly distributed variance. The operational time can be interpreted as linear representation of change such that the derivative dt is proportional to a deterministic constant. Hence the variance of g(t) scales proportionally to time intervals. Random evolution in variance may be obtained by considering representation of change that is random and nonlinear. Hence define a non-negative, non-decreasing stochastic process (W(1)) such that it randomly maps time instances while preserving their order, therefore changing the time. A time-changed stochastic process {2) is then defined as f(t) = g(W (t))
where the evolution of is governed by dW(1) instead of dt. In other words, the change in will have variance proportional to W (t) (s), instead of - 8 where & > s. Assuming that g(t) is Gaussian, this operation enables large deviations from Gaussian behaviour to occur when dW (t) large, while retaining conditionally Gaussian form.
The random evolution of W W(t) can be modelled as a subordinator that take values in 00, oc) such that it has independent and stationary increments with no fixed discontinuities (Feller (1966); |ertoon
(1997)) ) Thus, subordinator increases non-linearly with certain statistical distribution defined by the random number of discontinuities and their random magnitudes. A Lévy process W (t) in [0, oo) having no drift or Brownian motion defined through its characteristic function E [exp(iu (f)) =
(Kllenberg 2002), Corollary 15.8) where Q is a Lévy measure that satisfies fco..o) (1 A (dx) < oo (Bertoin (197), , p.72). The Lévy measure Q is defined on the random magnitudes of discontinuities, called jumps, and denotes the expected number of jumps per unit time whose magnitudes belong to some subset of the jump space (Tankov, P. and Cont, R. (015)).
exp [Samm (etww ()d())))
By the Lév decomposition, pure jump Lévy process (i.e. containing no Brownian motion) may be expressed using a stochastic integral as
W(t)= = / wN([0, t). dw.
(0,)0)
where N is a bivariate point process having mean measure Leb x Q on [0, x (c) which can be conveniently expressed using Poisson random measure as
N =
Sv.
where (Vi E [], are i.i.d. uniform random variables which give the times of arrival of jumps, (M.) are the sizes of the jumps and Sy. is Dirac measure centered at time VÚ2 and jump size M,. Substituting N into Eq. (1) leads to representation of Lévy jump process infinite series
3

---

Under review as a conference paper at ICLR 2023
w(t) Σ 1=1
M,Irgee a.s.
The almost sure convergence of this series to {W())} is proved in (Rosińiski 2001 Therefore, by sampling pairs of jump times and sizes [V, Mc a realisation of a Lévy process (t)) may be obtained.
The standard formulation of the time-change operation on [0, oo) can be extended to d-dimensional input spaces X' by considering the Poisson random measure representation N of Lévy process. A homogeneous Poisson process expressing the jump times can be generalised to any number di- mensions where we define arbitrary inputs x E R2d (Kingman, J.F.C. (1922). ). The independence properties of Lévy measure allow the definition on the unit time interval to be extended to unit d- dimensional volumes by appropriately scaling the rate of the process (Wolpert & lckstadt (T998b)) For multidimensional input transformations a subordination field on 1º is a d-dimensional stochas- tic process such that each of its dimensions is a subordinator. Thus the i-th dimension of an input vector (+() is mapped to (()) where (+( denotes the subordinator on i. Therefore distance d(æ', x) can be randomly transformed as dx'(,,,, W (z)). Hence, the choice of a Lévy measure characterise the distribution of the random distances over the input space. The notation introduced
for the multidimensional treatment of subordination is omitted for brevity in the following sections as it straightforward extend the model into multidimensional input spaces.
2.3 NON-GAUSSIAN PROCESSES
non- non-nonnnennn Gaussian process (NGP) prior on functions can be obtained by randomly transforming the inputs using subordinator and carrying out GP regression on the transformed input space. The resulting posterior distribution follows a non-Gaussian stochastic process. Given a set of input- output pairs (I,, 34 consider latent input transformation such that z; is mapped to W(r,) where (W x E X is a subordinator. The associated prior on the transformation function defined p(W) and the conditional prior over f is p(flW) ~ GP(mw (x), Kw (z', (x',x)) xx) where mw (x) = mW (x)),
then
(x',z) = ((W(x'), W(x) = ([[W(x') - Wx))) and ((-.-) .) is stationary kernel function e.g. squared exponential Matérn. The joint distribution over the product space of f and W, P(f. W y):)) characterises the NGP prior.
The conditional GP structure of NGP induces posterior mean mw(-) and kernel function (..) that can be evaluated analytically, i.e. p(flii,, W) ~ GPP(mw Kw). ). The conditional likelihood p(yln//W WW) of particular interest this framework since it is a measure of how well analytically.
the data is represented by the model given random transformation and it can also be evaluated
The NGP posterior distribution over the function space is found as
p(flyim) 1 p(flyums W)p(W ((WWi)))d
where p( ( (y).)) is the posterior distribution of the subordinator process. Inferring p( O(W] (Wy)).)) and hence p( //n))) is analytically intractable, however using approximate inference methods allow for straightforward extensions of the model and fully Bayesian inference discussed in the following.
3 SAMPLING AND INFERENCE
In this section, we review shot-noise simulation methods for simulating Lévy processes based on series representations (Rosiński (200T). We describe novel Metropolis-Hasting---
(MH-in-Gibbs) algorithm (Hastings) (1970); CCib Greenberg (1995)) to obtain samples from the posterior distribution of subordinator and estimate non-Gaussian process posterior (fll)):))
4

---

Under review as a conference paper at ICLR 2023
3.1 SHOT-NOISE SIMULATION METHODS
The jump magnitudes {M, shown in Eq. (3) of a Lévy process cannot be directly simulated because there may be infinite number of jumps in any finite interval. One way obtain approx- imate samples from such an infinite sequence is to consider ordering the jump magnitudes by size and simulating large jumps while ignoring or approximating the residual error as discussed in (Fer- guson & Klass (972); >: Rosinski (2001); Wolpert & Ickstadt ( (1998a b): Godsill & Kındap (2021). Once the ordered jump sizes have been obtained, the corresponding jump positions {Veft= may be simulated independently from a uniform distribution on (I1e. Fub)) where 216, Tub are some lower
and upper bounds, or sequentially in space from a homogeneous Poisson process if preferred.
Consider bivariate point process N' that has the same form as Eq. (2) where the jump magnitudes M, are expressed as the output of function (TT;) where { are the epochs of unit rate Poisson process, i.c. the cumulative sum of exponential random variables with unit rate, independent of [Vi] Similar to the standard inverse CDF method for random variate generation, the upper tail mass of a Lévy measure Q+ (xx Q[xx, 20)) can be inverted to produce jump magnitudes of a subordinator by passing epochs of a homogeneous Poisson process through the inverse Lévy measure Q (+) The corresponding function h(.) Q+ (-) is non -on-enccg thus {h(T,)} is an ordered sequence representing random jump sizes. Note that the epochs of a homogeneous Poisson process are analogous to uniformly distributed random variables in (0, o0) and the mapping theorem states that the resulting points {Vi, (Vi,h(T')) form a Poisson point process N = >i= óv,,K(r.) on (Img Fubb) x (0, oo) (ingmaan, J.F.C. ).92)). N' converges almost surely to N as the {T,} sequence is simulated indefinitely (Rosinski 2001)) and approximations of the point process may be obtained through finite samples.
The explicit evaluation of the inverse tail measure (-) in general is not possible. The Lévy measures considered in this paper possess a density function denoted as Q(x)) such that ((dz) = Q(r)dr. The approach taken in this work is to simulate from a tractable dominating point process No having Lévy measure Qo such that (o(dx)/Q((n) be explicitly evaluated. The resulting jump magnitudes belonging to No are then thinned with probability Q(z)/Qo(a) as in (Lewis & Shedler (1979)) obtain the desired approximate jump magnitudes } of subordinator.
2 1, Va (0,co) oo) for which hi)) can
As motivating example in this paper, we consider tempered stable (TS) processes which are com- monly used in mathematical finance to model stochastic volatility (Carr et al. (2003). We note that our methodology is equally applicable to other subordinator processes for which shot noise sim- ulation methods can be applied (Godsill & Kindap (202T); Godsill et al. ( A TS process exhibits both a-stable and Gaussian trends depending on the distance it travels. For short distances the stable characteristics prevail and the TS process produces larger jumps compared to Gaus- sian process. For longer distances the tempering causes a TS process to produce Gaussian trends (Rosiński (2007)) Thus, TS process is a natural extension to Gaussian processes.
The Lévy density for the subordinator TS process is defined as (Shephard & Barndorff- 2001 : Rosiński 2007))
(xx) Cx -1-11-
(4)
where a E (0, 1) is the tail parameter and / is the tempering parameter. The corresponding tail probability may be calculated in terms of gamma functions but it cannot be analytically inverted and numerical approximations are needed (Imai & Kawai (20TT). Instead, we adopt a thinning approach where the Lévy density is factorised into o-stable subordinator process with Lévy density [994):(Godsil Godsill et al. (019)) and tempering function e -az The tail mass of stable process can be found to be Q0 (z) = Sx- and inverting this
1-a Qo(x) = Cx-1-a (Samorodnitsky & Taqqu 1994);
function produces the simulation function h(7) (22) -1/0 Given points 2; from stable point process with density Qo(x), individually selecting (thinning) points with probability e - results in tempered stable process. The associated sampling algorithm is shown in Alg. (1) for reference.
Alg. generates the jumps that correspond to TS process in (0. 1). Since the jumps of a Lévy process are independent and stationary it is straightforward to adjust the interval. For instance, setting the rate of the underlying Poisson process produced in the second stage of Alg. (T) to the
5

---

Under review as a conference paper at ICLR 2023
Algorithm Generation of the jumps of tempered stable process with Lévy density Qrs(œ)
Cx-1- -- -2 where a is the tail parameter and B is the tempering parameter.
1. Assign Nrs = 0.
2. Generate the epochs of a unit rate Poisson process, (Ty; i 1,2,3, 3...
3. For 2..,....
1/a
Compute z;
With probability 8ur, accept 2; and assign Nrs Nrs U 24.
length of the interval (Iis, Fub) produces the associated TS process. Similarly, for d-dimensional input spaces the jumps on a n--dimensional hypercube can be simulated by setting the rate to the associated volume (Wolpert & lckstadt (1998b))
3.2 APPROXIMATE INFERENCE
Since a stochastic process is defined as an infinite collection of random variables, designing direct sampling methods from the posterior p(W p(Wjp (lyy)n) based batch Monte Carlo methods is difficult task. Instead a more appropriate approach to high dimensional problems is to use a Gibbs sampler which approximates samples from multivariate probability distribution or in this case a stochastic process.
Gibbs sampler approximating samples from p(W [1.))) can be implemented by simulating the associated bivariate random points that define the jump size and position small disjoint intervals == (1)) conditioned on the previous sample points - X' \ (z), and observations. Progressively simulating these points such that the whole input space is covered leads to approximate samples from the target distribution. While Gibbs sampling reduces the complexity of sampling a stochastic process for each small interval, direct sampling from the conditional posterior for each interval is still intractable in general. Thus for each interval a Mctropolis Haatopoliis-Hattnnss algorithm is used yielding MH-within Gibbbs sampling algorithm (Chib & Greenberg (1995). The proposal density for the MCMC sampling procedure p(W, W- which produces new bivariate points (jump sizes and times) some interval conditioned all points in as described in Alg. (2).
Algorithm 2 Simulating sample paths from the proposal density
Given a random length set Nw = (V((*) .r((k)) } and an interval (rz, m1) E X,
1. Simulate {V.",, _M.') with rate jzy - mi using Algorithm
2. Remove all points (V:(*) Mr(() from Nw such that v, v((h) < 37 and add (V:",
Nw Nww U (v/"), M
3. Substitute the points of Nw into Eq. to obtain the proposed sample path W
For each realisation of the subordinator (k),, the conditional likelihood p(31.n (K()) may be used as weight in Markov chain Monte Carlo sampler since it is proportional to the posterior distribu- tion p(W (1.m) and we are proposing from p PWW|WW_, As discussed in Section 2.3 the conditional likelihood p(z (W(k) may analytically found given the values of (*) Then given sample (K) and proposal W the acceptance probability for the proposal
aWW"), w(k) min (1, ((V../W(()))
(5)
p(yin (K))

---

Under review as a conference paper at ICLR 2023
Algorithm 3 MH-within Gibbs sampler for p(W (WWy).)).
1. Initialise W by simulating (V, Mi from the associated bivariate point process using
Alg. !
2. Analytically evaluate mwro, Kwo which define the conditional GP posterior p(f Vi:n, W and the conditional likelihood p(yi: W()).
J 3. For N times, iterate over Tj E X' where Ü T =X
j=1
(a) Using Ty and the points (v(k) 3 associated with W((), sample a proposed sam- ple path W using Alg. 2
(b) Evaluate mwo. K Kuro and p(y)=//W)))
(c) With probability a aWV(), ((k)) the proposal accepted and W(++1) W other- wise reject and ((++1) W(k).
The MH-within--ee Sibths sampling procedure is described in Alg. 3 The resulting samples (W()) are individually associated with conditional GP posterior functions p(flyi W()) that are com-
pletely defined through their mean Mwax and covariance Woo functions. Such collection forms
a Gaussian mixture distribution and the mean and covariance of the corresponding mixture density can be obtained
En mwa = ff
(6)
and
N N
Cov
[woo + (Mwwss mmy (mwro - =
(7)
where N is the number of samples and E Ef) Cov define the posterior mean and covari- ance of the random function f. It is straightforward to obtain the corresponding predictive density ((y"(y by adding the observation noise matrix Ω to each covariance matrix sample Ky Us- ing a constant noise matrix Ω corresponds the assumption that the observation likelihood model is Gaussian (Rasmussen & Williams (2066). This assumption can be relaxed by sampling a noise matrix (2) for each individual sample to consider non-Gaussian likelihood models such scale mixture of normals which includes the Student-t and Laplace distributions (Barndorff-Nielsee et al.) (1822). This results in doubly non-Gaussian behaviour which is highly expressive while retaining interpretation of individual components of the behaviour.
Following similar approach the hyperparameters C, a and 3 of the subordinator process may be included the sampling procedure by considering an appropriate prior distribution over their values. Hence these parameters may be marginalised out using the Monte Carlo procedure, which leaves the same number of kernel parameters that define a standard GP. This approach works successfully and will be reported in future publication. Furthermore, nonparametric kernel may be included in this framework by considering a prior distribution on stationary kernel functions and sampling a kernel function for each proposed sample. Some examples of nonparametric kernel design can be found in (Wilson& & Adams (2013): Tobar et al. (2015); Bruinsma et al. 2022 )
A straightffrward extension of Alg. suming that individual subordinator dimensions z((1) independent priori. The simulation steps defined by Alg. [I] and (2) can be independently applied to each dimension and the other steps remain unchanged, replacing step 3. (b) with the multidinensional GP likelihood.
multidimensional input spaces be achieved by
4 EXPERIMENTAL RESULTS AND DISCUSSION
In this section, we present the results of applying NGP regression to non-Gaussian data sets and compare the results with alternative standard GP regression settings. order to emphasise the differences in non-Gaussian and Gaussian processes, first synthetically generated data set
7

---

Under review as a conference paper at ICLR 2023
regression
dersity ubuovvllrs
goazd
interernees ceefernce searve!
24
03
-1
regression
regression
Samples posterior density (w-
-1
Figure 1: Regression analysis results for NGP and two alternative GP models where the observations are generated from NGP prior.
using tempered stable subordinator. Afterwards, in order to demonstrate the generalisation of NGP to multidimensional input spaces a two dimensional example problem is presented using a toy data set available in TensorFlow (Wickham (2066) on diamond prices.
Any stationary kernel function can be used to define the conditional kernel function in NGP re- gression. Here squared exponential (SH) kernel is selected for the conditional GP. In order to demonstrate the ability of a NGP in identifying local characteristics, a length-scale = 0.1 used and an observation set defined small region of the input space (0, 0.5) simulated from NGP prior. The latent transformation space in this example is generated as a TS subordinator with a = 0.8 and [ 5 and the observations have i.i.d. noise with standard deviation 0.1. The observations and results are shown in Fig. ①I.
The NGP predictive posterior density shown in the first column of Fig. (T) clearly identifies some local changes in the variance that match with large jumps observed in the latent transformation pos- terior samples above. Furthermore, for input differences on the order of 1 the density retains close to Gaussian behaviour. The posterior sample paths on the latent transformation space identify how fast the SE kernel decays to zero around different regions X. Large jump sizes break the correlation between local points and the associated observations are treated as statistically independent. From this perspective, if the model correctly identifies large jumps this shows that it discovers some obser- vations contain more information about their local region than stationary GP can encode. Hence, the sample paths provide an expressive probabilistic layer for interpreting non-Gaussian behaviour.
Two alternative GP regression results presented in the second column of Fig. 1 Firstly, using the ground truth value of = 0. a smooth approximation of the predictive posterior density can be obtained. Alternatively, a GP with a higher marginal likelihood can be obtained by optimising the length-scale. In order to adapt to large deviations in some regions of X' the optimisation results in smaller value of / 0.01. This can be considered as a trade-off between in-sample performance and the generalisation capability of the model. As gets smaller each observation is modelled as almost statistically independent and there are no long distance dependencies between inputs. Such a representation will likely be an overfitted model that does not have any generalisation capability NGP regression provides a more sparse representation
and interpolations will produce white-nois of the random function in the sense that it only defines local statistical independence assumptions if the observations show non-Gon non-Gaussian behaviour and otherwise retain long distance dependencies.
For the multidimensional case example regression model for diamond prices is presented. The features used in this task are the carat of a diamond which is a measure of its weight and the per- centage length of its table which is the largest flat facet of the diamond and affects how the diamond interacts with light. For case of visualisation, both input dimensions are linearly transformed to lie between [0, 1]. The experiment is designed such that 1000 randomly selected input-output pairs
8

---

Under review as a conference paper at ICLR 2023
Optimised GP surface (Angle 1)
Optimised GP surface (Anglo 2)
Carat space
Teadoorman manmma vallas
1943
7180
11 Garize
carat
tabi
Nomalizzdd values
NGP surface (Angle 1)
NGP surface (Angle 2)
Table space
Torsformed volues
2500
50
carat
Normalized values
Figure 2: Regression analysis results for NGP and GP models with for the diamond price data set using a TS subordinator. The posterior means are plotted as a surface and the =3 standard deviation surface is overlaid the mean wireframe plot. The right hand column shows the posterior subordinator samples for a TS subordinator.
chosen the training set and the learned posterior surface is compared against another randomly selected 1000 pairs.
Figure 2 show the results of NGP regression using TS subordinator. The main distinguishing property of the diamond price data set is that the price increases non-linearly with increases in the carat feature. The GP surface shown here forms smooth function that cannot model the rapid change in prices around these regions of non-linear increase. Furthermore, the predictive surface making non-zero predictions around out-of-sample regions which are undesirable in downstream decision-making tasks. The NGP surface is able to identify the non-linear increase in prices and the increased uncertainty prices for larger carat values. These properties can be most clearly identified in the posterior subordinator samples. Note that the mean log conditional likelihood of the MCMC samples are found as -97873.9 and the GP log likelihood is found to be 113731.8.
NGP regression with tempered stable subordinator presented in this work may be applied to datasets where there are local deviations from Gaussian behaviour tion closely follows a GP. Using alternative characterisations of the subordinator, varying degrees of non-Gaussian behaviour can be modelled in a NGP regression framework as briefly demonstrated in Section 4 Some practical examples of subordinator processes gamma processes (Rosiński (2001) and inverse Gaussian processes (Rydberg) (1997): (Barndoff--Nelkee
overall trend of the func-
(199)) which both lead to analytical probability density functions unlike the TS process. This fact can be utilised to design better proposal densities for a MCMC procedure. However given any Lévy density a similar formulation Section 3.1 can be readily designed and used for inference studied in Section 3.2
A particularly interesting case is the generalised anverse-Gaussie process which can capture vari- ous degrees of semi- hemi-heevy-ttileed behaviour, including the gamma and inverse Gaussian processes (Eberlein & Hammerstein (2004) : Godsill & Kindap/ (2021)).
Using NGP regression produces a generative model conditioned on dataset where samples from both the posterior density over the input transformation and the predictive density over the observa- tions can be generated. A probabilistic representation a latent layer between the input and output spaces may lead to new insights about the underlying data generating mechanism. Furthermore, our construction of a NGP using the time-change operation may potentially be extended to any prob- abilistic setting for interpreting non non-Gounseee Gaussian behaviour. For example, Lévy fields be used to model the first layer deep GP architecture (Titsias & Lawrence (2000) Damianou & Lawrence 2013)

---

Under review as a conference paper at ICLR 2023
REFERENCES
Ryan Prescott Adams and Oliver Stegle. Gaussian Process Product Models for Nonparametric Nonstationarity. In Proceedings of the 25th International Conference on Machine Learning, ICML '08, pp. 1-8, New York, NY, USA, 2008. Association for Computing Machinery. ISBN 9781605582054. doi: 10.1145/139015..1 139111 0350156.139013
URL https:://cc..rgg///// 11445
O. Barndorff-Nielsen, Kent, and M. Sørensen. Normal Variance-Mean Mixtures and z Distribu- tions. International Statistical Review / Revue Internationate de Statistique, 5012))145-1555 1982.
ISSN 03067734, 17515823. URL/ttpp: ltttp://wet /wwww cestor.orge www.orrgggeeee022222222
O. E. Barndorfi-Nielsen. Processes of normal inverse Gaussian type. Stochastics, 2(1):41-68, 1997. URL https://ideas repec. a spr/ finsto/
Finance and
v2y1997i1p41 hhmmml
Bertoin. Lévy Processes. Cambridge Tracts Mathematics, 121. Cambridge University Press, 1997.
Wessel Bruinsma, Martin Tegnér, and Richard E. Turner. Modelling Non-Smooth Signals With Complex Spectral Structure. In Proceedings of the 25th International Conference on Artificial Intelligence and Statisticx, Proceedings Machine Learning Research. PMLR, 2022.
Peter Carr, Hélyette Geman, Dilip B. Madan, and Marc Yor. Stochastic Volatility for Lévy Processes. Mathematical Finance, 13(3):345-322 2003. URL /ttpttpsss//ccooPaapee.. repec. org/
RePEc: :ta.::::baaaaaaaaa.a..... 11
Siddhartha Chib and Edward Greenberg. Understanding the Metropolis-Hasting Algorithm. The
American Statisticiun, 49(4):327-335, 1995. ISSN 00031305. URL http::/wwet jstor.
brg/stable//6655
Andreas Damianou and Neil D. Lawrence. Deep Gaussian processes. In Carlos M. Carvalho and Pradeep Ravikumar (eds.), Proceedings of the Sixteenth International Conference on Artificial telligence and Statistics, volume 31 of Proceedings of Machine Learning Research, pp. 207-215, Scottsdale, Arizona, USA, 29 Apr-01 May 2013. PMLR. URL https: / /tttpp://proceedinggs mlr. wtt..ppres/vvl//aiaaaaaa.
3a .html
R. L. Dobrushin. Gaussian and their Subordinated Self-similar Random Generalized Fields. The Annals of Probabiliiy, 7(1):1 28, 1979. doi: 10. 10.1214/a 176995115. URL htttps:///ol
atq/10.1214/ /11.00/1122///////799944414
E Eberlein and E. A. Hammerstein. Generalized hyperbolic and inverse Gaussian distributions: Limiting cases and approximation of processes. In Robert C. Dalang, Marco Dozzi, and Francesco Russo (eds.), Seminar Stochastic Analysis, Random Fields and Applications IV, pp. 221-264, Basel, 2004. Birkhäuser Basel. ISBN 978-3-0348-7943.4
W. Feller. An introduction to probability theory and its applications. Number v. 2 in Wiley mathe- matical statistics series. Wiley, 1966.
T. S. Ferguson and M. J. Klass. A Representation of Independent Increment Processes without Gaussian Components. The Annals of Mathematical Statistics, 43(5):1634 - 1643, 1972. doi: 10. /
URL htth:tts: http:://dooo.ooooo org/10 1214 htt.://1111....o........ /tt.a//mmm//171955355
Eric Ghysels, Andrew C. Harvey, and Eric Renault. Stochastic volatility. In Statistical Meth- ods in Finance, volume 14 of Handbook of Statistics, pp. 119-191. Elsevier, 1996. doi:
URL https: http://wwwww. www. wcwencedirectttttttccceeeo
Simon Godsill and Yaman Kindap. Point process simulation of generalised inverse Gaussian pro- cesses and estimation of the Jaeger integral. Statistics and Computing, 32(1):13, Dec 2021. ISSN 1573- doi: 10. 10010077:11222-011110 1212222021110072--1
URL attps: wttp::....r/...................0
10

---

Under review as a conference paper at ICLR 2023
Simon Godsill, Marina Riabiz, and Ioannis Kontoyiannis. The Lévy State Space Model. In 2019 53nd Asilomar Conference on Signals, Systems, and Computers, pp. 487-494 2019. doi: 10. (10//EEECCONF-46440009999888881
W. K. Hastings. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1):97-109, 1970. doi: 10.1099/biommet55...99 URL http::/bi bimmtt.
org//giggggoooo.oooooo exford htt.:////ggggeeeee/eaaaaaeaattee
Junichi Imai and Reiichiro Kawai. On finite truncation of infinite shot noise series representation of tempered stable laws. Physica A A-sttisstical Mechanics and Its Applications, 390:4411-4425, 2011.
Phillip A Jang, Andrew Loeb, Matthew Davidow, and Andrew G Wilson. Scalable Levy Process Priors for Spectral Kernel Learning. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fer- gus, S. Vishwanathan, and Garnett (eds.), Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL https:://prcccedggggg nwwrips....... ca/
O. Kallenberg. Foundations of Modern Probability. Springer- Verlag, 2nd edition, 2002.
Kingman, J.F.C. Poisson Processes. Oxford Studies in Probability. Clarendon Press, 1992.
Miguel Lázaro-Gredilla. Bayesian Warped Gaussian Processes. In Pereira, C.J. Burges, L. Bot- tou, and K.Q. Weinberger (eds.), Advances in Neural Information Processing Systems, vol- ume 25. Curran Associates, Inc., 2012. URL https:://ppoocett ngs neurips..... cc/
A. W Lewis and G. S. Shedler. Simulation of nonhomogeneous Poisson processes by thinning. Naval Research Logistics Quarterly, 26(3):403-413 September 1979. doi: 10. 00.10022/aa.388022033444 URL https://ideas. wwwpe.er....eeeeeeeee. htt.://a.....ee...o..///aa/aaaaa v26y197913p403-413 hhmml
Miguel Lázaro-Gredilla and Michalis Titsias. Variational Heteroscedastic Gaussian Process Regres- sion. pp. 841-848, 01 2011.
David J. C. MacKay. Information Theory, Inference, and Learning Algorithms. Copyright Cam- bridge University Press, 2003.
Robin Merkle and Andrea Barth. On Some Distributional Properties of Subordinated Gaussian Random Fields. Methodology and Computing in Applied Probability, apr 2022a. doi: 10.1007/ si (009-002200958.
Robin Merkle and Andrea Barth. Subordinated Gaussian random fields in elliptic partial differential equations. Stochastics and Partial Differential Equations: Analysis and Computations, pp. 1-49, 2022b.
Bernt Øksendal. Stochastic Differential Equations: An Introduction with Applications. Springer, 6th edition, January 2014.
C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. MIT Press, 2006.
Carl Rasmussen and Zoubin Ghahramani. Infinite Mixtures of Gaussian Process Experts. T. Di- etterich, S. Becker, and Z. Ghahramani (eds.), Advances in Neural Information Processing Sys- tems, volume 14. MIT Press, 2001. URL htt paattps://procceenngg (ttG://eee.acag32eeecca222222223
wewripsssss wwweeuuisss...pppppeeeeeeeeeeeee cwppeper///ppppppppppppocp.e .pdf
Rosiński. Series Representations of Lévy Processes from the Perspective of Point Processes, pp. 401-415. Birkhäuser Boston, Boston, MA, 2001. ISBN 978-1-4612-0117-77 doi: 10.1007/ 978- 178-1-4612-0197-7 18. URL /tt10//00/000////////////////000
18

---

Under review as a conference paper at ICLR 2023
Jan Rosiñski. Tempering stable processes. Stochastic Processes and their Applications, 117(6):
677-707, 2007. ISSN 0304-4149. doi: http://ddbiarrr///00000//......0 10.003. URLpettps
www.sciencedirect. wwt.r//ttt.ccccceeeeeeeeeee http://www.....aooooc.o..ooaaaaa
T. H. Rydberg. The normal inverse Gaussian Lévy process: simulation and approxima- tion. Communications Statistics. Stochastic Models, 13(4):887-910, 1997. doi: 10.1080/ 15326349708807456. URL
ooo.ooo.o/oooooooooooooooooooo.. 10. htt.://///...110132333333333031 (Tt..000000013333777077777777770
Y-LK Samo and Stephen J Roberts. String and membrane Gaussian processes. Journal of Machine Learning Research, 17, 2016.
G. Samorodnitsky and M. S. Taqqu. Stable non-Gaussian random processes stochastic models with infinite variance. CRC Press, 1994.
Alexandra M. Schmidt and Anthony O' O'Hagan. Bayesian Inference for Non-Stationary Spatial Covariance Structure via Spatial Deformations. Journal of the Royal Statistical Society. Se- ISSN 13697412, 14679868. URL
ries B (Statistical Methodology), 65(3):743 758, 2003.
htttp:/wwat
/ 3647549
Amar Shah, Andrew Wilson, and Zoubin Ghahramani. Student-t Processes as Alternatives Gaus- sian Processes. In Samuel Kaski and Jukka Corander (eds.), Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, volume 33 of Proceedings of Machine Learning Research, pp. paress/v333 prrss//333//haha44444 .htmll
Reykjavik, Iceland, 22-25 Apr 2014. PMLR. URL
htttps:/prroeeddingg ml
N. Shephard and O. E. Barndorf-NNicsee. Normal Modified Stable Processes. Economies Series Working Papers 72, University of Oxford, Department of Economics, July 2001. URL https:
//ideas. repec
html
Edward Snelson, Zoubin Ghahramani, and Carl Rasmussen. Warped Gaussian Processes. In S. Thrun, L. Saul, and B. Schölkopf (eds.), Advances in Neural Information Processing Sys- tems, volume 16. MIT Press, 2003. URL 2003/ fi 1e/ 6b57 313355540000 01533447778811111100000000131110
www..enurrppp cc/paper/ pddf
Jasper Snoek, Kevin Swersky, Rich Zemel, and Ryan Adams. Input Warping for Bayesian Op- timization of Non-Stationary Functions. In Eric Xing and Tony Jebara (eds.), Proceed- ings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pp. 1674-1688 Bejing, China, 24-24 Jun 2014. PMLR. URL mlr. press / v32 snoek14 html
https::/prrceeet
Tankov, P. and Cont, R. Financial Modelling with Jump Processes, Second Edition. Chapman and Hall/CRC Financial Mathematics Series. Taylor Francis, 2015.
Michalis Titsias and Neil D. Lawrence. Bayesian Gaussian Process Latent Variable Model. In Yee Whye Teh and Mike Titterington (eds.). Proceedings of the Thirleenth International Con- ference on Artificial Intelligence and Statisties, volume 9 of Proceedings of Machine Learning Research. pp. 844-851, Chia Laguna Resort, Sardinia, Italy, 13-15 May 2010. PMLR. URL https://proceedine
mlr pttp///eseee/////t/ttitiiita.ae. 10t...aa..mmmmmmmm .htll
Felipe Tobar, Thang D Bui, and Richard E Turner. Learning Stationary Time Series using Gaussian Processes with Nonparametric Kernels. In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 28. Curran As- sociates, Inc., 2015. URL ht http:ppppp/pproceeaiigs....oosst hwwp::/wwc................oo.oo. httpppps...pppppp.eeeeeeee cc/pppper/01///ilee 955683440a a 8113
pdf
Michael Unser and Pouya D. Tafti. An Introduction to Sparse Stochastic Processes. Cambridge University Press, 2014. doi: 10. 10.1017//B09781 1077415800.
AED Veraart and M Winkel. Time change. In Fncyclopediia of Quantitative Finance, pp. 1812- 1816. Wiley, 2010. doi: 10.1000/9788047000660220010/1110 10. 100.000//00 0788470000016000 19026
URL httt:://dkd doi org/
Peter Whittle, Stochastic-processes in several dimensions. Bulletin of the International Statistical Institute, 40(2):974-994, 1963.
12

---

Under review as a conference paper at ICLR 2023
Hadley Wickham, g8plo12: Elegant Graphics for Data Analysis. Springer- Verlag New York, 2016.
ISBN 978-3-319-24277-4 URL
www.ttrrtttiiyyiesee.ee.e. orr
Andrew Wilson and Ryan Adams. Gaussian Process Kernels for Pattern Discovery and Ex- trapolation. In Sanjoy Dasgupta and David McAllester (eds.), Proceedings of the 30th Inter- national Conference Machine Learning, volume 28 of Proceedings of Machine Learning Research, pp. 1067-1075, Atlanta, Georgia, USA, 17-19 Jun 2013. PMLR. URL https: //proceed ml p.r.sspprssss.m wataas//22www11o//
html
Andrew Gordon Wilson and Zoubin Ghahramani. Generalised Wishart Processes, 2011. URL https://arxiv: org htcpgg/bb//11 (01110010010111
R.I. Wolpert and K. Ickstadt. Poisson/gamma random field models for spatial statistics. Biometrika,
85(2)251-267, 06 1998a. ISSN 0006-3444. doi: 00.1099/biomm//852.255.. URL
https:
Robert L. Wolpert and Katja Ickstadt. Simulation of Lévy Random Fields, pp. 227-242. Springer
New York, New York, NY, 1998b. ISBN 978-1-4412-1172--999 doi: 10. 10.10079988-14611-173322
A APPENDIX
A.1 ADDITIONAL DETAILS OF THE EXPERIMENTAL SETUP
For the first experimental setup, Gibbs sampler defines grid of 100 disjoint intervals and iterates over the whole space 50 times. After initial samples are discarded for burn-in the average log conditional likelihood of the remaining samples is found to be 51.46 with standard deviation of 3.17. In comparison the log conditional likelihood of the data generating process is found as 59.09 which suggests that the NGP does not overfit. The log marginal likelihood of the GPs with = 0.1 and 0.01 are found as -88.339 and -71.38, respectively. The mean acceptance probability for each step in Alg. 3 is found as 0.72. The confidence intervals in regression results show standard deviations. Lastly, on the right column we show 5 samples from the NGP posterior density and its smooth approximation.
Given the results above the two main aspects that require further attention are the design of sensible prior distribution on the latent transformation space and the initialisation of the Gibbs sampler. The tempered stable (TS) subordinator is characterised by three parameters, the tail heaviness, tempering and scale. The expected value and variance of the subordinator process on an input space X' is a function of these parameters and the length (or measure) of X. In order to produce results that are comparable with Gaussian process (GP) regression, the expected value of the subordinator process is set to the length of X, i.e. given value is [mmaa - Tminl As discussed in Section 2.2 if the observed input points are assumed to lie on Euclidean space, the change in the covariance scales linearly according to [Timar Imin| On the latent transformation space, this corresponds to an identity map in (Zmin Tmaz). Setting the expected value to the length of the input set expresses a preference towards regular Gaussian behaviour and the deviations from the identity map provide insight into the characteristics of the observed data set.
B and C, that represent
interval (#min: Zmax) the expected
Note again that the shot-noise simulation methods studied in Section 3.1 produce approximate sam- ple paths from a Lévy process since the infinite series described by Eq. (3) have to be truncated to a finite number of terms. The convergence of the series is found in practice to be faster for smaller values of a and 9 parameters. The number of terms that are required obtain a sufficiently close es- timate can be adaptively found using probabilistic asymptotic bounds and will be presented in future extensions of our work. However, a simple way to ensure the required convergence is to increase the number of terms produced /3 gets larger. In practice for 3 = 5, producing 1000 terms in Alg. is found to work well.
The initialisation of the Gibbs sampler in practice is of the central issues in designing Markov chain Monte Carlo methods. In Alg. B the sampling method depends on the number and size of small disjoint intervals т. As the size of T is decreased, the acceptance rates for the Gibbs
13

---

Under review as a conference paper at ICLR 2023
sampler increase and convergence be achieved in few number of iterations of the whole input space. However, this also results in an increased number of intervals and each iteration requires more computation and time. Our strategy to initialise the Gibbs sampler is to obtain crude estimate of the latent transformation function by running the same algorithm on a few number of intervals т and set the initial state of the Gibbs sampler the jump magnitudes and positions that correspond the last sample path. The chain can then be simulated starting from this state and only a few samples have be discarded as burn-in. Alternatively, a simpler initialisation is to generate linearly spaced points with equal magnitudes that represent the identity map corresponding to GP regression setting. This alternative initialisation can be particularly useful for dataset that display close to Gaussian behaviour.
0
Figure 3: Regression analysis results for NGP and GP models with the Matérn kernel where the observations are generated from a NGP prior.
It is straightforward to change the particular choice of kernel function without any changes in our presented algorithms. In Fig. 3 an example regression problem presented using the conditional kernel function specified as the Matérn kernel. The Matérn kernel decays slower than squared exponential kernel and therefore can be used to model long range dependencies, Particular parameter values of the Matérn kernel correspond to the covariance of certain stochastic differential equations as described in Whittle (1963). The results are again compared with a GP regression setting to highlight the differences. The smoothness parameter is chosen
5/2 and the length-scale is I = 0.1. The subordinator parameters are identical to the previous example in Fig. [ A dataset of size 500 is generated (0, 1) from the NGP prior and set of size 100 randomly selected as the observations. In this case, the unobserved elements in the dataset serve as the test set. For the NGP regression, most unobserved data points lie inside the confidence intervals while the GP regression misses these points around local deviations. The log conditional likelihood of the NGP and the marginal likelihood of the GP models are 18.25 and 114.997 respectively.
Additionally, in order to emphasise that our framework works for any choice of Lévy process with measure that possesses a density function, the multidimensional experiment is run separately using gamma subordinator: The required simulation algorithm for the gamma subordinator is presented in Rosiński (200T); Godsill & Kindap (202T). The training set is independently sampled and the results are presented Figures and 5 Similar conclusions can be made using the gamma subordinator. However note that the gamma subordinator is characterised by the flat intervals shown in Fig. and there may be other applications that are more appropriate for such posterior surface.
14

---

Under review as a conference paper at ICLR 2023
Optimised GP surface (Angle 1)
Optimised GP surface (Angle
observed points
unobserved points
20020
20000
2500
1900
13002
15020
s00
10002
Inoce
15000
10020
3500
5030
5003
carat
0,6 0.
9.
carat
table
NGP surface (Angle 1)
NGP surface (Angle 2)
37502
17520
15002
15030
15000
10200
s000
1500
1500
5500
5003
193
PA Carat
04 92
carat
table
Figure 4: Regression analysis results for NGP and GP models with for the diamond price data using a gamma subordinator. The posterior means are plotted as a surface and the +-3 standard deviation surface is overlaid on the mean wireframe plot.
Carat space
Table space
Tarsformed values
Parnformeed valuees
Normllsed values
Normaized values
Figure 5: Posterior subordinator samples for gamma subordinator.
15


Reviews:
[{"review_id": "zaAZjjNPNF", "reviewer": "Reviewer_ebws", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "Mechanisms for MCMC inference in NGP are introduced. NGPs are defined GPs in which the input has been warped according to a monotonically non-decreasing function with a pure jump Lévy process prior.", "summary_of_the_review": "An interesting theoretical derivation, but this paper looks more like a solution in search of a problem. The experiments are unconvincing for being in lower dimension and not try hard enough to make the standard GP look good (measuring the MSE for the best hyperparameter selection -- instead measuring only the log-likelihood for what looks like an inadequate lengthscale/kernel selection).", "strength_and_weaknesses": "Strengths:\n- Clear and well written\n- Technically correct\n- Effort in the development of an efficient, non-trivial MCMC\n\nWeaknesses:\n\n- Only tested in 1D and 2D regression\n- The hyperparameters of the NGP aren't tuned with this procedure\n- The 1D problem is a toy example in which the data follows exactly the model\nIn the 2D problem it is unknown what MSE is, and how it is affected by the choice of hyperparameters. The standard GP seems oversmoothed. If the lengthscales are tuned, is it not possible to get a similar MSE to that of the NGP? This very fundamental question is left unanswered. Given the dimensionality of the data, both lengthscales could be swept to properly check this conclusively. It says \"optimized GP\", so I'm assuming that some lengthscale selection is happening, but his might be falling in a local minimum. The surface is clearly oversmoothed. Also, other kernels should be tested to attempt best match to the data while still remaining within the tractability of the GP.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "mRy3Jrxvmc", "reviewer": "Reviewer_AecU", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "In the manuscript, \"Non-Gaussian Process Regression\", the authors present a Bayesian modelling strategy for semi-parametric regression in which the standard setting of Gaussian process regression is augmented with a transformation of the input space through a 'time change' operation.  In particular, the latter is based on a Levy process to allow the regression function to better follow apparent discontinuities than can the original Gaussian process. ", "summary_of_the_review": "Overall, I think the authors have the start of a useful piece of work here (namely, the formulation and implementation of an interesting model), but in the present exposition I do not find a compelling case for becoming interested in this model over any other type of 'fancy' Gaussian process.", "strength_and_weaknesses": "The construction of the model and the strategy described for its computational implementation (i.e., approximate posterior simulation) appear to be sound and well thought out.  Yet the motivation for developing this type of model and its inferential strengths and weaknesses over other approaches are not established to my satisfaction.  \n\nIt is certainly the case that non-smoothness is evident in many real world datasets, and at times this may be such a magnitude that the smooth functional form preferred (or even 'demanded', according to the nature of the kernel) under the Gaussian process prior becomes a hinderance to effective inference.  The consequences of this model misspecification and the pathway to appropriate model revisions is highly context dependant.  In a general sense this means that the objectives and conditions of the inference must be specified clearly---ideally through an appropriate utility function---and the nature of the robustness sought must be established. E.g. is the exercise one of smoothing over noisy observations along a well-sampled but noisily observed window of data with the objective of minimising mean squared error relative to the true latent function, or is the focus rather on interpolation over a missing segment of the curve with the utility being a Frequentist coverage performance of the pointwise posterior bounds?  Only once these aims are established can one then effectively investigate the performance of the proposed solution, both in absolute terms and relative to competing methods.  This investigation can reasonably be either via computational experiments as in the present manuscript or via theoretical study; with ideally the two aspects of the literature complementing each other.  In the former (computational experiments) strength is given by the demonstration of substantive computational and/or accuracy superiority over existing methods in real world problems, or by demonstration on mock datasets that the new method solves a meaningful failure mode of the existing solution.\n\nIn the present set of experiments it seems that the exercise is one of smoothing noisy data (more so than interpolation) with discontinuities such a those drawn from the proposed NGP process itself, and the utility lies in the reduced degree of oversmoothing that results compared with the base Gaussian process.  In this setting it would seem that the natural comparison points are the likes of Gramacy & Lee's Bayesian treed GP model (JASA 2008) and Kim, Mallick & Holmes' piecewise Gaussian processes (JASA 2005).  Models in which the mean function of the GP is treated as a piecewise linear process would also be a natural alternative: why does the discontinuity need to be put through the kernel?  Indeed, one also thinks of Frequentist non-parametric regression methods (e.g. piecewise linear regression) here too: what is the Bayesian aspect of this problem adding from an inferential point of view (because it's unlikely to be competitive in speed with the need to do full posterior sampling)?  How is the hyper-prior specification complicated by the 'degeneracy' between the Levy process parameterisation and the kernel parameterisation?\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "Lh_aytKAV3", "reviewer": "Reviewer_bFEN", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposed a non-Gaussian process regression which allows to model GPs with time changes and thus accounts for non-Gaussian behaviours such as heavy tails. The model is constructed using a latent transformed input space, and the Levy process is used to model the random evolution of the latent transformation. Further, MCMC (MH-within-Gibbs sampler) is used for the inference. ", "summary_of_the_review": "The paper proposes an interesting method for modelling potentially non-Gaussian data with heavy tales. However, the current evaluation is rather limited and does not fully demonstrate the usefulness of the method compared to standard GP regression with Matern kernel and GP regression with change point kernel with Matern kernels as base kernels. From the presented examples, some limitations of the method also cannot be identified, i.e. how many jumps are feasible in terms of the inference, how robust MCMC will be when there are more jumps in the process (for example, does initializing Gibbs sampler become harder?). ", "strength_and_weaknesses": "Strengths:\n-- The proposed method is interesting and promising for cases when one wants to use a GP, but the data are possibly non-Gaussian with non-stationary jumps. \n-- The paper is overall clear and well presented. \n\nWeaknesses: \n-- A more extensive comparison with possible competitors is missing.\n-- The examples demonstrating the performance of the method are quite limited. In particular, it would be nice to see examples with longer data sets containing more jumps. \n-- The limitations of the method (computational and conceptual are not completely clear; how many dimensions are feasible; how many jumps are identifiable; how robust is the method to initialization of the Gibbs sampler etc.). \n\nDetailed comments:\n•\tIs the MCMC algorithm feasible in higher dimensions? Can sample paths in different dimensions be dependent?\n•\tThere are only 2 points with jumps in Figure 1; what would happen with a longer data series with more jumps (say above 100)? Is that still feasible to identify them? Would the inference be possible in this case, and what effect would initialization in Gibbs sampler have here? \n•\tIf the interest is in the location of the jumps, what would be the advantage of NGP compared to a GP regression and, say, change point kernel with Matern as base kernels?\n•\tIf we assume some non-stationarity, how the fitted results would be different between NGP and GP with a Matern kernel? What would be the advantage of using NGP over a standard GP regression with the Matern kernel? I further noticed that this kind of experiment is presented in Appendix A1, but the advantages are still not completely clear. Why NGP fit is considered to be better?\n•\tIt would be interesting to see how NGP performs out of sample, i.e. what happens with the predictions if you extrapolate GP to the right in Figures 1 and 3. \n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject