Добавлен анализ ANOVA

analysis.py

@@ -1,7 +1,15 @@
 # analysis.py
 """
-Анализ экспериментальных данных: регрессия, проверка значимости, дисперсия в центре.
+Модуль анализа экспериментальных данных.
+
+Выполняет:
+  - множественную линейную регрессию (МНК),
+  - проверку значимости коэффициентов (t-критерий),
+  - дисперсионный анализ (ANOVA) с разложением остатка на чистую ошибку и неадекватность (если есть центр),
+  - оценку воспроизводимости по центральным точкам,
+  - исключение незначимых факторов с пересчётом модели.
 """
+
 import numpy as np
 from core.table_io import handle_table_output
 from core.statistics import (
@@ -10,15 +18,34 @@ from core.statistics import (
     compute_center_stats,
     chi2_confidence_interval,
     get_t_critical,
-    regression_summary
+    regression_summary,
+    anova_regression          # функция ANOVA, добавленная в statistics.py
 )
 
+
 def analysis_menu(mixture, design, responses):
+    """
+    Главное меню подсистемы анализа.
+
+    Аргументы:
+        mixture (Mixture): текущий состав смеси (нужен для имён факторов).
+        design (ExperimentDesign): объект с шагами и матрицей плана.
+        responses (list): список откликов, каждый как {'name': str, 'values': list}.
+
+    Действия:
+        - проверяет наличие плана и откликов;
+        - определяет активные факторы (шаг > 0);
+        - запрашивает уровень значимости alpha;
+        - строит матрицу X (единичный столбец + активные факторы);
+        - позволяет выбрать один отклик или все для последовательного анализа.
+    """
+    # Проверка наличия данных
     if not design.plan or not responses:
         print("Нет данных для анализа: план или отклики пусты.")
         input("Нажмите Enter для продолжения...")
         return
 
+    # Активные факторы – те, для которых задан шаг > 0
     active_indices = design._active_indices()
     k = len(active_indices)
     if k < 2:
@@ -26,6 +53,7 @@ def analysis_menu(mixture, design, responses):
         input("Нажмите Enter для продолжения...")
         return
 
+    # Имена активных факторов берём из смеси (индекс = idx + 1, т.к. 0 – растворитель)
     factor_names = []
     for idx in active_indices:
         ing_index = idx + 1
@@ -34,6 +62,7 @@ def analysis_menu(mixture, design, responses):
         else:
             factor_names.append(f"x{idx+1}")
 
+    # Ввод уровня значимости alpha
     print("\nВведите уровень значимости alpha (например, 0.05, 0.01, 0.001):")
     while True:
         try:
@@ -49,11 +78,13 @@ def analysis_menu(mixture, design, responses):
         except ValueError:
             print("Введите число.")
 
+    # Матрица X: столбец единиц + столбцы активных факторов (в кодировке ±1)
     plan_matrix = np.array([[row[j] for j in active_indices] for row in design.plan])
     N = plan_matrix.shape[0]
-    X = np.hstack([np.ones((N, 1)), plan_matrix])
+    X = np.hstack([np.ones((N, 1)), plan_matrix])   # размер N x (k+1)
-    df_res = N - (k + 1)
+    df_res = N - (k + 1)   # степени свободы остатка для полной модели
 
+    # Цикл выбора отклика
     while True:
         print("\n" + "=" * 50)
         print("АНАЛИЗ ЗНАЧЕНИЙ")
@@ -67,15 +98,19 @@ def analysis_menu(mixture, design, responses):
         choice = input("Выберите отклик для анализа: ").strip()
         if choice == "":
             break
+
         try:
             idx = int(choice) - 1
             if idx < 0:
                 print("Неверный номер.")
                 continue
+
             if idx < len(responses):
+                # Анализ одного отклика
                 analyze_single_response(responses[idx], X, k, df_res, design, alpha, factor_names)
                 input("Нажмите Enter для продолжения...")
             elif idx == len(responses):
+                # Анализ всех откликов подряд
                 for resp in responses:
                     print(f"\n--- Анализ отклика '{resp['name']}' ---")
                     analyze_single_response(resp, X, k, df_res, design, alpha, factor_names)
@@ -85,26 +120,57 @@ def analysis_menu(mixture, design, responses):
         except ValueError:
             print("Введите число.")
 
+
 def analyze_single_response(response, X, k, df_res, design, alpha, factor_names):
+    """
+    Проводит полный анализ одного отклика.
+
+    Этапы:
+      1. Расчёт регрессии методом наименьших квадратов (МНК).
+      2. Вывод регрессионной статистики: R², скорректированный R², F-критерий,
+         коэффициенты, их стандартные ошибки и t-статистики.
+      3. Дисперсионный анализ (ANOVA) с разложением сумм квадратов.
+         Если есть центральные точки (>=2), остаток разбивается на чистую ошибку
+         и неадекватность (lack-of-fit).
+      4. Анализ центральных точек – оценка воспроизводимости.
+      5. Опция исключения незначимых факторов (p > alpha) и пересчёт упрощённой модели.
+
+    Аргументы:
+        response (dict): {'name': str, 'values': list}
+        X (np.ndarray): матрица регрессоров (с единичным столбцом)
+        k (int): число активных факторов
+        df_res (int): степени свободы остатка (N - k - 1)
+        design (ExperimentDesign): объект с планом
+        alpha (float): уровень значимости
+        factor_names (list): имена факторов
+    """
+    # Преобразуем значения отклика в массив numpy
     y = np.array(response["values"])
     N = len(y)
     if N != X.shape[0]:
         print(f"Ошибка: число значений отклика ({N}) не соответствует числу опытов ({X.shape[0]})")
         return
 
+    # ---------- Шаг 1: Расчёт регрессии ----------
     reg = compute_regression(X, y)
-    beta = reg["beta"]
+    beta = reg["beta"]          # коэффициенты (включая свободный член)
-    se = reg["se"]
+    se = reg["se"]              # стандартные ошибки коэффициентов
-    t_stats = reg["t_stats"]
+    t_stats = reg["t_stats"]    # t-статистики
-    s2_res = reg["s2_res"]
+    s2_res = reg["s2_res"]      # дисперсия остатков
-    df_res_actual = reg["df_res"]
+    df_res_actual = reg["df_res"] # степени свободы остатка (N - p)
-    y_pred = reg["y_pred"]
+    y_pred = reg["y_pred"]      # предсказанные значения модели
 
+    # Проверка на вырожденность матрицы
     if np.isnan(beta).any():
         print("Ошибка при расчёте регрессии: матрица вырождена.")
         return
 
+    # ---------- Шаг 2: Вспомогательная функция для вывода регрессионной статистики ----------
     def print_and_get_tcrit(beta, se, t_stats, df_res, y_pred, y, k, factor_names, alpha):
+        """
+        Выводит основную регрессионную статистику и возвращает критическое t.
+        """
+        # Расчёт R², скорр. R² и F-критерия
         summary = regression_summary(y, y_pred, k, df_res)
         r2 = summary["r2"]
         r2_adj = summary["r2_adj"]
@@ -115,13 +181,16 @@ def analyze_single_response(response, X, k, df_res, design, alpha, factor_names)
         print(f"Число опытов: {len(y)}")
         print(f"Число факторов: {k}")
         print(f"Степени свободы остатков: {df_res}")
+        # Стандартная ошибка регрессии (корень из дисперсии остатков)
         print(f"Стандартная ошибка регрессии: {np.sqrt(s2_res) if not np.isnan(s2_res) else np.nan:.6f}")
         print(f"Множественный R: {np.sqrt(r2) if not np.isnan(r2) else np.nan:.6f}")
         print(f"R-квадрат: {r2:.6f}")
         print(f"Скорректированный R-квадрат: {r2_adj:.6f}")
         print(f"F-статистика: {f_stat:.4f}, p-value: {f_pvalue:.6f}")
 
+        # Критическое значение t для двустороннего критерия
         t_crit = get_t_critical(df_res, alpha=alpha)
+
         print("\nКоэффициенты регрессии:")
         headers = ["Коэффициент", "Значение", "Станд. ошибка", "t-статистика", f"|t| > t_крит({alpha})?"]
         rows = []
@@ -135,14 +204,77 @@ def analyze_single_response(response, X, k, df_res, design, alpha, factor_names)
                 f"{t_stats[i]:.4f}",
                 sig
             ])
+        # Вывод таблицы с возможностью сохранения/копирования
         handle_table_output(headers, rows, title="", ask_each_time=True)
+
         if not np.isnan(t_crit):
             print(f"Критическое значение t (alpha={alpha}, df={df_res}): {t_crit:.4f}")
         return t_crit
 
+    # ---------- Шаг 2 (выполнение): вывод статистики для полной модели ----------
     t_crit = print_and_get_tcrit(beta, se, t_stats, df_res_actual, y_pred, y, k, factor_names, alpha)
 
+    # ---------- Шаг 3: Дисперсионный анализ (ANOVA) ----------
+    # Определяем индексы строк плана, соответствующие центральным точкам (все активные факторы = 0)
     center_indices = get_center_indices(design.plan, design._active_indices())
+
+    # Выполняем ANOVA (если центр есть, разлагаем остаток)
+    anova = anova_regression(y, y_pred, X, center_indices if len(center_indices) >= 2 else None)
+
+    print("\n=== Дисперсионный анализ (ANOVA) ===")
+    anova_headers = ["Источник", "SS", "df", "MS", "F", "p-value"]
+    anova_rows = []
+
+    # Регрессия
+    anova_rows.append([
+        "Регрессия",
+        f"{anova['SSR']:.6f}",
+        f"{anova['df_reg']}",
+        f"{anova['MSR']:.6f}",
+        f"{anova['F_reg']:.4f}",
+        f"{anova['p_reg']:.6f}"
+    ])
+    # Остаток
+    anova_rows.append([
+        "Остаток",
+        f"{anova['SSE']:.6f}",
+        f"{anova['df_res']}",
+        f"{anova['MSE']:.6f}",
+        "-",
+        "-"
+    ])
+
+    # Если есть центральные точки (>=2), добавляем разложение остатка
+    if len(center_indices) >= 2:
+        anova_rows.append([
+            "  Чистая ошибка",
+            f"{anova['SSPE']:.6f}",
+            f"{anova['df_pe']}",
+            f"{anova['MS_pe']:.6f}",
+            "-",
+            "-"
+        ])
+        anova_rows.append([
+            "  Неадекватность",
+            f"{anova['SSLOF']:.6f}",
+            f"{anova['df_lof']}",
+            f"{anova['MS_lof']:.6f}",
+            f"{anova['F_lof']:.4f}",
+            f"{anova['p_lof']:.6f}"
+        ])
+    # Итого (общая сумма квадратов)
+    anova_rows.append([
+        "Итого",
+        f"{anova['SST']:.6f}",
+        f"{len(y)-1}",
+        "-",
+        "-",
+        "-"
+    ])
+
+    handle_table_output(anova_headers, anova_rows, title="", ask_each_time=True)
+
+    # ---------- Шаг 4: Анализ центральных точек (воспроизводимость) ----------
     if len(center_indices) >= 2:
         y_center = y[center_indices]
         center_mean, center_var = compute_center_stats(y_center)
@@ -151,32 +283,46 @@ def analyze_single_response(response, X, k, df_res, design, alpha, factor_names)
         print(f"Среднее значение отклика в центре: {center_mean:.6f}")
         print(f"Дисперсия в центре (несмещённая): {center_var:.6f}")
         print(f"Стандартное отклонение: {np.sqrt(center_var):.6f}")
+
+        # Доверительный интервал для дисперсии (хи-квадрат)
         lower, upper = chi2_confidence_interval(center_var, len(center_indices), alpha=alpha)
         if not np.isnan(lower):
             print(f"Доверительный интервал для дисперсии (уровень доверия {1-alpha}): [{lower:.6f}, {upper:.6f}]")
     else:
         print("\nНедостаточно центральных точек (нужно минимум 2).")
 
+    # ---------- Шаг 5: Опция исключения незначимых факторов ----------
     if k > 0 and not np.isnan(t_crit):
         print("\nХотите исключить незначимые факторы и пересчитать модель?")
         choice = input("(y/n, по умолчанию n): ").strip().lower()
         if choice == 'y':
+            # Маска: первый столбец (свободный член) всегда включаем
             significant = [True] + [abs(t_stats[i]) > t_crit for i in range(1, k+1)]
+
+            # Если все факторы значимы, упрощённая модель не отличается
             if all(significant):
                 print("Все факторы значимы. Упрощённая модель не отличается от полной.")
                 return
+
+            # Выбираем только значимые столбцы
             X_new = X[:, significant]
+            # Пересчитываем регрессию
             reg_new = compute_regression(X_new, y)
             beta_new = reg_new["beta"]
             se_new = reg_new["se"]
             t_stats_new = reg_new["t_stats"]
             df_res_new = reg_new["df_res"]
             y_pred_new = reg_new["y_pred"]
+
             if np.isnan(beta_new).any():
                 print("Ошибка при пересчёте упрощённой модели.")
                 return
+
+            # Обновляем список имён факторов
             factor_names_new = ["Свободный член"] + [factor_names[i-1] for i in range(1, k+1) if significant[i]]
             k_new = len(factor_names_new) - 1
+
             print("\n=== Упрощённая модель (исключены незначимые факторы) ===")
             print_and_get_tcrit(beta_new, se_new, t_stats_new, df_res_new, y_pred_new, y, k_new, factor_names_new[1:], alpha)
+
             print("\nСовет: для использования упрощённой модели в крутом восхождении установите шаги для исключённых факторов равными 0 в планировании эксперимента.")

core/statistics.py

@@ -128,3 +128,83 @@ def compute_gradient_points(beta, base_levels, step_factor, num_steps, factor_na
         point["Предсказанный отклик"] = predicted
         points.append(point)
     return points
+
+
+def anova_regression(y, y_pred, X, center_indices=None):
+    """
+    Выполняет полный дисперсионный анализ (ANOVA) для регрессионной модели.
+
+    Аргументы:
+        y – фактические значения отклика
+        y_pred – предсказанные моделью значения
+        X – матрица регрессоров (с единичным столбцом)
+        center_indices – индексы центральных точек (опционально)
+
+    Возвращает словарь с компонентами ANOVA:
+        SSR, SSE, SST, df_reg, df_res, df_pe, df_lof,
+        MSR, MSE, MS_pe, MS_lof,
+        F_regression, p_regression, F_lof, p_lof
+    """
+    n = len(y)
+    p = X.shape[1]          # число параметров (включая свободный член)
+    df_reg = p - 1          # степени свободы регрессии (без свободного члена)
+    df_res = n - p          # степени свободы остатка
+
+    # Общая сумма квадратов (SST)
+    y_mean = np.mean(y)
+    SST = np.sum((y - y_mean)**2)
+
+    # Регрессионная сумма квадратов (SSR)
+    SSR = np.sum((y_pred - y_mean)**2)
+
+    # Остаточная сумма квадратов (SSE)
+    SSE = np.sum((y - y_pred)**2)
+
+    # Средние квадраты
+    MSR = SSR / df_reg if df_reg > 0 else np.nan
+    MSE = SSE / df_res if df_res > 0 else np.nan
+
+    # F-статистика для регрессии (значимость всей модели)
+    F_reg = MSR / MSE if MSE > 0 else np.nan
+    p_reg = 1 - stats.f.cdf(F_reg, df_reg, df_res) if not np.isnan(F_reg) and SCIPY_AVAILABLE else np.nan
+
+    # ----- Разложение остатка на чистую ошибку и неадекватность (если есть центр) -----
+    if center_indices is not None and len(center_indices) >= 2:
+        y_center = y[center_indices]
+        y_pred_center = y_pred[center_indices]
+        n_center = len(center_indices)
+        # Среднее в центре
+        y_center_mean = np.mean(y_center)
+        # Сумма квадратов чистой ошибки (SSPE) = сумма квадратов отклонений от среднего в центре
+        SSPE = np.sum((y_center - y_center_mean)**2)
+        df_pe = n_center - 1
+        MS_pe = SSPE / df_pe if df_pe > 0 else np.nan
+
+        # Сумма квадратов неадекватности (SSLOF) = SSE - SSPE
+        SSLOF = SSE - SSPE
+        df_lof = df_res - df_pe
+        MS_lof = SSLOF / df_lof if df_lof > 0 else np.nan
+
+        # F-статистика для проверки адекватности (lack-of-fit)
+        F_lof = MS_lof / MS_pe if MS_pe > 0 and not np.isnan(MS_lof) else np.nan
+        p_lof = 1 - stats.f.cdf(F_lof, df_lof, df_pe) if not np.isnan(F_lof) and SCIPY_AVAILABLE else np.nan
+    else:
+        SSPE = np.nan
+        df_pe = 0
+        MS_pe = np.nan
+        SSLOF = np.nan
+        df_lof = 0
+        MS_lof = np.nan
+        F_lof = np.nan
+        p_lof = np.nan
+
+    return {
+        "SSR": SSR, "SSE": SSE, "SST": SST,
+        "df_reg": df_reg, "df_res": df_res,
+        "MSR": MSR, "MSE": MSE,
+        "F_reg": F_reg, "p_reg": p_reg,
+        "SSPE": SSPE, "SSLOF": SSLOF,
+        "df_pe": df_pe, "df_lof": df_lof,
+        "MS_pe": MS_pe, "MS_lof": MS_lof,
+        "F_lof": F_lof, "p_lof": p_lof
+    }