多项式拟合polyfit之测试

最新推荐文章于 2025-06-12 14:46:29 发布

原创最新推荐文章于 2025-06-12 14:46:29 发布 · 2.3k 阅读

2 ·

CC 4.0 BY-SA版权

文章标签：

#拟合曲面

图像处理同时被 2 个专栏收录

4 篇文章

订阅专栏

matlab闲玩

3 篇文章

订阅专栏

本文介绍最小二乘法多项式曲线拟合的基本原理及应用，对比不同拟合方法，详细推导了拟合过程，并通过Matlab代码验证了方法的有效性和计算效率。

部分内容转载自：https://blog.youkuaiyun.com/JairusChan/article/details/7517773

概念

最小二乘法多项式曲线拟合，根据给定的m个点,并不要求这条曲线精确地经过这些点，而是曲线y=f(x)的近似曲线y= φ(x)。

原理

[原理部分由个人根据互联网上的资料进行总结，希望对大家能有用]

给定数据点pi(xi,yi)，其中i=1,2,…,m。求近似曲线y= φ(x)。并且使得近似曲线与y=f(x)的偏差最小。近似曲线在点pi处的偏差δi= φ(xi)-y，i=1,2,…,m。

常见的曲线拟合方法:

1.使偏差绝对值之和最小

2.使偏差绝对值最大的最小

3.使偏差平方和最小

按偏差平方和最小的原则选取拟合曲线，并且采取二项式方程为拟合曲线的方法,称为最小二乘法。

推导过程：

1. 设拟合多项式为：

2. 各点到这条曲线的距离之和，即偏差平方和如下：

3. 为了求得符合条件的a值，对等式右边求ai偏导数，因而我们得到了：

…….

4. 将等式左边进行一下化简，然后应该可以得到下面的等式：

…….

5. 把这些等式表示成矩阵的形式，就可以得到下面的矩阵：

6. 将这个范德蒙得矩阵化简后可得到:

7. 也就是说X*A=Y，那么A = (X’*X)-1*X’*Y，便得到了系数矩阵A，同时，我们也就得到了拟合曲线。

以上内容引用自：https://blog.youkuaiyun.com/JairusChan/article/details/7517773
+++++++++++++++++++++++++++++++++++++
以下内容为博客作者hmmwjs添加
理论推导到步骤5，后面的化简不知是怎么推导出来的。
通过数据进行了测试，
步骤5因为double精度问题（在三次多项式的时候，已经达到了6次方）；
步骤6结果与matlab自带的polyfit计算结果完全一致，且计算时间上也比较接近；
至于步骤7，方程没错（‘表示转置，-1表示矩阵求逆），但运行时间上变慢。
下面附上，matlab测试的代码

% 用来测试，【计算时间】
% 转换过程是对的：5对称矩阵（关于主对角线对称的方阵）-->6范德蒙行列式（第一列都是1）
% https://blog.youkuaiyun.com/JairusChan/article/details/7517773
clc;clear all;close all
format long g
warning off
x = [405.0000 , 404.0000 , 403.0000 , 402.0000 , 401.0000 , 400.0000 , 399.0000 , 398.0000 , 373.0000 , 372.0000 , 371.0000 , 370.0000 , 369.0000 , 368.0000 , 367.0000 , 366.0000 , 365.0000 , 364.0000 , 356.0000 , 355.0000 , 354.0000 , 353.0000 , 352.0000 , 351.0000 , 350.0000 , 347.0000 , 346.0000 , 343.0000];
y = [-49.5, -47.5, -46.5, -45.5, -44.5, -42.5, -41.5, -39.5,  -4.5,  -3.5,  -1.5,   0.5,   1.5,   3.5,   4.5,   6.5,   7.5,   8.5,  25.5,  27.5,  29.5,  31.5,  32.5,  37.5,  37.5,  37.5,  37.5,  37.5];
%% 拟合测试——以三次方拟合为例——y=a0+a1*x+a2*x^2+a3*x^3;共四个系数
nTimes = 100
%% 利用自己推导的，网页中的步骤5
tic
t(1) = toc;
disp('自己推导：')
for times = 1:nTimes
for i = 0:6
    eval(['x',num2str(i),' = x.^',num2str(i),';']);
    eval(['xSum',num2str(i),' = sum(x',num2str(i),'(:));']);
    eval(['yx',num2str(i),'=y.*x.^',num2str(i),';']);
    eval(['yxSum',num2str(i),'=sum(yx',num2str(i),'(:));']);
end
XX = [
    xSum0 xSum1 xSum2 xSum3
    xSum1 xSum2 xSum3 xSum4
    xSum2 xSum3 xSum4 xSum5
    xSum3 xSum4 xSum5 xSum6
    ];
YY = [
    yxSum0
    yxSum1
    yxSum2
    yxSum3
    ];
% disp('自己推导：')
A = XX\YY;
end
% A1 = inv(XX)*YY%结果与上面的左除一样
t(2) = toc;
disp('Matlab自带拟合：')
for times = 1:nTimes
%% 利用matlab自带的拟合
% disp('Matlab自带拟合：')
xishu = polyfit(x,y,3);
end
t(3) = toc;
disp('步骤6：')
for times = 1:nTimes
%% 利用网页上提供的方法 步骤6
% 测试表明，人家的公式是对的
for i = 0:3
    eval(['x_',num2str(i),'=(x'').^',num2str(i),';']);
end
XX2 = [x_0 x_1 x_2 x_3 ];
YY2 = y';
% disp('步骤6：')
A2 = XX2\YY2;
end
t(4) = toc;
disp('步骤7：')
for times = 1:nTimes
%% 利用网页上提供的方法  步骤7
for i = 0:3
    eval(['x_',num2str(i),'=(x'').^',num2str(i),';']);
end
XX2 = [x_0 x_1 x_2 x_3 ];
YY2 = y';
A21 = inv(XX2'*XX2)*XX2'*YY2;
end
t(5) = toc;
T = t(2:end)-t(1:end-1);
T'

（运行一百次耗时:s）计算结果如下：

nTimes =
   100
自己推导：
Matlab自带拟合：
步骤6：
步骤7：
ans =
         0.439228266018257
        0.0657886020166998
        0.0709421622497514
         0.075328753289399

附带之前运行时的计算结果

% 用来测试[计算结果  testPolyfitEquation0为计算时间的测试]
% ，转换过程是对的：5对称矩阵（关于主对角线对称的方阵）-->6范德蒙行列式（第一列都是1）
% https://blog.youkuaiyun.com/JairusChan/article/details/7517773
clc;clear all;close all
format long g
x = [405.0000 , 404.0000 , 403.0000 , 402.0000 , 401.0000 , 400.0000 , 399.0000 , 398.0000 , 373.0000 , 372.0000 , 371.0000 , 370.0000 , 369.0000 , 368.0000 , 367.0000 , 366.0000 , 365.0000 , 364.0000 , 356.0000 , 355.0000 , 354.0000 , 353.0000 , 352.0000 , 351.0000 , 350.0000 , 347.0000 , 346.0000 , 343.0000];
y = [-49.5, -47.5, -46.5, -45.5, -44.5, -42.5, -41.5, -39.5,  -4.5,  -3.5,  -1.5,   0.5,   1.5,   3.5,   4.5,   6.5,   7.5,   8.5,  25.5,  27.5,  29.5,  31.5,  32.5,  37.5,  37.5,  37.5,  37.5,  37.5];
%% 拟合测试——以三次方拟合为例——y=a0+a1*x+a2*x^2+a3*x^3;共四个系数
%% 利用自己推导的，网页中的步骤5
for i = 0:6
    eval(['x',num2str(i),' = x.^',num2str(i),';']);
    eval(['xSum',num2str(i),' = sum(x',num2str(i),'(:));']);
    eval(['yx',num2str(i),'=y.*x.^',num2str(i),';']);
    eval(['yxSum',num2str(i),'=sum(yx',num2str(i),'(:));']);
end
XX = [
    xSum0 xSum1 xSum2 xSum3
    xSum1 xSum2 xSum3 xSum4
    xSum2 xSum3 xSum4 xSum5
    xSum3 xSum4 xSum5 xSum6
    ];
YY = [
    yxSum0
    yxSum1
    yxSum2
    yxSum3
    ];
disp('自己推导：')
A = XX\YY
% A1 = inv(XX)*YY%结果与上面的左除一样
%% 利用matlab自带的拟合
disp('Matlab自带拟合：')
xishu = polyfit(x,y,3)
%% 利用网页上提供的方法 步骤6
% 测试表明，人家的公式是对的
for i = 0:3
    eval(['x_',num2str(i),'=(x'').^',num2str(i),';']);
end
XX2 = [x_0 x_1 x_2 x_3 ];
YY2 = y';
disp('步骤6：')
A2 = XX2\YY2
%% 利用网页上提供的方法  步骤7
for i = 0:3
    eval(['x_',num2str(i),'=(x'').^',num2str(i),';']);
end
XX2 = [x_0 x_1 x_2 x_3 ];
YY2 = y';
A21 = inv(XX2'*XX2)*XX2'*YY2

计算结果为：

自己推导：
A =
         -21577.7657993395
          178.557127477387
        -0.487265577436154
      0.000438644959464213
Matlab自带拟合：
xishu =
      0.000438645071094348        -0.487265702300147           178.55717395217         -21577.7715556855
步骤6：
A2 =
         -21577.7715556855
           178.55717395217
        -0.487265702300147
      0.000438645071094348
A21 =
          -21577.761617817
          178.557093705516
        -0.487265486699101
      0.000438644878352712