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, group_by=None)

Bases: object

The single user-facing PPL type (immutable).

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

Parameters:

result (Any | None)

property certificate

The estimation certificate, when a fit attached one.

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

each(by=None)

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

Two data layouts are supported:

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

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

Parameters:

by (str | None)

Return type:

RandomVariable

noncentered()

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

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

Return type:

RandomVariable

property scope: str

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

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

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

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

Parameters:

• name (str | None)

• result (Any | None)

Return type:

RandomVariable

independent()

Return an explicitly independent copy of one atomic random variable.

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

Return type:

RandomVariable

exp()

Return the deterministic exponential transform of this random variable.

Return type:

RandomVariable

log()

Return the deterministic logarithm transform of this random variable.

Return type:

RandomVariable

eq(other)

Build an equality constraint against a value or another expression.

ne(other)

Build an inequality constraint against a value or another expression.

given(constraint)

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

Return type:

RandomVariable

property is_bound: bool

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

property has_free: bool

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

property name: str | None

Return the optional user-visible variable name.

property parameter_dimension: int | None

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

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

property group_index: dict | None

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

property group_labels: tuple | None

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

property columns: list

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

Type:

For a constrain(...) joint RV

property dist

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

property components

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

property params

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

property result: Any | None

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

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

Draw samples from the represented distribution or derived expression.

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

Parameters:

• n (int | None)

• seed (int | None)

• size (int | None)

• max_attempts (int)

log_prob(x)

Evaluate the log probability or log density at x.

log_density(x)

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

log_likelihood(data)

Total log-likelihood of data under the fitted model.

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

Return type:

float

aic(data, k=None)

Akaike information criterion (lower is better).

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

Parameters:

k (int | None)

Return type:

float

bic(data, k=None)

Bayesian information criterion (lower is better).

Parameters:

k (int | None)

Return type:

float

pointwise_log_likelihood(data)

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

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

Return type:

ndarray

plugin_log_likelihood(data)

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

Return type:

ndarray

waic(data)

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

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

Return type:

dict

loo(data)

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

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

Return type:

dict

summary()

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

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

mean(samples=20000, seed=0)

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

Parameters:

• samples (int)

• seed (int)

var(samples=20000, seed=0)

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

Parameters:

• samples (int)

• seed (int)

prob(samples=40000, seed=999)

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

Parameters:

• samples (int)

• seed (int)

prob_of_event(samples=40000, seed=999)

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

Parameters:

• samples (int)

• seed (int)

predict(n=1, rng=None)

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

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

Parameters:

n (int | Mapping[str, Any])

posterior(x)

Posterior over a latent or a parameter.

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

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

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

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

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

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

Return type:

dict

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

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

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

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

Parameters:

• data (Sequence[Any])

• how (str)

• max_its (int)

• delta (float)

• backend (str)

• num_workers (int | None)

• engine (Any)

• precision (Any)

• print_iter (int)

• missing (str)

Return type:

RandomVariable

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

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

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

Parameters:

• rv (RandomVariable)

• target (str)

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

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

holds_fixed names the conventional parameters this family’s EM estimator leaves alone, so a free written in one of those slots is disclosed rather than silently ignored.

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, validator=None, holds_fixed=())

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.

validate_args(args)

Validate fixed conventional parameters without rejecting symbolic parameter expressions.

Parameters:

args (tuple[Any, ...])

Return type:

None

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, reduction='all', margin=None)

Bases: object

A boolean relation over one or more random variables.

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

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

property rv

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

eval(env)

Evaluate the constraint predicate against an environment of RV values.

eval_rows(env, rows=None)

Evaluate and reduce to exactly one boolean per sampled row.

Parameters:

rows (int | None)

Return type:

ndarray

with_reduction(reduction)

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

Parameters:

reduction (str)

Return type:

Constraint

contains(x)

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

Event

alias of Constraint

class ProbabilityEstimate(hits, trials, seed)

Bases: float

A Monte Carlo probability plus the experiment needed to interpret it.

The numeric value remains backward-compatible with ordinary probability arithmetic. hits, trials, seed, and a Wilson 95% interval distinguish zero observed hits from a mathematical proof that the event has probability zero.

Parameters:

• hits (int)

• trials (int)

• seed (int)

as_dict()

Return a serialization-friendly experiment receipt.

Return type:

dict

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