最小二乘法曲面拟合

2017-09-10 20:45:19

标签：

最近上课讲到最⼩⼆乘的基本原理，⽼师给出问题⾃⼰去调研⼀下。本科期间我第⼀次接触这个东西实在⼤学物理这个课程中，当时记得⽼师给了⼀个数值例⼦，然后还觉着挺简单的，当然也是不知道这个东西具体什么⽤的，今天⾃⼰了不少的资料，把最⼩⼆乘法给彻底的理解⼀下。

最⼩⼆乘法做拟合看到的最好的⼀个⽂章就是来⾃博友PureSky_Memory的⼀篇博⽂。

他的⽂章⾥讲的关于最⼩⼆乘法拟合直线我个⼈觉着写的⾮常好；我这⾥也是借⽤他⽂中的部分思想和⾃⼰所学的知识⾃⼰写了⼀个⽤最⼩⼆乘法来曲⾯拟合并给出数值例⼦。

1、基本原理

⾸先是平⾯上的⼆次⽅程：

化简为：

对于等精度的N组数据（x i , y i），i=1,…,N;x i , y i 是精确的，假如预测值是z i 。⽤最⼩⼆乘法估计参数时，要求观测值z i 的偏差加权平⽅和最⼩。⽤泛函误差δ表⽰

a0 , a1 , a2 , a3 , a4 , a5的值影响着δ的⼤⼩，对δ求最⼩值，⾃然就是对a0 , a1 , a2 ，a3 , a4 , a5 分别求偏导，实际就是⼀个泛函求极值问题。

具体求导：

整理后得到⽅程组：

⽅程组形式

矩阵形式

2、数值例⼦

本实验是采⽤双曲抛物⾯作为模型，在三维直⾓坐标中通过⾃⼰⽣成x、y的坐标数据求解z值，显⽰得到精确解。通过最⼩⼆乘法得到所求的系数，⽣成拟合⽅程，代⼊10%噪⾳的原始数据得到拟合结果。

双曲抛物⾯（马鞍⾯）⼀般公式：

在程序中所使⽤的a=b=1。

拟合⽅程：

f unction least_square_method()

%% 程序说明%%

%b本程序是使⽤双曲抛物来进⾏

%验证最⼩⼆乘法，并且加⼊10%的噪⾳

%% figure 双曲抛物

x = -19:20;

y = -19:20;

[X,Y] =meshgrid(x,y);

Z = X.^2 - Y.^2;%准确解

figure

mesh(X,Y,Z);

title('准确解')

xlabel('X')

ylabel('Y')

zlabel('Z')

%% 加⼊10%噪声后进⾏拟合%%

A = zeros(6,6);%系数矩阵

b = zeros(1,6);

[row,col] =size(X);

拟合直线

for i = 1:row

for j = 1:col

A(1,1) = A(1,1) + 1;

A(1,2) = A(1,2) + X(i,j);

A(1,3) = A(1,3) + Y(i,j);

A(1,4) = A(1,4) + X(i,j)^2;

A(1,5) = A(1,5) +X(i,j)*Y(i,j);

A(1,6) = A(1,6) + Y(i,j)^2;

A(2,1) = A(2,1) + X(i,j);

A(2,2) = A(2,2) +X(i,j)*X(i,j)*X(i,j); A(2,3) = A(2,3) +Y(i,j)*X(i,j);

A(2,4) = A(2,4) + X(i,j)^2*X(i,j);

A(2,5) = A(2,5) +X(i,j)*Y(i,j)*X(i,j); A(2,6) = A(2,6) +Y(i,j)^2*X(i,j);

A(3,1) = A(3,1) + Y(i,j);

A(3,2) = A(3,2) + X(i,j)*Y(i,j);

A(3,3) = A(3,3) +Y(i,j)*Y(i,j);

A(3,4) = A(3,4) +X(i,j)^2*Y(i,j);

A(3,5) = A(3,5) +X(i,j)*Y(i,j)*Y(i,j); A(3,6) = A(3,6) +Y(i,j)^2*Y(i,j);

A(4,1) = A(4,1) + X(i,j)^2;

A(4,2) = A(4,2) +X(i,j)*X(i,j)^2;

A(4,3) = A(4,3) +Y(i,j)*X(i,j)^2;

A(4,4) = A(4,4) +X(i,j)^2*X(i,j)^2; A(4,5) = A(4,5) +X(i,j)*Y(i,j)*X(i,j)^2; A(4,6) = A(4,6) +Y(i,j)^2*X(i,j)^2;

A(5,1) = A(5,1) +X(i,j)*Y(i,j);

A(5,2) = A(5,2) +X(i,j)*X(i,j)*Y(i,j); A(5,3) = A(5,3) +Y(i,j)*X(i,j)*Y(i,j);

A(5,4) = A(5,4) + X(i,j)^2*X(i,j)*Y(i,j);

A(5,5) = A(5,5) +X(i,j)*Y(i,j)*X(i,j)*Y(i,j);

A(5,6) = A(5,6) +Y(i,j)^2*X(i,j)*Y(i,j);

A(6,1) = A(6,1) + Y(i,j)^2;

A(6,2) = A(6,2) +X(i,j)*Y(i,j)^2;

A(6,3) = A(6,3) +Y(i,j)*Y(i,j)^2;

A(6,4) = A(6,4) +X(i,j)^2*Y(i,j)^2;

A(6,5) = A(6,5) +X(i,j)*Y(i,j)*Y(i,j)^2;

A(6,6) = A(6,6) +Y(i,j)^2*Y(i,j)^2;

b(1) = b(1) + Z(i,j);

b(2) = b(2) + Z(i,j)*X(i,j);

b(3) = b(3) + Z(i,j)*Y(i,j);

b(4) = b(4) + Z(i,j)*X(i,j)^2;

b(5) = b(5) +Z(i,j)*X(i,j)*Y(i,j);

b(6) = b(6) + Z(i,j)*Y(i,j)^2;

end

a = b*inv(A);%求解系数

x = X + 0.3*rand(row,col);%加⼊10%的噪⾳

y = Y + 0.3*rand(row,col);%加⼊10%的噪⾳

z = a(1) + a(2)*x + a(3)*y + a(4)*x.^2+ a(5)*x.*y + a(6)*y.^2;%拟合解 figure

mesh(X,Y,z)

xlabel('X')

ylabel('Y')

zlabel('Z')

title('加⼊10%噪⾳拟合结果')

%% Writed by 王明⽂ 2017/9/10%%

本文发布于:2024-09-22 01:06:31，感谢您对本站的认可！

本文链接：https://www.17tex.com/tex/3/359931.html

上一篇：直线相关与回归分析

下一篇：MATLAB中利用cftool导出代码实现一张图中拟合多条平滑曲线

标签：拟合乘法双曲得到

留言与评论（共有 0 条评论）