inst/reproducevarselection/varselection.script.R

In pierrejacob/debiasedmcmc: Unbiased MCMC with couplings

# Script that sets an initial seed using rlecuyer's package
# and compute a number of estimators using that seed
#
arguments <- commandArgs(TRUE)
if (length(arguments)!=8){
  cat("needs 8 integers as arguments; JOB_ID, NRUNS, design, n, p, SNR, k, m")
  q()
}

###
# arguments <- c("1", "10", "0", 1000", "500", "2", "0", "0")
### means JOB_ID=1, NRUNS=10, design=0, n=1000, p=500, SNR=2, k=0, m=0

JOB_ID <- as.numeric(arguments[1])
NRUNS <- as.numeric(arguments[2]) # number of runs for this job
design <- as.numeric(arguments[3]) # 0 for independent design, 1 for correlated
n <- as.numeric(arguments[4])
p <- as.numeric(arguments[5])
SNR <- as.numeric(arguments[6])
k <- as.numeric(arguments[7])
m <- as.numeric(arguments[8])

# file paths
workingdirectory <- ""
resultpath <- paste0("varselection.design", design, ".SNR", SNR, ".n", n, ".p", p, ".k", k, ".m", m, ".job", JOB_ID, ".RData")
#
setwd(workingdirectory)
# load packages
library(rlecuyer)
library(unbiasedmcmc)
setmytheme()
# initial seed
.lec.SetPackageSeed(c(42, 66, 101, 123454, 7, 54321))
nstream <- 1000 # number larger than total # of processors expected to run this script
stream.names <- paste(1:nstream)
.lec.CreateStream(stream.names)
.lec.CurrentStream(paste(JOB_ID))

### beginning of job
### load data (generated by "varselection.generatedata.R")

if (design == 0){
  load(paste0("varselection.dataSNR", SNR, ".RData"))
} else {
  load(paste0("varselection.data.correlatedSNR", SNR, ".RData"))
}
Ysub <- Y[1:n]
Xsub <- X[1:n,1:p]
Y2 <- (t(Ysub) %*% Ysub)[1,1]
g <- p^3
s0 <- 100
kappa <- 2
proportion_singleflip <- 0.5

vs <- get_variableselection(Ysub,Xsub,g,kappa,s0,proportion_singleflip)

# load functions to generate meeting times
prior <- vs$prior
marginal_likelihood <- vs$marginal_likelihood
rinit <- vs$rinit
single_kernel <- vs$single_kernel
coupled_kernel <- vs$coupled_kernel
coupled_chains <- vs$coupled_chains
unbiasedestimator <- vs$unbiasedestimator

result <- list()
for (irun in 1:NRUNS){
  cat("Run #", irun, "\n")
  ue <- unbiasedestimator(single_kernel, coupled_kernel, rinit, h = function(x) x, k = k, m = m)
  result[[irun]] <- list(irun = irun, JOB_ID = JOB_ID,
                         mcmcestimator = ue$mcmcestimator, correction = ue$correction,
                         uestimator = ue$uestimator, meetingtime = ue$meetingtime, niterations = ue$iteration)
  save(result, file = resultpath)
}
.lec.CurrentStreamEnd()