MCPeSe

Example

library(huge)source("mcpese.R")########################################################### Initialize seed number:#seed = Sys.time()#seed = as.integer(seed)#seed = seed %% 100000seed=59040set.seed(seed)########################################################### Graphical model simulation:p=200n=210Model="hub"HugeData= huge.generator(n=n,d=p,graph=Model)nlambda=100########################################################### Compute Glasso solution path:L= huge(HugeData$data,nlambda=nlambda,method="glasso")########################################################### Run the A-R selection (uniform prior),ARSelect= mcpese(L,n=n,M=1000)names(ARSelect)rhos=ARSelect$rhosARSelect$accept.rate########################################################### Run the M-H selection (uniform prior),MHSelect= mcpese(L,n=n,nSteps=1000,method="M-H")names(MHSelect)rhos=MHSelect$rhosMHSelect$accept.rateplot(rhos,type="l")##########################################################

# Either use the mean of rho values...mean(rhos)#ThetaARSelect = huge(Y, lambda=mean(rhos), method="glasso")#ThetaAR = as.matrix(ThetaARSelect$icov[[1]])# ... or pick the smallest tuning parameter value from the solution path# which is larger or equal to the mean value:optARrhoIndx=ARSelect$opt.index# This is sup{i : rho[i] >= mean(rhos)}optMHrhoIndx=MHSelect$opt.indexhuge.plot(HugeData$theta)title("Ground truth")##########################################################

huge.plot(L$path[[optARrhoIndx]])title("Accept-Reject sampling")##########################################################

huge.plot(L$path[[optMHrhoIndx]])title("Metropolis-Hastings sampling")

Reference