library(cmdstanr) P <- 15 N <- 10000 ## N <- 1000 coefs <- abs(runif(P, .6, .99)) nu <- runif(P, -2, 2) re_sd <- runif(P, .3, 1.5) ones <- rep(1, N) latent_factors <- rnorm(N) Yhat <- ones %*% t(nu) + latent_factors %*% t(coefs) Y <- Yhat + matrix(rnorm(N*P, 0, re_sd), nrow = N, ncol = P, byrow = TRUE) mod <- cmdstan_model("./mvn_test_cfa.stan") fit_flag <- mod$sample(list(P = P, N = N, y = Y, flag = 1), parallel_chains = 4) fit_no_flag <- mod$sample(list(P = P, N = N, y = Y, flag = 0), parallel_chains = 4)