unconstrain_pars() sometimes returns infinite values

Summary:

unconstrain_pars() may return infinities for constrained samples that (possibly due to numerical problems) are on the bound of the parameter space. This can cause downstream problems, e.g., when attempting to use the unconstrained prameters for bridge sampling.

Description:

I fit a model using brms and tried to estimate the model's marginal likelihood using bridge_sampling(), but got the following error:

Error in out[!from@positive] <- -out[!from@positive] : 
  NAs are not allowed in subscripted assignments
In addition: Warning messages:
1: effective sample size cannot be calculated, has been replaced by number of samples. 
2: 1 of the 60000 log_prob() evaluations on the posterior draws produced -Inf/Inf. 
3: 100 of the 60000 log_prob() evaluations on the proposal draws produced -Inf/Inf. 
Error: Bridgesampling failed. Perhaps you did not set 'save_pars = save_pars(all = TRUE)' when fitting your model?

After some digging, I realized that the reason was an infinite value returned by unconstrain_pars(). Specifically, in one MCMC draw contains a correlation coefficient estimate of exactly 1 was backtransformed via atanh() to Inf.

Reproducible Steps:

The model was fit using cmdstanr, so this requires a recent version of brms. I use:

 package        * version    date       lib source        
 brms           * 2.16.1     2021-08-23 [1] CRAN (R 4.1.1)

Download the stanfit-object here.

The following code reproduces the error in the bridge sample:

library("brms")
brms_samples <- readRDS("brms_samples.rds")
bridge_sampler(brms_samples)

The following shows where the infinite unconstrained value arises:

library("rstan")

# Internal brms code
samples <- brms:::restructure(brms_samples)
samples <- brms:::update_misc_env(samples)
samples <- samples$fit
stanfit_model <- samples

# Internal bridgesampling code
ex <- extract(samples, permuted = FALSE)
skeleton <- rstan:::.create_skeleton(samples@sim$pars_oi,
                                              samples@par_dims[samples@sim$pars_oi])
upars <- apply(ex, 1:2, FUN = function(theta) {
  unconstrain_pars(stanfit_model, rstan:::.rstan_relist(theta, skeleton))
})

Current Output:

Because the object returned by unconstrain_pars() does not give parameter names, I had to look up the parameter name in the constrained samples:

> which(is.infinite(upars), arr.ind = T)
         iterations chains
[1,] 131      14398      3


# Unconstrained samples that include parameter names
ex2 <- extract(samples, c("b", "Intercept", "sd_1", "z_1", "L_1"), permuted = FALSE)

# The index is shifted by 1 becasue of the additional `cholesky_factor_corr`
# components in the constrained smaples
> ex2[14398, 3, 132, drop = FALSE]
, , parameters = L_1[2,1]

          chains
iterations chain:3
      [1,]       1

> atanh(ex2[14398, 3, 132])
[1] Inf

Expected Output:

An appropriate finite value.

RStan Version:

2.21.2

R Version:

4.1.1 (2021-08-10)

Operating System:

Ubuntu 18.04

Start with the reported rstan::unconstrain_pars() path and the rstan:::.rstan_relist() call in the reproduction, using the supplied brms_samples.rds example. Trace how a constrained correlation value of exactly 1 is transformed and add coverage for this boundary case. Done means the reproduced draw yields a finite unconstrained value without breaking existing transformations.