Множественный регрессионный Python

Question

Я действительно не понимаю, что не так с моим (простым) кодом ...

Я просто хочу проверить множественную линейную регрессию (....!).

 import pandas as pd
 import numpy as np
 import scipy.stats as st
 import sklearn
 n = 1000

 X1 = linspace(2, 8.5, n)
 X2 = linspace(-4, 2.9, n)
 X3 = linspace(-1, 16, n)

 X = np.transpose( [X1, X2, X3] )

 Y = 2*X1 + 3.2*X2 -1.2*X3  +  4  + st.norm.rvs(size = n, loc = 0, scale = 0.6)

 X = pd.DataFrame( X , columns = ["X1", "X2", "X3"])
 Y = pd.DataFrame(Y, columns = ["Y"])

 #Create linear regression object:
 my_reg = sklearn.linear_model.LinearRegression()

 #Train:
 my_reg.fit(X, Y)

 print('Coefficients: \n', my_reg.coef_)
 print('Constant: \n', my_reg.intercept_)

И я получаю некоторые глупые результаты, например, коэффициенты [0.25127347 0.26673645 0.65717676] ...

Я также попробовал способ OLS, но я все еще получаю бессмысленный коф (немного другой, но все еще глупый)

(Это работа с линейной регрессией с одной переменной, что-то вроде Y = 2 * X + 5, я бы получил коэф и перехват, очень близкий к реальному)

Answer 1

Тот факт, что коэффициенты совсем не похожи на «истинные», которые вы установили, указывает на то, что мультиколлинеарность может быть проблемой. Проблема с вашим кодом в том, что ваша матрица X почти единственная, что делает числовые результаты нестабильными. Как видно из графиков @ R.yan, ваши X1 и X2 практически идентичны, за исключением линейного сдвига. Это подтверждается тем фактом, что ваша матрица X, которая имеет 1000 строк и три столбца, имеет ранг только 2. См .:

np.linalg.matrix_rank(X)
Out[26]: 2

Попробуйте вместо этого:

import pandas as pd
import numpy as np
import scipy.stats as st
import sklearn
from sklearn.linear_model import LinearRegression
n = 1000

# adding noise to your data:
X1 = np.linspace(2, 8.5, n) + st.norm.rvs(size=n ,loc = 0, scale = 1)
X2 = np.linspace(-4, 2.9, n) + st.norm.rvs(size=n ,loc = 0, scale = 1)
X3 = np.linspace(-1, 16, n) + st.norm.rvs(size=n ,loc = 0, scale = 1)

X = np.transpose( [X1, X2, X3] )

Y = 2*X1 + 3.2*X2 -1.2*X3  +  4  + st.norm.rvs(size=1000 ,loc = 0, scale = 1)

X = pd.DataFrame( X , columns = ["X1", "X2", "X3"])
Y = pd.DataFrame(Y, columns = ["Y"])
#Create linear regression object:
my_reg = sklearn.linear_model.LinearRegression(fit_intercept = True)

#Train:
res = my_reg.fit(X, Y)

print('Coefficients: \n', my_reg.coef_)
print('Constant: \n', my_reg.intercept_)

Coefficients: 
 [[ 1.99273588  3.20068392 -1.19688422]]
Constant: 
 [ 4.02296003]

Теперь мы получаем правильные коэффициенты и матрицу полного ранга:

np.linalg.matrix_rank(X)
Out[32]: 3

Обратите внимание, что в линейной регрессии X должен иметь ранг, равный количеству столбцов (или строк, если оно меньше). Если это не так, это означает, что существует мультиколлинеарность, которая отображает числовые результаты для инверсии X'X нестабильной (в зависимости от того, какой алгоритм используется). См. это описание для получения дополнительной информации о мультиколлинеарности.

Answer 2

Спасибо всем!

Я некоторое время не выполнял линейную регрессию, и, конечно, это потому, что X не обратим (в R это дает мне 'nan').

Так что это не был умный вопрос ...

Еще раз спасибо!

Answer 3

Полагаю, код дает правильный ответ. Я строю предсказанную базу Y на coef_ и intercept_ из вашей регрессии и получаю следующий график.

import pandas as pd
import numpy as np
import scipy.stats as st
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt
n = 1000

X1 = np.linspace(2, 8.5, n)
X2 = np.linspace(-4, 2.9, n)
X3 = np.linspace(-1, 16, n)

X = np.transpose( [X1, X2, X3] )

Y = 2*X1 + 3.2*X2 -1.2*X3  +  4  + st.norm.rvs(size=1000 ,loc = 0, scale = 0.6)

X = pd.DataFrame( X , columns = ["X1", "X2", "X3"])
Y = pd.DataFrame(Y, columns = ["Y"])
#Create linear regression object:
my_reg = sklearn.linear_model.LinearRegression()
plt.plot(Y, color='blue', label='Y')
#Train:
res = my_reg.fit(X, Y)

print('Coefficients: \n', my_reg.coef_)
print('Constant: \n', my_reg.intercept_)

plt.scatter(X.index.values,X['X1'], c='black')
plt.scatter(X.index.values,X['X2'], c='black')
plt.scatter(X.index.values,X['X3'], c='black')


Y_pred = my_reg.coef_[0][0]*X['X1'] + my_reg.coef_[0][1]*X['X2'] +my_reg.coef_[0][2]*X['X3'] + my_reg.intercept_
plt.plot(Y_pred, color="red", label='predict')
plt.legend()

Out[]: ('Coefficients: \n', array([[  3.13842691e+12,   1.01316187e+13,  -5.31223199e+12]]))
('Constant: \n', array([  2.89373889e+13]))