Skip to content

Commit

Permalink
adding fractional intrinsic inla model
Browse files Browse the repository at this point in the history
  • Loading branch information
David Bolin authored and David Bolin committed Dec 18, 2024
1 parent 9f3dd0b commit 74e2e05
Show file tree
Hide file tree
Showing 21 changed files with 1,604 additions and 596 deletions.
5 changes: 5 additions & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,7 @@ export(gg_df)
export(glance)
export(graph_data_rspde)
export(intrinsic.matern.operators)
export(intrinsic.operators)
export(matern.covariance)
export(matern.operators)
export(matern.rational)
Expand Down Expand Up @@ -137,6 +138,10 @@ if (getRversion() >= "3.6.0") {
S3method(inlabru::ibm_n, bru_mapper_metric_graph)
S3method(inlabru::ibm_values, bru_mapper_metric_graph)
S3method(inlabru::ibm_jacobian, bru_mapper_metric_graph)
S3method(inlabru::bru_mapper, metric_graph_dirichlet)
S3method(inlabru::ibm_n, bru_mapper_metric_graph_dirichlet)
S3method(inlabru::ibm_values, bru_mapper_metric_graph_dirichlet)
S3method(inlabru::ibm_jacobian, bru_mapper_metric_graph_dirichlet)
}

import(Matrix)
Expand Down
887 changes: 447 additions & 440 deletions R/inla_rspde.R
View file

Large diffs are not rendered by default.

150 changes: 150 additions & 0 deletions R/inla_rspde_intrinsic.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#' The current options are "beta" or "lognormal". The default is "lognormal".
#' @param nu.prec.inc Amount to increase the precision in the beta prior
#' distribution. Check details below.
#' @param diagonal Number added to diagonal of Q for increased stability.
#' @param type.rational.approx Which type of rational approximation
#' should be used? The current types are "chebfun", "brasil" or "chebfunLB".
#' @param shared_lib String specifying which shared library to use for the Cgeneric
Expand Down Expand Up @@ -54,6 +55,7 @@ rspde.intrinsic <- function(mesh,
prior.nu = NULL,
prior.nu.dist = "lognormal",
nu.prec.inc = 0.01,
diagonal = 1e-5,
type.rational.approx = "chebfun",
shared_lib = "detect",
debug = FALSE,
Expand Down Expand Up @@ -119,6 +121,8 @@ rspde.intrinsic <- function(mesh,
return(out)
}
Cmatrix <- to.inla.matrix(op$Q)
D <- Diagonal(dim(op$Q)[1], diagonal)
scaling <- RSpectra::eigs(as(G, "CsparseMatrix"), 2, which = "SM")$values[1]

list_args <-
list(
Expand All @@ -138,11 +142,13 @@ rspde.intrinsic <- function(mesh,
prior.nu.prec = prior.nu$prec,
prior.nu.logscale = prior.nu$logscale,
nu_upper_bound = nu.upper.bound,
scaling = scaling,
prior.nu.dist = prior.nu.dist,
Q = Cmatrix,
C = C,
Ci = Ci,
G = G,
D = D,
rational_table = rational_table
)

Expand Down Expand Up @@ -444,3 +450,147 @@ rspde.intrinsic.matern <- function(mesh,
model$rspde_version <- as.character(packageVersion("rSPDE"))
return(model)
}


#' result summary for intrinsic models
#' @noRd
rspde.intrinsic.result <- function(inla, name, rspde,
compute.summary = TRUE,
n_samples = 5000,
n_density = 1024) {

nu.upper.bound <- rspde$nu_upper_bound
result <- list()

if (!rspde$est_nu) {
row_names <- c("tau")
} else {
row_names <- c("tau", "nu")
}

result$summary.values <- inla$summary.random[[name]]

if (!is.null(inla$marginals.random[[name]])) {
result$marginals.values <- inla$marginals.random[[name]]
}

name_theta1 <- "tau"
name_theta2 <- "nu"

name_theta1_model <- "tau"
name_theta2_model <- "nu"

result[[paste0("summary.log.", name_theta1_model)]] <- INLA::inla.extract.el(
inla$summary.hyperpar,
paste("Theta1 for ", name, "$", sep = "")
)
rownames(result[[paste0("summary.log.", name_theta1_model)]]) <- paste0("log(", name_theta1_model, ")")

if (rspde$est_nu) {
result$summary.logit.nu <- INLA::inla.extract.el(inla$summary.hyperpar,
paste("Theta2 for ", name, "$", sep = ""))
rownames(result$summary.logit.nu) <- "logit(nu)"
}

if (!is.null(inla$marginals.hyperpar[[paste0("Theta1 for ", name)]])) {
result[[paste0("marginals.log.", name_theta1_model)]] <- INLA::inla.extract.el(
inla$marginals.hyperpar,
paste("Theta1 for ", name, "$", sep = "")
)
names(result[[paste0("marginals.log.", name_theta1_model)]]) <- name_theta1_model

if (rspde$est_nu) {
result$marginals.logit.nu <- INLA::inla.extract.el(
inla$marginals.hyperpar,
paste("Theta2 for ", name, "$", sep = "")
)
names(result$marginals.logit.nu) <- "nu"
}


result[[paste0("marginals.", name_theta1)]] <- lapply(
result[[paste0("marginals.log.", name_theta1)]],
function(x) {
INLA::inla.tmarginal(function(y) exp(y), x)
}
)

if (rspde$est_nu) {
result$marginals.nu <- lapply(result$marginals.logit.nu,
function(x) {
INLA::inla.tmarginal(
function(y) {
nu.upper.bound * exp(y) / (1 + exp(y))
}, x)
})
}

if (compute.summary) {
norm_const <- function(density_df) {
min_x <- min(density_df[, "x"])
max_x <- max(density_df[, "x"])
denstemp <- function(x) {
dens <- sapply(x, function(z) {
if (z < min_x) {
return(0)
} else if (z > max_x) {
return(0)
} else {
return(approx(
x = density_df[, "x"],
y = density_df[, "y"], xout = z
)$y)
}
})
return(dens)
}
norm_const <- stats::integrate(
f = function(z) {
denstemp(z)
}, lower = min_x, upper = max_x,
subdivisions = nrow(density_df),
stop.on.error = FALSE
)$value
return(norm_const)
}

norm_const_theta1 <- norm_const(result[[paste0("marginals.", name_theta1)]][[name_theta1]])
result[[paste0("marginals.", name_theta1)]][[name_theta1]][, "y"] <-
result[[paste0("marginals.", name_theta1)]][[name_theta1]][, "y"] / norm_const_theta1

result[[paste0("summary.", name_theta1)]] <- create_summary_from_density(result[[paste0("marginals.", name_theta1)]][[name_theta1]],
name = name_theta1
)


if (rspde$est_nu) {
norm_const_nu <- norm_const(result$marginals.nu$nu)
result$marginals.nu$nu[, "y"] <-
result$marginals.nu$nu[, "y"] / norm_const_nu

result$summary.nu <- create_summary_from_density(result$marginals.nu$nu,
name = "nu"
)
}
}
}

result$n_par <- length(rspde$start.theta)

class(result) <- c("rspde_result", "rspde_intrinsic")
result$stationary <- TRUE
result$parameterization <- "spde"

result$params <- c(name_theta1)
if (rspde$est_nu) {
result$params <- c(result$params, "nu")
}

if (!is.null(result$summary.nu)) {
if(nu.upper.bound - result$summary.nu$mean < 0.1 || nu.upper.bound - result$summary.nu$mode < 0.1){
warning("the mean or mode of nu is very close to nu.upper.bound, please consider increasing nu.upper.bound, and refitting the model.")
}
}

return(result)
}
41 changes: 40 additions & 1 deletion R/inla_rspde_spacetime.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@
#' default values will be used. The `mean` value is also used as starting value for gamma.
#' @param prior.precision A precision matrix for \eqn{\log(\kappa), \log(\sigma), \log(\gamma), \rho}. This matrix replaces the precision
#' element from `prior.kappa`, `prior.sigma`, `prior.gamma`, and `prior.rho` respectively. For dimension 1 `prior.precision` must be a 4x4 matrix. For dimension 2, \eqn{\rho} is a vector of length 2, so in this case `prior.precision` must be a 5x5 matrix. If `NULL`, a diagonal precision matrix with default values will be used.
#' @param graph_dirichlet For models on metric graphs, use Dirichlet vertex conditions at vertices of degree 1?
#' @param bounded_rho Logical. Should `rho` be bounded to ensure the existence, uniqueness, and well-posedness of the solution? Defaults to `TRUE`.
#' Note that this bounding is not a strict condition; there may exist values of rho beyond the upper bound that still satisfy these properties.
#' When `bounded_rho = TRUE`, the `rspde_lme` models enforce bounded `rho` for consistency.
Expand Down Expand Up @@ -82,6 +83,7 @@ rspde.spacetime <- function(mesh_space = NULL,
prior.rho = NULL,
prior.gamma = NULL,
prior.precision = NULL,
graph_dirichlet = TRUE,
bounded_rho = TRUE,
shared_lib = "detect",
debug = FALSE,
Expand Down Expand Up @@ -119,7 +121,8 @@ rspde.spacetime <- function(mesh_space = NULL,
rho = prior.rho$mean,
alpha = alpha,
beta = beta,
bounded_rho = bounded_rho
bounded_rho = bounded_rho,
graph_dirichlet = graph_dirichlet
)

default_precision <- 0.1
Expand Down Expand Up @@ -234,6 +237,9 @@ rspde.spacetime <- function(mesh_space = NULL,
model$mesh <- op$mesh_space
} else{
model$mesh <- graph
if(graph_dirichlet) {
class(model$mesh) <- c("metric_graph_dirichlet", class(model$mesh))
}
}
model$time_mesh <- op$mesh_time
model$rspde_version <- as.character(packageVersion("rSPDE"))
Expand All @@ -256,6 +262,11 @@ rspde.spacetime <- function(mesh_space = NULL,
#' S3method(inlabru::ibm_n, bru_mapper_metric_graph)
#' S3method(inlabru::ibm_values, bru_mapper_metric_graph)
#' S3method(inlabru::ibm_jacobian, bru_mapper_metric_graph)
#' S3method(inlabru::bru_mapper, metric_graph_dirichlet)
#' S3method(inlabru::ibm_n, bru_mapper_metric_graph_dirichlet)
#' S3method(inlabru::ibm_values, bru_mapper_metric_graph_dirichlet)
#' S3method(inlabru::ibm_jacobian, bru_mapper_metric_graph_dirichlet)

#' }
#'

Expand Down Expand Up @@ -294,3 +305,31 @@ ibm_jacobian.bru_mapper_metric_graph <- function(mapper, input, ...) {
}
mapper[["mesh"]][["fem_basis"]](input)
}

## Dirichlet mappers
#' @noRd
bru_mapper.metric_graph_dirichlet <- function(mesh, ...) {
mapper <- list(mesh = mesh)
inlabru::bru_mapper_define(mapper, new_class = "bru_mapper_metric_graph_dirichlet")
}

#' @noRd
ibm_n.bru_mapper_metric_graph_dirichlet <- function(mapper, ...) {
nrow(mapper[["mesh"]][["mesh"]][["VtE"]]) - sum(mapper[["mesh"]]$get_degrees()==1)
}

#' @noRd
ibm_values.bru_mapper_metric_graph_dirichlet <- function(mapper, ...) {
seq_len(nrow(mapper[["mesh"]][["mesh"]][["VtE"]]) - sum(mapper[["mesh"]]$get_degrees()==1))
}

#' @noRd
ibm_jacobian.bru_mapper_metric_graph_dirichlet <- function(mapper, input, ...) {
if (is.null(input)) {
return(Matrix::Matrix(0, 0, inlabru::ibm_n(mapper)))
}

ind.field <- setdiff(1:nrow(mapper[["mesh"]][["mesh"]][["VtE"]]),
which(mapper[["mesh"]]$get_degrees()==1))
mapper[["mesh"]][["fem_basis"]](input)[,ind.field]
}
Loading

0 comments on commit 74e2e05

Please sign in to comment.