轉(zhuǎn)引自http://blog.csdn.net/meng4411yu/article/details/8851187
robustfit
Robust regression(穩(wěn)健回歸)
語法
b=robustfit(X,y)
b=robustfit(X,y,wfun,tune)
b=robustfit(X,y,wfun,tune,const)
[b,stats]=robustfit(...)
描述
b=robustfit(X,y)通過執(zhí)行穩(wěn)健回歸來估計線性模型y=Xb蚪战,并返回一個由回歸系數(shù)組成的向量b悍募。X是一個np預測變量矩陣,y是一個n1觀測向量条获。計算使用的方法是加上bisquare加權(quán)函數(shù)的迭代重加權(quán)最小二乘法饰剥。默認的情況下殊霞,robustfit函數(shù)把全1的一個列向量添加進X中,此列向量與b中第一個元素的常數(shù)項對應汰蓉。注意不能直接對X添加一個全1的列向量绷蹲。可以在下面的第三個描述中通過改變變量“const”來更改robustfit函數(shù)的操作顾孽。robustfit函數(shù)把X或y中的NaNs作為一個缺省值祝钢,接著把它刪除。
b=robustfit(X,y,wfun,tune)增加了一個加權(quán)函數(shù)“wfun”和常數(shù)“tune”若厚±褂ⅲ“tune”是一個調(diào)節(jié)常數(shù),其在計算權(quán)重之前被分成殘差向量测秸,如果“wfun”被指定為一個函數(shù)疤估,那么“tune”是必不可少的。權(quán)重函數(shù)“wfun”可以為下表中的任何一個權(quán)重函數(shù):
如果“tune”未被指定霎冯,那么其默認值為表中對應值铃拇。“wfun”也可以是一個把殘差向量作為輸入沈撞,并產(chǎn)生一個權(quán)重向量作為輸出的函數(shù)慷荔。通過標準誤差估計和調(diào)節(jié)參數(shù)來調(diào)整殘差〔常“wfun”可以用@(@wyfun)显晶。
b=robustfit(X,y,wfun,tune,const)增加一個“const”控制模式內(nèi)是否包含一個常數(shù)項,默認為包含(on)壹士。
[b,stats]=robustfit(...)返回一個包含一下域的STATS結(jié)構(gòu)磷雇。
'ols_s' sigma estimate (rmse) from least squares fit
'robust_s' robust estimate of sigma
'mad_s' MAD estimate of sigma; used for scaling
residuals during the iterative fitting
's' final estimate of sigma, the larger of robust_s
and a weighted average of ols_s and robust_s
'se' standard error of coefficient estimates
't' ratio of b to stats.se
'p' p-values for stats.t
'covb' estimated covariance matrix for coefficient estimates
'coeffcorr' estimated correlation of coefficient estimates
'w' vector of weights for robust fit
'h' vector of leverage values for least squares fit
'dfe' degrees of freedom for error
'R' R factor in QR decomposition of X matrix
The ROBUSTFIT function estimates the variance-covariance matrix of the coefficient estimates as V=inv(X'X)STATS.S^2. The standard errors and correlations are derived from V.
matlab例子:
x = (1:10)';
y = 10 - 2*x + randn(10,1); y(10) = 0;
使用原始最小二乘估計和穩(wěn)健回歸估計結(jié)果如下:
bls = regress(y,[ones(10,1) x])
bls =
7.2481
-1.3208
brob = robustfit(x,y)
brob =
9.1063
-1.8231
顯示結(jié)果如下:
scatter(x,y,'filled'); grid on; hold on
plot(x,bls(1)+bls(2)x,'r','LineWidth',2);
plot(x,brob(1)+brob(2)x,'g','LineWidth',2)
legend('Data','Ordinary Least Squares','Robust Regression')
一個來自網(wǎng)的例子的matlab實現(xiàn):
估計:(K>0)
matlab實現(xiàn)代碼:
[plain] view plaincopy
<embed id="ZeroClipboardMovie_1" src="http://static.blog.csdn.net/scripts/ZeroClipboard/ZeroClipboard.swf" loop="false" menu="false" quality="best" bgcolor="#ffffff" name="ZeroClipboardMovie_1" allowscriptaccess="always" allowfullscreen="false" type="application/x-shockwave-flash" pluginspage="http://www.macromedia.com/go/getflashplayer" flashvars="id=1&width=16&height=16" wmode="transparent" style="box-sizing: border-box; outline: 0px;" width="16" height="16" align="middle">
<span style="font-size:12px;">function wf=robust(x,y,k)
% find starting values using Ordinary Least Squares
w = x\y;
r = y - x*w;
scale=1;
%optional I can compute the scale using MED
% scale = median(abs(r - median(r)))/0.6745;
cvg = 1;%convergence
while (cvg > 1e-5)
r = r/scale;
wf = w; %save w
WH=wfun(r,k);%diff(rho(xc)/x)
% do weighted least-squares
yst = y.*sqrt(WH);
xst = matmul(x,sqrt(WH));
w = xst\yst;
%the new residual
r = y - x*w;
% compute the convergence
cvg = max(abs(w-wf)./abs(wf));
end;
function W=wfun(r,k)
W=zeros(length(r),1);
for i=1:length(r)
if (r(i)>=-(k-1)) && (r(i)<=k)
W(i)=1;
elseif r(i)<-(k-1)
W(i)=(k-1)4/(r(i)4);
else
W(i)=k4/(r(i)4);
end
end</span>
倦春。
另外户敬,http://blog.csdn.net/abcjennifer/article/details/7449435#(M-estimator M估計法 用于幾何模型建立)
博客中對M估計法有蠻好的解釋落剪。