最小二乘法曲面拟合
http://blog.csdn.net/liumangmao1314/article/details/54179526?locationNum=12&fps=1
最近上课讲到最小二乘的基本原理，老师给出问题自己去调研一下。本科期间我第一次接触这个东西实在大学物理这个课程中，当时记得老师给了一个数值例子，然后还觉着挺简单的，当然也是不知道这个东西具体什么用的，今天自己找了不少的资料，把最小二乘法给彻底的理解一下。
最小二乘法做拟合看到的最好的一个文章就是来自博友PureSky_Memory的一篇博文。
http://blog.sina.com.cn/s/blog_648868460100hevs.html
他的文章里讲的关于最小二乘法拟合直线我个人觉着写的非常好；我这里也是借用他文中的部分思想和自己所学的知识自己写了一个用最小二乘法来曲面拟合并给出数值例子。
1、基本原理
首先是平面上的二次方程：
化简为：
对于等精度的N组数据（xi , yi），i=1,…,N;xi , yi 是精确的，假如预测值是zi 。用最小二乘法估计参数时，要求观测值 zi 的偏差加权平方和最小。用泛函误差δ表示
a0 , a1 , a2 , a3 , a4 , a5的值影响着 δ 的大小，对 δ 求最小值，自然就是对a0 , a1 , a2 ，a3 , a4 , a5 分别求偏导，实际就是一个泛函求极值问题。
具体求导：
整理后得到方程组：
2、数值例子
本实验是采用双曲抛物面作为模型，在三维直角坐标中通过自己生成x、y的坐标数据求解z值，显示得到精确解。通过最小二乘法得到所求的系数，生成拟合方程，代入10%噪音的原始数据得到拟合结果。
双曲抛物面（马鞍面）一般公式：
在程序中所使用的a=b=1。
拟合方程：
function 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
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%%
