Brown-Mood估计、Theil估计 2、Siegel估计、线性分位回归的r语言代码

setwd("I:/研一课程/2.2回归分析/R/data")#设定当前的工作目录，重要！

#Brown-Mood估计

d<-read.table("faithful.txt",header=T)

x=d[,1];y=d[,2];

md=median(x);xx1=x[x<=md];xx2=x[x>md];yy1=y[x<=md];yy2=y[x>md]

md1=median(xx1);md2=median(xx2)

mw1=median(yy1);mw2=median(yy2)

beta=(mw2-mw1)/(md2-md1)

alpha=median(y-beta*x)

par(mfrow=c(1,1))

plot(x,y)

abline(c(alpha,beta),lty=2)

cyx=coef(lm(y~x))

abline(c(cyx))

#Theil估计 2

#install.packages("mblm")

library(mblm);

fit=mblm(y~x,repeated=F)

summary(fit)

anova(fit)

confint(fit)

plot(x,y)

abline(c(alpha,beta),lty=2)

cyx=coef(lm(y~x))

abline(c(cyx))

abline(fit,lty=3,col=2)

#Siegel估计

library(mblm);

fit2=mblm(y~x,repeated=T)

summary(fit2)

anova(fit2)

confint(fit2)

plot(x,y)

abline(c(alpha,beta),lty=2)

cyx=coef(lm(y~x))

abline(c(cyx))

abline(fit,lty=3)

abline(fit2,col=2)

#线性分位回归

#install.packages("quantreg")  #我安装不了这个包

#install.packages("SparseM")

library(quantreg);

library(SparseM);

d=read.table("D:/data/CPIGINI.txt",header=T);

x=d[,1];y=d[,2];

q05=rq(y~x,tau=0.5)

cyx=coef(lm(y~x))

par(mfrow=c(1,1));

plot(x,y)

abline(c(cyx))

abline(rq(y~x,tau=0.5),col="red")#color=red

abline(rq(y~x,tau=0.1),col="green")#color=green

abline(rq(y~x,tau=0.9),col="blue")#color=blue