Package 'proxymix'

Description

Packages a gmm_imputation as a mice::mids object so that an arbitrary model estimand can be fitted and pooled with the established mice workflow. The joint Gaussian-mixture imputations – including the multimodal and heteroscedastic shapes a univariate imputer cannot produce – flow through unchanged; mice supplies with(), mice::pool(), and the pooled diagnostics.

Usage

as_mids(object)

Arguments

Value

A mice::mids object with m imputations.

See Also

gmm_impute(), proxy_pool() for the closed-form column-mean pooling.

Other imputation: gmm_complete(), gmm_imputation(), gmm_impute(), mechanism, proxy_fmi(), proxy_mnar_sensitivity(), proxy_pool()

Examples

set.seed(1)
x1 <- rnorm(150); x2 <- x1 + rnorm(150)
x2[runif(150) < 0.3] <- NA
imp <- gmm_impute(cbind(x1, x2), N = 1L, m = 10L, seed = 1L)
if (requireNamespace("mice", quietly = TRUE)) {
  fit <- with(as_mids(imp), lm(x2 ~ x1))
  summary(mice::pool(fit))
}

Description

An autoplot() method for gmm_fit objects, rendering the fitted mixture with ggplot2. The displayed coordinates are reduced to the requested one or two dimensions through the closed-form marginal gmm_marginalise(), so the method works for a proxy of any ambient dimension p.

A one-dimensional request draws the marginal mixture density, optionally with the per-component densities underneath and a rug of the target's samples. A two-dimensional request draws the marginal density as a viridis raster with white contour lines, optionally overlaying each component's mean and a probability-contour ellipse.

Arguments

Details

The method is registered against the ggplot2::autoplot() generic only when ggplot2 is installed; call it as ggplot2::autoplot(fit) or load ggplot2 first. It returns the ggplot object, so the usual + layering applies for further customisation.

Value

A ggplot object.

See Also

Other classes: glance.gmm_fit, gmm(), gmm_counterfactual_law(), gmm_dim(), gmm_fit(), gmm_n_components(), gmm_target(), gmm_weights(), is_proposal(), tidy.gmm

Examples

samples <- matrix(stats::rnorm(200), ncol = 2)
tgt <- gmm_target_from_samples(samples)
fit <- fit_proxymix(tgt, N = 2L, regime = "sample", max_iter = 25L)
ggplot2::autoplot(fit)
ggplot2::autoplot(fit, dims = 1L)

Description

A 2-D "banana" density obtained by warping an isotropic Gaussian through the map $(z_1, z_2) \mapsto (z_1, z_2 + \tfrac{1}{2}(z_1^2 - 1))$. The map has unit Jacobian, so the resulting density is exactly normalised.

Usage

banana_target(with_samples = FALSE, n = 2000L, seed = 1L)

Arguments

Value

A gmm_target in dimension 2.

See Also

Other targets: donut_target(), epanechnikov_target(), gmm_target_from_samples(), maxent_target(), mixture_target()

Examples

b <- banana_target()
b
b@log_density(matrix(c(0, 0, 1, 0), ncol = 2, byrow = TRUE))

Description

Returns the Bayesian and Akaike information criteria of a regime-(ii) fit, together with the integrated completed likelihood (ICL). All three are computed against the empirical log-likelihood of the samples used to fit the model and are reported on the same scale (smaller is better). They are NA for regimes that do not have an empirical likelihood ("moment", "kld").

Usage

bic_aic(fit)

Arguments

Details

The ICL of Biernacki, Celeux and Govaert (2000) adds to the BIC twice the entropy of the fitted classification, $\mathrm{ICL} = \mathrm{BIC} + 2 E_N$, where $E_N = -\sum_{i,k} \gamma_{ik} \log \gamma_{ik} \ge 0$ is the entropy of the responsibilities $\gamma_{ik}$. It therefore penalises mixtures whose components overlap (uncertain assignments), and favours well-separated clustering solutions over the merely best-fitting ones. Because $E_N \ge 0$, the ICL is never smaller than the BIC, and the two coincide for a single component ($K = 1$), where every responsibility is one. The classification entropy itself is returned as classification_entropy.

Value

A list with bic, aic, icl, classification_entropy, and n_params.

References

Biernacki, C., Celeux, G. and Govaert, G. (2000) Assessing a mixture model for clustering with the integrated completed likelihood. IEEE Transactions on Pattern Analysis and Machine Intelligence 22(7), 719–725. doi:10.1109/34.865189

See Also

Other diagnostics: ess_summary(), ess_trace(), gmm_anneal_path(), gmm_conditional_entropy(), gmm_entropy(), gmm_evidence(), gmm_fit_ensemble(), gmm_fit_quality(), gmm_independence_graph(), gmm_mutual_information(), hellinger_mc(), kld_trace(), proxy_functional_ci()

Examples

x <- matrix(stats::rnorm(200), ncol = 2)
tgt <- gmm_target_from_samples(x)
fit <- fit_proxymix(tgt, N = 2L, regime = "sample", max_iter = 25L)
bic_aic(fit)

Description

Evaluates the density (or log-density) of a Gaussian mixture at one or more points.

Usage

dgmm(x, g, log = FALSE)

Arguments

Value

A numeric vector of length nrow(x).

See Also

Other ops: gmm_canonicalise(), gmm_conditionalise(), gmm_divergence(), gmm_kld(), gmm_marginalise(), gmm_mean(), gmm_modes(), pgmm(), rgmm()

Examples

g <- gmm(weights = c(0.5, 0.5),
         means = list(c(-1, 0), c(1, 0)),
         covariances = list(diag(2), diag(2)))
dgmm(c(0, 0), g)
dgmm(c(0, 0), g, log = TRUE)

Description

A rotationally symmetric annulus on $\mathbb{R}^2$, with density

$f(x) \propto \exp\!\left(-\tfrac{(\Vert x \Vert - r_0)^2}{2 \sigma^2}\right).$

Numerical integration in polar coordinates fixes the normaliser; the returned target exposes a normalised log_density.

Usage

donut_target(r0 = 2.5, sigma = 0.5, with_samples = FALSE, n = 2000L, seed = 1L)

Arguments

Value

A gmm_target in dimension 2.

See Also

Other targets: banana_target(), epanechnikov_target(), gmm_target_from_samples(), maxent_target(), mixture_target()

Examples

d <- donut_target()
d

Description

The canonical compact-support target: a product of one-dimensional Epanechnikov densities $K(u) = \tfrac{3}{4}(1 - u^2)\,\mathbf{1}_{|u| \le 1}$, rescaled to ⁠[center - half_width, center + half_width]⁠ in each coordinate. No mixture of full-support Gaussians can have compact support, so this target is the clean case where regime (iii) is the only viable fitting route. It declares its support, which makes fit_kld_em() (and fit_proxymix() with regime = "kld") select a support-matched uniform proposal automatically instead of the default multivariate-t, which would place importance mass where the log-density is -Inf.

Usage

epanechnikov_target(
  n_dim = 1L,
  center = 0,
  half_width = 1,
  with_samples = FALSE,
  n = 2000L,
  seed = 1L
)

Arguments

Value

A gmm_target in dimension n_dim with a declared compact support.

See Also

Other targets: banana_target(), donut_target(), gmm_target_from_samples(), maxent_target(), mixture_target()

Examples

e <- epanechnikov_target()
e
e@log_density(matrix(c(0, 0.5, 1.5), ncol = 1L)) # finite, finite, -Inf

Description

Convenience accessor returning the headline IS-quality numbers for a regime-(iii) fit: effective sample size and its ratio to is_size, the largest self-normalised weight, the support fraction (proportion of draws that received a finite weight), and the Monte-Carlo standard error of the final KLD estimate. Returns NA fields for regimes that do not use importance sampling.

Usage

ess_summary(fit)

Arguments

Details

Validation-side numbers (⁠validation_*⁠) are populated only when the fit was called with validation_size > 0.

Value

A list of numeric scalars (or NAs where not applicable).

See Also

Other diagnostics: bic_aic(), ess_trace(), gmm_anneal_path(), gmm_conditional_entropy(), gmm_entropy(), gmm_evidence(), gmm_fit_ensemble(), gmm_fit_quality(), gmm_independence_graph(), gmm_mutual_information(), hellinger_mc(), kld_trace(), proxy_functional_ci()

Examples

fit <- fit_proxymix(banana_target(), N = 3L, regime = "kld",
                    is_size = 1500L, max_iter = 20L, seed = 1L,
                    validation_size = 1500L)
ess_summary(fit)

Description

Returns the effective sample size (1 / sum(W^2)) of the self-normalised importance weights used by a regime-(iii) fit. NA for regimes that do not use importance sampling.

Usage

ess_trace(fit)

Arguments

Value

Numeric scalar (or NA_real_).

See Also

Other diagnostics: bic_aic(), ess_summary(), gmm_anneal_path(), gmm_conditional_entropy(), gmm_entropy(), gmm_evidence(), gmm_fit_ensemble(), gmm_fit_quality(), gmm_independence_graph(), gmm_mutual_information(), hellinger_mc(), kld_trace(), proxy_functional_ci()

Examples

fit <- fit_proxymix(banana_target(), N = 2L, regime = "kld",
                    is_size = 1000L, max_iter = 15L, seed = 1L)
ess_trace(fit)

Description

Implements regime (ii) of Hoek and Elliott (2024). Runs the textbook expectation-maximisation algorithm for Gaussian mixtures on the supplied samples, with diagonal ridge regularisation for numerical stability, optional multi-start, and monotone-log-likelihood checking.

Usage

fit_em_samples(
  target,
  N = 2L,
  init = NULL,
  max_iter = 100L,
  tol = 1e-06,
  ridge_eps = 1e-06,
  n_starts = 5L,
  anneal = FALSE,
  temp_schedule = NULL,
  seed = NULL,
  canonicalise = TRUE
)

Arguments

Value

A gmm_fit with regime = "sample". When anneal = TRUE the diagnostics list also carries annealed = TRUE and the temp_schedule used.

See Also

Other fitting: fit_kld_em(), fit_moment_match(), from_kde(), from_objective(), select_N()

Examples

x <- matrix(stats::rnorm(200), ncol = 2)
tgt <- gmm_target_from_samples(x)
fit <- fit_em_samples(tgt, N = 2L, max_iter = 30L, n_starts = 2L)
fit@diagnostics$loglik_final

Description

Implements regime (iii) of Hoek and Elliott (2024). Minimises KL(f || g_theta) where f is supplied as an evaluable log-density on the target, via expectation-maximisation against importance-sampled draws from a user-chosen proposal q.

Usage

fit_kld_em(
  target,
  N = 3L,
  proposal = NULL,
  is_size = 5000L,
  init = NULL,
  max_iter = 100L,
  tol = 1e-05,
  ridge_eps = 1e-06,
  min_ess = 50,
  on_low_ess = c("warn", "abort"),
  seed = NULL,
  validation_size = NULL,
  validation_proposal = NULL,
  validation_seed = NULL,
  support_warn = TRUE,
  adapt = c("none", "pmc"),
  refresh_every = 5L,
  defensive_gamma = 0.15,
  inflate = 1.5,
  anneal = FALSE,
  temp_schedule = NULL,
  canonicalise = TRUE
)

Arguments

Details

With adapt = "none" (the default) the Monte Carlo draws from q are computed once at the start and the resulting self-normalised importance-sampling weights are reused at every EM iteration. With adapt = "pmc" the proposal is refreshed every refresh_every iterations with a defensive mixture built from the current iterate – the population-Monte-Carlo scheme: the fitted mixture (covariances inflated by inflate) carries 1 - defensive_gamma of the proposal mass and the original proposal q keeps defensive_gamma as a heavy-tailed anchor, a fresh IS batch is drawn, and EM continues on the refreshed weights. Because the refreshed proposal tracks the target, the effective sample size recovers from a poor initial proposal and the usable dimension range extends well beyond what a fixed proposal reaches; the per-batch ESS trace is reported as diagnostics$ess_history. While a batch is degenerate (its effective sample size is below min_ess), the refresh fires every iteration with an escalating covariance inflation floored at a growing fraction of the batch's sample covariance, so a collapsed iterate walks back out toward the target instead of freezing; and convergence is only accepted on an adapted batch, so a run that stabilises on the original proposal's draw is refreshed at least once before it is allowed to stop. The scheme is the mixture population-Monte-Carlo idea of Cappé et al. (2008) with the defensive-mixture safeguard of Owen and Zhou (2000); it re-draws rather than recycles batches (compare the adaptive multiple importance sampling of Cornuet et al., 2012).

Since v0.1.1 the function also draws an independent validation IS sample when validation_size > 0 and reports its own KLD estimate, effective sample size, and largest weight share. This lets users tell the difference between in-sample EM overfit to one particular IS draw and a fit that generalises across independent IS draws.

When the target's normalised property is FALSE or NA, the importance-sampled kld_final and kld_trace measure $\widehat{KL}(f \Vert g) - \log Z(f)$ rather than the absolute divergence. The fit's diagnostics list records this via kld_is_shifted = TRUE and a kld_shift_explanation string. When the target also supplies a finite log_normalizer, a corrected absolute estimate is reported as kld_final_absolute.

Value

A gmm_fit with regime = "kld". The diagnostics list contains, among others, kld_trace, kld_final, kld_is_shifted, kld_final_absolute (when computable), ess, ess_relative (ess / is_size), max_weight, support_fraction, mc_se_kld, validation_kld, validation_ess, and validation_max_weight.

References

Cappé, O., Douc, R., Guillin, A., Marin, J.-M. and Robert, C. P. (2008) Adaptive importance sampling in general mixture classes. Statistics and Computing 18, 447–459. doi:10.1007/s11222-008-9059-x

Cornuet, J.-M., Marin, J.-M., Mira, A. and Robert, C. P. (2012) Adaptive multiple importance sampling. Scandinavian Journal of Statistics 39, 798–812. doi:10.1111/j.1467-9469.2011.00756.x

Owen, A. and Zhou, Y. (2000) Safe and effective importance sampling. Journal of the American Statistical Association 95(449), 135–143. doi:10.1080/01621459.2000.10473909

See Also

Other fitting: fit_em_samples(), fit_moment_match(), from_kde(), from_objective(), select_N()

Examples

tgt <- banana_target()
q <- is_mvt(n_dim = 2L, mean = c(0, 0),
            sigma = 4 * diag(2), df = 5)
fit <- fit_kld_em(tgt, N = 3L, proposal = q,
                  is_size = 1500L, max_iter = 25L, seed = 1L,
                  validation_size = 1500L)
fit@diagnostics$kld_final
fit@diagnostics$validation_kld

Description

Implements regime (i) of Hoek and Elliott (2024). When N == 1, this is the exact moment match: mu is the target mean and Sigma is the target covariance. When N > 1, the function returns the deterministic moment-seed of init_moment_seed() wrapped as a gmm_fit, without iterative refinement — useful as a starting point for the iterative regimes.

Usage

fit_moment_match(target, N = 1L, ridge_eps = 1e-06, canonicalise = TRUE)

Arguments

Details

Either the target must carry an n by p samples matrix, or its metadata slot must contain pre-computed moments of the form ⁠list(mean = <p-vec>, cov = <p-by-p>)⁠.

Value

A gmm_fit with regime = "moment".

See Also

Other fitting: fit_em_samples(), fit_kld_em(), from_kde(), from_objective(), select_N()

Examples

x <- matrix(stats::rnorm(200), ncol = 2)
tgt <- gmm_target_from_samples(x)
fit_moment_match(tgt, N = 1L)

Description

The unified front door of proxymix. Picks a fitting regime (or honours an explicit choice) and dispatches to the corresponding regime-specific fitter:

Usage

fit_proxymix(
  target,
  N = 1L,
  regime = c("auto", "moment", "sample", "kld"),
  ...
)

Arguments

Details

• "moment" - closed-form moment matching (fit_moment_match()).

• "sample" - classical EM on i.i.d. samples (fit_em_samples()).

• "kld" - importance-sampled KLD-EM (fit_kld_em()) for a target that can be evaluated but not sampled.

With regime = "auto" the choice is made from the shape of the supplied target:

• N == 1 and the target carries samples or moments: "moment".

• N >= 2 and the target carries samples: "sample".

• The target carries log_density only (no samples): "kld".

Value

A gmm_fit.

Examples

## auto: samples + N=2 -> classical EM.
x <- matrix(stats::rnorm(200), ncol = 2)
tgt_s <- gmm_target_from_samples(x)
fit_proxymix(tgt_s, N = 2L, max_iter = 25L)

## explicit "kld" on a log-density-only target.
fit_proxymix(banana_target(), N = 3L, regime = "kld",
             is_size = 1000L, max_iter = 20L, seed = 1L)

Description

Assembles a joint Gaussian-mixture proxy over the outcome, the treatment, and the covariates, and returns an uplift_model that the decision verbs read in closed form. One fit yields prediction, heterogeneous treatment effects, optimal actions, off-line policy value, and an identification audit – see proxy_cate(), proxy_decide(), proxy_policy_value() and proxy_identification_report().

Usage

fit_uplift(
  data,
  outcome,
  treatment,
  covariates,
  N = "auto",
  regime = "auto",
  assume = c("ignorability", "latent_confounder"),
  outcome_type = c("continuous", "binary", "count"),
  n_grid = 1:4,
  seed = NULL,
  ...
)

Arguments

Details

The component count N may be fixed or chosen automatically. With N = "auto" the function sweeps n_grid and keeps the K that minimises the joint BIC; the full BIC trace is stored in the model's metadata. The treatment is binary at this version (⁠{t0, t1}⁠); a continuous dose is a future extension.

Value

An uplift_model.

See Also

proxy_cate(), proxy_decide(), proxy_identification_report()

Other decision: gmm_cf_mean(), gmm_cf_tail_prob(), gmm_cf_variance(), gmm_counterfactual(), gmm_intervene(), proxy_cate(), proxy_confounding_gap(), proxy_decide(), proxy_identification_report(), proxy_overlap(), proxy_policy_value(), proxy_predict(), proxy_regime_segments(), proxy_retrospective_uplift(), proxy_uplift(), uplift_identification(), uplift_model()

Examples

set.seed(1)
n <- 400L
x <- stats::rnorm(n)
t <- stats::rbinom(n, 1L, 0.5)
y <- 1 + 0.5 * x + (1 + x) * t + stats::rnorm(n, sd = 0.5)
dat <- data.frame(y = y, t = t, x = x)
m <- fit_uplift(dat, outcome = "y", treatment = "t", covariates = "x",
                N = 2L, regime = "sample", max_iter = 50L, seed = 1L)
m

Description

Fits an N-component Gaussian-mixture proxy to a (Gaussian, diagonal- bandwidth) kernel-density estimate over samples, via regime (iii) KLD-EM. The proxy is closed-form marginalisable, conditionable, and samplable; the KDE is none of those things on its own.

Usage

from_kde(
  samples,
  N = 3L,
  bandwidth = "silverman",
  proposal = NULL,
  is_size = 5000L,
  max_iter = 100L,
  tol = 1e-05,
  ridge_eps = 1e-06,
  min_ess = 50L,
  seed = NULL,
  validation_size = 0L,
  validation_proposal = NULL,
  validation_seed = NULL,
  support_warn = TRUE,
  canonicalise = TRUE
)

Arguments

Details

This is a compression operation: take an n-sample KDE and replace it with the closest N-component mixture in the Kullback-Leibler sense (which is much smaller than n for typical use). Bias inherited from the KDE is reproduced in the proxy; the bandwidth controls the bias-variance trade-off.

Dimensional scope. The dimensional guard is p <= 5 (recommended), p <= 10 (allowed with warning), p > 10 (rejected). Regime-(iii) KLD-EM is driven by importance sampling, whose effective sample size collapses sharply in high dimensions.

Value

A gmm_fit with regime = "kld" and metadata recording the KDE inputs (kde_samples_n, bandwidth, bandwidth_method).

See Also

Other fitting: fit_em_samples(), fit_kld_em(), fit_moment_match(), from_objective(), select_N()

Examples

set.seed(1L)
x <- rbind(
  mvnfast::rmvn(120L, mu = c(-2, 0), sigma = diag(2)),
  mvnfast::rmvn(120L, mu = c( 2, 0), sigma = diag(2))
)
fit <- from_kde(x, N = 2L, is_size = 2000L, max_iter = 40L, seed = 1L)
fit
ess_summary(fit)

Description

Fits a Gaussian-mixture proxy to the Gibbs measure $\exp(-f(x) / T)$ of a user-supplied objective f over a bounded box, by cooling a short temperature ladder through regime-(iii) importance-sampled KLD-EM (fit_kld_em()). As the temperature falls the mixture mass concentrates on the low regions of f, so the fitted mixture is a closed-form map over the optima rather than a single point estimate. Pair it with gmm_modes() to read off the distinct optima.

Usage

from_objective(
  objective,
  lower,
  upper,
  N = NULL,
  minimise = TRUE,
  temperature = NULL,
  n_steps = 6L,
  exploration = 0.5,
  inflate = 1.8,
  is_size = 10000L,
  max_iter = 70L,
  ridge_eps = 1e-04,
  seed = NULL
)

Arguments

Details

The Gibbs measure can be evaluated point-wise but not directly sampled, which is precisely the setting of regime (iii): minimising the Kullback-Leibler divergence from a Gaussian mixture to a peaked target is the rank-weighted Gaussian update at the heart of fit_kld_em(). This function is that fit, driven against a sequence of cooling Gibbs targets and warm-started from the previous fit at each step. Because a multimodal f produces a multimodal target, the components spread across the basins and recover the optima together.

Recovery is most reliable with component headroom – a number of components N comfortably larger than the number of optima you expect, so the defensive proposal can keep a component on each basin. Symmetric landscapes (where several optima are exchangeable) need the most headroom.

Dimensional scope. As with the rest of regime (iii), the importance- sampling effective sample size falls sharply with dimension; the guard is p <= 5 (recommended), p <= 10 (allowed with a warning), p > 10 (rejected).

Value

A gmm_fit (the fitted proxy) carrying a from_objective metadata record with the temperature ladder and box. Pass it to gmm_modes() to extract the distinct optima.

See Also

gmm_modes() to resolve the fitted map into distinct optima.

Other fitting: fit_em_samples(), fit_kld_em(), fit_moment_match(), from_kde(), select_N()

Examples

## A bimodal 1-D objective with minima at +/- 2.
f <- function(v) (v[1]^2 - 4)^2
fit <- from_objective(f, lower = -5, upper = 5, N = 6L,
                      is_size = 2000L, n_steps = 5L, seed = 1L)
gmm_modes(fit)$modes

Description

A broom-style glance() method: a one-row summary of a gmm_fit with the regime, the component count and dimension, convergence, iteration count, and the regime's headline fit statistics. Available as generics::glance(fit) when the generics package is installed.

Arguments

Value

A one-row data frame.

See Also

Other classes: autoplot.gmm_fit, gmm(), gmm_counterfactual_law(), gmm_dim(), gmm_fit(), gmm_n_components(), gmm_target(), gmm_weights(), is_proposal(), tidy.gmm

Examples

fit <- fit_proxymix(banana_target(), N = 2L, regime = "kld",
                    is_size = 1000L, max_iter = 10L, seed = 1L)
generics::glance(fit)

Description

Lightweight S7 class representing an N-component multivariate Gaussian mixture on $\mathbb{R}^p$. Use gmm() to construct, dgmm() / rgmm() to evaluate or sample, and gmm_marginalise() / gmm_conditionalise() for closed-form operations.

Usage

gmm(
  weights = numeric(0),
  means = list(),
  covariances = list(),
  name = "gmm",
  metadata = list()
)

Arguments

Value

An S7 object inheriting from gmm.

See Also

Other classes: autoplot.gmm_fit, glance.gmm_fit, gmm_counterfactual_law(), gmm_dim(), gmm_fit(), gmm_n_components(), gmm_target(), gmm_weights(), is_proposal(), tidy.gmm

Examples

g <- gmm(
  weights = c(0.4, 0.6),
  means = list(c(-1, 0), c(1, 0)),
  covariances = list(diag(2), diag(2))
)
g

Description

Returns the (closed-form) distribution of $Y = A X + b + \epsilon$ when $X \sim g$ is a Gaussian mixture and $\epsilon \sim \mathcal{N}(0, R)$ is independent additive Gaussian noise.

Usage

gmm_affine(g, A, b = 0, noise_cov = NULL, ridge_eps = 1e-06)

Arguments

Details

For each component k, the pushed-forward parameters are

$\mu'_k = A \mu_k + b, \qquad \Sigma'_k = A \Sigma_k A^\top + R,$

and the mixture weights are unchanged. This is the finite-mixture analogue of a Kalman-style predict step.

The channel is required to be affine in x and the noise is required to be Gaussian. Non-linear channels are not closed form and are not silently approximated: push samples through the map instead (rgmm() then the transform) and refit with fit_em_samples() when a mixture of the image is needed.

Value

A gmm in R^m with the same number of components and the same weights as g.

See Also

Other operators: gmm_aggregate(), gmm_convolve(), gmm_counterfactual(), gmm_filter(), gmm_intervene(), gmm_missing(), gmm_mix(), gmm_observe(), gmm_product(), gmm_reduce()

Examples

g <- gmm(weights = c(0.5, 0.5),
         means = list(c(-1, 0), c(1, 0)),
         covariances = list(diag(2), diag(2)))
A <- matrix(c(1, 0, 0, 1, 1, -1), nrow = 3L, byrow = TRUE)
gmm_affine(g, A, b = c(0, 0, 0), noise_cov = 0.01 * diag(3))

Description

A named alias for gmm_affine() when A is a (row-wise) aggregation matrix — e.g. a block-sum, block-average, or unequal-weight aggregation used in downscaling pipelines. The mathematics is identical to gmm_affine(); the alias gives the public API a clearer hook for aggregation-specific diagnostics in later releases.

Usage

gmm_aggregate(g, A, noise_cov = NULL, ridge_eps = 1e-06)

Arguments

Value

A gmm in R^m.

See Also

Other operators: gmm_affine(), gmm_convolve(), gmm_counterfactual(), gmm_filter(), gmm_intervene(), gmm_missing(), gmm_mix(), gmm_observe(), gmm_product(), gmm_reduce()

Examples

g <- gmm(weights = c(0.5, 0.5),
         means = list(c(-1, 0, 1), c(1, 0, -1)),
         covariances = list(diag(3), diag(3)))
# Sum coordinates 1 and 2 into a single aggregate; pass coord 3 through.
A <- matrix(c(1, 1, 0,
              0, 0, 1), nrow = 2L, byrow = TRUE)
gmm_aggregate(g, A)

Description

Tracks the number of distinct mixture centroids as a function of temperature under mass-constrained deterministic annealing (Rose, Gurewitz and Fox 1990), a physics-derived alternative to information-criterion model selection. The system starts at a high temperature where all k_max centroids collapse to the data centroid (a single effective component) and is cooled along a geometric schedule; at each critical temperature a centroid bifurcates, so the number of distinct centroids grows in steps. The temperatures at which it grows are the phase transitions, and the count occupying the widest temperature range is the discovered component number.

Usage

gmm_anneal_path(
  x,
  k_max = 8L,
  sigma = NULL,
  t_high = NULL,
  t_low = NULL,
  n_steps = 80L,
  n_inner = 30L,
  w = NULL,
  perturb = 0.02,
  merge_tol = 0.1,
  ridge_eps = 1e-06,
  seed = 1L
)

Arguments

Details

The first bifurcation has a closed-form critical temperature $T_c = \lambda_{\max}(\Sigma^{-1} C)$, where $C$ is the (weighted) data covariance and $\Sigma = \sigma^2 I$ the shared reference covariance. This value is returned as t_critical_analytic and serves as an independent analytic check on the empirically detected first transition. Subsequent transitions have no comparably simple closed form, and the count is a diagnostic rather than a guarantee.

Annealing fixes the component covariance to the reference $\Sigma$ so the temperature is the only scale; this is the clean isotropic regime in which the critical temperature is exact. For robust fitting under free covariances, use anneal = TRUE on fit_em_samples() or fit_kld_em() instead.

Value

A list with elements path (a data frame of temperature, n_effective and free_energy), critical_temperatures (the temperatures at which the count increased), first_critical_temperature (the first such, or NA if none was detected), t_critical_analytic ($\lambda_{\max}(\Sigma^{-1} C)$), k_selected (the widest-plateau component count), lambda_max and sigma.

References

Rose, K., Gurewitz, E. and Fox, G. C. (1990) Statistical mechanics and phase transitions in clustering. Physical Review Letters 65(8), 945–948. doi:10.1103/PhysRevLett.65.945

See Also

Other diagnostics: bic_aic(), ess_summary(), ess_trace(), gmm_conditional_entropy(), gmm_entropy(), gmm_evidence(), gmm_fit_ensemble(), gmm_fit_quality(), gmm_independence_graph(), gmm_mutual_information(), hellinger_mc(), kld_trace(), proxy_functional_ci()

Examples

set.seed(1)
x <- rbind(
  matrix(stats::rnorm(120, mean = -4), ncol = 2),
  matrix(stats::rnorm(120, mean =  4), ncol = 2)
)
path <- gmm_anneal_path(x, k_max = 4L, n_steps = 40L)
path$k_selected
path$first_critical_temperature

Description

Returns a new gmm (or gmm_fit) with the components permuted into a canonical order: weight descending, then ⁠||mu||⁠ descending as a tiebreaker. The mixture distribution is unchanged — only the bookkeeping order is — but the canonical ordering removes the EM label-switching nuisance from snapshot tests, cross-run comparisons, and printed summaries.

Usage

gmm_canonicalise(g)

Arguments

Details

Applied automatically by the regime-specific fitters (fit_moment_match(), fit_em_samples(), fit_kld_em()) and by the top-level dispatcher fit_proxymix() when canonicalise = TRUE (the default).

Value

A gmm (or gmm_fit) of the same subclass as g, with the components permuted into canonical order.

See Also

Other ops: dgmm(), gmm_conditionalise(), gmm_divergence(), gmm_kld(), gmm_marginalise(), gmm_mean(), gmm_modes(), pgmm(), rgmm()

Examples

g <- gmm(weights = c(0.1, 0.6, 0.3),
         means = list(c(0, 0), c(3, 0), c(-1, 1)),
         covariances = list(diag(2), diag(2), diag(2)))
gmm_canonicalise(g)

Description

Returns the one identified summary of a gmm_counterfactual_law – its mean.

Usage

gmm_cf_mean(x)

Arguments

Value

Numeric scalar.

See Also

Other decision: fit_uplift(), gmm_cf_tail_prob(), gmm_cf_variance(), gmm_counterfactual(), gmm_intervene(), proxy_cate(), proxy_confounding_gap(), proxy_decide(), proxy_identification_report(), proxy_overlap(), proxy_policy_value(), proxy_predict(), proxy_regime_segments(), proxy_retrospective_uplift(), proxy_uplift(), uplift_identification(), uplift_model()

Examples

g <- gmm(weights = 1, means = list(c(0, 0, 0)),
         covariances = list(diag(3)))
cf <- gmm_counterfactual(g, evidence = c(1, 0, 0.2),
                         do = c(NA, 1, NA), query = 1L)
gmm_cf_mean(cf)

Description

Like gmm_cf_variance(), a per-unit counterfactual tail probability $P(Y_{t'} > c \mid y, t, x)$ is not identified – it is a functional of the unidentified counterfactual law, not of its identified mean. This accessor refuses rather than mislead.

Usage

gmm_cf_tail_prob(x, threshold)

Arguments

Value

Never returns; always raises a proxymix_not_identified error.

See Also

Other decision: fit_uplift(), gmm_cf_mean(), gmm_cf_variance(), gmm_counterfactual(), gmm_intervene(), proxy_cate(), proxy_confounding_gap(), proxy_decide(), proxy_identification_report(), proxy_overlap(), proxy_policy_value(), proxy_predict(), proxy_regime_segments(), proxy_retrospective_uplift(), proxy_uplift(), uplift_identification(), uplift_model()

Examples

g <- gmm(weights = 1, means = list(c(0, 0, 0)),
         covariances = list(diag(3)))
cf <- gmm_counterfactual(g, evidence = c(1, 0, 0.2),
                         do = c(NA, 1, NA), query = 1L)
try(gmm_cf_tail_prob(cf, threshold = 2))

Description

The per-unit counterfactual variance is not identified from the joint density: it depends on the cross-world coupling of the structural residuals under the factual and counterfactual treatments, which no goodness-of-fit can certify. This accessor therefore raises an error rather than return the (misleading) spread of the abduction atoms.

Usage

gmm_cf_variance(x)

Arguments

Value

Never returns; always raises a proxymix_not_identified error.

See Also

Other decision: fit_uplift(), gmm_cf_mean(), gmm_cf_tail_prob(), gmm_counterfactual(), gmm_intervene(), proxy_cate(), proxy_confounding_gap(), proxy_decide(), proxy_identification_report(), proxy_overlap(), proxy_policy_value(), proxy_predict(), proxy_regime_segments(), proxy_retrospective_uplift(), proxy_uplift(), uplift_identification(), uplift_model()

Examples

g <- gmm(weights = 1, means = list(c(0, 0, 0)),
         covariances = list(diag(3)))
cf <- gmm_counterfactual(g, evidence = c(1, 0, 0.2),
                         do = c(NA, 1, NA), query = 1L)
try(gmm_cf_variance(cf))

Description

Extract completed datasets from a gmm_imputation

Usage

gmm_complete(object, which = 1L)

Arguments

Value

When which selects one completion, a single completed dataset (matrix, or data frame if the input was one); otherwise a list of them.

See Also

Other imputation: as_mids(), gmm_imputation(), gmm_impute(), mechanism, proxy_fmi(), proxy_mnar_sensitivity(), proxy_pool()

Description

Returns the differential entropy of the conditional mixture $g_{Y \mid X = x}$ obtained from gmm_conditionalise() – the predictive uncertainty of the target coordinates given the conditioned ones. The order-2 Renyi entropy is closed-form; order = "shannon" falls back to Monte Carlo. Multiple conditioning configurations are evaluated row-by-row.

Usage

gmm_conditional_entropy(
  g,
  given,
  order = c("renyi2", "shannon"),
  n_mc = 5000L,
  seed = NULL
)

Arguments

Value

A numeric scalar for a single configuration, or a numeric vector with one entropy per row of given.

See Also

gmm_entropy(), gmm_conditionalise()

Other diagnostics: bic_aic(), ess_summary(), ess_trace(), gmm_anneal_path(), gmm_entropy(), gmm_evidence(), gmm_fit_ensemble(), gmm_fit_quality(), gmm_independence_graph(), gmm_mutual_information(), hellinger_mc(), kld_trace(), proxy_functional_ci()

Examples

## Joint over (Y, X); predictive entropy of Y at several X values.
s <- matrix(c(2, 0.8, 0.8, 1), 2, 2)
g <- gmm(weights = 1, means = list(c(0, 0)), covariances = list(s))
gmm_conditional_entropy(g, given = rbind(c(NA, 0), c(NA, 1)))

Description

Computes the conditional distribution of a Gaussian mixture given fixed values of a subset of coordinates, by the Schur-complement formula applied component-wise and re-weighted by the marginal evidence $p(\textit{x}_b)$ of each component.

Usage

gmm_conditionalise(g, given)

Arguments

Value

A gmm object in dimension equal to the number of free coordinates.

See Also

Other ops: dgmm(), gmm_canonicalise(), gmm_divergence(), gmm_kld(), gmm_marginalise(), gmm_mean(), gmm_modes(), pgmm(), rgmm()

Examples

g <- gmm(weights = c(0.5, 0.5),
         means = list(c(-1, 0), c(1, 0)),
         covariances = list(diag(2), diag(2)))
gmm_conditionalise(g, given = c(NA, 0.5))

Description

The exact distribution of $X + Y$ for independent $X \sim g_1$ and $Y \sim g_2$: a Gaussian mixture with $K_1 K_2$ components,

$g_1 * g_2 = \sum_{ij} w_i v_j\, \mathcal{N}\!\left(\mu_i + m_j,\ \Sigma_i + S_j\right).$

Usage

gmm_convolve(g1, g2)

Arguments

Details

For the affine special case $X + c$ with a constant c, or $A X + \epsilon$ with Gaussian $\epsilon$, use gmm_affine(); the convolution operator is the general mixture-plus-mixture case that gmm_affine() cannot express.

Value

A gmm with K1 * K2 components.

See Also

Other operators: gmm_affine(), gmm_aggregate(), gmm_counterfactual(), gmm_filter(), gmm_intervene(), gmm_missing(), gmm_mix(), gmm_observe(), gmm_product(), gmm_reduce()

Examples

g1 <- gmm(weights = c(0.5, 0.5), means = list(-1, 1),
          covariances = list(matrix(0.5), matrix(0.5)))
g2 <- gmm(weights = 1, means = list(2), covariances = list(matrix(1)))
gmm_convolve(g1, g2)

Description

Computes the per-unit counterfactual of a query coordinate (the outcome) for an observed unit, under a do() intervention – Pearl's third rung on the latent-class structural causal model read off a fitted Gaussian mixture.

Usage

gmm_counterfactual(g, evidence, do, query, ridge_eps = 1e-06)

Arguments

Details

The three steps are closed form. Abduction recovers the regime posterior $\pi_k(\text{evidence})$ of the observed unit and, within each component, its structural residual. Action sets the do coordinates. Prediction re-evaluates the query coordinate. For a binary treatment the result is a discrete law on K atoms,

$Y_{t'} \mid (y, t, x) = \sum_k \pi_k(y, t, x)\, \delta\!\big(y + \beta_k^{T}(t' - t)\big),$

where $\beta_k^T$ is component k's within-class treatment slope. Only the mean of this law is identified; its spread reflects regime uncertainty, not the (unidentified) cross-world coupling, so the variance and tail accessors refuse to answer (see gmm_cf_variance()).

Value

A gmm_counterfactual_law object carrying the K atoms, their abduction weights, and the identified counterfactual mean.

See Also

gmm_intervene(), gmm_cf_variance(), proxy_retrospective_uplift()

Other operators: gmm_affine(), gmm_aggregate(), gmm_convolve(), gmm_filter(), gmm_intervene(), gmm_missing(), gmm_mix(), gmm_observe(), gmm_product(), gmm_reduce()

Other decision: fit_uplift(), gmm_cf_mean(), gmm_cf_tail_prob(), gmm_cf_variance(), gmm_intervene(), proxy_cate(), proxy_confounding_gap(), proxy_decide(), proxy_identification_report(), proxy_overlap(), proxy_policy_value(), proxy_predict(), proxy_regime_segments(), proxy_retrospective_uplift(), proxy_uplift(), uplift_identification(), uplift_model()

Examples

## Observed unit (y = 1.2, t = 0, x = 0.5); imagine t = 1.
g <- gmm(weights = c(0.6, 0.4),
         means = list(c(0, 0, 0), c(2, 1, 1)),
         covariances = list(diag(3), diag(3)))
cf <- gmm_counterfactual(g, evidence = c(1.2, 0, 0.5),
                         do = c(NA, 1, NA), query = 1L)
cf@mean

Description

The return type of gmm_counterfactual(): a discrete law on K atoms (the per-component counterfactual outcomes) with abduction weights, plus the one identified summary – the counterfactual mean. The atom spread is regime (epistemic) uncertainty about which structural component generated the unit; it is not the counterfactual outcome's dispersion, which is unidentified. Accordingly gmm_cf_variance() and gmm_cf_tail_prob() refuse to answer.

Usage

gmm_counterfactual_law(
  atoms = numeric(0),
  weights = numeric(0),
  mean = numeric(0),
  query = integer(0),
  evidence = numeric(0),
  do = numeric(0),
  name = "gmm_counterfactual_law",
  metadata = list()
)

Arguments

Value

An S7 object of class gmm_counterfactual_law.

See Also

Other classes: autoplot.gmm_fit, glance.gmm_fit, gmm(), gmm_dim(), gmm_fit(), gmm_n_components(), gmm_target(), gmm_weights(), is_proposal(), tidy.gmm

Description

Convenience accessor returning the ambient dimension $p$.

Usage

gmm_dim(x)

Arguments

Value

Integer scalar.

See Also

Other classes: autoplot.gmm_fit, glance.gmm_fit, gmm(), gmm_counterfactual_law(), gmm_fit(), gmm_n_components(), gmm_target(), gmm_weights(), is_proposal(), tidy.gmm

Examples

g <- gmm(weights = 1, means = list(c(0, 0)),
         covariances = list(diag(2)))
gmm_dim(g)

Description

Computes a divergence between two Gaussian mixtures of the same ambient dimension. The Cauchy-Schwarz divergence

$D_{\mathrm{CS}}(p, q) = \tfrac{1}{2}\log V(p, p) + \tfrac{1}{2}\log V(q, q) - \log V(p, q),$

with $V(p, q) = \int p(x) q(x)\, dx$, is closed-form, symmetric, non-negative, and zero exactly when $p \propto q$. The "kl" option delegates to gmm_kld(), a Monte-Carlo estimate of the asymmetric Kullback-Leibler divergence $\mathrm{KL}(p \Vert q)$.

Usage

gmm_divergence(p, q, type = c("cs", "kl"), n_mc = 5000L)

Arguments

Value

For type = "cs", a non-negative numeric scalar. For type = "kl", the list returned by gmm_kld().

See Also

gmm_entropy(), gmm_kld()

Other ops: dgmm(), gmm_canonicalise(), gmm_conditionalise(), gmm_kld(), gmm_marginalise(), gmm_mean(), gmm_modes(), pgmm(), rgmm()

Examples

p <- gmm(weights = c(0.5, 0.5),
         means = list(c(-1, 0), c(1, 0)),
         covariances = list(diag(2), diag(2)))
q <- gmm(weights = 1, means = list(c(0, 0)),
         covariances = list(diag(2) * 2))
gmm_divergence(p, q)
gmm_divergence(p, p)

Description

Computes the differential entropy of a Gaussian mixture. The quadratic (order-2) Renyi entropy $H_2(g) = -\log \int g(x)^2 \, dx$ is available in closed form, because $\int g^2$ is a finite sum of Gaussian-density evaluations. Shannon entropy has no closed form for a mixture (the integrand carries the logarithm of a sum) and is estimated by Monte Carlo, reported with its standard error and an analytic upper bound that brackets it from above.

Usage

gmm_entropy(g, order = c("renyi2", "shannon"), n_mc = 5000L, seed = NULL)

Arguments

Value

For order = "renyi2", a numeric scalar. For order = "shannon", a list with components mc (the estimate), mc_se (its standard error), upper_bound (the analytic upper bound), and n_mc.

See Also

gmm_divergence(), gmm_kld()

Other diagnostics: bic_aic(), ess_summary(), ess_trace(), gmm_anneal_path(), gmm_conditional_entropy(), gmm_evidence(), gmm_fit_ensemble(), gmm_fit_quality(), gmm_independence_graph(), gmm_mutual_information(), hellinger_mc(), kld_trace(), proxy_functional_ci()

Examples

g <- gmm(weights = c(0.5, 0.5),
         means = list(c(-2, 0), c(2, 0)),
         covariances = list(diag(2), diag(2)))
gmm_entropy(g)
gmm_entropy(g, order = "shannon", n_mc = 2000L, seed = 1L)

Description

Tests whether the last m observations of a series are consistent with a linear-Gaussian state-space model fitted on the rest, in the regime where m is small (even m = 1) and smaller than the parameter count – where ordinary structural-break tests (Chow, sup-Wald) are undefined because the post-break parameters cannot be estimated. The statistic is the sum of the last m squared standardised one-step innovations from the filter; a break inflates it.

Usage

gmm_eos_test(
  prior,
  dynamics,
  measurement,
  y,
  m = 1L,
  method = c("chisq", "andrews"),
  alpha = 0.05
)

Arguments

Details

Two calibrations are offered. method = "chisq" refers the statistic to a chi-square distribution on m * ncol(y) degrees of freedom, which is exact when the standardised innovations are Gaussian. method = "andrews" is the distribution-free Andrews (2003) subsampling P-test: the statistic's rank among the in-sample overlapping m-blocks of the same innovations gives the p-value, so it stays calibrated when the innovations are non-Gaussian (heavy-tailed observation noise, say). The model is supplied exactly as for gmm_filter().

Two finite-sample cautions. The chi-square calibration is exact when the model is given; with parameters estimated on a short series it over-rejects (size 0.068 at n = 30 against a nominal 0.05 in the validation study, settling to 0.044 by n = 120), so prefer method = "andrews" when the model is estimated on little data. The subsampling p-value has a floor of ⁠1 / (n - 2m + 2)⁠, so it can reject at level 0.05 only when ⁠n > 2m + 18⁠.

Value

An object of class gmm_eos_test: a list with the statistic, the p_value, the logical reject, the method, m, alpha, and the in-sample block statistics used for the subsampling calibration.

References

Andrews, D. W. K. (2003). End-of-Sample Instability Tests. Econometrica, 71(6), 1661–1694.

See Also

gmm_filter()

Examples

prior <- gmm(weights = 1, means = list(0), covariances = list(matrix(10)))
dyn   <- list(A = matrix(1), Q = matrix(0.04))
meas  <- list(C = matrix(1), R = matrix(1))
set.seed(1)
y <- c(rnorm(119), 6)                # a stable series with a jump at the end
gmm_eos_test(prior, dyn, meas, y, m = 1L, method = "andrews")

Description

Importance-sampling estimate of $Z = \int f(x)\, dx$ for the fit's target $f$, using the fitted mixture $\hat g$ as the proposal:

$\widehat{Z} = \frac{1}{n} \sum_{i=1}^n \frac{f(x_i)}{\hat g(x_i)}, \qquad x_i \sim \hat g,$

computed in the log domain. For a Bayesian posterior handed over as ⁠likelihood x prior⁠, $\log Z$ is the log marginal likelihood, so a fitted proxy doubles as a model-comparison device.

Usage

gmm_evidence(fit, n = 4000L, seed = NULL)

Arguments

Details

The estimator is exact in expectation for any proposal that dominates $f$, and its Monte Carlo error is driven by how well $\hat g$ matches $f$ – which is precisely what the fit optimised. The variance is finite only when $\hat g$ has tails at least as heavy as $f$; a right-tail diagnostic is returned (the effective sample size and the share of the estimate carried by the largest ten percent of weights, which sits near 0.10 for a well-matched proxy), and a classed warning (proxymix_heavy_tail) is raised when it indicates an untrustworthy tail. Results also report the delta-method standard error of $\log \widehat{Z}$.

When the target declares itself normalised (normalised = TRUE), the true value is $\log Z = 0$ and the function still estimates it – a useful end-to-end diagnostic of the fit.

Value

A list of class proxymix_evidence with elements log_z, se_log_z, n, ess, max_weight_share, top_decile_share, and flagged (the heavy-tail indicator).

See Also

Other diagnostics: bic_aic(), ess_summary(), ess_trace(), gmm_anneal_path(), gmm_conditional_entropy(), gmm_entropy(), gmm_fit_ensemble(), gmm_fit_quality(), gmm_independence_graph(), gmm_mutual_information(), hellinger_mc(), kld_trace(), proxy_functional_ci()

Examples

## An unnormalised target with a known constant: log f = log N(., 0, I) + 3.
tgt <- gmm_target(
  n_dim = 2L,
  log_density = function(x) {
    if (is.null(dim(x))) x <- matrix(x, ncol = 2L)
    -0.5 * rowSums(x^2) - log(2 * pi) + 3
  },
  normalised = FALSE, name = "shifted_gaussian"
)
fit <- fit_kld_em(tgt, N = 1L, is_size = 2000L, max_iter = 40L, seed = 1L)
ev <- gmm_evidence(fit, n = 2000L, seed = 2L)
ev$log_z   # close to 3

Description

Runs a Gaussian-sum filter over a series of n observations by alternating the predict operator gmm_affine(), the update operator gmm_observe() and an optional reduction gmm_reduce(). At one component this is the classical Kalman filter; at several components, with Gaussian-mixture process or measurement noise, it is the Gaussian-sum filter of Alspach and Sorenson (1972). The filter is a pure composition of the existing closed-form operators; nothing here leaves the affine-Gaussian world.

Usage

gmm_filter(
  prior,
  dynamics,
  measurement,
  y,
  k_max = NULL,
  reduce = c("merge", "anneal"),
  ridge_eps = 1e-08
)

Arguments

Details

At step $t$ the belief is propagated through the linear-Gaussian dynamics $x_t = A x_{t-1} + b + w$, then conditioned on the observation $y_t = C x_t + d + v$, then (if k_max is set) reduced back to at most k_max components. With Gaussian noise the component count is constant and the recursion is the Kalman filter. Non-Gaussian noise is supplied as a Gaussian sum – a gmm in place of the covariance matrix Q or R. Each such mixture multiplies the component count by its own component count every step, so a long horizon is only runnable with reduction enabled. Reduction is moment-preserving (see gmm_reduce()), so the filtered mean and covariance are unaffected to the moment-matching order while the count stays bounded.

Value

A list with elements

filtered

a length-n list of the filtered gmm beliefs, one per step (after predict, update and any reduction).

mean

an n-by-p matrix of the filtered mixture means.

cov

a length-n list of the p-by-p filtered mixture covariances.

summary

a data frame with one row per step: the step index, the per-coordinate filtered mean (mean_1, ...) and standard deviation (sd_1, ...), the component count n_components, and the step log marginal evidence log_evidence (whose sum over steps is the log likelihood of the series).

References

Alspach, D. L. and Sorenson, H. W. (1972) Nonlinear Bayesian estimation using Gaussian sum approximations. IEEE Transactions on Automatic Control 17(4), 439–448. doi:10.1109/TAC.1972.1100034

See Also

Other operators: gmm_affine(), gmm_aggregate(), gmm_convolve(), gmm_counterfactual(), gmm_intervene(), gmm_missing(), gmm_mix(), gmm_observe(), gmm_product(), gmm_reduce()

Examples

## A one-dimensional local-level filter (random walk observed in noise).
prior <- gmm(weights = 1, means = list(0), covariances = list(matrix(1)))
truth <- cumsum(stats::rnorm(20, sd = 0.3))
y <- truth + stats::rnorm(20, sd = 0.5)
out <- gmm_filter(
  prior,
  dynamics    = list(A = matrix(1), Q = matrix(0.09)),
  measurement = list(C = matrix(1), R = matrix(0.25)),
  y = y
)
head(out$summary)

## Heavy-tailed process noise as a two-component Gaussian sum, capped at
## four components per step.
q_noise <- gmm(weights = c(0.9, 0.1),
               means = list(0, 0),
               covariances = list(matrix(0.05), matrix(1)))
out2 <- gmm_filter(
  prior,
  dynamics    = list(A = matrix(1), Q = q_noise),
  measurement = list(C = matrix(1), R = matrix(0.25)),
  y = y, k_max = 4L
)
max(out2$summary$n_components)

Description

A gmm_fit is the result of fit_proxymix() (or one of the regime-specific fitters). It inherits the mixture parameters of gmm and adds a record of the target it was fitted to, the regime used, and the iteration diagnostics.

Usage

gmm_fit(
  weights = numeric(0),
  means = list(),
  covariances = list(),
  name = "gmm",
  metadata = list(),
  target = NULL,
  regime = NA_character_,
  diagnostics = list(),
  converged = NA,
  iterations = NA_integer_,
  call = NULL
)

Arguments

Value

An S7 object inheriting from gmm_fit (and gmm).

See Also

Other classes: autoplot.gmm_fit, glance.gmm_fit, gmm(), gmm_counterfactual_law(), gmm_dim(), gmm_n_components(), gmm_target(), gmm_weights(), is_proposal(), tidy.gmm

Examples

samples <- matrix(stats::rnorm(200), ncol = 2)
tgt <- gmm_target_from_samples(samples)
fit <- fit_proxymix(tgt, N = 2L, regime = "sample", max_iter = 25L)
inherits(fit, "proxymix::gmm_fit")

Description

Quantifies the sampling variability of the fitted mixture itself by a Bayesian (weighted) bootstrap: each replicate re-weights the fit's own observations with Dirichlet(1, ..., 1) weights and refits by a warm-started weighted EM. In regime "kld" the observations are the fit's cached importance draws with their self-normalised weights, so a replicate costs zero new target evaluations; in regimes "sample" and "moment" the observations are the target's samples. Summaries of any functional of the proxy are then read off the ensemble with proxy_functional_ci() – functional-space intervals sidestep the label-switching that makes parameter-space intervals incoherent.

Usage

gmm_fit_ensemble(fit, B = 200L, max_iter = 25L, tol = 1e-05, seed = NULL)

Arguments

Value

A list of class gmm_ensemble: fit (the base fit), members (a length-B list of gmm), B, and regime.

References

Rubin, D. B. (1981) The Bayesian bootstrap. The Annals of Statistics 9(1), 130–134.

See Also

Other diagnostics: bic_aic(), ess_summary(), ess_trace(), gmm_anneal_path(), gmm_conditional_entropy(), gmm_entropy(), gmm_evidence(), gmm_fit_quality(), gmm_independence_graph(), gmm_mutual_information(), hellinger_mc(), kld_trace(), proxy_functional_ci()

Examples

fit <- fit_proxymix(banana_target(), N = 2L, regime = "kld",
                    is_size = 1500L, max_iter = 20L, seed = 1L)
ens <- gmm_fit_ensemble(fit, B = 30L, seed = 2L)
proxy_functional_ci(ens, gmm_mean)

Description

Returns the fit-quality certificate: a small list recording the fitting regime, convergence, degeneracy, the effective-sample-size profile of the importance weights, and (when a validation split was drawn) the held-out validation gap. The certificate is stamped into the object's metadata at fit time and carried unchanged through every closed-form operator, so it can be read off a marginal, a conditional, a filtered belief, or any other derived mixture – alongside the provenance vector recording the chain of operations that produced it.

Usage

gmm_fit_quality(g)

Arguments

Details

Downstream verbs read the same certificate and raise a one-shot advisory (class proxymix_low_quality) when the source fit is flagged.

Value

A list with elements regime, converged, degenerate, ess, ess_relative, min_component_ess, max_weight, support_fraction, kld_final, and validation_gap (fields not applicable to the regime are NA), or NULL for a mixture that was never fitted (e.g. built directly with gmm()).

See Also

Other diagnostics: bic_aic(), ess_summary(), ess_trace(), gmm_anneal_path(), gmm_conditional_entropy(), gmm_entropy(), gmm_evidence(), gmm_fit_ensemble(), gmm_independence_graph(), gmm_mutual_information(), hellinger_mc(), kld_trace(), proxy_functional_ci()

Examples

fit <- fit_proxymix(banana_target(), N = 2L, regime = "kld",
                    is_size = 1500L, max_iter = 15L, seed = 1L)
gmm_fit_quality(fit)
## The certificate survives the operator calculus.
gmm_fit_quality(gmm_marginalise(fit, keep = 1L))

Description

The object returned by gmm_impute(). It carries the m completed data matrices, the bootstrap-fitted mixtures behind them (used by the analytic pooling in proxy_pool()), the mixture fitted to the full data, and a record of the missingness. Pass it to gmm_complete() to extract the completed datasets and to proxy_pool() / proxy_fmi() for inference.

Usage

gmm_imputation(
  data = NULL,
  completions = list(),
  fits = list(),
  point_fit = NULL,
  n_components = integer(0),
  m = integer(0),
  mechanism = "mar",
  observed = NULL,
  var_names = character(0),
  is_data_frame = FALSE,
  diagnostics = list(),
  call = NULL
)

Arguments

Value

An S7 object of class gmm_imputation.

See Also

Other imputation: as_mids(), gmm_complete(), gmm_impute(), mechanism, proxy_fmi(), proxy_mnar_sensitivity(), proxy_pool()

Description

Fits a Gaussian mixture to a numeric dataset that contains missing values and draws m completed datasets from the mixture conditional $p(x_{\mathrm{missing}} \mid x_{\mathrm{observed}})$. Because the mixture can be multimodal and heteroscedastic, the imputations follow the shape of the joint distribution rather than a single Gaussian, which keeps downstream inference valid on data that a single-Gaussian or linear-Gaussian imputer mis-specifies.

Usage

gmm_impute(
  data,
  N = NULL,
  m = 20L,
  mechanism = mar(),
  seed = NULL,
  max_iter = 100L,
  tol = 1e-06,
  ridge_eps = 1e-06
)

Arguments

Details

Imputation is conditioning. For a row with observed coordinates the missing coordinates follow the closed-form mixture conditional (the same Schur-complement algebra as gmm_conditionalise()). The mixture is fitted to the incomplete data by expectation-maximisation whose E-step uses each row's observed margin and whose M-step restores the conditional covariance of the filled entries, so component variances are not under-estimated.

Proper multiple imputation requires the fitting parameters themselves to carry uncertainty, otherwise the pooled intervals are too narrow. Each of the m imputations is therefore drawn under a mixture fitted to an independent bootstrap resample of the rows, so proxy_pool() reflects both imputation and parameter uncertainty.

The mechanism says how an entry came to be missing, which sets the conditional the missing value is drawn from: mar() (the default) for missing at random, censored() for a known interval such as a detection limit, or mnar() for a value-dependent selection model. The interval and value-dependent gates act on a single coordinate, and a row missing that coordinate must have its other coordinates observed. Numeric data only; categorical variables are out of scope.

Value

A gmm_imputation object.

See Also

gmm_complete() to extract completions, proxy_pool() to pool an estimand across them, gmm_conditionalise() for the conditioning algebra.

Other imputation: as_mids(), gmm_complete(), gmm_imputation(), mechanism, proxy_fmi(), proxy_mnar_sensitivity(), proxy_pool()

Examples

set.seed(1)
x1 <- rnorm(200)
x2 <- x1 + rnorm(200)
x2[runif(200) < plogis(x1)] <- NA          # missing at random on x1
imp <- gmm_impute(cbind(x1, x2), N = 1L, m = 10L, seed = 1L)
proxy_pool(imp, "x2")$estimate             # pooled mean of x2

Description

Returns the undirected second-order conditional-independence graph of a fitted Gaussian mixture: the partial-correlation (Gaussian graphical model) structure of the mixture's overall covariance. An edge $i - j$ is present when the partial correlation of coordinates $i$ and $j$ given all the others exceeds threshold in magnitude, and absent when it does not – the latter is the Markov statement $x_i \perp x_j \mid x_{\mathrm{rest}}$ at second order. The overall covariance

$\mathrm{Cov}(X) = \sum_k w_k (\Sigma_k + \mu_k \mu_k^\top) - \bar\mu\,\bar\mu^\top, \qquad \bar\mu = \sum_k w_k \mu_k,$

is closed-form in the mixture parameters, so no sampling is required.

Usage

gmm_independence_graph(g, threshold = 0.05)

Arguments

Details

This is a graphical-model (dependency-structure) diagnostic, not a causal discovery method: it recovers the undirected Markov skeleton, not edge directions. Its distinctive use is regime (iii): composed with fit_kld_em(), it recovers the dependency structure of a target you can only evaluate (an unnormalised energy / Gibbs density), where no sample exists to hand a sampling-based estimator. Being second order, it sees dependencies that enter the covariance; a purely higher-order coupling (zero correlation, nonzero dependence) is not detected – raise the mixture's component count, or read the coordinate-block dependence with gmm_mutual_information() instead.

Value

A symmetric integer adjacency matrix (1 = edge, 0 = none) with the coordinate names of g, carrying the partial-correlation matrix as the "pcor" attribute.

See Also

gmm_mutual_information(), fit_kld_em(), gmm_conditionalise()

Other diagnostics: bic_aic(), ess_summary(), ess_trace(), gmm_anneal_path(), gmm_conditional_entropy(), gmm_entropy(), gmm_evidence(), gmm_fit_ensemble(), gmm_fit_quality(), gmm_mutual_information(), hellinger_mc(), kld_trace(), proxy_functional_ci()

Examples

## Regime (iii): the Markov structure of an evaluable-but-unsampleable density.
## A continuous chain field x1 - x2 - x3 (couplings only between neighbours).
energy <- function(X) {
  X <- matrix(X, ncol = 3)
  rowSums((X^2 - 1)^2) - 0.7 * (X[, 1] * X[, 2] + X[, 2] * X[, 3])
}
target <- gmm_target(n_dim = 3L, log_density = function(X) -energy(X))
g <- fit_kld_em(target, N = 8L, proposal = is_uniform(3L, -3, 3),
                is_size = 8000L, anneal = TRUE, seed = 1L, support_warn = FALSE)
gmm_independence_graph(g)            # recovers x1 - x2 - x3 (no x1 - x3 edge)

Description

Reads a fitted Gaussian mixture as a latent-class structural causal model and returns the interventional distribution of the free coordinates under do() of some coordinates, optionally conditioning on others.

Usage

gmm_intervene(g, do, given = NULL, ridge_eps = 1e-06)

Arguments

Details

Intervened (do) coordinates are set inside every component but do not re-weight the regime gate – this is the graph surgery that distinguishes $p(\cdot \mid do(T = t))$ from $p(\cdot \mid T = t)$. Conditioned (given) coordinates re-weight the gate in the usual Bayesian way. Writing the component prior as $\pi_k$, the within-component conditional mean as $\mu_k$, and the given-coordinate evidence as $e_k$, the returned mixture has weights $\pi_k(given) \propto \pi_k\, e_k$ and per-component parameters from the Schur conditional on the union of the do and given coordinates.