Robust Locally Weighted Regression 魯棒局部加權回歸 -R實現

魯棒局部加權回歸

【轉載時請注明來源】：http://www.cnblogs.com/runner-ljt/

    Ljt

    作為一個初學者，水平有限，歡迎交流指正。

 

 算法參考文獻：

 （1） Robust Locally Weighted Regression and Smoothing Scatterplots (Willism_S.Cleveland)

 （2） 數據挖掘中強局部加權回歸算法實現 (虞樂,肖基毅)

 

 

 

R實現

#Robust Locally Weighted Regression 魯棒局部加權回歸

# 一元樣本值x,y ;待預測樣本點xp ;f局部加權窗口大小（一般取1/3~2/3）;d局部加權回歸階數;
#time魯棒局部加權回歸次數(一般取2就幾乎可以滿足收斂);
#step梯度下降法固定步長;error梯度下降法終止誤差;maxiter最大迭代次數
RobustLWRegression<-function(x,y,xp,f,d,time,step,error,maxiter){
  
  m<-nrow(x)
  r<-floor(f*m) #窗口內的樣本量
  xl<-abs(x-xp)
  xll<-xl[order(xl)]
  hr<-xll[r]  #h為離xp第r個最近的樣本到xp的距離
  
  #三次權值函數(幾乎在所有情況下都能夠提供充分平滑)
  xk<-(x-xp)/hr
  w<-ifelse(abs(xk)<1,(1-abs(xk^3))^3,0)
  
  #d次回歸函數
  for(i in 2:d){
    x<-cbind(x,x^i)
  }
  x<-cbind(1,x)
  n<-ncol(x)
  
  #梯度下降法（固定步長）求局部加權回歸的系數
  theta<-matrix(0,n,1) #theta 初始值都設置為0
  iter<-0
  newerror<-1   
  while((newerror>error)|(iter<maxiter)){
    iter<-iter+1
    h<-x%*%theta   
    des<-t(t(w*(h-y))%*%x)      #局部加權梯度     
    new_theta<-theta-step*des   #直接設置固定步長   
    newerror<-t(theta-new_theta)%*%(theta-new_theta)
    theta<-new_theta      
  }
  
  #time次魯棒局部加權回歸
  for(i in 1:time){
    e<-y-x%*%theta   
    s<-median(e) 
    #四次權值函數
    xb<-e/(6*s)
    R_w<-ifelse(abs(xb)<1,(1-xb^2)^2,0)
    
    #梯度下降法求魯棒加權局部回歸
    R_theta<-matrix(0,n,1) #theta 初始值都設置為0
    R_iter<-0
    R_newerror<-1
    while((R_newerror>error)|(R_iter<maxiter)){
      R_iter<-R_iter+1
      R_h<-x%*%R_theta   
      R_des<-t(t(w*R_w*(R_h-y))%*%x)    #魯棒局部加權梯度     
      R_new_theta<-R_theta-step*R_des   #直接設置固定步長   
      R_newerror<-t(R_theta-R_new_theta)%*%(R_theta-R_new_theta)
      R_theta<-R_new_theta      
    }
    theta<-R_theta   
  }
  
  
  for(i in 2:d){
     xp<-cbind(xp,xp^i)
  }
  xp<-cbind(1,xp)
  yp<-xp%*%theta
  # costfunction<-t(x%*%theta-y)%*%(x%*%theta-y)
  # result<-list(yp,theta,iter,costfunction)
  # names(result)<-c('擬合值','系數','迭代次數','誤差')
  # result
  yp
  
}

　　

實例比較 線性回歸、局部加權線性回歸和魯棒局部加權線性回歸：

>
> t(x)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,]   58   59   60   61   62   63   64   65   66    67    68    69    70    71    72
> t(y)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,]  111  115  121  123  131  130  140  136  142   145   147   151   148   151   148
> 
> lm(y~x)

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
    -50.245        2.864  

> yy<--50.245+2.864*x
> 
> plot(x,y,col='green',pch=20,xlim=c(57,73),ylim=c(109,159))
> lines(x,y,col='green')
> lines(x,yy,col='black')
> 
> g<-apply(x,1,function(xp){LWLRegression(x,y,xp,3,1e-7,100000,stepmethod=F,step=0.00001,alpha=0.25,beta=0.8)})
>
> points(x,g,col='blue',pch=20)
> lines(x,g,col='blue')
>
> gg<-apply(x,1,function(xp){RobustLWRegression(x,y,xp,0.6,2,2,0.00000001,1e-7,10000)})
>
> points(x,gg,col='red',pch=20)
> lines(x,gg,col='red')

> legend('bottomright',legend=c('散點圖','擬合直線','局部加權散點圖','魯棒局部加權散點圖'),lwd=1,col=c('green','black','blue','red'))
>