Fix solver selection for one-dimensional Matérn models with free boundaries: their non-symmetric operators no longer use the symmetric solver, avoiding non-finite gradients in estimation and the cross-validation vignette.
Respect control_opt(max_num_threads = ) when it is smaller than the number
of chains, and avoid R console output from optimizer worker threads.
Stop estimation with an error if any optimizer step is non-finite, instead of silently returning unchanged parameters from a failed chain.
Replace the large fractional-model comparison in the automated tests with
small one- and two-dimensional fit-and-predict regression tests. The optional
INLA/inlabru comparison is retained in
tests/demo/fractional-inlabru-comparison.R.
Fix the random-intercept fixture in the factor-cache regression tests to build its design matrix from each group's observations.
Adding an option to select the reordering method for the solvers on OSX.
Fix the gradient of the fixed effects under non-Gaussian measurement noise. It
weighted observations by V / sigma^2 instead of 1 / (sigma^2 V), and
ignored the correlation of correlated measurement noise, biasing the
fixed-effect estimates.
Fix parameter names being out of order: par_names, the convergence messages
and attr(fit, "conv_diag") gave the measurement-noise parameters and the
fixed effects each other's names.
Fix the cap on the measurement noise nu, which set log(nu) to 1e4
instead of log(1e4).
precond_sgd(): the preconditioner for the fixed effects and the measurement
noise mu now integrates out the latent field, instead of using the
complete-data curvature, which overstated it many times over under a
correlated field.
precond_sgd() uses the marginal Fisher information for the sigma of
non-Gaussian measurement noise, which lets NIG measurement noise converge
where it previously stalled. New control_opt(precond_meas_sigma = , fisher_refresh_every = ) select this curvature ("auto", "fisher" or
"complete") and how often it is recomputed.
control_opt(nig_param_std = ) now applies to NIG measurement noise; it is
not applied to latent NIG noise. For such fits attr(fit, "chain_params")
and conv_diag hold the optimiser's coordinates for mu, sigma and nu.
Remove moving_window from ngme(). It was declared and documented ("return
the average estimation of last .. iterations") but never read anywhere in the
function, so it silently did nothing.
spacetime() gains cc_variance_free, DEFAULT TRUE. The operator now puts
the whole cc factor on the temporal term (K = cc*BtCs + Ls) instead of
splitting it sqrt(cc) between the temporal and spatial terms. This makes the
variance less dependent on cc and typically improves mixing.
cross_validation() accepts a NAMED list of transforms and scores the same
posterior draws on each scale. transform is applied after the draws exist
and never affects sampling, so calling once per scale re-drew identical
samples and doubled the cost of a cross-validation. Verified to reproduce the
two-call scores to 10 significant figures. An unnamed list keeps its previous
meaning of one entry per model.
Fix n_burnin being silently dropped in cross_validation().
Live diagnostics ([iter ...], [stationarity], [polish]) now go to
stderr via REprintf rather than stdout.
The polish phase's geometric decay now starts from the step size the search
was using. polish_stepsize_factor defaults to 1.0.
Trajectory continuation now requires the same model.
Parallelise the gradient computation over replicates.
Fix the advection term of the non-separable spacetime() model, which was
quadratic in the advection field and simply wrong. This changes fitted results
for every non-separable model with free advection.
Fix spacetime() assembling a different operator in C++ than in R.
spacetime() now computes its operator traces from the diagonal blocks,
needing neither Hutchinson probes nor a factorization of the full
(nt*ns) x (nt*ns) operator.
The spacetime() operator now starts from its approximatestationary distribution
instead of in zero. spacetime(stationary_init = FALSE) restores the old first block for
reproducing earlier fits.
Convergence is now decided in two phases, which separates two questions the
old rule ran together. control_opt(conv_criterion = ) selects between them
and defaults to the new "stationarity".
mc_se_lim * se_stat, i.e. when the Monte Carlo
error is a fraction of the statistical error the data can support.Very short runs now report non-convergence rather than claiming it. The stationarity test needs enough sub-batch history to split into two halves, so a fit given only a few dozen iterations warns instead of declaring success.
The stationarity test's comparison window is derived from the measured
correlation length of the chains rather than being a fixed count
(control_opt(stationarity_window = 0, stationarity_eff = )). tau is
estimated as the ratio of the batch-means variance to the naive variance of
the same mean and the window is 2 * stationarity_eff * tau.
The optional trend-triggered step-size decay
(stepsize_decay(method = "on_trend")) used the same |t| <= trend_lim test
to decide the chains had stopped moving, and so would decay soonest when the
window was noisiest. It now uses the drift rate as well.
ngme(start = previous_fit) is now a true continuation. start_sd is no
longer applied on a warm start at all. On top of that, each chain resumes
from its own final parameters instead of from the chain-averaged fit.
Each chain's endpoint is now stored on the fit (attr(fit, "chain_params"))
and handed back to the optimizer, and the storedoptimisation trajectory is
concatenated across the restarts rather than replaced, so traceplot() and any
convergence diagnostic see the whole path.
Each chain's own latent state (W, and the mixing variables V) is carried
over too: start holds the chain average of those.
control_opt(continue_chains = FALSE) restores the previous behaviour.
Convergence gained a third clause: the estimate must be precise, not merely
stationary. R-hat and the trend test say the chains reached the same stationary
distribution and stopped drifting; neither says the answer is pinned down, because
constant-step SGD settles into a distribution whose width is set by the step size
rather than by the data. control_opt(mc_se_conv_check = , mc_se_lim = ) requires the
Monte Carlo error of the reported estimate to be at most mc_se_lim (default 0.5) of
se_stat = sqrt(diag(P^-1)), the standard error the data supports, read from the
preconditioner.
Measurement-noise parameter names are prefixed meas_. A latent noise and the
measurement noise both have mu/sigma/nu, so par_names could contain sigma_1
twice and anything keyed by name silently merged two unrelated parameters.
R-hat is now computed over the same window as the trend test (sub-batch means
spanning the second half of the run) instead of over the current batch alone, and
n_conv_batch defaults to 2 rather than 1.
The convergence trend test gained a magnitude floor, control_opt(trend_rel_lim = )
(default 0.01). A parameter now fails the test only when its drift is both
statistically detectable (|slope|/se > trend_lim) and materially fast, moving more
than trend_rel_lim of its own scale (|mean| + sd) per 100 iterations. It is a rate,
not a movement per regression window.
The trend window is now built from sub-batch means spread evenly over the second half
of the run so far, instead of the last n_slope_check whole checkpoints. It is
therefore available at the first checkpoint rather than the n_slope_check-th.
n_min_batch now defaults to 1 rather than min(n_batch, 3), since it is no
longer what holds a converged fit back.
With verbose = TRUE and print_check_info = FALSE, each convergence checkpoint now
prints a progress summary: how many parameters have converged, and the worst
R-hat, t-statistic and drift rate with the parameter responsible.
predict() gained a replicate argument. Every prediction previously came
from the first replicate whatever the caller asked for, since the latent field
W is per-replicate. The argument takes a replicate level (the name used in
the fit) or a position, and an unmatched value is an error listing what the
fit actually has.
A namespace-qualified ngme2::f(...) in a model formula is now recognised as
a latent term. terms.formula(specials = ) matches specials by name, so the
qualified call fell through to the fixed-effect path and the latent field was
dropped from the model silently. The qualification is stripped from the
function position only, and the bare name is bound for evaluation, so this
keeps working when ngme2 is not attached.
Fixed f(model = ngme2::<constructor>(...), noise = noise_normal()), which
failed with "the condition has length > 1". f() read the constructor name as
as.character(model_expr[[1]]), which is c("::", "ngme2", "tp") for a
qualified call, so the %in% test that follows was not a single logical.
Removed the dead preconditioner plumbing from control_opt()
(precond_strategy, precond_by_diff_chain, compute_precond_each_iter).
precond_sgd() has not carried a preconditioner field since the
preconditioner became the analytic second-order Hessian, so control_opt()
overwrote its own default with NULL and passed precond_strategy = numeric(0) to C++ -- which never read it. numerical_eps is unaffected.
The default nu prior for stationary NIG latent noise is now the penalised-complexity
prior of Cabral, Bolin and Rue (2023) calibrated as in their Section 3.3, from
P(1/nu > U) = alpha with U = 2.5 and alpha = 0.01, i.e. P(nu < 0.4) = 0.01.
A failed Cholesky of the measurement precision QQ no longer aborts the fit outright.
Such a failure is nearly always a transient SGD excursion so the diagonal is now
nudged by an escalating amount (from 1e-10 of the diagonal scale, six attempts) and the
factorization retried. Only if that fails does estimation stop, and the error now reports
the iteration, the smallest diagonal entry and the asymmetry instead of printing them to
the console from inside the parallel SGD region.
sf and inlabru are declared in Suggests; spacetime() now checks for sf
up front rather than failing inside sf::st_coordinates().
The Parana section of the introductory vignette is skipped when INLA cannot be loaded,
instead of failing the whole vignette.
Rao-Blackwellisation is now on by default (control_opt(rao_blackwellization = )).
Fixed the Gibbs sweep order that made Rao-Blackwellisation unusable with a non-Gaussian
family. The measurement V was redrawn after W, which left QQ and cond_W
conditioned on the previous V while the gradient used the new one; the fit diverged until
QQ lost positive-definiteness. Both V blocks are now drawn before W. Gaussian
families were unaffected, since there V is fixed.
The Rao-Blackwell traces are taken from an exact selected (Takahashi) inverse when the
Cholesky factor of QQ is cheap enough, and from Hutchinson probes otherwise, controlled
by control_opt(selinv_max_fill = ). Triangular operators (ar1, ou, arma) and 1d
meshes qualify; a 2d mesh does not. For those models the traces are now exact and the fit
no longer depends on n_trace_iter at all.
The Hutchinson probe count adapts during the run (control_opt(trace_adapt = )) so that
the trace estimator contributes a target share of the gradient variance instead of being
fixed a priori.
The convergence trend test is now scale-free: it compares the fitted slope to its own
standard error rather than to an absolute bound, so trend_lim is a t-value (default 2)
and no longer depends on the parameter's units or on n_batch.
A run that exhausts its iteration budget without converging now warns and names the
parameters still failing, instead of returning silently. Set
control_opt(warn_no_convergence = FALSE) for short runs where convergence is not expected.
Removed the Pflug convergence diagnostic (pflug_conv_check, pflug_alpha). It could
never fire: the rule looks for sign changes in the running sum of <g_t, g_{t-1}>, but the
gradient is evaluated on a Gibbs sample of the latent field, so consecutive gradients stay
positively correlated however close the optimum is.
src/Makevars now tracks header dependencies, so editing a header rebuilds every object
that includes it. Previously only each source file was compared against its own object,
which could leave objects compiled against two different layouts of the same class.
Fixed single_V gradients for the NIG and t noise types, which did not match the score of
the corresponding density and had the wrong sign over part of the range.
The "ws" replicate sampling strategy now applies the importance weight it was missing, so
it no longer differs in scale from "all" by roughly the number of replicates.
Update estimation code to only reassemble and refactorize the latent precision
matrix QQ when it can actually have changed. The matrix is also split into
K' diag(1/SV) K and its measurement block (AZ)' D (AZ), and these are
caches separately. Results are unchanged by the caching and setting the
environment variable NGME2_DISABLE_FACTOR_CACHE restores the uncached behaviour.
Fixed the trace estimators for non-symmetric operators, which made the
preconditioner Hessian incorrect and the first-order trace far noisier than
necessary. control_opt(nonsym_solver = ) chooses the factorization for this case.
"normal_equations" (default) for a Cholesky of t(K) K, or "lu" for a sparse LU
of K. The Cholesky is cheaper for every non-symmetric operator in the package.
Use "lu" if K is so ill-conditioned that the Cholesky of t(K) K fails.
A failed factorization of K for a non-symmetric operator is now reported for
both routes, not only for the LU.
Speed-ups for space-time models. spacetime() builds its spatial block once per
assembly of K rather than once per time node: with the default B_gamma_x /
B_gamma_y the advection field is the same at every time node, so all nt - 1
blocks are the same matrix. The Rao-Blackwell traces no longer assemble the
matrix they trace either: tr(QQ^-1 A' diag(d) B) is now taken as
A' (d * (B QU)) on the probe block instead of forming the triple product,
which is as dense as QQ and was built once per parameter per iteration.
NGME2_NO_FACTORED_TRACE restores the assembled trace for comparison.
Increased the burnin for the Gibbs chain for the posterior samples of W and V
that ngme() returns from 1 to 100, and before a posterior simulate() draw.
Default random seeds are now drawn from the ambient R random number stream
instead of the system clock, so set.seed() makes simulate(),
control_opt(), ngme(), predict(), cross_validation(),
compute_ngme_ci(), compute_ngme_sgld_samples() and test_ngme()
reproducible in the usual R way.
Restore AR1 model to use stationary initial conditions and fix all models with triangular operators to use exact trace calculations.
pm25_quarterly_2022, a processed example dataset of quarterly
PM2.5 monitor averages from EPA AQS/AirData for the contiguous United States
in 2022. The dataset includes monitor location, quarter, quarterly average
PM2.5, contributing day counts, and a replicate index for quarterly models.noise_moments() for stationary ngme_noise objects, with direct
helpers noise_nig_moments() and noise_gal_moments() returning variance,
standard deviation, skewness, kurtosis, excess kurtosis, and central moments.normal_nig, normal_gal,
and nig_gal) by expanding prediction projection matrices when required and
using the output-sample prediction path for non-mean estimators.ngme() now respects
control_opt(verbose = ...) more consistently and uses messages instead of
unconditional cat() output.traceplot() no longer prints last estimates as a
side effect.simulate.ngme() compatible with the S3 simulate() interface by
accepting posterior and m_noise through ....get_data_from_formula() and improve diagnostics that depend on
captured model/design output.var1() — a Vector Autoregressive order-1 latent field for bivariate
time series modelling.RCallback in
src/latents/rcallback.cpp).ngme2 noise distributions (NIG, GAL, normal) as a single
shared innovation noise.print() method displays the recovered $A$ matrix, its spectral radius, and
the raw $(p_1, p_2, p_3, p_4)$ values.vignettes/var1-model.Rmd)
covering: model specification, Cayley reparameterization, simulation study
with parameter recovery, convergence trace plots, and NIG vs Gaussian
model comparison.posterior_plot().plot() support for ngme_sgld_ci objects, reusing stored SGLD
samples to visualize marginal posterior distributions.~ 0 + ...): skip fixed-effect
centering when no intercept is present.fe() centering with structural zeros: grouped fe() columns are
centered using in-group rows only, so out-of-group structural zeros remain
zero.data_idx) instead of always
taking the first n rows.start = previous_fit) across
standardization settings by remapping fixed effects through the current model
parameterization."(Intercept)"*) now default to prior_none(), while non-intercept
columns keep the default N(0,10) prior.standardize_fixed = TRUE, add prior compatibility handling:
isotropic normal priors on standardized columns are transformed to the SVD
basis; incompatible custom prior_beta specifications now automatically
disable fixed-effect standardization with a warning.prior_inv_exponential(lambda, lower) for nu, implementing
kappa = 1 / nu ~ Exp(lambda) as a first-class prior option.prior_inv_exp(...) for the same prior.calibrate_inv_exp_lambda_driven_nig() and
calibrate_inv_exp_lambda() for choosing lambda from a driven-noise
tail-inflation target.R_c(nu) curves: the
helper now scans for crossings and reports observed R_c range when the
requested target is unattainable.nu prior in f() for NIG-driven noise: when nu prior is
not explicitly set and nu is stationary, use
prior_inv_exp(lambda = log(2)/median(h), lower = nu_lower_bound).
For non-stationary nu, keep the legacy N(0,10) default prior.ngme() estimation/sampling path: C++ exceptions
are now propagated as R errors (including OpenMP parallel regions) instead of
potentially terminating the R session.nu initialization in noise helper constructors to respect
nu_lower_bound, using theta_nu = log(nu - nu_lower_bound) and validating
nu > nu_lower_bound.normal_nig conversion, printing, and plotting with effective
parameterization nu = nu_lower_bound + exp(theta_nu).prior_normal(), prior_pc_sd(), prior_half_cauchy(), prior_none(), and priors(...).f() and ngme_noise() to accept unified prior = ... inputs
(remove prior_theta_K and prior_mu/prior_sigma/prior_nu arguments).ngme_prior() interface and its documentation entry.coef/field) for noise parameter priors and
per-parameter operator prior compilation.ngme(..., prior_beta = ...), using the
same prior_*()/priors(...) API.Prior Templates for Stationary and Non-Stationary Models.control_opt(stepsize_decay = "grad_norm_plateau") (epoch-level, synchronized across chains)stepsize_decay() helper for configuring decay optionscross_validation(data = ...) model rebuild for refit-on-new-data workflows:
it now resolves external formula symbols (for example mesh, B, n_basis) from the fitted object when needed, and falls back to rebuild-without-start plus hyperparameter transplant if start state dimensions differ.chain_combine = "predictive_average" in predict() and cross_validation(), which averages predictions across optimization chains instead of averaging parameters first.control_opt to specify the solver