Python: Разработка модели множественной линейной регрессии с нуля

Question

Я пытаюсь создать модель множественной линейной регрессии с нуля в python. Используемый набор данных: Boston Housing Dataset из Sklearn. Поскольку мое внимание было сосредоточено на построении модели, я не выполнял никаких шагов предварительной обработки данных. Однако я использовал модель OLS для вычисления p-значений и удалил 3 функции из данных. После этого я использовал модель линейной регрессии, чтобы узнать веса для каждой функции.

import pandas as pd
from sklearn.datasets import load_boston
from sklearn.linear_model import LinearRegression

X=load_boston()
data=pd.DataFrame(X.data,columns=X.feature_names)
y=X.target
data.head()

#dropping three features
data=data.drop(['INDUS','NOX','AGE'],axis=1)

#new shape of the data (506,10) not including the target variable

#Passed the whole dataset to Linear Regression Model
model_lr=LinearRegression()
model_lr.fit(data,y)

model_lr.score(data,y)
0.7278959820021539

model_lr.intercept_
22.60536462807957   #----- intercept value

model_lr.coef_
array([-0.09649731,  0.05281081,  2.3802989 ,  3.94059598, -1.05476566,
        0.28259531, -0.01572265, -0.75651996,  0.01023922, -0.57069861]) #--- coefficients

Теперь я хотел вычислить коэффициенты вручную в Excel перед созданием модели в python. Чтобы вычислить вес каждой функции, я использовал эту формулу:

Расчет весов элементов

Для вычисления точки пересечения я использовал формулу b0 = mean (y) -b1 * mean (x1) -b2 * (mean (x2) ....- bn * mean (xn)

Значение точки пересечения по моим расчетам было 22,63551387 (почти то же, что и в модели)

Проблема в том, что веса характеристик из моих вычислений далеки от веса линейной модели sklearn.

-0.002528644 #-- CRIM
-0.001028914 #-- Zn
-0.038663314 #-- CHAS
-0.035026972 #-- RM
-0.014275311 #-- DIS
-0.004058291 #-- RAD
-0.000241103 #-- TAX
-0.015035534 #-- PTRATIO
-0.000318376 #-- B
-0.006411897 #-- LSTAT

Используя первую строку в качестве тестовых данных для проверки моих расчетов, Я получаю 22,73167044199992, в то время как модель линейной регрессии предсказывает 30,42657776. Исходное значение - 24.

Но как только я проверяю другие строки, модель sklearn имеет больше вариаций, в то время как прогнозы, сделанные на основе весов из моих вычислений, все показывают значения, близкие к 22.

Я думаю, что делаю ошибку при вычислении весов, но не уверен, в чем проблема? Есть ли ошибка в моих расчетах? Почему все мои коэффициенты из расчетов так близки к 0?

Вот мой код для расчета коэффициентов: (здесь новичок)

 x_1=[]
 x_2=[]
 for i,j in zip(data['CRIM'],y):
      mean_x=data['CRIM'].mean()
      mean_y=np.mean(y)
      c=i-mean_x*(j-mean_y)
      d=(i-mean_x)**2
      x_1.append(c)
      x_2.append(d)
 print(sum(x_1)/sum(x_2))

Спасибо, что прочитали этот длинный пост, я ценю это.

Answer 1

Похоже, проблема в расчете коэффициента. Приведенная вами формула для вычисления коэффициентов представлена ​​в скалярной форме и используется для простейшего случая линейной регрессии, а именно с одним признаком x.

EDIT

Теперь, увидев ваш код для расчета коэффициента, проблема прояснилась. Вы не можете использовать это уравнение для вычисления коэффициентов каждой характеристики независимо друг от друга, поскольку каждый коэффициент будет зависеть от всех функций. Я предлагаю вам взглянуть на вывод решения этой задачи оптимизации методом наименьших квадратов в простом случае здесь и в общем случае здесь . И в качестве общего совета по возможности придерживайтесь матричной реализации, поскольку это значительно эффективнее.

Однако в этом случае у нас есть 10-мерный вектор признаков, и поэтому в матричной записи он становится.

См. Вывод здесь

Я подозреваю, что вы допустили некоторую вычислительную ошибку здесь, поскольку реализовали это в python с использованием скалярной формулы более утомительна и неопрятна, чем матричный эквивалент. Но так как вы не делились этим миром своего кода, это трудно понять.

Вот пример того, как вы бы это реализовали:

def calc_coefficients(X,Y):
    X=np.mat(X)
    Y = np.mat(Y)
    return np.dot((np.dot(np.transpose(X),X))**(-1),np.transpose(np.dot(Y,X)))

def score_r2(y_pred,y_true):
    ss_tot=np.power(y_true-y_true.mean(),2).sum()
    ss_res = np.power(y_true -y_pred,2).sum()
    return 1 -ss_res/ss_tot


X = np.ones(shape=(506,11))
X[:,1:] = data.values

B=calc_coefficients(X,y)
##### Coeffcients 
B[:]
matrix([[ 2.26053646e+01],
        [-9.64973063e-02],
        [ 5.28108077e-02],
        [ 2.38029890e+00],
        [ 3.94059598e+00],
        [-1.05476566e+00],
        [ 2.82595310e-01],
        [-1.57226536e-02],
        [-7.56519964e-01],
        [ 1.02392192e-02],
        [-5.70698610e-01]])

#### Intercept 
B[0]
matrix([[22.60536463]])

y_pred = np.dot(np.transpose(B),np.transpose(X))
##### First 5 rows predicted
np.array(y_pred)[0][:5]

array([30.42657776, 24.80818347, 30.69339701, 29.35761397, 28.6004966 ])

##### First 5 rows Ground Truth
y[:5]

array([24. , 21.6, 34.7, 33.4, 36.2])

### R^2 score
score_r2(y_pred,y)

0.7278959820021539