部分内容转载自: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.使偏差绝对值最大的最小
![](https://img-my.csdn.net/uploads/201204/27/1335513871_2493.png)
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