Latent Models

Latent models provide a structured way to define commonly used Gaussian Markov Random Fields (GMRFs) by specifying their key components: hyperparameters, precision matrices, mean vectors, and constraints. This abstraction enables easy construction and combination of standard models used in spatial statistics, time series analysis, and Bayesian modeling.

Overview

The LatentModel interface standardizes GMRF construction through a simple pattern: 2. Define the model structure (e.g., temporal, spatial, independence)

3. Specify hyperparameters and their validation

4. Construct GMRFs with automatic constraint handling

using GaussianMarkovRandomFields

# Define a temporal AR1 model
ar1 = AR1Model(100)
hyperparameters(ar1)  # (τ = Real, ρ = Real)

# Construct GMRF with specific parameters
gmrf = ar1(τ=2.0, ρ=0.8)  # Returns GMRF or ConstrainedGMRF automatically

Latent prior hierarchy

AbstractLatentPrior
   ├── LatentModel              ← Gaussian: precision_matrix, mean, callable
   └── NonGaussianLatentPrior   ← non-Gaussian: requires local_quadratic

LatentModel is the abstract type for Gaussian latent priors — every model that ships with this package (AR, RW, IID, Matérn, Besag, BYM2, …) is a LatentModel. NonGaussianLatentPrior is a separate sibling for priors whose log-density is not quadratic in x; instead of a materialised GMRF, they expose a local_quadratic method that re-linearises the prior around any reference point. The Gaussian Approximation reference page covers the iterated-linearisation Newton driven by these priors.

GaussianMarkovRandomFields.AbstractLatentPrior Type

AbstractLatentPrior

Top of the latent-prior hierarchy. Anything that can be the prior side of gaussian_approximation is an AbstractLatentPrior. Two children:

• LatentModel: a Gaussian latent prior. (Q, μ) are determined by the hyperparameters alone, so the prior materialises as a GMRF.

• NonGaussianLatentPrior: a latent prior whose log-density is not quadratic in x. There is no canonical materialised GMRF; the prior is only meaningful at a chosen reference point through local_quadratic.

Every concrete subtype must implement:

• length(prior) — number of latent variables.

• hyperparameters(prior) — names/types of hyperparameters.

• model_name(prior) — Symbol used for CombinedModel parameter prefixing.

• constraints(prior; θ...) — nothing or (A, e).

source

GaussianMarkovRandomFields.LatentModel Type

LatentModel <: AbstractLatentPrior

Abstract type for Gaussian latent priors. (Q, μ) depend only on the hyperparameters θ, so the prior materialises as a GMRF / ConstrainedGMRF. In addition to the AbstractLatentPrior interface, every concrete subtype must implement:

• precision_matrix(model; θ...)

• mean(model; θ...)

• (model)(; θ...) -> AbstractGMRF

local_quadratic(::LatentModel, x_ref; θ...) has a default that uses these methods together with logpdf, so subtypes don't need to override it for the iterated path.

Usage

model = AR1Model(100)
gmrf = model(; τ = 2.0, ρ = 0.8)   # materialise a GMRF

source

GaussianMarkovRandomFields.NonGaussianLatentPrior Type

NonGaussianLatentPrior <: AbstractLatentPrior

Abstract type for latent priors whose log p(x | θ) is not quadratic in x. There is no canonical materialised GMRF; the prior is queried through local_quadratic at a reference point x_ref. In addition to the AbstractLatentPrior interface, every concrete subtype must implement:

• local_quadratic(prior, x_ref; θ...) -> LocalLatentQuadratic

Models of this type cannot be called as prior(; θ...) because materialisation is x-dependent. Use gaussian_approximation(prior, obs_lik; θ...) to drive the iterated-linearisation Newton; pass x0 explicitly if zeros aren't an appropriate starting point.

source

GaussianMarkovRandomFields.hyperparameters Method

hyperparameters(prior::AbstractLatentPrior) -> NamedTuple

NamedTuple describing the hyperparameter names and their expected types.

source

GaussianMarkovRandomFields.precision_matrix Method

precision_matrix(model::LatentModel; θ...) -> AbstractMatrix

Precision matrix of the Gaussian latent prior at hyperparameters θ.

source

Statistics.mean Method

mean(model::LatentModel; θ...) -> AbstractVector

Mean vector of the Gaussian latent prior at hyperparameters θ.

source

GaussianMarkovRandomFields.constraints Method

constraints(prior::AbstractLatentPrior; θ...) -> Union{Nothing, Tuple}

Linear-equality constraint information for the prior at hyperparameters θ. Either nothing (unconstrained) or a tuple (A, e) such that A x = e is enforced.

source

GaussianMarkovRandomFields.model_name Method

model_name(prior::AbstractLatentPrior) -> Symbol

Symbol used as a parameter-name suffix when this prior is composed with others in a CombinedModel (so e.g. τ from two priors becomes τ_ar1 and τ_besag).

source

Local quadratic interface

GaussianMarkovRandomFields.LocalLatentQuadratic Type

LocalLatentQuadratic{TQ, Th, T, Tx}

Local quadratic approximation of log p(x | θ) around a reference point x_ref, in natural form

log p(x | θ) ≈ logp_ref + (h - Q · x_ref)' (x - x_ref)
                        - ½ (x - x_ref)' Q (x - x_ref)

equivalently c + h' x - ½ x' Q x for an implicit c.

Fields

• Q::TQ: -∇²ₓ log p(x | θ) evaluated at x_ref.

• h::Th: ∇ₓ log p(x | θ)|_{x_ref} + Q · x_ref. The natural-form linear coefficient. Well-defined even for rank-deficient (intrinsic) Q where μ = Q⁻¹h is not.

• logp_ref::T: Exact log p(x_ref | θ) (not the quadratic's value at x_ref, which agrees only at x_ref itself). Carrying the exact value lets the Laplace marginal-likelihood formula use it at convergence without re-evaluating the prior.

• x_ref::Tx: The reference point.

source

GaussianMarkovRandomFields.local_quadratic Function

local_quadratic(m::AbstractLatentPrior, x_ref::AbstractVector; θ...) -> LocalLatentQuadratic

Local quadratic approximation of log p(x | θ) at x_ref for the latent prior m.

For LatentModel (Gaussian) the default method uses precision_matrix, mean, and the logpdf of the materialised GMRF — (Q, h) are independent of x_ref so the Newton iteration reduces to fixed-Q.

For NonGaussianLatentPrior subtypes there is no default; each concrete subtype must implement this method to specify the natural-form quadratic at x_ref.

source

local_quadratic(m::AutoDiffLatentPrior, x_ref; θ...) -> LocalLatentQuadratic

Natural-form quadratic of the AD-defined log-density at x_ref, via one gradient and one (sparse) Hessian evaluation: Q = -∇²f, h = ∇f + Q·x_ref, logp_ref = f(x_ref).

source

GaussianMarkovRandomFields.prior_logdensity Function

prior_logdensity(m::AbstractLatentPrior, x::AbstractVector; θ...) -> Real

Exact log p(x | θ) of the latent prior at x.

The iterated-linearisation line search evaluates the prior log-density at trial points where it needs only this scalar — not the full local quadratic. The default delegates to local_quadratic(m, x; θ...).logp_ref, so models that don't override it keep working unchanged. A NonGaussianLatentPrior whose local_quadratic is expensive (e.g. it assembles a sparse Hessian) can override prior_logdensity with a direct, allocation-light evaluation of log p(x | θ) to make the line search cheaper.

source

prior_logdensity(m::AutoDiffLatentPrior, x; θ...) -> Real

Direct primal evaluation f(x; θ...) — no AD, no Hessian. This is the cheap line-search hook: backtracking calls it per trial point instead of rebuilding the full local quadratic.

source

AD-defined (TMB-style) priors

AutoDiffLatentPrior defines a non-Gaussian latent prior implicitly through a scalar joint log-density f(x; θ) = log p(x | θ); the natural-form local quadratic is obtained by automatic differentiation (gradient + sparse Hessian), so you never hand-derive local_quadratic. It is the latent-side analogue of AutoDiffObservationModel and reuses the same DifferentiationInterface backends and sparse-Hessian machinery.

GaussianMarkovRandomFields.AutoDiffLatentPrior Type

AutoDiffLatentPrior{F, B, SB, HN, C} <: NonGaussianLatentPrior

A non-Gaussian latent prior defined implicitly by a scalar joint log-density f(x; θ...) = log p(x | θ). The natural-form local quadratic (local_quadratic) is obtained by automatic differentiation: the gradient ∇f and a (sparse) Hessian ∇²f give

Q = -∇²f(x_ref),   h = ∇f(x_ref) + Q · x_ref,   logp_ref = f(x_ref).

This is the latent-side analogue of AutoDiffObservationModel: it reuses the same DifferentiationInterface backends, the eltype-keyed preparation cache, and — with the SparseADLikelihoods extension loaded (SparseConnectivityTracer + SparseMatrixColorings) — sparse-Hessian detection and colouring. The Hessian sparsity is detected structurally (value-agnostic), so the pattern is stable across Newton iterates, which is what lets the GMRFWorkspace path reuse a single symbolic factorisation.

A monolithic TMB-style joint f(x, θ; data) is just an AutoDiffLatentPrior carrying the entire energy, combined with a ZeroLikelihood — marginal_loglikelihood then reduces to the Laplace marginal log ∫ exp f(x, θ) dx. The structured case (an AD-defined latent prior plus an exact closed-form likelihood) works equally well by composing this prior with any ObservationLikelihood.

Write logp_func in tracer-friendly form

With the default sparse Hessian backend, the sparsity pattern is detected by running logp_func on SparseConnectivityTracer tracers. Express the density as plain arithmetic on x; constructs that branch on a tracer value — e.g. Distributions.Normal(x[i], σ) (whose check_args tests σ > 0) — will error. Inline the closed form (-0.5*((y-x)/σ)^2 - log(σ) - 0.5log(2π)), or pass a dense / KnownHessianSparsityDetector hessian_backend.

Capturing data with the Enzyme default backend

grad_backend defaults to Enzyme when it's loaded. Enzyme can't prove a logp_func that captures a mutable array (e.g. observation data folded into a monolithic joint) is read-only, and raises EnzymeMutabilityException. For such a logp_func, pass grad_backend = AutoForwardDiff() (its hessian_backend follows automatically). Functions that capture only scalars are unaffected.

Hyperparameter gradients

With the ForwardDiff extension loaded, the Laplace marginal likelihood is differentiable in the hyperparameters: ForwardDiff.gradient(θ -> marginal_loglikelihood( prior, obs_lik, gaussian_approximation(prior, obs_lik; θ...); θ...), θ₀) is exact. When θ carries Dual partials, gaussian_approximation runs primal Newton and applies the Implicit Function Theorem (CHOLMOD can't factorize Dual matrices, and the sparsity tracer can't see Dual θ), so the returned posterior carries dx*/dθ and dQ/dθ.

Fields

• logp_func::F: log-density with signature (x; θ...) -> Real.

• n::Int: number of latent variables.

• grad_backend::B / hess_backend::SB: DI backends for ∇f and ∇²f.

• prep_cache::_ADPrepCache: eltype-keyed DI preparation cache.

• hyperparams::HN: tuple of hyperparameter names (Symbols).

• name::Symbol: model_name, used for CombinedModel parameter prefixing.

• constraints::C: nothing, or a fixed (A, e) linear-equality constraint A x = e.

source

Factor-graph (structured) priors

StructuredLatentPrior represents a non-Gaussian latent prior as a factor graph: a tuple of LatentFactorGroups, each a group of conditional factors that share one small per-factor log-density. The gradient and sparse Hessian are assembled by differentiating those small factor functions and scattering the contributions, so AD specialises only on the tiny factor functions rather than one opaque whole-model closure. It carries the same interface as AutoDiffLatentPrior (local_quadratic, prior_logdensity), so the rest of the pipeline treats them identically.

GaussianMarkovRandomFields.StructuredLatentPrior Type

StructuredLatentPrior{G, P, HN, C} <: NonGaussianLatentPrior

A non-Gaussian latent prior expressed as a factor graph: a tuple of LatentFactorGroups whose per-factor log-densities sum to log p(x | θ). The gradient and sparse Hessian are assembled by differentiating each group's small factor function and scattering the local contributions directly into the global structures — so the heavy machinery is generic, the only model-specific code AD specialises on is the handful of tiny factor functions, and the work per Newton iterate is O(nnz), not O(n²).

The factor-structured counterpart of AutoDiffLatentPrior (which differentiates one opaque whole-model closure). It carries the same NonGaussianLatentPrior interface (local_quadratic, prior_logdensity, the IFT hooks), so the GA / IFT / marginal-likelihood pipeline treats it identically.

Fields

• n: number of latent variables.

• groups: tuple of LatentFactorGroups.

• pattern: structural Hessian sparsity (typically the prior∪observation union), so the per-iterate Q aligns positionally with a reused GMRFWorkspace factorisation.

• posmaps: per-group, per-factor nzval positions of the factor's K×K block in pattern (precomputed once so assembly scatters into the sparse Q with no dense intermediate).

• hyperparams: tuple of hyperparameter names.

• name: model name.

• constraints: nothing, or a fixed (A, e) linear-equality constraint.

source

GaussianMarkovRandomFields.LatentFactorGroup Type

LatentFactorGroup{K, F}

A group of conditional-prior factors that share one log-density form. Each factor touches K latent variables; vars[i] lists the global latent indices of the i-th factor (by convention the factor's own variable first, then its parents). logp(vals, θ::NamedTuple) -> Real is the factor log-density evaluated at the local values vals = x[vars[i]] and hyperparameters θ.

Because every factor in a group calls the same logp, automatic differentiation specialises on one small (K-input) function per group — not on the whole-model log-density. This is what keeps the per-model compile cost small.

source

Available Models

Temporal Models

GaussianMarkovRandomFields.ARModel Type

ARModel{P}(n::Int; alg=<auto>, constraint=nothing)

A stationary autoregressive model of order P using the PACF parameterization.

The AR(P) model represents a temporal process: x[t] = φ₁x[t-1] + ... + φₚx[t-P] + ε[t]. Instead of parameterizing directly by the AR coefficients φ, the model uses partial autocorrelation (PACF) values θ₁, ..., θₚ which guarantee stationarity when |θₖ| < 1. The k-th PACF value θₖ is the correlation between x[t] and x[t-k] after removing the linear effect of the intervening lags x[t-1], ..., x[t-k+1]. The PACF values are converted to AR coefficients internally via the Durbin-Levinson recursion.

Type Aliases

• AR1Model = ARModel{1} — first-order autoregressive (backward compatible)

Hyperparameters

• P=1: τ (precision, τ > 0) and ρ (correlation, |ρ| < 1)

• P≥2: τ (precision) and pacf1, pacf2, ..., pacfP (PACF values, each in (-1, 1))

Precision Matrix

Constructed as Q = τ · L' · D · L where L is unit lower triangular with bandwidth P and D is diagonal, encoding the stationary AR(P) structure.

• P=1: Returns SymTridiagonal (efficient tridiagonal solver)

• P≥2: Returns SparseMatrixCSC (banded with bandwidth P)

Example

# AR(1) — backward compatible
model = AR1Model(100)
gmrf = model(τ=2.0, ρ=0.8)

# AR(2) via PACF
model2 = ARModel{2}(100)
gmrf2 = model2(τ=1.0, pacf1=0.7, pacf2=-0.3)

# AR(3)
model3 = ARModel{3}(100)
gmrf3 = model3(τ=1.0, pacf1=0.5, pacf2=-0.3, pacf3=0.2)

source

GaussianMarkovRandomFields.AR1Model Type

Backward-compatible alias for ARModel{1}.

source

GaussianMarkovRandomFields.RWModel Type

RWModel{Order}(n::Int; regularization=1e-5, alg=<auto>, additional_constraints=nothing, scale_model=false)

A random walk latent model of arbitrary order for constructing intrinsic GMRFs.

The RW model of order k represents a process where k-th order differences are i.i.d. Gaussian: Δᵏx[i] ~ N(0, τ⁻¹). This creates a singular precision matrix Q = τ * Dₖ'Dₖ with rank n-k, where Dₖ is the k-th order difference operator.

Variance scaling

With scale_model = true, the intrinsic precision is normalized following Sørbye & Rue (2014): the precision is multiplied by a constant so the geometric mean of the marginal variances (of the constrained model) is 1. This makes τ interpretable as a precision and comparable across n and across orders, so a prior on τ (e.g. PCPrior.Precision) transfers between models. Unscaled (scale_model = false, the default) reproduces the previous behavior, where τ = 1 implies an n- and order-dependent marginal variance. This is the same normalization BesagModel/BYM2Model apply.

Orders

• Order 1 (RW1): First differences x[i+1] - x[i] ~ N(0, τ⁻¹). Tridiagonal precision.

• Order 2 (RW2): Second differences x[i+2] - 2x[i+1] + x[i] ~ N(0, τ⁻¹). Pentadiagonal precision.

• Order k: k-th differences. Precision has bandwidth k.

Since the precision matrix is singular (rank n-k), the model is handled as an intrinsic GMRF with k constraints from the polynomial null space (degrees 0, 1, ..., k-1).

Type Aliases

• RW1Model = RWModel{1} — first-order random walk

• RW2Model = RWModel{2} — second-order random walk

Hyperparameters

• τ: Precision parameter (τ > 0)

Fields

• n::Int: Length of the process

• regularization::Float64: Small value added to diagonal after scaling (default 1e-5). Keep it small when scale_model = true: the variance normalization is computed at a negligible internal regularization, so a large regularization would shift the realized marginal variances away from the normalized scale.

• alg::Alg: LinearSolve algorithm (default: LDLtFactorization() for Order=1, CHOLMODFactorization() for Order≥2)

• additional_constraints::C: Optional additional constraints beyond the required null space constraints

• scale_factor::Float64: Sørbye–Rue variance-normalization factor multiplying the intrinsic precision (1.0 when scale_model = false)

Example

# First-order random walk
model = RW1Model(100)
gmrf = model(τ=1.0)

# Second-order random walk (smoother)
model2 = RW2Model(100)
gmrf2 = model2(τ=1.0)

# Third-order random walk
model3 = RWModel{3}(100)
gmrf3 = model3(τ=1.0)

source

GaussianMarkovRandomFields.RW1Model Type

Backward-compatible alias for RWModel{1}.

source

GaussianMarkovRandomFields.RW2Model Type

Convenience alias for RWModel{2}.

source

Spatial Models

GaussianMarkovRandomFields.MaternModel Type

MaternModel{F, S, Alg, C, P, M} <: LatentModel

A Matérn latent model for constructing spatial GMRFs from discretized Matérn SPDEs. Validated user-facing constructors live in the GaussianMarkovRandomFieldsFEM extension.

The κ-independent FEM matrices (the lumped mass matrix C and the stiffness matrix G) are assembled once at construction and cached in fem_matrices, so that repeated precision_matrix calls only redo the κ-dependent work.

source

GaussianMarkovRandomFields.BarrierModel Type

BarrierModel{F, Alg, C, P, M} <: LatentModel

A barrier Matérn latent model (Bakka et al., 2019): a non-stationary variant of the ν = 1 (α = 2) Matérn SPDE in which spatial correlation does not flow across designated barrier triangles of the mesh. Barrier triangles are given a small fixed range (a fraction of the normal range), so the field decorrelates sharply across them — useful for e.g. coastlines/islands blocking a marine field.

The precision keeps the same sparsity and cost as the stationary model and is assembled per region. With no barrier triangles it reduces exactly to MaternModel(disc; smoothness = 0) (the stationary ν = 1 Matérn). Its sparsity pattern is padded to a precomputed structural pattern, so it is identical for all hyperparameter values — as required by fixed-pattern workspaces.

Validated user-facing constructors live in the GaussianMarkovRandomFieldsFEM extension. The range-independent per-region FEM matrices are assembled once at construction and cached in fem_matrices.

source

GaussianMarkovRandomFields.barrier_triangles Function

barrier_triangles(disc::FEMDiscretization{2}, polygon) -> Vector{Int}

Return the ids of mesh triangles whose centroid lies inside polygon, the standard rule for tagging barrier triangles (Bakka et al., 2019). polygon is a closed simple polygon given either as an N×2 matrix of vertices or as a vector of points ((x, y) tuples or Vecs). The result can be passed straight to BarrierModel as barrier_cells.

source

GaussianMarkovRandomFields.BesagModel Type

BesagModel(adjacency::AbstractMatrix; regularization::Float64 = 1e-5, normalize_var::Bool = true, singleton_policy::Symbol = :gaussian, alg=CHOLMODFactorization())

A Besag model for spatial latent effects on graphs using intrinsic Conditional Autoregressive (CAR) structure.

The Besag model represents spatial dependence where each node's precision depends on its graph neighbors. This creates a graph Laplacian precision structure that's widely used for spatial smoothing on irregular lattices.

Mathematical Description

For a graph with adjacency matrix W, the precision matrix follows:

• Q[i,j] = -τ if nodes i and j are neighbors (W[i,j] = 1)

• Q[i,i] = τ * degree[i] where degree[i] = sum(W[i,:])

• All other entries are 0

Since this matrix is singular (rank n-1), we handle it as an intrinsic GMRF by: 2. Scaling by τ first, then adding small regularization (1e-5) to diagonal for numerical stability

3. Adding sum-to-zero constraint: sum(x) = 0

Hyperparameters

• τ: Precision parameter (τ > 0)

Fields

• adjacency::M: Adjacency matrix W (preserves input structure - sparse, SymTridiagonal, etc.)

• regularization::Float64: Small value added to diagonal after scaling (default 1e-5)

• alg::Alg: LinearSolve algorithm for solving linear systems

Example

# 4-node cycle graph - can use sparse, SymTridiagonal, or Matrix
W = sparse(Bool[0 1 0 1; 1 0 1 0; 0 1 0 1; 1 0 1 0])
model = BesagModel(W)
gmrf = model(τ=1.0)  # Returns ConstrainedGMRF with sum-to-zero constraint using CHOLMODFactorization

# Or specify custom algorithm
model = BesagModel(W, alg=LDLtFactorization())
gmrf = model(τ=1.0)

source

GaussianMarkovRandomFields.BYM2Model Type

BYM2Model(adjacency::AbstractMatrix; regularization::Float64 = 1e-5, normalize_var::Bool = true, singleton_policy::Symbol = :gaussian, alg=CHOLMODFactorization(), additional_constraints=nothing, iid_constraint=nothing)

Besag-York-Mollié model with improved BYM2 parameterization (Riebler et al. 2016).

The BYM2 model is a reparameterization of the classic BYM model for disease mapping that improves interpretability and facilitates prior specification. It combines spatial (ICAR) and unstructured (IID) random effects with intuitive mixing and precision parameters.

Mathematical Description

The BYM2 model is a 2n-dimensional latent field stacking: 2. A scaled spatial component: u* ~ ICAR with variance-normalized precision Q*

3. An unstructured component: v* ~ N(0, I)

The block-diagonal precision matrix is:

Q = [ τ/(1-φ) * Q*    0           ]
    [ 0               τ/φ * I     ]

In the linear predictor, effects are added: η[i] = ... + u_[i] + v_[i]

This parameterization ensures:

• Var(u*[i]) ≈ (1-φ)/τ (spatial variance)

• Var(v*[i]) = φ/τ (unstructured variance)

• Var(u_[i] + v_[i]) = 1/τ (total variance)

• φ controls the proportion of variance that is unstructured

• φ = 0: pure spatial model (like scaled Besag)

• φ = 1: pure unstructured model (like scaled IID)

Hyperparameters

• τ: Overall precision parameter (τ > 0)

• φ: Mixing parameter (0 < φ < 1), proportion of unstructured variance

Fields

• besag::BesagModel: The spatial component (variance-normalized)

• iid::IIDModel: The unstructured component

• n::Int: Number of spatial units

• alg::Alg: LinearSolve algorithm for solving linear systems

Identifiability and Constraints

When using BYM2 with a fixed intercept, the model can be unidentifiable because the unstructured component v* can absorb any constant shift. To handle this: 2. No intercept (recommended): Use y ~ 0 + bym2(region) in formulas

3. Constrain IID: Pass iid_constraint=:sumtozero to ensure identifiability

The spatial component u* always has a sum-to-zero constraint (one per connected component).

Example

# 4-node cycle graph
W = sparse(Bool[0 1 0 1; 1 0 1 0; 0 1 0 1; 1 0 1 0])

# Standard BYM2 (use with y ~ 0 + ...)
model = BYM2Model(W)
gmrf = model(τ=1.0, φ=0.5)  # Equal spatial and unstructured variance

# BYM2 with IID constraint (use with y ~ 1 + ... to include intercept)
model_constrained = BYM2Model(W; iid_constraint=:sumtozero)
gmrf = model_constrained(τ=1.0, φ=0.5)  # Both components sum to zero

# More spatial smoothing (90% spatial, 10% unstructured)
gmrf = model(τ=1.0, φ=0.1)

# More unstructured variation (10% spatial, 90% unstructured)
gmrf = model(τ=1.0, φ=0.9)

References

Riebler, A., Sørbye, S. H., Simpson, D., & Rue, H. (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research, 25(4), 1145-1165.

source

Barrier Matérn model

The barrier model ([5]) is a non-stationary variant of the Matérn SPDE in which spatial correlation does not flow across designated barrier regions of the domain — for example land barriers (coastlines, islands) that should block correlation in a marine field. The stationary Matérn correlates points purely by Euclidean distance, so it incorrectly couples points that are close in straight-line distance but physically separated by a barrier.

BarrierModel tags a subset of mesh triangles as barriers and gives them a small fixed range (range_fraction times the normal range), so the field decorrelates sharply across them. It keeps the same sparsity and computational cost as the stationary model, exposes the same precision_matrix(model; τ, range) interface, and reduces exactly to MaternModel(disc; smoothness = 0) when there are no barrier triangles.

using Ferrite, FerriteGmsh, Gmsh, LibGEOS   # activate the FEM extension

disc = FEMDiscretization(grid, ip, qr)

# Tag barrier triangles by their centroid lying inside a barrier polygon
# (e.g. land). Pass an N×2 matrix of vertices or a vector of points.
land  = [(-0.1, -1.0), (0.1, -1.0), (0.1, 1.0), (-0.1, 1.0)]
cells = barrier_triangles(disc, land)

model = BarrierModel(disc; barrier_cells = cells, range_fraction = 0.1)
x = model(τ = 1.0, range = 0.5)   # a GMRF whose correlation does not cross the barrier

You can also pass an explicit Vector{Int} of triangle ids as barrier_cells. The model is defined for the ν = 1 (α = 2) Matérn on 2D discretizations.

Independence Models

GaussianMarkovRandomFields.IIDModel Type

IIDModel(n::Int; alg=DiagonalFactorization(), constraint=nothing)

An independent and identically distributed (IID) latent model for constructing simple diagonal GMRFs.

The IID model represents independent Gaussian random variables with identical precision τ. This is the simplest possible latent model, equivalent to a scaled identity precision matrix.

Mathematical Description

Each element is independent: x[i] ~ N(0, τ⁻¹) for i = 1,...,n. The precision matrix is simply: Q = τ * I(n)

This model is useful for:

• Modeling independent effects or noise

• Baseline comparisons with structured models

• Teaching/demonstration purposes

Hyperparameters

• τ: Precision parameter (τ > 0)

Fields

• n::Int: Length of the IID process

• alg::Alg: LinearSolve algorithm for solving linear systems

• constraint::C: Optional constraint, either nothing or (A, e) tuple

Example

model = IIDModel(100)
gmrf = model(τ=2.0)  # Returns unconstrained GMRF

# With sum-to-zero constraint (common in INLA for separating global offset)
model = IIDModel(100, constraint=:sumtozero)
gmrf = model(τ=2.0)  # Returns ConstrainedGMRF with sum-to-zero constraint

source

GaussianMarkovRandomFields.FixedEffectsModel Type

FixedEffectsModel(n::Int; λ::Real = 1e-6, alg=DiagonalFactorization())

Weakly regularized fixed-effects latent component.

Encodes standard GLM-style fixed effects inside the latent vector with a small ridge precision λ * I(n).

• No hyperparameters (returns NamedTuple())

• No constraints

• Name: :fixed (used for parameter prefixing in CombinedModel)

Fields

• n::Int: Length of the fixed effects vector

• λ::Float64: Ridge regularization parameter

• alg::Alg: LinearSolve algorithm for solving linear systems

Example

model = FixedEffectsModel(10)
gmrf = model()  # Returns GMRF with precision λ * I(10) using DiagonalFactorization

# Or specify custom algorithm
model = FixedEffectsModel(10, alg=CHOLMODFactorization())
gmrf = model()

source

Model Composition

GaussianMarkovRandomFields.CombinedModel Type

CombinedModel(components::Vector{<:LatentModel}; alg=CHOLMODFactorization())

A combination of multiple LatentModel instances into a single block-structured GMRF.

This enables modeling with multiple latent components, such as the popular BYM model (Besag-York-Mollié) which combines spatial (Besag) and independent (IID) effects.

Mathematical Description

Given k component models with sizes n₁, n₂, ..., nₖ:

• Combined precision matrix: Q = blockdiag(Q₁, Q₂, ..., Qₖ)

• Combined mean vector: μ = vcat(μ₁, μ₂, ..., μₖ)

• Combined constraints: Block-diagonal constraint structure preserving individual constraints

Parameter Naming

To avoid conflicts when multiple models have the same hyperparameters (e.g., multiple τ), parameters are automatically prefixed with model names:

• Single occurrence: τ_besag, ρ_ar1

• Multiple occurrences: τ_besag, τ_besag_2, τ_besag_3

Fields

• components::Vector{<:LatentModel}: The individual latent models

• component_sizes::Vector{Int}: Cached sizes of each component

• total_size::Int: Total size of the combined model

• alg::Alg: LinearSolve algorithm for solving linear systems

Example - BYM Model

# BYM model: spatial Besag + independent IID effects
W = sparse(adjacency_matrix)  # Spatial adjacency
besag = BesagModel(W)         # Spatial component
iid = IIDModel(n)            # Independent component

# Vector constructor
bym = CombinedModel([besag, iid])
# Or variadic constructor (syntactic sugar)
bym = CombinedModel(besag, iid)

# With custom algorithm
bym = CombinedModel([besag, iid], alg=LDLtFactorization())

# Usage with automatically prefixed parameters
gmrf = bym(τ_besag=1.0, τ_iid=2.0)

source

GaussianMarkovRandomFields.component_model Function

component_model(model::CombinedModel, name::Symbol)

Retrieve a component latent model by name. Names follow the same convention as hyperparameter prefixes (e.g., :matern, :besag, :iid, :iid_2).

Also accessible via property syntax: model.matern, model.iid, etc.

Example

combined = CombinedModel(MaternModel(pts; smoothness=1), IIDModel(10))
matern = combined.matern          # or component_model(combined, :matern)
iid    = combined.iid             # or component_model(combined, :iid)

source

GaussianMarkovRandomFields.SeparableModel Type

SeparableModel{T<:Tuple{Vararg{LatentModel}}, Alg} <: LatentModel

Represents a separable (Kronecker product) latent model for multi-dimensional processes.

The precision matrix is the Kronecker product of the component precision matrices:

Q = Q_1 ⊗ Q_2 ⊗ ... ⊗ Q_N

where the rightmost component varies fastest in the vectorized representation.

For example, SeparableModel(temporal, spatial) produces Q_time ⊗ Q_space, yielding a block-tridiagonal structure when Q_time is tridiagonal (e.g., RW1).

Fields

• components::T: Tuple of LatentModel components

• alg: Linear solver algorithm

Example

# Space-time separable model (standalone)
temporal = RW1Model(n_time)
spatial = BesagModel(adjacency_matrix)
st_model = SeparableModel(temporal, spatial)  # Q = Q_time ⊗ Q_space

# Instantiate with hyperparameters
# Standalone: τ_rw1, τ_besag
# In CombinedModel: τ_rw1_separable, τ_besag_separable
gmrf = st_model(τ_rw1=1.0, τ_besag=2.0)

Notes

• Requires at least 2 components

• Component ordering: rightmost component varies fastest (e.g., space in time×space model)

• Follows R-INLA convention: Q = Q_group ⊗ Q_main

• Hyperparameters are suffixed with component model names (e.g., τ_rw1, τ_besag)

• When used in CombinedModel, an additional _separable suffix is added by the parent

• Constraints from each component are composed via Kronecker products with identity matrices

• Redundant constraints are automatically removed to ensure full row rank

source

Usage Examples

Basic Model Construction

# Temporal AR1 process
ar1 = AR1Model(100)
gmrf = ar1(τ=2.0, ρ=0.8)

# Spatial Matérn model from points (stores observation coordinates)
points = [0.0 0.0; 1.0 0.0; 0.5 1.0]  # N×2 matrix
matern = MaternModel(points; smoothness = 2)
gmrf = matern(τ=1.0, range=1.5)

# Convenience: evaluation matrix and observation model from stored points
A = evaluation_matrix(matern)
obs_model = PointEvaluationObsModel(matern, Normal)

# Spatial Besag model
W = sparse_adjacency_matrix
besag = BesagModel(W)
gmrf = besag(τ=1.0)  # Returns ConstrainedGMRF with sum-to-zero constraint

Model Composition

# Classic BYM model: spatial + independent effects
W = spatial_adjacency_matrix
bym = CombinedModel(BesagModel(W), IIDModel(n))

# Check combined hyperparameters
hyperparameters(bym)  # (τ_besag = Real, τ_iid = Real)

# Construct combined GMRF
gmrf = bym(τ_besag=1.0, τ_iid=2.0)

# BYM2 model: improved parameterization (recommended)
# Uses mixing parameter φ instead of separate precisions
bym2 = BYM2Model(W)
hyperparameters(bym2)  # (τ = Real, φ = Real)
gmrf = bym2(τ=1.0, φ=0.5)  # φ ∈ (0,1): proportion unstructured

# Continuous spatial field + independent effects
points = generate_observation_points()
spatial_matern = MaternModel(points; smoothness = 1)
independent_effects = IIDModel(length(points))
combined = CombinedModel(spatial_matern, independent_effects)
gmrf = combined(τ_matern=1.0, range_matern=2.0, τ_iid=0.1)

Smart Parameter Naming

CombinedModel automatically handles parameter naming conflicts:

# Single occurrence: no suffix
CombinedModel(AR1Model(10), BesagModel(W))  # τ_ar1, ρ_ar1, τ_besag

# Multiple occurrences: numbered suffixes  
CombinedModel(IIDModel(5), IIDModel(10), IIDModel(15))  # τ_iid, τ_iid_2, τ_iid_3

Advanced Compositions

# Spatiotemporal model with three components
W = spatial_adjacency(n_regions)
model = CombinedModel(
    BesagModel(W),           # Spatial structure
    RW1Model(n_time),        # Temporal trend
    IIDModel(n_regions * n_time)  # Independent effects
)

gmrf = model(τ_besag=1.0, τ_rw1=2.0, τ_iid=0.1)

Separable (Kronecker Product) Models

For multi-dimensional data with separable structure, use SeparableModel to compose components via Kronecker products. This is particularly efficient for space-time models:

# Space-time separable model: Q = Q_time ⊗ Q_space
# (space component varies fastest in vectorized form)
temporal = RW1Model(n_time)
spatial = BesagModel(W)
spacetime = SeparableModel(temporal, spatial)

# Hyperparameters use suffix naming
hyperparameters(spacetime)  # (τ_rw1 = Real, τ_besag = Real)

# Construct GMRF with block-tridiagonal structure
gmrf = spacetime(τ_rw1=1.0, τ_besag=2.0)

# 3-way separable model: Q = Q_time ⊗ Q_space ⊗ Q_age
three_way = SeparableModel(
    RW1Model(n_time),
    BesagModel(W),
    IIDModel(n_age)
)
gmrf = three_way(τ_rw1=1.0, τ_besag=2.0, τ_iid=0.5)

Key Features

Automatic Constraint Handling

Models automatically determine GMRF type based on constraints:

ar1 = AR1Model(10)     # No constraints → GMRF
rw1 = RW1Model(10)     # Sum-to-zero constraint → ConstrainedGMRF
mixed = CombinedModel(ar1, rw1)  # Inherits constraints → ConstrainedGMRF

Custom Constraints

Latent models support adding custom linear equality constraints Ax = e. This is particularly useful for INLA-style modeling where you want to learn a global offset separately from the latent field.

Models with Optional Constraints

For models without built-in constraints (AR1, IID, FixedEffects, Matern), use the constraint parameter:

# Sum-to-zero constraint (shorthand)
ar1 = AR1Model(100, constraint=:sumtozero)
gmrf = ar1(τ=1.0, ρ=0.8)  # Returns ConstrainedGMRF

# Custom constraint: fix first two elements to sum to 1
A = [1.0 1.0 0.0 0.0 0.0]  # 1×n constraint matrix
e = [1.0]                   # Constraint value
iid = IIDModel(5, constraint=(A, e))
gmrf = iid(τ=2.0)          # x[1] + x[2] = 1.0

Intrinsic Models with Additional Constraints

Intrinsic models (RW1, Besag) already have built-in constraints. Use additional_constraints to add extra constraints without removing the built-in ones:

# RW1 already has sum-to-zero; add constraint to fix first element
A_extra = [1.0 0.0 0.0 0.0 0.0]  # Fix x[1] = 0
e_extra = [0.0]
rw1 = RW1Model(5, additional_constraints=(A_extra, e_extra))
gmrf = rw1(τ=1.0)  # Both sum(x) = 0 AND x[1] = 0

# Besag with additional constraint
besag = BesagModel(W, additional_constraints=(A_extra, e_extra))
gmrf = besag(τ=1.0)  # Component-wise sum-to-zero + extra constraint

Don't Remove Required Constraints

Intrinsic models like RW1 and Besag require their built-in constraints for identifiability. Attempting to use constraint=:sumtozero on RW1Model will throw an error since it already has this constraint. Always use additional_constraints for these models.

Parameter Validation

All models include built-in parameter validation:

ar1 = AR1Model(100)
ar1(τ=-1.0, ρ=0.8)   # ArgumentError: τ must be positive
ar1(τ=1.0, ρ=1.5)    # ArgumentError: |ρ| must be < 1

Efficient Sparse Operations

• Block-diagonal precision matrices are always sparse

• Input matrix structures (sparse, SymTridiagonal) are preserved

• Minimal memory allocations through size caching

See Also

• GMRF and ConstrainedGMRF for the underlying GMRF types

• Hard Constraints for constraint handling details

• Observation Models for linking GMRFs to data

Formula Terms

Prefer writing models with formulas? The latent components discussed here can be constructed via formula terms and assembled automatically into a combined model and design matrix. See the Formula Interface reference for details:

• Formula Interface — terms (IID, RandomWalk, AR1, Besag, Matern) and build_formula_components.

• For a worked example combining Besag + IID + fixed effects under a Poisson likelihood with offset, see the tutorial: Advanced GMRF modelling for disease mapping