`statspai.doe`¶

doe ¶

Design of experiments: building designs and analysing them.

statspai.experimental covers the randomised controlled trial (assignment, balance, attrition, sample size). This package covers the design of the runs of an experiment or of a computer model: which combinations of factor levels to try.

factorial_design / design_aberration / factorial_effects: full and fractional two-level factorials, their alias structure, and the analysis of unreplicated experiments.
mixture_design: factors that are shares summing to one.
doe_optimal: D-, A- and I-optimal designs for linear, nonlinear and generalised linear models.
space_filling / design_augment / design_criteria: Latin hypercube, maximum projection, maximin and uniform designs.
support_points / split_data: representative points of a distribution and train / test splits built from them.
sobol_indices / morris_screening: global sensitivity analysis of a model. factor_importance: the same question asked of a data set, without a model.
sequential_design: Bayesian optimisation and active learning with a Gaussian process.

DesignResult `dataclass` ¶

Bases: ResultProtocolMixin

An experimental design: the runs, how it was built and how good it is.

Attributes:

Name	Type	Description
`design`	`DataFrame`	One row per run, one column per factor, in the units of the factors.
`unit`	`ndarray`	The same runs scaled to the unit cube.
`method`	`str`
`criteria`	`dict`	Space-filling measures of the design (see `sp.design_criteria`).
`model_info`	`dict`

Examples:

>>> import statspai as sp
>>> d = sp.space_filling(8, 2, method="lhs", seed=1)
>>> d.design.shape
(8, 2)

to_frame ¶

to_frame() -> DataFrame

The runs as a DataFrame.

plot ¶

plot(ax: Any = None, **kwargs: Any) -> Any

Scatter of the first two factors, or a scatter matrix for more.

FactorialEffectsResult `dataclass` ¶

Bases: ResultProtocolMixin

Effects of a two-level factorial experiment.

Attributes:

Name	Type	Description
`effects`	`DataFrame`	Indexed by term. `effect` (high minus low), `coef` (half of it, the regression coefficient on the -1 / +1 column), `half_normal_score`, `lenth_t`, `lenth_p`, `lenth_p_simultaneous`, `active` (beyond the margin of error); and `se`, `t`, `p` when the design leaves residual degrees of freedom.
`intercept`	`float`
`pse, me, sme`	`float`	Lenth's pseudo standard error, margin of error and simultaneous margin of error.
`df_resid`	`int`
`aliases`	`dict`	Terms that were dropped because their column repeats that of an earlier term.
`model_info`	`dict`

Examples:

>>> import statspai as sp
>>> d = sp.factorial_design(3).design
>>> d["y"] = [504, 984, 928, 808, 992, 784, 464, 976]
>>> fit = sp.factorial_effects(d, "y")
>>> fit.effects.shape[0]
7

active `property` ¶

active: List[str]

Terms beyond the margin of error.

plot ¶

plot(ax: Any = None, label: Optional[int] = None, **kwargs: Any) -> Any

Half-normal plot of the absolute effects (Daniel 1959).

Noise effects fall on a line through the origin; active effects lie to the right of it. label is the number of largest effects to name (default: the active ones, at least three).

FactorialDesignResult `dataclass` ¶

Bases: ResultProtocolMixin

A factorial design and its alias structure.

Attributes:

Name	Type	Description
`design`	`DataFrame`	One row per run in standard order (or run order when randomised). Two-level factors are coded -1 / +1 unless level values were given.
`kind`	`str`
`generators`	`list of str`	For a fractional design, e.g. `['E = ABCD']`.
`defining_relation`	`list of str`	Words equal to the identity, e.g. `['ABCDE']`; a leading minus sign marks a word that equals minus the identity. Reported for regular two-level fractions only: it is empty for a full factorial and also for the three-, five- and seven-level fractions and the orthogonal arrays, see the notes of :func:`factorial_design`.
`resolution`	`int or None`	Length of the shortest word; `None` for a full factorial.
`word_length_pattern`	`list of int`	Number of words of length 1, 2, 3, ...
`aliases`	`dict`	Each main effect and two-factor interaction with the effects it cannot be separated from (up to `alias_order` factors). A leading minus sign means the two columns are opposite: the estimate is the first effect minus the second. An interaction that is itself a defining word is constant over the fraction, hence not an effect, and is not listed.
`model_info`	`dict`

Examples:

>>> import statspai as sp
>>> d = sp.factorial_design(5, n_runs=16)
>>> d.generators, d.resolution
(['E = ABCD'], 5)

to_frame ¶

to_frame() -> DataFrame

The runs as a DataFrame.

FactorImportanceResult `dataclass` ¶

Bases: ResultProtocolMixin

Model-free importance of the columns of a data set.

Attributes:

Name	Type	Description
`importance`	`Series`	Total-index importance per factor; exactly zero for a factor that was not selected.
`selected`	`list of str`	Factors with positive importance, most important first.
`noise_share`	`float`	Share of the variance of the outcome that the selected factors leave unexplained.
`model_info`	`dict`

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame(rng.normal(size=(500, 4)), columns=list("abcd"))
>>> df["y"] = 2 * df["a"] + np.sin(2 * df["b"]) + 0.1 * rng.normal(size=500)
>>> res = sp.factor_importance(df, "y")
>>> sorted(res.selected)
['a', 'b']

plot ¶

plot(ax: Any = None, **kwargs: Any) -> Any

Bars of the importance, largest first.

ModelDesignResult `dataclass` ¶

Bases: ResultProtocolMixin

An optimal design for a model.

Attributes:

Name	Type	Description
`design`	`DataFrame`	The support points with a `weight` column (approximate design) or an `n` column of replicate counts (exact design).
`criterion`	`str`
`criterion_value`	`float`	`det(M)^(1/k)` for D (larger is better), `trace(M^-1) / k` for A and the average prediction variance for I (smaller is better), with `M` the information matrix per run.
`efficiency`	`float`	For an approximate design, a lower bound on its efficiency from the equivalence theorem (1 certifies optimality on the candidate set). For an exact design, its efficiency relative to the optimal approximate design.
`information`	`DataFrame`	Information matrix per run (at the first parameter vector).
`terms`	`list of str`
`model_info`	`dict`

Examples:

>>> import statspai as sp
>>> d = sp.doe_optimal("x + I(x**2)", {"x": (0, 2)})
>>> d.design.round(3).values.tolist()
[[0.0, 0.333], [1.0, 0.333], [2.0, 0.333]]

to_frame ¶

to_frame() -> DataFrame

One row per run for an exact design, the support points otherwise.

MorrisResult `dataclass` ¶

Bases: ResultProtocolMixin

Elementary-effects screening.

Attributes:

Name	Type	Description
`effects`	`DataFrame`	Indexed by input: `mu` (mean elementary effect), `mu_star` (mean absolute effect, the importance measure), `sigma` (standard deviation: non-linearity or interaction), and `share` (`mu_star` as a share of the total).
`elementary_effects`	`DataFrame`	One row per trajectory, one column per input.
`design`	`DataFrame`	The runs, in the units of the inputs.
`n_evaluations`	`int`
`model_info`	`dict`

Examples:

>>> import statspai as sp
>>> res = sp.morris_screening(lambda X: 3 * X[:, 0] + 0 * X[:, 1], 2,
...                           r=4, pass_as="array", seed=0)
>>> [round(float(v), 6) for v in res.effects["mu_star"]]
[3.0, 0.0]

plot ¶

plot(ax: Any = None, **kwargs: Any) -> Any

sigma against mu_star, one labelled point per input.

SobolResult `dataclass` ¶

Bases: ResultProtocolMixin

Variance-based sensitivity indices.

Attributes:

Name	Type	Description
`indices`	`DataFrame`	Indexed by input: `first` and `total` indices, with `_lower` / `_upper` bootstrap limits when `n_boot > 0`.
`variance`	`float`	Variance of the output.
`mean`	`float`
`n_evaluations`	`int`
`model_info`	`dict`

Examples:

>>> import statspai as sp
>>> res = sp.sobol_indices(lambda X: X[:, 0] + 2 * X[:, 1], 3,
...                        n=1024, pass_as="array", seed=0, n_boot=0)
>>> [abs(round(float(v), 1)) for v in res.indices["first"]]
[0.2, 0.8, 0.0]

plot ¶

plot(ax: Any = None, **kwargs: Any) -> Any

Bars of the first-order index and of the interaction part.

SequentialDesignResult `dataclass` ¶

Bases: ResultProtocolMixin

Runs chosen one at a time with a Gaussian process.

Attributes:

Name	Type	Description
`design`	`DataFrame`	All runs: the factors, the response `y` and `stage` (`'initial'` or `'sequential'`).
`best`	`dict`	The best run (for `goal='minimize'` / `'maximize'`).
`trace`	`DataFrame`	One row per added run: the response, the value of the criterion that selected it and the best response so far.
`fit`	`GPResult`	The Gaussian process fitted to all runs, on factors scaled to [0, 1]. Use `result.predict(newdata)` for predictions in the units of the factors.
`goal`	`str`
`model_info`	`dict`

Examples:

>>> import statspai as sp
>>> f = lambda d: (d["x"] - 0.3) ** 2
>>> res = sp.sequential_design(f, {"x": (0, 1)}, n_new=6, seed=1)
>>> abs(res.best["x"] - 0.3) < 0.05
True

predict ¶

predict(newdata: DataFrame, **kwargs: Any) -> DataFrame

Surrogate prediction (mean, sd, lower, upper) at new points given in the units of the factors.

plot ¶

plot(ax: Any = None, **kwargs: Any) -> Any

Best response so far against the number of runs.

DataSplitResult `dataclass` ¶

Bases: ResultProtocolMixin

A split of a data set into two parts with the same distribution.

Unpacks as train, test = sp.split_data(df).

Attributes:

Name	Type	Description
`train, test`	`DataFrame`
`test_index`	`ndarray`	Positions (0-based) of the test rows in the data.
`energy`	`float`	Energy distance of the test set to the full data, up to a constant; lower is more representative.
`energy_random`	`float`	The same for random splits of that size, averaged.
`model_info`	`dict`

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=200)})
>>> df["y"] = df["x"] ** 2 + rng.normal(size=200)
>>> train, test = sp.split_data(df, test_size=0.2, seed=1)
>>> train.shape[0], test.shape[0]
(160, 40)

SupportPointsResult `dataclass` ¶

Bases: ResultProtocolMixin

Points that represent a distribution.

Attributes:

Name	Type	Description
`points`	`DataFrame`	The `n` points, in the units of the data.
`index`	`array or None`	With `subsample=True`, the rows of the data that were picked.
`energy`	`float`	Energy distance between the points and the sample, up to the term that involves the sample alone (so it can be negative); lower is better.
`energy_random`	`float`	The same for random subsamples of that size, averaged.
`model_info`	`dict`

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame(rng.normal(size=(2000, 2)), columns=["a", "b"])
>>> rep = sp.support_points(df, 20, seed=1)
>>> rep.points.shape
(20, 2)
>>> bool(rep.energy < rep.energy_random)
True

design_criteria ¶

design_criteria(design: Any, bounds: Optional[Mapping[str, Tuple[float, float]]] = None, delta: float = 0.0, r: Optional[float] = None, target: Any = None) -> Series

How well a design fills the experimental region.

Parameters:

Name	Type	Description	Default
`design`	`(DataFrame, array or DesignResult)`	One row per run. Arrays and frames are taken to be on the unit cube unless `bounds` is given.	required
`bounds`	`dict`	`{factor: (lower, upper)}`; the design is scaled to the unit cube with them before anything is computed.	`None`
`delta`	`float`	Added to each squared difference in the MaxPro criterion. With the default a design that repeats a level of a factor has an infinite criterion; a small positive value makes designs with repeated levels comparable.	`0`
`r`	`float`	Power of the reciprocal distance. Default: twice the number of factors.	`None`
`target`	`DataFrame or array`	A sample from the distribution the runs should represent. When given, `energy_distance` to that sample is added. This one is computed in the units of the sample, not on the unit cube.	`None`

Returns:

Type Description

Series

maximin (smallest distance between two runs, larger is better), reciprocal_distance (its smooth stand-in, smaller is better), maxpro (smaller is better; penalises runs that are close in any projection), wraparound_discrepancy and centered_discrepancy (L2 discrepancies, smaller is more uniform), min_projected_distance (smallest gap between two runs along a single factor), fill_distance (largest distance from a point of the region to its nearest run, smaller is better), max_abs_correlation.

Notes

maximin, reciprocal_distance, maxpro and wraparound_discrepancy reproduce maximin.crit, maxpro.crit and uniform.crit of the R package SFDesign; energy_distance reproduces twinning::energy.

fill_distance is evaluated on a finite set of points and is therefore a lower bound of the true value.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> grid = np.array([[0.25, 0.25], [0.25, 0.75], [0.75, 0.25], [0.75, 0.75]])
>>> crit = sp.design_criteria(grid)
>>> float(crit["maximin"])
0.5
>>> bool(np.isinf(crit["maxpro"]))  # two runs share a level
True

References

johnson1990minimax; morris1995exploratory; joseph2015maximum; hickernell1998generalized

factorial_effects ¶

factorial_effects(data: DataFrame, y: str, factors: Optional[Sequence[str]] = None, order: Optional[int] = None, alpha: float = 0.05, reference: str = 'simulated', n_sim: int = 20000, seed: Optional[int] = 0) -> FactorialEffectsResult

Main effects and interactions of a two-level factorial experiment.

Estimates every effect up to a chosen order and tells which of them stand out, including in a design with no replication, where a regression has no residual degrees of freedom and reports no standard errors.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`	One row per run.	required
`y`	`str`	The response.	required
`factors`	`list of str`	The factor columns. Each must take exactly two values; the smaller (in sort order) is coded -1 and the larger +1. Default: every other column with exactly two values.	`None`
`order`	`int`	Largest interaction fitted. Default: all orders, so that a full factorial is fitted with a saturated model.	`None`
`alpha`	`float`	Level of the margins of error.	`0.05`
`reference`	`('simulated', 't')`	Null distribution of `effect / PSE`. `'simulated'` draws it (all effects null, normal errors), as Lenth's own R package `unrepx` does; `'t'` is the t distribution with `m / 3` degrees of freedom of the 1989 paper, which is conservative: with seven effects its margin of error is 1.6 times the simulated one and flags 2.1% of null effects at a nominal 5%, where the simulated margin flags 5.1%.	`'simulated'`
`n_sim`	`int`	Simulated sets of effects.	`20000`
`seed`	`int`	For the simulation. Fixed by default so that repeated calls agree; `None` is read as 0.	`0`

Returns:

Type	Description
`FactorialEffectsResult`	`effects` (one row per term), `pse`, `me`, `sme`, `active`, `aliases`, `summary()`, `plot()` (half-normal).

Notes

effect is the mean response at the high level minus that at the low level, twice the regression coefficient coef.

Lenth (1989): with m effects, s0 = 1.5 median |effect| and the pseudo standard error is 1.5 times the median of the absolute effects below 2.5 s0. The margin of error (ME) is the PSE times the 1 - alpha quantile of |effect| / PSE under the null, the simultaneous margin (SME) uses the largest of the m ratios instead. lenth_p and lenth_p_simultaneous are the matching tail probabilities. With reference='t' the quantiles are those of t(m / 3) at 1 - alpha / 2 and at gamma = (1 + (1 - alpha)^(1/m)) / 2. The method assumes that only a minority of the effects are real; when half of them or more are, the PSE is inflated and nothing is found.

A term whose column equals (up to sign) that of an earlier term is not estimable separately: it is dropped and listed in aliases under the term that stays. Terms are ordered by interaction order, so in a fraction the lower-order member of each alias set is kept.

When there are residual degrees of freedom (replication, or order below the number of factors) the usual se, t and p are reported next to Lenth's. Centre points are not used: drop those runs before calling.

Examples:

An unreplicated 2^3 experiment:

>>> import statspai as sp
>>> d = sp.factorial_design(3).design
>>> d["y"] = [60, 72, 54, 68, 52, 83, 45, 80]
>>> fit = sp.factorial_effects(d, "y")
>>> round(float(fit.effects.loc["C", "effect"]), 2)
23.0
>>> fit.df_resid
0

References

lenth1989quick; daniel1959use; box2005statistics

design_aberration ¶

design_aberration(design: DataFrame, factors: Optional[Sequence[str]] = None, max_order: Optional[int] = None) -> Series

Generalized word length pattern of any factorial design.

Measures how badly effects are confounded in a design with factors at any number of levels, regular or not: A_j is the total squared aliasing between the constant and all j-factor interaction contrasts. A_1 = A_2 = 0 means every pair of factors is balanced (an orthogonal array of strength two); A_3 then measures the confounding of main effects with two-factor interactions. Between two designs of the same size the one with the smaller pattern, read from A_1 on, is preferred.

Parameters:

Name	Type	Description	Default
`design`	`DataFrame`	One row per run. Each column is treated as a categorical factor with the levels that appear in it.	required
`factors`	`list of str`	Columns to use. Default: all.	`None`
`max_order`	`int`	Largest `j`. Default: `min(number of factors, 6)`.	`None`

Returns:

Type	Description
`Series`	`A1` ... `A<max_order>`, plus `resolution` (the first `j` with `A_j > 0`; NaN when all are zero) and `strength` (`resolution - 1`).

Notes

Computed from the coincidences between runs (Xu and Wu 2001), which costs runs^2 * factors * max_order and never forms the contrast columns. For a regular two-level fraction it is the ordinary word length pattern. Reproduces DoE.base::GWLP in R.

Examples:

>>> import statspai as sp
>>> d = sp.factorial_design(4, generators=["D = ABC"]).design
>>> a = sp.design_aberration(d)
>>> [round(float(a[f"A{j}"]), 10) for j in (1, 2, 3, 4)]
[0.0, 0.0, 0.0, 1.0]
>>> int(a["resolution"])
4

References

xu2001generalized

factorial_design ¶

factorial_design(factors: Any, levels: Any = 2, n_runs: Optional[int] = None, generators: Optional[Sequence[str]] = None, center_points: int = 0, replicates: int = 1, randomize: bool = False, seed: Optional[int] = None, alias_order: int = 2) -> FactorialDesignResult

A full or two-level fractional factorial design.

For an experiment that varies several factors at once: every combination of levels (a full factorial), or a fraction of the combinations chosen so that main effects and low-order interactions stay separable.

Parameters:

Name	Type	Description	Default
`factors`	`int, list of str, or dict`	The number of factors (named `A`, `B`, ...), their names, or `{name: [level values]}`. With a dict the design is written in those values and `levels` is ignored.	required
`levels`	`int or list of int`	Levels per factor.	`2`
`n_runs`	`int`	Size of a fraction. Omitted, the full factorial. Available: two-level factors in a power of two (the minimum-aberration regular fraction) or in 12, 20 or 24 runs (Plackett-Burman); factors that all have 3, 5 or 7 levels in a power of that number (a regular fraction, e.g. four three-level factors in 9 runs); and 18 runs for one two-level and up to seven three-level factors (the orthogonal array L18).	`None`
`generators`	`list of str`	Generators of a two-level fraction, e.g. `['E = ABC', 'F = BCD']` (or `'E = ABC'` when names are longer than one character). The factors on the left are the added ones; the others form the full factorial. A leading minus sign on the right (`'E = -ABC'`) picks the other half.	`None`
`center_points`	`int`	Runs at the centre added to a two-level design with numeric levels; they give a test for curvature and a pure-error variance.	`0`
`replicates`	`int`	Copies of the whole design.	`1`
`randomize`	`bool`	Put the runs in random order (the order in which to run them). A `std_order` column keeps the standard order.	`False`
`seed`	`int`		`None`
`alias_order`	`int`	Largest interaction listed in `aliases`.	`2`

Returns:

Type	Description
`FactorialDesignResult`	`design`, `generators`, `defining_relation`, `resolution`, `word_length_pattern`, `aliases`, `summary()`.

Notes

The minimum-aberration search enumerates every choice of generators. It covers up to 15 factors in 32 runs and 12 in 64 (a few seconds at the upper end); beyond that give generators= from a published table. Designs found this way have the word length pattern of the designs returned by FrF2::FrF2 in R; the generators themselves may differ, as minimum-aberration designs are unique only up to relabelling.

defining_relation and aliases are reported for regular two-level fractions only. Regular three-, five- and seven-level fractions do have a defining relation (in words such as ABC^2), but it is not worked out here and defining_relation is left empty for them, as it is for Plackett-Burman designs and the L18 array, which have none; word_length_pattern then holds the generalized pattern A1, A2, ... of sp.design_aberration and resolution the first non-zero entry. In a Plackett-Burman design every main effect is partially aliased with many two-factor interactions: it screens main effects and cannot separate interactions.

Examples:

A half fraction of five factors in 16 runs:

>>> import statspai as sp
>>> d = sp.factorial_design(5, n_runs=16)
>>> d.generators
['E = ABCD']
>>> d.resolution
5

Seven factors in 8 runs, saturated, every two-factor interaction aliased with a main effect:

>>> d = sp.factorial_design(7, n_runs=8)
>>> d.resolution
3

A 2 x 3 full factorial in the units of the factors:

>>> d = sp.factorial_design({"price": [9, 12], "ad": ["a", "b", "c"]})
>>> d.design.shape
(6, 2)

References

box2005statistics; wu2021experiments; plackett1946design; joseph2025experimental

factor_importance ¶

factor_importance(data: DataFrame, y: str, factors: Optional[Sequence[str]] = None, n_neighbors: Optional[int] = None, n_forward: int = 2, standardize: bool = True, n_mc: Optional[int] = None, seed: Optional[int] = None) -> FactorImportanceResult

Rank and select regressors by total Sobol' index, without a model.

A screen to run before modelling: which columns carry information about the outcome at all, through whatever functional form and whatever interactions. Unlike a correlation it sees non-linear effects and interactions, and unlike a variable importance from a forest it does not depend on a fitted model.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`		required
`y`	`str`	The outcome. Numeric, or a 0 / 1 indicator.	required
`factors`	`list of str`	Candidate columns. Default: every other column. Categorical columns (strings, `category`, bool) are one-hot encoded for the neighbour search and count as one factor; numeric columns enter the distance as they are, ties included (see Notes).	`None`
`n_neighbors`	`int`	Neighbours (the observation included) used for each conditional variance. Default 2, or 3 when the outcome takes two values.	`None`
`n_forward`	`int`	Passes of forward selection. In a pass, the factor that lowers the conditional variance most is added, and candidates that do not lower it are dropped from that pass. A second pass gives the dropped ones another chance once more factors are in. Set it to the number of factors for a complete forward selection.	`2`
`standardize`	`bool`	Scale numeric factors to unit variance before distances are measured.	`True`
`n_mc`	`int`	Evaluate the conditional variances at a random subset of this many observations (neighbours are still sought among all). Default: all observations.	`None`
`seed`	`int`	For the subset when `n_mc` is given.	`None`

Returns:

Type	Description
`FactorImportanceResult`	`importance`, `selected`, `noise_share`, `summary()`, `plot()`.

Notes

With v(S) the nearest-neighbour estimate of E[Var(y | X_S)] and S the selected set after forward selection and backward elimination, the importance of a selected factor i is (v(S without i) - v(S)) / (Var(y) - v(S)): the share of the explained variance that is lost without it. Factors outside S get zero. Without subsampling the result equals that of the R package first.

One deliberate difference from first: for a set made of categorical factors only, all observations of a cell are tied nearest neighbours, and which k of them a tree search returns is arbitrary. Here the conditional variance is then the variance within the cell. On data with y = 2 a + 3 [g = v] + noise (the share of a in the explained variance is 0.67) this gives 0.62 for a where first gives 0.55.

The same rule covers numeric factors with repeated values (the coded columns of a factorial design, integer levels): the neighbourhood of an observation is its k nearest neighbours and every observation tied with the k-th, so exact duplicates of a row form one cell and the result does not depend on the order of the rows or on whether levels are coded as numbers or as labels. Data without ties are unaffected.

The selection is greedy. Two factors that matter only through their product, with no effect of either on its own, can be missed: in the population neither lowers the conditional variance when added alone, so whether the first of them enters is decided by sampling noise.

Read the numbers as a ranking and a selection. With correlated factors a total index is small for each of two factors that carry the same information, since either can stand in for the other. Nearest neighbours degrade as the number of selected factors grows; with more than about ten relevant factors the conditional variances are biased upward and weak factors are missed. The measure is about prediction: a factor with a large index is not thereby a cause of the outcome.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> df = pd.DataFrame(rng.uniform(-1, 1, size=(800, 5)),
...                   columns=["x1", "x2", "x3", "x4", "x5"])
>>> df["y"] = (df["x1"] * (1 + df["x2"]) + df["x3"] ** 2
...            + 0.05 * rng.normal(size=800))
>>> res = sp.factor_importance(df, "y")
>>> sorted(res.selected)
['x1', 'x2', 'x3']

References

huang2025factor; sobol2001global

mixture_design ¶

mixture_design(components: Any, kind: str = 'simplex_lattice', degree: int = 2, n: Optional[int] = None, lower: Optional[Sequence[float]] = None, constraint: Optional[Callable[..., Any]] = None, total: float = 1.0, seed: Optional[int] = None) -> DesignResult

A design for an experiment whose factors are shares that sum to one.

Parameters:

Name	Type	Description	Default
`components`	`int or list of str`	The number of components or their names.	required
`kind`	`('simplex_lattice', 'simplex_centroid', 'space_filling')`	`'simplex_lattice'` Every blend whose shares are multiples of `1 / degree` (Scheffe 1958); supports a polynomial of that degree. `'simplex_centroid'` The pure components, the equal blends of every two, of every three, ... up to all (Scheffe 1963); `2^q - 1` runs. `'space_filling'` `n` runs spread evenly over the simplex, or over the part of it that satisfies `constraint`: the support points of the uniform distribution on that region. The choice when the region is irregular or no polynomial is assumed.	`'simplex_lattice'`
`degree`	`int`	For the lattice.	`2`
`n`	`int`	Number of runs, for `'space_filling'`.	`None`
`lower`	`list of float`	Lower bounds on the shares. The design is built on the sub-simplex they leave (in pseudo-components) and mapped back.	`None`
`constraint`	`callable`	For `'space_filling'`: takes a DataFrame of candidate blends and returns one boolean per row.	`None`
`total`	`float`	What the components sum to (100 for percentages).	`1.0`
`seed`	`int`		`None`

Returns:

Type	Description
`DesignResult`	`design` holds the shares; rows sum to `total`.

Notes

Space-filling criteria in criteria are omitted: they are defined on a cube, not a simplex.

Examples:

>>> import statspai as sp
>>> d = sp.mixture_design(3, degree=2)
>>> d.design.shape
(6, 3)
>>> d.design.sum(axis=1).round(12).unique().tolist()
[1.0]
>>> sp.mixture_design(3, kind="simplex_centroid").n_runs
7

References

scheffe1958experiments; scheffe1963simplex; mak2018support

doe_optimal ¶

doe_optimal(model: str, factors: Optional[Mapping[str, Tuple[float, float]]] = None, n: Optional[int] = None, criterion: str = 'D', candidates: Optional[DataFrame] = None, params: Optional[Mapping[str, float]] = None, prior: Optional[Mapping[str, Any]] = None, family: str = 'gaussian', grid: Optional[int] = None, n_prior: int = 32, n_starts: Optional[int] = None, max_iter: int = 3000, seed: Optional[int] = None) -> ModelDesignResult

The design that estimates a given model most precisely.

For choosing the runs of an experiment when the model to be fitted is known in advance: which prices to test for a demand curve, which doses for a dose-response curve, which attribute combinations for a regression with interactions.

Parameters:

Name	Type	Description	Default
`model`	`str`	A linear model as a formula in the factors, `"x1 + x2 + x1:x2 + I(x1*2)"`; or a nonlinear model with its parameters in braces, `"{a} + exp(-{b} x)"` (the syntax of `sp.nls`).	required
`factors`	`dict`	`{name: (lower, upper)}`. The candidates are a grid over this box; the support points found on it are then moved continuously.	`None`
`n`	`int`	Number of runs of an exact design. Omitted, the approximate design: support points and weights.	`None`
`criterion`	`('D', 'A', 'I')`	For `'I'` the prediction variance is averaged uniformly over the box of `factors`, or over the rows of `candidates`. With a binomial or Poisson response it is the variance of the estimated linear predictor.	`'D'`
`candidates`	`DataFrame`	The admissible runs, one column per factor, instead of `factors`. Use it for an irregular region or for qualitative factors (which a formula handles with `C(name)`).	`None`
`params`	`dict`	Values of the parameters of a nonlinear model at which the design is optimal (a locally optimal design).	`None`
`prior`	`dict`	`{parameter: (lower, upper)}` for a uniform prior, or a frozen `scipy.stats` distribution. The criterion is averaged over the prior (a pseudo-Bayesian design), which protects against a wrong guess of the parameters. Parameters left out are fixed at their value in `params`.	`None`
`family`	`('gaussian', 'binomial', 'poisson')`	For a binomial (logit) or Poisson (log) response `model` is the linear predictor, written with braces.	`'gaussian'`
`grid`	`int`	Grid points per factor. Default 201, 41, 15 for one, two, three factors; 4,096 quasi-random points plus the corners for four to ten factors. Beyond ten, pass `candidates=`.	`None`
`n_prior`	`int`	Quasi-random draws from the prior.	`32`
`n_starts`	`int`	Random starts of the exchange algorithm for an exact design. Default: between 2 and 20, fewer for large problems.	`None`
`max_iter`	`int`	Iterations of the weight algorithm.	`3000`
`seed`	`int`		`None`

Returns:

Type	Description
`ModelDesignResult`	`design`, `criterion_value`, `efficiency`, `information`, `summary()`, `to_frame()`.

Notes

Approximate designs are computed with the multiplicative algorithm on the candidate set (monotone for D; Yu 2010), then the support is merged and moved continuously off the grid. efficiency is the equivalence-theorem bound k / max_x d(x) for D, with d the standardised prediction variance and k the number of parameters (Kiefer and Wolfowitz 1960), and 2 - max_x phi(x) / trace for A and I. It is a statement about the candidate set.

Exact designs are found by an exchange algorithm (one run at a time is replaced by the best candidate; at a local optimum a few runs are replaced at random and the search resumes) from several random starts. That gives a good design without a guarantee: the problem has many local optima of nearly equal value. efficiency compares the design with the approximate optimum, an upper bound on what n runs can achieve. On the examples of AlgDesign::optFederov in R the criterion values coincide.

The information matrix of a nonlinear model depends on the unknown parameters. A locally optimal design is only as good as the guess; it typically has as many support points as parameters and leaves no way to check the model. Averaging over a prior spreads the runs.

An optimal design is optimal for the model given. It does not guard against that model being wrong: with real doubt about the functional form, add space-filling runs (sp.design_augment).

Examples:

A quadratic in one factor: a third of the runs at each end and at the centre.

>>> import statspai as sp
>>> d = sp.doe_optimal("x + I(x**2)", {"x": (0, 2)})
>>> d.design.round(3).values.tolist()
[[0.0, 0.333], [1.0, 0.333], [2.0, 0.333]]

Twelve runs for a two-factor model with interaction:

>>> d = sp.doe_optimal("a * b", {"a": (0, 1), "b": (0, 1)}, n=12, seed=1)
>>> d.design["n"].tolist()
[3, 3, 3, 3]

Exponential decay with a guess of the rate: observe at 0 and at 1 / rate.

>>> d = sp.doe_optimal("{a} * exp(-{b} * t)", {"t": (0, 10)},
...                    params={"a": 1, "b": 0.5})
>>> d.design["t"].round(2).tolist()
[0.0, 2.0]

References

kiefer1960equivalence; fedorov1972theory; yu2010monotonic; chaloner1995bayesian

morris_screening ¶

morris_screening(func: Optional[Callable[..., Any]], factors: Any, r: int = 10, levels: int = 4, grid_jump: Optional[int] = None, pass_as: str = 'frame', seed: Optional[int] = None, design: Any = None, y: Any = None) -> MorrisResult

Screen the inputs of a model with one-factor-at-a-time trajectories.

A cheap first pass over a model with many inputs: r (p + 1) runs are enough to sort the inputs into those that do nothing, those with a roughly linear effect, and those whose effect is non-linear or depends on other inputs. Sobol' indices (sp.sobol_indices) need hundreds of runs per input for the same question.

Parameters:

Name	Type	Description	Default
`func`	`callable or None`	The model; see `sp.sobol_indices` for what it receives. `None` when `design` and `y` are given.	required
`factors`	`int, list of str, or dict`	The inputs: their number, their names, or `{name: (lower, upper)}`. Effects are measured per full range of an input.	required
`r`	`int`	Number of trajectories. Each takes `p + 1` runs.	`10`
`levels`	`int`	Grid points per input.	`4`
`grid_jump`	`int`	Grid steps moved at a time. Default `levels // 2`; with an even number of levels every level is then visited equally often (Morris 1991). With an odd number the middle levels are visited more often.	`None`
`pass_as`	`('frame', 'array', 'rows')`		`'frame'`
`seed`	`int`		`None`
`design`	`DataFrame or array`	Runs made already, in the units of the inputs: `p + 1` consecutive rows per trajectory, consecutive rows differing in exactly one input.	`None`
`y`	`array`	The outputs of `design`.	`None`

Returns:

Type	Description
`MorrisResult`	`effects` (`mu`, `mu_star`, `sigma`, `share`), `elementary_effects`, `design`, `summary()`, `plot()`.

Notes

The elementary effect of input i on a trajectory is the change in output over the change in that input (on the unit scale) for the step that moves it. mu_star is the mean absolute effect (Campolongo, Cariboni and Saltelli 2007), sigma the standard deviation over trajectories.

With few trajectories the ranking of unimportant inputs is noisy; what is reliable is the separation of the clearly important from the rest. Given the same design and outputs the statistics equal those of sensitivity::morris in R.

Examples:

>>> import statspai as sp
>>> f = lambda d: 5 * d["a"] + d["b"] ** 2 + 0 * d["c"]
>>> res = sp.morris_screening(f, ["a", "b", "c"], r=8, seed=1)
>>> res.effects["mu_star"].idxmax(), float(res.effects.loc["c", "mu_star"])
('a', 0.0)
>>> res.n_evaluations
32

References

morris1991factorial; campolongo2007effective

sobol_indices ¶

sobol_indices(func: Callable[..., Any], factors: Any, n: int = 1024, estimator: str = 'saltelli', sampling: str = 'sobol', pass_as: str = 'frame', n_boot: int = 200, level: float = 0.95, seed: Optional[int] = None, A: Any = None, B: Any = None) -> SobolResult

First-order and total Sobol' indices of a model's inputs.

Tells which inputs (parameters, calibrated values, assumptions) the output of a model is sensitive to over the whole range of the inputs, interactions included. Use it on a structural or simulation model before estimation to see which parameters the moments can identify, after estimation to report which assumptions drive a counterfactual, or on a fitted flexible model (sp.gp_regress, a forest) to rank the regressors.

Parameters:

Name	Type	Description	Default
`func`	`callable`	The model. By default it receives a DataFrame with one row per run and one column per input and returns one number per run.	required
`factors`	`int, list of str, or dict`	The inputs: their number (uniform on [0, 1]), their names, or `{name: (lower, upper)}` for uniform inputs and `{name: scipy.stats frozen distribution}` for others. Inputs are independent.	required
`n`	`int`	Base sample size. The model is evaluated `n (p + 2)` times. With `sampling='sobol'` it is rounded up to a power of two.	`1024`
`estimator`	`('saltelli', 'jansen')`	Estimator of the first-order index: that of Saltelli et al. (2010) on centred outputs, or Jansen (1999); see Notes for which to prefer. The total index is Jansen's in both cases.	`'saltelli'`
`sampling`	`('sobol', 'random')`	Scrambled Sobol' points (smaller error for a given `n`) or independent draws.	`'sobol'`
`pass_as`	`('frame', 'array', 'rows')`	What `func` receives: a DataFrame, a 2-D array, or one 1-D array per call (slow; for a function of a single run).	`'frame'`
`n_boot`	`int`	Bootstrap replications for the intervals; 0 for none.	`200`
`level`	`float`		`0.95`
`seed`	`int`		`None`
`A`	`DataFrame or array`	Two independent samples of the inputs to use instead of drawing them, in the units of the inputs (then `n`, `sampling` and any distributions in `factors` are not used).	`None`
`B`	`DataFrame or array`	Two independent samples of the inputs to use instead of drawing them, in the units of the inputs (then `n`, `sampling` and any distributions in `factors` are not used).	`None`

Returns:

Type	Description
`SobolResult`	`indices` (`first`, `total` and limits per input), `variance`, `summary()`, `plot()`.

Notes

With A and B two samples and AB_i the sample A with column i taken from B,

total_i = mean (f(A) - f(AB_i))^2 / (2 V) and first_i = mean (f(B) - m) (f(AB_i) - f(A)) / V (Saltelli et al., with m and V the mean and variance of the pooled outputs), or first_i = 1 - mean (f(B) - f(AB_i))^2 / (2 V) (Jansen, with V the variance of f(A)).

Neither first-order estimator dominates. On the eight-input borehole function with n = 512 Sobol' points the root mean squared errors are 0.008 (Saltelli) against 0.005 (Jansen) for the input that explains 83% of the variance, 0.005 against 0.011 for inputs that explain 4%, and 0.0001 against 0.009 for the three inert inputs. The default is the one with the smaller worst case, and the one that tells an inert input from a weak one. With sampling='random' Jansen's is three times more accurate for a dominant input.

estimator='jansen' reproduces sensitivity::soboljansen in R, which computes the same quantities with the sums divided by 2 n - 1 instead of 2 n; after that rescaling the two agree to rounding error.

Estimates can fall slightly below zero or sum to more than one by sampling error; that is the Monte Carlo noise to read the indices against. The intervals are percentile bootstrap intervals over the base sample; with Sobol' points they are conservative, since the bootstrap treats the points as independent draws.

The decomposition assumes independent inputs. With dependent inputs the indices no longer add up to shares of the variance.

Examples:

y = x1 + 2 x2 + x1 x3 with independent uniform inputs:

>>> import statspai as sp
>>> f = lambda d: d["x1"] + 2 * d["x2"] + d["x1"] * d["x3"]
>>> res = sp.sobol_indices(f, ["x1", "x2", "x3"], n=2048, seed=1, n_boot=0)
>>> bool(res.indices.loc["x2", "first"] > res.indices.loc["x1", "first"])
True
>>> bool(res.indices.loc["x3", "total"] > res.indices.loc["x3", "first"])
True

References

sobol2001global; jansen1999analysis; saltelli2010variance; harenberg2019uncertainty

sequential_design ¶

sequential_design(func: Callable[..., Any], factors: Any, n_new: int = 20, goal: str = 'minimize', n_init: Optional[int] = None, design: Optional[DataFrame] = None, y: Any = None, kernel: str = 'rbf', noisy: bool = False, n_candidates: Optional[int] = None, pass_as: str = 'frame', seed: Optional[int] = None, criterion: str = 'variance') -> SequentialDesignResult

Choose the runs of an expensive function one at a time.

Bayesian optimisation and active learning in one function: minimise (or maximise) a function that is costly to evaluate with as few evaluations as possible, or learn it everywhere to build a fast surrogate. Typical uses are the objective of a simulated-moments or indirect-inference estimator, a tuning parameter chosen by an expensive cross-validation, or a simulation model that is to be replaced by an emulator.

Parameters:

Name	Type	Description	Default
`func`	`callable`	The function. By default it receives a DataFrame with one row per run and returns one number per run.	required
`factors`	`dict`	`{name: (lower, upper)}`.	required
`n_new`	`int`	Runs to add after the initial design.	`20`
`goal`	`('minimize', 'maximize', 'emulate')`	`'minimize'` / `'maximize'` add the run with the largest expected improvement (Jones, Schonlau and Welch 1998); `'emulate'` the run that most reduces the uncertainty of the surrogate over the region (see `criterion`).	`'minimize'`
`n_init`	`int`	Size of the initial maximum projection design. Default `max(5 p, 6)`.	`None`
`design`	`DataFrame and array`	Runs already made and their responses, used instead of an initial design.	`None`
`y`	`DataFrame and array`	Runs already made and their responses, used instead of an initial design.	`None`
`kernel`	`('rbf', 'matern52', 'matern32')`		`'rbf'`
`noisy`	`bool`	`False`: the function returns the same value when called twice at the same point, and the surrogate interpolates. `True`: the output has noise (a simulation with fresh random numbers), a noise variance is estimated and the improvement is measured against the best predicted value at the runs.	`False`
`n_candidates`	`int`	Candidate points scored per step. Default `min(20000, 2000 p)`.	`None`
`pass_as`	`('frame', 'array', 'rows')`	See `sp.sobol_indices`.	`'frame'`
`seed`	`int`		`None`
`criterion`	`('variance', 'alc')`	For `goal='emulate'`. `'variance'` adds the run where the predictive variance is largest; `'alc'` the run that lowers the predictive variance averaged over the region the most. Neither is better throughout: with 10 initial and 20 added runs the root mean squared prediction error on the two-input Branin function was 0.18 (variance) against 0.36 (alc), and on a three-input test function 0.19 against 0.05, the first figure driven by one poor run in four.	`'variance'`

Returns:

Type	Description
`SequentialDesignResult`	`design`, `best`, `trace`, `fit`, `predict()`, `summary()`, `plot()`.

Notes

Each step refits the Gaussian process (sp.gp_regress), scores a fresh quasi-random candidate set, plus points near the best runs when optimising, and evaluates the function at the winner. The candidates are finite, so the optimum is located up to their spacing; finish with a local optimiser started from best when more digits are needed.

Expected improvement balances looking near good values with looking where little is known. It finds the neighbourhood of a global optimum in few runs; it is not a convergence proof. A Gaussian process is a poor surrogate for a function with jumps or with more than about ten inputs.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> f = lambda d: np.sin(3 * d["a"]) + (d["b"] - 0.5) ** 2
>>> res = sp.sequential_design(f, {"a": (0, 3), "b": (0, 1)},
...                            n_new=12, seed=1)
>>> bool(res.best["y"] < -0.98)
True

References

jones1998efficient; sacks1989design; joseph2025experimental

design_augment ¶

design_augment(design: Any, n_new: int, candidates: Any = None, criterion: str = 'maxpro', bounds: Optional[Mapping[str, Tuple[float, float]]] = None, constraint: Optional[Callable[..., Any]] = None, n_candidates: Optional[int] = None, delta: float = 0.0, seed: Optional[int] = None) -> DesignResult

Add runs to an existing design so that the whole stays space-filling.

For a second batch after the first has been run, for validation runs that keep away from the training runs, or for a design on a region of irregular shape (pass only feasible candidates).

Parameters:

Name	Type	Description	Default
`design`	`(DataFrame, array or DesignResult)`	The runs already made. They are kept as they are.	required
`n_new`	`int`	Runs to add.	required
`candidates`	`DataFrame or array`	The points the new runs are chosen from, in the units of the design and inside its region (the unit cube, the bounds of a DesignResult, or `bounds`); points outside are an error. Default: a scrambled Sobol' set over the region.	`None`
`criterion`	`('maxpro', 'maximin')`	`'maxpro'` adds, each time, the candidate with the smallest sum over existing runs of `1 / prod_l (x_l - c_l)^2`; `'maximin'` the candidate farthest from its nearest run.	`'maxpro'`
`bounds`	`dict`	`{factor: (lower, upper)}`. Needed when the design is not on the unit cube and is not a DesignResult.	`None`
`constraint`	`callable`	Filters the default candidates; see `sp.space_filling`.	`None`
`n_candidates`	`int`	Size of the default candidate set. Default `max(8192, 400 (n + n_new))` rounded up to a power of two.	`None`
`delta`	`float`	Added to squared differences in the MaxPro criterion. With 0 a candidate that repeats a level of an existing run is never chosen.	`0`
`seed`	`int`	For the default candidate set.	`None`

Returns:

Type	Description
`DesignResult`	The augmented design, old runs first. `model_info['new_rows']` lists the positions of the added runs, `model_info ['candidate_index']` their rows in `candidates`.

Notes

The search is greedy: each run is the best addition given those before it, which is what makes the result usable in batches. Given the same candidates the 'maxpro' choice is the one made by MaxPro::MaxProAugment in R.

Examples:

>>> import statspai as sp
>>> first = sp.space_filling(10, 2, seed=1)
>>> both = sp.design_augment(first, 10, seed=2)
>>> both.design.shape
(20, 2)
>>> bool((both.design.iloc[:10].to_numpy() == first.design.to_numpy()).all())
True

References

joseph2015maximum; johnson1990minimax

space_filling ¶

space_filling(n: int, factors: Any, method: str = 'maxpro', seed: Optional[int] = None, n_starts: int = 3, iterations: Optional[int] = None, polish: Optional[bool] = None, constraint: Optional[Callable[..., Any]] = None, n_candidates: Optional[int] = None, delta: float = 0.0, r: Optional[float] = None, qualitative: Optional[Mapping[str, Any]] = None) -> DesignResult

A space-filling design: runs spread evenly over the factor region.

The design to use when little is known about the response surface: a simulation study over a parameter region, the training runs of a surrogate model, starting values for a multi-start optimiser, or the draws of a simulated-moments estimator.

Parameters:

Name	Type	Description	Default
`n`	`int`	Number of runs.	required
`factors`	`int, list of str, or dict`	The number of factors, their names, or `{name: (lower, upper)}`. Without bounds the design is on the unit cube.	required
`method`	`str`	`'maxpro'` Maximum projection design (Joseph, Gul and Ba 2015): fills the full region and every lower-dimensional projection. The default, because usually only some factors matter. `'maximin'` Latin hypercube with the largest minimum distance between runs (Morris and Mitchell 1995). `'uniform'` Latin hypercube with the smallest wrap-around discrepancy; the one to use for numerical integration. `'lhs'` A random Latin hypercube (McKay, Beckman and Conover 1979). `'sobol'`, `'halton'` Scrambled low-discrepancy sequences. Extensible: the first `m` runs of a longer sequence are a sequence themselves. `'random'` Independent uniform draws, the baseline.	`'maxpro'`
`seed`	`int`		`None`
`n_starts`	`int`	Independent searches; the best design is returned.	`3`
`iterations`	`int`	Swaps tried per search. Default `min(400000, max(20000, 1000 n p))`.	`None`
`polish`	`bool`	Move the levels continuously after the search. Default `True` for `'maxpro'` (the design is then no longer on the Latin hypercube midpoints but still has `n` distinct levels per factor) and `False` otherwise.	`None`
`constraint`	`callable`	Takes a DataFrame of candidate runs in the units of the factors and returns one boolean per row (or takes one row and returns a boolean). When given, the design is built by adding feasible candidates one at a time (`sp.design_augment`), since a Latin hypercube does not exist on an irregular region; `method` must be `'maxpro'` or `'maximin'`.	`None`
`n_candidates`	`int`	Size of the candidate set under `constraint`.	`None`
`delta`	`float`	Added to squared differences in the MaxPro criterion.	`0`
`r`	`float`	Power of the reciprocal-distance criterion that stands in for maximin during the search. Default `2 p`.	`None`
`qualitative`	`dict`	`{name: [levels]}` for factors without an order (a model variant, a solver, a region). Each level gets the same number of runs, as nearly as `n` allows: when `n` is a multiple of the number of level combinations every combination occurs equally often, otherwise the counts of the levels of each factor differ by at most one. The quantitative factors are arranged so that the runs fill the space within every level as well as overall (Joseph, Gul and Ba 2020). `method` must be `'maxpro'`. `criteria['maxpro_qq']` is the criterion that was minimised; the other measures refer to the quantitative factors alone.	`None`

Returns:

Type	Description
`DesignResult`	`design` (the runs as a DataFrame), `criteria`, `summary()`, `plot()`.

Notes

The search is stochastic and finds a good design, not a certified optimum. On the criterion value the designs are on par with those of the R packages SFDesign and MaxPro.

Examples:

>>> import statspai as sp
>>> d = sp.space_filling(20, {"beta": (0.9, 0.99), "sigma": (1, 5)}, seed=1)
>>> d.design.shape
(20, 2)
>>> bool(d.design["beta"].between(0.9, 0.99).all())
True
>>> lhs = sp.space_filling(20, 2, method="lhs", seed=1)
>>> bool(d.criteria["maxpro"] < lhs.criteria["maxpro"])
True

References

joseph2015maximum; joseph2020designing; morris1995exploratory; mckay1979comparison; joseph2025experimental

split_data ¶

split_data(data: DataFrame, test_size: float = 0.2, columns: Optional[Sequence[str]] = None, method: str = 'auto', start: Optional[int] = None, seed: Optional[int] = None, max_iter: int = 500, tol: float = 0.0001) -> DataSplitResult

Split data into a training and a test set that resemble each other.

A random split leaves the test set unrepresentative by chance: a region of the covariates, or a tail of the outcome, can be missing from it, and the measured test error then depends on the luck of the draw. This split is chosen so that the test set has, as nearly as possible, the joint distribution of the whole data set, outcome included.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`	One row per observation, outcome and features together.	required
`test_size`	`float`	Share of the rows in the test set.	`0.2`
`columns`	`list of str`	Columns that define the distribution. Default: all. Categorical columns are expanded into Helmert contrasts.	`None`
`method`	`('auto', 'support', 'twinning')`	`'support'` (SPlit, Joseph and Vakayil 2022) computes support points of the data and takes the nearest rows. `'twinning'` (Vakayil and Joseph 2022) walks through the data taking one row out of each group of `r = 1 / test_size` neighbours; it is far faster on large data and needs `1 / test_size` to be a whole number. `'auto'` uses twinning when it applies and there are more than 2,000 rows, otherwise support points. Support points need more distinct rows than test rows; on data with fewer (all columns discrete, many repeated rows) `'auto'` uses twinning, which handles ties, and says so in `model_info['notes']`, while an explicit `'support'` is an error.	`'auto'`
`start`	`int`	Twinning only: the row (position) where the walk starts. With it the split is deterministic.	`None`
`seed`	`int`		`None`
`max_iter`	`int`	Passed to the support-points iteration.	`500`
`tol`	`int`	Passed to the support-points iteration.	`500`

Returns:

Type	Description
`DataSplitResult`	`train`, `test`, `test_index`, `energy`, `energy_random`; unpacks into `train, test`.

Notes

The split does not depend on any model. It is meant for assessing prediction error of a model fitted once. Do not use it where the independence of the two parts is what a procedure relies on (sample splitting for honest inference, cross-fitting): there the split must be random.

With method='twinning' and the same start the test rows are those returned by twinning::twin in R.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(3)
>>> df = pd.DataFrame({"x": rng.normal(size=300)})
>>> df["y"] = np.sin(2 * df["x"]) + 0.3 * rng.normal(size=300)
>>> split = sp.split_data(df, test_size=0.2, seed=1)
>>> split.test.shape
(60, 2)
>>> bool(split.energy < split.energy_random)
True

References

joseph2022split; vakayil2022data; mak2018support

support_points ¶

support_points(data: Any, n: int, columns: Optional[Sequence[str]] = None, weights: Any = None, subsample: bool = False, standardize: bool = True, seed: Optional[int] = None, n_starts: int = 1, max_iter: int = 500, tol: float = 0.0001, max_sample: int = 50000) -> SupportPointsResult

The n points that best represent a sample or a distribution.

A compact stand-in for a distribution: where a random subsample of n points reproduces expectations with error of order n^(-1/2), support points do markedly better, so far fewer model evaluations are needed to propagate input uncertainty, and a thinned posterior sample keeps the shape of the full one.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame, array, or dict of distributions`	A sample from the target distribution (data, MCMC draws, draws of model inputs), one row per draw. Or `{name: scipy.stats frozen distribution}` for independent inputs with known distributions; a large quasi-random sample is then drawn internally.	required
`n`	`int`	Number of points.	required
`columns`	`list of str`	Columns to use. Default: all, which must be numeric.	`None`
`weights`	`str or array`	Weights of the draws (importance weights, survey weights).	`None`
`subsample`	`bool`	Return rows of the data: each support point is replaced by the nearest row not yet taken. Use it when the points must be real observations.	`False`
`standardize`	`bool`	Scale every column by its standard deviation (under `weights`, that of the weighted distribution) before measuring distances. Energy distance is not invariant to the units.	`True`
`seed`	`int`		`None`
`n_starts`	`int`	Random starts; the lowest energy is kept.	`1`
`max_iter`	`int`		`500`
`tol`	`float`	Stop when no coordinate (in standard deviations) moves by more than `tol * sqrt(dimension)`.	`1e-4`
`max_sample`	`int`	Larger samples are thinned at random to this size for the optimisation; the energy is still evaluated on all rows.	`50000`

Returns:

Type	Description
`SupportPointsResult`	`points`, `index`, `energy`, `energy_random`, `summary()`.

Notes

The points minimise 2 mean|x_i - Y| - mean|x_i - x_j| over the sample Y by the convex-concave iteration of Mak and Joseph (2018). The problem is not convex, so different starts give different, almost equally good, point sets. On the energy criterion the result is on par with support::sp in R.

An average of a function over the points estimates its expectation. The points carry equal weight.

Examples:

Twenty points for a bivariate normal with correlation 0.5, from a sample:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> z = rng.multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]], size=5000)
>>> rep = sp.support_points(pd.DataFrame(z, columns=["x1", "x2"]), 20, seed=1)
>>> rep.points.shape
(20, 2)

From known input distributions:

>>> from scipy import stats
>>> rep = sp.support_points({"r": stats.norm(0.1, 0.016),
...                          "k": stats.uniform(990, 120)}, 30, seed=1)
>>> list(rep.points.columns)
['r', 'k']

References

mak2018support; szekely2013energy

statspai.doe¶

doe ¶

DesignResult dataclass ¶

to_frame ¶

plot ¶

FactorialEffectsResult dataclass ¶

active property ¶

plot ¶

FactorialDesignResult dataclass ¶

to_frame ¶

FactorImportanceResult dataclass ¶

plot ¶

ModelDesignResult dataclass ¶

to_frame ¶

MorrisResult dataclass ¶

plot ¶

SobolResult dataclass ¶

plot ¶

SequentialDesignResult dataclass ¶

predict ¶

plot ¶

DataSplitResult dataclass ¶

SupportPointsResult dataclass ¶

design_criteria ¶

factorial_effects ¶

design_aberration ¶

factorial_design ¶

factor_importance ¶

mixture_design ¶

doe_optimal ¶

morris_screening ¶

sobol_indices ¶

sequential_design ¶

design_augment ¶

space_filling ¶

split_data ¶

support_points ¶

`statspai.doe`¶

DesignResult `dataclass` ¶

FactorialEffectsResult `dataclass` ¶

active `property` ¶

FactorialDesignResult `dataclass` ¶

FactorImportanceResult `dataclass` ¶

ModelDesignResult `dataclass` ¶

MorrisResult `dataclass` ¶

SobolResult `dataclass` ¶

SequentialDesignResult `dataclass` ¶

DataSplitResult `dataclass` ¶

SupportPointsResult `dataclass` ¶