Package 'rSPDE'

Description

rSPDE is used for approximating fractional elliptic SPDEs

$L^\beta (\tau u(s)) = W,$

where $L$ is a differential operator and $\beta>0$ is a general fractional power.

Details

The approximation is based on a rational approximation of the fractional operator, and allows for computationally efficient inference and simulation.

The main functions for computing rational approximation objects are:

fractional.operators()

works for general rational operators

matern.operators()

works for random fields with stationary Matern covariance functions

spde.matern.operators()

works for random fields with defined as solutions to a possibly non-stationary Matern-type SPDE model.

rspde.matern()

R-INLA implementation of the covariance-based rational approximation for random fields with stationary Matern covariance functions

Basic statistical operations such as likelihood evaluations (see ⁠[rSPDE.loglike], [rSPDE.matern.loglike]⁠) and kriging predictions (see ⁠[predict.rSPDEobj], [predict.CBrSPDEobj]⁠) using the rational approximations are also implemented.

For illustration purposes, the package contains a simple FEM implementation for models on R. For spatial models, the FEM implementation in the R-INLA package is recommended.

For a more detailed introduction to the package, see the rSPDE Vignettes.

Author(s)

Maintainer: David Bolin [email protected]

Authors:

• Alexandre Simas [email protected]

Other contributors:

• Finn Lindgren [email protected] [contributor]

See Also

Useful links:

• https://davidbolin.github.io/rSPDE/

• Report bugs at https://github.com/davidbolin/rSPDE/issues

Description

Augment accepts a model object and a dataset and adds information about each observation in the dataset. It includes predicted values in the .fitted column, residuals in the .resid column, and standard errors for the fitted values in a .se.fit column. It also contains the New columns always begin with a . prefix to avoid overwriting columns in the original dataset.

Usage

## S3 method for class 'rspde_lme'
augment(
  x,
  newdata = NULL,
  loc = NULL,
  mesh = FALSE,
  which_repl = NULL,
  se_fit = FALSE,
  conf_int = FALSE,
  pred_int = FALSE,
  level = 0.95,
  n_samples = 100,
  ...
)

Arguments

Value

A tidyr::tibble() with columns:

• .fitted Fitted or predicted value.

• .fittedlwrconf Lower bound of the confidence interval, if conf_int = TRUE

• .fitteduprconf Upper bound of the confidence interval, if conf_int = TRUE

• .fittedlwrpred Lower bound of the prediction interval, if pred_int = TRUE

• .fitteduprpred Upper bound of the prediction interval, if pred_int = TRUE

• .fixed Prediction of the fixed effects.

• .random Prediction of the random effects.

• .resid The ordinary residuals, that is, the difference between observed and fitted values.

• .se_fit Standard errors of fitted values, if se_fit = TRUE.

See Also

glance.rspde_lme

Description

rSPDE inlabru mapper

Usage

## S3 method for class 'inla_rspde'
bru_get_mapper(model, ...)

## S3 method for class 'bru_mapper_inla_rspde'
ibm_n(mapper, ...)

## S3 method for class 'bru_mapper_inla_rspde'
ibm_values(mapper, ...)

## S3 method for class 'bru_mapper_inla_rspde'
ibm_jacobian(mapper, input, ...)

Arguments

Examples

# devel version
if (requireNamespace("INLA", quietly = TRUE) &&
  requireNamespace("inlabru", quietly = TRUE)) {
  library(INLA)
  library(inlabru)

  set.seed(123)
  m <- 100
  loc_2d_mesh <- matrix(runif(m * 2), m, 2)
  mesh_2d <- inla.mesh.2d(
    loc = loc_2d_mesh,
    cutoff = 0.05,
    max.edge = c(0.1, 0.5)
  )
  sigma <- 1
  range <- 0.2
  nu <- 0.8
  kappa <- sqrt(8 * nu) / range
  op <- matern.operators(
    mesh = mesh_2d, nu = nu,
    range = range, sigma = sigma, m = 2,
    parameterization = "matern"
  )
  u <- simulate(op)
  A <- inla.spde.make.A(
    mesh = mesh_2d,
    loc = loc_2d_mesh
  )
  sigma.e <- 0.1
  y <- A %*% u + rnorm(m) * sigma.e
  y <- as.vector(y)

  data_df <- data.frame(
    y = y, x1 = loc_2d_mesh[, 1],
    x2 = loc_2d_mesh[, 2]
  )
  rspde_model <- rspde.matern(
    mesh = mesh_2d,
    nu_upper_bound = 2
  )

  cmp <- y ~ Intercept(1) +
    field(cbind(x1,x2), model = rspde_model)


  rspde_fit <- bru(cmp, data = data_df)
  summary(rspde_fit)
}
# devel.tag

Description

rSPDE anisotropic inlabru mapper

Usage

## S3 method for class 'inla_rspde_anisotropic2d'
bru_get_mapper(model, ...)

Arguments

Description

rSPDE space time inlabru mapper

Usage

## S3 method for class 'inla_rspde_spacetime'
bru_get_mapper(model, ...)

Arguments

Description

This function evaluates the log-likelihood function for observations of a non-stationary Gaussian process defined as the solution to the SPDE

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

The observations are assumed to be generated as $Y_i = u(s_i) + \epsilon_i$, where $\epsilon_i$ are iid mean-zero Gaussian variables. The latent model is approximated using a rational approximation of the fractional SPDE model.

Usage

construct.spde.matern.loglike(
  object,
  Y,
  A,
  sigma.e = NULL,
  mu = 0,
  user_nu = NULL,
  user_m = NULL,
  log_scale = TRUE,
  return_negative_likelihood = TRUE
)

Arguments

Value

The log-likelihood function. The parameters of the returned function are given in the order theta, nu, sigma.e, whenever they are available.

See Also

matern.operators(), predict.CBrSPDEobj()

Examples

# this example illustrates how the function can be used for maximum
# likelihood estimation
# Sample a Gaussian Matern process on R using a rational approximation
set.seed(123)
sigma.e <- 0.1
n.rep <- 10
n.obs <- 100
n.x <- 51
# create mass and stiffness matrices for a FEM discretization
x <- seq(from = 0, to = 1, length.out = n.x)
fem <- rSPDE.fem1d(x)
tau <- rep(0.5, n.x)
nu <- 0.8
alpha <- nu + 0.5
kappa <- rep(1, n.x)
# Matern parameterization
# compute rational approximation
op <- spde.matern.operators(
  loc_mesh = x,
  kappa = kappa, tau = tau, alpha = alpha,
  parameterization = "spde", d = 1
)
# Sample the model
u <- simulate(op, n.rep)
# Create some data
obs.loc <- runif(n = n.obs, min = 0, max = 1)
A <- rSPDE.A1d(x, obs.loc)
noise <- rnorm(n.obs * n.rep)
dim(noise) <- c(n.obs, n.rep)
Y <- as.matrix(A %*% u + sigma.e * noise)
# define negative likelihood function for optimization using matern.loglike
mlik <- construct.spde.matern.loglike(op, Y, A)
#' #The parameters can now be estimated by minimizing mlik with optim

# Choose some reasonable starting values depending on the size of the domain
theta0 <- log(c(1 / sqrt(var(c(Y))), sqrt(8), 0.9, 0.01))
# run estimation and display the results
theta <- optim(theta0, mlik)
print(data.frame(
  tau = c(tau[1], exp(theta$par[1])), kappa = c(kappa[1], exp(theta$par[2])),
  nu = c(nu, exp(theta$par[3])), sigma.e = c(sigma.e, exp(theta$par[4])),
  row.names = c("Truth", "Estimates")
))


# SPDE parameterization
# compute rational approximation
op <- spde.matern.operators(
  kappa = kappa, tau = tau, alpha = alpha,
  loc_mesh = x, d = 1,
  parameterization = "spde"
)
# Sample the model
u <- simulate(op, n.rep)
# Create some data
obs.loc <- runif(n = n.obs, min = 0, max = 1)
A <- rSPDE.A1d(x, obs.loc)
noise <- rnorm(n.obs * n.rep)
dim(noise) <- c(n.obs, n.rep)
Y <- as.matrix(A %*% u + sigma.e * noise)
# define negative likelihood function for optimization using matern.loglike
mlik <- construct.spde.matern.loglike(op, Y, A)
#' #The parameters can now be estimated by minimizing mlik with optim

# Choose some reasonable starting values depending on the size of the domain
theta0 <- log(c(1 / sqrt(var(c(Y))), sqrt(8), 0.9, 0.01))
# run estimation and display the results
theta <- optim(theta0, mlik)
print(data.frame(
  tau = c(tau[1], exp(theta$par[1])), kappa = c(kappa[1], exp(theta$par[2])),
  nu = c(nu, exp(theta$par[3])), sigma.e = c(sigma.e, exp(theta$par[4])),
  row.names = c("Truth", "Estimates")
))

Description

Train and test splits

Usage

create_train_test_indices(
  data,
  cv_type = c("k-fold", "loo", "lpo"),
  k = 5,
  percentage = 20,
  number_folds = 10
)

Arguments

Value

A list with two elements, train containing the training indices and test containing indices.

Description

Obtain several scores for a list of fitted models according to a folding scheme.

Usage

cross_validation(
  models,
  model_names = NULL,
  scores = c("mse", "crps", "scrps", "dss"),
  cv_type = c("k-fold", "loo", "lpo"),
  k = 5,
  percentage = 20,
  number_folds = 10,
  n_samples = 1000,
  return_scores_folds = FALSE,
  orientation_results = c("negative", "positive"),
  include_best = TRUE,
  train_test_indexes = NULL,
  return_train_test = FALSE,
  return_post_samples = FALSE,
  return_true_test_values = FALSE,
  parallelize_RP = FALSE,
  n_cores_RP = parallel::detectCores() - 1,
  true_CV = TRUE,
  save_settings = FALSE,
  print = TRUE,
  fit_verbose = FALSE
)

Arguments

Value

A data.frame with the fitted models and the corresponding scores.

Description

folded.matern.covariance.1d evaluates the 1d folded Matern covariance function over an interval $[0,L]$.

Usage

folded.matern.covariance.1d(
  h,
  m,
  kappa,
  nu,
  sigma,
  L = 1,
  N = 10,
  boundary = c("neumann", "dirichlet", "periodic")
)

Arguments

Details

folded.matern.covariance.1d evaluates the 1d folded Matern covariance function over an interval $[0,L]$ under different boundary conditions. For periodic boundary conditions

$C_{\mathcal{P}}(h,m) = \sum_{k=-\infty}^{\infty} (C(h-m+2kL),$

for Neumann boundary conditions

$C_{\mathcal{N}}(h,m) = \sum_{k=-\infty}^{\infty} (C(h-m+2kL)+C(h+m+2kL)),$

and for Dirichlet boundary conditions:

$C_{\mathcal{D}}(h,m) = \sum_{k=-\infty}^{\infty} (C(h-m+2kL)-C(h+m+2kL)),$

where $C(\cdot)$ is the Matern covariance function:

$C(h) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)}(\kappa h)^\nu K_\nu(\kappa h).$

We consider the truncation:

$C_{{\mathcal{P}},N}(h,m) = \sum_{k=-N}^{N} C(h-m+2kL), C_{\mathcal{N},N}(h,m) = \sum_{k=-\infty}^{\infty} (C(h-m+2kL)+C(h+m+2kL)),$

and

$C_{\mathcal{D},N}(h,m) = \sum_{k=-N}^{N} (C(h-m+2kL)-C(h+m+2kL)).$

Value

A matrix with the corresponding covariance values.

Examples

x <- seq(from = 0, to = 1, length.out = 101)
plot(x, folded.matern.covariance.1d(rep(0.5, length(x)), x,
  kappa = 10, nu = 1 / 5, sigma = 1
),
type = "l", ylab = "C(h)", xlab = "h"
)

Description

folded.matern.covariance.2d evaluates the 2d folded Matern covariance function over an interval $[0,L]\times [0,L]$.

Usage

folded.matern.covariance.2d(
  h,
  m,
  kappa,
  nu,
  sigma,
  L = 1,
  N = 10,
  boundary = c("neumann", "dirichlet", "periodic", "R2")
)

Arguments

Details

folded.matern.covariance.2d evaluates the 1d folded Matern covariance function over an interval $[0,L]\times [0,L]$ under different boundary conditions. For periodic boundary conditions

$C_{\mathcal{P}}((h_1,h_2),(m_1,m_2)) = \sum_{k_2=-\infty}^\infty \sum_{k_1=-\infty}^{\infty} (C(\|(h_1-m_1+2k_1L,h_2-m_2+2k_2L)\|),$

for Neumann boundary conditions

$C_{\mathcal{N}}((h_1,h_2),(m_1,m_2)) = \sum_{k_2=-\infty}^\infty \sum_{k_1=-\infty}^{\infty} (C(\|(h_1-m_1+2k_1L,h_2-m_2+2k_2L)\|)+C(\|(h_1-m_1+2k_1L, h_2+m_2+2k_2L)\|)+C(\|(h_1+m_1+2k_1L,h_2-m_2+2k_2L)\|)+ C(\|(h_1+m_1+2k_1L,h_2+m_2+2k_2L)\|)),$

and for Dirichlet boundary conditions:

$C_{\mathcal{D}}((h_1,h_2),(m_1,m_2)) = \sum_{k_2=-\infty}^\infty \sum_{k_1=-\infty}^{\infty} (C(\|(h_1-m_1+2k_1L,h_2-m_2+2k_2L)\|)- C(\|(h_1-m_1+2k_1L,h_2+m_2+2k_2L)\|)-C(\|(h_1+m_1+2k_1L, h_2-m_2+2k_2L)\|)+C(\|(h_1+m_1+2k_1L,h_2+m_2+2k_2L)\|)),$

where $C(\cdot)$ is the Matern covariance function:

$C(h) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)} (\kappa h)^\nu K_\nu(\kappa h).$

We consider the truncation for $k_1,k_2$ from $-N$ to $N$.

Value

The correspoding covariance.

Examples

h <- c(0.5, 0.5)
m <- c(0.5, 0.5)
folded.matern.covariance.2d(h, m, kappa = 10, nu = 1 / 5, sigma = 1)

Description

fractional.operators is used for computing an approximation, which can be used for inference and simulation, of the fractional SPDE

$L^\beta (\tau u(s)) = W.$

Here $L$ is a differential operator, $\beta>0$ is the fractional power, $\tau$ is a positive scalar or vector that scales the variance of the solution $u$, and $W$ is white noise.

Usage

fractional.operators(L, beta, C, scale.factor, m = 1, tau = 1)

Arguments

Details

The approximation is based on a rational approximation of the fractional operator, resulting 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 methods in operator.operations() should be used for operations involving the matrices, since these methods are more numerically stable.

Value

fractional.operators returns an object of class "rSPDEobj". This object contains the following quantities:

See Also

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

Examples

# Compute rational approximation of a Gaussian process with a
# Matern covariance function on R
kappa <- 10
sigma <- 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)

# compute rational approximation of covariance function at 0.5
tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
  (4 * pi)^(1 / 2) * gamma(nu + 1 / 2)))
op <- fractional.operators(
  L = fem$G + kappa^2 * fem$C, beta = (nu + 1 / 2) / 2,
  C = fem$C, scale.factor = kappa^2, tau = tau
)

v <- t(rSPDE.A1d(x, 0.5))
c.approx <- Sigma.mult(op, v)

# plot the result and compare with the true Matern covariance
plot(x, matern.covariance(abs(x - 0.5), kappa, nu, sigma),
  type = "l", ylab = "C(h)",
  xlab = "h", main = "Matern covariance and rational approximations"
)
lines(x, c.approx, col = 2)

Description

Auxiliar function to obtain domain-based initial values for log-likelihood optimization in rSPDE models with a latent stationary Gaussian Matern model

Usage

get.initial.values.rSPDE(
  mesh = NULL,
  mesh.range = NULL,
  graph.obj = NULL,
  n.spde = 1,
  dim = NULL,
  B.tau = NULL,
  B.kappa = NULL,
  B.sigma = NULL,
  B.range = NULL,
  nu = NULL,
  parameterization = c("matern", "spde"),
  include.nu = TRUE,
  log.scale = TRUE,
  nu.upper.bound = NULL
)

Arguments

Value

A vector of the form (theta_1,theta_2,theta_3) or where theta_1 is the initial guess for tau, theta_2 is the initial guess for kappa and theta_3 is the initial guess for nu.

Description

Data frame for result objects from R-INLA fitted models to be used in ggplot2

Usage

gg_df(result, ...)

Arguments

Value

A data frame containing the posterior densities.

Description

Returns a ggplot-friendly data-frame with the marginal posterior densities.

Usage

## S3 method for class 'rspde_result'
gg_df(
  result,
  parameter = result$params,
  transform = TRUE,
  restrict_x_axis = NULL,
  restrict_quantiles = NULL,
  ...
)

Arguments

Value

A data frame containing the posterior densities.

Description

Glance accepts a rspde_lme object and returns a tidyr::tibble() with exactly one row of model summaries. The summaries are the square root of the estimated variance of the measurement error, residual degrees of freedom, AIC, BIC, log-likelihood, the type of latent model used in the fit and the total number of observations.

Usage

## S3 method for class 'rspde_lme'
glance(x, ...)

Arguments

Value

A tidyr::tibble() with exactly one row and columns:

• nobs Number of observations used.

• sigma the square root of the estimated residual variance

• logLik The log-likelihood of the model.

• AIC Akaike's Information Criterion for the model.

• BIC Bayesian Information Criterion for the model.

• deviance Deviance of the model.

• df.residual Residual degrees of freedom.

• model.type Type of latent model fitted.

See Also

augment.rspde_lme

Description

Extracts data from metric graphs to be used by 'INLA' and 'inlabru'.

Usage

graph_data_rspde(
  graph_rspde,
  name = "field",
  repl = NULL,
  repl_col = NULL,
  group = NULL,
  group_col = NULL,
  only_pred = FALSE,
  loc = NULL,
  bru = FALSE,
  tibble = FALSE,
  drop_na = FALSE,
  drop_all_na = TRUE
)

Arguments

Value

An 'INLA' and 'inlabru' friendly list with the data.

Description

Compute prediction of a formula-based expression on a testing set based on a training set.

Usage

group_predict(
  models,
  model_names = NULL,
  formula = NULL,
  train_indices,
  test_indices,
  n_samples = 1000,
  pseudo_predict = TRUE,
  return_samples = FALSE,
  return_hyper_samples = FALSE,
  n_hyper_samples = 1,
  compute_posterior_means = TRUE,
  print = TRUE,
  fit_verbose = FALSE
)

Arguments

Value

A data.frame with the fitted models and the corresponding scores.

Description

intrinsic.matern.operators is used for computing a covariance-based rational SPDE approximation of intrinsic fields on $R^d$ defined through the SPDE

$(-\Delta)^{\beta/2}(\kappa^2-\Delta)^{\alpha/2} (\tau u) = \mathcal{W}$

Usage

intrinsic.matern.operators(
  kappa,
  tau,
  alpha,
  beta = 1,
  G = NULL,
  C = NULL,
  d = NULL,
  mesh = NULL,
  graph = NULL,
  loc_mesh = NULL,
  m_alpha = 2,
  m_beta = 2,
  compute_higher_order = FALSE,
  return_block_list = FALSE,
  type_rational_approximation = c("chebfun", "brasil", "chebfunLB"),
  fem_mesh_matrices = NULL,
  scaling = NULL
)

Arguments

Details

The covariance operator

$\tau^{-2}(-\Delta)^{\beta}(\kappa^2-\Delta)^{\alpha}$

is approximated based on rational approximations of the two fractional components. The Laplacians are equipped with homogeneous Neumann boundary conditions and a zero-mean constraint is additionally imposed to obtained a non-intrinsic model.

Value

intrinsic.matern.operators returns an object of class "intrinsicCBrSPDEobj". This object is a list containing the following quantities:

Examples

if (requireNamespace("RSpectra", quietly = TRUE)) {
  x <- seq(from = 0, to = 10, length.out = 201)
  beta <- 1
  alpha <- 1
  kappa <- 1
  op <- intrinsic.matern.operators(
    kappa = kappa, tau = 1, alpha = alpha,
    beta = beta, loc_mesh = x, d = 1
  )
  # Compute and plot the variogram of the model
  Sigma <- op$A[,-1] %*% solve(op$Q[-1,-1], t(op$A[,-1]))
  One <- rep(1, times = ncol(Sigma))
  D <- diag(Sigma)
  Gamma <- 0.5 * (One %*% t(D) + D %*% t(One) - 2 * Sigma)
  k <- 100
  plot(x, Gamma[k, ], type = "l")
  lines(x,
    variogram.intrinsic.spde(x[k], x, kappa, alpha, beta, L = 10, d = 1),
    col = 2, lty = 2
  )
}

Description

matern.covariance evaluates the Matern covariance function

$C(h) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)}(\kappa h)^\nu K_\nu(\kappa h).$

Usage

matern.covariance(h, kappa, nu, sigma)

Arguments

Value

A vector with the values C(h).

Examples

x <- seq(from = 0, to = 1, length.out = 101)
plot(x, matern.covariance(abs(x - 0.5), kappa = 10, nu = 1 / 5, sigma = 1),
  type = "l", ylab = "C(h)", xlab = "h"
)

Description

matern.operators is used for computing a rational SPDE approximation of a stationary Gaussian random fields on $R^d$ with a Matern covariance function

$C(h) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)} (\kappa h)^\nu K_\nu(\kappa h)$

Usage

matern.operators(
  kappa = NULL,
  tau = NULL,
  alpha = NULL,
  sigma = NULL,
  range = NULL,
  nu = NULL,
  G = NULL,
  C = NULL,
  d = NULL,
  mesh = NULL,
  graph = NULL,
  range_mesh = NULL,
  loc_mesh = NULL,
  m = 1,
  type = c("covariance", "operator"),
  parameterization = c("spde", "matern"),
  compute_higher_order = FALSE,
  return_block_list = FALSE,
  type_rational_approximation = c("chebfun", "brasil", "chebfunLB"),
  compute_logdet = FALSE
)

Arguments

Details

If type is "covariance", we use the covariance-based rational approximation of the fractional operator. In the SPDE approach, we model $u$ as the solution of the following SPDE:

$L^{\alpha/2}(\tau u) = \mathcal{W},$

where $L = -\Delta +\kappa^2 I$ and $\mathcal{W}$ is the standard Gaussian white noise. The covariance operator of $u$ is given by $L^{-\alpha}$. Now, let $L_h$ be a finite-element approximation of $L$. We can use a rational approximation of order $m$ on $L_h^{-\alpha}$ to obtain the following approximation:

$L_{h,m}^{-\alpha} = L_h^{-m_\alpha} p(L_h^{-1})q(L_h^{-1})^{-1},$

where $m_\alpha = \lfloor \alpha\rfloor$, $p$ and $q$ are polynomials arising from such rational approximation. From this approximation we construct an approximate precision matrix for $u$.

If type is "operator", the approximation is based on a rational approximation of the fractional operator $(\kappa^2 -\Delta)^\beta$, where $\beta = (\nu + d/2)/2$. This results in an approximate model of 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 methods in operator.operations() should be used for operations involving the matrices, since these methods are more numerically stable.

Value

If type is "covariance", then matern.operators returns an object of class "CBrSPDEobj". This object is a list containing the following quantities:

If type is "operator", then matern.operators returns an object of class "rSPDEobj". This object contains the quantities listed in the output of fractional.operators(), the G matrix, the dimension of the domain, as well as the parameters of the covariance function.

See Also

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

Examples

# Compute the covariance-based rational approximation of a
# Gaussian process with a Matern covariance function on R
kappa <- 10
sigma <- 1
nu <- 0.8
range <- sqrt(8 * nu) / kappa

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

# compute rational approximation of covariance function at 0.5
op_cov <- matern.operators(
  loc_mesh = x, nu = nu,
  range = range, sigma = sigma, d = 1, m = 2,
  parameterization = "matern"
)

v <- t(rSPDE.A1d(x, 0.5))
# Compute the precision matrix
Q <- op_cov$Q
# A matrix here is the identity matrix
A <- Diagonal(nobs)
# We need to concatenate 3 A's since we are doing a covariance-based rational
# approximation of order 2
Abar <- cbind(A, A, A)
w <- rbind(v, v, v)
# The approximate covariance function:
c_cov.approx <- (Abar) %*% solve(Q, w)
c.true <- folded.matern.covariance.1d(
  rep(0.5, length(x)),
  abs(x), kappa, nu, sigma
)

# plot the result and compare with the true Matern covariance
plot(x, c.true,
  type = "l", ylab = "C(h)",
  xlab = "h", main = "Matern covariance and rational approximations"
)
lines(x, c_cov.approx, col = 2)


# Compute the operator-based rational approximation of a Gaussian
# process with a Matern covariance function on R
kappa <- 10
sigma <- 1
nu <- 0.8
range <- sqrt(8 * nu) / kappa

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

# compute rational approximation of covariance function at 0.5
op <- matern.operators(
  range = range, sigma = sigma, nu = nu,
  loc_mesh = x, d = 1,
  type = "operator",
  parameterization = "matern"
)

v <- t(rSPDE.A1d(x, 0.5))
c.approx <- Sigma.mult(op, v)
c.true <- folded.matern.covariance.1d(
  rep(0.5, length(x)),
  abs(x), kappa, nu, sigma
)

# plot the result and compare with the true Matern covariance
plot(x, c.true,
  type = "l", ylab = "C(h)",
  xlab = "h", main = "Matern covariance and rational approximation"
)
lines(x, c.approx, col = 2)

Description

The function is used for computing an approximation, which can be used for inference and simulation, of the fractional SPDE

$(\kappa^2 - \Delta)^{\alpha/2} (\tau u(s)) = W$

on intervals or metric graphs. Here $W$ is Gaussian white noise, $\kappa$ controls the range, $\alpha = \nu + 1/2$ with $\nu>0$ controls the smoothness and $\tau$ is related to the marginal variances through

$\sigma^2 = \frac{\Gamma(\nu)}{\tau^2\Gamma(\alpha)2\sqrt{\pi}\kappa^{2\nu}}.$

Usage

matern.rational(
  graph = NULL,
  loc = NULL,
  bc = c("free", "Neumann", "Dirichlet"),
  kappa = NULL,
  range = NULL,
  nu = NULL,
  sigma = NULL,
  tau = NULL,
  alpha = NULL,
  m = 2,
  parameterization = c("matern", "spde"),
  type_rational_approximation = "brasil",
  type_interp = "spline"
)

Arguments

Value

A model object for the the approximation

Examples

s <- seq(from = 0, to = 1, length.out = 101)
kappa <- 20
sigma <- 2
nu <- 0.8
r <- sqrt(8*nu)/kappa #range parameter
op_cov <- matern.rational(loc = s, nu = nu, range = r, sigma = sigma, m = 2, 
parameterization = "matern")
cov.true <- matern.covariance(abs(s-s[1]), kappa = kappa, sigma = sigma, nu = nu)
cov.approx <- op_cov$covariance(ind = 1)

plot(s, cov.true)
lines(s, cov.approx, col = 2)

Description

Computes a rational approximation of the Matern covariance function on intervals.

Usage

matern.rational.cov(
  h,
  order,
  kappa,
  nu,
  sigma,
  type_rational = "brasil",
  type_interp = "linear"
)

Arguments

Value

The covariance matrix of the approximation

Examples

h <- seq(from = 0, to = 1, length.out = 100)
cov.true <- matern.covariance(h, kappa = 10, sigma = 1, nu = 0.8)
cov.approx <- matern.rational.cov(h, kappa = 10, sigma = 1, nu = 0.8, order = 2)

plot(h, cov.true)
lines(h, cov.approx, col = 2)

Description

matern2d.operators is used for computing a rational SPDE approximation of a stationary Gaussian random fields on $R^d$ with a Matern covariance function

$C(h) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)}(\sqrt{h^T H^{-1}h})^\nu K_\nu(\sqrt{h^T H^{-1}h})$

, based on a SPDE representation of the form

$(I - \nabla\cdot(H\nabla))^{(\nu+1)/2}u = c\sigma W$

, where $c>0$ is a constant. The matrix $H$ is defined as

$\begin{bmatrix} h_x^2 & h_xh_yh_{xy} \\ h_xh_yh_{xy} & h_y^2 \end{bmatrix}$

Usage

matern2d.operators(
  hx = NULL,
  hy = NULL,
  hxy = NULL,
  nu = NULL,
  sigma = NULL,
  mesh = NULL,
  fem = NULL,
  m = 1,
  type_rational_approximation = c("chebfun", "brasil", "chebfunLB"),
  return_fem_matrices = FALSE
)

Arguments

Value

An object of type CBrSPDEobj2d

See Also

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

Examples

library(fmesher)
n_loc <- 2000
loc_2d_mesh <- matrix(runif(n_loc * 2), n_loc, 2)
mesh_2d <- fm_mesh_2d(loc = loc_2d_mesh, cutoff = 0.03, max.edge = c(0.1, 0.5))
op <- matern2d.operators(mesh = mesh_2d)

Description

Functions for multiplying and solving with the $P_r$ and $P_l$ operators as well as the latent precision matrix $Q = P_l C^{-1}P_l$ and covariance matrix $\Sigma = P_r Q^{-1} P_r^T$. These operations are done without first assembling $P_r$, $P_l$ in order to avoid numerical problems caused by ill-conditioned matrices.

Usage

Pr.mult(obj, v, transpose = FALSE)

Pr.solve(obj, v, transpose = FALSE)

Pl.mult(obj, v, transpose = FALSE)

Pl.solve(obj, v, transpose = FALSE)

Q.mult(obj, v)

Q.solve(obj, v)

Qsqrt.mult(obj, v, transpose = FALSE)

Qsqrt.solve(obj, v, transpose = FALSE)

Sigma.mult(obj, v)

Sigma.solve(obj, v)

Arguments

Details

Pl.mult, Pr.mult, and Q.mult multiplies the vector with the respective object. Changing mult to solve in the function names multiplies the vector with the inverse of the object. Qsqrt.mult and Qsqrt.solve performs the operations with the square-root type object $Q_r = C^{-1/2}P_l$ defined so that $Q = Q_r^T Q_r$.

Value

A vector with the values of the operation

Description

Function to get the precision matrix of a CBrSPDEobj object

Usage

precision(object, ...)

## S3 method for class 'CBrSPDEobj'
precision(
  object,
  user_nu = NULL,
  user_kappa = NULL,
  user_sigma = NULL,
  user_range = NULL,
  user_tau = NULL,
  user_m = NULL,
  ...
)

Arguments

Value

The precision matrix.

See Also

simulate.CBrSPDEobj(), matern.operators()

Examples

# Compute the covariance-based rational approximation of a
# Gaussian process with a Matern covariance function on R
kappa <- 10
sigma <- 1
nu <- 0.8
range <- 0.2

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

# compute rational approximation of covariance function at 0.5
tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
  (4 * pi)^(1 / 2) * gamma(nu + 1 / 2)))
op_cov <- matern.operators(
  loc_mesh = x, nu = nu,
  range = range, sigma = sigma, d = 1, m = 2,
  parameterization = "matern"
)

# Get the precision matrix:
prec_matrix <- precision(op_cov)

Description

Function to get the precision matrix of a CBrSPDEobj2d object

Usage

## S3 method for class 'CBrSPDEobj2d'
precision(
  object,
  user_nu = NULL,
  user_hx = NULL,
  user_hy = NULL,
  user_hxy = NULL,
  user_sigma = NULL,
  user_m = NULL,
  ...
)

Arguments

Value

The precision matrix.

See Also

simulate.CBrSPDEobj2d(), matern2d.operators()

Examples

library(fmesher)
n_loc <- 2000
loc_2d_mesh <- matrix(runif(n_loc * 2), n_loc, 2)
mesh_2d <- fm_mesh_2d(loc = loc_2d_mesh, cutoff = 0.03, max.edge = c(0.1, 0.5))
op <- matern2d.operators(mesh = mesh_2d)
Q <- precision(op)

Description

Function to get the precision matrix of an inla_rspde object created with the rspde.matern() function.

Usage

## S3 method for class 'inla_rspde'
precision(object, theta = NULL, ...)

Arguments

Value

The precision matrix.

See Also

precision.CBrSPDEobj(), matern.operators()

Description

Function to get the precision matrix of a rSPDEobj1d object

Usage

## S3 method for class 'rSPDEobj1d'
precision(
  object,
  user_loc = NULL,
  user_nu = NULL,
  user_kappa = NULL,
  user_sigma = NULL,
  user_range = NULL,
  user_tau = NULL,
  user_m = NULL,
  ordering = c("field", "location"),
  ldl = FALSE,
  ...
)

Arguments

Value

A list containing the precision matrix Q of the process and its derivatives if they exist, and a matrix A that extracts the elements corresponding to the process. If ldl=TRUE, the LDL factorization is returned instead of Q. If the locations are not ordered, the precision matrix is given for the ordered locations, but the A matrix returns to the original order.

See Also

simulate.rSPDEobj1d(), matern.rational()

Examples

# Compute the covariance-based rational approximation of a
# Gaussian process with a Matern covariance function on R
sigma <- 1
nu <- 0.8
range <- 0.2

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

op_cov <- matern.rational(
  loc = x, nu = nu,
  range = range, sigma = sigma, m = 2,
  parameterization = "matern"
)

# Get the precision matrix:
prec_matrix <- precision(op_cov)

Description

Function to get the precision matrix of a spacetimeobj object

Usage

## S3 method for class 'spacetimeobj'
precision(
  object,
  user_kappa = NULL,
  user_sigma = NULL,
  user_gamma = NULL,
  user_rho = NULL,
  ...
)

Arguments

Value

The precision matrix.

See Also

simulate.spacetimeobj(), spacetime.operators()

Examples

s <- seq(from = 0, to = 20, length.out = 101)
t <- seq(from = 0, to = 20, length.out = 31)

op_cov <- spacetime.operators(space_loc = s, time_loc = t,
                             kappa = 5, sigma = 10, alpha = 1,
                             beta = 2, rho = 1, gamma = 0.05)
prec_matrix <- precision(op_cov)

Description

The function is used for computing kriging predictions based on data $Y_i = u(s_i) + \epsilon_i$, where $\epsilon$ is mean-zero Gaussian measurement noise and $u(s)$ is defined by a fractional SPDE $(\kappa^2 I - \Delta)^{\alpha/2} (\tau u(s)) = W$, where $W$ is Gaussian white noise and $\alpha = \nu + d/2$, where $d$ is the dimension of the domain.

Usage

## S3 method for class 'CBrSPDEobj'
predict(
  object,
  A,
  Aprd,
  Y,
  sigma.e,
  mu = 0,
  compute.variances = FALSE,
  posterior_samples = FALSE,
  n_samples = 100,
  only_latent = FALSE,
  ...
)

Arguments

Value

A list with elements

Examples

set.seed(123)
# Sample a Gaussian Matern process on R using a rational approximation
kappa <- 10
sigma <- 1
nu <- 0.8
sigma.e <- 0.3
range <- 0.2

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

tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
  (4 * pi)^(1 / 2) * gamma(nu + 1 / 2)))

# Compute the covariance-based rational approximation
op_cov <- matern.operators(
  loc_mesh = x, nu = nu,
  range = range, sigma = sigma, d = 1, m = 2,
  parameterization = "matern"
)

# Sample the model
u <- simulate(op_cov)

# Create some data
obs.loc <- runif(n = 10, min = 0, max = 1)
A <- rSPDE.A1d(x, obs.loc)
Y <- as.vector(A %*% u + sigma.e * rnorm(10))

# compute kriging predictions at the FEM grid
A.krig <- rSPDE.A1d(x, x)
u.krig <- predict(op_cov,
  A = A, Aprd = A.krig, Y = Y, sigma.e = sigma.e,
  compute.variances = TRUE
)

plot(obs.loc, Y,
  ylab = "u(x)", xlab = "x", main = "Data and prediction",
  ylim = c(
    min(u.krig$mean - 2 * sqrt(u.krig$variance)),
    max(u.krig$mean + 2 * sqrt(u.krig$variance))
  )
)
lines(x, u.krig$mean)
lines(x, u.krig$mean + 2 * sqrt(u.krig$variance), col = 2)
lines(x, u.krig$mean - 2 * sqrt(u.krig$variance), col = 2)

Description

The function is used for computing kriging predictions based on data $Y_i = u(s_i) + \epsilon_i$, where $\epsilon$ is mean-zero Gaussian measurement noise and $u(s)$ is defined by a SPDE as described in matern2d.operators().

Usage

## S3 method for class 'CBrSPDEobj2d'
predict(
  object,
  A,
  Aprd,
  Y,
  sigma.e,
  mu = 0,
  compute.variances = FALSE,
  posterior_samples = FALSE,
  n_samples = 100,
  only_latent = FALSE,
  ...
)

Arguments

Value

A list with elements

Examples

library(fmesher)
 n_loc <- 2000
 loc_2d_mesh <- matrix(runif(n_loc * 2), n_loc, 2)
 mesh_2d <- fm_mesh_2d(loc = loc_2d_mesh, cutoff = 0.01, max.edge = c(0.1, 0.5))
 op <- matern2d.operators(hx = 0.08, hy = 0.08, hxy = 0.5, nu = 0.5, 
 sigma = 1, mesh = mesh_2d)
 u <- simulate(op)
 n.obs <- 2000
 obs.loc <- cbind(runif(n.obs),runif(n.obs))
 A <- fm_basis(mesh_2d,obs.loc)
 sigma.e <- 0.1
 Y <- as.vector(A%*%u + sigma.e*rnorm(n.obs))
 A <- op$make_A(obs.loc)
 proj <- fm_evaluator(mesh_2d, dims = c(100, 100),
             xlim = c(0,1), ylim = c(0,1))
 Aprd <- op$make_A(proj$lattice$loc)
 u.krig <- predict(op, A = A, Aprd = Aprd, Y = Y, sigma.e = sigma.e)

Description

Prediction of a mixed effects regression model on a metric graph.

Usage

## S3 method for class 'rspde_lme'
predict(
  object,
  newdata = NULL,
  loc = NULL,
  time = NULL,
  mesh = FALSE,
  which_repl = NULL,
  compute_variances = FALSE,
  posterior_samples = FALSE,
  n_samples = 100,
  sample_latent = FALSE,
  return_as_list = FALSE,
  return_original_order = TRUE,
  ...,
  data = deprecated()
)

Arguments