起本篇题目还是比较纠结的，原因是我本意打算寻找这样一个算法：在测量数据有比较大离群点时如何估计原始模型。
上一篇曲面拟合是假设测量数据基本符合均匀分布，没有特别大的离群点的情况下，我们使用最小二乘得到了不错的拟合结果。
但是当我加入比如10个大的离群点时，该方法得到的模型就很难看了。所以我就在网上搜了一下，有没有能够剔除离群点的方法。
结果找到了如下名词：加权最小二乘、迭代最小二乘、抗差最小二乘、稳健最小二乘。
他们细节的区别我就不过分研究了，不过这些最小二乘似乎表达的是一个意思：
构造权重函数，给不同测量值不同的权重，偏差大的值权重小，偏差小的权重大，采用迭代最小二乘的方式最优化目标函数。
下面是matlab中robustfit函数权重函数，可以参考一下：
权重函数( Weight Function )
等式( Equation )
默认调节常数( Default Tuning Constant )
'andrews'
w = (abs(r)
1.339
'bisquare' (default)
w = (abs(r)<1) .* (1 - r.^2).^2
4.685
'cauchy'
w = 1 ./ (1 + r.^2)
2.385
'fair'
w = 1 ./ (1 + abs(r))
1.400
'huber'
w = 1 ./ max(1, abs(r))
1.345
'logistic'
w = tanh(r) ./ r
1.205
'ols'
传统最小二乘估计 (无权重函数)
无
'talwar'
w = 1 * (abs(r)<1)
2.795
'welsch'
w = exp(-(r.^2))
2.985
代码如下：
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; %生成待拟合点，加个噪声
Z(1:10)=Z(1:10)+400; %加入离群点
hold on;
plot3(X,Y,Z,'o');
XX=[X.^2 Y.^2 X.*Y X Y ones(100,1)];
YY=Z;
C=inv(XX'*XX)*XX'*YY; %最小二乘
z=C(1)*x.^2+C(2)*y.^2+C(3)*x.*y+C(4)*x+C(5)*y +C(6); %拟合结果
Cm=C;
mesh(x,y,z)
z=C(1)*X.^2+C(2)*Y.^2+C(3)*X.*Y+C(4)*X+C(5)*Y +C(6);
C0=C;
while 1
r = z-Z;
w = tanh(r)./r; %权重函数
W=diag(w);
C=inv(XX'*W*XX)*XX'*W*YY; %加权最小二乘
z=C(1)*X.^2+C(2)*Y.^2+C(3)*X.*Y+C(4)*X+C(5)*Y +C(6); %拟合结果
if norm(C-C0)<1e-10
break;
end
C0=C;
end
z=C(1)*x.^2+C(2)*y.^2+C(3)*x.*y+C(4)*x+C(5)*y +C(6); %拟合结果
mesh(x,y,z)
结果如下：
图中一共三个曲面，最下层是原模型，最上层是普通最小二乘拟合模型，中间层是加权最小二乘拟合模型。
可以看出，加权最小二乘效果要好一些。
参考：
如果觉得《matlab fair matlab练习程序（加权最小二乘）》对你有帮助，请点赞、收藏，并留下你的观点哦！
