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

1 插值

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

1.1 插值的定义

设函数$y=f(x)$在区间$[a,b]$上有定义，且已知在点$a\le x_0<x_1<...<x_n\le 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)$的方法称为插值法。

常见的插值法有三种：

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

1.2 多项式插值

定理：设有$n+1$个不同的节点，则存在唯一的$n$次多项式$L_n(x)=a_0+a_{1}x+a_x{x}^2+...+a_n{x}^n$，使得其过所有节点。
简单证明如下：

1.2.1 拉格朗日插值法

对于多项式插值，拉格朗日给出了一个公式：
$L_n(x)={\sum }_{i=0}^n({\Pi }_{j=0,j\ne 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是内部使用的算法：

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

2 拟合

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

2.1 最小二乘法

举个例子，首先我们根据函数$y=5x+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=kx+b$，现在要求的就是$k$和$b$使样本点和拟合曲线最接近。

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

下面就是计算使得$\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2$最小的$k$和$b$了，记作$\hat{k},\hat{b}$，具体计算过程大家可以自行查看相关资料，这里从简，只给出结果：
$\hat{k}=\frac{n\sum _{i=1}^n{x}_i{y}_i-\sum _{i=1}^n{y}_i\sum _{i=1}^n{x}_i}{n\sum _{i=1}^n{x}_i^2-\sum _{i=1}^n{x}_i\sum _{i=1}^n{x}_i}$
$\hat{b}=\frac{\sum _{i=1}^n{x}_i^2\sum _{i=1}^n{y}_i-\sum _{i=1}^n{x}_i\sum _{i=1}^n{x}_i{y}_i}{n\sum _{i=1}^n{x}_i^2-\sum _{i=1}^n{x}_i\sum _{i=1}^n{x}_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 对于拟合的评价

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

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

对于线性拟合函数，有$SST=SSE+SSR$，拟合优度
$0\le R^2=\frac{SSR}{SST}=\frac{SST-SSE}{SST}=1-\frac{SSE}{SST}\le 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 对于线性函数的理解

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

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

3 回归

回归是数据分析中最基础、最终要的分析工具，很多数据分析问题都可以用回归的思想来解决。回归分析的任务是通过研究自变量X和因变量Y的相关关系，尝试去解释Y的形成机制，进而达到通过X去预测Y的目的。其中，相关性是由于我们无法解释事物复杂的因果性，退而求其次改成通过回归分析，研究相关关系。

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

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

本文主要介绍最常用的线性回归。

在预测模型中，面对不同的数据，往往有不同的处理方法：

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

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

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

【待更新】