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估计
#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估计