genData <- function(n,p,beta)
  {
    X <- matrix(runif(n*p),nrow=n,ncol=p)
    y <- rnorm(n,X%*%beta)
    return(list(X=X,y=y))
  }
ridge <- function(Data,n.lam=301)
  {
    require(MASS)
    lam <- c(0,exp(seq(log(1e-2),log(1e7),len=n.lam)))
    fit <- lm.ridge(y~X,Data,lambda=lam)
    ##plot(fit$lam,fit$GCV,type="o",pch=19,log="x")
    ##abline(v=fit$lam[which.min(fit$GCV)])
    return(coef(fit)[which.min(fit$GCV),])
  }
bestsub <- function(Data)
  {
    require(leaps)
    fit <- regsubsets(y~X,Data)
    b <- numeric(ncol(Data$X)+1)
    names(b) <- fit$xnames
    bb <- coef(fit,which.min(summary(fit)[["bic"]]))
    b[names(bb)] <- bb
    return(b)
  }
lasso <- function(Data)
  {
    require(glmnet)
    cvfit <- cv.glmnet(Data$X,Data$y)
    return(as.numeric(coef(cvfit,s=cvfit$lambda.min)))
  }
N <- 100
beta <- c(2,-2,rep(0,8))
p <- length(beta)
results <- array(NA,dim=c(N,p,4),
                 dimnames=list(1:N,1:p,c("Subset","Lasso","Ridge","OLS")))
for (i in 1:N)
  {
    Data <- genData(100,p,beta)
    results[i,,1] <- bestsub(Data)[-1]
    results[i,,2] <- lasso(Data)[-1]
    results[i,,3] <- ridge(Data)[-1]
    results[i,,4] <- coef(lm(y~X,Data))[-1]
    displayProgressBar(i,N)
  }
B <- apply(results,2:3,mean)-beta
V <- apply(results,2:3,var)
MSE <- B^2+V
apply(MSE,2,sum)

par(mfrow=c(1,2))
ylim <- range(results)
boxplot(results[,,1],col="gray",ylim=ylim,main="Subset")
boxplot(results[,,2],col="gray",ylim=ylim,main="Lasso")