set.seed(1)
p=100
n=1000
max.rep=1000
x=matrix(ncol=p,nrow=n)
coefi=rep(NA,p)
for (i in 1:p){
x[,i]=rnorm(n)
coefi[i]=rnorm(1)100
}
y=x%%coefi+rnorm(n)
beta=rep(0,p)
error=rep(0,max.rep)
for (j in 1:max.rep){
for (i in 1:p){
a=y-x%*%beta+beta[i]x[,i]
beta[i]=lm(a~x[,i])$coef[2]
}
error[j]=sum((y-x%%beta)^2)
}
error
plot(1:max.rep,error,ylim=c(0,1))
plot(1:10,error[1:10])