generate_posterior: Generate Posterior Samples from a Fitted lgspline

In lgspline: Lagrangian Multiplier Smoothing Splines for Smooth Function Estimation

Generate Posterior Samples from a Fitted lgspline

Description

Draws from the posterior distribution of model coefficients, with optional dispersion sampling, posterior predictive draws, and propagation of uncertainty in estimated correlation parameters.

Usage

generate_posterior(
  object,
  new_sigmasq_tilde = object$sigmasq_tilde,
  new_predictors = NULL,
  theta_1 = 0,
  theta_2 = 0,
  posterior_predictive_draw = function(N, mean, sqrt_dispersion, ...) {
     rnorm(N,
    mean, sqrt_dispersion)
 },
  draw_dispersion = TRUE,
  include_posterior_predictive = FALSE,
  num_draws = 1,
  enforce_qp_constraints = TRUE,
  draw_correlation = FALSE,
  correlation_param_mean = NULL,
  correlation_param_vcov = NULL,
  correlation_VhalfInv_fxn = NULL,
  correlation_Vhalf_fxn = NULL,
  correlation_param_vcov_scale = NULL,
  include_warnings = TRUE,
  ...
)

Arguments

object	A fitted lgspline model object.
new_sigmasq_tilde	Numeric; dispersion \tilde{\sigma}^2 used as the point estimate when draw_dispersion = FALSE. Default object$sigmasq_tilde.
new_predictors	Matrix; predictor matrix for posterior predictive sampling. Default uses in-sample predictors.
theta_1	Numeric; shape increment for the inverse-gamma prior on \sigma^2. Default 0.
theta_2	Numeric; rate increment for the inverse-gamma prior. Default 0.
posterior_predictive_draw	Function; sampler for posterior predictive realisations. Must accept N, mean, sqrt_dispersion, .... Default rnorm.
draw_dispersion	Logical; sample \sigma^2 from its posterior. Default TRUE.
include_posterior_predictive	Logical; generate posterior predictive draws at new_predictors. Default FALSE.
num_draws	Positive integer; number of draws. Default 1.
enforce_qp_constraints	Logical; if TRUE, enforce active QP inequality constraints during posterior sampling via the stored elliptical-slice constrained sampler. Default TRUE.
draw_correlation	Logical; propagate correlation parameter uncertainty. Requires VhalfInv_fxn and VhalfInv_params_estimates in the fitted object. Default FALSE.
correlation_param_mean	Numeric vector; mean of the approximate normal posterior for correlation parameters on the unbounded (working) scale. Default: object$VhalfInv_params_estimates.
correlation_param_vcov	Matrix; variance-covariance for correlation parameter draws. Default: inverse Hessian from BFGS (object$VhalfInv_params_vcov).
correlation_VhalfInv_fxn	Function; maps correlation parameter vector to \mathbf{V}^{-1/2}. Default object$VhalfInv_fxn.
correlation_Vhalf_fxn	Function or NULL; maps to \mathbf{V}^{1/2}. Passed through to generate_posterior_correlation; the current correlation-aware posterior path only requires correlation_VhalfInv_fxn.
correlation_param_vcov_scale	NULL or numeric; if supplied, divides a user-supplied correlation_param_vcov by this value before passing it to generate_posterior_correlation. When NULL, no additional scaling is applied.
include_warnings	Logical; emit warnings for degenerate draws, constraint violations, etc. Default TRUE.
...	Additional arguments forwarded to the GLM weight function, dispersion function, and posterior_predictive_draw.

Details

Uses a Laplace approximation centred at the MAP estimate for non-Gaussian responses.

Dispersion posterior. When draw_dispersion = TRUE, \sigma^2 is drawn from

\sigma^2 \mid \mathbf{y} \sim \mathrm{InvGamma}(\alpha_1, \alpha_2),

where

\alpha_1 = \theta_1 + \tfrac{1}{2}(N - s \cdot \mathrm{tr}(\mathbf{H})), \quad \alpha_2 = \theta_2 + \tfrac{1}{2}(N - s \cdot \mathrm{tr}(\mathbf{H})) \tilde{\sigma}^2,

\mathbf{H} = \mathbf{XUGX}^\top is the hat matrix, s = 1 when unbias_dispersion = TRUE (else s = 0), and \theta_1 = -1, \theta_2 = 0 recovers an improper uniform prior on \sigma^2, while \theta_1 = 0, \theta_2 = 0 gives the usual scale-invariant improper prior p(\sigma^2) \propto 1/\sigma^2.

Correlation parameter posterior. When draw_correlation = TRUE and the fitted model contains an estimated correlation structure, each draw first samples \boldsymbol{\rho} from

\boldsymbol{\rho}^{(m)} \sim \mathcal{N}(\hat{\boldsymbol{\rho}}_{\mathrm{REML}}, \mathbf{H}^{-1}_{\mathrm{BFGS}}),

rebuilds the posterior covariance under the drawn correlation structure (reusing all pre-computed design matrices, constraints, and penalty matrices) and then draws coefficients from the updated posterior. Knot placement, partitioning, coefficient re-estimation, and GCV tuning are skipped entirely. Draws producing non-positive-definite correlation matrices are rejected and redrawn (up to 50 attempts).

When draw_correlation = FALSE (default), correlation parameters are fixed at their estimated values.

Inequality constraints. Active QP inequalities can be enforced during posterior sampling via elliptical slice sampling, producing draws from the corresponding truncated multivariate normal posterior on the coefficient scale. The public enforce_qp_constraints argument is forwarded to the stored sampler for both the standard and correlation-aware posterior paths.

Value

When num_draws = 1, a named list:

post_draw_coefficients

List of length K+1; per-partition coefficient vectors on the original scale.

post_draw_sigmasq

Drawn (or fixed) dispersion.

post_pred_draw

Posterior predictive vector (only when include_posterior_predictive = TRUE).

post_draw_correlation_params

Drawn correlation parameters on the working scale (only when draw_correlation = TRUE).

When num_draws > 1, each element becomes a list of length num_draws, and post_pred_draw (if requested) is an N_{\mathrm{new}} \times M matrix, where M = \mathrm{num\_draws}.

See Also

lgspline, generate_posterior_correlation, wald_univariate

Examples

set.seed(1234)
n_blocks <- 100; block_size <- 5; N <- n_blocks * block_size
rho_true <- 0.3
t <- seq(-5, 5, length.out = N)
true_mean <- sin(t)
errors <- Reduce("rbind",
  lapply(1:n_blocks, function(i) {
    sigma <- diag(block_size) + rho_true *
      (matrix(1, block_size, block_size) - diag(block_size))
    matsqrt(sigma) %*% rnorm(block_size)
  })
)
y <- true_mean + errors * 0.5

model_fit <- lgspline(t, y,
  K = 3,
  correlation_id = rep(1:n_blocks, each = block_size),
  correlation_structure = "exchangeable",
  include_warnings = FALSE
)

## Propagate correlation uncertainty across 50 draws
post <- generate_posterior(model_fit,
  draw_correlation = TRUE, num_draws = 50,
  include_warnings = FALSE
)

## Fixed correlation parameters for comparison
post_fixed <- generate_posterior(model_fit, num_draws = 50)

corr_draws <- unlist(post$post_draw_correlation_params)
rho_draws <- exp(-exp(corr_draws))
print(summary(rho_draws))

lgspline documentation built on May 8, 2026, 5:07 p.m.