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:

• David Bolin [email protected]

• 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 anisotropic inlabru mapper

Usage

bru_get_mapper.inla_rspde_anisotropic2d(model, ...)

Arguments

Description

rSPDE inlabru mapper

Usage

bru_get_mapper.inla_rspde_fintrinsic(model, ...)

Arguments

Description

rSPDE stationary inlabru mapper

Usage

bru_get_mapper.inla_rspde_matern1d(model, ...)

ibm_n.bru_mapper_inla_rspde_matern1d(mapper, ...)

ibm_values.bru_mapper_inla_rspde_matern1d(mapper, ...)

ibm_jacobian.bru_mapper_inla_rspde_matern1d(mapper, input, ...)

Arguments

Description

rSPDE space time inlabru mapper

Usage

bru_get_mapper.inla_rspde_spacetime(model, ...)

Arguments

Description

rSPDE inlabru mapper

Usage

bru_get_mapper.intrinsic_matern(model, ...)

ibm_n.bru_mapper_intrinsic_matern(mapper, ...)

ibm_values.bru_mapper_intrinsic_matern(mapper, ...)

ibm_jacobian.bru_mapper_intrinsic_matern(mapper, input, ...)

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,
  nu = NULL,
  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

Generic for computing the covariance between mesh nodes and locations for model objects.

Usage

cov_function_mesh(object, ...)

Arguments

Description

Generic for computing the covariance between mesh nodes for model objects.

Usage

covariance_mesh(object, ...)

Arguments

Description

Creates train and test splits for cross-validation by handling multiple data types and supporting k-fold, leave-one-out (LOO), and leave-percentage-out (LPO) methods. Handles missing values and maintains data structure across multiple datasets.

Usage

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

Arguments

Details

The function handles NA values by removing rows with any missing values before creating splits. For multiple datasets, indices are mapped back to their original positions in each dataset.

Value

A list where each element contains:

Description

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

Usage

cross_validation(
  models,
  model_names = NULL,
  scores = c("mae", "mse", "crps", "scrps", "dss"),
  cv_type = c("k-fold", "loo", "lpo"),
  weight_thr = NULL,
  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,
  model_options_bru = list(),
  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,
  type_rational_approximation = "chebfunLB",
  d = NULL,
  x_min = NULL,
  wl2_table = NULL
)

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,
  time = NULL,
  bru = FALSE,
  tibble = FALSE,
  drop_na = FALSE,
  drop_all_na = TRUE
)

Arguments

Value

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

Description

hybrid.spde constructs a finite element/rational SPDE approximation of the model

$Y(s) = \beta L^{-\alpha/2} X(s) + \tau^{-1} L^{-\alpha/2} W(s),$

where $L = \kappa(s)^2 - \Delta$ is the (possibly non-stationary) Whittle-Matern differential operator, $\tau$ is a positive scale that affects only the random component, $X$ is a (matrix of) covariate field(s) supplied at the mesh nodes, $\beta$ is a vector of regression coefficients and $W$ is Gaussian white noise.

Usage

hybrid.spde(
  X = NULL,
  beta_X = NULL,
  kappa = NULL,
  tau = NULL,
  kappa_mu = NULL,
  theta = NULL,
  B.tau = NULL,
  B.kappa = NULL,
  B.sigma = NULL,
  B.range = NULL,
  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("brasil", "chebfun", "chebfunLB"),
  check_stationarity = TRUE
)

Arguments

Value

hybrid.spde returns an object of class "hybrid_spde" which also inherits from the class of the object returned by spde.matern.operators() (spde_matern_operator plus either rSPDEobj or CBrSPDEobj depending on type). The object stores the covariate matrix X, the regression coefficients beta_X, the cached operator object op_mean used to evaluate $L^{-\alpha/2}$, and the evaluated mean function mu at the mesh nodes.

See Also

spde.matern.operators(), simulate.hybrid_spde(), predict.hybrid_spde()

Examples

set.seed(1)
x <- seq(from = 0, to = 1, length.out = 51)
kappa <- 10
tau <- 0.5
alpha <- 1.3
X <- cbind(sin(2 * pi * x), cos(2 * pi * x))
beta_X <- c(1.5, -0.5)
op <- hybrid.spde(
  X = X, beta_X = beta_X,
  kappa = kappa, tau = tau, alpha = alpha,
  loc_mesh = x, d = 1, parameterization = "spde"
)
Y <- simulate(op, nsim = 1)

Description

rSPDE inlabru mapper

Usage

ibm_n.bru_mapper_inla_rspde_fintrinsic(mapper, ...)

ibm_values.bru_mapper_inla_rspde_fintrinsic(mapper, ...)

ibm_jacobian.bru_mapper_inla_rspde_fintrinsic(mapper, input, ...)

bru_get_mapper.inla_rspde(model, ...)

ibm_n.bru_mapper_inla_rspde(mapper, ...)

ibm_values.bru_mapper_inla_rspde(mapper, ...)

ibm_jacobian.bru_mapper_inla_rspde(mapper, input, ...)

Arguments

Examples

if (rspde_safe_inla() &&
  requireNamespace("inlabru", quietly = TRUE)) {
  library(INLA)
  library(inlabru)

  set.seed(123)
  m <- 50
  loc_2d_mesh <- matrix(runif(m * 2), m, 2)
  mesh_2d <- inla.mesh.2d(
    loc = loc_2d_mesh,
    cutoff = 0.1,
    max.edge = c(0.5, 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 = 1,
    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 = 1
  )

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


  rspde_fit <- bru(cmp, data = data_df)
  summary(rspde_fit)
}

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("brasil", "chebfun", "chebfunLB"),
  fem_mesh_matrices = NULL,
  scaling = NULL,
  opts = 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. Unless supplied, the scaling is computed as the lowest positive eigenvalue of sqrt(solve(c0))%*%g1%*%sqrt(solve(c0)). opts provides a list of options for the numerical calculation of the scaling factor: tol, the relative tolerance on the eigenvalue, and maxitr, the maximum number of Lanczos steps.

Value

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

Examples

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

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} (\tau u) = \mathcal{W}$

Usage

intrinsic.operators(
  tau = NULL,
  beta = NULL,
  G = NULL,
  C = NULL,
  d = NULL,
  mesh = NULL,
  graph = NULL,
  loc_mesh = NULL,
  m = 1,
  compute_higher_order = FALSE,
  return_block_list = FALSE,
  type_rational_approximation = c("brasil", "chebfun", "chebfunLB"),
  fem_mesh_matrices = NULL,
  scaling = NULL,
  opts = NULL
)

Arguments

Details

The covariance operator

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

is approximated based on a rational approximation. The Laplacian is equipped with homogeneous Neumann boundary conditions and a zero-mean constraint is additionally imposed to obtained a non-intrinsic model. The scaling is computed as the lowest positive eigenvalue of sqrt(solve(c0))%*%g1sqrt(solve(c0)). opts provides a list of options for the numerical calculation of the scaling factor: tol, the relative tolerance on the eigenvalue, and maxitr, the maximum number of Lanczos steps.

Value

intrinsic.operators returns an object of class "intrinsicCBrSPDEobj".

Examples

x <- seq(from = 0, to = 10, length.out = 201)
beta <- 1
alpha <- 1
op <- intrinsic.operators(tau = 1, 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 = 0, alpha = 0, 
  beta = beta, L = 10, d = 1),
  col = 2, lty = 2
)

Description

Generic for computing the projection matrix that links observation locations to model mesh nodes.

Usage

make_A(object, ...)

Arguments

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("brasil", "chebfun", "chebfunLB", "wl2"),
  x_min = NULL,
  kappa_ref = NULL,
  wl2_table = NULL,
  variance_correction = c("none", "nodal"),
  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",
  wl2_table = NULL
)

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("brasil", "chebfun", "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

Generic function for posterior cross-validation.

Usage

posterior_crossvalidation(object, ...)

## S3 method for class 'list'
posterior_crossvalidation(object, ..., tibble = TRUE, return_indices = FALSE)

## Default S3 method:
posterior_crossvalidation(object, ...)

Arguments

Value

A list with elements mu, var and scores, and indices if requested. For a list of models, mu and var are lists with one element per model, and scores has one row per model.

Description

Performs cross-validation on objects fitted with rspde_lme and computes common predictive scores. For true_CV = FALSE (default, "pseudo" cross-validation) the model parameters are kept fixed at the estimates from the full fit and only the held-out points are masked at prediction time. When use_precomputed = TRUE the parameter-dependent structures (updated rSPDE operator and the precision matrix Q) are computed once and reused across folds, mirroring the advanced_options$precompute_data / na_test_idx path of predict.rspde_lme. For true_CV = TRUE the model is refit on each training fold (with the held-out response values set to NA). The previous fit is forwarded via previous_fit so that the optimisation starts from the full-data estimates.

Usage

## S3 method for class 'rspde_lme'
posterior_crossvalidation(
  object,
  scores = c("logscore", "crps", "scrps", "mae", "rmse"),
  mode = "k-fold",
  k = 10,
  percentage = 20,
  number_folds = 10,
  train_test_indices = NULL,
  true_CV = FALSE,
  factor = 1,
  tibble = TRUE,
  parallel_folds = FALSE,
  parallel_fitting = FALSE,
  n_cores = parallel::detectCores() - 1,
  print = FALSE,
  seed = NULL,
  return_indices = FALSE,
  use_precomputed = TRUE,
  data = NULL,
  nelder_mead_init = FALSE,
  ...
)