D <- function(a=13,b=32,n=200,sd1=1,sde=1) {
  xi <- rnorm(n, mean=0, sd=sd1)
  ei <- rnorm(n, mean=0, sd=sde)
  yi<- a+b*xi+e
  data.frame(x1,y1)
}

statfun <- function(d) {
  coef(lm(yi~xi,data=d))
}

library(plyr)

raply(500,statfun(D()))