mixle.ppl package

mixle.ppl — an elegant probabilistic-programming surface over mixle.

A model is plain mixle construction with two new things allowed in a parameter slot: the token free (estimate it) or another distribution (make it random). Fit with .fit(data); query with .sample / .log_prob / .posterior.

from mixle.ppl import Normal, free m = Normal(free, free).fit(data) m.sample(100)

The mixle.stats distribution classes are untouched; this is a thin, optional dialect over the existing model, estimator, and inference contracts.

Layout: the modeling surface and constraints are in mixle.ppl.core; the dialect constructors (Normal/Gamma/Markov/…) in mixle.ppl.distributions; the family/composite lowering registry in mixle.ppl._lowering (imported for its registration side effects); and fields, regression, conformal prediction, survival, diagnostics, and experimental scaling-law research in their own modules. Physics, simulation, and PDE-inverse systems are intentionally outside the Mixle Core PPL boundary. This package init only re-exports and labels stable versus experimental APIs.

class RandomVariable(kind, *, family=None, args=(), name=None, keys=None, dist=None, result=None, scope='shared', reparam=None, group_by=None)

Bases: object

The single user-facing PPL type (immutable).

Two states: sample (a symbolic draw: a family + argument expressions, some of which may be free) and bound (wraps a concrete fitted distribution). The verb surface is fixed and state-independent; validity depends on state. Construct via the family functions in mixle.ppl or fit.

Parameters:

result (Any | None)

property certificate

The estimation certificate, when a fit attached one.

Penalized fits downgrade the certificate because the optimum is for a surrogate objective, not an unpenalized likelihood.

each(by=None)

Mark this prior as per-group (a random effect / local latent). Used in a parameter slot: Normal(Normal(m, t).each(), s) is a hierarchical model.

Two data layouts are supported:

• nested – each() with no argument: .fit(groups) where groups is a list of per-group observation lists (one list per group).

• indexed-flat – each(by="g"): .fit(y, given={"g": labels}) where y is one flat observation array and labels[i] is observation i’s group. This is the varying-intercepts / 8-schools idiom; groups are taken in sorted order of the unique labels.

Parameters:

by (str | None)

Return type:

RandomVariable

noncentered()

Sample this location-scale prior in non-centered form (offset/multiplier).

For mu = Normal(loc, scale) with a random scale (a hierarchical prior), the centered parameterization couples mu’s range to scale and creates Neal’s funnel – a geometry HMC/NUTS samples badly. Normal(loc, scale).noncentered() instead samples a standard normal z and sets mu = loc + scale * z, whose geometry is independent of scale. Mathematically identical posterior, far better mixing (fewer divergences) when the data are weakly informative. Applies to Normal priors; a no-op marker on others.

Return type:

RandomVariable

property scope: str

Return the variable scope, such as scalar, grouped, or global.

classmethod fitted(dist, *, name=None, result=None)

Wrap a concrete distribution as the FITTED random variable a fitter returns.

This is what every @register_fitter function has to hand back, and the only way to build one was a private classmethod – so the shipped extensibility-seam example could not be written against the public surface at all (P09-F13). result is the optional posterior object a sampling fitter attaches; name defaults to the distribution’s own.

Parameters:

• name (str | None)

• result (Any | None)

Return type:

RandomVariable

independent()

Return an explicitly independent copy of one atomic random variable.

Reusing the same object in an expression means reusing the same draw (x - x == 0). Use x.independent() when an identically distributed but independent draw is intended.

Return type:

RandomVariable

exp()

Return the deterministic exponential transform of this random variable.

Return type:

RandomVariable

log()

Return the deterministic logarithm transform of this random variable.

Return type:

RandomVariable

eq(other)

Build an equality constraint against a value or another expression.

ne(other)

Build an inequality constraint against a value or another expression.

given(constraint)

Condition this RV on a constraint over itself (e.g. x.given(x > 0) -> truncation). The result samples by rejection and scores with the renormalized density. For relations among several RVs (a < b) use constrain(...).

Return type:

RandomVariable

property is_bound: bool

Whether this variable has already been lowered to a fitted concrete object.

property has_free: bool

Whether this parameter/expression graph contains one or more inferable unknowns.

property name: str | None

Return the optional user-visible variable name.

property parameter_dimension: int | None

Declared inferable degrees of freedom, or None when a backend owns an unknown-size model.

A fitted artifact returns the count recorded before execution; an unfitted structural model is counted directly. Fixed constructor values and derived posterior summaries are excluded.

property group_index: dict | None

Stable label-to-posterior-index mapping for an indexed-group fit.

property group_labels: tuple | None

Group labels in stable first-observation order for an indexed-group fit.

property columns: list

the variable names, in sample-column order. A vector-valued variable expands to one name per entry (v[0], v[1], …).

Type:

For a constrain(...) joint RV

property dist

The lowered concrete distribution — the full original mixle API (escape hatch).

property components

Fitted sub-models of a composite (mixture components, HMM state emissions, sequence element) as RandomVariables — query each with the same verbs (.params, .sample, .log_prob). Raises for non-composite models.

property params

Fitted parameters in the same parameterization used to construct the model (e.g. {'mean': 5.0, 'sd': 2.0} for Normal — not the internal sigma2). Falls back to .dist for families without a registered reader.

property result: Any | None

Inference metadata (EM history / MCMC chain) when present; else None.

sample(n=None, seed=None, size=None, *, max_attempts=100)

Draw samples from the represented distribution or derived expression.

n and size are aliases (size matches the stats-layer samplers); pass at most one. None returns a single draw.

Parameters:

• n (int | None)

• seed (int | None)

• size (int | None)

• max_attempts (int)

log_prob(x)

Evaluate the log probability or log density at x.

log_density(x)

Alias of log_prob() – the mixle.stats density verb, so a random variable answers the same call a fitted distribution does.

log_likelihood(data)

Total log-likelihood of data under the fitted model.

For an ordinary fit this is the sum of per-observation log_prob. For a conditional (regression / mixed-effects) fit it is the MARGINAL log-likelihood recorded by the fit – for a mixed model the random effects are integrated out, so the likelihood factors per GROUP, not per observation, and no per-observation sum would be correct.

Return type:

float

aic(data, k=None)

Akaike information criterion (lower is better).

k defaults to the declared inferable dimension recorded by fit(); fixed constructor values and derived summaries are never counted. Models without such a record must pass k.

Parameters:

k (int | None)

Return type:

float

bic(data, k=None)

Bayesian information criterion (lower is better).

Parameters:

k (int | None)

Return type:

float

pointwise_log_likelihood(data)

Return the (n_draws, n_obs) log-likelihood matrix used by WAIC / PSIS-LOO.

For a Bayesian fit (how='mcmc'|'hmc'|'ensemble'|'vi') each row is the log-likelihood of the data under one posterior draw. Point estimates have no posterior-draw matrix; use plugin_log_likelihood() for their ordinary per-observation plug-in score.

Return type:

ndarray

plugin_log_likelihood(data)

Per-observation log likelihood under the fitted point estimate (no posterior integration).

Return type:

ndarray

waic(data)

Widely Applicable Information Criterion from the posterior (lower waic is better).

Returns {elpd_waic, p_waic, waic, se, n_draws, pointwise}. Estimates out-of-sample predictive accuracy by integrating over parameter uncertainty – the Bayesian analogue of aic/bic. Point-estimate fits are rejected.

Return type:

dict

loo(data)

Pareto-Smoothed Importance-Sampling Leave-One-Out cross-validation (lower loo better).

Returns {elpd_loo, p_loo, loo, se, khat_max, n_draws, pointwise}. khat_max above ~0.7 signals an unreliable estimate (refit with more posterior draws or prefer waic).

Return type:

dict

summary()

Posterior summary of a Bayesian fit, or the fitted params for a point estimate.

For how='mcmc'|'hmc'|'ensemble'|'vi' returns a per-parameter dict of {mean, std, q2.5, q97.5, mcse} (the 95% credible interval plus each mean’s own Monte Carlo standard error, STAT-RR19-11) plus _acceptance_rate and, for multi-chain runs, _rhat / _ess / _n_chains. For map/em it returns .params.

mean(samples=20000, seed=0)

Expected value, using an analytic distribution moment before bounded Monte Carlo.

Parameters:

• samples (int)

• seed (int)

var(samples=20000, seed=0)

Variance, using an analytic distribution moment before bounded Monte Carlo.

Parameters:

• samples (int)

• seed (int)

prob(samples=40000, seed=999)

Monte Carlo probability receipt for a constrain(...) relation.

Parameters:

• samples (int)

• seed (int)

prob_of_event(samples=40000, seed=999)

Monte Carlo probability receipt for the event underlying a conditioned RV.

Parameters:

• samples (int)

• seed (int)

predict(n=1, rng=None)

Posterior-predictive draws. For a Bayesian fit (conjugate/mcmc/hmc) this integrates over parameter uncertainty (draw params from the posterior, then data); for a point fit (EM/MAP) it is the plug-in predictive (sample from the fitted distribution).

For a fitted REGRESSION, pass the covariates to predict at – a mapping of field name to column, or a DataFrame – and one draw per row comes back. A fitted regression is a conditional law with no single distribution to sample from, so predict(count) on one had nothing to do and a caller who passed covariates got “predictive sample count must be an exact positive integer” (P08-F15).

Parameters:

n (int | Mapping[str, Any])

posterior(x)

Posterior over a latent or a parameter.

• posterior(data) -> latent-state posterior (the E-step; e.g. mixture responsibilities), routed to the lowered distribution’s seq_posterior.

• posterior(handle | name | index) -> parameter posterior draws, when this RV was fit with how='mcmc' (read from .result).

explain_fit(*, how='auto', constraints=None, potentials=None, **_)

Report which inference route .fit(how=...) took (or would take, before fitting).

Returns {'route', 'reason', 'caveats'}. This is the inspection surface for Mixle’s automatic cross-family inference selection: rv.explain_fit() answers how the expression will be fit, what result type it returns, and which diagnostics or limitations apply. The route mirrors fit() exactly by sharing _resolve_auto() for the flat tree and re-checking the same structural short-circuits.

Called on a bound RV (the result of .fit(...)), this reports what that fit actually did – fit() stashes its own answer to this question, computed while the pre-fit expression still carried its priors, since a bound RV’s _args is always empty and cannot be re-derived from. The record travels with the model through pickling, so a reloaded artifact still explains how it was fit. Raises only when no record exists at all – a bound model built directly rather than through .fit(), or one written before the record was carried – in which case call explain_fit() on the pre-fit expression instead.

Return type:

dict

fit(data, *, how='auto', max_its=100, delta=1e-8, backend='local', num_workers=None, engine=None, precision=None, print_iter=0, missing='error', **kw)

Estimate / infer parameters from data and return a bound RV.

how: 'em' (EM/MLE, default for plain free models), 'map' (maximize the joint with priors), 'mcmc' (posterior samples over parameters with priors), 'auto' picks map when the model has priors else em. EM threads mixle’s parallel/distributed backends (backend='mp'|'mpi'|'dask').

missing: 'error' (default) rejects non-finite entries; 'marginalize' integrates a missing entry (NaN in the data) out of the likelihood instead of imputing it – each leaf is fit from its present rows only, so you get a well-defined mode/posterior over the present data (no fabricated values). Supported on the EM path (the default for free models, i.e. the posterior mode under flat priors); for how='map'/'mcmc' with missing data build the model with mixle.stats.marginalized() leaves directly.

Parameters:

• data (Sequence[Any])

• how (str)

• max_its (int)

• delta (float)

• backend (str)

• num_workers (int | None)

• engine (Any)

• precision (Any)

• print_iter (int)

• missing (str)

Return type:

RandomVariable

lower(rv, *, target='dist')

The one routing site: symbolic RandomVariable -> existing mixle object.

target='dist' returns a concrete *Distribution (needs no free holes); target='estimator' returns a *Estimator. Results are cached per random variable.

Parameters:

• rv (RandomVariable)

• target (str)

class Guide(**latents)

Bases: object

A declared mean-field variational approximation: named latents, each an independent q-factor.

Build it with Guide(name=handle, ...) where each handle is the same RandomVariable object used in the model (latents are matched by object identity). To pin and validate the variational family, pass name=(handle, 'gaussian'|'gamma'|'dirichlet'); otherwise the conjugate factor implied by the latent’s prior is used.

Parameters:

latents (Any)

names()

Return the declared latent names in guide order.

Return type:

tuple[str, …]

class StructuredVIPosterior(gres, guide)

Bases: object

Posterior from structured_vi(): per-latent variational factors + the ELBO trace.

posterior(name) returns the factor’s hyperparameters (e.g. {'mean','sd'} for a Gaussian, {'alpha','mean'} for a Dirichlet, {'shape','rate','mean'} for a Gamma); mean / samples give the posterior mean / exact draws of that latent. Names are the guide’s keys (the latent handle itself is also accepted).

Parameters:

guide (Guide)

posterior(name)

Return posterior parameters for a declared latent name or raw handle.

Return type:

dict

mean(name)

Return the posterior mean for a declared latent.

samples(name, n=4000, rng=None)

Draw samples from the variational factor for a declared latent.

Parameters:

n (int)

summary()

Return per-latent posterior parameters and ELBO metadata.

Return type:

dict

structured_vi(observations, guide, *, max_its=300, tol=1e-8)

Fit a structured model by mean-field VMP / coordinate-ascent VI under a declared Guide.

Parameters:

• observations – a list of (model, data) pairs – each model is a PPL RandomVariable observation factor, data its observed values. Factors that reuse the same latent handle share that latent (their evidence is combined). A single factor may be passed as (model, data) directly.

• guide (Guide) – the Guide naming the latents to approximate and (optionally) their q-families.

• tol (float) – CAVI sweep budget and ELBO convergence tolerance.

• max_its (int)

• tol

Returns:

A StructuredVIPosterior over the guide’s latents (monotone ELBO).

Raises:

ValueError – if a guide latent is not actually an inferred latent of the model, or its declared q-family does not match the model’s conjugate factor (the projection constraint is checked).

Return type:

StructuredVIPosterior

admixture(docs, topics, *, alpha=1.0, max_its=100, inner_its=40, tol=1e-5, seed=0)

Fit an admixture (LDA-class) model by mean-field VI built from this surface’s Dirichlet primitives.

docs is a corpus – a list of documents, each a sequence of integer word ids over the vocabulary. topics are the Dirichlet RandomVariable handles you declare as the per-topic variational factors q(beta_k); their prior alpha vectors set the vocabulary size V and the topic-word prior eta. alpha is the per-document topic Dirichlet prior. Returns an AdmixturePosterior.

This is “LDA via the guide”: LDA = an admixture whose emission is Categorical over words. The same coordinate-ascent (q(theta_d), q(beta_k) Dirichlet + categorical responsibilities q(z)) fits any admixture over a categorical vocabulary – no LDADistribution involved.

Parameters:

• max_its (int)

• inner_its (int)

• tol (float)

• seed (int)

class AdmixturePosterior(lam, gamma, topic_handles, ll_trace)

Bases: object

Posterior from admixture() (LDA-class mean-field VI). posterior(topic_handle) returns that topic’s fitted Dirichlet {'alpha','mean'}; topics() is E[beta] (K x V); doc_topics is the per-document mixing E[theta_d]; log_likelihood is the fitted-model corpus LL trace.

topics()

Return the posterior mean topic-word matrix E[beta].

doc_topics(d=None)

Return posterior mean document-topic weights for one document or all documents.

posterior(topic)

Return the fitted Dirichlet posterior for a topic handle or topic index.

summary()

Return corpus-level topic count, vocabulary size, and likelihood metadata.

Return type:

dict

Normal(mean, sd, *, name=None, keys=None)

Symbolic normal distribution parameterized by mean and standard deviation.

The user-facing scale is sd. Lowering converts it to the variance parameter used by GaussianDistribution while fitted artifacts should remain interpretable in the constructor’s scale parameterization.

Core-catalog mapping: Normal(mean, sd) is GaussianDistribution(mu=mean, sigma2=sd**2) – the catalog’s second argument is the VARIANCE, not the standard deviation (see the module-level parameterization map).

Parameters:

• mean (Any)

• sd (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Poisson(rate, *, name=None, keys=None)

Symbolic Poisson count distribution parameterized by non-negative rate.

Parameters:

• rate (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Gamma(shape, rate, *, name=None, keys=None)

Symbolic Gamma distribution parameterized by shape and rate.

Core-catalog mapping: the catalog’s GammaDistribution(k, theta) uses SCALE theta; lowering maps rate to theta = 1 / rate (see the module-level parameterization map).

Parameters:

• shape (Any)

• rate (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Exponential(rate, *, name=None, keys=None)

Symbolic exponential distribution parameterized by event rate.

Core-catalog mapping: the catalog’s ExponentialDistribution(beta) uses the MEAN/scale beta = 1 / rate (see the module-level parameterization map).

Parameters:

• rate (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Categorical(probs=None, *, logits=None, dim=None, name=None, keys=None)

Categorical from a probability dict {value: p} or a list of probabilities. The probability vector is also an inferable parameter. Categorical(free) learns the category probabilities by maximum likelihood, discovering the categories (and their count) from the data – no dim= needed; pass dim=K to request the explicit simplex-parameter treatment for how='mcmc'|'ensemble'|'map'.

Categorical(logits=Net(out=K)) is neural classification: p(y|x) = softmax(Net(x)), the softmax-link sibling of logistic regression. Fit with the conditional verb .fit(y, given={"x": X}).

Parameters:

• probs (Any)

• logits (Any)

• dim (int | None)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Bernoulli(p, *, name=None, keys=None)

Bernoulli distribution for binary outcomes with success probability p.

Parameters:

• p (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Geometric(p, *, name=None, keys=None)

Symbolic geometric count distribution with success probability p.

Parameters:

• p (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Binomial(n, p, *, name=None, keys=None)

Symbolic binomial distribution with n trials and success probability p.

Core-catalog mapping: the catalog’s BinomialDistribution(p, n) takes the SAME two parameters in the opposite positional order (see the module-level parameterization map); keyword arguments are safe on both surfaces.

Parameters:

• n (Any)

• p (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Weibull(shape, scale, *, name=None, keys=None)

Symbolic Weibull distribution with positive shape and scale parameters.

Parameters:

• shape (Any)

• scale (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Laplace(loc, scale, *, name=None, keys=None)

Symbolic Laplace distribution with location and positive scale.

Parameters:

• loc (Any)

• scale (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Logistic(loc, scale, *, name=None, keys=None)

Symbolic logistic distribution with location and positive scale.

Parameters:

• loc (Any)

• scale (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Uniform(low, high, *, name=None, keys=None)

Symbolic continuous uniform distribution on the closed interval [low, high].

Parameters:

• low (Any)

• high (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Rayleigh(sigma, *, name=None, keys=None)

Symbolic Rayleigh distribution with positive scale sigma.

Parameters:

• sigma (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Pareto(scale, shape, *, name=None, keys=None)

Symbolic Pareto distribution with minimum value scale and tail index shape.

Parameters:

• scale (Any)

• shape (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Beta(a, b, *, name=None, keys=None)

Beta distribution over probabilities with concentration parameters a and b.

Parameters:

• a (Any)

• b (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

StudentT(df, loc, scale, *, name=None, keys=None)

Symbolic Student-t distribution with degrees of freedom, location, and scale.

Parameters:

• df (Any)

• loc (Any)

• scale (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

LogNormal(mu, sigma, *, name=None, keys=None)

Symbolic log-normal distribution where log(X) is normal.

Naming: this family is called LogGaussian in the core catalog (mixle.dist.LogGaussianDistribution, alias mixle.dist.LogNormalDistribution); LogGaussian() is this constructor’s same-object alias, so either name works on either surface (T2-06). Core-catalog mapping: sigma here is the log-space STANDARD DEVIATION; the catalog class takes the log-space variance sigma2 = sigma**2 (T2-05).

Parameters:

• mu (Any)

• sigma (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

LogGaussian(mu, sigma, *, name=None, keys=None)

Symbolic log-normal distribution where log(X) is normal.

Naming: this family is called LogGaussian in the core catalog (mixle.dist.LogGaussianDistribution, alias mixle.dist.LogNormalDistribution); LogGaussian() is this constructor’s same-object alias, so either name works on either surface (T2-06). Core-catalog mapping: sigma here is the log-space STANDARD DEVIATION; the catalog class takes the log-space variance sigma2 = sigma**2 (T2-05).

Parameters:

• mu (Any)

• sigma (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

EMG(mu, sigma, rate, *, name=None, keys=None)

Exponentially-modified Gaussian: X = Normal(mu, sigma) + Exponential(rate) (right-skewed).

Lowers to ExponentiallyModifiedGaussianDistribution(mu, sigma**2, lam=rate); rate is the exponential component’s rate (its mean is 1/rate). Note the catalog class takes the VARIANCE sigma2 = sigma**2 where this constructor takes the standard deviation sigma (see the module-level parameterization map). The MLE is iterative with no closed form, so EMG(free, free, free).fit(data) uses a consistent method-of-moments estimate.

Parameters:

• mu (Any)

• sigma (Any)

• rate (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

NegativeBinomial(r, p, *, name=None, keys=None)

Symbolic negative-binomial distribution with shape r and success probability p.

Parameters:

• r (Any)

• p (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

HalfNormal(sigma, *, name=None, keys=None)

Half-normal on [0, inf) with scale sigma – the standard weakly-informative scale prior.

Parameters:

• sigma (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

InverseGamma(alpha, beta, *, name=None, keys=None)

Inverse-gamma(alpha, beta) – the classic conjugate prior for a variance.

Parameters:

• alpha (Any)

• beta (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

InverseGaussian(mu, lam, *, name=None, keys=None)

Inverse-Gaussian (Wald) with mean mu and shape lam – a positive, right-skewed law.

Parameters:

• mu (Any)

• lam (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Gumbel(loc, scale, *, name=None, keys=None)

Gumbel (type-I extreme-value) with loc and scale – for maxima / extremes.

Parameters:

• loc (Any)

• scale (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

SkewNormal(loc, scale, shape, *, name=None, keys=None)

Skew-normal with loc, scale, and shape (skewness; shape=0 recovers the Normal).

Parameters:

• loc (Any)

• scale (Any)

• shape (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Skellam(mu1, mu2, *, name=None, keys=None)

Skellam: the difference of two independent Poisson(mu1) and Poisson(mu2) counts.

Parameters:

• mu1 (Any)

• mu2 (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

LogSeries(p, *, name=None, keys=None)

Logarithmic (log-series) distribution on {1, 2, ...} with parameter p in (0, 1).

Parameters:

• p (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

VonMises(mu, kappa, *, name=None, keys=None)

Von Mises (circular normal) on the angle (-pi, pi] with mean mu and concentration kappa.

Parameters:

• mu (Any)

• kappa (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

GEV(loc, scale, shape, *, name=None, keys=None)

Generalized extreme value with loc, scale, shape (the block-maxima limit law).

Parameters:

• loc (Any)

• scale (Any)

• shape (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Tweedie(mu, phi, *, name=None, keys=None)

Tweedie compound Poisson-Gamma (power p=1.5) with mean mu and dispersion phi.

A positive distribution with an atom at zero – the standard model for non-negative data that is part-zero, part-continuous (insurance claims, rainfall, ecological biomass).

Parameters:

• mu (Any)

• phi (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

GeneralizedGaussian(mu, alpha, beta, *, name=None, keys=None)

Generalized Gaussian (exponential power) with mu, scale alpha, shape beta.

beta=2 is the Normal and beta=1 is the Laplace, so it interpolates between light and heavy tails – a flexible symmetric error model.

Parameters:

• mu (Any)

• alpha (Any)

• beta (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

GeneralizedPareto(scale, shape, loc=0.0, *, name=None, keys=None)

Generalized Pareto (scale, tail shape, threshold loc) – the peaks-over-threshold tail law.

Parameters:

• scale (Any)

• shape (Any)

• loc (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Nakagami(m, omega, *, name=None, keys=None)

Nakagami-m distribution with shape m and spread omega (signal-fading amplitudes).

Parameters:

• m (Any)

• omega (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Rician(nu, sigma, *, name=None, keys=None)

Rician (Rice) distribution with non-centrality nu and scale sigma (signal-plus-noise magnitude).

Parameters:

• nu (Any)

• sigma (Any)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

Dirichlet(alpha, *, dim=None, name=None, keys=None)

Dirichlet over a simplex; used as a prior on Categorical probabilities (VMP). The concentration alpha is also an inferable parameter. Dirichlet(free) learns the concentration by maximum likelihood, inferring the dimension K from the observed simplex data (no dim= needed); pass dim=K to request the explicit positive-K-vector parameter treatment for how='mcmc'|'ensemble'|'map'.

Parameters:

• alpha (Any)

• dim (int | None)

• name (str | None)

• keys (str | None)

Return type:

RandomVariable

class Conv(field='x', *, channels=(64, 128, 256), out=10)

Bases: _NeuralPredictor

A conv-net predictor over image covariates (C, H, W). Conv(channels=[64,128,256], out=10) is a VGG-style stack (two 3x3 convs + BatchNorm + max-pool per channel stage), global-pooled into a linear head. Use it exactly like Net: Categorical(logits=Conv(out=10)).fit(y, given={"x": images}).

Parameters:

• field (Any)

• channels (Any)

• out (int)

build(in_shape)

Build a Torch convolutional module for image-shaped covariates.

Parameters:

in_shape (Any)

Return type:

Any

Embedding

alias of CategoricalEmbedding

Mix(components, weights=None, *, name=None)

Symbolic finite mixture over PPL variables or concrete distributions.

Mix([Normal(free, free), Normal(free, free)]).fit(data) fits a 2-component Gaussian mixture. The fitted object exposes mixture responsibilities through posterior helpers and inherits the numeric stability contract of the stats mixture implementation.

Parameters:

name (str | None)

Return type:

RandomVariable

class Net(field='x', *, hidden=(64,), out=1)

Bases: _NeuralPredictor

An MLP predictor over vector covariates. Net(hidden=[256], out=10) is a one-hidden-layer ReLU net.

Parameters:

• field (Any)

• hidden (Any)

• out (int)

build(in_shape)

Build a Torch MLP module for the inferred covariate shape.

Parameters:

in_shape (Any)

Return type:

Any

class Transformer(field='x', *, out, d_model=128, n_layer=3, n_head=4, embedding=None)

Bases: _NeuralPredictor

A causal decoder-only Transformer predictor over a (block,) context of token ids.

Categorical(logits=Transformer(out=vocab, d_model=256, n_layer=6, n_head=8)) is autoregressive next-token prediction p(token | context) – it lowers to the same SoftmaxNeuralLeaf as Net/Conv (cross-entropy = next-token NLL) and fits through the unchanged estimate() loop with .fit(next_tokens, given={"x": context_windows}). The context width is the block inferred from the data.

Parameters:

• field (Any)

• out (int)

• d_model (int)

• n_layer (int)

• n_head (int)

• embedding (Any)

build(in_shape)

Build a causal language-model module for the inferred token context width.

Parameters:

in_shape (Any)

Return type:

Any

Flow(dim, *, hidden=32, layers=4, m_steps=80, lr=5e-3)

An exact p(x) over R^dim via a RealNVP coupling flow. Fit with .fit(x).

Parameters:

• dim (int)

• hidden (int)

• layers (int)

• m_steps (int)

• lr (float)

Return type:

RandomVariable

MAF(dim, *, hidden=64, blocks=3, m_steps=80, lr=5e-3)

An exact p(x) over R^dim via a masked autoregressive flow (richer dependence). Fit with .fit(x).

Parameters:

• dim (int)

• hidden (int)

• blocks (int)

• m_steps (int)

• lr (float)

Return type:

RandomVariable

VAE(dim, *, latent=2, hidden=32, eval_samples=16, m_steps=120, lr=5e-3)

A latent-variable p(x) via a VAE whose score is a reproducible ELBO estimate; fit with .fit(x).

Parameters:

• dim (int)

• latent (int)

• hidden (int)

• eval_samples (int)

• m_steps (int)

• lr (float)

Return type:

RandomVariable

DiscreteAR(dim, cats, *, hidden=64, m_steps=100, lr=5e-3)

An exact p(x) over discrete vectors x in {0..cats-1}^dim (autoregressive). Fit with .fit(x).

Parameters:

• dim (int)

• cats (int)

• hidden (int)

• m_steps (int)

• lr (float)

Return type:

RandomVariable

EBM(dim, *, hidden=64, layers=3, noise_ratio=2, m_steps=250, lr=5e-3)

An energy-based p(x) ∝ exp(-E(x)) over R^dim, trained by NCE (approximately normalized). Fit .fit(x).

Parameters:

• dim (int)

• hidden (int)

• layers (int)

• noise_ratio (int)

• m_steps (int)

• lr (float)

Return type:

RandomVariable

MDN(x_dim, y_dim, *, k=5, hidden=32, field='x', m_steps=120, lr=5e-3)

A multimodal, heteroscedastic p(y | x) via a mixture density network. Fit .fit(y, given={"x": X}).

Parameters:

• x_dim (int)

• y_dim (int)

• k (int)

• hidden (int)

• field (str)

• m_steps (int)

• lr (float)

Return type:

RandomVariable

CondFlow(x_dim, y_dim, *, hidden=32, layers=4, field='x', m_steps=100, lr=5e-3)

An exact conditional p(y | x) via a conditional coupling flow (needs y_dim >= 2). Fit with covariates.

Parameters:

• x_dim (int)

• y_dim (int)

• hidden (int)

• layers (int)

• field (str)

• m_steps (int)

• lr (float)

Return type:

RandomVariable

CondDiscreteAR(x_dim, y_dim, cats, *, hidden=64, field='x', m_steps=120, lr=5e-3)

An exact conditional p(y | x) over discrete y (autoregressive, conditioned on x). Fit w/ covariates.

Parameters:

• x_dim (int)

• y_dim (int)

• cats (int)

• hidden (int)

• field (str)

• m_steps (int)

• lr (float)

Return type:

RandomVariable

SemiMix(components, weights=None, *, name=None)

Semi-supervised finite mixture over component RandomVariables (or concrete distributions).

Like Mix(), but each observation is a (value, prior) pair where prior is either None (unlabeled) or a sequence of (component_index, probability) pairs giving a partial label. Labeled rows restrict/re-weight the responsibilities to the listed components, so a few labels can anchor the components. SemiMix([Normal(free, free), Normal(free, free)]).fit(data) fits a 2-component Gaussian mixture from a mix of labeled and unlabeled rows.

Parameters:

name (str | None)

Return type:

RandomVariable

Seq(element, *, name=None)

Symbolic iid sequence whose observations are lists or arrays of element.

Parameters:

name (str | None)

Return type:

RandomVariable

Markov(emission, states=None, *, transitions=None, initial=None, name=None)

Hidden Markov model over latent states emitting emission.

Markov(Normal(free, free), states=2).fit(sequences) fits a 2-state Gaussian HMM by EM (emissions k-means++ seeded so states separate); .posterior(sequences) gives state posteriors. For per-state priors pass a list of emissions, one per state: Markov([Normal(m0, 1), Normal(m1, 1)]). The transition matrix and initial distribution are inferable parameters too: pass transitions=free / transitions=Dirichlet(alpha) (each row a simplex) and/or initial=free / initial=Dirichlet(alpha) and fit with how='mcmc'|'ensemble'|'map' (typically with an ordered-emission constraint for identifiability).

Parameters:

• states (int | None)

• name (str | None)

Return type:

RandomVariable

LDA(num_topics, vocab_size, *, alpha=1.0, name=None)

Latent Dirichlet allocation. Fit on a list of documents, each a bag of (word_id, count) pairs over word ids 0..vocab_size-1. Topics are recovered as word distributions; alpha (the document-topic Dirichlet) is fixed by default.

Parameters:

• num_topics (int)

• vocab_size (int)

• alpha (float)

• name (str | None)

Return type:

RandomVariable

MVN(dim, *, mean=None, cov=None, name=None)

Multivariate Gaussian of dimension dim (full covariance). Fit on a list of length-dim vectors; MVN(dim).fit(X) recovers mean and covariance by EM.

The mean vector and covariance matrix are also inferable parameters: pass mean=free (a dim-vector on the real line) or mean=ordered (increasing entries, for identifiability) and/or cov=free (a full SPD covariance via its Cholesky factor) and fit with how='mcmc'|'ensemble'|'map'.

Parameters:

• dim (int)

• name (str | None)

Return type:

RandomVariable

DiagGaussian(dim, *, mean=None, var=None, name=None)

Diagonal-covariance multivariate Gaussian of dimension dim. DiagGaussian(dim).fit(X) recovers mean and per-axis variance by EM; the mean vector (mean=free / ordered) and diagonal variances (var=free, a positive vector) are also inferable parameters via how='mcmc'|'ensemble'|'map'.

Parameters:

• dim (int)

• name (str | None)

Return type:

RandomVariable

LocalLevel(*, name=None)

Local-level state-space model (random walk + noise) for a time series. Fit on a 1-D series; recovers level/observation noise and smoothed states (Kalman/RTS + EM).

Parameters:

name (str | None)

Return type:

RandomVariable

AR1(*, name=None)

AR(1)-plus-noise state-space model; estimates the autoregressive coefficient phi.

Parameters:

name (str | None)

Return type:

RandomVariable

Graph()

A VMP factor graph for arbitrary conjugate-Gaussian DAGs with shared variables. See mixle.ppl.vmp.Graph.

class Field(name)

Bases: object

A named covariate (data column) for regression: a * Field("x") + b.

Parameters:

name (str)

class Group(name, slopes=())

Bases: object

A by-group random-effects term for mixed-effects models. Group("subject") is a random intercept (lme4’s (1|subject)); Group("subject", slopes=["x"]) adds a correlated random slope on x ((1 + x | subject)).

Parameters:

name (str)

compare(models, data, *, by='aic')

Compare fitted models on data. Returns rows sorted best-first by by (‘aic’ | ‘bic’ | ‘loglik’ | ‘waic’ | ‘loo’).

'waic' and 'loo' are the Bayesian predictive criteria (integrating over parameter uncertainty via the posterior draws of a Bayesian fit); 'aic'/'bic' use the point estimate. Each row also reports elpd differences from the best model (d_elpd) for waic/loo.

Parameters:

by (str)

posterior_predictive_check(fitted, data, *, statistics=None, n_rep=1000, seed=0)

Posterior predictive check of a fitted PPL model against data.

Draws n_rep replicate datasets (each the size of data) from fitted.predict – which integrates over parameter uncertainty for a Bayesian fit (conjugate/mcmc/hmc) and is the plug-in predictive for a point fit (em/map) – evaluates each named statistic on every replicate and on the observed data, and returns the Bayesian p-value per statistic.

Returns {'observed', 'replicated', 'p_value', 'n_rep'}: observed[name] the statistic on the data, replicated[name] its (n_rep,) replicate values, p_value[name] = P(T_rep >= T_obs).

Parameters:

• fitted (RandomVariable)

• data (Sequence[Any])

• statistics (dict[str, Callable[[ndarray], float]] | None)

• n_rep (int)

• seed (int)

Return type:

dict[str, Any]

prior_predictive(model, size, *, n_rep=1000, statistics=None, seed=0)

Prior predictive simulation: n_rep datasets of size drawn from model’s prior.

For each replicate it draws every prior parameter (and hyperparameter) and then size data points, so the result reflects what the model believes before seeing data – the check that a prior is neither absurdly tight nor absurdly diffuse. Returns {'replicated': {stat: (n_rep,)}, 'samples': (n_rep, size), 'n_rep'} with the per-replicate statistics and the raw simulated datasets.

Parameters:

• model (RandomVariable)

• size (int)

• n_rep (int)

• statistics (dict[str, Callable[[ndarray], float]] | None)

• seed (int)

Return type:

dict[str, Any]

prior_predictive_check(model, data, *, statistics=None, n_rep=1000, seed=0)

Prior predictive check: where the observed statistics sit in the prior predictive distribution.

Like posterior_predictive_check() but the replicates come from the prior (via prior_predictive()), so a p-value near 0 or 1 flags a prior that is inconsistent with the data before any fitting – often a sign the prior is mis-scaled.

Parameters:

• model (RandomVariable)

• data (Sequence[Any])

• statistics (dict[str, Callable[[ndarray], float]] | None)

• n_rep (int)

• seed (int)

Return type:

dict[str, Any]

fit_censored(model, time, *, event=None, lower=None, upper=None, seed=0, max_its=4000, n_starts=2)

Fit a distribution’s free parameters to right-censored and/or truncated data by ML.

model is a flat PPL distribution with free parameter slots, e.g. Weibull(free, free) or Exponential(free). time are the (possibly censored) values; event flags which are observed events vs right-censored (default all observed); lower/upper mark truncation of the sampling window. Maximizes censored_loglik() over the free slots (Nelder-Mead in the unconstrained space, respecting each slot’s positivity/unit support) and returns the fitted model as a bound RandomVariable (with .summary()).

Parameters:

• model (RandomVariable)

• time (Sequence[float])

• event (Sequence[Any] | None)

• lower (float | None)

• upper (float | None)

• seed (int)

• max_its (int)

• n_starts (int)

Return type:

RandomVariable

fit_with_provenance(rv, data, *, seed=None, **fit_kw)

Fit a PPL random variable and return an explicit provenance result.

One pre-fit snapshot (see _replayable_data()) is used for both the fit and the hash, and the snapshot is fingerprinted before the fit runs and again when the header is built. Those two fingerprints must agree: if the fit edits its input in place the header would otherwise report a hash of records the fit never saw, and a provenance header that silently describes the wrong dataset is worse than one that refuses to be built. The extra hash is one SHA-256 pass over the snapshot and no extra memory. Mutation inside a record object is deliberately still detected by this check but not prevented – the snapshot is shallow (MXR-080-1897).

seed=None has the deterministic effective seed 0; otherwise the exact non-negative integer seed initializes the RandomState passed to rv.fit. Supplying a separate rng is rejected so the recorded state cannot disagree with the state used.

Parameters:

• rv (Any)

• data (Any)

• seed (int | None)

• fit_kw (Any)

Return type:

PPLFitResult

class PPLFitResult(fitted, header)

Bases: NamedTuple

Explicit fitted-object and provenance pair.

A named tuple preserves the historical fitted, header = ... unpacking contract while avoiding optional, failure-prone mutation of fitted.

Parameters:

• fitted (Any)

• header (Any)

fitted: Any

Alias for field number 0

header: Any

Alias for field number 1

class PointwiseLogLikelihood(*args, **kwargs)

Bases: Protocol

A result that can score observations under each posterior draw.

pointwise_log_likelihood(data)

Return a (n_draws, n_observations) log-likelihood matrix.

Parameters:

data (Any)

Return type:

ndarray

class Predictive(*args, **kwargs)

Bases: Protocol

A result carrying an executable posterior-predictive callable.

class ResultCapabilities(summarizable, sampleable, predictive, pointwise_log_likelihood)

Bases: object

Normalized runtime capability view for any public fit result.

Parameters:

• summarizable (bool)

• sampleable (bool)

• predictive (bool)

• pointwise_log_likelihood (bool)

class Sampleable(*args, **kwargs)

Bases: Protocol

A result that can return parameter / latent samples(...) (read by RandomVariable.posterior).

supports(r, Sampleable) is the runtime capability check.

samples(param=..., *args, **kwargs)

Return posterior samples for a parameter, latent, or default result target.

Parameters:

• param (Any)

• args (Any)

• kwargs (Any)

Return type:

ndarray

class Summarizable(*args, **kwargs)

Bases: Protocol

A result that can report a posterior / fit summary() (read by RandomVariable.summary).

supports(r, Summarizable) is the runtime capability check.

summary()

Return a dictionary summary of the fitted result.

Return type:

dict

result_capabilities(result)

Return capability flags without assuming a particular fitter result class.

Parameters:

result (Any)

Return type:

ResultCapabilities

censored_loglik(dist, time, *, event=None, lower=None, upper=None)

Total log-likelihood of right-censored and/or truncated time under a fitted dist.

event[i] true (default all true) means time[i] is an observed event contributing log f(time[i]); false means it is right-censored, contributing the log-survival log(1 - F(time[i])). lower/upper truncate the support: every point then also subtracts log(F(upper) - F(lower)) (use None for an open end). Requires dist.cdf.

Parameters:

time (Sequence[float])

Return type:

float

kaplan_meier(time, event=None)

Kaplan-Meier nonparametric survival estimate S(t) from right-censored data.

Returns {'time', 'survival', 'at_risk', 'events'} over the distinct event times – the standard model-free survival curve to plot against, or compare a fitted parametric model to.

Parameters:

• time (Sequence[float])

• event (Sequence[Any] | None)

Return type:

dict[str, ndarray]

hdi(samples, prob=0.94)

Highest-density interval: the narrowest interval containing prob of the posterior mass.

For a unimodal posterior this is the shortest (low, high) such that P(low <= x <= high) = prob; unlike an equal-tailed interval it tracks an asymmetric or bounded posterior correctly.

Parameters:

• samples (Sequence[float])

• prob (float)

Return type:

tuple[float, float]

posterior_summary(fitted, *, hdi_prob=0.94)

Per-parameter posterior summary table for a fitted PPL model (best after how='mcmc').

Returns a fixed per-parameter schema with mean, sd, hdi_low, hdi_high, ess, ess_tail, r_hat, mcse, diagnostic_status, and diagnostic_error. mean/sd come from the fit’s own summary; the HDI is computed from the posterior draws (when the fit exposes them); ess (effective sample size) and r_hat (Gelman-Rubin, multi-chain) come from the sampler’s diagnostics when present; mcse is the Monte Carlo standard error of the mean, sd / sqrt(ess) – the parameter number’s own noise floor, published WITH the number instead of dropped (STAT-RR17-11). A point fit (em/map) yields just mean/sd.

diagnostic_status semantics (STAT-RR17-11 – “ok” used to mean only “the ESS call did not raise”, so a one-chain fit with r_hat=None or NaN split-R-hat still published parameter numbers under “ok”):

• "ok" – ess, ess_tail, mcse are finite, a finite multi-chain r_hat exists, AND the values clear the standard convergence thresholds: split R-hat <= 1.01, bulk and tail ESS >= 100 (STAT-RR22-13 – finiteness alone labeled R-hat-1.363 fits “ok”). The ONLY promotable state.

• "unconverged-by-diagnostics" – finite diagnostics that FAIL a threshold; the error names the failing number. Not promotable: run longer or with more chains.

• "divergent-transitions" – the fit recorded post-warmup NUTS divergences (STAT-RR23-12: 15- and 16-divergence funnel fits read ok on clean R-hat/ESS). Not promotable: divergences are direct evidence of unexplored posterior geometry; reparameterize or raise target_accept.

• "single-chain-mixing-unassessable" – ESS computed, but one chain cannot assess mixing (no R-hat); run >= 2 chains before treating the summary as converged evidence.

• "unusable" – a diagnostic evaluated non-finite (NaN/inf ESS or R-hat).

• "unavailable" / "failed" – draws or diagnostics missing / raised.

Parameters:

• fitted (RandomVariable)

• hdi_prob (float)

Return type:

dict[str, dict[str, Any]]

constrain(*constraints)

A joint random variable formed by conditioning several RVs on a relation among them.

constrain(a < b) is the pair (a, b) restricted to a < b; pass several constraints (or combine with & | ~) for richer regions, e.g. constrain(a < b, b < c) orders three variables. The result samples by joint rejection and answers .sample(n) (an (n, k) array, columns in .columns order), .mean()/.var() (per-variable), .prob() (probability the relation holds), and .log_prob(x) (renormalized joint density of the independent variables on the region).

Return type:

RandomVariable

class Constraint(leaves, pred, desc, residual=None, soft=False, reduction='all', margin=None)

Bases: object

A boolean relation over one or more random variables.

Produced by comparisons on RVs — x > 0 (RV vs constant), a < b (RV vs RV), or 2 * a - b >= 1 (linear/transformed expressions on either side) — and combined with & (and), | (or), ~ (not). A constraint over a single RV is consumed by rv.given(c) (truncation); a constraint over several RVs is consumed by constrain(c) (joint conditioning) or fit(..., constraints=c) (feasible region).

leaves are the distinct leaf RVs the relation depends on; pred(env) evaluates the relation given env, a dict mapping each leaf RV to its value(s).

property rv

The single RV this constraint restricts (back-compat for one-variable events).

eval(env)

Evaluate the constraint predicate against an environment of RV values.

eval_rows(env, rows=None)

Evaluate and reduce to exactly one boolean per sampled row.

Parameters:

rows (int | None)

Return type:

ndarray

with_reduction(reduction)

Return the same event with vector entries reduced by all or any per row.

Parameters:

reduction (str)

Return type:

Constraint

contains(x)

Evaluate a single-variable constraint directly on that variable’s value(s).

Event

alias of Constraint

eq(lhs, rhs)

Build an equality relation lhs == rhs over RVs/expressions/constants.

== is not overloaded on RandomVariable (RVs are used as dict keys by identity), so build equalities with this function or rv.eq(...). Equalities have measure zero and cannot be honored by rejection, so consume them with the soft-penalty inference path, e.g. model.fit(data, constraints=eq(a + b, 1.0), penalty=100.0).

Return type:

Constraint

equal(lhs, rhs)

Build an equality relation lhs == rhs over RVs/expressions/constants.

== is not overloaded on RandomVariable (RVs are used as dict keys by identity), so build equalities with this function or rv.eq(...). Equalities have measure zero and cannot be honored by rejection, so consume them with the soft-penalty inference path, e.g. model.fit(data, constraints=eq(a + b, 1.0), penalty=100.0).

Return type:

Constraint

ne(lhs, rhs)

Build an inequality relation lhs != rhs (boolean only; no smooth penalty surface).

Return type:

Constraint

potential(fn, *vars, name=None)

Add a custom log-factor fn(*values) to a model’s joint log-density.

vars are random variables: parameters of the model (named priors, or param(...) vector/matrix handles, exactly the references a constrain()/eq() constraint may use), or named prior RVs that are not parameters of the model, which a potential introduces as auxiliary latents: each one becomes an extra inference slot whose own prior joins the joint target (a constraint on such an RV is rejected; a potential is how a latent that no distribution slot reads enters the model). At each inference evaluation they are resolved to their current values and passed to fn positionally; fn returns a scalar log-weight that is added to log p(data, theta). Use it for anything the distribution slots can’t say directly: a soft coupling between two latents, a bespoke log-prior, a penalty/regularizer:

a = Normal(0, 10, name="a"); b = Normal(0, 10, name="b")
m = Normal(a, 1.0).fit(data, potentials=potential(lambda av, bv: -0.5 * (av - bv) ** 2, a, b))

Pass one potential or a list via fit(..., potentials=...). Potentials route inference through the numerical target (how in map / mcmc / hmc / nuts / ensemble; auto picks map); every referenced variable must be a named prior RV or a param(...) handle.

Return type:

_Potential

increasing(v, *, strict=False)

The entries of a vector RV/expression are non-decreasing (strict -> strictly increasing).

Parameters:

strict (bool)

Return type:

Constraint

decreasing(v, *, strict=False)

The entries of a vector RV/expression are non-increasing (strict -> strictly decreasing).

Parameters:

strict (bool)

Return type:

Constraint

monotone(v)

The entries are monotone — non-decreasing or non-increasing.

Return type:

Constraint

convex(v)

The entries are convex: the second difference is non-negative everywhere.

Return type:

Constraint

concave(v)

The entries are concave: the second difference is non-positive everywhere.

Return type:

Constraint

lipschitz(v, bound)

Bounded first difference: |v[i+1] - v[i]| <= bound (a discrete smoothness constraint).

Parameters:

bound (float)

Return type:

Constraint

ode_residual(v, f, dt=1.0, *, tol=1e-2)

A differential-equation constraint: v (a function sampled on a uniform grid of spacing dt) satisfies dv/dt = f(v). The signed residual diff(v)/dt - f(v[:-1]) feeds the soft-penalty inference path — fit(..., constraints=ode_residual(y, f), penalty=w) fits a physics-informed curve. Like an equality it is measure-zero, so consume it with penalty= rather than by rejection.

Parameters:

• dt (float)

• tol (float)

Return type:

Constraint

class GaussianField(index, kernel, name='field')

Bases: object

A latent field: an index grid plus a Gaussian (GP/GMRF) prior over its node values.

Parameters:

• index (ndarray)

• kernel (FieldKernel)

• name (str)

class FieldSystem(fields, coregion=None)

Bases: object

Several named latent fields fit jointly – the spine for coupled multiphysics / multivariate earth-science models (e.g. an ore-grade field constrained by several geophysical surveys, or coupled paleo-environment fields like temperature + salinity + pCO2 read through many proxies).

By default the prior is block-diagonal over the fields (each keeps its own kernel precision) and cross-field dependence is expressed in the proxy likelihoods: a proxy loglik(field_t, params, torch) receives the full params dict, so it can read any field by name and couple them (a coupling PDE residual, a grade that depends on density + susceptibility, …). Attach a proxy to a particular field with proxy.on('name'); an unattached proxy defaults to the first field.

Pass coregion=B (a K x K between-field covariance for the K fields) for prior-level coregionalization – the intrinsic coregionalization model, joint prior precision B^-1 (x) Lambda, where the fields share one spatial structure Lambda (the first field’s kernel) and are correlated a priori through B (e.g. temperature and salinity covary before any data). Requires all fields to share one index/dim. B must be symmetric positive-definite.

Parameters:

• fields (Sequence[GaussianField])

• coregion (ndarray | None)

class FieldKernel

Bases: object

A Gaussian field prior over an index grid, expressed as a precision matrix.

precision(index)

Return the prior precision matrix for the supplied field index grid.

Parameters:

index (ndarray)

Return type:

ndarray

covariance(index)

The prior covariance matrix, when available directly (a kernel defined by its covariance). Used by the low-rank Gauss-Newton posterior to get field marginals without a dense precision inverse; None (the default) falls back to the dense path.

Parameters:

index (ndarray)

Return type:

ndarray | None

class GreatCircleRBF(lengthscale=1.0, amplitude=1.0, jitter=1e-06, radius=1.0)

Bases: FieldKernel

Squared-exponential GP prior on the sphere for a (lat, lon) index (degrees).

K[i,j] = amplitude^2 exp(-0.5 chord_ij^2 / lengthscale^2) where chord is the straight-line distance through the R^3 embedding – so the kernel is positive-definite on the sphere at every lengthscale (a geodesic-distance RBF is not). For nearby points the chord matches the great-circle distance, so lengthscale reads as a spatial correlation length in the same units as radius (radius=6371.0088 -> kilometres on Earth; default unit sphere -> radians of arc).

Parameters:

• lengthscale (float)

• amplitude (float)

• jitter (float)

• radius (float)

covariance(index)

Return the chordal-distance RBF covariance on spherical coordinates.

Parameters:

index (ndarray)

Return type:

ndarray

precision(index)

Return the inverse of the spherical RBF covariance matrix.

Parameters:

index (ndarray)

Return type:

ndarray

class GreatCircleMatern(lengthscale=1.0, amplitude=1.0, nu=1.5, jitter=1e-06, radius=1.0)

Bases: FieldKernel

Matern GP prior on the sphere for a (lat, lon) index (degrees) – the geostatistics default.

A Matern covariance of the R^3 chordal distance (so PD on the sphere at every lengthscale), giving rougher, more realistic spatial fields than the infinitely-smooth RBF. nu in {0.5, 1.5, 2.5} (0.5 = exponential / Ornstein-Uhlenbeck, 1.5 and 2.5 are the common once/twice-differentiable choices). lengthscale and radius carry the same spatial-units meaning as in GreatCircleRBF.

Parameters:

• lengthscale (float)

• amplitude (float)

• nu (float)

• jitter (float)

• radius (float)

covariance(index)

Return the chordal-distance Matern covariance on spherical coordinates.

Parameters:

index (ndarray)

Return type:

ndarray

precision(index)

Return the inverse of the spherical Matern covariance matrix.

Parameters:

index (ndarray)

Return type:

ndarray

great_circle_distance(latlon_a, latlon_b=None, *, radius=1.0)

Great-circle (geodesic) distance between (lat, lon) points in degrees.

Returns radius * theta where theta is the central angle (radians). With radius=1 the result is the angle itself; radius=6371.0088 gives kilometres on Earth. latlon_b=None returns the full pairwise (n, n) matrix for latlon_a; otherwise the (na, nb) cross matrix (a float when both are single points). Uses the haversine formula for accuracy at small angles.

Parameters:

• latlon_a (ndarray)

• latlon_b (ndarray | None)

• radius (float)

class RandomWalk(scale=1.0, order=1, ridge=None)

Bases: FieldKernel

A Gaussian-Markov-random-field smoothness prior via finite differences.

order=1 penalizes sum (f[i+1]-f[i])^2 / scale^2 (a random-walk / integrated-noise prior); order=2 penalizes the discrete curvature (an integrated-Wiener / thin-plate prior). ridge adds I / ridge^2 so the otherwise-improper prior is proper (it anchors the level/trend).

Parameters:

• scale (float)

• order (int)

• ridge (float | None)

precision(index)

Build a spacing-aware finite-difference precision on a 1-D physical grid.

Parameters:

index (ndarray)

Return type:

ndarray

class RBF(lengthscale=1.0, amplitude=1.0, jitter=1e-06)

Bases: FieldKernel

A squared-exponential (RBF) GP prior: K[i,j] = amplitude^2 exp(-0.5 (d_ij/lengthscale)^2).

The precision is inv(K + jitter*I). index may be 1-D (a grid) or 2-D (coordinates, one row per node) – distances are Euclidean, so this is the spatial-field prior as well.

Parameters:

• lengthscale (float)

• amplitude (float)

• jitter (float)

covariance(index)

Return the RBF covariance matrix plus diagonal jitter for index.

Parameters:

index (ndarray)

Return type:

ndarray

precision(index)

Return the inverse of the jittered RBF covariance matrix.

Parameters:

index (ndarray)

Return type:

ndarray

class AnisotropicRBF(ranges=(1.0, 1.0), angle=0.0, amplitude=1.0, jitter=1e-06, metric=None)

Bases: FieldKernel

A geometrically-anisotropic squared-exponential GP prior – correlation that is longer along one direction than another (geological layering / bedding, faulted fabric, flow channels).

The squared distance is the Mahalanobis form (x-y)^T M (x-y) with M built from per-axis correlation ranges after rotating the coordinates: in 2-D, angle (radians, counter-clockwise from the x-axis) sets the principal direction and ranges=(major, minor) the correlation length along/across it. For >2-D, supply ranges per axis (axis-aligned) or a full metric matrix M directly. Provably positive-definite – it is an ordinary RBF on linearly-transformed coordinates.

Parameters:

• ranges (Sequence[float])

• angle (float)

• amplitude (float)

• jitter (float)

• metric (ndarray | None)

covariance(index)

Return the anisotropic RBF covariance matrix plus diagonal jitter.

Parameters:

index (ndarray)

Return type:

ndarray

precision(index)

Return the inverse of the anisotropic RBF covariance matrix.

Parameters:

index (ndarray)

Return type:

ndarray

class Proxy

Bases: object

A likelihood hung off the shared field. Subclasses declare params and score the data in torch.

on(field_name)

Attach this proxy to a named field of a FieldSystem; returns self for chaining.

Parameters:

field_name (str)

Return type:

Proxy

validate(field_dim)

Validate observation shapes against the attached field dimension.

Parameters:

field_dim (int)

Return type:

None

params()

Return free parameter specifications used by the joint field fit.

Return type:

list[_ParamSpec]

loglik(field_t, params, torch)

Return this proxy’s torch log-likelihood contribution.

Parameters:

• field_t (Any)

• params (dict)

Return type:

Any

residual(field_t, params, torch)

The standardized Gaussian residual (y - prediction) / scale (1-D), or None if this proxy is not a Gaussian-misfit observation. Used by the Gauss-Newton posterior (how='gauss_newton').

Parameters:

• field_t (Any)

• params (dict)

Return type:

Any

class GaussianProxy(y, index=None, slope=1.0, intercept=0.0, scale=1.0, prefix='gauss')

Bases: Proxy

A linear-Gaussian forward model: y_j ~ N(intercept + slope * field[index_j], scale).

slope, intercept and scale are each a fixed float or the free token (estimated). With everything fixed this is the exact linear-Gaussian observation that makes the joint posterior exactly Gaussian – the closed-form check on the Laplace covariance.

Parameters:

• y (ndarray)

• index (ndarray | None)

• slope (Any)

• intercept (Any)

• scale (Any)

• prefix (str)

validate(field_dim)

Validate observation shapes against the attached field dimension.

Parameters:

field_dim (int)

Return type:

None

params()

Return parameter specs for any free slope, intercept, or scale.

Return type:

list[_ParamSpec]

loglik(field_t, params, torch)

Return the linear-Gaussian observation log-likelihood.

residual(field_t, params, torch)

Return standardized Gaussian residuals for Gauss-Newton covariance estimation.

class LogisticNicheProxy(presence, mu_scale=2.0, prefix='niche')

Bases: Proxy

Thermal-niche occupancy: presence[i,t] ~ Bernoulli(sigmoid(b - 0.5 kappa_i (field[t] - mu_i)^2)).

Each row of presence is one taxon’s presence/absence across the field’s index bins; the latent niche location mu_i and precision kappa_i are co-estimated, and b is a shared baseline. The assemblage acts as a community thermometer – a unimodal response peaked at the niche optimum.

Parameters:

• presence (ndarray)

• mu_scale (float)

• prefix (str)

validate(field_dim)

Validate observation shapes against the attached field dimension.

Parameters:

field_dim (int)

Return type:

None

params()

Return niche-location, niche-precision, and baseline parameter specifications.

Return type:

list[_ParamSpec]

loglik(field_t, params, torch)

Return the Bernoulli niche-response log-likelihood with weak location regularization.

class PoissonProxy(counts, index=None, offset=0.0, prefix='cox')

Bases: Proxy

A log-Gaussian Cox process: counts[j] ~ Poisson(exp(offset_j + field[index_j])).

The field is the latent log-intensity; counts are the point/aggregate observations. This is the canonical Cox-process proxy – the namesake of the foundation.

Parameters:

• counts (ndarray)

• index (ndarray | None)

• offset (Any)

• prefix (str)

validate(field_dim)

Validate observation shapes against the attached field dimension.

Parameters:

field_dim (int)

Return type:

None

loglik(field_t, params, torch)

Return the Cox-process Poisson log-likelihood up to the count-factorial constant.

class CustomProxy(loglik_fn, param_specs=(), prefix='custom')

Bases: Proxy

An arbitrary proxy: supply a torch log-likelihood loglik_fn(field_t, params, torch) and the parameter specs it reads from params (each (name, shape, support, init)).

Parameters:

• loglik_fn (Callable)

• param_specs (Sequence[tuple])

• prefix (str)

params()

Return the custom parameter specifications declared by param_specs.

Return type:

list[_ParamSpec]

loglik(field_t, params, torch)

Delegate log-likelihood evaluation to loglik_fn.

fit_field(field, proxies, *, how='laplace', max_iter=1000, lr=0.4, init=None, vi_steps=400, vi_lr=0.05, vi_samples=4, rng=None, tolerance_grad=1.0e-6, tolerance_change=1.0e-9)

Fit a latent field jointly to a list of proxy likelihoods.

field=None runs pure-parameter inference (the proxies carry all the latents, e.g. ODE/PDE coefficients) with no shared field. how='map' returns the joint MAP (no covariance); how='laplace' adds the Gaussian posterior (the inverse-Hessian covariance, exact when every factor is Gaussian; needs a twice-differentiable forward); how='gauss_newton' builds the posterior from J^T J + prior with J the Jacobian of the standardized residual – first-order only, so it works with the sparse adjoint solve (where how='laplace' cannot run) and is exact only for a linear Gaussian forward with fixed scale; how='vi' fits a mean-field Gaussian variational posterior (reparameterized ADVI on the unconstrained vector), an explicit approximation for genuinely non-Gaussian posteriors (e.g. with total-variation / Potts priors), and it too is first-order so it works through the sparse solve. Every result reports convergence and uncertainty-approximation metadata; posteriors over any node are read off the returned FieldPosterior.

Parameters:

• field (GaussianField | None)

• proxies (Sequence[Proxy])

• how (str)

• max_iter (int)

• lr (float)

• init (dict | None)

• vi_steps (int)

• vi_lr (float)

• vi_samples (int)

• tolerance_grad (float)

• tolerance_change (float)

Return type:

FieldPosterior

class FieldPosterior(map_values, _cov, _layout, _field_name, _hessian, objective, _field_prior=<factory>, _proxy_info=<factory>, _field_coordinate_count=0, _supports=<factory>, _marg_var=<factory>, _variational_params=<factory>, _parameter_ids=<factory>, initialization_space='natural', converged=True, iterations=0, termination_reason='not-reported', curvature_method='none', curvature_rank=None, curvature_dim=None, regularization=0.0, uncertainty_status='unavailable', approximation='none', multistart_receipt=())

Bases: object

The joint posterior. posterior(node) returns (mean, sd) for the field or any proxy param.

The Laplace covariance is the inverse of the joint negative-log-posterior Hessian at the MAP – the prior precision plus every proxy’s observed local curvature. posterior(field, coupling=False) instead inverts only the field’s own Hessian block (the posterior conditional on the other nodes at their MAP). The subset diagnostic in field_posterior() retains all cross-field and nuisance-parameter curvature contributed by the selected proxies.

Parameters:

• map_values (dict)

• _cov (ndarray)

• _layout (dict)

• _field_name (str)

• _hessian (ndarray)

• objective (float)

• _field_prior (ndarray)

• _proxy_info (dict)

• _field_coordinate_count (int)

• _supports (dict)

• _marg_var (dict)

• _variational_params (dict)

• _parameter_ids (dict)

• initialization_space (str)

• converged (bool)

• iterations (int)

• termination_reason (str)

• curvature_method (str)

• curvature_rank (int | None)

• curvature_dim (int | None)

• regularization (float)

• uncertainty_status (str)

• approximation (str)

• multistart_receipt (tuple)

mean(node)

Return the posterior mean/MAP value for node in its natural parameter space.

Parameters:

node (str)

Return type:

ndarray

parameter_id(node)

Return the stable optimization-layout ID for a public node alias.

Parameters:

node (str)

Return type:

str

cov(node)

Posterior covariance of node in its natural space (delta method applied for positive nodes).

Parameters:

node (str)

Return type:

ndarray

sd(node)

Return posterior standard deviations for node in its natural parameter space.

Parameters:

node (str)

Return type:

ndarray

posterior(node, *, coupling=True)

(mean, sd) for node in its natural space. coupling=True (default) marginalizes the other nodes; coupling=False fixes them at the MAP for an additive-information diagnostic.

Parameters:

• node (str)

• coupling (bool)

Return type:

tuple[ndarray, ndarray]

field_posterior(include=None)

The field posterior (mean, sd) under a subset of proxies, evaluated at the one joint MAP.

Because observed curvature is additive, the joint local precision under a subset of proxies is the prior curvature plus those proxies’ full observed-curvature blocks. include=None uses every proxy (the joint posterior); include=["gauss"] uses only that proxy – so you can read how much each proxy sharpens the shared field without re-fitting (which would be ill-posed for a single proxy with a free forward-model gain). Cross-field and nuisance-parameter blocks are retained before the primary-field marginal is extracted. This remains a local-curvature ablation at the joint MAP, not a separately normalized posterior. The mean is the joint MAP for every subset.

Parameters:

include (Sequence[str] | None)

Return type:

tuple[ndarray, ndarray]

sample(size=1, rng=None, *, nodes=None, given=None)

Draw joint samples from the Gaussian posterior, returned in each node’s natural space.

The Laplace / Gauss-Newton posterior is the Gaussian N(map, _cov) in the unconstrained parameter space; this draws via a Cholesky factor of that joint covariance (so cross-node and within-field correlations are preserved) and maps each draw back through the node’s support transform (exp for positive nodes). That back-transform is exact – unlike cov(), which linearizes it with the delta method – so a positive node’s draws are properly lognormal. A mean-field how='vi' (or low-rank Woodbury) fit exposes only per-node marginal variances, so there the draws are independent across nodes (which is exactly the mean-field assumption).

given (a {node: value} dict) conditions the draw on fixed node values – the closed-form Gaussian conditional of the remaining coordinates given the observed ones (e.g. draw the field consistent with a pinned measurement or a fixed proxy parameter). Conditioning couples nodes through the joint covariance, so it requires the full-covariance fit (how='laplace'/'gauss_newton'); the fixed nodes retain their given values.

Parameters:

• size (int) – number of joint draws.

• rng – a numpy.random.RandomState, an integer seed, or None.

• nodes (Sequence[str] | None) – which nodes to return (default: all). The draw is always joint over the full vector.

• given (dict | None) – optional {node: value} to condition on (values in each node’s natural space).

Returns:

{node: ndarray} with shape (size,) for scalar nodes and (size, dim) otherwise.

Return type:

dict

summary()

Return posterior means and standard deviations for every non-empty node.

Return type:

dict

class GP(name, index, kernel)

Bases: object

A latent field node for the equation-style surface: T = GP("T", index=grid, kernel=...).

Supports affine algebra (c0 - c1*T, T + b) so a linear forward model reads as math; the result carries the field, the gain and the offset to joint().

Parameters:

• name (str)

• index (np.ndarray)

• kernel (FieldKernel)

Gaussian(y, *, mean, sd)

A linear-Gaussian observation y ~ N(gain*field + offset, sd); mean is an affine in a GP.

This is an observation helper for joint() – it returns a (field, proxy) pair, not a distribution. For the Gaussian distribution as a random variable (the stats Gaussian analogue), use mixle.ppl.Normal(). Also exported as GaussianObs to keep the distinction visible at the call site.

Return type:

tuple

GaussianObs(y, *, mean, sd)

A linear-Gaussian observation y ~ N(gain*field + offset, sd); mean is an affine in a GP.

This is an observation helper for joint() – it returns a (field, proxy) pair, not a distribution. For the Gaussian distribution as a random variable (the stats Gaussian analogue), use mixle.ppl.Normal(). Also exported as GaussianObs to keep the distinction visible at the call site.

Return type:

tuple

Niche(presence, *, over, mu_scale=2.0)

Logistic thermal-niche occupancy of presence over the field over (a community thermometer).

Parameters:

• over (GP)

• mu_scale (float)

Return type:

tuple

Cox(counts, *, log_intensity, offset=0.0)

A log-Gaussian Cox process counts ~ Poisson(exp(offset + field)); log_intensity is a GP.

Return type:

tuple

joint(observations)

Assemble a latent-field model from equation-style observations; call .fit(how=...) to fit.

Each item is the (field, proxy) pair returned by Gaussian(), Niche(), Cox() or mixle.ppl.Differential(). Observations sharing a field must name the same one (this surface targets one shared field); field-free observations (a pure-parameter ODE inverse problem) are allowed.

Parameters:

observations (Sequence[tuple])

Return type:

FieldModel

class FieldModel(field, proxies)

Bases: object

A latent-field model built from equation-style observations; fit with .fit(how=...).

The one fit verb, matching the rest of mixle.ppl: joint([...]).fit(how='map'|'laplace') returns a FieldPosterior with a posterior over any node. Delegates to fit_field().

Parameters:

• field (GaussianField | None)

• proxies (list)

fit(*, how='laplace', **kw)

Fit the assembled field/proxy model with fit_field() and return its posterior.

Parameters:

how (str)

Return type:

FieldPosterior

TotalVariation(over, shape, *, weight=1.0, eps=1e-3)

A smoothed total-variation prior on the field over a structured shape grid: weight * sum over neighbour pairs sqrt((f_a - f_b)^2 + eps^2). Edge-preserving (it does not penalize a jump as harshly as the squared GMRF prior). Returns the (field, proxy) pair for joint().

Parameters:

• weight (float)

• eps (float)

Return type:

tuple

Potts(over, levels, *, weight=1.0)

A discrete-composition prior: weight * sum_i prod_k (f_i - level_k)^2 – a multi-well potential whose minima are the given levels, pulling the field toward a few discrete material values (a smooth relaxation of the Potts model). Combine with TotalVariation() for piecewise-constant regions. Returns the (field, proxy) pair for joint().

Parameters:

weight (float)

Return type:

tuple

multistart(model, inits, *, how='map', **kw)

Fit model from several initializations and keep the best (lowest-objective) fit.

For the multimodal posteriors of nonlinear inverse problems (e.g. scattering / FWI cycle-skipping), a single optimization can land in a poor local mode; inits is a list of {node: value} dicts and the fit with the smallest objective is returned. (Frequency continuation – fit a coarse/low- frequency model, then seed a finer one via fit(init=...) – is the complementary strategy.)

Parameters:

• model (FieldModel)

• inits (Sequence[dict])

• how (str)

Return type:

FieldPosterior

conformal(result, y_cal, *, given, alpha=0.1)

Split-conformal calibration of a fitted regression result into prediction intervals.

Mirrors fit’s convention (labels positional, given= keyword), a one-liner over ConformalRegressor:

m = Normal(free * Field("x") + free, free).fit(y_tr, given={"x": x_tr})
cp = conformal(m.result, y_cal, given={"x": x_cal}, alpha=0.1)
lo, hi = cp.interval({"x": x_te})
cp.covers(y_te, given={"x": x_te}).mean()   # ~ 0.9

Parameters:

• result (Any)

• y_cal (Any)

• given (dict)

• alpha (float)

Return type:

ConformalRegressor

waic(loglik)

Return the WAIC of a (n_draws, n_obs) pointwise log-likelihood matrix.

Parameters:

loglik (ndarray)

Return type:

dict

loo(loglik)

PSIS-LOO under its conventional short name (see psis_loo()).

Parameters:

loglik (ndarray)

Return type:

dict

psis_loo(loglik)

Return PSIS-LOO of a (n_draws, n_obs) pointwise log-likelihood matrix.

Parameters:

loglik (ndarray)

Return type:

dict

loo_stack(logliks, *, iters=2000, tol=1.0e-10)

Stack K candidate models by LOO predictive performance.

logliks is a sequence of (n_draws_k, n_obs) pointwise log-likelihood matrices over the same, aligned observations. Returns the stacking weights, the (n_obs, K) per-model pointwise LOO densities, each model’s elpd_loo, the actually achieved stacked_elpd_loo, and an optimizer receipt. A finite iteration limit need not reach the exact optimum.

Parameters:

• logliks (Sequence[ndarray])

• iters (int)

• tol (float)

Return type:

dict

loo_stacking_weights(pointwise_lpd, iters=2000, tol=1.0e-10, *, return_result=False)

Return LOO stacking weights (Yao, Vehtari, Simpson & Gelman, 2018).

pointwise_lpd is an (n_obs, K) matrix of per-model pointwise LOO log-predictive densities (each column is psis_loo(model_k)["pointwise"]). The returned simplex weights w maximize the LOO log-score of the weighted predictive distribution, sum_i log(sum_k w_k * exp(lpd_ik)). This is concave in w and solved here by the standard mixture-weight EM update (no external optimizer), which respects the simplex by construction. By default this compatibility surface returns only the weights. Set return_result=True to receive the weights together with convergence, iteration, objective, and termination metadata.

Parameters:

• pointwise_lpd (ndarray)

• iters (int)

• tol (float)

• return_result (bool)

Return type:

ndarray | dict

split_rhat(draws)

Maximum rank-normalized and folded split-R-hat for one parameter.

draws must be a finite (n_chains, n_draws) matrix with at least two chains and four draws per chain. Splitting catches within-chain non-stationarity; folding about the pooled median also catches scale mismatch. Values within ~0.01 of 1.0 indicate convergence; > 1.01 is a warning. Constant chains have no estimable scale and return NaN rather than a healthy-looking value.

Parameters:

draws (ndarray)

Return type:

float

bulk_ess(draws)

Bulk ESS of split rank-normalized draws; constant chains return NaN.

Parameters:

draws (ndarray)

Return type:

float

tail_ess(draws)

Tail ESS from split 5% and 95% quantile indicators; unavailable cases return NaN.

Parameters:

draws (ndarray)

Return type:

float

convergence_diagnostics(draws)

Return modern convergence metrics plus an explicit availability receipt.

Every unavailability is a RECEIPT, not an exception: a metric that cannot be computed for these draws comes back NaN, is listed under unavailable, and gets a plain-language entry in unavailable_because. In particular a single chain – the default configuration of every sampling route – receipts split_rhat (which compares independent chains) while still returning the finite single-chain bulk_ess/tail_ess estimates, instead of raising the 2-chain ValueError for one cause and receipting the other (t5 wave-3). Malformed input (wrong shape, NaN/Inf draws, too few draws to split) still raises.

Parameters:

draws (ndarray)

Return type:

dict

class ConformalRegressor(result, y_cal, *, given, alpha=0.1)

Bases: object

Split-conformal prediction intervals around a fitted regression result.

Calibrates the absolute-residual nonconformity score on held-out (given, y_cal) and produces symmetric intervals predict(x) +/- qhat with marginal coverage at least 1 - alpha. result is any object with a predict(given) method returning the fitted mean (a RegressionResult, a location-scale result, or a GP regressor).

Parameters:

• result (Any)

• y_cal (Any)

• given (dict)

• alpha (float)

interval(given)

Return (lower, upper) arrays of the conformal interval at covariates given.

Parameters:

given (dict)

Return type:

tuple[ndarray, ndarray]

covers(y, *, given)

Boolean array: does the interval at given contain each observed y.

Parameters:

• y (Any)

• given (dict)

Return type:

ndarray

class ConformalClassifier(proba_cal, y_cal, *, alpha=0.1)

Bases: object

Split-conformal label sets from per-class probabilities.

proba_cal is an (n, K) matrix of calibration probabilities p(y | x) (any proper classifier — a mixle generative classifier’s class posterior, a softmax, …) and y_cal the integer labels. The nonconformity score is 1 - p(true) and the calibrated set keeps every label whose score is within the conformal quantile, so it covers the true label with probability at least 1 - alpha and grows from one label (confident) to several (hedging) as the model is unsure.

Parameters:

• proba_cal (Any)

• y_cal (Any)

• alpha (float)

predict_set(proba)

Boolean (n, K) label-inclusion matrix at probabilities proba.

Parameters:

proba (Any)

Return type:

ndarray

covers(proba, y)

Boolean array: is each true label y in the predicted set.

Parameters:

• proba (Any)

• y (Any)

Return type:

ndarray

set_sizes(proba)

Number of labels in the predicted set for each row of proba.

Parameters:

proba (Any)

Return type:

ndarray

class ConformalQuantileRegressor(lo, hi, y_cal, *, given, alpha=0.1)

Bases: object

Conformalized quantile regression (Romano, Patterson, Candes 2019).

Combines two fitted quantile regressions (a lower and an upper conditional quantile) with a split-conformal calibration so the band has exact marginal coverage and the adaptive, heteroscedastic width of quantile regression — wide where the data is noisy, narrow where it is tight, unlike the constant-width absolute-residual band of ConformalRegressor.

The nonconformity score is the signed distance outside the predicted band, E_i = max(qlo(x_i) - y_i, y_i - qhi(x_i)) (negative when y_i is comfortably inside), and the calibrated band is [qlo(x) - qhat, qhi(x) + qhat] with qhat the conformal quantile of the calibration scores. lo and hi are fitted quantile-regression results (from ...fit(..., quantile=tau)), typically at tau = alpha/2 and 1 - alpha/2.

Parameters:

• lo (Any)

• hi (Any)

• y_cal (Any)

• given (dict)

• alpha (float)

interval(given)

Return (lower, upper) arrays of the calibrated adaptive band at covariates given.

Parameters:

given (dict)

Return type:

tuple[ndarray, ndarray]

covers(y, *, given)

Boolean array: does the adaptive band at given contain each observed y.

Parameters:

• y (Any)

• given (dict)

Return type:

ndarray

class ConformalStructure(dist, calibration, *, alpha=0.1)

Bases: object

Split-conformal credible sets over combinatorial structures (rankings, matchings, spanning trees, permutations, …) from a fitted mixle distribution’s exact log-density.

The nonconformity of a structure s is -log p(s): the lower its model probability, the more surprising it is. Calibrating on held-out true structures yields a log-probability threshold, and the conformal set is {s : log p(s) >= threshold} — it contains the true structure with probability at least 1 - alpha whenever the calibration and test structures are exchangeable (for example iid draws), with no assumption on the model being correct.

dist is any structure distribution exposing log_density (PlackettLuceDistribution, MallowsDistribution, MatchingDistribution, SpanningTreeDistribution, …); calibration is a sequence of observed structures. Membership is always available; listing or counting the set additionally needs the distribution’s exact enumerator().

Parameters:

• dist (Any)

• calibration (Any)

• alpha (float)

property log_prob_threshold: float

Structures with log p(s) at or above this value are in the conformal set.

contains(structure)

Is structure in the conformal set (its log-probability above the threshold).

Parameters:

structure (Any)

Return type:

bool

covers(structures)

Boolean array: membership of each structure (use on held-out truths to check coverage).

Parameters:

structures (Any)

Return type:

ndarray

members()

List the structures in the conformal set, highest-probability first.

Requires the distribution’s exact enumerator() (raises EnumerationError otherwise). The enumerator yields structures in descending log-probability, so the scan stops at the threshold.

Return type:

list

size()

Number of structures in the conformal set (needs the exact enumerator()).

Return type:

int

class ConformalLinkPredictor(edge_prob, cal_edges, *, alpha=0.1, allow_self_loops=False)

Bases: object

Split-conformal candidate-neighbor sets for a random-graph model from its edge-probability matrix P (P[i, j] = p(edge i--j), e.g. X @ X.T from a fitted RDPG, or an Erdos-Renyi / stochastic-block-model edge probability).

The nonconformity of a present edge (i, j) is 1 - P[i, j]. Calibrating on held-out true edges gives a threshold; the predicted neighbor set of a node keeps every candidate j with 1 - P[i, j] <= tau, so it contains a true neighbor with probability at least 1 - alpha over exchangeable held-out edges (a random split of the observed edges).

Parameters:

• edge_prob (Any)

• cal_edges (Any)

• alpha (float)

• allow_self_loops (bool)

neighbor_set(i, candidates=None)

Candidate nodes j in node i’s conformal neighbor set.

Parameters:

• i (int)

• candidates (Any)

Return type:

ndarray

covers(edges)

Boolean array: is each held-out true edge’s endpoint in the predicted neighbor set.

Parameters:

edges (Any)

Return type:

ndarray

set_sizes(nodes=None)

Neighbor-set size per node (defaults to all nodes).

Parameters:

nodes (Any)

Return type:

ndarray

class ConformalKnowledgeGraph(kg, calibration, *, slot='tail', alpha=0.1)

Bases: object

Split-conformal completion sets for a knowledge-graph model (any-slot UQ).

Calibrating on held-out true triples turns the model’s completion posterior into a set of candidate fillers that contains the true one with probability at least 1 - alpha over exchangeable held-out triples. slot selects which slot is predicted – 'tail' for (h, r, ?), 'head' for (?, r, t), 'relation' for (h, ?, t) – using the model’s tail_log_posterior / head_log_posterior / relation_log_posterior. The nonconformity of a triple is 1 - p(true filler), so a confident model gives small completion sets and a recommended completion carries a coverage guarantee.

Parameters:

• kg (Any)

• calibration (Any)

• slot (str)

• alpha (float)

completion_set(h=None, r=None, t=None)

Candidate fillers in the conformal set for the missing slot of a query.

Parameters:

• h (int | None)

• r (int | None)

• t (int | None)

Return type:

ndarray

covers(triples)

Boolean array: is each held-out true triple’s filler in the completion set.

Parameters:

triples (Any)

Return type:

ndarray

set_sizes(triples)

Completion-set size for each query (the slot of each triple is treated as missing).

Parameters:

triples (Any)

Return type:

ndarray

class ScalingLawDiagnostics(usable, status, acceptance_rate, posterior_draws, burn, n_chains, posterior_finite, caveats=())

Bases: object

Quality receipt required before a fitted law may drive an allocation.

Parameters:

• usable (bool)

• status (str)

• acceptance_rate (float)

• posterior_draws (int)

• burn (int)

• n_chains (int)

• posterior_finite (bool)

• caveats (tuple[str, ...])

class ScalingLawFit(fitted, n0, d0, diagnostics=None)

Bases: object

A fitted loss = E + A/N**alpha + B/D**beta scaling law with a genuine posterior.

fitted is the mixle.ppl RandomVariable returned by .fit(..., how="mcmc") – the “power-law leaf” – carrying MCMC draws over (E, A, alpha, B, beta, sigma) (the “uncertainty receipts”). n0/d0 are the normalization constants N/D were divided by before fitting (numerical conditioning only; predictions are in the original units).

Parameters:

• fitted (Any)

• n0 (float)

• d0 (float)

• diagnostics (ScalingLawDiagnostics | None)

samples(name)

Posterior draws for parameter name (one of E, A, alpha, B, beta, sigma).

Parameters:

name (str)

Return type:

ndarray

hdi(name, prob=0.9)

Highest-density credible interval for parameter name at coverage prob.

Parameters:

• name (str)

• prob (float)

Return type:

tuple[float, float]

require_usable()

Refuse scientific decisions from a fit that lacks a passing quality receipt.

Return type:

None

predict_mean(n, d)

Posterior expectation of mean loss, computed from correlated joint draws.

Parameters:

• n (float)

• d (float)

Return type:

float

predict_samples(n, d)

Posterior-predictive DRAWS of the mean loss at (n, d) – integrates over parameter uncertainty (no observation noise added), for building a predictive credible interval.

Parameters:

• n (float)

• d (float)

Return type:

ndarray

fit_scaling_law(records, *, draws=4000, burn=4000, scale=0.02, seed=None, rng=None)

Fit loss = E + A/N**alpha + B/D**beta to records (a list of (N, D, loss)).

Uses mixle’s OWN PPL fitting machinery, not scipy’s curve_fit: the nonlinear power-law mean is expressed as a potential (custom log-likelihood term) over free Normal-prior parameters, and the whole thing is fit as an ordinary mixle.ppl model with how="mcmc" – exactly the pattern examples/flagship_physics_inverse.py (the D-track’s physics-inverse-problem flagship) uses to turn an arbitrary nonlinear forward model into PPL evidence. Returns a ScalingLawFit carrying full MCMC posterior draws (real uncertainty, not just a point estimate).

scale is passed straight through to how="mcmc"’s adaptive random-walk proposal as its initial per-coordinate step size. The carrier observation used here is a single vacuous point (the potential IS the evidence, exactly as in the physics-inverse flagship), so mixle’s default proposal-scale heuristic (~ data_std / sqrt(n_data)) sees n_data=1 and starts far too wide relative to this tightly-peaked 6-parameter likelihood – left at its default it needs many thousands of extra burn-in draws for the adaptive proposal to shrink into a workable acceptance-rate regime. Setting a small explicit scale up front (tuned to this problem’s posterior widths, ~1e-2 on the log-parameter scale) restores good mixing (~15-25% acceptance) at the draws/burn defaults below.

Parameters:

• records (Sequence[tuple[float, float, float]])

• draws (int)

• burn (int)

• scale (float | None)

• seed (int | None)

• rng (RandomState | None)

Return type:

ScalingLawFit

generate_synthetic_chinchilla_data(n_points=60, *, seed=0, noise_sd=0.015, n_range=(1.0e7, 1.0e11), d_range=(1.0e8, 1.0e12), e=CHINCHILLA_E, a=CHINCHILLA_A, b=CHINCHILLA_B, alpha=CHINCHILLA_ALPHA, beta=CHINCHILLA_BETA)

SYNTHETIC (N, D, loss) triples generated from the published Chinchilla functional form.

N/D are drawn log-uniformly over n_range/d_range (spanning several orders of magnitude, the way a real training-run sweep would), the mean loss is the exact published power law E + A/N**alpha + B/D**beta, and i.i.d. Gaussian noise of scale noise_sd (in loss units) is added – a realistic per-run measurement/optimization-noise floor. See the module docstring for why this is synthetic-from-known-exponents rather than a real published per-run table (no network access in this environment).

Parameters:

• n_points (int)

• seed (int)

• noise_sd (float)

• n_range (tuple[float, float])

• d_range (tuple[float, float])

• e (float)

• a (float)

• b (float)

• alpha (float)

• beta (float)

Return type:

SyntheticScalingRecords

allocate_compute(fit, compute_budget, *, n_bounds=(1.0e7, 1.0e12), n_init=8, n_iter=20, seed=0, flops_per_token_param=FLOPS_PER_TOKEN_PARAM)

Find the (N, D) split minimizing fit’s predicted loss under C ~= 6*N*D.

Reuses mixle.doe.bayesopt.minimize() (GP-surrogate expected-improvement Bayesian optimization) – this codebase’s existing DOE machinery – rather than a bespoke optimizer. The compute constraint C = 6*N*D is an exact algebraic EQUALITY, not a black-box inequality, so instead of routing through mixle.doe.constrained (built for black-box inequality constraints, which this is not), it is eliminated by substitution: for any candidate N, D is set to exactly satisfy the constraint, collapsing the 2-D allocation problem to a 1-D search over log10(N) that bayesopt.minimize drives directly.

Parameters:

• fit (ScalingLawFit)

• compute_budget (float)

• n_bounds (tuple[float, float])

• n_init (int)

• n_iter (int)

• seed (int)

• flops_per_token_param (float)

Return type:

tuple[float, float]

class ScalingLawAllocationController(*, n_bounds=(1.0e7, 1.0e12), design=None, seed=None)

Bases: object

D5-pattern controller for the compute-allocation decision – shares D5’s controller brain.

D5 (mixle/inference/conditional_jit_controller.py, PR #163) defines a generic LearnedController[StateT, ActionT] base (select_action(state) -> action / update(state, action, gain, cost) -> None) plus a concrete DesignModelController that wraps mixle.task.edge.DesignModel – a GP-surrogate design space model, warm-startable across DIFFERENT tasks via a fingerprint vector – to propose a continuous 1-D knob (there, budget_fraction) from logged (state, action, gain, cost) rows, falling back to a fixed default before at least two rows are logged. D5’s own docstring explicitly anticipates this reuse: “a future F5 … item could subclass LearnedController directly for its own state/action types … reusing the bandit/ DesignModel wiring pattern without needing block-EM’s ControllerState/ControllerAction dataclasses at all.”

This class does exactly that: it is NOT a subclass of D5’s concrete DesignModelController (that class’s state/action types are block-EM-specific), but it reuses the SAME mixle.task.edge.DesignModel wiring D5’s DesignModelController uses, against F5’s own ScalingLawState (fitted-law-parameters + compute-budget fingerprint) and AllocationAction (log10(N)) types, with the identical cold-start-fallback and fingerprint-conditioned-proposal shape. It does not subclass D5’s mixle.inference.conditional_jit_controller.LearnedController ABC directly (importing D5’s module here would pull an unrelated inference-internals dependency into mixle.doe for a class whose only job is to satisfy the abstract two-method surface); the class shape below is deliberately identical to it so the substitutability the roadmap asks for (“shares the controller brain with D5”) is structural, not merely nominal.

Use allocate_compute_learned() to create a pending proposal, then call record_outcome() with a measured gain, measured cost, and provenance. Only accepted measured outcomes enter the warm-start design history.

Parameters:

• n_bounds (tuple[float, float])

• design (Any)

• seed (int | None)

record_outcome(proposal, *, realized_gain, realized_cost, provenance)

Attach one measured outcome to a pending proposal exactly once.

Parameters:

• proposal (AllocationProposal)

• realized_gain (float)

• realized_cost (float)

• provenance (str)

Return type:

AllocationObservationReceipt

allocate_compute_learned(fit, compute_budget, *, controller=None, flops_per_token_param=FLOPS_PER_TOKEN_PARAM)

Propose (N, D) and return a pending receipt without fabricating an outcome.

This is the OPTIONAL learned path (see the roadmap item’s “optionally wire in D5’s LearnedController pattern”); allocate_compute() (plain GP-BO via mixle.doe.bayesopt) is the primary, required allocator and is what the acceptance test compares against the fixed heuristic. With a fresh (cold) controller and no logged history this falls back to the bounds midpoint, so it is not expected to beat the heuristic on a single cold call – its payoff is warm-starting across many allocation decisions, the same story as D5’s DesignModelController.

Parameters:

• fit (ScalingLawFit)

• compute_budget (float)

• controller (ScalingLawAllocationController | None)

• flops_per_token_param (float)

Return type:

tuple[float, float, ScalingLawAllocationController, AllocationProposal]

Submodules

• mixle.ppl.autograd module
• mixle.ppl.conformal module
• mixle.ppl.core module
• mixle.ppl.density module
• mixle.ppl.diagnostics module
• mixle.ppl.distributions module
• mixle.ppl.field module
• mixle.ppl.guide module
• mixle.ppl.inference module
• mixle.ppl.neural module
• mixle.ppl.predictive module
• mixle.ppl.priors module
• mixle.ppl.provenance module
• mixle.ppl.regression module
• mixle.ppl.rough_paths module
• mixle.ppl.scaling_laws module
• mixle.ppl.statespace module
• mixle.ppl.summarize module
• mixle.ppl.survival module
• mixle.ppl.vmp module