这里用到的还是最小二乘方法，和上一次这篇文章原理差不多。
就是首先构造最小二乘函数，然后对每一个系数计算偏导，构造矩阵乘法形式，最后解方程组。
比如有一个二次曲面：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)
结果如下，深色曲面是原模型，浅色曲面是用噪声数据拟合的模型：
注：加权最小二乘可以参考我后来的这篇文章。
请发表评论