`cholesky_corr_free` rejects `cholesky_factor_corr` draws read back from CmdStan CSV with the default `sig_figs`

Abierto
  • ¿Algo no está bien?
#3,486 0 comentarios 0 reacciones 0 asignados Ver en GitHub

Los mantenedores suelen responder en 1 día

Nadie ha tomado este issue todavía.

Evaluación

Dificultad
4/5
Tiempo estimado
3-5 días
Aptitud para principiantes
54/100
Tipo de issue
Error
Claridad
Bastante claro
Estado de actividad
Activo
Stack tecnológico
cpp
Área
backend, databases

Línea de trabajo

Read cholesky_corr_free and check_cholesky_factor_corr, especially the change in #3387, then trace how rounded CSV draws reach these checks. Reproduce with the provided R example using CmdStan 2.40.0 and the default sig_figs; compare behavior at 18 digits. Done means rounded draws are handled without weakening rejection of invalid factors, with regression coverage for both cases.

Escrito por el modelo de indexación a partir del texto del issue.

Descripción

Description

With CmdStan 2.40.0, standalone generated quantities and cmdstanr $unconstrain_draws() fail for models with a cholesky_factor_corr parameter. With CmdStan 2.39.0, the same code works.

cholesky_corr_free: x is not a valid unit vector. The sum of the squares of the elements should be 1, but is 1

Since #3387, cholesky_corr_free calls check_cholesky_factor_corr, which requires the sum of squares of each row to be within CONSTRAINT_TOLERANCE (1e-8) of 1. CmdStan writes the draws with 8 significant digits by default (stan-dev/cmdstan#1308, since 2.37.0). In the example below, this rounding moves the sum of squares up to 2.0e-8 away from 1. The message shows the sum as "1" because it prints only 6 significant digits.

Example
library(cmdstanr)
code <- "
parameters {
  cholesky_factor_corr[3] L;
}
model {
  L ~ lkj_corr_cholesky(2);
}
generated quantities {
  real rho = multiply_lower_tri_self_transpose(L)[2, 1];
}
"
mod <- cmdstan_model(write_stan_file(code), compile_model_methods = TRUE)
fit <- mod$sample(chains = 1, iter_sampling = 1000, seed = 1)
mod$generate_quantities(fit, seed = 1)$draws()
fit$unconstrain_draws()
CmdStan sig_figs generate_quantities() unconstrain_draws()
2.39.0 default (8) OK OK
2.39.0 18 OK OK
2.40.0 default (8) error error
2.40.0 18 OK OK

To see how often the check fails, I rounded random K x K correlation Cholesky factors (cov2cor of a Wishart draw) with signif(). The table gives the share of factors that fail the check, and the largest deviation at 8 digits:

K 6 digits 8 digits 9 digits max |1 - sum of squares| at 8 digits
3 100% 10.8% 0% 1.4e-8
10 100% 45.1% 0% 1.9e-8
40 100% 94.2% 0% 2.0e-8

If 1000 draws failed independently at the rate of 10.8% for K = 3, a run would contain a failing draw with a probability above 99.99%. The example run above has K = 3 and fails.

set.seed(1)
rl <- function(K) { A <- matrix(rnorm(K * (K + 5)), K + 5); t(chol(cov2cor(crossprod(A)))) }
fails <- function(K, sf, n = 4000) mean(replicate(n, any(abs(1 - rowSums(signif(rl(K), sf)^2)) > 1e-8)))
fails(3, 8)
Possible fixes
  1. In cholesky_corr_free, keep the checks for a square, lower-triangular matrix with a positive diagonal. Reject a row only if its sum of squares differs from 1 by more than a looser tolerance, for example 1e-4. Then scale each row to unit length before the transform. Rounded draws pass, and the example in #3380 (a row with a sum of squares of 0.5) still fails. This would accept also csv files written by older CmdStan versions with sig_figs = 6. A looser tolerance without the scaling would also remove the error.
  2. Increase the CmdStan default sig_figs to 9 (or bigger value). In the simulation, 9 digits gave no failures for K = 3, 10 and 40. This helps only CSV files written after the change, and each value gets one character longer.
Current Version:

v5.4.0 (CmdStan 2.40.0). Tested with R 4.6.1 and cmdstanr 0.9.0.9001 on Linux (CSC Roihu).

Lenguaje dominante
C++
Estrellas
845
Forks
225
Merge medio
3 d 9 h
PR fusionados (30 d)
12

Preparar el entorno

Primeros pasos

  1. Lee el issue completo y luego la guía de contribución del proyecto.
  2. Comenta en el issue que vas a ocuparte — evita que dos personas hagan lo mismo.
  3. Haz un fork del repositorio y trabaja en una rama.
  4. Abre un pull request que haga referencia al número del issue.

Más de stan-dev/math

Todos los issues de stan-dev/math

Issues similares

Más issues de C++

Recibe los nuevos issues en tu correo

Un resumen breve de issues de GitHub para principiantes.