options(repr.plot.width=5, repr.plot.height=5) mu = 0.5 sigma = 5 nsample = 100 ntrial = 1000 MSE = function(mu.hat, mu) { return(sum((mu.hat - mu)^2) / length(mu)) } alpha = seq(0,1,length=20) mse = numeric(length(alpha)) bias = (1 - alpha) * mu variance = alpha^2 * 25 / 100 for (i in 1:ntrial) { Z = rnorm(nsample) * sigma + mu for (j in 1:length(alpha)) { mse[j] = mse[j] + MSE(alpha[j] * mean(Z) * rep(1, nsample), mu * rep(1, nsample)) } } mse = mse / ntrial plot(alpha, mse, type='l', lwd=2, col='red', ylim=c(0, max(mse)), xlab=expression(paste('Shrinkage factor,', alpha)), ylab=expression(paste('MSE(', alpha, ')')), cex.lab=1.2) plot(alpha, mse, type='l', lwd=2, col='red', ylim=c(0, max(mse)), xlab=expression(paste('Shrinkage factor,', alpha)), ylab=expression(paste('MSE(', alpha, ')')), cex.lab=1.2) lines(alpha, bias^2, col='green', lwd=2) lines(alpha, variance, col='blue', lwd=2) lam = nsample / alpha - nsample plot(lam, mse, type='l', lwd=2, col='red', ylim=c(0, max(mse)), xlab=expression(paste('Penalty parameter,', lambda)), ylab=expression(paste('MSE(', lambda, ')'))) plot(lam, mse, type='l', lwd=2, col='red', ylim=c(0, max(mse)), xlab=expression(paste('Penalty parameter,', lambda)), ylab=expression(paste('MSE(', lambda, ')'))) lines(lam, bias^2, col='green', lwd=2) lines(lam, variance, col='blue', lwd=2) plot(alpha, mse, type='l', lwd=2, col='red', ylim=c(0, max(mse)), xlab=expression(paste('Shrinkage parameter ', alpha)), ylab=expression(paste('MSE(', alpha, ')'))) abline(v=mu^2/(mu^2+sigma^2/nsample), col='blue', lty=2) library(lars) data(diabetes) library(MASS) diabetes.ridge = lm.ridge(diabetes$y ~ diabetes$x, lambda=seq(0, 100, 0.5)) plot(diabetes.ridge, lwd=3) CV = function(Z, alpha, K=5) { cve = numeric(K) n = length(Z) for (i in 1:K) { g = seq(as.integer((i-1)*n/K)+1,as.integer((i*n/K))) mu.hat = mean(Z[-g]) * alpha cve[i] = sum((Z[g]-mu.hat)^2) } return(c(sum(cve) / n, sd(cve) / sqrt(n))) } alpha = seq(0.0,1,length=20) mse = numeric(length(alpha)) avg.cv = numeric(length(alpha)) for (i in 1:ntrial) { Z = rnorm(nsample) * sigma + mu for (j in 1:length(alpha)) { current_cv = CV(Z, alpha[j]) avg.cv[j] = avg.cv[j] + current_cv[1] } } avg.cv = avg.cv / ntrial plot(alpha, avg.cv, type='l', lwd=2, col='green', xlab='Shrinkage parameter, alpha', ylab='Average CV(alpha)') abline(v=mu^2/(mu^2+sigma^2/nsample), col='blue', lty=2) abline(v=min(alpha[avg.cv == min(avg.cv)]), col='red', lty=2) cv = numeric(length(alpha)) cv.sd = numeric(length(alpha)) Z = rnorm(nsample) * sigma + mu for (j in 1:length(alpha)) { current_cv = CV(Z, alpha[j]) cv[j] = current_cv[1] cv.sd[j] = current_cv[2] } cv = numeric(length(alpha)) cv.sd = numeric(length(alpha)) nsample = 1000 Z = rnorm(nsample) * sigma + mu for (j in 1:length(alpha)) { current_cv = CV(Z, alpha[j]) cv[j] = current_cv[1] cv.sd[j] = current_cv[2] } plot(alpha, cv, type='l', lwd=2, col='green', xlab='Shrinkage parameter, alpha', ylab='CV(alpha)', xlim=c(0,1)) abline(v=mu^2/(mu^2+sigma^2/nsample), col='blue', lty=2) abline(v=min(alpha[cv == min(cv)]), col='red', lty=2) par(cex.lab=1.5) plot(diabetes.ridge$lambda, diabetes.ridge$GCV, xlab='Lambda', ylab='GCV', type='l', lwd=3, col='orange') select(diabetes.ridge) Z = 2 lam = 3 beta = seq(-4, 4, length=801) value = (beta - Z)^2/2 + lam * abs(beta) plot(beta, value, type='l', lwd=2, col='red', main="Z=2, lam=3") Z = 4 lam = 3 beta = seq(-4, 4, length=801) value = (beta - Z)^2/2 + lam * abs(beta) plot(beta, value, type='l', lwd=2, col='red', main="Z=4, lam=3") Z = -4 lam = 3 beta = seq(-4, 4, length=801) value = (beta - Z)^2/2 + lam * abs(beta) plot(beta, value, type='l', lwd=2, col='red', main="Z=-4, lam=3") library(lars) data(diabetes) diabetes.lasso = lars(diabetes$x, diabetes$y, type='lasso') plot(diabetes.lasso) par(cex.lab=1.5) cv.lars(diabetes$x, diabetes$y, K=10, type='lasso') library(glmnet) G = glmnet(diabetes$x, diabetes$y) plot(G) plot(cv.glmnet(diabetes$x, diabetes$y)) cv.glmnet(diabetes$x, diabetes$y)$lambda.1se X_HIV = read.table('http://stats191.stanford.edu/data/NRTI_X.csv', header=FALSE, sep=',') Y_HIV = read.table('http://stats191.stanford.edu/data/NRTI_Y.txt', header=FALSE, sep=',') set.seed(0) Y_HIV = as.matrix(Y_HIV)[,1] X_HIV = as.matrix(X_HIV) library(glmnet) G = glmnet(X_HIV, Y_HIV) plot(G) CV = cv.glmnet(X_HIV, Y_HIV) plot(CV) beta.hat = coef(G, s=CV$lambda.1se) beta.hat # might want to use as.numeric(beta.hat) instead of a sparse vector plot(glmnet(X_HIV, Y_HIV, alpha=0.25)) plot(glmnet(X_HIV, Y_HIV, alpha=0))