多项式拟合polyfit之测试

最新推荐文章于 2024-05-30 09:42:44 发布

hmmwjs

最新推荐文章于 2024-05-30 09:42:44 发布

阅读量2.2k

点赞数

分类专栏：图像处理 matlab闲玩文章标签：拟合曲面

本文链接：https://blog.csdn.net/hmmwjs/article/details/79758006

版权

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

4 篇文章 0 订阅

订阅专栏

matlab闲玩

3 篇文章 0 订阅

订阅专栏

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

% 用来测试，【计算时间】
% 转换过程是对的：5对称矩阵（关于主对角线对称的方阵）-->6范德蒙行列式（第一列都是1）
% https://blog.csdn.net/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.csdn.net/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

hmmwjs

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
多项式拟合polyfit之测试

部分内容转载自：https://blog.csdn.net/JairusChan/article/details/7517773概念最小二乘法多项式曲线拟合，根据给定的m个点,并不要求这条曲线精确地经过这些点，而是曲线y=f(x)的近似曲线y= φ(x)。原理[原理部分由个人根据互联网上的资料进行总结，希望对大家能有用]&nbsp; &nbsp; &nbsp;给定数据点pi(x...
复制链接

扫一扫