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

现在我想在用 python 创建模型之前在 excel 中手动计算系数。为了计算每个特征的权重，我使用了这个公式：

Calculating the Weights of the Features

为了计算截距，我使用了公式 b0 = 均值(y)-b1*均值(x1)-b2*(均值(x2)....-bn*均值(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。

编辑

现在看到你的系数计算代码后，问题就更清楚了。 您不能使用这个方程来计算每个特征的系数，因为每个系数都取决于所有的特征。我建议你看看这个最小二乘优化问题在简单情况here 和一般情况here 中的解的推导。作为一般提示，只要有可能就坚持使用矩阵实现，因为这从根本上更有效。

但是，在这种情况下，我们有一个 10 维的特征向量，因此在矩阵表示法中它变成了。

见推导here

我怀疑你在这里犯了一些计算错误，因为在 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

Answer 2

完整的解决方案 - 2020 - 波士顿数据集

正如另一位所说，要计算线性回归的系数，您必须计算

β = (X^T X)^-1 X^T y

这会给你系数（特征的所有 B + 截距）。

一定要在 X 中添加一个全为 1 的列来计算截距（更多代码见）

Main.py

from sklearn.datasets import load_boston
import numpy as np

from CustomLibrary import CustomLinearRegression
from CustomLibrary import CustomMeanSquaredError

boston = load_boston()
X = np.array(boston.data, dtype="f")
Y = np.array(boston.target, dtype="f")

regression = CustomLinearRegression()
regression.fit(X, Y)
print("Projection matrix sk:", regression.coefficients, "\n")
print("bias sk:", regression.intercept, "\n")

Y_pred = regression.predict(X)
loss_sk = CustomMeanSquaredError(Y, Y_pred)
print("Model performance:")
print("--------------------------------------")
print("MSE is {}".format(loss_sk))
print("\n")

CustomLibrary.py

import numpy as np

class CustomLinearRegression():
    def __init__(self):
        self.coefficients = None
        self.intercept = None

    def fit(self, x , y):
        x = self.add_one_column(x)

        x_T = np.transpose(x)
        inverse = np.linalg.inv(np.dot(x_T, x))
        pseudo_inverse = inverse.dot(x_T)
        coef = pseudo_inverse.dot(y)

        self.intercept = coef[0]
        self.coefficients = coef[1:]

        return coef

    def add_one_column(self, x):
        '''
        the fit method with x feature return x coefficients ( include the intercept)
        so for have the intercept + x feature coefficients we have to add one column ( in the beginning )
        with all 1ones
        '''

        X = np.ones(shape=(x.shape[0], x.shape[1] +1))
        X[:, 1:] = x
        return X


    def predict(self, x):
        predicted = np.array([])
        for sample in x:
            result = self.intercept
            for idx, feature_value_in_sample in enumerate(sample):
                result += feature_value_in_sample * self.coefficients[idx]
            predicted = np.append(predicted, result)

        return predicted


def CustomMeanSquaredError(Y, Y_pred):
    mse = 0
    for idx,data in enumerate(Y):
        mse += (data - Y_pred[idx])**2

    return  mse * (1 / len(Y))