import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.linear_model import LinearRegression
from statsmodels.stats.stattools import durbin_watson
from scipy import stats
import matplotlib.pyplot as plt

# === Датасет Ирисы Фишера ===
iris = load_iris()
iris_pd = pd.DataFrame(data=np.c_[iris['data'], iris['target']],
                       columns=iris['feature_names'] + ['target'])

# === Признаки и целевая переменная ===
x1 = iris_pd['sepal length (cm)'].values
x2 = iris_pd['sepal width (cm)'].values
x3 = iris_pd['petal length (cm)'].values
x4 = iris_pd['petal width (cm)'].values  # целевая переменная

# === Парная регрессия: x3 -> x4 ===
lr = LinearRegression().fit(x3.reshape(-1,1), x4)
y_pred1 = lr.predict(x3.reshape(-1,1))
residuals1 = x4 - y_pred1

# === Множественная регрессия: x1,x2,x3 -> x4 ===
X = np.vstack((x1, x2, x3)).T
mr = LinearRegression().fit(X, x4)
y_pred2 = mr.predict(X)
residuals2 = x4 - y_pred2

# === Графики остатков ===
fig, axes = plt.subplots(1, 2, figsize=(12,5))
axes[0].hist(residuals1, bins=15, color='lightblue', edgecolor='black')
axes[0].set_title('Остатки парной регрессии')
axes[1].hist(residuals2, bins=15, color='lightgreen', edgecolor='black')
axes[1].set_title('Остатки множественной регрессии')
plt.show()

# === Проверка нормальности остатков (Shapiro-Wilk) ===
print("=== Проверка нормальности остатков ===")
for i, res in enumerate([residuals1, residuals2], start=1):
    stat, p = stats.shapiro(res)
    print(f"Регрессия {i}: Статистика Шапиро–Уилка = {stat:.4f}, p-value = {p:.4f}")
    if p > 0.05:
        print("→ Остатки распределены нормально (H0 не отвергается).")
    else:
        print("→ Остатки не распределены нормально (H0 отвергается).")

# === Статистика Дурбина–Уотсона ===
from statsmodels.stats.stattools import durbin_watson
print("\n=== Статистика Дурбина–Уотсона ===")
dw1 = durbin_watson(residuals1)
dw2 = durbin_watson(residuals2)
print(f"Парная регрессия: DW = {dw1:.4f}")
print(f"Множественная регрессия: DW = {dw2:.4f}")

# Интерпретация
def interpret_dw(dw):
    if dw < 1.5:
        return "→ Положительная автокорреляция остатков."
    elif dw > 2.5:
        return "→ Отрицательная автокорреляция остатков."
    else:
        return "→ Остатки независимы (автокорреляции нет)."

print("Парная регрессия:", interpret_dw(dw1))
print("Множественная регрессия:", interpret_dw(dw2))