matlab练习程序（加权最小二乘）

起本篇题目还是比较纠结的，原因是我本意打算寻找这样一个算法：在测量数据有比较大离群点时如何估计原始模型。

上一篇曲面拟合是假设测量数据基本符合均匀分布，没有特别大的离群点的情况下，我们使用最小二乘得到了不错的拟合结果。

但是当我加入比如10个大的离群点时，该方法得到的模型就很难看了。所以我就在网上搜了一下，有没有能够剔除离群点的方法。

结果找到了如下名词：加权最小二乘、迭代最小二乘、抗差最小二乘、稳健最小二乘。

他们细节的区别我就不过分研究了，不过这些最小二乘似乎表达的是一个意思：

构造权重函数，给不同测量值不同的权重，偏差大的值权重小，偏差小的权重大，采用迭代最小二乘的方式最优化目标函数。

下面是matlab中robustfit函数权重函数，可以参考一下：

'andrews'	w = (abs(r)<pi) .* sin(r) ./ 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=;b=;c=-;d=;e=;f=;   %系数

n=:0.2:;

x=repmat(n,,);

y=repmat(n',1,96);

z=a*x.^+b*y.^+c*x.*y+d*x+e*y +f;      %原始模型

surf(x,y,z)

N=;

ind=int8(rand(N,)*+);

X=x(sub2ind(size(x),ind(:,),ind(:,)));

Y=y(sub2ind(size(y),ind(:,),ind(:,)));

Z=z(sub2ind(size(z),ind(:,),ind(:,)))+rand(N,)*;       %生成待拟合点，加个噪声

Z(:)=Z(:)+;                    %加入离群点

hold on;

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

XX=[X.^ Y.^ X.*Y X Y ones(,)];

YY=Z;

C=inv(XX'*XX)*XX'*YY;                                          %最小二乘

z=C()*x.^+C()*y.^+C()*x.*y+C()*x+C()*y +C();           %拟合结果

Cm=C;

mesh(x,y,z)

z=C()*X.^+C()*Y.^+C()*X.*Y+C()*X+C()*Y +C();

C0=C;

while

    r = z-Z;

    w = tanh(r)./r;                                             %权重函数

    W=diag(w);

    C=inv(XX'*W*XX)*XX'*W*YY;                                   %加权最小二乘

    z=C()*X.^+C()*Y.^+C()*X.*Y+C()*X+C()*Y +C();        %拟合结果

    if norm(C-C0)<1e-10

        break;

    end

    C0=C;

end

z=C()*x.^+C()*y.^+C()*x.*y+C()*x+C()*y +C();           %拟合结果

mesh(x,y,z)

结果如下：

图中一共三个曲面，最下层是原模型，最上层是普通最小二乘拟合模型，中间层是加权最小二乘拟合模型。

可以看出，加权最小二乘效果要好一些。

参考：

https://www.cnblogs.com/xiongyunqi/p/3737323.html

https://blog.csdn.net/baidu_35570545/article/details/55212241