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 86 mixle.stats distribution classes are untouched; this is a thin, optional dialect. See notes/ppl-syntax-spec.md.

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); fields/regression/conformal and the PDE-inverse stack in their own modules. This package init only re-exports.

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

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 (invariant I3); validity depends on state. Construct via the family functions in mixle.ppl or fit.

Parameters:

result (PosteriorResult | None)

property certificate

The estimation certificate, when a fit attached one (penalized fits downgrade honestly; E2).

each(by=None)[source]

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

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
exp()[source]
Return type:

RandomVariable

log()[source]
Return type:

RandomVariable

eq(other)[source]
ne(other)[source]
given(constraint)[source]

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
property has_free: bool
property name: str | None
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: PosteriorResult | None

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

sample(n=None, seed=None)[source]
Parameters:

  • n (int | None)

  • seed (int | None)

log_prob(x)[source]
log_likelihood(data)[source]

Total log-likelihood of data under the fitted model (sum of log_prob).

Return type:

float

aic(data, k=None)[source]

Akaike information criterion (lower is better). k defaults to a heuristic parameter count from .params.

Parameters:

k (int | None)

Return type:

float

bic(data, k=None)[source]

Bayesian information criterion (lower is better).

Parameters:

k (int | None)

Return type:

float

pointwise_log_likelihood(data)[source]

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; for a point-estimate fit it is a single row.

Return type:

ndarray

waic(data)[source]

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 – and falls back to a point estimate for non-Bayesian fits.

Return type:

dict

loo(data)[source]

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

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} (the 95% credible interval) plus _acceptance_rate and, for multi-chain runs, _rhat / _ess / _n_chains. For map/em it returns .params.

mean(samples=20000, seed=0)[source]

Expected value of the random variable (Monte-Carlo; works for any RV — concrete, transformed, convolved, or conditioned). For a joint constrain(...) RV this is the per-variable mean vector.

Parameters:
var(samples=20000, seed=0)[source]

Variance of the random variable (Monte-Carlo); per-variable for a joint RV.

Parameters:
prob(samples=40000, seed=999)[source]

Probability that the relation holds (Monte-Carlo), for a constrain(...) RV.

Parameters:
prob_of_event()[source]

P(event) under the base distribution (Monte-Carlo), for a conditioned RV.

predict(n=1, rng=None)[source]

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

Parameters:

n (int)

posterior(x)[source]

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

Report which inference route .fit(how=...) will take, and why – without fitting.

Returns {'route', 'reason', 'caveats'}. This is the teaching surface for mixle’s automatic cross-family inference selection: rv.explain_fit() answers “how will this be fit, and what are the honest limits of that choice?”. The route mirrors fit() exactly (it shares _resolve_auto() for the flat tree and re-checks the same structural short-circuits).

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

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:
Return type:

RandomVariable

lower(rv, *, target='dist')[source]

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 RV (I7).

Parameters:

  • rv (RandomVariable)

  • target (str)

class Guide(**latents)[source]

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

tuple[str, …]

class StructuredVIPosterior(gres, guide)[source]

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

dict

mean(name)[source]
samples(name, n=4000, rng=None)[source]
Parameters:

n (int)

summary()[source]
Return type:

dict

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

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

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:
class AdmixturePosterior(lam, gamma, topic_handles, ll_trace)[source]

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()[source]
doc_topics(d=None)[source]
posterior(topic)[source]
summary()[source]
Return type:

dict

Normal(mean, sd, *, name=None, keys=None)[source]

Normal with mean and standard deviation (lowers to GaussianDistribution(mu, sd**2)).

Parameters:
Return type:

RandomVariable

Poisson(rate, *, name=None, keys=None)[source]
Parameters:

  • rate (Any)

  • name (str | None)

  • keys (str | None)

Return type:

RandomVariable

Gamma(shape, rate, *, name=None, keys=None)[source]

Gamma with shape and rate (lowers to GammaDistribution(k=shape, theta=1/rate)).

Parameters:

  • shape (Any)

  • rate (Any)

  • name (str | None)

  • keys (str | None)

Return type:

RandomVariable

Exponential(rate, *, name=None, keys=None)[source]

Exponential with rate (mean 1/rate; lowers to ExponentialDistribution(beta=1/rate)).

Parameters:

  • rate (Any)

  • name (str | None)

  • keys (str | None)

Return type:

RandomVariable

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

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

  • p (Any)

  • name (str | None)

  • keys (str | None)

Return type:

RandomVariable

Geometric(p, *, name=None, keys=None)[source]
Parameters:

  • p (Any)

  • name (str | None)

  • keys (str | None)

Return type:

RandomVariable

Binomial(n, p, *, name=None, keys=None)[source]

Binomial with n trials and success probability p (n is fixed/known).

Parameters:
Return type:

RandomVariable

Weibull(shape, scale, *, name=None, keys=None)[source]

Weibull with shape (k) and scale (lambda).

Parameters:

  • shape (Any)

  • scale (Any)

  • name (str | None)

  • keys (str | None)

Return type:

RandomVariable

Laplace(loc, scale, *, name=None, keys=None)[source]

Laplace (double-exponential) with location and scale (b).

Parameters:
Return type:

RandomVariable

Logistic(loc, scale, *, name=None, keys=None)[source]

Logistic with location and scale.

Parameters:
Return type:

RandomVariable

Uniform(low, high, *, name=None, keys=None)[source]

Continuous uniform on [low, high].

Parameters:
Return type:

RandomVariable

Rayleigh(sigma, *, name=None, keys=None)[source]

Rayleigh with scale sigma.

Parameters:

  • sigma (Any)

  • name (str | None)

  • keys (str | None)

Return type:

RandomVariable

Pareto(scale, shape, *, name=None, keys=None)[source]

Pareto with minimum value xm (scale) and tail index alpha (shape).

Parameters:

  • scale (Any)

  • shape (Any)

  • name (str | None)

  • keys (str | None)

Return type:

RandomVariable

Beta(a, b, *, name=None, keys=None)[source]
Parameters:
Return type:

RandomVariable

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

Student-t with degrees of freedom, location, scale (heavy-tailed Normal).

Parameters:
Return type:

RandomVariable

LogNormal(mu, sigma, *, name=None, keys=None)[source]

Log-normal: log(X) ~ Normal(mu, sigma).

Parameters:
Return type:

RandomVariable

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

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). The MLE is iterative with no closed form, so EMG(free, free, free).fit(data) uses a consistent method-of-moments estimate.

Parameters:
Return type:

RandomVariable

NegativeBinomial(r, p, *, name=None, keys=None)[source]

Negative binomial with r failures and success probability p.

Parameters:
Return type:

RandomVariable

HalfNormal(sigma, *, name=None, keys=None)[source]

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

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

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

Parameters:
Return type:

RandomVariable

Gumbel(loc, scale, *, name=None, keys=None)[source]

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

Parameters:
Return type:

RandomVariable

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

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

Parameters:
Return type:

RandomVariable

Skellam(mu1, mu2, *, name=None, keys=None)[source]

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

Parameters:
Return type:

RandomVariable

LogSeries(p, *, name=None, keys=None)[source]

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

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

Parameters:
Return type:

RandomVariable

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

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

Parameters:
Return type:

RandomVariable

Tweedie(mu, phi, *, name=None, keys=None)[source]

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:
Return type:

RandomVariable

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

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:
Return type:

RandomVariable

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

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

Parameters:
Return type:

RandomVariable

Nakagami(m, omega, *, name=None, keys=None)[source]

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

Parameters:
Return type:

RandomVariable

Rician(nu, sigma, *, name=None, keys=None)[source]

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

Parameters:
Return type:

RandomVariable

Dirichlet(alpha, *, dim=None, name=None, keys=None)[source]

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

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)

field
channels
out
build(in_shape)[source]
Parameters:

in_shape (Any)

Return type:

Any

Embedding

alias of CategoricalEmbedding

Mix(components, weights=None, *, name=None)[source]

Finite mixture over component RandomVariables (or concrete distributions).

Mix([Normal(free, free), Normal(free, free)]).fit(data) fits a 2-component Gaussian mixture; .posterior(data) returns the responsibilities.

Parameters:

name (str | None)

Return type:

RandomVariable

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

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)

field
hidden
out
build(in_shape)[source]
Parameters:

in_shape (Any)

Return type:

Any

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

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)

field
out
d_model
n_layer
n_head
embedding
build(in_shape)[source]
Parameters:

in_shape (Any)

Return type:

Any

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

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

Parameters:
Return type:

RandomVariable

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

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

Parameters:
Return type:

RandomVariable

VAE(dim, *, latent=2, hidden=32, m_steps=120, lr=5e-3)[source]

A latent-variable p(x) over R^dim via a VAE. log_density is the ELBO (a lower bound); fit .fit(x).

Parameters:
Return type:

RandomVariable

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

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

Parameters:
Return type:

RandomVariable

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

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

Parameters:
Return type:

RandomVariable

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

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

Parameters:
Return type:

RandomVariable

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

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

Parameters:
Return type:

RandomVariable

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

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

Parameters:
Return type:

RandomVariable

SemiMix(components, weights=None, *, name=None)[source]

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

IID sequence of element. Fit on a list of sequences (each a list/array).

Parameters:

name (str | None)

Return type:

RandomVariable

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

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

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:
Return type:

RandomVariable

MVN(dim, *, mean=None, cov=None, name=None)[source]

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:
Return type:

RandomVariable

DiagGaussian(dim, *, mean=None, var=None, name=None)[source]

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:
Return type:

RandomVariable

LocalLevel(*, name=None)[source]

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

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

Parameters:

name (str | None)

Return type:

RandomVariable

Graph()[source]

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

class Field(name)[source]

Bases: object

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

Parameters:

name (str)

name
class Group(name, slopes=())[source]

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)

name
slopes
compare(models, data, *, by='aic')[source]

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

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:
Return type:

dict[str, Any]

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

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:
Return type:

dict[str, Any]

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

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:
Return type:

dict[str, Any]

fit_censored(model, time, *, event=None, lower=None, upper=None, seed=0)[source]

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:
Return type:

RandomVariable

fit_with_provenance(rv, data, *, seed=None, **fit_kw)[source]

Fit a PPL RandomVariable on data and return (fitted_rv, header) with full provenance.

fit_kw is passed through to rv.fit (how=, max_its=, delta=, backend=, …). The header is built from the fitted model’s lowered distribution so it carries schema + final log-likelihood where available; method records the requested how. The header is returned regardless; it is also attached as fitted.header when the RandomVariable allows it (RVs with __slots__ do not).

Parameters:
censored_loglik(dist, time, *, event=None, lower=None, upper=None)[source]

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

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:
Return type:

dict[str, ndarray]

hdi(samples, prob=0.94)[source]

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:
Return type:

tuple[float, float]

posterior_summary(fitted, *, hdi_prob=0.94)[source]

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

Returns {param_name: {'mean', 'sd', 'hdi_low', 'hdi_high', 'ess', 'r_hat'}}. 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. A point fit (em/map) yields just mean/sd.

Parameters:

  • fitted (RandomVariable)

  • hdi_prob (float)

Return type:

dict[str, dict[str, Any]]

constrain(*constraints)[source]

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

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

leaves
pred
desc
residual
soft
property rv

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

eval(env)[source]
contains(x)[source]

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

Event

alias of Constraint

eq(lhs, rhs)[source]

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

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

Return type:

Constraint

potential(fn, *vars, name=None)[source]

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

vars are random-variable parameters of the model (named priors, or param(...) vector/matrix handles) – exactly the references a constrain()/eq() constraint may use. 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); like constraints, every referenced variable must be a parameter of the fitted model.

Return type:

_Potential

increasing(v, *, strict=False)[source]

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

Parameters:

strict (bool)

Return type:

Constraint

decreasing(v, *, strict=False)[source]

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

Parameters:

strict (bool)

Return type:

Constraint

monotone(v)[source]

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

Return type:

Constraint

convex(v)[source]

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

Return type:

Constraint

concave(v)[source]

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

Return type:

Constraint

lipschitz(v, bound)[source]

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

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:
Return type:

Constraint

class GaussianField(index, kernel, name='field')[source]

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 = 'field'
class FieldSystem(fields, coregion=None)[source]

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 = None
class FieldKernel[source]

Bases: object

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

precision(index)[source]
Parameters:

index (ndarray)

Return type:

ndarray

covariance(index)[source]

The prior covariance matrix, when available cheaply (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)[source]

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 = 1.0
amplitude: float = 1.0
jitter: float = 1e-06
radius: float = 1.0
covariance(index)[source]

The prior covariance matrix, when available cheaply (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

precision(index)[source]
Parameters:

index (ndarray)

Return type:

ndarray

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

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 = 1.0
amplitude: float = 1.0
nu: float = 1.5
jitter: float = 1e-06
radius: float = 1.0
covariance(index)[source]

The prior covariance matrix, when available cheaply (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

precision(index)[source]
Parameters:

index (ndarray)

Return type:

ndarray

great_circle_distance(latlon_a, latlon_b=None, *, radius=1.0)[source]

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:
class RandomWalk(scale=1.0, order=1, ridge=None)[source]

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 = 1.0
order: int = 1
ridge: float | None = None
precision(index)[source]
Parameters:

index (ndarray)

Return type:

ndarray

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

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 = 1.0
amplitude: float = 1.0
jitter: float = 1e-06
covariance(index)[source]

The prior covariance matrix, when available cheaply (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

precision(index)[source]
Parameters:

index (ndarray)

Return type:

ndarray

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

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] = (1.0, 1.0)
angle: float = 0.0
amplitude: float = 1.0
jitter: float = 1e-06
metric: ndarray | None = None
covariance(index)[source]

The prior covariance matrix, when available cheaply (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

precision(index)[source]
Parameters:

index (ndarray)

Return type:

ndarray

class Proxy[source]

Bases: object

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

prefix: str = 'proxy'
field: str | None = None
on(field_name)[source]

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

Parameters:

field_name (str)

Return type:

Proxy

params()[source]
Return type:

list[_ParamSpec]

loglik(field_t, params, torch)[source]
Parameters:
Return type:

Any

residual(field_t, params, torch)[source]

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:
Return type:

Any

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

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 = None
slope: Any = 1.0
intercept: Any = 0.0
scale: Any = 1.0
prefix: str = 'gauss'
params()[source]
Return type:

list[_ParamSpec]

loglik(field_t, params, torch)[source]
residual(field_t, params, torch)[source]

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

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

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 = 2.0
prefix: str = 'niche'
params()[source]
Return type:

list[_ParamSpec]

loglik(field_t, params, torch)[source]
class PoissonProxy(counts, index=None, offset=0.0, prefix='cox')[source]

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 = None
offset: Any = 0.0
prefix: str = 'cox'
loglik(field_t, params, torch)[source]
class CustomProxy(loglik_fn, param_specs=(), prefix='custom')[source]

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 = 'custom'
params()[source]
Return type:

list[_ParamSpec]

loglik(field_t, params, torch)[source]
fit_field(field, proxies, *, how='laplace', max_iter=500, lr=0.4, init=None, vi_steps=400, vi_lr=0.05, vi_samples=4)[source]

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 is the posterior for the sparse adjoint solve (where how='laplace' cannot run) and is exact for a linear forward; how='vi' fits a mean-field Gaussian variational posterior (reparameterized ADVI on the unconstrained vector), the calibrated 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. Posteriors over any node are read off the returned FieldPosterior.

Parameters:
Return type:

FieldPosterior

class FieldPosterior(map_values, _cov, _layout, _field_name, _hessian, objective, _field_prior=<factory>, _proxy_info=<factory>, _supports=<factory>, _marg_var=<factory>)[source]

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 Fisher information. posterior(field, coupling=False) instead inverts only the field’s own Hessian block (the posterior conditional on the other nodes at their MAP), which is the per-proxy additive-information picture.

Parameters:
map_values: dict
objective: float
mean(node)[source]
Parameters:

node (str)

Return type:

ndarray

cov(node)[source]

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

Parameters:

node (str)

Return type:

ndarray

sd(node)[source]
Parameters:

node (str)

Return type:

ndarray

posterior(node, *, coupling=True)[source]

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

Parameters:
Return type:

tuple[ndarray, ndarray]

field_posterior(include=None)[source]

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

Because information is additive, the field posterior precision under any subset of proxies is the prior precision plus those proxies’ Fisher-information 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). 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)[source]

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 come back at 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()[source]
Return type:

dict

class GP(name, index, kernel)[source]

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

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

Return type:

tuple

Niche(presence, *, over, mu_scale=2.0)[source]

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

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

Return type:

tuple

joint(observations)[source]

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

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)

field: GaussianField | None
proxies: list
fit(*, how='laplace', **kw)[source]
Parameters:

how (str)

Return type:

FieldPosterior

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

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:
Return type:

tuple

Potts(over, levels, *, weight=1.0)[source]

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

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:
Return type:

FieldPosterior

conformal(result, y_cal, *, given, alpha=0.1)[source]

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:
Return type:

ConformalRegressor

loo_stack(logliks)[source]

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, and the stacked_elpd_loo of the weighted predictive (which is >= the best single-model elpd_loo, since a one-hot weight is feasible).

Parameters:

logliks (Sequence[ndarray])

Return type:

dict

loo_stacking_weights(pointwise_lpd, iters=2000, tol=1.0e-10)[source]

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.

Parameters:
Return type:

ndarray

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

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

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

Parameters:

given (dict)

Return type:

tuple[ndarray, ndarray]

covers(y, *, given)[source]

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

Parameters:
Return type:

ndarray

class ConformalClassifier(proba_cal, y_cal, *, alpha=0.1)[source]

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

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

Parameters:

proba (Any)

Return type:

ndarray

covers(proba, y)[source]

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

Parameters:
Return type:

ndarray

set_sizes(proba)[source]

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

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

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

Parameters:

given (dict)

Return type:

tuple[ndarray, ndarray]

covers(y, *, given)[source]

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

Parameters:
Return type:

ndarray

class ConformalStructure(dist, calibration, *, alpha=0.1)[source]

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

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

Parameters:

structure (Any)

Return type:

bool

covers(structures)[source]

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

Parameters:

structures (Any)

Return type:

ndarray

members()[source]

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

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

Return type:

int

class ConformalLinkPredictor(edge_prob, cal_edges, *, alpha=0.1)[source]

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)

neighbor_set(i, candidates=None)[source]

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

Parameters:
Return type:

ndarray

covers(edges)[source]

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

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

Parameters:

nodes (Any)

Return type:

ndarray

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

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

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

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

Parameters:

triples (Any)

Return type:

ndarray

set_sizes(triples)[source]

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

Parameters:

triples (Any)

Return type:

ndarray

Submodules