diff --git a/inst/examples/pinferunemjmcmc_example.R b/inst/examples/pinferunemjmcmc_example.R
index 8ca9cd5..7d5b9ba 100644
--- a/inst/examples/pinferunemjmcmc_example.R
+++ b/inst/examples/pinferunemjmcmc_example.R
@@ -10,7 +10,7 @@ formula1 <- as.formula(
 )
 
 # define the number or cpus
-M <- 4
+M <- 1L
 # define the size of the simulated samples
 NM <- 1000
 # define \k_{max} + 1 from the paper
@@ -35,7 +35,7 @@ thf <- 0.05
       interact.param = list(allow_offsprings = 3, mutation_rate = 250,
       last.mutation = 10000, max.tree.size = 5, Nvars.max = 15,
       p.allow.replace = 0.9, p.allow.tree = 0.01, p.nor = 0.9, p.and = 0.9),
-      n.models = 10000, unique = TRUE, max.cpu = 4, max.cpu.glob = 4,
+      n.models = 10000, unique = TRUE, max.cpu = M, max.cpu.glob = M,
       create.table = FALSE, create.hash = TRUE, pseudo.paral = TRUE,
       burn.in = 100, print.freq = 1000,
       advanced.param = list(
@@ -81,8 +81,34 @@ formula1 <- as.formula(
 )
 
 \donttest{
+  # Defining a custom estimator function
+  estimate.bas.glm.cpen <- function(
+    formula, data, family, prior, logn, r = 0.1, yid=1,
+    relat =c("cosi","sigmoid","tanh","atan","erf","m(")
+  ) {
+    #only poisson and binomial families are currently adopted
+    X <- model.matrix(object = formula,data = data)
+    capture.output({out <- BAS::bayesglm.fit(x = X, y = data[,yid], family=family,coefprior=prior)})
+    fmla.proc<-as.character(formula)[2:3]
+    fobserved <- fmla.proc[1]
+    fmla.proc[2]<- stringi::stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
+    fmla.proc[2]<- stringi::stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
+    sj<-2*(stringi::stri_count_fixed(str = fmla.proc[2], pattern = "*"))
+    sj<-sj+1*(stringi::stri_count_fixed(str = fmla.proc[2], pattern = "+"))
+    for(rel in relat) {
+      sj<-sj+2*(stringi::stri_count_fixed(str = fmla.proc[2], pattern = rel))
+    }
+    mlik = ((-out$deviance +2*log(r)*sum(sj)))/2
+    return(
+      list(
+        mlik = mlik, waic = -(out$deviance + 2*out$rank),
+        dic = -(out$deviance + logn*out$rank),
+        summary.fixed = list(mean = coefficients(out))
+      )
+    )
+  }
   res <- pinferunemjmcmc(
-    n.cores = 30, report.level = 0.5, num.mod.best = NM, simplify = TRUE,
+    n.cores = M, report.level = 0.5, num.mod.best = NM, simplify = TRUE,
     predict = TRUE, test.data = as.data.frame(test), link.function = g,
     runemjmcmc.params = list(
       formula = formula1, data = data.example, gen.prob = c(1, 1, 1, 1, 0),
@@ -97,7 +123,7 @@ formula1 <- as.formula(
         allow_offsprings = 4, mutation_rate = 100, last.mutation = 1000,
         max.tree.size = 6, Nvars.max = 20, p.allow.replace = 0.5,
         p.allow.tree = 0.4, p.nor = 0.3, p.and = 0.9
-      ), n.models = 7000, unique = TRUE, max.cpu = 4, max.cpu.glob = 4,
+      ), n.models = 7000, unique = TRUE, max.cpu = M, max.cpu.glob = M,
       create.table = FALSE, create.hash = TRUE, pseudo.paral = TRUE,
       burn.in = 100, print.freq = 1000,
       advanced.param = list(
@@ -107,12 +133,6 @@ formula1 <- as.formula(
     )
   )
 
-  print(auc(prob = res$predictions, y = test$X))
-  for (i in 1:M) {
-    print(auc(prob = res$threads.stats[[i]]$preds, y = test$X))
-    print(res$threads.stats[[i]]$post.populi)
-  }
-
   for (jjjj in 1:10)
   {
     resw <- as.integer(res$predictions >= 0.1 * jjjj)
diff --git a/man/pinferunemjmcmc.Rd b/man/pinferunemjmcmc.Rd
index 01be789..75b7699 100644
--- a/man/pinferunemjmcmc.Rd
+++ b/man/pinferunemjmcmc.Rd
@@ -72,7 +72,7 @@ formula1 <- as.formula(
 )
 
 # define the number or cpus
-M <- 4
+M <- 1L
 # define the size of the simulated samples
 NM <- 1000
 # define \k_{max} + 1 from the paper
@@ -97,7 +97,7 @@ thf <- 0.05
       interact.param = list(allow_offsprings = 3, mutation_rate = 250,
       last.mutation = 10000, max.tree.size = 5, Nvars.max = 15,
       p.allow.replace = 0.9, p.allow.tree = 0.01, p.nor = 0.9, p.and = 0.9),
-      n.models = 10000, unique = TRUE, max.cpu = 4, max.cpu.glob = 4,
+      n.models = 10000, unique = TRUE, max.cpu = M, max.cpu.glob = M,
       create.table = FALSE, create.hash = TRUE, pseudo.paral = TRUE,
       burn.in = 100, print.freq = 1000,
       advanced.param = list(
@@ -143,8 +143,34 @@ formula1 <- as.formula(
 )
 
 \donttest{
+  # Defining a custom estimator function
+  estimate.bas.glm.cpen <- function(
+    formula, data, family, prior, logn, r = 0.1, yid=1,
+    relat =c("cosi","sigmoid","tanh","atan","erf","m(")
+  ) {
+    #only poisson and binomial families are currently adopted
+    X <- model.matrix(object = formula,data = data)
+    capture.output({out <- BAS::bayesglm.fit(x = X, y = data[,yid], family=family,coefprior=prior)})
+    fmla.proc<-as.character(formula)[2:3]
+    fobserved <- fmla.proc[1]
+    fmla.proc[2]<- stringi::stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
+    fmla.proc[2]<- stringi::stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
+    sj<-2*(stringi::stri_count_fixed(str = fmla.proc[2], pattern = "*"))
+    sj<-sj+1*(stringi::stri_count_fixed(str = fmla.proc[2], pattern = "+"))
+    for(rel in relat) {
+      sj<-sj+2*(stringi::stri_count_fixed(str = fmla.proc[2], pattern = rel))
+    }
+    mlik = ((-out$deviance +2*log(r)*sum(sj)))/2
+    return(
+      list(
+        mlik = mlik, waic = -(out$deviance + 2*out$rank),
+        dic = -(out$deviance + logn*out$rank),
+        summary.fixed = list(mean = coefficients(out))
+      )
+    )
+  }
   res <- pinferunemjmcmc(
-    n.cores = 30, report.level = 0.5, num.mod.best = NM, simplify = TRUE,
+    n.cores = M, report.level = 0.5, num.mod.best = NM, simplify = TRUE,
     predict = TRUE, test.data = as.data.frame(test), link.function = g,
     runemjmcmc.params = list(
       formula = formula1, data = data.example, gen.prob = c(1, 1, 1, 1, 0),
@@ -159,7 +185,7 @@ formula1 <- as.formula(
         allow_offsprings = 4, mutation_rate = 100, last.mutation = 1000,
         max.tree.size = 6, Nvars.max = 20, p.allow.replace = 0.5,
         p.allow.tree = 0.4, p.nor = 0.3, p.and = 0.9
-      ), n.models = 7000, unique = TRUE, max.cpu = 4, max.cpu.glob = 4,
+      ), n.models = 7000, unique = TRUE, max.cpu = M, max.cpu.glob = M,
       create.table = FALSE, create.hash = TRUE, pseudo.paral = TRUE,
       burn.in = 100, print.freq = 1000,
       advanced.param = list(
@@ -169,12 +195,6 @@ formula1 <- as.formula(
     )
   )
 
-  print(auc(prob = res$predictions, y = test$X))
-  for (i in 1:M) {
-    print(auc(prob = res$threads.stats[[i]]$preds, y = test$X))
-    print(res$threads.stats[[i]]$post.populi)
-  }
-
   for (jjjj in 1:10)
   {
     resw <- as.integer(res$predictions >= 0.1 * jjjj)