--- a
+++ b/R/calc.lBF.R
@@ -0,0 +1,38 @@
+#' calculate BF
+#'
+#' @param pars A list object as returned by the function set.parameters
+#' @param data.sub A 3-dimension array (N. SNPs x N. phenotypes x 3) with state likelihoods.
+#' @param w0 An integer with column index for the null state.
+#' @param log10 Return Bayes Factor in log 10 base.
+#'
+#' @return lBF Tree Bayes Factor
+#'
+calc.lBF <- function(
+    pars = NULL,
+    data.sub = NULL,
+    w0 = 2,
+    log10 = TRUE
+) {
+
+    if( is.null(pars) | is.null(data.sub) ) {
+        stop("Missing input data.\n")
+    }
+    
+    if (ncol(data.sub)==2) w0 <- 1;
+
+    llk.full <- calc.llk.tree(pars, data.sub);
+
+    q00 <- pars$p.stay + pars$p.switch * (1-pars$pi1);
+    lq00 <- log(q00);
+    n.trans <- nrow(data.sub)-1;
+    l.p.null <- log(1-pars$pi1)+n.trans*lq00;
+    l.lk.null <- sum(data.sub[,w0]);
+    tmp <- c(llk.full, l.p.null+l.lk.null);
+    mx <- max(tmp);
+    tmp <- exp(tmp-mx);
+    lBF <- mx+log(tmp[1]-tmp[2])-l.lk.null-log(1-exp(l.p.null));
+    if (log10) lBF <- lBF/log(10);
+    
+    return(lBF);
+}
+