simr2 <- function(n,p,nreps) {
   r2s <- vector(length=nreps)
   for (i in 1:nreps) {
      x <- matrix(rnorm(n*p),ncol=p)
      y <- x %*% rep(1,p) + rnorm(n,sd=sqrt(p))
      r2s[i] <- getr2(x,y)
   }
   hist(r2s)
}

getr2 <- function(x,y) {
   smm <- summary(lm(y ~ x))
   smm$r.squared
}