Polynomial Regression 多项式回归模型

该文章作为机器学习的第五篇文章，主要介绍多项式回归模型的原理和实现方法。这个算法是建立在线性回归基础之上的，所以需要对线性回归有一定的了解。如果大家不是很了解前置知识请移步机器学习原理到Python代码实现之LinearRegression补番。

难度系数：⭐⭐⭐

更多相关工作请参考：Github

算法介绍

多项式回归是线性回归的扩展，它使用多项式函数来拟合数据。多项式回归的原理和线性回归类似，都是通过最小二乘法来拟合数据，只不过多项式回归的模型函数是多项式函数。当自变量与因变量之间的关系不是线性关系时，线性回归可能无法提供准确的预测。在这种情况下，可以使用多项式回归来处理非线性关系。通过增加自变量的高次方（例如平方项或立方项），多项式回归可以更好地捕获数据中的非线性变化。

算法原理解析

当我们想要通过数学模型来描述两个变量之间的关系时，有时候我们会发现，这两个变量之间的关系并不是一条直线，而是更接近于一条曲线。这时，线性回归就不再适用了。多项式回归算法就是用来处理这种情况的一种方法。

简单来说，多项式回归就是通过增加自变量的次数来让数学模型更好地拟合数据。比如说，如果原始的自变量是x，我们可以尝试把x变成x^2、x3等等，这样模型就能更好地描述数据的走势。

举个例子，假设我们有一组数据，自变量是x（比如气温），因变量是y（比如气压）。如果我们用线性回归来拟合这组数据，可能发现结果并不理想，因为气温和气压之间的关系并不是线性的。这时，我们就可以用多项式回归，把x变成x^2、x3等等，然后选择一个最佳的次数使得模型最能拟合数据。

对于初学者来说，理解多项式回归的原理可能有些复杂，但其实它就是一种通过改变自变量的形式来让模型更好地描述数据的方法。当你熟悉了这种方法后，你会发现它在很多场合都非常有用，比如在数据分析、机器学习等领域。

算法的数学原理和推导

在数学推导过程中，我们可以优先选择最小二乘法的算法求解，这是一个比较简单的思路。而我们之前已经学习过了线性回归的数学原理和推导，这里我们需要注意的是，多项式回归实际上是可以转换成线性回归任务来求解的。具体思路如下：

首先关于多项式回归，其公式如下：
$$y=w_0+w_{1}x+w_2{x}^2+...+w_n{x}^n$$

其中，$w_0,w_1,w_2,...,w_n$是模型的参数。这里我们可以假设$x1=x^1,x2=x^2,...,xn=x^n$，那么多项式回归可以转换成如下形式：
$$y=w_0+w_1{x}_1+w_2{x}_2+...+w_n{x}_n$$

在这样的前提下，问题便转换成的线性回归的问题，即：

$$w=({\mathbf{X}}^{\mathrm{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{T}}\mathbf{y}$$

而问题只到了这里吗？我们需要注意的是多项式回归中也会存在多元的问题，我们以二元多项式回归为例，其公式如下：
$$y=w_0+w_1{x}_1+w_2{x}_2+w_3{x}_1^2+w_4{x}_2^2+w_5{x}_1{x}_2$$

同样地，我们可以把$x_1^2,x_2^2,x_1{x}_2$这些变量转换成新的变量。但在多项式中，如何生成这些联合变量是需要注意的。那么如何生成N个变量的联合变量呢？简单来说就是将N次幂下的组会全部列举出来即可，也就是说3元3次幂的多项式，我们可以生成9个变量。

数据集介绍

数据集我们将继续采用波士顿房价数据集 。该数据集包含506个样本，13个特征，以及一个目标变量——房屋价格中位数。

波士顿房价数据集是一个非常经典的数据集，被广泛用于机器学习和数据分析领域。这个数据集包含了波士顿地区不同社区的房价信息:

参数	属性
CRIM–城镇人均犯罪率	城镇人均犯罪率
ZN - 占地面积超过25,000平方英尺的住宅用地比例。	住宅用地所占比例
INDUS - 每个城镇非零售业务的比例。	城镇中非商业用地占比例
CHAS - Charles River虚拟变量（如果是河道，则为1;否则为0	查尔斯河虚拟变量，用于回归分析
NOX - 一氧化氮浓度（每千万份）	环保指标
RM - 每间住宅的平均房间数	每栋住宅房间数
AGE - 1940年以前建造的自住单位比例	1940年以前建造的自住单位比例
DIS -波士顿的五个就业中心加权距离	与波士顿的五个就业中心加权距离
RAD - 径向高速公路的可达性指数	距离高速公路的便利指数
TAX - 每10,000美元的全额物业税率	每一万美元的不动产税率
PTRATIO - 城镇的学生与教师比例	城镇中教师学生比例
B - 1000（Bk - 0.63）^ 2其中Bk是城镇黑人的比例	城镇中黑人比例
LSTAT - 人口状况下降％	房东属于低等收入阶层比例
MEDV - 自有住房的中位数报价, 单位1000美元	自住房屋房价中位数

这个数据集的主要目的是通过机器学习算法，利用这14个特征预测房价中位数。在数据集中，每个样本包含一个社区的房价信息和相关的特征变量，例如社区的犯罪率、住宅用地比例、非商业用地比例、是否临河、房间数等。机器学习算法将根据这些特征变量预测房价中位数，从而帮助房地产经纪人、投资者或购房者更好地了解市场趋势和预测房价。

波士顿房价数据集是一个非常有价值的数据集，因为它包含了多个与房价相关的特征变量，并且数据来源于一个实际的房地产市场。这个数据集被广泛用于机器学习和数据分析的教学和实践，是入门机器学习和数据分析领域的经典案例之一。

接下来我们将对原始数据集进行处理，并对其进行特征工程，最终得到一个更加适合线性回归模型的数据集。数据集的地址在dataset\housing.data，大家可以直接使用。

代码实现

python"># 准备好我们需要使用的第三方包
import os
import numpy as np
import pandas as pandas
import matplotlib.pyplot as plt

关于数据处理部分和线性回归的相关内容已经在LinearRegression完成，大家可以链接直达查看，这里我们直接完成多项式回归的部分。

python">def load_data(file_path):
    # 读取数据文件
    names = ["CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE", "DIS", "RAD", "TAX", "PTRATIO", "B", "LSTAT", "MEDV"]
    data = pandas.read_csv(file_path, names=names, delim_whitespace=True)
    # 删除包含缺失值的数据行
    data = data.dropna()
    return data

def preprocess_data(data, func="del"):
    # 删除有缺失的数据
    if func == "del":
        data = data.dropna()
    # 通过均值的方式填充确实的数据
    elif func == "fill":
        data = data.fillna(data.mean())
    return data

# 将波士顿数据集按照8:2的比例划分成训练集和验证集
def split_data(data, test_ratio):
    np.random.seed(42)
    shuffled_indices = np.random.permutation(len(data))
    test_set_size = int(len(data) * test_ratio)
    test_indices = shuffled_indices[:test_set_size]
    train_indices = shuffled_indices[test_set_size:]
    return data.iloc[train_indices], data.iloc[test_indices]

# 划分训练集和验证集
data = load_data("dataset\\housing.data")
data = preprocess_data(data)                 # 该数据不存在缺失值
train_set, test_set = split_data(data, 0.2)

python"># 该部分是线性回归的代码
# 通过最小二乘法求解线性回归

class MyLinearRegression:
    def __init__(self):
        self.mean, self.std = None, None
        self.w, self.b = None, None
    
    def fit(self, X, y):
        X = self.data_preprocess(X)
        self.w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
        self.b = np.mean(y - X.dot(self.w))

    def data_preprocess(self, X):
        if self.mean is None:
            self.mean = np.mean(X, axis=0)
            self.std = np.std(X, axis=0)
        return (X - self.mean) / self.std
    
    def loss(self, y, y_pred):
        return np.mean((y - y_pred) ** 2)
    
    def predict(self, X):
        X = (X - self.mean) / self.std
        return np.dot(X, self.w) + self.b

为了实现多项式回归，我们在该代码的基础上补充一个特征扩展的函数polynomial_features。

python">class MyPolynomialRegression:
    def __init__(self, degree):
        self.mean, self.std = None, None
        self.w, self.b = None, None
        self.degree = degree

    def polynomial_features(self, X):
        n_samples, n_features = X.shape
        X_poly = np.zeros((n_samples, (self.degree + 1) * n_features))

        for i in range(n_samples):
            for j in range(n_features):
                X_poly[i, j] = X[i, j]
                for d in range(1, self.degree + 1):
                    X_poly[i, j + d * n_features] = X[i, j] ** d
        return X_poly

    def fit(self, X, y):
        X = self.data_preprocess(X)
        # 换成广义逆矩阵
        self.w = np.linalg.pinv(X.T.dot(X)).dot(X.T).dot(y)
        self.b = np.mean(y - X.dot(self.w))

    def data_preprocess(self, X):
        if self.mean is None:
            self.mean = np.mean(X, axis=0)
            self.std = np.std(X, axis=0)
        return (X - self.mean) / self.std
    
    def loss(self, y, y_pred):
        return np.mean((y - y_pred) ** 2)
    
    def predict(self, X):
        X = (X - self.mean) / self.std
        return np.dot(X, self.w) + self.b

python">def main(train_set, test_set):
    X_train = train_set.drop("MEDV", axis=1).to_numpy()
    y_train = train_set["MEDV"].to_numpy()
    X_test = test_set.drop("MEDV", axis=1).to_numpy()
    y_test = test_set["MEDV"].to_numpy()
    
    model = MyPolynomialRegression(1)
    new_X_train = model.polynomial_features(X_train)
    model.fit(new_X_train, y_train)
    w, b = model.w, model.b
    new_X_test = model.polynomial_features(X_test)
    y_pred = model.predict(new_X_test)
    mse = model.loss(y_test, y_pred)

    print("w:%s b:%7.5f 均方误差：%7.5f" % (w, b, mse))

main(train_set, test_set)

w:[-0.50053464  0.35526556  0.13197518  0.35921674 -0.99961199  1.56969408
 -0.08369355 -1.53935366  1.13469792 -0.89562242 -1.01521842  0.56380198
 -1.80394376 -0.50053464  0.35526556  0.13197518  0.35921674 -0.99961199
  1.56969408 -0.08369355 -1.53935366  1.13469792 -0.89562242 -1.01521842
  0.56380198 -1.80394376] b:22.79309 均方误差：24.39683

python">def main(train_set, test_set):
    X_train = train_set.drop("MEDV", axis=1).to_numpy()
    y_train = train_set["MEDV"].to_numpy()
    X_test = test_set.drop("MEDV", axis=1).to_numpy()
    y_test = test_set["MEDV"].to_numpy()
    
    model = MyPolynomialRegression(2)
    new_X_train = model.polynomial_features(X_train)
    model.fit(new_X_train, y_train)
    w, b = model.w, model.b
    new_X_test = model.polynomial_features(X_test)
    y_pred = model.predict(new_X_test)
    mse = model.loss(y_test, y_pred)

    print("w:%s b:%7.5f 均方误差：%7.5f" % (w, b, mse))

main(train_set, test_set)

w:[-1.516094   -0.5809663  -0.73752631  0.22598272 -1.10093148 -6.67993547
 -0.51577348 -2.72712851  2.48993965 -2.54340527 -5.46159781  1.09655201
 -4.70326324 -1.516094   -0.5809663  -0.73752631  0.22598272 -1.10093148
 -6.67993547 -0.51577348 -2.72712851  2.48993965 -2.54340527 -5.46159781
  1.09655201 -4.70326324  1.47804288  0.95756322  1.53839073  0.22598272
 -0.34211681 15.85790143  0.9004084   2.79459599 -2.4950686   3.77310186
  9.17541831 -1.54631387  5.17081091] b:22.79309 均方误差：13.56476

python">def main(train_set, test_set):
    X_train = train_set.drop("MEDV", axis=1).to_numpy()
    y_train = train_set["MEDV"].to_numpy()
    X_test = test_set.drop("MEDV", axis=1).to_numpy()
    y_test = test_set["MEDV"].to_numpy()
    
    model = MyPolynomialRegression(3)
    new_X_train = model.polynomial_features(X_train)
    model.fit(new_X_train, y_train)
    w, b = model.w, model.b
    new_X_test = model.polynomial_features(X_test)
    y_pred = model.predict(new_X_test)
    mse = model.loss(y_test, y_pred)

    print("w:%s b:%7.5f 均方误差：%7.5f" % (w, b, mse))

main(train_set, test_set)

w:[ -3.12181614   0.69113051  -1.80279254   0.15202838  20.12708751
 -18.44771319   0.66207189  -5.25538413   5.22454548 -16.43490115
 -10.61926454  -0.12098772  -7.4680723   -3.12181614   0.69113051
  -1.80279254   0.15202838  20.12708751 -18.44771319   0.66207189
  -5.25538413   5.22454548 -16.43490115 -10.61926454  -0.12098772
  -7.4680723    7.95556202  -5.22530726   8.94110781   0.15202838
 -82.97740985  64.10630895  -3.94501293  13.05204952 -26.03785295
  67.4394229   30.27624233   5.38400794  18.48837252  -4.08551448
   3.98928579  -5.460344     0.15202838  40.67452607 -24.86669419
   2.30035344  -5.33391462  17.36990835 -34.68778302 -10.66518554
  -4.50598064  -7.93109814] b:22.79309 均方误差：13.44853

python">def main(train_set, test_set):
    X_train = train_set.drop("MEDV", axis=1).to_numpy()
    y_train = train_set["MEDV"].to_numpy()
    X_test = test_set.drop("MEDV", axis=1).to_numpy()
    y_test = test_set["MEDV"].to_numpy()
    
    model = MyPolynomialRegression(4)
    new_X_train = model.polynomial_features(X_train)
    model.fit(new_X_train, y_train)
    w, b = model.w, model.b
    new_X_test = model.polynomial_features(X_test)
    y_pred = model.predict(new_X_test)
    mse = model.loss(y_test, y_pred)

    print("w:%s b:%7.5f 均方误差：%7.5f" % (w, b, mse))

main(train_set, test_set)

w:[-3.13576217e+00  1.33795475e+00 -3.86671528e+00  1.10910762e-01
 -5.57027945e+01  1.96576133e+02  3.36709006e+00 -8.92384650e+00
  2.99095537e+01 -8.26461108e+01  1.36463387e+02 -7.90368466e-01
 -1.19095303e+01 -3.13575390e+00  1.33794116e+00 -3.86672215e+00
  1.10916585e-01 -5.57027901e+01  1.96576132e+02  3.36709862e+00
 -8.92382988e+00  2.99095545e+01 -8.26461126e+01  1.36463391e+02
 -7.90369213e-01 -1.19095310e+01  6.92491190e+00 -8.10996305e+00
  2.65727235e+01  1.10916141e-01  3.65061430e+02 -1.30470777e+03
 -2.84979051e+01  3.37356419e+01 -4.33856483e+02  5.94089716e+02
 -8.65133452e+02  1.61789535e+01  5.67380632e+01 -1.47426425e+00
  6.67566748e+00 -2.92213051e+01  1.10916144e-01 -4.06259928e+02
  1.43897908e+03  3.65053918e+01 -2.62669788e+01  1.02734907e+03
 -7.36140307e+02  9.06662245e+02 -2.47580733e+01 -6.19611741e+01
 -1.71383536e+00 -1.28440661e+00  9.96183919e+00  1.10916144e-01
  1.51032515e+02 -5.24888329e+02 -1.53365285e+01  7.60764390e+00
 -6.48326371e+02  3.03847059e+02 -3.15819847e+02  1.07407988e+01
  2.47860487e+01] b:22.79309 均方误差：12.90421

同样我们给出Sklearn的实现方式

python">from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression


def sklearn_main(train_set, test_set):
    X_train = train_set.drop("MEDV", axis=1).to_numpy()
    y_train = train_set["MEDV"].to_numpy()
    X_test = test_set.drop("MEDV", axis=1).to_numpy()
    y_test = test_set["MEDV"].to_numpy()
    poly = PolynomialFeatures(degree=2)
    X_train_poly = poly.fit_transform(X_train)
    X_test_poly = poly.transform(X_test)
    model = LinearRegression()
    model.fit(X_train_poly, y_train)
    y_pred = model.predict(X_test_poly)
    mse = np.mean((y_test - y_pred) ** 2)
    print("w:%s b:%7.5f 均方误差：%7.5f" % (model.coef_, model.intercept_, mse))

sklearn_main(train_set, test_set)

w:[ 5.50377850e+08 -7.25155954e+00  6.43056690e-01 -4.72286165e+00
  3.67693992e+01  2.85812138e+02  1.78769612e+01  6.98451634e-01
 -4.06639487e+00  2.63397019e+00  1.05224457e-02  7.65707397e+00
  1.63069648e-01 -2.01526636e-01  9.59578493e-05  7.81777220e-02
  6.29149427e-01  2.57872491e+00 -2.01595602e+00  2.17379223e-01
 -2.33140101e-03 -9.32989554e-02  2.60837521e-01 -3.47855999e-02
  6.19073750e-01 -4.81331723e-04  3.32912130e-02 -3.99060606e-04
 -5.51134937e-03 -3.54563187e-03 -1.35090986e+00 -3.15040131e-02
  1.41763713e-03 -1.34509936e-02 -2.04447839e-02  8.08488035e-04
  3.82531779e-04  3.64327939e-04 -1.14836963e-02  4.27863945e-02
 -5.34629140e-02 -2.04285326e-01  3.03266276e-01  3.31412539e-03
  1.65597030e-01 -7.52610044e-02  3.37154323e-04  3.31622015e-02
  2.23400004e-03 -1.55677424e-02  3.67693992e+01 -3.66681092e+01
 -5.26721223e+00 -3.86728128e-02 -9.06706235e-01  1.41939708e-01
 -4.25350723e-03 -8.56691474e-01  1.05400141e-02 -2.90967558e-01
 -9.88009749e+01 -8.83019108e+00  1.26581842e-01  1.48123555e+01
 -2.07155861e+00  2.66731755e-01 -1.45194402e+01 -1.67237935e-02
  1.00515923e+00  1.16590564e+00 -8.84064023e-02 -3.63283726e-01
 -2.23549284e-01 -1.19088508e-02 -5.10726196e-01 -5.24156067e-03
 -2.35689833e-02 -2.29435884e-04 -5.62095379e-04  1.38004970e-02
 -4.79761383e-04  4.07525858e-03 -4.61000170e-04 -9.78001291e-03
  4.35919277e-01  3.15534696e-03 -5.44676149e-03 -2.00588708e-01
 -4.32624813e-03  4.75137857e-02 -1.20917920e-01  8.57579379e-03
 -1.35700032e-01  3.42304727e-03 -2.32734725e-02 -4.68271257e-05
  5.11291048e-03 -3.20676019e-04 -1.28660473e-03  1.58204620e-02
  2.46381159e-03  2.31443378e-02 -3.53936035e-05 -1.57841718e-04
  1.67967389e-02] b:-550378077.60859 均方误差：14.46105

以上便是本次关于多项式回归的求解算法，希望大家可以有所收获。