R语言Theil法估计回归系数

 

代码：
x=c(33,45,30,20,39,34,34,21,27,38,30)
y=c(76,103,69,50,86,85,74,58,62,88,210)
r=rle(x)
if(length(r$lengths)==length(x))    #不打结
{y=y}
if(length(r$lengths)<length(x))    #打结
{
 m=r$lengths>1
 t1=which(m==1)
 for(i in 1:length(t1))
  {
    n=r$values[t1[i]]
    t3=which(x==n)          #确定结点位置 
    x=c(x[-t3],n)  
    y=c(y[-t3],median(y[t3]))     #替换原y值
   }
}
k=c(1:length(x))
a=(y[k>1]-y[1])/(x[k>1]-x[1]) 
for(i in 2:(length(x)-1))
{
 b=(y[k>i]-y[i])/(x[k>i]-x[i])  #任意两点斜率
 a=c(a,b)
}
beta=a
beta.e=median(beta)             #斜率估计值
alpha=y-beta.e*x
alpha.e=median(alpha)           #截距估计值
beta.e
alpha.e
plot(x,y)
abline(alpha.e,beta.e)
运行结果：
> beta.e
[1] 2.111111
> alpha.e
[1] 7.75