simhet <- function(n,p,nreps,sdpow) {
   bh1s <- vector(length=nreps)
   ses <- vector(length=nreps)
   for (i in 1:nreps) {
      x <- matrix(rnorm(n*p),ncol=p)
      meany <- x %*% rep(1,p) 
      sds <- abs(meany)^sdpow
      y <- meany + rnorm(n,sd=sds)
      lmout <- lm(y ~ x)
      bh1s[i] <- coef(lmout)[2]
      ses[i] <- sqrt(vcov(lmout)[2,2])
   }
   mean(abs(bh1s - 1.0) < 1.96*ses)
}