diff --git a/NAMESPACE b/NAMESPACE
index 0b16688..d9d5186 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,6 +4,7 @@ S3method(as.data.frame,iprobitData)
 S3method(fitted,iprobitMod)
 S3method(iprobit,default)
 S3method(iprobit,formula)
+S3method(iprobit,ipriorKernel)
 S3method(iprobit,iprobitMod)
 S3method(logLik,iprobitMod)
 S3method(plot,iprobitData)
@@ -12,39 +13,31 @@ S3method(predict,iprobitMod)
 S3method(print,iprobitLowerBound)
 S3method(print,iprobitMod)
 S3method(print,iprobitMod_summary)
-S3method(print,iprobitPredict)
-S3method(print,iprobitPredict_quant)
 S3method(summary,iprobitMod)
 S3method(update,iprobitMod)
-export(convert_prob)
 export(expit)
 export(gen_circle)
 export(gen_mixture)
 export(gen_spiral)
 export(get_brier_score)
 export(get_brier_scores)
-export(get_coef_se_mult)
 export(get_error_rate)
 export(get_error_rates)
-export(get_kernel)
+export(get_hyperparam)
+export(get_intercept)
 export(get_lbs)
-export(get_one.lam)
 export(iplot_dec_bound)
 export(iplot_error)
 export(iplot_fitted)
 export(iplot_lb)
 export(iplot_predict)
 export(iprobit)
-export(iprobit_bin)
-export(iprobit_mult)
 export(iprobit_parallel)
 export(is.iprobitData)
+export(is.iprobitMod)
 export(is.iprobitMod_bin)
 export(is.iprobitMod_mult)
 export(logit)
-export(predict_quant)
-export(quantile_prob)
-export(sample_prob_mult)
 import(ggplot2)
 import(iprior)
 importFrom(stats,coef)
diff --git a/R/Accessor_functions.R b/R/Accessor_functions.R
new file mode 100644
index 0000000..1618fce
--- /dev/null
+++ b/R/Accessor_functions.R
@@ -0,0 +1,141 @@
+################################################################################
+#
+#   iprobit: Binary and Multinomial Probit Regression with I-priors
+#   Copyright (C) 2017  Haziq Jamil
+#
+#   This program is free software: you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation, either version 3 of the License, or
+#   (at your option) any later version.
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License
+#   along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+#' Accessor functions for \code{ipriorMod} objects.
+#'
+#' @param object An \code{ipriorMod} object.
+#'
+#' @name Accessors
+NULL
+
+#' @describeIn Accessors Obtain all of the hyperparameters.
+#' @export
+get_hyperparam <- function(object) {
+  check_and_get_iprobitMod(object)
+  object$param.full
+}
+
+#' @describeIn Accessors Obtain the intercept.
+#' @export
+get_intercept <- function(object, by.class = FALSE) {
+  res <- get_hyperparam(object)
+  res <- res[grep("Intercept", rownames(res)), , drop = FALSE]
+  if (!isTRUE(by.class)) {
+    res <- c(res)
+    if (is.common.intercept(object)) {
+      res <- res[1]
+      names(res) <- "Intercept"
+    } else {
+      names(res) <- get_names(object, "intercept")
+    }
+  }
+  res
+}
+
+get_alpha <- get_intercept
+
+get_lambda <- function(object, by.class = FALSE) {
+  res <- get_hyperparam(object)
+  res <- res[grep("lambda", rownames(res)), , drop = FALSE]
+  if (!isTRUE(by.class)) {
+    res <- c(res)
+    if (is.common.intercept(object)) {
+      res <- res[1]
+      names(res) <- get_names(object, "lambda", FALSE)
+    } else {
+      names(res) <- get_names(object, "lambda", TRUE)
+    }
+  }
+  res
+}
+
+get_sd <- function(object) {
+  setNames(object$param.summ$S.D., rownames(object$param.summ))
+}
+
+get_sd_alpha <- function(object) {
+  res <- get_sd(object)
+  res[grep("Intercept", names(res))]
+}
+
+get_sd_lambda <- function(object) {
+  res <- get_sd(object)
+  res[grep("lambda", names(res))]
+}
+
+
+
+#' @export
+get_error_rate <- function(x) x$fitted.values$train.error
+
+#' @export
+get_error_rates <- function(x) {
+  res <- x$error
+  names(res) <- seq_along(res)
+  res
+}
+
+#' @export
+get_brier_score <- function(x) x$fitted.values$brier.score
+
+#' @export
+get_brier_scores <- function(x) {
+  res <- x$brier
+  names(res) <- seq_along(res)
+  res
+}
+
+get_m <- function(object) {
+  if (is.iprobitMod(object)) object <- object$ipriorKernel
+  length(object$y.levels)
+}
+
+
+#' Extract the variational lower bound
+#'
+#' @param object An object of class \code{ipriorProbit}.
+#' @param ... This is not used here.
+#'
+#' @return The variational lower bound.
+#' @export
+logLik.iprobitMod <- function(object, ...) {
+  lb <- object$lower.bound[!is.na(object$lower.bound)]
+  lb <- lb[length(lb)]
+  class(lb) <- "iprobitLowerBound"
+  lb
+}
+
+#' @export
+get_lbs <- function(x) x$lower.bound
+
+#' @export
+print.iprobitLowerBound <- function(x, ...) {
+  cat("Lower bound =", x)
+}
+
+get_theta <- function(object) object$theta
+
+get_w <- function(object) object$w
+
+get_y <- function(object) {
+  res <- factor(object$ipriorKernel$y)
+  levels(res) <- object$ipriorKernel$y.levels
+  res
+}
diff --git a/R/Plots.R b/R/Plots.R
index 5d84802..2bc191e 100644
--- a/R/Plots.R
+++ b/R/Plots.R
@@ -18,6 +18,13 @@
 #
 ################################################################################
 
+# iplot_fitted()
+# iplot_dec_bound()
+# iplot_predict()
+# iplot_lb()
+# iplot_error()
+# iplot_lb_and_error()
+
 #' @export
 plot.iprobitMod <- function(x, niter.plot = NULL, levels = NULL, ...) {
   iplot_fitted(x)
@@ -27,10 +34,9 @@ plot.iprobitMod <- function(x, niter.plot = NULL, levels = NULL, ...) {
 iplot_fitted <- function(object) {
   list2env(object, environment())
   list2env(ipriorKernel, environment())
-  list2env(model, environment())
 
   probs <- fitted(object)$prob
-  if (isNystrom(ipriorKernel)) probs <- probs[order(Nystrom$Nys.samp), ]
+  # if (isNystrom(ipriorKernel)) probs <- probs[order(Nystrom$Nys.samp), ]
   df.plot <- data.frame(probs, i = 1:n)
   colnames(df.plot) <- c(y.levels, "i")
   df.plot <- reshape2::melt(df.plot, id.vars = "i")
@@ -44,7 +50,7 @@ iplot_fitted <- function(object) {
 }
 
 #' @export
-iplot_lb <- function(x, niter.plot = NULL, lab.pos = c("up", "down"), ...) {
+iplot_lb <- function(x, niter.plot, lab.pos = c("up", "down"), ...) {
   lae.check <- FALSE
   extra.opt <- list(...)
   if (isTRUE(extra.opt$lb.and.error)) {
@@ -61,7 +67,7 @@ iplot_lb <- function(x, niter.plot = NULL, lab.pos = c("up", "down"), ...) {
 
   lb.original <- x$lower.bound
   if (missing(niter.plot)) niter.plot <- seq_along(lb.original)
-  else if (length(niter.plot) == 1) niter.plot <- c(1, niter.plot)
+  else if (length(niter.plot) == 1) niter.plot <- seq_len(niter.plot)
   lb <- lb.original[niter.plot]
   plot.df <- data.frame(Iteration = niter.plot, lb = lb)
   time.per.iter <- x$time$time / x$niter
@@ -191,7 +197,7 @@ iplot_dec_bound <- function(object, X.var = c(1, 2), col = "grey35", size = 0.8,
     geom_point(data = points.df, aes(X1, X2, col = Class)) +
     geom_contour(data = plot.df, aes(X1, X2, z = value, group = variable,
                                      size = "Decision\nboundary"),
-                 bins = 2, col = col, ...) +
+                 binwidth = 0.5 + 1e-12, col = col, ...) +
     coord_cartesian(xlim = mm[1, ], ylim = mm[2, ]) +
     scale_colour_manual(values = c(iprior::gg_col_hue(m), "grey30")) +
     scale_size_manual(values = size, name = NULL) +
@@ -230,18 +236,44 @@ iplot_predict <- function(object, X.var = c(1, 2), grid.len = 50,
 }
 
 iplot_predict_bin <- function(plot.df, points.df, x, y, m, dec.bound) {
-  ggplot() +
-    geom_raster(data = plot.df, aes(X1, X2, fill = class2), alpha = 0.5) +
+  p <- ggplot() +
+    geom_raster(data = plot.df, aes(X1, X2, fill = `2`), alpha = 0.5) +
     scale_fill_gradient(low = "#F8766D", high = "#00BFC4", limits = c(0, 1)) +
     # annotate(geom = "raster", x = plot.df[, 1], y = plot.df[, 2],
     #          alpha = 0.6 * plot.df[, 3], fill = "#F8766D") +
     # annotate(geom = "raster", x = plot.df[, 1], y = plot.df[, 2],
     #          alpha = 0.6 * plot.df[, 4], fill = "#00BFC4") +
-    geom_point(data = points.df, aes(X1, X2, col = Class)) +
-               # col = "black", shape = 21, stroke = 0.8) +
-    coord_cartesian(xlim = x, ylim = y) +
-    guides(fill = FALSE) +
     theme_bw()
+
+  # Add decision boundary ------------------------------------------------------
+  if (isTRUE(dec.bound)) {
+    p <- p +
+      geom_contour(data = plot.df, aes(X1, X2, z = `2`,
+                                       size = "Decision\nboundary"),
+                   binwidth = 0.5 + 1e-12, col = "grey35",
+                   linetype = "dashed") +
+      scale_size_manual(values = 0.8, name = NULL) +
+      scale_color_discrete(name = " Class") +
+      guides(fill = FALSE,
+             size = guide_legend(override.aes = list(linetype = 2, col = "grey35")),
+             col = guide_legend(order = 1, override.aes = list(linetype = 0))) +
+      theme(legend.key.width = unit(2, "line"))
+  } else {
+    p <- p + guides(fill = FALSE)
+  }
+
+  # Add points -----------------------------------------------------------------
+  p <- p +
+    ggrepel::geom_label_repel(data = points.df, segment.colour = "grey25",
+                              box.padding = 0.9, show.legend = FALSE,
+                              aes(X1, X2, col = Class, label = prob)) +
+    geom_point(data = points.df, aes(X1, X2, col = Class)) +
+    geom_point(data = subset(points.df, points.df$prob != ""), aes(X1, X2),
+               shape = 1, col = "grey25")
+
+  # Touch up plot and return ---------------------------------------------------
+  p + coord_cartesian(xlim = x, ylim = y)
+
 }
 
 iplot_predict_mult <- function(plot.df, points.df, x, y, m, dec.bound) {
@@ -262,15 +294,15 @@ iplot_predict_mult <- function(plot.df, points.df, x, y, m, dec.bound) {
                alpha = alpha * plot.df[, class.ind[j]], fill = fill.col[j])
   }
 
-  # Add points and decision boundary, and touch up remaining plot --------------
+  # Add decision boundary ------------------------------------------------------
   if (isTRUE(dec.bound)) {
     decbound.df <- reshape2::melt(plot.df, id.vars = c("X1", "X2"))
     p <- p +
       geom_contour(data = decbound.df, aes(X1, X2, z = value, group = variable,
                                            col = variable,
                                            size = "Decision\nboundary"),
-                   # binwidth = 0.5 + 1e-12,
-                   bins = 2,
+                   binwidth = 0.5 + 1e-12,
+                   # bins = 2,
                    linetype = "dashed") +
       scale_size_manual(values = 0.8, name = NULL) +
       scale_color_discrete(name = " Class") +
@@ -303,7 +335,7 @@ iplot_predict_mult <- function(plot.df, points.df, x, y, m, dec.bound) {
 iplot_error <- function(x, niter.plot, plot.test = TRUE) {
   if (x$niter < 2) stop("Nothing to plot.")
   if (missing(niter.plot)) niter.plot <- seq_along(x$train.error)
-  else if (length(niter.plot) == 1) niter.plot <- c(1, niter.plot)
+  else if (length(niter.plot) == 1) niter.plot <- seq_len(niter.plot)
 
   # Prepare plotting data frame ------------------------------------------------
   plot.df <- data.frame(Iteration = niter.plot,
@@ -360,9 +392,9 @@ iplot_error <- function(x, niter.plot, plot.test = TRUE) {
     theme(legend.position = "none")
 }
 
-iplot_lb_and_error <- function(x, niter.plot, lab.pos) {
+iplot_lb_and_error <- function(x, niter.plot, lab.pos, plot.test = TRUE) {
   suppressMessages(
-    p2 <- iplot_error(x, niter.plot) +
+    p2 <- iplot_error(x, niter.plot, plot.test) +
       scale_x_continuous(
         breaks = scales::pretty_breaks(n = min(5, ifelse(x$niter == 2, 1, x$niter)))
       )
diff --git a/R/Utilities.R b/R/Utilities.R
index 5778b10..7c2b984 100644
--- a/R/Utilities.R
+++ b/R/Utilities.R
@@ -78,187 +78,38 @@ is.iprobitMod_mult <- function(x) inherits(x, "iprobitMod_mult")
 #' @export
 is.iprobitData <- function(x) inherits(x, "iprobitData")
 
-isNystrom <- function(x) {
-  if (inherits(x, "ipriorKernel_old")) {
-    if (!is.list(x$Nystrom)) res <- x$Nystrom
-    else res <- TRUE
-  } else {
-    if (!is.list(x$ipriorKernel$Nystrom)) res <- x$ipriorKernel$Nystrom
-    else res <- TRUE
-  }
-  res
-}
-
-#' Extract the variational lower bound
-#'
-#' @param object An object of class \code{ipriorProbit}.
-#' @param ... This is not used here.
-#'
-#' @return The variational lower bound.
-#' @export
-logLik.iprobitMod <- function(object, ...) {
-  lb <- object$lower.bound[!is.na(object$lower.bound)]
-  lb <- lb[length(lb)]
-  class(lb) <- "iprobitLowerBound"
-  lb
-}
-
-#' @export
-get_lbs <- function(x) x$lower.bound
-
-#' @export
-print.iprobitLowerBound <- function(x, ...) {
-  cat("Lower bound =", x)
-}
-
-get_theta <- function(object) object$theta
-
-get_w <- function(object) object$w
-
-get_y <- function(object) {
-  res <- factor(object$ipriorKernel$y)
-  levels(res) <- object$ipriorKernel$y.levels
-  res
-}
-
-
-#' @export
-get_one.lam <- function(object) {
-  object$ipriorKernel$model$one.lam
-}
-
-#' @export
-get_kernel <- function(object, collapse = TRUE) {
-  kernel.used <- object$ipriorKernel$model$kernel
-  Hurst.used <- get_Hurst(object)
-  for (i in seq_along(kernel.used)) {
-    if (kernel.used[i] == "FBM")
-      kernel.used[i] <- paste0(kernel.used[i], ",", Hurst.used[i])
-  }
-  kernel.used
-}
-
-#' get_lambda <- function(object) {
-#'   lambda <- object$lambda
-#'   if (is.iprobitMod_bin(object)) {
-#'     if (length(lambda) > 1)
-#'       names(lambda) <- paste0("lambda[", seq_along(lambda), "]")
-#'     else
-#'       names(lambda) <- "lambda"
-#'   } else if (is.iprobitMod_mult(object)) {
-#'     if (nrow(lambda) > 1)
-#'       rownames(lambda) <- paste0("lambda[", seq_along(lambda[, 1]), ",]")
-#'     else
-#'       rownames(lambda) <- "lambda"
-#'     colnames(lambda) <- paste0("Class = ", seq_along(object$y.levels))
-#'   } else {
-#'     stop("Input iprobitMod objects only.")
-#'   }
-#'
-#'   lambda
-#' }
-
-
-# get_psi <- function(object) {
-#   alpha <- object$alpha
-#   if (is.iprobitMod_bin(object)) {
-#     # if (length(alpha) > 1)
-#     #   names(alpha) <- paste0("alpha[", seq_along(alpha), "]")
-#     # else
-#     names(alpha) <- "alpha"
-#   } else if (is.iprobitMod_mult(object)) {
-#     alpha <- matrix(alpha, nrow = 1)
-#     rownames(alpha) <- "alpha"
-#     colnames(alpha) <- paste0("Class = ", seq_along(object$y.levels))
+# isNystrom <- function(x) {
+#   if (inherits(x, "ipriorKernel_old")) {
+#     if (!is.list(x$Nystrom)) res <- x$Nystrom
+#     else res <- TRUE
 #   } else {
-#     stop("Input iprobitMod objects only.")
+#     if (!is.list(x$ipriorKernel$Nystrom)) res <- x$ipriorKernel$Nystrom
+#     else res <- TRUE
 #   }
-#   alpha
+#   res
 # }
 
-get_alpha <- function(object) {
-  param.full <- param.summ_to_param.full(object$param.summ)
-  param.full[grep("Intercept", names(param.full))]
+is.common.intercept <- function(x) {
+  check_and_get_iprobitMod(x)
+  return(x$common$intercept)
 }
 
-get_lambda <- function(object) {
-  param.full <- param.summ_to_param.full(object$param.summ)
-  param.full[grep("lambda", names(param.full))]
+is.common.RKHS.param <- function(x) {
+  check_and_get_iprobitMod(x)
+  return(x$common$RKHS.param)
 }
 
-get_sd <- function(object) {
-  setNames(object$param.summ$S.D., rownames(object$param.summ))
-}
-
-get_sd_alpha <- function(object) {
-  res <- get_sd(object)
-  res[grep("Intercept", names(res))]
-}
-
-get_sd_lambda <- function(object) {
-  res <- get_sd(object)
-  res[grep("lambda", names(res))]
-}
-
-#' @export
-get_coef_se_mult <- function(object) {
-  theta <- coef(object)
-  m <- ncol(theta)
-  l <- nrow(theta) - 1
-
-  if (isTRUE(object$control$common.intercept)) {
-    alpha <- theta[1, 1]
-    names(alpha) <- "alpha"
-    alpha.se <- object$se.alpha
-  } else {
-    alpha <- theta[1, ]
-    names(alpha) <- paste0("alpha[", 1:m, "]")
-    alpha.se <- rep(object$se.alpha, m)
-  }
-
-  if (isTRUE(object$control$common.RKHS.scale)) {
-    lambda <- theta[-1, 1]
-    if (length(lambda) > 1)
-      names(lambda) <- paste0("lambda[", seq_along(lambda), ",]")
-    else
-      names(lambda) <- "lambda"
-    lambda.se <- object$se.lambda[, 1]
+check_and_get_iprobitMod <- function(object, assign.to.env = FALSE) {
+  # Helper function to check whether object is of ipriorMod class.
+  #
+  # Args: An ipriorMod or ipriorKernel object; logical assign.to.env.
+  #
+  # Returns: Nothing - just checks. Unless assign.to.env is TRUE.
+  if (is.iprobitMod(object)) {
+    if (isTRUE(assign.to.env)) list2env(object$ipriorKernel, parent.frame())
   } else {
-    lambda <- c(t(theta[-1, ]))
-    lambda.names <- NULL
-    for (k in 1:l) {
-      lambda.names <- c(lambda.names, paste0("lambda[", k, ",", 1:m, "]"))
-    }
-    names(lambda) <- lambda.names
-    lambda.se <- c(t(object$se.lambda))
+    stop("Input an iprobitMod object.", call. = FALSE)
   }
-
-  list(theta = c(alpha, lambda), se = c(alpha.se, lambda.se))
-}
-
-#' @export
-get_error_rate <- function(x) x$fitted.values$train.error
-
-#' @export
-get_error_rates <- function(x) {
-  res <- x$error
-  names(res) <- seq_along(res)
-  res
-}
-
-#' @export
-get_brier_score <- function(x) x$fitted.values$brier.score
-
-#' @export
-get_brier_scores <- function(x) {
-  res <- x$brier
-  names(res) <- seq_along(res)
-  res
-}
-
-get_m <- function(object) {
-  if (is.iprobitMod(object)) object <- object$ipriorKernel
-  length(object$y.levels)
 }
 
 all.same <- function(v) {
@@ -266,196 +117,6 @@ all.same <- function(v) {
   all(sapply(as.list(v[-1]), FUN = function(z) identical(z, v[1])))
 }
 
-# lambda expansion and Hlam calculation for binary models ----------------------
-
-expand_lambda <- function(x, intr, intr.3plus = NULL) {
-  # Helper function to expand lambda or lambda.sq (scale parameters) according
-  # to any interactions specification.
-  #
-  # Args: lambda
-  #
-  # Returns: Expanded lambda.
-  if (is.vector(x)) {
-    # Binary models ------------------------------------------------------------
-    return(iprior::.expand_Hl_and_lambda(x, x, intr, intr.3plus)$lambda)
-  } else {
-    # Multinomial models -------------------------------------------------------
-    res <- NULL
-    m <- ncol(x)
-    for (j in seq_len(m)) {
-      res[[j]] <- iprior::.expand_Hl_and_lambda(x[, j], x[, j], intr,
-                                                intr.3plus)$lambda
-    }
-    return(matrix(unlist(res), ncol = m))
-  }
-}
-
-get_Hlam <- function(object, theta, theta.is.lambda = FALSE) {
-  # Obtains the kernel matrix Hlam.
-  #
-  # Args:
-  #
-  # Returns: For binary models, this calculate Hlam. For multinomial models,
-  # this calculates Hlam for every class---so a list is returned.
-  if (is.iprobit_bin(object)) {
-    if (is.matrix(theta)) theta <- theta[, 1]
-    return(iprior::.get_Hlam(object = object, theta = theta,
-                             theta.is.lambda = theta.is.lambda))
-  } else {
-    res <- NULL
-    m <- get_m(object)
-    for (j in seq_len(m)) {
-      res[[j]] <- iprior::.get_Hlam(object = object, theta = as.numeric(theta[, j]),
-                                    theta.is.lambda = theta.is.lambda)
-    }
-    return(res)
-  }
-}
-
-get_Htildelam <- function(object, theta, xstar) {
-  if (is.iprobit_bin(object)) {
-    if (is.matrix(theta)) theta <- theta[, 1]
-    return(iprior::.get_Htildelam(object, theta, xstar))
-  } else {
-    res <- NULL
-    m <- get_m(object)
-    for (j in seq_len(m)) {
-      res[[j]] <- iprior::.get_Htildelam(object, as.numeric(theta[, j]), xstar)
-    }
-    return(res)
-  }
-}
-
-
-get_Hlamsq <- function() {
-  # Calculate Hlamsq for closed-form VB.
-  #
-  # Args: mod is an ipriorKernel object defined in parent environment, along
-  # with lambda, lambdasq, Psql, Hsql, H2l, ind1, ind2, p, and no.int.
-  #
-  # Returns: Hlamsq mat for binary models, and list of Hlamsq matrices for
-  # multinomial models.
-  environment(Hlam_two_way_index) <- environment()
-  q <- p + no.int
-  m <- get_m(mod)
-
-  if (is.iprobit_bin(mod)) { # BINARY MODEL
-    # Calculate square terms of Hlamsq -----------------------------------------
-    if (is.null(Hsql))
-      square.terms <- Reduce("+", mapply("*", Psql[1:q], lambdasq[1:q],
-                                         SIMPLIFY = FALSE))
-    else
-      square.terms <- Reduce("+", mapply("*", Hsql[1:q], lambdasq[1:q],
-                                         SIMPLIFY = FALSE))
-
-    # Calculate two-way terms of Hlamsq ----------------------------------------
-    if (is.null(ind1) && is.null(ind2)) {
-      two.way.terms <- 0
-    } else {
-      lambda.two.way <- Hlam_two_way_index(lambda, lambdasq)
-      two.way.terms <- Reduce("+", mapply("*", H2l, lambda.two.way,
-                                          SIMPLIFY = FALSE))
-    }
-
-    return(square.terms + two.way.terms)
-  } else {
-    res <- NULL
-    for (j in seq_len(m)) { # MULTINOMIAL MODEL
-      # Calculate square terms of Hlamsq -----------------------------------------
-      if (is.null(Hsql))
-        square.terms <- Reduce("+", mapply("*", Psql[1:q], lambdasq[1:q, j],
-                                           SIMPLIFY = FALSE))
-      else
-        square.terms <- Reduce("+", mapply("*", Hsql[1:q], lambdasq[1:q, j],
-                                           SIMPLIFY = FALSE))
-
-      # Calculate two-way terms of Hlamsq ----------------------------------------
-      if (is.null(ind1) && is.null(ind2))
-        two.way.terms <- 0
-      else {
-        lambda.two.way <- Hlam_two_way_index(lambda[1:q, j], lambdasq[1:q, j])
-        two.way.terms <-
-          Reduce("+", mapply("*", H2l, lambda.two.way, SIMPLIFY = FALSE))
-      }
-
-      res[[j]] <- square.terms + two.way.terms
-    }
-    return(res)
-  }
-}
-
-
-Hlam_two_way_index <- function(lam, lamsq) {
-  # mod <- iprior::.kernL(Species ~ . ^ 2, iris)
-  # iprobit.env <- environment()
-  # list2env(mod, iprobit.env)
-  # list2env(BlockBstuff, iprobit.env)
-  # list2env(model, iprobit.env)
-
-  comb.ind12 <- cbind(ind1, ind2)
-  comb.ind12 <- split(comb.ind12, row(comb.ind12))
-
-  replace_ind <- function(x) {
-    if (any(x > p)) {
-      here <- which(x > p)
-      what <- x[here] - p
-      res <- c(x[-here], intr[, what])
-      sort(res)
-    } else {
-      x
-    }
-  }
-
-  lam_ind <- function(x) {
-    here <- which(duplicated(x))
-    if (length(here > 0)) {
-      what <- x[here]
-      what.not <- x[x != what]
-      prod(lamsq[what]) * prod(lam[what.not])
-    } else {
-      prod(lam[x])
-    }
-  }
-
-  lapply(lapply(comb.ind12, replace_ind), lam_ind)
-}
-
-# myfun <- function() {
-#   # Function to test lambdaExpand_mult() and the Hlam*_mult functions
-#   set.seed(123)
-#   dat <- gen_mixture(n = 4, m = 3)
-#   mod <- iprior::.kernL(y ~ . ^ 2, dat)
-#   iprobit.env <- environment()
-#   list2env(mod, iprobit.env)
-#   list2env(BlockBstuff, iprobit.env)
-#   list2env(model, iprobit.env)
-#   environment(lambdaExpand_mult) <- iprobit.env
-#   environment(HlamFn_mult) <- iprobit.env
-#   environment(HlamsqFn_mult) <- iprobit.env
-#   environment(Hlam_two_way_index) <- iprobit.env
-#   m <- length(y.levels)
-#   lambda <- matrix(2:3, ncol = m, nrow = l)
-#   lambda.sq <- lambda ^ 2
-#   lambdaExpand_mult(x = lambda, y = lambda.sq, env = iprobit.env)
-#   HlamFn_mult(env = iprobit.env)
-#   HlamsqFn_mult(env = iprobit.env)
-#   list(H = Hl, Hsq = Hsql, lambda = lambda, lambda.sq = lambda.sq,
-#        Hlam.mat = Hlam.mat[[1]], Hlam.matsq = Hlam.matsq[[1]])
-#
-#   # CHECK
-#   # 2 3 6 (lambda)
-#   # 4 9 36 (lambda.sq)
-#   #
-#   # H1 <- H[[1]]
-#   # H2 <- H[[2]]
-#   # H12 <- H[[1]] * H[[2]]
-#   # H1sq <- H1 %*% H1
-#   # H2sq <- H2 %*% H2
-#   # H12sq <- H12 %*% H12
-#   #
-#   # 4 * H1sq + 9 * H2sq + 36 * H12sq + 2 * 3 * (H1 %*% H2 + H2 %*% H1) + 4 * 3 * (H1 %*% H12 + H12 %*% H1) + 9 * 2 * (H2 %*% H12 + H12 %*% H2)
-# }
-
 hyperparam_to_theta <- function(lambda, hurst, lengthscale, offset) {
   # In order to be compatible with the iprior helper functions, need theta (the
   # unconstrained versions of the hyperparameters, excluding alpha).
@@ -538,7 +199,6 @@ param.full_to_coef <- function(param.full, object) {
     iprior::.reduce_theta(x, est.list = object$estl)$theta.reduced
   })
   res <- rbind("Intercept" = param.full[1, ], res)
-  rownames(res) <- get_names(object, expand = FALSE)
   res
 }
 
diff --git a/R/fitted_and_predict.R b/R/fitted_and_predict.R
index 9ec9b6b..505747e 100644
--- a/R/fitted_and_predict.R
+++ b/R/fitted_and_predict.R
@@ -44,7 +44,10 @@
 fitted.iprobitMod <- function(object, quantiles = FALSE, n.samp = 100,
                               transform = identity, raw = FALSE, ...) {
   if (isTRUE(quantiles)) {
-    return(predict_quant(object, n.samp, transform, raw = raw, type = "train"))
+    y <- as.numeric(factor(object$ipriorKernel$y))
+
+    return(predict_quant(object, n.samp, transform, raw = raw, type = "train",
+                         y = y))
   } else {
     res <- object$fitted.values
     res$prob <- transform(res$prob)
@@ -101,17 +104,20 @@ predict_iprobit <- function(object, xstar, y.test, quantiles = FALSE,
   res
 }
 
-calc_ystar <- function(object, xstar, alpha, theta, w) {
+calc_ystar <- function(object, xstar, alpha, theta, w, theta.is.lambda = FALSE) {
   # Args: An ipriorKernel object.
   #
   # Returns: For binary models, a vector of ystar. For multinomial models, this
   # is a matrix with columns representing the number of classes.
-  if (is.null(xstar)) Hlam.new <- get_Hlam(object, theta)
-  else Hlam.new <- get_Htildelam(object, theta, xstar)
+  if (is.null(xstar)) Hlam.new <- get_Hlam(object, theta, theta.is.lambda)
+  else Hlam.new <- get_Htildelam(object, theta, xstar, theta.is.lambda)
   if (is.iprobit_bin(object)) {
+    if (length(alpha) > 1) {
+      stopifnot(all.same(as.numeric(alpha)))
+      alpha <- alpha[1]
+    }
     return(as.numeric(alpha + Hlam.new %*% w))
   } else {
-    res <- NULL
     m <- get_m(object)
     if (length(alpha) == 1) alpha <- rep(alpha, m)
     return(rep(alpha, each = nrow(Hlam.new[[1]])) +
@@ -120,6 +126,7 @@ calc_ystar <- function(object, xstar, alpha, theta, w) {
 }
 
 probs_yhat_error <- function(y, y.levels, ystar, type = c("train", "test")) {
+  # Args: y is for comparison (test or training responses).
   m <- length(y.levels)
 
   if (m == 2) {
@@ -150,7 +157,8 @@ probs_yhat_error <- function(y, y.levels, ystar, type = c("train", "test")) {
   levels(y.hat) <- y.levels
 
   error.rate <- mean(as.numeric(y.hat) != as.numeric(y)) * 100
-  brier.score <- brier_score(as.numeric(y), as.numeric(y.hat), as.data.frame(p.hat))
+  brier.score <- brier_score(as.numeric(y), as.numeric(y.hat),
+                             as.data.frame(p.hat))
 
   structure(list(y = y.hat, prob = as.data.frame(p.hat),
                  type = match.arg(type, c("train", "test")),
@@ -158,7 +166,6 @@ probs_yhat_error <- function(y, y.levels, ystar, type = c("train", "test")) {
             class = "iprobitPredict")
 }
 
-#' @export
 predict_quant <- function(object, n.samp, transform = identity, raw = FALSE,
                           xstar = NULL, y = NULL, type = c("train", "test")) {
   # Helper function to get the quantiles of the probabilities, error rate and
@@ -170,49 +177,73 @@ predict_quant <- function(object, n.samp, transform = identity, raw = FALSE,
     tmp <- sample_prob_bin(object, n.samp, xstar, y)
   }
   if (is.iprobitMod_mult(object)) {
-    # tmp <- sample_prob_mult(object, n.samp, Hl, y)
+    tmp <- sample_prob_mult(object, n.samp, xstar, y)
   }
 
   phats <- convert_prob(tmp$phat.samp, transform)
   names(phats) <- object$ipriorKernel$y.levels
+
+  if (all(is.nan(tmp$error.samp))) {
+    error.rate <- NaN
+  } else {
+    error.rate <- stats::quantile(tmp$error.samp,
+                                  probs = c(0.025, 0.25, 0.5, 0.75, 0.975))
+  }
+
+  if (all(is.nan(tmp$brier.samp))) {
+    brier.score <- NaN
+  } else {
+    brier.score <- stats::quantile(tmp$brier.samp,
+                                   probs = c(0.025, 0.25, 0.5, 0.75, 0.975))
+  }
+
   if (isTRUE(raw)) {
     return(list(
-      prob = phats, train.error = tmp$error.samp, brier.score = tmp$brier.samp
+      prob = phats, error.rate = tmp$error.samp, brier.score = tmp$brier.samp
     ))
   } else {
     return(structure(list(
       prob        = quantile_prob(phats),
-      error.rate  = stats::quantile(tmp$error.samp, probs = c(0.025, 0.25, 0.5, 0.75, 0.975)),
-      brier.score = stats::quantile(tmp$brier.samp, probs = c(0.025, 0.25, 0.5, 0.75, 0.975)),
+      error.rate  = error.rate,
+      brier.score = brier.score,
       type        = match.arg(type, c("train", "test"))
     ), class = "iprobitPredict_quant"))
   }
 }
 
-# calc_ystar <- function(object, xstar, alpha, theta, w)
-
 sample_prob_bin <- function(object, n.samp, xstar = NULL, y = NULL) {
   # Helper function to sample probabilities, error rates and brier scores for
-  # binary models.
+  # binary models. For the closed-form VB-EM, then the hyperparameters all have
+  # normal distributions. For the VB-EM with Metropolis sampler, then the
+  # samples for theta are Bootstrapped from object$theta.samp, a remnant of the
+  # fitting procedure. Note that results will not be so good if n.samp >>
+  # nrow(object$theta.samp).
   #
   # Args: An iprobitMod_x object.
+  #
+  # Returns: A list of 1) a list of class probability samples; 2) a vector of
+  # samples of  misclassification rates; and 3) a vector of sample Brier scores.
+  # This is needed in predict_quant().
   alpha.samp <- sample_alpha(n.samp, get_alpha(object), get_sd_alpha(object))
+  til <- FALSE
   if (grepl("Closed-form", object$est.method)) {
     theta.samp <- sample_lambda(n.samp, get_lambda(object),
                                 get_sd_lambda(object))
+    til <- TRUE
   } else {
-    stop("Not implemented yet.")
+    theta.samp <- object$theta.samp
+    theta.samp <- theta.samp[sample(seq_len(nrow(theta.samp)), size = n.samp), ]
   }
   w.samp <- sample_w(n.samp, object$w, object$Varw)  # n.samp x n matrix
+
   phat.samp <- list()
   error.samp <- brier.samp <- rep(NA, n.samp)
-  if (is.null(y)) y <- as.numeric(factor(object$ipriorKernel$y))
 
   for (i in seq_len(n.samp)) {
     alpha <- alpha.samp[i]
-    theta <- hyperparam_to_theta(theta.samp[i, ])
+    theta <- theta.samp[i, ]
     w <- w.samp[i, ]
-    ystar.samp <- calc_ystar(object$ipriorKernel, xstar, alpha, theta, w)
+    ystar.samp <- calc_ystar(object$ipriorKernel, xstar, alpha, theta, w, til)
     tmp <- probs_yhat_error(y, object$ipriorKernel$y.levels, ystar.samp)
     phat.samp[[i]] <- tmp$prob
     error.samp[i] <- tmp$error
@@ -222,144 +253,98 @@ sample_prob_bin <- function(object, n.samp, xstar = NULL, y = NULL) {
   list(phat.samp = phat.samp, error.samp = error.samp, brier.samp = brier.samp)
 }
 
-sample_alpha <- function(n.samp, mean, sd) {
-  rnorm(n.samp, mean, sd)
-}
-
-sample_lambda <- function(n.samp, mean, sd) {
-  if (length(mean) > 1) {
-    return(mvtnorm::rmvnorm(n.samp, mean, diag(sd ^ 2)))
-  } else {
-    return(matrix(log(rnorm(n.samp, mean, sd))))  # note the log(lambda)
-  }
-}
-
-sample_w <- function(n.samp, mean, variance) {
-  mvtnorm::rmvnorm(n.samp, mean, variance)
-}
-
-#' calc_ystar <- function(object, xstar) {
-#'   # Helper function to calculate ystar values from Hl, lambda, alpha and w. If
-#'   # values are not supplied in the argument, these are obtained from object.
-#'   #
-#'   # Args: iprobitMod object.
-#'   #
-#'   # Returns:
-#'   list2env(object, environment())
-#'   list2env(ipriorKernel, environment())
-#'   list2env(model, environment())
-#'   if (!is.null(Hl.new)) Hl <- Hl.new
-#'   if (!is.null(alpha.new)) alpha <- alpha.new
-#'   if (!is.null(lambda.new)) lambda <- lambda.new
-#'   if (!is.null(w.new)) w <- w.new
-#'
-#'   if (is.iprobitMod_bin(object)) {
-#'     environment(lambdaExpand_bin) <- environment()
-#'     environment(HlamFn) <- environment()
-#'     lambdaExpand_bin(env = environment(), y = NULL)
-#'     HlamFn(env = environment())
-#'     return(as.vector(alpha + (Hlam.mat %*% w)))
-#'   }
-#'   if (is.iprobitMod_mult(object)) {
-#'     # environment(lambdaExpand_mult) <- environment()
-#'     # environment(HlamFn_mult) <- environment()
-#'     # lambdaExpand_mult(env = environment(), y = NULL)
-#'     # HlamFn_mult(env = environment())
-#'     # return(rep(alpha, each = nrow(Hlam.mat[[1]])) +
-#'     #          mapply("%*%", Hlam.mat, split(w, col(w))))
-#'   }
-#' }
-
-
-
-
-
-
-
-
-
-
-# predict_iprobit_mult <- function(y, y.levels, ystar, type = c("train", "test")) {
-#   m <- length(y.levels)
-#
-#   p.hat <- ystar
-#   for (i in seq_len(nrow(ystar))) {
-#     for (j in seq_along(y.levels)) {
-#       p.hat[i, j] <- EprodPhiZ(ystar[i, j] - ystar[i, seq_along(y.levels)[-j]])
-#     }
-#     # p.hat[i, m] <- 1 - sum(p.hat[i, 1:(m - 1)])
-#   }
-#   p.hat <- p.hat / matrix(rep(apply(p.hat, 1, sum), m), ncol = m)  # normalise
-#   colnames(p.hat) <- y.levels
-#
-#   tmp <- p.hat
-#   y.hat <- factor(
-#     apply(tmp, 1, function(x) which(x == max(x))),
-#     levels = seq_len(m)
-#   )
-#   levels(y.hat) <- y.levels
-#
-#   error.rate <- mean(as.numeric(y.hat) != as.numeric(y)) * 100
-#   brier.score <- brier_score(as.numeric(y), as.numeric(y.hat), as.data.frame(p.hat))
-#
-#   type <- match.arg(type, c("train", "test"))
-#
-#   structure(list(y = y.hat, prob = as.data.frame(p.hat), type = type,
-#                  train.error = error.rate, brier.score = brier.score),
-#             class = "iprobit_predict")
-# }
-
+sample_prob_mult <- function(object, n.samp, xstar = NULL, y = NULL) {
+  # Helper function to sample probabilities, error rates and brier scores for
+  # multinomial models.
+  alpha.samp <- theta.samp <- w.samp <- phat.samp <- list()
+  m <- get_m(object)
+  til <- FALSE
+  vb.closed <- grepl("Closed-form", object$est.method)
 
+  alpha <- c(get_alpha(object, by.class = TRUE))
+  sd.alpha <- get_sd_alpha(object)
 
+  lambda <- get_lambda(object, by.class = TRUE)
+  sd.lambda <- matrix(get_sd_lambda(object), ncol = m, byrow = TRUE)
+  if (isTRUE(vb.closed)) til <- TRUE
 
-#' @export
-sample_prob_mult <- function(object, n.samp, Hl.new = NULL, y = NULL) {
-  # Helper function to sample probabilities, error rates and brier scores for
-  # multinomial models.
-  alpha.samp <- lambda.samp <- w.samp <- phat.samp <- list()
-  m <- object$ipriorKernel$m
+  # First, sample per class ----------------------------------------------------
   for (j in seq_len(m)) {
-    # First sample per class
-    alpha.samp[[j]] <- rnorm(n.samp, object$alpha[j], object$se.alpha)
-    if (nrow(object$lambda) > 1) {
-      lambda.samp[[j]] <- mvtnorm::rmvnorm(n.samp, object$lambda[, j],
-                                           diag(object$se.lambda[, j] ^ 2))
-    } else {
-      lambda.samp[[j]] <- matrix(rnorm(n.samp, object$lambda[, j],
-                                       object$se.lambda[, j]))
+    alpha.samp[[j]] <- sample_alpha(n.samp, alpha[j], sd.alpha[j])
+    if (isTRUE(vb.closed)) {
+      theta.samp[[j]] <- sample_lambda(n.samp, lambda[, j], sd.lambda[, j])
     }
-    w.samp[[j]] <- mvtnorm::rmvnorm(n.samp, object$w[, j], object$Varw[[j]])
+    w.samp[[j]] <- sample_w(n.samp, object$w[, j], object$Varw[[j]])
   }
-  # Then combine each sample into a list of length n.samp (makes it easier in
-  # the calculation of the probabilities etc.)
+
+  # Then combine each sample into a list of length n.samp (makes it easier in --
+  # the calculation of the probabilities etc.) ---------------------------------
   alpha.samp <- mapply(function(x) unlist(lapply(alpha.samp, `[`, x)),
                        seq_len(n.samp), SIMPLIFY = FALSE)
-  lambda.samp <- mapply(function(x) matrix(unlist(lapply(lambda.samp, `[`, x, )), ncol = m),
-                        seq_len(n.samp), SIMPLIFY = FALSE)
+  if (isTRUE(vb.closed)) {
+    theta.samp <- mapply(function(x) matrix(unlist(lapply(theta.samp, `[`, x, )),
+                                            ncol = m),
+                         seq_len(n.samp), SIMPLIFY = FALSE)
+  } else {
+    theta.samp <- object$theta.samp
+    theta.samp <- theta.samp[sample(seq_along(theta.samp), size = n.samp)]
+  }
   w.samp <- mapply(function(x) matrix(unlist(lapply(w.samp, `[`, x, )), ncol = m),
                    seq_len(n.samp), SIMPLIFY = FALSE)
-  phat.samp <- list()
+
   error.samp <- brier.samp <- rep(NA, n.samp)
-  if (is.null(y)) y <- object$ipriorKernel$Y
 
   for (i in seq_len(n.samp)) {
     alpha <- alpha.samp[[i]]
-    lambda <- lambda.samp[[i]]
+    theta <- theta.samp[[i]]
     w <- w.samp[[i]]
-    ystar.samp <- calc_ystar(object, Hl.new = Hl.new, alpha.new = alpha,
-                             lambda.new = lambda, w.new = w)
-    tmp <- predict_iprobit_mult(y, object$ipriorKernel$y.levels, ystar.samp)
+    ystar.samp <- calc_ystar(object$ipriorKernel, xstar, alpha, theta, w, til)
+    tmp <- probs_yhat_error(y, object$ipriorKernel$y.levels, ystar.samp)
     phat.samp[[i]] <- tmp$prob
-    error.samp[i] <- tmp$train.error
-    brier.samp[i] <- tmp$brier.score
+    error.samp[i] <- tmp$error
+    brier.samp[i] <- tmp$brier
   }
 
   list(phat.samp = phat.samp, error.samp = error.samp, brier.samp = brier.samp)
 }
 
+sample_alpha <- function(n.samp, mean, sd) {
+  # Helper function for sample_prob_bin(). Samples the intercept from the
+  # (approximated) posterior N(mean, sd^2).
+  #
+  # Args: The number of samples n.samp, the mean and sd of the normal
+  # distribution.
+  #
+  # Returns: A vector of length n.samp containing the intercept samples.
+  rnorm(n.samp, mean, sd)
+}
+
+sample_lambda <- function(n.samp, mean, sd) {
+  # Helper function for sample_prob_bin(). Samples the intercept from the
+  # (approximated) posterior N(mean, sd^2).
+  #
+  # Args: The number of samples n.samp, the mean and sd of the normal
+  # distribution.
+  #
+  # Returns: A matrix with n.samp rows containing the lambda samples.
+  if (length(mean) > 1) {
+    return(mvtnorm::rmvnorm(n.samp, mean, diag(sd ^ 2)))
+  } else {
+    return(matrix(rnorm(n.samp, mean, sd)))
+  }
+}
 
+sample_w <- function(n.samp, mean, variance) {
+  # Helper function for sample_prob_bin(). Samples the intercept from the
+  # (approximated) posterior N(mean, variance).
+  #
+  # Args: The number of samples n.samp, the mean and variance of the
+  # multivariate normal distribution.
+  #
+  # Returns: A matrix with n.samp rows containing the w samples.
+  mvtnorm::rmvnorm(n.samp, mean, variance)
+}
 
-#' @export
 convert_prob <- function(phat, transform = identity) {
   # Converts list of length n.samp containing random samples of the
   # probabilities. This function converts it into a list of length no.classes so
@@ -373,7 +358,6 @@ convert_prob <- function(phat, transform = identity) {
   res
 }
 
-#' @export
 quantile_prob <- function(phats) {
   # Get quantiles for lists returned from convert_prob()
   for (i in seq_along(phats)) {
@@ -383,7 +367,6 @@ quantile_prob <- function(phats) {
   lapply(phats, as.data.frame)
 }
 
-#' @export
 print.iprobitPredict <- function(x, rows = 10, digits = 3, ...) {
   if (x$type == "train") {
     cat("Training error:", paste0(decimal_place(x$error.rate, digits), "%\n"))
@@ -413,7 +396,6 @@ print.iprobitPredict <- function(x, rows = 10, digits = 3, ...) {
   if (nrow(x$p) > rows) cat("# ... with", nrow(x$p) - rows, "more rows")
 }
 
-#' @export
 print.iprobitPredict_quant <- function(x, rows = 5, digits = 3, ...) {
   rows.act <- nrow(x$prob[[1]])
   rows <- min(rows.act, rows)
diff --git a/R/iprobit.R b/R/iprobit.R
index c1c11c2..9dec774 100644
--- a/R/iprobit.R
+++ b/R/iprobit.R
@@ -30,7 +30,7 @@ iprobit.default <- function(y, ..., kernel = "linear", interactions = NULL,
                             train.samp, control = list()) {
   # Load the I-prior model -----------------------------------------------------
   if (iprior::is.ipriorKernel(y)) {
-    mod <- y
+    mod <- remove_psi(y)  # remove psi from estimation procedure if not already
   } else {
     mod <- iprior::kernL(y, ..., kernel = kernel, interactions = interactions,
                          est.lambda = TRUE, est.psi = FALSE, psi = 1,
@@ -42,16 +42,17 @@ iprobit.default <- function(y, ..., kernel = "linear", interactions = NULL,
 
   # Set up controls ------------------------------------------------------------
   control_ <- list(
-    maxit          = 10,
+    maxit          = 100,
     stop.crit      = 1e-5,
     silent         = FALSE,
     alpha0         = NULL,  # if NULL, parameters are
-    lambda0        = NULL,  # initialised in VB
+    # lambda0        = NULL,  # initialised in VB
     w0             = NULL,  # routine
     theta0         = NULL,
     n.samp         = 100,  # settings for
     sd.samp        = 0.15,   # the metropolis
     thin.samp      = 2,    # sampler
+    seed           = NULL,
     restarts       = 0,
     restart.method = c("lb", "error", "brier")
   )
@@ -74,11 +75,11 @@ iprobit.default <- function(y, ..., kernel = "linear", interactions = NULL,
     if (m == 2) {
       # Binary models ----------------------------------------------------------
       if (est.method["em.closed"]) {  # VB CLOSED-FORM
-        res <- iprobit_bin(mod, maxit, stop.crit, silent, alpha0, lambda0, w0)
+        res <- iprobit_bin(mod, maxit, stop.crit, silent, alpha0, theta0, w0)
         res$est.method <- "Closed-form VB-EM algorithm."
       } else {
         res <- iprobit_bin_metr(mod, maxit, stop.crit, silent, alpha0, theta0,
-                                w0, n.samp, sd.samp, thin.samp)
+                                w0, n.samp, sd.samp, thin.samp, seed)
         res$est.method <- paste0("VB-EM with Metropolis sampler (",
                                  iprior::dec_plac(mean(res$acc.rate) * 100, 1),
                                  "% acc.).")
@@ -87,12 +88,16 @@ iprobit.default <- function(y, ..., kernel = "linear", interactions = NULL,
     } else {
       # Multinomial models -----------------------------------------------------
       if (est.method["em.closed"]) {  # VB CLOSED-FORM
-        res <- iprobit_mult(mod, maxit, stop.crit, silent, alpha0, lambda0, w0,
+        res <- iprobit_mult(mod, maxit, stop.crit, silent, alpha0, theta0, w0,
                             common.intercept, common.RKHS.scale)
         res$est.method <- "Closed-form VB-EM algorithm."
-        # res$coefficients <- rbind(get_alpha(res), get_lambda(res))
       } else {
-        stop("Not implemented yet.")
+        res <- iprobit_mult_metr(mod, maxit, stop.crit, silent, alpha0, theta0,
+                                 w0, n.samp, sd.samp, thin.samp, seed,
+                                 common.intercept, common.RKHS.scale)
+        res$est.method <- paste0("VB-EM with Metropolis sampler (",
+                                 iprior::dec_plac(mean(res$acc.rate) * 100, 1),
+                                 "% acc.).")
       }
       class(res) <- c("iprobitMod", "iprobitMod_mult")
     }
@@ -106,6 +111,7 @@ iprobit.default <- function(y, ..., kernel = "linear", interactions = NULL,
     res$ipriorKernel <- mod
   }
   res$coefficients <- param.full_to_coef(res$param.full, mod)
+  # rownames(res$coefficients) <- get_names(mod, expand = FALSE)
 
   # Change the call to "iprobit" -----------------------------------------------
   res$call <- iprior::.fix_call_default(match.call(), "iprobit")
@@ -146,7 +152,6 @@ iprobit.formula <- function(formula, data, kernel = "linear", one.lam = FALSE,
 iprobit.ipriorKernel <- function(object, control = list(),
                                  common.intercept = FALSE,
                                  common.RKHS.scale = FALSE, ...) {
-  object <- remove_psi(object)
   res <- iprobit.default(y = object, method = method, control = control,
                          common.intercept = common.intercept,
                          common.RKHS.scale = common.RKHS.scale)
diff --git a/R/iprobit_bin.R b/R/iprobit_bin.R
index 0e6e833..bbb25cb 100644
--- a/R/iprobit_bin.R
+++ b/R/iprobit_bin.R
@@ -18,9 +18,8 @@
 #
 ################################################################################
 
-#' @export
 iprobit_bin <- function(mod, maxit = 100, stop.crit = 1e-5, silent = FALSE,
-                        alpha0 = NULL, lambda0 = NULL, w0 = NULL) {
+                        alpha0 = NULL, theta0 = NULL, w0 = NULL) {
   # Declare all variables and functions to be used into environment ------------
   iprobit.env <- environment()
   list2env(mod, iprobit.env)
@@ -33,8 +32,12 @@ iprobit_bin <- function(mod, maxit = 100, stop.crit = 1e-5, silent = FALSE,
 
   # Initialise -----------------------------------------------------------------
   if (is.null(alpha0)) alpha0 <- rnorm(1)
-  if (is.null(lambda0)) lambda0 <- abs(rnorm(p))
+  if (is.null(theta0)) theta0 <- rnorm(p)
   if (is.null(w0)) w0 <- rep(0, n)
+  if (p == 1) lambda0 <- exp(theta0)  # ensures positive values
+  else lambda0 <- theta0
+
+  alpha <- alpha0
   lambda <- ct <- lambda0
   lambdasq <- lambda0 ^ 2
   Hl <- iprior::.expand_Hl_and_lambda(Hl, rep(1, p), intr, intr.3plus)$Hl  # expand Hl
@@ -42,7 +45,6 @@ iprobit_bin <- function(mod, maxit = 100, stop.crit = 1e-5, silent = FALSE,
   lambdasq <- expand_lambda(lambdasq[1:p], intr, intr.3plus)
   Hlam   <- get_Hlam(mod, lambda, theta.is.lambda = TRUE)
   Hlamsq <- get_Hlamsq()
-  alpha <- alpha0
   w <- w0
   Varw <- NA
   f.tmp <- as.numeric(alpha + Hlam %*% w)
@@ -105,7 +107,7 @@ iprobit_bin <- function(mod, maxit = 100, stop.crit = 1e-5, silent = FALSE,
       (sum(diag(W)) + determinant(A)$modulus + sum(log(ct))) / 2
 
     # theta --------------------------------------------------------------------
-    if (length(lambda) == 1) theta <- log(lambda[1:p])
+    if (p == 1) theta <- log(lambda[1])
     else theta <- lambda[1:p]
     theta <- matrix(theta, ncol = 2, nrow = length(theta))
 
@@ -138,7 +140,7 @@ iprobit_bin <- function(mod, maxit = 100, stop.crit = 1e-5, silent = FALSE,
     `97.5%` = c(alpha, lambda[1:p]) + qnorm(0.975) * sqrt(c(1 / n, 1 / ct[1:p]))
   )
   colnames(param.summ) <- c("Mean", "S.D.", "2.5%", "97.5%")
-  rownames(param.summ) <- get_names(mod, expand = FALSE)
+  rownames(param.summ) <- get_names(mod, c("intercept", "lambda"), FALSE)
 
   # Clean up and close ---------------------------------------------------------
   convergence <- niter == maxit
diff --git a/R/iprobit_bin_metr.R b/R/iprobit_bin_metr.R
index c25ea6a..17763a3 100644
--- a/R/iprobit_bin_metr.R
+++ b/R/iprobit_bin_metr.R
@@ -1,6 +1,40 @@
+################################################################################
+#
+#   iprobit: Binary and Multinomial Probit Regression with I-priors
+#   Copyright (C) 2017  Haziq Jamil
+#
+#   This program is free software: you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation, either version 3 of the License, or
+#   (at your option) any later version.
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License
+#   along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+# The updating equation for theta depends on a density that cannot be
+# recognised. Therefore, in order to obtain posterior samples from this density,
+# we must emebed a Metropolis step within the VB-EM. After getting samples
+# theta.samp, we can then calculate any statistic of interest, including
+# E[H_theta], E[H_theta^2], the unconstrained version of the parameters
+# g^{-1})(theta), and so on. Uncertainty of these estimates can be summarised
+# through the standard deviation of the samples, and we also calculate the
+# quantiles c(0.025, 0.975) for the relevant samples.
+#
+# In order to obtain E[H_theta] and E[H_theta^2], we create a list Hlaml (c.f.
+# Hlamsql) of length n.samp. For each value of theta.samp, calculate H_theta and
+# store it in the list. The estimate would then be the sample mean.
+
 iprobit_bin_metr <- function(mod, maxit = 20, stop.crit = 1e-5, silent = FALSE,
                              alpha0 = NULL, theta0 = NULL, w0 = NULL,
-                             n.samp = 1000, sd.samp = 0.1, thin.samp = 10) {
+                             n.samp = 1000, sd.samp = 0.15, thin.samp = 10,
+                             seed = NULL) {
   # Declare all variables and functions to be used into environment ------------
   iprobit.env <- environment()
   list2env(mod, iprobit.env)
@@ -10,6 +44,7 @@ iprobit_bin_metr <- function(mod, maxit = 20, stop.crit = 1e-5, silent = FALSE,
   thin.seq <- seq_len(n.samp)[seq_len(n.samp) %% thin.samp == 0]
 
   # Initialise -----------------------------------------------------------------
+  if (!is.null(seed)) set.seed(seed)
   if (is.null(alpha0)) alpha0 <- rnorm(1)
   if (is.null(theta0)) theta0 <- rnorm(thetal$n.theta)
   if (is.null(w0)) w0 <- rep(0, n)
@@ -83,6 +118,7 @@ iprobit_bin_metr <- function(mod, maxit = 20, stop.crit = 1e-5, silent = FALSE,
     }
     acc.rate[niter] <- count / n.samp
     theta <- apply(theta.samp[thin.seq, ], 2, mean)
+    theta <- matrix(theta, ncol = 2, nrow = length(theta))
     Hlam <- Reduce("+", Hlaml[thin.seq]) / length(Hlaml[thin.seq])
     Hlamsq <- Reduce("+", Hlamsql[thin.seq]) / length(Hlamsql[thin.seq])
 
@@ -116,10 +152,12 @@ iprobit_bin_metr <- function(mod, maxit = 20, stop.crit = 1e-5, silent = FALSE,
   time.taken <- iprior::as.time(end.time - start.time)
 
   # Calculate posterior s.d. and prepare param table ---------------------------
+  param.full <- theta_to_param.full(theta, alpha, mod)
   param.summ <- rbind(
-    alpha = c(alpha, 1 / n, alpha - qnorm(0.975) / n, alpha + qnorm(0.975) / n),
-    theta.samp_to_param.summ(theta.samp[thin.seq, ], mod)
+    c(alpha, 1 / n, alpha - qnorm(0.975) / n, alpha + qnorm(0.975) / n),
+    theta.samp_to_summ(theta.samp[thin.seq, ], mod)
   )
+  rownames(param.summ) <- get_names(mod, expand = FALSE)
 
   # Clean up and close ---------------------------------------------------------
   convergence <- niter == maxit
@@ -129,8 +167,8 @@ iprobit_bin_metr <- function(mod, maxit = 20, stop.crit = 1e-5, silent = FALSE,
     else cat("Converged after", niter, "iterations.\n")
   }
 
-  list(theta = theta, param.summ = param.summ, w = w, Varw = Varw,
-       lower.bound = as.numeric(na.omit(lb)), niter = niter,
+  list(theta = theta, param.full = param.full, param.summ = param.summ, w = w,
+       Varw = Varw, lower.bound = as.numeric(na.omit(lb)), niter = niter,
        start.time = start.time, end.time = end.time, time = time.taken,
        fitted.values = fitted.values, test = fitted.test,
        train.error = as.numeric(na.omit(train.error)),
@@ -141,7 +179,7 @@ iprobit_bin_metr <- function(mod, maxit = 20, stop.crit = 1e-5, silent = FALSE,
        theta.samp = theta.samp[thin.seq, ])
 }
 
-theta.samp_to_param.summ <- function(theta.samp, object) {
+theta.samp_to_summ <- function(theta.samp, object) {
   # Args: theta.samp is the sample from the (approximate) posterior q(theta),
   # and object is an ipriorKernel object.
   theta <- apply(theta.samp, 2, mean)
@@ -151,14 +189,14 @@ theta.samp_to_param.summ <- function(theta.samp, object) {
   ginv.theta.sd <- apply(ginv.theta.samp, 2, sd)
   ginv.theta.0025 <- apply(ginv.theta.samp, 2, stats::quantile, probs = 0.025)
   ginv.theta.0975 <- apply(ginv.theta.samp, 2, stats::quantile, probs = 0.975)
-  param.summ <- data.frame(
+  res <- data.frame(
     Mean    = ginv.theta,
     S.D.    = ginv.theta.sd,
     `2.5%`  = ginv.theta.0025,
     `97.5%` = ginv.theta.0975
   )
-  colnames(param.summ) <- c("Mean", "S.D.", "2.5%", "97.5%")
-  param.summ
+  colnames(res) <- c("Mean", "S.D.", "2.5%", "97.5%")
+  res
 }
 
 
diff --git a/R/iprobit_bin_metr (backup).R b/R/iprobit_bin_metr_backup.R
similarity index 100%
rename from R/iprobit_bin_metr (backup).R
rename to R/iprobit_bin_metr_backup.R
diff --git a/R/iprobit_helper.R b/R/iprobit_helper.R
index 116b415..8f8a619 100644
--- a/R/iprobit_helper.R
+++ b/R/iprobit_helper.R
@@ -18,6 +18,192 @@
 #
 ################################################################################
 
+expand_lambda <- function(x, intr, intr.3plus = NULL) {
+  # Helper function to expand lambda or lambda.sq (scale parameters) according
+  # to any interactions specification.
+  #
+  # Args: lambda
+  #
+  # Returns: Expanded lambda.
+  if (is.vector(x)) {
+    # Binary models ------------------------------------------------------------
+    return(iprior::.expand_Hl_and_lambda(x, x, intr, intr.3plus)$lambda)
+  } else {
+    # Multinomial models -------------------------------------------------------
+    res <- NULL
+    m <- ncol(x)
+    for (j in seq_len(m)) {
+      res[[j]] <- iprior::.expand_Hl_and_lambda(x[, j], x[, j], intr,
+                                                intr.3plus)$lambda
+    }
+    return(matrix(unlist(res), ncol = m))
+  }
+}
+
+get_Hlam <- function(object, theta, theta.is.lambda = FALSE) {
+  # Obtains the kernel matrix Hlam.
+  #
+  # Args:
+  #
+  # Returns: For binary models, this calculate Hlam. For multinomial models,
+  # this calculates Hlam for every class---so a list is returned.
+  if (is.iprobit_bin(object)) {
+    if (is.matrix(theta)) theta <- theta[, 1]
+    return(iprior::.get_Hlam(object, theta, theta.is.lambda))
+  } else {
+    res <- NULL
+    m <- get_m(object)
+    for (j in seq_len(m)) {
+      res[[j]] <- iprior::.get_Hlam(object, as.numeric(theta[, j]),
+                                    theta.is.lambda)
+    }
+    return(res)
+  }
+}
+
+get_Htildelam <- function(object, theta, xstar, theta.is.lambda = FALSE) {
+  if (is.iprobit_bin(object)) {
+    if (is.matrix(theta)) theta <- theta[, 1]
+    return(iprior::.get_Htildelam(object, theta, xstar, theta.is.lambda))
+  } else {
+    res <- NULL
+    m <- get_m(object)
+    for (j in seq_len(m)) {
+      res[[j]] <- iprior::.get_Htildelam(object, as.numeric(theta[, j]), xstar,
+                                         theta.is.lambda)
+    }
+    return(res)
+  }
+}
+
+get_Hlamsq <- function() {
+  # Calculate Hlamsq for closed-form VB.
+  #
+  # Args: mod is an ipriorKernel object defined in parent environment, along
+  # with lambda, lambdasq, Psql, Hsql, H2l, ind1, ind2, p, and no.int.
+  #
+  # Returns: Hlamsq mat for binary models, and list of Hlamsq matrices for
+  # multinomial models.
+  environment(Hlam_two_way_index) <- environment()
+  q <- p + no.int
+  m <- get_m(mod)
+
+  if (is.iprobit_bin(mod)) { # BINARY MODEL
+    # Calculate square terms of Hlamsq -----------------------------------------
+    if (is.null(Hsql))
+      square.terms <- Reduce("+", mapply("*", Psql[1:q], lambdasq[1:q],
+                                         SIMPLIFY = FALSE))
+    else
+      square.terms <- Reduce("+", mapply("*", Hsql[1:q], lambdasq[1:q],
+                                         SIMPLIFY = FALSE))
+
+    # Calculate two-way terms of Hlamsq ----------------------------------------
+    if (is.null(ind1) && is.null(ind2)) {
+      two.way.terms <- 0
+    } else {
+      lambda.two.way <- Hlam_two_way_index(lambda, lambdasq)
+      two.way.terms <- Reduce("+", mapply("*", H2l, lambda.two.way,
+                                          SIMPLIFY = FALSE))
+    }
+
+    return(square.terms + two.way.terms)
+  } else {
+    res <- NULL
+    for (j in seq_len(m)) { # MULTINOMIAL MODEL
+      # Calculate square terms of Hlamsq -----------------------------------------
+      if (is.null(Hsql))
+        square.terms <- Reduce("+", mapply("*", Psql[1:q], lambdasq[1:q, j],
+                                           SIMPLIFY = FALSE))
+      else
+        square.terms <- Reduce("+", mapply("*", Hsql[1:q], lambdasq[1:q, j],
+                                           SIMPLIFY = FALSE))
+
+      # Calculate two-way terms of Hlamsq ----------------------------------------
+      if (is.null(ind1) && is.null(ind2))
+        two.way.terms <- 0
+      else {
+        lambda.two.way <- Hlam_two_way_index(lambda[1:q, j], lambdasq[1:q, j])
+        two.way.terms <-
+          Reduce("+", mapply("*", H2l, lambda.two.way, SIMPLIFY = FALSE))
+      }
+
+      res[[j]] <- square.terms + two.way.terms
+    }
+    return(res)
+  }
+}
+
+Hlam_two_way_index <- function(lam, lamsq) {
+  # mod <- iprior::.kernL(Species ~ . ^ 2, iris)
+  # iprobit.env <- environment()
+  # list2env(mod, iprobit.env)
+  # list2env(BlockBstuff, iprobit.env)
+  # list2env(model, iprobit.env)
+
+  comb.ind12 <- cbind(ind1, ind2)
+  comb.ind12 <- split(comb.ind12, row(comb.ind12))
+
+  replace_ind <- function(x) {
+    if (any(x > p)) {
+      here <- which(x > p)
+      what <- x[here] - p
+      res <- c(x[-here], intr[, what])
+      sort(res)
+    } else {
+      x
+    }
+  }
+
+  lam_ind <- function(x) {
+    here <- which(duplicated(x))
+    if (length(here > 0)) {
+      what <- x[here]
+      what.not <- x[x != what]
+      prod(lamsq[what]) * prod(lam[what.not])
+    } else {
+      prod(lam[x])
+    }
+  }
+
+  lapply(lapply(comb.ind12, replace_ind), lam_ind)
+}
+
+# myfun <- function() {
+#   # Function to test lambdaExpand_mult() and the Hlam*_mult functions
+#   set.seed(123)
+#   dat <- gen_mixture(n = 4, m = 3)
+#   mod <- iprior::.kernL(y ~ . ^ 2, dat)
+#   iprobit.env <- environment()
+#   list2env(mod, iprobit.env)
+#   list2env(BlockBstuff, iprobit.env)
+#   list2env(model, iprobit.env)
+#   environment(lambdaExpand_mult) <- iprobit.env
+#   environment(HlamFn_mult) <- iprobit.env
+#   environment(HlamsqFn_mult) <- iprobit.env
+#   environment(Hlam_two_way_index) <- iprobit.env
+#   m <- length(y.levels)
+#   lambda <- matrix(2:3, ncol = m, nrow = l)
+#   lambda.sq <- lambda ^ 2
+#   lambdaExpand_mult(x = lambda, y = lambda.sq, env = iprobit.env)
+#   HlamFn_mult(env = iprobit.env)
+#   HlamsqFn_mult(env = iprobit.env)
+#   list(H = Hl, Hsq = Hsql, lambda = lambda, lambda.sq = lambda.sq,
+#        Hlam.mat = Hlam.mat[[1]], Hlam.matsq = Hlam.matsq[[1]])
+#
+#   # CHECK
+#   # 2 3 6 (lambda)
+#   # 4 9 36 (lambda.sq)
+#   #
+#   # H1 <- H[[1]]
+#   # H2 <- H[[2]]
+#   # H12 <- H[[1]] * H[[2]]
+#   # H1sq <- H1 %*% H1
+#   # H2sq <- H2 %*% H2
+#   # H12sq <- H12 %*% H12
+#   #
+#   # 4 * H1sq + 9 * H2sq + 36 * H12sq + 2 * 3 * (H1 %*% H2 + H2 %*% H1) + 4 * 3 * (H1 %*% H12 + H12 %*% H1) + 9 * 2 * (H2 %*% H12 + H12 %*% H2)
+# }
+
 # Helper function to determine when to stop the while loop for the VB-EM
 # algorithm.
 # If niter == 0 return TRUE because must complete 1 iteration.
diff --git a/R/iprobit_mult.R b/R/iprobit_mult.R
index 74d53dc..1494f47 100644
--- a/R/iprobit_mult.R
+++ b/R/iprobit_mult.R
@@ -18,9 +18,8 @@
 #
 ################################################################################
 
-#' @export
 iprobit_mult <- function(mod, maxit = 10, stop.crit = 1e-5, silent = FALSE,
-                         alpha0 = NULL, lambda0 = NULL, w0 = NULL,
+                         alpha0 = NULL, theta0 = NULL, w0 = NULL,
                          common.intercept = FALSE, common.RKHS.scale = FALSE) {
   # Declare all variables and functions to be used into environment ------------
   iprobit.env <- environment()
@@ -39,12 +38,15 @@ iprobit_mult <- function(mod, maxit = 10, stop.crit = 1e-5, silent = FALSE,
     if (isTRUE(common.intercept)) alpha0 <- rep(rnorm(1), m)
     else alpha0 <- rnorm(m)
   }
-  if (is.null(lambda0)) {
-    if (isTRUE(common.RKHS.scale)) lambda0 <- rep(abs(rnorm(p)), m)
-    else lambda0 <- abs(rnorm(p * m))
+  if (is.null(theta0)) {
+    if (isTRUE(common.RKHS.scale)) theta0 <- rep(rnorm(p), m)
+    else theta0 <- rnorm(p * m)
   }
+  if (p == 1) lambda0 <- exp(theta0)
+  else lambda0 <- theta0
   if (is.null(w0)) w0 <- matrix(0, ncol = m, nrow = n)
-  alpha <- rep(1, m); alpha[1:m] <- alpha0
+
+  alpha <- alpha0
   theta <- lambda <- ct <- dt <- matrix(lambda0, ncol = m, nrow = p)
   lambdasq <- lambda ^ 2
   Hl <- iprior::.expand_Hl_and_lambda(Hl, rep(1, p), intr, intr.3plus)$Hl  # expand Hl
@@ -141,8 +143,8 @@ iprobit_mult <- function(mod, maxit = 10, stop.crit = 1e-5, silent = FALSE,
     if (isTRUE(common.intercept)) alpha <- rep(mean(alpha), m)
 
     # theta --------------------------------------------------------------------
-    if (nrow(lambda) == 1) theta <- log(lambda)
-    else theta <- lambda
+    if (p == 1) theta <- log(lambda[1, , drop = FALSE])
+    else theta <- lambda[1:p, , drop = FALSE]
 
     # Calculate lower bound ----------------------------------------------------
     lb.ystar <- sum(logClb)
@@ -221,11 +223,6 @@ iprobit_mult <- function(mod, maxit = 10, stop.crit = 1e-5, silent = FALSE,
        test.brier = as.numeric(na.omit(test.brier)), convergence = convergence)
 }
 
-# set.seed(123)
-# dat <- gen_mixture(n = 400, m = 4)
-# mod <- kernL(y ~ ., dat, est.psi = FALSE)
-# set.seed(123); tmp <- iprobit_mult(mod)
-
 # else {
 #   # Nystrom approximation
 #   K.mm <- Hlamsq[[j]][1:Nystrom$m, 1:Nystrom$m]
diff --git a/R/iprobit_mult_metr.R b/R/iprobit_mult_metr.R
new file mode 100644
index 0000000..6587966
--- /dev/null
+++ b/R/iprobit_mult_metr.R
@@ -0,0 +1,306 @@
+################################################################################
+#
+#   iprobit: Binary and Multinomial Probit Regression with I-priors
+#   Copyright (C) 2017  Haziq Jamil
+#
+#   This program is free software: you can redistribute it and/or modify
+#   it under the terms of the GNU General Public License as published by
+#   the Free Software Foundation, either version 3 of the License, or
+#   (at your option) any later version.
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#   GNU General Public License for more details.
+#
+#   You should have received a copy of the GNU General Public License
+#   along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+# theta.samp would now be a list of length n.samp and each component of this
+# list is the matrix theta.samp.
+#
+# Similarly Hlaml and Hlamsql are lists of length n.samp, and each component of
+# the list is another list of length  m (no. of classes) where the values
+# Hlam_theta for each theta.samp are stored.
+
+iprobit_mult_metr <- function(mod, maxit = 5, stop.crit = 1e-5, silent = FALSE,
+                              alpha0 = NULL, theta0 = NULL, w0 = NULL,
+                              n.samp = 100, sd.samp = 0.15, thin.samp = 1,
+                              seed = NULL, common.intercept = FALSE,
+                              common.RKHS.scale = FALSE) {
+  # Declare all variables and functions to be used into environment ------------
+  iprobit.env <- environment()
+  list2env(mod, iprobit.env)
+  environment(loop_logical) <- iprobit.env
+  y <- as.numeric(factor(y))  # as.factor then as.numeric to get y = 1, 2, ...
+  m <- length(y.levels)
+  nm <- n * m
+  maxit <- max(1, maxit)
+  thin.seq <- seq_len(n.samp)[seq_len(n.samp) %% thin.samp == 0]
+
+  # Initialise -----------------------------------------------------------------
+  if (!is.null(seed)) set.seed(seed)
+  if (is.null(alpha0)) {
+    if (isTRUE(common.intercept)) alpha0 <- rep(rnorm(1), m)
+    else alpha0 <- rnorm(m)
+  }
+  if (is.null(theta0)) {
+    if (isTRUE(common.RKHS.scale)) theta0 <- rep(abs(rnorm(p)), m)
+    else theta0 <- abs(rnorm(p * m))
+  }
+  if (is.null(w0)) w0 <- matrix(0, ncol = m, nrow = n)
+
+  alpha <- alpha0
+  theta.samp <- Hlaml <- Hlamsql <- list(NULL)
+  theta <- matrix(theta0, ncol = m, nrow = thetal$n.theta)
+  Hlaml[[1]] <- Hlam <- get_Hlam(mod, theta)
+  Hlamsql[[1]] <- Hlamsq <- lapply(Hlaml[[1]], iprior::fastSquare)
+  theta.samp[[1]] <- theta
+
+  w <- f.tmp <- ystar <- w0
+  Varw <- W <- list(NULL)
+  f.tmp <- rep(alpha, each = n) + mapply("%*%", Hlam, split(w, col(w)))
+  logdetA <- rep(NA, m)
+  w.ind <- seq_len(
+    ifelse(isTRUE(common.intercept) && isTRUE(common.RKHS.scale), 1, m)
+  )
+  niter <- 0
+  lb <- train.error <- train.brier <- test.error <- test.brier <- rep(NA, maxit)
+  logClb <- rep(NA, n)
+  log.q.star <- rep(NA, m)
+  log.q <- matrix(NA, ncol = m, nrow = n.samp)
+  log.q[1, ] <- sapply(Hlamsq, function(x) -0.5 * sum(diag(x)))
+  acc.rate <- matrix(NA, nrow = maxit, ncol = m)
+
+  # The variational EM loop ----------------------------------------------------
+  if (!silent) pb <- txtProgressBar(min = 0, max = maxit, style = 1)
+  start.time <- Sys.time()
+
+  while (loop_logical()) {  # see loop_logical() function in iprobit_helper.R
+    # Update ystar -------------------------------------------------------------
+    for (i in 1:n) {
+      j <- as.numeric(y[i])
+      fi <- f.tmp[i, ]
+      fik <- fi[-j]; fij <- fi[j]
+      logClb[i] <- logC <- EprodPhiZ(fij - fik, log = TRUE)
+      for (k in seq_len(m)[-j]) {
+        logD <- log(integrate(
+          function(z) {
+            fij.minus.fil <- fij - fi[-c(k, j)]
+            logPhi.l <- 0
+            for (kk in seq_len(length(fij.minus.fil)))
+              logPhi.l <- logPhi.l + pnorm(z + fij.minus.fil[kk], log.p = TRUE)
+            logphi.k <- dnorm(z + fij - fi[k], log = TRUE)
+            exp(logphi.k + logPhi.l) * dnorm(z)
+          }, lower = -Inf, upper = Inf
+        )$value)
+        ystar[i, k] <- fi[k] - exp(logD - logC)
+      }
+      ystar[i, j] <- fi[j] - sum(ystar[i, -j] - fi[-j])
+    }
+
+    # Update w -----------------------------------------------------------------
+    for (j in w.ind) {
+      A <- Hlamsq[[j]] + diag(1, n)
+      a <- as.numeric(crossprod(Hlam[[j]], ystar[, j] - alpha[j]))
+      eigenA <- iprior::eigenCpp(A)
+      V <- eigenA$vec
+      u <- eigenA$val + 1e-8  # ensure positive eigenvalues
+      uinv.Vt <- t(V) / u
+      w[, j] <- as.numeric(crossprod(a, V) %*% uinv.Vt)
+      Varw[[j]] <- iprior::fastVDiag(V, 1 / u)  # V %*% uinv.Vt
+      W[[j]] <- Varw[[j]] + tcrossprod(w[, j])
+      logdetA[j] <- determinant(A)$mod
+    }
+    if (isTRUE(common.intercept) && isTRUE(common.RKHS.scale)) {
+      w <- matrix(w[, 1], nrow = n, ncol = m)
+      W <- rep(list(W[[1]]), m)
+      logdetA <- rep(logdetA[1], m)
+    }
+
+    # Update theta -------------------------------------------------------------
+    count <- rep(0, m)
+    if (isTRUE(common.RKHS.scale)) {
+      # SAME RKHS PARAMETERS, SO JUST NEED TO SAMPLE ONCE AND REPEAT IN EACH CLASS
+      for (t in seq_len(n.samp - 1)) {
+        theta.star <- theta.current <- theta.samp[[t]]
+        theta.star[] <- rnorm(thetal$n.theta, mean = theta.current[, 1], sd = sd.samp)
+        Hlam.star <- get_Hlam(mod, theta.star)
+        Hlamsq.star <- lapply(Hlam.star, iprior::fastSquare)
+        log.q.star[] <- as.numeric(-0.5 * (
+          sum(Hlamsq.star[[1]] * W[[1]]) -
+            2 * crossprod(ystar[, 1] - alpha[1], Hlam.star[[1]]) %*% w[, 1]
+        ))
+        log.prob.acc <- log.q.star - log.q[t, ]
+        prob.acc <- exp(log.prob.acc)
+        prob.acc[prob.acc > 1] <- 1
+        if (runif(1) < prob.acc[j]) {
+          theta.samp[[t + 1]] <- theta.star
+          Hlaml[[t + 1]] <- Hlam.star
+          Hlamsql[[t + 1]] <- Hlamsq.star
+          log.q[t + 1, ] <- log.q.star
+          count <- count + 1
+        } else {
+          theta.samp[[t + 1]] <- theta.current
+          Hlaml[[t + 1]] <- Hlaml[[t]]
+          Hlamsql[[t + 1]] <- Hlamsql[[t]]
+          log.q[t + 1, ] <- log.q[t, ]
+        }
+      }
+    } else {
+      # DIFFERENT RKHS SCALE, SO INDEPENDENT SAMPLER FOR EACH CLASS
+      for (t in seq_len(n.samp - 1)) {
+        theta.star <- theta.current <- theta.samp[[t]]
+        for (j in 1:m) {
+          theta.star[, j] <- rnorm(thetal$n.theta, mean = theta.current[, j], sd = sd.samp)
+        }
+        Hlam.star <- get_Hlam(mod, theta.star)
+        Hlamsq.star <- lapply(Hlam.star, iprior::fastSquare)
+        for (j in 1:m) {
+          log.q.star[j] <- as.numeric(-0.5 * (
+            sum(Hlamsq.star[[j]] * W[[j]]) -
+              2 * crossprod(ystar[, j] - alpha[j], Hlam.star[[j]]) %*% w[, j]
+          ))
+        }
+        log.prob.acc <- log.q.star - log.q[t, ]
+        prob.acc <- exp(log.prob.acc)
+        prob.acc[prob.acc > 1] <- 1
+        Hlam.tmp <- Hlaml[[t]]
+        Hlamsq.tmp <- Hlamsql[[t]]
+        theta.tmp <- theta.current
+        for (j in 1:m) {
+          if (runif(1) < prob.acc[j]) {
+            theta.tmp[, j] <- theta.star[, j]
+            Hlam.tmp[[j]] <- Hlam.star[[j]]
+            Hlamsq.tmp[[j]] <- Hlamsq.star[[j]]
+            log.q[t + 1, j] <- log.q.star[j]
+            count[j] <- count[j] + 1
+          } else {
+            log.q[t + 1, j] <- log.q[t, j]
+          }
+        }
+        theta.samp[[t + 1]] <- theta.tmp
+        Hlaml[[t + 1]] <- Hlam.tmp
+        Hlamsql[[t + 1]] <- Hlamsq.tmp
+      }
+    }
+    acc.rate[niter, ] <- count / n.samp
+    theta <- Reduce("+", theta.samp[thin.seq]) / length(theta.samp[thin.seq])
+
+    # Update Hlam and Hlamsq ---------------------------------------------------
+    for (j in 1:m) {
+      tmp <- lapply(Hlaml, function(x) x[[j]])
+      Hlam[[j]] <- Reduce("+", tmp[thin.seq]) / length(tmp[thin.seq])
+      tmp <- lapply(Hlamsql, function(x) x[[j]])
+      Hlamsq[[j]] <- Reduce("+", tmp[thin.seq]) / length(tmp[thin.seq])
+    }
+
+    # Update alpha -------------------------------------------------------------
+    alpha <- apply(ystar - mapply("%*%", Hlam, split(w, col(w))), 2, mean)
+    if (isTRUE(common.intercept)) alpha <- rep(mean(alpha), m)
+
+    # Calculate lower bound ----------------------------------------------------
+    lb.ystar <- sum(logClb)
+    lb.w <- 0.5 * (nm - sum(sapply(W, function(x) sum(diag(x)))) - sum(logdetA))
+    lb.lambda <- -sum(apply(log.q, 2, mean))
+    if (isTRUE(common.intercept))
+      lb.alpha <- 0.5 * (1 + log(2 * pi) - log(nm))
+    else
+      lb.alpha <- (m / 2) * (1 + log(2 * pi) - log(n))
+    lb[niter + 1] <- lb.ystar + lb.w + lb.lambda + lb.alpha
+
+    # Calculate fitted values and error rate -----------------------------------
+    f.tmp <- rep(alpha, each = n) + mapply("%*%", Hlam, split(w, col(w)))
+    fitted.values <- probs_yhat_error(y, y.levels, f.tmp)
+    train.error[niter + 1] <- fitted.values$error
+    train.brier[niter + 1] <- fitted.values$brier
+    fitted.test <- NULL
+    if (iprior::.is.ipriorKernel_cv(mod)) {
+      ystar.test <- calc_ystar(mod, mod$Xl.test, alpha, theta, w)
+      fitted.test <- probs_yhat_error(y.test, y.levels, ystar.test)
+      test.error[niter + 1] <- fitted.test$error
+      test.brier[niter + 1] <- fitted.test$brier
+    }
+
+    niter <- niter + 1
+    if (!silent) setTxtProgressBar(pb, niter)
+  }
+
+  end.time <- Sys.time()
+  time.taken <- iprior::as.time(end.time - start.time)
+
+  # Calculate posterior s.d. and prepare param table ---------------------------
+  param.full <- theta_to_param.full(theta, alpha, mod)
+  if (isTRUE(common.intercept)) {
+    alpha <- unique(alpha)  # since they're all the same
+    se.alpha <- sqrt(1 / nm)
+  } else {
+    se.alpha <- rep(sqrt(1 / n), m)
+  }
+  alpha.summ <- cbind(alpha, se.alpha, alpha - qnorm(0.975) * se.alpha,
+                      alpha + qnorm(0.975) * se.alpha)
+  theta.summ <- theta.samp_to_summ_mult(theta.samp, mod)
+  colnames(alpha.summ) <- colnames(theta.summ)
+  param.summ <- rbind(alpha.summ, theta.summ)
+  rownames(param.summ) <- c(get_names(mod, "intercept", !common.intercept),
+                            get_names(mod, c("lambda", "hurst", "lengthscale",
+                                             "offset"), !common.RKHS.scale))
+
+  # Clean up and close ---------------------------------------------------------
+  convergence <- niter == maxit
+  if (!silent) {
+    close(pb)
+    if (convergence) cat("Convergence criterion not met.\n")
+    else cat("Converged after", niter, "iterations.\n")
+  }
+
+  list(theta = theta, param.full = param.full, param.summ = param.summ, w = w,
+       Varw = Varw, lower.bound = as.numeric(na.omit(lb)), niter = niter,
+       start.time = start.time, end.time = end.time, time = time.taken,
+       fitted.values = fitted.values, test = fitted.test,
+       train.error = as.numeric(na.omit(train.error)),
+       train.brier = as.numeric(na.omit(train.brier)),
+       test.error = as.numeric(na.omit(test.error)),
+       test.brier = as.numeric(na.omit(test.brier)), convergence = convergence,
+       acc.rate = as.numeric(na.omit(acc.rate)),
+       theta.samp = theta.samp[thin.seq])
+}
+
+theta.samp_to_summ_mult <- function(theta.samp, object) {
+  m <- ncol(theta.samp[[1]])
+  theta.samp_col <- function(x, j) x[, j]
+  res <- NULL
+  for (j in seq_len(m)) {
+    tmp <- t(sapply(theta.samp, theta.samp_col, j = j))
+    res[[j]] <- theta.samp_to_summ(tmp, object)
+  }
+  res.lambda <- do.call(rbind, lapply(res, function(x) x[grep("lambda", rownames(x)), ]))
+  res.hurst <- do.call(rbind, lapply(res, function(x) x[grep("hurst", rownames(x)), ]))
+  res.lengthscale <- do.call(rbind, lapply(res, function(x) x[grep("lengthscale", rownames(x)), ]))
+  res.offset <- do.call(rbind, lapply(res, function(x) x[grep("offset", rownames(x)), ]))
+
+  res <- rbind(
+    res.lambda, res.hurst, res.lengthscale, res.offset
+  )
+  unq.ind <- match(unique(res[, 1]), res[, 1])
+  # rownames(res) <- get_names(mod, c("lambda", "hurst", "lengthscale", "offset"))
+  res[unq.ind, ]
+}
+
+# else {
+#   # Nystrom approximation
+#   K.mm <- Hlamsq[[j]][1:Nystrom$m, 1:Nystrom$m]
+#   eigenK.mm <- iprior::eigenCpp(K.mm)
+#   V <- Hlamsq[[j]][, 1:Nystrom$m] %*% eigenK.mm$vec
+#   u <- eigenK.mm$val
+#   u.Vt <- t(V) * u
+#   D <- u.Vt %*% V + diag(1, Nystrom$m)
+#   E <- solve(D, u.Vt, tol = 1e-18)
+#   # see https://stackoverflow.com/questions/22134398/mahalonobis-distance-in-r-error-system-is-computationally-singular
+#   # see https://stackoverflow.com/questions/21451664/system-is-computationally-singular-error
+#   w[, j] <- as.numeric(a - V %*% (E %*% a))
+#   W[[j]] <- (diag(1, n) - V %*% E) + tcrossprod(w[, j])
+#   logdetA[j] <- determinant(A)$mod
+# }
diff --git a/man/Accessors.Rd b/man/Accessors.Rd
new file mode 100644
index 0000000..a296596
--- /dev/null
+++ b/man/Accessors.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Accessor_functions.R
+\name{Accessors}
+\alias{Accessors}
+\alias{get_hyperparam}
+\alias{get_intercept}
+\title{Accessor functions for \code{ipriorMod} objects.}
+\usage{
+get_hyperparam(object)
+
+get_intercept(object, by.class = FALSE)
+}
+\arguments{
+\item{object}{An \code{ipriorMod} object.}
+}
+\description{
+Accessor functions for \code{ipriorMod} objects.
+}
+\section{Functions}{
+\itemize{
+\item \code{get_hyperparam}: Obtain all of the hyperparameters.
+
+\item \code{get_intercept}: Obtain the intercept.
+}}
+
diff --git a/man/logLik.iprobitMod.Rd b/man/logLik.iprobitMod.Rd
index 77381f4..63e7eec 100644
--- a/man/logLik.iprobitMod.Rd
+++ b/man/logLik.iprobitMod.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/Utilities.R
+% Please edit documentation in R/Accessor_functions.R
 \name{logLik.iprobitMod}
 \alias{logLik.iprobitMod}
 \title{Extract the variational lower bound}
diff --git a/tests/testthat/test-5kernels.R b/tests/testthat/OLDtest-5kernels.R
similarity index 100%
rename from tests/testthat/test-5kernels.R
rename to tests/testthat/OLDtest-5kernels.R
diff --git a/tests/testthat/test-1gen_data.R b/tests/testthat/test-01-gen_data.R
similarity index 100%
rename from tests/testthat/test-1gen_data.R
rename to tests/testthat/test-01-gen_data.R
diff --git a/tests/testthat/test-02-iprobit_bin.R b/tests/testthat/test-02-iprobit_bin.R
new file mode 100644
index 0000000..05d9db8
--- /dev/null
+++ b/tests/testthat/test-02-iprobit_bin.R
@@ -0,0 +1,146 @@
+context("Binary estimation methods")
+
+test_that("VB-EM (Closed-form)", {
+
+  mod <- iprior::kernL(y ~ ., gen_mixture(n = 4, seed = 123), est.psi = FALSE)
+  mod1 <- iprobit(mod, control = list(silent = TRUE, maxit = 20, theta0 = 1:2,
+                                      alpha0 = 1))
+  mod2 <- iprobit_bin(mod, 20, silent = TRUE, alpha0 = 1, theta0 = 1:2)
+  expect_equal(as.numeric(coef(mod1)[, 1]),
+               c(-0.060995395, 0.975918045, 0.001124226), tol = 1e-5)
+  expect_equal(mod1$param.full, mod2$param.full, tol = 1e-5)
+
+})
+
+test_that("VB-EM with Metropolis sampler", {
+
+  mod <- iprior::kernL(y ~ ., gen_circle(n = 4, seed = 123), kernel = "se",
+                       est.lengthscale = TRUE, est.psi = FALSE)
+  mod1 <- iprobit(mod, control = list(silent = TRUE, maxit = 2, theta0 = 1:4,
+                                      alpha0 = 1, n.samp = 10, thin.samp = 1,
+                                      seed = 123))
+  mod2 <- iprobit_bin_metr(mod, 2, silent = TRUE, alpha0 = 1, theta0 = 1:4,
+                           n.samp = 10, thin.samp = 1, seed = 123)
+  expect_equal(as.numeric(coef(mod1)[, 1]),
+               c(0.1395819, 1.0103551, 2.0785975, 18.2429103, 83.1534269),
+               tol = 1e-5)
+  expect_equal(mod1$param.full, mod2$param.full, tol = 1e-5)
+
+})
+
+context("Binary models")
+
+# Generate data and model fit --------------------------------------------------
+dat <- gen_mixture(n = 10)
+dat.test <- gen_mixture(n = 5)
+
+# Regular
+mod <- iprobit(dat$y, dat$X, control = list(maxit = 5, silent = TRUE))
+# Formula (with training samples)
+modf <- iprobit(y ~ ., dat, one.lam = TRUE, train.samp = c(1, 5),
+                control = list(maxit = 5, silent = TRUE))
+# Metropolis
+modmetr <- iprobit(y ~ ., dat, kernel = "fbm", est.hurst = TRUE,
+                   control = list(maxit = 2, silent = TRUE))
+
+# Predict without quantiles ----------------------------------------------------
+mod.predict <- predict(mod, newdata = list(dat.test$X))
+modf.predict <- predict(modf, newdata = dat.test)
+modmetr.predict <- predict(modmetr, newdata = dat.test)
+
+# Predict with quantiles -------------------------------------------------------
+mod.predict.q <- predict(mod, list(dat.test$X), quantiles = TRUE, n.samp = 3)
+modf.predict.q <- predict(modf, dat.test, quantiles = TRUE, n.samp = 3)
+modmetr.predict.q <- predict(modmetr, dat.test,  quantiles = TRUE, n.samp = 3)
+
+test_that("Fitting binary models", {
+
+	expect_s3_class(mod, "iprobitMod")
+	expect_s3_class(mod, "iprobitMod_bin")
+	expect_s3_class(modf, "iprobitMod")
+	expect_s3_class(modf, "iprobitMod_bin")
+
+})
+
+test_that("Print", {
+
+  expect_that(print(mod), prints_text("Training error rate"))
+  expect_that(print(modf), prints_text("Training error rate"))
+
+})
+
+test_that("Summary", {
+
+  expect_s3_class(summary(mod), "iprobitMod_summary")
+  expect_s3_class(summary(modf), "iprobitMod_summary")
+
+})
+
+test_that("Fitted", {
+
+  expect_that(fitted(mod), is_a("iprobitPredict"))
+  expect_that(fitted(modf), is_a("iprobitPredict"))
+
+})
+
+test_that("Fitted quantiles", {
+
+  expect_that(fitted(mod, TRUE, n.samp = 3), is_a("iprobitPredict_quant"))
+  expect_that(fitted(modf, TRUE, n.samp = 3), is_a("iprobitPredict_quant"))
+
+})
+
+test_that("Predict (without test error rate)", {
+
+  expect_s3_class(mod.predict, "iprobitPredict")
+  expect_s3_class(modf.predict, "iprobitPredict")
+  expect_s3_class(modmetr.predict, "iprobitPredict")
+
+  expect_s3_class(mod.predict.q, "iprobitPredict_quant")
+  expect_s3_class(modf.predict.q, "iprobitPredict_quant")
+  expect_s3_class(modmetr.predict.q, "iprobitPredict_quant")
+
+  expect_that(print(mod.predict), prints_text("Test data not provided."))
+
+})
+
+test_that("Predict (with test error rate)", {
+
+  mod.predict <- predict(mod, newdata = list(dat.test$X), y = dat.test$y)
+
+  expect_that(print(mod.predict), prints_text("Test error"))
+  expect_that(print(modf.predict), prints_text("Test error"))
+  expect_that(print(modmetr.predict), prints_text("Test error"))
+
+})
+
+# test_that("Convergence", {
+#
+#   set.seed(123)
+#   dat <- gen_mixture(n = 10)
+#   mod <- iprobit(dat$y, dat$X, control = list(maxit = 500, silent = TRUE))
+#   modf <- iprobit(y ~ ., dat, one.lam = TRUE,
+#                   control = list(maxit = 500, silent = TRUE))
+#   expect_equal(as.numeric(get_lambda(mod)), 0.11646, tolerance = 1e-3)
+#   expect_equal(as.numeric(get_lambda(mod)), 0.11646, tolerance = 1e-3)
+#
+#   # > summary(mod)
+#   # Call:
+#   # iprobit(y = dat$y, X1 = dat$X, control = list(maxit = 500, silent = TRUE))
+#   #
+#   # Classes:
+#   # RKHS used:
+#   # Linear (X1)
+#   #
+#   # Hyperparameters:
+#   #           Mean   S.D.    2.5%  97.5%
+#   # alpha  -0.0090 0.3162 -0.6288 0.6108
+#   # lambda  0.1165 0.0201  0.0770 0.1559
+#   # ---
+#   #
+#   # Closed-form VB-EM algorithm. Iterations: 71/500
+#   # Converged to within 1e-05 tolerance. Time taken: 0.08387208 secs
+#   # Variational lower bound: -4.897003
+#   # Training error: 0%. Brier score: 0.004980669
+#
+# })
diff --git a/tests/testthat/test-03-iprobit_mult.R b/tests/testthat/test-03-iprobit_mult.R
new file mode 100644
index 0000000..40ca6d9
--- /dev/null
+++ b/tests/testthat/test-03-iprobit_mult.R
@@ -0,0 +1,222 @@
+context("Multinomial estimation methods")
+
+test_that("VB-EM (Closed-form)", {
+
+  mod <- iprior::kernL(y ~ ., gen_mixture(n = 6, m = 3, seed = 123),
+                       est.psi = FALSE)
+  theta0 <- matrix(rep(1, 6), ncol = 3)
+  alpha0 <- rep(1, 3)
+
+  # Different hyperparameters in each class ------------------------------------
+  mod1 <- iprobit(mod, control = list(silent = TRUE, maxit = 20, theta0 = theta0,
+                                      alpha0 = alpha0))
+  mod2 <- iprobit_mult(mod, 20, silent = TRUE, alpha0 = alpha0, theta0 = theta0)
+  expect_equal(as.numeric(coef(mod1)[, 1]),
+               c(6.502061e-01, 5.067932e-01, -6.667536e-05), tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 2]),
+               c(0.998459716, 0.001371398, 0.760741027), tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 3]),
+               c(1.3513342, 0.1209030, 0.6443124), tol = 1e-5)
+  expect_equal(mod1$param.full, mod2$param.full, tol = 1e-5)
+
+  # Only different RKHS parameters in each class -------------------------------
+  mod1 <- iprobit(mod, common.intercept = TRUE,
+                  control = list(silent = TRUE, maxit = 20, theta0 = theta0,
+                                 alpha0 = alpha0))
+  mod2 <- iprobit_mult(mod, 20, silent = TRUE, alpha0 = alpha0, theta0 = theta0,
+                       common.intercept = TRUE)
+  expect_equal(as.numeric(coef(mod1)[, 1]),
+               c(1, 0.4390309548, -0.0001031719), tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 2]),
+               c(1, 0.0009601449, 0.7745982109), tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 3]),
+               c(1, 0.1150476, 0.6924896), tol = 1e-5)
+  expect_equal(mod1$param.full, mod2$param.full, tol = 1e-5)
+
+  # Only different intercepts in each class ------------------------------------
+  mod1 <- iprobit(mod, common.RKHS.scale = TRUE,
+                  control = list(silent = TRUE, maxit = 20, theta0 = theta0,
+                                 alpha0 = alpha0))
+  mod2 <- iprobit_mult(mod, 20, silent = TRUE, alpha0 = alpha0, theta0 = theta0,
+                       common.RKHS.scale = TRUE)
+  expect_equal(as.numeric(coef(mod1)[, 1]),
+               c(0.61014, -0.005820041, 0.048168828), tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 2]),
+               c(0.99696, -0.005820041, 0.048168828), tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 3]),
+               c(1.39291, -0.005820041, 0.048168828), tol = 1e-5)
+  expect_equal(mod1$param.full, mod2$param.full, tol = 1e-5)
+
+  # Same hyperparameters in each class -----------------------------------------
+  mod1 <- iprobit(mod, common.intercept = TRUE, common.RKHS.scale = TRUE,
+                  control = list(silent = TRUE, maxit = 20, theta0 = theta0,
+                                 alpha0 = alpha0))
+  mod2 <- iprobit_mult(mod, 20, silent = TRUE, alpha0 = alpha0, theta0 = theta0,
+                       common.intercept = TRUE, common.RKHS.scale = TRUE)
+  expect_equal(as.numeric(coef(mod1)[, 1]),
+               c(1, 0 , 0), tol = 1e-5)
+  expect_equal(mod1$param.full, mod2$param.full, tol = 1e-5)
+
+})
+
+test_that("VB-EM with Metropolis sampler", {
+
+  mod <- iprior::kernL(y ~ ., gen_circle(n = 6, m = 3, seed = 123), kernel = "fbm",
+                       est.hurst = TRUE, est.psi = FALSE)
+  theta0 <- matrix(c(1, 1, 0, 0), ncol = 3, nrow = 4)
+  alpha0 <- rep(1, 3)
+
+  # Different hyperparameters in each class ------------------------------------
+  mod1 <- iprobit(mod, control = list(silent = TRUE, maxit = 20, theta0 = theta0,
+                                      alpha0 = alpha0, n.samp = 5,
+                                      thin.samp = 1, seed = 123))
+  mod2 <- iprobit_mult_metr(mod, 20, silent = TRUE, alpha0 = alpha0,
+                            theta0 = theta0, n.samp = 5, thin.samp = 1,
+                            seed = 123)
+  expect_equal(as.numeric(coef(mod1)[, 1]),
+               c(0.9985691, 1.0150100, 1.1014431, 0.5292299, 0.6026026),
+               tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 2]),
+               c(1.0055504, 0.6488293, 0.9076343, 0.5466233, 0.4670734),
+               tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 3]),
+               c(0.9958806, 0.9745103, 0.9592181, 0.4971903, 0.4359529),
+               tol = 1e-5)
+  expect_equal(mod1$param.full, mod2$param.full, tol = 1e-5)
+
+  # Only different RKHS parameters in each class -------------------------------
+  mod1 <- iprobit(mod, common.intercept = TRUE,
+                  control = list(silent = TRUE, maxit = 20, theta0 = theta0,
+                                 alpha0 = alpha0, n.samp = 5, thin.samp = 1,
+                                 seed = 123))
+  mod2 <- iprobit_mult_metr(mod, 20, silent = TRUE, alpha0 = alpha0,
+                            theta0 = theta0, n.samp = 5, thin.samp = 1,
+                            seed = 123, common.intercept = TRUE)
+  expect_equal(as.numeric(coef(mod1)[, 1]),
+               c(1, 1.0150100, 1.1014431, 0.5292299, 0.6026026), tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 2]),
+               c(1, 0.6488293, 0.9076343, 0.5466233, 0.4670734), tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 3]),
+               c(1, 0.9745103, 0.9592181, 0.4971903, 0.4359529), tol = 1e-5)
+  expect_equal(mod1$param.full, mod2$param.full, tol = 1e-5)
+
+  # Only different intercepts in each class ------------------------------------
+  mod1 <- iprobit(mod, common.RKHS.scale = TRUE,
+                  control = list(silent = TRUE, maxit = 20, theta0 = theta0,
+                                 alpha0 = alpha0, n.samp = 5, thin.samp = 1,
+                                 seed = 123))
+  mod2 <- iprobit_mult_metr(mod, 20, silent = TRUE, alpha0 = alpha0,
+                            theta0 = theta0, n.samp = 5, thin.samp = 1,
+                            seed = 123, common.RKHS.scale = TRUE)
+  expect_equal(as.numeric(coef(mod1)[, 1]),
+               c(0.9984867, 1.0390585, 1.0274495, 0.4885230, 0.5595608),
+               tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 2]),
+               c(1.0051488, 1.0390585, 1.0274495, 0.4885230, 0.5595608),
+               tol = 1e-5)
+  expect_equal(as.numeric(coef(mod1)[, 3]),
+               c(0.9963646, 1.0390585, 1.0274495, 0.4885230, 0.5595608),
+               tol = 1e-5)
+  expect_equal(mod1$param.full, mod2$param.full, tol = 1e-5)
+
+  # Same hyperparameters in each class -----------------------------------------
+  mod1 <- iprobit(mod, common.intercept = TRUE, common.RKHS.scale = TRUE,
+                  control = list(silent = TRUE, maxit = 20, theta0 = theta0,
+                                 alpha0 = alpha0, n.samp = 5, thin.samp = 1,
+                                 seed = 123))
+  mod2 <- iprobit_mult_metr(mod, 20, silent = TRUE, alpha0 = alpha0,
+                            theta0 = theta0, n.samp = 5, thin.samp = 1,
+                            seed = 123, common.RKHS.scale = TRUE,
+                            common.intercept = TRUE)
+  expect_equal(as.numeric(coef(mod1)[, 1]),
+               c(1, 1.0390585, 1.0274495, 0.4885230, 0.5595608), tol = 1e-5)
+  expect_equal(mod1$param.full, mod2$param.full, tol = 1e-5)
+
+})
+
+context("Multinomial models")
+
+# Generate data and model fit --------------------------------------------------
+set.seed(123)
+dat <- gen_circle(n = 15, m = 3)
+dat.test <- gen_circle(n = 6, m = 3)
+
+# Regular
+mod <- iprobit(dat$y, dat$X, control = list(maxit = 5, silent = TRUE))
+# Formula (with training samples)
+modf <- iprobit(y ~ ., dat, one.lam = TRUE, train.samp = c(1, 5),
+                control = list(maxit = 5, silent = TRUE))
+# Metropolis
+modmetr <- iprobit(y ~ ., dat, kernel = "fbm", est.hurst = TRUE,
+                   control = list(maxit = 2, silent = TRUE))
+
+# Predict without quantiles ----------------------------------------------------
+mod.predict <- predict(mod, newdata = list(dat.test$X))
+modf.predict <- predict(modf, newdata = dat.test)
+modmetr.predict <- predict(modmetr, newdata = dat.test)
+
+# Predict with quantiles -------------------------------------------------------
+mod.predict.q <- predict(mod, list(dat.test$X), quantiles = TRUE, n.samp = 3)
+modf.predict.q <- predict(modf, dat.test, quantiles = TRUE, n.samp = 3)
+modmetr.predict.q <- predict(modmetr, dat.test,  quantiles = TRUE, n.samp = 3)
+
+test_that("Fitting binary models", {
+
+  expect_s3_class(mod, "iprobitMod")
+  expect_s3_class(mod, "iprobitMod_mult")
+  expect_s3_class(modf, "iprobitMod")
+  expect_s3_class(modf, "iprobitMod_mult")
+
+})
+
+test_that("Print", {
+
+  expect_that(print(mod), prints_text("Training error rate"))
+  expect_that(print(modf), prints_text("Training error rate"))
+
+})
+
+test_that("Summary", {
+
+  expect_s3_class(summary(mod), "iprobitMod_summary")
+  expect_s3_class(summary(modf), "iprobitMod_summary")
+
+})
+
+test_that("Fitted", {
+
+  expect_that(fitted(mod), is_a("iprobitPredict"))
+  expect_that(fitted(modf), is_a("iprobitPredict"))
+
+})
+
+test_that("Fitted quantiles", {
+
+  expect_that(fitted(mod, TRUE, n.samp = 3), is_a("iprobitPredict_quant"))
+  expect_that(fitted(modf, TRUE, n.samp = 3), is_a("iprobitPredict_quant"))
+
+})
+
+test_that("Predict (without test error rate)", {
+
+  expect_s3_class(mod.predict, "iprobitPredict")
+  expect_s3_class(modf.predict, "iprobitPredict")
+  expect_s3_class(modmetr.predict, "iprobitPredict")
+
+  expect_s3_class(mod.predict.q, "iprobitPredict_quant")
+  expect_s3_class(modf.predict.q, "iprobitPredict_quant")
+  expect_s3_class(modmetr.predict.q, "iprobitPredict_quant")
+
+  expect_that(print(mod.predict), prints_text("Test data not provided."))
+
+})
+
+test_that("Predict (with test error rate)", {
+
+  mod.predict <- predict(mod, newdata = list(dat.test$X), y = dat.test$y)
+
+  expect_that(print(mod.predict), prints_text("Test error"))
+  expect_that(print(modf.predict), prints_text("Test error"))
+  expect_that(print(modmetr.predict), prints_text("Test error"))
+
+})
diff --git a/tests/testthat/test-iprobit_helpers.R b/tests/testthat/test-04-iprobit_helpers.R
similarity index 100%
rename from tests/testthat/test-iprobit_helpers.R
rename to tests/testthat/test-04-iprobit_helpers.R
diff --git a/tests/testthat/test-05-plots.R b/tests/testthat/test-05-plots.R
new file mode 100644
index 0000000..f1c58f7
--- /dev/null
+++ b/tests/testthat/test-05-plots.R
@@ -0,0 +1,164 @@
+context("Plots")
+
+# Binary models ----------------------------------------------------------------
+
+set.seed(123)
+n <- 50
+dat <- gen_circle(n)
+mod <- iprobit(y ~ X1 + X2, dat, one.lam = TRUE, kernel = "fbm",
+               train.samp = sort(sample(seq_len(n), size = 40)),
+               control = list(silent = TRUE, maxit = 10))
+
+test_that("iplot_fitted()", {
+
+  expect_silent(p <- iplot_fitted(mod))
+
+})
+
+test_that("iplot_dec_bound()", {
+
+  expect_silent(p <- iplot_dec_bound(mod, grid.len = 5))
+
+})
+
+test_that("iplot_predict()", {
+
+  expect_silent(p <- iplot_predict(mod, grid.len = 5))
+  expect_silent(p <- iplot_predict(mod, grid.len = 5, dec.bound = FALSE))
+  expect_silent(p <- iplot_predict(mod, grid.len = 5, plot.test = FALSE))
+  expect_silent(p <- iplot_predict(mod, grid.len = 5, dec.bound = FALSE,
+                                   plot.test = FALSE))
+
+})
+
+test_that("iplot_lb()", {
+
+  expect_silent(p <- iplot_lb(mod))
+
+})
+
+test_that("iplot_error()", {
+
+  expect_silent(p <- iplot_error(mod, 3))
+  expect_silent(p <- iplot_error(mod, plot.test = FALSE))
+
+})
+
+test_that("iplot_lb_and_error()", {
+
+  expect_silent(p <- iplot_lb_and_error(mod, niter.plot = 3, plot.test = FALSE))
+
+})
+
+# Binary models ----------------------------------------------------------------
+
+set.seed(123)
+n <- 100
+dat <- gen_mixture(n, m = 3)
+mod <- iprobit(y ~ ., dat, kernel = "fbm", train.samp = sort(sample(seq_len(n),
+                                                                    size = 90)),
+               control = list(silent = TRUE, maxit = 4))
+
+test_that("iplot_fitted()", {
+
+  expect_silent(p <- iplot_fitted(mod))
+
+})
+
+test_that("iplot_dec_bound()", {
+
+  expect_silent(p <- iplot_dec_bound(mod, grid.len = 5))
+
+})
+
+test_that("iplot_predict()", {
+
+  expect_silent(p <- iplot_predict(mod, grid.len = 5))
+  expect_silent(p <- iplot_predict(mod, grid.len = 5, dec.bound = FALSE))
+  expect_silent(p <- iplot_predict(mod, grid.len = 5, plot.test = FALSE))
+  expect_silent(p <- iplot_predict(mod, grid.len = 5, dec.bound = FALSE,
+                                   plot.test = FALSE))
+
+})
+
+test_that("iplot_lb()", {
+
+  expect_silent(p <- iplot_lb(mod))
+
+})
+
+test_that("iplot_error()", {
+
+  expect_silent(p <- iplot_error(mod, 3))
+  expect_silent(p <- iplot_error(mod, plot.test = FALSE))
+
+})
+
+test_that("iplot_lb_and_error()", {
+
+  expect_silent(p <- iplot_lb_and_error(mod, niter.plot = 3, plot.test = FALSE))
+
+})
+
+# context("Methods for iprobitMod objects")
+
+# Update/refit -----------------------------------------------------------------
+#
+# test_that("Update from formula", {
+#
+#   dat <- gen_mixture(10)
+#   mod.original <- mod <- iprobit(y ~ ., dat, control = list(maxit = 1),
+#                                  silent = TRUE)
+#   mod <- iprobit(mod, maxit = 5)
+#
+#   expect_s3_class(mod, "iprobitMod")
+#   expect_equal(mod$niter, 6)
+#   expect_equal(length(mod$lower.bound), 6)
+#   expect_equal(length(mod$error), 6)
+#   expect_equal(length(mod$brier), 6)
+#   expect_equal(mod$formula, mod.original$formula)
+#
+# })
+#
+# test_that("Update from non-formula", {
+#
+#   dat <- gen_mixture(10)
+#   mod.original <- mod <- iprobit(dat$y, dat$X, control = list(maxit = 1),
+#                                  silent = TRUE)
+#   mod <- iprobit(mod, maxit = 5)
+#
+#   expect_s3_class(mod, "iprobitMod")
+#   expect_equal(mod$niter, 6)
+#   expect_equal(length(mod$lower.bound), 6)
+#   expect_equal(length(mod$error), 6)
+#   expect_equal(length(mod$brier), 6)
+#   expect_equal(mod$formula, mod.original$formula)
+#
+# })
+#
+# test_that("Update options", {
+#
+#   dat <- gen_mixture(10)
+#   mod.original <- mod <- iprobit(dat$y, dat$X, control = list(maxit = 1),
+#                                  silent = TRUE)
+#
+#   expect_message(mod <- iprobit(mod))
+#   expect_output(mod <- iprobit(mod, stop.crit = 1e-1, maxit = 5, silent = FALSE))
+#
+# })
+#
+# test_that("The function update()", {
+#
+#   dat <- gen_mixture(10)
+#   mod.original <- mod <- iprobit(dat$y, dat$X, control = list(maxit = 1),
+#                                  silent = TRUE)
+#   update(mod, maxit = 5)
+#
+#   expect_s3_class(mod, "iprobitMod")
+#   expect_equal(mod$niter, 6)
+#   expect_equal(length(mod$lower.bound), 6)
+#   expect_equal(length(mod$error), 6)
+#   expect_equal(length(mod$brier), 6)
+#   expect_equal(mod$formula, mod.original$formula)
+#
+# })
diff --git a/tests/testthat/test-7methods.R b/tests/testthat/test-06-accessors.R
similarity index 95%
rename from tests/testthat/test-7methods.R
rename to tests/testthat/test-06-accessors.R
index 6dcf980..a5951db 100644
--- a/tests/testthat/test-7methods.R
+++ b/tests/testthat/test-06-accessors.R
@@ -1,6 +1,6 @@
 # context("Methods for iprobitMod objects")
-#
-# # Update/refit -----------------------------------------------------------------
+
+# Update/refit -----------------------------------------------------------------
 #
 # test_that("Update from formula", {
 #
diff --git a/tests/testthat/test-6utilities.R b/tests/testthat/test-07-utilities.R
similarity index 100%
rename from tests/testthat/test-6utilities.R
rename to tests/testthat/test-07-utilities.R
diff --git a/tests/testthat/test-4errors.R b/tests/testthat/test-08-errors.R
similarity index 100%
rename from tests/testthat/test-4errors.R
rename to tests/testthat/test-08-errors.R
diff --git a/tests/testthat/test-8parallel.R b/tests/testthat/test-09-parallel.R
similarity index 100%
rename from tests/testthat/test-8parallel.R
rename to tests/testthat/test-09-parallel.R
diff --git a/tests/testthat/test-2iprobit_bin.R b/tests/testthat/test-2iprobit_bin.R
deleted file mode 100644
index e162211..0000000
--- a/tests/testthat/test-2iprobit_bin.R
+++ /dev/null
@@ -1,105 +0,0 @@
-context("Binary models")
-
-dat <- gen_mixture(n = 10)
-dat.test <- gen_mixture(n = 5)
-
-mod <- iprobit(dat$y, dat$X, control = list(maxit = 5, silent = TRUE))
-modf <- iprobit(y ~ ., dat, one.lam = TRUE, control = list(maxit = 5,
-                                                           silent = TRUE))
-
-dat.test <- gen_mixture(n = 5)
-mod.predict <- predict(mod, newdata = list(dat.test$X))
-modf.predict <- predict(modf, newdata = dat.test)
-
-mod.predict.q <- predict(mod, newdata = list(dat.test$X), quantiles = TRUE,
-                         n.samp = 3)
-modf.predict.q <- predict(modf, newdata = dat.test, quantiles = TRUE,
-                          n.samp = 3)
-
-
-test_that("Fitting binary models", {
-
-	expect_s3_class(mod, "iprobitMod")
-	expect_s3_class(mod, "iprobitMod_bin")
-	expect_s3_class(modf, "iprobitMod")
-	expect_s3_class(modf, "iprobitMod_bin")
-
-})
-
-test_that("Print", {
-
-  expect_that(print(mod), prints_text("Training error rate"))
-  expect_that(print(modf), prints_text("Training error rate"))
-
-})
-
-test_that("Summary", {
-
-  expect_s3_class(summary(mod), "iprobitMod_summary")
-  expect_s3_class(summary(modf), "iprobitMod_summary")
-
-})
-
-test_that("Fitted", {
-
-  expect_that(fitted(mod), is_a("iprobitPredict"))
-  expect_that(fitted(modf), is_a("iprobitPredict"))
-
-})
-
-test_that("Fitted quantiles", {
-
-  expect_that(fitted(mod, TRUE, n.samp = 3), is_a("iprobitPredict_quant"))
-  expect_that(fitted(modf, TRUE, n.samp = 3), is_a("iprobitPredict_quant"))
-
-})
-
-test_that("Predict (without test error rate)", {
-
-  expect_s3_class(mod.predict, "iprobitPredict")
-  expect_s3_class(modf.predict, "iprobitPredict")
-  expect_s3_class(mod.predict.q, "iprobitPredict_quant")
-  expect_s3_class(modf.predict.q, "iprobitPredict_quant")
-  expect_that(print(mod.predict), prints_text("Test data not provided."))
-
-})
-
-test_that("Predict (with test error rate)", {
-
-  mod.predict <- predict(mod, newdata = list(dat.test$X), y = dat.test$y)
-
-  expect_that(print(mod.predict), prints_text("Test error"))
-  expect_that(print(modf.predict), prints_text("Test error"))
-
-})
-
-test_that("Convergence", {
-
-  set.seed(123)
-  dat <- gen_mixture(n = 10)
-  mod <- iprobit(dat$y, dat$X, control = list(maxit = 500, silent = TRUE))
-  modf <- iprobit(y ~ ., dat, one.lam = TRUE,
-                  control = list(maxit = 500, silent = TRUE))
-  expect_equal(as.numeric(get_lambda(mod)), 0.11646, tolerance = 1e-3)
-  expect_equal(as.numeric(get_lambda(mod)), 0.11646, tolerance = 1e-3)
-
-  # > summary(mod)
-  # Call:
-  # iprobit(y = dat$y, X1 = dat$X, control = list(maxit = 500, silent = TRUE))
-  #
-  # Classes:
-  # RKHS used:
-  # Linear (X1)
-  #
-  # Hyperparameters:
-  #           Mean   S.D.    2.5%  97.5%
-  # alpha  -0.0090 0.3162 -0.6288 0.6108
-  # lambda  0.1165 0.0201  0.0770 0.1559
-  # ---
-  #
-  # Closed-form VB-EM algorithm. Iterations: 71/500
-  # Converged to within 1e-05 tolerance. Time taken: 0.08387208 secs
-  # Variational lower bound: -4.897003
-  # Training error: 0%. Brier score: 0.004980669
-
-})
diff --git a/tests/testthat/test-3iprobit_mult.R b/tests/testthat/test-3iprobit_mult.R
deleted file mode 100644
index 27e4e7c..0000000
--- a/tests/testthat/test-3iprobit_mult.R
+++ /dev/null
@@ -1,142 +0,0 @@
-# context("Multinomial models")
-#
-# test_that("Fitting multinomial models (IIA)", {
-#
-#   m <- 4
-# 	dat <- gen_mixture(n = 10, m = m)
-# 	mod1 <- iprobit(dat$y, dat$X, silent = TRUE, control = list(maxit = 5))
-# 	modf <- iprobit(y ~ ., dat, silent = TRUE, one.lam = TRUE,
-# 	                control = list(maxit = 5))
-# 	expect_s3_class(mod1, "iprobitMod")
-# 	expect_s3_class(mod1, "iprobitMod_mult")
-# 	expect_s3_class(modf, "iprobitMod")
-# 	expect_s3_class(modf, "iprobitMod_mult")
-#
-# 	# Common intercept
-# 	mod2 <- iprobit(dat$y, dat$X, silent = TRUE,
-# 	                control = list(maxit = 5, common.intercept = TRUE))
-# 	expect_true(all.same(get_alpha(mod2)))
-#
-# 	# Common RKHS scale parameter
-# 	mod3 <- iprobit(dat$y, dat$X, silent = TRUE,
-# 	                control = list(maxit = 5, common.RKHS.scale = TRUE))
-# 	expect_true(all.same(get_lambda(mod3)))
-#
-# 	# Common intercept and RKHS scale parameter
-# 	mod4 <- iprobit(dat$y, dat$X, silent = TRUE,
-# 	                control = list(maxit = 5,
-# 	                               common.intercept = TRUE,
-# 	                               common.RKHS.scale = TRUE))
-# 	expect_true(all.same(get_alpha(mod4)))
-# 	expect_true(all.same(get_lambda(mod4)))
-#
-# })
-#
-# test_that("Print (IIA)", {
-#
-#   m <- 4
-#   dat <- gen_mixture(n = 10, m = m)
-#   mod <- iprobit(dat$y, dat$X, silent = TRUE, control = list(maxit = 5))
-#   modf <- iprobit(y ~ ., dat, silent = TRUE, one.lam = TRUE, control = list(maxit = 5))
-#   expect_that(print(mod), prints_text("Training error rate"))
-#   expect_that(print(modf), prints_text("Training error rate"))
-#
-# })
-#
-# test_that("Summary", {
-#
-#   m <- 4
-#   dat <- gen_mixture(n = 10, m = m)
-#   mod <- iprobit(dat$y, dat$X, silent = TRUE, control = list(maxit = 5))
-#   modf <- iprobit(y ~ ., dat, silent = TRUE, one.lam = TRUE, control = list(maxit = 5))
-#   mod.summary <- summary(mod)
-#   modf.summary <- summary(modf)
-#   expect_s3_class(mod.summary, "iprobitSummary")
-#   expect_s3_class(modf.summary, "iprobitSummary")
-#
-# })
-#
-# test_that("Fitted", {
-#
-#   m <- 4
-#   dat <- gen_mixture(n = 10, m = m)
-#   mod <- iprobit(dat$y, dat$X, silent = TRUE, control = list(maxit = 5))
-#   modf <- iprobit(y ~ ., dat, silent = TRUE, one.lam = TRUE, control = list(maxit = 5))
-#   expect_that(fitted(mod), is_a("iprobit_predict"))
-#   expect_that(fitted(modf), is_a("iprobit_predict"))
-#
-# })
-#
-# test_that("Fitted quantiles", {
-#
-#   m <- 4
-#   dat <- gen_mixture(n = 10, m = m)
-#   mod <- iprobit(dat$y, dat$X, silent = TRUE, control = list(maxit = 5))
-#   modf <- iprobit(y ~ ., dat, silent = TRUE, one.lam = TRUE, control = list(maxit = 5))
-#   expect_that(fitted(mod, TRUE, n.samp = 3), is_a("iprobit_predict_quant"))
-#   expect_that(fitted(modf, TRUE, n.samp = 3), is_a("iprobit_predict_quant"))
-#
-# })
-#
-# test_that("Predict (without test error rate)", {
-#
-#   m <- 4
-#   dat <- gen_mixture(n = 10, m = m)
-#   dat.test <- gen_mixture(n = 5, m = m)
-#   mod <- iprobit(dat$y, dat$X, silent = TRUE, control = list(maxit = 5))
-#   modf <- iprobit(y ~ ., dat, silent = TRUE, one.lam = TRUE,
-#                   control = list(maxit = 5))
-#
-#   mod.predict <- predict(mod, newdata = list(dat.test$X))
-#   modf.predict <- predict(modf, newdata = dat.test)
-#
-#   expect_s3_class(mod.predict, "iprobit_predict")
-#   expect_s3_class(modf.predict, "iprobit_predict")
-#   expect_that(print(mod.predict), prints_text("Test data not provided."))
-#
-# })
-#
-# test_that("Predict (with test error rate)", {
-#
-#   m <- 4
-#   dat <- gen_mixture(n = 10, m = m)
-#   dat.test <- gen_mixture(n = 5, m = m)
-#   mod <- iprobit(dat$y, dat$X, silent = TRUE, control = list(maxit = 5))
-#   modf <- iprobit(y ~ ., dat, silent = TRUE, one.lam = TRUE,
-#                   control = list(maxit = 5))
-#
-#   mod.predict <- predict(mod, newdata = list(dat.test$X), y = dat.test$y)
-#   modf.predict <- predict(modf, newdata = dat.test)
-#
-#   expect_that(print(mod.predict), prints_text("Test error"))
-#   expect_that(print(modf.predict), prints_text("Test error"))
-#
-# })
-#
-# test_that("Convergence", {
-#
-#   set.seed(123)
-#   m <- 3
-#   dat <- gen_mixture(n = 10, m = m)
-#   mod <- iprobit(dat$y, dat$X,
-#                  control = list(maxit = 500, silent = TRUE,
-#                                 alpha0 = rep(1, m),
-#                                 lambda0 = rep(1, m)))
-#   modf <- iprobit(y ~ ., dat, one.lam = TRUE,
-#                   control = list(maxit = 500, silent = TRUE,
-#                                  alpha0 = rep(1, m),
-#                                  lambda0 = rep(1, m)))
-#   expect_equal(get_lambda(mod)[2], 0.1743832, tolerance = 1e-3)
-#   expect_equal(get_lambda(modf)[2], 0.1743832, tolerance = 1e-3)
-#
-#   # Single lambda
-#   # > modf
-#   # Training error rate: 20.00 %
-#   # Lower bound value: -12.87285
-#   #
-#   # Class = 1 Class = 2 Class = 3
-#   # alpha    0.60748   1.76979   0.62274
-#   # lambda   0.00000   0.17438   0.00000
-#
-# })
-#