Package 'ngme2'

Title: Linear Latent Non-Gaussian Models with Flexible Distributions
Description: Fits and analyzes linear latent non-Gaussian models for temporal, spatial, and space-time data. The package provides model components for autoregressive and Ornstein-Uhlenbeck processes, random walks, Matern fields based on stochastic partial differential equations, separable and non-separable space-time models, graph-based Matern models, bivariate type-G fields, and user-defined sparse operators. Latent fields and observation models can use Gaussian and non-Gaussian noise distributions, including normal inverse Gaussian, generalized asymmetric Laplace, and skew-t distributions. Functions are included for simulation, likelihood-based estimation, prediction, cross-validation, convergence diagnostics, stochastic gradient optimization, batch-means confidence intervals, and posterior-like sampling. The modeling framework is described in Bolin, Jin, Simas and Wallin (2026) "A Unified and Computationally Efficient Non-Gaussian Statistical Modeling Framework" <doi:10.48550/arXiv.2602.23987>.
Authors: David Bolin [aut, cph], Xiaotian Jin [aut, cre], Alexandre Simas [aut], Jonas Wallin [aut], Andrea V. Rocha [ctb] (contributed to many vignettes and examples), Timothy A. Davis [ctb, cph] (SuiteSparse components), Patrick R. Amestoy [ctb, cph] (SuiteSparse AMD and CAMD components), Iain S. Duff [ctb, cph] (SuiteSparse AMD and CAMD components), John K. Reid [ctb, cph] (SuiteSparse AMD-derived code), Yanqing Chen [ctb, cph] (SuiteSparse CAMD components), Sivasankaran Rajamanickam [ctb, cph] (SuiteSparse CCOLAMD components), Stefan Larimore [ctb, cph] (SuiteSparse CCOLAMD and COLAMD components), William W. Hager [ctb, cph] (SuiteSparse CHOLMOD Modify components), University of Florida [cph] (SuiteSparse CHOLMOD and CCOLAMD components), Regents of the University of Minnesota [cph] (METIS/GKlib components bundled with SuiteSparse), Free Software Foundation, Inc. [cph] (GKlib regex/getopt code bundled with SuiteSparse), Makoto Matsumoto [ctb, cph] (Mersenne Twister code bundled with SuiteSparse), Takuji Nishimura [ctb, cph] (Mersenne Twister code bundled with SuiteSparse), Alexander Chemeris [ctb, cph] (MS integer headers bundled with SuiteSparse)
Maintainer: Xiaotian Jin <[email protected]>
License: GPL (>= 3)
Version: 0.9.8
Built: 2026-05-19 08:17:09 UTC
Source: https://github.com/davidbolin/ngme2

Help Index

AdaGrad SGD optimization

Description

AdaGrad SGD optimization

Usage

adagrad(stepsize = 0.05, epsilon = 1e-08)

Arguments

stepsize

stepsize for SGD
epsilon

epsilon for numerical stability

Details

The update rule for AdaGrad is:

vt=vt−1+gt2v_t = v_{t-1} + g_t^2

xt+1=xt−stepsize∗gtvt+ϵx_{t+1} = x_t - \text{stepsize} * \frac{g_t}{\sqrt{v_t} + \epsilon}

Value

a list of control variables for optimization (used in control_opt function)

Adam SGD optimization

Description

Adam SGD optimization

Usage

adam(stepsize = 0.05, beta1 = 0.9, beta2 = 0.999, epsilon = 1e-08)

Arguments

stepsize

stepsize for SGD
beta1

beta1 for Adam
beta2

beta2 for Adam
epsilon

epsilon for numerical stability

Details

The update rule for Adam is:

mt=β1mt−1+(1−β1)gtm_t = \beta_1 m_{t-1} + (1 - \beta_1) g_t

vt=β2vt−1+(1−β2)gt2v_t = \beta_2 v_{t-1} + (1 - \beta_2) g_t^2

mt^=mt/(1−β1t)\hat{m_t} = m_t / (1 - \beta_1^t)

vt^=vt/(1−β2t)\hat{v_t} = v_t / (1 - \beta_2^t)

xt+1=xt−stepsize∗mt^vt^+ϵx_{t+1} = x_t - \text{stepsize} * \frac{\hat{m_t}}{\sqrt{\hat{v_t}} + \epsilon}

Value

a list of control variables for optimization (used in control_opt function)

AdamW SGD optimization

Description

AdamW SGD optimization

Usage

adamW(
  stepsize = 0.05,
  beta1 = 0.9,
  beta2 = 0.999,
  lambda = 0.01,
  epsilon = 1e-08
)

Arguments

stepsize

stepsize for SGD
beta1

beta1 for AdamW
beta2

beta2 for AdamW
lambda

lambda (weight decay) for AdamW
epsilon

epsilon for numerical stability

Details

The update rule for AdamW is:

mt=β1mt−1+(1−β1)gtm_t = \beta_1 m_{t-1} + (1 - \beta_1) g_t

vt=β2vt−1+(1−β2)gt2v_t = \beta_2 v_{t-1} + (1 - \beta_2) g_t^2

mt^=mt/(1−β1t)\hat{m_t} = m_t / (1 - \beta_1^t)

vt^=vt/(1−β2t)\hat{v_t} = v_t / (1 - \beta_2^t)

xt+1=xt−stepsize∗(λxt+mt^vt^+ϵ)x_{t+1} = x_t - \text{stepsize} * \left( \lambda x_t + \frac{\hat{m_t}}{\sqrt{\hat{v_t}} + \epsilon} \right)

Value

a list of control variables for optimization (used in control_opt function)

Adaptive gradient descent

Description

Adaptive gradient descent optimizer.

Usage

adaptive_gd(stepsize = 0.01)

Arguments

stepsize

initial stepsize for SGD

Details

Based on the method described in https://arxiv.org/pdf/1910.09529. The update rule for adaptive gradient descent is:

λk=min⁡(1+θk−1λk−1,∣∣xk−xk−1∣∣2∣∣∇f(xk)−∇f(xk−1)∣∣)\lambda_k = \min(\sqrt{1 + \theta_{k-1}} \lambda_{k-1}, \frac{||x_k - x_{k-1}||}{2 ||\nabla f(x_k) - \nabla f(x_{k-1})||} )

xk+1=xk−λk∇f(xk)x_{k+1} = x_k - \lambda_k \nabla f(x_k)

θk=λk/λk−1\theta_k = \lambda_k / \lambda_{k-1}

Value

a list of control variables for optimization (used in control_opt function)

ngme AR(p) model specification

Description

ngme AR(p) model specification

Usage

ar(mesh = NULL, rho = NULL, order = NULL)

Arguments

mesh

integer vector or inla.mesh.1d object, index to build the mesh
rho

vector of AR coefficients of length p (should satisfy stationarity conditions)
order

integer, the AR order p. If provided and rho is not specified, rho will be initialized as rep(0, p)

Value

ngme_operator object

Examples

# AR(2) model with specified coefficients
ar(c(1:10), rho = c(0.5, -0.3))
# AR(3) model with specified coefficients
ar(c(1:10), rho = c(0.4, 0.2, -0.1))
# AR(2) model with default coefficients (zeros)
ar(c(1:10), order = 2)
# AR(3) model with specified order and coefficients
ar(c(1:10), order = 3, rho = c(0.4, 0.2, -0.1))

ngme AR(1) model specification

Description

ngme AR(1) model specification

Usage

ar1(mesh = NULL, rho = 0)

Arguments

mesh

integer vector or inla.mesh.1d object, index to build the mesh
rho

the correlation parameter (between -1 and 1)

Value

ngme_operator object

Argo float dataset

Description

Argo floats measurements.

Usage

data("argo_float")

Format

Data frame containing 274 observations on 4 variables.

lat

Latitude.

lon

Longitude.

sal

Salinity.

temp

Temperature.

Details

The floats have a pressure case made of aluminium that is about 1.3m long and about 20cm diameter. They weigh about 40kg. On the top is an antenna to communicate with the satellites that fix the float's position and receive the data. Also on the top are the temperature, salinity and pressure sensors.

Source

Data can be obtained from Argo float website.

ngme ARMA(p, q) model specification

Description

General ARMA(p, q) latent operator under AZW: K W = Z eps, where K encodes the AR part and Z encodes the MA part.

Usage

arma(
  mesh = NULL,
  ar_order = 1,
  ma_order = 1,
  ar = NULL,
  ma = NULL,
  fix_ar = FALSE,
  fix_ma = FALSE
)

Arguments

mesh

Integer vector or 'inla.mesh.1d' object.
ar_order

Integer p (0, 1, or 2).
ma_order

Integer q (0, 1, or 2).
ar

Numeric length-p vector of AR coefficients in (-1, 1). Defaults to zeros.
ma

Numeric length-q vector of MA coefficients in (-1, 1). Defaults to zeros.
fix_ar

Logical or length-p logical vector to fix AR coefficients.
fix_ma

Logical or length-q logical vector to fix MA coefficients.

Details

- Supports p, q <= 2. The user provides standard AR/MA coefficients (ar, ma). For estimation, these are transformed in R into an unbounded parameter vector 'theta_K' via PACF-style mappings so optimizers (e.g., Adam/SGD) operate on unconstrained variables while K and Z remain valid: - AR(1): phi1 = tanh(raw1/2). AR(2): pi1=tanh(raw1/2), pi2=tanh(raw2/2), phi2 = pi2, phi1 = pi1 * (1 - pi2). - MA(1): theta1 = tanh(raw1/2). MA(2): psi1=tanh(raw1/2), psi2=tanh(raw2/2), theta2 = psi2, theta1 = psi1 * (1 - psi2). - The returned operator includes both K (AR part) and Z (MA part). During estimation the backend updates Z synchronously every time 'theta_K' changes.

Value

An 'ngme_operator' describing the ARMA(p, q) latent model.

Examples

mesh <- 1:10
op <- arma(mesh, ar_order = 2, ma_order = 2, ar = c(0.5, 0.3), ma = c(0.6, 0.4))

Convenience wrapper for ARMA(1,1)

Description

Convenience wrapper for ARMA(1,1)

Usage

arma11(mesh = NULL, ar = 0, ma = 0, fix_ar = FALSE, fix_ma = FALSE)

Arguments

mesh

Integer vector or 'inla.mesh.1d' object.
ar

Numeric length-p vector of AR coefficients in (-1, 1). Defaults to zeros.
ma

Numeric length-q vector of MA coefficients in (-1, 1). Defaults to zeros.
fix_ar

Logical or length-p logical vector to fix AR coefficients.
fix_ma

Logical or length-q logical vector to fix MA coefficients.

Value

An ngme_operator object describing an ARMA(1,1) latent operator. The object contains the AR precision component K, the MA filter component Z, integration weights h, and the unconstrained AR and MA parameters used by the optimizer. If mesh is NULL, returns an ngme_operator_def specification for delayed construction inside f().

Batch/checkpoint decay helper

Description

Batch/checkpoint decay helper

Usage

batch_decay(
  patience = 3,
  gamma = 0.5,
  min_delta = 0,
  warmup = 0,
  min_stepsize = 0
)

Arguments

patience

number of consecutive checkpoints without improvement before decay.
gamma

decay factor in (0,1)(0,1).
min_delta

minimum required decrease to count as improvement.
warmup

number of initial checkpoints to skip.
min_stepsize

lower bound for absolute stepsize.

Details

Convenience helper that enables grad-norm plateau decay and keeps iteration schedule constant.

Value

a stepsize_control() object.

Pooled Batch-Means Confidence Intervals from Multiple Chains

Description

Compute pooled point estimates, covariance, and Wald confidence intervals from multiple chain trajectories.

Usage

batch_means_ci(
  chain_trajectories,
  level = 0.95,
  alpha = 0.501,
  M = NULL,
  N = NULL,
  drop_burnin = TRUE,
  burnin_iter = 0
)

Arguments

chain_trajectories

list of numeric matrices, each with rows = iterations and columns = parameters.
level

confidence level.
alpha

stepsize decay exponent in (1/2,1)(1/2, 1).
M

number of retained batches per chain (excluding burn-in batch 0).
N

decorrelation constant used in the batch boundary formula.
drop_burnin

logical; if 'TRUE', discard batch 0.
burnin_iter

non-negative integer. Explicitly discard the first 'burnin_iter' iterations of each chain before Xi-style batching.

Value

A list with pooled estimates, covariance, standard errors, and confidence intervals.

Batch-Means Covariance Estimator for SGD Trajectories

Description

Estimate the asymptotic covariance from a single SGD trajectory using the increasing-batch construction from Xi et al. (2020).

Usage

batch_means_estimator(
  trajectory,
  alpha = 0.501,
  M = NULL,
  N = NULL,
  drop_burnin = TRUE,
  burnin_iter = 0
)

Arguments

trajectory

numeric matrix with rows = iterations and columns = parameters.
alpha

stepsize decay exponent in ηi=ηi−α\eta_i = \eta i^{-\alpha}. Must satisfy 1/2<α<11/2 < \alpha < 1.
M

number of retained batches (excluding burn-in batch 0). If 'NULL', use ⌊n(1−α)/2⌋\lfloor n^{(1-\alpha)/2} \rfloor.
N

decorrelation constant in ek=⌊((k+1)N)1/(1−α)⌋e_k = \lfloor ((k+1)N)^{1/(1-\alpha)} \rfloor. If 'NULL', use N=n1−α/(M+1)N = n^{1-\alpha}/(M+1).
drop_burnin

logical; if 'TRUE', discard batch 0.
burnin_iter

non-negative integer. Explicitly discard the first 'burnin_iter' iterations before building Xi-style batches. After trimming, batch boundaries are rebuilt from iteration 1 of the retained trajectory.

Value

A list containing the covariance estimate, pooled mean, batch sizes, batch boundaries, and chain-level metadata.

BFGS optimization

Description

BFGS optimization

Usage

bfgs(line_search = "wolfe")

Arguments

line_search

line search method, can be c("backtracking", "wofle") "bfgs" means use BFGS preconditioner

Value

a list of control variables for optimization (used in control_opt function)

Ngme bivariate model specification

Description

Giving 2 sub_models, build a correlated bivariate operator based on K = D(theta, rho)

D(θ,ρ)=(cos⁡(θ)+ρsin⁡(θ)−sin⁡(θ)1+ρ2sin⁡(θ)−ρcos⁡(θ)cos⁡(θ)1+ρ2)D(\theta, \rho) = \begin{pmatrix} \cos(\theta) + \rho \sin(\theta) & -\sin(\theta) \sqrt{1+\rho^2} \\ \sin(\theta) - \rho \cos(\theta) & \cos(\theta) \sqrt{1+\rho^2} \end{pmatrix}

Usage

bv(
  mesh,
  sub_models,
  theta = 0,
  rho = 0,
  c1 = 1,
  c2 = 1,
  share_param = FALSE,
  fix_theta = FALSE,
  use_c_param = FALSE
)

Arguments

mesh

the mesh where model is defined
sub_models

a list of sub_models (total 2 sub_models)
theta

the rotation parameter (starting point should be from -pi/4 to pi/4)
rho

the parameter related to correlation
c1

scaling factor for the first sub-model (default: 1)
c2

scaling factor for the second sub-model (default: 1)
share_param

TRUE if share the same parameter for 2 sub_models (of same type)
fix_theta

TRUE if fix the theta of bv model
use_c_param

TRUE to estimate c1 and c2 as parameters, FALSE to fix them at specified values (default: FALSE)

Details

The bivariate model combines two sub-models with a rotation matrix D and optional scaling factors c1 and c2. When use_c_param = TRUE, c1 and c2 are estimated as parameters (on log scale). When use_c_param = FALSE (default), c1 and c2 are fixed at their specified values (both default to 1).

Value

a ngme_operator of bivariate model

Ngme bivariate model with Matern sub_models

Description

Giving 2 sub_models, build a correlated bivaraite operator based on K = D(theta, eta)

D(θ,ρ)=(cos⁡(θ)+ρsin⁡(θ)−sin⁡(θ)1+ρ2sin⁡(θ)−ρcos⁡(θ)cos⁡(θ)1+ρ2)D(\theta, \rho) = \begin{pmatrix} \cos(\theta) + \rho \sin(\theta) & -\sin(\theta) \sqrt{1+\rho^2} \\ \sin(\theta) - \rho \cos(\theta) & \cos(\theta) \sqrt{1+\rho^2} \end{pmatrix}

Usage

bv_matern(
  mesh,
  sub_models,
  theta = 0,
  rho = 0,
  sd1 = 1,
  sd2 = 1,
  group = NULL,
  share_param = FALSE,
  fix_theta = FALSE,
  ...
)

Arguments

mesh

the mesh where model is defined
sub_models

a list of sub_models (should be two matern models)
theta

the parameter related to rotation
rho

the parameter related to correlation
sd1

scaling parameter for the first sub-model
sd2

scaling parameter for the second sub-model
group

group vector, can be inherited from ngme() function
share_param

TRUE if share the same parameter for 2 sub_models (of same type)
fix_theta

TRUE if the rotation parameter is fixed
...

ignore

Value

a list of specification of model

Calibrate Inverse-Exponential Prior for NIG Driven Noise

Description

Calibrate λ\lambda in the prior κ=1/ν∼Exp(λ)\kappa = 1/\nu \sim \mathrm{Exp}(\lambda) using a tail-inflation target for driven NIG noise.

The calibration target is

Pr⁡(Rc(ν)>rtarget)=α,\Pr(R_c(\nu) > r_\text{target}) = \alpha,

where

Rc(ν)=Pr⁡(∣U∣>c∣ν)Pr⁡(∣Z∣>c),Z∼N(0,1)R_c(\nu) = \frac{\Pr(|U| > c \mid \nu)}{\Pr(|Z| > c)},\quad Z\sim N(0,1)

and

U=μ(V−h)+σVZσh,V∼GIG(−1/2,ν,νh2).U = \frac{\mu(V-h) + \sigma\sqrt{V}Z}{\sigma\sqrt{h}},\quad V\sim\mathrm{GIG}(-1/2,\nu,\nu h^2).

The function solves Rc(νr)=rtargetR_c(\nu_r) = r_\text{target} using log-scale bisection, then returns

λ=−νrlog⁡(α).\lambda = -\nu_r \log(\alpha).

Usage

calibrate_inv_exp_lambda_driven_nig(
  r_target = 2,
  alpha = 0.1,
  c = 3,
  mu = 0,
  sigma = 1,
  h = 1,
  n_samples = 1e+05,
  nu_lower = 0.1,
  nu_upper = 100,
  tol = 0.05,
  max_iter = 30,
  max_expand = 30,
  seed = NULL
)

calibrate_inv_exp_lambda(...)

Arguments

r_target

target tail inflation level rtarget>1r_\text{target} > 1
alpha

target prior probability in (0,1)(0,1)
c

tail threshold for ∣U∣>c|U|>c, typically '2.5' or '3'
mu

NIG drift parameter in the driven noise term
sigma

NIG scale parameter (> 0)
h

positive increment scaling (> 0)
n_samples

Monte Carlo sample size used per evaluation of Rc(ν)R_c(\nu)
nu_lower

initial lower bracket for ν\nu
nu_upper

initial upper bracket for ν\nu
tol

relative tolerance for solving Rc(νr)=rtargetR_c(\nu_r)=r_\text{target}
max_iter

maximum number of bisection iterations
max_expand

maximum bracket expansion steps on each side
seed

optional integer seed for reproducible calibration
...

arguments forwarded to calibrate_inv_exp_lambda_driven_nig()

Value

A list with calibrated 'lambda', solved 'nu_r', achieved 'rc_nu_r', and diagnostics.

The swamp of Cienaga Grande in Santa Marta, Colombia

Description

There is a total of 114 locations where some properties of the swamp were measured. Those measurements were taken twice, however there is no information available about their chronological order so this data cannot be treated as spatiotemporal, despite that, the multiple measurements can be treated as replicates.

Usage

cienaga

Format

A data frame with 218 rows and 6 columns.

East, North

location

depth

depth of the swamp

temp

temperature

oxyg

oxygen

measurement

1 means the first measurement, 2 the second

Source

..

The x y location of the border of the swamp of Cienaga Grande in Santa Marta, Colombia

Description

The data is of dimension 472 * 2. It contains the x and y coordinates of the border of the swamp.

Usage

cienaga.border

Format

A data frame with 472 rows and 2 columns.

East, North

location

Source

..

Compare noise objects using Kullback-Leibler divergence

Description

This function compares multiple noise objects by calculating the KLD of each noise against the first noise object (reference).

Usage

compare_noise_kld(x = NULL, ..., xlim = c(-10, 10), n_points = 1000)

Arguments

x

first noise object (reference)
...

additional noise objects to compare, and optional parameters
xlim

x-axis range for evaluation (default: c(-10, 10))
n_points

number of evaluation points (default: 1000)

Value

named vector of KLD values

Examples

n1 <- noise_nig(mu = 0, sigma = 1, nu = 1)
n2 <- noise_nig(mu = 0.5, sigma = 1.2, nu = 0.8)
compare_noise_kld(n1, method2 = n2)

Helper function to compute the index_corr vector

Description

Helper function to compute the index_corr vector

Usage

compute_index_corr_from_map(map, eps = 0.1)

Arguments

map

used as location to compute distance, can be 1d (numerical) or 2d (data.frame)
eps

threshold to determine if two points are close (if close, we consider them as the same point)

Value

the index_corr vector for ngme correlated measurement noise

Examples

x_coord <- c(1.11, 1.12, 2, 1.3, 1.3)
y_coord <- c(2.11, 2.11, 2, 3.3, 3.3)
coord <- data.frame(x_coord, y_coord)
compute_index_corr_from_map(map = coord, 0.1)

Compute Gaussian log-likelihood

Description

Creates the underlying C++ block models and evaluates their Gaussian log-likelihood when applicable.

Usage

compute_log_like(ngme_model)

Arguments

ngme_model

An object created by 'ngme()'.

Value

A numeric scalar. Returns '0' if any replicate is non-Gaussian.

Refit an Existing ngme Object with SGD and Compute Batch-Means CI

Description

Run one additional SGD stage (warm-started from an existing 'ngme' fit), then compute Xi-style batch-means confidence intervals from the newly stored trajectories.

Usage

compute_ngme_ci(
  fit,
  iterations = 4000,
  alpha = 0.501,
  optimizer = sgd(stepsize = 0.03),
  burnin = 100,
  n_batch = 20,
  n_parallel_chain = 4,
  t0 = 1,
  start_sd = 0.2,
  seed = Sys.time(),
  verbose = FALSE,
  level = 0.95,
  name = "all",
  M = NULL,
  N = NULL,
  apply_transform = TRUE,
  drop_burnin = TRUE,
  burnin_iter = 0,
  control_opt = NULL,
  ...
)

compute_ngme_CI(
  fit,
  iterations = 4000,
  alpha = 0.501,
  optimizer = sgd(stepsize = 0.03),
  burnin = 100,
  n_batch = 20,
  n_parallel_chain = 4,
  t0 = 1,
  start_sd = 0.2,
  seed = Sys.time(),
  verbose = FALSE,
  level = 0.95,
  name = "all",
  M = NULL,
  N = NULL,
  apply_transform = TRUE,
  drop_burnin = TRUE,
  burnin_iter = 0,
  control_opt = NULL,
  ...
)

Arguments

fit

existing fitted 'ngme' object (can be obtained with any optimizer).
iterations

optimization iterations for the SGD CI stage.
alpha

Xi-style exponent used for both polynomial schedule and batch-means inference.
optimizer

optimizer for the CI stage; must be 'sgd(...)' for Xi theory.
burnin

burn-in iterations before SGD optimization.
n_batch

number of optimization checkpoints.
n_parallel_chain

number of parallel chains for SGD.
t0

non-negative schedule offset in ηt=η0(t+t0)−α\eta_t = \eta_0 (t + t0)^{-\alpha}.
start_sd

standard deviation for randomized chain initialization.
seed

random seed for CI-stage optimization.
verbose

logical; print optimization progress.
level

confidence level for CI.
name

parameter block for CI ('"all"', latent name, or '"general"').
M

number of retained Xi batches (excluding batch 0).
N

decorrelation constant in Xi batch boundaries.
apply_transform

logical; if 'TRUE', apply Post-delta transform after raw-space inference.
drop_burnin

logical; if 'TRUE', discard Xi batch 0.
burnin_iter

non-negative integer used both as optimizer schedule warmup ('stepsize_schedule_burnin_iter') and as explicit CI trajectory trimming before Xi-style batching.
control_opt

optional pre-built 'control_opt' object for the CI stage. If supplied, it is used directly (with 'store_traj' forced to 'TRUE').
...

additional arguments forwarded to [control_opt_batch_ci()] when 'control_opt' is 'NULL'.

Value

An object of class 'ngme_batch_ci' with extra fields: 'refit' (the SGD-refit model), 'refit_control_opt', and 'refit_source'.

Refit an Existing ngme Object with SGLD and Extract Samples

Description

Run one additional SGLD stage (warm-started from an existing 'ngme' fit), then extract posterior-like samples from stored trajectories.

Usage

compute_ngme_sgld_samples(
  fit,
  iterations = 4000,
  optimizer = sgld(stepsize = 0.01, temperature = 1),
  burnin = 100,
  n_batch = 20,
  n_parallel_chain = 4,
  alpha = 0.6,
  t0 = 10,
  start_sd = 0.2,
  seed = Sys.time(),
  verbose = FALSE,
  name = "all",
  burnin_iter = 0,
  thinning = 1,
  n_gibbs_samples = NULL,
  apply_transform = TRUE,
  combine_chains = TRUE,
  control_opt = NULL,
  ...
)

Arguments

fit

existing fitted 'ngme' object (can be obtained with any optimizer).
iterations

optimization iterations for the SGLD stage.
optimizer

optimizer for the sampling stage; must be 'sgld(...)'.
burnin

burn-in iterations before optimization.
n_batch

number of optimization checkpoints.
n_parallel_chain

number of parallel chains.
alpha

polynomial schedule exponent used by 'poly_decay(alpha, t0)'.
t0

non-negative schedule offset.
start_sd

standard deviation for randomized chain initialization.
seed

random seed for SGLD stage.
verbose

logical; print optimization progress.
name

parameter block to extract: '"all"' (default), latent model name, or '"general"'.
burnin_iter

non-negative integer used both as optimizer schedule warmup ('stepsize_schedule_burnin_iter') and as explicit trajectory trimming before sampling.
thinning

positive integer thinning interval.
n_gibbs_samples

optional positive integer. If supplied, override 'control_ngme$n_gibbs_samples' for the SGLD refit stage only; otherwise inherit the value from 'fit'.
apply_transform

logical; apply parameter transforms to user scale.
combine_chains

logical; if 'TRUE', return one combined data.frame, otherwise return one data.frame per chain.
control_opt

optional pre-built 'control_opt' object for the SGLD stage. If supplied, it is used directly (with 'store_traj' forced to 'TRUE').
...

additional arguments forwarded to [control_opt()] when 'control_opt' is 'NULL'.

Value

A data.frame (or list of data.frames) from [ngme_sgld_samples()] with extra attributes 'refit', 'refit_control_opt', and 'refit_source'.

Compute the scores given the prediction

Description

Compute the scores given the prediction

Usage

compute_score_given_pred(
  Y_N_1_thin,
  Y_N_2_thin,
  y_data,
  group_data,
  merge_groups = FALSE,
  merged_group_name = NULL,
  metric = NULL
)

Arguments

Y_N_1_thin

posterior predictive draws (rows = observations, columns = samples)
Y_N_2_thin

posterior predictive draws (rows = observations, columns = samples)
y_data

a vector of length n_obs
group_data

a vector of length n_obs
merge_groups

logical, if TRUE and there are exactly two groups, combine them into a vector-valued score
merged_group_name

optional label for the merged group
metric

optional custom metric function used to combine group-wise values before scoring

Generate control specifications for the ngme model

Description

Generate control specifications for the ngme model

Usage

control_ngme(
  n_gibbs_samples = 5,
  fix_beta = FALSE,
  n_post_samples = 100,
  beta_init = NULL,
  debug = FALSE,
  ...
)

Arguments

n_gibbs_samples

number of gibbs samples at each iteration
fix_beta

logical, fix fixed effect estimation
n_post_samples

number of posterior samples, see ?ngme_post_samples()
beta_init

fixed effect initial value on the original design scale. If control_opt(standardize_fixed = TRUE) (the default), the vector is internally rotated/scaled to match the SVD-standardized design matrix before being passed to the optimizer. Provide values on the raw design scale; no manual rescaling is required.
debug

debug mode
...

additional arguments. Legacy aliases feff and fix_feff are still recognized and mapped to beta_init and fix_beta.

Value

a list of control variables for ngme

Generate control specifications for ngme() function.

Description

These are configurations for ngme optimization process.

Usage

control_opt(
  seed = Sys.time(),
  burnin = 100,
  iterations = 500,
  estimation = TRUE,
  standardize_fixed = TRUE,
  n_batch = 10,
  iters_per_check = iterations/n_batch,
  optimizer = adam(),
  start = NULL,
  start_sd = 0.5,
  n_parallel_chain = 4,
  max_num_threads = n_parallel_chain,
  print_check_info = FALSE,
  max_relative_step = 0.5,
  max_absolute_step = 0.5,
  rao_blackwellization = FALSE,
  n_trace_iter = 10,
  sampling_strategy = "all",
  solver_backend = if (Sys.info()["sysname"] == "Darwin") "accelerate" else "cholmod",
  solver_type = "llt",
  verbose = FALSE,
  store_traj = TRUE,
  robust = FALSE,
  stepsize_control = NULL,
  n_min_batch = min(n_batch, 3),
  n_slope_check = min(n_batch, 3),
  trend_std_conv_check = TRUE,
  std_lim = 0.01,
  trend_lim = 0.01,
  R_hat_conv_check = TRUE,
  max_R_hat = 1.1,
  pflug_conv_check = TRUE,
  pflug_alpha = 0.9
)

Arguments

seed

set the seed for pesudo random number generator
burnin

interations for burn-in periods (before optimization)
iterations

optimization iterations
estimation

run the estimation process (call C++ in backend)
standardize_fixed

whether or not standardize the fixed effect. When TRUE (default), the design matrix is SVD-standardized internally and any user-supplied control_ngme(beta_init = ...) is automatically mapped from the original design scale to that standardized basis. Set to FALSE to keep both the design matrix and any provided beta_init on their original scale.
n_batch

number of checkpoints; optimization is split into n_batch equal batches
iters_per_check

run how many iterations between each check point (or specify n_batch)
optimizer

choose different sgd optimization method, currently support "sgd", "sgld", "precond_sgd", "momentum", "adagrad", "rmsprop", "adam", "adamW" see ?sgd, ?precond_sgd, ?momentum, ?adagrad, ?rmsprop, ?adam, ?adamW
start

deprecated guard argument. Do not pass model starts through control_opt(); use ngme(..., start = previous_fit) instead.
start_sd

standard deviation of the initial parameter (1st chain fixed, other chains random), set 0 to be fixed for all chains
n_parallel_chain

number of parallel chains
max_num_threads

maximum number of threads used for parallel computing, by default will be set same as n_parallel_chain. If it is more than n_parallel_chain, the rest will be used to parallel different replicates of the model.
print_check_info

print the convergence information
max_relative_step

max relative step allowed in 1 iteration
max_absolute_step

max absolute step allowed in 1 iteration
rao_blackwellization

use rao_blackwellization
n_trace_iter

use how many iterations to approximate the trace (Hutchinson’s trick)
sampling_strategy

subsampling method of replicates of model, c("all", "ws") "all" means using all replicates in each iteration, "ws" means weighted sampling (each iteration use 1 replicate to compute the gradient, the sample probability is proption to its number of observations)
solver_backend

backend in ("eigen", "cholmod", "accelerate", "pardiso")
solver_type

factorization type: "llt" or "ldlt"
verbose

print estimation
store_traj

store the optimizer trajectory for diagnostics (set FALSE to reduce memory)
robust

use robust mode in the backend optimizer/model updates
stepsize_control

unified stepsize configuration created by stepsize_control(), poly_decay(), or batch_decay(). For polynomial schedule, poly_decay(..., burnin_iter = B) keeps schedule scale at 1 for the first B iterations, then starts decay with reset local time index.
n_min_batch

minimum number of checkpoints before any convergence diagnostic is attempted
n_slope_check

number of checkpoints used as the regression window for the trend test
trend_std_conv_check

enable the trend/std diagnostic (uses std_lim, trend_lim, n_slope_check)
std_lim

maximum allowed standard deviation
trend_lim

maximum allowed slope
R_hat_conv_check

use the R-hat diagnostic for convergence checking
max_R_hat

maximum allowed R_hat
pflug_conv_check

use Pflug diagnostic for convergence check
pflug_alpha

scaling factor (0-1] for Pflug criterion: require pflug_sum < pflug_alpha * max_pflug_sum

Details

Convergence diagnostics (multi-chain): * R-hat: per-parameter Gelman–Rubin statistic; passes if R_hat <= max_R_hat. * Trend/Std: uses the last n_slope_check checkpoints after at least n_min_batch batches. Passes when both the relative std (sqrt(var)/|mean| <= std_lim) and linear trend of the means (|slope| <= trend_lim) satisfy their thresholds. * Pflug: per-chain criterion pflug_sum < pflug_alpha * max_pflug_sum in the latest batch; if all chains satisfy it, overall convergence is declared. Checks are evaluated every iters_per_check = iterations / n_batch. A parameter is marked converged if any enabled parameter-level diagnostic (R-hat or Trend/Std) passes; the run stops when all parameters converge or when the Pflug diagnostic triggers.

Value

list of control variables

Generate CI-focused control settings for batch-means inference

Description

Convenience wrapper for control_opt() with defaults tailored for Xi-style batch-means confidence intervals.

Usage

control_opt_batch_ci(
  optimizer = sgd(stepsize = 0.03),
  burnin = 100,
  iterations = 2000,
  n_batch = 20,
  n_parallel_chain = 4,
  alpha = 0.501,
  t0 = 1,
  schedule_burnin_iter = 0,
  ...
)

Arguments

optimizer

optimizer object, default sgd(stepsize = 0.03).
burnin

burn-in iterations before optimization.
iterations

optimization iterations.
n_batch

number of checkpoints.
n_parallel_chain

number of parallel chains.
alpha

polynomial stepsize exponent for poly_decay(alpha, t0).
t0

non-negative schedule offset.
schedule_burnin_iter

non-negative integer. Initial optimization iterations where polynomial schedule scaling is disabled.
...

additional arguments forwarded to control_opt() to override defaults.

Details

This helper enforces trajectory-friendly defaults:

  • store_traj = TRUE

  • trend_std_conv_check = FALSE

  • R_hat_conv_check = FALSE

  • pflug_conv_check = FALSE

  • stepsize_control = poly_decay(alpha, t0, schedule_burnin_iter)

Any of these can still be overridden through ....

Value

object of class control_opt.

Create paired indices for bivariate cross-validation Ensures that paired observations (e.g., u_wind and v_wind at same location) are kept together in the same fold

Description

Create paired indices for bivariate cross-validation Ensures that paired observations (e.g., u_wind and v_wind at same location) are kept together in the same fold

Usage

create_paired_cv_splits(data, loc_col, group, k, seed = NULL)

Arguments

data

data frame containing the observations
loc_col

character vector of column names defining locations (e.g., c("lon", "lat") or c("x", "y"))
group

character, name of the group column (e.g., "direction" for wind data)
k

number of folds
seed

random seed

Value

list with test_idx and train_idx, each containing k lists of indices

Examples

# Example with wind data
data_long <- data.frame(
  lon = rep(c(1, 2, 3, 4), 2),
  lat = rep(c(1, 1, 2, 2), 2),
  direction = rep(c("u_wind", "v_wind"), each = 4),
  wind = rnorm(8)
)
splits <- create_paired_cv_splits(data_long, c("lon", "lat"), "direction", k = 2)

Compute the cross-validation for the ngme model Perform cross-validation for ngme model first into sub_groups (a list of target, and train data)

Description

Compute the cross-validation for the ngme model Perform cross-validation for ngme model first into sub_groups (a list of target, and train data)

Usage

cross_validation(
  ngme,
  type = "k-fold",
  seed = NULL,
  print = TRUE,
  N_sim = 5,
  n_gibbs_samples = 500,
  n_burnin = 100,
  k = 5,
  percent = 0.2,
  times = 10,
  metric = NULL,
  test_idx = NULL,
  train_idx = NULL,
  keep_pred = FALSE,
  parallel = FALSE,
  thining_gap = 1,
  cores_layer1 = if (parallel) min(parallel::detectCores(), 2) else 1,
  cores_layer2 = if (parallel) min(parallel::detectCores(), 2) else 1,
  merge_groups = FALSE,
  merged_group_name = NULL,
  data = NULL,
  chain_combine = c("param_mean", "predictive_average")
)

Arguments

ngme

a ngme object, or a list of ngme object (if comparing multiple models)
type

character, in c("k-fold", "loo", "lpo", "custom") k-fold is k-fold cross-validation, provide k loo is leave-one-out, lpo is leave-percent-out, provide percent from 1 to 100 custom is user-defined group, provide target and data
seed

random seed
print

print information during computation
N_sim

integer, number of simulations (e.g., estimate MAE, MSE, .. N times)
n_gibbs_samples

number of gibbs samples of latent process, used for computing CRPS, sCRPS
n_burnin

number of burnin
k

integer (only for k-fold type)
percent

how many percent for testing? from 0 to 1 (for lpo type)
times

how many test cases (only for lpo type)
metric

Optional function or list of functions (one per model) that maps the group-wise observations/predictions for a single location to the quantity that should be scored. The function receives a list containing at least y (a named numeric vector of group values) and may optionally use samples1 and samples2 (matrices with rows named by group and columns indexing posterior draws). The function must return either a numeric scalar (optionally named), or a list with component y and optional components samples1, samples2, and label. When NULL, the original per-group scores are computed. Example: metric = function(data) 2 * data$y["A"] + data$y["B"]. To sum all groups, return sum(data$y).
test_idx

a list of indices of the data (which data points to be predicted) (only for custom type)
train_idx

a list of indices of the data (which data points to be used for re-sampling (not re-estimation)) (only for custom type)
keep_pred

logical, keep test information (pred_1, pred_2) in the return (as attributes), pred_1 and pred_2 are the prediction of the two chains
parallel

logical, run in parallel mode
thining_gap

integer, the gap between samples for thinning, if 0, then no thinning, if 1, then keep 50% of the samples for CRPS, sCRPS, etc.
cores_layer1

integer, number of cores for the first layer (over testing samples)
cores_layer2

integer, number of cores for the second layer (over computing scores for N_sim simulations)
merge_groups

logical, if TRUE, merge groups as vector components (e.g., for vector-valued wind data with north_wind, east_wind). MAE becomes Euclidean distance, MSE becomes squared Euclidean distance, etc.
merged_group_name

character, name for the merged group when merge_groups=TRUE. If NULL, uses "group1_group2" format (default: NULL)
data

optional data.frame used to replace the original fitting data before running CV. If 'NULL', the data stored in 'ngme' is used. If provided, the model is rebuilt on 'data' while reusing fitted parameters from 'ngme'.
chain_combine

how to combine multiple optimization chains: '"param_mean"' uses the fitted object directly (default), while '"predictive_average"' computes predictions from each optimization chain and averages at the predictive level.

Value

A list with components:

mean.scores

Mean of the N_sim estimates for MSE, MAE, CRPS, and sCRPS.

sd.scores

Standard deviation of the N_sim estimates for MSE, MAE, CRPS, and sCRPS.

Specifying a latent process model (wrapper function for each model)

Description

Function used for defining of smooth and spatial terms within ngme model formulae. The function is a wrapper function for specific submodels. (see ngme_models_types() for available models).

Usage

f(
  map,
  model,
  noise = noise_normal(),
  name = "field",
  data = NULL,
  group = NULL,
  which_group = NULL,
  replicate = NULL,
  A = NULL,
  W = NULL,
  fix_W = FALSE,
  fix_K = FALSE,
  prior = NULL,
  subset = rep(TRUE, length_map(map)),
  debug = FALSE
)

Arguments

map

symbol or numerical value: index or covariates to build index
model

an ngme_operator or ngme_operator_def created by model constructors such as ar1(), matern(), bv(), etc.
noise

ngme_noise object, noise_nig() or noise_gal()
name

name of the field, for later use, if not provided, will be "field1" etc.
data

specifed or inherit from ngme() function
group

group factor indicate resposne variable, can be inherited from ngme() function, (used for bivariate model)
which_group

belong to which group
replicate

factor indicating replicate structure for INLA-style replicates. When provided, operators and A matrices will be block-diagonalized across replicates.
A

observation matrix, automatically computed given map and model
W

starting value of the process
fix_W

stop sampling for W
fix_K

fix the estimation of parameter of K
prior

prior specification created by prior_*() or priors(...) for operator parameters (theta_K).
subset

subset of the model
debug

debug mode

Details

When using different meshes for different replicates, provide the mesh parameter as a list of mesh objects. The number of meshes should match the number of replicates. For example: mesh = list(mesh1, mesh2, mesh3) for 3 replicates.

When using the 'replicate' argument, the operator matrices will be block-diagonalized. For a generic operator K = param1 * Matrix1 + param2 * Matrix2, with 2 replicates, the resulting operator becomes: K = param1 * bdiag(Matrix1_rep1, Matrix1_rep2) + param2 * bdiag(Matrix2_rep1, Matrix2_rep2) Similarly, the A matrix will be block-diagonalized as bdiag(A_rep1, A_rep2, ...).

Value

a list for constructing latent model, e.g. A, h, C, G, which also has 1. Information about K matrix 2. Information about noise 3. Control variables

The Generalized Asymmetric Laplace (GAL) Distribution

Description

Density, distribution function, quantile function and random generation for the generalized asymmetric Laplace distribution with parameters mu, sigma and nu, delta.

Usage

dgal(x, delta, mu, nu, sigma, log = FALSE)

rgal(n, delta, mu, nu, sigma, seed = 0)

pgal(q, delta, mu, nu, sigma, lower.tail = TRUE, log.p = FALSE)

qgal(p, delta, mu, nu, sigma, lower.tail = TRUE, log.p = FALSE)

Arguments

x, q

vector of quantiles.
delta

A numeric value for the location parameter.
mu

A numeric value for the shift parameter.
nu

A numeric value for the shape parameter.
sigma

A numeric value for the scaling parameter.
log, log.p

logical; if TRUE, probabilities/densities pp are returned as log(p)log(p).
n

number of observations.
seed

Seed for the random generation.
lower.tail

logical; if TRUE, probabilities are P[X≤x]P[X\leq x], otherwise, P[X>x]P[X>x].
p

vector of probabilities.

Details

The generalized asymmetric Laplace distribution has density given by

f(x;p,a,b)=eν+μ(x−δ)/σ2νμ2/σ2+ν2πνσ2+(x−δ)2K1((νσ2+(x−δ)2)(μ2/σ4+ν/σ2)),f(x; p, a, b) = \frac{e^{\nu+\mu(x-\delta)/\sigma^2}\sqrt{\nu\mu^2/\sigma^2+\nu^2}}{\pi\sqrt{\nu\sigma^2+(x-\delta)^2}} K_1(\sqrt{(\nu\sigma^2+(x-\delta)^2)(\mu^2/\sigma^4+\nu/\sigma^2)}),

where KpK_p is modified Bessel function of the second kind of order pp, x>0x>0, ν>0\nu>0 and μ,δ,σ∈R\mu,\delta, \sigma\in\mathbb{R}. See Barndorff-Nielsen (1977, 1978 and 1997) for further details.

If the mixing variable VV follows a Gamma distribution (same parameterization in R):

V∼Γ(hν,ν),V \sim \Gamma(h \nu, \nu),

then the poserior follows the GAL distribution (a special case of GIG distribution):

−μ+μV+σVZ∼GIG(hν−0.5,2ν+(μσ)2,0)-\mu +\mu V + \sigma \sqrt{V} Z \sim GIG(h \nu - 0.5, 2 \nu + (\frac{\mu}{\sigma})^{2}, 0)

Value

dgal gives the density, pgal gives the distribution function, qgal gives the quantile function, and rgal generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rgal.

References

Barndorff-Nielsen, O. (1977) Exponentially decreasing distributions for the logarithm of particle size. Proceedings of the Royal Society of London.

Series A, Mathematical and Physical Sciences. The Royal Society. 353, 401–409. doi:10.1098/rspa.1977.0041

Barndorff-Nielsen, O. (1978) Hyperbolic Distributions and Distributions on Hyperbolae, Scandinavian Journal of Statistics. 5, 151–157.

See Also

dgig, dig, digam

Examples

rgal(100, delta = 0, mu = 5, sigma = 1, nu = 1)
pgal(0.4, delta = 0, mu = 5, sigma = 1, nu = 1)
qgal(0.8, delta = 0, mu = 5, sigma = 1, nu = 1)
plot(function(x){dgal(x, delta = 0, mu = 5, sigma = 1, nu = 1)}, main =
"generalized asymmetric Laplace density", ylab = "Probability density",
xlim = c(0,10))

Generic precision matrix operator

Description

Creates a flexible precision matrix K by allowing linear combinations of base matrices with parameter-dependent coefficients. The model follows the form: K = sum_i f_i(theta) * matrices_i, where f_i depends on the parameter transformation.

Usage

generic(
  theta_K,
  matrices,
  h,
  trans = NULL,
  mesh = NULL,
  zero_trace = FALSE,
  model = "generic",
  param_name = NULL,
  param_trans = NULL,
  ...
)

Arguments

theta_K

The parameter vector (in real scale). If missing, initialized as double(0)
matrices

A list of matrices to be combined with parameter-dependent coefficients
h

The integration weights vector
trans

Transformation for each parameter. Must be a list where each element is a character vector defining transformations for each matrix. For example, ‘trans=list(kappa=c("exp2", "identity", "null"))' means parameter ’kappa' affects matrix 1 with exp2 transformation, matrix 2 with identity, and doesn't affect matrix 3. Available transformations: "identity", "exp", "exp2", "exp4", "sqrt", "square", "log", "tanh"
mesh

The mesh object
zero_trace

Whether the trace of K should be zero
model

The model type name
param_name

Optional parameter names stored on the returned operator for diagnostics
param_trans

Optional parameter transformations stored on the returned operator for diagnostics
...

Additional arguments (ignored)

Details

This function can be used to represent various models like AR1, Matern, and RW1 in a unified framework.

The 'generic' function provides a flexible way to construct precision matrices by combining existing matrices with transformed parameters. Here's how it works:

1. Each parameter in 'theta_K' can affect one or more matrices in the 'matrices' list 2. The 'trans' list defines how each parameter transforms before multiplying each matrix 3. For each matrix, the coefficients from all parameters are multiplied together 4. The final K is the sum of all transformed matrices

Common use cases:

- **AR1 model**: K = rho * C + G (where rho is between -1 and 1) - **Matern (alpha=2)**: K = kappa^2 * C + G (where kappa > 0) - **Matern (alpha=4)**: K = kappa^4 * C + 2*kappa^2 * G + G * Cinv * G

Value

An object of class ngme_operator. The object stores the sparse precision matrix K, integration weights h, parameter vector theta_K, base matrices, transformation specification trans, and an update_K function that rebuilds K for new parameter values. It is intended for use as the latent model inside f().

Examples

# AR1 model with rho = 0.5
ar1_obj <- ar1(1:10, rho = 0.5)
g <- name2fun("tanh", inv = TRUE)
generic_ar1 <- generic(
  theta_K = c(rho = g(0.5)),  # Transform from real scale
  trans = list(rho = c("tanh", "null")),  # Apply tanh to first matrix
  matrices = list(ar1_obj$C, ar1_obj$G),
  h = ar1_obj$h
)

# Matern model with kappa = 2
mesh <- fmesher::fm_mesh_1d(seq(0, 1, length.out = 10))
matern_obj <- matern(mesh)
log_kappa <- log(2)
generic_matern <- generic(
  theta_K = c(kappa = log_kappa),
  trans = list(kappa = c("exp2", "null")),  # exp(2*kappa) for C, nothing for G
  matrices = list(matern_obj$C, matern_obj$G),
  h = matern_obj$h
)

# Matern model with alpha = 4 and kappa = 2
set.seed(1)
mesh <- fmesher::fm_mesh_2d(cbind(x = runif(20), y = runif(20)))
matern_obj <- matern(mesh, alpha = 4)
C <- matern_obj$C
G <- matern_obj$G
Cinv <- C; diag(Cinv) <- 1 / Matrix::diag(C)

generic_matern4 <- generic(
  theta_K = c(kappa = log(2)),
  trans = list(kappa = c("exp4", "exp2", "null")),
  matrices = list(C, 2*G, G %*% Cinv %*% G),
  h = matern_obj$h
)

Non-stationary precision matrix operator with custom matrix combinations

Description

Creates a flexible non-stationary precision matrix K by allowing both linear combinations and matrix multiplications with parameter-dependent diagonal matrices. This enables modeling spatially varying parameters through basis expansions.

Usage

generic_ns(
  theta_K,
  matrices,
  position,
  h,
  model = "generic_ns",
  B_theta_K = NULL,
  trans = NULL,
  mesh = NULL,
  zero_trace = FALSE,
  param_name = NULL,
  param_trans = NULL,
  ...
)

Arguments

theta_K

List of parameter vectors (in real scale)
matrices

List of fixed matrices to be used in the model
position

List of vectors specifying matrix combinations (required)
h

The integration weights vector
model

The model type name stored on the returned operator
B_theta_K

List of basis matrices for parameters, defaults to matrices of 1s if not provided
trans

List of transformations for each parameter (one transformation per parameter)
mesh

The mesh object
zero_trace

Whether the trace of K should be zero
param_name

Optional parameter names stored on the returned operator for diagnostics
param_trans

Optional parameter transformations stored on the returned operator for diagnostics
...

Additional arguments (ignored)

Details

The key distinction from 'generic' is that 'generic_ns': 1. Requires explicit position specification for matrix combinations 2. Converts parameters to diagonal matrices for multiplication with fixed matrices 3. Allows for basis expansions via B_theta_K

The 'generic_ns' operator constructs K through the following steps:

1. Create diagonal matrices for each parameter using basis expansions: D_param = diag(B_param * theta_param) 2. Combine parameter matrices with fixed matrices according to position specifications 3. Sum all the resulting combinations

The 'position' parameter defines how matrices are combined: - Each element of the list represents a term in the sum - Each vector contains indices of matrices to multiply (parameter matrices first, then fixed matrices) - For example: 'position = list(c(1, 3), c(2, 4))' with one parameter means: Multiply the parameter diagonal matrix (index 1) with fixed matrix 1 (index 3), then add parameter diagonal matrix 2 (index 2) multiplied by fixed matrix 2 (index 4)

Spatially-varying parameters can be created using basis matrices in B_theta_K.

Value

An object of class ngme_operator. The object stores the sparse non-stationary precision matrix K, integration weights h, the flattened parameter vector theta_K, basis matrices B_theta_K, parameter grouping metadata param_map, matrix combination rules position, and an update function for rebuilding K. It is intended for use as the latent model inside f().

Examples

# Simple example with one parameter
n <- 5
A <- matrix(1, n, n)
B <- matrix(2, n, n)
alpha <- 0.4

# Create a simple generic_ns model: D_alpha * A + B
model <- generic_ns(
  theta_K = list(alpha = c(alpha)),
  matrices = list(A, B),
  position = list(c(1, 2), c(3)), # D_alpha * A + B
  h = rep(1, n)
)

# AR1 model representation: rho * C + G
ar1_obj <- ar1(1:10, rho = 0.5)
g <- name2fun("tanh", inv = TRUE)
generic_ar1 <- generic_ns(
  theta_K = list(x = c(g(0.5))), # Transform from real scale
  trans = list(x = "tanh"), # Apply tanh transformation
  matrices = list(ar1_obj$C, ar1_obj$G),
  position = list(c(1, 2), c(3)), # D_rho * C + G
  h = ar1_obj$h,
  mesh = 1:10
)

# Spatially varying Matern model
# Create a simple 1D mesh
mesh <- fmesher::fm_mesh_1d(seq(0, 1, length.out = 10))

# Create basis for spatially-varying kappa
n <- 10
B_kappa <- matrix(0, n, 2)
B_kappa[1:(n / 2), 1] <- 1 # First half of the domain
B_kappa[(n / 2 + 1):n, 2] <- 1 # Second half of the domain

# Create standard Matern components
matern_model <- matern(mesh)

# Create model with space-varying kappa
ns_model <- generic_ns(
  theta_K = list(kappa = c(log(1), log(2))), # Different values for each region
  matrices = list(matern_model$C, matern_model$G),
  B_theta_K = list(kappa = B_kappa),
  trans = list(kappa = "exp2"), # kappa^2 transformation for C
  h = matern_model$h,
  position = list(c(1, 2), c(3)), # D_kappa * C + G
  mesh = mesh
)

Extracts design matrix from a formula and data.

Description

This function takes a formula and a data frame, and returns a design matrix based on the right-hand side of the formula. If the formula is '~.', it treats all columns in the data as the design matrix. Otherwise, it constructs a design matrix without an intercept based on the specified formula.

Usage

get_data_from_formula(form, data)

Arguments

form

A formula object, e.g., '~ x1 + x2' or '~.'.
data

A data frame or matrix from which to extract the design matrix.

Value

A numeric matrix representing the design matrix.

Examples

data(mtcars)
# Extract a design matrix for 'mpg' and 'cyl'
X1 <- get_data_from_formula(~ mpg + cyl, mtcars)
# Extract a design matrix for all columns
X2 <- get_data_from_formula(~., mtcars)
# Extract a design matrix for 'hp'
X3 <- get_data_from_formula(~hp, mtcars)

Calculate parameter distance from true values

Description

Compute ||theta - theta_hat|| distance for all parameters across iterations. This provides a single measure of convergence combining all parameters.

Usage

get_parameter_distance(
  trajectories,
  true_values,
  norm_type = "euclidean",
  chain_summary = "mean"
)

Arguments

trajectories

ngme_trajectories object from get_trace_trajectories()
true_values

named list or vector of true parameter values
norm_type

character, type of norm to use: "euclidean" (default), "manhattan", "max"
chain_summary

character, how to summarize across chains: "mean" (default), "median", "chain1"

Value

numeric vector of distances across iterations

Get trace trajectories from ngme fitting

Description

Extract numerical trajectories of parameters during ngme optimization without creating plots. This function provides the underlying data used by traceplot().

Usage

get_trace_trajectories(ngme, name = "general", apply_transform = TRUE)

Arguments

ngme

ngme object
name

name of latent models, otherwise get fixed effects and measurement noise should be in names(ngme$models) or other
apply_transform

logical, whether to apply parameter transformations (default: TRUE)

Value

list containing:

trajectories

named list of trajectory matrices for each parameter

parameter_names

character vector of parameter names

transformations

list of transformation functions used

n_chains

number of chains

n_iterations

number of iterations

get the trajectories of parameters of the model

Description

Get the trajectories of the parameters of the model

Usage

get_trajectories(ngme_object, model_name)

Arguments

ngme_object

ngme object
model_name

name of the model, if NULL, return the trajectories of the parameters of the measurement noise and fixed effects

Value

a list of trajectories, each element is a matrix of trajectories (n_samples * n_trajectories)

The Generalised Inverse-Gaussian (GIG) Distribution

Description

Density, distribution function, quantile function and random generation for the generalised inverse-Gaussian distribution with parameters p, a and b.

Usage

dgig(x, p, a, b, log = FALSE)

rgig(n, p, a, b, seed = 0)

pgig(q, p, a, b, lower.tail = TRUE, log.p = FALSE)

qgig(prob, p, a, b, lower.tail = TRUE, log.p = FALSE)

Arguments

x, q

vector of quantiles.
p

parameter p.
a, b

parameters a and b. Must be positive.
log, log.p

logical; if TRUE, probabilities/densities pp are returned as log(p)log(p).
n

number of observations.
seed

Seed for the random generation.
lower.tail

logical; if TRUE, probabilities are P[X≤x]P[X\leq x], otherwise, P[X>x]P[X>x].
prob

vector of probabilities.

Details

The generalised inverse-Gaussian distribution has density given by

f(x;p,a,b)=((a/b)p/2)/(2Kp(ab))xp−1exp⁡{−(a/2)x−(b/2)/x},f(x; p, a, b) = ((a/b)^{p/2})/(2K_p(\sqrt{ab})) x^{p-1} \exp\{-(a/2)x - (b/2)/x\},

where KpK_p is modified Bessel function of the second kind of order pp, x>0x>0, a,b>0a,b>0 and p∈Rp\in\mathbb{R}. See Jørgensen (1982) for further details.

Value

dgig gives the density, pgig gives the distribution function, qgig gives the quantile function, and rgig generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rgig.

References

Jørgensen, Bent (1982). Statistical Properties of the Generalized Inverse Gaussian Distribution. Lecture Notes in Statistics. 9. New York–Berlin: Springer-Verlag. doi:10.1007/978-1-4612-5698-4

See Also

dnig, dig, digam

Examples

rgig(20, p = 1, a = 1, b = 1)
pgig(0.4, p = 1, a = 1, b = 1)
qgig(0.8, p = 1, a = 1, b = 1)
plot(function(x){dgig(x, p = 1, a = 1, b = 1)}, main =
"Generalised inverse-Gaussian density", ylab = "Probability density",
xlim = c(0,10))

The Inverse-Gaussian (IG) Distribution

Description

Density, distribution function, quantile function and random generation for the inverse-Gaussian distribution with parameters a and b.

Usage

dig(x, a, b, log = FALSE)

rig(n, a, b, seed = 0)

pig(q, a, b, lower.tail = TRUE, log.p = FALSE)

qig(p, a, b, lower.tail = TRUE, log.p = FALSE)

Arguments

x, q

vector of quantiles.
a, b

parameters a and b. Must be positive.
log, log.p

logical; if TRUE, probabilities/densities pp are returned as log(p)log(p).
n

number of observations.
seed

Seed for the random generation.
lower.tail

logical; if TRUE, probabilities are P[X≤x]P[X\leq x], otherwise, P[X>x]P[X>x].
p

vector of probabilities.

Details

The inverse-Gaussian distribution has density given by

f(x;a,b)=b2πx3exp⁡(−a2x−b2x+ab),f(x; a, b) = \frac{\sqrt{b}}{\sqrt{2\pi x^3}}\exp( -\frac{a}{2}x -\frac{b}{2x} + \sqrt{ab}),

where x>0x>0 and a,b>0a,b>0. In this parameterization, E(X)=b/aE(X) = \sqrt{b}/\sqrt{a}. See Tweedie (1957a, 1957b) for further details.

Value

dig gives the density, pig gives the distribution function, qig gives the quantile function, and rig generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rig.

References

Tweedie, M. C. K. (1957a). "Statistical Properties of Inverse Gaussian Distributions I". Annals of Mathematical Statistics. 28 (2): 362–377. doi:10.1214/aoms/1177706964

Tweedie, M. C. K. (1957b). "Statistical Properties of Inverse Gaussian Distributions II". Annals of Mathematical Statistics. 28 (3): 696–705. doi:10.1214/aoms/1177706881

See Also

dnig, dgig, digam

Examples

rig(100, a = 1, b = 1)
pig(0.4, a = 1, b = 1)
qig(0.8, a = 1, b = 1)
plot(function(x){dig(x, a = 1, b = 1)}, main =
"Inverse-Gaussian density", ylab = "Probability density",
xlim = c(0,10))

The Inverse-Gamma (IGam) Distribution

Description

Density, distribution function, quantile function and random generation for the inverse-Gamma distribution with parameters a and b.

Usage

digam(x, a, b, log = FALSE)

rigam(n, a, b)

pigam(q, a, b, lower.tail = TRUE, log.p = FALSE)

qigam(p, a, b, lower.tail = TRUE, log.p = FALSE)

Arguments

x, q

vector of quantiles.
a, b

parameters a and b. Must be positive.
log, log.p

logical; if TRUE, probabilities/densities pp are returned as log(p)log(p).
n

number of observations.
lower.tail

logical; if TRUE, probabilities are P[X≤x]P[X\leq x], otherwise, P[X>x]P[X>x].
p

vector of probabilities.

Details

The inverse-Gamma distribution has density given by

f(x;a,b)=baΓ(a)x−a−1exp⁡(−bx),f(x; a, b) = \frac{b^a}{\Gamma(a)}x^{-a-1}\exp( -\frac{b}{x}),

where x>0x>0 and a,b>0a,b>0.

Value

digam gives the density, pigam gives the distribution function, qigam gives the quantile function, and rigam generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rig.

See Also

dnig, dgig

Examples

rigam(100, a = 1, b = 1)
pigam(0.4, a = 1, b = 1)
qigam(0.8, a = 1, b = 1)
plot(function(x){digam(x, a = 1, b = 1)}, main =
"Inverse-Gamma density", ylab = "Probability density",
xlim = c(0,10))

ngme iid model specification

Description

ngme iid model specification

Usage

iid(mesh = NULL)

Arguments

mesh

integer or factor, index to build the mesh

Value

ngme_operator object

Create Time Series Cross-Validation Indices

Description

Creates indices for time series cross-validation with options for expanding window or sliding window approaches. Supports both single-step and multi-step forecasting, and can handle replicated observations.

Usage

make_time_series_cv_index(
  time_idx,
  train_length = NULL,
  test_length = 1,
  replicate = time_idx,
  gap = 0
)

Arguments

time_idx

A numeric vector of time indices in ascending order
train_length

An integer specifying the fixed length of training sets. If NULL (default), an expanding window approach is used.
test_length

An integer specifying the number of observations to include in each test set. Default is 1 (single-step forecasting).
replicate

An optional vector of the same length as time_idx, indicating which observations belong to the same replicate group. When provided, ensures that all observations with the same replicate value are either entirely in the training set or entirely in the test set.
gap

An integer specifying the gap between the training set and test set. Default is 0 (no gap, test set starts immediately after training set). For example, gap=1 means skip one time point between training and test (useful for 2-step ahead forecasting).

Details

Time series cross-validation requires respecting the temporal order of observations. This function implements two common approaches:

1. Expanding window (when train_length = NULL): The training set grows with each fold, starting with a minimal set and expanding to include all but the test data.

2. Sliding window (when train_length is specified): Uses a fixed-length window that slides through the time series, maintaining the same training size across folds.

The test_length parameter allows for multi-step forecasting evaluation.

When replicate is provided, the function ensures that all observations with the same replicate value are kept together, either all in the training set or all in the test set. This is useful for scenarios where multiple observations at the same time point should be treated as a group.

The gap parameter creates a separation between training and test sets, which is useful for multi-step ahead forecasting validation.

Value

A list with two components:

train

A list of numeric vectors, where each vector contains the time indices for training in that fold

test

A list of numeric vectors, where each vector contains the time indices for testing in that fold

Examples

# Expanding window approach with single-step forecasting
cv_expanding <- make_time_series_cv_index(1:10)

# Sliding window approach with window size 3 and single-step forecasting
cv_sliding <- make_time_series_cv_index(1:10, train_length = 3)

# Sliding window with multi-step forecasting (predict 2 steps ahead)
cv_multistep <- make_time_series_cv_index(1:10, train_length = 3, test_length = 2)

# Working with replicates
time_idx <- c(1, 1, 1, 2, 2, 3, 3)
replicates <- c(1, 1, 1, 2, 2, 3, 3)
cv_with_replicates <- make_time_series_cv_index(time_idx, replicate = replicates)

# 2-step ahead forecasting with a gap
cv_with_gap <- make_time_series_cv_index(1:10, gap = 1)

ngme Matern SPDE model specification

Description

ngme Matern SPDE model specification

Usage

matern(
  mesh = NULL,
  alpha = 2,
  fix_alpha = TRUE,
  kappa = NULL,
  theta_kappa = NULL,
  rational_order = 2,
  B_kappa = NULL
)

Arguments

mesh

an fmesher::fm_mesh_2d object, mesh for build the SPDE model
alpha

SPDE smoothness parameter (alpha = 2beta) 2 or 4 for integer case, otherwise for fractional case
fix_alpha

if FALSE, enable fractional modeling
kappa

the range parameter, for the stationary model, it is the range parameter
theta_kappa

kappa = exp(B_kappa * theta_kappa), for the non-stationary model, it is the range parameter
rational_order

order of rational approximation for fractional operators (default: 2)
B_kappa

bases for theta_kappa, ignore if use the stationary model

Value

ngme_operator object

taking mean over a list of nested lists

Description

taking mean over a list of nested lists

Usage

mean_list(lls, weights = NULL)

Arguments

lls

a list
weights

weights of each list

Value

a list of nested lists

Examples

ls <- list(
  list(a = 1, b = 2, t = "nig", ll = list(a = 1, b = 2, w = "ab")),
  list(a = 3, b = 5, t = "nig", ll = list(a = 1, b = 6, w = "ab")),
  list(a = 5, b = 5, t = "nig", ll = list(a = 4, b = 2, w = "ab"))
)
mean_list(ls)

Merge 2 noise into 1 noise

Description

Merge 2 noise into 1 noise

Usage

merge_noise(noise1, noise2)

Arguments

noise1

noise 1
noise2

noise 2

Value

merged noise

Merge model of replicates into model of 1 replicate given train_idx and test_idx, the merged model contains all the information of train_idx from different replicates.

Description

Merge model of replicates into model of 1 replicate given train_idx and test_idx, the merged model contains all the information of train_idx from different replicates.

Usage

merge_replicates(ngme, train_idx, test_idx)

Arguments

ngme

a ngme object
train_idx

a vector of indices of train data
test_idx

a vector of indices of test data

Momentum SGD optimization

Description

Momentum SGD optimization

Usage

momentum(stepsize = 0.05, beta1 = 0.9, beta2 = 1 - beta1)

Arguments

stepsize

stepsize for SGD
beta1

beta1 for momentum
beta2

beta2 for momentum

Details

The update rule for momentum is:

vt=β1vt−1+β2gtv_t = \beta_1 v_{t-1} + \beta_2 g_t

xt+1=xt−stepsize∗vtx_{t+1} = x_t - \text{stepsize} * v_t

Value

a list of control variables for optimization (used in control_opt function)

Convert transformation name to function

Description

This function converts a transformation name to the corresponding function. Convertion is from original scale to real scale by default. Available transformations:

  • exp4: exp(4x)exp(4x), inverse is log(x)/4log(x)/4

  • exp2: exp(2x)exp(2x), inverse is log(x)/2log(x)/2

  • tanh: Hyperbolic tangent transformation used for AR1 parameter, uses ar1_th2a and ar1_a2th

  • identity: Identity function, no transformation

  • exp: exp(x)exp(x), inverse is log(x)log(x)

  • sqrt: sqrt(x)sqrt(x), inverse is x2x^2

  • square: x2x^2, inverse is sqrt(x)sqrt(x)

  • log: log(x)log(x), inverse is exp(x)exp(x)

Usage

name2fun(trans, inv = FALSE)

Arguments

trans

Character string specifying the transformation type from original scale to real scale
inv

Logical, if TRUE returns the transformation function from real scale to original scale

Value

A function that applies the specified transformation

Fit an additive linear mixed effect model over replicates

Description

ngme function performs an analysis of non-gaussian additive models. It does the maximum likelihood estimation via stochastic gradient descent. The prediction of unknown location can be performed by leaving the response variable to be NA. The likelihood is specified by family. The model estimation control can be setted in control using control_opt() function, see ?control_opt for details. See ngme_model_types() for available models.

Usage

ngme(
  formula,
  data,
  family = "normal",
  control_opt = NULL,
  control_ngme = NULL,
  group = NULL,
  replicate = NULL,
  start = NULL,
  moving_window = 1,
  prior_beta = NULL,
  debug = FALSE
)

Arguments

formula

formula
data

a dataframe or a list providing data (Only response variable can contain NA value, NA value in other columns will cause problem)
family

likelihood type, same as measurement noise specification. It can be provided as:

  1. a string, e.g., "normal", "nig", "t".

  2. an ngme noise object, e.g., noise_normal(), noise_nig(mu = 0, sigma = 1, nu = 1), noise_t(nu=5).

control_opt

control for optimizer. by default it is control_opt(). See ?control_opt for details.
control_ngme

control for ngme model. by default it is control_ngme(). See ?control_ngme for details.
group

factor, used for bivariate model, indicating which group the observation belongs to
replicate

factor, used for divide data into different replicates
start

starting ngme object (usually object from last fit)
moving_window

number of iterations to average the estimation
prior_beta

prior specification for fixed effects ('beta'), created by prior_*() or priors(...).
debug

toggle debug mode

Value

random effects (for different replicate) + models(fixed effects, measuremnt noise, and latent process)

Examples

ngme(
  formula = Y ~ x1 + f(
    x2,
    model = ar1(rho = 0.5),
    noise = noise_nig()
  ) + f(x1,
    model = rw1(),
    noise = noise_normal(),
  ),
  family = noise_normal(sd = 0.5),
  data = data.frame(Y = 1:5, x1 = 2:6, x2 = 3:7),
  control_opt = control_opt(
    estimation = FALSE
  )
)

Convert sparse matrix into sparse dgCMatrix

Description

Convert sparse matrix into sparse dgCMatrix

Usage

ngme_as_sparse(G)

Arguments

G

matrix

Value

sparse dgCMatrix

Batch-Means Confidence Intervals from an ngme Fit

Description

Compute Xi-style batch-means confidence intervals directly from trajectories stored in an 'ngme' object.

Usage

ngme_batch_ci(
  ngme,
  name = "all",
  level = 0.95,
  alpha = 0.501,
  M = NULL,
  N = NULL,
  apply_transform = TRUE,
  drop_burnin = TRUE,
  burnin_iter = 0
)

Arguments

ngme

fitted 'ngme' object with 'store_traj = TRUE' during optimization.
name

either '"all"' (default) to use all parameters jointly, a latent model name, or '"general"' for measurement-noise/fixed-effect block.
level

confidence level.
alpha

stepsize decay exponent in (1/2,1)(1/2, 1).
M

number of retained batches (excluding burn-in batch 0).
N

decorrelation constant in the batch boundary formula.
apply_transform

logical; if 'TRUE', inference is first computed in the raw optimization space and then mapped to user scale via post-hoc Delta method (Post-delta).
drop_burnin

logical; if 'TRUE', discard batch 0.
burnin_iter

non-negative integer. Explicitly discard the first 'burnin_iter' optimization iterations before Xi-style batching.

Value

Same structure as [batch_means_ci()] with extra metadata fields 'name' and 'apply_transform'.

variance of the data or the latent field

Description

Compute the variance of the data or the latent field

Usage

ngme_cov_matrix(ngme_object, model_name = "data", replicate = 1)

Arguments

ngme_object

ngme_model
model_name

if model_name = "data", then return the covariance matrix of the data (without measurement noise) if the model_name is the name or index of the latent, then return the covariance matrix of the latent field
replicate

which replicate (default = 1)

Value

a data.frame of posterior samples (mesh_size * n_post_samples)

ngme make mesh for different replicates

Description

Make different mesh for different replicates

Usage

ngme_make_mesh_repls(data, map, replicate, mesh_type = "regular")

Arguments

data

provide the data.frame
map

provide the map to make mesh, i.g. ~x+y, x and y will be extracted from data to make a 2d mesh
replicate

provide the replicate information, i.g. ~id
mesh_type

type of mesh, "regular" means use all the point from same replicate to make a mesh

Value

a list of mesh of length of different replicates

Show ngme model types

Description

Show ngme model types

Usage

ngme_model_types()

Value

available types for models

ngme noise specification

Description

Function for specifying ngme noise. Please use noise_nig and noise_normal for simpler usage. Use ngme_noise_types() to check all the available types.

Usage

ngme_noise(
  noise_type,
  mu = 0,
  sigma = 1,
  nu = 1,
  B_mu = NULL,
  theta_mu = NULL,
  B_sigma = NULL,
  theta_sigma = NULL,
  B_nu = NULL,
  theta_nu = NULL,
  theta_sigma_normal = NULL,
  B_sigma_normal = NULL,
  fix_theta_mu = FALSE,
  fix_theta_sigma = FALSE,
  fix_rho = FALSE,
  fix_theta_sigma_normal = FALSE,
  fix_theta_nu = FALSE,
  V = NULL,
  fix_V = FALSE,
  single_V = FALSE,
  share_V = FALSE,
  corr_measurement = FALSE,
  index_corr = NULL,
  map_corr = NULL,
  nu_lower_bound = 0,
  rho = double(0),
  prior = NULL,
  ...
)

noise_normal(
  sigma = NULL,
  theta_sigma = NULL,
  B_sigma = matrix(1),
  corr_measurement = FALSE,
  index_corr = NULL,
  ...
)

noise_nig(
  mu = NULL,
  sigma = NULL,
  nu = NULL,
  V = NULL,
  theta_mu = NULL,
  theta_sigma = NULL,
  theta_nu = NULL,
  nu_lower_bound = 0,
  B_mu = matrix(1),
  B_sigma = matrix(1),
  B_nu = matrix(1),
  corr_measurement = FALSE,
  index_corr = NULL,
  ...
)

noise_gal(
  mu = NULL,
  sigma = NULL,
  nu = NULL,
  V = NULL,
  theta_mu = NULL,
  theta_sigma = NULL,
  theta_nu = NULL,
  nu_lower_bound = 0.01,
  B_mu = matrix(1),
  B_sigma = matrix(1),
  B_nu = matrix(1),
  corr_measurement = FALSE,
  index_corr = NULL,
  ...
)

noise_skew_t(
  mu = NULL,
  sigma = NULL,
  nu = NULL,
  theta_mu = NULL,
  theta_sigma = NULL,
  theta_nu = NULL,
  nu_lower_bound = 0.01,
  B_mu = matrix(1),
  B_sigma = matrix(1),
  B_nu = matrix(1),
  corr_measurement = FALSE,
  index_corr = NULL,
  ...
)

noise_t(
  nu = NULL,
  theta_nu = NULL,
  nu_lower_bound = 0,
  B_nu = matrix(1),
  corr_measurement = FALSE,
  index_corr = NULL,
  ...
)

noise_normal_nig(
  sigma_normal = NULL,
  mu = NULL,
  sigma_nig = NULL,
  nu = NULL,
  V = NULL,
  theta_mu = NULL,
  theta_sigma_nig = NULL,
  theta_sigma_normal = NULL,
  theta_nu = NULL,
  B_mu = matrix(1),
  B_sigma_nig = matrix(1),
  B_sigma_normal = matrix(1),
  B_nu = matrix(1),
  corr_measurement = FALSE,
  index_corr = NULL,
  ...
)

Arguments

noise_type

type of noise, "normal", "nig", "gal", "t", "skew_t"
mu

specify the NIG noise parameter mu, see ?nig
sigma

specify the noise parameter sigma, see ?nig
nu

specify the noise parameter nu, see ?nig and ?gal
B_mu

Basis matrix for mu (if non-stationary)
theta_mu

specify a non-stationary noise using theta_mu
B_sigma

Basis matrix for sigma (if non-stationary)
theta_sigma

specify a non-stationary noise using theta_sigma
B_nu

Basis matrix for nu (if non-stationary)
theta_nu

specify a non-stationary noise using theta_nu
theta_sigma_normal

for normal nosie with nig noise sharing same parameter
B_sigma_normal

for normal nosie with nig noise sharing same parameter
fix_theta_mu

fix the parameter of theta_mu
fix_theta_sigma

fix the parameter of theta_sigma, can be a single logical value or a vector of logical values with length equal to length(theta_sigma)
fix_rho

fix the parameter of rho
fix_theta_sigma_normal

fix the parameter of sigma_normal, used in noise_normal_nig()
fix_theta_nu

fix the parameter of nu
V

start value for V
fix_V

fix the sampling of V gives Y|W ~ N(mean * sigma, sigma^2)
single_V

TRUE if V is a single number
share_V

used only for bivariate model
corr_measurement

TRUE if we use correlated measurement noise
index_corr

used when corr_measurement=TRUE, indicate which observation has correlation
map_corr

1d, 2d, or formula, used when corr_measurement=TRUE, specify use which covariate to infer the index_corr.
nu_lower_bound

specify the lower bound of parameter nu; effective parametrization is nu=nulowerbound+exp⁡(Bνθν)nu = nu_lower_bound + \exp(B_\nu \theta_\nu) so theta_nu remains unconstrained. Default 0.
rho

used when corr_measurement=TRUE, starting point for correlation
prior

prior specification created by priors(...). Supported keys are mu, sigma, and nu.
...

additional arguments
sigma_normal

for normal nosie with nig noise sharing same parameter
sigma_nig

similar to sigma_normal
theta_sigma_nig

similar to theta_sigma_normal
B_sigma_nig

similar to B_sigma_nig

Details

The parameterization is given in ?nig and ?gal. Moreover, for specifying non-stationary mu and sigma, nu

μ=Bμθμ,\mu = B_{\mu} \theta_{\mu},

and

σ=exp⁡(Bσθσ),\sigma = \exp (B_{\sigma} \theta_{\sigma}),

ν=νlower+exp⁡(Bνθν).\nu = \nu_{\mathrm{lower}} + \exp (B_{\nu} \theta_{\nu}).

Value

a list of specification of noise

Examples

noise_normal(sigma = 2)
noise_nig(mu = 1, sigma = 2, nu = 1)
noise_gal(mu = 1, sigma = 2, nu = 1)
noise_skew_t(mu = 0, sigma = 1, nu = 5)
noise_t(nu = 5)

Show ngme noise types

Description

Show ngme noise types

Usage

ngme_noise_types()

Value

available types for noise

List supported optimizers

Description

This function returns a list of supported optimizers in the ngme package. The optimizers are categorized into three groups:

  • Gradient descent: "sgd", "sgld", "momentum", "adaptive_gd"

  • Adaptive learning rate: "adagrad", "rmsprop", "adam", "adamW"

  • Preconditioner: "precond_sgd", "bfgs"

Usage

ngme_optimizers()

Value

a character vector of supported optimizers

Parse the formula for ngme function

Description

Parse the formula for ngme function

Usage

ngme_parse_formula(
  fm,
  data,
  control_ngme,
  noise,
  group,
  replicate,
  prior_beta,
  standardize
)

Arguments

fm

formula
data

data.frame
control_ngme

control_ngme
noise

noise
group

group factor
replicate

replicate vector
prior_beta

prior specification for fixed effects
standardize

logical, whether fixed effects are standardized

Value

a list (replicate) of ngme_replicate models

posterior samples of different latent models

Description

Extract the posterior samples of different latent models

Usage

ngme_post_samples(ngme_object, model_name = 1, type = "W", replicate = 1)

Arguments

ngme_object

ngme object
model_name

name of the model, or index of the model
type

type of samples, "W" or "V"
replicate

which replicate

Value

a data.frame of posterior samples (mesh_size * n_post_samples)

Show ngme priors

Description

Show ngme priors

Usage

ngme_prior_types()

Value

available types of priors

Access the result of a ngme fitted model

Description

Access the result of a ngme fitted model

Usage

ngme_result(ngme_object, model = NULL, transformed = TRUE)

Arguments

ngme_object

a ngme fitted model object
model

latent model name to filter by (e.g., "my_ar1", "my_matern", "my_rw1", etc.), "data" for fixed effects and measurement noise, if NULL, return all models
transformed

logical, if TRUE (default) return transformed parameters, if FALSE return raw parameters

Value

a list of parameters for the specified model, or all models if model is NULL

Examples

# Fit a simple AR(1) model
set.seed(1)
Y <- 1:10
n_obs <- length(Y)
x1 <- rnorm(n_obs)
x2 <- runif(n_obs)

ngme_out <- ngme(
  Y ~ x1 + x2 + f(
    1:n_obs,
    name = "my_ar",
    model = ar1(rho = 0.5),
    noise = noise_nig(mu = 2, sigma = 3, nu = 1)
  ),
  data = data.frame(x1 = x1, x2 = x2, Y = Y),
  control_opt = control_opt(estimation = FALSE)
)

# Get all model parameters (transformed)
all_params <- ngme_result(ngme_out)
# Returns: list(my_ar = list(rho = 0.5, mu = 2, sigma = 3, nu = 1),
#               data = list(fixed_effects = c(...), sigma = 0.5))

# Get parameters for specific latent model
ar_params <- ngme_result(ngme_out, model = "my_ar")
# Returns: list(rho = 0.5, mu = 2, sigma = 3, nu = 1)

# Get raw (untransformed) parameters
ar_raw <- ngme_result(ngme_out, model = "my_ar", transformed = FALSE)
# Returns: list(theta_rho = 1.099, mu = 2, sigma = 3, nu = 1)

# Get fixed effects and measurement noise
data_params <- ngme_result(ngme_out, model = "data")
# Returns: list(fixed_effects = c(...), sigma = 0.5, ...)

# For models with multiple latent processes
ngme_out2 <- ngme(
  Y ~ x1 + x2 + f(
    1:n_obs,
    name = "my_ar",
    model = ar1(rho = 0.5),
    noise = noise_nig(mu = 2, sigma = 3, nu = 1)
  ) + f(
    1:n_obs,
    name = "my_ou",
    model = ou(),
    noise = noise_normal(sigma = 1)
  ),
  data = data.frame(x1 = x1, x2 = x2, Y = Y),
  control_opt = control_opt(estimation = FALSE)
)

# Get all models
all_models <- ngme_result(ngme_out2)
# Returns: list(my_ar = list(...), my_ou = list(...), data = list(...))

# Get specific model
ou_params <- ngme_result(ngme_out2, model = "my_ou")
# Returns: list(theta_K1 = 0.5, sigma = 1)

Quantile Confidence Intervals from SGLD Samples

Description

Compute parameter-wise quantile confidence intervals from posterior-like SGLD samples returned by [ngme_sgld_samples()].

Usage

ngme_sgld_ci(samples, lower = 0.025, upper = 0.975)

Arguments

samples

data.frame (or list of data.frames) returned by [ngme_sgld_samples()].
lower

lower quantile probability (e.g. 0.025).
upper

upper quantile probability (e.g. 0.975).

Value

A list with posterior-like point estimates ('estimates'), quantile intervals ('ci'), and sample covariance matrix ('covariance'). Class: 'ngme_sgld_ci'.

Extract Posterior-like Samples from Stored SGLD Trajectories

Description

Build posterior-like samples from optimizer trajectories by dropping an initial burn-in segment and applying thinning.

Usage

ngme_sgld_samples(
  ngme,
  name = "all",
  burnin_iter = 0,
  thinning = 1,
  apply_transform = TRUE,
  combine_chains = TRUE
)

Arguments

ngme

fitted 'ngme' object with 'store_traj = TRUE'.
name

parameter block to extract: '"all"' (default), latent model name, or '"general"'.
burnin_iter

non-negative integer. Number of initial iterations to discard before sampling.
thinning

positive integer thinning interval.
apply_transform

logical; apply parameter transforms to user scale.
combine_chains

logical; if 'TRUE', return one combined data.frame, otherwise return one data.frame per chain.

Value

A data.frame (or list of data.frames when 'combine_chains = FALSE') with columns '.chain', '.draw', '.iter', and one column per parameter.

Make observation matrix for time series

Description

Make observation matrix for time series

Usage

ngme_ts_make_A(loc, replicate = NULL, range = c(min(loc), max(loc)))

Arguments

loc

integers (after sorting, no gaps > 1)
replicate

indicating replicate measure at same location
range

range for the mesh by default range=(min(loc), max(loc))

Value

A matrix (length(loc) * length(unique(loc)))

Examples

ngme_ts_make_A(c(1, 2, 2), replicate = c(1, 1, 2))
ngme_ts_make_A(c(1, 2, 2), range = c(1, 5))

Update ngme2 to the latest stable version

Description

Update ngme2 to the latest stable version

Usage

ngme_update()

Value

void

The Normal Inverse-Gaussian (NIG) Distribution

Description

Density, distribution function, quantile function and random generation for the normal inverse-Gaussian distribution with parameters p, a and b.

Usage

dnig(x, delta, mu, nu, sigma, h = NULL, log = FALSE)

rnig(n, delta, mu, nu, sigma, h = NULL, seed = 0)

pnig(q, delta, mu, nu, sigma, h = NULL, lower.tail = TRUE, log.p = FALSE)

qnig(p, delta, mu, nu, sigma, h = NULL, lower.tail = TRUE, log.p = FALSE)

Arguments

x, q

vector of quantiles.
delta

A numeric value for the location parameter.
mu

A numeric value for the shift parameter.
nu

A numeric value for the shape parameter.
sigma

A numeric value for the scaling parameter.
h

A numeric value for the additional parameter, see details.
log, log.p

logical; if TRUE, probabilities/densities pp are returned as log(p)log(p).
n

number of observations.
seed

Seed for the random generation.
lower.tail

logical; if TRUE, probabilities are P[X≤x]P[X\leq x], otherwise, P[X>x]P[X>x].
p

vector of probabilities.

Details

The normal inverse-Gaussian distribution has density given by

f(x;δ,μ,σ,ν)=eν+μ(x−δ)/σ2νμ2/σ2+ν2πνσ2+(x−δ)2K1((νσ2+(x−δ)2)(μ2/σ4+ν/σ2)),f(x; \delta, \mu, \sigma, \nu) = \frac{e^{\nu+\mu(x-\delta)/\sigma^2}\sqrt{\nu\mu^2/\sigma^2+\nu^2}}{\pi\sqrt{\nu\sigma^2+(x-\delta)^2}} K_1(\sqrt{(\nu\sigma^2+(x-\delta)^2)(\mu^2/\sigma^4+\nu/\sigma^2)}),

where KpK_p is modified Bessel function of the second kind of order pp, x>0x>0, ν>0\nu>0 and μ,δ,σ∈R\mu,\delta, \sigma\in\mathbb{R}. See Barndorff-Nielsen (1977, 1978 and 1997) for further details.

The additional parameter h is used when

V∼IG(ν,νh2)V\sim IG(\nu,\nu h^{2})

. By the infinite divisibility,

1hV∼IG(νh,νh)\frac{1}{h} V \sim IG(\nu h, \nu h)

. Then

δ+μV+σVZ\delta+\mu V + \sigma \sqrt{V} Z

has the distribution of

NIG(δ=−μh,μ=μh,σ=σh,ν=νh).NIG(\delta=-\mu h,\mu= \mu h, \sigma=\sigma \sqrt{h}, \nu=\nu h).

Value

dnig gives the density, pnig gives the distribution function, qnig gives the quantile function, and rnig generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rnig.

References

Barndorff-Nielsen, O. (1977) Exponentially decreasing distributions for the logarithm of particle size. Proceedings of the Royal Society of London.

Series A, Mathematical and Physical Sciences. The Royal Society. 353, 401–409. doi:10.1098/rspa.1977.0041

Barndorff-Nielsen, O. (1978) Hyperbolic Distributions and Distributions on Hyperbolae, Scandinavian Journal of Statistics. 5, 151–157.

Barndorff-Nielsen, O. (1997) Normal Inverse Gaussian Distributions and Stochastic Volatility Modelling, Scandinavian Journal of Statistics. 24, 1-13. doi:10.1111/1467-9469.00045

See Also

dgig, dig, digam

Examples

rnig(100, delta = 0, mu = 5, sigma = 1, nu = 1)
pnig(0.4, delta = 0, mu = 5, sigma = 1, nu = 1)
qnig(0.8, delta = 0, mu = 5, sigma = 1, nu = 1)
plot(function(x){dnig(x, delta = 0, mu = 5, sigma = 1, nu = 1)}, main =
"Normal inverse-Gaussian density", ylab = "Probability density",
xlim = c(0,10))

Moments of stationary NIG and GAL noise

Description

Compute the mean, variance, skewness, kurtosis, and excess kurtosis for stationary NIG and GAL noise. The distribution is evaluated through the normal mean-variance mixture

X=δ+μV+σVZ,X = \delta + \mu V + \sigma \sqrt{V} Z,

where Z∼N(0,1)Z \sim N(0, 1). For the package's stationary noise parameterization, delta = -mu, so the mean is zero.

Usage

noise_moments(noise, ...)

## S3 method for class 'ngme_noise'
noise_moments(noise, ...)

noise_nig_moments(delta = -mu, mu = 0, sigma = 1, nu = 1)

noise_gal_moments(delta = -mu, mu = 0, sigma = 1, nu = 1)

Arguments

noise

An ngme_noise object created by noise_nig() or noise_gal().
...

Currently unused.
delta

Location parameter. Defaults to -mu, matching the stationary ngme noise convention.
mu

Shift/skewness parameter in the mixture mean.
sigma

Positive Gaussian scale parameter.
nu

Positive mixing-distribution shape parameter.

Value

A named numeric vector with entries skewness, kurtosis, excess_kurtosis, variance, sd, central_moment3, and central_moment4.

Examples

noise_nig_moments(mu = 1, sigma = 2, nu = 3)
noise_gal_moments(mu = 1, sigma = 2, nu = 3)
noise_moments(noise_nig(mu = 1, sigma = 2, nu = 3))

Test OpenMP availability and report the number of threads.

Description

This function checks if OpenMP is available in the current R environment and, if so, reports the number of OpenMP threads detected. If OpenMP is not available, it reports a corresponding message. It relies on an internal or external function 'get_openmp_threads()' which is assumed to return the number of threads or 0 if unavailable.

Usage

openmp_test()

Value

Invisibly returns the detected number of OpenMP threads. A value of zero means OpenMP is not available.

Ornstein-Uhlenbeck Process Model

Description

Implements the exact discrete-time representation of the continuous Ornstein-Uhlenbeck (OU) process using a matrix formulation KX=δ\mathbf{K}\mathbf{X} = \boldsymbol{\delta}.

Usage

ou(mesh = NULL, theta = 1)

Arguments

mesh

numerical vector or 'inla.mesh.1d' object giving the ordered time locations. Must be strictly increasing.
theta

positive mean-reversion rate (stiffness parameter). Internally stored on the log scale so optimization remains unconstrained.

Details

The OU process is defined by the stochastic differential equation (SDE):

dXt=−θXtdt+σdWtdX_t = -\theta X_t dt + \sigma dW_t

where θ>0\theta > 0 is the mean-reversion rate (stiffness) and σ\sigma is the volatility (diffusion coefficient).

## Exact Discretization

Unlike Euler-Maruyama approximations, the OU process has an exact solution between time points tit_i and ti+1t_{i+1} with step size Δti\Delta t_i:

Xi+1=Xie−θΔti+ϵiX_{i+1} = X_i e^{-\theta \Delta t_i} + \epsilon_i

This is an AR(1) form with autoregressive coefficient:

ρi=e−θΔti\rho_i = e^{-\theta \Delta t_i}

and Gaussian noise:

ϵi∼N(0,σ22θ(1−ρi2))\epsilon_i \sim \mathcal{N}\left(0, \frac{\sigma^2}{2\theta}(1 - \rho_i^2)\right)

## Matrix Representation

The precision matrix K\mathbf{K} for the linear system KX=δ\mathbf{K}\mathbf{X} = \boldsymbol{\delta} is constructed as:

K=[1−ρ1200⋯0−ρ110⋯00−ρ21⋯0⋮⋮⋱⋱⋮00⋯−ρn−11]\mathbf{K} = \begin{bmatrix} \sqrt{1-\rho_1^2} & 0 & 0 & \cdots & 0 \\ -\rho_1 & 1 & 0 & \cdots & 0 \\ 0 & -\rho_2 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & -\rho_{n-1} & 1 \end{bmatrix}

where ρi=exp⁡(−θΔti)\rho_i = \exp(-\theta \Delta t_i) with Δti=ti+1−ti\Delta t_i = t_{i+1} - t_i.

**Why the 1−ρ12\sqrt{1-\rho_1^2} term?** The first row ensures that X1X_1 is drawn from the stationary distribution N(0,σ2/(2θ))\mathcal{N}(0, \sigma^2/(2\theta)). This scaling factor matches the variance of the first component of δ\boldsymbol{\delta} with the transition noise variance in subsequent steps.

## Non-uniform Mesh

For non-uniform meshes, each ρi\rho_i is computed locally based on the varying time step Δti\Delta t_i, allowing the model to naturally handle irregularly spaced observations.

## The hh vector

The hh vector represents integration weights for the mass matrix. For the OU model, it is set as:

h=[Δt1,Δt2,…,Δtn−1,Δtn−1]h = [\Delta t_1, \Delta t_2, \ldots, \Delta t_{n-1}, \Delta t_{n-1}]

where the last element is duplicated. This ensures proper weighting in the finite element formulation.

## The δ\boldsymbol{\delta} vector

In the equation KX=δ\mathbf{K}\mathbf{X} = \boldsymbol{\delta}, the vector δ\boldsymbol{\delta} contains independent identically distributed (i.i.d.) random variables with unit variance. The precision matrix K\mathbf{K} is constructed such that the resulting process X\mathbf{X} has the correct OU covariance structure.

Value

An 'ngme_operator' object containing:

  • 'K': the precision matrix

  • 'h': integration weights

  • 'theta_K': log-transformed theta for unconstrained optimization

  • 'mesh': the mesh object

References

Uhlenbeck, G. E., & Ornstein, L. S. (1930). On the theory of the Brownian motion. Physical review, 36(5), 823.

Examples

# Uniform mesh
mesh_uniform <- seq(0, 10, by = 0.5)
ou_uniform <- ou(mesh = mesh_uniform, theta = 0.5)
print(ou_uniform)

# Non-uniform mesh
mesh_nonuniform <- c(0, 1, 2.5, 3.5, 5, 8, 11)
ou_nonunif <- ou(mesh = mesh_nonuniform, theta = 0.8)
print(ou_nonunif)

Plot the density of one or more stationary noise objects

Description

This function plots the probability density function for one or more stationary noise objects (e.g., NIG, GAL, or normal noise). Multiple noise objects can be compared on the same plot.

Usage

## S3 method for class 'ngme_noise'
plot(x = NULL, ...)

Arguments

x

An ngme_noise object (required).
...

Additional ngme_noise objects to plot, or plotting parameters such as xlim. Named arguments will be used as legend labels.

Value

A ggplot object showing the density curves for the provided noise objects.

Examples

plot(noise_nig(mu = 1, sigma = 2, nu = 1))
plot(n1 = noise_nig(mu = 0, sigma = 1, nu = 1), n2 = noise_nig(mu = 1, sigma = 1.5, nu = 0.5))

Plot method for parameter_distance

Description

Plot method for parameter_distance

Usage

## S3 method for class 'parameter_distance'
plot(x, ...)

Arguments

x

parameter_distance object
...

additional arguments for ggplot

Value

A ggplot object showing the parameter-distance trajectory over optimization iterations. The y-axis is the stored distance ∣∣theta−theta^∣∣||theta - \hat{theta}||; plot annotations report the norm and across-chain summary used when x was computed.

Polynomial schedule helper

Description

Polynomial schedule helper

Usage

poly_decay(alpha = 0.501, t0 = 1, burnin_iter = 0)

Arguments

alpha

polynomial exponent in (1/2,1)(1/2,1).
t0

non-negative schedule offset.
burnin_iter

non-negative integer. Initial iterations without polynomial schedule scaling.

Details

Convenience helper that enables polynomial schedule and disables checkpoint-based decay.

Value

a stepsize_control() object.

Plot Posterior Distributions from SGLD Samples

Description

Visualize marginal posterior distributions of parameters extracted from [ngme_sgld_samples()] or summarized by [ngme_sgld_ci()].

Usage

posterior_plot(
  x,
  parameters = NULL,
  type = c("density", "histogram"),
  by_chain = FALSE,
  show_estimate = TRUE,
  show_interval = TRUE,
  lower = NULL,
  upper = NULL,
  bins = 30,
  ncol = NULL,
  ...
)

## S3 method for class 'ngme_sgld_ci'
plot(x, ...)

Arguments

x

a data.frame (or list of data.frames) returned by [ngme_sgld_samples()], or an object of class 'ngme_sgld_ci'.
parameters

optional character vector of parameter names to plot. Defaults to all parameter columns.
type

character; '"density"' (default) or '"histogram"'.
by_chain

logical; if 'TRUE', color distributions by chain when a '.chain' column is available.
show_estimate

logical; add a vertical line at the posterior mean.
show_interval

logical; add vertical lines for the equal-tail interval.
lower

lower quantile probability used for intervals. If 'x' is an 'ngme_sgld_ci' object, the stored value is used by default.
upper

upper quantile probability used for intervals. If 'x' is an 'ngme_sgld_ci' object, the stored value is used by default.
bins

number of bins for histograms.
ncol

number of facet columns. If 'NULL', use 2 columns for multiple parameters and 1 for a single parameter.
...

unused.

Value

A faceted ggplot object.

Compute the precision matrix for multivariate model

Description

Compute the precision matrix for multivariate model

Usage

precision_matrix_multivariate(
  p,
  operator_list,
  rho,
  theta = NULL,
  Q = NULL,
  scale = NULL
)

Arguments

p

dimension, should be integer and greater than 1
operator_list

a list of ngme_operator object (length should be p)
rho

vector with the p(p-1)/2 correlation parameters rho_11, rho_21, rho_22, ... rho_p1, rho_p2, ... rho_p(p-1)
theta

parameter for Q matrix (length of 1 when p=2, length of 3 when p=3)
Q

orthogonal matrix of dim p*p (provide when p > 3)
scale

A vector of length p with constants to multiply each operator matrix with

Details

The general model is defined as $D diag(L_1, ..., L_p) x = M$. D is the dependence matrix, it is paramterized by $D = Q(theta) * D_l(cor_mat)$, where $Q$ is the orthogonal matrix, and $D_l$ is matrix controls the cross-correlation. See the section 2.2 of Bolin and Wallin (2020) for exact parameterization of Dependence matrix.

Value

the precision matrix of the multivariate model

References

Bolin, D. and Wallin, J. (2020), Multivariate type G Matérn stochastic partial differential equation random fields. J. R. Stat. Soc. B, 82: 215-239. https://doi.org/10.1111/rssb.12351

Examples

rho <- c(-0.5, 0.5, -0.25) # correlation parameters
operator_list <- list(ar1(1:5, rho = 0.4), ar1(1:5, rho = 0.5), ar1(1:5, rho = 0.6))
precision_matrix_multivariate(3, operator_list, rho, theta = c(1, 2, 3))

Compute the precision matrix for multivariate spde Matern model

Description

Compute the precision matrix for multivariate spde Matern model

Usage

precision_matrix_multivariate_spde(
  p,
  mesh,
  rho,
  alpha_list = NULL,
  theta_K_list = NULL,
  variance_list = NULL,
  B_K_list = NULL,
  theta = NULL,
  Q = NULL
)

Arguments

p

dimension, should be integer and greater than 1
mesh

an fmesher::fm_mesh_2d object, mesh for build the SPDE model
rho

vector with the p(p-1)/2 correlation parameters rho_11, rho_21, rho_22, ... rho_p1, rho_p2, ... rho_p(p-1)
alpha_list

a list of SPDE smoothness parameter
theta_K_list

a list (length is p) of theta_K
variance_list

If provided, it should be a vector of length p, where the kth element corresponds to a desired variance of the kth field. The kth operator is then scaled by a constant c so that this variance is achieved in the stationary case (default no scaling)
B_K_list

a list (length is p) of B_K (non-stationary case)
theta

parameter for Q matrix (length of 1 when p=2, length of 3 when p=3)
Q

orthogonal matrix of dim p*p (provide when p > 3)

Details

The general model is defined as $D diag(L_1, ..., L_p) x = M$. D is the dependence matrix, it is paramterized by $D = Q(theta) * D_l(cor_mat)$, where $Q$ is the orthogonal matrix, and $D_l$ is matrix controls the cross-correlation. See the section 2.2 of Bolin and Wallin (2020) for exact parameterization of Dependence matrix.

Value

the precision matrix of the multivariate model

References

Bolin, D. and Wallin, J. (2020), Multivariate type G Matérn stochastic partial differential equation random fields. J. R. Stat. Soc. B, 82: 215-239. https://doi.org/10.1111/rssb.12351

Examples

library(fmesher)
library(fields)
# Define mesh
x <- seq(from = 0, to = 1, length.out = 40)
mesh <- fm_rcdt_2d_inla(lattice = fm_lattice_2d(x, x), extend = FALSE)
# Set parameters
p <- 3 # number of fields
rho <- c(-0.5, 0.5, -0.25) # correlation parameters
log_kappa <- list(2, 2, 2) # log(kappa)
variances <- list(1, 1, 1) # set marginal variances to 1
alpha <- list(2, 2, 2) # smoothness parameters
# Compute precision
Q <- precision_matrix_multivariate_spde(p,
  mesh = mesh, rho = rho,
  alpha = alpha, theta_K_list = log_kappa,
  variance_list = variances
)
# Plot the cross covariances
A <- as.vector(fm_basis(mesh, loc = matrix(c(0.5, 0.5), 1, 2)))
Sigma <- as.vector(solve(Q, c(A, rep(0, 2 * mesh$n))))
r11 <- Sigma[1:mesh$n]
r12 <- Sigma[(mesh$n + 1):(2 * mesh$n)]
r13 <- Sigma[(2 * mesh$n + 1):(3 * mesh$n)]
Sigma <- as.vector(solve(Q, c(rep(0, mesh$n), A, rep(0, mesh$n))))
r22 <- Sigma[(mesh$n + 1):(2 * mesh$n)]
r23 <- Sigma[(2 * mesh$n + 1):(3 * mesh$n)]
Sigma <- as.vector(solve(Q, v <- c(rep(0, 2 * mesh$n), A)))
r33 <- Sigma[(2 * mesh$n + 1):(3 * mesh$n)]

proj <- fm_evaluator(mesh)

oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 3))
image.plot(fm_evaluate(proj, r11), main = "Cov(X_1(s0),X_1(s)")
plot.new()
plot.new()
image.plot(fm_evaluate(proj, r12), main = "Cov(X_1(s0),X_2(s)")
image.plot(fm_evaluate(proj, r22), main = "Cov(X_2(s0),X_2(s)")
plot.new()
image.plot(fm_evaluate(proj, r13), main = "Cov(X_1(s0),X_3(s)")
image.plot(fm_evaluate(proj, r23), main = "Cov(X_2(s0),X_3(s)")
image.plot(fm_evaluate(proj, r33), main = "Cov(X_3(s0),X_3(s)")
par(oldpar)

Preconditioner SGD optimization

Description

Preconditioner SGD optimization

Usage

precond_sgd(stepsize = 1, numerical_eps = 1e-05)

Arguments

stepsize

stepsize for SGD
numerical_eps

numerical, the gap used for estimate preconditioner, default is 1e-5

Details

Preconditioner SGD is using fisher information matrix as preconditioner (natural gradient descent).

Value

a list of control variables for optimization (used in control_opt function)

Predict function of ngme2 predict using ngme after estimation

Description

Predict function of ngme2 predict using ngme after estimation

Usage

## S3 method for class 'ngme'
predict(
  object,
  map,
  data = NULL,
  type = "lp",
  group = NULL,
  estimator = c("mean", "sd", "0.05q", "0.95q", "median", "mode"),
  sampling_size = 500,
  burnin_size = 100,
  seed = Sys.time(),
  train_idx = NULL,
  chain_combine = c("param_mean", "predictive_average"),
  return_samples = FALSE,
  ...
)

Arguments

object

a ngme object
map

a named list (or dataframe) of the locations to make the prediction
data

a data.frame or matrix of covariates (used for fixed effects) names(loc) corresponding to the name each latent model vector or matrix (n * 2) for spatial coords
type

what type of prediction, c("fe", "lp", <model_name>) "fe" is fixed effect prediction <model_name> is prediction of a specific model "lp" is linear predictor (including fixed effect and all sub-models) "response" is the linear predictor plus a fresh measurement-noise draw
group

which filed to predict (used for bivariate model, should be of same length as map)
estimator

what type of estimator. Options include: - "mean", "median", "mode", "sd": standard estimators - "0.XXXq": any quantile specified as probability (e.g., "0.025q", "0.5q", "0.975q")
sampling_size

size of posterior sampling
burnin_size

size of posterior burnin
seed

random seed
train_idx

optional vector of training indices to use for posterior sampling. If provided, only these indices from the original data will be used for training, similar to cross-validation. If NULL, uses all original training data.
chain_combine

how to combine multiple optimization chains:

  • "param_mean": default behavior using the fitted object parameters.

  • "predictive_average": run prediction for each optimization chain and average predictions across chains.

return_samples

logical; when 'TRUE', attach sample draws for the requested output in 'attr(ret, "samples")'. For 'type = "response"', the attached samples are response predictive draws.
...

additional arguments (currently unused)

Value

a list of outputs contains estimation of operator paramters, noise parameters

Print an ngme model

Description

Print an ngme model

Usage

## S3 method for class 'ngme'
print(x, ...)

Arguments

x

ngme model object
...

...

Value

Invisibly returns x, an ngme object. The method is called for its side effect of printing a model summary.

Print ngme model

Description

Print ngme model

Usage

## S3 method for class 'ngme_model'
print(x, padding = 0, ...)

Arguments

x

ngme model object
padding

number of white space padding in front
...

...

Value

a list (model specifications)

Print ngme noise

Description

Print ngme noise

Usage

## S3 method for class 'ngme_noise'
print(x, padding = 0, prefix = "Noise type", model_type = NULL, ...)

Arguments

x

noise object
padding

number of white space padding in front
prefix

prefix
model_type

model type
...

...

Value

a list (noise specifications)

Print ngme operator

Description

Print ngme operator

Usage

## S3 method for class 'ngme_operator'
print(x, padding = 0, prefix = "Model type", ...)

Arguments

x

ngme operator object
padding

number of white space padding in front
prefix

prefix string
...

...

Value

a list (operator specifications)

Print ngme object

Description

Print ngme object

Usage

## S3 method for class 'ngme_replicate'
print(x, ...)

Arguments

x

ngme object
...

ignored

Value

Invisibly returns x, an ngme_replicate object. The method is called for its side effect of printing a replicate summary.

Print method for ngme_trajectories

Description

Print method for ngme_trajectories

Usage

## S3 method for class 'ngme_trajectories'
print(x, ...)

Arguments

x

ngme_trajectories object
...

additional arguments

Value

Invisibly returns x, an ngme_trajectories object. The method is called for its side effect of printing a summary of stored parameter trajectory matrices, including parameter names, chain count, and iteration count.

Print method for noise_kld_comparison

Description

Print method for noise_kld_comparison

Usage

## S3 method for class 'noise_kld_comparison'
print(x, ...)

Arguments

x

noise_kld_comparison object
...

additional arguments

Value

Invisibly returns x, a noise_kld_comparison object. The method is called for its side effect of printing the reference noise object, the Kullback-Leibler divergence values for each comparison noise, and the closest comparison.

Print method for parameter_distance

Description

Print method for parameter_distance

Usage

## S3 method for class 'parameter_distance'
print(x, ...)

Arguments

x

parameter_distance object
...

additional arguments

Value

Invisibly returns x, a numeric vector with class parameter_distance. The method is called for its side effect of printing the norm type, chain summary, number of iterations, final distance, and true parameter values.

Prior Half-Cauchy

Description

Prior Half-Cauchy

Usage

prior_half_cauchy(scale = 1, target = "coef")

Arguments

scale

half-Cauchy scale
target

apply prior on coefficient scale ('"coef"') or field scale ('"field"')

Value

prior specification

Prior Inverse-Exponential

Description

Prior induced by κ=1/ν∼Exp(λ)\kappa = 1 / \nu \sim \mathrm{Exp}(\lambda), giving p(ν)=λexp⁡(−λ/ν)ν−2p(\nu) = \lambda \exp(-\lambda / \nu)\nu^{-2} for ν>0\nu > 0. Internally this prior is applied to ν=lower+exp⁡(θ)\nu = \mathrm{lower} + \exp(\theta).

Usage

prior_inv_exponential(lambda = 1, lower = 0, target = "coef")

prior_inv_exp(lambda = 1, lower = 0, target = "coef")

Arguments

lambda

exponential rate on κ=1/ν\kappa = 1/\nu
lower

lower shift used in ν=lower+exp⁡(θ)\nu = \mathrm{lower} + \exp(\theta)
target

apply prior on coefficient scale ('"coef"') or field scale ('"field"')

Value

prior specification

Prior None

Description

Prior None

Usage

prior_none(target = "coef")

Arguments

target

apply prior on coefficient scale ('"coef"') or field scale ('"field"')

Value

prior specification

Prior Normal

Description

Prior Normal

Usage

prior_normal(mean = 0, sd = 1, target = "coef")

Arguments

mean

prior mean
sd

prior standard deviation
target

apply prior on coefficient scale ('"coef"') or field scale ('"field"')

Value

prior specification

Prior PC-SD

Description

Prior PC-SD

Usage

prior_pc_sd(u, alpha, target = "coef")

Arguments

u

tail event threshold
alpha

tail probability
target

apply prior on coefficient scale ('"coef"') or field scale ('"field"')

Value

prior specification

Prior Container

Description

Prior Container

Usage

priors(...)

Arguments

...

named prior specifications

Value

named prior container

ngme random effect model

Description

ngme random effect model

Usage

re(map, theta_K = NULL, ...)

Arguments

map

numerical vector, covariates to build index for the process (can be formula, provided data)
theta_K

initial value for theta_K (build covariance matrix)
...

currently ignored

Value

ngme_operator object

Root Mean Square Propagation (RMSProp) SGD optimization

Description

Root Mean Square Propagation (RMSProp) SGD optimization

Usage

rmsprop(stepsize = 0.05, beta1 = 0.9, epsilon = 1e-08)

Arguments

stepsize

stepsize for SGD
beta1

beta1 for momentum
epsilon

epsilon for numerical stability

Details

The update rule for RMSProp is:

vt=β1vt−1+(1−β1)gt2v_t = \beta_1 v_{t-1} + (1 - \beta_1) g_t^2

xt+1=xt−stepsize∗gtvt+ϵx_{t+1} = x_t - \text{stepsize} * \frac{g_t}{\sqrt{v_t} + \epsilon}

Value

a list of control variables for optimization (used in control_opt function)

Random Walk Model of Order 1 (RW1)

Description

Constructs a first-order random walk model for spatial or temporal processes. The RW1 model assumes that first-order differences ΔWi=Wi−Wi−1\Delta W_i = W_i - W_{i-1} are independent and identically distributed Gaussian variables.

Usage

rw1(mesh = NULL, cyclic = FALSE, constr = TRUE)

Arguments

mesh

numerical vector or inla.mesh.1d object, locations to build the mesh. For numerical vectors, assumes equally spaced locations.
cyclic

logical, whether the mesh is circular. If TRUE, the first and last locations are treated as neighbors. Cannot be FALSE when constr = TRUE.
constr

logical, whether to enforce the sum-to-zero constraint ∑i=1nhiWi=0\sum_{i=1}^n h_i W_i = 0. If FALSE, fixes the first element W1=0W_1 = 0. Must be TRUE for cyclic models.

Details

The RW1 model is defined by the precision matrix KK that penalizes first-order differences. The model has different structures depending on the constraints:

**Non-cyclic, constrained (default)**: The first row of KK enforces the constraint ∑i=1nhiWi=0\sum_{i=1}^n h_i W_i = 0, where hih_i are the mesh weights.

**Non-cyclic, unconstrained**: The first element is fixed at 0 (W1=0W_1 = 0).

**Cyclic**: Treats the domain as circular, connecting the first and last locations as neighbors. The constraint ∑i=1nhiWi=0\sum_{i=1}^n h_i W_i = 0 is always enforced. The precision matrix is expanded to size (n+1)×(n+1)(n+1) \times (n+1) to handle the constraint properly.

Value

An 'ngme_operator' object containing the precision matrix and related components for the RW1 model.

Examples

# Non-cyclic constrained RW1 (default)
rw1_default <- rw1(1:5)
print(rw1_default$K)

# Non-cyclic unconstrained RW1 (fixes first element to 0)
rw1_fixed <- rw1(1:5, constr = FALSE)
print(rw1_fixed$K)

# Cyclic RW1 (connects first and last locations)
rw1_cyclic <- rw1(1:5, cyclic = TRUE)
print(rw1_cyclic$K)

# Using with unequally spaced locations
locations <- c(0, 1, 3, 6, 10)
rw1_unequal <- rw1(locations)
print(rw1_unequal$K)

Random Walk Model of Order 2 (RW2)

Description

Constructs a second-order random walk model for spatial or temporal processes. The RW2 model assumes that second-order differences Δ2Wi=Wi−2Wi−1+Wi−2\Delta^2 W_i = W_i - 2W_{i-1} + W_{i-2} are independent and identically distributed Gaussian variables.

Usage

rw2(mesh = NULL, cyclic = FALSE, constr = TRUE)

Arguments

mesh

numerical vector or inla.mesh.1d object, locations to build the mesh. For numerical vectors, assumes equally spaced locations. Must have at least 3 locations.
cyclic

logical, whether the mesh is circular. If TRUE, the first and last locations are treated as neighbors with second-order differences computed across the boundary.
constr

logical, whether to enforce the sum-to-zero constraint ∑i=1nhiWi=0\sum_{i=1}^n h_i W_i = 0 and the linear trend constraint ∑i=1n(∑j=1ihj−h1)⋅Wi=0\sum_{i=1}^n (\sum_{j=1}^i h_j - h_1) \cdot W_i = 0. If FALSE, fixes the first and second elements W1=0W_1 = 0 and W2=0W_2 = 0.

Details

The RW2 model is defined by the precision matrix KK that penalizes second-order differences, making it smoother than RW1. The model enforces constraints to ensure identifiability:

**Non-cyclic (default)**: The first two rows of KK enforce constraints: the first row implements ∑i=1nhiWi=0\sum_{i=1}^n h_i W_i = 0 (sum-to-zero), and the second row implements ∑i=1nhi⋅i⋅Wi=0\sum_{i=1}^n h_i \cdot i \cdot W_i = 0 (removes linear trend). The remaining rows penalize second-order differences.

**Cyclic**: Treats the domain as circular, connecting the first and last locations as neighbors. The constraint ∑i=1nhiWi=0\sum_{i=1}^n h_i W_i = 0 is enforced. The precision matrix is expanded to size (n+1)×(n+1)(n+1) \times (n+1) to handle the constraint and circular structure properly.

Value

An 'ngme_operator' object containing the precision matrix and related components for the RW2 model.

Examples

# Non-cyclic RW2 with constraints (default)
rw2_default <- rw2(1:6)
print(rw2_default$K)

# Cyclic RW2 (connects first and last locations)
rw2_cyclic <- rw2(1:6, cyclic = TRUE)
print(rw2_cyclic$K)

# Using with unequally spaced locations
locations <- c(0, 1, 3, 6, 10, 15)
rw2_unequal <- rw2(locations)
print(rw2_unequal$K)

Vanilla SGD optimization

Description

Simple stochastic gradient descent without momentum or preconditioning.

Usage

sgd(stepsize = 0.001)

Arguments

stepsize

stepsize for SGD

Value

a list of control variables for optimization (used in control_opt function)

Stochastic Gradient Langevin Dynamics (SGLD) optimization

Description

Stochastic Gradient Langevin Dynamics (SGLD) optimization

Usage

sgld(stepsize = 0.001, temperature = 1)

Arguments

stepsize

base stepsize
temperature

non-negative Langevin temperature

Details

SGLD adds Gaussian noise to vanilla SGD updates:

xt+1=xt−ηt∇U(xt)+2Tηt ξt,ξt∼N(0,I)x_{t+1} = x_t - \eta_t \nabla U(x_t) + \sqrt{2 T \eta_t}\,\xi_t, \quad \xi_t \sim \mathcal{N}(0, I)

where TT is temperature. The implementation applies this component-wise using the current effective stepsize.

Value

a list of control variables for optimization (used in control_opt function)

Simulate from a ngme object (possibly with replicates)

Description

Simulate from a ngme object (possibly with replicates)

Usage

## S3 method for class 'ngme'
simulate(object, nsim = 1, seed = NULL, ...)

Arguments

object

ngme object
nsim

number of simulations
seed

seed
...

optional arguments. Supported names are posterior (whether to simulate from posterior sampling of latent fields) and m_noise (whether to add the measurement noise).

Value

a realization of ngme object

Simulate latent process with noise

Description

Simulate latent process with noise

Usage

## S3 method for class 'ngme_model'
simulate(object, nsim = 1, seed = NULL, ...)

Arguments

object

ngme model specified by f() function
nsim

number of simulations
seed

seed
...

ignored

Value

a realization of latent model

Examples

simulate(f(1:10, model = ar1(rho = 0.4), noise = noise_nig()))
simulate(f(rnorm(10), model = rw1(), noise = noise_normal()))
simulate(f(1:10, model = ar1(rho = 0.4), noise = noise_t(nu = 5)))

Simulate ngme noise object

Description

Simulate ngme noise object

Usage

## S3 method for class 'ngme_noise'
simulate(object, nsim = 1, seed = NULL, h = NULL, ...)

Arguments

object

ngme noise object
nsim

number of simulations
seed

seed
h

should be of same length as nsim
...

ignored

Value

data.frame (each col is a realization)

Ngme space-time non-separable model specification

Description

Implements a non-separable space-time model based on the advection-diffusion SPDE with first-order derivative in time. The model combines temporal and spatial components through a finite difference method (implicit Euler) for temporal discretization and finite element method (continuous Galerkin) for spatial discretization. When advection dominates diffusion, the "Streamline Diffusion" stabilization technique is applied.

Usage

spacetime(
  mesh,
  lambda = 1,
  alpha = 2,
  cc = 1,
  kappa = 1,
  fix_gamma = FALSE,
  theta_gamma_x = 0,
  theta_gamma_y = 0,
  shared_theta_gamma = FALSE,
  B_gamma_x = matrix(1, nrow = mesh[[2]]$n, ncol = 1),
  B_gamma_y = matrix(1, nrow = mesh[[2]]$n, ncol = 1),
  B_gamma_x_list = NULL,
  B_gamma_y_list = NULL,
  stabilization = TRUE
)

Arguments

mesh

A list of two objects:

  • mesh_t - The temporal mesh

  • mesh_s - The spatial mesh

lambda

The spatial damping parameter.
alpha

2 or 4, SPDE smoothness parameter.
cc

Parameter c in the SPDE.
kappa

Kappa parameter from Matern SPDE.
fix_gamma

TRUE if fix gamma (advection term), FALSE if estimate gamma.
theta_gamma_x

The x component of the advection term: gamma_x = B_gamma_x %*% theta_gamma_x.
theta_gamma_y

The y component of the advection term: gamma_y = B_gamma_y %*% theta_gamma_y.
shared_theta_gamma

TRUE if share the same theta_gamma for all time nodes. (theta_gamma_x and theta_gamma_y will be the same)
B_gamma_x

The design matrix for the x component of the advection term.
B_gamma_y

The design matrix for the y component of the advection term.
B_gamma_x_list

A list of design matrices for the x component of the advection term on every time node, length(B_gamma_x_list) == nt-1.
B_gamma_y_list

A list of design matrices for the y component of the advection term on every time node, length(B_gamma_y_list) == nt-1.
stabilization

TRUE if using a stabilization term (for implicit Euler).

Details

The model is particularly useful for large space-time datasets in environmental science, offering computationally efficient methods for parameter estimation, kriging prediction, and conditional simulations.

For details, see doi:10.1016/j.spasta.2024.100847

Value

ngme_operator object.

Unified stepsize control

Description

Unified stepsize control

Usage

stepsize_control(
  schedule = stepsize_schedule(method = "constant"),
  decay = stepsize_decay(method = "none")
)

Arguments

schedule

schedule component. Either an object from stepsize_schedule() or a method string accepted by stepsize_schedule().
decay

decay component. Either an object from stepsize_decay() or a method string accepted by stepsize_decay().

Details

Bundle an iteration schedule (stepsize_schedule) and an optional checkpoint-based decay rule (stepsize_decay) into one object for control_opt().

Value

a bundled stepsize-control object for control_opt().

Stepsize decay schedule

Description

Stepsize decay schedule

Usage

stepsize_decay(
  method = c("none", "grad_norm_plateau"),
  patience = 3,
  gamma = 0.5,
  min_delta = 0,
  warmup = 0,
  min_stepsize = 0
)

Arguments

method

decay strategy. "none" disables decay; "grad_norm_plateau" decays when mean grad.norm() across chains fails to decrease for a number of epochs.
patience

number of consecutive epochs without grad.norm() improvement before decaying stepsize.
gamma

decay factor applied when triggered (0 < gamma < 1).
min_delta

minimum required decrease in grad.norm() to be counted as improvement.
warmup

number of initial epochs to skip decay checks.
min_stepsize

lower bound for stepsize after decay (absolute value).

Details

Define a stepsize decay strategy for use in control_opt(). Currently supports the plateau-based schedule using grad.norm().

Value

a list of control variables for stepsize decay (used in control_opt).

Stepsize schedule

Description

Stepsize schedule

Usage

stepsize_schedule(
  method = c("constant", "poly"),
  alpha = 0.501,
  t0 = 1,
  burnin_iter = 0
)

Arguments

method

schedule type. "constant" keeps ηt=η0\eta_t = \eta_0 and "poly" uses polynomial decay.
alpha

polynomial exponent used when method = "poly". Must satisfy 1/2<α<11/2 < \alpha < 1.
t0

non-negative offset in iteration index.
burnin_iter

non-negative integer. During these initial iterations, schedule scaling is fixed to 1. Afterward, polynomial decay is applied with reset local time index.

Details

Define an iteration-dependent stepsize schedule for use in control_opt(). This schedule is independent of stepsize_decay() and applies a direct multiplicative scaling by iteration index:

ηt=η0(t+t0)−α\eta_t = \eta_0 (t + t_0)^{-\alpha}

.

Value

a list of control variables for stepsize schedule (used in control_opt).

Summary of ngme fit result

Description

Summary of ngme fit result

Usage

## S3 method for class 'ngme'
summary(object, name = NULL, ...)

Arguments

object

an object of class ngme
name

name of the latent model to be summarized (if NULL, will print all)
...

other arguments

Value

a list of summary

Summary for Batch-Means CI Results

Description

Summarize an object of class 'ngme_batch_ci', including parameter means, standard errors, confidence intervals, and covariance matrix.

Usage

## S3 method for class 'ngme_batch_ci'
summary(object, digits = 4, ...)

Arguments

object

an object of class 'ngme_batch_ci'.
digits

number of digits used for printing.
...

currently unused.

Value

A list with summary table and covariance matrix.

ngme tensor-product model specification

Description

Given 2 operator (first and second), build a tensor-product operator based on K = K_first x K_second (here x is Kronecker product)

Usage

tp(first, second)

Arguments

first

ngme_operator, left side of kronecker model (usually a temporal or iid model)
second

ngme_operator, right side of kronecker model (ususally a temporal or spatial model)

Value

a list of specification of model

Trace plot of ngme fitting

Description

Trace plot of ngme fitting

Usage

traceplot(
  ngme,
  name = "all",
  moving_window = 1,
  hline = NULL,
  combine = TRUE,
  ncol = NULL
)

Arguments

ngme

ngme object
name

name of latent models, otherwise plot fixed effects and measurement noise. Use "all" to plot all latent models plus data parameters in one figure. should be in names(ngme$models) or other
moving_window

moving window for the traceplot
hline

vector, add hline to each plot
combine

bool, if TRUE, return a single faceted ggplot; otherwise return a list of ggplot objects and print them one by one.
ncol

number of facet columns when combine = TRUE. If NULL (default), use 2 columns for multi-parameter plots and 1 for single-parameter plots.

Value

A ggplot object when combine = TRUE; otherwise a list of ggplot objects. The mean parameter trajectories are stored in the avg_lines attribute of the returned object.

ngme VAR(1) bivariate model specification (Cayley re-parameterization)

Description

Builds a Vector Autoregressive order-1 (VAR(1)) operator for bivariate time-series data. The model follows:

Yt=A Yt−1+εt,A=(a11a12a21a22)Y_t = A\,Y_{t-1} + \varepsilon_t, \quad A = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix}

Usage

var1(mesh = NULL, p1 = 0, p2 = 1, p3 = 0, p4 = 1)

Arguments

mesh

Integer vector or inla.mesh.1d object for the time mesh (length TT). Pass NULL for deferred construction.
p1

Unconstrained skew-symmetric parameter (∈R\in \mathbb{R}); controls the rotation part of the Cayley mapping. Default 0.
p2, p3, p4

Unconstrained lower-triangular Cholesky parameters (∈R\in \mathbb{R}); control the contraction (positive-definite) part. Defaults p2 = 1, p3 = 0, p4 = 1.

Value

An ngme_operator (or ngme_operator_def when mesh = NULL) suitable for use inside f().

Cayley re-parameterization

To guarantee stationarity (ρ(A)<1\rho(A) < 1) at every SGD step, the 2×22 \times 2 coefficient matrix AA is obtained via the Cayley transform:

A=(I+S)(I−S)−1A = (I + S)(I - S)^{-1}

where SS is constructed from four unconstrained real parameters (p1,p2,p3,p4)∈R4(p_1, p_2, p_3, p_4) \in \mathbb{R}^4 as S=J−RS = J - R:

Skew-symmetric part (rotation):

J=(0p1−p10)J = \begin{pmatrix} 0 & p_1 \\ -p_1 & 0 \end{pmatrix}

Positive-definite part (contraction):

R=LL⊤+εIR = L L^{\top} + \varepsilon I, where L=(p20p3p4)L = \begin{pmatrix} p_2 & 0 \\ p_3 & p_4 \end{pmatrix} and ε=10−5\varepsilon = 10^{-5}.

Because S+S⊤=−2RS + S^{\top} = -2R is negative definite, all eigenvalues of SS have strictly negative real parts, and the Cayley transform then guarantees ∣λi(A)∣<1|\lambda_i(A)| < 1 for every ii.

The default values (p1,p2,p3,p4)=(0,1,0,1)(p_1, p_2, p_3, p_4) = (0, 1, 0, 1) give A≈0A \approx 0 (no auto-regression at initialization).

Precision operator

The 2T x 2T precision operator is:

K=M0+a11 M11+a22 M22+a12 M12+a21 M21K = M_0 + a_{11}\,M_{11} + a_{22}\,M_{22} + a_{12}\,M_{12} + a_{21}\,M_{21}

where the MijM_{ij} are fixed sparse block matrices constructed from the T×TT \times T first-order difference matrix CTC_T.

Examples

set.seed(1)
n_obs <- 10
dat <- data.frame(
  y     = rnorm(2 * n_obs),
  idx   = rep(seq_len(n_obs), 2),
  group = factor(rep(c("y1", "y2"), each = n_obs))
)

fit <- ngme(
  y ~ 0 + f(idx, model = var1(mesh = 1:n_obs), group = group,
            noise = noise_nig()),
  data    = dat,
  family  = noise_normal(sigma = 0.01, fix_sigma = TRUE),
  control_opt = control_opt(estimation = FALSE)
)
print(fit)