131 lines
5.1 KiB
Python
131 lines
5.1 KiB
Python
# 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
|