mixle.stats.matrix.lkj module

LKJ distribution over correlation matrices.

The Lewandowski-Kurowicka-Joe (LKJ) law places density f(R) = c_d(eta) * det(R)^(eta-1) on d x d correlation matrices (symmetric, unit diagonal, positive definite). The concentration eta > 0 tilts mass toward the identity: eta = 1 is uniform over correlation matrices, eta > 1 favours weak correlations (R near I), and eta < 1 favours strong ones. It is the standard prior on correlation matrices in hierarchical Bayesian models (Stan’s default), separating a covariance into scales times a correlation.

Normalizer (C-vine derivation, verified to high precision against arbitrary-precision integration over the correlation elliptope for d = 2, 3):

Z_d(eta) = prod_{k=1}^{d-1} B(eta + (d-1-k)/2, 1/2)^(d-k), c_d(eta) = 1 / Z_d(eta).

It samples exactly by the onion method (Lewandowski et al. 2009, sec. 3.2): each off-diagonal entry then has the exact marginal (r + 1)/2 ~ Beta(eta + (d-2)/2, eta + (d-2)/2). Because the density depends on R only through det(R), observations are encoded as log det(R), and eta is fit by maximum likelihood (a 1-D root find: the mean log-determinant equals sum_k (d-k)[psi(eta+e_k) - psi(eta+e_k+1/2)] with e_k = (d-1-k)/2).

Reference: Lewandowski, Kurowicka & Joe, “Generating random correlation matrices based on vines and extended onion method”, J. Multivariate Analysis 100 (2009).

class LKJDistribution(dim, eta, name=None, keys=None)

Bases: SequenceEncodableProbabilityDistribution

LKJ distribution over dim x dim correlation matrices with concentration eta > 0.

Parameters:

• dim (int)

• eta (float)

• name (str | None)

• keys (str | None)

density(x)

Return the probability density at a correlation matrix x.

Parameters:

x (Any)

Return type:

float

log_density(x)

Return the log-density at a dim x dim correlation matrix (-inf if not positive definite).

Parameters:

x (Any)

Return type:

float

seq_log_density(x)

Return vectorized log-density for a sequence-encoded array of log det(R) values.

Parameters:

x (ndarray)

Return type:

ndarray

sampler(seed=None)

Return an onion-method sampler for correlation matrices.

Parameters:

seed (int | None)

Return type:

LKJSampler

estimator(pseudo_count=None)

Return a maximum-likelihood estimator for the concentration eta (dim fixed).

Parameters:

pseudo_count (float | None)

Return type:

LKJEstimator

dist_to_encoder()

Return the data encoder (a correlation matrix is encoded as its log-determinant).

Return type:

LKJDataEncoder

class LKJSampler(dist, seed=None)

Bases: DistributionSampler

Sample correlation matrices by the onion method (Lewandowski-Kurowicka-Joe 2009).

Parameters:

• dist (LKJDistribution)

• seed (int | None)

sample(size=None)

Draw observations.

Combinator samplers (mixture/sequence/…) accept batched. With batched=True (the default) each child stream is drawn in one vectorized call instead of a per-draw Python loop – far faster. Because every child sampler owns an independent RandomState, batching consumes each stream in the same order as the loop, so the draws are identical to the legacy path. batched=False forces that legacy per-draw loop as a guaranteed- stable reference. Leaf samplers are already vectorized and ignore the flag.

Parameters:

size (int | None)

Return type:

Any

class LKJAccumulator(name=None, keys=None)

Bases: SequenceEncodableStatisticAccumulator

Accumulate (count, sum of log det(R)) – the sufficient statistics for the eta-MLE.

Parameters:

• name (str | None)

• keys (str | None)

update(x, weight, estimate)

Parameters:

• x (Any)

• weight (float)

• estimate (LKJDistribution | None)

Return type:

None

initialize(x, weight, rng)

Parameters:

• x (Any)

• weight (float)

• rng (RandomState | None)

Return type:

None

seq_update(x, weights, estimate)

Parameters:

• x (ndarray)

• weights (ndarray)

• estimate (Any)

Return type:

None

seq_initialize(x, weights, rng)

Parameters:

• x (ndarray)

• weights (ndarray)

• rng (RandomState | None)

Return type:

None

combine(suff_stat)

Parameters:

suff_stat (tuple[float, float])

Return type:

LKJAccumulator

value()

Return type:

tuple[float, float]

from_value(x)

Parameters:

x (tuple[float, float])

Return type:

LKJAccumulator

key_merge(stats_dict)

Pool this accumulator’s statistics into stats_dict under its merge key.

The structural default implements the common single-key pattern: store the accumulator under self.keys the first time the key is seen, else combine into the one already there. Accumulators with several named keys (e.g. an HMM’s init/trans/state keys) or a non-accumulator stats payload override this. A keys of None (the default) is a no-op.

Parameters:

stats_dict (dict[str, Any])

Return type:

None

key_replace(stats_dict)

Replace this accumulator’s statistics from the pooled stats_dict entry (see key_merge).

Parameters:

stats_dict (dict[str, Any])

Return type:

None

acc_to_encoder()

Return type:

LKJDataEncoder

class LKJAccumulatorFactory(name=None, keys=None)

Bases: StatisticAccumulatorFactory

Factory for LKJAccumulator.

Parameters:

• name (str | None)

• keys (str | None)

make()

Return type:

LKJAccumulator

class LKJEstimator(dim, eta_bounds=(0.05, 1.0e4), name=None, keys=None)

Bases: ParameterEstimator

Maximum-likelihood estimator for the concentration eta at fixed dimension dim.

Parameters:

• dim (int)

• eta_bounds (tuple[float, float])

• name (str | None)

• keys (str | None)

accumulator_factory()

Return type:

LKJAccumulatorFactory

estimate(nobs, suff_stat)

Parameters:

• nobs (float | None)

• suff_stat (tuple[float, float])

Return type:

LKJDistribution

class LKJDataEncoder

Bases: DataSequenceEncoder

Encode each correlation matrix as its log-determinant (the only data the density depends on).

seq_encode(x)

Encode the iid observation sequence x for vectorized evaluation.

Parameters:

x (Sequence[Any])

Return type:

ndarray