多项式回归实现与应用

5556 字 · 620 阅读 · 2023 年 05 月 03 日

本文已更新,你可以访问 AI By Doing 以获得更好的阅读体验。
本篇文章需 特别授权许可,内容版权归作者所有,未经授权,禁止转载。

介绍

前面的文章中,相信你已经对线性回归有了充分的了解。掌握一元和多元线性回归之后,我们就能针对一些有线性分布趋势的数据进行回归预测。但是,生活中还常常会遇到一些分布不那么「线性」的数据,例如像股市的波动、交通流量等。那么对于这类非线性分布的数据,就需要通过这篇文章介绍的方法来处理。

知识点

  • 多项式回归介绍
  • 多项式回归基础
  • 多项式回归预测

多项式回归介绍

在线性回归中,我们通过建立自变量 $x$ 的一次方程来拟合数据。而非线性回归中,则需要建立因变量和自变量之间的非线性关系。从直观上讲,也就是拟合的直线变成了「曲线」。

对于非线性回归问题而言,最简单也是最常见的方法就是这篇文章要讲解的「多项式回归」。多项式是中学时期就会接触到的概念,这里引用 维基百科 的定义如下:

多项式(Polynomial)是代数学中的基础概念,是由称为未知数的变量和称为系数的常量通过有限次加法、加减法、乘法以及自然数幂次的乘方运算得到的代数表达式。多项式是整式的一种。未知数只有一个的多项式称为一元多项式;例如 $x^2-3x+4$ 就是一个一元多项式。未知数不止一个的多项式称为多元多项式,例如 $x^3-2xyz^2+2yz+1$ 就是一个三元多项式。

多项式回归基础

首先,我们通过一组示例数据来认识多项式回归

# 加载示例数据
x = [4, 8, 12, 25, 32, 43, 58, 63, 69, 79]
y = [20, 33, 50, 56, 42, 31, 33, 46, 65, 75]

示例数据一共有 10 组,分别对应着横坐标和纵坐标。接下来,通过 Matplotlib 绘制数据,查看其变化趋势。

from matplotlib import pyplot as plt
%matplotlib inline

plt.scatter(x, y)

实现 2 次多项式拟合

接下来,通过多项式来拟合上面的散点数据。首先,一个标准的一元高阶多项式函数如下所示:

$$ y(x, w) = w*0 + w_1x + w_2x^2 +...+w_mx^m = \sum\limits*{j=0}^{m}w_jx^j \tag{1}$$

其中,$m$ 表示多项式的阶数,$x^j$ 表示 $x$ 的 $j$ 次幂,$w$ 则代表该多项式的系数。

当我们使用上面的多项式去拟合散点时,需要确定两个要素,分别是:多项式系数 $w$ 以及多项式阶数 $m$,这也是多项式的两个基本要素。

如果通过手动指定多项式阶数 $m$ 的大小,那么就只需要确定多项式系数 $w$ 的值是多少。例如,这里首先指定 $m=2$,多项式就变成了:

$$ y(x, w) = w*0 + w_1x + w_2x^2= \sum\limits*{j=0}^{2}w_jx^j \tag{2}$$

当我们确定 $w$ 的值的大小时,就回到了前面线性回归中学习到的内容。

首先,我们构造两个函数,分别是用于拟合的多项式函数,以及误差函数。

def func(p, x):
    # 根据公式,定义 2 次多项式函数
    w0, w1, w2 = p
    f = w0 + w1*x + w2*x*x
    return f


def err_func(p, x, y):
    # 残差函数(观测值与拟合值之间的差距)
    ret = func(p, x) - y
    return ret

接下来,就是使用最小二乘法求解最优参数的过程。这里为了方便,我们直接使用 Scipy 提供的最小二乘法类,得到最佳拟合参数。当然,你完全可以按照线性回归文章中最小二乘法公式自行求解参数。不过,实际工作中为了快速实现,往往会使用像 Scipy 这样现成的函数,这里也是为了给大家多介绍一种方法。

import numpy as np
from scipy.optimize import leastsq

p_init = np.random.randn(3)  # 生成 3 个随机数
# 使用 Scipy 提供的最小二乘法函数得到最佳拟合参数
parameters = leastsq(err_func, p_init, args=(np.array(x), np.array(y)))

print('Fitting Parameters: ', parameters[0])

小贴士

关于 scipy.optimize.leastsq() 的具体使用介绍,可以阅读 ↗ 官方文档。特别注意的是,上面我们定义的 p_init 并不是官方文档中最小化初始参数的意思,因为最小二乘法是解析解,其不会涉及到从某个参数开始迭代的过程。实际上,这里 p_init 的具体取值不会影响求解结果,所以我们使用了随机值,但是其值的个数决定了最终多项式的次数。具体来说,$n$ 个值则最终求解出的是 $n-1$ 次多项式。

我们这里得到的最佳拟合参数 $w_0$, $w_1$, $w_2$ 依次为 3.76893117e+01, -2.60474147e-018.00078082e-03。也就是说,我们拟合后的函数(保留两位有效数字)为:

$$ y(x) = 37 - 0.26x + 0.0080x^2 \tag{3}$$

然后,我们尝试绘制出拟合后的图像。

# 绘制拟合图像时需要的临时点
x_temp = np.linspace(0, 80, 10000)

# 绘制拟合函数曲线
plt.plot(x_temp, func(parameters[0], x_temp), 'r')

# 绘制原数据点
plt.scatter(x, y)

实现 N 次多项式拟合

你会发现,上面采用 2 次多项式拟合的结果也不能恰当地反映散点的变化趋势。此时,我们可以尝试 3 次及更高次多项式拟合。接下来的代码中,我们将针对上面 2 次多项式拟合的代码稍作修改,实现一个 N 次多项式拟合的方法。

def fit_func(p, x):
    """根据公式,定义 n 次多项式函数
    """
    f = np.poly1d(p)
    return f(x)


def err_func(p, x, y):
    """残差函数(观测值与拟合值之间的差距)
    """
    ret = fit_func(p, x) - y
    return ret


def n_poly(n):
    """n 次多项式拟合
    """
    p_init = np.random.randn(n)  # 生成 n 个随机数
    parameters = leastsq(err_func, p_init, args=(np.array(x), np.array(y)))
    return parameters[0]

可以使用 $n=3$(2 次多项式) 验证一下上面的代码是否可用。

n_poly(3)

此时得到的参数结果和公式 $(3)$ 的结果一致,只是顺序有出入。这是因为 NumPy 中的多项式函数 np.poly1d(3) 默认的样式是:

$$ y(x) = 0.0080x^2 - 0.26x + 37\tag{4}$$

接下来,我们绘制出 3,4,5,6,7, 8 次多项式的拟合结果。

# 绘制拟合图像时需要的临时点
x_temp = np.linspace(0, 80, 10000)

# 绘制子图
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

axes[0, 0].plot(x_temp, fit_func(n_poly(4), x_temp), 'r')
axes[0, 0].scatter(x, y)
axes[0, 0].set_title("m = 3")

axes[0, 1].plot(x_temp, fit_func(n_poly(5), x_temp), 'r')
axes[0, 1].scatter(x, y)
axes[0, 1].set_title("m = 4")

axes[0, 2].plot(x_temp, fit_func(n_poly(6), x_temp), 'r')
axes[0, 2].scatter(x, y)
axes[0, 2].set_title("m = 5")

axes[1, 0].plot(x_temp, fit_func(n_poly(7), x_temp), 'r')
axes[1, 0].scatter(x, y)
axes[1, 0].set_title("m = 6")

axes[1, 1].plot(x_temp, fit_func(n_poly(8), x_temp), 'r')
axes[1, 1].scatter(x, y)
axes[1, 1].set_title("m = 7")

axes[1, 2].plot(x_temp, fit_func(n_poly(9), x_temp), 'r')
axes[1, 2].scatter(x, y)
axes[1, 2].set_title("m = 8")

从上面的 6 张图可以看出,当 $m=4$(4 次多项式) 时,图像拟合的效果已经明显优于 $m=3$ 的结果。但是随着 $m$ 次数的增加,当 $m=8$ 时,曲线呈现出明显的震荡,这也就是线性回归文章中所讲到的过拟和(Overfitting)现象。

使用 scikit-learn 进行多项式拟合

除了像上面我们自己去定义多项式及实现多项式回归拟合过程,也可以使用 scikit-learn 提供的多项式回归方法来完成。这里,我们会用到sklearn.preprocessing.PolynomialFeatures() 这个类。PolynomialFeatures() 主要的作用是产生多项式特征矩阵。如果你第一次接触这个概念,可能需要仔细理解下面的内容。

对于一个二次多项式而言,我们知道它的标准形式为:$ y(x, w) = w_0 + w_1x + w_2x^2 $。但是,多项式回归却相当于线性回归的特殊形式。例如,我们这里令 $x = x_1$, $x^2 = x_2$ ,那么原方程就转换为:$ y(x, w) = w_0 + w_1x_1 + w_2x_2 $,这也就变成了多元线性回归。这就完成了一元高次多项式到多元一次多项式之间的转换。这就类似于高中数学学习过的「换元法」。

举例说明,对于自变量向量 $X$ 和因变量 $y$,如果 $X$:

$$ X=\left[ \begin{array}{c}{2} \\ {-1} \\ {3}\end{array}\right] \tag{5a}$$

我们可以通过 $ y = w_1 x + w_0$ 线性回归模型进行拟合。同样,如果对于一元二次多项式 $ y(x, w) = w_0 + w_1x + w_2x^2 $,如果能得到由 $x = x_1$, $x^2 = x_2$ 构成的特征矩阵,即:

$$ X=\left[X, X^{2}\right]=\left[ \begin{array}{cc}{2} & {4} \\ {-1} & {1} \\ {3} & {9}\end{array}\right] \tag{5b}$$

那么也就可以通过线性回归进行拟合了。

你可以手动计算上面的结果,但是当多项式为一元高次或者多元高次时,特征矩阵的表达和计算过程就变得比较复杂了。例如,下面是二元二次多项式的特征矩阵表达式。

$$ X = \left [ X_{1}, X_{2}, X_{1}^2, X_{1}X_{2}, X_{2}^2 \right ] \tag{5c}$$

还好,在 scikit-learn 中,我们可以通过 PolynomialFeatures() 类自动产生多项式特征矩阵,PolynomialFeatures() 类的默认参数及常用参数定义如下:

sklearn.preprocessing.PolynomialFeatures(degree=2, interaction_only=False, include_bias=True)


- degree: 多项式次数,默认为 2 次多项式

- interaction_only: 默认为 False,如果为 True 则产生相互影响的特征集。

- include_bias: 默认为 True,包含多项式中的截距项。


对应上面的特征向量,我们使用 PolynomialFeatures() 的主要作用是产生 2 次多项式对应的特征矩阵,如下所示:

from sklearn.preprocessing import PolynomialFeatures

X = [2, -1, 3]
X_reshape = np.array(X).reshape(len(X), 1)  # 转换为列向量
# 使用 PolynomialFeatures 自动生成特征矩阵
PolynomialFeatures(degree=2, include_bias=False).fit_transform(X_reshape)

对于上方单元格中的矩阵,第 1 列为 $X^1$,第 2 列为 $X^2$。我们就可以通过多元线性方程 $ y(x, w) = w_0 + w_1x_1 + w_2x_2 $ 对数据进行拟合。

小贴士

本节内容中,你会看到大量的 reshape 操作,它们的目的都是为了满足某些类或函数传参的数组形状。这些操作在这篇文章中是必须的,因为数据原始形状(如上面的一维数组)可能无法直接传入某些特定类或函数中。但在实际工作中并不是必须的,因为你手中的原始数据集形状可能支持直接传入。所以,不必为这些 reshape 操作感到疑惑,也不要死记硬背。

前面小节中的示例数据,其自变量应该是 $x$,而因变量是 $y$。如果我们使用 2 次多项式拟合,那么首先使用 PolynomialFeatures() 得到特征矩阵。

x = np.array(x).reshape(len(x), 1)  # 转换为列向量
y = np.array(y).reshape(len(y), 1)

# 使用 sklearn 得到 2 次多项式回归特征矩阵
poly_features = PolynomialFeatures(degree=2, include_bias=False)
poly_x = poly_features.fit_transform(x)

poly_x

可以看到,输出结果正好对应一元二次多项式特征矩阵公式 $\left [ X, X^2 \right ]$。然后,我们使用 scikit-learn 训练线性回归模型。

from sklearn.linear_model import LinearRegression

# 定义线性回归模型
model = LinearRegression()
model.fit(poly_x, y)  # 训练

# 得到模型拟合参数
model.intercept_, model.coef_

你会发现,这里得到的参数值和公式 $(3,4)$ 一致。为了更加直观,这里同样绘制出拟合后的图像。

# 绘制拟合图像
x_temp = np.array(x_temp).reshape(len(x_temp), 1)
poly_x_temp = poly_features.fit_transform(x_temp)

plt.plot(x_temp, model.predict(poly_x_temp), 'r')
plt.scatter(x, y)

你会发现,上图似曾相识。它和公式 $(3)$ 下方的图其实是一致的。

多项式回归预测

上面的内容中,我们了解了如何使用多项式去拟合数据。那么在接下来的内容中,就使用多项式回归去解决实际的预测问题。本次预测文章中,我们使用到由世界卫生组织和联合国儿童基金会提供的「世界麻疹疫苗接种率」数据集。而目标则是预测相应年份的麻疹疫苗接种率。

首先,我们导入并预览「世界麻疹疫苗接种率」数据集。

import pandas as pd

# 加载数据集
df = pd.read_csv(
    "https://labfile.oss.aliyuncs.com/courses/1081/course-6-vaccine.csv", header=0)
df

可以看出,该数据集由两列组成。其中 Year 表示年份,而 Values 则表示当年世界麻疹疫苗接种率,这里只取百分比的数值部分。我们将数据绘制成图表,查看变化趋势。

# 定义 x, y 的取值
x = df['Year']
y = df['Values']
# 绘图
plt.plot(x, y, 'r')
plt.scatter(x, y)

对于上图呈现出来的变化趋势,我们可能会认为多项式回归会优于线性回归。到底是不是这样呢?试一试便知。

线性回归与 2 次多项式回归对比

根据线性回归中学到的内容,在机器学习任务中,我们一般会将数据集划分为训练集和测试集。所以,这里将 70% 的数据划分为训练集,而另外 30% 则归为测试集。代码如下:

# 首先划分 dateframe 为训练集和测试集
train_df = df[:int(len(df)*0.7)]
test_df = df[int(len(df)*0.7):]

# 定义训练和测试使用的自变量和因变量
X_train = train_df['Year'].values
y_train = train_df['Values'].values

X_test = test_df['Year'].values
y_test = test_df['Values'].values

接下来,我们使用 scikit-learn 提供的多项式回归预测方法来训练模型。首先,我们先解决上面的问题,那就是:多项式回归会不会优于线性回归?

首先,训练线性回归模型,并进行预测。

# 建立线性回归模型
model = LinearRegression()
model.fit(X_train.reshape(len(X_train), 1), y_train.reshape(len(y_train), 1))
results = model.predict(X_test.reshape(len(X_test), 1))
results  # 线性回归模型在测试集上的预测结果

有了预测结果,我们就可以将其同真实的结果进行比较。这里,我们使用到平均绝对误差和均方误差两个指标。如果你对这两个指标仍不太熟悉,它们的定义如下:

平均绝对误差(MAE)就是绝对误差的平均值,它的计算公式 $(6)$ 如下:

$$ \textrm{MAE}(y, \hat{y} ) = \frac{1}{n}\sum_{i=1}^{n}{|y_{i}-\hat y_{i}|}\tag{6}$$

其中,$y_{i}$ 表示真实值,$\hat y_{i}$ 表示预测值,$n$ 则表示值的个数。MAE 的值越小,说明预测模型拥有更好的精确度。

均方误差(MSE)它表示误差的平方的期望值,它的计算公式 $(7)$ 如下:

$$ \textrm{MSE}(y, \hat{y} ) = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y})^{2}\tag{7}$$

其中,$y_{i}$ 表示真实值,$\hat y_{i}$ 表示预测值,$n$ 则表示值的个数。MSE 的值越小,说明预测模型拥有更好的精确度。

这里,我们直接使用 scikit-learn 提供的 MAE 和 MSE 计算方法。

from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error

print("线性回归平均绝对误差: ", mean_absolute_error(y_test, results.flatten()))
print("线性回归均方误差: ", mean_squared_error(y_test, results.flatten()))

接下来,开始训练 2 次多项式回归模型,并进行预测。

# 2 次多项式回归特征矩阵
poly_features_2 = PolynomialFeatures(degree=2, include_bias=False)
poly_X_train_2 = poly_features_2.fit_transform(
    X_train.reshape(len(X_train), 1))
poly_X_test_2 = poly_features_2.fit_transform(X_test.reshape(len(X_test), 1))

# 2 次多项式回归模型训练与预测
model = LinearRegression()
model.fit(poly_X_train_2, y_train.reshape(len(X_train), 1))  # 训练模型

results_2 = model.predict(poly_X_test_2)  # 预测结果

results_2.flatten()  # 打印扁平化后的预测结果
print("2 次多项式回归平均绝对误差: ", mean_absolute_error(y_test, results_2.flatten()))
print("2 次多项式均方误差: ", mean_squared_error(y_test, results_2.flatten()))

根据上面平均绝对误差和均方误差的定义,你已经知道这两个取值越小,代表模型的预测准确度越高。也就是说,线性回归模型的预测结果要优于 2 次多项式回归模型的预测结果。

更高次多项式回归预测

不必惊讶,这种情况是非常常见的。但这并不代表,这节中所讲的多项式回归就会比线性回归更差。下面,我们就试一试 3,4,5次多项式回归的结果。为了缩减代码量,我们重构代码,并一次性得到 3 个预测结果。

这里将通过实例化 make_pipeline 管道类,实现调用一次 fitpredict 方法即可应用于所有预测器。make_pipeline 是使用 sklearn 过程中的技巧创新,其可以将一个处理流程封装起来使用。

具体来讲,例如上面的多项式回归中,我们需要先使用 PolynomialFeatures 完成特征矩阵转换,再放入 LinearRegression 中。那么,PolynomialFeatures + LinearRegression 这一个处理流程,就可以通过 make_pipeline 封装起来使用。

from sklearn.pipeline import make_pipeline

X_train = X_train.reshape(len(X_train), 1)
X_test = X_test.reshape(len(X_test), 1)
y_train = y_train.reshape(len(y_train), 1)

for m in [3, 4, 5]:
    model = make_pipeline(PolynomialFeatures(
        m, include_bias=False), LinearRegression())
    model.fit(X_train, y_train)
    pre_y = model.predict(X_test)
    print("{} 次多项式回归平均绝对误差: ".format(m),
          mean_absolute_error(y_test, pre_y.flatten()))
    print("{} 次多项式均方误差: ".format(m), mean_squared_error(y_test, pre_y.flatten()))
    print("---")

从上面的结果可以得出,3,4,5 次多项式回归的结果均优于线性回归模型。所以,多项式回归还是有其优越性的。

小贴士

深入了解 make_pipeline 的使用,你可以阅读 ↗ 官方文档,或一篇不错的 ↗ 文章


多项式回归预测次数选择



实验进行到现在,你可能会有一个疑问:在选择多项式进行回归预测的过程中,到底几次多项式是最优呢?

对于上面的问题,其实答案很简单。我们可以选择一个误差指标,例如这里选择 MSE,然后计算出该指标随多项式次数增加而变化的图像,结果不就一目了然了吗?试一试。

# 计算 m 次多项式回归预测结果的 MSE 评价指标并绘图
mse = []  # 用于存储各最高次多项式 MSE 值
m = 1  # 初始 m 值
m_max = 10  # 设定最高次数
while m <= m_max:
    model = make_pipeline(PolynomialFeatures(
        m, include_bias=False), LinearRegression())
    model.fit(X_train, y_train)  # 训练模型
    pre_y = model.predict(X_test)  # 测试模型
    mse.append(mean_squared_error(y_test, pre_y.flatten()))  # 计算 MSE
    m = m + 1

print("MSE 计算结果: ", mse)
# 绘图
plt.plot([i for i in range(1, m_max + 1)], mse, 'r')
plt.scatter([i for i in range(1, m_max + 1)], mse)

# 绘制图名称等
plt.title("MSE of m degree of polynomial regression")
plt.xlabel("m")
plt.ylabel("MSE")

如上图所示,MSE 值在 2 次多项式回归预测时达到最高点,之后迅速下降。而 3 次之后的结果虽然依旧呈现逐步下降的趋势,但趋于平稳。一般情况下,我们考虑到模型的泛化能力,避免出现过拟合,这里就可以选择 3 次多项式为最优回归预测模型。

小结

这篇文章中,我们了解了什么是多项式回归,以及多项式回归与线性回归之间的联系与区别。同时,实验探索了动手实现多项式回归拟合,以及运用 scikit-learn 在真实数据集下构建多项式回归预测模型。


系列文章

本篇文章需 特别授权许可,内容版权归作者所有,未经授权,禁止转载。