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%%