365 lines
18 KiB
Python
365 lines
18 KiB
Python
# analysis.py
|
|
"""
|
|
Модуль анализа экспериментальных данных.
|
|
|
|
Выполняет:
|
|
- множественную линейную регрессию (МНК),
|
|
- проверку значимости коэффициентов (t-критерий),
|
|
- дисперсионный анализ (ANOVA) с разложением остатка на чистую ошибку и неадекватность (если есть центр),
|
|
- оценку воспроизводимости по центральным точкам,
|
|
- исключение незначимых факторов с пересчётом модели,
|
|
- интерпретацию результатов на естественном языке.
|
|
"""
|
|
|
|
import numpy as np
|
|
from core.table_io import handle_table_output
|
|
from core.interpreter import full_interpretation
|
|
from core.statistics import (
|
|
compute_regression,
|
|
get_center_indices,
|
|
compute_center_stats,
|
|
chi2_confidence_interval,
|
|
get_t_critical,
|
|
regression_summary,
|
|
anova_regression
|
|
)
|
|
|
|
|
|
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:
|
|
print("Для анализа необходимо минимум 2 активных фактора.")
|
|
input("Нажмите Enter для продолжения...")
|
|
return
|
|
|
|
# Имена активных факторов берём из смеси (индекс = idx + 1, т.к. 0 – растворитель)
|
|
factor_names = []
|
|
for idx in active_indices:
|
|
ing_index = idx + 1
|
|
if ing_index < len(mixture.ings):
|
|
factor_names.append(mixture.ings[ing_index].name)
|
|
else:
|
|
factor_names.append(f"x{idx+1}")
|
|
|
|
# Ввод уровня значимости alpha
|
|
print("\nВведите уровень значимости alpha (например, 0.05, 0.01, 0.001):")
|
|
while True:
|
|
try:
|
|
alpha_input = input("alpha (по умолчанию 0.05): ").strip()
|
|
if alpha_input == "":
|
|
alpha = 0.05
|
|
else:
|
|
alpha = float(alpha_input)
|
|
if alpha <= 0 or alpha >= 1:
|
|
print("alpha должна быть в интервале (0, 1).")
|
|
continue
|
|
break
|
|
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]) # размер N x (k+1)
|
|
df_res = N - (k + 1) # степени свободы остатка для полной модели
|
|
|
|
# Цикл выбора отклика
|
|
while True:
|
|
print("\n" + "=" * 50)
|
|
print("АНАЛИЗ ЗНАЧЕНИЙ")
|
|
print(f"Уровень значимости α = {alpha}")
|
|
print("Доступные отклики:")
|
|
for idx, resp in enumerate(responses, start=1):
|
|
print(f"{idx}. {resp['name']}")
|
|
print(f"{len(responses)+1}. Анализировать все отклики")
|
|
print("(Enter - выход в главное меню)")
|
|
print("=" * 50)
|
|
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)
|
|
input("Нажмите Enter для продолжения...")
|
|
else:
|
|
print("Неверный номер.")
|
|
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) и пересчёт упрощённой модели.
|
|
6. Интерпретация результатов на естественном языке.
|
|
|
|
Аргументы:
|
|
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"] # 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"]
|
|
f_stat = summary["f_stat"]
|
|
f_pvalue = summary["f_pvalue"]
|
|
|
|
print("\n=== Регрессионный анализ ===")
|
|
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 = []
|
|
coef_names = ["Свободный член"] + factor_names
|
|
for i, name in enumerate(coef_names):
|
|
sig = "Да" if (not np.isnan(t_stats[i]) and abs(t_stats[i]) > t_crit and not np.isnan(t_crit)) else "Нет"
|
|
rows.append([
|
|
name,
|
|
f"{beta[i]:.6f}",
|
|
f"{se[i]:.6f}",
|
|
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}")
|
|
|
|
# Возвращаем summary для дальнейшего использования
|
|
return t_crit, summary
|
|
|
|
# ---------- Шаг 2 (выполнение): вывод статистики для полной модели ----------
|
|
t_crit, summary_full = print_and_get_tcrit(beta, se, t_stats, df_res_actual, y_pred, y, k, factor_names, alpha)
|
|
|
|
# Извлекаем показатели качества
|
|
r2 = summary_full["r2"]
|
|
r2_adj = summary_full["r2_adj"]
|
|
f_stat = summary_full["f_stat"]
|
|
f_pvalue = summary_full["f_pvalue"]
|
|
|
|
# ---------- Шаг 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: Анализ центральных точек (воспроизводимость) ----------
|
|
# Инициализируем переменные для центра
|
|
center_mean = None
|
|
center_var = None
|
|
|
|
if len(center_indices) >= 2:
|
|
y_center = y[center_indices]
|
|
center_mean, center_var = compute_center_stats(y_center)
|
|
print("\n=== Анализ центральных точек ===")
|
|
print(f"Количество центральных точек: {len(center_indices)}")
|
|
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("Все факторы значимы. Упрощённая модель не отличается от полной.")
|
|
else:
|
|
# Выбираем только значимые столбцы
|
|
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("Ошибка при пересчёте упрощённой модели.")
|
|
else:
|
|
# Обновляем список имён факторов
|
|
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 в планировании эксперимента.")
|
|
|
|
# ---------- Шаг 6: Интерпретация результатов ----------
|
|
# Создаём полный список имён (включая свободный член) для интерпретатора
|
|
factor_names_full = ["Свободный член"] + factor_names
|
|
|
|
# Формируем полное заключение с помощью интерпретатора
|
|
interpretation = full_interpretation(
|
|
r2=r2,
|
|
r2_adj=r2_adj,
|
|
f_stat=f_stat,
|
|
f_pvalue=f_pvalue,
|
|
beta=beta,
|
|
t_stats=t_stats,
|
|
p_values=None, # пока не вычисляем p-values для коэффициентов
|
|
factor_names=factor_names_full, # <-- передаём полный список
|
|
anova=anova,
|
|
center_mean=center_mean,
|
|
center_var=center_var,
|
|
n_center=len(center_indices) if len(center_indices) >= 2 else 0,
|
|
alpha=alpha
|
|
)
|
|
|
|
print("\n" + interpretation)
|