#------------------------------------------------------------------------------------------- file <- url("http://www.trutschnig.net/geo_reg2.RData") load(file) head(geo_reg2) A<-geo_reg2 model<-lm(data=A,y ~ poly(x,2,raw=TRUE)) summary(model) ND<-data.frame(x=c(1.5,2.5)) p<-predict(model,new=ND) p plot(A) xg<-seq(-3,3,length=100) ND<-data.frame(x=xg) p<-predict(model,new=ND) lines(xg,p,col="red") #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 A<-data.frame(x=x,y=y) model<-lm(y ~ poly(x,2,raw=TRUE)) 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)