mixle.doe.analysis module

Analysis of designed experiments: factorial effects and second-order response surfaces.

These turn the runs of a design (see mixle.doe.factorial) and their measured responses into the quantities a practitioner reads off: the effect of each factor and interaction in a two-level design, and – for a response-surface design – the fitted second-order model, its stationary point, and the canonical (eigenvalue) analysis that says whether that point is a maximum, minimum, or saddle.

class FactorialEffects(terms, coef, effects, intercept, residual_std, se, aliases, estimable_contrasts)[source]

Bases: object

Estimated effects from a two-level factorial / fractional-factorial / Plackett-Burman design.

Only estimable quantities are reported as effects. When the design cannot separate two or more terms – classical aliasing/confounding, including a factor that never varies – those individual terms have no unique estimate: their entries in coef / effects are nan. aliases says which terms are aliased with which; estimable_contrasts gives the one well-defined combined quantity the data actually supports for each such group.

Parameters:
terms

term names ("intercept", "x0", "x0:x1", …).

Type:

list[str]

coef

least-squares regression coefficients in coded +/-1 units (nan for a term that is aliased with another and so has no individual estimate).

Type:

numpy.ndarray

effects

the classical effect per term – the change in mean response as a factor moves from its low to its high level, i.e. 2 * coef (the intercept entry is just the grand mean).

Type:

numpy.ndarray

intercept

the grand mean of the response (nan if the intercept itself is aliased, e.g. with a factor that never varies).

Type:

float

residual_std

residual standard deviation when the design has spare runs (else None).

Type:

float | None

se

standard error of each entry of effects. None throughout when the design has no spare degrees of freedom to estimate a residual variance (e.g. a saturated design); individual entries are nan wherever effects is (aliased terms).

Type:

numpy.ndarray | None

aliases

for each term, the other term names it is aliased with (empty list if the term is uniquely estimable).

Type:

dict[str, list[str]]

estimable_contrasts

for each alias group of two or more terms, keyed by its member term names joined by their signed defining relation (for example "x0-x1" for opposite columns), the (effect, se) of the exact combined quantity that is estimable – empty when the design has no aliasing.

Type:

dict[str, tuple[float, float | None]]

as_dict()[source]

Map each non-intercept term to its effect (nan where aliased – see estimable_contrasts).

Return type:

dict[str, float]

factorial_effects(design, y, *, interactions=True, coded=False)[source]

Estimate main effects and two-factor interactions from a two-level design.

Fits the linear model y ~ 1 + x_i (+ x_i x_j) in coded +/-1 units by least squares; the coefficients are half the classical effects. design is the (n, d) run matrix (the real factor levels, coded to +/-1 automatically – or pass coded=True if it is already coded, in which case every entry must be exactly -1.0 or 1.0), y the measured responses. Set interactions=False for a main-effects-only (e.g. screening) fit.

Rejects designs that cannot support the requested model outright: non-finite values, fewer runs than model parameters, and (for coded=True) levels other than +/-1. A design that has enough runs but still cannot separate every term – aliasing, including a factor that never varies – is not rejected: the aliased terms come back as nan in coef / effects, with the one combined quantity the data does support reported in estimable_contrasts (see FactorialEffects). That automatic resolution only covers the classical case where the aliasing is exact term-for-term confounding, as in any standard fractional-factorial / Plackett-Burman alias structure; a design that is rank-deficient for some other reason (e.g. too few runs and not a clean fractional design) raises instead of guessing which contrasts might be estimable.

Parameters:
Return type:

FactorialEffects

class ResponseSurface(coef, terms, b, B, stationary_point, eigenvalues, kind, residual_std, model_rank, n_parameters, degrees_of_freedom, estimable)[source]

Bases: object

A fitted second-order response surface y = b0 + b'x + x'Bx and its canonical analysis.

Parameters:
coef

full coefficient vector (intercept, linears, then the upper-triangular second-order terms).

Type:

numpy.ndarray

terms

matching term names.

Type:

list[str]

b

linear coefficient vector (d,).

Type:

numpy.ndarray

B

symmetric (d, d) matrix of quadratic coefficients (cross terms split onto both halves).

Type:

numpy.ndarray

stationary_point

x* solving the stationarity equation B x* = -b/2, or None when no such point exists – i.e. the gradient b + 2 B x is never zero anywhere, which can only happen when B is singular (kind == "no_stationary_point"; see below).

Type:

numpy.ndarray | None

eigenvalues

eigenvalues of B – all negative => a maximum, all positive => a minimum, mixed signs => a saddle; any eigenvalue ~0 (with a genuine stationary point) => a ridge, a whole line/subspace of equally-stationary points, of which stationary_point is one representative (the minimum-norm) point.

Type:

numpy.ndarray

kind

"maximum" / "minimum" / "saddle" / "ridge" / "no_stationary_point". The last case is distinct from a ridge: a ridge is a genuine (degenerate) stationary subspace, while "no_stationary_point" means the stationarity equation has no solution at all (e.g. the exact surface y = x0**2 + x1, whose gradient’s second component is always exactly 1).

Type:

str

residual_std

residual standard deviation when the design has spare runs (else None).

Type:

float | None

model_rank

rank of the fitted quadratic model matrix.

Type:

int

n_parameters

number of quadratic coefficients requested.

Type:

int

degrees_of_freedom

residual degrees of freedom, n_runs - n_parameters.

Type:

int

estimable

whether every requested coefficient is identified. Construction refuses a non-estimable fit, so returned surfaces always report True.

Type:

bool

predict(x)[source]

Predict the response at points x (m, d) from the fitted surface.

Return type:

ndarray

gradient(x)[source]

The response gradient b + 2 B x at x – its direction is the path of steepest ascent.

Return type:

ndarray

response_surface(x, y)[source]

Fit a full second-order (quadratic) response surface and analyse its stationary point.

Least-squares-fits y = b0 + sum b_i x_i + sum_{i<=j} b_{ij} x_i x_j to the design runs x (n, d) and responses y, then solves the stationarity equation B x* = -b/2 for the stationary point and classifies it from the eigenvalues of the quadratic matrix B. When B is singular a solution only exists if -b/2 lies in its range, which is checked explicitly via the residual of the stationarity equation – a least-squares “solution” that does not actually zero the gradient is reported as kind="no_stationary_point" rather than a fake stationary point (see ResponseSurface). Fit on the coded design for a well-conditioned model (the classic central-composite / Box-Behnken workflow).

Return type:

ResponseSurface

design_diagnostics(design, model, *, ref=None)[source]

Quality diagnostics for a design under a model – “is this a good design to run?”.

Builds the model matrix F = model(design) and reports, relative to a hypothetical perfectly orthogonal design (where each is 1.0):

  • d_efficiency – det(M)**(1/p) / n: overall coefficient-estimation precision;

  • a_efficiency – p / (n * trace(M^-1)): average coefficient variance;

  • g_efficiency – p / (n * max prediction variance) over ref (or the design itself);

  • condition_number of M (large => near-collinear / fragile to fit; None when the requested model is rank-deficient), plus a finite effective_condition_number within the estimable subspace and explicit rank / full_rank fields;

  • max_correlation – the largest absolute pairwise correlation among the non-intercept model columns (the aliasing check; 0 for an orthogonal design).

model is a model-matrix function such as mixle.doe.optimal.polynomial_features(). Use it on the coded design for meaningful efficiencies.

ref – when given – is a reference design: raw candidate points in the same (m, d) raw-factor-column form as design (e.g. a finer grid over the same region, or a wider region to probe extrapolation risk), not a pre-built model matrix. It is passed through the same model used to build the design’s own model matrix before its prediction variance is computed, so g_efficiency compares like with like; a ref whose raw column count does not match design’s raises. g_efficiency is guaranteed to fall on the documented (0, 1] scale only for the self-referential case (ref=None, or ref equal to design itself) – a reference set concentrated in a lower-variance region than the design’s own points (e.g. just its centre) can legitimately score above 1.0.

The raw design and every transformed model matrix must be nonempty, two-dimensional, finite, and preserve the run axis. A rank-zero model is rejected. A partially estimable model reports zero full-model efficiencies, condition_number=None, and its finite effective condition number rather than publishing NaN/Inf quality evidence.

Return type:

dict