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

Offen
  • Stimmt etwas nicht?
#3,486 0 Kommentare 0 Reaktionen 0 zugewiesene Personen Auf GitHub ansehen

Maintainer antworten meist innerhalb von 1 Tag

Dieses Issue hat noch niemand übernommen.

Bewertung

Schwierigkeit
4/5
Geschätzter Aufwand
3-5 Tage
Anfängerfreundlichkeit
54/100
Issue-Typ
Bug
Klarheit
Größtenteils klar
Aktivitätsstatus
Aktiv
Tech-Stack
cpp
Bereich
backend, databases

Rechercherichtung

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.

Vom Indexierungsmodell aus dem Issue-Text verfasst.

Beschreibung

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

Vorherrschende Sprache
C++
Sterne
845
Forks
225
Ø Merge
4 T. 21 Std.
Gemergte PRs (30 T.)
13

Entwicklungsumgebung

Erste Schritte

  1. Lesen Sie das ganze Issue und danach den Beitragsleitfaden des Projekts.
  2. Schreiben Sie ins Issue, dass Sie es übernehmen — das erspart doppelte Arbeit.
  3. Forken Sie das Repository und arbeiten Sie in einem Branch.
  4. Öffnen Sie einen Pull Request, der die Issue-Nummer nennt.

Mehr aus stan-dev/math

Alle Issues in stan-dev/math

Ähnliche Issues

Weitere Issues zu C++

Neue Issues direkt in Ihr Postfach

Eine kurze Übersicht über anfängerfreundliche GitHub-Issues.