matlab練習程序（曲面擬合）

這里用到的還是最小二乘方法，和上一次這篇文章原理差不多。

就是首先構造最小二乘函數，然后對每一個系數計算偏導，構造矩陣乘法形式，最后解方程組。

比如有一個二次曲面：z=ax^2+by^2+cxy+dx+ey+f

首先構造最小二乘函數，然后計算系數偏導（我直接手寫了）：

解方程組（下圖中A矩陣后面求和符號我就沒寫了啊），然后計算C：

代碼如下：

clear all;
close all;
clc;

a=2;b=2;c=-3;d=1;e=2;f=30;   %系數         
n=1:0.2:20;
x=repmat(n,96,1);
y=repmat(n',1,96);
z=a*x.^2+b*y.^2+c*x.*y+d*x+e*y +f;      %原始模型     
surf(x,y,z)

N=100;
ind=int8(rand(N,2)*95+1);

X=x(sub2ind(size(x),ind(:,1),ind(:,2)));
Y=y(sub2ind(size(y),ind(:,1),ind(:,2)));
Z=z(sub2ind(size(z),ind(:,1),ind(:,2)))+rand(N,1)*20;       %生成待擬合點，加個噪聲

hold on;
plot3(X,Y,Z,'o');

A=[N sum(Y) sum(X) sum(X.*Y) sum(Y.^2) sum(X.^2);
   sum(Y) sum(Y.^2) sum(X.*Y) sum(X.*Y.^2) sum(Y.^3) sum(X.^2.*Y);
   sum(X) sum(X.*Y) sum(X.^2) sum(X.^2.*Y) sum(X.*Y.^2) sum(X.^3);
   sum(X.*Y) sum(X.*Y.^2) sum(X.^2.*Y) sum(X.^2.*Y.^2) sum(X.*Y.^3) sum(X.^3.*Y);
   sum(Y.^2) sum(Y.^3) sum(X.*Y.^2) sum(X.*Y.^3) sum(Y.^4) sum(X.^2.*Y.^2);
   sum(X.^2) sum(X.^2.*Y) sum(X.^3) sum(X.^3.*Y) sum(X.^2.*Y.^2) sum(X.^4)];

B=[sum(Z) sum(Z.*Y) sum(Z.*X) sum(Z.*X.*Y) sum(Z.*Y.^2) sum(Z.*X.^2)]';

C=inv(A)*B;

z=C(6)*x.^2+C(5)*y.^2+C(4)*x.*y+C(3)*x+C(2)*y +C(1);           %擬合結果

mesh(x,y,z)

結果如下，深色曲面是原模型，淺色曲面是用噪聲數據擬合的模型：

注：加權最小二乘可以參考我后來的這篇文章。