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)
Run the code above in your browser using DataLab