Designs

Design generators. Every design is constructed from Factor objects and produces a pandas design matrix from generate_design().

Factorial aliasing and blocking

Source: factorial_aliasing_blocking.dot

Regular two-level fractions and regular factorial blocks both use treatment words and algebra over $\mathrm{GF}(2)$, but they solve different problems. A fractional factorial deliberately omits treatment combinations, so its defining relation determines which treatment effects are aliased. Resolution and the word-length pattern summarize how severe those aliases are, and foldover can add runs when a particular ambiguity matters.

Blocking keeps the full factorial treatment set but partitions it into restricted groups. The independent block generators produce a block-defining subgroup; treatment effects in that subgroup are intentionally confounded with blocks. FactorialDesign.block_structure() exposes that choice explicitly, and run-order randomization then occurs within blocks rather than destroying the restriction.

Base

industrialstats.designs.base

Base class for all experimental designs.

Factor dataclass

Factor(name: str, levels: list[str | float | int], factor_type: str = 'categorical')

Represents an experimental factor.

Attributes:

Name	Type	Description
name	str	Name of the factor as it appears in the design matrix.
levels	list[str | float | int]	Supported discrete levels for the factor.
factor_type	str	Either "categorical" or "continuous".

ExperimentalDesign

ExperimentalDesign(name: str)

Bases: ABC

Abstract base class for experimental designs.

Subclasses must implement :meth:generate_design and :meth:validate_design to provide bespoke construction and integrity checks for the experiment.

Initialise the design container.

Parameters:

Name	Type	Description	Default
name	str	Human-readable label describing the design.	required

Source code in src/industrialstats/designs/base.py

def __init__(self, name: str) -> None:
    """Initialise the design container.

    Parameters
    ----------
    name
        Human-readable label describing the design.
    """

    # Persist metadata and runtime state for the design definition.
    self.name = name
    self.factors: list[Factor] = []
    self.design_matrix: pd.DataFrame | None = None
    self.randomized: bool = False
    self.seed: int | None = None

run_count property

run_count: int

Number of experimental runs currently in the design.

factor_names property

factor_names: list[str]

List of factor names present in the design.

is_balanced property

is_balanced: bool

Indicate whether the design is balanced across categorical factors.

Returns:

Type	Description
bool	True when every combination of categorical factor levels is represented equally often without missing values. If the design has not been generated the property returns False.

Examples:

>>> from industrialstats.designs.factorial import FactorialDesign
>>> factors = [Factor("A", [0, 1]), Factor("B", [0, 1])]
>>> design = FactorialDesign(factors, randomize=False)
>>> design.generate_design()
>>> design.is_balanced
True

design_efficiency property

design_efficiency: dict[str, float]

Compute baseline efficiency metrics for the design.

Returns:

Type	Description
dict[str, float]	Dictionary containing efficiency statistics such as the run fraction, replication factor, balance index, and the proportion of missing data.

Examples:

>>> from industrialstats.designs.factorial import FactorialDesign
>>> design = FactorialDesign([Factor("A", [0, 1])], randomize=False)
>>> design.generate_design()
>>> metrics = design.design_efficiency
>>> sorted(metrics.keys())
['balance_index', 'missing_rate', 'replication_factor', 'run_fraction']

generate_design abstractmethod

generate_design() -> DataFrame

Generate the experimental design matrix.

Source code in src/industrialstats/designs/base.py

@abstractmethod
def generate_design(self) -> pd.DataFrame:
    """Generate the experimental design matrix."""
    pass

validate_design abstractmethod

validate_design() -> bool

Validate the experimental design.

Source code in src/industrialstats/designs/base.py

@abstractmethod
def validate_design(self) -> bool:
    """Validate the experimental design."""
    pass

add_factor

add_factor(factor: Factor) -> None

Add a factor to the design.

Parameters:

Name	Type	Description	Default
factor	Factor	Factor description to register.	required

Raises:

Type	Description
TypeError	If factor is not a :class:Factor instance.

Source code in src/industrialstats/designs/base.py

def add_factor(self, factor: Factor) -> None:
    """Add a factor to the design.

    Parameters
    ----------
    factor
        Factor description to register.

    Raises
    ------
    TypeError
        If ``factor`` is not a :class:`Factor` instance.
    """

    # Enforce a consistent factor representation for downstream logic.
    if not isinstance(factor, Factor):
        raise TypeError("factor must be a Factor instance")

    # Append the validated factor to the ordered factor list.
    self.factors.append(factor)

randomize

randomize(seed: int | None = None) -> None

Randomise the run order of the experiment.

Parameters:

Name	Type	Description	Default
seed	int | None	Optional random seed used to create deterministic shuffles.	None

Raises:

Type	Description
ValueError	If the design matrix has not been generated.

References

Montgomery, D.C. (2017). Design and Analysis of Experiments, 9th ed., Wiley.

Examples:

>>> from industrialstats.designs.factorial import FactorialDesign
>>> design = FactorialDesign({"A": [1, -1], "B": [1, -1]})
>>> design.generate_design()
>>> design.randomize(seed=42)
>>> design.design_matrix[["RunOrder", "A", "B"]].head()
   RunOrder  A  B
0         1  1 -1
1         2 -1  1

Source code in src/industrialstats/designs/base.py

def randomize(self, seed: int | None = None) -> None:
    """Randomise the run order of the experiment.

    Parameters
    ----------
    seed
        Optional random seed used to create deterministic shuffles.

    Raises
    ------
    ValueError
        If the design matrix has not been generated.

    References
    ----------
    Montgomery, D.C. (2017). *Design and Analysis of Experiments*, 9th
    ed., Wiley.

    Examples
    --------
    >>> from industrialstats.designs.factorial import FactorialDesign
    >>> design = FactorialDesign({"A": [1, -1], "B": [1, -1]})
    >>> design.generate_design()
    >>> design.randomize(seed=42)
    >>> design.design_matrix[["RunOrder", "A", "B"]].head()
       RunOrder  A  B
    0         1  1 -1
    1         2 -1  1
    """

    # Guard against missing design matrices prior to shuffling runs.
    if self.design_matrix is None:
        raise ValueError(
            "Design matrix not generated yet. Call generate_design() first."
        )

    # Persist the seed so exported metadata reflects the randomisation state.
    self.seed = seed

    # Shuffle the design matrix using a reproducible random state when provided.
    self.design_matrix = self.design_matrix.sample(
        frac=1, random_state=seed
    ).reset_index(drop=True)

    # Insert the sequential run order column for traceability of randomisation.
    self.design_matrix.insert(0, "RunOrder", range(1, len(self.design_matrix) + 1))
    self.randomized = True

to_csv

to_csv(filename: str) -> None

Export design to a CSV file.

Parameters:

Name	Type	Description	Default
filename	str	Destination file path.	required

Raises:

Type	Description
ValueError	If no design matrix is available.

Source code in src/industrialstats/designs/base.py

def to_csv(self, filename: str) -> None:
    """Export design to a CSV file.

    Parameters
    ----------
    filename
        Destination file path.

    Raises
    ------
    ValueError
        If no design matrix is available.
    """

    # Ensure that the design matrix exists before writing to disk.
    if self.design_matrix is None:
        raise ValueError("No design matrix to export.")

    # Persist the design matrix in a simple comma-separated format.
    self.design_matrix.to_csv(filename, index=False)

to_excel

to_excel(filename: str, include_metadata: bool = True) -> None

Export the design to an Excel workbook with optional metadata.

Parameters:

Name	Type	Description	Default
filename	str	Path to the output workbook (.xlsx or .xlsm).	required
include_metadata	bool	Whether to include a summary worksheet.	True

Raises:

Type	Description
ValueError	If no design matrix exists or the extension is invalid.
OSError	If the workbook cannot be written to disk.

Examples:

>>> from industrialstats.designs.factorial import FactorialDesign
>>> design = FactorialDesign([Factor("A", ["Low", "High"])], randomize=False)
>>> design.generate_design()
>>> design.to_excel("factorial_design.xlsx")

Source code in src/industrialstats/designs/base.py

def to_excel(self, filename: str, include_metadata: bool = True) -> None:
    """Export the design to an Excel workbook with optional metadata.

    Parameters
    ----------
    filename
        Path to the output workbook (``.xlsx`` or ``.xlsm``).
    include_metadata
        Whether to include a summary worksheet.

    Raises
    ------
    ValueError
        If no design matrix exists or the extension is invalid.
    OSError
        If the workbook cannot be written to disk.

    Examples
    --------
    >>> from industrialstats.designs.factorial import FactorialDesign
    >>> design = FactorialDesign([Factor("A", ["Low", "High"])], randomize=False)
    >>> design.generate_design()
    >>> design.to_excel("factorial_design.xlsx")
    """

    # Disallow export when the design has not yet been generated.
    if self.design_matrix is None:
        raise ValueError("No design matrix available for export.")

    # Resolve and validate the requested output path and extension.
    path = Path(filename)
    if path.suffix.lower() not in {".xlsx", ".xlsm"}:
        raise ValueError("Excel export only supports .xlsx or .xlsm files.")

    # Make a defensive copy so that formatting mutations do not affect the source.
    design_df = self.design_matrix.copy()

    try:
        # Open an Excel writer context for structured output.
        with pd.ExcelWriter(path) as writer:
            # Persist the design matrix on the primary sheet.
            design_df.to_excel(writer, index=False, sheet_name="Design")
            self._format_excel_sheet(writer, "Design", design_df)

            # Append supplementary metadata when requested.
            if include_metadata:
                metadata = self._build_metadata_frame()
                metadata.to_excel(writer, sheet_name="Summary")
                self._format_excel_sheet(writer, "Summary", metadata, autofit=False)
    except ModuleNotFoundError:  # pragma: no cover
        # Fallback gracefully by emitting CSV if Excel engines are missing.
        design_df.to_csv(path.with_suffix(".csv"), index=False)
    except OSError as exc:  # pragma: no cover
        raise OSError(f"Failed to write Excel file '{path}': {exc}") from exc

to_json

to_json(filename: str) -> None

Serialise the design to a JSON document suitable for APIs.

Parameters:

Name	Type	Description	Default
filename	str	Destination path ending with .json.	required

Raises:

Type	Description
ValueError	If the design matrix is missing or the extension is invalid.
OSError	If writing the JSON file fails.

Examples:

>>> from industrialstats.designs.factorial import FactorialDesign
>>> design = FactorialDesign([Factor("A", [1, -1])], randomize=False)
>>> design.generate_design()
>>> design.to_json("design.json")

Source code in src/industrialstats/designs/base.py

def to_json(self, filename: str) -> None:
    """Serialise the design to a JSON document suitable for APIs.

    Parameters
    ----------
    filename
        Destination path ending with ``.json``.

    Raises
    ------
    ValueError
        If the design matrix is missing or the extension is invalid.
    OSError
        If writing the JSON file fails.

    Examples
    --------
    >>> from industrialstats.designs.factorial import FactorialDesign
    >>> design = FactorialDesign([Factor("A", [1, -1])], randomize=False)
    >>> design.generate_design()
    >>> design.to_json("design.json")
    """

    # Confirm that a design matrix exists before constructing the payload.
    if self.design_matrix is None:
        raise ValueError("No design matrix available for export.")

    # Ensure the destination path targets a JSON extension for interoperability.
    path = Path(filename)
    if path.suffix.lower() != ".json":
        raise ValueError("JSON export requires a '.json' file extension.")

    # Build a serialisable payload containing metadata and the design matrix.
    payload = {
        "name": self.name,
        "exported_at": datetime.now(UTC).isoformat().replace("+00:00", "Z"),
        "metadata": self._metadata_payload(),
        "design_matrix": self.design_matrix.to_dict(orient="records"),
    }

    try:
        # Persist the JSON payload with indentation for readability.
        with path.open("w", encoding="utf-8") as handle:
            json.dump(payload, handle, indent=2)
    except OSError as exc:  # pragma: no cover
        raise OSError(f"Failed to write JSON file '{path}': {exc}") from exc

clone

clone() -> ExperimentalDesign

Create and return a deep copy of the design instance.

Returns:

Type	Description
ExperimentalDesign	Deep copy of self with independent factor and matrix structures.

Examples:

>>> from industrialstats.designs.factorial import FactorialDesign
>>> design = FactorialDesign([Factor("A", [0, 1])], randomize=False)
>>> design.generate_design()
>>> clone = design.clone()
>>> clone is design
False
>>> clone.design_matrix.equals(design.design_matrix)
True

Source code in src/industrialstats/designs/base.py

def clone(self) -> "ExperimentalDesign":
    """Create and return a deep copy of the design instance.

    Returns
    -------
    ExperimentalDesign
        Deep copy of ``self`` with independent factor
        and matrix structures.

    Examples
    --------
    >>> from industrialstats.designs.factorial import FactorialDesign
    >>> design = FactorialDesign([Factor("A", [0, 1])], randomize=False)
    >>> design.generate_design()
    >>> clone = design.clone()
    >>> clone is design
    False
    >>> clone.design_matrix.equals(design.design_matrix)
    True
    """

    # Delegate to :func:`copy.deepcopy` to replicate nested containers safely.
    return deepcopy(self)

merge_with

merge_with(other_design: ExperimentalDesign) -> ExperimentalDesign

Merge the design with another compatible design.

Parameters:

Name	Type	Description	Default
other_design	ExperimentalDesign	Design whose runs will be appended to this design.	required

Returns:

Type	Description
ExperimentalDesign	Deep copy of the current design containing runs from both designs.

Raises:

Type	Description
TypeError	If other_design is not an :class:ExperimentalDesign.
ValueError	If either design lacks a generated design matrix or the factor metadata is incompatible.

Examples:

>>> from industrialstats.designs.factorial import FactorialDesign
>>> factors = [Factor("A", [0, 1])]
>>> d1 = FactorialDesign(factors, randomize=False)
>>> _ = d1.generate_design()
>>> d2 = FactorialDesign(factors, randomize=False)
>>> _ = d2.generate_design()
>>> merged = d1.merge_with(d2)
>>> merged.run_count
4

Source code in src/industrialstats/designs/base.py

def merge_with(self, other_design: "ExperimentalDesign") -> "ExperimentalDesign":
    """Merge the design with another compatible design.

    Parameters
    ----------
    other_design
        Design whose runs will be appended to this design.

    Returns
    -------
    ExperimentalDesign
        Deep copy of the current design containing runs
        from both designs.

    Raises
    ------
    TypeError
        If ``other_design`` is not an :class:`ExperimentalDesign`.
    ValueError
        If either design lacks a generated design matrix or the
        factor metadata is incompatible.

    Examples
    --------
    >>> from industrialstats.designs.factorial import FactorialDesign
    >>> factors = [Factor("A", [0, 1])]
    >>> d1 = FactorialDesign(factors, randomize=False)
    >>> _ = d1.generate_design()
    >>> d2 = FactorialDesign(factors, randomize=False)
    >>> _ = d2.generate_design()
    >>> merged = d1.merge_with(d2)
    >>> merged.run_count
    4
    """

    # Validate type compatibility to avoid merging unrelated implementations.
    if not isinstance(other_design, ExperimentalDesign):
        raise TypeError("other_design must be an ExperimentalDesign instance")

    # Ensure both designs are generated prior to concatenation.
    if self.design_matrix is None or other_design.design_matrix is None:
        raise ValueError("Both designs must have generated matrices to merge.")

    # Clone so the current design remains immutable to callers.
    merged_design = self.clone()

    # Merge factor metadata to reconcile level information.
    merged_design.factors = self._merge_factor_metadata(other_design)

    # Concatenate design matrices while respecting differing column orders.
    merged_design.design_matrix = self._aligned_concat(other_design)
    merged_design.design_matrix = merged_design.design_matrix.reset_index(drop=True)

    # Reset randomisation metadata after the merge.
    merged_design.randomized = False
    merged_design.seed = None
    return merged_design

compare_to

compare_to(other_design: ExperimentalDesign) -> dict[str, Any]

Compare this design with another design.

Parameters:

Name	Type	Description	Default
other_design	ExperimentalDesign	Design to compare against.	required

Returns:

Type	Description
Dict[str, Any]	Summary of differences such as run count and factor sets.

Raises:

Type	Description
ValueError	If either design lacks a generated matrix.

Source code in src/industrialstats/designs/base.py

def compare_to(self, other_design: "ExperimentalDesign") -> dict[str, Any]:
    """Compare this design with another design.

    Parameters
    ----------
    other_design
        Design to compare against.

    Returns
    -------
    Dict[str, Any]
        Summary of differences such as run count and factor
        sets.

    Raises
    ------
    ValueError
        If either design lacks a generated matrix.
    """

    # Ensure both designs have been generated before comparing characteristics.
    if self.design_matrix is None or other_design.design_matrix is None:
        raise ValueError("Both designs must have generated matrices to compare")

    # Identify differing factor names between both designs.
    factor_diff = list(
        set(self.factor_names).symmetric_difference(other_design.factor_names)
    )
    # Compute the run-count difference for quick diagnostics.
    run_diff = other_design.run_count - self.run_count
    return {"factor_diff": factor_diff, "run_diff": run_diff}

summary

summary() -> dict[str, Any]

Return summary information about the design.

Returns:

Type	Description
dict[str, Any]	Key characteristics of the design.

Source code in src/industrialstats/designs/base.py

def summary(self) -> dict[str, Any]:
    """Return summary information about the design.

    Returns
    -------
    dict[str, Any]
        Key characteristics of the design.
    """

    # Provide a status indicator when the design has not been generated yet.
    if self.design_matrix is None:
        return {"status": "Design not generated"}

    # Report headline metadata for quick inspection.
    return {
        "design_name": self.name,
        "n_factors": len(self.factors),
        "n_runs": len(self.design_matrix),
        "randomized": self.randomized,
        "factors": [f.name for f in self.factors],
        "factor_levels": {f.name: f.levels for f in self.factors},
        "design_matrix_shape": self.design_matrix.shape,
    }

Full factorial

industrialstats.designs.factorial

Full factorial experimental designs with statistically defined blocking.

FactorialDesign

FactorialDesign(factors: list[Factor], replicates: int = 1, center_points: int = 0, randomize: bool = True, blocks: int | None = None, seed: int | None = None, block_generators: list[str] | None = None, allow_main_effect_confounding: bool = False)

Bases: FactorialDesign

Full factorial design with regular blocking for two-level experiments.

Create a full factorial design.

Source code in src/industrialstats/designs/factorial.py

def __init__(
    self,
    factors: list[Factor],
    replicates: int = 1,
    center_points: int = 0,
    randomize: bool = True,
    blocks: int | None = None,
    seed: int | None = None,
    block_generators: list[str] | None = None,
    allow_main_effect_confounding: bool = False,
) -> None:
    """Create a full factorial design."""
    super().__init__(
        factors=factors,
        replicates=replicates,
        center_points=center_points,
        randomize=randomize,
        blocks=None,
        seed=seed,
    )
    self.blocks = blocks
    self.block_generators = (
        list(block_generators) if block_generators is not None else None
    )
    self.allow_main_effect_confounding = allow_main_effect_confounding

    if not isinstance(allow_main_effect_confounding, bool):
        raise ValueError("allow_main_effect_confounding must be a boolean")
    if self.block_generators is not None and self.blocks is None:
        self.blocks = 2 ** len(self.block_generators)
    self._validate_blocking_configuration()

generate_design

generate_design() -> DataFrame

Generate the full factorial design and apply regular blocking.

Source code in src/industrialstats/designs/factorial.py

def generate_design(self) -> pd.DataFrame:
    """Generate the full factorial design and apply regular blocking."""
    if not self.validate_design():
        raise ValueError("Invalid design configuration")

    combinations_list = list(product(*(factor.levels for factor in self.factors)))
    design_data: list[dict[str, Any]] = []
    run_id = 1

    for replicate in range(1, self.replicates + 1):
        for combination in combinations_list:
            row: dict[str, Any] = {
                "RunID": run_id,
                "Replicate": replicate,
                "DesignPoint": "Factorial",
                "StdOrder": run_id,
            }
            for index, factor in enumerate(self.factors):
                row[factor.name] = combination[index]
            design_data.append(row)
            run_id += 1

    if self.center_points > 0:
        centers = self._calculate_center_points()
        for _ in range(self.center_points):
            row = {
                "RunID": run_id,
                "Replicate": 1,
                "DesignPoint": "Center",
                "StdOrder": run_id,
            }
            for index, factor in enumerate(self.factors):
                row[factor.name] = centers[index]
            design_data.append(row)
            run_id += 1

    self.design_matrix = pd.DataFrame(design_data)
    if self._blocking_requested():
        self._apply_regular_blocking()

    if self.randomize_flag:
        if self._blocking_requested():
            self._randomize_within_blocks()
        else:
            self.randomize(self.seed)

    self._store_factor_level_orders()
    return self.design_matrix

validate_design

validate_design() -> bool

Validate the factorial design and its optional block structure.

Source code in src/industrialstats/designs/factorial.py

def validate_design(self) -> bool:
    """Validate the factorial design and its optional block structure."""
    if not self.factors or self.replicates < 1:
        return False
    if any(len(factor.levels) < 2 for factor in self.factors):
        return False
    try:
        self._validate_blocking_configuration()
    except ValueError:
        return False
    return True

model_terms

model_terms(max_order: int | None = None) -> list[str]

Return hierarchical factorial terms through max_order.

Parameters:

Name	Type	Description	Default
max_order	int | None	Highest interaction order to include. None returns the saturated factorial hierarchy through the interaction involving every factor.	None

Returns:

Type	Description
list[str]	Ordered effect names such as A, A*B and A*B*C.

Source code in src/industrialstats/designs/factorial.py

def model_terms(self, max_order: int | None = None) -> list[str]:
    """Return hierarchical factorial terms through ``max_order``.

    Parameters
    ----------
    max_order
        Highest interaction order to include. ``None`` returns the saturated
        factorial hierarchy through the interaction involving every factor.

    Returns
    -------
    list[str]
        Ordered effect names such as ``A``, ``A*B`` and ``A*B*C``.
    """
    order = self._resolve_model_order(max_order)
    names = [factor.name for factor in self.factors]
    return [
        "*".join(term)
        for interaction_order in range(1, order + 1)
        for term in combinations(names, interaction_order)
    ]

degrees_of_freedom

degrees_of_freedom(max_order: int | None = None) -> dict[str, int]

Calculate factorial degrees of freedom for a hierarchical model.

For a term involving factors in a set S, the term degrees of freedom are the product prod(len(levels_j) - 1 for j in S). With max_order=None the model is saturated over the factorial treatment combinations. A finite max_order gives a truncated hierarchical model; omitted higher-order treatment variation remains in Error together with replication and center-point residual degrees of freedom.

Source code in src/industrialstats/designs/factorial.py

def degrees_of_freedom(self, max_order: int | None = None) -> dict[str, int]:
    """Calculate factorial degrees of freedom for a hierarchical model.

    For a term involving factors in a set ``S``, the term degrees of freedom
    are the product ``prod(len(levels_j) - 1 for j in S)``. With
    ``max_order=None`` the model is saturated over the factorial treatment
    combinations. A finite ``max_order`` gives a truncated hierarchical model;
    omitted higher-order treatment variation remains in ``Error`` together
    with replication and center-point residual degrees of freedom.
    """
    order = self._resolve_model_order(max_order)
    factor_by_name = {factor.name: factor for factor in self.factors}
    dof: dict[str, int] = {}

    for term_name in self.model_terms(order):
        term_factors = term_name.split("*")
        dof[term_name] = prod(
            len(factor_by_name[name].levels) - 1 for name in term_factors
        )

    total_runs = int(self.n_runs())
    model_dof = sum(dof.values())
    dof["Error"] = total_runs - model_dof - 1
    dof["Total"] = total_runs - 1
    return dof

model_structure

model_structure(max_order: int | None = None) -> dict[str, Any]

Describe a saturated or truncated hierarchical factorial model.

Source code in src/industrialstats/designs/factorial.py

def model_structure(self, max_order: int | None = None) -> dict[str, Any]:
    """Describe a saturated or truncated hierarchical factorial model."""
    order = self._resolve_model_order(max_order)
    dof = self.degrees_of_freedom(order)
    terms = self.model_terms(order)
    return {
        "max_order": order,
        "saturated": order == len(self.factors),
        "terms": terms,
        "model_degrees_of_freedom": sum(dof[term] for term in terms),
        "error_degrees_of_freedom": dof["Error"],
        "total_degrees_of_freedom": dof["Total"],
    }

calculate_effects

calculate_effects(response_data: list[float], max_order: int = 2) -> dict[str, float]

Calculate canonical two-level factorial effects.

Factor levels are coded according to their declared order: the first level is -1 and the second is +1. The returned factorial effect is twice the corresponding coefficient in a regression using these coded columns and their products.

Source code in src/industrialstats/designs/factorial.py

def calculate_effects(
    self,
    response_data: list[float],
    max_order: int = 2,
) -> dict[str, float]:
    """Calculate canonical two-level factorial effects.

    Factor levels are coded according to their declared order: the first
    level is ``-1`` and the second is ``+1``. The returned factorial effect
    is twice the corresponding coefficient in a regression using these
    coded columns and their products.
    """
    if self.design_matrix is None:
        raise ValueError("Design matrix not generated")
    if not self._is_two_level_design():
        raise ValueError("Effect calculation only supported for 2-level designs")
    if not (self.design_matrix["DesignPoint"] == "Factorial").all():
        raise ValueError(
            "Canonical factorial effects require factorial treatment points only"
        )

    return calculate_two_level_factorial_effects(
        self.design_matrix,
        response_data,
        [factor.name for factor in self.factors],
        max_order=max_order,
        level_orders={factor.name: factor.levels for factor in self.factors},
    )

block_structure

block_structure() -> dict[str, Any]

Return generators, defining contrasts, and block diagnostics.

Source code in src/industrialstats/designs/factorial.py

def block_structure(self) -> dict[str, Any]:
    """Return generators, defining contrasts, and block diagnostics."""
    if not self._blocking_requested():
        return {
            "n_blocks": 1,
            "generators": [],
            "defining_contrasts": [],
            "confounded_main_effects": [],
            "runs_per_block": int(self.n_factorial_runs()),
        }

    assert self.blocks is not None
    masks = self._resolve_block_masks()
    defining = self._defining_masks(masks)
    return {
        "n_blocks": self.blocks,
        "generators": [self._mask_to_word(mask) for mask in masks],
        "defining_contrasts": [self._mask_to_word(mask) for mask in defining],
        "confounded_main_effects": [
            self._mask_to_word(mask) for mask in defining if mask.bit_count() == 1
        ],
        "runs_per_block": int(self.n_factorial_runs() // self.blocks),
    }

blocking_scheme

blocking_scheme(block_size: int) -> DataFrame

Reassign an existing design using regular treatment contrasts.

Source code in src/industrialstats/designs/factorial.py

def blocking_scheme(self, block_size: int) -> pd.DataFrame:
    """Reassign an existing design using regular treatment contrasts."""
    if self.design_matrix is None:
        raise ValueError("Design matrix not generated")
    if isinstance(block_size, bool) or not isinstance(block_size, int):
        raise ValueError("block_size must be a positive integer")
    if block_size <= 0:
        raise ValueError("block_size must be a positive integer")
    if not (self.design_matrix["DesignPoint"] == "Factorial").all():
        raise ValueError("blocking_scheme requires factorial treatment points only")

    n_runs = len(self.design_matrix)
    if n_runs % block_size != 0:
        raise ValueError("block_size must divide the number of factorial runs")
    n_blocks = n_runs // block_size
    if n_blocks < 2:
        raise ValueError("block_size must define at least two blocks")

    self.blocks = n_blocks
    self.block_generators = None
    self._validate_blocking_configuration()

    if "RunOrder" in self.design_matrix.columns:
        self.design_matrix = self.design_matrix.drop(columns="RunOrder")
    self.design_matrix = self.design_matrix.sort_values(
        "StdOrder", kind="stable"
    ).reset_index(drop=True)
    self._apply_regular_blocking()
    if self.randomize_flag:
        self._randomize_within_blocks()
    self._store_factor_level_orders()
    return self.design_matrix

Fractional factorial

industrialstats.designs.fractional_factorial

Fractional factorial design implementation.

FractionalFactorialDesign

FractionalFactorialDesign(factors: list[Factor], fraction: str = '1/2', generators: list[str] | None = None, resolution: int | None = None, replicates: int = 1, randomize: bool = True)

Bases: ExperimentalDesign

Two-level fractional factorial design.

This class generates regular two-level fractional factorial designs for three to fifteen factors and automatically selects minimum aberration generator strings for the most common fractions (:math:1/2, :math:1/4, :math:1/8, and :math:1/16). The implementation mirrors the confounding analysis strategy described in :mod:industrialstats.designs.README by constructing the defining relation of the design and deriving alias chains through linear algebra over :math:\mathrm{GF}(2).

Notes

• All factors must be two-level factors with exactly two entries in Factor.levels.
• Fraction denominators must be powers of two; the numerator is assumed to be one.
• When generators are not supplied, the class performs a constrained search for generator sets that maximise the design resolution and, as a tie-breaker, minimise the word-length pattern (the minimum aberration criterion of Montgomery and Wu & Hamada).
• Alias structures reported by :meth:alias_structure include all factorial effects, not only main effects and two-factor interactions.

References

Montgomery, D. C. (2017). Design and Analysis of Experiments (9th ed.). Wiley. Wu, C. F. J., & Hamada, M. S. (2009). Experiments: Planning, Analysis, and Optimization (2nd ed.). Wiley. Groemping, U. (2014). "R package FrF2 for creating and analysing 2-level factorial designs". Journal of Statistical Software, 56(1). Xu, H. (2005). "A catalogue of three-level and four-level fractional factorial designs with minimum aberration". Technometrics, 47(3). NIST/SEMATECH (2012). e-Handbook of Statistical Methods, Chapter 5.3. "FrF2" R package documentation for catlg52 (minimum aberration tables). "Montgomery" Chapter 8 examples on foldover strategies. "NIST foldover guidance" on mitigating aliasing. "Xu & Wu (2001)" alias chain construction using defining relations.

Examples:

>>> from industrialstats.designs.base import Factor
>>> factors = [Factor(name, [-1, 1]) for name in "ABCDEFG"]
>>> design = FractionalFactorialDesign(factors, fraction="1/8")
>>> design.generators
['A*B*C', 'A*B*D', 'A*C*D']
>>> metrics = design.resolution_analysis()
>>> metrics["resolution"]
4
>>> alias = design.alias_structure()
>>> alias["A"][:3]
['A', 'A:B:C:F:G', 'A:B:D:E:G']
>>> design.foldover_options()[0]["type"]
'full'

Initialize fractional factorial design.

Parameters:

Name	Type	Description	Default
factors	list[Factor]	Factors in the experiment. Must all have two levels.	required
fraction	str	Fraction of the full design, e.g. "1/2" or "1/4".	'1/2'
generators	list[str]	Generator strings using factor names.	None
resolution	int	Desired design resolution (for reference only).	None
replicates	int	Number of replicates.	1
randomize	bool	Whether to randomize run order.	True

Source code in src/industrialstats/designs/fractional_factorial.py

def __init__(
    self,
    factors: list[Factor],
    fraction: str = "1/2",
    generators: list[str] | None = None,
    resolution: int | None = None,
    replicates: int = 1,
    randomize: bool = True,
) -> None:
    """Initialize fractional factorial design.

    Parameters
    ----------
    factors : list[Factor]
        Factors in the experiment. Must all have two levels.
    fraction : str, optional
        Fraction of the full design, e.g. ``"1/2"`` or ``"1/4"``.
    generators : list[str], optional
        Generator strings using factor names.
    resolution : int, optional
        Desired design resolution (for reference only).
    replicates : int, optional
        Number of replicates.
    randomize : bool, optional
        Whether to randomize run order.
    """
    super().__init__("Fractional Factorial Design")
    self.factors = factors
    self.fraction = fraction
    self.generators = generators or []
    self.requested_resolution = resolution
    self.replicates = replicates
    self.randomize_flag = randomize
    self._coded_matrix: pd.DataFrame | None = None

    factor_count = len(factors)
    if factor_count < 3 or factor_count > 15:
        raise ValueError(
            "Fractional factorial designs support between 3 and 15 factors"
        )
    if not all(len(f.levels) == 2 for f in factors):
        raise ValueError("All factors must have exactly two levels")

    try:
        numerator, denominator = fraction.split("/")
    except ValueError as exc:  # pragma: no cover - defensive
        raise ValueError(
            "Fraction must be provided as 'numerator/denominator'"
        ) from exc
    if numerator.strip() != "1":
        raise ValueError("Only regular fractions with numerator 1 are supported")
    denom = int(denominator)
    p = int(np.log2(denom))
    if 2**p != denom:
        raise ValueError("Fraction denominator must be a power of 2")
    self.p = p

    self._base_count = factor_count - p
    if self._base_count <= 0:
        raise ValueError(
            "Number of generators exceeds number of available base factors"
        )

    if not self.generators:
        self.generators = self._auto_generators()

calculate_resolution

calculate_resolution() -> tuple[int | None, dict[int, int]]

Return the design resolution and its word-length pattern.

Source code in src/industrialstats/designs/fractional_factorial.py

def calculate_resolution(self) -> tuple[int | None, dict[int, int]]:
    """Return the design resolution and its word-length pattern."""

    pattern = _word_length_pattern(self._defining_words())
    if not pattern:
        return None, {}
    return min(pattern), pattern

verify_resolution

verify_resolution(minimum: int) -> bool

Check that the design resolution meets minimum.

Raises:

Type	Description
ValueError	If the design resolution is undefined or below minimum.

Source code in src/industrialstats/designs/fractional_factorial.py

def verify_resolution(self, minimum: int) -> bool:
    """Check that the design resolution meets ``minimum``.

    Raises
    ------
    ValueError
        If the design resolution is undefined or below ``minimum``.
    """

    resolution, _ = self.calculate_resolution()
    if resolution is None:
        raise ValueError("Resolution is undefined for a design without generators")
    if resolution < minimum:
        raise ValueError(
            f"Design resolution {resolution} is below the required minimum of {minimum}."
        )
    return True

generate_design

generate_design() -> DataFrame

Generate fractional factorial design matrix.

Source code in src/industrialstats/designs/fractional_factorial.py

def generate_design(self) -> pd.DataFrame:
    """Generate fractional factorial design matrix."""
    base_factors = self.factors[: self._base_count]
    alias_factors = self.factors[self._base_count :]

    base_names = [factor.name for factor in base_factors]
    parsed_generators = self._parsed_generators()
    runs = list(product([-1, 1], repeat=len(base_names)))
    data = []
    run_id = 1
    for rep in range(1, self.replicates + 1):
        for run in runs:
            row = dict(zip(base_names, run, strict=True))
            for factor, terms in zip(alias_factors, parsed_generators, strict=True):
                row[factor.name] = self._evaluate_generator(terms, row)
            data.append({"RunID": run_id, "Replicate": rep, **row})
            run_id += 1

    coded_df = pd.DataFrame(data)
    df = coded_df.copy()
    for factor in self.factors:
        mapping = {-1: factor.levels[0], 1: factor.levels[1]}
        df[factor.name] = df[factor.name].map(mapping)

    if self.randomize_flag:
        df = df.sample(frac=1, random_state=None).reset_index(drop=True)
        df.insert(0, "RunOrder", range(1, len(df) + 1))
        self.randomized = True
    self.design_matrix = df
    self._coded_matrix = coded_df
    return df

validate_design

validate_design() -> bool

Validate fractional factorial parameters.

Source code in src/industrialstats/designs/fractional_factorial.py

def validate_design(self) -> bool:
    """Validate fractional factorial parameters."""
    return (
        3 <= len(self.factors) <= 15
        and all(len(f.levels) == 2 for f in self.factors)
        and self.replicates > 0
    )

alias_structure

alias_structure() -> dict[str, list[str]]

Return alias chains for all factorial effects.

The alias chains are computed by applying the defining relation of the design to every factorial effect. Each entry contains the canonical effect (lexicographically smallest mask) followed by the remaining members of its alias class.

Source code in src/industrialstats/designs/fractional_factorial.py

def alias_structure(self) -> dict[str, list[str]]:
    """Return alias chains for all factorial effects.

    The alias chains are computed by applying the defining relation of the
    design to every factorial effect.  Each entry contains the canonical
    effect (lexicographically smallest mask) followed by the remaining
    members of its alias class.
    """

    words = self._defining_words()
    total_factors = len(self.factors)
    alias_map: dict[str, list[str]] = {}
    processed: set[int] = set()
    for mask in range(1, 1 << total_factors):
        if mask in processed:
            continue
        alias_class = {mask}
        for word in words:
            alias_class.add(mask ^ word)
        canonical = min(alias_class)
        if canonical != mask:
            continue
        canonical_name = self._mask_to_effect(canonical)
        effect_names = sorted(
            self._mask_to_effect(alias_mask) for alias_mask in alias_class
        )
        ordered = [canonical_name] + [
            name for name in effect_names if name != canonical_name
        ]
        alias_map[canonical_name] = ordered
        processed.update(alias_class)
    return alias_map

resolution_analysis

resolution_analysis() -> dict[str, Any]

Analyze design resolution and clarity.

Source code in src/industrialstats/designs/fractional_factorial.py

def resolution_analysis(self) -> dict[str, Any]:
    """Analyze design resolution and clarity."""
    resolution, pattern = self.calculate_resolution()
    if resolution is None:
        return {
            "resolution": None,
            "word_length_pattern": {},
            "minimum_aberration": [],
            "meets_requested_resolution": None,
        }
    max_length = max(pattern)
    aberration = [
        (length, pattern.get(length, 0))
        for length in range(resolution, max_length + 1)
    ]
    meets = None
    if self.requested_resolution is not None:
        meets = resolution >= self.requested_resolution
    return {
        "resolution": resolution,
        "word_length_pattern": pattern,
        "minimum_aberration": aberration,
        "meets_requested_resolution": meets,
    }

foldover_options

foldover_options() -> list[dict[str, Any]]

Suggest foldover strategies with supporting metadata.

Foldover proposals follow Montgomery's guidance: a full foldover adds a replicate with all signs reversed, while partial foldovers target problematic columns. The suggestions focus on breaking aliases for main effects with the longest alias chains.

Source code in src/industrialstats/designs/fractional_factorial.py

def foldover_options(self) -> list[dict[str, Any]]:
    """Suggest foldover strategies with supporting metadata.

    Foldover proposals follow Montgomery's guidance: a *full* foldover adds
    a replicate with all signs reversed, while *partial* foldovers target
    problematic columns.  The suggestions focus on breaking aliases for main
    effects with the longest alias chains.
    """

    resolution_info = self.resolution_analysis()
    resolution = resolution_info["resolution"]
    alias_map = self.alias_structure()
    alias_names = [factor.name for factor in self.factors[self._base_count :]]
    parsed = self._parsed_generators()

    severity: list[tuple[str, int, list[str]]] = []
    for factor in self.factors:
        chain = alias_map.get(factor.name, [factor.name])
        severity.append((factor.name, len(chain) - 1, chain))
    severity.sort(key=lambda item: (-item[1], item[0]))

    # Full foldover suggestion.
    expected_resolution = (
        None if resolution is None else min(resolution + 1, len(self.factors))
    )
    options: list[dict[str, Any]] = [
        {
            "type": "full",
            "description": (
                "Add a full foldover by reversing the signs of all generators. "
                "This mitigates main-effect aliasing and typically improves the resolution by one."
            ),
            "generators_to_reverse": self.generators,
            "expected_resolution": expected_resolution,
        }
    ]

    # Partial foldovers for the most aliased main effects.
    for factor_name, _, chain in severity[: min(3, len(severity))]:
        impacted_generators = []
        for gen, terms, alias_name in zip(
            self.generators, parsed, alias_names, strict=True
        ):
            if factor_name == alias_name or factor_name in terms:
                impacted_generators.append(gen)
        options.append(
            {
                "type": "partial",
                "factor": factor_name,
                "description": (
                    f"Fold over factor {factor_name} (reverse signs of runs where it is high) "
                    "to separate it from its aliases."
                ),
                "generators_to_reverse": impacted_generators,
                "confounded_with": [
                    effect for effect in chain if effect != factor_name
                ],
            }
        )

    return options

Completely randomized design

industrialstats.designs.crd

Completely Randomized Design (CRD) implementation.

CompletelyRandomizedDesign

CompletelyRandomizedDesign(treatments: list[str], replicates: int, seed: int | None = None, response_variables: list[str] | None = None)

Bases: ExperimentalDesign

Completely Randomized Design (CRD).

This is the simplest experimental design where treatments are randomly assigned to experimental units without any restrictions or blocking.

Initialize CRD.

Parameters:

Name	Type	Description	Default
treatments	list of str	Names of treatment levels.	required
replicates	int	Number of replicates per treatment.	required
seed	int	Random seed for reproducible run ordering.	None
response_variables	list of str	Names of response variables measured in the experiment.	None

Source code in src/industrialstats/designs/crd.py

def __init__(
    self,
    treatments: list[str],
    replicates: int,
    seed: int | None = None,
    response_variables: list[str] | None = None,
) -> None:
    """Initialize CRD.

    Parameters
    ----------
    treatments : list of str
        Names of treatment levels.
    replicates : int
        Number of replicates per treatment.
    seed : int, optional
        Random seed for reproducible run ordering.
    response_variables : list of str, optional
        Names of response variables measured in the experiment.
    """
    super().__init__("Completely Randomized Design")

    if len(treatments) < 2:
        raise ValueError("Must have at least 2 treatments")
    if replicates < 1:
        raise ValueError("Must have at least 1 replicate")

    # Copy the caller's list so a later mutation cannot desynchronise
    # the design from the factor levels derived from it.
    self.treatments = list(treatments)
    self.replicates = replicates
    self.seed = seed
    self.response_variables = response_variables or []

    # Create a single factor with treatment levels
    treatment_levels: list[str | float | int] = list(self.treatments)
    treatment_factor = Factor("Treatment", treatment_levels, "categorical")
    self.factors = [treatment_factor]

generate_design

generate_design() -> DataFrame

Generate the CRD design matrix.

Returns:

Type	Description
DataFrame	Design matrix with randomized run order.

Source code in src/industrialstats/designs/crd.py

def generate_design(self) -> pd.DataFrame:
    """Generate the CRD design matrix.

    Returns
    -------
    pandas.DataFrame
        Design matrix with randomized run order.
    """
    if not self.validate_design():
        raise ValueError("Invalid design configuration")

    # Create all treatment-replicate combinations
    design_data = []
    run_id = 1

    for rep in range(1, self.replicates + 1):
        for treatment in self.treatments:
            design_data.append(
                {"RunID": run_id, "Treatment": treatment, "Replicate": rep}
            )
            run_id += 1

    self.design_matrix = pd.DataFrame(design_data)

    # Always randomize CRD (that's the point!)
    self.randomize(seed=self.seed)

    return self.design_matrix

validate_design

validate_design() -> bool

Validate CRD parameters.

Source code in src/industrialstats/designs/crd.py

def validate_design(self) -> bool:
    """Validate CRD parameters."""
    if len(self.treatments) < 2:
        return False
    if self.replicates < 1:
        return False
    if not all(isinstance(r, str) for r in self.response_variables):
        return False
    return len(self.response_variables) == len(set(self.response_variables))

n_runs

n_runs() -> int

Calculate total number of runs.

Source code in src/industrialstats/designs/crd.py

def n_runs(self) -> int:
    """Calculate total number of runs."""
    return len(self.treatments) * self.replicates

degrees_of_freedom

degrees_of_freedom() -> dict[str, int]

Calculate degrees of freedom for CRD analysis.

Source code in src/industrialstats/designs/crd.py

def degrees_of_freedom(self) -> dict[str, int]:
    """Calculate degrees of freedom for CRD analysis."""
    n_treatments = len(self.treatments)
    total_runs = self.n_runs()

    return {
        "Treatment": n_treatments - 1,
        "Error": total_runs - n_treatments,
        "Total": total_runs - 1,
    }

expected_mean_squares

expected_mean_squares() -> dict[str, str]

Return expected mean squares for CRD.

Source code in src/industrialstats/designs/crd.py

def expected_mean_squares(self) -> dict[str, str]:
    """Return expected mean squares for CRD."""
    return {"Treatment": "σ² + r·σ²ₜ", "Error": "σ²"}

efficiency_vs_rcbd

efficiency_vs_rcbd(block_variance: float, error_variance: float) -> float

Calculate relative efficiency compared to RCBD.

Parameters:

Name	Type	Description	Default
block_variance	float	Estimated variance between blocks.	required
error_variance	float	Estimated experimental error variance.	required

Returns:

Type	Description
float	Relative efficiency (> 1 means CRD is more efficient).

Source code in src/industrialstats/designs/crd.py

def efficiency_vs_rcbd(self, block_variance: float, error_variance: float) -> float:
    """Calculate relative efficiency compared to RCBD.

    Parameters
    ----------
    block_variance : float
        Estimated variance between blocks.
    error_variance : float
        Estimated experimental error variance.

    Returns
    -------
    float
        Relative efficiency (> 1 means CRD is more efficient).
    """
    # Relative efficiency = (RCBD error MS) / (CRD error MS)
    rcbd_error_ms = error_variance
    crd_error_ms = error_variance + block_variance

    return crd_error_ms / rcbd_error_ms

sample_size_calculation

sample_size_calculation(effect_size: float, alpha: float = 0.05, power: float = 0.8) -> int

Calculate required sample size per treatment.

Parameters:

Name	Type	Description	Default
effect_size	float	Expected effect size (Cohen's f).	required
alpha	float	Type I error rate. Defaults to 0.05.	0.05
power	float	Desired statistical power. Defaults to 0.8.	0.8

Returns:

Type	Description
int	Required number of replicates per treatment.

Source code in src/industrialstats/designs/crd.py

def sample_size_calculation(
    self, effect_size: float, alpha: float = 0.05, power: float = 0.8
) -> int:
    """Calculate required sample size per treatment.

    Parameters
    ----------
    effect_size : float
        Expected effect size (Cohen's ``f``).
    alpha : float, optional
        Type I error rate. Defaults to 0.05.
    power : float, optional
        Desired statistical power. Defaults to 0.8.

    Returns
    -------
    int
        Required number of replicates per treatment.
    """

    from scipy.stats import f

    k = len(self.treatments)  # Number of treatments

    # Degrees of freedom
    df1 = k - 1

    # Non-centrality parameter for desired power
    from scipy.stats import ncf

    # Iterative search for required sample size
    for n_per_group in range(2, 1000):
        df2 = k * (n_per_group - 1)
        lambda_nc = effect_size**2 * n_per_group * k

        f_crit = f.ppf(1 - alpha, df1, df2)
        calculated_power = 1 - ncf.cdf(f_crit, df1, df2, lambda_nc)

        if calculated_power >= power:
            return n_per_group

    return -1  # Could not find required sample size

create_data_collection_sheet

create_data_collection_sheet(response_variables: list[str] | None = None) -> DataFrame

Create a data collection sheet for the experiment.

Parameters:

Name	Type	Description	Default
response_variables	list of str	Names of response variables to include. Defaults to the design's stored response_variables or ["Response"] if none were specified.	None

Returns:

Type	Description
DataFrame	Data collection sheet with empty response columns.

Source code in src/industrialstats/designs/crd.py

def create_data_collection_sheet(
    self, response_variables: list[str] | None = None
) -> pd.DataFrame:
    """Create a data collection sheet for the experiment.

    Parameters
    ----------
    response_variables : list of str, optional
        Names of response variables to include. Defaults to the design's
        stored ``response_variables`` or ``["Response"]`` if none were
        specified.

    Returns
    -------
    pandas.DataFrame
        Data collection sheet with empty response columns.
    """
    design_matrix = self.design_matrix
    if design_matrix is None:
        design_matrix = self.generate_design()

    data_sheet = design_matrix.copy()

    responses = response_variables or self.response_variables or ["Response"]

    for response in responses:
        data_sheet[response] = np.nan

    # Add columns for data collection
    data_sheet["Date"] = ""
    data_sheet["Time"] = ""
    data_sheet["Observer"] = ""
    data_sheet["Notes"] = ""

    return data_sheet

summary_statistics

summary_statistics(data: DataFrame, response_columns: list[str]) -> dict[str, DataFrame]

Calculate summary statistics for multiple responses.

Parameters:

Name	Type	Description	Default
data	DataFrame	Experimental data with results.	required
response_columns	list of str	Names of response variable columns to analyze.	required

Returns:

Type	Description
dict of pandas.DataFrame	Mapping of response names to summary statistics by treatment.

Source code in src/industrialstats/designs/crd.py

def summary_statistics(
    self, data: pd.DataFrame, response_columns: list[str]
) -> dict[str, pd.DataFrame]:
    """Calculate summary statistics for multiple responses.

    Parameters
    ----------
    data : pandas.DataFrame
        Experimental data with results.
    response_columns : list of str
        Names of response variable columns to analyze.

    Returns
    -------
    dict of pandas.DataFrame
        Mapping of response names to summary statistics by treatment.
    """
    self._validate_response_data(data, response_columns)

    summaries: dict[str, pd.DataFrame] = {}
    from scipy.stats import t

    for column in response_columns:
        summary = (
            data.groupby("Treatment")[column]
            .agg(["count", "mean", "std", "min", "max", "median"])
            .round(3)
        )

        ci_lower: list[float] = []
        ci_upper: list[float] = []

        for treatment in summary.index:
            treatment_data = data[data["Treatment"] == treatment][column]
            n = len(treatment_data)
            mean = treatment_data.mean()
            std = treatment_data.std()

            t_critical = t.ppf(0.975, n - 1)
            margin_error = t_critical * std / np.sqrt(n)

            ci_lower.append(mean - margin_error)
            ci_upper.append(mean + margin_error)

        summary["CI_Lower"] = np.round(ci_lower, 3)
        summary["CI_Upper"] = np.round(ci_upper, 3)

        summaries[column] = summary

    return summaries

Randomized complete block design

industrialstats.designs.rcbd

Randomized Complete Block Design implementation.

RandomizedCompleteBlockDesign

RandomizedCompleteBlockDesign(treatments: list[str], blocks: list[str], blocking_factor: str = 'Block', seed: int | None = None)

Bases: ExperimentalDesign

Randomized Complete Block Design.

Initialize RCBD.

Parameters:

Name	Type	Description	Default
treatments	list of str	List of treatment names or levels.	required
blocks	list of str	Names of blocking levels.	required
blocking_factor	str	Column name for blocks in the design matrix. Defaults to "Block".	'Block'
seed	int	Random seed for reproducibility.	None

Raises:

Type	Description
ValueError	If fewer than two treatments or blocks are provided.

Source code in src/industrialstats/designs/rcbd.py

def __init__(
    self,
    treatments: list[str],
    blocks: list[str],
    blocking_factor: str = "Block",
    seed: int | None = None,
) -> None:
    """Initialize RCBD.

    Parameters
    ----------
    treatments : list of str
        List of treatment names or levels.
    blocks : list of str
        Names of blocking levels.
    blocking_factor : str, optional
        Column name for blocks in the design matrix. Defaults to ``"Block"``.
    seed : int, optional
        Random seed for reproducibility.

    Raises
    ------
    ValueError
        If fewer than two treatments or blocks are provided.
    """
    super().__init__("Randomized Complete Block Design")
    if len(treatments) < 2:
        raise ValueError("Must have at least 2 treatments")
    if len(blocks) < 2:
        raise ValueError("Must have at least 2 blocks")

    # Copy the caller's lists so a later mutation cannot desynchronise the
    # design from the factor levels derived from them.
    self.treatments = list(treatments)
    self.blocks = list(blocks)
    self.blocking_factor = blocking_factor
    self.seed = seed

    treatment_levels: list[str | float | int] = list(self.treatments)
    block_levels: list[str | float | int] = list(self.blocks)
    self.factors = [
        Factor("Treatment", treatment_levels, "categorical"),
        Factor(blocking_factor, block_levels, "categorical"),
    ]

generate_design

generate_design(seed: int | None = None) -> DataFrame

Generate RCBD matrix with proper randomization.

Each block contains all treatments exactly once. Blocks are internally randomized using :func:pandas.DataFrame.sample with an offset seed so that seed + i controls the shuffling of the i-th block. The combined design matrix is returned with a leading RunOrder column.

Parameters:

Name	Type	Description	Default
seed	int	Random seed for reproducible shuffling. If None (default), the method uses the seed provided during initialization.	None

Returns:

Type	Description
DataFrame	Randomized design matrix with columns RunOrder, blocking_factor and Treatment.

See Also

industrialstats.designs.base.ExperimentalDesign.randomize Generic randomization utility for arbitrary designs. efficiency_vs_crd Compute relative efficiency against a completely randomized design.

Examples:

>>> rcbd = RandomizedCompleteBlockDesign(
...     treatments=["A", "B", "C"], blocks=["B1", "B2"], seed=123
... )
>>> dm = rcbd.generate_design()
>>> dm.head()
   RunOrder Block Treatment
0         1    B1         B
1         2    B1         C
2         3    B1         A

References

.. [1] Montgomery, D.C. (2017). Design and Analysis of Experiments. 9th ed. Wiley.

Source code in src/industrialstats/designs/rcbd.py

def generate_design(self, seed: int | None = None) -> pd.DataFrame:
    """Generate RCBD matrix with proper randomization.

    Each block contains all treatments exactly once. Blocks are internally
    randomized using :func:`pandas.DataFrame.sample` with an offset seed so
    that ``seed + i`` controls the shuffling of the ``i``-th block. The
    combined design matrix is returned with a leading ``RunOrder`` column.

    Parameters
    ----------
    seed : int, optional
        Random seed for reproducible shuffling. If ``None`` (default), the
        method uses the seed provided during initialization.

    Returns
    -------
    pandas.DataFrame
        Randomized design matrix with columns ``RunOrder``,
        ``blocking_factor`` and ``Treatment``.

    See Also
    --------
    industrialstats.designs.base.ExperimentalDesign.randomize
        Generic randomization utility for arbitrary designs.
    efficiency_vs_crd
        Compute relative efficiency against a completely randomized design.

    Examples
    --------
    >>> rcbd = RandomizedCompleteBlockDesign(
    ...     treatments=["A", "B", "C"], blocks=["B1", "B2"], seed=123
    ... )
    >>> dm = rcbd.generate_design()
    >>> dm.head()
       RunOrder Block Treatment
    0         1    B1         B
    1         2    B1         C
    2         3    B1         A

    References
    ----------
    .. [1] Montgomery, D.C. (2017). *Design and Analysis of Experiments*.
           9th ed. Wiley.
    """
    if seed is not None:
        self.seed = seed

    design_rows = []
    run_id = 1
    for i, block in enumerate(self.blocks):
        block_rows = []
        for treatment in self.treatments:
            block_rows.append(
                {
                    "RunID": run_id,
                    self.blocking_factor: block,
                    "Treatment": treatment,
                }
            )
            run_id += 1
        # Randomize within block
        block_df = (
            pd.DataFrame(block_rows)
            .sample(
                frac=1,
                random_state=None if self.seed is None else self.seed + i,
            )
            .reset_index(drop=True)
        )
        design_rows.extend(block_df.to_dict(orient="records"))

    self.design_matrix = pd.DataFrame(design_rows)
    self.design_matrix.insert(0, "RunOrder", range(1, len(self.design_matrix) + 1))
    self.randomized = True
    return self.design_matrix

validate_design

validate_design() -> bool

Validate RCBD parameters.

Source code in src/industrialstats/designs/rcbd.py

def validate_design(self) -> bool:
    """Validate RCBD parameters."""
    return len(self.treatments) >= 2 and len(self.blocks) >= 2

efficiency_vs_crd

efficiency_vs_crd(block_variance: float, error_variance: float = 1.0) -> float

Calculate relative efficiency compared to CRD.

Source code in src/industrialstats/designs/rcbd.py

def efficiency_vs_crd(
    self, block_variance: float, error_variance: float = 1.0
) -> float:
    """Calculate relative efficiency compared to CRD."""
    if block_variance < 0 or error_variance <= 0:
        raise ValueError("variances must be positive")
    rcbd_error_ms = error_variance
    crd_error_ms = error_variance + block_variance
    return rcbd_error_ms / crd_error_ms

missing_plot_analysis

missing_plot_analysis(missing_positions: list[tuple[str, str]]) -> dict[str, Any]

Analyze impact of missing plots.

Source code in src/industrialstats/designs/rcbd.py

def missing_plot_analysis(
    self, missing_positions: list[tuple[str, str]]
) -> dict[str, Any]:
    """Analyze impact of missing plots."""
    design_matrix = self.design_matrix
    if design_matrix is None:
        design_matrix = self.generate_design()

    dm = design_matrix.copy()
    for block, treatment in missing_positions:
        mask = (dm[self.blocking_factor] == block) & (dm["Treatment"] == treatment)
        dm = dm.loc[~mask]

    balanced = (
        dm.groupby(self.blocking_factor)["Treatment"].nunique().nunique() == 1
    )
    return {
        "missing_count": len(missing_positions),
        "remaining_runs": len(dm),
        "balanced_after_missing": balanced,
    }

latin_square_option

latin_square_option() -> DataFrame | None

Generate Latin Square if conditions allow.

Source code in src/industrialstats/designs/rcbd.py

def latin_square_option(self) -> pd.DataFrame | None:
    """Generate Latin Square if conditions allow."""
    if len(self.treatments) != len(self.blocks) or len(self.treatments) < 3:
        return None

    treatments = self.treatments
    blocks = self.blocks
    n = len(treatments)
    ls_rows = []
    for i, row in enumerate(blocks):
        for j in range(n):
            treatment = treatments[(i + j) % n]
            ls_rows.append({"Row": row, "Column": j + 1, "Treatment": treatment})
    return pd.DataFrame(ls_rows)

Screening

industrialstats.designs.screening

Screening designs such as Plackett-Burman and definitive screening.

PlackettBurmanDesign

PlackettBurmanDesign(factors: list[Factor], randomize: bool = True, seed: int | None = None)

Bases: ExperimentalDesign

Plackett--Burman screening design for two-level factors.

The implemented Hadamard catalogue contains Sylvester orders and the standard 12- and 20-run Plackett--Burman base designs together with all powers-of-two doublings of those bases. For k factors, the constructor selects the smallest supported run size N satisfying N > k.

Parameters:

Name	Type	Description	Default
factors	list of Factor	Factors to include in the design. Each factor must have two levels.	required
randomize	bool	If True, randomize the run order. Defaults to True.	True
seed	int	Random seed for deterministic run-order shuffling.	None

Source code in src/industrialstats/designs/screening.py

def __init__(
    self, factors: list[Factor], randomize: bool = True, seed: int | None = None
) -> None:
    super().__init__("Plackett-Burman Design")
    self.factors = factors
    self.randomize_flag = randomize
    self.seed = seed

    if len(self.factors) < 2:
        raise ValueError("At least two factors are required")
    if not all(len(f.levels) == 2 for f in self.factors):
        raise ValueError("Plackett-Burman design requires 2-level factors")

is_supported_run_size classmethod

is_supported_run_size(n_runs: int) -> bool

Return whether n_runs belongs to the implemented PB catalogue.

Source code in src/industrialstats/designs/screening.py

@classmethod
def is_supported_run_size(cls, n_runs: int) -> bool:
    """Return whether ``n_runs`` belongs to the implemented PB catalogue."""
    if isinstance(n_runs, bool) or not isinstance(n_runs, int) or n_runs < 4:
        return False
    return any(
        n_runs % base == 0 and cls._is_power_of_two(n_runs // base)
        for base in cls._BASE_ORDERS
    )

supported_run_sizes classmethod

supported_run_sizes(max_runs: int) -> tuple[int, ...]

Return implemented Plackett--Burman run sizes up to max_runs.

Source code in src/industrialstats/designs/screening.py

@classmethod
def supported_run_sizes(cls, max_runs: int) -> tuple[int, ...]:
    """Return implemented Plackett--Burman run sizes up to ``max_runs``."""
    if isinstance(max_runs, bool) or not isinstance(max_runs, int):
        raise ValueError("max_runs must be an integer >= 4")
    if max_runs < 4:
        raise ValueError("max_runs must be >= 4")

    sizes: set[int] = set()
    for base in cls._BASE_ORDERS:
        order = base
        while order < 4:
            order *= 2
        while order <= max_runs:
            sizes.add(order)
            order *= 2
    return tuple(sorted(sizes))

run_size_for_factors classmethod

run_size_for_factors(n_factors: int) -> int

Return the smallest implemented run size able to hold n_factors.

Source code in src/industrialstats/designs/screening.py

@classmethod
def run_size_for_factors(cls, n_factors: int) -> int:
    """Return the smallest implemented run size able to hold ``n_factors``."""
    if isinstance(n_factors, bool) or not isinstance(n_factors, int):
        raise ValueError("n_factors must be an integer >= 2")
    if n_factors < 2:
        raise ValueError("n_factors must be >= 2")

    upper = 4
    while upper <= n_factors:
        upper *= 2
    candidates = cls.supported_run_sizes(upper)
    return min(order for order in candidates if order > n_factors)

run_size

run_size() -> int

Return the number of runs generated for this factor set.

Source code in src/industrialstats/designs/screening.py

def run_size(self) -> int:
    """Return the number of runs generated for this factor set."""
    return self.run_size_for_factors(len(self.factors))

generate_design

generate_design() -> DataFrame

Generate the design matrix.

Source code in src/industrialstats/designs/screening.py

def generate_design(self) -> pd.DataFrame:
    """Generate the design matrix."""
    design_matrix = self._pb_matrix(len(self.factors))
    df = pd.DataFrame(design_matrix, columns=[f.name for f in self.factors])

    self.design_matrix = df
    if self.randomize_flag:
        self.randomize(seed=self.seed)
    else:
        self.design_matrix.insert(
            0, "RunOrder", range(1, len(self.design_matrix) + 1)
        )
    return self.design_matrix

foldover

foldover() -> DataFrame

Create a foldover design to de-alias main effects.

Returns:

Type	Description
DataFrame	Foldover design matrix appended to the existing design.

Source code in src/industrialstats/designs/screening.py

def foldover(self) -> pd.DataFrame:
    """Create a foldover design to de-alias main effects.

    Returns
    -------
    pandas.DataFrame
        Foldover design matrix appended to the existing design.
    """
    design_matrix = self.design_matrix
    if design_matrix is None:
        design_matrix = self.generate_design()

    fold_df = design_matrix.copy()
    for col in self.factors:
        fold_df[col.name] = -fold_df[col.name]

    fold_df["RunOrder"] = range(len(design_matrix) + 1, 2 * len(design_matrix) + 1)
    self.design_matrix = pd.concat([design_matrix, fold_df], ignore_index=True)
    return fold_df

validate_design

validate_design() -> bool

Validate the design parameters and main-effect orthogonality.

Source code in src/industrialstats/designs/screening.py

def validate_design(self) -> bool:
    """Validate the design parameters and main-effect orthogonality."""
    if len(self.factors) < 2 or not all(len(f.levels) == 2 for f in self.factors):
        return False
    matrix = self._pb_matrix(len(self.factors)).astype(float)
    cross_product = matrix.T @ matrix
    return np.allclose(cross_product, len(matrix) * np.eye(len(self.factors)))

DefinitiveScreeningDesign

DefinitiveScreeningDesign(factors: list[Factor], randomize: bool = True, seed: int | None = None)

Bases: ExperimentalDesign

Conference-matrix definitive screening design for three-level factors.

This implementation follows the conference-matrix construction of Xiao, Lin, and Bai (2012). It uses the smallest Paley conference order large enough for the requested number of factors, takes the required columns, appends their foldover, and adds one center run.

For m factors, let q be the smallest supported value, either 1 or an odd prime, satisfying q + 1 >= m. The design then has 2 * (q + 1) + 1 runs. Consequently, designs are minimal 2m + 1 constructions when m = q + 1 and may contain additional runs when a larger conference order is required.

Factor columns are returned in coded levels -1, 0, 1. Three-level factors are interpreted as quantitative/continuous for the statistical construction even when legacy Factor metadata leaves factor_type at its default value. Mixed continuous/two-level categorical DSDs are not yet implemented.

Parameters:

Name	Type	Description	Default
factors	list of Factor	Factors to include; each must have exactly three levels.	required
randomize	bool	If True, shuffle the run order. Defaults to True.	True
seed	int	Random seed controlling the shuffle.	None

References

.. [1] Jones, B., Nachtsheim, C. J. (2011). A Class of Three-Level Designs for Definitive Screening in the Presence of Second-Order Effects. Journal of Quality Technology, 43(1), 1-15. .. [2] Xiao, L., Lin, D. K. J., Bai, F. (2012). Constructing Definitive Screening Designs Using Conference Matrices. Journal of Quality Technology, 44(1), 2-8.

Source code in src/industrialstats/designs/screening.py

def __init__(
    self, factors: list[Factor], randomize: bool = True, seed: int | None = None
) -> None:
    super().__init__("Definitive Screening Design")
    self.factors = factors
    self.randomize_flag = randomize
    self.seed = seed

    if len(self.factors) < 2:
        raise ValueError("At least two factors are required")
    if not all(len(f.levels) == 3 for f in self.factors):
        raise ValueError("Definitive screening requires 3-level factors")

generate_design

generate_design() -> DataFrame

Generate the coded definitive screening design matrix.

Source code in src/industrialstats/designs/screening.py

def generate_design(self) -> pd.DataFrame:
    """Generate the coded definitive screening design matrix."""
    coded = self._coded_matrix(len(self.factors))
    df = pd.DataFrame(coded, columns=[factor.name for factor in self.factors])

    self.design_matrix = df
    if self.randomize_flag:
        self.randomize(seed=self.seed)
    else:
        self.design_matrix.insert(
            0, "RunOrder", range(1, len(self.design_matrix) + 1)
        )
    return self.design_matrix

validate_design

validate_design() -> bool

Validate both inputs and defining algebraic DSD properties.

Source code in src/industrialstats/designs/screening.py

def validate_design(self) -> bool:
    """Validate both inputs and defining algebraic DSD properties."""
    if len(self.factors) < 2:
        return False
    if not all(len(f.levels) == 3 for f in self.factors):
        return False

    coded = self._coded_matrix(len(self.factors)).astype(float)
    n_runs, n_factors = coded.shape

    # Linear main effects are mutually orthogonal.
    cross_product = coded.T @ coded
    diagonal = np.diag(np.diag(cross_product))
    if not np.allclose(cross_product, diagonal):
        return False

    # Foldover symmetry makes linear effects orthogonal to pure quadratics.
    quadratics = coded**2
    if not np.allclose(coded.T @ quadratics, 0.0):
        return False

    # Main effects must be orthogonal to every two-factor interaction.
    interactions = np.column_stack(
        [
            coded[:, left] * coded[:, right]
            for left in range(n_factors)
            for right in range(left + 1, n_factors)
        ]
    )
    if interactions.size and not np.allclose(coded.T @ interactions, 0.0):
        return False

    # Intercept + all linear + all pure quadratic terms are estimable.
    second_order_main_model = np.column_stack((np.ones(n_runs), coded, quadratics))
    return np.linalg.matrix_rank(second_order_main_model) == 1 + 2 * n_factors

Response surface

industrialstats.designs.response_surface

Response surface methodology (RSM) designs.

This module implements the :class:ResponseSurfaceDesign class for constructing central composite and Box–Behnken designs used in process optimization.

ResponseSurfaceDesign

ResponseSurfaceDesign(factors: list[Factor], design_type: str = 'CCD', alpha: float | None = None, center_points: int = 5)

Bases: ExperimentalDesign

Response surface design for optimization studies.

Parameters:

Name	Type	Description	Default
factors	list[Factor]	Continuous factors (must be 2-level for coding).	required
design_type	str	"CCD" for Central Composite or "BBD" for Box-Behnken. Defaults to "CCD".	'CCD'
alpha	float	Alpha value for axial points (CCD only). If None, calculated for rotatability.	None
center_points	int	Number of center point replicates. Defaults to 5.	5

Examples:

Generate a Box-Behnken design with three factors and one center point::

>>> factors = [
...     Factor("A", [-1, 1]),
...     Factor("B", [-1, 1]),
...     Factor("C", [-1, 1]),
... ]
>>> design = ResponseSurfaceDesign(factors, design_type="BBD", center_points=1)
>>> design.generate_design().shape[0]
13

Initialize response surface design.

Parameters:

Name	Type	Description	Default
factors	list[Factor]	Continuous factors (must be 2-level for coding).	required
design_type	str	"CCD" for Central Composite or "BBD" for Box-Behnken. Defaults to "CCD".	'CCD'
alpha	float	Alpha value for axial points (CCD only). If None, calculated for rotatability.	None
center_points	int	Number of center point replicates. Defaults to 5.	5

Source code in src/industrialstats/designs/response_surface.py

def __init__(
    self,
    factors: list[Factor],
    design_type: str = "CCD",
    alpha: float | None = None,
    center_points: int = 5,
) -> None:
    """Initialize response surface design.

    Parameters
    ----------
    factors : list[Factor]
        Continuous factors (must be 2-level for coding).
    design_type : str, optional
        ``"CCD"`` for Central Composite or ``"BBD"`` for Box-Behnken. Defaults to ``"CCD"``.
    alpha : float, optional
        Alpha value for axial points (CCD only). If ``None``, calculated for rotatability.
    center_points : int, optional
        Number of center point replicates. Defaults to ``5``.
    """
    super().__init__(f"{design_type} Response Surface Design")

    if design_type not in ["CCD", "BBD"]:
        raise ValueError("design_type must be 'CCD' or 'BBD'")

    # Validate factors for RSM
    for factor in factors:
        if factor.factor_type != "continuous":
            raise ValueError(
                "Response surface designs require continuous factors only"
            )
        if len(factor.levels) != 2:
            raise ValueError("Factors must have exactly 2 levels for coding")

    self.factors = factors
    self.design_type = design_type
    self.center_points = center_points
    self.alpha = alpha

    if center_points < 1:
        raise ValueError("Must have at least 1 center point")

generate_design

generate_design() -> DataFrame

Generate response surface design matrix.

Source code in src/industrialstats/designs/response_surface.py

def generate_design(self) -> pd.DataFrame:
    """Generate response surface design matrix."""
    if not self.validate_design():
        raise ValueError("Invalid design configuration")

    if self.design_type == "CCD":
        return self._generate_ccd()
    if self.design_type == "BBD":
        return self._generate_bbd()
    raise ValueError(
        f"Unsupported design_type {self.design_type!r}; expected 'CCD' or 'BBD'"
    )

validate_design

validate_design() -> bool

Validate response surface design parameters.

Source code in src/industrialstats/designs/response_surface.py

def validate_design(self) -> bool:
    """Validate response surface design parameters."""
    if len(self.factors) < 2:
        return False

    if self.design_type == "BBD" and len(self.factors) < 3:
        return False

    for factor in self.factors:
        if factor.factor_type != "continuous":
            return False
        if len(factor.levels) != 2:
            return False

    return not self.center_points < 1

n_runs

n_runs() -> int

Calculate total number of runs.

Source code in src/industrialstats/designs/response_surface.py

def n_runs(self) -> int:
    """Calculate total number of runs."""
    k = len(self.factors)

    if self.design_type == "CCD":
        factorial_runs = 2**k
        axial_runs = 2 * k
        return factorial_runs + axial_runs + self.center_points
    if self.design_type == "BBD":
        # Box-Behnken: 2 * k * (k-1) + center points
        bbd_runs = 2 * k * (k - 1)
        return bbd_runs + self.center_points
    raise ValueError(
        f"Unsupported design_type {self.design_type!r}; expected 'CCD' or 'BBD'"
    )

design_properties

design_properties() -> dict[str, Any]

Calculate design properties (rotatability, orthogonality, etc.).

Source code in src/industrialstats/designs/response_surface.py

def design_properties(self) -> dict[str, Any]:
    """Calculate design properties (rotatability, orthogonality, etc.)."""
    if self.design_matrix is None:
        raise ValueError("Design not generated yet")

    k = len(self.factors)
    properties = {}

    # Basic properties
    properties["n_factors"] = k
    properties["n_runs"] = len(self.design_matrix)
    properties["design_type"] = self.design_type

    if self.design_type == "CCD":
        properties["alpha"] = self.alpha
        properties["rotatable"] = abs(self.alpha - (2**k) ** (1 / 4)) < 0.001

        # Efficiency calculations
        factorial_runs = 2**k
        axial_runs = 2 * k
        center_runs = self.center_points

        properties["factorial_fraction"] = factorial_runs / self.n_runs()
        properties["axial_fraction"] = axial_runs / self.n_runs()
        properties["center_fraction"] = center_runs / self.n_runs()

    elif self.design_type == "BBD":
        coded = self._get_design_matrix_coded()
        xtx = coded.T @ coded
        off_diag = xtx - np.diag(np.diag(xtx))
        properties["orthogonal"] = bool(np.allclose(off_diag, 0))
        properties["rotatable"] = False  # Box-Behnken is not rotatable

    return properties

prediction_variance

prediction_variance(prediction_points: list[list[float]]) -> ndarray

Calculate prediction variance at specified points.

Parameters:

Name	Type	Description	Default
prediction_points	list[list[float]]	Points in coded units where prediction variance is calculated.	required

Returns:

Type	Description
ndarray	Prediction variances.

Source code in src/industrialstats/designs/response_surface.py

def prediction_variance(self, prediction_points: list[list[float]]) -> np.ndarray:
    """Calculate prediction variance at specified points.

    Parameters
    ----------
    prediction_points : list[list[float]]
        Points in coded units where prediction variance is calculated.

    Returns
    -------
    np.ndarray
        Prediction variances.
    """
    if self.design_matrix is None:
        raise ValueError("Design not generated yet")

    # Convert design matrix to coded units for calculation
    X_coded = self._get_design_matrix_coded()

    # Add intercept column
    X_coded = np.column_stack([np.ones(len(X_coded)), X_coded])

    # Add quadratic terms
    k = len(self.factors)
    for i in range(k):
        X_coded = np.column_stack([X_coded, X_coded[:, i + 1] ** 2])

    # Add interaction terms
    for i in range(k):
        for j in range(i + 1, k):
            X_coded = np.column_stack(
                [X_coded, X_coded[:, i + 1] * X_coded[:, j + 1]]
            )

    # Calculate (X'X)^-1
    XtX_inv = np.linalg.inv(X_coded.T @ X_coded)

    # Calculate prediction variance for each point
    variances = []
    for point in prediction_points:
        # Create expanded point vector
        x_point = [1.0]  # intercept
        x_point.extend(point)  # linear terms
        x_point.extend([xi**2 for xi in point])  # quadratic terms

        # interaction terms
        for i in range(len(point)):
            for j in range(i + 1, len(point)):
                x_point.append(point[i] * point[j])

        x_point = np.array(x_point)
        variance = x_point.T @ XtX_inv @ x_point
        variances.append(variance)

    return np.array(variances)

response_surface_analysis

response_surface_analysis(response_data: list[float]) -> dict[str, Any]

Fit response surface model and analyze results.

Parameters:

Name	Type	Description	Default
response_data	list[float]	Response values for each design point.	required

Returns:

Type	Description
dict[str, Any]	Analysis results including coefficients, model fit, and optimum.

Source code in src/industrialstats/designs/response_surface.py

def response_surface_analysis(self, response_data: list[float]) -> dict[str, Any]:
    """Fit response surface model and analyze results.

    Parameters
    ----------
    response_data : list[float]
        Response values for each design point.

    Returns
    -------
    dict[str, Any]
        Analysis results including coefficients, model fit, and optimum.
    """
    if len(response_data) != len(self.design_matrix):
        raise ValueError("Response data length must match design matrix")

    # Get coded design matrix
    X_coded = self._get_design_matrix_coded()
    k = len(self.factors)

    # Build full quadratic model matrix
    # Intercept
    X_model = np.ones((len(X_coded), 1))
    term_names = ["Intercept"]

    # Linear terms
    X_model = np.column_stack([X_model, X_coded])
    term_names.extend([f.name for f in self.factors])

    # Quadratic terms
    for i, factor in enumerate(self.factors):
        X_model = np.column_stack([X_model, X_coded[:, i] ** 2])
        term_names.append(f"{factor.name}²")

    # Interaction terms
    for i in range(k):
        for j in range(i + 1, k):
            X_model = np.column_stack([X_model, X_coded[:, i] * X_coded[:, j]])
            term_names.append(f"{self.factors[i].name}*{self.factors[j].name}")

    # Fit model using least squares
    y = np.array(response_data)

    try:
        # Calculate coefficients
        XtX_inv = np.linalg.inv(X_model.T @ X_model)
        coefficients = XtX_inv @ X_model.T @ y

        # Calculate fitted values and residuals
        y_fitted = X_model @ coefficients
        residuals = y - y_fitted

        # Calculate R-squared
        ss_total = np.sum((y - np.mean(y)) ** 2)
        ss_residual = np.sum(residuals**2)
        r_squared = 1 - (ss_residual / ss_total)

        # Adjusted R-squared
        n = len(y)
        p = len(coefficients) - 1  # excluding intercept
        adj_r_squared = 1 - (1 - r_squared) * (n - 1) / (n - p - 1)

        # Standard errors
        mse = ss_residual / (n - len(coefficients))
        std_errors = np.sqrt(np.diag(XtX_inv) * mse)

        # t-statistics and p-values
        t_stats = coefficients / std_errors
        p_values = 2 * (1 - stats.t.cdf(np.abs(t_stats), n - len(coefficients)))

        # Find optimum (coded units)
        optimum_coded = self._find_optimum_coded(coefficients, term_names)

        # Convert optimum to actual units
        optimum_actual = {}
        if optimum_coded is not None:
            for i, factor in enumerate(self.factors):
                low_level = factor.levels[0]
                high_level = factor.levels[1]
                center = (low_level + high_level) / 2
                half_range = (high_level - low_level) / 2

                actual_value = center + optimum_coded[i] * half_range
                optimum_actual[factor.name] = actual_value

        # Prepare results
        results = {
            "coefficients": dict(zip(term_names, coefficients, strict=True)),
            "std_errors": dict(zip(term_names, std_errors, strict=True)),
            "t_statistics": dict(zip(term_names, t_stats, strict=True)),
            "p_values": dict(zip(term_names, p_values, strict=True)),
            "r_squared": r_squared,
            "adj_r_squared": adj_r_squared,
            "rmse": np.sqrt(mse),
            "fitted_values": y_fitted,
            "residuals": residuals,
            "optimum_coded": optimum_coded,
            "optimum_actual": optimum_actual,
            "information_matrix_inv": XtX_inv,
            "degrees_of_freedom": n - len(coefficients),
            "mse": mse,
        }

        return results

    except np.linalg.LinAlgError as e:
        raise ValueError("Unable to fit model - design matrix is singular") from e

steepest_ascent

steepest_ascent(model_results: dict[str, Any], start_point: dict[str, float] | None = None, step_length: float = 0.5, n_steps: int = 10, direction: str = 'ascent', visualize: bool = False, plot_factors: tuple[str, str] | None = None) -> dict[str, Any]

Compute the steepest ascent/descent path for the fitted surface.

This routine follows the gradient-based approach described by Box and Draper [1]_ to generate an ordered sequence of points in the direction of the largest increase (or decrease) in the response. Steps are taken in coded units and converted back to actual factor levels for reporting.

Parameters:

Name	Type	Description	Default
model_results	dict[str, Any]	Output from :meth:response_surface_analysis containing the fitted coefficients.	required
start_point	dict[str, float] | None	Starting location in actual units. If omitted, the design centre is used.	None
step_length	float	Length of each step in coded units once the gradient direction is normalised. Defaults to 0.5.	0.5
n_steps	int	Number of steps to compute along the path. Defaults to 10.	10
direction	str	Direction of movement relative to the gradient. Defaults to "ascent".	'ascent'
visualize	bool	Whether to return a contour plot overlay. Defaults to False.	False
plot_factors	tuple[str, str] | None	Pair of factor names to use on the contour plot. When None the first two factors are used.	None

Returns:

Type	Description
Dict[str, Any]	Dictionary with keys "path" (DataFrame of coded and actual points) and "figure" (Plotly figure when visualize is True).

Raises:

Type	Description
ValueError	If an invalid direction is provided or the gradient is zero.

Examples:

>>> factors = [Factor("x1", [-1, 1]), Factor("x2", [-1, 1])]
>>> design = ResponseSurfaceDesign(factors)
>>> dm = design.generate_design()
>>> y = dm["x1"] * -2 + dm["x2"]
>>> results = design.response_surface_analysis(y.tolist())
>>> path = design.steepest_ascent(results, n_steps=3)
>>> list(path["path"]["Step"])
[0, 1, 2, 3]

References

Box, G. E. P., & Draper, N. R. (2007). Response Surfaces, Mixtures, and Ridge Analyses (2nd ed.). Wiley.

Source code in src/industrialstats/designs/response_surface.py

def steepest_ascent(
    self,
    model_results: dict[str, Any],
    start_point: dict[str, float] | None = None,
    step_length: float = 0.5,
    n_steps: int = 10,
    direction: str = "ascent",
    visualize: bool = False,
    plot_factors: tuple[str, str] | None = None,
) -> dict[str, Any]:
    """Compute the steepest ascent/descent path for the fitted surface.

    This routine follows the gradient-based approach described by Box and
    Draper [1]_ to generate an ordered sequence of points in the direction of
    the largest increase (or decrease) in the response. Steps are taken in
    coded units and converted back to actual factor levels for reporting.

    Parameters
    ----------
    model_results
        Output from :meth:`response_surface_analysis`
        containing the fitted coefficients.
    start_point
        Starting location in actual units. If omitted, the
        design centre is used.
    step_length
        Length of each step in coded units once the gradient
        direction is normalised. Defaults to ``0.5``.
    n_steps
        Number of steps to compute along the path. Defaults to
        ``10``.
    direction
        Direction of movement relative to the gradient. Defaults
        to ``"ascent"``.
    visualize
        Whether to return a contour plot overlay. Defaults to
        ``False``.
    plot_factors
        Pair of factor names to use on the contour plot. When
        ``None`` the first two factors are used.

    Returns
    -------
    Dict[str, Any]
        Dictionary with keys ``"path"`` (DataFrame of coded
        and actual points) and ``"figure"`` (Plotly figure when
        ``visualize`` is ``True``).

    Raises
    ------
    ValueError
        If an invalid direction is provided or the gradient is
        zero.

    Examples
    --------
    >>> factors = [Factor("x1", [-1, 1]), Factor("x2", [-1, 1])]
    >>> design = ResponseSurfaceDesign(factors)
    >>> dm = design.generate_design()
    >>> y = dm["x1"] * -2 + dm["x2"]
    >>> results = design.response_surface_analysis(y.tolist())
    >>> path = design.steepest_ascent(results, n_steps=3)
    >>> list(path["path"]["Step"])
    [0, 1, 2, 3]

    References
    ----------
    Box, G. E. P., & Draper, N. R. (2007). *Response Surfaces, Mixtures,
    and Ridge Analyses* (2nd ed.). Wiley.
    """

    if direction not in {"ascent", "descent"}:
        raise ValueError("direction must be 'ascent' or 'descent'")

    coefficient_map = model_results.get("coefficients")
    if coefficient_map is None:
        raise ValueError("model_results must include 'coefficients'")

    _, b, _ = self._quadratic_components(coefficient_map)
    gradient = b if direction == "ascent" else -b

    norm = float(np.linalg.norm(gradient))
    if norm == 0:
        raise ValueError("Gradient is zero; steepest path is undefined")

    unit_direction = gradient / norm

    if start_point is None:
        coded_start = np.zeros(len(self.factors))
    else:
        missing = set(self.factor_names).difference(start_point)
        if missing:
            missing_str = ", ".join(sorted(missing))
            raise ValueError(
                f"start_point is missing levels for factors: {missing_str}"
            )
        ordered_actual = [start_point[f.name] for f in self.factors]
        coded_start = self._actual_to_coded_vector(ordered_actual)

    coded_points = [coded_start.copy()]
    for step in range(1, n_steps + 1):
        coded_points.append(coded_start + step_length * step * unit_direction)

    actual_points = [self._coded_to_actual_vector(point) for point in coded_points]
    responses = [
        self._evaluate_quadratic(coefficient_map, point) for point in coded_points
    ]

    index_lookup = {name: idx for idx, name in enumerate(self.factor_names)}

    data = {
        "Step": list(range(len(coded_points))),
        **{
            f"coded_{factor.name}": [point[i] for point in coded_points]
            for i, factor in enumerate(self.factors)
        },
        **{
            factor.name: [point[i] for point in actual_points]
            for i, factor in enumerate(self.factors)
        },
        "predicted_response": responses,
    }
    path_df = pd.DataFrame(data)

    figure = None
    if visualize:
        chosen = (
            plot_factors
            if plot_factors is not None
            else (self.factors[0].name, self.factors[1].name)
        )
        plotter = self._quadratic_plotter(coefficient_map)
        path_points = [
            (point[index_lookup[chosen[0]]], point[index_lookup[chosen[1]]])
            for point in actual_points
        ]
        figure = plotter.contour_plot(chosen[0], chosen[1], path=path_points)

    return {"path": path_df, "figure": figure}

ridge_analysis

ridge_analysis(model_results: dict[str, Any], radii: Sequence[float], constraints: Sequence[Callable[[ndarray], float]] | None = None, penalty_weight: float = 100.0, visualize: bool = False, plot_factors: tuple[str, str] | None = None) -> dict[str, Any]

Perform ridge analysis for constrained optimisation.

Ridge analysis seeks the best point on a hypersphere of radius r in coded units by solving the Lagrangian system (B + λI) x = -0.5 b for each candidate radius [1]_. Optional inequality constraints g_i(x) ≤ 0 are enforced through quadratic penalties to discourage infeasible solutions.

Parameters:

Name	Type	Description	Default
model_results	dict[str, Any]	Output from :meth:response_surface_analysis containing the fitted coefficients.	required
radii	sequence of float	Radii in coded units at which to evaluate the ridge solution.	required
constraints	sequence of callable	Functions mapping a coded point to a real value. Positive values are treated as violations and penalised. Defaults to None.	None
penalty_weight	float	Penalty scaling factor applied to squared constraint violations. Defaults to 100.0.	100.0
visualize	bool	When True returns a contour plot of the first two factors with the ridge path overlay. Defaults to False.	False
plot_factors	tuple[str, str]	Factor names to use for the visualisation. Defaults to the first two factors.	None

Returns:

Type	Description
dict[str, Any]	Dictionary containing the ridge solutions ("solutions" DataFrame) and "figure" with the optional Plotly contour.

References

.. [1] Box, G. E. P., & Draper, N. R. (2007). Response Surfaces, Mixtures, and Ridge Analyses (2nd ed.). Wiley.

Source code in src/industrialstats/designs/response_surface.py

def ridge_analysis(
    self,
    model_results: dict[str, Any],
    radii: Sequence[float],
    constraints: Sequence[Callable[[np.ndarray], float]] | None = None,
    penalty_weight: float = 100.0,
    visualize: bool = False,
    plot_factors: tuple[str, str] | None = None,
) -> dict[str, Any]:
    """Perform ridge analysis for constrained optimisation.

    Ridge analysis seeks the best point on a hypersphere of radius ``r`` in
    coded units by solving the Lagrangian system ``(B + λI) x = -0.5 b`` for
    each candidate radius [1]_. Optional inequality constraints ``g_i(x) ≤ 0``
    are enforced through quadratic penalties to discourage infeasible
    solutions.

    Parameters
    ----------
    model_results : dict[str, Any]
        Output from :meth:`response_surface_analysis` containing the fitted
        coefficients.
    radii : sequence of float
        Radii in coded units at which to evaluate the ridge solution.
    constraints : sequence of callable, optional
        Functions mapping a coded point to a real value. Positive values are
        treated as violations and penalised. Defaults to ``None``.
    penalty_weight : float, optional
        Penalty scaling factor applied to squared constraint violations.
        Defaults to ``100.0``.
    visualize : bool, optional
        When ``True`` returns a contour plot of the first two factors with the
        ridge path overlay. Defaults to ``False``.
    plot_factors : tuple[str, str], optional
        Factor names to use for the visualisation. Defaults to the first two
        factors.

    Returns
    -------
    dict[str, Any]
        Dictionary containing the ridge solutions (``"solutions"`` DataFrame)
        and ``"figure"`` with the optional Plotly contour.

    References
    ----------
    .. [1] Box, G. E. P., & Draper, N. R. (2007). *Response Surfaces,
       Mixtures, and Ridge Analyses* (2nd ed.). Wiley.
    """

    coefficient_map = model_results.get("coefficients")
    if coefficient_map is None:
        raise ValueError("model_results must include 'coefficients'")

    intercept, b, B = self._quadratic_components(coefficient_map)
    k = len(self.factors)
    identity = np.eye(k)
    penalties = constraints or []

    def _solve_radius(radius: float) -> tuple[np.ndarray, float, float]:
        def norm_difference(lmbda: float) -> float:
            matrix = B + lmbda * identity
            solution = np.linalg.solve(matrix, -0.5 * b)
            return float(solution @ solution - radius**2)

        # Identify a bracket for the root of norm_difference.
        bracket = None
        candidates = np.linspace(-50.0, 50.0, 400)
        previous_value = None
        previous_lambda = None
        for lmbda in candidates:
            try:
                value = norm_difference(lmbda)
            except np.linalg.LinAlgError:
                previous_value = None
                previous_lambda = None
                continue
            if previous_value is not None and value * previous_value <= 0:
                bracket = (previous_lambda, lmbda)
                break
            previous_value = value
            previous_lambda = lmbda

        if bracket is None:
            # Fall back to the best candidate if no sign change is found.
            feasible_values = []
            for lmbda in candidates:
                try:
                    difference = abs(norm_difference(lmbda))
                    feasible_values.append((difference, lmbda))
                except np.linalg.LinAlgError:
                    continue
            if not feasible_values:
                raise ValueError(
                    "Unable to bracket Lagrange multiplier for ridge analysis"
                )
            lmbda = min(feasible_values, key=lambda item: item[0])[1]
        else:
            lmbda = optimize.brentq(norm_difference, *bracket, maxiter=200)

        solution = np.linalg.solve(B + lmbda * identity, -0.5 * b)
        response = intercept + solution @ b + solution.T @ B @ solution
        penalty = 0.0
        if penalties:
            for constraint in penalties:
                violation = float(constraint(solution))
                penalty += penalty_weight * max(0.0, violation) ** 2
        return solution, response, penalty

    coded_solutions = []
    responses = []
    penalties_applied = []
    for radius in radii:
        if radius <= 0:
            raise ValueError("radii must contain positive values")
        point, response, penalty = _solve_radius(radius)
        coded_solutions.append(point)
        responses.append(response)
        penalties_applied.append(penalty)

    actual_points = [
        self._coded_to_actual_vector(point) for point in coded_solutions
    ]

    results_df = pd.DataFrame(
        {
            "radius": list(radii),
            "objective": responses,
            "penalty": penalties_applied,
            **{
                f"coded_{factor.name}": [point[i] for point in coded_solutions]
                for i, factor in enumerate(self.factors)
            },
            **{
                factor.name: [point[i] for point in actual_points]
                for i, factor in enumerate(self.factors)
            },
        }
    )

    figure = None
    if visualize:
        chosen = (
            plot_factors
            if plot_factors is not None
            else (self.factors[0].name, self.factors[1].name)
        )
        index_lookup = {name: idx for idx, name in enumerate(self.factor_names)}
        plotter = self._quadratic_plotter(coefficient_map)
        path_points = [
            (point[index_lookup[chosen[0]]], point[index_lookup[chosen[1]]])
            for point in actual_points
        ]
        figure = plotter.contour_plot(chosen[0], chosen[1], path=path_points)

    return {"solutions": results_df, "figure": figure}

canonical_analysis

canonical_analysis(model_results: dict[str, Any]) -> dict[str, Any]

Carry out canonical analysis of the fitted response surface.

Canonical analysis diagonalises the quadratic form to reveal the surface curvature and nature of the stationary point (minimum, maximum, saddle, or ridge) following Box and Draper [1]_. Eigenvectors define the canonical directions while eigenvalues quantify curvature along each axis.

Parameters:

Name	Type	Description	Default
model_results	dict[str, Any]	Output from :meth:response_surface_analysis containing model coefficients and residual degrees of freedom.	required

Returns:

Type	Description
dict[str, Any]	Summary including the stationary point in coded and actual units, eigenvalues/eigenvectors, surface classification, and 95 % confidence ellipsoid axes for the optimum when degrees of freedom are available.

References

.. [1] Box, G. E. P., & Draper, N. R. (2007). Response Surfaces, Mixtures, and Ridge Analyses (2nd ed.). Wiley.

Source code in src/industrialstats/designs/response_surface.py

def canonical_analysis(self, model_results: dict[str, Any]) -> dict[str, Any]:
    """Carry out canonical analysis of the fitted response surface.

    Canonical analysis diagonalises the quadratic form to reveal the surface
    curvature and nature of the stationary point (minimum, maximum, saddle,
    or ridge) following Box and Draper [1]_. Eigenvectors define the
    canonical directions while eigenvalues quantify curvature along each
    axis.

    Parameters
    ----------
    model_results : dict[str, Any]
        Output from :meth:`response_surface_analysis` containing model
        coefficients and residual degrees of freedom.

    Returns
    -------
    dict[str, Any]
        Summary including the stationary point in coded and actual units,
        eigenvalues/eigenvectors, surface classification, and 95 % confidence
        ellipsoid axes for the optimum when degrees of freedom are available.

    References
    ----------
    .. [1] Box, G. E. P., & Draper, N. R. (2007). *Response Surfaces,
       Mixtures, and Ridge Analyses* (2nd ed.). Wiley.
    """

    coefficient_map = model_results.get("coefficients")
    if coefficient_map is None:
        raise ValueError("model_results must include 'coefficients'")

    intercept, b, B = self._quadratic_components(coefficient_map)
    try:
        stationary_coded = -0.5 * np.linalg.solve(B, b)
    except np.linalg.LinAlgError:
        stationary_coded = None

    stationary_actual = (
        self._coded_to_actual_vector(stationary_coded)
        if stationary_coded is not None
        else None
    )

    eigenvalues, eigenvectors = np.linalg.eigh(B)

    if np.all(eigenvalues > 0):
        surface_type = "minimum"
    elif np.all(eigenvalues < 0):
        surface_type = "maximum"
    elif np.any(np.isclose(eigenvalues, 0)):
        surface_type = "ridge"
    else:
        surface_type = "saddle"

    df_resid = model_results.get("degrees_of_freedom")
    confidence = None
    if df_resid is not None and df_resid > 0:
        k = len(self.factors)
        f_value = stats.f.ppf(0.95, k, df_resid)
        axes = []
        for value in eigenvalues:
            if np.isclose(value, 0):
                axes.append(np.inf)
            else:
                axes.append(float(np.sqrt(f_value / abs(value))))
        confidence = {
            "alpha": 0.95,
            "axes": axes,
            "eigenvectors": eigenvectors,
        }

    stationary_response = None
    if stationary_coded is not None:
        stationary_response = self._evaluate_quadratic(
            coefficient_map, stationary_coded
        )

    return {
        "intercept": intercept,
        "stationary_point_coded": stationary_coded,
        "stationary_point_actual": stationary_actual,
        "stationary_response": stationary_response,
        "eigenvalues": eigenvalues,
        "eigenvectors": eigenvectors,
        "surface_type": surface_type,
        "confidence_region": confidence,
    }

multiple_response_optimization

multiple_response_optimization(response_models: dict[str, dict[str, Any]], weights: dict[str, float] | None = None, desirability_functions: dict[str, Callable[[float], float]] | None = None, constraint_functions: Sequence[Callable[[ndarray], float]] | None = None, grid_resolution: int = 25, search_radius: float = 1.5, penalty_weight: float = 50.0, weight_perturbation: float = 0.15) -> dict[str, Any]

Simultaneously optimise multiple responses.

The optimisation proceeds by evaluating fitted response models on a lattice in the coded factor space, computing desirability functions for each response, and aggregating them using weighted geometric means. A Pareto frontier of feasible points is additionally identified, and weight sensitivity analysis perturbs the supplied weights to study robustness.

Parameters:

Name	Type	Description	Default
response_models	dict[str, dict[str, Any]]	Mapping from response name to model results produced by :meth:response_surface_analysis. Each entry may include a "goal" field ("max" or "min").	required
weights	dict[str, float]	Importance weights for each response. If omitted, equal weights are assigned.	None
desirability_functions	dict[str, callable]	Custom desirability functions mapping response values to [0, 1]. Defaults to piecewise-linear ramps based on the response range.	None
constraint_functions	sequence of callable	Inequality constraints g_i(x) evaluated in coded space. Positive values are penalised quadratically.	None
grid_resolution	int	Number of grid points per factor in the coded space. Defaults to 25.	25
search_radius	float	Extent of the coded search space ([-radius, radius] per factor). Defaults to 1.5.	1.5
penalty_weight	float	Penalty multiplier for constraint violations. Defaults to 50.0.	50.0
weight_perturbation	float	Relative perturbation applied to weights during the sensitivity analysis. Defaults to 0.15.	0.15

Returns:

Type	Description
dict[str, Any]	Dictionary containing the best compromise solution, Pareto frontier, and weight sensitivity study.

References

.. [1] Box, G. E. P., & Draper, N. R. (2007). Response Surfaces, Mixtures, and Ridge Analyses (2nd ed.). Wiley.

Source code in src/industrialstats/designs/response_surface.py

def multiple_response_optimization(
    self,
    response_models: dict[str, dict[str, Any]],
    weights: dict[str, float] | None = None,
    desirability_functions: dict[str, Callable[[float], float]] | None = None,
    constraint_functions: Sequence[Callable[[np.ndarray], float]] | None = None,
    grid_resolution: int = 25,
    search_radius: float = 1.5,
    penalty_weight: float = 50.0,
    weight_perturbation: float = 0.15,
) -> dict[str, Any]:
    """Simultaneously optimise multiple responses.

    The optimisation proceeds by evaluating fitted response models on a
    lattice in the coded factor space, computing desirability functions for
    each response, and aggregating them using weighted geometric means. A
    Pareto frontier of feasible points is additionally identified, and
    weight sensitivity analysis perturbs the supplied weights to study
    robustness.

    Parameters
    ----------
    response_models : dict[str, dict[str, Any]]
        Mapping from response name to model results produced by
        :meth:`response_surface_analysis`. Each entry may include a
        ``"goal"`` field (``"max"`` or ``"min"``).
    weights : dict[str, float], optional
        Importance weights for each response. If omitted, equal weights are
        assigned.
    desirability_functions : dict[str, callable], optional
        Custom desirability functions mapping response values to ``[0, 1]``.
        Defaults to piecewise-linear ramps based on the response range.
    constraint_functions : sequence of callable, optional
        Inequality constraints ``g_i(x)`` evaluated in coded space. Positive
        values are penalised quadratically.
    grid_resolution : int, optional
        Number of grid points per factor in the coded space. Defaults to
        ``25``.
    search_radius : float, optional
        Extent of the coded search space (``[-radius, radius]`` per factor).
        Defaults to ``1.5``.
    penalty_weight : float, optional
        Penalty multiplier for constraint violations. Defaults to ``50.0``.
    weight_perturbation : float, optional
        Relative perturbation applied to weights during the sensitivity
        analysis. Defaults to ``0.15``.

    Returns
    -------
    dict[str, Any]
        Dictionary containing the best compromise solution, Pareto frontier,
        and weight sensitivity study.

    References
    ----------
    .. [1] Box, G. E. P., & Draper, N. R. (2007). *Response Surfaces,
       Mixtures, and Ridge Analyses* (2nd ed.). Wiley.
    """

    if not response_models:
        raise ValueError("response_models must not be empty")

    desirability_functions = desirability_functions or {}
    response_names = list(response_models)

    if weights is None:
        weights = dict.fromkeys(response_names, 1.0)
    missing_weights = set(response_names).difference(weights)
    if missing_weights:
        for name in missing_weights:
            weights[name] = 1.0

    weight_sum = float(sum(weights.values()))
    if weight_sum <= 0:
        raise ValueError("weights must sum to a positive value")
    normalised_weights = {k: v / weight_sum for k, v in weights.items()}

    axes = [
        np.linspace(-search_radius, search_radius, grid_resolution)
        for _ in self.factors
    ]
    grid_points = np.array(list(product(*axes)))

    predictions: dict[str, np.ndarray] = {}
    goals: dict[str, str] = {}
    for name, result in response_models.items():
        coeffs = result.get("coefficients")
        if coeffs is None:
            raise ValueError(
                f"response model '{name}' must include 'coefficients' from response_surface_analysis"
            )
        goals[name] = result.get("goal", "max").lower()
        preds = np.array(
            [self._evaluate_quadratic(coeffs, point) for point in grid_points]
        )
        if goals[name] not in {"max", "min"}:
            raise ValueError("goal must be 'max' or 'min'")
        predictions[name] = preds

    desirabilities: dict[str, np.ndarray] = {}
    for name in response_names:
        if name in desirability_functions:
            func = desirability_functions[name]
            desirabilities[name] = np.array(
                [func(value) for value in predictions[name]]
            )
        else:
            values = predictions[name]
            v_min = float(np.min(values))
            v_max = float(np.max(values))
            if np.isclose(v_max, v_min):
                desirabilities[name] = np.ones_like(values)
            else:
                if goals[name] == "max":
                    desirabilities[name] = np.clip(
                        (values - v_min) / (v_max - v_min), 0.0, 1.0
                    )
                else:
                    desirabilities[name] = np.clip(
                        (v_max - values) / (v_max - v_min), 0.0, 1.0
                    )

    penalties = np.zeros(len(grid_points), dtype=float)
    if constraint_functions:
        for idx, point in enumerate(grid_points):
            violation = 0.0
            for constraint in constraint_functions:
                violation_value = float(constraint(point))
                violation += max(0.0, violation_value) ** 2
            penalties[idx] = penalty_weight * violation

    def combined_desirability(weight_map: dict[str, float]) -> np.ndarray:
        overall = np.ones(len(grid_points), dtype=float)
        for name in response_names:
            overall *= desirabilities[name] ** weight_map[name]
        return overall

    base_scores = combined_desirability(normalised_weights) - penalties
    best_index = int(np.argmax(base_scores))

    best_coded = grid_points[best_index]
    best_actual = self._coded_to_actual_vector(best_coded)
    best_predictions = {
        name: predictions[name][best_index] for name in response_names
    }
    best_desirabilities = {
        name: desirabilities[name][best_index] for name in response_names
    }

    feasible_indices = (
        np.where(penalties <= 1e-8)[0]
        if constraint_functions
        else np.arange(len(grid_points))
    )
    pareto_indices: list[int] = []
    if feasible_indices.size > 0:
        signed = []
        for name in response_names:
            if goals[name] == "max":
                signed.append(predictions[name][feasible_indices])
            else:
                signed.append(-predictions[name][feasible_indices])
        signed_matrix = np.column_stack(signed)
        for idx, candidate in enumerate(signed_matrix):
            dominated = False
            for other_idx, other in enumerate(signed_matrix):
                if other_idx == idx:
                    continue
                if np.all(other >= candidate) and np.any(other > candidate):
                    dominated = True
                    break
            if not dominated:
                pareto_indices.append(int(feasible_indices[idx]))

    pareto_points = []
    for idx in pareto_indices:
        coded = grid_points[idx]
        actual = self._coded_to_actual_vector(coded)
        pareto_points.append(
            {
                "coded": coded,
                "actual": actual,
                "responses": {
                    name: predictions[name][idx] for name in response_names
                },
            }
        )

    def _weight_analysis(new_weights: dict[str, float]) -> dict[str, Any]:
        total = float(sum(new_weights.values()))
        if total <= 0:
            return {"weights": new_weights, "best_index": None}
        norm = {k: v / total for k, v in new_weights.items()}
        scores = combined_desirability(norm) - penalties
        idx = int(np.argmax(scores))
        return {
            "weights": norm,
            "coded_point": grid_points[idx],
            "actual_point": self._coded_to_actual_vector(grid_points[idx]),
            "responses": {name: predictions[name][idx] for name in response_names},
            "overall_desirability": float(scores[idx]),
        }

    sensitivity: list[dict[str, Any]] = []
    for name in response_names:
        base = dict(weights)
        base[name] *= 1 + weight_perturbation
        sensitivity.append(_weight_analysis(base))
        base = dict(weights)
        base[name] *= max(1 - weight_perturbation, 0)
        sensitivity.append(_weight_analysis(base))

    grid_df = pd.DataFrame(
        {
            **{
                f"coded_{factor.name}": grid_points[:, idx]
                for idx, factor in enumerate(self.factors)
            },
            **{f"pred_{name}": predictions[name] for name in response_names},
            "overall_desirability": combined_desirability(normalised_weights),
            "penalty": penalties,
        }
    )

    return {
        "optimum": {
            "coded": best_coded,
            "actual": best_actual,
            "responses": best_predictions,
            "desirabilities": best_desirabilities,
            "overall_desirability": float(base_scores[best_index]),
        },
        "pareto_frontier": pareto_points,
        "weight_sensitivity": sensitivity,
        "grid": grid_df,
    }

contour_data

contour_data(coefficients: dict[str, float], factor1: str, factor2: str, grid_size: int = 20) -> tuple[ndarray, ndarray, ndarray]

Generate contour plot data for two factors.

Parameters:

Name	Type	Description	Default
coefficients	dict[str, float]	Model coefficients from :func:response_surface_analysis.	required
factor1	str	Name of the first factor for the plot.	required
factor2	str	Name of the second factor for the plot.	required
grid_size	int	Grid resolution. Defaults to 20.	20

Returns:

Type	Description
tuple[ndarray, ndarray, ndarray]	X, Y, Z arrays for contour plotting.

Source code in src/industrialstats/designs/response_surface.py

def contour_data(
    self,
    coefficients: dict[str, float],
    factor1: str,
    factor2: str,
    grid_size: int = 20,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Generate contour plot data for two factors.

    Parameters
    ----------
    coefficients : dict[str, float]
        Model coefficients from :func:`response_surface_analysis`.
    factor1 : str
        Name of the first factor for the plot.
    factor2 : str
        Name of the second factor for the plot.
    grid_size : int, optional
        Grid resolution. Defaults to ``20``.

    Returns
    -------
    tuple[np.ndarray, np.ndarray, np.ndarray]
        ``X``, ``Y``, ``Z`` arrays for contour plotting.
    """
    # Find factor indices
    factor1_idx = next(i for i, f in enumerate(self.factors) if f.name == factor1)
    factor2_idx = next(i for i, f in enumerate(self.factors) if f.name == factor2)

    # Create grid in coded units
    x_coded = np.linspace(-2, 2, grid_size)
    y_coded = np.linspace(-2, 2, grid_size)
    X_coded, Y_coded = np.meshgrid(x_coded, y_coded)

    # Calculate response surface
    Z = np.zeros_like(X_coded)

    for i in range(grid_size):
        for j in range(grid_size):
            # Set other factors to center (0)
            point = np.zeros(len(self.factors))
            point[factor1_idx] = X_coded[i, j]
            point[factor2_idx] = Y_coded[i, j]

            # Calculate response using model
            response = coefficients["Intercept"]

            # Linear terms
            for k, factor in enumerate(self.factors):
                response += coefficients.get(factor.name, 0) * point[k]

            # Quadratic terms
            for k, factor in enumerate(self.factors):
                response += coefficients.get(f"{factor.name}²", 0) * point[k] ** 2

            # Interaction terms
            for k in range(len(self.factors)):
                for m in range(k + 1, len(self.factors)):
                    interaction_coef = coefficients.get(
                        f"{self.factors[k].name}*{self.factors[m].name}", 0
                    )
                    response += interaction_coef * point[k] * point[m]

            Z[i, j] = response

    # Convert grid to actual units
    factor1_obj = next(f for f in self.factors if f.name == factor1)
    factor2_obj = next(f for f in self.factors if f.name == factor2)

    # Factor 1 conversion
    center1 = (factor1_obj.levels[0] + factor1_obj.levels[1]) / 2
    range1 = (factor1_obj.levels[1] - factor1_obj.levels[0]) / 2
    X_actual = center1 + X_coded * range1

    # Factor 2 conversion
    center2 = (factor2_obj.levels[0] + factor2_obj.levels[1]) / 2
    range2 = (factor2_obj.levels[1] - factor2_obj.levels[0]) / 2
    Y_actual = center2 + Y_coded * range2

    return X_actual, Y_actual, Z

Optimal designs

industrialstats.designs.optimal

Optimal experimental designs using algorithmic approaches.

OptimalDesign

OptimalDesign(factors: list[Factor], n_runs: int, criterion: str = 'D', model_terms: list[str] | None = None)

Bases: ExperimentalDesign

Generate optimal experimental designs using exchange algorithms.

Supports D-optimal, A-optimal, G-optimal, and I-optimal criteria.

Initialize optimal design.

Parameters:

Name	Type	Description	Default
factors	list[Factor]	Experimental factors.	required
n_runs	int	Number of experimental runs.	required
criterion	str	Optimality criterion ("D", "A", "G", or "I"). Defaults to "D".	'D'
model_terms	list[str]	Model terms to include. Defaults to main effects and interactions.	None

Source code in src/industrialstats/designs/optimal.py

def __init__(
    self,
    factors: list[Factor],
    n_runs: int,
    criterion: str = "D",
    model_terms: list[str] | None = None,
) -> None:
    """Initialize optimal design.

    Parameters
    ----------
    factors : list[Factor]
        Experimental factors.
    n_runs : int
        Number of experimental runs.
    criterion : str, optional
        Optimality criterion (``"D"``, ``"A"``, ``"G"``, or ``"I"``). Defaults to ``"D"``.
    model_terms : list[str], optional
        Model terms to include. Defaults to main effects and interactions.
    """
    super().__init__(f"{criterion}-Optimal Design")

    if criterion not in ["D", "A", "G", "I"]:
        raise ValueError("criterion must be 'D', 'A', 'G', or 'I'")

    if n_runs < len(factors) + 1:
        raise ValueError("n_runs must be greater than number of factors")

    self.factors = factors
    self.n_runs = n_runs
    self.criterion = criterion
    self.model_terms = model_terms or self._default_model_terms()
    self.candidate_set: pd.DataFrame | None = None
    self.candidate_model_matrix: np.ndarray | None = None
    self.exchange_history: list[dict[str, Any]] = []

generate_candidate_set

generate_candidate_set(grid_density: int = 5) -> DataFrame

Generate candidate set of all possible design points.

Parameters:

Name	Type	Description	Default
grid_density	int	Number of levels for continuous factors. Defaults to 5.	5

Returns:

Type	Description
DataFrame	Candidate set of design points.

Source code in src/industrialstats/designs/optimal.py

def generate_candidate_set(self, grid_density: int = 5) -> pd.DataFrame:
    """Generate candidate set of all possible design points.

    Parameters
    ----------
    grid_density : int, optional
        Number of levels for continuous factors. Defaults to ``5``.

    Returns
    -------
    pd.DataFrame
        Candidate set of design points.
    """
    candidate_points = []

    # Generate levels for each factor
    factor_levels = []
    for factor in self.factors:
        if factor.factor_type == "categorical":
            factor_levels.append(factor.levels)
        else:
            # Create grid for continuous factors
            min_val = min(factor.levels)
            max_val = max(factor.levels)
            levels = np.linspace(min_val, max_val, grid_density)
            factor_levels.append(levels.tolist())

    # Generate all combinations
    from itertools import product

    for i, combination in enumerate(product(*factor_levels)):
        point = {"CandidateID": i + 1}
        for j, factor in enumerate(self.factors):
            point[factor.name] = combination[j]
        candidate_points.append(point)

    self.candidate_set = pd.DataFrame(candidate_points)
    return self.candidate_set

generate_design

generate_design(max_iterations: int = 1000, random_start: bool = True, n_random_starts: int = 5, improvement_threshold: float = 1e-06) -> DataFrame

Generate optimal design using coordinate exchange algorithm.

Parameters:

Name	Type	Description	Default
max_iterations	int	Maximum number of exchange iterations. Defaults to 1000.	1000
random_start	bool	Whether to use random starting design. Defaults to True.	True
n_random_starts	int	Number of random starts to try. Defaults to 5.	5
improvement_threshold	float	Minimum improvement in the criterion required to continue iterations. Defaults to 1e-6.	1e-06

Returns:

Type	Description
DataFrame	Optimal design matrix.

Source code in src/industrialstats/designs/optimal.py

def generate_design(
    self,
    max_iterations: int = 1000,
    random_start: bool = True,
    n_random_starts: int = 5,
    improvement_threshold: float = 1e-6,
) -> pd.DataFrame:
    """Generate optimal design using coordinate exchange algorithm.

    Parameters
    ----------
    max_iterations : int, optional
        Maximum number of exchange iterations. Defaults to ``1000``.
    random_start : bool, optional
        Whether to use random starting design. Defaults to ``True``.
    n_random_starts : int, optional
        Number of random starts to try. Defaults to ``5``.
    improvement_threshold : float, optional
        Minimum improvement in the criterion required to continue
        iterations. Defaults to ``1e-6``.

    Returns
    -------
    pd.DataFrame
        Optimal design matrix.
    """
    if not self.validate_design():
        raise ValueError("Invalid design configuration")

    if self.candidate_set is None:
        self.generate_candidate_set()
    # Precompute candidate model matrix for faster evaluation
    self.candidate_model_matrix = self._build_model_matrix(self.candidate_set)

    best_design = None
    best_criterion_value = float("-inf")

    # Try multiple random starts
    for _start in range(n_random_starts):
        design = self._coordinate_exchange(
            max_iterations, random_start, improvement_threshold
        )
        criterion_value = self._calculate_criterion(design)

        if self._is_better_criterion(criterion_value, best_criterion_value):
            best_design = design.copy()
            best_criterion_value = criterion_value

    self.design_matrix = best_design
    return self.design_matrix

validate_design

validate_design() -> bool

Validate optimal design parameters.

Source code in src/industrialstats/designs/optimal.py

def validate_design(self) -> bool:
    """Validate optimal design parameters."""
    if len(self.factors) == 0:
        return False

    if self.n_runs < len(self.model_terms):
        return False

    return self.criterion in ["D", "A", "G", "I"]

design_efficiency

design_efficiency(reference_design: DataFrame | None = None) -> dict[str, float]

Calculate design efficiency metrics.

Parameters:

Name	Type	Description	Default
reference_design	DataFrame	Reference design for comparison. Defaults to an orthogonal design.	None

Returns:

Type	Description
dict[str, float]	Efficiency metrics.

Source code in src/industrialstats/designs/optimal.py

def design_efficiency(
    self, reference_design: pd.DataFrame | None = None
) -> dict[str, float]:
    """Calculate design efficiency metrics.

    Parameters
    ----------
    reference_design : pd.DataFrame, optional
        Reference design for comparison. Defaults to an orthogonal design.

    Returns
    -------
    dict[str, float]
        Efficiency metrics.
    """
    if self.design_matrix is None:
        raise ValueError("Design not generated yet")

    # Calculate current design criterion
    current_X = self._build_model_matrix(self.design_matrix)
    current_XtX = current_X.T @ current_X

    efficiencies = {}

    try:
        if self.criterion == "D":
            current_det = np.linalg.det(current_XtX)

            # D-efficiency relative to orthogonal design
            p = len(self.model_terms)
            max_det = (self.n_runs / p) ** p  # Theoretical maximum
            efficiencies["D_efficiency"] = (current_det / max_det) ** (1 / p)

        elif self.criterion == "A":
            current_trace = np.trace(np.linalg.inv(current_XtX))

            # A-efficiency (simplified)
            p = len(self.model_terms)
            min_trace = p / self.n_runs  # Theoretical minimum
            efficiencies["A_efficiency"] = min_trace / current_trace

        # Relative efficiency compared to reference design
        if reference_design is not None:
            ref_X = self._build_model_matrix(reference_design)
            ref_XtX = ref_X.T @ ref_X

            if self.criterion == "D":
                ref_det = np.linalg.det(ref_XtX)
                efficiencies["Relative_D_efficiency"] = (current_det / ref_det) ** (
                    1 / p
                )

            elif self.criterion == "A":
                ref_trace = np.trace(np.linalg.inv(ref_XtX))
                efficiencies["Relative_A_efficiency"] = ref_trace / current_trace

    except np.linalg.LinAlgError:
        efficiencies["Error"] = "Singular information matrix"

    return efficiencies

prediction_variance_map

prediction_variance_map(factor1: str, factor2: str, grid_size: int = 20) -> tuple[ndarray, ndarray, ndarray]

Generate prediction variance map for two factors.

Parameters:

Name	Type	Description	Default
factor1	str	First factor name.	required
factor2	str	Second factor name.	required
grid_size	int	Grid resolution. Defaults to 20.	20

Returns:

Type	Description
tuple[ndarray, ndarray, ndarray]	X, Y, Z arrays for contour plotting.

Source code in src/industrialstats/designs/optimal.py

def prediction_variance_map(
    self, factor1: str, factor2: str, grid_size: int = 20
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Generate prediction variance map for two factors.

    Parameters
    ----------
    factor1 : str
        First factor name.
    factor2 : str
        Second factor name.
    grid_size : int, optional
        Grid resolution. Defaults to ``20``.

    Returns
    -------
    tuple[np.ndarray, np.ndarray, np.ndarray]
        X, Y, Z arrays for contour plotting.
    """
    if self.design_matrix is None:
        raise ValueError("Design not generated yet")

    # Get factor objects
    factor1_obj = next(f for f in self.factors if f.name == factor1)
    factor2_obj = next(f for f in self.factors if f.name == factor2)

    # Create grid
    if factor1_obj.factor_type == "continuous":
        x_range = np.linspace(
            min(factor1_obj.levels), max(factor1_obj.levels), grid_size
        )
    else:
        x_range = factor1_obj.levels

    if factor2_obj.factor_type == "continuous":
        y_range = np.linspace(
            min(factor2_obj.levels), max(factor2_obj.levels), grid_size
        )
    else:
        y_range = factor2_obj.levels

    X, Y = np.meshgrid(x_range, y_range)
    Z = np.zeros_like(X)

    # Calculate information matrix
    design_X = self._build_model_matrix(self.design_matrix)
    XtX_inv = np.linalg.inv(design_X.T @ design_X)

    # Calculate prediction variance at each grid point
    for i in range(len(y_range)):
        for j in range(len(x_range)):
            # Create point with other factors at center
            point_data = {}
            for factor in self.factors:
                if factor.name == factor1:
                    point_data[factor.name] = X[i, j]
                elif factor.name == factor2:
                    point_data[factor.name] = Y[i, j]
                else:
                    # Set other factors to center
                    if factor.factor_type == "continuous":
                        center = (min(factor.levels) + max(factor.levels)) / 2
                        point_data[factor.name] = center
                    else:
                        point_data[factor.name] = factor.levels[0]  # First level

            point = pd.Series(point_data)
            x_point = self._build_point_vector(point)
            variance = x_point.T @ XtX_inv @ x_point
            Z[i, j] = variance

    return X, Y, Z

augment_design

augment_design(additional_runs: int, current_data: DataFrame | None = None) -> DataFrame

Augment existing design with additional runs.

Parameters:

Name	Type	Description	Default
additional_runs	int	Number of additional runs to add.	required
current_data	DataFrame	Current experimental data. If None, uses the generated design.	None

Returns:

Type	Description
DataFrame	Augmented design.

Source code in src/industrialstats/designs/optimal.py

def augment_design(
    self, additional_runs: int, current_data: pd.DataFrame | None = None
) -> pd.DataFrame:
    """Augment existing design with additional runs.

    Parameters
    ----------
    additional_runs : int
        Number of additional runs to add.
    current_data : pd.DataFrame, optional
        Current experimental data. If ``None``, uses the generated design.

    Returns
    -------
    pd.DataFrame
        Augmented design.
    """
    if current_data is None:
        if self.design_matrix is None:
            raise ValueError("No current design available")
        current_data = self.design_matrix.copy()

    if self.candidate_set is None:
        self.generate_candidate_set()

    # Start with current design
    augmented_design = current_data.copy()
    current_n_runs = len(current_data)

    # Add runs one by one using exchange algorithm
    for new_run in range(additional_runs):
        best_addition = None
        best_criterion = (
            float("-inf") if self.criterion in ["D", "A"] else float("inf")
        )

        # Try each candidate point
        for _, candidate in self.candidate_set.iterrows():
            # Create trial design with new point
            trial_design = augmented_design.copy()
            new_point = {"RunID": current_n_runs + new_run + 1}
            for factor in self.factors:
                new_point[factor.name] = candidate[factor.name]

            trial_design = pd.concat(
                [trial_design, pd.DataFrame([new_point])], ignore_index=True
            )

            # Calculate criterion
            criterion_value = self._calculate_criterion(trial_design)

            if self._is_better_criterion(criterion_value, best_criterion):
                best_addition = new_point
                best_criterion = criterion_value

        # Add best point
        if best_addition is not None:
            augmented_design = pd.concat(
                [augmented_design, pd.DataFrame([best_addition])], ignore_index=True
            )

    return augmented_design

design_diagnostics

design_diagnostics() -> dict[str, Any]

Calculate design diagnostics and properties.

Returns:

Type	Description
dict[str, Any]	Diagnostic metrics for the current design.

Raises:

Type	Description
ValueError	If the design has not been generated.

Source code in src/industrialstats/designs/optimal.py

def design_diagnostics(self) -> dict[str, Any]:
    """Calculate design diagnostics and properties.

    Returns
    -------
    dict[str, Any]
        Diagnostic metrics for the current design.

    Raises
    ------
    ValueError
        If the design has not been generated.
    """
    if self.design_matrix is None:
        raise ValueError("Design not generated yet")

    diagnostics = {}

    # Basic properties
    diagnostics["n_runs"] = len(self.design_matrix)
    diagnostics["n_factors"] = len(self.factors)
    diagnostics["n_model_terms"] = len(self.model_terms)
    diagnostics["criterion"] = self.criterion

    # Model matrix properties
    X = self._build_model_matrix(self.design_matrix)
    XtX = X.T @ X

    try:
        # Condition number
        eigenvalues = np.linalg.eigvals(XtX)
        condition_number = np.max(eigenvalues) / np.min(eigenvalues)
        diagnostics["condition_number"] = condition_number

        # Determinant
        diagnostics["determinant"] = np.linalg.det(XtX)

        # Trace
        diagnostics["trace"] = np.trace(XtX)

        # Minimum eigenvalue
        diagnostics["min_eigenvalue"] = np.min(eigenvalues)

        # Design criterion value
        diagnostics["criterion_value"] = self._calculate_criterion(
            self.design_matrix
        )

        # Correlation matrix
        correlation_matrix = np.corrcoef(X.T)
        diagnostics["max_correlation"] = np.max(
            np.abs(correlation_matrix - np.eye(len(self.model_terms)))
        )

    except np.linalg.LinAlgError:
        diagnostics["error"] = "Singular information matrix"

    # Exchange algorithm convergence
    if self.exchange_history:
        diagnostics["exchange_iterations"] = len(self.exchange_history)
        diagnostics["final_improvement"] = self.exchange_history[-1]["improved"]

        # Convergence plot data
        criterion_values = [h["criterion_value"] for h in self.exchange_history]
        diagnostics["convergence_history"] = criterion_values

    return diagnostics

CustomOptimalDesign

CustomOptimalDesign(factors: list[Factor], n_runs: int, criterion_function: Callable[[ndarray], float], criterion_name: str = 'Custom')

Bases: OptimalDesign

Custom optimal design with user-defined criterion function.

Initialize custom optimal design.

Parameters:

Name	Type	Description	Default
factors	list[Factor]	Experimental factors.	required
n_runs	int	Number of experimental runs.	required
criterion_function	Callable[[ndarray], float]	Function that takes the model matrix X and returns a criterion value.	required
criterion_name	str	Name for the custom criterion. Defaults to "Custom".	'Custom'

Source code in src/industrialstats/designs/optimal.py

def __init__(
    self,
    factors: list[Factor],
    n_runs: int,
    criterion_function: Callable[[np.ndarray], float],
    criterion_name: str = "Custom",
) -> None:
    """Initialize custom optimal design.

    Parameters
    ----------
    factors : list[Factor]
        Experimental factors.
    n_runs : int
        Number of experimental runs.
    criterion_function : Callable[[np.ndarray], float]
        Function that takes the model matrix ``X`` and returns a criterion value.
    criterion_name : str, optional
        Name for the custom criterion. Defaults to ``"Custom"``.
    """
    super().__init__(factors, n_runs, criterion="D")  # Dummy criterion
    self.name = f"{criterion_name}-Optimal Design"
    self.criterion_function = criterion_function
    self.criterion_name = criterion_name

Advanced designs

industrialstats.designs.advanced

Advanced experimental designs.

SplitPlotDesign

SplitPlotDesign(whole_plot_factors: list[Factor], sub_plot_factors: list[Factor], replicates: int = 1, randomize: bool = True, seed: int | None = None)

Bases: ExperimentalDesign

Split-plot design with explicit whole-plot experimental units.

The design handles hard-to-change whole-plot factors and easier-to-change sub-plot factors. Each replicate of a whole-plot treatment combination is a distinct whole-plot experimental unit. Randomization is restricted so whole plots are randomized as units and sub-plots are randomized only within their parent whole plot.

Parameters:

Name	Type	Description	Default
whole_plot_factors	list[Factor]	Factors applied to whole plots (hard-to-change factors).	required
sub_plot_factors	list[Factor]	Factors applied within whole plots (easy-to-change factors).	required
replicates	int	Number of independent whole-plot replicates for each whole-plot treatment combination. Defaults to 1.	1
randomize	bool	Whether to randomize whole-plot order and sub-plot order within each whole plot. Defaults to True.	True
seed	int	Random seed for reproducible restricted randomization.	None

Examples:

Generate a split-plot design with one whole-plot factor and one sub-plot factor::

>>> from industrialstats.designs.base import Factor
>>> from industrialstats.designs.advanced import SplitPlotDesign
>>> wp = [Factor("Oven", [1, 2])]
>>> sp = [Factor("Temperature", [150, 200, 250])]
>>> design = SplitPlotDesign(wp, sp, seed=123)
>>> design.generate_design().head()
   RunOrder  StdOrder  Replicate  WholePlot  SubPlot  Oven  Temperature
0         1         4          1          2        1     2          150
1         2         6          1          2        3     2          250
2         3         5          1          2        2     2          200
3         4         1          1          1        1     1          150
4         5         3          1          1        3     1          250

Source code in src/industrialstats/designs/advanced.py

def __init__(
    self,
    whole_plot_factors: list[Factor],
    sub_plot_factors: list[Factor],
    replicates: int = 1,
    randomize: bool = True,
    seed: int | None = None,
) -> None:
    super().__init__("Split-Plot Design")
    if not whole_plot_factors:
        raise ValueError("At least one whole-plot factor is required")
    if not sub_plot_factors:
        raise ValueError("At least one sub-plot factor is required")
    if isinstance(replicates, bool) or not isinstance(replicates, int):
        raise ValueError("replicates must be an integer >= 1")
    if replicates < 1:
        raise ValueError("replicates must be >= 1")

    self.whole_plot_factors = whole_plot_factors
    self.sub_plot_factors = sub_plot_factors
    self.replicates = replicates
    self.randomize_flag = randomize
    self.seed = seed
    self.factors = whole_plot_factors + sub_plot_factors

generate_design

generate_design() -> DataFrame

Generate the split-plot design matrix.

Returns:

Type	Description
DataFrame	Generated design matrix with explicit Replicate, WholePlot, and SubPlot identifiers.

Source code in src/industrialstats/designs/advanced.py

def generate_design(self) -> pd.DataFrame:
    """Generate the split-plot design matrix.

    Returns
    -------
    pandas.DataFrame
        Generated design matrix with explicit ``Replicate``, ``WholePlot``,
        and ``SubPlot`` identifiers.
    """
    if not self.validate_design():
        raise ValueError("Invalid design configuration")

    design_rows: list[dict[str, Any]] = []
    run_id = 1
    whole_plot_id = 1
    wp_combinations = list(
        product(*(factor.levels for factor in self.whole_plot_factors))
    )
    sp_combinations = list(
        product(*(factor.levels for factor in self.sub_plot_factors))
    )

    for replicate in range(1, self.replicates + 1):
        for wp_combo in wp_combinations:
            for subplot_id, sp_combo in enumerate(sp_combinations, start=1):
                row: dict[str, Any] = {
                    "StdOrder": run_id,
                    "Replicate": replicate,
                    "WholePlot": whole_plot_id,
                    "SubPlot": subplot_id,
                }
                for index, factor in enumerate(self.whole_plot_factors):
                    row[factor.name] = wp_combo[index]
                for index, factor in enumerate(self.sub_plot_factors):
                    row[factor.name] = sp_combo[index]
                design_rows.append(row)
                run_id += 1
            whole_plot_id += 1

    self.design_matrix = pd.DataFrame(design_rows)

    if self.randomize_flag:
        self._randomize_restricted()

    return self.design_matrix

n_whole_plots

n_whole_plots() -> int

Return the number of independent whole-plot experimental units.

Source code in src/industrialstats/designs/advanced.py

def n_whole_plots(self) -> int:
    """Return the number of independent whole-plot experimental units."""
    return self.replicates * int(
        np.prod([len(factor.levels) for factor in self.whole_plot_factors])
    )

n_runs

n_runs() -> int

Return the total number of sub-plot experimental runs.

Source code in src/industrialstats/designs/advanced.py

def n_runs(self) -> int:
    """Return the total number of sub-plot experimental runs."""
    n_subplots_per_whole_plot = int(
        np.prod([len(factor.levels) for factor in self.sub_plot_factors])
    )
    return self.n_whole_plots() * n_subplots_per_whole_plot

validate_design

validate_design() -> bool

Validate split-plot design parameters.

Source code in src/industrialstats/designs/advanced.py

def validate_design(self) -> bool:
    """Validate split-plot design parameters."""
    return (
        bool(self.whole_plot_factors)
        and bool(self.sub_plot_factors)
        and isinstance(self.replicates, int)
        and not isinstance(self.replicates, bool)
        and self.replicates >= 1
        and all(factor.levels for factor in self.whole_plot_factors)
        and all(factor.levels for factor in self.sub_plot_factors)
    )

MixtureDesign

MixtureDesign(factors: list[Factor], order: int = 2, constraints: list[Callable[[ndarray], bool]] | None = None, randomize: bool = False, seed: int | None = None)

Bases: ExperimentalDesign

Simplex-lattice mixture design.

Parameters:

Name	Type	Description	Default
factors	list[Factor]	Mixture components. Must contain at least three continuous factors.	required
order	int	Simplex-lattice degree :math:m. Defaults to 2.	2
constraints	list[Callable[[ndarray], bool]]	Constraint functions applied to candidate mixtures. Each function receives an array of component proportions and returns True if the point satisfies the constraint. Defaults to None.	None
randomize	bool	Whether to randomize run order. Defaults to False.	False
seed	int	Random seed for reproducible randomization.	None

Examples:

Generate a simplex-lattice design for a three-component mixture::

>>> from industrialstats.designs.base import Factor
>>> from industrialstats.designs.advanced import MixtureDesign
>>> comps = [
...     Factor("A", [], "continuous"),
...     Factor("B", [], "continuous"),
...     Factor("C", [], "continuous"),
... ]
>>> design = MixtureDesign(comps, order=2)
>>> design.generate_design()
     A    B    C
0  1.0  0.0  0.0
1  0.0  1.0  0.0
2  0.0  0.0  1.0
3  0.5  0.5  0.0
4  0.5  0.0  0.5
5  0.0  0.5  0.5

References

.. [1] Cornell, J. A. (2011). Experiments with Mixtures.

Source code in src/industrialstats/designs/advanced.py

def __init__(
    self,
    factors: list[Factor],
    order: int = 2,
    constraints: list[Callable[[np.ndarray], bool]] | None = None,
    randomize: bool = False,
    seed: int | None = None,
) -> None:
    super().__init__("Mixture Design")
    if len(factors) < 3:
        raise ValueError("MixtureDesign requires at least three factors")
    self.factors = factors
    self.order = order
    self.constraints = constraints or []
    self.randomize_flag = randomize
    self.seed = seed

generate_design

generate_design() -> DataFrame

Generate the mixture design matrix.

Source code in src/industrialstats/designs/advanced.py

def generate_design(self) -> pd.DataFrame:
    """Generate the mixture design matrix."""
    if not self.validate_design():
        raise ValueError("Invalid mixture design configuration")

    points = self._generate_simplex_lattice()
    valid_points = []
    for pt in points:
        if all(constraint(pt) for constraint in self.constraints):
            if not np.isclose(pt.sum(), 1.0):
                raise ValueError("Mixture components must sum to 1")
            valid_points.append(pt)

    design = pd.DataFrame(valid_points, columns=[f.name for f in self.factors])

    if self.randomize_flag:
        rng = np.random.default_rng(self.seed)
        design = design.sample(
            frac=1,
            random_state=int(rng.integers(0, np.iinfo("int32").max)),
        ).reset_index(drop=True)
        design.insert(0, "RunOrder", range(1, len(design) + 1))
        self.randomized = True

    self.design_matrix = design
    return design

plot_simplex

plot_simplex(ax: Axes | None = None) -> Axes

Plot mixture design points for three components.

Parameters:

Name	Type	Description	Default
ax	Axes	Axes object to plot on. Created if None.	None

Returns:

Type	Description
Axes	Axes containing the simplex plot.

Raises:

Type	Description
ValueError	If the design has not been generated or the number of factors is not three.

Source code in src/industrialstats/designs/advanced.py

def plot_simplex(self, ax: Axes | None = None) -> Axes:
    """Plot mixture design points for three components.

    Parameters
    ----------
    ax : matplotlib.axes.Axes, optional
        Axes object to plot on. Created if ``None``.

    Returns
    -------
    matplotlib.axes.Axes
        Axes containing the simplex plot.

    Raises
    ------
    ValueError
        If the design has not been generated or the number of factors is
        not three.
    """
    if self.design_matrix is None:
        raise ValueError("Generate design before plotting")
    if len(self.factors) != 3:
        raise ValueError("Simplex plot currently supports three factors")

    import matplotlib.pyplot as plt

    if ax is None:
        _, ax = plt.subplots()

    data = self.design_matrix[[f.name for f in self.factors]].to_numpy()
    x = data[:, 1] + 0.5 * data[:, 2]
    y = (np.sqrt(3) / 2) * data[:, 2]
    ax.scatter(x, y)
    ax.set_xlabel(self.factors[1].name)
    ax.set_ylabel(self.factors[2].name)
    ax.set_title("Mixture Simplex")
    ax.set_aspect("equal")
    return ax

validate_design

validate_design() -> bool

Validate mixture design parameters.

Source code in src/industrialstats/designs/advanced.py

def validate_design(self) -> bool:
    """Validate mixture design parameters."""
    return len(self.factors) >= 3 and self.order >= 1