Вы не можете выбрать более 25 тем
Темы должны начинаться с буквы или цифры, могут содержать дефисы(-) и должны содержать не более 35 символов.
105 строки
3.5 KiB
Python
105 строки
3.5 KiB
Python
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) |