A Gaussian Process Regression Model for Distribution Inputs

TL;DR

This paper introduces a new family of Gaussian process models for distribution inputs using transportation-based kernels, enabling efficient forecasting of processes indexed by probability distributions.

Contribution

It develops positive definite kernels based on Wasserstein distances for Gaussian process modeling with distribution inputs, providing a probabilistic framework and forecasting methods.

Findings

New kernels based on Wasserstein distances for distribution inputs

Probabilistic characterization of Gaussian processes with these kernels

Efficient forecasting of distribution-indexed Gaussian processes

Abstract

Monge-Kantorovich distances, otherwise known as Wasserstein distances, have received a growing attention in statistics and machine learning as a powerful discrepancy measure for probability distributions. In this paper, we focus on forecasting a Gaussian process indexed by probability distributions. For this, we provide a family of positive definite kernels built using transportation based distances. We provide a probabilistic understanding of these kernels and characterize the corresponding stochastic processes. We prove that the Gaussian processes indexed by distributions corresponding to these kernels can be efficiently forecast, opening new perspectives in Gaussian process modeling.

Tables

Table 1. TABLE I: Values of different quality criteria for the “distribution”, “Legendre” and “PCA” Gaussian process models. The “distribution” Gaussian process model is based on covariance functions operating directly on the input distributions, while “Legendre” and “PCA” are based on linear projections of the input distributions on finite-dimensional spaces. For “Legendre” and “PCA”, the order value is the dimension of the projection space. The quality criteria are the root mean square error (RMSE) which should be minimal and the confidence interval ratio ( C ​ I ​ R 0.9 𝐶 𝐼 subscript 𝑅 0.9 CIR_{0.9} ) which should be close to 0.9 0.9 0.9 . The “distribution” Gaussian process model clearly outperforms the two other models.

model	RMSE	$CI{R}_{0.9}$
“distribution”	$0.094$	$0.92$
“Legendre” order 5	$0.49$	$0.92$
“Legendre” order 10	$0.34$	$0.89$
“Legendre” order 15	$0.29$	$0.91$
“PCA” order 5	$0.63$	$0.82$
“PCA” order 10	$0.52$	$0.87$
“PCA” order 15	$0.47$	$0.93$

Table 2. TABLE II: Same setting as in Table I , except that the input and predictand distributions are only observed indirectly, through sample values from them. The “distribution” model suggested in this paper clearly outperforms the “kernel regression” procedure.

model	RMSE	$CI{R}_{0.9}$
“distribution”	$0.21$	$0.91$
“kernel regression”	$0.93$

Equations

$${\mathcal{T}}_2(\mu ,\nu )=\underset{\pi \in \Pi (\mu ,\nu )}{\inf }\int {\left|x-y\right|}^{2}d\pi (x,y).$$

$$W_2(\mu ,\nu )={\mathcal{T}}_2(\mu ,\nu {)}^{1/2}.$$

$$F_{\mu }^{-1}(t)=\inf \{u,F_{\mu }(u)\ge t\}$$

$$M:x\mapsto E(X(x))$$

$$K:(x,y)\mapsto \mathrm{Cov}(X(x)X(y))$$

$$\sum _{i,j=1}^n{\lambda }_i{\lambda }_{j}K(x_i,x_j)\ge 0.$$

$$K(x,x)=0$$

$$\sum _{i,j=1}^n{c}_i{c}_{j}K(x_i,x_j)\le 0.$$

$$\mathrm{Cov}(X(g(x)),X(g(y)))=\mathrm{Cov}(X(x),X(y)).$$

$$\mathrm{Cov}\left(X(g(x))-X(g(o))\right)=\mathrm{Cov}\left(X(x)-X(o)\right).$$

$$g_{\#}(\mu ):A\mapsto \mu (g^{-1}(A)).$$

$$\begin{gathered} K^{H,{\mu }_0}(\mu ,\nu ) \\ =\frac{1}{2}\left(W_2^{2H}({\mu }_0,\mu )+W_2^{2H}({\mu }_0,\nu )-W_2^{2H}(\mu ,\nu )\right), \end{gathered}$$

$$\{\begin{array}{rl}EX(\mu )& =0,\\ \mathrm{Cov}(X(\mu ),X(\nu ))& =K^{H,{\mu }_0}(\mu ,\nu )\end{array}$$

$$E(X(\mu )-X(\nu ){)}^2=W_2^{2H}(\mu ,\nu ),$$

$$Y(\mu ):=X(\bar{\mu })+f(m(\mu )),$$

$$\mathrm{Var}(Y(\mu ))=E(X(\bar{\mu }){)}^2=E(X(\bar{\mu })-X({\delta }_0){)}^2$$

$$={\left(\mathrm{Var}(\mu )\right)}^H.$$

$$E(Y(\mu ))=E(F_{\mu }^{-1}(U))$$

$$\mathrm{Cov}(Y(\mu ),Y(\nu ))=\mathrm{Cov}(F_{\mu }^{-1}(U),F_{\mu }^{-1}(U)),$$

$$(\mu ,\nu )\mapsto F\left(W_2^{2H}(\mu ,\nu )\right)$$

$$(-1{)}^n{F}^{(n)}(x)\ge 0.$$

$$F(x)={\int }_{{\mathbb{R}}^{+}}{e}^{-\lambda x}d{\mu }_F(\lambda ).$$

$$e^{-\lambda W_2^{2H}(\mu ,\nu )},$$

$$\left\{K_{{\sigma }^2,\ell ,H}={\sigma }^2{e}^{-\frac{W_2^{2H}}{\ell }},({\sigma }^2,\ell ,H)\in C\times C^{\prime }\times [0,1]\right\},$$

$$(Z(\mu ){)}_{\mu \in {\mathcal{W}}_2(\mathbb{R})}:=(F_{\mu }^{-1}(U){)}_{\mu \in {\mathcal{W}}_2(\mathbb{R})},$$

$$(\mu ,\nu )\mapsto E(Z(\mu )-Z(\nu ){)}^2$$

$${\beta }_E:=\sup \left\{\beta >0|d^{\beta }\text{is negative definite}\right\}.$$

$$\{K_{\theta };\theta \in \Theta \},$$

$$L_{\theta }=\frac{1}{n}\ln (\det R_{\theta })+\frac{1}{n}{y}^t{R}_{\theta }^{-1}y,$$

$${\hat{Y}}_{\theta }(\mu )=r_{\theta }^t(\mu )R_{\theta }^{-1}y$$

$$r_{\theta}(\mu)=\left[\begin{array}[]{c}K_{\theta}(\mu,\mu_{1})\\ \vdots\\ K_{\theta}(\mu,\mu_{n})\end{array}\right].$$

Full text

I Introduction

Originally used in spatial statistics (see for instance [2] and references therein), Kriging has become very popular in many fields such as machine learning or computer experiment, as described in [3]. It consists in predicting the value of a function at some point by a linear combination of observed values at different points. The unknown function is modeled as the realization of a random process, usually Gaussian, and the Kriging forecast can be seen as the posterior mean, leading to the optimal linear unbiased predictor of the random process. Gaussian process models rely on the definition of a covariance function that characterizes the correlations between values of the process at different observation points. As the notion of similarity between data points is crucial, i.e. close location inputs are likely to have similar target values, covariance functions are the key ingredient in using Gaussian processes, since they define nearness or similarity. In order to obtain a satisfying model one need to chose a covariance function (i.e. a positive definite kernel) that respects the structure of the index space of the dataset. Continuity of the covariance is a minimal assumption, as one may ask for additional properties such as stationarity or stationary increments with respect to a distance. These stronger assumptions allow to obtain a model where the correlations between data points depend on the distance between them.

First used in Support Vector (see for instance [4]), positive definite kernels are nowadays used for a wide range of applications. There is a huge statistical literature dealing with the construction and properties of kernel functions over $\mathbb{R}^{d}$ for $d\geq 1$ (we refer for instance to [5] or [6] and references therein). Yet the construction of kernels with adequate properties on more complex spaces is still a growing field of research (see for example [7], [8], [9]).

Within this framework, we tackle the problem of forecasting a process indexed by one-dimensional distributions. Our motivations come from a variety of applied problems: in the classical ecological inference problem (see [10]), outputs are not known for individual inputs but for groups, for which the distribution of a covariate is known. This situation happens for instance in political studies, when one wants to infer the correlation between a vote and variables such as age, gender or wealth level, from the distributions of these covariates in different states (see for example [10]). The problem of causal inference can also be considered in a distribution learning setting (see [11]).

As [12] remarks, learning on distribution inputs offers two important advantages in the big data era. By bagging together individual inputs with similar outputs, one reduces the size of a dataset and anonymizes the data. Doing so results on learning on the distribution of the inputs in the bags.

Another application arises in numerical code experiments, when the prior knowledge of the input conditions may not be an exact value but rather a set of acceptable values that will be modeled using a prior distribution. Hence we observe output values for such probability distributions and want to forecast the process for other ones. A similar application of distribution inputs for numerical code experiments is given by non-negative functional inputs. We give a detailed example of this situation in Section II.

Several approaches already exist to deal with distribution inputs regression. An important class of methods relies on some notion of divergence between distributions (see [13, 14, 15]). Other methods have been proposed, such as kernel mean embedding [12] and kernel ridge regression methods [16]. In this paper we focus on Gaussian process regression method.

The first issue when considering Gaussian process regression for distribution inputs is to define a covariance function, which will allow to compare the similarity between probability distributions. Several approaches can be considered here. The simplest method is to compare a set of parametric features built from the probability distributions, such as the mean or the higher moments. This approach is limited as the effect of such parameters do not take into account the whole shape of the law. Specific kernels should be designed in order to map distributions into a reproducing kernel Hilbert space in which the whole arsenal of kernel methods can be extended to probability measures. This issue has recently been considered in [17] or [18].

In the past few years, transport based distances such as the Monge-Kantorovich or Wasserstein distance have become a growing way to assess similarity between probability measures and are used for numerous applications in learning and forecast problems. Since such distances are defined as a cost to transport one distribution to the other one, they appear to be a very relevant way to measure similarities between probability measures. Details on Wasserstein distances and their links with optimal transport problems can be found in  [19]. Applications in statistics are developed in [20], [21, 22] while kernels have been developed in [18] or [23].

In this paper, we construct covariance functions in order to obtain Gaussian processes indexed by probability measures. We provide a class of covariances which are functions of the Monge-Kantorovich distance, corresponding to stationary Gaussian processes. We also give covariances corresponding to the fractional Brownian processes indexed by probability distributions, which have stationary increments with respect to the Monge-Kantorovich distance. Furthermore we show original nondegeneracy results for these kernels. Then, in this framework, we focus on the selection of a stationary covariance kernel in a parametric model through maximum likelihood. We prove the consistency and asymptotic normality of the covariance parameter estimators. We then consider the Kriging of such Gaussian processes. We prove the asymptotic accuracy of the Kriging prediction under the estimated covariance parameters. In simulations, we show the strong benefit of the studied kernels, compared to more standard kernels operating on finite dimensional projections of the distributions. In addition, we show in the simulations that the Gaussian process model suggested in this article is significantly more accurate that the kernel smoothing based predictor of [24]. Our results consolidate the idea that the Monge-Kantorovich distance is an efficient tool to assess variability between distributions, leading to sharp predictions of the outcome of a Gaussian process with distribution-type inputs.

The paper falls into the following parts. In Section III we recall generalities on the Wasserstein space, covariance kernels and stationarity of Gaussian processes. Section IV is devoted to the construction and analysis of an appropriate kernel for probability measures on $\mathbb{R}$. Asymptotic results on the estimation of the covariance function and properties of the prediction of the associated Gaussian process are presented in Section V. Section VI is devoted to numerical applications while the proofs are postponed to the appendix.

II An applicative case from nuclear safety

The research that conducted to this article have been partially funded by CEA, and is motivated by a nuclear safety application, which we detail here.

A standard problem for used fissile storage process is the axial burn up analysis of fuel pins [25]. In this case study, fuel pins may be seen as one-dimensional curves $X:[0,1]\to\mathbb{R}^{+}$ [26]. These curves correspond to the axial irradiation profiles for fuel in transportation or storage packages which define the neutronic reactivity of the systems. From a curve $X$, corresponding to a given irradiation profile, it is then possible to compute the resulting neutron multiplication factor $k_{eff}(X)$ by numerical simulation [27]. It can be insightful, for profiles with a given total irradiation $\int_{0}^{1}X(t)dt$, to study the impact of the shape of the irradiation curve $X$ on the multiplication factor $k_{eff}(X)$. This type of study can be addressed by considering $k_{eff}$ as a realization of a Gaussian process indexed by one-dimensional distributions.

III Generalities

In this section we recall some basic definitions and properties of the Wasserstein spaces and of covariance kernels.

The Monge-Kantorovich distance

Let us consider the set $\mathcal{W}_{2}(\mathbb{R})$ of probability measures on $\mathbb{R}$ with a finite moment of order two. For two $\mu,\nu$ in $\mathcal{W}_{2}\left(\mathbb{R}\right)$ , we denote by $\Pi(\mu,\nu)$ the set of all probability measures $\pi$ over the product set $\mathbb{R}\times\mathbb{R}$ with first (resp. second) marginal $\mu$ (resp. $\nu$).

The transportation cost with quadratic cost function, or quadratic transportation cost, between these two measures $\mu$ and $\nu$ is defined as

[TABLE]

This transportation cost allows to endow the set $\mathcal{W}_{2}\left(\mathbb{R}\right)$ with a metric by defining the quadratic Monge-Kantorovich, or quadratic Wasserstein distance between $\mu$ and $\nu$ as

[TABLE]

A probability measure $\pi$ in $\Pi(\mu,\nu)$ realizing the infimum in (1) is called an optimal coupling. This vocabulary transfers to a random vector $(X_{1},X_{2})$ with distribution $\pi$. We will call $\mathcal{W}_{2}(\mathbb{R})$ endowed with the distance $W_{2}$ the Wasserstein space.

We will consider on several occasions the collection of random variables $(F^{-1}_{\mu}(U))_{\mu\in\mathcal{W}_{2}(\mathbb{R})}$, where $F^{-1}_{\mu}$ defined as

[TABLE]

denotes the quantile function of the distribution $\mu$, and $U$ is an uniform random variable on $[0,1]$. For every $\mu,\nu\in\mathcal{W}_{2}(\mathbb{R})$, the random vector $((F^{-1}_{\mu}(U)),(F^{-1}_{\nu}(U)))$ is an optimal coupling (see [19]). Notice that the random variable $F^{-1}_{\mu}(U)$ does not depend on $\nu$, so that $(F^{-1}_{\mu}(U))_{\mu\in\mathcal{W}_{2}(\mathbb{R})}$ is an optimal coupling between every distribution of $\mathcal{W}_{2}(\mathbb{R})$.

More details on Wasserstein distances and their links with optimal transport problems can be found in [28] or [19] for instance.

Covariance kernels

Let us recall that the law of a Gaussian random process $(X(x))_{x\in E}$ indexed by a set $E$ is entirely characterized by its mean and covariance functions

[TABLE]

and

[TABLE]

(see $e.g.$ [29]).

A function $K$ is actually the covariance of a random process if and only if it is a positive definite kernel, that is to say for every $x_{1},\cdots,x_{n}\in E$ and $\lambda_{1},\cdots,\lambda_{n}\in\mathbb{R},$

[TABLE]

In this case we say that $K$ is a * covariance kernel*.

On the other hand, any function can be chosen as the mean of a random process. Hence without loss of generality we focus on centered random processes in Section IV.

Positive definite kernels are closely related to negative definite kernels. A function $K:E\times E\rightarrow\mathbb{R}$ is said to be a negative definite kernel if for every $x\in E$,

[TABLE]

and for every $x_{1},\cdots,x_{n}\in E$ and $c_{1},\cdots,c_{n}\in\mathbb{R}$ such that $\sum_{i=1}^{n}c_{i}=0$,

[TABLE]

• Example

The variogram $(x,y)\mapsto\operatorname{\mathbb{E}}(X(x)-X(y))^{2}$ of any random field $X$ is a negative definite kernel.

If the inequality (3) (resp. (5)) is strict as soon as not every $\lambda_{i}$ (resp. $c_{i}$) is null and the $x_{i}$ are two by two distinct, a positive definite (resp. negative definite) kernel is said to be nondegenerate. Nondegeneracy of a covariance kernel is equivalent to the fact that every covariance matrix built with $K$ is invertible. We will say that a Gaussian random process is nondegenerate if its covariance function is a nondegenerate kernel. Nondegeneracy is is usually a desirable condition for Kriging, since the forecast is built using the inverse of the covariance matrix of the observations. In addition, Gaussian process models with degenerate kernels have structural restrictions that can prevent them for being flexible enough. We give the nondegeneracy of the fractional Brownian motion indexed by the Wasserstein space in Section IV.

Stationarity

Stationarity is a property of random processes that is standard in the Kriging literature. Roughly speaking, a stationary random process behaves in the same way at every point of the index space. It is also an enjoyable property for technical reasons. In particular it is a key assumption for the proofs of the properties we give in Section V.

We say that a random process $X$ indexed by a metric space $(E,d)$ is stationary if it has constant mean and for every isometry $g$ of the metric space we have

[TABLE]

Let us notice in particular that if the covariance of a random process is a function of the distance, equation (6) is verified. This is the assumption we make in Section V.

One can also find the assumption of stationarity for the increments of a random process. Many classical random processes have stationarity increments, such as the fractional Brownian motion. We prove the existence of fractional Brownian motion indexed by the Wasserstein space in Section IV.

We will say that $X$ has stationary increments starting in $o\in E$ if $X$ is centred, $X(o)=0$ almost surely, and for every isometry $g$ we have

[TABLE]

Notice that the variance of a random process with stationary increments increases as the input gets far from the origin point $o$.

Let us remark that the definitions we gave are usually called “in the wide sense”, in contrast with stationarity definitions “in the strict sense”, which asks for the law of the process (or its increments) to be invariant under the action of the isometries, and not only the first and second moments. Since we are only dealing with Gaussian processes those definitions coincide.

Isometries of the Wasserstein space

Since we are interested in processes indexed by the Wasserstein space with stationarity features, let us recall a few facts about isometries of the Wasserstein space $\mathcal{W}_{2}(\mathbb{R})$, that is to say maps $i:\mathcal{W}_{2}(\mathbb{R})\rightarrow\mathcal{W}_{2}(\mathbb{R})$ that preserve the Wasserstein distance.

Trivial isometries come from isometries of $\mathbb{R}$: to any isometry $g:\mathbb{R}\rightarrow\mathbb{R}$, we can associate an isometry $g_{\#}:\mathcal{W}_{2}(\mathbb{R})\rightarrow\mathcal{W}_{2}(\mathbb{R})$ that maps any measure $\mu\in\mathcal{W}_{2}(\mathbb{R})$ to the measure

[TABLE]

Stationarity of a random process with regard to these trivial isometries is an interesting feature, since it means that the statistical properties of the outputs do not change when we apply an isometry to the real line.

However let us mention that not every isometry of the Wasserstein space is trivial. In particular, mapping every distribution to its symmetric regarding its expectancy defines an isometry of $\mathcal{W}_{2}(\mathbb{R})$. We refer to [30] for a complete description of the isometries of the Wasserstein space.

IV Gaussian process models for distribution inputs

In this section we give covariance kernels on the space of probability distributions on the real line. This allows for modeling and Gaussian process regression of datasets with distribution inputs.

We start in Section IV-A by giving a generalization of the seminal fractional Brownian motion to distributions inputs endowed with the Wasserstein distance.

Then, in Section IV-B we give Gaussian processes that are stationary with respect to the Wasserstein distance on the inputs.

IV-A Fractional Brownian motion with distribution inputs

We first consider the family of fractional Brownian kernels

[TABLE]

where $0<H\leq 1$ and $\mu_{0}\in\mathcal{W}_{2}(\mathbb{R})$ are fixed.

Note that these kernels are obtained by taking the covariances of the classical fractional Brownian motions and replacing the distance $|t-s|$ between two times $s,t\in\mathbb{R}$ by the Wasserstein distance $W_{2}(\mu,\nu)$ between two distribution inputs. The measure $\mu_{0}\in\mathcal{W}_{2}(\mathbb{R})$ plays the role of the origin $0\in\mathbb{R}$.

Theorem IV.1**.**

For every $0\leq H\leq 1$ and a given $\mu_{0}\in\mathcal{W}_{2}(\mathbb{R})$ the function $K^{H,\mu_{0}}$ defined by (8) is a covariance function on $\mathcal{W}_{2}(\mathbb{R})$. Furthermore $K^{H,\mu_{0}}$ is nondegenerate if and only if $0<H<1$.

The Gaussian process $(X(\mu))_{\mu\in\mathcal{W}_{2}(\mathbb{R})}$ such that

[TABLE]

is the $H$-fractional Brownian motion with index space $\mathcal{W}_{2}(\mathbb{R})$ and origin in $\mu_{0}$. It inherits properties from the classical fractional Brownian motion.

It is easy to check that the output at the origin measure $\mu_{0}$ is zero, $X(\mu_{0})=0\text{ almost surely.}$ Furthermore

[TABLE]

from which we deduce that $(X(\mu))_{\mu\in\mathcal{W}_{2}(\mathbb{R})}$ has stationary increments, which means that the statistical properties of $X(\mu)-X(\nu)$ are the same as those of $X(g(\mu))-X(g(\nu))$ for every isometry $g$ of the Wasserstein space.

The fractional Brownion motion is well known for its parameter $H$ governing the regularity of the trajectories: small values of $H$ correspond to very irregular trajectories while greater values give steadier paths. Moreover for $H>1/2$ the process exhibits long-range dependence (see [31]).

From the modelling point of view, it is interesting to consider the following process: consider $(X(\mu))_{\mu\in\mathcal{W}_{2}(\mathbb{R})}$ the $H$-fractional Brownian motion with origin in $\delta_{0}$ the Dirac measure at [math], $f$ a real-valued function and define

[TABLE]

where $\bar{\mu}$ denotes the centred version of $\mu$. We then have, using $X(\delta_{0})=0$ almost surely and (10):

[TABLE]

Hence the mean of the output $Y(\mu)$ is a function of the mean of the input distribution $\mu$ and its dispersion is an increasing function of the dispersion of $\mu$. This is a valuable property when modeling a function $\mu\mapsto g(\mu)$ as a Gaussian Process realization $\mu\mapsto Y(\mu)$, when it is believed that the range of possible values for $g$ increases with the variance of the input $\mu$.

Let us further notice that for $f=id$ and $H=1$ we have

[TABLE]

and

[TABLE]

where $F^{-1}_{\mu}$ denotes the quantile function of the distribution $\mu$, and $U$ is an uniform random variable on $[0,1]$. In some sense, $Y$ is in this case the Gaussian process that mimics the statistical properties of the optimal coupling $(F^{-1}_{\mu}(U))_{\mu\in\mathcal{W}_{2}(\mathbb{R})}$ (see Section III a ).

From now on (with the exception of Section -A from the appendix where we prove Theorem IV.1) we will focus on stationary processes, which are more adapted to learning tasks on distributions where there is no a priori reason to associate different dispersion properties to the outputs corresponding to different distribution inputs.

IV-B Stationary processes

We now construct Gaussian processes which are stationary with respect to the Wasserstein distance.

Theorem IV.2**.**

For every completely monotone function $F$ and $0<H\leq 1$ the function

[TABLE]

is a covariance function on $\mathcal{W}_{2}(\mathbb{R})$. Furthermore a Gaussian random process with constant mean and covariance (12) is stationary with respect to the Wasserstein distance.

We recall that a $\mathcal{C}^{\infty}$ function $F:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ is said to be completely monotone if for every $n\in\mathbb{N}$ and $x\in\mathbb{R}^{+}$,

[TABLE]

Here $F^{(n)}$ denotes the derivative of order $n$ of $F$. The prototype of a completely monotone fuction is $x\mapsto e^{-\lambda x}$, for any positive $\lambda$. Furthermore $F$ is completely monotone if and only if it is the Laplace transform of a positive measure $\mu_{F}$ with finite mass on $\mathbb{R}^{+}$, that is to say

[TABLE]

Other examples of completely monotone functions include $x^{-\lambda}$ for positive values of $\lambda$ and $\log\left(1+\frac{1}{x}\right)$.

• Example

Applying theorem IV.2 with the completely monotone functions $e^{-\lambda x}$ we obtain the stationary covariance kernels

[TABLE]

for every $\lambda>0$ and $0<H\leq 1$.

These kernels are generalizations to distribution inputs of the kernels of the form $e^{-\lambda\|x-y\|^{2H}}$ on $\mathbb{R}^{d}$, which are classical in spatial statistics and machine learning. In particular setting $H=1/2$ gives the family of Laplace kernels, and $H=1$ the family of Gaussian kernels.

At this point we have obtained enough covariance functions to consider parametric models that fit practical datasets. Section V addresses the question of the selection of the best covariance kernel amongst a parametric family of stationary kernels, together with the prediction of the associated Gaussian process. In Section VI we carry out simulations with the following parametric model, which is directly derived from (13):

[TABLE]

where $C,C^{\prime}\subset(0,\infty)$ are two compact sets.

IV-C Ideas of proof

Theorems IV.1 and IV.2 are direct corollaries of the following result:

Theorem IV.3**.**

The function $W_{2}^{2H}$ is a negative definite kernel if and only if $0\leq H\leq 1$. Furthermore, it is nondegenerate if and only if $0<H<1$.

One can find in [18] a proof of the negative definiteness of the kernel $W_{2}^{2H}$ restricted to absolutely continuous distributions in $\mathcal{W}_{2}(\mathbb{R})$. The proof given here holds for any distribution of $\mathcal{W}_{2}(\mathbb{R})$, and we provide the nondegeneracy property of the kernel.

In short (see Appendix -A for a detailed proof), we consider $H=1$ and the optimal coupling (see Section III a)

[TABLE]

where $F^{-1}_{\mu}$ is the quantile function of the distribution $\mu$ and $U$ is an uniform random variable on $[0,1]$. This coupling can be seen as a (non-Gaussian!) random field indexed by $\mathcal{W}_{2}(\mathbb{R})$. As such, its variogram

[TABLE]

is a negative definite kernel. Furthermore it is equal to $W_{2}^{2}(\mu,\nu)$ since the coupling $(Z(\mu))$ is optimal (see (1)). The proof ends with the use of the following classical lemma:

Lemma IV.4**.**

If $K$ is a negative definite kernel then $K^{H}$ is a negative definite kernel for every $0\leq H\leq 1$.

See $e.g.$ [32] for a proof Lemma IV.4.

• Remark

In [8], Istas defines the fractional index of a metric space $E$ endowed with a distance $d$ by

[TABLE]

One of the interpretation of the fractional index is that $\beta_{E}/2$ it is the maximal regularity for a fractional Brownian motion indexed by $(E,d)$: indeed the $H$-fractional Brownian motion indexed by a metric space exists if and only if $H\leq\beta_{E}/2$. For instance, the fractional exponent of the Euclidean spaces $\mathbb{R}^{n}$ is equal to $2$, while the fractional index of the spheres $\mathbb{S}^{n}$ is only $1$. Recall that an $H$-fractional Brownian motion has more regular paths and exhibits long-distance correlation for large values of $H$. In a non-rigorous way, the fractional index can be seen as some measure of the difficulty to construct long-distance correlated random fields indexed by the space $(E,d)$.

It is in general a difficult problem to find the fractional index of a given space. Theorem IV.3 states that the fractional exponent $\beta_{\mathcal{W}_{2}(\mathbb{R})}$ of the Wasserstein space is equal to $2$.

V Model selection and Gaussian process regression

V-A Maximum Likelihood and prediction

Let us consider a Gaussian process $Y$ indexed by $\mathcal{W}_{2}(\mathbb{R})$, with zero mean function and unknown covariance function $K_{0}$. Most classically, it is assumed that the covariance function $K_{0}$ belongs to a parametric set of the form

[TABLE]

with $\Theta\subset\mathbb{R}^{p}$ and where $K_{\theta}$ is a covariance function and $\theta$ is called the covariance parameter. Hence we have $K_{0}=K_{\theta_{0}}$ for some true parameter $\theta_{0}\in\Theta$.

For instance, considering the fractional Brownian motion kernel given in (8), we can have $\theta=(\sigma^{2},H)$, $\Theta=(0,\infty)\times(0,1]$ and $K_{\theta}=\sigma^{2}K^{H,\eta}$, where $\eta$ is fixed in $\mathcal{W}_{2}(\mathbb{R})$. In this case, the covariance parameters are the order of magnitude parameter $\sigma^{2}$ and the regularity parameter $H$.

Typically, the covariance parameter $\theta$ is selected from a data set of the form $(\mu_{i},y_{i})_{i=1,...,n}$, with $y_{i}=Y(\mu_{i})$. Several techniques have been proposed for constructing an estimator $\hat{\theta}=\hat{\theta}(\mu_{1},y_{1},...,\mu_{n},y_{n})$, in particular maximum likelihood (see e.g. [33]) and cross validation [34, 35, 36]. In this paper, we shall focus on maximum likelihood, which is widely used in practice and has received a lot of theoretical attention.

Maximum Likelihood is based on maximizing the Gaussian likelihood of the vector of observations $(y_{1},...,y_{n})$. The estimator is $\hat{\theta}_{ML}\in\operatorname{argmin}L_{\theta}$ with

[TABLE]

where $R_{\theta}=[K_{\theta}(\mu_{i},\mu_{j})]_{1\leq i,j\leq n}$

Given the maximum likelihood estimator $\hat{\theta}_{ML}$, the value $Y(\mu)$, for any input $\mu\in\mathcal{W}_{2}(\mathbb{R})$, can be predicted by plugging (see for instance in [33]) $\hat{\theta}_{ML}$ in the conditional expectation (or posterior mean) expression for Gaussian processes. More precisely, $Y(\mu)$ is predicted by $\hat{Y}_{\hat{\theta}_{ML}}(\mu)$ with

[TABLE]

and

[TABLE]

Note that $\hat{Y}_{\theta}(\mu)$ is the conditional expectation of $Y(\mu)$ given $y_{1},...,y_{n}$, when assuming that $Y$ is a centered Gaussian process with covariance function $K_{\theta}$.

V-B Asymptotic properties

In this section, we aim at showing that some of the asymptotic results of the Gaussian process literature, which hold for Gaussian processes indexed by $\mathbb{R}^{d}$, can be extended to Gaussian processes indexed by $\mathcal{W}_{2}(\mathbb{R})$. To our knowledge, this extension has not been considered before.

For a Gaussian process indexed by $\mathbb{R}^{d}$, two main asymptotic frameworks are under consideration: fixed-domain and increasing-domain asymptotics [33]. Under increasing-domain asymptotics, as $n\to\infty$, the observation points $x_{1},...,x_{n}\in\mathbb{R}^{d}$ are so that $\min_{i\neq j}||x_{i}-x_{j}||$ is lower bounded. Under fixed-domain asymptotics, the sequence (or triangular array) of observation points $(x_{1},...,x_{n})$ becomes dense in a fixed bounded subset of $\mathbb{R}^{d}$. To be specific, for a Gaussian process indexed by $\mathbb{R}$, a standard increasing-domain framework would be given by $x_{i}=i$ for $i\in\mathbb{N}$, while a standard fixed-domain framework would be given by, for $n\in\mathbb{N}$, $x_{i}=i/n$ for $i=1,...,n$.

Let us now briefly review the existing results for Gaussian processes indexed by $\mathbb{R}^{d}$. Typically, under increasing-domain asymptotics, the true covariance parameter $\theta_{0}$ is estimated consistently by maximum likelihood, with asymptotic normality [37, 38, 39, 40, 41, 42]. Also, predicting with the estimated covariance parameter $\hat{\theta}$ is asymptotically as good as predicting with $\theta_{0}$ [41].

Under fixed-domain asymptotics, there are cases where some components of the true covariance parameter $\theta_{0}$ can not be consistently estimated [33, 43]. Nevertheless, these components which can not be estimated consistently do not have an asymptotic impact on prediction [44, 45, 46]. Some results on prediction with estimated covariance parameters are available in [47]. Also, asymptotic properties of maximum likelihood estimators are obtained in [48, 49, 50, 51, 52].

We remark, finally, that the above increasing-domain asymptotic results hold for fairly general classes of covariance functions, while fixed-domain asymptotic results currently have to be derived for specific covariance functions and on a case-by-case basis.

For this reason, in this paper, we focus on extending some of the above increasing-domain asymptotic results to Gaussian processes indexed by $\mathcal{W}_{2}(\mathbb{R})$. Indeed, this will enable us to obtain a fair amount of generality with respect to the type of covariance functions considered.

We thus extend the contributions of [41] in the case of Gaussian processes with probability distribution inputs. In the rest of the section, we first list and discuss technical conditions for the asymptotic results. Then, we show the consistency and asymptotic normality of maximum likelihood and show that predictions from the maximum likelihood estimator are asymptotically as good as those obtained from the true covariance parameter. In Section V-C, we study an explicit example, for which all the technical conditions can be satisfied. All the proofs are postponed to the appendix. At the end of Section V-C, we discuss the novelty of these proofs, compared to those of the literature, and especially those in [41].

The technical conditions for this section are listed below.

Condition V.1**.**

We consider a triangular array of observation points $\{\mu_{1},...,\mu_{n}\}=\{\mu_{1}^{(n)},...,\mu_{n}^{(n)}\}$ so that for all $n\in\mathbb{N}$ and $1\leq i\leq n$, $\mu_{i}$ has support in $[i,i+L]$ with a fixed $L<\infty$.

Condition V.2**.**

The model of covariance functions $\{K_{\theta},\theta\in\Theta\}$ satisfies

[TABLE]

with $F_{\theta}:\mathbb{R}^{+}\to\mathbb{R}$ and

[TABLE]

with a fixed $A<\infty$, $\tau>1$.

Condition V.3**.**

We have observations $y_{i}=Y(\mu_{i})$, $i=1,\cdots,n$ of the centered Gaussian process $Y$ with covariance function $K_{\theta_{0}}$ for some $\theta_{0}\in\Theta$.

Condition V.4**.**

The sequence of matrices $R_{\theta}=\left(K_{\theta}(\mu_{i},\mu_{j})\right)_{1\leq i,j\leq n}$ satisfies

[TABLE]

for a fixed $c>0$, where $\lambda_{\inf}(R_{\theta})$ denotes the smallest eigenvalue of $R_{\theta}$.

Condition V.5**.**

$\forall\alpha>0$,

[TABLE]

Condition V.6**.**

$\forall t\geq 0$, $F_{\theta}(t)$ is continuously differentiable with respect to $\theta$ and we have

[TABLE]

with $A,\tau$ as in Condition V.2.

Condition V.7**.**

$\forall t\geq 0$, $F_{\theta}(t)$ is three times continuously differentiable with respect to $\theta$ and we have, for $q\in\{2,3\}$, $i_{1}\cdots i_{q}\in\{1,\cdots p\}$,

[TABLE]

Condition V.8**.**

$\forall(\lambda_{1}\cdots,\lambda_{p})\neq(0,\cdots,0)$,

[TABLE]

Condition V.1 mimics the increasing-domain asymptotic framework discussed above for vectorial inputs. In particular, the observation measures $\mu_{i}$ and $\mu_{j}$ yield a large Wasserstein distance when $|i-j|$ is large.

Condition V.2 entails that all the covariance functions under consideration are stationary in the sense that the covariance between $\mu$ and $\nu$ depends only on the distance $W_{2}(\mu,\nu)$. Stationarity is also assumed when considering increasing-domain asymptotics for Gaussian processes indexed by $\mathbb{R}^{d}$ [37, 38, 39, 40, 41, 42]. Hence, we remark that the asymptotic results of the present section do not apply to the covariance functions of fractional Brownian motion in (8). On the other hand, these results apply to the power exponential covariance functions in (14).

Condition V.2 also imposes that the covariance functions in the parametric model decrease fast enough with the Wasserstein distance. This condition is standard in the case of vector inputs, and holds for instance for the covariance functions in (14).

Condition V.3 means that we address the well-specified case [34, 35], where there is a true covariance parameter $\theta_{0}$ to estimate.

Condition V.4 is technically necessary for the proof techniques of this paper. This condition holds whenever the covariance model satisfies, for all $\theta\in\Theta,w\geq 0$, $F_{\theta}(w)=\bar{F}_{\theta}(w)+\delta_{\theta}\mathbf{1}_{\{w=0\}}$, where $\bar{F}_{\theta}$ is a continuous covariance function and where $\inf_{\theta\in\Theta}\delta_{\theta}>0$. This situation corresponds to Gaussian processes observed with Gaussian measure errors, or to Gaussian processes with very small scale irregularities, and is thus representative of a significant range of practical applications.

In the case where $F_{\theta}$ is continuous (which usually means that we have exact observations of a Gaussian process with continuous realizations), then Condition V.4 implies that

[TABLE]

For a large class of Gaussian processes indexed by $\mathbb{R}^{d}$, it has been shown that the condition in (21) (with $W_{2}$ replaced by the Euclidean distance) is also sufficient for Condition V.4 [41, 53]. The proof relies on the Fourier transform on $\mathbb{R}^{d}$. For Gaussian processes indexed by $\mathcal{W}_{2}(\mathbb{R})$, one could expect the condition in (21) to be sufficient to guarantee V.4 in many cases, although, to our knowledge, obtaining rigorous proofs in this direction is an open problem.

Condition V.5 means that there is enough information in the triangular array $\{\mu_{1},...,\mu_{n}\}$ to differentiate between the covariance functions $K_{\theta_{0}}$ and $K_{\theta}$, when $\theta$ is bounded away from $\theta_{0}$. We believe that Condition V.5 can be checked for specific explicit instances of the triangular array $\{\mu_{1},...,\mu_{n}\}$, as it involves an explicit sum of covariance values.

Conditions V.6 and V.7 are standard regularity and asymptotic decorrelation conditions for the covariance model. They hold, in particular, for the power exponential covariance model of (14).

Finally, Condition V.8 is interpreted as an asymptotic local linear independence of the $p$ derivatives of the covariance function, around $\theta_{0}$. Since this condition involves an explicit sum of covariance function derivatives, we believe that it can be checked for specific instances of the triangular array $\{\mu_{1},...,\mu_{n}\}$.

We now provide the first result of this section, showing that the maximum likelihood estimator is asymptotically consistent.

Theorem V.9**.**

Let $\hat{\theta}_{ML}$ be as in (19). Under Conditions V.1 to V.5, we have as $n\rightarrow\infty$

[TABLE]

In the next theorem, we show that the maximum likelihood estimator is asymptotically Gaussian. In addition, the rate of convergence is $\sqrt{n}$, and the asymptotic covariance matrix $M_{ML}^{-1}$ of $\sqrt{n}(\hat{\theta}_{ML}-\theta_{0})$ (that may depend on $n$) is asymptotically bounded and invertible, see (22).

Theorem V.10**.**

Let $M_{ML}$ be the $p\times p$ matrix defined by

[TABLE]

with $R_{\theta}$ as in (19). Under Conditions V.1 to V.8 we have

[TABLE]

Furthermore,

[TABLE]

In the next theorem, we show that, when using the maximum likelihood estimator, the corresponding predictions of the values of $Y$ are asymptotically equal to the predictions using the true covariance parameter $\theta_{0}$. Note that, in the increasing-domain framework considered here, the mean square prediction error is typically lower-bounded, even when using the true covariance parameter. Indeed, this occurs in the case of Gaussian processes with vector inputs, see Proposition 5.2 in [41].

Theorem V.11**.**

Under Conditions V.1 to V.8 we have

[TABLE]

with $\hat{Y}_{\theta}(\mu)$ as in (20).

V-C An example

In this section, we provide an explicit example of triangular array of probability measures for which Conditions V.5 and V.8 are satisfied. We consider random probability measures $(\mu_{i})_{i\in\mathbb{N}}$ which are independent and identically distributed (up to support shifts to satisfy condition V.1). We then show that Conditions V.5 and V.8 are satisfied almost surely. The motivation for studying shifted independent and identically distributed random probability measures is that this this model is simple to describe and can generate a large range of sequences $\{\mu_{1},...,\mu_{n}\}$.

Proposition V.12**.**

Assume that Conditions V.2, V.6 and V.7 hold.

Assume that for $\theta\neq\theta_{0}$, $F_{\theta}$ and $F_{\theta_{0}}$ are not equal everywhere on $\mathbb{R}^{+}$. Assume that there does not exist $(\lambda_{1},...,\lambda_{p})\neq(0,...,0)$ so that $\sum_{i=1}^{p}(\partial/\partial\theta_{i})F_{\theta_{0}}$ is the zero function on $\mathbb{R}^{+}.$

Let $(Z_{i})_{i\in\mathbb{Z}}$ be independent and identically distributed Gaussian processes on $\mathbb{R}$ with continuous trajectories. Assume that $Z_{0}$ has mean function [math] and covariance function $C_{0}$. Assume that $C_{0}(u,v)=C_{0}(u^{\prime},v^{\prime})$ whenever $v-u=v^{\prime}-u^{\prime}$ and let $C_{0}(u,v)=C_{0}(u-v)$ for ease of notation. Let $\hat{C}_{0}(w)=\int_{\mathbb{R}}C_{0}(t)e^{-\mathrm{i}wt}dt$ with $\mathrm{i}^{2}=-1$. Assume that $\hat{C}_{0}(w)|w|^{2q}$ is bounded away from [math] and $\infty$ as $|w|\to\infty$, for some fixed $q\in(0,\infty)$.

Let $L>1$ be fixed. For $i\in\mathbb{Z}$, let $f_{i}:\mathbb{R}\to\mathbb{R}^{+}$ be defined by $f_{i}(t)=\exp(Z_{i}(t-i))/M_{i}$ if $t\in[i,i+L]$ and $f_{i}(t)=0$ else, where $M_{i}=\int_{i}^{i+L}\exp(Z_{i}(t-i))dt$. Let $\mu_{i}$ be the measure with probability density function $f_{i}$. Then, almost surely, with the sequence of random probability measures $\{\mu_{1},...,\mu_{n}\}$, Conditions V.5 and V.8 hold.

In Proposition V.12, the identifiability assumptions on $\{F_{\theta}\}$ are very mild, and hold for instance for the power exponential model in (14).

In Proposition V.12, the random probability measures have probability density functions obtained from exponentials of realizations of Gaussian processes. Hence, these measures have a non-parametric source of randomness, and can take flexible forms. Several standard covariance functions on $\mathbb{R}$ satisfy the conditions in Proposition V.12, in particular the Matérn covariance functions (see e.g. [33]).

We remark that, in the context of Proposition V.12, when $F_{\theta}(w)=\bar{F}_{\theta}(w)+\delta_{\theta}\mathbf{1}_{\{w=0\}}$, with $\bar{F}_{\theta}$ a continuous covariance function and $\inf_{\theta\in\Theta}\delta_{\theta}>0$, as described when discussing Condition V.4, then Conditions V.1 to V.8 hold so that Theorems V.9, V.10 and V.11 hold. If however $F_{\theta}$ is continuous, then Condition V.4 almost surely does not hold since $L>1$ (as there will almost surely be pairs of distributions $\mu_{i},\mu_{i}$, $i\neq j$, with arbitrarily small $W_{2}(\mu_{i},\mu_{j})$). Nevertheless, when $L<1$, it can be shown that Proposition V.12 still holds when, in the conditions of this proposition on $\{F_{\theta}\}$, $\mathbb{R}^{+}$ is replaced by $\cup_{i=1}^{\infty}[i-L,i+L]$. Also, as discussed above, when $L<1$, the condition in (21) is satisfied and one could expect Condition V.4 to hold.

We conclude this section by discussing the corresponding proofs (in the appendix). These proofs can be divided into two groups. In the first group (proofs of Theorems V.9, V.10 and V.11 and of Proposition .7) we show that the arguments in [41] can be adapted and extended to the setting of the present article. The main innovations in this first group compared to [41] are that we allow for triangular arrays of observation points, and are not restricted to the specific structure of observation points of [41].

The proofs of the second group (proofs of Lemma .4 and Proposition V.12) are specific to Gaussian processes with distribution inputs and are thus original for the most part. In particular, in the proof of Proposition V.12, we show that, for two measures obtained by taking exponentials of Gaussian processes, the corresponding random Wasserstein distance has maximal distribution support. In this aim, we use equivalence of Gaussian measure tools and specific technical manipulations of the Wasserstein distance.

VI Simulation study

We now compare the Gaussian process model suggested in the present paper, with various models for predicting scalar outputs corresponding to distributional inputs. Among the covariance functions introduced in this paper, we shall focus on the power-exponential model (14), since its covariance functions are stationary with respect to the Wasserstein distance. We will not consider the fractional Brownian motion model (8), since it imposes to choose a “zero distribution”, from which the variance increases with the distance. While this feature is relevant in some applications (for instance in finance), it is not natural in the simulation examples adressed here.

VI-A Comparison with projection-based covariance functions

In this section, we focus on Gaussian process models for prediction. We compare the covariance functions (14) of this paper, operating directly on the input probability distributions, to more classical covariance functions operating on projections of these probability measures on finite dimensional spaces.

VI-A1 Overview of the simulation procedure

We address the input-output map given by, for a distribution $\nu$ on $\mathbb{R}$,

[TABLE]

where $m_{k}(\nu)=\int_{\mathbb{R}}x^{k}d\nu(x)$.

We first simulate independently $n=100$ learning distributions $\nu_{1},...,\nu_{100}$ as follows. First, we sample uniformly $\mu_{i}\in[0.3,0.7]$ and $\sigma_{i}\in[0.001,0.2]$, and compute $f_{i}$, the density of the Gaussian distribution with mean $\mu_{i}$ and variance $\sigma_{i}^{2}$. Then, we generate the function $g_{i}$ with value $f_{i}(x)\exp(Z_{i}(x))$, $x\in[0,1]$, where $Z_{i}$ is a realization of a Gaussian process on $[0,1]$ with mean function [math] and Matérn $5/2$ covariance function with parameters $\sigma=1$ and $\ell=0.2$ (see e.g. [54] for the expression of this covariance function). Finally, $\nu_{i}$ is the distribution on $[0,1]$ having density $g_{i}/(\int_{0}^{1}g_{i})$. In Figure 1, we show the density functions of $10$ of these $n$ sampled distributions. From the figure, we see that the learning distributions keep a relatively strong underlying two dimensional structure, driven by the randomly generated means and standard deviations. At the same time, because of the random perturbations generated with the Gaussian processes $Z_{i}$, these distributions are not restricted in a finite-dimensional space, and can exhibit various degrees of asymmetries.

From the learning set $(\nu_{i},F(\nu_{i}))_{i=1,...,n}$, we fit three Gaussian process models, which we call “distribution”, “Legendre” and “PCA”, and for which we provide more details below. Each of these three Gaussian process models provide a conditional expectation function

[TABLE]

and a conditional variance function

[TABLE]

We then evaluate the quality of the three Gaussian process models on a test set of size $n_{t}=500$ of the form $(\nu_{t,i},F(\nu_{t,i}))_{i=1,...,n_{t}}$, where the $\nu_{t,i}$ are generated in the same way as the $\nu_{i}$ above. We consider the two following quality criteria. The first one is the root mean square error (RMSE),

[TABLE]

which should be minimal. The second one is the confidence interval ratio (CIR) at level $\alpha\in(0,1)$,

[TABLE]

with $q_{\alpha}$ the $\left(\frac{1}{2}+\frac{\alpha}{2}\right)$ quantile of the standard normal distribution. The $CIR_{\alpha}$ criterion should be close to $\alpha$.

VI-A2 Details on the Gaussian process models

The “distribution” Gaussian process model is based on the covariance functions discussed before, operating directly on probability distributions. In this model, the Gaussian process has mean function zero and a covariance function of the form

[TABLE]

We call the covariance parameters $\sigma^{2}>0$, $\ell>0$ and $H\in[0,1]$ the variance, correlation length and exponent. These parameters are estimated by maximum likelihood from the training set $(\nu_{i},F(\nu_{i}))_{i=1,...,n}$, which yields the estimates $\hat{\sigma}^{2},\hat{\ell},\hat{H}$. Finally, the Gaussian process model for which the conditional moments $\hat{F}(\nu)$ and $\hat{\sigma}^{2}(\nu)$ are computed is a Gaussian process with mean function zero and covariance function $K_{\hat{\sigma}^{2},\hat{\ell},\hat{H}}$.

The “Legendre” and “PCA” Gaussian process models are based on covariance functions operating on finite-dimensional linear projections of the distributions. These projection-based covariance functions are used in the literature, in the general framework of stochastic processes with functional inputs, see e.g. [55, 56]. For the “Legendre” covariance function, for a distribution $\nu$ with density $f_{\nu}$ and support $[0,1]$, we compute, for $i=0,...,o-1$

[TABLE]

where $p_{i}$ is the $i-th$ normalized Legendre polynomial, with $\int_{0}^{1}p_{i}^{2}(t)dt=1$. The integer $o$ is called the order of the decomposition. Then, the covariance function operates on the input vector $(a_{0}(\nu),...,a_{o-1}(\nu))$ and is of the form

[TABLE]

The covariance parameters $\sigma^{2}\geq 0,\ell_{0}>0,...,\ell_{o-1}>0,H\in(0,1]$ are estimated by maximum likelihood, from the learning set $(a_{0}(\nu_{i}),...,a_{o-1}(\nu_{i}),F(\nu_{i}))_{i=1,...,n}$. Finally, the conditional moments $\hat{F}(\nu)$ and $\hat{\sigma}^{2}(\nu)$ are computed as for the “distribution” Gaussian process model.

For the “PCA” covariance function, we discretize each of the $n$ probability density functions $f_{\nu_{i}}$ to obtain $n$ vectors $v_{i}=(f_{\nu_{i}}(j/(d-1)))_{j=0,...,d-1}$, with $d=100$. Then, we let $w_{1},...,w_{o}$ be the first $o$ principal component vectors of the set of vectors $(v_{1},...,n_{n})$. For any distribution $\nu$ with density $f_{\nu}$, we associate its projection vector $(a_{1}(\nu),...,a_{o}(\nu))$ defined as

[TABLE]

This procedure corresponds to the numerical implementation of functional principal component analysis presented in Section 2.3 of [57]. Then, the covariance function in the “PCA” case operates on the input vector $(a_{1}(\nu),...,a_{o}(\nu))$. Finally, the conditional moments $\hat{F}(\nu)$ and $\hat{\sigma}^{2}(\nu)$ are computed as for the “Legendre” Gaussian process model.

VI-A3 Results

In Table I we show the values of the RMSE and $CIR_{0.9}$ quality criteria for the “distribution”, “Legendre” and “PCA” Gaussian process models. From the values of the RMSE criterion, the “distribution” Gaussian process model clearly outperforms the two other models. The RMSE of the “Legendre” and “PCA” models slightly decreases when the order increases, and stay well above the RMSE of the “distribution” model. Note that with orders $10$ and $15$, despite being less accurate, the “Legendre” and “PCA” models are significantly more complex to fit and interpret than the “distribution” model. Indeed these two models necessitate to estimate $12$ and $17$ covariance parameters, against $3$ for the “distribution” model. The maximum likelihood estimation procedure thus takes more time for the “Legendre” and “PCA” models than for the “distribution” model. We also remark that all three models provide appropriate predictive confidence intervals, as the value of the $CIR_{0.9}$ criterion is close to $0.9$. Finally, “Legendre” performs slightly better than “PCA”.

Our interpretation for these results is that, because of the nature of the simulated data $(\nu_{i},F(\nu_{i}))$, working directly on distributions, and with the Wasserstein distance, is more appropriate than using linear projections. Indeed, in particular, two distributions with similar means and small variances are close to each other with respect to both the Wasserstein distance and the value of the output function $F$. However, if the ratio between the two variances is large, the probability density functions of the two distributions are very different from each other, with respect to the $L^{2}$ distance. Hence, linear projections based on probability density functions is inappropriate in the setting considered here.

VI-B Comparison with the kernel regression procedure of [24]

In this section, we compare the “distribution” method of Table I which is suggested in the present article, with the “kernel regression” procedure of [24]. This procedure consists in predicting $f(P)\in\mathbb{R}$, with $P\in\mathcal{W}_{2}(\mathbb{R})$, from $\hat{P},\hat{P}_{1},...,\hat{P}_{n},f(P_{1}),...,f(P_{n})$ where $\hat{P},\hat{P}_{1},...,\hat{P}_{n}$ are estimates of $P,P_{1},...,P_{n}\in\mathcal{W}_{2}(\mathbb{R})$ obtained from sample values of $P,P_{1},...,P_{n}$. In [24], $\hat{P},\hat{P}_{1},...,\hat{P}_{n}$ correspond to kernel smoothing estimates of probability density functions constructed from the sample values. Then, the prediction $\hat{f}(\hat{P})$ of $f(P)$ is obtained by a weighted average of $f(P_{1}),...,f(P_{n})$ where the weights are computed by applying a kernel to the distances $D(\hat{P},\hat{P}_{1}),...,D(\hat{P},\hat{P}_{n})$. The distances suggested in [24] are the $L^{1}$ distances between the estimated probability density functions. We remark that there is no estimate of the prediction error $f(P)-\hat{f}(\hat{P})$ in [24], which is a downside compared to the Gaussian process model considered in this paper.

An interesting feature of the setting of [24] is that the input $P$ of the function value $f(P)$ is not observed. Only a sample from $P$ is available (this is the “two-stage sampling” difficulty described in [16], which arises in various applications) We shall demonstrate in this section that Gaussian process models can accommodate with this constraint. The idea is that $f(\hat{P})$ differs from $f(P)$, and that this difference can be modeled by adding a nugget variance parameter to the Gaussian process model. More precisely, the covariance functions we shall study in this section are

[TABLE]

where $\delta\geq 0$ is an additional covariance parameter, which can also be estimated in the maximum likelihood procedure. Apart from this modification of the covariance model, we carry out the Gaussian process model computation as in Section VI-A, with always $W_{2}(P,Q)$ replaced by $W_{2}(\check{P},\check{Q})$, where $\check{P},\check{Q}$ are the empirical distributions corresponding to the available sample values from $P,Q$.

We first reproduce the “skewness of Beta” example of [24]. In this example $n=275$ distributions $P_{1},...,P_{n}$ are randomly and independently generated for the learning set. We have that $P_{i}=B_{a_{i}}$ is the Beta distribution with parameters $(a_{i},b)$ where $a_{i}$ is uniformly distributed on $[3,20]$ and $b=3$. The test set consists in $n_{t}=50$ distributions $P_{t,1},...,P_{t,n_{t}}$ generated independently in the same way. The function to predict is defined by $f(P_{a})=[2(b-a)(a+b+1)^{1/2}]/[(a+b+2)(ab)^{1/2}]$ and corresponds to the skewness of the Beta distribution. For each distribution, $500$ sample values are available. For the “kernel regression” procedure, we used the same settings (kernel, bandwidth selection, training and validation sets…) as in [24].

The predictions obtained by the “distribution” and “kernel regression” procedures are presented in Figure 2. We observe that both methods perform equally well. The prediction errors are small, and are essentially due to to the fact that we only observe random samples from the distributions. [We have repeated the simulation of Figure 2 with $5,000$ sample values instead of $500$, and the predicted values have become visually equal to the true values.] Our conclusion on this “skewness of Beta” example is that the setting is here very favourable (the input space of distributions is one-dimensional and $275$ observations of the function are available) so that both methods have similar good performances.

Next, we repeat the “distribution” and “kernel regression” procedures on the same setting as in Table I (except that each input and predictand distribution is only observed indirectly, through $500$ sample values from it). The prediction results, based on the same criteria as in Table I are given in Table II. We observe that the RMSE prediction criterion for the “distribution” model is deteriorated compared to Table I. This is due to the fact that the distributions are not observed exactly anymore. The $CIR_{0.9}$ criterion is equal to $0.91$ for the “distribution” Gaussian process model. Hence, thanks to the addition of the nugget variance parameter, the Gaussian process model is able to take into account the additional uncertainty due to the random samples of the unobserved distributions, and to yield appropriate conditional variances.

We also observe that the RMSE pediction criterion is much larger for the “kernel regression” procedure. Hence, in this more challenging scenario (the input-space of distributions is non-parametric and only $100$ learning function values are available), the “distribution” Gaussian process model become strongly preferable. In our opinion, this is because the Wasserstein distance is here more relevant than distances between probability density functions (as discussed for Table I). Also, Gaussian process prediction has benefits compared to prediction with weighted kernel averages. In particular, Gaussian process predictions come with a probabilistic model and have optimality properties under this model.

VI-C Numerical complexity of the method

Our method inherits the numerical complexity of Gaussian process regression in more classical settings. Given a learning dataset $(\mu_{i},y_{i})_{i=1}^{N}$ the complexity of the Kriging method is given by the inversion of the covariance matrix $K_{\theta}(\mu_{i},\mu_{j})_{i,j=1}^{N}$, which is in $O(N^{3})$ number of operations. The Wasserstein distances between every pair of $\mu_{i}$ need also to be evaluated, which costs $O(N^{2}q)$ operations, where $q$ is the size of the sampling of the distributions.

Each prediction is then obtained by a vector product in $O(N)$ operations, while the computation of the conditional variance at some outputs is obtained in $O(N^{2})$.

The $O(N^{3})$ cost of the overall method makes it challenging to use on very large datasets, however on moderately large datasets its good performances makes it an interesting choice, and in particular a preferable choice over the other methods it was compared to in this simulation study.

For the sake of illustration, we remark that it took around 9 seconds to carry out our whole suggested Gaussian process procedure, in the case of Table I, and around 30 seconds in the case of Table II.

See also [58] for a discussion of the covariance tapering method to reduce the numerical cost of Gaussian process regression.

VII Conclusion

We provided a new approach to learning with distribution inputs. Its strength relies on the existence of positive definite kernels on the distribution space, which enables the use of Gaussian process models and kernel learning methods. In particular, we generalized the seminal models that are the fractional Brownian motion and the power exponential stationary processes, to distribution inputs. The kernels we use are functions of the Wasserstein distance, which has proven its efficiency as a discrepancy measure between distributions in numerous applications. Our method requires only the distributions inputs to have a second order moment, which allows the simultaneous handling of very heterogeneous data, such as absolutely continuous distributions, deterministic inputs and empirical distributions, which is particularly important when only a sample of the input distributions is known.

Focusing on Gaussian process regression with stationary covariance functions, we proved that our method extends this classical tool to distribution inputs. In particular, we gave generalization of state of the art asymptotic results to our setting. As in vector input Kriging, the overall numerical complexity of the method is in $O(n^{3})$, where $n$ is the size of the dataset, which is more costly than other distribution regression methods (such as the kernel regression procedure from [24]), however our numerical simulations suggest that our method gives better prediction. Furthermore Kriging comes with an error estimation in the form of the conditional variance of the Gaussian process, which is an important guarantee in practice.

On the down side, the methods we use to prove the positive definiteness of our kernels are tightly related to the existence of an optimal coupling between every distribution, which existence is specific to dimension one. It is an important problem for numerous applications to give learning methods for multidimensional distributions. Hence, it would be valuable to obtain kernels based on the multidimensional Wasserstein space. This would require an other approach that the one used in the present paper, and constitutes an interesting problem for further research.

[Proofs]

-A Proofs for Section IV

Proof of Theorem IV.3.

We start with the negative definiteness. For any $\mu\in\mathcal{W}_{2}(\mathbb{R})$ we denote by $F_{\mu}^{-1}$ the quantile function associated to $\mu$. It is well known that given a uniform random variable $U$ on $[0,1]$, $F_{\mu}^{-1}(U)$ is a random variable with law $\mu$, and furthermore for every $\mu,\nu\in\mathcal{W}_{2}(\mathbb{R})$:

[TABLE]

that is to say the coupling of $\mu$ and $\nu$ given by the random vector $(F_{\mu}^{-1}(U),F_{\nu}^{-1}(U))$ is optimal. Consider now $\mu_{1},\cdots,\mu_{n}\in\mathcal{W}_{2}(\mathbb{R})$ and $c_{1},\cdots,c_{n}\in\mathbb{R}$ such that $\sum_{i=1}^{n}c_{i}=0$. We have

[TABLE]

Using $\sum_{i=1}^{n}c_{i}=0$ the first two sums vanish and we obtain

[TABLE]

which proves that $W^{2H}_{2}$ is a negative definite kernel for $0\leq H\leq 1$.

Let us now consider $H>1$. Using (1) it is clear that for every $x,y\in\mathbb{R}$, $W_{2}(\delta_{x},\delta_{y})=|x-y|$. It is well known (see e.g [8]) that $|x-y]^{2H}$ is not a negative definite kernel on $\mathbb{R}$ for $H>1$, hence the same is true for $W_{2}^{2H}$.

Let us now prove the nondegeneracy of the kernel: the idea of the proof is adapted from [59]: we consider $\mathcal{W}_{2}(\mathbb{R})\times\mathbb{R}$ endowed with the product distance

[TABLE]

We assume the degeneracy of the kernel $W^{2H}_{2}$ on $\mathcal{W}_{2}(\mathbb{R})$ and deduce that $d^{2H}$ is not negative definite on $\mathcal{W}_{2}(\mathbb{R})\times\mathbb{R}$, in contradiction with the following Lemma, from which we postpone the proof:

Lemma .1**.**

The function $d^{2H}$ is a negative definite kernel if and only if $0\leq H\leq 1$.

Let us fix $0<H<1$ and assume that $W_{2}^{2H}$ is degenerate. There exists $\mu_{1},\cdots,\mu_{n}\in\mathcal{W}_{2}(\mathbb{R})$ and $c_{1},\cdots,c_{n}\in\mathbb{R}$ such that $\sum_{i=1}^{n}c_{i}=0$ and

[TABLE]

In $\mathcal{W}_{2}(\mathbb{R})\times\mathbb{R}$ we now consider the points $P_{i}=(\mu_{i},0)$ for $1\leq i\leq n$ and $P_{n+1}=(\mu_{n},\varepsilon)$ with $\varepsilon>0$. We also set $c^{\prime}_{i}=c_{i}$ for every $1\leq i\leq n-1$ and $c^{\prime}_{n}=c^{\prime}_{n+1}=c_{n}/2$. Notice that we have

[TABLE]

Now

[TABLE]

We now use

[TABLE]

to obtain

[TABLE]

Finally using (25) and $H<1$ we obtain

[TABLE]

which is positive for $\varepsilon$ small enough. This shows that $d^{2H}$ is not negative definite, in contradiction with Lemma .1. In the end $W_{2}^{2H}$ is nondegenerate for every $0<H<1$.

We now use the same argument as in the end of the proof of Theorem IV.3. Since $W_{2}^{2H}(\delta_{x},\delta_{y})=|x-y|^{2H}$ and $|x-y|^{2}$ and $|x-y|^{0}$ are degenerate kernels on $\mathbb{R}$, $W_{2}^{0}$ and $W_{2}^{2}$ are degenerate kernels.

∎

Proof of Lemma .1.

For $H=1$ we have

[TABLE]

hence $d^{2}$ is negative definite as the sum of two negative definite kernels. From Lemma IV.4 we get that $d^{2H}$ is a negative definite kernel for every $0\leq H\leq 1$.

For $H>1$ we notice that $d^{2H}(\mu,x)(\mu,y)=|x-y|^{2H}$ and use again the fact that $|x-y|^{2H}$ is not a negative definite kernel to conclude that $d^{2H}$ is not negative definite. ∎

Proof of Theorem IV.1.

The fact that (8) are covariance kernels is a direct consequence of Theorem IV.3 and the following Schoenberg Theorem (which is proven in [32]):

Theorem .2** (Schoenberg).**

Given a set $X$, two functions $K,R:X\times X\rightarrow\mathbb{R}$, and $o\in X$ such that for every $x,y\in X$,

[TABLE]

and

[TABLE]

the function $R$ is a positive definite kernel if and only if $K$ is a negative definite kernel.

We now prove the degeneracy: let $X=(X(\mu))_{\mu\in\mathcal{W}_{2}(\mathbb{R})}$ denote the $H$-fractional Brownian field indexed by $\mathcal{W}_{2}(\mathbb{R})$ with origin in $\sigma$. Assume $X$ is degenerate: there exist $\lambda_{1},\cdots,\lambda_{n}\in\mathbb{R}$ and $\mu_{1},\cdots,\mu_{n}\in\mathcal{W}_{2}(\mathbb{R})$ such that

[TABLE]

Since $X(\sigma)=0$ almost surely, setting $\mu_{n+1}=\sigma$ and $\lambda_{n+1}=-\sum_{i=1}^{n}\lambda_{i}$, it is clear that

[TABLE]

which implies

[TABLE]

Since $\sum_{i=1}^{n+1}\lambda_{i}=0$ this shows that $W_{2}^{2H}$ is degenerate, in contradiction with Theorem IV.3. Therefore $X$ is nondegenerate for every $0<H<1$.

The degeneracy of the [math]-fractional and the $2$-fractional Brownian field indexed by $\mathcal{W}_{2}(\mathbb{R})$ is a direct consequence from the degeneracy of $W_{2}^{0}$ and $W_{2}^{2}$.

∎

Proof of Theorem IV.2.

The fact that (12) are covariance kernels is a direct consequence of Theorem IV.3 and the following Schoenberg Theorem (which proof can be found in [32]):

Theorem .3** (Schoenberg).**

Let $F:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ be a completely monotone function, and $K$ a negative definite kernel. Then $(x,y)\mapsto F(K(x,y))$ is a positive definite kernel.

Furthermore as a function of the distance $W_{2}$, (12) is obviously invariant under the action of any isometry of $\mathcal{W}_{2}(\mathbb{R})$, so that the second claim holds. ∎

-B Proofs for Section V-B

Proof of Theorem V.9.

We have $\hat{\theta}_{ML}\in\operatorname{argmin}L_{\theta}$ with

[TABLE]

From Lemma .5 we have that

[TABLE]

are bounded as $n\rightarrow\infty$. Hence we can proceed as in the beginning of the proof of Proposition 3.1 in [41] to obtain

[TABLE]

Following again the proof of Proposition 3.1 in [41] we obtain the existence of a positive $a$ such that

[TABLE]

with $|\Lambda|^{2}=(1/n)\sum_{i,j=1}^{n}\Lambda_{i,j}^{2}$.

Hence from Condition V.5 and (26) we have $\forall\alpha>0$,

[TABLE]

and so

[TABLE]

∎

Proof of Theorem V.10.

From Lemma .5 and Condition V.4 we have for every $n\in\mathbb{N}$, $\left|\left(M_{ML}\right)_{i,j}\right|\leq B$ for a fixed $B<\infty$.

In addition, for any $\lambda_{1},\cdots,\lambda_{p}\in\mathbb{R}$ such that $\sum_{i=1}^{p}\lambda_{i}^{2}=1$,

[TABLE]

with a fixed $C>0$, since for every $n$

[TABLE]

from Lemma .5. Hence from Condition V.8 we obtain

[TABLE]

Hence (22) is proved. Let us now assume that

[TABLE]

Then there exists a bounded measurable function $g:\mathbb{R}^{p}\rightarrow\mathbb{R}$, $\xi>0$ and a subsequence $n^{\prime}$ such that along $n^{\prime}$ we have

[TABLE]

with $U\sim\mathcal{N}(0,I_{p})$.

In addition, by compactness, up to extracting another subsequence we can assume that

[TABLE]

where $M_{\infty}$ is a symmetric positive definite matrix.

Now the remaining of the proof is similar to the proof of Proposition 3.2 in [41]. We have

[TABLE]

Hence, exactly as in the proof of Proposition D.9 in [41] we can show

[TABLE]

Let us compute

[TABLE]

We have

[TABLE]

and from Condition V.4 and Lemma .6,

[TABLE]

Hence

[TABLE]

Moreover, $\displaystyle\frac{\partial^{3}}{\partial\theta_{i}\partial\theta_{j}\partial\theta_{k}}L_{\theta}$ can be written as

[TABLE]

where $A_{\theta}$ and $B_{\theta}$ are sums of products of the matrices $R_{\theta}^{-1}$ or $\displaystyle\frac{\partial}{\partial\theta_{i_{1}}}\cdots\frac{\partial}{\partial\theta_{i_{q}}}R_{\theta}$ with $q\in\{0,\cdots,3\}$ and $i_{1},\cdots,i_{q}\in\{1,\cdots p\}$.

Hence from Condition V.4 and from Lemmas .5 and .6 we have

[TABLE]

Following exactly the proof of Proposition D.10 in [41] we can show that

[TABLE]

Moreover since $M_{ML}\underset{n\rightarrow\infty}{\rightarrow}M_{\infty}$ we have

[TABLE]

This is in contradiction with (27) and concludes the proof. ∎

Proof of Theorem V.11.

From Theorem V.9 it is enough to show for $i=1,\cdots,p$ that

[TABLE]

From a version of Sobolev embedding theorem (see Theorem 4.12, part I, case A in [60]), there exists a finite constant $A_{\Theta}$ depending only on $\Theta$ such that

[TABLE]

Therefore in order to prove the theorem it is sufficient to show that for $w_{\theta}(\mu)$ of the form $r_{\theta}(\mu)$ or $\frac{\partial}{\partial\theta_{i}}r_{\theta}(\mu)$ or $\frac{\partial}{\partial\theta_{i}}\frac{\partial}{\partial\theta_{j}}r_{\theta}(\mu)$, and for $W_{\theta}$ equal to a product of the matrices $R_{\theta}^{-1}$ or $\frac{\partial}{\partial\theta_{i}}R_{\theta}$ or $\frac{\partial}{\partial\theta_{i}}\frac{\partial}{\partial\theta_{j}}R_{\theta}$, we have

[TABLE]

From Fubini theorem for positive integrands we have

[TABLE]

Now there exists a constant $c_{p+1}$ so that for $X$ a centred Gaussian random variable,

[TABLE]

hence

[TABLE]

Now from Lemmas .5 and .6 there exists $B<\infty$ such that

[TABLE]

Thus

[TABLE]

Finally for some $q\in\{0,1,2\}$ and for $i_{1},\cdots,i_{q}\in\{1,\cdots p\}$ we have

[TABLE]

with $C<\infty$ coming from Condition V.2, V.6, and V.7.

Using the proof of Lemma .4 we see that this quantity is bounded, which finishes the proof of Theorem V.11. ∎

-C Technical lemmas for Section V-B

Lemma .4**.**

[TABLE]

is bounded as $n\rightarrow\infty$.

Proof.

Let $\mu\in\mathcal{W}_{2}(\mathbb{R})$ and $i^{*}\in\operatorname{argmin}\limits_{k\in\{1,\cdots n\}}W_{2}(\mu_{k},\mu)$. For every $j\in\{1,\cdots,n\},$ $W_{2}(\mu,\mu_{j})\geq W_{2}(\mu,\mu_{i^{*}}).$ Moreover from the triangle inequality we have

[TABLE]

hence

[TABLE]

Let us define

[TABLE]

From Condition V.2 we have

[TABLE]

Now

[TABLE]

where for every $t\in[0,1]$

[TABLE]

Note that from Condition V.1 for every $t\in[0,1]$,

[TABLE]

If $|j-i^{*}|\geq L$ we have

[TABLE]

so that

[TABLE]

Hence

[TABLE]

∎

Lemma .5**.**

Under Conditions V.1 to V.4,

[TABLE]

and

[TABLE]

are bounded as $n\rightarrow\infty$.

Proof.

[TABLE]

is bounded as $n\rightarrow\infty$ from Lemma .4. The proof is similar for

[TABLE]

∎

In a similar way we also obtain the following Lemma.

Lemma .6**.**

$\forall q\in\{2,3\}$, $\forall i_{1},\cdots,i_{q}\in\{1,\cdots p\}$,

[TABLE]

is bounded as $n\rightarrow\infty$.

-D Proofs for Section V-C

Proposition .7**.**

Under the setting of Proposition V.12, almost surely as $n\to\infty$,

[TABLE]

and the sum in the right-hand side of the above display is a continuous function of $\theta$.

Proof.

Let

[TABLE]

Let $(m_{n})_{n\in\mathbb{N}}$ be a sequence of integers so that as $n\to\infty$, $m_{n}\to\infty$ and $n/m_{n}\to\infty$. Let

[TABLE]

With the same proof as that of Lemma D.11 in [41], we can show (using (28)) that $|S_{\theta}-S_{\theta,m_{n}}|$ goes almost surely to zero as $n\to\infty$. Also

[TABLE]

say. From (28), one can show simply that $r\to 0$ almost surely as $n\to\infty$. Also, the $B_{k}$ are independent random variables with identical distribution, and they are bounded in absolute value by

[TABLE]

from (28) and Condition V.2. Hence, applying Theorem 2.1 in [61] yields

[TABLE]

Hence, finally he have obtained almost surely as $n\to\infty$

[TABLE]

Also, we have, for $|i-j|\geq L$

[TABLE]

from (28). Hence, we can simply show

[TABLE]

with

[TABLE]

From (28) and Condition V.6, we can show that there exists a deterministic finite constant $C$ so that

[TABLE]

Also, by dominated convergence $T_{\theta}$ is a continuously differentiable function of $\theta$ and

[TABLE]

where $C^{\prime}$ is also a deterministic finite constant. Hence $\sup_{\theta\in\Theta}|S_{\theta}-T_{\theta}|\to 0$ almost surely as $n\to\infty$.

∎

Proof of Proposition V.12.

Assume that

[TABLE]

Then from Proposition .7 and by compacity, there exists $\theta_{1}\neq\theta_{0}$ so that

[TABLE]

From the conditions on $\{F_{\theta}\}$, there exist $\beta>0$, $\delta>0$, $a\geq 0$ so that for $u\in[a-\delta,a+\delta]$ we have $|F_{\theta_{1}}(u)-F_{\theta_{0}}(u)|\geq\beta$. Hence, we have

[TABLE]

for $k$ so that $a\in(k-1,k]$.

Let now $g_{0}:[0,L]\to\mathbb{R}^{+}$ be defined by $g_{0}(u)=D_{0}\exp(-1/(1-u^{2}))\mathbf{1}_{\{u\in[-1,1]\}}$ where $0<D_{0}<\infty$ is so that $\int_{\mathbb{R}}g_{0}(u)du=1$. Then, $g_{0}$ is infinitely differentiable. Let $h_{0}(u)=(1/\sigma)g_{0}((u-\delta/4)/\sigma)$ and $h_{k-1}(u)=(1/\sigma)g_{0}((u-a)/\sigma)$, where $\sigma>0$ is chosen small enough so that, with $\nu_{0}$ and $\nu_{k-1}$ the distributions with probability density functions $h_{0}$ and $h_{k-1}$ we have $W_{2}(\nu_{0},\nu_{k-1})\in[a-\delta/2,a+\delta/2]$ and $\nu_{0},\nu_{k-1}$ have supports in $[0,L],[k-1,k-1+L]$.

Let now $P_{1},P_{2}$ be two distributions with support in $[0,L]$, with quantile functions $q_{1},q_{2}$, with cumulative distribution functions $F_{1},F_{2}$ and with probability density functions $f_{1},f_{2}$. Then we have

[TABLE]

Let $\tau>0$ be so that $L^{3/2}\tau^{1/2}\leq\delta/5$. Then, for any $f:[0,L]\to\mathbb{R}$ and $g:[k-1,k-1+L]\to\mathbb{R}$, we have that $|f/(\int_{0}^{L}f)-h_{0}|_{\infty}\leq\tau$ and $|g/(\int_{0}^{L}g)-h_{k-1}|_{\infty}\leq\tau$ imply $W_{2}(\nu_{f},\nu_{g})\in[a-\delta,a+\delta]$, where $\nu_{f},\nu_{g}$ are the measures with probability density functions $f$ and $g$. Since $h_{0}$ and $h_{k-1}$ are infinitely differentiable, have integral one, and have respective supports included in $[0,L]$ and $[k-1,k-1+L]$, it is easy to see that there exists $\epsilon>0$ so that $|f-h_{0}|_{\infty}\leq\epsilon$ implies $|f/(\int_{0}^{L}f)-h_{0}|_{\infty}\leq\tau$, and similarly for $g$ and $h_{k-1}$. Hence, if we can show that

[TABLE]

and

[TABLE]

we obtain a contradiction. The two probabilities above are shown to be non-zero similarly and we will address the first one only. It is sufficient to show that

[TABLE]

Since $h_{0}+\epsilon/2$ is continuous and bounded away from [math] and infinity on $[0,L]$, it is sufficient to show that for all $\kappa>0$,

[TABLE]

From e.g. Theorem 1.1 in [62], since $z_{0}$ has mean function zero, we have

[TABLE]

Consider now the Gaussian measures $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$, on the space of continuous functions from $[0,L]\to\mathbb{R}$, so that $\mathcal{G}_{1}$ is the measure of the Gaussian process $Z_{0}$ and $\mathcal{G}_{2}$ is that of $Z_{0}-\log(h_{0}+\epsilon/2)$. Then, from e.g. the discussion in (22) in Chapter 4.2 of [33], since $\log(h_{0}+\epsilon/2)$ is infinitely differentiable, and from the assumptions on the covariance function of $Z_{0}$, the Gaussian measures $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ are equivalent. Hence, since

[TABLE]

we also have

[TABLE]

which is exactly

[TABLE]

This concludes the proof that Condition V.5 holds.

The proof that Condition V.8 holds can be obtained in the same way. In particular, an analog of Proposition .7 can be obtained. We skip the details. ∎

Acknowledgements

We thank the anonymous reviewers, whose suggestions have greatly contributed to improve the manuscript.

We thank Yann Richet, from the French Radioprotection and Nuclear Safety Institute (IRSN), for introducing us to the problem of axial burn up analysis of fuel pins, which motivated the present work.