Utilities

Supporting utilities for validation, data handling, and export.

Validation

industrialstats.utils.validation

Validation utilities for experimental designs.

DesignValidator

Comprehensive design validation.

validate_factors staticmethod

validate_factors(factors: list[Factor]) -> list[str]

Validate factor specifications and return warnings.

Parameters:

Name	Type	Description	Default
factors	list of Factor	Factors to validate.	required

Returns:

Type	Description
list of str	Validation warnings, if any.

Source code in src/industrialstats/utils/validation.py

@staticmethod
def validate_factors(factors: list[Factor]) -> list[str]:
    """Validate factor specifications and return warnings.

    Parameters
    ----------
    factors : list of Factor
        Factors to validate.

    Returns
    -------
    list of str
        Validation warnings, if any.
    """
    warnings: list[str] = []
    factor_types = {f.factor_type for f in factors}
    if len(factor_types) > 1:
        warnings.append("Mixed factor types detected")

    for f in factors:
        if len(f.levels) < 2:
            warnings.append(f"Factor {f.name} has fewer than 2 levels")
        if len(set(f.levels)) != len(f.levels):
            warnings.append(f"Factor {f.name} has duplicate levels")
        if f.factor_type == "continuous" and not all(
            isinstance(level, (int, float)) for level in f.levels
        ):
            warnings.append(
                f"Factor {f.name} is continuous but has non-numeric levels"
            )
    return warnings

validate_design_matrix staticmethod

validate_design_matrix(design_matrix: DataFrame) -> dict[str, Any]

Validate a generated design matrix.

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Design matrix to inspect.	required

Returns:

Type	Description
dict	Validation summary including missing values, duplicates and single-level factors.

Source code in src/industrialstats/utils/validation.py

@staticmethod
def validate_design_matrix(design_matrix: pd.DataFrame) -> dict[str, Any]:
    """Validate a generated design matrix.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Design matrix to inspect.

    Returns
    -------
    dict
        Validation summary including missing values, duplicates and
        single-level factors.
    """
    result: dict[str, Any] = {
        "missing_values": design_matrix.isna().any().any(),
        "missing_counts": design_matrix.isna().sum().to_dict(),
        "duplicate_rows": design_matrix.duplicated().any(),
        "single_level_factors": [
            col
            for col in design_matrix.columns
            if design_matrix[col].nunique() <= 1
        ],
    }
    return result

check_confounding staticmethod

check_confounding(design_matrix: DataFrame) -> dict[str, Any]

Check for confounding patterns.

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Design matrix to analyze.	required

Returns:

Type	Description
dict	Dictionary containing high-correlation pairs, variance inflation factors (VIF), alias structures derived from the design matrix null space, and variance decomposition (:math:R^2) for each column.

References

.. [1] Box, G. E. P., Hunter, J. S., & Hunter, W. G. (2005). Statistics for Experimenters. .. [2] Montgomery, D. C. (2017). Design and Analysis of Experiments.

Source code in src/industrialstats/utils/validation.py

@staticmethod
def check_confounding(design_matrix: pd.DataFrame) -> dict[str, Any]:
    """Check for confounding patterns.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Design matrix to analyze.

    Returns
    -------
    dict
        Dictionary containing high-correlation pairs, variance inflation
        factors (VIF), alias structures derived from the design matrix null
        space, and variance decomposition (:math:`R^2`) for each column.

    References
    ----------
    .. [1] Box, G. E. P., Hunter, J. S., & Hunter, W. G. (2005). *Statistics
           for Experimenters*.
    .. [2] Montgomery, D. C. (2017). *Design and Analysis of Experiments*.
    """
    result: dict[str, Any] = {
        "high_correlation": {},
        "vif": {},
        "alias_structure": [],
        "variance_decomposition": {},
    }
    corr = design_matrix.corr(numeric_only=True).abs()
    for i, col in enumerate(corr.columns):
        for j in range(i + 1, len(corr.columns)):
            other = corr.columns[j]
            if corr.iloc[i, j] > 0.95:
                result["high_correlation"].setdefault(col, []).append(other)

    numeric = design_matrix.select_dtypes(include=[np.number])
    if numeric.shape[1] >= 2:
        X = numeric.values
        X_with_const = np.column_stack([np.ones(len(numeric)), X])
        # A perfectly confounded column has R^2 == 1, so the VIF
        # computation divides by zero and the design matrix is singular.
        # That is the defining case this validator exists to report, not an
        # anomaly, so let the arithmetic yield the mathematically correct
        # infinite inflation rather than surfacing numerical warnings that
        # the caller has already asked about by calling this function.
        with (
            warnings.catch_warnings(),
            np.errstate(divide="ignore", invalid="ignore"),
        ):
            # statsmodels signals rank deficiency and related model
            # problems through ModelWarning subclasses, and reports poor
            # conditioning as a plain UserWarning. Silence those two
            # specific channels rather than the whole UserWarning
            # category, so an unrelated diagnostic still reaches callers.
            warnings.filterwarnings("ignore", category=ModelWarning)
            warnings.filterwarnings(
                "ignore",
                message=".*poorly conditioned.*",
                category=UserWarning,
            )
            for i in range(1, X_with_const.shape[1]):
                result["vif"][numeric.columns[i - 1]] = variance_inflation_factor(
                    X_with_const, i
                )

        from scipy.linalg import null_space

        ns = null_space(X)
        for vec in ns.T:
            involved = [
                col
                for col, coeff in zip(numeric.columns, vec, strict=True)
                if abs(coeff) > 1e-10
            ]
            if involved:
                result["alias_structure"].append(involved)

        for i, col in enumerate(numeric.columns):
            y = X[:, i]
            X_other = np.delete(X, i, axis=1)
            if X_other.size == 0:
                continue
            beta, _, _, _ = np.linalg.lstsq(X_other, y, rcond=None)
            residuals = y - X_other @ beta
            ss_res = np.sum(residuals**2)
            ss_tot = np.sum((y - y.mean()) ** 2)
            r2 = 1 - ss_res / ss_tot if ss_tot > 0 else 1.0
            result["variance_decomposition"][col] = r2

    return result

estimate_power staticmethod

estimate_power(design_matrix: DataFrame, effect_size: float) -> float

Estimate design power for a given effect size.

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Design matrix.	required
effect_size	float	Expected effect size.	required

Returns:

Type	Description
float	Estimated statistical power.

Raises:

Type	Description
ValueError	If the design matrix is empty.

Source code in src/industrialstats/utils/validation.py

@staticmethod
def estimate_power(design_matrix: pd.DataFrame, effect_size: float) -> float:
    """Estimate design power for a given effect size.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Design matrix.
    effect_size : float
        Expected effect size.

    Returns
    -------
    float
        Estimated statistical power.

    Raises
    ------
    ValueError
        If the design matrix is empty.
    """
    from scipy.stats import f, ncf

    n = len(design_matrix)
    if n == 0:
        raise ValueError("Design matrix is empty")
    df_model = design_matrix.shape[1] - 1
    df_error = n - df_model - 1
    lambda_nc = effect_size**2 * n / 2
    f_crit = f.ppf(0.95, df_model, df_error)
    power = 1 - ncf.cdf(f_crit, df_model, df_error, lambda_nc)
    return power

Data generation

industrialstats.utils.data_generation

Utilities to simulate experimental data.

The :class:DataSimulator class centralizes routines for generating experimental and process-oriented data with rich noise structures. The implementation follows the mathematical guidelines in Montgomery [1] and Box & Jenkins [2] for factorial responses and stochastic process modelling, respectively.

DataSimulator

DataSimulator(seed: int | None = None)

Generate realistic experimental data.

Initialize the simulator.

Parameters:

Name	Type	Description	Default
seed	int	Random seed for reproducibility.	None

Source code in src/industrialstats/utils/data_generation.py

def __init__(self, seed: int | None = None) -> None:
    """Initialize the simulator.

    Parameters
    ----------
    seed : int, optional
        Random seed for reproducibility.
    """
    self.random_state = np.random.default_rng(seed)

simulate_factorial_response

simulate_factorial_response(design_matrix: DataFrame, main_effects: dict[str, float] | None = None, interactions: dict[tuple[str, str], float] | None = None, noise_level: float = 1.0, noise_dist: str = 'normal', noise_params: dict[str, float] | None = None, response_type: str = 'continuous', random_effects: dict[str, float] | None = None, corr: float = 0.0, heteroskedastic: Sequence[float] | Callable[[DataFrame], ndarray] | None = None, drift: float = 0.0, missing_rate: float = 0.0, missing_pattern: str = 'MCAR', measurement_error: dict[str, Any] | None = None) -> Series

Simulate response for a factorial design.

The deterministic part of the response follows the linear model

.. math:: y = Xeta + arepsilon,

where X is the encoded design matrix and arepsilon denotes the stochastic noise component. Interaction terms are formed by pairwise products of encoded factors. Optional random effects and AR(1) correlated noise may be superimposed on the deterministic structure.

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Design matrix.	required
main_effects	dict	Mapping of factor names to effect sizes. If None, all factors receive an effect size of 1.0.	None
interactions	dict	Mapping of (factor1, factor2) to interaction effect sizes.	None
noise_level	float	Scale of the random noise, by default 1.0.	1.0
noise_dist	('normal', 'laplace', 't', 'gamma', 'exponential')	Distribution for noise generation, by default 'normal'.	'normal'
noise_params	dict	Additional parameters for the selected distribution, e.g., degrees of freedom for 't'.	None
response_type	('continuous', 'binomial', 'poisson')	Type of response variable, by default 'continuous'.	'continuous'
random_effects	dict	Mapping of grouping column names to variance components for random intercepts.	None
corr	float	Correlation coefficient for AR(1) noise. A value of 0 implies independent errors.	0.0
heteroskedastic	sequence of float or callable	Observation-wise noise scales. Length must equal the number of design rows. Overrides noise_level when provided. If a callable is supplied, it receives the design matrix and must return a vector of scale factors.	None
drift	float	Linear drift coefficient applied in run order, by default 0.	0.0
missing_rate	float	Fraction of responses to set as missing. Must be in [0, 1].	0.0
missing_pattern	('MCAR', 'MAR', 'MNAR', 'block')	Missing-data mechanism. 'block' drops the last fraction of observations, 'MAR' and 'MNAR' implement missingness at random and not at random, respectively.	'MCAR'
measurement_error	dict	Parameters describing an additive measurement error model applied to the final response. Accepts {"scale": float, "distribution": str} following the same distribution names as noise_dist.	None

Returns:

Type	Description
Series	Simulated response values.

See Also

industrialstats.utils.validation.DesignValidator.check_confounding Assess correlation-based confounding in design matrices. industrialstats.analysis.power_analysis.factorial_power Power calculations for factorial designs.

Examples:

>>> import pandas as pd
>>> from industrialstats.utils.data_generation import DataSimulator
>>> dm = pd.DataFrame({"A": [1, -1, 1, -1], "B": [1, 1, -1, -1]})
>>> sim = DataSimulator(seed=1)
>>> sim.simulate_factorial_response(dm, main_effects={"A": 2, "B": 1}).round(2)
0    3.62
1    1.33
2    0.88
3   -3.53
Name: Response, dtype: float64

References

.. [1] Montgomery, D.C. (2017). Design and Analysis of Experiments. 9th ed. Wiley. .. [2] Box, G.E.P., Hunter, J.S., Hunter, W.G. (2005). Statistics for Experimenters, 2nd ed. Wiley. .. [3] Laird, N. M., & Ware, J. H. (1982). "Random-effects models for longitudinal data." Biometrics, 38(4), 963-974. .. [4] Carroll, R.J., Ruppert, D., Stefanski, L.A., & Crainiceanu, C.M. (2006). Measurement Error in Nonlinear Models, 2nd ed. Chapman & Hall/CRC.

Source code in src/industrialstats/utils/data_generation.py

def simulate_factorial_response(
    self,
    design_matrix: pd.DataFrame,
    main_effects: dict[str, float] | None = None,
    interactions: dict[tuple[str, str], float] | None = None,
    noise_level: float = 1.0,
    noise_dist: str = "normal",
    noise_params: dict[str, float] | None = None,
    response_type: str = "continuous",
    random_effects: dict[str, float] | None = None,
    corr: float = 0.0,
    heteroskedastic: Sequence[float]
    | Callable[[pd.DataFrame], np.ndarray]
    | None = None,
    drift: float = 0.0,
    missing_rate: float = 0.0,
    missing_pattern: str = "MCAR",
    measurement_error: dict[str, Any] | None = None,
) -> pd.Series:
    """Simulate response for a factorial design.

    The deterministic part of the response follows the linear model

    .. math:: y = X\beta + \varepsilon,

    where ``X`` is the encoded design matrix and ``\varepsilon`` denotes the
    stochastic noise component. Interaction terms are formed by pairwise
    products of encoded factors. Optional random effects and AR(1) correlated
    noise may be superimposed on the deterministic structure.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Design matrix.
    main_effects : dict, optional
        Mapping of factor names to effect sizes. If ``None``, all factors
        receive an effect size of 1.0.
    interactions : dict, optional
        Mapping of ``(factor1, factor2)`` to interaction effect sizes.
    noise_level : float, optional
        Scale of the random noise, by default 1.0.
    noise_dist : {'normal', 'laplace', 't', 'gamma', 'exponential'}, optional
        Distribution for noise generation, by default ``'normal'``.
    noise_params : dict, optional
        Additional parameters for the selected distribution, e.g., degrees of
        freedom for ``'t'``.
    response_type : {'continuous', 'binomial', 'poisson'}, optional
        Type of response variable, by default ``'continuous'``.
    random_effects : dict, optional
        Mapping of grouping column names to variance components for
        random intercepts.
    corr : float, optional
        Correlation coefficient for AR(1) noise. A value of ``0`` implies
        independent errors.
    heteroskedastic : sequence of float or callable, optional
        Observation-wise noise scales. Length must equal the number of
        design rows. Overrides ``noise_level`` when provided. If a callable is
        supplied, it receives the design matrix and must return a vector of
        scale factors.
    drift : float, optional
        Linear drift coefficient applied in run order, by default ``0``.
    missing_rate : float, optional
        Fraction of responses to set as missing. Must be in ``[0, 1]``.
    missing_pattern : {'MCAR', 'MAR', 'MNAR', 'block'}, optional
        Missing-data mechanism. ``'block'`` drops the last fraction of
        observations, ``'MAR'`` and ``'MNAR'`` implement missingness at random
        and not at random, respectively.
    measurement_error : dict, optional
        Parameters describing an additive measurement error model applied to
        the final response. Accepts ``{"scale": float, "distribution": str}``
        following the same distribution names as ``noise_dist``.

    Returns
    -------
    pandas.Series
        Simulated response values.

    See Also
    --------
    industrialstats.utils.validation.DesignValidator.check_confounding
        Assess correlation-based confounding in design matrices.
    industrialstats.analysis.power_analysis.factorial_power
        Power calculations for factorial designs.

    Examples
    --------
    >>> import pandas as pd
    >>> from industrialstats.utils.data_generation import DataSimulator
    >>> dm = pd.DataFrame({"A": [1, -1, 1, -1], "B": [1, 1, -1, -1]})
    >>> sim = DataSimulator(seed=1)
    >>> sim.simulate_factorial_response(dm, main_effects={"A": 2, "B": 1}).round(2)
    0    3.62
    1    1.33
    2    0.88
    3   -3.53
    Name: Response, dtype: float64

    References
    ----------
    .. [1] Montgomery, D.C. (2017). *Design and Analysis of Experiments*.
           9th ed. Wiley.
    .. [2] Box, G.E.P., Hunter, J.S., Hunter, W.G. (2005). *Statistics for
           Experimenters*, 2nd ed. Wiley.
    .. [3] Laird, N. M., & Ware, J. H. (1982). "Random-effects models for
           longitudinal data." *Biometrics*, 38(4), 963-974.
    .. [4] Carroll, R.J., Ruppert, D., Stefanski, L.A., & Crainiceanu, C.M.
           (2006). *Measurement Error in Nonlinear Models*, 2nd ed.
           Chapman & Hall/CRC.
    """
    if main_effects is None:
        main_effects = {
            c: 1.0
            for c in design_matrix.columns
            if c
            not in {"RunID", "Replicate", "DesignPoint", "StdOrder", "RunOrder"}
        }
    interactions = interactions or {}
    random_effects = random_effects or {}

    factor_cols = [
        c
        for c in design_matrix.columns
        if c not in {"RunID", "Replicate", "DesignPoint", "StdOrder", "RunOrder"}
    ]
    encoded_df = design_matrix[factor_cols].copy()
    for col in encoded_df.select_dtypes(exclude="number").columns:
        encoded_df[col] = encoded_df[col].astype("category").cat.codes

    encoded_array = encoded_df.to_numpy(dtype=float)
    coef_vector = np.array(
        [main_effects.get(col, 0.0) for col in encoded_df.columns], dtype=float
    )
    response = encoded_array @ coef_vector

    if interactions:
        cols = list(encoded_df.columns)
        term_matrix = [
            coef
            * encoded_array[:, cols.index(f1)]
            * encoded_array[:, cols.index(f2)]
            for (f1, f2), coef in interactions.items()
            if f1 in cols and f2 in cols
        ]
        if term_matrix:
            response += np.sum(np.stack(term_matrix, axis=0), axis=0)

    for col, var in random_effects.items():
        if col in design_matrix:
            groups = design_matrix[col]
            levels = groups.unique()
            re = self.random_state.normal(scale=np.sqrt(var), size=len(levels))
            mapping = dict(zip(levels, re, strict=True))
            response += groups.map(mapping).to_numpy()

    n = len(response)
    scales = self._resolve_scale(noise_level, heteroskedastic, design_matrix)
    if corr != 0.0:
        idx = np.arange(n)
        base = corr ** np.abs(np.subtract.outer(idx, idx))
        cov = np.outer(scales, scales) * base
        if noise_dist != "normal":
            raise ValueError(
                "Correlated noise currently supported only for normal distribution"
            )
        noise = self.random_state.multivariate_normal(np.zeros(n), cov)
    else:
        noise = self._draw_noise(noise_dist, n, scales, noise_params)

    response = response + noise

    if drift != 0.0:
        order = design_matrix.get("RunOrder", pd.Series(range(n))).to_numpy()
        response = response + drift * order

    if response_type == "continuous":
        final = response
    elif response_type == "binomial":
        p = 1 / (1 + np.exp(-response))
        final = self.random_state.binomial(1, p)
    elif response_type == "poisson":
        rate = np.exp(response)
        final = self.random_state.poisson(rate)
    else:
        raise ValueError("Unsupported response_type")

    final = pd.Series(final, name="Response", dtype=float)

    if measurement_error:
        err_scale = measurement_error.get("scale", noise_level)
        err_dist = measurement_error.get("distribution", "normal")
        meas_noise = self._draw_noise(
            err_dist,
            n,
            np.full(n, err_scale, dtype=float),
            measurement_error.get("params"),
        )
        final = final + meas_noise

    if missing_rate > 0:
        final = self._apply_missingness(
            final, missing_pattern, missing_rate, design_matrix
        )

    return final

simulate_correlated_responses

simulate_correlated_responses(design_matrix: DataFrame, main_effects_list: list[dict[str, float]], cov: ndarray, **kwargs) -> DataFrame

Simulate multiple correlated responses.

Each response uses simulate_factorial_response for its deterministic component. Correlated noise is then added using a multivariate normal distribution with covariance cov.

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Design matrix.	required
main_effects_list	list of dict	Main-effect specifications for each response.	required
cov	ndarray	Covariance matrix defining correlations between responses.	required
**kwargs		Additional arguments forwarded to :meth:simulate_factorial_response.	{}

Returns:

Type	Description
DataFrame	Simulated responses with one column per response.

Source code in src/industrialstats/utils/data_generation.py

def simulate_correlated_responses(
    self,
    design_matrix: pd.DataFrame,
    main_effects_list: list[dict[str, float]],
    cov: np.ndarray,
    **kwargs,
) -> pd.DataFrame:
    """Simulate multiple correlated responses.

    Each response uses ``simulate_factorial_response`` for its deterministic
    component. Correlated noise is then added using a multivariate normal
    distribution with covariance ``cov``.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Design matrix.
    main_effects_list : list of dict
        Main-effect specifications for each response.
    cov : numpy.ndarray
        Covariance matrix defining correlations between responses.
    **kwargs
        Additional arguments forwarded to
        :meth:`simulate_factorial_response`.

    Returns
    -------
    pandas.DataFrame
        Simulated responses with one column per response.
    """
    means = [
        self.simulate_factorial_response(
            design_matrix, main_effects=effects, noise_level=0.0, **kwargs
        ).to_numpy()
        for effects in main_effects_list
    ]
    mean_mat = np.column_stack(means)
    noise = self.random_state.multivariate_normal(
        np.zeros(len(main_effects_list)), cov, size=len(design_matrix)
    )
    result = mean_mat + noise
    cols = [f"Y{i + 1}" for i in range(len(main_effects_list))]
    return pd.DataFrame(result, columns=cols)

simulate_process_data

simulate_process_data(n_periods: int, model: Callable[[DataFrame], ndarray], covariates: DataFrame | None = None, freq: str = 'D', noise_dist: str = 'normal', noise_params: dict[str, float] | None = None, trend: dict[str, Any] | None = None, seasonality: dict[str, float] | None = None, ar_params: Sequence[float] | None = None, ma_params: Sequence[float] | None = None, heteroskedastic: Iterable[float] | Callable[[DataFrame], ndarray] | None = None, outliers: dict[str, dict[str, float | int]] | None = None, missing: dict[str, float | str] | None = None, measurement_error: dict[str, Any] | None = None, return_components: bool = False) -> DataFrame

Simulate a univariate or multivariate process response.

The function evaluates a custom model over a time-indexed frame and enriches the result with trend, seasonality, heteroskedasticity, and ARMA-style autocorrelation. Noise can be drawn from Gaussian, Student :math:t, Gamma, or Exponential families. Outlier generation follows the guidelines of Hawkins [3] whereas missingness mechanisms mirror Little & Rubin [4].

Parameters:

Name	Type	Description	Default
n_periods	int	Number of time points to simulate.	required
model	callable	Callable mapping a covariate frame to deterministic response values.	required
covariates	DataFrame	External drivers aligned with n_periods. A time column is appended automatically.	None
freq	str	Frequency string passed to :func:pandas.date_range for the time index, by default 'D'.	'D'
noise_dist	('normal', 'laplace', 't', 'gamma', 'exponential')	Distribution used for the innovation process.	'normal'
noise_params	dict	Distribution-specific parameters.	None
trend	dict	Specification of deterministic trend. Supports {'type': 'linear', 'slope': float, 'intercept': float} or {'type': 'poly', 'coeffs': Sequence[float]}.	None
seasonality	dict	Describes periodic fluctuations with keys 'period' and 'amplitude'. A 'phase' entry introduces a phase shift.	None
ar_params	sequence of float	Autoregressive parameters :math:\phi.	None
ma_params	sequence of float	Moving-average parameters :math:\theta.	None
heteroskedastic	iterable or callable	Observation-specific noise scales.	None
outliers	dict	Mapping of outlier types ('random', 'systematic', 'leverage') to configuration dictionaries.	None
missing	dict	Missing-data mechanism configuration with keys 'rate' and 'mechanism'.	None
measurement_error	dict	Additive measurement error with keys 'scale' and 'distribution' (defaults to Gaussian).	None
return_components	bool	If True, the returned frame includes deterministic and stochastic components for diagnostics.	False

Returns:

Type	Description
DataFrame	Simulated process data containing the response and metadata columns.

References

.. [1] Montgomery, D.C. (2017). Design and Analysis of Experiments, 9th ed. Wiley. .. [2] Box, G.E.P., Jenkins, G.M., Reinsel, G.C., & Ljung, G.M. (2015). Time Series Analysis: Forecasting and Control, 5th ed. Wiley. .. [3] Hawkins, D.M. (1980). Identification of Outliers. Chapman and Hall. .. [4] Little, R.J.A., & Rubin, D.B. (2002). Statistical Analysis with Missing Data, 2nd ed. Wiley.

Source code in src/industrialstats/utils/data_generation.py

def simulate_process_data(
    self,
    n_periods: int,
    model: Callable[[pd.DataFrame], np.ndarray],
    covariates: pd.DataFrame | None = None,
    freq: str = "D",
    noise_dist: str = "normal",
    noise_params: dict[str, float] | None = None,
    trend: dict[str, Any] | None = None,
    seasonality: dict[str, float] | None = None,
    ar_params: Sequence[float] | None = None,
    ma_params: Sequence[float] | None = None,
    heteroskedastic: Iterable[float]
    | Callable[[pd.DataFrame], np.ndarray]
    | None = None,
    outliers: dict[str, dict[str, float | int]] | None = None,
    missing: dict[str, float | str] | None = None,
    measurement_error: dict[str, Any] | None = None,
    return_components: bool = False,
) -> pd.DataFrame:
    r"""Simulate a univariate or multivariate process response.

    The function evaluates a custom ``model`` over a time-indexed frame and
    enriches the result with trend, seasonality, heteroskedasticity, and
    ARMA-style autocorrelation. Noise can be drawn from Gaussian, Student
    :math:`t`, Gamma, or Exponential families. Outlier generation follows the
    guidelines of Hawkins [3]_ whereas missingness mechanisms mirror Little &
    Rubin [4]_.

    Parameters
    ----------
    n_periods : int
        Number of time points to simulate.
    model : callable
        Callable mapping a covariate frame to deterministic response values.
    covariates : pandas.DataFrame, optional
        External drivers aligned with ``n_periods``. A ``time`` column is
        appended automatically.
    freq : str, optional
        Frequency string passed to :func:`pandas.date_range` for the time
        index, by default ``'D'``.
    noise_dist : {'normal', 'laplace', 't', 'gamma', 'exponential'}, optional
        Distribution used for the innovation process.
    noise_params : dict, optional
        Distribution-specific parameters.
    trend : dict, optional
        Specification of deterministic trend. Supports ``{'type': 'linear',
        'slope': float, 'intercept': float}`` or ``{'type': 'poly',
        'coeffs': Sequence[float]}``.
    seasonality : dict, optional
        Describes periodic fluctuations with keys ``'period'`` and
        ``'amplitude'``. A ``'phase'`` entry introduces a phase shift.
    ar_params : sequence of float, optional
        Autoregressive parameters :math:`\phi`.
    ma_params : sequence of float, optional
        Moving-average parameters :math:`\theta`.
    heteroskedastic : iterable or callable, optional
        Observation-specific noise scales.
    outliers : dict, optional
        Mapping of outlier types (``'random'``, ``'systematic'``,
        ``'leverage'``) to configuration dictionaries.
    missing : dict, optional
        Missing-data mechanism configuration with keys ``'rate'`` and
        ``'mechanism'``.
    measurement_error : dict, optional
        Additive measurement error with keys ``'scale'`` and
        ``'distribution'`` (defaults to Gaussian).
    return_components : bool, optional
        If ``True``, the returned frame includes deterministic and stochastic
        components for diagnostics.

    Returns
    -------
    pandas.DataFrame
        Simulated process data containing the response and metadata columns.

    References
    ----------
    .. [1] Montgomery, D.C. (2017). *Design and Analysis of Experiments*,
           9th ed. Wiley.
    .. [2] Box, G.E.P., Jenkins, G.M., Reinsel, G.C., & Ljung, G.M. (2015).
           *Time Series Analysis: Forecasting and Control*, 5th ed. Wiley.
    .. [3] Hawkins, D.M. (1980). *Identification of Outliers*. Chapman and
           Hall.
    .. [4] Little, R.J.A., & Rubin, D.B. (2002). *Statistical Analysis with
           Missing Data*, 2nd ed. Wiley.
    """

    if n_periods <= 0:
        raise ValueError("n_periods must be positive")

    time_index = pd.date_range("2000-01-01", periods=n_periods, freq=freq)
    base_df = covariates.copy() if covariates is not None else pd.DataFrame()
    base_df = base_df.reset_index(drop=True)
    if not base_df.empty and len(base_df) != n_periods:
        raise ValueError("covariates must have n_periods rows")
    base_df["time"] = time_index
    base_df["t"] = np.arange(n_periods, dtype=float)

    deterministic = np.asarray(model(base_df), dtype=float)
    if deterministic.shape not in {(n_periods,), (n_periods, 1)}:
        raise ValueError("model must return a vector of length n_periods")
    deterministic = deterministic.reshape(n_periods)

    if trend:
        kind = trend.get("type", "linear").lower()
        if kind == "linear":
            slope = float(trend.get("slope", 0.0))
            intercept = float(trend.get("intercept", 0.0))
            deterministic = (
                deterministic + intercept + slope * base_df["t"].to_numpy()
            )
        elif kind == "poly":
            coeffs = trend.get("coeffs", (0.0,))
            poly = np.poly1d(coeffs)
            deterministic = deterministic + poly(base_df["t"].to_numpy())
        else:
            raise ValueError("Unsupported trend specification")

    if seasonality:
        period = float(seasonality.get("period", 12.0))
        amplitude = float(seasonality.get("amplitude", 1.0))
        phase = float(seasonality.get("phase", 0.0))
        deterministic = deterministic + amplitude * np.sin(
            2 * np.pi * (base_df["t"].to_numpy() + phase) / max(period, 1e-6)
        )

    scale_vec = self._resolve_scale(1.0, heteroskedastic, base_df)
    base_noise = self._draw_noise(noise_dist, n_periods, scale_vec, noise_params)
    stochastic = self._arma_filter(base_noise, ar_params, ma_params)
    response = deterministic + stochastic
    response = self._inject_outliers(response, outliers, scale_vec)

    result = pd.DataFrame({"time": time_index, "response": response})
    result["deterministic"] = deterministic
    result["stochastic"] = stochastic

    if measurement_error:
        err_scale = float(measurement_error.get("scale", 0.1))
        err_dist = str(measurement_error.get("distribution", "normal"))
        meas = self._draw_noise(
            err_dist,
            n_periods,
            np.full(n_periods, err_scale),
            measurement_error.get("params"),
        )
        result["response"] = result["response"] + meas
        result["measurement_error"] = meas

    if missing:
        mechanism = str(missing.get("mechanism", "MCAR"))
        rate = float(missing.get("rate", 0.0))
        result["response"] = self._apply_missingness(
            result["response"], mechanism, rate, base_df
        )

    if not return_components:
        cols_to_keep = ["time", "response"]
        if "measurement_error" in result:
            cols_to_keep.append("measurement_error")
        result = result[cols_to_keep]

    return result

simulate_multi_response

simulate_multi_response(design_matrix: DataFrame, response_models: Sequence[Callable[[DataFrame], ndarray]], covariance: ndarray, response_types: Sequence[str] | None = None, noise_scales: Sequence[float] | None = None, noise_dist: str = 'normal', measurement_error: Sequence[dict[str, Any] | None] | None = None) -> DataFrame

Simulate correlated multi-response experimental outcomes.

Each response is constructed from a deterministic model augmented with a correlated latent noise term drawn from covariance. Response types may be continuous, categorical (binary logistic), or count (Poisson).

References

.. [1] Khuri, A.I., & Cornell, J.A. (1996). Response Surfaces: Design and Analyses. CRC Press. .. [2] Johnson, R.A., & Wichern, D.W. (2007). Applied Multivariate Statistical Analysis, 6th ed. Pearson.

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Design or feature matrix shared across responses.	required
response_models	sequence of callable	Deterministic response functions applied to the design matrix.	required
covariance	ndarray	Positive semi-definite covariance matrix governing latent noise.	required
response_types	sequence of {'continuous', 'categorical', 'count'}	Specifies the distribution of each response. Defaults to continuous.	None
noise_scales	sequence of float	Additional scale multipliers applied per response.	None
noise_dist	('normal', 't', 'gamma', 'exponential', 'laplace')	Distribution used to generate latent noise prior to correlating.	'normal'
measurement_error	sequence of dict	Optional measurement-error configuration for each response using the same schema as :meth:simulate_factorial_response.	None

Returns:

Type	Description
DataFrame	Multi-response dataset preserving the original design columns.

Source code in src/industrialstats/utils/data_generation.py

def simulate_multi_response(
    self,
    design_matrix: pd.DataFrame,
    response_models: Sequence[Callable[[pd.DataFrame], np.ndarray]],
    covariance: np.ndarray,
    response_types: Sequence[str] | None = None,
    noise_scales: Sequence[float] | None = None,
    noise_dist: str = "normal",
    measurement_error: Sequence[dict[str, Any] | None] | None = None,
) -> pd.DataFrame:
    """Simulate correlated multi-response experimental outcomes.

    Each response is constructed from a deterministic model augmented with a
    correlated latent noise term drawn from ``covariance``. Response types
    may be continuous, categorical (binary logistic), or count (Poisson).

    References
    ----------
    .. [1] Khuri, A.I., & Cornell, J.A. (1996). *Response Surfaces: Design
           and Analyses*. CRC Press.
    .. [2] Johnson, R.A., & Wichern, D.W. (2007). *Applied Multivariate
           Statistical Analysis*, 6th ed. Pearson.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Design or feature matrix shared across responses.
    response_models : sequence of callable
        Deterministic response functions applied to the design matrix.
    covariance : numpy.ndarray
        Positive semi-definite covariance matrix governing latent noise.
    response_types : sequence of {'continuous', 'categorical', 'count'}, optional
        Specifies the distribution of each response. Defaults to continuous.
    noise_scales : sequence of float, optional
        Additional scale multipliers applied per response.
    noise_dist : {'normal', 't', 'gamma', 'exponential', 'laplace'}, optional
        Distribution used to generate latent noise prior to correlating.
    measurement_error : sequence of dict, optional
        Optional measurement-error configuration for each response using the
        same schema as :meth:`simulate_factorial_response`.

    Returns
    -------
    pandas.DataFrame
        Multi-response dataset preserving the original design columns.
    """

    n = len(design_matrix)
    r = len(response_models)
    if covariance.shape != (r, r):
        raise ValueError(
            "covariance must be square with dimension equal to responses"
        )
    response_types = response_types or ["continuous"] * r
    if len(response_types) != r:
        raise ValueError("response_types length must match response_models")
    if noise_scales and len(noise_scales) != r:
        raise ValueError("noise_scales length must match response_models")

    deterministic_parts = []
    for model in response_models:
        deterministic_parts.append(np.asarray(model(design_matrix), dtype=float))
    deterministic = np.column_stack(deterministic_parts)

    noise_scales = noise_scales or [1.0] * r
    latent = self.random_state.multivariate_normal(np.zeros(r), covariance, size=n)
    if noise_dist != "normal":
        base = self._draw_noise(
            noise_dist,
            n * r,
            np.repeat(1.0, n * r),
            None,
        ).reshape(n, r)
        chol = np.linalg.cholesky(covariance + 1e-12 * np.eye(r))
        latent = base @ chol.T
    latent = latent * np.asarray(noise_scales, dtype=float)

    responses = np.empty_like(latent)
    columns = []
    meas_configs: list[dict[str, Any] | None] = (
        list(measurement_error) if measurement_error else [None] * r
    )
    for idx, (det_col, noise_col, kind, meas_cfg) in enumerate(
        zip(deterministic.T, latent.T, response_types, meas_configs, strict=True)
    ):
        column_name = f"Y{idx + 1}"
        columns.append(column_name)
        combined = det_col + noise_col
        if kind == "continuous":
            resp = combined
        elif kind == "categorical":
            prob = 1 / (1 + np.exp(-combined))
            resp = self.random_state.binomial(1, prob)
        elif kind == "count":
            rate = np.clip(np.exp(combined), 1e-9, None)
            resp = self.random_state.poisson(rate)
        else:
            raise ValueError("Unsupported response type")
        if meas_cfg:
            scale = float(meas_cfg.get("scale", 0.1))
            dist = str(meas_cfg.get("distribution", "normal"))
            resp = resp + self._draw_noise(
                dist, n, np.full(n, scale, dtype=float), meas_cfg.get("params")
            )
        responses[:, idx] = resp

    out_df = design_matrix.reset_index(drop=True).copy()
    for column, values in zip(columns, responses.T, strict=True):
        out_df[column] = values
    return out_df

validate_against_real_data

validate_against_real_data(simulated: DataFrame | Series, real_data: DataFrame | Series) -> dict[str, dict[str, float]]

Compare simulated data to real experimental measurements.

The function computes absolute differences in means and standard deviations for each variable, allowing users to gauge similarity between simulated and actual data sets.

Parameters:

Name	Type	Description	Default
simulated	Series or DataFrame	Simulated responses.	required
real_data	Series or DataFrame	Empirical measurements to compare against.	required

Returns:

Type	Description
dict of dict	Mapping each column name to {"mean_diff": float, "std_diff": float}.

Source code in src/industrialstats/utils/data_generation.py

def validate_against_real_data(
    self, simulated: pd.DataFrame | pd.Series, real_data: pd.DataFrame | pd.Series
) -> dict[str, dict[str, float]]:
    """Compare simulated data to real experimental measurements.

    The function computes absolute differences in means and standard
    deviations for each variable, allowing users to gauge similarity between
    simulated and actual data sets.

    Parameters
    ----------
    simulated : pandas.Series or pandas.DataFrame
        Simulated responses.
    real_data : pandas.Series or pandas.DataFrame
        Empirical measurements to compare against.

    Returns
    -------
    dict of dict
        Mapping each column name to ``{"mean_diff": float, "std_diff": float}``.
    """
    sim_df = pd.DataFrame(simulated)
    real_df = pd.DataFrame(real_data)[sim_df.columns]
    stats: dict[str, dict[str, float]] = {}
    for col in sim_df.columns:
        sim_col = sim_df[col].dropna()
        real_col = real_df[col].dropna()
        ks_stat = float(
            np.max(
                np.abs(
                    np.sort(sim_col.to_numpy())
                    - np.sort(real_col.to_numpy())[: len(sim_col)]
                )
            )
        )
        stats[col] = {
            "mean_diff": float(abs(sim_col.mean() - real_col.mean())),
            "std_diff": float(abs(sim_col.std() - real_col.std())),
            "ks_like": ks_stat,
        }
    return stats

Transforms

industrialstats.utils.transforms

Data transformation helpers.

center

center(df: DataFrame) -> DataFrame

Center numeric columns around zero.

Parameters:

Name	Type	Description	Default
df	DataFrame	Input DataFrame.	required

Returns:

Type	Description
DataFrame	Centered DataFrame.

Source code in src/industrialstats/utils/transforms.py

def center(df: pd.DataFrame) -> pd.DataFrame:
    """Center numeric columns around zero.

    Parameters
    ----------
    df : pandas.DataFrame
        Input DataFrame.

    Returns
    -------
    pandas.DataFrame
        Centered DataFrame.
    """
    centered = df.copy()
    for col in centered.select_dtypes(include=[np.number]).columns:
        centered[col] = centered[col] - centered[col].mean()
    return centered

standardize

standardize(df: DataFrame) -> DataFrame

Standardize numeric columns to unit variance.

Parameters:

Name	Type	Description	Default
df	DataFrame	Input DataFrame.	required

Returns:

Type	Description
DataFrame	Standardized DataFrame.

Source code in src/industrialstats/utils/transforms.py

def standardize(df: pd.DataFrame) -> pd.DataFrame:
    """Standardize numeric columns to unit variance.

    Parameters
    ----------
    df : pandas.DataFrame
        Input DataFrame.

    Returns
    -------
    pandas.DataFrame
        Standardized DataFrame.
    """
    standardized = df.copy()
    for col in standardized.select_dtypes(include=[np.number]).columns:
        std = standardized[col].std(ddof=0)
        if std != 0:
            standardized[col] = (standardized[col] - standardized[col].mean()) / std
    return standardized

log_transform

log_transform(df: DataFrame, columns: list[str]) -> DataFrame

Apply natural logarithm to specified columns.

Parameters:

Name	Type	Description	Default
df	DataFrame	Input DataFrame.	required
columns	list[str]	Columns to transform.	required

Returns:

Type	Description
DataFrame	DataFrame with transformed columns.

Source code in src/industrialstats/utils/transforms.py

def log_transform(df: pd.DataFrame, columns: list[str]) -> pd.DataFrame:
    """Apply natural logarithm to specified columns.

    Parameters
    ----------
    df : pandas.DataFrame
        Input DataFrame.
    columns : list[str]
        Columns to transform.

    Returns
    -------
    pandas.DataFrame
        DataFrame with transformed columns.
    """
    transformed = df.copy()
    for col in columns:
        transformed[col] = np.log(transformed[col])
    return transformed

Efficiency

industrialstats.utils.efficiency

Design efficiency metrics and visualization utilities.

d_efficiency

d_efficiency(design_matrix: DataFrame) -> float

Compute D-efficiency of a design.

D-efficiency is defined as (\det(X^T X)^{1/p}) / n where p is the number of parameters and n is the run count.

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Encoded design matrix including intercept.	required

Returns:

Type	Description
float	D-efficiency value.

References

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

Source code in src/industrialstats/utils/efficiency.py

def d_efficiency(design_matrix: pd.DataFrame) -> float:
    r"""Compute D-efficiency of a design.

    D-efficiency is defined as ``(\det(X^T X)^{1/p}) / n`` where ``p`` is the
    number of parameters and ``n`` is the run count.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Encoded design matrix including intercept.

    Returns
    -------
    float
        D-efficiency value.

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

    info = _information_matrix(design_matrix)
    n, p = design_matrix.shape
    det = np.linalg.det(info)
    return float(det ** (1 / p) / n)

a_efficiency

a_efficiency(design_matrix: DataFrame) -> float

Compute A-efficiency of a design.

A-efficiency is p / (\operatorname{trace}((X^T X)^{-1}) \cdot n).

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Encoded design matrix including intercept.	required

Returns:

Type	Description
float	A-efficiency value.

Source code in src/industrialstats/utils/efficiency.py

def a_efficiency(design_matrix: pd.DataFrame) -> float:
    r"""Compute A-efficiency of a design.

    A-efficiency is ``p / (\operatorname{trace}((X^T X)^{-1}) \cdot n)``.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Encoded design matrix including intercept.

    Returns
    -------
    float
        A-efficiency value.
    """

    info = _information_matrix(design_matrix)
    n, p = design_matrix.shape
    inv_trace = np.trace(np.linalg.inv(info))
    return float(p / (inv_trace * n))

g_efficiency

g_efficiency(design_matrix: DataFrame, candidate_points: DataFrame) -> float

Compute G-efficiency for a design.

G-efficiency is the reciprocal of the maximum scaled prediction variance over the candidate set.

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Design matrix used to fit the model.	required
candidate_points	DataFrame	Candidate matrix covering the region of interest.	required

Returns:

Type	Description
float	G-efficiency value.

Source code in src/industrialstats/utils/efficiency.py

def g_efficiency(
    design_matrix: pd.DataFrame,
    candidate_points: pd.DataFrame,
) -> float:
    """Compute G-efficiency for a design.

    G-efficiency is the reciprocal of the maximum scaled prediction variance
    over the candidate set.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Design matrix used to fit the model.
    candidate_points : pandas.DataFrame
        Candidate matrix covering the region of interest.

    Returns
    -------
    float
        G-efficiency value.
    """

    info_inv = np.linalg.inv(_information_matrix(design_matrix))
    Xc = np.asarray(candidate_points)
    pv = np.einsum("ij,jk,ik->i", Xc, info_inv, Xc)
    return float(1 / pv.max())

i_efficiency

i_efficiency(design_matrix: DataFrame, candidate_points: DataFrame) -> float

Compute I-efficiency for a design.

I-efficiency is the reciprocal of the average scaled prediction variance over the candidate set.

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Design matrix used to fit the model.	required
candidate_points	DataFrame	Candidate matrix covering the region of interest.	required

Returns:

Type	Description
float	I-efficiency value.

Source code in src/industrialstats/utils/efficiency.py

def i_efficiency(
    design_matrix: pd.DataFrame,
    candidate_points: pd.DataFrame,
) -> float:
    """Compute I-efficiency for a design.

    I-efficiency is the reciprocal of the average scaled prediction variance
    over the candidate set.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Design matrix used to fit the model.
    candidate_points : pandas.DataFrame
        Candidate matrix covering the region of interest.

    Returns
    -------
    float
        I-efficiency value.
    """

    info_inv = np.linalg.inv(_information_matrix(design_matrix))
    Xc = np.asarray(candidate_points)
    pv = np.einsum("ij,jk,ik->i", Xc, info_inv, Xc)
    return float(1 / pv.mean())

relative_efficiency

relative_efficiency(design_a: DataFrame, design_b: DataFrame, metric: str = 'D') -> float

Compare efficiencies of two designs.

Parameters:

Name	Type	Description	Default
design_a	DataFrame	Design matrices to compare.	required
design_b	DataFrame	Design matrices to compare.	required
metric	('D', 'A')	Efficiency measure for comparison, by default "D".	"D"

Returns:

Type	Description
float	Relative efficiency eff_a / eff_b.

Source code in src/industrialstats/utils/efficiency.py

def relative_efficiency(
    design_a: pd.DataFrame,
    design_b: pd.DataFrame,
    metric: str = "D",
) -> float:
    """Compare efficiencies of two designs.

    Parameters
    ----------
    design_a, design_b : pandas.DataFrame
        Design matrices to compare.
    metric : {"D", "A"}, optional
        Efficiency measure for comparison, by default ``"D"``.

    Returns
    -------
    float
        Relative efficiency ``eff_a / eff_b``.
    """

    metrics = {"D": d_efficiency, "A": a_efficiency}
    if metric not in metrics:
        raise ValueError("metric must be 'D' or 'A'")
    eff_a = metrics[metric](design_a)
    eff_b = metrics[metric](design_b)
    return float(eff_a / eff_b)

variance_inflation_factors

variance_inflation_factors(design_matrix: DataFrame) -> Series

Compute variance inflation factors (VIF) for regressors.

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Encoded design matrix including intercept.	required

Returns:

Type	Description
Series	VIF values indexed by column name.

Source code in src/industrialstats/utils/efficiency.py

def variance_inflation_factors(design_matrix: pd.DataFrame) -> pd.Series:
    """Compute variance inflation factors (VIF) for regressors.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Encoded design matrix including intercept.

    Returns
    -------
    pandas.Series
        VIF values indexed by column name.
    """

    X = np.asarray(design_matrix)
    info = X.T @ X
    info_inv = np.linalg.inv(info)
    vif_vals = np.diag(info_inv) * np.diag(info)
    return pd.Series(vif_vals, index=design_matrix.columns)

estimate_power

estimate_power(design_matrix: DataFrame, effect_contrast: ndarray, effect_size: float, sigma: float = 1.0, alpha: float = 0.05) -> float

Estimate power for detecting a specified contrast.

Parameters:

Name	Type	Description	Default
design_matrix	DataFrame	Encoded design matrix including intercept.	required
effect_contrast	ndarray	Contrast vector c specifying the linear combination of coefficients under test.	required
effect_size	float	Magnitude of the effect along c.	required
sigma	float	Residual standard deviation, by default 1.0.	1.0
alpha	float	Significance level, by default 0.05.	0.05

Returns:

Type	Description
float	Approximate statistical power.

Source code in src/industrialstats/utils/efficiency.py

def estimate_power(
    design_matrix: pd.DataFrame,
    effect_contrast: np.ndarray,
    effect_size: float,
    sigma: float = 1.0,
    alpha: float = 0.05,
) -> float:
    """Estimate power for detecting a specified contrast.

    Parameters
    ----------
    design_matrix : pandas.DataFrame
        Encoded design matrix including intercept.
    effect_contrast : numpy.ndarray
        Contrast vector ``c`` specifying the linear combination of coefficients
        under test.
    effect_size : float
        Magnitude of the effect along ``c``.
    sigma : float, optional
        Residual standard deviation, by default ``1.0``.
    alpha : float, optional
        Significance level, by default ``0.05``.

    Returns
    -------
    float
        Approximate statistical power.
    """

    X = np.asarray(design_matrix)
    n, p = X.shape
    info_inv = np.linalg.inv(X.T @ X)
    se = sigma * np.sqrt(effect_contrast @ info_inv @ effect_contrast)
    df = n - p
    tcrit = stats.t.ppf(1 - alpha / 2, df)
    ncp = effect_size / se
    power = 1 - stats.nct.cdf(tcrit, df, ncp) + stats.nct.cdf(-tcrit, df, ncp)
    return float(power)

plot_efficiencies

plot_efficiencies(efficiencies: dict[str, float]) -> Axes

Plot efficiency metrics for multiple designs.

Parameters:

Name	Type	Description	Default
efficiencies	dict	Mapping of design labels to efficiency values.	required

Returns:

Type	Description
Axes	Axes containing a bar chart of efficiencies.

Source code in src/industrialstats/utils/efficiency.py

def plot_efficiencies(efficiencies: dict[str, float]) -> plt.Axes:
    """Plot efficiency metrics for multiple designs.

    Parameters
    ----------
    efficiencies : dict
        Mapping of design labels to efficiency values.

    Returns
    -------
    matplotlib.axes.Axes
        Axes containing a bar chart of efficiencies.
    """

    labels = list(efficiencies)
    values = [efficiencies[k] for k in labels]
    _fig, ax = plt.subplots()
    ax.bar(labels, values, color="steelblue")
    ax.set_ylabel("Efficiency")
    ax.set_ylim(0, max(values) * 1.1)
    ax.set_title("Design Efficiency Comparison")
    return ax

Input and output

industrialstats.utils.io

Structured data-loading helpers for external tabular inputs.

load_csv

load_csv(path: str | Path, **kwargs: Any) -> DataFrame

Load a CSV file and expose operational failures through DataExcept.

Parameters:

Name	Type	Description	Default
path	str or Path	CSV source path.	required
**kwargs	Any	Additional arguments passed to :func:pandas.read_csv.	{}

Returns:

Type	Description
DataFrame	Loaded tabular data.

Raises:

Type	Description
DataLoadingError	If pandas or the filesystem cannot load the CSV source.

Source code in src/industrialstats/utils/io.py

def load_csv(path: str | Path, **kwargs: Any) -> pd.DataFrame:
    """Load a CSV file and expose operational failures through DataExcept.

    Parameters
    ----------
    path : str or pathlib.Path
        CSV source path.
    **kwargs
        Additional arguments passed to :func:`pandas.read_csv`.

    Returns
    -------
    pandas.DataFrame
        Loaded tabular data.

    Raises
    ------
    DataLoadingError
        If pandas or the filesystem cannot load the CSV source.
    """
    try:
        return pd.read_csv(path, **kwargs)
    except (OSError, ValueError, UnicodeError) as exc:
        raise DataLoadingError(str(path), exc) from exc

Export

industrialstats.utils.export

Data export utilities for industrialstats.

export_to_csv

export_to_csv(df: DataFrame, path: str | Path, include_index: bool = False, **kwargs: Any) -> None

Save a DataFrame to CSV.

Parameters:

Name	Type	Description	Default
df	DataFrame	DataFrame to export.	required
path	str or Path	Destination file path.	required
include_index	bool	Whether to include the DataFrame index. Defaults to False.	False
**kwargs	Any	Additional arguments passed to :func:pandas.DataFrame.to_csv.	{}

Raises:

Type	Description
FileWriteError	If pandas or the filesystem cannot write the destination.

Source code in src/industrialstats/utils/export.py

def export_to_csv(
    df: pd.DataFrame,
    path: str | Path,
    include_index: bool = False,
    **kwargs: Any,
) -> None:
    """Save a DataFrame to CSV.

    Parameters
    ----------
    df : pandas.DataFrame
        DataFrame to export.
    path : str or Path
        Destination file path.
    include_index : bool, optional
        Whether to include the DataFrame index. Defaults to ``False``.
    **kwargs
        Additional arguments passed to :func:`pandas.DataFrame.to_csv`.

    Raises
    ------
    FileWriteError
        If pandas or the filesystem cannot write the destination.
    """
    try:
        df.to_csv(path, index=include_index, **kwargs)
    except (OSError, ValueError, TypeError, UnicodeError) as exc:
        raise FileWriteError(str(path), exc) from exc

export_to_excel

export_to_excel(df: DataFrame, path: str | Path, include_index: bool = False, **kwargs: Any) -> None

Save a DataFrame to an Excel workbook.

Parameters:

Name	Type	Description	Default
df	DataFrame	DataFrame to export.	required
path	str or Path	Destination file path.	required
include_index	bool	Whether to include the index column. Defaults to False.	False
**kwargs	Any	Additional arguments passed to :func:pandas.DataFrame.to_excel.	{}

Raises:

Type	Description
FileWriteError	If pandas or the filesystem cannot write the destination.

Source code in src/industrialstats/utils/export.py

def export_to_excel(
    df: pd.DataFrame,
    path: str | Path,
    include_index: bool = False,
    **kwargs: Any,
) -> None:
    """Save a DataFrame to an Excel workbook.

    Parameters
    ----------
    df : pandas.DataFrame
        DataFrame to export.
    path : str or Path
        Destination file path.
    include_index : bool, optional
        Whether to include the index column. Defaults to ``False``.
    **kwargs
        Additional arguments passed to :func:`pandas.DataFrame.to_excel`.

    Raises
    ------
    FileWriteError
        If pandas or the filesystem cannot write the destination.
    """
    try:
        df.to_excel(path, index=include_index, **kwargs)
    except (OSError, ValueError, TypeError, UnicodeError) as exc:
        raise FileWriteError(str(path), exc) from exc

export_to_json

export_to_json(df: DataFrame, path: str | Path, **kwargs: Any) -> None

Save a DataFrame and metadata to JSON.

Parameters:

Name	Type	Description	Default
df	DataFrame	DataFrame to export.	required
path	str or Path	Destination file path.	required
**kwargs	Any	Additional JSON dump options.	{}

Raises:

Type	Description
FileWriteError	If serialization or filesystem writing fails.

Source code in src/industrialstats/utils/export.py

def export_to_json(df: pd.DataFrame, path: str | Path, **kwargs: Any) -> None:
    """Save a DataFrame and metadata to JSON.

    Parameters
    ----------
    df : pandas.DataFrame
        DataFrame to export.
    path : str or Path
        Destination file path.
    **kwargs
        Additional JSON ``dump`` options.

    Raises
    ------
    FileWriteError
        If serialization or filesystem writing fails.
    """
    data = {
        "data": df.to_dict(orient="records"),
        "columns": list(df.columns),
    }
    try:
        with Path(path).open("w", encoding="utf-8") as f:
            json.dump(data, f, indent=2, **kwargs)
    except (OSError, ValueError, TypeError, UnicodeError) as exc:
        raise FileWriteError(str(path), exc) from exc

Performance

industrialstats.utils.performance

Utilities for profiling code execution paths.

profile_function

profile_function(func: Callable[..., Any], *args: Any, **kwargs: Any) -> Stats

Profile a callable and return execution statistics.

Parameters:

Name	Type	Description	Default
func	Callable	Function or method to profile.	required
*args	Any	Positional arguments passed to func.	()
**kwargs	Any	Keyword arguments passed to func.	{}

Returns:

Type	Description
Stats	Profiling statistics sorted by cumulative time.

Source code in src/industrialstats/utils/performance.py

def profile_function(
    func: Callable[..., Any], *args: Any, **kwargs: Any
) -> pstats.Stats:
    """Profile a callable and return execution statistics.

    Parameters
    ----------
    func : Callable
        Function or method to profile.
    *args : Any
        Positional arguments passed to ``func``.
    **kwargs : Any
        Keyword arguments passed to ``func``.

    Returns
    -------
    pstats.Stats
        Profiling statistics sorted by cumulative time.
    """
    profiler = cProfile.Profile()
    profiler.enable()
    func(*args, **kwargs)
    profiler.disable()
    stats = pstats.Stats(profiler)
    stats.sort_stats("cumtime")
    return stats

Configuration

industrialstats.config

Global configuration management for industrialstats.

This module provides a simple configuration system that controls plotting preferences, numerical precision, and logging levels across the package. Configuration values can be loaded from a JSON or YAML file and are applied to relevant third-party libraries.

Config dataclass

Config(plot_style: str = 'ggplot', theme: str = 'whitegrid', precision: int = 4, log_level: str = 'INFO')

Package configuration settings.

Parameters:

Name	Type	Description	Default
plot_style	str	Matplotlib style to apply for plots.	"ggplot"
theme	str	Seaborn theme used to style figures.	"whitegrid"
precision	int	Number of decimal places for NumPy printing.	4
log_level	str	Logging level applied to the root logger.	"INFO"

apply

apply() -> None

Apply configuration to third-party libraries.

Source code in src/industrialstats/config.py

def apply(self) -> None:
    """Apply configuration to third-party libraries."""
    plt.style.use(self.plot_style)
    sns.set_theme(style=self.theme)
    np.set_printoptions(precision=self.precision)
    logging.getLogger().setLevel(self.log_level.upper())
    self._applied = True

update

update(**kwargs: Any) -> None

Update configuration values and reapply settings.

Parameters:

Name	Type	Description	Default
**kwargs	Any	Configuration fields to update.	{}

Source code in src/industrialstats/config.py

def update(self, **kwargs: Any) -> None:
    """Update configuration values and reapply settings.

    Parameters
    ----------
    **kwargs
        Configuration fields to update.
    """
    for key, value in kwargs.items():
        if hasattr(self, key):
            setattr(self, key, value)
    self.apply()

load_config

load_config(path: str | Path) -> None

Load configuration from a JSON or YAML file.

Parameters:

Name	Type	Description	Default
path	str or Path	Path to the configuration file.	required

Raises:

Type	Description
ValueError	If the file format is unsupported or PyYAML is required but not installed.

Source code in src/industrialstats/config.py

def load_config(path: str | Path) -> None:
    """Load configuration from a JSON or YAML file.

    Parameters
    ----------
    path : str or Path
        Path to the configuration file.

    Raises
    ------
    ValueError
        If the file format is unsupported or PyYAML is required but not
        installed.
    """
    path = Path(path)
    if path.suffix.lower() == ".json":
        data: dict[str, Any] = json.loads(path.read_text())
    elif path.suffix.lower() in {".yml", ".yaml"}:
        if yaml is None:  # pragma: no cover - handled above
            raise ValueError("PyYAML is required for YAML configuration files")
        data = yaml.safe_load(path.read_text())
    else:  # pragma: no cover - defensive
        raise ValueError("Unsupported configuration file format")

    config.update(**data)