diff --git a/NAMESPACE b/NAMESPACE
index 207bee5..94ef069 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -26,6 +26,7 @@ S3method(summary,group.lasso)
 S3method(summary,netReg.family)
 export(beta)
 export(binomial)
+export(bum)
 export(cv.edgenet)
 export(edgenet)
 export(family)
@@ -49,7 +50,10 @@ importFrom(tensorflow,install_tensorflow)
 importFrom(tensorflow,tf)
 importFrom(tfprobability,tfd_bernoulli)
 importFrom(tfprobability,tfd_beta)
+importFrom(tfprobability,tfd_categorical)
 importFrom(tfprobability,tfd_gamma)
 importFrom(tfprobability,tfd_inverse_gaussian)
+importFrom(tfprobability,tfd_mixture)
 importFrom(tfprobability,tfd_poisson)
+importFrom(tfprobability,tfd_uniform)
 useDynLib(netReg, .registration = TRUE)
diff --git a/R/family.R b/R/family.R
index 05b83f3..b8c07c9 100644
--- a/R/family.R
+++ b/R/family.R
@@ -161,6 +161,28 @@ beta <- function(link = c("logit", "probit", "log")) {
 }
 
 
+#' @export
+#' @rdname family-methods
+bum <- function(link = c("logit", "probit", "log")) {
+  warn.experimental("bum")
+  link <- match.arg(link)
+  linkinv <- switch(
+    link,
+    "logit" = logistic,
+    "log" = exp,
+    "probit" = gcdf,
+    stop("did not recognize link function", call. = FALSE)
+  )
+
+  .as.family(
+    "bum",
+    link,
+    linkinv,
+    bum.loss
+  )
+}
+
+
 #' @export
 #' @rdname family-methods
 inverse.gaussian <- function(link = c("1/mu^2")) {
diff --git a/R/likelihood.R b/R/likelihood.R
index 22ab4d2..f85a7d7 100644
--- a/R/likelihood.R
+++ b/R/likelihood.R
@@ -98,6 +98,54 @@ beta.loss <- function(y, mu.hat, ...) {
 }
 
 
+#' @noRd
+#' @importFrom tensorflow tf
+#' @importFrom tfprobability tfd_beta tfd_uniform tfd_mixture tfd_categorical
+bum.loss <- function(y, mu.hat, ...) {
+  obj <- 0
+  eps <- .Machine$double.eps * 1e9
+  N <- mu.hat$shape[[1]]
+  for (j in seq(ncol(y))) {
+    mu <- mu.hat[, j]
+    mixture_weight <- tf$constant(matrix(0.75, nrow = N, ncol = 2), dtype = tf$float32) # TODO: learn this parameter
+    phi <- 1 # TODO: replace this with tf$variable
+
+    # reparametrize
+    # concentration1 := alpha = mu * phi
+    p <- mu * phi
+    # concentration0 := beta = (1. - mu) * phi
+    q <- (1 - mu) * phi
+
+    # deal with numerical instabilities
+    p.trans <- tf$math$maximum(p, eps)
+    q.trans <- tf$math$maximum(q, eps)
+
+    # need correct batch dimensions for mixture
+    p.trans <- tf$stack(list(p.trans, p.trans), 0L)
+    q.trans <- tf$stack(list(q.trans, q.trans), 0L)
+
+    prob <- tfprobability::tfd_mixture(
+      cat = tfprobability::tfd_categorical(
+        probs = c(mixture_weight, 1 - mixture_weight)
+      ),
+      components = c(
+        tfprobability::tfd_beta(
+          concentration1 = p.trans, concentration0 = q.trans
+        ),
+        tfprobability::tfd_uniform(
+          low = tf$constant(tf$Variable(matrix(0, nrow = 2, ncol = N), dtype=tf$float32)),
+          high = tf$constant(tf$Variable(matrix(1, nrow = 2, ncol = N), dtype=tf$float32))
+        )
+      )
+    )
+
+    obj <- obj + tf$reduce_sum(prob$log_prob(y[, j]))
+  }
+
+  -obj
+}
+
+
 #' @noRd
 #' @importFrom tensorflow tf
 #' @importFrom tfprobability tfd_inverse_gaussian
diff --git a/man/family-methods.Rd b/man/family-methods.Rd
index a1bbae8..599d26b 100644
--- a/man/family-methods.Rd
+++ b/man/family-methods.Rd
@@ -8,6 +8,7 @@
 \alias{poisson}
 \alias{gamma}
 \alias{beta}
+\alias{bum}
 \alias{inverse.gaussian}
 \title{Family objects for models}
 \usage{
@@ -23,6 +24,8 @@ gamma(link = c("inverse"))
 
 beta(link = c("logit", "probit", "log"))
 
+bum(link = c("logit", "probit", "log"))
+
 inverse.gaussian(link = c("1/mu^2"))
 }
 \arguments{