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)

結果以下：

圖中一共三個曲面，最下層是原模型，最上層是普通最小二乘擬合模型，中間層是加權最小二乘擬合模型。

能夠看出，加權最小二乘效果要好一些。

參考：

http://www.javashuo.com/article/p-kuaphqmz-ck.html

http://www.javashuo.com/article/p-dcyjkqon-nr.html