# set dimension and sample size
p <- 10
n <- 10
# specify vector of penalty parameters
lambda <- c(2, 1)
# generate precision matrix
T1 <- matrix(0.7, p, p)
diag(T1) <- 1
T2 <- diag(rep(2, p))
# generate precision matrix
Omega <- matrix(0.4, p, p)
diag(Omega) <- 2
Sigma <- solve(Omega)
# data
Y <- mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=Sigma)
S <- cov(Y)
# find optimal penalty parameters through cross-validation
lambdaOpt <- optPenaltyPmultiT.kCVauto(Y, rep(10^(-10), 2),
rep(10^(10), 2), rep(1, 2),
targetList=list(T1=T1, T2=T2))
# unpenalized diagonal estimate
Phat <- ridgePmultiT(S, lambdaOpt, list(T1=T1, T2=T2))
Run the code above in your browser using DataLab