Getting around truncation limitation for mixture models in brms

Interfaces brms
3

Thanks, @paul.buerkner. Following the fix on the GitHub version, I was able to implement this. Here’s a working example:

library(dplyr)
library(brms) #using GitHub brms, following fix of #453

#make some data
dat <- data_frame(response = c(rnorm(200, 3, 1), rnorm(200, 6, 1))) %>%
  filter(response < 8)

# Define a distribution of a normal with a truncated upper bound
normal_ub <- custom_family("normal_ub", dpars = c("mu", "sigma", "upper"),
                           links = c("identity","identity","identity"),
                           type = "real",
                           lb = c(NA, 0, NA), ub = c(NA, NA, NA))
stan_funs <- "
real normal_ub_lpdf(real y, real mu, real sigma, real upper) {
return normal_lpdf(y | mu, sigma) - normal_lcdf(upper | mu, sigma);
}
real normal_ub_rng(real mu, real sigma, real upper) {
real p = normal_cdf(upper, mu, sigma);  // cdf for bounds
real u = uniform_rng(0, p);
return inv_Phi(u)*sigma + mu;  // inverse cdf for value
}
"

mix <- mixture(gaussian, normal_ub, order = 'mu')
prior <- c(
  prior(normal(3, 1), class="Intercept", dpar = "mu1"),
  prior(normal(6, 1), class="Intercept", dpar = "mu2"),
  prior(normal(8.1, 0.0001), class="upper2") #use get_prior() to find the name
)

trunc_mod <- brm(bf(response ~ 1), data = dat, family = mix,
                 stan_funs = stan_funs, prior = prior,
                 iter = 1000, warmup = 500, chains = 1,
                 control = list(adapt_delta = 0.999, max_treedepth=12),
                 save_model = "test_truncated.stan",
                 save_dso = TRUE,
                 seed = 89)

#Get our distribution functions to use post-fitting
expose_functions(trunc_mod, vectorize = TRUE)

log_lik_normal_ub <- function(i, draws) {
  mu <- draws$dpars$mu[, i]
  sigma <- draws$dpars$sigma
  upper <- draws$dpars$upper
  y <- draws$data$Y[i]
  normal_ub_lpmf(y, mu, sigma, upper)
}

predict_normal_ub <- function(i, draws, ...) {
  mu <- draws$dpars$mu[, i]
  sigma <- draws$dpars$sigma
  upper <- draws$dpars$upper
  normal_ub_rng(mu, sigma, upper)
}

pp_check(trunc_mod)
1 Like