# load libraries library(mvtnorm) library(rags2ridges) library(Matrix) library(porridge) library(ggtern) # set dimension and sample size p <- 101 n <- 25 # generate banded precision matrix betas <- seq(-2.5, 2.5, length.out=p) # generate banded precision matrix target1 <- betas/4 target2 <- rep(c(-1,1), p)[1:p] target3 <- rep(0, p) # put targets in a list targetMat <- cbind(target1, target2, target3) # set initial penalty parameter lambdas <- rep(1, ncol(targetMat)) # define storage variable lBalance <- numeric() # generate banded precision matrix diags <- list(rep(1, p), rep(0.5, p-1), rep(0.25, p-2), rep(0.1, p-3)) Sigma <- as.matrix(bandSparse(p, k = -c(0:3), diag = c(diags), symm=TRUE)) colnames(Sigma) <- rownames(Sigma) <- NULL # simulate for (it in 1:10000){ # generate data X <- rmvnorm(n, sigma=Sigma) X <- sweep(X, 2, runif(ncol(X), -10, 10), FUN="+") Y <- X %*% betas + rnorm(n) lambdasOpt <- optPenaltyGLMmultiT.kCVauto(Y, X, lambdasInit=lambdas, fold=10, stratified=TRUE, loss="loglik", targetMat, model="linear") lBalance <- rbind(lBalance, lambdasOpt) print(c(it, round(lambdasOpt/sum(lambdasOpt), 2))) } # process results fracs <- sweep(lBalance, 1, rowSums(lBalance), "/") colnames(fracs) <- c("target1", "target2", "target3") fracs <- as.data.frame(fracs) # make plot ggtern(fracs, aes(target1,target2,target3)) + geom_hex_tern(binwidth=0.02) + Tlab("-1/1") + Llab(expression(beta/4)) + Rlab(expression(0[p])) + labs(title = expression(paste("Density plot of ", (list(lambda[1],lambda[2],lambda[3]))/(lambda[1]+lambda[2]+lambda[3]), "; n=25, p=101", sep="")))