source('http://myweb.uiowa.edu/pbreheny/7600/s16/notes/fun.R')
require(ncvreg)
n <- 50
p <- 100
p0 <- 6
N <- 500

# Fixed gam
SNR <- seq(0.25, 4, 0.25)
MSE <- array(NA, dim=c(N, length(SNR), 2))
pb <- txtProgressBar(1, N, style=3)
for (i in 1:N) {
  for (j in 1:length(SNR)) {
    Data <- genData(n, p, J0=p0, SNR=SNR[j])
    pData <- genData(n, p, J0=p0, SNR=SNR[j])

    fit <- with(Data, ncvreg(X, y, penalty='lasso'))
    pred <- predict(fit, pData$X)
    lam <- fit$lambda[which.min(apply(pred-pData$y, 2, crossprod))]
    MSE[i,j,1] <- crossprod(Data$beta - coef(fit, lambda=lam)[-1])

    fit <- with(Data, ncvreg(X, y))
    pred <- predict(fit, pData$X)
    lam <- fit$lambda[which.min(apply(pred-pData$y, 2, crossprod))]
    MSE[i,j,2] <- crossprod(Data$beta - coef(fit, lambda=lam)[-1])
  }
  setTxtProgressBar(pb, i)
}
save(SNR, MSE, file='sim1.RData')