matlab练习程序(加权最小二乘)
起本篇题目还是比较纠结的,原因是我本意打算寻找这样一个算法:在测量数据有比较大离群点时如何估计原始模型。
上一篇曲面拟合是假设测量数据基本符合均匀分布,没有特别大的离群点的情况下,我们使用最小二乘得到了不错的拟合结果。
但是当我加入比如10个大的离群点时,该方法得到的模型就很难看了。所以我就在网上搜了一下,有没有能够剔除离群点的方法。
结果找到了如下名词:加权最小二乘、迭代最小二乘、抗差最小二乘、稳健最小二乘。
他们细节的区别我就不过分研究了,不过这些最小二乘似乎表达的是一个意思:
构造权重函数,给不同测量值不同的权重,偏差大的值权重小,偏差小的权重大,采用迭代最小二乘的方式最优化目标函数。
下面是matlab中robustfit函数权重函数,可以参考一下:
权重函数(Weight Function) | 等式(Equation) | 默认调节常数(Default Tuning Constant) |
---|---|---|
\'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=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)
结果如下:
图中一共三个曲面,最下层是原模型,最上层是普通最小二乘拟合模型,中间层是加权最小二乘拟合模型。
可以看出,加权最小二乘效果要好一些。
参考:
https://www.cnblogs.com/xiongyunqi/p/3737323.html
https://blog.csdn.net/baidu_35570545/article/details/55212241