Bayesian Spatial Inversion and Conjugate Selection Gaussian Prior Models

TL;DR

This paper develops conjugate prior models for Bayesian spatial inversion, introducing selection Gaussian priors that handle complex data features and demonstrate improved seismic inversion accuracy over traditional methods.

Contribution

It introduces selection Gaussian prior models that extend conjugate priors, enabling more flexible modeling of spatial data with multi-modality, skewness, and heavy tails.

Findings

Selection Gaussian priors are conjugate for Gauss-linear likelihoods.

Posterior models remain selection Gaussian with explicit parameter relationships.

Seismic inversion shows 20-40% improvement in mean-square-error.

Abstract

We introduce the concept of conjugate prior models for a given likelihood function in Bayesian spatial inversion. The conjugate class of prior models can be selection extended and still remain conjugate. We demonstrate the generality of selection Gaussian prior models, representing multi-modality, skewness and heavy-tailedness. For Gauss-linear likelihood functions, the posterior model is also selection Gaussian. The model parameters of the posterior pdf are explisite functions of the model parameters of the likelihood and prior models - and the actual observations, of course. Efficient algorithms for simulation of and prediction for the selection Gaussian posterior pdf are defined. Inference of the model parameters in the selection Gaussian prior pdf, based on one training image of the spatial variable, can be reliably made by a maximum likelihood criterion and numerical optimization. Lastly, a seismic inversion case study is presented, and improvements of $ 20$-$40\%$ in prediction mean-square-error, relative to traditional Gaussian inversion, are found.

Tables

Table 1. Table 1 : Model parameters for six cases, with μ = 0 𝜇 0 \mu=0 and σ 2 = 1 superscript 𝜎 2 1 \sigma^{2}=1 for all cases.

Case	$\gamma$	$d_h$	$d_v$	$A_i$	description
1	0.8000	2.0	2.0	$(\infty ,-0.3]\cup [0.3,\infty )$	sym. bimodal iso.
2	0.6500	6.0	0.85	$(\infty ,-0.3]\cup [0.3,\infty )$	sym. bimodal aniso.
3	0.9250	2.0	0.60	$(\infty ,-0.85]\cup [0.8,\infty )$	asym. bimodal aniso.
4	0.9995	3.0	3.0	$[-0.45,-0.2]\cup [-0.1,0.1]\cup [0.2,0.45]$	sym. trimodal iso.
5	0.7000	2.0	2.0	$(\infty ,-0.7]\cup [-0.1,2.5]$	asym. unimodal iso.
6	0.7000	2.0	2.0	$(\infty ,-1.75]\cup [-0.5,0.5]\cup [1.75,\infty )$	sym. heavy tailed iso.

Table 2. Table 2 : Model parameters for four posterior cases, with μ = 0 𝜇 0 \mu=0 and σ 2 = 1 superscript 𝜎 2 1 \sigma^{2}=1 for all cases.

Case	$\gamma$	$d_h$	$A_i$	description	cond. values
1	0.900	4	$(\infty ,-0.4]\cup [0.4,\infty )$	sym. bimodal	$2.5,-2.5$
2	0.999	4	$[-0.65,-0.4]\cup [0.12,0.12]\cup \left[\mathrm{0.40.65}\right]$	sym. trimodal	$0.55,-0.55$
3	0.600	4	$(\infty ,-1.5]\cup [-0.5,0.5)$	asym. unimodal	$1.0,-3.0$
4	0.700	4	$(\infty ,-1.75]\cup [-0.5,0.5]\cup [1.75,\infty )$	sym. heavy tailed	$3.0,-3.0$

Table 3. Table 3 : Summary of well predictions for the selection Gaussian and Gaussian model. Mean square error (MSE) of predictions, posterior and prior 80% coverage of prediction intervals.

	MSE		Prior 80% coverage		Posterior 80% coverage
	Selection	Gaussian	Selection	Gaussian	Selection	Gaussian
$\log v_p$	0.0034	0.0050	0.84	0.88	0.85	0.96
$\log v_s$	0.0112	0.0191	0.82	0.89	0.84	0.87
$\log \rho$	0.0009	0.0011	0.82	0.95	0.83	0.89

Table 4. Table 4 : Model parameters for six cases, with μ = 0 𝜇 0 \mu=0 and σ 2 = 1 superscript 𝜎 2 1 \sigma^{2}=1 for all cases.

Case	$\gamma$	$d_h$	$d_v$	$A_i$	description
1	0.8000	2.0	2.0	$(\infty ,-0.3]\cup [0.3,\infty )$	sym. bimodal iso.
2	0.6500	6.0	0.85	$(\infty ,-0.3]\cup [0.3,\infty )$	asym. bimodal aniso.
3	0.9250	2.0	0.60	$(\infty ,-0.85]\cup [0.8,\infty )$	sym. bimodal aniso.
4	0.9995	3.0	3.0	$[-0.45,-0.2]\cup [-0.1,0.1]\cup [0.2,0.45]$	sym. trimodal iso.
5	0.7000	2.0	2.0	$(\infty ,-0.7]\cup [-0.1,2.5]$	asym. unimodal iso.
6	0.7000	2.0	2.0	$(\infty ,-1.75]\cup [-0.5,0.5]\cup [1.75,\infty )$	sym. heavy tailed iso.

Table 5. Table 5 : Model parameters for four posterior cases, with μ = 0 𝜇 0 \mu=0 and σ 2 = 1 superscript 𝜎 2 1 \sigma^{2}=1 for all cases.

Case	$\gamma$	$d_h$	$A_i$	description	cond. values
1	0.900	4	$(\infty ,-0.4]\cup [0.4,\infty )$	sym. bimodal	$2.5,-2.5$
2	0.999	4	$[-0.65,-0.4]\cup [0.12,0.12]\cup \left[\mathrm{0.40.65}\right]$	sym. trimodal	$0.55,-0.55$
3	0.600	4	$(\infty ,-1.5]\cup [-0.5,0.5)$	asym. unimodal	$1.0,-3.0$
4	0.700	4	$(\infty ,-1.75]\cup [-0.5,0.5]\cup [1.75,\infty )$	sym. heavy tailed	$3.0,-3.0$

Equations

$$[\mathbf{r}|\mathbf{d}]\to f(\mathbf{r}|\mathbf{d})$$

$$=\text{mbox}const\times f(\mathbf{d}|\mathbf{r})f(\mathbf{r})$$

$$\mathbf{y}\to f(\mathbf{y})$$

$$=[2\pi {]}^{-n/2}|\mathbf{\Sigma }{|}^{-1/2}\exp \{-\frac{1}{2}[\mathbf{y}-\mu {]}^T{\mathbf{\Sigma }}^{-1}[\mathbf{y}-\mu ]\}$$

$$F(\mathbf{y}\in \mathcal{B})$$

$$f(\mathbf{r}|\mathbf{d};{\theta }_{r|d})=\text{mbox}const\times f(\mathbf{d}|\mathbf{r};{\psi }_d)f(\mathbf{r};{\theta }_r)$$

$$\displaystyle\left[\begin{array}[]{c}\boldsymbol{\mathbf{r}}\\ \boldsymbol{\mathbf{\nu}}\end{array}\right]\rightarrow f\left(\left[\begin{array}[]{c}\boldsymbol{\mathbf{r}}\\ \boldsymbol{\mathbf{\nu}}\end{array}\right]\right)=f(\boldsymbol{\mathbf{\nu}}|\boldsymbol{\mathbf{r}})f(\boldsymbol{\mathbf{r}})$$

$${\mathbf{r}}_A=[\mathbf{r}|\nu \in \mathcal{A}]\to f({\mathbf{r}}_A)$$

$$=[F(\nu \in \mathcal{A}){]}^{-1}\times F(\nu \in \mathcal{A}|\mathbf{r})f(\mathbf{r})$$

$${\mathbf{r}}_A=[\mathbf{r}|\nu \in \mathcal{A}]\to f({\mathbf{r}}_A)$$

$$=[F(\nu \in \mathcal{A};{\theta }_r,{\kappa }_{\nu }){]}^{-1}\times F(\nu \in \mathcal{A}|\mathbf{r};{\kappa }_{\nu })f(\mathbf{r};{\theta }_r)$$

$$\mathbf{r}\to f(\mathbf{r})={\varphi }_{n_r}(\mathbf{r};{\mu }_r,{\mathbf{\Sigma }}_r)$$

$$[\nu |\mathbf{r}]\to f(\nu |\mathbf{r})={\varphi }_{n_{\nu }}(\nu ;{\mu }_{\nu |r},{\mathbf{\Sigma }}_{\nu |r})$$

$${\mathbf{r}}_A=[\mathbf{r}|\nu \in \mathcal{A}]\to f({\mathbf{r}}_A)$$

$$=[{\Phi }_{n_{\nu }}(\mathcal{A};{\mu }_{\nu },{\mathbf{\Sigma }}_{\nu }){]}^{-1}\times {\Phi }_{n_{\nu }}(\mathcal{A};{\mu }_{\nu |r},{\mathbf{\Sigma }}_{\nu |r}){\varphi }_{n_r}(\mathbf{r};{\mu }_r,{\mathbf{\Sigma }}_r)$$

$$[{\mathbf{r}}_A|{\mathbf{d}}^o]\to f({\mathbf{r}}_A|{\mathbf{d}}^o)$$

$$=\text{mbox}{const}_1\times f({\mathbf{d}}^o|\mathbf{r})f(\mathbf{r}|\nu \in \mathcal{A})$$

$$=\text{mbox}{const}_2\times f({\mathbf{d}}^o|\mathbf{r})F(\nu \in \mathcal{A}|\mathbf{r})f(\mathbf{r})$$

$$=\text{mbox}{const}_3\times F(\nu \in \mathcal{A}|\mathbf{r},{\mathbf{d}}^o)f(\mathbf{r}|{\mathbf{d}}^o)$$

$$=[F(\nu \in \mathcal{A}|{\mathbf{d}}^o){]}^{-1}\times F(\nu \in \mathcal{A}|\mathbf{r},{\mathbf{d}}^o)f(\mathbf{r}|{\mathbf{d}}^o)$$

$$[\mathbf{d}|\mathbf{r}]=\mathbf{H}\mathbf{r}+{\mathbf{e}}_{d|r}\to f(\mathbf{d}|\mathbf{r})={\varphi }_{n_d}(\mathbf{d};\mathbf{H}\mathbf{r},{\mathbf{\Sigma }}_{d|r})$$

$$\mathbf{r}\to f(\mathbf{r})={\varphi }_{n_r}(\mathbf{r};\mu {\mathbf{i}}_{n_r},{\sigma }^2\mathbf{C})$$

$$[\nu |\mathbf{r}]=\gamma {\sigma }^{-1}[\mathbf{r}-\mu {\mathbf{i}}_{n_r}]+{\mathbf{e}}_{\nu |r}\to f(\nu |\mathbf{r})$$

$$=\prod _{i=1}^{n_r}{\varphi }_1({\nu }_i;\gamma {\sigma }^{-1}[r_i-\mu ],[1-{\gamma }^2])$$

$$\displaystyle\left[\begin{array}[]{c}\boldsymbol{\mathbf{r}}\\ \boldsymbol{\mathbf{\nu}}\end{array}\right]\rightarrow f\left(\left[\begin{array}[]{c}\boldsymbol{\mathbf{r}}\\ \boldsymbol{\mathbf{\nu}}\end{array}\right]\right)=\phi_{2n_{r}}\left(\left[\begin{array}[]{c}\boldsymbol{\mathbf{r}}\\ \boldsymbol{\mathbf{\nu}}\end{array}\right];\left[\begin{array}[]{c}\mu\boldsymbol{\mathbf{i}}_{n_{r}}\\ 0\boldsymbol{\mathbf{i}}_{n_{r}}\end{array}\right],\left[\begin{array}[]{cc}\sigma^{2}\boldsymbol{\mathbf{C}}&\gamma\sigma\boldsymbol{\mathbf{C}}\\ \gamma\sigma\boldsymbol{\mathbf{C}}&\gamma^{2}\boldsymbol{\mathbf{C}}+[1-\gamma^{2}]\boldsymbol{\mathbf{I}}_{n_{r}}\end{array}\right]\right)$$

$${\mathbf{r}}_A=[\mathbf{r}|\nu \in \mathcal{A}]\to$$

$$=[F(\nu \in \mathcal{A}){]}^{-1}\times F(\nu \in \mathcal{A}|\mathbf{r})f(\mathbf{r})$$

$$=[{\Phi }_{n_r}({\cup }_{i=1}^{n_r}{\mathcal{A}}_i;0{\mathbf{i}}_{n_r},{\gamma }^2\mathbf{C}+[1-{\gamma }^2]{\mathbf{I}}_{n_r}){]}^{-1}$$

$$\times \prod _{i=1}^{n_r}{\Phi }_1({\mathcal{A}}_i;\gamma {\sigma }^{-1}[r_i-\mu ],[1-{\gamma }^2])\text{mbox}{\varphi }_{n_r}(\mathbf{r};\mu {\mathbf{i}}_{n_r},{\sigma }^2\mathbf{C})$$

$$r_{Ai}=[r_i|\nu \in \mathcal{A}]\to$$

$$=[F(\nu \in \mathcal{A}){]}^{-1}\times F(\nu \in \mathcal{A}|r_i)f(r_i)$$

$$=[{\Phi }_{n_r}({\cup }_{i=1}^{n_r}{\mathcal{A}}_i;0{\mathbf{i}}_{n_r},{\gamma }^2\mathbf{C}+[1-{\gamma }^2]{\mathbf{I}}_{n_r}){]}^{-1}$$

$$\times {\Phi }_{n_r}({\cup }_{i=1}^{n_r}{\mathcal{A}}_i;\gamma {\sigma }^{-1}{\mathbf{c}}_i(r_i-\mu ),{\gamma }^2[\mathbf{C}-{\mathbf{c}}_i{\mathbf{c}}_i^T]+[1-{\gamma }^2]{\mathbf{I}}_{n_r})$$

$$\times {\varphi }_1(r_i;\mu ,{\sigma }^2)$$

$$\displaystyle\left[\begin{array}[]{c}\boldsymbol{\mathbf{r}}\\ \boldsymbol{\mathbf{\nu}}\\ \boldsymbol{\mathbf{d}}\end{array}\right]\rightarrow f\left(\left[\begin{array}[]{c}\boldsymbol{\mathbf{r}}\\ \boldsymbol{\mathbf{\nu}}\\ \boldsymbol{\mathbf{d}}\end{array}\right]\right)$$

$$\displaystyle=\phi_{2n_{r}+n_{d}}\left(\left[\begin{array}[]{c}\boldsymbol{\mathbf{r}}\\ \boldsymbol{\mathbf{\nu}}\\ \boldsymbol{\mathbf{d}}\end{array}\right];\left[\begin{array}[]{c}\mu\boldsymbol{\mathbf{i}}_{n_{r}}\\ 0\boldsymbol{\mathbf{i}}_{n_{r}}\\ \mu\boldsymbol{\mathbf{H}}\boldsymbol{\mathbf{i}}_{n_{r}}\end{array}\right],\left[\begin{array}[]{ccc}\sigma^{2}\boldsymbol{\mathbf{C}}&\gamma\sigma\boldsymbol{\mathbf{C}}&\sigma^{2}\boldsymbol{\mathbf{C}}\boldsymbol{\mathbf{H}}^{T}\\ \gamma\sigma\boldsymbol{\mathbf{C}}&\gamma^{2}\boldsymbol{\mathbf{C}}+[1-\gamma^{2}]\boldsymbol{\mathbf{I}}_{n_{r}}&\gamma\sigma\boldsymbol{\mathbf{C}}\boldsymbol{\mathbf{H}}^{T}\\ \sigma^{2}\boldsymbol{\mathbf{H}}\boldsymbol{\mathbf{C}}&\gamma\sigma\boldsymbol{\mathbf{H}}\boldsymbol{\mathbf{C}}&\sigma^{2}\boldsymbol{\mathbf{H}}\boldsymbol{\mathbf{C}}\boldsymbol{\mathbf{H}}^{T}+\boldsymbol{\mathbf{\Sigma}}_{d|r}\end{array}\right]\right)$$

$$[{\mathbf{r}}_A|\mathbf{d}]=[\mathbf{r}|\nu \in \mathcal{A},\mathbf{d}]\to$$

$$=[F(\nu \in \mathcal{A}|\mathbf{d}){]}^{-1}\times F(\nu \in \mathcal{A}|\mathbf{r},\mathbf{d})f(\mathbf{r}|\mathbf{d})$$

$$=[{\Phi }_{n_r}({\cup }_{i=1}^{n_r}{\mathcal{A}}_i;{\mu }_{\nu |d},{\mathbf{\Sigma }}_{\nu |d}){]}^{-1}$$

$$\times {\Phi }_{n_r}({\cup }_{i=1}^{n_r}{\mathcal{A}}_i;{\mu }_{\nu |rd},{\mathbf{\Sigma }}_{\nu |rd})$$

$$\times {\varphi }_{n_r}(\mathbf{r};{\mu }_{r|d},{\mathbf{\Sigma }}_{r|d})$$

$$\displaystyle\left[\begin{array}[]{c}\boldsymbol{\mathbf{\mu}}_{r|d}\\ \boldsymbol{\mathbf{\mu}}_{\nu|d}\end{array}\right]=\left[\begin{array}[]{c}\mu\boldsymbol{\mathbf{i}}_{n_{r}}\\ 0\boldsymbol{\mathbf{i}}_{n_{r}}\end{array}\right]+\left[\begin{array}[]{c}\sigma^{2}\boldsymbol{\mathbf{C}}\boldsymbol{\mathbf{H}}^{T}\\ \gamma\sigma\boldsymbol{\mathbf{C}}\boldsymbol{\mathbf{H}}^{T}\end{array}\right]\left[\sigma^{2}\boldsymbol{\mathbf{H}}\boldsymbol{\mathbf{C}}\boldsymbol{\mathbf{H}}^{T}+\boldsymbol{\mathbf{\Sigma}}_{d|r}\right]^{-1}\left[\boldsymbol{\mathbf{d}}-\mu\boldsymbol{\mathbf{H}}\boldsymbol{\mathbf{i}}_{n_{r}}\right]$$

Full text

Keywords: Inverse Problems, Spatial Prediction, Conditional Simulation, Spatial Inference, Seismic Inversion

1 Introduction

Inversion constitutes a challenge in many mathematical engineering problems. Observations from the variable of interest are often indirectly collected by some complex acquisition devise. The objective is naturally to predict the variable of interest based on the available observations. We consider spatial variables in this study and examples of inverse problems can be found in image analysis, remote sensing and geophysics. Inversion of seismic data is presented as a case study later in the paper.

The spatial variable of interest is $\{r(\boldsymbol{\mathbf{x}});\boldsymbol{\mathbf{x}}\in\mathcal{D}\subset\mathcal{R}^{m}\}$ with $r(\boldsymbol{\mathbf{x}})\in\mathcal{R}$ being the variable, having spatial reference $\boldsymbol{\mathbf{x}}$ running over the reference domain $\mathcal{D}\subset\mathcal{R}^{m}$ , which naturally has dimension $m$ equal to one, two or three. The variable is discretized $\{r(\boldsymbol{\mathbf{x}});\boldsymbol{\mathbf{x}}\in\mathcal{L}_{\mathcal{D}}\}$ where $\mathcal{L}_{\mathcal{D}}$ is a regular grid, of size $n_{r}$, covering $\mathcal{D}$ , and the spatial variable is represented by the $n_{r}$-vector $\boldsymbol{\mathbf{r}}\in\mathcal{R}^{n_{r}}$. Assume further that a $n_{d}$-vector of observations $\boldsymbol{\mathbf{d}}\in\mathcal{R}^{n_{d}}$ related to the variable of interest is collected. The focus is on assessing $\boldsymbol{\mathbf{r}}$ given $\boldsymbol{\mathbf{d}}$, $[\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}}]$.

We phrase the assessment in a probabilistic setting, by using Bayesian inversion, see Tarantola (2005),

[TABLE]

where $\boldsymbol{\mathbf{y}}\rightarrow f(\boldsymbol{\mathbf{y}})$ reads, the random variable $\boldsymbol{\mathbf{y}}$ is distributed according to the probability density function (pdf) $f(\boldsymbol{\mathbf{y}})$. The $f(\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}})$ is the posterior pdf being the ultimate solution of Bayesian inversion. The likelihood function $f(\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{r}})$ , being a function of $\boldsymbol{\mathbf{r}}$, represents the observation acquisition procedure, while the prior pdf $f(\boldsymbol{\mathbf{r}})$ summarizes prior information about the spatial variable of interest. The likelihood and prior models uniquely define the posterior model, although the integral in the normalizing constant usually is complicated to calculate.

Classical Bayesian inference, see Casella and Berger (2002), conserns estimation of model parameters, contrary to Bayesian inversion which is defined in a predictive setting. The classical approach focus on the posterior pdf of a low-dimensional vector of model parameters $\boldsymbol{\mathbf{\theta}}$ given a set of observations $\boldsymbol{\mathbf{y}}$ , hence on $f(\boldsymbol{\mathbf{\theta}}|\boldsymbol{\mathbf{y}})$. To assess this posterior pdf one must assume a likelihood function $f(\boldsymbol{\mathbf{y}}|\boldsymbol{\mathbf{\theta}})$ and specify a prior pdf $f(\boldsymbol{\mathbf{\theta}})$. In order to avoid to calculate complex integrals, classical Bayesian inference has introduced the concept of conjugate classes of parametric prior pdfs. For a given likelihood function $f(\boldsymbol{\mathbf{y}}|\boldsymbol{\mathbf{\theta}})$ , with a prior model $f(\boldsymbol{\mathbf{\theta}})$ from the corresponding conjugate class of pdfs, the resulting posterior pdf $f(\boldsymbol{\mathbf{\theta}}|\boldsymbol{\mathbf{y}})$ will belong to the same conjugate class of pdfs. Consequently, the hyper-parameters of the posterior pdf will depend only on the hyper-parameters of the likelihood function and the prior pdf, and the actual observations, hence complex integral calculations are avoided.

In Bayesian spatial inversion the integral calculations are even more challenging than in classical Bayesian inference, since the variable of interest is of much higher dimensions than the model parameters. Fortunately, the elicitation of the prior model is simpler in Bayesian inversion than in Bayesian inference since the prior is on observable variables in the former while it is on model parameters in the latter. Hence Bayesian inversion may naturally be cast in an empirical Bayesian setting, see Efron and Morris (1973) . In this study we consider Bayesian inversion and introduce the consept of conjugate classes of prior parametrised pdfs for a given class of likelihood functions. Further, we demonstrate that, for a given class of likelihood functions, the corresponding class of conjugate prior pdfs can be generalized by a selection mechanism. The resulting selection class of pdfs will also be conjugate with respect to the same class of likelihood functions. This construction makes it possible to define highly flexible classes of conjugate prior pdfs where the corresponding posterior pdfs can be determined by only the associated hyper-parameters and the actual observations. The assessment of the resulting posterior pdfs may require numerical or simulation schemes, but the closed form expressions for the posterior pdf makes it possible to design tailored efficient schemes. We demonstrate the selection conjugate consept on the familiar case of Gauss-linear likelihood functions and the conjugate class of spatial Gaussian pdfs.

Bayesian predictive inversion and Bayesian model inference can be combined to have hierarchical Bayesian inversion. It is, however, complicated to define general conjugate classes of prior models in this setting, see Røislien and Omre (2006), Arellano-Valle et al. (2009) and Branco et al. (2013). We present a brief discussion and evaluation of classical likelihood inference of the model parameters in the selection Gaussian prior pdf, from one available training image. Lastly, a case study of seismic inversion of real data is presented.

The developments of the selection Gaussian prior model is inspired by the early work on skewed pdfs in Azzalini (1985) and Azzalini and Dalla Valle (1996), see also Genton (2004) and Azzalini (2013). These models are generalized to spatial settings in Kim and Mallick (2004), Allard and Naveau (2007), and Rimstad and Omre (2014). In the current study we use generalized selection sets as discussed in Arellano-Valle and del Pino (2004) and Arellano-Valle et al. (2006), to model spatial prior pdfs with marginal multi-modality. Also marginal skewness and heavy-tailedness may be represented. In spatial modelling, multi-modal spatial histograms occur if latent categorical variables are present. These variables may be lithologies in subsurface reservoir modelling or tissue classes in medical image analysis. Modelling multi-modal marginal characteristics with spatial continuity is challenging. If the observations are collected with spatial convolution and errors, as in subsurface acquisition and medical imaging, the inversion challenge is even larger.

The major contribution of the current paper is, however, the discussion of conjugate prior models in Bayesian spatial inversion and the demonstration that this conjugate characteristic is closed under activation of a selection mechanism. This result holds for all types of conjugate prior models. Moreover, we demonstrate the large potential of this result in modelling continuous spatial variables with multi-modal marginal distributions.

In the presentation $f(\boldsymbol{\mathbf{y}})$ denotes a pdf of the random variable $\boldsymbol{\mathbf{y}}$, while $F(\boldsymbol{\mathbf{y}}\in\mathcal{B})$ denotes the probability that the random variable is in the sub-set $\mathcal{B}$ of its sample space. For the Gaussian random $n$-vector $\boldsymbol{\mathbf{y}}$ we write,

[TABLE]

Note that from a computational point of view $|\boldsymbol{\mathbf{\Sigma}}|^{-1}$ and $\boldsymbol{\mathbf{\Sigma}}^{-1}$ are demanding, and so is simulation from and calculation of the sub-set $\mathcal{B}$ probability of an arbitrary high-dimensional Gaussian pdf. The former two are widely studied, however, while the latter has drawn much less attention except for Genz and Bretz (2009). Later, we suggest efficient algorithms for simulation from and calculating of these sub-set $\mathcal{B}$ probabilities. We also use the notation $\boldsymbol{\mathbf{i}}_{n}$ for a unit $n$-vector and $\boldsymbol{\mathbf{I}}_{n}$ for a identity $(n\times n)$-matrix, while $I(A)$ is an indicator function taking value $1$ if $A$ is true and [math] otherwise.

In Section 2, Bayesian spatial inversion is discussed, and a selection extended conjugate class of prior pdfs for a given class of likelihood functions is defined. The conjugate class of selection Gaussian prior pdfs is developed and discussed in detail. Expressions for the model parameters of the corresponding posterior pdf are developed, and simulation algorithms for assessing the posterior pdf are presented. The flexibility of the class of prior pdfs is illustrated by several examples. Section 3 contains a discussion of model parameter inference and the development of maximum likelihood estimators for the model parameters of the conjugate class of selection Gaussian prior pdfs, based on one training realization of the spatial variable. A small empirical study demonstrates the consistency of the estimator as the size of the training variable increase. In Section 4, a case study based on real seismic data along a well profile from the Alvheim field in the North Sea is presented, see also Karimi et al. (2010). Comparisons with regular Bayesian Gaussian inversion are made. Lastly, in Section 5, the conclusions of the study are forwarded.

2 Bayesian Spatial Inversion

The focus is on prediction of a spatial variable discretized into the $n_{r}$-vector $\boldsymbol{\mathbf{r}}$, based on the observations represented in the $n_{d}$-vector $\boldsymbol{\mathbf{d}}$. We phrase the prediction as Bayesian inversion, see Expression 1, which requires that the likelihood function $f(\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{r}})$ is given, and that the prior pdf $f(\boldsymbol{\mathbf{r}})$ is specified. Hereby, the corresponding posterior pdf $f(\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}})$ is defined. Inspired by the consept of conjugate prior pdfs in traditional Bayesian inference we present the following definition,

Definition 1** (Conjugate class of prior pdfs)**

Consider Bayesian inversion,

[TABLE]

with likelihood function $f(\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{r}};\boldsymbol{\mathbf{\psi}}_{d})$ in a parametrized pdf class $\mathcal{L}_{\psi}$ and prior pdf $f(\boldsymbol{\mathbf{r}};\boldsymbol{\mathbf{\theta}}_{r})$ in a parametrized pdf class $\mathcal{P}_{\theta}$. If the associated posterior pdf $f(\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}};\boldsymbol{\mathbf{\theta}}_{r|d})$ also is in the pdf class $\mathcal{P}_{\theta}$, then the pdf class $\mathcal{P}_{\theta}$ is termed a conjugate class with respect to the likelihood function class $\mathcal{L}_{\psi}$. The posterior model parameters $\boldsymbol{\mathbf{\theta}}_{r|d}$ will be a function of $[\boldsymbol{\mathbf{\psi}}_{d},\boldsymbol{\mathbf{\theta}}_{r},\boldsymbol{\mathbf{d}}]$.

For continuous spatial variables, the class of Gaussian prior pdfs is known to be a conjugate class with respect to Gauss-linear likelihood functions. Hence, if the observations are collected through a linear forward model with additive Gaussian errors, and the prior pdf is specified to be Gaussian, then the posterior pdf will also be Gaussian. This characteristic is the basis for kriging prediction and conditional simulation in geostatistics, see Chiles and Delfiner (2012). Moreover, for event spatial variables the Poisson prior pdf is conjugate with respect to thinning likelihood functions, while for mosaic spatial variables the Markov prior pdf is conjugate with respect to conditionally independent single-site response likelihood functions. These conjugate characteristics are of course the major reason for the frequent use of these spatial models. In the next section we define an extended class of prior pdfs based on a selection consept, and demonstrate that this consept can be used to construct an extended conjugate class of prior pdfs.

2.1 Generalization by Selection

Consider the previously defined spatial variable represented by the $n_{r}$-vector $\boldsymbol{\mathbf{r}}$ with prior pdf $f(\boldsymbol{\mathbf{r}})$. Extend the dimension by an auxiliary random $n_{\nu}$-vector $\boldsymbol{\mathbf{\nu}}\in\mathcal{R}^{n_{\nu}}$, such that,

[TABLE]

with arbitrary chosen pdf $f(\boldsymbol{\mathbf{\nu}}|\boldsymbol{\mathbf{r}})$, and denote the pdf $f(\boldsymbol{\mathbf{r}})$ the basis-pdf. Consider an arbitrary sub-set $\mathcal{A}\subset\mathcal{R}^{n_{\nu}}$ and define the associated random selection $n_{r}$-vector $\boldsymbol{\mathbf{r}}_{A}$ by,

[TABLE]

Note in particular, that $f(\boldsymbol{\mathbf{r}}_{A})=f(\boldsymbol{\mathbf{r}})$ if we define $f(\boldsymbol{\mathbf{\nu}}|\boldsymbol{\mathbf{r}})=f(\boldsymbol{\mathbf{\nu}})$, which of course is the extreme choice of independence between $\boldsymbol{\mathbf{r}}$ and $\boldsymbol{\mathbf{\nu}}$. Based on this selection consept we define,

Definition 2** (Selection extension of prior pdfs)**

Consider a prior basis-pdf $f(\boldsymbol{\mathbf{r}};\boldsymbol{\mathbf{\theta}}_{r})$ in a parametrized pdf class $\mathcal{P}_{\theta}$ and define auxiliary variable $\boldsymbol{\mathbf{\nu}}\in\mathcal{R}^{n_{\nu}}$ related to $\boldsymbol{\mathbf{r}}$ by pdf $f(\boldsymbol{\mathbf{\nu}}|\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{\kappa}}_{\nu})$ in a parametrized pdf class $\mathcal{E}_{\kappa}$. Specify further a selection set $\mathcal{A}\subset\mathcal{R}^{n_{\nu}}$. Define the selection variable,

[TABLE]

with pdf $f(\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{\nu}}\in\mathcal{A};\boldsymbol{\mathbf{\theta}}_{r},\boldsymbol{\mathbf{\kappa}}_{\nu})$ in the parametrized selection extended pdf class $\mathcal{S}_{\mathcal{A}}[\mathcal{P}_{\theta}\times\mathcal{E}_{\kappa}]$.

The selection extension can be made for any basis-pdf class for arbitrary auxiliary variables with associated selection sets. The class of selection Gaussian pdfs with $f(\boldsymbol{\mathbf{r}})$ from the Gaussian class and $f(\boldsymbol{\mathbf{\nu}}|\boldsymbol{\mathbf{r}})$ being Gauss-linear with associated selection sets, hence $[\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{\nu}}]$ being jointly Gaussian, is thoroughly discussed in Arellano-Valle et al. (2006). We define this class of selection Gaussian pdf by a Gaussian basis-pdf

[TABLE]

with the expectation $n_{r}$-vector $\boldsymbol{\mathbf{\mu}}_{r}$ and the covariance $(n_{r}\times n_{r})$-matrix $\boldsymbol{\mathbf{\Sigma}}_{r}$. The auxiliary $n_{\nu}$-vector $\boldsymbol{\mathbf{\nu}}$ is defined as,

[TABLE]

with the conditional expectation $n_{\nu}$-vector being linear in $\boldsymbol{\mathbf{r}}$, $\boldsymbol{\mathbf{\mu}}_{\nu|r}=\boldsymbol{\mathbf{\mu}}_{\nu}+\boldsymbol{\mathbf{\Gamma}}_{\nu|r}(\boldsymbol{\mathbf{r}}-\boldsymbol{\mathbf{\mu}}_{r})$ with expectation $n_{\nu}$-vector $\boldsymbol{\mathbf{\mu}}_{\nu}$, coupling $(n_{\nu}\times n_{r})$-matrix $\boldsymbol{\mathbf{\Gamma}}_{\nu|r}$, and the conditional covariance $(n_{r}\times n_{\nu})$-matrix $\boldsymbol{\mathbf{\Sigma}}_{\nu|r}$. Hence $[\boldsymbol{\mathbf{\nu}}|\boldsymbol{\mathbf{r}}]$ is Gauss-linear, and since $\boldsymbol{\mathbf{r}}$ is Gaussian, the joint $(n_{r}+n_{\nu})$-vector $[\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{\nu}}]$ is Gaussian. By enforcing the selection $\boldsymbol{\mathbf{\nu}}\in\mathcal{A}\subset\mathcal{R}^{n_{\nu}}$, we obtain the selection Gaussian $n_{r}$-vector $\boldsymbol{\mathbf{r}}_{A}$,

[TABLE]

where the marginal covariance $(n_{\nu}\times n_{\nu})$-matrix is $\boldsymbol{\mathbf{\Sigma}}_{\nu}=\boldsymbol{\mathbf{\Gamma}}_{\nu|r}\boldsymbol{\mathbf{\Sigma}}_{r}\boldsymbol{\mathbf{\Gamma}}_{\nu|r}^{T}+\boldsymbol{\mathbf{\Sigma}}_{\nu|r}$. All valid sets of model parameters $(\boldsymbol{\mathbf{\mu}}_{r},\boldsymbol{\mathbf{\Sigma}}_{r},\boldsymbol{\mathbf{\mu}}_{\nu},\boldsymbol{\mathbf{\Gamma}}_{\nu|r},\boldsymbol{\mathbf{\Sigma}}_{\nu|r},\mathcal{A})$, define the class of selection Gaussian pdfs. Note, that by assigning $\boldsymbol{\mathbf{\Gamma}}_{\nu|r}$ a null-matrix or setting $\mathcal{A}=\mathcal{R}^{n_{\nu}}$, the selection Gaussian class of pdfs is identical to the Gaussian one.

The likelihood function of the actual observations $\boldsymbol{\mathbf{d}}^{o}$, $f(\boldsymbol{\mathbf{d}}^{o}|\boldsymbol{\mathbf{r}})$, is a function of $\boldsymbol{\mathbf{r}}$ and only dependent on the observation acquisition procedure, independent on choice of prior pdf. Hence, the posterior pdf based on a prior selection pdf, can be expressed as,

[TABLE]

which relies on the conditional independence relation $f(\boldsymbol{\mathbf{\nu}},\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{r}})=f(\boldsymbol{\mathbf{\nu}}|\boldsymbol{\mathbf{r}})f(\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{r}})$. Note that the posterior pdf corresponds to the selection pdf with basis-pdf being the conditional pdf $f(\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}})$, which provides the following theorem,

Theorem 1** (Selection extended conjugate class of prior pdfs)**

Consider a likelihood function in parametrized pdf class $\mathcal{L}_{\psi}$ and a prior pdf in the parametrized pdf class $\mathcal{P}_{\theta}$ . According to Definition 1 - let the pdf class $\mathcal{P}_{\theta}$ be a conjugate class with respect to likelihood function class $\mathcal{L}_{\psi}$. According to Definition 2 - define the associated selection extended pdf class $\mathcal{S}_{\mathcal{A}}[\mathcal{P}_{\theta}\times\mathcal{E}_{\kappa}]$ based on auxiliary pdf class $\mathcal{E}_{\kappa}$ and selection set $\mathcal{A}$. Then the pdf class $\mathcal{S}_{\mathcal{A}}[\mathcal{P}_{\theta}\times\mathcal{E}_{\kappa}]$ is a conjugate class with respect to the likelihood function class $\mathcal{L}_{\psi}$ for all pdf classes $\mathcal{E}_{\kappa}$ and selection sets $\mathcal{A}$. The conjugate characteristics of a prior pdf class $\mathcal{P}_{\theta}$ is closed under selection extension.

This closedness property for conjugate pdf classes is very general and applies to continuous, event and mosaic spatial variables. Moreover, it may be used in traditional Bayesian inference. For continuous spatial variables with a Gauss-linear likelihood function, prior pdfs from the class of Gaussian pdfs is known to be conjugate. According to the results above, also the selection Gaussian pdf for any arbitray auxiliary extension and selection set, will define a conjugate class of pdfs with repect to a Gauss-linear likelihood function. In the following sub-sections we will explore this opportunity to define more flexible prior pdfs in Bayesian spatial inversion with observations collected through a Gauss-likelihood function.

2.2 Likelihood Model

We limit the likelihood function to be from the Gauss-linear class,

[TABLE]

where $\boldsymbol{\mathbf{H}}$ is an observation acquisition $(n_{d}\times n_{r})$-matrix and $\boldsymbol{\mathbf{e}}_{d|r}$ is a centred Gaussian error $n_{d}$-vector with covariance $(n_{d}\times n_{d})$-matrix $\boldsymbol{\mathbf{\Sigma}}_{d|r}$, independent of $\boldsymbol{\mathbf{r}}$. Hence the model parameters are $\boldsymbol{\mathbf{\theta}}_{l}=(\boldsymbol{\mathbf{H}},\boldsymbol{\mathbf{\Sigma}}_{d|r})$. There are no constraint on the matrix $\boldsymbol{\mathbf{H}}$. It may be binary, selecting only a sub-set of elements in the vector $\boldsymbol{\mathbf{r}}$ or it could represent convolution by averaging over elements in the vector $\boldsymbol{\mathbf{r}}$. Also numerical differentiation or integration can be captured by $\boldsymbol{\mathbf{H}}$. In the case study we use a convolved, contrast-linearized approximation to the wave equation to model acquisition of seismic data, see Buland and Omre (2003).

2.3 Prior Model

The spatial variable of interest is represented in the vector $\boldsymbol{\mathbf{r}}$. We define the prior basis-pdf $f(\boldsymbol{\mathbf{r}})$ to be a spatially stationary Gaussian pdf,

[TABLE]

where the scalars $(\mu,\sigma^{2})$ are the stationary expectation and variance, respectively, while the spatial correlation $(n_{r}\times n_{r})$-matrix $\boldsymbol{\mathbf{C}}$ is defined by the spatial translation invariant correlation function $\rho(\boldsymbol{\mathbf{\tau}});\boldsymbol{\mathbf{\tau}}\in\mathcal{R}^{m}$ . This pdf is spatially stationary in the sense that the marginal pdfs $f(r_{i})=\phi_{1}(r_{i};\mu,\sigma^{2});i=1,\dots,n_{r}$ are all identical. Moreover, it exhibits ergodisity since $f(r_{i},r_{j})\rightarrow f(r_{i})f(r_{j})$ as $|\boldsymbol{\mathbf{x}}_{i}-\boldsymbol{\mathbf{x}}_{j}|\rightarrow\infty$, which entails that consistent estimates of the model parameters can be obtained.

We use this Gaussian pdf as basis-pdf to define a selection pdf, with the auxiliary $n_{r}$-vector $\boldsymbol{\mathbf{\nu}}$ extension,

[TABLE]

where $\gamma\in[-1,1]\subset\mathcal{R}$ is a coupling parameter while $\boldsymbol{\mathbf{e}}_{\nu|r}$ is a centred Gaussian $n_{r}$-vector with independent elements with variance $[1-\gamma^{2}]$, independent of $\boldsymbol{\mathbf{r}}$. The extended variable becomes jointly Gaussian,

[TABLE]

and consists of two grid-discretized spatial variables with variances $\sigma^{2}$ and $1$ respectively, with inter-correlation $\gamma\sigma$ and with identical spatial correlation function $\rho(\boldsymbol{\mathbf{\tau}})$.

Define the selection set $\mathcal{A}=\cup_{i=1}^{n_{r}}\mathcal{A}_{i}\subset\mathcal{R}^{n_{r}}$ with $\mathcal{A}_{i}=\mathcal{A}_{j};i,j=1,\dots,n_{r}$ , hence identical selection sets for each component in $\boldsymbol{\mathbf{r}}$. The corresponding spatial selection Gaussian pdf, which belong to the class of selection Gaussian pdfs, see Expression 9, is defined as,

[TABLE]

with the actual model parameters $\boldsymbol{\mathbf{\theta}}_{p}=(\mu,\sigma^{2},\gamma,\rho(\boldsymbol{\mathbf{\tau}}),\mathcal{A}_{i})$.

Note in particular that the prior model for the spatial variable of interest, represented by the selection Gaussian pdf on the vector $\boldsymbol{\mathbf{r}}$, is subject to the grid $\mathcal{L}_{\mathcal{D}}$. The selection Gaussian consept breaks down when the grid size tends to zero by infilling of the grid, see Minozzo and Ferracuti (2012) and Rimstad and Omre (2014). This lack of generality limits the use of the model, of course, but in many applications like image analysis, remote sensing and geophysics, there is a natural choice of grid $\mathcal{L}_{\mathcal{D}}$ due to the observation acquisition procedure. Moreover, this limitation is sheared by the categorical Markov random field model, see Besag (1974) and Kaiser et al. (2002), which has proven immensely useful in many applications.

The selection Gaussian class of pdfs can be shown to be closed under marginalization, see Arellano-Valle et al. (2006), and the uni-variate marginal pdf is,

[TABLE]

which is a selection Gaussian pdf with $\boldsymbol{\mathbf{c}}_{i}$ being the $i$’th column $n_{r}$-vector of the correlation matrix $\boldsymbol{\mathbf{C}}$. Note, however, that the two first moments of the marginal pdf, $\operatorname{E}(r_{Ai})$ and $\operatorname{Var}(r_{Ai})$ do not have nice closed form expressions. All lower dimensional marginal pdfs will also be selection Gaussian pdfs.

All marginal pdfs $f(r_{Ai});i=1,\dots,n_{r}$ are dependent on all $n_{r}$ elements of the auxiliary variable $\boldsymbol{\mathbf{\nu}}$, not only $\nu_{i}$ . This dependence causes the selection Gaussian pdf to be defined subjective to the grid $\mathcal{L}_{\mathcal{D}}$. This coupling is contrary to Gaussian marginal pdfs where all dependence of other dimensions are integrated out. The dependence will decline with distance, however, since the spatial correlation function $\rho(\boldsymbol{\mathbf{\tau}})$ will tend towards zero with increasing $|\boldsymbol{\mathbf{\tau}}|$. Consequently all marginal pdfs will be equal, due to symmetry of the regular grid $\mathcal{L}_{\mathcal{D}}$, except for edge effects of the grid. The decline of the spatial correlation with distance will also make the bi-variate $f(r_{Ai},r_{Aj})$ tend towards $f(r_{Ai})f(r_{Aj})$, hence independence, as the distance $|\boldsymbol{\mathbf{\tau}}_{ij}|=|\boldsymbol{\mathbf{x}}_{i}-\boldsymbol{\mathbf{x}}_{j}|$ increase. Consequently, the selection Gaussian pdf exhibit approximate stationarity and ergodicity in the sense defined above. These characteristics make the selection Gaussian pdf suitable as prior pdf in Bayesian spatial inversion.

The prior selection Gaussian pdf will naturally be inspected by simulation. Simulation is performed sequentially by first generating a realization of the auxiliary variable $\boldsymbol{\mathbf{\nu}}^{s}\in\mathcal{A}$ and thereafter generating the realization $\boldsymbol{\mathbf{r}}_{A}^{s}=[\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{\nu}}^{s}]^{s}$. Since the joint variable $[\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{\nu}}]$ is Gaussian, many efficient algorithms are available, see Supplement A.

The variety of the prior selection Gaussian model is exhibited in Figure 1 which is based on the model parameters listed in Table 1 with an anisotropic second-order exponential correlation function with anisotropy factor $(d_{h},d_{v})$. A more detailed discussion of the example is given in Supplement C. The spatial variable is represented on a $(64\times 64)$-grid. The anisotropy factors vary, and so do the correlation and the selection sets for the auxiliary variable. We observe a large variety of prior spatial models, all of them approximately stationary and ergodic in the sense discussed above. The prior models have marginal distributions that can be multi-modal, skewed or heavy-tailed, or a combination of these features. The computer demand for generating one such realization is typically a couple of minutes on a regular laptop computer.

2.4 Posterior Model

The posterior pdf is uniquely defined by the likelihood function and the prior pdf. With a likelihood function from the Gauss-linear class and a prior pdf from the selection Gaussian class, the posterior pdf will also, due to Theorem 1, be from the selection Gaussian class. The selection Gaussian class of prior pdfs is conjugate with respect to Gauss-linear likelihood functions. Hence the model parameters of the posterior pdf are analytically tractable based on the model parameters of the likelihood and prior models and the actual observations. The joint pdf is,

[TABLE]

and one may demonstrate that $[\boldsymbol{\mathbf{\nu}},\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{r}}]$ are conditionally independent. Note also that the joint $[\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{\nu}}\in\mathcal{A}]$ will be selection Gaussian, and so will the two marginals $[\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{\nu}}\in\mathcal{A}]$ and $[\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{\nu}}\in\mathcal{A}]$. Hence the marginal pdf of the observations will be dependent on the actual prior model, which the likelihood model will not. The focus of the study is on the posterior $[\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}},\boldsymbol{\mathbf{\nu}}\in\mathcal{A}]$ which will be selection Gaussian as well, see Appendix A.

From Theorem 1 and standard Gaussian theory one has,

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

This posterior pdf will of course be spatially non-stationary due to conditioning on the observations $\boldsymbol{\mathbf{d}}$. The pdf will, however, be in the class of selection Gaussian pdfs, see Expression 9, and hence be closed under marginalization and conditioning, with the corresponding model parameters analytically tractable. Assessment of the posterior pdf is usually made by simulation of realizations and locationwise prediction with associated precision intervals.

Simulation of realizations from the posterior pdf is made sequentially by first to generate a realization $[\boldsymbol{\mathbf{\nu}}^{s}\in\mathcal{A}|\boldsymbol{\mathbf{d}}]$ and thereafter to generate a realization $[\boldsymbol{\mathbf{r}}_{A}|\boldsymbol{\mathbf{d}}]^{s}=[\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{\nu}}^{s},\boldsymbol{\mathbf{d}}]^{s}$. Since the joint variable $[\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{\nu}}|\boldsymbol{\mathbf{d}}]$ is Gaussian, many efficient algorithms are available, and the algorithm actually used in the current study is specified in Supplement A with,

[TABLE]

where

[TABLE]

. Prediction of $[\boldsymbol{\mathbf{r}}_{A}|\boldsymbol{\mathbf{d}}]$ need to be carefully designed since we often define selection Gaussian prior models with multiple modes, and so will also the posterior pdf be. The traditional expectation (E) predictor based on a minimum locationwise squared error loss, denoted $\hat{\boldsymbol{\mathbf{r}}}_{E}=\operatorname{E}\{\boldsymbol{\mathbf{r}}_{A}|\boldsymbol{\mathbf{d}}\}$, will often appear in low-probability regions inbetween modes of the posterior pdf. The median (M) predictor based on a minimum locationwise absolute error criterion, denoted $\hat{\boldsymbol{\mathbf{r}}}_{M}=MED\{\boldsymbol{\mathbf{r}}_{A}|\boldsymbol{\mathbf{d}}\}$, shear the same tendency to appear in low-probability regions. The prefered predictor is the global maximum posterior predictor, but it is usually too computer demanding to determine since it requires optimization of a $n_{r}$-dimensional multi-modal function. Therefore we recommend the maximum posterior (MAP) predictor based on a maximum locationwise posterior criterion,

[TABLE]

with

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

which normally appear close to the dominant mode of the posterior pdf. This predictor is relatively simple to identify since the marginal posterior pdfs are known to be selection Gaussian pdfs with analytically assessable parameter values. The associated prediction $\alpha$-intervals will naturally be the interval between the upper/lower $\alpha/2$-quantiles of the marginal posterior pdfs, which usually must be assessed by simulation based inference. Note also, that these predictors and prediction intervals will correspond to kriging if the prior pdf is from the pure Gaussian class since MAP and E predictors coinside for uni-modal symmetrical pdfs.

The characteristics of the posterior selection Gaussian spatial model are exhibited in Figure 2 which is based on the parameter sets listed in Table 2. A more detailed discussion of the examples is presented in Supplement D. The spatial variable is represented on a $128$-grid with exact observations in grid nodes $16$ and $112$. The different prior models produces very different posterior realizations and predictions, all of them exactly honoring the observations of course. The MAP-predictor is particularly sensitive to multi-modal marginals in the prior model. The computer demand for this simple example is very modest since the posterior model is analytically tractable.

3 Model Parameter Inference

One challenge with this class of selection Gaussian pdfs is the lack of clear interpretation of the model parameters, even in the reduced parametrization used as spatial stationary prior pdf in this study. The fact that the model parameter values are dependent on the actual grid-design $\mathcal{L}_{\mathcal{D}}$ , complicates matters even more. In this section we discuss model parameter inference in some larger detail.

In order to perform Bayesian inversion, all model parameters of both the likelihood function and the prior pdf must be assigned values. The likelihood parameters $\boldsymbol{\mathbf{\theta}}_{l}$ are assumed to be known through studies of the observation acquisition procedure. The model parameters of the prior pdf $\boldsymbol{\mathbf{\theta}}_{p}$ are more complicated to elicit.

One may consider a hierarchical Bayesian inversion model, combining Bayesian inversion and Bayesian inference, and consider $\boldsymbol{\mathbf{\theta}}_{p}$ as a random variable with a suitable prior model $f(\boldsymbol{\mathbf{\theta}}_{p})$. Then, in principle, the posterior model for $\boldsymbol{\mathbf{\theta}}_{p}$ is available,

[TABLE]

Remember that $\boldsymbol{\mathbf{\theta}}_{p}=(\mu,\sigma^{2},\gamma,\rho(\boldsymbol{\mathbf{\tau}}),\mathcal{A}_{i})$ where $\mu\in\mathcal{R}$ , $\sigma^{2}\in\mathcal{R}_{+}$ and $\gamma\in[-1,1]\subset\mathcal{R}$ while $\rho(0)=1$ and $\rho(\boldsymbol{\mathbf{\tau}});\boldsymbol{\mathbf{\tau}}\in\mathcal{R}_{+}^{m}$ is a positive definite function, and $A_{i}\subset\mathcal{R}$. Hence, both assigning suitable prior models to $\boldsymbol{\mathbf{\theta}}_{p}$ and calculation of the normalizing constant appear as very complicated.

We recommend using training images of the spatial variable of interest and discretize them to a grid-design corresponding to $\mathcal{L}_{\mathcal{D}}$, ie. the same grid spacing along all dimensions. These discretized training images will in many applications be available since the grid-design often is defined by the observation acquisition procedure and hence used in several sites. Denote one such discretized training image by $n_{r}^{o}$-vector $\boldsymbol{\mathbf{r}}_{A}^{o}$. The corresponding Bayesian inference expression is,

[TABLE]

which also will be very complicated to assess for the full model parameter vector $\boldsymbol{\mathbf{\theta}}_{p}$. In Arellano-Valle et al. (2009) and Branco et al. (2013) this posterior model for the model parameter $\gamma$ given the other parameters in $\boldsymbol{\mathbf{\theta}}_{p}$ and with $\mathcal{A}_{i}=\mathcal{R}_{\oplus}$ is discussed. The authors of the former reference also provide guidelines for obtaining conjugate prior models for $\gamma$. The generalization of these results to cover the full prior model parameter vector $\boldsymbol{\mathbf{\theta}}_{p}$ appears as very complicated.

We choose to continue in a classical inference setting and develop the log-likelihood for observing $\boldsymbol{\mathbf{r}}_{A}^{o}$ as a function of $\boldsymbol{\mathbf{\theta}}_{p}$,

[TABLE]

In theory we may then define the maximum likelihood estimator for $\boldsymbol{\mathbf{\theta}}_{p}$ by,

[TABLE]

but in practice we need to parametrize also $\rho(\boldsymbol{\mathbf{\tau}})$ and $A_{i}$ in order to perform the optimization, and this parametrization will be problem specific. One major challenge in the optimization is that the probability $\Phi_{n_{r}^{o}}(\cdot)$ need to be recalculated for varying $\boldsymbol{\mathbf{\theta}}_{p}$, which may be extremely computer demanding. We calculate this probability by an importance blocking rejection algorithm, see Supplement B. Lastly, there is no guarantee that the object function $\log L(\boldsymbol{\mathbf{\theta}}_{p},\boldsymbol{\mathbf{r}}_{A}^{o})$ may not be multi-modal, which makes optimization notoriously complicated. It will be unfair to say that not many unresolved issues remain, but we present one encouraging example of prior model parameter elicitation based on training images below.

We evaluate the characteristics of the maximum likelihood estimator for $\boldsymbol{\mathbf{\theta}}_{p}=(\mu,\sigma^{2},d,\gamma,a)$, where $d=d_{h}=d_{v}$ and $\mathcal{A}_{i}:(-\infty,-a]\cup[a,\infty)$ , for case 1 in Figure 1 and Table 1. The results are exhibited in Figure 3, and a more detailed discussion is presented in Supplement E. The training images are subsets of the realization in Figure 1 of sizes $[8\times 8],[16\times 16],[24\times 24]$ and $[32\times 32]$. By repeating this inference on $1000$ realizations from the prior model we can assess the accuracy and precision of the estimator. We observe from Figure 3 that the estimator appears as biased, but consistent as the training image increases, which is as expected for maximum likelihood estimators for ergodic spatial models. Moreover, it appears as relatively reliable estimates can be obtained even for training images of size $[24\times 24]$. The computer demand for estimating $\boldsymbol{\mathbf{\theta}}_{p}$ for one training image of size $[32\times 32]$ is typically one minute on a regular laptop computer.

4 Case Study - Seismic inversion

The objective of seismic inversion is to predict the elastic material properties - pressure-wave velocity, shear-wave velocity and density - in the subsurface based on observed amplitude-versus-offset (AVO) seismic data collected at the surface. The data appear as time-laged, angle-dependent reflection intensities from the subsurface created by an air pulse generated at the surface. We model the log-transformed properties in order to have a linear likelihood function, see Buland and Omre (2003), $\boldsymbol{\mathbf{r}}=(\log\boldsymbol{\mathbf{v}}_{p},\log\boldsymbol{\mathbf{v}}_{s},\log\boldsymbol{\mathbf{\rho}})\in\mathcal{R}^{3n_{r}}$ and $\boldsymbol{\mathbf{d}}=(\boldsymbol{\mathbf{d}}^{1},\boldsymbol{\mathbf{d}}^{2},\boldsymbol{\mathbf{d}}^{3})\in\mathcal{R}^{3n_{d}}$ with upper-index representing three angles. Hence the objective is to assess $[\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}}]$, and we phrase the inversion in a Bayesian setting.

The case study is based on data from the Alvheim field in the North Sea ( Avseth et al. (2008) , Rimstad et al. (2012)). The subsurface contains a turbiditic oil and gas reservoir at about 2000 meters depth, but we use reflection time as depth reference with one meter (m) corresponding to approximately one milli-second (ms). We consider one vertical profile at the depth range $\mathcal{D}:[1935-2145]$ ms discretized to $\mathcal{L}_{\mathcal{D}}$ with $n=n_{r}=n_{d}=55$ grid nodes, where both AVO seismic data $\boldsymbol{\mathbf{d}}$ and exact observations of the elastic material properties $\boldsymbol{\mathbf{r}}^{o}$ in a well, are available, see Figure 4(a) and 5. Both $\boldsymbol{\mathbf{d}}$ and $\boldsymbol{\mathbf{r}}^{o}$ are used to infer the likelihood and prior model parameters $\boldsymbol{\mathbf{\theta}}_{l}$ and $\boldsymbol{\mathbf{\theta}}_{p}$, by considering $[\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{r}}^{o}]$ and $\boldsymbol{\mathbf{r}}^{o}$ as training images, respectively. In the Bayesian spatial inversion we consider $[\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}}]$, hence the posterior model is only conditioned on $\boldsymbol{\mathbf{d}}$. The training image $\boldsymbol{\mathbf{r}}^{o}$ is used to validate the results. In practical use, seismic inversion of profiles in the neighborhood of the well trace, without well observations, will be made, but then model validation is complicated. We perform Bayesian inversion based on two alternative prior models, one selection Gaussian and one traditional Gaussian prior model, and we compare the corresponding posterior models.

The likelihood model $f(\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{r}})$ link the seismic data $\boldsymbol{\mathbf{d}}$ and the elastic material properties of interest $\boldsymbol{\mathbf{r}}$. The model is based on a linearization of the wave equation as defined in Buland and Omre (2003),

[TABLE]

where $\boldsymbol{\mathbf{W}}$ is a convolution matrix defined by the kernels in Figure 4(b); matrix $\boldsymbol{\mathbf{A}}$ represents the angle-dependent linearized wave equation; $\boldsymbol{\mathbf{D}}$ is a diffentiation matrix; and $\boldsymbol{\mathbf{\epsilon}}_{d|r}$ is a centred Gaussian vector with covariance matrix $\boldsymbol{\mathbf{\Sigma}}_{d|r}=\sigma^{2}_{d|r}\boldsymbol{\mathbf{\Sigma}}_{d|r}^{o}$. The correlation matrix $\boldsymbol{\mathbf{\Sigma}}_{d|r}^{o}$ is defined by exponential correlation functions in angle and time with ranges $d^{a}$ and $d^{t}$, respectively. Hence the likelihood model is Gauss-linear with parameters $\boldsymbol{\mathbf{\theta}}_{l}=(\sigma^{2}_{d|r},d^{a},d^{t})$.

The prior model $f(\boldsymbol{\mathbf{r}})$ represents the general characteristics of the elastic material properties of interest. Figure 5 contain a plot of the exact observations of the properties $\boldsymbol{\mathbf{r}}^{o}$ along the profile plus Gaussian quantile-quantile plots and histograms of residuals after removing the linear vertical trend. The bi-modality of the histograms of $\log v_{p}$ and $\log v_{s}$ are caused by vertically varying rock types in the subsurface. By using a selection Gaussian prior model, this bi-modality in the marginal pdfs can be captured. The model as defined in Section 2.3 must be extended to represent the tri-variate $\boldsymbol{\mathbf{r}}$, and the parametrization become $\boldsymbol{\mathbf{\theta}}_{p}=(\boldsymbol{\mathbf{\mu}},\boldsymbol{\mathbf{\Sigma}},\boldsymbol{\mathbf{\gamma}},d^{r},\boldsymbol{\mathbf{a}})$. The spatial exponential correlation function with range $d^{r}$ is common for all three variables and the selection sets are parametrized as $\mathcal{A}_{i}:\{(-\infty,a],[a,\infty)\}$ with specific $a$ values for each variable. The alternative traditional Gaussian prior model $f^{G}(\boldsymbol{\mathbf{r}})$ has parametrization $\boldsymbol{\mathbf{\theta}}^{G}_{p}=(\boldsymbol{\mathbf{\mu}},\boldsymbol{\mathbf{\Sigma}},d^{r})$ and will not capture the bi-modality in the variables.

We infer the likelihood parameters $\boldsymbol{\mathbf{\theta}}_{l}$ from the available seismic data $\boldsymbol{\mathbf{d}}$ and the exact observations of the elastic material properties $\boldsymbol{\mathbf{r}}^{o}$ by using a maximum likelihood criterium,

[TABLE]

Likewise we infer the model parameters for the two alternative prior models, $\boldsymbol{\mathbf{\theta}}_{p}$ and $\boldsymbol{\mathbf{\theta}}_{p}^{G}$ from $\boldsymbol{\mathbf{r}}^{o}$. We set the location parameter $\boldsymbol{\mathbf{\mu}}$, for both models, equal to the vertical linear trend for each of the three variables, and estimate the remaining parameters by a maximum likelihood criterium,

[TABLE]

The optimizations of the likelihood functions all appear to converge to unique optima with computer demands of a few minutes on a regular lap-top computer. The actual estimates for the likelihood model parameters are,

[TABLE]

while the estimates for the two alternative prior model parameters are,

[TABLE]

[TABLE]

The parameter estimates for the two alternative prior models appear as consistent with comparable range lengths and dependence structures between the three variables. In Figure 6 the marginal pdfs of the two prior models are displayed together with the histograms of $\boldsymbol{\mathbf{r}}^{o}$. The selection Gaussian prior model captures the bi-modality of the histograms, without overfitting to the available well observations.

Based on the Gauss-linear likelihood model $f(\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{r}})$ with parameter values $\hat{\boldsymbol{\mathbf{\theta}}}_{l}$ and the selection Gaussian prior model $f(\boldsymbol{\mathbf{r}})$ with parameter values $\hat{\boldsymbol{\mathbf{\theta}}}_{p}$, we use Bayesian spatial inversion to assess the posterior model $f(\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}})$ which also will be selection Gaussian. By using the alternative traditional Gaussian prior model with associated parameter values we obtain a Gaussian posterior model $f^{G}(\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}})$.

Realizations from the two alternative posterior models, $f(\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}})$ and $f^{G}(\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}})$ are displayed in Figure 7. The realizations from the former can be generated sequentially as outlined in Section 2.4. The realizations from the selection Gaussian posterior model appear with abrupt changes between two levels defined by the two modes in the prior model, hence the time-integrated histograms are bi-modal. The realizations from the Gaussian posterior model are smoother and the corresponding histograms are uni-modal.

Predictions of $[\boldsymbol{\mathbf{r}}|\boldsymbol{\mathbf{d}}]$ based on the two alternative posterior models with associated 0.8-prediction intervals are displayed in Figure 8. Also the correct elastic material property profiles $\boldsymbol{\mathbf{r}}^{o}$ are presented. Moreover, predictions and 0.8-prediction intervals for the two alternative prior models are displayed. We use E-predictors for both models. In this case study, contrary to the example in Section 2.4, we have densely sampled data. The data consists of convolved, gradient observations in every node of $\mathcal{L}_{\mathcal{D}}$, and the marginal pdf of the selection Gaussian posterior model will appear as almost uni-modal, although flipping between modes vertically. Consequently the E, MED and MAP predictors will be almost identical and the former is used for computational convenience. The predictions based on both posterior models do reproduce the correct profiles relatively well with large improvements of the prior predictions. The predictions from the selection Gaussian posterior model appear with more abrupt changes than the Gaussian one, whenever the correct profiles has large steps. Moreover, the 0.8-prediction intervals are narrower for the former model than for the latter.

Table 3 contains summary statistics for Figure 8. The mean square error (MSE) of the predictions relative to the correct profiles for each variable for both alternative models are listed. Moreover, the coverage of the 0.8-prediction intervals for the correct profiles for the prior and posterior models are specified. The predictions from the selection Gaussian model appear as clearly superior to the predictions from the Gaussian one. The improvements in MSE are in the range of $20-40$ %. The coverage values for the selection Gaussian model is close to $0.8$ as they should be, while the coverages for the Gaussian model are far too large and more variable.

5 Concluding remarks

We study Bayesian spatial inversion and introduce the consept of conjugate classes of prior parametrized pdfs with respect to given classes of likelihood functions. For this class of prior pdfs the associated posterior pdfs will be in the same class. Such conjugate classes exist for continuous, event and mosaic spatial variables for frequently used likelihood functions. The conjugate class of prior models can be selection extended without loss of the conjugate characteristic.

We demonstrate the potential of the selection extension by introducing the class of selection Gaussian prior pdfs which is conjugate with respect to Gauss-linear likelihood functions. The flexibility of this selection Gaussian class is displayed in a variety of examples which represent multi-modality, skewness and heavy-tailedness in the marginal distributions.

By using a prior model from a conjugate class for a given likelihood function, the associated posterior model can be assessed exactly based only on the model parameters of the prior and likelihood models - and the actual observations, of course. The normalizing constant, which usually complicates Bayesian spatial inversion, will be available on parametric form. We demonstrate this favorable characteristic for the class of selection Gaussian prior pdfs, and show that the posterior selection Gaussian pdf is analytically tractable which is used to make efficient algorithms for simulation and prediction. Several examples presenting conditional simulations and predictions exposing multi-modality, skewness and heavy-tailedness are displayed.

The class of selection Gaussian prior pdfs is parametrized by a number of model parameters which are not easily interpretable. Based on one training image of the spatial variable we define maximum likelihood estimators for the model parameters. A limited simulation study is conducted and we conclude that the estimators appear as consistent and that even for relatively small training images reliable estimates can be obtained. These results are encouraging.

Lastly, a case study using the selection Gaussian prior model on seismic inversion of real data is presented. We demonstrate $20$-$40\%$ improvement in the mean-square-error of predictions compared to the traditional Gaussian inversion. Also the prediction intervals of the former model appear as more reliable than for the latter.

The selection extension of conjugate classes of prior pdfs in Bayesian spatial inversion appears to have a large potential. We have to some extent explored this potential for continuous spatial variables and the class of selection Gaussian prior pdfs. The challenge for these models appears in sampling from and calculation of sub-set probabilities in high-dimensional Gaussian pdfs. We have presented some relatively efficient algorithms for these purposes. Many improvements of these algorithms are definitely possible. For event and mosaic spatial variables, the class of Poisson and Markov pdfs are conjugate with respect to certain likelihood functions. Also these classes can be selection extended and still remain conjugate. We have not yet explored these possibilities.

Acknowledgments

The research is a part of the Uncertainty in Reservoir Evaluation (URE) activity at the Norwegian University of Science and Technology (NTNU).

Appendix A Selection Gaussian Model

The closedness properties of the selection Gaussian model is demonstrated.

Definition 3** (Selection Gaussian pdf)**

Consider the $n$-vector Gaussian basis-pdf,

[TABLE]

and Gauss-linear auxiliary $q$-vector variable,

[TABLE]

with $\boldsymbol{\mathbf{\mu}}_{\nu|r}=\boldsymbol{\mathbf{\mu}}_{\nu}+\boldsymbol{\mathbf{\Gamma}}_{\nu|r}(\boldsymbol{\mathbf{r}}-\boldsymbol{\mathbf{\mu}}_{r})$, where $\boldsymbol{\mathbf{\Gamma}}_{\nu|r}$ is denoted the coupling $(q\times n)$-matrix.

Define a selection set $\mathcal{A}_{\nu}\subset\mathcal{R}^{q}$, and the corresponding $n$-vector selection Gaussian pdf,

[TABLE]

with the covariance $(q\times q)$-matrix $\boldsymbol{\mathbf{\Sigma}}_{\nu}=\boldsymbol{\mathbf{\Gamma}}_{\nu|r}\boldsymbol{\mathbf{\Sigma}}_{r}\boldsymbol{\mathbf{\Gamma}}_{\nu|r}^{T}+\boldsymbol{\mathbf{\Sigma}}_{\nu|r}$.

The class of selection Gaussian pdfs is defined by all valid sets of parameters $(\boldsymbol{\mathbf{\mu}}_{r},\boldsymbol{\mathbf{\Sigma}}_{r},\boldsymbol{\mathbf{\mu}}_{\nu},\boldsymbol{\mathbf{\Gamma}}_{\nu|r},\boldsymbol{\mathbf{\Sigma}}_{\nu|r},\mathcal{A}_{\nu})$.

The following results are useful for later Proofs,

Result 1** (Conditional Probabilities)**

Consider the joint $(n+m)$-vectorial variable $(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$ with joint pdf $f(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$, then,

[TABLE]

and also for arbitrary subset $\mathcal{A}_{x}\subset\Omega_{x}$,

[TABLE]

For $f(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$ being a Gaussian pdf,

[TABLE]

The major statements are captured in the following Proposition,

Proposition 1** (Selection Gaussian Models)**

Consider the selection Gaussian prior model,

[TABLE]

and Gauss-linear $m$-vector likelihood model,

[TABLE]

with conditional expectation $\boldsymbol{\mathbf{\mu}}_{d|r}=\boldsymbol{\mathbf{H}}\boldsymbol{\mathbf{r}}$ where $\boldsymbol{\mathbf{H}}$ is an observation $(m\times n)$-matrix. Moreover, assume conditional independence of $[\boldsymbol{\mathbf{\nu}},\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{r}}]$.

Then the following holds:

A.

$[\boldsymbol{\mathbf{r}}_{A},\boldsymbol{\mathbf{d}}]$* is selection Gaussian*

B.

$\boldsymbol{\mathbf{d}}$* is selection Gaussian*

C.

$[\boldsymbol{\mathbf{r}}_{A}|\boldsymbol{\mathbf{d}}]$* is selection Gaussian*

The Proposition is justified by the following Proof,

Proof 1

The three proposition items are demonstrated sequentially.

The joint pdf in A is,

[TABLE]

with the last identity from conditional independence of $[\boldsymbol{\mathbf{\nu}},\boldsymbol{\mathbf{d}}|\boldsymbol{\mathbf{r}}]$.

Hence from Definition 3 , the joint $(n+m)$-vector $[\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{d}}]$ is selection Gaussian.

The marginal pdf in B is,

[TABLE]

with the last identity from Result R2G.

Hence from Definition 3, the marginal $m$-vector $\boldsymbol{\mathbf{d}}$ is selection Gaussian.

The conditional pdf in C is,

[TABLE]

Hence from Definition 3, the conditional $n$-vector $[\boldsymbol{\mathbf{r}}_{A}|\boldsymbol{\mathbf{d}}]$ is selection Gaussian.

Supplementary Material

Appendix A Sampling - truncated multivariate Gaussian distribution

Consider the problem of sampling from a $n$-dimensional truncated multivariate normal distribution with unnormalized density $I(\mathbf{x}\in A)\times\phi_{n}(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\Sigma})$, where $\mathbf{x},\boldsymbol{\mu}\in\mathcal{R}^{n}$, $\boldsymbol{\Sigma}\in\mathcal{R}^{n\times n}$, $A=A_{1}\times\ldots\times A_{n}$, $A_{i}\subseteq\mathcal{R}$, $I(\cdot)$ is the indicator function, and $\phi_{n}(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\Sigma})$ is the multivariate Gaussian density distribution with expectation vector $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$. In order to sample from this distribution we extend the Metropolis-Hastings algorithm in Robert (1995) with a block independent proposal distribution:

[TABLE]

where $n_{a}$ is the block size, $\mathbf{x}^{a}\in\mathcal{R}^{n_{a}},\mathbf{x}^{b}\in\mathbb{R}^{n-n_{a}}$, $\phi_{1}(x_{i}^{a}\mid\mathbf{x}_{1:i-1}^{a},\mathbf{x}^{b};\boldsymbol{\mu},\boldsymbol{\Sigma})$ the conditional Gaussian probability of $x_{i}^{a}$ given $\mathbf{x}_{1:i-1}^{a}$ and $\mathbf{x}^{b}$, and $\Phi_{1}(x^{a}_{i}\in A_{i}\mid\mathbf{x}_{1:i-1}^{a},\mathbf{x}^{b};\boldsymbol{\mu},\boldsymbol{\Sigma})$ is the probability of the set $A_{i}$ under the Gaussian probability distribution of $x_{i}$ given $\mathbf{x}_{1:i-1}^{a}$ and $\mathbf{x}^{b}$. We use the notation $\mathbf{x}_{1:i-1}=(x_{1},x_{2},\ldots x_{i-1})$. The approach is inspired by the importance sampler in Genz (1992). Note that $p^{*}(\mathbf{x}^{a}\mid\mathbf{x}^{b})$ is normalized and it is easy to sample from the distribution due to the sequential structure.

The acceptance probability in the accept/reject step is

[TABLE]

where ${\mathbf{x}^{a}}^{\prime}$ is the new proposed state. The Metropolis-Hastings algorithm is presented in Algorithm 1.

In practice we calculate the conditional distributions in Algorithm 1 in advance. To save memory and time we also limit the elements in $\mathbf{x}$, i.e. sets eligible for choice, such that all elements in $\mathbf{x}$ has approximately equal update probability. We normally use the block size $n_{a}=100$.

Appendix B Estimation - multivariate Gaussian probabilities

Consider the problem of estimating the multivariate Gaussian probability

[TABLE]

where $\mathbf{x},\boldsymbol{\mu}\in\mathcal{R}^{n}$, $\boldsymbol{\Sigma}\in\mathcal{R}^{n\times n}$, $A=A_{1}\times\ldots\times A_{n}$, $A_{i}\subset\mathcal{R}$, $I(\cdot)$ is the indicator function, and $\phi_{n}(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\Sigma})$ is the multivariate Gaussianl density distribution with expectation vector $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$. The usual importance sampling Monte Carlo approximation is

[TABLE]

with $\mathbf{x}^{j}\sim f_{n}(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\Sigma});\;j=1,\ldots N$ and $N$ is the number of Monte Carlo sampling points. We extend the approach presented in Genz (1992) by allowing $A_{i}$ to consist of several intervals, and use

[TABLE]

as importance function, where $\phi_{1}(x_{i}\mid\mathbf{x}_{1:i-1};\boldsymbol{\mu},\boldsymbol{\Sigma})$ the conditional Gaussian probability of $x_{i}$ given $\mathbf{x}_{1:i-1}$, and $\Phi_{1}(A_{i}\mid\mathbf{x}_{1:i-1};\boldsymbol{\mu},\boldsymbol{\Sigma})$ is the probability of the set $A_{i}$ under the Gaussian probability distribution of $x_{i}$ given $\mathbf{x}_{1:i-1}$. We use the notation $\mathbf{x}_{1:i-1}=(x_{1},x_{2},\ldots x_{i-1})$. However, we also introduce a mean shift parameter $\boldsymbol{\eta}$ in the importance function which is important for asymmetric sets $A_{i}$. Then the importance sampling approximation appear as

[TABLE]

with $\mathbf{x}^{j}\sim f_{n}(\mathbf{x};\boldsymbol{\mu}+\boldsymbol{\eta},\boldsymbol{\Sigma}),\;j=1,\ldots N$.

Appendix C Example: Prior Model

The flexibility of the selection Gaussian pdf as prior model for the spatial variable of interest is demonstrated by generating realizations with varying model parameter sets, $\boldsymbol{\mathbf{\theta}}_{p}=(\mu,\sigma^{2},\gamma,\rho(\boldsymbol{\mathbf{\tau}}),\mathcal{A}_{i})$. The spatial variable is represented on a $(64\times 64)$-grid $\mathcal{L}_{\mathcal{D}}$ covering $\mathcal{D}\subset\mathcal{R}^{2}$, hence a surface in two dimensions. The spatial correlation function is parametrized as $\rho(\boldsymbol{\mathbf{\tau}};(d_{h},d_{v}))=\exp\{-[[\tau_{h}^{2}/d_{h}^{2}]+[\tau_{v}^{2}/d_{v}^{2}]]\}$, hence to be a second-order exponential correlation function with anisotropy factors $(d_{h},d_{v})$. The selection set $\mathcal{A}_{i}\subset\mathcal{R}$ is parametrized as a number of line segments on $\mathcal{R}$. The example design is summarized in Table 4.

In order to simulate realizations from the selection Gaussian pdf we extend the Metropolis Hastings (MH) algorithm presented in Rimstad and Omre (2014) by allowing more general selection sets $\mathcal{A}$. The algorithm is summarized in Supplement A. The algorithm is a block proposal MH-algorithm, and we normally use block sizes about $100$ which in our examples give an acceptance rate of about $0.25$. The computer demand for generating one realization is a couple of minutes on a regular laptop computer. The burn-in and mixing appear as satisfactory and are not displayed.

Figure 9 displays the results from the six cases. The first column displays the marginal distribution in the $(32,32)$ location of the grid compared with a univariate Gaussian distribution with the same two first moments as the marginal selection Gaussian pdf. The selection set $\mathcal{A}$ is illustrated with a thick gray line at the bottom of the display. The selection set is comparable to the marginal of the selection Gaussian pdf because $\sigma^{2}=1$. The second column displays Gaussian quantile-quantile plots of the marginal distributions. The last column displays realizations from the selection Gaussian prior pdfs.

The first row in Figure 9, case 1, displays a symmetric bimodal spatially isotropic model. The selection region for the auxiliary variable is absolute values greater than $0.3$. The marginal pdf is symmetric and bimodal, and the quantile-quantile plot shows clear deviations from the Gaussian distribution. In the realization the two modes are visible as two separated levels with sharp transitions between them.

Case 2 is displayed in the second row in Figure 9 and this model is also symmetric and bimodal, but spatially anisotropic. The selection region for the auxiliary variable is absolute values greater than $0.3$, as in case 1. In this case the horizontal spatial correlation is increased and the vertical one decreased, while the coupling parameter is reduced. The resulting realization is clearly layered, with marginal pdf very similar to the one in case 1.

The third row in Figure 9, case 3, displays an asymmetric bimodal spatially anisotropic model. The selection is asymmetric and further out in the tails than for previous cases. The occurrence of two clearly separated modes is made possible by low spatial correlation which allows larger jumps. The asymmetric selection causes the mode to the left to be smaller than the mode to the right. In spite of the low correlation the realization appears with clear spatial anisotropy and with two distinct modes.

Case 4 is displayed in the forth row in Figure 9 and represents a symmetric trimodal spatial isotropic model. The selection contains three symmetric closed intervals, which provides a trimodal symmetric marginal pdf. The three modes are distinctly separated, and clearly visible in the realizations. The spatial transitions between the two outer modes seem always to pass through the middle mode.

The fifth row in Figure 9, case 5, displays a skewed model. The skewed models considered in Allard and Naveau (2007) and Rimstad and Omre (2014) only contain one-sided selection interval, hence the model formulation constrains the degree of skewness. In the current case one additional selection interval is introduced, which provides a more flexible skewness structure for the model. The skewness is evident in the marginal pdf, in the quantile-quantile plot, and in the realization of the selection Gaussian prior pdf.

The last row in Figure 9, case 6, displays a symmetric peaked heavy tailed model. Symmetric selection is used, which forces higher probability density in the centre and in the tails. The extreme tails still decay exponentially, as seen in the quantile-quantile plot, while the more visible effects is due to the heavy tails. The closest univariate Student-$t$ distribution, if we ignore the extreme tails, has about $2$ degrees of freedom.

Appendix D Example: Posterior Model

The posterior model and the effect of using different prediction criteria is demonstrated on a small example. The spatial variable is defined on a 128-grid $\mathcal{L}_{\mathcal{D}}$ covering $\mathcal{D}\subset\mathcal{R}$, hence a one-dimensional case. The prior from the selection Gaussian class is defined with the spatial correlation function $\rho(\tau)=\exp\{-\tau^{2}\}$, and a selection set $\mathcal{A}\subset\mathcal{R}$ consisting of several line segments. The example design is summarized in Table 5. The four cases shear about the same characteristics as case $1,4,5$ and $6$ in the example in Supplement C. The likelihood function provides exact observations at grid locations 16 and 112. The actual observed values are also specified in Table 5.

The posterior model is from the selection Gaussian class, and the exact observations are of course exactly reproduced. The corresponding conditional selection Gaussian pdf must be calculated and realizations from this pdf can be generated by the algorithm in Supplement A.

We compare Bayesian inversion based on selection Gaussian and Gaussian prior models, and the effect of using different predictors. The two first moments are set identical in the two alternative prior models for each case. The results are summarized in Figure 10. The first column in the figure contains the marginal pdfs at location $64$ for both the selection Gaussian and Gaussian prior models for the four cases. Selection sets for the former are plotted as grey line segments. The second column contain realizations from the posterior pdf of the selection Gaussian model given observations at locations $16$ and $112$. Lastly, the third column contains predictions based on different criteria for each of the two posterior models for each case. The E (black solid), MED (black hatched) and MAP (black hatch-dot) predictors for the selection Gaussian model and the E/MED/MAP (grey hatched) predictors, which coinsides, for the Gaussian model. We focus our discussion on the results from the selection Gaussian model, and compare them with the Gaussian predictor.

The first row, case 1, displays a symmetric bi-modal prior model. The marginal distribution in Figure 10 is clearly bimodal. We condition on the values $2.5$ and $-2.5$ at grid nodes $16$ and $112$, respectively. The realizations of the posterior model have a evident bi-modal structure. The E predictor is almost identical to the Gaussian predictor, while the MED and MAP predictors clearly deviate from the Gaussian predictor. The MAP predictor has a stepwise structure reflecting the bimodality of the posterior model and the MED predictor is somewhere between the MAP and E predictors, but closest to the E predictor.

Case 2, displayed in the second row, has a symmetric tri-modal prior model. We condition on the values $0.55$ and $-0.55$ at grid nodes $16$ and $112$, respectively. The three modes are clearly visible in both the marginal distribution of the prior and the posterior realizations. The E predictor is in this case also almost identical to the Gaussian predictor. The MAP predictor has a stepwise structure with three levels, and the MED predictor is in this case closest to the MAP predictor.

The third row, case 3, displays an asymmetric unimodal prior model. We condition on the values $1.0$ and $-3.0$ at grid nodes $16$ and $112$, respectively. The marginal distribution of the prior is obviously skewed and the posterior realizations have a skewed structure. The E predictor and the Gaussian predictor are again almost identical. The MAP and MED predictors are similar to the E predictor except that the stationary levels for the MAP and MED are somewhat shifted relative to the E predictor.

Case 4, displayed in the last row, has a symmetric heavy tailed prior model. We condition on the values $3.0$ and $-3.0$ at grid nodes $16$ and $112$, respectively. All the predictions have similar shapes, but the MAP predictor, followed by the MED predictor, decays faster toward the stationary level than the E predictor. Again the E predictor and the Gaussian predictor are almost identical. The fact that the E, MED and MAP predictors are not identical entails that the posterior distributions are asymmetric, in spite the prior distribution being unimodal and symmetric.

The E, MED and MAP predictors can be very different for selection Gaussian models, contrary to the Gaussian model where all the three predictors are identical. The predictors are particularly different for multi-modal prior models, where the MAP predictor appears as stepwise. The E predictors for the selection Gaussian model is almost identical to the predictor for the corresponding Gaussian prior model.

Appendix E Example: Model Parameter Inference

We consider case 1 in the example in Supplement C, which has a selection Gaussian prior model with symmetric, bi-modal marginal pdfs and isotropic spatial correlation. The prior model is parametrized by $\boldsymbol{\mathbf{\theta}}_{p}=(\mu,\sigma^{2},d,\gamma,a)$ , where $(\mu,\sigma^{2},d)$ defines the stationary, isotropic Gaussian basis-pdf with the two former parameters being expectation and variance respectively, while the latter is the range in an isotropic second-order exponential spatial correlation function. The parameter $\gamma$ is the coupling parameter and $a$ defines the selection set $\mathcal{A}_{i}:(-\infty,-a]\cup[a,\infty)\subset\mathcal{R}$.

Consider one training image $\boldsymbol{\mathbf{r}}_{A}^{o}$ discretized to grid $\mathcal{L}_{\mathcal{D}}$. We estimate the prior model parameters $\boldsymbol{\mathbf{\theta}}_{p}$ by a maximum likelihood approach. The challenging calculations are repeated assessment of $\log\Phi_{p}(\mathcal{A};\mathbf{0},(1-\gamma^{2})\mathbf{I}_{p}+\gamma^{2}\mathbf{C})$ for varying values of $\boldsymbol{\mathbf{\theta}}_{p}$, which we solve by Monte Carlo importance sampling approach inspired by Genz (1992) and Genz and Bretz (2009), see Supplement B. For a grid $\mathcal{L}_{\mathcal{D}}:[32\times 32]$ the computation of the likelihood function for one set of $\boldsymbol{\mathbf{\theta}}_{p}$ requires about one minute on a regular laptop computer. In order to ensure a smooth likelihood function which can be optimized by standard procedures we keep the Monte Carlo samples fixed during each optimization. We used $N=5000$ samples which according to results in Rimstad and Omre (2014) should ensure stable solutions. There may exist multiple local optima in the likelihood function, hence we initiate the optimization in multiple points, but we encountered few problems with multiple solutions.

We used the following experimental design: generate $1000$ realizations of the model with parameters $\boldsymbol{\mathbf{\theta}}_{p}=(0,1,2,0.8,0.3)$ on a grid $\mathcal{L}_{\mathcal{D}}:[32\times 32]$. For each realization select one subset on grid sizes: $[8\times 8],[16\times 16],[24\times 24]$ and $[32\times 32]$, and estimate the model parameter $\boldsymbol{\mathbf{\theta}}_{p}$ as previously described.

The results are summarized in Figure 11 and 12. From the Figure 11 we observe that the estimator for $\boldsymbol{\mathbf{\theta}}_{p}$ is relatively well centred at the correct values even for small grid sizes. The centering improves with increasing grid sizes, while the estimation variance decreases. Hence it appears as the estimator is consistent with increasing grid sizes, but not in general unbiased. These results are as expected for maximum likelihood estimators. In Figure 12 we display multivariate results for the grid size $[16\times 16]$. The negative correlation between the estimators for $\gamma$ and $a$, the coupling and selection set, is easy to understand since stronger coupling requires selections closer to zero. Moreover, there appears to be some positive correlation between estimators for $\sigma^{2}$ and $d$, variance and range.

It is encouraging that relatively stable parameter estimates can be obtained from small training images of grid size about $[24\times 24]$. These results make us believe that the parametrization of the selection Gaussian prior model is reasonable and that the model parameters can be robustly assessed from training images of reasonable size. Moreover, it makes us trust the numerical approximations used in the optimization procedure. Cases 2 through 6 in the example in Supplement C, require parametrizations with more parameters, which may complicate the evaluations of the likelihood functions and introduce ambiguities among parameters. These complications are not further considered in the current study.

We have based our inference study on a training image of exact, complete set of observations on a grid. Alternatively, maximum likelihood inference can be made from any set of observations from a Gauss-linear likelihood model since the posterior model will also be a selection Gaussian pdf.