Добавлены две функции для анализа ЦКП плана: build_quadratic_X и analyze_quadratic_model

core/statistics.py

@@ -208,3 +208,103 @@ def anova_regression(y, y_pred, X, center_indices=None):
         "MS_pe": MS_pe, "MS_lof": MS_lof,
         "F_lof": F_lof, "p_lof": p_lof
     }
+# core/statistics.py (дополнение)
+
+def build_quadratic_X(plan):
+    """
+    Строит расширенную матрицу регрессоров для квадратичной модели.
+    plan: numpy array (N x k) – матрица плана в кодированных переменных.
+    Возвращает: X (N x p), где p = 1 + 2k + k*(k-1)//2
+    (свободный член, линейные, квадраты, попарные произведения)
+    """
+    N, k = plan.shape
+    X = np.ones((N, 1))               # свободный член
+    X = np.hstack([X, plan])          # линейные (k столбцов)
+    X = np.hstack([X, plan**2])       # квадраты (k столбцов)
+    # взаимодействия (k*(k-1)//2 столбцов)
+    for i in range(k):
+        for j in range(i+1, k):
+            X = np.hstack([X, (plan[:, i] * plan[:, j]).reshape(-1, 1)])
+    return X
+
+def fit_quadratic_model(plan, y):
+    """
+    Выполняет регрессию для квадратичной модели по матрице плана.
+    Возвращает словарь с результатами (аналогично compute_regression, но с добавлением структуры).
+    """
+    X = build_quadratic_X(plan)
+    reg = compute_regression(X, y)   # используем существующую функцию
+    # Добавляем информацию о структуре модели
+    k = plan.shape[1]
+    p_total = X.shape[1]
+    reg["k"] = k
+    reg["X"] = X
+    reg["plan"] = plan
+    reg["p_total"] = p_total
+    # Можно добавить имена коэффициентов (для вывода)
+    coef_names = ["Свободный член"]
+    coef_names += [f"x{i+1}" for i in range(k)]
+    coef_names += [f"x{i+1}²" for i in range(k)]
+    for i in range(k):
+        for j in range(i+1, k):
+            coef_names.append(f"x{i+1}·x{j+1}")
+    reg["coef_names"] = coef_names
+    return reg
+
+def analyze_quadratic_model(beta, k, X, y):
+    """
+    Анализирует квадратичную модель: находит стационарную точку и тип экстремума.
+    beta – коэффициенты (включая все члены модели)
+    k – число факторов
+    X – матрица регрессоров (для расчёта предсказанных значений)
+    y – отклик
+    Возвращает словарь с результатами.
+    """
+    # Извлекаем коэффициенты
+    b0 = beta[0]
+    b_lin = beta[1:1+k]                # линейные
+    b_quad = beta[1+k:1+2*k]           # квадраты
+    b_inter = beta[1+2*k:]             # взаимодействия
+
+    # Строим матрицу B (гессиан) для квадратичной формы
+    B = np.diag(2 * b_quad)
+    idx = 0
+    for i in range(k):
+        for j in range(i+1, k):
+            B[i, j] = b_inter[idx]
+            B[j, i] = b_inter[idx]
+            idx += 1
+
+    # Стационарная точка: x* = -0.5 * inv(B) * b_lin
+    try:
+        x_star = -0.5 * np.linalg.solve(B, b_lin)
+    except np.linalg.LinAlgError:
+        x_star = None
+
+    # Предсказанное значение в стационарной точке
+    if x_star is not None:
+        y_star = b0 + np.dot(b_lin, x_star) + np.dot(x_star, B @ x_star) / 2
+    else:
+        y_star = None
+
+    # Тип экстремума
+    if x_star is not None:
+        eigvals = np.linalg.eigvalsh(B)
+        if np.all(eigvals < 0):
+            extremum_type = "Максимум"
+        elif np.all(eigvals > 0):
+            extremum_type = "Минимум"
+        elif np.any(eigvals < 0) and np.any(eigvals > 0):
+            extremum_type = "Седловая точка"
+        else:
+            extremum_type = "Неопределён (есть нулевые собственные значения)"
+    else:
+        extremum_type = "Не найден (матрица B вырождена)"
+
+    return {
+        "x_star": x_star,
+        "y_star": y_star,
+        "extremum_type": extremum_type,
+        "B": B,
+        "eigenvalues": np.linalg.eigvalsh(B) if x_star is not None else None
+    }