mixle.stats.multivariate.multivariate_gaussian module¶

Create, estimate, and sample from a multivariate normal distribution with mean vector ‘mu’ (length n), and covariance matrix ‘covar’ (n by n).

Defines the MultivariateGaussianDistribution, MultivariateGaussianSampler, MultivariateGaussianAccumulatorFactory, MultivariateGaussianAccumulator, MultivariateGaussianEstimator, and the MultivariateGaussianDataEncoder classes for use with mixle.

Data type: np.ndarray[float]

x = (x_1,x_2,..,x_n) ~ MVN(mu, covar), where mu is a length n numpy array, anc 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)[source]

Bases: FixedFisherView

Parameters:: dist (Any)

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

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()[source]

classmethod compute_declaration()[source]

static exp_family_sufficient_statistics(x, engine)[source]

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)[source]

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)[source]

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)[source]

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)[source]

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)[source]

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)[source]

Vectorized expected_log_density over sequence-encoded observations.

Parameters:: x (ndarray)
Return type:: ndarray

density(x)[source]

Evaluate the density at x.

Parameters:: x (np.ndarray) – Observation from multivariate Gaussian distribution.
Returns:: Density at x.
Return type:: float

log_density(x)[source]

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)[source]

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)[source]

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)[source]

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)[source]

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)[source]

Engine-neutral vectorized log-density for encoded data.

Parameters:

x (Any)
engine (Any)

Return type:

Any

classmethod backend_stacked_params(dists, engine)[source]

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)[source]

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)[source]

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)[source]: Return this distribution’s own Fisher view.

sampler(seed=None)[source]

Create a MultivariateGaussianSampler for sampling from this distribution.

Parameters:: seed (Optional[int]) – Seed to set for sampling with RandomState.
Returns:: MultivariateGaussianSampler object.

condition(observed)[source]

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)[source]

Return the marginal Gaussian over the dimensions keep (Gaussian marginals just drop rows).

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)[source]

Create a MultivariateGaussianEstimator for estimating this distribution.

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

Parameters:: pseudo_count (Optional[float]) – Used to inflate sufficient statistics in estimation.
Returns:: MultivariateGaussianEstimator object.

dist_to_encoder()[source]

Returns a MultivariateGaussianDataEncoder object for encoding sequences of iid observations.

Return type:: MultivariateGaussianDataEncoder

class MultivariateGaussianSampler(dist, seed=None)[source]

Bases: DistributionSampler

MultivariateGaussianSampler object for sampling from a MultivariateGaussianDistribution.

Parameters:

dist (MultivariateGaussianDistribution)
seed (int | None)

sample(size=None)[source]

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.
Returns:: Numpy array with shape (dim,) if size is None, else with shape (size, dim).
Return type:: ndarray

class MultivariateGaussianAccumulator(dim=None, keys=None, name=None)[source]

Bases: SequenceEncodableStatisticAccumulator

MultivariateGaussianAccumulator object for aggregating sufficient statistics from iid observations.

Parameters:

dim (int | None)
keys (str | None)
name (str | None)

update(x, weight, estimate)[source]

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)[source]

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)[source]

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)[source]

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)[source]

Merge the sufficient statistics of suff_stat 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:: MultivariateGaussianAccumulator object.
Return type:: MultivariateGaussianAccumulator

value()[source]

Returns the sufficient statistics (sum, sum of outer products, count) of the accumulator.

Return type:: tuple[ndarray, ndarray, float]

from_value(x)[source]

Set the sufficient statistics of the accumulator to x.

Parameters:: x (Tuple[np.ndarray, np.ndarray, float]) – Tuple of (weighted sum of observations, weighted sum of outer products, sum of weights).
Returns:: MultivariateGaussianAccumulator object.
Return type:: MultivariateGaussianAccumulator

key_merge(stats_dict)[source]

Combine sufficient statistics with other accumulators sharing a matching key.

Parameters:: stats_dict (Dict[str, Any]) – Dictionary mapping keys to aggregated statistics.
Returns:: None.
Return type:: None

key_replace(stats_dict)[source]

Replace sufficient statistics with values from stats_dict for a matching key.

Parameters:: stats_dict (Dict[str, Any]) – Dictionary mapping keys to aggregated statistics.
Returns:: None.
Return type:: None

acc_to_encoder()[source]

Returns a MultivariateGaussianDataEncoder object for encoding sequences of iid observations.

Return type:: MultivariateGaussianDataEncoder

class MultivariateGaussianAccumulatorFactory(dim, keys=None, name=None)[source]

Bases: StatisticAccumulatorFactory

MultivariateGaussianAccumulatorFactory object for creating MultivariateGaussianAccumulator objects.

Parameters:

dim (int | None)
keys (str | None)
name (str | None)

make()[source]

Returns a new MultivariateGaussianAccumulator with the factory’s dim, keys, and name.

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)[source]

Bases: ParameterEstimator

MultivariateGaussianEstimator object for estimating a multivariate normal distribution from aggregated sufficient statistics.

Parameters:

dim (int | None)
pseudo_count (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)

accumulator_factory()[source]

Returns a MultivariateGaussianAccumulatorFactory built from the estimator’s attributes.

Return type:: MultivariateGaussianAccumulatorFactory

model_log_density(model)[source]

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)[source]

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)[source]

Bases: DataSequenceEncoder

MultivariateGaussianDataEncoder object for encoding sequences of iid multivariate Gaussian observations.

Parameters:: dim (int | None)

seq_encode(x)[source]

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).