Optimal compromise between incompatible conditional probability distributions, with application to Objective Bayesian Kriging

TL;DR

This paper introduces a method to find the optimal joint distribution that best fits incompatible conditional distributions, enabling practical Bayesian Kriging analysis with improved computational tractability and reliable prediction intervals.

Contribution

It proposes a novel approach to approximate incompatible conditionals with an optimal compromise distribution, facilitating objective Bayesian inference in Kriging models.

Findings

The pseudo-Gibbs sampler converges to the optimal compromise distribution.

The method achieves near-optimal frequentist coverage in prediction intervals.

It simplifies Bayesian correlation parameter estimation in Kriging models.

Abstract

Models are often defined through conditional rather than joint distributions, but it can be difficult to check whether the conditional distributions are compatible, i.e. whether there exists a joint probability distribution which generates them. When they are compatible, a Gibbs sampler can be used to sample from this joint distribution. When they are not, the Gibbs sampling algorithm may still be applied, resulting in a "pseudo-Gibbs sampler". We show its stationary probability distribution to be the optimal compromise between the conditional distributions, in the sense that it minimizes a mean squared misfit between them and its own conditional distributions. This allows us to perform Objective Bayesian analysis of correlation parameters in Kriging models by using univariate conditional Jeffreys-rule posterior distributions instead of the widely used multivariate Jeffreys-rule posterior. This strategy makes the full-Bayesian procedure tractable. Numerical examples show it has near-optimal frequentist performance in terms of prediction interval coverage.

Figures

Tables

Table 1. ?tablename? 1 . RMSE (where the error is measured in terms of the distance in Definition 5.1 ) of the MLE and MAP estimators for several “true” vectors of correlation lengths. The last column displays in percents the decrease of the RMSE of the MAP estimator with respect to the MLE.

Corr. lengths	MLE	MAP	– (%)
0.4 – 0.8 – 0.2	3.49	2.97	15
0.5 – 0.5 – 0.5	4.00	3.46	13
0.7 – 1.3 – 0.4	4.02	3.64	9
0.8 – 0.3 – 0.6	3.75	3.26	13
0.8 – 1.0 – 0.9	4.65	4.18	10

Table 2. ?tablename? 2 . Average with respect to randomly sampled design sets and realizations of the Gaussian Process (with variance parameter 1 and smoothness parameter 5/2) of the coverage of 95% Prediction Intervals across the sample space. “True” stands for the prediction based on the knowledge of the true variance parameter and the true vector of correlation lengths.

Corr. lengths	True	MLE	MAP	FPD
0.4 – 0.8 – 0.2	0.95	0.88	0.91	0.95
0.5 – 0.5 – 0.5	0.95	0.89	0.90	0.94
0.7 – 1.3 – 0.4	0.95	0.90	0.92	0.95
0.8 – 0.3 – 0.6	0.95	0.89	0.91	0.95
0.8 – 1.0 – 0.9	0.95	0.90	0.92	0.94

Table 3. ?tablename? 3 . Average with respect to randomly sampled design sets and realizations of the Gaussian Process (with variance parameter 1 and smoothness parameter 5/2) of the mean length of 95% Prediction Intervals across the sample space. The numbers in parentheses represent in percents the increase when using the MLE/MAP/FPD instead of the “true” vector of correlation lengths and variance parameter.

Corr. lengths	True	MLE	MAP	FPD
0.4 – 0.8 – 0.2	2.23	2.05 (-8)	2.13 (-4)	2.59 (+16)
0.5 – 0.5 – 0.5	1.69	1.55 (-8)	1.58 (-6)	1.84   (+9)
0.7 – 1.3 – 0.4	1.09	1.02 (-6)	1.07 (-2)	1.21 (+11)
0.8 – 0.3 – 0.6	1.63	1.51 (-7)	1.56 (-4)	1.82 (+12)
0.8 – 1.0 – 0.9	0.71	0.66 (-7)	0.69 (-3)	0.76   (+8)

Table 4. ?tablename? 4 . Coverage of 95 % prediction intervals when emulating the Ackley function on the unit hypercube using a Gaussian Process with null mean function and a Matérn anisotropic geometric covariance kernel with smoothness 5/2, unknown variance parameter and unknown vector of correlation lengths. The design sets contain 100 points.

Design set type	MLE	MAP	FPD
Unoptimized LHS	0.89	0.92	0.93
Optimized LHS	0.74	0.76	0.80
Random design	0.87	0.88	0.91

Table 5. ?tablename? 5 . Mean length of 95 % prediction intervals when emulating the Ackley function on the unit hypercube using a Gaussian Process with null mean function and a Matérn anisotropic geometric covariance kernel with smoothness 5/2, unknown variance parameter and unknown vector of correlation lengths. The design sets contain 100 points.

Design set type	MLE	MAP	FPD
Unoptimized LHS	0.31	0.33	0.35
Optimized LHS	0.24	0.24	0.28
Random design	0.28	0.29	0.32

Equations

$$p_1(x|y):=\frac{d{\pi }_1(\cdot |y)}{d\lambda }(x)$$

$$p_2(y|x):=\frac{d{\pi }_2(\cdot |x)}{d\lambda }(y)$$

$$\frac{d{m}_{\ne 1}}{d\lambda }(y)$$

$$\frac{d{m}_{\ne 2}}{d\lambda }(x)$$

$$\frac{d{\pi }_1{m}_{\ne 1}}{d\lambda }(x,y)=\frac{d{\pi }_2{m}_{\ne 2}}{d\lambda }(x,y)=\frac{\exp (-xy-2x-3y)}{{\int }_0^{+\infty }(t+6{)}^{-1}\exp (-t)dt}.$$

$$\frac{d{\pi }_1(\cdot |y)}{d\lambda }(x)/\frac{d{\pi }_2(\cdot |x)}{d\lambda }(y)=u(x)v(y).$$

$$\frac{p_1(x|y)}{p_2(y|x)}=\frac{(x+3{)}^{-1}\exp (-2x)}{(y+2{)}^{-1}\exp (-3y)}.$$

$${\varphi }_{-i}(S_{-i})=\varphi ({\text{bigtimes}}_{j<i}{S}_j\times {\Omega }_i\times {\text{bigtimes}}_{k>i}{S}_k).$$

$${\pi }_i{m}_{\ne i}(S_{<i}\times S_i\times S_{>i})={\int }_{S_{<i}\times S_{>i}}{\pi }_i(S_i|{\omega }_{-i})d{m}_{\ne i}({\omega }_{-i}).$$

$${\left({\pi }_i{m}_{\ne i}\right)}_{-i}=m_{\ne i}.$$

$$m_{\ne i}=\frac{1}{r}\sum _{j=1}^r({\pi }_j{m}_{\ne j}{)}_{-i}.$$

$${\pi }_1{m}_{\ne 1}=\mathcal{N}\left(\left(\begin{array}{c}0\\ 0\end{array}\right),\left(\begin{array}{cc}1& 1/2\\ 1/2& 2\end{array}\right)\right)\mathrm{and}{\pi }_2{m}_{\ne 2}=\mathcal{N}\left(\left(\begin{array}{c}0\\ 0\end{array}\right),\left(\begin{array}{cc}1& 1\\ 1& 2\end{array}\right)\right)$$

$$E_{\lambda }(P)=\sum _{i=1}^r{\int }_{\mathcal{A}}{\left[\frac{d({\pi }_i{P}_{-i})}{d\lambda }(\omega )-\frac{dP}{d\lambda }(\omega )\right]}^{2}d\lambda (\omega ).$$

$$P_G=\frac{1}{r}\sum _{i=1}^r{\pi }_i(P_G{)}_{-i}.$$

$${\left({\pi }_1{m}_{\ne 1}\right)}_{-2}=m_{\ne 2}\text{ and}{\left({\pi }_2{m}_{\ne 2}\right)}_{-1}=m_{\ne 1}.$$

$$\forall i,j,k\in [[1,r]],\forall S_k\in {\mathcal{A}}_k,{\pi }_i{m}_{\ne i}\left({\text{bigtimes}}_{k^{\prime }<k}{\Omega }_{k^{\prime }}\times S_k\times {\text{bigtimes}}_{k^{\prime \prime }>k}{\Omega }_{k^{\prime \prime }}\right)={\pi }_j{m}_{\ne j}\left({\text{bigtimes}}_{k^{\prime }<k}{\Omega }_{k^{\prime }}\times S_k\times {\text{bigtimes}}_{k^{\prime \prime }>k}{\Omega }_{k^{\prime \prime }}\right).$$

$$P(\{{\omega }_1\}\times \{{\omega }_2\}|{\omega }_3)=\frac{1}{2}{\pi }_1(\{{\omega }_1\}|{\omega }_2,{\omega }_3)P({\Omega }_1\times \{{\omega }_2\}|{\omega }_3)+\frac{1}{2}{\pi }_1({\omega }_2|{\omega }_1,{\omega }_3)P(\{{\omega }_1\}\times {\Omega }_2|{\omega }_3).$$

$$\varphi (\{{\omega }_1\}\times \{{\omega }_2\}\times \{{\omega }_3\})$$

$$=Q(\{{\omega }_2\}\times \{{\omega }_3\}|{\omega }_1)\varphi (\{{\omega }_1\}\times {\Omega }_2\times {\Omega }_3)$$

$$=R(\{{\omega }_1\}\times \{{\omega }_3\}|{\omega }_2)\varphi ({\Omega }_2\times \{{\omega }_2\}\times {\Omega }_3).$$

$$\frac{\varphi (\{{\omega }_1\}\times \{{\omega }_2\}\times \{{\omega }_3\})}{\varphi ({\Omega }_1\times \{{\omega }_2\}\times \{{\omega }_3\})}$$

$$=\frac{R(\{{\omega }_1\}\times \{{\omega }_3\}|{\omega }_2)}{R({\Omega }_1\times \{{\omega }_2\}|{\omega }_3)}=\frac{1}{2}{\pi }_1(\{{\omega }_1\}|{\omega }_2,{\omega }_3)+\frac{1}{2}{\pi }_3(\{{\omega }_3\}|{\omega }_1,{\omega }_2)\frac{R(\{{\omega }_1\}\times {\Omega }_3|{\omega }_2)}{R({\Omega }_1\times \{{\omega }_3\}|{\omega }_2)}.$$

$$\frac{1}{2}{\pi }_2(\{{\omega }_2\}|{\omega }_1,{\omega }_3)\frac{\varphi (\{{\omega }_1\}\times {\Omega }_2\times \{{\omega }_3\})}{\varphi ({\Omega }_1\times \{{\omega }_2\}\times \{{\omega }_3\})}=\frac{1}{2}{\pi }_3(\{{\omega }_3\}|{\omega }_1,{\omega }_2)\frac{\varphi (\{{\omega }_1\}\times \{{\omega }_2\}\times {\Omega }_3)}{\varphi ({\Omega }_1\times \{{\omega }_2\}\times \{{\omega }_3\})}.$$

$$\sum _{{\omega }_3\in {\Omega }_3}P(\{{\omega }_1\}\times \{{\omega }_2\}|{\omega }_3)\varphi ({\Omega }_1\times {\Omega }_2\times \{{\omega }_3\})$$

$$=$$

$${\varphi }_{-3}(\{{\omega }_1\}\times \{{\omega }_2\})=\frac{1}{2}\sum _{{\omega }_3\in {\Omega }_3}{\pi }_1(\{{\omega }_1\}|{\omega }_2,{\omega }_3){\varphi }_{-1}(\{{\omega }_2\}\times \{{\omega }_3\})+\frac{1}{2}\sum _{{\omega }_3\in {\Omega }_3}{\pi }_2(\{{\omega }_2\}|{\omega }_1,{\omega }_3){\varphi }_{-2}(\{{\omega }_1\}\times \{{\omega }_3\}).$$

$${\pi }_1(X_1=1|{\omega }_2)=1_{\{{\omega }_2=0\}}+1/21_{\{{\omega }_2=1\}};$$

$${\pi }_2(X_2=1|{\omega }_1)=1/21_{\{{\omega }_1=0\}}+1_{\{{\omega }_1=1\}}.$$

$${\pi }_C(\{{\omega }_1\}\times \{{\omega }_2\})$$

$${\pi }_E(\{{\omega }_1\}\times \{{\omega }_2\})$$

$${\pi }_1(X_1=1|{\omega }_2,{\omega }_3)$$

$${\pi }_2(X_2=1|{\omega }_1,{\omega }_3)$$

$${\pi }_3(X_3=1|{\omega }_1,{\omega }_2)$$

$$P(\{{\omega }_1\}\times \{{\omega }_2\}|{\omega }_3)=1/4.$$

$${\pi }_C(\{{\omega }_1\}\times \{{\omega }_2\}\times \{{\omega }_3\})=\frac{1}{24}{1}_{\{{\omega }_1={\omega }_2,{\omega }_3=0\}}+\frac{1}{12}{1}_{\{{\omega }_1\ne {\omega }_2,{\omega }_3=0\}}+\frac{5}{24}{1}_{\{{\omega }_1={\omega }_2,{\omega }_3=1\}}+\frac{1}{6}{1}_{\{{\omega }_1\ne {\omega }_2,{\omega }_3=1\}}.$$

Full text

{resume}

Les modèles statistiques sont souvent définis par lois de probabilité conditionnelles plutôt que jointes. Il peut cependant être difficile de contrôler la compatibilité des lois conditionnelles, i.e. l’existence d’une loi jointe qui les engendre. Quand les lois conditionnelles sont compatibles, un échantillonneur de Gibbs peut être utilisé pour échantillonner selon la loi jointe dont elles procèdent. Dans le cas contraire, l’algorithme de l’échantillonneur de Gibbs peut parfois tout de même être appliqué. Une telle procédure est appelée « pseudo-échantillonneur de Gibbs ». Nous montrons que sa distribution stationnaire est le compromis optimal entre les lois conditionnelles, au sens où elle minimise un écart quadratique moyen entre elles et ses propres lois conditionnelles. Cela permet une analyse bayésienne objective des paramètres de corrélation de modèles de krigeage : nous utilisons des lois a posteriori de Jeffreys univariées conditionnelles plutôt que la très usitée loi a posteriori de Jeffreys multivariée. Cette stratégie rend la procédure pleinement bayésienne abordable. Des exemples numériques montrent que ses performances fréquentistes en termes de taux de couverture des intervalles prédictifs sont quasi-optimales.

1. Introduction

Generally speaking, there are two ways to create statistical models for multiple random variables. One can either consider them simultaneously and directly define their joint distribution, or one can define a system of conditional distributions. The first approach is conceptually easier and often (but not always) leads to models with well understood properties because a closed-form expression is available. The second one allows for more flexibility in modeling but makes theoretical analysis more difficult.

The main problem with the second approach is that conditional distributions may not be compatible. In this context, it means that there exists no joint distribution from which the conditional distributions can all be derived. Other definitions of compatibility exist in the literature. For example, in the context of a model with a given prior distribution, Dawid and Lauritzen [11] examine the problem of eliciting a compatible prior distribution for a submodel. In the domain of Bayesian Networks, a probability distribution can be compatible or not with a given Directed Acyclic Graph (DAG) [30]. Moreover, an abstraction (simplification) of a DAG can be compatible or not [34]. A concept of compatibility of two prior distributions also exists in the field of Bayesian model selection. It is based on the Kullback-Leibler divergence of the corresponding marginal (predictive) distributions [9]. In this paper however, the notion of compatibility concerns families of conditional distributions. A family of conditional distributions is called compatible if there exists a joint distribution that agrees with them all [16]. Uniqueness of this joint distribution is desirable but not included in the requirements of compatibility.

To illustrate this definition of compatibility, consider the following random-effect model ([16], cited by Robert [29] page 41), which shows that “in general, reasonable-seeming conditional models will not be compatible with any single joint distribution” [12]. Define $y_{ij}=\beta+u_{i}+\epsilon_{ij}$, $i\in[\![1,I]\!]$ and $j\in[\![1,J]\!]$ where $u_{i}\sim\mathcal{N}(0,\sigma^{2})$ and $\epsilon_{ij}\sim\mathcal{N}(0,\tau^{2})$. The parameters of the model being $(\beta,\sigma^{2},\tau^{2})$, let us consider the prior distribution $\pi(\beta,\sigma^{2},\tau^{2})\propto(\sigma^{2}\tau^{2})^{-1}$: the corresponding posterior is improper. The conditional posterior distributions cannot therefore be compatible. Nevertheless, they are all proper. With $\bar{y}$ and $\bar{u}$ being the empirical means of $\boldsymbol{y}=(y_{ij})_{i\in[\![1,I]\!],j\in[\![1,J]\!]}$ and $\boldsymbol{u}=(u_{i})_{i\in[\![1,I]\!]}$, $\beta|\boldsymbol{u},\boldsymbol{y},\sigma^{2},\tau^{2}\sim\mathcal{N}(\bar{y}-\bar{u},\tau^{2}/IJ)$, $\sigma^{2}|\boldsymbol{u},\boldsymbol{y},\beta,\tau^{2}\sim\mathcal{IG}\left(I/2,\sum_{i}u_{i}^{2}/2\right)$, $\tau^{2}|\boldsymbol{u},\beta,\boldsymbol{y},\sigma^{2}\sim\mathcal{IG}\left(IJ/2,\sum_{ij}(y_{ij}-u_{i}-\beta)^{2}/2\right)$. If the posterior distribution were presented only through these conditionals, one might miss the fact that Gibbs sampling is impossible in this case due to the corresponding Markov chain being null recurrent.

For finite state spaces, null recurrent Markov chains are impossible, but attempting to define a joint probability distribution through its conditional distributions may still lead to incompatibility. Accordingly, the problem of efficiently determining whether a given system of conditionals is compatible has received considerable attention over the years. Kuo et al. [22], after listing previous attempts, provide probably the best solution to date. Their idea relies on the Structural Ratio Matrix which contains ratios between conditional distributions.

However, even if a system contains incompatible conditional probability distributions, it does not follow that it is useless. Since Heckerman et al. [15], practitioners have been using systems of conditional probability distributions without reference to compatibility. Indeed, providing the Markov chain is positive recurrent, it is always possible to fire up Gibbs samplers to deal with a system of conditional distributions. Some authors use the colorful acronym PIGS for “Potentially Incompatible Gibbs Sampler” to describe such a procedure. When the conditionals are definitely known to be incompatible, the most widely used term seems to be “Pseudo-Gibbs Sampler” (PGS).

Behind the practice of PIGS is the intuition that the Gibbs sampler should converge to the joint distribution that best represents the system of conditionals. Kuo and Wang [21] provide a detailed analysis and geometrical interpretation of the behavior of Pseudo-Gibbs Samplers for discrete conditional distributions. In particular, they show how the scanning order determines its stationary distribution. In Section 2 of the present paper, we provide some theoretical foundation for the intuition that the stationary distribution of a PGS with random scanning order is, in case of uniqueness, the best “compromise” between incompatible conditionals. Section 3 provides further discussion of this theory by considering alterations to the main definitions and showing them to lead to undesirable results.

In Section 4, we use the theory of optimal compromise to derive Objective Bayesian inference on correlation parameters of Kriging models. Kriging models are widely used in spatial statistics, but they are more complex than standard models. Their parameters are numerous, and their interactions complex, which makes eliciting a joint Objective prior difficult. This makes an approach resting on conditional priors attractive despite the risk of incompatibility of the associated conditional posteriors. It is to deal with such situations that the theory of optimal compromise was developed.

The Objective Bayesian paradigm as explained by Berger [4] consists in eliciting for every model a “default”, reasonable prior distribution that could be used when no explicit prior information is available. In particular, the Berger-Bernardo reference prior [8], hereafter simply named “reference prior”, can be algorithmically computed with minimal user intervention.

For models with a single scalar parameter, the reference prior rewards parameter values that are easily discriminated by the likelihood function. Its definition is related to the Kullback-Leibler divergence between posterior and prior [8]. For usual continuous models – essentially models where the Fisher information matrix is equal to the opposite of the expectancy of the second derivative of the log-likelihood, see Clarke and Barron [10] for an exhaustive list of conditions – it coincides with the Jeffreys-rule prior.

For models with multiple parameters, the reference prior algorithm requires the user to specify an ordering on the parameters and then iteratively compute the reference prior on each parameter conditionally to all subsequent parameters. The only user input is therefore this ordering, and common sense arguments often make one more sensible than others. Yang and Berger [33] list different reference priors obtained with different parameter orderings for a large number of statistical models. Of course, one could also group several parameters and treat them as one single multi-dimensional parameter, but doing so tends to produce less satisfactory inference [5]. In particular, for usual continuous models Berger et al. [7] state “ We actually know of no multivariable example in which we would recommend the Jeffreys-rule prior. In higher dimensions, the prior always seems to be either ‘too diffuse’ […] or ‘too concentrated’ ”.

Berger et al. [6] were the first to derive a reference prior for the parameters of a Gaussian Process regression model. This model contained only one correlation parameter, however. When several correlation parameters are involved, there is no reasonable way to order them. Even if one were arbitrarily picked, computation of the prior would be analytically intractable. Several authors [25, 27, 19, 28, 13] have therefore resolved to treat all correlation parameters as a single multidimensional parameter. It is in order to avoid having to do this that we make use of PIGS.

The idea is simple: for every correlation parameter, it is possible to analytically derive the reference prior for this parameter conditionally to all others. Each of the corresponding posterior distributions can be seen as a conditional probability distribution on one correlation parameter when all others are known. These conditional distributions then serve as input to a PIGS.

Theorem 4.2 is the main result with respect to the application.

First, under reasonable assumptions, the PIGS admits one single stationary probability distribution. Second, the Markov kernel defined by the PIGS is uniformly ergodic. Since this Markov kernel is defined over an uncountable state space, the latter fact is significant. The stationary distribution, which we call the Gibbs reference posterior distribution, can be used to improve prediction of the value taken by the Gaussian process at unobserved points. Sections 5 and 6 illustrate the inferential and predictive performance of the stationary distribution respectively.

2. Optimal compromise: a general theory

2.1. Definitions and notations

In this section we introduce the concepts necessary to define the optimal compromise between potentially incompatible conditional distributions. For the sake of readability, all proofs are provided in Appendix A.

First, note that in this context, “conditional distribution” is really an informal way of referring to a Markov kernel.

Definition 2.1**.**

Let $(A,\mathcal{A})$ and $(B,\mathcal{B})$ be measurable sets. A mapping $\pi:A\times\mathcal{B}\rightarrow[0,1]$ is called a Markov kernel if:

1. (1)

for all $x\in A$, $\pi(x,\cdot):\mathcal{B}\rightarrow[0,1]$ is a probability distribution and 2. (2)

for all $S\in\mathcal{B}$, $\pi(\cdot,S):A\rightarrow[0,1]$ is $\mathcal{A}$-measurable.

We use the following notation: for every $(x,S)\in A\times\mathcal{B}$, $\pi(S|x):=\pi(x,S)$.

Let $r$ be a positive integer and let $(\Omega_{1},\mathcal{A}_{1})$,…,$(\Omega_{r},\mathcal{A}_{r})$ be measurable sets. Define $\Omega=\bigtimes_{i=1}^{r}\Omega_{i}=\Omega_{1}\times...\times\Omega_{r}$ and $\mathcal{A}:=\bigotimes_{i=1}^{r}\mathcal{A}_{r}=\mathcal{A}_{1}\otimes...\otimes\mathcal{A}_{r}$.

For every $i\in[\!|1,r]\!]$, let $\pi_{i}$ be a Markov kernel $\left(\bigtimes_{j\neq i}\Omega_{j}\right)\times\mathcal{A}_{i}\rightarrow[0,1]$.

Intuitively (we formalize this below), every $\pi_{i}$ should be assembled with a distribution $m_{\neq i}$ on $\bigotimes_{j\neq i}\mathcal{A}_{j}$ to create a “joint” distribution, that is a probability distribution on $\mathcal{A}$. We refer to every $m_{\neq i}$ ($i\in[\![1,r]\!])$ as an $(r-1)$-dimensional distribution. If the $m_{\neq i}$ can be chosen in such a way as to make all joint distributions equal, then the Markov kernels in the sequence $(\pi_{i})_{i\in[\![1,r]\!]}$ are called compatible. And if no choice of $(m_{\neq i})_{i\in[\![1,r]\!]}$ can make all joint distributions equal, we have to look for a “compromise” between the Markov kernels.

Consider the following example [3] with $r=2$ and $\Omega_{1}=\Omega_{2}=(0,+\infty)$ are endowed with their Borel $\sigma$-algebra: let $\pi_{1}$ and $\pi_{2}$ be Markov kernels such that for all $x,y>0$ $\pi_{1}(\cdot|y)$ and $\pi_{2}(\cdot|x)$ are absolutely continuous with respect to the Lebesgue measure on $(0,+\infty)$. Let $\lambda$ be this measure and let the densities of $\pi_{1}$ and $\pi_{2}$ respectively be

[TABLE]

Let us define $m_{\neq 1}$ and $m_{\neq 2}$ as the probability distributions with the following densities with respect to the Lebesgue measure on $(0,+\infty)$:

[TABLE]

Then the joint probability distributions $\pi_{1}m_{\neq 1}$ and $\pi_{2}m_{\neq 2}$ are equal. Now denoting by $\lambda$ the Lebesgue measure on $(0,+\infty)\times(0,+\infty)$, the density of $\pi_{1}m_{\neq 1}=\pi_{2}m_{\neq 2}$ is

[TABLE]

Remark* (Producing incompatibility is easy).*

Take $r=2$ and let $\Omega_{1}=\Omega_{2}$ be a Borel subset of $\mathbb{R}$. Assume that for all $x,y\in\Omega_{1}$, $\pi_{1}(\cdot|y)$ and $\pi_{2}(\cdot|x)$ are absolutely continuous with respect to $\lambda$, which denotes here the Lebesgue measure on $\Omega_{1}$. Let $p_{1}(\cdot|y)$ and $p_{2}(\cdot|x)$ be their respective density functions and further assume that for $\lambda$-almost all real numbers $x$ and $y$, $p_{1}(x|y)>0$ and $p_{2}(y|x)>0$. A necessary condition [3] for the compatibility of $\pi_{1}$ and $\pi_{2}$ can be derived from Bayes’ rule: there must exist two mappings $u$ and $v$ defined on $\Omega_{1}$ such that for $\lambda$-almost all real numbers $x$ and $y$

[TABLE]

In the previous example, $\Omega_{1}=\Omega_{2}=(0,+\infty)$ and this necessary condition is fulfilled: for all $x,y>0$,

[TABLE]

Taking $p_{1}(x|y)=(y+2)\exp(-(y+2)x)$ and $p_{2}(y|x)=\exp(-y)$ makes $\pi_{1}$ and $\pi_{2}$ fail the necessary condition for compatibility.

Definition 2.2**.**

Let $\phi$ be a probability distribution on $\mathcal{A}$. For every $i\in[\![1,r]\!]$, denote by $\phi_{-i}$ the probability distribution on $\bigotimes_{j\neq i}\mathcal{A}_{j}$ defined as follows. For every set $S_{-i}$ that can be decomposed as $S_{-i}=\bigtimes_{j\neq i}S_{j}$ (with $S_{j}\in\mathcal{A}_{j}$ for every $j\neq i$),

[TABLE]

$\phi_{-i}$* is called the $i$-th $(r-1)$-marginal distribution of $\phi$.*

Remark*.*

The above definition is valid because any probability distribution on $\mathcal{A}$ can be characterized by its values on “rectangles” $\bigtimes_{i=1}^{r}S_{i}$ (where for every $i\in[\![1,r]\!]$, $S_{i}\in\mathcal{A}_{i}$).

For every $i\in[\![1,r]\!]$ and every probability distribution $m_{\neq i}$ on $\bigotimes_{j\neq i}\mathcal{A}_{j}$, denote by $\pi_{i}m_{\neq i}$ the distribution on $\mathcal{A}$ defined as follows. For every $i\in[\![1,r]\!]$, for every set $S_{<i}\in\bigotimes_{j<i}\mathcal{A}_{j}$, every set $S_{>i}\in\bigotimes_{k>i}\mathcal{A}_{k}$ and every set $S_{i}\in\mathcal{A}_{i}$,

[TABLE]

Naturally, for $i=1$ (resp. $i=r$), remove $S_{<i}$ (resp. $S_{>i}$) from the formula above. In the following, do this kind of operation when $i=1$ or $i=r$.

Notice that for every $i\in[\![1,r]\!]$, $m_{\neq i}$ is the $i$-th $(r-1)$-marginal distribution of $\pi_{i}m_{\neq i}$:

[TABLE]

If there exists a sequence of $(r-1)$-dimensional distributions $(m_{\neq i})_{i\in[\![1,r]\!]}$ such that all distributions $\pi_{i}m_{\neq i}$ are equal, then the Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ are compatible. If no such sequence $(m_{\neq i})_{i\in[\![1,r]\!]}$ exists, then we wish to find a sequence $(m_{\neq i})_{i\in[\![1,r]\!]}$ that makes the $\pi_{i}m_{\neq i}$ share some “common ground”. The following definition expresses this constraint formally.

Definition 2.3**.**

A sequence of $(r-1)$-dimensional distributions $(m_{\neq i})_{i\in[\![1,r]\!]}$ (each $m_{\neq i}$ being a probability distribution on $\bigotimes_{j\neq i}\mathcal{A}_{j}$) is said to be compatible with the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ if for every $i\in[\![1,r]\!]$

[TABLE]

So the “common ground” we require for the sequence of distributions $(\pi_{i}m_{\neq i})_{i\in[\![1,r]_{!}]}$ is that their $(r-1)$-marginal distributions should be the same on average. Other constraints would have been possible, and we discuss some of them in Section 3.1 below.

Consider the case where $r=2$, $\Omega_{1}=\Omega_{2}=\mathbb{R}$ and for all $x,y\in\mathbb{R}$, $\pi_{1}(\cdot|y)=\mathcal{N}(y/4,1/8)$ and $\pi_{2}(\cdot|x)=\mathcal{N}(x,1)$, the second argument of $\mathcal{N}(\cdot,\cdot)$ being the variance. Now define $m_{\neq 1}=\mathcal{N}(0,2)$ and $m_{\neq 2}=\mathcal{N}(0,1)$. We have

[TABLE]

Since $(\pi_{1}m_{\neq 1})_{-2}=\mathcal{N}(0,1)=m_{\neq 2}$ and $(\pi_{2}m_{\neq 2})_{-1}=\mathcal{N}(0,2)=m_{\neq 1}$, we have a fortiori $m_{\neq 1}=1/2(\pi_{1}m_{\neq 1})_{-1}+1/2(\pi_{2}m_{\neq 2})_{-1}$ and $m_{\neq 2}=1/2\;(\pi_{1}m_{\neq 1})_{-2}+1/2\;(\pi_{2}m_{\neq 2})_{-2}$, so $(m_{\neq 1},m_{\neq 2})$ is compatible with $(\pi_{1},\pi_{2})$ in the sense of Definition 2.3.

The definition of a compromise follows from this new definition of compatibility.

Definition 2.4**.**

A probability distribution $P$ on $\mathcal{A}$ is called a compromise between the Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ if these two conditions are verified:

1. (1)

for every $i\in[\![1,r]\!]$, $\pi_{i}P_{-i}$ is absolutely continuous with respect to $P$; 2. (2)

the sequence $(P_{-i})_{i\in[\![1,r]\!]}$ of $P$’s $(r-1)$-marginal distributions is compatible with $(\pi_{i})_{i\in[\![1,r]\!]}$.

In the definition of a compromise, the first condition exists to give meaning to the definition of an optimal compromise below. It is reasonable on its own though: a compromise should not deem events impossible if they are considered possible by the Markov kernels.

Returning to the previous example, both $\pi_{1}m_{\neq 1}$ and $\pi_{2}m_{\neq 2}$ are compromises, and any convex combination of these two distributions is one as well. The following definition introduces a cost function for compromises in order to determine which compromises are optimal.

Definition 2.5**.**

Let $\lambda$ be a positive measure on $\mathcal{A}$. Let $P$ be a compromise between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ that is absolutely continuous with respect to $\lambda$. $P$ is called an optimal compromise with respect to $\lambda$ between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ if it minimizes the functional $E_{\lambda}$ over all compromises between $(\pi_{i})_{i\in[\![1,r]\!]}$ that are absolutely continuous with respect to $\lambda$. $E_{\lambda}$ is defined by:

[TABLE]

In the previous example, the optimal compromise with respect to the Lebesgue measure on $\mathbb{R}\times\mathbb{R}$ can be shown (cf. Section 2.2) to be $1/2\;\pi_{1}m_{-1}+1/2\;\pi_{2}m_{-2}$.

Proposition 2.6**.**

The set of all compromises between $(\pi_{i})_{i\in[\![1,r]\!]}$ is convex, as is the subset of all compromises absolutely continuous with respect to $\lambda$.

If the Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ are compatible and there exists a joint distribution $\pi$ on $\mathcal{A}$ that agrees with them all, then for every positive measure $\lambda$ on $\mathcal{A}$ such that $\pi$ is absolutely continuous with respect to $\lambda$, $E_{\lambda}(\pi)=0$ and $\pi$ is an optimal compromise with respect to $\lambda$.

Even though Definition 2.5 makes it seem like the notion of optimal compromise is tied to a reference measure $\lambda$, it turns out that in many situations there exists a compromise that is optimal with respect to all possible reference measures – cf. Theorem 2.10 below. In the previous example, it is given by $1/2\pi_{1}m_{-1}+1/2\pi_{2}m_{-2}$.

2.2. Deriving the optimal compromise

The concepts of compromise and optimal compromise defined in Section 2.1 are intimately linked to Gibbs sampling, or (because the Markov kernels are incompatible) pseudo-Gibbs sampling (PGS). Let us therefore recall the Gibbs sampling algorithm.

The method used to sample the index $i$ is called scanning order. A systematic scan moves through each index in turn using a deterministic pattern while an equiprobable random scan selects for each step an index within $[\![1,r]\!]$ with probability $1/r$.

If the Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ are compatible, then the scanning order influences only the rate of convergence of the algorithm. He et al. [14] provide theoretical results on the subject and Mitliagkas and Mackey [24] propose a measure for the quality of any scan order for Gibbs sampling on finite state spaces.

If the Markov kernels are not compatible, then every scan order may produce a different target distribution. Kuo and Wang [21] thoroughly study the links between all possible systematic scan orders for Markov kernels on finite state spaces. The equiprobable random scan order has yet another target distribution.

The notion of Gibbs compromise is central to this theory of compromises between incompatible Markov kernels.

Definition 2.7**.**

A probability distribution $P_{G}$ on $\mathcal{A}$ is called a Gibbs compromise between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ if it satisfies:

[TABLE]

If it exists, a Gibbs compromise is a stationary distribution for the Gibbs sampler with equiprobable random scan order. This fact makes it practical from a sampling standpoint. In the following, we show how it relates to the concepts of compromise and optimal compromise in the sense of Definitions 2.4 and 2.5.

Proposition 2.8**.**

A Gibbs compromise between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ is also a compromise between this sequence of Markov kernels in the sense of Definition 2.4.

The denomination “Gibbs compromise” is justified because it is a stationary distribution for the Gibbs sampler with random equiprobable scanning order.

The proposition below shows that all compromises are tied to Gibbs compromises.

Proposition 2.9**.**

If a sequence of $(r-1)$-dimensional probability distributions $(m_{\neq i})_{i\in[\![1,r]\!]}$ (each $m_{\neq i}$ being a probability distribution on $\bigotimes_{j\neq i}\mathcal{A}_{j}$) is compatible with the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$, then it is the sequence of $(r-1)$-dimensional distributions of a Gibbs compromise between the Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$.

Remark* (Equivalence relation and convexity).*

Let us say that two compromises are equivalent if they share the same sequence of $(r-1)$-marginal distributions. This is obviously an equivalence relation and every class of equivalence can be represented by a single Gibbs compromise. Moreover, each class of equivalence is a convex subset of the set of all compromises. And for any positive measure $\lambda$ on $\mathcal{A}$, its intersection with the set of all compromises absolutely continuous with respect to $\lambda$ is also convex. Finally, the functional $E_{\lambda}$ is convex over this intersection.

2.3. A theoretical justification of Pseudo-Gibbs sampling

At this stage, all necessary tools are available to derive the two most important results of this theory of compromise between incompatible Markov kernels.

Theorem 2.10**.**

If there exists a unique Gibbs compromise $P_{G}$ between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$, then the following statements holds:

1. (1)

$P_{G}$* is absolutely continuous with respect to any compromise between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$;* 2. (2)

for any positive measure $\lambda$ on $\mathcal{A}$ such that $P_{G}$ is absolutely continuous with respect to $\lambda$, $P_{G}$ is the unique optimal compromise with respect to $\lambda$.

Because of these two properties, we call $P_{G}$ the optimal compromise.

Remark* (Equivalence class and convexity – continued).*

The arguments used in the proof of Theorem 2.10 (cf. Appendix A) can also be used to deal with the case where there exist several different Gibbs compromises. First, one can show that each Gibbs compromise is absolutely continuous with respect to any equivalent compromise. Now let $\mathcal{C}$ be an equivalence class and let $\pi_{\mathcal{C}}$ be the unique Gibbs compromise in this class of equivalence. Let $\lambda$ be a positive measure on $\mathcal{A}$ such that $\pi_{\mathcal{C}}$ is absolutely continuous with respect to $\lambda$. $\pi_{C}$ (uniquely) minimizes $E_{\lambda}$ over the intersection of $\mathcal{C}$ and the set of all compromises absolutely continuous with respect to $\lambda$. This implies that for any positive measure $\lambda$ on $\mathcal{A}$, any optimal compromise with respect to $\lambda$ is a Gibbs compromise. Finally, note that the set of all Gibbs compromises is convex (however there is no reason $E_{\lambda}$ should be convex over this set!).

Theorem 2.10 has important practical implications. It opens the possibility of using Gibbs sampling to find the optimal compromise between incompatible Markov kernels and justifies using PIGS.

The next result shows that, under the conditions of Theorem 2.10, the optimal compromise remains the same for a fairly large class of reparametrizations. This result is key to the application of PIGS in an Objective Bayesian framework, where some degree of invariance by reparametrization of priors and posteriors is usually expected.

For every $i\in[\![1,r]\!]$, let $(\tilde{\Omega}_{i},\tilde{\mathcal{A}}_{i})$ be a measurable space and let $f_{i}$ be bijective measurable mapping $\Omega_{i}\rightarrow\tilde{\Omega}_{i}$ whose inverse $f_{i}^{-1}$ is also measurable. Define $f=(f_{1},...,f_{r}):\bigtimes_{i\in[\![1,r]\!]}\Omega_{i}\rightarrow\bigtimes_{i\in[\![1,r]\!]}\tilde{\Omega}_{i}$ and for every $i\in[\![1,r]\!]$ $f_{-i}=(f_{1},..,f_{i-1},f_{i+1},...,f_{r}):\bigtimes_{j\neq i}\Omega_{j}\rightarrow\bigtimes_{j\neq i}\tilde{\Omega}_{j}$.

Also let $\tilde{\pi}_{i}$ be the Markov kernel $\left(\bigtimes_{j\neq i}\tilde{\Omega}_{j}\right)\times\tilde{\mathcal{A}}_{i}\rightarrow[0,1]$ such that for every $\omega_{-i}\in\bigtimes_{j\neq i}\Omega_{j}$ and every $S_{i}\in\mathcal{A}_{i}$, $\tilde{\pi}_{i}(f_{i}(S_{i})|f_{-i}(\omega_{-i}))=\pi_{i}(S_{i}|\omega_{-i})$.

Proposition 2.11**.**

Assume there exists a unique Gibbs compromise $P_{G}$ between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$. Then the push-forward measure of $P_{G}$ by $f$ $\tilde{P}_{G}:=P_{G}\ast f$ is the unique Gibbs compromise between the sequence of Markov kernels $(\tilde{\pi}_{i})_{i\in[\![1,r]\!]}$.

3. Testing the definitions of compatibility

While the previous section presented the theory of compromise, this section aims to test its foundations, namely Definitions 2.3 and 2.4. Subsection 3.1 shows that strengthening their requirements is not possible since it would in many cases threaten the very existence of a compromise. Subsection 3.2 on the other hand shows that these requirements cannot be weakened at a small price.

3.1. Stronger definitions of compromises are not possible

While the definition of the optimal compromise is straightforward, as it involves minimizing some measure of distance between the “targeted” conditionals and the conditionals of the compromise, the definition of a compromise may seem arbitrary. To motivate this definition, let us focus on the two-dimensional case.

Suppose that $r=2$ and that $\pi_{1}$ and $\pi_{2}$ are incompatible. This means there exists no joint distribution $\pi$ which agrees with both Markov kernels. This being the case, it seems sensible to weaken the definition of compatibility by applying it to the “marginals” instead of the “joint” distribution. The following definition makes this idea precise.

Definition 3.1**.**

A pair of probability distributions $m_{\neq 1}$ (resp. $m_{\neq 2}$) on $\mathcal{A}_{2}$ (resp. $\mathcal{A}_{1}$) is compatible with the pair of Markov kernels $\pi_{1}$ and $\pi_{2}$ if the distributions $\pi_{1}m_{\neq 1}$ and $\pi_{2}m_{\neq 2}$ verify

[TABLE]

While this definition may seem more restrictive at first glance than Definition 2.3, both definitions are in fact equivalent when applied to a pair of Markov kernels, because $(r-1)$-dimensional distributions are simply $1$-dimensional distributions in this case. Indeed, following directly from Definition 2.3, we have this result which holds for any $r$:

Proposition 3.2**.**

If a sequence of $(r-1)$-dimensional probability distributions $(m_{\neq i})_{i\in[\![1,r]\!]}$ (each $m_{\neq i}$ being a probability distribution on $\bigotimes_{j\neq i}\mathcal{A}_{j}$) is compatible (in the sense of Definition 2.3) with the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$, then all joint distributions in the sequence $(\pi_{i}m_{\neq i})_{i\in[\![1,r]\!]}$ share the same marginals, that is

[TABLE]

Now let us consider the three-dimensional case ($r=3$). Because the aim of this section is merely to motivate the definitions of compromises and optimal compromises, there is no need for the discussion to be fully general. Let us therefore restrict the discussion to an important particular case. Assume that $\Omega_{1}$, $\Omega_{2}$ and $\Omega_{3}$ are finite sets and that $\mathcal{A}_{1}$, $\mathcal{A}_{2}$ and $\mathcal{A}_{3}$ are respectively the sets of all their subsets. This has an important consequence: any mapping from a subset of $\Omega$ to a subset of $\Omega$ is measurable.

Also assume that the Markov kernels $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$ are positive mappings. For $\pi_{1}$, this means that for every $(\omega_{1},\omega_{2},\omega_{3})\in\Omega_{1}\times\Omega_{2}\times\Omega_{3}$, $\pi_{1}(\{\omega_{1}\}|\omega_{2},\omega_{3})>0$.

If we consider $\omega_{3}$ known, then the situation is reduced to the two-dimensional case. Because $\pi_{1}$ and $\pi_{2}$ are positive mappings and both $\Omega_{1}$ and $\Omega_{2}$ are finite sets, Markov chain theory ensures there exists a unique Gibbs compromise $P(\cdot|\omega_{3})$. For every $(\omega_{1},\omega_{2})\in\Omega_{1}\times\Omega_{2}$,

[TABLE]

Thanks to Theorem 2.10, $P(\cdot|\omega_{3})$ is the optimal compromise. Moreover, notice that Equation (3.3) defines a Markov kernel on $\Omega_{3}\times\left(\mathcal{A}_{1}\otimes\mathcal{A}_{2}\right)$.

We may similarly derive Markov kernels $Q:\Omega_{1}\times\left(\mathcal{A}_{2}\otimes\mathcal{A}_{3}\right)$ and $R:\Omega_{2}\times\left(\mathcal{A}_{1}\otimes\mathcal{A}_{3}\right)$.

Once again, because $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$ are positive mappings, it follows from Markov chain theory that $P$, $Q$ and $R$ are also positive mappings. We now show it using only elementary arguments, because these arguments will be useful again later.

Assume $P$ is not a positive mapping. Then there exists $(\omega_{1}^{(0)},\omega_{2}^{(0)},\omega_{3}^{(0)})\in\Omega_{1}\times\Omega_{2}\times\Omega_{2}$ such that $P(\{\omega_{1}^{(0)}\}\times\{\omega_{2}^{(0)}\}|\omega_{3}^{(0)})=0$. Equation (3.3) then implies that $P(\Omega_{1}\times\{\omega_{2}^{(0)}\}|\omega_{3}^{(0)})=0$. So for every $\omega_{1}\in\Omega_{1}$, $P(\{\omega_{1}\}\times\{\omega_{2}^{(0)}\}|\omega_{3}^{(0)})=0$. But then Equation (3.3) implies that $P(\{\omega_{1}\}\times\Omega_{2}|\omega_{3}^{(0)})=0$. Since this holds for every $\omega_{1}\in\Omega_{1}$, $P(\Omega_{1}\times\Omega_{2}|\omega_{3}^{(0)})=0$, which is absurd since $P(\cdot|\omega_{3}^{(0)})$ is supposed to be a probability distribution. So $P$ is a positive mapping.

Ideally, we would wish to define the optimal compromise between $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$ as the joint distribution $\phi$ such that for every $(\omega_{1},\omega_{2},\omega_{3})\in\Omega_{1}\times\Omega_{2}\times\Omega_{3}$

[TABLE]

Unfortunately, the existence of such an “optimal compromise” $\phi$ implies that $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$ are compatible. First, one can show that for every $(\omega_{1},\omega_{2},\omega_{3})\in\Omega_{1}\times\Omega_{2}\times\Omega_{3}$, $\phi(\{\omega_{1}\}\times\{\omega_{2}\}\times\{\omega_{3}\})>0$. The arguments are similar to those used above to show that $P$ is a positive mapping. Rewriting Equations (3.4) and (3.6) yields

[TABLE]

Now combine Equations (3.7) and (3.8):

[TABLE]

As this holds for every $(\omega_{1},\omega_{2},\omega_{3})\in\Omega_{1}\times\Omega_{2}\times\Omega_{3}$, it implies that $\pi_{2}\phi_{-2}=\pi_{3}\phi_{-3}$. A similar proof then shows that $\pi_{3}\phi_{-3}=\pi_{1}\phi_{-1}$. This means that if an “optimal compromise” $\phi$ exists, then $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$ are compatible and no compromise was needed.

Similarly to what was done in the two-dimensional case, we avoid this difficulty by weakening the compatibility requirements: we no longer require $P(\cdot|\omega_{3})$ to be the optimal compromise between $\pi_{1}$ and $\pi_{2}$ for every $\omega_{3}\in\Omega_{3}$ (as expressed by Equation (3.3)), but only on average over $\omega_{3}\in\Omega_{3}$. So a compromise $\phi$ should still verify Equation (3.4), but it would only need to verify this weakened version of Equation (3.3) for every $(\omega_{1},\omega_{2})\in\Omega_{1}\times\Omega_{2}$:

[TABLE]

Because Equation (3.4) is still expected to hold, Equation (3.10) is equivalent to

[TABLE]

So the requirement boils down to $2\phi_{-3}=(\pi_{1}\phi_{-1})_{-3}+(\pi_{2}\phi_{-2})_{-3}$. Of course, we symmetrically require $2\phi_{-1}=(\pi_{2}\phi_{-2})_{-1}+(\pi_{3}\phi_{-3})_{-3}$ and $2\phi_{-2}=(\pi_{1}\phi_{-1})_{-2}+(\pi_{3}\phi_{-3})_{-2}$ as well.

To sum this part of the discussion up, the necessity of weakening our “ideal” requirements for “optimal compatibility” made us downgrade from a requirement about a “joint” $3$-dimensional distribution $\phi$ to requirements about its $(3-1)$-marginal distributions $\phi_{-1}$, $\phi_{-2}$ and $\phi_{-3}$. Definition 2.3 is just another formulation of these requirements, which are taken to define a compatible sequence of $(r-1)$-dimensional distributions.

Proposition 2.9 shows that in cases where there exists a unique Gibbs compromise, even with this weakened set of requirements for compatibility, there exists only one compatible sequence of $(r-1)$-dimensional distributions, so we may not strenghten it if there is to be any solution. Indeed, in cases where no Gibbs compromise exists, no compatible set of $(r-1)$-dimensional distributions exists either!

3.2. Weaker definitions of compromises are inconvenient

As was shown in the previous subsection, the requirements given by Definition 2.3 for the compatibility of a sequence of $(r-1)$-dimensional distributions with a given sequence of Markov kernels cannot be strengthened. The following shows they cannot be weakened either.

3.2.1. Why we need some notion of compatibility: example in 2 dimensions.

Why bother with the compatibility of $(r-1)$-dimensional distributions and not simply minimize the functional $E_{\lambda}$ of Definition 2.5 over all distributions absolutely continuous with respect to $\lambda$ ? The following 2-dimensional example shows that doing so yields unsatisfactory results.

Consider the following situation: $\Omega_{1}=\Omega_{2}=\{0,1\}$ and $\mathcal{A}_{1}=\mathcal{A}_{2}=\{\emptyset,\{0\},\{1\},\{0,1\}\}$. Let $X_{1}$ (resp. $X_{2}$) be the identity function on $\Omega_{1}$ (resp. $\Omega_{2}$). Both $X_{1}$ and $X_{2}$ are measurable functions (and thus random variables when $\mathcal{A}_{1}\otimes\mathcal{A}_{2}$ is endowed with a probability measure). Define the following Markov kernels:

[TABLE]

Let $\lambda$ be the counting measure. Denote by $\pi_{C}$ the optimal compromise with respect to $\lambda$ (minimizing $E_{\lambda}$ over all compromises) and $\pi_{E}$ the distribution on $\mathcal{A}_{1}\otimes\mathcal{A}_{2}$ that minimizes $E_{\lambda}$ over all distributions on $\mathcal{A}_{1}\otimes\mathcal{A}_{2}$. We have $E_{\lambda}(\pi_{C})=2/25>E_{\lambda}(\pi_{E})=1/15$ and

[TABLE]

As $\pi_{C}$ is the unique Gibbs compromise between $\pi_{1}$ and $\pi_{2}$, its marginals are the only compatible marginals: $(\pi_{C})_{-1}(X_{2}=1)=4/5$ and $(\pi_{C})_{-2}(X_{1}=1)=3/5$. The marginals of $\pi_{E}$ are noticeably different: $(\pi_{E})_{-1}(X_{2}=1)=13/15$ and $(\pi_{E})_{-2}(X_{1}=1)=8/15$.

Observe that according to $\pi_{1}$, $X_{2}=0$ implies $X_{1}=1$ but that according to $\pi_{2}$, $X_{1}=1$ implies that $X_{2}=1\neq 0$. This discrepancy is a major source of incompatibility between the two Markov kernels. So, as $\pi_{E}$ makes both $X_{1}=1$ and $X_{2}=0$ less likely than $\pi_{C}$, it “ignores the inconsistent parts” of $\pi_{1}$ and $\pi_{2}$ to some extent. Therefore, if the marginals are not set in advance (say, by imposing compatibility with the conditionals in the sense of Definition 2.3), one may “cheat” by having the marginals disadvantage inconvenient values for the parameters.

3.2.2. Why the compatibility requirements can hardly be weakened: example in 3 dimensions.

In the two-dimensional case, because of Proposition 3.2, Definitions 3.1 and 2.3 give the same meaning to the concept of compatibility of marginals, so Definition 2.3 may be thought of as a generalization of Definition 3.1 to cases with more than two dimensions. However, another generalization of the latter definition is possible. To avoid confusion, this other generalization will be called weak compatibility. In the following, the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ is defined as in Section 2.1.

Definition 3.3**.**

A sequence of $(r-1)$-dimensional distributions $(m_{\neq i})_{i\in[\![1,r]\!]}$ is weakly compatible with a sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ if Equation (3.2) holds.

Proposition 3.2 means that compatibility in the sense of Definition 2.3 implies weak compatibility in the sense of Definition 3.3, hence its denomination as “weak”.

Using the concept of weak compatibility of a sequence of $(r-1)$-marginal distributions, we define weak compromises and the optimal weak compromise as analogues to compromises and optimal compromises respectively.

Definition 3.4**.**

A probability distribution $P$ on $\bigotimes_{i\in[\![1,r]\!]}\mathcal{A}_{i}$ is called a weak compromise between the Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ if these two conditions are verified:

1. (1)

for every $i\in[\![1,r]\!]$, $\pi_{i}P_{-i}$ is absolutely continuous with respect to $P$; 2. (2)

the sequence $(P_{-i})_{i\in[\![1,r]\!]}$ of $P$’s $(r-1)$-marginal distributions is weakly compatible with $(\pi_{i})_{i\in[\![1,r]\!]}$.

Definition 3.5**.**

Let $\lambda$ be a positive measure on $\mathcal{A}$. Let $P$ be a weak compromise between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ that is absolutely continuous with respect to $\lambda$. $P$ is called an optimal weak compromise with respect to $\lambda$ between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$ if it minimizes the functional $E_{\lambda}$ over all compromises between $(\pi_{i})_{i\in[\![1,r]\!]}$ that are absolutely continuous with respect to $\lambda$. $E_{\lambda}$ is defined by Equation (2.13).

Because, for any positive measure $\lambda$ on $\mathcal{A}$, the set of all weak compromises absolutely continuous with respect to $\lambda$ includes the set of all compromises absolutely continuous with respect to $\lambda$, an optimal weak compromise with respect to $\lambda$ makes the functional $E_{\lambda}$ no greater than an optimal compromise with respect to $\lambda$. However, as shown in the following example with $r=3$, optimal weak compromises may have undesirable behavior.

Assume $\Omega_{1}=\Omega_{2}=\Omega_{3}=\{0,1\}$ and $\mathcal{A}_{1}=\mathcal{A}_{2}=\mathcal{A}_{3}=\{\emptyset,\{0\},\{1\},\{0,1\}\}$. Let $X_{1}$ (resp. $X_{2}$, $X_{3}$) be the identity function on $\Omega_{1}$ (resp. $\Omega_{2}$, $\Omega_{3}$).

Consider the following Markov kernels. For every $(\omega_{1},\omega_{2},\omega_{3})\in\Omega_{1}\times\Omega_{2}\times\Omega_{3}$

[TABLE]

Notice that, provided $\omega_{3}$ is known, $\pi_{1}$ and $\pi_{2}$ are compatible Markov kernels. The unique probability distribution on $\mathcal{A}_{1}\otimes\mathcal{A}_{2}$ that fits both $\pi_{1}$ and $\pi_{2}$ (conditional to $\omega_{3}$) verifies for every $(\omega_{1},\omega_{2})\in\Omega_{1}\times\Omega_{2}$

[TABLE]

Thus, any joint distribution fitting the Markov kernel $P$ would make $X_{1}$, $X_{2}$ and $X_{3}$ mutually independent. Unfortunately, no such joint distribution could fit $\pi_{3}$, but we may expect compromises between $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$ to retain the independence of $X_{1}$ and $X_{2}$.

Let $\lambda$ be the counting measure on $\mathcal{A}$.

Denote by $\pi_{C}$ the (unique) optimal compromise between $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$ with respect to $\lambda$: we have $E_{\lambda}(\pi_{C})=1/48\approx 0.021$. For every $(\omega_{1},\omega_{2},\omega_{3})\in\Omega_{1}\times\Omega_{2}\times\Omega_{3}$

[TABLE]

Notably, its third 2-marginal distribution $(\pi_{C})_{-3}$ verifies for every $(\omega_{1},\omega_{2})\in\Omega_{1}\times\Omega_{2}$

[TABLE]

So, as expected, $\pi_{C}$ retains the independence of $X_{1}$ and $X_{2}$. Because $\pi_{C}$ is the unique Gibbs compromise between $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$, Proposition 2.9 implies any other compromise between $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$ also retains this property.

Let us now consider an optimal weak compromise $\pi_{W}$ between $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$ with respect to $\lambda$. Numerical computation gives us the following approximation, with $E_{\lambda}(\pi_{W})\approx 0.019$. For every $(\omega_{1},\omega_{2},\omega_{3})\in\Omega_{1}\times\Omega_{2}\times\Omega_{3}$,

[TABLE]

Its third 2-marginal distribution $(\pi_{W})_{-3}$ is approximately

[TABLE]

Thus the independence of $X_{1}$ and $X_{2}$ is lost. Therefore, weak compatibility is no adequate notion of compatibility. As a matter of fact, although $(\pi_{W})_{-3}$ and $(\pi_{C})_{-3}$ share the same marginals, that is

[TABLE]

$(\pi_{W})_{-3}$ slightly disadvantages the event $X_{1}=X_{2}$, which is where the incompatibility between $\pi_{3}$ and the pair $(\pi_{1},\pi_{2})$ is most obvious: according to $\pi_{3}$, $X_{1}=X_{2}$ implies $X_{3}=1$, so conversely, $X_{3}=0$ should imply $X_{1}\neq X_{2}$, when in fact $\pi_{1}$ and $\pi_{2}$ state that even given $\omega_{3}=0$, $\{X_{1}\neq X_{2}\}$ only happens with probability 1/2. On the other hand, according to $\pi_{3}$, if $X_{1}\neq X_{2}$, then $X_{3}$ can with equal probability be 0 or 1, which matches $\pi_{1}$ and $\pi_{2}$ better.

4. Optimal compromise between Objective Posterior conditional distributions in Gaussian Process regression

Kriging is a surrogate model used to emulate a real-valued function on a spatial domain $\mathcal{D}$ when said function can only be evaluated on a finite subset of $\mathcal{D}$ called “design set”. The “Kriging prediction” is the mean function of the process taken conditionally to all known values of the emulated function, i.e. the values at the points in the design set. The main advantage of the framework is its natural way of representing uncertainty about the value of the function at unobserved points [31]. Prediction does not consist of a single value but of a complete Normal distribution. “Kriging” is the name given to the framework in the geostatistical literature [18], but is also frequently used in the context of computer experiments and machine learning under the label “Gaussian Process regression” [26]. In this work, we focus on Simple Kriging, where the Gaussian Process is assumed to be stationary with known mean, as opposed to Universal Kriging, which incorporates an unknown mean function.

The probability distribution of a stationary Gaussian Process is characterized by a variance parameter and a correlation function (also known as “correlation kernel”) which itself depends on parameters. So one should deal with uncertainty about model parameters.

The problem is “notoriously difficult”, as highlighted by Kennedy and O’Hagan [20], because the likelihood function may often be quite flat [23]. In a Bayesian framework, this uncertainty is represented by a prior distribution on the parameters.

4.1. Issues raised by objective prior elicitation for Gaussian processes

Let $Y(\boldsymbol{x})$, $\boldsymbol{x}\in\mathcal{D}$ be a real-valued random field on a bounded subset $\mathcal{D}$ of $\mathbb{R}^{r}$. We assume $Y$ is Gaussian with zero mean (or known mean) and with covariance of the form ${\rm Cov}(Y(\boldsymbol{x}),Y(\boldsymbol{x}^{\prime}))=\sigma^{2}K_{\boldsymbol{\theta}}(\boldsymbol{x}-\boldsymbol{x}^{\prime})$. $\sigma^{2}$ thus denotes the variance of the Gaussian Process and $\boldsymbol{\theta}\in(0,+\infty)^{r}$, hereafter named the “vector of correlation lengths”, is the vector of scaling parameters used by the chosen class of correlation kernels $K_{\boldsymbol{\theta}}$.

Consider a set of $n\in\mathbb{N}$ points $(\boldsymbol{x}^{(i)})_{i\in[\![1,n]\!]}$ belonging to the domain $\mathcal{D}$. This set is called the design set and $Y$ is observed at all points of this set. Let $\boldsymbol{Y}$ be the Gaussian vector $(Y(\boldsymbol{x}^{(i)}))_{i\in[\![1,n]\!]}$ and let $\boldsymbol{\Sigma}_{\boldsymbol{\theta}}$ be its correlation matrix: the distribution of $\boldsymbol{Y}$ is therefore $\mathcal{N}(\boldsymbol{0}_{n},\sigma^{2}\boldsymbol{\Sigma}_{\boldsymbol{\theta}})$.

Let $\boldsymbol{y}$ be the vector of observations. When applied to a matrix, $|\cdot|$ refers to its determinant.

With these notations, the likelihood of the parameters $\sigma^{2}$ and $\boldsymbol{\theta}$ is

[TABLE]

The reference prior with parameter ordering $\sigma^{2}\prec\boldsymbol{\theta}$ is given by Berger et al. [6]:

[TABLE]

The distribution $\pi(\sigma^{2}|\boldsymbol{\theta})$ has infinite mass: it is an improper prior.

Let us integrate $\sigma^{2}$ out of the likelihood (4.1):

[TABLE]

Ren et al. [27] provide the reference prior $\pi(\boldsymbol{\theta})$, where $\boldsymbol{\theta}$ is regarded as a single multidimensional parameter. It is proportional to the square root of the determinant of the $r\times r$ matrix with $(i,j)$-th entry

[TABLE]

However, this method has the disadvantage of requiring the use of a multidimensional Jeffreys-rule prior distribution, which may show the sort of undesirable behavior mentioned in the introduction.

Alternatively, we could draw inspiration from the one-dimensional case in the following way. Suppose that we know every entry of $\boldsymbol{\theta}$ except one, $\theta_{i}$. Then, according to Equation (4.4), the prior density on $\theta_{i}$ knowing all entries $\theta_{j}$ ($j\neq i$) would be

[TABLE]

The density functions $\pi_{i}(\theta_{i}\;|\;\theta_{j}\;\forall j\neq i)$ ($i\in[\![1,r]\!]$) define Markov kernels. Indeed, they are continuous with respect to the $\theta_{j}$ ($j\neq i$) and are probability densities with respect to the Lebesgue measure. They are unfortunately likely to violate the necessary condition for compatibility given by Equation (2.6).

Let us now consider the corresponding posterior conditional densities $\pi_{i}(\theta_{i}\;|\;\boldsymbol{y},\;\theta_{j}\;\forall j\neq i)$ ($i\in[\![1,r]\!]$). Just like their prior counterparts, they are likely to violate the necessary condition for compatibility. However, each of them represents our opinion about one parameter if all others were known. This is a setting where the results of Section 2 can be applied in order to find the optimal compromise between the Markov kernels $\mathbb{R}^{r-1}\times\mathcal{B}(\mathbb{R})$ they define. This optimal compromise will then be taken as posterior probability of the vector $\boldsymbol{\theta}$. In the following, we describe settings in which there exists a single Gibbs compromise between these Markov kernels. Theorem 2.10 then asserts it is the optimal compromise. We call this compromise the Gibbs reference posterior distribution because of its link to the reference posterior distribution in settings with a one-dimensional parameter $\boldsymbol{\theta}$.

However, even though we call it a “posterior” distribution, it is unclear whether there exists a prior distribution from which it could be derived using Bayes’ rule. Denote by $\pi_{G}(\boldsymbol{\theta}|\boldsymbol{y})$ the probability density with respect to the Lebesgue measure of the Gibbs reference posterior distribution. Bayes’ rule requires that in case a (proper or improper) prior density $\pi_{G}(\boldsymbol{\theta})$ exists, there should also exist a function $\tilde{L}(\boldsymbol{y})$ such that, for almost every $\boldsymbol{\theta}\in\mathbb{R}^{r}$ in the sense of the Lebesgue measure,

[TABLE]

As we have no explicit expression of $\pi_{G}(\boldsymbol{\theta}|\boldsymbol{y})$, we have no way to check whether Equation (4.6) holds or not.

In this section, we establish that, whenever Matérn anisotropic geometric or tensorized kernels with known smoothness parameter $\nu$ are used, under certain conditions to be detailed later, there exists a unique Gibbs compromise between the reference posterior conditionals, which thanks to Theorem 2.10 is the optimal compromise. Henceforth, it will be called “Gibbs reference posterior distribution”, even though this “posterior” has not been derived from a prior distribution using the Bayes rule.

All proofs for this section can be found in Appendix B.

4.2. Definitions

In this work, we use the following convention for the Fourier transform: the Fourier transform $\widehat{g}$ of a smooth function $g:\mathbb{R}^{r}\rightarrow\mathbb{R}$ verifies $g(\boldsymbol{x})=\int_{\mathbb{R}^{r}}\widehat{g}(\boldsymbol{\omega})e^{i\langle\boldsymbol{\omega}|\boldsymbol{x}\rangle}d\boldsymbol{\omega}$ and $\widehat{g}(\boldsymbol{\omega})=(2\pi)^{-r}\int_{\mathbb{R}^{r}}g(\boldsymbol{x})e^{-i\langle\boldsymbol{\omega}|\boldsymbol{x}\rangle}d\boldsymbol{x}$.

Let us set up a few notations.

1. (a)

$\mathcal{K}_{\nu}$ is the modified Bessel function of second kind with parameter $\nu$ ; 2. (b)

$K_{r,\nu}$ is the $r$-dimensional Matérn isotropic covariance kernel with variance 1, correlation length 1 and smoothness $\nu\in(0,+\infty)$ and $\widehat{K}_{r,\nu}$ is its Fourier transform:

1. (i)

$\forall\boldsymbol{x}\in\mathbb{R}^{r}$,

[TABLE] 2. (ii)

$\forall\boldsymbol{\omega}\in\mathbb{R}^{r},$

[TABLE] 3. (c)

$K_{r,\nu}^{tens}$ is the $r$-dimensional Matérn tensorized covariance kernel with variance 1, correlation length 1 and smoothness $\nu\in\mathbb{R}_{+}$ and $\widehat{K}_{r,\nu}^{tens}$ is its Fourier transform:

1. (i)

$\forall\boldsymbol{x}\in\mathbb{R}^{r}$,

[TABLE] 2. (ii)

$\forall\boldsymbol{\omega}\in\mathbb{R}^{r}$,

[TABLE] 4. (d)

if $\boldsymbol{t}\in\mathbb{R}^{r}$, $\frac{\boldsymbol{t}}{\boldsymbol{\theta}}=\left(\frac{t_{1}}{\theta_{1}},...,\frac{t_{r}}{\theta_{r}}\right)$ and $\boldsymbol{t}\,\boldsymbol{\mu}=\left(t_{1}\mu_{1},...,t_{r}\mu_{r}\right)$.

We define the Matérn geometric anisotropic covariance kernel with variance parameter $\sigma^{2}$, correlation lengths $\boldsymbol{\theta}$ (resp. inverse correlation lengths $\boldsymbol{\mu}$) and smoothness $\nu$ as the function $\boldsymbol{x}\mapsto\sigma^{2}K_{r,\nu}\left(\frac{\boldsymbol{x}}{\boldsymbol{\theta}}\right)$ (resp. $\boldsymbol{x}\mapsto\sigma^{2}K_{r,\nu}\left(\boldsymbol{x}\boldsymbol{\mu}\right)$).

Similarly, we define the Matérn tensorized covariance kernel with variance parameter $\sigma^{2}$, correlation lengths $\boldsymbol{\theta}$ (resp. inverse correlation lengths $\boldsymbol{\mu}$) and smoothness $\nu$ as the function $\boldsymbol{x}\mapsto\sigma^{2}K_{r,\nu}^{tens}\left(\frac{\boldsymbol{x}}{\boldsymbol{\theta}}\right)$ (resp. $\boldsymbol{x}\mapsto\sigma^{2}K_{r,\nu}^{tens}\left(\boldsymbol{x}\boldsymbol{\mu}\right)$).

Thanks to Proposition 2.11, we may choose any parametrization we wish for the Matérn correlation kernels. We have found that the parametrization involving inverse correlation lengths makes proofs easier.

Several key passages in the proofs (to be found in Appendix B) involve a technical assumption on the design set:

Definition 4.1**.**

A design set is said to have coordinate-distinct points, or simply to be coordinate-distinct, if for any distinct points in the set $\boldsymbol{x}$ and $\boldsymbol{x^{\prime}}$, every component of the vector $\boldsymbol{x}-\boldsymbol{x^{\prime}}$ differs from 0.

Most randomly sampled design sets almost surely have coordinate-distinct points – for instance Latin Hypercube Sampling. Cartesian product design sets, however, do not.

4.3. Main result

The result is valid for Simple Kriging models with the following characteristics:

1. (a)

the design set contains $n$ coordinate-distinct points in $\mathbb{R}^{r}$ ($n$ and $r$ are positive integers); 2. (b)

the covariance function is Matérn anisostropic geometric or tensorized with variance parameter $\sigma^{2}>0$, smoothness parameter $\nu$ and vector of correlation lengths (resp. inverse correlation lengths) $\boldsymbol{\theta}\in(0,+\infty)^{r}$ (resp. $\boldsymbol{\mu}\in(0,+\infty)^{r}$) ; 3. (c)

one of the following conditions is verified:

1. (i)

$\nu\in(0,1)$ and $n>1$ and the Matérn kernel is tensorized ; 2. (ii)

$\nu\in(1,2)$ and $n>r+2$ ; 3. (iii)

$\nu\in(2,3)$ and $n>r(r+1)/2+2r+3$.

Theorem 4.2**.**

In a Simple Kriging model with the characteristics described above, there exists a hyperplane $\mathcal{H}$ of $\mathbb{R}^{n}$ such that, for any $\boldsymbol{y}\in\mathbb{R}^{n}\setminus\mathcal{H}$, there exists a unique Gibbs compromise $\pi_{G}(\boldsymbol{\theta}|\boldsymbol{y})$ (resp. $\pi_{G}(\boldsymbol{\mu}|\boldsymbol{y})$) between the reference posterior conditionals $\pi_{i}(\theta_{i}|\boldsymbol{y},\boldsymbol{\theta}_{-i})$ (resp. $\pi_{i}(\mu_{i}|\boldsymbol{y},\boldsymbol{\mu}_{-i})$). It is the unique stationary distribution of the Markov kernel $P_{\boldsymbol{y}}:(0,+\infty)^{r}\times\mathcal{B}\left((0,+\infty)^{r}\right)\rightarrow[0,1]$ defined by

[TABLE]

The Markov kernel $P_{\boldsymbol{y}}$ is uniformly ergodic: $\lim_{n\to\infty}\sup_{\boldsymbol{\theta}^{(0)}\in(0,+\infty)^{r}}\|P_{\boldsymbol{y}}^{n}(\boldsymbol{\theta}^{(0)},\cdot)-\pi_{G}(\cdot|\boldsymbol{y})\|_{TV}=0$ (resp. $\lim_{n\to\infty}\sup_{\boldsymbol{\mu}^{(0)}\in(0,+\infty)^{r}}\|P_{\boldsymbol{y}}^{n}(\boldsymbol{\mu}^{(0)},\cdot)-\pi_{G}(\cdot|\boldsymbol{y})\|_{TV}=0$), where $\|\cdot\|_{TV}$ is the total variation norm.

Remark*.*

The reference posterior conditionals are invariant by reparametrization, so the Markov kernel $P_{\boldsymbol{y}}$ does not depend on whether the chosen parametrization uses correlation lengths $\boldsymbol{\theta}$ or inverse correlation lengths $\boldsymbol{\mu}$. Due to Proposition 2.11, the Gibbs compromise does not either. The parametrization using inverse correlation lenghts $\boldsymbol{\mu}$ is more convenient for proving this theorem, however.

Notice that in such a Kriging model, the vector of observations $\boldsymbol{y}$ almost surely belongs to $\mathbb{R}^{n}\setminus\mathcal{H}$, so this assumption is of no practical consequence. Theorem 4.2 therefore asserts that the Gibbs compromise between the incompatible conditionals $\pi_{i}(\mu_{i}|\boldsymbol{y},\boldsymbol{\mu}_{-i})$ exists, is unique, and can be sampled from using Potentially Incompatible Gibbs Sampling (PIGS). In the following, it is called “Gibbs reference posterior distribution”.

4.4. Using the Gibbs reference posterior distribution

Let $\boldsymbol{x}_{0}$ be a point in the domain $\mathcal{D}$ that does not belong to the design set. Denote by $\boldsymbol{\Sigma}_{\boldsymbol{\theta},0,\cdot}$ the correlation matrix between $Y(\boldsymbol{x}_{0})$ and $\boldsymbol{Y}$, and by $\boldsymbol{\Sigma}_{\boldsymbol{\theta},\cdot,0}$ its transpose the correlation matrix between $\boldsymbol{Y}$ and $Y(\boldsymbol{x}_{0})$.

Theorem 4.1.2. (case 4) of Santner et al. [31] provides this useful result for prediction:

Proposition 4.3**.**

Conditionally to $\boldsymbol{Y}=\boldsymbol{y}$ and assuming $\boldsymbol{\theta}$ is known, the random variable $Z_{0}$ defined below follows the Student t-distribution with $n$ degrees of freedom.

[TABLE]

Remark*.*

If $n$ exceeds 30, it is usually accepted that the Student t-distribution with $n$ degrees of freedom can be approximated by the standard Normal distribution. As this threshold should be exceeded in practical cases, we would recommend performing all computations as though the Student t-distribution were Normal.

This Proposition implicitely gives the distribution of $y_{0}:=Y(\boldsymbol{x}_{0})$ conditionally to $\boldsymbol{Y}=\boldsymbol{y}$ and $\boldsymbol{\theta}$. For later reference, denote it by $L^{1}(y_{0}|\boldsymbol{y},\boldsymbol{\theta})$. In practice, when $\boldsymbol{\theta}$ is unknown, the distribution of $y_{0}=Y(\boldsymbol{x}_{0})$ conditionally to $\boldsymbol{Y}=\boldsymbol{y}$ can be obtained once $\boldsymbol{\theta}$ has been sampled from the Gibbs reference posterior distribution:

[TABLE]

Its cdf can be approximated by averaging the cdfs of the Student t-distributions (or their Normal approximations) corresponding to every point in the sample.

5. Comparisons between the MLE and MAP estimators

To illustrate the inferential performance of the Gibbs reference posterior distribution, let us introduce the Maximum A Posteriori estimator (MAP). It takes the value $\hat{\boldsymbol{\theta}}_{MAP}$ of $\boldsymbol{\theta}$ where the density with respect to the Lebesgue measure of the Gibbs reference posterior distribution is largest. We contrast it with the Maximum Likelihood Estimator (MLE) $\hat{\boldsymbol{\theta}}_{MLE}$ which does the same with the likelihood function.

5.1. Methodology

In this section, we compare the MLE and MAP estimators for accuracy and robustness.

Our test cases are 3-dimensional Gaussian Processes with Matérn anisotropic geometric correlation kernels with smoothness 5/2. Their mean is the null function, which only leaves us with the matter of estimating their correlation length for each dimension.

We use uniform designs: our observation points are randomly generated according to the uniform distribution on a cube with side length 1.

In order to measure the performance of our estimators, we define a suitable distance between two vectors of correlation lengths. Then the error of an estimator is defined as its distance to the “true” vector of correlation lengths.

Let $g$ be the function such that for any $t$ in $(-1,1)$, $g(t)=\operatorname{argtanh}(t)$ and $g(-1)=g(1)=0$. We use the convention that, for any matrix $\boldsymbol{M}$ with elements in $[0,1]$, $g(\boldsymbol{M})$ is the matrix resulting from applying $g$ to every element of $\boldsymbol{M}$.

Definition 5.1**.**

For a given design set, the distance between two vectors of correlation lengths $\boldsymbol{\theta}^{1}$ and $\boldsymbol{\theta}^{2}$ is $\|g(\boldsymbol{\Sigma}_{\boldsymbol{\theta}^{1}})-g(\boldsymbol{\Sigma}_{\boldsymbol{\theta}^{2}})\|$, where $\|\cdot\|$ denotes the Frobenius norm.

This distance involves applying the Fisher transformation [17] (that is, the inverse hyperbolic tangent function) to every (non-unitary) correlation coefficient in both associated correlation matrices. This is a variance-stabilizing transformation. For any random variables $U$ and $V$ following the normal distribution with mean [math] and variance $1$, let $\rho$ denote the correlation coefficient between $U$ and $V$ ($-1<\rho<1$). If $(U_{i},V_{i})$ $(1\leqslant i\leqslant N)$ are independent copies of $(U,V)$, then $\hat{\rho}=\sum_{i=1}^{N}U_{i}V_{i}/n$ is a random variable and $\operatorname{argtanh}(\hat{\rho})$ follows the normal distribution with mean $\operatorname{argtanh}(\rho)$ and variance $1/(N-3)$. So the variance of $\operatorname{argtanh}(\hat{\rho})$ does not depend on $\rho$, whereas the variance of $\hat{\rho}$ does and goes to zero for $|\rho|\rightarrow 1$. Involving the Fisher transformation in the definition of the distance between two vectors of correlation lengths is therefore a way to assert that vectors of correlation lengths can be far apart even if they both lead to highly correlated observations.

This allows us to make sure errors made when estimating near-1 correlation coefficients are no less taken into account than errors made when estimating near-0 correlation coefficients.

Let us choose a “true” vector of correlation lengths (and also a variance parameter, but this parameter has no effect on either the MLE or the MAP). Then we need to:

1. (1)

Sample $n$ points randomly according to the uniform distribution on the unit cube (in the following, $n=30$). 2. (2)

Generate the observations of the Gaussian Process at the sampled points according to the selected “true” variance and correlation lengths. 3. (3)

Sample the vector of correlation lengths according to the Gibbs reference posterior distribution $\pi_{G}(\boldsymbol{\theta}|\boldsymbol{y})$ through PIGS. 4. (4)

Compute the MLE and the MAP of the vector of correlation lengths and their errors. 5. (5)

Repeat steps 1 to 4 $m-1$ times (in the following, $m=500$).

This method allows us to derive an approximate distribution of the errors of both estimators when both the realization of the Gaussian Process and the design set vary. Thus we get to test the robustness of both estimators versus the variability of both the Gaussian Process and the choice of design set.

5.2. Results

This subsection provides results obtained on 3-dimensional Gaussian Processes with null mean function and Matérn anisotropic geometric correlation kernels with smoothness 5/2. The results are divided by “true” vectors of correlation lengths. In each case, we give in Table 1 the empirical Root Mean Square Errors (RMSEs) of both MLE and MAP estimators as functions of varying instances of the Gaussian Process and uniform design sets on the unit cube.

Most of the “true” vectors of correlation lengths featured in Table 1 were selected in a way to showcase the behavior of both estimators in strongly anisotropic cases, but one (0.5 - 0.5 - 0.5) also showcases their behavior if the true kernel is actually isotropic. And the final one (0.8 - 1 - 0.9) is used to illustrate the performance in the case of a strongly correlated Gaussian Process: this case is fundamentally different from all others, because the Matérn anisotropic geometric family of correlation kernels is designed in such a way that the correlation length with greatest influence is the lowest. Informally speaking, it is enough for one correlation length to be near zero to make the whole process very uncorrelated, even should all other correlation lengths be very high.

In all studied cases, the MAP estimator was more robust than the MLE estimator: its RMSE (Root Mean Square Error) was between 9 and 15% lower, as showcased in Table 1.

To get a better sense of the distribution of the error when the design set and the realization of the Gaussian Process vary, we give in Figure 1 violin plots of the errors in the two most extreme case: very low correlation (0.4 – 0.8 – 0.2) and very high correlation (0.8 – 1.0 – 0.9)

6. Comparison of the MLE and MAP plug-in distributions and the full posterior predictive distribution

6.1. Methodology

We use the same test cases as before. In this section, our goal is to assess the accuracy of prediction intervals associated with both estimators and with the full posterior distribution. We consider 95% intervals: the lower bound is the 2.5% quantile and the upper bound the 97.5% quantile of plug-in distributions $\hat{P}_{MLE}(y_{0}|\boldsymbol{y})=L^{1}(y_{0}|\boldsymbol{y},\hat{\boldsymbol{\theta}}_{MLE})$, $\hat{P}_{MAP}(\boldsymbol{y}_{0}|\boldsymbol{y})=L^{1}(y_{0}|\boldsymbol{y},\hat{\boldsymbol{\theta}}_{MAP})$ and $P(\boldsymbol{y}_{0}|\boldsymbol{y})=\int L^{1}(y_{0}|\boldsymbol{y},\boldsymbol{\theta})\;\pi_{G}(\boldsymbol{\theta}|\boldsymbol{y})\;d\boldsymbol{\theta}$. For the sake of comprehensiveness, we also consider predictive intervals associated with the “true” predictive distribution $L(\boldsymbol{y}_{0}\;|\;\boldsymbol{y},\sigma^{2},\boldsymbol{\theta})$, which is the predictive distribution we would use if we knew the correct values of the parameters $\sigma^{2}$ and $\boldsymbol{\theta}$.

Let us choose a “true” vector of correlation lengths $\boldsymbol{\theta}$ (and also a variance parameter $\sigma^{2}$, but this parameter has no effect on predictive accuracy). Then we do the following:

1. (1)

Sample $n$ observation points randomly according to the uniform distribution on the unit cube (in the following, $n=30$). 2. (2)

Generate the observations of the Gaussian Process at the sampled points according to the selected “true” variance and correlation lengths. 3. (3)

Sample the vector of correlation lengths according to the Gibbs reference posterior distribution $\pi_{G}(\boldsymbol{\theta}|\boldsymbol{y})$ through PIGS. 4. (4)

Compute the MLE and the MAP of the vector of correlation lengths. 5. (5)

Sample $n_{0}$ test points randomly according to the uniform distribution on the unit cube (in the following, $n_{0}=100$). 6. (6)

At each point, determine the 95% prediction intervals derived from $L(\boldsymbol{y}_{0}\;|\;\boldsymbol{y},\sigma^{2},\boldsymbol{\theta})$ ($\sigma^{2}$ and $\boldsymbol{\theta}$ being the “true” parameters), $\hat{P}_{MLE}(\boldsymbol{y}_{0}\;|\;\boldsymbol{y})$, $\hat{P}_{MAP}(\boldsymbol{y}_{0}\;|\;\boldsymbol{y})$ and $P(\boldsymbol{y}_{0}\;|\;\boldsymbol{y})$. 7. (7)

Generate the values of the Gaussian Process at the newly sampled points (naturally, do this conditionally to the previously generated observations). 8. (8)

Count the number of points within the prediction intervals derived from each of the four distributions. Divide the counts by $n_{0}$: this yields four coverages corresponding to each type of predictive intervals. Also compute the mean length of every type of prediction interval. 9. (9)

Repeat steps 1 to 8 $m-1$ times (in the following, $m=500$).

6.2. Results

There is no reason for individual coverages of 95% predictive intervals given by the predictive distribution to be equal to 95%. Recall that any coverage is given for a unique realization of the Gaussian Process, and that the values of this process at different points are correlated. If the predictive interval at some point fails to cover the true value at this point, it is likely that predictive intervals at neighboring points will also fail to cover the true values at those points, even though the nominal value is 95% everywhere. Conversely, if it actually covers the true value, then prediction intervals at neighboring points are more than 95% likely to cover their true values.

In short, prediction intervals give information that is only valid if understood to refer to what can be guessed on the sole basis of the observations made at the design points, which is why coverages for individual realizations of the Gaussian Process are not necessarily 95% even if the predictive distribution is perfectly accurate (i.e. based on the true values of $\sigma^{2}$ and $\boldsymbol{\theta}$).

However, if the predictive distribution is perfectly accurate, then the average of the coverages is the nominal value: 95%. It is thus interesting to compute the average of the coverages for all distributions, whether they are plug-in distributions based on the MLE or MAP estimator, or the predictive distribution based on the full posterior distribution (hereafter noted FPD). In the above described methodology, the average was taken over the realizations of the Gaussian Process with the chosen true parameters and over all design sets with $n$ design points. The results below are obtained in this way.

The results given in Table 2 show that using the full posterior distribution (FPD) to derive the predictive distribution is the best possible choice from a frequentist point of view as the nominal value is nearly matched by the average coverage. Predictive intervals derived from the MAP estimator do not perform as well, and predictive intervals derived from the MLE perform even worse.

Let us now focus on the average (with respect to the uniform design sets and realizations of the Gaussian Process) of the mean (over the test set for a given realization of the Gaussian Process and a given uniform design set) length of prediction intervals. The results are given in Table 3, where the figures between parentheses give the increase or decrease (in percents) of the average mean length when compared to the average mean length of prediction intervals obtained using the true values of the parameters.

Predictive intervals derived from the full posterior distribution (FPD) are on average the largest, but not much larger than predictive intervals derived using the true parameters. In the tests we conducted, they seemed on average to be larger by about one fifth at worst. Predictive intervals derived from the MLE and MAP estimators are on average shorter than those derived from the true parameters. This can be interpreted as an under-estimation of the uncertainty of the prediction when fixing the vector of correlation lengths to the most likely value given the observations, and this can explain the low observed coverage in Table 2.

In Figure 2, we give violin plots of coverage and mean length of Prediction Intervals in the two most extreme cases: correlation lengths 0.4 – 0.8 – 0.2 (very low correlation) and 0.8 – 1.0 – 0.9 (very high correlation). The results are similar and illustrate the fact that the FPD gives larger intervals in order to reach the derived coverage value.

6.3. A higher-dimensional case

In this subsection, we emulate using Simple Kriging the 10-dimensional Ackley function:

[TABLE]

The goal in this section is to emulate the Ackley function on the unit hypercube $[0,1]^{10}$ using design sets with 100 observation points. Although the impact of the design set type is not the focus of this study, we present the results with a randomly chosen design according to the Uniform distribution on the domain $[0,1]^{10}$, a design obtained through Latin Hypercube Sampling (LHS), and a design obtained through LHS and subsequently optimized to maximize the minimum distance between two points. The Simple Kriging model uses the null function as mean function and the Matérn anisotropic geometric covariance kernel family with smoothness parameter 5/2. The Gibbs reference posterior distribution is accessed through a sample of 1000 points. The conditional densities are sampled using the Metropolis algorithm with normal instrumental density with standard deviation $0.4$ and a 100-step burn-in period.

To evaluate the performance of prediction intervals, we follow steps 3, 4, 5, 6 and 8 of the method presented in this section (step 7 is skipped as the “values of the Gaussian process” are naturally the values of the Ackley function) with $n_{0}=1000$. The results are presented in Tables 4 and 5.

As is shown in Table 4, prediction intervals derived using the Full Posterior Distribution perform better than those derived from the MAP, which themselves perform better than those derived from the MLE. This order of performance is the same regardless of the type of design set, although the optimized design set leads to much worse performances on average for prediction intervals than unoptimized designs. The latter fact is not surprising since space-filling designs ensure than no two points can be very close to each other, which makes it harder to determine the correlation lengths.

As expected, prediction intervals derived from the Full Posterior Distribution are on average longer than those derived from the MAP and a fortiori the MLE. Notice that prediction intervals are on average shorter with the optimized design set, which explains the poorer performances in terms of coverage.

7. Conclusion and Perspectives

We provided theoretical foundation to the claim that the stationary distribution of the Markov chain underlying PIGS with random scanning order is the optimal compromise between the potentially incompatible conditional distributions.

This theory is mainly derived from intuitive conceptions of what a compromise should be. In places where such conceptions were inconclusive, we relied on concrete examples to precisely determine what was acceptable and what was not in a compromise and used it to complete the definition. One strength of this theory is that it can be applied to continuous as well as discrete probability distributions, whereas previous studies focused on the discrete, or even finite, case.

A question that remains open outside the finite-state case is how compatibility of conditional distributions is to be checked in practice. Although compromises are useful, not needing them is better.

Further investigation is needed to fully understand the properties of the optimal compromise. Nevertheless, its invariance by reparametrization and its respect of pairwise independence show that it preserves important features of the conditional distributions.

The theory of optimal compromise suggests a framework for deriving a new objective posterior distribution based on the conditionals yielded by the reference prior theory on Simple Kriging parameters. Applying this framework to Matérn anisotropic kernels, we showed prediction to have good frequentist properties.

Future work should investigate whether this posterior distribution formally corresponds to some joint prior distribution. And if it does, how the joint prior distribution could be accessed, and how it relates to the conditional prior distributions.

Regarding the specific Kriging application presented in this paper, the next step is to extend the framework to Universal Kriging, where instead of being known, the mean function is only assumed to be a linear combination of known functions $f_{1},...,f_{p}$. The linear coefficients $\beta_{1},...,\beta_{p}$ are then considered parameters of the model. This extension is of practical relevance, because the mean function can rarely be considered known. It can probably be done in the same way Berger et al. [6] extended the reference prior from the Simple Kriging to the Universal Kriging framework: they used the flat improper prior as joint prior on $\beta_{1},...,\beta_{p}$ conditional to $\sigma^{2}$ and $\boldsymbol{\theta}$ and used it to integrate $\beta_{1},...,\beta_{p}$ out of the likelihood function, and then proceeded to derive the reference prior on $\sigma^{2}$ and $\boldsymbol{\theta}$ with respect to the integrated likelihood.

A further extension would involve deriving an objective prior on the smoothness parameter $\nu$. In this endeavor, one should take into account the relation between correlation length $\boldsymbol{\theta}$ and smoothness $\nu$. Unfortunately, asymptotic theory is not of much help in this regard, as Anderes [2] shows that provided the spatial domain $\mathcal{D}$ is of dimension at least 5, then all parameters of the Matérn anisotropic geometric kernel are microergodic (Zhang [35] shows this to be untrue for spatial domains of dimension 1, 2 or 3, but the non-microergodic parameters are $\sigma^{2}$ and $\boldsymbol{\theta}$, not $\nu$). This means that the Gaussian measures on $\mathcal{D}$ corresponding to Gaussian Processes with two different smoothness parameters are orthogonal, which suggests that there exists a consistent estimator (the MLE possibly). Stein [32] (section 6.6) considers the Fisher information on $\boldsymbol{\theta}$ and $\nu$, and gives examples (with a one-dimensional sample space $\mathcal{D}$) showing that the Fisher information on these parameters depends a lot on the design set. Fisher information relative to the smoothness parameter $\nu$ increases when design points are chosen to be close to one another (relative to the "true" correlation length $\boldsymbol{\theta}$), whereas Fisher information relative to correlation length $\boldsymbol{\theta}$ is maximized for design points that are farther apart. This, according to him, is coherent with the fact that $\boldsymbol{\theta}$ has greater influence on the low frequency behavior of the Matérn kernel while $\nu$ has greater influence on its high frequency behavior. This also suggests to us that the smoothness parameter $\nu$, like the variance parameter $\sigma^{2}$, can only be meaningfully estimated if the vector of correlation lengths $\boldsymbol{\theta}$ is known. Otherwise, the estimator could hardly tell which design points are close to each other, which intuitively seems a prerequisite to evaluating the smoothness of the process. If we wish to apply the reference prior algorithm to the case where $\nu$ is unknown, we should thus probably derive the reference prior on $\nu$ conditional to $\boldsymbol{\theta}$.

Acknowledgements

The author would like to thank his PhD advisor Professor Josselin Garnier (École Polytechnique, Centre de Mathématiques Appliquées) for his guidance, Loic Le Gratiet (EDF R&D, Chatou) and Anne Dutfoy (EDF R&D, Saclay) for their advice and helpful suggestions. I also thank the Associate Editor and both anonymous reviewers for their constructive criticism which considerably improved this paper.

The author acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-13-MONU-0005 (project CHORUS).

?appendixname? A Proofs of Section 2

Proof of Proposition 2.6.

Let $P^{(0)}$ and $P^{(1)}$ be two compromises between $(\pi_{i})_{i\in[\![1,r]\!]}$ and let $t\in(0,1)$. Let us check that $(1-t)P^{(0)}+tP^{(1)}$ verifies Conditions (1) and (2) from Definition 2.4.

For every $i\in[\![1,r]\!]$, $\pi_{i}\left((1-t)P^{(0)}+tP^{(1)}\right)_{-i}=(1-t)\pi_{i}P_{-i}^{(0)}+t\pi_{i}P_{-i}^{(1)}$.

Let $A$ be a measurable set such that $\left((1-t)P^{(0)}+tP^{(1)}\right)(A)=0$. Then $P^{(0)}(A)=0$ and $P^{(1)}(A)=0$. Because both $P^{(0)}$ and $P^{(1))}$ are compromises, for every $i\in[\![1,r]\!]$, $\pi_{i}P_{-i}^{(0)}(A)=0$ and $\pi_{i}P_{-i}^{(1)}(A)=0$. So $\left(\pi_{i}\left((1-t)P^{(0)}+tP^{(1)}\right)_{-i}\right)(A)=0$ and Condition (1) is verified.

Now, for every $i\in[\![1,r]\!]$,

[TABLE]

So Condition (2) is also verified. ∎

Proof of Proposition 2.8.

Let $P_{G}$ be a Gibbs compromise between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$. Then, for any measurable set $A$ such that $P_{G}(A)=0$,

[TABLE]

So for every integer $i\in[\![1,r]\!]$, $\pi_{i}(P_{G})_{-i}(A)=0$. This fulfills the first condition of Definition 2.4.

Its second condition is also fulfilled because for every integer $i\in[\![1,r]\!]$,

[TABLE]

∎

Proof of Proposition 2.9.

Define the probability distribution $P$ on $\mathcal{A}$ as follows:

[TABLE]

Then for every $i\in[\![1,r]\!]$ the $i$-th $(r-1)$-marginal distribution $P_{-i}$ of $P$ is given by

[TABLE]

where the last equality is due to $(m_{\neq i})_{i\in[\![1,r]\!]}$ being compatible with $(\pi_{i})_{i\in[\![1,r]\!]}$. Plugging this into Equation (A.4), we obtain that $P$ is a Gibbs compromise. Equation (A.5) then yields the result.

∎

Proof of Theorem 2.10.

If $P_{G}$ is the unique Gibbs compromise between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$, then Proposition 2.9 asserts that $\left((P_{G})_{-i}\right)_{i\in[\![1,r]\!]}$ is the only compatible sequence of $(r-1)$-dimensional distributions. So any compromise has the same sequence of $(r-1)$-marginal distributions. Let $P$ be such a compromise.

For every $i\in[\![1,r]\!]$, $\pi_{i}(P_{G})_{-i}=\pi_{i}P_{-i}$ is absolutely continuous with respect to $P$ . So the average $P_{G}$ is also absolutely continuous with respect to $P$.

Let $\lambda$ be a positive measure on $\mathcal{A}$ such that $P$ is absolutely continuous with respect to $\lambda$. Because $P_{G}$ and $P$ share the same sequence of $(r-1)$-marginal distributions, for $\lambda$-almost any $\omega\in\Omega$,

[TABLE]

Moreover, Equation (A.4) implies that for $\lambda$-almost any $\omega\in\Omega$, $\frac{dP_{G}}{d\lambda}(\omega)$ is the arithmetic average between the $\frac{d(\pi_{i}(P_{G})_{-i})}{d\lambda}(\omega)$ ($i\in[\![1,r]\!]$), so it minimizes the mean squared error. Together with Equation (A.6), this implies that for $\lambda$-almost any $\omega\in\Omega$,

[TABLE]

Consequently, $E_{\lambda}(P_{G})\leqslant E_{\lambda}(P)$. Moreover, if $P\neq P_{G}$, then there exists $S\in\mathcal{A}$ such that $\lambda(S)>0$ and for every $\omega\in S$,

[TABLE]

So for every $\omega\in S$, Equation (A.7) is a strict inequality and thus $E_{\lambda}(P_{G})<E_{\lambda}(P)$. $P_{G}$ is therefore the unique optimal compromise with respect to $\lambda$. ∎

Proof of Proposition 2.11.

For every $i\in[\![1,r]\!]$, for every $\tilde{S}_{i}\in\tilde{\mathcal{A}}_{i}$,

[TABLE]

Now, for every $i\in[\![1,r]\!]$, for every $\tilde{T}_{i}\in\tilde{\mathcal{A}}_{i}$,

[TABLE]

So $(P_{G})_{-i}\ast f_{-i}=(P_{G}\ast f)_{-i}=(\tilde{P}_{G})_{-i}$. Then, returning to Equation (A.9),

[TABLE]

Finally, we obtain

[TABLE]

$\tilde{P}_{G}$ is therefore a Gibbs compromise between the sequence of Markov kernels $(\tilde{\pi}_{i})_{i\in[\![1,r]\!]}$. We now prove its uniqueness. For every $i\in[\![1,r]\!]$, $f_{i}$ is bijective and $f_{i}^{-1}$ is measurable. So for any Gibbs compromise $\tilde{Q}_{G}$ between the sequence of Markov kernels $(\tilde{\pi}_{i})_{i\in[\![1,r]\!]}$, $\tilde{Q}_{G}\ast f^{-1}$ is a Gibbs compromise between the sequence of Markov kernels $(\pi_{i})_{i\in[\![1,r]\!]}$. Given $P_{G}$ is the unique Gibbs compromise between the Markov kernels in this sequence, $\tilde{Q}_{G}\ast f^{-1}=P_{G}$, so $\tilde{Q}_{G}=\tilde{Q}_{G}\ast f^{-1}\ast f=P_{G}\ast f=\tilde{P}_{G}$.

∎

?appendixname? B Proofs of Section 4

The following holds where there is no mention of the contrary. When applied to a vector, $\|\cdot\|$ denotes the Euclidean norm and when applied to a matrix, it denotes the Frobenius norm. The choice of norm does not matter much because in finite-dimensional vector spaces, all norms are equivalent.

B.1. Differentiating the Matérn correlation kernel

Lemma B.1**.**

The partial derivative with respect to $\mu_{i}$ of the Matérn tensorized kernel of variance $\sigma^{2}$, smoothness $\nu$ and inverse correlation length vector $\boldsymbol{\mu}$ is:

[TABLE]

This can be rewritten as:

[TABLE]

?proofname?

.

The first assertion is a simple matter of differentiating Equation (4.9). In the following calculation, the fourth line is given by formula 9.6.28 (page 376) in Abramowitz and Stegun [1].

[TABLE]

From there, Equation (B.1) follows immediately. Rewriting it in the form given in (B.2) only requires us to recall $\Gamma(\nu)=(\nu-1)\Gamma(\nu-1)$ (case $\nu>1$), $\Gamma(1)=1$ (case $\nu=1$) and $\mathcal{K}_{\nu-1}=\mathcal{K}_{1-\nu}$ (case $\nu<1$). ∎

Lemma B.2**.**

The partial derivative with respect to $\mu_{i}$ of the Matérn geometric anisotropic kernel of variance $\sigma^{2}$, smoothness $\nu$ and inverse correlation length vector $\boldsymbol{\mu}$ is:

[TABLE]

This can be rewritten as:

[TABLE]

?proofname?.

The first assertion is a simple matter of differentiating Equation (4.7). In the following calculation, the fourth line is given by formula 9.6.28 (page 376) in Abramowitz and Stegun [1].

[TABLE]

From there, Equation (B.4) follows immediately. Rewriting it in the form given in (B.5) only requires us to recall $\Gamma(\nu)=(\nu-1)\Gamma(\nu-1)$ (case $\nu>1$), $\Gamma(1)=1$ (case $\nu=1$) and $\mathcal{K}_{\nu-1}=\mathcal{K}_{1-\nu}$ (case $\nu<1$). ∎

B.2. Accounting for low correlation: $\|\boldsymbol{\mu}\|\to\infty$

In this subsection, we consider a fixed design set of $n$ coordinate-distinct points $\boldsymbol{x}^{(k)}$ ($k\in[\![1,n]\!]$) in $\mathbb{R}^{r}$.

Lemma B.3**.**

For any Matérn anisotropic geometric or tensorized correlation kernel with smoothness $\nu>0$, for all $b<2\min(1,\nu)-1$ and $c>1$ (and if $\nu\neq 1$, for all $b\leqslant 2\min(1,\nu)-1$),

1. (a)

$\forall\boldsymbol{\mu}\in(\mathbb{R}_{+})^{r}$, $\|\frac{\partial}{\partial\mu_{i}}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}\|\leqslant M_{i,1}\;\mu_{i}^{-c}$. 2. (b)

$\forall\boldsymbol{\mu}\in(\mathbb{R}_{+})^{r}$, $\|\frac{\partial}{\partial\mu_{i}}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}\|\leqslant M_{i,2}\;\mu_{i}^{b}$.

?proofname?.

This can be gathered from Lemma B.1 or B.2 after recalling that 1) a Matérn kernel is a bounded function, 2) $\forall\nu\geqslant 0$, as $z\to+\infty$, $\mathcal{K}_{\nu}(z)\sim\sqrt{\pi}\exp(-z)/\sqrt{2z}$ ([1] 9.7.2) and 3) as $z\to 0$, $\mathcal{K}_{0}(z)\sim-\log(z)$ ([1] 9.6.8). ∎

Let us define

[TABLE]

Proposition B.4**.**

For any Matérn anisotropic geometric or tensorized correlation kernel with smoothness $\nu>0$, for all $\mu_{i}\in(0,+\infty)$, $\pi(\mu_{i}|\boldsymbol{\mu}_{-i})$, seen as a function of $\boldsymbol{\mu}$, is well defined and continuous over $\{\boldsymbol{\mu}\in[0,+\infty)^{r}:\mu_{i}\neq 0,\,\boldsymbol{\mu}_{-i}\neq\boldsymbol{0}_{r-1}\}$.

?proofname?.

For any given $\tilde{\boldsymbol{\mu}}\in[0,+\infty)^{r}$ such that $\tilde{\mu}_{i}\neq 0$ and $\tilde{\boldsymbol{\mu}}_{-i}\neq\boldsymbol{0}_{r-1}$, we prove that $\pi(\mu_{i}|\boldsymbol{\mu}_{-i})$, seen as a function of $\boldsymbol{\mu}$, is well defined and continuous at $\boldsymbol{\mu}=\tilde{\boldsymbol{\mu}}$.

For a start, notice that if $\boldsymbol{\mu}$ is confined to a sufficiently small neighborhood of $\tilde{\boldsymbol{\mu}}$, then $\|\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\|$ remains bounded. Therefore, Lemma B.3 implies that $\int_{0}^{\infty}f_{i}(\mu_{i}=t\;|\;\boldsymbol{\mu}_{-i})dt$ is finite and, thanks to the dominated convergence theorem, that it is continuous at $\boldsymbol{\mu}_{-i}=\tilde{\boldsymbol{\mu}}_{-i}$.

∎

Definition B.5**.**

An anisotropic geometric or tensorized correlation kernel is said to be “well-behaved” if its one-dimensional version is, for any set of parameters, a positive decreasing function on $[0,+\infty)$ that vanishes in the neighborhood of $+\infty$.

Lemma B.6**.**

Provided a coordinate-distinct design set is used, a well-behaved anisotropic geometric or tensorized correlation kernel parametrized by $\boldsymbol{\mu}$ has the following properties:

1. (a)

for any fixed $\boldsymbol{\mu}_{-i}\in\left[0,+\infty\right)^{r-1}$, it is a decreasing function of $\mu_{i}$ ; 2. (b)

as $\|\boldsymbol{\mu}\|\rightarrow\infty$, $\|\boldsymbol{\Sigma}_{\boldsymbol{\mu}}-\boldsymbol{I}_{n}\|\rightarrow 0$.

Lemma B.7**.**

For any well-behaved correlation kernel, as $\|\boldsymbol{\mu}\|\rightarrow\infty$, $\operatorname{Tr}\left[\frac{\partial}{\partial\mu_{i}}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\right]=o\left(\left\|\frac{\partial}{\partial\mu_{i}}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}\right\|\right)$.

?proofname?.

This result is due to the fact that all $\frac{\partial}{\partial\mu_{i}}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}$’s diagonal coefficients are null and $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}$ goes to the identity matrix as $\|\boldsymbol{\mu}\|\rightarrow\infty$. ∎

Let us now define

[TABLE]

Lemma B.8**.**

For any well-behaved correlation kernel, as $\|\boldsymbol{\mu}\|\rightarrow\infty$, $f_{i}(\mu_{i}\;|\;\boldsymbol{\mu}_{-i})\sim h_{i}(\mu_{i}\;|\;\boldsymbol{\mu}_{-i})$.

?proofname?.

Because $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}$ goes to the identity matrix, this is a direct consequence of Lemma B.7. ∎

Corollary B.9**.**

For any well-behaved correlation kernel, there exist $S>0$, and $0<a<b$ such that, whenever $\|\boldsymbol{\mu}\|\geqslant S$,

[TABLE]

In the following, $\boldsymbol{\Sigma}_{\boldsymbol{\mu}_{-i}}$ is the correlation matrix that would be obtained if $\mu_{i}$ were replaced by 0. Moreover, if $\boldsymbol{M}$ is a matrix, $\boldsymbol{M}^{(kl)}$ is its element in the $k$-th row and $l$-th column.

Lemma B.10**.**

If a well-behaved correlation kernel is used, there exist real constants $S>0$ and $c>0$ such that, for all $\mu_{i}\in(0,+\infty)$ and whenever $\|\boldsymbol{\mu}_{-i}\|\geqslant S$,

[TABLE]

?proofname?.

If a well-behaved correlation kernel is used, then for any for any $\epsilon>0$, Corollary B.9 implies that

[TABLE]

The last inequality holds because the Frobenius norm of any matrix is smaller than or equal to the sum of the absolute values of its elements and the correlation kernel is a decreasing function of $\mu_{i}$. Now, for all $k\neq l$, when $\mu_{i}\rightarrow+\infty$, $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{(kl)}\rightarrow 0$ and when $\mu_{i}=0$, $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{(kl)}=\boldsymbol{\Sigma}_{\boldsymbol{\mu}_{-i}}^{(kl)}$. From this, we gather that

[TABLE]

From this, we deduce that

[TABLE]

∎

This Lemma has the following immediate consequence:

Proposition B.11**.**

If a well-behaved tensorized kernel is used, there exists $S>0$ and for every $i\in[\![1,r]\!]$, there exists a function $M_{i}:(0,+\infty)\rightarrow(0,+\infty)$ such that for all $\|\boldsymbol{\mu}_{-i}\|\geqslant S$, $\pi_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})\geqslant M_{i}(\mu_{i})$.

?proofname?.

If a tensorized correlation kernel is used, for every pair of integers $(k,l)\in[\![1,r]\!]^{2}$ such that $k\neq l$, define the function $M_{i}^{(kl)}:(0,+\infty)\rightarrow(0,+\infty)$ ; $t\mapsto\left|\frac{d}{dt}\boldsymbol{\Sigma}_{\mu_{i}=t,\boldsymbol{\mu}_{-i}=\boldsymbol{0}_{r-1}}^{(kl)}\right|$.

[TABLE]

This fact, joined with Lemma B.10, yields the result. ∎

Proposition B.12**.**

Assume a well-behaved anisotropic geometric correlation kernel is used. If the corresponding one-dimensional kernel $K$ has the properties (P1) and (P2), then for every $i\in[\![1,r]\!]$, there exist positive functions $s_{i}$ and $m_{i}$ defined on $(0,+\infty)$ such that, for all $\|\boldsymbol{\mu}_{-i}\|\geqslant s_{i}(\mu_{i})$, $\pi_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})\geqslant m_{i}(\mu_{i})$.

(P1) : There exist $S_{1}>0$ and $M_{1}>0$ such that, for all $t\geqslant S_{1}$, $|K^{\prime}(t)|\geqslant M_{1}tK(t)$.

(P2) : For any $a>0$, there exist $S_{2}(a)>0$ and $M_{2}(a)>0$ such that, whenever $t\geqslant S_{2}(a)$, $K(t+a)\geqslant M_{2}(a)K(t)$.

?proofname?.

From (P1), we gather that for all $a>0$ and $t\geqslant S_{1}$, $|K^{\prime}(\sqrt{t^{2}+a^{2}})|\geqslant M_{1}\sqrt{t^{2}+a^{2}}K(\sqrt{t^{2}+a^{2}})$. Now, because the correlation kernel is well-behaved, $K$ is a decreasing function. As $\sqrt{t^{2}+a^{2}}\leqslant t+a$, $K(\sqrt{t^{2}+a^{2}})\geqslant K(t+a)$.

Plugging this into the previous inequality, we get $|K^{\prime}(\sqrt{t^{2}+a^{2}})|\geqslant M_{1}\sqrt{t^{2}+a^{2}}K(t+a)$.

If $t\geqslant\max(S_{1},S_{2}(a))$, we can then use (P2) to obtain

[TABLE]

Independently from this, we have the following algebraic fact:

[TABLE]

Because we use a well-behaved anisotropic geometric kernel, defining the function $M_{i}^{(kl)}:(0,+\infty)\rightarrow(0,+\infty)$ ; $t\mapsto\left(x_{j}^{(k)}-x_{j}^{(l)}\right)^{2}$, we can write:

[TABLE]

Setting $a_{kl}:=|x_{i}^{(k)}-x_{i}^{(l)}|\mu_{i}$ and $t_{kl}:=\|(\boldsymbol{x}_{-i}^{(k)}-\boldsymbol{x}_{-i}^{(l)})\boldsymbol{\mu}_{-i}\|$ (and thus, naturally, $\sqrt{t_{kl}^{2}+a_{kl}^{2}}=\|(\boldsymbol{x}^{(k)}-\boldsymbol{x}^{(l)})\boldsymbol{\mu}\|$), and provided $\|\boldsymbol{\mu}_{-i}\|$ is sufficiently large to make all $t_{kl}$s meet the conditions necessary to apply (P1) and (P2) (that depend in the case of (P2) on the $a_{kl}$s), Equation (B.16) yields the existence of some number $m_{i}^{(kl)}(\mu_{i})>0$ such that

[TABLE]

Finally, setting $m_{i}(\mu_{i}):=\mu_{i}\min_{k\neq l}\left[\left(x_{i}^{(k)}-x_{i}^{(l)}\right)^{2}m_{i}^{(kl)}(\mu_{i})\right]$, we get

[TABLE]

Then, applying Lemma B.10 yields the result. ∎

Proposition B.13**.**

Matérn one-dimensional kernels with smoothness parameter $\nu>1$ have the properties (P1) and (P2) of Proposition B.12.

?proofname?.

(P1) is given by Lemma B.2, after noticing that, denoting by $K_{\nu}$ the Matérn one-dimensional kernel of smoothness $\nu>1$, provided $t$ is sufficiently large, $K_{\nu}(t)\leqslant K_{\nu-1}(t)$.

This inequality ensues from the fact that $\forall\nu\geqslant 0$, as $t\to+\infty$, $\mathcal{K}_{\nu}(t)\sim\sqrt{\pi}\exp(-t)/\sqrt{2t}$ ([1] 9.7.2) and thus $K_{\nu}(t)\sim 2/\Gamma(\nu)(\sqrt{\nu}t)^{\nu}\sqrt{\pi/(4\sqrt{\nu}t)}\exp(-2\sqrt{\nu}t)$. Moreover, this last equivalence relation also implies (P2). ∎

Proposition B.14**.**

For Matérn anisotropic geometric kernels with smoothness $\nu>1$ and Matérn tensorized correlation kernels with smoothness $\nu>0$,

for any $\delta>0$, $i\in[\![1,r]\!]$ and $\mu_{i}\in(0,+\infty)$, there exists $b_{i,\delta}(\mu_{i})>0$ such that, if $\|\boldsymbol{\mu}_{-i}\|\geqslant\delta$, then $\pi_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})\geqslant b_{i,\delta}(\mu_{i})$.

?proofname?.

Matérn correlation kernels with such smoothness parameters make Proposition B.11 or B.12 applicable. Therefore, there exist $s_{i}(\mu_{i})>0$ and $m_{i}(\mu_{i})>0$ such that, if $\|\boldsymbol{\mu}_{-i}\|\geqslant s_{i}(\mu_{i})$, $\pi_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})\geqslant m_{i}(\mu_{i})$. Besides, we know from Proposition B.4 that $\pi_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})$, seen as a function of $\boldsymbol{\mu}_{-i}$, is continuous and positive over the compact set $\{\boldsymbol{\mu}_{-i}:\delta\leqslant\|\boldsymbol{\mu}_{-i}\|\leqslant s_{i}(\mu_{i})\}$. Thus its minimum $\tilde{m}_{i,\delta}(\mu_{i})$ on this set is positive and we obtain the result by setting $b_{i,\delta}(\mu_{i}):=\min(m_{i}(\mu_{i}),\tilde{m}_{i,\delta}(\mu_{i}))$. ∎

Proposition B.15**.**

For Matérn anisotropic geometric correlation kernels with smoothness $\nu>1$ and for Matérn tensorized correlation kernels with smoothness $\nu>0$, for any $\boldsymbol{y}\in\mathbb{R}^{n}\setminus\{0\}^{n}$, any $\delta>0$ and any $\mu_{i}\in(0,+\infty)$, there exists $b_{i,\delta,\boldsymbol{y}}(\mu_{i})>0$ such that, if $\|\boldsymbol{\mu}_{-i}\|\geqslant\delta$, then $\pi_{i}(\mu_{i}|\boldsymbol{y},\boldsymbol{\mu}_{-i})\geqslant b_{i,\delta,\boldsymbol{y}}(\mu_{i})$.

?proofname?.

Set $\delta>0$ and $\boldsymbol{y}\in\mathbb{R}^{n}\setminus\{0\}^{n}$. There exist $m_{\delta}>0$ and $M_{\delta}>0$ s.t. $\forall\boldsymbol{\mu}\in(0,+\infty)^{r}$, $\|\boldsymbol{\mu}_{-i}\|\geqslant\delta\Rightarrow m_{\delta}\leqslant\|\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\|\leqslant M_{\delta}$, so there also exist $m_{\delta,\boldsymbol{y}}>0$ and $M_{\delta,\boldsymbol{y}}>0$ s.t. $m_{\delta,\boldsymbol{y}}\leqslant L(\boldsymbol{y}|\boldsymbol{\mu})\leqslant M_{\delta,\boldsymbol{y}}$. This, combined with Proposition B.14, yields the result. ∎

B.3. Accounting for high correlation: $\|\boldsymbol{\mu}\|\to 0$

This part of the proof relies on the combination of some spectral study of the Matérn kernels and on the study of the matrices that are part of the series expansion of the correlation matrix $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}$ when $\|\boldsymbol{\mu}\|\to 0$ for three types of Matérn kernels: isotropic, tensorized and anisotropic geometric.

Lemma B.16**.**

There exists a covariance kernel $\widetilde{K}_{r,\nu}$ such that for any design set $\boldsymbol{x}^{(1)},...,\boldsymbol{x}^{(n)}$, for all $\boldsymbol{\mu}\in(\mathbb{R}_{+})^{r}$ and $\boldsymbol{\xi}=(\xi_{1},...,\xi_{n})\in\mathbb{R}^{n}$,

[TABLE]

where $f_{r,\nu}(t)=(2\sqrt{\nu})^{-r-2\nu}t^{-r}$ if $t\geqslant\left(2\sqrt{\nu}\right)^{-1}$ and $f_{r,\nu}(t)=t^{2\nu}$ if $t\leqslant\left(2\sqrt{\nu}\right)^{-1}$.

?proofname?.

For all $\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{r}$, $K_{r,\nu}(\boldsymbol{x}-\boldsymbol{y})=\int_{\mathbb{R}^{r}}\widehat{K}_{r,\nu}(\boldsymbol{\omega})e^{i\langle\boldsymbol{\omega}|\boldsymbol{x}-\boldsymbol{y}\rangle}d\boldsymbol{\omega}$.

[TABLE]

Now, let $\widetilde{K}_{r,\nu}$ be the function with Fourier transform $\widehat{\widetilde{K}}_{r,\nu}(\boldsymbol{\omega})=\mathbf{1}_{\{\|\boldsymbol{\omega}\|\geqslant 1\}}\|\boldsymbol{\omega}\|^{-r-2\nu}$. According to Bochner’s theorem, $\widetilde{K}_{r,\nu}$ is a correlation kernel, which leads to the conclusion. ∎

Lemma B.17**.**

For every design set with coordinate-distinct points $\boldsymbol{x}^{(1)},...,\boldsymbol{x}^{(n)}$, there exists a constant $c_{\boldsymbol{x}}>0$ such that for all $\boldsymbol{\mu}\in(\mathbb{R}_{+})^{r}$,

[TABLE]

where $f_{r,\nu}(t)=(2\sqrt{\nu})^{-r-2\nu}t^{-r}$ if $t\geqslant\left(2\sqrt{\nu}\right)^{-1}$ and $f_{r,\nu}(t)=t^{2\nu}$ if $t\leqslant\left(2\sqrt{\nu}\right)^{-1}$.

?proofname?.

For every design set $\boldsymbol{x}^{(1)},...,\boldsymbol{x}^{(n)}$, the set of all design sets that can be written $\frac{\boldsymbol{\mu}}{\|\boldsymbol{\mu}\|_{\infty}}\boldsymbol{x}^{(1)},...,\frac{\boldsymbol{\mu}}{\|\boldsymbol{\mu}\|_{\infty}}\boldsymbol{x}^{(n)}$ ($\boldsymbol{\mu}\in(\mathbb{R}_{+})^{r}$) is compact. If the design set $\boldsymbol{x}^{(1)},...,\boldsymbol{x}^{(n)}$ has coordinate-distinct points, then every design set in the aforementioned compact set has no overlapping points. Thus the conclusion follows from Lemma B.16. ∎

Proposition B.18**.**

With Matérn anisotropic geometric or tensorized kernels, for every design set with coordinate-distinct points $\boldsymbol{x}^{(1)},...,\boldsymbol{x}^{(n)}$, as $\|\boldsymbol{\mu}\|\to 0$, $\left\|\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\right\|=O(\|\boldsymbol{\mu}\|^{-2\nu})$.

?proofname?.

For Matérn anisotropic geometric kernels, we need only apply Lemma B.17. In the case of tensorized Matérn kernels, analoguous results to Lemma B.16 and then Lemma B.17 may be used. ∎

Abramowitz and Stegun [1] give the following results on the modified Bessel function of second kind (usually noted $K_{\nu}$ and which we note $\mathcal{K}_{\nu}$ in order to avoid confusion with the Matérn correlation kernel). If $I_{\nu}$ is the modified Bessel function of first kind and $\psi$ is the function defined in (6.3.2) by $\psi:\mathbb{N}\setminus\{0\}\rightarrow\mathbb{R}$ ; $k\mapsto-\gamma+\sum_{i=1}^{k-1}i^{-1}$ :

[TABLE]

This gives us the series expansion of $K_{1,\nu}(z)$ ($\nu\in[0,+\infty)\setminus\mathbb{N}$) when $z\rightarrow 0$:

[TABLE]

In the remainder of this subsection, we consider a fixed design set with $n$ coordinate-distinct points $\boldsymbol{x}^{(k)}$ ($k\in[\![1,n]\!]$) in $\mathbb{R}^{r}$. Moreover, all Matérn kernels we consider are assumed to have non-integer smoothness parameter $\nu$.

Let us now define, for every nonnegative integer $k<\nu$ the matrix $\boldsymbol{D}^{k}$ whose $(i,j)$ element is

[TABLE]

Let us also define the matrix $\boldsymbol{D}^{\nu}$ whose $(i,j)$ element is

[TABLE]

If the correlation kernel is Matérn isotropic, $\boldsymbol{\Sigma}_{\mu}$ has the following series expansion if $\nu$ is not an integer when $\mu\rightarrow 0+$:

[TABLE]

In this expansion, $\mu^{-2\nu}\|\boldsymbol{R}_{\mu}\|\rightarrow 0$.

For any integer $i\in[\![1,r]\!]$ and any nonnegative integer $k<\nu$ define the matrix $\boldsymbol{D}_{i}^{k}$ whose $(m,p)$ element is

[TABLE]

and also the matrix $\boldsymbol{D}_{i}^{\nu}$ whose $(m,p)$ element is

[TABLE]

For every $i\in[\![1,r]\!]$, if the points in the design set differed only through their $i$-th coordinate, the series expansion of the correlation matrix (using a Matérn anisotropic geometric or tensorized kernel) when $\|\boldsymbol{\mu}\|\to 0$ (and thus when $\mu_{i}\to 0$) would be

[TABLE]

where $\mu_{i}^{-2\nu}$ $\|\boldsymbol{R}_{\mu_{i}}\|\rightarrow 0$ as $\mu_{i}\to 0$.

Note the following identities:

[TABLE]

If a tensorized correlation kernel is used, the correlation matrix $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{tens}$ may be written

[TABLE]

where the subscript $\circ$ above the symbol $\prod$ serves to denote the Hadamard product of matrices.

In case a Matérn anisotropic geometric kernel is used, then define for any nonnegative interger $k<\nu$ the matrix $\boldsymbol{D}^{k}(\boldsymbol{\mu})$ whose $(m,p)$ element is

[TABLE]

where $d_{m,p}(\boldsymbol{\mu})=\left\|\left(\boldsymbol{x}^{(m)}-\boldsymbol{x}^{(p)}\right)\boldsymbol{\mu}\right\|$.

And, similarly, we may define the matrix $\boldsymbol{D}^{\nu}(\boldsymbol{\mu})$ whose $(m,p)$ element is

[TABLE]

We thus have (if $\nu\in[0,+\infty)\setminus\mathbb{N}$)

[TABLE]

where $\|\boldsymbol{\mu}\|^{-2\nu}$ $\|\boldsymbol{R}_{\boldsymbol{\mu}}^{geom}\|\rightarrow 0$ as $\|\boldsymbol{\mu}\|\to 0$.

Similar identities to those of Equation (B.31) can be derived to make Equation (B.37) more explicit for small values of $\nu$.

[TABLE]

Fortunately, for small values of $\nu$, $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{tens}$ can also be simply written.

[TABLE]

In the three expressions above, $\|\boldsymbol{\mu}\|^{-2\nu}$ $\|\boldsymbol{R}_{\boldsymbol{\mu}}^{tens}\|\rightarrow 0$ as $\|\boldsymbol{\mu}\|\to 0$.

Define $k_{\nu}$ as the orthogonal complement in $\mathbb{R}^{n}$ of the vector space spanned by:

1. (1)

if $\nu\in(0,1)$: $\boldsymbol{1}$; 2. (2)

if $\nu\in(1,2)$: $\boldsymbol{1}$ and $\boldsymbol{X}_{i}$ ($i\in[\![1,r]\!]$); 3. (3)

if $\nu\in(2,3)$: $\boldsymbol{1}$ and $\boldsymbol{X}_{i}$ ($i\in[\![1,r]\!]$) and $\boldsymbol{X}_{i}\circ\boldsymbol{X}_{j}$ ($i,j\in[\![1,r]\!]$).

Clearly, for any $\nu\in(0,1)\cup(1,2)\cup(2,3)$, for any vector $\boldsymbol{v}\in k_{\nu}$,

[TABLE]

Since when $\mu\to 0$ $\|\boldsymbol{D}^{\nu}(\boldsymbol{\mu})\|=O(\|\boldsymbol{\mu}\|^{2\nu})$, $\|\boldsymbol{R}_{\boldsymbol{\mu}}^{geom}\|=o(\|\boldsymbol{\mu}\|^{2\nu})$ and $\|\boldsymbol{R}_{\boldsymbol{\mu}}^{tens}\|=o(\|\boldsymbol{\mu}\|^{2\nu})$, for any $\boldsymbol{\mu}\in(0,+\infty)^{r}$ such that $\|\boldsymbol{\mu}\|$ is small enough, there exists $c>0$ such that for any $\boldsymbol{v}\in k_{\nu}$,

[TABLE]

Proposition B.19**.**

For a Matérn anisotropic geometric or tensorized correlation kernel with smoothness parameter $\nu\in(0,1)\cup(1,2)\cup(2,3)$, for any vector $\boldsymbol{y}\in\mathbb{R}^{n}$ not orthogonal to $k_{\nu}$, when $\|\boldsymbol{\mu}\|\to 0$, $\|\boldsymbol{\mu}\|^{-2\nu}=O\left(\boldsymbol{y}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\boldsymbol{y}\right)$.

?proofname?.

Let $d_{\nu}$ be the dimension of $k_{\nu}$ and let $\boldsymbol{O}_{\nu}$ be an orthogonal $n\times n$ matrix whose first $n-d_{\nu}$ columns form an orthonormal basis of $k_{\nu}^{\perp}$ and whose last $d_{\nu}$ columns form an orthonormal basis of $k_{\nu}$.

Then $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}=\boldsymbol{O}_{\nu}\boldsymbol{O}_{\nu}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}\boldsymbol{O}_{\nu}\boldsymbol{O}_{\nu}^{\top}$. Consider the following decomposition of $\boldsymbol{O}_{\nu}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}\boldsymbol{O}_{\nu}$:

[TABLE]

where the blocks $\boldsymbol{A}_{\boldsymbol{\mu}}$, $\boldsymbol{B}_{\boldsymbol{\mu}}$ and $\boldsymbol{C}_{\boldsymbol{\mu}}$ are respectively $(n-d_{\nu})\times(n-d_{\mu})$, $(n-d_{\mu})\times d_{\mu}$ and $d_{\mu}\times d_{\mu}$ matrices. Note that $\boldsymbol{A}_{\boldsymbol{\mu}}$ and $\boldsymbol{C}_{\boldsymbol{\mu}}$ represent the restriction of the scalar product defined by $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}$ to $k_{\nu}^{\perp}$ and $k_{\nu}$ respectively. When $\|\boldsymbol{\mu}\|$ is small enough, defining $c>0$ as in Equation (B.46), $\|\boldsymbol{C}_{\boldsymbol{\mu}}\|\leqslant c\|\boldsymbol{\mu}\|^{2\nu}$.

[TABLE]

For any vector $\boldsymbol{y}\in\mathbb{R}^{n}$, there exist $\boldsymbol{y}_{1}\in\mathbb{R}^{n-d_{\nu}}$ and $\boldsymbol{y}_{2}\in\mathbb{R}^{d_{\nu}}$ such that

[TABLE]

[TABLE]

Given $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}$ is positive definite, the diagonal block $\left(\boldsymbol{A}-\boldsymbol{B}_{\boldsymbol{\mu}}\boldsymbol{C}_{\boldsymbol{\mu}}^{-1}\boldsymbol{B}_{\boldsymbol{\mu}}^{\top}\right)^{-1}$ is positive definite too. This implies $\boldsymbol{y}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\boldsymbol{y}\geqslant\boldsymbol{y}_{2}^{\top}\boldsymbol{C}_{\boldsymbol{\mu}}^{-1}\boldsymbol{y}_{2}.$ When $\|\boldsymbol{\mu}\|$ is small enough, $\boldsymbol{y}_{2}^{\top}\boldsymbol{C}_{\boldsymbol{\mu}}^{-1}\boldsymbol{y}_{2}\geqslant c^{-1}\|\boldsymbol{\mu}\|^{-2\nu}\|\boldsymbol{y}_{2}\|^{2}$.

If $\boldsymbol{y}$ is not orthogonal to $k_{\nu}$, then $\|\boldsymbol{y}_{2}\|\neq 0$ and thus $\|\boldsymbol{\mu}\|^{-2\nu}=O(\boldsymbol{y}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\boldsymbol{y})$.

∎

Proposition B.20**.**

Assume $\nu\in(1,2)\cup(2,3)$. For every $\boldsymbol{y}\in\mathbb{R}^{n}$ that is not orthogonal to the vector subspace $k_{\nu}$, $L(\boldsymbol{y}|\boldsymbol{\mu})f_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})$ is a bounded function of $\boldsymbol{\mu}$.

?proofname?.

Let $v_{1}(\boldsymbol{\mu})\geqslant v_{2}(\boldsymbol{\mu})\geqslant...\geqslant v_{n}(\boldsymbol{\mu})$ be the ordered eigenvalues of $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}$. We can now rewrite $L(\boldsymbol{y}|\boldsymbol{\mu})$ as

[TABLE]

Proposition B.19 asserts that for any $\boldsymbol{y}\in\mathbb{R}^{n}$ that is not orthogonal to $k_{\nu}$, $\left(\boldsymbol{y}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\boldsymbol{y}\right)^{-1}=O(\|\boldsymbol{\mu}\|^{2\nu})$ for $\|\boldsymbol{\mu}\|\to 0$. Besides, Proposition B.18 asserts that $\left\|\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\right\|=O(\|\boldsymbol{\mu}\|^{-2\nu})$, so $\left(\boldsymbol{y}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\boldsymbol{y}\right)^{-1}=O\left(\left\|\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\right\|^{-1}\right)$. This implies that for every integer $i\in[\![1,r]\!]$, $v_{k}(\boldsymbol{\mu})^{-1}\left(\boldsymbol{y}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\boldsymbol{y}\right)^{-1}=O(1)$.

Clearly, $\lim_{\|\boldsymbol{\mu}\|\to 0}|\boldsymbol{1}^{\top}\boldsymbol{v}_{1}(\boldsymbol{\mu})|=\|\boldsymbol{1}\|$ and $\lim_{\|\boldsymbol{\mu}\|\to 0}v_{1}(\boldsymbol{\mu})=n$. The latter implies $\lim_{\|\boldsymbol{\mu}\|\to 0}v_{1}(\boldsymbol{\mu})^{-1}=n^{-1}$ and $v_{1}(\boldsymbol{\mu})^{-1}\left(\boldsymbol{y}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\boldsymbol{y}\right)^{-1}=O\left(\|\boldsymbol{\mu}\|^{2\nu}\right)$. Now, for every $\boldsymbol{\mu}\in(0,+\infty)^{r}$, we may (thanks to the axiom of choice) choose a unit eigenvector $\boldsymbol{v}_{1}(\boldsymbol{\mu})$ corresponding to the largest eigenvalue $v_{1}(\boldsymbol{\mu})$ and $\boldsymbol{v}_{2}(\boldsymbol{\mu})$ corresponding to the second largest eigenvalue $v_{2}(\boldsymbol{\mu})$. Because $\boldsymbol{\Sigma}_{\boldsymbol{\mu}}$ is symmetric, $\boldsymbol{v}_{1}(\boldsymbol{\mu})^{\top}\boldsymbol{v}_{2}(\boldsymbol{\mu})=0$ for all $\boldsymbol{\mu}\in(0,+\infty)^{r}$, so $\lim_{\|\boldsymbol{\mu}\|\to 0}\boldsymbol{1}^{\top}\boldsymbol{v}_{2}(\boldsymbol{\mu})=0$.

[TABLE]

For all $\boldsymbol{\mu}\in(0,+\infty)^{r}$, let $i(\boldsymbol{\mu})$ be the smallest integer $i\in[\![1,r]\!]$ such that $\mu_{i}=\max_{j=1}^{r}\mu_{j}$. Now for every integer $i\in[\![1,r]\!]$ let $\boldsymbol{w}_{i}(\boldsymbol{\mu})$ be the unit vector that belongs to the space spanned by $\boldsymbol{v}_{1}(\boldsymbol{\mu})$ and $\boldsymbol{X}_{i}$ that verifies $\boldsymbol{v}_{1}(\boldsymbol{\mu})^{\top}\boldsymbol{w}_{i}(\boldsymbol{\mu})=0$ and $\boldsymbol{X}_{i}^{\top}\boldsymbol{w}_{i}(\boldsymbol{\mu})>0$.

[TABLE]

Because $\lim_{\|\boldsymbol{\mu}\|\to 0}|\boldsymbol{1}^{\top}\boldsymbol{v}_{1}(\boldsymbol{\mu})|=\|\boldsymbol{1}\|$, $\liminf_{\|\boldsymbol{\mu}\|\to 0}\boldsymbol{X}_{i(\boldsymbol{\mu})}^{\top}\boldsymbol{w}_{i(\boldsymbol{\mu})}(\boldsymbol{\mu})\geqslant\min_{i=1}^{r}\lim_{\|\boldsymbol{\mu}\|\to 0}\boldsymbol{X}_{i}^{\top}\boldsymbol{w}_{i}(\boldsymbol{\mu})>0$, so there exists a constant $c_{2}>0$ such that when $\|\boldsymbol{\mu}\|$ is small enough

[TABLE]

Recall $v_{2}(\boldsymbol{\mu})=\max\{\boldsymbol{\xi}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}\boldsymbol{\xi}|\boldsymbol{\xi}\in S^{n-1}\;\mathrm{and}\;\boldsymbol{\xi}^{\top}\boldsymbol{v}_{1}(\boldsymbol{\mu})\}$, so a fortiori $v_{2}(\boldsymbol{\mu})\geqslant c_{2}\|\boldsymbol{\mu}\|^{2}$.

This implies $v_{2}(\boldsymbol{\mu})^{-1}=O(\|\boldsymbol{\mu}\|^{-2})$ and therefore $v_{2}(\boldsymbol{\mu})^{-1}\left(\boldsymbol{y}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\boldsymbol{y}\right)^{-1}=O(\|\boldsymbol{\mu}\|^{2(\nu-1)})$ .

Finally, $L(\boldsymbol{y}|\boldsymbol{\mu})=O\left(\|\boldsymbol{\mu}\|^{\nu}\right)O\left(\|\boldsymbol{\mu}\|^{\nu-1}\right)=O\left(\|\boldsymbol{\mu}\|^{2\nu-1}\right)$. Given that $f_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})=O(\|\boldsymbol{\mu}\|^{1-2\nu})$, the product $L(\boldsymbol{y}|\boldsymbol{\mu})f_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})$ is bounded when $\|\boldsymbol{\mu}\|\to 0$. ∎

Proposition B.21**.**

Assume $\nu\in(0,1)$. For every $\boldsymbol{y}\in\mathbb{R}^{n}$ that is not collinear to $\boldsymbol{1}$, when $\|\boldsymbol{\mu}\|\to 0$, $L(\boldsymbol{y}|\boldsymbol{\mu})f_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})=O(\mu_{i}^{-1+\nu})$.

?proofname?.

This proof is similar to the previous one, so we use the same notations. $v_{1}(\boldsymbol{\mu})^{-1}=O(1)$, so $v_{1}(\boldsymbol{\mu})^{-1}\left(\boldsymbol{y}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\boldsymbol{y}\right)^{-1}=O(\|\boldsymbol{\mu}\|^{2\nu})$, which yields that $L(\boldsymbol{y}|\boldsymbol{\mu})=O(\|\boldsymbol{\mu}\|^{\nu})$.

This implies that $L(\boldsymbol{y}|\boldsymbol{\mu})\left\|\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}\right\|=O(\|\boldsymbol{\mu}\|^{-\nu})=O(\mu_{i}^{-\nu})$.

By Lemma B.3, $\left\|\frac{\partial}{\partial\mu_{i}}\boldsymbol{\Sigma}_{\boldsymbol{\mu}}\right\|=O(\mu_{i}^{-1+2\nu})$. Putting all this together, $L(\boldsymbol{y}|\boldsymbol{\mu})f_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})=O(\mu_{i}^{-1+\nu})$. ∎

Proposition B.22**.**

For Matérn anisotropic geometric or tensorized kernels with smoothness parameter $\nu\in(0,1)\cup(1,2)\cup(2,3)$, if $\boldsymbol{y}\in\mathbb{R}^{n}$ is not orthogonal to $k_{\nu}$, then the conditional posterior distribution $\pi_{i}(\mu_{i}|\boldsymbol{y},\boldsymbol{\mu}_{-i})$, seen as a function of $\boldsymbol{\mu}$, is continuous over $\{\boldsymbol{\mu}\in[0,+\infty)^{r}:\mu_{i}\neq 0\}$.

Moreover,

[TABLE]

?proofname?.

Given Proposition B.4 and the fact that $\forall\boldsymbol{y}\in\mathbb{R}^{n}$, $\forall i\in[\![1,r]\!]$ and $\forall\mu_{i}\in(0,+\infty)$, as $\|\boldsymbol{\mu}_{-i}\|\to 0$, $L(\boldsymbol{y}|\boldsymbol{\mu})f_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})$ converges pointwise to $L(\boldsymbol{y}|\mu_{i},\boldsymbol{\mu}_{-1}=\boldsymbol{0}_{r-1})f_{i}(\mu_{i}|\boldsymbol{\mu}_{-i}=\boldsymbol{0}_{r-1})$, we only need to show that

[TABLE]

Lemma B.3 implies that there exists $M_{i}>0$ such that

[TABLE]

Lemma B.6 then ensures $\left\|\boldsymbol{\Sigma}_{\boldsymbol{\mu}}^{-1}-\boldsymbol{I}_{n}\right\|\to 0$ as $\mu_{i}\to+\infty$, so the right member of the inequality is integrable in the neighborhood of $+\infty$. Let us now focus on the neighborhood of 0.

If $\nu>1$, Proposition B.20 asserts that $L(\boldsymbol{y}|\boldsymbol{\mu})f_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})$ is bounded in the neighborhood of 0.

If $\nu<1$, Proposition B.21 asserts that $L(\boldsymbol{y}|\boldsymbol{\mu})f_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})\mu_{i}^{1-\nu}$ is bounded in the neighborhood of 0.

Therefore, there exists a function independent of $\boldsymbol{\mu}_{-i}$ that is both greater than $L(\boldsymbol{y}|\boldsymbol{\mu})f_{i}(\mu_{i}|\boldsymbol{\mu}_{-i})$ and integrable over $\mu_{i}\in(0,+\infty)$, so the dominated convergence theorem is applicable.

∎

B.4. Lower bound for conditional reference posterior densities

The following Lemma provides the key to proving Theorem 4.2:

Lemma B.23**.**

In a Simple Kriging model with the characteristics described above, there exists a hyperplane $\mathcal{H}$ of $\mathbb{R}^{n}$ such that, for any $\boldsymbol{y}\in\mathbb{R}^{n}\setminus\mathcal{H}$ and any $i\in[\![1,r]\!]$, there exists a measurable function $m_{i,\boldsymbol{y}}:(0,+\infty)\rightarrow(0,+\infty)$ such that, for all $\boldsymbol{\mu}_{-i}\in(0,+\infty)^{r-1}$, the conditional reference posterior density verifies:

[TABLE]

?proofname?.

This proof consists in combining Proposition B.15 and Proposition B.22, which respectively deal with large and small values of $\|\boldsymbol{\mu}_{-i}\|$.

Proposition B.15 implies that for any $\boldsymbol{y}\in\mathbb{R}^{n}\setminus\{0\}^{n}$, for any $i\in[\![1,r]\!]$ and any $\mu_{i}\in(0,+\infty)$, there exists a compact neighborhood $N_{i}(\mu_{i})$ of $\boldsymbol{0}_{r-1}$ within $[0,+\infty)^{r}$ such that

[TABLE]

The vector space $k_{\nu}\subset\mathbb{R}^{n}$ has dimension greater or equal to

1. (a)

$n-1$ if $\nu\in(0,1)$; 2. (b)

$n-(r+1)$ if $\nu\in(1,2)$; 3. (c)

$n-(r+1)(r+2)/2$ if $\nu\in(2,3)$.

For all Simple Kriging models tackled by this lemma, the dimension of $k_{\nu}$ is therefore greater or equal to 1. Its orthogonal complement $k_{\nu}^{\perp}$ is then included within a hyperplane $\mathcal{H}$ of $\mathbb{R}^{n}$. Assuming $\boldsymbol{y}\in\mathbb{R}^{n}\setminus\mathcal{H}$, Proposition B.22 ensures that $\pi(\mu_{i}|\boldsymbol{y},\boldsymbol{\mu}_{-i})$ is a continuous and positive function of $\boldsymbol{\mu}$ on $\{\boldsymbol{\mu}\in[0,+\infty)^{r}:\mu_{i}\neq 0\}$. In particular, this implies that for any $\mu_{i}\in(0,+\infty)$ and any compact neighborhood $N_{i}(\mu_{i})$ of $\boldsymbol{0}_{r-1}$ within $[0,+\infty)^{r}$,

[TABLE]

Putting this together, if $\boldsymbol{y}\in\mathbb{R}^{n}\setminus\mathcal{H}$, for any $i\in[\![1,r]\!]$ and any $\mu_{i}\in(0,+\infty)$,

[TABLE]

The mapping $m_{i,\boldsymbol{y}}:\mu_{i}\mapsto\inf\{\pi(\mu_{i}|\boldsymbol{y},\boldsymbol{\mu}_{-i}):\boldsymbol{\mu}_{-i}\in[0,+\infty)^{r-1}\}$, which is measurable on $(0,+\infty)$ is therefore also positive on $(0,+\infty)$.

∎

Proof of Theorem 4.2.

Lemma B.23 implies that $\forall\boldsymbol{\mu}^{(0)}\in(0,+\infty)^{r}$,

[TABLE]

and thus $\forall\boldsymbol{\mu}^{(0)}\in(0,+\infty)^{r}$, $\forall n\geqslant r$,

[TABLE]

Defining $f_{\boldsymbol{y}}(\boldsymbol{\mu}):=r^{-r}\prod_{i=1}^{r}m_{i,\boldsymbol{y}}(\mu_{i})$, $f_{\boldsymbol{y}}$ is a measurable positive function. Therefore $f_{\boldsymbol{y}}$ is the density with respect to the Lebesgue measure of a positive measure with mass $\epsilon_{\boldsymbol{y}}>0$. So $\epsilon_{\boldsymbol{y}}^{-1}f_{\boldsymbol{y}}$ is a probability density with respect to the Lebesgue measure and the Markov kernel $P_{\boldsymbol{y}}$ thus satisfies the uniform $(n,\epsilon_{\boldsymbol{y}})$ Doeblin condition:

[TABLE]

This implies that $P_{\boldsymbol{y}}$ is uniformly ergodic: it has a unique invariant probability distribution $\pi_{G}(\cdot|\boldsymbol{y})$ and $\lim_{n\to\infty}\sup_{\boldsymbol{\mu}^{(0)}\in(0,+\infty)^{r}}\|P_{\boldsymbol{y}}^{n}(\boldsymbol{\mu}^{(0)},\cdot)-\pi_{G}(\cdot|\boldsymbol{y})\|_{TV}=0$, where $\|\cdot\|_{TV}$ is the total variation norm. By definition, $\pi_{G}(\cdot|\boldsymbol{y})$ is the Gibbs compromise between the incompatible posterior conditionals $\pi_{i}(\mu_{i}|\boldsymbol{y},\boldsymbol{\mu}_{-i})$.

∎

?refname?

• Abramowitz and Stegun [1964]

M. Abramowitz and I. A. Stegun.

Handbook of mathematical functions with formulas, graphs, and mathematical tables, volume 55 of Applied Mathematics Series.

National Bureau of Standards, 1964.

• Anderes [2010]

E. Anderes.

On the consistent separation of scale and variance for Gaussian random fields.

Annals of Statistics, 38(2):870–893, 2010.

• Arnold et al. [2001]

B. C. Arnold, E. Castillo, and J. M. Sarabia.

Conditionally specified distributions: an introduction.

Statistical Science, 16(3):268–269, 2001.

• Berger [2006]

J. Berger.

The case for objective bayesian analysis.

Bayesian analysis, 1(3):385–402, 2006.

• Berger and Bernardo [1992]

J. O. Berger and J. M. Bernardo.

On the Development of Reference Priors.

Bayesian statistics, 4(4):35–60, 1992.

• Berger et al. [2001]

J. O. Berger, V. De Oliveira, and B. Sansó.

Objective Bayesian analysis of spatially correlated data.

Journal of the American Statistical Association, 96(456):1361–1374, 2001.

• Berger et al. [2015]

James O Berger, Jose M Bernardo, and Dongchu Sun.

Overall objective priors.

Bayesian Analysis, 10(1):189–221, 2015.

• Bernardo [2005]

J. M. Bernardo.

Reference analysis.

In D. Dey and C. Rao, editors, Handbook of statistics, volume 25, pages 17–90. Elsevier, 2005.

• Celeux et al. [2005]

G. Celeux, J.-M. Marin, and C. Robert.

Sélection bayésienne de variables en régression linéaire.

Journal de la société française de statistique, 147(1):59–79, 2005.

• Clarke and Barron [1994]

B. S. Clarke and A. R. Barron.

Jeffreys’ prior is asymptotically least favorable under entropy risk.

Journal of Statistical planning and Inference, 41(1):37–60, 1994.

• Dawid and Lauritzen [2001]

A. P. Dawid and S. L. Lauritzen.

Compatible prior distributions.

In E. I. George, editor, Bayesian Methods with Applications to Science, Policy, and Official Statistics, Selected Papers from ISBA 2000: The Sixth World Meeting of the International Society for Bayesian Analysis, pages 109–118. Eurostat, Luxembourg, 2001.

• Gelman and Raghunathan [2001]

A. Gelman and T. E. Raghunathan.

Comment on “Conditionally specified distributions: an introduction” by B.C. Arnold, E. Castillo and J.M. Sarabia.

Statistical Science, 16(3):268–269, 2001.

• Gu et al. [2018]

M. Gu, X. Wang, and J. O. Berger.

Robust Gaussian stochastic process emulation.

The Annals of Statistics, 46(6A):3038–3066, 2018.

• He et al. [2016]

B. D. He, C. M. De Sa, I. Mitliagkas, and C. Ré.

Scan order in Gibbs sampling: Models in which it matters and bounds on how much.

In Advances in neural information processing systems, volume 29, pages 1–9, 2016.

• Heckerman et al. [2000]

D. Heckerman, D. M. Chickering, C. Meek, R. Rounthwaite, and C. Kadie.

Dependency networks for inference, collaborative filtering, and data visualization.

Journal of Machine Learning Research, 1:49–75, 2000.

• Hobert and Casella [1998]

J. P. Hobert and G. Casella.

Functional compatibility, Markov chains, and Gibbs sampling with improper posteriors.

Journal of Computational and Graphical Statistics, 7(1):42–60, 1998.

• Hotelling [1953]

H. Hotelling.

New light on the correlation coefficient and its transforms.

Journal of the Royal Statistical Society. Series B (Methodological), 15(2):193–232, 1953.

• Journel and Huijbregts [1978]

A. G. Journel and Ch. J. Huijbregts.

Mining geostatistics.

Academic press, New York, 1978.

• Kazianka and Pilz [2012]

H. Kazianka and J. Pilz.

Objective bayesian analysis of spatial data with uncertain nugget and range parameters.

Canadian Journal of Statistics, 40(2):304–327, 2012.

• Kennedy and O’Hagan [2001]

M. C. Kennedy and A. O’Hagan.

Bayesian calibration of computer models.

Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(3):425–464, 2001.

• Kuo and Wang [2017]

K. L. Kuo and Y. J. Wang.

Pseudo-Gibbs sampler for discrete conditional distributions.

Annals of the Institute of Statistical Mathematics, pages 1–13, 2017.

• Kuo et al. [2017]

K.-L. Kuo, C.-C. Song, and T. J. Jiang.

Exactly and almost compatible joint distributions for high-dimensional discrete conditional distributions.

Journal of Multivariate Analysis, 157:115–123, 2017.

• Li and Sudjianto [2005]

R. Li and A. Sudjianto.

Analysis of computer experiments using penalized likelihood in Gaussian Kriging models.

Technometrics, 47(2):111–120, 2005.

• Mitliagkas and Mackey [2017]

I. Mitliagkas and L. Mackey.

Improving Gibbs Sampler Scan Quality with DoGS.

In Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceeding of Machine Learning Research, pages 2469–2477, 2017.

• Paulo [2005]

R. Paulo.

Default priors for Gaussian processes.

Annals of Statistics, 33(2):556–582, 2005.

• Rasmussen and Williams [2006]

C. E. Rasmussen and C. K. I. Williams.

Gaussian processes for machine learning.

MIT Press, 2006.

• Ren et al. [2012]

C. Ren, D. Sun, and C. He.

Objective bayesian analysis for a spatial model with nugget effects.

Journal of Statistical Planning and Inference, 142(7):1933–1946, 2012.

• Ren et al. [2013]

C. Ren, D. Sun, and S. K. Sahu.

Objective bayesian analysis of spatial models with separable correlation functions.

Canadian Journal of Statistics, 41(3):488–507, 2013.

• Robert [2007]

C. P. Robert.

The Bayesian Choice : From Decision-Theoretic Foundations to Computational Implementation.

Springer-Verlag, New York, 2007.

• Roverato and Consonni [2004]

A. Roverato and G. Consonni.

Compatible prior distributions for directed acyclic graph models.

Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(1):47–61, 2004.

• Santner et al. [2003]

T. J. Santner, B. J. Williams, and W. I. Notz.

The Design and Analysis of Computer Experiments.

Springer-Verlag, New York, 2003.

• Stein [1999]

M. L. Stein.

Interpolation of Spatial Data. Some Theory for Kriging.

Springer Series in Statistics. Springer-Verlag, New York, 1999.

• Yang and Berger [1996]

R. Yang and J. O. Berger.

A catalog of noninformative priors.

Institute of Statistics and Decision Sciences, Duke University, 1996.

• Yet and Marsh [2014]

B. Yet and W. Marsh.

Compatible and incompatible abstractions in Bayesian networks.

Knowledge-Based Systems, 62:84–97, 2014.

• Zhang [2004]

H. Zhang.

Inconsistent Estimation and Asymptotically Equal Interpolations in Model-based Geostatistics.

Journal of the American Statistical Association, 99(465):250–261, 2004.