mixle.inference package

The inference concern — fit a model and quantify its parameters.

One home for turning data into a fitted/posterior model. Every entry point is the same idea — infer parameters from data — and differs only in what it requires of its input:

• closed-form conjugate Bayes (conjugate_posterior) — needs a ConjugateUpdatable family;

• MLE / EM / MAP (fit / optimize / run_em) — needs an Estimator for the model;

• sampling-based inference (nuts / advi) — needs a sampleable / differentiable target (a log-density callable, or a model that can be sampled). This is exactly why it belongs here and not in a separate package: it is inference under a capability precondition on the target, nothing more.

Everything physically lives in this package: the estimation / EM / fit / objectives / Fisher machinery (mixle.inference.{estimation,em,fit,objectives,fisher}), the MCMC samplers (mixle.inference.mcmc), the engine-agnostic NUTS/ADVI target facade (mixle.inference.target + .backends + .diagnostics). Conjugate Bayes is re-exported from its canonical home mixle.stats.bayes.

These imports are eager and cycle-free: the machinery’s only mixle.stats dependency is the compute layer (mixle.stats.compute.{pdist,sequence}), never the mixle.stats package surface — the vectorized seq_* drivers were moved out of mixle.stats.__init__ into compute.sequence for exactly this reason.

class PublicAPIEntry(name, source, status, lazy)

Bases: object

One name in the package export manifest.

Parameters:

• name (str)

• source (str)

• status (str)

• lazy (bool)

public_api_manifest()

Return the immutable, ordered package export manifest.

Return type:

tuple[PublicAPIEntry, …]

class Explanation(total, parts=<factory>, correction=0.0, responsibilities=None, component=None)

Bases: object

Exact additive attribution of a log-density (or a decision margin between two hypotheses).

correction is the explicitly named non-additive residual – 0.0 for purely additive structures (Composite, Bayesian network, and any margin comparison, since the normalizer cancels there); the logsumexp term for a Mixture’s absolute log p(x). sum(v for _, v in parts) + correction equals total exactly, for every supported structure – this is asserted to machine precision in tests.

Parameters:

• total (float)

• parts (list[tuple[str, float]])

• correction (float)

• responsibilities (ndarray | None)

• component (int | None)

most_anomalous(k=3)

Return the top contribution terms by anomaly score order.

Parameters:

k (int)

Return type:

list[tuple[str, float]]

ledger_sum()

sum(parts) + correction – should equal total to machine precision; see class docstring.

Return type:

float

is_exact(atol=1e-9)

Return whether ledger parts plus correction reconstruct the total.

Parameters:

atol (float)

Return type:

bool

summary()

Render a human-readable anomaly ledger summary.

Return type:

str

explain(model, x)

Exact per-part attribution of model.log_density(x) (see module docstring).

Parameters:

• model (Any)

• x (Any)

Return type:

Explanation

explain_margin(model, answer, runner_up)

Exact per-part attribution of the decision margin log p(answer) - log p(runner_up).

answer/runner_up are two full candidate records for Composite/Bayesian-network models, so a field that is identical between them contributes exactly 0 – the diagnostic for a single corrupted field. For a Mixture used as a generative classifier, use explain_margin_mixture() instead (the margin there needs the observed point plus two component indices, not two full records).

Parameters:

• model (Any)

• answer (Any)

• runner_up (Any)

Return type:

Explanation

explain_margin_mixture(model, x, answer, runner_up)

Exact per-part attribution of the decision margin between two Mixture components at one point x.

answer/runner_up are component indices (e.g. class labels) of a Mixture used as a generative classifier: the margin is (log_w[answer] + log p_answer(x)) - (log_w[runner_up] + log p_runner_up(x)). The mixture’s logsumexp normalizer – the same term explain() must name as a correction for the absolute log p(x) – cancels exactly in this subtraction, so the margin ledger needs no correction (computed and asserted at 0.0, not assumed away).

Both hypotheses must be exact, distinct, in-range integer component indices: a margin of a hypothesis against itself is 0 by construction and is a caller error, not evidence, so it is rejected rather than reported as a confident tie.

Raises:

• TypeError – if either hypothesis is not an exact non-Boolean integer.

• IndexError – if either hypothesis is outside range(len(model.components)).

• ValueError – if both hypotheses name the same component.

Parameters:

• model (Any)

• x (Any)

• answer (int)

• runner_up (int)

Return type:

Explanation

class FaultReport(dominant, evidence=<factory>, suggested_fix='', receipt=<factory>)

Bases: object

Structural diagnosis built from explain() ledgers over failing cases.

dominant names the structural element most responsible, when one rises above ordinary case-to-case variability. suggested_fix is one of _FIX_VOCAB or empty. evidence ranks elements by adverse contribution.

Parameters:

• dominant (str)

• evidence (list[tuple[str, float]])

• suggested_fix (str)

• receipt (dict[str, Any])

diagnose(model, cases, *, background=None, min_z=1.0, co_occurrence_threshold=0.5)

Aggregate explain() ledgers into a structural fault report.

Each case’s per-part contributions are compared to background (a reference sample of typical cases; defaults to cases themselves, though a real, separately-supplied background is needed to detect a fault that is systematic across every case, since self-baselining against the failing set cancels a shift common to all of it) via a robust z-score (median/MAD per part name), so “adverse” is relative to that part’s own normal variability, not a raw log-density magnitude.

A single part scoring far from baseline is not, by itself, evidence of a structural defect. The closed-vocabulary fault this function actively detects is a missing dependency: two parts that are adverse on the same cases more often than chance predicts, indicating a possible unmodeled edge. When no such co-anomalous pair is found, the report remains empty/low-severity rather than guessing a fix.

min_z is a robust z-score cutoff (finite, >= 0) and co_occurrence_threshold is a rate in [0, 1]; out-of-range values are rejected rather than silently turning both comparisons into tautologies that prescribe an add_edge backed by no co-anomalous case at all.

Parameters:

• model (Any)

• cases (Sequence[Any])

• background (Sequence[Any] | None)

• min_z (float)

• co_occurrence_threshold (float)

Return type:

FaultReport

class InterventionalNetwork(net, interventions, identification=None)

Bases: object

A graph-surgery scenario with an explicit causal-identification status.

Parameters:

• net (Any)

• interventions (dict[int, Any])

• identification (CausalIdentification | None)

sample(size=1, *, seed=None)

Ancestral sampling with the intervened fields clamped (their factors are never consulted).

Parameters:

• size (int)

• seed (int | None)

Return type:

list[tuple]

expectation(field, *, n=4000, seed=0)

Monte-Carlo expectation in this structural scenario.

Parameters:

• field (int)

• n (int)

• seed (int)

Return type:

float

distribution(field, *, n=4000, seed=0)

Structural-scenario marginal of a discrete field as {value: probability}.

Parameters:

• field (int)

• n (int)

• seed (int)

Return type:

dict[Any, float]

class CausalIdentification(graph_source, estimand, evidence, assumptions, exchangeability, positivity, consistency, no_interference, structural_counterfactuals=False)

Bases: object

Auditable assumptions under which a structural contrast is interpreted causally.

This record does not prove its assertions. It makes the assertions reviewable and prevents a learned observational arrow from silently becoming a causal result. graph_source must say how the direction was established, and evidence should identify the design, protocol, or domain artifact supporting the claim.

Parameters:

• graph_source (str)

• estimand (str)

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

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

• exchangeability (bool)

• positivity (bool)

• consistency (bool)

• no_interference (bool)

• structural_counterfactuals (bool)

property identified: bool

Whether all common treatment-effect identification assumptions were declared.

classmethod domain_asserted(evidence, *, estimand='average treatment effect', assumptions=('causal DAG supplied independently of fitted associations',), structural_counterfactuals=False)

Construct a receipt for a domain-asserted causal DAG.

This convenience constructor records that the caller accepts the standard identification assumptions. Applications needing qualified assumptions should instantiate the dataclass directly rather than use this shorthand.

Parameters:

• evidence (str)

• estimand (str)

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

• structural_counterfactuals (bool)

Return type:

CausalIdentification

to_dict()

Return a serializable copy of the receipt.

Return type:

dict[str, Any]

average_causal_effect(net, treatment, a, b, outcome, *, identification=None, n=4000, seed=0)

Return an identified average causal contrast under an explicit assumption receipt.

Parameters:

• net (Any)

• treatment (int)

• a (Any)

• b (Any)

• outcome (int)

• identification (CausalIdentification | None)

• n (int)

• seed (int)

Return type:

float

counterfactual(net, observed, interventions, *, identification=None)

What this observed record would have been under the intervention (abduction-action-prediction).

Per Pearl’s three steps, walked in topological order:

• abduction – a linear-Gaussian field’s exogenous noise is point-identified from the row: its residual eps = observed - coef @ parents_observed;

• action – intervened fields take their do values;

• prediction – the same residual replays through the counterfactual parents: cf = coef @ parents_cf + eps.

An identification receipt with structural_counterfactuals=True is required. Boundaries: (1) a field that is not linear-Gaussian keeps its observed value only while its parents are unchanged under the intervention (that much IS identified); if its parents change, its exogenous noise cannot be recovered from one observation and this raises — use average_causal_effect() for the population answer instead of a guessed individual one. (2) The counterfactual is relative to the network’s DAG as given: purely observational structure learning cannot orient Markov-equivalent edges (x -> y and y -> x fit equally well), so if the causal direction matters, assert it from domain knowledge rather than trusting the learned arrow.

Parameters:

• net (Any)

• observed (tuple)

• interventions (dict[int, Any])

• identification (CausalIdentification | None)

Return type:

tuple

bn_do(net, interventions, *, identification=None)

Apply graph surgery, labeling the result unidentified unless a valid receipt is supplied.

Parameters:

• net (Any)

• interventions (dict[int, Any])

• identification (CausalIdentification | None)

Return type:

InterventionalNetwork

class Forecast(mean, lo, hi, level, state_probs, samples=None, method='monte-carlo (chain-path simulation)', n_draws=0, mean_mcse=None)

Bases: object

Per-step predictive summaries plus the state-marginal trajectory.

mean, lo, and hi are MONTE CARLO estimates from n_draws predictive draws per step (STAT-RR17-10); only state_probs is exact. mean_mcse is the per-step Monte Carlo standard error of mean (draw std / sqrt(n_draws)) so two runs’ means can be read against their own noise instead of as reproducible constants – 25-draw means varied 4.39-4.72 across seeds on a fixed model. method names the estimator.

Parameters:

• mean (Any)

• lo (Any)

• hi (Any)

• level (float)

• state_probs (ndarray)

• samples (Any)

• method (str)

• n_draws (int)

• mean_mcse (Any)

forecast(model, history, horizon, *, level=0.9, n=4000, seed=0, keep_samples=False)

Forecast horizon steps beyond history under a fitted HMM.

Parameters:

• model (Any) – a fitted HiddenMarkovModelDistribution (any emission family with a sampler).

• history (Any) – the observed sequence to condition on (one sequence).

• horizon (int) – number of future steps to predict.

• level (float) – central-interval mass (0.9 -> the 5%..95% band).

• n (int) – Monte-Carlo draws per step for the emission quantiles (state marginals are exact).

• seed (int) – reproducibility.

• keep_samples (bool) – also return (H, n) JOINT sample paths – column j is one simulated trajectory of the chain, so cross-step dependence is the model’s own and summing along a column is a valid path functional (scalar emissions only).

Raises:

• TypeError – if model is not a fitted HMM.

• ValueError – for a non-exact/non-positive horizon or n, a level outside (0, 1), or an empty history.

Return type:

Forecast

Every control is checked up front, before any state is allocated or any draw is taken. Only the horizon used to be checked, and each remaining control failed somewhere else or not at all: level=-0.5 produced lo=6.75 above hi=2.25 – an interval whose lower endpoint exceeds its upper one, silently returned as an ordinary Forecast – while a NaN level surfaced as NumPy’s generic quantile error and n=0 warned about an empty sample and then raised an opaque IndexError partway through the loop.

class ParameterEstimator

Bases: ABC, Generic[SS]

Estimates a distribution from accumulated sufficient statistics.

accumulator_factory() supplies accumulators that gather sufficient statistics of type SS, and estimate(nobs, suff_stat) maps those statistics (plus optional regularization configured on the estimator) to a new distribution.

to_dict()

Return a safe JSON-compatible representation of this estimator.

Return type:

dict[str, Any]

classmethod from_dict(payload)

Reconstruct an estimator from to_dict output.

Parameters:

payload (dict[str, Any])

Return type:

ParameterEstimator

to_json(**kwargs)

Serialize this estimator as safe strict JSON.

Parameters:

kwargs (Any)

Return type:

str

classmethod from_json(text)

Deserialize an estimator from to_json output.

Parameters:

text (str)

Return type:

ParameterEstimator

abstractmethod estimate(nobs, suff_stat)

Estimate a distribution from accumulated sufficient statistics.

Parameters:

• nobs (float | None)

• suff_stat (SS)

Return type:

SequenceEncodableProbabilityDistribution

abstractmethod accumulator_factory()

Return the accumulator factory used to collect this estimator’s sufficient statistics.

Return type:

StatisticAccumulatorFactory

resident_accumulation_supported()

Return whether engine-resident (fixed-width) sufficient statistics suffice for estimate.

Most exponential-family M-steps consume only the resident sufficient statistics, so the default is True. Estimators whose M-step needs more than that (e.g. a full count histogram for the negative-binomial dispersion solve) override this to False so stacked/generated kernels fall back to the host accumulator, keeping every backend’s fixed point identical.

Return type:

bool

get_prior()

Return the parameter prior configured on this estimator, if any.

The unified estimation contract treats the prior as the single regularization concept: None gives maximum likelihood, a conjugate prior gives the Bayesian posterior update inside estimate. The default reads the prior attribute (None when unset).

Return type:

ProbabilityDistribution | None

model_log_density(model)

Return the prior log-density of model’s parameters (the ELBO global term).

Used by the variational/MAP objective in fit. The default is 0.0 (no prior); conjugate estimators override this to evaluate their prior at the model’s parameters, mapping to the prior’s parameterization first.

Parameters:

model (ProbabilityDistribution)

Return type:

float

estimate(data, estimator, prev_estimate=None)

Perform E-step in EM algorithm by iterating over all observations in ‘data’.

Arg estimator must be consistent with prev_estimate. That is, prev_estimate must be an estimate that could be obtained from estimator.

Data must type consistent with estimator and prev_estimate.

Returns the next iteration of EM algorithm by iterating over each observation of data. See seq_estimate() for a more computationally efficient implementation.

Parameters:

• data (Union[Sequence[T], pyspark.rdd.RDD]) – Sequence of iid observations of data type consistent with ‘estimator’ and/or ‘prev_estimate’.

• estimator (ParameterEstimator) – Model to be estimated from ‘data’.

• prev_estimate (Optional[SequenceEncodableProbabilityDistribution]) – Previous estimate of EM algorithm. Must be included for distributions that require initialization.

Returns:

SequenceEncodableProbabilityDistribution object.

Return type:

SequenceEncodableProbabilityDistribution

initialize(data, estimator, rng, p=0.1)

Randomly initialize a model corresponding to ParameterEstimator for iid observations data.

Note: ParameterEstimator must be of data type T, matching the input data.

This function sequentially iterates over the entire data set ‘data’, repeatedly calling initialize() method of the SequenceEncodableStatisticAccumulator object created from ‘estimator’. Data points are weighted 0 or 1 with probability p.

Seq_initialize() is much more efficient, and should produce the same initialized model for the same data sets.

Parameters:

• data (Union[Sequence[T], pyspark.rdd.RDD]) – Set of iid observations compatible with ‘estimator’.

• estimator (ParameterEstimator) – ParameterEstimator object for desired model to be estimated from data.

• rng (RandomState) – RandomState object for setting seed.

• p (float) – Proportion of data to randomly sample for initializing model.

Returns:

SequenceEncodableProbabilityDistribution object consistent with ‘estimator’.

Return type:

SequenceEncodableProbabilityDistribution

seq_estimate(enc_data, estimator, prev_estimate)

Perform vectorized E-step in EM algorithm for encoded sequence of observations in ‘enc_data’.

Arg estimator must be consistent with prev_estimate. That is, prev_estimate must be an estimate that could be obtained from estimator.

Arg enc_data must type consistent with estimator and prev_estimate (result of seq_encode() call).

Returns the next iteration of EM algorithm with vectorized calls to “seq_update()” of the corresponding SequenceEncodableStatsiticAccumulator objects.

Parameters:

• enc_data (Union[List[Tuple[int, T]], 'pyspark.rdd.RDD']) – Sequence encoded data of format matching output of seq_encode() function.

• estimator (ParameterEstimator) – Model to be estimated from ‘enc_data’.

• prev_estimate (SequenceEncodableProbabilityDistribution) – Previous estimate of EM algorithm.

Returns:

SequenceEncodableProbabilityDistribution object.

Return type:

T_D

seq_initialize(enc_data, estimator, rng, p=0.1)

Vectorized initialization of a model corresponding to ParameterEstimator for encoded sequences of iid data

observations.

Arg enc_data must type consistent with estimator (result of seq_encode() call). Arg estimator must be of data type consistent with encoded sequence data type in ‘enc_data’.

Vectorized initialization of SequenceEncodableProbabilityDistribution corresponding to ‘estimator’ from enc_data. Observations in the encoded sequence enc_data are kept with probability p.

This functions relies on calls to SequenceEncodableStatisticAccumulator.seq_initialize(), which is a vectorized initialization of the SequenceEncodableStatisticAccumulator object.

This method should produce the same initialized model as a call to initialize() if the data sets are the same.

Parameters:

• enc_data (Union[List[Tuple[int, T]], 'pyspark.rdd.RDD']) – Sequence encoded data of format matching output of seq_encode() function.

• estimator (ParameterEstimator) – Model to be estimated from ‘enc_data’.

• rng (RandomState) – RandomState object for setting seed.

• p (float) – Proportion of data to randomly sample for initializing model.

Returns:

SequenceEncodableProbabilityDistribution object consistent with ‘estimator’.

Return type:

SequenceEncodableProbabilityDistribution

optimize(data, estimator=None, max_its=10, delta=1.0e-9, init_estimator=None, init_p=0.1, rng=None, prev_estimate=None, vdata=None, enc_data=None, enc_vdata=None, out=None, print_iter=1, num_chunks=1, engine=None, precision=None, fields=None, resources=None, placement=None, sub_chunks=1, chunk_size=None, backend='local', num_workers=None, client=None, comm=None, root=0, root_only=False, strategy=None, reuse_estep_ll=True, objective='auto', on_step=None, structure='auto', schedule='full', monotone=None, track_best=None, seed=None, fused_options=None)

Fit estimator to data by a generalized-EM loop, for max_its iterations or until the

objective improves by less than delta.

Each iteration re-estimates every part of the model by whatever its structure calls for – closed-form for conjugate / exponential-family leaves, gradient descent for neural leaves, coordinate descent for GLMs, responsibility-weighted EM for latent structure (mixtures, HMMs) – so a single call fits a heterogeneous tree without the caller choosing an algorithm. (The convergence objective defaults to the family-defined maximum-likelihood objective: each leaf applies its own documented estimator update, which for some families is a closed-form moment update rather than an iterative likelihood maximization, and families may apply documented numerical floors/repairs – read them back via model.numerical_repairs(). A parameter prior switches the objective to penalized-LL / MAP, and a variational model to the ELBO – see objective.)

Missing values. With estimator=None (the default) auto-inference finds the gaps and wraps the affected leaf for you: data carrying None, NaN, pd.NA or pd.NaT fits an OptionalDistribution whose missingness rate is estimated alongside the base family. When you build the estimator yourself, that wrapper is yours to add, and it must name the sentinel your data actually carries – the wrapper matches its missing_value by identity, so the None default does not accept NaN:

optimize(rows, CompositeEstimator([OptionalEstimator(GaussianEstimator(),
                                                     missing_value=float("nan")),
                                   CategoricalEstimator()]))

NaN is the spelling float arrays and numeric pandas columns use, so missing_value=float('nan') is the one that works for tabular numeric data; missing_value=None (the default) is for records that carry None. mixle.stats.marginalized() is the same wrapper with no fitted rate, but it takes a DISTRIBUTION rather than an estimator and its default sentinel is not NaN, so NaN data must spell out marginalized(dist, missing_value=float('nan')).

A native +inf/-inf gets the same auto-inference treatment as the sentinels above – fitted as its own OptionalDistribution(..., missing_value=math.inf) (one wrapper per sign present) rather than rejected – because an infinity is not itself a value most continuous families can carry, so treating it as “a measurement that comes out infinite p of the time” keeps the auto-built model a proper generative law instead of refusing to fit. This is DIFFERENT from the None/NaN/pd.NA/pd.NaT case above: those are absences by convention; an inf is a legitimate (if unusual) numeric value, so estimator=None raises a UserWarning naming the affected field(s) when it reclassifies one – catch it with warnings.catch_warnings() if your pipeline needs to notice. Passing an estimator explicitly bypasses auto-inference and gets the family’s own contract instead: a family whose support excludes inf (e.g. GaussianEstimator) raises ValueError rather than silently wrapping it.

Identifier-like columns. estimator=None auto-inference (via mixle.utils.automatic.get_estimator) also freezes a scalar column it cannot otherwise type – an unrecognized scalar type (datetime, Decimal, …), an ambiguous bool/number mix, and an identifier-like/high-cardinality column (e.g. a user id, order id, or other near-unique key) – as the empirical categorical over the values observed at fit time, held fixed by IgnoredDistribution. A value not seen at fit time then scores -inf there, the same finite-support behavior every automatically fitted categorical has (CategoricalDistribution’s default_value docstring); this is silent at fit time – no warning, unlike the +inf/-inf reclassification above – because freezing is the correct, unsurprising answer for a column of this shape, not a defect being flagged. It matters at SCORING time: a held-out or production row carrying a fresh identifier scores that field at -inf, and because log-densities sum, one such field silently drives the WHOLE model’s log_density/Model.evaluate() score to -inf even though every other field scored normally (Model.evaluate() warns rather than raising on this, precisely because it is this documented, by-design behavior rather than a scoring bug – see its docstring). mixle.propose() fits the identical policy but additionally excludes just the affected field(s) from a candidate’s held-out verification score (with a note naming them, in Model.notes) rather than voiding the whole candidate – optimize() itself does not do this rescue, so an identifier-like field’s presence is worth checking for explicitly (e.g. via model.log_density per field, or mixle.utils.automatic.get_estimator(data) inspection) before trusting a held-out score built on a model this call produced.

DataFrames. optimize and fit take a pandas DataFrame or Series directly; the row-level stats API does not. mixle.stats.seq_encode consumes RECORDS and refuses a frame (“expected a sequence of 2-tuples, got DataFrame”) – converting is the caller’s job, in one call: mixle.data.dataframe_records(df) for the records, or mixle.data.seq_encode_dataframe(df, model=...), which is the frame-shaped counterpart of seq_encode. A frame’s gaps are canonicalized on the way in by what each column HOLDS – NaN for a column of numbers, None for anything else – so df and df.convert_dtypes() fit models carrying the same missing_value, and each model scores the other frame.

Parameters:

• data (Optional[List[T]]) – List of data type T containing observed data. Must be compatible with data type of estimator.

• estimator (ParameterEstimator | ProbabilityDistribution | None) – What to fit. A ParameterEstimator is used as-is; a distribution prototype (any ProbabilityDistribution) is coerced to its matching estimator via proto.estimator() so you build the model shape only once; None infers an estimator from raw data (mixle.utils.automatic.get_estimator). Inferring from rows that do not all carry the same number of fields is ambiguous between a table with a malformed row and variable-length sequence data, so it is decided on the arity evidence rather than guessed: an overwhelmingly dominant arity raises ContractError naming the offending row, a bare majority is fitted as a sequence and warns that it did, and widely spread arities are fitted as a sequence in silence. Pass mixle.utils.automatic.get_estimator(data, ragged='sequence') as the estimator to demand the sequence reading outright.

• max_its (int) – Maximum number of EM iterations to be performed. Default value is 10 iterations.

• delta (Optional[float]) – Stopping criteria for EM algorithm used if max_its is not set: Iterate until abs(old_loglikelihood - new_loglikelihood) < delta or iterations == max_its. delta=None requests a fixed max_its iteration count by disabling this convergence test – but it does NOT disable the separate acceptance gate every iteration still passes through (a proposed step that does not improve the objective, or produces a non-finite one, is rejected and ends the loop regardless of delta). That rejection can still cut a delta=None run short of max_its, from ordinary float noise near a fixed point; a UserWarning discloses it when it happens (compare fit_provenance().iterations against max_iterations to check yourself instead).

• init_estimator (Optional[ParameterEstimator]) – ParameterEstimator to used to initialize EM algorithm parameters. If None, estimator is used. Must be consistent with estimator.

• init_p (float) – Value in (0.0,1.0] for randomizing the proportion of data points used in initialization. This is a statistical knob, not an execution detail: the initialization is the estimator’s answer on the drawn subsample, and for latent models and for families whose documented update is approximate rather than an exact likelihood ascent (e.g. Weibull, GeneralizedExtremeValue, WrappedCauchy) the accepted EM trajectory – and therefore the returned parameters – depends on that start. In the extreme, such a family’s very first update can fail the monotone acceptance gate, and the fit then returns the initialization itself (fit_provenance() reports converged=False with iterations below the cap). init_p=1.0 initializes from the estimator’s answer on all the data.

• rng (RandomState) – RandomState used to set seed for initializing EM algorithm. None resolves to a FIXED seed, so the NumPy-driven parts of an un-seeded optimize/fit (initialization, EM, subsampling) are deterministic by default; pass your own RandomState when you WANT different initializations across calls (e.g. hand-rolled restarts). Torch-backed leaves are the deliberate exception: modules that consume torch’s global RNG (dropout, VAE reparameterization draws, minibatch shuffling) follow torch’s own default non-determinism – call torch.manual_seed yourself when a torch-backed fit must be exactly reproducible. An integer is accepted and coerced to RandomState(rng). Mutually exclusive with seed.

• vdata (Optional[Sequence[T]]) – Optional validation set.

• prev_estimate (Optional[SeqeuenceEncodableProbabilityDistribution]) – Optional model estimate used from prior fitting. Must be consistent with estimator.

• enc_data (Optional[List[Tuple[int, E]]]) – Optional encoded data of form List[Tuple[int, E]]. Formed from data if None.

• enc_vdata (Optional[List[Tuple[int, E0]]]) – Optional sequence encoded validation set.

• out (IO | None) – Stream for per-iteration EM progress lines. Defaults to None (quiet, so the library does not spam stdout in normal use); pass out=sys.stdout to watch convergence.

• print_iter (int) – Print the log-likelihood difference every print_iter iterations; the final converged iteration is always reported. Pass print_iter=0 to suppress the periodic lines (keeping only the converged line), or out=None to silence entirely.

• num_chunks (int) – Number of chunks for encoded data. For exact-sufficient-statistic leaves the chunk statistics combine exactly, so chunking changes only float summation order. Chunking also changes WHICH rows the init_p initialization subsample draws, however, so for the initialization-sensitive families described under init_p a different num_chunks can change the fitted parameters, not just the execution plan – measured at the percent level on Weibull/GEV/WrappedCauchy fits. When varying it, compare fit_provenance().final_objective (or hold the start fixed with init_p=1.0).

• engine (Optional[Any]) – Optional ComputeEngine for local kernel scoring/accumulation. Distributed engine placement is intentionally deferred to the orchestrator/planner layer.

• precision (Optional[Any]) – Optional floating-point precision such as 'float32' or np.float64. Pass 'auto' to let mixle.engines.auto_precision choose from the data and engine: float32 only on a GPU torch engine with well-conditioned numeric data, else float64. Pass 'minimal' for the data-aware CPU allocator (mixle.inference.precision_plan): it inspects the data magnitude and the model’s leaf families/conditioning and runs the reduced float32 fused kernel where verified safe (accumulation stays float64), else float64 – the “preserve accuracy with minimal compute” default for local fits.

• fields (Optional[Any]) – DataFrame column/field selection. A single field yields scalar observations; several fields yield tuple observations unless the estimator/model is record-shaped, in which case dict records are produced by source column name.

• resources (Optional[Any]) – Optional planner resources. When supplied with raw data, optimize encodes through the shared encoded-data factory so placement, sub-chunks, and per-shard engines use the orchestrator contract.

• placement (Optional[Any]) – Optional explicit placement produced by mixle.utils.parallel.planner.plan.

• sub_chunks (int) – Number of sub-chunks per placement shard when resources or placement is supplied.

• chunk_size (Optional[int]) – Approximate chunk size for ordinary local sequence encoding.

• backend (str) – Encoded-data backend for raw data. 'local' keeps the historical local encoding unless resources/placement are supplied; 'mp' and 'mpi' use the shared encoded-data factory.

• num_workers (Optional[int]) – Worker count for backend='mp' and optional partition count hint for backend='dask'. The default is the CPU count capped at the number of observations; an explicit value is used as given (capped only at the number of observations), so a value far above the core count spawns that many processes and fails slowly at the worker setup timeout rather than promptly. Keep it at or below the core count.

• client (Optional[Any]) – Existing dask.distributed client for backend='dask'. If omitted, the dask backend uses an active default client or starts a local threaded client.

• comm (Optional[Any]) – MPI communicator for backend='mpi'.

• root (int) – MPI root rank for backend='mpi'.

• root_only (bool) – MPI root-only data mode for backend='mpi'.

• strategy (Optional[Any]) – Optional EM strategy from mixle.inference.em (e.g. AnnealedEM, HardEM, MonteCarloEM, CompiledEM) or any callable (enc, estimator, model) -> model to use in place of the standard exact E/M step. None uses the standard step.

• reuse_estep_ll (bool) – Default True. Reuse the data log-likelihood computed during the E-step (the posterior normalizer / forward pass / variational ELBO) for convergence instead of running a separate scoring pass each iteration – typically ~1.5-2x faster per iteration for latent models (mixtures, HMMs and variants, topic models, associations, IBP, …) on the default local engine. Convergence then lags by one iteration (same fixed point) and the best-likelihood model is returned; fixed-iteration fits (delta=None) are identical to the standard loop. Automatically falls back to the standard loop for engines/strategies/ distributed backends or models that can’t report the LL (no slowdown there). Set False to force the exact historical per-iteration scoring behavior.

• objective (str) – Convergence/selection objective. 'auto' (default) makes the prior the single switch – a model exposing a variational ELBO (seq_local_elbo) is fit by variational Bayes ('vb'), an estimator carrying a parameter prior by penalized log-likelihood ('map'), and everything else by maximum likelihood ('mle'). 'mle' is the family-defined maximum-likelihood objective: what is guaranteed is each leaf family’s documented estimator update – exact MLE for exponential-family leaves, a documented closed-form moment update for families whose class docstring says so – and families may apply documented numerical floors/repairs, reported on the fitted model via numerical_repairs(). Pass 'mle' / 'map' / 'vb' to force a specific objective. fit accepts the same argument; both share this resolution so a Bayesian estimator is fit on the correct objective regardless of the verb used. (Only 'mle' is compatible with the fused E-step shortcut; reuse_estep_ll is ignored for 'map'/'vb'.)

• monotone (Optional[bool]) –

  Outer-objective acceptance policy. None (default) uses strict generalized-EM acceptance for immutable closed-form updates and best-seen selection for mutable neural or variational/approximate updates. In best-seen mode finite downhill steps may be traversed, but the returned model (including mutable module parameters) is restored to the best selected-objective value observed. Pass True to reject the first decreasing step, or False to permit a non-monotone trajectory explicitly.

  Convergence contract (worklist Q5.4; pinned by em_convergence_contract_test): under strict acceptance the accepted-round objectives are non-decreasing within tolerance and – bounded above by the estimator variance floors – the objective sequence converges; limit-point stationarity is the classical EM/GEM theory per family. Under best-seen selection the guarantee is the best visited iterate; a neural leaf fit with lr_decay in (0.5, 1] additionally follows a Robbins–Monro step schedule (the condition stochastic-approximation EM analyses require). A non-finite objective is never accepted and never becomes the convergence reference.

• track_best (Optional[bool]) – Whether to restore the best outer-objective state seen. None (default) does so except when an estimator explicitly declares a surrogate fitting objective, such as NCE; observed density is not a valid selector until such a model is normalized, so its final finite update is returned. Pass a boolean to override selection.

• on_step (Optional[Callable[[EMStep], None]]) – Optional per-iteration callback receiving an EMStep (iter, model, log_density, delta) for the accepted model. Use it to checkpoint a long run – e.g. on_step=registry.checkpointer('run', every=5) – and resume with prev_estimate=. Called on every iteration regardless of print_iter.

• structure (str) – 'auto' (default) makes the tagline literal for flat tuple records fit with no estimator: the cross-field dependency graph is discovered (mixle.inference.learn_bayesian_network()) and returned when it beats the independent composite by BIC — otherwise (no edges, non-record data, or any failure) the historical automatic-composite path proceeds untouched. A pandas DataFrame is converted to the same flat records first (honoring fields), so a table gets the same structure inference whether it arrives as a DataFrame or as a list of row tuples. 'off' restores the unconditional historical behavior. Only consulted when estimator is None and no prev_estimate/init_estimator/strategy/enc_data is supplied.

• schedule (str) – 'full' (default) performs exact full-tree EM every round and automatically uses whole-model or component-level compiled kernels when the local model is eligible. 'auto' additionally engages the block-coordinate-ascent scheduler (mixle.inference.block_em) when selective work is expected to win: after one full bootstrap sweep, components are ranked by their last observed complete-data Q gain per measured whole-block cost, and only the highest-value ones within a per-round budget are re-estimated. Inactive component scores are cached while their parameters remain unchanged. This is a scheduling choice: observed likelihood is transactionally gated non-decreasing every round, and when there is no useful ranking to do (e.g. every component looks equally worth updating) the scheduler degenerates to doing exactly what 'full' does. Only engaged when the model is a local-backend MixtureDistribution/ MixtureEstimator MLE/MAP fit with no explicit strategy/engine/resources/ placement – anything else silently falls back to 'full' (never an error, never a behavior change beyond scheduling).

• seed (Optional[int]) – Integer seed for initializing the EM algorithm – shorthand for rng=RandomState(seed), matching the seed= argument the samplers and the other entry points take. Mutually exclusive with rng (passing both raises TypeError).

• fused_options (Optional[dict]) – Tuning knobs for the fused compiled kernels, applied when the fit runs on a fused engine (an explicit fused engine= or the auto-fusion gate). Recognized keys: parallel (bool – force the chunk-parallel kernels on or off, overriding the observation-count auto-gate; results are bit-stable either way), lse_bits (int – opt SCORING into quantized log-sum-exp with a ~2**-bits relative bound; E-steps always stay exact), and lse_span (float, default 24.0 – the clipped LSE exponent range used with lse_bits). Unknown keys raise ValueError. Ignored (with the same validation) when the fit resolves to a non-fused path.

Returns:

the fitted model. The run behind it ended in one of

three ways – the objective gain fell below delta (a converged fit), the max_its cap was reached with the objective still improving (an unconverged fit; a UserWarning says so when delta is in force), or the family’s next update failed the monotone acceptance gate before the cap (no further progress is possible; see init_p). The receipt distinguishing them ships on the model: model.fit_provenance() reports converged, iterations, final_objective (the returned model’s objective value), and any numerical repairs.

Return type:

SequenceEncodableProbabilityDistribution

Model comparison. A fitted core distribution intentionally carries no aic()/bic() – the core cannot count free parameters exactly for every composed model. The supported paths: per-observation log-densities (model.log_density(x)) feed the paired comparison tests exported by mixle.inference (vuong_test(), clarke_test(), paired_score_difference(), compare_elpd()), and the mixle.ppl surface fits the same families with a recorded parameter dimension, giving m.aic(data) / m.bic(data) and best-first ranking via mixle.ppl.compare([m1, m2], data, by='bic').

jit_seq_log_density(model, engine=None)

Return a JittedScorer: the whole-tree log-density of model compiled to one XLA program.

jit_seq_log_density(model)(data) is bit-identical to model.seq_log_density(encode(data)) but runs as a single jax.jit XLA program over the entire composite tree – fast for repeated scoring of a fixed model. Requires the JAX optional extra and that every leaf in model declares JAX support (raises EngineNotSupportedError-style errors from the engine layer otherwise).

Parameters:

• model (Any)

• engine (Any)

Return type:

JittedScorer

jit_em_mixture(model, data, *, max_its=100, engine=None, return_receipt=False)

Fit a finite mixture of same-family scalar exponential-family leaves by EM, with the ENTIRE EM loop (every E-step + closed-form weighted M-step iteration) compiled to ONE jax.jit XLA program via lax.scan – the parameters are traced inputs threaded through the loop on-device, so there is no per-iteration recompile and no per-iteration host sync. This is roadmap A2 bullet 1 (“repeated EM iterations run as one XLA program”) realized literally.

model is the initial mixture (its components seed the EM); supported leaves: Gaussian, Poisson, Exponential. Runs a fixed max_its iterations (no host-side early stop – that is the point: the loop stays on-device). Returns a fitted MixtureDistribution, bit-close to the host EM from the same start (it is the same EM update). Raises NotImplementedError for unsupported structure.

SPEED – measured scope: the payoff is GPU/TPU and large scale, where XLA parallelizes the E-step over millions of points and many components. On an Apple M4 GPU (via jax-metal) this kernel runs ~21x faster than mixle’s vectorized NumPy EM (K=10, N=1e6, 50 iters: 156 ms vs 3254 ms) with identical estimates, and ~12x faster than the same jitted loop on CPU. On CPU it is *not* a speedup – mixle’s host EM is already vectorized NumPy (+ a fused path) and a long sequential scan of small steps loses to it (~0.1-0.8x for K up to 40). So: GPU -> big win, CPU -> use the host EM. (The CPU win from A2 is the single-pass scoring jit, jit_seq_log_density(), ~8x.) Note: jax-metal is version-pinned – the GPU result above used Python 3.11 + jax/jaxlib 0.4.34 + jax-metal 0.1.1; newer jaxlib emits StableHLO that jax-metal 0.1.1 cannot compile.

Parameters:

• model (Any)

• data (Any)

• max_its (int)

• engine (Any)

• return_receipt (bool)

class JITExecutionResult(model, receipt)

Bases: object

A fitted model paired with evidence of the JIT execution that actually ran.

Parameters:

• model (Any)

• receipt (dict[str, Any])

class JittedScorer(model, engine=None)

Bases: object

A jax.jit-compiled whole-tree log-density scorer for a fixed model.

Calling the scorer encodes data, runs the engine-neutral backend_seq_log_density over the whole model tree on the JAX engine under jax.jit, and returns host log-densities. The compiled program is cached and reused across calls with the same data shape.

Parameters:

• model (Any)

• engine (Any)

laplace_posterior(model, data, *, log_prior=None, eps=1e-4, ridge=1e-6, mode_tol=1e-3)

Build a validated Laplace approximation around model’s parameters.

log_prior is evaluated in the flattened unconstrained coordinates. If supplied, the target is the log posterior; otherwise the result is explicitly a likelihood-curvature approximation, despite the compatibility name of this function and result class.

Parameters:

• log_prior (Callable[[ndarray], float] | None)

• eps (float)

• ridge (float)

• mode_tol (float)

Return type:

LaplacePosterior

class LaplacePosterior(mode_model, u_mode, cov, rebuild, metadata)

Bases: object

Validated Gaussian Laplace approximation in unconstrained parameter space.

The compatibility name is retained, but is_posterior is true only when the target included an explicit prior.

sample(n=1, rng=None)

Draw model samples from the Laplace Gaussian approximation. rng may be None (unseeded), an int seed, or an existing RandomState.

Parameters:

n (int)

summary()

Return Laplace approximation metadata.

Return type:

dict

collapse_mixture(mixture)

Moment-match a mixture onto a single distribution, in closed form (exact for Gaussians).

Returns the Gaussian minimizing KL(mixture || Gaussian) – its mean and covariance are the mixture’s overall mean and covariance (law of total variance). This is the exact, sample-free counterpart to mixle.ops.project(mixture, GaussianDistribution(...).estimator()).

Parameters:

mixture (Any)

Return type:

Any

reduce_mixture(mixture, n_components, *, method='runnalls')

Reduce a Gaussian mixture to n_components by greedily merging the least-costly pair (Runnalls).

Repeatedly merges the two components whose merge costs the least KL (_runnalls_cost()) until n_components remain. Every merge is moment-preserving, so the reduced mixture has the same overall mean and covariance as the original – only higher moments are lost. Returns a mixture of the same kind as the input (a GaussianMixtureDistribution for multivariate input).

Parameters:

• mixture (Any) – a Gaussian mixture (GaussianMixtureDistribution or a mixture of Gaussian components).

• n_components (int) – target number of components (>= 1); no-op if already <= the current count.

• method (str) – only "runnalls" for now.

Return type:

Any

moment_project(teacher, target=None, *, exact=True, **sampling_kw)

Project teacher onto a smaller student – exactly when possible, else by sampling.

If teacher is a Gaussian mixture and target is None (or a single Gaussian family), the projection is the closed-form collapse_mixture() – no samples, machine-precision. Otherwise (or when exact=False) it delegates to mixle.ops.project(), the sampling M-projection onto target’s family. This gives one entry point that is exact where the structure allows and explicit about sampling where it does not.

Parameters:

• teacher (Any)

• target (Any)

• exact (bool)

• sampling_kw (Any)

Return type:

Any

gaussian_kl(p, q)

KL(p || q) between two Gaussians (uni- or multivariate), in nats. Accepts distributions.

0.5[ tr(Σq⁻¹ Σp) + (μq-μp)ᵀ Σq⁻¹ (μq-μp) - d + ln(det Σq / det Σp) ] – the analytic Gaussian KL.

Parameters:

• p (Any)

• q (Any)

Return type:

float

class TrialsToTarget(n_trials, final_model, history=<factory>)

Bases: object

How many trials a path needed to reach its target (None = never reached), plus its score history.

Parameters:

• n_trials (int | None)

• final_model (Any)

• history (list[float])

apply_add_edge_fix(model, fault, data)

Apply only the fix diagnose actually named: add a linear-Gaussian edge between the two fields in fault.dominant (parsed from its "field[i]|...field[j]|..." shape), refit that one factor on data, and return the corrected network – every other factor is left untouched.

Returns None (never guesses) when fault.suggested_fix is not "add_edge", or dominant does not name exactly two fields (nothing dominant, or a shape this translation doesn’t understand).

fault.dominant names an undirected pair – it carries no information about which field is the parent. Both orientations are fit and the better-fitting one on data is kept: field-index order (which field happens to have the smaller schema index) has no relationship to causal direction, so always taking parent, child = sorted(idx) risked silently adding the edge backwards.

“Add edge” means ADD. The candidate keeps the child’s existing parents, its categorical encoding and vector metadata, and is only built where doing so preserves the child’s factor family; an orientation that would instead rewrite the child’s model, or that would close a cycle, is reported as unsupported rather than performed under this function’s name. None is returned when neither orientation is a semantics-preserving edit.

Parameters:

• model (HeterogeneousBayesianNetwork)

• fault (FaultReport)

• data (Sequence[tuple])

Return type:

HeterogeneousBayesianNetwork | None

blind_search_trials_to_target(data, held_out, target_score, *, round_size=10, max_rounds=20, max_parents=1, seed=0)

The blind baseline: NO diagnosis. Re-run learn_bayesian_network from scratch on growing prefixes of data (round_size more examples each round, shuffled once up front) until the discovered structure’s held-out score reaches target_score (typically directed_correction()’s own held-out score) or max_rounds is exhausted.

Parameters:

• data (Sequence[tuple])

• held_out (Sequence[tuple])

• target_score (float)

• round_size (int)

• max_rounds (int)

• max_parents (int)

• seed (int)

Return type:

TrialsToTarget

directed_correction(model, cases, data, held_out, *, background=None)

One diagnosis, one targeted edit, one verification.

Costs exactly 1 trial if diagnose names an actionable fix and it verifiably improves held-out score over the original model; costs 0 (unreached) if the fix isn’t actionable or doesn’t verify – a correction that does not improve held-out is a failed diagnosis, logged as such, never silently kept anyway.

A candidate is only “verified” against a comparison that actually happened. held_out must be nonempty, and both totals must be comparable scores: a NaN (scoring failed) fails after <= before just as a real improvement does, so without this check a model that could not be scored at all is returned as a successful one-trial correction.

Raises:

ValueError – if held_out is empty – there is nothing to verify against.

Parameters:

• model (Any)

• cases (Sequence[tuple])

• data (Sequence[tuple])

• held_out (Sequence[tuple])

• background (Sequence[tuple] | None)

Return type:

TrialsToTarget

held_out_log_likelihood(model, data)

Total log-density of a fitted network (or anything exposing dist_to_encoder/seq_log_density) over held-out data – the metric both paths below are compared on.

Parameters:

• model (Any)

• data (Sequence[tuple])

Return type:

float

fisher_merge(estimates, fishers=None)

Fisher-weighted merge of parameter estimates – the closed-form Laplace-posterior combination.

Given estimates θ_i (each a flat parameter vector) and their Fisher information F_i, returns θ* = (Σ F_i)⁻¹ (Σ F_i θ_i) – the point that maximizes the sum of the local Laplace log-posteriors Σ_i -½(θ-θ_i)ᵀ F_i (θ-θ_i). This is Matena & Raffel Fisher merging (diagonal F) and, in general, the precision-weighted mean; for Gaussians it coincides with the product-of-experts mean. No gradient steps – a single linear solve.

Parameters:

• estimates (Any) – sequence of parameter vectors θ_i (each shape (p,)), or a stack (m, p).

• fishers (Any) – per-estimate Fisher information. Each may be a scalar/1-D vector (diagonal Fisher, per-coordinate precision) or a (p, p) matrix (full Fisher). None uses unit Fisher (a plain average). A single value is broadcast to every estimate.

Returns:

The merged parameter vector θ* of shape (p,).

Return type:

ndarray

fit(data, estimator=None, max_its=10, delta=1.0e-6, init_estimator=None, **kwargs)

Fit a model in the Bayesian (variational / MAP) sense, returning the posterior-bearing model.

This is the posterior-returning counterpart of optimize(). fit iterates the EM/VB update that maximizes the objective selected by objective (default 'auto'):

• 'auto' – the prior is the single switch: 'vb' when the model exposes seq_local_elbo, 'map' when the estimator carries a parameter prior, else 'mle';

• 'mle' – the family-defined maximum-likelihood objective on the data log-likelihood (ignores any prior in the objective). What is guaranteed is each leaf family’s documented estimator update – some families document a closed-form moment update rather than an iterative likelihood maximization – and families may apply documented numerical floors/repairs; read them back via model.numerical_repairs();

• 'map' / 'vb' – penalized log-likelihood / ELBO obj = data term + prior term, where the data term is the observed-data LL (MAP) or local-ELBO contributions (variational), and the prior term is estimator.model_log_density(model).

Convergence is checked on the chosen objective, so under 'map'/'vb' the prior is part of the stopping rule and conjugate updates never decrease it. The returned model carries its conjugate posterior forward as model.get_prior(). With no prior anywhere, every objective reduces to plain EM, so fit and optimize agree for frequentist estimators. optimize accepts the same objective argument; the two share this resolution so a Bayesian estimator is fit correctly regardless of which verb the caller reaches for.

Args otherwise mirror optimize() (local encoded path). Returns the model with the best validation log-likelihood seen during the run.

fit is a thin wrapper over optimize() – they share the one EM/objective loop. fit adds only the opt-in data-structure check, a Bayesian-leaning default delta (1e-6), and the exact per-iteration-scored loop (reuse_estep_ll=False). Every other optimize() keyword – engines, precision, distributed backend, on_step, the fused E-step – is accepted here too and forwarded verbatim, so reaching for a heavier knob never means switching verbs. estimator accepts the same three spellings as optimize() (estimator, distribution prototype, or None to infer from data).

Parameters:

• data (Sequence[T] | None)

• estimator (ParameterEstimator | ProbabilityDistribution | None)

• max_its (int)

• delta (float | None)

• init_estimator (ParameterEstimator | ProbabilityDistribution | None)

• kwargs (Any)

Return type:

SequenceEncodableProbabilityDistribution

class ObjectiveFitResult(model, value, iterations, history=(), converged=False, initial_value=None, final_delta=None, maximize=True, best_value=None, best_iteration=None, final_gradient_norm=None)

Bases: object

Optimization result for arbitrary differentiable objectives.

By default objective helpers restore the best trainable state seen during optimization. In that mode value and model refer to best_value / best_iteration while history still records every attempted iterate and final_delta describes the last attempted step.

Parameters:

• model (Any)

• value (float)

• iterations (int)

• history (tuple[float, ...])

• converged (bool)

• initial_value (float | None)

• final_delta (float | None)

• maximize (bool)

• best_value (float | None)

• best_iteration (int | None)

• final_gradient_norm (float | None)

as_tuple()

Return the historical (model, value) shape used by fit helpers.

Return type:

tuple[Any, float]

property objective_change: float | None

Return the signed objective change from the start of optimization.

property improvement: float | None

Return positive improvement in the requested optimization direction.

property best_improvement: float | None

Return best positive improvement seen during optimization.

class ObjectiveParameter(name, value, constraint='real')

Bases: object

Named trainable parameter for arbitrary differentiable objectives.

Supported constraints are real, positive / positive_vector / positive_matrix, unit_interval, simplex / simplex_vector, row_simplex_matrix, column_simplex_matrix, greater_than:<name>, and less_than:<name>. Coupled bound constraints refer to an earlier parameter in the same set. The optimizer stores unconstrained raw tensors and presents constrained values to the objective callable.

Parameters:

• name (str)

• value (Any)

• constraint (str)

class ObjectiveParameterSet(parameters, engine=None, precision=None, torch=None)

Bases: object

Engine-backed named parameters for user-supplied objectives.

Parameters:

• parameters (Any)

• engine (Any | None)

• precision (Any | None)

• torch (Any | None)

trainable_tensors()

Return raw tensors passed to the optimizer.

Return type:

Sequence[Any]

values()

Return constrained tensors keyed by parameter name.

Return type:

Mapping[str, Any]

detached_values()

Return plain Python/NumPy constrained values keyed by parameter name.

Return type:

Mapping[str, Any]

class ExpectedLogDensity(weights=None, normalize=False)

Bases: object

Objective sum_i w_i log q_model(x_i) for encoded observations.

Parameters:

• weights (Sequence[float] | None)

• normalize (bool)

class ObjectiveSum(*objectives)

Bases: object

Add several model objectives into one scalar objective.

Parameters:

objectives (ObjectiveCallable)

class UnnormalizedLogLikelihood(log_unnormalized, log_partition=None, partition_enc=None, reference_log_density=None, weights=None, normalize=False)

Bases: object

Objective for models specified by unnormalized log likelihoods.

The objective is

sum_i w_i * log f_theta(x_i) - sum_i w_i * log Z(theta)

where log f_theta is supplied by the caller. log Z(theta) can be supplied exactly with log_partition or estimated from reference samples using self-normalized importance form logmeanexp(log f_theta(y_j) - log q(y_j)).

Parameters:

• log_unnormalized (ObjectiveCallable)

• log_partition (Callable[[Any, Any], Any] | None)

• partition_enc (Any | None)

• reference_log_density (Callable[[Any, Any], Any] | None)

• weights (Sequence[float] | None)

• normalize (bool)

class CallableObjective(fn, name='callable_objective')

Bases: object

Small adapter naming arbitrary objective callables.

Parameters:

• fn (ObjectiveCallable)

• name (str)

fit_objective(enc, model, objective, engine=None, max_its=500, lr=0.05, optimizer='adam', tol=1.0e-7, maximize=True, out=None, print_iter=100, return_result=False, precision=None, restore_best=True)

Optimize a user-supplied differentiable objective over a distribution tree.

The objective is called as objective(shadow_model, enc, engine) and must return a scalar engine tensor. Distribution parameters are obtained from the same declaration-driven raw state used by fit_mle / fit_map. When restore_best is true, the returned model/value correspond to the best objective value seen, not necessarily the last attempted optimizer step.

Parameters:

• enc (Any)

• model (SequenceEncodableProbabilityDistribution)

• objective (Callable[[Any, Any, Any], Any])

• engine (Any | None)

• max_its (int)

• lr (float)

• optimizer (str)

• tol (float)

• maximize (bool)

• out (Any | None)

• print_iter (int)

• return_result (bool)

• precision (Any | None)

• restore_best (bool)

Return type:

Any

fit_parameter_objective(parameters, objective, enc=None, engine=None, max_its=500, lr=0.05, optimizer='adam', tol=1.0e-7, maximize=True, out=None, print_iter=100, precision=None, return_result=False, restore_best=True)

Optimize an arbitrary objective over named constrained parameters.

parameters may be a mapping of name -> initial_value for real parameters, a sequence of ObjectiveParameter objects, or a sequence of (name, initial_value, constraint) tuples. objective is called as objective(params, enc, engine), where params is a mapping of constrained engine tensors. When restore_best is true, returned detached values are the best seen named parameters.

Parameters:

• parameters (Any)

• objective (Callable[[Mapping[str, Any], Any, Any], Any])

• enc (Any | None)

• engine (Any | None)

• max_its (int)

• lr (float)

• optimizer (str)

• tol (float)

• maximize (bool)

• out (Any | None)

• print_iter (int)

• precision (Any | None)

• return_result (bool)

• restore_best (bool)

Return type:

Any

optimize_torch_objective(parameters, objective, engine=None, max_its=500, lr=0.05, optimizer='adam', tol=1.0e-7, maximize=True, out=None, print_iter=100, precision=None, return_result=False, restore_best=True)

Optimize an arbitrary Torch objective over supplied tensor parameters.

This is the escape hatch for models whose likelihood is not an iid distribution score, such as Gaussian-process marginal likelihoods or supervised neural-network losses. When restore_best is true, supplied tensors are copied back to their best seen values before returning.

Parameters:

• parameters (Iterable[Any])

• objective (Callable[[], Any])

• engine (Any | None)

• max_its (int)

• lr (float)

• optimizer (str)

• tol (float)

• maximize (bool)

• out (Any | None)

• print_iter (int)

• precision (Any | None)

• return_result (bool)

• restore_best (bool)

Return type:

Any

projection_samples(source, sample_size, seed=None)

Draw reusable samples for Monte-Carlo projection experiments.

Parameters:

• source (SequenceEncodableProbabilityDistribution)

• sample_size (int)

• seed (int | None)

Return type:

Sequence[Any]

variational_projection(source, target, data=None, enc=None, sample_size=1000, seed=None, engine=None, max_its=500, lr=0.05, optimizer='adam', tol=1.0e-7, out=None, print_iter=100, return_result=False, precision=None, restore_best=True)

Project source onto the family represented by target.

This minimizes a Monte-Carlo estimate of the forward KL KL(source || target) by maximizing E_source[log target(X)]. Pass data/enc to reuse common random numbers or a curated design set; otherwise samples are drawn from source.

Parameters:

• source (SequenceEncodableProbabilityDistribution)

• target (SequenceEncodableProbabilityDistribution)

• data (Sequence[Any] | None)

• enc (Any | None)

• sample_size (int)

• seed (int | None)

• engine (Any | None)

• max_its (int)

• lr (float)

• optimizer (str)

• tol (float)

• out (Any | None)

• print_iter (int)

• return_result (bool)

• precision (Any | None)

• restore_best (bool)

Return type:

Any

class EMStep(iter, model, log_density, delta)

Bases: NamedTuple

One accepted EM iteration, handed to an optimize(on_step=...) callback.

iter is 1-based; model is the current accepted model – snapshot it to checkpoint, and resume with optimize(prev_estimate=...); log_density is the training objective at this step; delta is its gain over the previous step (inf on the first iteration).

Parameters:

• iter (int)

• model (Any)

• log_density (float)

• delta (float)

iter: int

Alias for field number 0

model: Any

Alias for field number 1

log_density: float

Alias for field number 2

delta: float

Alias for field number 3

best_of(data, vdata, est, trials, max_its, init_p, delta, rng=None, init_estimator=None, enc_data=None, enc_vdata=None, out=None, print_iter=1, reuse_estep_ll=True, objective='auto', seed=None, fused_options=None)

Performs EM algorithm for trials-number of randomized initial conditions. Returns the best model fit in terms of

maximum log-likelihood value from validation data.

Parameters:

• data (Optional[List[T]]) – List of data of type T. If None is given, enc_data must be provided as List[Tuple[int, enc_data_type]].

• vdata (Optional[Sequence[T]]) – Optional validation set.

• est (ParameterEstimator) – ParameterEstimator for model to be estimated.

• trials (int) – Integer number >= 1, of randomized initial conditions to perform EM algorithm for.

• max_its (int) – Integer value >=1, sets the maximum number of iterations of EM to be performed as stopping criteria.

• init_p (float) – Value in (0.0,1.0] for randomizing the proportion of data points used in initialization.

• delta (float) – Stopping criteria for EM when abs(old-log-likelihood - new-log-likelihood) < delta.

• rng (RandomState) – RandomState for setting seed. An integer is coerced to RandomState(rng); None (default) resolves to the fixed default seed. Mutually exclusive with seed.

• init_estimator (Optional[ParameterEstimator]) – Optional ParameterEstimator used for fitting.

• enc_data (Optional[List[Tuple[int, E]]]) – Optional encoded data, if provided data need not be provided. If None, enc_data is set from data.

• enc_vdata (Optional[List[Tuple[int, E0]]]) – Optional sequence encoded validation set.

• out (I0) – Text output stream.

• print_iter (int) – Print iterations (i.e. log-likelihood difference) every print_iter-iterations.

• reuse_estep_ll (bool) – Default True. Forwarded to each trial’s optimize call – reuse the E-step likelihood for convergence instead of a separate scoring pass (see optimize). Set False to force the exact historical per-iteration scoring behavior.

• objective (str) – Convergence/selection objective forwarded to each trial’s optimize call; 'auto' (default) selects MLE / MAP / variational Bayes from the prior (see optimize).

• seed (Optional[int]) – Integer seed – shorthand for rng=RandomState(seed). Mutually exclusive with rng (passing both raises TypeError).

• fused_options (Optional[dict]) – Fused-kernel tuning knobs forwarded to every trial’s optimize call – see optimize(fused_options=...) for the recognized keys.

Returns:

Tuple of log-likelihood of best fitting model and the best fitting model from number of trials.

Return type:

tuple[float, SequenceEncodableProbabilityDistribution]

run_em(enc_data, estimator, initial_model, strategy=None, max_its=10, delta=1.0e-9, engine=None, objective=None, max_iter=None, monotone=True, tolerance=1.0e-12)

Run an EM-family strategy until convergence or max_its.

objective takes the same values as optimize(): None (MLE), a selection string ('auto' / 'mle' / 'map' / 'vb'), or a ready model -> float callable. max_its is the canonical iteration-cap spelling (matching optimize / fit / best_of); max_iter is accepted as a back-compat alias and overrides max_its when given. monotone=True transactionally rejects an objective decrease independently of whether delta is disabled. Set it to False only for a strategy whose changing surrogate objective intentionally permits temporary decreases in the supplied reporting objective.

The returned model carries a FitProvenance receipt (model.fit_provenance()): iterations actually run, whether the run converged on delta or stopped at the max_its cap (or on a rejected step), and the final objective. A run exiting on the cap reports converged=False, so a truncated fit is distinguishable from a finished one. A run whose very first step is rejected accepts nothing and returns initial_model unchanged; that case additionally raises a UserWarning, because the returned model is the initialization rather than a fit and the receipt alone does not say so.

Parameters:

• enc_data (Any)

• estimator (ParameterEstimator)

• initial_model (SequenceEncodableProbabilityDistribution)

• strategy (EMStrategy | None)

• max_its (int)

• delta (float | None)

• engine (Any | None)

• objective (str | Callable[[Any], float] | None)

• max_iter (int | None)

• monotone (bool)

• tolerance (float)

Return type:

SequenceEncodableProbabilityDistribution

class CompiledEM(*, compute_dtype=None)

Bases: object

Local full-mixture EM using the profiled fused component kernels.

This is the reusable form of the full-tree execution path used by block EM. It keeps the ordinary EM fixed point and update semantics, but can avoid generic accumulator traversal for fusible subtrees of a MixtureDistribution. Unsupported models or remote execution engines fall back to StandardEM.

Parameters:

compute_dtype (Any)

step(enc_data, estimator, model, engine=None, objective=None)

Run one compiled full-mixture E/M sweep.

Parameters:

• enc_data (Any)

• estimator (ParameterEstimator)

• model (SequenceEncodableProbabilityDistribution)

• engine (Any | None)

• objective (Callable[[Any], float] | None)

Return type:

EMStepResult

class SquaremEM(base_strategy=None, packer=None, tolerance=1.0e-12, step_growth=4.0)

Bases: object

SQUAREM acceleration (Varadhan & Roland 2008, SqS3) with an objective gate.

One step runs one SQUAREM cycle – THREE base-strategy sweeps: two ordinary EM sweeps give the secant pair r = theta1 - theta0, v = (theta2 - theta1) - r; the squared-extrapolation proposal theta0 - 2*alpha*r + alpha^2*v with alpha = min(-1, -sqrt(r.r/v.v)) is then stabilized by a third EM sweep. The gate keeps the stabilized proposal only if its objective is at least the plain two-sweep result’s; otherwise the cycle falls back to that plain result, so the accepted sequence is monotone whenever the base strategy is (StandardEM’s exact M-step is). An unpack that rejects the proposal (a constraint violation surfacing as ValueError / FloatingPointError / OverflowError) counts as a rejected cycle, never a crash.

Cost accounting is honest: max_its iterations of optimize(strategy=SquaremEM()) spend up to 3 * max_its E-step sweeps. The probe receipt for the win (overlapping 6-component GMM, n=100k): plain EM needed 200 sweeps to a target log-likelihood that SQUAREM reached in 33 sweeps, every cycle accepted, monotone throughout.

packer: (pack, unpack) between the model and an unconstrained parameter vector; default squarem_packer() (supported families documented there).

The strategy is STATEFUL across steps (the adaptive trust-region cap on alpha): use a fresh instance per fit, exactly as optimize(strategy=SquaremEM()) constructs one.

Parameters:

• base_strategy (Any | None)

• packer (tuple[Callable[[Any], np.ndarray], Callable[[np.ndarray], Any]] | None)

• tolerance (float)

• step_growth (float)

step(enc_data, estimator, model, engine=None, objective=None)

One SQUAREM cycle: two base sweeps, one gated extrapolation sweep.

Parameters:

• enc_data (Any)

• estimator (ParameterEstimator)

• model (SequenceEncodableProbabilityDistribution)

• engine (Any | None)

• objective (Callable[[Any], float] | None)

Return type:

EMStepResult

class SampledSufficientStatistics(suff_stat, nobs=None)

Bases: object

Unambiguous Monte-Carlo EM payload for an optional observation count.

Parameters:

• suff_stat (Any)

• nobs (float | None)

squarem_packer(model)

Build (pack, unpack) for SquaremEM over model’s primary parameters.

Supported out of the box: recursively nested MixtureDistribution and CompositeDistribution nodes with Gaussian / Laplace / Exponential / Poisson / Categorical leaves and no priors attached (a MAP fit changes the fixed point). Anything else raises NotImplementedError with the escape hatch named: pass an explicit packer=(pack, unpack) to SquaremEM for custom models.

pack(model) -> theta and unpack(theta) -> model round-trip losslessly (asserted in tests); mixture weights travel as log-weights re-normalized on unpack.

Parameters:

model (Any)

Return type:

tuple[Callable[[Any], ndarray], Callable[[ndarray], Any]]

class EMStrategy(*args, **kwargs)

Bases: Protocol

Structural contract for an EM-family strategy consumed by run_em().

Every strategy object in this module (StandardEM, PosteriorTransformEM, AnnealedEM, …) satisfies this Protocol structurally by exposing a step(...) -> EMStepResult method. run_em and _em_step_fn dispatch on it polymorphically; membership is decided by isinstance().

step(enc_data, estimator, model, engine=..., objective=...)

Run one EM-style update and return the resulting model and objective metadata.

Parameters:

• enc_data (Any)

• estimator (ParameterEstimator)

• model (SequenceEncodableProbabilityDistribution)

• engine (Any | None)

• objective (Callable[[Any], float] | None)

Return type:

EMStepResult

class StreamingEstimator(estimator, schedule=None, model=None, init_estimator=None, init_p=0.1, rng=None, encoder=None, num_chunks=1, dataset_size=None)

Bases: _StreamingBase

Decay-mode online estimator built from accumulator scaling and M-steps.

Parameters:

• estimator (ParameterEstimator)

• model (SequenceEncodableProbabilityDistribution | None)

• init_estimator (ParameterEstimator | None)

• init_p (float)

• rng (RandomState | None)

• num_chunks (int)

• dataset_size (float | None)

update(data=None, *, enc_data=None)

Consume one batch and return the updated model.

Every scale and count is validated finite before it touches the running accumulator, and the running state is committed only once the M-step has succeeded – so a failed update leaves the estimator exactly where it was and the caller can retry.

Parameters:

• data (Sequence[T] | None)

• enc_data (list[tuple[int, E0]] | None)

Return type:

SequenceEncodableProbabilityDistribution

class IncrementalEstimator(estimator, model=None, init_estimator=None, init_p=0.1, rng=None, encoder=None, num_chunks=1)

Bases: _StreamingBase

Neal-Hinton style incremental EM over replaceable data chunks.

Each chunk contributes a sufficient-statistic payload computed under the current model. Revisiting a chunk replaces that payload and rebuilds the pooled accumulator from the stored chunk summaries before running the ordinary M-step. Rebuilding is essential because mergeable statistics form a monoid, not necessarily a group: sums can be subtracted, but support minima/maxima, weighted medians, sketches, and similar summaries cannot.

Parameters:

• estimator (ParameterEstimator)

• model (SequenceEncodableProbabilityDistribution | None)

• init_estimator (ParameterEstimator | None)

• init_p (float)

• rng (RandomState | None)

• num_chunks (int)

update(data=None, *, enc_data=None, chunk_id=None)

Replace one chunk contribution and return the updated model.

chunk_id is keyword-only so this matches StreamingEstimator.update()’s (data, *, enc_data) shape across the streaming surface; it is required (a None chunk_id raises) because the Neal-Hinton update keys each batch’s contribution by it.

The chunk store, running statistics and model are committed only once the M-step has succeeded, so a failed update leaves the estimator exactly where it was – in particular the previous payload for a revisited chunk_id survives – and the caller can retry.

Parameters:

• data (Sequence[T] | None)

• enc_data (list[tuple[int, E0]] | None)

• chunk_id (Any)

Return type:

SequenceEncodableProbabilityDistribution

chunk_value(chunk_id)

Return a copy of one stored chunk contribution.

Parameters:

chunk_id (Any)

reset()

Drop all chunk contributions and fitted model state.

Return type:

None

class BayesianStreamingEstimator(estimator, mode='posterior_carry', schedule=None, model=None, init_estimator=None, init_p=0.1, rng=None, num_chunks=1)

Bases: object

Streaming / recursive-Bayes driver over the mixle.stats estimator protocol.

mode='posterior_carry' (the default) performs exact recursive conjugate updating: each fitted posterior is carried forward as the next batch’s prior by rebuilding the estimator from the fitted model (model.estimator() returns an estimator whose prior is the model’s posterior).

mode='forgetting' retains a running sufficient-statistic accumulator. Before each batch after the first, the retained OLD statistics and effective observation count are scaled by rho = schedule(step); the new batch is then added at full weight. Estimation always uses the original estimator/prior, avoiding both discounting new evidence and double-counting a carried posterior. Accumulator scaling preserves structural support metadata such as categorical bounds.

The public surface is BayesianStreamingEstimator(estimator, mode=..., schedule=...) plus .update(data=None, enc_data=None) and .reset().

Parameters:

• estimator (ParameterEstimator)

• mode (str | None)

• schedule (Any | None)

• model (SequenceEncodableProbabilityDistribution | None)

• init_estimator (ParameterEstimator | None)

• init_p (float)

• rng (RandomState | None)

• num_chunks (int)

update(data=None, enc_data=None)

Consume one batch and return the updated posterior-bearing model.

Parameters:

• data (Sequence[T] | None)

• enc_data (Any | None)

Return type:

SequenceEncodableProbabilityDistribution

reset()

Restore the original estimator, RNG state, and empty stream state.

Return type:

None

bayes_action(posterior, loss, actions, *, n=2000, seed=0, cvar_alpha=0.1, quantiles=(0.05, 0.5, 0.95), vectorized=None, report='selection')

Pick the Bayes action: argmin_a E_{draw ~ posterior}[ loss(a, draw) ].

The winner’s reported numbers carry selection bias unless you pay to remove it. The chosen action minimises K NOISY Monte-Carlo estimates, so its own expected_loss is optimistic post-argmin (the winner’s curse): measured bias about -2.5 Monte-Carlo standard errors at K = 100 candidates – and it does NOT shrink with n, because more draws shrink the SE and the bias together. expected_loss_mc_se is returned so the optimism has a scale; with many near-tied candidates, read the winner’s advertised loss as roughly that many SEs too good. report="fresh-draws" removes the bias instead: the choice is made exactly as before, then n NEW posterior draws are taken and the winner’s reported expected_loss / risk_profile are recomputed on them (selection-independent, hence conditionally unbiased for the chosen action). That costs one more samples(n, rng) call and n more loss evaluations (one vectorized call) FOR THE WINNER ONLY – deliberately not the default, because a loss may be metered or side-effecting and the exact-invocation-count contract below is itself a stated guarantee; selection_expected_loss then preserves the in-sample value the action won with. alternatives always carries selection-stage estimates.

Parameters:

• posterior (Any) – any object exposing samples(n, rng) – e.g. mixle.inference.posterior(model, data, over=...) (parameter, latent, or predictive).

• loss (Callable[[Any, Any], float]) – loss(action, draw) -> float (or a numpy-vectorized loss(action, draws) -> array).

• actions (Sequence[Any]) – the finite candidate-action set to minimise over.

• n (int) – number of posterior draws for the Monte-Carlo expectation.

• seed (int) – RNG seed for the posterior draw (reproducible).

• cvar_alpha (float) – tail mass for the CVaR / VaR of the chosen action (0.1 -> worst 10%).

• quantiles (Sequence[float]) – loss quantiles to report per action. Reported quantiles (and the CVaR/VaR fields) are plug-in Monte Carlo estimates and carry NO uncertainty quantification of their own – tail quantities converge far more slowly than the mean whose expected_loss_mc_se is reported, so raise n before leaning on a far tail (STAT-RR19-17).

• vectorized (bool | None) – the loss’s calling convention. True -> called once per action with the whole draw array; False -> called once per draw; None (default) -> read loss.vectorized if the loss declares it, else auto-detect by probing the array call once. A probe is a real invocation of the loss, so a loss that keeps state, records calls, or has side effects should declare its convention (argument or attribute) to be evaluated exactly the required number of times.

• report (str) – "selection" (default) reports the winner’s selection-stage estimates plus their Monte-Carlo SE; "fresh-draws" re-estimates the winner on n new draws (see above).

Returns:

{action, action_index, expected_loss, expected_loss_mc_se, risk_profile, alternatives, report, vectorized} – plus selection_expected_loss when report="fresh-draws".

Raises:

• ValueError – if actions is empty, n is not positive, cvar_alpha is outside (0, 1], a quantile level is outside [0, 1], the posterior returns an empty, ragged, or wrong-sized sample batch (samples(n, rng) must return exactly n draws), or a loss evaluates to NaN.

• TypeError – if n is not an exact integer, or cvar_alpha/a quantile level is a Boolean, a string, or anything else that is not already a real number.

• CapabilityError – if posterior does not expose the samples(n, rng) contract.

Return type:

dict[str, Any]

class RiskProfile(expected_loss, cvar, cvar_alpha, var, quantiles, std)

Bases: object

The tail-risk summary of a single action’s posterior loss distribution.

Parameters:

• expected_loss (float)

• cvar (float)

• cvar_alpha (float)

• var (float)

• quantiles (dict[float, float])

• std (float)

as_dict()

Return risk metrics and quantiles as JSON-compatible data.

Return type:

dict[str, Any]

certify(model, *, data=None, receipts=(), escape_tested=False, penalized=False)

Return the EstimationCertificate for a fitted model (or distribution prototype).

Model structure supplies only a candidate guarantee. data is checked by verify_estimation_conditions(); additional optimizer evidence may be supplied as structured receipts. A bare escape_tested=True is rejected because a caller assertion is not verification evidence.

Pass penalized (a reason string, or True) when the fit optimized a penalized objective – soft constraints, conservation/PINN residual factors, potentials (E2). The optimum is then of the penalized surrogate, not the likelihood, so no block may claim more than STATIONARY however clean its own solver is: every stronger block is downgraded with the penalty named in its reason.

Parameters:

• model (Any)

• data (Iterable[Any] | None)

• receipts (Iterable[VerificationReceipt])

• escape_tested (bool)

• penalized (str | bool)

Return type:

EstimationCertificate

plan_estimation(model)

Return a pre-fit plan whose guarantees remain explicitly unverified.

Parameters:

model (Any)

Return type:

EstimationCertificate

schedule(model, *, escape_tested=False, receipts=())

Plan the block-coordinate estimation schedule for model (planner v2, A3).

Built from the same block classification as certify(): EM blocks make the schedule a loop (E-step + per-block M-steps, repeated until convergence); without a latent block every block is one independent pass. Gradient blocks appear as explicit gradient passes with their pool placement, so the schedule is also the offload plan for the hybrid case.

receipts exists because escape_tested was otherwise unusable in both directions (MXR-080-1905). It forwards to certify(), which refuses escape_tested=True with “supply a VerificationReceipt” – and this function had no parameter through which to supply one, so the only non-default value raised an error naming a remedy the caller could not reach. Pass escape evidence here instead, bound to the model via receipt_subject() (MXR-080-1877).

Parameters:

• model (Any)

• escape_tested (bool)

• receipts (Iterable[VerificationReceipt])

Return type:

EstimationSchedule

class EstimationSchedule(passes=<factory>, latent=False)

Bases: object

The block-coordinate schedule planner v2 produces (A3): ordered passes + the loop structure.

A fully-factorized model schedules one independent pass per block (no loop). A latent model schedules the EM loop explicitly: an E-step over the latent, then one M-step pass PER BLOCK per round – each M-step named with its own method, so the schedule shows exactly where the closed forms live inside the iteration (and which pass, if any, is the gradient block a pool would take).

Parameters:

• passes (list[SchedulePass])

• latent (bool)

property per_round: list[SchedulePass]

Return schedule passes repeated inside each latent-variable round.

property gradient_passes: list[SchedulePass]

Return schedule passes assigned to gradient optimization.

describe()

Return a compact prose description of the estimation schedule.

Return type:

str

class SchedulePass(order, kind, block, method, placement, repeat)

Bases: object

One pass of the estimation schedule: what runs, on which block, how, where, and how often.

Parameters:

• order (int)

• kind (str)

• block (str)

• method (str)

• placement (str)

• repeat (str)

class EstimationCertificate(guarantee, blocks=<factory>, escape_tested=False)

Bases: object

The auditable proof of how a model was (or would be) estimated: per-block plans + the aggregate.

guarantee is the minimum over blocks. why_not_adam summarizes where gradient optimization was used and why those blocks could not use a stronger closed-form or convex route.

Parameters:

• guarantee (Guarantee)

• blocks (list[BlockPlan])

• escape_tested (bool)

property gradient_blocks: list[BlockPlan]

Return blocks that require gradient optimization.

property closed_form_blocks: list[BlockPlan]

Return blocks with global or stronger guarantees.

as_dict()

Return the certificate as JSON-compatible data.

Return type:

dict[str, Any]

why_not_adam()

The audit: which blocks needed gradient descent and why – everything else got a stronger method.

Return type:

str

table()

Render a human-readable block-level estimation table.

Return type:

str

class BlockPlan(name, kind, method, guarantee, gradient, placement, reason, candidate_guarantee=None, proof_obligations=(), verified_checks=(), receipt_ids=(), verified_by='', unbound_receipt_ids=())

Bases: object

One estimation block and the distinction between planned and verified strength.

Parameters:

• name (str)

• kind (str)

• method (str)

• guarantee (Guarantee)

• gradient (bool)

• placement (str)

• reason (str)

• candidate_guarantee (Guarantee | None)

• proof_obligations (tuple[ProofObligation, ...])

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

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

• verified_by (str)

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

class Guarantee(*values)

Bases: IntEnum

How strong the solution to an estimation block is, as an ordered ladder (higher = stronger).

property label: str

Return the enum name used in reports.

class ProofObligation(check, required_for, description)

Bases: object

A condition that must be evidenced before a block may claim a guarantee.

Parameters:

• check (str)

• required_for (Guarantee)

• description (str)

class VerificationReceipt(receipt_id, block, guarantee, checks, source, evidence, subject_hash='')

Bases: object

Evidence that named checks were performed for one estimation block.

Receipts are explicit audit inputs, not Boolean assertions. evidence must contain the measurements or result metadata supporting the checks.

Parameters:

• receipt_id (str)

• block (str)

• guarantee (Guarantee)

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

• source (str)

• evidence (dict[str, Any])

• subject_hash (str)

evidence_as_dict()

Return the evidence as plain JSON-compatible data.

Return type:

dict[str, Any]

receipt_subject(model)

Return the subject identifier a VerificationReceipt must carry to be about model.

A caller supplying optimizer evidence the library cannot produce itself – multi-start escape testing, a solver’s own convergence report – names the model it verified:

VerificationReceipt(..., subject_hash=receipt_subject(model))

Requiring this is the point (MXR-080-1877): before it, naming the right check strings was the whole test, so a receipt whose evidence read {"trust": "me"} upgraded an arbitrary GaussianDistribution(999, 1) to GLOBAL_UNIQUE. Binding does not make forgery impossible – anyone holding the model can compute this – but it means a receipt cannot be written without the model, cannot be reused against a different one, and stops applying the moment the parameters change. What a reader gains is in the certificate: BlockPlan.verified_by names who established each guarantee, so library-checked conditions no longer read identically to a caller’s assertion.

Parameters:

model (Any)

Return type:

str

verify_estimation_conditions(model, data)

Check data-dependent proof obligations for supported analytic estimators.

Unsupported or iterative estimators produce no receipt rather than a guessed guarantee. The returned receipts can be passed to certify().

Parameters:

• model (Any)

• data (Iterable[Any])

Return type:

list[VerificationReceipt]

uq(thing, data=None, *, alpha=0.1, equivalent=None)

Quantify the uncertainty of thing, choosing the method from what thing is.

Parameters:

• thing (Any) – a fitted mixle model, a torch module / point-predictor callable (or a list of them for a deep ensemble), or an LLM-style callable that maps a prompt to a generation.

• data (Any) – for a mixle model, the fitting data (builds a Laplace likelihood-curvature approximation); for a point predictor, (X_cal, y_cal) calibration data; for an LLM, optional example prompts used to calibrate an abstention threshold.

• alpha (float) – target miscoverage / abstention level (1 - alpha coverage).

• equivalent (Callable[[Any, Any], bool] | None) – for the LLM path, an optional meaning-equivalence predicate over generations (default: exact string match after stripping).

Returns:

A UQResult exposing the method-appropriate accessors and its own calibration numbers.

Return type:

UQResult

class UQResult(kind, method, payload)

Bases: object

The uncertainty of a predictor, with the method that produced it and receipts to check it.

Parameters:

• kind (str)

• method (str)

• payload (dict[str, Any])

sample_models(n=200, *, seed=None)

n fitted models drawn from the parameter approximation.

Parameters:

• n (int)

• seed (int | None)

Return type:

list[Any]

credible_interval(readout, alpha=0.1, *, n=400, seed=0)

A 1-alpha credible interval on readout(model) over the parameter posterior.

Parameters:

• readout (Callable[[Any], float])

• alpha (float)

• n (int)

• seed (int)

Return type:

tuple[float, float]

parameter_interval(readout, alpha=0.1, *, n=400, seed=0)

A central interval over draws from the disclosed parameter approximation.

Parameters:

• readout (Callable[[Any], float])

• alpha (float)

• n (int)

• seed (int)

Return type:

tuple[float, float]

interval(x, alpha=None)

Calibrated prediction interval(s) at x. alpha overrides the calibrated level.

Parameters:

• x (Any)

• alpha (float | None)

Return type:

tuple[ndarray, ndarray]

epistemic_std(x)

Ensemble disagreement (std across members) at x – 0.0 for a single predictor.

Parameters:

x (Any)

Return type:

ndarray

semantic_entropy(prompt, *, n=8)

Entropy (nats) over the meaning classes of n sampled generations for prompt.

Parameters:

• prompt (Any)

• n (int)

Return type:

float

confident(prompt, *, n=8, max_entropy=None)

True when BIAS-CORRECTED semantic entropy is below the threshold.

The plug-in entropy at small n is biased LOW by ~(K-1)/(2n) nats – measured -0.656 at the 8-draw default on a uniform ten-class generator – which makes a hallucinating model look confident, the exact failure this gate exists to catch (STAT-RR17-14). The decision therefore uses the Miller-Madow-corrected estimate: still an estimate, but its residual error is far smaller and the correction errs toward ABSTAINING. The calibrated threshold is a quantile of the SAME corrected statistic at the SAME draw count (STAT-RR19-12: a plug-in-scale threshold against a corrected-scale gate moved acceptance from 98.2% to 41.0% on a fixed generator), so serving a calibrated threshold at a different n is refused – the correction is (K-1)/(2n), and changing n changes the statistic the threshold was calibrated for. An explicit max_entropy override carries no calibration claim and may use any n. Use semantic_entropy_receipt() to see the plug-in value, the correction, n, and K.

Parameters:

• prompt (Any)

• n (int)

• max_entropy (float | None)

Return type:

bool

semantic_entropy_receipt(prompt, *, n=8)

The entropy estimate WITH its approximation receipt (n, K, bias, corrected value).

Parameters:

• prompt (Any)

• n (int)

Return type:

dict[str, Any]

report()

Return uncertainty-quantification metadata and scalar payload fields.

Return type:

dict[str, Any]

calibration_report(model, data, *, seed=0)

The calibration of model on held-out data (see module docstring).

data should be data the model was not fitted on – calibration measured on the training set is optimistic. Runs the PIT test when the model has a scalar predictive CDF; always reports the held-out mean log-density.

A predictive law with ATOMS is handled with the randomized PIT rather than the continuous one: F(Y) is uniform only for a continuous law, so routing a discrete model through the ordinary PIT reports a large error for an exactly correct model. seed fixes that randomization (it is recorded in note) and is unused for a continuous law, whose method is unchanged.

Scoring validity is checked first and carried in the report. The verdict used to consult only the PIT histogram, so a model whose CDF was exact but whose per-row log-densities were all NaN – or which returned one score for a thousand rows – passed the gate as "calibrated". Exactly one usable score per row is now required for a calibrated verdict; a cardinality mismatch raises and NaN scores are counted into n_invalid_scores, which forces is_calibrated() to False.

Raises:

ValueError – if the model’s batch scoring does not return exactly one log-density per row.

Parameters:

• model (Any)

• data (Any)

• seed (int)

Return type:

CalibrationReport

class CalibrationReport(n, mean_log_density, pit_error=None, pit_histogram=None, bins=10, method='', note='', n_invalid_scores=0)

Bases: object

Whether a fitted model’s uncertainty is calibrated on held-out data.

pit_error is the total-variation distance of the PIT histogram from uniform (0 = perfectly calibrated). It has a finite-sample floor ~``sqrt(bins/n)`` even for a perfect model, so is_calibrated() judges against that floor rather than a fixed constant.

Parameters:

• n (int)

• mean_log_density (float)

• pit_error (float | None)

• pit_histogram (dict[str, Any] | None)

• bins (int)

• method (str)

• note (str)

• n_invalid_scores (int)

noise_floor()

The PIT-error a perfectly calibrated model would show at this sample size (sampling noise).

Return type:

float

scoring_is_valid()

Whether every held-out row got a usable (finite or explicitly impossible) log-density.

Return type:

bool

max_pit_error()

The largest PIT error this histogram statistic can attain: all mass in one bin.

sum_b |freq_b - 1/bins| peaks at (1 - 1/bins) + (bins-1)/bins == 2(1 - 1/bins) – 1.8 for the default ten bins. A tolerance at or above this bound accepts every possible histogram, so the gate is not a test at all.

Return type:

float

has_enough_data()

Whether the holdout can support a per-bin uniformity judgement at all.

The histogram statistic needs a usable expected count per bin; below MIN_EXPECTED_PER_BIN the noise floor sqrt(bins/n) grows until 2.5x it exceeds max_pit_error() and the acceptance region becomes vacuous – at ten bins and a single observation the threshold was 7.91 against a statistic that cannot exceed 1.8, so even maximally clumped PIT values were reported calibrated.

Return type:

bool

is_calibrated(tol=None)

True when the PIT error is within tolerance AND the model actually scored the holdout.

Default tol = 2.5x the finite-sample noise floor (so genuine miscalibration, not sampling noise, is what fails). Unknown -> False, conservatively. An invalid predictive score on any row -> False: the PIT histogram alone cannot certify a model whose density is unavailable. Too little held-out data for has_enough_data() -> False: absence of evidence is not evidence of calibration.

Raises:

ValueError – if tol is not finite and nonnegative, or is at/above max_pit_error() (a tolerance that no histogram can exceed would approve any model).

Parameters:

tol (float | None)

Return type:

bool

as_dict()

Return rounded calibration metrics as JSON-compatible data.

Return type:

dict[str, Any]

to_dataframe()

Return this report as a pandas.DataFrame.

When a PIT histogram was computed (method == "PIT"), returns one row per bin – bin_left, bin_right, count, density, uniform – the natural columnar shape of a calibration report (a table you’d plot). Otherwise (no scalar predictive CDF, so no PIT histogram) returns a single summary row with the same unrounded fields as as_dict(). Requires the pandas extra (pip install mixle[pandas]).

Return type:

Any

to_parquet(path, **kwargs)

Write this report to a Parquet file; see to_dataframe().

kwargs forward to DataFrame.to_parquet (e.g. engine=, compression=). Needs a Parquet engine in addition to pandas – pip install mixle[arrow] (pyarrow) or fastparquet.

Parameters:

• path (Any)

• kwargs (Any)

Return type:

None

plan_placement(certificate, pool=None, *, telemetry=None, model_digest=None, data_digest=None, version_digest=None)

Decide local vs pool for every block of a certified fit (see module docstring).

Closed-form / convex / EM blocks always stay local. A gradient (pool_eligible) block goes to the pool only when pool.available and its estimated work clears pool.flop_threshold_tflop; otherwise it stays local too. Each decision carries a priced reason, and (with a pool) a placement telemetry event is emitted per block.

Parameters:

• certificate (EstimationCertificate)

• pool (PoolSpec | None)

• telemetry (Any)

• model_digest (str | None)

• data_digest (str | None)

• version_digest (str | None)

Return type:

PlacementPlan

class PlacementPlan(placements=<factory>, pool_spec=None, context_digest=None, bindings=<factory>, execution_status='not_executed')

Bases: object

An unexecuted local-vs-pool plan, never an execution receipt.

Parameters:

• placements (list[BlockPlacement])

• pool_spec (PoolSpec | None)

• context_digest (str | None)

• bindings (dict[str, str | None])

• execution_status (str)

property pool_blocks: list[BlockPlacement]

Return blocks assigned to pool execution.

property local_blocks: list[BlockPlacement]

Return blocks assigned to local execution.

property est_pool_cost: float

Return estimated total pool spend.

as_dict()

Return the placement plan as JSON-compatible data.

Return type:

dict[str, Any]

report()

Render a human-readable placement report.

Return type:

str

class BlockPlacement(name, kind, placement, reason, est_tflop=0.0, est_cost=0.0, estimate_source='none', execution_status='not_executed', observed_placement=None)

Bases: object

Where one block runs and why: ‘local’ or ‘pool’, with the priced justification.

Parameters:

• name (str)

• kind (str)

• placement (str)

• reason (str)

• est_tflop (float)

• est_cost (float)

• estimate_source (str)

• execution_status (str)

• observed_placement (str | None)

class PoolSpec(available=False, cost_per_hour=0.5, flop_threshold_tflop=1.0, pool_speedup=10.0)

Bases: object

What the shared GPU pool offers, and when offloading is worth it.

available: is a pool configured at all (else everything is local). cost_per_hour: the pool’s price. flop_threshold: the estimated work (in TFLOP) below which the round-trip is not worth it – small gradient blocks stay local. local_tflops / pool speedup drive the estimate.

Parameters:

• available (bool)

• cost_per_hour (float)

• flop_threshold_tflop (float)

• pool_speedup (float)

learn_placement_policy(rows, static_policy, *, cost_key='cost', k=8, min_neighbors=4)

Learn a placement policy from telemetry (features, choice, outcome) rows (see module docstring).

static_policy maps a feature dict to a choice and is the fall-back when history is too thin. cost_key names the outcome field to minimize (default "cost"). Feature standardization and the neighbor index are built from the rows; LearnedPolicy.decide() and evaluate follow.

Parameters:

• rows (list[tuple[dict[str, Any], str, dict[str, Any]]])

• static_policy (Callable[[dict[str, Any]], str])

• cost_key (str)

• k (int)

• min_neighbors (int)

Return type:

LearnedPolicy

class LearnedPolicy(keys, vecs, choices, costs, static, mean=<factory>, scale=<factory>, k=8, min_neighbors=4, margin=0.02)

Bases: object

A history-based placement policy that defers to a static teacher where it lacks evidence.

Parameters:

• keys (list[str])

• vecs (ndarray)

• choices (list[str])

• costs (ndarray)

• static (Callable[[dict[str, Any]], str])

• mean (ndarray)

• scale (ndarray)

• k (int)

• min_neighbors (int)

• margin (float)

decide(features)

Return (choice, learned) – the learned pick when confident, else the static fallback.

Parameters:

features (dict[str, Any])

Return type:

tuple[str, bool]

evaluate(rows, *, cost_key='cost')

Compare policies on the same contexts with actually observed action costs.

A context contributes only when both the learned and static selected actions have realized outcomes in the supplied holdout panel. No training-neighbor estimate is substituted for a missing outcome.

Parameters:

• rows (list[tuple[dict[str, Any], str, dict[str, Any]]])

• cost_key (str)

Return type:

dict[str, Any]

learn_action_policy(rows, static_scorer=None, *, value_key='value', k=8, min_neighbors=4)

Learn a reasoner acquisition policy from route telemetry (features, kind, outcome) rows.

static_scorer is the fall-back when history is thin (default mixle.substrate.act.score_action()). value_key names the outcome field to MAXIMIZE (default "value" – did the action yield evidence). Returns a LearnedAcquisition usable directly as investigate(..., scorer=policy).

Parameters:

• rows (list[tuple[dict[str, Any], str, dict[str, Any]]])

• static_scorer (Callable[[Any, str], float] | None)

• value_key (str)

• k (int)

• min_neighbors (int)

Return type:

LearnedAcquisition

class LearnedAcquisition(keys, vecs, values, static, mean=<factory>, scale=<factory>, k=8, min_neighbors=4)

Bases: object

A history-based action scorer for the reasoner: learns which actions pay off, else defers.

Drop-in for mixle.substrate.act.score_action() (call it as scorer=policy in investigate). From route telemetry – each row a fired action’s (features={kind,cost,overlap}, value) – it estimates the expected yield of an action in a query’s feature region and scores it yield / cost. Where nearby history is too thin, it falls back to the static lexical scorer. This is an abstention rule, not a guarantee that learned decisions are never worse.

Parameters:

• keys (list[str])

• vecs (ndarray)

• values (ndarray)

• static (Callable[[Any, str], float])

• mean (ndarray)

• scale (ndarray)

• k (int)

• min_neighbors (int)

expected_yield(features)

Estimated yield of an action with these features, or None when history is too thin to say.

Parameters:

features (dict[str, Any])

Return type:

float | None

learn_schedule_policy(rows, static_policy, *, latency_key='latency', k=8, min_neighbors=4)

Learned pool scheduling (J4): when work is pool-eligible, learn where and when it actually runs fastest.

The same evidence-gated shape as placement, keyed on realized latency instead of dollar cost: rows are (features, choice, outcome) where features describe the moment (queue depth, job size, local load), choice is the scheduling decision (“run_local” / “queue_pool” / “defer”), and the outcome’s latency is what the decision actually cost in wall-clock. Where nearby history is thin, the returned policy defers to the static scheduler.

Parameters:

• rows (list[tuple[dict[str, Any], str, dict[str, Any]]])

• static_policy (Callable[[dict[str, Any]], str])

• latency_key (str)

• k (int)

• min_neighbors (int)

Return type:

LearnedPolicy

meta_improve(rows, static_policy, *, holdout_rows=None, cost_key='cost', holdout_frac=0.3, seed=0, k=8, min_neighbors=4, propensity_key='propensities', min_propensity=0.05, min_overlap=0.8, confidence_level=0.95, min_improvement=0.0)

Learn from telemetry and promote only with isolated, supported holdout evidence.

Repeated outcomes for the same feature context are kept in one split. Evaluation first uses paired realized outcomes when both policies’ selected actions were observed for the same context. If a complete propensity map is logged instead, inverse-propensity evaluation is available. Promotion requires minimum overlap and the one-sided confidence bound on learned - static to beat the requested improvement threshold. Returns:

{promoted, policy, receipt: {learned_mean_cost, static_mean_cost,
 mean_cost_difference, upper_confidence_bound, overlap_fraction, n}}

policy is learned only when promoted; otherwise it is exactly static_policy. The receipt records either decision.

Parameters:

• rows (list[tuple[dict[str, Any], str, dict[str, Any]]])

• static_policy (Callable[[dict[str, Any]], str])

• holdout_rows (list[tuple[dict[str, Any], str, dict[str, Any]]] | None)

• cost_key (str)

• holdout_frac (float)

• seed (int)

• k (int)

• min_neighbors (int)

• propensity_key (str)

• min_propensity (float)

• min_overlap (float)

• confidence_level (float)

• min_improvement (float)

Return type:

dict[str, Any]

simulate(model)

Package a fitted model as a Simulator (see module docstring).

Parameters:

model (Any)

Return type:

Simulator

class Simulator(model)

Bases: object

A fitted model packaged as a data generator, runnable under a baseline or named scenarios.

Parameters:

model (Any)

scenario(name, interventions)

Register a named intervention scenario (requires a learned Bayesian network to apply).

Parameters:

• name (str)

• interventions (dict[int, Any])

Return type:

Simulator

run(n=100, *, scenario=None, interventions=None, seed=0)

Generate exactly n synthetic records under the baseline, a registered scenario, or ad-hoc interventions.

n must be an exact positive integer. Raises IncompleteSimulationError if the underlying generator yields a different number of records than requested.

Parameters:

• n (int)

• scenario (str | None)

• interventions (dict[int, Any] | None)

• seed (int)

Return type:

list[Any]

outcome_mean(field_index, *, scenario=None, n=2000, seed=0)

The mean of a numeric field over exactly n draws – the quantity to compare across conditions.

The run is size-checked first (see run()), so the returned mean is always over n observations rather than over however many the generator happened to yield.

Parameters:

• field_index (int)

• scenario (str | None)

• n (int)

• seed (int)

Return type:

float

compare(field_index, scenario_a, scenario_b, *, n=4000, seed=0)

mean(field | scenario_a) - mean(field | scenario_b) – the simulated effect of A vs B.

Parameters:

• field_index (int)

• scenario_a (str | None)

• scenario_b (str | None)

• n (int)

• seed (int)

Return type:

float

class Scenario(name, interventions=<factory>)

Bases: object

A named simulation condition: which fields are clamped to which values (an intervention).

Parameters:

• name (str)

• interventions (dict[int, Any])

synthesize(source, *, label=None, verify=None, n=100, max_tries=None, seed=0)

Build a verified dataset of n accepted rows from a generative source (see module docstring).

source is a fitted model (sampled), a list of real inputs (a generator is inferred over them), or a callable draw function. label (optional) is the teacher applied to each input. verify (optional) accepts verify(x) or verify(x, label) and gates each row – rejected rows are resampled up to max_tries total draws. The verifier is attached to the returned Dataset so consumers can recheck() independently.

n must be an exact non-negative integer (bool is not an integer here). A returned dataset always holds exactly n verified rows; if the attempt budget, an impossible verifier, a short sampler, or real-input deduplication leaves it short, IncompleteSynthesisError is raised with the partial dataset attached rather than a shortfall being returned as an ordinary result.

Parameters:

• source (Any)

• label (Callable[[Any], Any] | None)

• verify (Callable[[...], bool] | None)

• n (int)

• max_tries (int | None)

• seed (int)

Return type:

Dataset

class Dataset(inputs, labels=None, verify=None, acceptance_rate=1.0, n_rejected=0, provenance=<factory>)

Bases: object

A verified synthetic dataset: inputs, optional labels, and the verifier that vouched for them.

inputs and labels stay ordinary mutable lists, but a labeled dataset must always hold exactly one label per input. That alignment is checked at construction and again before every operation that walks the rows, so a post-construction edit that leaves the two out of step is reported instead of silently truncating the rows that get iterated, paired, or rechecked.

Parameters:

• inputs (list[Any])

• labels (list[Any] | None)

• verify (Callable[[...], bool] | None)

• acceptance_rate (float)

• n_rejected (int)

• provenance (dict[str, Any])

require_aligned()

Raise unless every input carries exactly one label (vacuous for an unlabeled dataset).

Return type:

None

pairs()

(input, label) pairs – raises if the dataset is unlabeled or misaligned.

Return type:

list[tuple[Any, Any]]

recheck()

Re-run the attached verifier over every row.

Returns True when every row still passes, or when no verifier is attached. Raises if the dataset is misaligned – a partial recheck that skipped unlabeled inputs would report a pass it never made.

Return type:

bool

class CyclicGroup(order, period=1.0)

Bases: object

Z_``order`` acting on a periodic real-valued coordinate of period period by rotation.

Each group element k in {0, ..., order - 1} is realized concretely as a rotation of the coordinate’s (cos, sin) embedding by angle 2*pi*k/order – an orthogonal (norm- and Jacobian-preserving) transform, so composing group elements is exactly addition mod order (compose()), and a density fit on the embedding is unaffected by which element aligned a point into it (the twist is undone before scoring, not baked into the density).

order must be an exact positive integer and period a finite positive number: a zero order divides by zero in the rotation law, a negative or fractional one describes no cyclic group at all, and a zero/non-finite period makes every embedding NaN. {0, ..., order - 1} is the whole element set, and it is checked – an out-of-range k used to be silently wrapped into a different element, so -1, order, order + 1 and True all acted as undeclared aliases.

Parameters:

• order (int)

• period (float)

element(k)

Validate k as a declared group element and return it as a plain int.

Parameters:

k (Any)

Return type:

int

embed(x)

The periodic coordinate’s (cos, sin) embedding, shape (..., 2).

Parameters:

x (Sequence[float])

Return type:

ndarray

act(embedded, k)

Rotate an (..., 2) embedding by group element k (the forward twist).

Parameters:

• embedded (ndarray)

• k (int)

Return type:

ndarray

inverse_act(embedded, k)

Undo group element k’s twist – act(inverse_act(v, k), k) == v.

Parameters:

• embedded (ndarray)

• k (int)

Return type:

ndarray

compose(k1, k2)

The group element equivalent to applying k1 then k2 – addition mod order.

Parameters:

• k1 (int)

• k2 (int)

Return type:

int

class TwistedMixtureResult(base_density, group)

Bases: object

A single shared base density plus the group whose elements twist it into each group’s local factor.

Parameters:

• base_density (Any)

• group (CyclicGroup)

log_normalizer()

The (cached, k-independent) log constant that normalizes the shared base density over a period.

Return type:

float

log_density(x, k)

log p(x | group=k): undo k’s twist, score under the shared base density, normalize.

A normalized density on the periodic coordinate – exp of this integrates to 1 over one period – so values are comparable against independent_log_density() and across fits.

Parameters:

• x (Sequence[float])

• k (int)

Return type:

ndarray

fit_twisted_mixture(group, data_by_group, *, n_components=2, seed=0, max_its=50)

Fit ONE base density on every group’s data pooled together after undoing each group’s twist.

data_by_group maps a group element k to that group’s (small) sample of the periodic coordinate. Every sample is embedded and rotated back by its own group’s inverse_act before pooling – so the fitted n_components-component density sees sum(len(v) for v in data_by_group.values()) points, not just one group’s slice, for the same parameter count as fit_independent_mixtures() spends on a SINGLE group.

Parameters:

• group (CyclicGroup)

• data_by_group (dict[int, Sequence[float]])

• n_components (int)

• seed (int)

• max_its (int)

Return type:

TwistedMixtureResult

fit_independent_mixtures(group, data_by_group, *, n_components=2, seed=0, max_its=50)

The untwisted baseline: one independently-fit n_components-component density per group.

Same per-group parameter count as the shared base density in fit_twisted_mixture(), but |groups| times the total parameters overall, and each fit sees only its own group’s (small) sample – the comparison fit_twisted_mixture() is measured against.

Parameters:

• group (CyclicGroup)

• data_by_group (dict[int, Sequence[float]])

• n_components (int)

• seed (int)

• max_its (int)

Return type:

dict[int, Any]

independent_log_density(models, group, x, k)

Score x under group k’s independently-fit density (the baseline sibling of TwistedMixtureResult.log_density()).

Normalized over one period by that model’s own log_circular_normalizer() – each independent fit has a different embedding normalizer, so without this the baseline comparison would be decided by normalization rather than by fit quality.

Parameters:

• models (dict[int, Any])

• group (CyclicGroup)

• x (Sequence[float])

• k (int)

Return type:

ndarray

create(data, *, budget=None, device=None, calibrate=None, quantify_uq=False, max_its=25, seed=0)

Create a certified model from data (see module docstring).

data is a list of records/scalars. calibrate (a fraction in (0,1)) reserves a holdout for a PIT calibration check. quantify_uq=True attaches an auto-selected UQ handle. budget / device (any object; recorded verbatim) constrain the fit toward a smaller model. Returns a CreatedModel bundling the fit, its certificate, and the requested post-conditions.

Parameters:

• data (Any)

• budget (Any | None)

• device (Any | None)

• calibrate (float | None)

• quantify_uq (bool)

• max_its (int)

• seed (int)

Return type:

CreatedModel

class CreatedModel(model, certificate, strategy, calibration=None, uq=None, provenance=<factory>, postconditions=<factory>)

Bases: object

A certified model artifact: the fitted model plus its guarantees and provenance.

postconditions records every post-fit condition the caller asked for and what became of it: {"calibration": {"requested": True, "performed": False, "error": "..."}, ...}. A requested postcondition that raised used to be replaced by a bare None field, so an artifact whose calibration and UQ had both failed was indistinguishable from one where neither was asked for – while still advertising the estimation certificate’s guarantee.

Parameters:

• model (Any)

• certificate (Any)

• strategy (str)

• calibration (Any | None)

• uq (Any | None)

• provenance (dict[str, Any])

• postconditions (dict[str, dict[str, Any]])

failed_postconditions()

Names of the postconditions the caller requested that did not complete.

Return type:

list[str]

is_certified()

Whether every requested postcondition actually completed.

The estimation certificate alone does not make an artifact certified: a caller who asked for calibration or UQ asked for it as part of the claim, so a silent failure invalidates the aggregate claim even though the fit itself succeeded.

Return type:

bool

property guarantee: Any

The aggregate estimation guarantee (MIN over blocks) – the artifact’s summary claim.

Reports Guarantee.UNVERIFIED when a requested postcondition failed: the artifact cannot expose a guarantee stronger than the weakest condition it was asked to establish.

why()

One line summarizing how the artifact was estimated.

Return type:

str

is_calibrated()

Whether the calibration holdout judged the model calibrated.

None means calibration was never requested. A requested calibration that failed is False, not None – the check was asked for and did not pass, which is a different statement from never having asked.

Return type:

bool | None

skill(name, obj, *, description='', tags=(), call=None, provenance=None, registry=None)

Package obj (a fitted model, a CreatedModel, or a function) as a reusable Skill and register it (see module docstring).

The estimation certificate is inherited from obj.certificate when present (so a certified model yields a certified skill). call overrides how the skill is invoked; otherwise a model verb (predict / __call__ / sampler) is used. The skill is added to registry (or the process default) and returned.

Parameters:

• name (str)

• obj (Any)

• description (str)

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

• call (Callable[[...], Any] | None)

• provenance (dict[str, Any] | None)

• registry (SkillRegistry | None)

Return type:

Skill

class Skill(name, call, description='', tags=(), certificate=None, provenance=<factory>, implementation_digest='', artifact_digest='', descriptor_digest='', _artifact=None)

Bases: object

A named, described, provenanced callable derived from a fitted artifact (or a plain function).

Parameters:

• name (str)

• call (Callable[[...], Any])

• description (str)

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

• certificate (Any | None)

• provenance (dict[str, Any])

• implementation_digest (str)

• artifact_digest (str)

• descriptor_digest (str)

• _artifact (Any)

property guarantee: Any | None

The inherited estimation guarantee, if this skill wraps a certified model.

score(query)

Lexical overlap of query with this skill’s name/description/tags, in [0, 1].

Parameters:

query (str)

Return type:

float

verify_identity()

Whether executable code, bound state, and descriptive policy still match this skill.

Return type:

bool

class SkillRegistry

Bases: object

A findable collection of skills – register once, retrieve by name or by query.

add(sk)

Register a skill and return it for chaining.

Parameters:

sk (Skill)

Return type:

Skill

get(name)

Return a registered skill by exact name.

Parameters:

name (str)

Return type:

Skill

all()

Return all registered skills in insertion order.

Return type:

list[Skill]

find(query, k=5)

The best k skills for query by lexical overlap (highest score first, ties by name).

Parameters:

• query (str)

• k (int)

Return type:

list[Skill]

best(query)

The single best-matching skill for query (or None if nothing overlaps).

Parameters:

query (str)

Return type:

Skill | None

index(substrate)

Mirror every skill into a Substrate for embedding-grade recall.

The registry’s own find is a lexical matcher; when a corpus grows past that, index the skills as artifact items and retrieve through the substrate instead. Returns the item ids.

Parameters:

substrate (Any)

Return type:

list[str]

default_registry()

The process-wide default registry skill() writes to when no registry is given.

Return type:

SkillRegistry

posterior(model, data=None, *, over='predictive', prior=None, method='auto', **kwargs)

Build the Posterior of model over the requested variables.

Parameters:

• model (Any) – a fitted mixle distribution (or a latent-variable model for over='latent').

• data (Any) – observations – required for over='params' and for over='latent' (the x the latent posterior conditions on); ignored for plug-in over='predictive'.

• over (str) – 'latent' -> q(z | x) (needs the latent_posterior capability); 'params' -> q(theta | data); 'predictive' -> draws of new data.

• prior (Any) – prior over parameters for over='params' (see conjugate_posterior / sample_parameter_posterior).

• method (str) – for over='params' – 'auto' (conjugate when the family supports it, else MCMC), 'conjugate', or 'mcmc'.

• **kwargs (Any) – forwarded to sample_parameter_posterior for the MCMC path (sampler, steps…).

Returns:

A Posterior – a LatentPosterior, ParameterPosterior, or PredictivePosterior.

Return type:

Posterior

class ParameterPosterior(draw_one, draw_many, *, mean_fn=None, chain=None, kind='')

Bases: Posterior

q(theta | data) over a model family’s parameters – exact (conjugate) or MCMC.

Built by posterior(model, data, over="params"); the exact-vs-MCMC distinction is hidden behind the shared Posterior interface. A single sample() returns one parameter set, and samples() returns n of them; mean() and interval() summarize the posterior.

Parameters:

• draw_one (Callable[[Any], Any])

• draw_many (Callable[[int, Any], Any])

• mean_fn (Callable[[], Any] | None)

• chain (np.ndarray | None)

• kind (str)

classmethod from_conjugate(cp)

Wrap a closed-form ConjugatePosterior (each draw is a parameter dict).

Parameters:

cp (ConjugatePosterior)

Return type:

ParameterPosterior

classmethod from_mcmc(result)

Wrap an MCMC MCMCResult; draws resample the retained chain, summaries use it directly.

Parameters:

result (Any)

Return type:

ParameterPosterior

sample(rng=None)

One parameter draw from the posterior.

Parameters:

rng (Any)

Return type:

Any

samples(n, rng=None)

n parameter draws (a dict of length-n arrays for conjugate; a list for MCMC).

Parameters:

• n (int)

• rng (Any)

Return type:

Any

mean()

The posterior mean of the parameters.

Return type:

Any

interval(level=0.9)

Central credible interval at level – [lo, hi] over the chain (MCMC) or 2000 draws.

Parameters:

level (float)

Return type:

Any

to_dataframe(n, rng=None)

Return n posterior parameter draws as a pandas.DataFrame, one row per draw.

Draws through samples(), so an MCMC-backed posterior resamples (with replacement) from its retained chain and a conjugate posterior draws exact iid samples. Columns are the parameter names when the sampler kept them (conjugate; MCMC over CategoricalDistribution); otherwise generic param_0..param_{d-1} (a single param_0 for a scalar family) since the other bridged MCMC families map to unnamed positional theta – see mixle.inference.mcmc.parameter_bridge. Requires the pandas extra (pip install mixle[pandas]).

Parameters:

• n (int)

• rng (Any)

Return type:

Any

to_parquet(path, n, rng=None, **kwargs)

Write n posterior parameter draws to a Parquet file; see to_dataframe().

kwargs forward to DataFrame.to_parquet (e.g. engine=, compression=). Needs a Parquet engine in addition to pandas – pip install mixle[arrow] (pyarrow) or fastparquet.

Parameters:

• path (Any)

• n (int)

• rng (Any)

• kwargs (Any)

Return type:

None

class PredictivePosterior(draw_one, draw_many)

Bases: Posterior

The posterior-predictive: draws of new data from a fitted model.

plug_in() wraps a single fitted model’s sampler. from_parameter_posterior() integrates parameter uncertainty – each predictive draw first samples theta ~ q(theta | data), rebuilds the model, then draws data – so the spread reflects both sampling and parameter uncertainty.

Parameters:

• draw_one (Callable[[Any], Any])

• draw_many (Callable[[int, Any], Any])

classmethod plug_in(model)

Plug-in predictive: draw new data from model at its fitted parameters.

Parameters:

model (Any)

Return type:

PredictivePosterior

classmethod from_parameter_posterior(param_post, build)

Posterior-predictive integrating parameter uncertainty.

build maps one parameter draw (the object ParameterPosterior.sample() returns) to a distribution; each predictive draw rebuilds the model from a fresh theta and samples it.

Parameters:

• param_post (ParameterPosterior)

• build (Callable[[Any], Any])

Return type:

PredictivePosterior

sample(rng=None)

One predictive draw of new data.

Parameters:

rng (Any)

Return type:

Any

samples(n, rng=None)

n predictive draws of new data.

Parameters:

• n (int)

• rng (Any)

Return type:

Any

class PriceForecast(mean, lo, hi, paths, level, calibration_count=0, interval_method='horizon_matched_split_conformal', coverage_assumptions=('held_out_exchangeability', 'stable_data_generating_process'), model_fit_length=0)

Bases: NamedTuple

A price/cost forecast with a horizon-matched interval calibration receipt.

Parameters:

• mean (ndarray)

• lo (ndarray)

• hi (ndarray)

• paths (ndarray)

• level (float)

• calibration_count (int)

• interval_method (str)

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

• model_fit_length (int)

mean: ndarray

Alias for field number 0

lo: ndarray

Alias for field number 1

hi: ndarray

Alias for field number 2

paths: ndarray

Alias for field number 3

level: float

Alias for field number 4

calibration_count: int

Alias for field number 5

interval_method: str

Alias for field number 6

coverage_assumptions: tuple[str, ...]

Alias for field number 7

model_fit_length: int

Alias for field number 8

forecast_price(model, history, horizon, *, level=0.9, cal_frac=0.3, seed=0, model_fit_length=None)

Forecast horizon steps of a price/cost series with a conformally-calibrated band.

Parameters:

• model (Any) – a fitted HiddenMarkovModelDistribution over the price/cost series (any scalar emission family with a sampler – see mixle.inference.forecast.forecast()).

• history (Any) – the observed series to condition on (one sequence, oldest first).

• horizon (int) – number of future steps to forecast.

• level (float) – central-interval mass (0.9 -> the calibrated 5%..95% band).

• cal_frac (float) – the fraction of history (most recent) reserved for calibration. Within that reserved window, a rolling-origin backtest of the same horizon-step-ahead forecast is scored against the real outcomes to build the calibration residuals.

• seed (int) – reproducibility for the calibration and requested-horizon Monte Carlo draws.

• model_fit_length (int | None) – REQUIRED DECLARATION (STAT-RR17-06) – how many leading observations of history the fitted model consumed (0 when it was built from prior knowledge or other data; len(history) when it was fitted on everything). The rolling-origin calibration reuses the model at every pseudo-origin without refitting, so when the model was fitted through ANY calibration-window outcome the residuals are partly in-sample and the band is optimistic: the receipt’s interval_method and coverage_assumptions then disclose exactly that instead of claiming a held-out guarantee (a spy measured a first origin at 13 with a model fitted through 19 – six future outcomes exposed). The clean alternative is refitting per origin, which this API deliberately does not fake by silently truncating your model.

Returns:

A PriceForecast with the calibrated (lo, hi) intervals, the point forecast mean, the raw per-step predictive draws paths (for Monte-Carlo DCF scenario analysis downstream), and level.

Return type:

PriceForecast

What (lo, hi) is: one MARGINAL interval PER LEAD, each with coverage AT LEAST level at its own horizon depth – NOT a simultaneous band, and not exact: split conformal is a lower-bound guarantee, so each lead can be conservative. Whole-path coverage is therefore NOT controlled in either direction by the per-lead guarantee: independent misses at exactly nominal would compound toward level**horizon, while conservative, strongly dependent leads can push joint coverage ABOVE level (an earlier version of this text said joint coverage “never exceeds” the requested level, which is false for a lower-bound guarantee – STAT-RR21-13). A plotted “tube” must not be read as containing the path with level probability in either direction; for path-level statements, work from paths directly.

Read calibration_count before trusting the exact level. At the defaults the reserved window yields on the order of 20 residuals per lead; split conformal’s finite-sample quantile moves in steps of 1/(n_cal + 1), so small counts make the interval coarse and typically conservative – and consecutive rolling origins are one step apart, so their horizon-step windows OVERLAP and the residuals are serially dependent: the exchangeability behind the conformal guarantee holds only approximately across overlapping origins. More history (or a larger cal_frac) buys both more residuals and less overlap distortion.

class Receipt(answer, produced_by='', ledger=None, trace=None, calibration=None, provenance=<factory>, policy=<factory>, executables=<factory>, bindings=<factory>, receipt_digest='', schema_version='mixle-answer-receipt-v1')

Bases: object

The bound artifact: an answer plus everything needed to re-verify it offline.

Parameters:

• answer (Any)

• produced_by (str)

• ledger (Explanation | None)

• trace (ExecutionTrace | None)

• calibration (dict[str, Any] | None)

• provenance (dict[str, Any])

• policy (dict[str, Any])

• executables (dict[str, str])

• bindings (dict[str, str])

• receipt_digest (str)

• schema_version (str)

to_json()

Return the answer receipt as JSON-compatible data.

Return type:

dict[str, Any]

to_knowledge_dict(*, id, project_id, task)

A plain dict shaped like mixle_knowledge.contracts.AnswerReceipt (id/project_id/task/ produced_by/answer/ledger/trace/calibration/provenance), aligned with the mixle-knowledge receipt contracts. Distinct from mixle_knowledge.contracts.ArtifactReceipt, which certifies a trained model/artifact, not one served answer – this is the per-answer evidence trail an offline consumer re-verifies (recompute the ledger, replay the trace, resolve the citations).

Stays a plain dict on purpose: mixle core carries no dependency on mixle-knowledge (platform contract packages depend on core, never the other way); constructing the validated pydantic object (AnswerReceipt(**receipt.to_knowledge_dict(...))) is the receiving side’s job.

Parameters:

• id (str)

• project_id (str)

• task (str)

Return type:

dict[str, Any]

class VerificationReport(checks=<factory>, required=())

Bases: object

Which claims were observed and which required observations determine the verdict.

Parameters:

• checks (dict[str, str])

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

property passed: bool

Whether at least one evidence claim was observed and every required check passed.

summary()

Return a compact comma-separated verification summary.

Return type:

str

verify_receipt(receipt, *, tools=None, evidence=None, tol=1e-9)

Re-check every claim the receipt actually makes, using only the receipt’s own bound data (plus the tools registry needed to re-execute a trace – the one piece that cannot be inlined into the receipt itself, since a tool is a function, not data).

Parameters:

• receipt (Receipt)

• tools (dict[str, Any] | None)

• evidence (dict[str, Any] | None)

• tol (float)

Return type:

VerificationReport

fit_and_record(data, estimator, *, seed, max_its=10, delta=1.0e-9, fit_options=None)

Execute one materialized fit and return the model with its complete replay receipt.

Parameters:

• data (Any)

• estimator (Any)

• seed (int)

• max_its (int)

• delta (float | None)

• fit_options (dict[str, Any] | None)

Return type:

tuple[Any, ReproReceipt]

record_fit(model, data, *, seed, estimator=None, max_its=10, delta=1.0e-9, fit_options=None)

Record a ReproReceipt for a model fitted on data with seed (see module docstring).

Pass the same max_its/delta actually used for the fit being recorded (they default to optimize()’s own defaults); verify_reproducible() replays the fit with these exact values.

Parameters:

• model (Any)

• data (Any)

• seed (int)

• estimator (Any)

• max_its (int)

• delta (float | None)

• fit_options (dict[str, Any] | None)

Return type:

ReproReceipt

verify_reproducible(estimator, data, receipt, *, seed=None, max_its=None)

Refit estimator on data and check the fit reproduces receipt (data + parameters).

Returns {reproducible, data_matches, params_match, refit_fingerprint}. reproducible is True iff BOTH the data fingerprint and the refit’s parameter fingerprint match the receipt – i.e. the exact fit can be recovered from the recorded recipe. seed defaults to the receipt’s seed. max_its defaults to the receipt’s own recorded value (the iteration budget the original fit actually used); delta always replays the receipt’s recorded value – a refit run with a different iteration budget or convergence tolerance than the original fit can land at a different (if nearby) optimum and spuriously report reproducible: False for a fit that would have reproduced exactly under its own original settings.

Parameters:

• estimator (Any)

• data (Any)

• receipt (ReproReceipt)

• seed (int | None)

• max_its (int | None)

Return type:

dict[str, Any]

class ReproReceipt(data_fingerprint, n, seed, estimator, param_fingerprint, max_its=10, delta=1e-09, schema_version='mixle-repro-v2', data_digest='', model_digest='', estimator_implementation_digest='', estimator_state_digest='', optimizer_implementation_digest='', dependency_digest='', environment=<factory>, fit_policy=<factory>, fit_policy_digest='', execution_digest='', receipt_digest='')

Bases: object

The complete executable recipe and observed result needed to re-derive a fit.

max_its/delta must match whatever was actually passed to the original optimize() call that produced the fingerprinted model (they default to optimize()’s own defaults). Without these, verify_reproducible() cannot replay the same fit – an iterative estimator refit with a different iteration budget or convergence tolerance can land at a different (if nearby) optimum and spuriously report reproducible: False for a fit that would have reproduced exactly under its own original settings.

Parameters:

• data_fingerprint (str)

• n (int)

• seed (int)

• estimator (str)

• param_fingerprint (str)

• max_its (int)

• delta (float | None)

• schema_version (str)

• data_digest (str)

• model_digest (str)

• estimator_implementation_digest (str)

• estimator_state_digest (str)

• optimizer_implementation_digest (str)

• dependency_digest (str)

• environment (dict[str, Any])

• fit_policy (dict[str, Any])

• fit_policy_digest (str)

• execution_digest (str)

• receipt_digest (str)

as_dict()

Return the receipt as JSON-compatible data.

Return type:

dict[str, Any]

matches_data(data)

Whether data is the exact dataset this fit was recorded on.

Parameters:

data (Any)

Return type:

bool

matches_model(model)

Whether model has the exact parameters this receipt fingerprinted.

Parameters:

model (Any)

Return type:

bool

data_fingerprint(data, *, ndigits=_NDIGITS)

A stable hash of a training dataset (order-sensitive; floats rounded) – identifies the exact input.

Parameters:

• data (Any)

• ndigits (int)

Return type:

str

param_fingerprint(model, *, ndigits=_NDIGITS)

A stable hash of a fitted model’s parameters, via its to_json state (floats rounded).

Parameters:

• model (Any)

• ndigits (int)

Return type:

str

class BeliefState

Bases: ABC

A distribution over a latent, exposing a uniform query + update interface.

Realizations answer where they are defined: mean(), cov(), var(), sd(), entropy(), interval(), sample(), marginal(), and – the point of a belief state – update(), which returns a new belief conditioned on fresh evidence.

abstractmethod mean()

The posterior mean of the latent.

Return type:

ndarray

abstractmethod entropy()

The differential/Shannon entropy H[q] (nats) – watch it shrink as evidence arrives.

Return type:

float

abstractmethod sample(n=1, rng=None)

Draw n latent samples from the belief.

Parameters:

• n (int)

• rng (Any)

Return type:

ndarray

abstractmethod update(*args, **kwargs)

Return a new belief conditioned on fresh evidence (the assimilation step).

Parameters:

• args (Any)

• kwargs (Any)

Return type:

BeliefState

cov()

The posterior covariance (not defined for every realization).

Return type:

ndarray

var()

Per-coordinate posterior variance.

Return type:

ndarray

sd()

Per-coordinate posterior standard deviation.

Return type:

ndarray

interval(level=0.9)

Per-coordinate central credible interval at level – an (d, 2) array of [lo, hi].

Parameters:

level (float)

Return type:

ndarray

marginal(indices)

The belief restricted to a subset of latent coordinates.

Parameters:

indices (Any)

Return type:

BeliefState

class GaussianBelief(mean, cov)

Bases: BeliefState

A multivariate-Gaussian belief N(mean, cov) over a continuous latent.

Evidence is a linear-Gaussian observation y = H z + noise, noise ~ N(0, R); update() applies the exact Kalman measurement update (Joseph form, so the covariance stays symmetric positive-definite). fuse() combines two beliefs about the same latent as a product of Gaussian experts only with an explicit fusion contract. condition() does noiseless Gaussian conditioning on a coordinate subset.

Parameters:

• mean (Any)

• cov (Any)

property dim: int

Return the latent dimensionality.

mean()

Return a copy of the Gaussian mean vector.

Return type:

ndarray

cov()

Return a copy of the Gaussian covariance matrix.

Return type:

ndarray

entropy()

Return differential entropy of the Gaussian belief in nats.

Return type:

float

interval(level=0.9)

Return marginal central intervals for each coordinate.

Parameters:

level (float)

Return type:

ndarray

sample(n=1, rng=None)

Draw samples from the Gaussian belief.

Parameters:

• n (int)

• rng (Any)

Return type:

ndarray

update(H, y, R)

Kalman measurement update: condition on y = H z + noise, noise ~ N(0, R).

Parameters:

• H (Any) – (k, d) observation matrix (or (d,) / scalar for a single linear readout).

• y (Any) – (k,) observed value (or scalar).

• R (Any) – (k, k) observation-noise covariance (or (k,) diagonal / scalar).

Return type:

GaussianBelief

fuse(other, *, shared_prior=None, other_is_likelihood=False)

Fuse a posterior or an explicitly identified Gaussian likelihood message.

For two posteriors derived from the same prior, pass shared_prior; its information is subtracted once. If other is instead a likelihood factor represented by a Gaussian kernel, set other_is_likelihood=True. Omitting both is ambiguous and rejected.

Parameters:

• other (GaussianBelief)

• shared_prior (GaussianBelief | None)

• other_is_likelihood (bool)

Return type:

GaussianBelief

condition(indices, values)

Noiseless Gaussian conditioning: fix latent coordinates indices to values.

Returns the belief over the remaining coordinates. This is the exact R -> 0 limit of an observation that reads off those coordinates.

Parameters:

• indices (Any)

• values (Any)

Return type:

GaussianBelief

marginal(indices)

Return the marginal Gaussian belief over selected coordinates.

Parameters:

indices (Any)

Return type:

GaussianBelief

as_belief(obj, node=None)

Adapt any object exposing mean/cov (a FieldPosterior node, a fitted Gaussian, a ParameterPosterior) into a GaussianBelief.

node is forwarded when the source is node-addressable (e.g. FieldPosterior.mean(node) / .cov(node)); otherwise mean/cov are called with no argument.

Parameters:

• obj (Any)

• node (Any)

Return type:

GaussianBelief

class UncertaintyDecomposition(total, aleatoric, epistemic, kind)

Bases: object

A predictive uncertainty split into aleatoric + epistemic (summing to total).

kind is "entropy" (values in nats) or "variance" (values in the outcome’s squared units). Each field is a scalar for a single query point, or an array over query points. Array fields are owned read-only copies, so the total = aleatoric + epistemic identity the record asserts cannot be broken after the fact (MXR-080-1899).

Parameters:

• total (ndarray)

• aleatoric (ndarray)

• epistemic (ndarray)

• kind (str)

property fraction_epistemic: ndarray

Share of the total uncertainty that is epistemic (reducible by more data), in [0, 1].

item()

Collapse size-1 arrays to Python floats (convenience for single-point decompositions).

Return type:

UncertaintyDecomposition

decompose_uncertainty(*, probs=None, means=None, variances=None)

Front door: decompose a predictive into aleatoric + epistemic uncertainty.

Pass exactly one representation of the per-member predictive:

• probs=(M, ..., K) – categorical predictives -> BALD entropy split (decompose_entropy());

• means=(M, ...) (with optional variances=(M, ...)) – continuous predictives -> law-of-total-variance split (decompose_variance()).

To decompose a fitted model’s parameter uncertainty, build the members first with posterior_ensemble() (then predictive_distribution() for the discrete case).

Parameters:

• probs (Any)

• means (Any)

• variances (Any)

Return type:

UncertaintyDecomposition

decompose_entropy(member_probs)

BALD entropy split of a discrete predictive.

Parameters:

member_probs (Any) – array (M, ..., K) – M posterior draws / ensemble members, each a categorical predictive over K outcomes (optionally batched over query points in the middle axes). Rows need not be normalized; each is renormalized over the last axis.

Returns:

An UncertaintyDecomposition with kind="entropy" (nats). epistemic is the PLUG-IN mutual information H(mean) - mean H, clamped to >= 0 (non-negative in exact arithmetic; the clamp removes floating-point negatives, and total is recomputed as aleatoric + epistemic afterward so the asserted identity holds exactly – audit U-8). ESTIMATOR LIMITS (audit U-2): the plug-in MI over M members is structurally capped at log M whatever the true disagreement – at the permitted minimum M = 2 a true 2.08-nat mutual information reads at most 0.693 – and is biased low below the cap; read it as a lower bound at small M, and prefer M large enough that log M comfortably exceeds the disagreement you need to resolve.

Return type:

UncertaintyDecomposition

decompose_variance(member_means, member_vars=None)

Law-of-total-variance split of a continuous predictive.

Parameters:

• member_means (Any) – array (M, ...) – each member’s predictive mean E[y | theta_m].

• member_vars (Any) – array (M, ...) – each member’s predictive variance Var[y | theta_m]. If None, aleatoric noise is taken as zero (members are point predictors) and the decomposition reports only the epistemic spread of the means.

Returns:

An UncertaintyDecomposition with kind="variance" (or kind="variance-epistemic-only" when member_vars is None – total then contains NO aleatoric noise and is not a predictive variance; audit U-9).

Return type:

UncertaintyDecomposition

predictive_distribution(members, support)

Evaluate an iterable of fitted distributions over a discrete support -> (M, K) probs.

Each member’s log_density is evaluated at every point of support and softmax-normalized over the support, giving one categorical row per member. Feed the result to decompose_entropy().

PRECONDITION (audit U-7): the softmax treats every support point as carrying EQUAL measure. That is exact for a discrete family evaluated on its full support, and a grid approximation for a continuous density ONLY when the grid is uniform – a non-uniform grid silently reweights the density by the missing cell widths, and any downstream entropy inherits the distortion. Restrict a discrete member’s support at your own risk: the row then represents the conditional given that subset.

Parameters:

• members (Iterable[Any])

• support (Sequence[Any])

Return type:

ndarray

posterior_ensemble(param_post, build, n=200, rng=None)

Materialize n models from a parameter posterior – an ensemble representing q(theta|data).

build maps one parameter draw (whatever ParameterPosterior.sample() returns) to a fitted distribution, mirroring from_parameter_posterior(). The returned list is the “members” the decomposition integrates over – so epistemic uncertainty here is genuine parameter uncertainty, not just ensemble disagreement.

Parameters:

• param_post (Any)

• build (Callable[[Any], Any])

• n (int)

• rng (Any)

Return type:

list[Any]

class Clustering(representatives, probs, labels)

Bases: object

Samples grouped into equivalence classes: representatives, class probs, per-sample labels.

The three fields are one interlocking statement – probs[c] is the mass of the class whose member representatives[c] stands for, and labels[i] says which class sample i landed in – so the record owns them: representatives is a tuple and the arrays are write-locked copies (MXR-080-1899). Consumers in this repository only read them (index, iterate, zip, argmax, tolist), which is why the sequence can be a tuple without breaking a caller.

Parameters:

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

• probs (ndarray)

• labels (ndarray)

cluster_samples(samples, equivalent=None)

Group samples into equivalence classes under equivalent (default exact ==).

For discrete draws whose surface form varies but meaning does not – e.g. LLM generations (“Paris”, “It’s Paris.”, “The capital is Paris”) – pass a semantic equivalent (embedding similarity, an entailment check, normalized match). True single linkage is computed as the connected components of all pairwise equivalence edges.

Parameters:

• samples (Sequence[Any])

• equivalent (Callable[[Any, Any], bool] | None)

Return type:

Clustering

marginalize_meaning(items, equivalent=None, *, log_probs=None, weights=None)

The distribution over meanings = the string distribution marginalized over each meaning class.

A generative model puts probability on strings; a meaning is an equivalence class of strings (equivalent decides sameness). The probability of a meaning c is the pushforward under the quotient – you sum the string probabilities over the class: P(c) = sum_{s in c} P(s). This returns that marginal (as a Clustering, probs = P(c)).

How the per-string probability enters:

• log_probs – the model’s sequence log-probabilities log P(s) for each DISTINCT item; classes are combined by logsumexp and renormalized over the items you passed, so the result is the meaning marginal CONDITIONAL ON the observed strings, not the model’s full marginal (audit U-5) – unobserved tail mass is excluded, and repeated strings would each add their own P(s) (estimating sum P(s)^2, a different quantity), so duplicates are refused in this branch. Use this when the items are distinct strings whose probabilities you know – it corrects the counting form’s hidden “every string in a class is equiprobable” assumption.

• weights – explicit non-negative masses per item (summed within class).

• neither – uniform mass, i.e. counting: P(c) = count_c / N. For i.i.d. samples from the model this is the unbiased Monte-Carlo estimate of the same marginal.

Parameters:

• items (Sequence[Any])

• equivalent (Callable[[Any, Any], bool] | None)

• log_probs (Any)

• weights (Any)

Return type:

Clustering

semantic_entropy(samples, equivalent=None, *, log_probs=None, weights=None)

PLUG-IN entropy (nats) over the meaning classes of samples.

Sample a stochastic generator (an LLM at temperature) n times, marginalize the string distribution over meaning classes (marginalize_meaning()), and take the entropy of that marginal. High semantic entropy means the model disagrees with itself about what the answer is (a hallucination signal), as opposed to merely phrasing one answer many ways (which collapses to low entropy).

APPROXIMATION, NOT A MEASUREMENT OF THE GENERATOR (STAT-RR17-14): with n counted samples this is the plug-in estimator, biased LOW by about (K - 1) / (2n) nats (Miller-Madow, K = observed meaning classes) – at the eight-draw default an abstention loop actually uses, a uniform ten-class generator’s true 2.303 nats measured 1.646 on average (bias -0.656): the bias direction makes a HALLUCINATING model look confident, the exact failure the signal exists to catch. Use semantic_entropy_receipt() for the estimate WITH its sample size, observed class count, bias estimate, and the bias-corrected value; abstention thresholds should use the corrected value (it errs toward abstaining). Pass log_probs (the sequence log-likelihoods) to marginalize with the actual string probabilities rather than by sample counting. Feed the clusters’ per-member distributions to decompose_entropy() for an epistemic/aleatoric split.

Parameters:

• samples (Sequence[Any])

• equivalent (Callable[[Any, Any], bool] | None)

• log_probs (Any)

• weights (Any)

Return type:

float

reliability_curve(prob, outcome, *, bins=10, strategy='uniform', ci=False, n_boot=1000, ci_level=0.95, seed=0)

Reliability diagram: observed frequency vs mean forecast probability, per bin.

Bins the forecasts, then within each bin compares the mean predicted probability to the observed event frequency. A perfectly calibrated forecaster lies on the diagonal observed == predicted.

Parameters:

• prob (ndarray) – (n,) predicted probabilities of the positive class (or top-label confidences).

• outcome (ndarray) – (n,) 0/1 outcomes (or correctness indicators).

• bins (int) – number of bins.

• strategy (str) – "uniform" (equal-width) or "quantile" (equal-count) bins.

• ci (bool) – if True attach a percentile bootstrap interval on the observed frequency in each bin.

• n_boot (int) – bootstrap resamples when ci is True.

• ci_level (float) – central probability of each PER-BIN interval (e.g. 0.95).

• seed (int | RandomState | None) – RNG seed for the bootstrap.

Returns:

{'mean_pred', 'obs_freq', 'count', 'bin_edges'} (one entry per non-empty bin), plus 'obs_lo' / 'obs_hi' / 'obs_ci_effective_boot' when ci is True.

Return type:

dict[str, ndarray]

The intervals are POINTWISE, one bin at a time – with 10 bins at 95% each, at least one bin escaping its interval under perfect calibration is a ~30-45% event, so judge single-bin excursions accordingly. A resample in which a bin comes up EMPTY contributes nothing to that bin’s interval; the quantiles are taken over the non-empty resamples only, which conditions the interval on the bin being observed (a real bias for sparse bins). obs_ci_effective_boot reports, per bin, how many of the n_boot resamples actually informed the interval – treat an interval resting on far fewer than n_boot replicates as unstable rather than narrow.

expected_calibration_error(prob, outcome, *, bins=10, strategy='uniform', norm='l1', ci=False, n_boot=1000, ci_level=0.95, seed=0)

Expected Calibration Error: count-weighted average gap between confidence and accuracy.

ECE = sum_b (n_b / n) |obs_b - pred_b| over bins (norm='l2' uses the squared gap, square- rooted). For multiclass classifiers reduce with top_label_confidence() first.

Zero is NOT the null value. The plug-in ECE takes an absolute value of binomial noise in every bin, so a PERFECTLY calibrated forecaster still scores roughly sqrt(bins / n) in expectation – and a percentile bootstrap interval around the plug-in estimate inherits that upward shift, so under perfect calibration it excludes the true value 0 essentially always. Judge the point estimate against calibration_null_expectation() (the measured perfect-calibration distribution of this same statistic at your prob profile and binning), and read the bootstrap interval as sampling uncertainty around the BIASED plug-in functional, not as a test of calibration.

Parameters:

• prob (ndarray) – (n,) predicted probabilities / confidences.

• outcome (ndarray) – (n,) 0/1 outcomes / correctness indicators.

• bins (int) – number of bins.

• strategy (str) – "uniform" or "quantile" binning.

• norm (str) – "l1" (mean absolute gap) or "l2" (root mean squared gap).

• ci (bool) – if True also return a percentile bootstrap interval.

• n_boot (int) – bootstrap controls.

• ci_level (float) – bootstrap controls.

• seed (int | RandomState | None) – bootstrap controls.

Returns:

The ECE (float), or (ece, lo, hi) when ci is True.

Return type:

float | tuple[float, float, float]

maximum_calibration_error(prob, outcome, *, bins=10, strategy='uniform')

Maximum Calibration Error: the worst per-bin gap max_b |obs_b - pred_b|.

Unlike expected_calibration_error() this is not count-weighted, so it surfaces a small but badly-miscalibrated region that the average would hide – and for exactly that reason its null value GROWS with bins: a max over more (and therefore sparser) bins of pure binomial noise rises from ~0.05 to ~0.6 as bins goes 10 -> 500 under perfect calibration, with the argmax typically landing on the sparsest bin. A large MCE from a thin bin is what noise looks like at that binning; compare to calibration_null_expectation() before reading it as a defect, and prefer more data per bin (fewer bins, or quantile binning) when the max matters.

Parameters:

• prob (ndarray)

• outcome (ndarray)

• bins (int)

• strategy (str)

Return type:

float

top_label_confidence(prob, labels)

Reduce a multiclass classifier to the (confidence, correct) top-label calibration problem.

This checks calibration OF THE ARGMAX CLASS ONLY – a deliberately narrow reduction, not THE multiclass calibration check. Probabilities the model puts on classes it never predicts are invisible to it: a forecaster that always assigns 0.4 to a class that never occurs (while its argmax confidences are accurate) passes top-label calibration with ECE 0. Classwise or full- vector calibration is a strictly stronger property and is not implemented here; when the non-argmax probabilities carry decisions (thresholded alerts on a minority class, expected-cost rankings), check those probabilities against their own outcomes directly.

Parameters:

• prob (ndarray) – (n, K) class-probability matrix.

• labels (ndarray) – (n,) integer true labels.

Returns:

the max predicted probability per row and a 0/1 indicator of whether that argmax class was the true label. Feed these to reliability_curve() / expected_calibration_error().

Return type:

(confidence, correct)

class ProbabilityCalibrator(method='isotonic')

Bases: object

Map raw scores to calibrated probabilities – fit against binary outcomes.

A raw score (a model’s confidence, a self-consistency fraction, a token likelihood) need not be a probability of anything: it can be monotone-but-miscalibrated, or have no relationship to the outcome at all. This learns the transform score -> P(outcome = 1 | score) from labeled data, so the output is a probability of the event you calibrated against.

• method="isotonic" – monotone, non-parametric (pool-adjacent-violators). Assumes higher score => not-lower probability; flexible, needs enough calibration points.

• method="platt" – logistic sigmoid(a * score + b). Two parameters, usable on little data, but assumes a sigmoidal relationship. Fit against Platt’s smoothed targets (N+ + 1)/(N+ + 2) and 1/(N- + 2) rather than raw 0/1: with raw targets, perfectly separable calibration data drives a to infinity and the “calibrated” outputs to exact 0/1 – maximal overconfidence from the method that was supposed to remove it. The smoothing keeps the map finite and the outputs strictly inside (0, 1).

The guarantee is conditional on HELD-OUT evaluation. The isotonic fit interpolates its own calibration outcomes, so re-checking calibration on the SAME data it was fit on passes by construction (a reliability curve of the fitted values against the fitting outcomes sits on the diagonal no matter how little signal the score carried); any honest check of a fitted calibrator – isotonic especially – uses data it never saw. Both methods assume the score/outcome relationship is stable between fitting and serving.

A near-flat fitted curve is itself the finding: it means the raw score carried little information about the outcome (its “likelihood” was unrelated to the event).

Parameters:

method (str)

fit(scores, outcomes)

Fit the score->probability map on scores with binary outcomes (0/1).

Parameters:

• scores (Any)

• outcomes (Any)

Return type:

ProbabilityCalibrator

predict(scores)

Calibrated probabilities for scores (clamped to [0, 1]).

Parameters:

scores (Any)

Return type:

ndarray

calibrate_probabilities(scores, outcomes, *, method='isotonic')

Fit a ProbabilityCalibrator mapping scores to P(outcome=1 | score).

Parameters:

• scores (Any)

• outcomes (Any)

• method (str)

Return type:

ProbabilityCalibrator

calibration_null_expectation(prob, *, bins=10, strategy='uniform', norm='l1', n_sim=500, seed=0)

What ECE and MCE look like for a PERFECTLY calibrated forecaster at this prob profile.

Simulates outcomes y*_i ~ Bernoulli(prob_i) – the definition of perfect calibration at the observed forecasts – and recomputes the binned statistics each time, with the same bins the real statistics would use. The plug-in ECE/MCE are strictly positive under this null (absolute values of binomial noise), so THESE numbers, not zero, are the calibrated baseline: an observed ECE near ece here is consistent with perfect calibration, and evidence of miscalibration starts around the *_q95 quantiles, which give a one-sided 5% test by construction.

Parameters:

• prob (ndarray) – (n,) the forecaster’s predicted probabilities (outcomes are not needed – the null generates them).

• bins (int) – the same binning controls you pass to the real statistics.

• strategy (str) – the same binning controls you pass to the real statistics.

• norm (str) – the same binning controls you pass to the real statistics.

• n_sim (int) – null simulations; at least 20. The *_q95 fields are the CONSERVATIVE order statistic X_(ceil(0.95 (n_sim+1))), calibrated for the STRICT comparison observed > q95: for a continuous statistic its fresh-null strict exceedance is exactly 1 - k/(n_sim + 1) <= 5% at every accepted n_sim (4.76% at the floor of 20). With DISCRETE ties the two operators split – a single forecast at p = 0.5 makes every null ECE equal the threshold, so strict exceedance is 0 while >= reads 1 (STAT-RR22-10) – so compare with > and read the bound as <= (ties only make it more conservative). The interpolated quantile realized 9.83% at the floor, and n_sim=1 once labeled a single draw a “q95” with 49.8% exceedance (STAT-RR19-10/STAT-RR21-09). The default 500 keeps the threshold’s own MC noise small.

• seed (int | RandomState | None) – RNG seed.

Returns:

{'ece', 'ece_q95', 'mce', 'mce_q95', 'n_sim'} – null means and 95th percentiles of the plug-in ECE (under norm) and MCE at this forecast profile and binning.

Return type:

dict[str, float]

pit_values(y, cdf)

Probability Integral Transform values u_i = F_i(y_i).

Under a calibrated continuous predictive distribution the PIT values are Uniform(0, 1). Pass either the precomputed CDF values F_i(y_i) or a callable cdf(y) -> F(y); a scalar-only callable such as a fitted univariate distribution’s .cdf is evaluated element-wise (mixtures and other structured families define no cdf; pass their CDF as your own callable).

Parameters:

• y (ndarray) – (n,) realised values.

• cdf (ndarray | Callable[[ndarray], ndarray]) – (n,) precomputed predictive-CDF values at y, or a callable applied to y.

Returns:

(n,) PIT values in [0, 1]. Invalid CDF values raise rather than being clipped.

Return type:

ndarray

pit_ensemble(y, forecasts, *, randomize=True, seed=0)

Rank-based PIT from a finite predictive ensemble.

u_i is the fraction of ensemble members <= y_i. With randomize=True ties are broken by a uniform jitter within the rank gap, which makes the PIT exactly Uniform(0, 1) under calibration even for discrete ensembles (the randomized PIT of Czado et al. 2009).

Parameters:

• y (ndarray) – (n,) realised values.

• forecasts (ndarray) – (n, m) ensemble (m draws per observation).

• randomize (bool) – jitter ties for an exactly-uniform PIT.

• seed (int | RandomState | None) – RNG seed when randomize is True.

Returns:

(n,) PIT values in [0, 1].

Return type:

ndarray

pit_histogram(pit, *, bins=10)

Histogram of PIT values with the uniform reference level.

Parameters:

• pit (ndarray) – (n,) PIT values.

• bins (int) – number of equal-width bins on [0, 1].

Returns:

{'counts', 'density', 'edges', 'uniform'} where density integrates to 1 and uniform is the flat reference density (1.0) a calibrated forecast would match.

Return type:

dict[str, ndarray]

pit_calibration_error(pit, *, bins=10)

Calibration error of a PIT histogram: TOTAL absolute deviation from uniform mass.

sum_b |count_b/n - 1/bins| – the SUM over bins (twice the total-variation distance to uniform, range [0, 2(1 - 1/bins)]), not a per-bin mean: a threshold sized for a mean would be off by a factor of bins. 0 when the PIT histogram is perfectly flat (calibrated), larger when it is U-shaped (under-dispersed) or humped (over-dispersed).

Parameters:

• pit (ndarray)

• bins (int)

Return type:

float

interval_coverage(lower, upper, y)

Empirical coverage and mean width of a set of prediction intervals.

Parameters:

• lower (ndarray) – (n,) lower endpoints.

• upper (ndarray) – (n,) upper endpoints.

• y (ndarray) – (n,) realised values.

Returns:

{'coverage', 'mean_width'} – the fraction of y inside [lower, upper] and the mean interval width. Compare coverage to the nominal level the interval was built for.

Return type:

dict[str, float]

coverage_curve(forecasts, y, *, levels=None)

Empirical-vs-nominal coverage of central intervals across a grid of nominal levels.

For each nominal central level c the per-observation central interval [quantile((1-c)/2), quantile((1+c)/2)] is read off the predictive ensemble and its empirical coverage is measured.

Judge empirical against null_expectation, NOT the diagonal. Plug-in quantiles of a finite ensemble cover BELOW nominal even when the ensemble is perfect: with m = 5 members, a perfectly calibrated forecast’s 95% interval covers only ~0.64 – the interval simply cannot reach past the sample extremes. null_expectation is the expected coverage of a perfectly calibrated m-member ensemble at each level, computed for a uniform reference with the same quantile rule. It is almost – not exactly – distribution-free: linear quantile interpolation leaves a small dependence on the ensemble’s shape at very small m (measured 0.632 / 0.644 / 0.636 for uniform / Gaussian / exponential at m = 5 and nominal 0.95, converging by m ~ 50), which is an order of magnitude smaller than the diagonal gap it corrects. empirical below the DIAGONAL at small m is expected; empirical meaningfully below null_expectation indicates over-confidence.

Parameters:

• forecasts (ndarray) – (n, m) predictive ensemble.

• y (ndarray) – (n,) realised values.

• levels (ndarray | None) – nominal central coverage levels in (0, 1); defaults to 0.05 .. 0.95 by 0.05.

Returns:

{'nominal', 'empirical', 'null_expectation'} arrays of equal length – null_expectation is the finite-m perfect-calibration benchmark per level.

Return type:

dict[str, ndarray]

bonferroni(pvals, *, alpha=0.05)

Bonferroni single-step FWER control: adjusted p_i = min(1, m p_i).

The simplest and most conservative correction; always valid regardless of dependence.

Parameters:

• pvals (ndarray) – (m,) raw p-values.

• alpha (float) – target family-wise error rate, in (0, 1).

Returns:

{'reject', 'pvals_adjusted', 'n_reject', 'alpha'}.

Return type:

dict

holm(pvals, *, alpha=0.05)

Holm step-down FWER control – uniformly more powerful than Bonferroni, equally valid.

Sorts p-values ascending and applies the running factor m - k with a cumulative maximum so the adjusted values stay monotone.

Parameters:

• pvals (ndarray)

• alpha (float)

Return type:

dict

hochberg(pvals, *, alpha=0.05)

Hochberg step-up FWER control (valid under independence / positive dependence).

Same per-step factor as Holm but applied step-up (from the largest p-value down), giving more rejections; requires the independence/PRDS assumption that Holm does not.

Parameters:

• pvals (ndarray)

• alpha (float)

Return type:

dict

benjamini_hochberg(pvals, *, alpha=0.05)

Benjamini–Hochberg FDR control; adjusted values are q-values.

Controls the expected false-discovery proportion at alpha under independence or positive regression dependence (PRDS). The standard choice for screening many hypotheses.

Parameters:

• pvals (ndarray)

• alpha (float)

Return type:

dict

benjamini_yekutieli(pvals, *, alpha=0.05)

Benjamini–Yekutieli FDR control, valid under arbitrary dependence.

Like benjamini_hochberg() but inflated by the harmonic factor c(m) = sum_{i=1}^m 1/i, so it holds for any dependence structure at the cost of power.

Parameters:

• pvals (ndarray)

• alpha (float)

Return type:

dict

adjust_pvalues(pvals, *, method='bh', alpha=0.05)

Unified dispatcher over the correction methods.

Parameters:

• pvals (ndarray) – (m,) raw p-values.

• method (str) – one of "bonferroni", "holm", "hochberg", "bh" (Benjamini–Hochberg), "by" (Benjamini–Yekutieli).

• alpha (float) – target error rate (FWER for the first three, FDR for the last two), in (0, 1).

Returns:

{'reject', 'pvals_adjusted', 'n_reject', 'alpha'}.

Return type:

dict

fisher_combine(pvals)

Fisher’s method: combine independent p-values via -2 sum log p ~ chi^2_{2k}.

Sensitive to a few very small p-values. For combining evidence for the same hypothesis across independent tests.

Returns:

{'statistic', 'pvalue', 'df'}.

Parameters:

pvals (ndarray)

Return type:

dict[str, float]

stouffer_combine(pvals, *, weights=None)

Stouffer’s Z method: combine p-values on the z-scale, optionally weighted.

Z = sum w_i Phi^{-1}(1 - p_i) / sqrt(sum w_i^2). Weights let more-precise studies (e.g. larger samples) count more; equal weights recover the unweighted combination. A weight of exactly 0 for one study while the others are positive is a legitimate “exclude this study” case (mirrors weighted_conformal()’s zero-weight convention) and stays valid; every weight being non-positive or non-finite together is degenerate and raises.

Parameters:

• pvals (ndarray) – (k,) raw p-values from independent studies/strata.

• weights (ndarray | None) – optional (k,) weights, non-negative with a positive total; None weights every study equally.

Returns:

{'z', 'pvalue'} (one-sided combined p-value).

Return type:

dict[str, float]

tippett_combine(pvals)

Tippett’s method (Sidak min-p): combined p = 1 - (1 - min p)^k.

Most powerful when a single strong signal exists among the tests.

Evaluated as -expm1(k * log1p(-min_p)) rather than literally (MXR-080-1604). Forming 1 - min_p first rounds the evidence away before it is ever exponentiated: below the float64 epsilon the subtraction returns exactly 1.0, so a minimum p-value of 1e-17, 1e-20 or 1e-100 combined to exactly 0.0 – an impossible result for finite nonzero inputs, and one that overstates certainty without bound. Even 1e-16 was inflated to 2.220446049250313e-16 (two epsilons) instead of the correct 2e-16. The stable form is accurate all the way down to the smallest subnormal.

Returns:

{'min_p', 'pvalue'}.

Parameters:

pvals (ndarray)

Return type:

dict[str, float]

bootstrap(data, statistic, *, n_boot=2000, method='bca', ci_level=0.95, seed=0, groups=None, clusters=None, block_length=None, m=None)

Bootstrap confidence interval for statistic(data).

Parameters:

• data (Any) – a single array (resampled along axis 0) or a tuple of arrays sharing their first axis (the statistic is then called as statistic(*parts)).

• statistic (Callable[[...], Any]) – maps the data to a scalar or fixed-length vector.

• n_boot (int) – number of bootstrap resamples.

• method (str) – "percentile", "basic" (pivotal), or "bca" (bias-corrected & accelerated). "bca"’s second-order accuracy holds for plain i.i.d. resampling of a SMOOTH (asymptotically normal, differentiable) functional; with groups / clusters / block_length / m set it falls back to "percentile". Non-smooth statistics – sample quantiles and the median above all – violate the smoothness the jackknife acceleration needs (the leave-one-out vector collapses to a couple of distinct values); when that degeneracy is detected the interval falls back to "percentile" and method records "percentile (bca-degenerate-jackknife)".

• ci_level (float) – central probability of the interval.

• seed (int | RandomState | None) – RNG seed.

• groups (ndarray | None) – (n,) labels for stratified resampling (resample within each group).

• clusters (ndarray | None) – (n,) labels for cluster resampling (resample whole clusters with replacement). Validity is asymptotic in the NUMBER OF CLUSTERS, not the number of rows: at G = 5 clusters, nominal-95% coverage measures ~0.5-0.7 however many rows each cluster holds, so treat few-cluster intervals as descriptive. Unequal cluster sizes additionally inflate replicate variance (each resample’s row count varies with which clusters were drawn).

• block_length (int | None) – moving-block length for serially dependent (time-series) data. Blocks are CIRCULAR (the series wraps), which removes the endpoint under-weighting of non-circular moving blocks. For picking the length, the coverage-optimal rate grows with the sample – on the order of n**(1/3) scaled by the dependence strength – so “the correlation length” alone is materially too short for long series: a length that ignores n leaves residual dependence between blocks uncounted.

• m (int | None) – subsample size for m-out-of-n subsampling (without replacement), m < n. Replicates are rescaled about the point estimate by sqrt(m/n) (Politis–Romano), which assumes the usual sqrt(n)-consistent statistic (a faster-rate statistic gets a conservative interval), and the interval takes the BASIC (pivotal) orientation that subsampling theory yields – method is overridden in this mode.

Returns:

A BootstrapResult.

Return type:

BootstrapResult

class BootstrapResult(estimate, ci_low, ci_high, distribution, method, ci_level, standard_error)

Bases: object

Result of a bootstrap() call.

Parameters:

• estimate (ndarray)

• ci_low (ndarray)

• ci_high (ndarray)

• distribution (ndarray)

• method (str)

• ci_level (float)

• standard_error (ndarray)

estimate

the statistic on the original data (scalar or vector).

Type:

numpy.ndarray

ci_low / ci_high

confidence-interval endpoints (same shape as estimate).

distribution

(n_boot, ...) array of bootstrap replicates.

Type:

numpy.ndarray

method

the interval method used.

Type:

str

ci_level

the central probability of the interval.

Type:

float

standard_error

bootstrap standard error (std of the replicates).

Type:

numpy.ndarray

block_bootstrap(data, statistic, block_length, *, n_boot=2000, ci_level=0.95, seed=0)

Moving-block bootstrap for serially dependent (time-series) data.

Convenience wrapper over bootstrap() with block_length set: resamples contiguous blocks so within-block autocorrelation is preserved. Choose block_length on the order of the series’ correlation length.

Parameters:

• data (Any)

• statistic (Callable[[...], Any])

• block_length (int)

• n_boot (int)

• ci_level (float)

• seed (int | RandomState | None)

Return type:

BootstrapResult

wild_bootstrap(fitted, residuals, statistic, *, n_boot=2000, kind='rademacher', ci_level=0.95, seed=0, leverage=None)

Wild (residual-multiplier) bootstrap for heteroscedastic regression errors.

Builds synthetic responses y* = fitted + residual * v where v are mean-zero, unit-variance two-point multipliers drawn independently per observation, then recomputes the statistic on each y*. Because each residual keeps its own magnitude, the procedure preserves heteroscedasticity that an i.i.d. residual resample would destroy.

The robustness claim is conditional on what the residuals carry. OLS residuals systematically UNDERSTATE the error at influential design points (E[e_i^2] = sigma_i^2 (1 - h_i) with h_i the hat-matrix diagonal), so multiplying raw residuals hands high-leverage observations too little noise exactly where mis-estimated variance hurts the most. Pass leverage (the hat diagonals) to apply the Davidson-Flachaire restoration e_i / (1 - h_i) before multiplication – the standard wild-bootstrap practice for regression coefficients. residuals must be RESPONSE-SCALE y - fitted: Pearson, deviance, or standardized residuals reconstruct a different response vector and silently change the estimand (the point estimate is computed from fitted + residuals, which equals the observed y only for raw residuals).

Parameters:

• fitted (ndarray) – (n,) fitted values from the model.

• residuals (ndarray) – (n,) raw response-scale residuals y - fitted.

• statistic (Callable[[ndarray], Any]) – maps a synthetic response vector y* to a scalar or vector (e.g. refit and return coefficients).

• n_boot (int) – number of resamples.

• kind (str) – "rademacher" (v in {-1, +1}) or "mammen" (Mammen’s two-point distribution).

• ci_level (float) – central probability of the percentile interval.

• seed (int | RandomState | None) – RNG seed.

• leverage (ndarray | None) – optional (n,) hat-matrix diagonals h_i in [0, 1); when given, residuals are inflated to e_i / (1 - h_i) inside the resampling loop (the point estimate still uses the raw residuals, i.e. the observed response).

Returns:

A BootstrapResult (percentile interval).

Return type:

BootstrapResult

permutation_test(x, y, *, statistic=None, n_perm=10000, alternative='two-sided', paired=False, stratify=None, seed=0, exact_max=10000, null_value=None)

Two-sample permutation test for an arbitrary statistic under a sharp null.

Under the null that the two samples are exchangeable, the labels can be shuffled freely; the statistic’s permutation distribution is the reference. For two-sided the statistic is centered at zero by construction (difference statistics) and compared on absolute value.

Parameters:

• x (ndarray) – the two samples (1-D), each with at least one observation. For paired=True they must have equal length and pairing is preserved by sign-flipping the within-pair differences.

• y (ndarray) – the two samples (1-D), each with at least one observation. For paired=True they must have equal length and pairing is preserved by sign-flipping the within-pair differences.

• statistic (Callable[[ndarray, ndarray], float] | None) – f(x, y) -> float; defaults to the difference in means. For paired it is applied to (differences, zeros) so the default reduces to the mean paired difference.

• n_perm (int) – number of random rearrangements (ignored if the exact set is enumerated).

• alternative (str) – "two-sided", "greater", or "less".

• paired (bool) – paired (sign-flip) permutation instead of label shuffling.

• stratify (ndarray | None) – (n,) group labels (concatenated x-then-y) for restricted permutation – labels are shuffled only within each group, preserving group structure.

• seed (int | RandomState | None) – RNG seed.

• exact_max (int) – if the number of distinct rearrangements is <= exact_max they are enumerated for an exact p-value.

• null_value (float | None) – reference value for a custom two-sided statistic. Required when statistic is supplied with alternative="two-sided" so the absolute tail is centered correctly.

Returns:

A PermutationResult.

Raises:

ValueError – if x or y is empty. A single observation per group (or a single paired difference) is still a legitimate, if weak, test – e.g. one-vs-one has 2 distinct label swaps and correctly returns pvalue=1.0, not an error. But an empty group has ZERO rearrangements to enumerate: exact enumeration falls back to comb(n, 0) == 1, silently “enumerating” one degenerate combination that recomputes stat on the same empty-vs- everything split, producing a null distribution of a single NaN. NaN >= NaN is False in numpy, so the observed-vs-null comparison in _pvalue() counts 0 exceedances out of 1 and returns an exact p-value of 0.0 – reading as maximal significance from zero evidence instead of failing loudly.

Return type:

PermutationResult

class PermutationResult(statistic, pvalue, null_distribution, n_perm, exact, alternative, null_hypothesis='')

Bases: object

Result of a permutation_test().

Parameters:

• statistic (float)

• pvalue (float)

• null_distribution (ndarray)

• n_perm (int)

• exact (bool)

• alternative (str)

• null_hypothesis (str)

statistic

the observed test statistic.

Type:

float

pvalue

the (one- or two-sided) p-value.

Type:

float

null_distribution

the statistic under each sampled/enumerated rearrangement.

Type:

numpy.ndarray

n_perm

number of rearrangements used.

Type:

int

exact

True if the full permutation set was enumerated.

Type:

bool

alternative

the alternative hypothesis.

Type:

str

sandwich_covariance(scores, bread, *, clusters=None, n_params=None, small_sample=True)

Generic sandwich covariance B M B' from per-observation scores and the bread.

Parameters:

• scores (ndarray) – (n, p) per-observation contributions to the estimating equation (for OLS the i-th row is x_i * e_i; for a GLM the score of observation i).

• bread (ndarray) – (p, p) inverse sensitivity B (e.g. (X'X)^{-1} for OLS, the inverse negative Hessian / inverse Fisher for an M-estimator).

• clusters (ndarray | None) – optional (n,) cluster labels; scores are summed within cluster before forming the meat (the robust-to-within-cluster-correlation meat). If None, observations are independent.

• n_params (int | None) – number of estimated parameters for the small-sample correction (defaults to p).

• small_sample (bool) – apply the usual finite-sample correction (n/(n-p) for independent data, G/(G-1) * (n-1)/(n-p) for clustered).

Returns:

The (p, p) robust covariance matrix.

Raises:

ValueError – if clusters is not aligned with scores, or holds fewer than 2 distinct clusters (across-cluster variation is what the clustered meat is estimated from), or n_params is not an exact nonnegative integer.

Return type:

ndarray

ols_robust_covariance(x, residuals, *, hc='hc1')

Heteroscedasticity-consistent (White) covariance for OLS coefficients.

Cov = (X'X)^{-1} ( sum_i w_i x_i x_i' ) (X'X)^{-1} where the per-observation weight w_i is the squared residual, optionally adjusted for leverage h_ii:

• hc0: e_i^2 (White’s original; biased down in small samples).

• hc1: e_i^2 * n/(n-p) (degrees-of-freedom correction; the common default).

• hc2: e_i^2 / (1 - h_ii).

• hc3: e_i^2 / (1 - h_ii)^2 (best small-sample behaviour; ~jackknife).

Parameters:

• x (ndarray) – (n, p) design matrix (include an intercept column if wanted).

• residuals (ndarray) – (n,) OLS residuals y - X beta_hat.

• hc (str) – one of "hc0", "hc1", "hc2", "hc3".

Returns:

The (p, p) robust covariance matrix.

Return type:

ndarray

cluster_robust_covariance(x, residuals, clusters, *, small_sample=True)

Cluster-robust (one-way or multi-way) covariance for OLS coefficients.

Allows arbitrary correlation within clusters while assuming independence across them. Pass one label array for one-way clustering, or a list/tuple of label arrays for multi-way clustering, which uses the Cameron–Gelbach–Miller inclusion–exclusion V_A + V_B - V_{A∩B} (two-way) and its higher-way generalisation.

Parameters:

• x (ndarray) – (n, p) design matrix.

• residuals (ndarray) – (n,) OLS residuals.

• clusters (ndarray | list) – (n,) labels, or a list of such arrays for multi-way clustering.

• small_sample (bool) – apply the G/(G-1) * (n-1)/(n-p) finite-sample correction per term.

Returns:

The (p, p) cluster-robust covariance matrix.

Raises:

ValueError – if clusters is an empty list of dimensions, is not aligned with x, or any clustering (or intersection clustering) has fewer than 2 distinct clusters.

Return type:

ndarray

newey_west_covariance(x, residuals, *, lags=None, small_sample=True)

Newey–West HAC covariance for OLS coefficients (serially correlated errors).

Heteroscedasticity- and autocorrelation-consistent: the meat is the long-run covariance of the scores s_t = x_t e_t, estimated with Bartlett (triangular) weights so it stays positive semi-definite:

S = Gamma_0 + sum_{l=1}^{L} (1 - l/(L+1)) (Gamma_l + Gamma_l’), Gamma_l = sum_t s_t s_{t-l}’.

Rows of x (and residuals) are assumed to be in time order.

Parameters:

• x (ndarray) – (n, p) design matrix, rows in time order.

• residuals (ndarray) – (n,) residuals in time order.

• lags (int | None) – truncation lag L; defaults to the Newey–West rule floor(4 (n/100)^{2/9}).

• small_sample (bool) – apply the n/(n-p) degrees-of-freedom correction.

Returns:

The (p, p) HAC covariance matrix.

Return type:

ndarray

robust_standard_errors(cov)

Standard errors sqrt(diag(cov)) from a covariance matrix.

A negative variance is not a small standard error – it means the covariance estimate is not a covariance. Clipping it to zero would turn an invalid, indefinite, or numerically failed estimate (an unreplicated multi-way cluster meat, a bread that did not invert) into the strongest possible claim the API can make: perfect precision. That is the one wrong answer a caller cannot detect, so a materially negative variance is rejected instead.

Only roundoff-scale negatives are healed: a variance whose magnitude is within 1e-10 of the covariance’s own diagonal scale is genuine floating-point noise around an exact zero (a perfectly-collinear direction, a zero-residual fit) and is floored to zero as before.

Raises:

ValueError – if the diagonal contains a non-finite value or a materially negative variance.

Parameters:

cov (ndarray)

Return type:

ndarray

glm(x, y, *, family='gaussian', link=None, offset=None, weights=None, weight_type='analytic', max_iter=100, tol=1e-8, robust=False, theta=1.0)

Fit a generalized linear model by iteratively reweighted least squares.

Parameters:

• x (ndarray) – (n, p) design matrix (include an intercept column explicitly if wanted).

• y (ndarray) – (n,) response (counts, 0/1 or proportions, positive reals, … per the family).

• family (str | Family) – "gaussian", "binomial", "poisson", "gamma", "inverse_gaussian", "negativebinomial", or a Family.

• link (str | Link | None) – a link name, or a Link instance to use directly (e.g. a custom link); defaults to the family’s default link (the canonical link, except gamma / inverse-Gaussian which default to log – their canonical inverse / inverse_squared links do not keep mu positive).

• offset (ndarray | None) – (n,) known additive term on the linear-predictor scale (e.g. log exposure).

• weights (ndarray | None) – (n,) non-negative prior weights, at least one positive. Zero-weight rows drop out of the fit entirely. What a weight MEANS is set by weight_type – the two readings agree on the coefficients and on nothing else.

• weight_type (str) –

  the convention weights are read under.

  • "analytic" (default; also called prior, precision or variance weights – the McCullagh–Nelder reading, and R’s for glm): row i is ONE observation measured with precision w_i, so Var(y_i) = phi V(mu_i) / w_i. The sample size is the number of rows: residual_df = rows - rank, dispersion divides by that, and the log-likelihood is that of n independent observations with dispersions phi/w_i. Where the dispersion is ESTIMATED, phi then absorbs the weights’ scale, so multiplying every weight by a constant leaves the whole fit unchanged – coefficients, standard errors, log-likelihood, AIC and BIC alike – and only the relative precisions matter. (Not so where phi is fixed at 1: for binomial / Poisson / negative-binomial the weight scale IS the amount of information, and halving every weight widens the standard errors by sqrt(2).) Use this for inverse-variance weighting, measurement precisions or sampling weights; fractional weights are accepted.

  • "frequency": w_i is a COUNT – row i stands for w_i identical observations – so it must be a whole number. Every reported quantity then equals what fitting the expanded sum(w)-row data set gives: residual_df = sum(w) - rank, dispersion divides by that, robust standard errors replicate the score, and AIC/BIC span sum(w) observations. Use it for aggregated frequency tables – survey cells, contingency tables, “distinct rows plus a count” extracts.

  The gap is not cosmetic. On a 6-row table whose counts sum to 114, the same call gives se = [0.369, 0.135] analytically and [0.070, 0.025] as frequencies – a factor sqrt((114 - 2) / (6 - 2)) = 5.29 – with p-values of 5e-04 against 6e-82. Fixed dispersion families (binomial, Poisson, negative-binomial) return identical coefficients AND standard errors either way, because phi is never estimated there; only residual_df and the BIC sample size move. Their likelihood is sum w log f under both conventions (a discrete density has no analytic-weight form), so ll_nobs is sum(w) for them regardless – see GLMResult.ll_nobs.

• max_iter (int) – IRLS controls (convergence on the relative deviance change).

• tol (float) – IRLS controls (convergence on the relative deviance change).

• robust (bool | str) –

  which coefficient-covariance estimator to report (recorded on the result as cov_type). False – the default – is the model-based phi (X'WX)^-1; True or "HC0" is the Huber–White sandwich; "HC1" | "HC2" | "HC3" apply the standard finite-sample leverage corrections to the same PSD-by-construction sandwich (MacKinnon–White 1985):

  • HC0: the plain sandwich. No finite-sample correction, so on small or leveraged designs it is biased LOW – measured 16-51% below HC3 on a 12-point design with one leveraged row, i.e. exactly where robust standard errors are reached for. (robust=True keeps meaning HC0 for compatibility.)

  • HC1: HC0 scaled by n / (n - rank) – the df correction Stata’s robust applies. Fixes the average small-sample bias but ignores WHERE the leverage sits.

  • HC2: per-observation scores scaled by 1 / sqrt(1 - h_i); exactly unbiased under homoskedasticity. Use when leverage is uneven but n is moderate.

  • HC3: scores scaled by 1 / (1 - h_i) (jackknife-like) – the standard recommendation for n <~ 50 or clearly leveraged designs (Long–Ervin 2000); the most conservative of the four.

  All four assume INDEPENDENT observations (none is cluster-robust) and still require the mean/link to be correct (audit G-7). h_i is the observation’s leverage in the converged working weighted regression; under weight_type="frequency" each replicate carries its own leverage, so every estimator equals the expanded data set’s. HC1-3 are refused on a saturated fit (their corrections divide by zero there), and HC2/HC3 on a unit-leverage observation (one fitted by itself alone).

• theta (float) – the negative-binomial dispersion parameter (family="negativebinomial" only; ignored otherwise). Not estimated from data – pass the value appropriate to your data (e.g. from a prior fit or a method-of-moments estimate); the default 1.0 is a plain placeholder, not a fitted value. WARNING (audit G-4): the MODEL-BASED standard errors assume theta is correct – with the placeholder against true theta = 10, measured SEs were 1.87x the true sampling spread (nominal 95% Wald intervals covering ~100%, tests losing essentially all power). Pass robust=True unless theta is trusted, and if theta came from the same data, remember AIC/BIC do not count it.

Returns:

A GLMResult.

Raises:

PerfectSeparationError – a binomial fit whose fitted probabilities reached exactly 0 or 1, leaving IRLS no finite step. Separation that stalls just short of that returns instead, with GLMResult.separated set, a UserWarning, and z_values / p_values refusing – so treat separated, not the absence of this exception, as the answer to “were my classes separated?” (the exception fires on a continuously coded separation and not on the dummy-coded form of the very same data).

Return type:

GLMResult

class GLMResult(coef, se, fitted, deviance, dispersion, log_likelihood, n_iter, family, link, cov, converged=True, rank=None, residual_df=None, dispersion_estimated=False, weight_type='analytic', ll_nobs=None, separated=False, cov_type='model', _link=None)

Bases: object

Fitted GLM.

Parameters:

• coef (ndarray)

• se (ndarray)

• fitted (ndarray)

• deviance (float)

• dispersion (float)

• log_likelihood (float | None)

• n_iter (int)

• family (str)

• link (str)

• cov (ndarray)

• converged (bool)

• rank (int | None)

• residual_df (int | None)

• dispersion_estimated (bool)

• weight_type (str)

• ll_nobs (float | None)

• separated (bool)

• cov_type (str)

• _link (Link)

coef

(p,) coefficient estimates.

Type:

numpy.ndarray

se

(p,) standard errors (model-based, or the robust variant named by cov_type). All NaN on a saturated dispersion-estimating fit – see dispersion.

Type:

numpy.ndarray

cov_type

which estimator se / cov came from: "model" for phi (X'WX)^-1, or "HC0" | "HC1" | "HC2" | "HC3" for the sandwich variants (see glm()’s robust argument). Recorded on the result for the same reason weight_type is: the meaning of a headline output must be readable off the result.

Type:

str

fitted

(n,) fitted means mu.

Type:

numpy.ndarray

deviance

residual deviance.

Type:

float

dispersion

estimated/assumed dispersion phi – the residual_df-corrected (Pearson) estimate the standard errors use, NOT the value the log-likelihood is evaluated at. NaN when a dispersion-estimating family is fit with residual_df <= 0 (a SATURATED fit – as many estimated parameters as observations): the coefficients, fitted values and deviance are still returned, but the dispersion, every standard error and the likelihood are undefined there, announced by a UserWarning at fit time, and z_values / p_values / aic / bic refuse with that reason.

Type:

float

residual_df

residual degrees of freedom, n - rank, where n counts rows under weight_type="analytic" and sum(weights) under weight_type="frequency". This is the divisor of dispersion and the t reference p_values uses, so it is where the weight convention shows up in the standard errors.

Type:

int | None

log_likelihood

maximised log-likelihood, or None when the supplied family/response does not define one – or when the fit is saturated (see dispersion), where the likelihood is unbounded in the dispersion and so has no maximised value. For dispersion-estimating families this is evaluated at the family’s own dispersion MLE (see _dispersion_mle), which differs from dispersion by design: the covariance wants the df-corrected estimate, the likelihood wants its maximiser. It is also computed under weight_type: analytic weights enter as per-observation precisions (Var = phi / w_i), frequency weights as replicate counts.

Type:

float | None

ll_nobs

how many observations log_likelihood is a sum over – the sample size bic penalises with. Row count for an analytic-weighted dispersion-estimating family; sum(weights) whenever the likelihood is a weighted sum of per-unit densities (every frequency-weighted fit, and every fixed-dispersion family, whose discrete densities have no analytic-weight form). Pairing a likelihood over sum(w) units with a log(rows) penalty is what made BIC agree with neither convention.

Type:

float | None

weight_type

the convention weights were read under; see glm().

Type:

str

separated

True when the binomial design separates the classes, so at least one coefficient has no finite maximum-likelihood estimate. deviance, fitted and log_likelihood are still correct (compare nested models by deviance), but coef and se along the separating direction record only where IRLS stopped, so z_values / p_values refuse. Announced by a UserWarning at fit time.

Type:

bool

n_iter

IRLS iterations to convergence.

Type:

int

converged

whether the IRLS convergence criterion was met. Public fits currently raise instead of returning this as false.

Type:

bool

rank

effective NUMERICAL rank of the weighted design sqrt(W) X at the solution, from the same SVD the covariance is computed from (cutoff at cond(X), not cond(X)^2; audit B4). rank < p – exact or near collinearity – is announced by a UserWarning at fit time and makes z_values / p_values refuse.

Type:

int | None

family / link

names.

cov

(p, p) coefficient covariance.

Type:

numpy.ndarray

Call forms: aic / bic are PROPERTIES (fit.aic, no parentheses) while z_values / p_values are METHODS (fit.z_values()). Both call forms are frozen for compatibility – changing either direction breaks every existing caller – so the split is stated here instead: mixing them up fails with TypeError: 'float' object is not callable (calling the property) or a bound-method object where an array was expected (reading the method), and neither traceback names the actual mistake. All four refuse with a stated reason (no likelihood, rank deficiency, separation, saturation) rather than return NaN.

predict(x, *, offset=None)

Predict the mean response mu at new design rows x.

Parameters:

• x (ndarray)

• offset (ndarray | None)

Return type:

ndarray

property aic: float

-2 ll + 2 k with k = rank + estimated dispersion.

Type:

Akaike information criterion

property bic: float

Bayesian information criterion over the observations the log-likelihood spans.

The penalty has to count the SAME observations the likelihood does (ll_nobs). Deriving it from residual_df instead pairs a likelihood summed over sum(weights) replicates with a log(rows) penalty on a frequency table, which is correct under neither weight convention – 107.7 where the two defensible answers were 116.5 and a likelihood of -3.5.

z_values()

Return Wald statistics for fitted coefficients (full-rank designs only).

Return type:

ndarray

p_values()

Two-sided Wald p-values: Student-t on the residual df when the dispersion was ESTIMATED (Gaussian/Gamma/inverse-Gaussian – the plug-in-dispersion normal reference rejected 9% at nominal 5% at n=8; audit G-1), normal otherwise (fixed-dispersion families), both asymptotic in the non-Gaussian mean model.

Return type:

ndarray

exception PerfectSeparationError

Bases: RuntimeError

The binomial classes are perfectly separated (or perfectly predicted): no finite MLE exists.

Raised by glm() when IRLS runs the fitted binomial probabilities all the way onto 0 or 1 and can no longer take a step – the terminal form of (quasi-)complete separation, where |coef| diverges and the working weights vanish. A subclass of RuntimeError, so pre-existing handlers keep working.

Separation that stops the iteration just SHORT of that point is the same statistical condition and is reported too, but as a UserWarning plus GLMResult.separated, because there the fit still carries a correct deviance / fitted / likelihood; only its Wald branch is undefined and refuses. Which of the two a given data set produces is a numerical accident, so never read the absence of this exception as evidence of no separation – read separated.

class Family(name, variance, canonical, unit_deviance, estimate_dispersion, extra=1.0, default_link=None)

Bases: object

An exponential-family error model: variance function, canonical link, deviance, dispersion.

canonical names the family’s mathematical canonical link; default_link (when set) is the link glm() fits with when none is requested. They are separate so families whose canonical link is numerically awkward (gamma’s inverse, inverse-Gaussian’s inverse_squared – neither keeps mu positive) can default to log without mislabeling the canonical link.

Parameters:

• name (str)

• variance (Callable[[ndarray], ndarray])

• canonical (str)

• unit_deviance (Callable[[ndarray, ndarray], ndarray])

• estimate_dispersion (bool)

• extra (float)

• default_link (str | None)

ridge_regression(x, y, alpha=1.0, *, fit_intercept=True)

Ridge (L2-penalised) linear regression in closed form.

Minimises ||y - X b||^2 + alpha ||b||^2; the intercept (if fitted) is not penalised.

Parameters:

• x (ndarray)

• y (ndarray)

• alpha (float)

• fit_intercept (bool)

Return type:

PenalizedResult

elastic_net(x, y, alpha=1.0, l1_ratio=0.5, *, fit_intercept=True, max_iter=1000, tol=1e-7)

Elastic-net linear regression by cyclic coordinate descent.

Minimises (1/2n) ||y - X b||^2 + alpha ( l1_ratio ||b||_1 + (1 - l1_ratio)/2 ||b||^2 ). l1_ratio = 1 is the lasso (sparse); l1_ratio = 0 is ridge-shaped, but note the 1/(2n) data term means elastic_net(X, y, a, 0.0) equals ridge_regression(X, y, n*a) – the penalties differ by a factor of n when cross-walking alpha (audit G-5).

Parameters:

• x (ndarray)

• y (ndarray)

• alpha (float)

• l1_ratio (float)

• fit_intercept (bool)

• max_iter (int)

• tol (float)

Return type:

PenalizedResult

lasso(x, y, alpha=1.0, **kw)

Lasso (L1) linear regression – elastic_net() with l1_ratio = 1.

Parameters:

• x (ndarray)

• y (ndarray)

• alpha (float)

Return type:

PenalizedResult

class PenalizedResult(coef, intercept, alpha, l1_ratio, n_iter, converged=True, rank=None)

Bases: object

Fitted penalized linear regression.

Parameters:

• coef (ndarray)

• intercept (float)

• alpha (float)

• l1_ratio (float)

• n_iter (int)

• converged (bool)

• rank (int | None)

coef

(p,) coefficients (excluding the intercept).

Type:

numpy.ndarray

intercept

fitted intercept.

Type:

float

alpha

overall penalty strength.

Type:

float

l1_ratio

elastic-net mixing (1 = lasso, 0 = ridge).

Type:

float

n_iter

coordinate-descent iterations (0 for the closed-form ridge).

Type:

int

predict(x)

Predict penalized-regression responses for design rows.

Parameters:

x (ndarray)

Return type:

ndarray

robust_regression(x, y, *, method='huber', c=None, max_iter=100, tol=1e-8)

Robust (M-estimator) linear regression by IRLS with a robust scale.

Down-weights observations with large residuals so a few outliers cannot dominate the fit. huber uses the Huber weight (tuning c = 1.345 for 95% Gaussian efficiency); tukey uses the redescending Tukey biweight (c = 4.685), which rejects gross outliers entirely.

The scale estimate (median absolute residual deviation / 0.6745) is only ever as robust as any MAD-based estimator is: past its own breakdown point – half or more of the rows sharing (near enough) the same residual under the current fit, an ordinary shape for zero-inflated or capped/rounded data and not a contrived one – it is exactly zero by construction and floors at 1e-8 (or, as of audit R-3, is treated the same way once it is merely far smaller than the response’s own spread, even without literally reaching that floor). That is an inherent limit of this estimator family, not something iterating further or switching method fixes: tukey’s redescending weight then hard-zeros essentially every row (raised below), and huber’s never-quite-zero weight instead settles – sometimes only after many iterations – on a near-zero-coefficient fit matching the majority tie that reports converged=True. A capped or censored measurement breaks down the same way without ever collapsing the scale estimate itself: a shared ceiling ties rows’ raw responses together but not their residuals once the fit carries any slope, so this is instead read directly off the response (see the module comment above _robust_weight_collapse in this file for the full two-signal detection design, including what audit R-3 closed and what it plainly did not: most notably, a majority and a minority that follow two genuinely DIFFERENT relationships – as opposed to sharing one tied value – is only ever PARTIALLY visible to this mechanism, and a location-only split with no covariate to anchor the response’s own spread is not visible to it at all). All of these are disclosed rather than left looking like an ordinary converged fit: see RegressionFit.degenerate_scale and the UserWarning issued at fit time. What is NOT, and cannot be, disclosed: WHICH of “a real minority signal was discarded,” “a minority of contamination was correctly rejected,” or “nothing was discarded at all, because the majority is genuinely this exact” produced a flagged pattern – all three leave the same fingerprint, so the warning names the observed condition, never a verdict on which case it is.

A sharper, method='huber'-specific gap (GAP C in the module comment, campaign nine, D-0209): huber’s weight never fully crushes, so a real minority following a different relationship can be systematically down-weighted enough to bias coef toward the majority WITHOUT tripping either disclosure signal at all – confirmed common (up to 100% of seeds at a 35% minority fraction, at everyday noise levels), not a rare edge case. Raising max_iter past a RuntimeError("... failed to converge ...") is what usually exposes it: that default failure is a real but incomplete guard, not a guarantee, so a non-convergence error under method='huber' is worth treating as a possible sign of exactly this shape rather than simply raising max_iter and moving on. Switching to method='tukey' is not a reliable workaround for this specific shape either outside a no-intercept design.

Raises:

RuntimeError – every observation lands at zero weight (the message names the breakdown pattern above when one was detected), coefficients become non-finite, or IRLS does not converge in max_iter iterations.

Parameters:

• x (ndarray)

• y (ndarray)

• method (str)

• c (float | None)

• max_iter (int)

• tol (float)

Return type:

RegressionFit

quantile_regression(x, y, tau=0.5, *, max_iter=200, tol=1e-7, eps=1e-6)

Linear quantile regression: the conditional tau-quantile by IRLS on the check loss.

Minimises the pinball loss sum rho_tau(y - X b) via iteratively reweighted least squares with weights tau / |r| for positive residuals and (1 - tau) / |r| for negative ones (a smoothed Newton scheme; eps floors |r| for stability).

Parameters:

• x (ndarray)

• y (ndarray)

• tau (float)

• max_iter (int)

• tol (float)

• eps (float)

Return type:

RegressionFit

class RegressionFit(coef, fitted, scale, n_iter, converged=True, rank=None, degenerate_scale=False)

Bases: object

Coefficients + fitted values from robust_regression() / quantile_regression().

Parameters:

• coef (ndarray)

• fitted (ndarray)

• scale (float)

• n_iter (int)

• converged (bool)

• rank (int | None)

• degenerate_scale (bool)

degenerate_scale

set by robust_regression() when it detects a sign of its own >50% breakdown point – EITHER the IRLS scale has collapsed (at or within a convergence step of its numerical floor, or far below the raw response’s own spread) with a real share of rows assigned (near) zero weight as a result, OR at least ~44% of the raw response is tied (near enough) to a single value regardless of the fit’s residuals (the shape a capped/censored/zero-inflated measurement produces, which need not collapse the scale at all; see the module comment above _robust_weight_collapse in this file for the full design and, just as importantly, what it still cannot see). Either signal firing means this is not a transient numerical accident, and coef MAY reflect a majority pattern rather than a fully recovered relationship – but may equally reflect a fit that is completely correct, with the flagged pattern arising because the majority genuinely fits (near-)exactly; see robust_regression(). This field, and the UserWarning announcing it at fit time, name the observed CONDITION only, never a verdict: neither can establish whether a real minority signal was discarded or a genuine minority of contamination was correctly rejected (both look identical from this pattern alone), and – a further, structural limitation, not merely an unresolved ambiguity – a majority and minority that follow two genuinely different relationships (as opposed to sharing one tied value) is largely invisible to this mechanism in the first place unless the majority’s own fit is exceptionally tight; see the module comment for exactly how far that reaches and where it stops. Always False for quantile_regression(), which has no comparable mechanism.

Type:

bool

predict(x)

Predict fitted regression values for design rows.

Parameters:

x (ndarray)

Return type:

ndarray

kaplan_meier(time, event=None, *, ci_level=0.95)

Kaplan–Meier product-limit estimate of the survival function S(t).

Parameters:

• time (ndarray) – (n,) observed times (event or censoring).

• event (ndarray | None) – (n,) 1 = event, 0 = right-censored (defaults to all events).

• ci_level (float) – confidence level for the POINTWISE log–log interval at each time point.

Returns:

{'time', 'survival', 'se', 'ci_low', 'ci_high', 'at_risk', 'n_events', 'median'}.

Return type:

dict[str, ndarray]

ci_low/ci_high are PER-TIME-POINT (pointwise) intervals, not a simultaneous confidence band: each is a MARGINAL, asymptotic (Greenwood log–log) interval for S(t) at its own t with probability ci_level under independent censoring, but the probability that the true curve stays inside the envelope at EVERY time point is far lower – measured simultaneous coverage of the pointwise 95% envelope is 0.60 / 0.51 / 0.43 at n = 50 / 100 / 200 (it FALLS with n, because more event times mean more chances to escape). Reading the plotted envelope as “the curve lives in here” is a simultaneous claim these intervals do not support; a Hall–Wellner band is the standard tool for that and is not implemented here.

nelson_aalen(time, event=None)

Nelson–Aalen estimate of the cumulative hazard H(t) = sum d_i / Y_i.

Returns:

{'time', 'cumhaz', 'se'} with the Poisson-type standard error of the cumulative hazard.

Parameters:

• time (ndarray)

• event (ndarray | None)

Return type:

dict[str, ndarray]

cox_ph(x, time, event, *, start=None, strata=None, subject=None, ties='efron', max_iter=50, tol=1e-9)

Cox proportional-hazards regression by Newton–Raphson on the partial likelihood.

The hazard is h(t | x) = h0(t) exp(x' beta); only the ordering of event times enters, so the baseline h0 is left unspecified (semi-parametric). Time-varying covariates are supported through the counting-process form: pass start so each row is an at-risk interval (start, stop] (a subject contributes several rows), and the risk set at an event time is every interval covering it.

Parameters:

• x (ndarray) – (n, p) covariates (no intercept – it is absorbed into the baseline).

• time (ndarray) – (n,) event/censoring times (the interval stop times).

• event (ndarray) – (n,) 1 = event, 0 = censored.

• start (ndarray | None) – optional (n,) interval start times for time-varying covariates / left truncation.

• strata (ndarray | None) – optional (n,) labels; each stratum gets its own baseline hazard (coefficients shared).

• subject (ndarray | None) – optional (n,) subject identifiers. Required to compute concordance for counting-process input with multiple rows per subject.

• ties (str) – "efron" (default, more accurate) or "breslow" tie handling.

• max_iter (int) – Newton controls.

• tol (float) – Newton controls.

Returns:

A CoxResult.

Return type:

CoxResult

What the standard errors assume – read this before using counting-process input for clustered or recurrent events. The SEs are the inverse observed information of the partial likelihood, which is valid when SUBJECTS are independent. The (start, stop] interface happily accepts multiple rows per subject (that is its purpose: time-varying covariates for one subject) and also LOOKS like it accepts recurrent events or cluster members as rows – but rows from the same subject/cluster are dependent, the information-based SEs do not know that, and they understate: on a simulated clustered design the mean reported SE was 0.68 of the true sampling SD. For clustered data use frailty_cox() (models the dependence) or bootstrap over SUBJECTS/CLUSTERS – resampling rows would repeat the same mistake. Proportional hazards itself (covariate effects constant over time) is assumed, not checked: no diagnostic ships in this module, so inspect Schoenfeld-style residuals or stratified log(-log) survival plots externally before reading coef as a time-constant log-hazard ratio.

class CoxResult(coef, se, cov, loglik, baseline_time, baseline_cumhaz, concordance, n_iter, converged, baseline_by_stratum=<factory>)

Bases: object

Fitted Cox proportional-hazards model.

Parameters:

• coef (ndarray)

• se (ndarray)

• cov (ndarray)

• loglik (float)

• baseline_time (ndarray)

• baseline_cumhaz (ndarray)

• concordance (float | None)

• n_iter (int)

• converged (bool)

• baseline_by_stratum (dict[object, dict[str, ndarray]])

coef

(p,) log-hazard-ratio coefficients.

Type:

numpy.ndarray

se

(p,) standard errors (inverse observed information).

Type:

numpy.ndarray

cov

(p, p) covariance.

Type:

numpy.ndarray

loglik

maximised partial log-likelihood.

Type:

float

baseline_time / baseline_cumhaz

Breslow baseline cumulative hazard for an unstratified fit. These arrays are empty for a multi-stratum fit.

baseline_by_stratum

mapping from each stratum label to its own time and cumhaz arrays.

Type:

dict[object, dict[str, numpy.ndarray]]

concordance

Harrell’s C-index, evaluated only within strata. This is None for counting-process input unless subject identifiers are supplied.

Type:

float | None

n_iter

Newton iterations.

Type:

int

converged

whether Newton’s method reached tol with a finite, well-identified fit. False means either it hit max_iter without reaching tol, a step produced a non-finite iterate (kept the last finite one instead), or the final standard errors are enormous – the classic symptom of (quasi-)complete separation, where the partial likelihood has no finite maximizer and Newton’s steps can shrink below tol while the coefficients themselves keep drifting toward infinity.

Type:

bool

hazard_ratios()

Return exponentiated Cox coefficients.

Return type:

ndarray

z_values()

Return Wald z statistics for Cox coefficients.

Return type:

ndarray

p_values()

Return two-sided normal-approximation p-values for Cox coefficients.

Return type:

ndarray

to_person_period(time, event, covariates=None)

Expand right-censored durations into a person-period (long) array for discrete-time models.

Each subject contributes one row per discrete period they were at risk; the binary outcome is 1 in the period the event occurred and 0 otherwise. Integer time is the number of periods observed.

Returns:

{'period', 'outcome', 'subject', 'covariates'} (covariates repeated per period if given).

Parameters:

• time (ndarray)

• event (ndarray)

• covariates (ndarray | None)

Return type:

dict[str, ndarray]

discrete_time_hazard(x, outcome, *, link='cloglog', offset=None)

Discrete-time hazard model: a binary GLM on the person-period array.

Fit on the long-format data from to_person_period() (the design x typically holds period indicators / a time trend plus covariates). cloglog gives the grouped-proportional-hazards (interval-censored Cox) interpretation; logit gives the proportional-odds hazard.

Returns:

a mixle.inference.glm.GLMResult (binomial family with the chosen link).

Parameters:

• x (ndarray)

• outcome (ndarray)

• link (str)

• offset (ndarray | None)

aalen_johansen(time, event, *, causes=None)

Aalen–Johansen cumulative incidence functions for competing risks.

With several mutually exclusive event types, the cause-specific CIF F_k(t) is the probability of failing from cause k by time t accounting for the competing causes (it is not 1 - KM on the cause, which overstates incidence).

Parameters:

• time (ndarray) – (n,) event/censoring times.

• event (ndarray) – (n,) integer cause label, 0 = censored, 1..K = causes.

• causes (ndarray | None) – optional explicit list of cause labels; inferred from event if None.

Returns:

{'time', 'cif': {cause: array}, 'overall_survival'}.

Return type:

dict

aalen_additive(x, time, event, *, intercept=True)

Aalen’s additive-hazards regression: cumulative regression functions B(t).

Models h(t | x) = b0(t) + sum_j x_j b_j(t) with time-varying additive effects. At each event time the increment dB is the least-squares solution over the risk set; the cumulative B(t) (returned) has interpretable slopes – a rising B_j means covariate j adds hazard.

The estimator is only defined while the risk-set design has full column rank, and it degrades BEFORE it breaks: as the risk set shrinks toward p subjects the least-squares increments are dominated by noise, and the accumulated tail of B(t) is numerical garbage that LOOKS like a dramatic late effect (measured: max abs(B) growing 1.4 -> 19 over the last few event times of a null design). Estimation therefore STOPS at the first event time whose risk-set design is rank-deficient or numerically ill-conditioned – the standard Aalen practice – instead of silently substituting a minimum-norm solution and accumulating it. truncated_at reports where (None if the full curve was estimable); the curve up to that point is unaffected by the truncation.

Returns:

{'time', 'cum_coef', 'truncated_at'} where cum_coef is (n_estimated_times, p[+1]) cumulative coefficients (the first column is the baseline when intercept is True) and truncated_at is the first event time excluded because the risk-set design lost full rank (None when no truncation occurred).

Parameters:

• x (ndarray)

• time (ndarray)

• event (ndarray)

• intercept (bool)

Return type:

dict

frailty_cox(x, time, event, groups, *, max_iter=50, tol=1e-5, ties='breslow', se_method='complete-data')

Shared gamma-frailty Cox model for clustered survival, by EM.

Subjects in the same group share an unobserved frailty w_g ~ Gamma(1/theta, 1/theta) (mean 1, variance theta) that multiplies the hazard, capturing within-group correlation. The E-step retains E[w_g], Var(w_g), and E[log(w_g)] under the conjugate gamma posterior. The coefficient/baseline M-step uses E[w_g] in the integrated-hazard term, while the dispersion M-step maximises the expected gamma-prior log likelihood using both E[w_g] and E[log(w_g)].

On theta: it is returned as a point estimate ONLY – see the honesty notes on FrailtyCoxResult before reading a small value as evidence for (or against) clustering.

Parameters:

• se_method (str) – how the fixed-effect SEs are computed. "complete-data" (default) inverts the partial-likelihood information with the posterior-mean frailties and theta PLUGGED IN as if known – it prices in none of their estimation uncertainty and understates: measured against the true sampling SD at theta = 1, se/SD was 0.79-0.88 and Wald-95 coverage 0.87-0.91 (clean at theta = 0). "jackknife" deletes one GROUP at a time and refits, pricing in frailty and theta estimation; it is honest but costs G extra EM fits, and its own validity is asymptotic in the NUMBER OF GROUPS (a handful of groups gives a noisy SE from few leave-outs).

• x (ndarray)

• time (ndarray)

• event (ndarray)

• groups (ndarray)

• max_iter (int)

• tol (float)

• ties (str)

Returns:

A FrailtyCoxResult (se_method records which SEs you got).

Return type:

FrailtyCoxResult

class FrailtyCoxResult(coef, se, theta, frailties, groups, frailty_variance=<factory>, frailty_log_mean=<factory>, n_iter=0, converged=False, ties='breslow', se_method='complete-data')

Bases: object

Shared gamma-frailty Cox result.

Parameters:

• coef (ndarray)

• se (ndarray)

• theta (float)

• frailties (ndarray)

• groups (ndarray)

• frailty_variance (ndarray)

• frailty_log_mean (ndarray)

• n_iter (int)

• converged (bool)

• ties (str)

• se_method (str)

coef / se

fixed-effect log-hazard-ratios and standard errors (see se_method for what the SEs account for – the default plug-in SEs understate).

theta

estimated frailty variance. A POINT estimate with no uncertainty attached: no SE or test accompanies it, the null theta = 0 sits on the parameter boundary (a likelihood ratio there is a 50:50 chi2_0/chi2_1 mixture, not implemented), and the EM estimate never actually reaches 0 under a true theta of 0 (measured mean 0.046 on null simulations) – so a small positive theta is what NO clustering looks like, not weak evidence of some.

Type:

float

frailties

posterior mean random effect per group.

Type:

numpy.ndarray

frailty_variance

posterior variance of each group random effect.

Type:

numpy.ndarray

frailty_log_mean

posterior expectation of log(w_g).

Type:

numpy.ndarray

groups

group labels aligned to frailties.

Type:

numpy.ndarray

n_iter

EM iterations.

Type:

int

converged

whether the EM convergence criterion was met.

Type:

bool

se_method

"complete-data" (plug-in information, understates) or "jackknife" (delete-one-group, accounts for frailty/theta estimation).

Type:

str

ordinal_regression(x, y, *, link='logit', max_iter=200)

Fit a cumulative-link ordinal regression (ordered logit / probit) by maximum likelihood.

Parameters:

• x (ndarray) – (n, p) covariates (no intercept – the thresholds play that role).

• y (ndarray) – (n,) integer category labels 0..K-1 (ordered).

• link (str) – "logit" (proportional odds) or "probit".

• max_iter (int) – optimiser iterations.

Returns:

An OrdinalResult.

Return type:

OrdinalResult

class OrdinalResult(coef, thresholds, se, log_likelihood, link, n_categories, n_iter=0, converged=True, rank=None)

Bases: object

Fitted ordinal (cumulative-link) regression.

Parameters:

• coef (ndarray)

• thresholds (ndarray)

• se (ndarray)

• log_likelihood (float)

• link (str)

• n_categories (int)

• n_iter (int)

• converged (bool)

• rank (int | None)

coef

(p,) slope coefficients (positive beta_j raises the latent score, shifting mass toward higher categories).

Type:

numpy.ndarray

thresholds

(K-1,) ordered cut points alpha.

Type:

numpy.ndarray

se

(p,) standard errors for coef.

Type:

numpy.ndarray

log_likelihood

maximised log-likelihood.

Type:

float

link

"logit" or "probit".

Type:

str

n_categories

number of ordered categories K.

Type:

int

predict_proba(x)

Per-category probabilities (n, K) at design rows x.

Parameters:

x (ndarray)

Return type:

ndarray

predict(x)

Most-probable ordered category per row.

Parameters:

x (ndarray)

Return type:

ndarray

concordance_summary(x, y)

All pairwise concordance measures between two ordinal variables.

Returns:

{'kendall_tau_b', 'gamma', 'somers_d_yx', 'somers_d_xy', 'concordant', 'discordant', 'tx', 'ty', 'txy'} – tx = pairs tied on x only, ty on y only, txy on both.

Parameters:

• x (ndarray)

• y (ndarray)

Return type:

dict[str, float]

kendall_tau(x, y)

Kendall’s tau-b rank correlation (tie-corrected) between two ordinal variables.

Parameters:

• x (ndarray)

• y (ndarray)

Return type:

float

goodman_kruskal_gamma(x, y)

Goodman–Kruskal gamma: (C - D) / (C + D) ignoring ties.

Parameters:

• x (ndarray)

• y (ndarray)

Return type:

float

somers_d(x, y, *, dependent='y')

Somers’ D, the asymmetric rank association treating dependent as the response.

Parameters:

• x (ndarray)

• y (ndarray)

• dependent (str)

Return type:

float

mann_whitney_u(x, y, *, alternative='two-sided', use_continuity=True)

Mann-Whitney U / Wilcoxon rank-sum test for two independent samples.

THE NULL IS FULL EXCHANGEABILITY (F = G), not mere stochastic equality: the reference variance is derived under identical distributions, so with equal P(X > Y) = 1/2 but UNEQUAL shapes/spreads and unbalanced sample sizes the level is not controlled – measured 12.5% rejection at nominal 5% one way and 1.5% with the sizes reversed (STAT-RR17-17). For the stochastic-equality null under unequal shapes use brunner_munzel(), whose variance estimates each group separately. Uses mid-ranks for ties, the tie-corrected normal approximation, and (default) a continuity correction. alternative is 'two-sided', 'greater' (x > y), or 'less'. The rank-biserial correlation 2*U1/(n1 n2) - 1 is reported as the effect size; null_hypothesis on the result names the null in force.

Parameters:

• x (Any)

• y (Any)

• alternative (str)

• use_continuity (bool)

Return type:

MannWhitneyResult

class MannWhitneyResult(statistic, statistic2, zscore, pvalue, rank_biserial, alternative, null_hypothesis='exchangeability (F = G); use brunner_munzel for unequal shapes')

Bases: object

Result of a two-sample Mann-Whitney U test.

Parameters:

• statistic (float)

• statistic2 (float)

• zscore (float)

• pvalue (float)

• rank_biserial (float)

• alternative (str)

• null_hypothesis (str)

wilcoxon_signed_rank(x, y=None, *, alternative='two-sided', zero_method='wilcox', correction=False)

Wilcoxon signed-rank test for paired samples (or one sample vs 0).

Ranks |d| for d = x - y (mid-ranks for ties), splits into positive / negative rank sums, and uses the EXACT subset-sum null when it is available – n <= 300 with no zero differences and no tied |d| – and the tie-corrected normal approximation otherwise (the same regime switch SciPy and R make; the normal approximation is level-violating at very small n, where its smallest attainable p sits below the exact one). zero_method='wilcox' drops zero differences (and their ranks); 'pratt' keeps them in the ranking but drops them from the sums, with the matching Pratt/Cureton zero corrections applied to the null mean and variance (as scipy does). correction applies only to the normal branch (exact tails need no continuity repair). The matched-pairs rank-biserial correlation is reported as the effect size.

Parameters:

• x (Any)

• y (Any)

• alternative (str)

• zero_method (str)

• correction (bool)

Return type:

WilcoxonResult

class WilcoxonResult(statistic, zscore, pvalue, rank_biserial, alternative, method='normal')

Bases: object

Result of a paired or one-sample Wilcoxon signed-rank test.

method names the null actually used (STAT-RR17-15): "exact" (full enumeration, the n <= 300 no-ties/no-zeros regime) or "normal" (tie/zero-corrected approximation). zscore is signed FOR THE STATED ALTERNATIVE – positive means the data lean toward it – so an all-positive sample under alternative='greater' reports a positive z next to its small exact p, where the old min(W+, W-)-based z reported -2.02 beside p = 0.03125 and a rank-biserial of +1.

Parameters:

• statistic (float)

• zscore (float)

• pvalue (float)

• rank_biserial (float)

• alternative (str)

• method (str)

sign_test(x, y=None, *, alternative='two-sided')

Sign test for paired samples (or one sample vs 0): exact binomial test on the signs of x - y.

Only the directions of the differences are used (ties dropped), so it is maximally robust but less powerful than the signed-rank test. extra carries n_positive and n (non-zero pairs).

Parameters:

• x (Any)

• y (Any)

• alternative (str)

Return type:

TestResult

kruskal_wallis(*samples)

Kruskal-Wallis H test: the rank-based k-sample generalization of Mann-Whitney (one-way ANOVA).

Tie-corrected H with a chi-square(k-1) reference. extra carries df and the epsilon_squared effect size (H - k + 1)/(N - k).

Parameters:

samples (Any)

Return type:

TestResult

friedman_test(*measurements)

Friedman test for k related samples (repeated measures): the rank-based repeated-measures ANOVA.

Pass each treatment as a separate equal-length array (one value per block). Ranks within each block, tie-corrects, and uses a chi-square(k-1) reference. extra carries df and Kendall’s W concordance effect size.

Parameters:

measurements (Any)

Return type:

TestResult

brunner_munzel(x, y, *, alternative='two-sided', distribution='t')

Brunner-Munzel test: the generalized Wilcoxon test that does not assume equal variances/shapes.

Tests the stochastic-equality null P(x < y) + 0.5 P(x = y) = 1/2. alternative is 'two-sided', 'greater' (x > y), or 'less' – the same direction convention as mann_whitney_u() (and scipy). distribution='t' uses a Satterthwaite t reference, 'normal' the normal approximation. The t reference is the better of the two but is still LIBERAL in small samples – below roughly 10 observations per group its true level runs above nominal (this is documented behavior of the approximation, not an implementation defect), so treat sub-10-per-group p-values near the threshold as optimistic; the studentized-permutation version of the test is the small-sample fix and is not implemented here. Reports the estimated relative effect p_hat = P(x < y) + 0.5 P(x = y) in extra.

Completely separated samples (all of one group beyond all of the other) make the variance estimate 0 and the studentized statistic undefined – and NO valid p-value exists under this test’s own null. The permutation tail 1/C(n1+n2, n1) an earlier revision reported is exact only under FULL exchangeability (F = G), a strictly stronger hypothesis than Brunner-Munzel’s stochastic equality: with X = +/-1 equiprobable and Y = 0, P(X < Y) = 0.5 holds exactly, yet all-X-below-Y has probability 2^(-n1) – at n1 = 4 the old report rejected 12.5% of that null at nominal 5% while quoting p = 0.0002 (STAT-RR19-04). The result therefore reports pvalue = nan, statistic = +/-inf, extra['method'] = 'separation-no-valid-p-under-stochastic-equality', the observed p_hat, and – for callers prepared to ASSERT the stronger F = G null – extra['p_exchangeability'], the exact permutation bound labeled as exactly that. Use mann_whitney_u() or a permutation test when exchangeability is the hypothesis you mean.

Parameters:

• x (Any)

• y (Any)

• alternative (str)

• distribution (str)

Return type:

TestResult

mood_median_test(*samples, ties='below')

Mood’s median test: contingency test that k samples sit alike around the POOLED median.

Cross-tabulates each observation as above / (at-or-below) the pooled grand median and tests independence of the resulting 2xk table. extra carries the grand_median, the reference used (method), and min_expected_count.

What rejection means – precisely: the groups differ in their probability of exceeding the pooled median. That is implied by unequal medians under a shift model, but it is NOT “the medians differ” unconditionally: groups with identical medians and different shapes near that median can reject, and the test has famously low power besides. Reference: the median split fixes BOTH margins of the table, so for k = 2 the conditional exact (Fisher) p-value is the calibrated one and is used automatically whenever any expected cell count is below 5 (the uncorrected chi-square over-rejects there); for k > 2 the chi-square approximation is all that is implemented – check min_expected_count and treat p-values from tables with expected counts below ~5 as approximate.

Parameters:

• samples (Any)

• ties (str)

Return type:

TestResult

dunn_test(*samples, p_adjust='holm')

Dunn’s post-hoc test: all pairwise rank-mean comparisons after a Kruskal-Wallis rejection.

Uses the pooled-rank z statistic with the shared tie-corrected variance, and adjusts the pairwise p-values by 'holm', 'bonferroni', or 'none'.

Each pairwise z is calibrated under the GLOBAL null (all k groups exchangeable), because both the pooled ranks and the shared variance come from the full pooled sample. That is exactly the null a post-hoc test runs under after rejecting it, so read the p-values as screening under partial nulls, not as k(k-1)/2 self-contained two-sample tests: when one group genuinely differs, its observations shift the pooled ranks and the shared variance for every OTHER pair too – with two identical groups plus one concentrated third, the identical pair’s nominal-5% rejection rate is ~24%. For a confirmatory conclusion about a specific pair, test that pair directly (brunner_munzel() uses only the two groups involved) with your own multiplicity adjustment; use Dunn to rank candidates, not to certify them.

Parameters:

• samples (Any)

• p_adjust (str)

Return type:

DunnResult

class DunnResult(comparisons, zscores, pvalues, p_adjust)

Bases: object

Post-hoc Dunn pairwise comparisons after Kruskal-Wallis.

Parameters:

• comparisons (list[tuple[int, int]])

• zscores (ndarray)

• pvalues (ndarray)

• p_adjust (str)

jonckheere_terpstra(*samples, alternative='increasing')

Jonckheere-Terpstra test for an ORDERED alternative across independent samples.

More powerful than Kruskal-Wallis when the groups are expected to shift monotonically in the given order. alternative='increasing' / 'decreasing' / 'two-sided'. Uses the tie-corrected normal approximation of the J statistic (sum of pairwise Mann-Whitney counts over ordered pairs).

The group ORDER is part of the hypothesis and must be chosen BEFORE looking at the data (the same pre-specification the sibling page_trend_test() states): the level calculation assumes one fixed ordering, and picking the ordering that best matches the observed group means – k!/2 candidate orderings – turns a nominal 5% test into a data-dredged one. If the ordering came from the data, this test’s p-value does not mean what it says.

The null variance is the published Jonckheere-Terpstra one (Lehmann 1975), whose leading term is [n(n-1)(2n+5) - sum n_i(n_i-1)(2n_i+5)] / 72 – algebraically the same as the equivalent [n^2(2n+3) - sum n_i^2(2n_i+3)] / 72 form. This implementation previously used n(n-1)(2n+3) - sum n_i(n_i-1)(2n_i+3), mixing the two (MXR-080-1599): a hybrid that matches neither published expression and, critically, does NOT reduce to the Mann-Whitney variance for two ordered groups. For two groups of five with complete separation it understated the variance as 21.5278 against the required n1*n2*(n+1)/12 = 22.9167, reporting z=2.69408, p=0.003529 where the rank-sum test gives z=2.61116, p=0.004512. The tied branch inherited the same altered base terms, including in its sum t_j(t_j-1)(2t_j+5) tie correction.

Reproducing that cross-check: the quoted rank-sum numbers use the normal approximation WITHOUT a continuity correction, because this J statistic’s normal reference carries none – call mann_whitney_u() with use_continuity=False to match (its DEFAULT is True, which reports a different, corrected p on the same data). The difference between the defaults is the correction convention, not a defect in either test.

Parameters:

• samples (Any)

• alternative (str)

Return type:

TestResult

page_trend_test(*measurements, decreasing=False)

Page’s trend test for an ORDERED alternative in repeated measures.

Like Friedman but for a pre-specified ordering of the k treatments (the columns, in order). Tests L = sum_j j * R_j against the normal approximation. Set decreasing=True to predict the reverse ordering. extra carries the z-score.

Ties are handled exactly, not ignored: within-block midranks enter L, and the null variance is the exact permutation variance of the observed (possibly tied) midranks – per block, Var = sum_j (w_j - w-bar)^2 * sum_i (r_i - r-bar)^2 / (k - 1) under uniform within-block permutation, summed over blocks. With no ties this reduces algebraically to the textbook n k^2 (k+1)(k^2-1) / 144; with ties the textbook constant overstates the variance (tied midranks vary less), which made the test conservative and its z biased toward 0.

Parameters:

• measurements (Any)

• decreasing (bool)

Return type:

TestResult

ks_1samp(x, cdf, *, alternative='two-sided')

One-sample Kolmogorov-Smirnov goodness-of-fit test against a fully-specified cdf callable.

cdf may be vectorized (ndarray -> ndarray) or scalar-only, such as the .cdf method of a fitted univariate mixle distribution; a scalar-only callable is evaluated element-wise. Mixtures, categoricals, composites, and sequence models define no cdf: for those, pass the CDF as your own callable (for a univariate mixture, the weighted sum of the component CDFs).

The p-value is ASYMPTOTIC in every mode: two-sided uses the limiting Kolmogorov distribution (scipy’s method='asymp') and the one-sided modes use the large-sample exponential bound exp(-2 n D^2). No small-sample (Marsaglia-Tsang-Wang exact) correction is applied, so at small n the reported p-value is LARGER than the exact one – it understates the evidence against the hypothesized cdf and overstates goodness-of-fit (e.g. ~26% relative near p=0.005 at n=50; the gap shrinks as n grows). For an exact small-sample p-value use scipy.stats.ks_1samp(x, cdf, method='exact'); the statistic here matches scipy’s to machine precision either way.

Parameters:

• x (Any)

• cdf (Callable[[ndarray], ndarray])

• alternative (str)

Return type:

TestResult

ks_2samp(x, y, *, alternative='two-sided')

Two-sample Kolmogorov-Smirnov test: max gap between the two empirical CDFs.

The p-value is exact at small samples and asymptotic otherwise (scipy method='auto'). The null distribution is distribution-free for CONTINUOUS data; with ties (discrete or rounded data) the p-value is conservative.

Parameters:

• x (Any)

• y (Any)

• alternative (str)

Return type:

TestResult

runs_test(x, *, cutoff='median')

Wald-Wolfowitz runs test for randomness of a binary/dichotomized sequence.

Dichotomizes x about its median (or a supplied numeric cutoff) and tests whether the run count departs from what independence predicts (too few runs => clustering/trend; too many => over-alternation). Two-sided. The null is EXACT (the closed-form run-count distribution given (n1, n2), big-integer arithmetic) for all n1 + n2 <= 5000 – the uncorrected normal approximation is level-violating at small and sparse sizes (exhaustive enumeration at n1 = n2 = 5 rejects 20/252 = 7.94%, and at n = 61, n1 = 3 9.83%, at nominal 5%; STAT-RR17-16/STAT-RR19-14). The normal approximation applies only above 5000 (STAT-RR22-16: this sentence previously still said 60 while execution switched at 5001). extra carries the run count, z-score, and method ("exact" / "normal").

Parameters:

• x (Any)

• cutoff (str | float)

Return type:

TestResult

cliffs_delta(x, y)

Cliff’s delta effect size in [-1, 1]: P(x > y) - P(x < y) (rank-based, ties count as 0).

Parameters:

• x (Any)

• y (Any)

Return type:

float

class TestResult(statistic, pvalue, extra=<factory>)

Bases: object

Generic statistic + p-value result; extra carries test-specific fields (effect size, df, …).

Parameters:

• statistic (float)

• pvalue (float)

• extra (dict[str, Any])

deming_regression(x, y, variance_ratio=1.0)

Errors-in-variables (Deming) regression of y on x when both are noisy.

Parameters:

• x – paired measurements; both may carry error.

• y – paired measurements; both may carry error.

• variance_ratio (float) – var(e_y) / var(e_x) – the ratio of output to input noise variance. 1.0 is total least squares (orthogonal regression); a large value -> ordinary least squares (no input error); a small value -> inverse regression (predictor dominated by error).

Returns:

A DemingFit with the unbiased slope / intercept and the recovered latent x*.

Return type:

DemingFit

class DemingFit(slope, intercept, variance_ratio, x, y)

Bases: object

Result of deming_regression(): slope/intercept plus the recovered latent predictor values.

conditional_mean(x_star)

The conditional mean E[y | x*] = a + b x* at true predictor values x*.

Parameters:

x_star (ndarray)

Return type:

ndarray

simex(fit_fn, x, y, sigma_u, *, lambdas=None, n_sims=100, extrapolation='quadratic', seed=0)

SIMEX: simulation–extrapolation correction for a predictor measured with known error.

When a predictor x is observed as x = x* + u with u ~ N(0, sigma_u^2), naive estimates are biased (attenuation). SIMEX adds further noise of variance lambda sigma_u^2 for a grid of lambda >= 0, refits at each level (averaging over n_sims noise draws), then extrapolates the estimate back to lambda = -1 (zero measurement error). Works for any estimator returning a parameter vector.

Parameters:

• fit_fn (Callable[[ndarray, ndarray], ndarray]) – f(x, y) -> theta returning the parameter vector for (possibly multi-column) x.

• x (ndarray) – (n,) or (n, p) error-prone predictor(s).

• y (ndarray) – (n,) response.

• sigma_u (float | ndarray) – measurement-error standard deviation (scalar, or per-column for matrix x).

• lambdas (ndarray | None) – extra-noise levels (defaults to 0, 0.5, 1.0, 1.5, 2.0).

• n_sims (int) – noise replications per level.

• extrapolation (str) – "quadratic" or "linear" extrapolant in lambda.

• seed (int | RandomState | None) – RNG seed.

Returns:

{'estimate', 'naive', 'lambdas', 'curve'} – the SIMEX-corrected parameter vector, the naive lambda=0 estimate, and the per-level averaged estimates.

Return type:

dict

propagate_uncertainty(func, samples, *, quantiles=(0.025, 0.5, 0.975))

Monte-Carlo propagation of an uncertainty set through an arbitrary functional.

Pushes input draws (a posterior sample, a bootstrap set, any uncertainty representation) through func and summarises the output distribution – the general “what is the uncertainty of g(theta)?” operation. func may be vectorised (accept the whole (n, ...) array) or act on a single draw; both are detected.

Parameters:

• func (Callable[[ndarray], ndarray]) – the functional to propagate. Returns a scalar or fixed-length vector per input draw.

• samples (ndarray) – (n, ...) input draws (rows are draws).

• quantiles (tuple[float, ...]) – output quantiles to report.

Returns:

{'mean', 'std', 'quantiles', 'levels', 'samples'} over the propagated outputs.

Return type:

dict

paired_score_difference(scores_a, scores_b, *, lower_is_better=True, ci_level=0.95)

Mean paired held-out score difference with a CI and a paired t-test.

Parameters:

• scores_a (ndarray) – (n,) per-observation held-out scores for the two models (same observations, same order).

• scores_b (ndarray) – (n,) per-observation held-out scores for the two models (same observations, same order).

• lower_is_better (bool) – True for losses/scores (CRPS, log loss, pinball); False for higher-is-better metrics.

• ci_level (float) – confidence level for the interval on the mean difference.

Returns:

{'mean_diff', 'se', 'ci_low', 'ci_high', 't', 'p_value', 'favored'} where mean_diff is mean(a - b) and favored is 'A' / 'B' / 'tie' at the given level. A zero- variance paired difference (every observation agrees exactly) is maximally significant when the shared value is nonzero, not an automatic tie – see the se == 0 handling below.

Return type:

dict

vuong_test(loglik_a, loglik_b, *, k_a=0, k_b=0, correction='none')

Vuong’s test for non-nested model selection.

Compares two models by their pointwise log-likelihoods. Under the null that both are equally close to the truth, the statistic sqrt(n) * mean(m) / sd(m) (with m_i = ll_a_i - ll_b_i, minus an optional complexity correction) is asymptotically standard normal. A large positive value favors A.

Parameters:

• loglik_a (ndarray) – (n,) pointwise log-likelihoods of the two (non-nested) models.

• loglik_b (ndarray) – (n,) pointwise log-likelihoods of the two (non-nested) models.

• k_a (int) – parameter counts, used only if correction is set.

• k_b (int) – parameter counts, used only if correction is set.

• correction (str) – "none", "aic" (subtract k_a - k_b), or "bic" (subtract (k_a - k_b) log n / 2) from the log-likelihood ratio.

Returns:

{'statistic', 'p_value', 'favored'}.

Return type:

dict

clarke_test(loglik_a, loglik_b, *, k_a=0, k_b=0, correction='none')

Clarke’s distribution-free paired sign test for non-nested models.

Counts how often model A’s pointwise log-likelihood beats B’s; under the null this count is Binomial(n, 0.5). More robust than vuong_test() when the per-observation log-ratio is heavy-tailed or skewed (where the normal approximation behind Vuong fails).

Returns:

{'statistic', 'p_value', 'favored', 'n'} – statistic is the number of points favoring A.

Parameters:

• loglik_a (ndarray)

• loglik_b (ndarray)

• k_a (int)

• k_b (int)

• correction (str)

Return type:

dict

compare_elpd(pointwise_a, pointwise_b)

Compare two models’ expected log pointwise predictive density (LOO/WAIC).

Takes the per-observation elpd contributions (the pointwise arrays returned by mixle.ppl.diagnostics.psis_loo() / waic) and returns the elpd difference with the standard error of the pointwise difference – the standard-error estimate for model comparison (a difference within ~2 SE of zero is not decisive).

Parameters:

• pointwise_a (ndarray) – (n,) per-observation elpd contributions (higher is better).

• pointwise_b (ndarray) – (n,) per-observation elpd contributions (higher is better).

Returns:

{'elpd_diff', 'se', 'z', 'favored'} – elpd_diff = sum(a - b). A zero-variance pointwise difference is maximally significant when elpd_diff is nonzero, not an automatic tie – see the se == 0 handling below (same degeneracy as paired_score_difference()).

Return type:

dict

split_conformal(cal_pred, cal_y, test_pred, *, alpha=0.1, side='two-sided')

Split (inductive) conformal interval from a calibration set.

Parameters:

• cal_pred (ndarray) – (n,) model predictions on the calibration set.

• cal_y (ndarray) – (n,) calibration responses.

• test_pred (ndarray) – (m,) predictions at the test points.

• alpha (float) – miscoverage level (1 - alpha coverage).

• side (str) – "two-sided" (|y - yhat| score), "upper" (one-sided upper bound), or "lower" (one-sided lower bound).

Returns:

(lower, upper) arrays of length m (an unbounded side is -inf / +inf).

Return type:

tuple[ndarray, ndarray]

jackknife_plus(x, y, fit_predict, x_test, *, alpha=0.1)

Jackknife+ intervals (leave-one-out), using all data for both fitting and calibration.

For each training point i the model is refit without i; R_i = |y_i - mu_{-i}(x_i)| is the LOO residual and mu_{-i}(x) the LOO prediction at a test point. The interval aggregates mu_{-i}(x) -/+ R_i across i (Barber et al. 2021), giving ~``1 - 2 alpha`` worst-case and ~``1 - alpha`` typical coverage without a data split. Costs n refits.

Returns:

(lower, upper) arrays of length len(x_test).

Parameters:

• x (ndarray)

• y (ndarray)

• fit_predict (Callable[[ndarray, ndarray, ndarray], ndarray])

• x_test (ndarray)

• alpha (float)

Return type:

tuple[ndarray, ndarray]

cv_plus(x, y, fit_predict, x_test, *, alpha=0.1, n_folds=10, seed=0)

CV+ intervals: the K-fold analogue of jackknife_plus() (only n_folds refits).

Each point’s residual uses the model trained on the other folds, and the test prediction uses the same out-of-fold model. Much cheaper than Jackknife+ with nearly the same guarantee.

Returns:

(lower, upper) arrays of length len(x_test).

Parameters:

• x (ndarray)

• y (ndarray)

• fit_predict (Callable[[ndarray, ndarray, ndarray], ndarray])

• x_test (ndarray)

• alpha (float)

• n_folds (int)

• seed (int)

Return type:

tuple[ndarray, ndarray]

mondrian_conformal(cal_pred, cal_y, cal_groups, test_pred, test_groups, *, alpha=0.1)

Mondrian (group-conditional) split conformal: a separate quantile per group.

Calibrates the conformal quantile within each group (taxonomy), so coverage holds conditional on the group rather than only marginally – the fix when error scale varies across known subpopulations.

Parameters:

• cal_pred (ndarray) – calibration predictions, responses, and group labels.

• cal_y (ndarray) – calibration predictions, responses, and group labels.

• cal_groups (ndarray) – calibration predictions, responses, and group labels.

• test_pred (ndarray) – test predictions and their group labels.

• test_groups (ndarray) – test predictions and their group labels.

• alpha (float) – miscoverage level.

Returns:

(lower, upper) arrays of length len(test_pred).

Return type:

tuple[ndarray, ndarray]

weighted_conformal(cal_pred, cal_y, test_pred, weights, *, alpha=0.1, test_weight=None)

Covariate-shift-weighted split conformal (Tibshirani et al. 2019).

Under covariate shift the calibration and test inputs follow different distributions; reweighting the calibration scores by the likelihood ratio w(x) = p_test(x)/p_train(x) restores coverage – PER QUERY. The theorem places EACH test point’s own w(x_test) alongside the calibration weights, so test_weight is REQUIRED and per-query: pass an (m,) array of each query’s likelihood ratio, or a scalar as your explicit assertion that w(x_test) is the same for every query (true only when the shift does not move the test points’ own ratios). An earlier signature defaulted to a single shared scalar and suggested “the mean test/train ratio”; on a two-point-covariate fixture with exact ratios w(0)=0.505 / w(1)=50 that default covered 0.589 of 20,000 trials at nominal 0.90 while the query-specific weights covered 1.000 (STAT-RR21-07) – a shared scalar cannot apply the heavy query’s required mass.

Parameters:

• cal_pred (ndarray) – calibration predictions and responses.

• cal_y (ndarray) – calibration predictions and responses.

• test_pred (ndarray) – (m,) test predictions.

• weights (ndarray) – (n,) likelihood-ratio weights for the calibration points (need not be normalised).

• alpha (float) – miscoverage level, in [0.0, 1.0] (0.0 and 1.0 are valid boundaries, as for _conformal_quantile()).

• test_weight (float | ndarray | None) – (m,) per-query likelihood ratios w(x_test_i) (or a scalar, asserting a constant ratio across these queries). Required.

Returns:

(lower, upper) arrays of length m (a symmetric interval per test point, each at its own weighted quantile).

Return type:

tuple[ndarray, ndarray]

conformal_label_threshold(cal_prob_true, *, alpha=0.1)

Calibrate the LAC (least-ambiguous set-valued classifier) score threshold for 1 - alpha coverage.

The nonconformity score of a calibration point is 1 - p_model[true_class] – which needs the model’s class scores to rank well, not to be a true probability (the whole point: a softmax over a ReLU net is not a describable random process, but conformal still gives a finite-sample coverage guarantee from how those scores behave on held-out, exchangeable data). Returns the conformal quantile qhat of the calibration scores; a class is admitted at test time iff 1 - p[c] <= qhat (see conformal_label_sets()).

Parameters:

• cal_prob_true (ndarray) – (n,) model score assigned to the true class of each calibration point.

• alpha (float) – miscoverage level (1 - alpha marginal coverage of the returned sets).

Returns:

qhat – the score threshold (+inf when n is too small for the requested alpha).

Return type:

float

conformal_label_sets(cal_prob_true, test_prob, *, alpha=0.1, qhat=None)

Split-conformal prediction sets for a classifier: distribution-free 1 - alpha label coverage.

Calibrates a LAC threshold (conformal_label_threshold()) on the held-out true-class scores, then admits every class whose score clears it. The returned boolean mask has guaranteed marginal coverage: the true label is in the set with probability >= 1 - alpha. A singleton set is a confident prediction; an empty or multi-label set is an explicit abstention – the signal a cost-aware cascade escalates on.

Parameters:

• cal_prob_true (ndarray) – (n,) score assigned to the true class of each calibration point.

• test_prob (ndarray) – (m, K) model class scores at the test points (rows need not sum to 1).

• alpha (float) – miscoverage level.

• qhat (float | None) – a precomputed threshold (e.g. from an earlier calibration); recomputed if None.

Returns:

(sets, qhat) – sets is an (m, K) boolean mask, qhat the threshold used.

Return type:

tuple[ndarray, float]

learn_structure(data, *, field_estimators=None, min_gain=0.0, max_levels=64, n_bins=4, max_its=30, rng=None)

Discover the dependency forest for heterogeneous data and return the fitted joint model.

Any field can be a parent: a discrete one conditions directly, a continuous one is quantile-binned into n_bins conditioning levels (so a real can drive a count, a category, or another real). Scores every (parent -> child) pair by dependency_gain(), greedily builds a maximum-gain acyclic forest (each field at most one parent), and fits each factor. Falls back to independent marginals where no dependence clears min_gain – never worse than a composite, much better when structure exists. This is “automatic inference for composable models of heterogeneous data” made real.

rng seeds every internal fit’s EM initialization; None resolves to a FIXED seed, so two calls on the same data return the same model. (Before this knob, fits whose detected family needs a randomized init – a Gaussian-mixture conditional, say – drew fresh OS entropy per call, so the learned model itself was nondeterministic.)

Parameters:

• data (Sequence[tuple])

• field_estimators (Sequence[Any] | None)

• min_gain (float)

• max_levels (int)

• n_bins (int)

• max_its (int)

• rng (RandomState | None)

Return type:

DependencyTreeDistribution

learn_mixture_structure(data, n_components, *, restarts=3, max_iter=15, seed=0, min_gain=0.0, n_bins=4, max_its=30, field_estimators=None)

Fit a MixtureOfDependencyTrees by hard EM: discover clusters and each cluster’s dependency graph.

Each iteration re-learns a dependency forest per cluster on its currently-assigned points (M-step), then reassigns every record to its most-probable cluster (E-step), until assignments stabilize. Runs restarts random initializations and returns the highest-likelihood fit. Empty/tiny clusters receive a rescue fit and temporary proposal mass, but the returned weights count only actual hard memberships and may therefore contain an honest zero. Deterministic given seed: the one RandomState drives initialization, rescue, and every per-cluster fit.

field_estimators pins each field’s family (forwarded to learn_structure()) instead of re-running the automatic detector per cluster per iteration. Beyond the speedup, pinning matters for identifiability: the detector models a multimodal column with a Gaussian MIXTURE, which lets ONE cluster absorb what the caller intended as two – with per-field families pinned to unimodal models, regimes that differ in level must separate into different components to score well.

Parameters:

• data (Sequence[tuple])

• n_components (int)

• restarts (int)

• max_iter (int)

• seed (int)

• min_gain (float)

• n_bins (int)

• max_its (int)

• field_estimators (Sequence[Any] | None)

Return type:

MixtureOfDependencyTrees

mixture_structure_health(mot, data, *, merged_sep_threshold=None)

The identifiability receipt for a fitted MixtureOfDependencyTrees.

The trap it names: a flexible per-field family (a Gaussian-mixture conditional, a heavy-tailed catch-all marginal) lets ONE component absorb what the caller intended as SEVERAL regimes – and because the absorbed fit scores well, likelihood-level receipts look healthy. Measured concretely: on two planted regimes, a one-component fit matches the two-component fit’s likelihood, so nothing downstream of the density can see the difference.

The check that can: STRUCTURE-CONDITIONAL multimodality. For every continuous field of every component, group that component’s dominated points by the field’s own learned conditioning (parent level for a binned conditional, regression residuals for a linear edge, the whole fiber for a root marginal) – after the tree has explained what it can, each group must be unimodal. A group that still splits (deterministic 2-means separation over the same 2.65 + 6/sqrt(n) finite-sample threshold the hvis merged-regime detector is calibrated to, minority share >= 20%) is a regime split the component is hiding, whichever family absorbed it.

Returns {"components": [...], "diagnosis": [str, ...]} – empty diagnosis means no component hides multimodal structure. Per component, multimodal_fields lists (field, where, separation) and mixture_factors lists factors the detector fitted as mixture families (supporting detail: where the absorbed structure went). The fix when it fires is more components or pinned unimodal field_estimators (see learn_mixture_structure()). Complementary to mixle.utils.hvis.model_fit_health, which audits density calibration and accepts this model directly.

Parameters:

• mot (MixtureOfDependencyTrees)

• data (Sequence[tuple])

• merged_sep_threshold (float | None)

Return type:

dict

learn_bayesian_network(data, *, max_parents=2, min_gain=0.0, max_its=30, weights=None, _routing=None)

Discover a heterogeneous DAG for data and return the fitted network.

Each field greedily gains up to max_parents parents by BIC-penalized conditional likelihood, keeping the graph acyclic. Continuous children become conditional-linear-Gaussian factors (regression on continuous + one-hot discrete parents); discrete/count children condition on the joint config of their discrete parents, or become GLM nodes when a driver is continuous. With weights (soft-EM responsibilities), every factor fit and the BIC search itself are responsibility-weighted, with effective sample size sum(weights).

A field the conditional-linear-Gaussian machinery cannot read as a real number – a string column with more than 64 distinct values, or a numeric column carrying missing values – is modeled by the same automatic detector mixle.utils.automatic.get_estimator() uses (a categorical/ignored leaf, or an Optional-wrapped one) and takes no part in the edge search; every other field still gets the full search. Such a column used to be handed a Gaussian template and crash inside the Gaussian encoder. See _column_routing(), which makes that decision; _routing is a private hook letting a caller that fits many components on SLICES of one dataset decide it once on the whole dataset instead.

Parameters:

• data (Sequence[tuple])

• max_parents (int)

• min_gain (float)

• max_its (int)

• weights (Sequence[float] | None)

• _routing (tuple[dict[int, int], list[bool], list[bool], list[Any], dict[int, list[Any]]] | None)

Return type:

HeterogeneousBayesianNetwork

learn_mixture_bayesian_network(data, n_components, *, restarts=3, max_iter=12, seed=0, max_parents=2, min_gain=0.0, max_its=30, em='hard')

Fit a MixtureOfBayesianNetworks by EM: discover clusters and each cluster’s DAG.

em="hard" (default): each iteration re-learns a network per cluster on its assigned points and reassigns every record to its most-probable cluster, until assignments stabilize. em="soft": proper EM – every record contributes to EVERY cluster with its responsibility, each component’s structure search and factor fits are responsibility-weighted (learn_bayesian_network(weights=)), and convergence is on the observed-data log-likelihood; boundary points shape both clusters instead of whipsawing between them. Both run from the same restarts (k-means-seeded + random); best final log-likelihood wins. Starved clusters are re-seeded (hard) / responsibility-floored (soft).

Parameters:

• data (Sequence[tuple])

• n_components (int)

• restarts (int)

• max_iter (int)

• seed (int)

• max_parents (int)

• min_gain (float)

• max_its (int)

• em (str)

Return type:

MixtureOfBayesianNetworks

select_mixture_components(data, k_values=(1, 2, 3, 4), *, em='hard', seed=0, **kwargs)

Model selection over the number of clusters by BIC.

Fits learn_bayesian_network() for k=1 and learn_mixture_bayesian_network() for each larger k, scores each by bayesian_network_bic(), and returns (best_model, report) where report = {"k": chosen, "bic": {k: score, ...}}. BIC’s n_params log n penalty is what stops a mixture of per-cluster DAGs from always preferring more clusters.

Parameters:

• data (Sequence[tuple])

• k_values (Sequence[int])

• em (str)

• seed (int)

• kwargs (Any)

Return type:

tuple[Any, dict[str, Any]]

bayesian_network_bic(model, data)

BIC (lower is better) for a fitted network or mixture: -2 LL + n_params log n.

Parameters:

• model (Any)

• data (Sequence[tuple])

Return type:

float

learn_structure_embedded(data, *, text_fields='auto', n_clusters=8, embed_dim=16, seed=0, max_parents=2, **embed_kw)

Discover cross-field structure where some fields are free text (see module docstring).

Parameters:

• data (Sequence[tuple]) – flat tuple records; text fields may hold arbitrary strings.

• text_fields (Sequence[int] | str) – which field indices to embed, or "auto" – every string-valued field with too many distinct values to be a categorical.

• n_clusters (int) – number of representative clusters surfaced by describe() (interpretability only; the model uses the full embedding vector, not a cluster code).

• embed_dim (int) – embedding dimension for mixle.represent.fit_embedder() – the vector node’s dim.

• **embed_kw (Any) – forwarded to fit_embedder (epochs, hidden, feature_dim, …).

• seed (int)

• max_parents (int)

• **embed_kw

Return type:

EmbeddedStructureModel

class EmbeddedStructureModel(net, codecs)

Bases: object

A discovered dependency graph over records whose text fields ride in as embedding VECTORS.

Parameters:

• net (Any)

• codecs (dict[int, EmbeddedFieldCodec])

encode_record(x)

The record with each text field replaced by its embedding vector.

Parameters:

x (tuple)

Return type:

tuple

encode_records(rows)

Replace configured text fields with embedding vectors for each record.

EmbeddedFieldCodec.vectors() guarantees one row per value, so the pairing below is total. It used to be a bare zip over an unchecked embedder result: an embedder that returned fewer vectors than it was given left the surplus records holding raw text in a slot the graph reads as numeric, and the batch silently shrank downstream instead of reporting that the embedder had underproduced.

Parameters:

rows (Sequence[tuple])

Return type:

list[tuple]

edges()

Return discovered graph edges.

Return type:

list[tuple[int, int]]

density_semantics()

Always LIKELIHOOD_FACTOR.

The score is exact for what it computes, but it is not a normalized law over records: text fields are replaced by embedding vectors before scoring, embedding is not a measure-preserving change of variables, and no Jacobian relates the two measures. Saying so in the docstring was not enough – the machine-readable contract (the ExactDensity capability, describe(), and every combinator that folds child semantics through join_density_semantics()) still read this as an exact density, so a model composed over one of these silently inherited an exactness claim.

log_density(x)

Log density of x in the embedded space, not of the original record.

Worth being exact about, because the name invites the stronger reading. Each configured text field is replaced by its embedding vector before scoring, and the graph then evaluates a density over that representation. Embedding is not a measure-preserving change of variables and no Jacobian relates the two, so this number is not the log density of the text record and does not integrate to one over any space of text.

What it is good for is comparison: between records whose non-text fields differ, and between records whose text embeds to genuinely different vectors. What it cannot support is an absolute probability claim about a record, or a comparison against a density from a model that scores text some other way. Two different texts that embed close together will score close together – that is the representation working as designed, not the model failing to tell them apart.

Parameters:

x (tuple)

Return type:

float

seq_log_density(rows)

Per-record log densities in the embedded space; see log_density() for what that means.

Parameters:

rows (Sequence[tuple])

Return type:

ndarray

describe()

The discovered structure with per-cluster representative examples for each text field.

Return type:

dict[str, Any]

class HeterogeneousBayesianNetwork(factors)

Bases: object

A DAG joint over a heterogeneous record: log p(x) = sum_i log P(x_i | parents(i)) over fitted factors.

Parameters:

factors (Sequence[Any])

edges()

Return DAG edges as (parent_field, child_field) pairs.

Return type:

list[tuple[int, int]]

describe(field_names=None)

A real, structured explanation of what the fit found – every field’s per-factor report (_MarginalFactor.describe() and friends: fitted regression coefficients, GLM weights, the conditional table, or the marginal summary), split into edges (fields with >=1 parent) and roots (independent fields). Pass field_names (e.g. the tuple position -> column name map) to label everything by name instead of integer index. Every value is read off this network’s own fitted parameters – nothing here is a placeholder or a re-derived approximation.

Parameters:

field_names (Sequence[str] | None)

Return type:

dict[str, Any]

log_density(x)

Evaluate the joint log density of one record.

Parameters:

x (tuple)

Return type:

float

seq_log_density(encoded)

Evaluate joint log density for encoded records.

Parameters:

encoded (Any)

Return type:

ndarray

dist_to_encoder()

Return the encoder for record batches consumed by seq_log_density.

Return type:

Any

sampler(seed=None)

Return a sampler for the fitted Bayesian network.

Parameters:

seed (int | None)

Return type:

Any

class MixtureOfBayesianNetworks(components, weights)

Bases: object

A latent mixture whose components each carry their own heterogeneous DAG (regression edges and all).

The fullest form of the moat: it discovers the clustering and each cluster’s cross-field graphical model, so the slope of one field on another (or the whole dependency graph) can differ between clusters – which a single network cannot represent. log p(x) = logsumexp_k(log w_k + log p_k(x)) over HeterogeneousBayesianNetwork components. Fit by learn_mixture_bayesian_network().

Parameters:

• components (Sequence[HeterogeneousBayesianNetwork])

• weights (Sequence[float])

log_density(x)

Evaluate the mixture joint log density of one record.

Parameters:

x (tuple)

Return type:

float

seq_log_density(encoded)

Evaluate mixture joint log density for encoded records.

Parameters:

encoded (Any)

Return type:

ndarray

dist_to_encoder()

Return the record encoder shared by all mixture components.

Return type:

Any

responsibilities(data)

Return posterior component probabilities for each record.

Parameters:

data (Sequence[tuple])

Return type:

ndarray

sampler(seed=None)

Return a sampler for the mixture of Bayesian networks.

Parameters:

seed (int | None)

Return type:

Any

property n_components: int

Return the number of mixture components.

dependency_gain(parent, child, child_estimator, *, max_its=30, penalty='bic', rng=None)

Description-length gain (nats) of modeling child conditioned on a discrete parent vs. independently.

Fits the marginal P(child) and the conditional P(child | parent) (a child model per parent value) on the same data and returns LL_cond - LL_marginal minus a complexity penalty for the extra parameters (BIC: 0.5 * (levels - 1) * k * ln n). Positive means the dependence is worth modeling. This is a model-based dependency test – it works across any pair of families, unlike a same-type MI estimate. rng seeds the fits’ EM initializations (None = a fixed seed: deterministic by default; matters when the child family needs a randomized init, e.g. a mixture).

Parameters:

• parent (Sequence[Any])

• child (Sequence[Any])

• child_estimator (Any)

• max_its (int)

• penalty (str)

• rng (RandomState | None)

Return type:

float

class DependencyTreeDistribution(parents, factors, binners=None)

Bases: FitProvenanceCarrier

A directed-forest joint over a heterogeneous record: each field is a marginal or a conditional on its parent.

log_density(record) = sum_root log P(f_root) + sum_child log P(f_child | f_parent). The dependence a CompositeDistribution assumes away is modeled here as per-parent-value conditionals – while it still scores, samples, and composes like any mixle distribution.

Parameters:

• parents (Sequence[int | None])

• factors (Sequence[Any])

• binners (Sequence[Any] | None)

log_density(x)

Evaluate the dependency-tree joint log density for one record.

Parameters:

x (tuple)

Return type:

float

seq_log_density(encoded)

Evaluate dependency-tree joint log density for encoded records.

Parameters:

encoded (Any)

Return type:

ndarray

dist_to_encoder()

Return the encoder for dependency-tree record batches.

Return type:

Any

sampler(seed=None)

Return a sampler for the dependency tree.

Parameters:

seed (int | None)

Return type:

Any

edges()

The learned dependency edges (parent_field, child_field).

Return type:

list[tuple[int, int]]

class MixtureOfDependencyTrees(components, weights)

Bases: object

A latent mixture whose components each carry their own discovered dependency structure.

log p(x) = logsumexp_k ( log w_k + log p_k(x) ) where each p_k is a DependencyTreeDistribution. This is the deep form of the tagline: it discovers both the clustering and the within-cluster cross-field dependence – so the same category can map to different reals in different clusters, which neither a single dependency tree (one relationship) nor a mixture of independent composites (no within-cluster dependence) can represent. Fit by learn_mixture_structure().

Parameters:

• components (Sequence[DependencyTreeDistribution])

• weights (Sequence[float])

property w: ndarray

Mixture-convention alias for weights – lets mixture-generic tooling (e.g. mixle.utils.hvis.model_fit_health / hvis_map) accept this model unchanged.

property log_w: ndarray

Mixture-convention alias for log_weights (see w).

log_density(x)

Evaluate mixture log density for one dependency-tree record.

Parameters:

x (tuple)

Return type:

float

seq_log_density(encoded)

Evaluate mixture log density for encoded dependency-tree records.

Parameters:

encoded (Any)

Return type:

ndarray

dist_to_encoder()

Return the record encoder shared by all tree components.

Return type:

Any

responsibilities(data)

Posterior p(component | record) for each record – the E-step and a soft cluster assignment.

Parameters:

data (Sequence[tuple])

Return type:

ndarray

sampler(seed=None)

Return a sampler for the mixture of dependency trees.

Parameters:

seed (int | None)

Return type:

Any

property n_components: int

Return the number of mixture components.

kfold(n, n_splits=5, *, shuffle=False, seed=0)

Standard k-fold split of n rows.

Parameters:

• n (int) – number of observations.

• n_splits (int) – number of folds k (each row is in exactly one test fold).

• shuffle (bool) – shuffle the row order before splitting (use for i.i.d. data; leave False to keep contiguous folds, i.e. blocked_kfold()).

• seed (int | RandomState | None) – RNG seed when shuffle is True.

Returns:

k (train_index, test_index) pairs.

Return type:

list[tuple[ndarray, ndarray]]

blocked_kfold(n, n_splits=5)

Contiguous-block k-fold (no shuffle) for serially dependent data.

Identical to kfold() with shuffle=False; named separately because keeping blocks contiguous is the point for time series, not an incidental default.

Parameters:

• n (int)

• n_splits (int)

Return type:

list[tuple[ndarray, ndarray]]

leave_one_out(n)

Leave-one-out CV: n folds, each holding out a single observation.

Parameters:

n (int)

Return type:

list[tuple[ndarray, ndarray]]

stratified_kfold(y, n_splits=5, *, shuffle=True, seed=0)

Stratified k-fold preserving class proportions in every fold.

Distributes each class’s indices across the folds so that class frequencies in each test fold mirror the overall frequencies – important for imbalanced classification.

Parameters:

• y (ndarray) – (n,) class labels.

• n_splits (int) – number of folds.

• shuffle (bool) – shuffle within each class before round-robin assignment.

• seed (int | RandomState | None) – RNG seed when shuffle is True.

Returns:

n_splits (train_index, test_index) pairs.

Return type:

list[tuple[ndarray, ndarray]]

leave_one_group_out(groups)

Leave-one-group-out CV: one fold per distinct group held out entirely.

Parameters:

groups (ndarray) – (n,) group labels (e.g. subject / site / genus id).

Returns:

One (train_index, test_index) pair per unique group.

Return type:

list[tuple[ndarray, ndarray]]

group_kfold(groups, n_splits=5)

Group k-fold: partition groups into k folds so no group spans train and test.

Greedily assigns whole groups (largest first) to the currently-smallest fold, balancing fold sizes while keeping each group intact – the right scheme for repeated measures when there are more groups than folds.

Parameters:

• groups (ndarray)

• n_splits (int)

Return type:

list[tuple[ndarray, ndarray]]

time_series_split(n, n_splits=5, *, gap=0, max_train_size=None)

Forward-chaining time-series CV: train always precedes test (expanding window).

The series is cut into n_splits + 1 contiguous blocks; fold i tests on block i+1 and trains on everything strictly before it. A gap (buffer) drops the gap points just before each test block from the training set so leakage through short-range autocorrelation is avoided.

Parameters:

• n (int) – number of time-ordered observations.

• n_splits (int) – number of train/test splits.

• gap (int) – number of observations to drop between the train and test segments.

• max_train_size (int | None) – cap on the training-window length (sliding window); None keeps it expanding.

Returns:

n_splits (train_index, test_index) pairs.

Return type:

list[tuple[ndarray, ndarray]]

purged_kfold(n, n_splits=5, *, embargo=0)

Purged/embargoed blocked k-fold (a.k.a. buffered CV).

Contiguous test blocks, but the embargo observations on each side of a test block are removed from the training set, so no training point sits within the autocorrelation reach of a test point. The standard guard for serially dependent data when you still want every block to serve as a test fold (cf. de Prado 2018).

Parameters:

• n (int) – number of time-ordered observations.

• n_splits (int) – number of contiguous folds.

• embargo (int) – buffer width removed from train on both sides of each test block.

Returns:

n_splits (train_index, test_index) pairs.

Return type:

list[tuple[ndarray, ndarray]]

spatial_block_kfold(coords, n_splits=5, *, block_size=None, n_side=None, seed=0)

Spatial-block k-fold: hold out contiguous spatial blocks, not scattered points.

Partitions space into a regular grid of blocks and assigns whole blocks to folds at random, so a test point’s spatial neighbours are held out with it rather than leaking into training. This is the right scheme for geostatistical / spatially autocorrelated data.

Parameters:

• coords (ndarray) – (n, d) spatial coordinates (typically d = 2).

• n_splits (int) – number of folds.

• block_size (float | None) – grid-cell width in coordinate units; if None it is derived from n_side.

• n_side (int | None) – number of grid cells per axis; defaults so there are comfortably more blocks than folds.

• seed (int | RandomState | None) – RNG seed for assigning blocks to folds.

Returns:

n_splits (train_index, test_index) pairs.

Return type:

list[tuple[ndarray, ndarray]]

nested_kfold(n, *, outer_splits=5, inner_splits=4, shuffle=False, seed=0)

Nested k-fold: an inner CV (for tuning) inside each outer fold (for evaluation).

The outer loop estimates generalisation; the inner loop selects hyper-parameters using only the outer-training data, so the reported score is not optimistically biased by tuning on the test set.

Parameters:

• n (int) – number of observations.

• outer_splits (int) – number of outer folds.

• inner_splits (int) – number of inner folds within each outer training set.

• shuffle (bool) – shuffle before splitting (i.i.d. data only).

• seed (int | RandomState | None) – RNG seed.

Returns:

A list of NestedFold; inner folds index the original array.

Return type:

list[NestedFold]

class NestedFold(train, test, inner)

Bases: object

One outer fold of nested_kfold().

Parameters:

• train (ndarray)

• test (ndarray)

• inner (list[tuple[ndarray, ndarray]])

train

outer training indices (used to build the inner CV).

Type:

numpy.ndarray

test

outer test indices (held out for the final, untuned-on evaluation).

Type:

numpy.ndarray

inner

inner (train_index, test_index) folds, indexing into the original array (not into train), so they can be applied directly.

Type:

list[tuple[numpy.ndarray, numpy.ndarray]]

log_score(prob, *, mean=True)

Logarithmic score (log loss) -log p(y) of the probability assigned to the outcome.

The log score is local (it depends only on the density at the realised value) and strictly proper. Pass the predictive probability/density evaluated at the observed outcome.

Parameters:

• prob (ndarray) – array of predictive probabilities (or densities) at the realised outcomes. Must be finite and non-negative – a probability or density can never be negative, NaN, or infinite. Zero is a legitimate (if poor) forecast, and is clipped away from zero before the log so it yields a large but finite penalty rather than inf.

• mean (bool) – if True (default) return the mean log loss; otherwise the per-observation vector.

Returns:

Mean log loss (float) or the per-observation array.

Raises:

ValueError – if any entry of prob is negative, NaN, or infinite. Clipping a negative value the same way as the legitimate zero-probability edge case would let an invalid input (e.g. -1) silently masquerade as a valid, if poor, forecast; and an infinite input would silently become the score -inf – which, this module being lower-is-better, looks like the single best score achievable, letting a broken upstream density silently win any comparison instead of being rejected.

Return type:

ndarray | float

brier_score(prob, outcome, *, mean=True)

Brier score (mean squared error of probabilistic classification).

For binary forecasts pass 1-D prob (predicted probability of the positive class) and 0/1 outcome. For K-class forecasts pass prob shaped (n, K) and outcome either as integer class labels (length n) or as a one-hot (n, K) matrix; the multiclass score sums the squared error across classes, so it ranges in [0, 2].

Parameters:

• prob (ndarray) – (n,) positive-class probabilities, or (n, K) class probabilities. Must be finite and lie in [0, 1] – these are genuine probabilities, not densities.

• outcome (ndarray) – (n,) 0/1 or integer labels, or (n, K) one-hot. Must be finite, and each entry a legitimate label for prob’s shape: 0/1 against a 1-D binary prob, an integer class index in [0, K) against an (n, K) prob, or 0/1 against an (n, K) one-hot matrix.

• mean (bool) – if True return the mean; otherwise the per-observation vector.

Returns:

Mean Brier score (float) or the per-observation array.

Raises:

ValueError – if prob is not finite or has entries outside [0, 1], or if outcome is not finite or has entries outside its valid domain. Without these checks an invalid prob/outcome pair can cancel out into a deceptively “perfect” score of 0.0 (e.g. brier_score([-1, 2], [-1, 2])), and an out-of-range integer label silently corrupts an unrelated one-hot slot instead of raising.

Return type:

ndarray | float

brier_decomposition(prob, outcome, *, bins=10)

Murphy’s reliability–resolution–uncertainty decomposition of the (binary) Brier score.

Bins the forecast probabilities into bins equal-width bins on [0, 1] and computes

Brier = reliability - resolution + uncertainty,

where reliability (lower is better) is the mean squared gap between forecast probability and observed frequency within a bin, resolution (higher is better) is how far bin frequencies move from the base rate, and uncertainty is the base-rate variance p_bar (1 - p_bar) (a property of the data, not the forecaster).

Parameters:

• prob (ndarray) – (n,) predicted probabilities of the positive class. Must be finite and lie in [0, 1].

• outcome (ndarray) – (n,) 0/1 outcomes. Must be finite and each entry exactly 0 or 1.

• bins (int) – number of equal-width probability bins.

Returns:

{'reliability', 'resolution', 'uncertainty', 'brier', 'brier_binned_reconstruction'}. brier is the ACTUAL Brier score mean((prob - outcome)^2) (audit SCORING-3: the key previously held the binned reconstruction, which differs from the real score by the within-bin variance term). brier_binned_reconstruction = reliability - resolution + uncertainty is the decomposition’s own sum. The components are PLUG-IN estimates whose finite-sample bias grows like bins/n (audit SCORING-2): with sparse bins, reliability is inflated by within-bin binomial noise and resolution by base-rate noise – read them against n/bins, and do not chase small differences at small samples.

Raises:

ValueError – if prob is not finite or has entries outside [0, 1], or if outcome is not finite or has entries outside {0, 1}. Without these checks an invalid input pair can silently produce a mathematically impossible negative “brier” value instead of raising (e.g. brier_decomposition([-1, 2], [-1, 2]) returns brier == -2.0).

Return type:

dict[str, float]

crps_ensemble(forecasts, y, *, fair=True, mean=True)

Continuous Ranked Probability Score from a finite predictive ensemble (sample).

Uses the energy form CRPS = E|X - y| - 1/2 E|X - X'| estimated from the ensemble draws. The CRPS generalises absolute error to distributional forecasts (a point forecast recovers |x-y|) and is reported in the same units as y.

Parameters:

• forecasts (ndarray) – predictive draws. Either (n, m) (m draws for each of n observations) or (m,) for a single observation. Ragged ensembles are not supported – pad or call per observation.

• y (ndarray) – (n,) realised values (or a scalar for the (m,) case).

• fair (bool) – if True (the DEFAULT – audit SCORING-1) use the unbiased 1/(m(m-1)) spread estimator (the “fair” CRPS), which is proper for finite ensembles: measured at m = 5, the expected fair score is minimized by reporting the true predictive (dispersion factor 1.000), while the 1/m^2 estimator rewards UNDER-dispersing to a factor of 0.700 – an improper score that trains and selects forecasters toward overconfidence. fair=False keeps the plain 1/m^2 estimator for comparability with sources that use it; treat those numbers as biased-low in ensemble spread.

• mean (bool) – if True return the mean CRPS; otherwise the per-observation vector.

Returns:

Mean CRPS (float) or the per-observation array.

Return type:

ndarray | float

crps_gaussian(mu, sigma, y, *, mean=True)

Closed-form CRPS for a Gaussian predictive distribution N(mu, sigma^2).

CRPS = sigma * [ z (2 Phi(z) - 1) + 2 phi(z) - 1/sqrt(pi) ] with z = (y - mu) / sigma (Gneiting & Raftery 2007). Exact, so it is the right reference when the forecast is Gaussian.

Parameters:

• mu (ndarray) – predictive means, broadcastable with y.

• sigma (ndarray) – predictive standard deviations (> 0), broadcastable with y.

• y (ndarray) – realised values.

• mean (bool) – if True return the mean; otherwise the per-observation vector.

Returns:

Mean CRPS (float) or the per-observation array.

Return type:

ndarray | float

interval_score(lower, upper, y, alpha, *, mean=True)

Winkler interval score for a central (1 - alpha) prediction interval.

IS = (u - l) + (2/alpha)(l - y) 1[y < l] + (2/alpha)(y - u) 1[y > u] (Gneiting & Raftery 2007). It rewards narrow intervals but penalises misses by 2/alpha times the shortfall, so at matched nominal coverage the tighter, better-placed interval scores lower. This is the proper way to rank interval methods (conformal, quantile-regression, GP credible bands) head to head.

Parameters:

• lower (ndarray) – (n,) interval lower endpoints.

• upper (ndarray) – (n,) interval upper endpoints.

• y (ndarray) – (n,) realised values.

• alpha (float) – miscoverage level (e.g. 0.1 for a 90% interval), in (0, 1).

• mean (bool) – if True return the mean; otherwise the per-observation vector.

Returns:

Mean interval score (float) or the per-observation array.

Return type:

ndarray | float

winkler_score(lower, upper, y, alpha, *, mean=True)

Alias for interval_score() (the score is due to Winkler 1972).

Parameters:

• lower (ndarray)

• upper (ndarray)

• y (ndarray)

• alpha (float)

• mean (bool)

Return type:

ndarray | float

pinball_loss(pred, y, tau, *, mean=True)

Pinball (quantile / check) loss for a quantile forecast.

rho_tau(u) = max(tau * u, (tau - 1) * u) with u = y - pred; its expectation is minimised by the true tau-quantile, so it is the proper score for quantile forecasts and the loss that quantile regression and conformalised quantiles optimise.

Parameters:

• pred (ndarray) – predicted tau-quantiles, shape (n,) (single level) or (n, q) (several levels, one per column of tau).

• y (ndarray) – (n,) realised values.

• tau (float | ndarray) – quantile level(s) in (0, 1) – a scalar for (n,) pred, or a length-q array matching the columns of a (n, q) pred.

• mean (bool) – if True return the mean over all entries; otherwise the per-observation (or per-observation-per-level) array.

Returns:

Mean pinball loss (float) or the per-observation array.

Return type:

ndarray | float

energy_score(forecasts, y, *, fair=True, mean=True)

Energy score: the multivariate generalisation of CRPS for vector-valued forecasts.

ES = E||X - y|| - 1/2 E||X - X'|| with Euclidean norms; for scalar y it equals crps_ensemble().

Parameters:

• forecasts (ndarray) – predictive draws shaped (n, m, d) (m d-vectors per observation) or (m, d) for a single observation.

• y (ndarray) – (n, d) realised vectors (or (d,) for the single-observation case).

• fair (bool) – if True (the DEFAULT – audit SCORING-1, same finite-ensemble impropriety as crps_ensemble()) use the unbiased 1/(m(m-1)) spread estimator; False keeps the biased 1/m^2 form for comparability, rewarding under-dispersion.

• mean (bool) – if True return the mean; otherwise the per-observation vector.

Returns:

Mean energy score (float) or the per-observation array.

Return type:

ndarray | float

skill_score(score, reference, *, perfect=0.0)

Skill score: fractional improvement of score over a reference score.

skill = (reference - score) / (reference - perfect) – 1 means the forecast hit the perfect score, 0 means it only matched the reference (e.g. climatology), and negative means it did worse than the reference. Works for any negatively-oriented score (lower is better) such as those above.

Parameters:

• score (float) – the forecast’s score (lower is better).

• reference (float) – the baseline/reference forecast’s score.

• perfect (float) – the score of a perfect forecast (0.0 for the rules in this module).

Returns:

The skill score (float); nan if the reference already equals the perfect score.

Return type:

float

select_best(candidates, *, score, lower_is_better=False, conformal_alpha=None, heuristic_alpha=None)

Score each candidate and return the best, the verifier-based test-time-compute selector.

Parameters:

• candidates (Any) – an iterable of candidate objects (anything score accepts).

• score (Callable[[Any], float]) – a verifier score(candidate) -> float; the winner maximizes it (or minimizes it when lower_is_better).

• lower_is_better (bool) – if True, the winner is the candidate with the smallest score.

• conformal_alpha (float | None) – deprecated compatibility alias for heuristic_alpha. It does not provide conformal coverage because candidate scores are neither held-out conformity scores nor exchangeable repeated measurements.

• heuristic_alpha (float | None) – optional normal score-spread threshold level in (0, 1). When given, confident is a heuristic ranking flag only; coverage_guarantee remains False. With exactly TWO candidates the heuristic is undefined – margin and spread are the same number up to sqrt(2), so the flag would be decided by alpha alone (never at 0.05, always at 0.20, whatever the scores) – and confident is returned False with confidence_method naming the degeneracy.

Returns:

A SelectionResult. It is also subscriptable (result["best"]), so callers may treat it as a small dict with keys best, best_index, scores, confident.

Raises:

ValueError – if candidates is empty, or conformal_alpha is outside (0, 1).

Return type:

SelectionResult

class SelectionResult(best, best_index, scores, confident=None, margin=None, band=None, confidence_method=None, coverage_guarantee=False, _extras=<factory>)

Bases: object

Result of a select_best() call.

Parameters:

• best (Any)

• best_index (int)

• scores (ndarray)

• confident (bool | None)

• margin (float | None)

• band (float | None)

• confidence_method (str | None)

• coverage_guarantee (bool)

• _extras (dict[str, Any])

best

the winning candidate (the actual object, not its index).

Type:

Any

best_index

the position of the winner in the input candidates.

Type:

int

scores

per-candidate scores in input order (a numpy float array).

Type:

numpy.ndarray

confident

whether the winner’s lead clears the optional score-spread heuristic — None when neither heuristic_alpha nor its deprecated alias was supplied.

Type:

bool | None

margin

the winner’s lead over the runner-up (in score units, always >= 0); None when there is a single candidate.

Type:

float | None

band

the heuristic band the margin was compared against — None when no threshold was requested or there is a single candidate.

Type:

float | None

confidence_method

the executed diagnostic, "normal_score_spread_heuristic", or None.

Type:

str | None

coverage_guarantee

always False; one verifier score per selected candidate cannot establish held-out conformal or bootstrap coverage.

Type:

bool

class ConjugateUpdatable

Bases: PredicateCapability

Has a closed-form conjugate Bayesian update (the top tier of the inference ladder).

supports(x, ConjugateUpdatable) is True iff conjugate_posterior(x, data) returns an exact closed-form posterior; otherwise Bayesian inference falls back to the numerical fitters (MAP / Laplace / MCMC / VI). Detection is the single conjugate_posterior registry.

classmethod check(obj)

Return whether obj has a registered closed-form conjugate update.

Parameters:

obj (Any)

Return type:

bool

conjugate_posterior(dist, data, prior=None, weights=None)

Closed-form conjugate posterior over the parameters of dist given data.

Every supported family returns an exact closed-form parameter posterior: exact parameter samples, marginal likelihood, and posterior mean / point estimate. A posterior predictive is exact where its compound law is implemented; otherwise the returned distribution carries is_exact_posterior_predictive=False and a machine-readable approximation_receipt. For families with a multi-parameter likelihood whose conjugate is conditional (Gamma, InverseGamma, InverseGaussian, Pareto, NegativeBinomial, vonMises), the non-target parameter (shape / location / scale / number-of-trials / concentration) is taken as known from dist – exactly as a Binomial’s number of trials is. Families with no closed-form conjugate (full Beta, full Gamma, LogSeries, …) raise a clear error rather than returning a partial answer.

Parameters:

• dist – A mixle likelihood distribution instance whose type selects the conjugate family.

• data – A sequence of observations of the kind dist scores.

• prior (dict | None) – Optional dict of conjugate-prior hyperparameters (family specific). None uses a weak proper prior.

• weights (ndarray | None) – Optional per-observation weights (e.g. EM responsibilities).

Returns:

A ConjugatePosterior exposing mean, sample, point_estimate, log_marginal_likelihood, posterior_predictive and summary.

Return type:

ConjugatePosterior

mixture_conjugate_posterior(dist, data, priors, prior_weights=None, weights=None)

Posterior under a prior that is a mixture of conjugate priors (Diaconis-Ylvisaker).

Parameters:

• dist – The likelihood distribution (must map to a closed-form conjugate realiser, since the reweighting needs each component’s marginal likelihood).

• data – The observations.

• priors (list[dict]) – One hyperparameter dict per mixture component (same keys conjugate_posterior() accepts for this family).

• prior_weights (ndarray | None) – Prior mixing weights w_m (default uniform); normalised internally.

• weights (ndarray | None) – Optional per-observation weights.

Returns:

the exact posterior, again a mixture of conjugate posteriors with weights w'_m proportional to w_m * Z_m.

Return type:

A MixtureConjugatePosterior

is_conjugate_family(dist)

Return whether dist (instance or type) has a closed-form conjugate posterior.

Single source of truth: membership in the conjugate_posterior builder registry. This is the family-level capability — “can this distribution be updated in closed form?” — distinct from the instance-level has_conj_prior flag (whether a conjugate prior is currently attached for the MAP path). Backs mixle.capability.ConjugateUpdatable and ProbabilityDistribution.has_conjugate_prior().

Parameters:

dist (Any)

Return type:

bool

class ConjugatePosterior

Bases: object

Base class for a closed-form conjugate posterior over a likelihood’s parameters.

Subclasses provide the family-specific hyperparameters and the four capabilities: mean (the posterior mean of the parameters), sample (exact draws of the parameters), point_estimate (a fitted distribution at the posterior mean), log_marginal_likelihood (the evidence of the observed data under the prior), and posterior_predictive (the distribution of a new draw).

mean()

Return posterior mean parameters as a family-specific dictionary.

Return type:

dict[str, Any]

sample(n=1, rng=None)

Draw n exact parameter samples from the posterior.

Parameters:

• n (int)

• rng (RandomState | None)

Return type:

dict[str, ndarray]

sampler(seed=None)

Return a sampler exposing the standard obj.sampler(seed).sample(size) API.

Mirrors the distribution sampling convention so conjugate posteriors read the same way as every other mixle object; here each draw is a parameter set from the posterior. size=None returns one parameter set (scalars), size=n a dict of length-n arrays. The explicit-rng form sample(n, rng) remains available.

Parameters:

seed (int | None)

Return type:

ConjugatePosteriorSampler

point_estimate()

Return a likelihood distribution built from a representative posterior parameter value.

log_marginal_likelihood()

Return the closed-form log evidence of the observed data.

Return type:

float

posterior_predictive()

Return the posterior predictive distribution for a new observation.

summary()

Return a compact summary of the posterior family, mean, and hyperparameters.

Return type:

dict[str, Any]

hyper()

Return family-specific posterior hyperparameters.

Return type:

dict[str, Any]

class MixtureConjugatePosterior(components, post_weights, prior_weights, comp_log_evidence)

Bases: ConjugatePosterior

Posterior under a prior that is itself a mixture of conjugate priors.

Diaconis-Ylvisaker (1979): if the prior is sum_m w_m * pi_m(theta) with each pi_m conjugate to the likelihood, the posterior is again a mixture of the component conjugate posteriors, sum_m w'_m * pi_m^post(theta), with the mixing weights reweighted by each component’s marginal likelihood Z_m:

w’_m proportional to w_m * Z_m.

Everything stays closed-form, so a prior can be multimodal (e.g. “the rate is near 1 OR near 10”) or robust (a heavy-tailed prior as a mixture of conjugates) without losing exact inference.

mean()

Return the posterior-weighted mean of component posterior means.

Return type:

dict[str, Any]

sample(n=1, rng=None)

Draw parameter samples from the posterior mixture.

Parameters:

• n (int)

• rng (RandomState | None)

Return type:

dict[str, ndarray]

point_estimate()

The maximum-a-posteriori component’s point estimate (the dominant conjugate mode).

posterior_predictive()

A mixle MixtureDistribution of the component predictives, weighted by the posterior weights.

log_marginal_likelihood()

Return the mixture-prior log evidence by log-summing component evidences.

Return type:

float

hyper()

Return posterior weights and component hyperparameters.

Return type:

dict[str, Any]

class FisherView(dist, estimator=None)

Bases: object

Accumulator-backed Fisher-geometry view of a distribution.

Parameters:

• dist (Any) – mixle.stats or mixle.bstats distribution.

• estimator (Any | None) – Optional estimator. When omitted, dist.estimator() is used.

Notes

The generic encoded-data path obtains per-row statistics by replaying seq_update with one-hot weights. That keeps the interface compatible with existing encoders, but it is a correctness fallback rather than a high-performance implementation. Important families should override to_fisher() with direct seq_expected_statistics kernels.

structured_statistics(x, estimate=None, weight=1.0)

Structured sufficient stats for one observation.

For latent-variable distributions, this is the posterior-expected complete-data sufficient statistic under estimate (or this view’s distribution when estimate is omitted).

Parameters:

• x (Any)

• estimate (Any | None)

• weight (float)

Return type:

Any

expected_structured_statistics(x, estimate=None, weight=1.0)

Alias for posterior-expected structured sufficient statistics.

Parameters:

• x (Any)

• estimate (Any | None)

• weight (float)

Return type:

Any

sufficient_statistics(x, estimate=None, vectorizer=None)

Return vectorized sufficient statistics for one observation.

Parameters:

• x (Any)

• estimate (Any | None)

• vectorizer (SufficientStatisticVectorizer | None)

Return type:

ndarray

expected_sufficient_statistics(x, estimate=None, vectorizer=None)

Alias for vectorized posterior-expected sufficient statistics.

Parameters:

• x (Any)

• estimate (Any | None)

• vectorizer (SufficientStatisticVectorizer | None)

Return type:

ndarray

seq_structured_statistics(enc_data, estimate=None)

Structured per-row stats from encoded data.

This generic implementation is intentionally conservative and may be slow for large data. It exists so every encoder-compatible model has a correct baseline.

Parameters:

• enc_data (Any)

• estimate (Any | None)

Return type:

list[Any]

statistics_matrix(data=None, enc_data=None, estimate=None, vectorizer=None, fit=True)

Return an n x d matrix of per-observation sufficient statistics.

Parameters:

• data (Sequence[Any] | None)

• enc_data (Any | None)

• estimate (Any | None)

• vectorizer (SufficientStatisticVectorizer | None)

• fit (bool)

Return type:

ndarray

expected_statistics_matrix(data=None, enc_data=None, estimate=None, vectorizer=None, fit=True)

Alias for a matrix of posterior-expected sufficient statistics.

Parameters:

• data (Sequence[Any] | None)

• enc_data (Any | None)

• estimate (Any | None)

• vectorizer (SufficientStatisticVectorizer | None)

• fit (bool)

Return type:

ndarray

seq_expected_statistics(enc_data, estimate=None, vectorizer=None, fit=True)

Return posterior-expected statistics for encoded observations.

Parameters:

• enc_data (Any)

• estimate (Any | None)

• vectorizer (SufficientStatisticVectorizer | None)

• fit (bool)

Return type:

ndarray

mean_statistics(stats=None, **kwargs)

Return empirical mean sufficient statistics.

Parameters:

• stats (ndarray | None)

• kwargs (Any)

Return type:

ndarray

score_center(stats=None, **kwargs)

Return the center used for score or Fisher-vector whitening.

Parameters:

• stats (ndarray | None)

• kwargs (Any)

Return type:

ndarray

fisher_information(stats=None, diagonal=False, ridge=1.0e-8, **kwargs)

Empirical Fisher approximation from per-observation statistic vectors.

Parameters:

• stats (ndarray | None)

• diagonal (bool)

• ridge (float)

• kwargs (Any)

Return type:

ndarray

fisher_vectors(stats=None, metric='diagonal', center=None, fisher=None, ridge=1.0e-8, **kwargs)

Return centered/whitened sufficient-statistic vectors.

metric=’identity’ returns centered statistics, ‘diagonal’ divides by per-coordinate Fisher standard deviations, and ‘full’ applies an empirical full-matrix whitening transform.

Parameters:

• stats (ndarray | None)

• metric (str)

• center (ndarray | None)

• fisher (ndarray | None)

• ridge (float)

• kwargs (Any)

Return type:

ndarray

observed_fisher_information(stats=None, diagonal=False, center=None, ridge=1.0e-8, **kwargs)

Observed Fisher estimate from score vectors for observed data.

For latent models the statistic rows are posterior-expected complete-data sufficient statistics, so centering them by their model expectation gives observed score vectors. Their covariance is the observed Fisher metric used by Fisher-vector embeddings. If a model center is unavailable, this falls back to the empirical mean.

Parameters:

• stats (ndarray | None)

• diagonal (bool)

• center (ndarray | None)

• ridge (float)

• kwargs (Any)

Return type:

ndarray

observed_fisher_vectors(stats=None, metric='diagonal', center=None, fisher=None, ridge=1.0e-8, **kwargs)

Return centered statistics whitened by an observed Fisher metric.

Parameters:

• stats (ndarray | None)

• metric (str)

• center (ndarray | None)

• fisher (ndarray | None)

• ridge (float)

• kwargs (Any)

Return type:

ndarray

fisher_vector(x, estimate=None, metric='diagonal', center=None, fisher=None, vectorizer=None, ridge=1.0e-8)

Return a single Fisher vector using a model or caller-supplied reference geometry.

A generic view cannot estimate either a center or Fisher metric from one observation: doing so centers the observation by itself and produces a zero vector by construction.

Parameters:

• x (Any)

• estimate (Any | None)

• metric (str)

• center (ndarray | None)

• fisher (ndarray | None)

• vectorizer (SufficientStatisticVectorizer | None)

• ridge (float)

Return type:

ndarray

natural_parameters()

Return natural parameters when a specialized view provides them.

Return type:

Any

class FixedFisherView(dist, labels)

Bases: FisherView

Distribution-specific Fisher view with fixed vector coordinates.

Parameters:

• dist (Any)

• labels (Sequence[Path])

structured_statistics(x, estimate=None, weight=1.0)

Return fixed-coordinate sufficient statistics for one observation.

Parameters:

• x (Any)

• estimate (Any | None)

• weight (float)

Return type:

Any

sufficient_statistics(x, estimate=None, vectorizer=None)

Return projected sufficient statistics for one observation.

Parameters:

• x (Any)

• estimate (Any | None)

• vectorizer (SufficientStatisticVectorizer | None)

Return type:

ndarray

seq_structured_statistics(enc_data, estimate=None)

Return fixed-coordinate statistics for encoded observations.

Parameters:

• enc_data (Any)

• estimate (Any | None)

Return type:

list[Any]

statistics_matrix(data=None, enc_data=None, estimate=None, vectorizer=None, fit=True)

Return projected fixed-coordinate statistics for raw or encoded data.

Parameters:

• data (Sequence[Any] | None)

• enc_data (Any | None)

• estimate (Any | None)

• vectorizer (SufficientStatisticVectorizer | None)

• fit (bool)

Return type:

ndarray

mean_statistics(stats=None, model=True, **kwargs)

Return model-mean or empirical sufficient statistics.

Parameters:

• stats (ndarray | None)

• model (bool)

• kwargs (Any)

Return type:

ndarray

fisher_information(stats=None, diagonal=False, ridge=1.0e-8, **kwargs)

Return model Fisher information in fixed coordinates.

Parameters:

• stats (ndarray | None)

• diagonal (bool)

• ridge (float)

• kwargs (Any)

Return type:

ndarray

fisher_vectors(stats=None, metric='diagonal', center=None, fisher=None, ridge=1.0e-8, **kwargs)

Return centered or whitened Fisher vectors in fixed coordinates.

Parameters:

• stats (ndarray | None)

• metric (str)

• center (ndarray | None)

• fisher (ndarray | None)

• ridge (float)

• kwargs (Any)

Return type:

ndarray

to_fisher(dist, **kwargs)

Return a FisherView for dist via its own to_fisher hook.

Fisher views are co-located with each distribution: a distribution owns its view by overriding ProbabilityDistribution.to_fisher. This module keeps only the shared base machinery (FisherView/FixedFisherView/SufficientStatisticVectorizer and the reusable CountFisherView / EmpiricalMetricFixedFisherView helpers) plus _legacy_to_fisher(), the type-name dispatch for families not yet migrated to a per-file hook (and the generic fallback).

Parameters:

• dist (Any)

• kwargs (Any)

Return type:

FisherView

nuts(target, *, backend='auto', dim=None, init=None, num_samples=1000, warmup=1000, chains=1, mass=1.0, target_accept=0.8, max_tree_depth=10, thin=1, rng=None, parallel=None, **backend_kwargs)

No-U-Turn Sampler over an arbitrary log-target, dispatched to a registered engine.

The target contract depends on the backend (the kinds cannot be auto-converted across autodiff systems):

• numpy / numba: a fused value_and_grad(theta) -> (logp, grad) (njit-jitted for numba). The caller supplies the (analytic) gradient.

• torch / jax: a scalar logp(theta); the backend builds value_and_grad by autodiff (torch.func / NumPyro).

Parameters:

• target (Callable[[...], Any]) – the log-target, per the contract above.

• backend (str) – "auto" (default), or one of available_backends(). "auto" honors an explicit choice elsewhere, otherwise prefers the always-present numpy path for a plain numpy value_and_grad; pass backend="numba" to run an @njit target, "jax" / "torch" for an autodiff scalar logp.

• dim (int | None) – parameter dimension. Required unless init is given.

• init (Any) – initial state, shape (dim,) or (chains, dim) for per-chain starts. Defaults to zeros. A 1-D init is reused (jittered after the first) across chains.

• num_samples (int) – retained post-warmup draws per chain.

• warmup (int) – step-size adaptation / burn-in iterations per chain.

• chains (int) – number of independent chains (>= 2 to get a meaningful R-hat).

• mass (Any) – forwarded to the sampler.

• target_accept (float) – forwarded to the sampler.

• max_tree_depth (int) – forwarded to the sampler.

• thin (int) – forwarded to the sampler.

• rng (RandomState | int | None) – seed / RandomState for reproducibility.

• parallel (bool | str | None) – run the chains independent chains concurrently. None/False (default) runs them in the backend’s usual single call; "thread" uses a thread pool (real speedup for the numba/torch backends, whose inner loops release the GIL); True or "process" uses a process pool (real speedup for the pure-numpy backend; requires a picklable target). Ignored for the jax backend, which already vectorizes chains internally. Each chain still gets an independent seed, so results are valid multi-chain draws (R-hat / ESS across chains).

• **backend_kwargs (Any) – forwarded to the chosen backend (e.g. compile=, device= for torch; chain_method= for jax).

Returns:

NutsResult with pooled samples (chains*draws, d), per-chain chains (chains, draws, d), per-dimension rhat and ess, the total target-evaluation count, the adapted step_size, and extra={"backend": name}.

Return type:

NutsResult

nuts_torch(logp, *, dim=None, init=None, num_samples=1000, warmup=1000, chains=1, mass=1.0, target_accept=0.8, max_tree_depth=10, thin=1, rng=None, compile=True, dtype=None, device=None)

Torch-native NUTS over a torch scalar logp(theta) -> Tensor[()] (GPU / large autodiff targets).

The whole leapfrog/tree runs in torch tensors with a torch.compile``d ``value_and_grad (no numpy round-trip / graph re-trace per gradient). Same multi-chain interface and NutsResult as nuts(). Also registered as the "torch" backend.

Performance note (be deliberate about when to use this): on CPU this is typically slower than the numpy nuts() (per-op torch dispatch + host syncs in the tree dominate when the target is low-cost). Its value is GPU (device=) and large autodiff targets. For CPU work prefer the numpy backend, numba (analytic gradient), or the jax/NumPyro backend (XLA, vectorized multi-chain). Chains run in a Python loop (per-chain latency).

Parameters:

• logp (Callable[[Any], Any])

• dim (int | None)

• init (Any)

• num_samples (int)

• warmup (int)

• chains (int)

• mass (Any)

• target_accept (float)

• max_tree_depth (int)

• thin (int)

• rng (RandomState | int | None)

• compile (bool)

• dtype (Any)

• device (Any)

Return type:

NutsResult

advi(target_batch, u0, s0, *, samples=1000, mc=16, steps=2000, lr=0.05, batch_size=None, family='meanfield', alpha=1.0, rng=None)

Automatic-differentiation VI over a user batched torch target.

Parameters:

• target_batch (Callable[[Any], Any]) – target_batch(U: Tensor(B, d)) -> Tensor(B,) returning the (unnormalized) joint log-target for each of B parameter draws. The caller owns any data minibatching/rescaling inside this callable; batch_size here is unused unless the caller wires it in (kept for signature parity with the internal ADVI).

• u0 – initial variational mean and (marginal) scale in the unconstrained space.

• s0 – initial variational mean and (marginal) scale in the unconstrained space.

• samples (int) – number of posterior draws to return.

• mc (int) – Monte-Carlo samples per step, optimizer steps, Adam learning rate.

• steps (int) – Monte-Carlo samples per step, optimizer steps, Adam learning rate.

• lr (float) – Monte-Carlo samples per step, optimizer steps, Adam learning rate.

• family (str) – 'meanfield' (diagonal) or 'fullrank' (Cholesky covariance).

• alpha (float) – Renyi/tilted objective exponent (1.0 = standard KL-ELBO).

• rng (RandomState | int | None) – seed / RandomState.

• batch_size (int | None)

Returns:

AdviResult with samples (samples, d) drawn from the fitted Gaussian q (in the same space target_batch consumes), plus the fitted mean and scale – and, for family='fullrank', the cholesky factor carrying the learned correlations.

Return type:

AdviResult

class NutsResult(samples, chains, rhat, ess, num_target_evals, step_size, extra=<factory>)

Bases: object

Draws and diagnostics from a multi-chain nuts() run.

samples is (chains * draws, d) pooled draws; chains is (n_chains, draws, d).

Parameters:

• samples (ndarray)

• chains (ndarray)

• rhat (ndarray)

• ess (ndarray)

• num_target_evals (int | None)

• step_size (float)

• extra (dict)

class AdviResult(samples, mean, scale, objective=None, objective_n_eval=None, cholesky=None)

Bases: object

Result of advi(): posterior draws plus the fitted variational parameters.

Parameters:

• samples (ndarray)

• mean (ndarray)

• scale (ndarray)

• objective (float | None)

• objective_n_eval (int | None)

• cholesky (ndarray | None)

objective: float | None = None

Final variational objective, a Monte Carlo ESTIMATE from objective_n_eval draws.

For |alpha - 1| <= 1e-6 this is the ELBO. Otherwise it is the K-sample tilted Renyi bound at K = objective_n_eval, whose expectation depends on K – compare two fits’ objectives only at equal objective_n_eval – and for alpha > 1 the finite-K estimate is NOT a lower bound on the evidence (it is upward-biased and can exceed log Z). The value comes from a FIXED number of Adam steps with no convergence criterion; finiteness, not stationarity, is what was checked.

cholesky: ndarray | None = None

Lower-triangular Cholesky factor of the fitted covariance, for family='fullrank' only.

None for meanfield, whose covariance is exactly diag(scale**2). A full-rank fit exists to learn the posterior correlations, and those live entirely in the off-diagonals: without this the returned parameters describe a diagonal Gaussian, so reconstructing q from (mean, scale) silently gives back the meanfield answer the caller explicitly asked not to have.

rhat(chains)

Gelman-Rubin potential scale reduction factor (R-hat) per parameter dimension.

R-hat compares the variance between chains to the variance within chains. Values near 1.0 indicate the chains have mixed and are sampling a common target; values above ~1.01-1.1 flag non-convergence. Standard multi-chain check (Gelman & Rubin 1992).

Parameters:

chains (Any) – array shaped (n_chains, n_draws) or (n_chains, n_draws, d). Needs at least two chains and two draws.

Returns:

A length-d array of R-hat values (one per parameter).

Return type:

ndarray

ess(samples, max_lag=None)

Effective sample size per parameter, n / tau with Geyer’s initial-monotone-sequence tau.

Accepts a single chain (n_draws,) / (n_draws, d) or a stack of chains (n_chains, n_draws, d) (the chains are pooled into one autocorrelation estimate after centering each chain by its own mean). Adjacent-lag autocorrelation pairs are summed until the first non-positive pair, independently per parameter, with the retained pairs forced non-increasing (Geyer 1992) – see _geyer_tau().

Returns:

A length-d array of ESS values.

Parameters:

• samples (Any)

• max_lag (int | None)

Return type:

ndarray

split_rhat(chains)

Rank-normalized split-R-hat (Vehtari et al. 2021) – the robust, recommended R-hat.

Splits each chain in half (doubling the chain count, so within-chain non-stationarity is caught), rank-normalizes the pooled draws (robust to heavy tails / non-normality), then applies the Gelman-Rubin R-hat. Convergence is typically declared at < 1.01.

Parameters:

chains (Any)

Return type:

ndarray

folded_split_rhat(chains)

Folded rank-normalized split-R-hat (Vehtari et al. 2021): split-R-hat on |x - median(x)|.

The plain split_rhat() compares chain locations; folding about the median makes it compare chain scales/tails, catching the case where chains share a mean but differ in spread. Use together with split_rhat() (or rhat_max()).

Parameters:

chains (Any)

Return type:

ndarray

rhat_max(chains)

The recommended convergence R-hat: max(split_rhat, folded_split_rhat) (Vehtari et al. 2021).

A single number that flags non-convergence in either the location (bulk) or the scale (folded) of the chains; declare convergence at < 1.01.

Parameters:

chains (Any)

Return type:

ndarray

mcse_mean(chains)

Monte Carlo standard error of the posterior mean: sd(x) / sqrt(ESS) per parameter.

The sampling error in the estimated mean from autocorrelated draws – the posterior standard deviation deflated by the (autocorrelation-based) effective sample size of the raw chains.

Parameters:

chains (Any)

Return type:

ndarray

geweke_z(chain, first=0.1, last=0.5)

Geweke convergence z-score comparing the start and end of a single chain (per parameter).

Tests stationarity within one chain (Geweke 1992): if the chain has converged, the mean of its first first fraction and its last last fraction agree, so the z-statistic (mean_a - mean_b) / sqrt(se_a^2 + se_b^2) – with autocorrelation-adjusted standard errors (the segment variance divided by its effective sample size) – is approximately standard normal. |z| < 2 indicates convergence; a large |z| flags a chain still drifting. Returns one z per parameter. Complements the multi-chain split_rhat().

Parameters:

• chain (Any)

• first (float)

• last (float)

Return type:

ndarray

mcmc_summary(chains)

Per-parameter posterior + convergence summary (the ArviZ-style summary table).

Returns one dict per parameter with the posterior mean/sd and 5/50/95% quantiles, plus the recommended convergence diagnostics: r_hat (= rhat_max(), the max of the bulk and folded rank-normalized split-R-hats), ess_bulk, ess_tail, and mcse_mean (Vehtari et al. 2021). Declare convergence at r_hat < 1.01 with ess_bulk/ess_tail comfortably in the hundreds.

Parameters:

chains (Any)

Return type:

list[dict[str, float]]

ess_bulk(chains)

Bulk effective sample size: ESS of the rank-normalized split chains (mixing in the distribution body).

Parameters:

chains (Any)

Return type:

ndarray

ess_tail(chains, prob=0.05)

Tail effective sample size: the smaller ESS of the lower-/upper-prob tail indicators.

Tail quantiles mix more slowly than the bulk; tail-ESS is min(ESS[1(x <= q_prob)], ESS[1(x >= q_{1-prob})]) over the split chains (Vehtari et al. 2021), surfacing poor tail exploration that bulk-ESS misses.

Parameters:

• chains (Any)

• prob (float)

Return type:

ndarray

available_backends()

Return the names of registered backends whose engine is importable, in registration order.

Each probe is isolated: one backend whose available() raises used to propagate out of this call and hide every healthy backend, so a single broken third-party registration made the always-present numpy path look unavailable too. A raising probe answers “unknown”, which is reported as a warning and treated as unavailable – it cannot be silently equivalent to a probe that returned False, because the caller may need to know its engine is misconfigured rather than absent.

Return type:

list[str]

class InferenceBackend(name, available, target_kind, nuts)

Bases: object

A registered inference engine: a name, an availability probe, a target contract, a sampler.

Parameters:

• name (str)

• available (Callable[[], bool])

• target_kind (str)

• nuts (Callable[..., NutsResult])

register_inference_backend(backend, *, replace=False)

Register an inference backend under backend.name.

The registry is global process state that select_backend() dispatches every nuts call through, so one malformed registration is not a local mistake – it corrupts dispatch for every caller. A backend with an empty name and an unrecognized target_kind used to register successfully and become selectable, even though select_backend() rejects that same kind string when a caller supplies it.

Replacing an already-registered name requires replace=True. Silent replacement let a second import or a typo’d name swap out a live sampler with no trace; an intentional override (a test double, a vendored engine) says so explicitly.

Raises:

• TypeError – if backend is not an InferenceBackend, or available/nuts are not callable.

• ValueError – for an empty/blank name, an unrecognized target_kind, or an unrequested replacement of a registered name.

Parameters:

• backend (InferenceBackend)

• replace (bool)

Return type:

None

class EventStudyResult(effect, se, z, p_value, ci, treated_mean, treated_se, control_mean, control_se, tau2_treated, n_treated, n_control, shrunk_treated, estimand, identified, interpretation, identification, ci_level=0.95, n_treated_population=None, n_control_population=None, df=None)

Bases: object

Pooled change contrast with explicit causal-identification metadata.

Parameters:

• effect (float)

• se (float)

• z (float)

• p_value (float)

• ci (tuple[float, float])

• treated_mean (float)

• treated_se (float)

• control_mean (float | None)

• control_se (float | None)

• tau2_treated (float)

• n_treated (int)

• n_control (int)

• shrunk_treated (ndarray)

• estimand (str)

• identified (bool)

• interpretation (str)

• identification (dict[str, Any] | None)

• ci_level (float)

• n_treated_population (int | None)

• n_control_population (int | None)

• df (float | None)

class EventStudyIdentification(design_evidence, parallel_trends_evidence, exchangeability, positivity, consistency, no_interference, no_anticipation, sensitivity_analysis=None)

Bases: object

Auditable assumptions for interpreting a treated/control contrast as causal DiD.

Parameters:

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

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

• exchangeability (bool)

• positivity (bool)

• consistency (bool)

• no_interference (bool)

• no_anticipation (bool)

• sensitivity_analysis (str | None)

to_dict()

Return a serializable copy of the receipt.

Return type:

dict[str, Any]

gaussian_effect(pre, post)

Per-subject mean shift and its (Welch) sampling variance from pre/post activity samples.

Parameters:

• pre (ndarray)

• post (ndarray)

Return type:

tuple[float, float]

poisson_lograte_effect(k_pre, t_pre, k_post, t_post)

Per-subject log activity-rate shift log(rate_post) - log(rate_pre) for event counts over windows.

k_* are event counts, t_* the window durations (or exposures). The estimate is the Haldane form log((k_post + 0.5)/t_post) - log((k_pre + 0.5)/t_pre) with the matching plug-in variance 1/(k_post + 0.5) + 1/(k_pre + 0.5) – STRICTLY POSITIVE counts only: a zero-count window is REFUSED with a pointer to poisson_pooled_rate_ratio() (STAT-RR17-09: at zero the Haldane offset is pure correction, not rate information, and pooling it drove a true null to p = 1.37e-12). An earlier version of this docstring claimed zero counts were made finite and quoted a 1/k variance; neither matched the executable contract (STAT-P20-03).

Parameters:

• k_pre (float)

• t_pre (float)

• k_post (float)

• t_post (float)

Return type:

tuple[float, float]

poisson_lograte_effects(k_pre, t_pre, k_post, t_post)

Per-subject log rate-ratio effects for Gaussian pooling, via the CONDITIONAL (binomial) route.

Conditioning is what makes the heterogeneous estimand reachable. Given a subject’s total n_i = k_pre,i + k_post,i and the common exposure ratio r = t_post/t_pre, k_post,i | n_i ~ Binomial(n_i, p_i) with logit(p_i) = log(theta_i) + log(r) – the subject’s BASELINE RATE lambda_i cancels exactly, so nothing here assumes rate homogeneity. Two earlier constructions failed measurably on that point (pass 19 and its pass-20 replay):

1. Count-estimated per-subject variances (1/(k+0.5) forms) let inverse-variance weights correlate with the effects (measured corr(1/v, y) = -0.72; weighted true-null mean -0.150; z = -7.8 at n = 1000, 100% false rejection, flat in n).

2. Debiasing at the pooled ARM means fixed the homogeneous fixture but assumed every subject shares one baseline rate: with half the subjects at rate 0.1 and half at 9.1 (exposures 1:3) and every true ratio exactly 1, it rejected 400/400 with mean z -20.75/-65.62 at n = 1e3/1e4 (STAT-RR19-03) – the nonlinear Haldane bias must be removed at each subject’s own law, which is unknowable at the rate scale and EXACT at the conditional scale.

Here each effect is the Haldane conditional logit log((k_post+0.5)/(k_pre+0.5)) - log(r), debiased by the exact pmf-summation bias of that statistic under Binomial(n_i, p_bar) (p_bar = pooled working proportion; under the null every subject’s p_i equals it exactly, whatever the baseline rates), with the exact conditional variance at (n_i, p_bar) as the pooling variance – a function of n_i alone given p_bar, so weights cannot correlate with the per-subject noise. Zero-total subjects (n_i = 0) carry no information about a rate ratio and are excluded from the returned arrays.

Measured on the pass-20 fixtures through the real DL pool: heterogeneous 0.1/9.1 (r = 3) true null rejects 0.055 at n = 1000 (400 reps) and 0.0417 at n = 10000 (1200 reps); the milder 2/8 (r = 2) mixture 0.018; the homogeneous 4.6 (r = 3) fixture 0.048; power at a true common ratio of 1.3 on the 2/8 mixture is 1.000. The 4.0 arm-mean floor remains as the pooled-z normal-approximation gate; below it use poisson_pooled_rate_ratio, exact at any sparsity for the common-ratio estimand.

Returns (effects, variances, selection_receipt). The receipt records the POPULATION this route can speak for (STAT-RR22-01): zero-total subjects are excluded because the conditional likelihood is empty for them, and that exclusion is OUTCOME-DEPENDENT – under extreme baseline heterogeneity (rates 0.01 and 10, true effects +1/-1 with all-subject ATT exactly zero), event-free low-rate subjects vanished differentially and the retained-subject contrast read -0.94 with 100% rejection while silently wearing the all-treated label. Hand the receipt to hierarchical_event_study() so the estimand names the population actually estimated.

Parameters:

• k_pre (Any)

• t_pre (float)

• k_post (Any)

• t_post (float)

Return type:

tuple[ndarray, ndarray, dict[str, int]]

poisson_pooled_rate_ratio(k_pre, t_pre, k_post, t_post, *, alpha=0.05)

EXACT conditional inference for a COMMON post/pre rate ratio – valid at any sparsity.

Sums the counts over subjects: with a common ratio theta and a constant per-subject exposure ratio r = t_post/t_pre, K_post | K_post + K_pre = N is exactly Binomial(N, theta*r / (1 + theta*r)), so the test of theta = 1 is an exact binomial test and the CI is the Clopper-Pearson interval transformed back to the ratio scale (STAT-RR17-09’s sparse-regime route). ESTIMAND: the common ratio – a deliberate reduction from the per-subject random-effects mean, which sparse counts cannot support; heterogeneity across subjects is NOT estimated here.

Parameters:

• k_pre (Any)

• t_pre (float)

• k_post (Any)

• t_post (float)

• alpha (float)

Return type:

dict[str, Any]

hierarchical_event_study(treated_effects, treated_vars, control_effects=None, control_vars=None, *, alpha=0.05, identification=None, treated_selection=None, control_selection=None)

Pool per-subject changes into an association or explicitly identified DiD estimate.

*_effects / *_vars are per-subject changes and sampling variances. A control group makes the numeric contrast treated_mean - control_mean. It is labeled an ATT only when a complete identification receipt is attached; otherwise it remains a treated/control change association.

THE WEIGHTING FOLLOWS THE ESTIMAND (STAT-RR21-01). The ATT is the treated-population ARITHMETIC mean of unit effects, so the identified path uses equal-subject-weight group means with their empirical standard errors (sd/sqrt(n) per group – per-subject sampling noise is inside the observed spread, so this SE covers noise and heterogeneity together). DerSimonian-Laird precision weighting estimates a DIFFERENT quantity whenever effect size and sampling variance are dependent – with treated unit effects evenly split between +1 and -1 (true ATT exactly 0), precision weights track the totals that track the effects, and the DL contrast reported +0.197 with 92-100% rejection at nominal 5%. DL pooling remains what the ASSOCIATION paths report (a precision-weighted summary, named as such), and its tau^2 / shrunk-effect diagnostics are attached in both modes.

THE POPULATION FOLLOWS THE RECEIPT (STAT-RR22-01). The conditional Poisson route excludes zero-total subjects – an OUTCOME-DEPENDENT exclusion, so its arithmetic contrast estimates the event-positive population, not all treated units (measured: with baseline rates 0.01/10 and true all-subject ATT exactly zero, differential exclusion drove the retained-subject contrast to -0.94 with 100% rejection). Pass each arm’s selection_receipt (the third return of poisson_lograte_effects()); when any subjects were excluded, the identified estimand names the event-positive population explicitly and the result carries the original denominators. Refusing to attach the all-treated label is the point: this route cannot identify it.

IDENTIFIED-PATH UNCERTAINTY (STAT-RR22-02): the reference is Student-t with Welch- Satterthwaite degrees of freedom (a normal quantile rejected 18.96% at nominal 5% with 2 subjects per arm), the empirical variance is FLOORED by the supplied per-subject sampling variances (mean(vars)/n per arm – a degenerate all-identical-estimates arm used to report SE 0 and the contradiction CI=[1,1] with p=1), the p-value and CI come from the same t reference, and subjects are assumed INDEPENDENT units: with 10 clusters repeated 10 times per arm, rejection measured 55.6% at nominal 5% – aggregate to cluster-level effects first.

Parameters:

• treated_effects (ndarray)

• treated_vars (ndarray)

• control_effects (ndarray | None)

• control_vars (ndarray | None)

• alpha (float)

• identification (EventStudyIdentification | None)

• treated_selection (dict[str, int] | None)

• control_selection (dict[str, int] | None)

Return type:

EventStudyResult

tipping_drift(result)

Sensitivity bound: the unmeasured differential drift that would explain the effect away.

Within-subject DiD is unbiased only if, absent treatment, treated and control would have drifted equally. This returns the differential drift delta (in effect units) that nullifies the estimate (= effect) and the value that pushes the result’s OWN confidence interval through zero – at the result’s own level and reference (STAT-RR23-03: a hard-coded Normal 1.96 disagreed with the identified path’s Welch-t interval edge, 0.1235 vs the actual 0.0604, and silently priced a 90% result at 95%). Larger = more robust.

Parameters:

result (EventStudyResult)

Return type:

dict

Subpackages

• mixle.inference.mcmc package
  • Submodules
    • mixle.inference.mcmc.conjugate module
    • mixle.inference.mcmc.gradients module
    • mixle.inference.mcmc.nuts_numba module
    • mixle.inference.mcmc.nuts_torch module
    • mixle.inference.mcmc.parameter_bridge module
    • mixle.inference.mcmc.proposals module
    • mixle.inference.mcmc.samplers module
• mixle.inference.production package
  • Submodules
    • mixle.inference.production.drift module
    • mixle.inference.production.monitor module
    • mixle.inference.production.provenance module
    • mixle.inference.production.registry module
    • mixle.inference.production.serving module

Submodules

• mixle.inference.backend_respecialization module
• mixle.inference.backends module
• mixle.inference.bayesian_network module
• mixle.inference.belief module
• mixle.inference.blackbox module
• mixle.inference.block_em module
• mixle.inference.block_gibbs module
• mixle.inference.calibrate_fit module
• mixle.inference.calibration module
• mixle.inference.calibration_gate module
• mixle.inference.causal module
• mixle.inference.condition module
• mixle.inference.conditional_jit_controller module
• mixle.inference.conformal module
• mixle.inference.copula_structure module
• mixle.inference.create module
• mixle.inference.cross_validation module
• mixle.inference.decision module
• mixle.inference.diagnostics module
• mixle.inference.em module
• mixle.inference.errors_in_variables module
• mixle.inference.estimation module
• mixle.inference.event_study module
• mixle.inference.explain module
• mixle.inference.fisher module
• mixle.inference.forecast module
• mixle.inference.freeze_rollup module
• mixle.inference.fusion_policy module
• mixle.inference.glm module
• mixle.inference.gradient_fit module
• mixle.inference.heterogeneous_executor module
• mixle.inference.integrity module
• mixle.inference.jit module
• mixle.inference.leaf_hotswap module
• mixle.inference.model_comparison module
• mixle.inference.multiple_testing module
• mixle.inference.node_precision_plan module
• mixle.inference.node_report module
• mixle.inference.nonparametric module
• mixle.inference.objectives module
• mixle.inference.orchestration module
• mixle.inference.ordinal module
• mixle.inference.placement module
• mixle.inference.planning module
• mixle.inference.posterior module
• mixle.inference.precision_plan module
• mixle.inference.price_forecast module
• mixle.inference.priors module
• mixle.inference.project module
• mixle.inference.receipt module
• mixle.inference.refine module
• mixle.inference.reproduce module
• mixle.inference.resampling module
• mixle.inference.risk module
• mixle.inference.robust module
• mixle.inference.scenario module
• mixle.inference.scoring module
• mixle.inference.select module
• mixle.inference.simulate module
• mixle.inference.skill module
• mixle.inference.spark_executor module
• mixle.inference.streaming module
• mixle.inference.structure module
• mixle.inference.structure_embedded module
• mixle.inference.survival module
• mixle.inference.synthesize module
• mixle.inference.target module
• mixle.inference.torsion module
• mixle.inference.transaction module
• mixle.inference.uncertainty module
• mixle.inference.uq module