插值与拟合
插值与拟合
插值与拟合是数学建模中的一种基本的数据分析手段，被公认为建模中常用的算法之一。
插值问题
已知函数在某区间（域）内若干点处的值，求函数在该区间（域）内其他点处的值，这种问题适宜用插值方法解决。
一维插值问题可描述为：已知函数在\(x_0,x_1,…,x_n\)处的值\(y_1,y_2,…,y_n\),求简单函数\(p(x)\)，使\(p(x_i)=y_i\)。
可以用范德蒙行列式和克莱姆法则证明，在\(x_0,x_1,…,x_n\)处的、取值\(y_1,y_2,…,y_n\)的多项式存在且唯一，即插值问题的解唯一存在。
常用的插值方法有Lagrange插值法和Newton插值法
拉格朗日插值法
拉格朗日插值公式指的是在节点上给出节点基函数，然后做基函数的2线性组合，组合系数为节点函数值的一种插值多项式。
高次插值的Runge现象
Runge通过对一个例子的研究发现，插值多项式次数越高，插值精度越高的结论仅仅在插值多项式的次数不超过七时成立；插值多项式的次数超过七时，插值多项式会出现严重的振荡现象，称之为Runge现象。
因此，在实际中不应使用七次以上的插值。
避免Runge现象的常用方法是：将插值区间分成若干小区间，在小区间内用低次（二次，三次）插值，即分段低次插值，如样条函数插值。
Matlab中的插值
1.一维插值
一维插值的命令是interp1
语法及说明
vq = interp1(x,v,xq)
%使用线性插值返回一维函数在特定查询点的插入值。向量 x 包含样本点，v 包含对应值 v(x)。向量 xq 包含查询点的坐标。如果您有多个在同一点坐标采样的数据集，则可以将 v 以数组的形式进行传递。数组 v 的每一列都包含一组不同的一维样本值。
vq = interp1(x,v,xq,method)
% 指定备选插值方法：'linear'、'nearest'、'next'、'previous'、'pchip'、'cubic'、'v5cubic'、'makima' 或 'spline'。默认方法为 'linear'。
vq = interp1(x,v,xq,method,extrapolation)
% 用于指定外插策略，来计算落在 x 域范围外的点。如果希望使用 method 算法进行外插，可将 extrapolation 设置为 'extrap'。您也可以指定一个标量值，这种情况下，interp1 将为所有落在 x 域范围外的点返回该标量值。
vq = interp1(v,xq)
%返回插入的值，并假定一个样本点坐标默认集。默认点是从 1 到 n 的数字序列，其中 n 取决于 v 的形状：
%当 v 是向量时，默认点是 1:length(v)。
%当 v 是数组时，默认点是 1:size(v,1)。
%如果您不在意点之间的绝对距离，则可使用此语法。
vq = interp1(v,xq,method)
%指定备选插值方法中的任意一种，并使用默认样本点。
vq = interp1(v,xq,method,extrapolation)
%指定外插策略，并使用默认样本点。
pp = interp1(x,v,method,'pp')
%使用 method 算法返回分段多项式形
示例
定义样本点 x 及其对应样本值 v。
x = 0:pi/4:2*pi;
v = sin(x);
将查询点定义为 x 范围内更精细的采样点。
xq = 0:pi/16:2*pi;
在查询点插入函数并绘制结果。
figure
vq1 = interp1(x,v,xq);
plot(x,v,'o',xq,vq1,':.');
xlim([0 2*pi]);
title('(Default) Linear Interpolation');
现在使用 'spline' 方法计算相同点处的 v。
figure
vq2 = interp1(x,v,xq,'spline');
plot(x,v,'o',xq,vq2,':.');
xlim([0 2*pi]);
title('Spline Interpolation');
定义一组函数值。
v = [0 1.41 2 1.41 0 -1.41 -2 -1.41 0];
定义一组介于默认点 1:9 之间的查询点。在这种情况下，默认点为 1:9，因为 v 包含 9 个值。
xq = 1.5:8.5;
计算 xq 处的 v。
vq = interp1(v,xq);
绘制结果。
figure
plot((1:9),v,'o',xq,vq,'*');
legend('v','vq');
2.二维插值
二维插值的命令是interp2。
语法及说明
Vq = interp2(X,Y,V,Xq,Yq)
%使用线性插值返回双变量函数在特定查询点的插入值。结果始终穿过函数的原始采样。X 和 Y 包含样本点的坐标。V 包含各样本点处的对应函数值。Xq 和 Yq 包含查询点的坐标。
Vq = interp2(V,Xq,Yq)
%假定一个默认的样本点网格。默认网格点覆盖矩形区域 X=1:n 和 Y=1:m，其中 [m,n] = size(V)。如果您希望节省内存且不在意点之间的绝对距离，则可使用此语法。
Vq = interp2(V)
%将每个维度上样本值之间的间隔分割一次，形成优化网格，并在这些网格上返回插入值。
Vq = interp2(V,k)
%将每个维度上样本值之间的间隔反复分割 k 次，形成优化网格，并在这些网格上返回插入值。这将在样本值之间生成 2^k-1 个插值点。
Vq = interp2(___,method)
%指定备选插值方法：'linear'、'nearest'、'cubic'、'makima' 或 'spline'。默认方法为 'linear'。
Vq = interp2(___,method,extrapval)
% 还指定标量值 extrapval，此参数会为处于样本点域范围外的所有查询点赋予该标量值。
%如果您为样本点域范围外的查询省略 extrapval 参数，则基于 method 参数，interp2 返回下列值之一：
%对于 'spline' 和 'makima' 方法，返回外插值
%对于其他内插方法，返回 NaN 值
示例
对 peaks 函数进行粗略采样。
[X,Y] = meshgrid(-3:3);
V = peaks(X,Y);
绘制粗略采样。
figure
surf(X,Y,V)
title('Original Sampling');
创建间距为 0.25 的查询网格。
[Xq,Yq] = meshgrid(-3:0.25:3);
对查询点插值。
Vq = interp2(X,Y,V,Xq,Yq);
绘制结果。
figure
surf(Xq,Yq,Vq);
title('Linear Interpolation Using Finer Grid');
对 peaks 函数进行粗略采样。
[X,Y] = meshgrid(-3:3);
V = peaks(7);
绘制粗略采样。
figure
surf(X,Y,V)
title('Original Sampling');
创建间距为 0.25 的查询网格。
[Xq,Yq] = meshgrid(-3:0.25:3);
对查询点插值，并指定三次插值。
Vq = interp2(X,Y,V,Xq,Yq,'cubic');
绘制结果。
figure
surf(Xq,Yq,Vq);
title('Cubic Interpolation Over Finer Grid');
拟合问题
根据离散数据求数据间近似函数关系的问题称为曲线拟合问题。
拟合问题与插值问题的区别
- 插值函数过已知点，而拟合函数不一定过已知点；
- 插值主要用于求函数值，而拟合的主要目的是求函数关系，从而进行预测等进一步的分析。
拟合的计算
曲线拟合需解决如下两个问题：
- 线型的选择；
- 线型中参数的计算。
线型的选择是拟合计算的关键和难点。通常主要根据专业知识和散点图确定线型。
线性拟合中参数的计算可采用最小二乘法，而非线性拟合参数的计算则要应用Gauss-Newton迭代法。
Matlab中的拟合
1.多项式拟合——polyfit
语法及说明
p = polyfit(x,y,n)
%返回次数为 n 的多项式 p(x) 的系数，该阶数是 y 中数据的最佳拟合（在最小二乘方式中）。p 中的系数按降幂排列，p 的长度为 n+1
[p,S] = polyfit(x,y,n)
%还返回一个结构体 S，后者可用作 polyval 的输入来获取误差估计值。
[p,S,mu] = polyfit(x,y,n)
%还返回 mu，后者是一个二元素向量，包含中心化值和缩放值。mu(1) 是 mean(x)，mu(2) 是 std(x)。使用这些值时，polyfit 将 x 的中心置于零值处并缩放为具有单位标准差
示例
在区间 [0,4*pi] 中沿正弦曲线生成 10 个等间距的点。
x = linspace(0,4*pi,10);
y = sin(x);
使用 polyfit 将一个 7 次多项式与这些点拟合。
p = polyfit(x,y,7);
在更精细的网格上计算多项式并绘制结果图。
x1 = linspace(0,4*pi);
y1 = polyval(p,x1);
figure
plot(x,y,'o')
hold on
plot(x1,y1)
hold off
创建一个由区间 [0,1] 中的 5 个等间距点组成的向量，并计算这些点处的 y(x)=(1+x)−1。
x = linspace(0,1,5);
y = 1./(1+x);
将 4 次多项式与 5 个点拟合。通常，对于 n 个点，可以拟合 n-1 次多项式以便完全通过这些点。
p = polyfit(x,y,4);
在由 0 和 2 之间的点组成的更精细网格上计算原始函数和多项式拟合。
x1 = linspace(0,2);
y1 = 1./(1+x1);
f1 = polyval(p,x1);
在更大的区间 [0,2] 中绘制函数值和多项式拟合，其中包含用于获取以圆形突出显示的多项式拟合的点。多项式拟合在原始 [0,1] 区间中的效果较好，但在该区间外部很快与拟合函数出现差异。
figure
plot(x,y,'o')
hold on
plot(x1,y1)
plot(x1,f1,'r--')
legend('y','y1','f1')
将一个简单线性回归模型与一组离散二维数据点拟合。
创建几个由样本数据点 (x,y) 组成的向量。对数据进行一次多项式拟合。
x = 1:50;
y = -0.3*x + 2*randn(1,50);
p = polyfit(x,y,1);
计算在 x 中的点处拟合的多项式 p。用这些数据绘制得到的线性回归模型。
f = polyval(p,x);
plot(x,y,'o',x,f,'-')
legend('data','linear fit')
将一个线性模型拟合到一组数据点并绘制结果，其中包含预测区间为 95% 的估计值。
创建几个由样本数据点 (x,y) 组成的向量。使用 polyfit 对数据进行一次多项式拟合。指定两个输出以返回线性拟合的系数以及误差估计结构体。
x = 1:100;
y = -0.3*x + 2*randn(1,100);
[p,S] = polyfit(x,y,1);
计算以 p 为系数的一次多项式在 x 中各点处的拟合值。将误差估计结构体指定为第三个输入，以便 polyval 计算标准误差的估计值。标准误差估计值在 delta 中返回。
[y_fit,delta] = polyval(p,x,S);
绘制原始数据、线性拟合和 95% 预测区间 y±2Δ。
plot(x,y,'bo')
hold on
plot(x,y_fit,'r-')
plot(x,y_fit+2*delta,'m--',x,y_fit-2*delta,'m--')
title('Linear Fit of Data with 95% Prediction Interval')
legend('Data','Linear Fit','95% Prediction Interval')
首先生成 x 点的向量，在区间 [0,2.5] 内等间距分布；然后计算这些点处的 erf(x)。
x = (0:0.1:2.5)';
y = erf(x);
确定 6 次逼近多项式的系数。
p = polyfit(x,y,6)
p = 1×7
0.0084 -0.0983 0.4217 -0.7435 0.1471 1.1064 0.0004
为了查看拟合情况如何，在各数据点处计算多项式，并生成说明数据、拟合和误差的一个表。
f = polyval(p,x);
T = table(x,y,f,y-f,'VariableNames',{'X','Y','Fit','FitError'})
T=26×4 table
X Y Fit FitError
___ _______ __________ ___________
0 0 0.00044117 -0.00044117
0.1 0.11246 0.11185 0.00060836
0.2 0.2227 0.22231 0.00039189
0.3 0.32863 0.32872 -9.7429e-05
0.4 0.42839 0.4288 -0.00040661
0.5 0.5205 0.52093 -0.00042568
0.6 0.60386 0.60408 -0.00022824
0.7 0.6778 0.67775 4.6383e-05
0.8 0.7421 0.74183 0.00026992
0.9 0.79691 0.79654 0.00036515
1 0.8427 0.84238 0.0003164
1.1 0.88021 0.88005 0.00015948
1.2 0.91031 0.91035 -3.9919e-05
1.3 0.93401 0.93422 -0.000211
1.4 0.95229 0.95258 -0.00029933
1.5 0.96611 0.96639 -0.00028097
⋮
在该区间中，插值与实际值非常符合。创建一个绘图，以显示在该区间以外，外插值与实际数据值如何快速偏离。
x1 = (0:0.1:5)';
y1 = erf(x1);
f1 = polyval(p,x1);
figure
plot(x,y,'o')
hold on
plot(x1,y1,'-')
plot(x1,f1,'r--')
axis([0 5 0 2])
hold off
局限性
- 在涉及很多点的问题中，使用
polyfit增加多项式拟合的次数并不总能得到较好的拟合。高次多项式可以在数据点之间振动，导致与数据之间的拟合较差。在这些情况下，可使用低次多项式拟合（点之间倾向于更平滑）或不同的方法，具体取决于该问题。 - 多项式在本质上是无边界的振荡函数。所以，它们并不非常适合外插有界的数据或单调（递增或递减）的数据。
2.Matlab拟合工具箱
为了更好、更便捷地进行拟合，Matlab提供了拟合工具箱。
介绍
Curve Fitting 应用程序提供了一个灵活的界面，您可以在其中以交互方式将曲线和曲面拟合到数据并查看图表。您可以： 创建、绘制和比较多个拟合。使用线性或非线性回归、插值、平滑和自定义方程。查看拟合优度统计数据、显示置信区间和残差、移除异常值并使用验证数据评估拟合。自动生成代码以拟合和绘制曲线和曲面，或将拟合导出到工作区以进行进一步分析。
使用方法
（1）.工具箱的启动
在命令窗口输入cftool即可启动工具箱。
（2）.数据的录入
在命令窗口录入自变量x和函数y的数据，然后在Data菜单中即可选中上述数据，并产生Data sets。
此时工具箱会自动画出散点图。
（3）.拟合
点击Fitting->New fit，可以修改Fit name，选择Data sets（自动）和Type of Fit。
Apply后即可完成拟合。
