mixle.stats.multivariate.multivariate_gaussian module

Multivariate Gaussian distributions over real-valued vectors.

Data type: np.ndarray[float]

x = (x_1,x_2,..,x_n) ~ MVN(mu, covar), where mu is a length n numpy array, and covar is an n by n positive definite covariance matrix.

The log-density is given by

log(p(x)) = -0.5*k*log(2*pi) - 0.5*log|covar| - 0.5*(x-mu)’ covar^{-1} (x-mu).

Reference: Mardia, Kent & Bibby, Multivariate Analysis (Academic Press, 1979).

class MultivariateGaussianFisherView(dist)

Bases: FixedFisherView

Fisher view over first and upper-triangular second moments for a full Gaussian.

Parameters:

dist (Any)

class MultivariateGaussianDistribution(mu, covar=MISSING, name=None, keys=None, covariance=MISSING, prior=None)

Bases: SequenceEncodableProbabilityDistribution

Multivariate normal distribution with mean vector mu and full covariance matrix covar.

Parameters:

• mu (list[float] | ndarray)

• covar (list[list[float]] | ndarray)

• name (str | None)

• keys (str | None)

• covariance (list[list[float]] | ndarray)

• prior (SequenceEncodableProbabilityDistribution | None)

classmethod compute_capabilities()

Declare backend support for multivariate Gaussian generated kernels.

classmethod compute_declaration()

Return the generated-compute declaration for the multivariate Gaussian.

static exp_family_sufficient_statistics(x, engine)

Return vector/matrix sufficient statistics for generated MVN scoring.

Parameters:

• x (Any)

• engine (Any)

Return type:

tuple[Any, …]

static exp_family_natural_parameters(params, engine)

Return natural parameters for generated MVN scoring.

Parameters:

• params (dict[str, Any])

• engine (Any)

Return type:

tuple[Any, …]

static exp_family_log_partition(params, engine)

Return the full-covariance Gaussian log partition.

Parameters:

• params (dict[str, Any])

• engine (Any)

Return type:

Any

static backend_legacy_sufficient_statistics(x, params, engine)

Return row-wise legacy accumulator statistics for generated resident reductions.

Parameters:

• x (Any)

• params (dict[str, Any])

• engine (Any)

Return type:

tuple[Any, …]

set_prior(prior)

Attach a parameter prior and precompute conjugate-prior expectations.

With a NormalWishart(m0, kappa, W, nu) prior over (mu, Lambda=covar^-1) this caches the prior parameters and E[ln|Lambda|], the quantities needed by expected_log_density. Any other prior (including None) leaves the distribution a plain point model.

Parameters:

prior (SequenceEncodableProbabilityDistribution | None)

Return type:

None

expected_log_density(x)

Variational expectation E_q[log p(x | mu, Lambda)] under the NormalWishart prior.

Falls back to the plug-in log_density(x) when no conjugate prior is attached.

Return type:

float

seq_expected_log_density(x)

Vectorized expected_log_density over sequence-encoded observations.

Parameters:

x (ndarray)

Return type:

ndarray

density(x)

Evaluate the density at x.

Parameters:

x (np.ndarray) – Observation from multivariate Gaussian distribution.

Returns:

Density at x.

Return type:

float

log_density(x)

Evaluate the log-density at x.

The log-density is given by

log(p(x)) = -0.5*k*log(2*pi) - 0.5*log|covar| - 0.5*(x-mu)’ covar^{-1} (x-mu).

Parameters:

x (np.ndarray) – Observation from multivariate Gaussian distribution.

Returns:

Log-density at x.

Return type:

float

density_cumulative(x)

Exact probability-ordered cumulative G(x) = P(p(Y) >= p(x)) – the highest-density-region mass whose boundary passes through x (the multivariate analogue of a CDF; a coordinate-wise CDF is undefined without a total order on R^d).

For a multivariate Gaussian p(y) >= p(x) iff the squared Mahalanobis distance is no larger, and that distance is chi-square with dim degrees of freedom, so G(x) = chi2.cdf(maha2, dim). Used by mixle.enumeration.density_rank.density_rank() to return an EXACT cumulative for the MVN.

Parameters:

x (ndarray)

Return type:

float

density_quantile(q)

Inverse of density_cumulative(): a representative point at cumulative-density index q.

The multivariate analogue of a quantile / inverse-CDF: a coordinate-wise quantile is undefined without a total order on R^d, but the density ordering gives one. q is the highest-density region mass, so the boundary is the squared-Mahalanobis level chi2.ppf(q, dim); we return a representative point on that contour, mu + sqrt(level) * L[:, 0] where covar = L L^T (so its Mahalanobis distance is exactly the level). Sweeping q enumerates the support in descending density.

Parameters:

q (float)

Return type:

ndarray

seq_log_density(x)

Vectorized evaluation of the log-density at a sequence-encoded input x.

Parameters:

x (np.ndarray) – Encoded data matrix with shape (sz, dim) from MultivariateGaussianDataEncoder.seq_encode().

Returns:

Numpy array of length sz containing the log-density of each encoded observation.

Return type:

ndarray

static backend_log_density_from_params(x, mu, inv_covar, log_det, engine)

Engine-neutral multivariate Gaussian log-density from inverse covariance.

Parameters:

• x (Any)

• mu (Any)

• inv_covar (Any)

• log_det (Any)

• engine (Any)

Return type:

Any

backend_seq_log_density(x, engine)

Engine-neutral vectorized log-density for encoded data.

Parameters:

• x (Any)

• engine (Any)

Return type:

Any

classmethod backend_stacked_params(dists, engine)

Return stacked full-covariance Gaussian parameters for a homogeneous mixture kernel.

Parameters:

• dists (Sequence[MultivariateGaussianDistribution])

• engine (Any)

Return type:

dict[str, Any]

classmethod backend_stacked_log_density(x, params, engine)

Return an (n, k) matrix of full-covariance Gaussian log densities.

Parameters:

• x (Any)

• params (dict[str, Any])

• engine (Any)

Return type:

Any

classmethod backend_stacked_sufficient_statistics(x, weights, params, engine)

Return component-stacked legacy sufficient statistics on the active engine.

The weighted second moment sum_n w[n,k] x_n x_n^T is accumulated per component as a gemm (x * w[:, k]).T @ x instead of reducing a per-sample, per-component (n, k, dim, dim) outer-product tensor. That intermediate is N*K*dim*dim (~20 GB at n=2e4, k=8, dim=128 — it OOMs a GPU); the per-component gemm holds only an (n, dim) temporary and hands the reduction to BLAS/cuBLAS. The first moment is likewise w.T @ x rather than a reduced (n, k, dim) tensor.

Parameters:

• x (Any)

• weights (Any)

• params (dict[str, Any])

• engine (Any)

Return type:

tuple[Any, …]

to_fisher(**kwargs)

Return this distribution’s own Fisher view.

sampler(seed=None)

Return a sampler for iid draws from this distribution.

Parameters:

seed (int | None) – Optional seed for the sampler’s random state.

Returns:

A configured MultivariateGaussianSampler.

condition(observed)

Return the conditional distribution over the unobserved dimensions given observed.

observed maps dimension index to its fixed value; the result is the closed-form Gaussian conditional over the remaining dimensions (in increasing index order):

mu_{u|o} = mu_u + Sigma_uo Sigma_oo^{-1} (x_o - mu_o), Sigma_{u|o} = Sigma_uu - Sigma_uo Sigma_oo^{-1} Sigma_ou.

Sampling the result is given=-style conditional sampling (draw the unobserved coordinates consistent with the observed ones). Raises if no dimension is left unobserved.

Parameters:

observed (dict[int, float])

Return type:

MultivariateGaussianDistribution

marginal(keep)

Return the marginal Gaussian over dimensions keep.

N(mu, Sigma) marginalized to index set keep is N(mu[keep], Sigma[keep, keep]). The order of keep is preserved, so the result’s dimensions follow the given order.

Parameters:

keep (Sequence[int])

Return type:

MultivariateGaussianDistribution

estimator(pseudo_count=None)

Return an estimator initialized from this distribution’s shape.

If pseudo_count is passed, the current mean and covariance are used to regularize the estimate.

Parameters:

pseudo_count (float | None) – Optional smoothing count applied to mean and covariance estimates.

Returns:

A MultivariateGaussianEstimator.

dist_to_encoder()

Return an encoder for iid multivariate Gaussian observations.

Return type:

MultivariateGaussianDataEncoder

class MultivariateGaussianSampler(dist, seed=None)

Bases: DistributionSampler

Sampler for iid multivariate Gaussian observations.

Parameters:

• dist (MultivariateGaussianDistribution)

• seed (int | None)

sample(size=None, *, batched=True)

Draw iid samples from the multivariate Gaussian distribution.

Parameters:

• size (Optional[int]) – Number of iid samples to draw. If None, a single sample is drawn.

• batched (bool)

Returns:

Numpy array with shape (dim,) if size is None, else with shape (size, dim).

Return type:

ndarray

class MultivariateGaussianSuffStat(sum_, sum2_, count_, anchored=None)

Bases: tuple

A (sum, sum_outer, count) sufficient statistic that also carries shift-anchored moments.

Behaves exactly like the plain 3-tuple everywhere it is indexed, unpacked, or iterated (it is one), so generic consumers – scale_suff_stat, the declared StatisticSpec reader, engine kernels, _assert_suff_close – see nothing new. anchored is extra payload: (anchor, sum_i w_i*(x_i - anchor), sum_i w_i*(x_i - anchor)(x_i - anchor)^T), which MultivariateGaussianAccumulator.combine() folds in and MultivariateGaussianEstimator.estimate() uses to compute a shift-invariant covariance. The raw moments stay the exchange format – (sum, sum2, count) is the declared statistic tuple – so the anchored track rides along beside them rather than replacing them. Mirrors mixle.stats.univariate.continuous.gaussian.GaussianSuffStat.

Parameters:

• sum_ (ndarray)

• sum2_ (ndarray)

• count_ (float)

• anchored (tuple | None)

class MultivariateGaussianAccumulator(dim=None, keys=None, name=None)

Bases: SequenceEncodableStatisticAccumulator

Accumulator for multivariate Gaussian sufficient statistics.

Alongside the declared raw moments (sum, sum2, count) this keeps a SHIFT-ANCHORED moment track, for the reason the univariate GaussianAccumulator keeps one: the covariance computed from raw reduced moments is the classic cancellation-prone E[xx^T] - mu mu^T form, which loses ~2*log2(abs(mean)/sd) bits, so data with unit spread at offset 1e7 fits a visibly wrong correlation and data at offset 1e8 fits a matrix so far from positive-definite that _robust_cho_factor refuses it – on data whose covariance a two-pass computation returns exactly (epoch seconds are ~1.7e9). Anchoring at the first value seen keeps every term of the scatter O(count * spread^2), making the M-step covariance shift-invariant. The track is CONDITIONING-GATED (see _needs_anchor()): a chunk the raw form handles to ~1e-9 relative error accumulates exactly the historical single-pass way – bit-identical statistics, no second pass – and the anchor activates only when a chunk (or a scalar update) would corrupt the covariance.

Parameters:

• dim (int | None)

• keys (str | None)

• name (str | None)

update(x, weight, estimate)

Update sufficient statistics with a single weighted observation.

Parameters:

• x (np.ndarray) – Length-dim observation vector.

• weight (float) – Weight for the observation.

• estimate (Optional[MultivariateGaussianDistribution]) – Kept for consistency with SequenceEncodableStatisticAccumulator (not used).

Returns:

None.

Return type:

None

initialize(x, weight, rng)

Initialize the accumulator with a weighted observation. Calls update().

Parameters:

• x (np.ndarray) – Length-dim observation vector.

• weight (float) – Weight for the observation.

• rng (Optional[RandomState]) – Kept for consistency with SequenceEncodableStatisticAccumulator.

Returns:

None.

Return type:

None

seq_update(x, weights, estimate)

Vectorized update of sufficient statistics with an encoded sequence of observations.

Parameters:

• x (np.ndarray) – Encoded data matrix with shape (sz, dim).

• weights (np.ndarray) – Numpy array of sz observation weights.

• estimate (Optional[MultivariateGaussianDistribution]) – Kept for consistency (not used).

Returns:

None.

Return type:

None

seq_initialize(x, weights, rng)

Vectorized initialization of the accumulator. Calls seq_update().

Parameters:

• x (np.ndarray) – Encoded data matrix with shape (sz, dim).

• weights (np.ndarray) – Numpy array of sz observation weights.

• rng (Optional[RandomState]) – Kept for consistency with SequenceEncodableStatisticAccumulator.

Returns:

None.

Return type:

None

combine(suff_stat)

Merge sufficient statistics into this accumulator.

Parameters:

suff_stat (Tuple[np.ndarray, np.ndarray, float]) – Tuple of (weighted sum of observations, weighted sum of outer products, sum of weights).

Returns:

This accumulator.

Return type:

MultivariateGaussianAccumulator

value()

Return (sum, sum_outer, count) sufficient statistics.

Once the shift-anchored moment track is live the returned value is a MultivariateGaussianSuffStat – a drop-in 3-tuple (indexing/unpacking/iteration all behave identically) that additionally carries those moments in its anchored attribute, so combine() can fold them in and MultivariateGaussianEstimator.estimate() can compute a shift-invariant covariance. Well-conditioned data never activates the track and gets the historical plain tuple.

Return type:

tuple[ndarray, ndarray, float]

from_value(x)

Replace this accumulator’s sufficient statistics.

Parameters:

x (Tuple[np.ndarray, np.ndarray, float]) – Tuple of (weighted sum of observations, weighted sum of outer products, sum of weights).

Returns:

This accumulator.

Return type:

MultivariateGaussianAccumulator

scale(c)

Scale the accumulated statistics in-place by c, anchored track included.

The structural default routes through from_value(scale_suff_stat(self.value(), c)), which rebuilds from a PLAIN tuple and would therefore drop the anchored payload – turning a scaled large-offset accumulator back into a cancellation-prone one. Scaling every weight by c scales both anchored moments by c and leaves the anchor (a data value) alone.

Parameters:

c (float)

Return type:

MultivariateGaussianAccumulator

acc_to_encoder()

Return an encoder compatible with this accumulator’s dimension.

Return type:

MultivariateGaussianDataEncoder

class MultivariateGaussianAccumulatorFactory(dim, keys=None, name=None)

Bases: StatisticAccumulatorFactory

Factory for multivariate Gaussian accumulators.

Parameters:

• dim (int | None)

• keys (str | None)

• name (str | None)

make()

Return a fresh accumulator with the factory configuration.

Return type:

MultivariateGaussianAccumulator

class MultivariateGaussianEstimator(dim=None, pseudo_count=(None, None), suff_stat=(None, None), name=None, keys=None, prior=None, min_covar=None, ridge=None, track_conditioning=False, degenerate_ratio=1.0e-6)

Bases: ParameterEstimator

Estimator for multivariate Gaussian distributions.

The fitted covariance is NOT the raw maximum-likelihood estimate unless ridge=0.0 is passed: by default the M-step regularizes each diagonal entry of the empirical covariance as cov_ii + eps_i with eps_i = ridge * cov_ii (default ridge=1e-6) so a singular empirical covariance – e.g. an EM component holding fewer than d points – remains factorizable. Each eps_i is RELATIVE to its own coordinate’s variance, so the ridge is a negligible ~1e-6 relative perturbation whatever the units, per column: the fit is equivariant under any per-coordinate change of units, and heterogeneous-unit data is no longer penalized (the former trace-scaled ridge priced a microscale column’s jitter off the macroscale columns, inflating the smallest variances – and with them the log-likelihood and AIC – by up to ~5 orders of magnitude, T1-05). A coordinate with no positive variance of its own gets its jitter from the mean of the coordinates that have one; the absolute min_covar (default 1e-8) is the last resort for a covariance with no scale at all (all-zero diagonal, or ridge=0.0) and becomes a hard per-coordinate floor as soon as it is passed explicitly.

A materially binding ridge is recorded on the fitted distribution as a covariance-ridged(...) entry in numerical_repairs() – and therefore in fit_provenance().repairs. “Materially” means the jitter exceeded _RIDGE_MATERIAL_RATIO of the smallest raw variance (an explicit min_covar dominating a small variance, or manufactured jitter on a variance-free coordinate – never the default per-coordinate ridge on a full-rank fit), or the raw covariance was not positive-definite on its own and the ridge is what makes it factorable – the rank-deficient case this ridge exists for. Pass ridge=0.0 to recover the exact MLE up to the min_covar absolute jitter floor.

The fit is SHIFT-EQUIVARIANT: fit(x + c) returns the covariance of x + c to within a few ulps for any constant c the data can carry, epoch seconds (~1.7e9) included. That is not free from the declared (sum, sum2, count) statistics – their covariance is the cancellation-prone E[xx^T] - mu mu^T, which at offset 1e7 fitted a visibly wrong correlation and at 1e8 a matrix so far from positive-definite that the constructor refused it – so MultivariateGaussianAccumulator carries a conditioning-gated shift-anchored moment track and this estimator reads it. Statistics that arrive already reduced and WITHOUT that track (an engine kernel’s stacked moments, a hand-built tuple) cannot be corrected here and take the historical path; when they are too ill-conditioned for it, estimate warns rather than returning a covariance it cannot stand behind.

Parameters:

• dim (int | None)

• pseudo_count (float | tuple[float | None, float | None] | None)

• suff_stat (tuple[ndarray | None, ndarray | None] | None)

• name (str | None)

• keys (str | None)

• prior (SequenceEncodableProbabilityDistribution | None)

• min_covar (float | None)

• ridge (float | None)

• track_conditioning (bool)

• degenerate_ratio (float)

accumulator_factory()

Return an accumulator factory matching this estimator.

Return type:

MultivariateGaussianAccumulatorFactory

model_log_density(model)

Log-density of the model parameters under the NormalWishart prior (ELBO global term).

The prior is over (mu, Lambda=covar^-1), so the model’s covariance is inverted before scoring.

Parameters:

model (MultivariateGaussianDistribution)

Return type:

float

estimate(nobs, suff_stat)

Estimate a multivariate normal distribution with from aggregated sufficient statistics.

Suff_stat is a Tuple of size 3 containing:

suff_stat[0] (np.ndarray): Component-wise sum of weighted observation values. suff_stat[1] (np.ndarray): Component-wise sum of weighted squared observation values. suff_stat[2] (float): Sum of weights for each observation.

Parameters:

• nobs (Optional[float]) – Weighted number of observations used in aggregation of suff stats.

• suff_stat (Tuple[np.ndarray, np.ndarray, float]) – See above for details.

Returns:

MultivariateGaussianDistribution

Return type:

MultivariateGaussianDistribution

class MultivariateGaussianDataEncoder(dim=None)

Bases: DataSequenceEncoder

Encoder for iid multivariate Gaussian observations.

Parameters:

dim (int | None)

seq_encode(x)

Encode a sequence of iid length-dim observations for vectorized seq_ calls.

Parameters:

x (Union[Sequence[List[float]], Sequence[List[np.ndarray]], np.ndarray]) – Sequence of length-dim observation vectors.

Returns:

Encoded data matrix with shape (len(x), dim).