mixle.ppl.core module

Core of the mixle.ppl probabilistic-programming surface.

One immutable wrapper type, RandomVariable, sits over mixle’s existing distribution / estimator / sampler objects. It adds no inference engine: every call lowers (one routing site, lower()) to machinery that already exists and then dispatches.

The core surface covers immutable random variables, free parameter holes, fitting, and lowering into existing Mixle estimators and distributions. Unsupported combinations should fail explicitly rather than returning a partially lowered model.

class RandomVariable(kind, *, family=None, args=(), name=None, keys=None, dist=None, result=None, scope='shared', reparam=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 (PosteriorResult | 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.

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 sample expression contains one or more free parameters.

property name: str | None

Return the optional user-visible variable name.

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)

Draw samples from the represented distribution or derived expression.

Parameters:

• n (int | None)

• seed (int | None)

log_prob(x)

Evaluate the log probability or log density at x.

log_likelihood(data)

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

Return type:

float

aic(data, k=None)

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)

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

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

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} (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)

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:

• samples (int)

• seed (int)

var(samples=20000, seed=0)

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

Parameters:

• samples (int)

• seed (int)

prob(samples=40000, seed=999)

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

Parameters:

• samples (int)

• seed (int)

prob_of_event()

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

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

Parameters:

n (int)

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. Raises if that record is unavailable (e.g. the model was reloaded from a saved artifact, which does not round-trip it) – 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)

register_family(name, dist_cls, est_cls, to_dist, arity, seed_at=None, positive=None, init_fit=None, read=None, support=None)

Register a flat PPL family and its distribution/estimator lowering rules.

Return type:

Family

class Family(name, dist_cls, est_cls, to_dist, arity, seed_at=None, positive=None, init_fit=None, read=None, support=None)

Bases: object

Lowering recipe for one distribution family.

Keeps the alias namespace and the engine objects in one place so the wrapper never hard-codes a distribution. to_dist maps user-facing (conventional) arguments to the underlying *Distribution kwargs; make_estimator builds the paired *Estimator for the all-free case.

make_dist(args, name)

Construct the concrete distribution for conventional PPL arguments.

Parameters:

• args (tuple[Any, ...])

• name (str | None)

make_estimator(name, keys)

Construct the estimator associated with this family.

Parameters:

• name (str | None)

• keys (str | None)

class Constraint(leaves, pred, desc, residual=None, soft=False)

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.

contains(x)

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

Event

alias of Constraint

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

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

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