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. See notes/ppl-syntax-spec.md for the design charter and invariants.

Slice 1-2 of the build order: the wrapper + free holes + fit (EM core). Algebra, conditioning, latent-arg Compound, and MCMC routing land in later slices.

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

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

exp()

Return type:

RandomVariable

log()

Return type:

RandomVariable

eq(other)

ne(other)

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

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)

Parameters:

• n (int | None)

• seed (int | None)

log_prob(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=...) 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)

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

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)

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.

name

dist_cls

est_cls

to_dist

arity

seed_at

positive

support

init_fit

read

make_dist(args, name)

Parameters:

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

• name (str | None)

make_estimator(name, keys)

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

leaves

pred

desc

residual

soft

property rv

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

eval(env)

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