Introduction

The objective of this notebook is the define the multivariate Matern’s covariance models to be implemented in gstlearn. The space is the euclidean model \(\mathbb{R}^d\) with \(d \geq 0\).

Here, we use the Matern’s model defined in Bourotte et al. (2016).

Definition of the models

The univariate model

The univariate Matern’s covariance model on \(\mathbb{R}^d\) is:

\[ \mathcal{M}(h; \nu) = \frac{2^{1-\nu}}{\Gamma(\nu)} (|h|)^{\nu} K_{\nu}(|h|) \] where \(\nu > 0\) is a smoothness parameter, and \(K_{\nu}\) denotes the modified Bessel function of the second kind of order \(\nu\). This is the convention used in gstlearn.

Its spectral density function is: \[ \mathcal{S}(\boldsymbol{\omega}; \nu) = \frac{\Gamma(\nu + \frac{d}{2})}{\Gamma(\nu)\, \pi^{\frac{d}{2}}} \frac{1}{(1+|\boldsymbol{\omega}|^2)^{\nu + \frac{d}{2}}} \]

The space anisotropy and scaling are first applied to the increments to convert them into normalized space distances: \(\boldsymbol{h} \mapsto d_s = \sqrt{\boldsymbol{h}^t Q \boldsymbol{h}}\). The space metric tensor \(Q = R^t D(1/\boldsymbol{a}^2) R\) is defined by the rotation \(R\) and the scaling diagonal matrix \(D(1/\boldsymbol{a}^2)\) defined by its spatial ranges (\(\boldsymbol{a} = (a_i)_{i \in 1:d} > 0\)).

The multivariate model

The multivariate Matern’s model is defined from

  1. The marginal models. Each variable \(i\) is characterized by \(r_{i} > 0\) its relative scale factor, and \(\nu_i > 0\) its smoothness parameter.

  2. A covariance matrix \(\boldsymbol{\Sigma} = [\sigma_{ij}]_{1 \leq i,j \leq p}\).

The family of simple and cross covariance functions is: \[ C_{ij}^{\mathcal{M}}(\boldsymbol{h}) = \sigma_{ij} \times \mathcal{M}(r_{ij} \, d_s(\boldsymbol{h}); \nu_{ij}). \]

with, the cross-parameters are computed from the marginal parameters:

  1. \(\nu_{ij} = (\nu_{i} + \nu_{j})/2\), and
  2. \(r_{ij} = \sqrt{(r_{i}^2 + r_{j}^2)/2}\)

The family of simple and cross spectrum functions is: \[ S_{ij}^{\mathcal{M}}(\boldsymbol{\omega}) = \sigma_{ij} \times \mathcal{S}(\frac{1}{r_{ij}} \, d^s(\boldsymbol{h}); \nu_{ij}) \, \frac{\Pi_{k=1}^d a_k}{r_{ij}^d}. \] where the \(\mathcal{S}\) the spectrum of the standard covariance function and the anisotropic spectral distance \(d^s\) defined by \[ d^s(\boldsymbol{\omega}) = \sqrt{\boldsymbol{\omega}^t R^t D(\boldsymbol{a}^2)R\boldsymbol{\omega} } \]

Note:

  • The covariance matrix must verify the condition that the matrix \([\sigma_{ij}/\tau_{ij}]_{i,j}\) be positive definite, with the maximum correlation matrix defined by \[ \tau_{ij} = \frac{\Gamma(\nu_{ij})}{\Gamma(\nu_{i})^{1/2}\Gamma(\nu_{j})^{1/2}}\frac{r_i^{\nu_{i}}r_j^{\nu_{j}}}{r_{ij}^{2\nu_{ij}}} \]

  • The intrinsic case corresponds to \(\nu_i = \nu\) and \(r_i = r = 1\) (if all variable ranges are equal, they can be taken into account by the metric defined on space-time): \[ C_{ij}(\boldsymbol{h}) = C_{ij}(0) \times \mathcal{M}(d_s; \nu). \]

Initialisation and parameters

Some functions are defined for the QC of the simulations.

The multivariate Matern’model

# parameters for the anisotropy
ani_scales = c(2, 1/2)
ani_angles = c(30, 0)
d = 2
err = gstlearn::defineDefaultSpace(gstlearn::ESpaceType_RN(), d)

# parameters of the variables
nu = c(1/2,  3/4, 2)[1:3]
rr = c(1, 2, 3)[1:3]
stopifnot(length(nu) == length(rr))
p = length(nu)

# parameters for the multivariate correlations
Nu = outer(X = nu, Y = nu, FUN = function(u,v){(u+v)/2})
Rr = outer(X = rr, Y = rr, FUN = function(u,v){sqrt((u^2+v^2)/2)})
Tau = matrix(NaN, p, p)
for (i in 1:p){
  for (j in 1:p) {
    Tau[i,j] = gamma(Nu[i,j])/Rr[i,j]^(2*Nu[i,j])/sqrt(
    (gamma(nu[i])/rr[i]^(2*nu[i]))*(gamma(nu[j])/rr[j]^(2*nu[j])))
  }
}
sigma0 = matrix(rnorm(p*p), p, p)
sigma0 = sigma0 %*% t(sigma0)
sigma  = sigma0 * Tau
stopifnot(all(eigen(sigma)$values >0))

# Multi-variate Matern's models
gsSigma = MatrixSymmetric(p)
for (i in 1:p) {
  for (j in 1:p) {
    gsSigma$setValue(i-1, j-1, sigma[i, j])
  }
}

# the covariance matrix is checked
stopifnot(gsSigma$isDefinitePositive())

# the Matern's covariance in gstlearn
ctxt  = CovContext(nvar = p, ndim = d)
gsMatern = CorMatern_create(
  type = ECov_MATERN(), ctxt = ctxt, 
  params = nu[1:p],
  kappas = rr[1:p],
  ranges = ani_scales[1:d] , 
  angles = ani_angles[1:d], 
  flagRange = FALSE)

# the model is defined by the covariance
gsMod = ModelGeneric(ctxt = ctxt) 
err   = gsMod$setCov(cova = gsMatern)

Covariance and spectrum

Covariances of the multi-variate model

Spectums of the multi-variate model

Simulations

err = grid$deleteColumns(names = paste(prefix, "*", sep = "."))

# simulation parameters
seed = 258425
set.seed(seed)
ns = 5000 # number of spectral waves
nbsimu = 3

err = simuSpectral(dbin = NULL, dbout = grid, cova = gsMatern, nbsimu = nbsimu, 
                   seed = seed, ns = ns, verbose = FALSE, 
                   namconv = NamingConvention(prefix))

for (isim in 1:nbsimu) {
  nm_var = paste(prefix, paste0("V", 1:p), paste0("S", isim), sep = ".")
  knitr::kable(dbStatisticsMono(db = grid, names = nm_var, opers = opers)$toTL(), digits = 4)
  plot_hist_map_vario(grid = grid, names = nm_var, model = gsMatern,
                      main = paste0("Matern: d = ", gsMatern$getNDim(),
                                    ", p = ", gsMatern$getNVar()),
                      flag.hist = TRUE, flag.map = TRUE, flag.vario = TRUE)
}

References