diff --git a/RBBGCMuso/R/multiSite.R b/RBBGCMuso/R/multiSite.R
index 9c5c0b0..c76ce86 100644
--- a/RBBGCMuso/R/multiSite.R
+++ b/RBBGCMuso/R/multiSite.R
@@ -33,6 +33,7 @@ multiSiteCalib <- function(measurements,
                            parameters,
                            dataVar, 
                            iterations = 100,
+                           burnin =ifelse(iterations < 3000, 3000, NULL),
                            likelihood,
                            execPath,
                            thread_prefix="thread",
@@ -247,12 +248,7 @@ multiSiteThread <- function(measuredData, parameters = NULL, startDate = NULL,
         }}
 
     print("optiMuso is randomizing the epc parameters now...",quote = FALSE)
-    if(iterations < 3000){
-        randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = 3000)
-        randVals[[2]]<- randVals[[2]][sample(1:3000,iterations),]
-    } else {
-        randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = iterations)
-    }
+    randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = iterations, burnin = burnin)
 
     origEpc <- readValuesFromFile(epcFile, randVals[[1]])
     partialResult <- matrix(ncol=length(randVals[[1]])+2*length(dataVar))
diff --git a/RBBGCMuso/R/musoRand.R b/RBBGCMuso/R/musoRand.R
index ea03a8a..416a693 100644
--- a/RBBGCMuso/R/musoRand.R
+++ b/RBBGCMuso/R/musoRand.R
@@ -8,7 +8,7 @@
 #' @importFrom limSolve xsample
 #' @export
 
-musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = NULL){
+musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = NULL, burnin = ifelse(iterations < 3000, 3000, NULL)){
     if(is.null(constrains)){
         constMatrix <- constrains
         constMatrix <- getOption("RMuso_constMatrix")[[fileType]][[as.character(getOption("RMuso_version"))]]
@@ -176,7 +176,7 @@ musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = N
         E <- do.call(rbind,lapply(Ef,function(x){x$E}))
         f <- do.call(c,lapply(Ef,function(x){x$f}))
          # browser()
-        randVal <- suppressWarnings(limSolve::xsample(G=G,H=h,E=E,F=f,iter = iterations))$X
+        randVal <- suppressWarnings(limSolve::xsample(G=G,H=h,E=E,F=f,burninlength=burnin, jiter = iterations))$X
     } else{
         Gh0<-genMat0(dependences)
         randVal <- suppressWarnings(xsample(G=Gh0$G,H=Gh0$h, iter = iterations))$X