diff --git a/R/bcf.R b/R/bcf.R
index 3a6a247c..5cb3226a 100644
--- a/R/bcf.R
+++ b/R/bcf.R
@@ -96,28 +96,29 @@
 #'
 #' @examples
 #' n <- 500
-#' x1 <- rnorm(n)
-#' x2 <- rnorm(n)
-#' x3 <- rnorm(n)
-#' x4 <- as.numeric(rbinom(n,1,0.5))
-#' x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-#' X <- cbind(x1,x2,x3,x4,x5)
-#' p <- ncol(X)
-#' g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-#' tau1 <- function(x) {rep(3,nrow(x))}
-#' tau2 <- function(x) {1+2*x[,2]*x[,4]}
-#' mu_x <- mu1(X)
-#' tau_x <- tau2(X)
-#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-#' Z <- rbinom(n,1,pi_x)
-#' E_XZ <- mu_x + Z*tau_x
-#' snr <- 4
-#' y <- E_XZ + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-#' X <- as.data.frame(X)
-#' X$x4 <- factor(X$x4, ordered = TRUE)
-#' X$x5 <- factor(X$x5, ordered = TRUE)
+#' p <- 5
+#' X <- matrix(runif(n*p), ncol = p)
+#' mu_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+#' )
+#' pi_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+#' )
+#' tau_x <- (
+#'     ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+#'     ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+#'     ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+#'     ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+#' )
+#' Z <- rbinom(n, 1, pi_x)
+#' noise_sd <- 1
+#' y <- mu_x + tau_x*Z + rnorm(n, 0, noise_sd)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
@@ -1399,28 +1400,29 @@ bcf <- function(X_train, Z_train, y_train, propensity_train = NULL, rfx_group_id
 #'
 #' @examples
 #' n <- 500
-#' x1 <- rnorm(n)
-#' x2 <- rnorm(n)
-#' x3 <- rnorm(n)
-#' x4 <- as.numeric(rbinom(n,1,0.5))
-#' x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-#' X <- cbind(x1,x2,x3,x4,x5)
-#' p <- ncol(X)
-#' g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-#' tau1 <- function(x) {rep(3,nrow(x))}
-#' tau2 <- function(x) {1+2*x[,2]*x[,4]}
-#' mu_x <- mu1(X)
-#' tau_x <- tau2(X)
-#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-#' Z <- rbinom(n,1,pi_x)
-#' E_XZ <- mu_x + Z*tau_x
-#' snr <- 4
-#' y <- E_XZ + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-#' X <- as.data.frame(X)
-#' X$x4 <- factor(X$x4, ordered = TRUE)
-#' X$x5 <- factor(X$x5, ordered = TRUE)
+#' p <- 5
+#' X <- matrix(runif(n*p), ncol = p)
+#' mu_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+#' )
+#' pi_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+#' )
+#' tau_x <- (
+#'     ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+#'     ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+#'     ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+#'     ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+#' )
+#' Z <- rbinom(n, 1, pi_x)
+#' noise_sd <- 1
+#' y <- mu_x + tau_x*Z + rnorm(n, 0, noise_sd)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
@@ -1576,22 +1578,27 @@ predict.bcfmodel <- function(object, X, Z, propensity = NULL, rfx_group_ids = NU
 #'
 #' @examples
 #' n <- 500
-#' x1 <- rnorm(n)
-#' x2 <- rnorm(n)
-#' x3 <- rnorm(n)
-#' x4 <- as.numeric(rbinom(n,1,0.5))
-#' x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-#' X <- cbind(x1,x2,x3,x4,x5)
-#' p <- ncol(X)
-#' g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-#' tau1 <- function(x) {rep(3,nrow(x))}
-#' tau2 <- function(x) {1+2*x[,2]*x[,4]}
-#' mu_x <- mu1(X)
-#' tau_x <- tau2(X)
-#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-#' Z <- rbinom(n,1,pi_x)
+#' p <- 5
+#' X <- matrix(runif(n*p), ncol = p)
+#' mu_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+#' )
+#' pi_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+#' )
+#' tau_x <- (
+#'     ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+#'     ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+#'     ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+#'     ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+#' )
+#' Z <- rbinom(n, 1, pi_x)
 #' E_XZ <- mu_x + Z*tau_x
 #' snr <- 3
 #' rfx_group_ids <- rep(c(1,2), n %/% 2)
@@ -1599,9 +1606,6 @@ predict.bcfmodel <- function(object, X, Z, propensity = NULL, rfx_group_ids = NU
 #' rfx_basis <- cbind(1, runif(n, -1, 1))
 #' rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 #' y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-#' X <- as.data.frame(X)
-#' X$x4 <- factor(X$x4, ordered = TRUE)
-#' X$x5 <- factor(X$x5, ordered = TRUE)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
@@ -1667,22 +1671,27 @@ getRandomEffectSamples.bcfmodel <- function(object, ...){
 #'
 #' @examples
 #' n <- 500
-#' x1 <- rnorm(n)
-#' x2 <- rnorm(n)
-#' x3 <- rnorm(n)
-#' x4 <- as.numeric(rbinom(n,1,0.5))
-#' x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-#' X <- cbind(x1,x2,x3,x4,x5)
-#' p <- ncol(X)
-#' g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-#' tau1 <- function(x) {rep(3,nrow(x))}
-#' tau2 <- function(x) {1+2*x[,2]*x[,4]}
-#' mu_x <- mu1(X)
-#' tau_x <- tau2(X)
-#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-#' Z <- rbinom(n,1,pi_x)
+#' p <- 5
+#' X <- matrix(runif(n*p), ncol = p)
+#' mu_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+#' )
+#' pi_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+#' )
+#' tau_x <- (
+#'     ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+#'     ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+#'     ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+#'     ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+#' )
+#' Z <- rbinom(n, 1, pi_x)
 #' E_XZ <- mu_x + Z*tau_x
 #' snr <- 3
 #' rfx_group_ids <- rep(c(1,2), n %/% 2)
@@ -1690,9 +1699,6 @@ getRandomEffectSamples.bcfmodel <- function(object, ...){
 #' rfx_basis <- cbind(1, runif(n, -1, 1))
 #' rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 #' y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-#' X <- as.data.frame(X)
-#' X$x4 <- factor(X$x4, ordered = TRUE)
-#' X$x5 <- factor(X$x5, ordered = TRUE)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
@@ -1832,22 +1838,27 @@ saveBCFModelToJson <- function(object){
 #'
 #' @examples
 #' n <- 500
-#' x1 <- rnorm(n)
-#' x2 <- rnorm(n)
-#' x3 <- rnorm(n)
-#' x4 <- as.numeric(rbinom(n,1,0.5))
-#' x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-#' X <- cbind(x1,x2,x3,x4,x5)
-#' p <- ncol(X)
-#' g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-#' tau1 <- function(x) {rep(3,nrow(x))}
-#' tau2 <- function(x) {1+2*x[,2]*x[,4]}
-#' mu_x <- mu1(X)
-#' tau_x <- tau2(X)
-#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-#' Z <- rbinom(n,1,pi_x)
+#' p <- 5
+#' X <- matrix(runif(n*p), ncol = p)
+#' mu_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+#' )
+#' pi_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+#' )
+#' tau_x <- (
+#'     ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+#'     ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+#'     ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+#'     ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+#' )
+#' Z <- rbinom(n, 1, pi_x)
 #' E_XZ <- mu_x + Z*tau_x
 #' snr <- 3
 #' rfx_group_ids <- rep(c(1,2), n %/% 2)
@@ -1855,9 +1866,6 @@ saveBCFModelToJson <- function(object){
 #' rfx_basis <- cbind(1, runif(n, -1, 1))
 #' rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 #' y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-#' X <- as.data.frame(X)
-#' X$x4 <- factor(X$x4, ordered = TRUE)
-#' X$x5 <- factor(X$x5, ordered = TRUE)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
@@ -1910,22 +1918,27 @@ saveBCFModelToJsonFile <- function(object, filename){
 #'
 #' @examples
 #' n <- 500
-#' x1 <- rnorm(n)
-#' x2 <- rnorm(n)
-#' x3 <- rnorm(n)
-#' x4 <- as.numeric(rbinom(n,1,0.5))
-#' x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-#' X <- cbind(x1,x2,x3,x4,x5)
-#' p <- ncol(X)
-#' g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-#' tau1 <- function(x) {rep(3,nrow(x))}
-#' tau2 <- function(x) {1+2*x[,2]*x[,4]}
-#' mu_x <- mu1(X)
-#' tau_x <- tau2(X)
-#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-#' Z <- rbinom(n,1,pi_x)
+#' p <- 5
+#' X <- matrix(runif(n*p), ncol = p)
+#' mu_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+#' )
+#' pi_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+#' )
+#' tau_x <- (
+#'     ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+#'     ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+#'     ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+#'     ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+#' )
+#' Z <- rbinom(n, 1, pi_x)
 #' E_XZ <- mu_x + Z*tau_x
 #' snr <- 3
 #' rfx_group_ids <- rep(c(1,2), n %/% 2)
@@ -1933,9 +1946,6 @@ saveBCFModelToJsonFile <- function(object, filename){
 #' rfx_basis <- cbind(1, runif(n, -1, 1))
 #' rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 #' y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-#' X <- as.data.frame(X)
-#' X$x4 <- factor(X$x4, ordered = TRUE)
-#' X$x5 <- factor(X$x5, ordered = TRUE)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
@@ -1990,22 +2000,27 @@ saveBCFModelToJsonString <- function(object){
 #'
 #' @examples
 #' n <- 500
-#' x1 <- rnorm(n)
-#' x2 <- rnorm(n)
-#' x3 <- rnorm(n)
-#' x4 <- as.numeric(rbinom(n,1,0.5))
-#' x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-#' X <- cbind(x1,x2,x3,x4,x5)
-#' p <- ncol(X)
-#' g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-#' tau1 <- function(x) {rep(3,nrow(x))}
-#' tau2 <- function(x) {1+2*x[,2]*x[,4]}
-#' mu_x <- mu1(X)
-#' tau_x <- tau2(X)
-#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-#' Z <- rbinom(n,1,pi_x)
+#' p <- 5
+#' X <- matrix(runif(n*p), ncol = p)
+#' mu_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+#' )
+#' pi_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+#' )
+#' tau_x <- (
+#'     ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+#'     ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+#'     ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+#'     ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+#' )
+#' Z <- rbinom(n, 1, pi_x)
 #' E_XZ <- mu_x + Z*tau_x
 #' snr <- 3
 #' rfx_group_ids <- rep(c(1,2), n %/% 2)
@@ -2013,9 +2028,6 @@ saveBCFModelToJsonString <- function(object){
 #' rfx_basis <- cbind(1, runif(n, -1, 1))
 #' rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 #' y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-#' X <- as.data.frame(X)
-#' X$x4 <- factor(X$x4, ordered = TRUE)
-#' X$x5 <- factor(X$x5, ordered = TRUE)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
@@ -2155,22 +2167,27 @@ createBCFModelFromJson <- function(json_object){
 #'
 #' @examples
 #' n <- 500
-#' x1 <- rnorm(n)
-#' x2 <- rnorm(n)
-#' x3 <- rnorm(n)
-#' x4 <- as.numeric(rbinom(n,1,0.5))
-#' x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-#' X <- cbind(x1,x2,x3,x4,x5)
-#' p <- ncol(X)
-#' g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-#' tau1 <- function(x) {rep(3,nrow(x))}
-#' tau2 <- function(x) {1+2*x[,2]*x[,4]}
-#' mu_x <- mu1(X)
-#' tau_x <- tau2(X)
-#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-#' Z <- rbinom(n,1,pi_x)
+#' p <- 5
+#' X <- matrix(runif(n*p), ncol = p)
+#' mu_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+#' )
+#' pi_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+#' )
+#' tau_x <- (
+#'     ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+#'     ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+#'     ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+#'     ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+#' )
+#' Z <- rbinom(n, 1, pi_x)
 #' E_XZ <- mu_x + Z*tau_x
 #' snr <- 3
 #' rfx_group_ids <- rep(c(1,2), n %/% 2)
@@ -2178,9 +2195,6 @@ createBCFModelFromJson <- function(json_object){
 #' rfx_basis <- cbind(1, runif(n, -1, 1))
 #' rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 #' y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-#' X <- as.data.frame(X)
-#' X$x4 <- factor(X$x4, ordered = TRUE)
-#' X$x5 <- factor(X$x5, ordered = TRUE)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
@@ -2238,22 +2252,27 @@ createBCFModelFromJsonFile <- function(json_filename){
 #'
 #' @examples
 #' n <- 500
-#' x1 <- rnorm(n)
-#' x2 <- rnorm(n)
-#' x3 <- rnorm(n)
-#' x4 <- as.numeric(rbinom(n,1,0.5))
-#' x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-#' X <- cbind(x1,x2,x3,x4,x5)
-#' p <- ncol(X)
-#' g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-#' tau1 <- function(x) {rep(3,nrow(x))}
-#' tau2 <- function(x) {1+2*x[,2]*x[,4]}
-#' mu_x <- mu1(X)
-#' tau_x <- tau2(X)
-#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-#' Z <- rbinom(n,1,pi_x)
+#' p <- 5
+#' X <- matrix(runif(n*p), ncol = p)
+#' mu_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+#' )
+#' pi_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+#' )
+#' tau_x <- (
+#'     ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+#'     ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+#'     ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+#'     ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+#' )
+#' Z <- rbinom(n, 1, pi_x)
 #' E_XZ <- mu_x + Z*tau_x
 #' snr <- 3
 #' rfx_group_ids <- rep(c(1,2), n %/% 2)
@@ -2261,9 +2280,6 @@ createBCFModelFromJsonFile <- function(json_filename){
 #' rfx_basis <- cbind(1, runif(n, -1, 1))
 #' rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 #' y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-#' X <- as.data.frame(X)
-#' X$x4 <- factor(X$x4, ordered = TRUE)
-#' X$x5 <- factor(X$x5, ordered = TRUE)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
@@ -2316,24 +2332,28 @@ createBCFModelFromJsonString <- function(json_string){
 #' @export
 #'
 #' @examples
-#' n <- 100
+#' n <- 500
 #' p <- 5
-#' x1 <- rnorm(n)
-#' x2 <- rnorm(n)
-#' x3 <- rnorm(n)
-#' x4 <- rnorm(n)
-#' x5 <- rnorm(n)
-#' X <- cbind(x1,x2,x3,x4,x5)
-#' p <- ncol(X)
-#' g <- function(x) {ifelse(x[,5] < -0.44,2,ifelse(x[,5] < 0.44,-1,4))}
-#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-#' tau1 <- function(x) {rep(3,nrow(x))}
-#' tau2 <- function(x) {1+2*x[,2]*(x[,4] > 0)}
-#' mu_x <- mu1(X)
-#' tau_x <- tau2(X)
-#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-#' Z <- rbinom(n,1,pi_x)
+#' X <- matrix(runif(n*p), ncol = p)
+#' mu_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+#' )
+#' pi_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+#' )
+#' tau_x <- (
+#'     ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+#'     ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+#'     ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+#'     ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+#' )
+#' Z <- rbinom(n, 1, pi_x)
 #' E_XZ <- mu_x + Z*tau_x
 #' snr <- 3
 #' rfx_group_ids <- rep(c(1,2), n %/% 2)
@@ -2341,9 +2361,6 @@ createBCFModelFromJsonString <- function(json_string){
 #' rfx_basis <- cbind(1, runif(n, -1, 1))
 #' rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 #' y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-#' X <- as.data.frame(X)
-#' X$x4 <- factor(X$x4, ordered = TRUE)
-#' X$x5 <- factor(X$x5, ordered = TRUE)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
@@ -2375,8 +2392,8 @@ createBCFModelFromJsonString <- function(json_string){
 #'                  rfx_group_ids_test = rfx_group_ids_test,
 #'                  rfx_basis_test = rfx_basis_test, 
 #'                  num_gfr = 10, num_burnin = 0, num_mcmc = 10)
-#' # bcf_json_list <- list(saveBCFModelToJson(bcf_model))
-#' # bcf_model_roundtrip <- createBCFModelFromCombinedJson(bcf_json_list)
+#' bcf_json_list <- list(saveBCFModelToJson(bcf_model))
+#' bcf_model_roundtrip <- createBCFModelFromCombinedJson(bcf_json_list)
 createBCFModelFromCombinedJson <- function(json_object_list){
     # Initialize the BCF model
     output <- list()
@@ -2528,24 +2545,28 @@ createBCFModelFromCombinedJson <- function(json_object_list){
 #' @export
 #'
 #' @examples
-#' n <- 100
+#' n <- 500
 #' p <- 5
-#' x1 <- rnorm(n)
-#' x2 <- rnorm(n)
-#' x3 <- rnorm(n)
-#' x4 <- rnorm(n)
-#' x5 <- rnorm(n)
-#' X <- cbind(x1,x2,x3,x4,x5)
-#' p <- ncol(X)
-#' g <- function(x) {ifelse(x[,5] < -0.44,2,ifelse(x[,5] < 0.44,-1,4))}
-#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-#' tau1 <- function(x) {rep(3,nrow(x))}
-#' tau2 <- function(x) {1+2*x[,2]*(x[,4] > 0)}
-#' mu_x <- mu1(X)
-#' tau_x <- tau2(X)
-#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-#' Z <- rbinom(n,1,pi_x)
+#' X <- matrix(runif(n*p), ncol = p)
+#' mu_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+#' )
+#' pi_x <- (
+#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+#' )
+#' tau_x <- (
+#'     ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+#'     ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+#'     ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+#'     ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+#' )
+#' Z <- rbinom(n, 1, pi_x)
 #' E_XZ <- mu_x + Z*tau_x
 #' snr <- 3
 #' rfx_group_ids <- rep(c(1,2), n %/% 2)
@@ -2553,9 +2574,6 @@ createBCFModelFromCombinedJson <- function(json_object_list){
 #' rfx_basis <- cbind(1, runif(n, -1, 1))
 #' rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 #' y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-#' X <- as.data.frame(X)
-#' X$x4 <- factor(X$x4, ordered = TRUE)
-#' X$x5 <- factor(X$x5, ordered = TRUE)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
@@ -2587,8 +2605,8 @@ createBCFModelFromCombinedJson <- function(json_object_list){
 #'                  rfx_group_ids_test = rfx_group_ids_test,
 #'                  rfx_basis_test = rfx_basis_test, 
 #'                  num_gfr = 10, num_burnin = 0, num_mcmc = 10)
-#' # bcf_json_string_list <- list(saveBCFModelToJsonString(bcf_model))
-#' # bcf_model_roundtrip <- createBCFModelFromCombinedJsonString(bcf_json_string_list)
+#' bcf_json_string_list <- list(saveBCFModelToJsonString(bcf_model))
+#' bcf_model_roundtrip <- createBCFModelFromCombinedJsonString(bcf_json_string_list)
 createBCFModelFromCombinedJsonString <- function(json_string_list){
     # Initialize the BCF model
     output <- list()
diff --git a/man/bart.Rd b/man/bart.Rd
index 3cf6f211..83e48409 100644
--- a/man/bart.Rd
+++ b/man/bart.Rd
@@ -142,7 +142,8 @@ X_test <- X[test_inds,]
 X_train <- X[train_inds,]
 y_test <- y[test_inds]
 y_train <- y[train_inds]
-bart_model <- bart(X_train = X_train, y_train = y_train, X_test = X_test)
+bart_model <- bart(X_train = X_train, y_train = y_train, X_test = X_test, 
+                   num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 plot(rowMeans(bart_model$y_hat_test), y_test, xlab = "predicted", ylab = "actual")
 abline(0,1,col="red",lty=3,lwd=3)
 }
diff --git a/man/bcf.Rd b/man/bcf.Rd
index 2fcb3993..76152c9e 100644
--- a/man/bcf.Rd
+++ b/man/bcf.Rd
@@ -147,28 +147,29 @@ Run the Bayesian Causal Forest (BCF) algorithm for regularized causal effect est
 }
 \examples{
 n <- 500
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- as.numeric(rbinom(n,1,0.5))
-x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-X <- cbind(x1,x2,x3,x4,x5)
-p <- ncol(X)
-g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-tau1 <- function(x) {rep(3,nrow(x))}
-tau2 <- function(x) {1+2*x[,2]*x[,4]}
-mu_x <- mu1(X)
-tau_x <- tau2(X)
-pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-Z <- rbinom(n,1,pi_x)
-E_XZ <- mu_x + Z*tau_x
-snr <- 4
-y <- E_XZ + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-X <- as.data.frame(X)
-X$x4 <- factor(X$x4, ordered = TRUE)
-X$x5 <- factor(X$x5, ordered = TRUE)
+p <- 5
+X <- matrix(runif(n*p), ncol = p)
+mu_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+)
+pi_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+)
+tau_x <- (
+    ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+    ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+    ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+    ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+)
+Z <- rbinom(n, 1, pi_x)
+noise_sd <- 1
+y <- mu_x + tau_x*Z + rnorm(n, 0, noise_sd)
 test_set_pct <- 0.2
 n_test <- round(test_set_pct*n)
 n_train <- n - n_test
@@ -188,7 +189,8 @@ tau_test <- tau_x[test_inds]
 tau_train <- tau_x[train_inds]
 bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train, 
                  propensity_train = pi_train, X_test = X_test, Z_test = Z_test, 
-                 propensity_test = pi_test)
+                 propensity_test = pi_test, num_gfr = 10, 
+                 num_burnin = 0, num_mcmc = 10)
 plot(rowMeans(bcf_model$mu_hat_test), mu_test, xlab = "predicted", 
      ylab = "actual", main = "Prognostic function")
 abline(0,1,col="red",lty=3,lwd=3)
diff --git a/man/createBARTModelFromCombinedJson.Rd b/man/createBARTModelFromCombinedJson.Rd
index b39421f3..35d185c3 100644
--- a/man/createBARTModelFromCombinedJson.Rd
+++ b/man/createBARTModelFromCombinedJson.Rd
@@ -38,7 +38,8 @@ X_test <- X[test_inds,]
 X_train <- X[train_inds,]
 y_test <- y[test_inds]
 y_train <- y[train_inds]
-bart_model <- bart(X_train = X_train, y_train = y_train)
+bart_model <- bart(X_train = X_train, y_train = y_train, 
+                   num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 bart_json <- list(saveBARTModelToJson(bart_model))
 bart_model_roundtrip <- createBARTModelFromCombinedJson(bart_json)
 }
diff --git a/man/createBARTModelFromCombinedJsonString.Rd b/man/createBARTModelFromCombinedJsonString.Rd
index 1ae404b3..a8470dee 100644
--- a/man/createBARTModelFromCombinedJsonString.Rd
+++ b/man/createBARTModelFromCombinedJsonString.Rd
@@ -38,7 +38,8 @@ X_test <- X[test_inds,]
 X_train <- X[train_inds,]
 y_test <- y[test_inds]
 y_train <- y[train_inds]
-bart_model <- bart(X_train = X_train, y_train = y_train)
+bart_model <- bart(X_train = X_train, y_train = y_train, 
+                   num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 bart_json_string_list <- list(saveBARTModelToJsonString(bart_model))
 bart_model_roundtrip <- createBARTModelFromCombinedJsonString(bart_json_string_list)
 }
diff --git a/man/createBARTModelFromJson.Rd b/man/createBARTModelFromJson.Rd
index bfd62f91..57686122 100644
--- a/man/createBARTModelFromJson.Rd
+++ b/man/createBARTModelFromJson.Rd
@@ -38,7 +38,8 @@ X_test <- X[test_inds,]
 X_train <- X[train_inds,]
 y_test <- y[test_inds]
 y_train <- y[train_inds]
-bart_model <- bart(X_train = X_train, y_train = y_train)
+bart_model <- bart(X_train = X_train, y_train = y_train, 
+                   num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 bart_json <- saveBARTModelToJson(bart_model)
 bart_model_roundtrip <- createBARTModelFromJson(bart_json)
 }
diff --git a/man/createBARTModelFromJsonFile.Rd b/man/createBARTModelFromJsonFile.Rd
index 5ed802cb..f714a94a 100644
--- a/man/createBARTModelFromJsonFile.Rd
+++ b/man/createBARTModelFromJsonFile.Rd
@@ -38,7 +38,8 @@ X_test <- X[test_inds,]
 X_train <- X[train_inds,]
 y_test <- y[test_inds]
 y_train <- y[train_inds]
-bart_model <- bart(X_train = X_train, y_train = y_train)
+bart_model <- bart(X_train = X_train, y_train = y_train, 
+                   num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 tmpjson <- tempfile(fileext = ".json")
 saveBARTModelToJsonFile(bart_model, file.path(tmpjson))
 bart_model_roundtrip <- createBARTModelFromJsonFile(file.path(tmpjson))
diff --git a/man/createBARTModelFromJsonString.Rd b/man/createBARTModelFromJsonString.Rd
index 50651393..9e77083a 100644
--- a/man/createBARTModelFromJsonString.Rd
+++ b/man/createBARTModelFromJsonString.Rd
@@ -38,7 +38,8 @@ X_test <- X[test_inds,]
 X_train <- X[train_inds,]
 y_test <- y[test_inds]
 y_train <- y[train_inds]
-bart_model <- bart(X_train = X_train, y_train = y_train)
+bart_model <- bart(X_train = X_train, y_train = y_train, 
+                   num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 bart_json <- saveBARTModelToJsonString(bart_model)
 bart_model_roundtrip <- createBARTModelFromJsonString(bart_json)
 y_hat_mean_roundtrip <- rowMeans(predict(bart_model_roundtrip, X_train)$y_hat)
diff --git a/man/createBCFModelFromCombinedJson.Rd b/man/createBCFModelFromCombinedJson.Rd
index b1fb9ac9..6f29569e 100644
--- a/man/createBCFModelFromCombinedJson.Rd
+++ b/man/createBCFModelFromCombinedJson.Rd
@@ -18,24 +18,28 @@ Convert a list of (in-memory) JSON strings that represent BCF models to a single
 which can be used for prediction, etc...
 }
 \examples{
-n <- 100
+n <- 500
 p <- 5
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- rnorm(n)
-x5 <- rnorm(n)
-X <- cbind(x1,x2,x3,x4,x5)
-p <- ncol(X)
-g <- function(x) {ifelse(x[,5] < -0.44,2,ifelse(x[,5] < 0.44,-1,4))}
-mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-tau1 <- function(x) {rep(3,nrow(x))}
-tau2 <- function(x) {1+2*x[,2]*(x[,4] > 0)}
-mu_x <- mu1(X)
-tau_x <- tau2(X)
-pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-Z <- rbinom(n,1,pi_x)
+X <- matrix(runif(n*p), ncol = p)
+mu_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+)
+pi_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+)
+tau_x <- (
+    ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+    ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+    ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+    ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+)
+Z <- rbinom(n, 1, pi_x)
 E_XZ <- mu_x + Z*tau_x
 snr <- 3
 rfx_group_ids <- rep(c(1,2), n \%/\% 2)
@@ -43,9 +47,6 @@ rfx_coefs <- matrix(c(-1, -1, 1, 1), nrow=2, byrow=TRUE)
 rfx_basis <- cbind(1, runif(n, -1, 1))
 rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-X <- as.data.frame(X)
-X$x4 <- factor(X$x4, ordered = TRUE)
-X$x5 <- factor(X$x5, ordered = TRUE)
 test_set_pct <- 0.2
 n_test <- round(test_set_pct*n)
 n_train <- n - n_test
@@ -76,7 +77,7 @@ bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train,
                  Z_test = Z_test, propensity_test = pi_test, 
                  rfx_group_ids_test = rfx_group_ids_test,
                  rfx_basis_test = rfx_basis_test, 
-                 num_gfr = 100, num_burnin = 0, num_mcmc = 100)
-# bcf_json_list <- list(saveBCFModelToJson(bcf_model))
-# bcf_model_roundtrip <- createBCFModelFromCombinedJson(bcf_json_list)
+                 num_gfr = 10, num_burnin = 0, num_mcmc = 10)
+bcf_json_list <- list(saveBCFModelToJson(bcf_model))
+bcf_model_roundtrip <- createBCFModelFromCombinedJson(bcf_json_list)
 }
diff --git a/man/createBCFModelFromCombinedJsonString.Rd b/man/createBCFModelFromCombinedJsonString.Rd
index 988c7346..bd7e63f2 100644
--- a/man/createBCFModelFromCombinedJsonString.Rd
+++ b/man/createBCFModelFromCombinedJsonString.Rd
@@ -18,24 +18,28 @@ Convert a list of (in-memory) JSON strings that represent BCF models to a single
 which can be used for prediction, etc...
 }
 \examples{
-n <- 100
+n <- 500
 p <- 5
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- rnorm(n)
-x5 <- rnorm(n)
-X <- cbind(x1,x2,x3,x4,x5)
-p <- ncol(X)
-g <- function(x) {ifelse(x[,5] < -0.44,2,ifelse(x[,5] < 0.44,-1,4))}
-mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-tau1 <- function(x) {rep(3,nrow(x))}
-tau2 <- function(x) {1+2*x[,2]*(x[,4] > 0)}
-mu_x <- mu1(X)
-tau_x <- tau2(X)
-pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-Z <- rbinom(n,1,pi_x)
+X <- matrix(runif(n*p), ncol = p)
+mu_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+)
+pi_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+)
+tau_x <- (
+    ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+    ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+    ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+    ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+)
+Z <- rbinom(n, 1, pi_x)
 E_XZ <- mu_x + Z*tau_x
 snr <- 3
 rfx_group_ids <- rep(c(1,2), n \%/\% 2)
@@ -43,9 +47,6 @@ rfx_coefs <- matrix(c(-1, -1, 1, 1), nrow=2, byrow=TRUE)
 rfx_basis <- cbind(1, runif(n, -1, 1))
 rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-X <- as.data.frame(X)
-X$x4 <- factor(X$x4, ordered = TRUE)
-X$x5 <- factor(X$x5, ordered = TRUE)
 test_set_pct <- 0.2
 n_test <- round(test_set_pct*n)
 n_train <- n - n_test
@@ -76,7 +77,7 @@ bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train,
                  Z_test = Z_test, propensity_test = pi_test, 
                  rfx_group_ids_test = rfx_group_ids_test,
                  rfx_basis_test = rfx_basis_test, 
-                 num_gfr = 100, num_burnin = 0, num_mcmc = 100)
-# bcf_json_string_list <- list(saveBCFModelToJsonString(bcf_model))
-# bcf_model_roundtrip <- createBCFModelFromCombinedJsonString(bcf_json_string_list)
+                 num_gfr = 10, num_burnin = 0, num_mcmc = 10)
+bcf_json_string_list <- list(saveBCFModelToJsonString(bcf_model))
+bcf_model_roundtrip <- createBCFModelFromCombinedJsonString(bcf_json_string_list)
 }
diff --git a/man/createBCFModelFromJson.Rd b/man/createBCFModelFromJson.Rd
index 2cde726a..bb20ac43 100644
--- a/man/createBCFModelFromJson.Rd
+++ b/man/createBCFModelFromJson.Rd
@@ -19,22 +19,27 @@ which can be used for prediction, etc...
 }
 \examples{
 n <- 500
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- as.numeric(rbinom(n,1,0.5))
-x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-X <- cbind(x1,x2,x3,x4,x5)
-p <- ncol(X)
-g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-tau1 <- function(x) {rep(3,nrow(x))}
-tau2 <- function(x) {1+2*x[,2]*x[,4]}
-mu_x <- mu1(X)
-tau_x <- tau2(X)
-pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-Z <- rbinom(n,1,pi_x)
+p <- 5
+X <- matrix(runif(n*p), ncol = p)
+mu_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+)
+pi_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+)
+tau_x <- (
+    ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+    ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+    ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+    ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+)
+Z <- rbinom(n, 1, pi_x)
 E_XZ <- mu_x + Z*tau_x
 snr <- 3
 rfx_group_ids <- rep(c(1,2), n \%/\% 2)
@@ -42,9 +47,6 @@ rfx_coefs <- matrix(c(-1, -1, 1, 1), nrow=2, byrow=TRUE)
 rfx_basis <- cbind(1, runif(n, -1, 1))
 rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-X <- as.data.frame(X)
-X$x4 <- factor(X$x4, ordered = TRUE)
-X$x5 <- factor(X$x5, ordered = TRUE)
 test_set_pct <- 0.2
 n_test <- round(test_set_pct*n)
 n_train <- n - n_test
@@ -77,7 +79,7 @@ bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train,
                  Z_test = Z_test, propensity_test = pi_test, 
                  rfx_group_ids_test = rfx_group_ids_test,
                  rfx_basis_test = rfx_basis_test, 
-                 num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
+                 num_gfr = 10, num_burnin = 0, num_mcmc = 10, 
                  mu_forest_params = mu_params, 
                  tau_forest_params = tau_params)
 bcf_json <- saveBCFModelToJson(bcf_model)
diff --git a/man/createBCFModelFromJsonFile.Rd b/man/createBCFModelFromJsonFile.Rd
index cb83403e..58056e11 100644
--- a/man/createBCFModelFromJsonFile.Rd
+++ b/man/createBCFModelFromJsonFile.Rd
@@ -19,22 +19,27 @@ to a BCF model object which can be used for prediction, etc...
 }
 \examples{
 n <- 500
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- as.numeric(rbinom(n,1,0.5))
-x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-X <- cbind(x1,x2,x3,x4,x5)
-p <- ncol(X)
-g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-tau1 <- function(x) {rep(3,nrow(x))}
-tau2 <- function(x) {1+2*x[,2]*x[,4]}
-mu_x <- mu1(X)
-tau_x <- tau2(X)
-pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-Z <- rbinom(n,1,pi_x)
+p <- 5
+X <- matrix(runif(n*p), ncol = p)
+mu_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+)
+pi_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+)
+tau_x <- (
+    ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+    ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+    ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+    ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+)
+Z <- rbinom(n, 1, pi_x)
 E_XZ <- mu_x + Z*tau_x
 snr <- 3
 rfx_group_ids <- rep(c(1,2), n \%/\% 2)
@@ -42,9 +47,6 @@ rfx_coefs <- matrix(c(-1, -1, 1, 1), nrow=2, byrow=TRUE)
 rfx_basis <- cbind(1, runif(n, -1, 1))
 rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-X <- as.data.frame(X)
-X$x4 <- factor(X$x4, ordered = TRUE)
-X$x5 <- factor(X$x5, ordered = TRUE)
 test_set_pct <- 0.2
 n_test <- round(test_set_pct*n)
 n_train <- n - n_test
@@ -77,7 +79,7 @@ bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train,
                  Z_test = Z_test, propensity_test = pi_test, 
                  rfx_group_ids_test = rfx_group_ids_test,
                  rfx_basis_test = rfx_basis_test, 
-                 num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
+                 num_gfr = 10, num_burnin = 0, num_mcmc = 10, 
                  mu_forest_params = mu_params, 
                  tau_forest_params = tau_params)
 # saveBCFModelToJsonFile(bcf_model, "test.json")
diff --git a/man/createBCFModelFromJsonString.Rd b/man/createBCFModelFromJsonString.Rd
index 1cd567ca..cd00edca 100644
--- a/man/createBCFModelFromJsonString.Rd
+++ b/man/createBCFModelFromJsonString.Rd
@@ -19,22 +19,27 @@ to a BCF model object which can be used for prediction, etc...
 }
 \examples{
 n <- 500
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- as.numeric(rbinom(n,1,0.5))
-x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-X <- cbind(x1,x2,x3,x4,x5)
-p <- ncol(X)
-g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-tau1 <- function(x) {rep(3,nrow(x))}
-tau2 <- function(x) {1+2*x[,2]*x[,4]}
-mu_x <- mu1(X)
-tau_x <- tau2(X)
-pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-Z <- rbinom(n,1,pi_x)
+p <- 5
+X <- matrix(runif(n*p), ncol = p)
+mu_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+)
+pi_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+)
+tau_x <- (
+    ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+    ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+    ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+    ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+)
+Z <- rbinom(n, 1, pi_x)
 E_XZ <- mu_x + Z*tau_x
 snr <- 3
 rfx_group_ids <- rep(c(1,2), n \%/\% 2)
@@ -42,9 +47,6 @@ rfx_coefs <- matrix(c(-1, -1, 1, 1), nrow=2, byrow=TRUE)
 rfx_basis <- cbind(1, runif(n, -1, 1))
 rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-X <- as.data.frame(X)
-X$x4 <- factor(X$x4, ordered = TRUE)
-X$x5 <- factor(X$x5, ordered = TRUE)
 test_set_pct <- 0.2
 n_test <- round(test_set_pct*n)
 n_train <- n - n_test
@@ -75,7 +77,7 @@ bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train,
                  Z_test = Z_test, propensity_test = pi_test, 
                  rfx_group_ids_test = rfx_group_ids_test,
                  rfx_basis_test = rfx_basis_test, 
-                 num_gfr = 100, num_burnin = 0, num_mcmc = 100)
+                 num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 # bcf_json <- saveBCFModelToJsonString(bcf_model)
 # bcf_model_roundtrip <- createBCFModelFromJsonString(bcf_json)
 }
diff --git a/man/getRandomEffectSamples.bartmodel.Rd b/man/getRandomEffectSamples.bartmodel.Rd
index 9f273732..0da1eb98 100644
--- a/man/getRandomEffectSamples.bartmodel.Rd
+++ b/man/getRandomEffectSamples.bartmodel.Rd
@@ -56,6 +56,6 @@ bart_model <- bart(X_train = X_train, y_train = y_train, X_test = X_test,
                    rfx_group_ids_test = rfx_group_ids_test, 
                    rfx_basis_train = rfx_basis_train, 
                    rfx_basis_test = rfx_basis_test, 
-                   num_gfr = 100, num_burnin = 0, num_mcmc = 100)
+                   num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 rfx_samples <- getRandomEffectSamples(bart_model)
 }
diff --git a/man/getRandomEffectSamples.bcfmodel.Rd b/man/getRandomEffectSamples.bcfmodel.Rd
index 410f44c4..5bc97da7 100644
--- a/man/getRandomEffectSamples.bcfmodel.Rd
+++ b/man/getRandomEffectSamples.bcfmodel.Rd
@@ -21,22 +21,27 @@ Extract raw sample values for each of the random effect parameter terms.
 }
 \examples{
 n <- 500
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- as.numeric(rbinom(n,1,0.5))
-x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-X <- cbind(x1,x2,x3,x4,x5)
-p <- ncol(X)
-g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-tau1 <- function(x) {rep(3,nrow(x))}
-tau2 <- function(x) {1+2*x[,2]*x[,4]}
-mu_x <- mu1(X)
-tau_x <- tau2(X)
-pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-Z <- rbinom(n,1,pi_x)
+p <- 5
+X <- matrix(runif(n*p), ncol = p)
+mu_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+)
+pi_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+)
+tau_x <- (
+    ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+    ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+    ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+    ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+)
+Z <- rbinom(n, 1, pi_x)
 E_XZ <- mu_x + Z*tau_x
 snr <- 3
 rfx_group_ids <- rep(c(1,2), n \%/\% 2)
@@ -44,9 +49,6 @@ rfx_coefs <- matrix(c(-1, -1, 1, 1), nrow=2, byrow=TRUE)
 rfx_basis <- cbind(1, runif(n, -1, 1))
 rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-X <- as.data.frame(X)
-X$x4 <- factor(X$x4, ordered = TRUE)
-X$x5 <- factor(X$x5, ordered = TRUE)
 test_set_pct <- 0.2
 n_test <- round(test_set_pct*n)
 n_train <- n - n_test
@@ -79,7 +81,7 @@ bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train,
                  Z_test = Z_test, propensity_test = pi_test, 
                  rfx_group_ids_test = rfx_group_ids_test,
                  rfx_basis_test = rfx_basis_test, 
-                 num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
+                 num_gfr = 10, num_burnin = 0, num_mcmc = 10, 
                  mu_forest_params = mu_params, 
                  tau_forest_params = tau_params)
 rfx_samples <- getRandomEffectSamples(bcf_model)
diff --git a/man/predict.bartmodel.Rd b/man/predict.bartmodel.Rd
index a6fbe26c..bd15bc42 100644
--- a/man/predict.bartmodel.Rd
+++ b/man/predict.bartmodel.Rd
@@ -56,7 +56,8 @@ X_test <- X[test_inds,]
 X_train <- X[train_inds,]
 y_test <- y[test_inds]
 y_train <- y[train_inds]
-bart_model <- bart(X_train = X_train, y_train = y_train)
+bart_model <- bart(X_train = X_train, y_train = y_train, 
+                   num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 y_hat_test <- predict(bart_model, X_test)$y_hat
 plot(rowMeans(y_hat_test), y_test, xlab = "predicted", ylab = "actual")
 abline(0,1,col="red",lty=3,lwd=3)
diff --git a/man/predict.bcfmodel.Rd b/man/predict.bcfmodel.Rd
index c0b14eb5..a535654a 100644
--- a/man/predict.bcfmodel.Rd
+++ b/man/predict.bcfmodel.Rd
@@ -39,28 +39,29 @@ Predict from a sampled BCF model on new data
 }
 \examples{
 n <- 500
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- as.numeric(rbinom(n,1,0.5))
-x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-X <- cbind(x1,x2,x3,x4,x5)
-p <- ncol(X)
-g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-tau1 <- function(x) {rep(3,nrow(x))}
-tau2 <- function(x) {1+2*x[,2]*x[,4]}
-mu_x <- mu1(X)
-tau_x <- tau2(X)
-pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-Z <- rbinom(n,1,pi_x)
-E_XZ <- mu_x + Z*tau_x
-snr <- 4
-y <- E_XZ + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-X <- as.data.frame(X)
-X$x4 <- factor(X$x4, ordered = TRUE)
-X$x5 <- factor(X$x5, ordered = TRUE)
+p <- 5
+X <- matrix(runif(n*p), ncol = p)
+mu_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+)
+pi_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+)
+tau_x <- (
+    ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+    ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+    ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+    ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+)
+Z <- rbinom(n, 1, pi_x)
+noise_sd <- 1
+y <- mu_x + tau_x*Z + rnorm(n, 0, noise_sd)
 test_set_pct <- 0.2
 n_test <- round(test_set_pct*n)
 n_train <- n - n_test
@@ -79,7 +80,8 @@ mu_train <- mu_x[train_inds]
 tau_test <- tau_x[test_inds]
 tau_train <- tau_x[train_inds]
 bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train, 
-                 propensity_train = pi_train)
+                 propensity_train = pi_train, num_gfr = 10, 
+                 num_burnin = 0, num_mcmc = 10)
 preds <- predict(bcf_model, X_test, Z_test, pi_test)
 plot(rowMeans(preds$mu_hat), mu_test, xlab = "predicted",
      ylab = "actual", main = "Prognostic function")
diff --git a/man/saveBARTModelToJson.Rd b/man/saveBARTModelToJson.Rd
index d06323d3..a617532e 100644
--- a/man/saveBARTModelToJson.Rd
+++ b/man/saveBARTModelToJson.Rd
@@ -36,6 +36,7 @@ X_test <- X[test_inds,]
 X_train <- X[train_inds,]
 y_test <- y[test_inds]
 y_train <- y[train_inds]
-bart_model <- bart(X_train = X_train, y_train = y_train)
+bart_model <- bart(X_train = X_train, y_train = y_train, 
+                   num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 bart_json <- saveBARTModelToJson(bart_model)
 }
diff --git a/man/saveBARTModelToJsonFile.Rd b/man/saveBARTModelToJsonFile.Rd
index 3b24e210..46a3110e 100644
--- a/man/saveBARTModelToJsonFile.Rd
+++ b/man/saveBARTModelToJsonFile.Rd
@@ -38,7 +38,8 @@ X_test <- X[test_inds,]
 X_train <- X[train_inds,]
 y_test <- y[test_inds]
 y_train <- y[train_inds]
-bart_model <- bart(X_train = X_train, y_train = y_train)
+bart_model <- bart(X_train = X_train, y_train = y_train, 
+                   num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 tmpjson <- tempfile(fileext = ".json")
 saveBARTModelToJsonFile(bart_model, file.path(tmpjson))
 unlink(tmpjson)
diff --git a/man/saveBARTModelToJsonString.Rd b/man/saveBARTModelToJsonString.Rd
index ec954f7b..c83f9e5d 100644
--- a/man/saveBARTModelToJsonString.Rd
+++ b/man/saveBARTModelToJsonString.Rd
@@ -36,6 +36,7 @@ X_test <- X[test_inds,]
 X_train <- X[train_inds,]
 y_test <- y[test_inds]
 y_train <- y[train_inds]
-bart_model <- bart(X_train = X_train, y_train = y_train)
+bart_model <- bart(X_train = X_train, y_train = y_train, 
+                   num_gfr = 10, num_burnin = 0, num_mcmc = 10)
 bart_json_string <- saveBARTModelToJsonString(bart_model)
 }
diff --git a/man/saveBCFModelToJson.Rd b/man/saveBCFModelToJson.Rd
index 171d0c53..300fd28a 100644
--- a/man/saveBCFModelToJson.Rd
+++ b/man/saveBCFModelToJson.Rd
@@ -17,22 +17,27 @@ Convert the persistent aspects of a BCF model to (in-memory) JSON
 }
 \examples{
 n <- 500
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- as.numeric(rbinom(n,1,0.5))
-x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-X <- cbind(x1,x2,x3,x4,x5)
-p <- ncol(X)
-g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-tau1 <- function(x) {rep(3,nrow(x))}
-tau2 <- function(x) {1+2*x[,2]*x[,4]}
-mu_x <- mu1(X)
-tau_x <- tau2(X)
-pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-Z <- rbinom(n,1,pi_x)
+p <- 5
+X <- matrix(runif(n*p), ncol = p)
+mu_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+)
+pi_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+)
+tau_x <- (
+    ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+    ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+    ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+    ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+)
+Z <- rbinom(n, 1, pi_x)
 E_XZ <- mu_x + Z*tau_x
 snr <- 3
 rfx_group_ids <- rep(c(1,2), n \%/\% 2)
@@ -40,9 +45,6 @@ rfx_coefs <- matrix(c(-1, -1, 1, 1), nrow=2, byrow=TRUE)
 rfx_basis <- cbind(1, runif(n, -1, 1))
 rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-X <- as.data.frame(X)
-X$x4 <- factor(X$x4, ordered = TRUE)
-X$x5 <- factor(X$x5, ordered = TRUE)
 test_set_pct <- 0.2
 n_test <- round(test_set_pct*n)
 n_train <- n - n_test
@@ -75,7 +77,7 @@ bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train,
                  Z_test = Z_test, propensity_test = pi_test, 
                  rfx_group_ids_test = rfx_group_ids_test,
                  rfx_basis_test = rfx_basis_test, 
-                 num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
+                 num_gfr = 10, num_burnin = 0, num_mcmc = 10, 
                  mu_forest_params = mu_params, 
                  tau_forest_params = tau_params)
 # bcf_json <- saveBCFModelToJson(bcf_model)
diff --git a/man/saveBCFModelToJsonFile.Rd b/man/saveBCFModelToJsonFile.Rd
index 2c8ea980..f57c7a7d 100644
--- a/man/saveBCFModelToJsonFile.Rd
+++ b/man/saveBCFModelToJsonFile.Rd
@@ -19,22 +19,27 @@ Convert the persistent aspects of a BCF model to (in-memory) JSON and save to a
 }
 \examples{
 n <- 500
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- as.numeric(rbinom(n,1,0.5))
-x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-X <- cbind(x1,x2,x3,x4,x5)
-p <- ncol(X)
-g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-tau1 <- function(x) {rep(3,nrow(x))}
-tau2 <- function(x) {1+2*x[,2]*x[,4]}
-mu_x <- mu1(X)
-tau_x <- tau2(X)
-pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-Z <- rbinom(n,1,pi_x)
+p <- 5
+X <- matrix(runif(n*p), ncol = p)
+mu_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+)
+pi_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+)
+tau_x <- (
+    ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+    ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+    ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+    ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+)
+Z <- rbinom(n, 1, pi_x)
 E_XZ <- mu_x + Z*tau_x
 snr <- 3
 rfx_group_ids <- rep(c(1,2), n \%/\% 2)
@@ -42,9 +47,6 @@ rfx_coefs <- matrix(c(-1, -1, 1, 1), nrow=2, byrow=TRUE)
 rfx_basis <- cbind(1, runif(n, -1, 1))
 rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-X <- as.data.frame(X)
-X$x4 <- factor(X$x4, ordered = TRUE)
-X$x5 <- factor(X$x5, ordered = TRUE)
 test_set_pct <- 0.2
 n_test <- round(test_set_pct*n)
 n_train <- n - n_test
@@ -77,7 +79,7 @@ bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train,
                  Z_test = Z_test, propensity_test = pi_test, 
                  rfx_group_ids_test = rfx_group_ids_test,
                  rfx_basis_test = rfx_basis_test, 
-                 num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
+                 num_gfr = 10, num_burnin = 0, num_mcmc = 10, 
                  mu_forest_params = mu_params, 
                  tau_forest_params = tau_params)
 # saveBCFModelToJsonFile(bcf_model, "test.json")
diff --git a/man/saveBCFModelToJsonString.Rd b/man/saveBCFModelToJsonString.Rd
index 3c0bdee1..b07a763f 100644
--- a/man/saveBCFModelToJsonString.Rd
+++ b/man/saveBCFModelToJsonString.Rd
@@ -17,22 +17,27 @@ Convert the persistent aspects of a BCF model to (in-memory) JSON string
 }
 \examples{
 n <- 500
-x1 <- rnorm(n)
-x2 <- rnorm(n)
-x3 <- rnorm(n)
-x4 <- as.numeric(rbinom(n,1,0.5))
-x5 <- as.numeric(sample(1:3,n,replace=TRUE))
-X <- cbind(x1,x2,x3,x4,x5)
-p <- ncol(X)
-g <- function(x) {ifelse(x[,5]==1,2,ifelse(x[,5]==2,-1,4))}
-mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
-mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
-tau1 <- function(x) {rep(3,nrow(x))}
-tau2 <- function(x) {1+2*x[,2]*x[,4]}
-mu_x <- mu1(X)
-tau_x <- tau2(X)
-pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
-Z <- rbinom(n,1,pi_x)
+p <- 5
+X <- matrix(runif(n*p), ncol = p)
+mu_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
+)
+pi_x <- (
+    ((0 <= X[,1]) & (0.25 > X[,1])) * (0.2) + 
+    ((0.25 <= X[,1]) & (0.5 > X[,1])) * (0.4) + 
+    ((0.5 <= X[,1]) & (0.75 > X[,1])) * (0.6) + 
+    ((0.75 <= X[,1]) & (1 > X[,1])) * (0.8)
+)
+tau_x <- (
+    ((0 <= X[,2]) & (0.25 > X[,2])) * (0.5) + 
+    ((0.25 <= X[,2]) & (0.5 > X[,2])) * (1.0) + 
+    ((0.5 <= X[,2]) & (0.75 > X[,2])) * (1.5) + 
+    ((0.75 <= X[,2]) & (1 > X[,2])) * (2.0)
+)
+Z <- rbinom(n, 1, pi_x)
 E_XZ <- mu_x + Z*tau_x
 snr <- 3
 rfx_group_ids <- rep(c(1,2), n \%/\% 2)
@@ -40,9 +45,6 @@ rfx_coefs <- matrix(c(-1, -1, 1, 1), nrow=2, byrow=TRUE)
 rfx_basis <- cbind(1, runif(n, -1, 1))
 rfx_term <- rowSums(rfx_coefs[rfx_group_ids,] * rfx_basis)
 y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
-X <- as.data.frame(X)
-X$x4 <- factor(X$x4, ordered = TRUE)
-X$x5 <- factor(X$x5, ordered = TRUE)
 test_set_pct <- 0.2
 n_test <- round(test_set_pct*n)
 n_train <- n - n_test
@@ -75,7 +77,7 @@ bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train,
                  Z_test = Z_test, propensity_test = pi_test, 
                  rfx_group_ids_test = rfx_group_ids_test,
                  rfx_basis_test = rfx_basis_test, 
-                 num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
+                 num_gfr = 10, num_burnin = 0, num_mcmc = 10, 
                  mu_forest_params = mu_params, 
                  tau_forest_params = tau_params)
 # saveBCFModelToJsonString(bcf_model)
diff --git a/man/stochtree-package.Rd b/man/stochtree-package.Rd
index 6d82b32a..81cd068c 100644
--- a/man/stochtree-package.Rd
+++ b/man/stochtree-package.Rd
@@ -4,9 +4,9 @@
 \name{stochtree-package}
 \alias{stochtree}
 \alias{stochtree-package}
-\title{stochtree: Stochastic tree Ensembles (XBART and BART) for Supervised Learning and Causal Inference}
+\title{stochtree: Stochastic Tree Ensembles (XBART and BART) for Supervised Learning and Causal Inference}
 \description{
-Flexible stochastic tree ensemble software. Robust implementations of Bayesian Additive Regression Trees (Chipman, George, McCulloch (2010) \doi{10.1214/09-AOAS285}) for supervised learning and (Bayesian Causal Forests (BCF) Hahn, Murray, Carvalho (2020) \doi{10.1214/19-BA1195}) for causal inference. Enables model serialization and parallel sampling and provides a low-level interface for custom stochastic forest samplers.
+Flexible stochastic tree ensemble software. Robust implementations of Bayesian Additive Regression Trees (BART) Chipman, George, McCulloch (2010) \doi{10.1214/09-AOAS285} for supervised learning and Bayesian Causal Forests (BCF) Hahn, Murray, Carvalho (2020) \doi{10.1214/19-BA1195} for causal inference. Enables model serialization and parallel sampling and provides a low-level interface for custom stochastic forest samplers.
 }
 \seealso{
 Useful links: