import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from statsmodels.stats.stattools import durbin_watson
from scipy import stats
import seaborn as sns
import matplotlib.pyplot as plt
import statsmodels.api as sm

# === Загружаем данные ===
columns = ['Position', 'Height(inches)', 'Weight(pounds)', 'Age']
df = pd.read_csv('SOCR_Data_MLB_HeightsWeights.txt', sep='\t', nrows=115, usecols=columns)

# Кодируем категориальный признак "Position"
df['Position_num'] = pd.factorize(df['Position'])[0]

print("=" * 60)
print("ПАРНАЯ РЕГРЕССИЯ: Weight = f(Height)")
print("=" * 60)

# === ПАРНАЯ РЕГРЕССИЯ (Weight = f(Height)) ===
X_simple = df[['Height(inches)']]
y_simple = df['Weight(pounds)']

model_simple = LinearRegression().fit(X_simple, y_simple)
y_pred_simple = model_simple.predict(X_simple)

# Параметры модели
r_squared_simple = model_simple.score(X_simple, y_simple)
n_simple = len(y_simple)
adj_r_squared_simple = 1 - (1 - r_squared_simple) * (n_simple - 1) / (n_simple - 2)
std_error_simple = np.std(y_simple - y_pred_simple)
multiple_r_simple = np.sqrt(r_squared_simple)

print('ПАРАМЕТРЫ ПАРНОЙ РЕГРЕССИИ:')
print(f'Коэффициент детерминации R²: {r_squared_simple:.4f}')
print(f'Скорректированный R²: {adj_r_squared_simple:.4f}')
print(f'Множественный коэффициент корреляции R: {multiple_r_simple:.4f}')
print(f'Стандартная ошибка регрессии: {std_error_simple:.4f}')
print(f'Коэффициенты модели (b0 + b1*x): {model_simple.intercept_:.4f} [{model_simple.coef_[0]:.4f}]')

# === Остатки парной регрессии ===
residuals_simple = y_simple - y_pred_simple

# === ПРОВЕРКА НОРМАЛЬНОСТИ ОСТАТКОВ (Парная регрессия) ===
print('\nПРОВЕРКА НОРМАЛЬНОСТИ ОСТАТКОВ (Парная регрессия):')

# Тест Шапиро-Уилка
shapiro_stat_simple, shapiro_p_simple = stats.shapiro(residuals_simple)
print(f'Тест Шапиро-Уилка:')
print(f'Статистика: {shapiro_stat_simple:.4f}, p-value = {shapiro_p_simple:.4f}')
if shapiro_p_simple > 0.05:
    print('→ Остатки распределены нормально (H₀ не отвергается)')
else:
    print('→ Остатки не распределены нормально (H₀ отвергается)')

# Тест Колмогорова-Смирнова
ks_stat_simple, ks_p_simple = stats.kstest(residuals_simple, 'norm', args=(np.mean(residuals_simple), np.std(residuals_simple)))
print(f'\nТест Колмогорова-Смирнова:')
print(f'Статистика: {ks_stat_simple:.4f}, p-value = {ks_p_simple:.4f}')
if ks_p_simple > 0.05:
    print('→ Остатки распределены нормально (H₀ не отвергается)')
else:
    print('→ Остатки не распределены нормально (H₀ отвергается)')

# === СТАТИСТИКА ДУРБИНА-УОТСОНА (Парная регрессия) ===
dw_stat_simple = durbin_watson(residuals_simple)
print(f'\nСтатистика Дурбина–Уотсона: {dw_stat_simple:.4f}')

if dw_stat_simple < 1.5:
    print('→ Наблюдается ПОЛОЖИТЕЛЬНАЯ автокорреляция остатков')
elif dw_stat_simple > 2.5:
    print('→ Наблюдается ОТРИЦАТЕЛЬНАЯ автокорреляция остатков')
else:
    print('→ Автокорреляция остатков ОТСУТСТВУЕТ (остатки независимы)')

# === ГРАФИКИ ДЛЯ ПАРНОЙ РЕГРЕССИИ ===
fig1, axes1 = plt.subplots(2, 2, figsize=(15, 12))
fig1.suptitle('АНАЛИЗ ОСТАТКОВ: ПАРНАЯ РЕГРЕССИЯ (Weight = f(Height))', fontsize=16)

# 1. Гистограмма
axes1[0, 0].hist(residuals_simple, bins=15, density=True, alpha=0.7, color='lightblue', edgecolor='black')
xmin, xmax = axes1[0, 0].get_xlim()
x = np.linspace(xmin, xmax, 100)
p = stats.norm.pdf(x, np.mean(residuals_simple), np.std(residuals_simple))
axes1[0, 0].plot(x, p, 'k', linewidth=2, label='Нормальное распределение')
axes1[0, 0].set_title('Гистограмма остатков')
axes1[0, 0].set_xlabel('Остатки')
axes1[0, 0].set_ylabel('Плотность')
axes1[0, 0].legend()
axes1[0, 0].grid(True, alpha=0.3)

# 2. QQ-plot
sm.qqplot(residuals_simple, line='s', ax=axes1[0, 1])
axes1[0, 1].set_title('QQ-график остатков')
axes1[0, 1].grid(True, alpha=0.3)

# 3. Остатки vs Предсказанные значения
axes1[1, 0].scatter(y_pred_simple, residuals_simple, alpha=0.7, color='blue')
axes1[1, 0].axhline(y=0, color='red', linestyle='--')
axes1[1, 0].set_xlabel('Предсказанные значения')
axes1[1, 0].set_ylabel('Остатки')
axes1[1, 0].set_title('Остатки vs Предсказанные значения')
axes1[1, 0].grid(True, alpha=0.3)

# 4. Фактические vs Предсказанные значения
axes1[1, 1].scatter(y_simple, y_pred_simple, alpha=0.7, color='green')
min_val = min(y_simple.min(), y_pred_simple.min())
max_val = max(y_simple.max(), y_pred_simple.max())
axes1[1, 1].plot([min_val, max_val], [min_val, max_val], 'r--', label='Идеальная линия')
axes1[1, 1].set_xlabel('Фактические значения')
axes1[1, 1].set_ylabel('Предсказанные значения')
axes1[1, 1].set_title('Фактические vs Предсказанные значения')
axes1[1, 1].legend()
axes1[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("\n" + "=" * 60)
print("МНОЖЕСТВЕННАЯ РЕГРЕССИЯ: Weight = f(Position, Height, Age)")
print("=" * 60)

# === МНОЖЕСТВЕННАЯ РЕГРЕССИЯ (Weight = f(Position, Height, Age)) ===
X_multi = df[['Position_num', 'Height(inches)', 'Age']]
y_multi = df['Weight(pounds)']

model_multi = LinearRegression().fit(X_multi, y_multi)
y_pred_multi = model_multi.predict(X_multi)

# Параметры модели
r_squared_multi = model_multi.score(X_multi, y_multi)
n_multi = len(y_multi)
p_multi = X_multi.shape[1]
adj_r_squared_multi = 1 - (1 - r_squared_multi) * (n_multi - 1) / (n_multi - p_multi - 1)
std_error_multi = np.std(y_multi - y_pred_multi)
multiple_r_multi = np.sqrt(r_squared_multi)

print('ПАРАМЕТРЫ МНОЖЕСТВЕННОЙ РЕГРЕССИИ:')
print(f'Коэффициент детерминации R²: {r_squared_multi:.4f}')
print(f'Скорректированный R²: {adj_r_squared_multi:.4f}')
print(f'Множественный коэффициент корреляции R: {multiple_r_multi:.4f}')
print(f'Стандартная ошибка регрессии: {std_error_multi:.4f}')
print(f'Коэффициенты модели: {model_multi.intercept_:.4f} {model_multi.coef_}')

# === Остатки множественной регрессии ===
residuals_multi = y_multi - y_pred_multi

# === ПРОВЕРКА НОРМАЛЬНОСТИ ОСТАТКОВ (Множественная регрессия) ===
print('\nПРОВЕРКА НОРМАЛЬНОСТИ ОСТАТКОВ (Множественная регрессия):')

# Тест Шапиро-Уилка
shapiro_stat_multi, shapiro_p_multi = stats.shapiro(residuals_multi)
print(f'Тест Шапиро-Уилка:')
print(f'Статистика: {shapiro_stat_multi:.4f}, p-value = {shapiro_p_multi:.4f}')
if shapiro_p_multi > 0.05:
    print('→ Остатки распределены нормально (H₀ не отвергается)')
else:
    print('→ Остатки не распределены нормально (H₀ отвергается)')

# Тест Колмогорова-Смирнова
ks_stat_multi, ks_p_multi = stats.kstest(residuals_multi, 'norm', args=(np.mean(residuals_multi), np.std(residuals_multi)))
print(f'\nТест Колмогорова-Смирнова:')
print(f'Статистика: {ks_stat_multi:.4f}, p-value = {ks_p_multi:.4f}')
if ks_p_multi > 0.05:
    print('→ Остатки распределены нормально (H₀ не отвергается)')
else:
    print('→ Остатки не распределены нормально (H₀ отвергается)')

# === СТАТИСТИКА ДУРБИНА-УОТСОНА (Множественная регрессия) ===
dw_stat_multi = durbin_watson(residuals_multi)
print(f'\nСтатистика Дурбина–Уотсона: {dw_stat_multi:.4f}')

if dw_stat_multi < 1.5:
    print('→ Наблюдается ПОЛОЖИТЕЛЬНАЯ автокорреляция остатков')
elif dw_stat_multi > 2.5:
    print('→ Наблюдается ОТРИЦАТЕЛЬНАЯ автокорреляция остатков')
else:
    print('→ Автокорреляция остатков ОТСУТСТВУЕТ (остатки независимы)')

# === ГРАФИКИ ДЛЯ МНОЖЕСТВЕННОЙ РЕГРЕССИИ ===
fig2, axes2 = plt.subplots(2, 2, figsize=(15, 12))
fig2.suptitle('АНАЛИЗ ОСТАТКОВ: МНОЖЕСТВЕННАЯ РЕГРЕССИЯ (Weight = f(Position, Height, Age))', fontsize=16)

# 1. Гистограмма
axes2[0, 0].hist(residuals_multi, bins=15, density=True, alpha=0.7, color='lightcoral', edgecolor='black')
xmin, xmax = axes2[0, 0].get_xlim()
x = np.linspace(xmin, xmax, 100)
p = stats.norm.pdf(x, np.mean(residuals_multi), np.std(residuals_multi))
axes2[0, 0].plot(x, p, 'k', linewidth=2, label='Нормальное распределение')
axes2[0, 0].set_title('Гистограмма остатков')
axes2[0, 0].set_xlabel('Остатки')
axes2[0, 0].set_ylabel('Плотность')
axes2[0, 0].legend()
axes2[0, 0].grid(True, alpha=0.3)

# 2. QQ-plot
sm.qqplot(residuals_multi, line='s', ax=axes2[0, 1])
axes2[0, 1].set_title('QQ-график остатков')
axes2[0, 1].grid(True, alpha=0.3)

# 3. Остатки vs Предсказанные значения
axes2[1, 0].scatter(y_pred_multi, residuals_multi, alpha=0.7, color='red')
axes2[1, 0].axhline(y=0, color='red', linestyle='--')
axes2[1, 0].set_xlabel('Предсказанные значения')
axes2[1, 0].set_ylabel('Остатки')
axes2[1, 0].set_title('Остатки vs Предсказанные значения')
axes2[1, 0].grid(True, alpha=0.3)

# 4. Фактические vs Предсказанные значения
axes2[1, 1].scatter(y_multi, y_pred_multi, alpha=0.7, color='orange')
min_val = min(y_multi.min(), y_pred_multi.min())
max_val = max(y_multi.max(), y_pred_multi.max())
axes2[1, 1].plot([min_val, max_val], [min_val, max_val], 'r--', label='Идеальная линия')
axes2[1, 1].set_xlabel('Фактические значения')
axes2[1, 1].set_ylabel('Предсказанные значения')
axes2[1, 1].set_title('Фактические vs Предсказанные значения')
axes2[1, 1].legend()
axes2[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# === СРАВНЕНИЕ МОДЕЛЕЙ ===
print("\n" + "=" * 60)
print("СРАВНЕНИЕ МОДЕЛЕЙ")
print("=" * 60)

print(f"{'Параметр':<30} {'Парная':<15} {'Множественная':<15}")
print(f"{'-'*60}")
print(f"{'R²':<30} {r_squared_simple:<15.4f} {r_squared_multi:<15.4f}")
print(f"{'Скорректированный R²':<30} {adj_r_squared_simple:<15.4f} {adj_r_squared_multi:<15.4f}")
print(f"{'Стандартная ошибка':<30} {std_error_simple:<15.4f} {std_error_multi:<15.4f}")
print(f"{'Shapiro-Wilk p-value':<30} {shapiro_p_simple:<15.4f} {shapiro_p_multi:<15.4f}")
print(f"{'Kolmogorov-Smirnov p-value':<30} {ks_p_simple:<15.4f} {ks_p_multi:<15.4f}")
print(f"{'Durbin-Watson':<30} {dw_stat_simple:<15.4f} {dw_stat_multi:<15.4f}")

# Вывод о лучшей модели
print(f"\nВЫВОД:")
if adj_r_squared_multi > adj_r_squared_simple and std_error_multi < std_error_simple:
    print("→ МНОЖЕСТВЕННАЯ РЕГРЕССИЯ показывает лучшие результаты")
else:
    print("→ ПАРНАЯ РЕГРЕССИЯ показывает сравнимые или лучшие результаты")