Skip to contents

This function returns a log-likelihood function for a Gaussian process with a Matern covariance function, that is observed under Gaussian measurement noise: \(Y_i = u(s_i) + \epsilon_i\), where \(\epsilon_i\) are iid mean-zero Gaussian variables. The latent model is approximated using the a rational approximation of the fractional SPDE model corresponding to the Gaussian process.


  sigma.e = NULL,
  mu = 0,
  user_nu = NULL,
  user_tau = NULL,
  user_kappa = NULL,
  user_sigma = NULL,
  user_range = NULL,
  parameterization = c("spde", "matern"),
  user_m = NULL,
  log_scale = TRUE,
  return_negative_likelihood = TRUE



The rational SPDE approximation, computed using matern.operators()


The observations, either a vector or a matrix where the columns correspond to independent replicates of observations.


An observation matrix that links the measurement location to the finite element basis.


IF non-null, the standard deviation of the measurement noise will be kept fixed in the returned likelihood.


Expectation vector of the latent field (default = 0).


If non-null, the shape parameter will be kept fixed in the returned likelihood.


If non-null, the tau parameter will be kept fixed in the returned likelihood. (Replaces sigma)


If non-null, the range parameter will be kept fixed in the returned likelihood.


If non-null, the standard deviation will be kept fixed in the returned likelihood.


If non-null, the range parameter will be kept fixed in the returned likelihood. (Replaces kappa)


If spde, then one will use the parameters tau and kappa. If matern, then one will use the parameters sigma and range.


If non-null, update the order of the rational approximation, which needs to be a positive integer.


Should the parameters be evaluated in log-scale?


Return minus the likelihood to turn the maximization into a minimization?


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


# this example illustrates how the function can be used for maximum
# likelihood estimation

# Sample a Gaussian Matern process on R using a rational approximation
nu <- 0.8
sigma <- 1
sigma.e <- 0.1
n.rep <- 10
n.obs <- 200
n.x <- 51
range <- 0.2
# create mass and stiffness matrices for a FEM discretization
x <- seq(from = 0, to = 1, length.out = n.x)
# 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, 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)
# \donttest{
# Define the negative likelihood function for optimization
# using CBrSPDE.matern.loglike
# Matern parameterization
loglike <- rSPDE.construct.matern.loglike(op_cov, Y, A, parameterization = "matern")

# The parameters can now be estimated by minimizing mlik with optim

# Choose some reasonable starting values depending on the size of the domain
theta0 <- c(
  get.initial.values.rSPDE(mesh.range = 1, dim = 1),
  log(0.1 * sd(as.vector(Y)))
# run estimation and display the results
theta <- optim(theta0, loglike,
  method = "L-BFGS-B"
  sigma = c(sigma, exp(theta$par[1])), range = c(range, 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")
#>               sigma    range       nu  sigma.e
#> Truth     1.0000000 0.200000 0.800000 0.100000
#> Estimates 0.9973174 0.180028 0.978009 0.100262
# }