mixle.epistemic.discrepancy module

Distance/divergence between two distributions, or between a predicted and an observed sample.

This is the “compare predicted vs. observed” hinge an epistemic loop’s UPDATE arrow needs: given a hypothesis’s predicted observation and the real one, how far apart are they? Nothing under mixle.stats/mixle.inference computed this directly before this module – proper scoring rules (mixle.inference.scoring) score a single outcome against a predictive distribution, which is a related but different question (they answer “how good was this one call”, not “how far apart are these two whole distributions”).

Every function here is generic over any object exposing log_density/sample (the same duck-typed surface mixle.capability already dispatches on) or .sampler(seed).sample(n) (the concrete shape every mixle.stats distribution actually has). A closed-form fast path is used only where one is exact and unambiguous (currently: two univariate Gaussians); everything else falls back to a Monte Carlo estimate, and discrepancy_report() says plainly which path was taken via its degraded flag – an honest signal, never a silently approximated number presented as exact.

class DiscrepancyResult(value, metric, degraded, seed=None, n_samples=None, reproducible=True)[source]

Bases: object

One discrepancy evaluation: the value, which metric computed it, and whether it was exact.

seed and n_samples record what a Monte Carlo estimate needs for exact reproduction: pass seed back into discrepancy_report() and the same value comes back out. Both are None when nothing was sampled – an exact closed form, or a metric applied directly to caller-supplied arrays – since recording a seed there would misleadingly imply randomness that was never used. seed is also None when the caller passed their own np.random.RandomState instance rather than an int: that object’s reproducibility is already the caller’s to manage, and there is no single integer that reconstructs it.

reproducible is False when a sampled side draws from randomness this module cannot control – a bare .sample(n) that takes no RNG/seed argument and reads NumPy’s global state. A seed passed to such a call never reaches the sampler, so two identically seeded reports return genuinely different values; seed is therefore also None in that case rather than recording an integer that does not reproduce anything. Give the sampler an rng / random_state / generator / seed keyword, or expose mixle’s .sampler(seed) shape, to make a path reproducible.

Parameters:
discrepancy_report(predicted, observed, *, metric='auto', seed=None)[source]

The actual delta_m(o_hat, o) entry point: compare a predicted and an observed value/distribution.

metric="auto" picks kl_divergence when both sides look like distributions (expose log_density), else mmd over raw arrays (the “predicted is a distribution, observed is a concrete measurement” case reduces to comparing observed against samples drawn from predicted). degraded=True whenever the underlying computation fell back to a Monte Carlo / sample-based estimate rather than an exact closed form – callers that need to know whether a number is exact or estimated read this field rather than guessing from the metric name.

Every path that draws samples does so from a seeded RNG, never from unrecorded OS entropy: pass seed for a run you choose to make reproducible, or leave it unset and read the seed this function picked back off result.seed – either way, calling again with that same seed reproduces the exact same value, PROVIDED the sampled objects let this module control their randomness. When one does not, result.reproducible is False and result.seed is None: the guarantee is reported as unavailable rather than quietly broken. See DiscrepancyResult for that contract and for when n_samples comes back None.

Parameters:
Return type:

DiscrepancyResult

kl_divergence(p, q, *, n=10_000, seed=None)[source]

KL(p || q) in nats: exact closed form when a known pair matches, else a Monte Carlo estimate.

The one closed-form entry in the dispatch table today is two univariate Gaussians (the exact formula, not an approximation); every other pair falls back to mean_{x ~ p}[log p(x) - log q(x)] using n samples drawn from p. Extending the closed-form table to more conjugate pairs (Categorical-Categorical, Dirichlet-Dirichlet, …) is legitimate future work – it was deliberately left at one entry here rather than half-built across several families with incompatible parameterizations (mixle’s categorical distribution keys its simplex by a pmap over arbitrary hashable labels, not a fixed-order probability vector, which is a real complication left to a dedicated follow-up rather than papered over).

Parameters:
Return type:

float

js_divergence(p, q, *, n=10_000, seed=None)[source]

Jensen-Shannon divergence: symmetric, bounded, computed via the sample-mixture estimator.

0.5 * KL(p || m) + 0.5 * KL(q || m) where m is the equal mixture of p and q; each term is estimated by sampling from the corresponding side and evaluating log m(x) = log(0.5 p(x) + 0.5 q(x)) via logaddexp for numerical stability. Symmetric by construction up to Monte Carlo noise (both halves use independent sample draws).

Parameters:
Return type:

float

wasserstein_distance(p, q, *, n=10_000, seed=None)[source]

1-Wasserstein (earth-mover) distance between two 1D distributions, via sorted sample matching.

Draws n samples from each side; the empirical 1D optimal transport cost is the mean absolute difference between the two sorted sample sequences (exact for the empirical distributions, a consistent estimator of the true distance as n grows). Samples are squeezed to plain 1D before sorting (_as_1d_samples()) rather than sorted-then-raveled – sorting a (n, 1) column vector in its original shape is a silent no-op (see that function’s docstring), which would silently return the mean absolute difference between the two unsorted sample sequences instead of the true empirical Wasserstein distance. Raises NotImplementedError for genuinely multivariate input rather than silently computing a coordinate-wise number that isn’t the true multivariate Wasserstein distance – there is no cheap exact estimator for that case, and returning a wrong-but-plausible-looking number would be worse than refusing.

Parameters:
Return type:

float

mmd(samples_p, samples_q, *, kernel='rbf', bandwidth=None)[source]

Maximum Mean Discrepancy between two raw sample sets: sqrt(max(mmd_squared(...), 0)).

mmd_squared() estimates squared MMD and is deliberately signed, so returning it under the name mmd published a “distance” that could be negative – identical inputs scored below zero, and because a discrepancy feeds DiscrepancyLikelihood’s exp(-discrepancy / temperature), a negative value produced a likelihood above 1 and reversed the intended penalty ordering right around zero, where the two distributions are most alike.

Taking the square root of the clipped estimate restores the contract the public discrepancy API is written against: non-negative, zero when the samples are indistinguishable, growing with separation, and in the same units as the data. The clip means this is not unbiased near zero (it can only overstate a true zero, never understate it); when unbiasedness is what matters – two-sample testing, permutation nulls – call mmd_squared() directly and handle its sign yourself.

Parameters:
Return type:

float

mmd_squared(samples_p, samples_q, *, kernel='rbf', bandwidth=None)[source]

The unbiased estimator of squared MMD between two raw sample sets – a SIGNED quantity.

This is Gretton et al.’s U-statistic, and unbiasedness is exactly why it can come out negative: under the null (both sample sets from the same distribution) the true value is 0, so an unbiased estimator of it must scatter to both sides. Identical arrays [0, 1] give -0.3934693403; identical singletons give -2.0. Those are correct outputs of this estimator, not errors – which is why it is named for what it estimates. Use it when you want the unbiased test statistic (a two-sample test threshold, a permutation null); use mmd() when you want a discrepancy that behaves like a distance.

Unlike the other functions here, this takes samples directly rather than distribution objects – it works even when neither side is a mixle.stats distribution (e.g. a real observation array vs. a synthesized/predicted one). bandwidth defaults to the median pairwise distance heuristic over the pooled samples. Only the RBF kernel is implemented; other kernel names raise NotImplementedError.

Parameters:
Return type:

float