spde.matern.operators {rSPDE}

Rational approximations of non-stationary Gaussian SPDE Matern random fields

Description

spde.matern.operators is used for computing a rational SPDE approximation of a Gaussian random fields on R^d defined as a solution to the SPDE

(\kappa(s) - \Delta)^\beta (\tau(s)u(s)) = W.

Usage

spde.matern.operators(
  kappa = NULL,
  tau = NULL,
  theta = NULL,
  B.tau = matrix(c(0, 1, 0), 1, 3),
  B.kappa = matrix(c(0, 0, 1), 1, 3),
  B.sigma = matrix(c(0, 1, 0), 1, 3),
  B.range = matrix(c(0, 0, 1), 1, 3),
  alpha = NULL,
  nu = NULL,
  parameterization = c("spde", "matern"),
  G = NULL,
  C = NULL,
  d = NULL,
  graph = NULL,
  mesh = NULL,
  range_mesh = NULL,
  loc_mesh = NULL,
  m = 1,
  type = c("covariance", "operator"),
  type_rational_approximation = c("chebfun", "brasil", "chebfunLB")
)

Arguments

kappa

Vector with the, possibly spatially varying, range parameter evaluated at the locations of the mesh used for the finite element discretization of the SPDE.

tau

Vector with the, possibly spatially varying, precision parameter evaluated at the locations of the mesh used for the finite element discretization of the SPDE.

theta

Theta parameter that connects B.tau and B.kappa to tau and kappa through a log-linear regression, in case the parameterization is spde, and that connects B.sigma and B.range to tau and kappa in case the parameterization is matern.

B.tau

Matrix with specification of log-linear model for \tau. Will be used if parameterization = 'spde'.

B.kappa

Matrix with specification of log-linear model for \kappa. Will be used if parameterization = 'spde'.

B.sigma

Matrix with specification of log-linear model for \sigma. Will be used if parameterization = 'matern'.

B.range

Matrix with specification of log-linear model for \rho, which is a range-like parameter (it is exactly the range parameter in the stationary case). Will be used if parameterization = 'matern'.

alpha

smoothness parameter. Will be used if the parameterization is 'spde'.

nu

Shape parameter of the covariance function. Will be used if the parameterization is 'matern'.

parameterization

Which parameterization to use? matern uses range, std. deviation and nu (smoothness). spde uses kappa, tau and nu (smoothness). The default is matern.

G

The stiffness matrix of a finite element discretization of the domain of interest.

C

The mass matrix of a finite element discretization of the domain of interest.

d

The dimension of the domain. Does not need to be given if mesh is used.

graph

An optional metric_graph object. Replaces d, C and G.

mesh

An optional inla mesh. d, C and G must be given if mesh is not given.

range_mesh

The range of the mesh. Will be used to provide starting values for the parameters. Will be used if mesh and graph are NULL, and if one of the parameters (kappa or tau for spde parameterization, or sigma or range for matern parameterization) are not provided.

loc_mesh

The mesh locations used to construct the matrices C and G. This option should be provided if one wants to use the rspde_lme() function and will not provide neither graph nor mesh. Only works for 1d data. Does not work for metric graphs. For metric graphs you should supply the graph using the graph argument.

m

The order of the rational approximation, which needs to be a positive integer. The default value is 1.

type

The type of the rational approximation. The options are "covariance" and "operator". The default is "covariance".

type_rational_approximation

Which type of rational approximation should be used? The current types are "chebfun", "brasil" or "chebfunLB".

Details

The approximation is based on a rational approximation of the fractional operator (\kappa(s)^2 -\Delta)^\beta, where \beta = (\nu + d/2)/2. This results in an approximate model on the form

P_l u(s) = P_r W,

where P_j = p_j(L) are non-fractional operators defined in terms of polynomials p_j for j=l,r. The order of p_r is given by m and the order of p_l is m + m_\beta where m_\beta is the integer part of \beta if \beta>1 and m_\beta = 1 otherwise.

The discrete approximation can be written as u = P_r x where x \sim N(0,Q^{-1}) and Q = P_l^T C^{-1} P_l. Note that the matrices P_r and Q may be be ill-conditioned for m>1. In this case, the metehods in operator.operations() should be used for operations involving the matrices, since these methods are more numerically stable.

Value

spde.matern.operators returns an object of class "rSPDEobj. This object contains the quantities listed in the output of fractional.operators() as well as the smoothness parameter \nu.

See Also

fractional.operators(), spde.matern.operators(), matern.operators()

Examples

# Sample non-stationary Matern field on R
tau <- 1
nu <- 0.8

# create mass and stiffness matrices for a FEM discretization
x <- seq(from = 0, to = 1, length.out = 101)
fem <- rSPDE.fem1d(x)

# define a non-stationary range parameter
kappa <- seq(from = 2, to = 20, length.out = length(x))
alpha <- nu + 1/2
# compute rational approximation
op <- spde.matern.operators(
  kappa = kappa, tau = tau, alpha = alpha,
  G = fem$G, C = fem$C, d = 1
)

# sample the field
u <- simulate(op)

# plot the sample
plot(x, u, type = "l", ylab = "u(s)", xlab = "s")