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)) 

axes[0,0].scatter(x3, x4)
axes[0,0].set_title('Зависимость X4 от X3')

axes[0,1].scatter(x3, x4)
axes[0,1].set_title('Предсказанные значения')



# Парная регрессия
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(x3.reshape(-1,1))
axes[0,0].plot(x3, line_y1, color='red', linewidth=2)

line_x3_pred = np.linspace(min(x3), 2*max(x3), 150)
line_y1_pred = lr.predict(line_x3_pred.reshape(-1,1))
axes[0,1].plot(line_x3_pred, line_y1_pred, color='red', linewidth=2)

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

print(np.dot(x4 - line_y1, line_y1))

print('Парная регрессия. Коэффициент детерминации: ', r2_score(x4,line_y1))
print('Скорр. коэффициент детерминации: ', 1 - (1 - r2_score(x4,line_y1)*(149/147)))
print('Множ. коэффициент корреляции: ', np.sqrt(r2_score(x4,line_y1)))
print('Стандартная ошибка регрессии: ', np.sqrt(Qint/147))

# Множественная регрессия
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_
axes[0,0].plot(line_x3, line_y2, color='green', linewidth=2)

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_
axes[0,1].plot(line_x3_pred, line_y2_pred, color='green', linewidth=2)

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

print('Множ. регресиия. Коэффициент детерминации: ', r2_score(x4,line_y2))
print('Скорр. коэффициент детерминации: ', 1 - (1 - r2_score(x4,line_y2)*(149/147)))
print('Множ. коэффициент корреляции: ', np.sqrt(r2_score(x4,line_y2)))
print('Стандартная ошибка регрессии: ', np.sqrt(Qint/147))

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

axes[1,0].scatter([i for i in range(1,151)], residuals1)
axes[1,0].set_title('Парная регрессия. График остатков')

axes[1,1].scatter([i for i in range(1,151)], residuals2)
axes[1,1].set_title('Множ. регрессия. График остатков')


plt.show()