rSPDE.construct.matern.loglike: Constructor of Matern loglikelihood functions.

In rSPDE: Rational Approximations of Fractional Stochastic Partial Differential Equations

Constructor of Matern loglikelihood functions.

Description

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.

Usage

rSPDE.construct.matern.loglike(
  object,
  Y,
  A,
  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
)

Arguments

object	The rational SPDE approximation, computed using matern.operators()
Y	The observations, either a vector or a matrix where the columns correspond to independent replicates of observations.
A	An observation matrix that links the measurement location to the finite element basis.
sigma.e	IF non-null, the standard deviation of the measurement noise will be kept fixed in the returned likelihood.
mu	Expectation vector of the latent field (default = 0).
user_nu	If non-null, the shape parameter will be kept fixed in the returned likelihood.
user_tau	If non-null, the tau parameter will be kept fixed in the returned likelihood. (Replaces sigma)
user_kappa	If non-null, the range parameter will be kept fixed in the returned likelihood.
user_sigma	If non-null, the standard deviation will be kept fixed in the returned likelihood.
user_range	If non-null, the range parameter will be kept fixed in the returned likelihood. (Replaces kappa)
parameterization	If spde, then one will use the parameters tau and kappa. If matern, then one will use the parameters sigma and range.
user_m	If non-null, update the order of the rational approximation, which needs to be a positive integer.
log_scale	Should the parameters be evaluated in log-scale?
return_negative_likelihood	Return minus the likelihood to turn the maximization into a minimization?

Value

The log-likelihood function. The parameters of the returned function are given in the order sigma, kappa, 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

set.seed(123)
# 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)

# 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"
)
print(data.frame(
  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")
))

rSPDE documentation built on Nov. 6, 2023, 1:06 a.m.