Log likelihood for `brms` custom family that depends on latent variables?

I’m implementing a brms custom family and would like to add a log_lik method as described in the “Post-Processing Custom Family Models” section of the vignette “Define Custom Response Distributions with brms”.

The challenge I am having is that the log likelihood \log p(y | \theta) depends on some “latent variables” (not only the distributional parameters of the model, which can be obtained via brms::get_dpar). I am unsure as to how I can extract these latent variable parameters from the model in order to compute the log likelihood.

Here is what I have written for the log_lik_latent_lognormal function so far (where latent_lognormal is the name of the family):

#' Calculate the pointwise log likelihood
#'
#' See [brms::log_lik()].
#'
#' @param i The index of the observation to calculate the log likelihood of
#' @param prep The result of a call to [`brms::prepare_predictions`]
#' @family postprocess
#' @autoglobal
#' @export
log_lik_latent_lognormal <- function(i, prep) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  sigma <- brms::get_dpar(prep, "sigma", i = i)
  y <- prep$data$Y[i]
  obs_t <- prep$data$vreal1[i]
  pwindow_width <- prep$data$vreal2[i]
  swindow_width <- prep$data$vreal3[i]

  # Generate values of the swindow_raw and pwindow_raw
  # Really these should be extracted from prep
  # How to do this?
  swindow_raw <- runif(prep$ndraws)
  pwindow_raw <- runif(prep$ndraws)

  pwindow <- pwindow_raw * pwindow_width
  swindow <- swindow_raw * swindow_width

  d <- y - pwindow + swindow
  obs_time <- obs_t - pwindow
  lpdf <- dlnorm(d, meanlog = mu, sdlog = sigma, log = TRUE)
  lcdf <- plnorm(obs_time, meanlog = mu, sdlog = sigma, log.p = TRUE)
  return(lpdf - lcdf)
}

The section with the comment is where I’d like to be getting the pwindow_raw and swindow_raw draws from the fitted model (via prep somehow).

I’ve looked a bit a the brms code (see issue for a log of my thinking [including at one point I think wrongly thinking about integrating these variables out to get a marginal likelihood]). Anyway, in brms:::prepare_predictions.bframel there are these prepare_predictions_xx functions for a range of xx. Are these like special cases where more than the distributional parameters are required? Would my situation be like one of these special cases?

    out$fe <- prepare_predictions_fe(x, draws, sdata, ...)
    out$sp <- prepare_predictions_sp(x, draws, sdata, ...)
    out$cs <- prepare_predictions_cs(x, draws, sdata, ...)
    out$sm <- prepare_predictions_sm(x, draws, sdata, ...)
    out$gp <- prepare_predictions_gp(x, draws, sdata, ...)
    out$re <- prepare_predictions_re(x, sdata, ...)
    out$ac <- prepare_predictions_ac(x, draws, sdata, nat_cov = FALSE, 

Another idea: is it possible it put within ... the fit object, and then rely on the user providing this to any functions requiring the log_lik?

Thank you! Let me know if I could make any part of this more clear.

Session information

  • Operating System: Mac OS
  • brms version: 2.21.6