import numpy as np
from sklearn.datasets import load_iris
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sb
import pingouin as pg
from scipy import stats

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

# Формирование выборок по признакам
x1 = np.array(iris_pd['sepal length (cm)']) # первый признак - длина чашелистика
x2 = np.array(iris_pd['sepal width (cm)']) # второй признак - ширина чашелистика
x3 = np.array(iris_pd['petal length (cm)']) # третий признак - длина лепестка
x4 = np.array(iris_pd['petal width (cm)']) # четвертый признак - ширина лепестка

# Пусть X4 - выходные переменные, а X3 - входные переменные, так как они коррелируют сильнее остальных комбинаций
fig, axes = plt.subplots(2, 2, figsize=(10, 10)) 


# Парная регрессия
N = 150
K = 2
lr = LinearRegression().fit(x3.reshape(-1,1), x4)

line_x3 = np.linspace(min(x3), max(x3),150)
line_y1 = lr.predict(line_x3.reshape(-1,1))

line_x3_pred = np.linspace(min(x3), 2*max(x3), 150)
line_y1_pred = lr.predict(line_x3_pred.reshape(-1,1))

meanY = np.mean(x4)
Qreg = np.sum((line_y1 - meanY)**2)
Qtotal = np.sum((x4 - meanY)**2)
Qint = np.sum((x4 - line_y1)**2)
Sreg = (1/(K-1))*Qreg
Stotal = (1/(N-1))*Qtotal
Sint = (1/(N-K))*Qint


# Множественная регрессия
N = 150
K = 4
X = np.vstack((x1, x2, x3)).T
mr = LinearRegression().fit(X, x4)

line_x1 = np.linspace(min(x1), max(x1), 150)
line_x2 = np.linspace(min(x2), max(x2), 150)
line_x3 = np.linspace(min(x3), max(x3), 150)

line_y2 = mr.coef_[0] * line_x1 + mr.coef_[1] * line_x2 + mr.coef_[2] * line_x3 + mr.intercept_

line_x1_pred = np.linspace(min(x1), 2*max(x1), 150)
line_x2_pred = np.linspace(min(x2), 2*max(x2), 150)
line_x3_pred = np.linspace(min(x3), 2*max(x3), 150)
line_y2_pred = mr.coef_[0] * line_x1_pred + mr.coef_[1] * line_x2_pred + mr.coef_[2] * line_x3_pred + mr.intercept_

Qreg = np.sum((line_y2 - meanY)**2)
Qint = np.sum((x4 - line_y2)**2)
Sreg = (1/(K-1))*Qreg
Stotal = (1/(N-1))*Qtotal
Sint = (1/(N-K))*Qint


# График остатков
residuals1 = x4 - line_y1
residuals2 = x4 - line_y2

axes[0,0].hist(residuals1)
axes[0,1].hist(residuals2)


plt.show()

# Проверка остатков на нормальность по критерию К-С
print('Тест Колмогорова-Смирнова. Уровень значимости a = 0.05')

for x in [residuals1,residuals2]:
    x_mean, std = np.mean(x), np.std(x, ddof = 1)
    x_std = (x - x_mean)/std
    pvalue = stats.kstest(x_std, "norm", alternative='less').pvalue
    if pvalue > 0.05: result = 'H0 не должна быть отвергнута'
    else: result = 'H0 должна быть отвергнута'
    print('p-value: {}, так что {}'.format(pvalue, result))
    
# Статистика Дурбина-Уотсона
# для парной регрессии
a1 = 0
gl = 1.706
gu = 1.760
for i in range(N-1): a1 += (residuals1[i] - residuals1[i+1])**2
gamma = a1/sum(residuals1**2)
print(gamma, gl, gu)

# для множ. регрессии
a2 = 0
gl = 1.679
gu = 1.788
for i in range(N-1): a2 += (residuals2[i] - residuals2[i+1])**2
gamma = a2/sum(residuals2**2)
print(gamma, gl, gu)