mixle.inference.resampling module

Bootstrap and permutation inference for arbitrary statistics.

Distribution-free uncertainty estimates by resampling the data itself:

  • bootstrap() – a confidence interval for any statistic T(data), with the resampling scheme matched to the data’s dependence structure: plain i.i.d., stratified (resample within groups), cluster/hierarchical (resample whole clusters – the unit of independence), moving block (preserve autocorrelation in a series), or m-out-of-n subsampling. Interval types: percentile, basic (pivotal), and bca (bias-corrected and accelerated – the second-order-accurate default for i.i.d. data).

  • wild_bootstrap() – residual bootstrap for regression that is robust to heteroscedasticity (Rademacher or Mammen two-point multipliers on the residuals).

  • permutation_test() – an exact/Monte-Carlo test for an arbitrary statistic under a sharp null, with stratified / restricted (within-group) shuffling and a paired (sign-flip) mode; when the number of distinct rearrangements is small it enumerates them for an exact p-value. “Exact” always means exact FOR THE STATED SHARP NULL, carried on the result as null_hypothesis (STAT-RR17-12): the sign-flip null is SYMMETRY of the paired differences about the null value, not mean equality (a mean-zero asymmetric alternative rejected 43% at nominal 5%), and the two-sample null is full exchangeability, not equal means under unequal variances.

Everything is pure NumPy. data may be a single array (resampled along axis 0) or a tuple of arrays sharing their first axis (e.g. (X, y)); the statistic is then called as statistic(*parts).

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

Bases: object

Result of a bootstrap() call.

Parameters:
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

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

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

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

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

BootstrapResult

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

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

class PermutationResult(statistic, pvalue, null_distribution, n_perm, exact, alternative, null_hypothesis='')[source]

Bases: object

Result of a permutation_test().

Parameters:
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

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

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