MATLAB 数据拟合（使用 polyfit 多项式曲线拟合、polyval）

文章共11,265字 · 阅读需要大约38分钟

一键AI生成摘要，助你高效阅读

问答

TXYGoodluck

35163人浏览 · 2020-08-03 10:25:30

TXYGoodluck · 2020-08-03 10:25:30 发布

解决数据拟合问题最重要方法是最小二乘法和回归分析。如，我们需要从一组测定的数据（例如N个点（xi，yi）（i=0,1，…，m））去求得自变量 x 和因变量 y 的一个近似解表达式 y=f（x），这就是由给定的 N 个点（xi，yi）（i=0,1，…，m）求数据拟合的问题。（注意数据拟合和数据插值是不同的，举个例子：因为测量数据往往不可避免地带有测试误差，而插值多项式又通过所有的点（xi，yi），这样就使插值多项式保留了这些误差，从而影响逼近精度，使得插值效果不理想）
所以使用最小二乘法曲线拟合法：即寻求已知函数的一个逼近函数y=f（x），使得逼近函数从总体与已知函数的偏差按某种方法度量能达到最小，而又不一定通过全部的点（xi，yi），这个时候就需要使用最小二乘法曲线拟合法。
数据拟合的具体做法是：对给定的数据（xi，yi）（i=0,1，…，m），在取定的函数类 $\phi$ 中使误差 $r_{i}=p\left(x_{i}\right)-y_{i}(i=0,1, \ldots, m)$ 的平方和最小，即
$\left[\sum_{i=0}^{m} r_{i}^{2}=\sum_{i=0}^{m}\left[p\left(x_{i}-y_{i}\right)\right]^{2}\right]_{\min }$
从几何意义讲，即寻求与给定点 $x_{i}-y_{i}(i=0,1, \ldots, m)$ 的距离平方和为最小的曲线y=p（x）。函数p（x）称为拟合函数或最小二乘解，求拟合函数p（x）的方法称为曲线拟合的最小二乘法。
在曲线拟合中，函数类 $\phi$ 可有不同的选取方法。
MATLAB工具箱中提供了最小二乘拟合函数 polyfit() -->多项式曲线拟合
具体调用格式有三种：

P = polyfit(X,Y,N)
[P,S] = polyfit(X,Y,N)
[P,S,MU] = polyfit(X,Y,N)

（1）P = polyfit(X,Y,N) 返回次数为 n 的多项式 p(x) 的系数，该阶数是 y 中数据的最佳拟合（在最小二乘方式中）。p 中的系数按降幂排列，p 的长度为 n+1.
其中 X 为输入的向量x，Y 为得到的函数值，N 表示拟合的最高次数，返回的P值为拟合的多项式：

$\times X^{n}+P(2) \times X^{n-1}+\ldots+P(N) \times X+P(N+1)$

（2） [P,S] = polyfit(X,Y,N)还返回一个结构体 S，S可用作 polyval 的输入来获取误差估计值（S为由范德蒙矩阵的QR分解的R分量）。
其中 X 为输入的向量x，Y 为得到的函数值，N 表示拟合的最高次数，返回的P值为拟合的多项式：

$\times X^{n}+P(2) \times X^{n-1}+\ldots+P(N) \times X+P(N+1)$

（3）[P,S,MU] = polyfit(X,Y,N)还返回 mu，mu是一个二元素向量，包含中心化值和缩放值。mu(1) 是 mean(x)，mu(2) 是 std(x)。使用这些值时，polyfit 将 x 的中心置于零值处并缩放为具有单位标准差：

$\hat{x}=\frac{x-\bar{x}}{\sigma_{x}}$

这种中心化和缩放变换可同时改善多项式和拟合算法的数值属性。

下面来看一些具体的例子：（来源于帮助文档，改编）

将多项式与三角函数拟合：
在区间 [0,4*pi] 中沿余弦曲线生成 10 个等间距的点。

>> x = linspace(0,4*pi,10);
>> y = cos(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,'m')
>> hold off

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

将多项式与点集拟合
创建一个由区间 [0,1] 中的 5 个等间距点组成的向量，并计算这些点处的 y（x）= (2+x)^-1 。

>> x=0:0.2:1;   //或者写成  >> x=0:0.2:1;  也可以实现相同的作用
>> y = 1./(2+x);

将 4 次多项式与 5 个点拟合。通常，对于 n 个点，可以拟合 n-1 次多项式以便完全通过这些点。

>> p = polyfit(x,y,4);

在由 0 和 2 之间的点组成的更精细网格上计算原始函数和多项式拟合。

>> x1 = linspace(0,2);
>> y1 = 1./(2+x1);
>> f1 = polyval(p,x1);

在更大的区间 [0,2] 中绘制函数值和多项式拟合，其中包含用于获取以圆形突出显示的多项式拟合的点。多项式拟合在原始 [0,1] 区间中的效果较好，但在该区间外部很快与拟合函数出现差异。

>> figure
>> plot(x,y,'o')
>> hold on
>> plot(x1,y1)
>> plot(x1,f1,'b--')
>> legend('y','y1','f1')

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

对误差函数进行多项式拟合：
首先生成 x 点的向量，在区间 [0,1.5*pi] 内等间距分布；然后计算这些点处的 erf(x)。
注：erf 是误差函数，也称高斯误差函数。

x = (0:0.10:1.5*pi)';
 y=erf(x);

确定6次逼近多项式的系数。

p=polyfit(x,y,6)

p =

0.0019   -0.0259    0.1242   -0.1785   -0.3442    1.2684   -0.0095

为了查看拟合情况如何，在各数据点处计算多项式，并生成说明数据、拟合和误差的一个表。

f=polyval(p,x);
T=table(x,y,f,y-f,'VariableNames',{'X','Y','Fit','FitError'})

T =

48×4 table

 X        Y          Fit         FitError  
___    _______    __________    ___________

  0          0    -0.0094811      0.0094811
0.1    0.11246       0.11375     -0.0012876
0.2     0.2227       0.22919     -0.0064911
0.3    0.32863       0.33619     -0.0075599
0.4    0.42839       0.43431      -0.005914
0.5     0.5205       0.52334     -0.0028408
0.6    0.60386       0.60326     0.00059392
0.7     0.6778        0.6742      0.0035981
0.8     0.7421       0.73643      0.0056695
0.9    0.79691       0.79033      0.0065794
  1     0.8427       0.83637      0.0063316
1.1    0.88021        0.8751      0.0051058
1.2    0.91031       0.90712       0.003194
1.3    0.93401       0.93307     0.00093826
1.4    0.95229       0.95361      -0.001323
1.5    0.96611        0.9694     -0.0032971
...

在该区间中，插值与实际值非常符合。创建一个绘图，以显示在该区间以外，外插值与实际数据值如何快速偏离。

 x1=(0:0.10:2*pi)';
 y1=erf(x1);
 f1=polyval(p,x1);
 plot(x,y,'o')
 plot(x1,y1,'-')
 hold on
 plot(x1,y1,'-')
 plot(x1,f1,'g-')
 axis([0  5.5  0  3])

在这里插入图片描述

使用中心化和缩放改善数值属性
创建一个由 1750 - 2000 年的人口数据组成的表，并绘制数据点。

 year = (1750:25:2000)';
 pop = 1e6*[791 856 978 1050 1262 1544 1650 2532 6122 8170 11560]';
 T = table(year, pop)

T =

11×2 table

year       pop   
____    _________

1750     7.91e+08
1775     8.56e+08
1800     9.78e+08
1825     1.05e+09
1850    1.262e+09
1875    1.544e+09
1900     1.65e+09
1925    2.532e+09
1950    6.122e+09
1975     8.17e+09
2000    1.156e+10

 plot(year,pop,'square')

在这里插入图片描述
使用带三个输入的 polyfit 拟合一个使用中心化和缩放的 5 次多项式，这将改善问题的数值属性。polyfit 将 year 中的数据以 0 为进行中心化，并缩放为具有标准差 1，这可避免在拟合计算中出现病态的 Vandermonde （范德蒙）矩阵。

[p,~,mu] = polyfit(T.year, T.pop, 5);

使用带四个输入的 polyval，根据缩放后的年份 (year-mu(1))/mu(2) 计算 p。绘制结果对原始年份的图。

f = polyval(p,year,[],mu);
hold on
plot(year,f)

在这里插入图片描述

简单线性回归
将一个简单线性回归模型与一组离散二维数据点拟合。
创建几个由样本数据点 (x,y) 组成的向量。对数据进行一次多项式拟合。

>> x=1:40;
>> y=0.4*x-1.5*randn(1,60);
注意：这里 矩阵维度必须一致。 所以第二行代码的输入是错误的
 
 y=0.4*x-1.5*randn(1,40);
 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% 预测区间 $\pm 2 \Delta$ 。

plot(x,y,'go')
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 – 查询点（向量）
查询点，指定为一个向量。x 中的点对应于 y 中包含的拟合函数值。如果 x 不是向量，则 polyfit 将其转换为列向量 x( : )。😃
x 具有重复（或接近重复）的点或者如果 x 可能需要中心化和缩放时的警告消息结果。
**数据类型：**single | double
复数支持： 是
y – 查询点位置的拟合值（向量）
查询点位置的拟合值，指定为向量。y 中的值对应于 x 中包含的查询点。如果 y 不是向量，则 polyfit 将其转换为列向量 y ( : ) 。😃
**数据类型：**single | double
复数支持： 是
n --多项式拟合的次数（正整数标量）
多项式拟合的次数，指定为正整数标量。n 指定 p 中最左侧系数的多项式幂。
输出参数解释：
p – 最小二乘拟合多项式系数（向量）

最小二乘拟合多项式系数，以向量的形式返回。p 的长度为 n+1，包含按降幂排列的多项式系数，最高幂为 n。如果 x 或 y 包含 NaN 值且 n < length(x)，则 p 的所有元素均为 NaN。

使用 polyval 计算 p 在查询点处的解。

S – 误差估计结构体（结构体）
误差估计结构体。此可选输出结构体主要用作 polyval 函数的输入，以获取误差估计值。S 包含以下字段：

字段	说明
R	Vandermonde 矩阵 x 的 QR 分解的三角因子
df	自由度
normr	残差的范数

如果 y 中的数据是随机的，则 p 的估计协方差矩阵是 (Rinv*Rinv’)*normr^2/df，其中 Rinv 是 R 的逆矩阵。

如果 y 中数据的误差呈独立正态分布，并具有常量方差，则 [y,delta] = polyval(…) 可生成至少包含 50% 的预测值的误差边界。即 y ± delta 至少包含 50% 对 x 处的未来观测值的预测值。

mu – 中心化值和缩放值

中心化值和缩放值，以二元素向量形式返回。mu(1) 为 mean(x)，mu(2) 为 std(x)。这些值以单位标准差将 x 中的查询点的中心置于零值处。

使用 mu 作为 polyval 的第四个输入以计算 p 在缩放点 (x - mu(1))/mu(2) 处的解。

局限性

在涉及很多点的问题中，使用 polyfit 增加多项式拟合的次数并不总能得到较好的拟合。高次多项式可以在数据点之间振动，导致与数据之间的拟合较差。在这些情况下，可使用低次多项式拟合（点之间倾向于更平滑）或不同的方法，具体取决于该问题。
多项式在本质上是无边界的振荡函数。所以，它们并不非常适合外插有界的数据或单调（递增或递减）的数据。

算法：
polyfit 使用 x 构造具有 n+1 列和 m = length(x) 行的 Vandermonde 矩阵 V 并生成线性方程组
$\left(\begin{array}{cccc} x_{1}^{n} & x_{1}^{n-1} & \ldots & 1 \\ n_{2} & x_{2}-1 & \ldots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ n_{m} & x_{m}-1 & \ldots & 1 \end{array}\right)\left(\begin{array}{c} p_{1} \\ p_{2} \\ \vdots \\ p_{n+1} \end{array}\right)=\left(\begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{m} \end{array}\right)$