Errors when right-censoring: bug in normal_lccdf?

I've been running into issues fitting models with a right-censored normal distribution, using brms. I've submitted an issue with brms here: https://github.com/paul-buerkner/brms/issues/1423 However, it appears that the error does not arise from brms. @paul-buerkner suggested it might instead arise from a bug in normal_lccdf, so I thought I'd post it here.

I appears that the issue occurs only

• using the gaussian family; others such as student or weibull work fine;
• when right-censoring, whereas left-censoring works fine (suggesting normal_lccdf as the culprit);
• quite recently (as others report running such models a few months back without issue).

Below is my original post to the brms issue tracker, containing a MWE. Apologies for not testing this with other Stan interfaces.

# Censoring example
set.seed(25)
n <- 5000
t1 <- 106 # Right censoring cutoff
t2 <- 94 # Left censoring cutoff

d <- tibble(y = rnorm(n, mean = 100, sd = 15),
            x = rnorm(n, 0, 1)) %>% 
  # Right censored variables
  mutate(y1   = if_else(y > t1, t1, y),
         cen1 = if_else(y > t1, "right", "none")) %>%
  # Left censored variables
  mutate(y2 = if_else(y < t2, t2, y),
         cen2 = if_else(y < t2, "left", "none"))

# Right censoring
fit <- brm(bf(y1 | cens(cen1) ~ x),
    data = d)

Compiling Stan program...
Start sampling
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
...
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
[1] "error occurred during calling the sampler; sampling not done"
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: y1 | cens(cen1) ~ x 
   Data: d (Number of observations: 5000) 

This occurs across multiple datasets and different kinds of models (linear mixed model with one random factor; mixed-effects location scale model), and when init = 0.

Stan code for this model (generated with brms 2.17.0):

functions {
}
data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
  int<lower=-1,upper=2> cens[N];  // indicates censoring
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  int Kc = K - 1;
  matrix[N, Kc] Xc;  // centered version of X without an intercept
  vector[Kc] means_X;  // column means of X before centering
  for (i in 2:K) {
    means_X[i - 1] = mean(X[, i]);
    Xc[, i - 1] = X[, i] - means_X[i - 1];
  }
}
parameters {
  vector[Kc] b;  // population-level effects
  real Intercept;  // temporary intercept for centered predictors
  real<lower=0> sigma;  // dispersion parameter
}
transformed parameters {
  real lprior = 0;  // prior contributions to the log posterior
  lprior += student_t_lpdf(Intercept | 3, 99.9, 9);
  lprior += student_t_lpdf(sigma | 3, 0, 9)
    - 1 * student_t_lccdf(0 | 3, 0, 9);
}
model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N] mu = Intercept + Xc * b;
    for (n in 1:N) {
    // special treatment of censored data
      if (cens[n] == 0) {
        target += normal_lpdf(Y[n] | mu[n], sigma);
      } else if (cens[n] == 1) {
        target += normal_lccdf(Y[n] | mu[n], sigma);
      } else if (cens[n] == -1) {
        target += normal_lcdf(Y[n] | mu[n], sigma);
      }
    }
  }
  // priors including constants
  target += lprior;
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept - dot_product(means_X, b);
}

In contrast, left-censoring works:

# Left censoring
brm(bf(y2 | cens(cen2) ~ x),
    data = d)

I tested this on two machines.

# Machine 1
> sessionInfo()
R version 4.0.4 (2021-02-15)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] brms_2.17.0     Rcpp_1.0.9      forcats_0.5.1   stringr_1.4.0   dplyr_1.0.6     purrr_0.3.4     readr_1.4.0     tidyr_1.1.3     tibble_3.1.8   
[10] ggplot2_3.3.6   tidyverse_1.3.1

# Machine 2
> sessionInfo()
R version 4.2.1 (2022-06-23 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.utf8  LC_CTYPE=English_United Kingdom.utf8    LC_MONETARY=English_United Kingdom.utf8 LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.utf8    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] brms_2.17.0     Rcpp_1.0.9      forcats_0.5.2   stringr_1.4.1   dplyr_1.0.9     purrr_0.3.4     readr_2.1.2     tidyr_1.2.0     tibble_3.1.8   
[10] ggplot2_3.3.6   tidyverse_1.3.2

Start by reproducing the supplied right-censoring example and inspect the generated Stan code, focusing on normal_lccdf and comparing it with the working normal_lcdf path. Trace the corresponding Stan Math implementation and add a regression test showing that right-censored normal models initialize and sample correctly.