if (FALSE) {
if (require.nowarnings("INLA")) {
# Generate mesh and SPDE model
n.lattice <- 10 # increase for more interesting, but slower, examples
x <- seq(from = 0, to = 10, length.out = n.lattice)
lattice <- fmesher::fm_lattice_2d(x = x, y = x)
mesh <- fmesher::fm_rcdt_2d_inla(lattice = lattice, extend = FALSE, refine = FALSE)
spde <- inla.spde2.matern(mesh, alpha = 2)
# Generate an artificial sample
sigma2.e <- 0.01
n.obs <- 100
obs.loc <- cbind(
runif(n.obs) * diff(range(x)) + min(x),
runif(n.obs) * diff(range(x)) + min(x)
)
Q <- inla.spde2.precision(spde, theta = c(log(sqrt(0.5)), log(sqrt(1))))
x <- inla.qsample(Q = Q)
A <- fmesher::fm_basis(mesh = mesh, loc = obs.loc)
Y <- as.vector(A %*% x + rnorm(n.obs) * sqrt(sigma2.e))
## Estimate the parameters using INLA
mesh.index <- inla.spde.make.index(name = "field", n.spde = spde$n.spde)
ef <- list(c(mesh.index, list(Intercept = 1)))
s.obs <- inla.stack(data = list(y = Y), A = list(A), effects = ef, tag = "obs")
s.pre <- inla.stack(data = list(y = NA), A = list(1), effects = ef, tag = "pred")
stack <- inla.stack(s.obs, s.pre)
formula <- y ~ -1 + Intercept + f(field, model = spde)
result <- inla(
formula = formula, family = "normal", data = inla.stack.data(stack),
control.predictor = list(
A = inla.stack.A(stack),
compute = TRUE
),
control.compute = list(
config = TRUE,
return.marginals.predictor = TRUE
),
num.threads = 1
)
## Calculate contour map with two levels
map <- contourmap.inla(result,
stack = stack, tag = "pred",
n.levels = 2, alpha = 0.1, F.limit = 0.1,
max.threads = 1
)
## Plot the results
cols <- contourmap.colors(map,
col = heat.colors(100, 1),
credible.col = grey(0.5, 1)
)
image(matrix(map$M[mesh$idx$lattice], n.lattice, n.lattice), col = cols)
}
}
Run the code above in your browser using DataLab