From 33fd9f1177627d82fa7d244af4645b730ddde93a Mon Sep 17 00:00:00 2001 From: Artemiy Date: Fri, 24 Jul 2026 00:35:08 +0500 Subject: [PATCH] =?UTF-8?q?=D0=94=D0=BE=D0=B1=D0=B0=D0=B2=D0=BB=D0=B5?= =?UTF-8?q?=D0=BD=20=D0=B0=D0=BD=D0=B0=D0=BB=D0=B8=D0=B7=20ANOVA?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- analysis.py | 166 ++++++++++++++++++++++++++++++++++++++++++--- core/statistics.py | 80 ++++++++++++++++++++++ 2 files changed, 236 insertions(+), 10 deletions(-) diff --git a/analysis.py b/analysis.py index 92733c1..96d4959 100644 --- a/analysis.py +++ b/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]) - df_res = N - (k + 1) + X = np.hstack([np.ones((N, 1)), plan_matrix]) # размер N x (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"] - se = reg["se"] - t_stats = reg["t_stats"] - s2_res = reg["s2_res"] - df_res_actual = reg["df_res"] - y_pred = reg["y_pred"] + beta = reg["beta"] # коэффициенты (включая свободный член) + se = reg["se"] # стандартные ошибки коэффициентов + t_stats = reg["t_stats"] # t-статистики + s2_res = reg["s2_res"] # дисперсия остатков + df_res_actual = reg["df_res"] # степени свободы остатка (N - p) + 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 в планировании эксперимента.") diff --git a/core/statistics.py b/core/statistics.py index d6d8663..27bdc53 100644 --- a/core/statistics.py +++ b/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 + }