[轉(zhuǎn)載]robustfith函數(shù)-最小二乘估計-M估計-Robust regression

轉(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ù):


圖片.png

如果“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')

image

一個來自網(wǎng)的例子的matlab實現(xiàn):

估計:
image

(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">

  1. <span style="font-size:12px;">function wf=robust(x,y,k)

  2. % find starting values using Ordinary Least Squares

  3. w = x\y;

  4. r = y - x*w;

  5. scale=1;

  6. %optional I can compute the scale using MED

  7. % scale = median(abs(r - median(r)))/0.6745;

  8. cvg = 1;%convergence

  9. while (cvg > 1e-5)

  10. r = r/scale;

  11. wf = w; %save w

  12. WH=wfun(r,k);%diff(rho(xc)/x)

  13. % do weighted least-squares

  14. yst = y.*sqrt(WH);

  15. xst = matmul(x,sqrt(WH));

  16. w = xst\yst;

  17. %the new residual

  18. r = y - x*w;

  19. % compute the convergence

  20. cvg = max(abs(w-wf)./abs(wf));

  21. end;

  22. function W=wfun(r,k)

  23. W=zeros(length(r),1);

  24. for i=1:length(r)

  25. if (r(i)>=-(k-1)) && (r(i)<=k)

  26. W(i)=1;

  27. elseif r(i)<-(k-1)

  28. W(i)=(k-1)4/(r(i)4);

  29. else

  30. W(i)=k4/(r(i)4);

  31. end

  32. end</span>

"wfun"函數(shù)為
image

,其中
image

墓卦。并且
image

倦春。

另外户敬,http://blog.csdn.net/abcjennifer/article/details/7449435#M-estimator M估計法 用于幾何模型建立

博客中對M估計法有蠻好的解釋落剪。

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個濱河市尿庐,隨后出現(xiàn)的幾起案子忠怖,更是在濱河造成了極大的恐慌,老刑警劉巖抄瑟,帶你破解...
    沈念sama閱讀 211,884評論 6 492
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件凡泣,死亡現(xiàn)場離奇詭異,居然都是意外死亡,警方通過查閱死者的電腦和手機鞋拟,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 90,347評論 3 385
  • 文/潘曉璐 我一進店門骂维,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人贺纲,你說我怎么就攤上這事航闺。” “怎么了猴誊?”我有些...
    開封第一講書人閱讀 157,435評論 0 348
  • 文/不壞的土叔 我叫張陵潦刃,是天一觀的道長。 經(jīng)常有香客問我懈叹,道長乖杠,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 56,509評論 1 284
  • 正文 為了忘掉前任澄成,我火速辦了婚禮胧洒,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘环揽。我一直安慰自己略荡,他們只是感情好,可當我...
    茶點故事閱讀 65,611評論 6 386
  • 文/花漫 我一把揭開白布歉胶。 她就那樣靜靜地躺著汛兜,像睡著了一般。 火紅的嫁衣襯著肌膚如雪通今。 梳的紋絲不亂的頭發(fā)上粥谬,一...
    開封第一講書人閱讀 49,837評論 1 290
  • 那天,我揣著相機與錄音辫塌,去河邊找鬼漏策。 笑死,一個胖子當著我的面吹牛臼氨,可吹牛的內(nèi)容都是我干的掺喻。 我是一名探鬼主播,決...
    沈念sama閱讀 38,987評論 3 408
  • 文/蒼蘭香墨 我猛地睜開眼储矩,長吁一口氣:“原來是場噩夢啊……” “哼感耙!你這毒婦竟也來了?” 一聲冷哼從身側(cè)響起持隧,我...
    開封第一講書人閱讀 37,730評論 0 267
  • 序言:老撾萬榮一對情侶失蹤即硼,失蹤者是張志新(化名)和其女友劉穎,沒想到半個月后屡拨,有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體只酥,經(jīng)...
    沈念sama閱讀 44,194評論 1 303
  • 正文 獨居荒郊野嶺守林人離奇死亡褥实,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 36,525評論 2 327
  • 正文 我和宋清朗相戀三年,在試婚紗的時候發(fā)現(xiàn)自己被綠了裂允。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片损离。...
    茶點故事閱讀 38,664評論 1 340
  • 序言:一個原本活蹦亂跳的男人離奇死亡,死狀恐怖绝编,靈堂內(nèi)的尸體忽然破棺而出草冈,到底是詐尸還是另有隱情,我是刑警寧澤瓮增,帶...
    沈念sama閱讀 34,334評論 4 330
  • 正文 年R本政府宣布怎棱,位于F島的核電站,受9級特大地震影響绷跑,放射性物質(zhì)發(fā)生泄漏拳恋。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點故事閱讀 39,944評論 3 313
  • 文/蒙蒙 一砸捏、第九天 我趴在偏房一處隱蔽的房頂上張望谬运。 院中可真熱鬧,春花似錦垦藏、人聲如沸梆暖。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,764評論 0 21
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽轰驳。三九已至,卻和暖如春弟灼,著一層夾襖步出監(jiān)牢的瞬間级解,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 31,997評論 1 266
  • 我被黑心中介騙來泰國打工田绑, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留勤哗,地道東北人。 一個月前我還...
    沈念sama閱讀 46,389評論 2 360
  • 正文 我出身青樓掩驱,卻偏偏與公主長得像芒划,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子欧穴,可洞房花燭夜當晚...
    茶點故事閱讀 43,554評論 2 349