Inference and prediction for mixed effects models with flexible non-Gaussian and Gaussian 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.

Author(s)

Maintainer: Xiaotian Jin xiaotian.jin@kaust.edu.sa

Authors:

• David Bolin davidbolin@gmail.com [copyright holder]

• Alexandre Simas alexandre.impa@gmail.com

• Jonas Wallin jonas.wallin81@gmail.com

Other contributors:

• Andrea V. Rocha dearocha@yahoo.com (contributed to many vignettes and examples) [contributor]

• Timothy A. Davis (SuiteSparse components) [contributor, copyright holder]

• Patrick R. Amestoy (SuiteSparse AMD and CAMD components) [contributor, copyright holder]

• Iain S. Duff (SuiteSparse AMD and CAMD components) [contributor, copyright holder]

• John K. Reid (SuiteSparse AMD-derived code) [contributor, copyright holder]

• Yanqing Chen (SuiteSparse CAMD components) [contributor, copyright holder]

• Sivasankaran Rajamanickam (SuiteSparse CCOLAMD components) [contributor, copyright holder]

• Stefan Larimore (SuiteSparse CCOLAMD and COLAMD components) [contributor, copyright holder]

• William W. Hager (SuiteSparse CHOLMOD Modify components) [contributor, copyright holder]

• University of Florida (SuiteSparse CHOLMOD and CCOLAMD components) [copyright holder]

• Regents of the University of Minnesota (METIS/GKlib components bundled with SuiteSparse) [copyright holder]

• Free Software Foundation, Inc. (GKlib regex/getopt code bundled with SuiteSparse) [copyright holder]

• Makoto Matsumoto (Mersenne Twister code bundled with SuiteSparse) [contributor, copyright holder]

• Takuji Nishimura (Mersenne Twister code bundled with SuiteSparse) [contributor, copyright holder]

• Alexander Chemeris (MS integer headers bundled with SuiteSparse) [contributor, copyright holder]

See Also

Useful links:

• https://davidbolin.github.io/ngme2/

• https://github.com/davidbolin/ngme2

• Report bugs at https://github.com/davidbolin/ngme2/issues

Helper function to compute AR(p) autocovariance

Description

Helper function to compute AR(p) autocovariance

Usage

ARcov(rho, p)

Arguments

Value

vector of autocovariances from lag 0 to lag p-1

AdaGrad SGD optimization

Description

AdaGrad SGD optimization

Usage

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

Arguments

Details

The update rule for AdaGrad is:

v_t = v_{t-1} + g_t^2

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

Details

The update rule for Adam is:

m_t = \beta_1 m_{t-1} + (1 - \beta_1) g_t

v_t = \beta_2 v_{t-1} + (1 - \beta_2) g_t^2

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

\hat{v_t} = v_t / (1 - \beta_2^t)

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

Details

The update rule for AdamW is:

m_t = \beta_1 m_{t-1} + (1 - \beta_1) g_t

v_t = \beta_2 v_{t-1} + (1 - \beta_2) g_t^2

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

\hat{v_t} = v_t / (1 - \beta_2^t)

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

Details

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

\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})||} )

x_{k+1} = x_k - \lambda_k \nabla f(x_k)

\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

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

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

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

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

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

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

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

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(\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

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(\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

Value

a list of specification of model

Calibrate Inverse-Exponential Prior for NIG Driven Noise

Description

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

The calibration target is

\Pr(R_c(\nu) > r_\text{target}) = \alpha,

where

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

and

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 R_c(\nu_r) = r_\text{target} using log-scale bisection, then returns

\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

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

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

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

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

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

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

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

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 = "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

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

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

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

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

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

Details

The generalized asymmetric Laplace distribution has density given by

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 K_p is modified Bessel function of the second kind of order p, x>0, \nu>0 and \mu,\delta, \sigma\in\mathbb{R}. See Barndorff-Nielsen (1977, 1978 and 1997) for further details.

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

V \sim \Gamma(h \nu, \nu),

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

-\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

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

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

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)

get the parameters of the noise

Description

Get the parameters of the noise

Usage

get_noise_param(noise, prefix)

Arguments

Value

a list with parameter name as key and parameter value as value

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

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

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

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

Details

The generalised inverse-Gaussian distribution has density given by

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

where K_p is modified Bessel function of the second kind of order p, x>0, a,b>0 and p\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

Details

The inverse-Gaussian distribution has density given by

f(x; a, b) = \frac{\sqrt{b}}{\sqrt{2\pi x^3}}\exp( -\frac{a}{2}x -\frac{b}{2x} + \sqrt{ab}),

where x>0 and a,b>0. In this parameterization, E(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

Details

The inverse-Gamma distribution has density given by

f(x; a, b) = \frac{b^a}{\Gamma(a)}x^{-a-1}\exp( -\frac{b}{x}),

where x>0 and a,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

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

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

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

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

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

Momentum SGD optimization

Description

Momentum SGD optimization

Usage

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

Arguments

Details

The update rule for momentum is:

v_t = \beta_1 v_{t-1} + \beta_2 g_t

x_{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), inverse is log(x)/4

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

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

• identity: Identity function, no transformation

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

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

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

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

Usage

name2fun(trans, inv = FALSE)

Arguments

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

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

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

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

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

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

Details

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

\mu = B_{\mu} \theta_{\mu},

and

\sigma = \exp (B_{\sigma} \theta_{\sigma}),

\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

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

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

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

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

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

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))

Check whether a newer stable version of ngme2 is available

Description

This function checks the package repository for the latest available ngme2 version. It does not install or update packages.

Usage

ngme_update()

Value

Invisibly returns a list with the local version, remote version, repository URL, and a logical update_available flag.

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

Details

The normal inverse-Gaussian distribution has density given by

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 K_p is modified Bessel function of the second kind of order p, x>0, \nu>0 and \mu,\delta, \sigma\in\mathbb{R}. See Barndorff-Nielsen (1977, 1978 and 1997) for further details.

The additional parameter h is used when

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

. By the infinite divisibility,

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

. Then

\delta+\mu V + \sigma \sqrt{V} Z

has the distribution of

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))

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 \mathbf{K}\mathbf{X} = \boldsymbol{\delta}.

Usage

ou(mesh = NULL, theta = 1)

Arguments

Details

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

dX_t = -\theta X_t dt + \sigma dW_t

where \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 t_i and t_{i+1} with step size \Delta t_i:

X_{i+1} = X_i e^{-\theta \Delta t_i} + \epsilon_i

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

\rho_i = e^{-\theta \Delta t_i}

and Gaussian noise:

\epsilon_i \sim \mathcal{N}\left(0, \frac{\sigma^2}{2\theta}(1 - \rho_i^2)\right)

## Matrix Representation

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

\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 \rho_i = \exp(-\theta \Delta t_i) with \Delta t_i = t_{i+1} - t_i.

**Why the \sqrt{1-\rho_1^2} term?** The first row ensures that X_1 is drawn from the stationary distribution \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 \rho_i is computed locally based on the varying time step \Delta t_i, allowing the model to naturally handle irregularly spaced observations.

## The h vector

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

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 \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 \mathbf{K} is constructed such that the resulting process \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

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

Value

A ggplot object showing the parameter-distance trajectory over optimization iterations. The y-axis is the stored distance ||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

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

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

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

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

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

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

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

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

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

Value

a list (operator specifications)

Print ngme object

Description

Print ngme object

Usage

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

Arguments

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

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

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

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

Value

prior specification

Prior Inverse-Exponential

Description

Prior induced by \kappa = 1 / \nu \sim \mathrm{Exp}(\lambda), giving p(\nu) = \lambda \exp(-\lambda / \nu)\nu^{-2} for \nu > 0. Internally this prior is applied to \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

Value

prior specification

Prior None

Description

Prior None

Usage

prior_none(target = "coef")

Arguments

Value

prior specification

Prior Normal

Description

Prior Normal

Usage

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

Arguments

Value

prior specification

Prior PC-SD

Description

Prior PC-SD

Usage

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

Arguments

Value

prior specification

Prior Container

Description

Prior Container

Usage

priors(...)

Arguments

Value

named prior container

ngme random effect model

Description

ngme random effect model

Usage

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

Arguments

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

Details

The update rule for RMSProp is:

v_t = \beta_1 v_{t-1} + (1 - \beta_1) g_t^2

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 \Delta W_i = W_i - W_{i-1} are independent and identically distributed Gaussian variables.

Usage

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

Arguments

Details

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

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

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

**Cyclic**: Treats the domain as circular, connecting the first and last locations as neighbors. The constraint \sum_{i=1}^n h_i W_i = 0 is always enforced. The precision matrix is expanded to size (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 \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

Details

The RW2 model is defined by the precision matrix K 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 K enforce constraints: the first row implements \sum_{i=1}^n h_i W_i = 0 (sum-to-zero), and the second row implements \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 \sum_{i=1}^n h_i W_i = 0 is enforced. The precision matrix is expanded to size (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

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

Details

SGLD adds Gaussian noise to vanilla SGD updates:

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 T 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

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

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

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

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

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

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

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:

\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

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

Value

A list with summary table and covariance matrix.

Test ngme function

Description

Test ngme function for different models

Usage

test_ngme(
  model,
  n_obs_per_rep,
  n_replicate = 1,
  control_opt = NULL,
  numer_grad = TRUE,
  f_noise = noise_nig(mu = -2, sigma = 1.5, nu = 0.5),
  n_gibbs_samples = 5,
  family = "normal",
  max.n = 1000,
  debug = FALSE,
  debug_f = FALSE,
  start = NULL,
  fix_theta_K = FALSE,
  fix_theta_mu = FALSE,
  fix_theta_sigma = FALSE,
  fix_theta_nu = FALSE,
  seed = Sys.time()
)

Arguments

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

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

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:

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

Value

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

Cayley re-parameterization

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

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

where S is constructed from four unconstrained real parameters (p_1, p_2, p_3, p_4) \in \mathbb{R}^4 as S = J - R:

Skew-symmetric part (rotation):

J = \begin{pmatrix} 0 & p_1 \\ -p_1 & 0 \end{pmatrix}

Positive-definite part (contraction):

R = L L^{\top} + \varepsilon I, where L = \begin{pmatrix} p_2 & 0 \\ p_3 & p_4 \end{pmatrix} and \varepsilon = 10^{-5}.

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

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

Precision operator

The 2T x 2T precision operator is:

K = M_0 + a_{11}\,M_{11} + a_{22}\,M_{22} + a_{12}\,M_{12} + a_{21}\,M_{21}

where the M_{ij} are fixed sparse block matrices constructed from the T \times T first-order difference matrix C_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)