实现多元线性回归

多元线性回归的实现
下面我们来使用python代码实现多元线性回归：

import numpy as np
from .metrics import r2_score

class LinearRegression:

    def __init__(self):
        """初始化Linear Regression模型"""
        self.coef_ = None    # 系数（theta0~1 向量）
        self.interception_ = None   # 截距（theta0 数）
        self._theta = None  # 整体计算出的向量theta

    def fit_normal(self, X_train, y_train):
        """根据训练数据X_train，y_train训练Linear Regression模型"""
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"
        # 正规化方程求解
        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)

        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]
        return self

    def predict(self, X_predict):
        """给定待预测的数据集X_predict，返回表示X_predict的结果向量"""
        assert self.interception_ is not None and self.coef_ is not None, \
            "must fit before predict"
        assert X_predict.shape[1] == len(self.coef_), \
            "the feature number of X_predict must be equal to X_train"
        X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
        y_predict = X_b.dot(self._theta)
        return y_predict

    def score(self, X_test, y_test):
        """很倔测试机X_test和y_test确定当前模型的准确率"""
        y_predict = self.predict(self, X_test)
        return r2_score(y_test, y_predict)
    

    def __repr__(self):
        return "LinearRegression()"

其实在代码中，思想很简单，就是使用公式即可。其中有一些知识点：
1、np.hstack(tup)：参数tup可以是元组，列表，或者numpy数组，返回结果为numpy的数组。按列顺序把数组给堆叠起来（加一个新列）。
2、np.ones()：返回一个全1的n维数组，有三个参数：shape（用来指定返回数组的大小）、dtype（数组元素的类型）、order（是否以内存中的C或Fortran连续（行或列）顺序存储多维数据）。后两个参数都是可选的，一般只需设定第一个参数。（类似的还有np.zeros()返回一个全0数组）
3、numpy.linalg模块包含线性代数的函数。使用这个模块，可以计算逆矩阵、求特征值、解线性方程组以及求解行列式等。inv函数计算逆矩阵
4、T：array的方法，对矩阵进行转置。
5、dot：点乘

调用
下面我们可以在jupyter notebook中调用我们的算法：

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

boston = datasets.load_boston()

X = boston.data
y = boston.target

X = X[y<50.0]
y = y[y<50.0]

X.shape
输出：(490, 13)

y.shape
输出：(490, )

from myAlgorithm.model_selection import train_test_split
from myAlgorithm.LinearRegression import LinearRegression

X_train, X_test, y_train, y_test = train_test_split(X, y, seed = 666)

reg = LinearRegression()
reg.fit_normal(X_train, y_train)

reg.coef_
输出：
array([-1.18919477e-01,  3.63991462e-02, -3.56494193e-02,  5.66737830e-02,
       -1.16195486e+01,  3.42022185e+00, -2.31470282e-02, -1.19509560e+00,
        2.59339091e-01, -1.40112724e-02, -8.36521175e-01,  7.92283639e-03,
       -3.81966137e-01])

reg.interception_
输出：
34.16143549622471

reg.score(X_test, y_test)
输出：
0.81298026026584658

我们看到，reg.coef_这一项的结果是13个系数，这13个系数有正有负。正负代表着该系数所乘的特征与预测目标是正相关还是负相关。正相关，特征越大房价越高；负相关，特征越大，房价越低。而系数绝对值的大小决定了影响程度。
下面我们对所有的系数按照数值由小到大进行排序：

np.argsort(reg.coef_)
输出：
array([ 4,  7, 10, 12,  0,  2,  6,  9, 11,  1,  3,  8,  5])

将这个返回结果作为索引，返回排序后索引所对应的特征名：

boston.feature_names[np.argsort(reg.coef_)]
输出：
array(['NOX', 'DIS', 'PTRATIO', 'LSTAT', 'CRIM', 'INDUS', 'AGE', 'TAX',
       'B', 'ZN', 'CHAS', 'RAD', 'RM'], dtype='<U7')

这也说明了线性回归算法，具有可解释性。