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

Aperta
  • Qualcosa non va?
#3,486 0 commenti 0 reazioni 0 assegnatari Vedi su GitHub

I maintainer di solito rispondono entro 1 giorno

Nessuno ha ancora preso questa issue.

Valutazione

Difficoltà
4/5
Tempo stimato
3-5 giorni
Idoneità per principianti
54/100
Tipo di issue
Bug
Chiarezza
Abbastanza chiara
Stato di attività
Attiva
Stack tecnologico
cpp
Ambito
backend, databases

Direzione di ricerca

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.

Scritto dal modello di indicizzazione a partire dal testo della issue.

Descrizione

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).

Lingua principale
C++
Stelle
845
Fork
225
Merge medio
3g 9h
PR unite (30g)
12

Preparare l'ambiente

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di stan-dev/math

Tutte le issue di stan-dev/math

Issue simili

Altre issue su C++

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.