library(devtools)
library(MASS)
library(gRim)
library(igraph)
install_github("pzwiernik/golazo")
library(golazo)
We first set-up the parameters of the simulation.
d <- 60
iters <- 10 # to average out the data effects
nrhos <- 20 # for each choice of other parameters 20 different rhos between 0 and 1 will be checked.
ns <- c(d/2,d,2*d) # set possible sample sizes
rs <- c(0.5,0.8) # set possible values of the signal strength r
gammas <- c(0,0.5) # set possible values of gamma (BIC corresponds to gamma=0)
You can also embed plots, for example:
pebic.res <- array(0,c(length(ns),length(rs),length(gammas),nrhos))
gebic.res <- array(0,c(length(ns),length(rs),length(gammas),nrhos))
for (n in ns){
for (r in rs){
for (gamma in gammas){
### Generate a data generating covariance matrix
## Option 1: over a star tree (comment out if neeeded)
#Sig <- matrix(r^2,d,d)
#Sig[1,] <- Sig[,1] <- r
#diag(Sig) <- 1
## Option 2: over a chain (comment out if needed)
Sig <- r^toeplitz(0:(d-1))
# set-up the simulation
pebices <- matrix(0,iters,nrhos)
gebices <- matrix(0,iters,nrhos)
for (it in 1:iters){
#generate data
X <- mvrnorm(n,mu=rep(0,d),Sigma = Sig)
S.sim <- cov.wt(X, method ="ML")$cov
R.sim <- cov2cor(S.sim)
## for PGLASSO
L <- matrix(0,d,d)
U <- matrix(1,d,d)
diag(U) <- 0
ebicres <- EBICgolazo(R.sim,n=n,L=L,U=U,rhomin=0.01,rhomax=1,nrhos=nrhos,gamma=gamma)
pebices[it,] <- ebicres$ebic
## for GLASSO
L <- matrix(-1,d,d)
U <- matrix(1,d,d)
diag(L) <- diag(U) <- 0
ebicres2 <- EBICgolazo(R.sim,n=n,L=L,U=U,rhomin=0.01,rhomax=1,nrhos=nrhos,gamma=gamma)
gebices[it,] <- ebicres2$ebic
#res2 <- golazo(R.sim,L=ebicres2$opt.rho*L,U=ebicres2$opt.rho*U)
}
pebic.res[which(ns==n),which(rs==r),which(gammas==gamma),] <- apply(pebices,2,mean)
gebic.res[which(ns==n),which(rs==r),which(gammas==gamma),] <- apply(gebices,2,mean)
}
}
}
The other plots can be obtained by changing xx
.
xx <- c(1,1,1)
main <- bquote(gamma==.(gammas[xx[3]]) ~", "~ r==.(rs[xx[2]]) ~", "~ d==60 ~", "~ n==.(ns[xx[1]]))
pmain <- bquote("PGLasso: "~ gamma==.(gammas[xx[3]]) ~", "~ r==.(rs[xx[2]]) ~", "~ d==60 ~", "~ n==.(ns[xx[1]]))
gmain <- bquote("GLasso: "~gamma==.(gammas[xx[3]]) ~", "~ r==.(rs[xx[2]]) ~", "~ d==60 ~", "~ n==.(ns[xx[1]]))
ymin <- min(pebic.res[xx[1],xx[2],xx[3],],gebic.res[xx[1],xx[2],xx[3],])
ymax <- max(pebic.res[xx[1],xx[2],xx[3],],gebic.res[xx[1],xx[2],xx[3],])
plot(ebicres2$rhos,pebic.res[xx[1],xx[2],xx[3],],type="l",lwd=2,xlab="rho",ylab="EBIC",main=main,ylim=c(ymin,ymax),col="red")
lines(ebicres2$rhos,gebic.res[xx[1],xx[2],xx[3],],type="l",lty="dashed",lwd=2,xlab="rho",ylab="EBIC")