【数学建模】插值、拟合与回归

HuaweiMian_

已于 2024-07-18 13:30:12 修改

阅读量1.5k

点赞数 24

文章标签：数学建模回归

于 2024-01-25 15:00:44 首次发布

本文链接：https://blog.csdn.net/qq_73176808/article/details/135832294

版权

1 插值

当我们需要根据已知的函数点进行数据处理时，有的时候现有数据非常少，不足以支撑分析，这时就需要使用一些数学方法来**“模拟产生”一些新的、可信度高的值来满足需求，这就是插值**。

1.1 插值的定义

设函数 $y = f (x)$ 在区间 $[a, b]$ 上有定义，且已知在点 $a≤x_0<x_1<...<x_n≤b$ 上的值分别为 $y_0,y_1,...,y_n$ ，若存在以简单函数 $P (x)$ ，使 $P(x_i)=y_i,(i=1,2,...,n)$ ，则称 $P (x)$ 为 $f (x)$ 的插值函数，点 $x_0,x_1,...,x_n$ 称为插值节点，区间 $[a, b]$ 称为插值区间，求 $P (x)$ 的方法称为插值法。

常见的插值法有三种：

多项式插值： $P (x)$ 是次数不超过n的代数多项式，即 $P(x)=a_0+a_1x+...+a_nx^n$ ；
分段插值： $P (x)$ 为分段多项式；
三角插值： $P (x)$ 为三角多项式。
其中，三角插值非常复杂，基本上用不到，因此本文不提及；多项式插值在次数较高的情况下，会在两端处产生极大的不稳定（Runge现象），因此最常用的是分段插值。

1.2 多项式插值

定理：设有 $n + 1$ 个不同的节点，则存在唯一的 $n$ 次多项式 $L_n(x)=a_0+a_1x+a_xx^2+...+a_nx^n$ ，使得其过所有节点。
简单证明如下：
在这里插入图片描述

1.2.1 拉格朗日插值法

对于多项式插值，拉格朗日给出了一个公式：
$L_n(x)=∑_{i=0}^n(Π_{j=0,j≠i}^n\frac{x-x_j}{x_i-x_j})$
但应用拉格朗日插值在高次的情况下，会产生Runge现象，即在两端极其不稳定，震荡明显。因此在不确定曲线运动趋势的情况下，不要使用高次插值。

1.2.2 牛顿插值法

对于多项式插值，牛顿给出的公式为：
$f(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+...+f[x_0,x_1,...,x_{n-1},x_n](x-x_0)(x-x_1)...(x-x_{n-2})(x-x_{n-1})$
其中， $f[x_0,x_k]=\frac{f(x_k)-f(x_0)}{x_k-x_0}$ 为函数 $f (x)$ 关于点 $x_0,x_k$ 的一阶差商（均差），二阶差商 $f[x_0,x_1,x_2]=\frac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}$ ，k阶段差商 $f[x_0,x_1,...,x_k]=\frac{f[x_1,...,x_{k-1},x_k]-f[x_0,x_1,...,x_{k-1}]}{x_k-x_0}$

这两种插值法相比来说，牛顿插值法具有继承性，即它是一个递推的过程，但牛顿插值法也存在Runge现象。同时，两种插值法不能全面反映被插值函数的性态，它们只满足了插值多项式在插值节点处和被插值函数有相同的函数值。在许多实际问题中，更是要求在一个或全部节点上插值多项式与被插值函数有相同的低阶甚至高阶导数值。两种插值法都不能满足这种要求。满足这种要求的插值多项式是埃尔米特插值多项式。

1.3 分段插值

分段插值即选取距离待插值点最近的若干点，使用多项式插值。因此对于整条曲线，使用了多次多项式插值，最终得到一个分段函数。

1.3.1 分段三次埃尔米特插值

埃尔米特插值的具体原理非常复杂，这里就不做过多赘述，对于数学建模，只需要会使用即可。直接使用埃尔米特插值得到的多项式次数较高，也存在Runge现象，因此实际应用中，往往使用分段三次Hermite插值多项式（PCHIP）。
在matlab中，我们直接调用内置函数pchip即可实现，其函数原型为：

p = pchip(x, y, new_x)

其中，x是已知样本点的横坐标，y是已知样本点的纵坐标，new_x是要插入处对应的横坐标。

1.3.2 三次样条插值

三次样条插值也是一种分段插值方法，同样，原理不做过多赘述，我们只需要会使用即可。在matlab中，我们直接调用内置函数spline即可。其函数原型为：

p = spline(x, y, new_x)

其参数和pchip一样，下面我们来实际应用一下这两中插值。

x = -pi: pi;
y = sin(x);
new_x = -pi: 0.1: pi;
p1 = pchip(x, y, new_x);   % 分段三次埃尔米特插值
p2 = spline(x, y, new_x);  % 三次样条插值
plot(x, y, 'o', new_x, p1, '-', new_x, p2, '-');
legend('插值节点', '分段三次埃尔米特插值', '三次样条插值')

得到的结果为：
在这里插入图片描述
可以看到，三次样条插值更加光滑一些。实际中，由于我们不知道数据的生成过程，因此两种插值法都可以使用。

1.4 n维数据的插值

n维数据插值使用较少，这里不做过多赘述，了解如何使用函数接口即可。matlab中内置函数interpn，原型为：

p = interpn(x1, x2, ..., xn, y, new_x1, new_x2, ..., new_xn, method)

其中，x1,x2,…,xn是已知样本点的横坐标，y是已知样本点的纵坐标，new_x1,new_x2,…,new_xn是要插入点的坐标，method是内部使用的算法：

‘linear’：线性插值（默认参数）
‘cubic’：三次插值
'spline：三次样条插值（最精确）
‘nearest’：最邻近插值算法

2 拟合

与插值不同的是，在拟合问题中，曲线不一定过给定的点。拟合的目标是找到一个函数，使得该曲线在某种准则下和所有数据点最近，即拟合的最好，也即最小化损失函数。

2.1 最小二乘法

举个例子，首先我们根据函数 $y = 5 x + 8$ 随机生成一些带有扰动项的点：

clear; clc
% 函数y=5x+8
x = 10 * rand(1, 10);
y = 5 * x + 8 + normrnd(0, 1, 1, 10);
f = @(x) 5 * x + 8;
plot(x, y, 'o')
hold on  % 继续作图
grid on  % 显示网格线
fplot(f, [0, 10])
legend('随机生成的数据', 'y=5x+8')

在这里插入图片描述
在我们不知道原始函数，只知道这些已知点的情况下，我们设置拟合曲线为 $y = k x + b$ ，现在要求的就是 $k$ 和 $b$ 使样本点和拟合曲线最接近。

那么，如何定义最接近呢？首先我们假设 $\hat{y}_i=kx_i+b$ ：
第一种定义是用一次绝对值，即使得 $\sum\limits_{i=1}^n|y_i-\hat{y}_i|$ 最小的 $k$ 和 $b$ ，但是这种算法含有绝对值，不容易求导，因此计算起来比较复杂；
第二种定义用差值的平方，即使得 $\sum\limits_{i=1}^n(y_i-\hat{y}_i)^2$ 最小的 $k$ 和 $b$ ；
对于更高次，首先奇数次会正负抵消，自然是不合理的；偶数高次对于异常值，会对拟合产生极大的影响，因此也不建议使用。因此使用平方最合理，这也就是最小二乘法。

下面就是计算使得 $\sum\limits_{i=1}^n(y_i-\hat{y}_i)^2$ 最小的 $k$ 和 $b$ 了，记作 $\hat{k},\hat{b}$ ，具体计算过程大家可以自行查看相关资料，这里从简，只给出结果：
$\hat{k}=\frac{n\sum\limits_{i=1}^nx_iy_i-\sum\limits_{i=1}^ny_i\sum\limits_{i=1}^nx_i}{n\sum\limits_{i=1}^nx_i^2-\sum\limits_{i=1}^nx_i\sum\limits_{i=1}^nx_i}$
$\hat{b}=\frac{\sum\limits_{i=1}^nx_i^2\sum\limits_{i=1}^ny_i-\sum\limits_{i=1}^nx_i\sum\limits_{i=1}^nx_iy_i}{n\sum\limits_{i=1}^nx_i^2-\sum\limits_{i=1}^nx_i\sum\limits_{i=1}^nx_i}$
其中 $n$ 是已知点的个数。
下面根据计算公式，给出matlab实现代码：

clear; clc
% 函数y=5x+8
x = 10 * rand(10, 1);
y = 5 * x + 8 + normrnd(0, 1, 10, 1);
f = @(x) 5 * x + 8;
plot(x, y, 'o')
hold on
grid on
fplot(f, [0, 10])
xlabel('x')
ylabel('y')

n = size(x, 1);
k = (n * sum(x .* y) - sum(x) * sum(y)) / (n * sum(x .* x) - sum(x) * sum(x))
b = (sum(x .* x) * sum(y) - sum(x) * sum(x .* y)) / (n * sum(x .* x) - sum(x) * sum(x))
f_new = @(x) k * x + b;
hold on
fplot(f_new, [0, 10])
legend('随机生成的数据', 'y=5x+8', '拟合曲线')

下面是运行结果：
在这里插入图片描述

2.2 对于拟合的评价

那么如何评价拟合的好坏呢？
这里首先定义一些概念：

拟合优度（可决系数） $R^2$ ；
总体平方和SST（Total sum of squares）： $SST=\sum\limits_{i=1}^n(y_i-\overline{y})^2$ ；
误差平方和SSE（Sum of squares due to error）： $SSE=\sum\limits_{i=1}^n(y_i-\hat{y})^2$ ；
回归平方和SSR（Sum of squares of the regression）： $SSR=\sum\limits_{i=1}^n(\hat{y}_i-\overline{y})^2$ ；

对于线性拟合函数，有 $SST = SSE + SSR$ ，拟合优度
$0≤R^2=\frac{SSR}{SST}=\frac{SST-SSE}{SST}=1-\frac{SSE}{SST}≤1$
$R^2$ 越接近1，说明误差越小，说明拟合越好。
对于其他函数，直接比较SSE即可。

下面给出计算拟合优度的matlab代码：

y_hat = k * x + b;
SSR = sum((y_hat - mean(y)) .^ 2)
SSE = sum((y_hat - y) .^ 2)
SST = sum((y - mean(y)) .^ 2)
R_2 = SSR / SST

对于本例，我们得到R_2为0.9955（随机数据可能值不相同，但都几乎为1，说明拟合效果好）

2.3 对于线性函数的理解

线性函数分为对变量为线性以及对参数为线性。线性拟合所针对的是对参数为线性的函数。对变量为线性的函数想必大家都了解，那么何为对参数为线性的函数？这里举几个例子：

$y=a+bx^2$ 是对参数为线性的；
$y=e^{a+bx}$ 是对参数为线性的（两侧取对数）；
$y = s in (b + c x)$ 不是对参数为线性的；
$y = ab x$ 不是对参数为线性的。
判断标准概括下来就是：在函数中，参数仅以一次方出现，且不能乘以或除以其他参数，且不能出现参数的复合函数形式（例如例子3）。

3 回归

回归分析的任务是：通过研究自变量X和因变量Y的相关关系，尝试去解释Y的形成机制，进而达到通过X去预测Y的目的。
要完成的任务主要有以下三个：

识别重要变量；
判断相关性的方向；
估计权重（回归系数）。

依据因变量Y的类型，我们将回归分为以下几类：

类型	模型	Y的特点	例子
线性回归	OLS、GLS（最小二乘）	连续数值型变量	GDP、产量、收入
0-1回归	Logistichuigui	二值变量	是否怎样
定序回归	probit定序回归	定序变量	等级评定
计数回归	泊松回归	计数变量	每分钟行人通过数
生存回归	Cox等比例风险回归	生存变量（截断数据）	企业、产品寿命

补充：
不同数据的处理方法：

数据类型	常见建模方法
横截面数据	多元线性回归
时间序列数据	移动平均、指数平滑、ARIMA、GARCH、VAR、协积
面板数据	固定效应和随机效应、静态面板和动态面板

其中前两种数据类型最常见，面板类数据往往需要深入的计量经济学知识。
解释一下这几种数据类型：

横截面数据：在某一时间收集的不同对象的数据；
时间序列数据：对同一对象在不同时间连续观察所取得的数据；
面板数据：横截面数据与时间序列数据综合起来的一种数据资源；

3.1、线性回归

![[一元线性回归.png]]

这里的线性指的是线性于系数。使用线性回归模型建模时，需要对数据进行预处理。也即计算出变量的对数、平方、交叉项等。

3.1.1、变量的内生性

内生性现象：引入了新的自变量后，对回归系数的影响非常大。
内生性解释：假设模型为 $\hat{y}=β_0+β_1x_1+β_2x_2+...+β_kx_k+μ$ ，μ为无法观测且满足一定条件的扰动项，包含了所有与y相关，但未添加到回归模型中的变量。如果满足==误差项μ和所有的自变量x均不相关，则称该回归模型具有外生性==。如果相关，则存在内生性，内生性会导致回归系数估计的不准确，不满足无偏和一致性。

无内生性要求太强，其中解释变量可以分为核心解释变量（关键的变量，当样本容量无限增大时，收敛于待估计参数的真值）和控制变量（不太关键的变量）。
实际中，只要保证核心解释变量和μ不相关即可。

3.1.2、回归系数的解释

![[回归系数的解释.png]]

几种情况：

一元线性回归： $y = a + b x + μ$ ，x每增加一个单位，y平均变化b个单位；
双对数模型： $l n y = a + b l n x + μ$ ，x每增加%1，y平均变化%b；
半对数模型： $y = a + b l n x + μ$ ，x每增加%1，y平均变化b/100个单位； $l n y = a + b x + μ$ ，x每增加1个单位，y平均变化(100b)%。

什么时候取对数：

与市场价值相关的，例如价格、销售额、工资等都可以取；
以年度量的变量，如受教育年限、工作经历等通常不取；
比例变量，例如失业率、参与率等，取不取都行；
变量取值必须是非负数，如果包含0，可以对y取对数 $l n (1 + y)$ 。
取对数的好处：
减弱数据的异方差性；
如果变量本身不服从正态分布，取对数后可能渐进服从正态分布；
可以让模型更具有经济学意义。

含有交互项的自变量
![[交叉项解释.png]]

3.1.3、虚拟变量

![[回归系数解释举例.png]]

![[举例续.png]]

多分类虚拟变量：例如不同省份，可以使用独热编码，模型可以写作 $\sum\limits_{i=1}^{n}β_i×VirtualVariable_i$
注意，为避免完全多重共线性，通常引入的虚拟变量的个数一般是分类数减1。减去的一类作为对照组，其独热编码全部为0。

3.1.4、扰动项满足的条件

扰动项要是球形扰动项，即满足同方差和无自相关两个条件。
![[球型扰动项.png]]

3.1.4.1 异方差

如果扰动项存在异方差，那么：

OLS估计出来的回归系数是无偏的、一致的；
假设检验无法使用（因为构造的统计量无效）；
OLS估计量不再是最优线性无偏估计量。
解决方法：
使用OLS+稳健的标准误（常用）；
使用广义最小二乘GLS（原理：方差较小的数据包含的信息较多，我们可以给予信息量大的数据更大的权重）。

3.1.4.2 检验异方差

可以绘制残差图，大体观察一下残差变化趋势。如果波动非常大，那么说明存在明显异方差，反之则不明显。
回归后运行命令：

rvfplot   // 画残差与拟合值散点图
graph export a1.png, replace  // 保存图片
rvpplot x // 画残差与自变量x散点图
graph export a2.png, replace

拟合值有可能出现负数的原因：因变量分布极其不均衡。

异方差的假设检验（推荐使用怀特）

BP检验
```
// 回归后使用
estat hettest, rhs iid
```
其中，原假设为扰动项不存在异方差，p值小于0.05说明在95%的置信水平下拒绝原假设，即存在异方差。
怀特检验
```
// 回归后使用
estat imtest, white
```

3.1.4.3 处理异方差（推荐使用第一种）

使用OLS+稳健标准误。只要样本容量较大，即使存在异方差，使用稳健标准误，那么所有的参数估计、假设检验均可照常进行。
广义最小二乘GLS。缺点：不知道扰动项的真实协方差矩阵，只能用样本数据来估计，得到的结果不稳健，存在偶然性。

OLS+稳健标准误操作方法：

// 或reg y x1 x2 ..., r
regress y x1 x2 ..., robust

3.1.5、多重共线性

![[多重共线性.png]]

3.1.5.1 检验多重共线性

通过计算方差膨胀因子VIF（Variance Inflation Factor）
假设有k个自变量，那么第m个自变量的 $VIF_m=\frac{1}{1-R^2_{1-k/m}}$ ， $R^2_{1-k/m}$ 是将第m个自变量作为因变量，对剩下k-1个自变量回归得到的拟合优度。
$VIF_m$ 越大，说明第m个自变量和其他变量的相关性越大。
定义回归模型的VIF为 $max\{VIF_1,VIF_2,...VIF_k\}$ 。一般VIF大于10，则认为该回归方程存在严重的多重共线性。

// 回归后使用
setat vif

![[结果.png]]

3.1.5.2 处理多重共线性

如果不关心具体的回归系数，而只关心整个方程预测被解释变量的能力，通常可以不用理会多重共线性（假设整个方程是显著的）。因为多重共线性的主要后果是使得对单个变量的贡献估计不准，但所有变量的整体效果仍可以较准确估计。
如果关心具体回归系数，但多重共线性不影响所关心变量的显著性，也可以忽略。即使在有方差膨胀的情况下，这些系数依然显著，没有多重共线性只会更加显著。
如果多重共线性严重影响所关心变量的显著性，则需要增大样本容量，删除导致严重共线性的变量，不要轻易删除（可能会导致内生性）；或者修改模型。

3.1.5.3 逐步回归

向前逐步回归Forward selection：将自变量逐个引入模型，每引入一个自变量后都进行检验，显著时才加入模型。
缺点：后来引入的自变量可能导致原来显著的自变量不显著了。

stepwise regress y x1 x2 ..., pe(#1)
// 筛选后变量太多，减小#1；反之则增加#1

向后逐步回归Backward elimination（推荐使用）：先将所有自变量放入模型，之后尝试将其中一个自变量剔除，看整个模型解释因变量的能力是否有显著变化，之后将最没有解释力的那个自变量剔除。此过程逐渐迭代，直到没有自变量符合条件。
缺点：计算量较大。

stepwise regress y x1 x2 ..., pr(#2)
// 筛选后变量太多，减小#2；反之则增加#2

注意：

x1，x2，…之间不能有完全多重共线性，如果有，需要自己手动剔除；
可以加上参数b（标准化回归系数）和r（稳健标准误）。
迫不得已的情况下再使用逐步回归，因为很容易导致内生性。

HuaweiMian_

关注

24
点赞
踩
28

收藏

觉得还不错? 一键收藏
1
评论
【数学建模】插值、拟合与回归

本文主要介绍数学建模中常用的插值、拟合与回归算法。
复制链接

扫一扫