#Datensatz 1 file <- url("http://www.trutschnig.net/geo_reg1.RData") load(file) head(geo_reg1) A<-geo_reg1 plot(A) model<-lm(data=A,y~x) model abline(model) summary(model) summary(model)$r.squared new<-data.frame(x=1.5) predict(model,newdata=new) new<-data.frame(x=c(1.25,1.5,1.75)) predict(model,newdata=new) #------------------------------------------------------------------------------- #Datensatz 2 file <- url("http://www.trutschnig.net/geo_reg2.RData") load(file) head(geo_reg2) B<-geo_reg2 plot(B) model<-lm(data=B,y~x+I(x^2)) model summary(model) summary(model)$r.squared new<-data.frame(x=seq(-3,4,length=101)) #plot using predict ynew<-predict(model,newdata=new) lines(new$x,ynew) #-------------------------------------------------------------------------------- #Überprüfung der Qualität der Schätzer mittels simulierter Daten n<-100 a<-3;b<-(-1) x<-runif(n,0,4) error<-rnorm(n,0,1) y<-a*x+b+error A<-data.frame(x=x,y=y) plot(x,y) model<-lm(data=A,y~x) model #systematisch: R<-1000 Schaetz<-data.frame(a=rep(0,R),b=rep(0,R)) n<-100 a<-1;b<-3 for(i in 1:R){ x<-runif(n,-4,4) error<-rnorm(n,0,1) y<-a*x+b+error A<-data.frame(x=x,y=y) model<-lm(data=A,y~x) Schaetz[i,1:2]<-as.numeric(model$coefficients)[2:1] } plot(Schaetz) #Was passiert für größeres oder kleineres n ? #Was passiert wenn die Fehler größere oder kleinere Varianz haben ? #---------------------------------------------------------------------------------- #analoge Überlegungen für Polynomregression #systematic order 2 plus rgl plot - install rgl package first !!! n<-1000 x<-runif(n,-3,3) error<-rnorm(n,0,1) y<-1-2*x+3*x^2+error A<-data.frame(x=x,y=y) plot(A,col="gray") model<-lm(data=A,y ~ x + I(x^2)) #model<-lm(y ~ poly(x,2,raw=TRUE)) #nice alternative new<-data.frame(x=seq(-3,3,length=101)) #plot using predict ynew<-predict(model,newdata=new) lines(new$x,ynew) ##### model summary(model) #systematic n<-100 R<-1000 E<-data.frame(a1=rep(0,R),a2=rep(0,R),b=rep(0,R)) for(i in 1:R){ x<-seq(-3,3,length=n) error<-rnorm(n,0,1) y<-1-2*x+3*x^2+error model<-lm(y ~ x + I(x^2)) E[i,3]<-as.numeric(coefficients(model))[1] E[i,1:2]<-as.numeric(coefficients(model))[2:3] } summary(E) library(rgl) open3d() plot3d(E$a1, E$a2, E$b, col="blue",xlab = "a1", ylab = "a2", zlab = "b",size=3) plot3d(-2,3,1,add=TRUE,col="red",size=10)