Deep Multivariate Skew t-Parsimonious Mixture Density Network (MST-PMDN)

Alex J. Cannon (alex.cannon@ec.gc.ca)

MST.PMDN is a 'torch for R' implementation of a distributional regression model based on a deep Multivariate Skew t-Parsimonious Mixture Density Network (MST-PMDN). The MST-PMDN framework represents complicated joint output distributions as mixtures of MST ('sn') components. A volume (L)-shape (A)-orientation (D) (LAD) eigenvalue decomposition parameterization provides a tractable, interpretable, and parsimonious representation of the scale matrices of the MST components, while explicit modeling of skewness and heavy tails can represent asymmetric behavior and tail dependence observed in real-world data (e.g., compound events and extremes). Overall, this provides an MST likelihood-based deep generative model.

In an MST-PMDN model, parameters of a mixture of multivariate skew t distributions that describe a multivariate output are estimated by training a deep learning model with two multi-modal input branches, one for tabular inputs and the other for (optional) image inputs. The two branches are provided as user-defined 'torch' modules. Outputs from each are concatenated and passed through a dense fusion network unless a custom fusion module is supplied, which then leads to the MST-PMDN head. In the absence of both branches, the tabular inputs are fed directly into the dense network. The overall network architecture is shown here.

Following the approach used in model-based clustering ('mclust'), scale matrices in the MST-PMDN head are represented using an LAD eigen-decomposition parameterization. LAD attributes, the nu (or degrees of freedom) parameter (n), and the alpha (or skewness) parameter (s) can be forced to be "V"ariable or "E"qual between mixture components (plus "I"dentity for A and D). An "N" constraint on n selects the exact infinite nu Gaussian-limit kernel; it is Gaussian when s is also "N" and skew-normal otherwise. An "N" constraint on s sets alpha to zero. Different model types are specified by setting constraint = "EIINN", "VEVEV", etc., where each letter position corresponds to an LADns attribute. With an "F" constraint on n, fixed_nu entries can be positive finite values, Inf for an exact Gaussian/skew-normal component, or NA for a component whose nu is learned over range_nu (default c(3, 50)). Fixed values are not clipped to range_nu; the range applies only to learned degrees of freedom. Furthermore, values of mu (the skew-t location parameter) (m), pi (or mixing coefficients) (x), volume-shape-orientation attributes (LAD), nu (n), and skewness (s) for the mixtures can be made independent of inputs by specifying any combination of constant_attr = "m", "mx", ..., "LADmxns". When skewness is nonzero, mu is not the component mean.

By combining appropriate values of constraint and constant_attr, MST-PMDN implements the Gaussian finite mixture models provided by 'mclust', i.e., for unconditional density estimation or model-based clustering:

mclust model	Description	MST-PMDN constraint =	MST-PMDN constant_attr =
EII	spherical, equal volume	"EIINN"	"LADmx"
VII	spherical, unequal volume	"VIINN"	"LADmx"
EEI	diagonal, equal volume and shape	"EEINN"	"LADmx"
VEI	diagonal, varying volume, equal shape	"VEINN"	"LADmx"
EVI	diagonal, equal volume, varying shape	"EVINN"	"LADmx"
VVI	diagonal, varying volume and shape	"VVINN"	"LADmx"
EEE	ellipsoidal, equal volume, shape, and orientation	"EEENN"	"LADmx"
EEV	ellipsoidal, equal volume and equal shape	"EEVNN"	"LADmx"
EVE	ellipsoidal, equal volume and orientation	"EVENN"	"LADmx"
VEE	ellipsoidal, equal shape and orientation	"VEENN"	"LADmx"
VEV	ellipsoidal, equal shape	"VEVNN"	"LADmx"
VVE	ellipsoidal, equal orientation	"VVENN"	"LADmx"
EVV	ellipsoidal, equal volume	"EVVNN"	"LADmx"
VVV	ellipsoidal, varying volume, shape, and orientation	"VVVNN"	"LADmx"

A comparison between 'mclust' and MST-PMDN with the constraints in the table above is shown here for the 'iris' dataset. Similarly, if the constraint on the nu parameter (n) is loosened (e.g., constraint = "VVVEN" with constant_attr = "LADmxn"), MST-PMDN can emulate model-based multivariate t clustering models provided by 'teigen'. Going one step further, removing the constraint on the skewness parameter (s) (e.g., constraint = "VVVEE" with constant_attr = "LADmxns") implements a restricted form of model-based multivariate skew t clustering ('EMMIXcskew' and previous packages).

While it can be used for model-based density estimation and clustering tasks, the primary purpose of the MST.PMDN package is to implement likelihood-based deep generative models. With unconstrained or partially constrained constant_attr, the MST-PMDN framework allows parameters of the mixture of multivariate Gaussian, t, or skew t distributions to depend on tabular and image covariates via user-specified torch modules. The interpretation layer evaluates scientifically meaningful distribution functionals, estimates tabular accumulated local effects and centred ICE, maps spatial image-occlusion effects, decomposes one-component contrasts among physical parameter channels, and attributes mixture exceedance risk to components. An example of the modelling use case, here demonstrated through simultaneous prediction of significant wave height and storm surge, is provided below.

Deep MST-PMDN Architecture

Installation

remotes::install_github("aljaca/MST.PMDN")

Example

library(MST.PMDN)

device <- ifelse(cuda_is_available(), "cuda", "cpu")
set.seed(1)
torch_manual_seed(1)

# Daily maximum significant wave height, storm surge, and covariates
# from Roberts Bank. Rows of date, x, x_image, and y are aligned. The lagged
# responses in x, all image channels, and y are already model-ready as
# documented in ?wave_surge.
date <- as.Date(wave_surge$date)
x <- wave_surge$x
x_image <- wave_surge$x_image  # 896 × 3 × 32 × 32: psl, uas, vas
y <- wave_surge$y

# The object also retains the image and response normalization statistics,
# the wave softplus scale, and both wave transformation functions.
x_image_mean <- wave_surge$x_image_mean
x_image_sd <- wave_surge$x_image_sd
y_mean <- wave_surge$y_mean
y_sd <- wave_surge$y_sd

# The TabularModule takes an input vector of length input_dim, runs it
# through two dense layers (input_dim→32 and 32→16) each with
# batch-norm (BN), ReLU and 50 %  dropout, then applies a final 16→16
# linear layer plus ReLU to produce a 16-dimensional output.
tabular_module <- nn_module(
  "TabularModule",
  initialize = function(
    input_dim,
    hidden_dims,
    output_dim,
    dropout_rate
  ) {
    # Number of hidden layers
    if (is.null(hidden_dims) || length(hidden_dims) == 0) {
      # No hidden layers
      self$n_hidden_layers <- 0
      self$hidden_dims <- c()
    } else if (!is.vector(hidden_dims) && !is.list(hidden_dims)) {
      # Single hidden size passed, wrap into vector
      self$hidden_dims <- c(hidden_dims)
      self$n_hidden_layers <- length(self$hidden_dims)
    } else {
      # Vector or list of hidden sizes
      self$hidden_dims <- hidden_dims
      self$n_hidden_layers <- length(self$hidden_dims)
    }
    # Store output size and dropout rate
    self$output_dim <- output_dim
    self$dropout_rate <- dropout_rate
    # Module lists for linear layers, batch-norms, dropouts
    self$layers <- nn_module_list()
    self$bns <- nn_module_list()
    if (self$dropout_rate > 0) {
      self$dropouts <- nn_module_list()
    }
    # Build hidden layers
    current_dim <- input_dim
    if (self$n_hidden_layers > 0) {
      for (i in seq_len(self$n_hidden_layers)) {
        # Linear transform
        self$layers$append(
          nn_linear(current_dim, self$hidden_dims[[i]])
        )
        # Batch normalization on hidden size
        self$bns$append(
          nn_batch_norm1d(self$hidden_dims[[i]])
        )
        # Optional dropout after activation
        if (self$dropout_rate > 0) {
          self$dropouts$append(
            nn_dropout(p = self$dropout_rate)
          )
        }
        # Update input size for next layer
        current_dim <- self$hidden_dims[[i]]
      }
    }
    # Final linear layer: last hidden (or input) → output_dim
    self$layers$append(
      nn_linear(current_dim, output_dim)
    )
  },
  forward = function(x) {
    # Pass through each hidden layer
    for (i in seq_len(self$n_hidden_layers)) {
      x <- self$layers[[i]](x)  # linear
      x <- self$bns[[i]](x)     # batch-norm
      x <- nnf_relu(x)          # activation
      # Apply dropout if configured
      if (self$dropout_rate > 0 && !is.null(self$dropouts[[i]])) {
        x <- self$dropouts[[i]](x)
      }
    }
    # Final projection and activation
    x <- self$layers[[length(self$layers)]](x)
    x <- nnf_relu(x)
    x
  }
)
tabular_mod <- tabular_module(
  input_dim = ncol(x),
  hidden_dims = c(32, 16),
  output_dim = 16,
  dropout_rate = 0.5
)

# The ImageModule accepts a 3×32×32 image, applies a 3×3 conv (3→16)
# with BN, ReLU  and 2×2 max-pool (→16×16), repeats with a 16→32 conv
# + BN, ReLU and max-pool (→8×8), flattens the 32×8×8 tensor to 2048
# units, and then projects it to 32 features via a linear layer, BN,
# and ReLU. Weight penalty (wd_image) is applied during training.
image_module <- nn_module(
  "ImageModule",
  initialize = function(
    in_channels,
    img_size,
    conv_channels,
    kernel_size = 3,
    pool_kernel = 2,
    output_dim = 32
  ) {
    # Store output dim
    self$output_dim <- output_dim
    # Build conv stack
    self$n_conv <- length(conv_channels)
    self$convs <- nn_module_list()
    self$bn_conv <- nn_module_list()
    # Track spatial dim through conv+pool
    spatial <- img_size
    pad <- floor(kernel_size / 2)
    for (i in seq_along(conv_channels)) {
      in_ch <- if (i == 1) in_channels else conv_channels[i-1]
      out_ch <- conv_channels[i]
      # conv keeps spatial size (with padding)
      self$convs$append(
        nn_conv2d(
          in_channels = in_ch,
          out_channels = out_ch,
          kernel_size = kernel_size,
          padding = pad
        )
      )
      self$bn_conv$append(nn_batch_norm2d(out_ch))
      # Pooling halves spatial dims
      spatial <- floor(spatial / pool_kernel)
    }
    # Store pooling layer and computed flatten_dim
    self$pool <- nn_max_pool2d(kernel_size = pool_kernel)
    self$flatten_dim <- tail(conv_channels, 1) * spatial * spatial
    # Final head: linear( flatten_dim → output_dim ) + BN
    self$fc    <- nn_linear(self$flatten_dim, output_dim)
    self$bn_fc <- nn_batch_norm1d(output_dim)
  },
  forward = function(x) {
    # conv → BN → ReLU → pool
    for (i in seq_len(self$n_conv)) {
      x <- self$convs[[i]](x)
      x <- self$bn_conv[[i]](x)
      x <- nnf_relu(x)
      x <- self$pool(x)
    }
    # Flatten and head
    x <- torch_flatten(x, start_dim = 2)
    x <- self$fc(x)
    nnf_relu(self$bn_fc(x))
  }
)
image_mod <- image_module(
  in_channels = dim(x_image)[2],
  img_size = dim(x_image)[3],
  conv_channels = c(16, 32),
  kernel_size = 3,
  pool_kernel = 2,
  output_dim = 32
)

# Define the fusion network, MST-PMDN head, and training setup
# Note: hyperparameters and number of epochs are not optimized
hidden_dim <- c(64, 32) # Layer widths in the default dense MLP
drop_hidden <- 0.1      # Dropout between non-final MLP layers
n_mixtures <- 2         # 2 components in the MST mixture model
constraint <- "VVIFN"   # LAD = "V"ariable-"V"ariable-"I"dentity; nu = 1 component "F"ixed; skewness = "N"ormal
fixed_nu <- c(Inf, NA)  # exact Gaussian 1st component; learned nu for 2nd
constant_attr <- ""     # All component attributes are free to vary with covariates
wd_tabular <- 0         # Weight decay for tabular module
wd_image <- 0.01        # Weight decay for image module
epochs <- 20            # Number of training epochs
lr <- 1e-3              # Initial Adam learning rate
batch_size <- 32        # Batch size

# Model training
fit <- train_mst_pmdn(
  inputs = x,
  outputs = y,
  hidden_dim = hidden_dim,
  drop_hidden = drop_hidden,
  n_mixtures = n_mixtures,
  constraint = constraint,
  fixed_nu = fixed_nu,
  constant_attr = constant_attr,
  epochs = epochs,
  lr = lr,
  batch_size = batch_size,
  wd_tabular = wd_tabular,
  wd_image = wd_image,
  image_inputs = x_image,
  image_module = image_mod,
  tabular_module = tabular_mod,
  checkpoint_path = "wave_surge_checkpoint.pt",
  device = device
)

# Model inference
pred <- predict_mst_pmdn(
  fit$model,
  new_inputs = x,
  image_inputs = x_image,
  device = device
)
print(names(pred))
print(pred$pi[1:3, ])
print(pred$mu[1:3, , ])
print(pred$nu[1:3, ])

# Draw samples
samples <- sample_mst_pmdn(
  pred,
  num_samples = 1000,
  device = device
)
print(str(samples))

# Evaluate CDF and quantiles
tt <- cdf_marginal_mst_pmdn(pred, y, draws = samples)
print(head(tt))

qq <- quantile_marginal_mst_pmdn(pred, tt, draws = samples)
print(head(qq))

# Functional-first interpretation of the complete two-component mixture
wave_q99 <- mst_functional("quantile", responses = 1L, prob = 0.99)
wave_q99_values <- functional_mst_pmdn(
  pred,
  wave_q99,
  num_samples = 4096L,
  seed = 20260806,
  device = device
)
print(head(wave_q99_values$data))

# ALE remains valid for a mixture. The matching image is held fixed whenever
# the selected tabular predictor changes.
wave_ale <- ale_mst_pmdn(
  fit$model,
  inputs = x,
  image_inputs = x_image,
  feature = 1L,
  functional = wave_q99,
  n_bins = 15L,
  num_samples = 4096L,
  seed = 20260806,
  device = device
)
plot(wave_ale)

# Mixture exceedance accounting separates component prevalence, within-
# component tail propensity, and contribution to total exceedance probability.
wave_tail_sources <- tail_components_mst_pmdn(
  pred,
  response = 1L,
  threshold = 2,
  num_samples = 4096L,
  seed = 20260806,
  device = device
)
print(head(wave_tail_sources$data))

Function Summaries

The deep MST-PMDN implementation consists of the following key functions and modules:

Function: t_cdf(z, nu)

• Purpose: Calculates a differentiable approximation of the univariate Student's t cumulative distribution function (CDF).
• Method: Uses Hill's transformed-normal approximation for finite nu >= 3, with exact closed-form CDFs for the Cauchy (nu = 1) and nu = 2 cases and the exact normal CDF for nu = Inf. A zero-safe factorization preserves gradients at the origin, while a direct, stable log-CDF path avoids lower-tail cancellation and probability flooring. The implementation is fully torch-compatible for use in autograd graphs.
• Context: Used within the loss function's skewness calculation to provide a fast, differentiable Student t log-CDF without switching between multiple implementations.

Function: sample_gamma(shape, scale, device)

• Purpose: Generates random samples from a Gamma distribution using torch.
• Method: Uses one vectorized R rgamma call and converts the output to a torch tensor with the requested dtype and device.
• Context: Used within the sample_mst_pmdn function to generate the scaling variable needed for sampling from the t-distribution component of the skew-t.

Function: build_orthogonal_matrix(params, dim)

• Purpose: Constructs a batch of orthogonal matrices (representing rotation/orientation D).
• Method: Uses the matrix exponential of a skew-symmetric matrix, where the input params parameterize the upper triangle of the skew-symmetric matrix.
• Context: Used in the main model (define_mst_pmdn) to generate the orientation component D of the LAD decomposition when orientation is not fixed to the identity matrix.

Function: init_mu_kmeans(model, outputs_train, ...)

• Purpose: Initializes the component location parameters (mu) using k-means clustering.
• Method: Applies k-means to the training output data to find initial centroids. These centroids initialize either the model$mu parameters (if constant) or the bias of the model$fc_mu layer (if network-dependent), setting initial weights to zero.
• Context: A heuristic to provide a potentially better starting point for training compared to random initialization, aiming for faster convergence.

Module: weight_norm_linear (nn_module)

• Purpose: Implements a linear layer with weight normalization.
• Method: Decomposes the weight matrix W into a direction V and a magnitude g, learning these instead of W directly.
• Context: Used for the MST-PMDN parameter-prediction heads. The default hidden MLP uses nn_linear layers; each non-final layer is followed by batch normalization, activation, and dropout.

Function: init_weight_norm(module)

• Purpose: Initializes the parameters (V, g) of a weight_norm_linear layer.
• Method: Uses Kaiming (He) normal initialization for the direction V and sets the initial magnitude g accordingly.
• Context: Applied recursively to the model to ensure proper initialization of all weight-normalized layers.

Module: define_mst_pmdn(...) (nn_module)

• Purpose: Defines the main MST-PMDN neural network architecture.
• Method:
  • Processes optional image and tabular inputs through dedicated modules or uses raw inputs.
  • When a fusion module is supplied, its output receives the optional drop_hidden dropout and is passed directly to the parameter-prediction heads. Otherwise, extracted features are concatenated and processed by a default hidden MLP constructed from nn_linear layers. In the custom-fusion case, hidden_dim creates no hidden layers and its final value is used only as a fallback for determining the fusion output dimension.
  • Predicts mixture parameters (pi, mu, L, A, D, nu, alpha) using separate output heads (mostly weight_norm_linear or nn_parameter if constant).
  • Applies Variable, Equal, Identity, exact Normal-limit, and Fixed constraints. Fixed nu = Inf marks an exact Gaussian/skew-normal component, while NA marks a learned finite-t component.
  • Constructs the full scale matrix Sigma = L * D * diag(A) * D^T and computes its Cholesky decomposition (scale_chol) for each component.
• Output: Returns a list containing all mixture parameters (pi, mu, scale_chol, nu, alpha) and LAD components (L, A, D), batched appropriately.

Function: loss_mst_pmdn(output, target)

• Purpose: Computes the negative log-likelihood (NLL) loss.
• Method:
  • For each data point and mixture component k:
    • Calculates residuals: diff = target - mu_k.
    • Standardizes residuals: v = scale_chol_k^{-1} * diff.
    • Calculates squared Mahalanobis distance: maha = ||v||^2.
    • For finite nu, calculates the symmetric multivariate t log-PDF using a cancellation-resistant normalizing constant and the skewness adjustment log(2 * T_CDF(alpha_k^T w, df=nu_k+d)).
    • For nu = Inf, evaluates the exact Gaussian or skew-normal limit using the multivariate normal log-PDF and log(2 * Phi(alpha_k^T v)).
  • Combines component log-densities using mixture weights pi via logsumexp.
  • Returns the mean NLL over the batch.
  • Optionally adds an L2 penalty on the final alpha values via lambda_alpha, and on (1/nu)^2 via lambda_nu_inv.

Function: sample_mst_pmdn(mdn_output, num_samples, ...)

• Purpose: On-device generation of random samples from the predicted mixture distribution.
• Method:
  • Samples component indices based on pi.
  • Gathers parameters for the selected components.
  • Generates t-distribution scaling factors W using sample_gamma for finite-t components and sets W = 1 exactly for Gaussian/skew-normal components.
  • Generates a standard multivariate skew-normal sample X based on the component's alpha (via delta).
  • Transforms the standard sample X to the output space: Y = mu_s + W * (scale_chol_s @ X).
• Output: Returns a list with
  • samples - a torch tensor of shape [S, B, d], where S is num_samples, B is the batch size (rows of the predictor matrix), and d is the response dimension.
  • components - a torch tensor of shape [S, B] giving the 1-based component label (1..G) used for each draw.

Function: sample_mst_pmdn_df(mdn_output, num_samples, ...)

• Purpose: Generates random samples from the predicted mixture distribution and returns a formatted R data frame.
• Method:
  • Samples component indices based on pi.
  • Gathers parameters for the selected components.
  • Generates t-distribution scaling factors W using sample_gamma for finite-t components and sets W = 1 exactly for Gaussian/skew-normal components.
  • Generates a standard multivariate skew-normal sample X based on the component's alpha (via delta).
  • Transforms the standard sample X to the output space: Y = mu_s + W * (scale_chol_s @ X).
• Output: A data frame with num_samples * batch_size rows containing
  • simulated response variables in columns V1 ... Vd;
  • row - the index (1..B) of the predictor row that generated the draw;
  • draw - the draw number (1..num_samples) for that predictor row;
  • comp - a factor giving the 1-based component label (1..G).

Function: cdf_marginal_mst_pmdn(mdn_output, y, var_index = NULL, ...)

• Purpose: Estimates the marginal CDF for one or more response dimensions.
• Method: Uses Monte Carlo sampling from the mixture to approximate (F(y_j) = sum_{k=1}^G pi_k F_{k,j}(y_j)), because the component skew-t marginal CDFs have no closed form in general.
• Context: Use when you need marginal probabilities for a fitted mixture; see quantile_marginal_mst_pmdn() for inverse-CDF summaries.

Function: quantile_marginal_mst_pmdn(mdn_output, probs, var_index = NULL, ...)

• Purpose: Estimates marginal quantiles for one of more response dimensions.
• Method: Uses Monte Carlo sampling to invert the mixture marginal CDF; the mixture quantile is not the weighted sum of component quantiles.
• Context: Use when summarizing predictive distributions; pairs naturally with cdf_marginal_mst_pmdn() for probability and quantile summaries.

Function: train_mst_pmdn(...)

• Purpose: Manages the model training process.
• Method: Includes data loading, model/optimizer setup (with k-means init), training loop (loss calculation, backpropagation, optimization), complete-case validation, learning rate scheduling, checkpointing, and early stopping. The latest resumable state is stored at checkpoint_path; the best model is stored separately at a derived _best path. Resumption restores the saved split, optimizer, histories, counters, and available R/torch RNG states. Handles optional image inputs and allows weighting an L2 penalty on alpha through lambda_alpha and on (1/nu)^2 through lambda_nu_inv.
• Output: Trained model, loss history, and training/validation indices.

Function: predict_mst_pmdn(model, new_inputs, ...)

• Purpose: Performs inference using the trained model.
• Method: Runs a forward pass on new inputs in evaluation mode (torch_no_grad()).
• Output: Raw model output list containing mixture parameters for the new inputs.

Function: scov_mst_pmdn(pred, type = c("scale", "scale_chol", "cov"), ...)

• Purpose: Returns component scale matrices, their Cholesky factors, or actual skew-t/skew-normal covariance matrices.
• Method: Uses the scale_chol factor employed by the likelihood. For type = "cov", the scale is transformed using both nu and alpha; covariance is undefined for finite nu <= 2, while nu = Inf uses the exact skew-normal limit.
• Output: A 4D tensor (or R array if as_array = TRUE) of shape [batch_size, M, d, d].

Distribution-functional interpretation

Functions: mst_functional(...) and functional_mst_pmdn(...)

• Purpose: Define and evaluate one scalar scientific summary per prediction row.
• Method: Means, variances, standard deviations, covariance, and correlation use exact component and mixture moments, including between-component covariance. Quantiles, marginal and joint exceedances, tail spread, and normalized generalized Bowley tail asymmetry use a parameter-independent latent bank created by latent_draws_mst_pmdn().
• Diagnostics: Tail summaries report the expected number of Monte Carlo draws in the relevant tail and flag inadequate resolution. Every public interpretation call aggregates its internal evaluations and emits at most one classed tail-resolution warning; the per-evaluation distribution remains available in diagnostics. Mixture results separately retain expected component draw counts; because component uniforms are shared across rows, component-selection Monte Carlo error does not average away along an effect curve. Exact moments return undefined values when finite degrees of freedom do not support them.

Functions: ale_mst_pmdn(...) and ice_mst_pmdn(...)

• Purpose: Explain how a tabular covariate changes a selected distribution functional.
• Method: One-dimensional ALE uses non-empty empirical bins and locally supported lower-to-upper contrasts. Effects are displayed at bin midpoints by subtracting half the current bin effect, which assumes within-bin linearity rather than using the boundary-valued Apley-Zhu convention. Centred ICE shows case-level heterogeneity with an optional or precomputed ALE overlay. A case's image remains aligned and fixed during tabular perturbations.
• Scope: A named feature is accepted when input columns are named; otherwise use R's 1-based column indices. Deterministically linked predictors, such as sine/cosine seasonal harmonics, should be perturbed as a scientifically coherent group outside this one-dimensional API.

Functions: image_contrast_mst_pmdn(...) and image_occlusion_mst_pmdn(...)

• Purpose: Measure whole-image and spatial patch effects on the same scalar functionals used by the tabular layer.
• Method: Observed image regions are replaced by a required reference field, optionally with cosine tapering and overlapping patches. Named channel_groups can attribute fields separately; NULL retains a joint synoptic-state perturbation. When an input channel is deterministically derived from another, a rebuild_channels callback should perturb the fundamental physical field, recompute every linked channel such as a pressure gradient, apply the original preprocessing, and return a complete model-ready tensor. Callback masks have one channel for a joint perturbation and the full model channel count under grouping.
• Monte Carlo stability: image_occlusion_mst_pmdn() accepts multiple independent latent banks, computes every CNN prediction once, and repeats only functional sampling. It returns bank-specific effects, their spread, and a strict same-sign indicator. The indicator is a sign-stability heuristic rather than a confidence interval.
• Interpretation: Patch effects are spatial sensitivity measures. Cosine taper weights are reported as metadata but do not define a linear normalization through the nonlinear network and functional. Raw effects are retained, and they do not generally add across overlapping patches to the whole-image contrast.

Function: decompose_mst_pmdn(...)

• Purpose: For one-component models, allocate a known functional contrast among location, complete scale, skewness, and degrees-of-freedom channels.
• Method: Exact Shapley averaging evaluates all hybrid parameter states—16 when all four channels are active—and reports the numerical sum-to-total residual. Structurally inactive or exactly unchanged channels disappear automatically; structural symmetry and an exactly zero skew vector are treated as the same endpoint. Diagnostics report the maximum location, Cholesky-factor, skew-vector, and inverse-df change. The complete Cholesky factor is the scale block; because it maps standardized skew direction into response space, its contribution includes the induced rotation of that direction.
• Scope: This is parameter-channel attribution rather than feature SHAP. Full channel decomposition is disabled for mixtures because component labels and compensating component changes are not identified.

Function: tail_components_mst_pmdn(...)

• Purpose: Provide a mixture-safe explanation of exceedance probability.
• Output: For each component, reports its weight, within-component exceedance probability, contribution to total exceedance probability, tail share, and contribution rank. Probability is estimated directly within each component and combined analytically with its mixture weight, avoiding extra mixture-selection noise. Component labels remain descriptive indices rather than assumed physical regimes.

References

Ambrogioni, L., Güçlü, U., van Gerven, M. A., & Maris, E. (2017). The kernel mixture network: A nonparametric method for conditional density estimation of continuous random variables. arXiv:1705.07111.

Andrews, J. L., & McNicholas, P. D. (2012). Model-based clustering, classification, and discriminant analysis via mixtures of multivariate t-distributions: the t EIGEN family. Statistics and Computing, 22, 1021-1029.

Azzalini, A., & Capitanio, A. (2003). Distributions generated by perturbation of symmetry with emphasis on a multivariate skew t-distribution. Journal of the Royal Statistical Society Series B: Statistical Methodology, 65(2), 367-389.

Andrews, J. L., Wickins, J. R., Boers, N. M., & McNicholas, P. D. (2018). teigen: An R package for model-based clustering and classification via the multivariate t distribution. Journal of Statistical Software, 83, 1-32.

Banfield, J. D., & Raftery, A. E. (1993). Model-based Gaussian and non-Gaussian clustering. Biometrics, 803-821.

Celeux, G., & Govaert, G. (1995). Gaussian parsimonious clustering models. Pattern Recognition, 28(5), 781-793.

Falbel D., & Luraschi, J. (2026). torch: Tensors and Neural Networks with 'GPU' Acceleration. R package version 0.17.0, doi:10.32614/CRAN.package.torch

Fraley, C., & Raftery, A. E. (2002). Model-based clustering, discriminant analysis, and density estimation. Journal of the American Statistical Association, 97(458), 611-631.

Fraley, C., & Raftery, A. E. (1998). How many clusters? Which clustering method? Answers via model-based cluster analysis. The Computer Journal, 41(8), 578-588.

Lee, S., & McLachlan, G. J. (2018). EMMIXcskew: An R Package for the Fitting of a Mixture of Canonical Fundamental Skew t-Distributions. Journal of Statistical Software, 83(3), 1-32.

Kingma, D. P., & Ba, J. (2015). Adam: a method for stochastic optimization. Proceedings of the 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA. arXiv:1412.6980

Klein, N. (2024). Distributional regression for data analysis. Annual Review of Statistics and Its Application, 11:321-346.

Peel, D., & McLachlan, G.J. (2000). Robust mixture modelling using the t distribution. Statistics and Computing 10, 339–348.

Srucca, L., Fop, M., Murphy, T. B., & Raftery, A. E. (2016). mclust 5: Clustering, classification and density estimation using Gaussian finite mixture models. The R Journal, 8(1), 289-317.

Williams, P. M. (1996). Using neural networks to model conditional multivariate densities. Neural Computation, 8(4), 843-854.