set.seed(42) n <- 100 p <- 50 library(MASS) Sigma <- matrix(NA,p,p) for (i in 1:p) for (j in 1:p) Sigma[i,j] <- 0^abs(i-j) x <- mvrnorm(n,rep(0,p),Sigma) y <- 2*x[,5] - 2*x[,10] + 1*x[,15] - 1*x[,20] + 0.5*x[,25] - 0.5*x[,30] + rnorm(n) x.new <- mvrnorm(n,rep(0,p),Sigma) y.new <- 2*x.new[,5] - 2*x.new[,10] + 1*x.new[,15] - 1*x.new[,20] + 0.5*x.new[,25] - 0.5*x.new[,30] + rnorm(n) x <- scale(x,center=TRUE,scale=FALSE) x <- scale(x,center=FALSE,scale=sqrt(colSums(x^2)/n)) y <- y - mean(y) beta.uni <- colSums(x*y)/n beta.ml <- drop(solve(t(x)%*%x) %*% t(x) %*% y) actual.lambda <- 0.1 actual.beta <- rep(0,p) actual.nz <- rep(FALSE,p) actual.resid <- y betamat <- actual.beta thresh <- function(beta) ifelse(abs(beta) < actual.lambda,0,sign(beta)*(abs(beta)-actual.lambda)) for (i in 1:1000) { print(i) old.beta <- actual.beta for (j in 1:p) { new.beta <- thresh(actual.beta[j] + sum(x[,j]*actual.resid)/n) if (new.beta != 0) { actual.beta[j] <- new.beta actual.resid <- y - drop(x %*% actual.beta) } } old.nz <- actual.nz actual.nz <- actual.beta != 0 print(sum(actual.nz != old.nz)) print(max(abs(old.beta - actual.beta))) print(actual.beta) betamat <- rbind(betamat,actual.beta) if (max(abs(old.beta - actual.beta)) < 0.001) break } matplot(betamat,type="l",xlab="iteration",ylab=expression(hat(beta))) mean((y.new - drop(x.new %*% actual.beta))^2)