A parameterization of anisotropic Gaussian fields with penalized complexity priors
Liam Llamazares-Elias Jonas Latz Finn Lindgren


apalike

#1 1





1
        
        1]L. Llamazares-EliasCorresponding author: L.S.Llamazares-Elias@sms.ed.ac.uk
        [1]School of Mathematics, University of Edinburgh, United Kingdom
        2]J. Latz
        [2]Department of Mathematics, University of Manchester, United Kingdom
        1]F. Lindgren
        A parameterization of anisotropic Gaussian fields with penalized complexity priors
    [
    
==================================================================================

 

0
        

            Anisotropic Gaussian random fields with identifiable parameters and penalized-complexity priors
        

         







§ ABSTRACT 
    Gaussian random fields (GFs) are fundamental tools in spatial modeling and can be represented flexibly and efficiently as solutions to stochastic partial differential equations (SPDEs). The SPDEs depend on specific parameters, which enforce various field behaviors and can be estimated using Bayesian inference. However, the likelihood typically only provides limited insights into the covariance structure under in-fill asymptotics. In response, it is essential to leverage priors to achieve appropriate, meaningful covariance structures in the posterior. This study introduces a smooth, invertible parameterization of the correlation length and diffusion matrix of an anisotropic GF and constructs penalized complexity (PC) priors for the model when the parameters are constant in space. The formulated prior is weakly informative, effectively penalizing complexity by pushing the correlation range toward infinity and the anisotropy to zero.



Keywords: Anisotropy, Bayesian, Penalized Complexity, Prior, Stochastic Partial Differential Equation.


1.9 


§ INTRODUCTION


Gaussian random fields (GF) are widely used to model spatial phenomena <cit.> while accounting for the uncertainty that may arise due to measurement error, model misspecification, or incomplete information. The prevalence of GFs is due to the fact that they are well understood theoretically, verify desirable properties, and are easily characterized – they are entirely specified by their mean and covariance <cit.>.

A convenient way of representing certain GFs is as solutions to stochastic partial differential equations (SPDEs). This representation allows for physical interpretation to be assigned to the parameters of the equation. Furthermore, it allows for computationally efficient inference, prediction, and uncertainty quantification using a finite element (FEM) approximation of the field <cit.>.

In the literature, a common choice is to model using stationary fields. That is, the correlation of the field at two locations only depends on the Euclidean distance between said locations. While this may be an appropriate assumption in some cases, in others, it is inadvisable. This limitation can be overcome by introducing additional parameters to model the anisotropy present in the field as in  <cit.>. In the following, we consider the semi-parametric estimation of the random field and its anisotropy parameters. The existing work leaves us with two significant challenges:


    
  * The anisotropy parameterization from <cit.> is non-identifiable as it has two parameter combinations for each anisotropy matrix, leading to a multi-modal likelihood, making it unsuitable as a general parameterization.
    
  * Given that not all parameters can be recovered under in-fill asymptotic <cit.>,  the choice of prior distribution on the parameters of the model may significantly impact the posterior distribution. As a result, suitable priors need to be defined.

To address these issues, we make the following contributions:

    
  * Alternative parameterization: We present an alternative parameterization of the anisotropy that preserves parameter interpretability and which is smooth and invertible.
    
  * Prior definition: We construct penalized complexity (PC) priors <cit.>  for the parameters in the model. An additional benefit of this construction is that it avoids overfitting by favoring simpler base models.
    
  * Validation and prediction: We conduct a simulation study that shows that PC priors outperform “non-informative” priors. We then use the derived model to study precipitation in Norway and show that the anisotropic model outperforms the isotropic one in the presence of limited information.


The outline of the work is as follows. In <Ref>, we introduce and motivate our anisotropic model. In <Ref>, we address Contribution <ref> and also construct a transformation that renders the parameters a Gaussian vector, providing the convenience of working with Gaussian random variables. In <Ref>, we present Contribution <ref>. Next, in <Ref>, we conduct a simulation study of the designed PC priors and compare the results to other possible priors on the parameters. In <Ref>, we study the performance of the model and priors on a real data set. Finally, in <Ref>, we synthesize the obtained results and discuss future avenues of research.




§ ANISOTROPIC MODEL 



 §.§ Formulation

In this section, we address the first main contribution. How can we introduce anisotropy in our model and parameterize it such that each parameter has an interpretable effect? We work in 2 dimensions and within the framework of the SPDE approach <cit.>. In this approach, the GF u is defined as the solution to a SPDE. Most commonly, the SPDE chosen is

    (κ^2-Δ)^α/2u=.

Here,  is Gaussian white noise on L^2(^2), and is defined such that  given f,g ∈ L^2(^2) the random measure

    (f):=⟨ f,⟩

verifies

    [(f)(g)]=f,g_L^2(^2).

The resulting field u has mean zero and is chosen to be isotropic to guarantee uniqueness. We recall that, by definition, the field u is isotropic  if there exists some function r such that, for all x,y∈^2,

    K(x,y):=Cov[u(x),u(y)]=r(y-x).

The isotropy of u may be an appropriate assumption in some cases. However, it is inadvisable if the correlation of the field is not equal in all spatial directions. To model this anisotropy, we will consider the model given by

    (κ^2-∇·H∇)u/σ_u =√(4 π)κ,

where the parameters are κ, σ_u ∈ (0,∞) which are positive and bounded away from zero, and a symmetric positive definite matrix  H∈^2× 2 with eigenvalues bounded away from zero. The parameters control the length scale, marginal variance, and anisotropy, respectively, as we will explain in the next subsections.

The formulation in (<ref>) preserves the advantages of the SPDE approach. Namely, representing u as the solution to a SPDE gives it a physical interpretation. The term κ^2 u represents reaction whereas -∇H·∇ u represents diffusion <cit.>.

Furthermore, using a finite element method, u can be projected onto the finite-dimensional Hilbert space _n spanned by the basis functions ψ_1, …, ψ_n linked to a mesh M_n of the domain. This gives a sequence of Gaussian Markov random fields u_n with sparse precision matrices (the precision matrix is the inverse of the covariance matrix), which converges in distribution to u as the mesh becomes finer and finer. This sparsity allows for a significant speed-up in computations (Kriging, posterior simulation, likelihood evaluation, etc.) <cit.>.


 §.§ Derivation

We now motivate the choice of the SPDE (<ref>). By definition, a field u on ^d  is stationary if its covariance K depends only on the relative position between two points. That is, for some function r called the covariance function

    K(x,y):=Cov[u(x),u(y)]=r(y-x),   ∀x,y∈^d.

Equation (<ref>), which is an extension of (<ref>), requires that the Euclidean translation y-x captures all relevant information related to the covariance of the field between any two spatial locations x and y. However, Euclidean geometry may not fit with the underlying properties of the field. For example, suppose that our field describes the geological properties of some homogeneous terrain. Even if the field was initially stationary, stationarity is lost if the terrain underwent a geological deformation ψ. Instead, the relationship

    K(x,y)=r(ψ^-1(x)-ψ^-1(y))

would be more appropriate. This
is the deformation method <cit.> and is, for instance, sometimes used to connect different layers in deep Gaussian processes. In deep Gaussian processes, this transformation is then random and not diffeomorphic <cit.>.


The deformation method provides an easy way to construct non-stationary fields starting from stationary ones. Furthermore, it can be used in a layered approach, where a stationary field is first built through some appropriate method and then deformed into a non-stationary field.

As discussed previously, SPDEs provide a convenient framework for the construction of stationary fields, for example, through (<ref>). Let us see what happens when we deform such a field. To this aim, consider an (unknown) d-dimensional manifold 𝒟 and our target manifold 𝒟 obtained through a diffeomorphic transformation

    ψ:𝒟→𝒟.

Further consider the solutions u to 
    (1-Δ)u=√(4π),  on 𝒟.

Here, either =^d, in which case, to obtain uniqueness, we impose a stationarity condition. Otherwise,  is some bounded subset of ℝ^d, in which case Neumann conditions or Dirichlet conditions are imposed on the boundary.
Then, a change of variables yields that u:= σ_u ·u∘ψ^-1 verifies the non-stationary SPDE

    1/Ψ[1-Ψ∇·ΨΨ^T/Ψ∇] u/σ_u=√(4π)/Ψ^1 / 2𝒲̇  on ,

where  Ψ is the Jacobian of  ψ, and we denote its determinant by Ψ <cit.>.

Let us write γ for the geometric mean of the eigenvalues of Ψ. Then, Ψ=γ^d and we obtain from (<ref>) that

    1/γ^d[1-γ^d ∇·γ^2-dH∇] u(x)/σ_u=√(4π)/γ^d / 2  on ,

where we defined H:=γ^-2ΨΨ^T. Note that, if we write κ:=γ^-1 and take the dimension d to be 2 we recover our model in (<ref>). For general d, it follows from the definition of H that

    
  * H is symmetric.
    
  * H has determinant 1.
    
  * H is positive definite.

From now on, we will impose these three assumptions on H. Furthermore, we will restrict ourselves to the stationary case by imposing that H is constant in space, or equivalently, we impose that  ψ is a linear deformation. That is, ψ(x)=Ψx.
By construction, the solution field u is thus geometrically anisotropic <cit.>. That is, if we replace the Euclidean metric with the deformed metric x_Ψ^-1:= Ψ^-1x, we obtain analogously to (<ref>) that, for some function r:^d → called the  covariance function of u,

    K(x,y):=Cov[u(x),u(y)] =r(y-x_Ψ^-1),   ∀x,y∈^d,

In the case =^d, the marginal variance of u is 𝔼[u (x)^2]=1 for all x∈ <cit.> (we recall that the solution to SPDEs of the form (<ref>) have mean 0). As a result, the marginal variance of u is  𝔼[u(x)^2]=σ_u^2 for all x∈.

In the case where  is a bounded domain, there is a boundary effect that affects the marginal variance of u. However, at a distance larger than twice the correlation length from the boundary, the bounded domain model is almost indistinguishable from the unbounded domain model. As a result, the boundary effect can be made negligible by embedding the region of interest in a sufficiently large domain.




§ PARAMETERIZATION

In this section, we parameterize H so that the parameters convey intrinsic geometric meaning about the field u. Recall conditions <ref>-<ref> and suppose for example that κ is fixed to 1 so that Ψ= √(H). Write 
    {(v,λ^2),(v_⊥,λ^-2)}
 for the eigensystem of H, where v:=(v_1,v_2)∈ℝ^2 and v_⊥=(-v_2,v_1) and we can suppose λ≥ 1 by reordering if λ< 1. Then, the previous discussion shows that  u corresponds to deforming and rescaling the stationary field u through

    u(x)= σ_u ·u(Ψx),  Ψx= λx,vv+λ^-1x,v_⊥v_⊥.

The rescaling by σ _u  corresponds to changing the variance of the field. The deformation corresponds to stretching the initial domain by a factor of λ in the direction of  v and contracting, also by a factor of λ, in the orthogonal direction v_⊥. The above shows that the eigensystem of H carries fundamental geometric information and motivates a parameterization of H in terms of its eigensystem.

Equation (<ref>) was also considered in <cit.>. Here, the authors defined v(α):=(cos(α),sin(α)) and parameterized H as

    H_v(α)=γI+βv(α)v(α)^T.

However, the map α↦H_v(α) is not injective as H_v(α)=H_-v(α). Because of this, it is impossible to recover the sign of v and, thus, this parameterization is not identifiable and leads to a bimodal likelihood.

The crucial step to obtain an identifiable parameterization is to consider the “half-angle” version v of v as an eigenvector of H.

    Given v=(v_1,v_2) ∈^2 define
    
    v:= vexp(i α /2 ),    where α := (v) ∈ [0,2 π).

    Write v=(v_1,v_2) and v_⊥=(-v_2,v_1).
    Then, the following defines a smooth, invertible parameterization on the space of symmetric positive definite matrices of determinant 1.
    
    H_v    =exp(v)/v^2 vv^T+exp(-v)/v^2 v_⊥v_⊥^T
            =cosh( |v|) I + sinh( |v|)/|v|[  v_1  v_2;  v_2 -v_1 ].


The link with the parameterization in <cit.> is the following

    Let H_v(α),H_v be as in (<ref>), (<ref>), then if we set
    
    v(α)=±v/v,   γ = exp(-v),  β =(1-γ^2)/γ.

    We obtain H_v(α)=H_v.


The idea behind the usage of the half angle version v of v is that it avoids any issues of identifiability in the sign as v and -v do not simultaneously belong to the parameter space. Parameterization using a Cholesky decomposition of H is also possible and is more readily generalized to higher dimensions. However, it is not as easily interpretable in terms of intrinsic properties of u.

In <Ref>, we show the half-angle vector field together with the parameterized diffusion matrices H_v. Here, each H_v is represented by the ellipse centred at v  whose main axis is exp(v)v/v and whose secondary axis is exp(-v)v_⊥/v. That is, the axes of the ellipses correspond to the eigenvectors of  H_v scaled by their respective eigenvalues. <Ref> shows visually how the anisotropy increases with v and is directed towards v. It can also be seen how the parameterization is injective (no two ellipses are the same) and smooth (the ellipses vary smoothly with v).


In <Ref> (a), we show the plot of the covariance K(x,0)=𝔼[u(x)u(0)] for κ=1 constant and as v is rotated around the X-axis by 90^∘  in each plot. As we can see, the vector field is most correlated in the direction of v, which is rotated by 45^∘ in each image, half the speed of rotation of v, and is least correlated in the direction of v_⊥. In Figure <Ref> (b), we show simulations of the field where again we leave κ=1 constant and rotate v by 90^∘ in each plot. The figure shows that the field diffuses the most in the direction of v and the least in the direction of v_⊥. The realizations of u are obtained using a FEM to solve the SPDE as detailed in <Ref>.


Another method that could be used to simulate the field u  and obtain the variance plots is the spectral method. The spectral method uses the spectral density of the field (the Fourier transform of the covariance function r if it exists)

    S(ξ):= ∫_exp(-2π iξ·x)r(x),   ξ∈,

to sample from the field, where  is Fourier domain (that is, =^d if ⊂^d is a box and =^d if =^d.). The field u is then sampled using the stochastic integral.
The spectrum for stationary fields u solving (<ref>) is given by (see <Ref> and <cit.>)

    S(ξ)=4πκ^2σ_u^2/(κ^2+4π^2ξ^THξ)^2.

Then, the covariance can be calculated as the Fourier transform of the spectrum, whereas the field u can be sampled using the stochastic integral

    u(x)=∫_𝒟u(x) 2π iξ·xdZ(ξ),

where Z(ξ) is called the spectral process. Using the fast Fourier transform, high-resolution samples of u can be obtained. See <cit.> for the details.

Thus far, we have worked with constant v. The same parameterization goes through when H(x) is allowed to be spatially varying. In this case, v(x) is a spatial vector field. In <Ref>, we take κ= σ_u=1 and show the covariance K(x,(2,2)) and field u  when v(x) is chosen to be the “twice-angle” field of the rotational field v(x)=(-x_2,x_1) (left of each subfigure), and from when v(x)=(x_1,x_2) (right of each subfigure). The figures show how the information of the field diffuses infinitesimally in the direction of v. The covariance and samples are once more obtained using the finite element method. In this non-stationary setting, the spectral method cannot be used as the spectrum is only defined for stationary fields.



In summary, we have parameterized the anisotropic field u using parameters (κ, v_1,v_2,σ_u), where u solves (<ref>) and H:= H_v is given by (<ref>). The parameterization is identifiable and smooth, and the parameters have an intrinsic geometric interpretation. This will be the parameterization used throughout the paper.


§ PENALIZED COMPLEXITY PRIORS

PC priors were originally developed in <cit.> to construct weakly informative priors while adhering to certain principles. The main idea is that one has a parametric family
of models M_θ with parameter θ and a base model M_0 (by convention corresponding to θ=0).

One views M_0 as the most suitable in the absence of any information. The larger the distance ζ(M_θ, M_0) between a model M_θ and M_0, the smaller the prior density for  θ should be, and the decrease is set to be exponential. That is, we choose the prior distribution for θ, such that:

    d(θ ):=ζ  (M_θ,M_0)∼Exp(λ_θ ).

Here, the rate λ_θ > 0 of the exponential distribution serves as a hyperparameter selected by the user, which governs the model's flexibility. A smaller value of λ_θ increases the model's flexibility, allowing for greater deviations of M_θ from the base model M_0. Conversely, a larger value of λ_θ imposes stricter penalties on these deviations, thereby reducing the model's flexibility. In the previously cited  <cit.>, this notion of “distance” was taken to be

    η (M_θ,M_0):=√(2 KLD(M_θ  | |M _0)),

where  KLD is the Kullback-Leibler divergence

    KLD(M_θ | | M_0 ):=∫_E log(Ṃ_θ/Ṃ_0 )Ṃ_θ .

One possible complication of using the KL divergence to define the distance is that the Radon–Nikodym derivative Ṃ_θ/Ṃ_0 does not exist. In fact, in the case of our model, if the variances of M_θ and M_0 are equal, for any possible base model u_0 and θ≠0 the measures M_θ and M_0 are singular. That is, with the usual convention,

    KLD(M_θ | | M_0 )=∞ .

The proof of this fact is given in <Ref>.
Equation (<ref>) makes an exact adherence to the previous steps impossible. As a result, we will merely adopt the principles of the PC prior construction (exponential penalization of complexity) while altering how we measure complexity
between models. The idea of modifying the metric used to measure complexity has also been carried out in different settings, such as in <cit.>, where a Wasserstein distance was used. The possibility of using the Wasserstein distance to measure the complexity of the model (<ref>) was also considered. However, the Wasserstein distance is bounded in this setting (see <Ref> in the supplementary material) and, as a result, was discarded. Thus, one of the main challenges is to find a computationally feasible distance which captures relevant information about our model. To this aim, we give the following definition.

    Given two sufficiently regular models u_A∼ M_A,u_B ∼ M_B, with respective spectral densities S_A, S_B and variances σ_A,σ_B, we define the pseudometric
    
    D_2(M_A ,M_B):=∫_^22πξ^4S_A(ξ)/σ_A^2-S_B(ξ)/σ_B^2^2 ξ^1/2 .


The above definition uses the rescaled spectral densities of each field to define a Sobolev seminorm on the difference of the correlations of M_A, M_B. For more details on how this distance was chosen and other distances that were considered, see <Ref>.

Write M_κ,v for the model given by SPDE (<ref>) with parameters (κ,v) and with variance fixed to σ _u=1. Due to the rescaling of the spectral densities in <Ref>, the choice of σ _u  does not affect the distance between models and can be set to any other positive value. We define the base model as the limit in distribution of M_κ,v as κ,v go to zero

    M_0:= lim_(κ_0,v_0)  →0M_κ_0,v_0 =0,1⊗1,

where 1: 𝒟→ℝ is the constant function , z ∈↦ 1 and (f⊗ g)(x,y):=f(x)g(y).
That is, u_0 is constant over space and follows a Gaussian distribution with variance 1. We choose this u_0 as our base model as it is simple, and (κ,v) =0 is the only distinguished point in the parameter space.

To reflect the dependency of this distance on the parameters, we use the notation

    d(κ,v):=D_2(M_κ,v,M_0) := lim_(κ_0,v_0)  →0D_2(M_κ,v,M_κ_0,v_0).



The exponential penalization in (<ref>) imposes one condition on the prior of (κ,v) whereas, since (κ,v) is three dimensional, two more conditions are necessary to determine the prior distribution uniquely. In <cit.>, this issue was circumvented by working iteratively, fixing one parameter while letting the other vary, building PC priors on each parameter, and then multiplying them together to get a joint prior. However, we prefer to work jointly from the start. This approach is made possible by the structure of d, which is calculated to have the form

    d(κ,v)=f(v)g(κ)

Since the angle α:= (v) does not affect (<ref>), we impose that α is uniformly distributed in [0,2π). This choice guarantees that we do not favor the alignment of the covariance of u in any direction of the plane. Next, since the contribution of f and  g to the distance is symmetric, we impose that  f-f(0) and  g knowing  f are exponentially distributed. The translation is necessary as f takes a nonzero minimum f(0)  at  v= 0, whereas  g takes a minimum of  0 at  κ =0.

This construction leads to the distance d(κ,v ) being exponentially distributed and gives us three conditions that uniquely determine (κ,v).


    A PC prior for (κ,v) with base model (κ,v) =0 is
    
    π(κ,v)     =λ_θλ_vf'(v ) f(v ) /2πvexp(-λ_v( f(v )-f(0) ) -λ_θ f(v ) κ ),

    where λ_θ>0 ,λ_v>0 are hyperparameters and
    
    f(r):=(π/3  (3 cosh (2 r )+1))^1/2,    f'(r)=√(π)sinh (2 r)/√(cosh (2 r)+1/3).


See <Ref> for the proof of this and other results of this section.
We plot the marginal density of the prior on κ and v  for λ_θ=λ_v=16π^2 in <Ref>. The marginal prior densities take a maximum at κ =0 and v=0  respectively, and by construction, the prior on v is radially symmetric.



The hyperparameter λ_θ determines the flexibility of the model (how much we penalize large values of d(κ,v )), whereas λ_v controls the degree of anisotropy (how much we penalize large values of v). Their values can be set to agree with desired quantiles using the following two results.

    The prior for r:=v satisfies ℙ[r  >r_0]=β if and only if
    
    λ_v =-log(β )/f(r_0)-f(0) .



 The prior for κ satisfies ℙ[κ >κ _0]=α if and only if
    
    λ_θ = 1/κ_0(1/f(0) W_0(λ_v f(0)λ_vf(0)/α ) -λ_v),

    where W_0 is the principal branch of the Lambert function. That is, W_0(x) is the real-valued inverse of  xx for x≥ 0,
    
    x=W_0(x)W_0(x),   ∀ x ≥ 0.


When specifying values of λ_θ,λ_v, it is useful to consider the ratio between the eigenvalues of H_v and the empirical correlation range. These measure respectively how much more correlated the field is in one direction in space and the distance at which the field becomes essentially uncorrelated (see Section <ref>). We denote these by

    a:= exp(v),   ρ:=√(8)κ^-1 .

<Ref> and  <Ref>  can then be rewritten as follows.

    The PC priors  in satisfy that ℙ[a>a_0]=β and ℙ[ρ <ρ_0]=α if and only if
    
    λ_v =-log(β )/f(log(a_0)) -f(0),   λ_θ = ρ _0/√(8)(1/f(0) W_0(λ_v f(0)λ_vf(0)/α ) -λ_v)
            .


In a practical application, α,β can be chosen to be small (for example 0.01), a_0 can be chosen to be an unexpectedly large amount of anisotropy, and  ρ_0 can be chosen as a surprisingly small correlation range for the field.

The parameters (κ,v) can be written as a joint transformation of a three-dimensional vector with standard multivariate Gaussian distribution. Thus, it is possible to efficiently generate independent samples (κ,v) through sampling multivariate Gaussian random variables, which is straightforward.


    Let Y∼(0, 𝐈_3) and write respectively Φ, R for the CDFs of a univariate standard Gaussian and Rayleigh distribution with shape parameter 1. Define
    
    A := √(Y_1^2 + Y_2^2),    B := f(0) - log(1 - R(A))/λ_v.

    Then, it holds that
    
    (v_1, v_2, κ) d=φ(Y_1, Y_2, Y_3) := ( f^-1(B) Y_1/A, f^-1(B) Y_2/A, -log(1 - Φ(Y_3))/λ_θ B),

    where f^-1(x) = 1/2cosh^-1( 4x^2 - 1/3).


The transformation φ only involves standard functions and can be evaluated efficiently. For the proof, see <Ref>.


§ SIMULATION STUDY



 §.§ Framework

This section aims to study the performance of the PC priors. To do so, we set different priors π_κ, v  on κ, v, observe a noisy realization of the field, and compare the behavior of the posteriors, resulting from each prior. We consider the random field u that solves our model (<ref>) on the square domain =[0,10]^2 and define the
observation process

    y= A u  +ε,

where u∈ℝ^n is a discrete approximation of the solution to (<ref>) obtained through a FEM on a mesh with n nodes, and A∈^m × n linearly interpolates to observe u at m=15 locations x_j_j=1^m   which are obtained by sampling uniformly from , and ε∈^m is a noise vector with

    ε∼(0,Q_ε ^-1),  Q_ε:=σ_ε^-2I_n.

The details of how u is sampled from and more detailed results can be found in <Ref>.
We set independent priors

    (κ , v) ∼π_κ , v,   σ_u∼Exp(λ_σ_u),   σ_ε∼Exp(λ_σ_ε).

The priors on σ_u, σ_ε correspond to their respective PC priors (see respectively <cit.> Section 3.3, <cit.> Theorem 2.1) with   hyperparameters λ_σ_u, λ_ε chosen so that

    ℙ[σ_u>σ_0]=0.01,   ℙ[σ_ε>σ_1]=0.01,

where we take σ_0=10,σ_1=1.5.

The prior π_κ, v is one of the priors to be compared (among them the PC prior). We consider the following options for π_κ , v:

    
  *  The PC priors in (<ref>) where the anisotropy hyperparameters λ_θ , λ_v are chosen so that the anisotropy ratio a=exp(|v|) and the correlation range ρ = √(8)κ ^-1 satisfy
          
    ℙ[a> a_0]= 0.01,   ℙ[ρ < ρ_0]=0.01,

          where we take a_0=10, ρ_0=1. This choice of hyperparameters corresponds to allowing with probability 0.01 that the field is 10 times more correlated in any given direction and, with the same probability, that the field has a correlation range smaller than 1.     The values of λ_θ,λ_v  can then be calculated using <Ref>.
    
  *  Independent priors κ∼Exp(λ_κ ) and v∼(0,σ_v^2 I_2 ). Under these priors, κ,v have the same mode as the PC priors (0 in each case). Additionally, λ_κ ,σ_v ^2 are chosen such that, under these priors, (<ref>) also holds. We denote this prior by π_EG.
    
  * Independent (improper) uniform priors for  log(κ),v with infinite support. We denote this prior by π_U.
    
  * Independent linear transformations of beta priors on log(κ),v_1,v_2 with shape parameters 1.1 such that the correlation range is supported in [ρ_0/w, wL]  and v_1,v_2  are supported in [-wa_0,w a_0]. The shape parameter is chosen so that the distribution is approximately uniform while having smooth density, which is relevant for the optimization. The purpose of w>1 is to extend the support of the parameters past ρ_0,a_0. The same value of ρ_0=1 is taken, L=10 is taken to be the length of , and w is set to 20. We denote this prior by π_β.

We simulate  θ^(j)^true from  π_sim∈{π_PC, π_EG, π_U, π_β}, use the FEM to simulate u^(j) from (<ref>) and then simulate y^(j) from (<ref>). Then, for each of the four priors π_κ,v∈{π_PC, π_EG, π_U, π_β}  we approximate the posterior distribution of the parameters θ=(κ, v, σ_u, σ_ε) given the data y^(j) and evaluate the performance of each of the four priors. This process is repeated J=600 times and repeated for each of the four possible values of π_sim. Since the posterior is not available in closed form, it is necessary to approximate it. We considered three choices: approximating the posterior by the Gaussian Z with the same median and whose precision is minus the Hessian of the posterior at its median, using importance sampling with Z as the proposal distribution, and using smoothed importance sampling with Z as the proposal distribution. Of these three approximations, the best performing method was smoothed importance sampling. As a result, this is the method we use in the following section to approximate integrals against the posterior. The details of the different approximations can be found in <Ref> in the supplementary material.



 §.§ Results

The focus of our study is the anisotropy parameters (κ, v_1, v_2). As a result, in this section, we will focus on the performance of the different priors on these parameters and relegate the results for the remaining parameters to <Ref>.

We first show in <Ref> the empirical cumulative distribution function (eCDF) of the vector of distances of the true anisotropy parameters
(log(κ^true),v_1^true,v_2^true) to the MAP estimates (log(κ),v_1,v_2) for each of the four different distributions on θ^true.

We observe that π_PC and π_EG perform the best and give almost identical results, to the point where it is difficult to distinguish them from the plots. This behavior is to be expected as both these priors are similar and lead to almost identical posteriors, as is discussed further in  <Ref>.

The difference in performance between π_PC, π_EG and π_U, π_β   is clearest for v_1,v_2  whereas for log(κ), when θ^true∼π_U or θ^true∼π_β all four models give comparable results. The figures relative to the parameters v_1 and v_2 are almost identical. This is expected by the symmetry in the priors and the likelihood, as there is no preferred direction of anisotropy


Next, in <Ref>, we show the length of the symmetric equal-tailed 0.95 credible interval for each parameter. As can be seen from the plots, both π_PC and π_EG limit the length of the credible intervals to a similar extent while π_U,π_β give much wider credible intervals. This connects with the motivating factor for the construction of the PC priors which is to penalize the complexity of the model. In <Ref>, we show the eCDF of the posterior mean complexity _π _θ|yd(κ , v) , where the complexity d  is defined in (<ref>) for each of the four different true distributions on θ. As can be seen, in all four cases, the PC and exponential-Gaussian priors reduce model complexity as compared to the uniform and beta priors.







§ AN APPLICATION TO RAINFALL DATA



 §.§ Framework

In this section, we analyze the performance of the anisotropic model and the PC priors on a data set for total annual precipitation in southern Norway between September 1, 2008, and August 31, 2009. This data set was studied in <cit.>, <cit.> using a linear model.

    y_i = β _0 + β _1 h_i + u(x_i)+ ε _i,    i=1,…,m= 233,

where y_i is the total annual precipitation at location x_i, h_i is the altitude at location x_i, u(x_i) is a spatially correlated random effect, and ε∼(0,σ_ε^2I) is a noise term.
In the articles above, u was modeled using a non-stationary Matérn process

    (κ^2(x)-Δ)(u(x)/σ_u(x))=√(4π)κ(x).

The non-stationarity κ^2(x),σ_u(x)  was parameterized by a log-linear model with covariates elevation and gradient of elevation and, in <cit.>, priors motivated through the PC prior approach were constructed for these extra parameters.




We will only study the stationary anisotropic setting where u| θ is a solution to (<ref>) with spatially constant parameters.  By incorporating β:=(β _0,β _1) into u, our linear model (<ref>) fits into the framework of <Ref>, where now

    A_β := (1_m, h,A),   u_β := (β, u)

take the place of A, u in (<ref>). We will consider

    β∼ (0, Q_β^-1),   Q_β = τ_βI_2

independent from u, θ and where we set the precision parameter to be τ_β= 10^-4.


 §.§ Maximum a posteriori estimates

We compare the model in (<ref>) under anisotropic PC and EG priors (see <Ref> and <Ref>) on (κ, v) and isotropic PC priors where v is set to 0. To do so, we must first derive an isotropic PC prior on κ using the distance metric (<ref>) restricted to the case where v=0. Using equation (<ref>) we obtain the following result.

    D_2 (M_κ , 0,M_0)= 1/√(12 π)κ.

So, by the principle of exponential penalization, κ∼Exp(λ_iso ).
As in <cit.>, we choose the hyperparameters in the three priors so that
    P(ρ<10)=0.05,   P(σ_u>3)=0.05,  P(σ_ε>3)=0.05.

For the anisotropy parameter v we choose the hyperparameter λ_v so that ℙ[a>10]=0.05. This choice of hyperparameters imposes that with probability 0.05, the field has a correlation range that is 10 times larger in any given direction.

In <Ref>, we plot the density of the priors on ρ. Since the marginal density of ρ is the same for the isotropic PC and EG priors, we do not plot it. We also plot the density of the PC and (anisotropic) EG priors on r:=v.

The decay of the marginal PC prior on κ and of the EG prior are both exponential. The decay rate of the marginal PC prior on v is exp(- c_1 exp(v )), whereas for the EG prior it decays slower, as exp(-c_2 v^2) for some constants c_1,c_2.

The MAP estimates, and symmetric 95% credible intervals for the anisotropic PC, anisotropic EG, and isotropic PC models are shown in <Ref>.



The credible intervals for v_1 in the anisotropic models do not contain 0, indicating that with high confidence, anisotropy is present in the precipitation field.

The half angle vector of the MAP for the anisotropic model with PC priors v  is  v=(0.02,0.48), and indicates that the precipitation is a=1.64 times more correlated in the North-South direction than in the East-West direction.
In <Ref>, we plot the posterior prediction, latent field, and the covariance function of the anisotropic field u with parameters κ,v.




 §.§ Model performance


To assess the performance of the models, we calculate a variety of scores for each model. We recall that, given a (predictive) distribution P and an (observed) point y, a score is a function S(P,y) that measures how closely the prediction P matches the observation y. If y follows the distribution Q, then the expected score is

    S(P,Q) := _Q[S(P,y)].

The score should be minimized when the true distribution Q matches the predictive distribution P. In this case, the score is said to be proper. It is strictly proper if it is minimized if and only if P=Q. We will consider the following scores

    
  * Given a distribution F̣ and an observation x ∈ the squared error (SE) is defined as

    SE(F,x)     := (x-∫_-∞^∞ tF̣(t))^2.

The expression above defines a proper scoring rule (its value is minimized when x follows the predictive distribution). However, it is not strictly proper (the minimum in F  is not unique).

  * Given a CDF F and an observation x ∈, the continuous ranked probability score (CRPS) is defined as  <cit.>

    CRPS(F,x)     := ∫_ (F(t)-1{x ≤ t})^2 d t.

The CRPS is a strictly proper scoring rule on distributions with finite expectations.

  * Given a random variable X with mean μ and precision Q, the Dawid-Sebastiani score (DSS) is defined as   <cit.>

    DSS(X;x)     := -log(Q ) + (x - μ)^TQ(x - μ).

The DSS is similar to a log Gaussian density score and is a strictly proper scoring rule for Gaussian random variables. Still, it also defines a proper score for non-Gaussian random variables.

Given a sample sample y=(y_1,...,y_n), a vector of predictive distributions F=(F_1,...,F_n) and a score S , the mean score is defined as

    S(F, y):= 1/n∑_i=1^n S(F_i,y_i).

By the linearity of the expectations, if S is a (strictly) proper scoring rule, then so is S.


 §.§ Score results

In this section, we calculate the mean scores

    RMSE:=(SE(F,y))^1/2,   CRPS(F,y),   DSS(F,y),

where F_i:= π(y_i| y_-i) is the LOO predictive distribution under each of the three models in Subsection <ref>, and y is the observation.


An initial calculation showed that the scores of the previous sections are very similar under all three priors. From this, we conclude that the model is very informative, and the influence of the priors is limited. As a result, to better compare the priors, we uniformly sub-sample y  and observe only y'  ∈^n_y' with n_y'< n_y=233. We then calculate the scores from the previous sections. Due to the extra variability introduced by sub-sampling the observations, we repeat this process 10 times. The resulting mean scores are shown in <Ref>. The anisotropic models perform better with lesser data, with the PC prior performing slightly better than the EG prior. Whereas with more data the results are almost equivalent.


From the trend of the above results it is reasonable to check whether the isotropic model outperforms the other two models with a larger amount of observations. To check this hypothesis we conducted a simulation study where a larger amount of observations are generated, up to n_y=1000. The results are included in <Ref> in the supplementary material and show that this is not the case. The isotropic model does not outperform the anisotropic models with a larger amount of observations.



§ DISCUSSION

In this study, we constructed a smooth, identifiable, and geometrically interpretable parameterization for a 2D anisotropic spatial Gaussian model. We developed penalized complexity (PC) priors for these parameters, defining model complexity as a Sobolev seminorm of the correlation function. This distance was calculated in a closed form using the spectral density.

Our performance comparison of the PC prior against other priors demonstrated its effectiveness in penalizing model complexity. We found that priors designed to match the quantiles of the PC prior yielded similar results. Both these priors significantly outperformed the other non-informative priors considered, highlighting the necessity of incorporating penalization in prior information to prevent overfitting.

Applying the anisotropic model to a real dataset of precipitation in Southern Norway, we observed that the anisotropic model outperformed the isotropic model when the number of observations was small, with the PC prior slightly outperforming the other choices. However, with a larger number of observations, the isotropic model performed similarly to the anisotropic models. These results indicate that the anisotropic model is more informative when the data is scarce, but as the number of observations increases, the isotropic model becomes more competitive. The results also suggest that a nonstationary model could be better suited to capture spatially varying patterns in the data.

In conclusion, we advocate for the use of informative priors for the anisotropy parameters in spatial Gaussian models. The PC prior is highly effective, but other priors designed to match desired quantiles can also be effective and are simpler to construct.

Looking forward, we aim to extend this study to the non-stationary setting where the model parameters vary spatially. This presents a challenge as there is no agreed-upon definition of correlation function or spectral density, necessitating a different definition of complexity. We are also interested in extending the parameterization to higher dimensions. The current construction relies on a “half-angle” parameterization of the anisotropy vector, which is not easily extendable to higher dimensions. Finally, we are interested in extending the construction to different orders of regularity.






§ SUPPLEMENTARY MATERIAL

tocsectionSupplementary Material 






§ PROOF OF <REF>


    By definition, H_v is symmetric with determinant 1. Furthermore, it is positive definite as its eigenvalues are positive. Given a positive definite symmetric matrix A with determinant 1, A has an eigensystem of the form
    
    (w,λ ),(w^⊥,λ ^-1),

    where λ≥ 1.
    By normalizing and taking -w if necessary we may suppose that w=α i with α∈ [0,π). Then, we obtain A = H_v for
    
    v=log(λ)2α i.

    This proves that the parameterization is surjective. Suppose now that H_v=H_v' then their eigenvalues and eigenvectors must be equal so
    
    v=v',   v=av'.

    For some a ∈. From the first condition, we deduce that v=v'. Taking absolute values in the second condition thus gives a=1. By construction v≠ -v' for all v,v' so necessarily  v=v'. This proves that the parameterization is invertible.

    Next, to derive the second equality in (<ref>), we use the half-angle formula
    
    cos(α/2)= √(1+cos(α ) /2)sign(sin(α))  ,   sin(α/2 )= √(1-cos(α) /2).

    This gives that
    
    vv^T= v^2/2I+ v/2[  v_1  v_2;  v_2 -v_1 ],   v_⊥v^T_⊥= v^2/2I-v/2[  v_1  v_2;  v_2 -v_1 ]    .

    The proof follows immediately by using the definition of H_v (first line of (<ref>)).



§ SELECTION OF A DISTANCE

In this section, we discuss the choice of the Sobolev distance in <Ref>. We first examine other possible choices, such as the Kullback-Leibler divergence, the L^2 distance, the Wasserstein distance, and the Hellinger distance, and show that they are unsuitable for our purposes. We then derive the exact form of the Sobolev distance in <ref> and show a possible alternative definition in <ref>.

Without pretending full generality we will use the following notation. Firstly, given a domain =[a,b] ⊂^d we denote the orthonormal Fourier basis on  by

    e_k(x):= vol( )^-1/2exp(ω_k·x/b-a),   ω_k:= 2 π i k  k∈^d,

where we used the elementwise division [x/(b-a)]_i:= x_i/(b_i-a_i). Given a stationary field u on , we denote its covariance function r and spectral density, if they exist, by

    r(x):= Cov[u(x),u(0)],    S(k):=∫_  r(x)e_k(x)x.

Additionally, we will denote the covariance operator of u by K. For stationary u ∈ L^2([a,b]) this is given by

    Kf,g_L^2( ):= ∫_∫_ r(y-x)f(x)g(y) xy,    f,g ∈ L^2( ).

For the calculations, the following lemma will be very useful.

    Let u be a stationary field in L^2() with spectral density S. Then, the covariance operator of u diagonalizes in the orthonormal Fourier basis e_k_k∈^d  with eigenvalues S(k).


    By the Bochner theorem, the covariance function r of u is the Fourier transform of its spectral density S. As a result, for all f,g ∈ L^2( ) we have
    
    Kf,g_L^2( )    =∫_∫_ r(y-x)f(x)g(y) xy=∫_∫_∑_k∈^d S(k)e_k(y-x)f(x)g(y) xy
        =∑_k∈^d S(k)∫_ f(x)e_k(x)x∫_ g(y)e_k(y) y=∑_k∈^d S(k)f(k)g(k),

    Applying this to f = e_j and g = e_k shows that
    
    K e_j,e_k_L^2( ) = S(k)δ_jk.

    This concludes the proof.

We will also repeatedly use the expression of the spectral density of the solution u  to (<ref>). To calculate this, note that the covariance operator of u is given by K=  ^-2 where

    := κ^2- ∇·H∇/√(4 π)κσ _u.

Let us consider first the case where = a,b⊂^d.
Using that ∇ e_k= i ω_k e_k we obtain that  diagonalizes in the Fourier basis

    [ ]_jk:=  e_j,e_k_L^2( ) = κ^2- ω_k·Hω_k/√(4 π)κσ _uδ_jk,    [ ^-1]_jk= √(4 π)κσ _u/κ^2- ω_k·Hω_kδ_jk.

As a result, by <Ref>, the spectral density of u is given by

    S(k)= 4 πκ ^2 σ _u^2/(κ^2- ω_k^T Hω_k)^2,   k∈^d.

If now =^d, the defining property of the spectral density is that, analogously to (<ref>), it holds that <cit.>

    Kf,g_L^2(^d )= ∫_^d S(ξ)f(ξ)g(ξ)ξ,   ∀ f,g ∈ L^2(^d).

Using that the Fourier transform is an isometry on L^2(^d),

    Kf,g_L^2(^d )    =  ^-1 f,  ^-1g_L^2(^d )= ^-1f,^-1g_L^2(^d )
       =∫_^d4πκ ^2 σ _u^2/(κ^2- 4π^2ξ^T Hξ)^2f(ξ)g(ξ)ξ.

So the stationary solution to (<ref>) on ^d  has spectral density

    S(ξ)= 4πκ ^2 σ _u^2/(κ^2- 4π^2ξ^T Hξ)^2.







 §.§ Kullback-Leibler divergence

As previously discussed, the KLD is unsuitable as a notion of distance between models as it is infinite for every possible parameter value. As a result, a different way to measure distance between models is necessary. This section discusses some of the options considered and how <Ref> was eventually chosen.
We begin by showing that the KLD is infinite for every possible parameter value with the same variance. To do so, we first give a sufficient condition for two stationary measures to be mutually singular. This is essentially the converse direction of  <cit.>. See also <cit.>.

    Write =[a,b]⊂^d and let u_A,u_B be two stationary Gaussian fields in L^2( ) with spectral densities S_A,S_B. Suppose that
    
    ∑_k ∈^d(S_A(k)/S_B(k)-1)^2 < ∞.

    Then, the Gaussian measures defined by u_A,u_B are mutually singular.


    Denote the covariance operators of u_A,u_B by K_A,K_B respectively and let I be the identity on L^2( ). By the Feldman–Hájek theorem, a necessary condition for u_A,u_B to not be mutually singular is that the following operator is Hilbert Schmidt,
    
    T:=(K_B^-1/2K_A^1/2)(K_B^-1/2K_A^1/2)^*-I.

    is a Hilbert-Schmidt operator where K_A,K_B are the covariance operators of u_A,u_B and I is the identity operator on L^2( ). By <Ref>, K_A, K_B both diagonalize in the orthonormal Fourier basis with eigenvalues given respectively by S_A(k) and S_B(k). As a result, the Hilbert-Schmidt norm of T  is given by
    
    T_HS^2= ∑_k∈^d(S_A(k)/S_B(k)-1)^2.

    The result follows by the necessary condition for mutual singularity.

Using the just proved <ref>, we show that the KLD between two different solutions to (<ref>) is infinite if they have the same variance parameter. In particular, if we were to renormalize the spectral densities of the two solutions to have the same variance, the KLD would be infinite for all different parameter values.

    Let u_A,u_B be the solutions to (<ref>) with Neumann or Dirichlet boundary conditions and parameters (κ _A,v_A,σ _A),(κ _B,v_B,σ _B) respectively. Then, the measures defined by u_A,u_B are mutually singular unless κ _A/κ _B=σ _A/σ _B and v_A=v_B. In particular, if u_A ≠ u_B both share the same variance then
    
    KLD(μ _A||μ _B)=∞.



    We have by (<ref>) that
    
    ∑_k∈^2(S_A(k)/S_B(k)-1)^2  = ∑_k∈^2(κ _A^2 σ _A^2/κ _B^2 σ _B^2(κ _B^2+4π^2k·H_v_Bk)^2/(κ _A^2+4π^2k·H_v_Ak)^2-1)^2 .

    Since H_v_A, H_v_B are positive definite, if their eigenvectors and eigenvalues are different, the term
    k·H_v_Bk/k·H_v_Ak
    can be made arbitrarily large by choosing k large enough. As a result, the sum diverges unless v_A=v_B and
    
    ∑_k∈^2(κ _A^2 σ _A^2/κ _B^2 σ _B^2(κ _B^2+4π^2k·H_v_Ak)^2/(κ _A^2+4π^2k·H_v_Ak)^2-1)^2 < ∞.

    The limit of the terms in the sum above as k→∞ is
    
    lim_k→∞κ _A^2 σ _A^2/κ _B^2 σ _B^2(κ _B^2+4π^2k·H_v_Ak)^2/(κ _A^2+4π^2k·H_v_Ak)^2-1= κ _A^2 σ _A^2/κ _B^2 σ _B^2-1.

    This limit is zero if and only if κ _A/κ _B=σ _A/σ _B. We deduce that the sum in (<ref>), and so also the sum in (<ref>), diverges unless κ _A/κ _B=σ _A/σ _B and v_A=v_B. The result follows by <Ref>.





 §.§ The L^2 distance

One possible option is to consider the L^2 distance between  u,u_0. An application of Fubini shows that given two Gaussian fields u_A∼(0,K_A),u_B∼(0,K_B)

    [u_A-u_B^2_L^2()]=Tr (K_A+K_B-2K_AB),

where K_AB is the covariance between u_A,u_B and Tr is the trace. It is unclear what choice of K_AB would be the most appropriate. For example, choosing  K_AB=0 would be too coarse a measure as it would not consider any of the non-stationarity that occurs off the diagonal of K_A, K_B (see Lemma <ref>). One can also choose K_AB as to minimize the L^2 norm while keeping (u_A,u_B) jointly Gaussian  by taking

    K_AB=(√(K_B)K_A√(K_B))^1/2.

This leads to the Wasserstein distance, which we discuss in the following subsection.



 §.§ Wasserstein distance between Gaussian measures

Consider a bounded domain =[0, T]^d. Each solution to SPDE (<ref>) induces a measure on L^2(). As a result, one possible way to measure the distance between u and u_0  is to take the Wasserstein distance between the induced Gaussian measures. It is known that given two Gaussian measures M _A∼(m_A, K_A),μ _B∼(m_B, K_B) on a separable Hilbert space, the Wasserstein distance between them is given by <cit.>

    W_2(μ _A,μ _B)=m_A-m_B_L^2()+Tr(K_A+K_B-2(√(K_B)K_A√(K_B) )^1/2),

where Tr is the trace. However, this approach poses some difficulties. Let us write K_A, K_B for the covariance operators of two stationary solutions u_A∼ M_A, u_B∼ M_B to (<ref>) with parameters (κ _A,v_A,σ _A),(κ _B,v_B,σ _B) respectively. By <Ref>, K_A, K_B diagonalize on the same basis, and using the expression of the spectral density in (<ref>), we obtain that

    W_2(M_A,M_B)     = (√(K)-√(K_0))^2= ∑_j,k∈^d[√(K)-√(K_0)]_jk^2                                                                                                                                     
        =∑_j ∈^d(√(4 π)κ_A σ_A /κ_A ^2+ 4 π^2ξ_j^T H_v_Aξ_j-√(4 π)κ_B σ_B /κ_B ^2+ 4 π^2ξ_j^T H_v_Bξ_j)^2 .

Making the domain  go to infinity shows that, by definition of ξ_j our discrete sum becomes an integral with

    W_2(M _A,M_B)=vol() /π∫_^2(κ_A /κ_A ^2+ ξ^T H_v_Aξ-κ _B/κ _B^2+ ξ^T H_v_Bξ)^2ξ +O(1).

The Wasserstein distance between μ _0,μ _1 scales as vol() times the Hellinger distance of the spectral densities of u_A,u_B. Thus, a reasonable option could be to define the scaling as the distance

    D_W_2 (M_A,M_B):= lim_T →∞W_2(M_A,M_B)/vol() =Hell(S_A,S_B).

However, a calculation in this simplified case shows that this will be bounded. For example, in the case where, as in the base model, v_A=v_B=0, and σ_A = σ _B =1, we have the following.

    D_W_2 (u_A,u_B)=2(1-κ_Aκ_B (log(κ _A^2)-log(κ_B ^2))/κ _A^2-κ_B ^2).

The above expression takes a maximum of 2. This makes putting an exponential distribution on the Wasserstein distance impossible, and as a result, the Wasserstein distance cannot be used to define PC priors for our model (<ref>).







 §.§ Sobolev distance between fields

Given a stationary field u on  with s-times differentiable covariance function r(x):= Cov[u(x,u(0))], we define the seminorm ·_s as the Sobolev seminorm of order s of r. That is,

    u_s:= ∇^s r_L^2()=2πξ^sS(ξ)_L^2( ):=∫_2πξ^2sS(ξ)^2 ξ^1/2,

where  is the Fourier domain of . That is, =^d if  is a box and =^d if =^d. If we think of the L^2() norm of r as a measure of the total “random energy” present in u. Then, ·_s measures the oscillation size in this energy to order s.

We now discuss which value of s is appropriate. In our case, as we saw in (<ref>), the solution u to (<ref>) with parameters κ, v, σ_u  has spectral density

    S_κ, v, σ_u(ξ)= 4πκ ^2σ_u^2/(κ ^2+4π^2ξ·H_vξ )^2 .

A change of variables ξ→κξ shows that, for dimension d=2,  the behavior in κ is

    u_s∼κ ^s-1 .

If we take s=0, then  u_0 becomes infinite, whereas if we choose s=1, we obtain that  u_1 is bounded in κ and thus cannot be exponentially distributed. This leads to the choice s=2 used in <Ref>.


    The seminorm u_s is different to the Sobolev  seminorm of u as we have that, given a multi-index α∈^d, the derivative  D^α u of a stationary field u  exists as long as
    
    ∫_(2πξ)^2α S(ξ)  ξ < ∞ .


Thus, one possible choice of metric, corresponding to s=1 /2, could also be to set

    u_*:=∫_2πξ S(ξ)  ξ   .

This distance is unbounded in κ and finite at the base model. A calculation gives a distance of a similar form to the one derived for the distance used in this paper in <Ref>.

    d_*(κ , v)^2     :=lim_κ_0,v→ 0 ∫_2πξ(S_κ, v, σ_u/σ_u^2-S_κ, v, σ_0/σ _0^2)  ξ=2π E(1-2v)-v/2κ,

where E is the complete elliptic integral

    E(m):=∫_0^π /2(1-msin^2(α ))^1/2α̣.

This distance is unbounded in κ and v and 0 at the base model. As a result, using d_* instead of d  could also provide a valid alternative.



§ DERIVATION OF PC PRIORS


    Let M_κ ,v,σ_u,M_κ _0,0,σ_0 be the stationary models given by (<ref>) with parameters  (κ,v,σ_u) and  (κ_0,0,σ_0) respectively. Then,
    
    d(κ ,v):=lim_κ_0 → 0D_2(M_κ ,v,σ_u,M_κ_0,0,σ_0)^2=π/3κ ^2 (3 cosh (2 | v| )+1) .



    For any (κ ,v,σ_u) the spectral density of u_κ ,v,σ_u is
    
    S_κ ,v,σ_u(ξ)=4πκ^2 σ _u^2/(κ ^2+4π^2ξ^T Hξ)^2 .

    As a result,  in a distributional sense
    
    lim_κ _0 → 0S_κ_0, 0,σ_0/σ_0^2(ξ)= lim_κ _0 → 04πκ_0^2 /(κ_0 ^2+4π^2ξ^2 )^2=δ0 (ξ),

    where δ0 is the Dirac delta at 0.
    In consequence, expanding the square in (<ref>) and using the change of variable ξ→κξ/(2π), we obtain that
    
    lim_κ_0 → 0D_2(M_κ ,v,σ_u,M_κ_0,0,σ_0)^2=∫_^22πξ^4 S_κ ,v, σ _u(ξ)^2/σ _u^4ξ  =4κ^2∫_^2ξ^4/(1+ξ^T H_vξ)^4ξ.

    Let α:= (v) and B_α /2 be a rotation by α /2 radians. Then
    
    B_α /2H_vB_α /2 ^-1= [  -v   0;   0 exv ] .

    Thus, by rotating and then changing to polar coordinates, we obtain
    
    lim_κ_0 → 0D_2(M_κ ,v,σ_u,M_κ_0,0,σ_0)^2=4κ^2 ∫_0^2π∫_0^∞r^5/(1+r^2 (| v| sin ^2(ϕ )+-| v| cos ^2(ϕ )))^4ṛϕ̣
        =∫_0^2π4κ ^2/6 (vsin ^2(ϕ )+-vcos
            ^2(ϕ ))^3ϕ̣=π/3  (3 cosh (2 | v| )+1)κ ^2 .

    This concludes the proof.

As we see from <Ref>, d(κ ,v) is independent of the argument α of v  and depends only on κ and r:=v. For this reason, we consider α to be independent of κ,v and, since a priori, there is no preferred direction for the anisotropy, we set a uniform prior on α

    π_α (α)=1/2π1_[0,2π ](α)

We will now define priors on κ,r following a sequential construction. To do so, we will use the fact that the distance to the base model ρ (κ,v ) decomposes as a product

    d(v,κ)  =f(r)g(κ),  
        f(r)            :=(π/3 (3 cosh (2 | v| )+1) )^1/2,   g(κ):=κ.

Together with the following lemma,

    Let X,Y be two positive random variables such that Y|X has density
    
    f_Y|X(y,x)=λ x exp(-λ x y )1_[0,+∞)(y).

    Then, the product Z:=XY follows an exponential distribution with parameter λ.


    By the law of total expectation. the cumulative density function of Z verifies
    
    F_Z(z)     =ℙ[XY≤ z]=[[1_XY≤ z|X]]=[∫_0^z/Xf_Y|X(y|X) dy] 
       =∫_0^∞∫_0^z /x f_Y|X(y|x)dℙ_X(x),

    where ℙ_X is the push-forward of ℙ by X.
    Applying the fundamental theorem of calculus shows that the density of  Z is
    
    F'_Z(z)=∫_0^∞1/xf_Y|X(z /x|x)dℙ_X(x)=λ-λ z∫_0^∞ dℙ_X(x)=λ-λ z.

    This concludes the proof.

The above lemma states that if we want the distance d=f g to be exponentially distributed, it suffices to let f follow  any distribution and have g conditional on  f be exponentially distributed with parameter proportional to  f. Since f takes a minimum of

    f(0)=√(4π/3 ),

the conditional distribution of f given g cannot follow an exponential distribution. Therefore, we will instead select the conditional distribution of g given f to be exponential.

While <Ref> allows complete freedom in choosing a distribution for f, the symmetry in the distance fg implies there is no intrinsic reason to prefer penalizing an increase in f either more or less than an increase in g. Consequently, we will also impose an exponential penalty on f by setting its marginal density to be equal to the following,

    π_f(f)= λ_v-λ_v (f- f(0))1_[f(0),∞).



We can set priors in the stationary setting in the same way as in the previous section.



 §.§ Proof of <Ref>
 

    Let us write
    
    r=v,   α =(v) .

    We have already determined the prior distributions
    
    π_f(f)= λ_v-λ_v(f- f(0))1_[f(0),∞)(f),  π_κ |f(κ |f)=λ_θ f -λ_θ f κ1_[0,∞)(κ )

    Since f:[0,∞) → [f(0), ∞) is bijective, we may apply a change variables to obtain that the marginal prior for r is
    
    π _r(r)=π_f(f(r))f'(r)= λ_v-λ_v( f(r)-f(0))f'(r)  .

    Again, by the bijectivity of f, we obtain that
    
    π_κ|r=π_κ|f=f(r)=λ_θ f(r) -λ_θ f(r) κ1_[0,∞)(κ ).

    Since α is uniformly distributed on [0,2π ] independently of r,κ the joint prior for (κ,α,r ) is
    
    π_κ,α,r (κ,α,r )=1_[0,2π](α)/2ππ_r (r)π_κ|r  (κ|r ) .

    Changing variables from polar to Cartesian coordinates gives
    
    π_κ,v(κ,v )=1/2πvπ_r (v )π_κ|r  (κ|v) .

    Using equations (<ref>) and (<ref>) in (<ref>) and using the expression for f derived in (<ref>)  concludes the proof. 



 §.§ Proof of <Ref>


    The prior for f is defined in (<ref>). Since f is increasing, we obtain that the cumulative distribution function of r is for all r_0≥ 0
    
    F_r(r_0)=ℙ(r≤ r_0)=ℙ(f(r)≤ f(r_0))=1--λ_θ (f(r_0)-f(0)).

    The theorem follows by a straightforward algebraic manipulation.




 §.§ Proof of <Ref>


    We begin by calculating the prior π_κ  for κ. A calculation shows that for all κ  >0
    
    π_κ(κ )     =∫_π_f(f)π_κ |f(κ  |f)df=∫_f(0)^∞λ_vλ_θ f  exp(-λ_v(f- f(0))-λ_θ  f κ  )df
            
       = λ_θλ_v-f(0) κλ_θ(f(0) κλ_θ +f(0) λ_v+1)/(κλ_θ +λ_v)^2,

    whereas for κ <0 we have π_κ (κ )=0.
    A further calculation shows that the CDF of κ is
    
    F_κ  (κ _0 )=∫_-∞^κ_0π_κ (κ)  κ̣= 1-λ_v-f(0) λ_θκ _0/λ_θκ _0+λ_v1_[0,∞](κ _0).

    Now solving for λ_θ gives that for any κ _0>0
    
    ℙ[κ >κ _0]=α⟺λ_θ = 1/f(0)κ _0 W_0(λ_v f(0)λ_vf(0)/α ) -λ_v/κ _0

    This concludes the proof.



§ JOINT TRANSFORMATION

In this section, we show how to write (κ,v) jointly as a function of a three-dimensional standard Gaussian Yd=(0, 𝐈_3). Our idea will be to generalize the method of inverse sampling. This is done by building an invertible generalization of the cumulative distribution function.

    Given a random variable X=(X_1,X_2) ∈^2 we define the generalized CDF of X as
    
    φ(X,x):=φ_X(x):=(ℙ[X_1≤ x_1],ℙ[X_2≤ x_2|X_1=x_1]).



    In the case where X_1, X_2 are independent we obtain that
    
    φ_X(x)=(F_X_1(x_1),F_X_2(x_2)).




By construction φ_X:^2→ [0,1]^2. The definition can be extended without difficulty to dimensions larger than 2. However, the notation becomes more cumbersome. The motivation for the above definition is given by the following properties, which mimic that of the 1-dimensional CDF.

    Let X be a two-dimensional random variable, then the following hold
    
        
  *  The distribution of X is determined by  φ. That is if Y is a two dimensional random variable such that φ_X=φ_Y, then Xd=Y.
        
  *  Consider a function ϕ:A→^2 defined on some subset A ⊂^2 of the form
              
    ϕ(x_1,x_2)=(ϕ_1(x_1),ϕ_2(x_1,x_2)),

              where ϕ_1 is invertible and  ϕ_1,ϕ_2(x_1,·) are monotone functions.
              Then,
              
    φ(ϕ(X),ϕ(x))=φ(X,x),   ∀x∈ A.

        
  * Let X be a  ^2 valued random variable. Then φ_X has a generalized inverse and φ_X, φ_X^-1 are both of the form  (<ref>).
    


    We prove the above points in order.
    
        
  * Suppose  φ_X=φ_Y, then by definition of φ
              we have that for all x=(x_1,x_2) ∈^2
              
    ℙ[X_1≤ x_1]=ℙ[Y_1≤ x_1],  ℙ[X_2≤ x_2|X_1=x_1]=ℙ[Y_2≤ x_2|Y_1=x_1] .

              From the basic theory of random variables, the left-hand side of the above implies that X_1d=Y_1. As a result, we obtain that
              
    ℙ[X_2≤ x_2]     =∫_ℙ[X_2 ≤ x_2|X_1=x_1]dℙ_X_1(x_1) 
       =∫_ℙ[Y_2 ≤ x_2|Y_1=x_1]dℙ_Y_1(x_1)=ℙ[Y_2≤ x_2].

              Since this holds for all x_2 ∈ we deduce as previously that X_2d= Y_2. Thus we obtain that Xd=Y, as desired.

        
  * To prove   <ref> we work with  φ componentwise. By the monotonicity of ϕ_1, it is clear that for the first component,               
    φ_1(ϕ(X),ϕ(x)):=ℙ[ϕ_1(X_1)≤ϕ_1(x_1)]=ℙ[X_1≤ x_1]=φ_1(X,x).

              The second part is proved similarly. We have that
              
    φ_2(ϕ(X),ϕ(x)):=ℙ[ϕ_2(X_1,X_2)≤ϕ_2(x_1,x_2)|ϕ_1(X_1)=ϕ_1(x_1)] 
       =ℙ[ϕ_2(X_1,X_2)≤ϕ_2(x_1,x_2)|X_1=x_1]=ℙ[ϕ_2(x_1,X_2)≤ϕ_2(x_1,x_2)|X_1=x_1]
                      
       = ℙ[X_2≤ x_2|X_1=x_1]=φ_2(X,x),

              where in the second equality, we used that ϕ_1 is invertible, so the σ-algebra generated by X_1 is the same as that generated by  ϕ(X_1) and ϕ(X_1)=ϕ(x_1) if and only if X_1=x_1. And in the fourth equality, we used that ϕ_2(x_1,· ) is monotone.
        
  * We first prove φ_X is invertible. Let us take p=(p_1,p_2) ∈ [0,1]^2. Then, the first component of φ_X is the univariate CDF F_X_1, and we can define
              
    x_1:=F_X_1^-1(p_1).

              Now, since φ_2(X,(x_1,· )) is increasing (as a function of the dot) for each x_1, it has a generalized inverse
              
    x_2:=φ_2(X,(x_1,· ))^-1(p_2).

              An algebraic verification now shows that φ_X has generalized inverse
              
    φ_X^-1(p)=(x_1,x_2).

              Furthermore, since φ_X is monotone in each component, so is φ^-1_X.
    
    This completes the proof.



    Given  random variables X,Y valued in ^2 it holds that
    
    φ_X^-1∘φ_Y(Y)d=X



    By <ref> of <Ref>, it suffices to prove that the generalized CDFs of X and  Y coincide. To this aim, we have that, by <ref> of the previous lemma
    
    φ(φ_X^-1∘φ_Y(Y),x)=φ(Y,φ_Y^-1∘φ_X(x))=φ_Y(φ_Y^-1(φ_X(x)))=φ_X(x).

    This concludes the proof.


We now calculate the transformation which makes (κ,r) Gaussian.

    Let (v,κ) be distributed according to the PC priors π (v,κ ) and define r:=v. Then, given U∼(0,𝐈_2)  it holds that
    
    (r,κ)d=( f^-1(f(0)-log(1-Φ(U_1))/λ_v ), (log(1-Φ(U_1)))/λ_v-f(0)  )^-1log(1-Φ(U_2))/λ_θ
            ),

    where Φ is the CDF of a univariate standard Gaussian and r:=v.


    The CDF F_r of r was calculated in (<ref>). An inversion gives
    
    x_1:=F_r^-1(p_1)                 = f^-1(f(0)-log(1-p_1)/λ_v  )

    To calculate the CDF of κ |r, we work with the density in (<ref>). This gives
    
    φ_2((r ,κ),(x_1,x_2))= 1--λ_θ f(x_1)x_2.

    For each fixed x_1, the above has the inverse
    
    x_2:=φ_2((r,κ),(x_1,· ))^-1(p_2)=- log(1-p_2)/λ_θ f(x_1)

    This, combined with (<ref>), the preceding <Ref> (with Y= U ), and <Ref> (applied to U), completes the proof.



 §.§ Proof of <Ref>


    Let α∼_[0,2π ] be independent from r. Then, by the definition r= v,
    
    vd=(rcos(α ),rsin(α )),

    The proof then follows from Lemma <ref> and by observing that R(√(Y_1^2+Y_2^2) ),Φ(Y_3) are i.i.d. uniformly distributed on  [0,1], and that
    
    Y_1/√(Y_1^2+Y_2^2)d=cos(α),  Y_1/√(Y_1^2+Y_2^2)d=sin(α)  .




 §.§ Proof of <Ref>



    The determinant of H_v(α) is (β+γ)γ. Setting this equal to 1 gives β=(1-γ^2)/γ. Next, we note that given any vector e, it holds that
    
    𝐈=ee^T+e_⊥e^T_⊥.

    Substituting this into (<ref>) with e=v(α)
    
    H_v(α)=(γ+β)v(α)v(α)^T+γv(α)_⊥v(α)_⊥ ^T=1/γv(α)v(α)^T+γv(α)_⊥v(α)_⊥ ^T.

    Setting H_v=H_v(α), we obtain the condition
    
    1/γv(α)v(α)^T+γv(α)_⊥v(α)^T_⊥=v/v^2 vv^T+-v/v^2 v_⊥v_⊥^T.

    Thus, equality holds in the conditions of the theorem.



§ SIMULATION STUDY DETAILS

In this section, we expand on the simulation study of <Ref>. We will discuss the details of the simulation of the SPDE, compare the priors, and show some additional results.


 §.§ Sampling from the SPDE 

To simulate the Gaussian field resulting from (<ref>), we use the  package. This package uses a FEM method to calculate to approximate the precision matrix of u. That is, we approximate u to be in the Hilbert space H spanned by our finite element basis ϕ_i_i=1^n and impose

    [ϕ_i,  u, i∈{1,...,n}]= [ϕ_i, κ, i∈{1,…,n}]

The condition that u∈ H leaves us with

    u(x)= ∑_i=1^n ϕ_i(x) u_i,

where one can show that the vector of weights u=(u_1,…,u_n) is Gaussian. Let us write Q_u for its precision matrix, u∼(0,Q_u^-1). Then, from
(<ref>), knowing Q_u determines  u. The value of Q_u is determined in turn by (<ref>), as substituting in  (<ref>) and applying integration by parts gives

    Lu∼(0,C_κ),

where

    L=C_κ+G_H,   [C_κ]_ij=ϕ_i,κ ^2ϕ_j_L^2( ),    [G_H]_ij=∇ϕ_i, H∇ϕ_j_L^2( ) .

From (<ref>) we deduce that the precision Q of the weight vector w is

    Q_u=(C_κ +G_H)^†C_κ^-1 (C_κ +G_H)=C_κ +2G_H+G_HC_κ ^-1G_H,

where we used that since H is symmetric, and we are considering real-valued functions, C_κ,G_H are symmetric.

The finite element basis is chosen so that C_κ, G_H are sparse. Since C_κ is not sparse, it is approximated by the mass lumped version

    [C_κ]_ij:=δ_ij∑_k=1^n [C_κ]_ik.

After this approximation we obtain a sparse precision Q_u which allows us to sample efficiently from u (or rather, its approximation u∼(0,Q_u)) by using a Cholesky decomposition Q_u=LL^T and solving L^Tz=u where z∼(0,𝐈).

It remains to discuss how the integrals defining C_κ,G_H are calculated. Let  be the mesh of the domain. For each triangle T∈ we denote the average value of κ , H on the nodes of T as κ^2(T),H(T) and approximate the integrals in (<ref>) by

    [C_κ]_ij≈∑_T∈→κ ^2(T)∫_Tϕ_i ϕ_j  dx,   
        [G_H]_ij≈∑_T∈→
     k,l=1,2[H(T)]_kl∫_T∂_kϕ_i ∂_lϕ_j  dx.

Finally, we take ϕ_i to be piecewise linear equal to 1 on node i and equal to 0 on every other node. An explicit formula for the integrals in (<ref>) can be found in <cit.>.


 §.§ Comparison of priors for the simulation study

To better understand the four priors π_PC,π_EG,π_U,π_β and their corresponding posteriors, in  <Ref> and <Ref>, we fix σ_u, σ_ε and plot each of the marginal prior and posterior densities mentioned above (up to a multiplicative constant) in log(κ) and v respectively, for a given observation y. As we can see, π_PC and π_EG give similar prior and (as a result) posterior densities, to the point where it is impossible to distinguish them in the plots. The uniform and beta priors and posteriors are similar to each other and differ mainly in the size of support. The beta priors and posteriors resemble a cutoff version of their uniform analogs. This indicates that choosing π_PC and π_EG as priors will give similar results. Likewise, π_U, π_β will also give similar results to each other. This is borne out in the results.






 §.§ Posterior sampling


Returning to the simulation, our parameters are  θ =(logκ,v, log(σ_u),log(σ_ε))∈ℝ^5. Bayes formula gives for any prior π(θ) on θ that

    π(θ  |y)     =π(θ , y)/π(y)π(u | θ ,y)/π(u | θ ,y)=π(θ ) π(u | θ ) π(y | u,θ)/π(y) π(u | θ  ,y).

We note that the right-hand side of (<ref>) involves u. However, the resulting expression is independent of the value of u chosen.
As a result, the log posterior ℓ (θ  |y) is up to an additive constant given by

    ℓ (θ  |y) = ℓ (θ  )+ℓ (u|θ  )+ℓ(y|u,θ)- ℓ (u|θ  ,y)


    The maximum of the right-hand side of (<ref>) is independent of u. However, in the case where ℓ (u|θ,y) were not known, the value of u would affect the result. If u|θ,y were approximated to be Gaussian, it would make the most sense to take u as the mode m_u 	of this Gaussian approximation.

In our case, all terms in (<ref>) are known exactly. The log prior ℓ(θ) is determined by (<ref>) and the remaining terms are given by

    ℓ(u|θ  )            =1/2(logQ_u-u-m_u_Q_u^2-nlog(2π)) 
    ℓ(y|u,θ)       =1/2(logQ_ε -y-Au_Q_ε^2-mlog(2π))
    ℓ(u|θ  ,y)     =1/2(logQ_u|y,θ -u- m_u|y,θ_Q_u|y,θ^2-nlog(2π)),

where we defined x_Q^2:= x^TQx and m is the dimension of y. The precision  Q_u of u is calculated by FEM and is given by  (<ref>) divided by 4πσ_u^2. The remaining precision matrices are

    Q_u|y,θ:= Q_u+ A^T Q_εA,  m_u|y,θ:=Q_u | y,θ^-1(Q_u m_u+A^T Q_εy).

Maximizing the expression in (<ref>) is made possible by the fact that the precision matrices are sparse. We obtain the MAP estimate

    θ:= max_θ∈^3ℓ (θ  |y)=max_θ∈^3ℓ (θ  |y).

Using this, we do the following simulations. We fix as our domain the square [-1,1]^2 from which we sample uniformly and independently m=15 locations x_j_j=1^m and use a mesh with n=2062 degrees of freedom. Then, we do the following.


For π_sim∈{π_PC, π_EG, π_U, π_β}:

    
  * For j= 1,...,N=200:
          
              
  *  We simulate  θ^(j)^true from  π_sim, use the FEM to simulate u^(j) from (<ref>) and then simulate y^(j) from (<ref>).
              
  * For prior π_est∈{π_PC, π_EG, π_U, π_β} on θ:
                    
                        
  * We calculate a maximum a posteriori estimate θ^(j) using (<ref>).
                        
  * We calculate θ^(j)^true_i -θ_i.
                        
  * We calculate a Gaussian approximation to the posterior distribution. To do so, we  write
                              
    Z∼(θ^(j), Σ_θ|y),

                              where the precision, covariance, and marginal standard deviations of Z are, respectively
                              
    [M_θ^(j)|y^(j)]_ij=-∂^2ℓ(θ|y^(j))/∂θ_i∂θ_j,  Σ_θ^(j)|y^(j)=M_θ^(j)|y^(j)^-1.

                              We write π_Z(θ) for the density of Z.
                              If the posterior distribution is Gaussian, then π_Z=π (θ|y^(j)).
                        
  * We approximate the posterior measure using self-normalized importance sampling. This step is necessary as we only have access to the unnormalized ℓ (see (<ref>)).

                              We do this step both with and without Pareto smoothing <cit.>. In both cases, using Z as our proposal distribution. Let {θ^(s)∼Z,s=1,..., S} be i.i.d. We denote the normalized, unnormalized and self-normalized importance ratios by
                              
    r(θ)= π(θ|y^(j))/π_ Z(θ),   r(θ)= exp(ℓ (θ))/π_ Z(θ),   r(θ^(s)):= r (θ^(s))/∑_s=1^S r (θ^(s)).

                              Since we do not have access to the normalizing constant of π (θ|y^(j)), it is the last two ratios that we use. Using r(θ^(s)) calculate the normalized smoothed weights w(θ^(s)) and approximate
                              
    π (θ|y^(j)) ≈π_IS(θ|y^(j)):=  ∑_s=1^S w(θ^(s)) δ_θ^(s)(θ) .

                              In the case where Pareto smoothing is not used, we directly take w(θ^(s))=r(θ^(s)).
                              The reason for using smoothing and no smoothing is to see if Pareto smoothing improves the approximation (for example, by comparing the frequency of times the true parameter is within the 0.95 confidence intervals). With this framework, the expect value of a function f(θ) can be approximated as
                              
    𝔼[f(θ)|y^(j)] ≈𝔼_IS[f(θ)|y^(j)]:= ∑_s=1^S w(θ^(s)) f(θ^(s)).

                              Additionally, we will be interested in the Kullback-Leibler divergence between the true posterior and the Gaussian approximation to the posterior. Here, f(θ)=log(r(θ|y^(j))) is known up to a normalizing constant, but we can approximate it using the self-normalized importance ratios as follows
                              
    KLD(π(·|y^(j))||π_Z )     = ∫_^5 π(θ|y^(j)) log(r (θ))θ∼∑_s=1^S r(θ^(s))  log(Sr (θ^(s))) 
       ∼∑_s=1^S w(θ^(s)) log(Sw (θ^(s))),

                              where the first equality is the definition of the Kullback-Leibler divergence and the first and second approximations are the self-normalized importance sampling step (r ∼ S r ) and smoothing step, respectively. This can be seen to converge to the true Kullback-Leibler divergence as S→∞ by using that r=_Zr r and the central limit theorem.
                        
  * Using π_IS (θ|y^(j)) we approximate the mean, covariance, and KL divergence and build a 0.95 credible interval I for θ. We also check whether the true parameter θ^(j)^true is in I and compare these results against the Gaussian approximation to the posterior in (<ref>).
                              
                        
  * We calculate a_i=ℙ_est[θ _i≤θ_i^(j)^true|y^(j)] for i=1,…,5 both in the case where θ is sampled from π_IS (θ|y^(j)) and from the Gaussian approximation to the posterior using prior π_est. We then plot the empirical cumulative distribution function of a_i and compare it to the uniform distribution. If the implementation is well calibrated, then a_i∼ U(0,1). See <cit.>.
                    
          

The above was tried using the Gaussian distribution Z to replace the posterior and not using importance sampling. The method using no smoothing in the importance sampling was also tested. But both the methods gave worse results than Pareto smoothed importance sampling and we only display this approximation in the results.
Using <Ref>, we can check whether the posterior distribution for θ  is accurate, and using <Ref>, we check whether the implementation is well-calibrated.


    Since the uniform priors give support to extreme values of κ,v, sampling the true parameter θ^true from the uniform prior often leads to singularities in the calculation of the importance samples. These cases made up a high percentage of the simulations (in one case, 90%), for which reason the improper uniform priors were not used to sample from θ^true in the simulations below. Rather, the width of the beta priors was also reduced when simulating these parameters to avoid singularities.

The results for the MAP estimate and the CI lengths were already presented for the anisotropy parameters log(κ),v in <Ref> <Ref> and <Ref>. Here, we begin by presenting these results for the remaining parameters log(σ_u),log(σ_ε).







In <Ref>, we show the frequency of the true parameter being in the 0.95 credible interval. As can be seen in the graphs, the performance depends on the distribution from which θ^true is sampled. If θ^true is sampled from π_PC or π_EG then these two priors perform better, whereas otherwise π_U,π_β perform better. We also observe that the confidence intervals for σ_ε are not wide enough. This may be due to the importance sampling.
The results for the remaining parameters log(σ_u),log(σ_ε) are more similar, which is explained by the fact that the same priors for these parameters are used in all cases. This is a common theme in the results, and we will not comment on it further.

In <Ref>, we show the eCDF of probabilities a_i=ℙ[θ_i≤θ_i^true|y] for each of the five parameters. If the implementation is well calibrated, then a_i should be uniform when π_sim=π. Once more, the performance divides into two cases. If θ^true is sampled from π_PC or π_EG then these two priors perform better, whereas otherwise π_U,π_β perform better.



In  <Ref>, we show the p-values for the KS test between a_i=ℙ_est[θ_i≤θ_i^true|y] and the CDF of the uniform distribution.



In <Ref>, we show the empirical difference against F(a_i) together with a 95% confidence interval. As can be seen from the plots, the implementation is well-calibrated.

In <Ref>, we show the eCDF of the Kullback-Leibler divergence between the posterior and the Gaussian approximation to the posterior around the median. As can be seen from the plots, in all cases π_PC(· |y) and π_EG(· |y) are closer to there respective Gaussian approximations than π_U(· |y) and π_β(· |y). This is because the posterior using π_U,π_β are very much not Gaussian.


To better compare PC and exponential-Gaussian priors, we also plot these separately. As can be seen, once more, they behave quite similarly, showing that both formulations lead to a practically identical penalization of complexity.



§ PRECIPITATION STUDY DETAILS



 §.§ Calculations

In this section, we show how to calculate scores of <Ref>. For what follows, we first define the leave out one mean of y_i having observed all locations except y_i

    m_y_i|y_-i    :=  ∫_^5(  ∫_ y_i π (y_i|y_-i,θ)ỵ_i) π(θ| y) θ = ∫_^5 m_y_i|y_-i,θπ (θ| y) θ.

The expression of m_y_i|y_-i is approximately that of y_i|y_-i, with the exception that the integral in (<ref>) is against π(θ| y) as opposed to π(θ| y_-i). This is a common approach to calculating the LOO predictions as it avoids the need to approximate π (θ| y_-i) for each i. The expression is justified by the fact that the posterior distribution of θ does not depend much on any single observation. That is, π (θ| y)≈π (θ| y_-i). Similarly, we will need the posterior leave out one variance

    σ^2_y_i|y_-i := ∫_^5σ^2_y_i|y_-i,θπ (θ| y) θ,

where σ^2_y_i|y_-i,θ is the posterior variance of y_i having observed all locations except y_i. This will be used later in the computations.

Because of the independence of u, β from θ  we have that

    . u_β| θ   ∼(0,Q_u_β^-1),   where  Q_u_β:=[     Q_β 0_2 × n;        ; 0_n × 2     Q_u ].

We can then calculate θ to maximize the (unnormalized) log-posterior as in (<ref>)

    ℓ (θ   |y) = ℓ (θ  )+ℓ (u_β |θ   )+ℓ(y|u_β ,θ )- ℓ (u_β |θ  ,y).

To calculate the RMSE, CRPS, and DSS, it is sufficient to calculate the posterior mean and variance of y_i given y_-i and θ. We show how to do this efficiently.



Due to the independence of u and ε and of ε _i and ε_-i knowing θ,  we deduce that y_-i = A_-iu+ ε_-i and ε _i are independent knowing θ. As a result,

    y_i | y_-i,θ    = (A_i  u) | y_-i,θ + ε _i | θ∼ (A_i m_u|y_-i,θ, A_i Σ_u|y_-i,θA_i^T + σ_ε^2)
        =  (m_y_i|y_-i,θ, σ_y_i|y_-i,θ^2),

where A_i ∈^1 × N is the i-th row of A. Furthermore, we have that, see (<ref>)

    Q_u|y_-i, θ    = Q_u + σ_ε ^-2A_-i^TA_-i= Q_u|y, θ - σ_ε ^-2A_i^T A_i 
    m_u|y_-i,θ    =σ_ε^-2Σ_u | y_-i,θA_-i^T y_-i,

where we wrote A_-i∈^n-1 × N for the matrix A with the i-th row removed. To avoid calculating an inverse for each i, we can calculate  Σ_u| y, θ once and then use the rank 1 correction provided by the Sherman-Morrison formula to obtain

    Σ_u|y_-i,θ    = Σ_u|y,θ + Σ_u|y,θA_i^TA_iΣ_u|y,θ/σ_ε^2 -A_iΣ_u|y,θA_i^T.


Using (<ref>) in (<ref>) and writing V_i :=A_i Σ_u|y,θA_i^T and q_ε:= σ_ε^-2,
gives

    σ_y_i|y_-i,θ^2     = V_i  + V_i^2/σ_ε^2 - V_i + σ_ε^2   = V_i/1-q_εV_i    +σ_ε^2  =   σ_ε^2 /1-q_εV_i.

Furthermore, for the mean, we have from (<ref>) that

    m_y_i|y_-i,θ    =q_εA_i Σ_u | y_-i,θA_-i^T y_-i= q_εA_i Σ_u | y_-i,θA^T y - q_εA_i Σ_u | y_-i,θA_i^T y_i

Let us write

    η_i := A_i u| y,θ= A_i m_u|y,θ= q_εA_i Σ_u|y,θA^T y,

Then, using the rank 1 correction in (<ref>) we obtain for the first term in (<ref>) that

    q_εA_i Σ_u | y_-i,θA^T y    =  q_εA_i Σ_u | y,θA^T y +q_εA_iΣ_u|y,θA_i^TA_iΣ_u|y,θ/σ_ε^2 -A_iΣ_u|y,θA_i^TA^Ty
        = η_i +  V_i/σ_ε^2 -V_iη_i= η_i /1-q_εV_i.

Whereas, for the second term in (<ref>), we have that

    q_εA_i Σ_u | y_-i,θA_i^T y_i     =  q_εA_i Σ_u | y,θA_i^T y_i +q_εA_iΣ_u|y,θA_i^TA_iΣ_u|y,θ/σ_ε^2 -A_iΣ_u|y,θA_i^TA_i^Ty_i        
        = q_ε V_i y_i + q_εV_i^2/σ_ε^2 - V_iy_i = q_ε  V_i y_i(1+ q_εV_i/1-q_εV_i)= q_ε  V_i y_i/1-q_εV_i

Using (<ref>) and (<ref>) in (<ref>) we obtain

    m_y_i|y_-i,θ    = η_i /1-q_εV_i - q_ε  V_i y_i/1-q_εV_i = η_i - q_ε  V_i y_i/1-q_εV_i= y_i + η_i-y_i/1- q_εV_i.




The term V_i can be calculated efficiently using a Takahashi recursion on the Cholesky factor of the posterior precision Q_u|y,θ without the need to calculate a dense matrix inverse, as implemented by . To calculate η_i, we use a matrix-vector solve.

Next, since the expression of π (θ| y_i) is known up to a normalizing constant, we can use importance sampling and Riemann integration together with (<ref>), (<ref>) to calculate (<ref>).
To estimate the CRPS, we use importance sampling and the expression for the posterior predictive distribution (<ref>). We have

    F_i(t) := ∫_^5(∫_-∞^t π(y_i | y_-i, θ)ỵ_i)p(θ| y)θ≈∑_j=1^J   Φ(t-m_y_i|y_-i,θ_j/σ_y_i|y_-i,θ_j ) w_j,

where Φ is the cumulative distribution function of a standard Gaussian variable, and θ_j,w_j are the importance samples and smoothed self-normalized weights.
Using this expression and the fact that there is an exact expression for the CRPS of a Gaussian mixture <cit.>, we obtain the CRPS.


 §.§ Spatial analysis of the scores

We hypothesize that perhaps some of the models perform better than others in different spatial areas. To confirm this hypothesis, we begin by plotting the difference in the scores between the isotropic and anisotropic PC models. We plot the difference of the scores for the RMSE, CRPS, and DSS in <Ref>.



We also do the same for the difference in the scores between the isotropic and anisotropic EG models in <Ref>. The results are very similar. We observe that the anisotropic models perform better towards the western coast and worse in the remaining region. This is to be expected as the correlation structure is different in the region with high elevation than the lower elevation parts of the domain and indicates that a non-stationary model might be more appropriate as the anisotropy varies in space.




 §.§ Precipitation simulation study

In this section, we repeat the structure of the simulations in <Ref>, but now, to investigate how the quality of the prediction changes with increased data, we use synthetic data. To do so, we simulate precipitation data on 4000 uniformly distributed locations in Norway from the model in (<ref>) with parameters

    ρ    = 201,    v    = (-0.45,0.03),    σ_u     =0.63,    σ_ε    = 0.14,    β_0     = 0.96,    β_1     =0.67  
    ρ    = 132,    v    = (0.00,0.00),     σ_u     =0.65,    σ_ε    = 0.13,    β_0     = 0.93,    β_1     =0.66.

These correspond to data generated from anisotropic and isotropic models, respectively, using the MAP obtained for the precipitation data in <Ref> using anisotropic and isotropic PC priors, respectively. The data is shown in <Ref>.


We then fit the model to the data and calculate the RMSE, CRPS, and DSS scores. We repeat this process for n_y=25,50,100,150,200,400, 600,800,1000 observations sampled uniformly from the simulated data. We repeat this process 100 times for each n_y and each prior. The results are shown in <Ref>. As can be seen from the plots, the anisotropic model performs better than the isotropic model for less data. However, as the number of observations increases, the scores become almost equal. The difference between the anisotropic PC and EG models is less pronounced than in the previous section.


Finally, we also plot the interval score for the parameters (log(κ ), v_1, v_2, log(σ _u)). Given a credible interval (L_F,U_F) with confidence level α , the interval score is defined by

    S_INT (F, y)=U_F-L_F+2/α(L_F-y) 𝕀(y<L_F)+2/α(y-U_F) 𝕀(y>U_F).

The interval score is a proper scoring rule consistent for equal-tail error probability intervals: S(F, G) is minimized for the narrowest P I that has expected coverage 1-α. The results are shown in <Ref>. As can be seen from the plots, the interval score for the parameters (log(κ ),  σ _u) is generally better for the anisotropic model. This is especially pronounced for log(κ ) when there is less data.  The interval score for v_1 and v_2 is not shown for the isotropic model as these parameters are not estimated in the isotropic model.




