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

Ouverte
  • Quelque chose ne va pas ?
#3,486 0 commentaires 0 réactions 0 personnes assignées Voir sur GitHub

Les mainteneurs répondent en général sous 1 jour

Personne n'a encore pris cette issue.

Évaluation

Difficulté
4/5
Temps estimé
3-5 jours
Accessibilité débutants
54/100
Type d'issue
Bug
Clarté
Plutôt claire
Activité
Active
Stack technique
cpp
Domaine
backend, databases

Piste de recherche

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.

Rédigé par le modèle d'indexation à partir du texte de l'issue.

Description

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

Langage dominant
C++
Étoiles
845
Forks
225
Merge moyen
4 j 21 h
PR mergées (30 j)
13

Préparer son environnement

Par où commencer

  1. Lisez l'issue en entier, puis le guide de contribution du projet.
  2. Signalez en commentaire que vous la prenez — cela évite que deux personnes fassent le même travail.
  3. Forkez le dépôt et travaillez sur une branche.
  4. Ouvrez une pull request qui référence le numéro de l'issue.

Autres issues de stan-dev/math

Toutes les issues de stan-dev/math

Issues similaires

Plus d'issues C++

Recevez les nouvelles issues par e-mail

Un résumé court des issues GitHub adaptées aux débutants.