library(devtools)
library(MASS)
library(gRim)
library(igraph)
install_github("pzwiernik/golazo")
library(golazo)

Set-up the parameters

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)

Run the simulations

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)
    }
  }
}

Plot the results

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")