diff --git a/code.R b/code.R
index 4f18595..085ac76 100644
--- a/code.R
+++ b/code.R
@@ -196,3 +196,26 @@ delta <- function(trait, tree,lambda0,se,sim,thin,burn) {
   return(deltaA)
 }
 
+
+# A SINGLE DELTA ITERATION
+delta_rep <- function(i, trait, tree, lambda0, se, sim, thin, burn){
+	rtrait <- sample(trait)
+	return(delta(rtrait, tree, lambda0, se, sim, thin, burn))
+}
+
+# COMPUTE DELTA AND P-VALUE USING SEVERAL CPUS
+delta_pvalue_parallel <- function(trait, tree, lambda0=0.1, se=0.0589, sim=10000, thin=10, burn=100, reps=100, cpus=1){
+
+    library(parallel)
+
+    # Compute delta
+    deltaA <- delta(trait, tree, lambda0, se, sim, thin, burn)
+
+    # Compute pvalue
+    random_delta <- mclapply(seq(1:reps), delta_rep, trait, tree, lambda0, se, sim, thin, burn, mc.cores = cpus)
+    p_value <- sum(random_delta>deltaA)/length(random_delta)
+    
+    message(paste("delta", "\t", "p-value"))
+    message(paste(deltaA, "\t", p_value))
+    return(c(deltaA, p_value))
+}
\ No newline at end of file
diff --git a/example.R b/example.R
index 96ee1d3..7f02e77 100644
--- a/example.R
+++ b/example.R
@@ -59,3 +59,16 @@ for ( i in 1:ns) {  mtrait[i,trait[i]] <- 1 }
 tiplabels(pie=mtrait,cex=0.5)
 
 
+###############################
+# EXAMPLE USING PARALLELIZATION
+###############################
+trait <- c(2,1,3,1,1,3,1,3,2,1,1,2,2,2,2,1,1,3,1,1)
+library(parallel)
+
+pvalue_reps=100
+cpus=detectCores()
+message(paste("CPUs:", cpus))
+
+res <- deltaA <- delta_pvalue_parallel(trait, tree, lambda0, se, sim, thin, burn, pvalue_reps, cpus)
+message(paste("delta:", res[1]))
+message(paste("pvalue:", res[2]))