## In this example, we calculate the excursion function
## for a partially observed AR process.
if (FALSE) {
if (require.nowarnings("INLA")) {
## Sample the process:
rho <- 0.9
tau <- 15
tau.e <- 1
n <- 100
x <- 1:n
mu <- 10 * ((x < n / 2) * (x - n / 2) + (x >= n / 2) * (n / 2 - x) + n / 4) / n
Q <- tau * sparseMatrix(
i = c(1:n, 2:n), j = c(1:n, 1:(n - 1)),
x = c(1, rep(1 + rho^2, n - 2), 1, rep(-rho, n - 1)),
dims = c(n, n), symmetric = TRUE
)
X <- mu + solve(chol(Q), rnorm(n))
## measure the sampled process at n.obs random locations
## under Gaussian measurement noise.
n.obs <- 50
obs.loc <- sample(1:n, n.obs)
A <- sparseMatrix(
i = 1:n.obs, j = obs.loc, x = rep(1, n.obs),
dims = c(n.obs, n)
)
Y <- as.vector(A %*% X + rnorm(n.obs) / sqrt(tau.e))
## Estimate the parameters using INLA
ef <- list(c(list(ar = x), list(cov = mu)))
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 + cov + f(ar, model = "ar1")
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
)
)
## calculate the level 0 positive excursion function
res.qc <- excursions.inla(result,
stack = stack, tag = "pred", alpha = 0.99, u = 0,
method = "QC", type = ">", max.threads = 2
)
## plot the excursion function and marginal probabilities
plot(res.qc$rho,
type = "l",
main = "marginal probabilities (black) and excursion function (red)"
)
lines(res.qc$F, col = 2)
}
}
Run the code above in your browser using DataLab