mixle.analysis package ¶

hill_estimator(data, k)[source]

Hill estimator of the tail index xi = 1/alpha from the top k order statistics.

xi_hat = (1/k) sum_{i=1}^{k} log(X_(n-i+1) / X_(n-k)) – consistent for heavy (Pareto-type, xi > 0) tails. For a Pareto tail with exponent alpha this estimates 1/alpha.

Parameters:

data (ndarray)
k (int)

Return type:

moment_estimator(data, k)[source]

Dekkers–Einmahl–de Haan moment estimator of the extreme-value index (any tail sign).

Generalises Hill to xi of either sign by combining the first two log-moments of the top k exceedances; works for heavy, light, and bounded (xi < 0) tails.

Parameters:

data (ndarray)
k (int)

Return type:

mean_residual_life(data, thresholds)[source]

Mean-excess (mean-residual-life) function for POT threshold selection.

e(u) = mean(X - u | X > u). Over a range where the GPD fits, e(u) is approximately linear in u (slope xi/(1-xi)); the lowest threshold from which the plot is linear is the choice.

Returns:

{'threshold', 'mean_excess', 'n_exceed'}.

Parameters:

data (ndarray)
thresholds (ndarray)

Return type:

endpoint_estimator(data, k, *, method='gpd')[source]

Estimate the finite right endpoint of a bounded support (frontier / boundary estimation).

Fits a GPD to the top k exceedances over X_(n-k); when the tail index xi is negative the support is bounded and the endpoint is x+ = X_(n-k) - beta/xi (which, by the GPD support constraint, necessarily exceeds the observed maximum). Generic to econometric frontier analysis, reliability limits, and image-edge localisation. Returns inf if the estimated tail is unbounded (xi >= 0).

Parameters:

data (ndarray) – the sample.
k (int) – number of upper order statistics (exceedances) used.
method (str) – "gpd" – GPD-MLE endpoint.

Returns:

The estimated right endpoint (inf if unbounded).

Return type:

n_records(data)[source]

Number of upper records. For an i.i.d. sequence of length n the expectation is H_n.

Parameters:: data (ndarray)
Return type:: int

record_times(data)[source]

Indices at which a new running maximum (upper record) occurs, including the first observation.

Parameters:: data (ndarray)
Return type:: ndarray

class KDE(data, *, bandwidth='silverman', bounds=None, adaptive=False)[source]

Bases: object

A fitted kernel density estimate.

Use kde() to construct. Evaluate with evaluate() (or call the instance). Supports a Gaussian kernel, reflection boundary correction (bounds), and adaptive bandwidths.

Parameters:

data (np.ndarray)
bounds (tuple[float | None, float | None] | None)
adaptive (bool)

evaluate(x)[source]

Density at points x (with reflection boundary correction if bounds was set).

Parameters:: x (ndarray)
Return type:: ndarray

kde(data, *, bandwidth='silverman', bounds=None, adaptive=False)[source]

Construct a kernel density estimate (Gaussian kernel).

Parameters:

data (ndarray) – (n,) sample.
bandwidth – "silverman", "scott", or a positive float.
bounds – (lo, hi) support limits for reflection boundary correction; either may be None for an unbounded side (e.g. (0.0, None) for a positive variable).
adaptive (bool) – use Abramson variable bandwidths (wider where the pilot density is low).

Returns:

A KDE.

Return type:

KDE

kde_mode(data, *, bandwidth='silverman', bounds=None, grid=None, ci=False, n_boot=500, ci_level=0.95, seed=0)[source]

Estimate the mode (peak location) of a density, optionally with a bootstrap CI.

Parameters:

data (ndarray) – (n,) sample.
bandwidth – passed to kde().
bounds – passed to kde().
grid (ndarray | None) – evaluation grid; defaults to 512 points spanning the data range.
ci (bool) – if True return a percentile bootstrap interval for the mode.
n_boot (int) – bootstrap controls.
ci_level (float) – bootstrap controls.
seed (int | RandomState | None) – bootstrap controls.

Returns:

The mode (float), or {'mode', 'ci_low', 'ci_high'} when ci is True.

Return type:

float | dict

intensity(events, grid, *, bandwidth='silverman', domain=None, edge_correct=True)[source]

Kernel intensity lambda(t) of an inhomogeneous Poisson / point process.

Unlike a density (which integrates to 1), the intensity integrates to the expected number of events: lambda_hat(t) = sum_i K_h(t - t_i). With edge_correct the estimate is divided by the fraction of the kernel falling inside domain, removing the downward bias near the boundary.

Parameters:

events (ndarray) – (m,) event locations.
grid (ndarray) – points t at which to evaluate the intensity.
bandwidth – "silverman", "scott", or a float.
domain (tuple[float, float] | None) – (lo, hi) observation window (defaults to the event range); used for edge correction.
edge_correct (bool) – divide by the in-window kernel mass at each t.

Returns:

The intensity evaluated on grid.

Return type:

ndarray

silverman_bandwidth(data)[source]

Silverman’s rule-of-thumb bandwidth 0.9 min(sd, IQR/1.34) n^{-1/5} (1-D).

Parameters:: data (ndarray)
Return type:: float

scott_bandwidth(data)[source]

Scott’s rule-of-thumb bandwidth sd * n^{-1/(d+4)}.

Parameters:: data (ndarray)
Return type:: float

turing_coverage(counts)[source]

Turing’s sample-coverage estimate and the complementary unseen probability mass.

C = 1 - f1 / n where f1 is the number of singletons and n the total count: the estimated probability that the next observation is a previously seen species. 1 - C = f1/n is the Good–Turing estimate of the total probability of all unseen species.

Returns:: {'coverage', 'unseen_mass', 'n', 'f1'}.
Parameters:: counts (ndarray)
Return type:: dict[str, float]

good_turing(counts)[source]

Simple Good–Turing smoothed probabilities (Gale & Sampson 1995).

Reallocates probability from seen to unseen items using the frequencies of frequencies. Empirical Turing estimates r* = (r+1) N_{r+1}/N_r are used for small r and a smoothed log-linear fit S(r) takes over once the two diverge (the Gale switch), giving stable discounts in the sparse tail.

Parameters:: counts (ndarray) – per-species abundances (zeros ignored).
Returns:: {'p0', 'proba', 'r_star', 'r'} – p0 is the total probability assigned to unseen species; proba are the smoothed probabilities of the input species (aligned to the positive entries of counts, summing to 1 - p0); r_star / r are the discounted and raw frequencies for the distinct abundance classes.
Return type:: dict[str, ndarray | float]

chao1(counts, *, ci_level=0.95)[source]

Chao1 nonparametric richness estimator from abundance data (bias-corrected).

S_chao1 = S_obs + f1 (f1 - 1) / (2 (f2 + 1)) (Chao 1984, bias-corrected form), a lower bound on total richness driven by the singleton (f1) and doubleton (f2) counts. Returns a standard error and a log-normal confidence interval for the number of undetected species (Chao 1987), so the interval respects S_chao1 >= S_obs.

Returns:

{'estimate', 'observed', 'f1', 'f2', 'se', 'ci_low', 'ci_high'}.

Parameters:

counts (ndarray)
ci_level (float)

Return type:

chao2(incidence, *, ci_level=0.95)[source]

Chao2 richness estimator from replicated incidence (presence/absence) data.

Parameters:

incidence (ndarray) – (n_species, n_sites) 0/1 matrix (or per-species counts of sites occupied, passed as a 1-D array together with ``… `` – a 2-D matrix is expected here).
ci_level (float) – confidence level for the log-normal interval.

Returns:

{'estimate', 'observed', 'q1', 'q2', 'se', 'ci_low', 'ci_high', 'sites'} where q1/q2 are the numbers of species found in exactly one / two sites.

Return type:

ace(counts, *, rare_threshold=10)[source]

ACE: Abundance-based Coverage Estimator of richness (Chao & Lee 1992).

Splits species into abundant (> rare_threshold) and rare, estimates sample coverage from the rare group’s singletons, and corrects for the coefficient of variation of the rare abundances.

Returns:

{'estimate', 'observed', 's_rare', 's_abund', 'c_ace'}.

Parameters:

counts (ndarray)
rare_threshold (int)

Return type:

ice(incidence, *, rare_threshold=10)[source]

ICE: Incidence-based Coverage Estimator of richness (the Chao–Lee estimator for incidence data).

Parameters:

incidence (ndarray) – (n_species, n_sites) 0/1 matrix.
rare_threshold (int) – species found in <= rare_threshold sites are treated as infrequent.

Returns:

{'estimate', 'observed', 's_infreq', 's_freq', 'c_ice'}.

Return type:

hill_numbers(counts, q=(0.0, 1.0, 2.0))[source]

Hill numbers (effective number of species) of order q.

The unified diversity profile: q=0 is observed richness, q=1 is the exponential of Shannon entropy, q=2 is the inverse Simpson concentration. Larger q weights common species more, so the profile D(q) summarises evenness as well as richness.

Parameters:

counts (ndarray) – per-species abundances.
q (float | ndarray) – a scalar order or an array of orders.

Returns:

Array of Hill numbers, one per requested order (scalar input still returns a length-1 array).

Return type:

ndarray

rarefaction_curve(counts, sizes=None)[source]

Individual-based rarefaction: expected richness when subsampling m individuals (Hurlbert).

E[S(m)] = sum_i (1 - C(n - x_i, m) / C(n, m)) – the expected number of species seen in a random subsample of m of the n individuals. Used to compare richness between samples at a common sample size (or coverage).

Parameters:

counts (ndarray) – per-species abundances.
sizes (ndarray | None) – subsample sizes m to evaluate; defaults to 1 .. n.

Returns:

{'sizes', 'expected_richness'}.

Return type:

class Variogram(model, nugget, psill, rng, nu=1.5, anisotropy=None)[source]

Bases: object

A fitted variogram model gamma(h) = nugget + psill * shape(h).

Parameters:

model (str)
nugget (float)
psill (float)
rng (float)
nu (float)
anisotropy (tuple[float, float] | None)

model

"spherical", "exponential", "gaussian" (aka "squared_exponential" / "rbf" – covariance psill * exp(-(h/rng)**2)), or "matern".

Type:: str

nugget

discontinuity at h=0 (measurement error / micro-scale variance).

Type:: float

psill

partial sill (correlated variance); nugget + psill is the total sill.

Type:: float

rng

range (correlation length).

Type:: float

nu

Matern smoothness (ignored by other models).

Type:: float

anisotropy

optional (angle_rad, ratio) geometric anisotropy – coordinates are rotated by angle and the minor axis scaled by 1/ratio before distances are taken.

Type:: tuple[float, float] | None

model: str

nugget: float

psill: float

rng: float

nu: float = 1.5

anisotropy: tuple[float, float] | None = None

gamma(h)[source]

Parameters:: h (ndarray)
Return type:: ndarray

cov_field(h)[source]

Covariance of the correlated field part (excludes the nugget): psill (1 - shape).

Parameters:: h (ndarray)
Return type:: ndarray

empirical_variogram(coords, values, *, n_bins=15, max_dist=None)[source]

Binned empirical (semi-)variogram: mean 0.5 (z_i - z_j)^2 by separation distance.

Returns:

{'lag', 'semivariance', 'count'} for each non-empty distance bin.

Parameters:

coords (ndarray)
values (ndarray)
n_bins (int)
max_dist (float | None)

Return type:

fit_variogram(coords, values, *, model='spherical', n_bins=15, nu=1.5)[source]

Fit a variogram model to data by least squares on the empirical variogram.

Returns:

A fitted Variogram (nugget, partial sill, range).

Parameters:

coords (ndarray)
values (ndarray)
model (str)
n_bins (int)
nu (float)

Return type:

Variogram

ordinary_kriging(coords, values, variogram, query, *, noise=None)[source]

Ordinary kriging: BLUP of the field at query under an unknown constant mean.

Parameters:

coords (ndarray) – (n, d) data locations.
values (ndarray) – (n,) measured field.
variogram (Variogram) – a fitted Variogram.
query (ndarray) – (q, d) prediction locations.
noise (ndarray | None) – optional (n,) per-observation measurement variance (heteroscedastic nugget); if None the homoscedastic variogram.nugget is used on the diagonal.

Returns:

{'prediction', 'variance'} arrays of length q.

Return type:

universal_kriging(coords, values, variogram, query, *, degree=1, noise=None)[source]

Universal kriging: kriging with a polynomial spatial trend (drift) of the given degree.

degree=1 removes a linear trend, degree=2 a quadratic one. Use when the field has a large-scale drift on top of the stationary residual the variogram describes.

Returns:

{'prediction', 'variance'}.

Parameters:

coords (ndarray)
values (ndarray)
variogram (Variogram)
query (ndarray)
degree (int)
noise (ndarray | None)

Return type:

calibrate_variance(predicted_var, residuals, *, target=0.9)[source]

Scale factor that makes kriging predictive intervals hit a target coverage.

Finds c so that standardised residuals residual / sqrt(c * predicted_var) achieve the target central coverage under a Gaussian predictive. Returns the variance multiplier c; multiply predicted_var by it to recalibrate (generic GP/kriging variance recalibration).

Parameters:

predicted_var (ndarray) – (m,) held-out kriging variances.
residuals (ndarray) – (m,) held-out actual - predicted.
target (float) – desired central coverage (e.g. 0.9).

Returns:

The variance multiplier c (> 0).

Return type:

borda_count(rankings)[source]

Borda positional aggregation: each item scores (m - 1 - position) summed over voters.

Returns:: {'consensus', 'scores'} – the consensus ordering (best first) and per-item Borda scores.
Parameters:: rankings (ndarray)
Return type:: dict[str, ndarray]

copeland(rankings)[source]

Copeland pairwise-majority aggregation.

Each ordered pair contributes a win/loss by majority across voters; an item’s score is wins minus losses. Closely tracks the Condorcet winner when one exists.

Returns:: {'consensus', 'scores', 'wins'}.
Parameters:: rankings (ndarray)
Return type:: dict[str, ndarray]

kemeny_consensus(rankings, *, exact_max_items=8)[source]

Kemeny median ranking: minimise the total Kendall-tau distance to all input orderings.

The Kemeny consensus is the maximum-likelihood aggregation under the Mallows–Kendall model and the Condorcet-consistent choice. Exact by enumeration when m <= exact_max_items; otherwise a local search (adjacent transpositions) from the Borda ordering.

Returns:

{'consensus', 'distance', 'exact'} – the consensus, its total Kendall distance, and whether the result is exact.

Parameters:

rankings (ndarray)
exact_max_items (int)

Return type:

dict

mallows_fit(rankings, *, exact_max_items=8)[source]

Fit a Mallows model (Kendall): central ranking + dispersion theta.

The central ranking is the kemeny_consensus(); theta is the MLE concentration, found by matching the observed mean Kendall distance to its expectation under the model. Larger theta means tighter agreement (theta -> 0 is uniform/no-consensus).

Returns:

{'center', 'theta', 'mean_distance'}.

Parameters:

rankings (ndarray)
exact_max_items (int)

Return type:

dict

kendall_distance(a, b)[source]

Kendall-tau distance: the number of item pairs ordered oppositely by a and b.

Parameters:

a (ndarray)
b (ndarray)

Return type:

int

spearman_footrule(a, b)[source]

Spearman footrule distance: the sum of absolute position differences across items.

Parameters:

a (ndarray)
b (ndarray)

Return type:

int

cayley_distance(a, b)[source]

Cayley distance: the minimum number of transpositions turning a into b (m - #cycles).

Parameters:

a (ndarray)
b (ndarray)

Return type:

int

class SpatialMixture(shape, n_components, emission, beta=1.0)[source]

Bases: object

A grid-structured mixture with a Potts prior on the latent labels and pluggable mixle emissions.

Parameters:

shape – grid shape, e.g. (nx, ny) or (nx, ny, nz) – defines the neighbour structure.
n_components (int) – number of mixture components (latent classes).
emission – a mixle ParameterEstimator for the per-component family, e.g. MultivariateGaussianEstimator() – this is what makes the class domain-agnostic.
beta (float) – Potts coupling (>= 0); larger smooths the labels more. 0 is an ordinary mixture.

fit(observations, *, max_iter=40, mf_iter=3, seed=0)[source]

Fit by mean-field variational EM. observations is a length-prod(shape) sequence of per-cell observations (row order matches shape.ravel()); each is a single emission datum.

Robustness: components are initialized by a short hard-assignment pass and the Potts coupling is annealed from 0 to beta over the first iterations, so components form before the smoothness prior is applied (a strong prior on a degenerate init otherwise collapses every cell into one).

Parameters:

max_iter (int)
mf_iter (int)
seed (int)

Return type:

SpatialMixture

responsibilities()[source]

The posterior label probabilities, (prod(shape), n_components) – a simplex per cell.

Return type:: ndarray

labels()[source]

The MAP label field, reshaped to shape.

Return type:: ndarray

entropy()[source]

Per-cell posterior entropy (label uncertainty), reshaped to shape.

Return type:: ndarray

component(j)[source]

The fitted mixle emission distribution of component j.

Parameters:: j (int)
Return type:: Any

class SmithMaxStable(sigma)[source]

Bases: object

The Smith max-stable process with Gaussian storm-profile covariance sigma (d x d, SPD).

A spatial process, not an i.i.d. leaf distribution: its full likelihood is intractable, so it is not a SequenceEncodableProbabilityDistribution – it exposes the things that do have closed forms (extremal_coefficient, bivariate_cdf) plus a sampler, and is fitted by the module-level fit_smith_maxstable() (composite/madogram estimation), mirroring the functional fit style of the other non-leaf spatial models. Margins are unit Frechet; spatial dependence grows with sigma.

Parameters:: sigma (np.ndarray)

extremal_coefficient(h)[source]

theta(h) = 2 * Phi(a/2) with a the Mahalanobis lag length – 1 at h=0 (full dependence) rising to 2 as the lag grows (independence).

Parameters:: h (ndarray)
Return type:: float

bivariate_cdf(z1, z2, h)[source]

P(Z(s) <= z1, Z(s+h) <= z2) = exp(-V(z1, z2)) – the Smith bivariate distribution.

Parameters:

z1 (float)
z2 (float)
h (ndarray)

Return type: