Initial commit: планирование эксперимента и оптимизация смесей

core/statistics.py

@@ -0,0 +1,130 @@
+# core/statistics.py
+"""
+Статистические функции для регрессии, проверки значимости, дисперсии и градиента.
+"""
+import numpy as np
+try:
+    from scipy import stats
+    SCIPY_AVAILABLE = True
+except ImportError:
+    SCIPY_AVAILABLE = False
+
+def compute_regression(X, y):
+    """МНК-регрессия, возвращает коэффициенты, ошибки, t-статистики."""
+    N, p = X.shape
+    try:
+        XTX = X.T @ X
+        XTX_inv = np.linalg.inv(XTX)
+        beta = XTX_inv @ (X.T @ y)
+        y_pred = X @ beta
+        residuals = y - y_pred
+        df_res = N - p
+        if df_res > 0:
+            s2_res = np.sum(residuals**2) / df_res
+            se = np.sqrt(np.diag(XTX_inv) * s2_res)
+        else:
+            s2_res = np.nan
+            se = [np.nan] * p
+        t_stats = beta / se if not np.isnan(se).any() else [np.nan] * p
+    except np.linalg.LinAlgError:
+        beta = [np.nan] * p
+        se = [np.nan] * p
+        t_stats = [np.nan] * p
+        df_res = 0
+        s2_res = np.nan
+        residuals = np.full(N, np.nan)
+        y_pred = np.full(N, np.nan)
+        XTX_inv = None
+    return {
+        "beta": beta,
+        "se": se,
+        "t_stats": t_stats,
+        "df_res": df_res,
+        "s2_res": s2_res,
+        "residuals": residuals,
+        "y_pred": y_pred,
+        "XTX_inv": XTX_inv
+    }
+
+def regression_summary(y, y_pred, k, df_res):
+    """Вычисляет R², скорректированный R², F-статистику и p-value."""
+    n = len(y)
+    ss_total = np.sum((y - np.mean(y))**2)
+    ss_res = np.sum((y - y_pred)**2)
+    r2 = 1 - ss_res / ss_total if ss_total > 0 else np.nan
+    if df_res > 0:
+        r2_adj = 1 - (1 - r2) * (n - 1) / df_res
+    else:
+        r2_adj = np.nan
+    if df_res > 0 and k > 0:
+        f_stat = (r2 / k) / ((1 - r2) / df_res) if r2 < 1 else np.inf
+    else:
+        f_stat = np.nan
+    if SCIPY_AVAILABLE and not np.isnan(f_stat) and f_stat != np.inf:
+        p_value = 1 - stats.f.cdf(f_stat, k, df_res)
+    else:
+        p_value = np.nan
+    return {"r2": r2, "r2_adj": r2_adj, "f_stat": f_stat, "f_pvalue": p_value}
+
+def get_center_indices(plan, active_indices):
+    """Возвращает индексы строк плана, где активные факторы равны 0."""
+    if not plan:
+        return []
+    plan_arr = np.array(plan)
+    if active_indices:
+        active_cols = plan_arr[:, active_indices]
+        center_mask = np.all(active_cols == 0, axis=1)
+        return np.where(center_mask)[0].tolist()
+    return []
+
+def compute_center_stats(y_center):
+    """Среднее и несмещённая дисперсия для центральных точек."""
+    if len(y_center) < 2:
+        return None, None
+    return np.mean(y_center), np.var(y_center, ddof=1)
+
+def chi2_confidence_interval(var, n, alpha=0.05):
+    """Доверительный интервал для дисперсии на основе χ²-распределения."""
+    if not SCIPY_AVAILABLE or n < 2:
+        return np.nan, np.nan
+    df = n - 1
+    chi2_lower = stats.chi2.ppf(alpha / 2, df)
+    chi2_upper = stats.chi2.ppf(1 - alpha / 2, df)
+    return df * var / chi2_upper, df * var / chi2_lower
+
+def get_t_critical(df, alpha=0.05):
+    """Критическое значение t (двустороннее)."""
+    if not SCIPY_AVAILABLE or df <= 0:
+        return np.nan
+    return stats.t.ppf(1 - alpha/2, df)
+
+def compute_gradient_points(beta, base_levels, step_factor, num_steps, factor_names=None, maximize=True):
+    """
+    Генерирует точки вдоль градиента (крутое восхождение или спуск).
+    beta – коэффициенты модели: [b0, b1, b2, ...]
+    base_levels – список базовых уровней для каждого фактора
+    step_factor – шаг по первому фактору
+    num_steps – количество шагов
+    factor_names – список названий факторов
+    maximize – True для максимизации, False для минимизации
+    """
+    if len(beta) - 1 != len(base_levels):
+        raise ValueError("Число коэффициентов (без b0) должно совпадать с числом факторов")
+    direction = 1 if maximize else -1
+    coefs = np.array(beta[1:]) * direction
+    max_abs = np.max(np.abs(coefs))
+    if max_abs == 0:
+        raise ValueError("Все коэффициенты равны нулю, направление не определено.")
+    norm_coefs = coefs / max_abs
+    steps = step_factor * norm_coefs
+
+    points = []
+    for step_num in range(1, num_steps + 1):
+        point = {}
+        for i, name in enumerate(factor_names or [f"x{i+1}" for i in range(len(base_levels))]):
+            point[name] = base_levels[i] + steps[i] * step_num
+        predicted = beta[0] + np.sum(beta[1:] * np.array([point[name] for name in factor_names or range(len(base_levels))]))
+        point["Шаг"] = step_num
+        point["Предсказанный отклик"] = predicted
+        points.append(point)
+    return points