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)